[MatL] Extract calculateDegradedStress().

The multiple inheritance for the phase field extension
is not actually needed and can be completely avoided.

MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.cpp

-/**
- * \copyright
- * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#include "LinearElasticIsotropicPhaseField.h"
-
-namespace MaterialLib
-{
-namespace Solids
-{
-
-template class LinearElasticIsotropicPhaseField<2>;
-template class LinearElasticIsotropicPhaseField<3>;
-
-}   // namespace Solids
-}   // namespace MaterialLib

MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.h

@@ -11,135 +11,84 @@
 
 #include "MathLib/KelvinVector.h"
 
-#include "LinearElasticIsotropic.h"
-#include "PhaseFieldExtension.h"
-
 namespace MaterialLib
 {
 namespace Solids
 {
+namespace Phasefield
+{
+/** Decompose the stiffness into tensile and compressive part.
+ * Judging by the physical observations, compression perpendicular
+ * to a crack does not cause crack propagation. Thus,
+ * the phase-field parameter is only involved into the tensile part
+ * to degrade the elastic strain energy.
+ */
 template <int DisplacementDim>
-class LinearElasticIsotropicPhaseField final
-    : public LinearElasticIsotropic<DisplacementDim>,
-      public PhaseFieldExtension<DisplacementDim>
+std::tuple<MathLib::KelvinVector::KelvinVectorType<
+               DisplacementDim> /* sigma_real */,
+           MathLib::KelvinVector::KelvinVectorType<
+               DisplacementDim> /* sigma_tensile */,
+           MathLib::KelvinVector::KelvinMatrixType<
+               DisplacementDim> /* C_tensile */,
+           MathLib::KelvinVector::KelvinMatrixType<
+               DisplacementDim> /* C_compressive */,
+           double /* strain_energy_tensile */,
+           double /* elastic_energy */
+           >
+calculateDegradedStress(
+    double const degradation,
+    double const bulk_modulus,
+    double const mu,
+    MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const& eps)
 {
-public:
     static int const KelvinVectorSize =
         MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
     using KelvinVector =
         MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
     using KelvinMatrix =
         MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>;
+    using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
+    // calculation of deviatoric parts
+    auto const& P_dev = Invariants::deviatoric_projection;
+    KelvinVector const epsd_curr = P_dev * eps;
 
-    explicit LinearElasticIsotropicPhaseField(
-        typename LinearElasticIsotropic<DisplacementDim>::MaterialProperties&&
-            material_properties)
-        : LinearElasticIsotropic<DisplacementDim>(
-              std::move(material_properties))
-    {
-    }
-
-    std::unique_ptr<
-        typename MechanicsBase<DisplacementDim>::MaterialStateVariables>
-    createMaterialStateVariables() const override
-    {
-        return LinearElasticIsotropic<
-            DisplacementDim>::createMaterialStateVariables();
-    }
-
-    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) const override
-    {
-        return LinearElasticIsotropic<DisplacementDim>::integrateStress(
-            t, x, dt, eps_prev, eps, sigma_prev, material_state_variables, T);
-    }
+    // Hydrostatic part for the stress and the tangent.
+    double const eps_curr_trace = Invariants::trace(eps);
 
-    double computeFreeEnergyDensity(
-        double const t,
-        ProcessLib::SpatialPosition const& x,
-        double const dt,
-        KelvinVector const& eps,
-        KelvinVector const& sigma,
-        typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
-            material_state_variables) const override
-    {
-        return LinearElasticIsotropic<DisplacementDim>::
-            computeFreeEnergyDensity(
-                t, x, dt, eps, sigma, material_state_variables);
-    }
+    KelvinMatrix C_tensile = KelvinMatrix::Zero();
+    KelvinMatrix C_compressive = KelvinMatrix::Zero();
 
-    /** Decompose the stiffness into tensile and compressive part.
-     * Judging by the physical observations, compression perpendicular
-     * to a crack does not cause crack propagation. Thus,
-     * the phase-field parameter is only involved into the tensile part
-     * to degrade the elastic strain energy.
-     */
-    bool calculateDegradedStress(double const t,
-                                 ProcessLib::SpatialPosition const& x,
-                                 KelvinVector const& eps,
-                                 double& strain_energy_tensile,
-                                 KelvinVector& sigma_tensile,
-                                 KelvinVector& sigma_compressive,
-                                 KelvinMatrix& C_tensile,
-                                 KelvinMatrix& C_compressive,
-                                 KelvinVector& sigma_real,
-                                 double const degradation,
-                                 double& elastic_energy) const override
+    if (eps_curr_trace >= 0)
     {
-        using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
-        // calculation of deviatoric parts
-        auto const& P_dev = Invariants::deviatoric_projection;
-        KelvinVector const epsd_curr = P_dev * eps;
-
-        // Hydrostatic part for the stress and the tangent.
-        double const eps_curr_trace = Invariants::trace(eps);
-
-        auto const& K =
-            LinearElasticIsotropic<DisplacementDim>::_mp.bulk_modulus(t, x);
-        auto const& mu = LinearElasticIsotropic<DisplacementDim>::_mp.mu(t, x);
-
-        C_tensile = KelvinMatrix::Zero();
-        C_compressive = KelvinMatrix::Zero();
-
-        if (eps_curr_trace >= 0)
-        {
-            strain_energy_tensile = K / 2 * eps_curr_trace * eps_curr_trace +
-                                    mu * epsd_curr.transpose() * epsd_curr;
-            sigma_tensile.noalias() =
-                K * eps_curr_trace * Invariants::identity2 + 2 * mu * epsd_curr;
-            sigma_compressive.noalias() = KelvinVector::Zero();
-            C_tensile.template topLeftCorner<3, 3>().setConstant(K);
-            C_tensile.noalias() += 2 * mu * P_dev * KelvinMatrix::Identity();
-            elastic_energy = degradation * strain_energy_tensile;
-        }
-        else
-        {
-            strain_energy_tensile = mu * epsd_curr.transpose() * epsd_curr;
-            sigma_tensile.noalias() = 2 * mu * epsd_curr;
-            sigma_compressive.noalias() =
-                K * eps_curr_trace * Invariants::identity2;
-            C_tensile.noalias() = 2 * mu * P_dev * KelvinMatrix::Identity();
-            C_compressive.template topLeftCorner<3, 3>().setConstant(K);
-            elastic_energy = K / 2 * eps_curr_trace * eps_curr_trace +
-                             mu * epsd_curr.transpose() * epsd_curr;
-        }
-
-        sigma_real.noalias() = degradation * sigma_tensile + sigma_compressive;
-        return true;
+        double const strain_energy_tensile =
+            bulk_modulus / 2 * eps_curr_trace * eps_curr_trace +
+            mu * epsd_curr.transpose() * epsd_curr;
+        KelvinVector const sigma_tensile =
+            bulk_modulus * eps_curr_trace * Invariants::identity2 +
+            2 * mu * epsd_curr;
+        C_tensile.template topLeftCorner<3, 3>().setConstant(bulk_modulus);
+        C_tensile.noalias() += 2 * mu * P_dev * KelvinMatrix::Identity();
+        double const elastic_energy = degradation * strain_energy_tensile;
+        KelvinVector const sigma_real =
+            degradation * sigma_tensile;  // + sigma_compressive, which is zero;
+        return std::make_tuple(sigma_real, sigma_tensile, C_tensile,
+                               C_compressive, strain_energy_tensile,
+                               elastic_energy);
     }
-};
-
-extern template class LinearElasticIsotropicPhaseField<2>;
-extern template class LinearElasticIsotropicPhaseField<3>;
-
+    double const strain_energy_tensile = mu * epsd_curr.transpose() * epsd_curr;
+    KelvinVector const sigma_tensile = 2 * mu * epsd_curr;
+    KelvinVector const sigma_compressive =
+        bulk_modulus * eps_curr_trace * Invariants::identity2;
+    C_tensile.noalias() = 2 * mu * P_dev * KelvinMatrix::Identity();
+    C_compressive.template topLeftCorner<3, 3>().setConstant(bulk_modulus);
+    double const elastic_energy =
+        bulk_modulus / 2 * eps_curr_trace * eps_curr_trace +
+        mu * epsd_curr.transpose() * epsd_curr;
+    KelvinVector const sigma_real =
+        degradation * sigma_tensile + sigma_compressive;
+    return std::make_tuple(sigma_real, sigma_tensile, C_tensile, C_compressive,
+                           strain_energy_tensile, elastic_energy);
+}
+}  // namespace Phasefield
 }  // namespace Solids
 }  // namespace MaterialLib

MaterialLib/SolidModels/PhaseFieldExtension.h

-/**
- * \file
- * \copyright
- * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#pragma once
-
-#include "MechanicsBase.h"
-
-namespace MaterialLib
-{
-namespace Solids
-{
-template <int DisplacementDim>
-struct PhaseFieldExtension : public MechanicsBase<DisplacementDim>
-{
-    using KelvinVector =
-        MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
-    using KelvinMatrix =
-        MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>;
-    virtual bool calculateDegradedStress(double const t,
-                                         ProcessLib::SpatialPosition const& x,
-                                         KelvinVector const& eps,
-                                         double& strain_energy_tensile,
-                                         KelvinVector& sigma_tensile,
-                                         KelvinVector& sigma_compressive,
-                                         KelvinMatrix& C_tensile,
-                                         KelvinMatrix& C_compressive,
-                                         KelvinVector& sigma_real,
-                                         double const degradation,
-                                         double& elastic_energy) const = 0;
-
-    /// Dynamic size Kelvin vector and matrix wrapper for the polymorphic
-    /// constitutive relation compute function.
-    bool calculateDegradedStress(
-        double const t,
-        ProcessLib::SpatialPosition const& x,
-        Eigen::Matrix<double, Eigen::Dynamic, 1> const& eps,
-        double& strain_energy_tensile,
-        Eigen::Matrix<double, Eigen::Dynamic, 1>& sigma_tensile,
-        Eigen::Matrix<double, Eigen::Dynamic, 1>& sigma_compressive,
-        Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
-            C_tensile,
-        Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
-            C_compressive,
-        Eigen::Matrix<double, Eigen::Dynamic, 1>& sigma_real,
-        double const degradation,
-        double& elastic_energy) const
-    {
-        // TODO Avoid copies of data:
-        // Using MatrixBase<Derived> not possible because template functions
-        // cannot be virtual. Maybe there is a workaround for this.  Using
-        // Map<Matrix<double, ...>> makes the interface (for the material model
-        // implementation) unnecessary difficult.
-        KelvinVector const eps_{eps};
-        KelvinVector sigma_tensile_{sigma_tensile};
-        KelvinVector sigma_compressive_{sigma_compressive};
-        KelvinMatrix C_tensile_{C_tensile};
-        KelvinMatrix C_compressive_{C_compressive};
-        KelvinVector sigma_real_{sigma_real};
-
-        bool const result = calculateDegradedStress(t,
-                                                    x,
-                                                    eps_,
-                                                    strain_energy_tensile,
-                                                    sigma_tensile_,
-                                                    sigma_compressive_,
-                                                    C_tensile_,
-                                                    C_compressive_,
-                                                    sigma_real_,
-                                                    degradation,
-                                                    elastic_energy);
-
-        sigma_tensile = sigma_tensile_;
-        sigma_compressive = sigma_compressive_;
-        C_tensile = C_tensile_;
-        C_compressive = C_compressive_;
-        sigma_real = sigma_real_;
-        return result;
-    }
-};
-
-}  // namespace Solids
-}  // namespace MaterialLib

ProcessLib/PhaseField/PhaseFieldFEM.h

@@ -70,20 +70,28 @@ struct IntegrationPointData final
 
     template <typename DisplacementVectorType>
     void updateConstitutiveRelation(double const t,
-                                    SpatialPosition const& x_position,
+                                    SpatialPosition const& x,
                                     double const /*dt*/,
                                     DisplacementVectorType const& /*u*/,
                                     double const degradation)
     {
-        static_cast<
-            MaterialLib::Solids::PhaseFieldExtension<DisplacementDim> const&>(
-            solid_material)
-            .calculateDegradedStress(t, x_position, eps, strain_energy_tensile,
-                                     sigma_tensile, sigma_compressive,
-                                     C_tensile, C_compressive, sigma,
-                                     degradation, elastic_energy);
+        auto linear_elastic_mp =
+            static_cast<MaterialLib::Solids::LinearElasticIsotropic<
+                DisplacementDim> const&>(solid_material)
+                .getMaterialProperties();
+
+        auto const bulk_modulus = linear_elastic_mp.bulk_modulus(t, x);
+        auto const mu = linear_elastic_mp.mu(t, x);
+
+        std::tie(sigma, sigma_tensile, C_tensile, C_compressive,
+                 strain_energy_tensile, elastic_energy) =
+            MaterialLib::Solids::Phasefield::calculateDegradedStress<
+                DisplacementDim>(degradation, bulk_modulus, mu, eps);
     }
     EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
+
+    static constexpr int kelvin_vector_size =
+        MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
 };
 
 /// Used by for extrapolation of the integration point values. It is ordered