Mercurial > repos > public > sbplib

changeset 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 8b1110385ee2
children b40359c9faed
files +scheme/ViscoElastic2d.m
diffstat 1 files changed, 22 insertions(+), 8 deletions(-) [+]
line wrap: on
line diff
diff -r 8b1110385ee2 -r cd0bfcaef0ba +scheme/ViscoElastic2d.m
--- a/+scheme/ViscoElastic2d.m	Tue Jul 28 21:58:22 2020 -0700
+++ b/+scheme/ViscoElastic2d.m	Tue Jul 28 21:59:41 2020 -0700
@@ -382,20 +382,22 @@
         %          Fields:
         %          -- tuning:           penalty strength, defaults to 1.0
         %          -- interpolation:    type of interpolation, default 'none'
-        function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
+        function [closure, penalty, forcingPenalties] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
 
             defaultType.tuning = 1.0;
             defaultType.interpolation = 'none';
             defaultType.type = 'standard';
             default_struct('type', defaultType);
 
+            forcingPenalties = [];
+
             switch type.type
             case 'standard'
                 [closure, penalty] = obj.interfaceStandard(boundary,neighbour_scheme,neighbour_boundary,type);
             case 'frictionalFault'
                 [closure, penalty] = obj.interfaceFrictionalFault(boundary,neighbour_scheme,neighbour_boundary,type);
             case 'normalTangential'
-                [closure, penalty] = obj.interfaceNormalTangential(boundary,neighbour_scheme,neighbour_boundary,type);
+                [closure, penalty, forcingPenalties] = obj.interfaceNormalTangential(boundary,neighbour_scheme,neighbour_boundary,type);
             end
 
         end
@@ -548,7 +550,7 @@
 
         end
 
-        function [closure, penalty] = interfaceNormalTangential(obj,boundary,neighbour_scheme,neighbour_boundary,type)
+        function [closure, penalty, forcingPenalties] = interfaceNormalTangential(obj,boundary,neighbour_scheme,neighbour_boundary,type)
             tuning = type.tuning;
 
             % u denotes the solution in the own domain
@@ -594,10 +596,16 @@
             end
 
             % Get elastic closure and penalty
-            [closure, penalty] = obj.elasticObj.interface(boundary, v.elasticObj, neighbour_boundary, type);
+            [closure, penalty, forcingPenalties] = obj.elasticObj.interface(boundary, v.elasticObj, neighbour_boundary, type);
             closure = Eu*closure*Eu';
             penalty = Eu*penalty*Ev';
 
+            for i = 1:numel(forcingPenalties)
+                forcingPenalties{i} = Eu*forcingPenalties{i};
+            end
+            forcing_u_n = forcingPenalties{1};
+            forcing_u_t = forcingPenalties{2};
+
             % ------ Traction coupling, viscous part -----------
             for i = 1:dim
                 for j = 1:dim
@@ -626,19 +634,25 @@
                         for l = 1:dim
                             for m = 1:dim
                                 % Normal component
-                                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}') );
-                                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}') );
+                                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}') );
+                                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}') );
+
 
                                 % Tangential component
-                                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}') );
-                                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}') );
+                                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}') );
+                                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}') );
                             end
+                            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) );
+                            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) );
                         end
                     end
                 end
             end
             %-------------------------------------------------
 
+            forcingPenalties{1} = forcing_u_n;
+            forcingPenalties{2} = forcing_u_t;
+
         end
 
         % Returns the outward unit normal vector for the boundary specified by the string boundary.