--- a/+scheme/LaplaceCurvilinear.m	Fri Jun 21 17:28:33 2019 +0200
+++ b/+scheme/LaplaceCurvilinear.m	Sun Sep 22 19:05:17 2019 -0700
@@ -9,6 +9,7 @@
         order % Order of accuracy for the approximation
 
         a,b % Parameters of the operator
+        weight % Parameter in front of time derivative (e.g. u_tt in wave equation) here: 1/a.
 
 
         % Inner products and operators for physical coordinates
@@ -62,6 +63,11 @@
             if isa(a, 'function_handle')
                 a = grid.evalOn(g, a);
             end
+
+            % If a is scalar
+            if length(a) == 1
+                a = a*ones(g.N(), 1);
+            end
             a = spdiag(a);
 
             if isa(b, 'function_handle')
@@ -239,6 +245,7 @@
             obj.dim = dim;
 
             obj.a = a;
+            obj.weight = inv(a);
             obj.b = b;
             obj.a11 = a11;
             obj.a12 = a12;
@@ -327,19 +334,25 @@
         %          -- interpolation:    type of interpolation, default 'none'
         function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
 
-            % error('Not implemented')
-
+            defaultType.coupling = 'sat';
             defaultType.tuning = 1.0;
             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);
+            switch type.coupling
+            case {'cg', 'CG'}
+                [closure, penalty] = interfaceCG(obj,boundary,neighbour_scheme,neighbour_boundary,type);
+            case {'sat', 'SAT'}
+                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
             otherwise
-                error('Unknown type of interpolation: %s ', type.interpolation);
+                error('Unknown type of coupling: %s ', type.coupling);
             end
         end
 
@@ -416,6 +429,29 @@
                           +a*Hi*e_u*tau*H_b*e_v';
         end
 
+        function [closure, penalty] = interfaceCG(obj,boundary,neighbour_scheme,neighbour_boundary,type)
+
+            % There is no penalty, only a closure. And the closure is the same as for Neumann BC
+            e               = obj.getBoundaryOperator('e', boundary);
+            d               = obj.getBoundaryOperator('d', boundary);
+            H_b             = obj.getBoundaryQuadrature(boundary);
+
+            e_v             = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary);
+
+            Hi = obj.Hi;
+            a = obj.a;
+            b_b = e'*obj.b*e;
+
+            tau1 = -1;
+            tau2 = 0;
+            tau = (tau1*e + tau2*d)*H_b;
+
+            closure =  a*Hi*tau*b_b*d';
+
+            % Zero penalty of correct dimensions
+            penalty = 0*e*e_v';
+        end
+
         function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
 
             % TODO: Make this work for curvilinear grids