Mercurial > repos > public > sbplib

comparison +multiblock/LaplaceSquared.m @ 1076:c7508ba70527

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge in feature/grids/LaplaceSquared
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 13 Feb 2019 20:52:14 +0100
parents 2aaced07d1e5
children
comparison
equal deleted inserted replaced
1061:c7b619cf5e34 1076:c7508ba70527
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('d3', 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