Mercurial > repos > public > sbplib

comparison +scheme/Wave2dCurve.m @ 82:df642750e8f7

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added Dirichelt BC to Wave2dCurve.
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 26 Nov 2015 16:24:25 +0100
parents a8ed986fcf57
children 08fbc284718f
comparison
equal deleted inserted replaced
81:9c0192cf099f 82:df642750e8f7
173 [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv] = obj.get_boundary_ops(boundary); 173 [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv] = obj.get_boundary_ops(boundary);
174 174
175 switch type 175 switch type
176 % Dirichlet boundary condition 176 % Dirichlet boundary condition
177 case {'D','d','dirichlet'} 177 case {'D','d','dirichlet'}
178 error('not implemented') 178 % v denotes the solution in the neighbour domain
179 alpha = obj.alpha; 179 tuning = 1.2;
180 180 % tuning = 20.2;
181 % tau1 < -alpha^2/gamma 181 [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t] = obj.get_boundary_ops(boundary);
182 tuning = 1.1; 182
183 tau1 = -tuning*alpha/gamm; 183 a_n = spdiag(coeff_n);
184 tau2 = s*alpha; 184 a_t = spdiag(coeff_t);
185 185
186 p = tau1*e + tau2*d; 186 F = (s * a_n * d_n' + s * a_t*d_t')';
187 187
188 closure = halfnorm_inv*p*e'; 188 u = obj;
189 189
190 pp = halfnorm_inv*p; 190 b1 = gamm*u.lambda./u.a11.^2;
191 b2 = gamm*u.lambda./u.a22.^2;
192
193 tau = -1./b1 - 1./b2;
194 tau = tuning * spdiag(tau(:));
195 sig1 = 1/2;
196
197 penalty_parameter_1 = halfnorm_inv_n*(tau + sig1*halfnorm_inv_t*F*e'*halfnorm_t)*e;
198
199 closure = obj.Ji*obj.c^2 * penalty_parameter_1*e';
200 pp = -obj.Ji*obj.c^2 * penalty_parameter_1;
191 switch class(data) 201 switch class(data)
192 case 'double' 202 case 'double'
193 penalty = pp*data; 203 penalty = pp*data;
194 case 'function_handle' 204 case 'function_handle'
195 penalty = @(t)pp*data(t); 205 penalty = @(t)pp*data(t);
252 b2_u = gamm_u*u.lambda./u.a22.^2; 262 b2_u = gamm_u*u.lambda./u.a22.^2;
253 b1_v = gamm_v*v.lambda./v.a11.^2; 263 b1_v = gamm_v*v.lambda./v.a11.^2;
254 b2_v = gamm_v*v.lambda./v.a22.^2; 264 b2_v = gamm_v*v.lambda./v.a22.^2;
255 265
256 266
257 tau = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v); 267 tau = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v);
258 m_tot = obj.m(1)*obj.m(2);
259 tau = tuning * spdiag(tau(:)); 268 tau = tuning * spdiag(tau(:));
260 sig1 = 1/2; 269 sig1 = 1/2;
261 sig2 = -1/2*s_u; 270 sig2 = -1/2*s_u;
262 271
263 penalty_parameter_1 = halfnorm_inv_u_n*(tau + sig1*halfnorm_inv_u_t*F_u*e_u'*halfnorm_u_t)*e_u; 272 penalty_parameter_1 = halfnorm_inv_u_n*(tau + sig1*halfnorm_inv_u_t*F_u*e_u'*halfnorm_u_t)*e_u;