Mercurial > repos > public > sbplib
changeset 1294:fe712a1fca3f feature/poroelastic
Add frictional fault interface in Elastic2dAnisotropicCurve
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Thu, 02 Jul 2020 15:00:21 -0700 |
| parents | 3e2c1df740df |
| children | cb053fabbedc |
| files | +scheme/Elastic2dCurvilinearAnisotropic.m |
| diffstat | 1 files changed, 171 insertions(+), 2 deletions(-) [+] |
line wrap: on
line diff
--- a/+scheme/Elastic2dCurvilinearAnisotropic.m Mon Jun 29 15:26:36 2020 -0700 +++ b/+scheme/Elastic2dCurvilinearAnisotropic.m Thu Jul 02 15:00:21 2020 -0700 @@ -35,6 +35,7 @@ tau_w, tau_e, tau_s, tau_n % Return vector field tau1_w, tau1_e, tau1_s, tau1_n % Return scalar field tau2_w, tau2_e, tau2_s, tau2_n % Return scalar field + tau_n_w, tau_n_e, tau_n_s, tau_n_n % Return scalar field % Inner products H @@ -319,6 +320,16 @@ obj.et_s = (obj.tangent_s{1}*obj.e1_s' + obj.tangent_s{2}*obj.e2_s')'; obj.et_n = (obj.tangent_n{1}*obj.e1_n' + obj.tangent_n{2}*obj.e2_n')'; + % obj.tau_n_w = (obj.en_w'*obj.e_w*obj.tau_w')'; + % obj.tau_n_e = (obj.en_e'*obj.e_e*obj.tau_e')'; + % obj.tau_n_s = (obj.en_s'*obj.e_s*obj.tau_s')'; + % obj.tau_n_n = (obj.en_n'*obj.e_n*obj.tau_n')'; + + obj.tau_n_w = (obj.n_w{1}*obj.tau1_w' + obj.n_w{2}*obj.tau2_w')'; + obj.tau_n_e = (obj.n_e{1}*obj.tau1_e' + obj.n_e{2}*obj.tau2_e')'; + obj.tau_n_s = (obj.n_s{1}*obj.tau1_s' + obj.n_s{2}*obj.tau2_s')'; + obj.tau_n_n = (obj.n_n{1}*obj.tau1_n' + obj.n_n{2}*obj.tau2_n')'; + % Stress operators sigma = cell(dim, dim); D1 = {obj.Dx, obj.Dy}; @@ -412,6 +423,75 @@ end end + function [closure, penalty] = displacementBCNormalTangential(boundary, bc, tuning) + error('not implemented'); + u = obj; + + component = bc{1}; + type = bc{2}; + + switch component + case 'n' + en_u = u.getBoundaryOperator('en', boundary); + tau_n_u = u.getBoundaryOperator('tau_n', boundary); + case 't' + en_u = u.getBoundaryOperator('et', boundary); + tau_n_u = u.getBoundaryOperator('tau_t', boundary); + end + + % Operators + e_u = u.getBoundaryOperatorForScalarField('e', boundary); + h11_u = u.getBorrowing(boundary); + n_u = u.getNormal(boundary); + + C_u = u.C; + Ji_u = u.Ji; + s_u = spdiag(u.(['s_' boundary])); + m_tot_u = u.grid.N(); + + 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 normal displacement, term 1: The symmetric term + Z_u = sparse(m_int, m_int); + Z_v = sparse(m_int, m_int); + for i = 1:dim + for j = 1:dim + for l = 1:dim + for k = 1:dim + Z_u = Z_u + n_u{i}*n_u{j}*e_u'*Ji_u*C_u{l,i,k,j}*e_u*n_u{k}*n_u{l}; + Z_v = Z_v + n_v{i}*n_v{j}*e_v'*Ji_v*C_v{l,i,k,j}*e_v*n_v{k}*n_v{l}; + end + end + end + end + + Z = -tuning*dim*( 1/(4*h11_u)*s_u*Z_u + 1/(4*h11_v)*s_v*Z_v ); + closure = closure + en_u*H_gamma*Z*en_u'; + penalty = penalty + en_u*H_gamma*Z*en_v'; + + % Continuity of normal displacement, term 2: The symmetrizing term + closure = closure + 1/2*tau_n_u*H_gamma*en_u'; + penalty = penalty + 1/2*tau_n_u*H_gamma*en_v'; + + % 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'; + + % Multiply all normal component terms by inverse of density x quadraure + closure = RHOi*Hi*closure; + penalty = RHOi*Hi*penalty; + end + % type Struct that specifies the interface coupling. % Fields: % -- tuning: penalty strength, defaults to 1.0 @@ -420,11 +500,100 @@ defaultType.tuning = 1.0; defaultType.interpolation = 'none'; + defaultType.type = 'standard'; default_struct('type', defaultType); - [closure, penalty] = obj.refObj.interface(boundary,neighbour_scheme.refObj,neighbour_boundary,type); + switch type.type + case 'standard' + [closure, penalty] = obj.refObj.interface(boundary,neighbour_scheme.refObj,neighbour_boundary,type); + case 'frictionalFault' + [closure, penalty] = obj.interfaceFrictionalFault(boundary,neighbour_scheme,neighbour_boundary,type); + end + end + function [closure, penalty] = interfaceFrictionalFault(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); + tau_n_u = u.getBoundaryOperator('tau_n', boundary); + h11_u = u.getBorrowing(boundary); + n_u = u.getNormal(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); + tau_n_v = v.getBoundaryOperator('tau_n', neighbour_boundary); + h11_v = v.getBorrowing(neighbour_boundary); + n_v = v.getNormal(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 normal displacement, term 1: The symmetric term + Z_u = sparse(m_int, m_int); + Z_v = sparse(m_int, m_int); + for i = 1:dim + for j = 1:dim + for l = 1:dim + for k = 1:dim + Z_u = Z_u + n_u{i}*n_u{j}*e_u'*Ji_u*C_u{l,i,k,j}*e_u*n_u{k}*n_u{l}; + Z_v = Z_v + n_v{i}*n_v{j}*e_v'*Ji_v*C_v{l,i,k,j}*e_v*n_v{k}*n_v{l}; + end + end + end + end + + Z = -tuning*dim*( 1/(4*h11_u)*s_u*Z_u + 1/(4*h11_v)*s_v*Z_v ); + closure = closure + en_u*H_gamma*Z*en_u'; + penalty = penalty + en_u*H_gamma*Z*en_v'; + + % Continuity of normal displacement, term 2: The symmetrizing term + closure = closure + 1/2*tau_n_u*H_gamma*en_u'; + penalty = penalty + 1/2*tau_n_u*H_gamma*en_v'; + + % 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'; + + % Multiply all normal component terms by inverse of density x quadraure + closure = RHOi*Hi*closure; + penalty = RHOi*Hi*penalty; + + % ---- Tangential tractions are imposed just like traction BC ------ + closure = closure + obj.boundary_condition(boundary, {'t', 't'}); + + end + + % Returns h11 for the boundary specified by the string boundary. % op -- string function h11 = getBorrowing(obj, boundary) @@ -450,7 +619,7 @@ % op -- string function o = getBoundaryOperator(obj, op, boundary) assertIsMember(boundary, {'w', 'e', 's', 'n'}) - assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en', 'et'}) + assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en', 'et', 'tau_n'}) o = obj.([op, '_', boundary]);