diff --git a/dumux/implicit/1p2c/1p2cfluxvariables.hh b/dumux/implicit/1p2c/1p2cfluxvariables.hh
index ff2a651a2916a7dd51cac9e6449adbeb08ef2009..7f4eefda5e1ea88886728e120b13f22d89cfadbd 100644
--- a/dumux/implicit/1p2c/1p2cfluxvariables.hh
+++ b/dumux/implicit/1p2c/1p2cfluxvariables.hh
@@ -293,58 +293,37 @@ protected:
         const VolumeVariables &volVarsI = elemVolVars[face().i];
         const VolumeVariables &volVarsJ = elemVolVars[face().j];
 
-        GlobalPosition tmp;
-        if (!problem.spatialParams().useTwoPointGradient(element, face().i, face().j)) {
-            // use finite-element gradients
-            tmp = 0.0;
-            for (unsigned int idx = 0;
-                    idx < face().numFap;
-                    idx++) // loop over adjacent vertices
-            {
-                // FE gradient at vertex idx
-                const GlobalPosition &feGrad = face().grad[idx];
-
-                // index for the element volume variables 
-                int volVarsIdx = face().fapIndices[idx];
-
-                // the pressure gradient
-                tmp = feGrad;
-                tmp *= elemVolVars[volVarsIdx].pressure();
-                potentialGrad_ += tmp;
-
-                // the mole-fraction gradient
-                tmp = feGrad;
-                tmp *= elemVolVars[volVarsIdx].moleFraction(transportCompIdx);
-                moleFractionGrad_ += tmp;
-
-                // phase viscosity
-                viscosity_ += elemVolVars[volVarsIdx].viscosity()*face().shapeValue[idx];
-
-                //phase molar density
-                molarDensity_ += elemVolVars[volVarsIdx].molarDensity()*face().shapeValue[idx];
-
-                //phase density
-                density_ += elemVolVars[volVarsIdx].density()*face().shapeValue[idx];
-            }
-        }
-        else {
-            // use two-point gradients
-            const auto geometry = element.geometry();
-            const GlobalPosition &globalPosI = geometry.corner(face().i);
-            const GlobalPosition &globalPosJ = geometry.corner(face().j);
-            tmp = globalPosI;
-            tmp -= globalPosJ;
-            Scalar dist = tmp.two_norm();
-
-            tmp = face().normal;
-            tmp /= face().normal.two_norm()*dist;
-
-            potentialGrad_ = tmp;
-            potentialGrad_ *= volVarsJ.pressure() - volVarsI.pressure();
-            moleFractionGrad_ = tmp;
-            moleFractionGrad_ *= volVarsJ.moleFraction(transportCompIdx) - volVarsI.moleFraction(transportCompIdx);
+        GlobalPosition tmp(0.0);
+        // use finite-element gradients
+        for (unsigned int idx = 0; idx < face().numFap; idx++) // loop over adjacent vertices
+        {
+            // FE gradient at vertex idx
+            const GlobalPosition &feGrad = face().grad[idx];
+
+            // index for the element volume variables 
+            int volVarsIdx = face().fapIndices[idx];
+
+            // the pressure gradient
+            tmp = feGrad;
+            tmp *= elemVolVars[volVarsIdx].pressure();
+            potentialGrad_ += tmp;
+
+            // the mole-fraction gradient
+            tmp = feGrad;
+            tmp *= elemVolVars[volVarsIdx].moleFraction(transportCompIdx);
+            moleFractionGrad_ += tmp;
+
+            // phase viscosity
+            viscosity_ += elemVolVars[volVarsIdx].viscosity()*face().shapeValue[idx];
+
+            //phase molar density
+            molarDensity_ += elemVolVars[volVarsIdx].molarDensity()*face().shapeValue[idx];
+
+            //phase density
+            density_ += elemVolVars[volVarsIdx].density()*face().shapeValue[idx];
         }
 
+
         ///////////////
         // correct the pressure gradients by the gravitational acceleration
         ///////////////
diff --git a/test/geomechanics/el1p2c/el1p2cspatialparams.hh b/test/geomechanics/el1p2c/el1p2cspatialparams.hh
index 3d491ba0ce64e0bfb6e2f1fc3c8b9bbd618a17ea..6d65645a96cb1e9ef088b5d948bd70a2f805585b 100644
--- a/test/geomechanics/el1p2c/el1p2cspatialparams.hh
+++ b/test/geomechanics/el1p2c/el1p2cspatialparams.hh
@@ -185,21 +185,6 @@ public:
         return param;
     }
 
-    /*!
-     * \brief Define if the model should apply two-point approximation 
-     * instead of a box approximation for the fluxes
-     *
-     * \param element The finite element
-     * \param vertexI first point for the two-point flux approximation
-     * \param vertexJ second point for the two-point flux approximation
-     */
-    bool useTwoPointGradient(const Element &element,
-                             int vertexI,
-                             int vertexJ) const
-    {
-        return false;
-    }
-
     /*!
      * \brief Return dispersivity (not needed here
      *
diff --git a/test/implicit/1p/1pnispatialparams.hh b/test/implicit/1p/1pnispatialparams.hh
index 42f4824c11c7c32007275604f9b45d9daf254c83..19895fe0cda9cc32735448bef4cfa83d118a4401 100644
--- a/test/implicit/1p/1pnispatialparams.hh
+++ b/test/implicit/1p/1pnispatialparams.hh
@@ -116,13 +116,6 @@ public:
         return 0;
     }
 
-    bool useTwoPointGradient(const Element &element,
-                             const int vertexI,
-                             const int vertexJ) const
-    {
-        return false;
-    }
-
     /*!
      * \brief Returns the heat capacity \f$[J / (kg K)]\f$ of the rock matrix.
      *
diff --git a/test/implicit/1p2c/1p2cnispatialparams.hh b/test/implicit/1p2c/1p2cnispatialparams.hh
index e771c7ade9c97f2212e4b4b7156f7e46b5f8e8d1..07206f29c12824c0b0844c5e46e1ec047edd1bf9 100644
--- a/test/implicit/1p2c/1p2cnispatialparams.hh
+++ b/test/implicit/1p2c/1p2cnispatialparams.hh
@@ -116,13 +116,6 @@ public:
         return 0;
     }
 
-    bool useTwoPointGradient(const Element &element,
-                             const int vertexI,
-                             const int vertexJ) const
-    {
-        return false;
-    }
-
     /*!
      * \brief Returns the heat capacity \f$[J / (kg K)]\f$ of the rock matrix.
      *
diff --git a/test/implicit/1p2c/1p2coutflowspatialparams.hh b/test/implicit/1p2c/1p2coutflowspatialparams.hh
index deefb643f57c167d8d2dabd0be1fa67f93ca5868..e329be316ae8ccb938f1147c04002d0ae8dbae9c 100644
--- a/test/implicit/1p2c/1p2coutflowspatialparams.hh
+++ b/test/implicit/1p2c/1p2coutflowspatialparams.hh
@@ -118,13 +118,6 @@ public:
         return 0;
     }
 
-    bool useTwoPointGradient(const Element &element,
-                             const int vertexI,
-                             const int vertexJ) const
-    {
-        return false;
-    }
-
     /*!
      * \brief Returns the heat capacity \f$[J / (kg K)]\f$ of the rock matrix.
      *
diff --git a/test/implicit/3p/3pnispatialparams.hh b/test/implicit/3p/3pnispatialparams.hh
index 17c6695fcd030ed0b640af11ba0e475bb08217b2..a1e4e6ab9c48841868c3199afcab435a851c2708 100644
--- a/test/implicit/3p/3pnispatialparams.hh
+++ b/test/implicit/3p/3pnispatialparams.hh
@@ -158,14 +158,6 @@ public:
 		return materialParams_;
     }
 
-
-    bool useTwoPointGradient(const Element &element,
-                             const int vertexI,
-                             const int vertexJ) const
-    {
-        return false;
-    }
-
     /*!
      * \brief Returns the heat capacity \f$[J / (kg K)]\f$ of the rock matrix.
      *
diff --git a/test/implicit/richards/richardsnispatialparams.hh b/test/implicit/richards/richardsnispatialparams.hh
index 017e393693afb0efc4232024d745eb5d62a2fb36..be724a4b1165d63f613099342a6945692d222855 100644
--- a/test/implicit/richards/richardsnispatialparams.hh
+++ b/test/implicit/richards/richardsnispatialparams.hh
@@ -165,14 +165,6 @@ public:
         return materialParams_;
     }
 
-
-    bool useTwoPointGradient(const Element &element,
-                             const int vertexI,
-                             const int vertexJ) const
-    {
-        return false;
-    }
-
     /*!
      * \brief Returns the heat capacity \f$[J / (kg K)]\f$ of the rock matrix.
      *