Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariable.m @ 958:72cd29107a9a feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
author Martin Almquist <malmquist@stanford.edu>
date Wed, 05 Dec 2018 18:58:10 -0800
parents e30aaa4a3e09
children c226fb8c2b8a
comparison
equal deleted inserted replaced
919:e30aaa4a3e09 958:72cd29107a9a
236 T_l{j} = cell(dim,dim); 236 T_l{j} = cell(dim,dim);
237 T_r{j} = cell(dim,dim); 237 T_r{j} = cell(dim,dim);
238 tau_l{j} = cell(dim,1); 238 tau_l{j} = cell(dim,1);
239 tau_r{j} = cell(dim,1); 239 tau_r{j} = cell(dim,1);
240 240
241 LAMBDA_l = e_l{j}'*LAMBDA*e_l{j};
242 LAMBDA_r = e_r{j}'*LAMBDA*e_r{j};
243 MU_l = e_l{j}'*MU*e_l{j};
244 MU_r = e_r{j}'*MU*e_r{j};
245
246 [~, n_l] = size(e_l{j});
247 [~, n_r] = size(e_r{j});
248
241 % Loop over components 249 % Loop over components
242 for i = 1:dim 250 for i = 1:dim
243 tau_l{j}{i} = sparse(m_tot,dim*m_tot); 251 tau_l{j}{i} = sparse(n_l, dim*m_tot);
244 tau_r{j}{i} = sparse(m_tot,dim*m_tot); 252 tau_r{j}{i} = sparse(n_r, dim*m_tot);
245 for k = 1:dim 253 for k = 1:dim
246 T_l{j}{i,k} = ... 254 T_l{j}{i,k} = ...
247 -d(i,j)*LAMBDA*(d(i,k)*e_l{k}*d1_l{k}' + db(i,k)*D1{k})... 255 -d(i,j)*LAMBDA_l*(d(i,k)*d1_l{k}' + db(i,k)*e_l{j}'*D1{k})...
248 -d(j,k)*MU*(d(i,j)*e_l{i}*d1_l{i}' + db(i,j)*D1{i})... 256 -d(j,k)*MU_l*(d(i,j)*d1_l{i}' + db(i,j)*e_l{j}'*D1{i})...
249 -d(i,k)*MU*e_l{j}*d1_l{j}'; 257 -d(i,k)*MU_l*d1_l{j}';
250 258
251 T_r{j}{i,k} = ... 259 T_r{j}{i,k} = ...
252 d(i,j)*LAMBDA*(d(i,k)*e_r{k}*d1_r{k}' + db(i,k)*D1{k})... 260 d(i,j)*LAMBDA_r*(d(i,k)*d1_r{k}' + db(i,k)*e_r{j}'*D1{k})...
253 +d(j,k)*MU*(d(i,j)*e_r{i}*d1_r{i}' + db(i,j)*D1{i})... 261 +d(j,k)*MU_r*(d(i,j)*d1_r{i}' + db(i,j)*e_r{j}'*D1{i})...
254 +d(i,k)*MU*e_r{j}*d1_r{j}'; 262 +d(i,k)*MU_r*d1_r{j}';
255 263
256 tau_l{j}{i} = tau_l{j}{i} + T_l{j}{i,k}*E{k}'; 264 tau_l{j}{i} = tau_l{j}{i} + T_l{j}{i,k}*E{k}';
257 tau_r{j}{i} = tau_r{j}{i} + T_r{j}{i,k}*E{k}'; 265 tau_r{j}{i} = tau_r{j}{i} + T_r{j}{i,k}*E{k}';
258 end 266 end
259 267
260 end 268 end
261 end 269 end
270
271 % Transpose T and tau to match boundary operator convention
272 for i = 1:dim
273 for j = 1:dim
274 tau_l{i}{j} = transpose(tau_l{i}{j});
275 tau_r{i}{j} = transpose(tau_r{i}{j});
276 for k = 1:dim
277 T_l{i}{j,k} = transpose(T_l{i}{j,k});
278 T_r{i}{j,k} = transpose(T_r{i}{j,k});
279 end
280 end
281 end
282
262 obj.T_l = T_l; 283 obj.T_l = T_l;
263 obj.T_r = T_r; 284 obj.T_r = T_r;
264 obj.tau_l = tau_l; 285 obj.tau_l = tau_l;
265 obj.tau_r = tau_r; 286 obj.tau_r = tau_r;
266 287
342 + d(i,j)* a_mu_i*MU ... 363 + d(i,j)* a_mu_i*MU ...
343 + db(i,j)*a_mu_ij*MU ); 364 + db(i,j)*a_mu_ij*MU );
344 365
345 % Loop over components that Dirichlet penalties end up on 366 % Loop over components that Dirichlet penalties end up on
346 for i = 1:dim 367 for i = 1:dim
347 C = T{k,i}; 368 C = transpose(T{k,i});
348 A = -d(i,k)*alpha(i,j); 369 A = -d(i,k)*alpha(i,j);
349 B = A + C; 370 B = A + e*C;
350 closure = closure + E{i}*RHOi*Hi*B'*e*H_gamma*(e'*E{k}' ); 371 closure = closure + E{i}*RHOi*Hi*B'*e*H_gamma*(e'*E{k}' );
351 penalty = penalty - E{i}*RHOi*Hi*B'*e*H_gamma; 372 penalty = penalty - E{i}*RHOi*Hi*B'*e*H_gamma;
352 end 373 end
353 374
354 % Free boundary condition 375 % Free boundary condition
355 case {'F','f','Free','free','traction','Traction','t','T'} 376 case {'F','f','Free','free','traction','Traction','t','T'}
356 closure = closure - E{k}*RHOi*Hi*e*H_gamma* (e'*tau{k} ); 377 closure = closure - E{k}*RHOi*Hi*e*H_gamma*tau{k}';
357 penalty = penalty + E{k}*RHOi*Hi*e*H_gamma; 378 penalty = penalty + E{k}*RHOi*Hi*e*H_gamma;
358 379
359 % Unknown boundary condition 380 % Unknown boundary condition
360 otherwise 381 otherwise
361 error('No such boundary condition: type = %s',type); 382 error('No such boundary condition: type = %s',type);
432 closure = closure - 1/2*E{i}*RHOi*Hi*e*H_gamma*e'*tau{i}; 453 closure = closure - 1/2*E{i}*RHOi*Hi*e*H_gamma*e'*tau{i};
433 penalty = penalty - 1/2*E{i}*RHOi*Hi*e*H_gamma*e_v'*tau_v{i}; 454 penalty = penalty - 1/2*E{i}*RHOi*Hi*e*H_gamma*e_v'*tau_v{i};
434 455
435 % Loop over components that we have interface conditions on 456 % Loop over components that we have interface conditions on
436 for k = 1:dim 457 for k = 1:dim
437 closure = closure + 1/2*E{i}*RHOi*Hi*T{k,i}'*e*H_gamma*e'*E{k}'; 458 closure = closure + 1/2*E{i}*RHOi*Hi*T{k,i}*H_gamma*e'*E{k}';
438 penalty = penalty - 1/2*E{i}*RHOi*Hi*T{k,i}'*e*H_gamma*e_v'*E_v{k}'; 459 penalty = penalty - 1/2*E{i}*RHOi*Hi*T{k,i}*H_gamma*e_v'*E_v{k}';
439 end 460 end
440 end 461 end
441 end 462 end
442 463
443 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. 464 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
475 496
476 if ~iscell(op) 497 if ~iscell(op)
477 op = {op}; 498 op = {op};
478 end 499 end
479 500
480 for i = 1:length(op) 501 for k = 1:length(op)
481 switch op{i} 502 switch op{k}
482 case 'e' 503 case 'e'
483 switch boundary 504 switch boundary
484 case {'w','W','west','West','s','S','south','South'} 505 case {'w','W','west','West','s','S','south','South'}
485 varargout{i} = obj.e_l{j}; 506 varargout{k} = obj.e_l{j};
486 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 507 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
487 varargout{i} = obj.e_r{j}; 508 varargout{k} = obj.e_r{j};
488 end 509 end
489 case 'd' 510 case 'd'
490 switch boundary 511 switch boundary
491 case {'w','W','west','West','s','S','south','South'} 512 case {'w','W','west','West','s','S','south','South'}
492 varargout{i} = obj.d1_l{j}; 513 varargout{k} = obj.d1_l{j};
493 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 514 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
494 varargout{i} = obj.d1_r{j}; 515 varargout{k} = obj.d1_r{j};
495 end 516 end
496 case 'H' 517 case 'H'
497 varargout{i} = obj.H_boundary{j}; 518 varargout{k} = obj.H_boundary{j};
498 case 'T' 519 case 'T'
499 switch boundary 520 switch boundary
500 case {'w','W','west','West','s','S','south','South'} 521 case {'w','W','west','West','s','S','south','South'}
501 varargout{i} = obj.T_l{j}; 522 varargout{k} = obj.T_l{j};
502 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 523 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
503 varargout{i} = obj.T_r{j}; 524 varargout{k} = obj.T_r{j};
504 end 525 end
505 case 'tau' 526 case 'tau'
506 switch boundary 527 switch boundary
507 case {'w','W','west','West','s','S','south','South'} 528 case {'w','W','west','West','s','S','south','South'}
508 varargout{i} = obj.tau_l{j}; 529 varargout{k} = obj.tau_l{j};
509 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 530 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
510 varargout{i} = obj.tau_r{j}; 531 varargout{k} = obj.tau_r{j};
511 end 532 end
512 case 'alpha' 533 case 'alpha'
513 % alpha = alpha(i,j) is the penalty strength for displacement BC. 534 % alpha = alpha(i,j) is the penalty strength for displacement BC.
514 tuning = 1.2; 535 tuning = 1.2;
515 LAMBDA = obj.LAMBDA; 536 LAMBDA = obj.LAMBDA;
524 a_mu_i = 2/theta_M; 545 a_mu_i = 2/theta_M;
525 a_mu_ij = 2/theta_H + 1/theta_R; 546 a_mu_ij = 2/theta_H + 1/theta_R;
526 547
527 d = @kroneckerDelta; % Kronecker delta 548 d = @kroneckerDelta; % Kronecker delta
528 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta 549 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
529 alpha = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ... 550 alpha = cell(obj.dim, obj.dim);
551
552 for i = 1:obj.dim
553 for j = 1:obj.dim
554 alpha{i,j} = tuning*( d(i,j)* a_lambda*LAMBDA ...
530 + d(i,j)* a_mu_i*MU ... 555 + d(i,j)* a_mu_i*MU ...
531 + db(i,j)*a_mu_ij*MU ); 556 + db(i,j)*a_mu_ij*MU );
532 557 end
533 varargout{i} = alpha; 558 end
559
560 varargout{k} = alpha;
534 otherwise 561 otherwise
535 error(['No such operator: operator = ' op{i}]); 562 error(['No such operator: operator = ' op{k}]);
536 end 563 end
537 end 564 end
538 565
539 end 566 end
540 567