Mercurial > repos > public > sbplib
comparison +scheme/Wave2dCurve.m @ 27:97a638f91fb8
Added function spdiag(). Fixed a bunch of bugs in the Wave2dCurve scheme.
author | Jonatan Werpers <jonatan@werpers.com> |
---|---|
date | Sun, 27 Sep 2015 22:15:37 +0200 |
parents | 5f6b0b6a012b |
children | a8ed986fcf57 |
comparison
equal
deleted
inserted
replaced
26:ed6a704b028d | 27:97a638f91fb8 |
---|---|
23 du_w, dv_w | 23 du_w, dv_w |
24 du_e, dv_e | 24 du_e, dv_e |
25 du_s, dv_s | 25 du_s, dv_s |
26 du_n, dv_n | 26 du_n, dv_n |
27 gamm_u, gamm_v | 27 gamm_u, gamm_v |
28 lambda | |
28 end | 29 end |
29 | 30 |
30 methods | 31 methods |
31 function obj = Wave2dCurve(m,ti,order,c,opSet) | 32 function obj = Wave2dCurve(m,ti,order,c,opSet) |
32 default_arg('opSet',@sbp.Variable); | 33 default_arg('opSet',@sbp.Variable); |
34 default_arg('c', 1); | |
33 | 35 |
34 if length(m) == 1 | 36 if length(m) == 1 |
35 m = [m m]; | 37 m = [m m]; |
36 end | 38 end |
37 | 39 |
148 obj.v = v; | 150 obj.v = v; |
149 obj.X = X; | 151 obj.X = X; |
150 obj.Y = Y; | 152 obj.Y = Y; |
151 obj.x = X(:); | 153 obj.x = X(:); |
152 obj.y = Y(:); | 154 obj.y = Y(:); |
155 obj.lambda = lambda; | |
153 | 156 |
154 obj.gamm_u = h_u*ops_u.borrowing.M.S; | 157 obj.gamm_u = h_u*ops_u.borrowing.M.S; |
155 obj.gamm_v = h_v*ops_v.borrowing.M.S; | 158 obj.gamm_v = h_v*ops_v.borrowing.M.S; |
156 end | 159 end |
157 | 160 |
228 | 231 |
229 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) | 232 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) |
230 % u denotes the solution in the own domain | 233 % u denotes the solution in the own domain |
231 % v denotes the solution in the neighbour domain | 234 % v denotes the solution in the neighbour domain |
232 tuning = 1.2; | 235 tuning = 1.2; |
233 [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); | 236 % tuning = 20.2; |
234 [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); | 237 [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, a_n_u, a_t_u] = obj.get_boundary_ops(boundary); |
235 | 238 [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, a_n_v, a_t_v] = neighbour_scheme.get_boundary_ops(boundary); |
236 F_u = s_u * a_n_u * d_n_u' + s_u * a_t_u*d_t_u'; | 239 |
237 F_v = s_v * a_n_v * d_n_v' + s_v * a_t_v*d_t_v'; | 240 a_n_u = spdiag(coeff_n_u); |
241 a_t_u = spdiag(coeff_t_u); | |
242 a_n_v = spdiag(coeff_n_v); | |
243 a_t_v = spdiag(coeff_t_v); | |
244 | |
245 F_u = (s_u * a_n_u * d_n_u' + s_u * a_t_u*d_t_u')'; | |
246 F_v = (s_v * a_n_v * d_n_v' + s_v * a_t_v*d_t_v')'; | |
238 | 247 |
239 u = obj; | 248 u = obj; |
240 v = neighbour_scheme; | 249 v = neighbour_scheme; |
241 | 250 |
242 b1_u = gamm_u*u.lambda./u.a11.^2; | 251 b1_u = gamm_u*u.lambda./u.a11.^2; |
243 b2_u = gamm_u*u.lambda./u.a22.^2; | 252 b2_u = gamm_u*u.lambda./u.a22.^2; |
244 b1_v = gamm_v*v.lambda./v.a11.^2; | 253 b1_v = gamm_v*v.lambda./v.a11.^2; |
245 b2_v = gamm_v*v.lambda./v.a22.^2; | 254 b2_v = gamm_v*v.lambda./v.a22.^2; |
246 | 255 |
247 | 256 |
248 tau = -1./(4*b1_u) -1/(4*b1_v) -1/(4*b2_u) -1/(4*b2_v); | 257 tau = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v); |
249 m_tot = obj.m(1)*obj.m(2); | 258 m_tot = obj.m(1)*obj.m(2); |
250 tau = tuning * spdiags(tau(:),0,m_tot,m_tot); | 259 tau = tuning * spdiag(tau(:)); |
251 sig1 = 1/2; | 260 sig1 = 1/2; |
252 sig2 = -1/2; | 261 sig2 = -1/2*s_u; |
253 | 262 |
254 penalty_parameter_1 = s_u*halfnorm_inv_u_n*(tau + sig1*halfnorm_inv_u_t*F_u'*halfnorm_u_t)*e_u; | 263 penalty_parameter_1 = halfnorm_inv_u_n*(tau + sig1*halfnorm_inv_u_t*F_u*e_u'*halfnorm_u_t)*e_u; |
255 penalty_parameter_2 = halfnorm_inv_u_n * sig2 * e_u; | 264 penalty_parameter_2 = halfnorm_inv_u_n * sig2 * e_u; |
256 | 265 |
257 | 266 |
258 closure = obj.Ji*obj.c^2 * ( penalty_parameter_1*e_u' + penalty_parameter_2*F_u'); | 267 closure = obj.Ji*obj.c^2 * ( penalty_parameter_1*e_u' + penalty_parameter_2*F_u'); |
259 penalty = obj.Ji*obj.c^2 * (-penalty_parameter_1*e_v' - penalty_parameter_2*F_v'); | 268 penalty = obj.Ji*obj.c^2 * (-penalty_parameter_1*e_v' - penalty_parameter_2*F_v'); |
260 end | 269 end |
261 | 270 |
262 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. | 271 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. |
263 % The right boundary is considered the positive boundary | 272 % The right boundary is considered the positive boundary |
264 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) | 273 function [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,a_n, a_t] = get_boundary_ops(obj,boundary) |
265 switch boundary | 274 switch boundary |
266 case 'w' | 275 case 'w' |
267 e = obj.e_w; | 276 e = obj.e_w; |
268 d_n = obj.du_w; | 277 d_n = obj.du_w; |
269 d_t = obj.dv_w; | 278 d_t = obj.dv_w; |
303 case {'w','e'} | 312 case {'w','e'} |
304 halfnorm_inv_n = obj.Hiu; | 313 halfnorm_inv_n = obj.Hiu; |
305 halfnorm_inv_t = obj.Hiv; | 314 halfnorm_inv_t = obj.Hiv; |
306 halfnorm_t = obj.Hv; | 315 halfnorm_t = obj.Hv; |
307 gamm = obj.gamm_u; | 316 gamm = obj.gamm_u; |
317 a_n = obj.a11; | |
318 a_t = obj.a12; | |
308 case {'s','n'} | 319 case {'s','n'} |
309 halfnorm_inv_n = obj.Hiv; | 320 halfnorm_inv_n = obj.Hiv; |
310 halfnorm_inv_t = obj.Hiu; | 321 halfnorm_inv_t = obj.Hiu; |
311 halfnorm_t = obj.Hu; | 322 halfnorm_t = obj.Hu; |
312 gamm = obj.gamm_v; | 323 gamm = obj.gamm_v; |
324 a_n = obj.a22; | |
325 a_t = obj.a12; | |
313 end | 326 end |
314 end | 327 end |
315 | 328 |
316 function N = size(obj) | 329 function N = size(obj) |
317 N = prod(obj.m); | 330 N = prod(obj.m); |