Mercurial > repos > public > sbplib

comparison +scheme/Wave2dCurve.m @ 5:5f6b0b6a012b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
First try at interface implementation in WaveCurve2s
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 21 Sep 2015 09:42:15 +0200
parents 48b6fb693025
children 97a638f91fb8
comparison
equal deleted inserted replaced
4:8e14b5a577a6 5:5f6b0b6a012b
227 end 227 end
228 228
229 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 229 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
230 % u denotes the solution in the own domain 230 % u denotes the solution in the own domain
231 % v denotes the solution in the neighbour domain 231 % v denotes the solution in the neighbour domain
232 tuning = 1.2;
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_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 [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
235 F_u = s_u * a_n_u * d_n_u' + s_u * a_t_u*d_t_u'; 236 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 F_v = s_v * a_n_v * d_n_v' + s_v * a_t_v*d_t_v';
237 238
238 tau = 111111111111111111111111111111; 239 u = obj;
240 v = neighbour_scheme;
241
242 b1_u = gamm_u*u.lambda./u.a11.^2;
243 b2_u = gamm_u*u.lambda./u.a22.^2;
244 b1_v = gamm_v*v.lambda./v.a11.^2;
245 b2_v = gamm_v*v.lambda./v.a22.^2;
246
247
248 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);
250 tau = tuning * spdiags(tau(:),0,m_tot,m_tot);
239 sig1 = 1/2; 251 sig1 = 1/2;
240 sig2 = -1/2; 252 sig2 = -1/2;
241 253
242 penalty_parameter_1 = s_u*halfnorm_inv_u_n*(tau + sig1*halfnorm_inv_u_t*F_u'*halfnorm_u_t)*e_u; 254 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; 255 penalty_parameter_2 = halfnorm_inv_u_n * sig2 * e_u;
244 256
245 tuning = 1.2; 257
246 258 closure = obj.Ji*obj.c^2 * ( penalty_parameter_1*e_u' + penalty_parameter_2*F_u');
247 alpha_u = obj.alpha; 259 penalty = obj.Ji*obj.c^2 * (-penalty_parameter_1*e_v' - penalty_parameter_2*F_v');
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 260 end
263 261
264 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. 262 % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
265 % The right boundary is considered the positive boundary 263 % 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) 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)