--- a/+scheme/LaplaceCurvilinearVirtaMin.m	Mon Jun 10 10:43:12 2019 +0200
+++ b/+scheme/LaplaceCurvilinearVirtaMin.m	Mon Jun 10 13:27:29 2019 +0200
@@ -37,8 +37,12 @@
         du_s, dv_s
         du_n, dv_n
         gamm_u, gamm_v
+        h11_u, h11_v
         lambda
 
+        % Number of boundary points in minumum check
+        bp
+
     end
 
     methods
@@ -54,6 +58,16 @@
                 error('Not implemented yet')
             end
 
+            % Number of boundary points in minimum check
+            switch order
+            case 2
+                obj.bp = 2;
+            case 4
+                obj.bp = 4;
+            case 6
+                obj.bp = 7;
+            end
+
             % assert(isa(g, 'grid.Curvilinear'))
             if isa(a, 'function_handle')
                 a = grid.evalOn(g, a);
@@ -224,6 +238,9 @@
 
             obj.gamm_u = h_u*ops_u.borrowing.M.d1;
             obj.gamm_v = h_v*ops_v.borrowing.M.d1;
+
+            obj.h11_u = h_u*ops_u.borrowing.H11;
+            obj.h11_v = h_v*ops_v.borrowing.H11;
         end
 
 
@@ -240,15 +257,43 @@
 
             [e, d] = obj.getBoundaryOperator({'e', 'd'}, boundary);
             H_b = obj.getBoundaryQuadrature(boundary);
-            gamm = obj.getBoundaryBorrowing(boundary);
+            [gamm, h11] = obj.getBoundaryBorrowing(boundary);
+
+            lambda = obj.lambda;
+            mx = obj.m(1);
+            my = obj.m(2);
+            bp = obj.bp;
+
+            lambda_mat = reshape(lambda, my, mx);
+            switch boundary
+            case 'w'
+                lambda_min = min(lambda_mat(:,1:bp), [], 2);
+                lambda_mat = repmat(lambda_min, 1, mx);
+            case 'e'
+                lambda_min = min(lambda_mat(:,end-bp+1:end), [], 2);
+                lambda_mat = repmat(lambda_min, 1, mx);
+            case 's'
+                lambda_min = min(lambda_mat(1:bp,:), [], 1);
+                lambda_mat = repmat(lambda_min, my, 1);
+            case 'n'
+                lambda_min = min(lambda_mat(end-bp+1:end,:), [], 1);
+                lambda_mat = repmat(lambda_min, my, 1);
+            end
+            lambda_min = lambda_mat(:);
 
             switch type
                 % Dirichlet boundary condition
                 case {'D','d','dirichlet'}
                     tuning = 1.0;
 
-                    b1 = gamm*obj.lambda./obj.a11.^2;
-                    b2 = gamm*obj.lambda./obj.a22.^2;
+                    switch boundary
+                    case {'w', 'e'}
+                        b1 = gamm*lambda_min./obj.a11.^2;
+                        b2 = h11 *lambda./obj.a22.^2;
+                    case {'s', 'n'}
+                        b1 = h11 *lambda./obj.a11.^2;
+                        b2 = gamm*lambda_min./obj.a22.^2;
+                    end
 
                     tau1 = tuning * spdiag(-1./b1 - 1./b2);
                     tau2 = 1;
@@ -281,7 +326,7 @@
         %          -- 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);
 
@@ -303,20 +348,79 @@
             [e_u, d_u] = obj.getBoundaryOperator({'e', 'd'}, boundary);
             H_b_u = obj.getBoundaryQuadrature(boundary);
             I_u = obj.getBoundaryIndices(boundary);
-            gamm_u = obj.getBoundaryBorrowing(boundary);
+            [gamm_u, h11_u] = obj.getBoundaryBorrowing(boundary);
 
             [e_v, d_v] = neighbour_scheme.getBoundaryOperator({'e', 'd'}, neighbour_boundary);
             H_b_v = neighbour_scheme.getBoundaryQuadrature(neighbour_boundary);
             I_v = neighbour_scheme.getBoundaryIndices(neighbour_boundary);
-            gamm_v = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary);
+            [gamm_v, h11_v] = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary);
 
             u = obj;
             v = neighbour_scheme;
 
-            b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2;
-            b2_u = gamm_u*u.lambda(I_u)./u.a22(I_u).^2;
-            b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2;
-            b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2;
+            % Minimum check for u
+            lambda_u = obj.lambda;
+            mx = obj.m(1);
+            my = obj.m(2);
+            bp = obj.bp;
+
+            lambda_mat = reshape(lambda_u, my, mx);
+            switch boundary
+            case 'w'
+                lambda_min = min(lambda_mat(:,1:bp), [], 2);
+                lambda_mat = repmat(lambda_min, 1, mx);
+            case 'e'
+                lambda_min = min(lambda_mat(:,end-bp+1:end), [], 2);
+                lambda_mat = repmat(lambda_min, 1, mx);
+            case 's'
+                lambda_min = min(lambda_mat(1:bp,:), [], 1);
+                lambda_mat = repmat(lambda_min, my, 1);
+            case 'n'
+                lambda_min = min(lambda_mat(end-bp+1:end,:), [], 1);
+                lambda_mat = repmat(lambda_min, my, 1);
+            end
+            lambda_min_u = lambda_mat(:);
+
+            % Minimum check for v
+            lambda_v = v.lambda;
+            mx = v.m(1);
+            my = v.m(2);
+            bp = v.bp;
+
+            lambda_mat = reshape(lambda_v, my, mx);
+            switch boundary
+            case 'w'
+                lambda_min = min(lambda_mat(:,1:bp), [], 2);
+                lambda_mat = repmat(lambda_min, 1, mx);
+            case 'e'
+                lambda_min = min(lambda_mat(:,end-bp+1:end), [], 2);
+                lambda_mat = repmat(lambda_min, 1, mx);
+            case 's'
+                lambda_min = min(lambda_mat(1:bp,:), [], 1);
+                lambda_mat = repmat(lambda_min, my, 1);
+            case 'n'
+                lambda_min = min(lambda_mat(end-bp+1:end,:), [], 1);
+                lambda_mat = repmat(lambda_min, my, 1);
+            end
+            lambda_min_v = lambda_mat(:);
+
+            switch boundary
+            case {'w', 'e'}
+                b1_u = gamm_u*u.lambda_min_u(I_u)./u.a11(I_u).^2;
+                b2_u = h11_u *u.lambda(I_u)./u.a22(I_u).^2;
+            case {'s', 'n'}
+                b1_u = h11_u *u.lambda(I_u)./u.a11(I_u).^2;
+                b2_u = gamm_u*u.lambda_min(I_u)./u.a22(I_u).^2;
+            end
+
+            switch neighbour_boundary
+            case {'w', 'e'}
+                b1_v = gamm_v*v.lambda_min_v(I_v)./v.a11(I_v).^2;
+                b2_v = h11_v *v.lambda(I_v)./v.a22(I_v).^2;
+            case {'s', 'n'}
+                b1_v = h11_v *v.lambda(I_v)./v.a11(I_v).^2;
+                b2_v = gamm_v*v.lambda_min(I_v)./v.a22(I_v).^2;
+            end
 
             tau1 = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v);
             tau1 = tuning * spdiag(tau1);
@@ -479,12 +583,14 @@
 
         % Returns borrowing constant gamma
         % boundary -- string
-        function gamm = getBoundaryBorrowing(obj, boundary)
+        function [gamm, h11] = getBoundaryBorrowing(obj, boundary)
             switch boundary
                 case {'w','e'}
                     gamm = obj.gamm_u;
+                    h11 = obj.h11_u;
                 case {'s','n'}
                     gamm = obj.gamm_v;
+                    h11 = obj.h11_v;
                 otherwise
                     error('No such boundary: boundary = %s',boundary);
             end