Mercurial > repos > public > sbplib
changeset 1342:d1dad4fbfe22 feature/D2_boundary_opt
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 | 663eb90a4559 |
| children | 09a5783a3d37 |
| files | +scheme/LaplaceCurvilinear.m |
| diffstat | 1 files changed, 58 insertions(+), 13 deletions(-) [+] |
line wrap: on
line diff
--- a/+scheme/LaplaceCurvilinear.m Fri Jul 22 16:37:49 2022 +0200 +++ b/+scheme/LaplaceCurvilinear.m Fri Aug 26 14:19:29 2022 +0200 @@ -103,6 +103,8 @@ obj.Hv = kr(I_u,H_v); obj.Hiu = kr(Hi_u,I_v); obj.Hiv = kr(I_u,Hi_v); + obj.H_u = H_u; + obj.H_v = H_v; e_w = kr(e_l_u,I_v); e_e = kr(e_r_u,I_v); @@ -289,7 +291,7 @@ switch type.interpolation case {'none', ''} [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type); - case {'op','OP'} + case {'op','OP','glue'} [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type); otherwise error('Unknown type of interpolation: %s ', type.interpolation); @@ -366,11 +368,13 @@ m_v = size(e_v, 2); Hi = obj.Hi; - a = obj.a; u = obj; v = neighbour_scheme; + a_u = u.a; + a_v = v.a; + b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2; b2_u = gamm_u*u.lambda(I_u)./u.a22(I_u).^2; b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2; @@ -384,19 +388,33 @@ beta_u = tau_v; % Build interpolation operators - intOps = interpOpSet(m_u, m_v, obj.order, neighbour_scheme.order); - Iu2v = intOps.Iu2v; - Iv2u = intOps.Iv2u; + switch type.interpolation + case {'op','OP','MC','mc'} + intOps = interpOpSet(m_u, m_v, obj.order, neighbour_scheme.order); + Iu2v = intOps.Iu2v; + Iv2u = intOps.Iv2u; + case 'glue' + optype = type.optype; % TODO: Should be stored by the scheme? + xu = obj.getBoundaryGrid(boundary); + xv = neighbour_scheme.getBoundaryGrid(neighbour_boundary); + hu = obj.getBoundaryGridSpacing(boundary); + hv = neighbour_scheme.getBoundaryGridSpacing(neighbour_boundary); + glueOps = interpOpSet(xu, xv, hu, hv, obj.order, neighbour_scheme.order, optype, optype); + Iu2v = glueOps.Iu2v; + Iv2u = glueOps.Iv2u; + otherwise + error(); + end - closure = a*Hi*e_u*tau_u*H_b_u*e_u' + ... - a*Hi*e_u*H_b_u*Iv2u.bad*beta_u*Iu2v.good*e_u' + ... - a*1/2*Hi*d_u*H_b_u*e_u' + ... - -a*1/2*Hi*e_u*H_b_u*d_u'; + closure = a_u*Hi*e_u*tau_u*H_b_u*e_u' + ... + a_v*Hi*e_u*H_b_u*Iv2u.bad*beta_u*Iu2v.good*e_u' + ... + a_u*1/2*Hi*d_u*H_b_u*e_u' + ... + -a_u*1/2*Hi*e_u*H_b_u*d_u'; - penalty = -a*Hi*e_u*tau_u*H_b_u*Iv2u.good*e_v' + ... - -a*Hi*e_u*H_b_u*Iv2u.bad*beta_u*e_v' + ... - -a*1/2*Hi*d_u*H_b_u*Iv2u.good*e_v' + ... - -a*1/2*Hi*e_u*H_b_u*Iv2u.bad*d_v'; + penalty = -a_u*Hi*e_u*tau_u*H_b_u*Iv2u.good*e_v' + ... + -a_v*Hi*e_u*H_b_u*Iv2u.bad*beta_u*e_v' + ... + -a_u*1/2*Hi*d_u*H_b_u*Iv2u.good*e_v' + ... + -a_v*1/2*Hi*e_u*H_b_u*Iv2u.bad*d_v'; end @@ -451,6 +469,33 @@ end end + function xb = getBoundaryGrid(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + if isa(obj.grid,'grid.Cartesian') + lg = obj.grid; + elseif isa(obj.grid,'grid.Curvilinear') + lg = obj.grid.logicalGrid(); + else + error('Grid type not supported'); + end + switch boundary + case {'w','e'} + xb = lg.x{2}; + case {'s','n'} + xb = lg.x{1}; + end + end + + function hb = getBoundaryGridSpacing(obj, boundary) + h = obj.grid.scaling(); + switch boundary + case {'w','e'} + hb = h(2); + case {'s','n'} + hb = h(1); + end + end + function N = size(obj) N = prod(obj.m); end