Exact 3D Delaunay tetrahedralization (Bowyer–Watson, ghost hull, BRIO/Hilbert, verifier) (#14)

* mesh: exact Bowyer-Watson Delaunay tetrahedralization with ghost tets, BRIO/Hilbert order and a verifier

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

* mesh: Delaunay benchmark, hash-paired cavity faces, mesh.delaunay.points debug command

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

* mesh: simplify Delaunay helpers, share predicates and parsers, halve verifier work

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

---------

Co-authored-by: jyviko <2501648+jyviko@users.noreply.github.com>
Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>

ROADMAP.md

@@ -18,7 +18,7 @@ hex meshes today and it gives the best elements where stresses are read.
 |---|---|---|
 | M0 Foundation | Workspace, quality gates, FE model, command engine with rollback and panic guard, CLI, geometry boundary on SolidCraft | done |
 | M1 Geometry in and surface preparation | Open `.solidcraft` designs with persistent face names, assemblies placed (done); face-tagged 3MF (done); test corpus, `cargo xtask corpus` (done); exact predicates (done); quality metrics and certified validity gate (done); weld, closed/manifold/orientation checks, feature edges | in progress |
-| M2 Tet meshing and size field | Size field (global, per face and edge seeds, curvature and proximity); Bowyer–Watson Delaunay; robust CDT boundary recovery (Diazzi et al. 2023); Delaunay refinement, ODT smoothing, sliver removal; tet4 and tet10 with mid-side nodes on the B-rep; mesh verify per set | in progress |
+| M2 Tet meshing and size field | Size field (global, per face and edge seeds, curvature and proximity); Bowyer–Watson Delaunay (exact, symbolic perturbation, verified; `modelcraft-mesh::delaunay`, done); robust CDT boundary recovery (Diazzi et al. 2023); Delaunay refinement, ODT smoothing, sliver removal; tet4 and tet10 with mid-side nodes on the B-rep; mesh verify per set | in progress |
 | 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; writers for Abaqus input, Nastran SOL 101/103 and LS-DYNA implicit; the `.modelcraft` file | planned |
 | M4 Own solver and first results | Linear static and modal solver (tet4/10, hex8, hex20; sparse direct, Lanczos), verified against analytic solutions; deformed shape, contours, mode shapes, probes, reaction sums | planned |
 | M5 Structured hex | 2D quad meshing of faces; mapped hex of six-sided regions; sweep meshing; conforming interfaces and tet–hex transition; automatic decomposition into sweepable and mappable regions, proposed as SolidCraft splits | planned |

crates/engine/src/cmd/mesh.rs

@@ -1,8 +1,9 @@
+use modelcraft_mesh::delaunay::Triangulation;
 use modelcraft_mesh::quality::{self, QualityReport};
 use modelcraft_model::{ElementKind, MAX_ELEMENTS};
 use serde_json::{Value, json};
 
-use super::{label, labels, xyz};
+use super::{label, labels, points, xyz};
 use crate::{CommandSpec, EngineError, Result, Session};
 
 pub const NODE_ADD: CommandSpec = CommandSpec {
@@ -32,6 +33,30 @@ pub const QUALITY: CommandSpec = CommandSpec {
     run: quality_report,
 };
 
+pub const DELAUNAY_POINTS: CommandSpec = CommandSpec {
+    id: "mesh.delaunay.points",
+    label: "Delaunay of Points",
+    params: r#"{"points": [[x, y, z], ...] (at most 20000)} — the exact Delaunay tetrahedralization (a debugging report; the model is not changed): {"dimension", "vertices", "duplicates": [[point, kept point], ...], "tets": [[i, j, k, l], ...]} with 0-based point indices in tet4 order; collinear or coplanar points give no tets"#,
+    run: delaunay_points,
+};
+
+/// The most points `mesh.delaunay.points` takes (its report lists every tetrahedron).
+const DELAUNAY_MAX_POINTS: usize = 20_000;
+
+fn delaunay_points(_: &mut Session, p: &Value) -> Result<Value> {
+    let points = points(p, "points", DELAUNAY_MAX_POINTS)?;
+    let tri = Triangulation::new(&points)?;
+    tri.verify()?;
+    let duplicates: Vec<[u32; 2]> = (0..points.len() as u32).filter_map(|v| tri.duplicate_of(v).map(|d| [v, d])).collect();
+    let tets: Vec<[u32; 4]> = tri.finite_tets().map(|(_, vs)| vs).collect();
+    Ok(json!({
+        "dimension": tri.dimension(),
+        "vertices": points.len() - duplicates.len(),
+        "duplicates": duplicates,
+        "tets": tets,
+    }))
+}
+
 fn node_add(s: &mut Session, p: &Value) -> Result<Value> {
     let l = label(p, "label")?;
     s.model.set_node(l, xyz(p, "xyz")?)?;
@@ -127,6 +152,32 @@ mod tests {
         s.execute("mesh.element.add", &json!({"label": label, "kind": "hex8", "nodes": nodes})).unwrap();
     }
 
+    #[test]
+    fn delaunay_of_points_reports_tets_and_duplicates() {
+        let mut s = Session::default();
+        let cube = json!({"points": [[0, 0, 0], [1, 0, 0], [0, 1, 0], [1, 1, 0], [0, 0, 1], [1, 0, 1], [0, 1, 1], [1, 1, 1], [1, 1, 1]]});
+        let r = s.execute("mesh.delaunay.points", &cube).unwrap();
+        assert_eq!(r["dimension"], 3);
+        assert_eq!(r["vertices"], 8);
+        assert_eq!(r["duplicates"], json!([[8, 7]]));
+        let tets = r["tets"].as_array().unwrap();
+        assert!((5..=6).contains(&tets.len()), "{tets:?}");
+        assert!(s.journal.is_empty(), "a report does not change the model");
+
+        let flat = s.execute("mesh.delaunay.points", &json!({"points": [[0, 0, 0], [1, 0, 0], [0, 1, 0], [1, 1, 0]]})).unwrap();
+        assert_eq!((flat["dimension"].clone(), flat["tets"].clone()), (json!(2), json!([])));
+        for bad in [
+            json!({}),
+            json!({"points": [[0, 0]]}),
+            json!({"points": [[0, 0, "x"]]}),
+            json!({"points": [[f64::MAX, 0, 0]]}),
+            json!({"points": [[1e-300, 0, 0]]}),
+            json!({"points": vec![[0, 0, 0]; 20_001]}),
+        ] {
+            assert!(s.execute("mesh.delaunay.points", &bad).is_err(), "{bad}");
+        }
+    }
+
     #[test]
     fn reports_quality_and_inverted_elements() {
         let mut s = Session::default();

crates/engine/src/cmd/mod.rs

@@ -20,6 +20,7 @@ pub static ALL: &[CommandSpec] = &[
     mesh::ELEMENT_ADD,
     mesh::ELEMENT_REMOVE,
     mesh::QUALITY,
+    mesh::DELAUNAY_POINTS,
 ];
 
 /// A required label parameter (1..=u32::MAX; the model checks its own range).
@@ -32,13 +33,22 @@ pub(crate) fn label(p: &Value, key: &str) -> Result<u32> {
 
 /// A required `[x, y, z]` parameter.
 pub(crate) fn xyz(p: &Value, key: &str) -> Result<[f64; 3]> {
-    let bad = || EngineError::Params(format!("`{key}` must be [x, y, z]"));
+    p.get(key).and_then(xyz_value).ok_or_else(|| EngineError::Params(format!("`{key}` must be [x, y, z]")))
-    let a = p.get(key).and_then(Value::as_array).filter(|a| a.len() == 3).ok_or_else(bad)?;
+}
-    let mut out = [0.0; 3];
+
-    for (o, v) in out.iter_mut().zip(a) {
+/// A JSON `[x, y, z]`.
-        *o = v.as_f64().ok_or_else(bad)?;
+fn xyz_value(v: &Value) -> Option<[f64; 3]> {
+    match v.as_array()?.as_slice() {
+        [x, y, z] => Some([x.as_f64()?, y.as_f64()?, z.as_f64()?]),
+        _ => None,
     }
-    Ok(out)
+}
+
+/// A required list of `[x, y, z]` points, at most `max_len` long.
+pub(crate) fn points(p: &Value, key: &str, max_len: usize) -> Result<Vec<[f64; 3]>> {
+    let bad = |what: String| EngineError::Params(format!("`{key}` must be a list of at most {max_len} [x, y, z]{what}"));
+    let raw = p.get(key).and_then(Value::as_array).filter(|a| a.len() <= max_len).ok_or_else(|| bad(String::new()))?;
+    raw.iter().enumerate().map(|(i, v)| xyz_value(v).ok_or_else(|| bad(format!(" (entry {i} is not)")))).collect()
 }
 
 /// A required list of labels (each in 1..=MAX_LABEL), at most `max_len` long.

crates/engine/src/lib.rs

@@ -179,6 +179,11 @@ mod tests {
             json!({"elements": [u32::MAX, -1, 1.5]}),
             json!({"elements": {}, "bins": -5}),
             json!({"elements": [1], "bins": 1e300}),
+            json!({"points": null}),
+            json!({"points": [[1, 2]]}),
+            json!({"points": [[1e308, -1e308, 0], [0, 0, 0]]}),
+            json!({"points": [[0, 0, 0], [0, 0, 0], [0, 0, 0], [0, 0, 0], [0, 0, 0]]}),
+            json!({"points": vec![[1, 1, 1]; 30_000]}),
         ];
         for spec in commands() {
             for p in &hostile {

crates/mesh/src/delaunay/mod.rs

@@ -0,0 +1,689 @@
+//! Exact, deterministic 3D Delaunay tetrahedralization of a point set.
+//!
+//! # Method
+//! - **Incremental cavity insertion** (A. Bowyer, *Computing Dirichlet tessellations*, Comput. J.
+//!   24(2):162–166, 1981; D. F. Watson, *Computing the n-dimensional Delaunay tessellation with
+//!   application to Voronoi polytopes*, Comput. J. 24(2):167–172, 1981). Each point finds the
+//!   tetrahedra whose circumsphere contains it (the cavity, grown by adjacency from one of them),
+//!   deletes them and joins the cavity's boundary faces to the point.
+//! - **Insertion order**: BRIO rounds with Hilbert order inside each round (Amenta, Choi, Rote,
+//!   SoCG 2003; Skilling 2004). See [`order`].
+//! - **Point location**: the stochastic visibility walk with memory (O. Devillers, S. Pion,
+//!   M. Teillaud, *Walking in a triangulation*, IJFCS 13(2):181–199, 2002, §2–§3), starting at
+//!   the last tetrahedron created. In a Delaunay triangulation the walk cannot cycle
+//!   (H. Edelsbrunner, *An acyclicity theorem for cell complexes in d dimensions*, Combinatorica
+//!   10(3):251–260, 1990); it is still capped at the number of tetrahedra, after which a linear
+//!   scan locates the point.
+//! - **The infinite vertex.** The convex hull is closed by *ghost* tetrahedra: each hull face is
+//!   joined to the vertex [`GHOST`], so every face has exactly two tetrahedra and the adjacency is
+//!   that of a triangulated 3-sphere (no boundary cases). A ghost tetrahedron's "circumsphere" is
+//!   the open half-space beyond its hull face; a point on the plane of the face is decided by the
+//!   finite tetrahedron on the other side of that face (the limit of spheres through the face
+//!   whose centre runs off to infinity), so the decision is the same perturbed one
+//!   (O. Devillers, M. Teillaud, *Perturbations for Delaunay and weighted Delaunay 3D
+//!   triangulations*, CGTA 44(3):160–168, 2011, §3 and the infinite-vertex discussion).
+//! - **Compact adjacency arrays**: per tetrahedron four vertex ids and four neighbour ids, where
+//!   neighbour `i` is across the face opposite vertex `i`. Deleted slots are recycled LIFO.
+//!
+//! # Exactness and degeneracies
+//! Every decision is an exact predicate from [`crate::predicates`]: `orient3d` for the walk and
+//! for the hull, `insphere_sos` (Devillers–Teillaud symbolic perturbation) for the cavity. There
+//! are no tolerances. The perturbation is a function of the coordinates only (it lifts each point
+//! by an infinitesimal that depends on its lexicographic rank), so the result is *the* Delaunay
+//! triangulation of the perturbed points: unique, independent of the insertion order, and with
+//! no flat tetrahedra, even on grids where every eight neighbours are cospherical.
+//!
+//! - **Duplicates are merged.** Points with equal coordinates (`-0.0 == 0.0`) become one vertex:
+//!   the one with the lowest input index. The others are reported by
+//!   [`Triangulation::duplicate_of`] and appear in no tetrahedron.
+//! - **Lower-dimensional inputs** (all points collinear or coplanar, a single point, none) have
+//!   no tetrahedra: [`Triangulation::dimension`] reports their affine dimension and the
+//!   tetrahedron lists are empty. Points arriving later that leave the plane start the 3D
+//!   triangulation.
+//! - **Hostile coordinates** are rejected at the boundary by [`Point3::new`] (NaN, infinities,
+//!   magnitudes outside `[2⁻¹²⁰, 2¹⁶⁰]`); at most [`MAX_POINTS`] points.
+//!
+//! # Determinism
+//! Same points in the same order → same tetrahedra in the same slots, bit for bit: the order is
+//! a fixed function of the input, the walk's "random" face choice comes from a fixed-seed
+//! generator, sorting has total keys, and nothing iterates a hash map.
+//!
+//! # Verification
+//! [`Triangulation::verify`] checks the adjacency invariants, the orientation of every
+//! tetrahedron, the convexity of the hull and the empty-circumsphere property of every interior
+//! face (exact and perturbed), which by Delaunay's lemma makes every circumsphere empty. Tests
+//! also use [`Triangulation::verify_empty_spheres`], the brute-force check of every point against
+//! every circumsphere.
+
+mod order;
+mod verify;
+
+#[cfg(test)]
+mod tests;
+
+use crate::predicates::{Axis, GenericPoint, Point3, Sign, generic_orient2d, insphere_sos, orient3d, tet_volume_sign};
+use crate::{Error, Result};
+
+/// The infinite vertex of the ghost tetrahedra that close the convex hull.
+pub const GHOST: u32 = u32::MAX;
+/// Marks a deleted tetrahedron slot (first vertex).
+const DEAD: u32 = u32::MAX - 1;
+/// No tetrahedron / no vertex.
+const NONE: u32 = u32::MAX;
+/// An empty slot of the edge table (no edge has two equal endpoints).
+const EMPTY_EDGE: u64 = u64::MAX;
+/// The most points a triangulation takes (vertex ids stay far below [`GHOST`] and the ~6.7
+/// tetrahedra per point fit `u32` slot ids).
+pub const MAX_POINTS: usize = 1 << 24;
+
+/// Where a point lies relative to the triangulation ([`Triangulation::locate`]).
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub enum Location {
+    /// In the closed finite tetrahedron `t` and not one of its vertices.
+    Inside(u32),
+    /// Strictly beyond the hull face of the ghost tetrahedron `t`.
+    Outside(u32),
+    /// Equal to this vertex.
+    Vertex(u32),
+}
+
+/// The outcome of [`Triangulation::insert`].
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub enum Insertion {
+    /// A new vertex.
+    Inserted(u32),
+    /// The point equals this existing vertex; it was recorded as a duplicate (its own id is the
+    /// last point index).
+    Duplicate(u32),
+}
+
+/// A Delaunay tetrahedralization with ghost tetrahedra closing the hull.
+///
+/// Vertex ids are input point indices. Tetrahedron ids are slots: [`Triangulation::slots`]
+/// counts them, and [`Triangulation::tet`] returns `None` for a recycled slot. A tetrahedron's
+/// vertices `[a, b, c, d]` have positive volume, `det[b − a, c − a, d − a] > 0` (the tet4 node
+/// order). A ghost tetrahedron has [`GHOST`] in one slot; putting a point in place of the ghost
+/// gives positive volume exactly when that point is beyond the hull face.
+#[derive(Debug, Clone)]
+pub struct Triangulation {
+    points: Vec<Point3>,
+    /// `NONE`, or the vertex this point duplicates.
+    duplicate_of: Vec<u32>,
+    tets: Vec<[u32; 4]>,
+    /// Neighbour `i` is across the face opposite vertex `i`.
+    adj: Vec<[u32; 4]>,
+    free: Vec<u32>,
+    /// A live tetrahedron incident to each vertex (`NONE` until it is in one).
+    vertex_tet: Vec<u32>,
+    /// While the points span less than 3D: affinely independent vertices, and the others.
+    basis: Vec<u32>,
+    pending: Vec<u32>,
+    /// Cavity marks: `epoch` = in the cavity, `epoch + 1` = tested and outside.
+    mark: Vec<u32>,
+    epoch: u32,
+    /// Where the next walk starts.
+    hint: u32,
+    /// State of the walk's face-order generator.
+    rng: u64,
+    scratch: Scratch,
+}
+
+#[derive(Debug, Clone, Default)]
+struct Scratch {
+    cavity: Vec<u32>,
+    stack: Vec<u32>,
+    /// `(cavity tet, face slot, outside neighbour)`.
+    boundary: Vec<(u32, usize, u32)>,
+    /// `(new tet vertices, slot of the outside neighbour, neighbour, the neighbour's slot back)`.
+    new_tets: Vec<([u32; 4], usize, u32, usize)>,
+    /// Open-addressing table `(edge key, tet, slot)` linking the new tetrahedra to each other.
+    edges: Vec<(u64, u32, usize)>,
+}
+
+/// A Delaunay triangulation is corrupt: this never happens with exact predicates, and is
+/// reported instead of panicking.
+fn corrupt(what: &str) -> Error {
+    Error::Failed(format!("Delaunay triangulation is inconsistent: {what}"))
+}
+
+/// The sign of the volume `det[b − a, c − a, d − a]` (the tet4 convention).
+#[inline]
+fn volume_sign(p: [&Point3; 4]) -> Sign {
+    let [a, b, c, d] = p;
+    tet_volume_sign(a, b, c, d)
+}
+
+/// Whether three points are collinear: all three axis projections are (exactly).
+fn collinear(a: &Point3, b: &Point3, c: &Point3) -> bool {
+    let [a, b, c] = [a, b, c].map(|p| GenericPoint::from(*p));
+    [Axis::X, Axis::Y, Axis::Z].into_iter().all(|axis| generic_orient2d(&a, &b, &c, axis).is_zero())
+}
+
+/// A coordinate triple as a [`Point3`], with `-0.0` made `+0.0` so equal points have equal bits.
+fn admit([x, y, z]: [f64; 3]) -> Result<Point3> {
+    Point3::new(x + 0.0, y + 0.0, z + 0.0)
+}
+
+impl Triangulation {
+    fn empty() -> Self {
+        Triangulation {
+            points: Vec::new(),
+            duplicate_of: Vec::new(),
+            tets: Vec::new(),
+            adj: Vec::new(),
+            free: Vec::new(),
+            vertex_tet: Vec::new(),
+            basis: Vec::new(),
+            pending: Vec::new(),
+            mark: Vec::new(),
+            epoch: 0,
+            hint: NONE,
+            rng: 0x2545_F491_4F6C_DD1D,
+            scratch: Scratch::default(),
+        }
+    }
+
+    /// The Delaunay tetrahedralization of `points`. Duplicates are merged into the lowest
+    /// index; fewer than four affinely independent points give no tetrahedra.
+    pub fn new(points: &[[f64; 3]]) -> Result<Self> {
+        if points.len() > MAX_POINTS {
+            return Err(Error::Limit(format!("{} points; a triangulation takes at most {MAX_POINTS}", points.len())));
+        }
+        let mut tri = Self::empty();
+        tri.points =
+            points.iter().enumerate().map(|(i, &c)| admit(c).map_err(|e| Error::Invalid(format!("point {i}: {e}")))).collect::<Result<_>>()?;
+        let n = tri.points.len();
+        tri.duplicate_of = vec![NONE; n];
+        tri.vertex_tet = vec![NONE; n];
+
+        // Merge duplicates: sort by coordinates, ties by index, keep the first of each run.
+        let mut by_coords: Vec<u32> = (0..n as u32).collect();
+        let key = |v: u32| tri.points.get(v as usize).map(|p| p.coords().map(f64::to_bits)).unwrap_or_default();
+        // Equal values have equal bits after `admit`; the bit order is total and deterministic.
+        by_coords.sort_unstable_by_key(|&v| (key(v), v));
+        for run in by_coords.chunk_by(|&a, &b| key(a) == key(b)) {
+            if let [keep, rest @ ..] = run {
+                for &v in rest {
+                    if let Some(d) = tri.duplicate_of.get_mut(v as usize) {
+                        *d = *keep;
+                    }
+                }
+            }
+        }
+        let unique: Vec<u32> = (0..n as u32).filter(|&v| tri.duplicate_of.get(v as usize) == Some(&NONE)).collect();
+
+        let expected = unique.len().saturating_mul(7);
+        tri.tets.reserve(expected);
+        tri.adj.reserve(expected);
+        tri.mark.reserve(expected);
+        for v in order::brio_hilbert(&tri.points, &unique) {
+            tri.add(v, true)?;
+        }
+        Ok(tri)
+    }
+
+    /// Insert one more point (a Steiner point, say). Its vertex id is its index in
+    /// [`Triangulation::points`]. A point equal to an existing vertex is recorded as a duplicate.
+    pub fn insert(&mut self, p: [f64; 3]) -> Result<Insertion> {
+        if self.points.len() >= MAX_POINTS {
+            return Err(Error::Limit(format!("a triangulation takes at most {MAX_POINTS} points")));
+        }
+        let v = self.points.len() as u32;
+        self.points.push(admit(p)?);
+        self.duplicate_of.push(NONE);
+        self.vertex_tet.push(NONE);
+        self.add(v, false)
+    }
+
+    /// The points, by vertex id (including duplicates).
+    pub fn points(&self) -> &[Point3] {
+        &self.points
+    }
+
+    /// The vertex a duplicate point was merged into.
+    pub fn duplicate_of(&self, v: u32) -> Option<u32> {
+        self.duplicate_of.get(v as usize).copied().filter(|&d| d != NONE)
+    }
+
+    /// The affine dimension of the distinct points: 3 when there are tetrahedra, otherwise 0
+    /// (none or one point), 1 (collinear) or 2 (coplanar).
+    pub fn dimension(&self) -> usize {
+        if self.tets.is_empty() { self.basis.len().saturating_sub(1) } else { 3 }
+    }
+
+    /// Tetrahedron slots, live or recycled, finite or ghost.
+    pub fn slots(&self) -> usize {
+        self.tets.len()
+    }
+
+    /// The vertices of the tetrahedron in slot `t` (one may be [`GHOST`]); `None` for a recycled
+    /// or out-of-range slot.
+    pub fn tet(&self, t: u32) -> Option<[u32; 4]> {
+        self.tets.get(t as usize).copied().filter(|v| v[0] != DEAD)
+    }
+
+    /// The neighbours of tetrahedron `t`: neighbour `i` shares the face opposite vertex `i`.
+    pub fn neighbors(&self, t: u32) -> Option<[u32; 4]> {
+        self.tet(t)?;
+        self.adj.get(t as usize).copied()
+    }
+
+    /// Whether `t` is a ghost tetrahedron (one vertex is [`GHOST`]).
+    pub fn is_ghost(&self, t: u32) -> bool {
+        self.tet(t).is_some_and(|v| v.contains(&GHOST))
+    }
+
+    /// A live tetrahedron incident to vertex `v` (finite or ghost), once `v` is triangulated.
+    pub fn incident_tet(&self, v: u32) -> Option<u32> {
+        self.vertex_tet.get(v as usize).copied().filter(|&t| t != NONE)
+    }
+
+    /// The finite tetrahedra with their slot ids, in slot order; vertices in tet4 order.
+    pub fn finite_tets(&self) -> impl Iterator<Item = (u32, [u32; 4])> + '_ {
+        self.tets.iter().enumerate().filter(|(_, v)| v[0] != DEAD && !v.contains(&GHOST)).map(|(t, v)| (t as u32, *v))
+    }
+
+    /// The number of finite tetrahedra.
+    pub fn finite_tet_count(&self) -> usize {
+        self.finite_tets().count()
+    }
+
+    /// Locate `p`: the walk from the last tetrahedron created.
+    pub fn locate(&self, p: [f64; 3]) -> Result<Option<Location>> {
+        if self.tets.is_empty() {
+            return Ok(None);
+        }
+        let point = admit(p)?;
+        let mut rng = self.rng;
+        self.walk(&point, self.hint, &mut rng).map(Some)
+    }
+
+    // ---- internals ----
+
+    fn pt(&self, v: u32) -> Result<&Point3> {
+        self.points.get(v as usize).ok_or_else(|| corrupt("vertex id"))
+    }
+
+    fn tet_at(&self, t: u32) -> Result<[u32; 4]> {
+        self.tet(t).ok_or_else(|| corrupt("tetrahedron id"))
+    }
+
+    fn adj_at(&self, t: u32) -> Result<[u32; 4]> {
+        self.adj.get(t as usize).copied().ok_or_else(|| corrupt("tetrahedron id"))
+    }
+
+    /// The four points of a finite tetrahedron.
+    fn pts(&self, [a, b, c, d]: [u32; 4]) -> Result<[&Point3; 4]> {
+        Ok([self.pt(a)?, self.pt(b)?, self.pt(c)?, self.pt(d)?])
+    }
+
+    /// The four points of `vs` with slot `k` replaced by `p` (`k` may hold the ghost).
+    fn with_point<'a>(&'a self, vs: [u32; 4], k: usize, p: &'a Point3) -> Result<[&'a Point3; 4]> {
+        let mut out = [p; 4];
+        for (i, (o, &v)) in out.iter_mut().zip(&vs).enumerate() {
+            if i != k {
+                *o = self.pt(v)?;
+            }
+        }
+        Ok(out)
+    }
+
+    /// Add vertex `v` (already in `points`). A duplicate is detected and recorded: in 3D by the
+    /// walk, below 3D by comparing with the vertices so far unless the caller knows `v` is
+    /// `unique` (the sort-based merge of [`Triangulation::new`]).
+    fn add(&mut self, v: u32, unique: bool) -> Result<Insertion> {
+        if !self.tets.is_empty() {
+            return self.insert_vertex(v);
+        }
+        let p = *self.pt(v)?;
+        if !unique && let Some(&u) = self.basis.iter().chain(&self.pending).find(|&&u| self.points.get(u as usize) == Some(&p)) {
+            if let Some(d) = self.duplicate_of.get_mut(v as usize) {
+                *d = u;
+            }
+            return Ok(Insertion::Duplicate(u));
+        }
+        match self.basis.as_slice() {
+            [] | [_] => self.basis.push(v),
+            &[a, b] => {
+                if collinear(self.pt(a)?, self.pt(b)?, &p) {
+                    self.pending.push(v);
+                } else {
+                    self.basis.push(v);
+                }
+            }
+            &[a, b, c, ..] => {
+                if orient3d(self.pt(a)?, self.pt(b)?, self.pt(c)?, &p).is_zero() {
+                    self.pending.push(v);
+                } else {
+                    self.start([a, b, c, v])?;
+                    for u in std::mem::take(&mut self.pending) {
+                        self.insert_vertex(u)?;
+                    }
+                }
+            }
+        }
+        Ok(Insertion::Inserted(v))
+    }
+
+    /// The first tetrahedron and its four ghosts.
+    fn start(&mut self, [a, b, c, d]: [u32; 4]) -> Result<()> {
+        let t0 = if volume_sign(self.pts([a, b, c, d])?) == Sign::Positive { [a, b, c, d] } else { [b, a, c, d] };
+        let mut all = vec![t0];
+        for j in 0..4 {
+            let mut g = t0;
+            g[j] = GHOST;
+            // Swapping two other slots turns the ghost to face away from the hull.
+            g.swap((j + 1) % 4, (j + 2) % 4);
+            all.push(g);
+        }
+        let mut faces: Vec<([u32; 3], u32, usize)> = Vec::with_capacity(20);
+        for vs in all {
+            let t = self.alloc(vs)?;
+            for k in 0..4 {
+                faces.push((verify::face(vs, k), t, k));
+            }
+            self.hint = t;
+        }
+        faces.sort_unstable();
+        for pair in faces.chunks(2) {
+            match pair {
+                [(f, t, i), (g, u, j)] if f == g => {
+                    self.set_adj(*t, *i, *u)?;
+                    self.set_adj(*u, *j, *t)?;
+                }
+                _ => return Err(corrupt("initial tetrahedron faces")),
+            }
+        }
+        self.basis.clear();
+        Ok(())
+    }
+
+    fn set_adj(&mut self, t: u32, i: usize, n: u32) -> Result<()> {
+        let slot = self.adj.get_mut(t as usize).and_then(|a| a.get_mut(i)).ok_or_else(|| corrupt("tetrahedron id"))?;
+        *slot = n;
+        Ok(())
+    }
+
+    fn alloc(&mut self, vs: [u32; 4]) -> Result<u32> {
+        let t = match self.free.pop() {
+            Some(t) => {
+                *self.tets.get_mut(t as usize).ok_or_else(|| corrupt("free slot"))? = vs;
+                *self.adj.get_mut(t as usize).ok_or_else(|| corrupt("free slot"))? = [NONE; 4];
+                t
+            }
+            None => {
+                if self.tets.len() >= (DEAD - 1) as usize {
+                    return Err(Error::Limit("too many tetrahedra".into()));
+                }
+                self.tets.push(vs);
+                self.adj.push([NONE; 4]);
+                self.mark.push(0);
+                (self.tets.len() - 1) as u32
+            }
+        };
+        for v in vs {
+            if let Some(s) = self.vertex_tet.get_mut(v as usize) {
+                *s = t;
+            }
+        }
+        Ok(t)
+    }
+
+    fn kill(&mut self, t: u32) -> Result<()> {
+        *self.tets.get_mut(t as usize).ok_or_else(|| corrupt("tetrahedron id"))? = [DEAD; 4];
+        self.free.push(t);
+        Ok(())
+    }
+
+    /// Whether `p` is in conflict with tetrahedron `t` (inside its perturbed circumsphere). Only
+    /// `dig` calls it: it reads the current cavity marks.
+    fn conflict(&self, t: u32, p: &Point3) -> Result<bool> {
+        let vs = self.tet_at(t)?;
+        let Some(g) = vs.iter().position(|&v| v == GHOST) else {
+            let [a, b, c, d] = self.pts(vs)?;
+            // `insphere` wants Shewchuk-positive orientation: swap the first two vertices.
+            return Ok(insphere_sos(b, a, c, d, p) == Sign::Positive);
+        };
+        Ok(match volume_sign(self.with_point(vs, g, p)?) {
+            Sign::Positive => true,
+            Sign::Negative => false,
+            // On the hull face's plane: the finite tetrahedron across that face decides.
+            Sign::Zero => {
+                let o = self.adj_at(t)?.get(g).copied().ok_or_else(|| corrupt("slot"))?;
+                // During `dig` the finite tetrahedron may already be decided.
+                match self.mark_of(o) {
+                    m if m == self.epoch => return Ok(true),
+                    m if m == self.epoch + 1 => return Ok(false),
+                    _ => {}
+                }
+                let ws = self.tet_at(o)?;
+                if ws.contains(&GHOST) {
+                    return Err(corrupt("two ghost tetrahedra share a hull face"));
+                }
+                let [a, b, c, d] = self.pts(ws)?;
+                insphere_sos(b, a, c, d, p) == Sign::Positive
+            }
+        })
+    }
+
+    /// The visibility walk towards `p` from `start` (Devillers, Pion, Teillaud 2002).
+    fn walk(&self, p: &Point3, start: u32, rng: &mut u64) -> Result<Location> {
+        let mut t = start;
+        let vs = self.tet_at(t)?;
+        if let Some(g) = vs.iter().position(|&v| v == GHOST) {
+            t = self.adj_at(t)?.get(g).copied().ok_or_else(|| corrupt("slot"))?;
+        }
+        let mut prev = NONE;
+        for _ in 0..self.tets.len().saturating_add(16) {
+            let vs = self.tet_at(t)?;
+            if vs.contains(&GHOST) {
+                return Ok(Location::Outside(t));
+            }
+            let adj = self.adj_at(t)?;
+            let ps = self.pts(vs)?;
+            // xorshift64 (Marsaglia 2003): which face to try first.
+            *rng ^= *rng << 13;
+            *rng ^= *rng >> 7;
+            *rng ^= *rng << 17;
+            let first = (*rng >> 61) as usize;
+            let mut next = NONE;
+            for i in 0..4 {
+                let k = (first + i) % 4;
+                let n = adj.get(k).copied().unwrap_or(NONE);
+                if n == prev {
+                    continue;
+                }
+                let mut q = ps;
+                if let Some(slot) = q.get_mut(k) {
+                    *slot = p;
+                }
+                if volume_sign(q) == Sign::Negative {
+                    next = n;
+                    break;
+                }
+            }
+            if next == NONE {
+                return self.inside(vs, t, p);
+            }
+            prev = t;
+            t = next;
+        }
+        self.locate_by_scan(p)
+    }
+
+    /// `p` lies in the closed finite tetrahedron `t`: one of its vertices, or inside.
+    fn inside(&self, vs: [u32; 4], t: u32, p: &Point3) -> Result<Location> {
+        for v in vs {
+            if self.pt(v)? == p {
+                return Ok(Location::Vertex(v));
+            }
+        }
+        Ok(Location::Inside(t))
+    }
+
+    /// Point location by testing every tetrahedron (the walk's fallback).
+    fn locate_by_scan(&self, p: &Point3) -> Result<Location> {
+        let mut outside = None;
+        for (t, &vs) in self.tets.iter().enumerate() {
+            if vs[0] == DEAD {
+                continue;
+            }
+            match vs.iter().position(|&v| v == GHOST) {
+                Some(g) => {
+                    if outside.is_none() && volume_sign(self.with_point(vs, g, p)?) == Sign::Positive {
+                        outside = Some(t as u32);
+                    }
+                }
+                None => {
+                    let mut inside = true;
+                    for k in 0..4 {
+                        if volume_sign(self.with_point(vs, k, p)?) == Sign::Negative {
+                            inside = false;
+                            break;
+                        }
+                    }
+                    if inside {
+                        return self.inside(vs, t as u32, p);
+                    }
+                }
+            }
+        }
+        outside.map(Location::Outside).ok_or_else(|| corrupt("a point is neither inside nor outside the hull"))
+    }
+
+    fn next_epoch(&mut self) {
+        if self.epoch >= u32::MAX - 3 {
+            self.mark.fill(0);
+            self.epoch = 0;
+        }
+        self.epoch += 2;
+    }
+
+    fn mark_of(&self, t: u32) -> u32 {
+        self.mark.get(t as usize).copied().unwrap_or(0)
+    }
+
+    fn set_mark(&mut self, t: u32, m: u32) {
+        if let Some(s) = self.mark.get_mut(t as usize) {
+            *s = m;
+        }
+    }
+
+    /// Insert vertex `v` into the 3D triangulation.
+    fn insert_vertex(&mut self, v: u32) -> Result<Insertion> {
+        let p = *self.pt(v)?;
+        let mut rng = self.rng;
+        let located = self.walk(&p, self.hint, &mut rng)?;
+        self.rng = rng;
+        let seed = match located {
+            Location::Vertex(u) => {
+                if let Some(d) = self.duplicate_of.get_mut(v as usize) {
+                    *d = u;
+                }
+                return Ok(Insertion::Duplicate(u));
+            }
+            Location::Inside(t) | Location::Outside(t) => t,
+        };
+        let mut s = std::mem::take(&mut self.scratch);
+        let result = self.dig(seed, &p, &mut s).and_then(|()| self.fill(v, &mut s));
+        self.scratch = s;
+        result.map(|()| Insertion::Inserted(v))
+    }
+
+    /// The cavity of `p`: every tetrahedron in conflict, grown by adjacency from `seed`, and the
+    /// faces between it and the rest.
+    fn dig(&mut self, seed: u32, p: &Point3, s: &mut Scratch) -> Result<()> {
+        self.next_epoch();
+        let (inside, outside) = (self.epoch, self.epoch + 1);
+        s.cavity.clear();
+        s.stack.clear();
+        s.boundary.clear();
+        self.set_mark(seed, inside);
+        s.cavity.push(seed);
+        s.stack.push(seed);
+        while let Some(t) = s.stack.pop() {
+            for (j, o) in self.adj_at(t)?.into_iter().enumerate() {
+                let m = self.mark_of(o);
+                if m == inside {
+                    continue;
+                }
+                if m != outside && self.conflict(o, p)? {
+                    self.set_mark(o, inside);
+                    s.cavity.push(o);
+                    s.stack.push(o);
+                } else {
+                    self.set_mark(o, outside);
+                    s.boundary.push((t, j, o));
+                }
+            }
+        }
+        Ok(())
+    }
+
+    /// Replace the cavity by the tetrahedra joining its boundary faces to `v`.
+    fn fill(&mut self, v: u32, s: &mut Scratch) -> Result<()> {
+        s.new_tets.clear();
+        // Everything is read before any slot is recycled: a new tetrahedron may reuse the id of
+        // a cavity tetrahedron that another outside neighbour still points to.
+        for &(t, j, o) in &s.boundary {
+            let mut vs = self.tet_at(t)?;
+            *vs.get_mut(j).ok_or_else(|| corrupt("slot"))? = v;
+            let back = self.adj_at(o)?.iter().position(|&n| n == t).ok_or_else(|| corrupt("one-sided adjacency"))?;
+            s.new_tets.push((vs, j, o, back));
+        }
+        for &t in &s.cavity {
+            self.kill(t)?;
+        }
+        // Each edge of the cavity boundary is shared by exactly two boundary faces; their new
+        // tetrahedra meet across the face spanned by that edge and `v`. Pair them through a
+        // small open-addressing table keyed by the edge (lookups only, so no order leaks out).
+        // 1.5 edges per new tetrahedron: at most 75% load.
+        let size = (s.new_tets.len() * 2).next_power_of_two().max(16);
+        s.edges.clear();
+        s.edges.resize(size, (EMPTY_EDGE, NONE, 0));
+        let mask = size - 1;
+        let mut open = 0usize;
+        for &(vs, j, o, back) in &s.new_tets {
+            let t = self.alloc(vs)?;
+            self.set_adj(t, j, o)?;
+            self.set_adj(o, back, t)?;
+            for k in (0..4).filter(|&k| k != j) {
+                let mut e = vs.iter().enumerate().filter(|&(i, _)| i != j && i != k).map(|(_, &u)| u);
+                let (a, b) = (e.next().unwrap_or(NONE), e.next().unwrap_or(NONE));
+                let key = (u64::from(a.min(b)) << 32) | u64::from(a.max(b));
+                let mut h = (key.wrapping_mul(0x9E37_79B9_7F4A_7C15) >> 32) as usize & mask;
+                let mut placed = false;
+                for _ in 0..size {
+                    let entry = s.edges.get_mut(h).ok_or_else(|| corrupt("edge table"))?;
+                    if entry.0 == EMPTY_EDGE {
+                        *entry = (key, t, k);
+                        open += 1;
+                        placed = true;
+                        break;
+                    }
+                    if entry.0 == key {
+                        let (_, u, l) = *entry;
+                        if u == NONE {
+                            return Err(corrupt("an edge of the cavity boundary has more than two faces"));
+                        }
+                        entry.1 = NONE;
+                        open -= 1;
+                        self.set_adj(t, k, u)?;
+                        self.set_adj(u, l, t)?;
+                        placed = true;
+                        break;
+                    }
+                    h = (h + 1) & mask;
+                }
+                if !placed {
+                    return Err(corrupt("edge table full"));
+                }
+            }
+            self.hint = t;
+        }
+        if open != 0 {
+            return Err(corrupt("the cavity boundary is not a closed surface"));
+        }
+        Ok(())
+    }
+}

crates/mesh/src/delaunay/order.rs

@@ -0,0 +1,138 @@
+//! Insertion order: biased randomized insertion order (BRIO) with Hilbert-curve order inside
+//! each round.
+//!
+//! - N. Amenta, S. Choi, G. Rote, *Incremental constructions con BRIO*, SoCG 2003, 211–219:
+//!   points go into rounds of geometrically growing size (each point independently lands in round
+//!   `k` from the end with probability 2⁻⁽ᵏ⁺¹⁾). Randomness between rounds keeps the expected
+//!   cost of a randomized incremental construction; order inside a round is free.
+//! - Inside a round the points follow a Hilbert curve, so consecutive points are close and the
+//!   point-location walk from the previous insertion is short. The index is computed with
+//!   J. Skilling, *Programming the Hilbert curve*, AIP Conf. Proc. 707:381–387, 2004
+//!   (axes → transposed index, then bit interleaving).
+//!
+//! The "random" round of a point is a fixed hash of its index, so the order, and with it every
+//! later step, is a deterministic function of the input. Floating point is used only to quantize
+//! coordinates for the curve; it never decides any topology.
+
+use crate::predicates::Point3;
+
+/// Bits per axis of the Hilbert grid (3 × 21 = 63-bit keys).
+const BITS: u32 = 21;
+/// The deepest BRIO round; deeper hashes are clamped (2⁻²⁰ of the points share it).
+const MAX_LEVEL: u32 = 20;
+
+/// SplitMix64 finalizer (Steele, Lea, Flood, *Fast splittable pseudorandom number generators*,
+/// OOPSLA 2014): a fixed, well-mixed hash of the vertex index.
+pub(super) fn mix(mut z: u64) -> u64 {
+    z = z.wrapping_add(0x9E37_79B9_7F4A_7C15);
+    z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
+    z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
+    z ^ (z >> 31)
+}
+
+/// Hilbert index of a point on the `2^BITS`-per-axis grid (Skilling 2004).
+fn hilbert(mut x: [u32; 3]) -> u64 {
+    let m = 1u32 << (BITS - 1);
+    // Inverse undo excess work.
+    let mut q = m;
+    while q > 1 {
+        let p = q - 1;
+        for i in 0..3 {
+            let xi = x[i];
+            if xi & q != 0 {
+                x[0] ^= p;
+            } else {
+                let t = (x[0] ^ xi) & p;
+                x[0] ^= t;
+                x[i] = xi ^ t;
+            }
+        }
+        q >>= 1;
+    }
+    // Gray encode.
+    x[1] ^= x[0];
+    x[2] ^= x[1];
+    let mut t = 0;
+    q = m;
+    while q > 1 {
+        if x[2] & q != 0 {
+            t ^= q - 1;
+        }
+        q >>= 1;
+    }
+    for xi in &mut x {
+        *xi ^= t;
+    }
+    // Interleave the transposed bits, most significant first.
+    let mut key = 0u64;
+    for b in (0..BITS).rev() {
+        for xi in x {
+            key = (key << 1) | u64::from((xi >> b) & 1);
+        }
+    }
+    key
+}
+
+/// The vertices in BRIO order, Hilbert order inside each round, ties by vertex index.
+pub(super) fn brio_hilbert(points: &[Point3], vertices: &[u32]) -> Vec<u32> {
+    let coords = |v: u32| points.get(v as usize).map(Point3::coords);
+    let (mut lo, mut hi) = ([f64::INFINITY; 3], [f64::NEG_INFINITY; 3]);
+    for c in vertices.iter().filter_map(|&v| coords(v)) {
+        for k in 0..3 {
+            lo[k] = lo[k].min(c[k]);
+            hi[k] = hi[k].max(c[k]);
+        }
+    }
+    let extent = (0..3).map(|k| hi[k] - lo[k]).fold(0.0f64, f64::max);
+    let top = f64::from((1u32 << BITS) - 1);
+    let scale = if extent > 0.0 && extent.is_finite() { top / extent } else { 0.0 };
+    let mut keyed: Vec<(u32, u64, u32)> = vertices
+        .iter()
+        .map(|&v| {
+            let c = coords(v).unwrap_or_default();
+            // `as` saturates; the clamp guards the last rounding step.
+            let q = std::array::from_fn(|k| (((c[k] - lo[k]) * scale) as u32).min((1 << BITS) - 1));
+            let level = mix(u64::from(v)).trailing_ones().min(MAX_LEVEL);
+            (MAX_LEVEL - level, hilbert(q), v)
+        })
+        .collect();
+    keyed.sort_unstable();
+    keyed.into_iter().map(|(_, _, v)| v).collect()
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn hilbert_visits_neighbouring_cells_in_order() {
+        // On a 2^BITS grid the curve over the first 8³ block is a walk of unit steps.
+        let mut cells: Vec<(u64, [u32; 3])> = Vec::new();
+        for x in 0..8 {
+            for y in 0..8 {
+                for z in 0..8 {
+                    cells.push((hilbert([x, y, z]), [x, y, z]));
+                }
+            }
+        }
+        cells.sort_unstable();
+        for w in cells.windows(2) {
+            assert_eq!(w[1].0, w[0].0 + 1, "the 8³ block is a prefix of the curve");
+            let d: u32 = (0..3).map(|k| w[0].1[k].abs_diff(w[1].1[k])).sum();
+            assert_eq!(d, 1, "{:?} -> {:?}", w[0].1, w[1].1);
+        }
+    }
+
+    #[test]
+    fn order_is_a_permutation_with_geometric_rounds() {
+        let points: Vec<Point3> = (0..10_000).map(|i| Point3::new(f64::from(i % 97), f64::from(i % 89), f64::from(i % 83)).unwrap()).collect();
+        let vertices: Vec<u32> = (0..10_000).collect();
+        let order = brio_hilbert(&points, &vertices);
+        let mut sorted = order.clone();
+        sorted.sort_unstable();
+        assert_eq!(sorted, vertices);
+        let last_round = vertices.iter().filter(|&&v| mix(u64::from(v)).trailing_ones() == 0).count();
+        assert!((4_500..5_500).contains(&last_round), "{last_round}");
+        assert_eq!(order, brio_hilbert(&points, &vertices), "deterministic");
+    }
+}

crates/mesh/src/delaunay/tests.rs

@@ -0,0 +1,326 @@
+use modelcraft_model::ElementKind;
+
+use super::*;
+use crate::quality::{GateOptions, enclosed_volume, validity_gate_elements};
+
+struct Rng(u64);
+
+impl Rng {
+    fn next_u64(&mut self) -> u64 {
+        self.0 = self.0.wrapping_add(1);
+        order::mix(self.0)
+    }
+    fn unit(&mut self) -> f64 {
+        (self.next_u64() >> 11) as f64 / (1u64 << 53) as f64
+    }
+    fn shuffle<T>(&mut self, v: &mut [T]) {
+        for i in (1..v.len()).rev() {
+            v.swap(i, (self.next_u64() % (i as u64 + 1)) as usize);
+        }
+    }
+}
+
+fn random(seed: u64, n: usize) -> Vec<[f64; 3]> {
+    let mut r = Rng(seed);
+    (0..n).map(|_| [r.unit(), r.unit(), r.unit()]).collect()
+}
+
+fn grid(nx: usize, ny: usize, nz: usize) -> Vec<[f64; 3]> {
+    let mut out = Vec::new();
+    for x in 0..nx {
+        for y in 0..ny {
+            for z in 0..nz {
+                out.push([x as f64, y as f64, z as f64]);
+            }
+        }
+    }
+    out
+}
+
+/// The finite tetrahedra as sorted coordinate quadruples: comparable across vertex numberings.
+fn tet_set(tri: &Triangulation) -> Vec<[[u64; 3]; 4]> {
+    let mut out: Vec<[[u64; 3]; 4]> = tri
+        .finite_tets()
+        .map(|(_, vs)| {
+            let mut q = vs.map(|v| tri.points()[v as usize].coords().map(f64::to_bits));
+            q.sort_unstable();
+            q
+        })
+        .collect();
+    out.sort_unstable();
+    out
+}
+
+/// The validity gate on the finite tetrahedra, with the hull volume expected; returns the volume.
+fn gate(tri: &Triangulation, expected: f64) -> f64 {
+    let elements = tri.finite_tets().map(|(t, vs)| (t + 1, ElementKind::Tet4, vs.map(|v| tri.points()[v as usize].coords())));
+    let report = validity_gate_elements(elements, &GateOptions { expected_volume: Some(expected), volume_tolerance: 1e-9 }).unwrap();
+    assert_eq!(report.checked, tri.finite_tet_count());
+    report.into_result().unwrap().volume
+}
+
+/// The volume enclosed by the hull faces, outward oriented, through the quality module.
+fn hull_volume(tri: &Triangulation) -> f64 {
+    let points: Vec<[f64; 3]> = tri.points().iter().map(Point3::coords).collect();
+    let mut triangles = Vec::new();
+    for t in 0..tri.slots() as u32 {
+        let Some(vs) = tri.tet(t) else { continue };
+        let Some(g) = vs.iter().position(|&v| v == GHOST) else { continue };
+        // Moving the ghost's slot to the end takes 3 - g transpositions.
+        let mut f: Vec<usize> = vs.iter().filter(|&&v| v != GHOST).map(|&v| v as usize).collect();
+        if (3 - g) % 2 == 1 {
+            f.swap(0, 1);
+        }
+        triangles.push([f[0], f[1], f[2]]);
+    }
+    enclosed_volume(&points, &triangles).unwrap()
+}
+
+/// Full check: invariants, brute-force empty spheres, the tet4 validity gate.
+fn check(tri: &Triangulation, hull_volume: f64) {
+    tri.verify().unwrap();
+    tri.verify_empty_spheres().unwrap();
+    gate(tri, hull_volume);
+}
+
+/// Delaunay tetrahedra by definition: every positively oriented 4-subset whose perturbed
+/// circumsphere excludes every other point.
+fn brute_force(points: &[[f64; 3]]) -> Vec<[[u64; 3]; 4]> {
+    let p: Vec<Point3> = points.iter().map(|&c| Point3::try_from(c).unwrap()).collect();
+    let n = p.len();
+    let mut out = Vec::new();
+    for i in 0..n {
+        for j in i + 1..n {
+            for k in j + 1..n {
+                for l in k + 1..n {
+                    let (a, b, c, d) = match orient3d(&p[i], &p[j], &p[k], &p[l]) {
+                        Sign::Zero => continue,
+                        Sign::Positive => (&p[i], &p[j], &p[k], &p[l]),
+                        Sign::Negative => (&p[j], &p[i], &p[k], &p[l]),
+                    };
+                    let empty = (0..n).filter(|m| ![i, j, k, l].contains(m)).all(|m| insphere_sos(a, b, c, d, &p[m]) == Sign::Negative);
+                    if empty {
+                        let mut q = [i, j, k, l].map(|v| points[v].map(f64::to_bits));
+                        q.sort_unstable();
+                        out.push(q);
+                    }
+                }
+            }
+        }
+    }
+    out.sort_unstable();
+    out
+}
+
+#[test]
+fn one_tetrahedron_and_its_ghosts() {
+    let tri = Triangulation::new(&[[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]).unwrap();
+    assert_eq!(tri.dimension(), 3);
+    assert_eq!(tri.finite_tet_count(), 1);
+    assert_eq!((0..tri.slots() as u32).filter(|&t| tri.is_ghost(t)).count(), 4);
+    check(&tri, 1.0 / 6.0);
+}
+
+#[test]
+fn random_points_are_delaunay_and_valid() {
+    let points = random(1, 1500);
+    let tri = Triangulation::new(&points).unwrap();
+    tri.verify().unwrap();
+    tri.verify_empty_spheres().unwrap();
+    // The tetrahedra fill the volume the hull faces (the ghosts' finite faces) enclose.
+    let vol = hull_volume(&tri);
+    assert!(vol > 0.9 && vol < 1.0, "{vol}");
+    gate(&tri, vol);
+}
+
+#[test]
+fn matches_the_definition_on_small_sets() {
+    let sets: Vec<Vec<[f64; 3]>> = vec![
+        random(7, 9),
+        random(8, 14),
+        grid(3, 3, 2),
+        grid(3, 3, 3),
+        // A cube's corners plus its centre: every face cospherical.
+        {
+            let mut g = grid(2, 2, 2);
+            g.push([0.5, 0.5, 0.5]);
+            g
+        },
+        // Coplanar on two planes.
+        vec![[0.0, 0.0, 0.0], [2.0, 0.0, 0.0], [0.0, 2.0, 0.0], [2.0, 2.0, 0.0], [1.0, 1.0, 0.0], [1.0, 1.0, 1.0], [0.0, 0.0, 1.0], [2.0, 2.0, 1.0]],
+    ];
+    for points in sets {
+        let tri = Triangulation::new(&points).unwrap();
+        tri.verify().unwrap();
+        assert_eq!(tet_set(&tri), brute_force(&points), "{points:?}");
+    }
+}
+
+#[test]
+fn the_result_does_not_depend_on_the_input_order() {
+    for points in [random(2, 400), grid(5, 4, 6)] {
+        let reference = tet_set(&Triangulation::new(&points).unwrap());
+        let mut r = Rng(99);
+        for _ in 0..3 {
+            let mut shuffled = points.clone();
+            r.shuffle(&mut shuffled);
+            let tri = Triangulation::new(&shuffled).unwrap();
+            tri.verify().unwrap();
+            assert_eq!(tet_set(&tri), reference);
+        }
+    }
+}
+
+#[test]
+fn deterministic_bit_for_bit() {
+    let points = random(3, 3000);
+    let a = Triangulation::new(&points).unwrap();
+    let b = Triangulation::new(&points).unwrap();
+    assert_eq!(a.tets, b.tets);
+    assert_eq!(a.adj, b.adj);
+}
+
+#[test]
+fn grids_are_cospherical_everywhere() {
+    let tri = Triangulation::new(&grid(7, 6, 5)).unwrap();
+    check(&tri, 6.0 * 5.0 * 4.0);
+    assert_eq!(hull_volume(&tri), 120.0);
+    // A rotated-in-place grid with non-integer spacing.
+    let skew: Vec<[f64; 3]> = grid(6, 6, 6).into_iter().map(|[x, y, z]| [0.1 * x, 0.3 * y + 1e3, 0.7 * z - 5.0]).collect();
+    let tri = Triangulation::new(&skew).unwrap();
+    tri.verify().unwrap();
+    tri.verify_empty_spheres().unwrap();
+}
+
+#[test]
+fn cospherical_and_cocircular_points() {
+    // Every integer point on the sphere x² + y² + z² = 9 (all cospherical).
+    let mut sphere = Vec::new();
+    for x in -3i32..=3 {
+        for y in -3i32..=3 {
+            for z in -3i32..=3 {
+                if x * x + y * y + z * z == 9 {
+                    sphere.push([f64::from(x), f64::from(y), f64::from(z)]);
+                }
+            }
+        }
+    }
+    assert_eq!(sphere.len(), 30);
+    let tri = Triangulation::new(&sphere).unwrap();
+    tri.verify().unwrap();
+    tri.verify_empty_spheres().unwrap();
+    assert_eq!(tet_set(&tri), brute_force(&sphere));
+
+    // Two regular 12-gons on parallel planes (a CAD cylinder's rings).
+    let mut cylinder = Vec::new();
+    for z in [0.0, 4.0] {
+        for i in 0..12 {
+            let a = f64::from(i) * std::f64::consts::TAU / 12.0;
+            cylinder.push([3.0 * a.cos(), 3.0 * a.sin(), z]);
+        }
+    }
+    let tri = Triangulation::new(&cylinder).unwrap();
+    tri.verify().unwrap();
+    tri.verify_empty_spheres().unwrap();
+    assert_eq!(tet_set(&tri), brute_force(&cylinder));
+}
+
+#[test]
+fn lower_dimensional_inputs_have_no_tetrahedra() {
+    assert_eq!(Triangulation::new(&[]).unwrap().dimension(), 0);
+    let one = Triangulation::new(&[[1.0, 2.0, 3.0]]).unwrap();
+    assert_eq!((one.dimension(), one.finite_tet_count()), (0, 0));
+    let line: Vec<[f64; 3]> = (0..50).map(|i| [f64::from(i), 2.0 * f64::from(i), -f64::from(i)]).collect();
+    let tri = Triangulation::new(&line).unwrap();
+    assert_eq!((tri.dimension(), tri.finite_tet_count()), (1, 0));
+    tri.verify().unwrap();
+
+    let plane: Vec<[f64; 3]> = grid(10, 10, 1).into_iter().map(|[x, y, _]| [x, y, x + 2.0 * y]).collect();
+    let mut tri = Triangulation::new(&plane).unwrap();
+    assert_eq!((tri.dimension(), tri.finite_tet_count()), (2, 0));
+    assert_eq!(tri.locate([0.0, 0.0, 0.0]).unwrap(), None);
+    // Duplicates of the planar points are still found.
+    assert_eq!(tri.insert([3.0, 4.0, 11.0]).unwrap(), Insertion::Duplicate(34));
+    // One point off the plane makes it 3D: every tetrahedron has that apex.
+    let apex = match tri.insert([4.5, 4.5, 40.0]).unwrap() {
+        Insertion::Inserted(v) => v,
+        other => panic!("{other:?}"),
+    };
+    assert_eq!(tri.dimension(), 3);
+    assert!(tri.finite_tets().all(|(_, vs)| vs.contains(&apex)));
+    tri.verify().unwrap();
+    tri.verify_empty_spheres().unwrap();
+}
+
+#[test]
+fn duplicates_merge_into_the_lowest_index() {
+    let base = grid(4, 4, 4);
+    let mut points = base.clone();
+    points.extend(base.iter().rev());
+    points.push([-0.0, 0.0, -0.0]);
+    let tri = Triangulation::new(&points).unwrap();
+    tri.verify().unwrap();
+    for i in 0..base.len() {
+        assert_eq!(tri.duplicate_of(i as u32), None);
+        assert_eq!(tri.duplicate_of((2 * base.len() - 1 - i) as u32), Some(i as u32));
+    }
+    assert_eq!(tri.duplicate_of(points.len() as u32 - 1), Some(0));
+    assert_eq!(tet_set(&tri), tet_set(&Triangulation::new(&base).unwrap()));
+}
+
+#[test]
+fn steiner_points_insert_into_a_built_triangulation() {
+    let mut tri = Triangulation::new(&grid(4, 4, 4)).unwrap();
+    let mut r = Rng(5);
+    for _ in 0..200 {
+        let p = [4.0 * r.unit() - 0.5, 4.0 * r.unit() - 0.5, 4.0 * r.unit() - 0.5];
+        assert!(matches!(tri.insert(p).unwrap(), Insertion::Inserted(_)));
+    }
+    // On a vertex, on an edge, on a face, on the hull's plane outside it.
+    assert_eq!(tri.insert([1.0, 2.0, 3.0]).unwrap(), Insertion::Duplicate(27));
+    assert_eq!(tri.duplicate_of(tri.points().len() as u32 - 1), Some(27));
+    for p in [[1.5, 2.0, 3.0], [1.5, 1.5, 0.0], [5.0, 5.0, 0.0], [-1.0, 1.0, 3.0]] {
+        assert!(matches!(tri.insert(p).unwrap(), Insertion::Inserted(_)), "{p:?}");
+    }
+    tri.verify().unwrap();
+    tri.verify_empty_spheres().unwrap();
+}
+
+#[test]
+fn locate_finds_the_containing_tetrahedron() {
+    let tri = Triangulation::new(&grid(3, 3, 3)).unwrap();
+    match tri.locate([0.5, 0.5, 0.5]).unwrap() {
+        Some(Location::Inside(t)) => assert!(!tri.is_ghost(t)),
+        other => panic!("{other:?}"),
+    }
+    assert!(matches!(tri.locate([9.0, 0.5, 0.5]).unwrap(), Some(Location::Outside(t)) if tri.is_ghost(t)));
+    assert_eq!(tri.locate([2.0, 2.0, 2.0]).unwrap(), Some(Location::Vertex(26)));
+    assert!(tri.locate([f64::NAN, 0.0, 0.0]).is_err());
+}
+
+#[test]
+fn hostile_coordinates_are_rejected() {
+    for bad in [[f64::NAN, 0.0, 0.0], [f64::INFINITY, 0.0, 0.0], [1e300, 0.0, 0.0], [1e-300, 0.0, 0.0]] {
+        assert!(matches!(Triangulation::new(&[[0.0; 3], bad]), Err(Error::Invalid(_))));
+    }
+    let mut tri = Triangulation::new(&grid(2, 2, 2)).unwrap();
+    assert!(tri.insert([0.0, f64::NEG_INFINITY, 0.0]).is_err());
+    tri.verify().unwrap();
+    // Extreme but admissible magnitudes still triangulate exactly.
+    let big: Vec<[f64; 3]> = random(4, 200).into_iter().map(|[x, y, z]| [x * 1e40, y * 1e-30 + 1.0, z]).collect();
+    Triangulation::new(&big).unwrap().verify().unwrap();
+}
+
+#[test]
+fn incident_tets_and_neighbors_are_consistent() {
+    let tri = Triangulation::new(&random(6, 300)).unwrap();
+    for v in 0..300 {
+        let t = tri.incident_tet(v).unwrap();
+        assert!(tri.tet(t).unwrap().contains(&v));
+    }
+    for (t, _) in tri.finite_tets() {
+        for n in tri.neighbors(t).unwrap() {
+            assert!(tri.neighbors(n).unwrap().contains(&t));
+        }
+    }
+}

crates/mesh/src/delaunay/verify.rs

@@ -0,0 +1,213 @@
+//! Independent checks of a triangulation: they re-derive every property from the arrays with
+//! the exact predicates; nothing here relies on how the insertion built them.
+
+use super::{DEAD, GHOST, NONE, Triangulation, volume_sign};
+use crate::predicates::{Point3, Sign, insphere, insphere_sos};
+use crate::quality::shape::{det, sub};
+use crate::{Error, Result};
+
+fn failed<T>(what: String) -> Result<T> {
+    Err(Error::Failed(format!("Delaunay verification: {what}")))
+}
+
+/// `det[b − a, c − a, d − a] / 6` in floating point (a measure for the volume balance only).
+pub(super) fn volume(p: [&Point3; 4]) -> f64 {
+    let [a, b, c, d] = p.map(Point3::coords);
+    det(sub(b, a), sub(c, a), sub(d, a)) / 6.0
+}
+
+/// Whether `e` is outside the circumsphere of the positive-volume tetrahedron `[a, b, c, d]`:
+/// never inside exactly, and outside with the perturbation.
+fn outside_sphere([a, b, c, d]: [&Point3; 4], e: &Point3) -> bool {
+    // `insphere` wants Shewchuk-positive orientation: swap the first two vertices.
+    match insphere(b, a, c, d, e) {
+        Sign::Positive => false,
+        Sign::Negative => true,
+        Sign::Zero => insphere_sos(b, a, c, d, e) == Sign::Negative,
+    }
+}
+
+/// The face of `vs` opposite slot `k`, sorted.
+pub(super) fn face(vs: [u32; 4], k: usize) -> [u32; 3] {
+    let mut f = [NONE; 3];
+    for (o, (_, &v)) in f.iter_mut().zip(vs.iter().enumerate().filter(|&(i, _)| i != k)) {
+        *o = v;
+    }
+    f.sort_unstable();
+    f
+}
+
+impl Triangulation {
+    fn point_of(&self, v: u32) -> Result<&Point3> {
+        match self.points.get(v as usize) {
+            Some(p) => Ok(p),
+            None => failed(format!("vertex {v} does not exist")),
+        }
+    }
+
+    fn points_of(&self, [a, b, c, d]: [u32; 4]) -> Result<[&Point3; 4]> {
+        Ok([self.point_of(a)?, self.point_of(b)?, self.point_of(c)?, self.point_of(d)?])
+    }
+
+    fn points_with<'a>(&'a self, vs: [u32; 4], k: usize, p: &'a Point3) -> Result<[&'a Point3; 4]> {
+        let mut out = [p; 4];
+        for (i, o) in out.iter_mut().enumerate() {
+            if i != k {
+                *o = self.point_of(vs.get(i).copied().unwrap_or(NONE))?;
+            }
+        }
+        Ok(out)
+    }
+
+    /// The vertex of `o` across its shared face with `t`, and the slot of that face in `o`.
+    fn across(&self, t: u32, o: u32) -> Result<([u32; 4], usize)> {
+        let Some(ws) = self.tet(o) else { return failed(format!("tetrahedron {t} has a deleted neighbour {o}")) };
+        let adj = self.adj.get(o as usize).copied().unwrap_or([NONE; 4]);
+        match (adj.iter().position(|&n| n == t), adj.iter().rposition(|&n| n == t)) {
+            (Some(k), Some(l)) if k == l => Ok((ws, k)),
+            _ => failed(format!("tetrahedra {t} and {o} are not mutual neighbours exactly once")),
+        }
+    }
+
+    /// Check the triangulation: adjacency is mutual and matches shared faces, every finite
+    /// tetrahedron has positive volume, the hull is convex, every interior face is locally
+    /// Delaunay (the opposite vertex is not inside the circumsphere, exactly, and is outside
+    /// with the perturbation), and every vertex is in a tetrahedron. As a floating-point sanity
+    /// check against a hull covered twice, the finite volume must match the volume the hull faces
+    /// enclose to 1e-9. By Delaunay's lemma, local Delaunay faces in a
+    /// triangulation of the hull make every circumsphere empty. `O(tetrahedra)`.
+    pub fn verify(&self) -> Result<()> {
+        let n = self.points.len();
+        let mut seen = vec![false; n];
+        let (mut finite_volume, mut hull_volume, mut scale) = (0.0f64, 0.0f64, 0.0f64);
+        let reference = self.tets.iter().flatten().copied().find(|&v| v < DEAD);
+        for (t, &vs) in self.tets.iter().enumerate() {
+            if vs[0] == DEAD {
+                continue;
+            }
+            let t = t as u32;
+            let adj = self.adj.get(t as usize).copied().unwrap_or([NONE; 4]);
+            for (i, &v) in vs.iter().enumerate() {
+                if vs.iter().skip(i + 1).any(|&w| w == v) {
+                    return failed(format!("tetrahedron {t} repeats vertex {v}"));
+                }
+                if v == GHOST {
+                    continue;
+                }
+                if self.duplicate_of(v).is_some() {
+                    return failed(format!("duplicate point {v} is a vertex"));
+                }
+                match seen.get_mut(v as usize) {
+                    Some(s) => *s = true,
+                    None => return failed(format!("tetrahedron {t} has vertex {v}, which does not exist")),
+                }
+            }
+            // Per face: the neighbour's vertices and the slot of the shared face in it.
+            let mut across = [([NONE; 4], 0); 4];
+            for (j, (&o, slot)) in adj.iter().zip(&mut across).enumerate() {
+                if o == t {
+                    return failed(format!("tetrahedron {t} is its own neighbour"));
+                }
+                let (ws, k) = self.across(t, o)?;
+                if face(vs, j) != face(ws, k) {
+                    return failed(format!("tetrahedra {t} and {o} are neighbours without sharing a face"));
+                }
+                *slot = (ws, k);
+            }
+            let apex = |j: usize| across.get(j).and_then(|&(ws, k)| ws.get(k).copied()).unwrap_or(NONE);
+            match vs.iter().position(|&v| v == GHOST) {
+                None => {
+                    let p = self.points_of(vs)?;
+                    if volume_sign(p) != Sign::Positive {
+                        return failed(format!("tetrahedron {t} is flat or inverted"));
+                    }
+                    let vol = volume(p);
+                    finite_volume += vol;
+                    scale += vol.abs();
+                    // The test is symmetric in the two tetrahedra: check each face once.
+                    for (j, &o) in adj.iter().enumerate() {
+                        let u = apex(j);
+                        if u != GHOST && t < o && !outside_sphere(p, self.point_of(u)?) {
+                            return failed(format!("vertex {u} is inside the circumsphere of tetrahedron {t}"));
+                        }
+                    }
+                }
+                Some(g) => {
+                    let fs = across.get(g).map(|&(ws, _)| ws).unwrap_or([NONE; 4]);
+                    if fs.contains(&GHOST) {
+                        return failed(format!("hull face of {t} borders another ghost"));
+                    }
+                    let inner = self.point_of(apex(g))?;
+                    if volume_sign(self.points_with(vs, g, inner)?) != Sign::Negative {
+                        return failed(format!("ghost tetrahedron {t} faces into the hull"));
+                    }
+                    for j in (0..4).filter(|&j| j != g) {
+                        let apex = self.point_of(apex(j))?;
+                        match volume_sign(self.points_with(vs, g, apex)?) {
+                            Sign::Positive => return failed(format!("the hull is not convex at ghost tetrahedron {t}")),
+                            Sign::Negative => {}
+                            Sign::Zero => {
+                                // Coplanar hull faces: the perturbed decision of the finite
+                                // tetrahedron under this face.
+                                if !outside_sphere(self.points_of(fs)?, apex) {
+                                    return failed(format!("coplanar hull faces at ghost tetrahedron {t} are not perturbed-Delaunay"));
+                                }
+                            }
+                        }
+                    }
+                    if let Some(r) = reference {
+                        hull_volume -= volume(self.points_with(vs, g, self.point_of(r)?)?);
+                    }
+                }
+            }
+        }
+        if self.dimension() == 3 {
+            for (v, s) in seen.iter().enumerate() {
+                let v = v as u32;
+                let triangulated = self.duplicate_of(v).is_none();
+                if *s != triangulated {
+                    return failed(format!("vertex {v} is {}in a tetrahedron", if *s { "" } else { "not " }));
+                }
+                if triangulated && !self.incident_tet(v).and_then(|t| self.tet(t)).is_some_and(|vs| vs.contains(&v)) {
+                    return failed(format!("vertex {v} has no valid incident tetrahedron"));
+                }
+            }
+        } else if self.tets.iter().any(|vs| vs[0] != DEAD) {
+            return failed("a lower-dimensional triangulation has tetrahedra".into());
+        }
+        if (finite_volume - hull_volume).abs() > 1e-9 * scale {
+            return failed(format!("the tetrahedra fill {finite_volume}, the hull encloses {hull_volume}"));
+        }
+        Ok(())
+    }
+
+    /// Brute force: no triangulated vertex is inside any circumsphere (exactly) or the perturbed
+    /// circumsphere, and none is beyond a hull face. `O(points × tetrahedra)`; for tests.
+    pub fn verify_empty_spheres(&self) -> Result<()> {
+        let vertices: Vec<u32> = (0..self.points.len() as u32).filter(|&v| self.duplicate_of(v).is_none()).collect();
+        for (t, &vs) in self.tets.iter().enumerate() {
+            if vs[0] == DEAD {
+                continue;
+            }
+            for &u in &vertices {
+                if vs.contains(&u) {
+                    continue;
+                }
+                let e = self.point_of(u)?;
+                match vs.iter().position(|&v| v == GHOST) {
+                    None => {
+                        if !outside_sphere(self.points_of(vs)?, e) {
+                            return failed(format!("vertex {u} is inside the circumsphere of tetrahedron {t}"));
+                        }
+                    }
+                    Some(g) => {
+                        if volume_sign(self.points_with(vs, g, e)?) == Sign::Positive {
+                            return failed(format!("vertex {u} is beyond the hull face of ghost tetrahedron {t}"));
+                        }
+                    }
+                }
+            }
+        }
+        Ok(())
+    }
+}

crates/mesh/src/lib.rs

@@ -7,6 +7,7 @@
 #![deny(clippy::unwrap_used, clippy::expect_used, clippy::panic, clippy::unimplemented, clippy::todo, clippy::unreachable)]
 #![forbid(unsafe_code)]
 
+pub mod delaunay;
 pub mod quality;
 
 #[derive(Debug, Clone, PartialEq, thiserror::Error)]

crates/mesh/src/predicates/mod.rs

@@ -104,6 +104,13 @@ impl std::ops::Mul for Sign {
     }
 }
 
+/// The sign of the volume `det[b − a, c − a, d − a]`: positive when `d` lies on the side of
+/// face `a, b, c`'s right-hand normal (the tet4 node order). It is [`orient3d`] negated.
+#[inline]
+pub fn tet_volume_sign(a: &Point3, b: &Point3, c: &Point3, d: &Point3) -> Sign {
+    -orient3d(a, b, c, d)
+}
+
 /// Whether `x` may be a point coordinate (see the module docs).
 #[inline]
 pub fn is_admissible(x: f64) -> bool {

crates/mesh/src/quality/mod.rs

@@ -21,14 +21,14 @@
 mod bernstein;
 pub mod jacobian;
 mod report;
-mod shape;
+pub(crate) mod shape;
 
 use modelcraft_model::{ElementKind, check_coordinates};
 
 pub use jacobian::{JacobianCheck, Validity};
 pub use report::{GateOptions, GateReport, MetricStats, QualityReport, enclosed_volume, quality_report, validity_gate, validity_gate_elements};
 
-use crate::predicates::{Point3, Sign, orient3d};
+use crate::predicates::{Point3, Sign, tet_volume_sign};
 use crate::{Error, Result};
 
 /// A quality metric. Angles are in degrees.
@@ -164,13 +164,7 @@ fn checked(kind: ElementKind, xyz: &[[f64; 3]], p: &[[f64; 3]]) -> Result<(Optio
 /// right-hand normal), or `None` when a coordinate is outside the predicates' domain.
 fn exact_tet_sign(xyz: &[[f64; 3]]) -> Option<Sign> {
     let p = |i: usize| xyz.get(i).and_then(|c| Point3::try_from(*c).ok());
-    // orient3d(a, b, c, d) is positive when d lies below the plane of a, b, c seen
+    Some(tet_volume_sign(&p(0)?, &p(1)?, &p(2)?, &p(3)?))
-    // counterclockwise, i.e. opposite to our convention.
-    Some(match orient3d(&p(0)?, &p(1)?, &p(2)?, &p(3)?) {
-        Sign::Positive => Sign::Negative,
-        Sign::Negative => Sign::Positive,
-        Sign::Zero => Sign::Zero,
-    })
 }
 
 /// The certified Jacobian check and exact volume of a solid element (tet4, tet10, wedge6, hex8,

crates/mesh/tests/delaunay_bench.rs

@@ -0,0 +1,50 @@
+//! Throughput of the Delaunay tetrahedralization on 1M random points and on a 100³ grid (every
+//! cell cospherical). Ignored by default; run with
+//! `cargo test --release -p modelcraft-mesh --test delaunay_bench -- --ignored --nocapture`.
+
+use std::time::Instant;
+
+use modelcraft_mesh::delaunay::Triangulation;
+
+struct Rng(u64);
+
+impl Rng {
+    fn next_u64(&mut self) -> u64 {
+        self.0 = self.0.wrapping_add(0x9E37_79B9_7F4A_7C15);
+        let mut z = self.0;
+        z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
+        z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
+        z ^ (z >> 31)
+    }
+    fn unit(&mut self) -> f64 {
+        (self.next_u64() >> 11) as f64 / (1u64 << 53) as f64
+    }
+}
+
+fn run(name: &str, points: &[[f64; 3]]) {
+    let start = Instant::now();
+    let tri = Triangulation::new(points).unwrap();
+    let build = start.elapsed().as_secs_f64();
+    let start = Instant::now();
+    tri.verify().unwrap();
+    let verify = start.elapsed().as_secs_f64();
+    let tets = tri.finite_tet_count();
+    println!(
+        "{name}: {} points -> {tets} tets in {build:.2} s ({:.0} points/s, {:.0} tets/s); verify {verify:.2} s",
+        points.len(),
+        points.len() as f64 / build,
+        tets as f64 / build
+    );
+}
+
+#[test]
+#[ignore]
+fn delaunay_throughput() {
+    let mut r = Rng(1);
+    let random: Vec<[f64; 3]> = (0..1_000_000).map(|_| [r.unit(), r.unit(), r.unit()]).collect();
+    run("random 1M", &random);
+
+    let n = 100;
+    let grid: Vec<[f64; 3]> = (0..n * n * n).map(|i| [(i % n) as f64, ((i / n) % n) as f64, (i / (n * n)) as f64]).collect();
+    run("grid 100^3", &grid);
+}

docs/meshing/literature.md

@@ -128,7 +128,14 @@ Delaunay methods with Voronoi/ODT improvement, then (2) all-hex and hex-dominant
   - Visibility or straight walking for point location: O. Devillers, S. Pion, M. Teillaud, *Walking in a
     triangulation*, IJFCS 13(2):181–199, 2002, DOI [10.1142/S0129054102001047](https://doi.org/10.1142/S0129054102001047).
   - A compact tet-adjacency array.
+  - Termination of the visibility walk in a Delaunay triangulation: H. Edelsbrunner, *An acyclicity theorem for cell
+    complexes in d dimensions*, Combinatorica 10(3):251–260, 1990, DOI
+    [10.1007/BF02122779](https://doi.org/10.1007/BF02122779).
+  - Hilbert index: J. Skilling, *Programming the Hilbert curve*, AIP Conf. Proc. 707:381–387, 2004, DOI
+    [10.1063/1.1751381](https://doi.org/10.1063/1.1751381).
 - **Code.** For a permissive reference, see geogram `delaunay_3d.cpp` (BSD).
+- **Implemented** in `modelcraft-mesh::delaunay` (ghost vertex for the hull, Devillers–Teillaud perturbation for
+  cospherical points, exact verifier).
 
 **A1.2 Marot, Pellerin, Remacle 2019, "One machine, one minute, three billion tetrahedra".** **Verdict: reference (later).**
 - **Citation.** IJNME 117(9):967–990, 2019. DOI [10.1002/nme.5987](https://doi.org/10.1002/nme.5987).