Mercurial > repos > public > sbplib
changeset 1318:6650b111742d feature/poroelastic
Implement normal-tangential interface treatment in viscoElastic2d
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Sun, 26 Jul 2020 20:42:06 -0700 |
| parents | 34997aced843 |
| children | 8b1110385ee2 |
| files | +scheme/ViscoElastic2d.m |
| diffstat | 1 files changed, 95 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- a/+scheme/ViscoElastic2d.m Sun Jul 26 20:06:06 2020 -0700 +++ b/+scheme/ViscoElastic2d.m Sun Jul 26 20:42:06 2020 -0700 @@ -394,6 +394,8 @@ [closure, penalty] = obj.interfaceStandard(boundary,neighbour_scheme,neighbour_boundary,type); case 'frictionalFault' [closure, penalty] = obj.interfaceFrictionalFault(boundary,neighbour_scheme,neighbour_boundary,type); + case 'normalTangential' + [closure, penalty] = obj.interfaceNormalTangential(boundary,neighbour_scheme,neighbour_boundary,type); end end @@ -546,6 +548,99 @@ end + 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; + + dim = obj.dim; + + n = u.getNormal(boundary); + t = u.getTangent(boundary); + H_gamma = u.getBoundaryQuadratureForScalarField(boundary); + e = u.getBoundaryOperatorForScalarField('e', boundary); + + ev = v.getBoundaryOperatorForScalarField('e', neighbour_boundary); + + H = u.H; + RHO = u.RHO; + ETA = u.ETA; + C = u.C; + Eu = u.Eu; + eU = u.eU; + eGamma = u.eGamma; + Egamma = u.Egamma; + + CV = v.C; + Ev = v.Eu; + eV = v.eU; + eGammaV = v.eGamma; + nV = v.getNormal(neighbour_boundary); + tV = v.getTangent(neighbour_boundary); + + % Reduce stiffness tensors to boundary size + for i = 1:dim + for j = 1:dim + for k = 1:dim + for l = 1:dim + C{i,j,k,l} = e'*C{i,j,k,l}*e; + CV{i,j,k,l} = ev'*CV{i,j,k,l}*ev; + end + end + end + end + + % Get elastic closure and penalty + [closure, penalty] = obj.elasticObj.interface(boundary, v.elasticObj, neighbour_boundary, type); + closure = Eu*closure*Eu'; + penalty = Eu*penalty*Ev'; + + % ------ Traction coupling, viscous part ----------- + for i = 1:dim + for j = 1:dim + for k = 1:dim + for l = 1:dim + for m = 1:dim + % Normal component + 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}') ); + 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}') ); + + % Tangential component + closure = closure + 1/2*eU{m}*( (RHO*H)\(e*t{m}*H_gamma*t{j}*n{i}*C{i,j,k,l}*e'*eGamma{k,l}') ); + penalty = penalty - 1/2*eU{m}*( (RHO*H)\(e*t{m}*H_gamma*tV{j}*nV{i}*CV{i,j,k,l}*ev'*eGammaV{k,l}') ); + end + end + end + end + end + %------------------------------------------------- + + % --- Displacement coupling ---------------------- + % 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 + % Normal component + 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}') ); + + % Tangential component + closure = closure - 1/2*eGamma{i,j}*( (H*ETA)\(e*t{l}*n{k}*C{i,j,k,l}*H_gamma*t{m}*e'*eU{m}') ); + penalty = penalty - 1/2*eGamma{i,j}*( (H*ETA)\(e*t{l}*n{k}*C{i,j,k,l}*H_gamma*tV{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)