question about implementation of elastoplasticity simulation

[RWTHLHK]

I'd like to implement my elastoplasticity constitutive model and perform simulations using Gridap. This must be done by a nonlinear solver incrementally. I just checked the tutorial about low level API in possion equation, and I can define the mechanical behavior for each quadrature point with low-level API like in the UMAT subroutine in ABAQUS, right? Or is there some high-level API so I can access the incrementation of solution values in each substep?

[amartinhuertas]

Did you take a look at this tutorial? https://gridap.github.io/Tutorials/dev/pages/t010_isotropic_damage/#isotropic_damage.jl-1

I think it might be helpful for your purpose.

[RWTHLHK]

Thanks a lot for your reply! I tried to implement a linear isotropic hardening model. But one major difference between this isotropic damage and isotropic hardening is when defining stress, not only the strain_in is required,, but also the dstrain_in. Is this ok? The following are my codes:
`
const E = 3.0e10 # Pa
const v = 0.3
const lamda = (Ev)/(1+v)(1-2v)
const mu = E/(2(1+v))
const sigma0 = 4.0e8 #Pa
const ytol = 1e-9 # yield tolerence
const h = 3.0e5 #Pa
sigmae(epsilon) = lamda * tr(epsilon) * one(epsilon) + 2.0muepsilon

function dev(m) # compute deviatoric stress
m - (tr(m)/3.0)*one(m)
end

function mises_stress(sigma_dev) #compute mises stress
sqrt(3.0/2.0 * sum(sigma_dev.*sigma_dev))
end

function hard_func(p) #linear isotropic hardening
sigma0 + h*p
end

function hard_derivative(p)
h
end
function check_yield(mises, p)
mises - hard_func(p)
end

function new_state(epsilon, deps, p_in)
sigma_trial = simgae(epsilon + deps)
sigma_tdev = dev(sigma_trial) # compute deviatoric trial stress
mises_sigmat = mises_stress(sigma_tdev) # compute trial mises stress
flow_n = sigma_tdev / mises_sigmat # flow direction
if check_yield(mises_sigmat, p_in) < sigma0*ytol
yield = false
delp = 0
p_out = p_in
else
yield = true
delp = solve_p(mises_sigmat, p_in) #newton raphson method
p_out = p_in + delp
end
yield, delp, p_out, flow_n
end

function sigma(epsion_in, deps, p_in)
_,delp,p_out,flow_n = new_state(epsion_in, deps, p_in)
depsp = 3.0/2.0 * delp * flow_n # plastic strain Incremental
depse = deps - depsp #elastic strain increment
sigmae(epsion_in + depse)
end

function dsigma(epsilon, deps, state)
yield, delp, p_out, flow_n = state
if yield
depsp = 3.0/2.0 * delp * flow_n # plastic strain Incremental
depse = deps - depsp #elastic strain increment
sigmae(depse)
else
sigmae(deps)
end
end

function solve_p(mises, p_in)
delp = 0
kewton = 1000
epoch = 0
yield_residual(delp) = mises - hard_func(p_in+delp) - 3mudelp
yield_jac(delp) = -3mu - h
res = yield_residual(delp)
jac = yield_jac(delp)
while res > sigma0ytol && epoch < kewton
delp -= res/jac
res = yield_residual(delp)
jac = yield_jac(delp)
epoch += 1
end
delp
`

[RWTHLHK]

Hi I just figured out how to implement my constitutive model. But to do that I need to store a tensor state variable. How can I initiliaze a tensor state variable using CellState