# Compute reaction forces for nonlinear material problem with dirichlet BCs

Hi!
I want to plot the load-displacement curve for a somewhat complex problem which I load using dirichlet BCs, but can’t make it work. I found this nice tutorial from @bleyerj for computing the internal forces. I’ve used it to create an MWE for a simple beam under tension with dirichlet BCs, with a linear material.
Although the response seems correct, I was wondering if this is the proper way to do it, since both my load and this reaction force method overwrite the same boundary vector `u_right`?
The code for this MWE:

``````import numpy as np
from mpi4py import MPI
import ufl
from dolfinx import fem, mesh, nls, io
import petsc4py.PETSc as PETSc
import dolfinx.fem.petsc
import dolfinx.nls.petsc
from dolfinx import nls

L = 5.0
H = 1.0
Nx = 20
Ny = 5
domain = mesh.create_rectangle(
MPI.COMM_WORLD,
[(0.0, -H / 2), (L, H / 2)],
(Nx, Ny),
diagonal=mesh.DiagonalType.crossed,
)

E = fem.Constant(domain, 1e5)
nu = fem.Constant(domain, 0.3)
mu = E / 2 / (1 + nu)
lmbda = E * nu / (1 + nu) / (1 - 2 * nu)

def eps(v):

def sigma(v):
return lmbda * ufl.tr(eps(v)) * ufl.Identity(2) + 2.0 * mu * eps(v)

V = fem.functionspace(domain, ("P", 1, (2,)))
u = fem.Function(V, name="Displacement")
u_ = ufl.TestFunction(V)
residual = ufl.inner(sigma(u), eps(u_)) * ufl.dx

# Dirichlet BC
def left(x):
return np.isclose(x[0], 0.0)

def right(x):
return np.isclose(x[0], L)

fdim = domain.topology.dim - 1
left_facets = mesh.locate_entities_boundary(domain, fdim, left)
right_facets = mesh.locate_entities_boundary(domain, fdim, right)

left_dofs = fem.locate_dofs_geometrical(V, left)
u_bc = fem.Function(V)
u_bc.x.array[:] = 0.0

# Right moving boundary
right_dofs = fem.locate_dofs_topological(V.sub(0), fdim, right_facets)
u_right = fem.Constant(domain, PETSc.ScalarType(1.))

bcs = [fem.dirichletbc(u_bc, left_dofs), fem.dirichletbc(u_right, right_dofs, V.sub(0))]

# Define problem
problem = fem.petsc.NonlinearProblem(residual, u, bcs)

# Newton solver
solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)

solver.rtol = 1e-6
solver.atol = 1e-7  # I use an additional absolute tolerance to avoid convergence issues when having loadsteps with no change in loading: the GP numerical inaccuracies keep giving slight tangent changes.
solver.convergence_criterion = "incremental"
solver.report = True

# Prepare internal forces
reac = fem.Function(V)
virtual_work_form = fem.form(ufl.action(residual, reac))

U_x = np.linspace(0, 0.2 , 11)[1:]

# Plot
file_results = io.XDMFFile(
domain.comm,
f"out_MWE.xdmf",
"w",
)
file_results.write_mesh(domain)

for i in range(len(U_x)):
u_right.value = U_x[i]
print(U_x[i])

its, converged = solver.solve(u)
assert converged
u.x.scatter_forward()

file_results.write_function(u, i)

# Compute internal forces
u_right.value = 1.0
fem.set_bc(reac.vector, bcs)
print(f"Horizontal reaction Rx = {fem.assemble_scalar(virtual_work_form):.6f}")

``````

However, when I apply it to my more complex problem (nonlinear material), the code runs but the reaction forces stay 0. What is the problem here?
I’ve provided the code below, but since it won’t run unless I provide a bunch more code. I hope that seeing the relevant parts of the code is sufficient to identify what my mistake is.

``````import numpy as np
from mpi4py import MPI
import ufl
from dolfinx import io, fem, mesh
from dolfinx.cpp.nls.petsc import NewtonSolver
from df_mat_utils.solvers import NonlinearMaterialProblem
from petsc4py.PETSc import ScalarType
import petsc4py.PETSc as PETSc

from constitutive_models import CustomHyperElastic_jax

"""
Mesh creation
"""
params = {
"mesh_element_size": 0.05,
"E": 3.13e3,
"nu": 0.37,
...
}

filename = f"dogbone_mesh.xdmf"

domain.topology.create_connectivity(domain.topology.dim - 1, domain.topology.dim)

material = CustomHyperElastic_jax(mesh=domain, params=params)

order = 1
deg_quad = 2 * (order - 1)
shape = (2,)

V = fem.functionspace(domain, ("P", order, shape))

def strain(u):
return ufl.as_vector(
[
u[0].dx(0),
u[1].dx(1),
1 / np.sqrt(2) * (u[1].dx(0) + u[0].dx(1)),
]
)

du = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
u = fem.Function(V, name = "u")     # total displacement

"""
Boundary conditions
"""
def bot_left_corner(x):
return np.logical_and(np.isclose(x[1], 0.0), np.isclose(x[0], 0.0))

def bot_right_corner(x):
return np.logical_and(np.isclose(x[1], 0.0), np.isclose(x[0], 1.0))

def left(x):
return np.isclose(x[0], 0)

def right(x):
return np.isclose(x[0], 1)

# Doing it like this, rather than locating geometrically directly, is necessary afaik to get only 1 dof per side.
fdim = domain.topology.dim-1
left_facets = mesh.locate_entities_boundary(domain, fdim, left)
right_facets = mesh.locate_entities_boundary(domain, fdim, right)
left_dofs = fem.locate_dofs_topological(V.sub(0), fdim, left_facets)
right_dofs = fem.locate_dofs_topological(V.sub(0), fdim, right_facets)

# bottom left:
bottom_value1 = fem.Function(V)
bottom_value1.x.array[:] = 0
left_corner_vertex = mesh.locate_entities_boundary(domain, 0, bot_left_corner)
dofs = fem.locate_dofs_topological((V, V), 0, left_corner_vertex)

# bottom right:
ME0y, _ = V.sub(1).collapse()
bottom_value2 = fem.Function(ME0y)
bottom_value2.x.array[:] = 0
right_corner_vertex = mesh.locate_entities_boundary(domain, 0, bot_right_corner)
right_dof = fem.locate_dofs_topological((V.sub(1), ME0y), 0, right_corner_vertex)

u_right = fem.Constant(domain,PETSc.ScalarType(1.))
#               bottom left corner                             bottom right corner                          left x                                                  right x
bcs = [fem.dirichletbc(bottom_value1, dofs, V), fem.dirichletbc(bottom_value2, right_dof, V.sub(1)), fem.dirichletbc(0.0, left_dofs, V.sub(0)), fem.dirichletbc(u_right, right_dofs, V.sub(0))]

"""
Problem setup
"""

sig = qmap.fluxes["Stress"]
Res = ufl.dot(sig, strain(v)) * qmap.dx
Jac = qmap.derivative(Res, u, du)

problem = NonlinearMaterialProblem(qmap, Res, Jac, u, bcs)

newton = NewtonSolver(MPI.COMM_WORLD)
newton.rtol = 1e-6
newton.atol = 1e-7  # I use an additional absolute tolerance to avoid convergence issues when having loadsteps with no change in loading: the GP numerical inaccuracies keep giving slight tangent changes.
newton.convergence_criterion = "incremental"
newton.report = True

U_x = np.linspace(0, 0.2 , 11)[1:]

# Prep to obtain the load value
# Based on: https://bleyerj.github.io/comet-fenicsx/tips/computing_reactions/computing_reactions.html#
reac = fem.Function(V)
virtual_work_form = fem.form(ufl.action(Res, reac))

for i in range(len(U_x)):
u_right.value = U_x[i]

converged, nr_iters = problem.solve(newton)

# Get internal forces, to plot load-displacement curve
u_right.value = 1.0
fem.set_bc(reac.vector, bcs)
print(f"Horizontal reaction Rx = {fem.assemble_scalar(virtual_work_form):.6f}")

``````

Thank you!

I’ve solved this issue by foregoing the method completely, and instead using a custom problem setup. I’ve modified the CustomNewtonProblem by BleyerJ, here is the full code:

``````class CustomNewtonProblem:
def __init__(self, quadrature_map, F, J, u, bcs, max_it=50, rtol=1e-8, atol=1e-8):
if isinstance(F, list):
self.L = [form(f) for f in F]
else:
self.L = form(F)
if isinstance(J, list):
self.a = [form(j) for j in J]
else:
self.a = form(J)
self.bcs = bcs
self._F, self._J = None, None
self.u = u
self.du = Function(self.u.function_space, name="Increment")
self.it = 0

if isinstance(self.a, list):
self.A = create_matrix(self.a[0])
else:
self.A = create_matrix(self.a)  # preallocating sparse jacobian matrix
if isinstance(self.L, list):
self.b = create_vector(self.L[0])  # preallocating residual vector
else:
self.b = create_vector(self.L)  # preallocating residual vector
self.max_it = max_it
self.rtol = rtol
self.atol = atol

self.internal_forces = None

def solve(self, solver, print_steps=True, print_solution=True):
i = 0  # number of iterations of the Newton solver
converged = False
while i < self.max_it:
with Timer("Constitutive update"):

# Assemble Jacobian and residual
with self.b.localForm() as loc_b:
loc_b.set(0)
self.A.zeroEntries()
if isinstance(self.a, list):
self.A.assemble()
for ai in self.a:
Ai = assemble_matrix(ai, bcs=self.bcs)
Ai.assemble()
self.A.axpy(1.0, Ai)
else:
assemble_matrix(self.A, self.a, bcs=self.bcs)
self.A.assemble()

if isinstance(self.L, list):
for Li in self.L:
bi = assemble_vector(Li)
self.b.axpy(1.0, bi)
else:
assemble_vector(self.b, self.L)
self.b.ghostUpdate(
)
self.b.scale(-1)

# Store internal_forces before applying BCs
self.internal_forces = self.b.copy()

# Compute b - J(u_D-u_(i-1))
if isinstance(self.a, list):
for ai in self.a:
apply_lifting(self.b, [ai], [self.bcs], x0=[self.u.vector], scale=1)
else:
apply_lifting(self.b, [self.a], [self.bcs], x0=[self.u.vector], scale=1)
# Set dx|_bc = u_{i-1}-u_D
set_bc(self.b, self.bcs, self.u.vector, 1.0)
self.b.ghostUpdate(
)

# Solve linear problem
solver.setOperators(self.A)
with Timer("Linear solve"):
solver.solve(self.b, self.du.vector)
self.du.x.scatter_forward()

# Update u_{i+1} = u_i + relaxation_param * delta x_i
self.u.x.array[:] += self.du.x.array[:]
i += 1
# Compute norm of update
correction_norm = self.du.vector.norm(0)
error_norm = self.b.norm(0)
if i == 1:
error_norm0 = error_norm
relative_norm = error_norm / error_norm0
if print_steps:
print(
f"    Iteration {i}:  Residual abs: {error_norm} rel: {relative_norm}"
)

if relative_norm < self.rtol or error_norm < self.atol:
converged = True
break

if print_solution:
if converged:
print(f"Solution reached in {i} iterations.")
else:
print(
f"No solution found after {i} iterations. Revert to previous solution and adjust solver parameters."
)

return converged, i

def get_reaction_forces(self, boundary_dofs):
# Return the *negative* of the internal forces at the specified dofs
if self.internal_forces is None:
raise ValueError("Solve the problem first to compute internal forces.")
return - self.internal_forces.getValues(boundary_dofs)
``````
1 Like