Mercurial > repos > public > sbplib

comparison +scheme/LaplaceCurvilinear.m @ 1342:d1dad4fbfe22 feature/D2_boundary_opt

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add support for glue interpolation operators, and separate wave speeds across interfaces
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 26 Aug 2022 14:19:29 +0200
parents 84200bbae101
children
comparison
equal deleted inserted replaced
1341:663eb90a4559 1342:d1dad4fbfe22
101 Dv = kr(I_u,D1_v); 101 Dv = kr(I_u,D1_v);
102 obj.Hu = kr(H_u,I_v); 102 obj.Hu = kr(H_u,I_v);
103 obj.Hv = kr(I_u,H_v); 103 obj.Hv = kr(I_u,H_v);
104 obj.Hiu = kr(Hi_u,I_v); 104 obj.Hiu = kr(Hi_u,I_v);
105 obj.Hiv = kr(I_u,Hi_v); 105 obj.Hiv = kr(I_u,Hi_v);
106 obj.H_u = H_u;
107 obj.H_v = H_v;
106 108
107 e_w = kr(e_l_u,I_v); 109 e_w = kr(e_l_u,I_v);
108 e_e = kr(e_r_u,I_v); 110 e_e = kr(e_r_u,I_v);
109 e_s = kr(I_u,e_l_v); 111 e_s = kr(I_u,e_l_v);
110 e_n = kr(I_u,e_r_v); 112 e_n = kr(I_u,e_r_v);
287 default_struct('type', defaultType); 289 default_struct('type', defaultType);
288 290
289 switch type.interpolation 291 switch type.interpolation
290 case {'none', ''} 292 case {'none', ''}
291 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type); 293 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type);
292 case {'op','OP'} 294 case {'op','OP','glue'}
293 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type); 295 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type);
294 otherwise 296 otherwise
295 error('Unknown type of interpolation: %s ', type.interpolation); 297 error('Unknown type of interpolation: %s ', type.interpolation);
296 end 298 end
297 end 299 end
364 % Find the number of grid points along the interface 366 % Find the number of grid points along the interface
365 m_u = size(e_u, 2); 367 m_u = size(e_u, 2);
366 m_v = size(e_v, 2); 368 m_v = size(e_v, 2);
367 369
368 Hi = obj.Hi; 370 Hi = obj.Hi;
369 a = obj.a;
370 371
371 u = obj; 372 u = obj;
372 v = neighbour_scheme; 373 v = neighbour_scheme;
374
375 a_u = u.a;
376 a_v = v.a;
373 377
374 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2; 378 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2;
375 b2_u = gamm_u*u.lambda(I_u)./u.a22(I_u).^2; 379 b2_u = gamm_u*u.lambda(I_u)./u.a22(I_u).^2;
376 b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2; 380 b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2;
377 b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2; 381 b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2;
382 tau_u = tuning * spdiag(tau_u); 386 tau_u = tuning * spdiag(tau_u);
383 tau_v = tuning * spdiag(tau_v); 387 tau_v = tuning * spdiag(tau_v);
384 beta_u = tau_v; 388 beta_u = tau_v;
385 389
386 % Build interpolation operators 390 % Build interpolation operators
387 intOps = interpOpSet(m_u, m_v, obj.order, neighbour_scheme.order); 391 switch type.interpolation
388 Iu2v = intOps.Iu2v; 392 case {'op','OP','MC','mc'}
389 Iv2u = intOps.Iv2u; 393 intOps = interpOpSet(m_u, m_v, obj.order, neighbour_scheme.order);
390 394 Iu2v = intOps.Iu2v;
391 closure = a*Hi*e_u*tau_u*H_b_u*e_u' + ... 395 Iv2u = intOps.Iv2u;
392 a*Hi*e_u*H_b_u*Iv2u.bad*beta_u*Iu2v.good*e_u' + ... 396 case 'glue'
393 a*1/2*Hi*d_u*H_b_u*e_u' + ... 397 optype = type.optype; % TODO: Should be stored by the scheme?
394 -a*1/2*Hi*e_u*H_b_u*d_u'; 398 xu = obj.getBoundaryGrid(boundary);
395 399 xv = neighbour_scheme.getBoundaryGrid(neighbour_boundary);
396 penalty = -a*Hi*e_u*tau_u*H_b_u*Iv2u.good*e_v' + ... 400 hu = obj.getBoundaryGridSpacing(boundary);
397 -a*Hi*e_u*H_b_u*Iv2u.bad*beta_u*e_v' + ... 401 hv = neighbour_scheme.getBoundaryGridSpacing(neighbour_boundary);
398 -a*1/2*Hi*d_u*H_b_u*Iv2u.good*e_v' + ... 402 glueOps = interpOpSet(xu, xv, hu, hv, obj.order, neighbour_scheme.order, optype, optype);
399 -a*1/2*Hi*e_u*H_b_u*Iv2u.bad*d_v'; 403 Iu2v = glueOps.Iu2v;
404 Iv2u = glueOps.Iv2u;
405 otherwise
406 error();
407 end
408
409 closure = a_u*Hi*e_u*tau_u*H_b_u*e_u' + ...
410 a_v*Hi*e_u*H_b_u*Iv2u.bad*beta_u*Iu2v.good*e_u' + ...
411 a_u*1/2*Hi*d_u*H_b_u*e_u' + ...
412 -a_u*1/2*Hi*e_u*H_b_u*d_u';
413
414 penalty = -a_u*Hi*e_u*tau_u*H_b_u*Iv2u.good*e_v' + ...
415 -a_v*Hi*e_u*H_b_u*Iv2u.bad*beta_u*e_v' + ...
416 -a_u*1/2*Hi*d_u*H_b_u*Iv2u.good*e_v' + ...
417 -a_v*1/2*Hi*e_u*H_b_u*Iv2u.bad*d_v';
400 418
401 end 419 end
402 420
403 % Returns the boundary operator op for the boundary specified by the string boundary. 421 % Returns the boundary operator op for the boundary specified by the string boundary.
404 % op -- string 422 % op -- string
449 case {'s','n'} 467 case {'s','n'}
450 gamm = obj.gamm_v; 468 gamm = obj.gamm_v;
451 end 469 end
452 end 470 end
453 471
472 function xb = getBoundaryGrid(obj, boundary)
473 assertIsMember(boundary, {'w', 'e', 's', 'n'})
474 if isa(obj.grid,'grid.Cartesian')
475 lg = obj.grid;
476 elseif isa(obj.grid,'grid.Curvilinear')
477 lg = obj.grid.logicalGrid();
478 else
479 error('Grid type not supported');
480 end
481 switch boundary
482 case {'w','e'}
483 xb = lg.x{2};
484 case {'s','n'}
485 xb = lg.x{1};
486 end
487 end
488
489 function hb = getBoundaryGridSpacing(obj, boundary)
490 h = obj.grid.scaling();
491 switch boundary
492 case {'w','e'}
493 hb = h(2);
494 case {'s','n'}
495 hb = h(1);
496 end
497 end
498
454 function N = size(obj) 499 function N = size(obj)
455 N = prod(obj.m); 500 N = prod(obj.m);
456 end 501 end
457 end 502 end
458 end 503 end