Mixed elements boundary conditions for switching from plane strain to plane stress

I’ve used an adapted version of this demo, which works. However, when using plane stress instead of plane strain, the NewtonSolver no longer finds a solution.

For plane stress I switch to mixed elements, and I believe the problem is likely related to the definition of the boundary conditions (I’ve seen this fix which is already used here).

Since I’m using an adapted version of dolfinx_materials, the underlying simplified code won’t work for you, but since the plane_strain version works and the differences between plane strain and plane stress are minimal (search for “hypothesis” below) and only use dolfinx, I hope it is sufficient to find the problem.

import numpy as np
from mpi4py import MPI
import ufl
from dolfinx import mesh, fem
from dolfinx.cpp.nls.petsc import NewtonSolver
from dolfinx_materials.quadrature_map import QuadratureMap
from dolfinx_materials.solvers import NonlinearMaterialProblem

from const_model_hybrid import PF_mixture_model
from petsc4py.PETSc import ScalarType

domain = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)

# hypothesis = "plane_strain"
hypothesis = "plane_stress"

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

if hypothesis == "plane_strain":
    V = fem.functionspace(domain, ("P", order, shape))

    def strain(u):
        return ufl.as_vector( [ u[0].dx(0), u[1].dx(1), 0.0, 1 / np.sqrt(2) * (u[1].dx(0) + u[0].dx(1)), 0.0, 0.0, ] )
elif hypothesis == "plane_stress":
    Ue = ufl.VectorElement("CG", domain.ufl_cell(), order)
    Ee = ufl.FiniteElement("DG", domain.ufl_cell(), order - 1)
    V = fem.FunctionSpace(domain, ufl.MixedElement([Ue, Ee]))
    V_u, _ = V.sub(0).collapse()

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

du = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
v_u, v_e = ufl.split(v)
u = fem.Function(V, name = "u")

# Boundary conditions
ds_m = ufl.Measure("ds", domain=domain)
n = ufl.FacetNormal(domain)

# Dirichlet
def top(x):
    return np.isclose(x[1], 1)

if hypothesis == "plane_strain":
    dofs = fem.locate_dofs_geometrical(V, top)
    bc_top = fem.dirichletbc(fem.Function(V), dofs)
elif hypothesis == "plane_stress":
    dofs = fem.locate_dofs_geometrical((V.sub(0), V_u), top)
    bc_top = fem.dirichletbc(fem.Function(V_u), dofs)
bcs = [bc_top]

# Neumann
q_lim = 40
loading = fem.Constant(domain, ScalarType(0. * q_lim))

dataloc = "../data/neo_hookean_100_1.0"
material = PF_mixture_model(domain, dataloc, E=3.13e3, nu=0.37)

qmap = QuadratureMap(domain, deg_quad, material)
qmap.register_gradient(material.gradient_names[0], strain(u))

sig = qmap.fluxes["Stress"]
if hypothesis == "plane_strain":
    Res = ufl.dot(sig, strain(v)) * qmap.dx - loading * ufl.inner(ufl.as_vector([0, -1]), v) * ds_m
elif hypothesis == "plane_stress":
    Res = ufl.dot(sig, strain(v)) * qmap.dx - loading * ufl.inner(ufl.as_vector([0, -1]), v_u) * ds_m
Jac = qmap.derivative(Res, u, du)

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

newton = NewtonSolver(MPI.COMM_WORLD)
newton.rtol = 1e-6
newton.convergence_criterion = "incremental"
newton.report = True

load_steps = np.linspace(0, 1, 20)[1:] ** 0.5

for i, load_step in enumerate(load_steps):
    print(f"Iteration {i}/{len(load_steps)}")
    loading.value = load_step * q_lim
    converged, it = problem.solve(newton)

When using locate_dofs_topological I get:

Process finished with exit code 139 (interrupted by signal 11:SIGSEGV)

and using locate_dofs_geometrical:

RuntimeError: Newton solver did not converge because maximum number of iterations reached

Hi surrogater,

As in lines 96 and 97 of the demo plane_elastoplasticity.py you linked to, when using mixed elements, you need to specify which subspace of the mixed finite element space you are setting the Dirichlet BCs in.

Try changing

bc_top = fem.dirichletbc(fem.Function(V_u), dofs)


bc_top = fem.dirichletbc(fem.Function(V_u), dofs, V.sub(0))

and see if this helps.



Thank you for the reponse hherlyng.
That was a mistake I made while creating this more minimal example, thank you for pointing it out!
It does not resolve my problem unfortunately.

For this simplified example, I can find a solution to the first load step in 40 newton raphson steps, and the second load step in 49 NR steps. Plane strain takes 4 and 3 newton raphson steps for these respectively.

I was wondering if my use of v and v_u here looks correct?

if hypothesis == "plane_strain":
    Res = ufl.dot(sig, strain(v)) * qmap.dx - loading * ufl.inner(ufl.as_vector([0, -1]), v) * ds_m
elif hypothesis == "plane_stress":
    Res = ufl.dot(sig, strain(v)) * qmap.dx - loading * ufl.inner(ufl.as_vector([0, -1]), v_u) * ds_m

At a first glance that looks correct, yes, since you have extracted v_u from the mixed trial function.

When using locate_dofs_topological you said you get an error with signal 11:SIGSEGV, this looks like a segmentation fault to me. Do you get the error on the exact line where you used locate_dofs_topological, or at some later instance?


I get the SIGSEGV error on this line of code when calling the petsc NewtonSolver’s solve method as:

it, converged = solver.solve(self.u.vector)

I think I’ve found the problem, which is a wrong tangent computation in my constitutive model. This then explains why the convergence was either super slow or not happening.

My apologies for troubling you with this.

1 Like

No worries, great to hear that you’ve solved your problem!

1 Like