--- a/+scheme/ViscoElastic2d.m	Tue Jul 21 21:07:37 2020 -0700
+++ b/+scheme/ViscoElastic2d.m	Wed Jul 22 15:35:25 2020 -0700
@@ -411,94 +411,85 @@
 
             % u denotes the solution in the own domain
             % v denotes the solution in the neighbour domain
-
             u = obj;
             v = neighbour_scheme;
 
-            % Operators, u side
-            e_u       = u.getBoundaryOperatorForScalarField('e', boundary);
-            en_u       = u.getBoundaryOperator('en', boundary);
-            tau_n_u     = u.getBoundaryOperator('tau_n', boundary);
-            h11_u     = u.getBorrowing(boundary);
-            n_u      = u.getNormal(boundary);
+            dim = obj.dim;
 
-            C_u = u.C;
-            Ji_u = u.Ji;
-            s_u = spdiag(u.(['s_' boundary]));
-            m_tot_u = u.grid.N();
+            n       = u.getNormal(boundary);
+            H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
+            e       = u.getBoundaryOperatorForScalarField('e', boundary);
 
-            % Operators, v side
-            e_v       = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
-            en_v       = v.getBoundaryOperator('en', neighbour_boundary);
-            tau_n_v     = v.getBoundaryOperator('tau_n', neighbour_boundary);
-            h11_v     = v.getBorrowing(neighbour_boundary);
-            n_v      = v.getNormal(neighbour_boundary);
+            ev      = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
 
-            C_v = v.C;
-            Ji_v = v.Ji;
-            s_v = spdiag(v.(['s_' neighbour_boundary]));
-            m_tot_v = v.grid.N();
-
-            % Operators that are only required for own domain
-            Hi      = u.E{1}*inv(u.H)*u.E{1}' + u.E{2}*inv(u.H)*u.E{2}';
-            RHOi    = u.E{1}*inv(u.RHO)*u.E{1}' + u.E{2}*inv(u.RHO)*u.E{2}';
+            H       = u.H;
+            RHO     = u.RHO;
+            ETA     = u.ETA;
+            C       = u.C;
+            Eu      = u.Eu;
+            eU      = u.eU;
+            eGamma  = u.eGamma;
+            Egamma  = u.Egamma;
 
-            % Shared operators
-            H_gamma         = u.getBoundaryQuadratureForScalarField(boundary);
-            dim             = u.dim;
+            CV       = v.C;
+            Ev      = v.Eu;
+            eV      = v.eU;
+            eGammaV = v.eGamma;
+            nV      = v.getNormal(neighbour_boundary);
 
-            % Preallocate
-            [~, m_int] = size(H_gamma);
-            closure = sparse(dim*m_tot_u, dim*m_tot_u);
-            penalty = sparse(dim*m_tot_u, dim*m_tot_v);
-
-            % Continuity of normal displacement, term 1: The symmetric term
-            Z_u = sparse(m_int, m_int);
-            Z_v = sparse(m_int, m_int);
+            % Reduce stiffness tensors to boundary size
             for i = 1:dim
                 for j = 1:dim
-                    for l = 1:dim
-                        for k = 1:dim
-                            Z_u = Z_u + n_u{i}*n_u{j}*e_u'*Ji_u*C_u{i,j,k,l}*e_u*n_u{k}*n_u{l};
-                            Z_v = Z_v + n_v{i}*n_v{j}*e_v'*Ji_v*C_v{i,j,k,l}*e_v*n_v{k}*n_v{l};
+                    for k = 1:dim
+                        for l = 1:dim
+                            C{i,j,k,l} = e'*C{i,j,k,l}*e;
+                            CV{i,j,k,l} = ev'*CV{i,j,k,l}*ev;
                         end
                     end
                 end
             end
 
-            Z = -tuning*dim*( 1/(4*h11_u)*s_u*Z_u + 1/(4*h11_v)*s_v*Z_v );
-            closure = closure + en_u*H_gamma*Z*en_u';
-            penalty = penalty + en_u*H_gamma*Z*en_v';
-
-            % Continuity of normal displacement, term 2: The symmetrizing term
-            closure = closure + 1/2*tau_n_u*H_gamma*en_u';
-            penalty = penalty + 1/2*tau_n_u*H_gamma*en_v';
-
-            % 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';
-
-            % Multiply all normal component terms by inverse of density x quadraure
-            closure = RHOi*Hi*closure;
-            penalty = RHOi*Hi*penalty;
+            % Get elastic closure and penalty
+            [closure, penalty] = obj.elasticObj.interface(boundary, v.elasticObj, neighbour_boundary, type);
+            closure = Eu*closure*Eu';
+            penalty = Eu*penalty*Ev';
 
             % ---- Tangential tractions are imposed just like traction BC ------
-            closure = closure + obj.boundary_condition(boundary, {'t', 't'});
-
-        end
+            % We only need the viscous part
+            closure_tangential = obj.boundary_condition(boundary, {'t', 't'});
+            closure = closure + closure_tangential*Egamma*Egamma';
 
 
-        % Returns h11 for the boundary specified by the string boundary.
-        % op -- string
-        function h11 = getBorrowing(obj, boundary)
-            assertIsMember(boundary, {'w', 'e', 's', 'n'})
+            % ------ Coupling of normal component -----------
+            % Add viscous part of traction coupling
+            for i = 1:dim
+                for j = 1:dim
+                    for k = 1:dim
+                        for l = 1:dim
+                            for m = 1:dim
+                                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}') );
+                                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}') );
+                            end
+                        end
+                    end
+                end
+            end
 
-            switch boundary
-            case {'w','e'}
-                h11 = obj.refObj.h11{1};
-            case {'s', 'n'}
-                h11 = obj.refObj.h11{2};
+            % Add penalty to strain rate eq
+            for i = 1:dim
+                for j = 1:dim
+                    for k = 1:dim
+                        for l = 1:dim
+                            for m = 1:dim
+                                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}') );
+                            end
+                        end
+                    end
+                end
             end
+            %-------------------------------------------------
+
         end
 
         % Returns the outward unit normal vector for the boundary specified by the string boundary.
@@ -523,7 +514,7 @@
             assertIsMember(boundary, {'w', 'e', 's', 'n'})
             assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en', 'et', 'tau_n', 'tau_t'})
 
-            o = obj.([op, '_', boundary]);
+            o = obj.elasticObj.([op, '_', boundary]);
 
         end