diff --git a/ProcessLib/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
index b2e870bab5e63aa533d29c4f0c7fe4b49f5b3b48..0ef541b0c6d5a2e8bb9ba9c0e33e33a6555ab63b 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
@@ -285,6 +285,8 @@ public:
             auto const& N = ip_data.N;
             auto const& dNdx = ip_data.dNdx;
 
+            pos.setIntegrationPoint(ip);
+
             auto const& K =
                 _process_data.porous_media_properties.getIntrinsicPermeability(
                     t, pos);
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
index 3748d45dc23261d3214766959d59018219bce645..00da123f9aba0d585eeb6a032a99bf4a34967e90 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
@@ -15,6 +15,7 @@
 #include "MaterialLib/SolidModels/KelvinVector.h"
 #include "MaterialLib/SolidModels/LinearElasticIsotropic.h"
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "NumLib/DOF/DOFTableUtil.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "ProcessLib/Deformation/BMatrixPolicy.h"
 #include "ProcessLib/Deformation/LinearBMatrix.h"
@@ -364,43 +365,6 @@ public:
         }
     }
 
-    void postTimestepConcrete(std::vector<double> const& local_x) override
-    {
-        double const& t = _process_data.t;
-
-        auto p =
-            Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
-                pressure_size> const>(local_x.data() + pressure_index,
-                                      pressure_size);
-        using GlobalDimVectorType =
-            typename ShapeMatricesTypePressure::GlobalDimVectorType;
-
-        unsigned const n_integration_points =
-            _integration_method.getNumberOfPoints();
-
-        SpatialPosition x_position;
-        x_position.setElementID(_element.getID());
-        for (unsigned ip = 0; ip < n_integration_points; ip++)
-        {
-            x_position.setIntegrationPoint(ip);
-            double const K_over_mu =
-                _process_data.intrinsic_permeability(t, x_position)[0] /
-                _process_data.fluid_viscosity(t, x_position)[0];
-
-            auto const rho_fr = _process_data.fluid_density(t, x_position)[0];
-            auto const& b = _process_data.specific_body_force;
-
-            // Compute the velocity
-            auto const& dNdx_p = _ip_data[ip].dNdx_p;
-            GlobalDimVectorType const darcy_velocity =
-                -K_over_mu * dNdx_p * p - K_over_mu * rho_fr * b;
-            for (unsigned d = 0; d < DisplacementDim; ++d)
-            {
-                _darcy_velocities[d][ip] = darcy_velocity[d];
-            }
-        }
-    }
-
     Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
         const unsigned integration_point) const override
     {
@@ -522,34 +486,54 @@ public:
         return getIntPtEpsilon(cache, 5);
     }
 
-    std::vector<double> const& getIntPtDarcyVelocityX(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& /*cache*/) const override
+    std::vector<double> const& getIntPtDarcyVelocity(
+        const double t,
+        GlobalVector const& current_solution,
+        NumLib::LocalToGlobalIndexMap const& dof_table,
+        std::vector<double>& cache) const override
     {
-        assert(!_darcy_velocities.empty());
-        return _darcy_velocities[0];
-    }
+        auto const num_intpts = _ip_data.size();
 
-    std::vector<double> const& getIntPtDarcyVelocityY(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& /*cache*/) const override
-    {
-        assert(_darcy_velocities.size() > 1);
-        return _darcy_velocities[1];
-    }
+        auto const indices = NumLib::getIndices(_element.getID(), dof_table);
+        assert(!indices.empty());
+        auto const local_x = current_solution.get(indices);
 
-    std::vector<double> const& getIntPtDarcyVelocityZ(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& /*cache*/) const override
-    {
-        assert(_darcy_velocities.size() > 2);
-        return _darcy_velocities[2];
+        cache.clear();
+        auto cache_vec = MathLib::createZeroedMatrix<Eigen::Matrix<
+            double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
+            cache, DisplacementDim, num_intpts);
+
+        SpatialPosition pos;
+        pos.setElementID(_element.getID());
+
+        auto p =
+            Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
+                pressure_size> const>(local_x.data() + pressure_index,
+                                      pressure_size);
+        using GlobalDimVectorType =
+            typename ShapeMatricesTypePressure::GlobalDimVectorType;
+
+        unsigned const n_integration_points =
+            _integration_method.getNumberOfPoints();
+
+        SpatialPosition x_position;
+        x_position.setElementID(_element.getID());
+        for (unsigned ip = 0; ip < n_integration_points; ip++) {
+            x_position.setIntegrationPoint(ip);
+            double const K_over_mu =
+                _process_data.intrinsic_permeability(t, x_position)[0] /
+                _process_data.fluid_viscosity(t, x_position)[0];
+
+            auto const rho_fr = _process_data.fluid_density(t, x_position)[0];
+            auto const& b = _process_data.specific_body_force;
+
+            // Compute the velocity
+            auto const& dNdx_p = _ip_data[ip].dNdx_p;
+            cache_vec.col(ip).noalias() =
+                -K_over_mu * dNdx_p * p - K_over_mu * rho_fr * b;
+        }
+
+        return cache;
     }
 
 private:
@@ -602,11 +586,6 @@ private:
         typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType>
         _secondary_data;
 
-    std::vector<std::vector<double>> _darcy_velocities =
-        std::vector<std::vector<double>>(
-            DisplacementDim,
-            std::vector<double>(_integration_method.getNumberOfPoints()));
-
     static const int pressure_index = 0;
     static const int pressure_size = ShapeFunctionPressure::NPOINTS;
     static const int displacement_index = ShapeFunctionPressure::NPOINTS;
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsProcess-impl.h b/ProcessLib/HydroMechanics/HydroMechanicsProcess-impl.h
index ee8477b123b1e1fdbd73e12cce1670e5daa4e90f..40d1565e622eba6439beb9de21a1b3ded2c22e9c 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsProcess-impl.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsProcess-impl.h
@@ -164,22 +164,10 @@ void HydroMechanicsProcess<DisplacementDim>::initializeConcreteProcess(
             &LocalAssemblerInterface::getIntPtEpsilonXY));
 
     Base::_secondary_variables.addSecondaryVariable(
-        "velocity_x",
+        "velocity",
         makeExtrapolator(
             1, getExtrapolator(), _local_assemblers,
-            &LocalAssemblerInterface::getIntPtDarcyVelocityX));
-
-    Base::_secondary_variables.addSecondaryVariable(
-        "velocity_y",
-        makeExtrapolator(
-            1, getExtrapolator(), _local_assemblers,
-            &LocalAssemblerInterface::getIntPtDarcyVelocityY));
-
-    Base::_secondary_variables.addSecondaryVariable(
-        "velocity_z",
-        makeExtrapolator(
-            1, getExtrapolator(), _local_assemblers,
-            &LocalAssemblerInterface::getIntPtDarcyVelocityZ));
+            &LocalAssemblerInterface::getIntPtDarcyVelocity));
 }
 
 template <int DisplacementDim>
diff --git a/ProcessLib/HydroMechanics/LocalAssemblerInterface.h b/ProcessLib/HydroMechanics/LocalAssemblerInterface.h
index 95b48de589c1fa372eaa82dd0c4db754e808d480..f3dd1f9a6e40532db269ebed5eb7f8d46e838633 100644
--- a/ProcessLib/HydroMechanics/LocalAssemblerInterface.h
+++ b/ProcessLib/HydroMechanics/LocalAssemblerInterface.h
@@ -91,19 +91,7 @@ struct LocalAssemblerInterface
         NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
         std::vector<double>& cache) const = 0;
 
-    virtual std::vector<double> const& getIntPtDarcyVelocityX(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& cache) const = 0;
-
-    virtual std::vector<double> const& getIntPtDarcyVelocityY(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& cache) const = 0;
-
-    virtual std::vector<double> const& getIntPtDarcyVelocityZ(
+    virtual std::vector<double> const& getIntPtDarcyVelocity(
         const double /*t*/,
         GlobalVector const& /*current_solution*/,
         NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
index b2c89ddc270d28aa48e2b0b5195f5f7ba0b5160f..d726b76bd1999936036e40aa674d6f774fdfe3a1 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
@@ -19,6 +19,7 @@
 
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "NumLib/Extrapolation/ExtrapolatableElement.h"
+#include "NumLib/DOF/DOFTableUtil.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "ProcessLib/LocalAssemblerInterface.h"
@@ -39,19 +40,7 @@ class LiquidFlowLocalAssemblerInterface
       public NumLib::ExtrapolatableElement
 {
 public:
-    virtual std::vector<double> const& getIntPtDarcyVelocityX(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& /*cache*/) const = 0;
-
-    virtual std::vector<double> const& getIntPtDarcyVelocityY(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& /*cache*/) const = 0;
-
-    virtual std::vector<double> const& getIntPtDarcyVelocityZ(
+    virtual std::vector<double> const& getIntPtDarcyVelocity(
         const double /*t*/,
         GlobalVector const& /*current_solution*/,
         NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
@@ -122,34 +111,43 @@ public:
         return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
     }
 
-    std::vector<double> const& getIntPtDarcyVelocityX(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& /*cache*/) const override
-    {
-        assert(!_darcy_velocities.empty());
-        return _darcy_velocities[0];
-    }
-
-    std::vector<double> const& getIntPtDarcyVelocityY(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& /*cache*/) const override
-    {
-        assert(_darcy_velocities.size() > 1);
-        return _darcy_velocities[1];
-    }
-
-    std::vector<double> const& getIntPtDarcyVelocityZ(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& /*cache*/) const override
+    std::vector<double> const& getIntPtDarcyVelocity(
+        const double t,
+        GlobalVector const& current_solution,
+        NumLib::LocalToGlobalIndexMap const& dof_table,
+        std::vector<double>& cache) const override
     {
-        assert(_darcy_velocities.size() > 2);
-        return _darcy_velocities[2];
+        // auto const num_nodes = ShapeFunction_::NPOINTS;
+        auto const num_intpts = _shape_matrices.size();
+
+        auto const indices = NumLib::getIndices(
+            _element.getID(), dof_table);
+        assert(!indices.empty());
+        auto const local_x = current_solution.get(indices);
+        auto const local_x_vec =
+            MathLib::toVector<Eigen::Matrix<double, ShapeFunction::NPOINTS, 1>>(
+                local_x, ShapeFunction::NPOINTS);
+
+        cache.clear();
+        auto cache_vec = MathLib::createZeroedMatrix<
+            Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
+            cache, GlobalDim, num_intpts);
+
+        SpatialPosition pos;
+        pos.setElementID(_element.getID());
+
+        // TODO fix
+#if 0
+        for (unsigned i = 0; i < num_intpts; ++i) {
+            pos.setIntegrationPoint(i);
+            auto const k = _process_data.hydraulic_conductivity(t, pos)[0];
+            // dimensions: (d x 1) = (d x n) * (n x 1)
+            cache_vec.col(i).noalias() =
+                -k * _shape_matrices[i].dNdx * local_x_vec;
+        }
+#endif
+
+        return cache;
     }
 
 private:
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
index 68a7e6552d5330857b98bc750fc66abfd17a004c..beacb9f51b7e3f5f22b9df9af282af69ab1682ab 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
@@ -69,27 +69,10 @@ void LiquidFlowProcess::initializeConcreteProcess(
         *_material_properties);
 
     _secondary_variables.addSecondaryVariable(
-        "darcy_velocity_x",
-        makeExtrapolator(
-            1, getExtrapolator(), _local_assemblers,
-            &LiquidFlowLocalAssemblerInterface::getIntPtDarcyVelocityX));
-
-    if (mesh.getDimension() > 1)
-    {
-        _secondary_variables.addSecondaryVariable(
-            "darcy_velocity_y",
-            makeExtrapolator(
-                1, getExtrapolator(), _local_assemblers,
-                &LiquidFlowLocalAssemblerInterface::getIntPtDarcyVelocityY));
-    }
-    if (mesh.getDimension() > 2)
-    {
-        _secondary_variables.addSecondaryVariable(
-            "darcy_velocity_z",
-            makeExtrapolator(
-                1, getExtrapolator(), _local_assemblers,
-                &LiquidFlowLocalAssemblerInterface::getIntPtDarcyVelocityZ));
-    }
+        "darcy_velocity",
+        makeExtrapolator(mesh.getDimension(),
+            getExtrapolator(), _local_assemblers,
+            &LiquidFlowLocalAssemblerInterface::getIntPtDarcyVelocity));
 }
 
 void LiquidFlowProcess::assembleConcreteProcess(
diff --git a/ProcessLib/RichardsFlow/RichardsFlowFEM.h b/ProcessLib/RichardsFlow/RichardsFlowFEM.h
index c21934444c9b5c666940d7924f2cecfdda40f9ce..53c3b1c5eb78e3a1e49559326496e49e2af16b8c 100644
--- a/ProcessLib/RichardsFlow/RichardsFlowFEM.h
+++ b/ProcessLib/RichardsFlow/RichardsFlowFEM.h
@@ -14,6 +14,7 @@
 #include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "MeshLib/CoordinateSystem.h"
+#include "NumLib/DOF/DOFTableUtil.h"
 #include "NumLib/Extrapolation/ExtrapolatableElement.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
@@ -63,19 +64,7 @@ public:
         NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
         std::vector<double>& /*cache*/) const = 0;
 
-    virtual std::vector<double> const& getIntPtDarcyVelocityX(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& /*cache*/) const = 0;
-
-    virtual std::vector<double> const& getIntPtDarcyVelocityY(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& /*cache*/) const = 0;
-
-    virtual std::vector<double> const& getIntPtDarcyVelocityZ(
+    virtual std::vector<double> const& getIntPtDarcyVelocity(
         const double /*t*/,
         GlobalVector const& /*current_solution*/,
         NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
@@ -110,9 +99,6 @@ public:
           _process_data(process_data),
           _integration_method(integration_order),
           _saturation(
-              std::vector<double>(_integration_method.getNumberOfPoints())),
-          _darcy_velocities(
-              GlobalDim,
               std::vector<double>(_integration_method.getNumberOfPoints()))
     {
         // This assertion is valid only if all nodal d.o.f. use the same shape
@@ -240,9 +226,44 @@ public:
         }  // end of mass lumping
     }
 
-    void computeSecondaryVariableConcrete(
-        double const t, std::vector<double> const& local_x) override
+    Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
+        const unsigned integration_point) const override
     {
+        auto const& N = _ip_data[integration_point].N;
+
+        // assumes N is stored contiguously in memory
+        return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
+    }
+
+    std::vector<double> const& getIntPtSaturation(
+        const double /*t*/,
+        GlobalVector const& /*current_solution*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::vector<double>& /*cache*/) const override
+    {
+        assert(!_saturation.empty());
+        return _saturation;
+    }
+
+    std::vector<double> const& getIntPtDarcyVelocity(
+        const double t,
+        GlobalVector const& current_solution,
+        NumLib::LocalToGlobalIndexMap const& dof_table,
+        std::vector<double>& cache) const override
+    {
+        // auto const num_nodes = ShapeFunction_::NPOINTS;
+        auto const num_intpts = _shape_matrices.size();
+
+        auto const indices = NumLib::getIndices(
+            _element.getID(), dof_table);
+        assert(!indices.empty());
+        auto const local_x = current_solution.get(indices);
+
+        cache.clear();
+        auto cache_vec = MathLib::createZeroedMatrix<
+            Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
+            cache, GlobalDim, num_intpts);
+
         SpatialPosition pos;
         pos.setElementID(_element.getID());
         const int material_id =
@@ -278,7 +299,7 @@ public:
             auto const mu = _process_data.material->getFluidViscosity(
                 p_int_pt, temperature);
             auto const K_mat_coeff = permeability * (k_rel / mu);
-            GlobalDimVectorType velocity =
+            cache_vec.col(ip).noalias() =
                 -K_mat_coeff * _ip_data[ip].dNdx * p_nodal_values;
             if (_process_data.has_gravity)
             {
@@ -288,63 +309,11 @@ public:
                 assert(body_force.size() == GlobalDim);
                 // here it is assumed that the vector body_force is directed
                 // 'downwards'
-                velocity += K_mat_coeff * rho_w * body_force;
-            }
-
-            for (unsigned d = 0; d < GlobalDim; ++d)
-            {
-                _darcy_velocities[d][ip] = velocity[d];
+                cache_vec.col(ip).noalias() += K_mat_coeff * rho_w * body_force;
             }
         }
-    }
-
-    Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
-        const unsigned integration_point) const override
-    {
-        auto const& N = _ip_data[integration_point].N;
-
-        // assumes N is stored contiguously in memory
-        return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
-    }
 
-    std::vector<double> const& getIntPtSaturation(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& /*cache*/) const override
-    {
-        assert(!_saturation.empty());
-        return _saturation;
-    }
-
-    std::vector<double> const& getIntPtDarcyVelocityX(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& /*cache*/) const override
-    {
-        assert(!_darcy_velocities.empty());
-        return _darcy_velocities[0];
-    }
-
-    std::vector<double> const& getIntPtDarcyVelocityY(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& /*cache*/) const override
-    {
-        assert(_darcy_velocities.size() > 1);
-        return _darcy_velocities[1];
-    }
-
-    std::vector<double> const& getIntPtDarcyVelocityZ(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& /*cache*/) const override
-    {
-        assert(_darcy_velocities.size() > 2);
-        return _darcy_velocities[2];
+        return cache;
     }
 
 private:
@@ -361,7 +330,6 @@ private:
             NodalRowVectorType, GlobalDimNodalMatrixType, NodalMatrixType>>>
         _ip_data;
     std::vector<double> _saturation;
-    std::vector<std::vector<double>> _darcy_velocities;
 };
 
 }  // namespace RichardsFlow
diff --git a/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp b/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
index 52372d654e5f211e23dd0d85647b2816f785ea82..612dc11b1251953f324334d80dd488af8304d6f9 100644
--- a/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
+++ b/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
@@ -55,27 +55,10 @@ void RichardsFlowProcess::initializeConcreteProcess(
             &RichardsFlowLocalAssemblerInterface::getIntPtSaturation));
 
     _secondary_variables.addSecondaryVariable(
-        "darcy_velocity_x",
+        "darcy_velocity",
         makeExtrapolator(
-            1, getExtrapolator(), _local_assemblers,
-            &RichardsFlowLocalAssemblerInterface::getIntPtDarcyVelocityX));
-
-    if (mesh.getDimension() > 1)
-    {
-        _secondary_variables.addSecondaryVariable(
-            "darcy_velocity_y",
-            makeExtrapolator(
-                1, getExtrapolator(), _local_assemblers,
-                &RichardsFlowLocalAssemblerInterface::getIntPtDarcyVelocityY));
-    }
-    if (mesh.getDimension() > 2)
-    {
-        _secondary_variables.addSecondaryVariable(
-            "darcy_velocity_z",
-            makeExtrapolator(
-                1, getExtrapolator(), _local_assemblers,
-                &RichardsFlowLocalAssemblerInterface::getIntPtDarcyVelocityZ));
-    }
+            mesh.getDimension(), getExtrapolator(), _local_assemblers,
+            &RichardsFlowLocalAssemblerInterface::getIntPtDarcyVelocity));
 }
 
 void RichardsFlowProcess::assembleConcreteProcess(
@@ -107,15 +90,5 @@ void RichardsFlowProcess::assembleWithJacobianConcreteProcess(
         dx_dx, M, K, b, Jac, _coupling_term);
 }
 
-void RichardsFlowProcess::computeSecondaryVariableConcrete(
-    double const t, GlobalVector const& x)
-{
-    DBUG("Compute the Darcy velocity for RichardsFlowProcess");
-    GlobalExecutor::executeMemberOnDereferenced(
-        &RichardsFlowLocalAssemblerInterface::computeSecondaryVariable,
-        _local_assemblers, *_local_to_global_index_map, t, x,
-        _coupling_term);
-}
-
 }  // namespace RichardsFlow
 }  // namespace ProcessLib
diff --git a/ProcessLib/RichardsFlow/RichardsFlowProcess.h b/ProcessLib/RichardsFlow/RichardsFlowProcess.h
index 939b9132692d004dbed1ff971620d6f7f797bae8..9f59e799854d7c6554c2704a0c578571ee788070 100644
--- a/ProcessLib/RichardsFlow/RichardsFlowProcess.h
+++ b/ProcessLib/RichardsFlow/RichardsFlowProcess.h
@@ -46,9 +46,6 @@ public:
     bool isLinear() const override { return false; }
     //! @}
 
-    void computeSecondaryVariableConcrete(
-        double const t, GlobalVector const& x) override;
-
 private:
     void initializeConcreteProcess(
         NumLib::LocalToGlobalIndexMap const& dof_table,