[GeoMechanicsApplication] Implement hardening process for Mohr-Coulomb, based on the modified class structure #13965
Conversation
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
…tructure-for-hardening
avdg81
reviewed
Nov 12, 2025
Contributor
avdg81
left a comment
avdg81
left a comment
There was a problem hiding this comment.
Hi Mohamed,
Thank you very much for preparing this clear and well-scoped PR. It was very much what I had expected. I have a few suggestions about the iteration loop as well as moving some functions around, but all of them are relatively minor things. Well done!
applications/GeoMechanicsApplication/custom_constitutive/coulomb_with_tension_cut_off_impl.cpp
Show resolved
Hide resolved
applications/GeoMechanicsApplication/custom_constitutive/coulomb_with_tension_cut_off_impl.cpp
Outdated
Show resolved
Hide resolved
applications/GeoMechanicsApplication/custom_constitutive/coulomb_with_tension_cut_off_impl.cpp
Outdated
Show resolved
Hide resolved
applications/GeoMechanicsApplication/custom_constitutive/coulomb_with_tension_cut_off_impl.cpp
Show resolved
Hide resolved
applications/GeoMechanicsApplication/custom_constitutive/coulomb_with_tension_cut_off_impl.cpp
Outdated
Show resolved
Hide resolved
...icsApplication/tests/cpp_tests/custom_constitutive/test_mohr_coulomb_with_tension_cutoff.cpp
Outdated
Show resolved
Hide resolved
...icsApplication/tests/cpp_tests/custom_constitutive/test_mohr_coulomb_with_tension_cutoff.cpp
Outdated
Show resolved
Hide resolved
...icsApplication/tests/cpp_tests/custom_constitutive/test_mohr_coulomb_with_tension_cutoff.cpp
Outdated
Show resolved
Hide resolved
...icsApplication/tests/cpp_tests/custom_constitutive/test_mohr_coulomb_with_tension_cutoff.cpp
Outdated
Show resolved
Hide resolved
...icsApplication/tests/cpp_tests/custom_constitutive/test_mohr_coulomb_with_tension_cutoff.cpp
Show resolved
Hide resolved
…TensionCutOffImpl`
…thTensionCutOffImpl`
…WithTensionCutOffImpl`
…mbWithTensionCutOffImpl`
…hTensionCutOffImpl`
… multiplier rather than passing it through a function parameter
It now uses proper mathematical structures rather than working things out at the component level.
… Coulomb yield surface to the Mohr-Coulomb implementation class
- Simplified some code by using a utility function that creates vectors. - Made member `CheckMaterialProperties` `const`. - Extracted a private member function that checks hardening coefficients (to reduce duplication).
avdg81
reviewed
Nov 28, 2025
Contributor
avdg81
left a comment
avdg81
left a comment
There was a problem hiding this comment.
Hi Mohamed,
Thanks a lot for processing the review suggestions. After reviewing your latest changes, I have a few more (minor) suggestions. I believe we're almost there. Thanks for all the work that you have done.
Comment on lines
88
to
99
| // Check whether the tension cut-off lies beyond the apex | ||
| auto result = Vector{ZeroVector(2)}; | ||
| result[0] = mCoulombYieldSurface.CalculateApex(); | ||
| const auto tensile_strength = mTensionCutOff.GetTensileStrength(); | ||
| if (tensile_strength > result[0]) return result; | ||
|
|
||
| const auto cohesion = mCoulombYieldSurface.GetCohesion(); | ||
| const auto sin_phi = std::sin(mCoulombYieldSurface.GetFrictionAngleInRadians()); | ||
| const auto cos_phi = std::cos(mCoulombYieldSurface.GetFrictionAngleInRadians()); | ||
| result[0] = (tensile_strength - cohesion * cos_phi) / (1.0 - sin_phi); | ||
| result[1] = (cohesion * cos_phi - tensile_strength * sin_phi) / (1.0 - sin_phi); | ||
| return result; |
Contributor
There was a problem hiding this comment.
By using UblasUtilities::CreateVector we can make the code slightly more elegant:
Suggested change
| // Check whether the tension cut-off lies beyond the apex | |
| auto result = Vector{ZeroVector(2)}; | |
| result[0] = mCoulombYieldSurface.CalculateApex(); | |
| const auto tensile_strength = mTensionCutOff.GetTensileStrength(); | |
| if (tensile_strength > result[0]) return result; | |
| const auto cohesion = mCoulombYieldSurface.GetCohesion(); | |
| const auto sin_phi = std::sin(mCoulombYieldSurface.GetFrictionAngleInRadians()); | |
| const auto cos_phi = std::cos(mCoulombYieldSurface.GetFrictionAngleInRadians()); | |
| result[0] = (tensile_strength - cohesion * cos_phi) / (1.0 - sin_phi); | |
| result[1] = (cohesion * cos_phi - tensile_strength * sin_phi) / (1.0 - sin_phi); | |
| return result; | |
| const auto tensile_strength = mTensionCutOff.GetTensileStrength(); | |
| if (const auto apex = mCoulombYieldSurface.CalculateApex(); tensile_strength > apex) | |
| return UblasUtilities::CreateVector({apex, 0.0}); | |
| const auto cohesion = mCoulombYieldSurface.GetCohesion(); | |
| const auto sin_phi = std::sin(mCoulombYieldSurface.GetFrictionAngleInRadians()); | |
| const auto cos_phi = std::cos(mCoulombYieldSurface.GetFrictionAngleInRadians()); | |
| return UblasUtilities::CreateVector({(tensile_strength - cohesion * cos_phi) / (1.0 - sin_phi), | |
| (cohesion * cos_phi - tensile_strength * sin_phi) / (1.0 - sin_phi)}); |
Note that the early exit now looks consistent with the result returned at line 129.
| CoulombYieldSurface::KappaDependentFunction MakeLinearFunction(double Value, double Coefficient) | ||
| { | ||
| return [Value, Coefficient](double kappa) { return std::max(Value + Coefficient * kappa, 0.0); }; | ||
| return [Value, Coefficient](double kappa) { return Value + Coefficient * kappa; }; |
Contributor
There was a problem hiding this comment.
Nitpicking (to follow the Kratos Style Guide):
Suggested change
| return [Value, Coefficient](double kappa) { return Value + Coefficient * kappa; }; | |
| return [Value, Coefficient](double Kappa) { return Value + Coefficient * Kappa; }; |
| mMaterialProperties[GEO_COULOMB_HARDENING_CONVERGENCE_TOLERANCE] = 1.0e-8; | ||
| } | ||
|
|
||
| this->CheckMaterialProperties(); |
Contributor
There was a problem hiding this comment.
No this pointer required here, so please leave it out.
Suggested change
| this->CheckMaterialProperties(); | |
| CheckMaterialProperties(); |
Comment on lines
193
to
194
| const auto mean = | ||
| std::accumulate(principal_stress_vector.begin(), principal_stress_vector.end(), 0.0) / 3.0; |
Contributor
There was a problem hiding this comment.
We could consider to replace the hard-coded size by asking for the vector's size:
Suggested change
| const auto mean = | |
| std::accumulate(principal_stress_vector.begin(), principal_stress_vector.end(), 0.0) / 3.0; | |
| const auto mean = std::accumulate(principal_stress_vector.begin(), principal_stress_vector.end(), 0.0) / | |
| static_cast<double>(principal_stress_vector.size()); |
What do you think? And perhaps @WPK4FEM may also want to weigh in?
Comment on lines
210
to
232
| check_properties.CheckAvailability(GEO_COHESION_FUNCTION_COEFFICIENTS); | ||
| const auto& cohesion_vec = mMaterialProperties[GEO_COHESION_FUNCTION_COEFFICIENTS]; | ||
| for (unsigned int i = 0; i < cohesion_vec.size(); ++i) { | ||
| KRATOS_ERROR_IF(cohesion_vec[i] < 0.0) | ||
| << "Entry " << i << " in " << GEO_COHESION_FUNCTION_COEFFICIENTS.Name() | ||
| << " out of range. Value: " << cohesion_vec[i]; | ||
| } | ||
|
|
||
| check_properties.CheckAvailability(GEO_FRICTION_ANGLE_FUNCTION_COEFFICIENTS); | ||
| const auto& friction_vec = mMaterialProperties[GEO_FRICTION_ANGLE_FUNCTION_COEFFICIENTS]; | ||
| for (unsigned int i = 0; i < friction_vec.size(); ++i) { | ||
| KRATOS_ERROR_IF(friction_vec[i] < 0.0) | ||
| << "Entry " << i << " in " << GEO_FRICTION_ANGLE_FUNCTION_COEFFICIENTS.Name() | ||
| << " out of range. Value: " << friction_vec[i]; | ||
| } | ||
|
|
||
| check_properties.CheckAvailability(GEO_DILATANCY_ANGLE_FUNCTION_COEFFICIENTS); | ||
| const auto& dilatancy_vec = mMaterialProperties[GEO_DILATANCY_ANGLE_FUNCTION_COEFFICIENTS]; | ||
| for (unsigned int i = 0; i < dilatancy_vec.size(); ++i) { | ||
| KRATOS_ERROR_IF(dilatancy_vec[i] < 0.0) | ||
| << "Entry " << i << " in " << GEO_DILATANCY_ANGLE_FUNCTION_COEFFICIENTS.Name() | ||
| << " out of range. Value: " << dilatancy_vec[i]; | ||
| } |
Contributor
There was a problem hiding this comment.
The checks of the three coefficient vectors look very similar. Perhaps we can extract it to a private member? I was thinking of something like:
void CoulombYieldSurface::CheckHardeningCoefficients(const Variable<Vector>& rCoefficientsVariable,
const CheckProperties& rChecker) const
{
rChecker.CheckAvailability(rCoefficientsVariable);
const auto& coefficients = mMaterialProperties[rCoefficientsVariable];
for (auto i = std::size_t{0}; i < coefficients.size(); ++i) {
KRATOS_ERROR_IF(coefficients[i] < 0.0) << "Entry " << i << " in " << rCoefficientsVariable.Name()
<< " out of range. Value: " << coefficients[i];
}
}Once that member function is available, this piece of code can be reduced to three function calls:
Suggested change
| check_properties.CheckAvailability(GEO_COHESION_FUNCTION_COEFFICIENTS); | |
| const auto& cohesion_vec = mMaterialProperties[GEO_COHESION_FUNCTION_COEFFICIENTS]; | |
| for (unsigned int i = 0; i < cohesion_vec.size(); ++i) { | |
| KRATOS_ERROR_IF(cohesion_vec[i] < 0.0) | |
| << "Entry " << i << " in " << GEO_COHESION_FUNCTION_COEFFICIENTS.Name() | |
| << " out of range. Value: " << cohesion_vec[i]; | |
| } | |
| check_properties.CheckAvailability(GEO_FRICTION_ANGLE_FUNCTION_COEFFICIENTS); | |
| const auto& friction_vec = mMaterialProperties[GEO_FRICTION_ANGLE_FUNCTION_COEFFICIENTS]; | |
| for (unsigned int i = 0; i < friction_vec.size(); ++i) { | |
| KRATOS_ERROR_IF(friction_vec[i] < 0.0) | |
| << "Entry " << i << " in " << GEO_FRICTION_ANGLE_FUNCTION_COEFFICIENTS.Name() | |
| << " out of range. Value: " << friction_vec[i]; | |
| } | |
| check_properties.CheckAvailability(GEO_DILATANCY_ANGLE_FUNCTION_COEFFICIENTS); | |
| const auto& dilatancy_vec = mMaterialProperties[GEO_DILATANCY_ANGLE_FUNCTION_COEFFICIENTS]; | |
| for (unsigned int i = 0; i < dilatancy_vec.size(); ++i) { | |
| KRATOS_ERROR_IF(dilatancy_vec[i] < 0.0) | |
| << "Entry " << i << " in " << GEO_DILATANCY_ANGLE_FUNCTION_COEFFICIENTS.Name() | |
| << " out of range. Value: " << dilatancy_vec[i]; | |
| } | |
| CheckHardeningCoefficients(GEO_COHESION_FUNCTION_COEFFICIENTS, check_properties); | |
| CheckHardeningCoefficients(GEO_FRICTION_ANGLE_FUNCTION_COEFFICIENTS, check_properties); | |
| CheckHardeningCoefficients(GEO_DILATANCY_ANGLE_FUNCTION_COEFFICIENTS, check_properties); |
| CalculatePlasticMultiplier(rSigmaTau, DerivativeOfFlowFunction(rSigmaTau, AveragingType)); | ||
| } | ||
|
|
||
| void CoulombYieldSurface::CheckMaterialProperties() |
Contributor
There was a problem hiding this comment.
When member functions can be made const, please do so:
Suggested change
| void CoulombYieldSurface::CheckMaterialProperties() | |
| void CoulombYieldSurface::CheckMaterialProperties() const |
|
|
||
| private: | ||
| void InitializeKappaDependentFunctions(); | ||
| void CheckMaterialProperties(); |
Contributor
There was a problem hiding this comment.
As a result of my previous suggestion:
Suggested change
| void CheckMaterialProperties(); | |
| void CheckMaterialProperties() const; |
Contributor
There was a problem hiding this comment.
Important note: this suggestion is unrelated to your PR, but I just happened to see it: member function InterfaceCoulombWithTensionCutOff::WorkingSpaceDimension apparently assumes a line interface element (see the comment at line 63). This needs to be fixed (including a (unit?) test) in a separate PR.
Comment on lines
141
to
143
| double kappa_start = -1.0; | ||
| mapped_sigma_tau = mCoulombWithTensionCutOffImpl.DoReturnMapping( | ||
| trial_sigma_tau, CoulombYieldSurface::CoulombAveragingType::NO_AVERAGING, kappa_start); |
Contributor
There was a problem hiding this comment.
I would like to propose a slightly different approach where we can tell mCoulombWithTensionCutOffImpl to save geo/experimental-changes-for-coulomb-hardening.
Comment on lines
-72
to
-77
| KRATOS_EXPECT_EXCEPTION_IS_THROWN( | ||
| [[maybe_unused]] const auto unused = law.Check(properties, element_geometry, process_info), | ||
| "GEO_COHESION does not exist in the property with Id 3.") | ||
| properties.SetValue(GEO_COHESION, -1.0); | ||
| KRATOS_EXPECT_EXCEPTION_IS_THROWN( | ||
| [[maybe_unused]] const auto unused = law.Check(properties, element_geometry, process_info), "GEO_COHESION in the property with Id 3 has an invalid value: -1 is out of the range [0, -).") |
Contributor
There was a problem hiding this comment.
Could you please try to reproduce these tests for the Coulomb yield surface?
avdg81
approved these changes
Dec 3, 2025
Contributor
avdg81
left a comment
avdg81
left a comment
There was a problem hiding this comment.
Hi Mohamed,
Thanks for taking care of the review suggestions. I feel this PR is ready to go.
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
📝 Description
The code is copied and modified to match the recently implemented structure for Mohr-Coulomb
🆕 Changelog