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 4454512..8809f2c 100644 --- a/docs/src/finite_strains.md +++ b/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). diff --git a/docs/src/small_strains.md b/docs/src/small_strains.md index 24c1ea3..501d953 100644 --- a/docs/src/small_strains.md +++ b/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. 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/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/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/runtests.jl b/test/runtests.jl index a0e18f8..5c0ccea 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_stress_from_state.jl") # Test finite strain behaviors include("test_hyperelastic.jl") diff --git a/test/test_stress_from_state.jl b/test/test_stress_from_state.jl new file mode 100644 index 0000000..2d34801 --- /dev/null +++ b/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