diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h b/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
index b32a8d6b70f5c06d4ed612174b9470d9743c08d8..b5f1abbfbc7ca93b68af676feb1ff7f974f51c13 100644
--- a/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
@@ -62,31 +62,6 @@ struct IntegrationPointData final
         material_state_variables->pushBackState();
     }
 
-    static const int kelvin_vector_size =
-        KelvinVectorDimensions<DisplacementDim>::value;
-    using Invariants = MaterialLib::SolidModels::Invariants<kelvin_vector_size>;
-
-    typename BMatricesType::KelvinMatrixType updateConstitutiveRelation(
-        double const t,
-        SpatialPosition const& x_position,
-        double const dt,
-        double const linear_thermal_strain)
-    {
-        // assume isotropic thermal expansion
-        eps_m.noalias() = eps - linear_thermal_strain * Invariants::identity2;
-        auto&& solution = solid_material.integrateStress(
-            t, x_position, dt, eps_m_prev, eps_m, sigma_prev,
-            *material_state_variables);
-
-        if (!solution)
-            OGS_FATAL("Computation of local constitutive relation failed.");
-
-        KelvinMatrixType<DisplacementDim> C;
-        std::tie(sigma, material_state_variables, C) = std::move(*solution);
-
-        return C;
-    }
-
     EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
 };
 
@@ -240,9 +215,15 @@ public:
                 DisplacementDim, ShapeFunction::NPOINTS,
                 typename BMatricesType::BMatrixType>(dNdx, N, x_coord,
                                                      _is_axially_symmetric);
-            auto const& sigma = _ip_data[ip].sigma;
+
+            auto& sigma = _ip_data[ip].sigma;
+            auto const& sigma_prev = _ip_data[ip].sigma_prev;
             auto& eps = _ip_data[ip].eps;
 
+            auto& eps_m = _ip_data[ip].eps_m;
+            auto& eps_m_prev = _ip_data[ip].eps_m_prev;
+            auto& state = _ip_data[ip].material_state_variables;
+
             double const delta_T =
                 N.dot(T) - _process_data.reference_temperature;
             // calculate thermally induced strain
@@ -256,8 +237,21 @@ public:
             // displacement equation, displacement part
             //
             eps.noalias() = B * u;
-            auto C = _ip_data[ip].updateConstitutiveRelation(
-                t, x_position, dt, linear_thermal_strain);
+
+            using Invariants =
+                MaterialLib::SolidModels::Invariants<kelvin_vector_size>;
+
+            // assume isotropic thermal expansion
+            eps_m.noalias() =
+                eps - linear_thermal_strain * Invariants::identity2;
+            auto&& solution = _ip_data[ip].solid_material.integrateStress(
+                t, x_position, dt, eps_m_prev, eps_m, sigma_prev, *state);
+
+            if (!solution)
+                OGS_FATAL("Computation of local constitutive relation failed.");
+
+            KelvinMatrixType<DisplacementDim> C;
+            std::tie(sigma, state, C) = std::move(*solution);
 
             local_Jac
                 .template block<displacement_size, displacement_size>(
@@ -276,9 +270,6 @@ public:
                        i, i * displacement_size / DisplacementDim)
                     .noalias() = N;
 
-            using Invariants =
-                MaterialLib::SolidModels::Invariants<kelvin_vector_size>;
-
             // calculate real density
             auto const rho_sr =
                 _process_data.reference_solid_density(t, x_position)[0];