Skip to content

[Mohr-Coulomb] Add corner return for simultaneous tension-shear yielding#116

Open
Curiim wants to merge 9 commits into
masterfrom
constitutive/mohr_coulomb
Open

[Mohr-Coulomb] Add corner return for simultaneous tension-shear yielding#116
Curiim wants to merge 9 commits into
masterfrom
constitutive/mohr_coulomb

Conversation

@Curiim
Copy link
Copy Markdown
Collaborator

@Curiim Curiim commented Feb 25, 2026

Summary

This PR extracts only the theory-critical changes in Mohr-Coulomb stress integration.

Changes

  • Fix yield admissibility tolerance in compute_yield_state (1e-7)
  • Refactor compute_stress into:
    • elastic predictor
    • single-surface return
    • corner return (multi-surface)
  • Add dedicated helpers:
    • compute_single_surface_return
    • compute_corner_return
    • update_softening_parameters

@Curiim Curiim added the enhancement New feature or request label Feb 25, 2026
@bodhinandach bodhinandach changed the title Mohr-Coulomb: add corner return for simultaneous tension-shear yielding [Mohr-Coulomb] Add corner return for simultaneous tension-shear yielding Feb 25, 2026
Copy link
Copy Markdown
Member

@bodhinandach bodhinandach left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm blocking this for now, but I think we should take a closer look at this, @ningyiw01.

Copy link
Copy Markdown
Member

@bodhinandach bodhinandach left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Something to check

yield_function(0) <= -Tolerance) {
(*state_vars).at("yield_state") = 1; // Shear
} else {
(*state_vars).at("yield_state") = 0; // Elastic
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should this not be possible?

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Talked this through with Ningyi a bit. I think the purpose is just to finalize the yield state on that final iteration of the loop depending on which surface we're closest too. Arguably this should be accomplished as a part of the compute_yield_state functionality. Regardless, I'd think that we shouldn't ever be able to return an elastic yield state if we need to enter this iteration loop.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Better to change this to a corner state

// A. Re-evaluate invariants and active surface at current stress.
// ---------------------------------------------------------------------
this->compute_stress_invariants(current_stress, state_vars);
yield_type = this->compute_yield_state(&yield_function, (*state_vars));
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In the first step, this will be the same as Line 373-377, right? Maybe we should just evaluate this at the very end of the iteration instead?

Comment on lines +443 to +447
if (yield_function(0) > yield_function(1)) {
yield_type = mpm::mohrcoulomb::FailureState::Tensile;
} else {
yield_type = mpm::mohrcoulomb::FailureState::Shear;
}
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do we need to assign (*state_vars).at("yield_state") here?

Comment on lines +437 to +439
if (corner_success) {
stress_correction = corner_correction;
dpdstrain_inc = corner_dpdstrain;
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

what is the suitable (*state_vars).at("yield_state") here?

Copy link
Copy Markdown
Collaborator

@cgeudeker cgeudeker left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Quite a discrepancy in evaluated stresses in the test cases right now. Often on the order of 25% increases in magnitude, so it seems likely that there may be some bug here. May also be best to limit the determination of the yield state to the associated compute_yield_state function and have any following if else statements to be determined by that evaluation.

const mpm::dense_map& state_vars) {
// Tolerance for yield function
const double Tolerance = -1E-1;
const double Tolerance = 1E-7;
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This reduction leads to issues with identifying the expected yield state for some of the test cases. Reverting back to the previous value fixes the issue.

yield_function(0) <= -Tolerance) {
(*state_vars).at("yield_state") = 1; // Shear
} else {
(*state_vars).at("yield_state") = 0; // Elastic
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Talked this through with Ningyi a bit. I think the purpose is just to finalize the yield state on that final iteration of the loop depending on which surface we're closest too. Arguably this should be accomplished as a part of the compute_yield_state functionality. Regardless, I'd think that we shouldn't ever be able to return an elastic yield state if we need to enter this iteration loop.

// ---------------------------------------------------------------------
// B. Convergence check
// ---------------------------------------------------------------------
if (yield_function(0) <= Tolerance && yield_function(1) <= Tolerance) {
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should be an absolute value to avoid overshooting.

yield_function(0) <= -Tolerance) {
(*state_vars).at("yield_state") = 1; // Shear
} else {
(*state_vars).at("yield_state") = 0; // Elastic
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Better to change this to a corner state

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

enhancement New feature or request

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants