tensors_material.jl

# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE

using Tensors
import Parameters

Parameters.@with_kw struct IdealPlasticMaterialParameters{dim}
    # Material parameters
    E :: Float64
    ν :: Float64
    σy :: Float64
    λ ::Float64 = E*ν / ((1+ν) * (1 - 2ν))
    μ ::Float64 = E / (2(1+ν))
end

mutable struct IdealPlasticMaterialState{dim, M}
    # Internal variables
    σ :: SymmetricTensor{2, dim, Float64, M}
    ϵp :: SymmetricTensor{2, dim, Float64, M}
end

function IdealPlasticMaterialState{dim}() where {dim}
    ϵp = zero(SymmetricTensor{2, dim})
    σ = zero(SymmetricTensor{2, dim})
    return IdealPlasticMaterialState(ϵp, σ)
end

# Could cache This
function compute_elastic_stiffness(mp::IdealPlasticMaterialParameters{dim}) where {dim}
    δ(i,j) = i == j ? 1.0 : 0.0
    g(i,j,k,l) = mp.λ*δ(i,j)*δ(k,l) + mp.μ*(δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k))
    return SymmetricTensor{4, dim}(g)
end

# Returns incremental stress and consistent tangent stiffness
function integrate_material(mp::IdealPlasticMaterialParameters{dim}, ms::IdealPlasticMaterialState{dim},
                            Δϵ::SymmetricTensor{2, dim}, Δt::Number) where {dim}

    D = compute_elastic_stiffness(mp)
    
    Δσ_tr = D ⊡ Δϵ

    σ_tr = ms.σ + Δσ_tr
    σ_tr_dev = dev(σ_tr)
    σ_tr_v = √(3/2) * norm(σ_tr_dev)
    if σ_tr_v <= mp.σy
        return Δσ_tr, D
    else
        n = 3/2 * σ_tr_dev/σ_tr_v
        Δλ = 1/(3*mp.μ)*(σ_tr_v - mp.σy)
        Δϵ_p = Δλ * n
        Δσ = D ⊡ (Δϵ - Δϵ_p)
        Dn = (D ⊡ n)
        tangent_stiffness = Dn ⊗ Dn / (n ⊡ Dn)
        return Δσ, tangent_stiffness
    end
end

function bench_ideal_plastic()
    mp = IdealPlasticMaterialParameters{3}(E=200e9, ν=0.3, σy = 100e6)
    ms = IdealPlasticMaterialState{3}()
    Δt = 0.25
    Δϵ = 0.25 * 1e-5 *  diagm(SymmetricTensor{2, 3}, [-0.3, -0.3, 1.0])
    σ_devs = Float64[]
    for i in 1:1000
        Δσ, D = integrate_material(mp, ms, Δϵ, Δt)
        ms.σ += Δσ
        push!(σ_devs, norm(dev(ms.σ)))
    end
    return σ_devs
end

using BenchmarkTools
@btime bench_ideal_plastic()
s = bench_ideal_plastic()

#=
using UnicodePlots
lineplot(s)
            ┌────────────────────────────────────────┐
   90000000 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
            │⠀⠀⠀⠀⠀⠀⠀⢰⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒│
            │⠀⠀⠀⠀⠀⠀⠀⡎⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
            │⠀⠀⠀⠀⠀⠀⣸⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
            │⠀⠀⠀⠀⠀⢀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
            │⠀⠀⠀⠀⠀⡼⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
            │⠀⠀⠀⠀⢠⠃⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
            │⠀⠀⠀⠀⡞⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
            │⠀⠀⠀⢰⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
            │⠀⠀⢀⡏⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
            │⠀⠀⡸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
            │⠀⢠⠇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
            │⠀⡞⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
            │⢰⠃⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
          0 │⡏⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
            └────────────────────────────────────────┘
            0                                     1000
=#