--- a/+scheme/Elastic2dVariable.m	Wed Dec 19 06:54:47 2018 +0100
+++ b/+scheme/Elastic2dVariable.m	Wed Dec 19 06:58:10 2018 +0100
@@ -32,10 +32,14 @@
 
         H, Hi, H_1D % Inner products
         e_l, e_r
+        e_w, e_e, e_s, e_n
+
+
         d1_l, d1_r % Normal derivatives at the boundary
         E % E{i}^T picks out component i
 
         H_boundary % Boundary inner products
+        H_w, H_e, H_s, H_n
 
         % Kroneckered norms and coefficients
         RHOi_kron
@@ -146,6 +150,10 @@
             obj.e_l{2} = kron(I{1},e_l{2});
             obj.e_r{1} = kron(e_r{1},I{2});
             obj.e_r{2} = kron(I{1},e_r{2});
+            obj.e_w = obj.e_l{1};
+            obj.e_e = obj.e_r{1};
+            obj.e_s = obj.e_l{2};
+            obj.e_n = obj.e_r{2};
 
             obj.d1_l{1} = kron(d1_l{1},I{2});
             obj.d1_l{2} = kron(I{1},d1_l{2});
@@ -184,6 +192,10 @@
             obj.H_boundary{1} = H{2};
             obj.H_boundary{2} = H{1};
             obj.H_1D = {H{1}, H{2}};
+            obj.H_w = H{2};
+            obj.H_e = H{2};
+            obj.H_s = H{1};
+            obj.H_n = H{1};
 
             % E{i}^T picks out component i.
             E = cell(dim,1);
@@ -330,6 +342,7 @@
             j = obj.get_boundary_number(boundary);
             [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary);
 
+
             E = obj.E;
             Hi = obj.Hi;
             LAMBDA = obj.LAMBDA;
@@ -354,7 +367,7 @@
                 % Loop over components that Dirichlet penalties end up on
                 for i = 1:dim
                     C = transpose(T{k,i});
-                    A = -alpha{i,k};
+                    A = -e*transpose(alpha{i,k});
                     B = A + e*C;
                     closure = closure + E{i}*RHOi*Hi*B'*e*H_gamma*(e'*E{k}' );
                     penalty = penalty - E{i}*RHOi*Hi*B'*e*H_gamma;
@@ -392,6 +405,8 @@
         end
 
         function [closure, penalty] = interfaceStandard(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
             % Operators without subscripts are from the own domain.
@@ -457,8 +472,8 @@
                 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*E{i}*RHOi*Hi*e*H_gamma*e'*tau{i};
-                penalty = penalty - 1/2*E{i}*RHOi*Hi*e*H_gamma*e_v'*tau_v{i};
+                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}';
 
                 % Loop over components that we have interface conditions on
                 for k = 1:dim
@@ -525,8 +540,6 @@
                             case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
                                 varargout{k} = obj.d1_r{j};
                         end
-                    case 'H'
-                        varargout{k} = obj.H_boundary{j};
                     case 'T'
                         switch boundary
                             case {'w','W','west','West','s','S','south','South'}
@@ -541,8 +554,17 @@
                             case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
                                 varargout{k} = obj.tau_r{j};
                         end
+                    case 'H'
+                        varargout{k} = obj.H_boundary{j};
                     case 'alpha'
                         % alpha = alpha(i,j) is the penalty strength for displacement BC.
+                        switch boundary
+                            case {'w','W','west','West','s','S','south','South'}
+                                e = obj.e_l{j};
+                            case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
+                                e = obj.e_r{j};
+                        end
+
                         tuning = 1.2;
                         LAMBDA = obj.LAMBDA;
                         MU = obj.MU;
@@ -565,7 +587,7 @@
                                               + db(i,j)*a_mu_ij*MU );
                         for i = 1:obj.dim
                             for l = 1:obj.dim
-                                alpha{i,l} = d(i,l)*alpha_func(i,j);
+                                alpha{i,l} = d(i,l)*alpha_func(i,j)*e;
                             end
                         end