diff --git a/ROADMAP.md b/ROADMAP.md index 4781f18..d674d0c 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -20,7 +20,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 and self-intersection checks, feature edges (`modelcraft-mesh::surface`, `mesh.surface.check`, done) | in progress | -| M2 Tet meshing and size field | Size field (global, per face and edge seeds, curvature and proximity; the `modelcraft_mesh::size::SizeField` trait and edge-division helpers are in); 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 | +| M2 Tet meshing and size field | Size field (global, per face and edge seeds, curvature and proximity; the `modelcraft_mesh::size::SizeField` trait and edge-division helpers are in); Bowyer–Watson Delaunay (exact, symbolic perturbation, verified; `modelcraft-mesh::delaunay`, done); robust CDT boundary recovery (Diazzi et al. 2023: implicit LNC Steiner points on segments, gift-wrapping facet recovery, inside/outside labelling, B-rep face tags on boundary tet faces, few-ulp rounding repair checked by the validity gate; `modelcraft-mesh::cdt` and `mesh.tet.generate`, done: all 20 valid corpus parts); 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; `model.check` naming every missing or bad item; `model.*` and `mesh.tag` commands (model definition: done); writers for Abaqus input (done: `deck.write`), Nastran SOL 101/103 (done: `deck.write`) and LS-DYNA implicit (done); the `.modelcraft` file (text definition, binary mesh, geometry references reopened on load, journal; atomic saves, autosave and recovery; undo/redo on the journal: done) | in progress | | M4 Own solver and first results | Element library per docs/elements.md (1D and 2D first); Linear static and modal solver (tet4/10, wedge6, hex8, hex20; sparse direct, Lanczos), verified against analytic solutions (solver: done — `modelcraft-fem` kernels, `modelcraft-solver` with nested dissection, sparse Cholesky and block Lanczos, `solve.static`/`solve.frequency`/`solve.probe`, patch, cantilever, Lamé, Kirsch and beam-frequency tests; reference points with kinematic and distributing couplings and moments by master–slave constraint elimination, verified by rigid motion, clamping, Saint-Venant tension, bending and tube torsion); deformed shape, contours, mode shapes, probes, reaction sums | started | | M5 Structured hex | 2D quad meshing of faces (`modelcraft-mesh::quad`: transfinite mapped meshes, free all-quad meshes of polygons with holes on prescribed boundary nodes, exact planar quad validity, shared `SizeField` and even edge divisions; done for planar faces); mapped hex of six-sided regions; sweep meshing; conforming interfaces and tet–hex transition; automatic decomposition into sweepable and mappable regions, proposed as SolidCraft splits | started | diff --git a/crates/engine/src/catalogue.rs b/crates/engine/src/catalogue.rs index ebafb67..c9d087b 100644 --- a/crates/engine/src/catalogue.rs +++ b/crates/engine/src/catalogue.rs @@ -159,7 +159,13 @@ pub static CATALOGUE: &[Capability] = &[ soon("mesh.size.global", "Global Size", Mesh, "Size", "M2", "One element size for the whole part"), soon("mesh.size.seed", "Seed Edges", Mesh, "Size", "M2", "Element counts, sizes and bias per edge and face"), soon("mesh.size.adapt", "Curvature", Mesh, "Size", "M2", "Refine by curvature and proximity"), - soon("mesh.tet.generate", "Tet Mesh", Mesh, "Tet", "M2", "Delaunay tetrahedral mesh of a body, boundary recovered and refined"), + ready( + "mesh.tet.generate", + "Tet Mesh", + Mesh, + "Tet", + "Constrained Delaunay tet mesh of a body: every boundary triangle kept, faces tagged (refinement comes next)", + ), soon("mesh.order.quadratic", "Quadratic", Mesh, "Tet", "M2", "Second-order elements with mid-side nodes on the geometry"), ready("mesh.delaunay.points", "Delaunay Points", Mesh, "Tet", "Delaunay tetrahedralization of a point list (debugging aid)"), soon("mesh.quad.faces", "Quad Faces", Mesh, "Structured Hex", "M5", "All-quad meshes of faces on shared edge nodes"), diff --git a/crates/engine/src/cmd/mesh.rs b/crates/engine/src/cmd/mesh.rs index f41b064..8235a60 100644 --- a/crates/engine/src/cmd/mesh.rs +++ b/crates/engine/src/cmd/mesh.rs @@ -1,10 +1,11 @@ use std::borrow::Cow; use modelcraft_geometry::SurfaceMesh; +use modelcraft_mesh::cdt::{CdtOptions, tetrahedralize}; use modelcraft_mesh::delaunay::Triangulation; use modelcraft_mesh::quality::{self, QualityReport}; use modelcraft_mesh::surface::{BoundarySurface, Defect, FeatureKind, SurfaceOptions}; -use modelcraft_model::{ElementKind, MAX_ELEMENTS}; +use modelcraft_model::{ElementKind, GeometryRef, MAX_ELEMENTS, MAX_LABEL}; use serde_json::{Value, json}; use super::{label, labels, listed, optional_number, points, xyz}; @@ -78,6 +79,102 @@ pub const SURFACE_CHECK: CommandSpec = CommandSpec { run: surface_check, }; +pub const TET_GENERATE: CommandSpec = CommandSpec { + id: "mesh.tet.generate", + label: "Generate Tet Mesh", + params: r#"{"body": "Body1", "tolerance": 0.05} or {"surface": 0}, the boundary options of mesh.surface.check ("weld", "weld_distance", "feature_angle"), optional "max_steiner": 50000 (0..=10000000), "max_tets": 20000000 (1..=50000000) — the constrained Delaunay tetrahedralization of the welded boundary: every boundary triangle is a union of tet faces, Steiner points only on boundary edges, no quality refinement yet. Adds tet4 elements and their nodes (labels after the largest existing ones), tags the elements with the body and each boundary element face with its geometry face, and fails, changing nothing, if the surface is invalid or the mesh rounded to doubles does not pass the validity gate"#, + kind: CommandKind::Model, + run: tet_generate, +}; + +/// The most Steiner points and tetrahedra `mesh.tet.generate` may be asked for. +const TET_MAX_STEINER: u64 = 10_000_000; +const TET_MAX_TETS: u64 = 50_000_000; + +/// An optional whole-number parameter in `range`. +fn optional_count(p: &Value, key: &str, default: u64, range: std::ops::RangeInclusive) -> Result { + let v = match p.get(key) { + None => default, + Some(v) => v + .as_u64() + .filter(|n| range.contains(n)) + .ok_or_else(|| EngineError::Params(format!("`{key}` must be a whole number in {}..={}", range.start(), range.end())))?, + }; + usize::try_from(v).map_err(|_| EngineError::Params(format!("`{key}` is too large"))) +} + +/// The boundary options of `mesh.surface.check` and `mesh.tet.generate`. +fn surface_options(p: &Value, mesh: &SurfaceMesh) -> Result { + let defaults = SurfaceOptions::default(); + Ok(SurfaceOptions { + weld_tolerance: optional_number(p, "weld", defaults.weld_tolerance)?, + weld_distance: optional_number(p, "weld_distance", mesh.weld_distance)?, + feature_angle: optional_number(p, "feature_angle", defaults.feature_angle)?, + }) +} + +fn tet_generate(s: &mut Session, p: &Value) -> Result { + let defaults = CdtOptions::default(); + let cdt_options = CdtOptions { + max_steiner: optional_count(p, "max_steiner", defaults.max_steiner as u64, 0..=TET_MAX_STEINER)?, + max_tets: optional_count(p, "max_tets", defaults.max_tets as u64, 1..=TET_MAX_TETS)?, + }; + let (surface, body, face_names) = { + let Input { mesh, body, faces } = surface_input(s, p)?; + let options = surface_options(p, &mesh)?; + (BoundarySurface::new(&mesh.positions, &mesh.triangles, &mesh.tri_face, &options)?, body, faces) + }; + let cdt = tetrahedralize(&surface, &cdt_options)?; + let gate = cdt.gate(surface.volume())?.into_result()?; + + let next = |last: Option| last.map_or(Some(1), |l| l.checked_add(1)); + let first_node = next(s.model.nodes().keys().next_back().copied()); + let first_element = next(s.model.elements().keys().next_back().copied()); + let fits = |first: Option, count: usize| first.filter(|&f| u64::from(f) + count as u64 - 1 <= u64::from(MAX_LABEL)); + let (Some(first_node), Some(first_element)) = (fits(first_node, cdt.points.len()), fits(first_element, cdt.tets.len())) else { + return Err(EngineError::Params(format!("the mesh needs labels beyond {MAX_LABEL}"))); + }; + for (i, xyz) in cdt.points.iter().enumerate() { + s.model.set_node(first_node + i as u32, *xyz)?; + } + let body_tag = s.model.geometry_tag(&GeometryRef::body(body.clone())?)?; + for (i, t) in cdt.tets.iter().enumerate() { + let label = first_element + i as u32; + s.model.set_element(label, ElementKind::Tet4, t.iter().map(|&v| first_node + v).collect())?; + s.model.tag_element(label, body_tag)?; + } + let mut face_tags: Vec> = vec![None; face_names.len()]; + for b in &cdt.boundary { + let face = surface.tri_face().get(b.triangle as usize).copied().unwrap_or(u32::MAX) as usize; + let tag = match face_tags.get(face).copied().flatten() { + Some(t) => t, + None => { + let name = + face_names.get(face).ok_or_else(|| EngineError::Params(format!("boundary triangle {} has no geometry face", b.triangle)))?; + let t = s.model.geometry_tag(&GeometryRef::face(body.clone(), name.clone())?)?; + if let Some(slot) = face_tags.get_mut(face) { + *slot = Some(t); + } + t + } + }; + let number = b.tet4_face().ok_or_else(|| EngineError::Params(format!("boundary face {} has no tet4 face number", b.opposite)))?; + s.model.tag_element_face(first_element + b.tet, number, tag)?; + } + Ok(json!({ + "body": body, + "nodes": cdt.points.len(), + "elements": cdt.tets.len(), + "first_node": first_node, + "first_element": first_element, + "steiner": cdt.steiner.len(), + "cavities": cdt.stats.cavities, + "boundary_faces": cdt.boundary.len(), + "volume": gate.volume, + "boundary_volume": surface.volume(), + })) +} + fn node_add(s: &mut Session, p: &Value) -> Result { let l = label(p, "label")?; s.model.set_node(l, xyz(p, "xyz")?)?; @@ -150,17 +247,29 @@ fn quality_report(s: &mut Session, p: &Value) -> Result { })) } +/// A boundary to mesh: the surface, with the body name and the geometry name of each face index +/// for tagging. +struct Input<'a> { + mesh: Cow<'a, SurfaceMesh>, + body: String, + faces: Vec, +} + /// The surface named by `body` (tessellated at `tolerance`) or `surface` (an imported index). -fn surface_input<'a>(s: &'a Session, p: &Value) -> Result> { +fn surface_input<'a>(s: &'a Session, p: &Value) -> Result> { match (p.get("body"), p.get("surface")) { - (Some(_), None) => Ok(Cow::Owned(super::geometry::body_mesh(s, p)?.1)), + (Some(_), None) => { + let (shape, mesh) = super::geometry::body_mesh(s, p)?; + Ok(Input { mesh: Cow::Owned(mesh), body: shape.name.clone(), faces: shape.faces().iter().map(|f| f.name.clone()).collect() }) + } (None, Some(i)) => { let n = s.surfaces.len(); - i.as_u64() + let t = i + .as_u64() .and_then(|i| usize::try_from(i).ok()) .and_then(|i| s.surfaces.get(i)) - .map(|t| Cow::Borrowed(&t.mesh)) - .ok_or_else(|| EngineError::Params(format!("`surface` must be the index of an imported surface (there are {n})"))) + .ok_or_else(|| EngineError::Params(format!("`surface` must be the index of an imported surface (there are {n})")))?; + Ok(Input { mesh: Cow::Borrowed(&t.mesh), body: t.name.clone(), faces: (0..t.faces).map(|f| f.to_string()).collect() }) } _ => Err(EngineError::Params("give either `body` (with `tolerance`) or `surface`".into())), } @@ -171,13 +280,8 @@ fn defect_json(d: &Defect) -> Value { } fn surface_check(s: &mut Session, p: &Value) -> Result { - let mesh = surface_input(s, p)?; - let defaults = SurfaceOptions::default(); - let options = SurfaceOptions { - weld_tolerance: optional_number(p, "weld", defaults.weld_tolerance)?, - weld_distance: optional_number(p, "weld_distance", mesh.weld_distance)?, - feature_angle: optional_number(p, "feature_angle", defaults.feature_angle)?, - }; + let mesh = surface_input(s, p)?.mesh; + let options = surface_options(p, &mesh)?; let surface = match BoundarySurface::new(&mesh.positions, &mesh.triangles, &mesh.tri_face, &options) { Ok(surface) => surface, Err(modelcraft_mesh::Error::Surface(d)) => { @@ -350,4 +454,77 @@ mod tests { assert!(s.execute("mesh.surface.check", &bad).is_err(), "{bad}"); } } + + #[test] + fn tet_generate_meshes_a_body_with_tags_and_appends() { + use std::sync::Arc; + + use modelcraft_geometry::open_design; + use modelcraft_geometry::sample::{Plate, plate_design}; + + let plate = Plate::default(); + let mut s = Session { design: Some(Arc::new(open_design(&plate_design(&plate).unwrap()).unwrap())), ..Session::default() }; + let r = s.execute("mesh.tet.generate", &json!({"body": "Body1", "tolerance": 0.5})).unwrap(); + let elements = r["elements"].as_u64().unwrap() as usize; + assert!(elements > 0, "{r}"); + assert_eq!(s.model.elements().len(), elements); + assert!(s.model.elements().values().all(|e| e.kind == modelcraft_model::ElementKind::Tet4)); + let (volume, boundary) = (r["volume"].as_f64().unwrap(), r["boundary_volume"].as_f64().unwrap()); + assert!((volume - boundary).abs() <= 1e-9 * boundary, "{r}"); + // Every element is on the body; boundary faces carry the 8 face names. + let c = s.model.classification(); + assert_eq!(c.elements().len(), elements); + let faces: std::collections::BTreeSet = c.elements().values().flat_map(|t| t.faces.values().copied()).collect(); + assert_eq!(faces.len(), 8); + let tagged: usize = c.elements().values().map(|t| t.faces.len()).sum(); + assert_eq!(tagged as u64, r["boundary_faces"].as_u64().unwrap()); + // The certified quality report agrees. + let q = s.execute("mesh.quality", &json!({})).unwrap(); + assert_eq!(q["valid"].as_u64().unwrap() as usize, elements); + // A second run appends after the existing labels. + let again = s.execute("mesh.tet.generate", &json!({"body": "Body1", "tolerance": 0.5})).unwrap(); + assert_eq!(again["first_element"].as_u64().unwrap() as usize, elements + 1); + assert_eq!(s.model.elements().len(), 2 * elements); + assert_eq!(s.journal.len(), 2); + + let journal = s.journal.len(); + let model = s.model.clone(); + for bad in [ + json!({"body": "Body1", "tolerance": 0.5, "max_steiner": -1}), + json!({"body": "Body1", "tolerance": 0.5, "max_steiner": 1e300}), + json!({"body": "Body1", "tolerance": 0.5, "max_tets": 0}), + json!({"body": "Body1", "tolerance": 0.5, "max_tets": 10}), + json!({"body": "Body1", "tolerance": 0.5, "max_tets": "many"}), + json!({"body": "Body1"}), + json!({"body": "Nope", "tolerance": 0.5}), + json!({"surface": 0}), + json!({}), + ] { + assert!(s.execute("mesh.tet.generate", &bad).is_err(), "{bad}"); + } + assert_eq!((s.journal.len(), &s.model), (journal, &model), "failures change nothing"); + } + + #[test] + fn tet_generate_on_an_imported_surface_and_an_open_one() { + use std::sync::Arc; + + use modelcraft_geometry::SurfaceMesh; + use modelcraft_geometry::tagged_mesh::TaggedSurface; + + let p = vec![[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 = vec![[0, 2, 1], [0, 1, 3], [1, 2, 3], [0, 3, 2]]; + let tet = SurfaceMesh { positions: p.clone(), triangles: t.clone(), tri_face: vec![0, 1, 1, 1], weld_distance: 0.0 }; + let open = SurfaceMesh { positions: p, triangles: t[..3].to_vec(), tri_face: vec![0, 1, 1], weld_distance: 0.0 }; + let mut s = Session::default(); + for mesh in [tet, open] { + s.surfaces.push(Arc::new(TaggedSurface { name: "Tet".into(), mesh, faces: 2 })); + } + let r = s.execute("mesh.tet.generate", &json!({"surface": 0})).unwrap(); + assert_eq!((r["elements"].clone(), r["nodes"].clone(), r["boundary_faces"].clone()), (json!(1), json!(4), json!(4))); + let element = s.model.element(1).unwrap(); + assert_eq!(element.nodes.len(), 4); + assert!(s.execute("mesh.tet.generate", &json!({"surface": 1})).is_err(), "an open surface bounds no volume"); + assert_eq!(s.model.elements().len(), 1); + } } diff --git a/crates/engine/src/cmd/mod.rs b/crates/engine/src/cmd/mod.rs index 41c0e08..58fc1f7 100644 --- a/crates/engine/src/cmd/mod.rs +++ b/crates/engine/src/cmd/mod.rs @@ -29,6 +29,7 @@ pub static ALL: &[CommandSpec] = &[ mesh::QUALITY, mesh::DELAUNAY_POINTS, mesh::SURFACE_CHECK, + mesh::TET_GENERATE, tags::TAG, tags::TAG_CLEAR, model::UNITS, diff --git a/crates/mesh/src/cdt/complex.rs b/crates/mesh/src/cdt/complex.rs new file mode 100644 index 0000000..3c055bb --- /dev/null +++ b/crates/mesh/src/cdt/complex.rs @@ -0,0 +1,590 @@ +//! The working tetrahedralization of the boundary recovery: tetrahedra over explicit and implicit +//! ([`Lnc`](crate::predicates::Lnc)) vertices, with the operations the recovery needs: point +//! location, Bowyer–Watson insertion of Steiner points, vertex stars, edge lookup, and the +//! replacement of a set of tetrahedra by another filling the same region. +//! +//! Layout: per tetrahedron four vertex ids in tet4 order (`det[b − a, c − a, d − a] > 0`, decided +//! exactly) and four neighbour ids, neighbour `i` across the face opposite vertex `i`; [`NONE`] on +//! the hull. The hull is the bounding box the recovery adds around the surface, so every Steiner +//! point is strictly inside it and no ghost tetrahedra are needed. + +use std::collections::HashMap; +use std::hash::{BuildHasherDefault, Hasher}; + +use crate::delaunay::Triangulation; +use crate::predicates::{GenericPoint, Sign, generic_insphere_sos, generic_lex, generic_orient3d}; +use crate::{Error, Result}; + +/// A multiplicative hasher for small integer keys (vertex triples). The maps it serves are only +/// looked up, never iterated, so its order cannot leak into the output. +#[derive(Default)] +pub(super) struct IdHasher(u64); + +impl Hasher for IdHasher { + fn finish(&self) -> u64 { + self.0 + } + + fn write(&mut self, bytes: &[u8]) { + for &b in bytes { + self.write_u64(u64::from(b)); + } + } + + fn write_u32(&mut self, x: u32) { + self.write_u64(u64::from(x)); + } + + fn write_u64(&mut self, x: u64) { + self.0 = (self.0.rotate_left(5) ^ x).wrapping_mul(0x51_7c_c1_b7_27_22_0a_95); + } +} + +/// A hash map keyed by vertex ids. +pub(super) type IdMap = HashMap>; + +/// No tetrahedron / no vertex. +pub(super) const NONE: u32 = u32::MAX; +/// Marks a recycled tetrahedron slot. +const DEAD: u32 = u32::MAX - 1; + +pub(super) fn corrupt(what: &str) -> Error { + Error::Failed(format!("constrained tetrahedralization is inconsistent: {what}")) +} + +/// The face of `vs` opposite slot `k`, ordered so that vertex `k` is on its positive side +/// ([`Complex::volume`] of the face plus that vertex is positive). +pub(super) fn face_toward(vs: [u32; 4], k: usize) -> [u32; 3] { + let [a, b, c, d] = vs; + match k { + 0 => [b, d, c], + 1 => [a, c, d], + 2 => [a, d, b], + _ => [a, b, c], + } +} + +/// A triangle's vertices sorted: the key of an unoriented face. +pub(super) fn sorted3([a, b, c]: [u32; 3]) -> [u32; 3] { + let mut f = [a, b, c]; + f.sort_unstable(); + f +} + +/// The working tetrahedralization. +#[derive(Debug, Clone)] +pub(super) struct Complex { + points: Vec, + /// A box certainly containing each point (implicit points are not exactly representable). + lo: Vec<[f64; 3]>, + hi: Vec<[f64; 3]>, + /// Rounded coordinates (floating-point estimates only, never decisions). + approx: Vec<[f64; 3]>, + tets: Vec<[u32; 4]>, + adj: Vec<[u32; 4]>, + free: Vec, + vertex_tet: Vec, + /// Traversal marks: a tetrahedron is visited in the current traversal when its mark equals + /// `epoch`. + mark: Vec, + epoch: u32, + rng: u64, + /// Most tetrahedron slots allowed. + max_tets: usize, + live: usize, +} + +/// Where a point lies ([`Complex::locate`]). +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub(super) enum Located { + /// In the closed tetrahedron and not one of its vertices. + Inside(u32), + /// Equal to this vertex. + Vertex(u32), +} + +impl Complex { + /// The finite tetrahedra of a Delaunay triangulation of `points` (in the same vertex order), + /// whose hull faces get no neighbour. + pub(super) fn from_delaunay(tri: &Triangulation, points: Vec, max_tets: usize) -> Result { + let mut id = vec![NONE; tri.slots()]; + let mut tets = Vec::new(); + for (t, vs) in tri.finite_tets() { + if let Some(slot) = id.get_mut(t as usize) { + *slot = tets.len() as u32; + } + tets.push(vs); + } + if tets.len() > max_tets { + return Err(Error::Limit(format!("{} tetrahedra; the limit is {max_tets}", tets.len()))); + } + let mut adj = Vec::with_capacity(tets.len()); + for (t, _) in tri.finite_tets() { + let ns = tri.neighbors(t).ok_or_else(|| corrupt("Delaunay neighbours"))?; + adj.push(ns.map(|n| id.get(n as usize).copied().unwrap_or(NONE))); + } + let (lo, hi) = points.iter().map(GenericPoint::bounds).unzip(); + let approx = points.iter().map(GenericPoint::approx).collect(); + let mut cx = Complex { + vertex_tet: vec![NONE; points.len()], + points, + lo, + hi, + approx, + mark: vec![0; tets.len()], + live: tets.len(), + tets, + adj, + free: Vec::new(), + epoch: 0, + rng: 0x9E37_79B9_7F4A_7C15, + max_tets, + }; + for (t, vs) in cx.tets.iter().enumerate() { + for &v in vs { + if let Some(s) = cx.vertex_tet.get_mut(v as usize) { + *s = t as u32; + } + } + } + Ok(cx) + } + + // ---- read access ---- + + pub(super) fn points(&self) -> &[GenericPoint] { + &self.points + } + + pub(super) fn pt(&self, v: u32) -> Result<&GenericPoint> { + self.points.get(v as usize).ok_or_else(|| corrupt("vertex id")) + } + + /// Vertex `v`'s coordinates rounded to doubles: for estimates and messages only. + pub(super) fn approx(&self, v: u32) -> Result<[f64; 3]> { + self.approx.get(v as usize).copied().ok_or_else(|| corrupt("vertex id")) + } + + /// The box certainly containing vertex `v`. + pub(super) fn bounds(&self, v: u32) -> ([f64; 3], [f64; 3]) { + match (self.lo.get(v as usize), self.hi.get(v as usize)) { + (Some(&l), Some(&h)) => (l, h), + _ => ([f64::NEG_INFINITY; 3], [f64::INFINITY; 3]), + } + } + + /// The box certainly containing all of `vs`. + pub(super) fn bounds_of(&self, vs: &[u32]) -> ([f64; 3], [f64; 3]) { + let mut lo = [f64::INFINITY; 3]; + let mut hi = [f64::NEG_INFINITY; 3]; + for &v in vs { + let (l, h) = self.bounds(v); + for k in 0..3 { + lo[k] = lo[k].min(l[k]); + hi[k] = hi[k].max(h[k]); + } + } + (lo, hi) + } + + pub(super) fn slots(&self) -> usize { + self.tets.len() + } + + pub(super) fn live_count(&self) -> usize { + self.live + } + + /// The vertices of a live tetrahedron. + pub(super) fn tet(&self, t: u32) -> Result<[u32; 4]> { + self.tets.get(t as usize).copied().filter(|v| v[0] != DEAD).ok_or_else(|| corrupt("tetrahedron id")) + } + + pub(super) fn is_live(&self, t: u32) -> bool { + self.tets.get(t as usize).is_some_and(|v| v[0] != DEAD) + } + + pub(super) fn adj(&self, t: u32) -> Result<[u32; 4]> { + self.tet(t)?; + self.adj.get(t as usize).copied().ok_or_else(|| corrupt("tetrahedron id")) + } + + /// The live tetrahedra in slot order. + pub(super) fn live_tets(&self) -> impl Iterator + '_ { + self.tets.iter().enumerate().filter(|(_, v)| v[0] != DEAD).map(|(t, v)| (t as u32, *v)) + } + + /// The sign of `det[b − a, c − a, d − a]` (tet4 orientation), exact. + pub(super) fn volume(&self, [a, b, c, d]: [u32; 4]) -> Result { + Ok(-generic_orient3d(self.pt(a)?, self.pt(b)?, self.pt(c)?, self.pt(d)?)) + } + + /// Whether `x` is inside the perturbed circumsphere of the positive tetrahedron `vs`. + pub(super) fn in_sphere(&self, [a, b, c, d]: [u32; 4], x: u32) -> Result { + // `insphere` wants Shewchuk-positive orientation (tet4-negative): swap the first two. + Ok(generic_insphere_sos(self.pt(b)?, self.pt(a)?, self.pt(c)?, self.pt(d)?, self.pt(x)?) == Sign::Positive) + } + + /// Whether two vertices are the same point (exact). + fn same_point(&self, u: u32, v: u32) -> Result { + Ok(generic_lex(self.pt(u)?, self.pt(v)?).is_eq()) + } + + // ---- traversal marks ---- + + /// A new traversal: marks `epoch` (visited) and `epoch + 1` (rejected) are fresh. + fn next_epoch(&mut self) { + if self.epoch >= u32::MAX - 4 { + self.mark.fill(0); + self.epoch = 0; + } + self.epoch += 2; + } + + fn visited(&self, t: u32) -> bool { + self.mark.get(t as usize) == Some(&self.epoch) + } + + fn visit(&mut self, t: u32) { + if let Some(m) = self.mark.get_mut(t as usize) { + *m = self.epoch; + } + } + + fn rejected(&self, t: u32) -> bool { + self.mark.get(t as usize) == Some(&(self.epoch + 1)) + } + + fn reject(&mut self, t: u32) { + if let Some(m) = self.mark.get_mut(t as usize) { + *m = self.epoch + 1; + } + } + + // ---- stars and edges ---- + + /// A live tetrahedron incident to `v`. + fn incident(&self, v: u32) -> Result { + let t = self.vertex_tet.get(v as usize).copied().unwrap_or(NONE); + if self.tets.get(t as usize).is_some_and(|vs| vs[0] != DEAD && vs.contains(&v)) { + return Ok(t); + } + // The hint is maintained on every allocation, so this scan is a safety net only. + self.live_tets().find(|(_, vs)| vs.contains(&v)).map(|(t, _)| t).ok_or_else(|| corrupt("a vertex is in no tetrahedron")) + } + + /// The tetrahedra incident to `v` (its star), into `out`. + pub(super) fn star(&mut self, v: u32, out: &mut Vec) -> Result<()> { + out.clear(); + let start = self.incident(v)?; + self.next_epoch(); + self.visit(start); + out.push(start); + let mut i = 0; + while let Some(&t) = out.get(i) { + i += 1; + let vs = self.tet(t)?; + let ns = self.adj(t)?; + for (&u, &n) in vs.iter().zip(&ns) { + if u == v || n == NONE || self.visited(n) { + continue; + } + self.visit(n); + out.push(n); + } + if out.len() > self.live { + return Err(corrupt("a vertex star larger than the mesh")); + } + } + Ok(()) + } + + /// Whether `u v` is an edge. + pub(super) fn has_edge(&mut self, u: u32, v: u32, scratch: &mut Vec) -> Result { + self.star(u, scratch)?; + for &t in scratch.iter() { + if self.tet(t)?.contains(&v) { + return Ok(true); + } + } + Ok(false) + } + + /// The tetrahedra whose (conservative) box meets the box `lo..hi`, connected to the star of + /// `v` through such tetrahedra; at most `cap` of them. + pub(super) fn tets_near(&mut self, v: u32, lo: [f64; 3], hi: [f64; 3], cap: usize) -> Result> { + let mut out = Vec::new(); + self.star(v, &mut out)?; + self.next_epoch(); + for &t in &out { + self.visit(t); + } + let mut i = 0; + while let Some(&t) = out.get(i) { + i += 1; + if out.len() >= cap { + break; + } + for n in self.adj(t)? { + if n == NONE || self.visited(n) { + continue; + } + self.visit(n); + let (l, h) = self.bounds_of(&self.tet(n)?); + if (0..3).all(|k| l[k] <= hi[k] && h[k] >= lo[k]) { + out.push(n); + } + } + } + Ok(out) + } + + // ---- point location ---- + + /// The visibility walk towards `p` from `start` (Devillers, Pion, Teillaud 2002), capped at + /// the number of tetrahedra, then a scan. + pub(super) fn locate(&mut self, p: u32, start: u32) -> Result { + let mut t = if self.is_live(start) { start } else { self.live_tets().next().map(|(t, _)| t).ok_or_else(|| corrupt("empty mesh"))? }; + let mut prev = NONE; + for _ in 0..self.live.saturating_add(16) { + let vs = self.tet(t)?; + let ns = self.adj(t)?; + self.rng ^= self.rng << 13; + self.rng ^= self.rng >> 7; + self.rng ^= self.rng << 17; + let first = (self.rng >> 61) as usize; + let mut next = None; + for i in 0..4 { + let k = (first + i) % 4; + let n = ns.get(k).copied().unwrap_or(NONE); + if (n != prev || n == NONE) && self.beyond(vs, k, p)? { + next = Some(n); + break; + } + } + match next { + None => return self.inside(vs, t, p), + Some(NONE) => return Err(Error::Failed("a Steiner point lies outside the bounding box".into())), + Some(n) => { + prev = t; + t = n; + } + } + } + // The walk may cycle in a non-Delaunay mesh: scan. + let live: Vec<(u32, [u32; 4])> = self.live_tets().collect(); + for (t, vs) in live { + let mut inside = true; + for k in 0..4 { + if self.beyond(vs, k, p)? { + inside = false; + break; + } + } + if inside { + return self.inside(vs, t, p); + } + } + Err(Error::Failed("a Steiner point is in no tetrahedron".into())) + } + + /// Whether `p` is strictly beyond the face of `vs` opposite slot `k`. + fn beyond(&self, mut vs: [u32; 4], k: usize, p: u32) -> Result { + if let Some(s) = vs.get_mut(k) { + *s = p; + } + Ok(self.volume(vs)? == Sign::Negative) + } + + fn inside(&self, vs: [u32; 4], t: u32, p: u32) -> Result { + for v in vs { + if self.same_point(v, p)? { + return Ok(Located::Vertex(v)); + } + } + Ok(Located::Inside(t)) + } + + // ---- modification ---- + + /// Add a point (not yet in any tetrahedron); its vertex id. + pub(super) fn push_point(&mut self, p: GenericPoint) -> Result { + if self.points.len() >= (NONE - 2) as usize { + return Err(Error::Limit("too many vertices".into())); + } + let (l, h) = p.bounds(); + self.approx.push(p.approx()); + self.points.push(p); + self.lo.push(l); + self.hi.push(h); + self.vertex_tet.push(NONE); + Ok((self.points.len() - 1) as u32) + } + + /// Insert vertex `v` (from [`Complex::push_point`]) by Bowyer–Watson with the perturbed + /// in-sphere test: the mesh must be Delaunay for that predicate, which it is until the facet + /// recovery starts. The walk starts at a tetrahedron incident to `near`. Returns the vertices + /// of the deleted tetrahedra (the only ones whose edges can have changed). + pub(super) fn insert(&mut self, v: u32, near: u32) -> Result> { + let start = self.incident(near)?; + let seed = match self.locate(v, start)? { + Located::Vertex(u) => { + let at = self.approx(u)?; + return Err(Error::Failed(format!("a Steiner point coincides with vertex {u} at {at:?}"))); + } + Located::Inside(t) => t, + }; + // The cavity: every tetrahedron whose perturbed circumsphere contains v, grown from the + // one containing it (which always conflicts: v is in the closed tetrahedron, so strictly + // inside its ball). + self.next_epoch(); + let mut cavity = vec![seed]; + self.visit(seed); + let mut boundary: Vec<(u32, usize)> = Vec::new(); + let mut i = 0; + while let Some(&t) = cavity.get(i) { + i += 1; + for (j, n) in self.adj(t)?.into_iter().enumerate() { + if n == NONE || self.rejected(n) { + boundary.push((t, j)); + continue; + } + if self.visited(n) { + continue; + } + if self.in_sphere(self.tet(n)?, v)? { + self.visit(n); + cavity.push(n); + } else { + self.reject(n); + boundary.push((t, j)); + } + } + if cavity.len() > self.live { + return Err(corrupt("a cavity larger than the mesh")); + } + } + let mut new = Vec::with_capacity(boundary.len()); + for (t, j) in boundary { + let mut vs = self.tet(t)?; + if let Some(s) = vs.get_mut(j) { + *s = v; + } + if self.volume(vs)? != Sign::Positive { + return Err(corrupt("a Bowyer–Watson cavity is not star-shaped from the new point")); + } + new.push(vs); + } + let mut touched: Vec = Vec::with_capacity(4 * cavity.len()); + for &t in &cavity { + touched.extend(self.tet(t)?); + } + self.replace(&cavity, &new)?; + touched.sort_unstable(); + touched.dedup(); + Ok(touched) + } + + 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() >= self.max_tets { + return Err(Error::Limit(format!("more than {} tetrahedra", self.max_tets))); + } + 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; + } + } + self.live += 1; + 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); + self.live = self.live.saturating_sub(1); + Ok(()) + } + + fn set_adj(&mut self, t: u32, i: usize, n: u32) -> Result<()> { + *self.adj.get_mut(t as usize).and_then(|a| a.get_mut(i)).ok_or_else(|| corrupt("tetrahedron id"))? = n; + Ok(()) + } + + /// Replace the tetrahedra `old` by `new` (positive tet4 tuples) filling the same region: each + /// face of a new tetrahedron is either shared with exactly one other new tetrahedron or is a + /// face of the region's boundary, where it takes over the outside neighbour. Anything else is + /// reported, never patched. The new slot ids. + pub(super) fn replace(&mut self, old: &[u32], new: &[[u32; 4]]) -> Result> { + // Read the outside links before any slot is recycled. + // Boundary face → (outside neighbour, its slot back, covered yet). + let mut external: IdMap<[u32; 3], (u32, usize, bool)> = IdMap::with_capacity_and_hasher(old.len() * 2, Default::default()); + self.next_epoch(); + for &t in old { + self.visit(t); + } + for &t in old { + let vs = self.tet(t)?; + for (j, n) in self.adj(t)?.into_iter().enumerate() { + if n != NONE && self.visited(n) { + continue; + } + let back = if n == NONE { 0 } else { self.adj(n)?.iter().position(|&m| m == t).ok_or_else(|| corrupt("one-sided adjacency"))? }; + if external.insert(sorted3(face_toward(vs, j)), (n, back, false)).is_some() { + return Err(corrupt("a region boundary face is listed twice")); + } + } + } + for &t in old { + self.kill(t)?; + } + let mut ids = Vec::with_capacity(new.len()); + let mut open: IdMap<[u32; 3], (u32, usize)> = IdMap::with_capacity_and_hasher(new.len() * 2, Default::default()); + let mut matched_external = 0usize; + for &vs in new { + let t = self.alloc(vs)?; + ids.push(t); + for k in 0..4 { + let key = sorted3(face_toward(vs, k)); + if let Some(entry) = external.get_mut(&key) { + let (n, back, covered) = *entry; + if covered { + return Err(corrupt("a region boundary face is covered twice")); + } + entry.2 = true; + self.set_adj(t, k, n)?; + if n != NONE { + self.set_adj(n, back, t)?; + } + matched_external += 1; + } else if let Some((u, l)) = open.remove(&key) { + self.set_adj(t, k, u)?; + self.set_adj(u, l, t)?; + } else { + open.insert(key, (t, k)); + } + } + } + if !open.is_empty() || matched_external != external.len() { + return Err(corrupt(&format!( + "a retriangulated region does not close: {} unmatched inner faces, {} of {} boundary faces covered", + open.len(), + matched_external, + external.len() + ))); + } + Ok(ids) + } +} diff --git a/crates/mesh/src/cdt/facets.rs b/crates/mesh/src/cdt/facets.rs new file mode 100644 index 0000000..c564302 --- /dev/null +++ b/crates/mesh/src/cdt/facets.rs @@ -0,0 +1,414 @@ +//! Facet recovery (Diazzi et al. 2023, §3.4 and §4.4). Every input triangle is a facet; once all +//! segments are edges, a facet is *missing* when an edge of the tetrahedralization pierces its +//! interior. The tetrahedra with such edges form face-connected *cavities* (each meets the +//! facet's plane inside the facet only, because the facet's edges are mesh edges). Each cavity is +//! deleted and split by the facet's plane into two half-cavities, which are refilled one after +//! the other by the modified gift-wrapping of §4.4 (after Shewchuk 2000, *Sweep algorithms for +//! constructing higher-dimensional constrained Delaunay triangulations*): repeatedly take an +//! open triangle `σ` of the half-cavity's front and join it to the vertex `w` of the half-cavity +//! such that +//! - (i) `w` is on the open side of `σ` (positive volume), +//! - (ii) the tetrahedron meets every cavity, front and constraint triangle properly +//! ([`super::meet`]), and contains no other cavity vertex, +//! - (iii) among those, the perturbed circumsphere of `σ w` contains no other candidate: the +//! candidate that is first in the pencil of spheres through `σ`, with the same symbolic +//! perturbation as the Delaunay tetrahedralization, so both half-cavities induce the same +//! triangulation on the facet. +//! +//! The first half-cavity leaves its triangles on the facet's plane open; the second closes them. +//! Constraint triangles already recovered that lie inside a cavity stay constraints: both of +//! their sides enter the front. +//! +//! Si's cavity expansion with local Delaunay tetrahedralizations (§3.4) is a speed-up the paper +//! uses before falling back to gift-wrapping; gift-wrapping alone is the always-correct path +//! (§4.4: it builds the constrained Delaunay tetrahedralization of any half-cavity that admits +//! one), so it is the one implemented here. + +use std::collections::{HashMap, HashSet}; + +use crate::predicates::{GenericPoint, Sign, generic_orient3d}; +use crate::{Error, Result}; + +use super::Recovery; +use super::complex::{NONE, corrupt, face_toward, sorted3}; +use super::meet::{in_closed_tet, meet_properly}; + +/// Most rounds of finding and refilling the cavities of one facet (one suffices in exact +/// arithmetic). +const MAX_ROUNDS_PER_FACET: usize = 16; + +/// A box: `(lo, hi)`. +type Bounds = ([f64; 3], [f64; 3]); + +/// The open front of a half-cavity: faces, their boxes, and whether each is still open. +type Front<'a> = (&'a [[u32; 3]], &'a [Bounds], &'a [bool]); + +/// An oriented triangle rotated so its smallest id comes first: the key of an oriented face. +fn canon([a, b, c]: [u32; 3]) -> [u32; 3] { + if a <= b && a <= c { + [a, b, c] + } else if b <= a && b <= c { + [b, c, a] + } else { + [c, a, b] + } +} + +fn boxes_meet((l1, h1): Bounds, (l2, h2): Bounds) -> bool { + (0..3).all(|k| l1[k] <= h2[k] && l2[k] <= h1[k]) +} + +/// The side of the facet's plane each vertex is on, cached for one facet. +struct Plane { + corners: [GenericPoint; 3], + side: HashMap, +} + +impl Recovery { + fn plane(&self, f: usize) -> Result { + let tri = *self.triangles.get(f).ok_or_else(|| corrupt("facet id"))?; + let corner = |v: u32| self.cx.pt(v).cloned(); + Ok(Plane { corners: [corner(tri[0])?, corner(tri[1])?, corner(tri[2])?], side: HashMap::new() }) + } + + fn facet_side(&self, plane: &mut Plane, v: u32) -> Result { + if let Some(&s) = plane.side.get(&v) { + return Ok(s); + } + let [a, b, c] = &plane.corners; + // tet4 volume sign of (a, b, c, v): positive on the side the facet's normal points to. + let s = -generic_orient3d(a, b, c, self.cx.pt(v)?); + plane.side.insert(v, s); + Ok(s) + } + + /// Whether edge `p q` crosses the facet's open interior. + fn pierces(&self, plane: &mut Plane, p: u32, q: u32) -> Result { + if self.facet_side(plane, p)? * self.facet_side(plane, q)? != Sign::Negative { + return Ok(false); + } + let [a, b, c] = &plane.corners; + let (pp, pq) = (self.cx.pt(p)?, self.cx.pt(q)?); + let s0 = generic_orient3d(pp, pq, a, b); + if s0.is_zero() { + return Ok(false); + } + Ok(generic_orient3d(pp, pq, b, c) == s0 && generic_orient3d(pp, pq, c, a) == s0) + } + + fn crosses(&self, plane: &mut Plane, t: u32) -> Result { + let vs = self.cx.tet(t)?; + for (i, j) in [(0, 1), (0, 2), (0, 3), (1, 2), (1, 3), (2, 3)] { + if let (Some(&p), Some(&q)) = (vs.get(i), vs.get(j)) + && self.pierces(plane, p, q)? + { + return Ok(true); + } + } + Ok(false) + } + + /// The vertices on facet `f`'s boundary: its corners and the Steiner points on its edges. + pub(super) fn facet_ring(&self, f: usize) -> Result> { + let tri = *self.triangles.get(f).ok_or_else(|| corrupt("facet id"))?; + let mut ring = Vec::new(); + for k in 0..3 { + let (a, b) = (tri[k], tri[(k + 1) % 3]); + let s = *self.edge_segment.get(&(a.min(b), a.max(b))).ok_or_else(|| corrupt("facet edge without a segment"))?; + let chain = self.segment(s)?.chain_from(a); + ring.extend(chain.iter().take(chain.len().saturating_sub(1))); + } + Ok(ring) + } + + /// The face-connected sets of tetrahedra crossing facet `f`'s interior. Each contains a + /// tetrahedron incident to a vertex of the facet's boundary `ring`. + fn crossing_cavities(&mut self, plane: &mut Plane, ring: &[u32]) -> Result>> { + let mut cavities = Vec::new(); + let mut seen: HashSet = HashSet::new(); + let mut star = Vec::new(); + for &v in ring { + self.cx.star(v, &mut star)?; + for &t in &star { + if seen.contains(&t) || !self.crosses(plane, t)? { + continue; + } + seen.insert(t); + let mut cavity = vec![t]; + let mut i = 0; + while let Some(&c) = cavity.get(i) { + i += 1; + for n in self.cx.adj(c)? { + if n != NONE && !seen.contains(&n) && self.crosses(plane, n)? { + seen.insert(n); + cavity.push(n); + } + } + } + cavities.push(cavity); + } + } + Ok(cavities) + } + + /// Whether facet `f` is a union of mesh triangles: no tetrahedron crosses it. + pub(super) fn facet_present(&mut self, f: usize) -> Result { + let ring = self.facet_ring(f)?; + let mut plane = self.plane(f)?; + Ok(self.crossing_cavities(&mut plane, &ring)?.is_empty()) + } + + /// Recover facet `f`, then register its triangles as constraints. Filling one cavity leaves + /// the others untouched (they share no tetrahedron), so one scan finds them all; the scan is + /// repeated until it finds none. + pub(super) fn recover_facet(&mut self, f: usize) -> Result<()> { + let ring = self.facet_ring(f)?; + let mut plane = self.plane(f)?; + for _ in 0..MAX_ROUNDS_PER_FACET { + let cavities = self.crossing_cavities(&mut plane, &ring)?; + if cavities.is_empty() { + return self.register_facet(f, &ring); + } + for cavity in cavities { + self.stats.cavities += 1; + self.fill_cavity(f, &mut plane, &cavity)?; + } + } + Err(Error::Limit(format!("facet {f} was still crossed after {MAX_ROUNDS_PER_FACET} rounds of cavity filling"))) + } + + /// Record the mesh triangles that make up facet `f` (all three vertices on its boundary). + fn register_facet(&mut self, f: usize, ring: &[u32]) -> Result<()> { + let on: HashSet = ring.iter().copied().collect(); + let mut star = Vec::new(); + for &v in ring { + self.cx.star(v, &mut star)?; + for &t in &star { + let vs = self.cx.tet(t)?; + for k in 0..4 { + let face = face_toward(vs, k); + if face.contains(&v) && face.iter().all(|u| on.contains(u)) { + let key = sorted3(face); + match self.constraint.insert(key, f as u32) { + Some(g) if g != f as u32 => { + return Err(Error::Failed(format!("input triangles {g} and {f} overlap"))); + } + _ => {} + } + } + } + } + } + Ok(()) + } + + /// Delete the cavity and refill its two halves by gift-wrapping. + fn fill_cavity(&mut self, f: usize, plane: &mut Plane, cavity: &[u32]) -> Result<()> { + let inside: HashSet = cavity.iter().copied().collect(); + // Per half: the boundary faces (cavity on their positive side) and the constraint faces + // inside it, and the other sides of those constraints. + let (mut fixed_up, mut fixed_down) = (Vec::new(), Vec::new()); + let (mut back_up, mut back_down) = (Vec::new(), Vec::new()); + let mut vertices: Vec = Vec::new(); + for &t in cavity { + let vs = self.cx.tet(t)?; + vertices.extend(vs); + for (j, n) in self.cx.adj(t)?.into_iter().enumerate() { + let face = face_toward(vs, j); + let interior = n != NONE && inside.contains(&n); + if interior && (n < t || !self.constraint.contains_key(&sorted3(face))) { + continue; + } + let mut pos = false; + let mut neg = false; + for v in face { + match self.facet_side(plane, v)? { + Sign::Positive => pos = true, + Sign::Negative => neg = true, + Sign::Zero => {} + } + } + let (fixed, back) = match (pos, neg) { + (true, false) => (&mut fixed_up, &mut back_up), + (false, true) => (&mut fixed_down, &mut back_down), + _ => { + let at = self.cx.approx(face[0])?; + return Err(Error::Failed(format!("facet {f}: a cavity face near {at:?} is not on one side of the facet's plane"))); + } + }; + fixed.push(face); + if interior { + // A recovered constraint inside the cavity: both sides are open. + let [a, b, c] = face; + back.push([a, c, b]); + } + } + } + vertices.sort_unstable(); + vertices.dedup(); + let mut sides = Vec::with_capacity(vertices.len()); + for &v in &vertices { + sides.push(self.facet_side(plane, v)?); + } + let above: Vec = vertices.iter().zip(&sides).filter(|(_, s)| **s != Sign::Negative).map(|(v, _)| *v).collect(); + let below: Vec = vertices.iter().zip(&sides).filter(|(_, s)| **s != Sign::Positive).map(|(v, _)| *v).collect(); + let on_plane = |v: u32| plane.side.get(&v) == Some(&Sign::Zero); + + let front_up = fixed_up.iter().chain(&back_up).copied().collect(); + let (mut new, open) = self.gift_wrap(f, front_up, &fixed_up, &above, Some(&on_plane))?; + if open.iter().any(|face| !face.iter().all(|&v| on_plane(v))) { + return Err(Error::Failed(format!("facet {f}: the upper half-cavity was left open off the facet"))); + } + // The facet's triangles made by the upper half bound the lower one. + fixed_down.extend(open.iter().copied()); + let front_down = fixed_down.iter().chain(&back_down).copied().collect(); + let (lower, rest) = self.gift_wrap(f, front_down, &fixed_down, &below, None)?; + if !rest.is_empty() { + return Err(Error::Failed(format!("facet {f}: the lower half-cavity did not close"))); + } + new.extend(lower); + let used: HashSet = new.iter().flatten().copied().collect(); + if let Some(&lost) = vertices.iter().find(|v| !used.contains(v)) { + let at = self.cx.approx(lost)?; + return Err(Error::Failed(format!("facet {f}: refilling a cavity lost vertex {lost} at {at:?}"))); + } + self.stats.gift_wrapped_tets += new.len(); + self.cx.replace(cavity, &new)?; + Ok(()) + } + + /// Fill a half-cavity from its open `front` faces (each with the unfilled region on its + /// positive side), using only `vertices`. `fixed` are the triangles bounding it, which no new + /// tetrahedron may cross. Faces on the facet's plane are left open when `skip` says so. The new + /// tetrahedra and the faces left open. + fn gift_wrap( + &mut self, + f: usize, + front: Vec<[u32; 3]>, + fixed: &[[u32; 3]], + vertices: &[u32], + skip: Option<&dyn Fn(u32) -> bool>, + ) -> Result<(Vec<[u32; 4]>, Vec<[u32; 3]>)> { + let mut faces: Vec<[u32; 3]> = Vec::with_capacity(front.len() * 2); + let mut boxes: Vec<([f64; 3], [f64; 3])> = Vec::with_capacity(front.len() * 2); + let mut alive: Vec = Vec::new(); + let mut index: HashMap<[u32; 3], usize> = HashMap::new(); + let mut queue: Vec = Vec::new(); + for face in front { + if index.insert(canon(face), faces.len()).is_some() { + return Err(corrupt("a cavity face is listed twice")); + } + queue.push(faces.len()); + faces.push(face); + boxes.push(self.cx.bounds_of(&face)); + alive.push(true); + } + queue.reverse(); + let fixed_boxes: Vec<([u32; 3], Bounds)> = fixed.iter().map(|&t| (t, self.cx.bounds_of(&t))).collect(); + let cap = 64 * (faces.len() + vertices.len()) + 64; + let mut tets: Vec<[u32; 4]> = Vec::new(); + while let Some(i) = queue.pop() { + if !alive.get(i).copied().unwrap_or(false) { + continue; + } + let sigma = *faces.get(i).ok_or_else(|| corrupt("front index"))?; + if let Some(on) = skip + && sigma.iter().all(|&v| on(v)) + { + continue; + } + let w = self.apex(f, sigma, (&faces, &boxes, &alive), &fixed_boxes, vertices)?; + let tet = [sigma[0], sigma[1], sigma[2], w]; + tets.push(tet); + if tets.len() > cap { + return Err(Error::Limit(format!("facet {f}: gift-wrapping a cavity made more than {cap} tetrahedra"))); + } + if let Some(a) = alive.get_mut(i) { + *a = false; + } + index.remove(&canon(sigma)); + for k in 0..3 { + let g = face_toward(tet, k); + if let Some(j) = index.remove(&canon(g)) { + if let Some(a) = alive.get_mut(j) { + *a = false; + } + continue; + } + let r = [g[0], g[2], g[1]]; + if index.insert(canon(r), faces.len()).is_some() { + let at = self.cx.approx(r[0])?; + return Err(Error::Failed(format!("facet {f}: gift-wrapping overlapped itself near {at:?}"))); + } + queue.push(faces.len()); + faces.push(r); + boxes.push(self.cx.bounds_of(&r)); + alive.push(true); + } + } + let open: Vec<[u32; 3]> = faces.iter().zip(&alive).filter(|(_, a)| **a).map(|(f, _)| *f).collect(); + Ok((tets, open)) + } + + /// The apex for front face `sigma`: conditions (i)–(iii) of the module docs. + fn apex(&mut self, f: usize, sigma: [u32; 3], front: Front<'_>, fixed: &[([u32; 3], Bounds)], vertices: &[u32]) -> Result { + let mut candidates = Vec::new(); + for &w in vertices { + if !sigma.contains(&w) && self.cx.volume([sigma[0], sigma[1], sigma[2], w])? == Sign::Positive { + candidates.push(w); + } + } + let mut excluded = vec![false; candidates.len()]; + for _ in 0..candidates.len() { + let mut best: Option<(usize, u32)> = None; + for (k, &w) in candidates.iter().enumerate() { + if excluded.get(k).copied().unwrap_or(true) { + continue; + } + best = match best { + None => Some((k, w)), + Some((bk, b)) => { + if self.cx.in_sphere([sigma[0], sigma[1], sigma[2], b], w)? { + Some((k, w)) + } else { + Some((bk, b)) + } + } + }; + } + let Some((k, w)) = best else { break }; + if self.valid(sigma, w, front, fixed, vertices)? { + return Ok(w); + } + if let Some(e) = excluded.get_mut(k) { + *e = true; + } + } + let at = self.cx.approx(sigma[0])?; + Err(Error::Failed(format!("facet {f}: gift-wrapping found no valid apex for a cavity face near {at:?}"))) + } + + /// Condition (ii): the tetrahedron `sigma w` meets every front and fixed triangle properly + /// and contains no other cavity vertex. + fn valid(&self, sigma: [u32; 3], w: u32, (faces, boxes, alive): Front<'_>, fixed: &[([u32; 3], Bounds)], vertices: &[u32]) -> Result { + let tet = [sigma[0], sigma[1], sigma[2], w]; + let bx = self.cx.bounds_of(&tet); + let points = self.cx.points(); + for &u in vertices { + if !tet.contains(&u) && boxes_meet(bx, self.cx.bounds(u)) && in_closed_tet(points, tet, u)? { + return Ok(false); + } + } + for ((face, fb), _) in faces.iter().zip(boxes).zip(alive).filter(|(_, a)| **a) { + if boxes_meet(bx, *fb) && !meet_properly(points, &tet, face)? { + return Ok(false); + } + } + for (face, fb) in fixed { + if boxes_meet(bx, *fb) && !meet_properly(points, &tet, face)? { + return Ok(false); + } + } + Ok(true) + } +} diff --git a/crates/mesh/src/cdt/meet.rs b/crates/mesh/src/cdt/meet.rs new file mode 100644 index 0000000..51ff77c --- /dev/null +++ b/crates/mesh/src/cdt/meet.rs @@ -0,0 +1,159 @@ +//! Whether two closed simplices meet properly: their intersection is exactly the simplex spanned +//! by their shared vertices (empty when they share none). The gift-wrapping step of the facet +//! recovery accepts a new tetrahedron only if it meets every cavity and front triangle properly +//! (Diazzi et al. 2023, §4.4 condition ii and App. B), and it must be exact on the coplanar +//! configurations CAD boundaries are made of. +//! +//! # Method +//! Two closed convex polytopes `A`, `B` whose relative interiors are disjoint are weakly separated +//! by a hyperplane `H` (`A` on one closed side, `B` on the other), and then `A ∩ B = (A ∩ H) ∩ +//! (B ∩ H)`. The separating hyperplanes form a pointed polyhedral cone (in `(n, c)` space) whose +//! extreme rays are hyperplanes through `d` affinely independent points of `A ∪ B` (`d` the +//! dimension), so trying the hyperplanes through `d` of the vertices finds one whenever one +//! exists. When none exists the relative interiors overlap and the intersection is not a common +//! face. When one exists, every one gives the same `A ∩ B`, and the question recurses into `H` +//! (one dimension lower) with the vertices on it; a point is the base case. Before that, if all +//! of one simplex's vertices are shared, that simplex is a face of the other and they meet +//! properly. +//! +//! Every test is an exact predicate: `orient3d` in space, `orient2d` projected along an axis the +//! plane is not parallel to, and coordinate comparison along an axis the line is not +//! perpendicular to. + +use crate::predicates::{Axis, GenericPoint, Sign, generic_compare, generic_orient2d, generic_orient3d}; +use crate::{Error, Result}; + +use super::complex::face_toward; + +#[derive(Debug, Clone, Copy)] +enum Level { + Space, + /// A plane, compared by orientation in its projection along this axis. + Plane(Axis), + /// A line, compared along this axis. + Line(Axis), + Point, +} + +fn pt(p: &[GenericPoint], v: u32) -> Result<&GenericPoint> { + p.get(v as usize).ok_or_else(|| Error::Failed("intersection test: vertex id out of range".into())) +} + +/// Whether the closed simplices with vertex ids `a` and `b` (each affinely independent, ids of +/// distinct points) intersect exactly in the hull of their shared vertices. +pub(super) fn meet_properly(p: &[GenericPoint], a: &[u32], b: &[u32]) -> Result { + proper(p, Level::Space, a, b) +} + +/// Whether `x` is in the closed positive tetrahedron `t` (Diazzi et al. 2023, Alg. 3). +pub(super) fn in_closed_tet(p: &[GenericPoint], t: [u32; 4], x: u32) -> Result { + for k in 0..4 { + let [u, v, w] = face_toward(t, k); + // The tetrahedron is on the face's positive side (tet4 volume, i.e. orient3d negative). + if generic_orient3d(pt(p, u)?, pt(p, v)?, pt(p, w)?, pt(p, x)?) == Sign::Positive { + return Ok(false); + } + } + Ok(true) +} + +fn subset(x: &[u32], y: &[u32]) -> bool { + x.iter().all(|v| y.contains(v)) +} + +fn proper(p: &[GenericPoint], level: Level, a: &[u32], b: &[u32]) -> Result { + if subset(a, b) || subset(b, a) { + return Ok(true); + } + // Beyond the subset rule, a point holds one vertex, which both sides would share. + let k = match level { + Level::Space => 3, + Level::Plane(_) => 2, + Level::Line(_) => 1, + Level::Point => return Ok(true), + }; + let mut all: Vec = a.to_vec(); + all.extend(b.iter().filter(|v| !a.contains(v))); + // Candidate hyperplanes through k of the points, those spanned by one simplex first (they + // separate most pairs). + let mut candidates: Vec<[u32; 3]> = Vec::new(); + combinations(&all, k, &mut [0; 3], 0, 0, &mut candidates); + candidates.sort_by_key(|h| !(subset(&h[..k], a) || subset(&h[..k], b))); + let mut spans = false; + for h in &candidates { + let h = &h[..k]; + if let (Level::Space, &[u, v, w]) = (level, h) + && plane_axis(p, [u, v, w])?.is_none() + { + continue; // collinear: no plane + } + spans = true; + let mut signs: Vec<(u32, Sign)> = Vec::with_capacity(all.len()); + for &x in &all { + signs.push((x, if h.contains(&x) { Sign::Zero } else { side(p, level, h, x)? })); + } + let sign_of = |v: &u32| signs.iter().find(|(x, _)| x == v).map_or(Sign::Zero, |(_, s)| *s); + let has = |set: &[u32], s: Sign| set.iter().any(|v| sign_of(v) == s); + let separated = (!has(a, Sign::Negative) && !has(b, Sign::Positive)) || (!has(a, Sign::Positive) && !has(b, Sign::Negative)); + if !separated { + continue; + } + let on = |set: &[u32]| -> Vec { set.iter().copied().filter(|v| sign_of(v) == Sign::Zero).collect() }; + let lower = match (level, h) { + (Level::Space, &[u, v, w]) => { + Level::Plane(plane_axis(p, [u, v, w])?.ok_or_else(|| Error::Failed("intersection test: degenerate plane".into()))?) + } + (Level::Plane(_), &[u, v]) => Level::Line(line_axis(p, u, v)?), + _ => Level::Point, + }; + return proper(p, lower, &on(a), &on(b)); + } + if !spans && let (Level::Space, &[u, v, ..]) = (level, all.as_slice()) { + // All points collinear: compare along their line. + return proper(p, Level::Line(line_axis(p, u, v)?), a, b); + } + Ok(false) +} + +/// An axis the plane through `u v w` is not parallel to (`None`: the points are collinear). +fn plane_axis(p: &[GenericPoint], [u, v, w]: [u32; 3]) -> Result> { + let (a, b, c) = (pt(p, u)?, pt(p, v)?, pt(p, w)?); + Ok(Axis::ALL.into_iter().find(|&ax| !generic_orient2d(a, b, c, ax).is_zero())) +} + +/// An axis along which two distinct points differ. +fn line_axis(p: &[GenericPoint], u: u32, v: u32) -> Result { + let (a, b) = (pt(p, u)?, pt(p, v)?); + Axis::ALL + .into_iter() + .find(|&ax| generic_compare(a, b, ax).is_ne()) + .ok_or_else(|| Error::Failed("intersection test: two vertices are the same point".into())) +} + +fn side(p: &[GenericPoint], level: Level, h: &[u32], x: u32) -> Result { + let px = pt(p, x)?; + Ok(match (level, h) { + (Level::Space, &[u, v, w]) => generic_orient3d(pt(p, u)?, pt(p, v)?, pt(p, w)?, px), + (Level::Plane(ax), &[u, v]) => generic_orient2d(pt(p, u)?, pt(p, v)?, px, ax), + (Level::Line(ax), &[u]) => match generic_compare(px, pt(p, u)?, ax) { + std::cmp::Ordering::Less => Sign::Negative, + std::cmp::Ordering::Equal => Sign::Zero, + std::cmp::Ordering::Greater => Sign::Positive, + }, + _ => return Err(Error::Failed("intersection test: malformed hyperplane".into())), + }) +} + +/// The `k`-subsets of `all` (`k ≤ 3`), each in the first `k` slots. +fn combinations(all: &[u32], k: usize, cur: &mut [u32; 3], len: usize, from: usize, out: &mut Vec<[u32; 3]>) { + if len == k { + out.push(*cur); + return; + } + for i in from..all.len() { + if let (Some(&v), Some(slot)) = (all.get(i), cur.get_mut(len)) { + *slot = v; + combinations(all, k, cur, len + 1, i + 1, out); + } + } +} diff --git a/crates/mesh/src/cdt/mod.rs b/crates/mesh/src/cdt/mod.rs new file mode 100644 index 0000000..a653085 --- /dev/null +++ b/crates/mesh/src/cdt/mod.rs @@ -0,0 +1,361 @@ +//! Constrained Delaunay tetrahedralization of a closed boundary surface: every input triangle is +//! a union of tetrahedron faces, with no tolerance anywhere in the decisions. +//! +//! # Method +//! L. Diazzi, D. Panozzo, A. Vaxman, M. Attene, *Constrained Delaunay Tetrahedrization: A Robust +//! and Practical Approach*, ACM TOG 42(6):181, 2023 (arXiv 2309.09805; paper only, the authors' +//! code is GPL/LGPL and is not read), on the theory of H. Si, K. Gärtner, *Meshing Piecewise +//! Linear Complexes by Constrained Delaunay Tetrahedralizations*, IMR 2005, and J. R. Shewchuk, +//! *Constrained Delaunay Tetrahedralizations and Provably Good Boundary Recovery*, IMR 2002 (also +//! papers only; TetGen is AGPL). +//! +//! 1. The Delaunay tetrahedralization ([`crate::delaunay`]) of the surface's vertices and the +//! eight corners of a box around them, so every Steiner point lies strictly inside the hull. +//! 2. **Segment recovery** ([`segments`], §3.3): every edge of the surface is a segment; a missing +//! one is split at an implicit Steiner point (an [`Lnc`](crate::predicates::Lnc) point on the +//! original segment, §4.3) placed by Si's rule and inserted by Bowyer–Watson with the same +//! perturbed in-sphere test as the Delaunay step, until all segments are edges. +//! 3. **Facet recovery** ([`facets`], §3.4, §4.4): every input triangle is a facet; the +//! tetrahedra crossing a missing facet are deleted and the two half-cavities refilled by the +//! modified gift-wrapping, which builds the constrained Delaunay tetrahedralization of each +//! half with exact predicates on implicit points. +//! 4. **Inside/outside** (§4.5): a flood fill from the box corners that flips at every facet +//! triangle; the orientation of every facet triangle against its inside tetrahedron is then +//! checked exactly. +//! +//! Every decision is an exact predicate ([`crate::predicates`]): `orient3d`, `orient2d` and +//! coordinate comparison on explicit and LNC points, and `insphere` with the Devillers–Teillaud +//! symbolic perturbation. The result is checked before it is returned: every segment is an edge, +//! no edge crosses any facet, every tetrahedron is positive (exactly), and the labelling is +//! consistent with the surface's orientation. +//! +//! # Output +//! [`TetMesh`]: the inside tetrahedra in tet4 order, the surface's vertices first (same ids) and +//! then the Steiner points, each with its segment and parameter, rounded to doubles; and for each +//! boundary tetrahedron face the input triangle it lies on (so the B-rep face tag carries over). +//! Rounding a Steiner point to the nearest doubles moves it off its segment by under an ulp and +//! can, on a nearly flat tetrahedron, flip it (§5.1). [`round`] then tries the doubles a few ulps +//! around such a point for one that keeps every incident tetrahedron positive; [`TetMesh::gate`] +//! runs the validity gate on the rounded mesh and any failure left is reported, never hidden. +//! +//! # Determinism and bounds +//! Same surface → same mesh: the Delaunay step is deterministic, segments and facets are +//! processed in a fixed order, and hash maps are only looked up. Steiner points, tetrahedra, +//! passes and cavities are capped ([`CdtOptions`]); hitting a cap is an [`Error::Limit`]. + +mod complex; +mod facets; +mod meet; +mod round; +mod segments; + +#[cfg(test)] +mod tests; + +use std::collections::HashMap; + +use modelcraft_model::ElementKind; + +use crate::delaunay::Triangulation; +use crate::predicates::{GenericPoint, Point3, Sign}; +use crate::quality::{GateOptions, GateReport, validity_gate_elements}; +use crate::surface::{BoundarySurface, Box3}; +use crate::{Error, Result}; + +use complex::{Complex, NONE, corrupt, face_toward, sorted3}; +use segments::Segment; + +/// Caps on the work of [`tetrahedralize`]. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub struct CdtOptions { + /// Most Steiner points added on segments. A boundary whose triangles are long next to + /// features much closer than their length (needle triangles) needs Steiner points at that + /// spacing all along them; the default is well above what the corpus needs (30k at most) and + /// stops such input in seconds. + pub max_steiner: usize, + /// Most tetrahedra (inside and outside, while building). + pub max_tets: usize, +} + +impl Default for CdtOptions { + fn default() -> Self { + CdtOptions { max_steiner: 50_000, max_tets: 20_000_000 } + } +} + +/// A Steiner point: `points[a] + t (points[b] − points[a])` of the surface edge `[a, b]`. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct SteinerPoint { + pub segment: [u32; 2], + pub t: f64, +} + +/// A tetrahedron face on the boundary: tetrahedron `tet`'s face opposite its vertex +/// `opposite` (0..4) lies on surface triangle `triangle`. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub struct BoundaryFace { + pub tet: u32, + pub opposite: u8, + pub triangle: u32, +} + +impl BoundaryFace { + /// The face's 1-based number in the model's tet4 face order ([`ElementKind::faces`]): the + /// face without the opposite vertex. + pub fn tet4_face(&self) -> Option { + let k = usize::from(self.opposite); + let n = ElementKind::Tet4.faces().iter().position(|f| !f.contains(&k))?; + u8::try_from(n + 1).ok() + } +} + +/// What the recovery did. +#[derive(Debug, Clone, Copy, Default, PartialEq, Eq)] +pub struct CdtStats { + /// Tetrahedra of the initial Delaunay tetrahedralization (with the box). + pub delaunay_tets: usize, + /// Missing subsegments split. + pub missing_subsegments: usize, + /// Cavities refilled by gift-wrapping. + pub cavities: usize, + /// Tetrahedra made by gift-wrapping. + pub gift_wrapped_tets: usize, + /// Tetrahedra in the box at the end. + pub box_tets: usize, + /// Steiner points moved by a few ulps from their nearest doubles to keep every incident + /// tetrahedron positive ([`round`]). + pub rounding_adjusted: usize, +} + +/// The constrained Delaunay tetrahedralization of a closed surface (see the module docs). +#[derive(Debug, Clone, PartialEq)] +pub struct TetMesh { + /// The surface's vertices (same ids as [`BoundarySurface::points`]), then the Steiner + /// points rounded to doubles. + pub points: Vec<[f64; 3]>, + /// The Steiner points, for `points[surface vertices..]`. + pub steiner: Vec, + /// The inside tetrahedra, tet4 order (positive volume in exact arithmetic). + pub tets: Vec<[u32; 4]>, + /// Every tetrahedron face on the surface, with its input triangle. + pub boundary: Vec, + pub stats: CdtStats, +} + +impl TetMesh { + /// The validity gate on the rounded coordinates: every tetrahedron positive, and the volume + /// equal to `expected_volume` (the surface's, [`BoundarySurface::volume`]) within 1e-9. + pub fn gate(&self, expected_volume: f64) -> Result { + let coords = |t: &[u32; 4]| -> Result<[[f64; 3]; 4]> { + let p = |v: u32| self.points.get(v as usize).copied().ok_or_else(|| corrupt("output vertex id")); + Ok([p(t[0])?, p(t[1])?, p(t[2])?, p(t[3])?]) + }; + let elements = self.tets.iter().enumerate().map(|(i, t)| Ok((i as u32 + 1, ElementKind::Tet4, coords(t)?))).collect::>>()?; + validity_gate_elements(elements, &GateOptions { expected_volume: Some(expected_volume), volume_tolerance: 1e-9 }) + } +} + +/// The state of one recovery. +pub(crate) struct Recovery { + cx: Complex, + /// The surface's triangles (the facets), vertex ids as in the complex. + triangles: Vec<[u32; 3]>, + segments: Vec, + /// Sorted edge → segment. + edge_segment: HashMap<(u32, u32), usize>, + /// Recovered facet triangles (sorted vertex ids) → facet. + constraint: HashMap<[u32; 3], u32>, + /// The surface's vertices are the complex's first `surface` vertices, then 8 box corners, + /// then the Steiner points. + surface: usize, + /// Steiner points in creation order: (segment, t). + steiner: Vec<(usize, f64)>, + options: CdtOptions, + stats: CdtStats, +} + +/// The constrained Delaunay tetrahedralization of `surface` (see the module docs). +pub fn tetrahedralize(surface: &BoundarySurface, options: &CdtOptions) -> Result { + let mut rec = Recovery::new(surface, options)?; + rec.classify_acute()?; + rec.recover_segments()?; + for f in 0..rec.triangles.len() { + rec.recover_facet(f)?; + } + rec.verify()?; + let inside = rec.label()?; + rec.output(&inside) +} + +impl Recovery { + fn new(surface: &BoundarySurface, options: &CdtOptions) -> Result { + let n = surface.points().len(); + let triangles = surface.triangles().to_vec(); + if triangles.is_empty() { + return Err(Error::Invalid("the surface has no triangles".into())); + } + let mut points: Vec<[f64; 3]> = surface.points().iter().map(|p| p.map(|x| x + 0.0)).collect(); + let mut bbox = Box3::EMPTY; + for p in &points { + bbox.add_point(*p); + } + let size = (0..3).map(|k| bbox.hi[k] - bbox.lo[k]).fold(0.0, f64::max).max(1.0); + for corner in 0..8 { + points.push(std::array::from_fn(|k| { + let x = if corner >> k & 1 == 0 { bbox.lo[k] - size } else { bbox.hi[k] + size }; + crate::predicates::snap_to_domain(x) + 0.0 + })); + } + let tri = Triangulation::new(&points)?; + if tri.dimension() != 3 || (0..points.len() as u32).any(|v| tri.duplicate_of(v).is_some()) { + return Err(Error::Invalid("the surface's vertices must be distinct".into())); + } + let generic = points.iter().map(|&p| Point3::try_from(p).map(GenericPoint::Explicit)).collect::>>()?; + let cx = Complex::from_delaunay(&tri, generic, options.max_tets)?; + let mut edges: Vec<(u32, u32)> = + triangles.iter().flat_map(|t| (0..3).map(move |k| (t[k].min(t[(k + 1) % 3]), t[k].max(t[(k + 1) % 3])))).collect(); + edges.sort_unstable(); + edges.dedup(); + if let Some(&(a, b)) = edges.iter().find(|(a, b)| a == b || *b as usize >= n) { + return Err(Error::Invalid(format!("surface edge {a}–{b} is degenerate or has no vertex"))); + } + let segments: Vec = edges.iter().map(|&(a, b)| Segment::new(a, b)).collect(); + let edge_segment = edges.iter().enumerate().map(|(i, &e)| (e, i)).collect(); + let stats = CdtStats { delaunay_tets: cx.live_count(), ..CdtStats::default() }; + Ok(Recovery { cx, triangles, segments, edge_segment, constraint: HashMap::new(), surface: n, steiner: Vec::new(), options: *options, stats }) + } + + /// Every segment is an edge, no edge crosses a facet, every facet has triangles, every + /// tetrahedron is positive. + fn verify(&mut self) -> Result<()> { + let mut scratch = Vec::new(); + for s in 0..self.segments.len() { + let chain: Vec = self.segment(s)?.splits.iter().map(|&(_, v)| v).collect(); + for w in chain.windows(2) { + if let &[u, v] = w + && !self.cx.has_edge(u, v, &mut scratch)? + { + let at = self.cx.approx(u)?; + return Err(Error::Failed(format!("a boundary edge near {at:?} was lost during facet recovery"))); + } + } + } + let mut has_faces = vec![false; self.triangles.len()]; + for &f in self.constraint.values() { + if let Some(h) = has_faces.get_mut(f as usize) { + *h = true; + } + } + for f in 0..self.triangles.len() { + if !has_faces.get(f).copied().unwrap_or(false) || !self.facet_present(f)? { + return Err(Error::Failed(format!("input triangle {f} is not a union of tetrahedron faces"))); + } + } + let live: Vec<[u32; 4]> = self.cx.live_tets().map(|(_, vs)| vs).collect(); + for vs in live { + if self.cx.volume(vs)? != Sign::Positive { + let at = self.cx.approx(vs[0])?; + return Err(corrupt(&format!("a tetrahedron near {at:?} is not positive"))); + } + } + Ok(()) + } + + /// Inside (`true`) or outside, per tetrahedron slot: a flood fill from the box corners that + /// flips across every facet triangle (Diazzi et al. §4.5), then the exact check that each + /// facet triangle has its inside tetrahedron on the side its normal points away from. + fn label(&mut self) -> Result> { + let mut label: Vec> = vec![None; self.cx.slots()]; + let corners = self.surface..self.surface + 8; + let corner = self.cx.live_tets().find(|(_, vs)| vs.iter().any(|&v| corners.contains(&(v as usize)))); + let (start, _) = corner.ok_or_else(|| corrupt("no tetrahedron touches the bounding box"))?; + if let Some(l) = label.get_mut(start as usize) { + *l = Some(false); + } + let mut queue = vec![start]; + while let Some(t) = queue.pop() { + let inside = label.get(t as usize).copied().flatten().unwrap_or(false); + let vs = self.cx.tet(t)?; + for (k, nb) in self.cx.adj(t)?.into_iter().enumerate() { + let want = inside != self.constraint.contains_key(&sorted3(face_toward(vs, k))); + match label.get_mut(nb as usize) { + None => {} // the hull + Some(l @ None) => { + *l = Some(want); + queue.push(nb); + } + Some(Some(got)) if *got != want => { + let at = self.cx.approx(vs[k])?; + return Err(Error::Failed(format!( + "inside/outside labelling is inconsistent near {at:?}: the surface does not bound a volume" + ))); + } + Some(_) => {} + } + } + } + let live: Vec<(u32, [u32; 4])> = self.cx.live_tets().collect(); + for (t, vs) in live { + match label.get(t as usize).copied().flatten() { + None => return Err(corrupt("a tetrahedron the flood fill did not reach")), + Some(false) => {} + Some(true) => { + for k in 0..4 { + let Some(&f) = self.constraint.get(&sorted3(face_toward(vs, k))) else { continue }; + let [a, b, c] = *self.triangles.get(f as usize).ok_or_else(|| corrupt("facet id"))?; + let x = vs.get(k).copied().unwrap_or(NONE); + if self.cx.volume([a, b, c, x])? != Sign::Negative { + let at = self.cx.approx(a)?; + return Err(Error::Failed(format!("input triangle {f} near {at:?} faces into the solid: the surface is inverted"))); + } + } + } + } + } + Ok(label.into_iter().map(|l| l == Some(true)).collect()) + } + + fn output(&self, inside: &[bool]) -> Result { + let n = self.surface; + // Surface vertices keep their ids; box corners disappear; Steiner points follow. + let total = self.cx.points().len(); + let mut id = vec![NONE; total]; + let mut points = Vec::with_capacity(n + self.steiner.len()); + for v in 0..n { + if let Some(s) = id.get_mut(v) { + *s = v as u32; + } + points.push(self.cx.approx(v as u32)?); + } + let mut steiner = Vec::with_capacity(self.steiner.len()); + for (k, &(s, t)) in self.steiner.iter().enumerate() { + let v = n + 8 + k; + if let Some(slot) = id.get_mut(v) { + *slot = points.len() as u32; + } + points.push(self.cx.approx(v as u32)?); + steiner.push(SteinerPoint { segment: self.segment(s)?.ends, t }); + } + let mut tets = Vec::new(); + let mut boundary = Vec::new(); + for (t, vs) in self.cx.live_tets() { + if !inside.get(t as usize).copied().unwrap_or(false) { + continue; + } + let mapped = vs.map(|v| id.get(v as usize).copied().unwrap_or(NONE)); + if mapped.contains(&NONE) { + return Err(corrupt("an inside tetrahedron uses a bounding-box corner")); + } + for k in 0..4 { + if let Some(&f) = self.constraint.get(&sorted3(face_toward(vs, k))) { + boundary.push(BoundaryFace { tet: tets.len() as u32, opposite: k as u8, triangle: f }); + } + } + tets.push(mapped); + } + let mut stats = self.stats; + stats.box_tets = self.cx.live_count(); + stats.rounding_adjusted = round::repair(&mut points, &tets, n)?; + Ok(TetMesh { points, steiner, tets, boundary, stats }) + } +} diff --git a/crates/mesh/src/cdt/round.rs b/crates/mesh/src/cdt/round.rs new file mode 100644 index 0000000..1dddce3 --- /dev/null +++ b/crates/mesh/src/cdt/round.rs @@ -0,0 +1,126 @@ +//! Rounding the Steiner points to doubles without inverting tetrahedra. +//! +//! The exact tetrahedralization is valid, but a Steiner point rounded to its nearest doubles is +//! off its segment by up to half an ulp per coordinate, which flips a tetrahedron whose volume is +//! below that scale (Diazzi et al. 2023, §5.1, §6.2: the nearly flat ones, typically four Steiner +//! points on nearly parallel boundary edges). For each tetrahedron that is not positive after +//! rounding, each of its Steiner points tries the doubles within [`REACH`] ulps of its rounded +//! coordinates, nearest first, and takes the first with which every tetrahedron around it is +//! positive (exact `orient3d` on the doubles). Surface vertices never move. The moves are a few +//! ulps, so boundary faces stay on their triangles to rounding precision. The number of passes is +//! capped; what is left is for the validity gate to report. +//! +//! This is the interim answer. The paper's own (§5.1) improves the exact mesh by flips that never +//! cross a facet before rounding, which removes most nearly flat tetrahedra instead of nudging +//! their vertices; it belongs with the flips of the quality stage. + +use crate::predicates::{Point3, Sign, tet_volume_sign}; +use crate::{Error, Result}; + +/// How far, in ulps per coordinate, a rounded Steiner point may move. +const REACH: i32 = 2; +/// Passes over the tetrahedra that are not positive. +const MAX_PASSES: usize = 4; + +fn positive(points: &[[f64; 3]], t: &[u32; 4]) -> Result { + let p = |v: u32| -> Result { + let c = points.get(v as usize).copied().ok_or_else(|| Error::Failed("rounding: vertex id out of range".into()))?; + Point3::try_from(c) + }; + Ok(tet_volume_sign(&p(t[0])?, &p(t[1])?, &p(t[2])?, &p(t[3])?) == Sign::Positive) +} + +/// `x` moved by `k` ulps. +fn step(mut x: f64, k: i32) -> f64 { + for _ in 0..k.unsigned_abs() { + x = if k > 0 { x.next_up() } else { x.next_down() }; + } + x +} + +/// Offsets within [`REACH`], nearest (smallest total) first, in a fixed order. +fn offsets() -> Vec<[i32; 3]> { + let r = -REACH..=REACH; + let mut o: Vec<[i32; 3]> = r.clone().flat_map(|a| r.clone().flat_map(move |b| (-REACH..=REACH).map(move |c| [a, b, c]))).collect(); + o.retain(|d| d != &[0, 0, 0]); + o.sort_by_key(|d| (d.iter().map(|k| k.abs()).sum::(), *d)); + o +} + +/// Move Steiner points (ids `>= surface`) so every tetrahedron is positive where a few-ulp move +/// does it. The number of points moved. +pub(super) fn repair(points: &mut [[f64; 3]], tets: &[[u32; 4]], surface: usize) -> Result { + let mut bad = Vec::new(); + for (i, t) in tets.iter().enumerate() { + if !positive(points, t)? { + bad.push(i); + } + } + if bad.is_empty() { + return Ok(0); + } + let mut around: Vec> = vec![Vec::new(); points.len()]; + for (i, t) in tets.iter().enumerate() { + for &v in t { + if (v as usize) >= surface + && let Some(a) = around.get_mut(v as usize) + { + a.push(i as u32); + } + } + } + let offsets = offsets(); + let mut moved = vec![false; points.len()]; + for _ in 0..MAX_PASSES { + let mut still = Vec::new(); + for &i in &bad { + let t = *tets.get(i).ok_or_else(|| Error::Failed("rounding: tetrahedron id".into()))?; + if positive(points, &t)? { + continue; + } + let mut fixed = false; + for v in t.into_iter().filter(|&v| (v as usize) >= surface) { + let Some(origin) = points.get(v as usize).copied() else { continue }; + let incident = around.get(v as usize).cloned().unwrap_or_default(); + for d in &offsets { + let candidate: [f64; 3] = std::array::from_fn(|k| step(origin[k], d[k])); + if candidate.iter().any(|x| !crate::predicates::is_admissible(*x)) { + continue; + } + if let Some(p) = points.get_mut(v as usize) { + *p = candidate; + } + let mut ok = true; + for &j in &incident { + let u = tets.get(j as usize).ok_or_else(|| Error::Failed("rounding: tetrahedron id".into()))?; + if !positive(points, u)? { + ok = false; + break; + } + } + if ok { + fixed = true; + if let Some(m) = moved.get_mut(v as usize) { + *m = true; + } + break; + } + if let Some(p) = points.get_mut(v as usize) { + *p = origin; + } + } + if fixed { + break; + } + } + if !fixed { + still.push(i); + } + } + if still.is_empty() || still.len() == bad.len() { + break; + } + bad = still; + } + Ok(moved.iter().filter(|m| **m).count()) +} diff --git a/crates/mesh/src/cdt/segments.rs b/crates/mesh/src/cdt/segments.rs new file mode 100644 index 0000000..e5eb9fe --- /dev/null +++ b/crates/mesh/src/cdt/segments.rs @@ -0,0 +1,279 @@ +//! Segment recovery (Si & Gärtner 2005, as restated by Diazzi et al. 2023, §3.3): while an edge +//! of the input surface (a *segment*) is not an edge of the Delaunay tetrahedralization, split +//! it at a Steiner point and insert the point (Bowyer–Watson). Each Steiner point is the LNC +//! point `a + t (b − a)` of its original segment `a b` (Diazzi et al. §4.2–§4.3): it lies on the +//! segment exactly, and only `t` is computed in floating point, so the placement rule may round +//! but never moves a point off its segment. +//! +//! The placement follows Si's reference-point rule. The reference point `r` of a missing +//! subsegment `v₁ v₂` is the vertex inside its diametral sphere for which the circle through +//! `v₁, v₂, r` is largest. A vertex is *acute* when two segments meet there at less than 90°; +//! the rule depends on how many of the original segment's endpoints are acute: +//! - none: split where the smaller of the spheres centred at `v₁`, `v₂` through `r` cuts the +//! segment, or at the midpoint when both radii exceed half the length; +//! - one (`w`): `p` is where the sphere centred at `w` through `r` cuts the segment; keep `p` if +//! it is closer to `r` than to the far endpoint `v_o`, else shift it towards the near endpoint +//! `v_w` by `|r p|` when that is under half of `|v_w p|`, else take the midpoint of `v_w p`; +//! - two: split at the midpoint first, after which each half has one acute endpoint. +//! +//! Any value that falls outside the open subsegment becomes the midpoint. Termination is Si's +//! (the spheres around acute vertices keep splits from cascading); the number of Steiner points +//! and of passes is still capped. + +use crate::predicates::{Axis, GenericPoint, Lnc, Point3, generic_orient2d}; +use crate::vec3::{cross, dot, norm, sub}; +use crate::{Error, Result}; + +use super::Recovery; + +/// Most tetrahedra searched for a reference point. +const NEAR_CAP: usize = 20_000; +/// Most passes over the segments. +const MAX_PASSES: usize = 10_000; + +/// An input segment and its split points. +#[derive(Debug, Clone)] +pub(super) struct Segment { + /// The original endpoints (explicit vertices), `ends[0] < ends[1]`. + pub ends: [u32; 2], + /// `(t, vertex)` along the segment from `ends[0]` (`t = 0`) to `ends[1]` (`t = 1`). + pub splits: Vec<(f64, u32)>, + /// Whether each original endpoint is acute. + pub acute: [bool; 2], +} + +impl Segment { + pub(super) fn new(a: u32, b: u32) -> Self { + let ends = [a.min(b), a.max(b)]; + Segment { ends, splits: vec![(0.0, ends[0]), (1.0, ends[1])], acute: [false; 2] } + } + + /// The vertices from `from` to the other end, in order. + pub(super) fn chain_from(&self, from: u32) -> Vec { + let mut v: Vec = self.splits.iter().map(|&(_, v)| v).collect(); + if from != self.ends[0] { + v.reverse(); + } + v + } +} + +fn dist(a: [f64; 3], b: [f64; 3]) -> f64 { + norm(sub(a, b)) +} + +impl Recovery { + /// Mark acute vertices: two segments at a vertex make an angle under 90° (floating point: + /// this only steers where Steiner points go). + pub(super) fn classify_acute(&mut self) -> Result<()> { + let n = self.cx.points().len(); + let mut dirs: Vec> = vec![Vec::new(); n]; + for s in &self.segments { + let [a, b] = s.ends; + let (pa, pb) = (self.cx.approx(a)?, self.cx.approx(b)?); + if let Some(d) = dirs.get_mut(a as usize) { + d.push(sub(pb, pa)); + } + if let Some(d) = dirs.get_mut(b as usize) { + d.push(sub(pa, pb)); + } + } + let acute: Vec = dirs.iter().map(|d| d.iter().enumerate().any(|(i, x)| d.iter().skip(i + 1).any(|y| dot(*x, *y) > 0.0))).collect(); + for s in &mut self.segments { + s.acute = s.ends.map(|v| acute.get(v as usize).copied().unwrap_or(false)); + } + Ok(()) + } + + /// Split missing segments until every one is an edge of the (Delaunay) tetrahedralization. + /// Alg. 2 sweeps all segments until none is missing; an insertion only changes edges between + /// vertices of the tetrahedra it deletes, so after the first sweep only the segments through + /// those vertices are checked again. + pub(super) fn recover_segments(&mut self) -> Result<()> { + let mut scratch = Vec::new(); + let mut through: Vec> = vec![Vec::new(); self.cx.points().len()]; + for (s, seg) in self.segments.iter().enumerate() { + for v in seg.ends { + if let Some(l) = through.get_mut(v as usize) { + l.push(s as u32); + } + } + } + let mut queued = vec![false; self.segments.len()]; + let mut work: Vec = (0..self.segments.len()).collect(); + for _ in 0..MAX_PASSES { + if work.is_empty() { + return Ok(()); + } + queued.fill(false); + let mut next: Vec = Vec::new(); + let mut requeue = |s: usize, next: &mut Vec| { + if let Some(q) = queued.get_mut(s) + && !*q + { + *q = true; + next.push(s); + } + }; + for &s in &work { + let mut i = 0; + while let (Some(&(tu, u)), Some(&(tv, v))) = (self.segment(s)?.splits.get(i), self.segment(s)?.splits.get(i + 1)) { + if self.cx.has_edge(u, v, &mut scratch)? { + i += 1; + continue; + } + self.stats.missing_subsegments += 1; + let t = self.split_parameter(s, (tu, u), (tv, v))?; + let [a, b] = self.segment(s)?.ends; + let (pa, pb) = (explicit(self.cx.pt(a)?)?, explicit(self.cx.pt(b)?)?); + let w = self.cx.push_point(GenericPoint::Lnc(Lnc::new(&pa, &pb, t)?))?; + through.push(vec![s as u32]); + self.steiner.push((s, t)); + if self.steiner.len() > self.options.max_steiner { + let at = self.cx.approx(u)?; + return Err(Error::Limit(format!( + "more than {} Steiner points on boundary edges (the last near {at:?}): the boundary has features much thinner than its triangles are long; raise the cap or tessellate coarser", + self.options.max_steiner + ))); + } + let touched = self.cx.insert(w, u)?; + if let Some(seg) = self.segments.get_mut(s) { + seg.splits.insert(i + 1, (t, w)); + } + requeue(s, &mut next); + for x in touched { + for &r in through.get(x as usize).map(Vec::as_slice).unwrap_or_default() { + requeue(r as usize, &mut next); + } + } + // Both halves wait for the next sweep (Alg. 2 splits each missing segment + // once per sweep). + i += 2; + } + } + work = next; + } + Err(Error::Limit(format!("segment recovery did not finish in {MAX_PASSES} passes"))) + } + + pub(super) fn segment(&self, s: usize) -> Result<&Segment> { + self.segments.get(s).ok_or_else(|| super::complex::corrupt("segment id")) + } + + /// Si's split parameter for the missing subsegment `(tu, u)`–`(tv, v)` of segment `s`. + fn split_parameter(&mut self, s: usize, (tu, u): (f64, u32), (tv, v): (f64, u32)) -> Result { + let (ends, acute) = { + let seg = self.segment(s)?; + (seg.ends, seg.acute) + }; + let [a, b] = ends; + let (pa, pb) = (self.cx.approx(a)?, self.cx.approx(b)?); + let length = dist(pa, pb); + let mid = 0.5 * (tu + tv); + let r = self.reference_point(u, v)?; + let chosen = match (r, acute) { + (None, _) => mid, + // Two acute ends: the first split is the midpoint. + (_, [true, true]) if tu < 0.5 && tv > 0.5 => 0.5, + (Some(r), acute) => { + let pr = self.cx.approx(r)?; + let (pu, pv) = (self.cx.approx(u)?, self.cx.approx(v)?); + match acute { + [false, false] => { + let (ru, rv) = (dist(pu, pr), dist(pv, pr)); + let half = 0.5 * dist(pu, pv); + if ru > half && rv > half { + mid + } else if ru <= rv { + tu + ru / length + } else { + tv - rv / length + } + } + _ => { + // The acute end w: the only one, or the end of the half this + // subsegment lies in. + let w_is_a = match acute { + [true, false] => true, + [false, true] => false, + _ => tv <= 0.5, + }; + let pw = if w_is_a { pa } else { pb }; + let radius = dist(pw, pr); + let tp = if w_is_a { radius / length } else { 1.0 - radius / length }; + let (t_near, t_far) = if w_is_a { (tu, tv) } else { (tv, tu) }; + let at = |t: f64| -> [f64; 3] { std::array::from_fn(|k| pa[k] + t * (pb[k] - pa[k])) }; + let p = at(tp); + if dist(p, pr) < dist(p, at(t_far)) { + tp + } else { + let d = dist(pr, p); + let near = dist(at(t_near), p); + if d < 0.5 * near { if t_near < tp { tp - d / length } else { tp + d / length } } else { 0.5 * (t_near + tp) } + } + } + } + } + }; + let t = if chosen.is_finite() && chosen > tu && chosen < tv { chosen } else { mid }; + if !(t > tu && t < tv) { + let at = self.cx.approx(u)?; + return Err(Error::Failed(format!("a missing segment near {at:?} is too short to split further in double precision"))); + } + Ok(t) + } + + /// The reference point of the missing subsegment `u v`: the vertex strictly inside its + /// diametral sphere with the largest circumcircle through `u, v` (floating point; `None` if + /// rounding hides every encroaching vertex). Candidates are the vertices of the tetrahedra + /// whose boxes meet the subsegment's box, connected to `u`'s star: they include every + /// tetrahedron the subsegment crosses, whose vertices are the ones keeping it from being an + /// edge. (On the corpus this picks the same point as searching the whole diametral ball, at a + /// fraction of the cost.) A vertex exactly on the open subsegment makes the + /// input invalid (a segment through a vertex), which is reported. + fn reference_point(&mut self, u: u32, v: u32) -> Result> { + let (pu, pv) = (self.cx.approx(u)?, self.cx.approx(v)?); + let c: [f64; 3] = std::array::from_fn(|k| 0.5 * (pu[k] + pv[k])); + let r2 = 0.25 * dot(sub(pu, pv), sub(pu, pv)); + let lo: [f64; 3] = std::array::from_fn(|k| pu[k].min(pv[k])); + let hi: [f64; 3] = std::array::from_fn(|k| pu[k].max(pv[k])); + let near = self.cx.tets_near(u, lo, hi, NEAR_CAP)?; + let mut seen: std::collections::HashSet = std::collections::HashSet::new(); + let mut best: Option<(f64, u32)> = None; + for t in near { + for x in self.cx.tet(t)? { + if x == u || x == v || !seen.insert(x) { + continue; + } + let px = self.cx.approx(x)?; + let d = sub(px, c); + if dot(d, d) >= r2 { + continue; + } + let (ea, eb) = (sub(pu, px), sub(pv, px)); + let area2 = norm(cross(ea, eb)); + let radius = if area2 > 0.0 { norm(ea) * norm(eb) * norm(sub(ea, eb)) / (2.0 * area2) } else { f64::INFINITY }; + if !radius.is_finite() && self.on_segment(u, v, x)? { + let at = self.cx.approx(x)?; + return Err(Error::Invalid(format!("vertex {x} at {at:?} lies on a boundary edge: the surface is not a valid complex"))); + } + if best.is_none_or(|(b, _)| radius > b) { + best = Some((radius, x)); + } + } + } + Ok(best.map(|(_, x)| x)) + } + + /// Whether `x` is exactly collinear with `u v` (it is inside their diametral ball, so then on + /// the open segment). + fn on_segment(&self, u: u32, v: u32, x: u32) -> Result { + let (a, b, c) = (self.cx.pt(u)?, self.cx.pt(v)?, self.cx.pt(x)?); + Ok(Axis::ALL.into_iter().all(|ax| generic_orient2d(a, b, c, ax).is_zero())) + } +} + +/// The explicit point of an original segment endpoint. +fn explicit(p: &GenericPoint) -> Result { + p.explicit().copied().ok_or_else(|| super::complex::corrupt("a segment endpoint is not an input vertex")) +} diff --git a/crates/mesh/src/cdt/tests.rs b/crates/mesh/src/cdt/tests.rs new file mode 100644 index 0000000..fdbb4c8 --- /dev/null +++ b/crates/mesh/src/cdt/tests.rs @@ -0,0 +1,320 @@ +//! Tests of the constrained Delaunay tetrahedralization: exact conformity to every input +//! triangle, the enclosed volume, the validity gate, and the classic polyhedra that have no +//! tetrahedralization without Steiner points. + +use super::meet::{in_closed_tet, meet_properly}; +use super::*; +use crate::predicates::{GenericPoint, Lnc, Point3}; +use crate::surface::{BoundarySurface, SurfaceOptions}; +use crate::vec3::{cross, dot, norm, sub}; + +/// A closed surface from points and triangles, each triangle oriented outward from `centre` +/// (the shapes here are star-shaped about it). +fn surface(points: &[[f64; 3]], triangles: &[[u32; 3]], centre: [f64; 3]) -> BoundarySurface { + let oriented: Vec<[u32; 3]> = triangles + .iter() + .map(|&[a, b, c]| { + let (pa, pb, pc) = (points[a as usize], points[b as usize], points[c as usize]); + let n = cross(sub(pb, pa), sub(pc, pa)); + if dot(n, sub(pa, centre)) > 0.0 { [a, b, c] } else { [a, c, b] } + }) + .collect(); + let faces: Vec = (0..oriented.len() as u32).collect(); + BoundarySurface::new(points, &oriented, &faces, &SurfaceOptions::default()).unwrap() +} + +/// The checks every result must pass: the gate (validity, volume), every input triangle covered +/// by boundary faces of exactly its area, each boundary face on its triangle's plane. +fn check(s: &BoundarySurface, m: &TetMesh) { + let report = m.gate(s.volume()).unwrap(); + assert!(report.passed(), "{report:?}"); + let mut area = vec![0.0; s.triangles().len()]; + for b in &m.boundary { + let t = m.tets[b.tet as usize]; + let face = face_toward(t, b.opposite as usize).map(|v| m.points[v as usize]); + area[b.triangle as usize] += 0.5 * norm(cross(sub(face[1], face[0]), sub(face[2], face[0]))); + let tri = s.triangles()[b.triangle as usize].map(|v| s.points()[v as usize]); + let tn = cross(sub(tri[1], tri[0]), sub(tri[2], tri[0])); + for p in face { + let d = dot(tn, sub(p, tri[0])) / norm(tn); + assert!(d.abs() < 1e-9, "{d}"); + } + } + for (i, t) in s.triangles().iter().enumerate() { + let p = t.map(|v| s.points()[v as usize]); + let want = 0.5 * norm(cross(sub(p[1], p[0]), sub(p[2], p[0]))); + assert!((area[i] - want).abs() <= 1e-9 * want.max(1e-300), "triangle {i}: {} vs {want}", area[i]); + } +} + +fn cube() -> (Vec<[f64; 3]>, Vec<[u32; 3]>) { + let p: Vec<[f64; 3]> = (0..8).map(|i| [(i & 1) as f64, (i >> 1 & 1) as f64, (i >> 2 & 1) as f64]).collect(); + let quads = [[0, 1, 3, 2], [4, 5, 7, 6], [0, 1, 5, 4], [2, 3, 7, 6], [0, 2, 6, 4], [1, 3, 7, 5]]; + let t = quads.iter().flat_map(|q| [[q[0], q[1], q[2]], [q[0], q[2], q[3]]]).collect(); + (p, t) +} + +#[test] +fn a_cube_conforms() { + let (p, t) = cube(); + let s = surface(&p, &t, [0.5; 3]); + let m = tetrahedralize(&s, &CdtOptions::default()).unwrap(); + check(&s, &m); + // A square's diagonal is cocircular with the other two corners, so it is not strongly + // Delaunay; at most one split per face (its midpoint) recovers it. + assert!(m.steiner.len() <= 6, "{:?}", m.steiner); + assert_eq!(m.boundary.len(), 12 + 2 * m.steiner.len()); +} + +/// Schönhardt's twisted prism (Schönhardt 1928): no tetrahedralization without Steiner points. +fn schonhardt(twist_deg: f64) -> (Vec<[f64; 3]>, Vec<[u32; 3]>) { + let mut p = Vec::new(); + for z in 0..2 { + for i in 0..3 { + let a = (i as f64) * 2.0 * std::f64::consts::PI / 3.0 + if z == 1 { twist_deg.to_radians() } else { 0.0 }; + p.push([a.cos(), a.sin(), z as f64 * 1.5]); + } + } + let mut t = vec![[0, 1, 2], [3, 4, 5]]; + for i in 0..3u32 { + let j = (i + 1) % 3; + // The reflex diagonal: from bottom i to top j (the top is turned towards j). + t.push([i, j, 3 + j]); + t.push([i, 3 + j, 3 + i]); + } + (p, t) +} + +#[test] +fn schonhardt_prism_gets_steiner_points_and_conforms() { + let (p, t) = schonhardt(30.0); + let s = surface(&p, &t, [0.0, 0.0, 0.75]); + let m = tetrahedralize(&s, &CdtOptions::default()).unwrap(); + check(&s, &m); + assert!(!m.steiner.is_empty(), "Schönhardt's prism has no tetrahedralization without Steiner points"); + for sp in &m.steiner { + assert!(sp.t > 0.0 && sp.t < 1.0); + } +} + +/// A UV sphere: many cospherical vertices. +fn uv_sphere(rings: usize, sectors: usize) -> (Vec<[f64; 3]>, Vec<[u32; 3]>) { + let mut p = vec![[0.0, 0.0, 1.0], [0.0, 0.0, -1.0]]; + for r in 1..rings { + let th = std::f64::consts::PI * r as f64 / rings as f64; + for s in 0..sectors { + let ph = 2.0 * std::f64::consts::PI * s as f64 / sectors as f64; + p.push([th.sin() * ph.cos(), th.sin() * ph.sin(), th.cos()]); + } + } + let id = |r: usize, s: usize| (2 + (r - 1) * sectors + s % sectors) as u32; + let mut t = Vec::new(); + for s in 0..sectors { + t.push([0, id(1, s), id(1, s + 1)]); + t.push([1, id(rings - 1, s + 1), id(rings - 1, s)]); + for r in 1..rings - 1 { + t.push([id(r, s), id(r + 1, s), id(r + 1, s + 1)]); + t.push([id(r, s), id(r + 1, s + 1), id(r, s + 1)]); + } + } + (p, t) +} + +#[test] +fn a_sphere_conforms_and_is_deterministic() { + let (p, t) = uv_sphere(12, 24); + let s = surface(&p, &t, [0.0; 3]); + let a = tetrahedralize(&s, &CdtOptions::default()).unwrap(); + check(&s, &a); + let b = tetrahedralize(&s, &CdtOptions::default()).unwrap(); + assert_eq!(a, b); +} + +#[test] +fn caps_are_reported() { + let (p, t) = schonhardt(30.0); + let s = surface(&p, &t, [0.0, 0.0, 0.75]); + let none = CdtOptions { max_steiner: 0, max_tets: 1_000_000 }; + assert!(matches!(tetrahedralize(&s, &none), Err(Error::Limit(_)))); + let few = CdtOptions { max_steiner: 1000, max_tets: 10 }; + assert!(matches!(tetrahedralize(&s, &few), Err(Error::Limit(_)))); +} + +fn e(x: f64, y: f64, z: f64) -> GenericPoint { + GenericPoint::Explicit(Point3::new(x, y, z).unwrap()) +} + +#[test] +fn proper_meeting_handles_coplanar_and_shared_cases() { + let p = vec![ + e(0.0, 0.0, 0.0), // 0 + e(1.0, 0.0, 0.0), // 1 + e(0.0, 1.0, 0.0), // 2 + e(0.0, 0.0, 1.0), // 3: tet 0 1 2 3 + e(1.0, 1.0, 0.0), // 4: beside face 0 1 2, in its plane + e(0.2, 0.2, 0.0), // 5: inside face 0 1 2 + e(0.0, 0.0, -1.0), // 6: below + e(0.5, 0.5, 2.0), // 7: above, off the tet + e(2.0, 2.0, 2.0), // 8: far + e(0.1, 0.1, 0.1), // 9: inside the tet + ]; + let tet = [0, 1, 2, 3]; + let yes = |b: &[u32]| assert!(meet_properly(&p, &tet, b).unwrap(), "{b:?}"); + let no = |b: &[u32]| assert!(!meet_properly(&p, &tet, b).unwrap(), "{b:?}"); + yes(&[0, 1, 2]); // a face + yes(&[1, 2, 4]); // coplanar, sharing edge 1 2, on the other side + yes(&[0, 1, 6]); // shares an edge, below + yes(&[3, 7, 8]); // shares a vertex, away + yes(&[6, 8, 4]); // disjoint + no(&[0, 1, 5]); // coplanar, overlapping the face 0 1 2 + no(&[1, 2, 5]); // coplanar, overlapping + no(&[6, 7, 4]); // passes through the tet + no(&[3, 9, 7]); // shares a vertex, enters the tet + no(&[5, 6, 8]); // touches the face 0 1 2 at an interior point + assert!(in_closed_tet(&p, tet, 9).unwrap()); + assert!(in_closed_tet(&p, tet, 5).unwrap()); + assert!(!in_closed_tet(&p, tet, 8).unwrap()); + // An LNC point on edge 0 3 is on the tet's boundary. + let mut q = p.clone(); + q.push(GenericPoint::Lnc(Lnc::new(&Point3::new(0.0, 0.0, 0.0).unwrap(), &Point3::new(0.0, 0.0, 1.0).unwrap(), 0.3).unwrap())); + assert!(in_closed_tet(&q, tet, 10).unwrap()); + assert!(!meet_properly(&q, &tet, &[10, 6, 8]).unwrap(), "touches edge 0 3 at a non-shared point"); + // Tetrahedra sharing a face meet properly; overlapping ones do not. + assert!(meet_properly(&p, &tet, &[0, 1, 2, 6]).unwrap()); + assert!(!meet_properly(&p, &tet, &[0, 1, 2, 9]).unwrap()); +} + +/// Two solids: a wide tetrahedron standing on the big triangle z = 0, and a small one whose top +/// vertex is just under that triangle's centre. The Delaunay edge from the big apex to the small +/// top pierces the big triangle while all its edges are present: only facet recovery (the +/// gift-wrapping of a cavity) restores it. +fn pierced() -> (Vec<[f64; 3]>, Vec<[u32; 3]>, Vec<[f64; 3]>) { + let p = vec![ + [-10.0, -10.0, 0.0], + [10.0, -10.0, 0.0], + [0.0, 10.0, 0.0], + [0.0, 0.0, 5.0], + [-1.0, -1.0, -1.0], + [1.0, -1.0, -1.0], + [0.0, 1.0, -1.0], + [0.0, 0.0, -0.01], + ]; + let t = vec![[0, 1, 2], [0, 1, 3], [1, 2, 3], [2, 0, 3], [4, 5, 6], [4, 5, 7], [5, 6, 7], [6, 4, 7]]; + (p, t, vec![[0.0, -3.0, 1.0], [0.0, -0.3, -0.7]]) +} + +/// Orients each triangle outward from the centre of the shell it belongs to. +fn shells(points: &[[f64; 3]], triangles: &[[u32; 3]], centres: &[[f64; 3]], per_shell: usize) -> BoundarySurface { + let oriented: Vec<[u32; 3]> = triangles + .iter() + .enumerate() + .map(|(i, &[a, b, c])| { + let centre = centres[i / per_shell]; + let (pa, pb, pc) = (points[a as usize], points[b as usize], points[c as usize]); + let n = cross(sub(pb, pa), sub(pc, pa)); + if dot(n, sub(pa, centre)) > 0.0 { [a, b, c] } else { [a, c, b] } + }) + .collect(); + let faces: Vec = (0..oriented.len() as u32).collect(); + BoundarySurface::new(points, &oriented, &faces, &SurfaceOptions::default()).unwrap() +} + +#[test] +fn a_pierced_facet_is_recovered_by_gift_wrapping() { + let (p, t, centres) = pierced(); + let s = shells(&p, &t, ¢res, 4); + assert_eq!(s.shells().len(), 2); + let m = tetrahedralize(&s, &CdtOptions::default()).unwrap(); + check(&s, &m); + assert!(m.stats.cavities > 0, "{:?}", m.stats); + assert!(m.stats.gift_wrapped_tets > 0); +} + +/// A Chazelle-style polyhedron (Chazelle 1984): a block whose top is cut by V-notches along `x` +/// and whose bottom is cut by V-notches along `y`, the notch edges almost touching. Its +/// tetrahedralizations need many Steiner points. +fn chazelle(notches: usize, gap: f64) -> (Vec<[f64; 3]>, Vec<[u32; 3]>) { + let n = 2 * notches; + let at = |i: usize| i as f64 / n as f64; + // Top height by y index: ridges at even, valleys at odd. Bottom height by x index. + let top = |j: usize| if j.is_multiple_of(2) { 1.0 } else { 0.5 + gap / 2.0 }; + let bottom = |i: usize| if i.is_multiple_of(2) { 0.0 } else { 0.5 - gap / 2.0 }; + let mut p = Vec::new(); + let mut id = std::collections::HashMap::new(); + let mut vertex = |x: usize, y: usize, up: bool, p: &mut Vec<[f64; 3]>| -> u32 { + *id.entry((x, y, up)).or_insert_with(|| { + p.push([at(x), at(y), if up { top(y) } else { bottom(x) }]); + (p.len() - 1) as u32 + }) + }; + let mut t = Vec::new(); + let quad = |q: [u32; 4], out: [f64; 3], p: &Vec<[f64; 3]>, t: &mut Vec<[u32; 3]>| { + for tri in [[q[0], q[1], q[2]], [q[0], q[2], q[3]]] { + let (a, b, c) = (p[tri[0] as usize], p[tri[1] as usize], p[tri[2] as usize]); + let nrm = cross(sub(b, a), sub(c, a)); + t.push(if dot(nrm, out) > 0.0 { tri } else { [tri[0], tri[2], tri[1]] }); + } + }; + for i in 0..n { + for j in 0..n { + let q = [vertex(i, j, true, &mut p), vertex(i + 1, j, true, &mut p), vertex(i + 1, j + 1, true, &mut p), vertex(i, j + 1, true, &mut p)]; + quad(q, [0.0, 0.0, 1.0], &p, &mut t); + let q = + [vertex(i, j, false, &mut p), vertex(i + 1, j, false, &mut p), vertex(i + 1, j + 1, false, &mut p), vertex(i, j + 1, false, &mut p)]; + quad(q, [0.0, 0.0, -1.0], &p, &mut t); + } + } + for k in 0..n { + for (fixed, out) in [(0, -1.0), (n, 1.0)] { + // Walls x = fixed (along y) and y = fixed (along x). + let q = [ + vertex(fixed, k, false, &mut p), + vertex(fixed, k + 1, false, &mut p), + vertex(fixed, k + 1, true, &mut p), + vertex(fixed, k, true, &mut p), + ]; + quad(q, [out, 0.0, 0.0], &p, &mut t); + let q = [ + vertex(k, fixed, false, &mut p), + vertex(k + 1, fixed, false, &mut p), + vertex(k + 1, fixed, true, &mut p), + vertex(k, fixed, true, &mut p), + ]; + quad(q, [0.0, out, 0.0], &p, &mut t); + } + } + (p, t) +} + +#[test] +fn chazelle_polyhedron_conforms() { + let (p, t) = chazelle(3, 0.02); + let faces: Vec = (0..t.len() as u32).collect(); + let s = BoundarySurface::new(&p, &t, &faces, &SurfaceOptions::default()).unwrap(); + assert_eq!(s.shells().len(), 1); + let m = tetrahedralize(&s, &CdtOptions::default()).unwrap(); + check(&s, &m); + assert!(!m.steiner.is_empty()); +} + +#[test] +fn hostile_and_degenerate_surfaces_are_errors() { + // A flat "solid" (all points coplanar) bounds no volume. + let p = vec![[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [1.0, 1.0, 0.0]]; + let t = vec![[0, 1, 2], [1, 3, 2], [0, 2, 1], [1, 2, 3]]; + let faces = vec![0; 4]; + if let Ok(s) = BoundarySurface::new(&p, &t, &faces, &SurfaceOptions::default()) { + assert!(tetrahedralize(&s, &CdtOptions::default()).is_err()); + } +} + +#[test] +fn the_gate_reports_an_inverted_tetrahedron() { + let (p, t) = cube(); + let s = surface(&p, &t, [0.5; 3]); + let mut m = tetrahedralize(&s, &CdtOptions::default()).unwrap(); + m.tets[0].swap(0, 1); + let report = m.gate(s.volume()).unwrap(); + assert_eq!(report.invalid, vec![1]); + assert!(report.into_result().is_err()); +} diff --git a/crates/mesh/src/lib.rs b/crates/mesh/src/lib.rs index 486bad1..b4e0cd4 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 cdt; pub mod delaunay; pub mod quad; pub mod quality; diff --git a/crates/mesh/src/predicates/indirect.rs b/crates/mesh/src/predicates/indirect.rs index 68bdc73..8908849 100644 --- a/crates/mesh/src/predicates/indirect.rs +++ b/crates/mesh/src/predicates/indirect.rs @@ -1,7 +1,10 @@ //! Indirect predicates on implicit points (Attene 2020, §4–§5), from the paper only. //! //! An implicit point is kept as the explicit points that define it and is never rounded: a -//! line–plane intersection ([`Lpi`], §4.2) or a three-plane intersection ([`Tpi`]). Each has +//! line–plane intersection ([`Lpi`], §4.2), a three-plane intersection ([`Tpi`]), or a linear +//! combination of two points ([`Lnc`], the Steiner points of Diazzi, Panozzo, Vaxman, Attene, +//! *Constrained Delaunay Tetrahedrization: A Robust and Practical Approach*, ACM TOG 42(6), 2023, +//! §4.3). Each has //! homogeneous coordinates `(X, W)` that are polynomials in the defining coordinates, so a //! predicate on implicit points is the sign of a polynomial times the signs of the `W`s //! (Attene §4.3–§4.6, §5.1: "count how many of them have a negative sign"). @@ -19,9 +22,10 @@ use std::cmp::Ordering; use super::dyadic::Dyadic; +use super::eft::{two_diff, two_product, two_sum}; use super::interval::Interval; -use super::ring::{Ring, Signed, det2, det3}; -use super::{Point2, Point3, Sign, orient2d, orient3d}; +use super::ring::{Ring, Signed, det2, det3, sq3}; +use super::{Point2, Point3, Sign, insphere, orient2d, orient3d}; use crate::{Error, Result}; /// A coordinate axis: the axis compared by [`generic_compare`], or dropped by @@ -34,6 +38,9 @@ pub enum Axis { } impl Axis { + /// The three axes, in order. + pub const ALL: [Axis; 3] = [Axis::X, Axis::Y, Axis::Z]; + pub(crate) fn index(self) -> usize { match self { Axis::X => 0, @@ -127,6 +134,15 @@ fn tpi_hom(planes: &[[[f64; 3]; 3]; 3]) -> Hom { Hom { x, w: Some(det3(&[n1, n2, n3])) } } +/// `a + t (b - a)` (Diazzi et al. 2023, §4.3 write `t v₁ + (1 - t) v₂`; this form keeps every +/// operand a double, so it is exact in the ring). `W = 1`. +fn lnc_hom(a: &[f64; 3], b: &[f64; 3], t: f64) -> Hom { + let (a, b, t) = (lift::(a), lift::(b), T::from_f64(t)); + let d = sub3(&b, &a); + let x = std::array::from_fn(|k| coord(&a, k) + t.clone() * coord(&d, k)); + Hom { x, w: None } +} + /// Decides the sign of `W` exactly: filter first, integers if needed. fn exact_w_sign(filter: &Hom, exact: impl FnOnce() -> Hom) -> Sign { w_sign(&filter.w).or_else(|| w_sign(&exact().w)).unwrap_or(Sign::Zero) @@ -172,12 +188,61 @@ impl Tpi { } } +/// The point `a + t (b - a)` of segment `a b`, `0 < t < 1`: an LNC point (LiNear Combination, +/// Diazzi et al. 2023, §4.3), the Steiner points of constrained Delaunay boundary recovery. It +/// lies exactly on the segment, and points on one segment are ordered by `t`. +#[derive(Debug, Clone)] +pub struct Lnc { + a: [f64; 3], + b: [f64; 3], + t: f64, + filter: [Interval; 3], +} + +impl Lnc { + /// Fails unless `a ≠ b` and `t` is an admissible double strictly between 0 and 1. + pub fn new(a: &Point3, b: &Point3, t: f64) -> Result { + if !(t > 0.0 && t < 1.0 && super::is_admissible(t)) { + return Err(Error::Invalid(format!("segment parameter {t:e} is not an admissible number strictly between 0 and 1"))); + } + if a == b { + return Err(Error::Invalid("a point on a segment needs two distinct endpoints".into())); + } + let filter = lnc_hom::(&a.0, &b.0, t).x; + Ok(Self { a: a.0, b: b.0, t, filter }) + } + + /// The segment parameter `t`. + pub fn t(&self) -> f64 { + self.t + } + + /// The segment's endpoints `a` (at `t = 0`) and `b` (at `t = 1`). + pub fn endpoints(&self) -> [[f64; 3]; 2] { + [self.a, self.b] + } + + /// The coordinates rounded to doubles, within about one unit in the last place of each + /// (error-free transforms: `b - a` and `t (b - a)` are kept as two-term sums). + pub fn rounded(&self) -> [f64; 3] { + std::array::from_fn(|k| { + let (a, b) = (self.a[k], self.b[k]); + let (d, d_lo) = two_diff(b, a); + let (p, p_lo) = two_product(self.t, d); + let tail = p_lo + self.t * d_lo; + let (s, s_lo) = two_sum(a, p); + s + (s_lo + tail) + }) + } +} + /// An explicit point or an implicit one, for the `generic_*` predicates. #[derive(Debug, Clone)] pub enum GenericPoint { Explicit(Point3), Lpi(Lpi), Tpi(Tpi), + Lnc(Lnc), } impl From for GenericPoint { @@ -187,7 +252,8 @@ impl From for GenericPoint { } impl GenericPoint { - fn explicit(&self) -> Option<&Point3> { + /// The point, when it is explicit. + pub fn explicit(&self) -> Option<&Point3> { match self { GenericPoint::Explicit(p) => Some(p), _ => None, @@ -199,6 +265,7 @@ impl GenericPoint { GenericPoint::Explicit(p) => Hom { x: lift(&p.0), w: None }, GenericPoint::Lpi(l) => l.filter.clone(), GenericPoint::Tpi(t) => t.filter.clone(), + GenericPoint::Lnc(l) => Hom { x: l.filter, w: None }, } } @@ -207,6 +274,7 @@ impl GenericPoint { GenericPoint::Explicit(p) => Hom { x: lift(&p.0), w: None }, GenericPoint::Lpi(l) => lpi_hom(&l.points), GenericPoint::Tpi(t) => tpi_hom(&t.planes), + GenericPoint::Lnc(l) => lnc_hom(&l.a, &l.b, l.t), } } @@ -215,6 +283,7 @@ impl GenericPoint { pub fn approx(&self) -> [f64; 3] { match self { GenericPoint::Explicit(p) => p.0, + GenericPoint::Lnc(l) => l.rounded(), _ => { let h = self.filter_hom(); let w = h.w.map_or(1.0, |w| w.mid()); @@ -222,6 +291,30 @@ impl GenericPoint { } } } + + /// A box that certainly contains the point: `(lo, hi)` per axis. Unbounded (`±∞`) when the + /// filter cannot bound the division of an LPI or TPI. + pub fn bounds(&self) -> ([f64; 3], [f64; 3]) { + match self { + GenericPoint::Explicit(p) => (p.0, p.0), + GenericPoint::Lnc(l) => (l.filter.map(|x| x.lo()), l.filter.map(|x| x.hi())), + _ => { + let h = self.filter_hom(); + let mut lo = [f64::NEG_INFINITY; 3]; + let mut hi = [f64::INFINITY; 3]; + if let Some(w) = h.w.filter(|w| w.lo() > 0.0 || w.hi() < 0.0) { + for (k, x) in h.x.iter().enumerate() { + let q = [x.lo() / w.lo(), x.lo() / w.hi(), x.hi() / w.lo(), x.hi() / w.hi()]; + if q.iter().all(|v| v.is_finite()) { + lo[k] = q.iter().copied().fold(f64::INFINITY, f64::min).next_down(); + hi[k] = q.iter().copied().fold(f64::NEG_INFINITY, f64::max).next_up(); + } + } + } + (lo, hi) + } + } + } } /// Runs a formula through the interval filter, then exactly. @@ -264,6 +357,41 @@ pub fn generic_orient3d(a: &GenericPoint, b: &GenericPoint, c: &GenericPoint, d: if flipped { -sign } else { sign } } +/// `det[p_i - e, |p_i - e|²]` (rows `a, b, c, d`), Shewchuk's in-sphere determinant, in +/// homogeneous form: `R_i = W_e X_i - W_i X_e = W_i W_e (p_i - e)`, and row `i` scaled by the +/// positive `(W_i W_e)²` becomes `[W_i W_e R_i, |R_i|²]`, so the sign needs no division. +fn insphere_formula([a, b, c, d, e]: [Hom; 5]) -> Option { + let row = |p: &Hom| -> ([T; 3], T) { + let r: [T; 3] = std::array::from_fn(|k| mulw(coord(&p.x, k), &e.w) - mulw(coord(&e.x, k), &p.w)); + let [x, y, z] = &r; + let lift = sq3(x, y, z); + let scale = match (&p.w, &e.w) { + (None, None) => None, + (Some(w), None) | (None, Some(w)) => Some(w.clone()), + (Some(w), Some(v)) => Some(w.clone() * v.clone()), + }; + (r.map(|v| mulw(v, &scale)), lift) + }; + let [(ra, la), (rb, lb), (rc, lc), (rd, ld)] = [row(&a), row(&b), row(&c), row(&d)]; + // Expansion along the lift column: -la·M_a + lb·M_b - lc·M_c + ld·M_d. + let m = |p: &[T; 3], q: &[T; 3], r: &[T; 3]| det3(&[p.clone(), q.clone(), r.clone()]); + let det = ld * m(&ra, &rb, &rc) - lc * m(&ra, &rb, &rd) + lb * m(&ra, &rc, &rd) - la * m(&rb, &rc, &rd); + det.sign() +} + +/// [`insphere`] on explicit or implicit points, same sign convention (positive when `e` is +/// inside the sphere through `a, b, c, d`, for `orient3d(a, b, c, d) > 0`). Exact. +pub fn generic_insphere(a: &GenericPoint, b: &GenericPoint, c: &GenericPoint, d: &GenericPoint, e: &GenericPoint) -> Sign { + if let (Some(a), Some(b), Some(c), Some(d), Some(e)) = (a.explicit(), b.explicit(), c.explicit(), d.explicit(), e.explicit()) { + return insphere(a, b, c, d, e); + } + // The 5×5 form with rows `[p, |p|², 1]` is alternating in its arguments, so moving an + // explicit point to the pivot only flips the sign. + let (points, flipped) = pivot_last([a, b, c, d, e]); + let sign = decide(points, insphere_formula, insphere_formula); + if flipped { -sign } else { sign } +} + /// `det[a' - c'; b' - c']` of the projections: `det[R_a; R_b] / (W_a W_b W_c²)`. fn orient2d_formula([a, b, c]: [Hom; 3], (u, v): (usize, usize)) -> Option { let row = |p: &Hom, k: usize| mulw(coord(&p.x, k), &c.w) - mulw(coord(&c.x, k), &p.w); diff --git a/crates/mesh/src/predicates/interval.rs b/crates/mesh/src/predicates/interval.rs index 665f575..48fdd3d 100644 --- a/crates/mesh/src/predicates/interval.rs +++ b/crates/mesh/src/predicates/interval.rs @@ -27,6 +27,16 @@ impl Interval { if lo.is_nan() || hi.is_nan() { Self::NAN } else { Self { lo: lo.next_down(), hi: hi.next_up() } } } + /// The lower bound. + pub(crate) fn lo(&self) -> f64 { + self.lo + } + + /// The upper bound. + pub(crate) fn hi(&self) -> f64 { + self.hi + } + /// The midpoint, an approximation of the value. pub(crate) fn mid(&self) -> f64 { 0.5 * self.lo + 0.5 * self.hi diff --git a/crates/mesh/src/predicates/mod.rs b/crates/mesh/src/predicates/mod.rs index c38c8d1..1a8275f 100644 --- a/crates/mesh/src/predicates/mod.rs +++ b/crates/mesh/src/predicates/mod.rs @@ -13,6 +13,11 @@ //! - M. Attene, *Indirect Predicates for Geometric Constructions*, CAD 126:102856, 2020 (paper //! only; his code is LGPL and is not read): [`Lpi`], [`Tpi`], [`GenericPoint`] and the //! `generic_*` predicates. +//! - L. Diazzi, D. Panozzo, A. Vaxman, M. Attene, *Constrained Delaunay Tetrahedrization: A +//! Robust and Practical Approach*, ACM TOG 42(6):181, 2023, §4.3 (paper only; the authors' +//! code is GPL/LGPL and is not read): the [`Lnc`] point on a segment, [`generic_insphere`], +//! and its perturbed form [`generic_insphere_sos`] with the same lexicographic lift as +//! [`insphere_sos`]. //! //! # Arithmetic assumptions //! IEEE 754 binary64, round-to-nearest-even, no extended intermediate precision, no FMA @@ -40,8 +45,8 @@ mod ring; mod sos; pub use adaptive::{incircle, incircle_filter, insphere, insphere_filter, orient2d, orient2d_filter, orient3d, orient3d_filter}; -pub use indirect::{Axis, GenericPoint, Lpi, Tpi, generic_compare, generic_orient2d, generic_orient3d, orient2d_projected}; -pub use sos::{incircle_sos, insphere_sos}; +pub use indirect::{Axis, GenericPoint, Lnc, Lpi, Tpi, generic_compare, generic_insphere, generic_orient2d, generic_orient3d, orient2d_projected}; +pub use sos::{generic_insphere_sos, generic_lex, incircle_sos, insphere_sos}; use crate::{Error, Result}; diff --git a/crates/mesh/src/predicates/sos.rs b/crates/mesh/src/predicates/sos.rs index 91d6923..69e9ced 100644 --- a/crates/mesh/src/predicates/sos.rs +++ b/crates/mesh/src/predicates/sos.rs @@ -20,6 +20,7 @@ use std::cmp::Ordering; +use super::indirect::{Axis, GenericPoint, generic_compare, generic_insphere, generic_orient3d}; use super::{Point2, Point3, Sign, incircle, insphere, orient2d, orient3d}; /// Lexicographic order on coordinates (total: `f64::total_cmp`). @@ -50,7 +51,14 @@ pub fn insphere_sos(a: &Point3, b: &Point3, c: &Point3, d: &Point3, e: &Point3) #[cold] fn insphere_perturbed(a: &Point3, b: &Point3, c: &Point3, d: &Point3, e: &Point3) -> Sign { - for i in ranked(&[a.0, b.0, c.0, d.0, e.0]) { + perturbed_cofactors([a, b, c, d, e], ranked(&[a.0, b.0, c.0, d.0, e.0]), orient3d) +} + +/// The sign of the first nonzero lift cofactor, in rank order (`order`: argument positions by +/// decreasing rank). The single statement of the perturbation rule for every point type. +fn perturbed_cofactors

(points: [&P; 5], order: [usize; 5], orient: impl Fn(&P, &P, &P, &P) -> Sign) -> Sign { + let [a, b, c, d, e] = points; + for i in order { let [p, q, r, s] = match i { 0 => [b, c, d, e], 1 => [a, c, d, e], @@ -58,7 +66,7 @@ fn insphere_perturbed(a: &Point3, b: &Point3, c: &Point3, d: &Point3, e: &Point3 3 => [a, b, c, e], _ => [a, b, c, d], }; - let minor = orient3d(p, q, r, s); + let minor = orient(p, q, r, s); if !minor.is_zero() { // Row i + 1 (1-based), lift column 4: cofactor sign (-1)^(i + 5). return if i % 2 == 1 { minor } else { -minor }; @@ -67,6 +75,46 @@ fn insphere_perturbed(a: &Point3, b: &Point3, c: &Point3, d: &Point3, e: &Point3 Sign::Zero } +/// The lexicographic order of two points' exact coordinates (`x`, then `y`, then `z`): the rank +/// order of the perturbation, for explicit and implicit points alike. On explicit points it is +/// the order [`insphere_sos`] uses (callers admit `+0.0`, never `-0.0`). +pub fn generic_lex(p: &GenericPoint, q: &GenericPoint) -> Ordering { + if let (Some(p), Some(q)) = (p.explicit(), q.explicit()) { + return lex(&p.coords(), &q.coords()); + } + Axis::ALL.into_iter().map(|a| generic_compare(p, q, a)).find(|o| o.is_ne()).unwrap_or(Ordering::Equal) +} + +/// [`insphere_sos`] on explicit or implicit points: the same lift by lexicographic rank, so a +/// triangulation built with [`insphere_sos`] stays consistent when implicit points are added. +/// Never [`Sign::Zero`] when `a, b, c, d` is not coplanar; points must be pairwise distinct. +pub fn generic_insphere_sos(a: &GenericPoint, b: &GenericPoint, c: &GenericPoint, d: &GenericPoint, e: &GenericPoint) -> Sign { + if let (Some(a), Some(b), Some(c), Some(d), Some(e)) = (a.explicit(), b.explicit(), c.explicit(), d.explicit(), e.explicit()) { + return insphere_sos(a, b, c, d, e); + } + let exact = generic_insphere(a, b, c, d, e); + if !exact.is_zero() { + return exact; + } + let points = [a, b, c, d, e]; + // Positions by decreasing rank, ties by position: an insertion sort, which needs no total + // order guarantee from the comparator to stay in bounds. + let mut order = [0usize, 1, 2, 3, 4]; + for i in 1..order.len() { + let mut j = i; + while j > 0 { + let (Some(&x), Some(&y)) = (order.get(j - 1), order.get(j)) else { break }; + let (Some(px), Some(py)) = (points.get(x), points.get(y)) else { break }; + if generic_lex(py, px) != Ordering::Greater { + break; + } + order.swap(j - 1, j); + j -= 1; + } + } + perturbed_cofactors(points, order, generic_orient3d) +} + /// [`incircle`] with symbolic perturbation: never [`Sign::Zero`] when `a, b, c` is not /// collinear. Same conventions as [`insphere_sos`]. #[inline] diff --git a/crates/mesh/src/predicates/tests.rs b/crates/mesh/src/predicates/tests.rs index c9afd53..697aafc 100644 --- a/crates/mesh/src/predicates/tests.rs +++ b/crates/mesh/src/predicates/tests.rs @@ -605,3 +605,138 @@ fn undefined_constructions_are_rejected() { assert_eq!(generic_compare(&t, &GenericPoint::Explicit(o), axis), Ordering::Equal); } } + +// --------------------------------------------------------------------------------------------- +// LNC points and the generic in-sphere (Diazzi et al. 2023, §4.3) +// --------------------------------------------------------------------------------------------- + +fn r_lnc(a: &Point3, b: &Point3, t: f64) -> RPoint { + let (a, b) = (rp(a), rp(b)); + std::array::from_fn(|k| a[k].clone() + q(t) * (b[k].clone() - a[k].clone())) +} + +fn r_insphere(p: &[RPoint; 5]) -> Sign { + let rows: Vec> = (0..4).map(|i| lifted(rsub(&p[i], &p[4]).to_vec())).collect(); + rsign(&rdet(&rows)) +} + +#[test] +fn lnc_points_lie_on_their_segment_in_parameter_order() { + let mut r = Rng(11); + for _ in 0..200 { + let (a, b) = (rand_pt(&mut r, 10.0), rand_pt(&mut r, 10.0)); + let (t0, t1) = (r.unit().clamp(1e-9, 0.5), 0.5 + r.unit().clamp(1e-9, 0.49)); + let (p0, p1) = (GenericPoint::Lnc(Lnc::new(&a, &b, t0).unwrap()), GenericPoint::Lnc(Lnc::new(&a, &b, t1).unwrap())); + let (ga, gb) = (GenericPoint::Explicit(a), GenericPoint::Explicit(b)); + for axis in [Axis::X, Axis::Y, Axis::Z] { + assert_eq!(generic_orient2d(&ga, &gb, &p0, axis), Sign::Zero); + assert_eq!(generic_orient2d(&p0, &p1, &gb, axis), Sign::Zero); + let k = axis.index(); + let want = r_lnc(&a, &b, t0)[k].cmp(&r_lnc(&a, &b, t1)[k]); + assert_eq!(generic_compare(&p0, &p1, axis), want); + } + let exact = r_lnc(&a, &b, t0); + let rounded = p0.approx(); + let (lo, hi) = p0.bounds(); + for k in 0..3 { + let x = exact[k].to_f64().unwrap(); + assert!((rounded[k] - x).abs() <= 2.0 * f64::EPSILON * x.abs().max(1e-300), "{} vs {x}", rounded[k]); + assert!(q(lo[k]) <= exact[k] && exact[k] <= q(hi[k])); + } + } + let (a, b) = (p3(0.0, 0.0, 0.0), p3(1.0, 1.0, 1.0)); + for bad in [0.0, 1.0, -0.5, 1.5, f64::NAN, f64::INFINITY, 1e-300] { + assert!(Lnc::new(&a, &b, bad).is_err(), "{bad}"); + } + assert!(Lnc::new(&a, &a, 0.5).is_err()); +} + +#[test] +fn generic_orient3d_and_insphere_match_rationals_on_lnc_points() { + let mut r = Rng(12); + let mut checked = 0; + for i in 0..120 { + let pts: Vec = (0..9).map(|_| rand_pt(&mut r, 1.0)).collect(); + let t = r.unit().clamp(1e-6, 1.0 - 1e-6); + let lnc = Lnc::new(&pts[0], &pts[1], t).unwrap(); + let Ok(lpi) = Lpi::new(&pts[2], &pts[3], &pts[4], &pts[5], &pts[6]) else { continue }; + let g = [GenericPoint::Lnc(lnc), GenericPoint::Lpi(lpi)]; + let rg = [r_lnc(&pts[0], &pts[1], t), r_lpi(&pts[2], &pts[3], &pts[4], &pts[5], &pts[6])]; + let ex: Vec = (0..5).map(|_| rand_pt(&mut r, 1.0)).collect(); + for mask in 0u32..243 { + if !(i + mask as usize).is_multiple_of(5) { + continue; + } + let (mut gp, mut rpv) = (Vec::new(), Vec::new()); + let mut m = mask; + for e in &ex { + match m % 3 { + 0 => { + gp.push(GenericPoint::Explicit(*e)); + rpv.push(rp(e)); + } + k => { + gp.push(g[k as usize - 1].clone()); + rpv.push(rg[k as usize - 1].clone()); + } + } + m /= 3; + } + let want = r_orient3d(&[rpv[0].clone(), rpv[1].clone(), rpv[2].clone(), rpv[3].clone()]); + assert_eq!(generic_orient3d(&gp[0], &gp[1], &gp[2], &gp[3]), want, "orient mask {mask}"); + let want = r_insphere(&[rpv[0].clone(), rpv[1].clone(), rpv[2].clone(), rpv[3].clone(), rpv[4].clone()]); + assert_eq!(generic_insphere(&gp[0], &gp[1], &gp[2], &gp[3], &gp[4]), want, "insphere mask {mask}"); + checked += 1; + } + } + assert!(checked > 3000, "{checked}"); + // The explicit fast path and the formula agree on the sign convention. + let (a, b, c, d) = (p3(0.0, 0.0, 0.0), p3(1.0, 0.0, 0.0), p3(0.0, 1.0, 0.0), p3(0.0, 0.0, 1.0)); + let inside = GenericPoint::Lnc(Lnc::new(&a, &p3(1.0, 1.0, 1.0), 0.25).unwrap()); + let e = |p: &Point3| GenericPoint::Explicit(*p); + let shewchuk_positive = orient3d(&b, &a, &c, &d) == Sign::Positive; + let (x, y) = if shewchuk_positive { (b, a) } else { (a, b) }; + assert_eq!(insphere(&x, &y, &c, &d, &p3(0.25, 0.25, 0.25)), Sign::Positive); + assert_eq!(generic_insphere(&e(&x), &e(&y), &e(&c), &e(&d), &inside), Sign::Positive); +} + +/// On a grid, a segment's midpoint (an LNC point) is a grid point: the perturbed in-sphere of +/// the implicit form must give the explicit form's answer, including every cospherical tie. +#[test] +fn perturbed_generic_insphere_agrees_with_explicit_points() { + let mut r = Rng(13); + let mut ties = 0; + for _ in 0..800 { + let v: Vec<[f64; 3]> = (0..5).map(|_| [r.int(1) as f64, r.int(1) as f64, r.int(1) as f64]).collect(); + if (0..5).any(|i| (0..i).any(|j| v[i] == v[j])) { + continue; + } + let p: Vec = v.iter().map(|c| p3(c[0], c[1], c[2])).collect(); + if orient3d(&p[0], &p[1], &p[2], &p[3]).is_zero() { + continue; + } + // Each point as the midpoint of a segment through it (direction from the generator). + let g: Vec = p + .iter() + .map(|c| { + if r.int(1) == 0 { + return GenericPoint::Explicit(*c); + } + let d = [r.int(2) as f64, r.int(2) as f64, (r.int(2) as f64) + 7.0]; + let x = c.coords(); + let a = p3(x[0] - d[0], x[1] - d[1], x[2] - d[2]); + let b = p3(x[0] + d[0], x[1] + d[1], x[2] + d[2]); + GenericPoint::Lnc(Lnc::new(&a, &b, 0.5).unwrap()) + }) + .collect(); + let want = insphere_sos(&p[0], &p[1], &p[2], &p[3], &p[4]); + assert_eq!(generic_insphere_sos(&g[0], &g[1], &g[2], &g[3], &g[4]), want, "{v:?}"); + ties += insphere(&p[0], &p[1], &p[2], &p[3], &p[4]).is_zero() as usize; + for i in 0..5 { + for j in 0..5 { + assert_eq!(generic_lex(&g[i], &g[j]), v[i].partial_cmp(&v[j]).unwrap()); + } + } + } + assert!(ties > 20, "{ties}"); +} diff --git a/crates/testkit/src/runner.rs b/crates/testkit/src/runner.rs index 9a62684..8c8f1f2 100644 --- a/crates/testkit/src/runner.rs +++ b/crates/testkit/src/runner.rs @@ -14,6 +14,7 @@ use std::time::Instant; use crate::corpus::{self, CATEGORIES, Part}; use modelcraft_geometry::{Design, Measure, SurfaceMesh}; +use modelcraft_mesh::cdt::{CdtOptions, tetrahedralize}; use modelcraft_mesh::surface::{BoundarySurface, SurfaceOptions}; pub const FILE: &str = "docs/meshing/corpus.md"; @@ -41,7 +42,7 @@ pub trait Stage { /// The stages, in column order. Meshers add theirs here. fn stages() -> Vec> { - vec![Box::new(SurfaceVolume), Box::new(BoundaryCheck)] + vec![Box::new(SurfaceVolume), Box::new(BoundaryCheck), Box::new(TetCdt)] } /// The boundary surface's enclosed volume against SolidCraft's measured volume and the closed @@ -86,6 +87,33 @@ impl Stage for BoundaryCheck { } } +/// The constrained Delaunay tetrahedralization (`modelcraft_mesh::cdt`) of each body's boundary: +/// tetrahedra, Steiner points on boundary edges, cavities refilled by gift-wrapping, and the +/// validity gate on the rounded mesh (every tetrahedron positive, volume equal to the boundary's). +struct TetCdt; + +impl Stage for TetCdt { + fn headers(&self) -> &'static [&'static str] { + &["CDT tets", "Steiner", "cavities", "tet gate", "CDT ms"] + } + + fn run(&self, l: &Loaded) -> Result, String> { + let start = Instant::now(); + let (mut tets, mut steiner, mut cavities) = (0, 0, 0); + for (body, mesh) in l.design.bodies.iter().zip(&l.surfaces) { + let options = SurfaceOptions { weld_distance: mesh.weld_distance, ..SurfaceOptions::default() }; + let s = BoundarySurface::new(&mesh.positions, &mesh.triangles, &mesh.tri_face, &options).map_err(|e| format!("{}: {e}", body.name))?; + let m = tetrahedralize(&s, &CdtOptions::default()).map_err(|e| format!("{}: {e}", body.name))?; + m.gate(s.volume()).and_then(|r| r.into_result()).map_err(|e| format!("{}: rounded mesh: {e}", body.name))?; + tets += m.tets.len(); + steiner += m.steiner.len(); + cavities += m.stats.cavities; + } + let ms = start.elapsed().as_secs_f64() * 1e3; + Ok(vec![tets.to_string(), steiner.to_string(), cavities.to_string(), "pass".into(), format!("{ms:.0}")]) + } +} + fn load(part: &Part) -> Result<(Loaded<'_>, f64), String> { let start = Instant::now(); let design = corpus::load(part).map_err(|e| e.to_string())?; @@ -224,4 +252,20 @@ mod tests { assert_eq!(cells[1], genus, "{name}: {cells:?}"); } } + + /// The constrained Delaunay tetrahedralization passes the gate on parts that exercise each + /// of its stages: an assembly whose rounding needs a few-ulp repair, cavities refilled by + /// gift-wrapping, and a holed plate. + #[test] + fn cdt_meshes_corpus_parts() { + let mut parts = corpus::parts().unwrap(); + parts.retain(|p| ["bolted-assembly", "thin-wall-enclosure", "analytic-plate-hole"].contains(&p.name.as_str())); + assert_eq!(parts.len(), 3); + for part in &parts { + let (loaded, _) = load(part).unwrap(); + let cells = TetCdt.run(&loaded).unwrap_or_else(|e| panic!("{}: {e}", part.name)); + assert_eq!(cells[3], "pass", "{}: {cells:?}", part.name); + assert!(cells[0].parse::().unwrap() > 0); + } + } } diff --git a/crates/ui-egui/src/control.rs b/crates/ui-egui/src/control.rs index 3b1f079..780897d 100644 --- a/crates/ui-egui/src/control.rs +++ b/crates/ui-egui/src/control.rs @@ -404,11 +404,11 @@ mod tests { assert_eq!(r["result"]["tab"], "mesh", "{r}"); let entries = r["result"]["entries"].as_array().unwrap(); assert_eq!(entries.len(), modelcraft_engine::catalogue::CATALOGUE.len()); - let tet = entries.iter().find(|e| e["id"] == "mesh.tet.generate").unwrap(); - assert_eq!((tet["coming_soon"].as_str(), tet["command"].is_null()), (Some("M2"), true)); + let soon = entries.iter().find(|e| e["id"] == "mesh.order.quadratic").unwrap(); + assert_eq!((soon["coming_soon"].as_str(), soon["command"].is_null()), (Some("M2"), true)); assert_eq!(call(&mut app, "ui.set", json!({"tab": "nope"}))["ok"], false); - let r = call(&mut app, "ui.start", json!({"command": "mesh.tet.generate"})); - assert_eq!((r["ok"].as_bool(), r["error"].as_str()), (Some(false), Some("Tet Mesh is coming soon (M2)"))); + let r = call(&mut app, "ui.start", json!({"command": "mesh.order.quadratic"})); + assert_eq!((r["ok"].as_bool(), r["error"].as_str()), (Some(false), Some("Quadratic is coming soon (M2)"))); } #[test] diff --git a/crates/ui-egui/src/lib.rs b/crates/ui-egui/src/lib.rs index c42e30b..3d814d1 100644 --- a/crates/ui-egui/src/lib.rs +++ b/crates/ui-egui/src/lib.rs @@ -520,8 +520,8 @@ mod tests { assert_eq!(app.dialog.as_ref().map(|d| d.command), Some("model.material.create")); app.dialog = None; let journal = app.session.journal.len(); - let err = app.start("mesh.tet.generate").unwrap_err(); - assert_eq!(err, "Tet Mesh is coming soon (M2)"); + let err = app.start("mesh.order.quadratic").unwrap_err(); + assert_eq!(err, "Quadratic is coming soon (M2)"); assert!(app.dialog.is_none()); assert_eq!(app.session.journal.len(), journal); assert_eq!(app.status.as_ref().map(|s| s.error), Some(false)); diff --git a/docs/architecture.md b/docs/architecture.md index 6512e09..c61ef75 100644 --- a/docs/architecture.md +++ b/docs/architecture.md @@ -61,8 +61,9 @@ Rules beyond the layer order: position. 3. **Tetrahedral meshing.** Bowyer–Watson Delaunay on Shewchuk's adaptive predicates (symbolic perturbation for degeneracies), boundary recovery by the robust constrained Delaunay method of - Diazzi et al. 2023 (implicit Steiner points on indirect predicates, gift-wrapping fallback), so - every boundary triangle and its B-rep face tag survives. Then constrained Delaunay refinement to a + Diazzi et al. 2023 (implicit Steiner points on indirect predicates, gift-wrapping facet + recovery; `modelcraft-mesh::cdt`, command `mesh.tet.generate`), so every boundary triangle and + its B-rep face tag survives. Then constrained Delaunay refinement to a size field, ODT smoothing interleaved as in Tournois et al. 2009, and the Klingner–Shewchuk 2007 improvement schedule against slivers. 4. **Hexahedral meshing (the core).** Octree/grid dual all-hex first, the one family that meshes diff --git a/docs/features.md b/docs/features.md index 7b76de9..2e67add 100644 --- a/docs/features.md +++ b/docs/features.md @@ -27,7 +27,7 @@ shown greyed with the roadmap milestone that brings them (ROADMAP.md). | Size | Global Size — One element size for the whole part | coming soon (M2) | | | Size | Seed Edges — Element counts, sizes and bias per edge and face | coming soon (M2) | | | Size | Curvature — Refine by curvature and proximity | coming soon (M2) | | -| Tet | Tet Mesh — Delaunay tetrahedral mesh of a body, boundary recovered and refined | coming soon (M2) | | +| Tet | Tet Mesh — Constrained Delaunay tet mesh of a body: every boundary triangle kept, faces tagged (refinement comes next) | ready | `mesh.tet.generate` | | Tet | Quadratic — Second-order elements with mid-side nodes on the geometry | coming soon (M2) | | | Tet | Delaunay Points — Delaunay tetrahedralization of a point list (debugging aid) | ready | `mesh.delaunay.points` | | Structured Hex | Quad Faces — All-quad meshes of faces on shared edge nodes | coming soon (M5) | | diff --git a/docs/meshing/corpus.md b/docs/meshing/corpus.md index 57e6bd7..4a09133 100644 --- a/docs/meshing/corpus.md +++ b/docs/meshing/corpus.md @@ -8,29 +8,29 @@ clock on the machine that ran it; every other column is deterministic. ## Results -| part | categories | bodies | faces | triangles | tol. mm | build ms | enclosed mm³ | SolidCraft mm³ | rel. diff | exact mm³ | boundary | genus | feature edges | check ms | -|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---| -| bent-tube | sweepable, thin-wall, genus | 1 | 18 | 4296 | 0.121 | 432 | 4213.1 | 4159.4 | 1.3e-2 | 4159.8 | ok, 2148 vertices | 1 | 268 | 11.4 | -| bolted-assembly | assembly, blocky | 3 | 22 | 284 | 0.117 | 57 | 33004.7 | 32999.0 | 1.7e-4 | 32998.9 | ok, 144 vertices | 0, 1, 1 | 104 | 0.7 | -| cross-cylinders | boolean, sharp | 1 | 10 | 1395 | 0.087 | 7548 | 32197.4 | 32357.8 | -5.0e-3 | 32365.8 | error: Cross: invalid boundary surface: 28 defects; triangles 142 and 184 are oriented inconsistently across edge 0–396 at (-15.488656021315903, 4.577998041166734, 7.853092194999726); triangles 184 and 239 are oriented inconsistently across edge 0–554 at (-18.75000000000002, 5.634596901641276, 7.632255089815496); edge 103–108 at (-0.16007784063084018, 0.03984236936862052, 9.998639270596865) is shared by 3 triangles; … | | | | -| filleted-bracket | fillet, genus, sweepable | 1 | 13 | 540 | 0.088 | 388 | 32162.3 | 32145.0 | 5.4e-4 | 32144.8 | ok, 268 vertices | 2 | 135 | 1.4 | -| i-beam | sweepable, thin-wall | 1 | 14 | 44 | 0.208 | 4 | 61600.0 | 61600.0 | 0.0e0 | 61600.0 | ok, 24 vertices | 0 | 36 | 0.1 | -| loft-transition | organic | 1 | 6 | 292 | 0.075 | 54 | 56390.2 | 56450.6 | -1.1e-3 | – | ok, 148 vertices | 0 | 40 | 0.7 | -| perforated-plate | genus, thin-wall, sweepable | 1 | 54 | 1644 | 0.144 | 270 | 26822.2 | 26764.6 | 2.2e-3 | 26764.2 | ok, 776 vertices | 24 | 828 | 4.9 | -| pinhole-block | small-feature, genus | 1 | 38 | 4044 | 0.073 | 206 | 49919.4 | 49909.6 | 2.0e-4 | 49909.5 | ok, 1992 vertices | 16 | 1548 | 12.9 | -| rounded-box | fillet | 1 | 26 | 2540 | 0.054 | 149 | 23325.7 | 23340.5 | -6.3e-4 | 23340.8 | ok, 1272 vertices | 0 | 216 | 6.7 | -| shallow-groove | near-tangent, boolean | 1 | 8 | 62 | 0.075 | 11 | 47795.2 | 47782.3 | 2.7e-4 | 47781.5 | ok, 33 vertices | 0 | 24 | 0.1 | -| sharp-wedge | sharp, thin-wall, sweepable | 1 | 5 | 8 | 0.083 | 2 | 3200.0 | 3200.0 | 0.0e0 | 3200.0 | ok, 6 vertices | 0 | 9 | 0.0 | -| sphere-on-plate | near-tangent, boolean | 1 | 12 | 2826 | 0.098 | 682 | 69353.3 | 69509.0 | -2.2e-3 | 69509.7 | ok, 1415 vertices | 0 | 124 | 6.1 | -| spline-vase | organic | 1 | 11 | 3890 | 0.106 | 217 | 94891.2 | 94948.5 | -6.0e-4 | – | ok, 1947 vertices | 0 | 243 | 8.9 | -| stepped-block | blocky, sharp | 1 | 14 | 96 | 0.078 | 35 | 47000.0 | 47000.0 | 0.0e0 | 47000.0 | ok, 48 vertices | 1 | 60 | 0.4 | -| thin-ribs | thin-wall, blocky | 1 | 26 | 92 | 0.119 | 69 | 32000.0 | 32000.0 | -4.5e-16 | 32000.0 | ok, 48 vertices | 0 | 72 | 0.5 | -| thin-wall-enclosure | thin-wall, blocky | 1 | 11 | 36 | 0.108 | 9 | 15564.0 | 15564.0 | 0.0e0 | 15564.0 | ok, 20 vertices | 0 | 28 | 0.1 | -| analytic-sphere | analytic | 1 | 6 | 2784 | 0.069 | 245 | 33353.4 | 33507.9 | -4.6e-3 | 33510.3 | ok, 1394 vertices | 0 | 96 | 6.1 | -| analytic-torus | analytic, genus | 1 | 6 | 7230 | 0.109 | 334 | 37644.1 | 37893.7 | -6.6e-3 | 37899.3 | ok, 3615 vertices | 1 | 180 | 16.1 | -| analytic-cylinder | analytic, sweepable | 1 | 4 | 1004 | 0.049 | 25 | 12533.1 | 12564.7 | -2.5e-3 | 12566.4 | ok, 504 vertices | 0 | 106 | 2.7 | -| analytic-pipe | analytic, genus, sweepable | 1 | 6 | 1448 | 0.104 | 44 | 94150.9 | 94238.0 | -9.3e-4 | 94247.8 | ok, 724 vertices | 1 | 220 | 3.9 | -| analytic-plate-hole | analytic, genus, sweepable | 1 | 8 | 112 | 0.224 | 11 | 98449.1 | 98429.5 | 2.0e-4 | 98429.2 | ok, 56 vertices | 1 | 62 | 0.4 | +| part | categories | bodies | faces | triangles | tol. mm | build ms | enclosed mm³ | SolidCraft mm³ | rel. diff | exact mm³ | boundary | genus | feature edges | check ms | CDT tets | Steiner | cavities | tet gate | CDT ms | +|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---| +| bent-tube | sweepable, thin-wall, genus | 1 | 18 | 4296 | 0.121 | 435 | 4213.1 | 4159.4 | 1.3e-2 | 4159.8 | ok, 2148 vertices | 1 | 268 | 10.9 | 15550 | 2338 | 0 | pass | 550 | +| bolted-assembly | assembly, blocky | 3 | 22 | 284 | 0.117 | 54 | 33004.7 | 32999.0 | 1.7e-4 | 32998.9 | ok, 144 vertices | 0, 1, 1 | 104 | 0.7 | 690 | 181 | 0 | pass | 31 | +| cross-cylinders | boolean, sharp | 1 | 10 | 1395 | 0.087 | 7682 | 32197.4 | 32357.8 | -5.0e-3 | 32365.8 | error: Cross: invalid boundary surface: 28 defects; triangles 142 and 184 are oriented inconsistently across edge 0–396 at (-15.488656021315903, 4.577998041166734, 7.853092194999726); triangles 184 and 239 are oriented inconsistently across edge 0–554 at (-18.75000000000002, 5.634596901641276, 7.632255089815496); edge 103–108 at (-0.16007784063084018, 0.03984236936862052, 9.998639270596865) is shared by 3 triangles; … | | | | error: Cross: invalid boundary surface: 28 defects; triangles 142 and 184 are oriented inconsistently across edge 0–396 at (-15.488656021315903, 4.577998041166734, 7.853092194999726); triangles 184 and 239 are oriented inconsistently across edge 0–554 at (-18.75000000000002, 5.634596901641276, 7.632255089815496); edge 103–108 at (-0.16007784063084018, 0.03984236936862052, 9.998639270596865) is shared by 3 triangles; … | | | | | +| filleted-bracket | fillet, genus, sweepable | 1 | 13 | 540 | 0.088 | 381 | 32162.3 | 32145.0 | 5.4e-4 | 32144.8 | ok, 268 vertices | 2 | 135 | 1.4 | 5582 | 1613 | 0 | pass | 187 | +| i-beam | sweepable, thin-wall | 1 | 14 | 44 | 0.208 | 4 | 61600.0 | 61600.0 | 0.0e0 | 61600.0 | ok, 24 vertices | 0 | 36 | 0.1 | 94 | 24 | 0 | pass | 3 | +| loft-transition | organic | 1 | 6 | 292 | 0.075 | 55 | 56390.2 | 56450.6 | -1.1e-3 | – | ok, 148 vertices | 0 | 40 | 0.7 | 1329 | 256 | 1 | pass | 33 | +| perforated-plate | genus, thin-wall, sweepable | 1 | 54 | 1644 | 0.144 | 265 | 26822.2 | 26764.6 | 2.2e-3 | 26764.2 | ok, 776 vertices | 24 | 828 | 4.6 | 33611 | 9798 | 0 | pass | 1195 | +| pinhole-block | small-feature, genus | 1 | 38 | 4044 | 0.073 | 202 | 49919.4 | 49909.6 | 2.0e-4 | 49909.5 | ok, 1992 vertices | 16 | 1548 | 12.4 | 96411 | 30221 | 0 | pass | 4408 | +| rounded-box | fillet | 1 | 26 | 2540 | 0.054 | 156 | 23325.7 | 23340.5 | -6.3e-4 | 23340.8 | ok, 1272 vertices | 0 | 216 | 6.3 | 11666 | 1913 | 0 | pass | 371 | +| shallow-groove | near-tangent, boolean | 1 | 8 | 62 | 0.075 | 11 | 47795.2 | 47782.3 | 2.7e-4 | 47781.5 | ok, 33 vertices | 0 | 24 | 0.2 | 267 | 87 | 0 | pass | 11 | +| sharp-wedge | sharp, thin-wall, sweepable | 1 | 5 | 8 | 0.083 | 2 | 3200.0 | 3200.0 | 0.0e0 | 3200.0 | ok, 6 vertices | 0 | 9 | 0.0 | 11 | 3 | 0 | pass | 0 | +| sphere-on-plate | near-tangent, boolean | 1 | 12 | 2826 | 0.098 | 703 | 69353.3 | 69509.0 | -2.2e-3 | 69509.7 | ok, 1415 vertices | 0 | 124 | 5.7 | 4898 | 309 | 9 | pass | 214 | +| spline-vase | organic | 1 | 11 | 3890 | 0.106 | 230 | 94891.2 | 94948.5 | -6.0e-4 | – | ok, 1947 vertices | 0 | 243 | 8.5 | 12461 | 1512 | 0 | pass | 419 | +| stepped-block | blocky, sharp | 1 | 14 | 96 | 0.078 | 30 | 47000.0 | 47000.0 | 0.0e0 | 47000.0 | ok, 48 vertices | 1 | 60 | 0.4 | 223 | 41 | 0 | pass | 6 | +| thin-ribs | thin-wall, blocky | 1 | 26 | 92 | 0.119 | 62 | 32000.0 | 32000.0 | -4.5e-16 | 32000.0 | ok, 48 vertices | 0 | 72 | 0.5 | 171 | 33 | 0 | pass | 6 | +| thin-wall-enclosure | thin-wall, blocky | 1 | 11 | 36 | 0.108 | 8 | 15564.0 | 15564.0 | 0.0e0 | 15564.0 | ok, 20 vertices | 0 | 28 | 0.1 | 191 | 47 | 2 | pass | 5 | +| analytic-sphere | analytic | 1 | 6 | 2784 | 0.069 | 263 | 33353.4 | 33507.9 | -4.6e-3 | 33510.3 | ok, 1394 vertices | 0 | 96 | 5.8 | 4211 | 58 | 0 | pass | 258 | +| analytic-torus | analytic, genus | 1 | 6 | 7230 | 0.109 | 345 | 37644.1 | 37893.7 | -6.6e-3 | 37899.3 | ok, 3615 vertices | 1 | 180 | 15.4 | 19790 | 466 | 0 | pass | 717 | +| analytic-cylinder | analytic, sweepable | 1 | 4 | 1004 | 0.049 | 25 | 12533.1 | 12564.7 | -2.5e-3 | 12566.4 | ok, 504 vertices | 0 | 106 | 2.7 | 7437 | 1568 | 0 | pass | 259 | +| analytic-pipe | analytic, genus, sweepable | 1 | 6 | 1448 | 0.104 | 43 | 94150.9 | 94238.0 | -9.3e-4 | 94247.8 | ok, 724 vertices | 1 | 220 | 3.7 | 6238 | 1340 | 0 | pass | 264 | +| analytic-plate-hole | analytic, genus, sweepable | 1 | 8 | 112 | 0.224 | 11 | 98449.1 | 98429.5 | 2.0e-4 | 98429.2 | ok, 56 vertices | 1 | 62 | 0.4 | 1868 | 543 | 0 | pass | 53 | ## Categories diff --git a/docs/meshing/literature.md b/docs/meshing/literature.md index f110ab7..f428ced 100644 --- a/docs/meshing/literature.md +++ b/docs/meshing/literature.md @@ -238,6 +238,21 @@ Delaunay methods with Voronoi/ODT improvement, then (2) all-hex and hex-dominant - **Complexity.** High: about 8–12k lines of Rust including the predicate generator. It is the single most valuable investment in the tet pipeline. - **Code.** `MarcoAttene/CDT` is **GPL, or LGPL by build option** (README): FORBIDDEN. +- **Implemented** in `modelcraft-mesh::cdt` from the paper only, with these choices: + - Steiner points are LNC points (§4.3) on the original segment, `a + t (b − a)` with `t` a double; the in-sphere test + on them is exact (interval filter, then exact dyadic integers) and perturbed with the same lexicographic lift as + the Delaunay step, so Bowyer–Watson insertion of Steiner points stays consistent with it. + - Facet recovery uses the modified gift-wrapping (§4.4) for every cavity, not only after Si's cavity expansion + fails: it is the always-correct path, and Si's expansion is a speed-up. Condition (ii) is an exact + proper-intersection test of closed simplices (separating plane through three of the vertices, then recursion + into the contact plane), which handles the coplanar configurations of CAD tessellations; recovered facets inside + a cavity stay constraints. + - Reference points (§3.3) are searched among the tetrahedra whose boxes meet the missing subsegment's box (they + contain the tetrahedra the subsegment crosses); on the corpus this picks the same point as searching the whole + diametral ball. + - Rounding (§5.1): a Steiner point whose rounding leaves an incident tetrahedron non-positive tries the doubles + within two ulps per coordinate, nearest first; the flip-based improvement of §5.1 is not implemented. All 20 + valid corpus parts pass the validity gate after rounding. **A2.4 Diazzi and Attene 2021, convex polyhedral meshing.** **Verdict: reference.** - **Citation.** L. Diazzi, M. Attene. *Convex polyhedral meshing for robust solid modeling.* ACM TOG 40(6):259, 2021.