Mercurial > repos > public > sbplib
changeset 1265:6696b216b982 feature/poroelastic
Start working on displacement bc, clean up traction.
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Thu, 30 Apr 2020 11:30:55 -0700 |
| parents | 066fdfaa2411 |
| children | ad31d9c4cec2 |
| files | +scheme/Elastic2dStaggeredAnisotropic.m |
| diffstat | 1 files changed, 67 insertions(+), 19 deletions(-) [+] |
line wrap: on
line diff
diff -r 066fdfaa2411 -r 6696b216b982 +scheme/Elastic2dStaggeredAnisotropic.m --- a/+scheme/Elastic2dStaggeredAnisotropic.m Wed Apr 29 21:05:01 2020 -0700 +++ b/+scheme/Elastic2dStaggeredAnisotropic.m Thu Apr 30 11:30:55 2020 -0700 @@ -30,6 +30,7 @@ % Traction operators tau_w, tau_e, tau_s, tau_n % Return scalar field + T_w, T_e, T_s, T_n % Act on scalar, return scalar % Inner products H, H_u, H_s @@ -345,12 +346,36 @@ tau_s = cell(nGrids, 1); tau_n = cell(nGrids, 1); + T_w = cell(nGrids, nGrids); + T_e = cell(nGrids, nGrids); + T_s = cell(nGrids, nGrids); + T_n = cell(nGrids, nGrids); + for b = 1:nGrids + [~, m_we] = size(e_w_s{b}); + [~, m_sn] = size(e_s_s{b}); + for c = 1:nGrids + T_w{b,c} = cell(dim, dim); + T_e{b,c} = cell(dim, dim); + T_s{b,c} = cell(dim, dim); + T_n{b,c} = cell(dim, dim); + + 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}); + end + end + end + end + for b = 1:nGrids tau_w{b} = cell(dim, 1); tau_e{b} = cell(dim, 1); tau_s{b} = cell(dim, 1); tau_n{b} = cell(dim, 1); - + for j = 1:dim tau_w{b}{j} = sparse(N, m_s{b}(2)); tau_e{b}{j} = sparse(N, m_s{b}(2)); @@ -369,6 +394,13 @@ tau_e{b}{j} = tau_e{b}{j} + (e_e_s{b}'*n_e(i)*sigma_b_ij)'; tau_s{b}{j} = tau_s{b}{j} + (e_s_s{b}'*n_s(i)*sigma_b_ij)'; tau_n{b}{j} = tau_n{b}{j} + (e_n_s{b}'*n_n(i)*sigma_b_ij)'; + + 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)'; end end end @@ -381,6 +413,11 @@ obj.tau_s = tau_s; obj.tau_n = tau_n; + obj.T_w = T_w; + obj.T_e = T_e; + obj.T_s = T_s; + obj.T_n = T_n; + % D1 = obj.D1; % Traction tensors, T_ij @@ -465,9 +502,10 @@ comp = obj.getComponent(comp, boundary); end - e = obj.getBoundaryOperatorForScalarField('e_u', boundary); + e_u = obj.getBoundaryOperatorForScalarField('e_u', boundary); + e_s = obj.getBoundaryOperatorForScalarField('e_s', boundary); tau = obj.getBoundaryOperator('tau', boundary); - % T = obj.getBoundaryTractionOperator(boundary); + T = obj.getBoundaryTractionOperator(boundary); % h11 = obj.getBorrowing(boundary); H_gamma = obj.getBoundaryQuadratureForScalarField(boundary); nu = obj.getNormal(boundary); @@ -478,13 +516,18 @@ RHO = obj.RHO; C = obj.C; + %---- Grid layout ------- + % gu1 = xp o yp; + % gu2 = xd o yd; + % gs1 = xd o yp; + % gs2 = xp o yd; + %------------------------ + switch boundary - case {'s', 'n'} - e = rot90(e,2); - G = rot90(G,2); - U = rot90(U,2); - H = rot90(H,2); - RHO = rot90(RHO,2); + case {'w', 'e'} + gridCombos = {{1,1}, {2,2}}; + case {'s', 'n'} + gridCombos = {{2,1}, {1,2}}; end dim = obj.dim; @@ -495,7 +538,6 @@ % Preallocate [~, col] = size(tau{1}{1}); closure = sparse(m_tot, m_tot); -% penalty = sparse(m_tot, col); penalty = cell(nGrids, 1); j = comp; @@ -503,7 +545,6 @@ % Dirichlet boundary condition case {'D','d','dirichlet','Dirichlet','displacement','Displacement'} - error('not implemented'); if numel(bc) >= 3 dComps = bc{3}; @@ -518,11 +559,15 @@ % 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; + for a = 1:nGrids + for b = 1:nGrids + 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}; + end + end end % Symmetric part only required on components with displacement BC. @@ -542,9 +587,12 @@ % Free boundary condition case {'F','f','Free','free','traction','Traction','t','T'} - for a = 1:nGrids - closure = closure - G{a}*U{a}{j}*(RHO{a}\(H{a}\(e{a}*H_gamma{a}*tau{a}{j}'))); - penalty{a} = G{a}*U{a}{j}*(RHO{a}\(H{a}\(e{a}*H_gamma{a}))); + for m = 1:numel(gridCombos) + gc = gridCombos{m}; + a = gc{1}; + b = gc{2}; + closure = closure - G{a}*U{a}{j}*(RHO{a}\(H{a}\(e_u{a}*H_gamma{b}*tau{b}{j}'))); + penalty{b} = G{a}*U{a}{j}*(RHO{a}\(H{a}\(e_u{a}*H_gamma{b}))); end % Unknown boundary condition