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

dolfinx
1

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
2

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)

to

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

and see if this helps.

Best,
Halvor

2 Likes
3

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
4

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?

Cheers

5

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)
6

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
7

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

1 Like