--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+rv/+diffops/constructSymmetricD2.m	Thu Jun 27 11:04:34 2019 +0200
@@ -0,0 +1,54 @@
+function D2 = constructSymmetricD2(g, order, opSet)
+    m = g.size();
+    ops = cell(g.D(),1);
+    I = cell(g.D(),1);
+    for i = 1:g.D()
+       lim = {g.x{i}(1), g.x{i}(end)};
+       ops{i} = opSet(m(i), lim, order);
+       I{i} = speye(m(i));
+    end
+
+    switch g.D()
+        case 1
+            e_r = ops{1}.e_r;
+            e_l = ops{1}.e_l;
+            Hi = ops{1}.HI;
+            B = e_r*e_r' - e_l*e_l';
+            if isequal(opSet,@sbp.D1Upwind)
+                Dm = ops{1}.Dm;
+                Dp = ops{1}.Dp;
+                M =  Dm - Hi*B;
+                D2 = @(Viscosity) M*Viscosity*Dp;
+            else
+                % TODO: Fix closure for D2Variable
+                % TODO: Fix Viscosity not being vector
+                D2 = @(Viscosity)ops{1}.D2(diag(Viscosity));
+            end
+        case 2
+            % TODO: 
+            % Currently only implemented for upwind operators.
+            % Remove this part once the time-dependent D2 operator is implemented for other opSets
+            % or if it is decided that it should only be supported for upwind operators.
+            assert(isequal(opSet,@sbp.D1Upwind))
+            e_e = kron(ops{1}.e_r,I{2});
+            e_w = kron(ops{1}.e_l,I{2});
+            Dm_x = kron(ops{1}.Dm,I{2});
+            Dp_x  = kron(ops{1}.Dp,I{2});
+            H_x  = kron(ops{1}.HI,I{2});
+            B_x = e_e*e_e' - e_w*e_w';
+            M_x = Dm_x-H_x*B_x;
+            D2_x = @(Viscosity) M_x*Viscosity*Dp_x;
+
+            e_n = kron(I{1},ops{2}.e_r);
+            e_s = kron(I{1},ops{2}.e_l);
+            Dm_y = kron(I{1},ops{2}.Dm);
+            Dp_y  = kron(I{1},ops{2}.Dp);
+            H_y = kron(I{1},ops{2}.HI);
+            B_y = e_n*e_n' - e_s*e_s';
+            M_y = Dm_y-H_y*B_y;
+            D2_y = @(Viscosity) M_y*Viscosity*Dp_y;
+            D2 = @(Viscosity) D2_x(Viscosity) + D2_y(Viscosity);
+        otherwise
+            error('3D not yet implemented')
+    end
+end
\ No newline at end of file