Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariableAnisotropic.m @ 1204:687515778437 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add anisotropic interface
author Martin Almquist <malmquist@stanford.edu>
date Fri, 06 Sep 2019 14:21:18 -0700
parents 25cadc69a589
children 6b203030fb37 05a01f77d0e3
comparison
equal deleted inserted replaced
1203:25cadc69a589 1204:687515778437
405 405
406 j = comp; 406 j = comp;
407 switch type 407 switch type
408 408
409 % Dirichlet boundary condition 409 % Dirichlet boundary condition
410 % OBS! Cannot yet set one component at a time unless one assumes Displacement for all components
411 case {'D','d','dirichlet','Dirichlet','displacement','Displacement'} 410 case {'D','d','dirichlet','Dirichlet','displacement','Displacement'}
412 411
413 if numel(bc) >= 3 412 if numel(bc) >= 3
414 dComps = bc{3}; 413 dComps = bc{3};
415 else 414 else
456 end 455 end
457 end 456 end
458 457
459 % type Struct that specifies the interface coupling. 458 % type Struct that specifies the interface coupling.
460 % Fields: 459 % Fields:
461 % -- tuning: penalty strength, defaults to 1.2 460 % -- tuning: penalty strength, defaults to 1.0
462 % -- interpolation: type of interpolation, default 'none' 461 % -- interpolation: type of interpolation, default 'none'
463 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) 462 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
464 463
465 defaultType.tuning = 1.2; 464 defaultType.tuning = 1.0;
466 defaultType.interpolation = 'none'; 465 defaultType.interpolation = 'none';
467 default_struct('type', defaultType); 466 default_struct('type', defaultType);
468 467
469 switch type.interpolation 468 switch type.interpolation
470 case {'none', ''} 469 case {'none', ''}
479 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type) 478 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
480 tuning = type.tuning; 479 tuning = type.tuning;
481 480
482 % u denotes the solution in the own domain 481 % u denotes the solution in the own domain
483 % v denotes the solution in the neighbour domain 482 % v denotes the solution in the neighbour domain
484 % Operators without subscripts are from the own domain. 483
485 484 u = obj;
486 % Get boundary operators 485 v = neighbour_scheme;
487 e = obj.getBoundaryOperator('e', boundary); 486
488 tau = obj.getBoundaryOperator('tau', boundary); 487 % Operators, u side
489 488 e_u = u.getBoundaryOperatorForScalarField('e', boundary);
490 e_v = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary); 489 tau_u = u.getBoundaryOperator('tau', boundary);
491 tau_v = neighbour_scheme.getBoundaryOperator('tau', neighbour_boundary); 490 h11_u = u.getBorrowing(boundary);
492 491 nu_u = u.getNormal(boundary);
493 H_gamma = obj.getBoundaryQuadrature(boundary); 492
494 493 E_u = u.E;
495 % Operators and quantities that correspond to the own domain only 494 C_u = u.C;
496 Hi = obj.Hi_kron; 495 m_tot_u = u.grid.N();
497 RHOi = obj.RHOi_kron; 496
498 497 % Operators, v side
499 % Penalty strength operators 498 e_v = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
500 alpha_u = 1/4*tuning*obj.getBoundaryOperator('alpha', boundary); 499 tau_v = v.getBoundaryOperator('tau', neighbour_boundary);
501 alpha_v = 1/4*tuning*neighbour_scheme.getBoundaryOperator('alpha', neighbour_boundary); 500 h11_v = v.getBorrowing(neighbour_boundary);
502 501 nu_v = v.getNormal(neighbour_boundary);
503 closure = -RHOi*Hi*e*H_gamma*(alpha_u' + alpha_v'*e_v*e'); 502
504 penalty = RHOi*Hi*e*H_gamma*(alpha_u'*e*e_v' + alpha_v'); 503 E_v = v.E;
505 504 C_v = v.C;
506 closure = closure - 1/2*RHOi*Hi*e*H_gamma*tau'; 505 m_tot_v = v.grid.N();
507 penalty = penalty - 1/2*RHOi*Hi*e*H_gamma*tau_v'; 506
508 507 % Operators that are only required for own domain
509 closure = closure + 1/2*RHOi*Hi*tau*H_gamma*e'; 508 Hi = u.Hi_kron;
510 penalty = penalty - 1/2*RHOi*Hi*tau*H_gamma*e_v'; 509 RHOi = u.RHOi_kron;
510 e_kron = u.getBoundaryOperator('e', boundary);
511 T_u = u.getBoundaryTractionOperator(boundary);
512
513 % Shared operators
514 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
515 H_gamma_kron = u.getBoundaryQuadrature(boundary);
516 dim = u.dim;
517
518 % Preallocate
519 [~, m_int] = size(H_gamma);
520 closure = sparse(dim*m_tot_u, dim*m_tot_u);
521 penalty = sparse(dim*m_tot_u, dim*m_tot_v);
522
523 % ---- Continuity of displacement ------
524
525 % Y: symmetrizing part of penalty
526 % Z: symmetric part of penalty
527 % X = Y + Z.
528
529 % Loop over components to couple across interface
530 for j = 1:dim
531
532 % Loop over components that penalties end up on
533 for i = 1:dim
534 Y = 1/2*T_u{j,i}';
535 Z_u = sparse(m_int, m_int);
536 Z_v = sparse(m_int, m_int);
537 for l = 1:dim
538 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_v = Z_v + e_v'*nu_v(l)*C_v{i,l,k,j}*nu_v(k)*e_v;
541 end
542 end
543 Z = -tuning*dim*( 1/(4*h11_u)*Z_u + 1/(4*h11_v)*Z_v );
544 X = Y + Z*e_u';
545 closure = closure + E_u{i}*X'*H_gamma*e_u'*E_u{j}';
546 penalty = penalty - E_u{i}*X'*H_gamma*e_v'*E_v{j}';
547 end
548 end
549
550 % ---- Continuity of traction ------
551 closure = closure - 1/2*e_kron*H_gamma_kron*tau_u';
552 penalty = penalty - 1/2*e_kron*H_gamma_kron*tau_v';
553
554 % ---- Multiply by inverse of density x quadraure ----
555 closure = RHOi*Hi*closure;
556 penalty = RHOi*Hi*penalty;
511 557
512 end 558 end
513 559
514 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type) 560 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
515 error('Non-conforming interfaces not implemented yet.'); 561 error('Non-conforming interfaces not implemented yet.');