diff -r bd8e607768ce -r 34997aced843 +scheme/Elastic2dCurvilinearAnisotropic.m
--- a/+scheme/Elastic2dCurvilinearAnisotropic.m	Sun Jul 26 20:05:10 2020 -0700
+++ b/+scheme/Elastic2dCurvilinearAnisotropic.m	Sun Jul 26 20:06:06 2020 -0700
@@ -495,30 +495,33 @@
         %          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.refObj.interface(boundary,neighbour_scheme.refObj,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);
             case 'frictionalFault'
-                [closure, penalty] = obj.interfaceFrictionalFault(boundary,neighbour_scheme,neighbour_boundary,type);
+                [closure, penalty, forcingPenalties] = obj.interfaceFrictionalFault(boundary,neighbour_scheme,neighbour_boundary,type);
             end
 
         end
 
-        function [closure, penalty] = interfaceFrictionalFault(obj,boundary,neighbour_scheme,neighbour_boundary,type)
+        function [closure, penalty, forcingPenalties] = interfaceFrictionalFault(obj,boundary,neighbour_scheme,neighbour_boundary,type)
             tuning = type.tuning;
 
             % u denotes the solution in the own domain
             % v denotes the solution in the neighbour domain
 
+            forcingPenalties = cell(1, 1);
             u = obj;
             v = neighbour_scheme;
 
@@ -584,25 +587,30 @@
             % Continuity of normal traction
             closure = closure - 1/2*en_u*H_gamma*tau_n_u';
             penalty = penalty + 1/2*en_u*H_gamma*tau_n_v';
+            forcing_tau_n = 1/2*en_u*H_gamma;
 
             % Multiply all normal component terms by inverse of density x quadraure
             closure = RHOi*Hi*closure;
             penalty = RHOi*Hi*penalty;
+            forcing_tau_n = RHOi*Hi*forcing_tau_n;
 
             % ---- Tangential tractions are imposed just like traction BC ------
             closure = closure + obj.boundary_condition(boundary, {'t', 't'});
 
+            forcingPenalties{1} = forcing_tau_n;
+
         end
 
         % Same interface conditions as in interfaceStandard, but imposed in the normal-tangential
         % coordinate system. For the isotropic case, the components decouple nicely.
         % The resulting scheme is not identical to that of interfaceStandard. This appears to be better.
-        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
             % v denotes the solution in the neighbour domain
 
+            forcingPenalties = cell(2, 1);
             u = obj;
             v = neighbour_scheme;
 
@@ -676,10 +684,12 @@
             % Continuity of normal component
             closure = closure + en_u*H_gamma*Zn*en_u';
             penalty = penalty + en_u*H_gamma*Zn*en_v';
+            forcing_u_n = -en_u*H_gamma*Zn;
 
             % Continuity of tangential component
             closure = closure + et_u*H_gamma*Zt*et_u';
             penalty = penalty + et_u*H_gamma*Zt*et_v';
+            forcing_u_t = -et_u*H_gamma*Zt;
             %------------------------------------------------------------------
 
             % --- Continuity of displacement, term 2: The symmetrizing term
@@ -687,6 +697,7 @@
             % Continuity of normal displacement
             closure = closure + 1/2*tau_n_u*H_gamma*en_u';
             penalty = penalty + 1/2*tau_n_u*H_gamma*en_v';
+            forcing_u_n = forcing_u_n - 1/2*tau_n_u*H_gamma;
 
             % Continuity of tangential displacement
             closure = closure + 1/2*tau_t_u*H_gamma*et_u';
@@ -707,6 +718,11 @@
             % Multiply all terms by inverse of density x quadraure
             closure = RHOi*Hi*closure;
             penalty = RHOi*Hi*penalty;
+            forcing_u_n = RHOi*Hi*forcing_u_n;
+            forcing_u_t = RHOi*Hi*forcing_u_t;
+
+            forcingPenalties{1} = forcing_u_n;
+            forcingPenalties{2} = forcing_u_t;
 
         end