Mercurial > repos > public > sbplib
changeset 1266:ad31d9c4cec2 feature/poroelastic
Add displacement BC in StaggeredAnisotropic
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Thu, 30 Apr 2020 20:55:26 -0700 |
| parents | 6696b216b982 |
| children | e61083f178be |
| files | +scheme/Elastic2dStaggeredAnisotropic.m |
| diffstat | 1 files changed, 48 insertions(+), 30 deletions(-) [+] |
line wrap: on
line diff
diff -r 6696b216b982 -r ad31d9c4cec2 +scheme/Elastic2dStaggeredAnisotropic.m --- a/+scheme/Elastic2dStaggeredAnisotropic.m Thu Apr 30 11:30:55 2020 -0700 +++ b/+scheme/Elastic2dStaggeredAnisotropic.m Thu Apr 30 20:55:26 2020 -0700 @@ -149,9 +149,13 @@ end % Borrowing constants - for i = 1:dim - obj.h11{i} = ops{i}.H_dual(1,1); - end + % for i = 1:dim + % obj.h11{i} = ops{i}.H_dual(1,1); + % end + obj.h11{1}{1} = ops{1}.H_dual(1,1); + obj.h11{1}{2} = ops{2}.H_primal(1,1); + obj.h11{2}{1} = ops{1}.H_primal(1,1); + obj.h11{2}{2} = ops{2}.H_dual(1,1); %---- Grid layout ------- % gu1 = xp o yp; @@ -361,10 +365,10 @@ for i = 1:dim for j = 1:dim - T_w{b,c}{i,j} = sparse(m_we, N_u{c}); - T_e{b,c}{i,j} = sparse(m_we, N_u{c}); - T_s{b,c}{i,j} = sparse(m_sn, N_u{c}); - T_n{b,c}{i,j} = sparse(m_sn, N_u{c}); + T_w{b,c}{i,j} = sparse(N_u{c}, m_we); + T_e{b,c}{i,j} = sparse(N_u{c}, m_we); + T_s{b,c}{i,j} = sparse(N_u{c}, m_sn); + T_n{b,c}{i,j} = sparse(N_u{c}, m_sn); end end end @@ -397,10 +401,10 @@ S_bc_ijl = C{b}{i,j,k,l}*D1_u2s{b,c}{k}; - % T_w{b,c}{i,l} = T_w{b,c}{i,l} + (e_w_s{b}'*n_w(i)*S_bc_ijl)'; - % T_e{b,c}{i,l} = T_e{b,c}{i,l} + (e_e_s{b}'*n_e(i)*S_bc_ijl)'; - % T_s{b,c}{i,l} = T_s{b,c}{i,l} + (e_s_s{b}'*n_s(i)*S_bc_ijl)'; - % T_n{b,c}{i,l} = T_n{b,c}{i,l} + (e_n_s{b}'*n_n(i)*S_bc_ijl)'; + T_w{b,c}{j,l} = T_w{b,c}{j,l} + (e_w_s{b}'*n_w(i)*S_bc_ijl)'; + T_e{b,c}{j,l} = T_e{b,c}{j,l} + (e_e_s{b}'*n_e(i)*S_bc_ijl)'; + T_s{b,c}{j,l} = T_s{b,c}{j,l} + (e_s_s{b}'*n_s(i)*S_bc_ijl)'; + T_n{b,c}{j,l} = T_n{b,c}{j,l} + (e_n_s{b}'*n_n(i)*S_bc_ijl)'; end end end @@ -506,7 +510,6 @@ e_s = obj.getBoundaryOperatorForScalarField('e_s', boundary); tau = obj.getBoundaryOperator('tau', boundary); T = obj.getBoundaryTractionOperator(boundary); - % h11 = obj.getBorrowing(boundary); H_gamma = obj.getBoundaryQuadratureForScalarField(boundary); nu = obj.getNormal(boundary); @@ -539,6 +542,10 @@ [~, col] = size(tau{1}{1}); closure = sparse(m_tot, m_tot); penalty = cell(nGrids, 1); + for a = 1:nGrids + [~, col] = size(e_u{a}); + penalty{a} = sparse(m_tot, col); + end j = comp; switch type @@ -559,30 +566,41 @@ % Nonsymmetric part goes on all components to % yield traction in discrete energy rate - for a = 1:nGrids - for b = 1:nGrids + for c = 1:nGrids + for m = 1:numel(gridCombos) + gc = gridCombos{m}; + a = gc{1}; + b = gc{2}; + for i = 1:dim - Y = T{a,b}{j,i}'; - X = e_s{a}*Y; - closure = closure + U{i}*RHOi*Hi*X'*e_s{a}*H_gamma{a}*(e_u{b}'*U{j}'*G{b}' ); - penalty = penalty - U{i}*RHOi*Hi*X'*e_s{a}*H_gamma{a}; + Y = T{a,c}{j,i}'; + closure = closure + G{c}*U{c}{i}*((RHO{c}*H{c})\(Y'*H_gamma{a}*(e_u{b}'*U{b}{j}'*G{b}') )); + penalty{b} = penalty{b} - G{c}*U{c}{i}*((RHO{c}*H{c})\(Y'*H_gamma{a}) ); end end end % Symmetric part only required on components with displacement BC. % (Otherwise it's not symmetric.) - for i = dComps - Z = sparse(m_tot, m_tot); - for l = 1:dim - for k = 1:dim - Z = Z + nu(l)*C{l,i,k,j}*nu(k); + for m = 1:numel(gridCombos) + gc = gridCombos{m}; + a = gc{1}; + b = gc{2}; + + h11 = obj.getBorrowing(b, boundary); + + for i = dComps + Z = 0*C{b}{1,1,1,1}; + for l = 1:dim + for k = 1:dim + Z = Z + nu(l)*C{b}{l,i,k,j}*nu(k); + end end + Z = -tuning*dim/h11*Z; + X = e_s{b}'*Z*e_s{b}; + closure = closure + G{a}*U{a}{i}*((RHO{a}*H{a})\(e_u{a}*X'*H_gamma{b}*(e_u{a}'*U{a}{j}'*G{a}' ) )); + penalty{a} = penalty{a} - G{a}*U{a}{i}*((RHO{a}*H{a})\(e_u{a}*X'*H_gamma{b} )); end - Z = -tuning*dim/h11*Z; - 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 % Free boundary condition @@ -764,14 +782,14 @@ % Returns h11 for the boundary specified by the string boundary. % op -- string - function h11 = getBorrowing(obj, boundary) + function h11 = getBorrowing(obj, stressGrid, boundary) assertIsMember(boundary, {'w', 'e', 's', 'n'}) switch boundary case {'w','e'} - h11 = obj.h11{1}; + h11 = obj.h11{stressGrid}{1}; case {'s', 'n'} - h11 = obj.h11{2}; + h11 = obj.h11{stressGrid}{2}; end end