diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
index c7d2fb5f04e0667c5a2b4ca9b7bea8c60b2c5eae..6d09a5688bdbe0036c91c10d66ea4f785a642bc4 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
@@ -15,6 +15,7 @@
 
 #include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
 #include "MathLib/KelvinVector.h"
+#include "NumLib/Fem/CoordinatesMapping/NaturalNodeCoordinates.h"
 #include "NumLib/Function/Interpolation.h"
 #include "ProcessLib/CoupledSolutionsForStaggeredScheme.h"
 
@@ -832,6 +833,47 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
     saturation_avg /= n_integration_points;
 
     (*_process_data.element_saturation)[_element.getID()] = saturation_avg;
+
+    // For each higher order node evaluate the shape matrices for the lower
+    // order element (the base nodes)
+    // TODO (naumov) Extract this method to be useful for other processes.
+    auto interpolate_p = [&]() {
+        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();
+
+        for (int n = 0; n < number_base_nodes; ++n)
+        {
+            std::size_t const global_index = _element.getNodeIndex(n);
+            (*_process_data.pressure_interpolated)[global_index] = p_L[n];
+        }
+
+        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.computeShapeFunctions(
+                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);
+            (*_process_data.pressure_interpolated)[global_index] = N_p * p_L;
+        }
+    };
+    interpolate_p();
 }
 
 }  // namespace RichardsMechanics
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp b/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
index d1ef1a3cd66cb61a4845b9773352bd416071d9a1..2765c51b04dc5bf4f5cac2ebbb8c5cda7b1186ae 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
@@ -195,6 +195,11 @@ void RichardsMechanicsProcess<DisplacementDim>::initializeConcreteProcess(
     _process_data.element_saturation = MeshLib::getOrCreateMeshProperty<double>(
         const_cast<MeshLib::Mesh&>(mesh), "saturation_avg",
         MeshLib::MeshItemType::Cell, 1);
+
+    _process_data.pressure_interpolated =
+        MeshLib::getOrCreateMeshProperty<double>(
+            const_cast<MeshLib::Mesh&>(mesh), "pressure_interpolated",
+            MeshLib::MeshItemType::Node, 1);
 }
 
 template <int DisplacementDim>
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h b/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h
index f737f7ce28ada2b4bbc29dcfc7e26832ddf99a16..c788c645b81f205237cf31d402f6c21a74248b87 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h
@@ -104,6 +104,7 @@ struct RichardsMechanicsProcessData
     double t = 0.0;
 
     MeshLib::PropertyVector<double>* element_saturation = nullptr;
+    MeshLib::PropertyVector<double>* pressure_interpolated = nullptr;
 
     EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
 };