diff --git a/MaterialLib/SolidModels/LinearElasticIsotropic.h b/MaterialLib/SolidModels/LinearElasticIsotropic.h
index 063f1f074333b6a20311f7799fcde58a468e54d2..bf2d729c2fc88c9bda61a805bb85b86438da4dcd 100644
--- a/MaterialLib/SolidModels/LinearElasticIsotropic.h
+++ b/MaterialLib/SolidModels/LinearElasticIsotropic.h
@@ -79,7 +79,7 @@ public:
     {
     }
 
-    void computeConstitutiveRelation(
+    bool computeConstitutiveRelation(
         double const t,
         ProcessLib::SpatialPosition const& x,
         double const /*dt*/,
@@ -97,6 +97,7 @@ public:
         C.noalias() += 2 * _mp.mu(t, x) * KelvinMatrix::Identity();
 
         sigma.noalias() = sigma_prev + C * (eps - eps_prev);
+        return true;
     }
 
 private:
diff --git a/MaterialLib/SolidModels/Lubby2-impl.h b/MaterialLib/SolidModels/Lubby2-impl.h
index 392fd020edec5ff549a75cc6161e98a41335de43..909029fd66c591b202cbad97f6a0dd75e00600ec 100644
--- a/MaterialLib/SolidModels/Lubby2-impl.h
+++ b/MaterialLib/SolidModels/Lubby2-impl.h
@@ -17,7 +17,7 @@ namespace MaterialLib
 namespace Solids
 {
 template <int DisplacementDim>
-void Lubby2<DisplacementDim>::computeConstitutiveRelation(
+bool Lubby2<DisplacementDim>::computeConstitutiveRelation(
     double const t,
     ProcessLib::SpatialPosition const& x,
     double const dt,
@@ -115,7 +115,7 @@ void Lubby2<DisplacementDim>::computeConstitutiveRelation(
         auto const success_iterations = newton_solver.solve(K_loc);
 
         if (!success_iterations)
-            OGS_FATAL("The local Newton method did not succeed.");
+            return false;
 
         // If the Newton loop didn't run, the linear solver will not be
         // initialized.
@@ -142,6 +142,8 @@ void Lubby2<DisplacementDim>::computeConstitutiveRelation(
 
     auto const& P_sph = Invariants::spherical_projection;
     C.noalias() = state.GM * dzdE * P_dev + 3. * state.KM * P_sph;
+
+    return true;
 }
 
 template <int DisplacementDim>
diff --git a/MaterialLib/SolidModels/Lubby2.h b/MaterialLib/SolidModels/Lubby2.h
index e2d8d50672e27cf8c3506a8c87add55b6e31da19..2d41b6fc954a9f14924baf46b0c0d44fec0ccc7c 100644
--- a/MaterialLib/SolidModels/Lubby2.h
+++ b/MaterialLib/SolidModels/Lubby2.h
@@ -127,7 +127,7 @@ public:
     {
     }
 
-    void computeConstitutiveRelation(
+    bool computeConstitutiveRelation(
         double const t,
         ProcessLib::SpatialPosition const& x_position,
         double const dt,
diff --git a/MaterialLib/SolidModels/MechanicsBase.h b/MaterialLib/SolidModels/MechanicsBase.h
index 8875291d6fc29d4f827f0f6bc2583876e7c85246..7040c53aed8e73ed38802e9ce047267a8483f8cb 100644
--- a/MaterialLib/SolidModels/MechanicsBase.h
+++ b/MaterialLib/SolidModels/MechanicsBase.h
@@ -53,7 +53,9 @@ struct MechanicsBase
 
     /// Dynamic size Kelvin vector and matrix wrapper for the polymorphic
     /// constitutive relation compute function.
-    void computeConstitutiveRelation(
+    /// Returns false in case of errors in the computation if Newton iterations
+    /// did not converge, for example.
+    bool computeConstitutiveRelation(
         double const t,
         ProcessLib::SpatialPosition const& x,
         double const dt,
@@ -76,25 +78,29 @@ struct MechanicsBase
         KelvinVector sigma_{sigma};
         KelvinMatrix C_{C};
 
-        computeConstitutiveRelation(t,
-                                    x,
-                                    dt,
-                                    eps_prev_,
-                                    eps_,
-                                    sigma_prev_,
-                                    sigma_,
-                                    C_,
-                                    material_state_variables);
+        bool const result =
+            computeConstitutiveRelation(t,
+                                        x,
+                                        dt,
+                                        eps_prev_,
+                                        eps_,
+                                        sigma_prev_,
+                                        sigma_,
+                                        C_,
+                                        material_state_variables);
 
         sigma = sigma_;
         C = C_;
+        return result;
     }
 
     /// Computation of the constitutive relation for specific material model.
     /// 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.
-    virtual void computeConstitutiveRelation(
+    /// Returns false in case of errors in the computation if Newton iterations
+    /// did not converge, for example.
+    virtual bool computeConstitutiveRelation(
         double const t,
         ProcessLib::SpatialPosition const& x,
         double const dt,
diff --git a/ProcessLib/SmallDeformation/SmallDeformationFEM.h b/ProcessLib/SmallDeformation/SmallDeformationFEM.h
index eed4503422d0e27fb16f0804fba7dfb5044519cc..01cb3ab818c3fc424ab74420f152ad1622342bed 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationFEM.h
+++ b/ProcessLib/SmallDeformation/SmallDeformationFEM.h
@@ -234,9 +234,10 @@ public:
                 Eigen::Map<typename BMatricesType::NodalForceVectorType const>(
                     local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);
 
-            _ip_data[ip]._solid_material.computeConstitutiveRelation(
-                t, x_position, _process_data.dt, eps_prev, eps, sigma_prev,
-                sigma, C, material_state_variables);
+            if (!_ip_data[ip]._solid_material.computeConstitutiveRelation(
+                    t, x_position, _process_data.dt, eps_prev, eps, sigma_prev,
+                    sigma, C, material_state_variables))
+                OGS_FATAL("Computation of local constitutive relation failed.");
 
             local_b.noalias() -= B.transpose() * sigma * detJ * wp.getWeight();
             local_Jac.noalias() +=