diff --git a/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h b/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h
index eaf661fd32f27264a226909763c03c99bc7dd4b8..c4864173aa7d4e9f623a7b0a3da1fa3910561438 100644
--- a/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h
+++ b/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h
@@ -10,6 +10,7 @@
 
 #pragma once
 
+#include "MaterialLib/SolidModels/MechanicsBase.h"
 #include "NumLib/Extrapolation/ExtrapolatableElement.h"
 #include "ProcessLib/LocalAssemblerInterface.h"
 
@@ -17,6 +18,7 @@ namespace ProcessLib
 {
 namespace RichardsMechanics
 {
+template <int DisplacementDim>
 struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
                                  public NumLib::ExtrapolatableElement
 {
@@ -43,6 +45,13 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
         std::vector<GlobalVector*> const& x,
         std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
         std::vector<double>& cache) const = 0;
+
+    // TODO move to NumLib::ExtrapolatableElement
+    virtual unsigned getNumberOfIntegrationPoints() const = 0;
+
+    virtual typename MaterialLib::Solids::MechanicsBase<
+        DisplacementDim>::MaterialStateVariables const&
+    getMaterialStateVariablesAt(unsigned /*integration_point*/) const = 0;
 };
 
 }  // namespace RichardsMechanics
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
index c3d7b3ecf5ec20967a65e083b343efa33dc2606a..d4450f4f295d16d82dca4e385f7ca26ef6c54f35 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
@@ -886,5 +886,25 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
                          *_process_data.pressure_interpolated);
 }
 
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+          typename IntegrationMethod, int DisplacementDim>
+unsigned RichardsMechanicsLocalAssembler<
+    ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
+    DisplacementDim>::getNumberOfIntegrationPoints() const
+{
+    return _integration_method.getNumberOfPoints();
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+          typename IntegrationMethod, int DisplacementDim>
+typename MaterialLib::Solids::MechanicsBase<
+    DisplacementDim>::MaterialStateVariables const&
+RichardsMechanicsLocalAssembler<
+    ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
+    DisplacementDim>::getMaterialStateVariablesAt(unsigned integration_point)
+    const
+{
+    return *_ip_data[integration_point].material_state_variables;
+}
 }  // namespace RichardsMechanics
 }  // namespace ProcessLib
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
index 45db4b3697186b755a912d5a47f4b156832c5b78..89ce9cfa5068b6afe727e64e5d9d17e395712217 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
@@ -42,7 +42,8 @@ struct SecondaryData
 
 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
           typename IntegrationMethod, int DisplacementDim>
-class RichardsMechanicsLocalAssembler : public LocalAssemblerInterface
+class RichardsMechanicsLocalAssembler
+    : public LocalAssemblerInterface<DisplacementDim>
 {
 public:
     using ShapeMatricesTypeDisplacement =
@@ -244,6 +245,12 @@ private:
         std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
         LocalCoupledSolutions const& local_coupled_solutions);
 
+    unsigned getNumberOfIntegrationPoints() const override;
+
+    typename MaterialLib::Solids::MechanicsBase<
+        DisplacementDim>::MaterialStateVariables const&
+    getMaterialStateVariablesAt(unsigned integration_point) const override;
+
 private:
     RichardsMechanicsProcessData<DisplacementDim>& _process_data;
 
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp b/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
index 33d9b5debec87dda9d1ba1830507ea5678006b1e..fd5ba4b532579916ef426a52254a19c634294e05 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
@@ -14,9 +14,9 @@
 
 #include "MeshLib/Elements/Utils.h"
 #include "NumLib/DOF/ComputeSparsityPattern.h"
+#include "ProcessLib/Deformation/SolidMaterialInternalToSecondaryVariables.h"
 #include "ProcessLib/Process.h"
 #include "ProcessLib/RichardsMechanics/CreateLocalAssemblers.h"
-
 #include "RichardsMechanicsFEM.h"
 #include "RichardsMechanicsProcessData.h"
 
@@ -168,29 +168,39 @@ void RichardsMechanicsProcess<DisplacementDim>::initializeConcreteProcess(
         _local_assemblers, mesh.isAxiallySymmetric(), integration_order,
         _process_data);
 
-    _secondary_variables.addSecondaryVariable(
-        "sigma",
-        makeExtrapolator(MathLib::KelvinVector::KelvinVectorType<
-                             DisplacementDim>::RowsAtCompileTime,
-                         getExtrapolator(), _local_assemblers,
-                         &LocalAssemblerInterface::getIntPtSigma));
-
-    _secondary_variables.addSecondaryVariable(
-        "epsilon",
-        makeExtrapolator(MathLib::KelvinVector::KelvinVectorType<
-                             DisplacementDim>::RowsAtCompileTime,
-                         getExtrapolator(), _local_assemblers,
-                         &LocalAssemblerInterface::getIntPtEpsilon));
-
-    _secondary_variables.addSecondaryVariable(
-        "velocity",
-        makeExtrapolator(DisplacementDim, getExtrapolator(), _local_assemblers,
-                         &LocalAssemblerInterface::getIntPtDarcyVelocity));
-
-    _secondary_variables.addSecondaryVariable(
-        "saturation",
-        makeExtrapolator(1, getExtrapolator(), _local_assemblers,
-                         &LocalAssemblerInterface::getIntPtSaturation));
+    auto add_secondary_variable = [&](std::string const& name,
+                                      int const num_components,
+                                      auto get_ip_values_function) {
+        _secondary_variables.addSecondaryVariable(
+            name,
+            makeExtrapolator(num_components, getExtrapolator(),
+                             _local_assemblers,
+                             std::move(get_ip_values_function)));
+    };
+
+    add_secondary_variable("sigma",
+                           MathLib::KelvinVector::KelvinVectorType<
+                               DisplacementDim>::RowsAtCompileTime,
+                           &LocalAssemblerIF::getIntPtSigma);
+
+    add_secondary_variable("epsilon",
+                           MathLib::KelvinVector::KelvinVectorType<
+                               DisplacementDim>::RowsAtCompileTime,
+                           &LocalAssemblerIF::getIntPtEpsilon);
+
+    add_secondary_variable("velocity",
+                           DisplacementDim,
+                           &LocalAssemblerIF::getIntPtDarcyVelocity);
+
+    add_secondary_variable("saturation", 1,
+                           &LocalAssemblerIF::getIntPtSaturation);
+
+    //
+    // enable output of internal variables defined by material models
+    //
+    ProcessLib::Deformation::solidMaterialInternalToSecondaryVariables<
+        LocalAssemblerIF>(_process_data.solid_materials,
+                                 add_secondary_variable);
 
     _process_data.element_saturation = MeshLib::getOrCreateMeshProperty<double>(
         const_cast<MeshLib::Mesh&>(mesh), "saturation_avg",
@@ -203,7 +213,7 @@ void RichardsMechanicsProcess<DisplacementDim>::initializeConcreteProcess(
 
     // Initialize local assemblers after all variables have been set.
     GlobalExecutor::executeMemberOnDereferenced(
-        &LocalAssemblerInterface::initialize, _local_assemblers,
+        &LocalAssemblerIF::initialize, _local_assemblers,
         *_local_to_global_index_map);
 }
 
@@ -332,7 +342,7 @@ void RichardsMechanicsProcess<DisplacementDim>::preTimestepConcreteProcess(
         ProcessLib::ProcessVariable const& pv =
             getProcessVariables(process_id)[0];
         GlobalExecutor::executeSelectedMemberOnDereferenced(
-            &LocalAssemblerInterface::preTimestep, _local_assemblers,
+            &LocalAssemblerIF::preTimestep, _local_assemblers,
             pv.getActiveElementIDs(), *_local_to_global_index_map,
             *x[process_id], t, dt);
     }
@@ -355,7 +365,7 @@ void RichardsMechanicsProcess<
 
     // Calculate strain, stress or other internal variables of mechanics.
     GlobalExecutor::executeSelectedMemberOnDereferenced(
-        &LocalAssemblerInterface::postNonLinearSolver, _local_assemblers,
+        &LocalAssemblerIF::postNonLinearSolver, _local_assemblers,
         pv.getActiveElementIDs(), getDOFTable(process_id), x, t, dt,
         _use_monolithic_scheme);
 }
@@ -370,7 +380,7 @@ void RichardsMechanicsProcess<
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     GlobalExecutor::executeSelectedMemberOnDereferenced(
-        &LocalAssemblerInterface::computeSecondaryVariable, _local_assemblers,
+        &LocalAssemblerIF::computeSecondaryVariable, _local_assemblers,
         pv.getActiveElementIDs(), getDOFTable(process_id), t, x,
         _coupled_solutions);
 }
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.h b/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.h
index f14bb2bec8acf27adb5754350edf435e004204d4..c50b236dc46ad951c7ff5a251d9977681539eccb 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.h
@@ -18,8 +18,6 @@ namespace ProcessLib
 {
 namespace RichardsMechanics
 {
-struct LocalAssemblerInterface;
-
 /// Linear kinematics poro-mechanical/biphasic (fluid-solid mixture) model.
 ///
 /// The mixture momentum balance and the mixture mass balance are solved under
@@ -64,6 +62,8 @@ public:
         const int process_id) const override;
 
 private:
+    using LocalAssemblerIF = LocalAssemblerInterface<DisplacementDim>;
+
     void constructDofTable() override;
 
     void initializeConcreteProcess(
@@ -100,7 +100,7 @@ private:
     std::unique_ptr<MeshLib::MeshSubset const> _mesh_subset_base_nodes;
     RichardsMechanicsProcessData<DisplacementDim> _process_data;
 
-    std::vector<std::unique_ptr<LocalAssemblerInterface>> _local_assemblers;
+    std::vector<std::unique_ptr<LocalAssemblerIF>> _local_assemblers;
 
     std::unique_ptr<NumLib::LocalToGlobalIndexMap>
         _local_to_global_index_map_single_component;