mesh: the tet mesher takes its input in a canonical order (#115)

SolidCraft lists a solid's faces in a different order from run to run
(the rounded box, for one), which reordered the triangle soup and the
face tags with every coordinate the same; welding keeps the first copy
it meets and every later stage visits ids in order, so the tet mesh
changed between runs. CanonicalSoup sorts the triangles by their
coordinates, numbers vertices by first use and face tags by first
appearance, and maps the tags back; mesh.tet.generate and the corpus
use it.

Tests: the washer's tet10 mesh has a fixed fingerprint (identical on
arm64 and x86-64 builds); a reordered, renumbered input gives the same
mesh; heavy: rounded-box and filleted-bracket mesh bit for bit the same
on repeated runs and with 1 and 8 threads.

The bracket oracle cases stay unfrozen: SolidCraft's tessellation
itself differs in the last bits between arm64 and x86-64. The coarse
case moves to 203.6 MPa with the canonical order; its band is 5%.

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

crates/mesh/src/tet/canonical.rs

@@ -0,0 +1,152 @@
+//! The canonical form of a boundary triangle soup: the tetrahedral mesh of a part depends on its
+//! geometry only, never on the order in which the CAD kernel lists its faces, triangles and
+//! vertices.
+//!
+//! A kernel may enumerate the faces of one solid in a different order from run to run (its
+//! internal maps are hash maps), which reorders the triangle soup and renumbers the face tags
+//! while every coordinate stays the same. Welding keeps the first copy of a point it meets, and
+//! every later stage visits vertices and faces in id order, so a reordered soup gives a different
+//! mesh. [`CanonicalSoup`] sorts the triangles by their coordinates (rotated to start at their
+//! smallest corner, which keeps their orientation), numbers the vertices by first use in that
+//! order and the face tags by first appearance, and maps the tags back afterwards.
+
+use std::collections::BTreeMap;
+
+use crate::{Error, Result};
+
+use super::QualityMesh;
+use super::complex::NONE;
+use super::geometry::AnalyticSurface;
+
+/// A triangle soup in canonical order, and the original face tag of each canonical one.
+#[derive(Debug, Clone, PartialEq)]
+pub struct CanonicalSoup {
+    pub positions: Vec<[f64; 3]>,
+    pub triangles: Vec<[u32; 3]>,
+    pub tri_face: Vec<u32>,
+    /// `faces[c]`: the original tag of canonical tag `c`.
+    pub faces: Vec<u32>,
+}
+
+/// A coordinate as an integer whose order is the numeric order of the coordinate (−0 before +0;
+/// the coordinates of a valid soup are finite).
+fn ordered(x: f64) -> u64 {
+    let b = x.to_bits();
+    if b >> 63 == 1 { !b } else { b | (1 << 63) }
+}
+
+fn point_key(p: [f64; 3]) -> [u64; 3] {
+    p.map(ordered)
+}
+
+impl CanonicalSoup {
+    /// The canonical form of `triangles` over `positions`, each triangle on face `tri_face[i]`
+    /// (all on face 0 when `tri_face` is empty).
+    pub fn new(positions: &[[f64; 3]], triangles: &[[u32; 3]], tri_face: &[u32]) -> Result<CanonicalSoup> {
+        if !tri_face.is_empty() && tri_face.len() != triangles.len() {
+            return Err(Error::Invalid(format!("{} face tags for {} triangles", tri_face.len(), triangles.len())));
+        }
+        let point = |v: u32| positions.get(v as usize).copied().ok_or_else(|| Error::Invalid(format!("triangle vertex {v} does not exist")));
+        // Each triangle rotated to start at its smallest corner, keyed by its corners' coordinates.
+        let mut order: Vec<([[u64; 3]; 3], [u32; 3], u32)> = Vec::with_capacity(triangles.len());
+        for (i, t) in triangles.iter().enumerate() {
+            let keys = [point_key(point(t[0])?), point_key(point(t[1])?), point_key(point(t[2])?)];
+            let first = (0..3).min_by(|&a, &b| keys[a].cmp(&keys[b])).unwrap_or(0);
+            let r = |k: usize| (first + k) % 3;
+            let tag = tri_face.get(i).copied().unwrap_or(0);
+            order.push(([keys[r(0)], keys[r(1)], keys[r(2)]], [t[r(0)], t[r(1)], t[r(2)]], tag));
+        }
+        // Equal keys are the same triangle twice; the original tag breaks the tie only between
+        // copies on different faces, which no valid soup has.
+        order.sort_unstable_by(|a, b| a.0.cmp(&b.0).then(a.2.cmp(&b.2)));
+        let mut id = vec![NONE; positions.len()];
+        let mut face_id: BTreeMap<u32, u32> = BTreeMap::new();
+        let mut out = CanonicalSoup {
+            positions: Vec::with_capacity(positions.len()),
+            triangles: Vec::with_capacity(triangles.len()),
+            tri_face: Vec::with_capacity(triangles.len()),
+            faces: Vec::new(),
+        };
+        for (_, t, tag) in order {
+            let mut nt = [0u32; 3];
+            for (slot, v) in nt.iter_mut().zip(t) {
+                let s = id.get_mut(v as usize).ok_or_else(|| Error::Invalid(format!("triangle vertex {v} does not exist")))?;
+                if *s == NONE {
+                    *s = out.positions.len() as u32;
+                    out.positions.push(point(v)?);
+                }
+                *slot = *s;
+            }
+            let c = *face_id.entry(tag).or_insert_with(|| {
+                out.faces.push(tag);
+                out.faces.len() as u32 - 1
+            });
+            out.triangles.push(nt);
+            out.tri_face.push(c);
+        }
+        Ok(out)
+    }
+
+    /// The canonical tag of an original face tag.
+    pub fn canonical_tag(&self, original: u32) -> Option<u32> {
+        self.faces.iter().position(|&f| f == original).map(|c| c as u32)
+    }
+
+    /// The exact surfaces by original tag, reordered by canonical tag.
+    pub fn analytic(&self, by_original: &[Option<AnalyticSurface>]) -> Vec<Option<AnalyticSurface>> {
+        self.faces.iter().map(|&f| by_original.get(f as usize).copied().flatten()).collect()
+    }
+
+    /// Give a mesh of the canonical soup its original face tags back.
+    pub fn restore_tags(&self, mesh: &mut QualityMesh) {
+        for b in &mut mesh.boundary {
+            if let Some(&f) = self.faces.get(b.face as usize) {
+                b.face = f;
+            }
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    /// A tetrahedron's surface listed in two orders, with its vertices and faces renumbered,
+    /// has one canonical form, and the tags map back.
+    #[test]
+    fn reordered_soups_have_one_canonical_form() {
+        let p = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]];
+        let t = [[0, 2, 1], [0, 1, 3], [0, 3, 2], [1, 2, 3]];
+        let a = CanonicalSoup::new(&p, &t, &[10, 11, 12, 13]).unwrap();
+        // Vertices reversed, triangles reversed and rotated, tags permuted.
+        let q: Vec<[f64; 3]> = p.iter().rev().copied().collect();
+        let map = |v: u32| 3 - v;
+        let u: Vec<[u32; 3]> = t.iter().rev().map(|f| [map(f[1]), map(f[2]), map(f[0])]).collect();
+        let b = CanonicalSoup::new(&q, &u, &[3, 2, 1, 0]).unwrap();
+        assert_eq!((&a.positions, &a.triangles, &a.tri_face), (&b.positions, &b.triangles, &b.tri_face));
+        for c in 0..4u32 {
+            // Triangle k of `t` is on face 10 + k in a and, listed as u[3 − k], on face k in b.
+            assert_eq!(a.faces[c as usize] - 10, b.faces[c as usize]);
+            assert_eq!(a.canonical_tag(a.faces[c as usize]), Some(c));
+        }
+        // Orientation is kept: the enclosed volume has the same sign.
+        let vol = |s: &CanonicalSoup| -> f64 {
+            s.triangles
+                .iter()
+                .map(|f| {
+                    let [x, y, z] = f.map(|v| s.positions[v as usize]);
+                    x[0] * (y[1] * z[2] - y[2] * z[1]) - x[1] * (y[0] * z[2] - y[2] * z[0]) + x[2] * (y[0] * z[1] - y[1] * z[0])
+                })
+                .sum()
+        };
+        assert!(vol(&a) != 0.0 && vol(&a) == vol(&b));
+    }
+
+    #[test]
+    fn hostile_soups_are_refused() {
+        let p = [[0.0, 0.0, 0.0]];
+        assert!(CanonicalSoup::new(&p, &[[0, 0, 7]], &[]).is_err());
+        assert!(CanonicalSoup::new(&p, &[[0, 0, 0]], &[1, 2]).is_err());
+        assert_eq!(CanonicalSoup::new(&p, &[], &[]).unwrap().triangles.len(), 0);
+    }
+}

crates/mesh/src/tet/mod.rs

@@ -33,6 +33,7 @@
 // clearer than iterator chains (as in the fem and solver crates).
 #![allow(clippy::needless_range_loop)]
 
+pub mod canonical;
 pub mod complex;
 pub mod geometry;
 pub mod improve;

crates/mesh/src/tet/tests.rs

@@ -349,3 +349,108 @@ mod pipeline {
         }
     }
 }
+
+/// Determinism (#115): the quality mesh is a function of the geometry, the same bit for bit on
+/// every platform, run and listing order of the input.
+mod determinism {
+    use crate::cdt::{CdtOptions, tetrahedralize};
+    use crate::surface::{BoundarySurface, SurfaceOptions};
+    use crate::tet::canonical::CanonicalSoup;
+    use crate::tet::{Order, QualityMesh, QualityOptions, quality_mesh};
+
+    /// cos and sin of k · 15° from square roots only (correctly rounded everywhere), so the
+    /// test input itself is the same on every platform.
+    fn unit_24(k: usize) -> [f64; 2] {
+        let (s2, s3, s6) = (2f64.sqrt(), 3f64.sqrt(), 6f64.sqrt());
+        let c = [1.0, (s6 + s2) / 4.0, s3 / 2.0, s2 / 2.0, 0.5, (s6 - s2) / 4.0, 0.0];
+        // Quadrant q, step j: cos(90° q + 15° j).
+        let (q, j) = (k / 6 % 4, k % 6);
+        let (x, y) = (c[j], c[6 - j]);
+        match q {
+            0 => [x, y],
+            1 => [-y, x],
+            2 => [-x, -y],
+            _ => [y, -x],
+        }
+    }
+
+    /// A washer (inner radius 6, outer 15, height 4) as a 24-sided prism ring: faces bottom 0,
+    /// top 1, inside 2, outside 3; triangles oriented outward.
+    pub(crate) fn washer() -> (Vec<[f64; 3]>, Vec<[u32; 3]>, Vec<u32>) {
+        const N: usize = 24;
+        let mut p = Vec::new();
+        for (r, z) in [(6.0, 0.0), (15.0, 0.0), (6.0, 4.0), (15.0, 4.0)] {
+            for k in 0..N {
+                let [c, s] = unit_24(k);
+                p.push([r * c, r * s, z]);
+            }
+        }
+        let id = |ring: usize, k: usize| (ring * N + k % N) as u32;
+        let (mut t, mut f) = (Vec::new(), Vec::new());
+        for k in 0..N {
+            // (ring a, ring b): quad a_k b_k b_k+1 a_k+1, outward when listed this way.
+            for (a, b, tag) in [(1, 0, 0), (2, 3, 1), (0, 2, 2), (3, 1, 3)] {
+                let q = [id(a, k), id(b, k), id(b, k + 1), id(a, k + 1)];
+                t.push([q[0], q[1], q[2]]);
+                t.push([q[0], q[2], q[3]]);
+                f.extend([tag, tag]);
+            }
+        }
+        (p, t, f)
+    }
+
+    fn mesh(p: &[[f64; 3]], t: &[[u32; 3]], f: &[u32]) -> QualityMesh {
+        let soup = CanonicalSoup::new(p, t, f).unwrap();
+        let s = BoundarySurface::new(&soup.positions, &soup.triangles, &soup.tri_face, &SurfaceOptions::default()).unwrap();
+        assert!(s.volume() > 0.0);
+        let cdt = tetrahedralize(&s, &CdtOptions::default()).unwrap();
+        let mut m =
+            quality_mesh(&s, &cdt, Vec::new(), &QualityOptions { size: Some(3.0), order: Order::Quadratic, ..QualityOptions::default() }).unwrap();
+        soup.restore_tags(&mut m);
+        m
+    }
+
+    /// FNV-1a over every coordinate's bits, the elements and the tagged boundary faces.
+    pub(crate) fn fingerprint(m: &QualityMesh) -> u64 {
+        let mut h = 0xcbf2_9ce4_8422_2325u64;
+        let mut eat = |x: u64| {
+            for b in x.to_le_bytes() {
+                h = (h ^ u64::from(b)).wrapping_mul(0x0100_0000_01b3);
+            }
+        };
+        m.points.iter().flatten().for_each(|c| eat(c.to_bits()));
+        m.elements.iter().flatten().for_each(|&v| eat(u64::from(v)));
+        m.boundary.iter().for_each(|b| eat(u64::from(b.tet) << 40 | u64::from(b.opposite) << 32 | u64::from(b.face)));
+        h
+    }
+
+    /// The washer's tet10 mesh has one fingerprint on every platform (CI checks it on Linux
+    /// x86-64, macOS arm64 and Windows): only correctly rounded arithmetic and the crate's own
+    /// `atan2` and `cbrt` make its decisions, and no hash map is iterated into it. A change that
+    /// moves the tet mesh on purpose updates the constant.
+    #[test]
+    fn the_quality_mesh_is_the_same_on_every_platform() {
+        let (p, t, f) = washer();
+        let m = mesh(&p, &t, &f);
+        assert!(m.elements.len() > 100, "{:?}", m.report);
+        assert_eq!(fingerprint(&m), 6_285_915_972_108_509_003, "{} tets, {} nodes", m.elements.len(), m.points.len());
+    }
+
+    /// Listing the triangles, their corners and the vertices in another order, and numbering
+    /// the faces differently, gives the same mesh with the faces' own tags.
+    #[test]
+    fn the_quality_mesh_does_not_depend_on_the_input_order() {
+        let (p, t, f) = washer();
+        let a = mesh(&p, &t, &f);
+        let n = p.len() as u32;
+        let q: Vec<[f64; 3]> = p.iter().rev().copied().collect();
+        let u: Vec<[u32; 3]> = t.iter().rev().map(|x| [n - 1 - x[1], n - 1 - x[2], n - 1 - x[0]]).collect();
+        let g: Vec<u32> = f.iter().rev().map(|&x| 10 + 3 - x).collect();
+        let b = mesh(&q, &u, &g);
+        assert_eq!((&a.points, &a.elements), (&b.points, &b.elements));
+        assert_eq!(a.boundary.len(), b.boundary.len());
+        for (x, y) in a.boundary.iter().zip(&b.boundary) {
+            assert_eq!((x.tet, x.opposite, 10 + 3 - x.face), (y.tet, y.opposite, y.face));
+        }
+    }
+}