--- a/+scheme/Wave2dCurve.m	Fri Nov 03 10:43:27 2017 +0100
+++ b/+scheme/Wave2dCurve.m	Fri Nov 03 10:53:15 2017 +0100
@@ -27,6 +27,8 @@
         gamm_u, gamm_v
         lambda
 
+        Dx, Dy % Physical derivatives
+
         x_u
         x_v
         y_u
@@ -38,6 +40,8 @@
             default_arg('opSet',@sbp.D2Variable);
             default_arg('c', 1);
 
+            warning('Use LaplaceCruveilinear instead')
+
             assert(isa(g, 'grid.Curvilinear'))
 
             m = g.size();
@@ -153,6 +157,9 @@
             obj.D = obj.Ji*c^2*(Duu + Duv + Dvu + Dvv);
             obj.lambda = lambda;
 
+            obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv;
+            obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv;
+
             obj.gamm_u = h_u*ops_u.borrowing.M.d1;
             obj.gamm_v = h_v*ops_v.borrowing.M.d1;
         end
@@ -225,6 +232,7 @@
 
                     tau = -c.^2 * 1/beta*obj.Ji*e;
 
+                    warning('is this right?! /c?')
                     closure{1} = halfnorm_inv*tau/c*spdiag(scale_factor)*e';
                     closure{2} = halfnorm_inv*tau*beta*d';
                     penalty = -halfnorm_inv*tau;