Mercurial > repos > public > sbplib
changeset 1308:5016f3f3badb feature/poroelastic
Add viscoElastic traction bc for normal-tangential
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Sun, 19 Jul 2020 21:24:48 -0700 |
| parents | fcca6ad8b102 |
| children | f59b849df30f |
| files | +scheme/ViscoElastic2d.m |
| diffstat | 1 files changed, 27 insertions(+), 5 deletions(-) [+] |
line wrap: on
line diff
diff -r fcca6ad8b102 -r 5016f3f3badb +scheme/ViscoElastic2d.m --- a/+scheme/ViscoElastic2d.m Sun Jul 19 20:30:16 2020 -0700 +++ b/+scheme/ViscoElastic2d.m Sun Jul 19 21:24:48 2020 -0700 @@ -40,6 +40,7 @@ e_scalar_w, e_scalar_e, e_scalar_s, e_scalar_n n_w, n_e, n_s, n_n % Physical normals + tangent_w, tangent_e, tangent_s, tangent_n % Physical tangents elasticObj @@ -208,6 +209,11 @@ obj.n_s = elasticObj.n_s; obj.n_n = elasticObj.n_n; + obj.tangent_w = elasticObj.tangent_w; + obj.tangent_e = elasticObj.tangent_e; + obj.tangent_s = elasticObj.tangent_s; + obj.tangent_n = elasticObj.tangent_n; + obj.H_w = elasticObj.H_w; obj.H_e = elasticObj.H_e; obj.H_s = elasticObj.H_s; @@ -267,18 +273,34 @@ switch type case {'F','f','Free','free','traction','Traction','t','T'} + + % Get elastic closure and penalty [closure, penalty] = obj.elasticObj.boundary_condition(boundary, bc, tuning); closure = Eu*closure*Eu'; penalty = Eu*penalty; - j = component; - for i = 1:dim - for k = 1:dim - for l = 1:dim - closure = closure + eU{j}*( (RHO*H)\(C{i,j,k,l}*e*H_gamma*n{i}*e'*eGamma{k,l}') ); + switch component + case 't' + dir = obj.getTangent(boundary); + case 'n' + dir = obj.getNormal(boundary); + case 1 + dir = {1, 0}; + case 2 + dir = {0, 1}; + end + + % Add viscous part of closure + for j = 1:dim + for i = 1:dim + for k = 1:dim + for l = 1:dim + closure = closure + eU{j}*( (RHO*H)\(C{i,j,k,l}*e*dir{j}^2*H_gamma*n{i}*e'*eGamma{k,l}') ); + end end end end + end end