diff --git a/Project.toml b/Project.toml index eb14fe8..034afdf 100644 --- a/Project.toml +++ b/Project.toml @@ -21,10 +21,12 @@ Newton = "0.2" FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +MaterialModelsTesting = "882b014b-b96c-4115-8629-e17fb35110d2" [targets] -test = ["Test", "ForwardDiff", "FiniteDiff"] +test = ["Test", "ForwardDiff", "FiniteDiff", "MaterialModelsTesting"] [sources] -MaterialModelsBase = {rev = "main", url = "https://github.com/knutam/MaterialModelsBase.jl"} -Newton = {rev = "main", url = "https://github.com/knutam/Newton.jl"} +MaterialModelsBase = {url = "https://github.com/knutam/MaterialModelsBase.jl"} +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 ce8cca2..6c5ef2e 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -8,5 +8,5 @@ Newton = "83aa5b51-0588-403c-85e4-434ec185aae7" Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4" [sources] -MaterialModelsBase = {rev = "main", url = "https://github.com/knutam/MaterialModelsBase.jl"} -Newton = {rev = "main", url = "https://github.com/knutam/Newton.jl"} +MaterialModelsBase = {url = "https://github.com/knutam/MaterialModelsBase.jl"} +Newton = {url = "https://github.com/knutam/Newton.jl"} diff --git a/docs/src/devdocs.md b/docs/src/devdocs.md index 7d48817..2857c75 100644 --- a/docs/src/devdocs.md +++ b/docs/src/devdocs.md @@ -15,6 +15,4 @@ compute_stress compute_tangent compute_stress_and_tangent static_vector -get_resid_eltype -get_num_unknowns ``` \ No newline at end of file diff --git a/src/CrystalPlasticity.jl b/src/CrystalPlasticity.jl index 045a480..e9ec0c9 100644 --- a/src/CrystalPlasticity.jl +++ b/src/CrystalPlasticity.jl @@ -20,7 +20,7 @@ The parameters for this model are, See [Meyer (2020)](https://doi.org/10.1016/j.ijsolstr.2020.04.037) for a description of this model in a finite strain framework. """ -@kwdef struct CrystalPlasticity{C, E, T, OS} <: AbstractMaterial +@kwdef struct CrystalPlasticity{C, E, T, OS, T_tol} <: AbstractMaterial crystal::C # E.g. BCC, FCC, etc. elastic::E # Elastic definition yield::T # Initial yield limit @@ -33,6 +33,8 @@ in a finite strain framework. Hkin::T # Kinematic hardening modulus β∞::T # Kinematic saturation stress overstress::OS # Overstress function + maxiter::Int = 100 + tolerance::T_tol = sqrt(eps(typeof(q))) end struct CrystalPlasticityState{T, nslip} <: MMB.AbstractMaterialState @@ -96,7 +98,7 @@ function MMB.material_response(mat::CrystalPlasticity, ϵ::SymmetricTensor{2,3}, Δγ0 = zero(typeof(τ_trial_red)) ν_trial = sign.(τ_trial_red) rf(x) = residual(mat, x, old, ϵ, ν_trial, Δt) - Δγ, ∂r∂x, converged = newtonsolve(rf, Δγ0) + Δγ, ∂r∂x, converged = newtonsolve(rf, Δγ0; maxiter = mat.maxiter, tol = mat.tolerance) if converged σ = calculate_stress(mat, Δγ, old, ϵ, ν_trial) # R(Δγ(ϵ), ϵ) = 0 => ∂R/∂x ∂γ/∂ϵ = - ∂R/∂ϵ diff --git a/src/Elastic.jl b/src/Elastic.jl index 541ea87..4f232d5 100644 --- a/src/Elastic.jl +++ b/src/Elastic.jl @@ -55,7 +55,7 @@ struct LinearElastic{T, case, N} <: AbstractMaterial p::SVector{N,T} end -MMB.get_params_eltype(::LinearElastic{T}) where T = T +MMB.get_vector_eltype(::LinearElastic{T}) where T = T # General symmetry LinearElastic{:general}(C::SymmetricTensor) = LinearElastic(C) @@ -94,7 +94,7 @@ end calculate_stress(m::LinearElastic, ϵ::SymmetricTensor) = m.C⊡ϵ # Functions for conversion between material and parameter vectors -MMB.get_num_params(::LinearElastic{<:Any,<:Any,N}) where{N} = N +MMB.get_vector_length(::LinearElastic{<:Any,<:Any,N}) where{N} = N function MMB.fromvector(v::AbstractVector, ::LinearElastic{<:Any, :isotropic}; offset=0) return LinearElastic{:isotropic}(;E=v[offset+1], ν=v[offset+2]) @@ -115,12 +115,10 @@ function MMB.differentiate_material!(deriv::MaterialDerivatives{T}, m::LinearEla tomandel!(deriv.dσdϵ, dσdϵ) p = tovector(m) σ_from_param(p_vector) = tomandel(calculate_stress(fromvector(p_vector, m), ϵ)) - ForwardDiff.jacobian!(deriv.dσdp,σ_from_param,p) + ForwardDiff.jacobian!(deriv.dσdp, σ_from_param, p) - fill!(deriv.dσdⁿs, zero(T)) - fill!(deriv.dsdϵ, zero(T)) - fill!(deriv.dsdp, zero(T)) - fill!(deriv.dsdⁿs, zero(T)) + fill!(deriv.dsdϵ, zero(T)) + fill!(deriv.dsdp, zero(T)) end # Show methods diff --git a/src/ExtraOutputs.jl b/src/ExtraOutputs.jl index 3be955b..1747be5 100644 --- a/src/ExtraOutputs.jl +++ b/src/ExtraOutputs.jl @@ -10,25 +10,27 @@ mutable struct DiffOutputHelper{TX, T} <: AbstractExtraOutput const ∂R∂ⁿsᴹ::Matrix{T} const ∂R∂pᴹ::Matrix{T} const ∂s∂Xᴹ::Matrix{T} - const ∂s∂ⁿsᴹ::Matrix{T} const ∂s∂pᴹ::Matrix{T} + const dσdⁿsᴹ::Matrix{T} + const dsdⁿsᴹ::Matrix{T} end function DiffOutputHelper(m::AbstractMaterial) - T = MMB.get_params_eltype(m) + T = MMB.get_vector_eltype(m) X = initial_guess(m, initial_material_state(m), zero(SymmetricTensor{2,3})) - NR = get_num_unknowns(X) + NR = get_vector_length(X) Nσ = get_num_tensorcomponents(m) Ns = get_num_statevars(m) - Np = get_num_params(m) + Np = get_vector_length(m) dRdX_invᴹ = zeros(T, NR, NR) ∂R∂ϵᴹ = zeros(T, NR, Nσ) ∂R∂ⁿsᴹ = zeros(T, NR, Ns) ∂R∂pᴹ = zeros(T, NR, Np) ∂s∂Xᴹ = zeros(T, Ns, NR) - ∂s∂ⁿsᴹ = zeros(T, Ns, Ns) ∂s∂pᴹ = zeros(T, Ns, Np) - return DiffOutputHelper(X, dRdX_invᴹ, false, ∂R∂ϵᴹ, ∂R∂ⁿsᴹ, ∂R∂pᴹ, ∂s∂Xᴹ, ∂s∂ⁿsᴹ, ∂s∂pᴹ) + dσdⁿsᴹ = zeros(T, Nσ, Ns) + dsdⁿsᴹ = zeros(T, Ns, Ns) + return DiffOutputHelper(X, dRdX_invᴹ, false, ∂R∂ϵᴹ, ∂R∂ⁿsᴹ, ∂R∂pᴹ, ∂s∂Xᴹ, ∂s∂pᴹ, dσdⁿsᴹ, dsdⁿsᴹ) end update_extras!(extras::DiffOutputHelper) = (extras.updated = false) @@ -37,4 +39,4 @@ function update_extras!(extras::DiffOutputHelper, X, dRdX_invᴹ) extras.X = X extras.dRdX_invᴹ .= dRdX_invᴹ extras.updated = true -end \ No newline at end of file +end diff --git a/src/FiniteStrainPlastic.jl b/src/FiniteStrainPlastic.jl index f30065a..4aa0f6f 100644 --- a/src/FiniteStrainPlastic.jl +++ b/src/FiniteStrainPlastic.jl @@ -28,24 +28,24 @@ m = FiniteStrainPlastic(elastic = CompressibleNeoHooke(G=80.e3, K=160.e3), ``` """ -struct FiniteStrainPlastic{E,YLD,IsoH,KinH,OS} <: AbstractMaterial +struct FiniteStrainPlastic{E,YLD,IsoH,KinH,OS,T} <: AbstractMaterial elastic::E # Elastic definition yield::YLD # Initial yield limit isotropic::IsoH # Tuple of isotropic hardening definitions kinematic::KinH # Tuple of kinematic hardening definitions overstress::OS # Overstress function + maxiter::Int # Maximum number of iterations for local problem (default 10) + tolerance::T # Absolute tolerance for local problem (default √eps(T)) end -function FiniteStrainPlastic(;elastic::AbstractHyperElastic, yield, isotropic=nothing, kinematic=nothing, overstress=RateIndependent()) +function FiniteStrainPlastic(;elastic::AbstractHyperElastic, yield, isotropic=nothing, kinematic=nothing, overstress=RateIndependent(), maxiter = 100, tolerance = sqrt(eps(MMB.get_vector_eltype(elastic)))) iso = maketuple_or_nothing(isotropic) kin = maketuple_or_nothing(kinematic) yld = default_yield_criteria(yield) - return FiniteStrainPlastic(elastic, yld, iso, kin, overstress) + return FiniteStrainPlastic(elastic, yld, iso, kin, overstress, maxiter, tolerance) end - -MMB.get_params_eltype(m::FiniteStrainPlastic) = typeof(initial_yield_limit(m.yield)) - -MMB.get_tensorbase(::FiniteStrainPlastic) = Tensor{2,3} +MMB.get_vector_eltype(m::FiniteStrainPlastic) = typeof(initial_yield_limit(m.yield)) +MMB.get_tensorbase(::FiniteStrainPlastic) = Tensor{2, 3} # Definition of material state struct FiniteStrainPlasticState{NKin,NIso,TFp,Tκ<:NTuple{NIso},TFk<:NTuple{NKin}} <: AbstractMaterialState @@ -55,7 +55,7 @@ struct FiniteStrainPlasticState{NKin,NIso,TFp,Tκ<:NTuple{NIso},TFk<:NTuple{NKin end function MMB.initial_material_state(m::FiniteStrainPlastic) - T = MMB.get_params_eltype(m) + T = MMB.get_vector_eltype(m) I2 = one(Tensor{2,3,T}) FiniteStrainPlasticState( I2, # Fp @@ -78,7 +78,7 @@ function FiniteStrainPlasticResidual(Mred, Δλ, κ::Tuple, Mk::Tuple) return FiniteStrainPlasticResidual(Mred, Δλ, map(x->convert(Tκ, x), κ), map(x->convert(TMk, x), Mk)) end -function get_num_unknowns(::FiniteStrainPlasticResidual{NKin_m1,NIso}) where{NKin_m1,NIso} +function MMB.get_vector_length(::FiniteStrainPlasticResidual{NKin_m1,NIso}) where{NKin_m1,NIso} return 10 + NIso + 9*NKin_m1 end @@ -98,22 +98,12 @@ function MMB.fromvector(v::AbstractVector{T}, ::FiniteStrainPlasticResidual{NKin return FiniteStrainPlasticResidual(Mred, Δλ, κ, Mk) end -function MMB.allocate_material_cache(m::FiniteStrainPlastic) - T = MMB.get_params_eltype(m) - s = initial_material_state(m) - F = one(Tensor{2,3}) - x = initial_guess(m, s, F) - xv = Vector{T}(undef, get_num_unknowns(x)) - rf!(r_vector, x_vector) = vector_residual!((x)->residual(x, m, old, F, zero(T)), r_vector, x_vector, x) - return NewtonCache(xv) -end - function mandel_stress(m_el::AbstractHyperElastic, Fe::Tensor{2,3}) Ce = tdot(Fe) return Ce ⋅ compute_stress(m_el, Ce) end -function MMB.material_response(m::FiniteStrainPlastic, F::Tensor{2,3}, old::FiniteStrainPlasticState, Δt, cache, extras) +function MMB.material_response(m::FiniteStrainPlastic, F::Tensor{2,3}, old::FiniteStrainPlasticState, Δt, _, extras) Fpinv = inv(old.Fp) Fe = F ⋅ Fpinv M = mandel_stress(m.elastic, Fe) @@ -127,7 +117,7 @@ function MMB.material_response(m::FiniteStrainPlastic, F::Tensor{2,3}, old::Fini else x0 = initial_guess(m, old, M) rf(x) = tovector(SVector, residual(fromvector(x, x0), m, old, F, Δt)) - x_vector, ∂R∂X, converged = newtonsolve(rf, tovector(SVector, x0)) + x_vector, ∂R∂X, converged = newtonsolve(rf, tovector(SVector, x0); maxiter = m.maxiter, tol = m.tolerance) if converged x_sol = fromvector(x_vector, x0) @@ -152,7 +142,7 @@ function MMB.material_response(m::FiniteStrainPlastic, F::Tensor{2,3}, old::Fini return P, dPdF, new else - throw(MMB.NoLocalConvergence("$(typeof(m)): newtonsolve! didn't converge, F = ", F)) + throw(MMB.NoLocalConvergence("$(typeof(m)): newtonsolve didn't converge, F = ", F)) end end end @@ -243,20 +233,20 @@ end #= # Functions for conversion between material and parameter vectors -function MMB.get_num_params(m::Plastic) +function MMB.get_vector_length(m::Plastic) return sum(( - get_num_params(m.elastic), - get_num_params(m.yield), + get_vector_length(m.elastic), + get_vector_length(m.yield), sum(get_num_params, m.isotropic), sum(get_num_params, m.kinematic), - get_num_params(m.overstress))) + get_vector_length(m.overstress))) end function fromvectortuple(v::AbstractVector, materialtuple; offset=0) n = 0 mout = map(materialtuple) do mt m = fromvector(v, mt; offset=(offset+n)) - n += get_num_params(m) + n += get_vector_length(m) m end return mout, n @@ -264,28 +254,28 @@ end function MMB.fromvector(v::AbstractVector, m::Plastic; offset=0) i = offset - elastic = fromvector(v, m.elastic, offset=i); i += get_num_params(m.elastic) - yield = fromvector(v, m.yield, offset=i); i += get_num_params(m.yield) + elastic = fromvector(v, m.elastic, offset=i); i += get_vector_length(m.elastic) + yield = fromvector(v, m.yield, offset=i); i += get_vector_length(m.yield) isotropic, nisoparam = fromvectortuple(v, m.isotropic, offset=i); i+=nisoparam kinematic, nkinparam = fromvectortuple(v, m.kinematic, offset=i); i+=nkinparam - overstress = fromvector(v, m.overstress, offset=i); # i += get_num_params(m.overstress) + overstress = fromvector(v, m.overstress, offset=i); # i += get_vector_length(m.overstress) return Plastic(;elastic, yield, isotropic, kinematic, overstress) end function MMB.tovector!(v::AbstractVector, m::Plastic; offset=0) i = offset - tovector!(v, m.elastic; offset=i); i += get_num_params(m.elastic) - tovector!(v, m.yield; offset=i); i += get_num_params(m.yield) + tovector!(v, m.elastic; offset=i); i += get_vector_length(m.elastic) + tovector!(v, m.yield; offset=i); i += get_vector_length(m.yield) for iso in m.isotropic - tovector!(v, iso;offset=i); i+= get_num_params(iso) + tovector!(v, iso;offset=i); i+= get_vector_length(iso) end for kin in m.kinematic - tovector!(v, kin; offset=i); i+= get_num_params(kin) + tovector!(v, kin; offset=i); i+= get_vector_length(kin) end - tovector!(v, m.overstress; offset=i); # i += get_num_params(m.overstress) + tovector!(v, m.overstress; offset=i); # i += get_vector_length(m.overstress) return v end diff --git a/src/Plastic.jl b/src/Plastic.jl index 7424bc7..abdd954 100644 --- a/src/Plastic.jl +++ b/src/Plastic.jl @@ -28,21 +28,23 @@ m = Plastic(elastic = LinearElastic(E=210.e3, ν=0.3), ``` """ -struct Plastic{E,YLD,IsoH,KinH,OS} <: AbstractMaterial +struct Plastic{E,YLD,IsoH,KinH,OS,T} <: AbstractMaterial elastic::E # Elastic definition yield::YLD # Initial yield limit isotropic::IsoH # Tuple of isotropic hardening definitions kinematic::KinH # Tuple of kinematic hardening definitions overstress::OS # Overstress function + maxiter::Int # Maximum number of iterations for local problem (default 100) + tolerance::T # Absolute tolerance for local problem (default √eps(T)) end -function Plastic(;elastic, yield, isotropic=nothing, kinematic=nothing, overstress=RateIndependent()) +function Plastic(;elastic, yield, isotropic=nothing, kinematic=nothing, overstress=RateIndependent(), maxiter = 100, tolerance = sqrt(eps(MMB.get_vector_eltype(elastic)))) iso = maketuple_or_nothing(isotropic) kin = maketuple_or_nothing(kinematic) yld = default_yield_criteria(yield) - return Plastic(elastic, yld, iso, kin, overstress) + return Plastic(elastic, yld, iso, kin, overstress, maxiter, tolerance) end -MMB.get_params_eltype(m::Plastic) = MMB.get_params_eltype(m.elastic) +MMB.get_vector_eltype(m::Plastic) = MMB.get_vector_eltype(m.elastic) # Definition of material state struct PlasticState{NKin,NIso,Tϵp,Tκ,Tβ} <: AbstractMaterialState @@ -53,15 +55,15 @@ struct PlasticState{NKin,NIso,Tϵp,Tκ,Tβ} <: AbstractMaterialState end function MMB.initial_material_state(m::Plastic) - T = MMB.get_params_eltype(m) + T = MMB.get_vector_eltype(m) PlasticState(zero(SymmetricTensor{2,3,T}), ntuple(i->get_initial_value(m.isotropic[i]), length(m.isotropic)), ntuple(i->zero(SymmetricTensor{2,3,T}), length(m.kinematic))) end -MMB.get_num_statevars(::PlasticState{NK, NI}) where {NK, NI} = NI + 6*(1+NK) +MMB.get_vector_length(::PlasticState{NK, NI}) where {NK, NI} = NI + 6*(1+NK) MMB.get_num_statevars(m::Plastic) = length(m.isotropic) + 6*(1 + length(m.kinematic)) -function MMB.get_statevar_eltype(::PlasticState{<:Any, <:Any, Tϵp, Tκ, Tβ}) where {Tϵp, Tκ, Tβ} +function MMB.get_vector_eltype(::PlasticState{<:Any, <:Any, Tϵp, Tκ, Tβ}) where {Tϵp, Tκ, Tβ} return promote_type(Tϵp, Tκ, Tβ) end @@ -89,11 +91,11 @@ struct PlasticResidual{NKin,NIso,Tσ,Tλ,Tκ,Tβ} <: AbstractResidual β::NTuple{NKin, SymmetricTensor{2,3,Tβ,6}} end -function get_resid_eltype(::PlasticResidual{<:Any, <:Any, Tσ, Tλ, Tκ, Tβ}) where {Tσ, Tλ, Tκ, Tβ} +function MMB.get_vector_eltype(::PlasticResidual{<:Any, <:Any, Tσ, Tλ, Tκ, Tβ}) where {Tσ, Tλ, Tκ, Tβ} return promote_type(Tσ, Tλ, Tκ, Tβ) end -get_num_unknowns(::PlasticResidual{NKin,NIso}) where {NKin, NIso} = 7 + NIso + 6 * NKin +MMB.get_vector_length(::PlasticResidual{NKin,NIso}) where {NKin, NIso} = 7 + NIso + 6 * NKin function MMB.tovector!(v::AbstractVector, r::PlasticResidual{NKin, NIso}; offset = 0) where {NKin, NIso} tomandel!(v, r.σ; offset=offset) @@ -113,17 +115,17 @@ function MMB.fromvector(v::AbstractVector, ::PlasticResidual{NKin,NIso}; offset= return PlasticResidual(σ, Δλ, κ, β) end -struct PlasticCache{TN,TR} +struct PlasticCache{TN,TR} <: AbstractMaterialCache newton::TN # NewtonCache if needed, otherwise nothing resid::TR # Cache used in residual function if needed, otherwise nothing end function get_newton_cache(m::Plastic, residual_cache) - T = MMB.get_params_eltype(m) + T = MMB.get_vector_eltype(m) s = initial_material_state(m) ϵ = zero(SymmetricTensor{2,3}) x = initial_guess(m, s, ϵ) - xv = Vector{T}(undef, get_num_unknowns(x)) + xv = Vector{T}(undef, get_vector_length(x)) rf!(r_vector, x_vector) = vector_residual!((x)->residual(x, m, old, ϵ, zero(T), residual_cache), r_vector, x_vector, x) return NewtonCache(xv) end @@ -142,7 +144,7 @@ function MMB.material_response(m::Plastic, ϵ::SymmetricTensor{2,3}, old::Plasti σ_trial, dσdϵ_elastic = elastic_response(m, ϵ, old) Φ_trial = yield_criterion(m.yield, σ_trial-sum(old.β), sum(old.κ)) - if Φ_trial < 0 + if real(Φ_trial) < 0 update_extras!(extras) return σ_trial, dσdϵ_elastic, old else @@ -150,7 +152,7 @@ function MMB.material_response(m::Plastic, ϵ::SymmetricTensor{2,3}, old::Plasti rf!(r_vector, x_vector) = vector_residual!(xx->residual(xx, m, old, ϵ, Δt, cache.resid), r_vector, x_vector, x) x_vector = getx(cache.newton) tovector!(x_vector, x) - x_vector, dRdx, converged = newtonsolve(rf!, x_vector, cache.newton) + x_vector, dRdx, converged = newtonsolve(rf!, x_vector, cache.newton; tol = m.tolerance, maxiter = m.maxiter) if converged x_sol = fromvector(x_vector, x) check_solution(x_sol) @@ -206,20 +208,20 @@ function calculate_elastic_strain(old::PlasticState, ϵ, ν, Δλ) end # Functions for conversion between material and parameter vectors -function MMB.get_num_params(m::Plastic) +function MMB.get_vector_length(m::Plastic) return sum(( - get_num_params(m.elastic), - get_num_params(m.yield), - sum(get_num_params, m.isotropic), - sum(get_num_params, m.kinematic), - get_num_params(m.overstress))) + get_vector_length(m.elastic), + get_vector_length(m.yield), + sum(get_vector_length, m.isotropic), + sum(get_vector_length, m.kinematic), + get_vector_length(m.overstress))) end function fromvectortuple(v::AbstractVector, materialtuple; offset=0) n = Ref(0) mout = map(materialtuple) do mt m = fromvector(v, mt; offset=(offset+n[])) - n[] += get_num_params(mt) + n[] += get_vector_length(mt) m end return mout, n[] @@ -227,32 +229,32 @@ end function MMB.fromvector(v::AbstractVector, m::Plastic; offset=0) i = offset - elastic = fromvector(v, m.elastic, offset=i); i += get_num_params(m.elastic) - yield = fromvector(v, m.yield, offset=i); i += get_num_params(m.yield) + elastic = fromvector(v, m.elastic, offset=i); i += get_vector_length(m.elastic) + yield = fromvector(v, m.yield, offset=i); i += get_vector_length(m.yield) isotropic, nisoparam = fromvectortuple(v, m.isotropic, offset=i) i += nisoparam kinematic, nkinparam = fromvectortuple(v, m.kinematic, offset=i) i += nkinparam - overstress = fromvector(v, m.overstress, offset=i); # i += get_num_params(m.overstress) + overstress = fromvector(v, m.overstress, offset=i); # i += get_vector_length(m.overstress) # The following constructor call doesn't seem to be type stable when used in e.g. differentiate_material. # Should be checked for (a) v and m same eltype and (b) v Dual and m Float64 eltypes. #return Plastic(;elastic, yield, isotropic, kinematic, overstress) - return Plastic(elastic, yield, isotropic, kinematic, overstress) + return Plastic(elastic, yield, isotropic, kinematic, overstress, m.maxiter, m.tolerance) end function MMB.tovector!(v::AbstractVector, m::Plastic; offset=0) i = offset - tovector!(v, m.elastic; offset=i); i += get_num_params(m.elastic) - tovector!(v, m.yield; offset=i); i += get_num_params(m.yield) + tovector!(v, m.elastic; offset=i); i += get_vector_length(m.elastic) + tovector!(v, m.yield; offset=i); i += get_vector_length(m.yield) for iso in m.isotropic - tovector!(v, iso;offset=i); i+= get_num_params(iso) + tovector!(v, iso;offset=i); i+= get_vector_length(iso) end for kin in m.kinematic - tovector!(v, kin; offset=i); i+= get_num_params(kin) + tovector!(v, kin; offset=i); i+= get_vector_length(kin) end - tovector!(v, m.overstress; offset=i); # i += get_num_params(m.overstress) + tovector!(v, m.overstress; offset=i); # i += get_vector_length(m.overstress) return v end diff --git a/src/PlasticDifferentiate.jl b/src/PlasticDifferentiate.jl index facc2d4..87528ae 100644 --- a/src/PlasticDifferentiate.jl +++ b/src/PlasticDifferentiate.jl @@ -1,15 +1,16 @@ MMB.allocate_differentiation_output(m::Plastic) = DiffOutputHelper(m) + function MMB.differentiate_material!(deriv::MaterialDerivatives, m::Plastic, ϵ, ⁿs, Δt, cache, diff_helper::DiffOutputHelper, dσdϵ) if diff_helper.updated # Plastic response differentiate_material_plastic!(deriv, m, ϵ, ⁿs, Δt, cache, diff_helper, dσdϵ) else # Elastic response - differentiate_material_elastic!(deriv, m, ϵ, ⁿs, dσdϵ) + differentiate_material_elastic!(deriv, m, ϵ, ⁿs, diff_helper, dσdϵ) end end -function differentiate_material_elastic!(deriv::MaterialDerivatives{T}, m, ϵ, ⁿs, dσdϵ) where{T} +function differentiate_material_elastic!(deriv::MaterialDerivatives{T}, m, ϵ, ⁿs, dh::DiffOutputHelper, dσdϵ) where{T} p = tovector(m) - sv = zeros(get_num_statevars(ⁿs)) + sv = tovector(ⁿs) # Differentiate stress tomandel!(deriv.dσdϵ, dσdϵ) # Already calculated @@ -18,22 +19,22 @@ function differentiate_material_elastic!(deriv::MaterialDerivatives{T}, m, ϵ, ## Differentiate stress by old state σ_from_state(old_vector) = tomandel(calculate_elastic_stress(m, ϵ, fromvector(old_vector, ⁿs))) - ForwardDiff.jacobian!(deriv.dσdⁿs, σ_from_state, tovector!(sv, ⁿs)) + ForwardDiff.jacobian!(dh.dσdⁿsᴹ, σ_from_state, sv) ## Differentiate stress by parameters σ_from_param(p_vector) = tomandel(calculate_elastic_stress(fromvector(p_vector, m), ϵ, ⁿs)) ForwardDiff.jacobian!(deriv.dσdp, σ_from_param, p) + deriv.dσdp .+= dh.dσdⁿsᴹ * deriv.dsdp # Differentiate state - fill!(deriv.dsdϵ, zero(T)) # Constant state - # new state equal to old state => derivative is the identity matrix - copyto!(deriv.dsdⁿs, LinearAlgebra.I) - fill!(deriv.dsdp, zero(T)) # Constant state + fill!(deriv.dsdϵ, zero(T)) # Constant state + # deriv.dsdp # Constant state => s depends equally on p as it did in the last time step end function differentiate_material_plastic!(deriv::MaterialDerivatives, m::Plastic, ϵ, ⁿs, Δt, cache, diff_helper::DiffOutputHelper, dσdϵ) # Extract from input p = tovector(m) + sv = tovector(ⁿs) X = diff_helper.X ∂R∂Xinvᴹ = diff_helper.dRdX_invᴹ @@ -42,17 +43,16 @@ function differentiate_material_plastic!(deriv::MaterialDerivatives, m::Plastic, ∂R∂ⁿsᴹ = diff_helper.∂R∂ⁿsᴹ ∂R∂pᴹ = diff_helper.∂R∂pᴹ ∂s∂Xᴹ = diff_helper.∂s∂Xᴹ - ∂s∂ⁿsᴹ = diff_helper.∂s∂ⁿsᴹ ∂s∂pᴹ = diff_helper.∂s∂pᴹ + dσdⁿsᴹ = diff_helper.dσdⁿsᴹ + dsdⁿsᴹ = diff_helper.dsdⁿsᴹ # Precalculations - ⁿs_vector = tovector(ⁿs) - R_from_strain(ϵ_vector) = tovector(residual(X, m, ⁿs, frommandel(baseof(ϵ), ϵ_vector), Δt, cache.resid)) ForwardDiff.jacobian!(∂R∂ϵᴹ, R_from_strain, tomandel(ϵ)) R_from_state(old_vector) = tovector(residual(X, m, fromvector(old_vector, ⁿs), ϵ, Δt, cache.resid)) - ForwardDiff.jacobian!(∂R∂ⁿsᴹ, R_from_state, ⁿs_vector) + ForwardDiff.jacobian!(∂R∂ⁿsᴹ, R_from_state, sv) R_from_param(p_vector) = tovector(residual(X, fromvector(p_vector, m), ⁿs, ϵ, Δt, cache.resid)) ForwardDiff.jacobian!(∂R∂pᴹ, R_from_param, p) @@ -60,27 +60,30 @@ function differentiate_material_plastic!(deriv::MaterialDerivatives, m::Plastic, # Differentiate the stress ## dσdXᴹ = [I 0] where I is 6x6 identity matrix ∂σ∂Xᴹ_times_∂R∂Xinvᴹ = ∂R∂Xinvᴹ[1:Nσ, :] - ## dσdϵ + ## dσdϵ tomandel!(deriv.dσdϵ, dσdϵ) # Already calculated - ## dσdⁿs - deriv.dσdⁿs .= -∂σ∂Xᴹ_times_∂R∂Xinvᴹ * ∂R∂ⁿsᴹ # ∂σ∂ⁿs=0 (σ=X.σ) + ## dσdⁿs + dσdⁿsᴹ .= -∂σ∂Xᴹ_times_∂R∂Xinvᴹ * ∂R∂ⁿsᴹ # ∂σ∂ⁿs=0 (σ=X.σ) ## dσdp deriv.dσdp .= -∂σ∂Xᴹ_times_∂R∂Xinvᴹ * ∂R∂pᴹ # ∂σ∂p=0 (σ=X.σ) + deriv.dσdp .+= dσdⁿsᴹ * deriv.dsdp # Differentiate the new state ## Calculate ∂s∂Xᴹ first s_from_X(X_vector) = tovector(get_plastic_result(m, fromvector(X_vector, X), ⁿs)[2]) ForwardDiff.jacobian!(∂s∂Xᴹ, s_from_X, tovector(X)) ∂s∂Xᴹ_times_∂R∂Xinvᴹ = ∂s∂Xᴹ*∂R∂Xinvᴹ + ## dsdϵ deriv.dsdϵ .= -∂s∂Xᴹ_times_∂R∂Xinvᴹ*∂R∂ϵᴹ # ∂s∂ϵ = 0 (get_plastic_result does not take the strain as argument) + ## dsdⁿs s_from_state(old_vector) = tovector(get_plastic_result(m, X, fromvector(old_vector, ⁿs))[2]) - ForwardDiff.jacobian!(∂s∂ⁿsᴹ, s_from_state, ⁿs_vector) - deriv.dsdⁿs .= -∂s∂Xᴹ_times_∂R∂Xinvᴹ*∂R∂ⁿsᴹ + ∂s∂ⁿsᴹ + ForwardDiff.jacobian!(dsdⁿsᴹ, s_from_state, sv) + dsdⁿsᴹ .-= ∂s∂Xᴹ_times_∂R∂Xinvᴹ*∂R∂ⁿsᴹ + ## dsdp s_from_p(p_vector) = tovector(get_plastic_result(fromvector(p_vector, m), X, ⁿs)[2]) ForwardDiff.jacobian!(∂s∂pᴹ, s_from_p, p) - deriv.dsdp .= -∂s∂Xᴹ_times_∂R∂Xinvᴹ*∂R∂pᴹ + ∂s∂pᴹ - -end \ No newline at end of file + deriv.dsdp .= -∂s∂Xᴹ_times_∂R∂Xinvᴹ*∂R∂pᴹ .+ ∂s∂pᴹ .+ dsdⁿsᴹ * deriv.dsdp +end diff --git a/src/ViscoElastic.jl b/src/ViscoElastic.jl index 2374c4e..bce9b2b 100644 --- a/src/ViscoElastic.jl +++ b/src/ViscoElastic.jl @@ -50,14 +50,14 @@ function GeneralizedMaxwell(base::LinearElastic, chains::Maxwell...) return GeneralizedMaxwell(base, chains) end -MMB.get_params_eltype(::GeneralizedMaxwell{<:Any, T}) where {T} = T +MMB.get_vector_eltype(::GeneralizedMaxwell{<:Any, T}) where {T} = T struct GeneralizedMaxwellState{T, num} <: AbstractMaterialState ϵv::NTuple{num, SymmetricTensor{2,3,T,6}} end function MMB.initial_material_state(m::GeneralizedMaxwell{<:Any, <:Any, num}) where {num} - T = MMB.get_params_eltype(m) + T = MMB.get_vector_eltype(m) return GeneralizedMaxwellState(ntuple(_ -> zero(SymmetricTensor{2, 3, T}), num)) end diff --git a/src/hyper_elasticity/HyperElastic.jl b/src/hyper_elasticity/HyperElastic.jl index cab4ea2..ef69b7e 100644 --- a/src/hyper_elasticity/HyperElastic.jl +++ b/src/hyper_elasticity/HyperElastic.jl @@ -36,13 +36,13 @@ function compute_stress_and_tangent(m::AbstractHyperElastic, C::SymmetricTensor) return S, 2*∂S∂C end -function MMB.material_response(m::AbstractHyperElastic, F::Tensor{2,3}, args...) +function MMB.material_response(m::AbstractHyperElastic, F::Tensor{2,3}, old::AbstractMaterialState, args...) C = tdot(F) S, ∂S∂E = compute_stress_and_tangent(m, C) P = F ⋅ S I2 = one(F) ∂P∂F = F ⋅ ∂S∂E ⊡ otimesu(F',I2) + otimesu(I2, S) - return P, ∂P∂F, NoMaterialState() + return P, ∂P∂F, old end MMB.get_tensorbase(::AbstractHyperElastic) = Tensor{2,3} diff --git a/src/hyper_elasticity/NeoHooke.jl b/src/hyper_elasticity/NeoHooke.jl index 0633bba..20618f8 100644 --- a/src/hyper_elasticity/NeoHooke.jl +++ b/src/hyper_elasticity/NeoHooke.jl @@ -16,6 +16,7 @@ which has no influence if `C` truly is incompressible, i.e. ``\\det(\\boldsymbol end compute_potential(m::NeoHooke, C::SymmetricTensor) = (m.G / 2) * (tr(C) / cbrt(det(C)) - 3) +MMB.get_vector_eltype(::NeoHooke{T}) where {T} = T """ CompressibleNeoHooke(; G, K) @@ -42,3 +43,5 @@ function compute_potential(m::CompressibleNeoHooke, C::SymmetricTensor) ΨK = (m.K / 2) * (sqrt(detC) - 1)^2 return ΨG + ΨK end + +MMB.get_vector_eltype(::CompressibleNeoHooke{T}) where {T} = T diff --git a/src/plasticity_components/IsotropicHardening.jl b/src/plasticity_components/IsotropicHardening.jl index 830082d..9f31d34 100644 --- a/src/plasticity_components/IsotropicHardening.jl +++ b/src/plasticity_components/IsotropicHardening.jl @@ -1,5 +1,6 @@ # Isotropic hardening abstract type AbstractIsotropicHardening{T} end +MMB.get_vector_eltype(::AbstractIsotropicHardening{T}) where {T} = T # Default initial value for κ get_initial_value(::AbstractIsotropicHardening{T}) where {T} = zero(T) @@ -59,8 +60,8 @@ get_initial_value(param::Swift) = param.K*(param.λ0^param.n) ## Conversion from vector to material -MMB.get_num_params(::Voce) = 2 -MMB.get_num_params(::Swift) = 3 +MMB.get_vector_length(::Voce) = 2 +MMB.get_vector_length(::Swift) = 3 MMB.fromvector(v::AbstractVector, ::Voce; offset=0) = Voce(Hiso=v[offset+1], κ∞=v[offset+2]) MMB.fromvector(v::AbstractVector, ::Swift; offset=0) = Swift(K=v[offset+1], λ0=v[offset+2], n=v[offset+3]) diff --git a/src/plasticity_components/KinematicHardening.jl b/src/plasticity_components/KinematicHardening.jl index 4cc6549..db4f866 100644 --- a/src/plasticity_components/KinematicHardening.jl +++ b/src/plasticity_components/KinematicHardening.jl @@ -1,5 +1,6 @@ # Kinematic hardening abstract type AbstractKinematicHardening{T} end +MMB.get_vector_eltype(::AbstractKinematicHardening{T}) where {T} = T """ ArmstrongFrederick(; Hkin, β∞) @@ -65,8 +66,8 @@ function get_evolution(param::OhnoWang{Tp}, ν::SecondOrderTensor, βᵢ::Second end ## Conversion from vector to material -MMB.get_num_params(::ArmstrongFrederick) = 2 -MMB.get_num_params(::Union{Delobelle,OhnoWang}) = 3 +MMB.get_vector_length(::ArmstrongFrederick) = 2 +MMB.get_vector_length(::Union{Delobelle,OhnoWang}) = 3 MMB.fromvector(v::AbstractVector, ::ArmstrongFrederick; offset=0) = ArmstrongFrederick(Hkin=v[offset+1], β∞=v[offset+2]) MMB.fromvector(v::AbstractVector, ::Delobelle; offset=0) = Delobelle(Hkin=v[offset+1], β∞=v[offset+2], δ=v[offset+3]) diff --git a/src/plasticity_components/OverstressFunctions.jl b/src/plasticity_components/OverstressFunctions.jl index 90925ab..d3dc36e 100644 --- a/src/plasticity_components/OverstressFunctions.jl +++ b/src/plasticity_components/OverstressFunctions.jl @@ -9,7 +9,7 @@ the so-called KKT loading/unloading conditions """ struct RateIndependent end -MMB.get_num_params(::RateIndependent) = 0 +MMB.get_vector_length(::RateIndependent) = 0 MMB.tovector!(v, ::RateIndependent; kwargs...) = nothing MMB.fromvector(v, ::RateIndependent; kwargs...) = RateIndependent() @@ -40,7 +40,8 @@ end overstress_function(ratelaw::NortonOverstress, Φ, σy) = (1/ratelaw.tstar) * macaulay(Φ/σy)^ratelaw.nexp -MMB.get_num_params(::NortonOverstress) = 2 +MMB.get_vector_length(::NortonOverstress) = 2 +MMB.get_vector_eltype(::NortonOverstress{TT, TN}) where {TT, TN} = promote_type(TT, TN) MMB.fromvector(v::AbstractVector, ::NortonOverstress; offset=0) = NortonOverstress(v[offset+1], v[offset+2]) function MMB.tovector!(v::Vector, o::NortonOverstress; offset=0) diff --git a/src/plasticity_components/YieldCriteria.jl b/src/plasticity_components/YieldCriteria.jl index 0aac4af..590e7b1 100644 --- a/src/plasticity_components/YieldCriteria.jl +++ b/src/plasticity_components/YieldCriteria.jl @@ -29,7 +29,8 @@ effective_stress(::VonMises, σred) = vonmises(σred) effective_stress_gradient(yc::VonMises, σred) = (3/2)*dev(σred)/effective_stress(yc, σred) # More efficient than using AD above -MMB.get_num_params(::VonMises) = 1 +MMB.get_vector_length(::VonMises) = 1 +MMB.get_vector_eltype(::VonMises{T}) where {T} = T MMB.fromvector(v::AbstractVector, ::VonMises; offset=0) = VonMises(v[offset+1]) function MMB.tovector!(v::Vector, yl::VonMises; offset=0) v[offset+1] = yl.Y0 @@ -60,7 +61,8 @@ function effective_stress(yc::DruckerPrager, σred::SymmetricTensor{2,3,T}) wher return sqrt(T(3)/2)*norm(dev(σred)) - yc.B*tr(σred) end -MMB.get_num_params(::DruckerPrager) = 2 +MMB.get_vector_length(::DruckerPrager) = 2 +MMB.get_vector_eltype(::DruckerPrager{T}) where {T} = T MMB.fromvector(v::AbstractVector, ::DruckerPrager; offset=0) = DruckerPrager(v[offset+1], v[offset+2]) function MMB.tovector!(v::Vector, yl::DruckerPrager; offset=0) v[offset+1] = yl.Y0 diff --git a/src/specialized_models/SimplePlastic.jl b/src/specialized_models/SimplePlastic.jl index 79c5310..b89cb7e 100644 --- a/src/specialized_models/SimplePlastic.jl +++ b/src/specialized_models/SimplePlastic.jl @@ -1,11 +1,13 @@ -@kwdef struct SimplePlastic{T} <: AbstractMaterial +@kwdef struct SimplePlastic{T, T_tol} <: AbstractMaterial G::T K::T Y0::T Hiso::T κ∞::T Hkin::T - β∞::T + β∞::T + maxiter::Int = 100 + tolerance::T_tol = sqrt(eps(typeof(G))) end function SimplePlastic(mp::Plastic{Ela, Yld, Iso, Kin, Rat}) where { @@ -20,7 +22,8 @@ function SimplePlastic(mp::Plastic{Ela, Yld, Iso, Kin, Rat}) where { return SimplePlastic(; G = E/(2 * (1 + ν)), K = E / (3 * (1 - 2 * ν)), Y0 = mp.yield.Y0, Hiso = only(mp.isotropic).Hiso, κ∞ = only(mp.isotropic).κ∞, - Hkin = only(mp.kinematic).Hkin, β∞ = only(mp.kinematic).β∞) + Hkin = only(mp.kinematic).Hkin, β∞ = only(mp.kinematic).β∞, + maxiter = mp.maxiter, tolerance = mp.tolerance) end struct SimplePlasticState{T} <: AbstractMaterialState @@ -44,7 +47,7 @@ function MMB.material_response(m::SimplePlastic, ϵ::SymmetricTensor{2,3}, old:: return σ_trial, E4, old else rf(x) = residual(x, m, ϵ, old) - Δλ, ∂r∂x, converged = newtonsolve(rf, 0.0) + Δλ, ∂r∂x, converged = newtonsolve(rf, 0.0; maxiter = m.maxiter, tol = m.tolerance) #= # Using bisection Δλ1, ∂r∂x1, converged = newtonsolve(rf, 0.0) diff --git a/src/utils.jl b/src/utils.jl index a8e2db6..f54c07d 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -58,26 +58,6 @@ baseof(::TT) where{TT} = Tensors.get_base(TT) abstract type AbstractResidual end -""" - get_num_unknowns(res::AbstractResidual) - -Return the number of unknowns for the residual `res` -""" -function get_num_unknowns end - -""" - get_resid_eltype(res::AbstractResidual) - -Return the element type used to store `res` as a vector -""" -function get_resid_eltype end - -function MMB.tovector(r::AbstractResidual) - v = zeros(get_resid_eltype(r), get_num_unknowns(r)) - MMB.tovector!(v, r) - return v -end - """ vector_residual!(rf::Function, r_vector::AbstractVector, x_vector::AbstractVector, xbase::RT) diff --git a/test/runtests.jl b/test/runtests.jl index 1cef967..a0e18f8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,15 +1,4 @@ -using Tensors, ForwardDiff, FiniteDiff, LinearAlgebra - -using MechanicalMaterialModels -import MechanicalMaterialModels as MechMat - -using MaterialModelsBase -import MaterialModelsBase as MMB - -using Test - -# Include files used for testing -include("utilities4testing.jl") +include("setup_tests.jl") # Test of utility features include("test_crystals.jl") diff --git a/test/setup_tests.jl b/test/setup_tests.jl new file mode 100644 index 0000000..63dafc1 --- /dev/null +++ b/test/setup_tests.jl @@ -0,0 +1,13 @@ +using Tensors, ForwardDiff, FiniteDiff, LinearAlgebra + +using MechanicalMaterialModels +import MechanicalMaterialModels as MechMat +import MaterialModelsTesting as MatTest + +using MaterialModelsBase +import MaterialModelsBase as MMB + +using Test + +# Include files used for testing +include("utilities4testing.jl") \ No newline at end of file diff --git a/test/test_differentiate.jl b/test/test_differentiate.jl index e33415a..fc9b7f6 100644 --- a/test/test_differentiate.jl +++ b/test/test_differentiate.jl @@ -1,5 +1,5 @@ function test_numerical_derivative(deriv, deriv_num, m, ϵ, state, Δt, cache, deriv_extra_output, dσdϵ) - obtain_numerical_material_derivative!(deriv_num, m, ϵ, state, Δt) + MatTest.obtain_numerical_material_derivative!(deriv_num, m, ϵ, state, Δt) differentiate_material!(deriv, m, ϵ, state, Δt, cache, deriv_extra_output, dσdϵ) for field in fieldnames(typeof(deriv)) @test isapprox( @@ -56,3 +56,84 @@ end end end +@testset "Numerical (3D material)" begin + E = 210.e3 + ν = 0.3 + e = LinearElastic(;E, ν) + σ_y0=150.0 + iso = (Voce(Hiso=20.e3, κ∞=200.0), Swift(K=10.0, λ0=1.e-4, n=0.4)) + kin = (ArmstrongFrederick(Hkin=5.e3, β∞=100.0), Delobelle(Hkin=3.e3, β∞=200.0, δ=0.4)) + m1 = Plastic(elastic=e, yield=σ_y0, isotropic=iso, kinematic=kin, overstress=RateIndependent()) + m2 = Plastic(elastic=e, yield=σ_y0, isotropic=iso, kinematic=kin, overstress=NortonOverstress(; tstar = 0.123, nexp = 2.5)) + materials = ["elastic" => e, "plastic" => m1, "viscoplastic" => m2] + for (key, m) in materials + @testset "$key" begin + @testset "Initial response" begin + ϵ = rand(SymmetricTensor{2,3}) * 0.1 * σ_y0 / E + Δt = 1e-2 + MatTest.test_derivative(m, ϵ, initial_material_state(m), Δt; + numdiffsettings = (fdtype = Val{:central},), + comparesettings = (atol_min = 1e-10, rtol_min = 1e-8), + ) + diff = MaterialDerivatives(m) + copy!(diff.dsdp, rand(size(diff.dsdp)...)) + MatTest.test_derivative(m, ϵ, initial_material_state(m), Δt; + numdiffsettings = (fdtype = Val{:central},), + comparesettings = (atol_min = 1e-10, rtol_min = 1e-8), + diff) + end + @testset "Just past yielding" begin + ϵ = SymmetricTensor{2,3}((i,j) -> i==j==1 ? 1.0 : (i == j ? -ν : 0.0)) * 1.01 * σ_y0 / E + Δt = 1e-2 + MatTest.test_derivative(m, ϵ, initial_material_state(m), Δt; + numdiffsettings = (fdtype = Val{:central}, relstep = 1e-8), + comparesettings = (atol_min = 1e-8, rtol_min = 1e-4, print_tol = false) + ) + end + @testset "After shear loading" begin + num_steps = 10; t_end = 0.01 + ϵ21 = num_steps * 0.1 * 3 * σ_y0 / E; + relstep = 1e-6 + stressfun(p) = MatTest.runstrain(fromvector(p, m), ϵ21, (2, 1), t_end, num_steps)[1] + dσ21_dp_num = FiniteDiff.finite_difference_jacobian(stressfun, BigFloat.(tovector(m)), Val{:central}; relstep) + σv, state, dσ21_dp, diff = MatTest.runstrain_diff(m, ϵ21, (2, 1), t_end, num_steps) + @test σv ≈ stressfun(tovector(m)) + #@test dσ21_dp ≈ dσ21_dp_num + isapprox(dσ21_dp_num, dσ21_dp; atol = 1000 * eps()) + scaled_error, maxtol = MatTest.compare_derivatives(dσ21_dp, dσ21_dp_num, σv, tovector(m) * relstep; atol_min = 1e-12, rtol_min = 1e-4, print_tol = false) + @test all(x -> x ≤ 1, scaled_error) + end + @testset "FullStressState" begin + ϵ21 = 3 * σ_y0 / E; num_steps = 10; t_end = 0.01 + ss = FullStressState() + # Check that we get the same result for MatTest.runstresstate and MatTest.runstrain + σ_ss, s_ss = MatTest.runstresstate(ss, m, ϵ21, (2, 1), t_end, num_steps) + σ, s = MatTest.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)), + (UniaxialNormalStress(), (1,2)), + (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 + num_steps = 10; t_end = 0.01 + ϵij = num_steps * 0.1 * 3 * σ_y0 / E + relstep = 1e-6 + stressfun(p) = MatTest.runstresstate(stress_state, fromvector(p, m), ϵij, ij, t_end, num_steps)[1] + dσij_dp_num = FiniteDiff.finite_difference_jacobian(stressfun, BigFloat.(tovector(m)), Val{:central}; relstep) + σv, state, dσij_dp, diff = MatTest.runstresstate_diff(stress_state, m, ϵij, ij, t_end, num_steps) + isapprox(dσij_dp_num, dσij_dp; atol = 1000 * eps()) + + scaled_error, maxtol = MatTest.compare_derivatives(dσij_dp, dσij_dp_num, σv, tovector(m) * relstep; atol_min = 1e-6, rtol_min = 1e-4, print_tol = false) + if !all(x -> x ≤ 1, scaled_error) + display(scaled_error) + end + @test all(x -> x ≤ 1, scaled_error) + end + end + end + end + end + \ No newline at end of file diff --git a/test/test_elastic.jl b/test/test_elastic.jl index f78d05d..c99f684 100644 --- a/test/test_elastic.jl +++ b/test/test_elastic.jl @@ -8,8 +8,8 @@ m2 = LinearElastic{:cubicsymmetry}(C1111=v2[1], C1122=v2[2], C1212=v2[3]) # Test that construction via vector gives the same material - @test m1 == MMB.fromvector(v1, m1) - @test m2 == MMB.fromvector(v2, m2) + @test m1 == fromvector(v1, m1) + @test m2 == fromvector(v2, m2) # Test get stress function E = 210.e3; ν = 0.3 diff --git a/test/test_hardening.jl b/test/test_hardening.jl index 22dd667..89243f6 100644 --- a/test/test_hardening.jl +++ b/test/test_hardening.jl @@ -22,10 +22,17 @@ @test !isapprox(af_h0, ow_h1) # Should not be equal # Check constructors from vectors - @test af == MMB.fromvector([Hkin, β∞], af) - @test af == MMB.fromvector([1.0, Hkin, β∞], af, offset=1) - @test db == MMB.fromvector([Hkin, β∞, δ], db) - @test ow == MMB.fromvector([Hkin, β∞, m], ow) + @test af == fromvector([Hkin, β∞], af) + @test af == fromvector([1.0, Hkin, β∞], af, offset=1) + @test db == fromvector([Hkin, β∞, δ], db) + @test ow == fromvector([Hkin, β∞, m], ow) + + # Test conversions + for T in (Float64, MatTest.DualT{Float32}) + for hardlaw in (af, db, ow) + MatTest.test_vectorconversion(T, hardlaw) + end + end # Check show methods @test contains(show_as_string(af), "ArmstrongFrederick with") @@ -41,7 +48,7 @@ end voce = Voce(Hiso=Hiso, κ∞=κ∞) # Test conversion from vector - @test voce == MMB.fromvector(vv, voce) + @test voce == fromvector(vv, voce) # Calculate Voce response using analytical function vocelaw(param, λ) = param.κ∞*(one(λ) - exp(-param.Hiso*λ/param.κ∞)) @@ -56,7 +63,7 @@ end swiftlaw(param, λ) = param.K * (param.λ0 + λ)^param.n # Test conversion from vector - @test swift == MMB.fromvector(sv, swift) + @test swift == fromvector(sv, swift) κ0 = swiftlaw(swift, 0.0) κ = swiftlaw(swift, λ) @@ -65,6 +72,13 @@ end @test dκdλ ≈ MechMat.get_evolution(swift, κ) @test κ0 ≈ MechMat.get_initial_value(swift) + # Test conversions + for T in (Float64, MatTest.DualT{Float32}) + for hardlaw in (voce, swift) + MatTest.test_vectorconversion(T, hardlaw) + end + end + # Show @test contains(show_as_string(voce), "Voce with") @test contains(show_as_string(swift), "Swift with") diff --git a/test/test_overstress.jl b/test/test_overstress.jl index 3e1fbbb..d7bb742 100644 --- a/test/test_overstress.jl +++ b/test/test_overstress.jl @@ -1,7 +1,7 @@ @testset "OverstressFunction" begin @testset "RateIndependent" begin of = RateIndependent() - @test MMB.get_num_params(of) == 0 + @test get_vector_length(of) == 0 @test of == fromvector(rand(10), of) @test contains(show_as_string(of), "Rate independent response") @@ -10,8 +10,11 @@ @testset "NortonOverstress" begin tstar, nexp = rand(2) of = NortonOverstress(;tstar, nexp) - @test MMB.get_num_params(of) == 2 - test_conversion(of) + @test get_vector_length(of) == 2 + + MatTest.test_vectorconversion(Float64, of) + MatTest.test_vectorconversion(Float32, of) + MatTest.test_vectorconversion(MatTest.DualT{Float32}, of) Φ, σy = rand(2) # Scaling of Φ and σy doesn't change values diff --git a/test/test_plastic.jl b/test/test_plastic.jl index ea166ae..a84583f 100644 --- a/test/test_plastic.jl +++ b/test/test_plastic.jl @@ -15,10 +15,10 @@ @test nisop == length(piso) # Test for the full material param = vcat(pel, σ_y0, piso, pkin) - @test m1 == MMB.fromvector(param, m1) - pvec = zeros(get_num_params(m1)) + @test m1 == fromvector(param, m1) + pvec = zeros(get_vector_length(m1)) @test length(pvec) == length(param) - MMB.tovector!(pvec, m1) + tovector!(pvec, m1) @test pvec == param # Run an example, not test of result, just checking that it runs diff --git a/test/test_viscoplastic.jl b/test/test_viscoplastic.jl index 503e1ec..c109122 100644 --- a/test/test_viscoplastic.jl +++ b/test/test_viscoplastic.jl @@ -29,7 +29,7 @@ # Conversions function check_conversions(m) - v = rand(get_num_params(m)) + v = rand(get_vector_length(m)) mr = fromvector(v, m) vc = tovector(mr) @test vc ≈ v diff --git a/test/test_yield.jl b/test/test_yield.jl index 30a0d12..2ef9639 100644 --- a/test/test_yield.jl +++ b/test/test_yield.jl @@ -19,8 +19,9 @@ @test MechMat.yield_criterion(yc, σ_shear, k) ≈ (s - (sy+k)) # Conversions - @test MMB.get_num_params(yc) == 1 - test_conversion(yc) + @test get_vector_length(yc) == 1 + MatTest.test_vectorconversion(Float64, yc) + MatTest.test_vectorconversion(Float32, yc) # Show @test contains(show_as_string(yc), "VonMises with") @@ -44,8 +45,9 @@ @test MechMat.yield_criterion(yc, σ_shear, k) ≈ (s - (Y0+k)) # Conversions - @test MMB.get_num_params(yc) == 2 - test_conversion(yc) + @test get_vector_length(yc) == 2 + MatTest.test_vectorconversion(Float64, yc) + MatTest.test_vectorconversion(MatTest.DualT{Float32}, yc) # Show @test contains(show_as_string(yc), "DruckerPrager with") diff --git a/test/utilities4testing.jl b/test/utilities4testing.jl index e64be59..4d7cb01 100644 --- a/test/utilities4testing.jl +++ b/test/utilities4testing.jl @@ -32,58 +32,6 @@ function run_normal(m, ϵ11_vec; Δt = NaN, stress_state=UniaxialNormalStress(), return s11, σ, dσdϵ, state, ϵ_full end -function obtain_numerical_material_derivative!(deriv, m, ϵ, old, Δt) - p = tovector(m) - # σ, dσdϵ, state = material_response(m, ϵ, old, Δt) - dmdz = (σ = (ϵ = deriv.dσdϵ, s = deriv.dσdⁿs, p = deriv.dσdp), - s = (ϵ = deriv.dsdϵ, s = deriv.dsdⁿs, p = deriv.dsdp)) - _tovector(t::AbstractTensor) = tomandel(t) - _tovector(s::AbstractMaterialState) = tovector(s) - funs = Tuple((ϵ=evec -> _tovector(material_response(m, frommandel(typeof(ϵ),evec), old, Δt)[i]), # f(ϵ) - s=svec -> _tovector(material_response(m, ϵ, MMB.fromvector(svec, old), Δt)[i]), # f(s) - p=pvec -> _tovector(material_response(fromvector(pvec, m), ϵ, old, Δt)[i])) # f(p) - for i in [1,3]) - x0 = (ϵ = tomandel(ϵ), s = tovector(old), p = p) - for (i, m) in enumerate((:σ,:s)) - tmp1 = getfield(dmdz,m) - for z in fieldnames(typeof(x0)) - if length(getfield(x0, z)) > 0 - tmp=getfield(tmp1,z) - fun(x) = getfield(funs[i], z)(x) - tmp .= FiniteDiff.finite_difference_jacobian(fun, getfield(x0, z)) - end - end - end -end - -function obtain_numerical_material_derivative!(deriv, m::LinearElastic, ϵ, old, Δt) - # Set all derivatives to zero - for field in fieldnames(typeof(deriv)) - fill!(getfield(deriv, field), 0) - end - dmdz = (σ = (ϵ = deriv.dσdϵ, s = deriv.dσdⁿs, p = deriv.dσdp), - s = (ϵ = deriv.dsdϵ, s = deriv.dsdⁿs, p = deriv.dsdp)) - # σ, dσdϵ, state = material_response(m, ϵ, old, Δt) - funs = (ϵ=evec -> tomandel(material_response(m, frommandel(typeof(ϵ),evec), old, Δt)[1]), # f(ϵ) - p=pvec -> tomandel(material_response(fromvector(pvec,m), ϵ, old, Δt)[1])) # f(p) - x0 = (ϵ=tomandel(ϵ), p=tovector(m)) - - dσdz = getfield(dmdz,:σ) - for z in fieldnames(typeof(x0)) - tmp=getfield(dσdz,z) - fun(x) = getfield(funs, z)(x) - tmp .= FiniteDiff.finite_difference_jacobian(fun, getfield(x0, z)) - end -end - -function test_conversion(mbase) - v0 = [rand() for _ in 1:MMB.get_num_params(mbase)] - v1 = similar(v0) - m = fromvector(v0, mbase) - tovector!(v1, m) - @test v0 ≈ v1 -end - function show_as_string(value, mime=MIME"text/plain"()) io = IOBuffer() show(IOContext(io), mime, value)