[GeoMechanicsApplication] Create a framework for a stress-dependent Youngs' modulus by markelov208 · Pull Request #14656 · KratosMultiphysics/Kratos

markelov208 · 2026-08-06T10:05:44Z

Added new GeoMechanicsApplication variables and registration for configuring Young’s modulus formulations (e.g., Eur) and related parameters
Updated GeoIncrementalLinearElasticLaw to compute the elastic matrix via utilities and to optionally adjust Young’s modulus based on the chosen formulation.
Added unit tests to cover new Eur functionalities
Added an integration test for Eur formulation
Updated README with Eur description

Copilot

Pull request overview

This PR (currently marked draft) introduces the initial framework to support a stress-dependent Young’s modulus formulation in the GeoMechanics incremental linear elastic constitutive law, including new application-level variables and test updates to accommodate the updated material-property requirements.

Changes:

Added new GeoMechanicsApplication variables and registration for configuring Young’s modulus formulations (e.g., Eur) and related parameters.
Updated GeoIncrementalLinearElasticLaw to compute the elastic matrix via utilities and to optionally adjust Young’s modulus based on the chosen formulation.
Updated incremental linear elastic law unit tests to set explicit material properties for initialization paths.

Reviewed changes

Copilot reviewed 7 out of 7 changed files in this pull request and generated 6 comments.

Show a summary per file

File	Description
applications/GeoMechanicsApplication/tests/cpp_tests/custom_constitutive/test_incremental_linear_elastic_plane_strain_2D_law.cpp	Updates tests to set material properties before initialization.
applications/GeoMechanicsApplication/tests/cpp_tests/custom_constitutive/test_incremental_linear_elastic_3D_law.cpp	Updates tests to set material properties before initialization; contains a small redundancy.
applications/GeoMechanicsApplication/geo_mechanics_application.cpp	Registers new variables for Young’s modulus formulation framework.
applications/GeoMechanicsApplication/geo_mechanics_application_variables.h	Defines new application variables for the formulation and its parameters.
applications/GeoMechanicsApplication/geo_mechanics_application_variables.cpp	Creates new application variables for runtime use/registration.
applications/GeoMechanicsApplication/custom_constitutive/incremental_linear_elastic_law.h	Introduces policy concept and new APIs/state to support stress-dependent modulus.
applications/GeoMechanicsApplication/custom_constitutive/incremental_linear_elastic_law.cpp	Implements formulation selection, Eur computation, and updated elastic-matrix computation path.

Suppressed comments (4)

applications/GeoMechanicsApplication/custom_constitutive/incremental_linear_elastic_law.cpp:271

mPolicy affects the computed Young’s modulus but is not serialized. After restart/load, the constructor defaults it to Constant, while mIsModelInitialized may load as true, so InitializePolicy() may never run again and Eur behavior would be silently lost. Consider serializing the active policy (e.g., by variant index) alongside the other state.

void GeoIncrementalLinearElasticLaw::save(Serializer& rSerializer) const
{
    KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, GeoLinearElasticLaw)
    rSerializer.save("ConstitutiveDimension"s, mpConstitutiveDimension);
    rSerializer.save("StressVector"s, mStressVector);
    rSerializer.save("StressVectorFinalized"s, mStressVectorFinalized);
    rSerializer.save("DeltaStrainVector"s, mDeltaStrainVector);
    rSerializer.save("StrainVectorFinalized"s, mStrainVectorFinalized);
    rSerializer.save("IsModelInitialized"s, mIsModelInitialized);
}

applications/GeoMechanicsApplication/tests/cpp_tests/custom_constitutive/test_incremental_linear_elastic_plane_strain_2D_law.cpp:178

The new material Properties set only YOUNG_MODULUS, but the law now requires POISSON_RATIO and expects a drainage type. Set POISSON_RATIO and GEO_DRAINAGE_TYPE here as well to avoid runtime errors.

    auto p_material_props              = Kratos::make_shared<Kratos::Properties>(0);
    (*p_material_props)[YOUNG_MODULUS] = 30000000.0;
    initial_parameters.SetMaterialProperties(*p_material_props);

applications/GeoMechanicsApplication/tests/cpp_tests/custom_constitutive/test_incremental_linear_elastic_3D_law.cpp:157

Same issue as above: the newly added Properties only set YOUNG_MODULUS, but POISSON_RATIO and GEO_DRAINAGE_TYPE are required by the updated law/utilities. Add them here to prevent runtime errors.

    auto p_material_props              = Kratos::make_shared<Kratos::Properties>(0);
    (*p_material_props)[YOUNG_MODULUS] = 30000000.0;
    initial_parameters.SetMaterialProperties(*p_material_props);

applications/GeoMechanicsApplication/tests/cpp_tests/custom_constitutive/test_incremental_linear_elastic_3D_law.cpp:187

In this test, the new Properties only set YOUNG_MODULUS, but the updated law expects POISSON_RATIO and GEO_DRAINAGE_TYPE to be present when computing elastic properties. Add them here as well.

    auto p_material_props              = Kratos::make_shared<Kratos::Properties>(0);
    (*p_material_props)[YOUNG_MODULUS] = 30000000.0;
    initial_parameters.SetMaterialProperties(*p_material_props);

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

+    auto p_material_props              = Kratos::make_shared<Kratos::Properties>(0);
+    (*p_material_props)[YOUNG_MODULUS] = 30000000.0;
+    parameters.SetMaterialProperties(*p_material_props);


+    auto p_material_props              = Kratos::make_shared<Kratos::Properties>(0);
+    (*p_material_props)[YOUNG_MODULUS] = 30000000.0;
+    parameters.SetMaterialProperties(*p_material_props);


…ents

Copilot

Pull request overview

Copilot reviewed 7 out of 7 changed files in this pull request and generated no new comments.

Suppressed comments (2)

applications/GeoMechanicsApplication/tests/cpp_tests/custom_constitutive/test_incremental_linear_elastic_3D_law.cpp:183

Duplicate SetStressVector call; the second call is redundant and could hide copy/paste mistakes in test setup.

    initial_parameters.SetStressVector(initial_stress);
    initial_parameters.SetStressVector(initial_stress);

applications/GeoMechanicsApplication/custom_constitutive/incremental_linear_elastic_law.cpp:21

This file uses std::decay_t and std::is_same_v in GetYoungsModulus() but does not include <type_traits>. Relying on transitive includes can break builds on different standard library implementations.

#include <cmath>
#include <limits>
#include <string>

rfaasse

Hi Gennady, thank you for the hard work on this very useful new framework! It is well tested and easy-to-read!

I do think these changes have some differences with respect to the design you proposed via e-mail, because of which we do not fully utilize the extensibility and generic applicability of the formulations + std::variant pattern. Most of the comments are regarding this point and are hopefully getting this point across clearly! If not, do not hesitate to ask for a meeting, because that might be more efficient 👍

rfaasse · 2026-08-13T12:44:51Z

    KRATOS_CATCH("")
 }

+double GeoIncrementalLinearElasticLaw::GetYoungsModulus(const Properties& rProperties, double YoungsModulus) const


This lambda function within the std::visit is different than I'd expect. In the option you sent us via mail when looking at the design, it looked like:

double GetYoungModulus(const Properties& rProps, const Vector& stress) const { return std::visit([&](auto&& policy){ return policy(rProps, stress); }, mPolicy); }

This only invokes the call operator on the policy (or formulation, which is now how it's called) without having to know the exact type of the formulation (i.e. the law is 'formulation-agnostic' and only cares about being able to call this function).

To do this (see also other comments) and fully utilize this pattern with the formulations and std::variant, I'd propose to move the implementation (this will need some small adjustments in the function itself) of CalculateYoungsModulusForEur to the call operator of the Eur formulation

rfaasse · 2026-08-17T07:07:04Z

+namespace Formulations
+{
+struct Constant {
+    static constexpr const char* Name = "Constant";


Although we can't make a std::string constexpr yet, I'd still prefer to use modern types over c-style char*, curious what you think!

char is old and it works. I tried std::string_view but serializer does not support it. It is getting unnecessary complex.

rfaasse · 2026-08-17T07:10:49Z

+{
+struct Constant {
+    static constexpr const char* Name = "Constant";
+    double                       operator()(const Properties&, double) const;


These operators now don't have an implementation and are not used. That needs to change if we want to make full use of this std::variant pattern. For the Constant variation, the implementation is probably simple (something like { return rProperties[YOUNG_MODULUS]; }). For the Eur formulation, I'd expect a function similar to the CalculateYoungsModulusForEur function which now is only available in the incremental linear elastic law.

rfaasse · 2026-08-17T07:29:02Z

+
+struct Eur {
+    static constexpr const char* Name = "Eur";
+    double                       operator()(const Properties&, double) const;


Since these formulations could be stress dependent (like this one), I would expect the call operator (also of the constant version) to have the following two arguments:

Suggested change

double operator()(const Properties&, double) const;

double operator()(const Properties&, const Vector&) const;

This also allows for things like position dependent youngs modulus later (since that could also be captured in a Vector). The Young's modulus, which we pass via the second argument now can be retrieved from the properties, so it's not needed to add a second argument there.

This signature is also more in line with the proposal you sent via e-mail (the re-name to formulation is not a problem in my opinion):

#include <variant> namespace Policies { struct Constant { double operator()(const Properties&, const Vector&) const; }; struct Eur { double operator()(const Properties&, const Vector&) const; }; struct E50{ double operator()(const Properties&, const Vector&) const; }; struct Esomething { double operator()(const Properties&, const Vector&) const; }; }

rfaasse · 2026-08-17T07:30:54Z

+
+struct Eur {
+    static constexpr const char* Name = "Eur";
+    double                       operator()(const Properties&, double) const;


When the implementation of this calculation is added here, I think it makes sense to move the 'formulations' namespace to its own files.

rfaasse · 2026-08-17T12:01:45Z

+
+E = E_ur^ref ((-sigma_3') / p_ref)^m
+
+using E_ur^ref = 1.0e7, p_ref = 50 and m = 1, the analytical axial stage strain for the stress change between -100 and -50 kPa is:


Could you add units to E and p_ref?

added Pa and changed kPa on Pa.

rfaasse · 2026-08-17T12:02:14Z

+
+using E_ur^ref = 1.0e7, p_ref = 50 and m = 1, the analytical axial stage strain for the stress change between -100 and -50 kPa is:
+
+epsilon_yy = (p_ref / E_ur^ref) ln(2)


Nice that we have an analytical solution here

rfaasse · 2026-08-17T12:04:15Z

+
+## Stage 1 - Unloading
+
+The top normal stress is reduced from -100 kPa to -50 kPa.


In this test we don't assert the displacements, only the stresses, so I think it's okay to leave everything in kPa instead of Pa. I'm not sure if it would work correctly otherwise (but I might be wrong there, maybe @WPK4FEM knows?)

Good catch! The test also asserts displacements. Changed kPa on Pa.

rfaasse · 2026-08-17T12:07:08Z

+                    "value": [0.0, 0.0],
+                    "table": [1, 0],


Since table 1 is constant, I think we could use the value and disregard the table here:

Suggested change

"value": [0.0, 0.0],

"table": [1, 0],

"value": [100.0, 0.0],

"table": [0, 0],

There could be an issue with this though (see #14642, at least we saw the issue for prescribed pressures, not sure if it'll influence prescribed stress), so if the results change we can just revert it and leave it for later.

removed table 1. This test is not sensitive to use a table or a value.

rfaasse · 2026-08-17T12:08:54Z

+                    "variable_name": "NORMAL_CONTACT_STRESS",
+                    "active": [true, false],
+                    "value": [0.0, 0.0],
+                    "table": [1, 0],


See previous comment about table 1

markelov208

Hi Richard, thank you very much for the review. Most of comments are answered except comments related to std::variant.

markelov208 · 2026-08-18T10:05:03Z

+    const auto friction_angle_rad = ConstitutiveLawUtilities::GetFrictionAngleInRadians(rProperties);
+    const auto stress_shift =
+        rProperties[GEO_COHESION] * std::cos(friction_angle_rad) / std::sin(friction_angle_rad);


Currently, they are required and checked. if they do not exist then the error is thrown.

markelov208 · 2026-08-18T10:07:24Z

 * $`C`$: The elastic constitutive tensor.
 * $`\Delta \Delta u`$: The incremental relative displacement vector.

+### 1.2 Incremental linear elastic E_ur law


re-phrased the title and beginning. now it is described as an option.

markelov208 · 2026-08-18T10:27:12Z

+
+where $c$ is `GEO_COHESION` and $\phi$ is `GEO_FRICTION_ANGLE`.
+
+For numerical robustness, the implemented formulation applies lower bounds to the confinement terms.


There is a check for a non-positive base for the power operation. I removed this sentence.

markelov208 · 2026-08-18T10:29:27Z

+    auto p_material_props              = Kratos::make_shared<Kratos::Properties>(0);
+    (*p_material_props)[YOUNG_MODULUS] = 30000000.0;
+    parameters.SetMaterialProperties(*p_material_props);


Now InitializeMaterialResponseCauchy needs any properties to set up the formulation. Currently, I extracted properties from Calculate3DStress to avoid this extra settings.

markelov208 · 2026-08-18T10:30:04Z


+### 1.2 Incremental linear elastic E_ur law
+
+This law is an incremental linear elastic continuum model with stress-dependent stiffness. It is an option in the linear elastic material model. To activate it, `GEO_YOUNGS_MODULUS_FORMULATION` shall be added to material properties. It has to have `Eur` value. This option is implemented using `std::variant` that combines the flexibility of Strategy with the performance of stack allocation. It avoids both deep inheritance trees and virtual calls during simulation. Another benefit is an easy extension with other formulations like `E50` and `Eoed`.


Added their full names: Secant Stiffness at 50% Strength and Oedometer Modulus

markelov208 · 2026-08-18T12:05:12Z

    parameters.SetStressVector(initial_stress);
+    auto p_material_props              = Kratos::make_shared<Kratos::Properties>(0);
+    (*p_material_props)[YOUNG_MODULUS] = 30000000.0;
+    parameters.SetMaterialProperties(*p_material_props);


after introduction of the formulation, InitializeMaterialResponseCauchy needs properties, which were not available here.

markelov208 · 2026-08-18T12:06:10Z

+        "displacement_relative_tolerance": 1.0e-06,
+        "displacement_absolute_tolerance": 1.0e-09,
+        "residual_relative_tolerance": 1.0e-04,
+        "residual_absolute_tolerance": 1.0e-09,


I use that script but it does not cover all possible cases yet. removed residual*

markelov208 · 2026-08-18T12:06:40Z

+
+E = E_ur^ref ((-sigma_3') / p_ref)^m
+
+using E_ur^ref = 1.0e7, p_ref = 50 and m = 1, the analytical axial stage strain for the stress change between -100 and -50 kPa is:


added Pa and changed kPa on Pa.

markelov208 · 2026-08-18T12:07:29Z

+
+## Stage 1 - Unloading
+
+The top normal stress is reduced from -100 kPa to -50 kPa.


Good catch! The test also asserts displacements. Changed kPa on Pa.

markelov208 · 2026-08-18T12:08:13Z

+                    "value": [0.0, 0.0],
+                    "table": [1, 0],


removed table 1. This test is not sensitive to use a table or a value.

Copilot

Pull request overview

Copilot reviewed 17 out of 17 changed files in this pull request and generated 1 comment.

Suppressed comments (3)

applications/GeoMechanicsApplication/custom_constitutive/README.md:77

Missing space in the required-inputs bullet makes the rendered docs harder to read.

- `GEO_YOUNGS_MODULUS_FORMULATION`with "Eur" value

applications/GeoMechanicsApplication/custom_constitutive/README.md:66

This paragraph says GEO_YOUNGS_MODULUS_FORMULATION can have only one value (Eur), but the implementation also supports Constant (and defaults to Constant when the property is absent). Please align the documentation with the supported values and clarify that Eur is the stress-dependent option.

The linear elastic law has an option with a stress-dependent Young's modulus. To activate this option, `GEO_YOUNGS_MODULUS_FORMULATION` shall be added to material properties. By now unloading/reloading Stiffness formulation is implemented and the keyword can have only one value `Eur`. This feature is implemented using `std::variant` that combines the flexibility of Strategy with the performance of stack allocation. It avoids both deep inheritance trees and virtual calls during simulation. Another benefit is an easy extension with other formulations like Secant Stiffness at 50% Strength and Oedometer Modulus.

applications/GeoMechanicsApplication/custom_constitutive/incremental_linear_elastic_law.cpp:150

The Eur formulation uses trigonometric functions of GEO_FRICTION_ANGLE, so the Check should also enforce an upper bound (<= 90°), consistent with other constitutive-law checks (e.g., InterfaceCoulombLaw). Currently only a lower bound is enforced via AllExclusive.

        const CheckProperties check_properties(rMaterialProperties, "parameters of material",
                                               CheckProperties::Bounds::AllExclusive);
        check_properties.Check(GEO_PRESSURE_REFERENCE);
        check_properties.Check(GEO_STRESS_DEPENDENCY_EXPONENT);
        check_properties.Check(GEO_COHESION);
        check_properties.Check(GEO_FRICTION_ANGLE);
    }

+YoungsModulusVariant InitializeFormulation(const std::string& rFormulation)
+{
+    if (rFormulation == Constant::Name) return Constant{};
+    if (rFormulation == Eur::Name) return Eur{};
+    KRATOS_ERROR << "Unknown GEO_YOUNGS_MODULUS_FORMULATION: " << rFormulation;
+}


markelov208 added 2 commits August 6, 2026 11:32

added new variables

10ebd09

added variant for const Youngs modulus

fbf7838

markelov208 self-assigned this Aug 6, 2026

markelov208 requested a review from a team as a code owner August 6, 2026 10:05

renames some variables

15e402e

markelov208 linked an issue Aug 6, 2026 that may be closed by this pull request

[GeoMechanicsApplication] Create a framework for a stress-dependent Youngs'modulus #14542

Open

fixed undrained case and 3 code smells

f5f296f

markelov208 requested a lite review from Copilot August 7, 2026 06:29

Copilot started reviewing on behalf of markelov208 August 7, 2026 06:30 View session

Copilot AI reviewed Aug 7, 2026

View reviewed changes

added mPolicy to constructors, serializer, fixed some of Copilot comm…

9d0067b

…ents

markelov208 requested a lite review from Copilot August 7, 2026 10:39

Copilot started reviewing on behalf of markelov208 August 7, 2026 10:40 View session

Copilot AI reviewed Aug 7, 2026

View reviewed changes

markelov208 added 14 commits August 7, 2026 12:48

added type_traits

ee7b63f

removed accidentally copied line

94eebaa

minor cleanup

082b402

added a description to README

db12fac

added unit tests, 5 still fail

b9ba93d

fixed Check unit test

b72a853

corrected expected value

26fa611

corrected expected stress

b99ee59

fixed reset test

e8507c2

commented copy test

13cf25b

used static CreateIncrementalLinearElastic3DLaw

fdd5463

used inline instead of static

315c236

reverted

671d10b

fixed GeoIncrementalLinearElasticEur3DLaw_CopyConstructorCop

9b3d031

markelov208 requested a lite review from Copilot August 9, 2026 06:47

Copilot started reviewing on behalf of markelov208 August 9, 2026 06:48 View session

markelov208 requested review from WPK4FEM and rfaasse August 10, 2026 05:55

markelov208 added 5 commits August 10, 2026 08:51

corrected README

1056681

corrected table numbers

19cc01d

used template instead of auto

1b9b71a

renamed Policy with Formulation

dfa9e2e

updated integration test

49aee46

rfaasse reviewed Aug 17, 2026

View reviewed changes

markelov208 added 8 commits August 18, 2026 04:33

cleaned integration test

a03a5f6

extracted SetBasicPropertiesFor3D

f4eb62e

extracted SetBasicPropertiesFor2D

e32be64

split MultipleScenarios test

085bdbc

removed 3 tests and moved 1 test to other file

a2d9776

renamed CreateConstantYoungsModulusProperties

b267997

corrected README

34995c7

split InputYoungsModulusFormulation test

36bd499

markelov208 commented Aug 18, 2026

View reviewed changes

markelov208 added 6 commits August 18, 2026 17:30

Merge branch 'master' into geo/14542-framework-stress-dependend-youngs

9751655

updated test_element_lab

93f58a6

corrected test_element_lab.py

77d25c2

moved variant staff into separate files

7def7bf

used operator ()

6c68788

added const

5a4cd7b

markelov208 requested a lite review from Copilot August 20, 2026 13:48

Copilot started reviewing on behalf of markelov208 August 20, 2026 13:48 View session

Copilot AI reviewed Aug 20, 2026

View reviewed changes

markelov208 added 4 commits August 20, 2026 16:10

corrected README

5a1954e

moved checks to youngs_modulus_formulations

2026f2a

removed optional

fb535f6

cleaned GetYoungsModulus

18c7681

	double operator()(const Properties&, double) const;
	double operator()(const Properties&, const Vector&) const;


		E = E_ur^ref ((-sigma_3') / p_ref)^m

		using E_ur^ref = 1.0e7, p_ref = 50 and m = 1, the analytical axial stage strain for the stress change between -100 and -50 kPa is:


		using E_ur^ref = 1.0e7, p_ref = 50 and m = 1, the analytical axial stage strain for the stress change between -100 and -50 kPa is:

		epsilon_yy = (p_ref / E_ur^ref) ln(2)


		## Stage 1 - Unloading

		The top normal stress is reduced from -100 kPa to -50 kPa.


		where $c$ is `GEO_COHESION` and $\phi$ is `GEO_FRICTION_ANGLE`.

		For numerical robustness, the implemented formulation applies lower bounds to the confinement terms.


		### 1.2 Incremental linear elastic E_ur law

		This law is an incremental linear elastic continuum model with stress-dependent stiffness. It is an option in the linear elastic material model. To activate it, `GEO_YOUNGS_MODULUS_FORMULATION` shall be added to material properties. It has to have `Eur` value. This option is implemented using `std::variant` that combines the flexibility of Strategy with the performance of stack allocation. It avoids both deep inheritance trees and virtual calls during simulation. Another benefit is an easy extension with other formulations like `E50` and `Eoed`.

Uh oh!

Conversation

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

Uh oh!

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Reviewed changes

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

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!

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!

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

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Uh oh!

Reviewers

Assignees

Labels

Projects

markelov208 commented Aug 6, 2026 •

edited

Loading