Differentiable Meshing

Dual contouring built bottom-up in JAX, so every continuous stage carries an exact derivative.

cadjoint.meshing turns a signed distance field into a surface mesh. It is not a wrapper around a black-box extractor with a gradient bolted on afterwards: every stage is written against the field callable and jax.grad of it, so the continuous half of the pipeline differentiates exactly with respect to the CAD parameters that produced the field.

The governing split is discrete topology, continuous motion. Which lattice edges cross the surface, which cells are active, and how cells connect are discrete choices — they are frozen per extraction and cannot be differentiated. Crossing positions, Hermite normals, and QEF vertices are continuous, and those carry derivatives. An optimization loop therefore holds topology fixed while gradients flow through node positions, and re-extracts every few steps.

Extract a mesh

GridSpec describes the sampling lattice; extract_mesh runs the whole chain and returns a Mesh.

import jax.numpy as jnp
from cadjoint.meshing import GridSpec, extract_mesh

def sphere_sdf(p):
    return jnp.sqrt(jnp.sum(p * p)) - 1.0

grid = GridSpec.from_bounds((-1.3, -1.3, -1.3), (2.6, 2.6, 2.6), 26)
mesh = extract_mesh(sphere_sdf, grid)

print(mesh.vertices.shape, mesh.faces.shape, mesh.quads.shape)
# (1832, 3) (3660, 3) (1830, 4)

from_bounds(bounds, size, resolution) takes a corner, an extent, and either one cell count or a per-axis triple. A Mesh is a NamedTuple of vertices and normals (JAX arrays, one per active cell), faces (triangles), quads (the native dual-contour output), and cells (the lattice index each vertex came from).

extract_mesh defaults to sharp=True, which places vertices with a rank-revealing SVD so creases and corners come out crisp. That path is forward-only — the SVD truncation is not differentiable. Pass sharp=False for the Tikhonov-regularized QEF, which is smooth and differentiable everywhere.

Differentiating through it

The differentiable recipe is explicit rather than hidden: sample the lattice, freeze the discrete structure, then rebuild only the continuous stages inside the traced function.

import jax
import jax.numpy as jnp
from cadjoint.meshing import (
    GridSpec, edge_hermite_data, find_crossing_edges,
    manifold_cell_incidence, qef_vertices, sample_grid,
)

def sphere(radius):
    return lambda p: jnp.sqrt(jnp.sum(p * p)) - radius

grid = GridSpec.from_bounds((-1.3, -1.3, -1.3), (2.6, 2.6, 2.6), 26)

# Discrete half — frozen once, outside the gradient.
values = sample_grid(sphere(1.0), grid)
edges = find_crossing_edges(values)
incidence = manifold_cell_incidence(edges, grid, values < 0.0)

# Continuous half — re-evaluated under the trace.
def loss(radius):
    hermite = edge_hermite_data(sphere(radius), grid, edges)
    vertices, _ = qef_vertices(hermite, incidence, grid)
    return jnp.sum(vertices ** 2)

print(jax.grad(loss)(jnp.float32(1.0)))  # 3638.59

The gradient is exact, not an approximation of the extractor. Root finding runs bisection and secant iterations on stop_gradient values and then applies one differentiable Newton correction, t = t0 - f(x(t0)) / (df/dt). At a converged root its derivative is precisely the implicit-function-theorem value dt/dθ = -(∂f/∂θ) / (∂f/∂t), so how the root was found does not affect the derivative that comes out.

Bisection also matters for correctness, not just speed. Hard CSG enters the field through jnp.minimum / jnp.maximum, which makes it piecewise smooth; a bracketing method keyed on sign cannot be fooled by the kink, and the one-sided Hermite normals at a seam are exactly the subgradients dual contouring needs to reconstruct a crease.

Stages

extract_mesh is a convenience over stages that are individually public, so a caller can substitute or instrument any of them.

Stage	Module	Key entry points
Hermite edge detection	`meshing.edge_detection`	`sample_grid`, `find_crossing_edges`, `edge_hermite_data`, `detect_edges`
Feature classification	`meshing.features`	`manifold_cell_incidence`, `classify_feature_cells`, `active_branches`, `detect_branch_changes`
Dual contouring	`meshing.dual_contouring`	`qef_vertices`, `sharp_qef_vertices`, `dual_faces`, `extract_mesh`
Adaptive octree	`meshing.adaptive`	`surface_cells`, `sparse_crossing_edges`
Simplification	`meshing.simplify`	`simplify_mesh`
Export	`meshing.export`	`save_obj`, `save_stl`, `save_step`, `merge_planar_faces`
Diagnostics	`meshing.diagnostics`	`mesh_report`, `surface_deviation`, `triangle_quality`

Features are classified two ways. classify_feature_cells measures normal spread with an SVD and reports each cell as FACE, CREASE, or CORNER. active_branches + detect_branch_changes are exact instead of thresholded: they record which CSG branch wins the min/max at each sample, and a cell where the winner changes is a seam, with no tolerance to tune.

Adaptive octree

Sampling the full lattice is wasteful when the surface occupies a thin shell of it. surface_cells and sparse_crossing_edges prune with an octree that descends only into blocks whose distance bound admits a crossing. The results are bit-identical to dense detection; lipschitz is the caller’s contract that the field does not grow faster than that bound (1.0 for a true distance function).

from cadjoint.meshing import GridSpec, surface_cells
import jax.numpy as jnp

grid = GridSpec.from_bounds((-1.3, -1.3, -1.3), (2.6, 2.6, 2.6), 26)
stats = {}
cells = surface_cells(lambda p: jnp.sqrt(jnp.sum(p * p)) - 1.0,
                      grid, lipschitz=1.0, stats=stats)
print(cells.shape[0], stats["evaluations"])  # 2216 5857

5857 field evaluations instead of the 19683 lattice points a dense sweep would touch, for the same set of crossings.

Export

Three writers, all taking the Mesh and a path:

import jax.numpy as jnp
from cadjoint.meshing import GridSpec, extract_mesh, save_obj, save_step, save_stl

grid = GridSpec.from_bounds((-1.3, -1.3, -1.3), (2.6, 2.6, 2.6), 26)
mesh = extract_mesh(lambda p: jnp.sqrt(jnp.sum(p * p)) - 1.0, grid)

save_obj(mesh, "part.obj")    # merge_planar=True by default
save_stl(mesh, "part.stl")    # binary=True by default
save_step(mesh, "part.step")  # AP214 faceted BREP, metres

save_obj merges coplanar triangles into n-gons first, so a flat wall exports as one face rather than a fan of triangles; merge_planar_faces exposes that directly with an angle and plane tolerance. save_step writes a CLOSED_SHELL inside a MANIFOLD_SOLID_BREP and is validated against the OCCT kernel in tests/meshing/test_step_kernel.py.

All three write the mesh as it is, and a faceted solid is the whole of what this package exports: one triangle per cell crossing, a 28-gon standing in for a bore. The STEP is real B-rep topology — an ADVANCED_FACE per merged polygon over its own PLANE, so a CAD kernel opens the file as one closed solid rather than as a triangle soup — but its faces are the mesh’s faces, and a bore arrives as the ring of narrow facets the lattice cut it into. A STEP in which a plate with a hole is seven faces instead of three thousand, the bore one CYLINDRICAL_SURFACE, is the step_export plugin kind; with nothing registered for it the writers on this page are what an export uses, and the report says so.

Checking a mesh

mesh_report bundles the four diagnostics used throughout the test suite:

import jax.numpy as jnp
from cadjoint.meshing import GridSpec, extract_mesh, mesh_report

sdf = lambda p: jnp.sqrt(jnp.sum(p * p)) - 1.0
grid = GridSpec.from_bounds((-1.3, -1.3, -1.3), (2.6, 2.6, 2.6), 26)
mesh = extract_mesh(sdf, grid)

report = mesh_report(sdf, mesh, samples=512, pairs=1024, seed=1)
print(report["watertight"], report["euler_characteristic"])
print(report["surface_deviation"])
# True 2
# {'max_abs': 0.00294, 'mean_abs': 0.00083, 'samples': 512}

Surface deviation samples the mesh and evaluates the true field there, so it measures how far the discretization strayed from the geometry it claims to represent — a max absolute error of 0.0029 on a unit sphere at 26 cells.

Tet meshes need a finer grid than hexes

The same declared resolution buys you three different things. method="hex" only has to decide in/out per lattice cell, so a rib 1.3 cells thick still produces a legal (if staircased) HEX8 mesh. method="tet4" / "tet10" hand the dual-contour surface to TetGen as a piecewise linear complex (a second route, mesher="gmsh", hands the same surface to Gmsh and sizes elements by the part rather than by the lattice), and a PLC has to be clean: watertight, manifold, and free of self-intersections. One QEF vertex per cell cannot represent two surface sheets passing through one cell, so where a wall is thinner than about two cells the extracted quads fold over each other and TetGen refuses the whole thing: method="cutfem" builds no volume mesh at all: the lattice’s cut cells are the discretisation, and the same dual-contour surface is extracted only to be drawn and selected on — see the simulation page.

RuntimeError: TetGen rejected the surface: The input surface mesh contain
self-intersections. ...

Blends make it worse, not better: a smooth union whose radius is comparable to the rib thickness folds the surface at the fillet even when the rib itself would have survived.

Reading the record

Every mesh that came through the ladder carries a refinement record, and a SimulationResult re-exports it as result.refinement (with a digest under describe()["refinement"], None when nothing was refined). A plate with a blended rib two-thirds of a cell thick is enough to see it fire:

import jax.numpy as jnp
from cadjoint.fem import SimMesh
from cadjoint.sdf.boolean import Union
from cadjoint.sdf.primitives import Box
from cadjoint.sdf.transforms import Translate

plate = Box([1.0, 1.0, 0.1])
rib = Translate(Box([0.12, 0.8, 0.4]), offset=jnp.array([0.0, 0.0, 0.4]))
part = Union((plate, rib), smoothness=0.1)

sim_mesh = SimMesh(name="ribbed", resolution=(10, 10, 6), method="tet4",
                   bounds=(-1.1, -1.1, -0.2), size=(2.2, 2.2, 1.2))
mesh = sim_mesh.build(part)

record = mesh.refinement
print(record["declared"], "->", record["used"], record["refined"])
# (10, 10, 6) -> (15, 15, 9) True

for attempt in record["attempts"]:
    print(attempt["resolution"], attempt["sharp"], attempt["outcome"])
# (10, 10, 6)  True  self-intersecting
# (10, 10, 6)  False self-intersecting
# (15, 15, 9)  True  meshed

outcome is one of "meshed", "self-intersecting" (the diagnostic fired, TetGen was never run) or "rejected" (TetGen ran and refused; the attempt also carries "error"). A refined mesh is a different mesh — more nodes, more elements, a different frozen topology — so if you are comparing solves across designs, read used rather than assuming the declared resolution.

Pass max_refinements=0 to restore the old behaviour of failing at the declared grid.

When the ladder does not help

If no rung works, the error names both ends of the ladder and the heuristic:

TetGen rejected the surface: ... The surface stays self-intersecting up to
(27, 27, 14) (declared (12, 12, 6), refined x1.5 and x2.25); the part likely
has features thinner than two cells at the declared resolution — raise the
declared resolution or use method='hex'.

Refinement fixes a surface that is merely under-resolved. It does not fix a feature that is small in absolute terms and stays small however fine the grid gets. scenes/end_cap.py’s housing is the standing example: measured over the ladder and well past it, it self-intersects at (26, 26, 13), (39, 39, 20), (59, 59, 30), (78, 78, 39), (104, 104, 52) and (130, 130, 65) alike. What changes with resolution is only which defect: at the declared (26, 26, 13) the surface is not even manifold (one edge shared by three triangles) and carries 48 folded triangle pairs; by (59, 59, 30) it is watertight and manifold, and still carries 16, all of them clustered on the snap-ring groove of the bearing seat (0.06 tall by 0.05 deep) and on the gusset-rib crests. Those two features are ~1.5 cells even at 2.25x. That part wants method="hex", a coarser blend, or a model change — not a finer grid.

Worth knowing when you read such a count: sdf_to_tet_mesh Newton-projects the QEF vertices onto the zero set before tetrahedralizing, and at a coarse grid that projection can create folds — the same end-cap surface at (26, 26, 13) has 1 self-intersecting pair before projection and 48 after, because the projection clamp is half the cell diagonal and a vertex can be dragged across a thin wall. Refining shrinks the clamp along with the cells, which is part of why the ladder works at all.

The Gmsh route

mesher picks who fills a tet4/tet10 mesh. The default, "tetgen", is everything above: the lattice sizes the elements and every node follows the design through recompute_tet_points.

mesh = SimMesh(name="body", resolution=48, method="tet10", mesher="gmsh")

mesher="gmsh" takes the same dual-contour surface, writes it as an STL, and hands it to Gmsh through the tet_mesher plugin. Gmsh’s classifySurfaces + createGeometry split the triangle soup into one reparametrised surface per smooth region — so what it meshes is a dozen faces of a part rather than three thousand facets — HXT fills the volume, and setOrder(2) puts the midside nodes on those surfaces rather than at chord midpoints. Three things change:

Element size is the part’s, not the lattice’s. The lattice is only how the surface was extracted; target_size (the grid’s smallest spacing by default) is what sizes the tets, and the refinement ladder does not run.
Quality is an order up. The plate measures a worst radius ratio of 0.293 against the DC/TetGen path’s 0.04, with the median at 0.74.
Every node is tagged with the patches that own it. assign_ownership evaluates |f_p| on each Gmsh surface entity’s nodes for every patch of the scene’s decomposition and keeps the ones that clear the bar; the mesh carries the resulting OwnedNodes record.

built = mesh.build(scene)
built.mesher                      # 'gmsh'
built.owned.arity_counts()        # {0: 1105, 1: 1494, 2: 194, 3: 8}
built.stats["blend_surfaces"]     # 0 — a hard CSG solid has no blend faces

Gmsh is GPL-2.0-or-later, so it is an optional extra (pip install 'cadjoint[gmsh]') or, in production, the cadjoint_tet_gmsh image where the licence boundary is a process boundary. Nothing under cadjoint/ imports it at module scope. See plugins.

Every vertex on every surface it touches

A boundary vertex follows the design by being re-solved onto the model, and which zero sets it is solved onto is the whole question. Solved onto the scene’s single field, a vertex sitting on a crease steps between the two faces that meet there — its position wanders along the crease and its derivative is one face’s, not the edge’s. Measured on the starter in the cut-cell study: 136 of 860 surface vertices, and 11-13 % on the objective’s gradient.

So a SimMesh lowers its scene to the node table (cadjoint.zeroset), classifies each boundary vertex onto the census surfaces it lies on — one for a face, two for an edge, three at a corner — and moves it with the minimum-norm Newton solve on all of them at once (cadjoint.zeroset.project), which is the same solve the zero-set protocol’s Solve performs. That kernel is the one projection in the repository: the mesh motion, the viewer’s feature-edge seams and the tangent-plane constructor all call it.

Where the scene cannot be lowered (a bare callable field), the single-field projection is used exactly as before.

Frozen geometry

A Gmsh mesh’s topology is decided once and held. Its positions follow the design through the node_map plugin kind — a per-arity Newton projection onto each node’s owning patches with the implicit-function adjoint, midsides re-solved on their own surfaces, interior nodes following by Laplacian relaxation — and nothing in this repository fills that kind (cadjoint.tier reports whether anything does).

So here the mesh is frozen geometry: it builds, solves, inspects and exports exactly as any other TetMesh, SimMesh.inspect() reports frozen_geometry: true, and an Optimization over a study on it is refused at declaration with one sentence naming the way out. The refusal is deliberate: cadjoint could slide every boundary node onto the scene’s zero set with a one-field projection and call that a gradient, but it would take crease nodes off their creases and put midsides back on chords, silently. mesher="tetgen" — the default — is differentiable throughout.

Program size

Every stage above evaluates the scene’s SDF under jax.jit, so the wall clock of a first mesh is mostly XLA compiling one program per stage — and how big those programs are is decided by functionalize, not by the mesher. Three things keep the tree from flattening as it is traced.

A sketch profile’s vertices are stacked into one (N, 2) array and its edge loop becomes two reductions, so a 168-vertex profile emits the same operations as a 12-vertex one.
A pattern traces its child once, over an array of instance transforms, and a node reachable from more than one parent is built once and outlined into its own StableHLO function.
Passing the parameter dicts to the jitted function rather than closing over them (functionalize_parametric) keeps values out of the program entirely, so a re-dimensioned design reuses the executable the compilation cache already holds.

Measured on scenes/end_cap.py — eight gusset ribs, three polar bolt circles, a mirrored dowel — the gradient program dropped from 1.62 MB / 18 090 operations to 0.58 MB / 5 771, and its XLA compile from 25.4 s to 0.68 s.

Known limitations

One QEF vertex per cell cannot represent two surface sheets crossing one cell face, so when a CSG seam grazes a lattice plane uniform dual contouring can emit nonmanifold edges. This is tracked as a strict xfail in tests/meshing/test_scenes.py and wants cell disambiguation. Multi-size octree leaves (2:1 balancing, transition stitching) are also open.