diff --git a/MaterialLib/SolidModels/Ehlers-impl.h b/MaterialLib/SolidModels/Ehlers-impl.h
index ec441ec48d6a828ad4adb93cdd4907ad0b11f9fb..e6e1eef452d418822c3d32622f00e527038ef9a6 100644
--- a/MaterialLib/SolidModels/Ehlers-impl.h
+++ b/MaterialLib/SolidModels/Ehlers-impl.h
@@ -655,10 +655,10 @@ newton(double const dt, MaterialProperties const& mp,
 }
 
 template <int DisplacementDim>
-std::tuple<typename SolidEhlers<DisplacementDim>::KelvinVector,
-           std::unique_ptr<
-               typename MechanicsBase<DisplacementDim>::MaterialStateVariables>,
-           typename SolidEhlers<DisplacementDim>::KelvinMatrix>
+boost::optional<std::tuple<typename SolidEhlers<DisplacementDim>::KelvinVector,
+                           std::unique_ptr<typename MechanicsBase<
+                               DisplacementDim>::MaterialStateVariables>,
+                           typename SolidEhlers<DisplacementDim>::KelvinMatrix>>
 SolidEhlers<DisplacementDim>::integrateStress(
     double const t,
     ProcessLib::SpatialPosition const& x,
@@ -727,7 +727,7 @@ SolidEhlers<DisplacementDim>::integrateStress(
                 state.eps_p_prev, s, sigma))
             std::tie(sigma, state.eps_p, linear_solver) = *solution;
         else
-                return {sigma, nullptr, tangentStiffness};
+            return {};
 
         if (_damage_properties)
         {
@@ -759,11 +759,12 @@ SolidEhlers<DisplacementDim>::integrateStress(
         sigma_final *= 1 - state.damage.value();
 
     return {
-        sigma_final,
-        std::unique_ptr<
-            typename MechanicsBase<DisplacementDim>::MaterialStateVariables>{
-            new StateVariables<DisplacementDim>{state}},
-        tangentStiffness};
+        {sigma_final,
+         std::unique_ptr<
+             typename MechanicsBase<DisplacementDim>::MaterialStateVariables>{
+             new StateVariables<DisplacementDim>{
+                 static_cast<StateVariables<DisplacementDim> const&>(state)}},
+         tangentStiffness}};
 }
 
 }  // namespace Ehlers
diff --git a/MaterialLib/SolidModels/Ehlers.h b/MaterialLib/SolidModels/Ehlers.h
index ff262fccd60ca48a142f0519cddf25816ce8db13..b0a712e42cccfa9d07442b177f3e811e46c266d7 100644
--- a/MaterialLib/SolidModels/Ehlers.h
+++ b/MaterialLib/SolidModels/Ehlers.h
@@ -223,9 +223,10 @@ public:
     {
     }
 
-    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,
diff --git a/MaterialLib/SolidModels/LinearElasticIsotropic.h b/MaterialLib/SolidModels/LinearElasticIsotropic.h
index 313b24c38b68eead815fec44a5a27e659df35c43..828230456bf1d20b9c81c8c5b07b2b8e7a7f2956 100644
--- a/MaterialLib/SolidModels/LinearElasticIsotropic.h
+++ b/MaterialLib/SolidModels/LinearElasticIsotropic.h
@@ -87,10 +87,10 @@ public:
     {
     }
 
-    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,
@@ -108,14 +108,13 @@ public:
 
         KelvinVector sigma = sigma_prev + C * (eps - eps_prev);
 
-        return std::make_tuple(
-            sigma,
-            std::unique_ptr<typename MechanicsBase<
-                DisplacementDim>::MaterialStateVariables>{
-                new MaterialStateVariables{
-                    static_cast<MaterialStateVariables const&>(
-                        material_state_variables)}},
-            C);
+        return {{sigma,
+                 std::unique_ptr<typename MechanicsBase<
+                     DisplacementDim>::MaterialStateVariables>{
+                     new MaterialStateVariables{
+                         static_cast<MaterialStateVariables const&>(
+                             material_state_variables)}},
+                 C}};
     }
 
 private:
diff --git a/MaterialLib/SolidModels/Lubby2-impl.h b/MaterialLib/SolidModels/Lubby2-impl.h
index 7ca157ca7d3e7433777e29828b2e2327120c81b6..fcd2739e1e2b1d31b94a81cea75b5adfc16f70e8 100644
--- a/MaterialLib/SolidModels/Lubby2-impl.h
+++ b/MaterialLib/SolidModels/Lubby2-impl.h
@@ -65,10 +65,10 @@ ProcessLib::KelvinMatrixType<DisplacementDim> tangentStiffnessA(
 };
 
 template <int DisplacementDim>
-std::tuple<typename Lubby2<DisplacementDim>::KelvinVector,
-           std::unique_ptr<
-               typename MechanicsBase<DisplacementDim>::MaterialStateVariables>,
-           typename Lubby2<DisplacementDim>::KelvinMatrix>
+boost::optional<std::tuple<typename Lubby2<DisplacementDim>::KelvinVector,
+                           std::unique_ptr<typename MechanicsBase<
+                               DisplacementDim>::MaterialStateVariables>,
+                           typename Lubby2<DisplacementDim>::KelvinMatrix>>
 Lubby2<DisplacementDim>::integrateStress(
     double const t,
     ProcessLib::SpatialPosition const& x,
@@ -164,7 +164,7 @@ Lubby2<DisplacementDim>::integrateStress(
         auto const success_iterations = newton_solver.solve(K_loc);
 
         if (!success_iterations)
-            return std::make_tuple(sigma_prev, nullptr, KelvinMatrix::Zero());
+            return {};
 
         // If the Newton loop didn't run, the linear solver will not be
         // initialized.
@@ -180,13 +180,13 @@ Lubby2<DisplacementDim>::integrateStress(
 
     // Hydrostatic part for the stress and the tangent.
     double const eps_i_trace = Invariants::trace(eps);
-    return std::make_tuple(
-        local_lubby2_properties.GM0 * sigd_j +
-            local_lubby2_properties.KM0 * eps_i_trace * Invariants::identity2,
-        std::unique_ptr<
-            typename MechanicsBase<DisplacementDim>::MaterialStateVariables>{
-            new MaterialStateVariables{state}},
-        C);
+    return {
+        {local_lubby2_properties.GM0 * sigd_j +
+             local_lubby2_properties.KM0 * eps_i_trace * Invariants::identity2,
+         std::unique_ptr<
+             typename MechanicsBase<DisplacementDim>::MaterialStateVariables>{
+             new MaterialStateVariables{state}},
+         C}};
 }
 
 template <int DisplacementDim>
diff --git a/MaterialLib/SolidModels/Lubby2.h b/MaterialLib/SolidModels/Lubby2.h
index 0747951c04723cd5e97524ce378a12694dfd3d8e..5c3516b7bfefc11294547b4b49228c91f312f83b 100644
--- a/MaterialLib/SolidModels/Lubby2.h
+++ b/MaterialLib/SolidModels/Lubby2.h
@@ -181,10 +181,10 @@ public:
     {
     }
 
-    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,
diff --git a/MaterialLib/SolidModels/MechanicsBase.h b/MaterialLib/SolidModels/MechanicsBase.h
index edcf6072952a83b22ee46f089e9bd8ebe4b33126..6f3135144976136beaeee78bd984e73d6b2cba36 100644
--- a/MaterialLib/SolidModels/MechanicsBase.h
+++ b/MaterialLib/SolidModels/MechanicsBase.h
@@ -9,6 +9,7 @@
 
 #pragma once
 
+#include <boost/optional.hpp>
 #include <memory>
 #include <tuple>
 
@@ -60,11 +61,11 @@ struct MechanicsBase
 
     /// Dynamic size Kelvin vector and matrix wrapper for the polymorphic
     /// constitutive relation compute function.
-    /// Returns false in case of errors in the computation if Newton iterations
-    /// did not converge, for example.
-    std::tuple<KelvinVector,
-               std::unique_ptr<MaterialStateVariables>,
-               KelvinMatrix>
+    /// 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>>
     integrateStress(double const t,
                     ProcessLib::SpatialPosition const& x,
                     double const dt,
@@ -90,11 +91,11 @@ struct MechanicsBase
     /// This should be implemented in the derived model. Fixed Kelvin vector and
     /// matrix size version; for dynamic size arguments there is an overloaded
     /// wrapper function.
-    /// Returns false in case of errors in the computation if Newton iterations
-    /// did not converge, for example.
-    virtual std::tuple<KelvinVector,
-                       std::unique_ptr<MaterialStateVariables>,
-                       KelvinMatrix>
+    /// 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>>
     integrateStress(double const t,
                     ProcessLib::SpatialPosition const& x,
                     double const dt,
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
index 3170269f1dcfda142221fd1f4b621b078683bbc5..1868619a5d489e5f9f872d0b4a018f0436bf5598 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
@@ -90,17 +90,15 @@ struct IntegrationPointData final
     {
         eps.noalias() = b_matrices * u;
 
-        KelvinMatrixType<DisplacementDim> C;
-        std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
-            DisplacementDim>::MaterialStateVariables>
-            new_state;
-        std::tie(sigma_eff, new_state, C) = solid_material.integrateStress(
+        auto&& solution = solid_material.integrateStress(
             t, x_position, dt, eps_prev, eps, sigma_eff_prev,
             *material_state_variables);
 
-        if (!new_state)
+        if (!solution)
             OGS_FATAL("Computation of local constitutive relation failed.");
-        material_state_variables = std::move(new_state);
+
+        KelvinMatrixType<DisplacementDim> C;
+        std::tie(sigma_eff, material_state_variables, C) = std::move(*solution);
 
         return C;
     }
diff --git a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix-impl.h b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix-impl.h
index ad5d9e2eb23d52a28b7943178a38db6d2e9fc131..9921426409e8d994c93a6d219e54cd4c06c537b2 100644
--- a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix-impl.h
+++ b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix-impl.h
@@ -224,18 +224,15 @@ assembleBlockMatricesWithJacobian(
 
         eps.noalias() = B * u;
 
-        KelvinMatrixType<GlobalDim> C;
-        std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
-            GlobalDim>::MaterialStateVariables>
-            new_state;
-        std::tie(sigma_eff, new_state, C) =
-            _ip_data[ip].solid_material.integrateStress(
-                t, x_position, _process_data.dt, eps_prev, eps, sigma_eff_prev,
-                *state);
-
-        if (!new_state)
+        auto&& solution = _ip_data[ip].solid_material.integrateStress(
+            t, x_position, _process_data.dt, eps_prev, eps, sigma_eff_prev,
+            *state);
+
+        if (!solution)
             OGS_FATAL("Computation of local constitutive relation failed.");
-        state = std::move(new_state);
+
+        KelvinMatrixType<GlobalDim> C;
+        std::tie(sigma_eff, state, C) = std::move(*solution);
 
         q.noalias() = - k_over_mu * (dNdx_p * p + rho_fr * gravity_vec);
 
@@ -337,17 +334,15 @@ computeSecondaryVariableConcreteWithBlockVectors(
 
         eps.noalias() = B * u;
 
-        std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
-            GlobalDim>::MaterialStateVariables>
-            new_state;
-        std::tie(sigma_eff, new_state, ip_data.C) =
-            _ip_data[ip].solid_material.integrateStress(
-                t, x_position, _process_data.dt, eps_prev, eps, sigma_eff_prev,
-                *state);
+        auto&& solution = _ip_data[ip].solid_material.integrateStress(
+            t, x_position, _process_data.dt, eps_prev, eps, sigma_eff_prev,
+            *state);
 
-        if (!new_state)
+        if (!solution)
             OGS_FATAL("Computation of local constitutive relation failed.");
-        state = std::move(new_state);
+
+        KelvinMatrixType<GlobalDim> C;
+        std::tie(sigma_eff, state, C) = std::move(*solution);
 
         q.noalias() = - k_over_mu * (dNdx_p * p + rho_fr * gravity_vec);
     }
diff --git a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrix-impl.h b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrix-impl.h
index e34f3933147626478a15f75106302295bad10a3d..3418b593a9969cb3c903dd5a82b0cbbb1a8fad59 100644
--- a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrix-impl.h
+++ b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrix-impl.h
@@ -144,18 +144,14 @@ assembleWithJacobian(
             Eigen::Map<typename BMatricesType::NodalForceVectorType const>(
                 local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);
 
-        KelvinMatrixType<DisplacementDim> C;
-        std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
-            DisplacementDim>::MaterialStateVariables>
-            new_state;
-        std::tie(sigma, new_state, C) =
-            _ip_data[ip]._solid_material.integrateStress(
-                t, x_position, _process_data.dt, eps_prev, eps, sigma_prev,
-                *state);
-
-        if (!new_state)
+        auto&& solution = _ip_data[ip]._solid_material.integrateStress(
+            t, x_position, _process_data.dt, eps_prev, eps, sigma_prev, *state);
+
+        if (!solution)
             OGS_FATAL("Computation of local constitutive relation failed.");
-        state = std::move(new_state);
+
+        KelvinMatrixType<DisplacementDim> C;
+        std::tie(sigma, state, C) = std::move(*solution);
 
         local_b.noalias() -=
             B.transpose() * sigma * detJ * wp.getWeight() * integralMeasure;
diff --git a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture-impl.h b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture-impl.h
index acaa2047718ada6915afbc487e32b12b29984f37..84b70ec46453964f98a6e7ab64a2e7cbe6f852a1 100644
--- a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture-impl.h
+++ b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture-impl.h
@@ -212,18 +212,14 @@ assembleWithJacobian(
 
         eps.noalias() = B * nodal_total_u;
 
-        KelvinMatrixType<DisplacementDim> C;
-        std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
-            DisplacementDim>::MaterialStateVariables>
-            new_state;
-        std::tie(sigma, new_state, C) =
-            _ip_data[ip]._solid_material.integrateStress(
-                t, x_position, _process_data.dt, eps_prev, eps, sigma_prev,
-                *state);
-
-        if (!new_state)
+        auto&& solution = _ip_data[ip]._solid_material.integrateStress(
+            t, x_position, _process_data.dt, eps_prev, eps, sigma_prev, *state);
+
+        if (!solution)
             OGS_FATAL("Computation of local constitutive relation failed.");
-        state = std::move(new_state);
+
+        KelvinMatrixType<DisplacementDim> C;
+        std::tie(sigma, state, C) = std::move(*solution);
 
         // r_u = B^T * Sigma = B^T * C * B * (u+phi*[u])
         // r_[u] = (phi*B)^T * Sigma = (phi*B)^T * C * B * (u+phi*[u])
diff --git a/ProcessLib/SmallDeformation/SmallDeformationFEM.h b/ProcessLib/SmallDeformation/SmallDeformationFEM.h
index dd857394f9e999d6a5892875a6f5ea70679fb9e1..e8493f15fce106d4ba992576c66eccd72ae319d1 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationFEM.h
+++ b/ProcessLib/SmallDeformation/SmallDeformationFEM.h
@@ -254,18 +254,15 @@ public:
                 Eigen::Map<typename BMatricesType::NodalForceVectorType const>(
                     local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);
 
-            KelvinMatrixType<DisplacementDim> C;
-            std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
-                DisplacementDim>::MaterialStateVariables>
-                new_state;
-            std::tie(sigma, new_state, C) =
-                _ip_data[ip].solid_material.integrateStress(
-                    t, x_position, _process_data.dt, eps_prev, eps, sigma_prev,
-                    *state);
-
-            if (!new_state)
+            auto&& solution = _ip_data[ip].solid_material.integrateStress(
+                t, x_position, _process_data.dt, eps_prev, eps, sigma_prev,
+                *state);
+
+            if (!solution)
                 OGS_FATAL("Computation of local constitutive relation failed.");
-            state = std::move(new_state);
+
+            KelvinMatrixType<DisplacementDim> C;
+            std::tie(sigma, state, C) = std::move(*solution);
 
             local_b.noalias() -= B.transpose() * sigma * w;
             local_Jac.noalias() += B.transpose() * C * B * w;