Shell analysis walkthrough
This tutorial builds one shell model end to end: a mesh, a domain, boundary conditions, a material, loads, the solve, and the outputs — then makes the thickness a design variable and takes a total derivative.
It assumes the hermit conda environment (DOLFINx 0.9 or 0.11) is active. Every
number Hermit returns is a csdl.Variable, so the whole thing has to run inside a
CSDL recorder.
import numpy as np
import csdl_alpha as csdl
import hermit as hm
recorder = csdl.Recorder(inline=True)
recorder.start()
1. A mesh
A 2 × 10 cantilever plate of quadrilaterals. Hermit takes an ordinary DOLFINx mesh, so
in a real workflow this comes from your geometry pipeline — hm.read_mesh("plate.xdmf")
reads one from XDMF. Here we build it in memory to keep the tutorial self-contained.
Two traps are worth knowing about, and both are silent:
create_meshswaps its last two positional arguments between DOLFINx 0.9 and 0.11, so dispatch on the signature rather than a version string.basix orders quadrilateral vertices by tensor product —
(0,0), (1,0), (0,1), (1,1)— not cyclically around the cell. Handcreate_mesha cyclic quad and you get bowtie cells, a solve that runs happily, and answers off by orders of magnitude.
import inspect
import basix.ufl
import ufl
from mpi4py import MPI
from dolfinx.mesh import create_mesh
LENGTH, WIDTH = 10.0, 2.0
nx, ny = 20, 4
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])
# tensor-product vertex order per cell, not cyclic
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
mesh = create_mesh(MPI.COMM_WORLD, cells, el, points)
else: # DOLFINx 0.9
mesh = create_mesh(MPI.COMM_WORLD, cells, points, el)
print(f"{points.shape[0]} vertices, {cells.shape[0]} cells")
105 vertices, 80 cells
2. The domain
ShellDomain is the one-time setup: the mixed state space, the mesh tags, and the
maps between mesh-file ordering and FE dof ordering. It carries no CSDL and no
per-solve cost, so build it once and reuse it across every load case, sweep and
optimizer iteration.
The default element is CG2CG1 — quadratic displacement, linear rotation. CG1CG1 is
cheaper; CG2CR1 uses a Crouzeix–Raviart rotation that resists shear locking on thin
shells, and is triangle-only.
The mesh must be serial: ShellDomain requires its vertex numbering to be a
permutation of the file order, and raises on an MPI-partitioned mesh.
domain = hm.ShellDomain(mesh, element="CG2CG1")
print(f"{domain.n_nodes} nodes, {domain.n_cells} cells, quadrature degree {domain.quadrature_degree}")
105 nodes, 80 cells, quadrature degree 4
3. Boundary conditions
A BC is a region predicate plus a dof mask. The predicate is the raw DOLFINx entity
locator — it receives a (3, N) coordinate array and returns an (N,) boolean mask —
and hm.near / hm.on_plane build the common ones.
hm.clamp fixes all six generalized dofs; hm.pin takes a named subset of
("ux", "uy", "uz", "rx", "ry", "rz"); hm.symmetry applies a plane mask; hm.gauge
pins a single point to remove leftover rigid-body modes. All of them compose with +,
and each takes value= for a non-zero prescribed displacement.
Like the domain, boundary conditions are compiled once and reused.
bcs = hm.clamp(domain, where=hm.near("x", 0.0))
print(f"{len(bcs.penalty_terms)} penalty region(s), {len(bcs.strong)} strong constraint(s)")
1 penalty region(s), 0 strong constraint(s)
4. A material
Every material input may be a scalar (broadcast to the mesh), a per-vertex or per-cell
array, or a Field on its own FE space — they need not agree with each other. Here
they are all constants.
E_VAL, NU_VAL, H_VAL, RHO_VAL = 4.32e8, 0.0, 0.2, 1.0
material = hm.isotropic(domain, E=E_VAL, nu=NU_VAL, thickness=H_VAL, density=RHO_VAL)
print("ABD assembled on", material.A.space)
ABD assembled on ('DG', 1, (3, 3))
5. Loads
hm.pressure is a scalar load along the shell normal. hm.traction and hm.moment
take explicit global vectors per unit area — on a curved surface those are very
different loads, and self weight is a traction, not a pressure.
hm.edge_traction / hm.edge_moment / hm.edge_pressure apply the same loads per
unit length on the exterior facets a where= predicate selects, and hm.point_load
places a consistent point force or moment at any physical coordinate. Loads compose
with +, and each term keeps its own FE space all the way into the residual.
PRESSURE_Z = 2.0
loads = hm.pressure(domain, PRESSURE_Z)
# composing is just addition -- e.g. adding a downward line load at the free tip:
# loads = loads + hm.edge_traction(domain, [0.0, 0.0, -1.0], where=hm.near("x", LENGTH))
print(f"{len(loads.pressure_terms)} pressure term(s)")
1 pressure term(s)
6. Solve
hm.solve is a CSDL implicit custom operation wrapping the assembly, the direct MUMPS
factorization, and the adjoint. It returns a ShellState that back-references
everything it was given, so the output functions need only the state.
material, loads and bcs must all be built against this exact domain — DOLFINx
silently drops a boundary condition located against a second, structurally identical
function space, so Hermit checks rather than assumes.
state = hm.solve(domain, material, loads, bcs)
print("state.disp_solid:", state.disp_solid.shape, "mixed-space dofs")
state.disp_solid: (1422,) mixed-space dofs
7. Scalar outputs
Each output is a free function of the state. They are all csdl.Variables, ready to
become an objective or a constraint.
compliance = hm.compliance(state)
mass = hm.mass(state)
cg = hm.center_of_gravity(state)
energy = hm.elastic_energy(state)
Ix = WIDTH * H_VAL**3 / 12.0
eb_tip = PRESSURE_Z * WIDTH * LENGTH**4 / (8.0 * E_VAL * Ix)
print(f"compliance : {float(np.ravel(compliance.value)[0]):.6e}")
print(f"elastic energy : {float(np.ravel(energy.value)[0]):.6e} (= compliance / 2)")
print(f"mass : {float(np.ravel(mass.value)[0]):.6f} (exact {RHO_VAL * H_VAL * WIDTH * LENGTH})")
print(f"centre of gravity: {np.round(np.ravel(cg.value), 6)}")
print(f"tip deflection : {np.abs(state.disp_solid.value).max():.6e} (Euler-Bernoulli {eb_tip:.6e})")
compliance : 1.387702e-01
elastic energy : 6.938509e-02 (= compliance / 2)
mass : 4.000000 (exact 4.0)
centre of gravity: [5. 1. 0.]
tip deflection : 8.676101e-03 (Euler-Bernoulli 8.680556e-03)
8. Field outputs, and the frame they live in
Field outputs come back as Field objects, which you can ask for on any DG or CG
space and recover by L2 project, interpolate, or DG0 average.
Frames matter here. A DG strain field defaults to the element-local in-plane
frame, which follows the mesh parametrisation and is therefore a per-cell artefact: on
an unstructured mesh, plotting the raw xx component gives a patchy picture of a
smooth physical field. A CG field defaults to global Cartesian instead, because
per-cell frames are ambiguous at a shared node — and because the components are then a
full 3-D tensor, it carries six Voigt components rather than three. Use to_global()
or to_frame(...) to move between representations.
This flat, axis-aligned plate is the one case where the distinction does not show:
every element frame already coincides with global \(x\)/\(y\). ex_orientation_fields.py
runs the same comparison on a triangle mesh, where the element frames alternate cell
to cell and the three pictures are visibly different.
eps, kappa, gamma = hm.strain_fields(state) # DG2, element-local frame
print(f"curvature: {kappa.values.shape} (n_cells, [xx, yy, 2xy]), kind={kappa.kind}")
print(f" root cell, element frame : {np.round(kappa.values[0], 8)}")
print(f" root cell, global frame : {np.round(kappa.to_global().values[0], 8)}")
print(" (identical here: the plate is flat and axis-aligned)")
_, kappa_cg, _ = hm.strain_fields(state, space=("Lagrange", 1)) # smooth nodal
print(f"\n CG1 curvature field: {kappa_cg.space}, kind={kappa_cg.kind} (global Voigt-6)")
u_nodal = hm.nodal_displacement(state).value # (n_nodes, 3), mesh-file order
tip = domain.node_coords[:, 0] > LENGTH - 1e-9
print(f"\nnodal displacement: {u_nodal.shape}; mean w at the free tip = {u_nodal[tip, 2].mean():.6e}")
print(" (positive: pressure acts along +n, and the plate normal is +z)")
curvature: (80, 3) (n_cells, [xx, yy, 2xy]), kind=strain2
root cell, element frame : [ 0.00033015 -0. 0. ]
root cell, global frame : [ 0.00033015 -0. 0. ]
(identical here: the plate is flat and axis-aligned)
CG1 curvature field: ('Lagrange', 1, (6,)), kind=tensor3 (global Voigt-6)
nodal displacement: (105, 3); mean w at the free tip = 8.676101e-03
(positive: pressure acts along +n, and the plate normal is +z)
9. Fields as inputs: spaces and ordering
Field is also how you give Hermit a spatially varying input. The builder you pick
declares which ordering your array is in — that is the whole point of having several:
builder |
input |
space |
|---|---|---|
|
one value |
any |
|
|
CG1 |
|
|
DG0 |
|
a callable of the dof coordinates |
any |
|
raw coefficients, FE dof order |
any |
Getting this wrong is silent — the two orderings coincide only when the mesh permutation happens to be the identity — so say which one you mean.
# a thickness that tapers from root to tip, given per mesh vertex
taper = 0.3 - 0.02 * domain.node_coords[:, 0]
tapered = hm.isotropic(domain, E=E_VAL, nu=NU_VAL,
thickness=hm.from_nodal(domain, taper), density=RHO_VAL)
tapered_state = hm.solve(domain, tapered, loads, bcs)
print(f"uniform t=0.2 : compliance {float(np.ravel(compliance.value)[0]):.6e}, "
f"mass {float(np.ravel(mass.value)[0]):.4f}")
print(f"tapered 0.3->0.1: compliance {float(np.ravel(hm.compliance(tapered_state).value)[0]):.6e}, "
f"mass {float(np.ravel(hm.mass(tapered_state).value)[0]):.4f}")
uniform t=0.2 : compliance 1.387702e-01, mass 4.0000
tapered 0.3->0.1: compliance 6.413965e-02, mass 4.0000
10. Derivatives
Everything above is one CSDL graph, so a total derivative is a compute_totals call.
Reverse-mode sensitivities are available with respect to thickness, density, ABD, \(E\),
\(\nu\), ply angles and heights, every load coefficient, the fibre orientation, and the
mesh coordinates.
Note that the design variable has to be a csdl.Variable before it is built into a
field, so that the graph records the dependence.
t = csdl.Variable(value=H_VAL * np.ones(domain.n_nodes), name="thickness")
dv_material = hm.isotropic(domain, E=E_VAL, nu=NU_VAL,
thickness=hm.from_nodal(domain, t), density=RHO_VAL)
dv_state = hm.solve(domain, dv_material, loads, bcs)
dv_compliance = hm.compliance(dv_state)
sim = csdl.experimental.PySimulator(recorder)
dC_dt = np.asarray(sim.compute_totals([dv_compliance], [t])[dv_compliance, t]).ravel()
print(f"dC/dt: {dC_dt.shape}, all negative = {bool((dC_dt < 0).all())}")
print(f" most sensitive node at x = {domain.node_coords[np.argmin(dC_dt), 0]:.2f} (the root)")
print(f" least sensitive node at x = {domain.node_coords[np.argmax(dC_dt), 0]:.2f} (the tip)")
dC/dt: (105,), all negative = True
most sensitive node at x = 0.50 (the root)
least sensitive node at x = 10.00 (the tip)
recorder.stop()
Where to go next
Thickness optimization and shape derivatives — driving an optimizer, and differentiating with respect to the mesh coordinates.
Examples —
ex_composite_plate.py(ply-angle sweep and its derivatives),ex_fiber_orientation.pyandex_orientation_fields.py(laminate orientation and strain-field frames),ex_thickness_opt.py(a full SLSQP run).Verification benchmarks — 28 cases with published or closed-form references, each runnable on its own.