[GeoMechanicsApplicaiton] Make tension cut off truly optional for Mohr Coulomb model by rfaasse · Pull Request #14467 · KratosMultiphysics/Kratos

rfaasse · 2026-06-03T07:18:11Z

📝 Description

This PR enables the user to specify a Mohr Coulomb model with or without tension cut-off, using the new GEO_ENABLE_TENSION_CUT_OFF parameter in the material parameters. When the tension cut-off is turned off, the Mohr Coulomb model collapses to a simplified version with just the regular failure zone, elastic zone and apex return.

Key changes

The GEO_ENABLE_TENSION_CUT_OFF is added as a material parameter for the Mohr Coulomb model. If it's not specified, but a tensile strength is specified, it is assumed the tension cut-off is enabled (for backward compatibility).
Unit tests and integration tests are adjusted (GEO_ENABLE_TENSION_CUT_OFF with the correct value is added). For some tests, the tensile strength was specified as cohesion/tan(phi), which is equivalent to disabling the tension cut-off. For the tests where this was the case, the GEO_ENABLE_TENSION_CUT_OFF is set to false
Added equivalent registered names for "GeoMohrCoulombWithTensionCutOff2D", "GeoMohrCoulombWithTensionCutOff3D", "GeoInterfaceCoulombWithTensionCutOff" -> "GeoMohrCoulomb2D", "GeoMohrCoulomb3D", "GeoInterfaceCoulomb". The old names are kept for now, but deprecated (to be removed after some time).

Note: A PR for the user documentation will be created separately

…on cutoff object based on that

…ion-cut-off-truly-optional

…-off

… values) in all crow materials

…ility)

rfaasse · 2026-06-04T13:25:27Z

+    // Keeping the old names for backward compatibility
+    KRATOS_REGISTER_CONSTITUTIVE_LAW("GeoMohrCoulombWithTensionCutOff2D", mMohrCoulomb2D)
+    KRATOS_REGISTER_CONSTITUTIVE_LAW("GeoMohrCoulombWithTensionCutOff3D", mMohrCoulomb3D)
+    KRATOS_REGISTER_CONSTITUTIVE_LAW("GeoInterfaceCoulombWithTensionCutOff", mInterfaceCoulomb)


I was doubting whether we should keep the old ones. Curious what others think

…enabled

WPK4FEM

HI Richard,
This is what I managed to read and comment on today, I will have another look.
Wijtze Pieter

WPK4FEM · 2026-06-04T15:07:31Z

    const auto coulomb_yield_function_value = mCoulombYieldSurface.YieldFunctionValue(rTrialStressState);
-    const auto tension_yield_function_value = mTensionCutOff.YieldFunctionValue(rTrialStressState);
-    constexpr auto tolerance                = 1.0e-10;
+    constexpr auto tolerance         = 1.0e-10;


Please relate there to the available relative or absolute numerical tolerances.

These are in test_utilities.h (so in the test folder) and therefore cannot be reached here. Also, the defaults are 1e-6 for relative and 1e-12 for absolute. What do you think we should use here? I wasn't that involved in the development of Coulomb, so I'm not sure why this exact tolerance was chosen at the time

WPK4FEM · 2026-06-04T15:11:12Z

+    auto tension_admissable = true;
+    if (mTensionCutOff) {
+        const auto tension_yield_function_value = mTensionCutOff->YieldFunctionValue(rTrialStressState);
+        const auto tension_tolerance = tolerance * (1.0 + std::abs(tension_yield_function_value));
+        tension_admissable           = tension_yield_function_value < tension_tolerance;
+    }
+
+    const auto admissible_state = coulomb_yield_function_value < coulomb_tolerance && tension_admissable;


The MC part is always there so lets start from that

Suggested change

auto tension_admissable = true;

if (mTensionCutOff) {

const auto tension_yield_function_value = mTensionCutOff->YieldFunctionValue(rTrialStressState);

const auto tension_tolerance = tolerance * (1.0 + std::abs(tension_yield_function_value));

tension_admissable = tension_yield_function_value < tension_tolerance;

}

const auto admissible_state = coulomb_yield_function_value < coulomb_tolerance && tension_admissable;

auto admissible_state = coulomb_yield_function_value < coulomb_tolerance;

if (admissible_state && mTensionCutOff) {

const auto tension_yield_function_value = mTensionCutOff->YieldFunctionValue(rTrialStressState);

const auto tension_tolerance = tolerance * (1.0 + std::abs(tension_yield_function_value));

admissible_state = tension_yield_function_value < tension_tolerance;

}

Nice suggestion, done!

WPK4FEM · 2026-06-04T15:15:44Z

+        if (mTensionCutOff && IsStressAtTensionApexReturnZone(trial_traction)) {
            mPlasticityStatus = PlasticityStatus::TENSION_APEX;
            return ReturnStressAtTensionApexReturnZone(rTrialStressState);
        }

-        if (IsStressAtTensionCutoffReturnZone(trial_traction)) {
+        if (mTensionCutOff && IsStressAtTensionCutoffReturnZone(trial_traction)) {
            mPlasticityStatus = PlasticityStatus::TENSION_CUT_OFF;
            return ReturnStressAtTensionCutoffReturnZone(rTrialStressState,
                                                         rElasticConstitutiveTensor, AveragingType);
        }


Both previous conditions need the TensionCutOff to be active so lets check that in an encompassing condition.

Suggested change

if (mTensionCutOff) {

if (IsStressAtTensionApexReturnZone(trial_traction)) {

mPlasticityStatus = PlasticityStatus::TENSION_APEX;

return ReturnStressAtTensionApexReturnZone(rTrialStressState);

}

if (IsStressAtTensionCutoffReturnZone(trial_traction)) {

mPlasticityStatus = PlasticityStatus::TENSION_CUT_OFF;

return ReturnStressAtTensionCutoffReturnZone(rTrialStressState,

rElasticConstitutiveTensor, AveragingType);

}

}

I agree with this suggestion. 👍

avdg81

Hi Richard,
Thanks a lot for making the tension cut-off of the Mohr-Coulomb model truly optional. With your changes, there's no doubt about when a tension cut-off was intended or not. It also helps to structure the code even more.
I have several questions and suggestions. Hope all of them will help to make your valuable changes stand out even more. 👍

avdg81 · 2026-06-08T09:09:27Z

+    if (rMaterialProperties.Has(GEO_ENABLE_TENSION_CUT_OFF) && rMaterialProperties[GEO_ENABLE_TENSION_CUT_OFF]) {
+        return std::make_optional<TensionCutoff>(rMaterialProperties[GEO_TENSILE_STRENGTH]);
+    }
+
+    // The following statement is to support backward compatibility (i.e. GEO_ENABLE_TENSION_CUT_OFF is not specified,
+    // but the GEO_TENSILE_STRENGTH is provided), which results in an enabled tension cutoff.
+    if (rMaterialProperties.Has(GEO_TENSILE_STRENGTH) && !rMaterialProperties.Has(GEO_ENABLE_TENSION_CUT_OFF)) {
+        return std::make_optional<TensionCutoff>(rMaterialProperties[GEO_TENSILE_STRENGTH]);
+    }


Since the return statements are identical, I would suggest to join the condition checks and put them in a local helper function named WantTensionCutOff (or any other name if you can think of a better one):

Suggested change

if (rMaterialProperties.Has(GEO_ENABLE_TENSION_CUT_OFF) && rMaterialProperties[GEO_ENABLE_TENSION_CUT_OFF]) {

return std::make_optional<TensionCutoff>(rMaterialProperties[GEO_TENSILE_STRENGTH]);

}

// The following statement is to support backward compatibility (i.e. GEO_ENABLE_TENSION_CUT_OFF is not specified,

// but the GEO_TENSILE_STRENGTH is provided), which results in an enabled tension cutoff.

if (rMaterialProperties.Has(GEO_TENSILE_STRENGTH) && !rMaterialProperties.Has(GEO_ENABLE_TENSION_CUT_OFF)) {

return std::make_optional<TensionCutoff>(rMaterialProperties[GEO_TENSILE_STRENGTH]);

}

if (WantTensionCutOff(rMaterialProperties)) {

return std::make_optional<TensionCutoff>(rMaterialProperties[GEO_TENSILE_STRENGTH]);

}

where

bool WantTensionCutOff(const Properties& rMaterialProperties) { if (rMaterialProperties.Has(GEO_ENABLE_TENSION_CUT_OFF) && rMaterialProperties[GEO_ENABLE_TENSION_CUT_OFF]) return true; // The following statement is to support backward compatibility (i.e. GEO_ENABLE_TENSION_CUT_OFF is not specified, // but the GEO_TENSILE_STRENGTH is provided), which results in an enabled tension cutoff. return !rMaterialProperties.Has(GEO_ENABLE_TENSION_CUT_OFF) && rMaterialProperties.Has(GEO_TENSILE_STRENGTH); }

As you can see, I have also swapped the two conditions in the backward compatibility check. In that way, the check whether the flag GEO_ENABLE_TENSION_CUT_OFF exists always comes first, which is a little bit clearer for me.

avdg81 · 2026-06-08T12:15:43Z


 template <typename StressStateType>
-bool CoulombWithTensionCutOffImpl::IsAdmissibleStressState(const StressStateType& rTrialStressState)
+bool CoulombImpl::IsAdmissibleStressState(const StressStateType& rTrialStressState)


Perhaps we can extract the checking whether the given trial stress state is admissible for a given yield surface to a new helper function, and then use that to implement this member function.

In the anonymous we can add the following function template:

template <typename YieldSurfaceType, typename StressStateType> bool IsAdmissibleStressState(const YieldSurfaceType& rYieldSurface, const StressStateType& rTrialStressState) { const auto yield_function_value = rYieldSurface.YieldFunctionValue(rTrialStressState); const auto tolerance = 1.0e-10 * (1.0 + std::abs(yield_function_value)); return yield_function_value < tolerance; }

And then the member function template can be simplified as follows:

template <typename StressStateType> bool CoulombImpl::IsAdmissibleStressState(const StressStateType& rTrialStressState) { if (!::IsAdmissibleStressState(mCoulombYieldSurface, rTrialStressState)) return false; if (mTensionCutOff && !::IsAdmissibleStressState(*mTensionCutOff, rTrialStressState)) return false; mPlasticityStatus = PlasticityStatus::ELASTIC; return true; }

Note that with this new structure in-place, adding the optional cap yield surface becomes straightforward.

avdg81 · 2026-06-08T12:32:55Z

+        if (mTensionCutOff && IsStressAtTensionApexReturnZone(trial_traction)) {
            mPlasticityStatus = PlasticityStatus::TENSION_APEX;
            return ReturnStressAtTensionApexReturnZone(rTrialStressState);
        }

-        if (IsStressAtTensionCutoffReturnZone(trial_traction)) {
+        if (mTensionCutOff && IsStressAtTensionCutoffReturnZone(trial_traction)) {
            mPlasticityStatus = PlasticityStatus::TENSION_CUT_OFF;
            return ReturnStressAtTensionCutoffReturnZone(rTrialStressState,
                                                         rElasticConstitutiveTensor, AveragingType);
        }


I agree with this suggestion. 👍

avdg81 · 2026-06-08T13:34:44Z

+Geo::SigmaTau CoulombImpl::CalculateCornerPoint() const
 {
-    const auto tensile_strength = mTensionCutOff.GetTensileStrength();
+    if (const auto apex = mCoulombYieldSurface.CalculateApex(); !mTensionCutOff) return apex;


Since apex is not being used by the check itself, I feel we'd better postpone its calculation to the latest possible moment:

Suggested change

if (const auto apex = mCoulombYieldSurface.CalculateApex(); !mTensionCutOff) return apex;

if (!mTensionCutOff) return mCoulombYieldSurface.CalculateApex();

(And let's add an additional blank line to clearly show that we've handled the first case.)

avdg81 · 2026-06-08T15:06:36Z

-    std::size_t         mMaxNumberOfPlasticIterations{100};
-    PlasticityStatus    mPlasticityStatus{PlasticityStatus::ELASTIC};
+    CoulombYieldSurface          mCoulombYieldSurface;
+    std::optional<TensionCutoff> mTensionCutOff;


Would a name like mOptionalTensionCutOff better reflect that the tension cut-off is optional? In other words, that it may or may not be there.

avdg81 · 2026-06-08T15:40:30Z

+                    "GEO_ENABLE_TENSION_CUT_OFF": true,
                    "GEO_TENSILE_STRENGTH": 0.0,


Suggested change

"GEO_ENABLE_TENSION_CUT_OFF": true,

"GEO_TENSILE_STRENGTH": 0.0,

"GEO_ENABLE_TENSION_CUT_OFF": false,

avdg81 · 2026-06-08T19:22:54Z

+    KRATOS_REGISTER_CONSTITUTIVE_LAW("GeoMohrCoulomb2D", mMohrCoulomb2D)
+    KRATOS_REGISTER_CONSTITUTIVE_LAW("GeoMohrCoulomb3D", mMohrCoulomb3D)
+    KRATOS_REGISTER_CONSTITUTIVE_LAW("GeoInterfaceCoulomb", mInterfaceCoulomb)


These new names are slightly different than what was proposed in the issue:

"GeoMohrCoulombWithTensionCutOff2D" => "GeoMohrCoulombLawPlaneStrain"
"GeoMohrCoulombWithTensionCutOff3D" => "GeoMohrCoulombLaw3D"
"GeoInterfaceCoulombWithTensionCutOff" => "GeoInterfaceCoulombLawPlaneStrain"

I'm wondering why the names are somewhat different?

avdg81 · 2026-06-08T19:24:40Z


-    const MohrCoulombWithTensionCutOff mMohrCoulombWithTensionCutOff2D{std::make_unique<PlaneStrain>()};
-    const MohrCoulombWithTensionCutOff mMohrCoulombWithTensionCutOff3D{std::make_unique<ThreeDimensional>()};
+    const MohrCoulomb mMohrCoulomb2D{std::make_unique<PlaneStrain>()};


I feel mMohrCoulombPlaneStrain is a better name (or perhaps even better: mMohrCoulombLawPlaneStrain). Perhaps @WPK4FEM wants to weigh in here, too?

avdg81 · 2026-06-08T19:32:41Z

-    const auto tensile_strength = cohesion / std::tan(MathUtils<>::DegreesToRadians(phi_in_degrees));
-    properties.SetValue(GEO_TENSILE_STRENGTH, tensile_strength);


Perhaps it's clearer to be explicit about not wanting a tension cut-off rather than relying on a default value:

properties.SetValue(GEO_ENABLE_TENSION_CUT_OFF, false);

avdg81 · 2026-06-08T19:34:07Z

    // Act
-    auto cauchy_stress_vector = UblasUtilities::CreateVector(
-        {tensile_strength + 20.0, tensile_strength + 10.0, tensile_strength, 0.0});
+    const auto sigma_corner = cohesion / std::tan(MathUtils<>::DegreesToRadians(phi_in_degrees));


Perhaps apex is a better name?

Suggested change

const auto sigma_corner = cohesion / std::tan(MathUtils<>::DegreesToRadians(phi_in_degrees));

const auto apex = cohesion / std::tan(MathUtils<>::DegreesToRadians(phi_in_degrees));

rfaasse added 6 commits June 1, 2026 11:23

Changed type to std::optional, but still assume it's always filled

2944d3f

Added the GEO_ENABLE_TENSION_CUT_OFF variable and construct the tensi…

b803901

…on cutoff object based on that

Fixed typo in README of Mohr Coulomb

2371c82

Made the usage of the tension cut-off optional

e959925

Removed redundant parentheses

60b720e

Clarified meaning of variable (hopefully correct)

78b6f0d

rfaasse self-assigned this Jun 3, 2026

rfaasse linked an issue Jun 3, 2026 that may be closed by this pull request

[GeoMechanicsApplication] Make tension cut-off truly optional #14463

Open

rfaasse added 11 commits June 3, 2026 09:34

Merge remote-tracking branch 'origin/master' into geo/14463-make-tens…

1a90355

…ion-cut-off-truly-optional

Attempt at fixing serialization

20a319f

Converted an integration test to use mohr coulomb without tension cut…

f6e480c

…-off

Renamed laws and impl

e16733c

Moved tensile bool above the value

3166bd6

fixed serialization tests and backward compatibility

b82c4cd

Removed tension cut-off explicitly (instead of implicitly through the…

73ad2ac

… values) in all crow materials

Renamed variable

4f608a5

Formatting

e60a27f

Added simplified registered names (but kept the old ones for compatib…

bdfdfbd

…ility)

Formatting

eb93996

rfaasse commented Jun 4, 2026

View reviewed changes

rfaasse added 3 commits June 4, 2026 16:01

Adjusted the readme after making tension cutoff optional

8e982fd

Added warning if tensile strength is set, but tension cut off is not …

d651da7

…enabled

Clarified comment about backward compatibility

587162c

rfaasse added the GeoMechanics Issues related to the GeoMechanicsApplication label Jun 4, 2026

rfaasse added 2 commits June 4, 2026 16:23

Clarified comment about deprecated names

d1d4ef0

Made function private and static

719b214

rfaasse marked this pull request as ready for review June 4, 2026 14:41

rfaasse requested a review from a team as a code owner June 4, 2026 14:41

WPK4FEM reviewed Jun 4, 2026

View reviewed changes

avdg81 reviewed Jun 8, 2026

View reviewed changes

	if (const auto apex = mCoulombYieldSurface.CalculateApex(); !mTensionCutOff) return apex;
	if (!mTensionCutOff) return mCoulombYieldSurface.CalculateApex();

		"GEO_ENABLE_TENSION_CUT_OFF": true,
		"GEO_TENSILE_STRENGTH": 0.0,

	"GEO_ENABLE_TENSION_CUT_OFF": true,
	"GEO_TENSILE_STRENGTH": 0.0,
	"GEO_ENABLE_TENSION_CUT_OFF": false,

		const auto tensile_strength = cohesion / std::tan(MathUtils<>::DegreesToRadians(phi_in_degrees));
		properties.SetValue(GEO_TENSILE_STRENGTH, tensile_strength);

	const auto sigma_corner = cohesion / std::tan(MathUtils<>::DegreesToRadians(phi_in_degrees));
	const auto apex = cohesion / std::tan(MathUtils<>::DegreesToRadians(phi_in_degrees));

Conversation

rfaasse commented Jun 3, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

📝 Description

Key changes

Uh oh!

Choose a reason for hiding this comment

Uh oh!

WPK4FEM left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

avdg81 left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

rfaasse commented Jun 3, 2026 •

edited

Loading