--- a/+scheme/ViscoElastic2d.m	Mon Jul 20 22:21:57 2020 -0700
+++ b/+scheme/ViscoElastic2d.m	Tue Jul 21 21:07:37 2020 -0700
@@ -337,13 +337,75 @@
 
             switch type.type
             case 'standard'
-                [closure, penalty] = obj.refObj.interface(boundary,neighbour_scheme.refObj,neighbour_boundary,type);
+                [closure, penalty] = obj.interfaceStandard(boundary,neighbour_scheme,neighbour_boundary,type);
             case 'frictionalFault'
                 [closure, penalty] = obj.interfaceFrictionalFault(boundary,neighbour_scheme,neighbour_boundary,type);
             end
 
         end
 
+        function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
+
+            % u denotes the solution in the own domain
+            % v denotes the solution in the neighbour domain
+            u = obj;
+            v = neighbour_scheme;
+
+            dim = obj.dim;
+
+            n       = u.getNormal(boundary);
+            H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
+            e       = u.getBoundaryOperatorForScalarField('e', boundary);
+
+            ev      = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
+
+            H       = u.H;
+            RHO     = u.RHO;
+            ETA     = u.ETA;
+            C       = u.C;
+            Eu      = u.Eu;
+            eU      = u.eU;
+            eGamma  = u.eGamma;
+
+            CV       = v.C;
+            Ev      = v.Eu;
+            eV      = v.eU;
+            eGammaV = v.eGamma;
+            nV      = v.getNormal(neighbour_boundary);
+
+
+            % Get elastic closure and penalty
+            [closure, penalty] = obj.elasticObj.interface(boundary, v.elasticObj, neighbour_boundary, type);
+            closure = Eu*closure*Eu';
+            penalty = Eu*penalty*Ev';
+
+            % Add viscous part of traction coupling
+            for i = 1:dim
+                for j = 1:dim
+                    for k = 1:dim
+                        for l = 1:dim
+                            closure = closure + 1/2*eU{j}*( (RHO*H)\(C{i,j,k,l}*e*H_gamma*n{i}*e'*eGamma{k,l}') );
+                            penalty = penalty + 1/2*eU{j}*( (RHO*H)\(e*H_gamma*nV{i}*(ev'*CV{i,j,k,l}*ev)*ev'*eGammaV{k,l}') );
+                        end
+                    end
+                end
+            end
+
+            % Add penalty to strain rate eq
+            for i = 1:dim
+                for j = 1:dim
+                    for k = 1:dim
+                        for l = 1:dim
+                            closure = closure - 1/2*eGamma{i,j}*( (H*ETA)\(C{i,j,k,l}*e*H_gamma*n{k}*e'*eU{l}') );
+                            penalty = penalty + 1/2*eGamma{i,j}*( (H*ETA)\(C{i,j,k,l}*e*H_gamma*n{k}*ev'*eV{l}') );
+                        end
+                    end
+                end
+            end
+
+
+        end
+
         function [closure, penalty] = interfaceFrictionalFault(obj,boundary,neighbour_scheme,neighbour_boundary,type)
             tuning = type.tuning;