Mercurial > repos > public > sbplib
changeset 1291:a8e730db76e9 feature/poroelastic
Add normal-tangential traction BC ElasticCurvilinearAnisotropic
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Mon, 29 Jun 2020 15:23:01 -0700 |
| parents | 01a0500de446 |
| children | 3c87e8d08076 |
| files | +scheme/Elastic2dCurvilinearAnisotropic.m |
| diffstat | 1 files changed, 49 insertions(+), 3 deletions(-) [+] |
line wrap: on
line diff
diff -r 01a0500de446 -r a8e730db76e9 +scheme/Elastic2dCurvilinearAnisotropic.m --- a/+scheme/Elastic2dCurvilinearAnisotropic.m Sun Jun 07 11:33:44 2020 -0700 +++ b/+scheme/Elastic2dCurvilinearAnisotropic.m Mon Jun 29 15:23:01 2020 -0700 @@ -26,6 +26,7 @@ Dx, Dy % Physical derivatives sigma % Cell matrix of physical stress operators n_w, n_e, n_s, n_n % Physical normals + tangent_w, tangent_e, tangent_s, tangent_n % Physical tangents % Boundary operators in cell format, used for BC T_w, T_e, T_s, T_n @@ -50,6 +51,7 @@ e2_w, e2_e, e2_s, e2_n % Act on vector field, return scalar field at boundary e_scalar_w, e_scalar_e, e_scalar_s, e_scalar_n; % Act on scalar field, return scalar field en_w, en_e, en_s, en_n % Act on vector field, return normal component + et_w, et_e, et_s, et_n % Act on vector field, return tangential component % E{i}^T picks out component i E @@ -280,25 +282,30 @@ obj.n_s = cell(2,1); obj.n_n = cell(2,1); + % Compute normal and rotate (exactly!) 90 degrees counter-clockwise to get tangent n_w_1 = (1./s_w).*e_w'*(J.*(K{1,1}*nu_w(1) + K{1,2}*nu_w(2))); n_w_2 = (1./s_w).*e_w'*(J.*(K{2,1}*nu_w(1) + K{2,2}*nu_w(2))); obj.n_w{1} = spdiag(n_w_1); obj.n_w{2} = spdiag(n_w_2); + obj.tangent_w = {-obj.n_w{2}, obj.n_w{1}}; n_e_1 = (1./s_e).*e_e'*(J.*(K{1,1}*nu_e(1) + K{1,2}*nu_e(2))); n_e_2 = (1./s_e).*e_e'*(J.*(K{2,1}*nu_e(1) + K{2,2}*nu_e(2))); obj.n_e{1} = spdiag(n_e_1); obj.n_e{2} = spdiag(n_e_2); + obj.tangent_e = {-obj.n_e{2}, obj.n_e{1}}; n_s_1 = (1./s_s).*e_s'*(J.*(K{1,1}*nu_s(1) + K{1,2}*nu_s(2))); n_s_2 = (1./s_s).*e_s'*(J.*(K{2,1}*nu_s(1) + K{2,2}*nu_s(2))); obj.n_s{1} = spdiag(n_s_1); obj.n_s{2} = spdiag(n_s_2); + obj.tangent_s = {-obj.n_s{2}, obj.n_s{1}}; n_n_1 = (1./s_n).*e_n'*(J.*(K{1,1}*nu_n(1) + K{1,2}*nu_n(2))); n_n_2 = (1./s_n).*e_n'*(J.*(K{2,1}*nu_n(1) + K{2,2}*nu_n(2))); obj.n_n{1} = spdiag(n_n_1); obj.n_n{2} = spdiag(n_n_2); + obj.tangent_n = {-obj.n_n{2}, obj.n_n{1}}; % Operators that extract the normal component obj.en_w = (obj.n_w{1}*obj.e1_w' + obj.n_w{2}*obj.e2_w')'; @@ -306,6 +313,12 @@ obj.en_s = (obj.n_s{1}*obj.e1_s' + obj.n_s{2}*obj.e2_s')'; obj.en_n = (obj.n_n{1}*obj.e1_n' + obj.n_n{2}*obj.e2_n')'; + % Operators that extract the tangential component + obj.et_w = (obj.tangent_w{1}*obj.e1_w' + obj.tangent_w{2}*obj.e2_w')'; + obj.et_e = (obj.tangent_e{1}*obj.e1_e' + obj.tangent_e{2}*obj.e2_e')'; + 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')'; + % Stress operators sigma = cell(dim, dim); D1 = {obj.Dx, obj.Dy}; @@ -354,15 +367,48 @@ default_arg('tuning', 1.0); assert( iscell(bc), 'The BC type must be a 2x1 or 3x1 cell array' ); - [closure, penalty] = obj.refObj.boundary_condition(boundary, bc, tuning); + component = bc{1}; + type = bc{2}; + + switch component + + % If conditions on Cartesian components + case {1,2} + [closure, penalty] = obj.refObj.boundary_condition(boundary, bc, tuning); + + % If conditions on normal or tangential components + case {'n', 't'} - type = bc{2}; + switch component + case 'n' + en = obj.getBoundaryOperator('en', boundary); + case 't' + en = obj.getBoundaryOperator('et', boundary); + end + e1 = obj.getBoundaryOperator('e1', boundary); + + bc = {1, type}; + [c1, p1] = obj.refObj.boundary_condition(boundary, bc, tuning); + bc = {2, type}; + c2 = obj.refObj.boundary_condition(boundary, bc, tuning); + + switch type + case {'F','f','Free','free','traction','Traction','t','T'} + closure = en*en'*(c1+c2); + penalty = en*e1'*p1; + case {'D','d','dirichlet','Dirichlet','displacement','Displacement'} + error('Not implemented'); + end + + end switch type case {'F','f','Free','free','traction','Traction','t','T'} + s = obj.(['s_' boundary]); s = spdiag(s); penalty = penalty*s; + end end @@ -404,7 +450,7 @@ % op -- string function o = getBoundaryOperator(obj, op, boundary) assertIsMember(boundary, {'w', 'e', 's', 'n'}) - assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en'}) + assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en', 'et'}) o = obj.([op, '_', boundary]);