diff --git a/Project.toml b/Project.toml index ff95639..f964898 100644 --- a/Project.toml +++ b/Project.toml @@ -1,25 +1,29 @@ name = "MaterialModelsBase" uuid = "af893363-701d-44dc-8b1e-d9a2c129bfc9" +version = "0.4" authors = ["Knut Andreas Meyer and contributors"] -version = "0.3.1" [deps] ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4" +[sources] +MaterialModelsTesting = {url = "https://github.com/KnutAM/MaterialModelsTesting.jl"} +MechanicalMaterialModels = {url = "https://github.com/KnutAM/MechanicalMaterialModels.jl"} +Newton = {url = "https://github.com/KnutAM/Newton.jl"} + [compat] julia = "1.11" [extras] -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" MaterialModelsTesting = "882b014b-b96c-4115-8629-e17fb35110d2" +MechanicalMaterialModels = "b3282f9b-607f-4337-ab95-e5488ab5652c" +Newton = "83aa5b51-0588-403c-85e4-434ec185aae7" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -[sources] -MaterialModelsTesting = {url = "https://github.com/KnutAM/MaterialModelsTesting.jl"} -#MaterialModelsTesting = {path = "/Users/knutan/.julia/dev/MaterialModelsTesting"} [targets] -test = ["Test", "FiniteDiff", "MaterialModelsTesting", "Random"] +test = ["FiniteDiff", "MaterialModelsTesting", "MechanicalMaterialModels", "Newton", "Random", "Test"] diff --git a/src/differentiation.jl b/src/differentiation.jl index e2e6902..c8b5d41 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -157,22 +157,53 @@ function reset_derivatives!(sd::StressStateDerivatives, sd_prev::StressStateDeri return sd end -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 +function differentiate_material!(ssd::StressStateDerivatives, stress_state::AbstractStressState, m::AbstractMaterial, strain::AbstractTensor, args::Vararg{Any,N}) where {N} + stress_3d, stiff_3d, state, strain_3d = stress_state_material_response(stress_state, m, strain, args...) + differentiate_material!(ssd.mderiv, m, strain_3d, args..., stiff_3d) + _update_derivatives!(ssd, stress_state, strain_3d, stress_3d, stiff_3d) + return reduce_tensordim(stress_state, stress_3d), reduce_stiffness(stress_state, stiff_3d, strain_3d, stress_3d), state, strain_3d +end + +function _update_derivatives!(ssd::StressStateDerivatives, ::NoIterationState, args...) + copy!(ssd.dσdp, ssd.mderiv.dσdp) + fill!(ssd.dϵdp, 0) + return ssd +end + +# IterationState +# drdp = 0 = ∂r∂p + ∂r∂x * dxdp => dxdp = -∂r∂x\∂r∂p +# r(x(p),P(x,p)) +# Calculated: ∂P∂p (for full strain, F) and ∂P∂x (∂P∂F) +# r = r̃(τ(P(x(p), p), x(p)), x(p)) +# drdp = ∂r̃∂τ ⋅ [∂τ∂P ⋅ [∂P∂x ⋅ dxdp + ∂P∂p] + ∂τ∂x ⋅ dxdp] + ∂r̃∂x ⋅ dxdp +# = [∂r̃∂τ ⋅ [∂τ∂P ⋅ ∂P∂x + ∂τ∂x] + ∂r̃∂x] ⋅ dxdp + ∂r̃∂τ ⋅ ∂τ∂P ⋅ ∂P∂p = 0 +# dxdp = -([∂r̃∂τ ⋅ [∂τ∂P ⋅ ∂P∂x + ∂τ∂x] + ∂r̃∂x])⁻¹ ⋅ [∂r̃∂τ ⋅ ∂τ∂P ⋅ ∂P∂p] +# r = r̂(x(p), p) = r̃(τ(P(x(p), p), x(p)), x(p)) +# ∂r̂∂x = get_drdx(...) # This is the function we already have +# ∂r̂∂p = ∂r̃∂τ ⋅ ∂τ∂P ⋅ ∂P∂p +# dxdp = -∂r̂∂x \ ∂r̂∂p +# Small strains, σ = τ = P & r = τ[sc] +# => ∂r̂∂p = ∂r̃∂τᶠ ⋅ ∂τᶠ∂Pᶠ ⋅ ∂Pᶠ∂p = [I ⋅ I ⋅ ∂Pᶠ∂p] = ∂Pᶠ∂p = ∂σᶠ∂p +# => dϵdp = - ∂r̂∂x \ ∂σᶠ∂p +# Finite strains +# ∂r̃∂τ depends on specific gauge constraints +# ∂τ∂P = I ⊗̄ F (τ = P F'; => P_ij F_jk' / dP_lm = δ_il δ_jm F_jk' = δ_il F_km) +# Not yet implemented + +# SmallStrains +function _update_derivatives!(ssd::StressStateDerivatives, stress_state::IterationState, ϵ::SymmetricTensor{2,3}, σ::SymmetricTensor{2,3}, dσdϵ::SymmetricTensor{4,3}) + sc = stress_controlled_indices(stress_state, ϵ) + ec = strain_controlled_indices(stress_state, ϵ) + minus_∂r∂x_inv = -inv(get_drdx(stress_state, dσdϵ, ϵ, σ)) # -[∂σᶠ∂ϵᶠ]⁻¹ + ssd.dϵdp[sc, :] .= minus_∂r∂x_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, :] + return ssd +end + +# FiniteStrains +function _update_derivatives!(ssd::StressStateDerivatives, stress_state::IterationState, F::Tensor{2,3}, P::Tensor{2,3}, dPdF::Tensor{4,3}) + error("StressStateDerivatives calculation for finite strains not implemented, see notes in source file for a start to implement") end """ diff --git a/src/stressiterations.jl b/src/stressiterations.jl index 969d2d4..5f889cb 100644 --- a/src/stressiterations.jl +++ b/src/stressiterations.jl @@ -25,8 +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 - σ, 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 + stress_3d, stiff_3d, state, strain_3d = stress_state_material_response(stress_state, m, args...) + return reduce_tensordim(stress_state, stress_3d), reduce_stiffness(stress_state, stiff_3d, strain_3d, stress_3d), state, strain_3d end update_stress_state!(::AbstractStressState, σ) = nothing @@ -98,17 +98,17 @@ struct UniaxialStrain <: AbstractStressState end # Cases with stress iterations """ - IterationSettings(;tolerance=1e-8, max_iter=10) + IterationSettings(;tolerance = 1e-8, maxiter = 10) Settings for stress iterations. Constructors for iterative stress states forwards given keyword arguments to this constructor and saves the result. """ @kwdef struct IterationSettings{T} - tolerance::T=1.e-8 - max_iter::Int=10 + tolerance::T = 1.e-8 + maxiter::Int = 10 end get_tolerance(is::IterationSettings) = is.tolerance -get_max_iter(is::IterationSettings) = is.max_iter +get_maxiter(is::IterationSettings) = is.maxiter """ UniaxialStress(; kwargs...) @@ -117,7 +117,12 @@ Uniaxial stress such that ``\\sigma_{ij}=0 \\forall (i,j)\\neq (1,1)`` The strain input can be 1d (`SecondOrderTensor{1}`). A 3d input is also accepted and used as an initial -guess for the unknown strain components. +guess for the unknown strain components. + +For finite strains, a uniaxial Kirchhoff stress, τ = P ⋅ F', is found, such that +``\\tau_{ij} = 0 \\forall (i,j) \\neq (1, 1)`` +Additionally, the deformation gradient is forced to be symmetric to the potential +arbitrary rotations, i.e. ``F_{ij} = F_{ji}``. The optional keyword arguments are forwarded to [`IterationSettings`](@ref). """ @@ -126,16 +131,21 @@ struct UniaxialStress{T} <: AbstractStressState end UniaxialStress(; kwargs...) = UniaxialStress(IterationSettings(; kwargs...)) get_tolerance(ss::UniaxialStress) = get_tolerance(ss.settings) -get_max_iter(ss::UniaxialStress) = get_max_iter(ss.settings) +get_maxiter(ss::UniaxialStress) = get_maxiter(ss.settings) """ PlaneStress(; kwargs...) -Plane stress such that -``\\sigma_{33}=\\sigma_{23}=\\sigma_{13}=\\sigma_{32}=\\sigma_{31}=0`` -The strain input should be at least 2d (components 11, 12, 21, and 22). -A 3d input is also accepted and used as an initial guess for the unknown -out-of-plane strain components. +For small strain, find the plane stress state such that +``\\sigma_{33}=\\sigma_{23}=\\sigma_{13}=0`` + +For finite strain, find the plane Kirchhoff stress, τ = P ⋅ F', such that +``\\tau{33}=\\tau{23}=\\tau{13}=0``. +The out-of-plane shear deformation gradient components are forced to be symmetric, +i.e. ``F_{23} = F_{32}`` and ``F_{13} = F_{31}``. + +The strain input should be at least 2d, but a 3d input is also accepted and +used as an initial guess for the unknown out-of-plane strain components. The optional keyword arguments are forwarded to [`IterationSettings`](@ref). """ @@ -145,15 +155,20 @@ end PlaneStress(; kwargs...) = PlaneStress(IterationSettings(; kwargs...)) get_tolerance(ss::PlaneStress) = get_tolerance(ss.settings) -get_max_iter(ss::PlaneStress) = get_max_iter(ss.settings) +get_maxiter(ss::PlaneStress) = get_maxiter(ss.settings) """ UniaxialNormalStress(; kwargs...) This is a variation of the uniaxial stress state, such that only -``\\sigma_{22}=\\sigma_{33}=0`` -The strain input must be 3d, and the components -``\\epsilon_{22}`` and ``\\epsilon_{33}`` are used as initial guesses. +``\\sigma_{22}=\\sigma_{33}=0`` for small strains. The strain input +must be 3d, and the components ``\\epsilon_{22}`` and ``\\epsilon_{33}`` +are used as initial guesses. + +For finite strains, the condition is on the Kirchhoff stress, τ = P ⋅ F', +``\\tau_{22}=\\tau_{33}=0``. The deformation gradient should be 3d, and +``F_{22}`` and ``F_{33}`` are used as initial guesses. + This case is useful when simulating strain-controlled axial-shear experiments. Note that the stress and stiffness outputs are the 3d tensors, and that the stiffness is **not** modified to account for the stress constraints. @@ -166,16 +181,23 @@ end UniaxialNormalStress(; kwargs...) = UniaxialNormalStress(IterationSettings(; kwargs...)) get_tolerance(ss::UniaxialNormalStress) = get_tolerance(ss.settings) -get_max_iter(ss::UniaxialNormalStress) = get_max_iter(ss.settings) +get_maxiter(ss::UniaxialNormalStress) = get_maxiter(ss.settings) """ GeneralStressState(σ_ctrl::AbstractTensor{2,3,Bool}, σ::AbstractTensor{2,3}; kwargs...) -Construct a general stress state controlled by `σ_ctrl` whose component is `true` if that +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, σ gives the value to which it is controlled. The current stress, for stress-controlled components -can be updated by calling `update_stress_state!(s::GeneralStressState, σ)`. Components in -σ that are not stress-controlled are ignored. +can be updated by calling `update_stress_state!(s::GeneralStressState, σ)`. Components in +σ that are not stress-controlled are ignored. + +For finite strains, this works directly on the first Piola-Kirchhoff stress, P, and +no gauge conditions (e.g. symmetry of deformation gradient) is enforced. Hence, it is the user's +responsibility to choose conditions that result in a valid equation system. Note that this might +change in the future to be consistent with other iteration states, but could also be generalized +further to allow more generic constraints. For finite strains this stress state should thus be +used for testing purposes only. Note that the stress and stiffness outputs are the 3d tensors, and that the stiffness is **not** modified to account for the stress constraints. @@ -224,7 +246,7 @@ end update_stress_state!(s::GeneralStressState, σ) = (s.σ = σ) get_tolerance(ss::GeneralStressState) = get_tolerance(ss.settings) -get_max_iter(ss::GeneralStressState) = get_max_iter(ss.settings) +get_maxiter(ss::GeneralStressState) = get_maxiter(ss.settings) const NoIterationState = Union{FullStressState,PlaneStrain,UniaxialStrain} const IterationState = Union{UniaxialStress, PlaneStress, UniaxialNormalStress, GeneralStressState} @@ -247,167 +269,266 @@ reduce_tensordim(::State3D, a::AbstractTensor) = a reduce_tensordim(::State2D, a::AbstractTensor) = reduce_tensordim(Val{2}(), a) reduce_tensordim(::State1D, a::AbstractTensor) = reduce_tensordim(Val{1}(), a) #reduce_tensordim(::Val{dim}, v::Vec{3}) = Vec{dim}(i->v[i]) -reduce_tensordim(::Val{dim}, a::Tensor{2}) where dim = Tensor{2,dim}((i,j)->a[i,j]) -reduce_tensordim(::Val{dim}, A::Tensor{4}) where dim = Tensor{4,dim}((i,j,k,l)->A[i,j,k,l]) -reduce_tensordim(::Val{dim}, a::SymmetricTensor{2}) where dim = SymmetricTensor{2,dim}((i,j)->a[i,j]) -reduce_tensordim(::Val{dim}, A::SymmetricTensor{4}) where dim = SymmetricTensor{4,dim}((i,j,k,l)->A[i,j,k,l]) +reduce_tensordim(::Val{dim}, P::Tensor{2}) where dim = Tensor{2,dim}((i, j) -> P[i, j]) +reduce_tensordim(::Val{dim}, dPdF::Tensor{4}) where dim = Tensor{4,dim}((i, j, k, l) -> dPdF[i, j, k, l]) +reduce_tensordim(::Val{dim}, σ::SymmetricTensor{2}) where dim = SymmetricTensor{2,dim}((i, j) -> σ[i, j]) +reduce_tensordim(::Val{dim}, dσdϵ::SymmetricTensor{4}) where dim = SymmetricTensor{4,dim}((i, j, k, l) -> dσdϵ[i, j, k, l]) function stress_state_material_response(stress_state::NoIterationState, - m::AbstractMaterial, ϵ::AbstractTensor, args::Vararg{Any,N}) where N + m::AbstractMaterial, strain::AbstractTensor, args::Vararg{Any,N}) where N - ϵ_full = expand_tensordim(stress_state, ϵ) - σ, dσdϵ, new_state = material_response(m, ϵ_full, args...) - return σ, dσdϵ, new_state, ϵ_full + strain_3d = expand_tensordim(stress_state, strain) + stress_3d, stiff_3d, new_state = material_response(m, strain_3d, args...) + return stress_3d, stiff_3d, new_state, strain_3d end function stress_state_material_response(stress_state::IterationState, - m::AbstractMaterial, ϵ::AbstractTensor, args::Vararg{Any,N}) where N + m::AbstractMaterial, strain::AbstractTensor, args::Vararg{Any,N}) where N # Newton options tol = get_tolerance(stress_state) - maxiter = get_max_iter(stress_state) + maxiter = get_maxiter(stress_state) - ϵ_full = expand_tensordim(stress_state, ϵ) + strain_3d = expand_tensordim(stress_state, strain) for _ in 1:maxiter - σ_full, dσdϵ_full, new_state = material_response(m, ϵ_full, args...) - σ_mandel = get_unknowns(stress_state, σ_full) + stress_3d, stiff_3d, new_state = material_response(m, strain_3d, args...) + r = get_residual(stress_state, stress_3d, strain_3d) - if norm(σ_mandel) < tol - return σ_full, dσdϵ_full, new_state, ϵ_full + if norm(r) < tol + return stress_3d, stiff_3d, new_state, strain_3d end - dσdϵ_mandel = get_unknowns(stress_state, dσdϵ_full) - ϵ_full -= get_full_tensor(stress_state, ϵ, dσdϵ_mandel\σ_mandel) + drdx = get_drdx(stress_state, stiff_3d, strain_3d, stress_3d) + strain_3d -= get_full_tensor(stress_state, strain, drdx\r) end throw(NoStressConvergence("Stress iterations with the NewtonSolver did not converge")) end -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) +reduce_stiffness(ss::State3D, stiff_3d::AbstractTensor{4,3}, strain_3d, stress_3d) = reduce_stiffness(Val(3), ss, stiff_3d, strain_3d, stress_3d) +reduce_stiffness(ss::State2D, stiff_3d::AbstractTensor{4,3}, strain_3d, stress_3d) = reduce_stiffness(Val(2), ss, stiff_3d, strain_3d, stress_3d) +reduce_stiffness(ss::State1D, stiff_3d::AbstractTensor{4,3}, strain_3d, stress_3d) = reduce_stiffness(Val(1), ss, stiff_3d, strain_3d, stress_3d) -reduce_stiffness(::Val{3}, ::State3D, dσdϵ_full::AbstractTensor{4,3}) = dσdϵ_full +reduce_stiffness(::Val{3}, ::State3D, stiff_3d::AbstractTensor{4,3}, args...) = stiff_3d -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) +function reduce_stiffness(::Union{Val{1}, Val{2}}, stress_state::IterationState, stiff_3d::AbstractTensor{4,3}, strain_3d, stress_3d) + # Using σ & ϵ, but for finite strain these are P & F + # •ᶜ: Strain components known (not part of x), stress components calculated + ∂r∂x, ∂r∂ϵᶜ, ∂σᶜ∂x, ∂σᶜ∂ϵᶜ = extract_substiffnesses(stress_state, stiff_3d, strain_3d, stress_3d) + dσᶜdϵᶜ = ∂σᶜ∂ϵᶜ - ∂σᶜ∂x * (∂r∂x \ ∂r∂ϵᶜ) + return _totensor(dσᶜdϵᶜ, stress_state, stiff_3d) # Convert SMatrix to Tensor{4} or SymmetricTensor{4} 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) +function reduce_stiffness(::Union{Val{1}, Val{2}}, stress_state::NoIterationState, stiff_3d::AbstractTensor{4,3}, args...) + return reduce_tensordim(stress_state, stiff_3d) 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ϵᶜ) -convert_stiffness(dσᶜdϵᶜ::SMatrix{4,4}, ::State2D, ::Tensor) = frommandel(Tensor{4,2}, dσᶜdϵᶜ) +_totensor(dσᶜdϵᶜ::SMatrix{1,1}, ::State1D, ::SymmetricTensor) = frommandel(SymmetricTensor{4,1}, dσᶜdϵᶜ) +_totensor(dσᶜdϵᶜ::SMatrix{1,1}, ::State1D, ::Tensor) = frommandel(Tensor{4,1}, dσᶜdϵᶜ) +_totensor(dσᶜdϵᶜ::SMatrix{3,3}, ::State2D, ::SymmetricTensor) = frommandel(SymmetricTensor{4,2}, dσᶜdϵᶜ) +_totensor(dσᶜdϵᶜ::SMatrix{4,4}, ::State2D, ::Tensor) = frommandel(Tensor{4,2}, dσᶜdϵᶜ) -# Conversions to mandel SArray for solving equation system +# =============================== # +# Small strains, SymmetricTensor =# +# =============================== # # Internal numbering in Tensors.jl -# Tensor: 11,21,31,12,22,32,13,23,33 # SymmetricTensor: 11,21,31,22,23,33 # UniaxialStress: only σ11 != 0 -# -SymmetricTensor # i: 1, 2, 3, 4, 5 # v contains 22, 33, 32, 31, 21 function get_full_tensor(::UniaxialStress, ::SymmetricTensor, v::SVector{5,T}) where T s = T(1/√2) # 11, 21, 31, 22, 32, 33 SymmetricTensor{2,3}((0, v[5]*s, v[4]*s, v[1], v[3]*s, v[2])) end -function get_unknowns(::UniaxialStress, a::SymmetricTensor{2,3}) - SVector{5}(a[2,2], a[3,3], a[3,2]*√2, a[3,1]*√2, a[2,1]*√2) -end -function get_unknowns(::UniaxialStress, A::SymmetricTensor{4, 3}) - @SMatrix [ A[2,2,2,2] A[2,2,3,3] √2*A[2,2,3,2] √2*A[2,2,3,1] √2*A[2,2,2,1]; - A[3,3,2,2] A[3,3,3,3] √2*A[3,3,3,2] √2*A[3,3,3,1] √2*A[3,3,2,1]; - √2*A[3,2,2,2] √2*A[3,2,3,3] 2*A[3,2,3,2] 2*A[3,2,3,1] 2*A[3,2,2,1]; - √2*A[3,1,2,2] √2*A[3,1,3,3] 2*A[3,1,3,2] 2*A[3,1,3,1] 2*A[3,1,2,1]; - √2*A[2,1,2,2] √2*A[2,1,3,3] 2*A[2,1,3,2] 2*A[2,1,3,1] 2*A[2,1,2,1]] +function get_residual(::UniaxialStress, σ::SymmetricTensor{2,3}, args...) + SVector{5}(σ[2,2], σ[3,3], σ[3,2]*√2, σ[3,1]*√2, σ[2,1]*√2) +end +function get_drdx(::UniaxialStress, dσdϵ::SymmetricTensor{4, 3}, args...) + @SMatrix [ dσdϵ[2,2,2,2] dσdϵ[2,2,3,3] √2*dσdϵ[2,2,3,2] √2*dσdϵ[2,2,3,1] √2*dσdϵ[2,2,2,1]; + dσdϵ[3,3,2,2] dσdϵ[3,3,3,3] √2*dσdϵ[3,3,3,2] √2*dσdϵ[3,3,3,1] √2*dσdϵ[3,3,2,1]; + √2*dσdϵ[3,2,2,2] √2*dσdϵ[3,2,3,3] 2*dσdϵ[3,2,3,2] 2*dσdϵ[3,2,3,1] 2*dσdϵ[3,2,2,1]; + √2*dσdϵ[3,1,2,2] √2*dσdϵ[3,1,3,3] 2*dσdϵ[3,1,3,2] 2*dσdϵ[3,1,3,1] 2*dσdϵ[3,1,2,1]; + √2*dσdϵ[2,1,2,2] √2*dσdϵ[2,1,3,3] 2*dσdϵ[2,1,3,2] 2*dσdϵ[2,1,3,1] 2*dσdϵ[2,1,2,1]] +end +function extract_substiffnesses(stress_state::UniaxialStress, D::SymmetricTensor{4,3}, ϵ, σ) + # •ᶜ: Strain components known (not part of x), stress components calculated + ∂r∂x = get_drdx(stress_state, D, ϵ, σ) + ∂r∂ϵᶜ = SMatrix{5,1}(D[2,2,1,1], D[3,3,1,1], √2*D[2,3,1,1], √2*D[1,3,1,1], √2*D[1,2,1,1]) + ∂σᶜ∂x = SMatrix{1,5}(D[1,1,2,2], D[1,1,3,3], √2*D[1,1,2,3], √2*D[1,1,1,3], √2*D[1,1,1,2]) + ∂σᶜ∂ϵᶜ = SMatrix{1,1}(D[1,1,1,1]) + return ∂r∂x, ∂r∂ϵᶜ, ∂σᶜ∂x, ∂σᶜ∂ϵᶜ end -# -Tensor -# i: 1, 2, 3, 4, 5, 6, 7, 8 -# v contains 22, 33, 23, 13, 12, 32, 31, 21 -function get_full_tensor(::UniaxialStress, ϵ::Tensor, v::SVector{8}) - # 11, 21, 31, 12, 22, 32, 13, 23, 33 - return Tensor{2,3}((0, v[8], v[7], v[5], v[1], v[6], v[4], v[3], v[2])) +# UniaxialNormalStress: σ22 = σ33 = 0 +# i: 1, 2 +# v contains 22, 33 +function get_full_tensor(::UniaxialNormalStress, ::SymmetricTensor, v::SVector{2}) + # 11,21,31, 22,32, 33 + SymmetricTensor{2,3}((0, 0, 0, v[1], 0, v[2])) end - -function get_unknowns(::UniaxialStress, a::Tensor{2,3}) - SVector{8}(a[2,2], a[3,3], a[2,3], a[1,3], a[1,2], a[3,2], a[3,1], a[2,1]) +function get_residual(::UniaxialNormalStress, σ::SymmetricTensor{2,3}, args...) + SVector{2}(σ[2,2], σ[3,3]) end - -function get_unknowns(::UniaxialStress, A::Tensor{4, 3}) - @SMatrix [ A[2,2,2,2] A[2,2,3,3] A[2,2,2,3] A[2,2,1,3] A[2,2,1,2] A[2,2,3,2] A[2,2,3,1] A[2,2,2,1]; - A[3,3,2,2] A[3,3,3,3] A[3,3,2,3] A[3,3,1,3] A[3,3,1,2] A[3,3,3,2] A[3,3,3,1] A[3,3,2,1]; - A[2,3,2,2] A[2,3,3,3] A[2,3,2,3] A[2,3,1,3] A[2,3,1,2] A[2,3,3,2] A[2,3,3,1] A[2,3,2,1]; - A[1,3,2,2] A[1,3,3,3] A[1,3,2,3] A[1,3,1,3] A[1,3,1,2] A[1,3,3,2] A[1,3,3,1] A[1,3,2,1]; - A[1,2,2,2] A[1,2,3,3] A[1,2,2,3] A[1,2,1,3] A[1,2,1,2] A[1,2,3,2] A[1,2,3,1] A[1,2,2,1]; - A[3,2,2,2] A[3,2,3,3] A[3,2,2,3] A[3,2,1,3] A[3,2,1,2] A[3,2,3,2] A[3,2,3,1] A[3,2,2,1]; - A[3,1,2,2] A[3,1,3,3] A[3,1,2,3] A[3,1,1,3] A[3,1,1,2] A[3,1,3,2] A[3,1,3,1] A[3,1,2,1]; - A[2,1,2,2] A[2,1,3,3] A[2,1,2,3] A[2,1,1,3] A[2,1,1,2] A[2,1,3,2] A[2,1,3,1] A[2,1,2,1]] +function get_drdx(::UniaxialNormalStress, dσdϵ::SymmetricTensor{4, 3}, args...) + @SMatrix [dσdϵ[2,2,2,2] dσdϵ[2,2,3,3]; + dσdϵ[3,3,2,2] dσdϵ[3,3,3,3]] end -# PlaneStress: σ33=σ23=σ13=0=σ31=σ32 -# -SymmetricTensor + +# PlaneStress: σ33 = σ23 = σ13 = 0 # i: 1, 2, 3 # v contains 33, 23, 13 function get_full_tensor(::PlaneStress, ::SymmetricTensor, v::SVector{3,T}) where T s = T(1/√2) # 11,21, 31, 22, 32, 33 SymmetricTensor{2,3}((0, 0, v[3]*s, 0, v[2]*s, v[1])) end -function get_unknowns(::PlaneStress, a::SymmetricTensor{2,3}) - SVector{3}(a[3,3], a[2,3]*√2, a[3,1]*√2) +function get_residual(::PlaneStress, σ::SymmetricTensor{2,3}, args...) + SVector{3}(σ[3,3], σ[2,3]*√2, σ[3,1]*√2) +end +function get_drdx(::PlaneStress, dσdϵ::SymmetricTensor{4, 3}, ϵ::SymmetricTensor{2,3}, σ::SymmetricTensor{2,3}) + @SMatrix [ dσdϵ[3,3,3,3] √2*dσdϵ[3,3,3,2] √2*dσdϵ[3,3,3,1]; + √2*dσdϵ[3,2,3,3] 2*dσdϵ[3,2,3,2] 2*dσdϵ[3,2,3,1]; + √2*dσdϵ[3,1,3,3] 2*dσdϵ[3,1,3,2] 2*dσdϵ[3,1,3,1]] +end +function extract_substiffnesses(stress_state::PlaneStress, D::SymmetricTensor{4,3}, ϵ, σ) + # f=strain free, stress constrained + # c=strain constrained, stress free + ∂σᶠ∂ϵᶠ = get_drdx(stress_state, D, ϵ, σ) + ∂σᶠ∂ϵᶜ = @SMatrix [ D[3,3,1,1] D[3,3,2,2] √2*D[3,3,2,1]; + √2*D[3,2,1,1] √2*D[3,2,2,2] 2*D[3,2,2,1]; + √2*D[3,1,1,1] √2*D[3,1,2,2] 2*D[3,1,2,1]] + + ∂σᶜ∂ϵᶠ = @SMatrix [ D[1,1,3,3] √2*D[1,1,3,2] √2*D[1,1,3,1]; + D[2,2,3,3] √2*D[2,2,3,2] √2*D[2,2,3,1]; + √2*D[2,1,3,3] 2*D[2,1,3,2] 2*D[2,1,3,1]] + + ∂σᶜ∂ϵᶜ = @SMatrix [ D[1,1,1,1] D[1,1,2,2] √2*D[1,1,2,1]; + D[2,2,1,1] D[2,2,2,2] √2*D[2,2,2,1]; + √2*D[2,1,1,1] √2*D[2,1,2,2] 2*D[2,1,2,1]] + + return ∂σᶠ∂ϵᶠ, ∂σᶠ∂ϵᶜ, ∂σᶜ∂ϵᶠ, ∂σᶜ∂ϵᶜ +end + +# =============================== # +# Finite strains, Tensor =# +# =============================== # +# Internal numbering in Tensors.jl +# Tensor: 11,21,31,12,22,32,13,23,33 + +function get_drdx(ss::AbstractStressState, dPdF::Tensor{4, 3, T}, F::Tensor{2, 3}, P::Tensor{2, 3}) where {T} + drdx, _ = get_drdx_dτdF(ss, dPdF, F, P) + return drdx end -function get_unknowns(::PlaneStress, A::SymmetricTensor{4, 3}) - @SMatrix [ A[3,3,3,3] √2*A[3,3,3,2] √2*A[3,3,3,1]; - √2*A[3,2,3,3] 2*A[3,2,3,2] 2*A[3,2,3,1]; - √2*A[3,1,3,3] 2*A[3,1,3,2] 2*A[3,1,3,1]] + +# Stiffness and stress conversions +function dPdF_to_dτdF(dPdF::Tensor{4,3}, F::Tensor, P::Tensor) + # TODO: If important, this can be optimized with Tensors.get_expression + # dτ_ij/dF_kl = d(P_im F_mj')/dF_kl = F_mj' * dP_im/dF_kl + P_im dF_mj'/dF_kl + # = δ_in F_jm dP_nm/dF_kl + P_il δ_jk + # : (I ⊗̄ F) : dPdF + P ⊗̲ I + I2 = one(F) + return otimesu(I2, F) ⊡ dPdF + otimesl(P, I2) # dτdF end +convert_P_to_τ(P::Tensor{2,3}, F::Tensor{2,3}) = P ⋅ F' -# -Tensor -# i: 1, 2, 3, 4, 5, 6 -# v contains 33, 23, 13, 32, 31 +# UniaxialStress +# i: 1, 2, 3, 4, 5, 6, 7, 8 +# x = [F22, F33, F23, F13, F12, F32, F31, F21] +function get_full_tensor(::UniaxialStress, ϵ::Tensor, v::SVector{8}) + # 11, 21, 31, 12, 22, 32, 13, 23, 33 + return Tensor{2,3}((0, v[8], v[7], v[5], v[1], v[6], v[4], v[3], v[2])) +end + +function get_residual(::UniaxialStress, P::Tensor{2,3}, F::Tensor{2,3}) + τ = convert_P_to_τ(P, F) + SVector{8}(τ[2,2], τ[3,3], τ[2,3], τ[1,3], τ[1,2], F[2,3] - F[3,2], F[1,3] - F[3,1], F[1,2] - F[2,1]) +end + +function get_drdx_dτdF(::UniaxialStress, dPdF::Tensor{4, 3, T}, F::Tensor{2, 3}, P::Tensor{2, 3}) where {T} + dτdF = dPdF_to_dτdF(dPdF, F, P) + drdx = @SMatrix [ + dτdF[2,2,2,2] dτdF[2,2,3,3] dτdF[2,2,2,3] dτdF[2,2,1,3] dτdF[2,2,1,2] dτdF[2,2,3,2] dτdF[2,2,3,1] dτdF[2,2,2,1]; + dτdF[3,3,2,2] dτdF[3,3,3,3] dτdF[3,3,2,3] dτdF[3,3,1,3] dτdF[3,3,1,2] dτdF[3,3,3,2] dτdF[3,3,3,1] dτdF[3,3,2,1]; + dτdF[2,3,2,2] dτdF[2,3,3,3] dτdF[2,3,2,3] dτdF[2,3,1,3] dτdF[2,3,1,2] dτdF[2,3,3,2] dτdF[2,3,3,1] dτdF[2,3,2,1]; + dτdF[1,3,2,2] dτdF[1,3,3,3] dτdF[1,3,2,3] dτdF[1,3,1,3] dτdF[1,3,1,2] dτdF[1,3,3,2] dτdF[1,3,3,1] dτdF[1,3,2,1]; + dτdF[1,2,2,2] dτdF[1,2,3,3] dτdF[1,2,2,3] dτdF[1,2,1,3] dτdF[1,2,1,2] dτdF[1,2,3,2] dτdF[1,2,3,1] dτdF[1,2,2,1]; + zero(T) zero(T) one(T) zero(T) zero(T) -one(T) zero(T) zero(T); # R6 = F[2,3] - F[3,2] + zero(T) zero(T) zero(T) one(T) zero(T) zero(T) -one(T) zero(T); # R7 = F[1,3] - F[3,1] + zero(T) zero(T) zero(T) zero(T) one(T) zero(T) zero(T) -one(T)] # R8 = F[1,2] - F[2,1] + return drdx, dτdF +end + +function extract_substiffnesses(stress_state::UniaxialStress, dPdF::Tensor{4, 3, T}, F::Tensor{2, 3}, P::Tensor{2, 3}) where {T} + # •ᶜ: Strain components known (not part of x), stress components calculated + ∂r∂x, dτdF = get_drdx_dτdF(stress_state, dPdF, F, P) + ∂r∂Fᶜ = SMatrix{8,1}(dτdF[2,2,1,1], dτdF[3,3,1,1], dτdF[2,3,1,1], dτdF[1,3,1,1], dτdF[1,2,1,1], zero(T), zero(T), zero(T)) + ∂Pᶜ∂x = SMatrix{1,8}(dPdF[1,1,2,2], dPdF[1,1,3,3], dPdF[1,1,2,3], dPdF[1,1,1,3], dPdF[1,1,1,2], dPdF[1,1,3,2], dPdF[1,1,3,1], dPdF[1,1,2,1]) + ∂Pᶜ∂Fᶜ = SMatrix{1,1}(dPdF[1,1,1,1]) + return ∂r∂x, ∂r∂Fᶜ, ∂Pᶜ∂x, ∂Pᶜ∂Fᶜ +end + + +# PlaneStress +# i: 1, 2, 3, 4, 5 +# x = [F33, F23, F13, F32, F31] function get_full_tensor(::PlaneStress, ϵ::Tensor, v::SVector{5}) # 11,21, 31,12,22, 32, 13, 23, 33 return Tensor{2,3}((0, 0, v[5], 0, 0, v[4], v[3], v[2], v[1])) end +function get_residual(::PlaneStress, P::Tensor{2,3}, F::Tensor) + τ = convert_P_to_τ(P, F) + SVector{5}(τ[3,3], τ[2,3], τ[1,3], F[2,3] - F[3,2], F[1,3] - F[3,1]) +end +function get_drdx_dτdF(::PlaneStress, dPdF::Tensor{4, 3, T}, F::Tensor{2,3}, P::Tensor{2,3}) where {T} + dτdF = dPdF_to_dτdF(dPdF, F, P) + drdx = @SMatrix [ + dτdF[3,3,3,3] dτdF[3,3,2,3] dτdF[3,3,1,3] dτdF[3,3,3,2] dτdF[3,3,3,1]; + dτdF[2,3,3,3] dτdF[2,3,2,3] dτdF[2,3,1,3] dτdF[2,3,3,2] dτdF[2,3,3,1]; + dτdF[1,3,3,3] dτdF[1,3,2,3] dτdF[1,3,1,3] dτdF[1,3,3,2] dτdF[1,3,3,1]; + zero(T) one(T) zero(T) -one(T) zero(T); # R4 = F23 - F32 + zero(T) zero(T) one(T) zero(T) -one(T)] # R5 = F13 - F31 + return drdx, dτdF +end +function extract_substiffnesses(stress_state::PlaneStress, dPdF::Tensor{4, 3, T}, F::Tensor{2, 3}, P::Tensor{2, 3}) where {T} + # •ᶜ: Strain components known (not part of x), stress components calculated + ∂r∂x, dτdF = get_drdx_dτdF(stress_state, dPdF, F, P) + ∂r∂Fᶜ = @SMatrix [dτdF[3,3,1,1] dτdF[3,3,2,2] dτdF[3,3,1,2] dτdF[3,3,2,1]; + dτdF[2,3,1,1] dτdF[2,3,2,2] dτdF[2,3,1,2] dτdF[2,3,2,1]; + dτdF[1,3,1,1] dτdF[1,3,2,2] dτdF[1,3,1,2] dτdF[1,3,2,1]; + zero(T) zero(T) zero(T) zero(T); + zero(T) zero(T) zero(T) zero(T)] + + ∂Pᶜ∂x = @SMatrix [dPdF[1,1,3,3] dPdF[1,1,2,3] dPdF[1,1,1,3] dPdF[1,1,3,2] dPdF[1,1,3,1]; + dPdF[2,2,3,3] dPdF[2,2,2,3] dPdF[2,2,1,3] dPdF[2,2,3,2] dPdF[2,2,3,1]; + dPdF[1,2,3,3] dPdF[1,2,2,3] dPdF[1,2,1,3] dPdF[1,2,3,2] dPdF[1,2,3,1]; + dPdF[2,1,3,3] dPdF[2,1,2,3] dPdF[2,1,1,3] dPdF[2,1,3,2] dPdF[2,1,3,1]] -function get_unknowns(::PlaneStress, a::Tensor{2,3}) - SVector{5}(a[3,3], a[2,3], a[1,3], a[3,2], a[3,1]) -end - -function get_unknowns(::PlaneStress, A::Tensor{4, 3}) - @SMatrix [ A[3,3,3,3] A[3,3,2,3] A[3,3,1,3] A[3,3,3,2] A[3,3,3,1]; - A[2,3,3,3] A[2,3,2,3] A[2,3,1,3] A[2,3,3,2] A[2,3,3,1]; - A[1,3,3,3] A[1,3,2,3] A[1,3,1,3] A[1,3,3,2] A[1,3,3,1]; - A[3,2,3,3] A[3,2,2,3] A[3,2,1,3] A[3,2,3,2] A[3,2,3,1]; - A[3,1,3,3] A[3,1,2,3] A[3,1,1,3] A[3,1,3,2] A[3,1,3,1]] + ∂Pᶜ∂Fᶜ = @SMatrix [dPdF[1,1,1,1] dPdF[1,1,2,2] dPdF[1,1,1,2] dPdF[1,1,2,1]; + dPdF[2,2,1,1] dPdF[2,2,2,2] dPdF[2,2,1,2] dPdF[2,2,2,1]; + dPdF[1,2,1,1] dPdF[1,2,2,2] dPdF[1,2,1,2] dPdF[1,2,2,1]; + dPdF[2,1,1,1] dPdF[2,1,2,2] dPdF[2,1,1,2] dPdF[2,1,2,1]] + return ∂r∂x, ∂r∂Fᶜ, ∂Pᶜ∂x, ∂Pᶜ∂Fᶜ end -# UniaxialNormalStress: σ22=σ33=0 -# -SymmetricTensor and Tensor +# UniaxialNormalStress: τ22=τ33=0 +# -Tensor # i: 1, 2 # v contains 22, 33 -function get_full_tensor(::UniaxialNormalStress, ::SymmetricTensor, v::SVector{2}) - # 11,21,31, 22,32, 33 - SymmetricTensor{2,3}((0, 0, 0, v[1], 0, v[2])) -end -function get_full_tensor(::UniaxialNormalStress, ϵ::Tensor, v::SVector{2}) - # 11,21,31,12, 22,32,13,23, 33 -return Tensor{2,3}((0, 0, 0, 0, v[1], 0, 0, 0, v[2])) +function get_full_tensor(::UniaxialNormalStress, ::Tensor, v::SVector{2}) + # 11,21,31,12, 22,32,13,23, 33 + return Tensor{2,3}((0, 0, 0, 0, v[1], 0, 0, 0, v[2])) end -function get_unknowns(::UniaxialNormalStress, a::AbstractTensor{2,3}) - SVector{2}(a[2,2], a[3,3]) +function get_residual(::UniaxialNormalStress, P::Tensor{2,3}, F::Tensor{2,3}) + τ = convert_P_to_τ(P, F) + return SVector{2}(τ[2,2], τ[3,3]) end -function get_unknowns(::UniaxialNormalStress, A::AbstractTensor{4, 3}) - @SMatrix [A[2,2,2,2] A[2,2,3,3]; - A[3,3,2,2] A[3,3,3,3]] +function get_drdx(::UniaxialNormalStress, dPdF::Tensor{4, 3, T}, F::Tensor{2, 3}, P::Tensor{2, 3}) where {T} + dτdF = dPdF_to_dτdF(dPdF, F, P) + @SMatrix [dτdF[2,2,2,2] dτdF[2,2,3,3]; + dτdF[3,3,2,2] dτdF[3,3,3,3]] end # GeneralStressState @@ -425,79 +546,17 @@ function get_full_tensor(state::GeneralStressState{Nσ}, ::TT, v::SVector{Nσ,T} return TB(f) end -function get_unknowns(state::GeneralStressState{Nσ}, a::AbstractTensor{2,3,T}) where {Nσ, T} - shear_factor = a isa SymmetricTensor ? sqrt(2 * one(T)) : one(T) +function get_residual(state::GeneralStressState{Nσ}, stress_3d::AbstractTensor{2,3,T}, args...) where {Nσ, T} + shear_factor = stress_3d isa SymmetricTensor ? sqrt(2 * one(T)) : 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])) + f(c) = ((i,j) = c; s(i,j)*(stress_3d[i,j] - state.σ[i,j])) 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 = a isa SymmetricTensor ? sqrt(2 * one(T)) : one(T) +function get_drdx(state::GeneralStressState{Nσ}, stiff_3d::AbstractTensor{4,3,T}, args...) where {Nσ,T} + shear_factor = stiff_3d isa SymmetricTensor ? sqrt(2 * one(T)) : 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)) + f(c1,c2) = ((i,j) = c1; (k,l) = c2; stiff_3d[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)) end -# Extract each part of the stiffness tensor as SMatrix -# UniaxialStress -function extract_substiffnesses(stress_state::UniaxialStress, D::SymmetricTensor{4,3}) - # f=strain free, stress constrained - # c=strain constrained, stress free - ∂σᶠ∂ϵᶠ = get_unknowns(stress_state, D) - ∂σᶠ∂ϵᶜ = SMatrix{5,1}(D[2,2,1,1], D[3,3,1,1], √2*D[2,3,1,1], √2*D[1,3,1,1], √2*D[1,2,1,1]) - ∂σᶜ∂ϵᶠ = SMatrix{1,5}(D[1,1,2,2], D[1,1,3,3], √2*D[1,1,2,3], √2*D[1,1,1,3], √2*D[1,1,1,2]) - ∂σᶜ∂ϵᶜ = SMatrix{1,1}(D[1,1,1,1]) - return ∂σᶠ∂ϵᶠ, ∂σᶠ∂ϵᶜ, ∂σᶜ∂ϵᶠ, ∂σᶜ∂ϵᶜ -end -function extract_substiffnesses(stress_state::UniaxialStress, D::Tensor{4,3}) - # f=strain free, stress constrained - # c=strain constrained, stress free - ∂σᶠ∂ϵᶠ = get_unknowns(stress_state, D) - ∂σᶠ∂ϵᶜ = SMatrix{8,1}(D[2,2,1,1], D[3,3,1,1], D[2,3,1,1], D[1,3,1,1], D[1,2,1,1], D[3,2,1,1], D[3,1,1,1], D[2,1,1,1]) - ∂σᶜ∂ϵᶠ = SMatrix{1,8}(D[1,1,2,2], D[1,1,3,3], D[1,1,2,3], D[1,1,1,3], D[1,1,1,2], D[1,1,3,2], D[1,1,3,1], D[1,1,2,1]) - ∂σᶜ∂ϵᶜ = SMatrix{1,1}(D[1,1,1,1]) - return ∂σᶠ∂ϵᶠ, ∂σᶠ∂ϵᶜ, ∂σᶜ∂ϵᶠ, ∂σᶜ∂ϵᶜ -end - -# PlaneStress -function extract_substiffnesses(stress_state::PlaneStress, D::SymmetricTensor{4,3}) - # f=strain free, stress constrained - # c=strain constrained, stress free - ∂σᶠ∂ϵᶠ = get_unknowns(stress_state, D) - ∂σᶠ∂ϵᶜ = @SMatrix [ D[3,3,1,1] D[3,3,2,2] √2*D[3,3,2,1]; - √2*D[3,2,1,1] √2*D[3,2,2,2] 2*D[3,2,2,1]; - √2*D[3,1,1,1] √2*D[3,1,2,2] 2*D[3,1,2,1]] - - ∂σᶜ∂ϵᶠ = @SMatrix [ D[1,1,3,3] √2*D[1,1,3,2] √2*D[1,1,3,1]; - D[2,2,3,3] √2*D[2,2,3,2] √2*D[2,2,3,1]; - √2*D[2,1,3,3] 2*D[2,1,3,2] 2*D[2,1,3,1]] - - ∂σᶜ∂ϵᶜ = @SMatrix [ D[1,1,1,1] D[1,1,2,2] √2*D[1,1,2,1]; - D[2,2,1,1] D[2,2,2,2] √2*D[2,2,2,1]; - √2*D[2,1,1,1] √2*D[2,1,2,2] 2*D[2,1,2,1]] - - return ∂σᶠ∂ϵᶠ, ∂σᶠ∂ϵᶜ, ∂σᶜ∂ϵᶠ, ∂σᶜ∂ϵᶜ -end - -function extract_substiffnesses(stress_state::PlaneStress, D::Tensor{4,3}) - # f=strain free, stress constrained - # c=strain constrained, stress free - ∂σᶠ∂ϵᶠ = get_unknowns(stress_state, D) - ∂σᶠ∂ϵᶜ = @SMatrix [D[3,3,1,1] D[3,3,2,2] D[3,3,1,2] D[3,3,2,1]; - D[2,3,1,1] D[2,3,2,2] D[2,3,1,2] D[2,3,2,1]; - D[1,3,1,1] D[1,3,2,2] D[1,3,1,2] D[1,3,2,1]; - D[3,2,1,1] D[3,2,2,2] D[3,2,1,2] D[3,2,2,1]; - D[3,1,1,1] D[3,1,2,2] D[3,1,1,2] D[3,1,2,1]] - - ∂σᶜ∂ϵᶠ = @SMatrix [D[1,1,3,3] D[1,1,2,3] D[1,1,1,3] D[1,1,3,2] D[1,1,3,1]; - D[2,2,3,3] D[2,2,2,3] D[2,2,1,3] D[2,2,3,2] D[2,2,3,1]; - D[1,2,3,3] D[1,2,2,3] D[1,2,1,3] D[1,2,3,2] D[1,2,3,1]; - D[2,1,3,3] D[2,1,2,3] D[2,1,1,3] D[2,1,3,2] D[2,1,3,1]] - - ∂σᶜ∂ϵᶜ = @SMatrix [D[1,1,1,1] D[1,1,2,2] D[1,1,1,2] D[1,1,2,1]; - D[2,2,1,1] D[2,2,2,2] D[2,2,1,2] D[2,2,2,1]; - D[1,2,1,1] D[1,2,2,2] D[1,2,1,2] D[1,2,2,1]; - D[2,1,1,1] D[2,1,2,2] D[2,1,1,2] D[2,1,2,1]] - return ∂σᶠ∂ϵᶠ, ∂σᶠ∂ϵᶜ, ∂σᶜ∂ϵᶠ, ∂σᶜ∂ϵᶜ -end diff --git a/test/differentiation.jl b/test/differentiation.jl index c617bb7..d8dcf86 100644 --- a/test/differentiation.jl +++ b/test/differentiation.jl @@ -8,7 +8,11 @@ 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) + if TT == Tensor && isa(stress_state, UniaxialStress) + @test sc == 2:9 + else + #@test tomandel(ϵmock)[sc] ≈ MMB.get_unknowns(stress_state, ϵmock) + end end for stress_state in (UniaxialStrain(), PlaneStrain(), FullStressState()) sc = MMB.stress_controlled_indices(stress_state, ϵmock) diff --git a/test/errors.jl b/test/errors.jl index 94eadd5..99c610d 100644 --- a/test/errors.jl +++ b/test/errors.jl @@ -8,7 +8,7 @@ # Test that "impossible" stress iteration throws # the correct error m = LinearElastic(80e3, 150e3) - ss = PlaneStress(;max_iter = 1, tolerance=1e-100) + ss = PlaneStress(;maxiter = 1, tolerance=1e-100) ϵ = rand(SymmetricTensor{2,2}) @test_throws NoStressConvergence material_response(ss, m, ϵ, NoMaterialState{Float64}()) -end \ No newline at end of file +end diff --git a/test/runtests.jl b/test/runtests.jl index 49971fc..87ae3a5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,9 +5,12 @@ using Tensors, StaticArrays import MaterialModelsBase as MMB using FiniteDiff: FiniteDiff using MaterialModelsTesting: - LinearElastic, ViscoElastic, NeoHooke, test_derivative, obtain_numerical_material_derivative!, + LinearElastic, ViscoElastic, NeoHooke, + test_derivative, obtain_numerical_material_derivative!, runstrain, runstrain_diff, runstresstate, runstresstate_diff +import MechanicalMaterialModels as MechMat + include("utils4testing.jl") include("vector_conversion.jl") diff --git a/test/stressiterations.jl b/test/stressiterations.jl index bc474bd..7501548 100644 --- a/test/stressiterations.jl +++ b/test/stressiterations.jl @@ -12,9 +12,7 @@ iter_mandel = ( ([2,3,4,5,6,7,8,9],[2,3,4,5,6]), ([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 + # Test that all conversions are invertible (convert back and forth) for _stress_state in all_states for TT in (Tensor, SymmetricTensor) stress_state = if isa(_stress_state, NamedTuple) @@ -34,48 +32,21 @@ iter_mandel = ( ([2,3,4,5,6,7,8,9],[2,3,4,5,6]), @test ared ≈ MMB.reduce_tensordim(stress_state, afull) end end - - 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) - Ax = MMB.get_unknowns(state, A) - # Corresponding to regular mandel - @test ax ≈ tomandel(a)[ii] - @test Ax ≈ tomandel(A)[ii,ii] - - am = tomandel(a) - am[setdiff(1:length(am), ii)] .= 0 - ac = frommandel(typeof(a), am) - # Test that other comps become zero, and the rest are converted correctly - @test ac ≈ MMB.get_full_tensor(state, a, MMB.get_unknowns(state, a)) - end - end end -@testset "stiffness_calculations" begin +@testset "small strain reduced stiffness" begin 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 - D_mandel = tomandel(dσdϵ) - jj = setdiff(1:size(D_mandel,1), ii) - @test dσᶜdϵᶜ ≈ frommandel(TT{4,_getdim(state)}, inv(inv(D_mandel)[jj,jj])) - else - @test dσdϵ == dσᶜdϵᶜ - end + ii = inds[2] + TT = SymmetricTensor + state = isa(_stress_state, NamedTuple) ? _stress_state[:SymmetricTensor] : _stress_state + dσdϵ = rand(TT{4,3}) + one(TT{4,3}) + dσᶜdϵᶜ = MMB.reduce_stiffness(state, dσdϵ, zero(TT{2,3}), zero(TT{2,3})) + if _getdim(state) < 3 # Otherwise, conversion is direct + D_mandel = tomandel(dσdϵ) + jj = setdiff(1:size(D_mandel,1), ii) + @test dσᶜdϵᶜ ≈ frommandel(TT{4,_getdim(state)}, inv(inv(D_mandel)[jj,jj])) + else + @test dσdϵ == dσᶜdϵᶜ end end end @@ -172,6 +143,13 @@ end m = NeoHooke(;G, K) old = initial_material_state(m) # UniaxialStress + # Zero strain, F = I2 + P, dPdF, state, Ffull = material_response(UniaxialStress(), m, one(Tensor{2,1}), old, 0.0) + @test P[1,1] ≈ 0 atol = eps(E) + @test dPdF[1,1,1,1] ≈ E + @test Ffull ≈ one(Tensor{2,3}) + + # Small strain 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) @@ -186,6 +164,15 @@ end @test isapprox(P_full[2,2], λ*Δϵ; rtol) # PlaneStress + # * Zero strain + 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 tovoigt(symmetric(dPdF)) ≈ Dvoigt + @test norm(P) ≈ 0 atol = eps(E) + @test Ffull ≈ one(Tensor{2,3}) + + # * Small strain ϵ = Δϵ * rand(SymmetricTensor{2,2}) F = one(Tensor{2,2}) + ϵ P, dPdF, state, Ffull = material_response(PlaneStress(), m, F, old, 0.0) @@ -194,6 +181,28 @@ end @test isapprox(P, fromvoigt(SymmetricTensor{4,2}, Dvoigt)⊡ϵ; rtol) end +@testset "finite strain plasticity" begin + # This test causes failure without preventing rotation + # around x-axis in uniaxial stress + m = MechMat.FiniteStrainPlastic(; + elastic=MechMat.CompressibleNeoHooke(;G = 80.e3, K = 155e3), + yield=100.0, + isotropic=(MechMat.Voce(;Hiso=100e3, κ∞=50.0),), + kinematic=(MechMat.ArmstrongFrederick(;Hkin=80e3, β∞=130.0),), + overstress=MechMat.NortonOverstress(; tstar=0.5, nexp=2.0) + ) + old = initial_material_state(m) + ss = UniaxialStress() + Δϵ = 1e-5 + for ϵ in range(0, 0.2, 100)[2:end] + P1, dPdF1, state, Ffull = material_response(ss, m, Tensor{2,1}((1 + ϵ,)), old, 1e-3) + P2, dPdF2, _ = material_response(ss, m, Tensor{2,1}((1 + ϵ + Δϵ,)), old, 1e-3) + dPdF = (dPdF1 + dPdF2) / 2 + @test dPdF[1,1,1,1] ≈ (P2[1,1] - P1[1,1]) / Δϵ rtol = 1e-4 + old = state + end +end + @testset "visco_elastic" begin # Test a nonlinear, state and rate dependent, material by running stress iterations. # Check that state variables and time dependence are handled correctly