Mercurial > repos > public > sbplib

annotate +multiblock/LaplaceSquared.m @ 1198:2924b3a9b921 feature/d2_compatible

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add OpSet for fully compatible D2Variable, created from regular D2Variable by replacing d1 by first row of D1. Formal reduction by one order of accuracy at the boundary point.
author Martin Almquist <malmquist@stanford.edu>
date Fri, 16 Aug 2019 14:30:28 -0700
parents 2aaced07d1e5
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
1070
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
1 classdef LaplaceSquared < scheme.Scheme
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
2 properties
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
3 grid
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
4 order
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
5 laplaceDiffOp
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
6
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
7 D
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
8 H
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
9 Hi
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
10
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
11 a,b
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
12 end
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
13
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
14 methods
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
15 % Discretisation of a*nabla*b*nabla
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
16 function obj = LaplaceSquared(g, order, a, b, opGen)
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
17 default_arg('order', 4);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
18 default_arg('a', 1);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
19 default_arg('b', 1);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
20 default_arg('opGen', @sbp.D4Variable);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
21
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
22 if isscalar(a)
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
23 a = grid.evalOn(g, a);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
24 end
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
25
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
26 if isscalar(b)
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
27 b = grid.evalOn(g, b);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
28 end
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
29
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
30 obj.grid = g;
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
31 obj.order = order;
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
32 obj.a = a;
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
33 obj.b = b;
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
34
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
35 obj.laplaceDiffOp = multiblock.Laplace(g, order, 1, 1, opGen);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
36
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
37 obj.H = obj.laplaceDiffOp.H;
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
38 obj.Hi = spdiag(1./diag(obj.H));
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
39
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
40 A = spdiag(a);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
41 B = spdiag(b);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
42
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
43 D_laplace = obj.laplaceDiffOp.D;
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
44 obj.D = A*D_laplace*B*D_laplace;
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
45 end
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
46
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
47 function s = size(obj)
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
48 s = size(obj.laplaceDiffOp);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
49 end
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
50
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
51 function op = getBoundaryOperator(obj, opName, boundary)
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
52 switch opName
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
53 case 'e'
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
54 op = getBoundaryOperator(obj.laplaceDiffOp, 'e', boundary);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
55 case 'd1'
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
56 op = getBoundaryOperator(obj.laplaceDiffOp, 'd', boundary);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
57 case 'd2'
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
58 e = getBoundaryOperator(obj.laplaceDiffOp, 'e', boundary);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
59 op = (e'*obj.laplaceDiffOp.D)';
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
60 case 'd3'
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
61 d1 = getBoundaryOperator(obj.laplaceDiffOp, 'd', boundary);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
62 op = (d1'*spdiag(obj.b)*obj.laplaceDiffOp.D)';
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
63 end
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
64 end
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
65
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
66 function op = getBoundaryQuadrature(obj, boundary)
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
67 op = getBoundaryQuadrature(obj.laplaceDiffOp, boundary);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
68 end
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
69
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
70 function [closure, penalty] = boundary_condition(obj,boundary,type) % TODO: Change name to boundaryCondition
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
71 switch type
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
72 case 'e'
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
73 error('Bc of type ''e'' not implemented')
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
74 case 'd1'
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
75 error('Bc of type ''d1'' not implemented')
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
76 case 'd2'
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
77 e = obj.getBoundaryOperator('e', boundary);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
78 d1 = obj.getBoundaryOperator('d1', boundary);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
79 d2 = obj.getBoundaryOperator('d2', boundary);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
80 H_b = obj.getBoundaryQuadrature(boundary);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
81
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
82 A = spdiag(obj.a);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
83 B_b = spdiag(e'*obj.b);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
84
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
85 tau = obj.Hi*A*d1*B_b*H_b;
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
86 closure = tau*d2';
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
87 penalty = -tau;
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
88 case 'd3'
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
89 e = obj.getBoundaryOperator('e', boundary);
1073
95113a592421 Fix bug in d3 boundary condition
Jonatan Werpers <jonatan@werpers.com>
parents: 1070
diff changeset
90 d3 = obj.getBoundaryOperator('d3', boundary);
1070
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
91 H_b = obj.getBoundaryQuadrature(boundary);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
92
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
93 A = spdiag(obj.a);
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
94
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
95 tau = -obj.Hi*A*e*H_b;
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
96 closure = tau*d3';
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
97 penalty = -tau;
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
98 end
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
99 end
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
100
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
101 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
102 error('Not implemented')
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
103 end
f6b3af6febf3 Start work on LaplaceSquared. FULL OF ERRORS?!
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
104 end
1074
2aaced07d1e5 Formatting fix
Jonatan Werpers <jonatan@werpers.com>
parents: 1073
diff changeset
105 end