diff --git a/ROADMAP.md b/ROADMAP.md index a4a1002..8eed8e9 100644 --- a/ROADMAP.md +++ b/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 | diff --git a/crates/engine/src/cmd/mesh.rs b/crates/engine/src/cmd/mesh.rs index da60020..e439388 100644 --- a/crates/engine/src/cmd/mesh.rs +++ b/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 { + 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 { 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(); diff --git a/crates/engine/src/cmd/mod.rs b/crates/engine/src/cmd/mod.rs index 06aa4f0..d1936c6 100644 --- a/crates/engine/src/cmd/mod.rs +++ b/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 { /// 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]")); - 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) { - *o = v.as_f64().ok_or_else(bad)?; + p.get(key).and_then(xyz_value).ok_or_else(|| EngineError::Params(format!("`{key}` must be [x, y, z]"))) +} + +/// A JSON `[x, y, z]`. +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> { + 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. diff --git a/crates/engine/src/lib.rs b/crates/engine/src/lib.rs index 20386d4..1122928 100644 --- a/crates/engine/src/lib.rs +++ b/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 { diff --git a/crates/mesh/src/delaunay/mod.rs b/crates/mesh/src/delaunay/mod.rs new file mode 100644 index 0000000..189c909 --- /dev/null +++ b/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, + /// `NONE`, or the vertex this point duplicates. + duplicate_of: Vec, + tets: Vec<[u32; 4]>, + /// Neighbour `i` is across the face opposite vertex `i`. + adj: Vec<[u32; 4]>, + free: Vec, + /// A live tetrahedron incident to each vertex (`NONE` until it is in one). + vertex_tet: Vec, + /// While the points span less than 3D: affinely independent vertices, and the others. + basis: Vec, + pending: Vec, + /// Cavity marks: `epoch` = in the cavity, `epoch + 1` = tested and outside. + mark: Vec, + 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, + stack: Vec, + /// `(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::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 { + 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::>()?; + 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 = (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 = (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 { + 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 { + 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 { + 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 + '_ { + 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> { + 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 { + 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 { + 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 { + 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 { + 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 { + 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 { + 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 { + 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(()) + } +} diff --git a/crates/mesh/src/delaunay/order.rs b/crates/mesh/src/delaunay/order.rs new file mode 100644 index 0000000..4e3fdf0 --- /dev/null +++ b/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 { + 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 = (0..10_000).map(|i| Point3::new(f64::from(i % 97), f64::from(i % 89), f64::from(i % 83)).unwrap()).collect(); + let vertices: Vec = (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"); + } +} diff --git a/crates/mesh/src/delaunay/tests.rs b/crates/mesh/src/delaunay/tests.rs new file mode 100644 index 0000000..60aec50 --- /dev/null +++ b/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(&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 = 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 = 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![ + 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)); + } + } +} diff --git a/crates/mesh/src/delaunay/verify.rs b/crates/mesh/src/delaunay/verify.rs new file mode 100644 index 0000000..723a795 --- /dev/null +++ b/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(what: String) -> Result { + 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 = (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(()) + } +} diff --git a/crates/mesh/src/lib.rs b/crates/mesh/src/lib.rs index 0185565..b1c0e78 100644 --- a/crates/mesh/src/lib.rs +++ b/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)] diff --git a/crates/mesh/src/predicates/mod.rs b/crates/mesh/src/predicates/mod.rs index d3d4f56..22848c4 100644 --- a/crates/mesh/src/predicates/mod.rs +++ b/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 { diff --git a/crates/mesh/src/quality/mod.rs b/crates/mesh/src/quality/mod.rs index e69751f..35a2431 100644 --- a/crates/mesh/src/quality/mod.rs +++ b/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 { 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 - // 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, - }) + Some(tet_volume_sign(&p(0)?, &p(1)?, &p(2)?, &p(3)?)) } /// The certified Jacobian check and exact volume of a solid element (tet4, tet10, wedge6, hex8, diff --git a/crates/mesh/tests/delaunay_bench.rs b/crates/mesh/tests/delaunay_bench.rs new file mode 100644 index 0000000..5e1c791 --- /dev/null +++ b/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); +} diff --git a/docs/meshing/literature.md b/docs/meshing/literature.md index 89821a7..8144ed2 100644 --- a/docs/meshing/literature.md +++ b/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).