Mercurial > repos > public > sbplib

comparison +scheme/Wave2dCurve.m @ 97:33dba20b4b9d feature/arclen-param

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merged with deafult.
author Martin Almquist <martin.almquist@it.uu.se>
date Wed, 02 Dec 2015 16:29:54 +0100
parents 19d0c9325a3e
children f18142c1530b
comparison
equal deleted inserted replaced
88:7b092f0fd620 97:33dba20b4b9d
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);
232 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 242 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
233 % u denotes the solution in the own domain 243 % u denotes the solution in the own domain
234 % v denotes the solution in the neighbour domain 244 % v denotes the solution in the neighbour domain
235 tuning = 1.2; 245 tuning = 1.2;
236 % tuning = 20.2; 246 % tuning = 20.2;
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] = obj.get_boundary_ops(boundary); 247 [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);
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] = neighbour_scheme.get_boundary_ops(boundary); 248 [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);
239 249
240 a_n_u = spdiag(coeff_n_u); 250 a_n_u = spdiag(coeff_n_u);
241 a_t_u = spdiag(coeff_t_u); 251 a_t_u = spdiag(coeff_t_u);
242 a_n_v = spdiag(coeff_n_v); 252 a_n_v = spdiag(coeff_n_v);
243 a_t_v = spdiag(coeff_t_v); 253 a_t_v = spdiag(coeff_t_v);
246 F_v = (s_v * a_n_v * d_n_v' + s_v * a_t_v*d_t_v')'; 256 F_v = (s_v * a_n_v * d_n_v' + s_v * a_t_v*d_t_v')';
247 257
248 u = obj; 258 u = obj;
249 v = neighbour_scheme; 259 v = neighbour_scheme;
250 260
251 b1_u = gamm_u*u.lambda./u.a11.^2; 261 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2;
252 b2_u = gamm_u*u.lambda./u.a22.^2; 262 b2_u = gamm_u*u.lambda(I_u)./u.a22(I_u).^2;
253 b1_v = gamm_v*v.lambda./v.a11.^2; 263 b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2;
254 b2_v = gamm_v*v.lambda./v.a22.^2; 264 b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2;
255 265
256 266 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);
258 m_tot = obj.m(1)*obj.m(2);
259 tau = tuning * spdiag(tau(:)); 267 tau = tuning * spdiag(tau(:));
260 sig1 = 1/2; 268 sig1 = 1/2;
261 sig2 = -1/2*s_u; 269 sig2 = -1/2;
262 270
263 penalty_parameter_1 = halfnorm_inv_u_n*(tau + sig1*halfnorm_inv_u_t*F_u*e_u'*halfnorm_u_t)*e_u; 271 penalty_parameter_1 = halfnorm_inv_u_n*(e_u*tau + sig1*halfnorm_inv_u_t*F_u*e_u'*halfnorm_u_t*e_u);
264 penalty_parameter_2 = halfnorm_inv_u_n * sig2 * e_u; 272 penalty_parameter_2 = halfnorm_inv_u_n * sig2 * e_u;
265 273
266 274
267 closure = obj.Ji*obj.c^2 * ( penalty_parameter_1*e_u' + penalty_parameter_2*F_u'); 275 closure = obj.Ji*obj.c^2 * ( penalty_parameter_1*e_u' + penalty_parameter_2*F_u');
268 penalty = obj.Ji*obj.c^2 * (-penalty_parameter_1*e_v' - penalty_parameter_2*F_v'); 276 penalty = obj.Ji*obj.c^2 * (-penalty_parameter_1*e_v' + penalty_parameter_2*F_v');
269 end 277 end
270 278
271 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. 279 % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
272 % The right boundary is considered the positive boundary 280 % The right boundary is considered the positive 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) 281 %
282 % I -- the indecies of the boundary points in the grid matrix
283 function [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t, I] = get_boundary_ops(obj,boundary)
284
285 gridMatrix = zeros(obj.m(2),obj.m(1));
286 gridMatrix(:) = 1:numel(gridMatrix);
287
274 switch boundary 288 switch boundary
275 case 'w' 289 case 'w'
276 e = obj.e_w; 290 e = obj.e_w;
277 d_n = obj.du_w; 291 d_n = obj.du_w;
278 d_t = obj.dv_w; 292 d_t = obj.dv_w;
279 s = -1; 293 s = -1;
280 294
281 coeff_n = obj.a11(:,1); 295 I = gridMatrix(:,1);
282 coeff_t = obj.a12(:,1); 296 coeff_n = obj.a11(I);
297 coeff_t = obj.a12(I);
283 case 'e' 298 case 'e'
284 e = obj.e_e; 299 e = obj.e_e;
285 d_n = obj.du_e; 300 d_n = obj.du_e;
286 d_t = obj.dv_e; 301 d_t = obj.dv_e;
287 s = 1; 302 s = 1;
288 303
289 coeff_n = obj.a11(:,end); 304 I = gridMatrix(:,end);
290 coeff_t = obj.a12(:,end); 305 coeff_n = obj.a11(I);
306 coeff_t = obj.a12(I);
291 case 's' 307 case 's'
292 e = obj.e_s; 308 e = obj.e_s;
293 d_n = obj.dv_s; 309 d_n = obj.dv_s;
294 d_t = obj.du_s; 310 d_t = obj.du_s;
295 s = -1; 311 s = -1;
296 312
297 coeff_n = obj.a22(1,:)'; 313 I = gridMatrix(1,:)';
298 coeff_t = obj.a12(1,:)'; 314 coeff_n = obj.a22(I);
315 coeff_t = obj.a12(I);
299 case 'n' 316 case 'n'
300 e = obj.e_n; 317 e = obj.e_n;
301 d_n = obj.dv_n; 318 d_n = obj.dv_n;
302 d_t = obj.du_n; 319 d_t = obj.du_n;
303 s = 1; 320 s = 1;
304 321
305 coeff_n = obj.a22(end,:)'; 322 I = gridMatrix(end,:)';
306 coeff_t = obj.a12(end,:)'; 323 coeff_n = obj.a22(I);
324 coeff_t = obj.a12(I);
307 otherwise 325 otherwise
308 error('No such boundary: boundary = %s',boundary); 326 error('No such boundary: boundary = %s',boundary);
309 end 327 end
310 328
311 switch boundary 329 switch boundary