[GeoMechanicsApplication] Improve the Coulomb model framework by avdg81 · Pull Request #14187 · KratosMultiphysics/Kratos

avdg81 · 2026-02-06T16:08:45Z

📝 Description
The Coulomb model now accepts principal stresses as well as sigma and tau. This also applies to the tension cut-off.

When these changes are in place, we can start to fix the incorrect return mapping of the tension cut-off as well as include the effect of the elastic constitutive tensor.

The new name makes clearer what it returns.

This resolves an issue found by SonarQube.

…RotationMatrix`

Also renamed some variables and made a few ones `const`.

…unction`

…he-coulomb-model-framework

WPK4FEM

Hi Anne,
In retrospect, my review comments are more about the future use than about the current state. I think the PR is a step in a good direction.
Regards, Wijtze Pieter

WPK4FEM · 2026-02-09T09:41:49Z

applications/GeoMechanicsApplication/custom_constitutive/coulomb_with_tension_cut_off_impl.cpp

+#include "custom_constitutive/geo_sigma_tau.hpp"
+#include "custom_constitutive/principal_stresses.hpp"


At first glance the naming is not consistent. Why geo_ preceding sigma_tau and omitting geo_ before principal_stresses in the names of the hpp files.

Well, both classes were prepared by two different developers who adopted different naming schemes. After a brief discussion with you, @rfaasse and myself, we decided to remove the geo_ prefix from geo_sigma_tau.{cpp,hpp}.

WPK4FEM · 2026-02-09T10:03:32Z

applications/GeoMechanicsApplication/custom_constitutive/coulomb_with_tension_cut_off_impl.cpp

+Geo::SigmaTau CoulombWithTensionCutOffImpl::ReturnStressAtTensionCutoffReturnZone(const Geo::SigmaTau& rSigmaTau) const
 {
    const auto derivative_of_flow_function = mTensionCutOff.DerivativeOfFlowFunction(rSigmaTau);
-    const auto lambda_tc = (mTensionCutOff.GetTensileStrength() - rSigmaTau[0] - rSigmaTau[1]) /
+    const auto lambda_tc = (mTensionCutOff.GetTensileStrength() - rSigmaTau.Sigma() - rSigmaTau.Tau()) /
                           (derivative_of_flow_function[0] + derivative_of_flow_function[1]);
-    return rSigmaTau + lambda_tc * derivative_of_flow_function;
+    return Geo::SigmaTau{rSigmaTau.Values() + lambda_tc * derivative_of_flow_function};
 }


now the lambda computation only depends on the content of the tension cutoff yield surface, but we will need a correction that includes the elastic constitutive tensor. After that correction should this lambda computation be here or a part of the tension cut off yield surface. The latter would make it necessary to know about the elastic constitutive tensor in the yield surfaces.

I would say that the answer depends on what is conceptually most natural. When the elastic constitutive tensor is naturally associated with the yield surface, the tensor can be part of it. If it's more of a case where the yield surface just uses the elastic constitutive tensor, we may simply pass it through the argument list. What do you think is most natural?

Pass it through the argument list, otherwise the yield surface may need to know about non-linear elasticity.

WPK4FEM · 2026-02-09T10:05:15Z

applications/GeoMechanicsApplication/custom_constitutive/coulomb_with_tension_cut_off_impl.cpp

+
+Geo::SigmaTau CoulombWithTensionCutOffImpl::ReturnStressAtRegularFailureZone(
+    const Geo::SigmaTau& rSigmaTau, CoulombYieldSurface::CoulombAveragingType AveragingType) const
 {


Also this lambda computation will need to know about the elastic constitutive tensor.

WPK4FEM · 2026-02-09T10:11:47Z

...cations/GeoMechanicsApplication/custom_constitutive/interface_coulomb_with_tension_cut_off.h

    Vector                       mTractionVector;
    Vector                       mTractionVectorFinalized;


Should we change these to Geo::SigmaTau?

Over time, yes. But for now, no. The reason I'm against doing it now, is simply scope creep. When we change these data members, we are also forced to make class Geo::SigmaTau serializable (i.e. we need to implement and test load and save functions), or else the serialization of InterfaceCoulombWithTensionCutOff will be broken. All of this is doable, but not with just a few extra lines. Is it okay for you if I add a backlog issue for this suggestion, just so we can pick it up later?

rfaasse

Thanks for implementing this improvement to the framework! I have some conceptual questions and some minor suggestions for improvements, but I don't expect anything major will come out of it!

rfaasse · 2026-02-09T13:35:49Z

applications/GeoMechanicsApplication/custom_constitutive/coulomb_with_tension_cut_off_impl.cpp

+                                                            CoulombYieldSurface::CoulombAveragingType AveragingType)
+{
+    auto sigma_tau_to_sigma_tau = [](const Geo::SigmaTau& rSigmaTau) { return rSigmaTau; };
+    return DoReturnMapping<>(rTrialSigmaTau, sigma_tau_to_sigma_tau, AveragingType);


When seeing it like this, it feels a bit weird to pass this trivial conversion function. Would it be possible to have the conversion function as an optional? Then we just skip the conversion if it's not passed

Hmm, I'm not sure. We cannot really skip the conversion, unless we're willing to assume that when no conversion function is supplied, the given stress state is to be taken as a $\sigma$ and $\tau$ stress state. With the current setup we make things very explicit, which I like more than an implicit assumption. Perhaps @WPK4FEM would like to share his point of view?

When first reading, my reaction was the same as Richards. Still I would keep the explicit passing of a trivial conversion. We probably will encounter other sets of stress invariants, the name of the convertor is then helping to make explicit from which to which perception of stress we are moving.

rfaasse · 2026-02-09T13:48:05Z

applications/GeoMechanicsApplication/custom_constitutive/coulomb_with_tension_cut_off_impl.h

-    [[nodiscard]] bool   IsAdmissibleSigmaTau(const Vector& rTrialSigmaTau) const;
-    [[nodiscard]] Vector DoReturnMapping(const Vector&                             rTrialSigmaTau,
-                                         CoulombYieldSurface::CoulombAveragingType AveragingType);
+    [[nodiscard]] bool IsAdmissibleStressState(const Geo::SigmaTau& rTrialSigmaTau) const;


You could argue that the name of this argument doesn't need to repeat 'sigma tau' in the name, because it's already it's type (which would also result in code that has a bit less repetition, e.g. rTrialStress.Sigma() instead of rTrialSigmaTau.Sigma()). I'm find with both ways, but am slightly in favour of not repeating the type in variable/argument names.

In a separate conversation Anne asked me about this. The use of sigma and tau here refers to the tractions on a chosen plane ( which also is the description we choose in an interface element ). Therefore I would find the following natural:
Geo::SigmaTau rTrialTraction
and for selecting the components of this variable you then get rTrialTraction.sigma and rTrialTracion.tau

Then I will rename the objects of type Geo::SigmaTau to something with "traction" in their names. 👍

rfaasse · 2026-02-09T13:49:37Z

applications/GeoMechanicsApplication/custom_constitutive/coulomb_with_tension_cut_off_impl.cpp

+#include "custom_constitutive/sigma_tau.hpp"
 #include "custom_utilities/constitutive_law_utilities.h"
+#include "custom_utilities/stress_strain_utilities.h"
 #include "custom_utilities/ublas_utilities.h"


Due to the changes ublas_utilities.h became redundant as an include here.

Good catch! Removed the now redundant #include.

rfaasse · 2026-02-09T13:55:28Z

applications/GeoMechanicsApplication/custom_constitutive/coulomb_yield_surface.cpp

+    return YieldFunctionValue(StressStrainUtilities::TransformPrincipalStressesToSigmaTau(rPrincipalStresses));
+}
+
 Vector CoulombYieldSurface::DerivativeOfFlowFunction(const Vector& rSigmaTau) const


Probably the same holds here (we only keep it for now, but remove it after the PQ stress state is added?)

No, I think we will continue to use this specific one. It is being called from CoulombWithTensionCutOffImpl::IsAdmissibleStressState, which accepts either a principal stress state or a $\sigma$ and $\tau$. Member function MohrCoulombWithTensionCutOff::CalculateMaterialResponseCauchy which calls this check is most naturally expressed in terms of principal stresses, so I wouldn't want to move the conversion from principal stresses to $\sigma$ and $\tau$ to there. Let me know what you think.

rfaasse · 2026-02-09T13:57:48Z

applications/GeoMechanicsApplication/custom_constitutive/coulomb_yield_surface.cpp


 void CoulombYieldSurface::SetKappa(double kappa) { mKappa = kappa; }

+// At some point in time we would like to get rid of this API. For now, just forward the request.


Just to validate my understanding: is this only kept for now because we don't have the PQ stress state yet (and so have to implement the Vector version to keep supporting building the CompressionCap)?

It's only there because we have one more yield surface that doesn't use specific stress states yet (i.e. the compression cap). Since all yield surfaces still derive from YieldSurface, we must implement this interface (since they are pure virtual functions). However, the Coulomb yield surface and the tension cut-off no longer use this API (as you can see from the code coverage analysis). As soon as we have expressed the compression cap in terms of $p$ and $q$, we no longer need the YieldSurface interface, and these overrides can be removed. That was what I was trying to say with the comment "At some point in time we would like to get rid of this API. For now, just forward the request.".

rfaasse · 2026-02-09T14:03:33Z

applications/GeoMechanicsApplication/custom_constitutive/mohr_coulomb_with_tension_cutoff.cpp

-    StressStrainUtilities::CalculatePrincipalStresses(
-        trial_stress_vector, principal_trial_stress_vector, rotation_matrix);
-
-    if (auto trial_sigma_tau = StressStrainUtilities::TransformPrincipalStressesToSigmaTau(principal_trial_stress_vector);


Nice that the conversion is no longer the callers responsibility!

rfaasse · 2026-02-09T14:07:32Z

applications/GeoMechanicsApplication/custom_constitutive/tension_cutoff.cpp

+    // Note that we attempt to calculate the principal stress vector (which consists of three values) from the given
+    // sigma and tau (two values). As a result, the second principal stress cannot be uniquely determined.
+    const auto principal_stress_vector = Geo::PrincipalStresses{};
+    return YieldFunctionValue(StressStrainUtilities::TransformSigmaTauToPrincipalStresses(


This means that the second component of the 'new' principal stress would be a 0.0. I don't know enough to determine whether that would lead to a problem or not, so I leave the question here without knowing if it's a redundant one (I see that it's the same as before, with an extra comment, so the situation certainly didn't get worse, just trying to understand).

rfaasse · 2026-02-09T14:11:04Z

applications/GeoMechanicsApplication/custom_utilities/stress_strain_utilities.cpp

+    auto eigen_vectors      = Matrix{};
+    CalculatePrincipalStresses(rStressVector, principal_stresses, eigen_vectors);
+
+    return std::make_pair(Geo::PrincipalStresses{principal_stresses}, eigen_vectors);


It's good to keep in mind we create an extra new object here (the PrincipalStresses object from the Vector), so it might impact performance (since this function is probably called for every integration point in a simulation). We don't need to change anything now, since we first want to get it working, but it's something to keep in mind

rfaasse · 2026-02-09T14:12:48Z

applications/GeoMechanicsApplication/tests/cpp_tests/custom_constitutive/test_yield_surface.cpp

+    ::testing::Values(std::make_tuple(Geo::PrincipalStresses{3.0, 2.0, 1.0}, 1.0),
+                      std::make_tuple(Geo::PrincipalStresses{1.7071067811865475, 1.0, 0.2928932188134525}, 0.0),
+                      std::make_tuple(Geo::PrincipalStresses{0.1715728752538099, -1.0, -1.8284271247461901}, -1.0)));


rfaasse · 2026-02-09T14:19:20Z

applications/GeoMechanicsApplication/tests/cpp_tests/custom_constitutive/test_yield_surface.cpp

-    auto p_coulomb_yield_surface =
-        std::unique_ptr<YieldSurface>{std::make_unique<CoulombYieldSurface>(material_properties)};
-    auto serializer = StreamSerializer{};
+    const auto coulomb_yield_surface         = CoulombYieldSurface{material_properties};


I think this test tests a bit less now (we test saving a concrete implementation, instead of saving a pointer to an abstract base with a concrete implementation behind it). I think this could be fine, since saving unique pointers is tested in core (probably?), just wanted to double-check if this was indeed your intention.

Or does it have to do with the idea to move to something less 'inheritance-based' in the future?

avdg81 added 17 commits February 5, 2026 10:43

Copied changes from prototype as starting point

8b1c95e

Removed a member function that was no longer needed

30b0eea

Avoid two unnecessary copies

27d69ca

Renamed a stress-strain utility function

0d02d54

The new name makes clearer what it returns.

Use move operations whenever we can

6c8579f

Apparently, Matrix objects cannot be moved

cb63692

This resolves an issue found by SonarQube.

Simplified a return statement by using a different constructor

aac5ee5

Added a comment that says some API will be removed soon

56db878

Renamed a unit test for consistency

a6d3493

Use functionality provided by UblasUtilities

b5c3136

Added a unit test for utility function `CalculatePrincipalStressesAnd…

40ea077

…RotationMatrix`

Cover the new overload of TensionCutoff::YieldFunctionValue

3955ada

Also renamed some variables and made a few ones `const`.

Parameterized a unit test

12288e8

Added a parameterized unit test for `TensionCutoff::DerivativeOfFlowF…

13fd63e

…unction`

Parameterized another unit test

0415f4a

Renamed a parameterized test fixture

3151ee1

Removed an overload that was no longer needed

e62d25c

avdg81 requested review from WPK4FEM, mnabideltares and rfaasse February 6, 2026 16:08

avdg81 self-assigned this Feb 6, 2026

avdg81 requested a review from a team as a code owner February 6, 2026 16:08

avdg81 added the GeoMechanics Issues related to the GeoMechanicsApplication label Feb 6, 2026

avdg81 added this to GeoMechanicsApplication (Deltares) Feb 6, 2026

avdg81 linked an issue Feb 6, 2026 that may be closed by this pull request

[GeoMechanicsApplication] Modify the Coulomb model such that it accepts sigma and tau as well as principal stresses #14180

Open

github-project-automation bot moved this to In progress in GeoMechanicsApplication (Deltares) Feb 6, 2026

Merge remote-tracking branch 'origin/master' into geo/14180-improve-t…

1adc97f

…he-coulomb-model-framework

WPK4FEM reviewed Feb 9, 2026

View reviewed changes

Removed the prefix geo_ from two files

d896d4f

rfaasse reviewed Feb 9, 2026

View reviewed changes

		#include "custom_constitutive/geo_sigma_tau.hpp"
		#include "custom_constitutive/principal_stresses.hpp"


		void CoulombYieldSurface::SetKappa(double kappa) { mKappa = kappa; }

		// At some point in time we would like to get rid of this API. For now, just forward the request.

Conversation

avdg81 commented Feb 6, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

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!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

rfaasse 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!

avdg81 Feb 9, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

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!

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

avdg81 commented Feb 6, 2026 •

edited

Loading

avdg81 Feb 9, 2026 •

edited

Loading