diff --git a/MaterialLib/SolidModels/CreepBGRa.cpp b/MaterialLib/SolidModels/CreepBGRa.cpp
index 78da877423ef33e63e2acc59b751cdf6f2612ad0..b59f6d85341dbb7558204c82f198616ba8243aa6 100644
--- a/MaterialLib/SolidModels/CreepBGRa.cpp
+++ b/MaterialLib/SolidModels/CreepBGRa.cpp
@@ -39,7 +39,7 @@ CreepBGRa<DisplacementDim>::integrateStress(
     KelvinVector const& sigma_prev,
     typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
     /*material_state_variables*/,
-    double const T)
+    double const T) const
 {
     using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
 
diff --git a/MaterialLib/SolidModels/CreepBGRa.h b/MaterialLib/SolidModels/CreepBGRa.h
index ee29d0e23c4d75249d39d2f34302534aea1f5772..affcbd6091f589ecd05d0c562c7e20e68262c516 100644
--- a/MaterialLib/SolidModels/CreepBGRa.h
+++ b/MaterialLib/SolidModels/CreepBGRa.h
@@ -55,7 +55,7 @@ public:
 
     std::unique_ptr<
         typename MechanicsBase<DisplacementDim>::MaterialStateVariables>
-    createMaterialStateVariables() override
+    createMaterialStateVariables() const override
     {
         return LinearElasticIsotropic<
             DisplacementDim>::createMaterialStateVariables();
@@ -75,17 +75,17 @@ public:
     {
     }
 
-    boost::optional<
-        std::tuple<KelvinVector, std::unique_ptr<typename MechanicsBase<
-                                     DisplacementDim>::MaterialStateVariables>,
-                   KelvinMatrix>>
+    boost::optional<std::tuple<KelvinVector,
+                               std::unique_ptr<typename MechanicsBase<
+                                   DisplacementDim>::MaterialStateVariables>,
+                               KelvinMatrix>>
     integrateStress(
         double const t, ProcessLib::SpatialPosition const& x, double const dt,
         KelvinVector const& eps_prev, KelvinVector const& eps,
         KelvinVector const& sigma_prev,
         typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
             material_state_variables,
-        double const T) override;
+        double const T) const override;
 
     ConstitutiveModel getConstitutiveModel() const override
     {
diff --git a/MaterialLib/SolidModels/Ehlers.cpp b/MaterialLib/SolidModels/Ehlers.cpp
index 23e0f9aef5a9ce77aeae9f40376f22891a01cfd2..680dcb246a3b803223a561ec6589eec011e2a207 100644
--- a/MaterialLib/SolidModels/Ehlers.cpp
+++ b/MaterialLib/SolidModels/Ehlers.cpp
@@ -473,14 +473,12 @@ boost::optional<std::tuple<typename SolidEhlers<DisplacementDim>::KelvinVector,
                                DisplacementDim>::MaterialStateVariables>,
                            typename SolidEhlers<DisplacementDim>::KelvinMatrix>>
 SolidEhlers<DisplacementDim>::integrateStress(
-    double const t,
-    ProcessLib::SpatialPosition const& x,
-    double const dt,
-    KelvinVector const& eps_prev,
-    KelvinVector const& eps,
+    double const t, ProcessLib::SpatialPosition const& x, double const dt,
+    KelvinVector const& eps_prev, KelvinVector const& eps,
     KelvinVector const& sigma_prev,
     typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
-        material_state_variables, double const /*T*/)
+        material_state_variables,
+    double const /*T*/) const
 {
     assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
                &material_state_variables) != nullptr);
diff --git a/MaterialLib/SolidModels/Ehlers.h b/MaterialLib/SolidModels/Ehlers.h
index 0be72d6005d89439e6a8db9adf1901415b1fe352..2958c1c94ed4150a3fb145091f9863c839d2538f 100644
--- a/MaterialLib/SolidModels/Ehlers.h
+++ b/MaterialLib/SolidModels/Ehlers.h
@@ -267,7 +267,7 @@ public:
 public:
     std::unique_ptr<
         typename MechanicsBase<DisplacementDim>::MaterialStateVariables>
-    createMaterialStateVariables() override
+    createMaterialStateVariables() const override
     {
         return std::make_unique<StateVariables<DisplacementDim>>();
     }
@@ -309,19 +309,17 @@ public:
         return (eps - eps_p.D - eps_p.V / 3 * identity2).dot(sigma) / 2;
     }
 
-    boost::optional<
-        std::tuple<KelvinVector, std::unique_ptr<typename MechanicsBase<
-                                     DisplacementDim>::MaterialStateVariables>,
-                   KelvinMatrix>>
+    boost::optional<std::tuple<KelvinVector,
+                               std::unique_ptr<typename MechanicsBase<
+                                   DisplacementDim>::MaterialStateVariables>,
+                               KelvinMatrix>>
     integrateStress(
-        double const t,
-        ProcessLib::SpatialPosition const& x,
-        double const dt,
-        KelvinVector const& eps_prev,
-        KelvinVector const& eps,
+        double const t, ProcessLib::SpatialPosition const& x, double const dt,
+        KelvinVector const& eps_prev, KelvinVector const& eps,
         KelvinVector const& sigma_prev,
         typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
-            material_state_variables, double const T) override;
+            material_state_variables,
+        double const T) const override;
 
     std::vector<typename MechanicsBase<DisplacementDim>::InternalVariable>
     getInternalVariables() const override;
diff --git a/MaterialLib/SolidModels/LinearElasticIsotropic.cpp b/MaterialLib/SolidModels/LinearElasticIsotropic.cpp
index 14d2c7bdcab8d92b5467f9131bebb43652e27799..e28133f1c425211b6c55bc8a9494df9d72286e1c 100644
--- a/MaterialLib/SolidModels/LinearElasticIsotropic.cpp
+++ b/MaterialLib/SolidModels/LinearElasticIsotropic.cpp
@@ -20,14 +20,12 @@ boost::optional<
                    DisplacementDim>::MaterialStateVariables>,
                typename LinearElasticIsotropic<DisplacementDim>::KelvinMatrix>>
 LinearElasticIsotropic<DisplacementDim>::integrateStress(
-    double const t,
-    ProcessLib::SpatialPosition const& x,
-    double const /*dt*/,
-    KelvinVector const& eps_prev,
-    KelvinVector const& eps,
+    double const t, ProcessLib::SpatialPosition const& x, double const /*dt*/,
+    KelvinVector const& eps_prev, KelvinVector const& eps,
     KelvinVector const& sigma_prev,
     typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
-        material_state_variables, double const T)
+        material_state_variables,
+    double const T) const
 {
     KelvinMatrix C = getElasticTensor(t, x, T);
 
diff --git a/MaterialLib/SolidModels/LinearElasticIsotropic.h b/MaterialLib/SolidModels/LinearElasticIsotropic.h
index 46bf4554f57416a38088cf7a693018f2d07112f3..2bd88582428d34a70668280d0816076ccb72ea17 100644
--- a/MaterialLib/SolidModels/LinearElasticIsotropic.h
+++ b/MaterialLib/SolidModels/LinearElasticIsotropic.h
@@ -75,7 +75,7 @@ public:
 
     std::unique_ptr<
         typename MechanicsBase<DisplacementDim>::MaterialStateVariables>
-    createMaterialStateVariables() override
+    createMaterialStateVariables() const override
     {
         return std::unique_ptr<
             typename MechanicsBase<DisplacementDim>::MaterialStateVariables>{
@@ -113,14 +113,12 @@ public:
                                    DisplacementDim>::MaterialStateVariables>,
                                KelvinMatrix>>
     integrateStress(
-        double const t,
-        ProcessLib::SpatialPosition const& x,
-        double const /*dt*/,
-        KelvinVector const& eps_prev,
-        KelvinVector const& eps,
-        KelvinVector const& sigma_prev,
+        double const t, ProcessLib::SpatialPosition const& x,
+        double const /*dt*/, KelvinVector const& eps_prev,
+        KelvinVector const& eps, KelvinVector const& sigma_prev,
         typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
-            material_state_variables, double const T) override;
+            material_state_variables,
+        double const T) const override;
 
     KelvinMatrix getElasticTensor(double const t,
                                   ProcessLib::SpatialPosition const& x,
diff --git a/MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.h b/MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.h
index 8edc0986ec478dd6ead0050245d25cf608a60c7e..003c68e70ffafcbb3a7413cf94d727e95ad5e825 100644
--- a/MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.h
+++ b/MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.h
@@ -41,7 +41,7 @@ public:
 
     std::unique_ptr<
         typename MechanicsBase<DisplacementDim>::MaterialStateVariables>
-    createMaterialStateVariables() override
+    createMaterialStateVariables() const override
     {
         return LinearElasticIsotropic<
             DisplacementDim>::createMaterialStateVariables();
@@ -52,14 +52,12 @@ public:
                                    DisplacementDim>::MaterialStateVariables>,
                                KelvinMatrix>>
     integrateStress(
-        double const t,
-        ProcessLib::SpatialPosition const& x,
-        double const dt,
-        KelvinVector const& eps_prev,
-        KelvinVector const& eps,
+        double const t, ProcessLib::SpatialPosition const& x, double const dt,
+        KelvinVector const& eps_prev, KelvinVector const& eps,
         KelvinVector const& sigma_prev,
         typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
-            material_state_variables, double const T) override
+            material_state_variables,
+        double const T) const override
     {
         return LinearElasticIsotropic<DisplacementDim>::integrateStress(
             t, x, dt, eps_prev, eps, sigma_prev, material_state_variables, T);
diff --git a/MaterialLib/SolidModels/Lubby2.cpp b/MaterialLib/SolidModels/Lubby2.cpp
index ac9160c34e8b749ca81b65c8c79266798ea7fde6..ab9024fd027e381e0fbe8d41183ba4d82ee88b2e 100644
--- a/MaterialLib/SolidModels/Lubby2.cpp
+++ b/MaterialLib/SolidModels/Lubby2.cpp
@@ -73,14 +73,12 @@ boost::optional<std::tuple<typename Lubby2<DisplacementDim>::KelvinVector,
                                DisplacementDim>::MaterialStateVariables>,
                            typename Lubby2<DisplacementDim>::KelvinMatrix>>
 Lubby2<DisplacementDim>::integrateStress(
-    double const t,
-    ProcessLib::SpatialPosition const& x,
-    double const dt,
-    KelvinVector const& /*eps_prev*/,
-    KelvinVector const& eps,
+    double const t, ProcessLib::SpatialPosition const& x, double const dt,
+    KelvinVector const& /*eps_prev*/, KelvinVector const& eps,
     KelvinVector const& /*sigma_prev*/,
     typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
-        material_state_variables, double const /*T*/)
+        material_state_variables,
+    double const /*T*/) const
 {
     using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
 
@@ -204,7 +202,7 @@ void Lubby2<DisplacementDim>::calculateResidualBurgers(
     KelvinVector& strain_Max_curr,
     const KelvinVector& strain_Max_t,
     ResidualVector& res,
-    detail::LocalLubby2Properties<DisplacementDim> const& properties)
+    detail::LocalLubby2Properties<DisplacementDim> const& properties) const
 {
     // calculate stress residual
     res.template segment<KelvinVectorSize>(0).noalias() =
@@ -231,7 +229,7 @@ void Lubby2<DisplacementDim>::calculateJacobianBurgers(
     double s_eff,
     const KelvinVector& sig_i,
     const KelvinVector& eps_K_i,
-    detail::LocalLubby2Properties<DisplacementDim> const& properties)
+    detail::LocalLubby2Properties<DisplacementDim> const& properties) const
 {
     Jac.setZero();
 
diff --git a/MaterialLib/SolidModels/Lubby2.h b/MaterialLib/SolidModels/Lubby2.h
index f8b6c5ed21ec395e18eb5c93cd694a26863fcdd2..e24b37d72b62c3dbbb18cc626104995f30109f6a 100644
--- a/MaterialLib/SolidModels/Lubby2.h
+++ b/MaterialLib/SolidModels/Lubby2.h
@@ -157,7 +157,7 @@ public:
 
     std::unique_ptr<
         typename MechanicsBase<DisplacementDim>::MaterialStateVariables>
-    createMaterialStateVariables() override
+    createMaterialStateVariables() const override
     {
         return std::unique_ptr<
             typename MechanicsBase<DisplacementDim>::MaterialStateVariables>{
@@ -224,14 +224,12 @@ public:
                                    DisplacementDim>::MaterialStateVariables>,
                                KelvinMatrix>>
     integrateStress(
-        double const t,
-        ProcessLib::SpatialPosition const& x,
-        double const dt,
-        KelvinVector const& eps_prev,
-        KelvinVector const& eps,
+        double const t, ProcessLib::SpatialPosition const& x, double const dt,
+        KelvinVector const& eps_prev, KelvinVector const& eps,
         KelvinVector const& sigma_prev,
         typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
-            material_state_variables, double const T) override;
+            material_state_variables,
+        double const T) const override;
 
 private:
     /// Calculates the 18x1 residual vector.
@@ -244,7 +242,7 @@ private:
         KelvinVector& strain_Max_curr,
         const KelvinVector& strain_Max_t,
         ResidualVector& res,
-        detail::LocalLubby2Properties<DisplacementDim> const& properties);
+        detail::LocalLubby2Properties<DisplacementDim> const& properties) const;
 
     /// Calculates the 18x18 Jacobian.
     void calculateJacobianBurgers(
@@ -255,7 +253,7 @@ private:
         double s_eff,
         const KelvinVector& sig_i,
         const KelvinVector& eps_K_i,
-        detail::LocalLubby2Properties<DisplacementDim> const& properties);
+        detail::LocalLubby2Properties<DisplacementDim> const& properties) const;
 
 private:
     NumLib::NewtonRaphsonSolverParameters const _nonlinear_solver_parameters;
diff --git a/MaterialLib/SolidModels/MechanicsBase.h b/MaterialLib/SolidModels/MechanicsBase.h
index a839aa3114088e442ddca7ca7bbff52361d7d10f..26f76b6622f1597aa2fc0c059525b252bf46cb19 100644
--- a/MaterialLib/SolidModels/MechanicsBase.h
+++ b/MaterialLib/SolidModels/MechanicsBase.h
@@ -66,7 +66,7 @@ struct MechanicsBase
     /// Polymorphic creator for MaterialStateVariables objects specific for a
     /// material model.
     virtual std::unique_ptr<MaterialStateVariables>
-    createMaterialStateVariables() = 0;
+    createMaterialStateVariables() const = 0;
 
     using KelvinVector =
         MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
@@ -77,9 +77,8 @@ struct MechanicsBase
     /// constitutive relation compute function.
     /// Returns nothing in case of errors in the computation if Newton
     /// iterations did not converge, for example.
-    boost::optional<std::tuple<KelvinVector,
-                               std::unique_ptr<MaterialStateVariables>,
-                               KelvinMatrix>>
+    boost::optional<std::tuple<
+        KelvinVector, std::unique_ptr<MaterialStateVariables>, KelvinMatrix>>
     integrateStress(double const t,
                     ProcessLib::SpatialPosition const& x,
                     double const dt,
@@ -87,7 +86,7 @@ struct MechanicsBase
                     Eigen::Matrix<double, Eigen::Dynamic, 1> const& eps,
                     Eigen::Matrix<double, Eigen::Dynamic, 1> const& sigma_prev,
                     MaterialStateVariables const& material_state_variables,
-                    double const T)
+                    double const T) const
     {
         // TODO Avoid copies of data:
         // Using MatrixBase<Derived> not possible because template functions
@@ -109,9 +108,8 @@ struct MechanicsBase
     /// wrapper function.
     /// Returns nothing in case of errors in the computation if Newton
     /// iterations did not converge, for example.
-    virtual boost::optional<std::tuple<KelvinVector,
-                                       std::unique_ptr<MaterialStateVariables>,
-                                       KelvinMatrix>>
+    virtual boost::optional<std::tuple<
+        KelvinVector, std::unique_ptr<MaterialStateVariables>, KelvinMatrix>>
     integrateStress(double const t,
                     ProcessLib::SpatialPosition const& x,
                     double const dt,
@@ -119,7 +117,7 @@ struct MechanicsBase
                     KelvinVector const& eps,
                     KelvinVector const& sigma_prev,
                     MaterialStateVariables const& material_state_variables,
-                    double const T) = 0;
+                    double const T) const = 0;
 
     /// Helper type for providing access to internal variables.
     struct InternalVariable