--- a/+rv/ResidualViscosity.m	Tue Jun 25 13:24:01 2019 +0200
+++ b/+rv/ResidualViscosity.m	Tue Jun 25 14:14:25 2019 +0200
@@ -1,70 +1,71 @@
 classdef ResidualViscosity < handle
     properties
+        grid % grid
         Df % Diff op approximating the gradient of the flux f(u)
         waveSpeed % Wave speed at each grid point, e.g f'(u). %TBD: Better naming?
         Cmax % Constant controlling relative amount of upwind dissipation
         Cres % Constant controling relative amount of upwind dissipation
         h % Length scale used for scaling the viscosity. Typically grid spacing.
         normalization % Function used to normalize the residual such that it is amplified in the
-                      % shocks.
-        boundaryIndices % indices of the boundary of the domain
+                      % shocks and suppressed elsewhere.
+        Mres % Coefficients for the residual viscosity
+        Mfirst % Coefficients for the first order viscosity
+        fRes % Function handle for computing the residual.
     end
 
     methods
         %  TBD: Decide on how to treat waveSpeed. It would be nice to just pass a constant value without
         %       wrapping it in a function.
-        function obj = ResidualViscosity(g, Df, waveSpeed, Cmax, Cres, h, normalization)
+        function obj = ResidualViscosity(grid, Df, waveSpeed, Cmax, Cres, h, normalization, postProcess)
             default_arg('normalization',@(v)abs(obj.minmaxDiffNeighborhood1d(v)-norm(v-mean(v),inf)));
+            default_arg('postProcess','maximum neighbors');
             obj.Df = Df;
             obj.waveSpeed = waveSpeed;
             obj.h = h;
             obj.Cmax = Cmax;
+            
             obj.Cres = Cres;
             obj.normalization = normalization;
-            obj.boundaryIndices = obj.getBoundaryIndices(g);
+            obj.grid = grid;
+            obj.Mres = obj.Cres*obj.h^2;
+            obj.Mfirst = obj.Cmax*obj.h;
+            switch postProcess
+                case 'none'
+                    obj.fRes = @(v,dvdt) obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v);
+                case 'maximum neighbors'
+                    obj.fRes = @obj.maxResidualNeighbors
+            end
         end
 
         function viscosity = evaluateViscosity(obj, v, dvdt)
-            viscosity = min(obj.Cmax*obj.h*abs(obj.waveSpeed(v)), obj.Cres*obj.h^2*abs(dvdt + obj.Df(v))./obj.normalization(v));
-            % TODO: If the boundary conditions are implemented correctly, this might not be needed.
-            viscosity(obj.boundaryIndices) = 0;
+            viscosity = min(obj.Mfirst*abs(obj.waveSpeed(v)), obj.fRes(v,dvdt));
         end
 
         function [viscosity, Df, firstOrderViscosity, residualViscosity] = evaluate(obj, v, dvdt)
             Df = obj.Df(v);
-            firstOrderViscosity = obj.Cmax*obj.h*abs(obj.waveSpeed(v));
-            residualViscosity = obj.Cres*obj.h^2*abs(dvdt + Df)./obj.normalization(v);
-            viscosity = min(firstOrderViscosity, residualViscosity);
-            % TODO: If the boundary conditions are implemented correctly, this might not be needed.
-            viscosity(obj.boundaryIndices) = 0;
+            firstOrderViscosity = obj.Mfirst*abs(obj.waveSpeed(v));
+            residualViscosity = obj.fRes(v,dvdt);
+            viscosity = min(firstOrderViscosity, obj.fRes(v,dvdt));
+        end
+        function res = maxResidualNeighbors(obj,v,dvdt)
+            res = obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v);
+            resMat = grid.funcToMatrix(obj.grid,res);
+            res = reshape(max(movmax(resMat,3,1),movmax(resMat,3,2))',obj.grid.N(),1);
         end
     end
     methods (Static)
-
-        function boundaryind = getBoundaryIndices(g)
-            ind = grid.funcToMatrix(g, 1:g.N()); 
-            switch g.D()
-            case 1
-                boundaryind = [ind(1) ind(end)];
-            case 2
-                boundaryind = [ind(1,:) ind(end,:), ind(:,1)' ind(:,end)'];
-            case 3
-                error;
-            end
-        end
-
         function minmaxDiff = minmaxDiffNeighborhood1d(u)
             umax = movmax(u,3);
             umin = movmin(u,3);
             minmaxDiff = umax - umin;
         end
-        function minmaxDiff = minmaxDiffNeighborhood2d(g, u)
-            m = g.m();
-            uMatrix = grid.funcToMatrix(g,u);
+        
+        function minmaxDiff = minmaxDiffNeighborhood2d(grid, u)
+            uMatrix = grid.funcToMatrix(grid,u);
             umax = max(movmax(uMatrix,3,1),movmax(uMatrix,3,2));
             umin = min(movmin(uMatrix,3,1),movmin(uMatrix,3,2));
             minmaxDiff = umax - umin;
-            minmaxDiff = reshape(minmaxDiff',m(1)*m(2),1);
+            minmaxDiff = reshape(minmaxDiff',g.N(),1);
         end
     end
 end
\ No newline at end of file