Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dCurvilinearAnisotropic.m @ 1315:6c308da9dcbc feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add new interface method interfaceNormalTangential in Elastic2dCurvilinearAnisotropic.
author Martin Almquist <malmquist@stanford.edu>
date Sun, 26 Jul 2020 17:38:28 -0700
parents 633757e582e5
children 34997aced843
comparison
equal deleted inserted replaced
1314:58df4a35fe43 1315:6c308da9dcbc
484 484
485 % Term 2: The symmetrizing term 485 % Term 2: The symmetrizing term
486 closure = closure + tau_n*H_gamma*en'; 486 closure = closure + tau_n*H_gamma*en';
487 penalty = penalty - tau_n*H_gamma; 487 penalty = penalty - tau_n*H_gamma;
488 488
489 % Multiply all normal component terms by inverse of density x quadraure 489 % Multiply all terms by inverse of density x quadraure
490 closure = RHOi*Hi*closure; 490 closure = RHOi*Hi*closure;
491 penalty = RHOi*Hi*penalty; 491 penalty = RHOi*Hi*penalty;
492 end 492 end
493 493
494 % type Struct that specifies the interface coupling. 494 % type Struct that specifies the interface coupling.
503 default_struct('type', defaultType); 503 default_struct('type', defaultType);
504 504
505 switch type.type 505 switch type.type
506 case 'standard' 506 case 'standard'
507 [closure, penalty] = obj.refObj.interface(boundary,neighbour_scheme.refObj,neighbour_boundary,type); 507 [closure, penalty] = obj.refObj.interface(boundary,neighbour_scheme.refObj,neighbour_boundary,type);
508 case 'normalTangential'
509 [closure, penalty] = obj.interfaceNormalTangential(boundary,neighbour_scheme,neighbour_boundary,type);
508 case 'frictionalFault' 510 case 'frictionalFault'
509 [closure, penalty] = obj.interfaceFrictionalFault(boundary,neighbour_scheme,neighbour_boundary,type); 511 [closure, penalty] = obj.interfaceFrictionalFault(boundary,neighbour_scheme,neighbour_boundary,type);
510 end 512 end
511 513
512 end 514 end
590 % ---- Tangential tractions are imposed just like traction BC ------ 592 % ---- Tangential tractions are imposed just like traction BC ------
591 closure = closure + obj.boundary_condition(boundary, {'t', 't'}); 593 closure = closure + obj.boundary_condition(boundary, {'t', 't'});
592 594
593 end 595 end
594 596
597 % Same interface conditions as in interfaceStandard, but imposed in the normal-tangential
598 % coordinate system. For the isotropic case, the components decouple nicely.
599 % The resulting scheme is not identical to that of interfaceStandard. This appears to be better.
600 function [closure, penalty] = interfaceNormalTangential(obj,boundary,neighbour_scheme,neighbour_boundary,type)
601 tuning = type.tuning;
602
603 % u denotes the solution in the own domain
604 % v denotes the solution in the neighbour domain
605
606 u = obj;
607 v = neighbour_scheme;
608
609 % Operators, u side
610 e_u = u.getBoundaryOperatorForScalarField('e', boundary);
611 en_u = u.getBoundaryOperator('en', boundary);
612 et_u = u.getBoundaryOperator('et', boundary);
613 tau_n_u = u.getBoundaryOperator('tau_n', boundary);
614 tau_t_u = u.getBoundaryOperator('tau_t', boundary);
615 h11_u = u.getBorrowing(boundary);
616 n_u = u.getNormal(boundary);
617 t_u = u.getTangent(boundary);
618
619 C_u = u.C;
620 Ji_u = u.Ji;
621 s_u = spdiag(u.(['s_' boundary]));
622 m_tot_u = u.grid.N();
623
624 % Operators, v side
625 e_v = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
626 en_v = v.getBoundaryOperator('en', neighbour_boundary);
627 et_v = v.getBoundaryOperator('et', neighbour_boundary);
628 tau_n_v = v.getBoundaryOperator('tau_n', neighbour_boundary);
629 tau_t_v = v.getBoundaryOperator('tau_t', neighbour_boundary);
630 h11_v = v.getBorrowing(neighbour_boundary);
631 n_v = v.getNormal(neighbour_boundary);
632 t_v = v.getTangent(neighbour_boundary);
633
634 C_v = v.C;
635 Ji_v = v.Ji;
636 s_v = spdiag(v.(['s_' neighbour_boundary]));
637 m_tot_v = v.grid.N();
638
639 % Operators that are only required for own domain
640 Hi = u.E{1}*inv(u.H)*u.E{1}' + u.E{2}*inv(u.H)*u.E{2}';
641 RHOi = u.E{1}*inv(u.RHO)*u.E{1}' + u.E{2}*inv(u.RHO)*u.E{2}';
642
643 % Shared operators
644 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
645 dim = u.dim;
646
647 % Preallocate
648 [~, m_int] = size(H_gamma);
649 closure = sparse(dim*m_tot_u, dim*m_tot_u);
650 penalty = sparse(dim*m_tot_u, dim*m_tot_v);
651
652 % -- Continuity of displacement, term 1: The symmetric term ---
653 Zn_u = sparse(m_int, m_int);
654 Zn_v = sparse(m_int, m_int);
655 Zt_u = sparse(m_int, m_int);
656 Zt_v = sparse(m_int, m_int);
657 for i = 1:dim
658 for j = 1:dim
659 for l = 1:dim
660 for k = 1:dim
661 % Penalty strength for normal component
662 Zn_u = Zn_u + n_u{i}*n_u{j}*e_u'*Ji_u*C_u{i,j,k,l}*e_u*n_u{k}*n_u{l};
663 Zn_v = Zn_v + n_v{i}*n_v{j}*e_v'*Ji_v*C_v{i,j,k,l}*e_v*n_v{k}*n_v{l};
664
665 % Penalty strength for tangential component
666 Zt_u = Zt_u + n_u{i}*t_u{j}*e_u'*Ji_u*C_u{i,j,k,l}*e_u*n_u{k}*t_u{l};
667 Zt_v = Zt_v + n_v{i}*t_v{j}*e_v'*Ji_v*C_v{i,j,k,l}*e_v*n_v{k}*t_v{l};
668 end
669 end
670 end
671 end
672
673 Zn = -tuning*dim*( 1/(4*h11_u)*s_u*Zn_u + 1/(4*h11_v)*s_v*Zn_v );
674 Zt = -tuning*dim*( 1/(4*h11_u)*s_u*Zt_u + 1/(4*h11_v)*s_v*Zt_v );
675
676 % Continuity of normal component
677 closure = closure + en_u*H_gamma*Zn*en_u';
678 penalty = penalty + en_u*H_gamma*Zn*en_v';
679
680 % Continuity of tangential component
681 closure = closure + et_u*H_gamma*Zt*et_u';
682 penalty = penalty + et_u*H_gamma*Zt*et_v';
683 %------------------------------------------------------------------
684
685 % --- Continuity of displacement, term 2: The symmetrizing term
686
687 % Continuity of normal displacement
688 closure = closure + 1/2*tau_n_u*H_gamma*en_u';
689 penalty = penalty + 1/2*tau_n_u*H_gamma*en_v';
690
691 % Continuity of tangential displacement
692 closure = closure + 1/2*tau_t_u*H_gamma*et_u';
693 penalty = penalty + 1/2*tau_t_u*H_gamma*et_v';
694 % ------------------------------------------------------------------
695
696 % --- Continuity of tractions -----------------------------
697
698 % Continuity of normal traction
699 closure = closure - 1/2*en_u*H_gamma*tau_n_u';
700 penalty = penalty + 1/2*en_u*H_gamma*tau_n_v';
701
702 % Continuity of tangential traction
703 closure = closure - 1/2*et_u*H_gamma*tau_t_u';
704 penalty = penalty + 1/2*et_u*H_gamma*tau_t_v';
705 %--------------------------------------------------------------------
706
707 % Multiply all terms by inverse of density x quadraure
708 closure = RHOi*Hi*closure;
709 penalty = RHOi*Hi*penalty;
710
711 end
712
595 713
596 % Returns h11 for the boundary specified by the string boundary. 714 % Returns h11 for the boundary specified by the string boundary.
597 % op -- string 715 % op -- string
598 function h11 = getBorrowing(obj, boundary) 716 function h11 = getBorrowing(obj, boundary)
599 assertIsMember(boundary, {'w', 'e', 's', 'n'}) 717 assertIsMember(boundary, {'w', 'e', 's', 'n'})