--- a/+scheme/Elastic2dVariableAnisotropicIncompatible.m	Tue Apr 30 14:18:47 2024 +0200
+++ b/+scheme/Elastic2dVariableAnisotropicIncompatible.m	Tue Apr 30 14:58:35 2024 +0200
@@ -532,7 +532,7 @@
         %          -- interpolation:    type of interpolation, default 'none'
         function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
 
-            defaultType.tuning = 1.0;
+            defaultType.tuning = 1.2;
             defaultType.interpolation = 'none';
             default_struct('type', defaultType);
 
@@ -558,22 +558,33 @@
             % Operators, u side
             e_u       = u.getBoundaryOperatorForScalarField('e', boundary);
             tau_u     = u.getBoundaryOperator('tau', boundary);
-            h11_u     = u.getBorrowing(boundary);
+            [h11_u, th_R_u] = u.getBorrowing(boundary);
             nu_u      = u.getNormal(boundary);
+            beta_u    = u.computeBorrowingFromC(boundary);
 
             E_u = u.E;
             C_u = u.C;
             m_tot_u = u.grid.N();
+            beta_inv_u = spdiag(1./beta_u);
 
             % Operators, v side
             e_v       = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
             tau_v     = v.getBoundaryOperator('tau', neighbour_boundary);
-            h11_v     = v.getBorrowing(neighbour_boundary);
+            [h11_v, th_R_v] = v.getBorrowing(neighbour_boundary);
             nu_v      = v.getNormal(neighbour_boundary);
+            beta_v    = v.computeBorrowingFromC(neighbour_boundary);
 
             E_v = v.E;
             C_v = v.C;
             m_tot_v = v.grid.N();
+            beta_inv_v = spdiag(1./beta_v);
+            
+            % Compensate at corners
+            [~, N] = size(e_u);
+            corner_weight = ones(N, 1);
+            corner_weight(1) = 2;
+            corner_weight(end) = 2;
+            corner_weight =  spdiag(corner_weight);
 
             % Fix {'e', 's'}, {'w', 'n'}, and {'x','x'} couplings
             flipFlag = false;
@@ -636,12 +647,12 @@
                             Z_v = Z_v + e_v'*nu_v(l)*C_v{l,i,k,j}*nu_v(k)*e_v;
                         end
                     end
-
+                    Z_u = (corner_weight/(4*h11_u) + beta_inv_u/(4*th_R_u))*Z_u; 
+                    Z_v = (corner_weight/(4*h11_v) + beta_inv_v/(4*th_R_v))*Z_v;
                     if flipFlag
                         Z_v = rot90(Z_v,2);
                     end
-
-                    Z = -tuning*dim*( 1/(4*h11_u)*Z_u + 1/(4*h11_v)*Z_v );
+                    Z = -tuning*(Z_u + 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_flip'*E_v{j}';