Add stress_from_state material-specific implementations for MaterialModelsBase.jl#12

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"}

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"}

docs/src/finite_strains.md

@@ -101,3 +101,9 @@ Otherwise, the overstress function, ``\eta(\varPhi)``, determines the evolution
 ```math
 \dot{\lambda} = \eta(\varPhi, Y_0 + \kappa)
 ```
+
+## [Postprocessing](@id finite_strain_postprocessing)
+`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).

docs/src/small_strains.md

@@ -86,3 +86,22 @@ BCC12
 GenericCrystallography
 CrystalPlasticity
 ```
+
+## [Postprocessing](@id small_strain_postprocessing)
+`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.

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
 

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. 

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.β)

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)

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

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))

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_stress_from_state.jl")
 
 # Test finite strain behaviors
 include("test_hyperelastic.jl")

test/test_stress_from_state.jl

@@ -0,0 +1,175 @@
+@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 stress_from_state(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 = 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
+        @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 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 = 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
+        @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 = stress_from_state(rss, ϵ1_red, state1_red)
+        @test σ1_red_current ≈ σ1_red
+
+        # 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)
+        @test state1_strain.ϵp[3, 3] != 0 # sanity: this test only matters if ϵp33 != 0
+        σ1_strain_current = stress_from_state(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 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 = 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)
+        @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 stress_from_state(rss, ϵ1_red, state1_red) ≈ σ1_red
+    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 = stress_from_state(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 = 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
+        # 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 stress_from_state(rss, ϵ1_red, state1_red) ≈ σ1_red
+    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 stress_from_state(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 stress_from_state(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 stress_from_state(m, F1, state1) ≈ P1
+        @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.
+        F2 = Tensor{2,3}((i, j) -> i == j ? (i == 1 ? 1.03 : 1.0) : 0.0)
+        σ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
+        @assert 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 stress_from_state(rss, F1_red, state1_red) ≈ P1_red
+    end
+end