# 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.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

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

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

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

to

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

and see if this helps.

Best,
Halvor

2 Likes

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?

Cheers

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