--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+rv/constructDiffOps.m	Fri Dec 14 18:33:10 2018 +0100
@@ -0,0 +1,97 @@
+function [D_rv, D_flux, DvDt, solutionPenalties, residualPenalties] = constructDiffOps(scheme, grid, order, opSet, waveSpeed, BCs)
+    assert(size(BCs,1) == grid.D());
+    
+    %% DiffOps for solution vector
+    [D, solutionPenalties] = constructTotalFluxDiffOp(scheme, grid, order, opSet, waveSpeed, BCs);
+    D2 = constructSymmetricD2Operator(grid, order, opSet);
+    D_rv = @(v,viscosity)(D + D2(viscosity))*v;
+
+    %% DiffOps for residual viscosity
+    [D_flux, residualPenalties] = constructTotalFluxDiffOp(scheme, grid, order-2, opSet, waveSpeed, BCs);
+    % DiffOp for flux in residual viscosity. Due to sign conventions of the implemnted schemes, we need to
+    % change the sign.
+    D_flux = -D_flux;
+    D_flux = @(v) D_flux*v;
+    % DiffOp for time derivative in residual viscosity
+    DvDt = @(v)D*v;
+end
+
+function [D, penalties] = constructTotalFluxDiffOp(scheme, grid, order, opSet, waveSpeed, BCs)
+    diffOp = scheme(grid, order, opSet, waveSpeed);
+    [D, penalties]  = addClosuresToDiffOp(diffOp, BCs);
+    if (isequal(opSet,@sbp.D1Upwind))
+        D = D + getUpwindDissOp(diffOp, waveSpeed);
+    end
+end
+
+function [D, penalties] = addClosuresToDiffOp(diffOp, BCs)
+    D = diffOp.D;
+    penalties = cell(size(BCs));
+    for i = 1:size(BCs,1)
+        for j = 1:size(BCs,2)
+            [closure, penalties{i,j}] = diffOp.boundary_condition(BCs{i,j}.boundary, BCs{i,j}.type);
+            D = D + closure;
+        end
+    end
+end
+
+% TODO: Pass flux splitting components here instead of hard-coding the splitting.
+% TBD: Decide if this should be done in the scheme instead of by the user.
+function  DissOp = getUpwindDissOp(diffOp, waveSpeed)
+    switch diffOp.grid.D()
+        case 1
+            DissOp = -(abs(waveSpeed(1))*diffOp.Diss);
+        case 2
+            DissOp = -(abs(waveSpeed(1))*diffOp.DissOpx + abs(waveSpeed(2))*diffOp.DissOpy);
+        otherwise
+            error('3D not yet implemented')
+    end
+    
+end
+
+function D2 = constructSymmetricD2Operator(grid, order, opSet)
+    % TODO: 
+    % Currently only implemented for upwind operators.
+    % Remove this part once the time-dependent D2 operator is implemented for other opSets
+    assert(isequal(opSet,@sbp.D1Upwind))
+
+    m = grid.size();
+    ops = cell(grid.D(),1);
+    I = cell(grid.D(),1);
+    for i = 1:grid.D()
+       lim = {grid.x{i}(1), grid.x{i}(end)};
+       ops{i} = opSet(m(i), lim, order);
+       I{i} = speye(m(i));
+    end
+
+    % TBD: How is this generalized to a loop over dimensions or similar?
+    switch grid.D()
+        case 1
+            e_r = ops{1}.e_r;
+            e_l = ops{1}.e_l;
+            Dm = ops{1}.Dm;
+            Dp = ops{1}.Dp;
+            Hi = ops{1}.HI;
+            B = e_r*e_r' - e_l*e_l';
+            D2 = @(viscosity) Dm*spdiag(viscosity)*Dp-Hi*(B*spdiag(viscosity)*Dp);
+        case 2
+            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';
+            D2_x = @(viscosity) Dm_x*spdiag(viscosity)*Dp_x-H_x*(B_x*spdiag(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';
+            D2_y = @(viscosity) Dm_y*spdiag(viscosity)*Dp_y-H_y*(B_y*spdiag(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