From 9c37afa45605005015a73b55a7f95ed9ea02eb12 Mon Sep 17 00:00:00 2001 From: ClaudeBot Date: Thu, 17 Sep 2026 11:34:43 -0400 Subject: [PATCH 1/6] Add calculate_current_stress prototype for MaterialModelsBase.jl#12 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Detailed change: Introduces `calculate_current_stress(m, ϵ, state)` as a prototype for MaterialModelsBase.jl#12: given a strain and an already-converged/current material state, it returns the corresponding stress without invoking any local Newton iteration that would advance history variables. This is meant for postprocessing, e.g. recomputing the stress at a slightly different strain than the one that produced `state`. Implemented for `LinearElastic` (via a generic `NoMaterialState` fallback), `Plastic`, `GeneralizedMaxwell`, and the `RotatedMaterial` wrapper. Also supports lower-dimensional stress/strain via `MaterialModelsBase.ReducedStressState` for `LinearElastic` and `Plastic`, replacing the ad hoc per-material `calculate_stress` dispatch hack used for plane-stress postprocessing in FerriteAssembly.jl#94's mixed_materials tutorial (that tutorial itself is not changed here; a follow-up PR should call this function instead). Went through the dual-review workflow with two rounds of independent Codex review (same thread). The first round caught two real bugs before implementation: (1) the reduced-dimensional `Plastic` formula originally reduced the plastic strain before subtracting it from the strain, which discards its out-of-plane component and gives the wrong answer for `PlaneStrain` (fixed by expanding the total strain to 3d first, then subtracting the full 3d plastic strain); (2) the wrapper methods for `RotatedMaterial`/`ReducedStressState` were ambiguous with the `NoMaterialState` fallback when wrapping a stateless material (fixed by adding explicit disambiguating overloads). The second round, run against the finished diff, reported no further findings. Test results: `Pkg.test()` passes (all pre-existing tests plus 16 new tests in test/test_current_stress.jl, including "frozen-state" checks that a different strain gives a purely elastic increment rather than triggering a fresh plastic/viscous correction). `docs/make.jl` builds cleanly (only pre-existing, unrelated warnings). This package has no literate tutorials/howtos, so that R2 step is not applicable here. Remaining/deferred: reduced-dimensional support for `GeneralizedMaxwell` and for `RotatedMaterial` wrapping a stateful material was not implemented in this prototype (calling `calculate_current_stress` on those combinations throws a clean MethodError rather than silently giving a wrong answer). FiniteStrainPlastic, CrystalPlasticity, and the hyperelastic models are also not covered, since the motivating FerriteAssembly.jl#94 use case and the issue's plasticity postprocessing scenario are small-strain. Co-Authored-By: Claude Sonnet 5 --- docs/src/small_strains.md | 5 ++ src/CurrentStress.jl | 81 ++++++++++++++++++++++++ src/MechanicalMaterialModels.jl | 3 + test/runtests.jl | 1 + test/test_current_stress.jl | 108 ++++++++++++++++++++++++++++++++ 5 files changed, 198 insertions(+) create mode 100644 src/CurrentStress.jl create mode 100644 test/test_current_stress.jl diff --git a/docs/src/small_strains.md b/docs/src/small_strains.md index 24c1ea3..33baa74 100644 --- a/docs/src/small_strains.md +++ b/docs/src/small_strains.md @@ -86,3 +86,8 @@ BCC12 GenericCrystallography CrystalPlasticity ``` + +## [Postprocessing](@id small_strain_postprocessing) +```@docs +calculate_current_stress +``` diff --git a/src/CurrentStress.jl b/src/CurrentStress.jl new file mode 100644 index 0000000..19c6373 --- /dev/null +++ b/src/CurrentStress.jl @@ -0,0 +1,81 @@ +""" + calculate_current_stress(m::AbstractMaterial, ϵ, state::AbstractMaterialState) + calculate_current_stress(rss::ReducedStressState, ϵ, state::AbstractMaterialState) + +Calculate the stress that is energy-conjugated to `ϵ`, consistent with the *given* +`state`, without invoking any local iteration that would advance history/internal +variables. `state` is normally the already-converged state obtained from a previous +call to `material_response` (e.g. during postprocessing, where `ϵ` may differ +slightly from the strain that produced `state`, such as an interpolated quadrature +point value). + +This is a prototype for [MaterialModelsBase.jl#12](https://github.com/KnutAM/MaterialModelsBase.jl/issues/12), +exploring how such an interface would work for different material models, +including support for a reduced-dimensional stress state +(see `MaterialModelsBase.ReducedStressState`) when only the +reduced-dimensional strain is supplied. This replaces the ad hoc, per-material +`calculate_stress` dispatch previously used for postprocessing in +[FerriteAssembly.jl#94](https://github.com/KnutAM/FerriteAssembly.jl/pull/94). + +Currently supported materials: [`LinearElastic`](@ref) (via the generic +`NoMaterialState` fallback), [`Plastic`](@ref), [`GeneralizedMaxwell`](@ref), +and [`RotatedMaterial`](@ref) wrapping any of these. Reduced-dimensional support +(via `ReducedStressState`) is currently implemented for `LinearElastic` and +`Plastic` only. +""" +function calculate_current_stress end + +# Generic fallback for genuinely stateless materials: since there is no history +# to accidentally advance, delegating to `material_response` is safe and exact. +function calculate_current_stress(m::AbstractMaterial, ϵ, state::MMB.NoMaterialState) + σ, _, _ = MMB.material_response(m, ϵ, state) + return σ +end + +function calculate_current_stress(stress_state::MMB.AbstractStressState, m::AbstractMaterial, ϵ, state::MMB.NoMaterialState) + σ, _, _, _ = MMB.material_response(stress_state, m, ϵ, state) + return σ +end + +# Plastic.jl +function calculate_current_stress(m::Plastic, ϵ::SymmetricTensor{2,3}, state::PlasticState) + return calculate_stress(m.elastic, ϵ - state.ϵp) +end + +function calculate_current_stress(stress_state::MMB.AbstractStressState, m::Plastic, ϵ, state::PlasticState) + # Expand the (possibly reduced) total strain to 3d before removing the plastic + # strain: for non-iterative states (e.g. PlaneStrain) the zero-padded + # out-of-plane *total* strain is exact by definition of the state, whereas + # reducing `state.ϵp` first would incorrectly discard its out-of-plane part. + ϵ_3d = MMB.expand_tensordim(stress_state, ϵ) + ϵₑ = ϵ_3d - state.ϵp + σ, _, _, _ = MMB.material_response(stress_state, m.elastic, ϵₑ, MMB.initial_material_state(m.elastic)) + return σ +end + +# ViscoElastic.jl +function calculate_current_stress(m::GeneralizedMaxwell, ϵ::SymmetricTensor{2,3}, state::GeneralizedMaxwellState) + σ0 = calculate_stress(m.base, ϵ) + return mapreduce((c, ϵv) -> 2 * c.G * (dev(ϵ) - ϵv), +, m.chains, state.ϵv; init=σ0) +end + +# RotatedMaterial.jl +function _calculate_current_stress_rotated(rm::RotatedMaterial, ϵ::SymmetricTensor{2,3}, state) + θ = norm(rm.rotation) + ϵ_rot = rotate(ϵ, rm.rotation, -θ) + σ_rot = calculate_current_stress(rm.material, ϵ_rot, state) + return rotate(σ_rot, rm.rotation, θ) +end +calculate_current_stress(rm::RotatedMaterial, ϵ::SymmetricTensor{2,3}, state) = _calculate_current_stress_rotated(rm, ϵ, state) +# Disambiguates against the `(AbstractMaterial, ϵ, ::NoMaterialState)` fallback above, +# which would otherwise be equally specific when `rm.material` is stateless. +calculate_current_stress(rm::RotatedMaterial, ϵ::SymmetricTensor{2,3}, state::MMB.NoMaterialState) = _calculate_current_stress_rotated(rm, ϵ, state) + +# ReducedStressState (MaterialModelsBase.jl) +function _calculate_current_stress_reduced(rss::MMB.ReducedStressState, ϵ, state) + return calculate_current_stress(rss.stress_state, rss.material, ϵ, state) +end +calculate_current_stress(rss::MMB.ReducedStressState, ϵ, state) = _calculate_current_stress_reduced(rss, ϵ, state) +# Disambiguates against the `(AbstractMaterial, ϵ, ::NoMaterialState)` fallback above, +# which would otherwise be equally specific when `rss.material` is stateless. +calculate_current_stress(rss::MMB.ReducedStressState, ϵ, state::MMB.NoMaterialState) = _calculate_current_stress_reduced(rss, ϵ, state) diff --git a/src/MechanicalMaterialModels.jl b/src/MechanicalMaterialModels.jl index 8a6ccff..d2155bb 100644 --- a/src/MechanicalMaterialModels.jl +++ b/src/MechanicalMaterialModels.jl @@ -57,4 +57,7 @@ export NeoHooke, CompressibleNeoHooke, SaintVenant include("FiniteStrainPlastic.jl") export FiniteStrainPlastic +include("CurrentStress.jl") +export calculate_current_stress + end diff --git a/test/runtests.jl b/test/runtests.jl index a0e18f8..772a575 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -17,6 +17,7 @@ include("test_viscoplastic.jl") include("test_viscoelastic.jl") include("test_differentiate.jl") include("test_crystal_plasticity.jl") +include("test_current_stress.jl") # Test finite strain behaviors include("test_hyperelastic.jl") diff --git a/test/test_current_stress.jl b/test/test_current_stress.jl new file mode 100644 index 0000000..e3f856c --- /dev/null +++ b/test/test_current_stress.jl @@ -0,0 +1,108 @@ +@testset "calculate_current_stress" begin + @testset "LinearElastic" begin + m = LinearElastic(E=210.e3, ν=0.3) + state = initial_material_state(m) + ϵ = rand(SymmetricTensor{2,3}) + @test calculate_current_stress(m, ϵ, state) ≈ m.C ⊡ ϵ + + # Reduced stress state: matches material_response's own (iterative) result, + # and the analytical plane-stress relation + E, ν = 210.e3, 0.3 + me = LinearElastic(; E, ν) + rss = ReducedStressState(PlaneStress(), me) + ϵ11 = 0.01 + ϵ_red = SymmetricTensor{2,2}((ϵ11, 0.0, 0.0)) + σ_direct = calculate_current_stress(rss, ϵ_red, initial_material_state(rss)) + σ_mr, _, _, _ = material_response(PlaneStress(), me, ϵ_red, initial_material_state(me)) + @test σ_direct ≈ σ_mr + @test σ_direct[1, 1] ≈ E / (1 - ν^2) * ϵ11 + @test σ_direct[2, 2] ≈ E / (1 - ν^2) * ν * ϵ11 + end + + @testset "Plastic" begin + E, ν = 210.e3, 0.3 + e = LinearElastic(; E, ν) + m = Plastic(elastic=e, yield=100.0, isotropic=Voce(Hiso=10.e3, κ∞=200.0), kinematic=ArmstrongFrederick(Hkin=1.e4, β∞=150.0)) + + # Load to a converged, plastically loaded state + state0 = initial_material_state(m) + ϵ1 = SymmetricTensor{2,3}((i, j) -> (i, j) == (1, 1) ? 0.01 : 0.0) + σ1, _, state1 = material_response(m, ϵ1, state0, nothing) + @test calculate_current_stress(m, ϵ1, state1) ≈ σ1 + + # Frozen-state postprocessing: a different strain should give a purely + # elastic increment from state1, NOT a fresh plastic correction. + ϵ2 = ϵ1 + SymmetricTensor{2,3}((i, j) -> (i, j) == (1, 1) ? 0.02 : 0.0) + σ2_frozen = calculate_current_stress(m, ϵ2, state1) + @test σ2_frozen ≈ σ1 + e.C ⊡ (ϵ2 - ϵ1) + σ2_true, _, state2_true = material_response(m, ϵ2, state1, nothing) + @test !(σ2_true ≈ σ2_frozen) # material_response would further evolve plastically + @test state2_true.ϵp != state1.ϵp + + # Reduced stress state (mirrors the FerriteAssembly#94 plane-stress postprocessing fix) + rss = ReducedStressState(PlaneStress(), m) + ϵ1_red = SymmetricTensor{2,2}((0.01, 0.0, 0.0)) + state0_red = initial_material_state(rss) + σ1_red, _, state1_red, _ = material_response(rss, ϵ1_red, state0_red, nothing) + σ1_red_current = calculate_current_stress(rss, ϵ1_red, state1_red) + @test σ1_red_current ≈ σ1_red + + # PlaneStrain: verify that the non-iterative shortcut retains the transverse + # plastic strain's elastic coupling (regression check for issue found in review) + rss_strain = ReducedStressState(PlaneStrain(), m) + state0_strain = initial_material_state(rss_strain) + σ1_strain, _, state1_strain, _ = material_response(rss_strain, ϵ1_red, state0_strain, nothing) + @test state1_strain.ϵp[3, 3] != 0 # sanity: this test only matters if ϵp33 != 0 + σ1_strain_current = calculate_current_stress(rss_strain, ϵ1_red, state1_strain) + @test σ1_strain_current ≈ σ1_strain + end + + @testset "GeneralizedMaxwell" begin + me = LinearElastic(E=210.e3, ν=0.3) + chain = Maxwell(G=1.e3, t=1.0) + m = GeneralizedMaxwell(me, chain) + + state0 = initial_material_state(m) + ϵ1 = rand(SymmetricTensor{2,3}) / 100 + σ1, _, state1 = material_response(m, ϵ1, state0, 0.5) + @test calculate_current_stress(m, ϵ1, state1) ≈ σ1 + + # Frozen-state: evaluating at a different strain must not re-solve the + # viscous strain evolution (which requires Δt); it should be a pure + # elastic-type increment using the given (fixed) viscous strain. + ϵ2 = ϵ1 + rand(SymmetricTensor{2,3}) / 100 + σ2_frozen = calculate_current_stress(m, ϵ2, state1) + σ2_expected = MechMat.calculate_stress(me, ϵ2) + 2 * chain.G * (dev(ϵ2) - state1.ϵv[1]) + @test σ2_frozen ≈ σ2_expected + σ2_true, _, _ = material_response(m, ϵ2, state1, 0.5) + @test !(σ2_true ≈ σ2_frozen) + end + + @testset "RotatedMaterial" begin + e = LinearElastic(E=210.e3, ν=0.3) + m_plastic = Plastic(elastic=e, yield=100.0, isotropic=Voce(Hiso=10.e3, κ∞=200.0), kinematic=ArmstrongFrederick(Hkin=1.e4, β∞=150.0)) + r = 2 * π * rand(Vec{3}) + rm = RotatedMaterial(m_plastic, r) + + state0 = initial_material_state(rm) + ϵ_global1 = SymmetricTensor{2,3}((i, j) -> (i, j) == (1, 1) ? 0.01 : 0.0) + _, _, state1 = material_response(rm, ϵ_global1, state0, nothing) + + ϵ_global2 = ϵ_global1 + SymmetricTensor{2,3}((i, j) -> (i, j) == (1, 1) ? 0.001 : 0.0) + σ_current = calculate_current_stress(rm, ϵ_global2, state1) + + θ = norm(r) + ϵ_local2 = rotate(ϵ_global2, r, -θ) + σ_local_expected = MechMat.calculate_stress(e, ϵ_local2 - state1.ϵp) + σ_expected = rotate(σ_local_expected, r, θ) + @test σ_current ≈ σ_expected + + # Stateless wrapped material: verify no dispatch ambiguity and correct rotation + m_el = LinearElastic{:cubicsymmetry}(C1111=1 + rand(), C1122=1 + rand(), C1212=1 + rand()) + rm_el = RotatedMaterial(m_el, r) + ϵ = rand(SymmetricTensor{2,3}) + σ_rm_el = calculate_current_stress(rm_el, ϵ, initial_material_state(rm_el)) + σ_local_el = calculate_current_stress(m_el, rotate(ϵ, r, -θ), initial_material_state(m_el)) + @test σ_rm_el ≈ rotate(σ_local_el, r, θ) + end +end From aae31dfa437b5602947e2a1af874c1406aa3f908 Mon Sep 17 00:00:00 2001 From: ClaudeBot Date: Fri, 18 Sep 2026 09:56:44 -0400 Subject: [PATCH 2/6] Extend calculate_current_stress to finite-strain models MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Detailed change: Extends the `calculate_current_stress` prototype (MaterialModelsBase.jl#12) from the previous commit to the finite-strain models in this package: - `NeoHooke`, `CompressibleNeoHooke`, `SaintVenant`: required no new code, for both full-3d and `ReducedStressState`-wrapped use (e.g. plane stress) — already covered by the existing generic `NoMaterialState` fallbacks, since MaterialModelsBase's stress-state iteration machinery already supports finite-strain (`Tensor{2,3}`-based) reduced states generically. Added tests to confirm this. - `FiniteStrainPlastic`: added a full-3d method that reuses the existing internal `calculate_PKstress(m, state, F)` helper, which already computes the frozen-`Fp` (converged state, no local Newton re-solve) 1st Piola-Kirchhoff stress. Added reduced-dimensional support via a small internal `FrozenStressMaterial <: AbstractMaterial` wrapper that packages the frozen PK-stress closure as a material, so it can ride MaterialModelsBase's existing stress-state Newton iteration (e.g. for `PlaneStress`); the tangent needed for that iteration is obtained via `Tensors.gradient` automatic differentiation. Not covered (documented as a limitation in the docstring): `CrystalPlasticity` (small-strain despite its docstring mentioning a finite-strain framework — a separate, pre-existing gap), `ReducedStressState` for `GeneralizedMaxwell`, and `RotatedMaterial` wrapping a finite-strain material (the latter already errors in `RotatedMaterial`'s own `material_response` due to a hard `::SymmetricTensor{2,3}` type assertion, independent of this change). Went through the dual-review workflow (same Codex thread as the prior change's plan review... actually a new thread this session). The plan review caught a real bug before implementation: the planned `MMB.NoMaterialState()` call has no zero-arg constructor (`NoMaterialState{T}` is parametric), which would throw a `MethodError` before any stress iteration starts. Fixed by removing the default argument from `FrozenStressMaterial`'s `material_response` (the state is now always passed explicitly) and constructing `MMB.NoMaterialState{eltype(F)}()` explicitly at the call site. The diff review after implementation reported no further findings. Test results: `Pkg.test()` passes (27 tests in the `calculate_current_stress` testset, up from 16; full suite green). `docs/make.jl` builds cleanly (only pre-existing, unrelated warnings). Co-Authored-By: Claude Sonnet 5 --- docs/src/finite_strains.md | 7 ++++++ docs/src/small_strains.md | 2 ++ src/CurrentStress.jl | 46 +++++++++++++++++++++++++++++++---- test/test_current_stress.jl | 48 +++++++++++++++++++++++++++++++++++++ 4 files changed, 98 insertions(+), 5 deletions(-) diff --git a/docs/src/finite_strains.md b/docs/src/finite_strains.md index 4454512..e576df6 100644 --- a/docs/src/finite_strains.md +++ b/docs/src/finite_strains.md @@ -101,3 +101,10 @@ Otherwise, the overstress function, ``\eta(\varPhi)``, determines the evolution ```math \dot{\lambda} = \eta(\varPhi, Y_0 + \kappa) ``` + +## [Postprocessing](@id finite_strain_postprocessing) +[`calculate_current_stress`](@ref) (documented under [Postprocessing](@ref small_strain_postprocessing)) +also supports the finite-strain models on this page: [`NeoHooke`](@ref), +[`CompressibleNeoHooke`](@ref), [`SaintVenant`](@ref), and +[`FiniteStrainPlastic`](@ref) (including reduced-dimensional stress states +for the latter). diff --git a/docs/src/small_strains.md b/docs/src/small_strains.md index 33baa74..17243bb 100644 --- a/docs/src/small_strains.md +++ b/docs/src/small_strains.md @@ -88,6 +88,8 @@ CrystalPlasticity ``` ## [Postprocessing](@id small_strain_postprocessing) +While `calculate_current_stress` itself is not specific to small strains, it +is documented here (see also [Finite Strains](@ref finite_strain_postprocessing)). ```@docs calculate_current_stress ``` diff --git a/src/CurrentStress.jl b/src/CurrentStress.jl index 19c6373..70b7e80 100644 --- a/src/CurrentStress.jl +++ b/src/CurrentStress.jl @@ -17,11 +17,21 @@ reduced-dimensional strain is supplied. This replaces the ad hoc, per-material `calculate_stress` dispatch previously used for postprocessing in [FerriteAssembly.jl#94](https://github.com/KnutAM/FerriteAssembly.jl/pull/94). -Currently supported materials: [`LinearElastic`](@ref) (via the generic -`NoMaterialState` fallback), [`Plastic`](@ref), [`GeneralizedMaxwell`](@ref), -and [`RotatedMaterial`](@ref) wrapping any of these. Reduced-dimensional support -(via `ReducedStressState`) is currently implemented for `LinearElastic` and -`Plastic` only. +Currently supported materials: [`LinearElastic`](@ref), [`NeoHooke`](@ref), +[`CompressibleNeoHooke`](@ref), and [`SaintVenant`](@ref) (all via the generic +`NoMaterialState` fallback), [`Plastic`](@ref), [`FiniteStrainPlastic`](@ref), +[`GeneralizedMaxwell`](@ref), and [`RotatedMaterial`](@ref) wrapping any of +the small-strain materials above. Reduced-dimensional support (via +`ReducedStressState`) is currently implemented for `LinearElastic`, +`NeoHooke`, `CompressibleNeoHooke`, `SaintVenant`, `Plastic`, and +`FiniteStrainPlastic`. + +!!! note "Not (yet) supported" + `CrystalPlasticity` (small-strain, despite referencing a finite-strain + framework in its docstring), `GeneralizedMaxwell`/`RotatedMaterial` under + `ReducedStressState`, and `RotatedMaterial` wrapping a finite-strain + material (the latter already errors in `RotatedMaterial`'s own + `material_response`, independently of this function). """ function calculate_current_stress end @@ -59,6 +69,32 @@ function calculate_current_stress(m::GeneralizedMaxwell, ϵ::SymmetricTensor{2,3 return mapreduce((c, ϵv) -> 2 * c.G * (dev(ϵ) - ϵv), +, m.chains, state.ϵv; init=σ0) end +# FiniteStrainPlastic.jl +# `calculate_PKstress(m, state, F)` already computes the frozen-state (converged +# `state.Fp`, no Newton re-solve) 1st Piola-Kirchhoff stress; it is used +# internally for the elastic-predictor branch of `material_response`. +calculate_current_stress(m::FiniteStrainPlastic, F::Tensor{2,3}, state::FiniteStrainPlasticState) = calculate_PKstress(m, state, F) + +# Wraps a frozen-state stress formula (strain -> stress, at fixed history +# variables) as an `AbstractMaterial`, so that it can ride MaterialModelsBase's +# existing stress-state Newton iteration (e.g. for `PlaneStress`). The tangent +# needed for that iteration is obtained via automatic differentiation, exactly +# analogous to how `compute_stress_and_tangent` differentiates through +# `compute_stress` in `HyperElastic.jl`. +struct FrozenStressMaterial{F} <: AbstractMaterial + f::F # F::Tensor{2,3} -> P::Tensor{2,3} +end +function MMB.material_response(fm::FrozenStressMaterial, F::Tensor{2,3}, old::MMB.AbstractMaterialState, args...) + dPdF, P = Tensors.gradient(fm.f, F, :all) + return P, dPdF, old +end + +function calculate_current_stress(stress_state::MMB.AbstractStressState, m::FiniteStrainPlastic, F, state::FiniteStrainPlasticState) + frozen = FrozenStressMaterial(F_ -> calculate_PKstress(m, state, F_)) + σ, _, _, _ = MMB.material_response(stress_state, frozen, F, MMB.NoMaterialState{eltype(F)}()) + return σ +end + # RotatedMaterial.jl function _calculate_current_stress_rotated(rm::RotatedMaterial, ϵ::SymmetricTensor{2,3}, state) θ = norm(rm.rotation) diff --git a/test/test_current_stress.jl b/test/test_current_stress.jl index e3f856c..dcf7c6c 100644 --- a/test/test_current_stress.jl +++ b/test/test_current_stress.jl @@ -105,4 +105,52 @@ σ_local_el = calculate_current_stress(m_el, rotate(ϵ, r, -θ), initial_material_state(m_el)) @test σ_rm_el ≈ rotate(σ_local_el, r, θ) end + + @testset "HyperElastic" begin + models = (NeoHooke(G=1 + rand()), CompressibleNeoHooke(G=1 + rand(), K=10 + rand()), SaintVenant(LinearElastic(E=210.e3, ν=0.3))) + for m in models + state = initial_material_state(m) + F = one(Tensor{2,3}) + rand(Tensor{2,3}) / 20 + P, _, _ = material_response(m, F, state) + @test calculate_current_stress(m, F, state) ≈ P + + # Reduced stress state: this already works via the generic + # `NoMaterialState` fallback, since MaterialModelsBase's stress-state + # iteration machinery already supports finite-strain (Tensor{2,3}) + # reduced states generically. + rss = ReducedStressState(PlaneStress(), m) + F_red = one(Tensor{2,2}) + rand(Tensor{2,2}) / 20 + state_red = initial_material_state(rss) + P_red, _, _, _ = material_response(rss, F_red, state_red) + @test calculate_current_stress(rss, F_red, state_red) ≈ P_red + end + end + + @testset "FiniteStrainPlastic" begin + E, ν, Y0 = 210.e3, 0.3, 100.0 + nh = CompressibleNeoHooke(G=convert_hooke_param(:G; E, ν), K=convert_hooke_param(:K; E, ν)) + m = FiniteStrainPlastic(elastic=nh, yield=Y0, isotropic=Voce(Hiso=10.e3, κ∞=200.0), kinematic=ArmstrongFrederick(Hkin=1.e4, β∞=150.0)) + + # Load to a converged, plastically loaded state + state0 = initial_material_state(m) + F1 = Tensor{2,3}((i, j) -> i == j ? (i == 1 ? 1.02 : 1.0) : 0.0) + P1, _, state1 = material_response(m, F1, state0, nothing) + @test calculate_current_stress(m, F1, state1) ≈ P1 + @test state1.Fp != state0.Fp # sanity: this test only matters if plastic loading occurred + + # Frozen-state postprocessing: a different F should give the frozen-Fp + # elastic response, NOT a fresh plastic correction. + F2 = Tensor{2,3}((i, j) -> i == j ? (i == 1 ? 1.03 : 1.0) : 0.0) + σ2_frozen = calculate_current_stress(m, F2, state1) + σ2_true, _, state2_true = material_response(m, F2, state1, nothing) + @test !(σ2_true ≈ σ2_frozen) # material_response would further evolve plastically + @test state2_true.Fp != state1.Fp + + # Reduced stress state, via the FrozenStressMaterial + MMB stress-state iteration + rss = ReducedStressState(PlaneStress(), m) + F1_red = Tensor{2,2}((1.02, 0.0, 0.0, 1.0)) + state0_red = initial_material_state(rss) + P1_red, _, state1_red, _ = material_response(rss, F1_red, state0_red, nothing) + @test calculate_current_stress(rss, F1_red, state1_red) ≈ P1_red + end end From 06f737df8a241c82f04c3a393e2baeff48348244 Mon Sep 17 00:00:00 2001 From: ClaudeBot Date: Mon, 21 Sep 2026 07:20:44 -0400 Subject: [PATCH 3/6] Generalize calculate_current_stress's ReducedStressState support MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Detailed change: Addresses KnutAM's review feedback on PR #13: a material-model developer should only need to implement the full-dimensional calculate_current_stress(m, ϵ, state), and ReducedStressState support should then "work by itself" - this is intended to be upstreamed to MaterialModelsBase.jl later. - Added a generic reduced-dimensional fallback, calculate_current_stress(stress_state::AbstractStressState, m::AbstractMaterial, strain, state::AbstractMaterialState), which wraps the model's own full-dimensional calculate_current_stress in the existing FrozenStressMaterial helper and rides MaterialModelsBase's stress-state Newton iteration (e.g. PlaneStress), autodiff-ing through it for the tangent. - Removed the now-redundant FiniteStrainPlastic-specific reduced method (no efficiency loss, since it already used FrozenStressMaterial + autodiff internally). - Generalized FrozenStressMaterial's material_response from F::Tensor{2,3} to strain::SecondOrderTensor{3}, so the same wrapper serves both finite- and small-strain materials. - As a consequence, GeneralizedMaxwell and RotatedMaterial wrapping a stateful small-strain material now get reduced-dimensional support "for free" (previously unsupported); added tests for both. - Kept the reduced-dimensional NoMaterialState and Plastic-specific fast paths, since both avoid autodiff via a cheaper analytic alternative (Julia's dispatch picks the more specific method automatically, no ambiguity). Removed the now-unnecessary RotatedMaterial/ReducedStressState NoMaterialState disambiguation overloads, since the broad full-dimensional fallback they guarded against no longer exists (see next point). - Removed the full-3d generic NoMaterialState fallback entirely, replacing it with two dedicated, gradient-free methods: LinearElastic (calculate_stress(m, ϵ), no autodiff at all) and AbstractHyperElastic (F ⋅ compute_stress(m, tdot(F)), one unavoidable gradient instead of the two that material_response/compute_stress_and_tangent would otherwise compute and discard). Went through the dual-review workflow (new Codex thread). The plan review caught a real issue with my initial reading of the repo owner's 3rd review comment ("Can be removed, since we want to special case to avoid calculating the gradient in this case", on the full-3d NoMaterialState fallback): I initially assumed this was a misplaced GitHub comment anchor, since that fallback computes no gradient directly. Codex correctly identified that it computes one indirectly, for hyperelastic materials, and fixed the plan accordingly (see above). The diff review after implementation reported no further findings. Test results: Pkg.test() passes (29 tests in the calculate_current_stress testset, up from 27; full suite green). docs/make.jl builds cleanly. Co-Authored-By: Claude Sonnet 5 --- src/CurrentStress.jl | 101 +++++++++++++++++++++--------------- test/test_current_stress.jl | 17 ++++++ 2 files changed, 75 insertions(+), 43 deletions(-) diff --git a/src/CurrentStress.jl b/src/CurrentStress.jl index 70b7e80..93ec1e1 100644 --- a/src/CurrentStress.jl +++ b/src/CurrentStress.jl @@ -17,35 +17,43 @@ reduced-dimensional strain is supplied. This replaces the ad hoc, per-material `calculate_stress` dispatch previously used for postprocessing in [FerriteAssembly.jl#94](https://github.com/KnutAM/FerriteAssembly.jl/pull/94). +A material-model developer only needs to implement the full-dimensional method, +`calculate_current_stress(m::MyMaterial, ϵ, state::MyMaterialState)`. Support for +a reduced-dimensional stress state (via `ReducedStressState`) then follows +automatically from a generic fallback, which rides `MaterialModelsBase`'s +existing stress-state Newton iteration (e.g. `PlaneStress`) using an internal +`FrozenStressMaterial` wrapper, with the tangent obtained by automatic +differentiation. A specific reduced-dimensional method only needs to be added +when a cheaper, non-autodiff alternative exists (as done here for +`Plastic` and for stateless/`NoMaterialState` materials). + Currently supported materials: [`LinearElastic`](@ref), [`NeoHooke`](@ref), -[`CompressibleNeoHooke`](@ref), and [`SaintVenant`](@ref) (all via the generic -`NoMaterialState` fallback), [`Plastic`](@ref), [`FiniteStrainPlastic`](@ref), -[`GeneralizedMaxwell`](@ref), and [`RotatedMaterial`](@ref) wrapping any of -the small-strain materials above. Reduced-dimensional support (via -`ReducedStressState`) is currently implemented for `LinearElastic`, -`NeoHooke`, `CompressibleNeoHooke`, `SaintVenant`, `Plastic`, and -`FiniteStrainPlastic`. +[`CompressibleNeoHooke`](@ref), and [`SaintVenant`](@ref) (each with a +dedicated, gradient-free implementation), [`Plastic`](@ref), +[`FiniteStrainPlastic`](@ref), [`GeneralizedMaxwell`](@ref), and +[`RotatedMaterial`](@ref) wrapping any of these. Reduced-dimensional support +(via `ReducedStressState`) works for all of the above, generically for +`GeneralizedMaxwell` and for `RotatedMaterial` wrapping a small-strain +material (via the generic fallback), and with a dedicated +non-autodiff implementation for stateless materials and for `Plastic`. !!! note "Not (yet) supported" `CrystalPlasticity` (small-strain, despite referencing a finite-strain - framework in its docstring), `GeneralizedMaxwell`/`RotatedMaterial` under - `ReducedStressState`, and `RotatedMaterial` wrapping a finite-strain - material (the latter already errors in `RotatedMaterial`'s own - `material_response`, independently of this function). + framework in its docstring) has no `calculate_current_stress` method at + all yet. `RotatedMaterial` wrapping a finite-strain material already + errors in `RotatedMaterial`'s own `material_response` (a hard + `::SymmetricTensor{2,3}` type assertion), independently of this function. """ function calculate_current_stress end -# Generic fallback for genuinely stateless materials: since there is no history -# to accidentally advance, delegating to `material_response` is safe and exact. -function calculate_current_stress(m::AbstractMaterial, ϵ, state::MMB.NoMaterialState) - σ, _, _ = MMB.material_response(m, ϵ, state) - return σ -end +# LinearElastic.jl: stress-only, no gradient (material_response would compute one, +# via `m.C`, that `calculate_current_stress` doesn't need). +calculate_current_stress(m::LinearElastic, ϵ::SymmetricTensor{2,3}, ::MMB.NoMaterialState) = calculate_stress(m, ϵ) -function calculate_current_stress(stress_state::MMB.AbstractStressState, m::AbstractMaterial, ϵ, state::MMB.NoMaterialState) - σ, _, _, _ = MMB.material_response(stress_state, m, ϵ, state) - return σ -end +# HyperElastic.jl: stress-only. Computing `S = 2 ∂Ψ/∂C` (once) is unavoidable to get +# the stress at all, but `material_response` additionally differentiates through +# that once more to get the tangent, which `calculate_current_stress` doesn't need. +calculate_current_stress(m::AbstractHyperElastic, F::Tensor{2,3}, ::MMB.NoMaterialState) = F ⋅ compute_stress(m, tdot(F)) # Plastic.jl function calculate_current_stress(m::Plastic, ϵ::SymmetricTensor{2,3}, state::PlasticState) @@ -57,6 +65,8 @@ function calculate_current_stress(stress_state::MMB.AbstractStressState, m::Plas # strain: for non-iterative states (e.g. PlaneStrain) the zero-padded # out-of-plane *total* strain is exact by definition of the state, whereas # reducing `state.ϵp` first would incorrectly discard its out-of-plane part. + # This avoids autodiff entirely, by delegating to `m.elastic`'s own analytic + # stress-state response. ϵ_3d = MMB.expand_tensordim(stress_state, ϵ) ϵₑ = ϵ_3d - state.ϵp σ, _, _, _ = MMB.material_response(stress_state, m.elastic, ϵₑ, MMB.initial_material_state(m.elastic)) @@ -75,43 +85,48 @@ end # internally for the elastic-predictor branch of `material_response`. calculate_current_stress(m::FiniteStrainPlastic, F::Tensor{2,3}, state::FiniteStrainPlasticState) = calculate_PKstress(m, state, F) -# Wraps a frozen-state stress formula (strain -> stress, at fixed history +# Wraps a frozen-state stress formula, `f`, mapping a strain (`SecondOrderTensor{3}`, +# i.e. `Tensor{2,3}` or `SymmetricTensor{2,3}`) to a stress (at fixed history/internal # variables) as an `AbstractMaterial`, so that it can ride MaterialModelsBase's -# existing stress-state Newton iteration (e.g. for `PlaneStress`). The tangent -# needed for that iteration is obtained via automatic differentiation, exactly -# analogous to how `compute_stress_and_tangent` differentiates through -# `compute_stress` in `HyperElastic.jl`. +# existing stress-state Newton iteration (e.g. for `PlaneStress`). The tangent needed +# for that iteration is obtained via automatic differentiation. This is what powers +# the generic reduced-dimensional fallback of `calculate_current_stress` below. struct FrozenStressMaterial{F} <: AbstractMaterial - f::F # F::Tensor{2,3} -> P::Tensor{2,3} + f::F +end +function MMB.material_response(fm::FrozenStressMaterial, strain::SecondOrderTensor{3}, old::MMB.AbstractMaterialState, args::Vararg{Any,N}) where {N} + dσdϵ, σ = Tensors.gradient(fm.f, strain, :all) + return σ, dσdϵ, old end -function MMB.material_response(fm::FrozenStressMaterial, F::Tensor{2,3}, old::MMB.AbstractMaterialState, args...) - dPdF, P = Tensors.gradient(fm.f, F, :all) - return P, dPdF, old + +# Generic reduced-dimensional fallback: as long as `calculate_current_stress(m, ϵ, +# state)` (full-dimensional) is implemented for `m`, this makes `ReducedStressState` +# support "just work", by autodiff-ing through it. More specific methods above/below +# (e.g. for `Plastic` or `NoMaterialState`) take precedence when a cheaper, +# non-autodiff alternative exists. +function calculate_current_stress(stress_state::MMB.AbstractStressState, m::AbstractMaterial, strain, state::MMB.AbstractMaterialState) + frozen = FrozenStressMaterial(e -> calculate_current_stress(m, e, state)) + σ, _, _, _ = MMB.material_response(stress_state, frozen, strain, MMB.NoMaterialState{eltype(strain)}()) + return σ end -function calculate_current_stress(stress_state::MMB.AbstractStressState, m::FiniteStrainPlastic, F, state::FiniteStrainPlasticState) - frozen = FrozenStressMaterial(F_ -> calculate_PKstress(m, state, F_)) - σ, _, _, _ = MMB.material_response(stress_state, frozen, F, MMB.NoMaterialState{eltype(F)}()) +# Reduced-dimensional fast path for stateless materials: avoids the autodiff in the +# generic fallback above by delegating directly to `material_response`'s own +# (analytic, for `LinearElastic`) stress-state handling. +function calculate_current_stress(stress_state::MMB.AbstractStressState, m::AbstractMaterial, strain, state::MMB.NoMaterialState) + σ, _, _, _ = MMB.material_response(stress_state, m, strain, state) return σ end # RotatedMaterial.jl -function _calculate_current_stress_rotated(rm::RotatedMaterial, ϵ::SymmetricTensor{2,3}, state) +function calculate_current_stress(rm::RotatedMaterial, ϵ::SymmetricTensor{2,3}, state) θ = norm(rm.rotation) ϵ_rot = rotate(ϵ, rm.rotation, -θ) σ_rot = calculate_current_stress(rm.material, ϵ_rot, state) return rotate(σ_rot, rm.rotation, θ) end -calculate_current_stress(rm::RotatedMaterial, ϵ::SymmetricTensor{2,3}, state) = _calculate_current_stress_rotated(rm, ϵ, state) -# Disambiguates against the `(AbstractMaterial, ϵ, ::NoMaterialState)` fallback above, -# which would otherwise be equally specific when `rm.material` is stateless. -calculate_current_stress(rm::RotatedMaterial, ϵ::SymmetricTensor{2,3}, state::MMB.NoMaterialState) = _calculate_current_stress_rotated(rm, ϵ, state) # ReducedStressState (MaterialModelsBase.jl) -function _calculate_current_stress_reduced(rss::MMB.ReducedStressState, ϵ, state) +function calculate_current_stress(rss::MMB.ReducedStressState, ϵ, state) return calculate_current_stress(rss.stress_state, rss.material, ϵ, state) end -calculate_current_stress(rss::MMB.ReducedStressState, ϵ, state) = _calculate_current_stress_reduced(rss, ϵ, state) -# Disambiguates against the `(AbstractMaterial, ϵ, ::NoMaterialState)` fallback above, -# which would otherwise be equally specific when `rss.material` is stateless. -calculate_current_stress(rss::MMB.ReducedStressState, ϵ, state::MMB.NoMaterialState) = _calculate_current_stress_reduced(rss, ϵ, state) diff --git a/test/test_current_stress.jl b/test/test_current_stress.jl index dcf7c6c..8998067 100644 --- a/test/test_current_stress.jl +++ b/test/test_current_stress.jl @@ -76,6 +76,15 @@ @test σ2_frozen ≈ σ2_expected σ2_true, _, _ = material_response(m, ϵ2, state1, 0.5) @test !(σ2_true ≈ σ2_frozen) + + # Reduced stress state: previously unsupported, now works automatically + # via the generic fallback (GeneralizedMaxwell has no dedicated reduced + # method, only the full-dimensional one used above). + rss = ReducedStressState(PlaneStress(), m) + ϵ1_red = SymmetricTensor{2,2}((0.01, 0.0, 0.0)) + state0_red = initial_material_state(rss) + σ1_red, _, state1_red, _ = material_response(rss, ϵ1_red, state0_red, 0.5) + @test calculate_current_stress(rss, ϵ1_red, state1_red) ≈ σ1_red end @testset "RotatedMaterial" begin @@ -104,6 +113,14 @@ σ_rm_el = calculate_current_stress(rm_el, ϵ, initial_material_state(rm_el)) σ_local_el = calculate_current_stress(m_el, rotate(ϵ, r, -θ), initial_material_state(m_el)) @test σ_rm_el ≈ rotate(σ_local_el, r, θ) + + # Reduced stress state wrapping a rotated, stateful material: previously + # unsupported, now works automatically via the generic fallback. + rss = ReducedStressState(PlaneStress(), rm) + ϵ1_red = SymmetricTensor{2,2}((0.01, 0.0, 0.0)) + state0_red = initial_material_state(rss) + σ1_red, _, state1_red, _ = material_response(rss, ϵ1_red, state0_red, nothing) + @test calculate_current_stress(rss, ϵ1_red, state1_red) ≈ σ1_red end @testset "HyperElastic" begin From be9fa3357aea08daa8458ac79ca1befcab89f987 Mon Sep 17 00:00:00 2001 From: ClaudeBot Date: Mon, 21 Sep 2026 10:46:22 -0400 Subject: [PATCH 4/6] Migrate to upstreamed MaterialModelsBase.stress_from_state Detailed change: The postprocessing interface prototyped in this repo (calculate_current_stress) has been upstreamed into MaterialModelsBase.jl itself as stress_from_state (KnutAM/MaterialModelsBase.jl#21, currently open, on branch knutambot/MaterialModelsBase.jl#cb/calculate_current_stress). This keeps only the material-specific implementations here, extending MMB.stress_from_state: - src/CurrentStress.jl renamed to src/StressFromState.jl, trimmed to only: LinearElastic and AbstractHyperElastic (gradient-free NoMaterialState overrides - upstream's own generic NoMaterialState fallback deliberately discards a tangent for stateless materials, matching this repo's earlier design choice; a material-specific override remains the recommended, opt-in way to avoid that), Plastic (full-dimensional formula plus a reduced-dimensional fast path, both cheaper than upstream's autodiff-based generic fallback), GeneralizedMaxwell and FiniteStrainPlastic (full-dimensional only now - their reduced-dimensional support comes entirely from upstream's generic fallback, which uses its own FrozenStressMaterial), and RotatedMaterial (full-dimensional wrapper, keeping the generic + NoMaterialState-disambiguation two-method pattern, since upstream's own NoMaterialState fallback would otherwise be ambiguous against a single generic method here). - Deleted: the local FrozenStressMaterial struct, the generic reduced-dimensional fallback, and ReducedStressState delegation - all now live in MaterialModelsBase.jl. - Project.toml and docs/Project.toml: temporarily pin [sources] for MaterialModelsBase to the cb/calculate_current_stress branch, with an inline comment marking this temporary and noting to revert (and tighten [compat]) once MaterialModelsBase.jl#21 merges and releases. - Renamed calculate_current_stress -> stress_from_state throughout (source, tests, docs); dropped the export line, since stress_from_state is not re-exported by this package, matching how material_response itself isn't re-exported either. - docs/src/small_strains.md and docs/src/finite_strains.md: replaced the local @docs stress_from_state blocks with prose, since the canonical docstring now lives upstream and this package's makedocs(modules=[MechanicalMaterialModels]) call wouldn't surface it anyway. Went through the dual-review workflow (new Codex thread). The plan review caught 4 real issues before implementation: (1) RotatedMaterial becomes ambiguous against upstream's own NoMaterialState fallback for a stateless wrapped material - fixed by restoring the two-method disambiguation pattern; (2) docs/Project.toml needs its own [sources] override too, since dependency [sources] entries aren't inherited into sibling environments - fixed; (3) removing the local docstring would leave @docs stress_from_state unable to find any docstring, since makedocs filters docstrings to the listed modules - fixed by using prose instead of @docs; (4) the unchanged [compat] range ("0.3, 0.4") doesn't protect downstream users once this merges for real - documented as a follow-up requirement, since there's no real release to reference yet. The diff review after implementation reported no further findings (it noted one pre-existing, unrelated upstream dispatch ambiguity in get_drdx for finite-strain GeneralStressState, not touched by or affecting this change). Test results: Pkg.test() passes against the upstream branch (29 tests in the stress_from_state testset, full suite green). docs/make.jl builds cleanly (only pre-existing, unrelated warnings). Remaining risk: this branch temporarily depends on an unmerged, force-pushable fork branch. If that branch is rewritten or deleted, this repo's CI/tests will break until the [sources] override is reverted or updated - explicitly acknowledged as intentional/temporary, per instruction, not a bug. This must be reverted before merging MaterialModelsBase.jl#21. Co-Authored-By: Claude Sonnet 5 --- Project.toml | 6 +- docs/Project.toml | 3 +- docs/src/finite_strains.md | 9 +- docs/src/small_strains.md | 22 ++- src/CurrentStress.jl | 132 ------------------ src/MechanicalMaterialModels.jl | 3 +- src/StressFromState.jl | 72 ++++++++++ test/runtests.jl | 2 +- ...nt_stress.jl => test_stress_from_state.jl} | 38 ++--- 9 files changed, 121 insertions(+), 166 deletions(-) delete mode 100644 src/CurrentStress.jl create mode 100644 src/StressFromState.jl rename test/{test_current_stress.jl => test_stress_from_state.jl} (85%) diff --git a/Project.toml b/Project.toml index 90b42f3..f2de29f 100644 --- a/Project.toml +++ b/Project.toml @@ -27,6 +27,10 @@ MaterialModelsTesting = "882b014b-b96c-4115-8629-e17fb35110d2" test = ["Test", "ForwardDiff", "FiniteDiff", "MaterialModelsTesting"] [sources] -MaterialModelsBase = {url = "https://github.com/knutam/MaterialModelsBase.jl"} +# TEMPORARY: pins to the branch implementing `stress_from_state` +# (https://github.com/KnutAM/MaterialModelsBase.jl/pull/21). Revert to +# `{url = "https://github.com/knutam/MaterialModelsBase.jl"}` once that PR +# merges and is released (and tighten the `[compat]` bound above accordingly). +MaterialModelsBase = {url = "https://github.com/knutambot/MaterialModelsBase.jl", rev = "cb/calculate_current_stress"} Newton = {url = "https://github.com/knutam/Newton.jl"} MaterialModelsTesting = {url = "https://github.com/KnutAM/MaterialModelsTesting.jl"} diff --git a/docs/Project.toml b/docs/Project.toml index 6c5ef2e..fe959e2 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -8,5 +8,6 @@ Newton = "83aa5b51-0588-403c-85e4-434ec185aae7" Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4" [sources] -MaterialModelsBase = {url = "https://github.com/knutam/MaterialModelsBase.jl"} +# TEMPORARY: see the matching note in ../Project.toml - revert together. +MaterialModelsBase = {url = "https://github.com/knutambot/MaterialModelsBase.jl", rev = "cb/calculate_current_stress"} Newton = {url = "https://github.com/knutam/Newton.jl"} diff --git a/docs/src/finite_strains.md b/docs/src/finite_strains.md index e576df6..8809f2c 100644 --- a/docs/src/finite_strains.md +++ b/docs/src/finite_strains.md @@ -103,8 +103,7 @@ Otherwise, the overstress function, ``\eta(\varPhi)``, determines the evolution ``` ## [Postprocessing](@id finite_strain_postprocessing) -[`calculate_current_stress`](@ref) (documented under [Postprocessing](@ref small_strain_postprocessing)) -also supports the finite-strain models on this page: [`NeoHooke`](@ref), -[`CompressibleNeoHooke`](@ref), [`SaintVenant`](@ref), and -[`FiniteStrainPlastic`](@ref) (including reduced-dimensional stress states -for the latter). +`MaterialModelsBase.stress_from_state` (see [Postprocessing](@ref small_strain_postprocessing)) +also has methods here for [`NeoHooke`](@ref), [`CompressibleNeoHooke`](@ref), +[`SaintVenant`](@ref), and [`FiniteStrainPlastic`](@ref) (including +lower-dimensional stress states for the latter). diff --git a/docs/src/small_strains.md b/docs/src/small_strains.md index 17243bb..501d953 100644 --- a/docs/src/small_strains.md +++ b/docs/src/small_strains.md @@ -88,8 +88,20 @@ CrystalPlasticity ``` ## [Postprocessing](@id small_strain_postprocessing) -While `calculate_current_stress` itself is not specific to small strains, it -is documented here (see also [Finite Strains](@ref finite_strain_postprocessing)). -```@docs -calculate_current_stress -``` +`MaterialModelsBase.stress_from_state` ([documented there](https://github.com/KnutAM/MaterialModelsBase.jl/pull/21), +pending release) provides a generic postprocessing interface: given a strain +and an already-converged state, it returns the corresponding stress without +advancing any history/internal variables. This package implements it (i.e. +extends `stress_from_state` with a method) for [`LinearElastic`](@ref), +[`Plastic`](@ref), and [`GeneralizedMaxwell`](@ref), including +lower-dimensional stress states (via `MaterialModelsBase.ReducedStressState`) +for `LinearElastic` and `Plastic` (see also +[Finite Strains](@ref finite_strain_postprocessing)). + +!!! note "Not (yet) supported" + `CrystalPlasticity` (small-strain, despite referencing a finite-strain + framework in its docstring) has no `stress_from_state` method. + `RotatedMaterial` wrapping a finite-strain material errors in + `RotatedMaterial`'s own `material_response` (a hard + `::SymmetricTensor{2,3}` type assertion), independently of this + interface. diff --git a/src/CurrentStress.jl b/src/CurrentStress.jl deleted file mode 100644 index 93ec1e1..0000000 --- a/src/CurrentStress.jl +++ /dev/null @@ -1,132 +0,0 @@ -""" - calculate_current_stress(m::AbstractMaterial, ϵ, state::AbstractMaterialState) - calculate_current_stress(rss::ReducedStressState, ϵ, state::AbstractMaterialState) - -Calculate the stress that is energy-conjugated to `ϵ`, consistent with the *given* -`state`, without invoking any local iteration that would advance history/internal -variables. `state` is normally the already-converged state obtained from a previous -call to `material_response` (e.g. during postprocessing, where `ϵ` may differ -slightly from the strain that produced `state`, such as an interpolated quadrature -point value). - -This is a prototype for [MaterialModelsBase.jl#12](https://github.com/KnutAM/MaterialModelsBase.jl/issues/12), -exploring how such an interface would work for different material models, -including support for a reduced-dimensional stress state -(see `MaterialModelsBase.ReducedStressState`) when only the -reduced-dimensional strain is supplied. This replaces the ad hoc, per-material -`calculate_stress` dispatch previously used for postprocessing in -[FerriteAssembly.jl#94](https://github.com/KnutAM/FerriteAssembly.jl/pull/94). - -A material-model developer only needs to implement the full-dimensional method, -`calculate_current_stress(m::MyMaterial, ϵ, state::MyMaterialState)`. Support for -a reduced-dimensional stress state (via `ReducedStressState`) then follows -automatically from a generic fallback, which rides `MaterialModelsBase`'s -existing stress-state Newton iteration (e.g. `PlaneStress`) using an internal -`FrozenStressMaterial` wrapper, with the tangent obtained by automatic -differentiation. A specific reduced-dimensional method only needs to be added -when a cheaper, non-autodiff alternative exists (as done here for -`Plastic` and for stateless/`NoMaterialState` materials). - -Currently supported materials: [`LinearElastic`](@ref), [`NeoHooke`](@ref), -[`CompressibleNeoHooke`](@ref), and [`SaintVenant`](@ref) (each with a -dedicated, gradient-free implementation), [`Plastic`](@ref), -[`FiniteStrainPlastic`](@ref), [`GeneralizedMaxwell`](@ref), and -[`RotatedMaterial`](@ref) wrapping any of these. Reduced-dimensional support -(via `ReducedStressState`) works for all of the above, generically for -`GeneralizedMaxwell` and for `RotatedMaterial` wrapping a small-strain -material (via the generic fallback), and with a dedicated -non-autodiff implementation for stateless materials and for `Plastic`. - -!!! note "Not (yet) supported" - `CrystalPlasticity` (small-strain, despite referencing a finite-strain - framework in its docstring) has no `calculate_current_stress` method at - all yet. `RotatedMaterial` wrapping a finite-strain material already - errors in `RotatedMaterial`'s own `material_response` (a hard - `::SymmetricTensor{2,3}` type assertion), independently of this function. -""" -function calculate_current_stress end - -# LinearElastic.jl: stress-only, no gradient (material_response would compute one, -# via `m.C`, that `calculate_current_stress` doesn't need). -calculate_current_stress(m::LinearElastic, ϵ::SymmetricTensor{2,3}, ::MMB.NoMaterialState) = calculate_stress(m, ϵ) - -# HyperElastic.jl: stress-only. Computing `S = 2 ∂Ψ/∂C` (once) is unavoidable to get -# the stress at all, but `material_response` additionally differentiates through -# that once more to get the tangent, which `calculate_current_stress` doesn't need. -calculate_current_stress(m::AbstractHyperElastic, F::Tensor{2,3}, ::MMB.NoMaterialState) = F ⋅ compute_stress(m, tdot(F)) - -# Plastic.jl -function calculate_current_stress(m::Plastic, ϵ::SymmetricTensor{2,3}, state::PlasticState) - return calculate_stress(m.elastic, ϵ - state.ϵp) -end - -function calculate_current_stress(stress_state::MMB.AbstractStressState, m::Plastic, ϵ, state::PlasticState) - # Expand the (possibly reduced) total strain to 3d before removing the plastic - # strain: for non-iterative states (e.g. PlaneStrain) the zero-padded - # out-of-plane *total* strain is exact by definition of the state, whereas - # reducing `state.ϵp` first would incorrectly discard its out-of-plane part. - # This avoids autodiff entirely, by delegating to `m.elastic`'s own analytic - # stress-state response. - ϵ_3d = MMB.expand_tensordim(stress_state, ϵ) - ϵₑ = ϵ_3d - state.ϵp - σ, _, _, _ = MMB.material_response(stress_state, m.elastic, ϵₑ, MMB.initial_material_state(m.elastic)) - return σ -end - -# ViscoElastic.jl -function calculate_current_stress(m::GeneralizedMaxwell, ϵ::SymmetricTensor{2,3}, state::GeneralizedMaxwellState) - σ0 = calculate_stress(m.base, ϵ) - return mapreduce((c, ϵv) -> 2 * c.G * (dev(ϵ) - ϵv), +, m.chains, state.ϵv; init=σ0) -end - -# FiniteStrainPlastic.jl -# `calculate_PKstress(m, state, F)` already computes the frozen-state (converged -# `state.Fp`, no Newton re-solve) 1st Piola-Kirchhoff stress; it is used -# internally for the elastic-predictor branch of `material_response`. -calculate_current_stress(m::FiniteStrainPlastic, F::Tensor{2,3}, state::FiniteStrainPlasticState) = calculate_PKstress(m, state, F) - -# Wraps a frozen-state stress formula, `f`, mapping a strain (`SecondOrderTensor{3}`, -# i.e. `Tensor{2,3}` or `SymmetricTensor{2,3}`) to a stress (at fixed history/internal -# variables) as an `AbstractMaterial`, so that it can ride MaterialModelsBase's -# existing stress-state Newton iteration (e.g. for `PlaneStress`). The tangent needed -# for that iteration is obtained via automatic differentiation. This is what powers -# the generic reduced-dimensional fallback of `calculate_current_stress` below. -struct FrozenStressMaterial{F} <: AbstractMaterial - f::F -end -function MMB.material_response(fm::FrozenStressMaterial, strain::SecondOrderTensor{3}, old::MMB.AbstractMaterialState, args::Vararg{Any,N}) where {N} - dσdϵ, σ = Tensors.gradient(fm.f, strain, :all) - return σ, dσdϵ, old -end - -# Generic reduced-dimensional fallback: as long as `calculate_current_stress(m, ϵ, -# state)` (full-dimensional) is implemented for `m`, this makes `ReducedStressState` -# support "just work", by autodiff-ing through it. More specific methods above/below -# (e.g. for `Plastic` or `NoMaterialState`) take precedence when a cheaper, -# non-autodiff alternative exists. -function calculate_current_stress(stress_state::MMB.AbstractStressState, m::AbstractMaterial, strain, state::MMB.AbstractMaterialState) - frozen = FrozenStressMaterial(e -> calculate_current_stress(m, e, state)) - σ, _, _, _ = MMB.material_response(stress_state, frozen, strain, MMB.NoMaterialState{eltype(strain)}()) - return σ -end - -# Reduced-dimensional fast path for stateless materials: avoids the autodiff in the -# generic fallback above by delegating directly to `material_response`'s own -# (analytic, for `LinearElastic`) stress-state handling. -function calculate_current_stress(stress_state::MMB.AbstractStressState, m::AbstractMaterial, strain, state::MMB.NoMaterialState) - σ, _, _, _ = MMB.material_response(stress_state, m, strain, state) - return σ -end - -# RotatedMaterial.jl -function calculate_current_stress(rm::RotatedMaterial, ϵ::SymmetricTensor{2,3}, state) - θ = norm(rm.rotation) - ϵ_rot = rotate(ϵ, rm.rotation, -θ) - σ_rot = calculate_current_stress(rm.material, ϵ_rot, state) - return rotate(σ_rot, rm.rotation, θ) -end - -# ReducedStressState (MaterialModelsBase.jl) -function calculate_current_stress(rss::MMB.ReducedStressState, ϵ, state) - return calculate_current_stress(rss.stress_state, rss.material, ϵ, state) -end diff --git a/src/MechanicalMaterialModels.jl b/src/MechanicalMaterialModels.jl index d2155bb..8b01ddc 100644 --- a/src/MechanicalMaterialModels.jl +++ b/src/MechanicalMaterialModels.jl @@ -57,7 +57,6 @@ export NeoHooke, CompressibleNeoHooke, SaintVenant include("FiniteStrainPlastic.jl") export FiniteStrainPlastic -include("CurrentStress.jl") -export calculate_current_stress +include("StressFromState.jl") end diff --git a/src/StressFromState.jl b/src/StressFromState.jl new file mode 100644 index 0000000..84ba99a --- /dev/null +++ b/src/StressFromState.jl @@ -0,0 +1,72 @@ +# Material-specific implementations of `MaterialModelsBase.stress_from_state` +# (https://github.com/KnutAM/MaterialModelsBase.jl/pull/21), which upstreamed the +# generic postprocessing machinery originally prototyped in this file (see +# https://github.com/KnutAM/MechanicalMaterialModels.jl/pull/13): the generic +# `NoMaterialState` fallbacks, `FrozenStressMaterial`, the generic +# reduced-dimensional fallback, and `ReducedStressState` delegation now all live +# in `MaterialModelsBase.jl` itself. Only the material-specific methods below - +# either required (no `NoMaterialState`-based default exists) or a cheaper, +# non-autodiff alternative to the generic fallback - remain here. +# +# NOTE: while MaterialModelsBase.jl#21 is not yet merged, Project.toml and +# docs/Project.toml temporarily point MaterialModelsBase at the branch +# implementing it (knutambot/MaterialModelsBase.jl#cb/calculate_current_stress). +# Revert both `[sources]` entries (and tighten the `MaterialModelsBase` `[compat]` +# bound to the first release containing `stress_from_state`) once that PR merges +# and is released. + +# LinearElastic.jl: stress-only, no gradient (material_response would compute one, +# via `m.C`, that `stress_from_state`'s generic `NoMaterialState` fallback would +# otherwise compute via `material_response` and discard). +MMB.stress_from_state(m::LinearElastic, ϵ::SymmetricTensor{2,3}, ::MMB.NoMaterialState) = calculate_stress(m, ϵ) + +# HyperElastic.jl: stress-only. Computing `S = 2 ∂Ψ/∂C` (once) is unavoidable to get +# the stress at all, but the generic `NoMaterialState` fallback's `material_response` +# call additionally differentiates through that once more to get the tangent, which +# `stress_from_state` doesn't need. +MMB.stress_from_state(m::AbstractHyperElastic, F::Tensor{2,3}, ::MMB.NoMaterialState) = F ⋅ compute_stress(m, tdot(F)) + +# Plastic.jl +function MMB.stress_from_state(m::Plastic, ϵ::SymmetricTensor{2,3}, state::PlasticState) + return calculate_stress(m.elastic, ϵ - state.ϵp) +end + +function MMB.stress_from_state(stress_state::MMB.AbstractStressState, m::Plastic, ϵ, state::PlasticState) + # Expand the (possibly reduced) total strain to 3d before removing the plastic + # strain: for non-iterative states (e.g. PlaneStrain) the zero-padded + # out-of-plane *total* strain is exact by definition of the state, whereas + # reducing `state.ϵp` first would incorrectly discard its out-of-plane part. + # This avoids autodiff entirely, by delegating to `m.elastic`'s own analytic + # stress-state response. + ϵ_3d = MMB.expand_tensordim(stress_state, ϵ) + ϵₑ = ϵ_3d - state.ϵp + σ, _, _, _ = MMB.material_response(stress_state, m.elastic, ϵₑ, MMB.initial_material_state(m.elastic)) + return σ +end + +# ViscoElastic.jl +function MMB.stress_from_state(m::GeneralizedMaxwell, ϵ::SymmetricTensor{2,3}, state::GeneralizedMaxwellState) + σ0 = calculate_stress(m.base, ϵ) + return mapreduce((c, ϵv) -> 2 * c.G * (dev(ϵ) - ϵv), +, m.chains, state.ϵv; init=σ0) +end + +# FiniteStrainPlastic.jl +# `calculate_PKstress(m, state, F)` already computes the frozen-state (converged +# `state.Fp`, no Newton re-solve) 1st Piola-Kirchhoff stress; it is used +# internally for the elastic-predictor branch of `material_response`. No +# reduced-dimensional method is needed: MaterialModelsBase's generic fallback +# (autodiff-ing through this method via its own `FrozenStressMaterial`) covers it. +MMB.stress_from_state(m::FiniteStrainPlastic, F::Tensor{2,3}, state::FiniteStrainPlasticState) = calculate_PKstress(m, state, F) + +# RotatedMaterial.jl +function _stress_from_state_rotated(rm::RotatedMaterial, ϵ::SymmetricTensor{2,3}, state) + θ = norm(rm.rotation) + ϵ_rot = rotate(ϵ, rm.rotation, -θ) + σ_rot = MMB.stress_from_state(rm.material, ϵ_rot, state) + return rotate(σ_rot, rm.rotation, θ) +end +MMB.stress_from_state(rm::RotatedMaterial, ϵ::SymmetricTensor{2,3}, state) = _stress_from_state_rotated(rm, ϵ, state) +# Disambiguates against MaterialModelsBase's own `(AbstractMaterial, ϵ, +# ::NoMaterialState)` fallback, which would otherwise be equally specific when +# `rm.material` is stateless. +MMB.stress_from_state(rm::RotatedMaterial, ϵ::SymmetricTensor{2,3}, state::MMB.NoMaterialState) = _stress_from_state_rotated(rm, ϵ, state) diff --git a/test/runtests.jl b/test/runtests.jl index 772a575..5c0ccea 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -17,7 +17,7 @@ include("test_viscoplastic.jl") include("test_viscoelastic.jl") include("test_differentiate.jl") include("test_crystal_plasticity.jl") -include("test_current_stress.jl") +include("test_stress_from_state.jl") # Test finite strain behaviors include("test_hyperelastic.jl") diff --git a/test/test_current_stress.jl b/test/test_stress_from_state.jl similarity index 85% rename from test/test_current_stress.jl rename to test/test_stress_from_state.jl index 8998067..3cefaeb 100644 --- a/test/test_current_stress.jl +++ b/test/test_stress_from_state.jl @@ -1,9 +1,9 @@ -@testset "calculate_current_stress" begin +@testset "stress_from_state" begin @testset "LinearElastic" begin m = LinearElastic(E=210.e3, ν=0.3) state = initial_material_state(m) ϵ = rand(SymmetricTensor{2,3}) - @test calculate_current_stress(m, ϵ, state) ≈ m.C ⊡ ϵ + @test stress_from_state(m, ϵ, state) ≈ m.C ⊡ ϵ # Reduced stress state: matches material_response's own (iterative) result, # and the analytical plane-stress relation @@ -12,7 +12,7 @@ rss = ReducedStressState(PlaneStress(), me) ϵ11 = 0.01 ϵ_red = SymmetricTensor{2,2}((ϵ11, 0.0, 0.0)) - σ_direct = calculate_current_stress(rss, ϵ_red, initial_material_state(rss)) + σ_direct = stress_from_state(rss, ϵ_red, initial_material_state(rss)) σ_mr, _, _, _ = material_response(PlaneStress(), me, ϵ_red, initial_material_state(me)) @test σ_direct ≈ σ_mr @test σ_direct[1, 1] ≈ E / (1 - ν^2) * ϵ11 @@ -28,12 +28,12 @@ state0 = initial_material_state(m) ϵ1 = SymmetricTensor{2,3}((i, j) -> (i, j) == (1, 1) ? 0.01 : 0.0) σ1, _, state1 = material_response(m, ϵ1, state0, nothing) - @test calculate_current_stress(m, ϵ1, state1) ≈ σ1 + @test stress_from_state(m, ϵ1, state1) ≈ σ1 # Frozen-state postprocessing: a different strain should give a purely # elastic increment from state1, NOT a fresh plastic correction. ϵ2 = ϵ1 + SymmetricTensor{2,3}((i, j) -> (i, j) == (1, 1) ? 0.02 : 0.0) - σ2_frozen = calculate_current_stress(m, ϵ2, state1) + σ2_frozen = stress_from_state(m, ϵ2, state1) @test σ2_frozen ≈ σ1 + e.C ⊡ (ϵ2 - ϵ1) σ2_true, _, state2_true = material_response(m, ϵ2, state1, nothing) @test !(σ2_true ≈ σ2_frozen) # material_response would further evolve plastically @@ -44,7 +44,7 @@ ϵ1_red = SymmetricTensor{2,2}((0.01, 0.0, 0.0)) state0_red = initial_material_state(rss) σ1_red, _, state1_red, _ = material_response(rss, ϵ1_red, state0_red, nothing) - σ1_red_current = calculate_current_stress(rss, ϵ1_red, state1_red) + σ1_red_current = stress_from_state(rss, ϵ1_red, state1_red) @test σ1_red_current ≈ σ1_red # PlaneStrain: verify that the non-iterative shortcut retains the transverse @@ -53,7 +53,7 @@ state0_strain = initial_material_state(rss_strain) σ1_strain, _, state1_strain, _ = material_response(rss_strain, ϵ1_red, state0_strain, nothing) @test state1_strain.ϵp[3, 3] != 0 # sanity: this test only matters if ϵp33 != 0 - σ1_strain_current = calculate_current_stress(rss_strain, ϵ1_red, state1_strain) + σ1_strain_current = stress_from_state(rss_strain, ϵ1_red, state1_strain) @test σ1_strain_current ≈ σ1_strain end @@ -65,13 +65,13 @@ state0 = initial_material_state(m) ϵ1 = rand(SymmetricTensor{2,3}) / 100 σ1, _, state1 = material_response(m, ϵ1, state0, 0.5) - @test calculate_current_stress(m, ϵ1, state1) ≈ σ1 + @test stress_from_state(m, ϵ1, state1) ≈ σ1 # Frozen-state: evaluating at a different strain must not re-solve the # viscous strain evolution (which requires Δt); it should be a pure # elastic-type increment using the given (fixed) viscous strain. ϵ2 = ϵ1 + rand(SymmetricTensor{2,3}) / 100 - σ2_frozen = calculate_current_stress(m, ϵ2, state1) + σ2_frozen = stress_from_state(m, ϵ2, state1) σ2_expected = MechMat.calculate_stress(me, ϵ2) + 2 * chain.G * (dev(ϵ2) - state1.ϵv[1]) @test σ2_frozen ≈ σ2_expected σ2_true, _, _ = material_response(m, ϵ2, state1, 0.5) @@ -84,7 +84,7 @@ ϵ1_red = SymmetricTensor{2,2}((0.01, 0.0, 0.0)) state0_red = initial_material_state(rss) σ1_red, _, state1_red, _ = material_response(rss, ϵ1_red, state0_red, 0.5) - @test calculate_current_stress(rss, ϵ1_red, state1_red) ≈ σ1_red + @test stress_from_state(rss, ϵ1_red, state1_red) ≈ σ1_red end @testset "RotatedMaterial" begin @@ -98,7 +98,7 @@ _, _, state1 = material_response(rm, ϵ_global1, state0, nothing) ϵ_global2 = ϵ_global1 + SymmetricTensor{2,3}((i, j) -> (i, j) == (1, 1) ? 0.001 : 0.0) - σ_current = calculate_current_stress(rm, ϵ_global2, state1) + σ_current = stress_from_state(rm, ϵ_global2, state1) θ = norm(r) ϵ_local2 = rotate(ϵ_global2, r, -θ) @@ -110,8 +110,8 @@ m_el = LinearElastic{:cubicsymmetry}(C1111=1 + rand(), C1122=1 + rand(), C1212=1 + rand()) rm_el = RotatedMaterial(m_el, r) ϵ = rand(SymmetricTensor{2,3}) - σ_rm_el = calculate_current_stress(rm_el, ϵ, initial_material_state(rm_el)) - σ_local_el = calculate_current_stress(m_el, rotate(ϵ, r, -θ), initial_material_state(m_el)) + σ_rm_el = stress_from_state(rm_el, ϵ, initial_material_state(rm_el)) + σ_local_el = stress_from_state(m_el, rotate(ϵ, r, -θ), initial_material_state(m_el)) @test σ_rm_el ≈ rotate(σ_local_el, r, θ) # Reduced stress state wrapping a rotated, stateful material: previously @@ -120,7 +120,7 @@ ϵ1_red = SymmetricTensor{2,2}((0.01, 0.0, 0.0)) state0_red = initial_material_state(rss) σ1_red, _, state1_red, _ = material_response(rss, ϵ1_red, state0_red, nothing) - @test calculate_current_stress(rss, ϵ1_red, state1_red) ≈ σ1_red + @test stress_from_state(rss, ϵ1_red, state1_red) ≈ σ1_red end @testset "HyperElastic" begin @@ -129,7 +129,7 @@ state = initial_material_state(m) F = one(Tensor{2,3}) + rand(Tensor{2,3}) / 20 P, _, _ = material_response(m, F, state) - @test calculate_current_stress(m, F, state) ≈ P + @test stress_from_state(m, F, state) ≈ P # Reduced stress state: this already works via the generic # `NoMaterialState` fallback, since MaterialModelsBase's stress-state @@ -139,7 +139,7 @@ F_red = one(Tensor{2,2}) + rand(Tensor{2,2}) / 20 state_red = initial_material_state(rss) P_red, _, _, _ = material_response(rss, F_red, state_red) - @test calculate_current_stress(rss, F_red, state_red) ≈ P_red + @test stress_from_state(rss, F_red, state_red) ≈ P_red end end @@ -152,13 +152,13 @@ state0 = initial_material_state(m) F1 = Tensor{2,3}((i, j) -> i == j ? (i == 1 ? 1.02 : 1.0) : 0.0) P1, _, state1 = material_response(m, F1, state0, nothing) - @test calculate_current_stress(m, F1, state1) ≈ P1 + @test stress_from_state(m, F1, state1) ≈ P1 @test state1.Fp != state0.Fp # sanity: this test only matters if plastic loading occurred # Frozen-state postprocessing: a different F should give the frozen-Fp # elastic response, NOT a fresh plastic correction. F2 = Tensor{2,3}((i, j) -> i == j ? (i == 1 ? 1.03 : 1.0) : 0.0) - σ2_frozen = calculate_current_stress(m, F2, state1) + σ2_frozen = stress_from_state(m, F2, state1) σ2_true, _, state2_true = material_response(m, F2, state1, nothing) @test !(σ2_true ≈ σ2_frozen) # material_response would further evolve plastically @test state2_true.Fp != state1.Fp @@ -168,6 +168,6 @@ F1_red = Tensor{2,2}((1.02, 0.0, 0.0, 1.0)) state0_red = initial_material_state(rss) P1_red, _, state1_red, _ = material_response(rss, F1_red, state0_red, nothing) - @test calculate_current_stress(rss, F1_red, state1_red) ≈ P1_red + @test stress_from_state(rss, F1_red, state1_red) ≈ P1_red end end From 501a024d8c17b6f9b60731cb3643ee63688c1b39 Mon Sep 17 00:00:00 2001 From: ClaudeBot Date: Mon, 21 Sep 2026 11:37:01 -0400 Subject: [PATCH 5/6] Remove redundant Plastic reduced-dimensional stress_from_state method MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Detailed change: MaterialModelsBase's generic reduced-dimensional fallback (autodiff-ing through the full-dimensional stress_from_state(m::Plastic, ϵ, state)) already gives the same result as the removed Plastic-specific fast path: the underlying formula, calculate_stress(m.elastic, ϵ - state.ϵp), is linear in ϵ, so the autodiff-derived tangent used in the generic fallback's Newton iteration is exact, same as m.elastic's own analytic stiffness. The dedicated method was therefore unnecessary duplication. Pkg.test() confirms all 29 stress_from_state tests still pass unchanged, including the PlaneStress/PlaneStrain reduced-dimensional Plastic tests, now exercising the generic fallback. docs/make.jl builds cleanly. Co-Authored-By: Claude Sonnet 5 --- src/StressFromState.jl | 17 ++++------------- test/test_stress_from_state.jl | 6 ++++-- 2 files changed, 8 insertions(+), 15 deletions(-) diff --git a/src/StressFromState.jl b/src/StressFromState.jl index 84ba99a..2cef5a7 100644 --- a/src/StressFromState.jl +++ b/src/StressFromState.jl @@ -27,23 +27,14 @@ MMB.stress_from_state(m::LinearElastic, ϵ::SymmetricTensor{2,3}, ::MMB.NoMateri MMB.stress_from_state(m::AbstractHyperElastic, F::Tensor{2,3}, ::MMB.NoMaterialState) = F ⋅ compute_stress(m, tdot(F)) # Plastic.jl +# Reduced-dimensional stress states (e.g. PlaneStress) need no dedicated method +# here: MaterialModelsBase's generic fallback autodiffs through this full-dim +# method and gives the same result (this formula is linear in ϵ, so the +# autodiff-derived tangent is exact, same as the elastic stiffness itself). function MMB.stress_from_state(m::Plastic, ϵ::SymmetricTensor{2,3}, state::PlasticState) return calculate_stress(m.elastic, ϵ - state.ϵp) end -function MMB.stress_from_state(stress_state::MMB.AbstractStressState, m::Plastic, ϵ, state::PlasticState) - # Expand the (possibly reduced) total strain to 3d before removing the plastic - # strain: for non-iterative states (e.g. PlaneStrain) the zero-padded - # out-of-plane *total* strain is exact by definition of the state, whereas - # reducing `state.ϵp` first would incorrectly discard its out-of-plane part. - # This avoids autodiff entirely, by delegating to `m.elastic`'s own analytic - # stress-state response. - ϵ_3d = MMB.expand_tensordim(stress_state, ϵ) - ϵₑ = ϵ_3d - state.ϵp - σ, _, _, _ = MMB.material_response(stress_state, m.elastic, ϵₑ, MMB.initial_material_state(m.elastic)) - return σ -end - # ViscoElastic.jl function MMB.stress_from_state(m::GeneralizedMaxwell, ϵ::SymmetricTensor{2,3}, state::GeneralizedMaxwellState) σ0 = calculate_stress(m.base, ϵ) diff --git a/test/test_stress_from_state.jl b/test/test_stress_from_state.jl index 3cefaeb..3476cbc 100644 --- a/test/test_stress_from_state.jl +++ b/test/test_stress_from_state.jl @@ -47,8 +47,10 @@ σ1_red_current = stress_from_state(rss, ϵ1_red, state1_red) @test σ1_red_current ≈ σ1_red - # PlaneStrain: verify that the non-iterative shortcut retains the transverse - # plastic strain's elastic coupling (regression check for issue found in review) + # PlaneStrain: verify the transverse plastic strain's elastic coupling is + # retained (regression check for an issue found in review of an earlier, + # since-removed Plastic-specific reduced-dimensional method; kept to + # confirm the generic fallback gets this right too). rss_strain = ReducedStressState(PlaneStrain(), m) state0_strain = initial_material_state(rss_strain) σ1_strain, _, state1_strain, _ = material_response(rss_strain, ϵ1_red, state0_strain, nothing) From 070572d2c7c14fcb7c947dcae7b18950121a6633 Mon Sep 17 00:00:00 2001 From: ClaudeBot Date: Mon, 21 Sep 2026 12:08:53 -0400 Subject: [PATCH 6/6] Relocate stress_from_state methods into their model files; fix test asserts MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Detailed change: Deletes src/StressFromState.jl and moves each material-specific MaterialModelsBase.stress_from_state method into its own model's existing source file, alongside its other AbstractMaterial interface implementations (material_response, initial_material_state, etc.): LinearElastic -> src/Elastic.jl, AbstractHyperElastic -> src/hyper_elasticity/HyperElastic.jl, Plastic -> src/Plastic.jl, GeneralizedMaxwell -> src/ViscoElastic.jl, FiniteStrainPlastic -> src/FiniteStrainPlastic.jl, RotatedMaterial (both methods plus the shared rotation helper) -> src/RotatedMaterial.jl. No behavioral change - same method bodies and signatures, just relocated, since implementing stress_from_state is part of the AbstractMaterial interface each model already implements, not a separate cross-cutting concern. Also addresses 2 minor review comments on test/test_stress_from_state.jl: changed two "sanity" checks (verifying the test setup itself produced plastic loading, not the feature under test) from @test x != y to @assert x ≉ y. @assert because a failure there means the test data is malformed, not a feature regression, so it shouldn't count toward the test summary; approximate inequality (not exact) since these compare floating-point tensors. The reviewer's literal suggested text used `!≈`, which is not valid Julia syntax (confirmed by a ParseError); used the correct, semantically equivalent operator `≉` (Base's negation of `≈`) instead. Went through the dual-review workflow (new Codex thread); both the plan review and the diff review after implementation reported no findings. Test results: Pkg.test() passes (27 tests in the stress_from_state testset, down from 29 as expected since two @tests became @asserts; full suite green). docs/make.jl builds cleanly. Co-Authored-By: Claude Sonnet 5 --- src/Elastic.jl | 5 +++ src/FiniteStrainPlastic.jl | 7 ++++ src/MechanicalMaterialModels.jl | 2 - src/Plastic.jl | 8 ++++ src/RotatedMaterial.jl | 12 ++++++ src/StressFromState.jl | 63 ---------------------------- src/ViscoElastic.jl | 5 +++ src/hyper_elasticity/HyperElastic.jl | 6 +++ test/test_stress_from_state.jl | 4 +- 9 files changed, 45 insertions(+), 67 deletions(-) delete mode 100644 src/StressFromState.jl diff --git a/src/Elastic.jl b/src/Elastic.jl index 7d8141b..d1c5888 100644 --- a/src/Elastic.jl +++ b/src/Elastic.jl @@ -93,6 +93,11 @@ end calculate_stress(m::LinearElastic, ϵ::SymmetricTensor) = m.C⊡ϵ +# Stress-only, no gradient (material_response would compute one, via `m.C`, +# that `stress_from_state`'s generic `NoMaterialState` fallback would +# otherwise compute via `material_response` and discard). +MMB.stress_from_state(m::LinearElastic, ϵ::SymmetricTensor{2,3}, ::MMB.NoMaterialState) = calculate_stress(m, ϵ) + # Functions for conversion between material and parameter vectors MMB.get_vector_length(::LinearElastic{<:Any,<:Any,N}) where{N} = N diff --git a/src/FiniteStrainPlastic.jl b/src/FiniteStrainPlastic.jl index 8d598a3..4e1c4d8 100644 --- a/src/FiniteStrainPlastic.jl +++ b/src/FiniteStrainPlastic.jl @@ -196,6 +196,13 @@ function calculate_PKstress(m::FiniteStrainPlastic, Fp::Tensor, F::Tensor) return P end +# `calculate_PKstress(m, state, F)` already computes the frozen-state (converged +# `state.Fp`, no Newton re-solve) 1st Piola-Kirchhoff stress; it is used above +# for the elastic-predictor branch of `material_response`. No +# reduced-dimensional method is needed: MaterialModelsBase's generic fallback +# (autodiff-ing through this method via its own `FrozenStressMaterial`) covers it. +MMB.stress_from_state(m::FiniteStrainPlastic, F::Tensor{2,3}, state::FiniteStrainPlasticState) = calculate_PKstress(m, state, F) + check_solution(x::FiniteStrainPlasticResidual) = x.Δλ < 0 ? throw(MMB.NoLocalConvergence("Plastic: Invalid solution, x.Δλ = ", x.Δλ, " < 0")) : nothing # TODO: Could be replaced by exponential map. diff --git a/src/MechanicalMaterialModels.jl b/src/MechanicalMaterialModels.jl index 8b01ddc..8a6ccff 100644 --- a/src/MechanicalMaterialModels.jl +++ b/src/MechanicalMaterialModels.jl @@ -57,6 +57,4 @@ export NeoHooke, CompressibleNeoHooke, SaintVenant include("FiniteStrainPlastic.jl") export FiniteStrainPlastic -include("StressFromState.jl") - end diff --git a/src/Plastic.jl b/src/Plastic.jl index 3136ec2..2131010 100644 --- a/src/Plastic.jl +++ b/src/Plastic.jl @@ -171,6 +171,14 @@ end check_solution(x::PlasticResidual) = x.Δλ < 0 ? throw(MMB.NoLocalConvergence("Plastic: Invalid solution, x.Δλ = ", x.Δλ, " < 0")) : nothing +# Reduced-dimensional stress states (e.g. PlaneStress) need no dedicated method +# here: MaterialModelsBase's generic fallback autodiffs through this full-dim +# method and gives the same result (this formula is linear in ϵ, so the +# autodiff-derived tangent is exact, same as the elastic stiffness itself). +function MMB.stress_from_state(m::Plastic, ϵ::SymmetricTensor{2,3}, state::PlasticState) + return calculate_stress(m.elastic, ϵ - state.ϵp) +end + # General residual function function residual(x::PlasticResidual{NKin,NIso}, m::Plastic, old::PlasticState, ϵ, Δt, cache) where{NKin,NIso} σ_red = x.σ - sum(x.β) diff --git a/src/RotatedMaterial.jl b/src/RotatedMaterial.jl index 650def6..f863ece 100644 --- a/src/RotatedMaterial.jl +++ b/src/RotatedMaterial.jl @@ -33,3 +33,15 @@ function MMB.material_response(rm::RotatedMaterial, strain::AbstractTensor, args stiff = rotate(stiff_rot, rm.rotation, θ) return stress, stiff, state end + +function _stress_from_state_rotated(rm::RotatedMaterial, ϵ::SymmetricTensor{2,3}, state) + θ = norm(rm.rotation) + ϵ_rot = rotate(ϵ, rm.rotation, -θ) + σ_rot = MMB.stress_from_state(rm.material, ϵ_rot, state) + return rotate(σ_rot, rm.rotation, θ) +end +MMB.stress_from_state(rm::RotatedMaterial, ϵ::SymmetricTensor{2,3}, state) = _stress_from_state_rotated(rm, ϵ, state) +# Disambiguates against MaterialModelsBase's own `(AbstractMaterial, ϵ, +# ::NoMaterialState)` fallback, which would otherwise be equally specific when +# `rm.material` is stateless. +MMB.stress_from_state(rm::RotatedMaterial, ϵ::SymmetricTensor{2,3}, state::MMB.NoMaterialState) = _stress_from_state_rotated(rm, ϵ, state) diff --git a/src/StressFromState.jl b/src/StressFromState.jl deleted file mode 100644 index 2cef5a7..0000000 --- a/src/StressFromState.jl +++ /dev/null @@ -1,63 +0,0 @@ -# Material-specific implementations of `MaterialModelsBase.stress_from_state` -# (https://github.com/KnutAM/MaterialModelsBase.jl/pull/21), which upstreamed the -# generic postprocessing machinery originally prototyped in this file (see -# https://github.com/KnutAM/MechanicalMaterialModels.jl/pull/13): the generic -# `NoMaterialState` fallbacks, `FrozenStressMaterial`, the generic -# reduced-dimensional fallback, and `ReducedStressState` delegation now all live -# in `MaterialModelsBase.jl` itself. Only the material-specific methods below - -# either required (no `NoMaterialState`-based default exists) or a cheaper, -# non-autodiff alternative to the generic fallback - remain here. -# -# NOTE: while MaterialModelsBase.jl#21 is not yet merged, Project.toml and -# docs/Project.toml temporarily point MaterialModelsBase at the branch -# implementing it (knutambot/MaterialModelsBase.jl#cb/calculate_current_stress). -# Revert both `[sources]` entries (and tighten the `MaterialModelsBase` `[compat]` -# bound to the first release containing `stress_from_state`) once that PR merges -# and is released. - -# LinearElastic.jl: stress-only, no gradient (material_response would compute one, -# via `m.C`, that `stress_from_state`'s generic `NoMaterialState` fallback would -# otherwise compute via `material_response` and discard). -MMB.stress_from_state(m::LinearElastic, ϵ::SymmetricTensor{2,3}, ::MMB.NoMaterialState) = calculate_stress(m, ϵ) - -# HyperElastic.jl: stress-only. Computing `S = 2 ∂Ψ/∂C` (once) is unavoidable to get -# the stress at all, but the generic `NoMaterialState` fallback's `material_response` -# call additionally differentiates through that once more to get the tangent, which -# `stress_from_state` doesn't need. -MMB.stress_from_state(m::AbstractHyperElastic, F::Tensor{2,3}, ::MMB.NoMaterialState) = F ⋅ compute_stress(m, tdot(F)) - -# Plastic.jl -# Reduced-dimensional stress states (e.g. PlaneStress) need no dedicated method -# here: MaterialModelsBase's generic fallback autodiffs through this full-dim -# method and gives the same result (this formula is linear in ϵ, so the -# autodiff-derived tangent is exact, same as the elastic stiffness itself). -function MMB.stress_from_state(m::Plastic, ϵ::SymmetricTensor{2,3}, state::PlasticState) - return calculate_stress(m.elastic, ϵ - state.ϵp) -end - -# ViscoElastic.jl -function MMB.stress_from_state(m::GeneralizedMaxwell, ϵ::SymmetricTensor{2,3}, state::GeneralizedMaxwellState) - σ0 = calculate_stress(m.base, ϵ) - return mapreduce((c, ϵv) -> 2 * c.G * (dev(ϵ) - ϵv), +, m.chains, state.ϵv; init=σ0) -end - -# FiniteStrainPlastic.jl -# `calculate_PKstress(m, state, F)` already computes the frozen-state (converged -# `state.Fp`, no Newton re-solve) 1st Piola-Kirchhoff stress; it is used -# internally for the elastic-predictor branch of `material_response`. No -# reduced-dimensional method is needed: MaterialModelsBase's generic fallback -# (autodiff-ing through this method via its own `FrozenStressMaterial`) covers it. -MMB.stress_from_state(m::FiniteStrainPlastic, F::Tensor{2,3}, state::FiniteStrainPlasticState) = calculate_PKstress(m, state, F) - -# RotatedMaterial.jl -function _stress_from_state_rotated(rm::RotatedMaterial, ϵ::SymmetricTensor{2,3}, state) - θ = norm(rm.rotation) - ϵ_rot = rotate(ϵ, rm.rotation, -θ) - σ_rot = MMB.stress_from_state(rm.material, ϵ_rot, state) - return rotate(σ_rot, rm.rotation, θ) -end -MMB.stress_from_state(rm::RotatedMaterial, ϵ::SymmetricTensor{2,3}, state) = _stress_from_state_rotated(rm, ϵ, state) -# Disambiguates against MaterialModelsBase's own `(AbstractMaterial, ϵ, -# ::NoMaterialState)` fallback, which would otherwise be equally specific when -# `rm.material` is stateless. -MMB.stress_from_state(rm::RotatedMaterial, ϵ::SymmetricTensor{2,3}, state::MMB.NoMaterialState) = _stress_from_state_rotated(rm, ϵ, state) diff --git a/src/ViscoElastic.jl b/src/ViscoElastic.jl index bce9b2b..cac874a 100644 --- a/src/ViscoElastic.jl +++ b/src/ViscoElastic.jl @@ -71,3 +71,8 @@ function MMB.material_response(m::GeneralizedMaxwell, ϵ::SymmetricTensor{2,3}, state = GeneralizedMaxwellState(map((c, ϵv_old) -> calculate_viscous_strain(c, ϵ, ϵv_old, Δt), m.chains, old.ϵv)) return σ, dσdϵ, state end + +function MMB.stress_from_state(m::GeneralizedMaxwell, ϵ::SymmetricTensor{2,3}, state::GeneralizedMaxwellState) + σ0 = calculate_stress(m.base, ϵ) + return mapreduce((c, ϵv) -> 2 * c.G * (dev(ϵ) - ϵv), +, m.chains, state.ϵv; init=σ0) +end diff --git a/src/hyper_elasticity/HyperElastic.jl b/src/hyper_elasticity/HyperElastic.jl index ef69b7e..106e433 100644 --- a/src/hyper_elasticity/HyperElastic.jl +++ b/src/hyper_elasticity/HyperElastic.jl @@ -46,3 +46,9 @@ function MMB.material_response(m::AbstractHyperElastic, F::Tensor{2,3}, old::Abs end MMB.get_tensorbase(::AbstractHyperElastic) = Tensor{2,3} + +# Stress-only. Computing `S = 2 ∂Ψ/∂C` (once) is unavoidable to get the stress +# at all, but the generic `NoMaterialState` fallback's `material_response` +# call additionally differentiates through that once more to get the tangent, +# which `stress_from_state` doesn't need. +MMB.stress_from_state(m::AbstractHyperElastic, F::Tensor{2,3}, ::MMB.NoMaterialState) = F ⋅ compute_stress(m, tdot(F)) diff --git a/test/test_stress_from_state.jl b/test/test_stress_from_state.jl index 3476cbc..2d34801 100644 --- a/test/test_stress_from_state.jl +++ b/test/test_stress_from_state.jl @@ -155,7 +155,7 @@ F1 = Tensor{2,3}((i, j) -> i == j ? (i == 1 ? 1.02 : 1.0) : 0.0) P1, _, state1 = material_response(m, F1, state0, nothing) @test stress_from_state(m, F1, state1) ≈ P1 - @test state1.Fp != state0.Fp # sanity: this test only matters if plastic loading occurred + @assert state1.Fp ≉ state0.Fp # sanity: this test only matters if plastic loading occurred # Frozen-state postprocessing: a different F should give the frozen-Fp # elastic response, NOT a fresh plastic correction. @@ -163,7 +163,7 @@ σ2_frozen = stress_from_state(m, F2, state1) σ2_true, _, state2_true = material_response(m, F2, state1, nothing) @test !(σ2_true ≈ σ2_frozen) # material_response would further evolve plastically - @test state2_true.Fp != state1.Fp + @assert state2_true.Fp ≉ state1.Fp # Reduced stress state, via the FrozenStressMaterial + MMB stress-state iteration rss = ReducedStressState(PlaneStress(), m)