Mercurial > repos > public > sbplib
changeset 1202:31d7288d0653 feature/poroelastic
Make Anisotropic work for mixed displacement/traction BC at the same boundary.
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Thu, 05 Sep 2019 16:41:49 -0700 |
| parents | 8f4e79aa32ba |
| children | 25cadc69a589 |
| files | +scheme/Elastic2dVariableAnisotropic.m |
| diffstat | 1 files changed, 18 insertions(+), 3 deletions(-) [+] |
line wrap: on
line diff
--- a/+scheme/Elastic2dVariableAnisotropic.m Thu Sep 05 14:13:00 2019 -0700 +++ b/+scheme/Elastic2dVariableAnisotropic.m Thu Sep 05 16:41:49 2019 -0700 @@ -370,7 +370,7 @@ function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning) default_arg('tuning', 1.0); - assert( iscell(bc), 'The BC type must be a 2x1 cell array' ); + assert( iscell(bc), 'The BC type must be a 2x1 or 3x1 cell array' ); comp = bc{1}; type = bc{2}; if ischar(comp) @@ -404,13 +404,28 @@ % OBS! Cannot yet set one component at a time unless one assumes Displacement for all components case {'D','d','dirichlet','Dirichlet','displacement','Displacement'} - % Loop over components that Dirichlet penalties end up on + if numel(bc) >= 3 + dComps = bc{3}; + else + dComps = 1:dim; + end + + % Loops over components that Dirichlet penalties end up on % Y: symmetrizing part of penalty % Z: symmetric part of penalty % X = Y + Z. + + % Nonsymmetric part goes on all components to + % Yield traction in discrete energy rate for i = 1:dim Y = T{j,i}'; + X = e*Y; + closure = closure + E{i}*RHOi*Hi*X'*e*H_gamma*(e'*E{j}' ); + penalty = penalty - E{i}*RHOi*Hi*X'*e*H_gamma; + end + % Symmetric part only required on components with displacement BC + for i = dComps Z = sparse(m_tot, m_tot); for l = 1:dim for k = 1:dim @@ -418,7 +433,7 @@ end end Z = -tuning*dim/h11*Z; - X = Z + e*Y; + X = Z; closure = closure + E{i}*RHOi*Hi*X'*e*H_gamma*(e'*E{j}' ); penalty = penalty - E{i}*RHOi*Hi*X'*e*H_gamma; end