--- a/+scheme/LaplaceCurvilinear.m	Fri Oct 12 08:50:25 2018 +0200
+++ b/+scheme/LaplaceCurvilinear.m	Thu Jan 17 10:44:12 2019 +0100
@@ -38,6 +38,7 @@
         du_n, dv_n
         gamm_u, gamm_v
         lambda
+
     end
 
     methods
@@ -53,7 +54,11 @@
                 error('Not implemented yet')
             end
 
-            assert(isa(g, 'grid.Curvilinear'))
+            % assert(isa(g, 'grid.Curvilinear'))
+            if isa(a, 'function_handle')
+                a = grid.evalOn(g, a);
+                a = spdiag(a);
+            end
 
             m = g.size();
             m_u = m(1);
@@ -268,11 +273,31 @@
             end
         end
 
-        function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
+        % type     Struct that specifies the interface coupling.
+        %          Fields:
+        %          -- tuning:           penalty strength, defaults to 1.2
+        %          -- interpolation:    type of interpolation, default 'none'
+        function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
+
+            defaultType.tuning = 1.2;
+            defaultType.interpolation = 'none';
+            default_struct('type', defaultType);
+
+            switch type.interpolation
+            case {'none', ''}
+                [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type);
+            case {'op','OP'}
+                [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type);
+            otherwise
+                error('Unknown type of interpolation: %s ', type.interpolation);
+            end
+        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
-            tuning = 1.2;
-            % tuning = 20.2;
             [e_u, d_u, gamm_u, H_b_u, I_u] = obj.get_boundary_ops(boundary);
             [e_v, d_v, gamm_v, H_b_v, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
 
@@ -298,15 +323,66 @@
             penalty = obj.a*obj.Hi*(-tau*e_v' + sig*d_v');
         end
 
-        % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
+        function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
+
+            % TODO: Make this work for curvilinear grids
+            warning('LaplaceCurvilinear: Non-conforming grid interpolation only works for Cartesian grids.');
+
+            % User can request special interpolation operators by specifying type.interpOpSet
+            default_field(type, 'interpOpSet', @sbp.InterpOpsOP);
+            interpOpSet = type.interpOpSet;
+            tuning = type.tuning;
+
+
+            % u denotes the solution in the own domain
+            % v denotes the solution in the neighbour domain
+            [e_u, d_u, gamm_u, H_b_u, I_u] = obj.get_boundary_ops(boundary);
+            [e_v, d_v, gamm_v, H_b_v, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
+
+            % Find the number of grid points along the interface
+            m_u = size(e_u, 2);
+            m_v = size(e_v, 2);
+
+            Hi = obj.Hi;
+            a = obj.a;
+
+            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;
+
+            tau_u = -1./(4*b1_u) -1./(4*b2_u);
+            tau_v = -1./(4*b1_v) -1./(4*b2_v);
+
+            tau_u = tuning * spdiag(tau_u);
+            tau_v = tuning * spdiag(tau_v);
+            beta_u = tau_v;
+
+            % Build interpolation operators
+            intOps = interpOpSet(m_u, m_v, obj.order, neighbour_scheme.order);
+            Iu2v = intOps.Iu2v;
+            Iv2u = intOps.Iv2u;
+
+            closure = a*Hi*e_u*tau_u*H_b_u*e_u' + ...
+                      a*Hi*e_u*H_b_u*Iv2u.bad*beta_u*Iu2v.good*e_u' + ...
+                      a*1/2*Hi*d_u*H_b_u*e_u' + ...
+                      -a*1/2*Hi*e_u*H_b_u*d_u';
+
+            penalty = -a*Hi*e_u*tau_u*H_b_u*Iv2u.good*e_v' + ...
+                      -a*Hi*e_u*H_b_u*Iv2u.bad*beta_u*e_v' + ...
+                      -a*1/2*Hi*d_u*H_b_u*Iv2u.good*e_v' + ...
+                      -a*1/2*Hi*e_u*H_b_u*Iv2u.bad*d_v';
+
+        end
+
+        % Returns the boundary ops and sign for the boundary specified by the string boundary.
         % The right boundary is considered the positive boundary
         %
-        %  I -- the indecies of the boundary points in the grid matrix
+        %  I -- the indices of the boundary points in the grid matrix
         function [e, d, gamm, H_b, I] = get_boundary_ops(obj, boundary)
-
-            % gridMatrix = zeros(obj.m(2),obj.m(1));
-            % gridMatrix(:) = 1:numel(gridMatrix);
-
             ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m));
 
             switch boundary