--- a/+rv/ResidualViscosity.m	Mon Feb 18 09:00:00 2019 +0100
+++ b/+rv/ResidualViscosity.m	Tue Mar 05 10:53:34 2019 +0100
@@ -24,6 +24,12 @@
             obj.boundaryIndices = obj.getBoundaryIndices(g);
         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;
+        end
+
         function [viscosity, Df, firstOrderViscosity, residualViscosity] = evaluate(obj, v, dvdt)
             Df = obj.Df(v);
             firstOrderViscosity = obj.Cmax*obj.h*abs(obj.waveSpeed(v));
@@ -47,15 +53,6 @@
             end
         end
 
-        % function f = minmaxDiffNeighborhood(g)
-        %     switch g.D()
-        %         case 1
-        %             f = @(u)minmaxDiffNeighborhood1d(u);
-        %         case 2
-        %             f = @(u)minmaxDiffNeighborhood2d(g,u);
-        %     end
-        % end
-
         function minmaxDiff = minmaxDiffNeighborhood1d(u)
             umax = movmax(u,3);
             umin = movmin(u,3);