Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariableAnisotropic.m @ 1207:05a01f77d0e3 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Swap indices and fix bug in ElasticVariableAnisotropic
author Martin Almquist <malmquist@stanford.edu>
date Fri, 20 Sep 2019 14:57:02 -0700
parents 687515778437
children 7f427275bc9c
comparison
equal deleted inserted replaced
1205:3258dca12af8 1207:05a01f77d0e3
138 d1_0{i} = ops{i}.d1_l; 138 d1_0{i} = ops{i}.d1_l;
139 d1_m{i} = ops{i}.d1_r; 139 d1_m{i} = ops{i}.d1_r;
140 end 140 end
141 141
142 %====== Assemble full operators ======== 142 %====== Assemble full operators ========
143 I_dim = speye(dim, dim);
143 RHO = spdiag(rho); 144 RHO = spdiag(rho);
144 obj.RHO = RHO; 145 obj.RHO = RHO;
145 obj.RHOi = inv(RHO); 146 obj.RHOi = inv(RHO);
147 obj.RHOi_kron = kron(obj.RHOi, I_dim);
146 148
147 obj.D1 = cell(dim,1); 149 obj.D1 = cell(dim,1);
148 obj.D2 = cell(dim,dim,dim); 150 obj.D2 = cell(dim,dim,dim);
149 151
150 % D1 152 % D1
162 e_scalar_w = e_l{1}; 164 e_scalar_w = e_l{1};
163 e_scalar_e = e_r{1}; 165 e_scalar_e = e_r{1};
164 e_scalar_s = e_l{2}; 166 e_scalar_s = e_l{2};
165 e_scalar_n = e_r{2}; 167 e_scalar_n = e_r{2};
166 168
167 I_dim = speye(dim, dim);
168 e_w = kron(e_scalar_w, I_dim); 169 e_w = kron(e_scalar_w, I_dim);
169 e_e = kron(e_scalar_e, I_dim); 170 e_e = kron(e_scalar_e, I_dim);
170 e_s = kron(e_scalar_s, I_dim); 171 e_s = kron(e_scalar_s, I_dim);
171 e_n = kron(e_scalar_n, I_dim); 172 e_n = kron(e_scalar_n, I_dim);
172 173
190 e2_s = (e_scalar_s'*E{2}')'; 191 e2_s = (e_scalar_s'*E{2}')';
191 e2_n = (e_scalar_n'*E{2}')'; 192 e2_n = (e_scalar_n'*E{2}')';
192 193
193 194
194 % D2 195 % D2
195 for i = 1:dim 196 for j = 1:dim
196 for k = 2:dim 197 for k = 1:dim
197 for l = 2:dim 198 for l = 1:dim
198 obj.D2{i,k,l} = sparse(m_tot); 199 obj.D2{j,k,l} = sparse(m_tot,m_tot);
199 end 200 end
200 end 201 end
201 end 202 end
202 ind = grid.funcToMatrix(g, 1:m_tot); 203 ind = grid.funcToMatrix(g, 1:m_tot);
203 204
204 k = 1; 205 k = 1;
205 for r = 1:m(2) 206 for r = 1:m(2)
206 p = ind(:,r); 207 p = ind(:,r);
207 for i = 1:dim 208 for j = 1:dim
208 for l = 1:dim 209 for l = 1:dim
209 coeff = C{i,k,k,l}; 210 coeff = C{k,j,k,l};
210 D_kk = D2{1}(coeff(p)); 211 D_kk = D2{1}(coeff(p));
211 obj.D2{i,k,l}(p,p) = D_kk; 212 obj.D2{j,k,l}(p,p) = D_kk;
212 end 213 end
213 end 214 end
214 end 215 end
215 216
216 k = 2; 217 k = 2;
217 for r = 1:m(1) 218 for r = 1:m(1)
218 p = ind(r,:); 219 p = ind(r,:);
219 for i = 1:dim 220 for j = 1:dim
220 for l = 1:dim 221 for l = 1:dim
221 coeff = C{i,k,k,l}; 222 coeff = C{k,j,k,l};
222 D_kk = D2{2}(coeff(p)); 223 D_kk = D2{2}(coeff(p));
223 obj.D2{i,k,l}(p,p) = D_kk; 224 obj.D2{j,k,l}(p,p) = D_kk;
224 end 225 end
225 end 226 end
226 end 227 end
227 228
228 % Quadratures 229 % Quadratures
236 237
237 % Differentiation matrix D (without SAT) 238 % Differentiation matrix D (without SAT)
238 D2 = obj.D2; 239 D2 = obj.D2;
239 D1 = obj.D1; 240 D1 = obj.D1;
240 D = sparse(dim*m_tot,dim*m_tot); 241 D = sparse(dim*m_tot,dim*m_tot);
241 d = @kroneckerDelta; % Kronecker delta
242 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
243 for i = 1:dim 242 for i = 1:dim
244 for j = 1:dim 243 for j = 1:dim
245 for k = 1:dim 244 for k = 1:dim
246 for l = 1:dim 245 for l = 1:dim
247 D = D + E{i}*inv(RHO)*( d(j,k)*D2{i,k,l}*E{l}' +... 246 if i == k
248 db(j,k)*D1{j}*C_mat{i,j,k,l}*D1{k}*E{l}' ... 247 D = D + E{j}*D2{j,k,l}*E{l}';
249 ); 248 else
249 D = D + E{j}*D1{i}*C_mat{i,j,k,l}*D1{k}*E{l}';
250 end
250 end 251 end
251 end 252 end
252 end 253 end
253 end 254 end
255 D = obj.RHOi_kron*D;
254 obj.D = D; 256 obj.D = D;
255 %=========================================%' 257 %=========================================%'
256 258
257 % Numerical traction operators for BC. 259 % Numerical traction operators for BC.
258 % 260 %
286 T_r{j}{i,l} = sparse(m_tot, n_r); 288 T_r{j}{i,l} = sparse(m_tot, n_r);
287 289
288 % Derivative direction k 290 % Derivative direction k
289 for k = 1:dim 291 for k = 1:dim
290 T_l{j}{i,l} = T_l{j}{i,l} ... 292 T_l{j}{i,l} = T_l{j}{i,l} ...
291 - (e_l{j}'*C_mat{i,j,k,l}*D1{k})'; 293 - (e_l{j}'*C_mat{j,i,k,l}*D1{k})';
292 T_r{j}{i,l} = T_r{j}{i,l} ... 294 T_r{j}{i,l} = T_r{j}{i,l} ...
293 + (e_r{j}'*C_mat{i,j,k,l}*D1{k})'; 295 + (e_r{j}'*C_mat{j,i,k,l}*D1{k})';
294 end 296 end
295 tau_l{j}{i} = tau_l{j}{i} + (T_l{j}{i,l}'*E{l}')'; 297 tau_l{j}{i} = tau_l{j}{i} + (T_l{j}{i,l}'*E{l}')';
296 tau_r{j}{i} = tau_r{j}{i} + (T_r{j}{i,l}'*E{l}')'; 298 tau_r{j}{i} = tau_r{j}{i} + (T_r{j}{i,l}'*E{l}')';
297 end 299 end
298 end 300 end
342 obj.tau_e = (e_e'*e1_e*obj.tau1_e')' + (e_e'*e2_e*obj.tau2_e')'; 344 obj.tau_e = (e_e'*e1_e*obj.tau1_e')' + (e_e'*e2_e*obj.tau2_e')';
343 obj.tau_s = (e_s'*e1_s*obj.tau1_s')' + (e_s'*e2_s*obj.tau2_s')'; 345 obj.tau_s = (e_s'*e1_s*obj.tau1_s')' + (e_s'*e2_s*obj.tau2_s')';
344 obj.tau_n = (e_n'*e1_n*obj.tau1_n')' + (e_n'*e2_n*obj.tau2_n')'; 346 obj.tau_n = (e_n'*e1_n*obj.tau1_n')' + (e_n'*e2_n*obj.tau2_n')';
345 347
346 % Kroneckered norms and coefficients 348 % Kroneckered norms and coefficients
347 obj.RHOi_kron = kron(obj.RHOi, I_dim);
348 obj.Hi_kron = kron(obj.Hi, I_dim); 349 obj.Hi_kron = kron(obj.Hi, I_dim);
349 350
350 % Misc. 351 % Misc.
351 obj.m = m; 352 obj.m = m;
352 obj.h = h; 353 obj.h = h;
433 % (Otherwise it's not symmetric.) 434 % (Otherwise it's not symmetric.)
434 for i = dComps 435 for i = dComps
435 Z = sparse(m_tot, m_tot); 436 Z = sparse(m_tot, m_tot);
436 for l = 1:dim 437 for l = 1:dim
437 for k = 1:dim 438 for k = 1:dim
438 Z = Z + nu(l)*C{i,l,k,j}*nu(k); 439 Z = Z + nu(l)*C{l,i,k,j}*nu(k);
439 end 440 end
440 end 441 end
441 Z = -tuning*dim/h11*Z; 442 Z = -tuning*dim/h11*Z;
442 X = Z; 443 X = Z;
443 closure = closure + E{i}*RHOi*Hi*X'*e*H_gamma*(e'*E{j}' ); 444 closure = closure + E{i}*RHOi*Hi*X'*e*H_gamma*(e'*E{j}' );
534 Y = 1/2*T_u{j,i}'; 535 Y = 1/2*T_u{j,i}';
535 Z_u = sparse(m_int, m_int); 536 Z_u = sparse(m_int, m_int);
536 Z_v = sparse(m_int, m_int); 537 Z_v = sparse(m_int, m_int);
537 for l = 1:dim 538 for l = 1:dim
538 for k = 1:dim 539 for k = 1:dim
539 Z_u = Z_u + e_u'*nu_u(l)*C_u{i,l,k,j}*nu_u(k)*e_u; 540 Z_u = Z_u + e_u'*nu_u(l)*C_u{l,i,k,j}*nu_u(k)*e_u;
540 Z_v = Z_v + e_v'*nu_v(l)*C_v{i,l,k,j}*nu_v(k)*e_v; 541 Z_v = Z_v + e_v'*nu_v(l)*C_v{l,i,k,j}*nu_v(k)*e_v;
541 end 542 end
542 end 543 end
543 Z = -tuning*dim*( 1/(4*h11_u)*Z_u + 1/(4*h11_v)*Z_v ); 544 Z = -tuning*dim*( 1/(4*h11_u)*Z_u + 1/(4*h11_v)*Z_v );
544 X = Y + Z*e_u'; 545 X = Y + Z*e_u';
545 closure = closure + E_u{i}*X'*H_gamma*e_u'*E_u{j}'; 546 closure = closure + E_u{i}*X'*H_gamma*e_u'*E_u{j}';