Mercurial > repos > public > sbplib

comparison +multiblock/LaplaceSquared.m @ 1070:f6b3af6febf3 feature/grids/LaplaceSquared

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Start work on LaplaceSquared. FULL OF ERRORS?!
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 24 Jul 2018 21:08:08 -0700
parents
children 95113a592421
comparison
equal deleted inserted replaced
784:085fc0fe537f 1070:f6b3af6febf3
1 classdef LaplaceSquared < scheme.Scheme
2 properties
3 grid
4 order
5 laplaceDiffOp
6
7 D
8 H
9 Hi
10
11 a,b
12 end
13
14 methods
15 % Discretisation of a*nabla*b*nabla
16 function obj = LaplaceSquared(g, order, a, b, opGen)
17 default_arg('order', 4);
18 default_arg('a', 1);
19 default_arg('b', 1);
20 default_arg('opGen', @sbp.D4Variable);
21
22 if isscalar(a)
23 a = grid.evalOn(g, a);
24 end
25
26 if isscalar(b)
27 b = grid.evalOn(g, b);
28 end
29
30 obj.grid = g;
31 obj.order = order;
32 obj.a = a;
33 obj.b = b;
34
35 obj.laplaceDiffOp = multiblock.Laplace(g, order, 1, 1, opGen);
36
37 obj.H = obj.laplaceDiffOp.H;
38 obj.Hi = spdiag(1./diag(obj.H));
39
40 A = spdiag(a);
41 B = spdiag(b);
42
43 D_laplace = obj.laplaceDiffOp.D;
44 obj.D = A*D_laplace*B*D_laplace;
45 end
46
47 function s = size(obj)
48 s = size(obj.laplaceDiffOp);
49 end
50
51 function op = getBoundaryOperator(obj, opName, boundary)
52 switch opName
53 case 'e'
54 op = getBoundaryOperator(obj.laplaceDiffOp, 'e', boundary);
55 case 'd1'
56 op = getBoundaryOperator(obj.laplaceDiffOp, 'd', boundary);
57 case 'd2'
58 e = getBoundaryOperator(obj.laplaceDiffOp, 'e', boundary);
59 op = (e'*obj.laplaceDiffOp.D)';
60 case 'd3'
61 d1 = getBoundaryOperator(obj.laplaceDiffOp, 'd', boundary);
62 op = (d1'*spdiag(obj.b)*obj.laplaceDiffOp.D)';
63 end
64 end
65
66 function op = getBoundaryQuadrature(obj, boundary)
67 op = getBoundaryQuadrature(obj.laplaceDiffOp, boundary);
68 end
69
70 function [closure, penalty] = boundary_condition(obj,boundary,type) % TODO: Change name to boundaryCondition
71 switch type
72 case 'e'
73 error('Bc of type ''e'' not implemented')
74 case 'd1'
75 error('Bc of type ''d1'' not implemented')
76 case 'd2'
77 e = obj.getBoundaryOperator('e', boundary);
78 d1 = obj.getBoundaryOperator('d1', boundary);
79 d2 = obj.getBoundaryOperator('d2', boundary);
80 H_b = obj.getBoundaryQuadrature(boundary);
81
82 A = spdiag(obj.a);
83 B_b = spdiag(e'*obj.b);
84
85 tau = obj.Hi*A*d1*B_b*H_b;
86 closure = tau*d2';
87 penalty = -tau;
88 case 'd3'
89 e = obj.getBoundaryOperator('e', boundary);
90 d3 = obj.getBoundaryOperator('d1', boundary);
91 H_b = obj.getBoundaryQuadrature(boundary);
92
93 A = spdiag(obj.a);
94
95 tau = -obj.Hi*A*e*H_b;
96 closure = tau*d3';
97 penalty = -tau;
98 end
99 end
100
101 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
102 error('Not implemented')
103 end
104 end
105 end