Mercurial > repos > public > sbplib
annotate +scheme/LaplaceCurvilinear.m @ 1031:2ef20d00b386 feature/advectionRV
For easier comparison, return both the first order and residual viscosity when evaluating the residual. Add the first order and residual viscosity to the state of the RungekuttaRV time steppers
author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
---|---|
date | Thu, 17 Jan 2019 10:25:06 +0100 |
parents | 706d1c2b4199 |
children | a0b3161e44f3 |
rev | line source |
---|---|
450
8d455e49364f
Copy Wave2dCurve to new scheme LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
449
diff
changeset
|
1 classdef LaplaceCurvilinear < scheme.Scheme |
0 | 2 properties |
3 m % Number of points in each direction, possibly a vector | |
4 h % Grid spacing | |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
5 |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
6 grid |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
7 |
0 | 8 order % Order accuracy for the approximation |
9 | |
552
ffaf13533c27
Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents:
551
diff
changeset
|
10 a,b % Parameters of the operator |
ffaf13533c27
Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents:
551
diff
changeset
|
11 |
ffaf13533c27
Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents:
551
diff
changeset
|
12 |
ffaf13533c27
Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents:
551
diff
changeset
|
13 % Inner products and operators for physical coordinates |
ffaf13533c27
Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents:
551
diff
changeset
|
14 D % Laplace operator |
ffaf13533c27
Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents:
551
diff
changeset
|
15 H, Hi % Inner product |
ffaf13533c27
Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents:
551
diff
changeset
|
16 e_w, e_e, e_s, e_n |
ffaf13533c27
Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents:
551
diff
changeset
|
17 d_w, d_e, d_s, d_n % Normal derivatives at the boundary |
ffaf13533c27
Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents:
551
diff
changeset
|
18 H_w, H_e, H_s, H_n % Boundary inner products |
ffaf13533c27
Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents:
551
diff
changeset
|
19 Dx, Dy % Physical derivatives |
ffaf13533c27
Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents:
551
diff
changeset
|
20 M % Gradient inner product |
ffaf13533c27
Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents:
551
diff
changeset
|
21 |
ffaf13533c27
Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents:
551
diff
changeset
|
22 % Metric coefficients |
0 | 23 J, Ji |
24 a11, a12, a22 | |
552
ffaf13533c27
Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents:
551
diff
changeset
|
25 x_u |
ffaf13533c27
Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents:
551
diff
changeset
|
26 x_v |
ffaf13533c27
Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents:
551
diff
changeset
|
27 y_u |
ffaf13533c27
Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents:
551
diff
changeset
|
28 y_v |
0 | 29 |
552
ffaf13533c27
Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents:
551
diff
changeset
|
30 % Inner product and operators for logical coordinates |
0 | 31 H_u, H_v % Norms in the x and y directions |
552
ffaf13533c27
Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents:
551
diff
changeset
|
32 Hi_u, Hi_v |
0 | 33 Hu,Hv % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. |
34 Hiu, Hiv | |
35 du_w, dv_w | |
36 du_e, dv_e | |
37 du_s, dv_s | |
38 du_n, dv_n | |
39 gamm_u, gamm_v | |
27
97a638f91fb8
Added function spdiag(). Fixed a bunch of bugs in the Wave2dCurve scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
5
diff
changeset
|
40 lambda |
720
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
41 |
0 | 42 end |
43 | |
44 methods | |
452
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
45 % Implements a*div(b*grad(u)) as a SBP scheme |
539 | 46 % TODO: Implement proper H, it should be the real physical quadrature, the logic quadrature may be but in a separate variable (H_logic?) |
47 | |
906
0499239496cf
Remove property interpolation_type in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
904
diff
changeset
|
48 function obj = LaplaceCurvilinear(g ,order, a, b, opSet) |
303
f18142c1530b
Fixed Wave2dCurve for new operators. Some cleanup in d2_variable_4.
Jonatan Werpers <jonatan@werpers.com>
parents:
96
diff
changeset
|
49 default_arg('opSet',@sbp.D2Variable); |
452
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
50 default_arg('a', 1); |
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
51 default_arg('b', 1); |
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
52 |
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
53 if b ~=1 |
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
54 error('Not implemented yet') |
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
55 end |
0 | 56 |
720
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
57 % assert(isa(g, 'grid.Curvilinear')) |
743
f4595f14d696
Change schemes to work for special coefficients.
Martin Almquist <malmquist@stanford.edu>
parents:
720
diff
changeset
|
58 if isa(a, 'function_handle') |
f4595f14d696
Change schemes to work for special coefficients.
Martin Almquist <malmquist@stanford.edu>
parents:
720
diff
changeset
|
59 a = grid.evalOn(g, a); |
f4595f14d696
Change schemes to work for special coefficients.
Martin Almquist <malmquist@stanford.edu>
parents:
720
diff
changeset
|
60 a = spdiag(a); |
f4595f14d696
Change schemes to work for special coefficients.
Martin Almquist <malmquist@stanford.edu>
parents:
720
diff
changeset
|
61 end |
0 | 62 |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
63 m = g.size(); |
0 | 64 m_u = m(1); |
65 m_v = m(2); | |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
66 m_tot = g.N(); |
0 | 67 |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
68 h = g.scaling(); |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
69 h_u = h(1); |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
70 h_v = h(2); |
0 | 71 |
553
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
72 |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
73 % 1D operators |
303
f18142c1530b
Fixed Wave2dCurve for new operators. Some cleanup in d2_variable_4.
Jonatan Werpers <jonatan@werpers.com>
parents:
96
diff
changeset
|
74 ops_u = opSet(m_u, {0, 1}, order); |
f18142c1530b
Fixed Wave2dCurve for new operators. Some cleanup in d2_variable_4.
Jonatan Werpers <jonatan@werpers.com>
parents:
96
diff
changeset
|
75 ops_v = opSet(m_v, {0, 1}, order); |
0 | 76 |
77 I_u = speye(m_u); | |
78 I_v = speye(m_v); | |
79 | |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
80 D1_u = ops_u.D1; |
303
f18142c1530b
Fixed Wave2dCurve for new operators. Some cleanup in d2_variable_4.
Jonatan Werpers <jonatan@werpers.com>
parents:
96
diff
changeset
|
81 D2_u = ops_u.D2; |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
82 H_u = ops_u.H; |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
83 Hi_u = ops_u.HI; |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
84 e_l_u = ops_u.e_l; |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
85 e_r_u = ops_u.e_r; |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
86 d1_l_u = ops_u.d1_l; |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
87 d1_r_u = ops_u.d1_r; |
0 | 88 |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
89 D1_v = ops_v.D1; |
303
f18142c1530b
Fixed Wave2dCurve for new operators. Some cleanup in d2_variable_4.
Jonatan Werpers <jonatan@werpers.com>
parents:
96
diff
changeset
|
90 D2_v = ops_v.D2; |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
91 H_v = ops_v.H; |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
92 Hi_v = ops_v.HI; |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
93 e_l_v = ops_v.e_l; |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
94 e_r_v = ops_v.e_r; |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
95 d1_l_v = ops_v.d1_l; |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
96 d1_r_v = ops_v.d1_r; |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
97 |
553
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
98 |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
99 % Logical operators |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
100 Du = kr(D1_u,I_v); |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
101 Dv = kr(I_u,D1_v); |
553
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
102 obj.Hu = kr(H_u,I_v); |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
103 obj.Hv = kr(I_u,H_v); |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
104 obj.Hiu = kr(Hi_u,I_v); |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
105 obj.Hiv = kr(I_u,Hi_v); |
0 | 106 |
553
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
107 e_w = kr(e_l_u,I_v); |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
108 e_e = kr(e_r_u,I_v); |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
109 e_s = kr(I_u,e_l_v); |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
110 e_n = kr(I_u,e_r_v); |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
111 obj.du_w = kr(d1_l_u,I_v); |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
112 obj.dv_w = (e_w'*Dv)'; |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
113 obj.du_e = kr(d1_r_u,I_v); |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
114 obj.dv_e = (e_e'*Dv)'; |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
115 obj.du_s = (e_s'*Du)'; |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
116 obj.dv_s = kr(I_u,d1_l_v); |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
117 obj.du_n = (e_n'*Du)'; |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
118 obj.dv_n = kr(I_u,d1_r_v); |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
119 |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
120 |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
121 % Metric coefficients |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
122 coords = g.points(); |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
123 x = coords(:,1); |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
124 y = coords(:,2); |
0 | 125 |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
126 x_u = Du*x; |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
127 x_v = Dv*x; |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
128 y_u = Du*y; |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
129 y_v = Dv*y; |
0 | 130 |
131 J = x_u.*y_v - x_v.*y_u; | |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
132 a11 = 1./J .* (x_v.^2 + y_v.^2); |
0 | 133 a12 = -1./J .* (x_u.*x_v + y_u.*y_v); |
134 a22 = 1./J .* (x_u.^2 + y_u.^2); | |
135 lambda = 1/2 * (a11 + a22 - sqrt((a11-a22).^2 + 4*a12.^2)); | |
136 | |
553
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
137 obj.x_u = x_u; |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
138 obj.x_v = x_v; |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
139 obj.y_u = y_u; |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
140 obj.y_v = y_v; |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
141 |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
142 |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
143 % Assemble full operators |
551
c5a7a13c03dc
Switch to using spdiag instead of spdiags
Jonatan Werpers <jonatan@werpers.com>
parents:
539
diff
changeset
|
144 L_12 = spdiag(a12); |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
145 Duv = Du*L_12*Dv; |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
146 Dvu = Dv*L_12*Du; |
0 | 147 |
148 Duu = sparse(m_tot); | |
149 Dvv = sparse(m_tot); | |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
150 ind = grid.funcToMatrix(g, 1:m_tot); |
0 | 151 |
152 for i = 1:m_v | |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
153 D = D2_u(a11(ind(:,i))); |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
154 p = ind(:,i); |
0 | 155 Duu(p,p) = D; |
156 end | |
157 | |
158 for i = 1:m_u | |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
159 D = D2_v(a22(ind(i,:))); |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
160 p = ind(i,:); |
0 | 161 Dvv(p,p) = D; |
162 end | |
163 | |
553
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
164 |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
165 % Physical operators |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
166 obj.J = spdiag(J); |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
167 obj.Ji = spdiag(1./J); |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
168 |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
169 obj.D = obj.Ji*a*(Duu + Duv + Dvu + Dvv); |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
170 obj.H = obj.J*kr(H_u,H_v); |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
171 obj.Hi = obj.Ji*kr(Hi_u,Hi_v); |
0 | 172 |
553
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
173 obj.e_w = e_w; |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
174 obj.e_e = e_e; |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
175 obj.e_s = e_s; |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
176 obj.e_n = e_n; |
376 | 177 |
557
2a856a589510
Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents:
554
diff
changeset
|
178 %% normal derivatives |
2a856a589510
Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents:
554
diff
changeset
|
179 I_w = ind(1,:); |
2a856a589510
Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents:
554
diff
changeset
|
180 I_e = ind(end,:); |
2a856a589510
Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents:
554
diff
changeset
|
181 I_s = ind(:,1); |
2a856a589510
Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents:
554
diff
changeset
|
182 I_n = ind(:,end); |
2a856a589510
Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents:
554
diff
changeset
|
183 |
2a856a589510
Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents:
554
diff
changeset
|
184 a11_w = spdiag(a11(I_w)); |
2a856a589510
Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents:
554
diff
changeset
|
185 a12_w = spdiag(a12(I_w)); |
2a856a589510
Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents:
554
diff
changeset
|
186 a11_e = spdiag(a11(I_e)); |
2a856a589510
Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents:
554
diff
changeset
|
187 a12_e = spdiag(a12(I_e)); |
2a856a589510
Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents:
554
diff
changeset
|
188 a22_s = spdiag(a22(I_s)); |
2a856a589510
Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents:
554
diff
changeset
|
189 a12_s = spdiag(a12(I_s)); |
2a856a589510
Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents:
554
diff
changeset
|
190 a22_n = spdiag(a22(I_n)); |
2a856a589510
Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents:
554
diff
changeset
|
191 a12_n = spdiag(a12(I_n)); |
2a856a589510
Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents:
554
diff
changeset
|
192 |
565 | 193 s_w = sqrt((e_w'*x_v).^2 + (e_w'*y_v).^2); |
194 s_e = sqrt((e_e'*x_v).^2 + (e_e'*y_v).^2); | |
195 s_s = sqrt((e_s'*x_u).^2 + (e_s'*y_u).^2); | |
196 s_n = sqrt((e_n'*x_u).^2 + (e_n'*y_u).^2); | |
562
11d8d6ccbcd7
Add boundary inner products
Jonatan Werpers <jonatan@werpers.com>
parents:
561
diff
changeset
|
197 |
566
9c98a0526afc
Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents:
565
diff
changeset
|
198 obj.d_w = -1*(spdiag(1./s_w)*(a11_w*obj.du_w' + a12_w*obj.dv_w'))'; |
9c98a0526afc
Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents:
565
diff
changeset
|
199 obj.d_e = (spdiag(1./s_e)*(a11_e*obj.du_e' + a12_e*obj.dv_e'))'; |
9c98a0526afc
Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents:
565
diff
changeset
|
200 obj.d_s = -1*(spdiag(1./s_s)*(a22_s*obj.dv_s' + a12_s*obj.du_s'))'; |
9c98a0526afc
Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents:
565
diff
changeset
|
201 obj.d_n = (spdiag(1./s_n)*(a22_n*obj.dv_n' + a12_n*obj.du_n'))'; |
554
6e71d4f8b155
Better names for the metric coefficients
Jonatan Werpers <jonatan@werpers.com>
parents:
553
diff
changeset
|
202 |
553
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
203 obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv; |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
204 obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv; |
384
32df00102268
Fix bug in Characteristic BC for Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
360
diff
changeset
|
205 |
562
11d8d6ccbcd7
Add boundary inner products
Jonatan Werpers <jonatan@werpers.com>
parents:
561
diff
changeset
|
206 %% Boundary inner products |
564
9bf49338f8e6
Fix bug in creation of boundary inner products
Jonatan Werpers <jonatan@werpers.com>
parents:
563
diff
changeset
|
207 obj.H_w = H_v*spdiag(s_w); |
9bf49338f8e6
Fix bug in creation of boundary inner products
Jonatan Werpers <jonatan@werpers.com>
parents:
563
diff
changeset
|
208 obj.H_e = H_v*spdiag(s_e); |
9bf49338f8e6
Fix bug in creation of boundary inner products
Jonatan Werpers <jonatan@werpers.com>
parents:
563
diff
changeset
|
209 obj.H_s = H_u*spdiag(s_s); |
9bf49338f8e6
Fix bug in creation of boundary inner products
Jonatan Werpers <jonatan@werpers.com>
parents:
563
diff
changeset
|
210 obj.H_n = H_u*spdiag(s_n); |
553
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
211 |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
212 % Misc. |
0 | 213 obj.m = m; |
214 obj.h = [h_u h_v]; | |
215 obj.order = order; | |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
216 obj.grid = g; |
0 | 217 |
452
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
218 obj.a = a; |
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
219 obj.b = b; |
0 | 220 obj.a11 = a11; |
221 obj.a12 = a12; | |
222 obj.a22 = a22; | |
27
97a638f91fb8
Added function spdiag(). Fixed a bunch of bugs in the Wave2dCurve scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
5
diff
changeset
|
223 obj.lambda = lambda; |
0 | 224 |
389
42c89b5eedc0
Add borrowing constants for D2 operators in D4Variable
Jonatan Werpers <jonatan@werpers.com>
parents:
384
diff
changeset
|
225 obj.gamm_u = h_u*ops_u.borrowing.M.d1; |
42c89b5eedc0
Add borrowing constants for D2 operators in D4Variable
Jonatan Werpers <jonatan@werpers.com>
parents:
384
diff
changeset
|
226 obj.gamm_v = h_v*ops_v.borrowing.M.d1; |
0 | 227 end |
228 | |
229 | |
230 % Closure functions return the opertors applied to the own doamin to close the boundary | |
231 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. | |
232 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. | |
233 % type is a string specifying the type of boundary condition if there are several. | |
234 % data is a function returning the data that should be applied at the boundary. | |
235 % neighbour_scheme is an instance of Scheme that should be interfaced to. | |
236 % neighbour_boundary is a string specifying which boundary to interface to. | |
347
85c2fe06d551
Implemented characteristic form BC in Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
337
diff
changeset
|
237 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter) |
0 | 238 default_arg('type','neumann'); |
347
85c2fe06d551
Implemented characteristic form BC in Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
337
diff
changeset
|
239 default_arg('parameter', []); |
0 | 240 |
567
33b962620e24
Remove unneeded sign information from get_boundary_ops
Jonatan Werpers <jonatan@werpers.com>
parents:
566
diff
changeset
|
241 [e, d, gamm, H_b, ~] = obj.get_boundary_ops(boundary); |
0 | 242 switch type |
243 % Dirichlet boundary condition | |
244 case {'D','d','dirichlet'} | |
82
df642750e8f7
Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
55
diff
changeset
|
245 tuning = 1.2; |
df642750e8f7
Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
55
diff
changeset
|
246 % tuning = 20.2; |
df642750e8f7
Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
55
diff
changeset
|
247 |
566
9c98a0526afc
Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents:
565
diff
changeset
|
248 b1 = gamm*obj.lambda./obj.a11.^2; |
9c98a0526afc
Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents:
565
diff
changeset
|
249 b2 = gamm*obj.lambda./obj.a22.^2; |
0 | 250 |
566
9c98a0526afc
Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents:
565
diff
changeset
|
251 tau1 = tuning * spdiag(-1./b1 - 1./b2); |
9c98a0526afc
Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents:
565
diff
changeset
|
252 tau2 = 1; |
0 | 253 |
566
9c98a0526afc
Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents:
565
diff
changeset
|
254 tau = (tau1*e + tau2*d)*H_b; |
0 | 255 |
566
9c98a0526afc
Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents:
565
diff
changeset
|
256 closure = obj.a*obj.Hi*tau*e'; |
9c98a0526afc
Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents:
565
diff
changeset
|
257 penalty = -obj.a*obj.Hi*tau; |
0 | 258 |
259 | |
260 % Neumann boundary condition | |
261 case {'N','n','neumann'} | |
558
54c775c3348a
Clean up use of the sign in boundary and interface routines
Jonatan Werpers <jonatan@werpers.com>
parents:
557
diff
changeset
|
262 tau1 = -1; |
0 | 263 tau2 = 0; |
566
9c98a0526afc
Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents:
565
diff
changeset
|
264 tau = (tau1*e + tau2*d)*H_b; |
0 | 265 |
566
9c98a0526afc
Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents:
565
diff
changeset
|
266 closure = obj.a*obj.Hi*tau*d'; |
9c98a0526afc
Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents:
565
diff
changeset
|
267 penalty = -obj.a*obj.Hi*tau; |
0 | 268 |
269 | |
270 % Unknown, boundary condition | |
271 otherwise | |
272 error('No such boundary condition: type = %s',type); | |
273 end | |
274 end | |
275 | |
928
1c61d8fa9903
Replace opts by type everywhere
Martin Almquist <malmquist@stanford.edu>
parents:
927
diff
changeset
|
276 % type Struct that specifies the interface coupling. |
904
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
277 % Fields: |
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
278 % -- tuning: penalty strength, defaults to 1.2 |
935
1bf200417957
Fix doc for LaplaceCurvilinear.interface
Jonatan Werpers <jonatan@werpers.com>
parents:
928
diff
changeset
|
279 % -- interpolation: type of interpolation, default 'none' |
928
1c61d8fa9903
Replace opts by type everywhere
Martin Almquist <malmquist@stanford.edu>
parents:
927
diff
changeset
|
280 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) |
927
4291731570bb
Rename AWW OP. Make Interpolation operator classes take grid points as arguments. Remove LaplCurv.interpolationOperators. Introduce default struct in LaplCurv.interface.
Martin Almquist <malmquist@stanford.edu>
parents:
925
diff
changeset
|
281 |
4291731570bb
Rename AWW OP. Make Interpolation operator classes take grid points as arguments. Remove LaplCurv.interpolationOperators. Introduce default struct in LaplCurv.interface.
Martin Almquist <malmquist@stanford.edu>
parents:
925
diff
changeset
|
282 defaultType.tuning = 1.2; |
4291731570bb
Rename AWW OP. Make Interpolation operator classes take grid points as arguments. Remove LaplCurv.interpolationOperators. Introduce default struct in LaplCurv.interface.
Martin Almquist <malmquist@stanford.edu>
parents:
925
diff
changeset
|
283 defaultType.interpolation = 'none'; |
928
1c61d8fa9903
Replace opts by type everywhere
Martin Almquist <malmquist@stanford.edu>
parents:
927
diff
changeset
|
284 default_struct('type', defaultType); |
927
4291731570bb
Rename AWW OP. Make Interpolation operator classes take grid points as arguments. Remove LaplCurv.interpolationOperators. Introduce default struct in LaplCurv.interface.
Martin Almquist <malmquist@stanford.edu>
parents:
925
diff
changeset
|
285 |
928
1c61d8fa9903
Replace opts by type everywhere
Martin Almquist <malmquist@stanford.edu>
parents:
927
diff
changeset
|
286 switch type.interpolation |
927
4291731570bb
Rename AWW OP. Make Interpolation operator classes take grid points as arguments. Remove LaplCurv.interpolationOperators. Introduce default struct in LaplCurv.interface.
Martin Almquist <malmquist@stanford.edu>
parents:
925
diff
changeset
|
287 case {'none', ''} |
940
31186559236d
Add forgotten type argument in call to interfaceStandard in Schrodinger2d and LaplCurv.
Martin Almquist <malmquist@stanford.edu>
parents:
938
diff
changeset
|
288 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type); |
927
4291731570bb
Rename AWW OP. Make Interpolation operator classes take grid points as arguments. Remove LaplCurv.interpolationOperators. Introduce default struct in LaplCurv.interface.
Martin Almquist <malmquist@stanford.edu>
parents:
925
diff
changeset
|
289 case {'op','OP'} |
928
1c61d8fa9903
Replace opts by type everywhere
Martin Almquist <malmquist@stanford.edu>
parents:
927
diff
changeset
|
290 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type); |
927
4291731570bb
Rename AWW OP. Make Interpolation operator classes take grid points as arguments. Remove LaplCurv.interpolationOperators. Introduce default struct in LaplCurv.interface.
Martin Almquist <malmquist@stanford.edu>
parents:
925
diff
changeset
|
291 otherwise |
4291731570bb
Rename AWW OP. Make Interpolation operator classes take grid points as arguments. Remove LaplCurv.interpolationOperators. Introduce default struct in LaplCurv.interface.
Martin Almquist <malmquist@stanford.edu>
parents:
925
diff
changeset
|
292 error('Unknown type of interpolation: %s ', type.interpolation); |
904
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
293 end |
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
294 end |
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
295 |
928
1c61d8fa9903
Replace opts by type everywhere
Martin Almquist <malmquist@stanford.edu>
parents:
927
diff
changeset
|
296 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type) |
1c61d8fa9903
Replace opts by type everywhere
Martin Almquist <malmquist@stanford.edu>
parents:
927
diff
changeset
|
297 tuning = type.tuning; |
904
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
298 |
906
0499239496cf
Remove property interpolation_type in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
904
diff
changeset
|
299 % u denotes the solution in the own domain |
0499239496cf
Remove property interpolation_type in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
904
diff
changeset
|
300 % v denotes the solution in the neighbour domain |
904
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
301 [e_u, d_u, gamm_u, H_b_u, I_u] = obj.get_boundary_ops(boundary); |
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
302 [e_v, d_v, gamm_v, H_b_v, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); |
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
303 |
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
304 u = obj; |
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
305 v = neighbour_scheme; |
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
306 |
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
307 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2; |
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
308 b2_u = gamm_u*u.lambda(I_u)./u.a22(I_u).^2; |
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
309 b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2; |
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
310 b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2; |
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
311 |
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
312 tau1 = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v); |
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
313 tau1 = tuning * spdiag(tau1); |
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
314 tau2 = 1/2; |
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
315 |
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
316 sig1 = -1/2; |
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
317 sig2 = 0; |
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
318 |
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
319 tau = (e_u*tau1 + tau2*d_u)*H_b_u; |
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
320 sig = (sig1*e_u + sig2*d_u)*H_b_u; |
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
321 |
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
322 closure = obj.a*obj.Hi*( tau*e_u' + sig*d_u'); |
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
323 penalty = obj.a*obj.Hi*(-tau*e_v' + sig*d_v'); |
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
324 end |
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
325 |
928
1c61d8fa9903
Replace opts by type everywhere
Martin Almquist <malmquist@stanford.edu>
parents:
927
diff
changeset
|
326 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type) |
904
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
327 |
921
19c05eefc8c6
Add warning that interfaceNonConforming in LaplaceCurve only works for Cartesian grids.
Martin Almquist <malmquist@stanford.edu>
parents:
910
diff
changeset
|
328 % TODO: Make this work for curvilinear grids |
19c05eefc8c6
Add warning that interfaceNonConforming in LaplaceCurve only works for Cartesian grids.
Martin Almquist <malmquist@stanford.edu>
parents:
910
diff
changeset
|
329 warning('LaplaceCurvilinear: Non-conforming grid interpolation only works for Cartesian grids.'); |
19c05eefc8c6
Add warning that interfaceNonConforming in LaplaceCurve only works for Cartesian grids.
Martin Almquist <malmquist@stanford.edu>
parents:
910
diff
changeset
|
330 |
927
4291731570bb
Rename AWW OP. Make Interpolation operator classes take grid points as arguments. Remove LaplCurv.interpolationOperators. Introduce default struct in LaplCurv.interface.
Martin Almquist <malmquist@stanford.edu>
parents:
925
diff
changeset
|
331 % User can request special interpolation operators by specifying type.interpOpSet |
928
1c61d8fa9903
Replace opts by type everywhere
Martin Almquist <malmquist@stanford.edu>
parents:
927
diff
changeset
|
332 default_field(type, 'interpOpSet', @sbp.InterpOpsOP); |
1c61d8fa9903
Replace opts by type everywhere
Martin Almquist <malmquist@stanford.edu>
parents:
927
diff
changeset
|
333 interpOpSet = type.interpOpSet; |
1c61d8fa9903
Replace opts by type everywhere
Martin Almquist <malmquist@stanford.edu>
parents:
927
diff
changeset
|
334 tuning = type.tuning; |
904
14b093a344eb
Outline new interface method in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
743
diff
changeset
|
335 |
927
4291731570bb
Rename AWW OP. Make Interpolation operator classes take grid points as arguments. Remove LaplCurv.interpolationOperators. Introduce default struct in LaplCurv.interface.
Martin Almquist <malmquist@stanford.edu>
parents:
925
diff
changeset
|
336 |
906
0499239496cf
Remove property interpolation_type in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
904
diff
changeset
|
337 % u denotes the solution in the own domain |
0499239496cf
Remove property interpolation_type in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
904
diff
changeset
|
338 % v denotes the solution in the neighbour domain |
938
97291e1bd57c
Make LaplCurv send only number of grid points to InterpOpsXX. Remove coordinates from get_boundary_ops method.
Martin Almquist <malmquist@stanford.edu>
parents:
936
diff
changeset
|
339 [e_u, d_u, gamm_u, H_b_u, I_u] = obj.get_boundary_ops(boundary); |
97291e1bd57c
Make LaplCurv send only number of grid points to InterpOpsXX. Remove coordinates from get_boundary_ops method.
Martin Almquist <malmquist@stanford.edu>
parents:
936
diff
changeset
|
340 [e_v, d_v, gamm_v, H_b_v, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); |
720
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
341 |
938
97291e1bd57c
Make LaplCurv send only number of grid points to InterpOpsXX. Remove coordinates from get_boundary_ops method.
Martin Almquist <malmquist@stanford.edu>
parents:
936
diff
changeset
|
342 % Find the number of grid points along the interface |
97291e1bd57c
Make LaplCurv send only number of grid points to InterpOpsXX. Remove coordinates from get_boundary_ops method.
Martin Almquist <malmquist@stanford.edu>
parents:
936
diff
changeset
|
343 m_u = size(e_u, 2); |
97291e1bd57c
Make LaplCurv send only number of grid points to InterpOpsXX. Remove coordinates from get_boundary_ops method.
Martin Almquist <malmquist@stanford.edu>
parents:
936
diff
changeset
|
344 m_v = size(e_v, 2); |
922
1d91c2a8aada
Refactor LaplCurv.interfaceNonConf. Add method for interpolation operator generation in LaplCurv
Martin Almquist <malmquist@stanford.edu>
parents:
921
diff
changeset
|
345 |
720
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
346 Hi = obj.Hi; |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
347 a = obj.a; |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
348 |
5
5f6b0b6a012b
First try at interface implementation in WaveCurve2s
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
349 u = obj; |
5f6b0b6a012b
First try at interface implementation in WaveCurve2s
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
350 v = neighbour_scheme; |
5f6b0b6a012b
First try at interface implementation in WaveCurve2s
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
351 |
85
643bc513b8b8
Wave2DCurve: Fixed error in tau. Also made random guess at fixing theory from paper.
Jonatan Werpers <jonatan@werpers.com>
parents:
83
diff
changeset
|
352 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2; |
643bc513b8b8
Wave2DCurve: Fixed error in tau. Also made random guess at fixing theory from paper.
Jonatan Werpers <jonatan@werpers.com>
parents:
83
diff
changeset
|
353 b2_u = gamm_u*u.lambda(I_u)./u.a22(I_u).^2; |
643bc513b8b8
Wave2DCurve: Fixed error in tau. Also made random guess at fixing theory from paper.
Jonatan Werpers <jonatan@werpers.com>
parents:
83
diff
changeset
|
354 b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2; |
643bc513b8b8
Wave2DCurve: Fixed error in tau. Also made random guess at fixing theory from paper.
Jonatan Werpers <jonatan@werpers.com>
parents:
83
diff
changeset
|
355 b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2; |
5
5f6b0b6a012b
First try at interface implementation in WaveCurve2s
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
356 |
720
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
357 tau_u = -1./(4*b1_u) -1./(4*b2_u); |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
358 tau_v = -1./(4*b1_v) -1./(4*b2_v); |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
359 |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
360 tau_u = tuning * spdiag(tau_u); |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
361 tau_v = tuning * spdiag(tau_v); |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
362 beta_u = tau_v; |
0 | 363 |
922
1d91c2a8aada
Refactor LaplCurv.interfaceNonConf. Add method for interpolation operator generation in LaplCurv
Martin Almquist <malmquist@stanford.edu>
parents:
921
diff
changeset
|
364 % Build interpolation operators |
938
97291e1bd57c
Make LaplCurv send only number of grid points to InterpOpsXX. Remove coordinates from get_boundary_ops method.
Martin Almquist <malmquist@stanford.edu>
parents:
936
diff
changeset
|
365 intOps = interpOpSet(m_u, m_v, obj.order, neighbour_scheme.order); |
927
4291731570bb
Rename AWW OP. Make Interpolation operator classes take grid points as arguments. Remove LaplCurv.interpolationOperators. Introduce default struct in LaplCurv.interface.
Martin Almquist <malmquist@stanford.edu>
parents:
925
diff
changeset
|
366 Iu2v = intOps.Iu2v; |
4291731570bb
Rename AWW OP. Make Interpolation operator classes take grid points as arguments. Remove LaplCurv.interpolationOperators. Introduce default struct in LaplCurv.interface.
Martin Almquist <malmquist@stanford.edu>
parents:
925
diff
changeset
|
367 Iv2u = intOps.Iv2u; |
922
1d91c2a8aada
Refactor LaplCurv.interfaceNonConf. Add method for interpolation operator generation in LaplCurv
Martin Almquist <malmquist@stanford.edu>
parents:
921
diff
changeset
|
368 |
720
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
369 closure = a*Hi*e_u*tau_u*H_b_u*e_u' + ... |
927
4291731570bb
Rename AWW OP. Make Interpolation operator classes take grid points as arguments. Remove LaplCurv.interpolationOperators. Introduce default struct in LaplCurv.interface.
Martin Almquist <malmquist@stanford.edu>
parents:
925
diff
changeset
|
370 a*Hi*e_u*H_b_u*Iv2u.bad*beta_u*Iu2v.good*e_u' + ... |
720
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
371 a*1/2*Hi*d_u*H_b_u*e_u' + ... |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
372 -a*1/2*Hi*e_u*H_b_u*d_u'; |
0 | 373 |
927
4291731570bb
Rename AWW OP. Make Interpolation operator classes take grid points as arguments. Remove LaplCurv.interpolationOperators. Introduce default struct in LaplCurv.interface.
Martin Almquist <malmquist@stanford.edu>
parents:
925
diff
changeset
|
374 penalty = -a*Hi*e_u*tau_u*H_b_u*Iv2u.good*e_v' + ... |
4291731570bb
Rename AWW OP. Make Interpolation operator classes take grid points as arguments. Remove LaplCurv.interpolationOperators. Introduce default struct in LaplCurv.interface.
Martin Almquist <malmquist@stanford.edu>
parents:
925
diff
changeset
|
375 -a*Hi*e_u*H_b_u*Iv2u.bad*beta_u*e_v' + ... |
4291731570bb
Rename AWW OP. Make Interpolation operator classes take grid points as arguments. Remove LaplCurv.interpolationOperators. Introduce default struct in LaplCurv.interface.
Martin Almquist <malmquist@stanford.edu>
parents:
925
diff
changeset
|
376 -a*1/2*Hi*d_u*H_b_u*Iv2u.good*e_v' + ... |
4291731570bb
Rename AWW OP. Make Interpolation operator classes take grid points as arguments. Remove LaplCurv.interpolationOperators. Introduce default struct in LaplCurv.interface.
Martin Almquist <malmquist@stanford.edu>
parents:
925
diff
changeset
|
377 -a*1/2*Hi*e_u*H_b_u*Iv2u.bad*d_v'; |
0 | 378 |
379 end | |
380 | |
906
0499239496cf
Remove property interpolation_type in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
904
diff
changeset
|
381 % Returns the boundary ops and sign for the boundary specified by the string boundary. |
0 | 382 % The right boundary is considered the positive boundary |
85
643bc513b8b8
Wave2DCurve: Fixed error in tau. Also made random guess at fixing theory from paper.
Jonatan Werpers <jonatan@werpers.com>
parents:
83
diff
changeset
|
383 % |
906
0499239496cf
Remove property interpolation_type in LaplaceCurvilinear
Martin Almquist <malmquist@stanford.edu>
parents:
904
diff
changeset
|
384 % I -- the indices of the boundary points in the grid matrix |
938
97291e1bd57c
Make LaplCurv send only number of grid points to InterpOpsXX. Remove coordinates from get_boundary_ops method.
Martin Almquist <malmquist@stanford.edu>
parents:
936
diff
changeset
|
385 function [e, d, gamm, H_b, I] = get_boundary_ops(obj, boundary) |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
386 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m)); |
85
643bc513b8b8
Wave2DCurve: Fixed error in tau. Also made random guess at fixing theory from paper.
Jonatan Werpers <jonatan@werpers.com>
parents:
83
diff
changeset
|
387 |
0 | 388 switch boundary |
389 case 'w' | |
390 e = obj.e_w; | |
561
3a13916f8ff0
Switch to using boundary ops for normal derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
560
diff
changeset
|
391 d = obj.d_w; |
566
9c98a0526afc
Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents:
565
diff
changeset
|
392 H_b = obj.H_w; |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
393 I = ind(1,:); |
0 | 394 case 'e' |
395 e = obj.e_e; | |
561
3a13916f8ff0
Switch to using boundary ops for normal derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
560
diff
changeset
|
396 d = obj.d_e; |
566
9c98a0526afc
Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents:
565
diff
changeset
|
397 H_b = obj.H_e; |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
398 I = ind(end,:); |
0 | 399 case 's' |
400 e = obj.e_s; | |
561
3a13916f8ff0
Switch to using boundary ops for normal derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
560
diff
changeset
|
401 d = obj.d_s; |
566
9c98a0526afc
Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents:
565
diff
changeset
|
402 H_b = obj.H_s; |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
403 I = ind(:,1)'; |
0 | 404 case 'n' |
405 e = obj.e_n; | |
561
3a13916f8ff0
Switch to using boundary ops for normal derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
560
diff
changeset
|
406 d = obj.d_n; |
566
9c98a0526afc
Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents:
565
diff
changeset
|
407 H_b = obj.H_n; |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
408 I = ind(:,end)'; |
0 | 409 otherwise |
410 error('No such boundary: boundary = %s',boundary); | |
411 end | |
412 | |
413 switch boundary | |
414 case {'w','e'} | |
415 gamm = obj.gamm_u; | |
416 case {'s','n'} | |
417 gamm = obj.gamm_v; | |
418 end | |
922
1d91c2a8aada
Refactor LaplCurv.interfaceNonConf. Add method for interpolation operator generation in LaplCurv
Martin Almquist <malmquist@stanford.edu>
parents:
921
diff
changeset
|
419 end |
1d91c2a8aada
Refactor LaplCurv.interfaceNonConf. Add method for interpolation operator generation in LaplCurv
Martin Almquist <malmquist@stanford.edu>
parents:
921
diff
changeset
|
420 |
0 | 421 function N = size(obj) |
422 N = prod(obj.m); | |
423 end | |
424 end | |
566
9c98a0526afc
Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents:
565
diff
changeset
|
425 end |