Compute reaction forces for nonlinear material problem with dirichlet BCs

dolfinx
1

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):
    return ufl.sym(ufl.grad(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))

# Load
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.quadrature_map import QuadratureMap
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"

mesh_reader = io.XDMFFile(MPI.COMM_WORLD,filename,"r")

domain = mesh_reader.read_mesh(name="dogbone")
domain.topology.create_connectivity(domain.topology.dim - 1, domain.topology.dim)
print(f"Finished loading mesh...")

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
"""
qmap = QuadratureMap(domain, deg_quad, material)
qmap.register_gradient(material.gradient_names[0], strain(u))

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!

2

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):
        self.quadrature_map = quadrature_map
        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"):
                self.quadrature_map.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(
                addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE
            )
            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(
                addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD
            )

            # 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.")
                self.quadrature_map.advance()
            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