--- a/+scheme/Elastic2dCurvilinearAnisotropicUpwind.m	Sun May 31 20:16:54 2020 -0700
+++ b/+scheme/Elastic2dCurvilinearAnisotropicUpwind.m	Sun May 31 20:18:45 2020 -0700
@@ -64,7 +64,7 @@
         % optFlag -- if true, extra computations are performed, which may be helpful for optimization.
         function obj = Elastic2dCurvilinearAnisotropicUpwind(g, order, rho, C, opSet, optFlag)
             default_arg('rho', @(x,y) 0*x+1);
-            default_arg('opSet',{@sbp.D2VariableCompatible, @sbp.D2VariableCompatible});
+            default_arg('opSet',{@sbp.D1Upwind, @sbp.D1Upwind});
             default_arg('optFlag', false);
             dim = 2;
 
@@ -105,10 +105,11 @@
             m_tot = g.N();
 
             % 1D operators
+            opSetMetric = {@sbp.D2VariableCompatible, @sbp.D2VariableCompatible};
             m_u = m(1);
             m_v = m(2);
-            ops_u = opSet{1}(m_u, {0, 1}, order);
-            ops_v = opSet{2}(m_v, {0, 1}, order);
+            ops_u = opSetMetric{1}(m_u, {0, 1}, order);
+            ops_v = opSetMetric{2}(m_v, {0, 1}, order);
 
             h_u = ops_u.h;
             h_v = ops_v.h;
@@ -184,7 +185,7 @@
             end
 
             gRef = grid.equidistant([m_u, m_v], {0,1}, {0,1});
-            refObj = scheme.Elastic2dVariableAnisotropicUpwind(gRef, order, rho_tilde, PHI, []);
+            refObj = scheme.Elastic2dVariableAnisotropicUpwind(gRef, order, rho_tilde, PHI, opSet);
 
             %---- Set object properties ------
             obj.RHO = spdiag(rho);
@@ -192,7 +193,7 @@
             % Volume quadrature
             obj.J = spdiag(J);
             obj.Ji = spdiag(1./J);
-            obj.H = obj.J*kr(H_u,H_v);
+            obj.H = obj.J*refObj.H;
 
             % Boundary quadratures
             s_w = sqrt((e_w'*x_v).^2 + (e_w'*y_v).^2);