Regarding the FEniCSx demo: Elasto-plastic analysis of a 2D von Mises material - Elasto-plastic analysis of a 2D von Mises material — Computational Mechanics Numerical Tours with FEniCSx, I have two questions:
- In the legacy demo (Elasto-plastic analysis of a 2D von Mises material — Numerical tours of continuum mechanics using FEniCS master documentation), it is said the plastic strain (in quadrature function space) must first be projected onto the previously defined DG FunctionSpace:
p_avg.assign(project(p, P0)); file_results.write(p_avg, t), how to do this in FEniCSx? - Following item1, I defined a project function, but when trying to save the vector function in VTK file, the kernel died. please see the following code.
import numpy as np
import matplotlib.pyplot as plt
import gmsh
from mpi4py import MPI
import ufl
import basix
from dolfinx import mesh, fem, io
import dolfinx.fem.petsc
from petsc4py import PETSc
hsize = 0.2
Re = 1.3
Ri = 1.0
gmsh.initialize()
gdim = 2
model_rank = 0
if MPI.COMM_WORLD.rank == 0:
gmsh.option.setNumber("General.Terminal", 0) # to disable meshing info
gmsh.model.add("Model")
geom = gmsh.model.geo
center = geom.add_point(0, 0, 0)
p1 = geom.add_point(Ri, 0, 0)
p2 = geom.add_point(Re, 0, 0)
p3 = geom.add_point(0, Re, 0)
p4 = geom.add_point(0, Ri, 0)
x_radius = geom.add_line(p1, p2)
outer_circ = geom.add_circle_arc(p2, center, p3)
y_radius = geom.add_line(p3, p4)
inner_circ = geom.add_circle_arc(p4, center, p1)
boundary = geom.add_curve_loop([x_radius, outer_circ, y_radius, inner_circ])
surf = geom.add_plane_surface([boundary])
geom.synchronize()
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", hsize)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", hsize)
gmsh.model.addPhysicalGroup(gdim, [surf], 1)
gmsh.model.addPhysicalGroup(gdim - 1, [x_radius], 1, name="bottom")
gmsh.model.addPhysicalGroup(gdim - 1, [y_radius], 2, name="left")
gmsh.model.addPhysicalGroup(gdim - 1, [inner_circ], 3, name="inner")
gmsh.model.mesh.generate(gdim)
domain, _, facets = io.gmshio.model_to_mesh(
gmsh.model, MPI.COMM_WORLD, model_rank, gdim=gdim
)
gmsh.finalize()
E = fem.Constant(domain, 70e3) # in MPa
nu = fem.Constant(domain, 0.3)
lmbda = E * nu / (1 + nu) / (1 - 2 * nu)
mu = E / 2.0 / (1 + nu)
sig0 = fem.Constant(domain, 250.0) # yield strength in MPa
Et = E / 100.0 # tangent modulus
H = E * Et / (E - Et) # hardening modulus
deg_u = 2
shape = (gdim,)
V = fem.functionspace(domain, ("P", deg_u, shape))
Vx, _ = V.sub(0).collapse()
Vy, _ = V.sub(1).collapse()
bottom_dofsy = fem.locate_dofs_topological((V.sub(1), Vy), gdim - 1, facets.find(1))
top_dofsx = fem.locate_dofs_topological((V.sub(0), Vx), gdim - 1, facets.find(2))
# used for post-processing
def bottom_inside(x):
return np.logical_and(np.isclose(x[0], Ri), np.isclose(x[1], 0))
bottom_inside_dof = fem.locate_dofs_geometrical((V.sub(0), Vx), bottom_inside)[0]
u0x = fem.Function(Vx)
u0y = fem.Function(Vy)
bcs = [
fem.dirichletbc(u0x, top_dofsx, V.sub(0)),
fem.dirichletbc(u0y, bottom_dofsy, V.sub(1)),
]
n = ufl.FacetNormal(domain)
q_lim = float(2 / np.sqrt(3) * np.log(Re / Ri) * sig0)
loading = fem.Constant(domain, 0.0)
deg_quad = 2 # quadrature degree for internal state variable representation
W0e = basix.ufl.quadrature_element(
domain.basix_cell(), value_shape=(), scheme="default", degree=deg_quad
)
We = basix.ufl.quadrature_element(
domain.basix_cell(), value_shape=(4,), scheme="default", degree=deg_quad
)
W = fem.functionspace(domain, We)
W0 = fem.functionspace(domain, W0e)
sig = fem.Function(W)
sig_old = fem.Function(W)
n_elas = fem.Function(W)
beta = fem.Function(W0)
p = fem.Function(W0, name="Cumulative_plastic_strain")
dp = fem.Function(W0)
u = fem.Function(V, name="Total_displacement")
du = fem.Function(V, name="Iteration_correction")
Du = fem.Function(V, name="Current_increment")
v = ufl.TrialFunction(V)
u_ = ufl.TestFunction(V)
P0 = fem.functionspace(domain, ("DG", 0))
p_avg = fem.Function(P0, name="Plastic_strain")
def eps(v):
e = ufl.sym(ufl.grad(v))
return ufl.as_tensor([[e[0, 0], e[0, 1], 0], [e[0, 1], e[1, 1], 0], [0, 0, 0]])
def elastic_behavior(eps_el):
return lmbda * ufl.tr(eps_el) * ufl.Identity(3) + 2 * mu * eps_el
def as_3D_tensor(X):
return ufl.as_tensor([[X[0], X[3], 0], [X[3], X[1], 0], [0, 0, X[2]]])
def to_vect(X):
return ufl.as_vector([X[0, 0], X[1, 1], X[2, 2], X[0, 1]])
ppos = lambda x: ufl.max_value(x, 0)
def constitutive_update(Δε, old_sig, old_p):
sig_n = as_3D_tensor(old_sig)
sig_elas = sig_n + elastic_behavior(Δε)
s = ufl.dev(sig_elas)
sig_eq = ufl.sqrt(3 / 2.0 * ufl.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 to_vect(new_sig), to_vect(n_elas), beta, dp
def sigma_tang(eps):
N_elas = as_3D_tensor(n_elas)
return (
elastic_behavior(eps)
- 3 * mu * (3 * mu / (3 * mu + H) - beta) * ufl.inner(N_elas, eps) * N_elas
- 2 * mu * beta * ufl.dev(eps)
)
ds = ufl.Measure("ds", domain=domain, subdomain_data=facets)
dx = ufl.Measure(
"dx",
domain=domain,
metadata={"quadrature_degree": deg_quad, "quadrature_scheme": "default"},
)
Residual = ufl.inner(eps(u_), as_3D_tensor(sig)) * dx - ufl.inner(
-loading * n, u_
) * ds(3)
tangent_form = ufl.inner(eps(v), sigma_tang(eps(u_))) * dx
basix_celltype = getattr(basix.CellType, domain.topology.cell_type.name)
quadrature_points, weights = basix.make_quadrature(basix_celltype, deg_quad)
map_c = domain.topology.index_map(domain.topology.dim)
num_cells = map_c.size_local + map_c.num_ghosts
cells = np.arange(0, num_cells, dtype=np.int32)
def interpolate_quadrature(ufl_expr, function):
expr_expr = fem.Expression(ufl_expr, quadrature_points)
expr_eval = expr_expr.eval(domain, cells)
function.x.array[:] = expr_eval.flatten()[:]
class CustomLinearProblem(fem.petsc.LinearProblem):
def assemble_rhs(self, u=None):
"""Assemble right-hand side and lift Dirichlet bcs.
Parameters
----------
u : dolfinx.fem.Function, optional
For non-zero Dirichlet bcs u_D, use this function to assemble rhs with the value u_D - u_{bc}
where u_{bc} is the value of the given u at the corresponding. Typically used for custom Newton methods
with non-zero Dirichlet bcs.
"""
# Assemble rhs
with self._b.localForm() as b_loc:
b_loc.set(0)
fem.petsc.assemble_vector(self._b, self._L)
# Apply boundary conditions to the rhs
x0 = [] if u is None else [u.vector]
fem.petsc.apply_lifting(self._b, [self._a], bcs=[self.bcs], x0=x0, scale=1.0)
self._b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
x0 = None if u is None else u.vector
fem.petsc.set_bc(self._b, self.bcs, x0, scale=1.0)
def assemble_lhs(self):
self._A.zeroEntries()
fem.petsc.assemble_matrix_mat(self._A, self._a, bcs=self.bcs)
self._A.assemble()
def solve_system(self):
# Solve linear system and update ghost values in the solution
self._solver.solve(self._b, self._x)
self.u.x.scatter_forward()
tangent_problem = CustomLinearProblem(
tangent_form,
-Residual,
u=du,
bcs=bcs,
petsc_options={
"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps",
},
)
def project(source_function, target_space):
"""
Project a function from its current function space onto a target function space.
Parameters:
- source_function: fem.Function
The function to be projected.
- target_space: fem.FunctionSpace
The function space onto which the source function should be projected.
Returns:
- fem.Function
The projected function in the target function space.
"""
# Create the projected function in the target space
projected_function = fem.Function(target_space)
# Define the test function and trial function in the target space
v = ufl.TestFunction(target_space)
u = ufl.TrialFunction(target_space)
# Define the variational problem (a is the bilinear form, L is the linear form)
a = ufl.inner(u, v) * ufl.dx
L = ufl.inner(source_function, v) * ufl.dx
# Solve the variational problem to project the source function onto the target space
problem = fem.petsc.LinearProblem(a, L, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
projected_function = problem.solve()
return projected_function
# alternative
def project1(v, V):
"""
Project a function or expression `v` onto a function space `V`.
Parameters:
v : dolfinx.fem.Function or ufl.Expression
The function or expression to project.
V : dolfinx.fem.FunctionSpace
The target function space to project onto.
Returns:
u_proj : dolfinx.fem.Function
The projected function in the function space `V`.
"""
# Define trial and test functions in the target function space
u = ufl.TrialFunction(V)
w = ufl.TestFunction(V)
# Define the bilinear form (left-hand side)
a = ufl.inner(u, w) * dx
# Define the linear form (right-hand side) using the function to project
L = ufl.inner(v, w) * dx
# Convert the forms to dolfinx.fem.Form objects
a_form = dolfinx.fem.form(a)
L_form = dolfinx.fem.form(L)
# Assemble the system matrix A and the right-hand side vector b
A = dolfinx.fem.petsc.assemble_matrix(a_form)
A.assemble()
b = dolfinx.fem.petsc.assemble_vector(L_form)
# Create a function in the target space to hold the projected result
u_proj = dolfinx.fem.Function(V)
# Solve the linear system to compute the projection
solver = PETSc.KSP().create(V.mesh.comm)
solver.setOperators(A)
solver.solve(b, u_proj.vector)
u_proj.x.scatter_forward()
return u_proj
SigOut = fem.functionspace(domain, ("DG", 0, (4, )))
sigmaout=fem.Function(SigOut, name="Stresses")
vtk = io.VTKFile(domain.comm, "results1/elastoplastic.pvd", "w")
Nitermax, tol = 200, 1e-6 # 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, 3))
# we set all functions to zero before entering the loop in case we would like to reexecute this code cell
sig.vector.set(0.0)
sig_old.vector.set(0.0)
p.vector.set(0.0)
u.vector.set(0.0)
n_elas.vector.set(0.0)
beta.vector.set(0.0)
Δε = eps(Du)
sig_, n_elas_, beta_, dp_ = constitutive_update(Δε, sig_old, p)
for i, t in enumerate(load_steps):
loading.value = t * q_lim
# compute the residual norm at the beginning of the load step
tangent_problem.assemble_rhs()
nRes0 = tangent_problem._b.norm()
nRes = nRes0
Du.x.array[:] = 0
niter = 0
while nRes / nRes0 > tol and niter < Nitermax:
# solve for the displacement correction
tangent_problem.assemble_lhs()
tangent_problem.solve_system()
# update the displacement increment with the current correction
Du.vector.axpy(1, du.vector) # Du = Du + 1*du
Du.x.scatter_forward()
# interpolate the new stresses and internal state variables
interpolate_quadrature(sig_, sig)
interpolate_quadrature(n_elas_, n_elas)
interpolate_quadrature(beta_, beta)
# compute the new residual
tangent_problem.assemble_rhs()
nRes = tangent_problem._b.norm()
niter += 1
# Update the displacement with the converged increment
u.vector.axpy(1, Du.vector) # u = u + 1*Du
u.x.scatter_forward()
# Update the previous plastic strain
interpolate_quadrature(dp_, dp)
p.vector.axpy(1, dp.vector)
# Update the previous stress
sig_old.x.array[:] = sig.x.array[:]
p_avg.interpolate(project(p, P0))
vtk.write_function(p_avg,t)
sigmaout.interpolate(project(sig,SigOut))
vtk.write_function(sigmaout,t)
vtk.close()