diff --git a/ProcessLib/LargeDeformation/LargeDeformationProcess.cpp b/ProcessLib/LargeDeformation/LargeDeformationProcess.cpp
index f19707a2dc24ff82466d0b01ece8fdb1d375b24a..724f100c1ce69c80879648152d56cd458516b711 100644
--- a/ProcessLib/LargeDeformation/LargeDeformationProcess.cpp
+++ b/ProcessLib/LargeDeformation/LargeDeformationProcess.cpp
@@ -18,6 +18,7 @@
 #include "MeshLib/Utils/getOrCreateMeshProperty.h"
 #include "NumLib/DOF/DOFTableUtil.h"
 #include "ProcessLib/Deformation/SolidMaterialInternalToSecondaryVariables.h"
+#include "ProcessLib/Output/CellAverageAlgorithm.h"
 #include "ProcessLib/Process.h"
 #include "ProcessLib/Reflection/ReflectionForExtrapolation.h"
 #include "ProcessLib/Reflection/ReflectionForIPWriters.h"
@@ -207,6 +208,8 @@ void LargeDeformationProcess<DisplacementDim>::computeSecondaryVariableConcrete(
     GlobalExecutor::executeSelectedMemberOnDereferenced(
         &LocalAssemblerInterface::computeSecondaryVariable, _local_assemblers,
         pv.getActiveElementIDs(), dof_tables, t, dt, x, x_prev, process_id);
+
+    computeCellAverages<DisplacementDim>(cell_average_data_, _local_assemblers);
 }
 template class LargeDeformationProcess<2>;
 template class LargeDeformationProcess<3>;
diff --git a/ProcessLib/Output/CellAverageAlgorithm.h b/ProcessLib/Output/CellAverageAlgorithm.h
new file mode 100644
index 0000000000000000000000000000000000000000..1bb82cc45f4dda0fce7957c924f2ed4292ce7f01
--- /dev/null
+++ b/ProcessLib/Output/CellAverageAlgorithm.h
@@ -0,0 +1,63 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "CellAverageData.h"
+#include "ProcessLib/Reflection/ReflectionIPData.h"
+
+namespace ProcessLib
+{
+namespace detail
+{
+void computeCellAverages(CellAverageData& cell_average_data,
+                         std::string const& name, unsigned const num_comp,
+                         auto&& flattened_ip_data_accessor,
+                         auto const& local_assemblers)
+{
+    auto& prop_vec =
+        cell_average_data.getOrCreatePropertyVector(name, num_comp);
+
+    for (std::size_t i = 0; i < local_assemblers.size(); ++i)
+    {
+        auto const& loc_asm = *local_assemblers[i];
+        auto const& ip_data = flattened_ip_data_accessor(loc_asm);
+        assert(ip_data.size() % num_comp == 0);
+        auto const num_ips =
+            static_cast<Eigen::Index>(ip_data.size() / num_comp);
+        Eigen::Map<const Eigen::MatrixXd> ip_data_mapped{ip_data.data(),
+                                                         num_comp, num_ips};
+
+        Eigen::Map<Eigen::VectorXd>{&prop_vec[i * num_comp], num_comp} =
+            ip_data_mapped.rowwise().mean();
+    }
+}
+}  // namespace detail
+
+template <int dim, typename LAIntf>
+void computeCellAverages(
+    CellAverageData& cell_average_data,
+    std::vector<std::unique_ptr<LAIntf>> const& local_assemblers)
+{
+    auto const callback = [&cell_average_data, &local_assemblers](
+                              std::string const& name,
+                              unsigned const num_comp,
+                              auto&& flattened_ip_data_accessor)
+    {
+        detail::computeCellAverages(cell_average_data, name, num_comp,
+                                    flattened_ip_data_accessor,
+                                    local_assemblers);
+    };
+
+    ProcessLib::Reflection::forEachReflectedFlattenedIPDataAccessor<dim,
+                                                                    LAIntf>(
+        LAIntf::getReflectionDataForOutput(), callback);
+}
+}  // namespace ProcessLib
diff --git a/ProcessLib/Output/CellAverageData.cpp b/ProcessLib/Output/CellAverageData.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1e46fa1a4cf53f272016dbb0475da1dd29c8a7b5
--- /dev/null
+++ b/ProcessLib/Output/CellAverageData.cpp
@@ -0,0 +1,59 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "CellAverageData.h"
+
+#include "MeshLib/Utils/getOrCreateMeshProperty.h"
+
+namespace ProcessLib
+{
+MeshLib::PropertyVector<double>& CellAverageData::getOrCreatePropertyVector(
+    std::string const& name, unsigned const num_comp)
+{
+    if (auto const it = cell_averages_.find(name); it != cell_averages_.end())
+    {
+        auto& prop_vec = *it->second;
+        auto const num_comp_mesh = prop_vec.getNumberOfGlobalComponents();
+        if (num_comp_mesh == static_cast<int>(num_comp))
+        {
+            return prop_vec;
+        }
+
+        OGS_FATAL(
+            "The requested property '{}' has {} components, but the one "
+            "present in the mesh has {} components.",
+            name, num_comp, num_comp_mesh);
+    }
+
+    auto const name_in_mesh = name + "_avg";
+    auto [it, emplaced] = cell_averages_.emplace(
+        name, MeshLib::getOrCreateMeshProperty<double>(
+                  const_cast<MeshLib::Mesh&>(mesh_), name_in_mesh,
+                  MeshLib::MeshItemType::Cell, num_comp));
+
+    if (!it->second)
+    {
+        OGS_FATAL("The cell property '{}' could not be added to the mesh.",
+                  name_in_mesh);
+    }
+
+    if (!emplaced)
+    {
+        OGS_FATAL(
+            "Internal logic error. Something very bad happened. The cell "
+            "property '{}' was not added to the list of cell averages to "
+            "compute. There is some very strange inconsistency in the "
+            "code. Trouble ahead!",
+            name_in_mesh);
+    }
+
+    return *it->second;
+}
+}  // namespace ProcessLib
diff --git a/ProcessLib/Output/CellAverageData.h b/ProcessLib/Output/CellAverageData.h
new file mode 100644
index 0000000000000000000000000000000000000000..19dda7699d12a75daadac234b60f7a411355ec2d
--- /dev/null
+++ b/ProcessLib/Output/CellAverageData.h
@@ -0,0 +1,29 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "MeshLib/Mesh.h"
+#include "MeshLib/PropertyVector.h"
+
+namespace ProcessLib
+{
+struct CellAverageData
+{
+    explicit CellAverageData(MeshLib::Mesh& mesh) : mesh_{mesh} {}
+
+    MeshLib::PropertyVector<double>& getOrCreatePropertyVector(
+        std::string const& name, unsigned const num_comp);
+
+private:
+    MeshLib::Mesh const& mesh_;
+    std::map<std::string, MeshLib::PropertyVector<double>*> cell_averages_;
+};
+}  // namespace ProcessLib
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index 98ed6566750c5a9f1d7583f0eb1beaf29e4d7053..9e9277e0169239101727d654193997de34c085e7 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -54,6 +54,7 @@ Process::Process(
     : name(std::move(name_)),
       _mesh(mesh),
       _secondary_variables(std::move(secondary_variables)),
+      cell_average_data_(mesh),
       _jacobian_assembler(std::move(jacobian_assembler)),
       _global_assembler(*_jacobian_assembler),
       _use_monolithic_scheme(use_monolithic_scheme),
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index 1816b836668da186e2062fa64c059408b627a3a9..be39d2bcb1369664303c9423abd2250578302833 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -21,6 +21,7 @@
 #include "ParameterLib/Parameter.h"
 #include "ProcessLib/BoundaryConditionAndSourceTerm/BoundaryConditionCollection.h"
 #include "ProcessLib/BoundaryConditionAndSourceTerm/SourceTermCollection.h"
+#include "ProcessLib/Output/CellAverageData.h"
 #include "ProcessLib/Output/ExtrapolatorData.h"
 #include "ProcessLib/Output/SecondaryVariable.h"
 #include "ProcessVariable.h"
@@ -361,6 +362,8 @@ protected:
 
     SecondaryVariableCollection _secondary_variables;
 
+    CellAverageData cell_average_data_;
+
     // TODO (CL) switch to the parallel vector matrix assembler here, once
     // feature complete, or use the assembly mixin and remove these two members.
     std::unique_ptr<ProcessLib::AbstractJacobianAssembler> _jacobian_assembler;
diff --git a/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp b/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp
index 2b6223c552456699dc8eb2744e25c912cd45c6d9..d414735825172d99ef31c13481107b0e92b983af 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp
+++ b/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp
@@ -15,6 +15,7 @@
 #include "MeshLib/Utils/IntegrationPointWriter.h"
 #include "MeshLib/Utils/getOrCreateMeshProperty.h"
 #include "ProcessLib/Deformation/SolidMaterialInternalToSecondaryVariables.h"
+#include "ProcessLib/Output/CellAverageAlgorithm.h"
 #include "ProcessLib/Process.h"
 #include "ProcessLib/Reflection/ReflectionForExtrapolation.h"
 #include "ProcessLib/Reflection/ReflectionForIPWriters.h"
@@ -216,6 +217,8 @@ void SmallDeformationProcess<DisplacementDim>::computeSecondaryVariableConcrete(
     GlobalExecutor::executeSelectedMemberOnDereferenced(
         &LocalAssemblerInterface::computeSecondaryVariable, _local_assemblers,
         pv.getActiveElementIDs(), dof_tables, t, dt, x, x_prev, process_id);
+
+    computeCellAverages<DisplacementDim>(cell_average_data_, _local_assemblers);
 }
 template class SmallDeformationProcess<2>;
 template class SmallDeformationProcess<3>;
diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index 154248a716b36a6b3576d6aab655ab79d156e17d..42d127276f9a65aaa8cf626ac7393123d3c2988d 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -1910,23 +1910,10 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
                          temperature,
                          *this->process_data_.temperature_interpolated);
 
-    unsigned const n_integration_points =
-        this->integration_method_.getNumberOfPoints();
-
-    double saturation_avg = 0;
-
     ConstitutiveRelations::ConstitutiveModels<DisplacementDim> const models{
         this->solid_material_, *this->process_data_.phase_transition_model_};
 
     updateConstitutiveVariables(local_x, local_x_prev, t, dt, models);
-
-    for (unsigned ip = 0; ip < n_integration_points; ip++)
-    {
-        saturation_avg += this->current_states_[ip].S_L_data.S_L;
-    }
-    saturation_avg /= n_integration_points;
-    (*this->process_data_.element_saturation)[this->element_.getID()] =
-        saturation_avg;
 }
 
 }  // namespace TH2M
diff --git a/ProcessLib/TH2M/TH2MProcess.cpp b/ProcessLib/TH2M/TH2MProcess.cpp
index b3289ca9d146d002cf4841f474efc27f61ee5d8d..376369b107e41a2beb98a765fb25e77da72b1648 100644
--- a/ProcessLib/TH2M/TH2MProcess.cpp
+++ b/ProcessLib/TH2M/TH2MProcess.cpp
@@ -19,6 +19,7 @@
 #include "MeshLib/Utils/getOrCreateMeshProperty.h"
 #include "NumLib/DOF/DOFTableUtil.h"
 #include "ProcessLib/Deformation/SolidMaterialInternalToSecondaryVariables.h"
+#include "ProcessLib/Output/CellAverageAlgorithm.h"
 #include "ProcessLib/Process.h"
 #include "ProcessLib/Reflection/ReflectionForExtrapolation.h"
 #include "ProcessLib/Reflection/ReflectionForIPWriters.h"
@@ -175,10 +176,6 @@ void TH2MProcess<DisplacementDim>::initializeConcreteProcess(
             _process_data.solid_materials, local_assemblers_,
             _integration_point_writer, integration_order);
 
-    _process_data.element_saturation = MeshLib::getOrCreateMeshProperty<double>(
-        const_cast<MeshLib::Mesh&>(mesh), "saturation_avg",
-        MeshLib::MeshItemType::Cell, 1);
-
     _process_data.gas_pressure_interpolated =
         MeshLib::getOrCreateMeshProperty<double>(
             const_cast<MeshLib::Mesh&>(mesh), "gas_pressure_interpolated",
@@ -320,6 +317,8 @@ void TH2MProcess<DisplacementDim>::computeSecondaryVariableConcrete(
         &LocalAssemblerInterface<DisplacementDim>::computeSecondaryVariable,
         local_assemblers_, pv.getActiveElementIDs(), getDOFTables(x.size()), t,
         dt, x, x_prev, process_id);
+
+    computeCellAverages<DisplacementDim>(cell_average_data_, local_assemblers_);
 }
 
 template <int DisplacementDim>
diff --git a/ProcessLib/TH2M/TH2MProcessData.h b/ProcessLib/TH2M/TH2MProcessData.h
index 8aa138d9758cbc32bc70827c3eba720b67e2d77e..e23f42d8cff647772698b2a9a1237d88c4fc4c9a 100644
--- a/ProcessLib/TH2M/TH2MProcessData.h
+++ b/ProcessLib/TH2M/TH2MProcessData.h
@@ -52,7 +52,6 @@ struct TH2MProcessData
 
     const bool use_TaylorHood_elements;
 
-    MeshLib::PropertyVector<double>* element_saturation = nullptr;
     MeshLib::PropertyVector<double>* gas_pressure_interpolated = nullptr;
     MeshLib::PropertyVector<double>* capillary_pressure_interpolated = nullptr;
     MeshLib::PropertyVector<double>* temperature_interpolated = nullptr;
diff --git a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h
index f76508b2bc1dae3577536257e34bfe3272e73ef8..3a9c9776320d3184d3037b359871fe0f2ee182c9 100644
--- a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h
+++ b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h
@@ -541,14 +541,6 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
     unsigned const n_integration_points =
         this->integration_method_.getNumberOfPoints();
 
-    double saturation_avg = 0;
-    double porosity_avg = 0;
-    double liquid_density_avg = 0;
-    double viscosity_avg = 0;
-
-    using KV = MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
-    KV sigma_avg = KV::Zero();
-
     typename ConstitutiveTraits::ConstitutiveSetting constitutive_setting;
 
     auto models = ConstitutiveTraits::createConstitutiveModels(
@@ -603,29 +595,7 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
                                   eps, current_state, this->prev_states_[ip],
                                   this->material_states_[ip], tmp, output_data,
                                   CD);
-
-        saturation_avg += std::get<SaturationData>(current_state).S_L;
-        porosity_avg += std::get<PorosityData>(current_state).phi;
-
-        liquid_density_avg += std::get<LiquidDensityData>(output_data).rho_LR;
-        viscosity_avg += std::get<LiquidViscosityData>(output_data).viscosity;
-        sigma_avg += ConstitutiveTraits::ConstitutiveSetting::statefulStress(
-            current_state);
     }
-    saturation_avg /= n_integration_points;
-    porosity_avg /= n_integration_points;
-    viscosity_avg /= n_integration_points;
-    liquid_density_avg /= n_integration_points;
-    sigma_avg /= n_integration_points;
-
-    (*process_data.element_saturation)[e_id] = saturation_avg;
-    (*process_data.element_porosity)[e_id] = porosity_avg;
-    (*process_data.element_liquid_density)[e_id] = liquid_density_avg;
-    (*process_data.element_viscosity)[e_id] = viscosity_avg;
-
-    Eigen::Map<KV>(
-        &(*process_data.element_stresses)[e_id * KV::RowsAtCompileTime]) =
-        MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_avg);
 
     NumLib::interpolateToHigherOrderNodes<
         ShapeFunction, typename ShapeFunctionDisplacement::MeshElement,
diff --git a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.cpp b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.cpp
index b2961f64854c1e46eb39059cf492537c02a444f7..f00e5791de941bd16e1d6370c1e5a51ecca93f6c 100644
--- a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.cpp
+++ b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.cpp
@@ -18,6 +18,7 @@
 #include "MeshLib/Utils/getOrCreateMeshProperty.h"
 #include "NumLib/DOF/DOFTableUtil.h"
 #include "ProcessLib/Deformation/SolidMaterialInternalToSecondaryVariables.h"
+#include "ProcessLib/Output/CellAverageAlgorithm.h"
 #include "ProcessLib/Process.h"
 #include "ProcessLib/Reflection/ReflectionForExtrapolation.h"
 #include "ProcessLib/Reflection/ReflectionForIPWriters.h"
@@ -156,29 +157,6 @@ void ThermoRichardsMechanicsProcess<DisplacementDim, ConstitutiveTraits>::
             process_data_.solid_materials, local_assemblers_,
             _integration_point_writer, integration_order);
 
-    process_data_.element_saturation = MeshLib::getOrCreateMeshProperty<double>(
-        const_cast<MeshLib::Mesh&>(mesh), "saturation_avg",
-        MeshLib::MeshItemType::Cell, 1);
-
-    process_data_.element_porosity = MeshLib::getOrCreateMeshProperty<double>(
-        const_cast<MeshLib::Mesh&>(mesh), "porosity_avg",
-        MeshLib::MeshItemType::Cell, 1);
-
-    process_data_.element_liquid_density =
-        MeshLib::getOrCreateMeshProperty<double>(
-            const_cast<MeshLib::Mesh&>(mesh), "liquid_density_avg",
-            MeshLib::MeshItemType::Cell, 1);
-
-    process_data_.element_viscosity = MeshLib::getOrCreateMeshProperty<double>(
-        const_cast<MeshLib::Mesh&>(mesh), "viscosity_avg",
-        MeshLib::MeshItemType::Cell, 1);
-
-    process_data_.element_stresses = MeshLib::getOrCreateMeshProperty<double>(
-        const_cast<MeshLib::Mesh&>(mesh), "stress_avg",
-        MeshLib::MeshItemType::Cell,
-        MathLib::KelvinVector::KelvinVectorType<
-            DisplacementDim>::RowsAtCompileTime);
-
     process_data_.pressure_interpolated =
         MeshLib::getOrCreateMeshProperty<double>(
             const_cast<MeshLib::Mesh&>(mesh), "pressure_interpolated",
@@ -309,6 +287,8 @@ void ThermoRichardsMechanicsProcess<DisplacementDim, ConstitutiveTraits>::
         &LocalAssemblerIF::computeSecondaryVariable, local_assemblers_,
         pv.getActiveElementIDs(), getDOFTables(x.size()), t, dt, x, x_prev,
         process_id);
+
+    computeCellAverages<DisplacementDim>(cell_average_data_, local_assemblers_);
 }
 
 template <int DisplacementDim, typename ConstitutiveTraits>
diff --git a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcessData.h b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcessData.h
index b7933599c8f2e59e167197f0fb5ba65709c9615f..b01f9d61280d1838bb0f3f03c96fbc652ff3bb2e 100644
--- a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcessData.h
+++ b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcessData.h
@@ -55,11 +55,6 @@ struct ThermoRichardsMechanicsProcessData
     InitializePorosityFromMediumProperty const
         initialize_porosity_from_medium_property;
 
-    MeshLib::PropertyVector<double>* element_saturation = nullptr;
-    MeshLib::PropertyVector<double>* element_porosity = nullptr;
-    MeshLib::PropertyVector<double>* element_liquid_density = nullptr;
-    MeshLib::PropertyVector<double>* element_viscosity = nullptr;
-    MeshLib::PropertyVector<double>* element_stresses = nullptr;
     MeshLib::PropertyVector<double>* temperature_interpolated = nullptr;
     MeshLib::PropertyVector<double>* pressure_interpolated = nullptr;
 
diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small.prj b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small.prj
index 02b646e0280fefaa383070b7ab3e5ecf70a0f4e9..ae7609e210976b5e6e387e0692be8a5b7fefe92a 100644
--- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small.prj
+++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small.prj
@@ -205,6 +205,7 @@
                 <variable>NodalForces</variable>
                 <variable>pressure_interpolated</variable>
                 <variable>saturation_avg</variable>
+                <variable>sigma_avg</variable>
             </variables>
             <suffix>_ts_{:timestep}_t_{:time}</suffix>
         </output>
@@ -354,6 +355,12 @@
             <absolute_tolerance>1e-8</absolute_tolerance>
             <relative_tolerance>0</relative_tolerance>
         </vtkdiff>
+        <vtkdiff>
+            <regex>RichardsFlow_2d_small_ts_.*.vtu</regex>
+            <field>sigma_avg</field>
+            <absolute_tolerance>2.8e-12</absolute_tolerance>
+            <relative_tolerance>0</relative_tolerance>
+        </vtkdiff>
         <vtkdiff>
             <regex>RichardsFlow_2d_small_ts_.*.vtu</regex>
             <field>epsilon</field>
diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_0_t_0.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_0_t_0.000000.vtu
index ce7d33c10c99f0e8f4d84829f3445c123dad298b..b9dd3ced86605adba1dd2e32593ced6d8f31dc2e 100644
--- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_0_t_0.000000.vtu
+++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_0_t_0.000000.vtu
@@ -27,7 +27,7 @@
         <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0"                    RangeMax="0"                    offset="2724"                />
         <DataArray type="Float64" Name="porosity_avg" format="appended" RangeMin="0.38"                 RangeMax="0.38"                 offset="2788"                />
         <DataArray type="Float64" Name="saturation_avg" format="appended" RangeMin="0.37116239774"        RangeMax="0.37116239774"        offset="2864"                />
-        <DataArray type="Float64" Name="stress_avg" NumberOfComponents="4" format="appended" RangeMin="0"                    RangeMax="0"                    offset="2940"                />
+        <DataArray type="Float64" Name="sigma_avg" NumberOfComponents="4" format="appended" RangeMin="0"                    RangeMax="0"                    offset="2940"                />
       </CellData>
       <Points>
         <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="2.00039996"           offset="3020"                />
diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_104_t_2000.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_104_t_2000.000000.vtu
index 91e4affb9790b53ca8cf1bc9a2c0ccdcea2fb8b9..784c259c3e415afc5a54505a106315f4ca8134a2 100644
--- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_104_t_2000.000000.vtu
+++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_104_t_2000.000000.vtu
@@ -27,7 +27,7 @@
         <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0"                    RangeMax="0"                    offset="75392"               />
         <DataArray type="Float64" Name="porosity_avg" format="appended" RangeMin="0.38"                 RangeMax="0.38"                 offset="75456"               />
         <DataArray type="Float64" Name="saturation_avg" format="appended" RangeMin="0.370119056"          RangeMax="0.95"                 offset="75532"               />
-        <DataArray type="Float64" Name="stress_avg" NumberOfComponents="4" format="appended" RangeMin="26.927710272"         RangeMax="882.33307117"         offset="75932"               />
+        <DataArray type="Float64" Name="sigma_avg" NumberOfComponents="4" format="appended" RangeMin="26.927710272"         RangeMax="882.33307117"         offset="75932"               />
       </CellData>
       <Points>
         <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="2.00039996"           offset="78984"               />
diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_19_t_318.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_19_t_318.000000.vtu
index a9f53bc77b87934f9129e9a6bacb496c2b8786b7..339fa819ca5b7cd771f00f5b59f6d492c6efdd4e 100644
--- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_19_t_318.000000.vtu
+++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_19_t_318.000000.vtu
@@ -27,7 +27,7 @@
         <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0"                    RangeMax="0"                    offset="75980"               />
         <DataArray type="Float64" Name="porosity_avg" format="appended" RangeMin="0.38"                 RangeMax="0.38"                 offset="76044"               />
         <DataArray type="Float64" Name="saturation_avg" format="appended" RangeMin="0.36527958308"        RangeMax="0.94979723169"        offset="76120"               />
-        <DataArray type="Float64" Name="stress_avg" NumberOfComponents="4" format="appended" RangeMin="5.3979639061"         RangeMax="679.33577683"         offset="76536"               />
+        <DataArray type="Float64" Name="sigma_avg" NumberOfComponents="4" format="appended" RangeMin="5.3979639061"         RangeMax="679.33577683"         offset="76536"               />
       </CellData>
       <Points>
         <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="2.00039996"           offset="79632"               />
diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_1_t_1.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_1_t_1.000000.vtu
index 703fec56f30d26feaea178f6af484011bd1fb8ae..83ab412ef9244b181d7c849d1dbd3a5ca5a62047 100644
--- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_1_t_1.000000.vtu
+++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_1_t_1.000000.vtu
@@ -27,7 +27,7 @@
         <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0"                    RangeMax="0"                    offset="76136"               />
         <DataArray type="Float64" Name="porosity_avg" format="appended" RangeMin="0.38"                 RangeMax="0.38"                 offset="76200"               />
         <DataArray type="Float64" Name="saturation_avg" format="appended" RangeMin="0.29769080019"        RangeMax="0.62977672082"        offset="76276"               />
-        <DataArray type="Float64" Name="stress_avg" NumberOfComponents="4" format="appended" RangeMin="0.13873837932"        RangeMax="135.9312842"          offset="76588"               />
+        <DataArray type="Float64" Name="sigma_avg" NumberOfComponents="4" format="appended" RangeMin="0.13873837932"        RangeMax="135.9312842"          offset="76588"               />
       </CellData>
       <Points>
         <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="2.00039996"           offset="79764"               />
diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_29_t_518.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_29_t_518.000000.vtu
index b6202102d148aff3be099e12a886f1f0c30d82a5..8cd523c95c1d55f474143683492dbbe1d01ce4da 100644
--- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_29_t_518.000000.vtu
+++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_29_t_518.000000.vtu
@@ -27,7 +27,7 @@
         <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0"                    RangeMax="0"                    offset="75708"               />
         <DataArray type="Float64" Name="porosity_avg" format="appended" RangeMin="0.38"                 RangeMax="0.38"                 offset="75772"               />
         <DataArray type="Float64" Name="saturation_avg" format="appended" RangeMin="0.34825015959"        RangeMax="0.94997120127"        offset="75848"               />
-        <DataArray type="Float64" Name="stress_avg" NumberOfComponents="4" format="appended" RangeMin="4.5539936988"         RangeMax="734.83684114"         offset="76232"               />
+        <DataArray type="Float64" Name="sigma_avg" NumberOfComponents="4" format="appended" RangeMin="4.5539936988"         RangeMax="734.83684114"         offset="76232"               />
       </CellData>
       <Points>
         <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="2.00039996"           offset="79300"               />
diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_2_t_3.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_2_t_3.000000.vtu
index 264abfd42822389170b5ca56846095d4f908c79b..8d6a382bdfa77bec33f9fe30f67bc3776fa18e0e 100644
--- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_2_t_3.000000.vtu
+++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_2_t_3.000000.vtu
@@ -27,7 +27,7 @@
         <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0"                    RangeMax="0"                    offset="76684"               />
         <DataArray type="Float64" Name="porosity_avg" format="appended" RangeMin="0.38"                 RangeMax="0.38"                 offset="76748"               />
         <DataArray type="Float64" Name="saturation_avg" format="appended" RangeMin="0.2993620397"         RangeMax="0.63504683156"        offset="76824"               />
-        <DataArray type="Float64" Name="stress_avg" NumberOfComponents="4" format="appended" RangeMin="0.14376166452"        RangeMax="136.18913915"         offset="77160"               />
+        <DataArray type="Float64" Name="sigma_avg" NumberOfComponents="4" format="appended" RangeMin="0.14376166452"        RangeMax="136.18913915"         offset="77160"               />
       </CellData>
       <Points>
         <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="2.00039996"           offset="80340"               />
diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_39_t_718.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_39_t_718.000000.vtu
index cce79dd592d7076174a838fd5fea7fb9ed50e6a6..39c90e840f4efec0379654679bb3861f4f311f70 100644
--- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_39_t_718.000000.vtu
+++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_39_t_718.000000.vtu
@@ -27,7 +27,7 @@
         <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0"                    RangeMax="0"                    offset="75512"               />
         <DataArray type="Float64" Name="porosity_avg" format="appended" RangeMin="0.38"                 RangeMax="0.38"                 offset="75576"               />
         <DataArray type="Float64" Name="saturation_avg" format="appended" RangeMin="0.36726703994"        RangeMax="0.94999968705"        offset="75652"               />
-        <DataArray type="Float64" Name="stress_avg" NumberOfComponents="4" format="appended" RangeMin="13.395372625"         RangeMax="800.64448175"         offset="76044"               />
+        <DataArray type="Float64" Name="sigma_avg" NumberOfComponents="4" format="appended" RangeMin="13.395372625"         RangeMax="800.64448175"         offset="76044"               />
       </CellData>
       <Points>
         <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="2.00039996"           offset="79096"               />
diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_3_t_8.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_3_t_8.000000.vtu
index 47475355ff2c6067ba7cd4c9ec4d48cab715432a..03367285fba66847888693cbb8e091bf37a32a47 100644
--- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_3_t_8.000000.vtu
+++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_3_t_8.000000.vtu
@@ -27,7 +27,7 @@
         <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0"                    RangeMax="0"                    offset="77176"               />
         <DataArray type="Float64" Name="porosity_avg" format="appended" RangeMin="0.38"                 RangeMax="0.38"                 offset="77240"               />
         <DataArray type="Float64" Name="saturation_avg" format="appended" RangeMin="0.30376757647"        RangeMax="0.64829678481"        offset="77316"               />
-        <DataArray type="Float64" Name="stress_avg" NumberOfComponents="4" format="appended" RangeMin="0.1376442936"         RangeMax="137.60530976"         offset="77680"               />
+        <DataArray type="Float64" Name="sigma_avg" NumberOfComponents="4" format="appended" RangeMin="0.1376442936"         RangeMax="137.60530976"         offset="77680"               />
       </CellData>
       <Points>
         <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="2.00039996"           offset="80876"               />
diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_49_t_918.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_49_t_918.000000.vtu
index 6f67c0b548634c9115a5fed6b20ebfcb01402d18..19d052c30121f83cea505bf78a9982c5943b1e19 100644
--- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_49_t_918.000000.vtu
+++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_49_t_918.000000.vtu
@@ -27,7 +27,7 @@
         <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0"                    RangeMax="0"                    offset="75492"               />
         <DataArray type="Float64" Name="porosity_avg" format="appended" RangeMin="0.38"                 RangeMax="0.38"                 offset="75556"               />
         <DataArray type="Float64" Name="saturation_avg" format="appended" RangeMin="0.35432299817"        RangeMax="0.94999995825"        offset="75632"               />
-        <DataArray type="Float64" Name="stress_avg" NumberOfComponents="4" format="appended" RangeMin="14.165778776"         RangeMax="815.19348278"         offset="76040"               />
+        <DataArray type="Float64" Name="sigma_avg" NumberOfComponents="4" format="appended" RangeMin="14.165778776"         RangeMax="815.19348278"         offset="76040"               />
       </CellData>
       <Points>
         <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="2.00039996"           offset="79108"               />
diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_4_t_18.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_4_t_18.000000.vtu
index a993960c70948092b1faaf01df3d1f6e08fa6d60..ea244a19a65992faa3fbde95ca93428d0e10d4df 100644
--- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_4_t_18.000000.vtu
+++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_4_t_18.000000.vtu
@@ -27,7 +27,7 @@
         <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0"                    RangeMax="0"                    offset="77412"               />
         <DataArray type="Float64" Name="porosity_avg" format="appended" RangeMin="0.38"                 RangeMax="0.38"                 offset="77476"               />
         <DataArray type="Float64" Name="saturation_avg" format="appended" RangeMin="0.3136966799"         RangeMax="0.67496680755"        offset="77552"               />
-        <DataArray type="Float64" Name="stress_avg" NumberOfComponents="4" format="appended" RangeMin="0.05750444578"        RangeMax="143.95094316"         offset="77920"               />
+        <DataArray type="Float64" Name="sigma_avg" NumberOfComponents="4" format="appended" RangeMin="0.05750444578"        RangeMax="143.95094316"         offset="77920"               />
       </CellData>
       <Points>
         <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="2.00039996"           offset="81096"               />
diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_59_t_1118.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_59_t_1118.000000.vtu
index 6a6c5a26731da1bbc6ed4b13f73355820708dfdf..22d5a98f6bd28f913bdb79ca20b83842aa85eabe 100644
--- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_59_t_1118.000000.vtu
+++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_59_t_1118.000000.vtu
@@ -27,7 +27,7 @@
         <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0"                    RangeMax="0"                    offset="75280"               />
         <DataArray type="Float64" Name="porosity_avg" format="appended" RangeMin="0.38"                 RangeMax="0.38"                 offset="75344"               />
         <DataArray type="Float64" Name="saturation_avg" format="appended" RangeMin="0.35576146526"        RangeMax="0.94999999999"        offset="75420"               />
-        <DataArray type="Float64" Name="stress_avg" NumberOfComponents="4" format="appended" RangeMin="11.48782448"          RangeMax="835.43491611"         offset="75816"               />
+        <DataArray type="Float64" Name="sigma_avg" NumberOfComponents="4" format="appended" RangeMin="11.48782448"          RangeMax="835.43491611"         offset="75816"               />
       </CellData>
       <Points>
         <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="2.00039996"           offset="78868"               />
diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_5_t_38.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_5_t_38.000000.vtu
index 89bef788e6863590eca82f17fd420bdec6eddb31..4ee3963156d1d1d9afa71b236f7e601a5c990728 100644
--- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_5_t_38.000000.vtu
+++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_5_t_38.000000.vtu
@@ -27,7 +27,7 @@
         <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0"                    RangeMax="0"                    offset="77228"               />
         <DataArray type="Float64" Name="porosity_avg" format="appended" RangeMin="0.38"                 RangeMax="0.38"                 offset="77292"               />
         <DataArray type="Float64" Name="saturation_avg" format="appended" RangeMin="0.33831365519"        RangeMax="0.72704656573"        offset="77368"               />
-        <DataArray type="Float64" Name="stress_avg" NumberOfComponents="4" format="appended" RangeMin="0.088480039377"       RangeMax="168.62458098"         offset="77724"               />
+        <DataArray type="Float64" Name="sigma_avg" NumberOfComponents="4" format="appended" RangeMin="0.088480039377"       RangeMax="168.62458098"         offset="77724"               />
       </CellData>
       <Points>
         <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="2.00039996"           offset="80888"               />
diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_69_t_1318.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_69_t_1318.000000.vtu
index 38a1c9d9cd6f7881f8eb55a5e5c00ebcb083318e..879a21c0febd456a54761f6f1e253f1f6dc14c2d 100644
--- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_69_t_1318.000000.vtu
+++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_69_t_1318.000000.vtu
@@ -27,7 +27,7 @@
         <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0"                    RangeMax="0"                    offset="75160"               />
         <DataArray type="Float64" Name="porosity_avg" format="appended" RangeMin="0.38"                 RangeMax="0.38"                 offset="75224"               />
         <DataArray type="Float64" Name="saturation_avg" format="appended" RangeMin="0.36713370976"        RangeMax="0.95"                 offset="75300"               />
-        <DataArray type="Float64" Name="stress_avg" NumberOfComponents="4" format="appended" RangeMin="20.824043994"         RangeMax="855.98600706"         offset="75684"               />
+        <DataArray type="Float64" Name="sigma_avg" NumberOfComponents="4" format="appended" RangeMin="20.824043994"         RangeMax="855.98600706"         offset="75684"               />
       </CellData>
       <Points>
         <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="2.00039996"           offset="78732"               />
diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_6_t_58.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_6_t_58.000000.vtu
index 99d0546d156b5829cdd4190ea931d540aaeb2e16..0848e3a65dd0f464bab7f1485d9c458e53752c85 100644
--- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_6_t_58.000000.vtu
+++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_6_t_58.000000.vtu
@@ -27,7 +27,7 @@
         <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0"                    RangeMax="0"                    offset="77184"               />
         <DataArray type="Float64" Name="porosity_avg" format="appended" RangeMin="0.38"                 RangeMax="0.38"                 offset="77248"               />
         <DataArray type="Float64" Name="saturation_avg" format="appended" RangeMin="0.36688449277"        RangeMax="0.77425673495"        offset="77324"               />
-        <DataArray type="Float64" Name="stress_avg" NumberOfComponents="4" format="appended" RangeMin="0.12160563802"        RangeMax="198.72934807"         offset="77688"               />
+        <DataArray type="Float64" Name="sigma_avg" NumberOfComponents="4" format="appended" RangeMin="0.12160563802"        RangeMax="198.72934807"         offset="77688"               />
       </CellData>
       <Points>
         <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="2.00039996"           offset="80860"               />
diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_79_t_1518.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_79_t_1518.000000.vtu
index f27991cbb38b44724edb341889e2c90d7bace39e..f4f39ef92a85ae77571f765d3aa848571f1ee713 100644
--- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_79_t_1518.000000.vtu
+++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_79_t_1518.000000.vtu
@@ -27,7 +27,7 @@
         <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0"                    RangeMax="0"                    offset="75500"               />
         <DataArray type="Float64" Name="porosity_avg" format="appended" RangeMin="0.38"                 RangeMax="0.38"                 offset="75564"               />
         <DataArray type="Float64" Name="saturation_avg" format="appended" RangeMin="0.35812247015"        RangeMax="0.95"                 offset="75640"               />
-        <DataArray type="Float64" Name="stress_avg" NumberOfComponents="4" format="appended" RangeMin="21.748341205"         RangeMax="864.27964012"         offset="76056"               />
+        <DataArray type="Float64" Name="sigma_avg" NumberOfComponents="4" format="appended" RangeMin="21.748341205"         RangeMax="864.27964012"         offset="76056"               />
       </CellData>
       <Points>
         <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="2.00039996"           offset="79124"               />
diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_7_t_78.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_7_t_78.000000.vtu
index 287a168431e5625802b3e1f6c2e803bec2f66f96..017f00e2fab97fdcd41d6fbf77cfbdfb680ca8d9 100644
--- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_7_t_78.000000.vtu
+++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_7_t_78.000000.vtu
@@ -27,7 +27,7 @@
         <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0"                    RangeMax="0"                    offset="77100"               />
         <DataArray type="Float64" Name="porosity_avg" format="appended" RangeMin="0.38"                 RangeMax="0.38"                 offset="77164"               />
         <DataArray type="Float64" Name="saturation_avg" format="appended" RangeMin="0.36448867031"        RangeMax="0.81710537713"        offset="77240"               />
-        <DataArray type="Float64" Name="stress_avg" NumberOfComponents="4" format="appended" RangeMin="0.12410710626"        RangeMax="229.0455663"          offset="77604"               />
+        <DataArray type="Float64" Name="sigma_avg" NumberOfComponents="4" format="appended" RangeMin="0.12410710626"        RangeMax="229.0455663"          offset="77604"               />
       </CellData>
       <Points>
         <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="2.00039996"           offset="80760"               />
diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_89_t_1718.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_89_t_1718.000000.vtu
index e71fddbc39c07dfe55b3a0d0168411ffe4344524..d012b085e7a1d613d3f195b7e99f643ef0fa115e 100644
--- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_89_t_1718.000000.vtu
+++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_89_t_1718.000000.vtu
@@ -27,7 +27,7 @@
         <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0"                    RangeMax="0"                    offset="75160"               />
         <DataArray type="Float64" Name="porosity_avg" format="appended" RangeMin="0.38"                 RangeMax="0.38"                 offset="75224"               />
         <DataArray type="Float64" Name="saturation_avg" format="appended" RangeMin="0.35485167754"        RangeMax="0.95"                 offset="75300"               />
-        <DataArray type="Float64" Name="stress_avg" NumberOfComponents="4" format="appended" RangeMin="23.098608117"         RangeMax="870.28500247"         offset="75708"               />
+        <DataArray type="Float64" Name="sigma_avg" NumberOfComponents="4" format="appended" RangeMin="23.098608117"         RangeMax="870.28500247"         offset="75708"               />
       </CellData>
       <Points>
         <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="2.00039996"           offset="78748"               />
diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_8_t_98.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_8_t_98.000000.vtu
index 0b00240c9165a58805bf52a5311ad06ee3d0f10f..bd511e3f76102d7a3f70c2833414856ec35385de 100644
--- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_8_t_98.000000.vtu
+++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_8_t_98.000000.vtu
@@ -27,7 +27,7 @@
         <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0"                    RangeMax="0"                    offset="77440"               />
         <DataArray type="Float64" Name="porosity_avg" format="appended" RangeMin="0.38"                 RangeMax="0.38"                 offset="77504"               />
         <DataArray type="Float64" Name="saturation_avg" format="appended" RangeMin="0.3577465938"         RangeMax="0.85566017945"        offset="77580"               />
-        <DataArray type="Float64" Name="stress_avg" NumberOfComponents="4" format="appended" RangeMin="0.14825098694"        RangeMax="262.70856268"         offset="77948"               />
+        <DataArray type="Float64" Name="sigma_avg" NumberOfComponents="4" format="appended" RangeMin="0.14825098694"        RangeMax="262.70856268"         offset="77948"               />
       </CellData>
       <Points>
         <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="2.00039996"           offset="81116"               />
diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_99_t_1918.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_99_t_1918.000000.vtu
index 69eacf9c7f67f01129c498f609ec6aff8f24c380..df0b16698e09284e7c50280ed86ce7199b6af490 100644
--- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_99_t_1918.000000.vtu
+++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_99_t_1918.000000.vtu
@@ -27,7 +27,7 @@
         <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0"                    RangeMax="0"                    offset="75312"               />
         <DataArray type="Float64" Name="porosity_avg" format="appended" RangeMin="0.38"                 RangeMax="0.38"                 offset="75376"               />
         <DataArray type="Float64" Name="saturation_avg" format="appended" RangeMin="0.36376927165"        RangeMax="0.95"                 offset="75452"               />
-        <DataArray type="Float64" Name="stress_avg" NumberOfComponents="4" format="appended" RangeMin="25.794361378"         RangeMax="878.71459569"         offset="75860"               />
+        <DataArray type="Float64" Name="sigma_avg" NumberOfComponents="4" format="appended" RangeMin="25.794361378"         RangeMax="878.71459569"         offset="75860"               />
       </CellData>
       <Points>
         <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="2.00039996"           offset="78904"               />
diff --git a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_9_t_118.000000.vtu b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_9_t_118.000000.vtu
index 8add72d94dbeee2a45d94a273fe2155587f01f73..47c3c4b5601087f7ec4a897adaef98861e5bb6d9 100644
--- a/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_9_t_118.000000.vtu
+++ b/Tests/Data/ThermoRichardsMechanics/RichardsFlow2D/RichardsFlow_2d_small_ts_9_t_118.000000.vtu
@@ -27,7 +27,7 @@
         <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0"                    RangeMax="0"                    offset="77608"               />
         <DataArray type="Float64" Name="porosity_avg" format="appended" RangeMin="0.38"                 RangeMax="0.38"                 offset="77672"               />
         <DataArray type="Float64" Name="saturation_avg" format="appended" RangeMin="0.35168479802"        RangeMax="0.88895963149"        offset="77748"               />
-        <DataArray type="Float64" Name="stress_avg" NumberOfComponents="4" format="appended" RangeMin="0.093494384573"       RangeMax="304.50768442"         offset="78128"               />
+        <DataArray type="Float64" Name="sigma_avg" NumberOfComponents="4" format="appended" RangeMin="0.093494384573"       RangeMax="304.50768442"         offset="78128"               />
       </CellData>
       <Points>
         <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="2.00039996"           offset="81280"               />
diff --git a/Tests/ProcessLib/TestReflectIPData.cpp b/Tests/ProcessLib/TestReflectIPData.cpp
index 5ace325219c359237ead1705635c6f1ac03bfd38..ad6fbe060b630023c3c2797620e445957a4f70ba 100644
--- a/Tests/ProcessLib/TestReflectIPData.cpp
+++ b/Tests/ProcessLib/TestReflectIPData.cpp
@@ -14,6 +14,9 @@
 #include <numeric>
 #include <random>
 
+#include "MeshToolsLib/MeshGenerators/MeshGenerator.h"
+#include "ProcessLib/Output/CellAverageAlgorithm.h"
+#include "ProcessLib/Reflection/ReflectionForIPWriters.h"
 #include "ProcessLib/Reflection/ReflectionIPData.h"
 #include "ProcessLib/Reflection/ReflectionSetIPData.h"
 
@@ -130,6 +133,8 @@ struct LocAsmIF
             makeReflectionData(&LocAsmIF::ip_data_level1),
             makeReflectionData(&LocAsmIF::ip_data_level1b)};
     }
+
+    static auto getReflectionDataForOutput() { return reflect(); }
 };
 
 template <int dim>
@@ -251,6 +256,7 @@ std::vector<double> initKelvin(LocAsmIF<dim>& loc_asm,
         for (std::size_t ip = 0; ip < num_int_pts; ++ip)
         {
             ip_data_accessor(loc_asm, ip) =
+                // the local assembler stores Kelvin vector data...
                 MathLib::KelvinVector::symmetricTensorToKelvinVector(
                     Eigen::Vector<double, kv_size>::LinSpaced(
                         kv_size, ip * kv_size + start_value,
@@ -268,6 +274,7 @@ std::vector<double> initKelvin(LocAsmIF<dim>& loc_asm,
     }
 
     // prepare reference data
+    // ... the reference data used in the tests is not Kelvin-mapped!
     std::vector<double> vector_expected(num_int_pts * kv_size);
     iota(begin(vector_expected), end(vector_expected), start_value);
     return vector_expected;
@@ -329,6 +336,10 @@ std::vector<double> initMatrix(LocAsmIF<dim>& loc_asm,
 template <int dim>
 struct ReferenceData
 {
+private:
+    ReferenceData() = default;
+
+public:
     std::vector<double> scalar;
     std::vector<double> vector;
     std::vector<double> kelvin;
@@ -348,6 +359,8 @@ struct ReferenceData
     std::vector<double> scalar2b;
     std::vector<double> scalar3b;
 
+    // Computes reference data and initializes the internal (integration point)
+    // data of the passed \c loc_asm.
     static ReferenceData<dim> create(LocAsmIF<dim>& loc_asm,
                                      bool const for_read_test)
     {
@@ -484,6 +497,45 @@ struct ReferenceData
     }
 };
 
+template <int dim>
+void checkAll(ReferenceData<dim> const& ref, auto&& checker,
+              bool const check_level_0_data = true)
+{
+    auto constexpr kv_size =
+        MathLib::KelvinVector::kelvin_vector_dimensions(dim);
+
+    // The storage of "level 0 data" (i.e., non-reflected data directly in the
+    // local assembler, e.g. a std::vector<double> of integration point data of
+    // a scalar value) is not compatible with the implemented IP data input
+    // logic. Therefore we don't test level 0 data in that case.
+    // This incompatibility is not a problem in practice, as such data is not
+    // used ATM in OGS's local assemblers.
+    if (check_level_0_data)
+    {
+        checker("scalar", 1, ref.scalar);
+        checker("vector", dim, ref.vector);
+        checker("kelvin", kv_size, ref.kelvin);
+    }
+
+    // level 1
+    checker("scalar1", 1, ref.scalar1);
+    checker("vector1", dim, ref.vector1);
+    checker("kelvin1", kv_size, ref.kelvin1);
+
+    // level 3
+    checker("scalar3", 1, ref.scalar3);
+    checker("vector3", dim, ref.vector3);
+    checker("kelvin3", kv_size, ref.kelvin3);
+    checker("matrix3", dim * 4, ref.matrix3);
+    checker("matrix3_1", dim * 4, ref.matrix3_1);
+    checker("matrix3_2", 4, ref.matrix3_2);
+
+    // b levels
+    checker("scalar1b", 1, ref.scalar1b);
+    checker("scalar2b", 1, ref.scalar2b);
+    checker("scalar3b", 1, ref.scalar3b);
+}
+
 template <class Dim>
 struct ProcessLib_ReflectIPData : ::testing::Test
 {
@@ -499,8 +551,6 @@ TYPED_TEST_SUITE(ProcessLib_ReflectIPData, ProcessLib_ReflectIPData_TestCases);
 TYPED_TEST(ProcessLib_ReflectIPData, ReadTest)
 {
     constexpr int dim = TypeParam::value;
-    auto constexpr kv_size =
-        MathLib::KelvinVector::kelvin_vector_dimensions(dim);
 
     using LocAsm = LocAsmIF<dim>;
 
@@ -548,28 +598,7 @@ TYPED_TEST(ProcessLib_ReflectIPData, ReadTest)
             << "Values differ for ip data with name '" << name << "'";
     };
 
-    // level 0
-    check("scalar", 1, ref.scalar);
-    check("vector", dim, ref.vector);
-    check("kelvin", kv_size, ref.kelvin);
-
-    // level 1
-    check("scalar1", 1, ref.scalar1);
-    check("vector1", dim, ref.vector1);
-    check("kelvin1", kv_size, ref.kelvin1);
-
-    // level 3
-    check("scalar3", 1, ref.scalar3);
-    check("vector3", dim, ref.vector3);
-    check("kelvin3", kv_size, ref.kelvin3);
-    check("matrix3", dim * 4, ref.matrix3);
-    check("matrix3_1", dim * 4, ref.matrix3_1);
-    check("matrix3_2", 4, ref.matrix3_2);
-
-    // b levels
-    check("scalar1b", 1, ref.scalar1b);
-    check("scalar2b", 1, ref.scalar2b);
-    check("scalar3b", 1, ref.scalar3b);
+    checkAll(ref, check);
 }
 
 TYPED_TEST(ProcessLib_ReflectIPData, WriteTest)
@@ -583,7 +612,8 @@ TYPED_TEST(ProcessLib_ReflectIPData, WriteTest)
 
     auto const ref = ReferenceData<dim>::create(loc_asm, false);
 
-    // data getters - used for checks //////////////////////////////////////////
+    // set up data getters - used for checks ///////////////////////////////////
+    // that critically relies on the read test above being successful!
 
     std::map<std::string, NumCompAndFunction<dim>>
         map_name_to_num_comp_and_function;
@@ -603,7 +633,7 @@ TYPED_TEST(ProcessLib_ReflectIPData, WriteTest)
     // checks //////////////////////////////////////////////////////////////////
 
     auto check = [&map_name_to_num_comp_and_function, &loc_asm](
-                     std::string const& name, std::size_t const size_expected,
+                     std::string const& name, unsigned const num_comp,
                      std::vector<double> const& values_plain)
     {
         auto const it = map_name_to_num_comp_and_function.find(name);
@@ -611,19 +641,34 @@ TYPED_TEST(ProcessLib_ReflectIPData, WriteTest)
         ASSERT_NE(map_name_to_num_comp_and_function.end(), it)
             << "No accessor found for ip data with name '" << name << "'";
 
-        auto const& [num_comp, fct] = it->second;
+        auto const& [num_comp_2, fct] = it->second;
+
+        // consistency checks
+        ASSERT_EQ(num_comp, num_comp_2);
+        ASSERT_EQ(0, values_plain.size() % num_comp);
+        auto const num_int_pts_expected = values_plain.size() / num_comp;
 
         EXPECT_THAT(fct(loc_asm), testing::Each(testing::IsNan()))
-            << "All values must be initialize to NaN in this unit test. Check "
+            << "All values must be initialized to NaN in this unit test. Check "
                "failed for ip data with name '"
             << name << "'";
 
-        auto const size = ProcessLib::Reflection::reflectSetIPData<dim>(
-            name, values_plain.data(), loc_asm.ip_data_level1);
+        // function under test /////////////////////////////////////////////////
 
-        EXPECT_EQ(size_expected, size)
-            << "Unexpected size obtained for ip data with name '" << name
-            << "'";
+        auto const num_int_pts_actual =
+            ProcessLib::Reflection::reflectSetIPData<dim>(
+                name, values_plain.data(), loc_asm.ip_data_level1);
+        auto const num_int_pts_actual_1 =
+            ProcessLib::Reflection::reflectSetIPData<dim>(
+                name, values_plain.data(), loc_asm.ip_data_level1b);
+
+        // end function under test /////////////////////////////////////////////
+
+        ASSERT_EQ(num_int_pts_expected, num_int_pts_actual)
+            << "Unexpected number of integration points obtained for ip data "
+               "with name '"
+            << name << "'";
+        ASSERT_EQ(num_int_pts_expected, num_int_pts_actual_1);
 
         // check set values via round-trip with getter tested in previous unit
         // test
@@ -633,19 +678,7 @@ TYPED_TEST(ProcessLib_ReflectIPData, WriteTest)
             << "'";
     };
 
-    check("scalar1", num_int_pts, ref.scalar1);
-    check("vector1", num_int_pts, ref.vector1);
-    check("kelvin1", num_int_pts, ref.kelvin1);
-
-    check("scalar3", num_int_pts, ref.scalar3);
-    check("vector3", num_int_pts, ref.vector3);
-    check("kelvin3", num_int_pts, ref.kelvin3);
-    check("matrix3", num_int_pts, ref.matrix3);
-    check("matrix3_1", num_int_pts, ref.matrix3_1);
-    check("matrix3_2", num_int_pts, ref.matrix3_2);
-
-    check("scalar2b", num_int_pts, ref.scalar2b);
-    check("scalar3b", num_int_pts, ref.scalar3b);
+    checkAll(ref, check, false);
 }
 
 TYPED_TEST(ProcessLib_ReflectIPData, RawDataTypes)
@@ -687,3 +720,147 @@ TYPED_TEST(ProcessLib_ReflectIPData, RawDataTypes)
     static_assert(
         !PRD::is_raw_data<Eigen::Matrix<double, kv_size, kv_size>>::value);
 }
+
+TYPED_TEST(ProcessLib_ReflectIPData, IPWriterTest)
+{
+    constexpr int dim = TypeParam::value;
+
+    using LocAsm = LocAsmIF<dim>;
+
+    std::size_t const num_int_pts = 8;
+    std::vector<std::unique_ptr<LocAsm>> loc_asms;
+    // emplace...new is a workaround for an MSVC compiler warning:
+    // C:\Program Files (x86)\Microsoft Visual
+    // Studio\2019\Community\VC\Tools\MSVC\14.29.30133\include\memory(3382):
+    // warning C4244: 'argument': conversion from 'const uintptr_t' to 'const
+    // unsigned int', possible loss of data
+    // C:/Users/gitlab/gitlab/_b/bg4d5s_d/0/ogs/ogs/Tests/ProcessLib/
+    // TestReflectIPData.cpp(792):
+    // note: see reference to function template instantiation
+    // 'std::unique_ptr<LocAsm,std::default_delete<LocAsm>>
+    // std::make_unique<LocAsm,const size_t&,0>(const size_t &)' being compiled
+    loc_asms.emplace_back(new LocAsm{num_int_pts});
+
+    std::unique_ptr<MeshLib::Mesh> mesh{
+        MeshToolsLib::MeshGenerator::generateLineMesh(1.0, 1)};
+
+    auto const ref = ReferenceData<dim>::create(*loc_asms.front(), true);
+
+    // function under test /////////////////////////////////////////////////////
+
+    std::vector<std::unique_ptr<MeshLib::IntegrationPointWriter>>
+        integration_point_writers;
+    constexpr int integration_order = 0xc0ffee;  // dummy value
+    ProcessLib::Reflection::addReflectedIntegrationPointWriters<dim>(
+        LocAsm::getReflectionDataForOutput(), integration_point_writers,
+        integration_order, loc_asms);
+
+    // this finally invokes the IP writers
+    MeshLib::addIntegrationPointDataToMesh(*mesh, integration_point_writers);
+
+    // checks //////////////////////////////////////////////////////////////////
+
+    auto check = [&mesh](std::string const& name,
+                         unsigned const num_comp,
+                         std::vector<double> const& values_expected)
+    {
+        auto const& props = mesh->getProperties();
+
+        ASSERT_TRUE(props.existsPropertyVector<double>(name + "_ip"))
+            << "Property vector '" << name + "_ip"
+            << "' does not exist in the mesh.";
+
+        ASSERT_EQ(0, values_expected.size() % num_comp);
+        auto const num_int_pts = values_expected.size() / num_comp;
+
+        auto const& ip_data_actual =
+            *props.getPropertyVector<double>(name + "_ip");
+
+        ASSERT_EQ(MeshLib::MeshItemType::IntegrationPoint,
+                  ip_data_actual.getMeshItemType());
+        ASSERT_EQ(mesh->getNumberOfElements() * num_int_pts,
+                  ip_data_actual.getNumberOfTuples());
+        ASSERT_EQ(num_comp, ip_data_actual.getNumberOfGlobalComponents());
+
+        EXPECT_THAT(ip_data_actual,
+                    testing::Pointwise(testing::DoubleEq(), values_expected))
+            << "Values differ for ip data with name '" << name << "'";
+    };
+
+    checkAll(ref, check);
+}
+
+TYPED_TEST(ProcessLib_ReflectIPData, CellAverageTest)
+{
+    constexpr int dim = TypeParam::value;
+
+    using LocAsm = LocAsmIF<dim>;
+
+    std::size_t const num_int_pts = 8;
+    std::vector<std::unique_ptr<LocAsm>> loc_asms;
+    // emplace...new is a workaround for an MSVC compiler warning:
+    // C:\Program Files (x86)\Microsoft Visual
+    // Studio\2019\Community\VC\Tools\MSVC\14.29.30133\include\memory(3382):
+    // warning C4244: 'argument': conversion from 'const uintptr_t' to 'const
+    // unsigned int', possible loss of data
+    // C:/Users/gitlab/gitlab/_b/bg4d5s_d/0/ogs/ogs/Tests/ProcessLib/
+    // TestReflectIPData.cpp(792):
+    // note: see reference to function template instantiation
+    // 'std::unique_ptr<LocAsm,std::default_delete<LocAsm>>
+    // std::make_unique<LocAsm,const size_t&,0>(const size_t &)' being compiled
+    loc_asms.emplace_back(new LocAsm{num_int_pts});
+
+    std::unique_ptr<MeshLib::Mesh> mesh{
+        MeshToolsLib::MeshGenerator::generateLineMesh(1.0, 1)};
+
+    auto const ref = ReferenceData<dim>::create(*loc_asms.front(), true);
+
+    // function under test /////////////////////////////////////////////////////
+
+    ProcessLib::CellAverageData cell_average_data{*mesh};
+    computeCellAverages<dim>(cell_average_data, loc_asms);
+
+    // checks //////////////////////////////////////////////////////////////////
+
+    auto check = [&mesh](std::string const& name,
+                         unsigned const num_comp_expected,
+                         std::vector<double> const& ip_data)
+    {
+        ASSERT_EQ(num_int_pts * num_comp_expected, ip_data.size());
+
+        // TODO the implementation here in the unit test computing cell average
+        // reference data might be too similar to the production code. In fact,
+        // it's almost the same code.
+
+        // each integration point corresponds to a column in the mapped
+        // matrix, vector components/matrix entries are stored contiguously
+        // in memory.
+        Eigen::Map<const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic,
+                                       Eigen::ColMajor>>
+            ip_data_mat{ip_data.data(), num_comp_expected, num_int_pts};
+
+        Eigen::VectorXd const cell_avg_expected = ip_data_mat.rowwise().mean();
+
+        auto const& props = mesh->getProperties();
+
+        ASSERT_TRUE(props.existsPropertyVector<double>(name + "_avg"))
+            << "Property vector '" << name + "_avg"
+            << "' does not exist in the mesh.";
+
+        auto const& cell_avg_actual =
+            *props.getPropertyVector<double>(name + "_avg");
+
+        ASSERT_EQ(MeshLib::MeshItemType::Cell,
+                  cell_avg_actual.getMeshItemType());
+        ASSERT_EQ(mesh->getNumberOfElements(),
+                  cell_avg_actual.getNumberOfTuples());
+        ASSERT_EQ(num_comp_expected,
+                  cell_avg_actual.getNumberOfGlobalComponents());
+
+        EXPECT_THAT(cell_avg_actual,
+                    testing::Pointwise(testing::DoubleEq(), cell_avg_expected))
+            << "Values differ for cell average data with name '" << name << "'";
+    };
+
+    checkAll(ref, check);
+}