From 583c6377a52359cc10fde41c6c2752bf52208fb4 Mon Sep 17 00:00:00 2001 From: Knut Andreas Date: Tue, 26 Aug 2025 23:07:56 +0200 Subject: [PATCH 01/14] Initializing improved differentiation interface Need to add more tests as no changes was required, specifically for diff with state variables + stress iterations --- Project.toml | 8 +- src/MaterialModelsBase.jl | 3 +- src/differentiation.jl | 92 +++++++- src/stressiterations.jl | 28 ++- src/vector_conversion.jl | 4 + test/TestMaterials/Manifest.toml | 216 ++++++++++++++++++ test/TestMaterials/Project.toml | 15 ++ test/{ => TestMaterials/src}/TestMaterials.jl | 7 +- test/differentiation.jl | 21 ++ test/runtests.jl | 6 +- 10 files changed, 376 insertions(+), 24 deletions(-) create mode 100644 test/TestMaterials/Manifest.toml create mode 100644 test/TestMaterials/Project.toml rename test/{ => TestMaterials/src}/TestMaterials.jl (95%) diff --git a/Project.toml b/Project.toml index 3b86dc8..146a130 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ authors = ["Knut Andreas Meyer and contributors"] version = "0.2.3" [deps] +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4" @@ -11,8 +12,11 @@ Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4" julia = "1.8" [extras] -ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +TestMaterials = "aa4f4413-326f-4ddd-a53e-c58801ebfeaf" + +[sources] +TestMaterials = {path = "test/TestMaterials"} [targets] -test = ["Test", "ForwardDiff"] +test = ["Test", "TestMaterials"] diff --git a/src/MaterialModelsBase.jl b/src/MaterialModelsBase.jl index 3d4a00b..0a01bc7 100644 --- a/src/MaterialModelsBase.jl +++ b/src/MaterialModelsBase.jl @@ -1,5 +1,6 @@ module MaterialModelsBase using Tensors, StaticArrays +using ForwardDiff: ForwardDiff import Base: @kwdef # General interface for <:AbstractMaterial @@ -153,8 +154,8 @@ abstract type AbstractExtraOutput end struct NoExtraOutput <: AbstractExtraOutput end include("vector_conversion.jl") -include("differentiation.jl") include("stressiterations.jl") +include("differentiation.jl") include("ErrorExceptions.jl") end \ No newline at end of file diff --git a/src/differentiation.jl b/src/differentiation.jl index 9e7db4a..e4702e9 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -12,8 +12,12 @@ end A struct that saves all derivative information using a `Matrix{T}` for each derivative, where `T=get_params_eltype(m)`. The dimensions are obtained from `get_num_tensorcomponents`, -`get_num_statevars`, and `get_num_params`. The values should be updated in `differentiate_material!` -by direct access of the fields, where `σ` is the stress, `ϵ` the strain, `s` and `ⁿs` are the current +`get_num_statevars`, and `get_num_params`. `m` must support `tovector` and `fromvector`, while +the output of `initial_material_state` must support `tovector`, and in addition the element type +of `tovector(initial_material_state(m))` must respect the element type in `tovector(m)` for any `m`. + +The values should be updated in `differentiate_material!` by direct access of the fields, +where `σ` is the stress, `ϵ` the strain, `s` and `ⁿs` are the current and old state variables, and `p` the material parameter vector. * `dσdϵ` @@ -29,13 +33,14 @@ function MaterialDerivatives(m::AbstractMaterial) n_tensor = get_num_tensorcomponents(m) n_state = get_num_statevars(m) n_params = get_num_params(m) + dsdp = ForwardDiff.jacobian(p -> tovector(initial_material_state(fromvector(p, m))), tovector(m)) return MaterialDerivatives( zeros(T, n_tensor, n_tensor), # dσdϵ zeros(T, n_tensor, n_state), # dσdⁿs zeros(T, n_tensor, n_params), # dσdp zeros(T, n_state, n_tensor), # dsdϵ zeros(T, n_state, n_state), # dsdⁿs - zeros(T, n_state, n_params) # dsdp + dsdp ) end @@ -66,4 +71,83 @@ allocate_differentiation_output(::AbstractMaterial) = NoExtraOutput() Calculate the derivatives and save them in `diff`, see [`MaterialDerivatives`](@ref) for a description of the fields in `diff`. """ -function differentiate_material! end \ No newline at end of file +function differentiate_material! end + +struct StressStateDerivatives{T} + mderiv::MaterialDerivatives{T} + dϵdp::Matrix{T} + dσdp::Matrix{T} + ϵindex::SMatrix{3, 3, Int} # To allow indexing by (i, j) into only + σindex::SMatrix{3, 3, Int} # saved values to avoid storing unused rows. + # TODO: Reduce the dimensions, for now all entries (even those that are zero) are stored. +end + +""" + differentiate_material!(ssd::StressStateDerivatives, stress_state, m, args...) + +For material models implementing `material_response(m, args...)` and `differentiate_material!(::MaterialDerivatives, m, args...)`, +this method will work automatically by +1) Calling `σ, dσdϵ, state = material_response(stress_state, m, args...)` (except that `dσdϵ::FourthOrderTensor{dim = 3}` is extracted) +2) Calling `differentiate_material!(ssd.mderiv::MaterialDerivatives, m, args..., dσdϵ::FourthOrderTensor{3})` +3) Updating `ssd` according to the constraints imposed by the `stress_state`. + +For material models that directly implement `material_response(stress_state, m, args...)`, this function should be overloaded directly +to calculate the derivatives in `ssd`. Here the user has full control and no modifications occur automatically, however, typically the +(total) derivatives `ssd.dσdp`, `ssd.dϵdp`, and `ssd.mderiv.dsdp` should be updated. +""" +function differentiate_material!(ssd::StressStateDerivatives, stress_state::AbstractStressState, m::AbstractMaterial, args::Vararg{Any,N}) where {N} + σ_full, dσdϵ_full, state, ϵ_full = stress_state_material_response(stress_state, m, args...) + differentiate_material!(ssd.mderiv, m, args..., dσdϵ_full) + dσᶠdϵᶠ_inv = inv(get_unknowns(stress_state, dσdϵ_full)) # f: unknown strain components solved for during stress iterations + ssd.dϵdp[sc, :] .= .-dσᶠdϵᶠ_inv * ssd.mderiv.dσdp[sc, :] + ssd.dσdp[ec, :] .= ssd.mderiv.dσdp .+ ssd.mderiv.dσdϵ[ec, sc] * ssd.dϵdp[sc, :] + ssd.mderiv.dsdp .+= ssd.mderiv.dsdϵ[:, sc] * ssd.dϵdp[sc, :] + return return reduce_tensordim(stress_state, σ_full), reduce_stiffness(stress_state, dσdϵ_full), state, ϵ_full +end + +""" + stress_controlled_indices(stress_state::AbstractStressState, ::AbstractTensor)::SVector{N, Int} + +Get the `N` indices that are stress-controlled in `stress_state`. The tensor input is used to +determine if a symmetric or full tensor is used. +""" +function stress_controlled_indices end + +""" + strain_controlled_indices(stress_state::AbstractStressState, ::AbstractTensor)::SVector{N, Int} + +Get the `N` indices that are strain-controlled in `stress_state`. The tensor input is used to +determine if a symmetric or full tensor is used. +""" +function strain_controlled_indices end + +# NoIterationState +stress_controlled_indices(::NoIterationState, ::AbstractTensor) = SVector{0,Int}() +strain_controlled_indices(::NoIterationState, ::SymmetricTensor) = @SVector([1, 2, 3, 4, 5, 6]) +strain_controlled_indices(::NoIterationState, ::Tensor) = @SVector([1, 2, 3, 4, 5, 6, 7, 8, 9]) + +# UniaxialStress +stress_controlled_indices(::UniaxialStress, ::SymmetricTensor) = @SVector([2, 3, 4, 5, 6]) +stress_controlled_indices(::UniaxialStress, ::Tensor) = @SVector([2, 3, 4, 5, 6, 7, 8, 9]) +strain_controlled_indices(::UniaxialStress, ::AbstractTensor) = @SVector([1]) + +# UniaxialNormalStress +stress_controlled_indices(::UniaxialNormalStress, ::AbstractTensor) = @SVector([2,3]) +strain_controlled_indices(::UniaxialNormalStress, ::SymmetricTensor) = @SVector([1, 4, 5, 6]) +strain_controlled_indices(::UniaxialNormalStress, ::Tensor) = @SVector([1, 4, 5, 6, 7, 8, 9]) + +# PlaneStress 12 -> 6, 21 -> 9 +stress_controlled_indices(::PlaneStress, ::SymmetricTensor) = @SVector([3, 4, 5]) +stress_controlled_indices(::PlaneStress, ::Tensor) = @SVector([3, 4, 5, 7, 8]) +strain_controlled_indices(::PlaneStress, ::SymmetricTensor) = @SVector([1, 2, 6]) +strain_controlled_indices(::PlaneStress, ::Tensor) = @SVector([1, 2, 6, 9]) + +# GeneralStressState +stress_controlled_indices(ss::GeneralStressState{Nσ}, ::AbstractTensor) where Nσ = controlled_indices_from_tensor(ss.σ_ctrl, true, Val(Nσ)) +function strain_controlled_indices(ss::GeneralStressState{Nσ,TT}, ::AbstractTensor) where {Nσ,TT} + N = Tensors.n_components(Tensors.get_base(TT)) - Nσ + return controlled_indices_from_tensor(ss.σ_ctrl, false, Val(N)) +end +function controlled_indices_from_tensor(ctrl::AbstractTensor, return_if, ::Val{N}) where N + return SVector{N}(i for (i, v) in pairs(tovoigt(SVector, ctrl)) if v == return_if) +end \ No newline at end of file diff --git a/src/stressiterations.jl b/src/stressiterations.jl index aaff31e..8de71e5 100644 --- a/src/stressiterations.jl +++ b/src/stressiterations.jl @@ -25,7 +25,8 @@ See also [`ReducedStressState`](@ref). material_response(::AbstractStressState, ::AbstractMaterial, args...) @inline function material_response(stress_state::AbstractStressState, m::AbstractMaterial, args::Vararg{Any,N}) where N - return reduced_material_response(stress_state, m, args...) + σ, dσdϵ, state, ϵ_full = stress_state_material_response(stress_state, m, args...) + return reduce_tensordim(stress_state, σ), reduce_stiffness(stress_state, dσdϵ), state, ϵ_full end update_stress_state!(::AbstractStressState, σ) = nothing @@ -168,7 +169,7 @@ get_tolerance(ss::UniaxialNormalStress) = get_tolerance(ss.settings) get_max_iter(ss::UniaxialNormalStress) = get_max_iter(ss.settings) """ - GeneralStressState(σ_ctrl::AbstractTensor{2,3,Bool}, σ::AbstractTensor{2,3,Bool}; kwargs...) + GeneralStressState(σ_ctrl::AbstractTensor{2,3,Bool}, σ::AbstractTensor{2,3}; kwargs...) Construct a general stress state controlled by `σ_ctrl` whose component is `true` if that component is stress-controlled and `false` if it is strain-controlled. If stress-controlled, @@ -252,15 +253,15 @@ reduce_tensordim(::Val{dim}, a::SymmetricTensor{2}) where dim = SymmetricTensor{ reduce_tensordim(::Val{dim}, A::SymmetricTensor{4}) where dim = SymmetricTensor{4,dim}((i,j,k,l)->A[i,j,k,l]) -function reduced_material_response(stress_state::NoIterationState, +function stress_state_material_response(stress_state::NoIterationState, m::AbstractMaterial, ϵ::AbstractTensor, args::Vararg{Any,N}) where N ϵ_full = expand_tensordim(stress_state, ϵ) σ, dσdϵ, new_state = material_response(m, ϵ_full, args...) - return reduce_tensordim(stress_state, σ), reduce_tensordim(stress_state, dσdϵ), new_state, ϵ_full + return σ, dσdϵ, new_state, ϵ_full end -function reduced_material_response(stress_state::IterationState, +function stress_state_material_response(stress_state::IterationState, m::AbstractMaterial, ϵ::AbstractTensor, args::Vararg{Any,N}) where N # Newton options @@ -274,8 +275,7 @@ function reduced_material_response(stress_state::IterationState, σ_mandel = get_unknowns(stress_state, σ_full) if norm(σ_mandel) < tol - dσdϵ = reduce_stiffness(stress_state, dσdϵ_full) - return reduce_tensordim(stress_state, σ_full), dσdϵ, new_state, ϵ_full + return σ_full, dσdϵ_full, new_state, ϵ_full end dσdϵ_mandel = get_unknowns(stress_state, dσdϵ_full) @@ -284,14 +284,22 @@ function reduced_material_response(stress_state::IterationState, throw(NoStressConvergence("Stress iterations with the NewtonSolver did not converge")) end -reduce_stiffness(::State3D, dσdϵ_full::AbstractTensor{4,3}) = dσdϵ_full +reduce_stiffness(ss::State3D, dσdϵ_full) = reduce_stiffness(Val(3), ss, dσdϵ_full) +reduce_stiffness(ss::State2D, dσdϵ_full) = reduce_stiffness(Val(2), ss, dσdϵ_full) +reduce_stiffness(ss::State1D, dσdϵ_full) = reduce_stiffness(Val(1), ss, dσdϵ_full) -function reduce_stiffness(stress_state, dσdϵ_full::AbstractTensor{4,3}) +reduce_stiffness(::Val{3}, ::State3D, dσdϵ_full::AbstractTensor{4,3}) = dσdϵ_full + +function reduce_stiffness(::Union{Val{1}, Val{2}}, stress_state::IterationState, dσdϵ_full::AbstractTensor{4,3}) ∂σᶠ∂ϵᶠ, ∂σᶠ∂ϵᶜ, ∂σᶜ∂ϵᶠ, ∂σᶜ∂ϵᶜ = extract_substiffnesses(stress_state, dσdϵ_full) dσᶜdϵᶜ = ∂σᶜ∂ϵᶜ - ∂σᶜ∂ϵᶠ * (∂σᶠ∂ϵᶠ \ ∂σᶠ∂ϵᶜ) return convert_stiffness(dσᶜdϵᶜ, stress_state, dσdϵ_full) end +function reduce_stiffness(::Union{Val{1}, Val{2}}, stress_state::NoIterationState, dσdϵ_full::AbstractTensor{4,3}) + return reduce_tensordim(stress_state, dσdϵ_full) +end + convert_stiffness(dσᶜdϵᶜ::SMatrix{1,1}, ::State1D, ::SymmetricTensor) = frommandel(SymmetricTensor{4,1}, dσᶜdϵᶜ) convert_stiffness(dσᶜdϵᶜ::SMatrix{1,1}, ::State1D, ::Tensor) = frommandel(Tensor{4,1}, dσᶜdϵᶜ) convert_stiffness(dσᶜdϵᶜ::SMatrix{3,3}, ::State2D, ::SymmetricTensor) = frommandel(SymmetricTensor{4,2}, dσᶜdϵᶜ) @@ -411,7 +419,7 @@ function get_full_tensor(state::GeneralStressState{Nσ}, ::TT, v::SVector{Nσ,T} end function get_unknowns(state::GeneralStressState{Nσ}, a::AbstractTensor{2,3,T}) where {Nσ, T} - shear_factor = T(√2) + shear_factor = a isa SymmetricTensor ? T(√2) : one(T) s(i,j) = i==j ? one(T) : shear_factor f(c) = ((i,j) = c; s(i,j)*(a[i,j]-state.σ[i,j])) return SVector{Nσ}((f(c) for c in state.σm_inds)) diff --git a/src/vector_conversion.jl b/src/vector_conversion.jl index 8d6cf7f..55c0b55 100644 --- a/src/vector_conversion.jl +++ b/src/vector_conversion.jl @@ -83,6 +83,8 @@ Puts the Mandel components of `a` into the vector `v`. """ function tovector! end +tovector!(v::AbstractVector, ::NoMaterialState; kwargs...) = v + # Tensors.jl implementation function tovector!(v::AbstractVector, a::SecondOrderTensor; offset = 0) return tomandel!(v, a; offset) @@ -107,6 +109,8 @@ Create a tensor with shape of `TT` with entries from the Mandel components in `v """ function fromvector end +fromvector(::AbstractVector, ::NoMaterialState; kwargs...) = NoMaterialState() + # Tensors.jl implementation function fromvector(v::AbstractVector, ::TT; offset = 0) where {TT <: SecondOrderTensor} return frommandel(Tensors.get_base(TT), v; offset) diff --git a/test/TestMaterials/Manifest.toml b/test/TestMaterials/Manifest.toml new file mode 100644 index 0000000..0cb032c --- /dev/null +++ b/test/TestMaterials/Manifest.toml @@ -0,0 +1,216 @@ +# This file is machine-generated - editing it directly is not advised + +julia_version = "1.11.6" +manifest_format = "2.0" +project_hash = "a195ace9399d41c37efb6635bf312fbd3e9045d6" + +[[deps.Artifacts]] +uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" +version = "1.11.0" + +[[deps.CommonSubexpressions]] +deps = ["MacroTools"] +git-tree-sha1 = "cda2cfaebb4be89c9084adaca7dd7333369715c5" +uuid = "bbf7d656-a473-5ed7-a52c-81e309532950" +version = "0.3.1" + +[[deps.CompilerSupportLibraries_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" +version = "1.1.1+0" + +[[deps.Dates]] +deps = ["Printf"] +uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" +version = "1.11.0" + +[[deps.DiffResults]] +deps = ["StaticArraysCore"] +git-tree-sha1 = "782dd5f4561f5d267313f23853baaaa4c52ea621" +uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" +version = "1.1.0" + +[[deps.DiffRules]] +deps = ["IrrationalConstants", "LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"] +git-tree-sha1 = "23163d55f885173722d1e4cf0f6110cdbaf7e272" +uuid = "b552c78f-8df3-52c6-915a-8e097449b14b" +version = "1.15.1" + +[[deps.DocStringExtensions]] +git-tree-sha1 = "7442a5dfe1ebb773c29cc2962a8980f47221d76c" +uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +version = "0.9.5" + +[[deps.ForwardDiff]] +deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions"] +git-tree-sha1 = "ce15956960057e9ff7f1f535400ffa14c92429a4" +uuid = "f6369f11-7733-5829-9624-2563aa707210" +version = "1.1.0" +weakdeps = ["StaticArrays"] + + [deps.ForwardDiff.extensions] + ForwardDiffStaticArraysExt = "StaticArrays" + +[[deps.IrrationalConstants]] +git-tree-sha1 = "e2222959fbc6c19554dc15174c81bf7bf3aa691c" +uuid = "92d709cd-6900-40b7-9082-c6be49f344b6" +version = "0.2.4" + +[[deps.JLLWrappers]] +deps = ["Artifacts", "Preferences"] +git-tree-sha1 = "0533e564aae234aff59ab625543145446d8b6ec2" +uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" +version = "1.7.1" + +[[deps.Libdl]] +uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" +version = "1.11.0" + +[[deps.LinearAlgebra]] +deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"] +uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +version = "1.11.0" + +[[deps.LogExpFunctions]] +deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"] +git-tree-sha1 = "13ca9e2586b89836fd20cccf56e57e2b9ae7f38f" +uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688" +version = "0.3.29" + + [deps.LogExpFunctions.extensions] + LogExpFunctionsChainRulesCoreExt = "ChainRulesCore" + LogExpFunctionsChangesOfVariablesExt = "ChangesOfVariables" + LogExpFunctionsInverseFunctionsExt = "InverseFunctions" + + [deps.LogExpFunctions.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + ChangesOfVariables = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0" + InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" + +[[deps.MacroTools]] +git-tree-sha1 = "1e0228a030642014fe5cfe68c2c0a818f9e3f522" +uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" +version = "0.5.16" + +[[deps.MaterialModelsBase]] +deps = ["ForwardDiff", "StaticArrays", "Tensors"] +path = "../.." +uuid = "af893363-701d-44dc-8b1e-d9a2c129bfc9" +version = "0.2.3" + +[[deps.NaNMath]] +deps = ["OpenLibm_jll"] +git-tree-sha1 = "9b8215b1ee9e78a293f99797cd31375471b2bcae" +uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" +version = "1.1.3" + +[[deps.OpenBLAS_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] +uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" +version = "0.3.27+1" + +[[deps.OpenLibm_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "05823500-19ac-5b8b-9628-191a04bc5112" +version = "0.8.5+0" + +[[deps.OpenSpecFun_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl"] +git-tree-sha1 = "1346c9208249809840c91b26703912dff463d335" +uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" +version = "0.5.6+0" + +[[deps.PrecompileTools]] +deps = ["Preferences"] +git-tree-sha1 = "5aa36f7049a63a1528fe8f7c3f2113413ffd4e1f" +uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a" +version = "1.2.1" + +[[deps.Preferences]] +deps = ["TOML"] +git-tree-sha1 = "0f27480397253da18fe2c12a4ba4eb9eb208bf3d" +uuid = "21216c6a-2e73-6563-6e65-726566657250" +version = "1.5.0" + +[[deps.Printf]] +deps = ["Unicode"] +uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" +version = "1.11.0" + +[[deps.Random]] +deps = ["SHA"] +uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +version = "1.11.0" + +[[deps.SHA]] +uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" +version = "0.7.0" + +[[deps.SIMD]] +deps = ["PrecompileTools"] +git-tree-sha1 = "fea870727142270bdf7624ad675901a1ee3b4c87" +uuid = "fdea26ae-647d-5447-a871-4b548cad5224" +version = "3.7.1" + +[[deps.SpecialFunctions]] +deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"] +git-tree-sha1 = "41852b8679f78c8d8961eeadc8f62cef861a52e3" +uuid = "276daf66-3868-5448-9aa4-cd146d93841b" +version = "2.5.1" + + [deps.SpecialFunctions.extensions] + SpecialFunctionsChainRulesCoreExt = "ChainRulesCore" + + [deps.SpecialFunctions.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + +[[deps.StaticArrays]] +deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"] +git-tree-sha1 = "cbea8a6bd7bed51b1619658dec70035e07b8502f" +uuid = "90137ffa-7385-5640-81b9-e52037218182" +version = "1.9.14" + + [deps.StaticArrays.extensions] + StaticArraysChainRulesCoreExt = "ChainRulesCore" + StaticArraysStatisticsExt = "Statistics" + + [deps.StaticArrays.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" + +[[deps.StaticArraysCore]] +git-tree-sha1 = "192954ef1208c7019899fbf8049e717f92959682" +uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" +version = "1.4.3" + +[[deps.Statistics]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "ae3bb1eb3bba077cd276bc5cfc337cc65c3075c0" +uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +version = "1.11.1" + + [deps.Statistics.extensions] + SparseArraysExt = ["SparseArrays"] + + [deps.Statistics.weakdeps] + SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" + +[[deps.TOML]] +deps = ["Dates"] +uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" +version = "1.0.3" + +[[deps.Tensors]] +deps = ["ForwardDiff", "LinearAlgebra", "PrecompileTools", "SIMD", "StaticArrays", "Statistics"] +git-tree-sha1 = "399ec4c46786c47380121f8a6497c65b89f0573f" +uuid = "48a634ad-e948-5137-8d70-aa71f2a747f4" +version = "1.16.2" + +[[deps.Unicode]] +uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" +version = "1.11.0" + +[[deps.libblastrampoline_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" +version = "5.11.0+0" diff --git a/test/TestMaterials/Project.toml b/test/TestMaterials/Project.toml new file mode 100644 index 0000000..9d2e072 --- /dev/null +++ b/test/TestMaterials/Project.toml @@ -0,0 +1,15 @@ +name = "TestMaterials" +uuid = "aa4f4413-326f-4ddd-a53e-c58801ebfeaf" +authors = ["Knut Andreas "] +version = "0.1.0" + +[deps] +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +MaterialModelsBase = "af893363-701d-44dc-8b1e-d9a2c129bfc9" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" +Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4" + +[compat] +ForwardDiff = "1.1.0" +StaticArrays = "1.9.14" +Tensors = "1.16.2" diff --git a/test/TestMaterials.jl b/test/TestMaterials/src/TestMaterials.jl similarity index 95% rename from test/TestMaterials.jl rename to test/TestMaterials/src/TestMaterials.jl index e90147e..85cf857 100644 --- a/test/TestMaterials.jl +++ b/test/TestMaterials/src/TestMaterials.jl @@ -5,9 +5,10 @@ using Tensors, StaticArrays, ForwardDiff # Overloaded functions import MaterialModelsBase as MMB -# Exported materials -export LinearElastic, ViscoElastic -# Todo: Also NeoHookean to test nonlinear geometric materials +# Test material implementations +# * LinearElastic +# * ViscoElastic +# TODO: NeoHookean (geometrically nonlinear materials) # LinearElastic material struct LinearElastic{T} <: AbstractMaterial diff --git a/test/differentiation.jl b/test/differentiation.jl index 70ba2f6..75d7bae 100644 --- a/test/differentiation.jl +++ b/test/differentiation.jl @@ -1,3 +1,24 @@ +@testset "differentiation_utils" begin + for TT in (Tensor, SymmetricTensor) + ϵmock = TT == Tensor ? ones(TT{2,3}) : TT{2,3}((i,j) -> i==j ? 1.0 : 1/√2) + for stress_state in ( + UniaxialStress(), UniaxialNormalStress(), PlaneStress(), + GeneralStressState(rand(TT{2, 3, Bool}), zero(TT{2, 3})), + ) + sc = MMB.stress_controlled_indices(stress_state, ϵmock) + ec = MMB.strain_controlled_indices(stress_state, ϵmock) + @test length(union(Set(sc), Set(ec))) == Tensors.n_components(TT{2, 3}) + @test tomandel(ϵmock)[sc] ≈ MMB.get_unknowns(stress_state, ϵmock) + end + for stress_state in (UniaxialStrain(), PlaneStrain(), FullStressState()) + sc = MMB.stress_controlled_indices(stress_state, ϵmock) + ec = MMB.strain_controlled_indices(stress_state, ϵmock) + @test length(sc) == 0 + @test length(ec) == Tensors.n_components(TT{2,3}) + end + end +end + @testset "differentiation" begin G = rand() K = G + rand() diff --git a/test/runtests.jl b/test/runtests.jl index ca289c3..823f7bf 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,11 +1,9 @@ using MaterialModelsBase using Test using Tensors, StaticArrays +import MaterialModelsBase as MMB +using TestMaterials: LinearElastic, ViscoElastic -const MMB = MaterialModelsBase - -include("TestMaterials.jl") -using .TestMaterials include("utils4testing.jl") include("vector_conversion.jl") From d6e221c0bb92a8f0a1d756b5958b509b641d1bea Mon Sep 17 00:00:00 2001 From: Knut Andreas Date: Wed, 27 Aug 2025 12:53:18 +0200 Subject: [PATCH 02/14] Fix get_unknowns for GeneralStressState with no unknowns --- src/stressiterations.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/stressiterations.jl b/src/stressiterations.jl index 8de71e5..932b43a 100644 --- a/src/stressiterations.jl +++ b/src/stressiterations.jl @@ -422,14 +422,14 @@ function get_unknowns(state::GeneralStressState{Nσ}, a::AbstractTensor{2,3,T}) shear_factor = a isa SymmetricTensor ? T(√2) : one(T) s(i,j) = i==j ? one(T) : shear_factor f(c) = ((i,j) = c; s(i,j)*(a[i,j]-state.σ[i,j])) - return SVector{Nσ}((f(c) for c in state.σm_inds)) + return SVector{Nσ,T}((f(c) for c in state.σm_inds)) end function get_unknowns(state::GeneralStressState{Nσ}, a::AbstractTensor{4,3,T}) where {Nσ,T} shear_factor = T(√2) s(i,j) = i==j ? one(T) : shear_factor f(c1,c2) = ((i,j) = c1; (k,l) = c2; a[i,j,k,l]*s(i,j)*s(k,l)) - return SMatrix{Nσ,Nσ}((f(c1,c2) for c1 in state.σm_inds, c2 in state.σm_inds)) + return SMatrix{Nσ,Nσ,T}((f(c1,c2) for c1 in state.σm_inds, c2 in state.σm_inds)) end # Extract each part of the stiffness tensor as SMatrix From c2ed48bdf3bbd5dab9f1bafde1835e1d65e558b6 Mon Sep 17 00:00:00 2001 From: Knut Andreas Date: Wed, 27 Aug 2025 12:54:38 +0200 Subject: [PATCH 03/14] Set min julia to 1.11 to allow using sources --- .github/workflows/CI.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 1d306e1..f7dd871 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -19,7 +19,7 @@ jobs: matrix: version: - '1' - - '1.8' + - '1.11' os: - ubuntu-latest - macOS-latest From 04e775b663df644204007d860ff4b59f842282f4 Mon Sep 17 00:00:00 2001 From: Knut Andreas Date: Wed, 27 Aug 2025 12:54:50 +0200 Subject: [PATCH 04/14] Set min julia to 1.11 to allow using sources (also in Project.toml) --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 146a130..e0bf6e3 100644 --- a/Project.toml +++ b/Project.toml @@ -9,7 +9,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4" [compat] -julia = "1.8" +julia = "1.11" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" From 7ae7b06262bb9531e381feb87da062e30415885c Mon Sep 17 00:00:00 2001 From: Knut Andreas Date: Wed, 27 Aug 2025 13:41:40 +0200 Subject: [PATCH 05/14] Add devdocs --- docs/make.jl | 1 + docs/src/devdocs.md | 8 ++++++++ 2 files changed, 9 insertions(+) create mode 100644 docs/src/devdocs.md diff --git a/docs/make.jl b/docs/make.jl index 443ef17..ca062e1 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -19,6 +19,7 @@ makedocs(; "Conversion" => "conversion.md", "Differentiation" => "differentiation.md", "Implementation" => "implementing.md", + "Devdocs" => "devdocs.md" ], ) diff --git a/docs/src/devdocs.md b/docs/src/devdocs.md new file mode 100644 index 0000000..db1ab35 --- /dev/null +++ b/docs/src/devdocs.md @@ -0,0 +1,8 @@ +# Internal documentation +The following documentation returns to internals of the package and may change +at any time. + +```@docs +MaterialModelsBase.stress_controlled_indices +MaterialModelsBase.strain_controlled_indices +``` From 31ecae874b8299d75983fa57debc364bba7dae75 Mon Sep 17 00:00:00 2001 From: Knut Andreas Date: Sat, 30 Aug 2025 22:51:27 +0200 Subject: [PATCH 06/14] Add proper tests for differentiation --- Project.toml | 3 +- src/MaterialModelsBase.jl | 3 +- src/differentiation.jl | 54 +++++++++--- test/TestMaterials/Manifest.toml | 103 +++++++++++++++++++++- test/TestMaterials/Project.toml | 2 + test/TestMaterials/src/TestMaterials.jl | 48 ++++++++++- test/TestMaterials/src/loadcases.jl | 85 ++++++++++++++++++ test/TestMaterials/src/numdiff.jl | 90 +++++++++++++++++++ test/differentiation.jl | 109 ++++++++++++++++++++---- test/runtests.jl | 3 +- 10 files changed, 466 insertions(+), 34 deletions(-) create mode 100644 test/TestMaterials/src/loadcases.jl create mode 100644 test/TestMaterials/src/numdiff.jl diff --git a/Project.toml b/Project.toml index e0bf6e3..338a28e 100644 --- a/Project.toml +++ b/Project.toml @@ -14,9 +14,10 @@ julia = "1.11" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TestMaterials = "aa4f4413-326f-4ddd-a53e-c58801ebfeaf" +FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" [sources] TestMaterials = {path = "test/TestMaterials"} [targets] -test = ["Test", "TestMaterials"] +test = ["Test", "FiniteDiff", "TestMaterials"] diff --git a/src/MaterialModelsBase.jl b/src/MaterialModelsBase.jl index 0a01bc7..03a10e8 100644 --- a/src/MaterialModelsBase.jl +++ b/src/MaterialModelsBase.jl @@ -26,7 +26,8 @@ export update_stress_state! # For nonzero stress export tovector, tovector!, fromvector # Convert to/from `AbstractVector`s export get_num_tensorcomponents, get_num_statevars # Information about the specific material export get_num_params, get_params_eltype # -export MaterialDerivatives, differentiate_material! # Differentiation routines +export MaterialDerivatives, StressStateDerivatives # Derivative collections +export differentiate_material! # Differentiation routines export allocate_differentiation_output # abstract type AbstractMaterial end diff --git a/src/differentiation.jl b/src/differentiation.jl index e4702e9..69b030d 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -1,12 +1,20 @@ struct MaterialDerivatives{T} dσdϵ::Matrix{T} - dσdⁿs::Matrix{T} + #dσdⁿs::Matrix{T} dσdp::Matrix{T} dsdϵ::Matrix{T} - dsdⁿs::Matrix{T} + #dsdⁿs::Matrix{T} dsdp::Matrix{T} end +function Base.getproperty(d::MaterialDerivatives, key::Symbol) + if key === :dσdⁿs || key === :dsdⁿs + error("You are probably assuming MaterialModelsBase v0.2 behavior for differentiation") + else + @inline getfield(d, key) + end +end + """ MaterialDerivatives(m::AbstractMaterial) @@ -36,10 +44,8 @@ function MaterialDerivatives(m::AbstractMaterial) dsdp = ForwardDiff.jacobian(p -> tovector(initial_material_state(fromvector(p, m))), tovector(m)) return MaterialDerivatives( zeros(T, n_tensor, n_tensor), # dσdϵ - zeros(T, n_tensor, n_state), # dσdⁿs zeros(T, n_tensor, n_params), # dσdp zeros(T, n_state, n_tensor), # dsdϵ - zeros(T, n_state, n_state), # dsdⁿs dsdp ) end @@ -82,6 +88,22 @@ struct StressStateDerivatives{T} # TODO: Reduce the dimensions, for now all entries (even those that are zero) are stored. end +function StressStateDerivatives(::AbstractStressState, m::AbstractMaterial) + mderiv = MaterialDerivatives(m) + np = get_num_params(m) + nt = get_num_tensorcomponents(m) # Should be changed to only save non-controlled entries + dϵdp = zeros(nt, np) + dσdp = zeros(nt, np) + # Should be changed to only save non-controlled entries + vo = Tensors.DEFAULT_VOIGT_ORDER[3] + if nt == 6 + index = SMatrix{3, 3, Int}(min(vo[i, j], vo[j, i]) for i in 1:3, j in 1:3) + else + index = SMatrix{3, 3, Int}(vo) + end + return StressStateDerivatives(mderiv, dϵdp, dσdp, index, index) +end + """ differentiate_material!(ssd::StressStateDerivatives, stress_state, m, args...) @@ -95,14 +117,22 @@ For material models that directly implement `material_response(stress_state, m, to calculate the derivatives in `ssd`. Here the user has full control and no modifications occur automatically, however, typically the (total) derivatives `ssd.dσdp`, `ssd.dϵdp`, and `ssd.mderiv.dsdp` should be updated. """ -function differentiate_material!(ssd::StressStateDerivatives, stress_state::AbstractStressState, m::AbstractMaterial, args::Vararg{Any,N}) where {N} - σ_full, dσdϵ_full, state, ϵ_full = stress_state_material_response(stress_state, m, args...) - differentiate_material!(ssd.mderiv, m, args..., dσdϵ_full) - dσᶠdϵᶠ_inv = inv(get_unknowns(stress_state, dσdϵ_full)) # f: unknown strain components solved for during stress iterations - ssd.dϵdp[sc, :] .= .-dσᶠdϵᶠ_inv * ssd.mderiv.dσdp[sc, :] - ssd.dσdp[ec, :] .= ssd.mderiv.dσdp .+ ssd.mderiv.dσdϵ[ec, sc] * ssd.dϵdp[sc, :] - ssd.mderiv.dsdp .+= ssd.mderiv.dsdϵ[:, sc] * ssd.dϵdp[sc, :] - return return reduce_tensordim(stress_state, σ_full), reduce_stiffness(stress_state, dσdϵ_full), state, ϵ_full +function differentiate_material!(ssd::StressStateDerivatives, stress_state::AbstractStressState, m::AbstractMaterial, ϵ::AbstractTensor, args::Vararg{Any,N}) where {N} + σ_full, dσdϵ_full, state, ϵ_full = stress_state_material_response(stress_state, m, ϵ, args...) + differentiate_material!(ssd.mderiv, m, ϵ_full, args..., dσdϵ_full) + + if isa(stress_state, NoIterationState) + copy!(ssd.dσdp, ssd.mderiv.dσdp) + fill!(ssd.dϵdp, 0) + else + sc = stress_controlled_indices(stress_state, ϵ) + ec = strain_controlled_indices(stress_state, ϵ) + dσᶠdϵᶠ_inv = inv(get_unknowns(stress_state, dσdϵ_full)) # f: unknown strain components solved for during stress iterations + ssd.dϵdp[sc, :] .= .-dσᶠdϵᶠ_inv * ssd.mderiv.dσdp[sc, :] + ssd.dσdp[ec, :] .= ssd.mderiv.dσdp[ec, :] .+ ssd.mderiv.dσdϵ[ec, sc] * ssd.dϵdp[sc, :] + ssd.mderiv.dsdp .+= ssd.mderiv.dsdϵ[:, sc] * ssd.dϵdp[sc, :] + end + return reduce_tensordim(stress_state, σ_full), reduce_stiffness(stress_state, dσdϵ_full), state, ϵ_full end """ diff --git a/test/TestMaterials/Manifest.toml b/test/TestMaterials/Manifest.toml index 0cb032c..328d424 100644 --- a/test/TestMaterials/Manifest.toml +++ b/test/TestMaterials/Manifest.toml @@ -2,7 +2,53 @@ julia_version = "1.11.6" manifest_format = "2.0" -project_hash = "a195ace9399d41c37efb6635bf312fbd3e9045d6" +project_hash = "d41ae232891445fa7506f236787c2504cd12d9df" + +[[deps.Adapt]] +deps = ["LinearAlgebra", "Requires"] +git-tree-sha1 = "f7817e2e585aa6d924fd714df1e2a84be7896c60" +uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" +version = "4.3.0" + + [deps.Adapt.extensions] + AdaptSparseArraysExt = "SparseArrays" + AdaptStaticArraysExt = "StaticArrays" + + [deps.Adapt.weakdeps] + SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" + StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" + +[[deps.ArrayInterface]] +deps = ["Adapt", "LinearAlgebra"] +git-tree-sha1 = "9606d7832795cbef89e06a550475be300364a8aa" +uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" +version = "7.19.0" + + [deps.ArrayInterface.extensions] + ArrayInterfaceBandedMatricesExt = "BandedMatrices" + ArrayInterfaceBlockBandedMatricesExt = "BlockBandedMatrices" + ArrayInterfaceCUDAExt = "CUDA" + ArrayInterfaceCUDSSExt = "CUDSS" + ArrayInterfaceChainRulesCoreExt = "ChainRulesCore" + ArrayInterfaceChainRulesExt = "ChainRules" + ArrayInterfaceGPUArraysCoreExt = "GPUArraysCore" + ArrayInterfaceReverseDiffExt = "ReverseDiff" + ArrayInterfaceSparseArraysExt = "SparseArrays" + ArrayInterfaceStaticArraysCoreExt = "StaticArraysCore" + ArrayInterfaceTrackerExt = "Tracker" + + [deps.ArrayInterface.weakdeps] + BandedMatrices = "aae01518-5342-5314-be14-df237901396f" + BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0" + CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" + CUDSS = "45b445bb-4962-46a0-9369-b4df9d0f772e" + ChainRules = "082447d4-558c-5d27-93f4-14fc19e9eca2" + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527" + ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" + SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" + StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" + Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" [[deps.Artifacts]] uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" @@ -19,6 +65,21 @@ deps = ["Artifacts", "Libdl"] uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" version = "1.1.1+0" +[[deps.ConstructionBase]] +git-tree-sha1 = "b4b092499347b18a015186eae3042f72267106cb" +uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9" +version = "1.6.0" + + [deps.ConstructionBase.extensions] + ConstructionBaseIntervalSetsExt = "IntervalSets" + ConstructionBaseLinearAlgebraExt = "LinearAlgebra" + ConstructionBaseStaticArraysExt = "StaticArrays" + + [deps.ConstructionBase.weakdeps] + IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953" + LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" + StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" + [[deps.Dates]] deps = ["Printf"] uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" @@ -41,6 +102,24 @@ git-tree-sha1 = "7442a5dfe1ebb773c29cc2962a8980f47221d76c" uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" version = "0.9.5" +[[deps.FiniteDiff]] +deps = ["ArrayInterface", "LinearAlgebra", "Setfield"] +git-tree-sha1 = "31fd32af86234b6b71add76229d53129aa1b87a9" +uuid = "6a86dc24-6348-571c-b903-95158fe2bd41" +version = "2.28.1" + + [deps.FiniteDiff.extensions] + FiniteDiffBandedMatricesExt = "BandedMatrices" + FiniteDiffBlockBandedMatricesExt = "BlockBandedMatrices" + FiniteDiffSparseArraysExt = "SparseArrays" + FiniteDiffStaticArraysExt = "StaticArrays" + + [deps.FiniteDiff.weakdeps] + BandedMatrices = "aae01518-5342-5314-be14-df237901396f" + BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0" + SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" + StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" + [[deps.ForwardDiff]] deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions"] git-tree-sha1 = "ce15956960057e9ff7f1f535400ffa14c92429a4" @@ -51,6 +130,11 @@ weakdeps = ["StaticArrays"] [deps.ForwardDiff.extensions] ForwardDiffStaticArraysExt = "StaticArrays" +[[deps.Future]] +deps = ["Random"] +uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820" +version = "1.11.0" + [[deps.IrrationalConstants]] git-tree-sha1 = "e2222959fbc6c19554dc15174c81bf7bf3aa691c" uuid = "92d709cd-6900-40b7-9082-c6be49f344b6" @@ -142,6 +226,12 @@ deps = ["SHA"] uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" version = "1.11.0" +[[deps.Requires]] +deps = ["UUIDs"] +git-tree-sha1 = "62389eeff14780bfe55195b7204c0d8738436d64" +uuid = "ae029012-a4dd-5104-9daa-d747884805df" +version = "1.3.1" + [[deps.SHA]] uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" version = "0.7.0" @@ -152,6 +242,12 @@ git-tree-sha1 = "fea870727142270bdf7624ad675901a1ee3b4c87" uuid = "fdea26ae-647d-5447-a871-4b548cad5224" version = "3.7.1" +[[deps.Setfield]] +deps = ["ConstructionBase", "Future", "MacroTools", "StaticArraysCore"] +git-tree-sha1 = "c5391c6ace3bc430ca630251d02ea9687169ca68" +uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46" +version = "1.1.2" + [[deps.SpecialFunctions]] deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"] git-tree-sha1 = "41852b8679f78c8d8961eeadc8f62cef861a52e3" @@ -206,6 +302,11 @@ git-tree-sha1 = "399ec4c46786c47380121f8a6497c65b89f0573f" uuid = "48a634ad-e948-5137-8d70-aa71f2a747f4" version = "1.16.2" +[[deps.UUIDs]] +deps = ["Random", "SHA"] +uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" +version = "1.11.0" + [[deps.Unicode]] uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" version = "1.11.0" diff --git a/test/TestMaterials/Project.toml b/test/TestMaterials/Project.toml index 9d2e072..956dc43 100644 --- a/test/TestMaterials/Project.toml +++ b/test/TestMaterials/Project.toml @@ -4,12 +4,14 @@ authors = ["Knut Andreas "] version = "0.1.0" [deps] +FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" MaterialModelsBase = "af893363-701d-44dc-8b1e-d9a2c129bfc9" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4" [compat] +FiniteDiff = "2.28.1" ForwardDiff = "1.1.0" StaticArrays = "1.9.14" Tensors = "1.16.2" diff --git a/test/TestMaterials/src/TestMaterials.jl b/test/TestMaterials/src/TestMaterials.jl index 85cf857..793fe79 100644 --- a/test/TestMaterials/src/TestMaterials.jl +++ b/test/TestMaterials/src/TestMaterials.jl @@ -1,10 +1,13 @@ module TestMaterials using MaterialModelsBase -using Tensors, StaticArrays, ForwardDiff +using Tensors, StaticArrays, ForwardDiff, FiniteDiff # Overloaded functions import MaterialModelsBase as MMB +include("numdiff.jl") +include("loadcases.jl") + # Test material implementations # * LinearElastic # * ViscoElastic @@ -57,11 +60,37 @@ struct ViscoElastic{T} <: AbstractMaterial E2::LinearElastic{T} η::T end + +function MMB.get_num_params(m::ViscoElastic) + return sum(get_num_params.((m.E1, m.E2))) + 1 +end + +function MMB.tovector!(v::AbstractVector, m::ViscoElastic; offset = 0) + i = offset + tovector!(v, m.E1; offset = i); i += get_num_params(m.E1) + tovector!(v, m.E2; offset = i); i += get_num_params(m.E2) + v[i + 1] = m.η + return v +end + +function MMB.fromvector(v::AbstractVector, m::ViscoElastic; offset = 0) + i = offset + E1 = fromvector(v, m.E1; offset = i); i += get_num_params(m.E1) + E2 = fromvector(v, m.E2; offset = i); i += get_num_params(m.E2) + η = v[i + 1] + return ViscoElastic(E1, E2, η) +end + struct ViscoElasticState{T} <: AbstractMaterialState ϵv::SymmetricTensor{2,3,T,6} end MMB.initial_material_state(::ViscoElastic) = ViscoElasticState(zero(SymmetricTensor{2,3})) +MMB.get_num_statevars(s::ViscoElasticState) = 6 +MMB.get_statevar_eltype(::ViscoElasticState{T}) where {T} = T +MMB.tovector!(v, s::ViscoElasticState; offset = 0) = tovector!(v, s.ϵv; offset) +MMB.fromvector(v, s::ViscoElasticState; offset = 0) = ViscoElasticState(fromvector(v, s.ϵv; offset)) + # To allow easy automatic differentiation function get_stress_ϵv(m::ViscoElastic, ϵ::SymmetricTensor{2,3}, old::ViscoElasticState, Δt) E1 = get_stiffness(m.E1) @@ -78,7 +107,7 @@ function MMB.material_response( if norm(ϵ) ≈ 0 return zero(ϵ), NaN*zero(get_stiffness(m.E1)), old else - throw(NoLocalConvergence("Δt=0 is not allowed for non-zero strain")) + throw(NoLocalConvergence("Δt=0 is not allowed")) end end σ, ϵv = get_stress_ϵv(m,ϵ,old,Δt) @@ -86,4 +115,19 @@ function MMB.material_response( return σ, dσdϵ, ViscoElasticState(ϵv) end +function MMB.differentiate_material!(diff::MMB.MaterialDerivatives, m::ViscoElastic, ϵ, old::ViscoElasticState, Δt, _, _, dσdϵ) + tomandel!(diff.dσdϵ, dσdϵ) + dσdⁿs = zeros(get_num_tensorcomponents(m), get_num_statevars(m)) + ForwardDiff.jacobian!(dσdⁿs, sv -> tovector(get_stress_ϵv(m, ϵ, fromvector(sv, old), Δt)[1]), tovector(old)) + ∂σ∂p = ForwardDiff.jacobian(p -> tovector(get_stress_ϵv(fromvector(p, m), ϵ, old, Δt)[1]), tovector(m)) + + dsdⁿs = zeros(get_num_statevars(m), get_num_statevars(m)) + ForwardDiff.jacobian!(diff.dsdϵ, e -> tovector(ViscoElasticState(get_stress_ϵv(m, fromvector(e, ϵ), old, Δt)[2])), tovector(ϵ)) + ForwardDiff.jacobian!(dsdⁿs, sv -> tovector(ViscoElasticState(get_stress_ϵv(m, ϵ, fromvector(sv, old), Δt)[2])), tovector(old)) + ∂s∂p = ForwardDiff.jacobian(p -> tovector(ViscoElasticState(get_stress_ϵv(fromvector(p, m), ϵ, old, Δt)[2])), tovector(m)) + + diff.dσdp .= ∂σ∂p .+ dσdⁿs * diff.dsdp + diff.dsdp .= ∂s∂p .+ dsdⁿs * diff.dsdp +end + end \ No newline at end of file diff --git a/test/TestMaterials/src/loadcases.jl b/test/TestMaterials/src/loadcases.jl new file mode 100644 index 0000000..783dfbd --- /dev/null +++ b/test/TestMaterials/src/loadcases.jl @@ -0,0 +1,85 @@ +function symmetric_components(components::NTuple{2, Int}) + return components[1] < components[2] ? reverse(components) : components +end + +construct_tensor(ϵ::AbstractTensor, _) = ϵ + +function construct_tensor(ϵ::Number, components) + i, j = symmetric_components(components) + return SymmetricTensor{2,3}((k, l) -> k == i && l == j ? ϵ : zero(ϵ)) +end + +function runstrain(m, ϵ_ij_end::Number, components, t_end, num_steps) + ϵ_end = construct_tensor(ϵ_ij_end, components) + i, j = symmetric_components(components) + state = initial_material_state(m) + cache = allocate_material_cache(m) + Δt = t_end / num_steps + σv = zeros(num_steps + 1) + for (k, scale) in enumerate(range(0, 1, num_steps + 1)[2:end]) + ϵ = ϵ_end * scale + σ, _, state = material_response(m, ϵ, state, Δt, cache) + σv[k + 1] = σ[i, j] + end + return σv, state +end + +function runstrain_diff(m, ϵ_ij_end::Number, components, t_end, num_steps) + ϵ_end = construct_tensor(ϵ_ij_end, components) + i, j = symmetric_components(components) + mind = Tensors.DEFAULT_VOIGT_ORDER[3][j, i] + mscale = i == j ? 1.0 : 1/√2 + state = initial_material_state(m) + cache = allocate_material_cache(m) + diff = MaterialDerivatives(m) + extras = allocate_differentiation_output(m) + Δt = t_end / num_steps + σv = zeros(num_steps + 1) + dσdp = zeros(num_steps + 1, get_num_params(m)) + for (k, scale) in enumerate(range(0, 1, num_steps + 1)[2:end]) + ϵ = ϵ_end * scale + old_state = state + σ, dσdϵ, state = material_response(m, ϵ, old_state, Δt, cache, extras) + differentiate_material!(diff, m, ϵ, old_state, Δt, cache, extras, dσdϵ) + σv[k + 1] = σ[i, j] + dσdp[k + 1, :] .= diff.dσdp[mind, :] .* mscale + end + return σv, state, dσdp, diff +end + +function runstresstate(stress_state, m, ϵend::Union{Number, AbstractTensor}, components, t_end, num_steps) + ϵt = construct_tensor(ϵend, components) + i, j = symmetric_components(components) + state = initial_material_state(m) + cache = allocate_material_cache(m) + Δt = t_end / num_steps + σv = zeros(num_steps + 1) + for (k, scale) in enumerate(range(0, 1, num_steps + 1)[2:end]) + ϵ = ϵt * scale + σ, _, state = material_response(stress_state, m, ϵ, state, Δt, cache) + σv[k + 1] = σ[i, j] + end + return σv, state +end + +function runstresstate_diff(stress_state, m, ϵend::Union{Number, AbstractTensor}, components, t_end, num_steps) + ϵt = construct_tensor(ϵend, components) + i, j = symmetric_components(components) + mind = Tensors.DEFAULT_VOIGT_ORDER[3][j, i] + mscale = i == j ? 1.0 : 1/√2 + state = initial_material_state(m) + cache = allocate_material_cache(m) + diff = StressStateDerivatives(stress_state, m) + extras = allocate_differentiation_output(m) + Δt = t_end / num_steps + σv = zeros(num_steps + 1) + dσdp = zeros(num_steps + 1, get_num_params(m)) + for (k, scale) in enumerate(range(0, 1, num_steps + 1)[2:end]) + ϵ = ϵt * scale + old_state = state + σ, _, state = differentiate_material!(diff, stress_state, m, ϵ, old_state, Δt, cache, extras) + σv[k + 1] = σ[i, j] + dσdp[k + 1, :] .= diff.dσdp[mind, :] .* mscale + end + return σv, state, dσdp, diff +end diff --git a/test/TestMaterials/src/numdiff.jl b/test/TestMaterials/src/numdiff.jl new file mode 100644 index 0000000..99fa264 --- /dev/null +++ b/test/TestMaterials/src/numdiff.jl @@ -0,0 +1,90 @@ +""" + obtain_numerical_material_derivative!(deriv, m, ϵ, old, Δt; fdtype = Val{:forward}, kwargs...) + +Obtain the numerical derivative of the material `m` at the strain `ϵ` and old state variables, `old`, +for a time step `Δt`. `fdtype` and `kwargs...` are passed to `FiniteDiff.finite_difference_jacobian`. +""" +function obtain_numerical_material_derivative!( + deriv::MMB.MaterialDerivatives, + m::AbstractMaterial, ϵ, old, Δt; + fdtype = Val{:forward}, kwargs... + ) + cache = allocate_material_cache(m) + p = tovector(m) + ⁿs = tovector(old) + e = tovector(ϵ) + + function numjac(f::F, x) where {F} + sz = (length(f(x)), length(x)) + prod(sz) == 0 && return zeros(sz...) + return FiniteDiff.finite_difference_jacobian(f, x, fdtype; kwargs...) + end + numjac!(J, f::F, x) where {F} = copy!(J, numjac(f, x)) + + funs = NamedTuple{(:σ, :s)}( + (s = svec -> tovector(material_response(m, ϵ, fromvector(svec, old), Δt, cache)[i]), # f(s) + p = pvec -> tovector(material_response(fromvector(pvec, m), ϵ, old, Δt, cache)[i]), # f(p) + ϵ = evec -> tovector(material_response(m, fromvector(evec, ϵ), old, Δt, cache)[i])) # f(ϵ) + for i in (1, 3) + ) + + dⁿsdp = copy(deriv.dsdp) # Copy not needed as dsdp field updated last, but safer in this test routine. + + # Fill all partial derivatives + numjac!(deriv.dσdϵ, funs.σ.ϵ, e) + dσdⁿs = numjac(funs.σ.s, ⁿs) + #numjac!(deriv. dσdⁿs, funs.σ.s, ⁿs) # not needed in deriv + numjac!(deriv.dσdp, funs.σ.p, p) + + numjac!(deriv.dsdϵ, funs.s.ϵ, e) + #numjac!(deriv. dsdⁿs, funs.s.s, ⁿs) # not needed in deriv + dsdⁿs = numjac(funs.s.s, ⁿs) + numjac!(deriv.dsdp, funs.s.p, p) + + # Adjust for dependence on previous state variable, i.e. ⁿs = ⁿs(p) + deriv.dσdp .+= dσdⁿs * dⁿsdp + deriv.dsdp .+= dsdⁿs * dⁿsdp + +end + +""" + obtain_numerical_material_derivative!(ssd, stress_state, m, ϵ, old, Δt; fdtype = Val{:forward}, kwargs...) + +Obtain the numerical derivative of the material `m` considering the `stress_state` iterations at the strain +`ϵ` and old state variables, `old`, for a time step `Δt`. `fdtype` and `kwargs...` are passed to `FiniteDiff.finite_difference_jacobian`. +""" +function obtain_numerical_material_derivative!( + ssd::MMB.StressStateDerivatives, + stress_state::AbstractStressState, + m::AbstractMaterial, ϵ, old, Δt; + fdtype = Val{:forward}, kwargs... + ) + cache = allocate_material_cache(m) + p = tovector(m) + ⁿs = tovector(old) + + function numjac(f::F, x) where {F} + sz = (length(f(x)), length(x)) + prod(sz) == 0 && return zeros(sz...) + return FiniteDiff.finite_difference_jacobian(f, x, fdtype; kwargs...) + end + numjac!(J, f::F, x) where {F} = copy!(J, numjac(f, x)) + + funs = NamedTuple{(:σ, :ϵ, :s)}( + (s = svec -> tovector(material_response(stress_state, m, ϵ, fromvector(svec, old), Δt, cache)[i]), # f(s) + p = pvec -> tovector(material_response(stress_state, fromvector(pvec, m), ϵ, old, Δt, cache)[i])) # f(p) + for i in (1, 4, 3) + ) + + dⁿsdp = copy(ssd.mderiv.dsdp) + + #numjac!(ssd.mderiv. dσdⁿs, funs.σ.s, ⁿs) + dσdⁿs = numjac(funs.σ.s, ⁿs) + numjac!(ssd.dσdp, funs.σ.p, p); ssd.dσdp .+= dσdⁿs * dⁿsdp + + dϵdⁿs = numjac(funs.ϵ.s, ⁿs); + numjac!(ssd.dϵdp, funs.ϵ.p, p); ssd.dϵdp .+= dϵdⁿs * dⁿsdp + + dsdⁿs = numjac(funs.s.s, sv) + numjac!(ssd.mderiv.dsdp, funs.s.p, p); ssd.mderiv.dsdp .+= dsdⁿs * dⁿsdp +end diff --git a/test/differentiation.jl b/test/differentiation.jl index 75d7bae..e3af6ad 100644 --- a/test/differentiation.jl +++ b/test/differentiation.jl @@ -20,20 +20,97 @@ end @testset "differentiation" begin - G = rand() - K = G + rand() - m = LinearElastic(G, K) - ϵ = rand(SymmetricTensor{2,3}) - diff = MaterialDerivatives(m) - IxI = one(SymmetricTensor{2,3}) ⊗ one(SymmetricTensor{2,3}) - dσdϵ = 2G*(one(SymmetricTensor{4,3}) - IxI/3) + K*IxI - dσdG = 2*(one(SymmetricTensor{4,3}) - IxI/3) ⊡ ϵ - dσdK = IxI ⊡ ϵ + @testset "Analytical for LinearElastic" begin + G = rand() + K = G + rand() + m = LinearElastic(G, K) + ϵ = rand(SymmetricTensor{2,3}) + diff = MaterialDerivatives(m) + IxI = one(SymmetricTensor{2,3}) ⊗ one(SymmetricTensor{2,3}) + dσdϵ = 2G*(one(SymmetricTensor{4,3}) - IxI/3) + K*IxI + dσdG = 2*(one(SymmetricTensor{4,3}) - IxI/3) ⊡ ϵ + dσdK = IxI ⊡ ϵ - old = NoMaterialState(); Δt = nothing - cache = NoMaterialCache(); extras = allocate_differentiation_output(m) - differentiate_material!(diff, m, ϵ, old, Δt, cache, extras, dσdϵ) - @test diff.dσdϵ ≈ tomandel(dσdϵ) - @test diff.dσdp[:,1] ≈ tomandel(dσdG) - @test diff.dσdp[:,2] ≈ tomandel(dσdK) -end \ No newline at end of file + old = NoMaterialState(); Δt = nothing + cache = NoMaterialCache(); extras = allocate_differentiation_output(m) + differentiate_material!(diff, m, ϵ, old, Δt, cache, extras, dσdϵ) + @test diff.dσdϵ ≈ tomandel(dσdϵ) + @test diff.dσdp[:,1] ≈ tomandel(dσdG) + @test diff.dσdp[:,2] ≈ tomandel(dσdK) + end + function test_derivative(m, ϵ, state, Δt; comparesettings = (), numdiffsettings = (), diff = MaterialDerivatives(m)) + diff_num = MaterialDerivatives(m) + copy!(diff_num.dsdp, diff.dsdp) + + cache = allocate_material_cache(m) + extras = allocate_differentiation_output(m) + _, dσdϵ, _ = material_response(m, ϵ, state, Δt, cache, extras) + differentiate_material!(diff, m, ϵ, state, Δt, cache, extras, dσdϵ) + TestMaterials.obtain_numerical_material_derivative!(diff_num, m, ϵ, state, Δt; numdiffsettings...) + for key in fieldnames(typeof(diff)) + @testset "$key" begin + check = isapprox(getfield(diff, key), getfield(diff_num, key); comparesettings...) + @test check + if !check + println("Printing derivative, numerical derivative, and relative diff of each term") + display(getfield(diff, key)) + display(getfield(diff_num, key)) + display((getfield(diff, key) .- getfield(diff_num, key)) ./ (1e-100 .+ abs.(getfield(diff, key)) + abs.(getfield(diff_num, key)))) + println("Done printing for failed test above") + end + end + end + end + + @testset "Numerical (3D material)" begin + elastic = LinearElastic(0.52, 0.77) + viscoelastic = ViscoElastic(elastic, LinearElastic(0.33, 0.54), 0.36) + for m in (elastic, viscoelastic) + @testset "Initial response" begin + ϵ = rand(SymmetricTensor{2,3}) * 1e-3 + Δt = 1e-2 + test_derivative(m, ϵ, initial_material_state(m), Δt; + numdiffsettings = (fdtype = Val{:central},), + comparesettings = () + ) + diff = MaterialDerivatives(m) + copy!(diff.dsdp, rand(size(diff.dsdp)...)) + test_derivative(m, ϵ, initial_material_state(m), Δt; + numdiffsettings = (fdtype = Val{:central},), + comparesettings = (), + diff) + end + @testset "After shear loading" begin + ϵ21 = 0.01; num_steps = 10; t_end = 0.01 + stressfun(p) = TestMaterials.runstrain(fromvector(p, m), ϵ21, (2, 1), t_end, num_steps)[1] + dσ21_dp_num = FiniteDiff.finite_difference_jacobian(stressfun, tovector(m), Val{:central}; relstep = 1e-6) + σv, state, dσ21_dp, diff = TestMaterials.runstrain_diff(m, ϵ21, (2, 1), t_end, num_steps) + @test σv ≈ stressfun(tovector(m)) + @test dσ21_dp ≈ dσ21_dp_num + end + @testset "FullStressState" begin + ϵ21 = 0.01; num_steps = 10; t_end = 0.01 + ss = FullStressState() + # Check that we get the same result for runstresstate and runstrain + σ_ss, s_ss = TestMaterials.runstresstate(ss, m, ϵ21, (2, 1), t_end, num_steps) + σ, s = TestMaterials.runstrain(m, ϵ21, (2, 1), t_end, num_steps) + @test σ_ss ≈ σ + @test tovector(s_ss) ≈ tovector(s) + end + for (stress_state, ij) in ( + (UniaxialStress(), (1,1)), (UniaxialStrain(), (1,1)), + (UniaxialNormalStress(), (1,1)), (UniaxialNormalStress(), (2,1)), + (PlaneStress(), (2, 2)), (PlaneStrain(), (2, 1)), + (GeneralStressState(SymmetricTensor{2,3,Bool}((true, false, false, false, true, true)), rand(SymmetricTensor{2,3})), (2,2)) + ) + @testset "$(nameof(typeof(stress_state))), (i,j) = ($(ij[1]), $(ij[2]))" begin + ϵij = 0.01; num_steps = 1; t_end = 0.01 + stressfun(p) = TestMaterials.runstresstate(stress_state, fromvector(p, m), ϵij, ij, t_end, num_steps)[1] + dσij_dp_num = FiniteDiff.finite_difference_jacobian(stressfun, tovector(m), Val{:central}; relstep = 1e-6) + σv, state, dσij_dp, diff = TestMaterials.runstresstate_diff(stress_state, m, ϵij, ij, t_end, num_steps) + @test dσij_dp ≈ dσij_dp_num + end + end + end + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 823f7bf..b84f8ed 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,7 +2,8 @@ using MaterialModelsBase using Test using Tensors, StaticArrays import MaterialModelsBase as MMB -using TestMaterials: LinearElastic, ViscoElastic +using TestMaterials: TestMaterials, LinearElastic, ViscoElastic +using FiniteDiff: FiniteDiff include("utils4testing.jl") From 44a6dd485e48c1ba95f5752c3880d827b798c735 Mon Sep 17 00:00:00 2001 From: Knut Andreas Date: Sun, 2 Nov 2025 12:40:20 +0100 Subject: [PATCH 07/14] Make interface stricter --- Project.toml | 6 +++--- src/MaterialModelsBase.jl | 7 ++++--- src/vector_conversion.jl | 3 ++- test/differentiation.jl | 37 +++++++------------------------------ test/runtests.jl | 4 +++- 5 files changed, 19 insertions(+), 38 deletions(-) diff --git a/Project.toml b/Project.toml index 338a28e..a50ec63 100644 --- a/Project.toml +++ b/Project.toml @@ -13,11 +13,11 @@ julia = "1.11" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -TestMaterials = "aa4f4413-326f-4ddd-a53e-c58801ebfeaf" FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" +MaterialModelsTesting = "882b014b-b96c-4115-8629-e17fb35110d2" [sources] -TestMaterials = {path = "test/TestMaterials"} +MaterialModelsTesting = {url = "https://github.com/KnutAM/MaterialModelsTesting.jl"} [targets] -test = ["Test", "FiniteDiff", "TestMaterials"] +test = ["Test", "FiniteDiff", "MaterialModelsTesting"] diff --git a/src/MaterialModelsBase.jl b/src/MaterialModelsBase.jl index 03a10e8..567096e 100644 --- a/src/MaterialModelsBase.jl +++ b/src/MaterialModelsBase.jl @@ -116,12 +116,13 @@ of the material `m`. Defaults to the empty `NoMaterialState()` """ -function initial_material_state(::AbstractMaterial) - return NoMaterialState() +function initial_material_state(m::AbstractMaterial) + T = get_params_eltype(m) + return NoMaterialState{T}() end abstract type AbstractMaterialState end -struct NoMaterialState <: AbstractMaterialState end +struct NoMaterialState{T} <: AbstractMaterialState end # Material cache """ diff --git a/src/vector_conversion.jl b/src/vector_conversion.jl index 55c0b55..22099d5 100644 --- a/src/vector_conversion.jl +++ b/src/vector_conversion.jl @@ -49,6 +49,7 @@ Get the type used to store each scalar component of the material state variables defaults to `Float64`. """ get_statevar_eltype(::AbstractMaterialState) = Float64 +get_statevar_eltype(::NoMaterialState{T}) where {T} = T """ get_num_params(m::AbstractMaterial) @@ -109,7 +110,7 @@ Create a tensor with shape of `TT` with entries from the Mandel components in `v """ function fromvector end -fromvector(::AbstractVector, ::NoMaterialState; kwargs...) = NoMaterialState() +fromvector(::AbstractVector{T}, ::NoMaterialState; kwargs...) where {T} = NoMaterialState{T}() # Tensors.jl implementation function fromvector(v::AbstractVector, ::TT; offset = 0) where {TT <: SecondOrderTensor} diff --git a/test/differentiation.jl b/test/differentiation.jl index e3af6ad..194da6a 100644 --- a/test/differentiation.jl +++ b/test/differentiation.jl @@ -38,29 +38,6 @@ end @test diff.dσdp[:,1] ≈ tomandel(dσdG) @test diff.dσdp[:,2] ≈ tomandel(dσdK) end - function test_derivative(m, ϵ, state, Δt; comparesettings = (), numdiffsettings = (), diff = MaterialDerivatives(m)) - diff_num = MaterialDerivatives(m) - copy!(diff_num.dsdp, diff.dsdp) - - cache = allocate_material_cache(m) - extras = allocate_differentiation_output(m) - _, dσdϵ, _ = material_response(m, ϵ, state, Δt, cache, extras) - differentiate_material!(diff, m, ϵ, state, Δt, cache, extras, dσdϵ) - TestMaterials.obtain_numerical_material_derivative!(diff_num, m, ϵ, state, Δt; numdiffsettings...) - for key in fieldnames(typeof(diff)) - @testset "$key" begin - check = isapprox(getfield(diff, key), getfield(diff_num, key); comparesettings...) - @test check - if !check - println("Printing derivative, numerical derivative, and relative diff of each term") - display(getfield(diff, key)) - display(getfield(diff_num, key)) - display((getfield(diff, key) .- getfield(diff_num, key)) ./ (1e-100 .+ abs.(getfield(diff, key)) + abs.(getfield(diff_num, key)))) - println("Done printing for failed test above") - end - end - end - end @testset "Numerical (3D material)" begin elastic = LinearElastic(0.52, 0.77) @@ -82,9 +59,9 @@ end end @testset "After shear loading" begin ϵ21 = 0.01; num_steps = 10; t_end = 0.01 - stressfun(p) = TestMaterials.runstrain(fromvector(p, m), ϵ21, (2, 1), t_end, num_steps)[1] + stressfun(p) = runstrain(fromvector(p, m), ϵ21, (2, 1), t_end, num_steps)[1] dσ21_dp_num = FiniteDiff.finite_difference_jacobian(stressfun, tovector(m), Val{:central}; relstep = 1e-6) - σv, state, dσ21_dp, diff = TestMaterials.runstrain_diff(m, ϵ21, (2, 1), t_end, num_steps) + σv, state, dσ21_dp, diff = runstrain_diff(m, ϵ21, (2, 1), t_end, num_steps) @test σv ≈ stressfun(tovector(m)) @test dσ21_dp ≈ dσ21_dp_num end @@ -92,8 +69,8 @@ end ϵ21 = 0.01; num_steps = 10; t_end = 0.01 ss = FullStressState() # Check that we get the same result for runstresstate and runstrain - σ_ss, s_ss = TestMaterials.runstresstate(ss, m, ϵ21, (2, 1), t_end, num_steps) - σ, s = TestMaterials.runstrain(m, ϵ21, (2, 1), t_end, num_steps) + σ_ss, s_ss = runstresstate(ss, m, ϵ21, (2, 1), t_end, num_steps) + σ, s = runstrain(m, ϵ21, (2, 1), t_end, num_steps) @test σ_ss ≈ σ @test tovector(s_ss) ≈ tovector(s) end @@ -104,10 +81,10 @@ end (GeneralStressState(SymmetricTensor{2,3,Bool}((true, false, false, false, true, true)), rand(SymmetricTensor{2,3})), (2,2)) ) @testset "$(nameof(typeof(stress_state))), (i,j) = ($(ij[1]), $(ij[2]))" begin - ϵij = 0.01; num_steps = 1; t_end = 0.01 - stressfun(p) = TestMaterials.runstresstate(stress_state, fromvector(p, m), ϵij, ij, t_end, num_steps)[1] + ϵij = 0.01; num_steps = 10; t_end = 0.01 + stressfun(p) = runstresstate(stress_state, fromvector(p, m), ϵij, ij, t_end, num_steps)[1] dσij_dp_num = FiniteDiff.finite_difference_jacobian(stressfun, tovector(m), Val{:central}; relstep = 1e-6) - σv, state, dσij_dp, diff = TestMaterials.runstresstate_diff(stress_state, m, ϵij, ij, t_end, num_steps) + σv, state, dσij_dp, diff = runstresstate_diff(stress_state, m, ϵij, ij, t_end, num_steps) @test dσij_dp ≈ dσij_dp_num end end diff --git a/test/runtests.jl b/test/runtests.jl index b84f8ed..9264bc2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,8 +2,10 @@ using MaterialModelsBase using Test using Tensors, StaticArrays import MaterialModelsBase as MMB -using TestMaterials: TestMaterials, LinearElastic, ViscoElastic using FiniteDiff: FiniteDiff +using MaterialModelsTesting: + LinearElastic, ViscoElastic, test_derivative, obtain_numerical_material_derivative!, + runstrain, runstrain_diff, runstresstate, runstresstate_diff include("utils4testing.jl") From 0e95001fb73f7ccf8003483ab877dce43bc7770b Mon Sep 17 00:00:00 2001 From: Knut Andreas Date: Sun, 2 Nov 2025 14:45:11 +0100 Subject: [PATCH 08/14] Provide type when creating material state --- src/MaterialModelsBase.jl | 3 ++- test/differentiation.jl | 2 +- test/errors.jl | 2 +- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/MaterialModelsBase.jl b/src/MaterialModelsBase.jl index 567096e..84a410a 100644 --- a/src/MaterialModelsBase.jl +++ b/src/MaterialModelsBase.jl @@ -114,7 +114,8 @@ end Return the (default) initial `state::AbstractMaterialState` of the material `m`. -Defaults to the empty `NoMaterialState()` +Defaults to the empty `NoMaterialState{T}()`, where +`T = get_params_eltype(m)` (which defaults to `Float64`). """ function initial_material_state(m::AbstractMaterial) T = get_params_eltype(m) diff --git a/test/differentiation.jl b/test/differentiation.jl index 194da6a..4543c5e 100644 --- a/test/differentiation.jl +++ b/test/differentiation.jl @@ -31,7 +31,7 @@ end dσdG = 2*(one(SymmetricTensor{4,3}) - IxI/3) ⊡ ϵ dσdK = IxI ⊡ ϵ - old = NoMaterialState(); Δt = nothing + old = NoMaterialState{Float64}(); Δt = nothing cache = NoMaterialCache(); extras = allocate_differentiation_output(m) differentiate_material!(diff, m, ϵ, old, Δt, cache, extras, dσdϵ) @test diff.dσdϵ ≈ tomandel(dσdϵ) diff --git a/test/errors.jl b/test/errors.jl index 9743442..94eadd5 100644 --- a/test/errors.jl +++ b/test/errors.jl @@ -10,5 +10,5 @@ m = LinearElastic(80e3, 150e3) ss = PlaneStress(;max_iter = 1, tolerance=1e-100) ϵ = rand(SymmetricTensor{2,2}) - @test_throws NoStressConvergence material_response(ss, m, ϵ, NoMaterialState()) + @test_throws NoStressConvergence material_response(ss, m, ϵ, NoMaterialState{Float64}()) end \ No newline at end of file From 5c39918509ad929c5d3a90e4c6be8d62efbcef2a Mon Sep 17 00:00:00 2001 From: Knut Andreas Date: Sun, 2 Nov 2025 21:49:59 +0100 Subject: [PATCH 09/14] Fix stress iterations --- src/stressiterations.jl | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/src/stressiterations.jl b/src/stressiterations.jl index 932b43a..a49e0d8 100644 --- a/src/stressiterations.jl +++ b/src/stressiterations.jl @@ -412,10 +412,17 @@ end # GeneralStressState function get_full_tensor(state::GeneralStressState{Nσ}, ::TT, v::SVector{Nσ,T}) where {Nσ,T,TT} - shear_factor = T(1/√2) + TB = Tensors.get_base(TT) + shear_factor = if TB == SymmetricTensor{2,3} + T(1/√2) + elseif TB == Tensor{2,3} + one(T) + else + error("GeneralStressState expects full dimension of strain input") + end s(i,j) = i==j ? one(T) : shear_factor f(i,j) = state.σ_ctrl[i,j] ? v[state.σ_minds[i,j]]*s(i,j) : zero(T) - return Tensors.get_base(TT)(f) + return TB(f) end function get_unknowns(state::GeneralStressState{Nσ}, a::AbstractTensor{2,3,T}) where {Nσ, T} @@ -426,7 +433,7 @@ function get_unknowns(state::GeneralStressState{Nσ}, a::AbstractTensor{2,3,T}) end function get_unknowns(state::GeneralStressState{Nσ}, a::AbstractTensor{4,3,T}) where {Nσ,T} - shear_factor = T(√2) + shear_factor = a isa SymmetricTensor ? T(√2) : one(T) s(i,j) = i==j ? one(T) : shear_factor f(c1,c2) = ((i,j) = c1; (k,l) = c2; a[i,j,k,l]*s(i,j)*s(k,l)) return SMatrix{Nσ,Nσ,T}((f(c1,c2) for c1 in state.σm_inds, c2 in state.σm_inds)) From 888af92586bb3b1c45c8c2c7dd6a97d44bb15474 Mon Sep 17 00:00:00 2001 From: Knut Andreas Date: Sun, 2 Nov 2025 22:50:31 +0100 Subject: [PATCH 10/14] Better test coverage for GeneralStressState, and fix following move to MaterialModelsTesting.jl of test features --- test/runtests.jl | 2 +- test/stressiterations.jl | 78 ++++++++++++++++++++++++---------------- 2 files changed, 49 insertions(+), 31 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 9264bc2..46e747c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,7 +4,7 @@ using Tensors, StaticArrays import MaterialModelsBase as MMB using FiniteDiff: FiniteDiff using MaterialModelsTesting: - LinearElastic, ViscoElastic, test_derivative, obtain_numerical_material_derivative!, + LinearElastic, ViscoElastic, NeoHooke, test_derivative, obtain_numerical_material_derivative!, runstrain, runstrain_diff, runstresstate, runstresstate_diff include("utils4testing.jl") diff --git a/test/stressiterations.jl b/test/stressiterations.jl index 5169f2f..bc474bd 100644 --- a/test/stressiterations.jl +++ b/test/stressiterations.jl @@ -1,18 +1,27 @@ all_states = ( FullStressState(), UniaxialStrain(), PlaneStrain(), UniaxialStress(), PlaneStress(), UniaxialNormalStress()) -iter_states = ( UniaxialStress(), PlaneStress(), UniaxialNormalStress()) +gss_ctrl = SymmetricTensor{2,3,Bool}((true, true, falses(4)...)) +iter_states = ( UniaxialStress(), PlaneStress(), UniaxialNormalStress(), + (SymmetricTensor = GeneralStressState(gss_ctrl, 0.0 * gss_ctrl), + Tensor = GeneralStressState(Tensor{2,3}(gss_ctrl), 0.0 * Tensor{2,3}(gss_ctrl))) + ) iter_mandel = ( ([2,3,4,5,6,7,8,9],[2,3,4,5,6]), ([3,4,5,7,8], [3,4,5]), - ([2,3], [2,3])) - + ([2,3], [2,3]), + ([1,6,9], [1,6])) @testset "conversions" begin # Test that # 1) All conversions are invertible (convert back and fourth) # 2) Are compatible with the expected mandel dynamic tensors - for stress_state in all_states + for _stress_state in all_states for TT in (Tensor, SymmetricTensor) + stress_state = if isa(_stress_state, NamedTuple) + _stress_state[nameof(TT)] + else + _stress_state + end ared = rand(TT{2,_getdim(stress_state)}) afull = MMB.expand_tensordim(stress_state, ared) @test _getdim(afull) == 3 # Full dimensional tensor @@ -26,8 +35,13 @@ iter_mandel = ( ([2,3,4,5,6,7,8,9],[2,3,4,5,6]), end end - for (state, inds) in zip(iter_states, iter_mandel) + for (_stress_state, inds) in zip(iter_states, iter_mandel) for (TT, ii) in zip((Tensor, SymmetricTensor), inds) + state = if isa(_stress_state, NamedTuple) + _stress_state[nameof(TT)] + else + _stress_state + end a = rand(TT{2,3}) A = rand(TT{4,3}) ax = MMB.get_unknowns(state, a) @@ -46,8 +60,13 @@ iter_mandel = ( ([2,3,4,5,6,7,8,9],[2,3,4,5,6]), end @testset "stiffness_calculations" begin - for (state, inds) in zip(iter_states, iter_mandel) + for (_stress_state, inds) in zip(iter_states, iter_mandel) for (TT, ii) in zip((Tensor, SymmetricTensor), inds) + state = if isa(_stress_state, NamedTuple) + _stress_state[nameof(TT)] + else + _stress_state + end dσdϵ = rand(TT{4,3}) + one(TT{4,3}) dσᶜdϵᶜ = MMB.reduce_stiffness(state, dσdϵ) if _getdim(state) < 3 # Otherwise, conversion is direct @@ -150,30 +169,29 @@ end Δϵ = 1e-6 # Small strain to compare with linear case rtol = 1e-5 # Relative tolerance to compare with linear case - @testset for m in (StVenant(;G, K), NeoHooke(;G, K)) - old = initial_material_state(m) - # UniaxialStress - P, dPdF, state, Ffull = material_response(UniaxialStress(), m, Tensor{2,1}((1 + Δϵ,)), old, 0.0) - @test isapprox(P[1,1], E*Δϵ; rtol) - @test isapprox(dPdF[1,1,1,1], E; rtol) - @test Ffull[2,2] ≈ Ffull[3,3] - @test Ffull[2,2] ≈ 1-ν*Δϵ - - # UniaxialStrain - P, dPdF, state, Ffull = material_response(UniaxialStrain(), m, Tensor{2,1}((1 + Δϵ,)), old, 0.0) - @test isapprox(P[1,1], (2G+λ)*Δϵ; rtol) - @test isapprox(dPdF[1,1,1,1], (2G+λ); rtol) - P_full, dPdF_full, _ = material_response(m, Ffull, old, 0.0) - @test isapprox(P_full[2,2], λ*Δϵ; rtol) - - # PlaneStress - ϵ = Δϵ * rand(SymmetricTensor{2,2}) - F = one(Tensor{2,2}) + ϵ - P, dPdF, state, Ffull = material_response(PlaneStress(), m, F, old, 0.0) - Dvoigt = (E/(1-ν^2))*[1 ν 0; ν 1 0; 0 0 (1-ν)/2] - @test isapprox(tovoigt(symmetric(dPdF)), Dvoigt; rtol) - @test isapprox(P, fromvoigt(SymmetricTensor{4,2}, Dvoigt)⊡ϵ; rtol) - end + m = NeoHooke(;G, K) + old = initial_material_state(m) + # UniaxialStress + P, dPdF, state, Ffull = material_response(UniaxialStress(), m, Tensor{2,1}((1 + Δϵ,)), old, 0.0) + @test isapprox(P[1,1], E*Δϵ; rtol) + @test isapprox(dPdF[1,1,1,1], E; rtol) + @test Ffull[2,2] ≈ Ffull[3,3] + @test Ffull[2,2] ≈ 1-ν*Δϵ + + # UniaxialStrain + P, dPdF, state, Ffull = material_response(UniaxialStrain(), m, Tensor{2,1}((1 + Δϵ,)), old, 0.0) + @test isapprox(P[1,1], (2G+λ)*Δϵ; rtol) + @test isapprox(dPdF[1,1,1,1], (2G+λ); rtol) + P_full, dPdF_full, _ = material_response(m, Ffull, old, 0.0) + @test isapprox(P_full[2,2], λ*Δϵ; rtol) + + # PlaneStress + ϵ = Δϵ * rand(SymmetricTensor{2,2}) + F = one(Tensor{2,2}) + ϵ + P, dPdF, state, Ffull = material_response(PlaneStress(), m, F, old, 0.0) + Dvoigt = (E/(1-ν^2))*[1 ν 0; ν 1 0; 0 0 (1-ν)/2] + @test isapprox(tovoigt(symmetric(dPdF)), Dvoigt; rtol) + @test isapprox(P, fromvoigt(SymmetricTensor{4,2}, Dvoigt)⊡ϵ; rtol) end @testset "visco_elastic" begin From 2f2f41fe0b3e97747a9b031fa5d492f701d4ff7c Mon Sep 17 00:00:00 2001 From: Knut Andreas Date: Sun, 2 Nov 2025 23:18:43 +0100 Subject: [PATCH 11/14] Make docs appear for stress state derivative updates --- src/differentiation.jl | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index 69b030d..bb03fce 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -76,6 +76,19 @@ allocate_differentiation_output(::AbstractMaterial) = NoExtraOutput() Calculate the derivatives and save them in `diff`, see [`MaterialDerivatives`](@ref) for a description of the fields in `diff`. + + differentiate_material!(ssd::StressStateDerivatives, stress_state, m, args...) + +For material models implementing `material_response(m, args...)` and `differentiate_material!(::MaterialDerivatives, m, args...)`, +this method will work automatically by +1) Calling `σ, dσdϵ, state = material_response(stress_state, m, args...)` (except that `dσdϵ::FourthOrderTensor{dim = 3}` is extracted) +2) Calling `differentiate_material!(ssd.mderiv::MaterialDerivatives, m, args..., dσdϵ::FourthOrderTensor{3})` +3) Updating `ssd` according to the constraints imposed by the `stress_state`. + +For material models that directly implement `material_response(stress_state, m, args...)`, this function should be overloaded directly +to calculate the derivatives in `ssd`. Here the user has full control and no modifications occur automatically, however, typically the +(total) derivatives `ssd.dσdp`, `ssd.dϵdp`, and `ssd.mderiv.dsdp` should be updated. +The exact interface of the `StressStateDerivatives` datastructure is still not fixed, and may change in the future. """ function differentiate_material! end @@ -104,19 +117,6 @@ function StressStateDerivatives(::AbstractStressState, m::AbstractMaterial) return StressStateDerivatives(mderiv, dϵdp, dσdp, index, index) end -""" - differentiate_material!(ssd::StressStateDerivatives, stress_state, m, args...) - -For material models implementing `material_response(m, args...)` and `differentiate_material!(::MaterialDerivatives, m, args...)`, -this method will work automatically by -1) Calling `σ, dσdϵ, state = material_response(stress_state, m, args...)` (except that `dσdϵ::FourthOrderTensor{dim = 3}` is extracted) -2) Calling `differentiate_material!(ssd.mderiv::MaterialDerivatives, m, args..., dσdϵ::FourthOrderTensor{3})` -3) Updating `ssd` according to the constraints imposed by the `stress_state`. - -For material models that directly implement `material_response(stress_state, m, args...)`, this function should be overloaded directly -to calculate the derivatives in `ssd`. Here the user has full control and no modifications occur automatically, however, typically the -(total) derivatives `ssd.dσdp`, `ssd.dϵdp`, and `ssd.mderiv.dsdp` should be updated. -""" function differentiate_material!(ssd::StressStateDerivatives, stress_state::AbstractStressState, m::AbstractMaterial, ϵ::AbstractTensor, args::Vararg{Any,N}) where {N} σ_full, dσdϵ_full, state, ϵ_full = stress_state_material_response(stress_state, m, ϵ, args...) differentiate_material!(ssd.mderiv, m, ϵ_full, args..., dσdϵ_full) From 8f06e7775fb735c16f7f44f230ca6bf189b3cc6a Mon Sep 17 00:00:00 2001 From: Knut Andreas Date: Tue, 4 Nov 2025 20:00:39 +0100 Subject: [PATCH 12/14] Define consistent conversion routine naming --- Project.toml | 1 + docs/src/conversion.md | 16 +++- src/MaterialModelsBase.jl | 6 +- src/differentiation.jl | 10 +-- src/stressiterations.jl | 8 +- src/vector_conversion.jl | 173 ++++++++++++++++---------------------- 6 files changed, 99 insertions(+), 115 deletions(-) diff --git a/Project.toml b/Project.toml index ba9b275..0a0f2ca 100644 --- a/Project.toml +++ b/Project.toml @@ -18,6 +18,7 @@ MaterialModelsTesting = "882b014b-b96c-4115-8629-e17fb35110d2" [sources] MaterialModelsTesting = {url = "https://github.com/KnutAM/MaterialModelsTesting.jl"} +#MaterialModelsTesting = {path = "/Users/knutan/.julia/dev/MaterialModelsTesting"} [targets] test = ["Test", "FiniteDiff", "MaterialModelsTesting"] diff --git a/docs/src/conversion.md b/docs/src/conversion.md index 0f841ee..3a550ed 100644 --- a/docs/src/conversion.md +++ b/docs/src/conversion.md @@ -2,25 +2,33 @@ CurrentModule = MaterialModelsBase ``` # Conversion -`MaterialModelsBase` defines an interface for converting parameters and state variables +`MaterialModelsBase` defines an interface for converting parameters, state variables, and other objects to and from `AbstractVector`s. This is useful when doing parameter identification, or when interfacing with different languages that require information to be passed as arrays. These function can be divided into those providing information about the material and state variables, and those doing the actual conversion. ## Information functions +The following should be defined for new materials, ```@docs get_tensorbase +get_vector_length +get_vector_eltype +``` + +Whereas the following already have default implementations that should work provided that the above are implemented. +```@docs get_num_statevars -get_statevar_eltype -get_num_params -get_params_eltype get_num_tensorcomponents ``` ## Conversion functions +The following should be defined for new materials, state variables, and alike ```@docs tovector! fromvector +``` +whereas the following already have default implemenations that should work provided that the above functions are implemented, +```@docs tovector ``` diff --git a/src/MaterialModelsBase.jl b/src/MaterialModelsBase.jl index 84a410a..68fd10b 100644 --- a/src/MaterialModelsBase.jl +++ b/src/MaterialModelsBase.jl @@ -25,7 +25,7 @@ export update_stress_state! # For nonzero stress # 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 -export get_num_params, get_params_eltype # +export get_vector_length, get_vector_eltype # Get the length and type when converting objects to vectors export MaterialDerivatives, StressStateDerivatives # Derivative collections export differentiate_material! # Differentiation routines export allocate_differentiation_output # @@ -115,10 +115,10 @@ Return the (default) initial `state::AbstractMaterialState` of the material `m`. Defaults to the empty `NoMaterialState{T}()`, where -`T = get_params_eltype(m)` (which defaults to `Float64`). +`T = get_vector_eltype(m)` (which defaults to `Float64`). """ function initial_material_state(m::AbstractMaterial) - T = get_params_eltype(m) + T = get_vector_eltype(m) return NoMaterialState{T}() end diff --git a/src/differentiation.jl b/src/differentiation.jl index bb03fce..fa3acfe 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -19,8 +19,8 @@ end MaterialDerivatives(m::AbstractMaterial) A struct that saves all derivative information using a `Matrix{T}` for each derivative, -where `T=get_params_eltype(m)`. The dimensions are obtained from `get_num_tensorcomponents`, -`get_num_statevars`, and `get_num_params`. `m` must support `tovector` and `fromvector`, while +where `T=get_vector_eltype(m)`. The dimensions are obtained from `get_num_tensorcomponents`, +`get_num_statevars`, and `get_vector_length`. `m` must support `tovector` and `fromvector`, while the output of `initial_material_state` must support `tovector`, and in addition the element type of `tovector(initial_material_state(m))` must respect the element type in `tovector(m)` for any `m`. @@ -37,10 +37,10 @@ and old state variables, and `p` the material parameter vector. """ function MaterialDerivatives(m::AbstractMaterial) - T = get_params_eltype(m) + T = get_vector_eltype(m) n_tensor = get_num_tensorcomponents(m) n_state = get_num_statevars(m) - n_params = get_num_params(m) + n_params = get_vector_length(m) dsdp = ForwardDiff.jacobian(p -> tovector(initial_material_state(fromvector(p, m))), tovector(m)) return MaterialDerivatives( zeros(T, n_tensor, n_tensor), # dσdϵ @@ -103,7 +103,7 @@ end function StressStateDerivatives(::AbstractStressState, m::AbstractMaterial) mderiv = MaterialDerivatives(m) - np = get_num_params(m) + np = get_vector_length(m) nt = get_num_tensorcomponents(m) # Should be changed to only save non-controlled entries dϵdp = zeros(nt, np) dσdp = zeros(nt, np) diff --git a/src/stressiterations.jl b/src/stressiterations.jl index b833ef7..dad4a2c 100644 --- a/src/stressiterations.jl +++ b/src/stressiterations.jl @@ -38,8 +38,8 @@ Creates a subtype of `AbstractMaterial` that wraps a stress state and a material calls to `material_response(w::ReducedStressState, args...)` gives the same result as `material_response(s, m, args...)`. Calls to `initial_material_state`, `allocate_material_cache`, -`get_num_tensorcomponents`, `get_num_statevars`, `get_num_params`, -`get_params_eltype`, `tovector!`, `tovector`, +`get_num_tensorcomponents`, `get_num_statevars`, `get_vector_length`, +`get_vector_eltype`, `tovector!`, `tovector`, and `allocate_differentiation_output` are forwarded with `m` as the argument. `fromvector` returns `ReducedStressState` and is supported as well. """ @@ -48,8 +48,8 @@ struct ReducedStressState{S<:AbstractStressState,M<:AbstractMaterial} <: Abstrac material::M end for op in ( :initial_material_state, :allocate_material_cache, - :get_num_tensorcomponents, :get_num_statevars, :get_num_params, - :get_params_eltype, :allocate_differentiation_output) + :get_num_tensorcomponents, :get_num_statevars, :get_vector_length, + :get_vector_eltype, :allocate_differentiation_output) @eval @inline $op(rss::ReducedStressState) = $op(rss.material) end function tovector!(v::AbstractVector, rss::ReducedStressState) diff --git a/src/vector_conversion.jl b/src/vector_conversion.jl index 65a089c..ff1f359 100644 --- a/src/vector_conversion.jl +++ b/src/vector_conversion.jl @@ -1,92 +1,115 @@ - """ - get_tensorbase(m::AbstractMaterial) + tovector!(v::AbstractVector, obj; offset = 0) -Return the type of the primary input (strain-like) and associated output (stress-like) -to the material model. The default is `SymmetricTensor{2, 3}` for small-strain material -models in 3D, but it can be -* Small strain material model: Small strain tensor, `ϵ::SymmetricTensor{2, dim}` -* Large strain material model: Deformation gradient, `F::Tensor{2, dim}` -* Traction-separation law: Displacement jump, `Δu::Vec{dim}` +Store the parameters in the object `obj` in `v`. This is typically used with `obj` as +an `AbstractMaterial`, an `AbstractMaterialState`, or an `AbstractTensor`. The `offset` input makes it possible +to write values starting at an offset location in `v`. -Where `dim = 3` is most common, but any value, 1, 2, or 3, is valid. +The implementation for `AbstractTensor`s is included in the `MaterialModelsBase.jl` package, and use the Mandel +notation for the conversion. """ -get_tensorbase(::AbstractMaterial) = SymmetricTensor{2, 3} +function tovector! end """ - get_num_tensorcomponents(::AbstractMaterial) + tovector([T::Type{<:AbstractArray}], obj)::T -Returns the number of independent components for the given material. - -!!! note - It is not required to implement this function, it is inferred by the - implementation of [`get_tensorbase`](@ref) +Out-of place version of `tovector!`. Relies on `get_vector_length` and +`get_vector_eltype` to be correctly defined, and defaults to `T = Array`. +Experimental support for `T = SArray` is also available, but performance may be suboptimal, +and can be improved by implementing a custom function for the given type. """ -get_num_tensorcomponents(m::AbstractMaterial) = Tensors.n_components(get_tensorbase(m)) +@inline tovector(obj) = tovector(Array, obj) -""" - get_num_statevars(s::AbstractMaterialState) - get_num_statevars(m::AbstractMaterial) +function tovector(::Type{<:Array}, obj) + T = get_vector_eltype(obj) + return tovector!(zeros(T, get_vector_length(obj)), obj) +end -Return the number of state variables. -A tensorial state variable should be counted by how many components it has. -E.g. if a state consists of one scalar and one symmetric 2nd order tensor in 3d, -`get_num_statevars` should return 7. +function tovector(::Type{<:SArray}, obj) + T = get_vector_eltype(obj) + N = get_vector_length(obj) + m = MVector{T, N}() + @inline tovector!(m, obj) + return SVector{T, N}(m) +end -It suffices to define for the state, `s`, only, but defining for material, `m`, -directly as well can improve performance. """ -function get_num_statevars end + fromvector(v::AbstractVector, ::OT; offset = 0) -get_num_statevars(::NoMaterialState) = 0 -get_num_statevars(m::AbstractMaterial) = get_num_statevars(initial_material_state(m)) +Output an object of similar type to `OT`, but with parameters according to `v`. This is typically used with +an `AbstractMaterial`, an `AbstractMaterialState`, or an `AbstractTensor`. The `offset` input makes it possible +to read values starting at an offset location in `v`. +The implementation for `AbstractTensor`s is included in the `MaterialModelsBase.jl` package, and use the Mandel +notation for the conversion. """ - get_statevar_eltype(s::AbstractMaterialState) +function fromvector end -Get the type used to store each scalar component of the material state variables, -defaults to `Float64`. """ -get_statevar_eltype(::AbstractMaterialState) = Float64 -get_statevar_eltype(::NoMaterialState{T}) where {T} = T + get_vector_length(obj) +Return the length of the vector representation of `obj` when using +`tovector` or `tovector!`. The default implementation of `tovector` relies +on this function being defined. """ - get_num_params(m::AbstractMaterial) +function get_vector_length end -Return the number of material parameters in `m`. No default value implemented. """ -function get_num_params end + get_vector_eltype(obj) +Return the element type of the vector representation of `obj` when using +`tovector`, i.e. with `T = get_vector_eltype(obj)` and `N = get_vector_length(obj)`, +we have `typeof(obj) == typeof(fromvector(zeros(T, N), obj))`. """ - get_params_eltype(m::AbstractMaterial) +function get_vector_eltype end -Return the number type for the scalar material parameters, defaults to `Float64` """ -get_params_eltype(::AbstractMaterial) = Float64 + get_tensorbase(m::AbstractMaterial) + +Return the type of the primary input (strain-like) and associated output (stress-like) +to the material model. The default is `SymmetricTensor{2, 3}` for small-strain material +models in 3D, but it can be +* Small strain material model: Small strain tensor, `ϵ::SymmetricTensor{2, dim}` +* Large strain material model: Deformation gradient, `F::Tensor{2, dim}` +* Traction-separation law: Displacement jump, `Δu::Vec{dim}` -# Conversion functions +Where `dim = 3` is most common, but any value, 1, 2, or 3, is valid. """ - tovector!(v::AbstractVector, m::AbstractMaterial; offset = 0) +get_tensorbase(::AbstractMaterial) = SymmetricTensor{2, 3} -Put the material parameters of `m` into the vector `v`. -This is typically used when the parameters should be fitted. +""" + get_num_tensorcomponents(::AbstractMaterial) - tovector!(v::AbstractVector, s::AbstractMaterialState; offset = 0) +Returns the number of independent components for the given material. -Put the state variables in `s` into the vector `v`. -This is typically used when differentiating the material -wrt. the the old state variables. +!!! note + It is not required to implement this function, it is inferred by the + implementation of [`get_tensorbase`](@ref) - tovector!(v::AbstractVector, a::Union{SecondOrderTensor, FourthOrderTensor}; offset = 0) +""" +get_num_tensorcomponents(m::AbstractMaterial) = Tensors.n_components(get_tensorbase(m)) -Puts the Mandel components of `a` into the vector `v`. """ -function tovector! end + get_num_statevars(m::AbstractMaterial) + +Return the number of scalar values required to store the state of `m`, i.e. `length(tovector(initial_material_state(m)))`. +The default implementation works provided that `s = initial_material_state(m)` and `get_vector_length(s)` is defined. +""" +function get_num_statevars(m::AbstractMaterial) + return get_vector_length(initial_material_state(m)) +end +# Implementations for types defined in MaterialModelsBase tovector!(v::AbstractVector, ::NoMaterialState; kwargs...) = v +fromvector(::AbstractVector{T}, ::NoMaterialState; kwargs...) where {T} = NoMaterialState{T}() +get_vector_length(::NoMaterialState) = 0 +get_vector_eltype(::NoMaterialState{T}) where {T} = T # Tensors.jl implementation +get_vector_length(::TT) where {TT <: AbstractTensor} = Tensors.n_components(Tensors.get_base(TT)) +get_vector_eltype(a::AbstractTensor) = eltype(a) + function tovector!(v::AbstractVector, a::SecondOrderTensor; offset = 0) return tomandel!(v, a; offset) end @@ -97,22 +120,6 @@ function tovector!(v::AbstractVector, a::Union{Tensor{4, <:Any, <:Any, M}, Symme return v end -""" - fromvector(v::AbstractVector, ::MT; offset = 0) where {MT<:AbstractMaterial} -Create a material of type `MT` with the parameters according to `v` - - fromvector(v::AbstractVector, ::ST; offset = 0) where {ST<:AbstractMaterialState} -Create a material state of type `ST` with the values according to `v` - - fromvector(v, ::TT; offset = 0) where {TT <: Union{SecondOrderTensor, FourthOrderTensor}} - -Create a tensor with shape of `TT` with entries from the Mandel components in `v`. -""" -function fromvector end - -fromvector(::AbstractVector{T}, ::NoMaterialState; kwargs...) where {T} = NoMaterialState{T}() - -# Tensors.jl implementation function fromvector(v::AbstractVector, ::TT; offset = 0) where {TT <: SecondOrderTensor} return frommandel(Tensors.get_base(TT), v; offset) end @@ -123,35 +130,3 @@ function fromvector(v::AbstractVector, ::TT; offset = 0) where {TT <: FourthOrde N = round(Int, sqrt(M)) return frommandel(TB, reshape(view(v, offset .+ (1:M)), (N, N))) end - -""" - tovector(m::AbstractMaterial) - -Out-of place version of `tovector!`. Relies on `get_num_params` and -`get_params_eltype` to be correctly defined -""" -function tovector(m::AbstractMaterial) - T = get_params_eltype(m) - return tovector!(zeros(T, get_num_params(m)), m) -end - -""" - tovector(m::AbstractMaterialState) - -Out-of place version of `tovector!`. Relies on `get_num_statevars` and -`get_statevar_eltype` to be correctly defined -""" -function tovector(s::AbstractMaterialState) - T = get_statevar_eltype(s) - return tovector!(zeros(T, get_num_statevars(s)), s) -end - -tovector(a::Union{SecondOrderTensor, FourthOrderTensor}) = (m = tomandel(a); reshape(m, length(m))) - -# Backwards compatibility -const get_parameter_type = get_params_eltype -const material2vector! = tovector! -const material2vector = tovector -const vector2material = fromvector - -export material2vector!, material2vector, vector2material From 5ca86519578db12d8265e83a03945076c846483e Mon Sep 17 00:00:00 2001 From: Knut Andreas Date: Tue, 4 Nov 2025 20:16:59 +0100 Subject: [PATCH 13/14] Fix spelling --- docs/src/conversion.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/conversion.md b/docs/src/conversion.md index 3a550ed..19a2a45 100644 --- a/docs/src/conversion.md +++ b/docs/src/conversion.md @@ -28,7 +28,7 @@ The following should be defined for new materials, state variables, and alike tovector! fromvector ``` -whereas the following already have default implemenations that should work provided that the above functions are implemented, +whereas the following already have default implementations that should work provided that the above functions are implemented, ```@docs tovector ``` From 0941ef2f399be9e4786884619766122b38bdfdf8 Mon Sep 17 00:00:00 2001 From: Knut Andreas Date: Tue, 4 Nov 2025 20:19:14 +0100 Subject: [PATCH 14/14] Try adding an ignore to OT --- .typos.toml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.typos.toml b/.typos.toml index 8a8d8b3..ac18a55 100644 --- a/.typos.toml +++ b/.typos.toml @@ -1,2 +1,5 @@ [files] extend-exclude = ["*.toml"] + +[default.extend-words] +OT = "OT" \ No newline at end of file