Thickness optimization and shape derivatives

Two things the analysis walkthrough only mentions: driving an optimizer with Hermit’s outputs, and differentiating with respect to the mesh coordinates.

Both run here for real, on a small cantilever plate. The optimization needs the optional solver stack (pip install 'hermit[opt]').

Setup

The same 2 × 10 cantilever plate as the analysis walkthrough, coarser so the optimizer finishes quickly.

import inspect

import numpy as np
import basix.ufl
import ufl
from mpi4py import MPI
from dolfinx.mesh import create_mesh
import csdl_alpha as csdl

import hermit as hm

LENGTH, WIDTH = 10.0, 2.0
E_VAL, NU_VAL, H_VAL, RHO_VAL = 4.32e8, 0.0, 0.2, 1.0
PRESSURE_Z = 2.0
nx, ny = 10, 2


def plate_mesh(nx, ny):
    xs = np.linspace(0.0, LENGTH, nx + 1)
    ys = np.linspace(0.0, WIDTH, ny + 1)
    points = np.array([[x, y, 0.0] for y in ys for x in xs])
    cells = np.array([[j * (nx + 1) + i, j * (nx + 1) + i + 1,
                       (j + 1) * (nx + 1) + i, (j + 1) * (nx + 1) + i + 1]
                      for j in range(ny) for i in range(nx)], dtype=np.int64)
    el = ufl.Mesh(basix.ufl.element("Lagrange", "quadrilateral", 1, shape=(3,)))
    if list(inspect.signature(create_mesh).parameters)[2] == "e":     # DOLFINx 0.11
        return create_mesh(MPI.COMM_WORLD, cells, el, points)
    return create_mesh(MPI.COMM_WORLD, cells, points, el)             # DOLFINx 0.9


mesh = plate_mesh(nx, ny)
print(f"{mesh.geometry.x.shape[0]} vertices")

33 vertices

Thickness optimization

Minimise compliance at fixed mass. The whole model is one CSDL graph, so the optimizer interface is three calls on the variables — set_as_design_variable, set_as_constraint, set_as_objective — and then handing the recorder to modopt.

Choose the design space deliberately. Giving every node its own free thickness is the obvious thing to do and it is ill-posed: bending stiffness goes as \(\int t^3\) while mass goes as \(\int t\), so at fixed mass an oscillating thickness genuinely raises the stiffness. A free nodal design space therefore has mesh-dependent checkerboard optima that really do beat any smooth profile in the discrete objective — the optimizer is right, the target is wrong. The usual remedies are a filter, a perimeter penalty, or a restricted design space.

We take the last, and it makes the example sharper. Restrict the thickness to \(t(x) = A\,\xi^{p}\) with \(\xi = (L-x)/L + \varepsilon\), and let the optimizer find both \(A\) and the exponent \(p\). In the beam limit \(C = \int M^2/t^3\,\mathrm{d}x\), and minimising that at fixed \(\int t\,\mathrm{d}x\) gives \(t \propto \sqrt{M}\); under uniform load \(M \propto (L-x)^2\), so the optimum is the linear wedge \(p = 1\). That is a real answer to recover, so we start \(p\) deliberately away from it.

recorder = csdl.Recorder(inline=True)
recorder.start()

domain = hm.ShellDomain(mesh, element="CG2CG1")
bcs = hm.clamp(domain, where=hm.near("x", 0.0))
loads = hm.pressure(domain, PRESSURE_Z)

EPS = 1e-3                                        # keeps the tip thickness positive
xi = (LENGTH - domain.node_coords[:, 0]) / LENGTH + EPS

# t(x) = A * xi**p. A variable exponent needs exp(p log xi); xi is a constant array,
# so this stays a clean CSDL expression in the two design variables.
scale = csdl.Variable(value=np.array([H_VAL]), name="scale")
exponent = csdl.Variable(value=np.array([0.3]), name="exponent")   # away from p = 1
log_xi = csdl.Variable(value=np.log(xi))
thickness = csdl.expand(scale, xi.shape) * csdl.exp(csdl.expand(exponent, xi.shape) * log_xi)

material = hm.isotropic(domain, E=E_VAL, nu=NU_VAL,
                        thickness=hm.from_nodal(domain, thickness), density=RHO_VAL)
state = hm.solve(domain, material, loads, bcs)
compliance, mass = hm.compliance(state), hm.mass(state)

c0 = float(np.ravel(compliance.value)[0])
mass_target = RHO_VAL * H_VAL * LENGTH * WIDTH
print(f"starting point: p = 0.3, compliance {c0:.6e}, mass {float(np.ravel(mass.value)[0]):.4f}")

starting point: p = 0.3, compliance 1.689510e-01, mass 3.0465

scale.set_as_design_variable(lower=1e-3, upper=1.0)
exponent.set_as_design_variable(lower=0.0, upper=3.0)
mass.set_as_constraint(lower=mass_target, upper=mass_target)
compliance.set_as_objective()

from modopt import CSDLAlphaProblem, PySLSQP

sim = csdl.experimental.PySimulator(recorder)
problem = CSDLAlphaProblem(problem_name="tutorial_taper", simulator=sim)
PySLSQP(problem, solver_options={"maxiter": 200, "acc": 1e-12}).solve()

recorder.stop()

Optimization terminated successfully    (Exit mode 0)
            Final objective value                : 4.205578e-02
            Final optimality                     : 1.892155e-14
            Final feasibility                    : 3.002043e-13
            Number of major iterations           : 8
            Number of function evaluations       : 8
            Number of derivative evaluations     : 8
            Average Function evaluation time     : 0.092107 s per evaluation
            Average Derivative evaluation time   : 0.077418 s per evaluation
            Total Function evaluation time       : 0.736854 s [ 54.26%]
            Total Derivative evaluation time     : 0.619346 s [ 45.61%]
            Optimizer time                       : 0.000200 s [  0.01%]
            Processing time                      : 0.001602 s [  0.12%]
            Visualization time                   : 0.000000 s [  0.00%]
            Total optimization time              : 1.358003 s [100.00%]
            Summary saved to                     : tutorial_taper_outputs/2026-09-09_14.40.22.510003/slsqp_summary.out

p_opt = float(np.ravel(exponent.value)[0])
c1 = float(np.ravel(compliance.value)[0])

print(f"recovered exponent : p = {p_opt:.4f}   (analytic 1.0)")
print(f"compliance         : {c0:.6e} -> {c1:.6e}   ({100 * (1 - c1 / c0):.1f}% lower)")
print(f"mass               : {float(np.ravel(mass.value)[0]):.6f}   (target {mass_target})")

recovered exponent : p = 1.0603   (analytic 1.0)
compliance         : 1.689510e-01 -> 4.205578e-02   (75.1% lower)
mass               : 4.000000   (target 4.0)

The optimizer recovers the analytic exponent to within about 6% on this deliberately tiny mesh, driving through the shell solve’s adjoint from a starting point well away from the answer. The remaining error is honest: the plate has finite width and transverse shear, and the tip thickness is \(\varepsilon\) rather than zero, so it is not exactly the beam the closed form assumes. ex_optimal_thickness_taper.py runs this as a gated benchmark on finer meshes and holds it to 5%.

That benchmark’s docstring also records the measurements behind the ill-posedness warning above — on its fixture, at equal mass, uniform gives \(C = 0.056\), the analytic wedge \(C = 0.017\), and a free-nodal optimum \(C = 0.0037\). The checkerboard beats the wedge by five times, and it is right to. examples/advanced_examples/ex_thickness_opt.py runs the free-nodal problem if you want to see that for yourself.

Shape (mesh-coordinate) derivatives

Hermit differentiates the residual and every output functional with respect to ufl.SpatialCoordinate, then scatters the sensitivity back to your node ordering. This is the gradient a shape optimizer, or a chain rule back to a geometry parameterisation, consumes.

Make the geometry a live input by passing geometry= a hm.geometry(...) built from a csdl.Variable — either node_disp= (a perturbation added to the reference coordinates) or nodes= (absolute coordinates). Leave it out and the geometry is a constant: no mesh-derivative forms are built, so there is no cost to not using this.

recorder = csdl.Recorder(inline=True)
recorder.start()

domain = hm.ShellDomain(plate_mesh(nx, ny), element="CG2CG1")
bcs = hm.clamp(domain, where=hm.near("x", 0.0))
material = hm.isotropic(domain, E=E_VAL, nu=NU_VAL, thickness=H_VAL, density=RHO_VAL)

node_disp = csdl.Variable(value=np.zeros((domain.n_nodes, 3)), name="node_disp")
state = hm.solve(domain, material, hm.pressure(domain, PRESSURE_Z), bcs,
                 geometry=hm.geometry(domain, node_disp=node_disp))
compliance = hm.compliance(state)

sim = csdl.experimental.PySimulator(recorder)
dC_dX = np.asarray(
    sim.compute_totals([compliance], [node_disp])[compliance, node_disp]
).reshape(domain.n_nodes, 3)

recorder.stop()

print(f"shape-derivative field: {dC_dX.shape}   (one gradient vector per mesh node)")

shape-derivative field: (33, 3)   (one gradient vector per mesh node)

dC_dX[k] is \(\partial(\text{compliance})/\partial\mathbf{x}_k\) at mesh node \(k\), in your node ordering. Two checks worth doing on a new model:

Stretching the plate along its axis makes it more compliant, so \(\partial C/\partial x\) should be positive at the tip and negative at the root — moving material away from the support in one case, towards it in the other.

The out-of-plane channel should be exactly zero on a flat plate: perturbing a node by \(+\delta z\) and by \(-\delta z\) gives mirror-image configurations with identical compliance, so the first-order term vanishes by symmetry. That is a genuine result, not a dead channel — give the plate a shallow curve and the same component is nonzero.

x = domain.node_coords[:, 0]
tip, root = x > LENGTH - 1e-9, x < 1e-9
print(f"mean dC/dx at the tip  : {dC_dX[tip, 0].mean():+.4e}")
print(f"mean dC/dx at the root : {dC_dX[root, 0].mean():+.4e}")
print(f"max |dC/dz|, flat plate: {np.abs(dC_dX[:, 2]).max():.4e}   (zero by reflection symmetry)")

mean dC/dx at the tip  : +2.3078e-02
mean dC/dx at the root : -2.2939e-02
max |dC/dz|, flat plate: 0.0000e+00   (zero by reflection symmetry)

# the same gradient about a shallow parabolic camber -- the z channel is live
recorder = csdl.Recorder(inline=True)
recorder.start()

curved = hm.ShellDomain(plate_mesh(nx, ny), element="CG2CG1")
camber = np.zeros((curved.n_nodes, 3))
camber[:, 2] = 0.5 * (curved.node_coords[:, 0] / LENGTH) ** 2
nd = csdl.Variable(value=camber, name="node_disp")

curved_state = hm.solve(
    curved,
    hm.isotropic(curved, E=E_VAL, nu=NU_VAL, thickness=H_VAL, density=RHO_VAL),
    hm.pressure(curved, PRESSURE_Z),
    hm.clamp(curved, where=hm.near("x", 0.0)),
    geometry=hm.geometry(curved, node_disp=nd),
)
C_curved = hm.compliance(curved_state)

sim = csdl.experimental.PySimulator(recorder)
g = np.asarray(sim.compute_totals([C_curved], [nd])[C_curved, nd]).reshape(curved.n_nodes, 3)
recorder.stop()

print(f"max |dC/dz|, cambered  : {np.abs(g[:, 2]).max():.4e}")

max |dC/dz|, cambered  : 2.0848e-03

Requirements and caveats

  • Serial mesh whose input_global_indices form a node permutation. ShellDomain enforces this for every use, not only for shape derivatives.

  • UFL patch. UFL’s CoordinateDerivative handler crashes on the shell bending term, in both UFL 2024.2 (DOLFINx 0.9) and 2026.1 (DOLFINx 0.11). Hermit applies a small runtime monkeypatch on import hermit; the installed UFL files are untouched.

  • Non-smooth stabilization scale. UFL drops the non-differentiable subgradient of CellDiameter, which appears only in the drilling / penalty stabilization scale. The effect on outputs is at the 1e-5 level.

  • Penalty conditioning. Under a very large geometry perturbation the penalty-BC system becomes ill-conditioned and finite-difference checks stop converging. That is a conditioning limit, not a bug — use the strong-BC path or a smaller step to verify.

The mesh-coordinate adjoints are finite-difference validated in tests/test_mesh_coord_deriv.py, which measures a few parts in \(10^{5}\) relative error and gates at \(3\times10^{-3}\); ex_shape_derivative_fd.py gates at \(5\times10^{-3}\). Those gates are set by finite-difference truncation and cancellation, not by the adjoint.