Mercurial > repos > public > sbplib
annotate +scheme/LaplaceCurvilinear.m @ 534:55a7777dfcd0 feature/boundaryGroup
Implement boundary_condition for boundary groups
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Thu, 03 Aug 2017 18:07:27 +0200 |
| parents | 56eb7c088bd4 |
| children | 08b6281ba2a9 |
| 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 | |
| 10 D % non-stabalized scheme operator | |
| 11 M % Derivative norm | |
|
452
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
12 a,b |
| 0 | 13 J, Ji |
| 14 a11, a12, a22 | |
| 15 | |
| 16 H % Discrete norm | |
| 17 Hi | |
| 18 H_u, H_v % Norms in the x and y directions | |
| 19 Hu,Hv % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. | |
| 20 Hi_u, Hi_v | |
| 21 Hiu, Hiv | |
| 22 e_w, e_e, e_s, e_n | |
| 23 du_w, dv_w | |
| 24 du_e, dv_e | |
| 25 du_s, dv_s | |
| 26 du_n, dv_n | |
| 27 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
|
28 lambda |
|
384
32df00102268
Fix bug in Characteristic BC for Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
360
diff
changeset
|
29 |
|
449
0707a7192bc3
Add first derivative operators to Wave2DCurve
Jonatan Werpers <jonatan@werpers.com>
parents:
395
diff
changeset
|
30 Dx, Dy % Physical derivatives |
|
0707a7192bc3
Add first derivative operators to Wave2DCurve
Jonatan Werpers <jonatan@werpers.com>
parents:
395
diff
changeset
|
31 |
|
384
32df00102268
Fix bug in Characteristic BC for Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
360
diff
changeset
|
32 x_u |
|
32df00102268
Fix bug in Characteristic BC for Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
360
diff
changeset
|
33 x_v |
|
32df00102268
Fix bug in Characteristic BC for Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
360
diff
changeset
|
34 y_u |
|
32df00102268
Fix bug in Characteristic BC for Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
360
diff
changeset
|
35 y_v |
| 0 | 36 end |
| 37 | |
| 38 methods | |
|
452
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
39 % Implements a*div(b*grad(u)) as a SBP scheme |
|
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
40 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
|
41 default_arg('opSet',@sbp.D2Variable); |
|
452
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
42 default_arg('a', 1); |
|
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
43 default_arg('b', 1); |
|
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
44 |
|
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
45 |
|
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
46 if b ~=1 |
|
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
47 error('Not implemented yet') |
|
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
48 end |
| 0 | 49 |
|
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
50 assert(isa(g, 'grid.Curvilinear')) |
| 0 | 51 |
|
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
52 m = g.size(); |
| 0 | 53 m_u = m(1); |
| 54 m_v = m(2); | |
|
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
55 m_tot = g.N(); |
| 0 | 56 |
|
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
57 h = g.scaling(); |
|
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
58 h_u = h(1); |
|
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
59 h_v = h(2); |
| 0 | 60 |
| 61 % Operators | |
|
303
f18142c1530b
Fixed Wave2dCurve for new operators. Some cleanup in d2_variable_4.
Jonatan Werpers <jonatan@werpers.com>
parents:
96
diff
changeset
|
62 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
|
63 ops_v = opSet(m_v, {0, 1}, order); |
| 0 | 64 |
| 65 I_u = speye(m_u); | |
| 66 I_v = speye(m_v); | |
| 67 | |
|
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
68 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
|
69 D2_u = ops_u.D2; |
|
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
70 H_u = ops_u.H; |
|
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
71 Hi_u = ops_u.HI; |
|
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
72 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
|
73 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
|
74 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
|
75 d1_r_u = ops_u.d1_r; |
| 0 | 76 |
|
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
77 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
|
78 D2_v = ops_v.D2; |
|
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
79 H_v = ops_v.H; |
|
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
80 Hi_v = ops_v.HI; |
|
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
81 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
|
82 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
|
83 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
|
84 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
|
85 |
|
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
86 Du = kr(D1_u,I_v); |
|
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
87 Dv = kr(I_u,D1_v); |
| 0 | 88 |
| 89 % Metric derivatives | |
|
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
90 coords = g.points(); |
|
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
91 x = coords(:,1); |
|
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
92 y = coords(:,2); |
| 0 | 93 |
|
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
94 x_u = Du*x; |
|
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
95 x_v = Dv*x; |
|
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
96 y_u = Du*y; |
|
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
97 y_v = Dv*y; |
| 0 | 98 |
| 99 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
|
100 a11 = 1./J .* (x_v.^2 + y_v.^2); |
| 0 | 101 a12 = -1./J .* (x_u.*x_v + y_u.*y_v); |
| 102 a22 = 1./J .* (x_u.^2 + y_u.^2); | |
| 103 lambda = 1/2 * (a11 + a22 - sqrt((a11-a22).^2 + 4*a12.^2)); | |
| 104 | |
|
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
105 % Assemble full operators |
|
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
106 L_12 = spdiags(a12, 0, m_tot, m_tot); |
|
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
107 Duv = Du*L_12*Dv; |
|
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
108 Dvu = Dv*L_12*Du; |
| 0 | 109 |
| 110 Duu = sparse(m_tot); | |
| 111 Dvv = sparse(m_tot); | |
|
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
112 ind = grid.funcToMatrix(g, 1:m_tot); |
| 0 | 113 |
| 114 for i = 1:m_v | |
|
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
115 D = D2_u(a11(ind(:,i))); |
|
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
116 p = ind(:,i); |
| 0 | 117 Duu(p,p) = D; |
| 118 end | |
| 119 | |
| 120 for i = 1:m_u | |
|
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
121 D = D2_v(a22(ind(i,:))); |
|
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
122 p = ind(i,:); |
| 0 | 123 Dvv(p,p) = D; |
| 124 end | |
| 125 | |
| 126 obj.H = kr(H_u,H_v); | |
| 127 obj.Hi = kr(Hi_u,Hi_v); | |
| 128 obj.Hu = kr(H_u,I_v); | |
| 129 obj.Hv = kr(I_u,H_v); | |
| 130 obj.Hiu = kr(Hi_u,I_v); | |
| 131 obj.Hiv = kr(I_u,Hi_v); | |
| 132 | |
| 133 obj.e_w = kr(e_l_u,I_v); | |
| 134 obj.e_e = kr(e_r_u,I_v); | |
| 135 obj.e_s = kr(I_u,e_l_v); | |
| 136 obj.e_n = kr(I_u,e_r_v); | |
| 137 obj.du_w = kr(d1_l_u,I_v); | |
| 138 obj.dv_w = (obj.e_w'*Dv)'; | |
| 139 obj.du_e = kr(d1_r_u,I_v); | |
| 140 obj.dv_e = (obj.e_e'*Dv)'; | |
| 141 obj.du_s = (obj.e_s'*Du)'; | |
| 142 obj.dv_s = kr(I_u,d1_l_v); | |
| 143 obj.du_n = (obj.e_n'*Du)'; | |
| 370 | 144 obj.dv_n = kr(I_u,d1_r_v); |
| 376 | 145 |
|
384
32df00102268
Fix bug in Characteristic BC for Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
360
diff
changeset
|
146 obj.x_u = x_u; |
|
32df00102268
Fix bug in Characteristic BC for Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
360
diff
changeset
|
147 obj.x_v = x_v; |
|
32df00102268
Fix bug in Characteristic BC for Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
360
diff
changeset
|
148 obj.y_u = y_u; |
|
32df00102268
Fix bug in Characteristic BC for Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
360
diff
changeset
|
149 obj.y_v = y_v; |
|
32df00102268
Fix bug in Characteristic BC for Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
360
diff
changeset
|
150 |
| 0 | 151 obj.m = m; |
| 152 obj.h = [h_u h_v]; | |
| 153 obj.order = order; | |
|
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
154 obj.grid = g; |
| 0 | 155 |
|
452
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
156 obj.a = a; |
|
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
157 obj.b = b; |
|
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
158 obj.J = spdiags(J, 0, m_tot, m_tot); |
|
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
159 obj.Ji = spdiags(1./J, 0, m_tot, m_tot); |
| 0 | 160 obj.a11 = a11; |
| 161 obj.a12 = a12; | |
| 162 obj.a22 = a22; | |
|
452
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
163 obj.D = obj.Ji*a*(Duu + Duv + Dvu + Dvv); |
|
27
97a638f91fb8
Added function spdiag(). Fixed a bunch of bugs in the Wave2dCurve scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
5
diff
changeset
|
164 obj.lambda = lambda; |
| 0 | 165 |
|
449
0707a7192bc3
Add first derivative operators to Wave2DCurve
Jonatan Werpers <jonatan@werpers.com>
parents:
395
diff
changeset
|
166 obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv; |
|
0707a7192bc3
Add first derivative operators to Wave2DCurve
Jonatan Werpers <jonatan@werpers.com>
parents:
395
diff
changeset
|
167 obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv; |
|
0707a7192bc3
Add first derivative operators to Wave2DCurve
Jonatan Werpers <jonatan@werpers.com>
parents:
395
diff
changeset
|
168 |
|
389
42c89b5eedc0
Add borrowing constants for D2 operators in D4Variable
Jonatan Werpers <jonatan@werpers.com>
parents:
384
diff
changeset
|
169 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
|
170 obj.gamm_v = h_v*ops_v.borrowing.M.d1; |
| 0 | 171 end |
| 172 | |
| 173 | |
| 174 % Closure functions return the opertors applied to the own doamin to close the boundary | |
| 175 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. | |
| 176 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. | |
| 177 % type is a string specifying the type of boundary condition if there are several. | |
| 178 % data is a function returning the data that should be applied at the boundary. | |
| 179 % neighbour_scheme is an instance of Scheme that should be interfaced to. | |
| 180 % 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
|
181 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter) |
| 0 | 182 default_arg('type','neumann'); |
|
347
85c2fe06d551
Implemented characteristic form BC in Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
337
diff
changeset
|
183 default_arg('parameter', []); |
| 0 | 184 |
|
384
32df00102268
Fix bug in Characteristic BC for Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
360
diff
changeset
|
185 [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv , ~, ~, ~, scale_factor] = obj.get_boundary_ops(boundary); |
| 0 | 186 switch type |
| 187 % Dirichlet boundary condition | |
| 188 case {'D','d','dirichlet'} | |
|
82
df642750e8f7
Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
55
diff
changeset
|
189 % v denotes the solution in the neighbour domain |
|
df642750e8f7
Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
55
diff
changeset
|
190 tuning = 1.2; |
|
df642750e8f7
Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
55
diff
changeset
|
191 % tuning = 20.2; |
|
df642750e8f7
Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
55
diff
changeset
|
192 [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t] = obj.get_boundary_ops(boundary); |
|
df642750e8f7
Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
55
diff
changeset
|
193 |
|
df642750e8f7
Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
55
diff
changeset
|
194 a_n = spdiag(coeff_n); |
|
df642750e8f7
Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
55
diff
changeset
|
195 a_t = spdiag(coeff_t); |
|
df642750e8f7
Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
55
diff
changeset
|
196 |
|
df642750e8f7
Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
55
diff
changeset
|
197 F = (s * a_n * d_n' + s * a_t*d_t')'; |
|
df642750e8f7
Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
55
diff
changeset
|
198 |
|
df642750e8f7
Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
55
diff
changeset
|
199 u = obj; |
| 0 | 200 |
|
82
df642750e8f7
Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
55
diff
changeset
|
201 b1 = gamm*u.lambda./u.a11.^2; |
|
df642750e8f7
Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
55
diff
changeset
|
202 b2 = gamm*u.lambda./u.a22.^2; |
| 0 | 203 |
|
82
df642750e8f7
Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
55
diff
changeset
|
204 tau = -1./b1 - 1./b2; |
|
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
205 tau = tuning * spdiag(tau); |
|
360
447ceb41fb65
Fixed penalty parameter in Wave2Dcurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
347
diff
changeset
|
206 sig1 = 1; |
| 0 | 207 |
|
82
df642750e8f7
Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
55
diff
changeset
|
208 penalty_parameter_1 = halfnorm_inv_n*(tau + sig1*halfnorm_inv_t*F*e'*halfnorm_t)*e; |
| 0 | 209 |
|
452
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
210 closure = obj.Ji*obj.a * penalty_parameter_1*e'; |
|
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
211 penalty = -obj.Ji*obj.a * penalty_parameter_1; |
| 0 | 212 |
| 213 | |
| 214 % Neumann boundary condition | |
| 215 case {'N','n','neumann'} | |
| 216 a_n = spdiags(coeff_n,0,length(coeff_n),length(coeff_n)); | |
| 217 a_t = spdiags(coeff_t,0,length(coeff_t),length(coeff_t)); | |
| 218 d = (a_n * d_n' + a_t*d_t')'; | |
| 219 | |
| 220 tau1 = -s; | |
| 221 tau2 = 0; | |
|
452
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
222 tau = obj.a * obj.Ji*(tau1*e + tau2*d); |
| 0 | 223 |
| 224 closure = halfnorm_inv*tau*d'; | |
|
347
85c2fe06d551
Implemented characteristic form BC in Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
337
diff
changeset
|
225 penalty = -halfnorm_inv*tau; |
| 0 | 226 |
|
347
85c2fe06d551
Implemented characteristic form BC in Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
337
diff
changeset
|
227 % Characteristic boundary condition |
|
85c2fe06d551
Implemented characteristic form BC in Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
337
diff
changeset
|
228 case {'characteristic', 'char', 'c'} |
|
85c2fe06d551
Implemented characteristic form BC in Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
337
diff
changeset
|
229 default_arg('parameter', 1); |
|
85c2fe06d551
Implemented characteristic form BC in Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
337
diff
changeset
|
230 beta = parameter; |
| 0 | 231 |
|
347
85c2fe06d551
Implemented characteristic form BC in Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
337
diff
changeset
|
232 a_n = spdiags(coeff_n,0,length(coeff_n),length(coeff_n)); |
|
85c2fe06d551
Implemented characteristic form BC in Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
337
diff
changeset
|
233 a_t = spdiags(coeff_t,0,length(coeff_t),length(coeff_t)); |
|
85c2fe06d551
Implemented characteristic form BC in Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
337
diff
changeset
|
234 d = s*(a_n * d_n' + a_t*d_t')'; % outward facing normal derivative |
|
85c2fe06d551
Implemented characteristic form BC in Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
337
diff
changeset
|
235 |
|
452
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
236 tau = -obj.a * 1/beta*obj.Ji*e; |
|
347
85c2fe06d551
Implemented characteristic form BC in Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
337
diff
changeset
|
237 |
|
452
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
238 closure{1} = halfnorm_inv*tau*spdiag(scale_factor)*e'; |
|
347
85c2fe06d551
Implemented characteristic form BC in Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
337
diff
changeset
|
239 closure{2} = halfnorm_inv*tau*beta*d'; |
|
85c2fe06d551
Implemented characteristic form BC in Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
337
diff
changeset
|
240 penalty = -halfnorm_inv*tau; |
| 0 | 241 |
| 242 % Unknown, boundary condition | |
| 243 otherwise | |
| 244 error('No such boundary condition: type = %s',type); | |
| 245 end | |
| 246 end | |
| 247 | |
| 248 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) | |
| 249 % u denotes the solution in the own domain | |
| 250 % v denotes the solution in the neighbour domain | |
|
5
5f6b0b6a012b
First try at interface implementation in WaveCurve2s
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
251 tuning = 1.2; |
|
27
97a638f91fb8
Added function spdiag(). Fixed a bunch of bugs in the Wave2dCurve scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
5
diff
changeset
|
252 % tuning = 20.2; |
|
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
|
253 [e_u, d_n_u, d_t_u, coeff_n_u, coeff_t_u, s_u, gamm_u, halfnorm_inv_u_n, halfnorm_inv_u_t, halfnorm_u_t, I_u] = obj.get_boundary_ops(boundary); |
|
643bc513b8b8
Wave2DCurve: Fixed error in tau. Also made random guess at fixing theory from paper.
Jonatan Werpers <jonatan@werpers.com>
parents:
83
diff
changeset
|
254 [e_v, d_n_v, d_t_v, coeff_n_v, coeff_t_v, s_v, gamm_v, halfnorm_inv_v_n, halfnorm_inv_v_t, halfnorm_v_t, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); |
| 0 | 255 |
|
27
97a638f91fb8
Added function spdiag(). Fixed a bunch of bugs in the Wave2dCurve scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
5
diff
changeset
|
256 a_n_u = spdiag(coeff_n_u); |
|
97a638f91fb8
Added function spdiag(). Fixed a bunch of bugs in the Wave2dCurve scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
5
diff
changeset
|
257 a_t_u = spdiag(coeff_t_u); |
|
97a638f91fb8
Added function spdiag(). Fixed a bunch of bugs in the Wave2dCurve scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
5
diff
changeset
|
258 a_n_v = spdiag(coeff_n_v); |
|
97a638f91fb8
Added function spdiag(). Fixed a bunch of bugs in the Wave2dCurve scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
5
diff
changeset
|
259 a_t_v = spdiag(coeff_t_v); |
|
97a638f91fb8
Added function spdiag(). Fixed a bunch of bugs in the Wave2dCurve scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
5
diff
changeset
|
260 |
|
97a638f91fb8
Added function spdiag(). Fixed a bunch of bugs in the Wave2dCurve scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
5
diff
changeset
|
261 F_u = (s_u * a_n_u * d_n_u' + s_u * a_t_u*d_t_u')'; |
|
97a638f91fb8
Added function spdiag(). Fixed a bunch of bugs in the Wave2dCurve scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
5
diff
changeset
|
262 F_v = (s_v * a_n_v * d_n_v' + s_v * a_t_v*d_t_v')'; |
| 0 | 263 |
|
5
5f6b0b6a012b
First try at interface implementation in WaveCurve2s
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
264 u = obj; |
|
5f6b0b6a012b
First try at interface implementation in WaveCurve2s
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
265 v = neighbour_scheme; |
|
5f6b0b6a012b
First try at interface implementation in WaveCurve2s
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
266 |
|
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
|
267 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
|
268 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
|
269 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
|
270 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
|
271 |
|
82
df642750e8f7
Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
55
diff
changeset
|
272 tau = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v); |
|
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
273 tau = tuning * spdiag(tau); |
| 0 | 274 sig1 = 1/2; |
| 83 | 275 sig2 = -1/2; |
| 0 | 276 |
|
96
19d0c9325a3e
Wave2dCurve: Removed some comments.
Jonatan Werpers <jonatan@werpers.com>
parents:
89
diff
changeset
|
277 penalty_parameter_1 = halfnorm_inv_u_n*(e_u*tau + sig1*halfnorm_inv_u_t*F_u*e_u'*halfnorm_u_t*e_u); |
| 0 | 278 penalty_parameter_2 = halfnorm_inv_u_n * sig2 * e_u; |
| 279 | |
| 280 | |
|
452
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
281 closure = obj.Ji*obj.a * ( penalty_parameter_1*e_u' + penalty_parameter_2*F_u'); |
|
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
282 penalty = obj.Ji*obj.a * (-penalty_parameter_1*e_v' + penalty_parameter_2*F_v'); |
| 0 | 283 end |
| 284 | |
| 285 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. | |
| 286 % 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
|
287 % |
|
643bc513b8b8
Wave2DCurve: Fixed error in tau. Also made random guess at fixing theory from paper.
Jonatan Werpers <jonatan@werpers.com>
parents:
83
diff
changeset
|
288 % I -- the indecies of the boundary points in the grid matrix |
|
384
32df00102268
Fix bug in Characteristic BC for Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
360
diff
changeset
|
289 function [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t, I, scale_factor] = get_boundary_ops(obj, 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
|
290 |
|
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
291 % gridMatrix = zeros(obj.m(2),obj.m(1)); |
|
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
292 % gridMatrix(:) = 1:numel(gridMatrix); |
|
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
293 |
|
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
294 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
|
295 |
| 0 | 296 switch boundary |
| 297 case 'w' | |
| 298 e = obj.e_w; | |
| 299 d_n = obj.du_w; | |
| 300 d_t = obj.dv_w; | |
| 301 s = -1; | |
| 302 | |
|
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
303 I = ind(1,:); |
|
89
53fe4b64f65e
Minor cleaning of Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
85
diff
changeset
|
304 coeff_n = obj.a11(I); |
|
53fe4b64f65e
Minor cleaning of Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
85
diff
changeset
|
305 coeff_t = obj.a12(I); |
|
384
32df00102268
Fix bug in Characteristic BC for Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
360
diff
changeset
|
306 scale_factor = sqrt(obj.x_v(I).^2 + obj.y_v(I).^2); |
| 0 | 307 case 'e' |
| 308 e = obj.e_e; | |
| 309 d_n = obj.du_e; | |
| 310 d_t = obj.dv_e; | |
| 311 s = 1; | |
| 312 | |
|
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
313 I = ind(end,:); |
|
89
53fe4b64f65e
Minor cleaning of Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
85
diff
changeset
|
314 coeff_n = obj.a11(I); |
|
53fe4b64f65e
Minor cleaning of Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
85
diff
changeset
|
315 coeff_t = obj.a12(I); |
|
384
32df00102268
Fix bug in Characteristic BC for Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
360
diff
changeset
|
316 scale_factor = sqrt(obj.x_v(I).^2 + obj.y_v(I).^2); |
| 0 | 317 case 's' |
| 318 e = obj.e_s; | |
| 319 d_n = obj.dv_s; | |
| 320 d_t = obj.du_s; | |
| 321 s = -1; | |
| 322 | |
|
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
323 I = ind(:,1)'; |
|
89
53fe4b64f65e
Minor cleaning of Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
85
diff
changeset
|
324 coeff_n = obj.a22(I); |
|
53fe4b64f65e
Minor cleaning of Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
85
diff
changeset
|
325 coeff_t = obj.a12(I); |
|
384
32df00102268
Fix bug in Characteristic BC for Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
360
diff
changeset
|
326 scale_factor = sqrt(obj.x_u(I).^2 + obj.y_u(I).^2); |
| 0 | 327 case 'n' |
| 328 e = obj.e_n; | |
| 329 d_n = obj.dv_n; | |
| 330 d_t = obj.du_n; | |
| 331 s = 1; | |
| 332 | |
|
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
333 I = ind(:,end)'; |
|
89
53fe4b64f65e
Minor cleaning of Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
85
diff
changeset
|
334 coeff_n = obj.a22(I); |
|
53fe4b64f65e
Minor cleaning of Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
85
diff
changeset
|
335 coeff_t = obj.a12(I); |
|
384
32df00102268
Fix bug in Characteristic BC for Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
360
diff
changeset
|
336 scale_factor = sqrt(obj.x_u(I).^2 + obj.y_u(I).^2); |
| 0 | 337 otherwise |
| 338 error('No such boundary: boundary = %s',boundary); | |
| 339 end | |
| 340 | |
| 341 switch boundary | |
| 342 case {'w','e'} | |
| 343 halfnorm_inv_n = obj.Hiu; | |
| 344 halfnorm_inv_t = obj.Hiv; | |
| 345 halfnorm_t = obj.Hv; | |
| 346 gamm = obj.gamm_u; | |
| 347 case {'s','n'} | |
| 348 halfnorm_inv_n = obj.Hiv; | |
| 349 halfnorm_inv_t = obj.Hiu; | |
| 350 halfnorm_t = obj.Hu; | |
| 351 gamm = obj.gamm_v; | |
| 352 end | |
| 353 end | |
| 354 | |
| 355 function N = size(obj) | |
| 356 N = prod(obj.m); | |
| 357 end | |
| 358 | |
| 359 | |
| 360 end | |
| 361 end |