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 by dual contouring, with an exact derivative through every continuous stage, and SimMesh turns that surface into a tet or hex volume mesh for the solvers.

Example

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)

Extraction

GridSpec describes the sampling lattice; extract_mesh runs the whole chain and returns a Mesh. GridSpec.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 stay sharp. That path is forward-only, because the SVD truncation is not differentiable. Pass sharp=False for the Tikhonov-regularized QEF, which is smooth and differentiable everywhere.

Gradients

Which lattice edges cross the surface, which cells are active and how cells connect are discrete choices, frozen per extraction and not differentiable. Crossing positions, Hermite normals and QEF vertices are continuous and carry derivatives. An optimization loop holds topology fixed while gradients flow through node positions, and re-extracts every few steps.

To differentiate, 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. Root finding runs bisection and secant iterations on stop_gradient values, then applies one differentiable Newton correction, t = t0 - f(x(t0)) / (df/dt). At a converged root its derivative is the implicit-function-theorem value dt/dθ = -(∂f/∂θ) / (∂f/∂t), independent of how the root was found.

Bisection is also a correctness choice. Hard CSG enters the field through jnp.minimum / jnp.maximum, which makes it piecewise smooth; a bracketing method keyed on sign is unaffected by the kink, and the one-sided Hermite normals at a seam are the subgradients dual contouring needs to reconstruct a crease.

Stages

extract_mesh composes stages that are individually public, so a caller can substitute or instrument any of them. The stage-to-test map and the benchmark policy are in the meshing pipeline note.

Stage	Module	Key entry points
Lattice fitting	`meshing.bounds`	`fit_bounds`, `boundary_inside`
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`, `remove_slivers`
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 and 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

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

Diagnostics

mesh_report bundles four diagnostics: watertightness, Euler characteristic, surface deviation and triangle quality.

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 evaluates the true field at points sampled on the mesh, measuring how far the discretization strays from the geometry: a maximum absolute error of 0.0029 on a unit sphere at 26 cells.

Export

Three writers, each taking a 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, so the export is a faceted solid: one triangle per cell crossing, a 28-gon standing in for a bore. The STEP has real B-rep topology (an ADVANCED_FACE per merged polygon over its own PLANE, so a CAD kernel opens it as one closed solid rather than 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 job of the step_export plugin kind; with nothing registered for it, exports use the writers on this page and the report says so.

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 the size of those programs is decided by functionalize, not by the mesher. Three measures 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.

On scenes/end_cap.py (eight gusset ribs, three polar bolt circles, a mirrored dowel) these reduce the gradient program 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.

Volume meshes

SimMesh builds the volume mesh a solver uses; see simulation. The same declared resolution means different things per method. method="hex" only decides in/out per lattice cell, so a rib 1.3 cells thick still produces a legal, staircased HEX8 mesh. method="tet4" / "tet10" hand the dual-contour surface to TetGen as a piecewise linear complex (or to Gmsh, with mesher="gmsh", which sizes elements by the part rather than the lattice), and a PLC must be 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 rejects the surface. method="cutfem" builds no volume mesh: the lattice’s cut cells are the discretization, and the dual-contour surface is extracted only to be drawn and selected on.

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

Blends make this worse: a smooth union whose radius is comparable to the rib thickness folds the surface at the fillet even where the rib alone would mesh.

Refinement 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 triggers it:

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 solves compared across designs should read used rather than assume the declared resolution.

max_refinements=0 (default 2) fails at the declared grid instead.

Refinement failures

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 only under-resolved. It does not fix a feature that is small in absolute terms. scenes/end_cap.py’s housing, 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. Only the kind of defect changes with resolution: 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.

sdf_to_tet_mesh Newton-projects the QEF vertices onto the zero set before tetrahedralizing, and on 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 with the cells, which is part of why the ladder works.

Boundary-vertex projection

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

A SimMesh therefore 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), the same solve the zero-set protocol’s Solve performs. The mesh motion, the viewer’s feature-edge seams and the tangent-plane constructor all call this one kernel.

A scene that cannot be lowered (a bare callable field) uses the single-field projection.

Gmsh

mesher selects what 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. The example below is a plate with a through bore, a hard Difference with sharp edges:

from cadjoint.sdf.boolean import Difference
from cadjoint.sdf.primitives import Cylinder

plate = Difference((Box([0.6, 0.6, 0.4]), Cylinder(radius=0.25, height=0.9)),
                   smoothness=0.0)
mesh = SimMesh(name="plate", resolution=20, method="tet10", mesher="gmsh",
               bounds=(-0.83, -0.83, -0.63), size=(1.66, 1.66, 1.26))

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 and createGeometry split the triangle soup into one reparametrised surface per smooth region, so it meshes 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. Differences from TetGen:

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 of magnitude better. At target_size=0.16 the plate’s worst radius ratio is 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(plate)
built.mesher                      # 'gmsh'
built.owned.arity_counts()        # {0: 18281, 1: 9060, 2: 504, 3: 8}: 8 corners
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 the tet_gmsh plugin.

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.

Without it the mesh is frozen geometry: it builds, solves, inspects and exports like 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 alternative. A one-field projection of every boundary node onto the scene’s zero set would produce a gradient, but it would take crease nodes off their creases and put midsides back on chords. mesher="tetgen", the default, is differentiable throughout.

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 needs cell disambiguation.
Multi-size octree leaves (2:1 balancing, transition stitching) are not implemented.
Exports are faceted; an exact-surface STEP needs a step_export plugin.
A Gmsh mesh is frozen geometry unless a node_map plugin is installed.