Mercurial > repos > public > sbplib
changeset 1207:05a01f77d0e3 feature/poroelastic
Swap indices and fix bug in ElasticVariableAnisotropic
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Fri, 20 Sep 2019 14:57:02 -0700 |
| parents | 3258dca12af8 |
| children | 7f427275bc9c |
| files | +scheme/Elastic2dVariableAnisotropic.m |
| diffstat | 1 files changed, 23 insertions(+), 22 deletions(-) [+] |
line wrap: on
line diff
--- a/+scheme/Elastic2dVariableAnisotropic.m Sat Sep 07 13:05:26 2019 -0700 +++ b/+scheme/Elastic2dVariableAnisotropic.m Fri Sep 20 14:57:02 2019 -0700 @@ -140,9 +140,11 @@ end %====== Assemble full operators ======== + I_dim = speye(dim, dim); RHO = spdiag(rho); obj.RHO = RHO; obj.RHOi = inv(RHO); + obj.RHOi_kron = kron(obj.RHOi, I_dim); obj.D1 = cell(dim,1); obj.D2 = cell(dim,dim,dim); @@ -164,7 +166,6 @@ e_scalar_s = e_l{2}; e_scalar_n = e_r{2}; - I_dim = speye(dim, dim); e_w = kron(e_scalar_w, I_dim); e_e = kron(e_scalar_e, I_dim); e_s = kron(e_scalar_s, I_dim); @@ -192,10 +193,10 @@ % D2 - for i = 1:dim - for k = 2:dim - for l = 2:dim - obj.D2{i,k,l} = sparse(m_tot); + for j = 1:dim + for k = 1:dim + for l = 1:dim + obj.D2{j,k,l} = sparse(m_tot,m_tot); end end end @@ -204,11 +205,11 @@ k = 1; for r = 1:m(2) p = ind(:,r); - for i = 1:dim + for j = 1:dim for l = 1:dim - coeff = C{i,k,k,l}; + coeff = C{k,j,k,l}; D_kk = D2{1}(coeff(p)); - obj.D2{i,k,l}(p,p) = D_kk; + obj.D2{j,k,l}(p,p) = D_kk; end end end @@ -216,11 +217,11 @@ k = 2; for r = 1:m(1) p = ind(r,:); - for i = 1:dim + for j = 1:dim for l = 1:dim - coeff = C{i,k,k,l}; + coeff = C{k,j,k,l}; D_kk = D2{2}(coeff(p)); - obj.D2{i,k,l}(p,p) = D_kk; + obj.D2{j,k,l}(p,p) = D_kk; end end end @@ -238,19 +239,20 @@ D2 = obj.D2; D1 = obj.D1; D = sparse(dim*m_tot,dim*m_tot); - d = @kroneckerDelta; % Kronecker delta - db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta for i = 1:dim for j = 1:dim for k = 1:dim for l = 1:dim - D = D + E{i}*inv(RHO)*( d(j,k)*D2{i,k,l}*E{l}' +... - db(j,k)*D1{j}*C_mat{i,j,k,l}*D1{k}*E{l}' ... - ); + if i == k + D = D + E{j}*D2{j,k,l}*E{l}'; + else + D = D + E{j}*D1{i}*C_mat{i,j,k,l}*D1{k}*E{l}'; + end end end end end + D = obj.RHOi_kron*D; obj.D = D; %=========================================%' @@ -288,9 +290,9 @@ % Derivative direction k for k = 1:dim T_l{j}{i,l} = T_l{j}{i,l} ... - - (e_l{j}'*C_mat{i,j,k,l}*D1{k})'; + - (e_l{j}'*C_mat{j,i,k,l}*D1{k})'; T_r{j}{i,l} = T_r{j}{i,l} ... - + (e_r{j}'*C_mat{i,j,k,l}*D1{k})'; + + (e_r{j}'*C_mat{j,i,k,l}*D1{k})'; end tau_l{j}{i} = tau_l{j}{i} + (T_l{j}{i,l}'*E{l}')'; tau_r{j}{i} = tau_r{j}{i} + (T_r{j}{i,l}'*E{l}')'; @@ -344,7 +346,6 @@ obj.tau_n = (e_n'*e1_n*obj.tau1_n')' + (e_n'*e2_n*obj.tau2_n')'; % Kroneckered norms and coefficients - obj.RHOi_kron = kron(obj.RHOi, I_dim); obj.Hi_kron = kron(obj.Hi, I_dim); % Misc. @@ -435,7 +436,7 @@ Z = sparse(m_tot, m_tot); for l = 1:dim for k = 1:dim - Z = Z + nu(l)*C{i,l,k,j}*nu(k); + Z = Z + nu(l)*C{l,i,k,j}*nu(k); end end Z = -tuning*dim/h11*Z; @@ -536,8 +537,8 @@ 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{i,l,k,j}*nu_u(k)*e_u; - Z_v = Z_v + e_v'*nu_v(l)*C_v{i,l,k,j}*nu_v(k)*e_v; + 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; end end Z = -tuning*dim*( 1/(4*h11_u)*Z_u + 1/(4*h11_v)*Z_v );