--- a/+scheme/Elastic2dVariableAnisotropic.m	Thu Sep 05 17:10:58 2019 -0700
+++ b/+scheme/Elastic2dVariableAnisotropic.m	Fri Sep 06 14:21:18 2019 -0700
@@ -407,7 +407,6 @@
             switch type
 
             % Dirichlet boundary condition
-            % OBS! Cannot yet set one component at a time unless one assumes Displacement for all components
             case {'D','d','dirichlet','Dirichlet','displacement','Displacement'}
 
                 if numel(bc) >= 3
@@ -458,11 +457,11 @@
 
         % type     Struct that specifies the interface coupling.
         %          Fields:
-        %          -- tuning:           penalty strength, defaults to 1.2
+        %          -- tuning:           penalty strength, defaults to 1.0
         %          -- interpolation:    type of interpolation, default 'none'
         function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
 
-            defaultType.tuning = 1.2;
+            defaultType.tuning = 1.0;
             defaultType.interpolation = 'none';
             default_struct('type', defaultType);
 
@@ -481,33 +480,80 @@
 
             % u denotes the solution in the own domain
             % v denotes the solution in the neighbour domain
-            % Operators without subscripts are from the own domain.
+
+            u = obj;
+            v = neighbour_scheme;
+
+            % Operators, u side
+            e_u       = u.getBoundaryOperatorForScalarField('e', boundary);
+            tau_u     = u.getBoundaryOperator('tau', boundary);
+            h11_u     = u.getBorrowing(boundary);
+            nu_u      = u.getNormal(boundary);
+
+            E_u = u.E;
+            C_u = u.C;
+            m_tot_u = u.grid.N();
 
-            % Get boundary operators
-            e   = obj.getBoundaryOperator('e', boundary);
-            tau = obj.getBoundaryOperator('tau', boundary);
+            % Operators, v side
+            e_v       = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
+            tau_v     = v.getBoundaryOperator('tau', neighbour_boundary);
+            h11_v     = v.getBorrowing(neighbour_boundary);
+            nu_v      = v.getNormal(neighbour_boundary);
+
+            E_v = v.E;
+            C_v = v.C;
+            m_tot_v = v.grid.N();
 
-            e_v   = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary);
-            tau_v = neighbour_scheme.getBoundaryOperator('tau', neighbour_boundary);
+            % Operators that are only required for own domain
+            Hi      = u.Hi_kron;
+            RHOi    = u.RHOi_kron;
+            e_kron  = u.getBoundaryOperator('e', boundary);
+            T_u     = u.getBoundaryTractionOperator(boundary);
 
-            H_gamma = obj.getBoundaryQuadrature(boundary);
+            % Shared operators
+            H_gamma         = u.getBoundaryQuadratureForScalarField(boundary);
+            H_gamma_kron    = u.getBoundaryQuadrature(boundary);
+            dim             = u.dim;
 
-            % Operators and quantities that correspond to the own domain only
-            Hi = obj.Hi_kron;
-            RHOi = obj.RHOi_kron;
+            % 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 displacement ------
 
-            % Penalty strength operators
-            alpha_u = 1/4*tuning*obj.getBoundaryOperator('alpha', boundary);
-            alpha_v = 1/4*tuning*neighbour_scheme.getBoundaryOperator('alpha', neighbour_boundary);
+            % Y: symmetrizing part of penalty
+            % Z: symmetric part of penalty
+            % X = Y + Z.
+
+            % Loop over components to couple across interface
+            for j = 1:dim
 
-            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');
+                % Loop over components that penalties end up on
+                for i = 1:dim
+                    Y = 1/2*T_u{j,i}';
+                    Z_u = sparse(m_int, m_int);
+                    Z_v = sparse(m_int, m_int);
+                    for l = 1:dim
+                        for k = 1:dim
+                            Z_u = Z_u + e_u'*nu_u(l)*C_u{i,l,k,j}*nu_u(k)*e_u;
+                            Z_v = Z_v + e_v'*nu_v(l)*C_v{i,l,k,j}*nu_v(k)*e_v;
+                        end
+                    end
+                    Z = -tuning*dim*( 1/(4*h11_u)*Z_u + 1/(4*h11_v)*Z_v );
+                    X = Y + Z*e_u';
+                    closure = closure + E_u{i}*X'*H_gamma*e_u'*E_u{j}';
+                    penalty = penalty - E_u{i}*X'*H_gamma*e_v'*E_v{j}';
+                end
+            end
 
-            closure = closure - 1/2*RHOi*Hi*e*H_gamma*tau';
-            penalty = penalty - 1/2*RHOi*Hi*e*H_gamma*tau_v';
+            % ---- Continuity of traction ------
+            closure = closure - 1/2*e_kron*H_gamma_kron*tau_u';
+            penalty = penalty - 1/2*e_kron*H_gamma_kron*tau_v';
 
-            closure = closure + 1/2*RHOi*Hi*tau*H_gamma*e';
-            penalty = penalty - 1/2*RHOi*Hi*tau*H_gamma*e_v';
+            % ---- Multiply by inverse of density x quadraure ----
+            closure = RHOi*Hi*closure;
+            penalty = RHOi*Hi*penalty;
 
         end