Mercurial > repos > public > sbplib
changeset 1279:546ee16760d5 feature/poroelastic
Add interfaces in Elastic2dStaggeredAnisotropic
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Sun, 07 Jun 2020 11:30:49 -0700 |
| parents | ce422ba8964e |
| children | c7f5b480e455 |
| files | +scheme/Elastic2dStaggeredAnisotropic.m |
| diffstat | 1 files changed, 103 insertions(+), 76 deletions(-) [+] |
line wrap: on
line diff
--- a/+scheme/Elastic2dStaggeredAnisotropic.m Sun Jun 07 11:28:22 2020 -0700 +++ b/+scheme/Elastic2dStaggeredAnisotropic.m Sun Jun 07 11:30:49 2020 -0700 @@ -103,7 +103,11 @@ for j = 1:dim for k = 1:dim for l = 1:dim - C_mat{a}{i,j,k,l} = spdiag(grid.evalOn(g_s{a}, C{i,j,k,l})); + if numel(C) == nGrids + C_mat{a}{i,j,k,l} = spdiag(C{a}{i,j,k,l}); + else + C_mat{a}{i,j,k,l} = spdiag(grid.evalOn(g_s{a}, C{i,j,k,l})); + end end end end @@ -230,7 +234,7 @@ % D1_u2s{a, b}{i} approximates ddi and - %takes from u grid number b to s grid number a + % takes from u grid number b to s grid number a % Some of D1_x2y{a, b} are 0. D1_u2s = cell(nGrids, nGrids); D1_s2u = cell(nGrids, nGrids); @@ -543,7 +547,7 @@ % Preallocate [~, col] = size(tau{1}{1}); closure = sparse(m_tot, m_tot); - penalty = cell(nGrids, 1); + penalty = cell(1, nGrids); for a = 1:nGrids [~, col] = size(e_u{a}); penalty{a} = sparse(m_tot, col); @@ -619,6 +623,8 @@ otherwise error('No such boundary condition: type = %s',type); end + + penalty = cell2mat(penalty); end % type Struct that specifies the interface coupling. @@ -651,106 +657,127 @@ v = neighbour_scheme; % Operators, u side - e_u = u.getBoundaryOperatorForScalarField('e', boundary); + eu_u = u.getBoundaryOperatorForScalarField('e_u', boundary); + es_u = u.getBoundaryOperatorForScalarField('e_s', boundary); tau_u = u.getBoundaryOperator('tau', boundary); - h11_u = u.getBorrowing(boundary); nu_u = u.getNormal(boundary); - E_u = u.E; + G_u = u.G; + U_u = u.U; C_u = u.C; - m_tot_u = u.grid.N(); + m_tot_u = u.N; % Operators, v side - e_v = v.getBoundaryOperatorForScalarField('e', neighbour_boundary); + eu_v = v.getBoundaryOperatorForScalarField('e_u', neighbour_boundary); + es_v = v.getBoundaryOperatorForScalarField('e_s', neighbour_boundary); tau_v = v.getBoundaryOperator('tau', neighbour_boundary); - h11_v = v.getBorrowing(neighbour_boundary); nu_v = v.getNormal(neighbour_boundary); - E_v = v.E; + G_v = v.G; + U_v = v.U; C_v = v.C; - m_tot_v = v.grid.N(); - - % Fix {'e', 's'}, {'w', 'n'}, and {'x','x'} couplings - flipFlag = false; - e_v_flip = e_v; - if (strcmp(boundary,'s') && strcmp(neighbour_boundary,'e')) || ... - (strcmp(boundary,'e') && strcmp(neighbour_boundary,'s')) || ... - (strcmp(boundary,'w') && strcmp(neighbour_boundary,'n')) || ... - (strcmp(boundary,'n') && strcmp(neighbour_boundary,'w')) || ... - (strcmp(boundary,'s') && strcmp(neighbour_boundary,'s')) || ... - (strcmp(boundary,'n') && strcmp(neighbour_boundary,'n')) || ... - (strcmp(boundary,'w') && strcmp(neighbour_boundary,'w')) || ... - (strcmp(boundary,'e') && strcmp(neighbour_boundary,'e')) - - flipFlag = true; - e_v_flip = fliplr(e_v); - - t1 = tau_v(:,1:2:end-1); - t2 = tau_v(:,2:2:end); - - t1 = fliplr(t1); - t2 = fliplr(t2); - - tau_v(:,1:2:end-1) = t1; - tau_v(:,2:2:end) = t2; - end + m_tot_v = v.N; % Operators that are only required for own domain - Hi = u.Hi_kron; - RHOi = u.RHOi_kron; - e_kron = u.getBoundaryOperator('e', boundary); + % Hi = u.Hi_kron; + % RHOi = u.RHOi_kron; + % e_kron = u.getBoundaryOperator('e', boundary); + H = u.H_u; + RHO = u.RHO; T_u = u.getBoundaryTractionOperator(boundary); % Shared operators H_gamma = u.getBoundaryQuadratureForScalarField(boundary); - H_gamma_kron = u.getBoundaryQuadrature(boundary); + % H_gamma_kron = u.getBoundaryQuadrature(boundary); dim = u.dim; + nGrids = obj.nGrids; % Preallocate - [~, m_int] = size(H_gamma); - closure = sparse(dim*m_tot_u, dim*m_tot_u); - penalty = sparse(dim*m_tot_u, dim*m_tot_v); + % [~, m_int] = size(H_gamma); + closure = sparse(m_tot_u, m_tot_u); + penalty = sparse(m_tot_u, m_tot_v); - % ---- Continuity of displacement ------ - - % Y: symmetrizing part of penalty - % Z: symmetric part of penalty - % X = Y + Z. + %---- Grid layout ------- + % gu1 = xp o yp; + % gu2 = xd o yd; + % gs1 = xd o yp; + % gs2 = xp o yd; + %------------------------ - % Loop over components to couple across interface - for j = 1:dim + switch boundary + case {'w', 'e'} + switch neighbour_boundary + case {'w', 'e'} + gridCombos = {{1,1,1}, {2,2,2}}; + case {'s', 'n'} + gridCombos = {{1,1,2}, {2,2,1}}; + end + case {'s', 'n'} + switch neighbour_boundary + case {'s', 'n'} + gridCombos = {{2,1,1}, {1,2,2}}; + case {'w', 'e'} + gridCombos = {{2,1,2}, {1,2,1}}; + end + end - % Loop over components that penalties end up on - for i = 1:dim - Y = 1/2*T_u{j,i}'; - Z_u = sparse(m_int, m_int); - Z_v = sparse(m_int, m_int); - for l = 1:dim - for k = 1:dim - Z_u = Z_u + e_u'*nu_u(l)*C_u{l,i,k,j}*nu_u(k)*e_u; - Z_v = Z_v + e_v'*nu_v(l)*C_v{l,i,k,j}*nu_v(k)*e_v; + % Symmetrizing part + for c = 1:nGrids + for m = 1:numel(gridCombos) + gc = gridCombos{m}; + a = gc{1}; + b = gc{2}; + + for i = 1:dim + for j = 1:dim + Y = 1/2*T_u{a,c}{j,i}'; + closure = closure + G_u{c}*U_u{c}{i}*((RHO{c}*H{c})\(Y'*H_gamma{a}*(eu_u{b}'*U_u{b}{j}'*G_u{b}') )); + penalty = penalty - G_u{c}*U_u{c}{i}*((RHO{c}*H{c})\(Y'*H_gamma{a}*(eu_v{b}'*U_v{b}{j}'*G_v{b}') )); end end - - if flipFlag - Z_v = rot90(Z_v,2); - end - - Z = -tuning*dim*( 1/(4*h11_u)*Z_u + 1/(4*h11_v)*Z_v ); - X = Y + Z*e_u'; - closure = closure + E_u{i}*X'*H_gamma*e_u'*E_u{j}'; - penalty = penalty - E_u{i}*X'*H_gamma*e_v_flip'*E_v{j}'; - end end - % ---- Continuity of traction ------ - closure = closure - 1/2*e_kron*H_gamma_kron*tau_u'; - penalty = penalty - 1/2*e_kron*H_gamma_kron*tau_v'; + % Symmetric part + for m = 1:numel(gridCombos) + gc = gridCombos{m}; + a = gc{1}; + b = gc{2}; + bv = gc{3}; + + h11_u = u.getBorrowing(b, boundary); + h11_v = v.getBorrowing(bv, neighbour_boundary); - % ---- Multiply by inverse of density x quadraure ---- - closure = RHOi*Hi*closure; - penalty = RHOi*Hi*penalty; + for i = 1:dim + for j = 1:dim + Z_u = 0*es_u{b}'*es_u{b}; + Z_v = 0*es_v{bv}'*es_v{bv}; + for l = 1:dim + for k = 1:dim + Z_u = Z_u + es_u{b}'*nu_u(l)*C_u{b}{l,i,k,j}*nu_u(k)*es_u{b}; + Z_v = Z_v + es_v{bv}'*nu_v(l)*C_v{bv}{l,i,k,j}*nu_v(k)*es_v{bv}; + end + end + Z = -tuning*dim*( 1/(4*h11_u)*Z_u + 1/(4*h11_v)*Z_v ); + % X = es_u{b}'*Z*es_u{b}; + % X = Z; + closure = closure + G_u{a}*U_u{a}{i}*((RHO{a}*H{a})\(eu_u{a}*Z'*H_gamma{b}*(eu_u{a}'*U_u{a}{j}'*G_u{a}' ) )); + penalty = penalty - G_u{a}*U_u{a}{i}*((RHO{a}*H{a})\(eu_u{a}*Z'*H_gamma{b}*(eu_v{a}'*U_v{a}{j}'*G_v{a}' ) )); + end + end + end + + % Continuity of traction + for j = 1:dim + for m = 1:numel(gridCombos) + gc = gridCombos{m}; + a = gc{1}; + b = gc{2}; + bv = gc{3}; + closure = closure - 1/2*G_u{a}*U_u{a}{j}*(RHO{a}\(H{a}\(eu_u{a}*H_gamma{b}*tau_u{b}{j}'))); + penalty = penalty - 1/2*G_u{a}*U_u{a}{j}*(RHO{a}\(H{a}\(eu_u{a}*H_gamma{b}*tau_v{bv}{j}'))); + end + end end @@ -871,7 +898,7 @@ end function N = size(obj) - N = obj.dim*prod(obj.m); + N = length(obj.D); end end end