diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/DrySolidDensity.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/DrySolidDensity.h
index a2af87eae9fbcad0d7102ab740b895387b279ca4..c021158782ccab86fc998c32233277b8532dd965 100644
--- a/ProcessLib/RichardsMechanics/ConstitutiveRelations/DrySolidDensity.h
+++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/DrySolidDensity.h
@@ -10,10 +10,17 @@
 
 #pragma once
 
+#include <string_view>
+
 #include "BaseLib/StrongType.h"
 
 namespace ProcessLib::RichardsMechanics
 {
 // Apparent dry solid density
 using DrySolidDensity = BaseLib::StrongType<double, struct DrySolidDensityTag>;
+
+constexpr std::string_view ioName(struct DrySolidDensityTag*)
+{
+    return "dry_density_solid";
+}
 }  // namespace ProcessLib::RichardsMechanics
diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/MicroPressure.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/MicroPressure.h
index 541be0780caef784f9bbf52f6db76e55bddaaa69..96ace800113e35d144b4299240cd281868268ef9 100644
--- a/ProcessLib/RichardsMechanics/ConstitutiveRelations/MicroPressure.h
+++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/MicroPressure.h
@@ -10,9 +10,16 @@
 
 #pragma once
 
+#include <string_view>
+
 #include "BaseLib/StrongType.h"
 
 namespace ProcessLib::RichardsMechanics
 {
 using MicroPressure = BaseLib::StrongType<double, struct MicroPressureTag>;
+
+constexpr std::string_view ioName(struct MicroPressureTag*)
+{
+    return "micro_pressure";
+}
 }  // namespace ProcessLib::RichardsMechanics
diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/MicroSaturation.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/MicroSaturation.h
index 2456caa8ed36239e66ca58cd8544b1617d589f6e..4048fdd58fec55d5ad277ea32914dfad44ada5ab 100644
--- a/ProcessLib/RichardsMechanics/ConstitutiveRelations/MicroSaturation.h
+++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/MicroSaturation.h
@@ -10,9 +10,16 @@
 
 #pragma once
 
+#include <string_view>
+
 #include "BaseLib/StrongType.h"
 
 namespace ProcessLib::RichardsMechanics
 {
 using MicroSaturation = BaseLib::StrongType<double, struct MicroSaturationTag>;
+
+constexpr std::string_view ioName(struct MicroSaturationTag*)
+{
+    return "micro_saturation";
+}
 }  // namespace ProcessLib::RichardsMechanics
diff --git a/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h b/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h
index 529002d1f651b7f59957c2d24c19c659abbe5ee6..a9a57f99fdfd4b7721c9d6b333ffb45544762d52 100644
--- a/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h
+++ b/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h
@@ -66,82 +66,6 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
         std::string_view const name, double const* values,
         int const integration_order) = 0;
 
-    virtual std::vector<double> getSigma() const = 0;
-
-    virtual std::vector<double> const& getIntPtSigma(
-        const double t,
-        std::vector<GlobalVector*> const& x,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
-        std::vector<double>& cache) const = 0;
-
-    virtual std::vector<double> getSwellingStress() const = 0;
-
-    virtual std::vector<double> const& getIntPtSwellingStress(
-        const double t,
-        std::vector<GlobalVector*> const& x,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
-        std::vector<double>& cache) const = 0;
-
-    virtual std::vector<double> getEpsilon() const = 0;
-
-    virtual std::vector<double> const& getIntPtEpsilon(
-        const double t,
-        std::vector<GlobalVector*> const& x,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
-        std::vector<double>& cache) const = 0;
-
-    virtual std::vector<double> const& getIntPtDarcyVelocity(
-        const double t,
-        std::vector<GlobalVector*> const& x,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
-        std::vector<double>& cache) const = 0;
-
-    virtual std::vector<double> getSaturation() const = 0;
-
-    virtual std::vector<double> const& getIntPtSaturation(
-        const double t,
-        std::vector<GlobalVector*> const& x,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
-        std::vector<double>& cache) const = 0;
-
-    virtual std::vector<double> getMicroSaturation() const = 0;
-
-    virtual std::vector<double> const& getIntPtMicroSaturation(
-        const double t,
-        std::vector<GlobalVector*> const& x,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
-        std::vector<double>& cache) const = 0;
-
-    virtual std::vector<double> getMicroPressure() const = 0;
-
-    virtual std::vector<double> const& getIntPtMicroPressure(
-        const double t,
-        std::vector<GlobalVector*> const& x,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
-        std::vector<double>& cache) const = 0;
-
-    virtual std::vector<double> getPorosity() const = 0;
-
-    virtual std::vector<double> const& getIntPtPorosity(
-        const double t,
-        std::vector<GlobalVector*> const& x,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
-        std::vector<double>& cache) const = 0;
-
-    virtual std::vector<double> getTransportPorosity() const = 0;
-
-    virtual std::vector<double> const& getIntPtTransportPorosity(
-        const double t,
-        std::vector<GlobalVector*> const& x,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
-        std::vector<double>& cache) const = 0;
-
-    virtual std::vector<double> const& getIntPtDryDensitySolid(
-        const double t,
-        std::vector<GlobalVector*> const& x,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
-        std::vector<double>& cache) const = 0;
-
     virtual std::vector<double> getMaterialStateVariableInternalState(
         std::function<std::span<double>(
             typename MaterialLib::Solids::MechanicsBase<DisplacementDim>::
@@ -157,6 +81,14 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
         DisplacementDim>::MaterialStateVariables const&
     getMaterialStateVariablesAt(unsigned /*integration_point*/) const = 0;
 
+    static auto getReflectionDataForOutput()
+    {
+        using Self = LocalAssemblerInterface<DisplacementDim>;
+
+        return ProcessLib::Reflection::reflectWithoutName(
+            &Self::current_states_, &Self::output_data_);
+    }
+
 protected:
     RichardsMechanicsProcessData<DisplacementDim>& process_data_;
     NumLib::GenericIntegrationMethod const& integration_method_;
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
index 4192ce317b52edc3d1a101893dc47d79bbffcdd5..f1bd168a25d130a85c5a59b129ccca10b0e0ca47 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
@@ -1557,121 +1557,6 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
         .noalias() += Kup * p_L;
 }
 
-template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
-          int DisplacementDim>
-std::vector<double> RichardsMechanicsLocalAssembler<
-    ShapeFunctionDisplacement, ShapeFunctionPressure,
-    DisplacementDim>::getSigma() const
-{
-    constexpr int kelvin_vector_size =
-        MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
-
-    return transposeInPlace<kelvin_vector_size>(
-        [this](std::vector<double>& values)
-        { return getIntPtSigma(0, {}, {}, values); });
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
-          int DisplacementDim>
-std::vector<double> const& RichardsMechanicsLocalAssembler<
-    ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
-    getIntPtSigma(
-        const double /*t*/,
-        std::vector<GlobalVector*> const& /*x*/,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
-        std::vector<double>& cache) const
-{
-    return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
-        this->current_states_,
-        [](auto& tuple) -> auto& {
-            return std::get<ProcessLib::ThermoRichardsMechanics::
-                                ConstitutiveStress_StrainTemperature::
-                                    EffectiveStressData<DisplacementDim>>(tuple)
-                .sigma_eff;
-        },
-        cache);
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
-          int DisplacementDim>
-std::vector<double> RichardsMechanicsLocalAssembler<
-    ShapeFunctionDisplacement, ShapeFunctionPressure,
-    DisplacementDim>::getSwellingStress() const
-{
-    constexpr int kelvin_vector_size =
-        MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
-
-    return transposeInPlace<kelvin_vector_size>(
-        [this](std::vector<double>& values)
-        { return getIntPtSwellingStress(0, {}, {}, values); });
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
-          int DisplacementDim>
-std::vector<double> const& RichardsMechanicsLocalAssembler<
-    ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
-    getIntPtSwellingStress(
-        const double /*t*/,
-        std::vector<GlobalVector*> const& /*x*/,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
-        std::vector<double>& cache) const
-{
-    constexpr int kelvin_vector_size =
-        MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
-    auto const n_integration_points = _ip_data.size();
-
-    cache.clear();
-    auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
-        double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
-        cache, kelvin_vector_size, n_integration_points);
-
-    for (unsigned ip = 0; ip < n_integration_points; ++ip)
-    {
-        auto const& sigma_sw =
-            std::get<ProcessLib::ThermoRichardsMechanics::
-                         ConstitutiveStress_StrainTemperature::
-                             SwellingDataStateful<DisplacementDim>>(
-                this->current_states_[ip])
-                .sigma_sw;
-        cache_mat.col(ip) =
-            MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_sw);
-    }
-
-    return cache;
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
-          int DisplacementDim>
-std::vector<double> RichardsMechanicsLocalAssembler<
-    ShapeFunctionDisplacement, ShapeFunctionPressure,
-    DisplacementDim>::getEpsilon() const
-{
-    constexpr int kelvin_vector_size =
-        MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
-
-    return transposeInPlace<kelvin_vector_size>(
-        [this](std::vector<double>& values)
-        { return getIntPtEpsilon(0, {}, {}, values); });
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
-          int DisplacementDim>
-std::vector<double> const& RichardsMechanicsLocalAssembler<
-    ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
-    getIntPtEpsilon(
-        const double /*t*/,
-        std::vector<GlobalVector*> const& /*x*/,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
-        std::vector<double>& cache) const
-{
-    return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
-        this->current_states_,
-        [](auto& tuple) -> auto& {
-            return std::get<StrainData<DisplacementDim>>(tuple).eps;
-        },
-        cache);
-}
-
 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
           int DisplacementDim>
 int RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
@@ -1698,202 +1583,6 @@ std::vector<double> RichardsMechanicsLocalAssembler<
         n_components);
 }
 
-template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
-          int DisplacementDim>
-std::vector<double> const& RichardsMechanicsLocalAssembler<
-    ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
-    getIntPtDarcyVelocity(
-        const double /*t*/,
-        std::vector<GlobalVector*> const& /*x*/,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
-        std::vector<double>& cache) const
-{
-    unsigned const n_integration_points =
-        this->integration_method_.getNumberOfPoints();
-
-    cache.clear();
-    auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
-        double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
-        cache, DisplacementDim, n_integration_points);
-
-    for (unsigned ip = 0; ip < n_integration_points; ip++)
-    {
-        cache_matrix.col(ip).noalias() = *std::get<
-            ProcessLib::ThermoRichardsMechanics::DarcyLawData<DisplacementDim>>(
-            this->output_data_[ip]);
-    }
-
-    return cache;
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
-          int DisplacementDim>
-std::vector<double> RichardsMechanicsLocalAssembler<
-    ShapeFunctionDisplacement, ShapeFunctionPressure,
-    DisplacementDim>::getSaturation() const
-{
-    std::vector<double> result;
-    getIntPtSaturation(0, {}, {}, result);
-    return result;
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
-          int DisplacementDim>
-std::vector<double> const& RichardsMechanicsLocalAssembler<
-    ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
-    getIntPtSaturation(
-        const double /*t*/,
-        std::vector<GlobalVector*> const& /*x*/,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
-        std::vector<double>& cache) const
-{
-    return ProcessLib::getIntegrationPointScalarData(
-        this->current_states_,
-        [](auto& tuple) -> auto&
-        {
-            return std::get<
-                       ProcessLib::ThermoRichardsMechanics::SaturationData>(
-                       tuple)
-                .S_L;
-        },
-        cache);
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
-          int DisplacementDim>
-std::vector<double> RichardsMechanicsLocalAssembler<
-    ShapeFunctionDisplacement, ShapeFunctionPressure,
-    DisplacementDim>::getMicroSaturation() const
-{
-    std::vector<double> result;
-    getIntPtMicroSaturation(0, {}, {}, result);
-    return result;
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
-          int DisplacementDim>
-std::vector<double> const& RichardsMechanicsLocalAssembler<
-    ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
-    getIntPtMicroSaturation(
-        const double /*t*/,
-        std::vector<GlobalVector*> const& /*x*/,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
-        std::vector<double>& cache) const
-{
-    return ProcessLib::getIntegrationPointScalarData(
-        this->current_states_,
-        [](auto& tuple) -> auto& { return *std::get<MicroSaturation>(tuple); },
-        cache);
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
-          int DisplacementDim>
-std::vector<double> RichardsMechanicsLocalAssembler<
-    ShapeFunctionDisplacement, ShapeFunctionPressure,
-    DisplacementDim>::getMicroPressure() const
-{
-    std::vector<double> result;
-    getIntPtMicroPressure(0, {}, {}, result);
-    return result;
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
-          int DisplacementDim>
-std::vector<double> const& RichardsMechanicsLocalAssembler<
-    ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
-    getIntPtMicroPressure(
-        const double /*t*/,
-        std::vector<GlobalVector*> const& /*x*/,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
-        std::vector<double>& cache) const
-{
-    return ProcessLib::getIntegrationPointScalarData(
-        this->current_states_,
-        [](auto& tuple) -> auto& { return *std::get<MicroPressure>(tuple); },
-        cache);
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
-          int DisplacementDim>
-std::vector<double> RichardsMechanicsLocalAssembler<
-    ShapeFunctionDisplacement, ShapeFunctionPressure,
-    DisplacementDim>::getPorosity() const
-{
-    std::vector<double> result;
-    getIntPtPorosity(0, {}, {}, result);
-    return result;
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
-          int DisplacementDim>
-std::vector<double> const& RichardsMechanicsLocalAssembler<
-    ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
-    getIntPtPorosity(
-        const double /*t*/,
-        std::vector<GlobalVector*> const& /*x*/,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
-        std::vector<double>& cache) const
-{
-    return ProcessLib::getIntegrationPointScalarData(
-        this->current_states_,
-        [](auto& tuple) -> auto&
-        {
-            return std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(
-                       tuple)
-                .phi;
-        },
-        cache);
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
-          int DisplacementDim>
-std::vector<double> RichardsMechanicsLocalAssembler<
-    ShapeFunctionDisplacement, ShapeFunctionPressure,
-    DisplacementDim>::getTransportPorosity() const
-{
-    std::vector<double> result;
-    getIntPtTransportPorosity(0, {}, {}, result);
-    return result;
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
-          int DisplacementDim>
-std::vector<double> const& RichardsMechanicsLocalAssembler<
-    ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
-    getIntPtTransportPorosity(
-        const double /*t*/,
-        std::vector<GlobalVector*> const& /*x*/,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
-        std::vector<double>& cache) const
-{
-    return ProcessLib::getIntegrationPointScalarData(
-        this->current_states_,
-        [](auto& tuple) -> auto&
-        {
-            return std::
-                get<ProcessLib::ThermoRichardsMechanics::TransportPorosityData>(
-                       tuple)
-                    .phi;
-        },
-        cache);
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
-          int DisplacementDim>
-std::vector<double> const& RichardsMechanicsLocalAssembler<
-    ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
-    getIntPtDryDensitySolid(
-        const double /*t*/,
-        std::vector<GlobalVector*> const& /*x*/,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
-        std::vector<double>& cache) const
-{
-    return ProcessLib::getIntegrationPointScalarData(
-        this->output_data_,
-        [](auto& tuple) -> auto& { return *std::get<DrySolidDensity>(tuple); },
-        cache);
-}
-
 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
           int DisplacementDim>
 void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
index d18799cba7d5b32213d8b36392c680bbad150e1f..e00e0c864c9cbecaa7b4ffe16f4c3124be256b60 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
@@ -196,69 +196,6 @@ public:
         return Eigen::Map<const Eigen::RowVectorXd>(N_u.data(), N_u.size());
     }
 
-    std::vector<double> getSigma() const override;
-
-    std::vector<double> const& getIntPtDarcyVelocity(
-        const double t,
-        std::vector<GlobalVector*> const& x,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
-        std::vector<double>& cache) const override;
-
-    std::vector<double> getSaturation() const override;
-    std::vector<double> const& getIntPtSaturation(
-        const double t,
-        std::vector<GlobalVector*> const& x,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
-        std::vector<double>& cache) const override;
-
-    std::vector<double> getMicroSaturation() const override;
-    std::vector<double> const& getIntPtMicroSaturation(
-        const double t,
-        std::vector<GlobalVector*> const& x,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
-        std::vector<double>& cache) const override;
-
-    std::vector<double> getMicroPressure() const override;
-    std::vector<double> const& getIntPtMicroPressure(
-        const double t,
-        std::vector<GlobalVector*> const& x,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
-        std::vector<double>& cache) const override;
-
-    std::vector<double> getPorosity() const override;
-    std::vector<double> const& getIntPtPorosity(
-        const double t,
-        std::vector<GlobalVector*> const& x,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
-        std::vector<double>& cache) const override;
-
-    std::vector<double> getTransportPorosity() const override;
-    std::vector<double> const& getIntPtTransportPorosity(
-        const double t,
-        std::vector<GlobalVector*> const& x,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
-        std::vector<double>& cache) const override;
-
-    std::vector<double> const& getIntPtSigma(
-        const double t,
-        std::vector<GlobalVector*> const& x,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
-        std::vector<double>& cache) const override;
-
-    std::vector<double> getSwellingStress() const override;
-    std::vector<double> const& getIntPtSwellingStress(
-        const double t,
-        std::vector<GlobalVector*> const& x,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
-        std::vector<double>& cache) const override;
-
-    std::vector<double> getEpsilon() const override;
-    std::vector<double> const& getIntPtEpsilon(
-        const double t,
-        std::vector<GlobalVector*> const& x,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
-        std::vector<double>& cache) const override;
-
     int getMaterialID() const override;
 
     std::vector<double> getMaterialStateVariableInternalState(
@@ -267,12 +204,6 @@ public:
                 MaterialStateVariables&)> const& get_values_span,
         int const& n_components) const override;
 
-    std::vector<double> const& getIntPtDryDensitySolid(
-        const double t,
-        std::vector<GlobalVector*> const& x,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
-        std::vector<double>& cache) const override;
-
 private:
     /**
      * Assemble local matrices and vectors arise from the linearized discretized
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp b/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
index ff36291ec52ea932d86ffe883b2a28fc29fc897c..4ab85ce74c2511533906b1711a132e2d080fb9a9 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
@@ -16,6 +16,8 @@
 #include "MeshLib/Utils/getOrCreateMeshProperty.h"
 #include "NumLib/DOF/ComputeSparsityPattern.h"
 #include "ProcessLib/Deformation/SolidMaterialInternalToSecondaryVariables.h"
+#include "ProcessLib/Reflection/ReflectionForExtrapolation.h"
+#include "ProcessLib/Reflection/ReflectionForIPWriters.h"
 #include "ProcessLib/Utils/CreateLocalAssemblersTaylorHood.h"
 #include "ProcessLib/Utils/SetIPDataInitialConditions.h"
 #include "RichardsMechanicsFEM.h"
@@ -48,43 +50,10 @@ RichardsMechanicsProcess<DisplacementDim>::RichardsMechanicsProcess(
     _hydraulic_flow = MeshLib::getOrCreateMeshProperty<double>(
         mesh, "MassFlowRate", MeshLib::MeshItemType::Node, 1);
 
-    // TODO (naumov) remove ip suffix. Probably needs modification of the mesh
-    // properties, s.t. there is no "overlapping" with cell/point data.
-    // See getOrCreateMeshProperty.
-    _integration_point_writer.emplace_back(
-        std::make_unique<MeshLib::IntegrationPointWriter>(
-            "sigma_ip",
-            static_cast<int>(mesh.getDimension() == 2 ? 4 : 6) /*n components*/,
-            integration_order, _local_assemblers, &LocalAssemblerIF::getSigma));
-
-    _integration_point_writer.emplace_back(
-        std::make_unique<MeshLib::IntegrationPointWriter>(
-            "saturation_ip", 1 /*n components*/, integration_order,
-            _local_assemblers, &LocalAssemblerIF::getSaturation));
-
-    _integration_point_writer.emplace_back(
-        std::make_unique<MeshLib::IntegrationPointWriter>(
-            "porosity_ip", 1 /*n components*/, integration_order,
-            _local_assemblers, &LocalAssemblerIF::getPorosity));
-
-    _integration_point_writer.emplace_back(
-        std::make_unique<MeshLib::IntegrationPointWriter>(
-            "transport_porosity_ip", 1 /*n components*/, integration_order,
-            _local_assemblers, &LocalAssemblerIF::getTransportPorosity));
-
-    _integration_point_writer.emplace_back(
-        std::make_unique<MeshLib::IntegrationPointWriter>(
-            "swelling_stress_ip",
-            static_cast<int>(mesh.getDimension() == 2 ? 4 : 6) /*n components*/,
-            integration_order, _local_assemblers,
-            &LocalAssemblerIF::getSwellingStress));
-
-    _integration_point_writer.emplace_back(
-        std::make_unique<MeshLib::IntegrationPointWriter>(
-            "epsilon_ip",
-            static_cast<int>(mesh.getDimension() == 2 ? 4 : 6) /*n components*/,
-            integration_order, _local_assemblers,
-            &LocalAssemblerIF::getEpsilon));
+    ProcessLib::Reflection::addReflectedIntegrationPointWriters<
+        DisplacementDim>(LocalAssemblerIF::getReflectionDataForOutput(),
+                         _integration_point_writer, integration_order,
+                         _local_assemblers);
 }
 
 template <int DisplacementDim>
@@ -202,6 +171,10 @@ void RichardsMechanicsProcess<DisplacementDim>::initializeConcreteProcess(
         NumLib::IntegrationOrder{integration_order}, mesh.isAxiallySymmetric(),
         _process_data);
 
+    ProcessLib::Reflection::addReflectedSecondaryVariables<DisplacementDim>(
+        LocalAssemblerIF::getReflectionDataForOutput(), _secondary_variables,
+        getExtrapolator(), _local_assemblers);
+
     auto add_secondary_variable = [&](std::string const& name,
                                       int const num_components,
                                       auto get_ip_values_function)
@@ -213,41 +186,6 @@ void RichardsMechanicsProcess<DisplacementDim>::initializeConcreteProcess(
                              std::move(get_ip_values_function)));
     };
 
-    add_secondary_variable("sigma",
-                           MathLib::KelvinVector::KelvinVectorType<
-                               DisplacementDim>::RowsAtCompileTime,
-                           &LocalAssemblerIF::getIntPtSigma);
-
-    add_secondary_variable("swelling_stress",
-                           MathLib::KelvinVector::KelvinVectorType<
-                               DisplacementDim>::RowsAtCompileTime,
-                           &LocalAssemblerIF::getIntPtSwellingStress);
-
-    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);
-
-    add_secondary_variable("micro_saturation", 1,
-                           &LocalAssemblerIF::getIntPtMicroSaturation);
-
-    add_secondary_variable("micro_pressure", 1,
-                           &LocalAssemblerIF::getIntPtMicroPressure);
-
-    add_secondary_variable("porosity", 1, &LocalAssemblerIF::getIntPtPorosity);
-
-    add_secondary_variable("transport_porosity", 1,
-                           &LocalAssemblerIF::getIntPtTransportPorosity);
-
-    add_secondary_variable("dry_density_solid", 1,
-                           &LocalAssemblerIF::getIntPtDryDensitySolid);
-
     //
     // enable output of internal variables defined by material models
     //