--- a/+rv/constructDiffOps.m	Fri Dec 14 18:33:10 2018 +0100
+++ b/+rv/constructDiffOps.m	Wed Dec 19 20:00:27 2018 +0100
@@ -1,13 +1,14 @@
-function [D_rv, D_flux, DvDt, solutionPenalties, residualPenalties] = constructDiffOps(scheme, grid, order, opSet, waveSpeed, BCs)
+function [D_rv, D_flux, DvDt, solutionPenalties, residualPenalties] = constructDiffOps(scheme, grid, order, opSet, waveSpeed, BCs, fluxSplitting)
+    default_arg('fluxSplitting',[]);
     assert(size(BCs,1) == grid.D());
-    
+
     %% DiffOps for solution vector
-    [D, solutionPenalties] = constructTotalFluxDiffOp(scheme, grid, order, opSet, waveSpeed, BCs);
+    [D, solutionPenalties] = constructTotalFluxDiffOp(scheme, grid, order, opSet, waveSpeed, BCs, fluxSplitting);
     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);
+    [D_flux, residualPenalties] = constructTotalFluxDiffOp(scheme, grid, order-2, opSet, waveSpeed, BCs, fluxSplitting);
     % DiffOp for flux in residual viscosity. Due to sign conventions of the implemnted schemes, we need to
     % change the sign.
     D_flux = -D_flux;
@@ -16,12 +17,14 @@
     DvDt = @(v)D*v;
 end
 
-function [D, penalties] = constructTotalFluxDiffOp(scheme, grid, order, opSet, waveSpeed, BCs)
-    diffOp = scheme(grid, order, opSet, waveSpeed);
+function [D, penalties] = constructTotalFluxDiffOp(scheme, grid, order, opSet, waveSpeed, BCs, fluxSplitting)
+    if isequal(opSet, @sbp.D1Upwind)
+        assert(size(fluxSplitting,1) == grid.D());
+        diffOp = scheme(grid, order, opSet, waveSpeed, fluxSplitting);
+    else
+        diffOp = scheme(grid, order, opSet, waveSpeed);
+    end
     [D, penalties]  = addClosuresToDiffOp(diffOp, BCs);
-    if (isequal(opSet,@sbp.D1Upwind))
-        D = D + getUpwindDissOp(diffOp, waveSpeed);
-    end
 end
 
 function [D, penalties] = addClosuresToDiffOp(diffOp, BCs)
@@ -35,24 +38,11 @@
     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
+    % or if it is decided that it should only be supported for upwind operators.
     assert(isequal(opSet,@sbp.D1Upwind))
 
     m = grid.size();