From a70a36a48c9d3b45a26eb06602c9a7ff362237a8 Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Tue, 23 May 2017 14:07:43 +0200
Subject: [PATCH] Return boost::optional from integrateStress().

---
 MaterialLib/SolidModels/Ehlers-impl.h         | 21 +++++------
 MaterialLib/SolidModels/Ehlers.h              |  7 ++--
 .../SolidModels/LinearElasticIsotropic.h      | 23 ++++++------
 MaterialLib/SolidModels/Lubby2-impl.h         | 24 ++++++-------
 MaterialLib/SolidModels/Lubby2.h              |  8 ++---
 MaterialLib/SolidModels/MechanicsBase.h       | 21 +++++------
 ProcessLib/HydroMechanics/HydroMechanicsFEM.h | 12 +++----
 .../HydroMechanicsLocalAssemblerMatrix-impl.h | 35 ++++++++-----------
 ...mallDeformationLocalAssemblerMatrix-impl.h | 18 ++++------
 ...ionLocalAssemblerMatrixNearFracture-impl.h | 18 ++++------
 .../SmallDeformation/SmallDeformationFEM.h    | 19 +++++-----
 11 files changed, 95 insertions(+), 111 deletions(-)

diff --git a/MaterialLib/SolidModels/Ehlers-impl.h b/MaterialLib/SolidModels/Ehlers-impl.h
index ec441ec48d6..e6e1eef452d 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 ff262fccd60..b0a712e42cc 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 313b24c38b6..828230456bf 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 7ca157ca7d3..fcd2739e1e2 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 0747951c047..5c3516b7bfe 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 edcf6072952..6f313514497 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 3170269f1dc..1868619a5d4 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 ad5d9e2eb23..9921426409e 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 e34f3933147..3418b593a99 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 acaa2047718..84b70ec4645 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 dd857394f9e..e8493f15fce 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;
-- 
GitLab