sbplib: +scheme/ViscoElastic2d.m comparison

comparison +scheme/ViscoElastic2d.m @ 1313:a4c8bf2371f3 feature/poroelastic

Add viscoelastic fault interface

author	Martin Almquist <malmquist@stanford.edu>
date	Wed, 22 Jul 2020 15:35:25 -0700
parents	6d64b57caf46
children	58df4a35fe43

comparison

equal deleted inserted replaced

-:6d64b57caf46
+:a4c8bf2371f3
 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
+dim = obj.dim;
-e_u       = u.getBoundaryOperatorForScalarField('e', boundary);
-en_u       = u.getBoundaryOperator('en', boundary);
+n       = u.getNormal(boundary);
-tau_n_u     = u.getBoundaryOperator('tau_n', boundary);
+H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
-h11_u     = u.getBorrowing(boundary);
+e       = u.getBoundaryOperatorForScalarField('e', boundary);
-n_u      = u.getNormal(boundary);
+ev      = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
-C_u = u.C;
-Ji_u = u.Ji;
+H       = u.H;
-s_u = spdiag(u.(['s_' boundary]));
+RHO     = u.RHO;
-m_tot_u = u.grid.N();
+ETA     = u.ETA;
+C       = u.C;
-% Operators, v side
+Eu      = u.Eu;
-e_v       = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
+eU      = u.eU;
-en_v       = v.getBoundaryOperator('en', neighbour_boundary);
+eGamma  = u.eGamma;
-tau_n_v     = v.getBoundaryOperator('tau_n', neighbour_boundary);
+Egamma  = u.Egamma;
-h11_v     = v.getBorrowing(neighbour_boundary);
-n_v      = v.getNormal(neighbour_boundary);
+CV       = v.C;
+Ev      = v.Eu;
-C_v = v.C;
+eV      = v.eU;
-Ji_v = v.Ji;
+eGammaV = v.eGamma;
-s_v = spdiag(v.(['s_' neighbour_boundary]));
+nV      = v.getNormal(neighbour_boundary);
-m_tot_v = v.grid.N();
+% Reduce stiffness tensors to boundary size
-% Operators that are only required for own domain
+for i = 1:dim
-Hi      = u.E{1}*inv(u.H)*u.E{1}' + u.E{2}*inv(u.H)*u.E{2}';
+for j = 1:dim
-RHOi    = u.E{1}*inv(u.RHO)*u.E{1}' + u.E{2}*inv(u.RHO)*u.E{2}';
+for k = 1:dim
+for l = 1:dim
-% Shared operators
+C{i,j,k,l} = e'*C{i,j,k,l}*e;
-H_gamma         = u.getBoundaryQuadratureForScalarField(boundary);
+CV{i,j,k,l} = ev'*CV{i,j,k,l}*ev;
-dim             = u.dim;
+end
+end
-% Preallocate
+end
-[~, m_int] = size(H_gamma);
+end
-closure = sparse(dim*m_tot_u, dim*m_tot_u);
-penalty = sparse(dim*m_tot_u, dim*m_tot_v);
+% Get elastic closure and penalty
+[closure, penalty] = obj.elasticObj.interface(boundary, v.elasticObj, neighbour_boundary, type);
-% Continuity of normal displacement, term 1: The symmetric term
+closure = Eu*closure*Eu';
-Z_u = sparse(m_int, m_int);
+penalty = Eu*penalty*Ev';
-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{i,j,k,l}*e_u*n_u{k}*n_u{l};
-Z_v = Z_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};
-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'});
+% We only need the viscous part
+closure_tangential = obj.boundary_condition(boundary, {'t', 't'});
-end
+closure = closure + closure_tangential*Egamma*Egamma';
-% Returns h11 for the boundary specified by the string boundary.
+% ------ Coupling of normal component -----------
-% op -- string
+% Add viscous part of traction coupling
-function h11 = getBorrowing(obj, boundary)
+for i = 1:dim
-assertIsMember(boundary, {'w', 'e', 's', 'n'})
+for j = 1:dim
+for k = 1:dim
-switch boundary
+for l = 1:dim
-case {'w','e'}
+for m = 1:dim
-h11 = obj.refObj.h11{1};
+closure = closure + 1/2*eU{m}*( (RHO*H)\(e*n{m}*H_gamma*n{j}*n{i}*C{i,j,k,l}*e'*eGamma{k,l}') );
-case {'s', 'n'}
+penalty = penalty - 1/2*eU{m}*( (RHO*H)\(e*n{m}*H_gamma*nV{j}*nV{i}*CV{i,j,k,l}*ev'*eGammaV{k,l}') );
-h11 = obj.refObj.h11{2};
+end
 end
+end
+end
+end
+% Add penalty to strain rate eq
+for i = 1:dim
+for j = 1:dim
+for k = 1:dim
+for l = 1:dim
+for m = 1:dim
+closure = closure - 1/2*eGamma{i,j}*( (H*ETA)\(e*n{l}*n{k}*C{i,j,k,l}*H_gamma*n{m}*e'*eU{m}') );
+penalty = penalty - 1/2*eGamma{i,j}*( (H*ETA)\(e*n{l}*n{k}*C{i,j,k,l}*H_gamma*nV{m}*ev'*eV{m}') );
+end
+end
+end
+end
+end
+%-------------------------------------------------
 end
 % Returns the outward unit normal vector for the boundary specified by the string boundary.
 % n is a cell of diagonal matrices for each normal component, n{1} = n_1, n{2} = n_2.
 function n = getNormal(obj, boundary)
 % op -- string
 function o = getBoundaryOperator(obj, op, boundary)
 assertIsMember(boundary, {'w', 'e', 's', 'n'})
 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en', 'et', 'tau_n', 'tau_t'})
-o = obj.([op, '_', boundary]);
+o = obj.elasticObj.([op, '_', boundary]);
 end
 % Returns the boundary operator op for the boundary specified by the string boundary.
 % op -- string

Mercurial > repos > public > sbplib

comparison +scheme/ViscoElastic2d.m @ 1313:a4c8bf2371f3 feature/poroelastic