diff --git a/NumLib/Function/Interpolation.h b/NumLib/Function/Interpolation.h
index 8999d16d0b49b5d9c93e7e3567872fd5acc0fcc2..876ac7d8ee0205534b69fa5d190277c49a4ed1ab 100644
--- a/NumLib/Function/Interpolation.h
+++ b/NumLib/Function/Interpolation.h
@@ -12,6 +12,9 @@
 #include<array>
 #include<cassert>
 
+#include "NumLib/Fem/CoordinatesMapping/NaturalNodeCoordinates.h"
+#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+
 namespace NumLib
 {
 
@@ -85,4 +88,52 @@ void shapeFunctionInterpolate(
                                         interpolated_values...);
 }
 
+/// Interpolates \c node_values given in lower order element nodes (e.g. the
+/// base nodes) to higher order element's nodes (e.g. quadratic nodes) and
+/// writes the result into the global property vector.
+///
+/// The base nodes' values are copied.  For each higher order node the shape
+/// matrices are evaluated for the lower order element (the base nodes), and
+/// used for the the scalar quantity interpolation.
+template <typename LowerOrderShapeFunction, typename HigherOrderMeshElementType,
+          int GlobalDim, typename EigenMatrixType>
+void interpolateToHigherOrderNodes(
+    MeshLib::Element const& element, bool const is_axially_symmetric,
+    Eigen::MatrixBase<EigenMatrixType> const& node_values,
+    MeshLib::PropertyVector<double>& interpolated_values_global_vector)
+{
+    using SF = LowerOrderShapeFunction;
+    using ShapeMatricesType = ShapeMatrixPolicyType<SF, GlobalDim>;
+    using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
+    using FemType = TemplateIsoparametric<SF, ShapeMatricesType>;
+
+    FemType fe(*static_cast<const typename SF::MeshElement*>(&element));
+    int const number_base_nodes = element.getNumberOfBaseNodes();
+    int const number_all_nodes = element.getNumberOfNodes();
+
+    // Copy the values for linear nodes.
+    for (int n = 0; n < number_base_nodes; ++n)
+    {
+        std::size_t const global_index = element.getNodeIndex(n);
+        interpolated_values_global_vector[global_index] = node_values[n];
+    }
+
+    // Shape matrices storage reused in the interpolation loop.
+    ShapeMatrices shape_matrices{SF::DIM, GlobalDim, SF::NPOINTS};
+
+    // Interpolate values for higher order nodes.
+    for (int n = number_base_nodes; n < number_all_nodes; ++n)
+    {
+        // Evaluated at higher order nodes' coordinates.
+        fe.template computeShapeFunctions<ShapeMatrixType::N>(
+            NaturalCoordinates<HigherOrderMeshElementType>::coordinates[n]
+                .data(),
+            shape_matrices, GlobalDim, is_axially_symmetric);
+
+        std::size_t const global_index = element.getNodeIndex(n);
+        interpolated_values_global_vector[global_index] =
+            shape_matrices.N * node_values;
+    }
+}
+
 } // namespace NumLib
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
index 0442978c6b54bf992141524549ff3ff6c41e0c40..f386dc0ddd57634c10b151ee9e1e543d91707e1e 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
@@ -15,7 +15,6 @@
 
 #include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
 #include "MathLib/KelvinVector.h"
-#include "NumLib/Fem/CoordinatesMapping/NaturalNodeCoordinates.h"
 #include "NumLib/Function/Interpolation.h"
 #include "ProcessLib/CoupledSolutionsForStaggeredScheme.h"
 
@@ -23,53 +22,6 @@ namespace ProcessLib
 {
 namespace RichardsMechanics
 {
-/// For each higher order node evaluate the shape matrices for the lower order
-/// element (the base nodes)
-template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
-          int DisplacementDim, typename EigenMatrixType>
-void interpolate_p(MeshLib::Element const& element,
-                   bool const is_axially_symmetric,
-                   Eigen::MatrixBase<EigenMatrixType> const& p_L,
-                   MeshLib::PropertyVector<double>& pressure_interpolated)
-{
-    using ShapeMatricesTypePressure =
-        ShapeMatrixPolicyType<ShapeFunctionPressure, DisplacementDim>;
-    using FemType = NumLib::TemplateIsoparametric<ShapeFunctionPressure,
-                                                  ShapeMatricesTypePressure>;
-
-    FemType fe(*static_cast<const typename ShapeFunctionPressure::MeshElement*>(
-        &element));
-    int const number_base_nodes = element.getNumberOfBaseNodes();
-    int const number_all_nodes = element.getNumberOfNodes();
-
-    // Copy the values for linear nodes.
-    for (int n = 0; n < number_base_nodes; ++n)
-    {
-        std::size_t const global_index = element.getNodeIndex(n);
-        pressure_interpolated[global_index] = p_L[n];
-    }
-
-    // Interpolate values for higher order nodes.
-    for (int n = number_base_nodes; n < number_all_nodes; ++n)
-    {
-        // Evaluated at higher order nodes' coordinates.
-        typename ShapeMatricesTypePressure::ShapeMatrices shape_matrices_p{
-            ShapeFunctionPressure::DIM, DisplacementDim,
-            ShapeFunctionPressure::NPOINTS};
-
-        fe.template computeShapeFunctions<NumLib::ShapeMatrixType::N>(
-            NumLib::NaturalCoordinates<
-                typename ShapeFunctionDisplacement::MeshElement>::coordinates[n]
-                .data(),
-            shape_matrices_p, DisplacementDim, is_axially_symmetric);
-
-        auto const& N_p = shape_matrices_p.N;
-
-        std::size_t const global_index = element.getNodeIndex(n);
-        pressure_interpolated[global_index] = N_p * p_L;
-    }
-}
-
 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
           typename IntegrationMethod, int DisplacementDim>
 RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
@@ -881,9 +833,10 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
 
     (*_process_data.element_saturation)[_element.getID()] = saturation_avg;
 
-    interpolate_p<ShapeFunctionDisplacement, ShapeFunctionPressure,
-                  DisplacementDim>(_element, _is_axially_symmetric, p_L,
-                                   *_process_data.pressure_interpolated);
+    NumLib::interpolateToHigherOrderNodes<
+        ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
+        DisplacementDim>(_element, _is_axially_symmetric, p_L,
+                         *_process_data.pressure_interpolated);
 }
 
 }  // namespace RichardsMechanics