Add stress_from_state interface for postprocessing

docs/Project.toml

@@ -3,3 +3,6 @@ CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 MaterialModelsBase = "af893363-701d-44dc-8b1e-d9a2c129bfc9"
 Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4"
+
+[sources]
+MaterialModelsBase = {path = ".."}

docs/src/index.md

@@ -104,6 +104,14 @@ given as
 AbstractExtraOutput
 ```
 
+### Postprocessing
+When postprocessing an already converged simulation, it is often useful to calculate
+the stress again without advancing any history variables.
+`stress_from_state` provides this and supports lower-dimensional stress states as well.
+```@docs
+stress_from_state
+```
+
 ### Exceptions
 Finally, the following exceptions are included
 ```@docs

docs/src/stressiterations.md

@@ -33,4 +33,4 @@ state and the material type into a single type that is passed to the element rou
 The wrapper `ReducedStressState` is provided for that purpose.
 ```@docs
 ReducedStressState
-```
+```

src/MaterialModelsBase.jl

@@ -22,6 +22,9 @@ export PlaneStress, UniaxialStress, UniaxialNormalStress    # Iterative stress s
 export GeneralStressState                                   # General iterative 3D non-zero stress state
 export update_stress_state!                                 # For nonzero stress-conditions
 
+# Postprocessing
+export stress_from_state                                    # Stress consistent with a frozen state
+
 # For parameter identification and differentiation of materials
 export tovector, tovector!, fromvector                      # Convert to/from `AbstractVector`s
 export get_num_tensorcomponents, get_num_statevars          # Information about the specific material
@@ -158,6 +161,7 @@ struct NoExtraOutput <: AbstractExtraOutput end
 
 include("vector_conversion.jl")
 include("stressiterations.jl")
+include("stress_from_state.jl")
 include("differentiation.jl")
 include("ErrorExceptions.jl")
 

src/stress_from_state.jl

@@ -0,0 +1,80 @@
+"""
+    stress_from_state(m::AbstractMaterial, strain, state::AbstractMaterialState)
+    stress_from_state(stress_state::AbstractStressState, m::AbstractMaterial, strain, state::AbstractMaterialState)
+    stress_from_state(rss::ReducedStressState, strain, state::AbstractMaterialState)
+
+## Using this interface
+Calculate the stress that is energy-conjugated to `strain`, 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.
+
+!!! warning
+    Differentiating this function wrt. `strain` while holding `state` fixed gives a
+    frozen-state tangent, which generally differs from `material_response`'s consistent
+    tangent whenever internal/history variables would evolve with `strain` (e.g. for a
+    plastic or viscous material). It is therefore not a general replacement for
+    `material_response` during, e.g., equilibrium iterations.
+
+## Implementing this interface
+A material-model developer only needs to implement the full-dimensional method,
+`stress_from_state(m::MyMaterial, strain, state::MyMaterialState)`. If `MyMaterial`
+has no state (i.e. `initial_material_state(m) isa NoMaterialState`), this is not
+required either, since [`material_response`](@ref) is then already frozen-state by
+definition and a generic fallback is provided. Support for a reduced-dimensional stress 
+state (e.g. via [`ReducedStressState`](@ref)) then follows automatically from a generic 
+fallback, using the tangent obtained by automatic differentiation via `Tensors.gradient`. 
+A specific reduced-dimensional method, 
+`stress_from_state(stress_state::AbstractStressState, m::MyMaterial, strain, state::MyMaterialState)`,
+can be added when a cheaper, non-autodiff alternative exists.
+"""
+function stress_from_state end
+
+# Fully generic: a material with no state has, by definition, nothing to freeze -
+# `material_response` already gives the frozen-state stress.
+function stress_from_state(m::AbstractMaterial, strain, state::NoMaterialState)
+    σ, _, _ = material_response(m, strain, state)
+    return σ
+end
+
+# Wraps a material `m` and a frozen state `s` as an `AbstractMaterial`, whose
+# `material_response` evaluates `stress_from_state(m, strain, s)` (at fixed
+# history/internal variables) so that it can ride the existing stress-state Newton
+# iteration (e.g. for `PlaneStress`). The tangent needed for that iteration is
+# obtained via automatic differentiation. This powers the generic reduced-dimensional
+# fallback of `stress_from_state` below.
+struct FrozenStressMaterial{MT <: AbstractMaterial, ST <: AbstractMaterialState} <: AbstractMaterial
+    m::MT
+    s::ST
+end
+function material_response(fm::FrozenStressMaterial, strain::SecondOrderTensor{3}, old::AbstractMaterialState, args...)
+    dσdϵ, σ = Tensors.gradient(e -> stress_from_state(fm.m, e, fm.s), strain, :all)
+    return σ, dσdϵ, old
+end
+
+# Generic reduced-dimensional fallback: as long as `stress_from_state(m, strain,
+# state)` (full-dimensional) is implemented for `m`, this makes `ReducedStressState`
+# support "just work", by autodiff-ing through it. The `NoMaterialState` fast path
+# below takes precedence when a cheaper, non-autodiff alternative exists.
+function stress_from_state(stress_state::AbstractStressState, m::AbstractMaterial, strain, state::AbstractMaterialState)
+    frozen = FrozenStressMaterial(m, state)
+    σ, _, _, _ = material_response(stress_state, frozen, strain, 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
+# (potentially analytic) stress-state handling.
+function stress_from_state(stress_state::AbstractStressState, m::AbstractMaterial, strain, state::NoMaterialState)
+    return first(material_response(stress_state, m, strain, state))
+end
+
+function stress_from_state(rss::ReducedStressState, strain, state::AbstractMaterialState)
+    return stress_from_state(rss.stress_state, rss.material, strain, state)
+end
+
+# Disambiguates the two 3-argument methods above for a `ReducedStressState` wrapping a
+# stateless material.
+function stress_from_state(rss::ReducedStressState, strain, state::NoMaterialState)
+    return stress_from_state(rss.stress_state, rss.material, strain, state)
+end

test/runtests.jl

@@ -15,6 +15,7 @@ include("utils4testing.jl")
 
 include("vector_conversion.jl")
 include("stressiterations.jl")
+include("stress_from_state.jl")
 include("differentiation.jl")
 include("errors.jl")
 include("performance.jl")

test/stress_from_state.jl

@@ -0,0 +1,188 @@
+module CurrentStressTestMaterials
+using MaterialModelsBase
+using Tensors
+import MaterialModelsBase as MMB
+
+# Well-conditioned isotropic elastic stiffness tensors, so that the `PlaneStress`
+# Newton iteration used below converges reliably (unlike an arbitrary random
+# 4th-order tensor, which need not be positive definite).
+function isotropic_C(G, K)
+    I2 = one(SymmetricTensor{2,3})
+    I4vol = otimes(I2, I2) / 3
+    I4dev = one(SymmetricTensor{4,3}) - I4vol
+    return 2G * I4dev + 3K * I4vol
+end
+function isotropic_C_finite(G, K)
+    I2 = one(Tensor{2,3})
+    return 2G * otimesu(I2, I2) + K * otimes(I2, I2)
+end
+
+# Stateless (NoMaterialState) toy material: exercises the fully generic
+# NoMaterialState fallbacks without any material-specific
+# `stress_from_state` method at all.
+struct ToyElastic{T} <: AbstractMaterial
+    C::SymmetricTensor{4,3,T}
+end
+MMB.initial_material_state(m::ToyElastic{T}) where {T} = MMB.NoMaterialState{T}()
+function MMB.material_response(m::ToyElastic, ϵ::SymmetricTensor{2,3}, state, args...)
+    return m.C ⊡ ϵ, m.C, state
+end
+
+# Same as `ToyElastic`, but with a specialized `PlaneStress` method that omits the
+# optional 4th (full-strain) output, as explicitly permitted by the
+# `material_response(::AbstractStressState, ...)` interface. Regression test for a
+# `stress_from_state` fast path that would otherwise assume 4 outputs.
+struct ToyElasticSpecialized{T} <: AbstractMaterial
+    C::SymmetricTensor{4,3,T}
+end
+MMB.initial_material_state(m::ToyElasticSpecialized{T}) where {T} = MMB.NoMaterialState{T}()
+function MMB.material_response(m::ToyElasticSpecialized, ϵ::SymmetricTensor{2,3}, state, args...)
+    return m.C ⊡ ϵ, m.C, state
+end
+function MMB.material_response(::PlaneStress, m::ToyElasticSpecialized, ϵ::SymmetricTensor{2,2}, state, args...)
+    C_red = SymmetricTensor{4,2}((i, j, k, l) -> m.C[i, j, k, l])
+    return C_red ⊡ ϵ, C_red, state
+end
+
+# Small-strain toy material with a state that unconditionally "evolves" every call
+# (unlike real plasticity with a yield surface), so that a frozen-state evaluation
+# at a different strain is guaranteed to differ from a fresh `material_response` call.
+struct ToyHistory{T} <: AbstractMaterial
+    C::SymmetricTensor{4,3,T}
+    H::T
+end
+struct ToyHistoryState{T} <: AbstractMaterialState
+    ϵp::SymmetricTensor{2,3,T}
+end
+MMB.initial_material_state(m::ToyHistory{T}) where {T} = ToyHistoryState(zero(SymmetricTensor{2,3,T}))
+function MMB.material_response(m::ToyHistory, ϵ::SymmetricTensor{2,3}, state::ToyHistoryState, args...)
+    ϵp_new = state.ϵp + m.H * (ϵ - state.ϵp)
+    σ = m.C ⊡ (ϵ - ϵp_new)
+    dσdϵ = (1 - m.H) * m.C # Consistent tangent: ϵp_new depends linearly on ϵ
+    return σ, dσdϵ, ToyHistoryState(ϵp_new)
+end
+MMB.stress_from_state(m::ToyHistory, ϵ::SymmetricTensor{2,3}, state::ToyHistoryState) = m.C ⊡ (ϵ - state.ϵp)
+
+# Finite-strain (nonsymmetric `Tensor{2,3}`) analogue, to check that
+# `FrozenStressMaterial`'s autodiff also works for the `Tensor` tensor family.
+struct ToyHistoryFinite{T} <: AbstractMaterial
+    C::Tensor{4,3,T}
+    H::T
+end
+struct ToyHistoryFiniteState{T} <: AbstractMaterialState
+    Fp::Tensor{2,3,T}
+end
+MMB.initial_material_state(m::ToyHistoryFinite{T}) where {T} = ToyHistoryFiniteState(zero(Tensor{2,3,T}))
+function MMB.material_response(m::ToyHistoryFinite, F::Tensor{2,3}, state::ToyHistoryFiniteState, args...)
+    Fp_new = state.Fp + m.H * (F - state.Fp)
+    P = m.C ⊡ (F - Fp_new)
+    dPdF = (1 - m.H) * m.C # Consistent tangent: Fp_new depends linearly on F
+    return P, dPdF, ToyHistoryFiniteState(Fp_new)
+end
+MMB.stress_from_state(m::ToyHistoryFinite, F::Tensor{2,3}, state::ToyHistoryFiniteState) = m.C ⊡ (F - state.Fp)
+
+# A plain, stateless material equivalent to `ToyHistory` frozen at a given `ϵp`, used
+# as an independent reference to check the generic reduced-dimensional fallback
+# (which goes through `FrozenStressMaterial` and autodiff) against MMB's own,
+# already-tested, `PlaneStress` Newton iteration on an explicit material.
+struct FrozenLinear{T} <: AbstractMaterial
+    C::SymmetricTensor{4,3,T}
+    ϵp::SymmetricTensor{2,3,T}
+end
+MMB.initial_material_state(m::FrozenLinear{T}) where {T} = MMB.NoMaterialState{T}()
+function MMB.material_response(m::FrozenLinear, ϵ::SymmetricTensor{2,3}, state, args...)
+    return m.C ⊡ (ϵ - m.ϵp), m.C, state
+end
+
+end # module
+
+import .CurrentStressTestMaterials as CT
+
+@testset "stress_from_state" begin
+    @testset "NoMaterialState generic fallback" begin
+        C = CT.isotropic_C(80.e3, 160.e3)
+        m = CT.ToyElastic(C)
+        state = initial_material_state(m)
+        @test state isa MaterialModelsBase.NoMaterialState
+        ϵ = rand(SymmetricTensor{2,3})
+        @test stress_from_state(m, ϵ, state) ≈ C ⊡ ϵ
+
+        # Reduced-dimensional fast path (no material-specific method exists at all)
+        rss = ReducedStressState(PlaneStress(), m)
+        ϵ_red = rand(SymmetricTensor{2,2})
+        state_red = initial_material_state(rss)
+        σ_direct = stress_from_state(rss, ϵ_red, state_red)
+        σ_mr, _, _, _ = material_response(rss, ϵ_red, state_red)
+        @test σ_direct ≈ σ_mr
+
+        # Regression test: a specialized `material_response(stress_state, m, ...)`
+        # method is allowed to omit the optional 4th (full-strain) output.
+        m_spec = CT.ToyElasticSpecialized(C)
+        state_spec = initial_material_state(m_spec)
+        σ_spec_direct = stress_from_state(PlaneStress(), m_spec, ϵ_red, state_spec)
+        σ_spec_mr, _, _ = material_response(PlaneStress(), m_spec, ϵ_red, state_spec)
+        @test σ_spec_direct ≈ σ_spec_mr
+    end
+
+    @testset "Stateful material, full dimension" begin
+        C = CT.isotropic_C(80.e3, 160.e3)
+        m = CT.ToyHistory(C, 0.5)
+
+        state0 = initial_material_state(m)
+        ϵ1 = rand(SymmetricTensor{2,3})
+        σ1, _, state1 = material_response(m, ϵ1, state0)
+        @test stress_from_state(m, ϵ1, state1) ≈ σ1
+
+        # Frozen-state postprocessing: a different strain should give the frozen-ϵp
+        # response, NOT a fresh history update.
+        ϵ2 = ϵ1 + rand(SymmetricTensor{2,3}) / 10
+        σ2_frozen = stress_from_state(m, ϵ2, state1)
+        @test σ2_frozen ≈ C ⊡ (ϵ2 - state1.ϵp)
+        σ2_true, _, state2_true = material_response(m, ϵ2, state1)
+        @test !(σ2_true ≈ σ2_frozen)
+        @test state2_true.ϵp != state1.ϵp
+    end
+
+    @testset "Stateful material, generic reduced-dimensional fallback" begin
+        C = CT.isotropic_C(80.e3, 160.e3)
+        m = CT.ToyHistory(C, 0.5)
+        rss = ReducedStressState(PlaneStress(), m)
+
+        state0 = initial_material_state(rss)
+        ϵ1 = rand(SymmetricTensor{2,2}) / 10
+        σ1, _, state1, _ = material_response(rss, ϵ1, state0)
+        @test stress_from_state(rss, ϵ1, state1) ≈ σ1
+
+        ϵ2 = ϵ1 + rand(SymmetricTensor{2,2}) / 10
+        σ2_frozen = stress_from_state(rss, ϵ2, state1)
+
+        # Independent reference: the frozen material is exactly linear elastic in
+        # (ϵ - state1.ϵp), so its plane-stress response can be obtained directly
+        # from MMB's own (already-tested) stress-state iteration on an explicit,
+        # equivalent material, bypassing `stress_from_state` entirely.
+        flin = CT.FrozenLinear(C, state1.ϵp)
+        σ2_expected, _, _, _ = material_response(PlaneStress(), flin, ϵ2, initial_material_state(flin))
+        @test σ2_frozen ≈ σ2_expected
+
+        σ2_true, _, state2_true, _ = material_response(rss, ϵ2, state1)
+        @test !(σ2_true ≈ σ2_frozen)
+        @test state2_true.ϵp != state1.ϵp
+    end
+
+    @testset "Finite-strain (Tensor) frozen-state fallback" begin
+        C = CT.isotropic_C_finite(80.e3, 160.e3)
+        m = CT.ToyHistoryFinite(C, 0.5)
+        rss = ReducedStressState(PlaneStress(), m)
+
+        state0 = initial_material_state(rss)
+        F1 = one(Tensor{2,2}) + rand(Tensor{2,2}) / 20
+        P1, _, state1, _ = material_response(rss, F1, state0)
+        @test stress_from_state(rss, F1, state1) ≈ P1
+
+        F2 = F1 + rand(Tensor{2,2}) / 20
+        P2_frozen = stress_from_state(rss, F2, state1)
+        P2_true, _, state2_true, _ = material_response(rss, F2, state1)
+        @test !(P2_true ≈ P2_frozen)
+        @test state2_true.Fp != state1.Fp
+    end
+end