# [legacy fenics2019] 3D Elasto-plastic analysis

Hello, following this 2D elastoplastic example, I am trying to modify the code to a 3D problem, I believe I have set the dimensions correctly to full 3D but the Newton Raphson doesn’t converge. in the following code I try to point out things that I had changed (other than domain and BCs, and loading force), while the rest follows the tutorial. I have tried to do a test by setting the yield stress `sigma0` to a very large number(` = dl.Constant(1e9.)` for example), to keep the simulation in the elastic region and still having the same problem of Newton Raphson not converging, can anyone help spotting any issue in the code?

(update: the same mesh and BCs and load were testing in a separate hyperplasticity code and were fine, so I don’t think the problem would be in there)

``````##### domain and BCs stuff
import dolfin as dl
import numpy as np
import warnings
import matplotlib.pyplot as plt
longdim=1000.
shortdim=150.
mesh = dl.BoxMesh(dl.Point(-longdim/2., shortdim/2., -shortdim/2.), dl.Point(longdim/2., -shortdim/2., shortdim/2.),     10, 4, 4)
class LeftSurf(dl.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and dl.near(x[0],   -longdim/2.)
class RightSurf(dl.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and dl.near(x[0],   longdim/2.)
boundary_parts = dl.MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
boundary_parts.set_all(0) #marks whole cube as domain 0
left = LeftSurf()
left.mark(boundary_parts, 1)
right = RightSurf()
right.mark(boundary_parts, 2)
ds=dl.Measure('ds', domain=mesh,subdomain_data=boundary_parts)
``````
``````### following tutorial
E = dl.Constant(70e3)
nu = dl.Constant(0.3)
lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2./(1+nu)
sig0 = dl.Constant(250.)
Et = E/100.
H = E*Et/(E-Et)
deg_u = 2
deg_stress = 2
V = dl.VectorFunctionSpace(mesh, "CG", deg_u)
u = dl.Function(V, name="Total displacement")
du = dl.Function(V, name="Iteration correction")
Du = dl.Function(V, name="Current increment")
v = dl.TrialFunction(V)
u_ = dl.TestFunction(V)
``````
``````##### here the dimension of We is updated to 6 corresponding
#####to the 6 unique components in Voigt notation in 3D
W = dl.FunctionSpace(mesh, We)
``````
``````### following tutorial
sig = dl.Function(W)
sig_old = dl.Function(W)
n_elas = dl.Function(W)
W0 = dl.FunctionSpace(mesh, W0e)
beta = dl.Function(W0)
p = dl.Function(W0, name="Cumulative plastic strain")
``````
``````##### some more BCs and loading stuff
bc = [dl.DirichletBC(V, dl.Constant((0.,0.,0.)), LeftSurf()) ,  dl.DirichletBC(V, dl.Constant((0.,0.,0.)), RightSurf())]
def __init__(self,vl, t, **kwargs):
super().__init__(**kwargs)
self.value = vl
self.tol = 1.
self.t = t
def eval(self, values, x):
if dl.near(x[0],0.,self.tol) and dl.near(x[1],0.,self.tol) and dl.near(x[2],-shortdim/2.,self.tol):
values[0] = 0.
values[1] = 0.
values[2] = self.value  * self.t
else:
values[0] = 0.
values[1] = 0.
values[2] = 0.
def value_shape(self):
return (3,)
def F_ext(v):
``````
``````#### updated vector to tensor conversion with 6 Voigt components
def as_3D_tensor(X):
return dl.as_tensor([[X[0], X[3], X[4]],
[X[3], X[1], X[5]],
[X[4], X[5], X[2]]])
def sigma(eps_el):
return lmbda*dl.tr(eps_el)*dl.Identity(3) + 2*mu*eps_el
ppos = lambda x: (x+abs(x))/2.

#### output Voigt vectors are updated to 6 components
def proj_sig(deps, old_sig, old_p):
sig_n = as_3D_tensor(old_sig)
sig_elas = sig_n + sigma(deps)
s = dl.dev(sig_elas)
sig_eq = dl.sqrt(3/2.*dl.inner(s, s))
f_elas = sig_eq - sig0 - H*old_p
dp = ppos(f_elas)/(3*mu+H)
n_elas = s/sig_eq*ppos(f_elas)/f_elas
beta = 3*mu*dp/sig_eq
new_sig = sig_elas-beta*s
return dl.as_vector([new_sig[0, 0], new_sig[1, 1], new_sig[2, 2], new_sig[0, 1], new_sig[0, 2], new_sig[1, 2]]), \
dl.as_vector([n_elas[0, 0], n_elas[1, 1], n_elas[2, 2], n_elas[0, 1], n_elas[0, 2], n_elas[1, 2]]), \
beta, dp
``````
``````#### the rest of the code follows tutorial
def sigma_tang(e):
N_elas = as_3D_tensor(n_elas)
return sigma(e) - 3*mu*(3*mu/(3*mu+H)-beta)*dl.inner(N_elas, e)*N_elas-2*mu*beta*dl.dev(e)
def local_project(v, V, u=None):
dv = dl.TrialFunction(V)
v_ = dl.TestFunction(V)
a_proj = dl.inner(dv, v_)*dxm
b_proj = dl.inner(v, v_)*dxm
solver = dl.LocalSolver(a_proj, b_proj)
solver.factorize()
if u is None:
u = dl.Function(V)
solver.solve_local_rhs(u)
return u
else:
solver.solve_local_rhs(u)
return
file_results = dl.XDMFFile("plasticity_results.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
P0 = dl.FunctionSpace(mesh, "DG", 0)
p_avg = dl.Function(P0, name="Plastic strain")
Nitermax, tol = 200, 1e-8  # parameters of the Newton-Raphson procedure
Nincr = 20
results = np.zeros((Nincr+1, 2))
a_Newton = dl.inner(eps(v), sigma_tang(eps(u_)))*dxm
res = -dl.inner(eps(u_), as_3D_tensor(sig))*dxm + F_ext(u_)
A, Res = dl.assemble_system(a_Newton, res, bc)
nRes0 = Res.norm("l2")
nRes = nRes0
Du.interpolate(dl.Constant((0., 0., 0.))).  ##### updated to 3 components
print("Increment:", str(i+1), ', t = ', t)
niter = 0
while nRes/nRes0 > tol and niter < Nitermax:
dl.solve(A, du.vector(), Res, "mumps")
Du.assign(Du+du)
deps = eps(Du)
sig_, n_elas_, beta_, dp_ = proj_sig(deps, sig_old, p)
local_project(sig_, W, sig)
local_project(n_elas_, W, n_elas)
local_project(beta_, W0, beta)
A, Res = dl.assemble_system(a_Newton, res, bc)
nRes = Res.norm("l2")
print("Iter #", str(niter) ,"    Residual:", nRes)
niter += 1
u.assign(u+Du)
p.assign(p+local_project(dp_, W0))
sig_old.assign(sig)
file_results.write(u, t)
p_avg.assign(dl.project(p, P0))
file_results.write(p_avg, t)
results[i+1, :] = (u(0.,0., -shortdim/2.)[3], t)
``````

Please make an effort to simplify as much as you can the code (if it possible), and make sure to post it in a single code section ````...````, otherwise it’s impossible for others to copy and run it.

Hi Francesco, here is a cleaned up and put together version:

``````import dolfin as dl
import numpy as np
mesh = dl.BoxMesh(dl.Point(-500,75,-75), dl.Point(500, -75, 75,     10, 4, 4)
class LeftSurf(dl.SubDomain):
def inside(self, x, on_boundary):return on_boundary and dl.near(x[0],   -longdim/2.)
class RightSurf(dl.SubDomain):
def inside(self, x, on_boundary):return on_boundary and dl.near(x[0],   longdim/2.)
boundary_parts = dl.MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
boundary_parts.set_all(0)
LeftSurf().mark(boundary_parts, 1)
RightSurf().mark(boundary_parts, 2)
ds=dl.Measure('ds', domain=mesh,subdomain_data=boundary_parts)
E = dl.Constant(70e3)
nu = dl.Constant(0.3)
lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2./(1+nu)
sig0 = dl.Constant(250.)
Et = E/100.
H = E*Et/(E-Et)
deg_u = 2
deg_stress = 2
V = dl.VectorFunctionSpace(mesh, "CG", deg_u)
u = dl.Function(V)
du = dl.Function(V)
Du = dl.Function(V)
v = dl.TrialFunction(V)
u_ = dl.TestFunction(V)
sig = dl.Function(W)
sig_old = dl.Function(W)
n_elas = dl.Function(W)
W0 = dl.FunctionSpace(mesh, W0e)
beta = dl.Function(W0)
p = dl.Function(W0)
bc = [dl.DirichletBC(V, dl.Constant((0.,0.,0.)), LeftSurf()) ,  dl.DirichletBC(V, dl.Constant((0.,0.,0.)), RightSurf())]
def __init__(self,vl, t, **kwargs):
super().__init__(**kwargs)
self.value = vl
self.tol = 1.
self.t = t
def eval(self, values, x):
if dl.near(x[0],0.,self.tol) and dl.near(x[1],0.,self.tol) and dl.near(x[2],-75.,self.tol):
values[0] = 0.
values[1] = 0.
values[2] = self.value  * self.t
else:
values[0] = 0.
values[1] = 0.
values[2] = 0.
def value_shape(self):
return (3,)
def as_3D_tensor(X):
return dl.as_tensor([[X[0], X[3], X[4]],
[X[3], X[1], X[5]],
[X[4], X[5], X[2]]])
def sigma(eps_el):return lmbda*dl.tr(eps_el)*dl.Identity(3) + 2*mu*eps_el
ppos = lambda x: (x+abs(x))/2.
def proj_sig(deps, old_sig, old_p):
sig_n = as_3D_tensor(old_sig)
sig_elas = sig_n + sigma(deps)
s = dl.dev(sig_elas)
sig_eq = dl.sqrt(3/2.*dl.inner(s, s))
f_elas = sig_eq - sig0 - H*old_p
dp = ppos(f_elas)/(3*mu+H)
n_elas = s/sig_eq*ppos(f_elas)/f_elas
beta = 3*mu*dp/sig_eq
new_sig = sig_elas-beta*s
return dl.as_vector([new_sig[0, 0], new_sig[1, 1], new_sig[2, 2], new_sig[0, 1], new_sig[0, 2], new_sig[1, 2]]), \
dl.as_vector([n_elas[0, 0], n_elas[1, 1], n_elas[2, 2], n_elas[0, 1], n_elas[0, 2], n_elas[1, 2]]), \
beta, dp
def sigma_tang(e):
N_elas = as_3D_tensor(n_elas)
return sigma(e) - 3*mu*(3*mu/(3*mu+H)-beta)*dl.inner(N_elas, e)*N_elas-2*mu*beta*dl.dev(e)
def local_project(v, V, u=None):
dv = dl.TrialFunction(V)
v_ = dl.TestFunction(V)
a_proj = dl.inner(dv, v_)*dxm
b_proj = dl.inner(v, v_)*dxm
solver = dl.LocalSolver(a_proj, b_proj)
solver.factorize()
if u is None:
u = dl.Function(V)
solver.solve_local_rhs(u)
return u
else:
solver.solve_local_rhs(u)
return
Nitermax, tol = 200, 1e-8
a_Newton = dl.inner(eps(v), sigma_tang(eps(u_)))*dxm
res = -dl.inner(eps(u_), as_3D_tensor(sig))*dxm + F_ext(u_)
A, Res = dl.assemble_system(a_Newton, res, bc)
nRes0 = Res.norm("l2")
nRes = nRes0
Du.interpolate(dl.Constant((0., 0., 0.)))
niter = 0
while nRes/nRes0 > tol and niter < Nitermax:
dl.solve(A, du.vector(), Res, "mumps")
Du.assign(Du+du)
deps = eps(Du)
sig_, n_elas_, beta_, dp_ = proj_sig(deps, sig_old, p)
local_project(sig_, W, sig)
local_project(n_elas_, W, n_elas)
local_project(beta_, W0, beta)
A, Res = dl.assemble_system(a_Newton, res, bc)
nRes = Res.norm("l2")
print("Iter #", str(niter) ,"    Residual:", nRes)
niter += 1
u.assign(u+Du)
p.assign(p+local_project(dp_, W0))
sig_old.assign(sig)
``````