diff --git a/dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh b/dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh
index 928f04308f3630bf2d21088e7326153b618a3b8b..089b699afd44fc9eabb692fe5e28e30707349df8 100644
--- a/dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh
+++ b/dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh
@@ -142,8 +142,8 @@ public:
             }
 
             Dune::FieldVector<Scalar, dim> refVelocity(0);
-            refVelocity[0] = 0.5 * (flux[1] - flux[0]);
-            refVelocity[1] = 0.5 * (flux[3] - flux[2]);
+            for (int i = 0; i < dim; i++)
+                refVelocity[i] = 0.5 * (flux[2*i + 1] - flux[2*i]);
 
             typedef Dune::GenericReferenceElements<Scalar, dim> ReferenceElements;
             const Dune::FieldVector<Scalar, dim>& localPos = ReferenceElements::general(eIt->geometry().type()).position(0,
diff --git a/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh b/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh
index 5b061207e9a1a98c25fb57d7792b15b130fc5d20..02b36c210ea7a9124063dc6d0a00f8c1651df121 100644
--- a/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh
+++ b/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh
@@ -207,11 +207,8 @@ public:
             }
 
             Dune::FieldVector < Scalar, dim > refVelocity(0);
-            refVelocity[0] = 0.5 * (fluxW[1] - fluxW[0]);
-            if (dim > 1)
-            refVelocity[1] = 0.5 * (fluxW[3] - fluxW[2]);
-            if (dim == 3)
-            refVelocity[2] = 0.5 * (fluxW[5] - fluxW[4]);
+            for (int i = 0; i < dim; i++)
+                refVelocity[i] = 0.5 * (fluxW[2*i + 1] - fluxW[2*i]);
 
             const Dune::FieldVector<Scalar, dim>& localPos =
             ReferenceElementContainer::general(eIt->geometry().type()).position(0, 0);
@@ -229,11 +226,8 @@ public:
             velocity[globalIdx] = elementVelocity;
 
             refVelocity = 0;
-            refVelocity[0] = 0.5 * (fluxNW[1] - fluxNW[0]);
-            if (dim > 1)
-            refVelocity[1] = 0.5 * (fluxNW[3] - fluxNW[2]);
-            if (dim == 3)
-            refVelocity[2] = 0.5 * (fluxNW[5] - fluxNW[4]);
+            for (int i = 0; i < dim; i++)
+                refVelocity[i] = 0.5 * (fluxNW[2*i + 1] - fluxNW[2*i]);
 
             // calculate the element velocity by the Piola transformation
             elementVelocity = 0;
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundvelocity2p.hh b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundvelocity2p.hh
index 242111aaee9b148cca2ed5a1307c55c6c92b7a58..0d31e740b01cd91c0baad73bfb36d0030d032043 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundvelocity2p.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundvelocity2p.hh
@@ -185,8 +185,8 @@ public:
             }
 
             DimVector refVelocity(0);
-            refVelocity[0] = 0.5 * (fluxW[1] - fluxW[0]);
-            refVelocity[1] = 0.5 * (fluxW[3] - fluxW[2]);
+            for (int i = 0; i < dim; i++)
+                refVelocity[i] = 0.5 * (fluxW[2*i + 1] - fluxW[2*i]);
 
             const DimVector& localPos =
                     ReferenceElementContainer::general(eIt->geometry().type()).position(0, 0);
@@ -202,8 +202,8 @@ public:
             velocityWetting[globalIdx] = elementVelocity;
 
             refVelocity = 0;
-            refVelocity[0] = 0.5 * (fluxNW[1] - fluxNW[0]);
-            refVelocity[1] = 0.5 * (fluxNW[3] - fluxNW[2]);
+            for (int i = 0; i < dim; i++)
+                refVelocity[i] = 0.5 * (fluxNW[2*i + 1] - fluxNW[2*i]);
 
             // calculate the element velocity by the Piola transformation
             elementVelocity = 0;
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundvelocity2padaptive.hh b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundvelocity2padaptive.hh
index 8c8261788100e67c8003189fc77b4f548e21ff80..387ed571ec183838b54bfbc1cfb2c5edcc84a88d 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundvelocity2padaptive.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundvelocity2padaptive.hh
@@ -191,8 +191,8 @@ public:
              }
 
              DimVector refVelocity(0);
-             refVelocity[0] = 0.5 * (fluxW[1] - fluxW[0]);
-             refVelocity[1] = 0.5 * (fluxW[3] - fluxW[2]);
+            for (int i = 0; i < dim; i++)
+                refVelocity[i] = 0.5 * (fluxW[2*i + 1] - fluxW[2*i]);
 
              const DimVector& localPos =
                      ReferenceElementContainer::general(eIt->geometry().type()).position(0, 0);
@@ -208,8 +208,8 @@ public:
              velocityWetting[globalIdx] = elementVelocity;
 
              refVelocity = 0;
-             refVelocity[0] = 0.5 * (fluxNW[1] - fluxNW[0]);
-             refVelocity[1] = 0.5 * (fluxNW[3] - fluxNW[2]);
+            for (int i = 0; i < dim; i++)
+                refVelocity[i] = 0.5 * (fluxNW[2*i + 1] - fluxNW[2*i]);
 
              // calculate the element velocity by the Piola transformation
              elementVelocity = 0;
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfao2pfaboundvelocity2p.hh b/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfao2pfaboundvelocity2p.hh
index fd14abb154da5e33c3406cc6aafa17d7a62d2699..787fe0eb4669dd12b8d3416a058c4d79156af696 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfao2pfaboundvelocity2p.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfao2pfaboundvelocity2p.hh
@@ -185,8 +185,8 @@ public:
             }
 
             DimVector refVelocity(0);
-            refVelocity[0] = 0.5 * (fluxW[1] - fluxW[0]);
-            refVelocity[1] = 0.5 * (fluxW[3] - fluxW[2]);
+            for (int i = 0; i < dim; i++)
+                refVelocity[i] = 0.5 * (fluxW[2*i + 1] - fluxW[2*i]);
 
             const DimVector& localPos =
                     ReferenceElementContainer::general(eIt->geometry().type()).position(0, 0);
@@ -202,8 +202,8 @@ public:
             velocityWetting[globalIdx] = elementVelocity;
 
             refVelocity = 0;
-            refVelocity[0] = 0.5 * (fluxNW[1] - fluxNW[0]);
-            refVelocity[1] = 0.5 * (fluxNW[3] - fluxNW[2]);
+            for (int i = 0; i < dim; i++)
+                refVelocity[i] = 0.5 * (fluxNW[2*i + 1] - fluxNW[2*i]);
 
             // calculate the element velocity by the Piola transformation
             elementVelocity = 0;