Mercurial > repos > public > sbplib

comparison +scheme/ViscoElastic2d.m @ 1320:cd0bfcaef0ba feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add some interface forcing penalties in ViscoElastic2d
author Martin Almquist <malmquist@stanford.edu>
date Tue, 28 Jul 2020 21:59:41 -0700
parents 6650b111742d
children
comparison
equal deleted inserted replaced
1319:8b1110385ee2 1320:cd0bfcaef0ba
380 380
381 % type Struct that specifies the interface coupling. 381 % type Struct that specifies the interface coupling.
382 % Fields: 382 % Fields:
383 % -- tuning: penalty strength, defaults to 1.0 383 % -- tuning: penalty strength, defaults to 1.0
384 % -- interpolation: type of interpolation, default 'none' 384 % -- interpolation: type of interpolation, default 'none'
385 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) 385 function [closure, penalty, forcingPenalties] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
386 386
387 defaultType.tuning = 1.0; 387 defaultType.tuning = 1.0;
388 defaultType.interpolation = 'none'; 388 defaultType.interpolation = 'none';
389 defaultType.type = 'standard'; 389 defaultType.type = 'standard';
390 default_struct('type', defaultType); 390 default_struct('type', defaultType);
391
392 forcingPenalties = [];
391 393
392 switch type.type 394 switch type.type
393 case 'standard' 395 case 'standard'
394 [closure, penalty] = obj.interfaceStandard(boundary,neighbour_scheme,neighbour_boundary,type); 396 [closure, penalty] = obj.interfaceStandard(boundary,neighbour_scheme,neighbour_boundary,type);
395 case 'frictionalFault' 397 case 'frictionalFault'
396 [closure, penalty] = obj.interfaceFrictionalFault(boundary,neighbour_scheme,neighbour_boundary,type); 398 [closure, penalty] = obj.interfaceFrictionalFault(boundary,neighbour_scheme,neighbour_boundary,type);
397 case 'normalTangential' 399 case 'normalTangential'
398 [closure, penalty] = obj.interfaceNormalTangential(boundary,neighbour_scheme,neighbour_boundary,type); 400 [closure, penalty, forcingPenalties] = obj.interfaceNormalTangential(boundary,neighbour_scheme,neighbour_boundary,type);
399 end 401 end
400 402
401 end 403 end
402 404
403 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type) 405 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
546 end 548 end
547 %------------------------------------------------- 549 %-------------------------------------------------
548 550
549 end 551 end
550 552
551 function [closure, penalty] = interfaceNormalTangential(obj,boundary,neighbour_scheme,neighbour_boundary,type) 553 function [closure, penalty, forcingPenalties] = interfaceNormalTangential(obj,boundary,neighbour_scheme,neighbour_boundary,type)
552 tuning = type.tuning; 554 tuning = type.tuning;
553 555
554 % u denotes the solution in the own domain 556 % u denotes the solution in the own domain
555 % v denotes the solution in the neighbour domain 557 % v denotes the solution in the neighbour domain
556 u = obj; 558 u = obj;
592 end 594 end
593 end 595 end
594 end 596 end
595 597
596 % Get elastic closure and penalty 598 % Get elastic closure and penalty
597 [closure, penalty] = obj.elasticObj.interface(boundary, v.elasticObj, neighbour_boundary, type); 599 [closure, penalty, forcingPenalties] = obj.elasticObj.interface(boundary, v.elasticObj, neighbour_boundary, type);
598 closure = Eu*closure*Eu'; 600 closure = Eu*closure*Eu';
599 penalty = Eu*penalty*Ev'; 601 penalty = Eu*penalty*Ev';
602
603 for i = 1:numel(forcingPenalties)
604 forcingPenalties{i} = Eu*forcingPenalties{i};
605 end
606 forcing_u_n = forcingPenalties{1};
607 forcing_u_t = forcingPenalties{2};
600 608
601 % ------ Traction coupling, viscous part ----------- 609 % ------ Traction coupling, viscous part -----------
602 for i = 1:dim 610 for i = 1:dim
603 for j = 1:dim 611 for j = 1:dim
604 for k = 1:dim 612 for k = 1:dim
624 for j = 1:dim 632 for j = 1:dim
625 for k = 1:dim 633 for k = 1:dim
626 for l = 1:dim 634 for l = 1:dim
627 for m = 1:dim 635 for m = 1:dim
628 % Normal component 636 % 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}') ); 637 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}') ); 638 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}') );
639
631 640
632 % Tangential component 641 % 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}') ); 642 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}') ); 643 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 644 end
645 forcing_u_n = forcing_u_n + 1/2*eGamma{i,j}*( (H*ETA)\(e*n{l}*n{k}*C{i,j,k,l}*H_gamma) );
646 forcing_u_t = forcing_u_t + 1/2*eGamma{i,j}*( (H*ETA)\(e*t{l}*n{k}*C{i,j,k,l}*H_gamma) );
636 end 647 end
637 end 648 end
638 end 649 end
639 end 650 end
640 %------------------------------------------------- 651 %-------------------------------------------------
652
653 forcingPenalties{1} = forcing_u_n;
654 forcingPenalties{2} = forcing_u_t;
641 655
642 end 656 end
643 657
644 % Returns the outward unit normal vector for the boundary specified by the string boundary. 658 % Returns the outward unit normal vector for the boundary specified by the string boundary.
645 % n is a cell of diagonal matrices for each normal component, n{1} = n_1, n{2} = n_2. 659 % n is a cell of diagonal matrices for each normal component, n{1} = n_1, n{2} = n_2.