Mercurial > repos > public > sbplib

comparison +scheme/ViscoElastic2d.m @ 1318:6650b111742d feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Implement normal-tangential interface treatment in viscoElastic2d
author Martin Almquist <malmquist@stanford.edu>
date Sun, 26 Jul 2020 20:42:06 -0700
parents 58df4a35fe43
children cd0bfcaef0ba
comparison
equal deleted inserted replaced
1317:34997aced843 1318:6650b111742d
392 switch type.type 392 switch type.type
393 case 'standard' 393 case 'standard'
394 [closure, penalty] = obj.interfaceStandard(boundary,neighbour_scheme,neighbour_boundary,type); 394 [closure, penalty] = obj.interfaceStandard(boundary,neighbour_scheme,neighbour_boundary,type);
395 case 'frictionalFault' 395 case 'frictionalFault'
396 [closure, penalty] = obj.interfaceFrictionalFault(boundary,neighbour_scheme,neighbour_boundary,type); 396 [closure, penalty] = obj.interfaceFrictionalFault(boundary,neighbour_scheme,neighbour_boundary,type);
397 case 'normalTangential'
398 [closure, penalty] = obj.interfaceNormalTangential(boundary,neighbour_scheme,neighbour_boundary,type);
397 end 399 end
398 400
399 end 401 end
400 402
401 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type) 403 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
544 end 546 end
545 %------------------------------------------------- 547 %-------------------------------------------------
546 548
547 end 549 end
548 550
551 function [closure, penalty] = interfaceNormalTangential(obj,boundary,neighbour_scheme,neighbour_boundary,type)
552 tuning = type.tuning;
553
554 % u denotes the solution in the own domain
555 % v denotes the solution in the neighbour domain
556 u = obj;
557 v = neighbour_scheme;
558
559 dim = obj.dim;
560
561 n = u.getNormal(boundary);
562 t = u.getTangent(boundary);
563 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
564 e = u.getBoundaryOperatorForScalarField('e', boundary);
565
566 ev = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
567
568 H = u.H;
569 RHO = u.RHO;
570 ETA = u.ETA;
571 C = u.C;
572 Eu = u.Eu;
573 eU = u.eU;
574 eGamma = u.eGamma;
575 Egamma = u.Egamma;
576
577 CV = v.C;
578 Ev = v.Eu;
579 eV = v.eU;
580 eGammaV = v.eGamma;
581 nV = v.getNormal(neighbour_boundary);
582 tV = v.getTangent(neighbour_boundary);
583
584 % Reduce stiffness tensors to boundary size
585 for i = 1:dim
586 for j = 1:dim
587 for k = 1:dim
588 for l = 1:dim
589 C{i,j,k,l} = e'*C{i,j,k,l}*e;
590 CV{i,j,k,l} = ev'*CV{i,j,k,l}*ev;
591 end
592 end
593 end
594 end
595
596 % Get elastic closure and penalty
597 [closure, penalty] = obj.elasticObj.interface(boundary, v.elasticObj, neighbour_boundary, type);
598 closure = Eu*closure*Eu';
599 penalty = Eu*penalty*Ev';
600
601 % ------ Traction coupling, viscous part -----------
602 for i = 1:dim
603 for j = 1:dim
604 for k = 1:dim
605 for l = 1:dim
606 for m = 1:dim
607 % Normal component
608 closure = closure + 1/2*eU{m}*( (RHO*H)\(e*n{m}*H_gamma*n{j}*n{i}*C{i,j,k,l}*e'*eGamma{k,l}') );
609 penalty = penalty - 1/2*eU{m}*( (RHO*H)\(e*n{m}*H_gamma*nV{j}*nV{i}*CV{i,j,k,l}*ev'*eGammaV{k,l}') );
610
611 % Tangential component
612 closure = closure + 1/2*eU{m}*( (RHO*H)\(e*t{m}*H_gamma*t{j}*n{i}*C{i,j,k,l}*e'*eGamma{k,l}') );
613 penalty = penalty - 1/2*eU{m}*( (RHO*H)\(e*t{m}*H_gamma*tV{j}*nV{i}*CV{i,j,k,l}*ev'*eGammaV{k,l}') );
614 end
615 end
616 end
617 end
618 end
619 %-------------------------------------------------
620
621 % --- Displacement coupling ----------------------
622 % Add penalty to strain rate eq
623 for i = 1:dim
624 for j = 1:dim
625 for k = 1:dim
626 for l = 1:dim
627 for m = 1:dim
628 % Normal component
629 closure = closure - 1/2*eGamma{i,j}*( (H*ETA)\(e*n{l}*n{k}*C{i,j,k,l}*H_gamma*n{m}*e'*eU{m}') );
630 penalty = penalty - 1/2*eGamma{i,j}*( (H*ETA)\(e*n{l}*n{k}*C{i,j,k,l}*H_gamma*nV{m}*ev'*eV{m}') );
631
632 % Tangential component
633 closure = closure - 1/2*eGamma{i,j}*( (H*ETA)\(e*t{l}*n{k}*C{i,j,k,l}*H_gamma*t{m}*e'*eU{m}') );
634 penalty = penalty - 1/2*eGamma{i,j}*( (H*ETA)\(e*t{l}*n{k}*C{i,j,k,l}*H_gamma*tV{m}*ev'*eV{m}') );
635 end
636 end
637 end
638 end
639 end
640 %-------------------------------------------------
641
642 end
643
549 % Returns the outward unit normal vector for the boundary specified by the string boundary. 644 % Returns the outward unit normal vector for the boundary specified by the string boundary.
550 % n is a cell of diagonal matrices for each normal component, n{1} = n_1, n{2} = n_2. 645 % n is a cell of diagonal matrices for each normal component, n{1} = n_1, n{2} = n_2.
551 function n = getNormal(obj, boundary) 646 function n = getNormal(obj, boundary)
552 assertIsMember(boundary, {'w', 'e', 's', 'n'}) 647 assertIsMember(boundary, {'w', 'e', 's', 'n'})
553 648