`dolfinx_mpc.NonlinearProblem` with periodic geometrical constraints diverges for mixed hyperelasticity

Hi Jorgen, many thanks for the fix in https://github.com/jorgensd/dolfinx_mpc/pull/237. I was trying to implement a homogenization example which fails on the nonlinear solve and wanted to check what could be going wrong in the problem setup (since the same material model with normal dirichlet BC's converges just fine).

## Summary

`dolfinx_mpc.NonlinearProblem` with `create_periodic_constraint_geometrical` on a
**mixed displacement–pressure** space diverges with `SNES reason = -6`
(`DIVERGED_LINE_SEARCH`) at every load step. The identical problem with standard
Dirichlet BCs converges in 2–3 Newton iterations per step.

---

## Problem description

### Kinematics

Unit cube $\Omega = [0,1]^3$, periodic homogenization. Total deformation gradient:

$$\mathbf{F} = \mathbf{I} + \overline{\mathbf{H}} + \nabla\tilde{\mathbf{u}},$$

where $\overline{\mathbf{H}}$ is a prescribed macro gradient (a `fem.Constant`) and
$\tilde{\mathbf{u}}$ is the periodic fluctuation.

### Constitutive model

Near-incompressible Neo-Hookean,
mixed space $[\mathbb{P}_2]^3 \times \mathbb{P}_1$:

$$\psi = p\,(J - J^\star) + \tfrac{\mu}{2}(I_1 - 3) - \mu\ln J^\star + \tfrac{\kappa}{2}(J^\star-1)^2,$$

$$I_1 = \text{tr}(\mathbf{F}^T\mathbf{F}), \quad J = \det\mathbf{F}, \quad
J^\star = \frac{(\kappa+p) + \sqrt{(\kappa+p)^2 + 4\kappa\mu}}{2\kappa}.$$

Parameters: $\mu = 1$, $\kappa = 200$.

### Periodic constraint [Abaqus `*EQUATION` style](http://130.149.89.49:2080/v2016/books/usb/default.htm?startat=pt08ch35s02aus131.html#usb-cni-pequation)

For every opposite-face node pair $(\mathbf{x}^+, \mathbf{x}^-)$ with
$\mathbf{x}^+ = \mathbf{x}^- + \mathbf{e}_k$:

$$\tilde{u}_i(\mathbf{x}^+) - \tilde{u}_i(\mathbf{x}^-) = 0.$$

In Abaqus CTENODS notation (with corner reference nodes $C_0, C_k$):

$$u_i^{\text{slave}} - u_i^{\text{master}} - u_i^{C_k} + u_i^{C_0} = 0,$$

```
*EQUATION
4
slave_node, i, 1.0
master_node, i, -1.0
Ck_node, i, -1.0
C0_node, i, 1.0
```

Here $\overline{\mathbf{H}}$ enters via `fem.Constant` directly in UFL, so no
macroscopic strain is carried by the corner nodes; the constraint reduces to
slave = master (zero RHS). The origin is pinned for rigid-body modes; pressure
is pinned at the origin.

```
Periodic face-pairing (axis 0):

y y
^ slave (x=1) ^ master (x=0)
| ●──●──● | ●──●──●
| │ │ │ ← ũᵢˢ = ũᵢᵐ → | │ │ │
| ●──●──● | ●──●──●
└─────────→ z └─────────→ z
```

### `dolfinx_mpc` implementation

```python
for comp in range(3):
for axis in range(3):
mpc.create_periodic_constraint_geometrical(
W.sub(0).sub(comp), # displacement component
_periodic_indicator(axis), # slave: high face of axis
_periodic_relation(axis), # mapping: x[axis] -= 1
[bc_origin],
1.0,
)
```

---

## Reproducer

issue_mpc_periodic_nonlinear.py

```python
from __future__ import annotations
import importlib.metadata
import basix.ufl
import dolfinx_mpc
import numpy as np
import ufl
from dolfinx import default_scalar_type, fem, mesh
from dolfinx.fem import petsc as fem_petsc
from mpi4py import MPI
from petsc4py import PETSc

EXX = 0.03 # applied macro strain
NCELLS = 3 # uniform NCELLS³ mesh
NSTEPS = 6 # load steps
MU = 1.0 # shear modulus
BULK = 200.0 # bulk modulus (bulk/mu = 200 → near-incompressible)
TOL = 1e-10

_SNES_OPTS = {
"snes_type": "newtonls",
"snes_linesearch_type": "bt",
"snes_max_it": 30,
"snes_rtol": 1e-8,
"snes_atol": 1e-8,
"ksp_type": "preonly",
"pc_type": "lu",
}

def _face(axis, value):
return lambda x: np.isclose(x[axis], value, atol=TOL)

def _origin(x):
return np.isclose(x[0], 0.) & np.isclose(x[1], 0.) & np.isclose(x[2], 0.)

def _periodic_indicator(axis):
def f(x):
mask = np.isclose(x[axis], 1., atol=TOL)
for later in range(axis + 1, 3):
mask &= x[later] < 1. - TOL
return mask
return f

def _periodic_relation(axis):
def f(x):
y = x.copy(); y[axis] -= 1.; return y
return f

def _zero_fn(W_sub):
V, _ = W_sub.collapse()
fn = fem.Function(V); fn.x.array[:] = 0.; fn.x.scatter_forward()
return V, fn

def _build_forms(msh, W, up, macro_H):
mu = fem.Constant(msh, default_scalar_type(MU))
kappa = fem.Constant(msh, default_scalar_type(BULK))
du, v = ufl.TrialFunction(W), ufl.TestFunction(W)
u, p = ufl.split(up)
I = ufl.Identity(3)
F = ufl.variable(I + macro_H + ufl.grad(u))
I1 = ufl.tr(F.T * F)
J = ufl.det(F)
J_star = ((kappa + p) + ufl.sqrt((kappa + p)**2 + 4.*kappa*mu)) / (2.*kappa)
psi = (p*(J - J_star) + .5*mu*(I1 - 3.)
- mu*ufl.ln(J_star) + .5*kappa*(J_star - 1.)**2)
R = ufl.derivative(psi * ufl.dx, up, v)
K = ufl.derivative(R, up, du)
return R, K

# ── Case 1: Dirichlet BCs (converges) ────────────────────────────────────────
def run_nonperiodic():
comm = MPI.COMM_WORLD
msh = mesh.create_unit_cube(comm, NCELLS, NCELLS, NCELLS)
u_el = basix.ufl.element(basix.ElementFamily.P, msh.basix_cell(), 2,
shape=(3,), dtype=default_scalar_type)
p_el = basix.ufl.element(basix.ElementFamily.P, msh.basix_cell(), 1,
dtype=default_scalar_type)
W = fem.functionspace(msh, basix.ufl.mixed_element([u_el, p_el]))
V_ux, ux0 = _zero_fn(W.sub(0).sub(0))
V_uy, uy0 = _zero_fn(W.sub(0).sub(1))
V_uz, uz0 = _zero_fn(W.sub(0).sub(2))
Q_p, p0 = _zero_fn(W.sub(1))
ux_load = fem.Function(V_ux)
bcs = [
fem.dirichletbc(ux0, fem.locate_dofs_geometrical((W.sub(0).sub(0), V_ux), _face(0, 0.)), W.sub(0).sub(0)),
fem.dirichletbc(uy0, fem.locate_dofs_geometrical((W.sub(0).sub(1), V_uy), _face(1, 0.)), W.sub(0).sub(1)),
fem.dirichletbc(uz0, fem.locate_dofs_geometrical((W.sub(0).sub(2), V_uz), _face(2, 0.)), W.sub(0).sub(2)),
fem.dirichletbc(ux_load, fem.locate_dofs_geometrical((W.sub(0).sub(0), V_ux), _face(0, 1.)), W.sub(0).sub(0)),
fem.dirichletbc(p0, fem.locate_dofs_geometrical((W.sub(1), Q_p), _origin), W.sub(1)),
]
up = fem.Function(W)
macro_H = fem.Constant(msh, default_scalar_type(np.zeros((3, 3))))
R, K = _build_forms(msh, W, up, macro_H)
problem = fem_petsc.NonlinearProblem(R, up, bcs=bcs, J=K,
petsc_options_prefix="nonperiodic_",
petsc_options=_SNES_OPTS)
ok = True
for step, alpha in enumerate(np.linspace(0., 1., NSTEPS + 1)):
ux_load.x.array[:] = alpha * EXX; ux_load.x.scatter_forward()
if step == 0:
continue
try:
problem.solve()
reason = problem.solver.getConvergedReason()
iters = problem.solver.getIterationNumber()
print(f"[nonperiodic] step={step} alpha={alpha:.4f} reason={reason} iters={iters}")
except PETSc.Error as exc:
print(f"[nonperiodic] step={step} alpha={alpha:.4f} FAILED: {exc}")
ok = False; break
return ok

# ── Case 2: dolfinx_mpc periodic BCs (diverges — this is the bug) ────────────
def run_periodic():
comm = MPI.COMM_WORLD
msh = mesh.create_unit_cube(comm, NCELLS, NCELLS, NCELLS)
u_el = basix.ufl.element(basix.ElementFamily.P, msh.basix_cell(), 2,
shape=(3,), dtype=default_scalar_type)
p_el = basix.ufl.element(basix.ElementFamily.P, msh.basix_cell(), 1,
dtype=default_scalar_type)
W = fem.functionspace(msh, basix.ufl.mixed_element([u_el, p_el]))
mpc = dolfinx_mpc.MultiPointConstraint(W)
bcs = []
for comp in range(3):
V_ui, ui0 = _zero_fn(W.sub(0).sub(comp))
dofs = fem.locate_dofs_geometrical((W.sub(0).sub(comp), V_ui), _origin)
bc = fem.dirichletbc(ui0, dofs, W.sub(0).sub(comp))
bcs.append(bc)
for axis in range(3):
mpc.create_periodic_constraint_geometrical(
W.sub(0).sub(comp),
_periodic_indicator(axis),
_periodic_relation(axis),
[bc], default_scalar_type(1.), default_scalar_type(TOL),
)
Q_p, p0 = _zero_fn(W.sub(1))
bcs.append(fem.dirichletbc(p0,
fem.locate_dofs_geometrical((W.sub(1), Q_p), _origin),
W.sub(1)))
mpc.finalize()
up = fem.Function(mpc.function_space)
macro_H = fem.Constant(msh, default_scalar_type(np.zeros((3, 3))))
R, K = _build_forms(msh, W, up, macro_H)
problem = dolfinx_mpc.NonlinearProblem(
R, up, mpc=mpc, bcs=bcs, J=K,
petsc_options={**_SNES_OPTS, "snes_error_if_not_converged": False},
)
ok = True
for step, alpha in enumerate(np.linspace(0., 1., NSTEPS + 1)):
macro_H.value[...] = default_scalar_type(np.diag([alpha * EXX, 0., 0.]))
if step == 0:
continue
_, reason, iters = problem.solve()
converged = reason > 0
print(f"[periodic] step={step} alpha={alpha:.4f} reason={reason} iters={iters} converged={converged}")
if not converged:
ok = False
return ok

def _print_versions():
import dolfinx, basix, ufl, petsc4py, mpi4py, numpy as np
libs = [
("dolfinx", getattr(dolfinx, "__version__", "n/a")),
("basix", getattr(basix, "__version__", "n/a")),
("ufl", getattr(ufl, "__version__", "n/a")),
("petsc4py", getattr(petsc4py, "__version__", "n/a")),
("mpi4py", getattr(mpi4py, "__version__", "n/a")),
("numpy", np.__version__),
]
for name in ("dolfinx_mpc",):
try:
libs.append((name, importlib.metadata.version(name)))
except Exception:
import dolfinx_mpc as _m
libs.append((name, getattr(_m, "__version__", "n/a")))
print("Library versions:")
for name, ver in libs:
print(f" {name:<14} {ver}")
print()

if __name__ == "__main__":
_print_versions()
print("=" * 60)
print("Case 1: non-periodic (Dirichlet BCs)")
print("=" * 60)
ok_np = run_nonperiodic()
print("=" * 60)
print("Case 2: periodic (dolfinx_mpc geometrical constraints)")
print("=" * 60)
ok_p = run_periodic()
print("=" * 60)
print(f"nonperiodic converged : {ok_np}")
print(f"periodic converged : {ok_p}")
```

---

## Output

```
Library versions:
dolfinx 0.11.0.dev0
basix 0.11.0.dev0
ufl 2025.3.0.dev0
petsc4py 3.24.5
mpi4py 4.1.1
numpy 2.4.3
dolfinx_mpc 0.10.0.dev0

============================================================
Case 1: non-periodic (Dirichlet BCs)
============================================================
[nonperiodic] step=1 alpha=0.1667 reason=2 iters=3
[nonperiodic] step=2 alpha=0.3333 reason=2 iters=2
[nonperiodic] step=3 alpha=0.5000 reason=2 iters=2
[nonperiodic] step=4 alpha=0.6667 reason=2 iters=2
[nonperiodic] step=5 alpha=0.8333 reason=2 iters=2
[nonperiodic] step=6 alpha=1.0000 reason=2 iters=2
============================================================
Case 2: periodic (dolfinx_mpc geometrical constraints)
============================================================
[periodic] step=1 alpha=0.1667 reason=-6 iters=1 converged=False
[periodic] step=2 alpha=0.3333 reason=-6 iters=0 converged=False
[periodic] step=3 alpha=0.5000 reason=-6 iters=0 converged=False
[periodic] step=4 alpha=0.6667 reason=-6 iters=0 converged=False
[periodic] step=5 alpha=0.8333 reason=-6 iters=0 converged=False
[periodic] step=6 alpha=1.0000 reason=-6 iters=0 converged=False
============================================================
nonperiodic converged : True
periodic converged : False
```

`reason = -6` → `SNES_DIVERGED_LINE_SEARCH` at every load step.

---

## Expected behaviour

The periodic and non-periodic formulations are mathematically equivalent for a
homogeneous single-material cube (no inclusions). Both should converge in a small
number of Newton iterations.

## Actual behaviour

`dolfinx_mpc.NonlinearProblem` diverges immediately with
`SNES_DIVERGED_LINE_SEARCH` (`reason = -6`) at every load step,
even at very small applied strains ($\overline{H}_{11} \approx 0.005$).

Any clues what is wrong with the setup ? Maybe the way the mpc is setup?