Mercurial > repos > public > sbplib

comparison +multiblock/LaplaceSquared.m @ 1108:5ec23b9bf360 feature/laplace_curvilinear_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default
author Martin Almquist <malmquist@stanford.edu>
date Wed, 10 Apr 2019 11:00:27 -0700
parents 2aaced07d1e5
children
comparison
equal deleted inserted replaced
1087:867307f4d80f 1108:5ec23b9bf360
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