Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dCurvilinearAnisotropic.m @ 1317:34997aced843 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add some interface forcing penalties in ElasticCurvilinearAnisotropic
author Martin Almquist <malmquist@stanford.edu>
date Sun, 26 Jul 2020 20:06:06 -0700
parents 6c308da9dcbc
children 8b1110385ee2
comparison
equal deleted inserted replaced
1316:bd8e607768ce 1317:34997aced843
493 493
494 % type Struct that specifies the interface coupling. 494 % type Struct that specifies the interface coupling.
495 % Fields: 495 % Fields:
496 % -- tuning: penalty strength, defaults to 1.0 496 % -- tuning: penalty strength, defaults to 1.0
497 % -- interpolation: type of interpolation, default 'none' 497 % -- interpolation: type of interpolation, default 'none'
498 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) 498 function [closure, penalty, forcingPenalties] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
499 499
500 defaultType.tuning = 1.0; 500 defaultType.tuning = 1.0;
501 defaultType.interpolation = 'none'; 501 defaultType.interpolation = 'none';
502 defaultType.type = 'standard'; 502 defaultType.type = 'standard';
503 default_struct('type', defaultType); 503 default_struct('type', defaultType);
504 504
505 forcingPenalties = [];
506
505 switch type.type 507 switch type.type
506 case 'standard' 508 case 'standard'
507 [closure, penalty] = obj.refObj.interface(boundary,neighbour_scheme.refObj,neighbour_boundary,type); 509 [closure, penalty] = obj.refObj.interface(boundary,neighbour_scheme.refObj,neighbour_boundary,type);
508 case 'normalTangential' 510 case 'normalTangential'
509 [closure, penalty] = obj.interfaceNormalTangential(boundary,neighbour_scheme,neighbour_boundary,type); 511 [closure, penalty, forcingPenalties] = obj.interfaceNormalTangential(boundary,neighbour_scheme,neighbour_boundary,type);
510 case 'frictionalFault' 512 case 'frictionalFault'
511 [closure, penalty] = obj.interfaceFrictionalFault(boundary,neighbour_scheme,neighbour_boundary,type); 513 [closure, penalty, forcingPenalties] = obj.interfaceFrictionalFault(boundary,neighbour_scheme,neighbour_boundary,type);
512 end 514 end
513 515
514 end 516 end
515 517
516 function [closure, penalty] = interfaceFrictionalFault(obj,boundary,neighbour_scheme,neighbour_boundary,type) 518 function [closure, penalty, forcingPenalties] = interfaceFrictionalFault(obj,boundary,neighbour_scheme,neighbour_boundary,type)
517 tuning = type.tuning; 519 tuning = type.tuning;
518 520
519 % u denotes the solution in the own domain 521 % u denotes the solution in the own domain
520 % v denotes the solution in the neighbour domain 522 % v denotes the solution in the neighbour domain
521 523
524 forcingPenalties = cell(1, 1);
522 u = obj; 525 u = obj;
523 v = neighbour_scheme; 526 v = neighbour_scheme;
524 527
525 % Operators, u side 528 % Operators, u side
526 e_u = u.getBoundaryOperatorForScalarField('e', boundary); 529 e_u = u.getBoundaryOperatorForScalarField('e', boundary);
582 penalty = penalty + 1/2*tau_n_u*H_gamma*en_v'; 585 penalty = penalty + 1/2*tau_n_u*H_gamma*en_v';
583 586
584 % Continuity of normal traction 587 % Continuity of normal traction
585 closure = closure - 1/2*en_u*H_gamma*tau_n_u'; 588 closure = closure - 1/2*en_u*H_gamma*tau_n_u';
586 penalty = penalty + 1/2*en_u*H_gamma*tau_n_v'; 589 penalty = penalty + 1/2*en_u*H_gamma*tau_n_v';
590 forcing_tau_n = 1/2*en_u*H_gamma;
587 591
588 % Multiply all normal component terms by inverse of density x quadraure 592 % Multiply all normal component terms by inverse of density x quadraure
589 closure = RHOi*Hi*closure; 593 closure = RHOi*Hi*closure;
590 penalty = RHOi*Hi*penalty; 594 penalty = RHOi*Hi*penalty;
595 forcing_tau_n = RHOi*Hi*forcing_tau_n;
591 596
592 % ---- Tangential tractions are imposed just like traction BC ------ 597 % ---- Tangential tractions are imposed just like traction BC ------
593 closure = closure + obj.boundary_condition(boundary, {'t', 't'}); 598 closure = closure + obj.boundary_condition(boundary, {'t', 't'});
599
600 forcingPenalties{1} = forcing_tau_n;
594 601
595 end 602 end
596 603
597 % Same interface conditions as in interfaceStandard, but imposed in the normal-tangential 604 % Same interface conditions as in interfaceStandard, but imposed in the normal-tangential
598 % coordinate system. For the isotropic case, the components decouple nicely. 605 % 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. 606 % 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) 607 function [closure, penalty, forcingPenalties] = interfaceNormalTangential(obj,boundary,neighbour_scheme,neighbour_boundary,type)
601 tuning = type.tuning; 608 tuning = type.tuning;
602 609
603 % u denotes the solution in the own domain 610 % u denotes the solution in the own domain
604 % v denotes the solution in the neighbour domain 611 % v denotes the solution in the neighbour domain
605 612
613 forcingPenalties = cell(2, 1);
606 u = obj; 614 u = obj;
607 v = neighbour_scheme; 615 v = neighbour_scheme;
608 616
609 % Operators, u side 617 % Operators, u side
610 e_u = u.getBoundaryOperatorForScalarField('e', boundary); 618 e_u = u.getBoundaryOperatorForScalarField('e', boundary);
674 Zt = -tuning*dim*( 1/(4*h11_u)*s_u*Zt_u + 1/(4*h11_v)*s_v*Zt_v ); 682 Zt = -tuning*dim*( 1/(4*h11_u)*s_u*Zt_u + 1/(4*h11_v)*s_v*Zt_v );
675 683
676 % Continuity of normal component 684 % Continuity of normal component
677 closure = closure + en_u*H_gamma*Zn*en_u'; 685 closure = closure + en_u*H_gamma*Zn*en_u';
678 penalty = penalty + en_u*H_gamma*Zn*en_v'; 686 penalty = penalty + en_u*H_gamma*Zn*en_v';
687 forcing_u_n = -en_u*H_gamma*Zn;
679 688
680 % Continuity of tangential component 689 % Continuity of tangential component
681 closure = closure + et_u*H_gamma*Zt*et_u'; 690 closure = closure + et_u*H_gamma*Zt*et_u';
682 penalty = penalty + et_u*H_gamma*Zt*et_v'; 691 penalty = penalty + et_u*H_gamma*Zt*et_v';
692 forcing_u_t = -et_u*H_gamma*Zt;
683 %------------------------------------------------------------------ 693 %------------------------------------------------------------------
684 694
685 % --- Continuity of displacement, term 2: The symmetrizing term 695 % --- Continuity of displacement, term 2: The symmetrizing term
686 696
687 % Continuity of normal displacement 697 % Continuity of normal displacement
688 closure = closure + 1/2*tau_n_u*H_gamma*en_u'; 698 closure = closure + 1/2*tau_n_u*H_gamma*en_u';
689 penalty = penalty + 1/2*tau_n_u*H_gamma*en_v'; 699 penalty = penalty + 1/2*tau_n_u*H_gamma*en_v';
700 forcing_u_n = forcing_u_n - 1/2*tau_n_u*H_gamma;
690 701
691 % Continuity of tangential displacement 702 % Continuity of tangential displacement
692 closure = closure + 1/2*tau_t_u*H_gamma*et_u'; 703 closure = closure + 1/2*tau_t_u*H_gamma*et_u';
693 penalty = penalty + 1/2*tau_t_u*H_gamma*et_v'; 704 penalty = penalty + 1/2*tau_t_u*H_gamma*et_v';
694 % ------------------------------------------------------------------ 705 % ------------------------------------------------------------------
705 %-------------------------------------------------------------------- 716 %--------------------------------------------------------------------
706 717
707 % Multiply all terms by inverse of density x quadraure 718 % Multiply all terms by inverse of density x quadraure
708 closure = RHOi*Hi*closure; 719 closure = RHOi*Hi*closure;
709 penalty = RHOi*Hi*penalty; 720 penalty = RHOi*Hi*penalty;
721 forcing_u_n = RHOi*Hi*forcing_u_n;
722 forcing_u_t = RHOi*Hi*forcing_u_t;
723
724 forcingPenalties{1} = forcing_u_n;
725 forcingPenalties{2} = forcing_u_t;
710 726
711 end 727 end
712 728
713 729
714 % Returns h11 for the boundary specified by the string boundary. 730 % Returns h11 for the boundary specified by the string boundary.