Mercurial > repos > public > sbplib
changeset 1315:6c308da9dcbc feature/poroelastic
Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Sun, 26 Jul 2020 17:38:28 -0700 |
| parents | 58df4a35fe43 |
| children | bd8e607768ce |
| files | +scheme/Elastic2dCurvilinearAnisotropic.m |
| diffstat | 1 files changed, 119 insertions(+), 1 deletions(-) [+] |
line wrap: on
line diff
diff -r 58df4a35fe43 -r 6c308da9dcbc +scheme/Elastic2dCurvilinearAnisotropic.m --- a/+scheme/Elastic2dCurvilinearAnisotropic.m Sat Jul 25 15:56:31 2020 -0700 +++ b/+scheme/Elastic2dCurvilinearAnisotropic.m Sun Jul 26 17:38:28 2020 -0700 @@ -486,7 +486,7 @@ closure = closure + tau_n*H_gamma*en'; penalty = penalty - tau_n*H_gamma; - % Multiply all normal component terms by inverse of density x quadraure + % Multiply all terms by inverse of density x quadraure closure = RHOi*Hi*closure; penalty = RHOi*Hi*penalty; end @@ -505,6 +505,8 @@ switch type.type case 'standard' [closure, penalty] = obj.refObj.interface(boundary,neighbour_scheme.refObj,neighbour_boundary,type); + case 'normalTangential' + [closure, penalty] = obj.interfaceNormalTangential(boundary,neighbour_scheme,neighbour_boundary,type); case 'frictionalFault' [closure, penalty] = obj.interfaceFrictionalFault(boundary,neighbour_scheme,neighbour_boundary,type); end @@ -592,6 +594,122 @@ end + % Same interface conditions as in interfaceStandard, but imposed in the normal-tangential + % coordinate system. For the isotropic case, the components decouple nicely. + % The resulting scheme is not identical to that of interfaceStandard. This appears to be better. + function [closure, penalty] = interfaceNormalTangential(obj,boundary,neighbour_scheme,neighbour_boundary,type) + tuning = type.tuning; + + % u denotes the solution in the own domain + % v denotes the solution in the neighbour domain + + u = obj; + v = neighbour_scheme; + + % Operators, u side + e_u = u.getBoundaryOperatorForScalarField('e', boundary); + en_u = u.getBoundaryOperator('en', boundary); + et_u = u.getBoundaryOperator('et', boundary); + tau_n_u = u.getBoundaryOperator('tau_n', boundary); + tau_t_u = u.getBoundaryOperator('tau_t', boundary); + h11_u = u.getBorrowing(boundary); + n_u = u.getNormal(boundary); + t_u = u.getTangent(boundary); + + C_u = u.C; + Ji_u = u.Ji; + s_u = spdiag(u.(['s_' boundary])); + m_tot_u = u.grid.N(); + + % Operators, v side + e_v = v.getBoundaryOperatorForScalarField('e', neighbour_boundary); + en_v = v.getBoundaryOperator('en', neighbour_boundary); + et_v = v.getBoundaryOperator('et', neighbour_boundary); + tau_n_v = v.getBoundaryOperator('tau_n', neighbour_boundary); + tau_t_v = v.getBoundaryOperator('tau_t', neighbour_boundary); + h11_v = v.getBorrowing(neighbour_boundary); + n_v = v.getNormal(neighbour_boundary); + t_v = v.getTangent(neighbour_boundary); + + C_v = v.C; + Ji_v = v.Ji; + s_v = spdiag(v.(['s_' neighbour_boundary])); + m_tot_v = v.grid.N(); + + % Operators that are only required for own domain + Hi = u.E{1}*inv(u.H)*u.E{1}' + u.E{2}*inv(u.H)*u.E{2}'; + RHOi = u.E{1}*inv(u.RHO)*u.E{1}' + u.E{2}*inv(u.RHO)*u.E{2}'; + + % Shared operators + H_gamma = u.getBoundaryQuadratureForScalarField(boundary); + dim = u.dim; + + % Preallocate + [~, m_int] = size(H_gamma); + closure = sparse(dim*m_tot_u, dim*m_tot_u); + penalty = sparse(dim*m_tot_u, dim*m_tot_v); + + % -- Continuity of displacement, term 1: The symmetric term --- + Zn_u = sparse(m_int, m_int); + Zn_v = sparse(m_int, m_int); + Zt_u = sparse(m_int, m_int); + Zt_v = sparse(m_int, m_int); + for i = 1:dim + for j = 1:dim + for l = 1:dim + for k = 1:dim + % Penalty strength for normal component + Zn_u = Zn_u + n_u{i}*n_u{j}*e_u'*Ji_u*C_u{i,j,k,l}*e_u*n_u{k}*n_u{l}; + Zn_v = Zn_v + n_v{i}*n_v{j}*e_v'*Ji_v*C_v{i,j,k,l}*e_v*n_v{k}*n_v{l}; + + % Penalty strength for tangential component + Zt_u = Zt_u + n_u{i}*t_u{j}*e_u'*Ji_u*C_u{i,j,k,l}*e_u*n_u{k}*t_u{l}; + Zt_v = Zt_v + n_v{i}*t_v{j}*e_v'*Ji_v*C_v{i,j,k,l}*e_v*n_v{k}*t_v{l}; + end + end + end + end + + Zn = -tuning*dim*( 1/(4*h11_u)*s_u*Zn_u + 1/(4*h11_v)*s_v*Zn_v ); + Zt = -tuning*dim*( 1/(4*h11_u)*s_u*Zt_u + 1/(4*h11_v)*s_v*Zt_v ); + + % Continuity of normal component + closure = closure + en_u*H_gamma*Zn*en_u'; + penalty = penalty + en_u*H_gamma*Zn*en_v'; + + % Continuity of tangential component + closure = closure + et_u*H_gamma*Zt*et_u'; + penalty = penalty + et_u*H_gamma*Zt*et_v'; + %------------------------------------------------------------------ + + % --- Continuity of displacement, term 2: The symmetrizing term + + % Continuity of normal displacement + closure = closure + 1/2*tau_n_u*H_gamma*en_u'; + penalty = penalty + 1/2*tau_n_u*H_gamma*en_v'; + + % Continuity of tangential displacement + closure = closure + 1/2*tau_t_u*H_gamma*et_u'; + penalty = penalty + 1/2*tau_t_u*H_gamma*et_v'; + % ------------------------------------------------------------------ + + % --- Continuity of tractions ----------------------------- + + % Continuity of normal traction + closure = closure - 1/2*en_u*H_gamma*tau_n_u'; + penalty = penalty + 1/2*en_u*H_gamma*tau_n_v'; + + % Continuity of tangential traction + closure = closure - 1/2*et_u*H_gamma*tau_t_u'; + penalty = penalty + 1/2*et_u*H_gamma*tau_t_v'; + %-------------------------------------------------------------------- + + % Multiply all terms by inverse of density x quadraure + closure = RHOi*Hi*closure; + penalty = RHOi*Hi*penalty; + + end + % Returns h11 for the boundary specified by the string boundary. % op -- string