From 0e097d7753a73d26fb95b90035890160a3184205 Mon Sep 17 00:00:00 2001
From: Elias <elias.isak.borjesson@gmail.com>
Date: Fri, 9 Jul 2021 17:16:49 +0200
Subject: [PATCH 1/3] extra output and new material model

---
 src/FiniteStrain/largedef_plastic.jl | 187 +++++++++++++++++++++++++++
 src/MaterialModels.jl                |  18 +++
 2 files changed, 205 insertions(+)
 create mode 100644 src/FiniteStrain/largedef_plastic.jl

diff --git a/src/FiniteStrain/largedef_plastic.jl b/src/FiniteStrain/largedef_plastic.jl
new file mode 100644
index 0000000..c1f4604
--- /dev/null
+++ b/src/FiniteStrain/largedef_plastic.jl
@@ -0,0 +1,187 @@
+export MatHyperElasticPlastic, MatHyperElasticPlasticState
+
+"""
+    MatHyperElasticPlastic <: AbstractMaterial
+
+Hyper elastic plastic material for large deformations
+"""
+struct MatHyperElasticPlastic{M <: AbstractMaterial} <: AbstractMaterial
+    density::Float64
+    elastic_material::M #Elastic material, Yeoh or Neo-hook
+    τ₀::Float64		    #Yield stress
+    H::Float64		    #Hardening
+end
+
+struct MatHyperElasticPlasticState <: AbstractMaterialState
+    ϵᵖ::Float64					       #Plastic strain
+    Fᵖ::Tensor{2,3,Float64,9}		   #Plastic deformation grad
+    ν::Tensor{2,3,Float64,9}	
+end
+
+struct MatHyperElasticPlasticExtras <: AbstractExtras
+    D  ::Float64 #Dissipation
+    ∂D ::SymmetricTensor{2,3,Float64,6}
+end
+
+# # # # # # #
+# Constructors
+# # # # # # #
+
+function init_material_state(::MatHyperElasticPlastic)
+    Fᵖ = one(Tensor{2,3,Float64})
+    ν = zero(Tensor{2,3,Float64})
+    return MatHyperElasticPlasticState(0.0, Fᵖ, ν)
+end
+
+function MatHyperElasticPlastic(; 
+    elastic_material::HyperElasticMaterial,
+    τ₀  ::T,
+    H   ::T,
+    density ::T = NaN, 
+    ) where T
+
+    return MatHyperElasticPlastic(density, elastic_material, τ₀, H)
+end
+
+# # # # # # #
+# Drivers
+# # # # # # #
+
+function _compute_2nd_PK(mp::MatHyperElasticPlastic, C::SymmetricTensor{2,dim,T}, state::MatHyperElasticPlasticState, compute_dissipation::Bool) where {dim,T}
+    emat, τ₀, H = (mp.elastic_material, mp.τ₀, mp.H)
+    
+    I = one(C)
+    Iᵈᵉᵛ = otimesu(I,I) - 1/3*otimes(I,I)
+    
+    ϵᵖ = state.ϵᵖ
+    Fᵖ = state.Fᵖ
+    ν = state.ν
+    
+    dγ = 0.0
+    
+    phi, Mᵈᵉᵛ, Cᵉ, S̃, ∂S̃∂Cₑ = hyper_yield_function(C, state.Fᵖ, state.ϵᵖ, emat, τ₀, H)
+    dFᵖdC = zero(Tensor{4,dim,T})
+    dΔγdC = zero(Tensor{2,dim,T})
+    dgdC = zero(Tensor{2,dim,T})
+    g = 0.0
+
+    if phi > 0 #Plastic step
+        #Setup newton varibles
+        TOL = 1e-9
+        newton_error, counter = TOL + 1, 0
+        dγ = 0.0
+        
+        local ∂ϕ∂Mᵈᵉᵛ, ∂Mᵈᵉᵛ∂M, ∂M∂Cₑ, dϕdΔγ, ∂ϕ∂M, ∂M∂S̃
+        while true; counter +=1 
+
+            Fᵖ = state.Fᵖ - dγ*state.Fᵖ⋅state.ν
+            ϵᵖ = state.ϵᵖ + dγ
+
+            R, Mᵈᵉᵛ, Cₑ, S̃, ∂S̃∂Cₑ = _hyper_yield_function(C, Fᵖ, ϵᵖ, emat, τ₀, H) 
+            
+            #Numerical diff
+            #h = 1e-6
+            #J = (yield_function(C, state.Fᵖ - (dγ+h)*state.Fᵖ⋅state.ν, state.ϵᵖ + (dγ+h), μ, λ, τ₀, H)[1] - yield_function(C, state.Fᵖ - dγ*state.Fᵖ⋅state.ν, state.ϵᵖ + dγ, μ, λ, τ₀, H)[1]) / h
+
+            ∂ϕ∂M = sqrt(3/2) * (Mᵈᵉᵛ/norm(Mᵈᵉᵛ)) #⊡ Iᵈᵉᵛ #  state.ν#
+            ∂M∂Cₑ = otimesu(I, S̃)
+            ∂Cₑ∂Fᵖ = otimesl(I, transpose(Fᵖ)⋅C) + otimesu(transpose(Fᵖ)⋅C, I)
+            ∂Fᵖ∂Δγ = -state.Fᵖ⋅state.ν
+            ∂M∂S̃ = otimesu(Cₑ, I)
+
+            dϕdΔγ = ∂ϕ∂M ⊡ ((∂M∂Cₑ + ∂M∂S̃ ⊡ ∂S̃∂Cₑ) ⊡ ∂Cₑ∂Fᵖ ⊡ ∂Fᵖ∂Δγ) - H
+             
+            ddγ = dϕdΔγ\-R
+            dγ += ddγ
+
+            newton_error = norm(R)
+
+            if newton_error < TOL
+                break
+            end
+
+            if counter > 6
+                #@warn("Could not find equilibrium in material")
+                break
+            end
+        end
+        
+        ∂Cₑ∂C = otimesu(transpose(Fᵖ),transpose(Fᵖ))
+        ∂ϕ∂C = ∂ϕ∂M ⊡ ((∂M∂Cₑ + ∂M∂S̃ ⊡ ∂S̃∂Cₑ) ⊡ ∂Cₑ∂C) 
+        dΔγdC = -inv(dϕdΔγ) * ∂ϕ∂C
+        dFᵖdC = -(state.Fᵖ⋅state.ν) ⊗ dΔγdC
+
+    end
+
+    
+    ν = norm(Mᵈᵉᵛ)==0.0 ? zero(Mᵈᵉᵛ) : sqrt(3/2)*Mᵈᵉᵛ/(norm(Mᵈᵉᵛ))
+
+    S = (Fᵖ ⋅ S̃ ⋅ transpose(Fᵖ))
+
+    ∂S∂Fᵖ = otimesu(I, (Fᵖ ⋅ S̃)) + otimesl((Fᵖ ⋅ S̃), I)
+    ∂S∂S̃ = otimesu(Fᵖ,Fᵖ)
+    ∂Cₑ∂Fᵖ = otimesl(I, transpose(Fᵖ) ⋅ C) + otimesu(transpose(Fᵖ) ⋅ C, I)
+    ∂Cₑ∂C = otimesu(transpose(Fᵖ),transpose(Fᵖ))
+    dCₑdC = ∂Cₑ∂C + ∂Cₑ∂Fᵖ ⊡ dFᵖdC
+    dSdC = ∂S∂Fᵖ ⊡ dFᵖdC + ∂S∂S̃ ⊡ ∂S̃∂Cₑ ⊡ dCₑdC
+
+    g = 0.0, 
+    dgdC = zero(SymmetricTensor{2,3,T,6})
+    if compute_dissipation
+        if phi > 0.0
+            M = Cᵉ ⋅ S̃
+            g, dgdC = _compute_dissipation(M, ν, dγ, dCₑdC, ∂S̃∂Cₑ, dΔγdC, ∂M∂Cₑ, ∂M∂S̃, Iᵈᵉᵛ, I)
+        end
+    end
+
+    return symmetric(S), dSdC, ϵᵖ, ν, Fᵖ, g, dgdC
+end
+
+function _hyper_yield_function(C::SymmetricTensor, Fᵖ::Tensor, ϵᵖ::T, emat::MatNeoHook, τ₀::Float64, H::Float64) where T
+    Cᵉ = symmetric(transpose(Fᵖ)⋅C⋅Fᵖ)
+
+    S̃, ∂S̃∂Cₑ = _constitutive_driver(emat, Cᵉ)
+
+    #Mandel stress
+    Mbar = Cᵉ⋅S̃
+    Mᵈᵉᵛ = dev(Mbar)
+
+    yield_func = sqrt(3/2)*norm(Mᵈᵉᵛ) - (τ₀ + H*ϵᵖ)
+
+    return yield_func, Mᵈᵉᵛ, Cᵉ, S̃, ∂S̃∂Cₑ
+end
+
+function _compute_dissipation(M, ν, Δγ, dCₑdC, ∂S̃∂Cₑ, dΔγdC, ∂M∂Cₑ, ∂M∂S̃, Iᵈᵉᵛ, I)
+
+    Mᵈᵉᵛ = dev(M)
+
+    ∂g∂M = ν * Δγ
+    ∂g∂ν = M * Δγ
+    ∂g∂Δγ = M ⊡ ν
+    
+    ∂ν∂Mᵈᵉᵛ = √(3/2) * (1/norm(Mᵈᵉᵛ)) * (Iᵈᵉᵛ - (Mᵈᵉᵛ ⊗ Mᵈᵉᵛ)/(norm(Mᵈᵉᵛ)^2) )
+    ∂Mᵈᵉᵛ∂M = Iᵈᵉᵛ
+    dMdC = (∂M∂Cₑ + ∂M∂S̃ ⊡ ∂S̃∂Cₑ) ⊡ dCₑdC
+
+    g = M ⊡ ν*Δγ
+    dgdC = (∂g∂M + ∂g∂ν ⊡ ∂ν∂Mᵈᵉᵛ ⊡ ∂Mᵈᵉᵛ∂M) ⊡ dMdC  +  ∂g∂Δγ * dΔγdC 
+
+    return g, dgdC
+end
+
+function material_response(mp::MatHyperElasticPlastic, C::SymmetricTensor{2,3,T,6}, state::MatHyperElasticPlasticState, Δt=0.0
+    ; nothing, options::Dict{Symbol, Any} = Dict{Symbol, Any}()) 
+
+    S, ∂S∂C, ϵᵖ, ν, Fᵖ, _, _ = _compute_2nd_PK(mp, C, state, false)
+
+    return S, ∂S∂C, MatHyperElasticPlasticState(ϵᵖ, Fᵖ, ν)
+end
+
+
+function material_response(mp::MatHyperElasticPlastic, C::SymmetricTensor{2,3,T,6}, state::MatHyperElasticPlasticState, Δt=0.0, ::Symbol
+    ; nothing, ::Dict{Symbol, Any} = Dict{Symbol, Any}()) 
+    
+    S, ∂S∂C, ϵᵖ, ν, Fᵖ, g, dgdC = _compute_2nd_PK(mp, C, state, true)
+
+    return S, ∂S∂C, MatHyperElasticPlasticState(ϵᵖ, Fᵖ, ν), MatHyperElasticPlasticExtras(g, dgdC)
+end
diff --git a/src/MaterialModels.jl b/src/MaterialModels.jl
index ca375f5..35b42c9 100644
--- a/src/MaterialModels.jl
+++ b/src/MaterialModels.jl
@@ -27,6 +27,14 @@ Store state variables here. For now, this should **not** be mutable, a new objec
 abstract type AbstractMaterialState end
 abstract type AbstractResiduals end
 
+"""
+    AbstractExtras
+
+Outputs from the material routine in addition to what is usual or strictly necessary (stress, tangent and state)
+"""
+abstract type AbstractExtras end
+struct EmptyExtras <: AbstractExtras end
+
 """
     material_response(m::AbstractMaterial, Δε::SymmetricTensor{2,3}, state::AbstractMaterialState, Δt; cache, options)
 
@@ -40,6 +48,16 @@ This function signature must be the same for all material models, even if they d
 """
 function material_response end
 
+"""
+    material_response(m::AbstractMaterial, Δε::SymmetricTensor{2,3}, state::AbstractMaterialState, Δt, ::Symbol; cache, options)
+
+Fallback for materials that has no extra outputs
+"""
+function material_response(m::AbstractMaterial, strain, state::AbstractMaterialState, Δt, ::Symbol; cache, options)
+    stress, tangent, state = material_response(m, strain, state, Δt; cache=cache, options=options)
+    return stress, tangent, state, EmptyExtras()
+end
+
 """
     initial_material_state(::AbstractMaterial)
 

From a2d298ba91e12c08d3f6179fe1712f078a428a4e Mon Sep 17 00:00:00 2001
From: Elias <elias.isak.borjesson@gmail.com>
Date: Tue, 24 Aug 2021 19:26:32 +0200
Subject: [PATCH 2/3] add some docs and tests

---
 src/FiniteStrain/largedef_plastic.jl         |  84 ++++++++++++++-----
 src/MaterialModels.jl                        |   2 +
 test/jld2_files/MatHyperElasticPlastic1.jld2 | Bin 0 -> 22783 bytes
 test/runtests.jl                             |   3 +-
 test/test_largedef_plastic.jl                |  73 ++++++++++++++++
 5 files changed, 140 insertions(+), 22 deletions(-)
 create mode 100644 test/jld2_files/MatHyperElasticPlastic1.jld2
 create mode 100644 test/test_largedef_plastic.jl

diff --git a/src/FiniteStrain/largedef_plastic.jl b/src/FiniteStrain/largedef_plastic.jl
index c1f4604..912d27f 100644
--- a/src/FiniteStrain/largedef_plastic.jl
+++ b/src/FiniteStrain/largedef_plastic.jl
@@ -1,15 +1,23 @@
 export MatHyperElasticPlastic, MatHyperElasticPlasticState
 
 """
-    MatHyperElasticPlastic <: AbstractMaterial
+    MatHyperElasticPlastic(elastic_material, σ_y, H)
 
-Hyper elastic plastic material for large deformations
+Large strain plasticity model with kinematic hardening.
+
+# Arguments
+- `elastic_material::AbstractMaterial`: Hyperelastic material, for exampel NeoHook
+- `σ_y:` yield limit
+- `H:` hardening modulus
+
+# Reference
+J.C. Simo, 1998, Numerical analysis and simulation of plasticity 
 """
 struct MatHyperElasticPlastic{M <: AbstractMaterial} <: AbstractMaterial
-    density::Float64
     elastic_material::M #Elastic material, Yeoh or Neo-hook
-    τ₀::Float64		    #Yield stress
+    σ_y::Float64		#Yield stress
     H::Float64		    #Hardening
+    density::Float64
 end
 
 struct MatHyperElasticPlasticState <: AbstractMaterialState
@@ -27,20 +35,20 @@ end
 # Constructors
 # # # # # # #
 
-function init_material_state(::MatHyperElasticPlastic)
+function initial_material_state(::MatHyperElasticPlastic)
     Fᵖ = one(Tensor{2,3,Float64})
     ν = zero(Tensor{2,3,Float64})
     return MatHyperElasticPlasticState(0.0, Fᵖ, ν)
 end
 
 function MatHyperElasticPlastic(; 
-    elastic_material::HyperElasticMaterial,
-    τ₀  ::T,
+    elastic_material::AbstractMaterial,
+    σ_y  ::T,
     H   ::T,
     density ::T = NaN, 
     ) where T
 
-    return MatHyperElasticPlastic(density, elastic_material, τ₀, H)
+    return MatHyperElasticPlastic(elastic_material, σ_y, H, density)
 end
 
 # # # # # # #
@@ -48,7 +56,7 @@ end
 # # # # # # #
 
 function _compute_2nd_PK(mp::MatHyperElasticPlastic, C::SymmetricTensor{2,dim,T}, state::MatHyperElasticPlasticState, compute_dissipation::Bool) where {dim,T}
-    emat, τ₀, H = (mp.elastic_material, mp.τ₀, mp.H)
+    emat, σ_y, H = (mp.elastic_material, mp.σ_y, mp.H)
     
     I = one(C)
     Iᵈᵉᵛ = otimesu(I,I) - 1/3*otimes(I,I)
@@ -59,7 +67,7 @@ function _compute_2nd_PK(mp::MatHyperElasticPlastic, C::SymmetricTensor{2,dim,T}
     
     dγ = 0.0
     
-    phi, Mᵈᵉᵛ, Cᵉ, S̃, ∂S̃∂Cₑ = hyper_yield_function(C, state.Fᵖ, state.ϵᵖ, emat, τ₀, H)
+    phi, Mᵈᵉᵛ, Cᵉ, S̃, ∂S̃∂Cₑ = _hyper_yield_function(C, state.Fᵖ, state.ϵᵖ, emat, σ_y, H)
     dFᵖdC = zero(Tensor{4,dim,T})
     dΔγdC = zero(Tensor{2,dim,T})
     dgdC = zero(Tensor{2,dim,T})
@@ -77,11 +85,11 @@ function _compute_2nd_PK(mp::MatHyperElasticPlastic, C::SymmetricTensor{2,dim,T}
             Fᵖ = state.Fᵖ - dγ*state.Fᵖ⋅state.ν
             ϵᵖ = state.ϵᵖ + dγ
 
-            R, Mᵈᵉᵛ, Cₑ, S̃, ∂S̃∂Cₑ = _hyper_yield_function(C, Fᵖ, ϵᵖ, emat, τ₀, H) 
+            R, Mᵈᵉᵛ, Cₑ, S̃, ∂S̃∂Cₑ = _hyper_yield_function(C, Fᵖ, ϵᵖ, emat, σ_y, H) 
             
             #Numerical diff
             #h = 1e-6
-            #J = (yield_function(C, state.Fᵖ - (dγ+h)*state.Fᵖ⋅state.ν, state.ϵᵖ + (dγ+h), μ, λ, τ₀, H)[1] - yield_function(C, state.Fᵖ - dγ*state.Fᵖ⋅state.ν, state.ϵᵖ + dγ, μ, λ, τ₀, H)[1]) / h
+            #J = (yield_function(C, state.Fᵖ - (dγ+h)*state.Fᵖ⋅state.ν, state.ϵᵖ + (dγ+h), μ, λ, σ_y, H)[1] - yield_function(C, state.Fᵖ - dγ*state.Fᵖ⋅state.ν, state.ϵᵖ + dγ, μ, λ, σ_y, H)[1]) / h
 
             ∂ϕ∂M = sqrt(3/2) * (Mᵈᵉᵛ/norm(Mᵈᵉᵛ)) #⊡ Iᵈᵉᵛ #  state.ν#
             ∂M∂Cₑ = otimesu(I, S̃)
@@ -115,7 +123,6 @@ function _compute_2nd_PK(mp::MatHyperElasticPlastic, C::SymmetricTensor{2,dim,T}
 
     
     ν = norm(Mᵈᵉᵛ)==0.0 ? zero(Mᵈᵉᵛ) : sqrt(3/2)*Mᵈᵉᵛ/(norm(Mᵈᵉᵛ))
-
     S = (Fᵖ ⋅ S̃ ⋅ transpose(Fᵖ))
 
     ∂S∂Fᵖ = otimesu(I, (Fᵖ ⋅ S̃)) + otimesl((Fᵖ ⋅ S̃), I)
@@ -125,7 +132,7 @@ function _compute_2nd_PK(mp::MatHyperElasticPlastic, C::SymmetricTensor{2,dim,T}
     dCₑdC = ∂Cₑ∂C + ∂Cₑ∂Fᵖ ⊡ dFᵖdC
     dSdC = ∂S∂Fᵖ ⊡ dFᵖdC + ∂S∂S̃ ⊡ ∂S̃∂Cₑ ⊡ dCₑdC
 
-    g = 0.0, 
+    g = 0.0 
     dgdC = zero(SymmetricTensor{2,3,T,6})
     if compute_dissipation
         if phi > 0.0
@@ -137,16 +144,16 @@ function _compute_2nd_PK(mp::MatHyperElasticPlastic, C::SymmetricTensor{2,dim,T}
     return symmetric(S), dSdC, ϵᵖ, ν, Fᵖ, g, dgdC
 end
 
-function _hyper_yield_function(C::SymmetricTensor, Fᵖ::Tensor, ϵᵖ::T, emat::MatNeoHook, τ₀::Float64, H::Float64) where T
+function _hyper_yield_function(C::SymmetricTensor, Fᵖ::Tensor, ϵᵖ::T, emat::AbstractMaterial, σ_y::Float64, H::Float64) where T
     Cᵉ = symmetric(transpose(Fᵖ)⋅C⋅Fᵖ)
 
-    S̃, ∂S̃∂Cₑ = _constitutive_driver(emat, Cᵉ)
+    S̃, ∂S̃∂Cₑ = material_response(emat, Cᵉ)
 
     #Mandel stress
     Mbar = Cᵉ⋅S̃
     Mᵈᵉᵛ = dev(Mbar)
 
-    yield_func = sqrt(3/2)*norm(Mᵈᵉᵛ) - (τ₀ + H*ϵᵖ)
+    yield_func = sqrt(3/2)*norm(Mᵈᵉᵛ) - (σ_y + H*ϵᵖ)
 
     return yield_func, Mᵈᵉᵛ, Cᵉ, S̃, ∂S̃∂Cₑ
 end
@@ -169,17 +176,52 @@ function _compute_dissipation(M, ν, Δγ, dCₑdC, ∂S̃∂Cₑ, dΔγdC, ∂M
     return g, dgdC
 end
 
-function material_response(mp::MatHyperElasticPlastic, C::SymmetricTensor{2,3,T,6}, state::MatHyperElasticPlasticState, Δt=0.0
-    ; nothing, options::Dict{Symbol, Any} = Dict{Symbol, Any}()) 
+"""
+    material_response(mp::MatHyperElasticPlastic, C::SymmetricTensor{2,3,T,6}, state::MatHyperElasticPlasticState, Δt=0.0;  <keyword arguments>)
+
+The material model assumes a multiplicative split of the deformation gradient
+
+```math
+\\boldsymbol{F} = \\boldsymbol{F}_e \\boldsymbol{F}_p
+```
+The yield function is of von Mises type is formulated in terms of the Mandel stress, \$\\boldsymbol{M}\$
+
+```math
+\\phi = \\sqrt{\\frac{3}{2}} \\left| \\boldsymbol{M}_{\\text{dev}} \\right| - \\sigma_y - H \\varepsilon_p
+```
+
+The evolution equation of associative type for the plastic velocity gradient is
+
+```math
+\\boldsymbol{L}_p = \\dot \\boldsymbol{F}_p \\boldsymbol{F}_p^{-1} = \\dot \\gamma \\frac{\\partial \\phi}{\\partial \\boldsymbol{M}} = \\dot \\gamma \\boldsymbol{\\nu}
+```
+
+As an simplification for the local backward Euler integration in the material routine, it is assumed that
+\$ \\boldsymbol{\\nu} = ^n\\boldsymbol{\\nu} \$. This reduces the system of equations to only one equation 
+(namely the yield function).
+    
+"""
+function material_response(mp::MatHyperElasticPlastic, C::SymmetricTensor{2,3,T,6}, state::MatHyperElasticPlasticState, Δt=0.0; 
+    cache=nothing, options=nothing) where T
 
     S, ∂S∂C, ϵᵖ, ν, Fᵖ, _, _ = _compute_2nd_PK(mp, C, state, false)
 
     return S, ∂S∂C, MatHyperElasticPlasticState(ϵᵖ, Fᵖ, ν)
 end
 
+"""
+    material_response(mp::MatHyperElasticPlastic, C::SymmetricTensor{2,3,T,6}, state::MatHyperElasticPlasticState, Δt, extras::Symbol;  <keyword arguments>)
+
+In addition to returning the stress, tangent and state, it also return extra data related to the amount 
+of dissipation energy, according to:
 
-function material_response(mp::MatHyperElasticPlastic, C::SymmetricTensor{2,3,T,6}, state::MatHyperElasticPlasticState, Δt=0.0, ::Symbol
-    ; nothing, ::Dict{Symbol, Any} = Dict{Symbol, Any}()) 
+```math
+D = \\dot \\boldsymbol{M} : \\boldsymbol{L}_p (\\geq 0)
+```
+
+"""
+function material_response(mp::MatHyperElasticPlastic, C::SymmetricTensor{2,3,T,6}, state::MatHyperElasticPlasticState, Δt, ::Symbol
+    ; cache=nothing, options=nothing) where T 
     
     S, ∂S∂C, ϵᵖ, ν, Fᵖ, g, dgdC = _compute_2nd_PK(mp, C, state, true)
 
diff --git a/src/MaterialModels.jl b/src/MaterialModels.jl
index 35b42c9..e4945de 100644
--- a/src/MaterialModels.jl
+++ b/src/MaterialModels.jl
@@ -86,6 +86,7 @@ function update_cache! end
 
 include("LinearElastic.jl")
 include("Plastic.jl")
+include("FiniteStrain/largedef_plastic.jl")
 include("CrystalViscoPlastic/slipsystems.jl")
 include("CrystalViscoPlastic/CrystalViscoPlastic.jl")
 include("CrystalViscoPlastic/CrystalViscoPlasticRed.jl")
@@ -103,6 +104,7 @@ export material_response
 export AbstractMaterial
 export LinearElastic, Plastic
 export LinearElasticState, PlasticState
+export MatHyperElasticPlastic, MatHyperElasticPlasticState
 export OneD, UniaxialStrain, UniaxialStress, PlaneStrain, PlaneStress
 
 export NeoHook, Yeoh, StVenant
diff --git a/test/jld2_files/MatHyperElasticPlastic1.jld2 b/test/jld2_files/MatHyperElasticPlastic1.jld2
new file mode 100644
index 0000000000000000000000000000000000000000..fce11b55f061c3181ad2d9649b938e9a6c824cae
GIT binary patch
literal 22783
zcmeI430RcX8pjU<%7B0<CgM8eo`N$lGk_Z(hKfqL5Q%0gC?<-6fP&_RWx1qT?xtyJ
zin&{;n9I1AW}4a!HIv5cE~$abEth*fzVn}<!<k@i-RIt(XTC>y`M<w&-t(UKJ0Bd*
zz(hrMG6W>VrzZ~9bRC_X6t9Vl&xqG_N=X|LpV3k?CNV8NDP^Qa8>kD^5lu6ti7wC>
zs0~OPq|+Gn0SQSNn&^(riR1cPR$zKmWG7Ge+U_c~nh=#Tg{QdU6y5J-#Xh{&eQEco
z$X+V)wsU(T8EJD3B_rcUBoY!<&53wZrKSI!=zl`gBsD&bvJAS2h<6zE$ULNk5aJO*
z$8)Me+=&YrCa=5!8DmouiL*)$AYXDap7j5hvY?j|RQVhqjZ563mCrjr-#xmdg**nm
zoLGA~f_y!nl9m`KukDxQQRL`IvIHwUV&qt7=66@zivOk2smXFDdbu#}p-6>yN=}K-
zFzRJjCB`v}%hi>J6}6RFe_s~=rNoixDQW3}F=I#2N~R?ZQXIsUwbQV6@imlBmN1L)
zPBX1$?QPA=jqz!04WPfnDs^?#;eS&HceXCQ*gD7az#aoqg?kQ+EtQ5T%;LV~FwMe{
z{3`g&VI384RKQUIM+F=ea8$ri0Y?QK6>wDG?_L4roFcLk{c+aW&B;YMPH=W7&ej74
zbtC!5c?~_MupTn7pMHNgd(W9>cI<JSg`+?GtLOnkRo2X-sT{rv4fMOC44ibxG4t86
zgfsm)9Den9dHgKm@r)8@II~BYZC+FHPh<8fvz^%4z$eVQ)Kc0FVs<~X{(g#oGP4(%
zZC_jQA3k#ME2T{5o!_U&IhB>e8yz7!O?XLF2$6`s3ei7xCM1ZcMebnYEpqEgEs@(m
zG$OZ=G!wZibVMW{XyAy5rb*~5EV>f?ymrQs;dU$rdTTA+>8dGn^^ebfu_T<6zmxL{
zS~Yr%B%k#RrbkPJ{7q{pFC^4ZuAe+{Xjl5V2R+(S(bGFM@gPHz5|alHkf$r3*QHy|
z^W+CYIZ`9!X*Jr3UY$putD)3Bkb4k4h?AB+ck~~<B%-EeccBW8-`i@_4|a8w@4S3{
z9<aTNT^;4=%2&*)v5%CFFL9}se9h4ND9zNg5-ue-T+tQ%xO?S>UUq~(T0eOGj1vQ;
z(|K$2n{4-zFihiJ;YauUBif?J#VoIJi0V~d2`%jiKg9Lwfy85}lA+JumP6Opw}cyI
zIy<o7Oh1bb$a`@67S+>E<o2~A{DPKyKVNz_Q##W%bikft4K3kfI|t`XxIM(8!yjJQ
zGLvdk>XTJ=gfH$ew_ugi1S$Ee+`iRzSi^61s?_kxr>U0kOA-oq-K9Fla5m%rC;X@v
zu5226CQG^=cfbGA{05e6|G77HL)NYl79A3o=Q@MxeAUr3JF?e4*?gkY>I`Y*#t%MQ
z->#M=T<2H*m2&;rVHW-K`7PV-QC(-&4~2Gwe?c30w@PNJ6dC_8rPp$AOStV@7f!7|
z-Fj<d9<-a(gz6b*ld9Pf{@oYW&U|myP-#V_o-1no=41)K^;o+rzPA%By6eI3`qZWR
z%X8P*x8@2xwg2Vtr@!fcvx85zlrS%2^}1!&g{jKSxjyy6XiN6TBW{xURA&ypJ=>1r
z6i-Q9KgKjoS~R%L?TKxy;rBHgJUMO06ifKZT_3$Mh3bc`H0+CLg^F|OPL=nMwVERZ
zEYDrL__Jp%;i_Iue*N&I85VtVx6je@R7W&$S#3x5)o*349Pn_ilyf2W#{J6H@R#Ec
z%pCCSdzSF~(-ZC_Q@w2Vy<1B7|IV)Kr<SDTTHH#pJ%beMW>D<D5XFvdtyo<<#a<6r
z>_;6G+oz*qt3@gHtFDU8enGKeJrsMlmtx0vY&)@JUQRl-^EFq=*pm?XMHjm5vbFCL
z`sv0?=%+g`p`UKOgnqjB68b}kR_G5Up+bKc{WdPz#ObsoN(^+QhZqErI57w&1I0j3
zhKhlKq>6!&ZoGx*hLFi(5K5+rK^U1M1VLJ|NDOpjg%|{p95Dzco5es+c8Gz2n8d(H
zw|1hL2ZfN6Vh~F5#UP9n2tlxxToVHwxg!Qa<e?Y@)7NrE=?3eGS_}-tTMUe(mKcN(
zjTnTIW?~RVv_hcQl29?w5lIY!NR$`^lOAHACvjq6AOpp~NQR0*2uT%#P?9AEVPvup
z7_?-X80g3xF$f}y#2}cg5Cc8Q5d#C+ECxohLkvQQNen{Cmtqh`P6~lhOY+4)M+(Fs
zh+GqcU~)$c^yHx!&~H~N(eq(6615nF5N|OECAGvLjA(=)L`#~9fsSa!Ac%yDK`@cT
zKu@B?z(9J4fsw?CK?oTr2BBoA7=)2jAqdryEHTiL$wFZFVTSy(G2I(jeh_z~pJ%TQ
zouF(UveeEG=;pym**qwx55%#Tkc%cCDB#pe=h#c|B?)ruCHRs!_7awpI6OFd$6kV8
zB?iY{f-i~Dv6tXW67v7DmymyB?#q&IhlKu;?e7ou<zIfY_fUUy9VOVs##%4sV5WkH
z`tqv@^0lV=N_JpW!H)m&p}xFrlCO_9&Fidem)5I9r}|qPD(%$osOYKwG*yi@TYSD4
zNg@<)*^l|tm*r1yaZ4&o6c5j7xHemkm?ec0!fW^|YRl_yFE6h>pBFzrp8whU@$vZl
zcwJt5dpX~p*ZlarF0Vf2an$eStZ)46S+^s4Mht&@Fo_tG7<gvk`PwGr4gJw~`kJwg
zY_HpzhckR+j`ksICS-NlYnBiPdUn6l$ysJ=d~jILHoKmD<>_jIhF!mK_vq5`;8$MW
zda@TB8vL)y2lMCihgQd!atf1eKlg`nkDV1g_k6y{`}w5SgK8c$vzPCl?1cvBE#}Yr
z5g#0yZwGCQV|;LE^z&u+-;29{(6;7}{`Kz;SR>bOws)=e_N?3DdY=Zp(w9VhzUZ^+
zg-4$=A)ZBdJ5sk{J=^PcT*2{bGDrJ7eV2mE+sqQ;K&yNFuxxd{H9k1(lH<Qsy6oy}
zf`-)}?fzZqc<?K)n|hhRp~3&MeAxWmxcs5-n%;k1Gu-yNeaj{Gf@t0H`6BOQ<9<_(
zH<{V%!d@n5aNc76ydUwwq4{>uwm8NIheki#|Ha(>gSIt)^sn32XN+9G?mN$%vuE9Q
zd2QYPBdv?2-^To`@BPmWCd7Nhbm6|qrw6qwmf+{@lj|N$Wxl@Zh~uX(Z7?0am2mjN
z46~*E5vE%`-u`8uq-7RAAIId?qNlX5EDx&e)c>uG)_CZrzA~!a!Qz@GaA@$2<>aF+
zC!c*RpQxLSYm8}MEai@MJxWrn`E6s`x5d?=_3|D|7?-!Y*E4-HA6_oQ`Fzo@(Cwzi
zyvr<sLxb}c^XL8O4-U<@gSN#nJ~%Y`ITyIZr_=iEycNFW!cuE}b}?P|>nruf3g(ag
zfj5>F)u;9QvHt!7d)DozPP_6?Jy$H9%UM3Y^=0e2J<oK5jt7qK39-FylT=x=c2fS1
z*J~4IT-jg>t<vbwdsEC3@(N>mPK|k0XGL1)1rGb_&oBDS&o0ZG6b^0LVzpU9J6P`z
z&W(IJtG)>w8hkd(Cxpuf`OPgSpLmu}Ve^3eIQm>SP{sG%-rbs?hUw##t8d<&`?ySg
z$m7GD8e0nm_pKVsAMIcR{e8cW$}>yg(BQnq{CPj>0S?W#gT{E!uoxd48oYUfdxq|`
zK6P$xoc+l=*80?AI&j^CYDI$iqrY>j1_KY$`dvwlu<PsXtfoJD`E?=@wNjq??ONN~
zgm}J8e?Fq>q$6Io_X+Ne2gb<UkF`H-s1;Lsot<I2{k4y_c^|UI$2hPjepr6a|4|bY
zG%V^44U78#4h@ccPI37lzl!FA`Sb0eE4vN2`_^l=_X(%I9d%Q*Pry1rzOYyq;LzZ_
z#r%0c>H!YTw}Zxb(6AUE92y+=pTB?5w&suieomv#%k@hdvbT;sUvDQ|{vy9gHxl9T
zlJV68uY{No4_eBP?)PeS+w1njh?A{l?#cS^(bmj}`_wEU4)oH)Z#~FOwZ;dB?Y1*{
zZmR~3Owh1(mn`m8Iv)JBpDcgA?f%*(aA@#9n@=>$r#+iLpFea?{*n$OhT2}Y?_YDf
zDq6RGEPvz+dt$_;$<KUYmcXIGd5ih;e$)dTnr{b<@t|QbJ~%Wu?mvJ3pl!_`{WGQ?
z*U9y(YFvKNo^?Cp(yDRB9wcJNp6kUstLROLm&f#wIzQbgdgc$-ZQJluGDrWQjdeaO
zy&u&v9ov0a-G+6n@i7i;z{h*P$d2&0j1T{?6Jw{AjtBn_IlqjM_qX8C;D0usmMkC4
z|79+J%;$-&?aG)W+v|4OnvhGPb<5|Ayr(pu5+(0%!J)x<i}~|@^aqFL+d<pn7#|!O
z{g$))&)+}j4{gjJ{bP>(v_P(3uP)c_+Ouv4+#PtYb1aE?H>hCmM|oNk;>}??;Pqdh
zo9SnJ-TL^R^OHIH$E?WuHg~I8LLBI|T|4CUK4^^(4jVPCT^sG;x+ZAYu(olg((&Mr
zQu#+OZ0utKhX((v^1=N1{6A&$Y4dVNpVq@{uiJ#E-!6*QEuSy)j_vy(VACG61P%?(
zTg;#Lqdzz_-wxUq$N1pT=!g5?nY(|`w&sui(K(yk<@!yX+D&<-!>MAgJH&?X=pvtk
zeelWH(9-8%(CTu|!7|x77~1!(Il=2nsdKPCOi$jHd*F}G!4`7oD45s&b#FMAjt76e
zCTUqi*Ls=2p~081`TW^@+}ZrIx%{Dn-`?^5ADx3`+c*cS!}3SIuw7j{>*RAVaA<Je
zV*b1z^#F(F+d*SIXjqI74h@d`zks`c(6;7}{=HRUv*h}H)nNEqd%oT})-5XddyD7D
zv2OqPIr?9%+v@+?nQ*sGu{7Vg_r*>>ThGCu*WBGv<R5JN9PFLl^)r1azxwFs2fvs%
zm{fh;7e|gWOIQWa&&&@}6=|&JB;c?OYaDgGom`d=n{D(uv&<}^9c-r)?k^g=8=AnO
z!I6(Mmk)GB^NBRBUb@3jEO{3WyRa?Pnjf_5td;Mc-}1Oje#j&6!%2oj!8s70FWRp;
zkmj;E*DQfUgYy>i=l$pp4$ZfNw#6|%I5hgTBd0!IOY763bw<kOH?8%7emv>qdeQvR
zf6do%H?GtA4T(%+e+yEveZs&OG%3EfilhZ@iBrG3YF)PrneJfppZb%(?RER8^0I>C
zlzYw`-ah&I29x^Cu?2JE%@Xo*WBSf3HCG01T59ow!@jdAuCcb6b>8q#=$zj1+M8wx
z8rCx|WcB3UfhKTh3;&DqaoZmM&6?)LQg~ib@9$4p^Ou+|P<uT6WsU81J7N6fpeDbT
zeci!}<&S(}i|(|1{yII60fz?XE#}YrQ4ervz8y5igNDWU;Lza1>bJc#k=AFh%Rg?X
zb+g{TAxuBr`s1KBg88F=o0IK&dFMQq_GvO>?Aa&W@6~RApQ%WCv~=A5{iR-afbQ+`
zb&r4=w%6^&VS|r9^(mGngg0C_pT6$!+~H>vR9(yx@`5%j>C@i-MeBVBhux)3kB%-X
z%lj4gZ{eS8me3A%i^qF&!%cb<I5haga`G{jlh1CJPs07LZ(X6^k79qUIxjlSnjdt7
zY4`j0TxxsW-dj<>f41Q3EuSy)R`-2%Lz5L|2^<=nx0pZgM}Kf=z8$nJj`6{v(eM1#
z_GAOC&%mrvp>2Jw>lV7vef^SIg88HW0l&JLZ_xUkk8ZQso^{*jU_t7Xc|}s=;nC}I
z==XWVgMRbd^zYxOY<t~K&v1?FNcoV3O@}O?uk)(vYML(4nkB@6{;XSm_PN3jEPim<
zh}bAk^B(KEfdAc&`zB6|F-y>}buaZfo7b_m2^<<6`8ad=K$n+K2bRyT2b^yhyox2Q
z_o@#2s#)`cZc@1KPVOAr>o)D#dD-9JFZ=bD&lh?3D$WTXI^QgTLxb}c^XL8O4-U<@
zgSN#nJ~%Y`U0K{~^L5W+De>fK|F_N7`an1Bx2@7;!Tiy`)x*4#4{7}#6sZoCtXq<!
zsxCvkK`X2~{R2|^))(}>KApUs)LP}De|f1~oIeO}&%pVE@g@hHUytAX<opKw`Xpg*
zZo-BuF|}{rXMcRW=nqq-?@2;a_Ds)6OH5BsOdl{Td2mqEYAkOCGvY_SmN+s)fy-($
dNDD!GRxWwl*z~orY;alv;U!4%MsS_({|2;vI!pin

literal 0
HcmV?d00001

diff --git a/test/runtests.jl b/test/runtests.jl
index 62d1536..806bc29 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -12,4 +12,5 @@ include("test_crystal_visco_plastic_red.jl")
 include("test_wrappers.jl")
 include("test_neohook.jl")
 include("test_yeoh.jl")
-include("test_stvenant.jl")
\ No newline at end of file
+include("test_stvenant.jl")
+include("test_largedef_plastic.jl")
\ No newline at end of file
diff --git a/test/test_largedef_plastic.jl b/test/test_largedef_plastic.jl
new file mode 100644
index 0000000..54d7812
--- /dev/null
+++ b/test/test_largedef_plastic.jl
@@ -0,0 +1,73 @@
+@testset "MatHyperElasticPlastic" begin
+    # constructor
+    E=200e3
+    ν=0.3
+    neohook = NeoHook(
+        μ = E / (2(1 + ν)),
+        λ = (E * ν) / ((1 + ν) * (1 - 2ν))
+    )
+
+    m = MatHyperElasticPlastic(
+        elastic_material = neohook,
+        σ_y = 200.0, 
+        H = 50.0,
+    )
+
+    # initial state
+    state = initial_material_state(m)
+
+    # elastic branch
+    ε11_yield = m.σ_y / E 
+    Δt = 0.0
+
+    F = Tensor{2,3}((0.5ε11_yield + 1, 0.0, 0.0, 0.0, (-0.5ε11_yield*ν + 1.0), 0.0, 0.0, 0.0, (-0.5ε11_yield*ν + 1.0)))
+    C = tdot(F)
+    S_neo, ∂S_neo, _ = material_response(neohook, C)
+    S, ∂S, temp_state = material_response(m, C, state)
+    _,_,_,extras = material_response(m, C, state, Δt , :extras)
+
+    @test S_neo ≈ S
+    @test ∂S == ∂S_neo
+    @test temp_state.ϵᵖ == 0.0
+    @test temp_state.Fᵖ == one(Tensor{2,3,Float64})
+    @test extras.D == 0.0 # No dissipation
+
+    # Plastic branch
+    F = Tensor{2,3}((1.1ε11_yield + 1, 0.0, 0.0, 0.0, (-1.1ε11_yield*ν + 1.0), 0.0, 0.0, 0.0, (-1.1ε11_yield*ν + 1.0)))
+    C = tdot(F)
+
+    S, ∂S, temp_state, extras = material_response(m, C, state, Δt, :extras)
+
+    @test temp_state.ϵᵖ > 0.0
+    @test extras.D > 0.0
+
+end
+
+
+function get_MatHyperElasticPlastic_loading()
+    strain1 = range(0.0,  0.002, length=5)
+    strain2 = range(0.002, 0.001, length=5)
+    strain3 = range(0.001, 0.007, length=5)
+
+    _x = [strain1[1:end]..., strain2[1:end]..., strain3[1:end]...]
+    _F = [Tensor{2,3}((x + 1.0, x/10, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0)) for x in _x]
+    C = [tdot(F) for F in _F]
+
+    return C
+end  
+
+@testset "MatHyperElasticPlastic jld2" begin
+
+    E=200e3
+    ν=0.3
+    m = MatHyperElasticPlastic(
+        elastic_material = MaterialModels.NeoHook(
+            μ = E / (2(1 + ν)),
+            λ = (E * ν) / ((1 + ν) * (1 - 2ν))
+        ),
+        σ_y = 200.0, 
+        H = 50.0,
+    )
+    loading = get_MatHyperElasticPlastic_loading()
+    check_jld2(m, loading, "MatHyperElasticPlastic1")#, debug_print=true, OVERWRITE_JLD2=true)
+end
\ No newline at end of file

From b1d29d57d8da16f92573ae8c4bc9ecf3a5b5525d Mon Sep 17 00:00:00 2001
From: Elias <elias.isak.borjesson@gmail.com>
Date: Wed, 25 Aug 2021 12:53:03 +0200
Subject: [PATCH 3/3] add docs

---
 docs/make.jl                          | 3 ++-
 docs/src/materials/LargeDefPlastic.md | 6 ++++++
 2 files changed, 8 insertions(+), 1 deletion(-)
 create mode 100644 docs/src/materials/LargeDefPlastic.md

diff --git a/docs/make.jl b/docs/make.jl
index 28b4a39..b4baa9f 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -8,7 +8,8 @@ makedocs(
         "Interface" => "interface.md",
         "Materials" => [
             "materials/LinearElastic.md",
-            "materials/Plastic.md"
+            "materials/Plastic.md",
+            "materials/LargeDefPlastic.md"
         ]
     ]
 )
diff --git a/docs/src/materials/LargeDefPlastic.md b/docs/src/materials/LargeDefPlastic.md
new file mode 100644
index 0000000..6a22a10
--- /dev/null
+++ b/docs/src/materials/LargeDefPlastic.md
@@ -0,0 +1,6 @@
+# LargeDefPlastic
+```@docs
+MatHyperElasticPlastic
+material_response(mp::MatHyperElasticPlastic, C::SymmetricTensor{2, 3, T, 6}, state::MatHyperElasticPlasticState) where T
+material_response(mp::MatHyperElasticPlastic, C::SymmetricTensor{2, 3, T, 6}, state::MatHyperElasticPlasticState, Δt, ::Symbol; kwargs...) where T
+```
\ No newline at end of file