[GeoMechanicsApplication] Make tension cut off truly optional for Mohr Coulomb model

applications/GeoMechanicsApplication/custom_constitutive/README.md

@@ -94,7 +94,7 @@ where:
 This criterion represents a linear envelope in the Mohr stress space, approximating the shear strength of a material under different stress states. For $`F_{MC} \le 0`$ we have elastic behavior, while $`F_{MC} \gt 0`$ corresponds to plastic behavior.
 
 
-Since the Mohr-Coulomb criterion primarily accounts for shear failure, it does not limit tensile stresses. In geomechanical applications, materials such as rocks and soils have a limited tensile strength . A tension cutoff is imposed as:
+Since the Mohr-Coulomb criterion primarily accounts for shear failure, it does not limit tensile stresses. In geomechanical applications, materials such as rocks and soils have a limited tensile strength . Therefore, the Mohr-Coulomb model in the geomechanics application supports an optional tension cutoff, which is imposed as:
 
 ```math
     F_{tc}(\sigma) = \sigma_1 - t_c = 0
@@ -106,6 +106,8 @@ The combination of these two, yields to the following figure:
 
 <img src="documentation_data/mohr-coulomb-with-tension-cutoff-zones.svg" alt="Mohr-Coulomb with tension cutoff" title="Mohr-Coulomb with tension cutoff" width="800">
 
+If the tension cutoff is not active, the figure can be simplified. In that case, the corner return coincides with the apex, resulting in the regular failure zone and the adjacent apex return zone, bound by the $\frac{\partial G_{MC}}{\partial \sigma}$ line and the horizontal $\sigma$ axis.
+
 
 ### 2.2 Implementation Mohr-Coulomb with tension cutoff
 
@@ -123,16 +125,16 @@ To incorporate the Mohr-Coulomb model with tension cutoff in numerical simulatio
 
 4. Evaluate the condition and mapping  
         4.1. If the trial stress falls in the elastic zone, it stays unchanged. No mapping is applied.  
-        4.2. If the trial stress falls in the tensile apex return zone. The trial stress then needs to be mapped back to the apex.  
-        4.3. If the trial stress falls in the tension cutoff zone. The trial stress then needs to be mapped back to the tension cutoff yield surface.  
-        4.4. If it falls in the tensile corner return zone, then it needs to be mapped to the corner point.  
+        4.2. If the trial stress falls in the (tensile) apex return zone. The trial stress then needs to be mapped back to the apex.  
+        4.3. If the trial stress falls in the tension cutoff zone. The trial stress then needs to be mapped back to the tension cutoff yield surface (optional).  
+        4.4. If it falls in the tensile corner return zone, then it needs to be mapped to the corner point (optional).  
         4.5. In the case of regular failure zone, then it is mapped back to the Mohr-Coulomb yield surface along the normal direction of flow function. The flow function is defined by  
 ```math
     G_{MC}(\sigma) = \frac{\sigma_1 - \sigma_3}{2} + \frac{\sigma_1 + \sigma_3}{2} \sin⁡{\psi}
 ```
   where $`\psi`$ is the dilatancy angle.
 
-5. If after mapping, the condidition $`\sigma_1 \ge \sigma_2 \ge \sigma_3`$ is not valid, average the principal stresses of stage 2 and the direction of the mapping and map the principal stresses again.
+5. If after mapping, the condition $`\sigma_1 \ge \sigma_2 \ge \sigma_3`$ is not valid, average the principal stresses of stage 2 and the direction of the mapping and map the principal stresses again.
 
 - If $`\sigma_1 \le \sigma_2`$ set:
 ```math

applications/GeoMechanicsApplication/custom_constitutive/coulomb_with_tension_cut_off_impl.h → applications/GeoMechanicsApplication/custom_constitutive/coulomb_impl.h

@@ -15,6 +15,8 @@
 
 #pragma once
 
+#include <optional>
+
 #include "coulomb_yield_surface.h"
 #include "custom_constitutive/principal_stresses.hpp"
 #include "geo_mechanics_application_constants.h"
@@ -31,11 +33,11 @@ namespace Geo
 class SigmaTau;
 } // namespace Geo
 
-class CoulombWithTensionCutOffImpl
+class CoulombImpl
 {
 public:
-    CoulombWithTensionCutOffImpl() = default;
-    explicit CoulombWithTensionCutOffImpl(const Properties& rMaterialProperties);
+    CoulombImpl() = default;
+    explicit CoulombImpl(const Properties& rMaterialProperties);
 
     [[nodiscard]] bool IsAdmissibleStressState(const Geo::SigmaTau& rTrialTraction);
     [[nodiscard]] bool IsAdmissibleStressState(const Geo::PrincipalStresses& rTrialPrincipalStresses);
@@ -52,12 +54,12 @@ class CoulombWithTensionCutOffImpl
     [[nodiscard]] PlasticityStatus GetPlasticityStatus() const;
 
 private:
-    CoulombYieldSurface mCoulombYieldSurface;
-    TensionCutoff       mTensionCutOff;
-    double              mSavedKappaOfCoulombYieldSurface{0.0};
-    double              mAbsoluteYieldFunctionValueTolerance{1.0e-8};
-    std::size_t         mMaxNumberOfPlasticIterations{100};
-    PlasticityStatus    mPlasticityStatus{PlasticityStatus::ELASTIC};
+    CoulombYieldSurface          mCoulombYieldSurface;
+    std::optional<TensionCutoff> mTensionCutOff;
+    double                       mSavedKappaOfCoulombYieldSurface{0.0};
+    double                       mAbsoluteYieldFunctionValueTolerance{1.0e-8};
+    std::size_t                  mMaxNumberOfPlasticIterations{100};
+    PlasticityStatus             mPlasticityStatus{PlasticityStatus::ELASTIC};
 
     template <typename StressStateType>
     [[nodiscard]] bool IsAdmissibleStressState(const StressStateType& rTrialStressState);

applications/GeoMechanicsApplication/custom_constitutive/interface_coulomb_with_tension_cut_off.cpp → applications/GeoMechanicsApplication/custom_constitutive/interface_coulomb_law.cpp

@@ -13,7 +13,7 @@
 //
 
 // Application includes
-#include "custom_constitutive/interface_coulomb_with_tension_cut_off.h"
+#include "custom_constitutive/interface_coulomb_law.h"
 #include "custom_constitutive/constitutive_law_dimension.h"
 #include "custom_constitutive/sigma_tau.hpp"
 #include "custom_utilities/check_utilities.hpp"
@@ -28,26 +28,26 @@
 namespace Kratos
 {
 
-InterfaceCoulombWithTensionCutOff::InterfaceCoulombWithTensionCutOff(std::unique_ptr<ConstitutiveLawDimension> pConstitutiveDimension)
+InterfaceCoulombLaw::InterfaceCoulombLaw(std::unique_ptr<ConstitutiveLawDimension> pConstitutiveDimension)
     : mpConstitutiveDimension(std::move(pConstitutiveDimension)),
       mTractionVector(ZeroVector(mpConstitutiveDimension->GetStrainSize())),
       mTractionVectorFinalized(ZeroVector(mpConstitutiveDimension->GetStrainSize())),
       mRelativeDisplacementVectorFinalized(ZeroVector(mpConstitutiveDimension->GetStrainSize()))
 {
 }
 
-ConstitutiveLaw::Pointer InterfaceCoulombWithTensionCutOff::Clone() const
+ConstitutiveLaw::Pointer InterfaceCoulombLaw::Clone() const
 {
-    auto p_result = std::make_shared<InterfaceCoulombWithTensionCutOff>(mpConstitutiveDimension->Clone());
+    auto p_result = std::make_shared<InterfaceCoulombLaw>(mpConstitutiveDimension->Clone());
     p_result->mTractionVector                      = mTractionVector;
     p_result->mTractionVectorFinalized             = mTractionVectorFinalized;
     p_result->mRelativeDisplacementVectorFinalized = mRelativeDisplacementVectorFinalized;
-    p_result->mCoulombWithTensionCutOffImpl        = mCoulombWithTensionCutOffImpl;
+    p_result->mCoulombImpl                         = mCoulombImpl;
     p_result->mIsModelInitialized                  = mIsModelInitialized;
     return p_result;
 }
 
-Vector& InterfaceCoulombWithTensionCutOff::GetValue(const Variable<Vector>& rVariable, Vector& rValue)
+Vector& InterfaceCoulombLaw::GetValue(const Variable<Vector>& rVariable, Vector& rValue)
 {
     if (rVariable == GEO_EFFECTIVE_TRACTION_VECTOR) {
         rValue = mTractionVector;
@@ -57,17 +57,15 @@ Vector& InterfaceCoulombWithTensionCutOff::GetValue(const Variable<Vector>& rVar
     return rValue;
 }
 
-int& InterfaceCoulombWithTensionCutOff::GetValue(const Variable<int>& rVariable, int& rValue)
+int& InterfaceCoulombLaw::GetValue(const Variable<int>& rVariable, int& rValue)
 {
     if (rVariable == GEO_PLASTICITY_STATUS) {
-        rValue = static_cast<int>(mCoulombWithTensionCutOffImpl.GetPlasticityStatus());
+        rValue = static_cast<int>(mCoulombImpl.GetPlasticityStatus());
     }
     return rValue;
 }
 
-void InterfaceCoulombWithTensionCutOff::SetValue(const Variable<Vector>& rVariable,
-                                                 const Vector&           rValue,
-                                                 const ProcessInfo&      rCurrentProcessInfo)
+void InterfaceCoulombLaw::SetValue(const Variable<Vector>& rVariable, const Vector& rValue, const ProcessInfo& rCurrentProcessInfo)
 {
     if (rVariable == GEO_EFFECTIVE_TRACTION_VECTOR) {
         mTractionVector = rValue;
@@ -76,15 +74,15 @@ void InterfaceCoulombWithTensionCutOff::SetValue(const Variable<Vector>& rVariab
     }
 }
 
-SizeType InterfaceCoulombWithTensionCutOff::WorkingSpaceDimension()
+SizeType InterfaceCoulombLaw::WorkingSpaceDimension()
 {
     // Note that this implementation assumes line interface elements. It needs to be modified when planar interface elements become available.
     return N_DIM_2D;
 }
 
-int InterfaceCoulombWithTensionCutOff::Check(const Properties&   rMaterialProperties,
-                                             const GeometryType& rElementGeometry,
-                                             const ProcessInfo&  rCurrentProcessInfo) const
+int InterfaceCoulombLaw::Check(const Properties&   rMaterialProperties,
+                               const GeometryType& rElementGeometry,
+                               const ProcessInfo&  rCurrentProcessInfo) const
 {
     const auto result = ConstitutiveLaw::Check(rMaterialProperties, rElementGeometry, rCurrentProcessInfo);
 
@@ -93,48 +91,48 @@ int InterfaceCoulombWithTensionCutOff::Check(const Properties&   rMaterialProper
     constexpr auto max_value_angle = 90.0;
     check_properties.SingleUseBounds(CheckProperties::Bounds::AllExclusive).Check(GEO_FRICTION_ANGLE, 0.0, max_value_angle);
     check_properties.Check(GEO_DILATANCY_ANGLE, rMaterialProperties[GEO_FRICTION_ANGLE]);
-    check_properties.Check(
-        GEO_TENSILE_STRENGTH,
-        rMaterialProperties[GEO_COHESION] /
-            std::tan(MathUtils<>::DegreesToRadians(rMaterialProperties[GEO_FRICTION_ANGLE])));
+    if (ConstitutiveLawUtilities::WantTensionCutOff(rMaterialProperties)) {
+        check_properties.Check(
+            GEO_TENSILE_STRENGTH,
+            rMaterialProperties[GEO_COHESION] /
+                std::tan(MathUtils<>::DegreesToRadians(rMaterialProperties[GEO_FRICTION_ANGLE])));
+    }
     check_properties.Check(INTERFACE_NORMAL_STIFFNESS);
     check_properties.Check(INTERFACE_SHEAR_STIFFNESS);
     return result;
 }
 
-ConstitutiveLaw::StressMeasure InterfaceCoulombWithTensionCutOff::GetStressMeasure()
+ConstitutiveLaw::StressMeasure InterfaceCoulombLaw::GetStressMeasure()
 {
     return ConstitutiveLaw::StressMeasure_Cauchy;
 }
 
-SizeType InterfaceCoulombWithTensionCutOff::GetStrainSize() const
+SizeType InterfaceCoulombLaw::GetStrainSize() const
 {
     // Note that this implementation assumes line interface elements. It needs to be modified when planar interface elements become available.
     return VOIGT_SIZE_2D_INTERFACE;
 }
 
-ConstitutiveLaw::StrainMeasure InterfaceCoulombWithTensionCutOff::GetStrainMeasure()
+ConstitutiveLaw::StrainMeasure InterfaceCoulombLaw::GetStrainMeasure()
 {
     return ConstitutiveLaw::StrainMeasure_Infinitesimal;
 }
 
-bool InterfaceCoulombWithTensionCutOff::IsIncremental() { return true; }
+bool InterfaceCoulombLaw::IsIncremental() { return true; }
 
-bool InterfaceCoulombWithTensionCutOff::RequiresInitializeMaterialResponse() { return true; }
+bool InterfaceCoulombLaw::RequiresInitializeMaterialResponse() { return true; }
 
-void InterfaceCoulombWithTensionCutOff::InitializeMaterial(const Properties& rMaterialProperties,
-                                                           const Geometry<Node>&,
-                                                           const Vector&)
+void InterfaceCoulombLaw::InitializeMaterial(const Properties& rMaterialProperties, const Geometry<Node>&, const Vector&)
 {
-    mCoulombWithTensionCutOffImpl = CoulombWithTensionCutOffImpl{rMaterialProperties};
+    mCoulombImpl = CoulombImpl{rMaterialProperties};
 
     mRelativeDisplacementVectorFinalized =
         HasInitialState() ? GetInitialState().GetInitialStrainVector() : ZeroVector{GetStrainSize()};
     mTractionVectorFinalized =
         HasInitialState() ? GetInitialState().GetInitialStressVector() : ZeroVector{GetStrainSize()};
 }
 
-void InterfaceCoulombWithTensionCutOff::InitializeMaterialResponseCauchy(Parameters& rConstitutiveLawParameters)
+void InterfaceCoulombLaw::InitializeMaterialResponseCauchy(Parameters& rConstitutiveLawParameters)
 {
     if (!mIsModelInitialized) {
         mTractionVectorFinalized             = rConstitutiveLawParameters.GetStressVector();
@@ -143,7 +141,7 @@ void InterfaceCoulombWithTensionCutOff::InitializeMaterialResponseCauchy(Paramet
     }
 }
 
-void InterfaceCoulombWithTensionCutOff::CalculateMaterialResponseCauchy(Parameters& rConstitutiveLawParameters)
+void InterfaceCoulombLaw::CalculateMaterialResponseCauchy(Parameters& rConstitutiveLawParameters)
 {
     if (!rConstitutiveLawParameters.GetOptions().Is(ConstitutiveLaw::COMPUTE_STRESS)) {
         return;
@@ -159,8 +157,8 @@ void InterfaceCoulombWithTensionCutOff::CalculateMaterialResponseCauchy(Paramete
     const auto negative   = std::signbit(trial_sigma_tau.Tau());
     trial_sigma_tau.Tau() = std::abs(trial_sigma_tau.Tau());
 
-    if (!mCoulombWithTensionCutOffImpl.IsAdmissibleStressState(trial_sigma_tau)) {
-        mapped_sigma_tau = mCoulombWithTensionCutOffImpl.DoReturnMapping(
+    if (!mCoulombImpl.IsAdmissibleStressState(trial_sigma_tau)) {
+        mapped_sigma_tau = mCoulombImpl.DoReturnMapping(
             trial_sigma_tau, mpConstitutiveDimension->CalculateElasticConstitutiveTensor(r_properties),
             Geo::PrincipalStresses::AveragingType::NO_AVERAGING);
         if (negative) mapped_sigma_tau.Tau() *= -1.0;
@@ -170,9 +168,9 @@ void InterfaceCoulombWithTensionCutOff::CalculateMaterialResponseCauchy(Paramete
     rConstitutiveLawParameters.GetStressVector() = mTractionVector;
 }
 
-Geo::SigmaTau InterfaceCoulombWithTensionCutOff::CalculateTrialTractionVector(const Vector& rRelativeDisplacementVector,
-                                                                              double NormalStiffness,
-                                                                              double ShearStiffness) const
+Geo::SigmaTau InterfaceCoulombLaw::CalculateTrialTractionVector(const Vector& rRelativeDisplacementVector,
+                                                                double NormalStiffness,
+                                                                double ShearStiffness) const
 {
     constexpr auto number_of_normal_components = std::size_t{1};
     return Geo::SigmaTau{mTractionVectorFinalized +
@@ -181,15 +179,15 @@ Geo::SigmaTau InterfaceCoulombWithTensionCutOff::CalculateTrialTractionVector(co
                               rRelativeDisplacementVector - mRelativeDisplacementVectorFinalized)};
 }
 
-void InterfaceCoulombWithTensionCutOff::FinalizeMaterialResponseCauchy(Parameters& rConstitutiveLawParameters)
+void InterfaceCoulombLaw::FinalizeMaterialResponseCauchy(Parameters& rConstitutiveLawParameters)
 {
     mRelativeDisplacementVectorFinalized = rConstitutiveLawParameters.GetStrainVector();
     mTractionVectorFinalized             = mTractionVector;
 }
 
-Matrix& InterfaceCoulombWithTensionCutOff::CalculateValue(Parameters& rConstitutiveLawParameters,
-                                                          const Variable<Matrix>& rVariable,
-                                                          Matrix&                 rValue)
+Matrix& InterfaceCoulombLaw::CalculateValue(Parameters&             rConstitutiveLawParameters,
+                                            const Variable<Matrix>& rVariable,
+                                            Matrix&                 rValue)
 {
     if (rVariable == CONSTITUTIVE_MATRIX) {
         const auto&    r_properties = rConstitutiveLawParameters.GetMaterialProperties();
@@ -204,25 +202,25 @@ Matrix& InterfaceCoulombWithTensionCutOff::CalculateValue(Parameters& rConstitut
     return rValue;
 }
 
-void InterfaceCoulombWithTensionCutOff::save(Serializer& rSerializer) const
+void InterfaceCoulombLaw::save(Serializer& rSerializer) const
 {
     KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, ConstitutiveLaw)
     rSerializer.save("ConstitutiveDimension", mpConstitutiveDimension);
     rSerializer.save("TractionVector", mTractionVector);
     rSerializer.save("TractionVectorFinalized", mTractionVectorFinalized);
     rSerializer.save("RelativeDisplacementVectorFinalized", mRelativeDisplacementVectorFinalized);
-    rSerializer.save("CoulombWithTensionCutOffImpl", mCoulombWithTensionCutOffImpl);
+    rSerializer.save("CoulombImpl", mCoulombImpl);
     rSerializer.save("IsModelInitialized", mIsModelInitialized);
 }
 
-void InterfaceCoulombWithTensionCutOff::load(Serializer& rSerializer)
+void InterfaceCoulombLaw::load(Serializer& rSerializer)
 {
     KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, ConstitutiveLaw)
     rSerializer.load("ConstitutiveDimension", mpConstitutiveDimension);
     rSerializer.load("TractionVector", mTractionVector);
     rSerializer.load("TractionVectorFinalized", mTractionVectorFinalized);
     rSerializer.load("RelativeDisplacementVectorFinalized", mRelativeDisplacementVectorFinalized);
-    rSerializer.load("CoulombWithTensionCutOffImpl", mCoulombWithTensionCutOffImpl);
+    rSerializer.load("CoulombImpl", mCoulombImpl);
     rSerializer.load("IsModelInitialized", mIsModelInitialized);
 }
 

applications/GeoMechanicsApplication/custom_constitutive/interface_coulomb_with_tension_cut_off.h → applications/GeoMechanicsApplication/custom_constitutive/interface_coulomb_law.h

@@ -14,29 +14,29 @@
 
 #pragma once
 
-#include "custom_constitutive/coulomb_with_tension_cut_off_impl.h"
+#include "custom_constitutive/coulomb_impl.h"
 #include "includes/constitutive_law.h"
 
 namespace Kratos
 {
 
 class ConstitutiveLawDimension;
 
-class KRATOS_API(GEO_MECHANICS_APPLICATION) InterfaceCoulombWithTensionCutOff : public ConstitutiveLaw
+class KRATOS_API(GEO_MECHANICS_APPLICATION) InterfaceCoulombLaw : public ConstitutiveLaw
 {
 public:
-    KRATOS_CLASS_POINTER_DEFINITION(InterfaceCoulombWithTensionCutOff);
+    KRATOS_CLASS_POINTER_DEFINITION(InterfaceCoulombLaw);
 
-    InterfaceCoulombWithTensionCutOff() = default;
-    explicit InterfaceCoulombWithTensionCutOff(std::unique_ptr<ConstitutiveLawDimension> pConstitutiveDimension);
+    InterfaceCoulombLaw() = default;
+    explicit InterfaceCoulombLaw(std::unique_ptr<ConstitutiveLawDimension> pConstitutiveDimension);
 
     // Copying is not allowed. Use member `Clone` instead.
-    InterfaceCoulombWithTensionCutOff(const InterfaceCoulombWithTensionCutOff&)            = delete;
-    InterfaceCoulombWithTensionCutOff& operator=(const InterfaceCoulombWithTensionCutOff&) = delete;
+    InterfaceCoulombLaw(const InterfaceCoulombLaw&)            = delete;
+    InterfaceCoulombLaw& operator=(const InterfaceCoulombLaw&) = delete;
 
     // Moving is supported
-    InterfaceCoulombWithTensionCutOff(InterfaceCoulombWithTensionCutOff&&) noexcept = default;
-    InterfaceCoulombWithTensionCutOff& operator=(InterfaceCoulombWithTensionCutOff&&) noexcept = default;
+    InterfaceCoulombLaw(InterfaceCoulombLaw&&) noexcept            = default;
+    InterfaceCoulombLaw& operator=(InterfaceCoulombLaw&&) noexcept = default;
 
     [[nodiscard]] ConstitutiveLaw::Pointer Clone() const override;
     SizeType                               WorkingSpaceDimension() override;
@@ -69,7 +69,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) InterfaceCoulombWithTensionCutOff :
     Vector                                    mTractionVector;
     Vector                                    mTractionVectorFinalized;
     Vector                                    mRelativeDisplacementVectorFinalized;
-    CoulombWithTensionCutOffImpl              mCoulombWithTensionCutOffImpl;
+    CoulombImpl                               mCoulombImpl;
     bool                                      mIsModelInitialized = false;
 
     [[nodiscard]] Geo::SigmaTau CalculateTrialTractionVector(const Vector& rRelativeDisplacementVector,
@@ -79,6 +79,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) InterfaceCoulombWithTensionCutOff :
     friend class Serializer;
     void save(Serializer& rSerializer) const override;
     void load(Serializer& rSerializer) override;
-}; // Class InterfaceMohrCoulombWithTensionCutOff
+}; // Class InterfaceMohrCoulomb
 
 } // namespace Kratos