Error Converting Elasto-plastic tutorial to 3D

General
1

Hi Everyone

I’m trying to convert the 2D Elasto-plastic tutorial Here for 3D models. However I’m getting the following error:

Generating mesh with CGAL 3D mesh generator
Index out of bounds.
Traceback (most recent call last):
File “elastoplasticity.py”, line 140, in
a_Newton = inner(eps(v), sigma_tang(eps(u_)))*dxm
File “elastoplasticity.py”, line 72, in sigma_tang
N_elas = vector2tensor(n_elas)
File “elastoplasticity.py”, line 49, in vector2tensor
return as_tensor([[X[0], X[3], X[5]], [X[3], X[1], X[4]], [X[5], X[4], X[2]]])
File “/home/numericalmining/anaconda3/envs/fenics2018/lib/python3.7/site-packages/ufl/exproperators.py”, line 449, in _getitem
all_indices, slice_indices, repeated_indices = create_slice_indices(component, shape, self.ufl_free_indices)
File “/home/numericalmining/anaconda3/envs/fenics2018/lib/python3.7/site-packages/ufl/index_combination_utils.py”, line 164, in create_slice_indices
error(“Index out of bounds.”)
File “/home/numericalmining/anaconda3/envs/fenics2018/lib/python3.7/site-packages/ufl/log.py”, line 172, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Index out of bounds.

Here is the full code:

from dolfin import *
from mshr import *
from fenics import *
from dolfin import *
import numpy as np
parameters["form_compiler"]["representation"] = 'quadrature'
import warnings
from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning
warnings.simplefilter("once", QuadratureRepresentationDeprecationWarning)

# Geometria y Meshing
L=20 #ancho de la base
z1=5 #alto base
w1=4 #ancho pilar direccion x
w2=4 #ancho pilar direccion y
h=5 #altura pilar

piso=Box(Point(0,0,0),Point(L,L,z1))
techo=Box(Point(0,0,z1+h),Point(L,L,2*z1+h))
pilar=Box(Point(0.5*L-0.5*w1,0.5*L-0.5*w2,z1),Point(0.5*L+0.5*w1,0.5*L+0.5*w2,z1+h))

g3d=piso+pilar+techo

mesh= generate_mesh(g3d,25)

#Valores de Entrada

P=Constant((0.,0.,-10.))
E=Constant(40000.)
nu = Constant(0.3)
sig0 = Constant(250.)  # yield strength
Et = E/100.  # tangent modulus
H = E*Et/(E-Et)  # hardening modulus
mu = E/2/(1+nu)
lmbda =E*nu/(1+nu)/(1-2*nu)
dim=3


def eps(v):
    return sym(grad(v))

def sigma(eps_el):
    return lmbda*tr(eps_el)*Identity(dim) + 2*mu*eps_el

def vector2tensor(X):
    if dim is 2:
        return as_tensor([[X[0], X[3]], [X[3], X[1]]])
    else:
        return as_tensor([[X[0], X[3], X[5]], [X[3], X[1], X[4]], [X[5], X[4], X[2]]])

def tensor2vector(X):
    if dim is 2:
        return as_vector([X[0, 0], X[1, 1], 0, X[0, 1]])
    else:
        return as_vector([X[0, 0], X[1, 1], X[2, 2], X[0, 1], X[1, 2], X[0, 2]])

ppos = lambda x: (x+abs(x))/2.

def proj_sig(deps, old_sig, old_p):
    sig_n = vector2tensor(old_sig)
    sig_elas = sig_n + sigma(deps)
    s = dev(sig_elas)
    sig_eq = sqrt(3/2.*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 tensor2vector(new_sigma),tensor2vector(nelas),beta, dp

def sigma_tang(e):
    N_elas = vector2tensor(n_elas)
    return sigma(e) - 3*mu*(3*mu/(3*mu+H)-beta)*inner(N_elas, e)*N_elas-2*mu*beta*dev(e)


deg_u = 2
deg_stress = 2
V = VectorFunctionSpace(mesh, "CG", deg_u)
We = VectorElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, dim=4, quad_scheme='default')
W = FunctionSpace(mesh, We)
W0e = FiniteElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, quad_scheme='default')
W0 = FunctionSpace(mesh, W0e)

H=2*z1+h
#Condiciones de Borde Exterior
def bottom(x, on_boundary):
  return near(x[2], 0)

class superior(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[2], H) and on_boundary

def norte(x, on_boundary):
  return near(x[1], L)

def sur(x, on_boundary):
  return near(x[1], 0)

def este(x, on_boundary):
  return near(x[0], L)

def oeste(x, on_boundary):
  return near(x[0], 0)

#Condiciones de borde de Dirichlet
bc1 = DirichletBC(V, Constant((0.,0.,0.)), bottom)
bc2 = DirichletBC(V.sub(1), Constant(0.), norte)
bc3 = DirichletBC(V.sub(1), Constant(0.), sur)
bc4 = DirichletBC(V.sub(0), Constant(0.), este)
bc5 = DirichletBC(V.sub(0), Constant(0.), oeste)
bc=[bc1,bc2,bc3,bc4,bc5]

#Condiciones de Borde de Newmann
facets = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
facets.set_all(0)
superior().mark(facets,1)
ds = Measure('ds', subdomain_data=facets)
q=10
loading = Expression("-q*t", q=q, t=0, degree=2)
n = FacetNormal(mesh)
def F_ext(v):
    return loading*dot(n, v)*ds(1)

sig = Function(W)
sig_old = Function(W)
n_elas = Function(W)
beta = Function(W0)
p = Function(W0, name="Cumulative plastic strain")
u = Function(V, name="Total displacement")
du = Function(V, name="Iteration correction")
Du = Function(V, name="Current increment")
v = TrialFunction(V)
u_ = TestFunction(V)
#
#

metadata = {"quadrature_degree": deg_stress, "quadrature_scheme": "default"}
dxm = dx(metadata=metadata)

a_Newton = inner(eps(v), sigma_tang(eps(u_)))*dxm
res = -inner(eps(u_), as_3D_tensor(sig))*dxm + F_ext(u_)

def local_project(v, V, u=None):
    dv = TrialFunction(V)
    v_ = TestFunction(V)
    a_proj = inner(dv, v_)*dxm
    b_proj = inner(v, v_)*dxm
    solver = LocalSolver(a_proj, b_proj)
    solver.factorize()
    if u is None:
        u = Function(V)
        solver.solve_local_rhs(u)
        return u
    else:
        solver.solve_local_rhs(u)
        return

file_results = XDMFFile("plasticity_results.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
P0 = FunctionSpace(mesh, "DG", 0)
p_avg = Function(P0, name="Plastic strain")

Nitermax, tol = 5, 1e-8  # parameters of the Newton-Raphson procedure
Nincr = 20
load_steps = np.linspace(0, 1.1, Nincr+1)[1:]**0.5
results = np.zeros((Nincr+1, 2))
for (i, t) in enumerate(load_steps):
    loading.t = t
    A, Res = assemble_system(a_Newton, res, bc)
    nRes0 = Res.norm("l2")
    nRes = nRes0
    Du.interpolate(Constant((0., 0.,0.)))
    print("Increment:", str(i+1))
    niter = 0
    while nRes/nRes0 > tol and niter < Nitermax:
        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 = assemble_system(a_Newton, res, bc)
        nRes = Res.norm("l2")
        print("    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(project(p, P0))
file_results.write(p_avg, t)

Some idea how to solve it?

Thank

2

As you have provided far from a minimal example, I can just give you some pointers:

eps(v)

and

eps(u_)

are tensors of shape (3, 3)

n_elas is a (4, ) vector.
Your vector2tensor conversion assumes that n_elas is alot longer:

i.e. that it has 6 entries.

Note that for further questions, I would recommend you to remove code not needed to reproduce the error message, see below for an example

from dolfin import *

mesh = UnitCubeMesh(10, 5, 3)

E=Constant(40000.)
nu = Constant(0.3)
Et = E/100.  # tangent modulus
H = E*Et/(E-Et)  # hardening modulus
mu = E/2/(1+nu)
lmbda =E*nu/(1+nu)/(1-2*nu)
dim=3


def eps(v):
    return sym(grad(v))

def sigma(eps_el):
    return lmbda*tr(eps_el)*Identity(dim) + 2*mu*eps_el

def vector2tensor(X):
    if dim is 2:
        return as_tensor([[X[0], X[3]], [X[3], X[1]]])
    else:
        return as_tensor([[X[0], X[3], X[5]], [X[3], X[1], X[4]], [X[5], X[4], X[2]]])

def tensor2vector(X):
    if dim is 2:
        return as_vector([X[0, 0], X[1, 1], 0, X[0, 1]])
    else:
        return as_vector([X[0, 0], X[1, 1], X[2, 2], X[0, 1], X[1, 2], X[0, 2]])

ppos = lambda x: (x+abs(x))/2.

def sigma_tang(e):
    N_elas = vector2tensor(n_elas)
    return sigma(e) - 3*mu*(3*mu/(3*mu+H)-beta)*inner(N_elas, e)*N_elas-2*mu*beta*dev(e)


deg_u = 2
deg_stress = 2
V = VectorFunctionSpace(mesh, "CG", deg_u)
We = VectorElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, dim=4, quad_scheme='default')
W = FunctionSpace(mesh, We)
W0e = FiniteElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, quad_scheme='default')
W0 = FunctionSpace(mesh, W0e)

sig = Function(W)
sig_old = Function(W)
n_elas = Function(W)
beta = Function(W0)
p = Function(W0, name="Cumulative plastic strain")
u = Function(V, name="Total displacement")
du = Function(V, name="Iteration correction")
Du = Function(V, name="Current increment")
v = TrialFunction(V)
u_ = TestFunction(V)
#
#

metadata = {"quadrature_degree": deg_stress, "quadrature_scheme": "default"}
dxm = dx(metadata=metadata)
a_Newton = inner(eps(v), sigma_tang(eps(u_)))*dxm