Dear Community,
I am currently constructing an elastoplastic model using FEniCSx, referencing the content in Elasto-plastic analysis of a 2D von Mises material. In this example, Dirichlet BCs are applied to the left and bottom edges of the quarter hollow cylinder model, with a time-varying loading applied to the inner wall of the hollow cylinder.
I now wish to modify the Dirichlet BCs such that the lower edge has fixed displacement of 0 in both the x and y directions, while the left edge gets a displacement increasing over time in the positive x direction up to 0.05mm (as a time-dependent Dirichlet BC). Additionally, I wish to remove the load applied to the inner wall of the hollow cylinder. So I changed the following part of the code (the rest remains same):
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
u_bottom = fem.Constant(domain, ScalarType((0.000, 0.000)))
bottom_dofs = fem.locate_dofs_topological(V, gdim - 1, facets.find(1))
bc_bottom = fem.dirichletbc(u_bottom, bottom_dofs, V)
u_left = fem.Constant(domain, ScalarType((0.000)))
left_dofs_x = fem.locate_dofs_topological(V.sub(0), gdim - 1, facets.find(2))
bc_left_x = fem.dirichletbc(u_left, left_dofs_x, V.sub(0))
bcs = [bc_left_x, bc_bottom]
Nitermax, tol = 200, 1e-6 # parameters of the Newton-Raphson procedure
Nincr = 20
load_steps = np.linspace(0, 2, 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)
Tensor_element = element("Lagrange", domain.basix_cell(), 1, shape=(2,))
V_Tensor = fem.functionspace(domain, Tensor_element)
Result_u = fem.Function(V_Tensor)
v_load = 0.05e-3
with io.XDMFFile(domain.comm, foldername + "/u.xdmf", "w") as xdmf_u:
xdmf_u.write_mesh(domain)
for i, t in enumerate(load_steps):
# loading.value = t * q_lim
u_left.value = ScalarType((v_load * t))
bc_left_x = fem.dirichletbc(u_left, left_dofs_x, V.sub(0))
bcs=[bc_left_x, bc_bottom]
# 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[:]
if len(bottom_inside_dof) > 0: # test if proc has dof
results[i + 1, :] = (u.x.array[bottom_inside_dof[0]], t, niter)
# -
Result_u.name = "u"
Result_u.interpolate(u)
xdmf_u.write_function(Result_u, t)
However, I have encountered the following error:
Traceback (most recent call last):
File "/root/shared/Code/2026_02_20/plasticity_original_example.py", line 271, in <module>
while nRes / nRes0 > tol and niter < Nitermax:
ZeroDivisionError: float division by zero
I believe I need to modify the CustomLinearProblem function to adjust the BCs, which is:
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",
},
)
But after searching through the FEniCSx tutorials, I remain unable to find a solution. Could someone kindly provide me some suggestions? If anyone needs, I can also provide the full code. Many thanks in advance!