Shape derivatives

[jorgensd]

Adds shape control: differentiate a functional with respect to the geometry of the mesh it is posed on.

S = dolfinx_adjoint.geometry_function_space(mesh)
s = dolfinx_adjoint.Function(S)
dolfinx_adjoint.move(mesh, s) # annotating counterpart of scifem.mesh.move
uh = dolfinx_adjoint.LinearProblem(a, L, bcs=bcs).solve()
Jhat = pyadjoint.ReducedFunctional(dolfinx_adjoint.assemble_scalar(j(uh)), pyadjoint.Control(s))
The control is the displacement, never the mesh — the mesh is the intermediate block variable linking it to every form, mirroring dolfin-adjoint's ALE.move. Geometry dependence is carried by ufl.SpatialCoordinate, and ufl.derivative with respect to it is the shape derivative.

Covered, each verified by value against a central difference of an independently rebuilt forward: assemble_scalar over dx/ds/dS and error_norm; LinearProblem and NonlinearProblem; blocked (Taylor-Hood) problems; time-dependent problems; higher-order (curved P2) meshes; multiple moves; interpolate of a coordinate expression and Dirichlet values built from one; and the shape Hessian of a functional with no PDE in between.

The adjoint doesn't use dolfin-adjoint's workaround
For a rank-1 form the symbolic route this codebase uses everywhere else — ufl.action(ufl.adjoint(dFdX), λ) — raises: UFL cannot expand a coordinate derivative of a coefficient in physical space. dolfin-adjoint assembles dF/dX as a rectangular matrix and transposes it with PETSc.

This contracts before differentiating instead. F is linear in its test function, so substituting the adjoint solution for it turns F into a functional whose coordinate derivative is a one-form on the geometry space:

d/dX_j F(u, λ) = Σ_i λ_i ∂F_i/∂X_j = [(∂F/∂X)*λ]_j

Verified against the matrix route to 1e-15. No PETSc Mat is created (so it can't reintroduce the collective-destruction hazard), and it extends to a mixed test space for free.

Bugs fixed
Every unblocked vector-valued problem failed at derivative(). compute_adjoint called ufl.extract_blocks unconditionally, which splits any vector element into scalar components DOLFINx can't compile. It survived because every vector-valued test problem in the suite is blocked, where the split is correct. Fixed using the blocked: bool parameter already written on dokken/topology-optimization — expect a conflict when that branch merges.
interpolate of a coordinate expression dropped its geometry term — 5.3% wrong silently; a Dirichlet value built from one, worse. Now 1.7e-11 and 2.7e-08.
Guarded rather than fixed — all raise, none silent
A recomputing checkpoint schedule (refused at move(), before the forward runs); TLM/Hessian across a solve; an expression mixing a coefficient with SpatialCoordinate; a form containing CellDiameter/Circumradius/MinCellEdgeLength; a displacement outside the geometry space.

Each would otherwise return a plausible, self-consistent gradient that a Taylor test reports as converging at rate 2 — the recomputing-schedule case was 17% wrong that way. Two guards are allowlists, deliberately: an unrecognised schedule or geometric quantity is refused rather than waved through.

Demos
shape_optimization.py — torsional rigidity, a square rounds towards a disk (+11.1%). stokes_shape_optimization.py — the dolfin-adjoint Stokes drag-minimization demo with its mesh generation folded in and a second-order mesh: Pironneau's rugby ball at aspect ratio 2.59, dissipation −11.7%, obstacle volume held to 0.19%.

Known gap
NonmatchingInterpolationBlock registers no mesh dependency and point location between meshes is geometry-dependent. Untested; presumably wrong in the same silent way. The one case that neither works nor raises.

18 files, +2582/−37. 157 passed, 1 xfailed, 2 xpassed serially; ruff and mypy clean. The two xpasses are pre-existing version-gated pyadjoint markers.

Two caveats on that last line, since you'll be the one merging: the mpirun -n 2 run against the final revision hadn't finished when I wrote this — earlier revisions all passed with the same counts, but it's worth confirming. And the xpass claim rests on boundary-control/branch-summary.md recording the same pair, not on my re-running main.