diff -r 23d9ca6755be -r 104f0af001e0 +scheme/Elastic2dVariable.m
--- a/+scheme/Elastic2dVariable.m	Tue Dec 25 07:23:38 2018 +0100
+++ b/+scheme/Elastic2dVariable.m	Tue Dec 25 08:33:35 2018 +0100
@@ -438,76 +438,28 @@
             % v denotes the solution in the neighbour domain
             % Operators without subscripts are from the own domain.
 
-            % j is the coordinate direction of the boundary
-            j = obj.get_boundary_number(boundary);
-            j_v = neighbour_scheme.get_boundary_number(neighbour_boundary);
-
             % Get boundary operators
-            [e, T, tau, H_gamma] = obj.getBoundaryOperator({'e','T','tau','H'}, boundary);
-            [e_v, tau_v] = neighbour_scheme.getBoundaryOperator({'e','tau'}, neighbour_boundary);
+            [e, tau] = obj.getBoundaryOperator({'e_tot','tau_tot'}, boundary);
+            [e_v, tau_v] = neighbour_scheme.getBoundaryOperator({'e_tot','tau_tot'}, neighbour_boundary);
+            H_gamma = obj.getBoundaryQuadrature(boundary);
 
             % Operators and quantities that correspond to the own domain only
-            Hi = obj.Hi;
-            RHOi = obj.RHOi;
-            dim = obj.dim;
-
-            %--- Other operators ----
-            m_tot_u = obj.grid.N();
-            E = obj.E;
-            LAMBDA_u = obj.LAMBDA;
-            MU_u = obj.MU;
-            lambda_u = e'*LAMBDA_u*e;
-            mu_u = e'*MU_u*e;
+            Hi = obj.Hi_kron;
+            RHOi = obj.RHOi_kron;
 
-            m_tot_v = neighbour_scheme.grid.N();
-            E_v = neighbour_scheme.E;
-            LAMBDA_v = neighbour_scheme.LAMBDA;
-            MU_v = neighbour_scheme.MU;
-            lambda_v = e_v'*LAMBDA_v*e_v;
-            mu_v = e_v'*MU_v*e_v;
-            %-------------------------
-
-            % Borrowing constants
-            theta_R_u = obj.theta_R{j};
-            theta_H_u = obj.theta_H{j};
-            theta_M_u = obj.theta_M{j};
-
-            theta_R_v = neighbour_scheme.theta_R{j_v};
-            theta_H_v = neighbour_scheme.theta_H{j_v};
-            theta_M_v = neighbour_scheme.theta_M{j_v};
+            % Penalty strength operators
+            alpha_u = 1/4*obj.getBoundaryOperator('alpha_tot', boundary);
+            alpha_v = 1/4*neighbour_scheme.getBoundaryOperator('alpha_tot', neighbour_boundary);
 
-            function [alpha_ii, alpha_ij] = computeAlpha(th_R, th_H, th_M, lambda, mu)
-                alpha_ii = dim*lambda/(4*th_H) + lambda/(4*th_R) + mu/(2*th_M);
-                alpha_ij = mu/(2*th_H) + mu/(4*th_R);
-            end
-
-            [alpha_ii_u, alpha_ij_u] = computeAlpha(theta_R_u, theta_H_u, theta_M_u, lambda_u, mu_u);
-            [alpha_ii_v, alpha_ij_v] = computeAlpha(theta_R_v, theta_H_v, theta_M_v, lambda_v, mu_v);
-            sigma_ii = alpha_ii_u + alpha_ii_v;
-            sigma_ij = alpha_ij_u + alpha_ij_v;
-
-            d = @kroneckerDelta;  % Kronecker delta
-            db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
-            sigma = @(i,j) tuning*(d(i,j)*sigma_ii + db(i,j)*sigma_ij);
+            closure = -RHOi*Hi*e*H_gamma*(alpha_u' + alpha_v'*e_v*e');
+            penalty = RHOi*Hi*e*H_gamma*(alpha_u'*e*e_v' + alpha_v');
 
-            % Preallocate
-            closure = sparse(dim*m_tot_u, dim*m_tot_u);
-            penalty = sparse(dim*m_tot_u, dim*m_tot_v);
-
-            % Loop over components that penalties end up on
-            for i = 1:dim
-                closure = closure - E{i}*RHOi*Hi*e*sigma(i,j)*H_gamma*e'*E{i}';
-                penalty = penalty + E{i}*RHOi*Hi*e*sigma(i,j)*H_gamma*e_v'*E_v{i}';
+            closure = closure - 1/2*RHOi*Hi*e*H_gamma*tau';
+            penalty = penalty - 1/2*RHOi*Hi*e*H_gamma*tau_v';
 
-                closure = closure - 1/2*E{i}*RHOi*Hi*e*H_gamma*tau{i}';
-                penalty = penalty - 1/2*E{i}*RHOi*Hi*e*H_gamma*tau_v{i}';
+            closure = closure + 1/2*RHOi*Hi*tau*H_gamma*e';
+            penalty = penalty - 1/2*RHOi*Hi*tau*H_gamma*e_v';
 
-                % Loop over components that we have interface conditions on
-                for k = 1:dim
-                    closure = closure + 1/2*E{i}*RHOi*Hi*T{k,i}*H_gamma*e'*E{k}';
-                    penalty = penalty - 1/2*E{i}*RHOi*Hi*T{k,i}*H_gamma*e_v'*E_v{k}';
-                end
-            end
         end
 
         function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)