testkit: octree hex column in the meshing corpus

The corpus runs mesh.hex.octree's mesher on every part at 1/25 of its diagonal (the size of
the structured hex column) and reports hexahedra, the smallest and mean scaled Jacobian, the
mesh volume against SolidCraft's and the time; a part it cannot mesh gets the reason. The
column caps the elements at 500 000 so a part whose grid explodes is refused in seconds.
Loaded parts carry their bounding-box diagonal for the stages. docs/meshing/corpus.md
regenerated: 15 of 21 parts mesh, including filleted-bracket, rounded-box, sphere-on-plate
and the analytic sphere.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm

ROADMAP.md

@@ -32,7 +32,7 @@ pressure) is a degree of freedom of the existing DOF framework, not a separate s
 | M3 First deck | Sets and surfaces bound to SolidCraft face, edge, vertex and body names, resolved after meshing; materials (isotropic and orthotropic elastic, density) with a unit-system label and consistency check; sections (solid, shell, beam profiles) assigned by set; element formulation per region mapped to each solver; BCs (fixed, displacement, symmetry); loads (force and moment on reference points, pressure, gravity); reference points with kinematic and distributing coupling; steps (linear static, frequency) with explicit per-step load and BC tables and output presets; `model.check` naming every missing or bad item; `model.*` and `mesh.tag` commands (model definition: done); writers for Abaqus input (done: `deck.write`), Nastran SOL 101/103 (done: `deck.write`) and LS-DYNA implicit (done); the `.fieldcraft` file (text definition, binary mesh, geometry references reopened on load, journal; atomic saves, autosave and recovery; undo/redo on the journal: done) | in progress |
 | M4 Own solver and first results | Element library per docs/elements.md (1D and 2D first); Linear static and modal solver (tet4/10, wedge6, hex8, hex20; sparse direct, Lanczos), verified against analytic solutions (solver: done — `fieldcraft-fem` kernels, `fieldcraft-solver` with nested dissection, sparse Cholesky and block Lanczos, `solve.static`/`solve.frequency`/`solve.probe`, patch, cantilever, Lamé, Kirsch and beam-frequency tests; reference points with kinematic and distributing couplings and moments by master–slave constraint elimination, verified by rigid motion, clamping, Saint-Venant tension, bending and tube torsion; 1D elements I: rod2/rod3 trusses, point masses, rotary inertia, springs and dashpots in both solvers and all three decks, verified by the method of joints, the three-bar truss, spring–mass, two-mass chain and torsional pendulum frequencies, and damped decay; 2D elements: plane stress, plane strain and axisymmetric tri3, tri6, quad4 (incompatible modes QM6), quad8 with two degrees of freedom per node, centrifugal loads, `mesh.planar` 2D meshes of polygons and planar faces, verified by patch tests, Kirsch, Lamé (plane strain and axisymmetric), rotating discs and cylinders, and a plane stress cantilever; deck mapping to Abaqus CPS/CPE/CAX, Nastran membranes, LS-DYNA ELFORM 12/13); shells and beams (done: quad4 MITC4 and tri3 MITC3 flat shells with a Hughes–Brezzi drilling rotation, two-node Timoshenko beams with profile properties and orientations, six degrees of freedom at shell and beam nodes shared with solids and trusses, section forces and surface stresses, verified by shell patch tests, Navier's plate, cantilever and portal-frame beams, and the Scordelis–Lo and pinched-cylinder benchmarks); results display (done: `results.*` commands, deformed shape to scale, banded contours of displacement, von Mises, principal and component stresses with a legend, mode shapes with their frequencies and animation, probing a node by clicking it, reaction force and moment sums per set, stale results marked after an edit; explicit steps: field frames by slider or animation with contours per frame, energy histories plotted) | started |
 | M5 Structured hex | 2D quad meshing of faces (`fieldcraft-mesh::quad`: transfinite mapped meshes, free all-quad meshes of polygons with holes on prescribed boundary nodes, exact planar quad validity, shared `SizeField` and even edge divisions; done for planar faces); mapped hex of six-sided regions (trivariate TFI, curved faces) and sweep meshing (multi-level extrusions with imprinted cross-sections, revolutions, one-to-one sweeps through linking faces; interval assignment, hex8/hex20 on the geometry, B-rep face tags, `mesh.hex.generate`/`sweep`/`mapped`: done, 16 of 21 corpus parts all-hex); conforming interfaces between cells (done: `mesh.hex.blocks` meshes a part split with SolidCraft's Split Body into mapped blocks as one conforming mesh — composite faces for seams and tangent faces, interval assignment over all cells, shared faces meshed once, grading off holes and fillets, free-internal-face check; the filleted bracket in 30 cells, hex20 within 0.1 % of a fine tet10 tip deflection); sweepable cells in a conforming mesh and tet–hex transition; automatic decomposition into sweepable and mappable regions, proposed as SolidCraft splits | started |
-| M6 Automatic hex | Octree/grid-dual all-hex for any part; hex-dominant by tet recombination; polycube for blocky parts; smoothing, untangling, boundary padding, layers through thin walls; mesh editing and orphan-mesh import | planned |
+| M6 Automatic hex | Octree/grid all-hex for any part (first version: done — `fieldcraft-mesh::octree`, `mesh.hex.octree`: 3-refinement grid with conforming templates refined by curvature, thickness and narrow faces, a core cleared to a manifold boundary, feature map onto it, padding and a pillow per patch capturing every sharp curve, foldover-free fitting, certified validity, B-rep face tags; the corpus hex column and the filleted-bracket oracle case; next: transitions for concave refined regions so thin and curved parts stop over-refining); hex-dominant by tet recombination; polycube for blocky parts; layers through thin walls; mesh editing and orphan-mesh import | started |
 | M7 Assemblies, contact and friction | Needs M15. Part instances from SolidCraft occurrences; tie between non-matching meshes (mortar); implicit contact: surface-to-surface mortar and node-to-surface, penalty and augmented Lagrange, self-contact, automatic pair detection with a bounding-volume search shared with explicit contact; friction: Coulomb with stick–slip, static/kinetic, anisotropic, pressure-, rate- and temperature-dependent coefficients; contact output (pressure, slip, status); thermal contact conductance and frictional heating; wear (Archard); bolt pretension, connectors with behaviour (elastic, damped, stops, locks); amplitudes, multiple steps, restart; nonlinear steps in all three decks. Verified by Hertz (sphere and cylinder), the Cattaneo–Mindlin partial-slip problem, the sliding block, NAFEMS contact benchmarks | planned |
 | M8 Materials, explicit, shells, composites | Plasticity, hyperelasticity, damage data per solver; explicit dynamics decks (Abaqus/Explicit, LS-DYNA); midsurface shells; composite layups and material orientations, continuum shells on swept hex; hybrid (mixed u–p) hex8/hex20/tet10 for nearly incompressible materials (done); heat transfer (solids done: steady and θ-method transient conduction with flux, convection and sources, thermal strain in static steps (solids and 2D) and explicit dynamics (with plasticity), coupled temperature–displacement, `solve.heat`; deck output, 2D/shell/beam conduction planned) | started |
 | M9 Post-processing and deck readers | Results readers for public formats (Nastran OP2, LS-DYNA d3plot, Abaqus `.fil`); XY and path plots, cuts, display groups, free-body forces; deck readers that keep unknown keywords; mesh convergence studies | planned |
@@ -60,7 +60,7 @@ In progress, in this order of landing:
 
 1. Tet quality: sliver repair and Delaunay refinement to a size field, so curved and filleted parts
    (the filleted bracket) solve on tets (M2).
-2. Octree all-hex for any part (M6), the project's main meshing goal.
+2. Octree all-hex (M6): concave transitions (Schneiders' modified templates or Livesu's dual schemes) so curved and thin-walled parts mesh without filling their hull with fine cells.
 3. The analysis video regenerated with the setup's symbols (done: support, load, reference-point
    and coupling symbols in the viewport and headless renders), captions and an orbit (screenshot
    camera and the ribbon following the command: done).

crates/engine/src/catalogue.rs

@@ -281,7 +281,13 @@ pub static CATALOGUE: &[Capability] = &[
     ),
     soon("mesh.hex.decompose", "Decompose", Mesh, "Structured Hex", "M5", "Find sweepable and mappable regions and propose the splits"),
     soon("mesh.transition", "Transition", Mesh, "Structured Hex", "M5", "Conforming tet–hex interfaces with pyramids"),
-    soon("mesh.hex.octree", "All-Hex", Mesh, "Automatic Hex", "M6", "Octree dual all-hex mesh of any part"),
+    ready(
+        "mesh.hex.octree",
+        "Octree Hex",
+        Mesh,
+        "Automatic Hex",
+        "Unstructured all-hex mesh of any part: adaptive grid, boundary layers on every sharp edge, untangled and proven valid",
+    ),
     soon("mesh.hex.dominant", "Hex-Dominant", Mesh, "Automatic Hex", "M6", "Hex-dominant mesh by recombining tets"),
     soon("mesh.hex.polycube", "Polycube", Mesh, "Automatic Hex", "M6", "Polycube hex mesh of blocky parts"),
     soon("mesh.smooth", "Smooth", Mesh, "Improve", "M6", "Smooth node positions, keeping every element valid"),

crates/engine/src/cmd/hex.rs

@@ -3,7 +3,7 @@
 use std::borrow::Cow;
 
 use fieldcraft_geometry::{FaceSurface, Shape, SurfaceMesh};
-use fieldcraft_mesh::hex::{Cell, HexMesh, HexMethod, HexOptions, MAX_CELLS, Surface, hex_mesh, hex_mesh_cells};
+use fieldcraft_mesh::hex::{BoundaryFace, Cell, HexMesh, HexMethod, HexOptions, MAX_CELLS, Surface, hex_mesh, hex_mesh_cells};
 use fieldcraft_mesh::size::UniformSize;
 use fieldcraft_mesh::surface::BoundarySurface;
 use fieldcraft_model::{GeometryRef, Model};
@@ -28,13 +28,13 @@ const MIN_RELATIVE_TOLERANCE: f64 = 1e-5;
 const DEFAULT_RELATIVE_SIZE: f64 = 1.0 / 25.0;
 const DEFAULT_RELATIVE_TOLERANCE: f64 = 1e-3;
 
-/// One body to mesh.
-struct Input<'a> {
-    name: String,
-    mesh: Cow<'a, SurfaceMesh>,
+/// One body to mesh (also the input of `mesh.hex.octree`).
+pub(super) struct Input<'a> {
+    pub(super) name: String,
+    pub(super) mesh: Cow<'a, SurfaceMesh>,
     surfaces: Vec<Option<Surface>>,
     face_names: Vec<String>,
-    diagonal: f64,
+    pub(super) diagonal: f64,
 }
 
 fn analytic(s: FaceSurface) -> Surface {
@@ -47,7 +47,7 @@ fn analytic(s: FaceSurface) -> Surface {
     }
 }
 
-fn number(p: &Value, key: &str) -> Result<Option<f64>> {
+pub(super) fn number(p: &Value, key: &str) -> Result<Option<f64>> {
     match p.get(key) {
         None => Ok(None),
         Some(v) => v
@@ -73,7 +73,7 @@ fn shape_input(shape: &Shape, tolerance: Option<f64>) -> Result<Input<'static>>
     Ok(Input { name: shape.name.clone(), mesh: Cow::Owned(shape.tessellate(tol)?), surfaces, face_names, diagonal: diag })
 }
 
-fn inputs<'a>(s: &'a Session, p: &Value) -> Result<Vec<Input<'a>>> {
+pub(super) fn inputs<'a>(s: &'a Session, p: &Value) -> Result<Vec<Input<'a>>> {
     let tolerance = number(p, "tolerance")?;
     match (p.get("body"), p.get("surface")) {
         (Some(b), None) => {
@@ -220,7 +220,7 @@ fn hex_blocks(s: &mut Session, p: &Value) -> Result<Value> {
         model.set_element(label(e)?, m.kind, nodes)?;
         model.tag_element(label(e)?, body_tags.get(cell).copied().ok_or_else(|| EngineError::Params("element without a cell".into()))?)?;
     }
-    tag_faces(&mut model, &inputs, |e| cm.element_cell.get(e as usize).map_or(0, |&c| c as usize), m, 1)?;
+    tag_faces(&mut model, &inputs, |e| cm.element_cell.get(e as usize).map_or(0, |&c| c as usize), &m.boundary, 1)?;
     s.model = model;
     let blocks: Vec<Value> =
         inputs.iter().zip(&cm.blocks).map(|(b, d)| json!({"body": b.name, "divisions": d, "elements": d.iter().product::<usize>()})).collect();
@@ -318,7 +318,7 @@ fn generate(s: &mut Session, p: &Value, default: HexMethod, allowed: &[HexMethod
             model.tag_element(element_label, body_tag)?;
             element_label = element_label.checked_add(1).ok_or_else(|| EngineError::Params("too many elements".into()))?;
         }
-        tag_faces(&mut model, std::slice::from_ref(b), |_| 0, m, first_element)?;
+        tag_faces(&mut model, std::slice::from_ref(b), |_| 0, &m.boundary, first_element)?;
         report.push(json!({
             "body": b.name,
             "method": m.report.method.name(),
@@ -338,9 +338,15 @@ fn generate(s: &mut Session, p: &Value, default: HexMethod, allowed: &[HexMethod
 
 /// Tag every boundary element face with the B-rep face it lies on, of the body
 /// (`bodies[cell(element)]`) its element meshes.
-fn tag_faces(model: &mut Model, bodies: &[Input<'_>], cell: impl Fn(u32) -> usize, m: &HexMesh, first_element: u32) -> Result<()> {
+pub(super) fn tag_faces(
+    model: &mut Model,
+    bodies: &[Input<'_>],
+    cell: impl Fn(u32) -> usize,
+    boundary: &[BoundaryFace],
+    first_element: u32,
+) -> Result<()> {
     let mut tags: std::collections::BTreeMap<(usize, u32), u32> = std::collections::BTreeMap::new();
-    for f in &m.boundary {
+    for f in boundary {
         let c = cell(f.element);
         let tag = match tags.get(&(c, f.brep_face)) {
             Some(&t) => t,

crates/engine/src/cmd/mesh.rs

@@ -93,7 +93,7 @@ const TET_MAX_STEINER: u64 = 10_000_000;
 const TET_MAX_TETS: u64 = 50_000_000;
 
 /// An optional whole-number parameter in `range`.
-fn optional_count(p: &Value, key: &str, default: u64, range: std::ops::RangeInclusive<u64>) -> Result<usize> {
+pub(super) fn optional_count(p: &Value, key: &str, default: u64, range: std::ops::RangeInclusive<u64>) -> Result<usize> {
     let v = match p.get(key) {
         None => default,
         Some(v) => v

crates/engine/src/cmd/mod.rs

@@ -9,6 +9,7 @@ mod geometry;
 mod hex;
 pub(crate) mod mesh;
 mod model;
+pub(crate) mod octree;
 mod results;
 pub(crate) mod solve;
 mod tags;
@@ -36,6 +37,7 @@ pub static ALL: &[CommandSpec] = &[
     mesh::SURFACE_CHECK,
     mesh::PLANAR,
     mesh::TET_GENERATE,
+    octree::HEX_OCTREE,
     block::BLOCK,
     hex::HEX_GENERATE,
     hex::HEX_SWEEP,

crates/engine/src/cmd/octree.rs

@@ -0,0 +1,177 @@
+//! Unstructured all-hexahedral meshing (`mesh.hex.octree`) of the open design's bodies or of
+//! imported surfaces.
+
+use fieldcraft_mesh::octree::{OctreeHexMesh, OctreeOptions, octree_hex_mesh};
+use fieldcraft_mesh::size::UniformSize;
+use fieldcraft_mesh::surface::BoundarySurface;
+use fieldcraft_model::{ElementKind, GeometryRef, MAX_ELEMENTS};
+use serde_json::{Value, json};
+
+use super::hex::{Input, inputs, number, tag_faces};
+use super::mesh::optional_count;
+use crate::{CommandKind, CommandSpec, EngineError, Result, Session};
+
+pub const HEX_OCTREE: CommandSpec = CommandSpec {
+    id: "mesh.hex.octree",
+    label: "Octree Hex Mesh",
+    params: r#"{"body": "Body1" (default: every body of the open design) or "surface": 0 (an imported surface), "size": 2.0 (element size in the interior, mm; default 1/25 of the bounding-box diagonal), "tolerance": 0.05 (chordal tolerance of the tessellation, mm; default 1e-3 of the diagonal), "max_depth": 3 (refinement levels below the size, each a third, 0..=5), "max_elements": 2000000 (1..=10000000; default the session limit), "attempts": 4 (builds, refining where one failed, 1..=8), "max_sweeps": 300 (optimiser sweeps, 0..=10000), and the boundary options of mesh.tet.generate ("weld", "weld_distance", "feature_angle")} — replaces the mesh with an unstructured all-hex8 mesh of any part: an adaptive 3-refinement grid with conforming transition templates, a core inside the part, two boundary layers that capture every B-rep face boundary and sharp edge, and untangling. Every element is proven valid; boundary element faces are tagged with their B-rep faces (named as in geometry.bodies; by index when a face has no name), so sets by geometry resolve. Reports per body the elements, quality, volume and grid statistics"#,
+    kind: CommandKind::Model,
+    run: hex_octree,
+};
+
+/// Smallest element size, relative to the bounding-box diagonal (the grid's cell count grows
+/// with the cube of its inverse).
+const MIN_RELATIVE_SIZE: f64 = 1e-3;
+/// Default element size, relative to the bounding-box diagonal.
+const DEFAULT_RELATIVE_SIZE: f64 = 1.0 / 25.0;
+/// Caps on the options.
+const MAX_DEPTH: u64 = 5;
+const MAX_ATTEMPTS: u64 = 8;
+const MAX_SWEEPS: u64 = 10_000;
+const MAX_HEX_ELEMENTS: u64 = 10_000_000;
+/// The session default of `max_elements` ([`crate::Limits::max_hex_elements`]).
+pub(crate) const DEFAULT_MAX_ELEMENTS: usize = 2_000_000;
+
+fn hex_octree(s: &mut Session, p: &Value) -> Result<Value> {
+    let size = number(p, "size")?;
+    let defaults = OctreeOptions::default();
+    let options = OctreeOptions {
+        max_depth: u8::try_from(optional_count(p, "max_depth", u64::from(defaults.max_depth), 0..=MAX_DEPTH)?).unwrap_or(defaults.max_depth),
+        max_elements: optional_count(p, "max_elements", s.limits.max_hex_elements as u64, 1..=MAX_HEX_ELEMENTS)?.min(MAX_ELEMENTS),
+        attempts: optional_count(p, "attempts", defaults.attempts as u64, 1..=MAX_ATTEMPTS)?,
+        max_sweeps: optional_count(p, "max_sweeps", defaults.max_sweeps as u64, 0..=MAX_SWEEPS)?,
+        ..defaults
+    };
+    let bodies = inputs(s, p)?;
+    let mut meshes: Vec<(Input, OctreeHexMesh)> = Vec::with_capacity(bodies.len());
+    let mut total = 0usize;
+    for b in bodies {
+        if !(b.diagonal > 0.0 && b.diagonal.is_finite()) {
+            return Err(EngineError::Params(format!("body `{}` has no extent", b.name)));
+        }
+        let h = size.unwrap_or(DEFAULT_RELATIVE_SIZE * b.diagonal);
+        if h < MIN_RELATIVE_SIZE * b.diagonal {
+            return Err(EngineError::Params(format!("`size` must be at least {} mm for body `{}`", MIN_RELATIVE_SIZE * b.diagonal, b.name)));
+        }
+        let surface_options = super::mesh::surface_options(p, &b.mesh)?;
+        let surface = BoundarySurface::new(&b.mesh.positions, &b.mesh.triangles, &b.mesh.tri_face, &surface_options)
+            .map_err(|e| EngineError::Params(format!("body `{}`: {e}", b.name)))?;
+        let left = OctreeOptions { max_elements: options.max_elements.saturating_sub(total).max(1), ..options.clone() };
+        let mesh = octree_hex_mesh(&surface, &UniformSize(h), &left).map_err(|e| EngineError::Params(format!("body `{}`: {e}", b.name)))?;
+        total += mesh.hexes.len();
+        meshes.push((b, mesh));
+    }
+    // Replace the mesh: definitions stay, nodes, elements and geometry tags are new.
+    let mut model = s.model.without_mesh();
+    let mut report = Vec::with_capacity(meshes.len());
+    let (mut node_label, mut element_label) = (1u32, 1u32);
+    for (b, m) in &meshes {
+        let first_node = node_label;
+        for xyz in &m.nodes {
+            model.set_node(node_label, *xyz)?;
+            node_label = node_label.checked_add(1).ok_or_else(|| EngineError::Params("too many nodes".into()))?;
+        }
+        let body_tag = model.geometry_tag(&GeometryRef::body(b.name.as_str())?)?;
+        let first_element = element_label;
+        for hx in &m.hexes {
+            model.set_element(element_label, ElementKind::Hex8, hx.iter().map(|&i| first_node + i).collect())?;
+            model.tag_element(element_label, body_tag)?;
+            element_label = element_label.checked_add(1).ok_or_else(|| EngineError::Params("too many elements".into()))?;
+        }
+        tag_faces(&mut model, std::slice::from_ref(b), |_| 0, &m.boundary, first_element)?;
+        let r = &m.report;
+        report.push(json!({
+            "body": b.name,
+            "kind": "hex8",
+            "elements": m.hexes.len(),
+            "nodes": m.nodes.len(),
+            "boundary_faces": m.boundary.len(),
+            "quality": {"min_scaled_jacobian": r.min_scaled_jacobian, "mean_scaled_jacobian": r.mean_scaled_jacobian},
+            "volume": r.volume,
+            "surface_volume": r.surface_volume,
+            "grid": {"cells": r.grid_cells, "core_elements": r.core_elements, "layer_elements": r.layer_elements, "frame": r.frame},
+            "features": {"corners": r.corners, "curves": r.curves, "patches": r.patches},
+            "attempts": r.attempts,
+            "retries": r.retries,
+            "sweeps": r.sweeps,
+        }));
+    }
+    s.model = model;
+    Ok(json!({"bodies": report, "nodes": node_label - 1, "elements": element_label - 1}))
+}
+
+#[cfg(test)]
+mod tests {
+    use std::sync::Arc;
+
+    use fieldcraft_geometry::open_design;
+    use fieldcraft_geometry::sample::{Plate, plate_design};
+    use serde_json::json;
+
+    use crate::Session;
+
+    fn opened() -> Session {
+        let plate = Plate::default();
+        Session { design: Some(Arc::new(open_design(&plate_design(&plate).unwrap()).unwrap())), ..Session::default() }
+    }
+
+    #[test]
+    fn plate_with_a_hole_meshes_all_hex_and_face_sets_resolve() {
+        let mut s = opened();
+        let r = s.execute("mesh.hex.octree", &json!({"body": "Body1", "size": 5.0, "max_depth": 2})).unwrap();
+        let body = &r["bodies"][0];
+        let n = body["elements"].as_u64().unwrap() as usize;
+        assert!(n > 50, "{r}");
+        assert!(body["quality"]["min_scaled_jacobian"].as_f64().unwrap() > 0.0, "{r}");
+        assert_eq!(s.model.element_counts(), vec![(fieldcraft_model::ElementKind::Hex8, n)]);
+        let (v, sv) = (body["volume"].as_f64().unwrap(), body["surface_volume"].as_f64().unwrap());
+        assert!((v - sv).abs() < 0.02 * sv, "{r}");
+        // The certified quality report agrees.
+        let q = s.execute("mesh.quality", &json!({})).unwrap();
+        assert_eq!(q["valid"].as_u64().unwrap() as usize, n);
+        // Every boundary face is tagged: the skin of the body set; each named face resolves.
+        s.execute("model.set.create", &json!({"name": "skin", "geometry": [{"body": "Body1"}]})).unwrap();
+        let skin = s.execute("model.set.resolve", &json!({"name": "skin"})).unwrap();
+        assert_eq!(skin["faces"]["count"], body["boundary_faces"], "{skin}");
+        let bodies = s.execute("geometry.bodies", &json!({})).unwrap();
+        for f in bodies[0]["faces"].as_array().unwrap() {
+            let name = f["name"].as_str().unwrap();
+            s.execute("model.set.create", &json!({"name": "f", "geometry": [{"body": "Body1", "face": name}]})).unwrap();
+            let set = s.execute("model.set.resolve", &json!({"name": "f"})).unwrap();
+            assert!(set["faces"]["count"].as_u64().unwrap() > 0, "face {name}: {set}");
+            s.execute("model.set.delete", &json!({"name": "f"})).unwrap();
+        }
+    }
+
+    #[test]
+    fn hostile_params_are_errors_and_change_nothing() {
+        let mut s = opened();
+        let before = s.model.clone();
+        for p in [
+            json!({"body": "Body1", "size": 0}),
+            json!({"body": "Body1", "size": -1}),
+            json!({"body": "Body1", "size": 1e-300}),
+            json!({"body": "Body1", "size": "big"}),
+            json!({"body": "Body1", "tolerance": 1e-300}),
+            json!({"body": "Body1", "tolerance": -1}),
+            json!({"body": "Body1", "max_depth": 9}),
+            json!({"body": "Body1", "max_depth": -1}),
+            json!({"body": "Body1", "max_elements": 0}),
+            json!({"body": "Body1", "max_elements": 1e300}),
+            json!({"body": "Body1", "attempts": 0}),
+            json!({"body": "Body1", "attempts": 100}),
+            json!({"body": "Body1", "max_sweeps": "many"}),
+            json!({"body": "Body1", "size": 5.0, "max_depth": 2, "max_elements": 10}),
+            json!({"body": "Nope"}),
+            json!({"body": 7}),
+            json!({"surface": 0}),
+            json!({"surface": -1}),
+            json!({"body": "Body1", "surface": 0}),
+        ] {
+            let r = s.execute("mesh.hex.octree", &p);
+            assert!(r.is_err(), "{p}: {r:?}");
+            assert_eq!(s.model, before, "{p}");
+        }
+        assert!(Session::default().execute("mesh.hex.octree", &json!({})).is_err());
+    }
+}

crates/engine/src/lib.rs

@@ -135,6 +135,8 @@ pub struct Limits {
     pub max_steiner: usize,
     /// `mesh.tet.generate`'s `max_tets`.
     pub max_tets: usize,
+    /// `mesh.hex.octree`'s `max_elements`.
+    pub max_hex_elements: usize,
     /// Threads the solvers use (the `solve.*` commands' `threads`): by default
     /// `FIELDCRAFT_THREADS` or the available cores, 1 to [`fieldcraft_parallel::MAX_THREADS`].
     /// Results do not depend on it.
@@ -144,7 +146,12 @@ pub struct Limits {
 impl Default for Limits {
     fn default() -> Self {
         let cdt = fieldcraft_mesh::cdt::CdtOptions::default();
-        Limits { max_steiner: cdt.max_steiner, max_tets: cdt.max_tets, threads: fieldcraft_parallel::default_threads() }
+        Limits {
+            max_steiner: cdt.max_steiner,
+            max_tets: cdt.max_tets,
+            max_hex_elements: cmd::octree::DEFAULT_MAX_ELEMENTS,
+            threads: fieldcraft_parallel::default_threads(),
+        }
     }
 }
 
@@ -471,9 +478,11 @@ mod tests {
             json!({"geometry": {"body": "Bar", "vertex": "v"}, "nodes": [1], "edges": [{"element": 1, "edge": 0}]}),
         ];
         // Small work caps, so input that would mesh up to the production caps (needle triangles
-        // need Steiner points all along them) is refused at once. The caps themselves are what
+        // need Steiner points all along them, an octree hex mesh is fitted and untangled) is
+        // refused at once. The caps themselves are what
         // `tet_generate_honours_the_session_caps` checks.
-        let capped = || Session { limits: Limits { max_steiner: 64, max_tets: 100_000, ..Limits::default() }, ..Session::default() };
+        let capped =
+            || Session { limits: Limits { max_steiner: 64, max_tets: 100_000, max_hex_elements: 500, ..Limits::default() }, ..Session::default() };
         for spec in commands() {
             for p in &hostile {
                 let mut s = capped();

crates/mesh/src/hex/mod.rs

@@ -37,7 +37,7 @@
 mod arrangement;
 mod brep;
 mod cells;
-mod geom;
+pub(crate) mod geom;
 mod interval;
 mod layered;
 mod mapped;
@@ -583,7 +583,7 @@ fn finish(cells: &Cells, raw: Raw, method: HexMethod) -> Result<HexMesh> {
 
 /// `∫ x n_x dA` over a quad4 or quad8 face (outward normal): its contribution to the enclosed
 /// volume. The 3 × 3 Gauss rule integrates it exactly (degree at most 5 per direction).
-fn face_flux(face: &[P]) -> f64 {
+pub(crate) fn face_flux(face: &[P]) -> f64 {
     const G: [(f64, f64); 3] = [(-0.774_596_669_241_483_4, 5.0 / 9.0), (0.0, 8.0 / 9.0), (0.774_596_669_241_483_4, 5.0 / 9.0)];
     let mut sum = 0.0;
     for (r, wr) in G {

crates/mesh/src/lib.rs

@@ -12,6 +12,7 @@ pub mod cdt;
 pub mod delaunay;
 pub mod hex;
 pub mod join;
+pub mod octree;
 pub mod planar;
 pub mod quad;
 pub mod quality;

crates/mesh/src/octree/geom.rs

@@ -0,0 +1,636 @@
+//! The boundary as the hex mesher sees it: a triangle tree for nearest-point, ray and box
+//! queries, the exact inside test, and the feature graph (corners, curves, patches).
+//!
+//! # Feature graph
+//! Every feature edge of the [`BoundarySurface`] (between two B-rep faces, or a crease inside a
+//! face) whose dihedral angle reaches the sharp angle is kept; smooth face boundaries (a fillet
+//! meeting its faces tangentially) are not features, the faces on either side share a patch.
+//! Neither are blades and slits, edges whose normals turn by more than the blade angle (a 3°
+//! wedge, a sphere touching a plane): no valid hexahedra can follow them, so the mesh rounds
+//! them within the boundary tolerance. *Corners* are vertices where the number of feature edges is not two, or where a
+//! curve turns by more than [`CORNER_TURN`]; *curves* are the chains of feature edges between
+//! corners, or closed loops; *patches* are the connected sets of triangles between feature edges,
+//! so every patch lies in one B-rep face and is bounded by curves. This is the partition the hex
+//! mesh's boundary must reproduce (Gao, Shen, Panozzo 2019, §3.3).
+
+use std::collections::BTreeMap;
+
+pub(crate) use crate::hex::geom::closest_on_triangle;
+use crate::predicates::{Axis, Point3, Sign, orient2d_projected, orient3d, snap_to_domain};
+use crate::surface::BoundarySurface;
+use crate::vec3::{P, cross, dot, norm, sub};
+use crate::{Error, Result};
+
+/// A degree-2 feature vertex whose curve turns by more than this (degrees) is a corner.
+pub(crate) const CORNER_TURN: f64 = 30.0;
+/// Triangles per tree leaf.
+const LEAF: usize = 4;
+
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub(crate) struct Aabb {
+    pub lo: P,
+    pub hi: P,
+}
+
+impl Aabb {
+    pub const EMPTY: Aabb = Aabb { lo: [f64::INFINITY; 3], hi: [f64::NEG_INFINITY; 3] };
+    pub fn add(&mut self, p: P) {
+        for ((lo, hi), x) in self.lo.iter_mut().zip(&mut self.hi).zip(p) {
+            *lo = lo.min(x);
+            *hi = hi.max(x);
+        }
+    }
+    fn join(&mut self, o: &Aabb) {
+        self.add(o.lo);
+        self.add(o.hi);
+    }
+    pub fn overlaps(&self, o: &Aabb) -> bool {
+        (0..3).all(|k| self.lo[k] <= o.hi[k] && o.lo[k] <= self.hi[k])
+    }
+    /// Squared distance from a point to the box (0 inside).
+    fn dist2(&self, p: P) -> f64 {
+        (0..3).map(|k| (self.lo[k] - p[k]).max(0.0).max(p[k] - self.hi[k])).map(|d| d * d).sum()
+    }
+    /// Whether the segment `a + t (b − a)`, `t ∈ [0, 1]`, may meet the box (slab test, widened by
+    /// a relative margin so rounding never drops a hit).
+    fn meets_segment(&self, a: P, b: P) -> bool {
+        let (mut t0, mut t1) = (0.0f64, 1.0f64);
+        for k in 0..3 {
+            let pad = 1e-9 * (self.hi[k] - self.lo[k]).abs().max(1e-300) + 1e-12 * (a[k].abs() + b[k].abs());
+            let (lo, hi) = (self.lo[k] - pad, self.hi[k] + pad);
+            let d = b[k] - a[k];
+            if d == 0.0 {
+                if a[k] < lo || a[k] > hi {
+                    return false;
+                }
+                continue;
+            }
+            let (mut u, mut v) = ((lo - a[k]) / d, (hi - a[k]) / d);
+            if u > v {
+                std::mem::swap(&mut u, &mut v);
+            }
+            t0 = t0.max(u);
+            t1 = t1.min(v);
+            if t0 > t1 {
+                return false;
+            }
+        }
+        true
+    }
+}
+
+#[derive(Debug, Clone, Copy)]
+struct Node {
+    bounds: Aabb,
+    children: Option<usize>,
+    start: usize,
+    end: usize,
+}
+
+/// A median-split tree over triangle boxes.
+#[derive(Debug, Clone)]
+struct Tree {
+    nodes: Vec<Node>,
+    order: Vec<u32>,
+}
+
+impl Tree {
+    fn build(boxes: &[Aabb]) -> Tree {
+        let mut order: Vec<u32> = (0..boxes.len() as u32).collect();
+        let empty = Node { bounds: Aabb::EMPTY, children: None, start: 0, end: 0 };
+        let mut nodes = vec![empty];
+        let mut work = vec![(0usize, 0usize, boxes.len())];
+        while let Some((ni, start, end)) = work.pop() {
+            let Some(slice) = order.get_mut(start..end) else { continue };
+            let (mut bounds, mut centres) = (Aabb::EMPTY, Aabb::EMPTY);
+            for &t in slice.iter() {
+                if let Some(b) = boxes.get(t as usize) {
+                    bounds.join(b);
+                    centres.add([0, 1, 2].map(|k| b.lo[k] + b.hi[k]));
+                }
+            }
+            let mut node = Node { bounds, children: None, start, end };
+            if end - start > LEAF {
+                let ext = [0, 1, 2].map(|k| centres.hi[k] - centres.lo[k]);
+                let axis = if ext[0] >= ext[1] && ext[0] >= ext[2] {
+                    0
+                } else if ext[1] >= ext[2] {
+                    1
+                } else {
+                    2
+                };
+                let key = |t: &u32| boxes.get(*t as usize).map_or(0.0, |b| b.lo[axis] + b.hi[axis]);
+                let mid = (end - start) / 2;
+                slice.select_nth_unstable_by(mid, |a, b| key(a).total_cmp(&key(b)).then(a.cmp(b)));
+                let l = nodes.len();
+                nodes.push(empty);
+                nodes.push(empty);
+                node.children = Some(l);
+                work.push((l, start, start + mid));
+                work.push((l + 1, start + mid, end));
+            }
+            if let Some(n) = nodes.get_mut(ni) {
+                *n = node;
+            }
+        }
+        Tree { nodes, order }
+    }
+
+    fn leaf(&self, n: &Node) -> &[u32] {
+        self.order.get(n.start..n.end).unwrap_or(&[])
+    }
+}
+
+pub(crate) fn lerp(a: P, b: P, t: f64) -> P {
+    [0, 1, 2].map(|k| a[k] + (b[k] - a[k]) * t)
+}
+
+/// The closest point of segment `ab` to `p` and its parameter.
+pub(crate) fn closest_on_segment(p: P, a: P, b: P) -> (P, f64) {
+    let ab = sub(b, a);
+    let l2 = dot(ab, ab);
+    let t = if l2 > 0.0 { (dot(sub(p, a), ab) / l2).clamp(0.0, 1.0) } else { 0.0 };
+    (lerp(a, b, t), t)
+}
+
+pub(crate) fn dist2(a: P, b: P) -> f64 {
+    let d = sub(a, b);
+    dot(d, d)
+}
+
+/// A chain of feature edges.
+#[derive(Debug, Clone, PartialEq)]
+pub(crate) struct Curve {
+    /// Surface vertices in order; a closed loop repeats its first vertex at the end.
+    pub vertices: Vec<u32>,
+    pub closed: bool,
+    /// The patches on either side (sorted, deduplicated: one entry for a crease seam).
+    pub patches: Vec<u32>,
+}
+
+#[derive(Debug, Clone, PartialEq)]
+pub(crate) struct Patch {
+    pub triangles: Vec<u32>,
+    /// The B-rep faces of its triangles (sorted, deduplicated; one for a well-formed boundary).
+    pub faces: Vec<u32>,
+    pub area: f64,
+    /// Twice the area over the length of the bounding curves: the width of a strip, the
+    /// diameter of a disc, `∞` for a closed surface without curves.
+    pub width: f64,
+}
+
+/// The boundary with its queries and feature graph.
+pub(crate) struct Geometry<'a> {
+    pub surface: &'a BoundarySurface,
+    tree: Tree,
+    pub boxes: Vec<Aabb>,
+    pub normals: Vec<P>,
+    pub tri_patch: Vec<u32>,
+    pub patches: Vec<Patch>,
+    pub curves: Vec<Curve>,
+    /// Corner vertices, ascending.
+    pub corners: Vec<u32>,
+    pub bounds: Aabb,
+    /// Far enough to leave the boundary from any point of its box.
+    far: f64,
+}
+
+impl<'a> Geometry<'a> {
+    pub fn new(surface: &'a BoundarySurface, sharp_angle: f64, blade_angle: f64) -> Result<Geometry<'a>> {
+        let pts = surface.points();
+        let tris = surface.triangles();
+        let xyz = |v: u32| pts.get(v as usize).copied().unwrap_or([0.0; 3]);
+        let mut bounds = Aabb::EMPTY;
+        let mut boxes = Vec::with_capacity(tris.len());
+        let mut normals = Vec::with_capacity(tris.len());
+        for t in tris {
+            let mut b = Aabb::EMPTY;
+            for &v in t {
+                b.add(xyz(v));
+            }
+            bounds.join(&b);
+            boxes.push(b);
+            let n = cross(sub(xyz(t[1]), xyz(t[0])), sub(xyz(t[2]), xyz(t[0])));
+            let l = norm(n);
+            normals.push(if l > 0.0 { n.map(|c| c / l) } else { [0.0; 3] });
+        }
+        if tris.is_empty() || !(0..3).all(|k| bounds.hi[k] >= bounds.lo[k]) {
+            return Err(Error::Invalid("the boundary has no triangles".into()));
+        }
+        let tree = Tree::build(&boxes);
+        let diag = norm(sub(bounds.hi, bounds.lo));
+        let far = 4.0 * diag.max(1e-9) + (0..3).map(|k| bounds.lo[k].abs().max(bounds.hi[k].abs())).fold(0.0, f64::max);
+
+        // Feature edges and corners.
+        // A boundary between two B-rep faces is kept or dropped as a whole, by the median angle
+        // along it, so tessellation noise never breaks a curve into pieces; creases inside a
+        // face go edge by edge.
+        let face_of = |t: u32| surface.tri_face().get(t as usize).copied().unwrap_or(0);
+        let mut pair_angles: BTreeMap<(u32, u32), Vec<f64>> = BTreeMap::new();
+        for e in surface.feature_edges() {
+            let (a, b) = (face_of(e.triangles[0]), face_of(e.triangles[1]));
+            if a != b {
+                pair_angles.entry((a.min(b), a.max(b))).or_default().push(e.angle);
+            }
+        }
+        let pair_kept: BTreeMap<(u32, u32), bool> = pair_angles
+            .into_iter()
+            .map(|(k, mut v)| {
+                v.sort_by(f64::total_cmp);
+                let m = v.get(v.len() / 2).copied().unwrap_or(0.0);
+                (k, m >= sharp_angle && m <= blade_angle)
+            })
+            .collect();
+        let kept = |e: &&crate::surface::FeatureEdge| {
+            let (a, b) = (face_of(e.triangles[0]), face_of(e.triangles[1]));
+            if a != b { pair_kept.get(&(a.min(b), a.max(b))).copied().unwrap_or(false) } else { e.angle >= sharp_angle && e.angle <= blade_angle }
+        };
+        let mut edge_set: BTreeMap<(u32, u32), [u32; 2]> = BTreeMap::new();
+        for e in surface.feature_edges().iter().filter(kept) {
+            edge_set.insert((e.vertices[0], e.vertices[1]), e.triangles);
+        }
+        let mut adj: BTreeMap<u32, Vec<u32>> = BTreeMap::new();
+        for &(a, b) in edge_set.keys() {
+            adj.entry(a).or_default().push(b);
+            adj.entry(b).or_default().push(a);
+        }
+        let mut corners: Vec<u32> = adj
+            .iter()
+            .filter(|(v, n)| {
+                if n.len() != 2 {
+                    return true;
+                }
+                let (p, a, b) = (xyz(**v), xyz(n[0]), xyz(n[1]));
+                let (u, w) = (sub(p, a), sub(b, p));
+                crate::vec3::angle_deg(u, w) > CORNER_TURN
+            })
+            .map(|(v, _)| *v)
+            .collect();
+        corners.sort_unstable();
+
+        // Patches: flood fill across non-feature edges.
+        let neighbors = surface.neighbors();
+        let mut tri_patch = vec![u32::MAX; tris.len()];
+        let mut patches: Vec<Patch> = Vec::new();
+        for seed in 0..tris.len() {
+            if tri_patch[seed] != u32::MAX {
+                continue;
+            }
+            let id = patches.len() as u32;
+            let mut patch = Patch { triangles: Vec::new(), faces: Vec::new(), area: 0.0, width: f64::INFINITY };
+            tri_patch[seed] = id;
+            let mut stack = vec![seed as u32];
+            while let Some(t) = stack.pop() {
+                patch.triangles.push(t);
+                let tri = tris.get(t as usize).copied().unwrap_or([0; 3]);
+                let n = neighbors.get(t as usize).copied().unwrap_or([u32::MAX; 3]);
+                for k in 0..3 {
+                    let (a, b) = (tri[k], tri[(k + 1) % 3]);
+                    if edge_set.contains_key(&(a.min(b), a.max(b))) {
+                        continue;
+                    }
+                    if let Some(slot) = tri_patch.get_mut(n[k] as usize)
+                        && *slot == u32::MAX
+                    {
+                        *slot = id;
+                        stack.push(n[k]);
+                    }
+                }
+            }
+            patch.triangles.sort_unstable();
+            for &t in &patch.triangles {
+                let tri = tris.get(t as usize).copied().unwrap_or([0; 3]);
+                patch.area += 0.5 * norm(cross(sub(xyz(tri[1]), xyz(tri[0])), sub(xyz(tri[2]), xyz(tri[0]))));
+                patch.faces.push(surface.tri_face().get(t as usize).copied().unwrap_or(0));
+            }
+            patch.faces.sort_unstable();
+            patch.faces.dedup();
+            let mut perimeter = 0.0;
+            for &t in &patch.triangles {
+                let tri = tris.get(t as usize).copied().unwrap_or([0; 3]);
+                for k in 0..3 {
+                    let (a, b) = (tri[k], tri[(k + 1) % 3]);
+                    if edge_set.contains_key(&(a.min(b), a.max(b))) {
+                        perimeter += norm(sub(xyz(b), xyz(a)));
+                    }
+                }
+            }
+            if perimeter > 0.0 {
+                patch.width = 2.0 * patch.area / perimeter;
+            }
+            patches.push(patch);
+        }
+
+        // Curves: chains from corners, then closed loops.
+        let mut used: BTreeMap<(u32, u32), bool> = edge_set.keys().map(|k| (*k, false)).collect();
+        let is_corner = |v: u32| corners.binary_search(&v).is_ok();
+        let mut curves: Vec<Curve> = Vec::new();
+        let trace = |start: u32, next: u32, used: &mut BTreeMap<(u32, u32), bool>, curves: &mut Vec<Curve>| {
+            let mut verts = vec![start];
+            let (mut prev, mut cur) = (start, next);
+            loop {
+                let key = (prev.min(cur), prev.max(cur));
+                if let Some(u) = used.get_mut(&key) {
+                    *u = true;
+                }
+                verts.push(cur);
+                if cur == start || is_corner(cur) || verts.len() > edge_set.len() + 2 {
+                    break;
+                }
+                let nxt = adj.get(&cur).and_then(|n| n.iter().copied().find(|&w| w != prev && used.get(&(cur.min(w), cur.max(w))) == Some(&false)));
+                match nxt {
+                    Some(w) => {
+                        prev = cur;
+                        cur = w;
+                    }
+                    None => break,
+                }
+            }
+            let closed = verts.first() == verts.last() && verts.len() > 2;
+            let mut ps: Vec<u32> = Vec::new();
+            for w in verts.windows(2) {
+                if let Some(t) = edge_set.get(&(w[0].min(w[1]), w[0].max(w[1]))) {
+                    ps.extend(t.iter().filter_map(|&t| tri_patch.get(t as usize).copied()));
+                }
+            }
+            ps.sort_unstable();
+            ps.dedup();
+            curves.push(Curve { vertices: verts, closed, patches: ps });
+        };
+        for &c in &corners {
+            let ns = adj.get(&c).cloned().unwrap_or_default();
+            for n in ns {
+                if used.get(&(c.min(n), c.max(n))) == Some(&false) {
+                    trace(c, n, &mut used, &mut curves);
+                }
+            }
+        }
+        let keys: Vec<(u32, u32)> = used.keys().copied().collect();
+        for k in keys {
+            if used.get(&k) == Some(&false) {
+                trace(k.0, k.1, &mut used, &mut curves);
+            }
+        }
+        Ok(Geometry { surface, tree, boxes, normals, tri_patch, patches, curves, corners, bounds, far })
+    }
+
+    pub fn xyz(&self, v: u32) -> P {
+        self.surface.points().get(v as usize).copied().unwrap_or([0.0; 3])
+    }
+
+    /// The centroid of triangle `t`.
+    pub fn triangle_centre(&self, t: u32) -> P {
+        let [a, b, c] = self.triangle(t);
+        [0, 1, 2].map(|k| (a[k] + b[k] + c[k]) / 3.0)
+    }
+
+    pub fn triangle(&self, t: u32) -> [P; 3] {
+        let tri = self.surface.triangles().get(t as usize).copied().unwrap_or([0; 3]);
+        tri.map(|v| self.xyz(v))
+    }
+
+    /// The closest boundary point to `p` among triangles that pass `keep`: (triangle, point,
+    /// squared distance).
+    pub fn nearest(&self, p: P, keep: impl Fn(u32) -> bool) -> Option<(u32, P, f64)> {
+        let mut best: Option<(u32, P, f64)> = None;
+        let mut stack = vec![0usize];
+        while let Some(n) = stack.pop() {
+            let Some(node) = self.tree.nodes.get(n) else { continue };
+            let bd = node.bounds.dist2(p);
+            if best.is_some_and(|b| bd >= b.2) {
+                continue;
+            }
+            match node.children {
+                Some(l) => {
+                    let (dl, dr) = (
+                        self.tree.nodes.get(l).map_or(f64::INFINITY, |x| x.bounds.dist2(p)),
+                        self.tree.nodes.get(l + 1).map_or(f64::INFINITY, |x| x.bounds.dist2(p)),
+                    );
+                    if dl <= dr {
+                        stack.extend([l + 1, l]);
+                    } else {
+                        stack.extend([l, l + 1]);
+                    }
+                }
+                None => {
+                    for &t in self.tree.leaf(node) {
+                        if !keep(t) {
+                            continue;
+                        }
+                        let [a, b, c] = self.triangle(t);
+                        let q = closest_on_triangle(p, a, b, c);
+                        let d = dist2(p, q);
+                        if best.is_none_or(|x| d < x.2 || (d == x.2 && t < x.0)) {
+                            best = Some((t, q, d));
+                        }
+                    }
+                }
+            }
+        }
+        best
+    }
+
+    /// The triangles whose boxes overlap `b`, ascending.
+    pub fn overlapping(&self, b: &Aabb) -> Vec<u32> {
+        let mut out = Vec::new();
+        let mut stack = vec![0usize];
+        while let Some(n) = stack.pop() {
+            let Some(node) = self.tree.nodes.get(n) else { continue };
+            if !node.bounds.overlaps(b) {
+                continue;
+            }
+            match node.children {
+                Some(l) => stack.extend([l, l + 1]),
+                None => out.extend(self.tree.leaf(node).iter().copied().filter(|&t| self.boxes.get(t as usize).is_some_and(|x| x.overlaps(b)))),
+            }
+        }
+        out.sort_unstable();
+        out
+    }
+
+    /// The first triangle hit by the ray from `p` along `d` (unit) within `max`, with its
+    /// distance (floating point; used to measure, never to decide).
+    pub fn ray(&self, p: P, d: P, max: f64, skip: u32) -> Option<(u32, f64)> {
+        let q = [0, 1, 2].map(|k| p[k] + d[k] * max);
+        let mut best: Option<(u32, f64)> = None;
+        let mut stack = vec![0usize];
+        while let Some(n) = stack.pop() {
+            let Some(node) = self.tree.nodes.get(n) else { continue };
+            if !node.bounds.meets_segment(p, q) {
+                continue;
+            }
+            match node.children {
+                Some(l) => stack.extend([l, l + 1]),
+                None => {
+                    for &t in self.tree.leaf(node) {
+                        if t == skip {
+                            continue;
+                        }
+                        let [a, b, c] = self.triangle(t);
+                        if let Some(s) = ray_triangle(p, d, a, b, c)
+                            && s > 0.0
+                            && s <= max
+                            && best.is_none_or(|x| s < x.1)
+                        {
+                            best = Some((t, s));
+                        }
+                    }
+                }
+            }
+        }
+        best
+    }
+
+    /// Whether `p` is strictly inside the solid, decided exactly: the parity of the boundary
+    /// crossings of a segment from `p` to a far point, tried along fixed directions until one
+    /// meets no edge or vertex and does not start on the boundary. `None` when every direction is
+    /// degenerate (p lies on the boundary) or `p` is not admissible.
+    pub fn inside(&self, p: P) -> Option<bool> {
+        const DIRS: [P; 7] = [
+            [0.5773, 0.6123, 0.5401],
+            [-0.6301, 0.4372, 0.6417],
+            [0.3989, -0.7068, 0.5839],
+            [0.7211, 0.4109, -0.5577],
+            [-0.4402, -0.6005, -0.6675],
+            [0.1351, 0.9013, -0.4116],
+            [-0.8803, 0.1237, 0.4583],
+        ];
+        let pp = Point3::new(snap_to_domain(p[0]), snap_to_domain(p[1]), snap_to_domain(p[2])).ok()?;
+        'dirs: for d in DIRS {
+            let q = [0, 1, 2].map(|k| snap_to_domain(p[k] + d[k] * self.far));
+            let Ok(qq) = Point3::new(q[0], q[1], q[2]) else { continue };
+            let mut parity = false;
+            let mut stack = vec![0usize];
+            while let Some(n) = stack.pop() {
+                let Some(node) = self.tree.nodes.get(n) else { continue };
+                if !node.bounds.meets_segment(p, q) {
+                    continue;
+                }
+                match node.children {
+                    Some(l) => stack.extend([l, l + 1]),
+                    None => {
+                        for &t in self.tree.leaf(node) {
+                            match self.crossing(&pp, &qq, t) {
+                                Some(true) => parity = !parity,
+                                Some(false) => {}
+                                None => continue 'dirs,
+                            }
+                        }
+                    }
+                }
+            }
+            return Some(parity);
+        }
+        None
+    }
+
+    /// Whether segment `pq` crosses triangle `t` properly (`Some(true)`), misses it
+    /// (`Some(false)`), or touches it degenerately (`None`): exact.
+    fn crossing(&self, p: &Point3, q: &Point3, t: u32) -> Option<bool> {
+        let tri = self.surface.triangles().get(t as usize)?;
+        let a = self.surface.point3(tri[0])?;
+        let b = self.surface.point3(tri[1])?;
+        let c = self.surface.point3(tri[2])?;
+        let (sp, sq) = (orient3d(&a, &b, &c, p), orient3d(&a, &b, &c, q));
+        if sq.is_zero() {
+            // The far end is outside the boundary's box, so never on a triangle; a segment in
+            // the plane is ambiguous.
+            return None;
+        }
+        if sp.is_zero() {
+            // p in the triangle's plane: the segment leaves the plane at p, so it touches the
+            // triangle only if p lies in it.
+            let axis = [Axis::Z, Axis::X, Axis::Y].into_iter().find(|&ax| !orient2d_projected(&a, &b, &c, ax).is_zero())?;
+            let o = orient2d_projected(&a, &b, &c, axis);
+            let inside = [(&a, &b), (&b, &c), (&c, &a)].into_iter().all(|(u, v)| orient2d_projected(u, v, p, axis) != -o);
+            return if inside { None } else { Some(false) };
+        }
+        if sp == sq {
+            return Some(false);
+        }
+        let s = [orient3d(p, q, &a, &b), orient3d(p, q, &b, &c), orient3d(p, q, &c, &a)];
+        if s.iter().any(|x| x.is_zero()) {
+            let nz: Vec<Sign> = s.iter().copied().filter(|x| !x.is_zero()).collect();
+            return if nz.windows(2).all(|w| w[0] == w[1]) { None } else { Some(false) };
+        }
+        Some(s[0] == s[1] && s[1] == s[2])
+    }
+
+    /// The closest point of curve `c` to `p`: (point, segment index, parameter on it).
+    pub fn nearest_on_curve(&self, c: u32, p: P) -> Option<(P, usize, f64)> {
+        let curve = self.curves.get(c as usize)?;
+        let mut best: Option<(P, usize, f64, f64)> = None;
+        for (i, w) in curve.vertices.windows(2).enumerate() {
+            let (q, t) = closest_on_segment(p, self.xyz(w[0]), self.xyz(w[1]));
+            let d = dist2(p, q);
+            if best.is_none_or(|b| d < b.3) {
+                best = Some((q, i, t, d));
+            }
+        }
+        best.map(|(q, i, t, _)| (q, i, t))
+    }
+}
+
+/// Möller–Trumbore ray–triangle distance (floating point).
+fn ray_triangle(o: P, d: P, a: P, b: P, c: P) -> Option<f64> {
+    let (e1, e2) = (sub(b, a), sub(c, a));
+    let h = cross(d, e2);
+    let det = dot(e1, h);
+    if det.abs() < 1e-300 {
+        return None;
+    }
+    let inv = 1.0 / det;
+    let s = sub(o, a);
+    let u = inv * dot(s, h);
+    if !(-1e-12..=1.0 + 1e-12).contains(&u) {
+        return None;
+    }
+    let q = cross(s, e1);
+    let v = inv * dot(d, q);
+    if v < -1e-12 || u + v > 1.0 + 1e-12 {
+        return None;
+    }
+    Some(inv * dot(e2, q))
+}
+
+#[cfg(test)]
+pub(crate) mod tests {
+    use super::*;
+    use crate::surface::SurfaceOptions;
+
+    /// A closed axis-aligned box surface, one B-rep face per side.
+    pub(crate) fn box_surface(lo: P, hi: P) -> BoundarySurface {
+        let c = |i: usize| [if i & 1 == 0 { lo[0] } else { hi[0] }, if i & 2 == 0 { lo[1] } else { hi[1] }, if i & 4 == 0 { lo[2] } else { hi[2] }];
+        let pts: Vec<P> = (0..8).map(c).collect();
+        let quads = [[0, 2, 3, 1], [4, 5, 7, 6], [0, 1, 5, 4], [2, 6, 7, 3], [0, 4, 6, 2], [1, 3, 7, 5]];
+        let mut tris = Vec::new();
+        let mut face = Vec::new();
+        for (f, q) in quads.iter().enumerate() {
+            tris.push([q[0], q[1], q[2]]);
+            tris.push([q[0], q[2], q[3]]);
+            face.extend([f as u32; 2]);
+        }
+        BoundarySurface::new(&pts, &tris, &face, &SurfaceOptions::default()).unwrap()
+    }
+
+    #[test]
+    fn box_feature_graph_inside_and_nearest() {
+        let s = box_surface([0.0; 3], [2.0, 3.0, 4.0]);
+        let g = Geometry::new(&s, 30.0, 160.0).unwrap();
+        assert_eq!(g.corners.len(), 8);
+        assert_eq!(g.curves.len(), 12);
+        assert_eq!(g.patches.len(), 6);
+        assert!(g.curves.iter().all(|c| c.vertices.len() == 2 && c.patches.len() == 2 && !c.closed));
+        assert_eq!(g.inside([1.0, 1.5, 2.0]), Some(true));
+        assert_eq!(g.inside([3.0, 1.5, 2.0]), Some(false));
+        assert_eq!(g.inside([1.0, 1.5, 0.0]), None, "on the boundary");
+        assert_eq!(g.inside([0.0, 0.0, 0.0]), None, "a corner");
+        let (_, q, d) = g.nearest([1.0, 1.0, -1.0], |_| true).unwrap();
+        assert!((d - 1.0).abs() < 1e-12 && q[2] == 0.0);
+        let top = g.tri_patch[2];
+        let (_, q, _) = g.nearest([1.0, 1.0, -1.0], |t| g.tri_patch[t as usize] == top).unwrap();
+        assert_eq!(q[2], 4.0, "restricted to the top patch");
+        let (t, s) = g.ray([1.0, 1.5, 2.0], [0.0, 0.0, 1.0], 10.0, u32::MAX).unwrap();
+        assert!((s - 2.0).abs() < 1e-12 && g.tri_patch[t as usize] == top);
+    }
+}

crates/mesh/src/octree/grid.rs

@@ -0,0 +1,376 @@
+//! The adaptive grid: a box of root cubes, each refined as a 27-tree (a leaf splits into
+//! 3 × 3 × 3), strongly balanced and closed under the template patterns, then cut into a
+//! conforming hexahedral mesh by [`super::templates`].
+//!
+//! # Balance and closure
+//! *Strong balance*: leaves that share a vertex differ by at most one level, so a leaf's finer
+//! neighbours are exactly one level finer and a leaf next to a coarser one never has a finer one
+//! on that side. A vertex of a leaf at level `L` is *marked* when a finer leaf touches it, that
+//! is, when one of the eight level-`L` cells around the vertex is split. *Closure*: a leaf whose
+//! marks fit no template (two diagonal corners of a face, three of a face, …) is split too.
+//! Balancing and closing repeat until neither splits anything; both only split, so this ends,
+//! and every round runs over the leaves in sorted order, so the result is the same on every run.
+//!
+//! # Coordinates
+//! Everything is integer. A cell at level `L` has index `i` (per axis, in units of its own edge)
+//! and spans `6 · 3^(D − L)` lattice units, `D` the deepest level allowed; template vertices are
+//! in sixths of a leaf, so every vertex is a lattice point and vertices shared by neighbours are
+//! found by exact lookup.
+
+use std::collections::HashMap;
+
+use super::templates::{self, Pattern, Sym};
+use crate::{Error, Result};
+
+/// A cell: its level and index (in units of its edge at that level).
+pub(crate) type Key = (u8, [i64; 3]);
+
+/// The deepest level a tree may have (`3^12` lattice steps per root cell edge, × 6).
+pub(crate) const MAX_LEVEL: u8 = 12;
+
+fn pow3(n: u8) -> i64 {
+    3i64.pow(u32::from(n))
+}
+
+/// What covers a cell-sized region at some level.
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+enum Cover {
+    /// A leaf (at that level or coarser).
+    Leaf(Key),
+    /// The region is split at that level (finer leaves inside).
+    Split,
+    /// Outside the root box.
+    Outside,
+}
+
+/// The tree. `nodes` maps every cell that exists to whether it is a leaf.
+#[derive(Debug, Clone)]
+pub(crate) struct Tree {
+    /// Root cubes per axis.
+    pub roots: [i64; 3],
+    /// The deepest level allowed.
+    pub depth: u8,
+    nodes: HashMap<Key, bool>,
+    leaves: usize,
+    max_leaves: usize,
+}
+
+/// The sort order of [`Tree::sorted_leaves`]: level, then z, y, x.
+fn order(k: &Key) -> (u8, i64, i64, i64) {
+    (k.0, k.1[2], k.1[1], k.1[0])
+}
+
+/// The 26 neighbour offsets, in a fixed order.
+fn offsets() -> impl Iterator<Item = [i64; 3]> {
+    (0..27).map(|k| [k % 3 - 1, (k / 3) % 3 - 1, k / 9 - 1]).filter(|d| *d != [0, 0, 0])
+}
+
+impl Tree {
+    /// One leaf per root cube.
+    pub fn new(roots: [i64; 3], depth: u8, max_leaves: usize) -> Result<Tree> {
+        if roots.iter().any(|&n| n < 1) || depth > MAX_LEVEL {
+            return Err(Error::Invalid(format!("grid of {roots:?} root cells to depth {depth}")));
+        }
+        let count = roots.iter().try_fold(1usize, |a, &n| usize::try_from(n).ok().and_then(|n| a.checked_mul(n)));
+        let count = count.filter(|&c| c <= max_leaves).ok_or_else(|| Error::Limit(format!("{roots:?} root cells exceed {max_leaves} cells")))?;
+        let mut nodes = HashMap::with_capacity(count);
+        for z in 0..roots[2] {
+            for y in 0..roots[1] {
+                for x in 0..roots[0] {
+                    nodes.insert((0, [x, y, z]), true);
+                }
+            }
+        }
+        Ok(Tree { roots, depth, nodes, leaves: count, max_leaves })
+    }
+
+    pub fn leaf_count(&self) -> usize {
+        self.leaves
+    }
+
+    /// Whether `key` is a leaf.
+    pub fn is_leaf(&self, key: &Key) -> bool {
+        self.nodes.get(key).copied().unwrap_or(false)
+    }
+
+    fn inside(&self, (level, i): Key) -> bool {
+        let n = pow3(level);
+        (0..3).all(|k| i[k] >= 0 && i[k] < self.roots[k] * n)
+    }
+
+    /// What covers the level-`L` cell `key`.
+    fn cover(&self, key: Key) -> Cover {
+        if !self.inside(key) {
+            return Cover::Outside;
+        }
+        let (level, i) = key;
+        for l in (0..=level).rev() {
+            let s = pow3(level - l);
+            let k = (l, [i[0].div_euclid(s), i[1].div_euclid(s), i[2].div_euclid(s)]);
+            match self.nodes.get(&k) {
+                Some(true) => return Cover::Leaf(k),
+                Some(false) => return Cover::Split,
+                None => {}
+            }
+        }
+        Cover::Outside
+    }
+
+    /// Split a leaf into its 27 children.
+    pub fn split(&mut self, key: Key) -> Result<()> {
+        let level = key.0;
+        if !self.is_leaf(&key) {
+            return Ok(());
+        }
+        if level >= self.depth {
+            return Err(Error::Limit(format!("refinement needs more than {} levels", self.depth)));
+        }
+        if self.leaves + 26 > self.max_leaves {
+            return Err(Error::Limit(format!("the grid needs more than {} cells", self.max_leaves)));
+        }
+        self.nodes.insert(key, false);
+        for child in Self::children(key) {
+            self.nodes.insert(child, true);
+        }
+        self.leaves += 26;
+        Ok(())
+    }
+
+    /// The leaves in sorted order (level, then z, y, x).
+    pub fn sorted_leaves(&self) -> Vec<Key> {
+        let mut v: Vec<Key> = self.nodes.iter().filter(|(_, leaf)| **leaf).map(|(k, _)| *k).collect();
+        v.sort_by_key(order);
+        v
+    }
+
+    /// The children of a split cell.
+    pub fn children(key: Key) -> impl Iterator<Item = Key> {
+        let (level, i) = key;
+        (0..27).map(move |k| (level + 1, [3 * i[0] + k % 3, 3 * i[1] + (k / 3) % 3, 3 * i[2] + k / 9]))
+    }
+
+    /// Split until leaves sharing a vertex differ by at most one level, checking only the leaves
+    /// `work` (in sorted order) and the children of the cells it splits: a split only makes
+    /// neighbours finer, so only new leaves can be too fine for a neighbour, and the result is
+    /// the unique coarsest balanced refinement whatever the order. Every split cell goes to `log`.
+    fn balance_from(&mut self, mut work: Vec<Key>, log: &mut Vec<Key>) -> Result<()> {
+        work.reverse();
+        while let Some(key) = work.pop() {
+            if !self.is_leaf(&key) || key.0 < 2 {
+                continue;
+            }
+            let (level, i) = key;
+            for d in offsets() {
+                if let Cover::Leaf(k) = self.cover((level, [i[0] + d[0], i[1] + d[1], i[2] + d[2]]))
+                    && k.0 + 1 < level
+                {
+                    self.split(k)?;
+                    log.push(k);
+                    let mut kids: Vec<Key> = Self::children(k).collect();
+                    kids.reverse();
+                    work.extend(kids);
+                    work.push(key);
+                    break;
+                }
+            }
+        }
+        Ok(())
+    }
+
+    /// The leaves whose marks a split of `key` may change: its children and the leaves beside it
+    /// on its level.
+    fn touched_by(&self, key: Key, out: &mut Vec<Key>) {
+        out.extend(Self::children(key));
+        let (level, i) = key;
+        out.extend(offsets().map(|d| (level, [i[0] + d[0], i[1] + d[1], i[2] + d[2]])).filter(|k| self.is_leaf(k)));
+    }
+
+    /// The marked corners of a leaf (bit `c` for corner `c`: bit 0 x, bit 1 y, bit 2 z).
+    pub fn marks(&self, (level, i): Key) -> u8 {
+        let mut mask = 0u8;
+        for c in 0..8usize {
+            let o = [(c & 1) as i64, ((c >> 1) & 1) as i64, ((c >> 2) & 1) as i64];
+            let split = (0..8usize).any(|e| {
+                let n = [i[0] + o[0] - 1 + (e & 1) as i64, i[1] + o[1] - 1 + ((e >> 1) & 1) as i64, i[2] + o[2] - 1 + ((e >> 2) & 1) as i64];
+                self.nodes.get(&(level, n)) == Some(&false)
+            });
+            if split {
+                mask |= 1 << c;
+            }
+        }
+        mask
+    }
+
+    /// Balance, then split every leaf whose marks fit no template (and every leaf marked at all
+    /// eight corners, whose template is the split), until nothing changes. `max_rounds` bounds the
+    /// rounds. Each round re-checks only the leaves the previous splits can have changed, so a
+    /// round costs what it splits, not the size of the tree.
+    pub fn close(&mut self, max_rounds: usize, skip: &dyn Fn(Key) -> bool) -> Result<()> {
+        let mut fresh = self.sorted_leaves();
+        let mut touched = fresh.clone();
+        for _ in 0..max_rounds {
+            let mut log = Vec::new();
+            self.balance_from(std::mem::take(&mut fresh), &mut log)?;
+            for &k in &log {
+                self.touched_by(k, &mut touched);
+            }
+            touched.retain(|k| self.is_leaf(k));
+            touched.sort_by_key(order);
+            touched.dedup();
+            let bad: Vec<Key> = touched
+                .drain(..)
+                .filter(|&k| !skip(k))
+                .filter(|&k| !matches!(templates::classify(self.marks(k)), Some(Pattern::None | Pattern::Corner | Pattern::Edge | Pattern::Face)))
+                .collect();
+            if bad.is_empty() {
+                return Ok(());
+            }
+            for &k in &bad {
+                self.split(k)?;
+            }
+            for &k in &bad {
+                fresh.extend(Self::children(k));
+                self.touched_by(k, &mut touched);
+            }
+            fresh.sort_by_key(order);
+        }
+        Err(Error::Limit(format!("grid closure did not settle in {max_rounds} rounds")))
+    }
+
+    /// Lattice units per edge of a level-`L` cell.
+    pub fn span(&self, level: u8) -> i64 {
+        6 * pow3(self.depth - level)
+    }
+
+    /// Cut every leaf by its template into a conforming mesh on the integer lattice. Fails if a
+    /// leaf's marks fit no template (call [`Tree::close`] first).
+    pub fn mesh(&self, skip: &dyn Fn(Key) -> bool) -> Result<LatticeMesh> {
+        let syms: Vec<Sym> = templates::symmetries();
+        let canon: Vec<(Pattern, Vec<templates::Hex>)> =
+            [Pattern::None, Pattern::Corner, Pattern::Edge, Pattern::Face].into_iter().map(|p| (p, templates::template(p))).collect();
+        let mut out = LatticeMesh::default();
+        let mut index: HashMap<[i64; 3], u32> = HashMap::new();
+        for key in self.sorted_leaves() {
+            if skip(key) {
+                continue;
+            }
+            let mask = self.marks(key);
+            let pattern = templates::classify(mask).filter(|p| *p != Pattern::Full);
+            let (pattern, sym) = pattern
+                .and_then(|p| templates::orientation(&syms, p, mask).map(|s| (p, s)))
+                .ok_or_else(|| Error::Failed(format!("grid cell {key:?} has marks {mask:#010b} that fit no template")))?;
+            let hexes = canon.iter().find(|(p, _)| *p == pattern).map(|(_, h)| h.as_slice()).unwrap_or(&[]);
+            let (level, i) = key;
+            let unit = pow3(self.depth - level);
+            for h in hexes {
+                let mut g = h.map(|p| {
+                    let q = sym.apply(p);
+                    [0, 1, 2].map(|k| (6 * i[k] + q[k]) * unit)
+                });
+                if templates::corner0_orientation(&g) < 0 {
+                    g = templates::reversed(g);
+                }
+                let mut ids = [0u32; 8];
+                for (id, p) in ids.iter_mut().zip(g) {
+                    *id = match index.get(&p) {
+                        Some(&v) => v,
+                        None => {
+                            let v = u32::try_from(out.points.len()).map_err(|_| Error::Limit("too many grid points".into()))?;
+                            index.insert(p, v);
+                            out.points.push(p);
+                            v
+                        }
+                    };
+                }
+                out.hexes.push(ids);
+                out.level.push(level);
+            }
+        }
+        Ok(out)
+    }
+}
+
+/// A conforming hexahedral mesh on the integer lattice.
+#[derive(Debug, Clone, Default)]
+pub(crate) struct LatticeMesh {
+    pub points: Vec<[i64; 3]>,
+    pub hexes: Vec<[u32; 8]>,
+    /// The tree level of the leaf each hexahedron came from.
+    pub level: Vec<u8>,
+}
+
+#[cfg(test)]
+mod tests {
+    use std::collections::HashMap;
+
+    use super::*;
+    use crate::quality::{Validity, element_validity};
+    use fieldcraft_model::ElementKind;
+
+    /// A tiny deterministic generator for test trees.
+    struct Lcg(u64);
+    impl Lcg {
+        fn next(&mut self) -> u64 {
+            self.0 = self.0.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407);
+            self.0 >> 33
+        }
+    }
+
+    fn check_conforming(tree: &Tree) {
+        let m = tree.mesh(&|_| false).unwrap();
+        let mut faces: HashMap<[u32; 4], usize> = HashMap::new();
+        let mut volume = 0.0;
+        for h in &m.hexes {
+            let xyz: Vec<[f64; 3]> = h.iter().map(|&v| m.points[v as usize].map(|c| c as f64)).collect();
+            let c = element_validity(ElementKind::Hex8, &xyz).unwrap();
+            assert_eq!(c.validity, Validity::Valid);
+            volume += c.volume;
+            for f in ElementKind::Hex8.faces() {
+                let mut q = [h[f[0]], h[f[1]], h[f[2]], h[f[3]]];
+                q.sort();
+                *faces.entry(q).or_default() += 1;
+            }
+        }
+        let side = tree.span(0) as f64;
+        let expected = side.powi(3) * (tree.roots[0] * tree.roots[1] * tree.roots[2]) as f64;
+        assert!((volume - expected).abs() < 1e-6 * expected, "volume {volume} vs {expected}");
+        let hi = tree.roots.map(|n| n * tree.span(0));
+        for (q, n) in faces {
+            let on_box = (0..3).any(|k| q.iter().all(|&v| m.points[v as usize][k] == 0) || q.iter().all(|&v| m.points[v as usize][k] == hi[k]));
+            assert_eq!(n, if on_box { 1 } else { 2 }, "face {q:?}");
+        }
+    }
+
+    #[test]
+    fn random_trees_close_into_conforming_valid_meshes() {
+        let mut rng = Lcg(7);
+        for case in 0..12 {
+            let roots = [1 + (rng.next() % 2) as i64, 1 + (rng.next() % 2) as i64, 1];
+            let mut tree = Tree::new(roots, 3, 1_000_000).unwrap();
+            for _ in 0..(1 + rng.next() % 6) {
+                let leaves = tree.sorted_leaves();
+                let k = leaves[(rng.next() as usize) % leaves.len()];
+                if k.0 < 3 {
+                    tree.split(k).unwrap();
+                }
+            }
+            tree.close(100, &|_| false).unwrap();
+            for k in tree.sorted_leaves() {
+                assert!(templates::classify(tree.marks(k)).is_some(), "case {case}");
+            }
+            check_conforming(&tree);
+        }
+    }
+
+    #[test]
+    fn caps_are_errors() {
+        assert!(Tree::new([0, 1, 1], 2, 10).is_err());
+        assert!(Tree::new([10, 10, 10], 2, 999).is_err());
+        let mut t = Tree::new([1, 1, 1], 1, 100).unwrap();
+        t.split((0, [0, 0, 0])).unwrap();
+        assert!(t.split((1, [0, 0, 0])).is_err(), "deeper than the depth");
+        let mut t = Tree::new([1, 1, 1], 3, 30).unwrap();
+        t.split((0, [0, 0, 0])).unwrap();
+        assert!(t.split((1, [0, 0, 0])).is_err(), "more cells than allowed");
+    }
+}

crates/mesh/src/octree/map.rs

@@ -0,0 +1,282 @@
+//! The topological map from the boundary's feature graph (corners, curves, patches) onto the
+//! core's boundary quads, after X. Gao, H. Shen, D. Panozzo, *Feature preserving octree-based
+//! hexahedral meshing*, CGF 38(5), 2019, §3.3 (paper only; the authors' code has no licence and
+//! is not read). Gao et al. trace the curves first and read the patches off them; here the
+//! patches come first and the curves are their borders, which keeps every region on the side of
+//! the part it faces:
+//!
+//! 1. **Label.** Each quad takes the patch its outward normal meets first (within two quad
+//!    sizes), or else the patch nearest its centre.
+//! 2. **Clean.** A quad whose label differs from all its edge neighbours takes their majority
+//!    label; then every patch keeps only its largest connected group of quads, the other
+//!    groups joining the label they border most. Both repeat a bounded number of times.
+//! 3. **Check.** The map must be a bijection: every patch has a region, and two regions border
+//!    each other exactly where the part has a curve between their patches.
+//! 4. **Features.** Region borders are the curves (each border edge takes the curve between its
+//!    two patches nearest to it); a vertex where three or more regions meet is a corner and takes
+//!    the nearest corner of the part around those patches. Corners that are not region junctions
+//!    (where one curve turns sharply) go to the nearest vertex on their curve's border.
+//!
+//! Any failure names where the grid is too coarse for the features; the caller refines there and
+//! tries again (§3.3, "Bijective map").
+
+use std::collections::{BTreeMap, BTreeSet};
+
+use super::geom::{Geometry, closest_on_segment, dist2};
+use super::smooth::Constraint;
+use super::topo::quad_edges;
+use crate::vec3::{P, cross, norm, sub};
+
+/// Rounds of label cleaning.
+const CLEAN_ROUNDS: usize = 16;
+
+/// Where and why the map failed.
+#[derive(Debug, Clone)]
+pub(crate) struct MapFailure {
+    pub at: Vec<P>,
+    pub why: String,
+}
+
+/// The map: a constraint per quad vertex and a patch per quad.
+#[derive(Debug, Clone)]
+pub(crate) struct FeatureMap {
+    pub role: BTreeMap<u32, Constraint>,
+    pub quad_patch: Vec<u32>,
+}
+
+/// Map the feature graph onto `quads` (the core's boundary, outward) over `xyz`.
+pub(crate) fn map_features(g: &Geometry, quads: &[[u32; 4]], xyz: &dyn Fn(u32) -> P) -> Result<FeatureMap, MapFailure> {
+    let n = quads.len();
+    let centre = |q: usize| {
+        let mut c = [0.0; 3];
+        for &v in &quads[q] {
+            let p = xyz(v);
+            for k in 0..3 {
+                c[k] += 0.25 * p[k];
+            }
+        }
+        c
+    };
+    let edges = quad_edges(quads);
+    let mut nbr: Vec<Vec<u32>> = vec![Vec::new(); n];
+    for qs in edges.values() {
+        if let [a, b] = qs.as_slice() {
+            nbr[*a as usize].push(*b);
+            nbr[*b as usize].push(*a);
+        }
+    }
+
+    // 1. Label: the patch the quad faces (a ray along its outward normal), else the nearest.
+    let mut label: Vec<u32> = (0..n)
+        .map(|q| {
+            let c = centre(q);
+            let p = quads[q].map(xyz);
+            let nrm = cross(sub(p[2], p[0]), sub(p[3], p[1]));
+            let l = norm(nrm);
+            let reach = 2.0 * (norm(sub(p[2], p[0])) + norm(sub(p[3], p[1])));
+            let hit = if l > 0.0 { g.ray(c, nrm.map(|x| x / l), reach, u32::MAX).map(|(t, _)| t) } else { None };
+            hit.or_else(|| g.nearest(c, |_| true).map(|(t, _, _)| t)).and_then(|t| g.tri_patch.get(t as usize).copied()).unwrap_or(0)
+        })
+        .collect();
+
+    // A patch no quad faces (a narrow face the core's boundary passes by) takes the quad
+    // nearest to it, which the cleaning keeps.
+    let mut seeded: BTreeSet<usize> = BTreeSet::new();
+    let have: BTreeSet<u32> = label.iter().copied().collect();
+    for p in 0..g.patches.len() as u32 {
+        if have.contains(&p) {
+            continue;
+        }
+        let best = (0..n)
+            .filter(|q| !seeded.contains(q))
+            .filter_map(|q| g.nearest(centre(q), |t| g.tri_patch.get(t as usize) == Some(&p)).map(|x| (q, x.2)))
+            .min_by(|a, b| a.1.total_cmp(&b.1).then(a.0.cmp(&b.0)));
+        if let Some((q, _)) = best {
+            label[q] = p;
+            seeded.insert(q);
+        }
+    }
+    // 2. Clean.
+    for _ in 0..CLEAN_ROUNDS {
+        let mut changed = false;
+        for q in 0..n {
+            let mut count: BTreeMap<u32, usize> = BTreeMap::new();
+            for &m in &nbr[q] {
+                *count.entry(label[m as usize]).or_default() += 1;
+            }
+            let same = count.get(&label[q]).copied().unwrap_or(0);
+            if same == 0
+                && !nbr[q].is_empty()
+                && !seeded.contains(&q)
+                && let Some((&best, _)) = count.iter().filter(|(l, _)| **l != label[q]).max_by(|a, b| a.1.cmp(b.1).then(b.0.cmp(a.0)))
+            {
+                label[q] = best;
+                changed = true;
+            }
+        }
+        let mut comp = vec![usize::MAX; n];
+        let mut comps: Vec<(u32, Vec<usize>)> = Vec::new();
+        for s in 0..n {
+            if comp[s] != usize::MAX {
+                continue;
+            }
+            let id = comps.len();
+            let mut members = vec![s];
+            comp[s] = id;
+            let mut i = 0;
+            while i < members.len() {
+                let q = members[i];
+                i += 1;
+                for &m in &nbr[q] {
+                    let m = m as usize;
+                    if comp[m] == usize::MAX && label[m] == label[s] {
+                        comp[m] = id;
+                        members.push(m);
+                    }
+                }
+            }
+            comps.push((label[s], members));
+        }
+        let mut largest: BTreeMap<u32, usize> = BTreeMap::new();
+        for (i, (l, m)) in comps.iter().enumerate() {
+            let e = largest.entry(*l).or_insert(i);
+            if comps[*e].1.len() < m.len() {
+                *e = i;
+            }
+        }
+        for (i, (l, members)) in comps.iter().enumerate() {
+            if largest.get(l) == Some(&i) {
+                continue;
+            }
+            let mut count: BTreeMap<u32, usize> = BTreeMap::new();
+            for &q in members {
+                for &m in &nbr[q] {
+                    if comp[m as usize] != i {
+                        *count.entry(label[m as usize]).or_default() += 1;
+                    }
+                }
+            }
+            if let Some((&best, _)) = count.iter().max_by(|a, b| a.1.cmp(b.1).then(b.0.cmp(a.0))) {
+                for &q in members {
+                    label[q] = best;
+                }
+                changed = true;
+            }
+        }
+        if !changed {
+            break;
+        }
+    }
+
+    // 3. Check: every patch has a region; borders match curves.
+    let present: BTreeSet<u32> = label.iter().copied().collect();
+    for (p, patch) in g.patches.iter().enumerate() {
+        if !present.contains(&(p as u32)) {
+            let at: Vec<P> = patch.triangles.iter().take(64).map(|&t| g.triangle_centre(t)).collect();
+            return Err(MapFailure { at, why: format!("patch {p} (B-rep faces {:?}) has no region", patch.faces) });
+        }
+    }
+    let mut curves_between: BTreeMap<(u32, u32), Vec<u32>> = BTreeMap::new();
+    for (i, c) in g.curves.iter().enumerate() {
+        match c.patches.as_slice() {
+            [a, b] => curves_between.entry((*a, *b)).or_default().push(i as u32),
+            [a] => curves_between.entry((*a, *a)).or_default().push(i as u32),
+            _ => {}
+        }
+    }
+    let mut role: BTreeMap<u32, Constraint> = BTreeMap::new();
+    let mut border_pairs: BTreeSet<(u32, u32)> = BTreeSet::new();
+    let mut vertex_labels: BTreeMap<u32, BTreeSet<u32>> = BTreeMap::new();
+    for (q, quad) in quads.iter().enumerate() {
+        for &v in quad {
+            vertex_labels.entry(v).or_default().insert(label[q]);
+        }
+    }
+    for ((a, b), qs) in &edges {
+        let (q0, q1) = match qs.as_slice() {
+            [x, y] => (*x as usize, *y as usize),
+            _ => continue,
+        };
+        let (la, lb) = (label[q0].min(label[q1]), label[q0].max(label[q1]));
+        if la == lb {
+            continue;
+        }
+        let Some(cs) = curves_between.get(&(la, lb)) else {
+            return Err(MapFailure { at: vec![xyz(*a), xyz(*b)], why: format!("regions of patches {la} and {lb} meet where the part has no curve") });
+        };
+        border_pairs.insert((la, lb));
+        let mid = [0, 1, 2].map(|k| 0.5 * (xyz(*a)[k] + xyz(*b)[k]));
+        let c = nearest_curve(g, cs, mid);
+        for v in [*a, *b] {
+            role.entry(v).or_insert(Constraint::Curve(c));
+        }
+    }
+    for (pair, cs) in &curves_between {
+        if pair.0 != pair.1 && !border_pairs.contains(pair) {
+            let at: Vec<P> = cs.iter().flat_map(|&c| g.curves[c as usize].vertices.iter().map(|&v| g.xyz(v)).collect::<Vec<_>>()).collect();
+            return Err(MapFailure { at, why: format!("the regions of patches {} and {} do not meet along their curve", pair.0, pair.1) });
+        }
+    }
+
+    // 4. Corners.
+    let mut taken: BTreeSet<u32> = BTreeSet::new();
+    let junctions: Vec<u32> = vertex_labels.iter().filter(|(_, l)| l.len() >= 3).map(|(v, _)| *v).collect();
+    for &v in &junctions {
+        let p = xyz(v);
+        let labels = &vertex_labels[&v];
+        let best = g
+            .corners
+            .iter()
+            .copied()
+            .filter(|c| !taken.contains(c) && corner_patches(g, *c).is_superset(labels))
+            .min_by(|&a, &b| dist2(g.xyz(a), p).total_cmp(&dist2(g.xyz(b), p)).then(a.cmp(&b)));
+        match best {
+            Some(c) => {
+                taken.insert(c);
+                role.insert(v, Constraint::Corner(c));
+            }
+            None => return Err(MapFailure { at: vec![p], why: format!("regions {labels:?} meet where the part has no corner") }),
+        }
+    }
+    for &c in &g.corners {
+        if taken.contains(&c) {
+            continue;
+        }
+        let p = g.xyz(c);
+        let cp = corner_patches(g, c);
+        let best = role
+            .iter()
+            .filter(|(_, r)| matches!(r, Constraint::Curve(_)))
+            .map(|(v, _)| *v)
+            .filter(|v| vertex_labels.get(v).is_some_and(|l| l.is_subset(&cp)))
+            .min_by(|&a, &b| dist2(xyz(a), p).total_cmp(&dist2(xyz(b), p)).then(a.cmp(&b)));
+        match best {
+            Some(v) => {
+                role.insert(v, Constraint::Corner(c));
+            }
+            None => return Err(MapFailure { at: vec![p], why: format!("the feature corner at {p:?} has no vertex") }),
+        }
+    }
+    for (q, quad) in quads.iter().enumerate() {
+        for &v in quad {
+            role.entry(v).or_insert(Constraint::Patch(label[q]));
+        }
+    }
+    Ok(FeatureMap { role, quad_patch: label })
+}
+
+/// The patches around a corner of the part (on either side of the curves that end there).
+fn corner_patches(g: &Geometry, c: u32) -> BTreeSet<u32> {
+    g.curves.iter().filter(|cv| cv.vertices.first() == Some(&c) || cv.vertices.last() == Some(&c)).flat_map(|cv| cv.patches.iter().copied()).collect()
+}
+
+/// The curve of `cs` nearest to `p`.
+fn nearest_curve(g: &Geometry, cs: &[u32], p: P) -> u32 {
+    cs.iter().copied().min_by(|&a, &b| curve_distance(g, a, p).total_cmp(&curve_distance(g, b, p)).then(a.cmp(&b))).unwrap_or(0)
+}
+
+fn curve_distance(g: &Geometry, c: u32, p: P) -> f64 {
+    g.curves.get(c as usize).map_or(f64::INFINITY, |cv| {
+        cv.vertices.windows(2).map(|w| dist2(p, closest_on_segment(p, g.xyz(w[0]), g.xyz(w[1])).0)).fold(f64::INFINITY, f64::min)
+    })
+}

crates/mesh/src/octree/mod.rs

@@ -0,0 +1,801 @@
+//! Unstructured all-hexahedral meshing of any closed boundary: an adaptive grid cut into
+//! conforming hexahedra, a core of grid elements inside the part, two boundary layers that carry
+//! the part's faces, sharp curves and corners, and an optimisation that fits them to the
+//! geometry while keeping every element valid.
+//!
+//! # Pipeline
+//! 1. **Frame.** The grid is aligned with the part: the normals of its two largest planar
+//!    patches when they are orthogonal and cover more area than the coordinate axes do, the
+//!    coordinate axes otherwise (octree methods are orientation dependent; N. Pietroni et al.,
+//!    *Hex-mesh generation and processing: a survey*, ACM TOG 42(2), 2022, §4.6).
+//! 2. **Grid** ([`grid`]). Root cubes of the target size over the part's box, refined as a
+//!    27-tree where a B-rep face curves (its normals turn by more than
+//!    [`OctreeOptions::curvature_angle`] inside a cell), where the part is thin (fewer than
+//!    [`OctreeOptions::thickness_cells`] cells across, measured along the inward normal of each
+//!    boundary triangle) and where a patch is narrow, then strongly balanced and closed under the
+//!    transition templates.
+//! 3. **Conforming hexahedra** ([`templates`]): the 3-refinement node-marking templates
+//!    (R. Schneiders, 1996; Y. Ito, A. M. Shih, B. K. Soni, IJNME 77(13), 2009), so every grid
+//!    leaf becomes 1, 4, 11, 22 or 27 hexahedra and the grid is conforming by construction.
+//! 4. **Core.** Elements whose eight corners are inside the part (exact ray parity on the
+//!    boundary triangles) and at least [`OctreeOptions::core_offset`] of their cell size from it.
+//!    Elements are dropped until the core's boundary is a 2-manifold and every boundary vertex can
+//!    move outward (its quads' normals lie in one half-space: the buffer rule of H. Tong,
+//!    E. Halilaj, Y. J. Zhang, *HybridOctree_Hex*, arXiv 2401.05984, 2024); its components and
+//!    genus must match the part's shells.
+//! 5. **Feature map** ([`map`]): the part's patches, sharp curves and corners mapped onto the
+//!    core's boundary quads (X. Gao, H. Shen, D. Panozzo, *Feature preserving octree-based
+//!    hexahedral meshing*, CGF 38(5), 2019, §3.3), patches first.
+//! 6. **Padding and pillowing** ([`topo`]): one layer of hexahedra on the core's boundary (so no
+//!    element keeps two faces on the boundary), then each patch's layer is pillowed, which inserts
+//!    a ring of quads along every sharp curve on both sides: a curve then runs between quads that
+//!    each have one edge on it, so straight and curved features are captured without degenerate
+//!    quads (Gao et al. 2019, §3.4, Fig. 9).
+//! 7. **Fit** ([`smooth`]): the foldover-free energy of V. Garanzha et al., *Foldover-free maps
+//!    in 50 lines of code*, ACM TOG 40(4), 2021, plus a fitting term that draws every boundary
+//!    vertex onto its corner, curve or patch with a growing weight (Gao et al. 2019, §3.5); the
+//!    barrier keeps every element valid while the boundary moves. Each boundary vertex then goes
+//!    exactly onto its feature where its elements stay valid; where they would not (a blade, a
+//!    near-tangent contact) the mesh stays off by [`OctreeReport::boundary_distance`].
+//! 8. **Gate.** The worst elements are improved by compass search, then every element is proven
+//!    valid by the certified Jacobian check (Bernstein bounds, [`crate::quality`]; corner
+//!    Jacobians alone do not prove a trilinear hexahedron valid, J. Marschner et al., CGF 39(5),
+//!    2020); elements that are not are repaired locally first. A failure in 4–8 names where the
+//!    grid is too coarse; the leaves there are split once more and the mesh is rebuilt, at most
+//!    [`OctreeOptions::attempts`] times, and a rebuild may not grow the grid more than
+//!    [`RETRY_GROWTH`] times (Gao et al. 2019, §3.3).
+//!
+//! # Why 3-refinement templates
+//! The dual schemes of M. Livesu, L. Pitzalis, G. Cherchi, *Optimal dual schemes for adaptive
+//! grid based hexmeshing* (ACM TOG 41(2), 2022) give a simpler singular structure, but their
+//! installation depends on a paired grid whose transition vertices the reference pipeline
+//! (L. Pitzalis et al., *Generalized adaptive refinement for grid-based hexahedral meshing*, ACM
+//! TOG 40(6), 2021) selects with an integer linear program, and on heuristics that are retried
+//! with other parameters when they fail. The node-marking templates need only the balance and
+//! closure rules and are proven conforming by an exhaustive test. Their limit is the one
+//! Schneiders states: they close convex refined regions only, so a leaf next to a concave corner
+//! of a finer region is split, and the split spreads to the region's orthogonal hull (H. Tong,
+//! Y. J. Zhang, *MCHex*, arXiv 2511.02064, 2025, meet the same limit with the same templates).
+//! Curvature therefore refines only [`OctreeOptions::curvature_levels`] levels; deeper levels
+//! follow thickness and narrow faces, whose refined regions are mostly convex. Schneiders'
+//! modified edge rule (every edge with a marked end split in three) or the dual schemes are the
+//! upgrade path.
+//!
+//! # Determinism and bounds
+//! The same boundary and options give the same mesh: grid cells are processed in sorted order,
+//! maps are only looked up, the optimiser sweeps vertices in index order. Grid cells, elements,
+//! refinement depth, closure rounds, optimiser sweeps and attempts are capped.
+
+mod geom;
+mod grid;
+mod map;
+mod smooth;
+mod templates;
+mod topo;
+
+#[cfg(test)]
+mod tests;
+
+use std::cell::{Cell, RefCell};
+use std::collections::{BTreeMap, BTreeSet, HashMap, HashSet};
+
+use fieldcraft_model::ElementKind;
+
+pub use crate::hex::BoundaryFace;
+use crate::hex::face_flux;
+use crate::quality::{GateOptions, Metric, Validity, element_quality, validity_gate_checks};
+use crate::size::{SizeField, checked_size};
+use crate::surface::BoundarySurface;
+use crate::vec3::{P, cross, dot, norm, sub};
+use crate::{Error, Result};
+
+use geom::{Aabb, Geometry};
+use grid::Tree;
+use smooth::{Constraint, Optimizer};
+
+/// Options and caps of [`octree_hex_mesh`].
+#[derive(Debug, Clone, PartialEq)]
+pub struct OctreeOptions {
+    /// Refinement levels below the root cells (each one a third of the size), at most 8.
+    pub max_depth: u8,
+    /// The most grid cells and the most output elements.
+    pub max_elements: usize,
+    /// A cell is split where the normals of one B-rep face inside it turn by more than this
+    /// (degrees).
+    pub curvature_angle: f64,
+    /// Levels below the root that curvature may refine (see the module docs).
+    pub curvature_levels: u8,
+    /// A cell is split where the part is less than this many cells thick.
+    pub thickness_cells: f64,
+    /// Core elements keep this fraction of their cell size from the boundary.
+    pub core_offset: f64,
+    /// Optimiser sweeps over all vertices.
+    pub max_sweeps: usize,
+    /// Builds, counting the rebuilds after refining where a stage failed (1 = no retry).
+    pub attempts: usize,
+    /// Boundary edges whose faces meet at more than this angle (degrees between normals) are
+    /// sharp features the mesh captures with its own edges; smoother face boundaries are not.
+    pub sharp_angle: f64,
+    /// Edges whose normals turn by more than this (degrees: thin blades, near-tangent contacts)
+    /// are not captured; the mesh rounds them.
+    pub blade_angle: f64,
+}
+
+impl Default for OctreeOptions {
+    fn default() -> Self {
+        OctreeOptions {
+            max_depth: 2,
+            max_elements: 5_000_000,
+            curvature_angle: 25.0,
+            curvature_levels: 1,
+            thickness_cells: 3.0,
+            core_offset: 0.3,
+            max_sweeps: 300,
+            attempts: 4,
+            sharp_angle: 30.0,
+            blade_angle: 160.0,
+        }
+    }
+}
+
+/// What [`octree_hex_mesh`] did.
+#[derive(Debug, Clone, PartialEq)]
+pub struct OctreeReport {
+    /// Builds run (1 + rebuilds after a failure).
+    pub attempts: usize,
+    /// Grid leaves of the last build.
+    pub grid_cells: usize,
+    /// Core elements (from the grid) and boundary-layer elements.
+    pub core_elements: usize,
+    pub layer_elements: usize,
+    /// The grid's axes.
+    pub frame: [[f64; 3]; 3],
+    /// Sharp corners, curves and patches of the boundary.
+    pub corners: usize,
+    pub curves: usize,
+    pub patches: usize,
+    /// Optimiser sweeps of the last build.
+    pub sweeps: usize,
+    /// Smallest and mean corner scaled Jacobian.
+    pub min_scaled_jacobian: f64,
+    pub mean_scaled_jacobian: f64,
+    /// Exact volume of the elements and the boundary's enclosed volume.
+    pub volume: f64,
+    pub surface_volume: f64,
+    /// The largest distance of a boundary vertex from its corner, curve or patch: zero when
+    /// every boundary vertex lies on the geometry; where the geometry is sharper than valid
+    /// elements can follow (blades, slits), the mesh stays this far off.
+    pub boundary_distance: f64,
+    /// Why earlier builds were refined.
+    pub retries: Vec<String>,
+}
+
+/// An all-hex8 mesh.
+#[derive(Debug, Clone, PartialEq)]
+pub struct OctreeHexMesh {
+    pub nodes: Vec<[f64; 3]>,
+    /// Model node order; every element is proven valid.
+    pub hexes: Vec<[u32; 8]>,
+    pub boundary: Vec<BoundaryFace>,
+    pub report: OctreeReport,
+}
+
+/// A rebuild after a failure may have at most this many times the grid cells of the build
+/// before…
+pub const RETRY_GROWTH: usize = 4;
+/// … or this many cells, whichever is more.
+const RETRY_MIN_CELLS: usize = 250_000;
+/// The fitting weights of the successive optimisation stages.
+const FIT_STAGES: &[f64] = &[1e-3, 1e-2, 1e-1, 1.0, 10.0, 100.0];
+/// A boundary vertex goes exactly onto its feature only if its elements keep at least this
+/// corner scaled Jacobian.
+const MIN_SNAP_QUALITY: f64 = 0.05;
+/// Elements below this corner scaled Jacobian are improved by compass search.
+const IMPROVE_BELOW: f64 = 0.3;
+/// Closure rounds of the grid (each splits at least one leaf; the cell cap bounds them too).
+const CLOSURE_ROUNDS: usize = 1000;
+/// Rounds of the core clean-up (each drops at least one element).
+const CORE_ROUNDS: usize = 1000;
+
+/// A failure that refining at `at` may cure.
+struct Retry {
+    at: Vec<P>,
+    why: String,
+}
+
+enum Attempt {
+    Done(Box<OctreeHexMesh>),
+    Retry(Retry),
+}
+
+/// Mesh the solid bounded by `surface` with hexahedra of about `field`'s size (the root cell
+/// size is the field at the centre of the part's box; finer cells follow the part's curvature,
+/// thickness and narrow faces). See the module docs.
+pub fn octree_hex_mesh(surface: &BoundarySurface, field: &dyn SizeField, options: &OctreeOptions) -> Result<OctreeHexMesh> {
+    check_options(options)?;
+    let g = Geometry::new(surface, options.sharp_angle, options.blade_angle)?;
+    let frame = choose_frame(&g);
+    let mut lo = [f64::INFINITY; 3];
+    let mut hi = [f64::NEG_INFINITY; 3];
+    for &p in surface.points() {
+        let u = to_frame(&frame, p);
+        for k in 0..3 {
+            lo[k] = lo[k].min(u[k]);
+            hi[k] = hi[k].max(u[k]);
+        }
+    }
+    let centre_u = [0, 1, 2].map(|k| 0.5 * (lo[k] + hi[k]));
+    let h = checked_size(field, from_frame(&frame, centre_u))?;
+    let mut roots = [0i64; 3];
+    for k in 0..3 {
+        let n = ((hi[k] - lo[k]) / h + 1.0).ceil();
+        if !(n.is_finite() && n >= 1.0 && n <= options.max_elements as f64) {
+            return Err(Error::Limit(format!("a size of {h} needs {n} root cells along an axis")));
+        }
+        roots[k] = n as i64;
+    }
+    let origin = [0, 1, 2].map(|k| centre_u[k] - 0.5 * roots[k] as f64 * h);
+    let ctx = Context {
+        g: &g,
+        field,
+        options,
+        frame,
+        origin,
+        h,
+        roots,
+        thickness: triangle_thickness(&g),
+        outside_memo: RefCell::new(HashMap::new()),
+        last_leaves: Cell::new(0),
+        too_thin: Cell::new(f64::INFINITY),
+    };
+    let mut hints: Vec<P> = Vec::new();
+    let mut retries: Vec<String> = Vec::new();
+    for attempt in 0..options.attempts {
+        let cap = if attempt == 0 {
+            options.max_elements
+        } else {
+            options.max_elements.min(RETRY_GROWTH.saturating_mul(ctx.last_leaves.get()).max(RETRY_MIN_CELLS))
+        };
+        let built = match ctx.build(&hints, cap) {
+            Err(Error::Limit(why)) if attempt > 0 => {
+                retries.push(format!("refining for another build needs more than {cap} grid cells ({why})"));
+                break;
+            }
+            other => other?,
+        };
+        match built {
+            Attempt::Done(mut m) => {
+                m.report.attempts = attempt + 1;
+                m.report.retries = retries;
+                return Ok(*m);
+            }
+            Attempt::Retry(r) => {
+                retries.push(r.why);
+                hints.extend(r.at);
+            }
+        }
+    }
+    Err(Error::Failed(format!("octree hex meshing failed after {} builds: {}", retries.len(), retries.join("; "))))
+}
+
+fn check_options(o: &OctreeOptions) -> Result<()> {
+    let ok = o.max_depth <= 8
+        && o.curvature_levels <= o.max_depth
+        && o.max_elements >= 1
+        && o.curvature_angle.is_finite()
+        && o.curvature_angle > 0.0
+        && o.thickness_cells.is_finite()
+        && o.thickness_cells >= 0.0
+        && o.core_offset.is_finite()
+        && (0.0..1.0).contains(&o.core_offset)
+        && (1..=16).contains(&o.attempts)
+        && o.sharp_angle.is_finite()
+        && o.sharp_angle > 0.0
+        && o.blade_angle.is_finite()
+        && o.blade_angle > o.sharp_angle
+        && o.blade_angle <= 180.0
+        && o.max_sweeps <= 100_000;
+    if ok { Ok(()) } else { Err(Error::Invalid(format!("octree options out of range: {o:?}"))) }
+}
+
+fn from_frame(f: &[[f64; 3]; 3], u: P) -> P {
+    [0, 1, 2].map(|k| u[0] * f[0][k] + u[1] * f[1][k] + u[2] * f[2][k])
+}
+
+fn to_frame(f: &[[f64; 3]; 3], p: P) -> P {
+    [dot(p, f[0]), dot(p, f[1]), dot(p, f[2])]
+}
+
+/// The grid axes (rows): the coordinate axes, or the normals of the two largest orthogonal planar
+/// patches when more boundary area is aligned with them.
+fn choose_frame(g: &Geometry) -> [[f64; 3]; 3] {
+    let tol = 1f64.to_radians().sin();
+    let normal = |t: u32| g.normals.get(t as usize).copied().unwrap_or([0.0; 3]);
+    let mut planar: Vec<(f64, usize, P)> = Vec::new();
+    for (i, p) in g.patches.iter().enumerate() {
+        let Some(&t0) = p.triangles.first() else { continue };
+        let n0 = normal(t0);
+        if p.triangles.iter().all(|&t| norm(cross(normal(t), n0)) < tol) {
+            planar.push((p.area, i, n0));
+        }
+    }
+    planar.sort_by(|a, b| b.0.total_cmp(&a.0).then(a.1.cmp(&b.1)));
+    let identity = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]];
+    let Some(&(_, _, a0)) = planar.first() else { return identity };
+    let Some(&(_, _, a1)) = planar.iter().find(|x| dot(x.2, a0).abs() < tol) else { return identity };
+    let a1 = {
+        let v = sub(a1, a0.map(|x| x * dot(a1, a0)));
+        let l = norm(v);
+        v.map(|x| x / l)
+    };
+    let candidate = [a0, a1, cross(a0, a1)];
+    let score = |f: &[[f64; 3]; 3]| -> f64 {
+        (0..g.surface.triangles().len() as u32)
+            .filter(|&t| f.iter().any(|a| norm(cross(normal(t), *a)) < tol))
+            .map(|t| {
+                let [a, b, c] = g.triangle(t);
+                0.5 * norm(cross(sub(b, a), sub(c, a)))
+            })
+            .sum()
+    };
+    if score(&candidate) > score(&identity) * (1.0 + 1e-9) { candidate } else { identity }
+}
+
+/// The distance through the part from each boundary triangle: a ray from its centre along the
+/// inward normal to the next triangle (`∞` when none).
+fn triangle_thickness(g: &Geometry) -> Vec<f64> {
+    let far = 2.0 * norm(sub(g.bounds.hi, g.bounds.lo)) + 1.0;
+    (0..g.surface.triangles().len() as u32)
+        .map(|t| {
+            let n = g.normals.get(t as usize).copied().unwrap_or([0.0; 3]);
+            g.ray(g.triangle_centre(t), n.map(|x| -x), far, t).map_or(f64::INFINITY, |(_, s)| s)
+        })
+        .collect()
+}
+
+fn centroid(points: &[P], hx: &[u32; 8]) -> P {
+    [0, 1, 2].map(|k| hx.iter().map(|&v| points.get(v as usize).map_or(0.0, |p| p[k])).sum::<f64>() / 8.0)
+}
+
+fn edge_mean(points: &[P], hx: &[u32; 8]) -> f64 {
+    let x = hx.map(|v| points.get(v as usize).copied().unwrap_or([0.0; 3]));
+    let edges = [(0, 1), (1, 2), (2, 3), (3, 0), (4, 5), (5, 6), (6, 7), (7, 4), (0, 4), (1, 5), (2, 6), (3, 7)];
+    let m = edges.iter().map(|&(a, b)| norm(sub(x[a], x[b]))).sum::<f64>() / 12.0;
+    if m > 0.0 { m } else { 1.0 }
+}
+
+/// The part being meshed and the state carried from one build to the next.
+struct Context<'a> {
+    g: &'a Geometry<'a>,
+    field: &'a dyn SizeField,
+    options: &'a OctreeOptions,
+    frame: [[f64; 3]; 3],
+    origin: P,
+    h: f64,
+    roots: [i64; 3],
+    thickness: Vec<f64>,
+    outside_memo: RefCell<HashMap<grid::Key, bool>>,
+    /// Grid leaves of the last build.
+    last_leaves: Cell<usize>,
+    /// The thinnest wall the deepest cells cannot cross in `thickness_cells` (∞ when none).
+    too_thin: Cell<f64>,
+}
+
+/// The mesh while its boundary layers are built: points, elements, their target sizes and how
+/// each point may move.
+struct Work {
+    points: Vec<P>,
+    hexes: Vec<[u32; 8]>,
+    size: Vec<f64>,
+    role: Vec<Constraint>,
+    core_count: usize,
+}
+
+impl Context<'_> {
+    fn edge(&self, level: u8) -> f64 {
+        self.h / 3f64.powi(i32::from(level))
+    }
+
+    /// The world box (inflated by `pad` cell edges) and the centre of a cell.
+    fn cell_box(&self, (level, i): grid::Key, pad: f64) -> (Aabb, P) {
+        let e = self.edge(level);
+        let mut b = Aabb::EMPTY;
+        for c in 0..8 {
+            b.add(from_frame(&self.frame, [0, 1, 2].map(|k| self.origin[k] + (i[k] as f64 + ((c >> k) & 1) as f64) * e)));
+        }
+        let centre = from_frame(&self.frame, [0, 1, 2].map(|k| self.origin[k] + (i[k] as f64 + 0.5) * e));
+        let p = pad * e;
+        (Aabb { lo: b.lo.map(|x| x - p), hi: b.hi.map(|x| x + p) }, centre)
+    }
+
+    fn needs_split(&self, key: grid::Key) -> Result<bool> {
+        let e = self.edge(key.0);
+        let (bx, centre) = self.cell_box(key, 0.05);
+        if e > checked_size(self.field, centre)? * (1.0 + 1e-9) {
+            return Ok(true);
+        }
+        let reach = 0.5 * 3f64.sqrt() * e * 1.05;
+        let cos = self.options.curvature_angle.to_radians().cos();
+        let mut first: BTreeMap<u32, P> = BTreeMap::new();
+        for t in self.g.overlapping(&bx) {
+            let [a, b, c] = self.g.triangle(t);
+            if geom::dist2(centre, geom::closest_on_triangle(centre, a, b, c)) > reach * reach {
+                continue;
+            }
+            let n = self.g.normals.get(t as usize).copied().unwrap_or([0.0; 3]);
+            let face = self.g.surface.tri_face().get(t as usize).copied().unwrap_or(0);
+            let f = *first.entry(face).or_insert(n);
+            if key.0 < self.options.curvature_levels && dot(f, n) < cos {
+                return Ok(true);
+            }
+            if let Some(&th) = self.thickness.get(t as usize)
+                && th < self.options.thickness_cells * e
+            {
+                if key.0 + 1 >= self.options.max_depth && th < self.options.thickness_cells * e / 3.0 {
+                    self.too_thin.set(self.too_thin.get().min(th));
+                }
+                return Ok(true);
+            }
+            // Narrow patches (a hole's wall, a thin rim) need cells across them.
+            if self.g.tri_patch.get(t as usize).and_then(|&p| self.g.patches.get(p as usize)).is_some_and(|pt| pt.width < 2.0 * e) {
+                return Ok(true);
+            }
+        }
+        Ok(false)
+    }
+
+    /// Whether a cell lies wholly outside the part (its centre outside and farther from the
+    /// boundary than its half diagonal): no core element can come from it, so its transition
+    /// pattern does not matter and it is not meshed. Memoised.
+    fn outside(&self, key: grid::Key) -> bool {
+        if let Some(&v) = self.outside_memo.borrow().get(&key) {
+            return v;
+        }
+        let (_, centre) = self.cell_box(key, 0.0);
+        let half = 0.5 * 3f64.sqrt() * self.edge(key.0);
+        let far = self.g.nearest(centre, |_| true).is_some_and(|x| x.2 > half * half * (1.0 + 1e-9));
+        let v = far && self.g.inside(centre) == Some(false);
+        self.outside_memo.borrow_mut().insert(key, v);
+        v
+    }
+
+    /// The level-`level` cells within one cell of a hint point.
+    fn hint_cells(&self, hints: &[P], level: u8) -> HashSet<grid::Key> {
+        let e = self.edge(level);
+        let mut out = HashSet::new();
+        for &p in hints {
+            let u = to_frame(&self.frame, p);
+            let i = [0, 1, 2].map(|k| ((u[k] - self.origin[k]) / e).floor());
+            if !i.iter().all(|x| x.is_finite() && x.abs() < 1e15) {
+                continue;
+            }
+            for d in 0..27i64 {
+                out.insert((level, [i[0] as i64 + d % 3 - 1, i[1] as i64 + (d / 3) % 3 - 1, i[2] as i64 + d / 9 - 1]));
+            }
+        }
+        out
+    }
+
+    /// The refined, balanced and closed grid; leaves near `hints` (where an earlier build
+    /// failed) are split once more.
+    fn tree(&self, hints: &[P], max_leaves: usize) -> Result<Tree> {
+        let mut tree = Tree::new(self.roots, self.options.max_depth, max_leaves)?;
+        for level in 0..self.options.max_depth {
+            for key in tree.sorted_leaves() {
+                if key.0 == level && self.needs_split(key)? {
+                    tree.split(key)?;
+                }
+            }
+        }
+        if !hints.is_empty() {
+            let near: Vec<HashSet<grid::Key>> = (0..self.options.max_depth).map(|l| self.hint_cells(hints, l)).collect();
+            let split: Vec<grid::Key> =
+                tree.sorted_leaves().into_iter().filter(|k| near.get(usize::from(k.0)).is_some_and(|s| s.contains(k))).collect();
+            for key in split {
+                tree.split(key)?;
+            }
+        }
+        tree.close(CLOSURE_ROUNDS, &|k| self.outside(k))?;
+        Ok(tree)
+    }
+
+    fn build(&self, hints: &[P], max_leaves: usize) -> Result<Attempt> {
+        let g = self.g;
+        let tree = self.tree(hints, max_leaves)?;
+        self.last_leaves.set(tree.leaf_count());
+        let lattice = tree.mesh(&|k| self.outside(k))?;
+        if lattice.hexes.len() > self.options.max_elements {
+            return Err(Error::Limit(format!("the grid has {} elements; the limit is {}", lattice.hexes.len(), self.options.max_elements)));
+        }
+        let unit = self.h / tree.span(0) as f64;
+        let points: Vec<P> = lattice.points.iter().map(|q| from_frame(&self.frame, [0, 1, 2].map(|k| self.origin[k] + q[k] as f64 * unit))).collect();
+
+        // Core.
+        let dist: Vec<f64> = points.iter().map(|&p| g.nearest(p, |_| true).map_or(0.0, |x| x.2.sqrt())).collect();
+        let keep = self.core(&lattice, &points, &dist);
+        let core: Vec<usize> = (0..lattice.hexes.len()).filter(|&i| keep[i]).collect();
+        if core.is_empty() {
+            let at = (0..g.surface.triangles().len() as u32).map(|t| g.triangle_centre(t)).collect();
+            return Ok(Attempt::Retry(Retry { at, why: "no grid element lies inside the part".into() }));
+        }
+        let mut w = Work {
+            hexes: core.iter().map(|&i| lattice.hexes[i]).collect(),
+            size: core.iter().map(|&i| self.edge(lattice.level[i])).collect(),
+            role: Vec::new(),
+            points,
+            core_count: core.len(),
+        };
+        let qb = topo::boundary(&w.hexes);
+        let q: Vec<[u32; 4]> = qb.iter().map(|x| x.1).collect();
+        if let Some(r) = check_topology(g, &q) {
+            let thin = self.too_thin.get();
+            if thin.is_finite() {
+                return Err(Error::Failed(format!(
+                    "{}: the part has walls {thin:.3} thick, thinner than {} cells of the deepest level ({:.3}); use a smaller size, more levels, or a structured (swept) hex mesh",
+                    r.why,
+                    self.options.thickness_cells,
+                    self.edge(self.options.max_depth)
+                )));
+            }
+            return Ok(Attempt::Retry(r));
+        }
+
+        // Feature map, padding and pillows.
+        let fmap = {
+            let pts = &w.points;
+            match map::map_features(g, &q, &|v| pts.get(v as usize).copied().unwrap_or([0.0; 3])) {
+                Ok(m) => m,
+                Err(f) => return Ok(Attempt::Retry(Retry { at: f.at, why: f.why })),
+            }
+        };
+        self.layers(&mut w, &qb, &fmap, &dist)?;
+        if w.hexes.len() > self.options.max_elements {
+            return Err(Error::Limit(format!("the mesh needs {} elements; the limit is {}", w.hexes.len(), self.options.max_elements)));
+        }
+
+        // Fit, then improve and repair.
+        let mut opt = Optimizer::new(g, std::mem::take(&mut w.points), &w.hexes, &w.size, w.role.clone());
+        let mut ran = 0;
+        let sweeps = (self.options.max_sweeps / FIT_STAGES.len()).max(1);
+        for &weight in FIT_STAGES {
+            opt.fit = weight;
+            ran = opt.run(sweeps, 1e-4);
+        }
+        let (_, boundary_distance) = opt.snap(MIN_SNAP_QUALITY);
+        let inverted: Vec<P> = w
+            .hexes
+            .iter()
+            .filter(|hx| smooth::min_scaled_jacobian(&hx.map(|v| opt.points[v as usize])) <= 0.0)
+            .map(|hx| centroid(&opt.points, hx))
+            .collect();
+        if !inverted.is_empty() {
+            let why = format!("{} elements stay inverted after untangling", inverted.len());
+            return Ok(Attempt::Retry(Retry { at: inverted, why }));
+        }
+        opt.improve_worst(IMPROVE_BELOW, 4);
+        let uncertified = opt.repair_certified(6);
+        if !uncertified.is_empty() {
+            let at: Vec<P> = uncertified.iter().filter_map(|&e| w.hexes.get(e as usize)).map(|hx| centroid(&opt.points, hx)).collect();
+            let why = format!("{} elements fail the certified validity check", at.len());
+            return Ok(Attempt::Retry(Retry { at, why }));
+        }
+        w.points = opt.points;
+        self.finish(&w, &tree, ran, boundary_distance)
+    }
+
+    /// The core: elements wholly inside the part, cleaned to a 2-manifold boundary whose
+    /// vertices can all move outward.
+    fn core(&self, lattice: &grid::LatticeMesh, points: &[P], dist: &[f64]) -> Vec<bool> {
+        let inside: Vec<bool> = points.iter().map(|&p| self.g.inside(p) == Some(true)).collect();
+        let mut keep: Vec<bool> = lattice
+            .hexes
+            .iter()
+            .zip(&lattice.level)
+            .map(|(hx, &l)| {
+                let floor = self.options.core_offset * self.edge(l);
+                hx.iter().all(|&v| inside.get(v as usize) == Some(&true) && dist.get(v as usize).is_some_and(|&d| d >= floor))
+            })
+            .collect();
+        for _ in 0..CORE_ROUNDS {
+            let mut drop = topo::nonmanifold(&lattice.hexes, &keep);
+            if drop.is_empty() {
+                drop = topo::hidden_vertices(&lattice.hexes, &keep, &|v| points.get(v as usize).copied().unwrap_or([0.0; 3]), 0.1);
+            }
+            if drop.is_empty() {
+                break;
+            }
+            for d in drop {
+                if let Some(k) = keep.get_mut(d as usize) {
+                    *k = false;
+                }
+            }
+        }
+        keep
+    }
+
+    /// The boundary layers: padding on the core's boundary quads `qb`, whose outer copies take the
+    /// map's roles and start halfway to the boundary along the core's normals, then a pillow
+    /// under each patch, whose copies start a quarter of the way into their elements. Every
+    /// element starts valid or nearly so; the fit untangles the rest.
+    fn layers(&self, w: &mut Work, qb: &[(topo::FaceRef, [u32; 4])], fmap: &map::FeatureMap, dist: &[f64]) -> Result<()> {
+        let q: Vec<[u32; 4]> = qb.iter().map(|x| x.1).collect();
+        let core_count = w.core_count;
+        let copy = topo::pad(&mut w.points, &mut w.hexes, &q)?;
+        w.role = vec![Constraint::Free; w.points.len()];
+        let mut normal: BTreeMap<u32, P> = BTreeMap::new();
+        for f in &q {
+            let p = f.map(|v| w.points[v as usize]);
+            let n = cross(sub(p[2], p[0]), sub(p[3], p[1]));
+            let l = norm(n);
+            if l > 0.0 {
+                for &v in f {
+                    let e = normal.entry(v).or_insert([0.0; 3]);
+                    for k in 0..3 {
+                        e[k] += n[k] / l;
+                    }
+                }
+            }
+        }
+        for (&v, &c) in &copy {
+            if let Some(r) = fmap.role.get(&v) {
+                w.role[c as usize] = *r;
+            }
+            let n = normal.get(&v).copied().unwrap_or([0.0; 3]);
+            let l = norm(n);
+            let d = 0.5 * dist.get(v as usize).copied().unwrap_or(0.0);
+            if l > 0.0 {
+                w.points[c as usize] = [0, 1, 2].map(|k| w.points[v as usize][k] + d * n[k] / l);
+            }
+        }
+        for ((owner, _), _) in qb {
+            w.size.push(0.5 * w.size.get(*owner as usize).copied().unwrap_or(self.h));
+        }
+        let patches: BTreeSet<u32> = fmap.quad_patch.iter().copied().collect();
+        for p in patches {
+            let inside: Vec<bool> =
+                (0..w.hexes.len()).map(|i| i >= core_count && i < core_count + q.len() && fmap.quad_patch.get(i - core_count) == Some(&p)).collect();
+            let before = w.hexes.len();
+            let centroids: Vec<P> = w.hexes.iter().map(|hx| centroid(&w.points, hx)).collect();
+            let copies = topo::pillow(&mut w.points, &mut w.hexes, &inside)?;
+            w.role.resize(w.points.len(), Constraint::Free);
+            for pc in &copies {
+                let orig = w.points[pc.original as usize];
+                let n = pc.elements.len().max(1) as f64;
+                let mut c = [0.0; 3];
+                for &e in &pc.elements {
+                    let ce = centroids.get(e as usize).copied().unwrap_or(orig);
+                    for k in 0..3 {
+                        c[k] += ce[k] / n;
+                    }
+                }
+                w.points[pc.copy as usize] = [0, 1, 2].map(|k| orig[k] + 0.25 * (c[k] - orig[k]));
+                w.role[pc.copy as usize] = match w.role.get(pc.original as usize).copied().unwrap_or(Constraint::Free) {
+                    Constraint::Free => Constraint::Free,
+                    _ => Constraint::Patch(p),
+                };
+            }
+            for hx in &w.hexes[before..] {
+                w.size.push(edge_mean(&w.points, hx));
+            }
+        }
+        Ok(())
+    }
+
+    /// The gate and the output: every element proven valid, nodes renumbered, boundary faces with
+    /// their B-rep faces.
+    fn finish(&self, w: &Work, tree: &Tree, sweeps: usize, boundary_distance: f64) -> Result<Attempt> {
+        let g = self.g;
+        // One certified pass: every element's quality (the same metric as mesh.quality and
+        // mesh.hex.generate report) and its Jacobian check for the gate.
+        let mut bad: Vec<P> = Vec::new();
+        let mut checks = Vec::with_capacity(w.hexes.len());
+        let (mut sj_min, mut sj_sum) = (f64::INFINITY, 0.0);
+        for (e, hx) in w.hexes.iter().enumerate() {
+            let q = element_quality(ElementKind::Hex8, &hx.map(|v| w.points[v as usize]))?;
+            let check = q.jacobian.ok_or_else(|| Error::Failed("an element was not certified".into()))?;
+            if check.validity != Validity::Valid {
+                bad.push(centroid(&w.points, hx));
+            }
+            let sj = q.get(Metric::ScaledJacobian).unwrap_or(-1.0);
+            sj_min = sj_min.min(sj);
+            sj_sum += sj;
+            checks.push((e as u32 + 1, check));
+        }
+        if !bad.is_empty() {
+            let why = format!("{} elements fail the certified validity check", bad.len());
+            return Ok(Attempt::Retry(Retry { at: bad, why }));
+        }
+        // The volume the mesh's own boundary encloses equals the elements' total exactly when they
+        // tile it without overlapping.
+        let faces = topo::boundary(&w.hexes);
+        let enclosed: f64 = faces.iter().map(|(_, quad)| face_flux(&quad.map(|v| w.points[v as usize]))).sum();
+        let gate =
+            validity_gate_checks(checks.iter().map(|(l, c)| (*l, c)), &GateOptions { expected_volume: Some(enclosed), volume_tolerance: 1e-9 })?
+                .into_result()?;
+        let volume = gate.volume;
+        let mut new_id = vec![u32::MAX; w.points.len()];
+        let mut nodes = Vec::new();
+        let mut hexes = Vec::with_capacity(w.hexes.len());
+        for hx in &w.hexes {
+            let mut o = [0u32; 8];
+            for (k, &v) in hx.iter().enumerate() {
+                if new_id[v as usize] == u32::MAX {
+                    new_id[v as usize] = u32::try_from(nodes.len()).map_err(|_| Error::Limit("too many nodes".into()))?;
+                    nodes.push(w.points[v as usize]);
+                }
+                o[k] = new_id[v as usize];
+            }
+            hexes.push(o);
+        }
+        let mut boundary = Vec::new();
+        for ((e, f), quad) in faces {
+            let patch = quad.iter().find_map(|&v| match w.role.get(v as usize) {
+                Some(Constraint::Patch(p)) => Some(*p),
+                _ => None,
+            });
+            let c = [0, 1, 2].map(|k| quad.iter().map(|&v| w.points[v as usize][k]).sum::<f64>() / 4.0);
+            let brep_face = match patch.and_then(|p| g.patches.get(p as usize)) {
+                Some(pt) if pt.faces.len() == 1 => pt.faces[0],
+                _ => {
+                    let keep = |t: u32| patch.is_none_or(|p| g.tri_patch.get(t as usize) == Some(&p));
+                    g.nearest(c, keep).and_then(|(t, _, _)| g.surface.tri_face().get(t as usize).copied()).unwrap_or(0)
+                }
+            };
+            boundary.push(BoundaryFace { element: e, face: f + 1, brep_face });
+        }
+        let report = OctreeReport {
+            attempts: 0,
+            grid_cells: tree.leaf_count(),
+            core_elements: w.core_count,
+            layer_elements: w.hexes.len() - w.core_count,
+            frame: self.frame,
+            corners: g.corners.len(),
+            curves: g.curves.len(),
+            patches: g.patches.len(),
+            sweeps,
+            min_scaled_jacobian: sj_min,
+            mean_scaled_jacobian: sj_sum / w.hexes.len().max(1) as f64,
+            volume,
+            surface_volume: g.surface.volume(),
+            boundary_distance,
+            retries: Vec::new(),
+        };
+        Ok(Attempt::Done(Box::new(OctreeHexMesh { nodes, hexes, boundary, report })))
+    }
+}
+
+/// The core boundary's components and genus against the part's shells; on a mismatch, a retry
+/// that refines near the whole boundary.
+fn check_topology(g: &Geometry, quads: &[[u32; 4]]) -> Option<Retry> {
+    let edges = topo::quad_edges(quads);
+    let mut comp: Vec<usize> = (0..quads.len()).collect();
+    for qs in edges.values() {
+        for w in qs.windows(2) {
+            let (a, b) = (topo::find(&mut comp, w[0] as usize), topo::find(&mut comp, w[1] as usize));
+            comp[a] = b;
+        }
+    }
+    // Per component: vertices, edges, faces.
+    let mut stats: BTreeMap<usize, (BTreeSet<u32>, i64, i64)> = BTreeMap::new();
+    for (i, q) in quads.iter().enumerate() {
+        let s = stats.entry(topo::find(&mut comp, i)).or_default();
+        s.2 += 1;
+        s.0.extend(q.iter().copied());
+    }
+    for qs in edges.values() {
+        if let Some(&q0) = qs.first()
+            && let Some(s) = stats.get_mut(&topo::find(&mut comp, q0 as usize))
+        {
+            s.1 += 1;
+        }
+    }
+    let mut got: Vec<i64> = stats.values().map(|s| (2 - (s.0.len() as i64 - s.1 + s.2)) / 2).collect();
+    let mut want: Vec<i64> = g.surface.shells().iter().map(|s| i64::from(s.genus)).collect();
+    got.sort_unstable();
+    want.sort_unstable();
+    if got == want {
+        return None;
+    }
+    let at = (0..g.surface.triangles().len() as u32).map(|t| g.triangle_centre(t)).collect();
+    Some(Retry { at, why: format!("the core's boundary has components of genus {got:?}, the part {want:?}") })
+}

crates/mesh/src/octree/smooth.rs

@@ -0,0 +1,641 @@
+//! Untangling, fitting and shape optimisation of the hexahedral mesh.
+//!
+//! # Energy
+//! The foldover-free mapping energy of V. Garanzha, I. Kaporin, L. Kudryavtseva, F. Protais,
+//! N. Ray, D. Sokolov, *Foldover-free maps in 50 lines of code*, ACM TOG 40(4):102, 2021 (paper
+//! only; the authors' code is AGPL and is not read). Each hexahedron is measured at nine points,
+//! its eight corners (the trapezoidal rule the paper uses for quads) and its centre, where a
+//! twisted element folds first: `J` is the sample's three edge directions against the ideal cube
+//! of the element's target edge `s`, `D = det J`, and the sample costs
+//!
+//! `f_ε + λ g_ε`, `f_ε = |J|³_F / χ(D, ε)`, `g_ε = (D² + 1) / χ(D, ε)`,
+//! `χ(D, ε) = (D + √(ε² + D²)) / 2`,
+//!
+//! weighted by `s³`. `f` is the shape (conformal) term, `g` the volume term. The paper's shape
+//! term is `|J|²/D^{2/3}`; its 3/2 power `|J|³/D`, the classical polyconvex distortion in three
+//! dimensions (J. M. Ball, 1977), is minimal at the same ideal shapes and needs only square
+//! roots, so the optimiser gives the same bits on every platform (no `pow` of the system maths
+//! library). `χ` keeps both finite for inverted samples, and `ε` follows the paper's rule (7),
+//! `ε = √(10⁻¹² + 0.04 min(0, min D)²)`, so it vanishes once every sample is positive and the
+//! energy becomes a barrier against inversion.
+//!
+//! A vertex on a corner, curve or patch adds `fit · 8 s̄ · dist²` (`s̄` the mean target edge of its
+//! elements) to its distance from that feature: with `fit` raised stage by stage the boundary is
+//! drawn onto the geometry while the barrier keeps every element valid (Gao, Shen, Panozzo 2019,
+//! §3.5). With `fit = 0` constrained vertices are held on their features instead (each step is
+//! restricted to the tangent and projected back).
+//!
+//! # Minimisation
+//! Nonlinear Gauss–Seidel: each vertex in turn takes a Newton step on its own energy with the
+//! paper's positive approximation of the Hessian (the second derivatives of `χ` and `det` are
+//! dropped, §2.3, App. A), accepted only if the energy decreased (halving line search). The paper
+//! minimises the same energy globally with L-BFGS; the vertex sweep is the local form of the same
+//! descent. Then a compass search raises the worst corner scaled Jacobians, and another repairs
+//! elements whose trilinear map the certified check does not prove positive. Sweeps, line-search
+//! steps, compass steps and rounds are capped.
+
+use fieldcraft_model::ElementKind;
+
+use crate::quality::{Validity, element_validity};
+use crate::vec3::{P, cross, dot, norm, sub};
+
+use super::geom::dist2;
+
+use super::geom::Geometry;
+
+/// How a vertex may move.
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub(crate) enum Constraint {
+    Free,
+    /// On a corner of the part (the boundary surface's vertex).
+    Corner(u32),
+    /// On a feature curve.
+    Curve(u32),
+    /// On a patch of the boundary.
+    Patch(u32),
+}
+
+/// Model-order corner `c` of a hex8: its three neighbours along x, y, z and the signs that make
+/// the reference corner frame right-handed.
+const CORNER_FRAMES: [([usize; 3], [f64; 3]); 8] = [
+    ([1, 3, 4], [1.0, 1.0, 1.0]),
+    ([0, 2, 5], [-1.0, 1.0, 1.0]),
+    ([3, 1, 6], [-1.0, -1.0, 1.0]),
+    ([2, 0, 7], [1.0, -1.0, 1.0]),
+    ([5, 7, 0], [1.0, 1.0, -1.0]),
+    ([4, 6, 1], [-1.0, 1.0, -1.0]),
+    ([7, 5, 2], [-1.0, -1.0, -1.0]),
+    ([6, 4, 3], [1.0, -1.0, -1.0]),
+];
+
+/// Weight of the volume term.
+const LAMBDA: f64 = 0.1;
+
+fn chi(d: f64, eps: f64) -> f64 {
+    0.5 * (d + (eps * eps + d * d).sqrt())
+}
+
+fn chi_d(d: f64, eps: f64) -> f64 {
+    0.5 * (1.0 + d / (eps * eps + d * d).sqrt())
+}
+
+/// The corner Jacobian (columns) of corner `c` of `x` against an ideal cube of edge `s`.
+fn corner_jacobian(x: &[P; 8], c: usize, s: f64) -> [P; 3] {
+    let (n, sg) = CORNER_FRAMES[c];
+    [0, 1, 2].map(|a| sub(x[n[a]], x[c]).map(|v| v * sg[a] / s))
+}
+
+/// The sample points of the energy: the eight corners and the centre. Column `a` of sample `s`'s
+/// Jacobian is `Σ_k W[s][a][k] x_k` (before dividing by the target edge). The centre keeps
+/// twisted elements (positive corners, folded inside) out of reach, as the nine-point measure of
+/// ParaView and of Zhang et al.'s meshers does.
+const SAMPLES: usize = 9;
+
+fn sample_weights(s: usize) -> [[f64; 8]; 3] {
+    let mut w = [[0.0; 8]; 3];
+    if s < 8 {
+        let (n, sg) = CORNER_FRAMES[s];
+        for a in 0..3 {
+            w[a][n[a]] += sg[a];
+            w[a][s] -= sg[a];
+        }
+    } else {
+        w = [
+            [-0.25, 0.25, 0.25, -0.25, -0.25, 0.25, 0.25, -0.25],
+            [-0.25, -0.25, 0.25, 0.25, -0.25, -0.25, 0.25, 0.25],
+            [-0.25, -0.25, -0.25, -0.25, 0.25, 0.25, 0.25, 0.25],
+        ];
+    }
+    w
+}
+
+/// The Jacobian (columns) of sample `smp` of `x` against an ideal cube of edge `s`.
+fn sample_jacobian(x: &[P; 8], smp: usize, s: f64) -> [P; 3] {
+    let w = sample_weights(smp);
+    [0, 1, 2].map(|a| {
+        let mut c = [0.0; 3];
+        for (k, p) in x.iter().enumerate() {
+            for r in 0..3 {
+                c[r] += w[a][k] * p[r];
+            }
+        }
+        c.map(|v| v / s)
+    })
+}
+
+fn det3(j: &[P; 3]) -> f64 {
+    dot(j[0], cross(j[1], j[2]))
+}
+
+/// The cofactor columns: `∂ det / ∂ column a`.
+fn cofactors(j: &[P; 3]) -> [P; 3] {
+    [cross(j[1], j[2]), cross(j[2], j[0]), cross(j[0], j[1])]
+}
+
+/// The energy of one corner.
+fn corner_energy(j: &[P; 3], eps: f64) -> f64 {
+    let d = det3(j);
+    let c = chi(d, eps);
+    let fro: f64 = j.iter().map(|v| dot(*v, *v)).sum();
+    fro * fro.sqrt() / c + LAMBDA * (d * d + 1.0) / c
+}
+
+/// The scaled Jacobian of each corner, smallest (`min over corners of det / product of edge
+/// lengths`), the usual hex quality measure.
+pub(crate) fn min_scaled_jacobian(x: &[P; 8]) -> f64 {
+    (0..8)
+        .map(|c| {
+            let j = corner_jacobian(x, c, 1.0);
+            let l = norm(j[0]) * norm(j[1]) * norm(j[2]);
+            if l > 0.0 { det3(&j) / l } else { -1.0 }
+        })
+        .fold(f64::INFINITY, f64::min)
+}
+
+/// The optimiser's state.
+pub(crate) struct Optimizer<'a> {
+    pub geom: &'a Geometry<'a>,
+    pub points: Vec<P>,
+    pub hexes: &'a [[u32; 8]],
+    /// Target edge per element.
+    pub size: &'a [f64],
+    pub constraint: Vec<Constraint>,
+    /// Elements around each vertex with the vertex's local position.
+    around: Vec<Vec<(u32, u8)>>,
+    pub eps: f64,
+    /// Weight of the fitting term; 0 holds boundary vertices on their features exactly.
+    pub fit: f64,
+}
+
+impl<'a> Optimizer<'a> {
+    pub fn new(geom: &'a Geometry<'a>, points: Vec<P>, hexes: &'a [[u32; 8]], size: &'a [f64], constraint: Vec<Constraint>) -> Optimizer<'a> {
+        let mut around: Vec<Vec<(u32, u8)>> = vec![Vec::new(); points.len()];
+        for (e, h) in hexes.iter().enumerate() {
+            for (k, &v) in h.iter().enumerate() {
+                if let Some(a) = around.get_mut(v as usize) {
+                    a.push((e as u32, k as u8));
+                }
+            }
+        }
+        Optimizer { geom, points, hexes, size, constraint, around, eps: 1.0, fit: 0.0 }
+    }
+
+    fn element(&self, e: u32) -> ([P; 8], f64) {
+        let h = self.hexes.get(e as usize).copied().unwrap_or([0; 8]);
+        (h.map(|v| self.points.get(v as usize).copied().unwrap_or([0.0; 3])), self.size.get(e as usize).copied().unwrap_or(1.0))
+    }
+
+    /// The smallest normalised corner determinant over the mesh, and the total energy.
+    pub fn measure(&self) -> (f64, f64) {
+        let (mut m, mut total) = (f64::INFINITY, 0.0);
+        for e in 0..self.hexes.len() as u32 {
+            let (x, s) = self.element(e);
+            for c in 0..SAMPLES {
+                let j = sample_jacobian(&x, c, s);
+                m = m.min(det3(&j));
+                total += s * s * s * corner_energy(&j, self.eps);
+            }
+        }
+        (m, total)
+    }
+
+    /// Where vertex `v`'s feature is closest to `p` (`None` for a free vertex).
+    pub fn target(&self, v: u32, p: P) -> Option<P> {
+        match self.constraint.get(v as usize).copied().unwrap_or(Constraint::Free) {
+            Constraint::Free => None,
+            Constraint::Corner(c) => Some(self.geom.xyz(c)),
+            Constraint::Curve(c) => self.geom.nearest_on_curve(c, p).map(|x| x.0),
+            Constraint::Patch(pt) => self.geom.nearest(p, |t| self.geom.tri_patch.get(t as usize) == Some(&pt)).map(|x| x.1),
+        }
+    }
+
+    /// The fitting weight of vertex `v`: the fit factor times its elements' size (so a distance
+    /// of one element weighs like the shape of its corners), ten times more for a corner.
+    fn fit_weight(&self, v: u32) -> f64 {
+        let a = self.around.get(v as usize).map_or(&[][..], |a| a.as_slice());
+        if a.is_empty() || self.fit == 0.0 {
+            return 0.0;
+        }
+        let s = a.iter().map(|&(e, _)| self.size.get(e as usize).copied().unwrap_or(1.0)).sum::<f64>() / a.len() as f64;
+        let corner = if matches!(self.constraint.get(v as usize), Some(Constraint::Corner(_))) { 10.0 } else { 1.0 };
+        8.0 * self.fit * corner * s
+    }
+
+    /// The energy of the corners that involve vertex `v`, with its fitting term.
+    fn local_energy(&self, v: u32) -> f64 {
+        let mut total = 0.0;
+        if self.fit > 0.0 {
+            let p = self.points.get(v as usize).copied().unwrap_or([0.0; 3]);
+            if let Some(t) = self.target(v, p) {
+                total += self.fit_weight(v) * dist2(p, t);
+            }
+        }
+        for &(e, k) in self.around.get(v as usize).map_or(&[][..], |a| a.as_slice()) {
+            let (x, s) = self.element(e);
+            let (n, _) = CORNER_FRAMES[k as usize];
+            for c in [k as usize, n[0], n[1], n[2], 8] {
+                total += s * s * s * corner_energy(&sample_jacobian(&x, c, s), self.eps);
+            }
+        }
+        total
+    }
+
+    /// The gradient and approximate Hessian of the local energy at vertex `v`.
+    fn derivatives(&self, v: u32) -> (P, [[f64; 3]; 3]) {
+        let mut g = [0.0; 3];
+        let mut h = [[0.0; 3]; 3];
+        for &(e, k) in self.around.get(v as usize).map_or(&[][..], |a| a.as_slice()) {
+            let (x, s) = self.element(e);
+            let k = k as usize;
+            let w = s * s * s;
+            let (nk, _) = CORNER_FRAMES[k];
+            for c in [k, nk[0], nk[1], nk[2], 8] {
+                // ∂ column a / ∂ x_v = coef[a] I.
+                let wt = sample_weights(c);
+                let coef: [f64; 3] = [0, 1, 2].map(|a| wt[a][k] / s);
+                let j = sample_jacobian(&x, c, s);
+                let d = det3(&j);
+                let ch = chi(d, self.eps);
+                let chd = chi_d(d, self.eps);
+                let fro: f64 = j.iter().map(|q| dot(*q, *q)).sum();
+                let cof = cofactors(&j);
+                let root = fro.sqrt();
+                // dE/dJ = A J + B cof(J).
+                let a_coef = 3.0 * root / ch;
+                let f_d = -fro * root * chd / (ch * ch);
+                let g_d = 2.0 * d / ch - (d * d + 1.0) * chd / (ch * ch);
+                let b_coef = f_d + LAMBDA * g_d;
+                let mut dd = [0.0; 3];
+                let mut sum_c2 = 0.0;
+                for a in 0..3 {
+                    if coef[a] == 0.0 {
+                        continue;
+                    }
+                    sum_c2 += coef[a] * coef[a];
+                    for r in 0..3 {
+                        g[r] += w * coef[a] * (a_coef * j[a][r] + b_coef * cof[a][r]);
+                        dd[r] += coef[a] * cof[a][r];
+                    }
+                }
+                // Positive part of the second derivative along D, χ'' and det'' dropped.
+                let f_dd = 2.0 * fro * root * chd * chd / (ch * ch * ch);
+                let g_dd = 2.0 / ch - 4.0 * d * chd / (ch * ch) + 2.0 * (d * d + 1.0) * chd * chd / (ch * ch * ch);
+                let k_dd = (f_dd + LAMBDA * g_dd).max(0.0);
+                for r in 0..3 {
+                    h[r][r] += w * a_coef * sum_c2;
+                    for q in 0..3 {
+                        h[r][q] += w * k_dd * dd[r] * dd[q];
+                    }
+                }
+            }
+        }
+        if self.fit > 0.0 {
+            let p = self.points.get(v as usize).copied().unwrap_or([0.0; 3]);
+            if let Some(t) = self.target(v, p) {
+                let w = self.fit_weight(v);
+                for r in 0..3 {
+                    g[r] += 2.0 * w * (p[r] - t[r]);
+                    h[r][r] += 2.0 * w;
+                }
+            }
+        }
+        (g, h)
+    }
+
+    /// Project a point onto vertex `v`'s constraint.
+    fn project(&self, v: u32, p: P) -> P {
+        self.target(v, p).unwrap_or(p)
+    }
+
+    /// The direction a constrained vertex may move along: the step restricted to its tangent.
+    fn restrict(&self, v: u32, d: P) -> P {
+        let p = self.points.get(v as usize).copied().unwrap_or([0.0; 3]);
+        match self.constraint.get(v as usize).copied().unwrap_or(Constraint::Free) {
+            Constraint::Free => d,
+            Constraint::Corner(_) => [0.0; 3],
+            Constraint::Curve(c) => {
+                let Some((_, i, _)) = self.geom.nearest_on_curve(c, p) else { return [0.0; 3] };
+                let Some(cv) = self.geom.curves.get(c as usize) else { return [0.0; 3] };
+                let (a, b) = (cv.vertices.get(i).copied().unwrap_or(0), cv.vertices.get(i + 1).copied().unwrap_or(0));
+                let t = sub(self.geom.xyz(b), self.geom.xyz(a));
+                let l2 = dot(t, t);
+                if l2 > 0.0 { t.map(|x| x * dot(d, t) / l2) } else { [0.0; 3] }
+            }
+            Constraint::Patch(pt) => {
+                let Some((t, _, _)) = self.geom.nearest(p, |t| self.geom.tri_patch.get(t as usize) == Some(&pt)) else { return [0.0; 3] };
+                let n = self.geom.normals.get(t as usize).copied().unwrap_or([0.0; 3]);
+                let dn = dot(d, n);
+                [0, 1, 2].map(|k| d[k] - dn * n[k])
+            }
+        }
+    }
+
+    /// One Newton step of vertex `v`; whether it moved.
+    fn relax(&mut self, v: u32) -> bool {
+        let soft = self.fit > 0.0;
+        if (!soft && matches!(self.constraint.get(v as usize), Some(Constraint::Corner(_))))
+            || self.around.get(v as usize).is_none_or(|a| a.is_empty())
+        {
+            return false;
+        }
+        let (g, h) = self.derivatives(v);
+        let Some(step) = solve3(h, g.map(|x| -x)) else { return false };
+        let step = if soft { step } else { self.restrict(v, step) };
+        if !step.iter().all(|x| x.is_finite()) || norm(step) == 0.0 {
+            return false;
+        }
+        let old = self.points.get(v as usize).copied().unwrap_or([0.0; 3]);
+        let e0 = self.local_energy(v);
+        let mut tau = 1.0;
+        for _ in 0..12 {
+            let moved = [0, 1, 2].map(|k| old[k] + tau * step[k]);
+            let trial = if soft { moved } else { self.project(v, moved) };
+            if let Some(p) = self.points.get_mut(v as usize) {
+                *p = trial;
+            }
+            let e1 = self.local_energy(v);
+            if e1 < e0 && e1.is_finite() {
+                return true;
+            }
+            tau *= 0.5;
+        }
+        if let Some(p) = self.points.get_mut(v as usize) {
+            *p = old;
+        }
+        false
+    }
+
+    /// Set the regularisation from the smallest corner determinant (Garanzha et al. 2021: ε
+    /// shrinks as the mesh untangles); that determinant and the energy at the new ε.
+    fn update_eps(&mut self) -> (f64, f64) {
+        let (min_d, energy) = self.measure();
+        let old = self.eps;
+        self.eps = (1e-12 + 0.04 * min_d.min(0.0).powi(2)).sqrt();
+        (min_d, if self.eps == old { energy } else { self.measure().1 })
+    }
+
+    /// Sweep all vertices until every corner is positive and the energy settles, or `max_sweeps`;
+    /// the sweeps run.
+    pub fn run(&mut self, max_sweeps: usize, tolerance: f64) -> usize {
+        let n = self.points.len() as u32;
+        let (_, mut energy) = self.update_eps();
+        let mut sweeps = 0;
+        for sweep in 0..max_sweeps {
+            let mut moved = 0usize;
+            for v in 0..n {
+                if self.relax(v) {
+                    moved += 1;
+                }
+            }
+            let old_eps = self.eps;
+            let (min_d, e) = self.update_eps();
+            sweeps = sweep + 1;
+            let settled = old_eps == self.eps && (energy - e).abs() <= tolerance * energy.abs();
+            energy = e;
+            if moved == 0 || (min_d > 0.0 && settled) {
+                break;
+            }
+        }
+        sweeps
+    }
+}
+
+impl Optimizer<'_> {
+    /// The smallest corner scaled Jacobian of the elements around vertex `v`.
+    fn local_quality(&self, v: u32) -> f64 {
+        self.around
+            .get(v as usize)
+            .map_or(f64::INFINITY, |a| a.iter().map(|&(e, _)| min_scaled_jacobian(&self.element(e).0)).fold(f64::INFINITY, f64::min))
+    }
+
+    /// Move every constrained vertex exactly onto its feature where its elements keep a corner
+    /// scaled Jacobian of at least `floor`; the others stay where the fit left them. Returns
+    /// how many stayed off and the largest distance of a constrained vertex from its feature.
+    pub fn snap(&mut self, floor: f64) -> (usize, f64) {
+        let (mut off, mut worst) = (0usize, 0.0f64);
+        for v in 0..self.points.len() as u32 {
+            let p = self.points.get(v as usize).copied().unwrap_or([0.0; 3]);
+            let Some(t) = self.target(v, p) else { continue };
+            if p == t {
+                continue;
+            }
+            if let Some(x) = self.points.get_mut(v as usize) {
+                *x = t;
+            }
+            if self.local_quality(v) < floor {
+                if let Some(x) = self.points.get_mut(v as usize) {
+                    *x = p;
+                }
+                off += 1;
+                worst = worst.max(dist2(p, t).sqrt());
+            }
+        }
+        (off, worst)
+    }
+
+    /// Raise the worst elements: every vertex of an element whose corner scaled Jacobian is
+    /// below `threshold` moves, within its constraint, by a compass search that maximises the
+    /// smallest scaled Jacobian around it (steps shrink by halves; at most `passes` rounds over
+    /// the worst elements, 12 step sizes per vertex).
+    pub fn improve_worst(&mut self, threshold: f64, passes: usize) {
+        for _ in 0..passes {
+            let mut verts: Vec<u32> = Vec::new();
+            for (e, h) in self.hexes.iter().enumerate() {
+                if min_scaled_jacobian(&self.element(e as u32).0) < threshold {
+                    verts.extend(h.iter().copied());
+                }
+            }
+            verts.sort_unstable();
+            verts.dedup();
+            if verts.is_empty() {
+                break;
+            }
+            let mut moved = false;
+            for &v in &verts {
+                if matches!(self.constraint.get(v as usize), Some(Constraint::Free)) {
+                    moved |= self.compass(v, &|o: &Self, v| o.local_quality(v));
+                }
+            }
+            if !moved {
+                break;
+            }
+        }
+    }
+
+    /// Compass search on vertex `v` (within its constraint) for the largest `quality`; steps
+    /// start at 0.15 of its largest element size and halve when no direction improves, 12 times
+    /// at most. Whether it moved.
+    fn compass(&mut self, v: u32, quality: &dyn Fn(&Self, u32) -> f64) -> bool {
+        let p0 = self.points.get(v as usize).copied().unwrap_or([0.0; 3]);
+        let reach = self
+            .around
+            .get(v as usize)
+            .map_or(0.0, |a| a.iter().map(|&(e, _)| self.size.get(e as usize).copied().unwrap_or(0.0)).fold(0.0, f64::max));
+        let mut step = 0.15 * reach;
+        let mut best = quality(self, v);
+        let mut cur = p0;
+        let mut moved = false;
+        for _ in 0..12 {
+            if step <= 1e-9 * reach.max(1e-300) {
+                break;
+            }
+            let mut improved = false;
+            for d in [[1.0, 0.0, 0.0], [-1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, 1.0], [0.0, 0.0, -1.0]] {
+                let dir = self.restrict(v, d.map(|x| x * step));
+                if norm(dir) == 0.0 {
+                    continue;
+                }
+                let trial = self.project(v, [0, 1, 2].map(|k| cur[k] + dir[k]));
+                if let Some(p) = self.points.get_mut(v as usize) {
+                    *p = trial;
+                }
+                let q = quality(self, v);
+                if q > best + 1e-12 {
+                    best = q;
+                    cur = trial;
+                    improved = true;
+                    moved = true;
+                } else if let Some(p) = self.points.get_mut(v as usize) {
+                    *p = cur;
+                }
+            }
+            if !improved {
+                step *= 0.5;
+            }
+        }
+        if let Some(p) = self.points.get_mut(v as usize) {
+            *p = cur;
+        }
+        moved
+    }
+
+    /// The smallest certified Jacobian ratio of the elements around `v` (negative where the
+    /// trilinear map folds somewhere, even with positive corners).
+    fn certified_quality(&self, v: u32) -> f64 {
+        self.around.get(v as usize).map_or(f64::INFINITY, |a| {
+            a.iter()
+                .map(|&(e, _)| {
+                    let x = self.element(e).0;
+                    if min_scaled_jacobian(&x) <= 0.0 {
+                        return -1e9;
+                    }
+                    element_validity(ElementKind::Hex8, &x).map_or(-1.0, |c| match c.validity {
+                        Validity::Valid => c.ratio.max(1e-300),
+                        _ => c.ratio.min(0.0) - 1e-12,
+                    })
+                })
+                .fold(f64::INFINITY, f64::min)
+        })
+    }
+
+    /// Repair elements whose trilinear map is not proven positive although their corners are:
+    /// the vertices of each such element (free ones, and boundary ones along their features)
+    /// move by compass search to raise the smallest certified bound around them; `rounds` rounds
+    /// at most. Returns the elements still not proven valid.
+    pub fn repair_certified(&mut self, rounds: usize) -> Vec<u32> {
+        let check = |o: &Self, e: u32| element_validity(ElementKind::Hex8, &o.element(e).0).is_ok_and(|c| c.validity == Validity::Valid);
+        let mut bad: Vec<u32> = (0..self.hexes.len() as u32).filter(|&e| !check(self, e)).collect();
+        for _ in 0..rounds {
+            if bad.is_empty() {
+                break;
+            }
+            let mut verts: Vec<u32> = bad.iter().flat_map(|&e| self.hexes.get(e as usize).copied().unwrap_or([0; 8])).collect();
+            verts.sort_unstable();
+            verts.dedup();
+            for v in verts {
+                if !matches!(self.constraint.get(v as usize), Some(Constraint::Corner(_)) | None) {
+                    self.compass(v, &|o: &Self, v| o.certified_quality(v));
+                }
+            }
+            let mut around: Vec<u32> = bad
+                .iter()
+                .flat_map(|&e| self.hexes.get(e as usize).copied().unwrap_or([0; 8]))
+                .flat_map(|v| self.around.get(v as usize).map_or(Vec::new(), |a| a.iter().map(|x| x.0).collect()))
+                .collect();
+            around.sort_unstable();
+            around.dedup();
+            bad = around.into_iter().filter(|&e| !check(self, e)).collect();
+        }
+        bad
+    }
+}
+
+/// Solve the symmetric 3 × 3 system `h x = b` (Cramer), `None` when singular.
+fn solve3(h: [[f64; 3]; 3], b: P) -> Option<P> {
+    let d = dot(h[0], cross(h[1], h[2]));
+    if !(d.abs() > 1e-300) || !d.is_finite() {
+        return None;
+    }
+    let col = |k: usize| [h[0][k], h[1][k], h[2][k]];
+    let (c0, c1, c2) = (col(0), col(1), col(2));
+    Some([dot(b, cross(c1, c2)) / d, dot(c0, cross(b, c2)) / d, dot(c0, cross(c1, b)) / d])
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::octree::geom::tests::box_surface;
+
+    #[test]
+    fn corner_frames_are_right_handed_on_the_unit_cube() {
+        let cube: [P; 8] = [[0., 0., 0.], [1., 0., 0.], [1., 1., 0.], [0., 1., 0.], [0., 0., 1.], [1., 0., 1.], [1., 1., 1.], [0., 1., 1.]];
+        for c in 0..8 {
+            let j = corner_jacobian(&cube, c, 1.0);
+            assert!((det3(&j) - 1.0).abs() < 1e-15, "corner {c}");
+        }
+        assert!((min_scaled_jacobian(&cube) - 1.0).abs() < 1e-15);
+    }
+
+    #[test]
+    fn a_tangled_interior_vertex_is_untangled() {
+        // A 2 × 2 × 2 block of cubes with the centre vertex pushed through a face; the boundary
+        // fixed.
+        let s = box_surface([0.0; 3], [2.0; 3]);
+        let g = Geometry::new(&s, 30.0, 160.0).unwrap();
+        let id = |x: u32, y: u32, z: u32| x + 3 * (y + 3 * z);
+        let mut p = Vec::new();
+        for z in 0..3 {
+            for y in 0..3 {
+                for x in 0..3 {
+                    p.push([x as f64, y as f64, z as f64]);
+                }
+            }
+        }
+        let mut h = Vec::new();
+        for z in 0..2 {
+            for y in 0..2 {
+                for x in 0..2 {
+                    h.push([
+                        id(x, y, z),
+                        id(x + 1, y, z),
+                        id(x + 1, y + 1, z),
+                        id(x, y + 1, z),
+                        id(x, y, z + 1),
+                        id(x + 1, y, z + 1),
+                        id(x + 1, y + 1, z + 1),
+                        id(x, y + 1, z + 1),
+                    ]);
+                }
+            }
+        }
+        let centre = id(1, 1, 1) as usize;
+        p[centre] = [1.7, 1.9, 0.2];
+        let size = vec![1.0; h.len()];
+        // Boundary vertices slide on the box's faces (hard constraints, patch of the nearest
+        // face), the centre is free.
+        let cons: Vec<Constraint> = (0..p.len())
+            .map(|i| {
+                if i == centre {
+                    return Constraint::Free;
+                }
+                let t = g.nearest(p[i], |_| true).unwrap().0;
+                Constraint::Patch(g.tri_patch[t as usize])
+            })
+            .collect();
+        let mut o = Optimizer::new(&g, p, &h, &size, cons);
+        assert!(o.measure().0 < 0.0);
+        let st = o.run(200, 1e-9);
+        assert!(o.measure().0 > 0.9, "{st:?}");
+        assert!(sub(o.points[centre], [1.0; 3]).iter().all(|c| c.abs() < 1e-3), "{:?}", o.points[centre]);
+    }
+}

crates/mesh/src/octree/templates.rs

@@ -0,0 +1,374 @@
+//! The 3-refinement transition templates: how a leaf of the balanced tree is cut into
+//! hexahedra so that it conforms to finer neighbours.
+//!
+//! Every vertex of a leaf is *marked* when a finer leaf touches it. Under strong balance the
+//! marked set of a leaf is one of five patterns (after the closure in [`super::grid`] has split
+//! every leaf that had another one), each with one template:
+//!
+//! | marks | template | hexahedra |
+//! |---|---|---|
+//! | none | the leaf itself | 1 |
+//! | one corner | corner cube and three blocks fanning to the far corner | 4 |
+//! | the two ends of an edge | three cubes along the edge and two four-hex slabs | 11 |
+//! | the four corners of a face | nine cubes on the face, nine above them, four on top | 22 |
+//! | all eight | the 3 × 3 × 3 split | 27 |
+//!
+//! A face of a leaf is cut according to its own marks only (none, one corner, one edge, all
+//! four), and an edge into thirds at its marked ends, so two leaves sharing a face cut it the same
+//! way: the mesh is conforming by construction. This is the node-marking form of Schneiders'
+//! 3-refinement (R. Schneiders, *Refining quadrilateral and hexahedral element meshes*, NUMGRID
+//! 1996; Y. Ito, A. M. Shih, B. K. Soni, *Octree-based reasonable-quality hexahedral mesh
+//! generation using a new set of refinement templates*, IJNME 77(13):1809–1833, 2009). The
+//! templates are our own construction from those papers' description: the face template's top
+//! block is the 3 × 3 → 1 transition traced from its dual sheets, the edge template is two
+//! copies of the planar 3 → 1 transition swept across the leaf, and the tests prove every
+//! template's faces, volume and validity.
+//!
+//! Coordinates are integers in sixths of the leaf's edge, so every vertex shared with a neighbour
+//! lies on the global integer lattice of [`super::grid`] and coincident vertices merge exactly.
+
+/// A point of a leaf in sixths of its edge, `0..=6` per axis.
+pub(crate) type Local = [i64; 3];
+/// A hexahedron in the model's node order (corners 1-2-3-4 then 5-6-7-8 above them).
+pub(crate) type Hex = [Local; 8];
+
+/// The marked-corner pattern of a leaf, up to symmetry.
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub(crate) enum Pattern {
+    None,
+    Corner,
+    Edge,
+    Face,
+    Full,
+}
+
+/// A symmetry of the cube: coordinate `k` of the image is coordinate `perm[k]` of the point,
+/// reflected (`6 − v`) when `flip[k]`.
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub(crate) struct Sym {
+    perm: [usize; 3],
+    flip: [bool; 3],
+}
+
+impl Sym {
+    pub(crate) fn apply(&self, p: Local) -> Local {
+        let mut q = [0; 3];
+        for (k, out) in q.iter_mut().enumerate() {
+            let v = p[self.perm[k]];
+            *out = if self.flip[k] { 6 - v } else { v };
+        }
+        q
+    }
+}
+
+/// The 48 symmetries in a fixed order.
+pub(crate) fn symmetries() -> Vec<Sym> {
+    let perms = [[0, 1, 2], [0, 2, 1], [1, 0, 2], [1, 2, 0], [2, 0, 1], [2, 1, 0]];
+    let mut out = Vec::with_capacity(48);
+    for perm in perms {
+        for bits in 0..8u8 {
+            out.push(Sym { perm, flip: [bits & 1 != 0, bits & 2 != 0, bits & 4 != 0] });
+        }
+    }
+    out
+}
+
+/// The leaf corner `c` (bit 0: x, bit 1: y, bit 2: z) in sixths.
+pub(crate) fn corner(c: usize) -> Local {
+    [(c & 1) as i64 * 6, ((c >> 1) & 1) as i64 * 6, ((c >> 2) & 1) as i64 * 6]
+}
+
+/// The corner index of a corner point.
+fn corner_index(p: Local) -> usize {
+    usize::from(p[0] == 6) | (usize::from(p[1] == 6) << 1) | (usize::from(p[2] == 6) << 2)
+}
+
+/// The image of a corner mask under a symmetry.
+fn map_mask(sym: &Sym, mask: u8) -> u8 {
+    (0..8).filter(|c| mask & (1 << c) != 0).fold(0u8, |m, c| m | (1 << corner_index(sym.apply(corner(c)))))
+}
+
+/// The pattern of a corner mask, or `None` when no template fits it.
+pub(crate) fn classify(mask: u8) -> Option<Pattern> {
+    let corners: Vec<Local> = (0..8).filter(|c| mask & (1 << c) != 0).map(corner).collect();
+    let shared = |n: usize| (0..3).filter(|&k| corners.iter().all(|p| p[k] == corners[0][k])).count() == n;
+    match corners.len() {
+        0 => Some(Pattern::None),
+        1 => Some(Pattern::Corner),
+        2 if shared(2) => Some(Pattern::Edge),
+        4 if shared(1) => Some(Pattern::Face),
+        8 => Some(Pattern::Full),
+        _ => None,
+    }
+}
+
+/// The canonical mask of a pattern: corner 0; edge 0–1 (along x); face z = 0.
+fn canonical_mask(p: Pattern) -> u8 {
+    match p {
+        Pattern::None => 0,
+        Pattern::Corner => 0b1,
+        Pattern::Edge => 0b11,
+        Pattern::Face => 0b1111,
+        Pattern::Full => 0xff,
+    }
+}
+
+/// The first symmetry (in [`symmetries`] order) taking the pattern's canonical mask to `mask`.
+pub(crate) fn orientation(syms: &[Sym], pattern: Pattern, mask: u8) -> Option<Sym> {
+    let c = canonical_mask(pattern);
+    syms.iter().find(|s| map_mask(s, c) == mask).copied()
+}
+
+/// Thirds to sixths.
+fn t(x: i64, y: i64, z: i64) -> Local {
+    [2 * x, 2 * y, 2 * z]
+}
+
+/// A thirds-aligned cube with lower corner `(x, y, z)` (in thirds) and edge `e` thirds.
+fn cube(x: i64, y: i64, z: i64, e: i64) -> Hex {
+    let (a, b) = (t(x, y, z), t(x + e, y + e, z + e));
+    [
+        [a[0], a[1], a[2]],
+        [b[0], a[1], a[2]],
+        [b[0], b[1], a[2]],
+        [a[0], b[1], a[2]],
+        [a[0], a[1], b[2]],
+        [b[0], a[1], b[2]],
+        [b[0], b[1], b[2]],
+        [a[0], b[1], b[2]],
+    ]
+}
+
+/// The planar 3 → 1 transition in `(x, s)`: x in thirds, `s` in halves of the transition height,
+/// three segments at `s = 0` and one at `s = 2`.
+const THREE_TO_ONE: [[(i64, i64); 4]; 4] =
+    [[(0, 0), (1, 0), (1, 1), (0, 2)], [(1, 0), (2, 0), (2, 1), (1, 1)], [(2, 0), (3, 0), (3, 2), (2, 1)], [(1, 1), (2, 1), (3, 2), (0, 2)]];
+
+/// The canonical template of a pattern, in sixths of the leaf.
+pub(crate) fn template(p: Pattern) -> Vec<Hex> {
+    match p {
+        Pattern::None => vec![cube(0, 0, 0, 3)],
+        Pattern::Full => (0..27).map(|k| cube(k % 3, (k / 3) % 3, k / 9, 1)).collect(),
+        Pattern::Corner => vec![
+            cube(0, 0, 0, 1),
+            [t(1, 0, 0), t(3, 0, 0), t(3, 3, 0), t(1, 1, 0), t(1, 0, 1), t(3, 0, 3), t(3, 3, 3), t(1, 1, 1)],
+            [t(0, 1, 0), t(1, 1, 0), t(3, 3, 0), t(0, 3, 0), t(0, 1, 1), t(1, 1, 1), t(3, 3, 3), t(0, 3, 3)],
+            [t(0, 0, 1), t(1, 0, 1), t(1, 1, 1), t(0, 1, 1), t(0, 0, 3), t(3, 0, 3), t(3, 3, 3), t(0, 3, 3)],
+        ],
+        Pattern::Edge => {
+            let mut out: Vec<Hex> = (0..3).map(|i| cube(i, 0, 0, 1)).collect();
+            // Two slabs over the edge column's faces y = 1 and z = 1, each the 3 → 1 transition
+            // swept from the leaf face to the diagonal plane y = z.
+            for q in THREE_TO_ONE {
+                let b = |(x, s): (i64, i64), w: i64| t(x, 1 + s, w * (1 + s));
+                let c = |(x, s): (i64, i64), w: i64| t(x, w * (1 + s), 1 + s);
+                out.push([b(q[0], 0), b(q[1], 0), b(q[2], 0), b(q[3], 0), b(q[0], 1), b(q[1], 1), b(q[2], 1), b(q[3], 1)]);
+                out.push([c(q[0], 0), c(q[1], 0), c(q[2], 0), c(q[3], 0), c(q[0], 1), c(q[1], 1), c(q[2], 1), c(q[3], 1)]);
+            }
+            out
+        }
+        Pattern::Face => {
+            let mut out: Vec<Hex> = (0..9).map(|k| cube(k % 3, k / 3, 0, 1)).collect();
+            // The layer above: its top is the corners raised to the leaf's top, the side nodes of
+            // the 3 → 1 transitions at two thirds, the inside nodes at half height.
+            let top = |a: i64, b: i64| -> Local {
+                let edge = |v: i64| v == 0 || v == 3;
+                match (edge(a), edge(b)) {
+                    (true, true) => [2 * a, 2 * b, 6],
+                    (false, false) => [2 * a, 2 * b, 3],
+                    _ => [2 * a, 2 * b, 4],
+                }
+            };
+            for k in 0..9 {
+                let (i, j) = (k % 3, k / 3);
+                out.push([
+                    t(i, j, 1),
+                    t(i + 1, j, 1),
+                    t(i + 1, j + 1, 1),
+                    t(i, j + 1, 1),
+                    top(i, j),
+                    top(i + 1, j),
+                    top(i + 1, j + 1),
+                    top(i, j + 1),
+                ]);
+            }
+            // The four on top: the y-arch (low), the x-arch (high) and the two side blocks.
+            out.push([[2, 2, 3], [4, 2, 3], [4, 4, 3], [2, 4, 3], [2, 0, 4], [4, 0, 4], [4, 6, 4], [2, 6, 4]]);
+            out.push([[2, 0, 4], [4, 0, 4], [4, 6, 4], [2, 6, 4], [0, 0, 6], [6, 0, 6], [6, 6, 6], [0, 6, 6]]);
+            out.push([[2, 2, 3], [2, 4, 3], [2, 6, 4], [2, 0, 4], [0, 2, 4], [0, 4, 4], [0, 6, 6], [0, 0, 6]]);
+            out.push([[4, 2, 3], [4, 4, 3], [4, 6, 4], [4, 0, 4], [6, 2, 4], [6, 4, 4], [6, 6, 6], [6, 0, 6]]);
+            out
+        }
+    }
+}
+
+/// Six times the signed volume of the corner tetrahedron at node 0 (exact on the lattice).
+pub(crate) fn corner0_orientation(h: &[[i64; 3]; 8]) -> i128 {
+    let e = |j: usize| [0, 1, 2].map(|k| i128::from(h[j][k] - h[0][k]));
+    let (a, b, c) = (e(1), e(3), e(4));
+    a[0] * (b[1] * c[2] - b[2] * c[1]) - a[1] * (b[0] * c[2] - b[2] * c[0]) + a[2] * (b[0] * c[1] - b[1] * c[0])
+}
+
+/// The hexahedron with its orientation reversed (the same cell, mirrored node order).
+pub(crate) fn reversed<T: Copy>(h: [T; 8]) -> [T; 8] {
+    [h[0], h[3], h[2], h[1], h[4], h[7], h[6], h[5]]
+}
+
+#[cfg(test)]
+mod tests {
+    use std::collections::HashMap;
+
+    use super::*;
+    use crate::quality::{Validity, element_validity};
+    use fieldcraft_model::ElementKind;
+
+    /// The face of the leaf a quad lies on (axis, value), if any.
+    fn on_leaf_face(q: &[Local; 4]) -> Option<(usize, i64)> {
+        (0..3).find_map(|k| [0, 6].into_iter().find(|&v| q.iter().all(|p| p[k] == v)).map(|v| (k, v)))
+    }
+
+    fn sorted(mut q: [Local; 4]) -> [Local; 4] {
+        q.sort();
+        q
+    }
+
+    /// The quads a leaf face must carry for its four corner marks, as sorted vertex sets, on the
+    /// face `(axis, value)`.
+    fn face_template(axis: usize, value: i64, marks: [bool; 4]) -> Vec<[Local; 4]> {
+        // Face coordinates (u, v) are the other two axes in increasing order; corners in order
+        // (0,0), (6,0), (0,6), (6,6).
+        let (ua, va) = match axis {
+            0 => (1, 2),
+            1 => (0, 2),
+            _ => (0, 1),
+        };
+        let lift = |u: i64, v: i64| {
+            let mut p = [0; 3];
+            p[axis] = value;
+            p[ua] = u;
+            p[va] = v;
+            p
+        };
+        // Build in a canonical frame, then map the 2D symmetry that takes the canonical marks to
+        // these marks.
+        let count = marks.iter().filter(|m| **m).count();
+        let canon: Vec<[(i64, i64); 4]> = match count {
+            0 => vec![[(0, 0), (6, 0), (6, 6), (0, 6)]],
+            1 => vec![[(0, 0), (2, 0), (2, 2), (0, 2)], [(2, 0), (6, 0), (6, 6), (2, 2)], [(0, 2), (2, 2), (6, 6), (0, 6)]],
+            2 => {
+                let mut v: Vec<[(i64, i64); 4]> = (0..3).map(|i| [(2 * i, 0), (2 * i + 2, 0), (2 * i + 2, 2), (2 * i, 2)]).collect();
+                v.extend([
+                    [(0, 2), (2, 2), (2, 4), (0, 6)],
+                    [(2, 2), (4, 2), (4, 4), (2, 4)],
+                    [(4, 2), (6, 2), (6, 6), (4, 4)],
+                    [(2, 4), (4, 4), (6, 6), (0, 6)],
+                ]);
+                v
+            }
+            _ => (0..9)
+                .map(|k| {
+                    [(2 * (k % 3), 2 * (k / 3)), (2 * (k % 3) + 2, 2 * (k / 3)), (2 * (k % 3) + 2, 2 * (k / 3) + 2), (2 * (k % 3), 2 * (k / 3) + 2)]
+                })
+                .collect(),
+        };
+        let canon_marks: [bool; 4] = match count {
+            0 => [false; 4],
+            1 => [true, false, false, false],
+            2 => [true, true, false, false],
+            _ => [true; 4],
+        };
+        let maps: [fn(i64, i64) -> (i64, i64); 8] = [
+            |u, v| (u, v),
+            |u, v| (6 - u, v),
+            |u, v| (u, 6 - v),
+            |u, v| (6 - u, 6 - v),
+            |u, v| (v, u),
+            |u, v| (6 - v, u),
+            |u, v| (v, 6 - u),
+            |u, v| (6 - v, 6 - u),
+        ];
+        let corners = [(0, 0), (6, 0), (0, 6), (6, 6)];
+        let idx = |c: (i64, i64)| corners.iter().position(|&x| x == c).unwrap();
+        let m =
+            maps.iter().find(|m| (0..4).all(|c| canon_marks[c] == marks[idx(m(corners[c].0, corners[c].1))])).expect("a face symmetry for the marks");
+        canon.iter().map(|q| sorted(q.map(|(u, v)| m(u, v)).map(|(u, v)| lift(u, v)))).collect()
+    }
+
+    fn face_corner_marks(axis: usize, value: i64, mask: u8) -> [bool; 4] {
+        let (ua, va) = match axis {
+            0 => (1, 2),
+            1 => (0, 2),
+            _ => (0, 1),
+        };
+        [(0, 0), (6, 0), (0, 6), (6, 6)].map(|(u, v)| {
+            let mut p = [0; 3];
+            p[axis] = value;
+            p[ua] = u;
+            p[va] = v;
+            mask & (1 << corner_index(p)) != 0
+        })
+    }
+
+    /// Every template under every symmetry: valid elements, the leaf's volume, closed (each inner
+    /// face shared by two elements) and, on each leaf face, exactly the face template of that
+    /// face's marks.
+    #[test]
+    fn templates_are_valid_closed_and_conform_to_their_face_marks() {
+        let syms = symmetries();
+        for pattern in [Pattern::None, Pattern::Corner, Pattern::Edge, Pattern::Face, Pattern::Full] {
+            for sym in &syms {
+                let mask = map_mask(sym, canonical_mask(pattern));
+                assert_eq!(classify(mask), Some(pattern));
+                let hexes: Vec<Hex> = template(pattern)
+                    .into_iter()
+                    .map(|h| {
+                        let h = h.map(|p| sym.apply(p));
+                        if corner0_orientation(&h) < 0 { reversed(h) } else { h }
+                    })
+                    .collect();
+                let mut volume = 0.0;
+                let mut faces: HashMap<[Local; 4], usize> = HashMap::new();
+                for h in &hexes {
+                    let xyz: Vec<[f64; 3]> = h.iter().map(|p| p.map(|v| v as f64)).collect();
+                    let check = element_validity(ElementKind::Hex8, &xyz).unwrap();
+                    assert_eq!(check.validity, Validity::Valid, "{pattern:?} {sym:?} {h:?}");
+                    volume += check.volume;
+                    for f in ElementKind::Hex8.faces() {
+                        *faces.entry(sorted([h[f[0]], h[f[1]], h[f[2]], h[f[3]]])).or_default() += 1;
+                    }
+                }
+                assert!((volume - 216.0).abs() < 1e-9, "{pattern:?}: volume {volume}");
+                let mut outer: HashMap<(usize, i64), Vec<[Local; 4]>> = HashMap::new();
+                for (q, n) in &faces {
+                    match on_leaf_face(q) {
+                        Some(f) => {
+                            assert_eq!(*n, 1);
+                            outer.entry(f).or_default().push(*q);
+                        }
+                        None => assert_eq!(*n, 2, "{pattern:?} {sym:?}: inner face {q:?} used {n} times"),
+                    }
+                }
+                for axis in 0..3 {
+                    for value in [0, 6] {
+                        let mut got = outer.remove(&(axis, value)).unwrap_or_default();
+                        let mut want = face_template(axis, value, face_corner_marks(axis, value, mask));
+                        got.sort();
+                        want.sort();
+                        assert_eq!(got, want, "{pattern:?} {sym:?} face {axis}={value}");
+                    }
+                }
+            }
+        }
+    }
+
+    #[test]
+    fn patterns_without_a_template_are_recognised() {
+        assert_eq!(classify(0b1001), None, "diagonal of a face");
+        assert_eq!(classify(0b1000_0001), None, "opposite corners");
+        assert_eq!(classify(0b0111), None, "three of a face");
+        assert_eq!(classify(0b0110_1001), None, "four not on a face");
+        assert_eq!(classify(0b0111_1111), None);
+        assert_eq!(classify(0b0011_0011), Some(Pattern::Face));
+        assert_eq!(classify(0b0101_0000), Some(Pattern::Edge));
+    }
+}

crates/mesh/src/octree/tests.rs

@@ -0,0 +1,98 @@
+use super::*;
+use crate::octree::geom::tests::box_surface;
+use crate::size::UniformSize;
+
+#[test]
+fn a_box_meshes_with_exact_volume_and_tagged_faces() {
+    let s = box_surface([0.0; 3], [10.0, 8.0, 6.0]);
+    let m = octree_hex_mesh(&s, &UniformSize(2.0), &OctreeOptions::default()).unwrap();
+    let r = &m.report;
+    assert!((r.volume - 480.0).abs() < 1e-9 * 480.0, "{r:?}");
+    assert!(r.min_scaled_jacobian > 0.3, "{r:?}");
+    let faces: BTreeSet<u32> = m.boundary.iter().map(|b| b.brep_face).collect();
+    assert_eq!(faces.len(), 6);
+    // Each boundary quad lies in the plane of its face.
+    for b in &m.boundary {
+        let hx = m.hexes[b.element as usize];
+        let q = topo::FACES[(b.face - 1) as usize].map(|i| m.nodes[hx[i] as usize]);
+        let (axis, value) = match b.brep_face {
+            0 => (2, 0.0),
+            1 => (2, 6.0),
+            2 => (1, 0.0),
+            3 => (1, 8.0),
+            4 => (0, 0.0),
+            _ => (0, 10.0),
+        };
+        assert!(q.iter().all(|p| (p[axis] - value).abs() < 1e-9), "{b:?} {q:?}");
+    }
+}
+
+/// A closed cylinder of radius `r` and height `h` on the z axis, `n` segments: faces 0 (bottom
+/// cap), 1 (top cap), 2 (side).
+pub(crate) fn cylinder_surface(r: f64, h: f64, n: usize) -> BoundarySurface {
+    let mut p = vec![[0.0, 0.0, 0.0], [0.0, 0.0, h]];
+    for z in [0.0, h] {
+        for i in 0..n {
+            let a = std::f64::consts::TAU * i as f64 / n as f64;
+            p.push([r * a.cos(), r * a.sin(), z]);
+        }
+    }
+    let (b, t) = (|i: usize| (2 + i % n) as u32, |i: usize| (2 + n + i % n) as u32);
+    let mut tris = Vec::new();
+    let mut face = Vec::new();
+    for i in 0..n {
+        tris.push([0, b(i + 1), b(i)]);
+        face.push(0);
+        tris.push([1, t(i), t(i + 1)]);
+        face.push(1);
+        tris.push([b(i), b(i + 1), t(i + 1)]);
+        tris.push([b(i), t(i + 1), t(i)]);
+        face.extend([2, 2]);
+    }
+    BoundarySurface::new(&p, &tris, &face, &crate::surface::SurfaceOptions::default()).unwrap()
+}
+
+/// A UV sphere of radius `r`, one face.
+pub(crate) fn sphere_surface(r: f64, n: usize) -> BoundarySurface {
+    let rings = n / 2;
+    let mut p = vec![[0.0, 0.0, -r], [0.0, 0.0, r]];
+    for j in 1..rings {
+        let th = std::f64::consts::PI * j as f64 / rings as f64 - std::f64::consts::FRAC_PI_2;
+        for i in 0..n {
+            let a = std::f64::consts::TAU * i as f64 / n as f64;
+            p.push([r * th.cos() * a.cos(), r * th.cos() * a.sin(), r * th.sin()]);
+        }
+    }
+    let id = |j: usize, i: usize| (2 + (j - 1) * n + i % n) as u32;
+    let mut tris = Vec::new();
+    for i in 0..n {
+        tris.push([0, id(1, i + 1), id(1, i)]);
+        tris.push([1, id(rings - 1, i), id(rings - 1, i + 1)]);
+    }
+    for j in 1..rings - 1 {
+        for i in 0..n {
+            tris.push([id(j, i), id(j, i + 1), id(j + 1, i + 1)]);
+            tris.push([id(j, i), id(j + 1, i + 1), id(j + 1, i)]);
+        }
+    }
+    let face = vec![0; tris.len()];
+    BoundarySurface::new(&p, &tris, &face, &crate::surface::SurfaceOptions::default()).unwrap()
+}
+
+#[test]
+fn a_cylinder_meshes_with_its_rims_captured() {
+    let s = cylinder_surface(5.0, 8.0, 64);
+    let m = octree_hex_mesh(&s, &UniformSize(2.0), &OctreeOptions::default()).unwrap();
+    let r = &m.report;
+    assert!(((r.volume - r.surface_volume) / r.surface_volume).abs() < 0.01, "{r:?}");
+    let faces: BTreeSet<u32> = m.boundary.iter().map(|b| b.brep_face).collect();
+    assert_eq!(faces.len(), 3);
+}
+
+#[test]
+fn a_sphere_meshes() {
+    let s = sphere_surface(5.0, 48);
+    let m = octree_hex_mesh(&s, &UniformSize(2.0), &OctreeOptions::default()).unwrap();
+    let r = &m.report;
+    assert!(((r.volume - r.surface_volume) / r.surface_volume).abs() < 0.01, "{r:?}");
+}

crates/mesh/src/octree/topo.rs

@@ -0,0 +1,404 @@
+//! Hexahedral mesh topology for the boundary stages: boundary faces, the manifold clean-up of
+//! the core, padding a boundary with a layer of hexahedra, and pillowing a shrink set.
+//!
+//! *Padding* puts one hexahedron on every boundary quad, between the quad and a copy of it, so
+//! no hexahedron keeps more than one face on the boundary. *Pillowing* a shrink set `S` inserts
+//! a layer of hexahedra along the faces between `S` and the rest of the mesh; the layer reaches
+//! the boundary along the loop where `S`'s boundary faces meet the others' (S. A. Mitchell,
+//! T. J. Tautges, *Pillowing doublets: refining a mesh to ensure that faces share at most one
+//! edge*, IMR 1995; Gao, Shen, Panozzo 2019, §3.4 use the same insertion around features).
+
+use std::collections::{BTreeMap, HashMap};
+
+use crate::{Error, Result};
+
+/// The six faces of a hex8 as local node positions, right-hand normals outward (the model's
+/// `ElementKind::Hex8.faces()` order).
+pub(crate) const FACES: [[usize; 4]; 6] = [[0, 3, 2, 1], [4, 5, 6, 7], [0, 1, 5, 4], [1, 2, 6, 5], [2, 3, 7, 6], [3, 0, 4, 7]];
+
+/// A face key: the vertices sorted.
+pub(crate) fn key(q: [u32; 4]) -> [u32; 4] {
+    let mut k = q;
+    k.sort_unstable();
+    k
+}
+
+/// The root of `i` in the union-find forest `c`, compressing the path from `i`.
+pub(crate) fn find(c: &mut [usize], i: usize) -> usize {
+    let mut r = i;
+    while let Some(&p) = c.get(r).filter(|&&p| p != r) {
+        r = p;
+    }
+    if let Some(x) = c.get_mut(i) {
+        *x = r;
+    }
+    r
+}
+
+/// A face of an element: (hexahedron, local face 0..6).
+pub(crate) type FaceRef = (u32, u8);
+
+/// The vertices of a local face, outward.
+pub(crate) fn face(h: &[u32; 8], f: u8) -> [u32; 4] {
+    FACES.get(f as usize).map_or([0; 4], |l| l.map(|i| h[i]))
+}
+
+/// Every face's elements, keyed by sorted vertices (deterministic map).
+pub(crate) fn face_map(hexes: &[[u32; 8]]) -> BTreeMap<[u32; 4], Vec<FaceRef>> {
+    let mut m: BTreeMap<[u32; 4], Vec<FaceRef>> = BTreeMap::new();
+    for (i, h) in hexes.iter().enumerate() {
+        for f in 0..6u8 {
+            m.entry(key(face(h, f))).or_default().push((i as u32, f));
+        }
+    }
+    m
+}
+
+/// The boundary faces (on one element only), outward, in face-map order.
+pub(crate) fn boundary(hexes: &[[u32; 8]]) -> Vec<(FaceRef, [u32; 4])> {
+    face_map(hexes)
+        .into_values()
+        .filter(|v| v.len() == 1)
+        .filter_map(|v| v.first().copied())
+        .map(|(h, f)| ((h, f), hexes.get(h as usize).map_or([0; 4], |x| face(x, f))))
+        .collect()
+}
+
+/// Edges of a quad list with the quads using each (ascending vertex pair → quad indices).
+pub(crate) fn quad_edges(quads: &[[u32; 4]]) -> BTreeMap<(u32, u32), Vec<u32>> {
+    let mut m: BTreeMap<(u32, u32), Vec<u32>> = BTreeMap::new();
+    for (i, q) in quads.iter().enumerate() {
+        for k in 0..4 {
+            let (a, b) = (q[k], q[(k + 1) % 4]);
+            m.entry((a.min(b), a.max(b))).or_default().push(i as u32);
+        }
+    }
+    m
+}
+
+/// The hexahedra to drop so that the boundary of `keep` is a 2-manifold: elements at edges with
+/// four boundary faces (elements touching along an edge only) and at vertices whose boundary
+/// faces form more than one fan. Returns the indices to drop (empty when manifold).
+pub(crate) fn nonmanifold(hexes: &[[u32; 8]], keep: &[bool]) -> Vec<u32> {
+    let kept: Vec<u32> = (0..hexes.len() as u32).filter(|&i| keep.get(i as usize).copied().unwrap_or(false)).collect();
+    let sub: Vec<[u32; 8]> = kept.iter().filter_map(|&i| hexes.get(i as usize).copied()).collect();
+    let b = boundary(&sub);
+    let quads: Vec<[u32; 4]> = b.iter().map(|(_, q)| *q).collect();
+    let mut bad_vertices: Vec<u32> = Vec::new();
+    for ((a, c), qs) in quad_edges(&quads) {
+        if qs.len() != 2 {
+            bad_vertices.extend([a, c]);
+        }
+    }
+    // Fans: around each vertex, boundary quads connected through shared edges at the vertex.
+    let mut at: BTreeMap<u32, Vec<u32>> = BTreeMap::new();
+    for (i, q) in quads.iter().enumerate() {
+        for &v in q {
+            at.entry(v).or_default().push(i as u32);
+        }
+    }
+    for (v, qs) in &at {
+        let mut comp: Vec<usize> = (0..qs.len()).collect();
+        for i in 0..qs.len() {
+            for j in i + 1..qs.len() {
+                let (qa, qb) = (quads[qs[i] as usize], quads[qs[j] as usize]);
+                // Share an edge through v: another common vertex adjacent to v in both.
+                let nb = |q: [u32; 4]| {
+                    let p = q.iter().position(|&x| x == *v).unwrap_or(0);
+                    [q[(p + 1) % 4], q[(p + 3) % 4]]
+                };
+                if nb(qa).iter().any(|x| nb(qb).contains(x)) {
+                    let (ri, rj) = (find(&mut comp, i), find(&mut comp, j));
+                    comp[ri] = rj;
+                }
+            }
+        }
+        let roots = (0..qs.len()).filter(|&i| find(&mut comp, i) == i).count();
+        if roots > 1 {
+            bad_vertices.push(*v);
+        }
+    }
+    bad_vertices.sort_unstable();
+    bad_vertices.dedup();
+    if bad_vertices.is_empty() {
+        return Vec::new();
+    }
+    kept.into_iter().filter(|&i| hexes.get(i as usize).is_some_and(|h| h.iter().any(|v| bad_vertices.binary_search(v).is_ok()))).collect()
+}
+
+/// Pad: a hexahedron on each listed boundary face (outward), between it and a copy; the copies
+/// of the faces' vertices are new points (appended, initially at the original positions).
+/// Returns the copy of each original vertex.
+pub(crate) fn pad(points: &mut Vec<[f64; 3]>, hexes: &mut Vec<[u32; 8]>, faces: &[[u32; 4]]) -> Result<BTreeMap<u32, u32>> {
+    let mut copy: BTreeMap<u32, u32> = BTreeMap::new();
+    for q in faces {
+        for &v in q {
+            if let std::collections::btree_map::Entry::Vacant(slot) = copy.entry(v) {
+                let n = u32::try_from(points.len()).map_err(|_| Error::Limit("too many points".into()))?;
+                points.push(points.get(v as usize).copied().unwrap_or([0.0; 3]));
+                slot.insert(n);
+            }
+        }
+    }
+    for q in faces {
+        let c = q.map(|v| copy.get(&v).copied().unwrap_or(v));
+        hexes.push([q[0], q[1], q[2], q[3], c[0], c[1], c[2], c[3]]);
+    }
+    Ok(copy)
+}
+
+/// One vertex copy made by [`pillow`]: the original, the copy, and the shrink-set elements
+/// around the original that now use the copy.
+#[derive(Debug, Clone)]
+pub(crate) struct PillowCopy {
+    pub original: u32,
+    pub copy: u32,
+    pub elements: Vec<u32>,
+}
+
+/// Pillow the shrink set `inside` (one flag per element): a new layer of hexahedra on every face
+/// between an element of the set and one outside it. Vertices of those faces get one copy per
+/// group of set elements around them (connected through faces inside the set), which the set's
+/// elements now use; the new layer joins the originals to the copies. Faces on the boundary are
+/// not pillowed: where the layer meets the boundary, the set's boundary faces move to the copies,
+/// so the boundary gains a ring of quads along the loop. New points are appended at the
+/// originals' positions. Returns the copies made.
+pub(crate) fn pillow(points: &mut Vec<[f64; 3]>, hexes: &mut Vec<[u32; 8]>, inside: &[bool]) -> Result<Vec<PillowCopy>> {
+    let fm = face_map(hexes);
+    let in_set = |h: u32| inside.get(h as usize).copied().unwrap_or(false);
+    // Faces to pillow, oriented from the outside element into the set.
+    let mut sheet: Vec<(u32, [u32; 4])> = Vec::new();
+    for refs in fm.values() {
+        if let [(a, fa), (b, fb)] = refs.as_slice() {
+            match (in_set(*a), in_set(*b)) {
+                (true, false) => sheet.push((*a, hexes.get(*a as usize).map_or([0; 4], |h| face(h, *fa)))),
+                (false, true) => sheet.push((*b, hexes.get(*b as usize).map_or([0; 4], |h| face(h, *fb)))),
+                _ => {}
+            }
+        }
+    }
+    if sheet.is_empty() {
+        return Ok(Vec::new());
+    }
+    let mut on_sheet: Vec<u32> = sheet.iter().flat_map(|(_, q)| *q).collect();
+    on_sheet.sort_unstable();
+    on_sheet.dedup();
+    // Set elements around each sheet vertex.
+    let mut around: HashMap<u32, Vec<u32>> = HashMap::new();
+    for (i, h) in hexes.iter().enumerate() {
+        if !in_set(i as u32) {
+            continue;
+        }
+        for &v in h {
+            if on_sheet.binary_search(&v).is_ok() {
+                let e = around.entry(v).or_default();
+                if e.last() != Some(&(i as u32)) {
+                    e.push(i as u32);
+                }
+            }
+        }
+    }
+    let mut copies: Vec<PillowCopy> = Vec::new();
+    // (vertex, element) → copy.
+    let mut copy_of: HashMap<(u32, u32), u32> = HashMap::new();
+    for &v in &on_sheet {
+        let mut els = around.get(&v).cloned().unwrap_or_default();
+        els.sort_unstable();
+        els.dedup();
+        // Group by faces inside the set that contain v.
+        let mut group: Vec<usize> = (0..els.len()).collect();
+        for i in 0..els.len() {
+            for j in i + 1..els.len() {
+                let (Some(hi), Some(hj)) = (hexes.get(els[i] as usize), hexes.get(els[j] as usize)) else { continue };
+                let shares = (0..6u8).any(|f| {
+                    let q = face(hi, f);
+                    q.contains(&v) && (0..6u8).any(|g| key(face(hj, g)) == key(q))
+                });
+                if shares {
+                    let (a, b) = (find(&mut group, i), find(&mut group, j));
+                    group[a] = b;
+                }
+            }
+        }
+        // Per group: its copy and that copy's index in `copies`.
+        let mut made: BTreeMap<usize, (u32, usize)> = BTreeMap::new();
+        for (i, &e) in els.iter().enumerate() {
+            let r = find(&mut group, i);
+            let (c, at) = match made.get(&r) {
+                Some(&m) => m,
+                None => {
+                    let c = u32::try_from(points.len()).map_err(|_| Error::Limit("too many points".into()))?;
+                    points.push(points.get(v as usize).copied().unwrap_or([0.0; 3]));
+                    copies.push(PillowCopy { original: v, copy: c, elements: Vec::new() });
+                    made.insert(r, (c, copies.len() - 1));
+                    (c, copies.len() - 1)
+                }
+            };
+            copy_of.insert((v, e), c);
+            if let Some(pc) = copies.get_mut(at) {
+                pc.elements.push(e);
+            }
+        }
+    }
+    // The new layer: original face (seen from outside, normal into the set) below, copies above.
+    let mut layer: Vec<[u32; 8]> = Vec::with_capacity(sheet.len());
+    for (h, q) in &sheet {
+        let r = [q[0], q[3], q[2], q[1]];
+        let c = r.map(|v| copy_of.get(&(v, *h)).copied().unwrap_or(v));
+        layer.push([r[0], r[1], r[2], r[3], c[0], c[1], c[2], c[3]]);
+    }
+    for (i, h) in hexes.iter_mut().enumerate() {
+        if in_set(i as u32) {
+            for v in h.iter_mut() {
+                if let Some(&c) = copy_of.get(&(*v, i as u32)) {
+                    *v = c;
+                }
+            }
+        }
+    }
+    hexes.extend(layer);
+    Ok(copies)
+}
+
+/// Hexahedra to drop so that every vertex of the boundary of `keep` can be pushed outward: the
+/// normals of its boundary quads must lie strictly inside one half-space (their mean has a
+/// positive product with each, at least `margin`). Where this fails, padding from the vertex
+/// to the part's surface inverts elements, so for each such vertex the kept element at it with
+/// the most boundary faces is dropped (the buffer-zone rule of H. Tong, E. Halilaj, Y. J. Zhang,
+/// *HybridOctree_Hex*, arXiv 2401.05984, 2024, §2.3, in a mean-normal form).
+pub(crate) fn hidden_vertices(hexes: &[[u32; 8]], keep: &[bool], xyz: &dyn Fn(u32) -> [f64; 3], margin: f64) -> Vec<u32> {
+    let kept: Vec<u32> = (0..hexes.len() as u32).filter(|&i| keep.get(i as usize).copied().unwrap_or(false)).collect();
+    let sub: Vec<[u32; 8]> = kept.iter().filter_map(|&i| hexes.get(i as usize).copied()).collect();
+    let b = boundary(&sub);
+    let mut normals: BTreeMap<u32, Vec<[f64; 3]>> = BTreeMap::new();
+    let mut faces_of: BTreeMap<u32, usize> = BTreeMap::new();
+    for ((h, _), q) in &b {
+        let p = q.map(xyz);
+        let d1 = [0, 1, 2].map(|k| p[2][k] - p[0][k]);
+        let d2 = [0, 1, 2].map(|k| p[3][k] - p[1][k]);
+        let n = [d1[1] * d2[2] - d1[2] * d2[1], d1[2] * d2[0] - d1[0] * d2[2], d1[0] * d2[1] - d1[1] * d2[0]];
+        let l = (n[0] * n[0] + n[1] * n[1] + n[2] * n[2]).sqrt();
+        if l > 0.0 {
+            for &v in q {
+                normals.entry(v).or_default().push(n.map(|c| c / l));
+            }
+        }
+        *faces_of.entry(*h).or_default() += 1;
+    }
+    let mut drop: Vec<u32> = Vec::new();
+    for (v, ns) in &normals {
+        let mut m = [0.0; 3];
+        for n in ns {
+            for k in 0..3 {
+                m[k] += n[k];
+            }
+        }
+        let l = (m[0] * m[0] + m[1] * m[1] + m[2] * m[2]).sqrt();
+        let ok = l > 0.0 && ns.iter().all(|n| (n[0] * m[0] + n[1] * m[1] + n[2] * m[2]) / l > margin);
+        if ok {
+            continue;
+        }
+        // The kept element at v with the most boundary faces (lowest index on ties).
+        let best = sub
+            .iter()
+            .enumerate()
+            .filter(|(_, h)| h.contains(v))
+            .max_by(|a, b| faces_of.get(&(a.0 as u32)).unwrap_or(&0).cmp(faces_of.get(&(b.0 as u32)).unwrap_or(&0)).then(b.0.cmp(&a.0)))
+            .and_then(|(i, _)| kept.get(i).copied());
+        if let Some(e) = best {
+            drop.push(e);
+        }
+    }
+    drop.sort_unstable();
+    drop.dedup();
+    drop
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn faces_follow_the_model_hex8_face_order() {
+        let model = fieldcraft_model::ElementKind::Hex8.faces();
+        assert_eq!(model.len(), FACES.len());
+        for (m, f) in model.iter().zip(FACES) {
+            assert_eq!(*m, &f[..]);
+        }
+    }
+
+    /// An `n × n × n` block of unit cubes.
+    pub(crate) fn block(n: u32) -> (Vec<[f64; 3]>, Vec<[u32; 8]>) {
+        let id = |x: u32, y: u32, z: u32| x + (n + 1) * (y + (n + 1) * z);
+        let mut p = Vec::new();
+        for z in 0..=n {
+            for y in 0..=n {
+                for x in 0..=n {
+                    p.push([x as f64, y as f64, z as f64]);
+                }
+            }
+        }
+        let mut h = Vec::new();
+        for z in 0..n {
+            for y in 0..n {
+                for x in 0..n {
+                    h.push([
+                        id(x, y, z),
+                        id(x + 1, y, z),
+                        id(x + 1, y + 1, z),
+                        id(x, y + 1, z),
+                        id(x, y, z + 1),
+                        id(x + 1, y, z + 1),
+                        id(x + 1, y + 1, z + 1),
+                        id(x, y + 1, z + 1),
+                    ]);
+                }
+            }
+        }
+        (p, h)
+    }
+
+    fn closed(hexes: &[[u32; 8]]) {
+        for (k, v) in face_map(hexes) {
+            assert!(v.len() <= 2, "face {k:?} on {} elements", v.len());
+        }
+        let b = boundary(hexes);
+        for (_, qs) in quad_edges(&b.iter().map(|x| x.1).collect::<Vec<_>>()) {
+            assert_eq!(qs.len(), 2);
+        }
+    }
+
+    #[test]
+    fn padding_and_pillowing_keep_a_closed_conforming_mesh() {
+        let (mut p, mut h) = block(3);
+        let b: Vec<[u32; 4]> = boundary(&h).into_iter().map(|x| x.1).collect();
+        assert_eq!(b.len(), 54);
+        let copy = pad(&mut p, &mut h, &b).unwrap();
+        assert_eq!(copy.len(), 56);
+        assert_eq!(h.len(), 27 + 54);
+        closed(&h);
+        // Pillow the padding hexes on the face z = 0 (outer faces down).
+        let pad_start = 27;
+        let inside: Vec<bool> = (0..h.len()).map(|i| i >= pad_start && (0..4).all(|k| p[h[i][k] as usize][2] == 0.0)).collect();
+        assert_eq!(inside.iter().filter(|x| **x).count(), 9);
+        let before = h.len();
+        let copies = pillow(&mut p, &mut h, &inside).unwrap();
+        closed(&h);
+        // Sheet: the 9 inner faces plus the 12 side faces to the neighbouring padding.
+        assert_eq!(h.len(), before + 9 + 12);
+        // The boundary gains the ring: 54 + 12 quads.
+        assert_eq!(boundary(&h).len(), 54 + 12);
+        assert!(!copies.is_empty());
+    }
+
+    #[test]
+    fn nonmanifold_edges_and_vertices_are_found() {
+        let (_, h) = block(2);
+        // Two cubes touching along an edge.
+        let keep = |ids: &[usize]| (0..8).map(|i| ids.contains(&i)).collect::<Vec<bool>>();
+        assert!(!nonmanifold(&h, &keep(&[0, 3])).is_empty(), "edge contact");
+        assert!(!nonmanifold(&h, &keep(&[0, 7])).is_empty(), "vertex contact");
+        assert!(nonmanifold(&h, &keep(&[0, 1, 3])).is_empty());
+        assert!(nonmanifold(&h, &keep(&[0, 1, 2, 3, 4, 5, 6, 7])).is_empty());
+    }
+}

crates/testkit/src/runner.rs

@@ -16,6 +16,7 @@ use crate::corpus::{self, CATEGORIES, Part};
 use fieldcraft_geometry::{Design, FaceSurface, Measure, SurfaceMesh};
 use fieldcraft_mesh::cdt::{CdtOptions, tetrahedralize};
 use fieldcraft_mesh::hex::{Cell, HexOptions, Surface, hex_mesh, hex_mesh_cells};
+use fieldcraft_mesh::octree::{OctreeOptions, octree_hex_mesh};
 use fieldcraft_mesh::size::UniformSize;
 use fieldcraft_mesh::surface::{BoundarySurface, SurfaceOptions};
 
@@ -33,6 +34,8 @@ pub struct Loaded<'a> {
     pub surfaces: Vec<SurfaceMesh>,
     /// The chordal tolerance the surfaces were made with (mm).
     pub tolerance: f64,
+    /// The bounding-box diagonal of all bodies (mm).
+    pub diagonal: f64,
 }
 
 /// A column group of the results table: a mesher, a check.
@@ -44,7 +47,7 @@ pub trait Stage {
 
 /// The stages, in column order. Meshers add theirs here.
 fn stages() -> Vec<Box<dyn Stage>> {
-    vec![Box::new(SurfaceVolume), Box::new(BoundaryCheck), Box::new(TetCdt), Box::new(StructuredHex)]
+    vec![Box::new(SurfaceVolume), Box::new(BoundaryCheck), Box::new(TetCdt), Box::new(StructuredHex), Box::new(OctreeHex)]
 }
 
 /// The boundary surface's enclosed volume against SolidCraft's measured volume and the closed
@@ -116,7 +119,7 @@ impl Stage for TetCdt {
     }
 }
 
-/// Element size of the structured hex stage, relative to the part's bounding-box diagonal.
+/// Element size of the hex stages, relative to the part's bounding-box diagonal.
 const HEX_RELATIVE_SIZE: f64 = 1.0 / 25.0;
 /// Characters of a "not meshed" reason kept in the table.
 const REASON_CHARS: usize = 240;
@@ -134,15 +137,7 @@ impl Stage for StructuredHex {
 
     fn run(&self, l: &Loaded) -> Result<Vec<String>, String> {
         let start = Instant::now();
-        let (lo, hi) = l.measures.iter().fold(([f64::INFINITY; 3], [f64::NEG_INFINITY; 3]), |(mut lo, mut hi), m| {
-            for i in 0..3 {
-                lo[i] = lo[i].min(m.bbox_min[i]);
-                hi[i] = hi[i].max(m.bbox_max[i]);
-            }
-            (lo, hi)
-        });
-        let diagonal = (0..3).map(|i| (hi[i] - lo[i]).powi(2)).sum::<f64>().sqrt();
-        let size = UniformSize(HEX_RELATIVE_SIZE * diagonal);
+        let size = UniformSize(HEX_RELATIVE_SIZE * l.diagonal);
         let linear = HexOptions { chordal_tolerance: l.tolerance, ..HexOptions::default() };
         let quadratic = HexOptions { quadratic: true, ..linear.clone() };
         if l.part.categories.iter().any(|c| c == "partitioned") {
@@ -231,6 +226,49 @@ pub fn analytic(s: FaceSurface) -> Surface {
     }
 }
 
+/// Element cap of the octree hex column: well above every part that meshes, and low enough that a
+/// part whose grid explodes (a thin spline wall refined everywhere) is refused in seconds instead
+/// of optimised for an hour.
+const OCTREE_MAX_ELEMENTS: usize = 500_000;
+
+/// The unstructured all-hex mesher (`fieldcraft_mesh::octree`) at a size of 1/25 of the part's
+/// diagonal: elements, the smallest and mean scaled Jacobian, and the mesh volume against
+/// SolidCraft's. Every element is proven valid before the mesher returns; a part it cannot mesh
+/// gets the reason, as in the structured hex column.
+struct OctreeHex;
+
+impl Stage for OctreeHex {
+    fn headers(&self) -> &'static [&'static str] {
+        &["octree hexes", "min SJ", "mean SJ", "hex vol. diff", "octree ms"]
+    }
+
+    fn run(&self, l: &Loaded) -> Result<Vec<String>, String> {
+        let start = Instant::now();
+        let (mut hexes, mut min_sj, mut sum_sj, mut volume) = (0usize, f64::INFINITY, 0.0, 0.0);
+        for (body, mesh) in l.design.bodies.iter().zip(&l.surfaces) {
+            let options = SurfaceOptions { weld_distance: mesh.weld_distance, ..SurfaceOptions::default() };
+            let s = BoundarySurface::new(&mesh.positions, &mesh.triangles, &mesh.tri_face, &options).map_err(|e| format!("{}: {e}", body.name))?;
+            let options = OctreeOptions { max_elements: OCTREE_MAX_ELEMENTS, ..OctreeOptions::default() };
+            let m = match octree_hex_mesh(&s, &UniformSize(HEX_RELATIVE_SIZE * l.diagonal), &options) {
+                Ok(m) => m,
+                Err(e) => {
+                    let reason: String = e.to_string().chars().take(REASON_CHARS).collect();
+                    let ms = start.elapsed().as_secs_f64() * 1e3;
+                    return Ok(vec![format!("no: {}: {reason}", body.name), String::new(), String::new(), String::new(), format!("{ms:.0}")]);
+                }
+            };
+            hexes += m.hexes.len();
+            min_sj = min_sj.min(m.report.min_scaled_jacobian);
+            sum_sj += m.report.mean_scaled_jacobian * m.hexes.len() as f64;
+            volume += m.report.volume;
+        }
+        let measured: f64 = l.measures.iter().map(|m| m.volume).sum();
+        let rel = if measured.abs() > 0.0 { (volume - measured) / measured } else { f64::NAN };
+        let ms = start.elapsed().as_secs_f64() * 1e3;
+        Ok(vec![hexes.to_string(), format!("{min_sj:.3}"), format!("{:.3}", sum_sj / hexes.max(1) as f64), format!("{rel:.1e}"), format!("{ms:.0}")])
+    }
+}
+
 fn load(part: &Part) -> Result<(Loaded<'_>, f64), String> {
     let start = Instant::now();
     let design = corpus::load(part).map_err(|e| e.to_string())?;
@@ -250,7 +288,7 @@ fn load(part: &Part) -> Result<(Loaded<'_>, f64), String> {
     let surfaces =
         design.bodies.iter().map(|b| b.tessellate(tolerance).map_err(|e| format!("tessellate {}: {e}", b.name))).collect::<Result<Vec<_>, _>>()?;
     let ms = start.elapsed().as_secs_f64() * 1e3;
-    Ok((Loaded { part, design, measures, surfaces, tolerance }, ms))
+    Ok((Loaded { part, design, measures, surfaces, tolerance, diagonal }, ms))
 }
 
 /// A failed part or stage: the error, then blanks to fill `width` cells.

docs/architecture.md

@@ -68,12 +68,21 @@ Rules beyond the layer order:
    its B-rep face tag survives. Then constrained Delaunay refinement to a
    size field, ODT smoothing interleaved as in Tournois et al. 2009, and the Klingner–Shewchuk 2007
    improvement schedule against slivers.
-4. **Hexahedral meshing (the core).** Octree/grid dual all-hex first, the one family that meshes
-   any input validly (Livesu et al. 2021 dual schemes, Pitzalis et al. 2021 refinement, Gao et al.
-   2019 feature projection and padding, Garanzha et al. 2021 untangling). Then hex-dominant by tet
+4. **Hexahedral meshing (the core).** Grid-based all-hex first, the one family that meshes any
+   input validly. `fieldcraft-mesh::octree` (command `mesh.hex.octree`): a 27-tree over the part,
+   refined by curvature, thickness and narrow faces, balanced and cut by the 3-refinement
+   node-marking templates (Schneiders 1996, Ito et al. 2009) into a conforming grid; the core of
+   grid elements inside the part (exact inside test) is cleaned to a manifold boundary; the part's
+   patches, sharp curves and corners are mapped onto it (Gao et al. 2019); one padding layer and a
+   pillow under each patch carry the boundary and capture every sharp curve; the foldover-free
+   energy (Garanzha et al. 2021) with a growing fitting term moves the boundary onto the geometry
+   while every element stays valid; the worst elements are improved, and every element is proven
+   valid by the certified Bernstein Jacobian check (Johnen et al. 2013; Marschner et al. 2020)
+   before the boundary quads are tagged with their B-rep faces. Where a stage fails the grid is
+   refined there and the mesh rebuilt. The dual schemes of Livesu et al. 2021 (fewer elements,
+   concave transitions) are the next upgrade of the transitions. Then hex-dominant by tet
    recombination (Pellerin et al. 2018), then polycube (DualCube labelling) for blocky parts;
-   frame-field/integer-grid-map methods stay research. Validity uses a certified trilinear Jacobian
-   check (Marschner et al. 2020), not corner Jacobians alone.
+   frame-field/integer-grid-map methods stay research.
 5. **Structured hex and face meshing.** Faces get all-quad meshes on shared, prescribed edge
    nodes (`fieldcraft-mesh::quad`): transfinite interpolation for four-sided faces (Gordon and
    Hall 1973), and for free faces a half-resolution constrained Delaunay triangulation,

docs/features.md

@@ -6,7 +6,7 @@ Every entry is a toolbar button. **Ready** entries run their command; **coming s
 shown greyed with the roadmap milestone that brings them (ROADMAP.md). Scripts, the command line,
 the control channel and AI agents over MCP (docs/mcp.md) run the same commands.
 
-55 ready, 49 coming soon.
+56 ready, 48 coming soon.
 
 ## Geometry
 
@@ -38,7 +38,7 @@ the control channel and AI agents over MCP (docs/mcp.md) run the same commands.
 | Structured Hex | Block Hex — One conforming hex mesh of a part split into six-sided blocks: shared faces meshed once | ready | `mesh.hex.blocks` |
 | Structured Hex | Decompose — Find sweepable and mappable regions and propose the splits | coming soon (M5) |  |
 | Structured Hex | Transition — Conforming tet–hex interfaces with pyramids | coming soon (M5) |  |
-| Automatic Hex | All-Hex — Octree dual all-hex mesh of any part | coming soon (M6) |  |
+| Automatic Hex | Octree Hex — Unstructured all-hex mesh of any part: adaptive grid, boundary layers on every sharp edge, untangled and proven valid | ready | `mesh.hex.octree` |
 | Automatic Hex | Hex-Dominant — Hex-dominant mesh by recombining tets | coming soon (M6) |  |
 | Automatic Hex | Polycube — Polycube hex mesh of blocky parts | coming soon (M6) |  |
 | Improve | Smooth — Smooth node positions, keeping every element valid | coming soon (M6) |  |

docs/meshing/corpus.md

@@ -8,31 +8,31 @@ clock on the machine that ran it; every other column is deterministic.
 
 ## Results
 
-| part | categories | bodies | faces | triangles | tol. mm | build ms | enclosed mm³ | SolidCraft mm³ | rel. diff | exact mm³ | boundary | genus | feature edges | check ms | CDT tets | Steiner | cavities | tet gate | CDT ms | hex | hex elements | hex8 min SJ | hex20 min SJ | hex ms |
-|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
-| bent-tube | sweepable, thin-wall, genus | 1 | 18 | 4296 | 0.121 | 428 | 4213.1 | 4159.4 | 1.3e-2 | 4159.8 | ok, 2148 vertices | 1 | 268 | 11.8 | 15550 | 2338 | 0 | pass | 598 | sweep | 192 | 0.859 | 0.997 | 62 |
-| bolted-assembly | assembly, blocky | 3 | 22 | 284 | 0.117 | 53 | 33004.7 | 32999.0 | 1.7e-4 | 32998.9 | ok, 144 vertices | 0, 1, 1 | 104 | 0.7 | 690 | 181 | 0 | pass | 33 | extrude | 357 | 0.398 | 0.222 | 40 |
-| cantilever-i-beam | sweepable, thin-wall | 1 | 14 | 44 | 0.208 | 4 | 61600.0 | 61600.0 | 0.0e0 | 61600.0 | ok, 24 vertices | 0 | 36 | 0.1 | 90 | 24 | 0 | pass | 3 | extrude | 336 | 0.707 | 0.707 | 33 |
-| cross-cylinders | boolean, sharp | 1 | 10 | 1395 | 0.087 | 10640 | 32197.4 | 32357.8 | -5.0e-3 | 32365.8 | error: Cross: invalid boundary surface: 28 defects; triangles 142 and 184 are oriented inconsistently across edge 0–396 at (-15.488656021315903, 4.577998041166734, 7.853092194999726); triangles 184 and 239 are oriented inconsistently across edge 0–554 at (-18.75000000000002, 5.634596901641276, 7.632255089815496); edge 103–108 at (-0.16007784063084018, 0.03984236936862052, 9.998639270596865) is shared by 3 triangles; … |  |  |  | error: Cross: invalid boundary surface: 28 defects; triangles 142 and 184 are oriented inconsistently across edge 0–396 at (-15.488656021315903, 4.577998041166734, 7.853092194999726); triangles 184 and 239 are oriented inconsistently across edge 0–554 at (-18.75000000000002, 5.634596901641276, 7.632255089815496); edge 103–108 at (-0.16007784063084018, 0.03984236936862052, 9.998639270596865) is shared by 3 triangles; … |  |  |  |  | error: Cross: invalid boundary surface: 28 defects; triangles 142 and 184 are oriented inconsistently across edge 0–396 at (-15.488656021315903, 4.577998041166734, 7.853092194999726); triangles 184 and 239 are oriented inconsistently across edge 0–554 at (-18.75000000000002, 5.634596901641276, 7.632255089815496); edge 103–108 at (-0.16007784063084018, 0.03984236936862052, 9.998639270596865) is shared by 3 triangles; … |  |  |  |  |
-| filleted-bracket | fillet, genus, sweepable | 1 | 13 | 540 | 0.088 | 494 | 32162.3 | 32145.0 | 5.4e-4 | 32144.8 | ok, 268 vertices | 2 | 135 | 1.7 | 5582 | 1613 | 0 | pass | 233 | no: Bracket: not supported yet: no structured hex method fits this body (extrude: cylindrical face 8 is not parallel to the direction; sweep: from face 1: linking face 6 has 2 boundary loops (periodic or holed linking faces are not swept); from face 6:  |  |  |  | 5 |
-| filleted-bracket-blocks | fillet, genus, partitioned | 30 | 190 | 712 | 0.088 | 1057 | 32169.8 | 32145.1 | 7.7e-4 | 32144.8 | ok, 416 vertices | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 | 474 | 1.7 | 1423 | 349 | 0 | pass | 64 | blocks (30 cells, 45 shared faces) | 4986 | 0.705 | 0.705 | 1440 |
-| i-beam | sweepable, thin-wall | 1 | 14 | 44 | 0.208 | 234 | 61600.0 | 61600.0 | 0.0e0 | 61600.0 | ok, 24 vertices | 0 | 36 | 0.3 | 94 | 24 | 0 | pass | 8 | extrude | 336 | 0.707 | 0.707 | 74 |
-| loft-transition | organic | 1 | 6 | 292 | 0.075 | 493 | 56390.2 | 56450.6 | -1.1e-3 | – | ok, 148 vertices | 0 | 40 | 1.5 | 1329 | 256 | 1 | pass | 80 | sweep | 2873 | 0.189 | 0.189 | 892 |
-| perforated-plate | genus, thin-wall, sweepable | 1 | 54 | 1644 | 0.144 | 2861 | 26822.2 | 26764.6 | 2.2e-3 | 26764.2 | ok, 776 vertices | 24 | 828 | 10.6 | 33611 | 9798 | 0 | pass | 7636 | extrude | 467 | 0.304 | 0.033 | 555 |
-| pinhole-block | small-feature, genus | 1 | 38 | 4044 | 0.073 | 4821 | 49919.4 | 49909.6 | 2.0e-4 | 49909.5 | ok, 1992 vertices | 16 | 1548 | 29.3 | 96411 | 30221 | 0 | pass | 23460 | extrude | 5236 | 0.408 | 0.392 | 1853 |
-| rounded-box | fillet | 1 | 26 | 2540 | 0.054 | 1822 | 23325.7 | 23340.5 | -6.3e-4 | 23340.8 | ok, 1272 vertices | 0 | 216 | 15.1 | 11687 | 1921 | 0 | pass | 1081 | no: Box: not supported yet: no structured hex method fits this body (extrude: cylindrical face 6 is not parallel to the direction; sweep: from face 0: 12 faces are not on the sweep path; from face 5: 12 faces are not on the sweep path; revolve: plan |  |  |  | 37 |
-| shallow-groove | near-tangent, boolean | 1 | 8 | 62 | 0.075 | 266 | 47795.2 | 47782.3 | 2.7e-4 | 47781.5 | ok, 33 vertices | 0 | 24 | 0.4 | 267 | 87 | 0 | pass | 27 | extrude | 1960 | 0.991 | 0.986 | 500 |
-| sharp-wedge | sharp, thin-wall, sweepable | 1 | 5 | 8 | 0.083 | 143 | 3200.0 | 3200.0 | 0.0e0 | 3200.0 | ok, 6 vertices | 0 | 9 | 0.0 | 11 | 3 | 0 | pass | 1 | extrude | 294 | 0.050 | 0.050 | 67 |
-| sphere-on-plate | near-tangent, boolean | 1 | 12 | 2826 | 0.098 | 3307 | 69353.3 | 69509.0 | -2.2e-3 | 69509.7 | ok, 1415 vertices | 0 | 124 | 15.4 | 4898 | 309 | 9 | pass | 736 | no: Plate: not supported yet: no structured hex method fits this body (extrude: face 6 is curved across the direction; sweep: from face 0: 6 faces are not on the sweep path; from face 5: linking face 11 has 2 boundary loops (periodic or holed linking  |  |  |  | 55 |
-| spline-vase | organic | 1 | 11 | 3890 | 0.106 | 1175 | 94891.2 | 94948.5 | -6.0e-4 | – | ok, 1947 vertices | 0 | 243 | 20.5 | 12461 | 1512 | 0 | pass | 2093 | sweep | 4554 | 0.649 | 0.649 | 846 |
-| stepped-block | blocky, sharp | 1 | 14 | 96 | 0.078 | 32 | 47000.0 | 47000.0 | 0.0e0 | 47000.0 | ok, 48 vertices | 1 | 60 | 0.5 | 223 | 41 | 0 | pass | 8 | extrude | 1416 | 0.999 | 0.999 | 158 |
-| thin-ribs | thin-wall, blocky | 1 | 26 | 92 | 0.119 | 68 | 32000.0 | 32000.0 | -4.5e-16 | 32000.0 | ok, 48 vertices | 0 | 72 | 0.6 | 171 | 33 | 0 | pass | 8 | extrude | 777 | 0.844 | 0.844 | 86 |
-| thin-wall-enclosure | thin-wall, blocky | 1 | 11 | 36 | 0.108 | 9 | 15564.0 | 15564.0 | 0.0e0 | 15564.0 | ok, 20 vertices | 0 | 28 | 0.1 | 191 | 47 | 2 | pass | 6 | extrude | 894 | 0.297 | 0.297 | 100 |
-| analytic-sphere | analytic | 1 | 6 | 2784 | 0.069 | 284 | 33353.4 | 33507.9 | -4.6e-3 | 33510.3 | ok, 1394 vertices | 0 | 96 | 7.0 | 4211 | 58 | 0 | pass | 358 | no: Sphere: not supported yet: no structured hex method fits this body (extrude: no direction to try; sweep: no planar face to sweep from; revolve: no axis to try; mapped: meshing failed: meshing failed: mesh validity gate: 8 inverted or degenerate ele |  |  |  | 35 |
-| analytic-torus | analytic, genus | 1 | 6 | 7230 | 0.109 | 627 | 37644.1 | 37893.7 | -6.6e-3 | 37899.3 | ok, 3615 vertices | 1 | 180 | 30.0 | 19790 | 466 | 0 | pass | 1091 | revolve | 1134 | 0.706 | 0.706 | 210 |
-| analytic-cylinder | analytic, sweepable | 1 | 4 | 1004 | 0.049 | 32 | 12533.1 | 12564.7 | -2.5e-3 | 12566.4 | ok, 504 vertices | 0 | 106 | 3.4 | 7437 | 1568 | 0 | pass | 344 | extrude | 3200 | 0.707 | 0.707 | 409 |
-| analytic-pipe | analytic, genus, sweepable | 1 | 6 | 1448 | 0.104 | 51 | 94150.9 | 94238.0 | -9.3e-4 | 94247.8 | ok, 724 vertices | 1 | 220 | 4.5 | 6238 | 1340 | 0 | pass | 327 | extrude | 1288 | 0.998 | 0.998 | 176 |
-| analytic-plate-hole | analytic, genus, sweepable | 1 | 8 | 112 | 0.224 | 13 | 98449.1 | 98429.5 | 2.0e-4 | 98429.2 | ok, 56 vertices | 1 | 62 | 0.5 | 1868 | 543 | 0 | pass | 69 | extrude | 321 | 0.320 | 0.142 | 49 |
+| part | categories | bodies | faces | triangles | tol. mm | build ms | enclosed mm³ | SolidCraft mm³ | rel. diff | exact mm³ | boundary | genus | feature edges | check ms | CDT tets | Steiner | cavities | tet gate | CDT ms | hex | hex elements | hex8 min SJ | hex20 min SJ | hex ms | octree hexes | min SJ | mean SJ | hex vol. diff | octree ms |
+|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
+| bent-tube | sweepable, thin-wall, genus | 1 | 18 | 4296 | 0.121 | 433 | 4213.1 | 4159.4 | 1.3e-2 | 4159.8 | ok, 2148 vertices | 1 | 268 | 11.4 | 15550 | 2338 | 0 | pass | 609 | sweep | 192 | 0.859 | 0.997 | 63 | no: Tube: meshing failed: the core's boundary has components of genus [0, 0, 0, 0], the part [1]: the part has walls 0.984 thick, thinner than 3 cells of the deepest level (0.536); use a smaller size, more levels, or a structured (swept) hex mesh |  |  |  | 572 |
+| bolted-assembly | assembly, blocky | 3 | 22 | 284 | 0.117 | 54 | 33004.7 | 32999.0 | 1.7e-4 | 32998.9 | ok, 144 vertices | 0, 1, 1 | 104 | 0.7 | 690 | 181 | 0 | pass | 35 | extrude | 357 | 0.398 | 0.222 | 42 | 41344 | 0.323 | 0.961 | 2.0e-4 | 11261 |
+| cantilever-i-beam | sweepable, thin-wall | 1 | 14 | 44 | 0.208 | 3 | 61600.0 | 61600.0 | 0.0e0 | 61600.0 | ok, 24 vertices | 0 | 36 | 0.2 | 90 | 24 | 0 | pass | 3 | extrude | 336 | 0.707 | 0.707 | 34 | no: Beam: meshing failed: the core's boundary has components of genus [0, 0], the part [0]: the part has walls 2.000 thick, thinner than 3 cells of the deepest level (0.924); use a smaller size, more levels, or a structured (swept) hex mesh |  |  |  | 908 |
+| cross-cylinders | boolean, sharp | 1 | 10 | 1395 | 0.087 | 7820 | 32197.4 | 32357.8 | -5.0e-3 | 32365.8 | error: Cross: invalid boundary surface: 28 defects; triangles 142 and 184 are oriented inconsistently across edge 0–396 at (-15.488656021315903, 4.577998041166734, 7.853092194999726); triangles 184 and 239 are oriented inconsistently across edge 0–554 at (-18.75000000000002, 5.634596901641276, 7.632255089815496); edge 103–108 at (-0.16007784063084018, 0.03984236936862052, 9.998639270596865) is shared by 3 triangles; … |  |  |  | error: Cross: invalid boundary surface: 28 defects; triangles 142 and 184 are oriented inconsistently across edge 0–396 at (-15.488656021315903, 4.577998041166734, 7.853092194999726); triangles 184 and 239 are oriented inconsistently across edge 0–554 at (-18.75000000000002, 5.634596901641276, 7.632255089815496); edge 103–108 at (-0.16007784063084018, 0.03984236936862052, 9.998639270596865) is shared by 3 triangles; … |  |  |  |  | error: Cross: invalid boundary surface: 28 defects; triangles 142 and 184 are oriented inconsistently across edge 0–396 at (-15.488656021315903, 4.577998041166734, 7.853092194999726); triangles 184 and 239 are oriented inconsistently across edge 0–554 at (-18.75000000000002, 5.634596901641276, 7.632255089815496); edge 103–108 at (-0.16007784063084018, 0.03984236936862052, 9.998639270596865) is shared by 3 triangles; … |  |  |  |  | error: Cross: invalid boundary surface: 28 defects; triangles 142 and 184 are oriented inconsistently across edge 0–396 at (-15.488656021315903, 4.577998041166734, 7.853092194999726); triangles 184 and 239 are oriented inconsistently across edge 0–554 at (-18.75000000000002, 5.634596901641276, 7.632255089815496); edge 103–108 at (-0.16007784063084018, 0.03984236936862052, 9.998639270596865) is shared by 3 triangles; … |  |  |  |  |
+| filleted-bracket | fillet, genus, sweepable | 1 | 13 | 540 | 0.088 | 388 | 32162.3 | 32145.0 | 5.4e-4 | 32144.8 | ok, 268 vertices | 2 | 135 | 1.4 | 5582 | 1613 | 0 | pass | 201 | no: Bracket: not supported yet: no structured hex method fits this body (extrude: cylindrical face 8 is not parallel to the direction; sweep: from face 1: linking face 6 has 2 boundary loops (periodic or holed linking faces are not swept); from face 6:  |  |  |  | 5 | 30676 | 0.305 | 0.963 | 8.3e-4 | 10213 |
+| filleted-bracket-blocks | fillet, genus, partitioned | 30 | 190 | 712 | 0.088 | 1006 | 32169.8 | 32145.1 | 7.7e-4 | 32144.8 | ok, 416 vertices | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 | 474 | 1.6 | 1423 | 349 | 0 | pass | 56 | blocks (30 cells, 45 shared faces) | 4986 | 0.705 | 0.705 | 542 | 149898 | 0.156 | 0.935 | 9.7e-4 | 36760 |
+| i-beam | sweepable, thin-wall | 1 | 14 | 44 | 0.208 | 3 | 61600.0 | 61600.0 | 0.0e0 | 61600.0 | ok, 24 vertices | 0 | 36 | 0.1 | 94 | 24 | 0 | pass | 4 | extrude | 336 | 0.707 | 0.707 | 33 | no: Beam: meshing failed: the core's boundary has components of genus [0, 0], the part [0]: the part has walls 2.000 thick, thinner than 3 cells of the deepest level (0.924); use a smaller size, more levels, or a structured (swept) hex mesh |  |  |  | 874 |
+| loft-transition | organic | 1 | 6 | 292 | 0.075 | 56 | 56390.2 | 56450.6 | -1.1e-3 | – | ok, 148 vertices | 0 | 40 | 0.7 | 1329 | 256 | 1 | pass | 36 | sweep | 2873 | 0.189 | 0.189 | 329 | 67188 | 0.210 | 0.959 | -1.2e-3 | 16559 |
+| perforated-plate | genus, thin-wall, sweepable | 1 | 54 | 1644 | 0.144 | 252 | 26822.2 | 26764.6 | 2.2e-3 | 26764.2 | ok, 776 vertices | 24 | 828 | 4.7 | 33611 | 9798 | 0 | pass | 1276 | extrude | 467 | 0.304 | 0.033 | 98 | 185128 | 0.300 | 0.958 | 2.5e-3 | 70289 |
+| pinhole-block | small-feature, genus | 1 | 38 | 4044 | 0.073 | 195 | 49919.4 | 49909.6 | 2.0e-4 | 49909.5 | ok, 1992 vertices | 16 | 1548 | 12.7 | 96411 | 30221 | 0 | pass | 4479 | extrude | 5236 | 0.408 | 0.392 | 587 | no: Block: meshing failed: octree hex meshing failed after 2 builds: the feature corner at [9.76, 10.180000000000001, 10.526315789473685] has no vertex; refining for another build needs more than 250000 grid cells (the grid needs more than 250000 cell |  |  |  | 2009 |
+| rounded-box | fillet | 1 | 26 | 2540 | 0.054 | 147 | 23325.7 | 23340.5 | -6.3e-4 | 23340.8 | ok, 1272 vertices | 0 | 216 | 6.3 | 11701 | 1928 | 0 | pass | 387 | no: Box: not supported yet: no structured hex method fits this body (extrude: cylindrical face 6 is not parallel to the direction; sweep: from face 0: 12 faces are not on the sweep path; from face 5: 12 faces are not on the sweep path; revolve: plan |  |  |  | 14 | 76430 | 0.561 | 0.985 | -1.3e-3 | 26673 |
+| shallow-groove | near-tangent, boolean | 1 | 8 | 62 | 0.075 | 11 | 47795.2 | 47782.3 | 2.7e-4 | 47781.5 | ok, 33 vertices | 0 | 24 | 0.2 | 267 | 87 | 0 | pass | 12 | extrude | 1960 | 0.991 | 0.986 | 196 | 2926 | 0.301 | 0.920 | 3.4e-4 | 982 |
+| sharp-wedge | sharp, thin-wall, sweepable | 1 | 5 | 8 | 0.083 | 2 | 3200.0 | 3200.0 | 0.0e0 | 3200.0 | ok, 6 vertices | 0 | 9 | 0.0 | 11 | 3 | 0 | pass | 0 | extrude | 294 | 0.050 | 0.050 | 30 | 84952 | 0.003 | 0.980 | -3.5e-2 | 23994 |
+| sphere-on-plate | near-tangent, boolean | 1 | 12 | 2826 | 0.098 | 686 | 69353.3 | 69509.0 | -2.2e-3 | 69509.7 | ok, 1415 vertices | 0 | 124 | 5.9 | 4898 | 309 | 9 | pass | 227 | no: Plate: not supported yet: no structured hex method fits this body (extrude: face 6 is curved across the direction; sweep: from face 0: 6 faces are not on the sweep path; from face 5: linking face 11 has 2 boundary loops (periodic or holed linking  |  |  |  | 10 | 25194 | 0.065 | 0.945 | 4.4e-2 | 20477 |
+| spline-vase | organic | 1 | 11 | 3890 | 0.106 | 217 | 94891.2 | 94948.5 | -6.0e-4 | – | ok, 1947 vertices | 0 | 243 | 8.6 | 12461 | 1512 | 0 | pass | 439 | sweep | 4554 | 0.649 | 0.649 | 522 | no: Vase: limit reached: the grid needs more than 500000 cells |  |  |  | 2308 |
+| stepped-block | blocky, sharp | 1 | 14 | 96 | 0.078 | 29 | 47000.0 | 47000.0 | 0.0e0 | 47000.0 | ok, 48 vertices | 1 | 60 | 0.4 | 223 | 41 | 0 | pass | 7 | extrude | 1416 | 0.999 | 0.999 | 137 | 52274 | 0.076 | 0.975 | -2.7e-6 | 14558 |
+| thin-ribs | thin-wall, blocky | 1 | 26 | 92 | 0.119 | 62 | 32000.0 | 32000.0 | -4.5e-16 | 32000.0 | ok, 48 vertices | 0 | 72 | 0.5 | 171 | 33 | 0 | pass | 7 | extrude | 777 | 0.844 | 0.844 | 76 | no: Plate: limit reached: the grid needs more than 500000 cells |  |  |  | 1111 |
+| thin-wall-enclosure | thin-wall, blocky | 1 | 11 | 36 | 0.108 | 8 | 15564.0 | 15564.0 | 0.0e0 | 15564.0 | ok, 20 vertices | 0 | 28 | 0.1 | 191 | 47 | 2 | pass | 5 | extrude | 894 | 0.297 | 0.297 | 89 | no: Enclosure: meshing failed: the core's boundary has components of genus [0, 0], the part [0]: the part has walls 1.000 thick, thinner than 3 cells of the deepest level (0.479); use a smaller size, more levels, or a structured (swept) hex mesh |  |  |  | 1623 |
+| analytic-sphere | analytic | 1 | 6 | 2784 | 0.069 | 247 | 33353.4 | 33507.9 | -4.6e-3 | 33510.3 | ok, 1394 vertices | 0 | 96 | 5.9 | 4211 | 58 | 0 | pass | 266 | no: Sphere: not supported yet: no structured hex method fits this body (extrude: no direction to try; sweep: no planar face to sweep from; revolve: no axis to try; mapped: meshing failed: meshing failed: mesh validity gate: 8 inverted or degenerate ele |  |  |  | 23 | 2400 | 0.538 | 0.866 | -1.3e-2 | 966 |
+| analytic-torus | analytic, genus | 1 | 6 | 7230 | 0.109 | 337 | 37644.1 | 37893.7 | -6.6e-3 | 37899.3 | ok, 3615 vertices | 1 | 180 | 15.8 | 19790 | 466 | 0 | pass | 761 | revolve | 1134 | 0.706 | 0.706 | 168 | 18814 | 0.460 | 0.881 | -1.2e-2 | 5646 |
+| analytic-cylinder | analytic, sweepable | 1 | 4 | 1004 | 0.049 | 26 | 12533.1 | 12564.7 | -2.5e-3 | 12566.4 | ok, 504 vertices | 0 | 106 | 2.7 | 7437 | 1568 | 0 | pass | 267 | extrude | 3200 | 0.707 | 0.707 | 331 | 2424 | 0.532 | 0.919 | -9.1e-3 | 699 |
+| analytic-pipe | analytic, genus, sweepable | 1 | 6 | 1448 | 0.104 | 44 | 94150.9 | 94238.0 | -9.3e-4 | 94247.8 | ok, 724 vertices | 1 | 220 | 3.8 | 6238 | 1340 | 0 | pass | 272 | extrude | 1288 | 0.998 | 0.998 | 145 | 53680 | 0.520 | 0.944 | -7.8e-4 | 14728 |
+| analytic-plate-hole | analytic, genus, sweepable | 1 | 8 | 112 | 0.224 | 11 | 98449.1 | 98429.5 | 2.0e-4 | 98429.2 | ok, 56 vertices | 1 | 62 | 0.4 | 1868 | 543 | 0 | pass | 53 | extrude | 321 | 0.320 | 0.142 | 38 | 163664 | 0.309 | 0.989 | 2.2e-4 | 51714 |
 
 ## Categories
 

docs/meshing/literature.md

@@ -531,6 +531,53 @@ Delaunay methods with Voronoi/ODT improvement, then (2) all-hex and hex-dominant
   with conforming templates and no dualization.
 - **Use.** Lets us refine a coarse aligned hex mesh, for example polycube output, as well as an octree.
 - **Code.** `cg3hci/Gen-Adapt-Ref-for-Hexmeshing` is **MIT**: READ OK.
+- **Finding (2026-10).** The reference pipeline chooses transition vertices with an integer linear
+  program (Gurobi) and installs CinoLib's schemes by heuristics that are retried with other
+  parameters when they fail; neither was adopted (see B1.10).
+
+**B1.7 Schneiders, Algorithms for quadrilateral and hexahedral mesh generation (VKI lecture notes,
+2000), §5.** **Verdict: built on (the transition templates of `fieldcraft-mesh::octree`).**
+- **Link.** https://www.robertschneiders.de/papers/vki.pdf (author's copy).
+- **Method.** 3-refinement driven by nodal marks: an element is split by the template of its
+  marked corners, and faces split by their own marks only, so the refinement is conforming.
+- **Limit (stated there).** The 2D templates of the original set make some 3D mark patterns
+  impossible: a hexahedron with seven marked corners has 51 boundary quads, an odd number, so no
+  hexahedral mesh of it exists (Mitchell's parity condition). Concave refined regions therefore
+  need their concave cells refined. The fix is a modified edge rule (an edge with a marked end
+  splits in three) with buffer layers in the face templates.
+
+**B1.8 Ito, Shih, Soni 2009, Octree-based reasonable-quality hexahedral mesh generation using a
+new set of refinement templates.** **Verdict: reference.**
+- **Citation.** IJNME 77(13):1809–1833, 2009. Templates for concave edges of the refined region
+  (paper not obtained; cited after Tong and Zhang 2025).
+
+**B1.9 Tong, Halilaj, Zhang 2024, HybridOctree_Hex.** **Verdict: reference (buffer-zone rule adopted).**
+- **Citation.** H. Tong, E. Halilaj, Y. J. Zhang, *HybridOctree_Hex: Hybrid octree-based adaptive
+  all-hexahedral mesh generation with Jacobian control*, arXiv 2401.05984, 2024.
+- **Method.** Dual of a strongly balanced, paired octree with five templates (a 13-hexahedron
+  frustum per transition face plus four transition-edge templates); a core cleared at the
+  boundary so that the normals of the quads around each core boundary vertex lie in one
+  half-space, then a buffer layer projected onto the surface and optimised.
+- **Code.** `CMU-CBML/HybridOctree_Hex` is **BSD-3-Clause**: READ OK (hard-coded template
+  geometry, not adopted).
+
+**B1.10 Tong, Zhang 2025, MCHex.** **Verdict: reference (next boundary stage to evaluate).**
+- **Citation.** H. Tong, Y. J. Zhang, *MCHex: Marching cubes based adaptive hexahedral mesh
+  generation with guaranteed positive Jacobian*, arXiv 2511.02064, 2025.
+- **Method.** The same 3-refinement templates (B1.7); the boundary cells are cut by marching
+  cubes on a signed distance field, every cut cell's pieces are 3-regular polyhedra, and their
+  midpoint subdivision gives hexahedra with a proven positive Jacobian (Bernstein coefficients) and
+  a boundary that is pillowed by construction; 17 of the 31 cube cases are refined locally first.
+- **Limits (stated there).** Sharp edges and corners are lost (marching cubes smooths them), and
+  the convex-only transitions over-refine. ModelCraft keeps the core-padding-pillow boundary with
+  feature capture for CAD parts; MCHex's cut cells are the candidate for organic shapes.
+
+**Implemented (`fieldcraft-mesh::octree`, `mesh.hex.octree`, 2026-10).** B1.7 templates on a
+27-tree refined by curvature, thickness and narrow faces; a core cleared by the B1.9 rule; the
+feature map of B1.4 (patches first); padding plus a pillow per patch (B1.4 §3.4); the B6.4 energy
+with a fitting term (B1.4 §3.5); the certified check of B6.6 with local repair. The dual schemes
+of B1.5 were not adopted (B1.6 finding); concave refined regions are closed by refinement (B1.7
+limit), so curvature refines one level only.
 
 ---
 

docs/oracle.md

@@ -2,7 +2,7 @@
 
 Generated by `cargo xtask oracle` from [`tests/oracle/`](../tests/oracle/); the policy is [verification.md](verification.md). Every case is our own model, built with FieldCraft commands from the problem's published description, and compared with the reference value its `expected.json` cites, within the stated tolerance and for the stated reason. Accepted results are frozen in each case's `fixture.json`: a later change that moves one beyond its band (10⁻⁹ of its size) fails CI even while the reference still passes. No reference value was produced with a commercial product.
 
-**163 of 169 quantities pass** (6 known failure(s)) in 49 active cases: 40 theory, 9 published benchmark, 0 experiment; 9 pending cases wait for capabilities still being built.
+**170 of 176 quantities pass** (6 known failure(s)) in 50 active cases: 41 theory, 9 published benchmark, 0 experiment; 9 pending cases wait for capabilities still being built.
 
 ## Theory: closed-form solutions, patch tests and mesh invariants
 
@@ -90,6 +90,13 @@ Generated by `cargo xtask oracle` from [`tests/oracle/`](../tests/oracle/); the
 | [nl-thermal-expansion](../tests/oracle/nl-thermal-expansion/) | stress of the held block, σ_xx (hydrostatic) | -3.5108 | -3.5108 | < 1e-7% | 1e-8 rel. | pass | matches |
 | [nl-yeoh-uniaxial](../tests/oracle/nl-yeoh-uniaxial/) | axial Cauchy stress, stretch λ = 2 | 3.1743 | 3.1743 | < 1e-7% | 1e-6 rel. | pass | matches |
 | [nl-yeoh-uniaxial](../tests/oracle/nl-yeoh-uniaxial/) | lateral displacement at (1, 1, 1): l − 1 | -0.28545 | -0.28545 | < 1e-7% | 1e-6 rel. | pass | matches |
+| [octree-filleted-bracket](../tests/oracle/octree-filleted-bracket/) | x reaction on the clamped face (5 mm mesh) | -2000 | -2000 | < 1e-7% | 1e-6 rel. | pass | matches |
+| [octree-filleted-bracket](../tests/oracle/octree-filleted-bracket/) | elements the certified Jacobian check does not prove valid (5 mm mesh) | 0 | 0 | ≈ 0 | exact | pass | matches |
+| [octree-filleted-bracket](../tests/oracle/octree-filleted-bracket/) | mesh volume against the CAD volume (5 mm mesh) | 32145 | 32177 | 0.1% | 0.3% | pass | matches |
+| [octree-filleted-bracket](../tests/oracle/octree-filleted-bracket/) | σ_yy on the outer face of the upright leg at y = 19, z = 6 (5 mm mesh) | 145.31 | 149.26 | 2.72% | 5% | pass | matches |
+| [octree-filleted-bracket](../tests/oracle/octree-filleted-bracket/) | σ_yy on the inner face of the upright leg at y = 19, z = 6 (5 mm mesh) | -145.31 | -147.28 | -1.35% | 5% | pass | matches |
+| [octree-filleted-bracket](../tests/oracle/octree-filleted-bracket/) | von Mises stress at the middle of the fillet, 7 mm mesh | 143.21 | 138.53 | -3.27% | 5% | pass | matches |
+| [octree-filleted-bracket](../tests/oracle/octree-filleted-bracket/) | von Mises stress at the middle of the fillet, 5 mm mesh | 143.21 | 144.53 | 0.921% | 2% | pass | matches |
 | [patch-hex20](../tests/oracle/patch-hex20/) | u_x at the corner (1, 1, 1) | -2.5000e-4 | -2.5000e-4 | < 1e-7% | 1e-9 rel. | pass | matches |
 | [patch-hex20](../tests/oracle/patch-hex20/) | u_y at the corner (1, 1, 1) | 0.001 | 0.001 | < 1e-7% | 1e-9 rel. | pass | matches |
 | [patch-hex20](../tests/oracle/patch-hex20/) | u_z at the corner (1, 1, 1) | 0.00225 | 0.00225 | < 1e-7% | 1e-9 rel. | pass | matches |
@@ -220,6 +227,7 @@ Generated by `cargo xtask oracle` from [`tests/oracle/`](../tests/oracle/); the
 - **[nl-rigid-rotation](../tests/oracle/nl-rigid-rotation/)** (theory): A quarter turn strains and stresses nothing (objectivity). Source: Objectivity of the total Lagrangian formulation and of the hyperelastic stress: a rigid motion x = R X + c gives C = I and σ = 0 whatever the rotation (Bonet & Wood 2008, §3.4 and §5.2; T. Belytschko, W. K. Liu, B. Moran, Nonlinear Finite Elements for Continua and Structures, 2nd ed., Wiley 2014, §5.5). Tolerances — ux: the inner nodes follow the rigid rotation exactly; uy: the inner nodes follow the rigid rotation exactly; uz: the inner nodes follow the rigid rotation exactly; von_mises: zero up to rounding of stresses of order E = 1000: a 90° turn leaves no strain
 - **[nl-thermal-expansion](../tests/oracle/nl-thermal-expansion/)** (theory): Free and fully constrained thermal expansion at finite strain. Source: The multiplicative split F = F_m F_θ with F_θ = (1 + α ΔT) I (V. A. Lubarda, Appl. Mech. Rev. 57(2):95–108, 2004): a free body follows F = F_θ stress-free; held on every face, F = I and F_m = ϑ⁻¹ I, so the compressible neo-Hookean material carries the pressure of its volume ratio, σ = κ (ϑ⁻³ − 1). Tolerances — free_u: a homogeneous stretch every element represents exactly; free_stress: rounding only; held_stress: a homogeneous state; E = 10, ν = 0.3 (κ = 8.3333), ϑ = 1.2
 - **[nl-yeoh-uniaxial](../tests/oracle/nl-yeoh-uniaxial/)** (theory): Yeoh block in uniaxial tension. Source: Uniaxial tension of the decoupled compressible Yeoh energy W = C₁₀ (Ī₁ − 3) + C₂₀ (Ī₁ − 3)² + C₃₀ (Ī₁ − 3)³ + κ/2 (J − 1)² (O. H. Yeoh, Rubber Chem. Technol. 66(5):754–771, 1993): σ = (λ/J) ∂W/∂λ with the lateral stretch found so σ_yy = 0 (bisection; E = 3, ν = 0.49, C₁₀ = μ/2, C₂₀ = −0.02, C₃₀ = 0.004). Tolerances — sigma: a homogeneous deformation the hybrid hex20 represents exactly: the tolerance admits the Newton tolerance and rounding; lateral: a homogeneous deformation the hybrid hex20 represents exactly: the tolerance admits the Newton tolerance and rounding
+- **[octree-filleted-bracket](../tests/oracle/octree-filleted-bracket/)** (theory): Filleted L-bracket on the octree all-hex mesh: equilibrium, beam bending, fillet stress convergence. Source: equilibrium; Euler–Bernoulli bending of the 8 × 40 mm leg, σ = M c / I with M = F (50 − y) (S. P. Timoshenko, J. N. Goodier, Theory of Elasticity, 3rd ed., 1970, §102); grid convergence of the fillet stress against the 4 mm mesh of this case (P. J. Roache, Verification and Validation in Computational Science and Engineering, 1998, ch. 5). Tolerances — reaction_x: the reactions balance the 2 kN load to solver precision; invalid: the octree mesher returns only elements proven valid; volume: bilinear boundary faces on the R6 fillet and the two Ø8 holes, and the 1e-3 chordal tessellation they are fitted to; syy_outer: beam theory ignores the corner 5 mm below and the hole 9 mm above; the probe takes the nearest node's extrapolated stress; syy_inner: as for the outer face; fillet_7: the coarsest mesh, within 5% of the 4 mm mesh; fillet_5: converging: within 2% of the 4 mm mesh
 - **[patch-hex20](../tests/oracle/patch-hex20/)** (theory): Constant-stress patch test, hex20. Source: Irons and Razzaque, 'Experience with the patch test for convergence of finite elements', in Aziz (ed.), The Mathematical Foundations of the Finite Element Method, Academic Press 1972, pp. 557–587; MacNeal and Harder, Finite Elements in Analysis and Design 1 (1985) 3–20, constant-stress loading, p. 14 Table 7. Tolerances — ux: a patch test is passed exactly: the tolerance only admits rounding; uy: a patch test is passed exactly: the tolerance only admits rounding; uz: a patch test is passed exactly: the tolerance only admits rounding; sxx: a patch test is passed exactly: the tolerance only admits rounding; syy: a patch test is passed exactly: the tolerance only admits rounding; szz: a patch test is passed exactly: the tolerance only admits rounding; sxy: a patch test is passed exactly: the tolerance only admits rounding (stresses of order 1)
 - **[patch-hex20-reduced](../tests/oracle/patch-hex20-reduced/)** (theory): Constant-stress patch test, hex20-reduced. Source: Irons and Razzaque, 'Experience with the patch test for convergence of finite elements', in Aziz (ed.), The Mathematical Foundations of the Finite Element Method, Academic Press 1972, pp. 557–587; MacNeal and Harder, Finite Elements in Analysis and Design 1 (1985) 3–20, constant-stress loading, p. 14 Table 7. Tolerances — ux: a patch test is passed exactly: the tolerance only admits rounding; uy: a patch test is passed exactly: the tolerance only admits rounding; uz: a patch test is passed exactly: the tolerance only admits rounding; sxx: a patch test is passed exactly: the tolerance only admits rounding; syy: a patch test is passed exactly: the tolerance only admits rounding; szz: a patch test is passed exactly: the tolerance only admits rounding; sxy: a patch test is passed exactly: the tolerance only admits rounding (stresses of order 1)
 - **[patch-hex8](../tests/oracle/patch-hex8/)** (theory): Constant-stress patch test, hex8. Source: Irons and Razzaque, 'Experience with the patch test for convergence of finite elements', in Aziz (ed.), The Mathematical Foundations of the Finite Element Method, Academic Press 1972, pp. 557–587; MacNeal and Harder, Finite Elements in Analysis and Design 1 (1985) 3–20, constant-stress loading, p. 14 Table 7. Tolerances — ux: a patch test is passed exactly: the tolerance only admits rounding; uy: a patch test is passed exactly: the tolerance only admits rounding; uz: a patch test is passed exactly: the tolerance only admits rounding; sxx: a patch test is passed exactly: the tolerance only admits rounding; syy: a patch test is passed exactly: the tolerance only admits rounding; szz: a patch test is passed exactly: the tolerance only admits rounding; sxy: a patch test is passed exactly: the tolerance only admits rounding (stresses of order 1)

tests/oracle/octree-filleted-bracket/README.md

@@ -0,0 +1,21 @@
+# Filleted L-bracket on the octree all-hex mesh
+
+The corpus part `filleted-bracket` (an 8 mm thick L, 60 × 50 × 40 mm, an R6 fillet inside the
+corner and a Ø8 hole through each leg), built with SolidCraft from `design.json` and meshed by
+`mesh.hex.octree` (hex8 with incompatible modes, steel E = 210 GPa, ν = 0.3). The bottom face of
+the horizontal leg (`F2:side:l1`) is clamped and the top face of the upright leg (`F2:side:l5`)
+carries 2 kN along +x.
+
+Checked on the 5 mm mesh: equilibrium; the certified validity of every element; the mesh volume
+against the CAD volume; the bending stress of the upright leg at y = 19 on both faces,
+σ_yy = M c / I = 2000 · 31 · 4 / 1706.7 = 145.3 MPa (M is exact by statics, the plane-section
+stress is the reference, the corner and a hole are near). The von Mises stress at the middle of
+the fillet, which has no closed form, is checked for convergence against the 4 mm mesh:
+
+| mesh size | hexahedra | fillet von Mises (MPa) | σ_yy outer / inner at y = 19 (MPa) |
+|---|---|---|---|
+| 7 mm | 4291 | 138.5 | 149.4 / −146.9 |
+| 5 mm | 13273 | 144.5 | 149.3 / −147.3 |
+| 4 mm (reference) | 21508 | 143.2 | 147.4 / −145.5 |
+
+A 3 mm mesh needs more than the direct solver's factor cap.

tests/oracle/octree-filleted-bracket/design.json

@@ -0,0 +1,89 @@
+{
+ "commands": [
+  {
+   "command": "sketch.create",
+   "params": {
+    "plane": "XY",
+    "name": "Profile"
+   }
+  },
+  {
+   "command": "sketch.line",
+   "params": {
+    "points": [
+     [
+      0,
+      0
+     ],
+     [
+      60,
+      0
+     ],
+     [
+      60,
+      8
+     ],
+     [
+      8,
+      8
+     ],
+     [
+      8,
+      50
+     ],
+     [
+      0,
+      50
+     ]
+    ],
+    "closed": true
+   }
+  },
+  {
+   "command": "sketch.finish",
+   "params": {}
+  },
+  {
+   "command": "solid.extrude",
+   "params": {
+    "distance": 40,
+    "body_name": "Bracket"
+   }
+  },
+  {
+   "command": "solid.fillet",
+   "params": {
+    "edges": [
+     [
+      8,
+      8,
+      20
+     ]
+    ],
+    "radius": 6
+   }
+  },
+  {
+   "command": "solid.hole",
+   "params": {
+    "position": [
+     40,
+     8,
+     20
+    ],
+    "diameter": 8
+   }
+  },
+  {
+   "command": "solid.hole",
+   "params": {
+    "position": [
+     8,
+     32,
+     20
+    ],
+    "diameter": 8
+   }
+  }
+ ]
+}

tests/oracle/octree-filleted-bracket/expected.json

@@ -0,0 +1,99 @@
+{
+ "title": "Filleted L-bracket on the octree all-hex mesh: equilibrium, beam bending, fillet stress convergence",
+ "kind": "theory",
+ "source": "equilibrium; Euler–Bernoulli bending of the 8 × 40 mm leg, σ = M c / I with M = F (50 − y) (S. P. Timoshenko, J. N. Goodier, Theory of Elasticity, 3rd ed., 1970, §102); grid convergence of the fillet stress against the 4 mm mesh of this case (P. J. Roache, Verification and Validation in Computational Science and Engineering, 1998, ch. 5)",
+ "slow": true,
+ "quantities": [
+  {
+   "id": "reaction_x",
+   "quantity": "x reaction on the clamped face (5 mm mesh)",
+   "measure": {
+    "from": "solve5",
+    "pointer": "/reaction_total/0"
+   },
+   "value": -2000.0,
+   "tolerance": {
+    "relative": 1e-06
+   },
+   "reason": "the reactions balance the 2 kN load to solver precision"
+  },
+  {
+   "id": "invalid",
+   "quantity": "elements the certified Jacobian check does not prove valid (5 mm mesh)",
+   "measure": {
+    "from": "quality5",
+    "pointer": "/invalid/count"
+   },
+   "value": 0.0,
+   "tolerance": {
+    "exact": true
+   },
+   "reason": "the octree mesher returns only elements proven valid"
+  },
+  {
+   "id": "volume",
+   "quantity": "mesh volume against the CAD volume (5 mm mesh)",
+   "measure": {
+    "from": "mesh5",
+    "pointer": "/bodies/0/volume"
+   },
+   "value": 32144.778925,
+   "tolerance": {
+    "relative": 0.003
+   },
+   "reason": "bilinear boundary faces on the R6 fillet and the two Ø8 holes, and the 1e-3 chordal tessellation they are fitted to"
+  },
+  {
+   "id": "syy_outer",
+   "quantity": "σ_yy on the outer face of the upright leg at y = 19, z = 6 (5 mm mesh)",
+   "measure": {
+    "from": "probe5",
+    "pointer": "/nodes/1/stress/1"
+   },
+   "value": 145.3125,
+   "tolerance": {
+    "relative": 0.05
+   },
+   "reason": "beam theory ignores the corner 5 mm below and the hole 9 mm above; the probe takes the nearest node's extrapolated stress"
+  },
+  {
+   "id": "syy_inner",
+   "quantity": "σ_yy on the inner face of the upright leg at y = 19, z = 6 (5 mm mesh)",
+   "measure": {
+    "from": "probe5",
+    "pointer": "/nodes/2/stress/1"
+   },
+   "value": -145.3125,
+   "tolerance": {
+    "relative": 0.05
+   },
+   "reason": "as for the outer face"
+  },
+  {
+   "id": "fillet_7",
+   "quantity": "von Mises stress at the middle of the fillet, 7 mm mesh",
+   "measure": {
+    "from": "probe7",
+    "pointer": "/nodes/0/von_mises"
+   },
+   "value": 143.21,
+   "tolerance": {
+    "relative": 0.05
+   },
+   "reason": "the coarsest mesh, within 5% of the 4 mm mesh"
+  },
+  {
+   "id": "fillet_5",
+   "quantity": "von Mises stress at the middle of the fillet, 5 mm mesh",
+   "measure": {
+    "from": "probe5",
+    "pointer": "/nodes/0/von_mises"
+   },
+   "value": 143.21,
+   "tolerance": {
+    "relative": 0.02
+   },
+   "reason": "converging: within 2% of the 4 mm mesh"
+  }
+ ]
+}

tests/oracle/octree-filleted-bracket/fixture.json

@@ -0,0 +1,33 @@
+{
+  "note": "FieldCraft's accepted results for this case, written by `cargo xtask oracle --bless`. A result outside its band fails CI.",
+  "values": {
+    "fillet_5": {
+      "band": 1.4452906175564113e-7,
+      "value": 144.52906175564112
+    },
+    "fillet_7": {
+      "band": 1.4321000000000002e-7,
+      "value": 138.5286423070889
+    },
+    "invalid": {
+      "band": 0.0,
+      "value": 0.0
+    },
+    "reaction_x": {
+      "band": 2.0000000000000003e-6,
+      "value": -1999.999999998785
+    },
+    "syy_inner": {
+      "band": 1.472793063397077e-7,
+      "value": -147.2793063397077
+    },
+    "syy_outer": {
+      "band": 1.4926027331710374e-7,
+      "value": 149.26027331710372
+    },
+    "volume": {
+      "band": 0.00003217691091849411,
+      "value": 32176.910918494108
+    }
+  }
+}

tests/oracle/octree-filleted-bracket/recipe.json

@@ -0,0 +1,228 @@
+[
+ {
+  "command": "model.material.create",
+  "params": {
+   "name": "Steel",
+   "elastic": {
+    "type": "isotropic",
+    "young": 210000.0,
+    "poisson": 0.3
+   }
+  }
+ },
+ {
+  "command": "mesh.hex.octree",
+  "params": {
+   "body": "Bracket",
+   "size": 7.0
+  },
+  "id": "mesh7"
+ },
+ {
+  "command": "model.set.create",
+  "params": {
+   "name": "all",
+   "geometry": [
+    {
+     "body": "Bracket"
+    }
+   ]
+  }
+ },
+ {
+  "command": "model.set.create",
+  "params": {
+   "name": "base",
+   "geometry": [
+    {
+     "body": "Bracket",
+     "face": "F2:side:l1"
+    }
+   ]
+  }
+ },
+ {
+  "command": "model.set.create",
+  "params": {
+   "name": "top",
+   "geometry": [
+    {
+     "body": "Bracket",
+     "face": "F2:side:l5"
+    }
+   ]
+  }
+ },
+ {
+  "command": "model.section.create",
+  "params": {
+   "name": "solid",
+   "set": "all",
+   "material": "Steel",
+   "kind": {
+    "type": "solid"
+   },
+   "formulation": {
+    "incompatible_modes": true
+   }
+  }
+ },
+ {
+  "command": "model.boundary_condition.create",
+  "params": {
+   "name": "clamp",
+   "type": "fixed",
+   "target": {
+    "set": "base"
+   }
+  }
+ },
+ {
+  "command": "model.load.create",
+  "params": {
+   "name": "F",
+   "type": "force",
+   "target": {
+    "set": "top"
+   },
+   "vector": [
+    2000.0,
+    0,
+    0
+   ]
+  }
+ },
+ {
+  "command": "model.step.create",
+  "params": {
+   "name": "Load",
+   "procedure": {
+    "type": "static"
+   },
+   "boundary_conditions": [
+    "clamp"
+   ],
+   "loads": [
+    "F"
+   ]
+  }
+ },
+ {
+  "command": "mesh.quality",
+  "params": {},
+  "id": "quality7"
+ },
+ {
+  "command": "solve.static",
+  "params": {},
+  "id": "solve7"
+ },
+ {
+  "command": "solve.probe",
+  "params": {
+   "step": "Load",
+   "at": [
+    [
+     9.757359312880714,
+     9.757359312880714,
+     10.0
+    ],
+    [
+     0.0,
+     19.0,
+     6.0
+    ],
+    [
+     8.0,
+     19.0,
+     6.0
+    ]
+   ]
+  },
+  "id": "probe7"
+ },
+ {
+  "command": "mesh.hex.octree",
+  "params": {
+   "body": "Bracket",
+   "size": 5.0
+  },
+  "id": "mesh5"
+ },
+ {
+  "command": "mesh.quality",
+  "params": {},
+  "id": "quality5"
+ },
+ {
+  "command": "solve.static",
+  "params": {},
+  "id": "solve5"
+ },
+ {
+  "command": "solve.probe",
+  "params": {
+   "step": "Load",
+   "at": [
+    [
+     9.757359312880714,
+     9.757359312880714,
+     10.0
+    ],
+    [
+     0.0,
+     19.0,
+     6.0
+    ],
+    [
+     8.0,
+     19.0,
+     6.0
+    ]
+   ]
+  },
+  "id": "probe5"
+ },
+ {
+  "command": "mesh.hex.octree",
+  "params": {
+   "body": "Bracket",
+   "size": 4.0
+  },
+  "id": "mesh4"
+ },
+ {
+  "command": "mesh.quality",
+  "params": {},
+  "id": "quality4"
+ },
+ {
+  "command": "solve.static",
+  "params": {},
+  "id": "solve4"
+ },
+ {
+  "command": "solve.probe",
+  "params": {
+   "step": "Load",
+   "at": [
+    [
+     9.757359312880714,
+     9.757359312880714,
+     10.0
+    ],
+    [
+     0.0,
+     19.0,
+     6.0
+    ],
+    [
+     8.0,
+     19.0,
+     6.0
+    ]
+   ]
+  },
+  "id": "probe4"
+ }
+]