|
0
|
1 classdef Wave2dCurve < scheme.Scheme
|
|
|
2 properties
|
|
|
3 m % Number of points in each direction, possibly a vector
|
|
|
4 h % Grid spacing
|
|
|
5 u,v % Grid
|
|
|
6 x,y % Values of x and y for each grid point
|
|
|
7 X,Y % Grid point locations as matrices
|
|
|
8 order % Order accuracy for the approximation
|
|
|
9
|
|
|
10 D % non-stabalized scheme operator
|
|
|
11 M % Derivative norm
|
|
|
12 c
|
|
|
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
|
|
|
28 end
|
|
|
29
|
|
|
30 methods
|
|
|
31 function obj = Wave2dCurve(m,ti,order,c,opSet)
|
|
|
32 default_arg('opSet',@sbp.Variable);
|
|
|
33
|
|
|
34 if length(m) == 1
|
|
|
35 m = [m m];
|
|
|
36 end
|
|
|
37
|
|
|
38 m_u = m(1);
|
|
|
39 m_v = m(2);
|
|
|
40 m_tot = m_u*m_v;
|
|
|
41
|
|
|
42 [u, h_u] = util.get_grid(0, 1, m_u);
|
|
|
43 [v, h_v] = util.get_grid(0, 1, m_v);
|
|
|
44
|
|
|
45
|
|
|
46 % Operators
|
|
|
47 ops_u = opSet(m_u,h_u,order);
|
|
|
48 ops_v = opSet(m_v,h_v,order);
|
|
|
49
|
|
|
50 I_u = speye(m_u);
|
|
|
51 I_v = speye(m_v);
|
|
|
52
|
|
|
53 D1_u = sparse(ops_u.derivatives.D1);
|
|
|
54 D2_u = ops_u.derivatives.D2;
|
|
|
55 H_u = sparse(ops_u.norms.H);
|
|
|
56 Hi_u = sparse(ops_u.norms.HI);
|
|
|
57 % M_u = sparse(ops_u.norms.M);
|
|
|
58 e_l_u = sparse(ops_u.boundary.e_1);
|
|
|
59 e_r_u = sparse(ops_u.boundary.e_m);
|
|
|
60 d1_l_u = sparse(ops_u.boundary.S_1);
|
|
|
61 d1_r_u = sparse(ops_u.boundary.S_m);
|
|
|
62
|
|
|
63 D1_v = sparse(ops_v.derivatives.D1);
|
|
|
64 D2_v = ops_v.derivatives.D2;
|
|
|
65 H_v = sparse(ops_v.norms.H);
|
|
|
66 Hi_v = sparse(ops_v.norms.HI);
|
|
|
67 % M_v = sparse(ops_v.norms.M);
|
|
|
68 e_l_v = sparse(ops_v.boundary.e_1);
|
|
|
69 e_r_v = sparse(ops_v.boundary.e_m);
|
|
|
70 d1_l_v = sparse(ops_v.boundary.S_1);
|
|
|
71 d1_r_v = sparse(ops_v.boundary.S_m);
|
|
|
72
|
|
|
73
|
|
|
74 % Metric derivatives
|
|
|
75 [X,Y] = ti.map(u,v);
|
|
|
76
|
|
|
77 [x_u,x_v] = gridDerivatives(X,D1_u,D1_v);
|
|
|
78 [y_u,y_v] = gridDerivatives(Y,D1_u,D1_v);
|
|
|
79
|
|
|
80
|
|
|
81
|
|
|
82 J = x_u.*y_v - x_v.*y_u;
|
|
|
83 a11 = 1./J .* (x_v.^2 + y_v.^2); %% GÖR SOM MATRISER
|
|
|
84 a12 = -1./J .* (x_u.*x_v + y_u.*y_v);
|
|
|
85 a22 = 1./J .* (x_u.^2 + y_u.^2);
|
|
|
86 lambda = 1/2 * (a11 + a22 - sqrt((a11-a22).^2 + 4*a12.^2));
|
|
|
87
|
|
|
88 dof_order = reshape(1:m_u*m_v,m_v,m_u);
|
|
|
89
|
|
|
90 Duu = sparse(m_tot);
|
|
|
91 Dvv = sparse(m_tot);
|
|
|
92
|
|
|
93 for i = 1:m_v
|
|
|
94 D = D2_u(a11(i,:));
|
|
|
95 p = dof_order(i,:);
|
|
|
96 Duu(p,p) = D;
|
|
|
97 end
|
|
|
98
|
|
|
99 for i = 1:m_u
|
|
|
100 D = D2_v(a22(:,i));
|
|
|
101 p = dof_order(:,i);
|
|
|
102 Dvv(p,p) = D;
|
|
|
103 end
|
|
|
104
|
|
|
105 L_12 = spdiags(a12(:),0,m_tot,m_tot);
|
|
|
106 Du = kr(D1_u,I_v);
|
|
|
107 Dv = kr(I_u,D1_v);
|
|
|
108
|
|
|
109 Duv = Du*L_12*Dv;
|
|
|
110 Dvu = Dv*L_12*Du;
|
|
|
111
|
|
|
112
|
|
|
113
|
|
|
114 obj.H = kr(H_u,H_v);
|
|
|
115 obj.Hi = kr(Hi_u,Hi_v);
|
|
|
116 obj.Hu = kr(H_u,I_v);
|
|
|
117 obj.Hv = kr(I_u,H_v);
|
|
|
118 obj.Hiu = kr(Hi_u,I_v);
|
|
|
119 obj.Hiv = kr(I_u,Hi_v);
|
|
|
120
|
|
|
121 % obj.M = kr(M_u,H_v)+kr(H_u,M_v);
|
|
|
122 obj.e_w = kr(e_l_u,I_v);
|
|
|
123 obj.e_e = kr(e_r_u,I_v);
|
|
|
124 obj.e_s = kr(I_u,e_l_v);
|
|
|
125 obj.e_n = kr(I_u,e_r_v);
|
|
|
126 obj.du_w = kr(d1_l_u,I_v);
|
|
|
127 obj.dv_w = (obj.e_w'*Dv)';
|
|
|
128 obj.du_e = kr(d1_r_u,I_v);
|
|
|
129 obj.dv_e = (obj.e_e'*Dv)';
|
|
|
130 obj.du_s = (obj.e_s'*Du)';
|
|
|
131 obj.dv_s = kr(I_u,d1_l_v);
|
|
|
132 obj.du_n = (obj.e_n'*Du)';
|
|
|
133 obj.dv_n = kr(I_u,d1_r_v);
|
|
|
134
|
|
|
135 obj.m = m;
|
|
|
136 obj.h = [h_u h_v];
|
|
|
137 obj.order = order;
|
|
|
138
|
|
|
139
|
|
|
140 obj.c = c;
|
|
|
141 obj.J = spdiags(J(:),0,m_tot,m_tot);
|
|
|
142 obj.Ji = spdiags(1./J(:),0,m_tot,m_tot);
|
|
|
143 obj.a11 = a11;
|
|
|
144 obj.a12 = a12;
|
|
|
145 obj.a22 = a22;
|
|
|
146 obj.D = obj.Ji*c^2*(Duu + Duv + Dvu + Dvv);
|
|
|
147 obj.u = u;
|
|
|
148 obj.v = v;
|
|
|
149 obj.X = X;
|
|
|
150 obj.Y = Y;
|
|
|
151 obj.x = X(:);
|
|
|
152 obj.y = Y(:);
|
|
|
153
|
|
|
154 obj.gamm_u = h_u*ops_u.borrowing.M.S;
|
|
|
155 obj.gamm_v = h_v*ops_v.borrowing.M.S;
|
|
|
156 end
|
|
|
157
|
|
|
158
|
|
|
159 % Closure functions return the opertors applied to the own doamin to close the boundary
|
|
|
160 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
|
|
|
161 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
|
|
|
162 % type is a string specifying the type of boundary condition if there are several.
|
|
|
163 % data is a function returning the data that should be applied at the boundary.
|
|
|
164 % neighbour_scheme is an instance of Scheme that should be interfaced to.
|
|
|
165 % neighbour_boundary is a string specifying which boundary to interface to.
|
|
|
166 function [closure, penalty] = boundary_condition(obj,boundary,type,data)
|
|
|
167 default_arg('type','neumann');
|
|
|
168 default_arg('data',0);
|
|
|
169
|
|
|
170 [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv] = obj.get_boundary_ops(boundary);
|
|
|
171
|
|
|
172 switch type
|
|
|
173 % Dirichlet boundary condition
|
|
|
174 case {'D','d','dirichlet'}
|
|
|
175 error('not implemented')
|
|
|
176 alpha = obj.alpha;
|
|
|
177
|
|
|
178 % tau1 < -alpha^2/gamma
|
|
|
179 tuning = 1.1;
|
|
|
180 tau1 = -tuning*alpha/gamm;
|
|
|
181 tau2 = s*alpha;
|
|
|
182
|
|
|
183 p = tau1*e + tau2*d;
|
|
|
184
|
|
|
185 closure = halfnorm_inv*p*e';
|
|
|
186
|
|
|
187 pp = halfnorm_inv*p;
|
|
|
188 switch class(data)
|
|
|
189 case 'double'
|
|
|
190 penalty = pp*data;
|
|
|
191 case 'function_handle'
|
|
|
192 penalty = @(t)pp*data(t);
|
|
|
193 otherwise
|
|
|
194 error('Weird data argument!')
|
|
|
195 end
|
|
|
196
|
|
|
197
|
|
|
198 % Neumann boundary condition
|
|
|
199 case {'N','n','neumann'}
|
|
|
200 c = obj.c;
|
|
|
201
|
|
|
202
|
|
|
203 a_n = spdiags(coeff_n,0,length(coeff_n),length(coeff_n));
|
|
|
204 a_t = spdiags(coeff_t,0,length(coeff_t),length(coeff_t));
|
|
|
205 d = (a_n * d_n' + a_t*d_t')';
|
|
|
206
|
|
|
207 tau1 = -s;
|
|
|
208 tau2 = 0;
|
|
|
209 tau = c.^2 * obj.Ji*(tau1*e + tau2*d);
|
|
|
210
|
|
|
211 closure = halfnorm_inv*tau*d';
|
|
|
212
|
|
|
213 pp = halfnorm_inv*tau;
|
|
|
214 switch class(data)
|
|
|
215 case 'double'
|
|
|
216 penalty = pp*data;
|
|
|
217 case 'function_handle'
|
|
|
218 penalty = @(t)pp*data(t);
|
|
|
219 otherwise
|
|
|
220 error('Weird data argument!')
|
|
|
221 end
|
|
|
222
|
|
|
223 % Unknown, boundary condition
|
|
|
224 otherwise
|
|
|
225 error('No such boundary condition: type = %s',type);
|
|
|
226 end
|
|
|
227 end
|
|
|
228
|
|
|
229 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
|
|
|
230 % u denotes the solution in the own domain
|
|
|
231 % v denotes the solution in the neighbour domain
|
|
|
232 [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] = obj.get_boundary_ops(boundary);
|
|
|
233 [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] = neighbour_scheme.get_boundary_ops(boundary);
|
|
|
234
|
|
|
235 F_u = s_u * a_n_u * d_n_u' + s_u * a_t_u*d_t_u';
|
|
|
236 F_v = s_v * a_n_v * d_n_v' + s_v * a_t_v*d_t_v';
|
|
|
237
|
|
|
238 tau = 111111111111111111111111111111;
|
|
|
239 sig1 = 1/2;
|
|
|
240 sig2 = -1/2;
|
|
|
241
|
|
|
242 penalty_parameter_1 = s_u*halfnorm_inv_u_n*(tau + sig1*halfnorm_inv_u_t*F_u'*halfnorm_u_t)*e_u;
|
|
|
243 penalty_parameter_2 = halfnorm_inv_u_n * sig2 * e_u;
|
|
|
244
|
|
|
245 tuning = 1.2;
|
|
|
246
|
|
|
247 alpha_u = obj.alpha;
|
|
|
248 alpha_v = neighbour_scheme.alpha;
|
|
|
249
|
|
|
250 % tau1 < -(alpha_u/gamm_u + alpha_v/gamm_v)
|
|
|
251
|
|
|
252 tau1 = -(alpha_u/gamm_u + alpha_v/gamm_v) * tuning;
|
|
|
253 tau2 = s_u*1/2*alpha_u;
|
|
|
254 sig1 = s_u*(-1/2);
|
|
|
255 sig2 = 0;
|
|
|
256
|
|
|
257 tau = tau1*e_u + tau2*d_u;
|
|
|
258 sig = sig1*e_u + sig2*d_u;
|
|
|
259
|
|
|
260 closure = halfnorm_inv*( tau*e_u' + sig*alpha_u*d_u');
|
|
|
261 penalty = halfnorm_inv*(-tau*e_v' - sig*alpha_v*d_v');
|
|
|
262 end
|
|
|
263
|
|
|
264 % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
|
|
|
265 % The right boundary is considered the positive boundary
|
|
|
266 function [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t] = get_boundary_ops(obj,boundary)
|
|
|
267 switch boundary
|
|
|
268 case 'w'
|
|
|
269 e = obj.e_w;
|
|
|
270 d_n = obj.du_w;
|
|
|
271 d_t = obj.dv_w;
|
|
|
272 s = -1;
|
|
|
273
|
|
|
274 coeff_n = obj.a11(:,1);
|
|
|
275 coeff_t = obj.a12(:,1);
|
|
|
276 case 'e'
|
|
|
277 e = obj.e_e;
|
|
|
278 d_n = obj.du_e;
|
|
|
279 d_t = obj.dv_e;
|
|
|
280 s = 1;
|
|
|
281
|
|
|
282 coeff_n = obj.a11(:,end);
|
|
|
283 coeff_t = obj.a12(:,end);
|
|
|
284 case 's'
|
|
|
285 e = obj.e_s;
|
|
|
286 d_n = obj.dv_s;
|
|
|
287 d_t = obj.du_s;
|
|
|
288 s = -1;
|
|
|
289
|
|
|
290 coeff_n = obj.a22(1,:)';
|
|
|
291 coeff_t = obj.a12(1,:)';
|
|
|
292 case 'n'
|
|
|
293 e = obj.e_n;
|
|
|
294 d_n = obj.dv_n;
|
|
|
295 d_t = obj.du_n;
|
|
|
296 s = 1;
|
|
|
297
|
|
|
298 coeff_n = obj.a22(end,:)';
|
|
|
299 coeff_t = obj.a12(end,:)';
|
|
|
300 otherwise
|
|
|
301 error('No such boundary: boundary = %s',boundary);
|
|
|
302 end
|
|
|
303
|
|
|
304 switch boundary
|
|
|
305 case {'w','e'}
|
|
|
306 halfnorm_inv_n = obj.Hiu;
|
|
|
307 halfnorm_inv_t = obj.Hiv;
|
|
|
308 halfnorm_t = obj.Hv;
|
|
|
309 gamm = obj.gamm_u;
|
|
|
310 case {'s','n'}
|
|
|
311 halfnorm_inv_n = obj.Hiv;
|
|
|
312 halfnorm_inv_t = obj.Hiu;
|
|
|
313 halfnorm_t = obj.Hu;
|
|
|
314 gamm = obj.gamm_v;
|
|
|
315 end
|
|
|
316 end
|
|
|
317
|
|
|
318 function N = size(obj)
|
|
|
319 N = prod(obj.m);
|
|
|
320 end
|
|
|
321
|
|
|
322 end
|
|
|
323
|
|
|
324 methods(Static)
|
|
|
325 % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u
|
|
|
326 % and bound_v of scheme schm_v.
|
|
|
327 % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l')
|
|
|
328 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v)
|
|
|
329 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v);
|
|
|
330 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
|
|
|
331 end
|
|
|
332 end
|
|
|
333 end |
