From 6c36f5b2ad9c0df6c4d05b4d09befbcc77ac2c8c Mon Sep 17 00:00:00 2001 From: Irakli Kourtis Date: Tue, 13 Oct 2026 16:15:54 -0400 Subject: [PATCH] Planar all-quad meshing: size fields, TFI mapped meshes, free quad mesher (M5 base) (#18) * Planar all-quad meshing: size fields, TFI mapped meshes, free quad mesher, planar quad validity - modelcraft_mesh::size: SizeField trait, division counts, graded polyline division, per-loop parity - modelcraft_mesh::quad: transfinite (Gordon-Hall) mapped meshes of four-sided regions; free all-quad meshes of polygons with holes on prescribed boundary nodes (half-resolution CDT with Ruppert refinement, Lee-Lo recombination, uniform split, smart Laplacian plus Freitag-Plassmann untangling), exact validity gate on every element - modelcraft_mesh::quality: Robinson aspect ratio and skew for quads, exact planar quad validity Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm * quad mesher: true-boundary guard for Steiner points, untangling by corner area, parity retry, recombination quality floor; I-beam and perforated plate tests Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm * docs: quad meshing papers (B7), roadmap M5 and M2, architecture pipeline step Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm * quad mesher: simplify (shared corner and validity helpers, vertex rotation helper, guard decoupled from the CDT, verify inside the CDT, non-allocating final gate) Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm --------- Co-authored-by: jyviko <2501648+jyviko@users.noreply.github.com> Co-authored-by: Claude Opus 5.5 --- ROADMAP.md | 4 +- crates/mesh/src/lib.rs | 2 + crates/mesh/src/predicates/mod.rs | 7 + crates/mesh/src/quad/cdt.rs | 775 ++++++++++++++++++++++++++++++ crates/mesh/src/quad/domain.rs | 309 ++++++++++++ crates/mesh/src/quad/free.rs | 361 ++++++++++++++ crates/mesh/src/quad/guard.rs | 121 +++++ crates/mesh/src/quad/mod.rs | 287 +++++++++++ crates/mesh/src/quad/smooth.rs | 250 ++++++++++ crates/mesh/src/quad/tests.rs | 539 +++++++++++++++++++++ crates/mesh/src/quad/tfi.rs | 108 +++++ crates/mesh/src/quality/mod.rs | 58 ++- crates/mesh/src/quality/shape.rs | 31 ++ crates/mesh/src/size.rs | 263 ++++++++++ docs/architecture.md | 9 +- docs/meshing/literature.md | 67 +++ 16 files changed, 3185 insertions(+), 6 deletions(-) create mode 100644 crates/mesh/src/quad/cdt.rs create mode 100644 crates/mesh/src/quad/domain.rs create mode 100644 crates/mesh/src/quad/free.rs create mode 100644 crates/mesh/src/quad/guard.rs create mode 100644 crates/mesh/src/quad/mod.rs create mode 100644 crates/mesh/src/quad/smooth.rs create mode 100644 crates/mesh/src/quad/tests.rs create mode 100644 crates/mesh/src/quad/tfi.rs create mode 100644 crates/mesh/src/size.rs diff --git a/ROADMAP.md b/ROADMAP.md index a3049dc..0dfdbb9 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -20,10 +20,10 @@ 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); 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); Delaunay refinement, ODT smoothing, sliver removal; tet4 and tet10 with mid-side nodes on the B-rep; mesh verify per set | in progress | | M3 First deck | Sets and surfaces bound to SolidCraft face, edge, vertex and body names, resolved after meshing; materials (isotropic and orthotropic elastic, density) with a unit-system label and consistency check; sections (solid, shell, beam profiles) assigned by set; element formulation per region mapped to each solver; BCs (fixed, displacement, symmetry); loads (force and moment on reference points, pressure, gravity); reference points with kinematic and distributing coupling; steps (linear static, frequency) with explicit per-step load and BC tables and output presets; writers for Abaqus input, Nastran SOL 101/103 and LS-DYNA implicit; the `.modelcraft` file | planned | | M4 Own solver and first results | Linear static and modal solver (tet4/10, hex8, hex20; sparse direct, Lanczos), verified against analytic solutions; deformed shape, contours, mode shapes, probes, reaction sums | planned | -| M5 Structured hex | 2D quad meshing of faces; mapped hex of six-sided regions; sweep meshing; conforming interfaces and tet–hex transition; automatic decomposition into sweepable and mappable regions, proposed as SolidCraft splits | planned | +| 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 | | M6 Automatic hex | Octree/grid-dual all-hex for any part; hex-dominant by tet recombination; polycube for blocky parts; smoothing, untangling, boundary padding, layers through thin walls; mesh editing and orphan-mesh import | planned | | M7 Assemblies and contact | Part instances from SolidCraft occurrences; tie between non-matching meshes; general and pairwise contact with automatic pair detection, friction; bolt pretension, springs, point masses; amplitudes, multiple steps; nonlinear static steps in decks | planned | | M8 Materials, explicit, shells, composites | Plasticity, hyperelasticity, damage data per solver; explicit dynamics decks (Abaqus/Explicit, LS-DYNA); midsurface shells; composite layups and material orientations, continuum shells on swept hex; heat transfer | planned | diff --git a/crates/mesh/src/lib.rs b/crates/mesh/src/lib.rs index 1682305..486bad1 100644 --- a/crates/mesh/src/lib.rs +++ b/crates/mesh/src/lib.rs @@ -8,7 +8,9 @@ #![forbid(unsafe_code)] pub mod delaunay; +pub mod quad; pub mod quality; +pub mod size; pub mod surface; mod vec3; diff --git a/crates/mesh/src/predicates/mod.rs b/crates/mesh/src/predicates/mod.rs index 7ff8fb4..c38c8d1 100644 --- a/crates/mesh/src/predicates/mod.rs +++ b/crates/mesh/src/predicates/mod.rs @@ -117,6 +117,13 @@ pub fn is_admissible(x: f64) -> bool { x == 0.0 || (x.abs() >= MIN_MAGNITUDE && x.abs() <= MAX_MAGNITUDE) } +/// `x` with a magnitude below [`MIN_MAGNITUDE`] rounded to zero (tessellation and construction +/// noise); other values unchanged. The result is admissible when `|x| <= MAX_MAGNITUDE`. +#[inline] +pub fn snap_to_domain(x: f64) -> f64 { + if x.abs() < MIN_MAGNITUDE { 0.0 } else { x } +} + fn check(c: [f64; N]) -> Result<[f64; N]> { match c.iter().find(|x| !is_admissible(**x)) { None => Ok(c), diff --git a/crates/mesh/src/quad/cdt.rs b/crates/mesh/src/quad/cdt.rs new file mode 100644 index 0000000..e042dfd --- /dev/null +++ b/crates/mesh/src/quad/cdt.rs @@ -0,0 +1,775 @@ +//! Planar constrained Delaunay triangulation with Delaunay refinement to a size field. +//! +//! - Insertion: Bowyer–Watson (Bowyer, Comput. J. 24(2), 1981; Watson, Comput. J. 24(2), 1981) +//! inside a super-triangle, with exact `orient2d`/`incircle` (Shewchuk 1997). The cavity grows +//! from the triangle containing the point and never crosses a constrained edge, which is the +//! constrained Delaunay insertion of Shewchuk, *Delaunay refinement algorithms for triangular +//! mesh generation*, CGTA 22:21–74, 2002, §2. +//! - Constraint recovery: Sloan, *A fast algorithm for generating constrained Delaunay +//! triangulations*, Computers & Structures 47(3):441–450, 1993: flip the edges crossing a +//! segment until it appears, then restore the Delaunay property with Lawson flips. +//! - Refinement: Ruppert, *A Delaunay refinement algorithm for quality 2-dimensional mesh +//! generation*, J. Algorithms 18(3):548–585, 1995 (circumcentres of triangles that are too +//! large for the size field or have a circumradius-to-shortest-edge ratio above √2), with one +//! change: boundary segments carry prescribed nodes and are never split, so a circumcentre that +//! falls outside the region, comes closer than a fraction of the local size to a vertex, or is +//! refused by the caller's guard (clearance from the true boundary) is not inserted (a +//! too-large triangle gets its centroid instead). Ruppert's encroachment rule (diametral +//! circles) would forbid any refinement next to a boundary coarser than the size field; the +//! clearances keep the packing argument for termination. +//! +//! Every loop is bounded (walk steps, cavity rounds, flips, refinement passes and the triangle +//! count), and a cavity is checked before it is retriangulated: the new triangles all have +//! positive exact orientation and no constrained edge is removed. + +use std::collections::VecDeque; + +use super::domain::Region; +use crate::predicates::{Point2, Sign, incircle, orient2d, snap_to_domain}; +use crate::size::{SizeField, checked_size}; +use crate::{Error, Result}; + +const NONE: u32 = u32::MAX; +/// Rounds of cavity shrinking before an insertion is rejected. +const MAX_CAVITY_ROUNDS: usize = 32; +/// Refinement passes over the triangles. +const MAX_REFINE_PASSES: usize = 64; + +#[derive(Debug, Clone, Copy)] +struct Tri { + /// Counterclockwise vertices. + v: [u32; 3], + /// `n[i]`: the neighbour across the edge opposite `v[i]` (from `v[i+1]` to `v[i+2]`). + n: [u32; 3], + /// Whether the edge opposite `v[i]` is a constrained (boundary) edge. + fixed: [bool; 3], + inside: bool, + alive: bool, +} + +/// The triangles of the region: counterclockwise vertex triples into `pts`. +#[derive(Debug, Clone)] +pub(crate) struct Triangulation { + pub pts: Vec<[f64; 2]>, + pub p2: Vec, + pub tris: Vec<[u32; 3]>, +} + +pub(crate) struct Cdt { + pts: Vec<[f64; 2]>, + p2: Vec, + tris: Vec, + free: Vec, + vt: Vec, + last: u32, + mark: Vec, + epoch: u32, + max_tris: usize, +} + +/// Steiner points keep this fraction of the local size from vertices (and, through the +/// caller's guard, from the boundary). +pub(crate) const CLEARANCE: f64 = 0.3; + +/// A constructed point snapped into the predicates' domain, or `None` when it is outside the +/// model's coordinate range. +fn steiner_point(c: [f64; 2]) -> Option { + let c = [snap_to_domain(c[0]), snap_to_domain(c[1])]; + if !(c[0].abs() <= modelcraft_model::MAX_COORD && c[1].abs() <= modelcraft_model::MAX_COORD) { + return None; + } + Point2::new(c[0], c[1]).ok() +} + +fn err(what: &str) -> Error { + Error::Failed(format!("triangulation: {what}")) +} + +/// The edge index `i` of `tri` that runs from `a` to `b` (opposite `v[i]`). +fn edge_index(tri: &Tri, a: u32, b: u32) -> Option { + (0..3).find(|&i| tri.v[(i + 1) % 3] == a && tri.v[(i + 2) % 3] == b) +} + +/// Every edge of `tri` keeps `p` on its left or on it (exact). +fn contains(tri: &Tri, p2: &[Point2], p: &Point2) -> bool { + (0..3).all(|i| match (p2.get(tri.v[(i + 1) % 3] as usize), p2.get(tri.v[(i + 2) % 3] as usize)) { + (Some(a), Some(b)) => orient2d(a, b, p) != Sign::Negative, + _ => false, + }) +} + +impl Cdt { + fn t(&self, t: u32) -> Result<&Tri> { + self.tris.get(t as usize).ok_or_else(|| err("triangle index out of range")) + } + fn tm(&mut self, t: u32) -> Result<&mut Tri> { + self.tris.get_mut(t as usize).ok_or_else(|| err("triangle index out of range")) + } + fn p(&self, v: u32) -> Result<&Point2> { + self.p2.get(v as usize).ok_or_else(|| err("vertex index out of range")) + } + fn xy(&self, v: u32) -> Result<[f64; 2]> { + self.pts.get(v as usize).copied().ok_or_else(|| err("vertex index out of range")) + } + + /// Triangulate `region`'s boundary points (ids preserved) and recover its loops as + /// constrained edges; triangles inside the region are marked. + pub fn build(region: &Region, max_tris: usize) -> Result { + let n = region.pts.len(); + let (lo, hi) = super::domain::bbox(region.pts.iter().copied()); + let c = [0.5 * (lo[0] + hi[0]), 0.5 * (lo[1] + hi[1])]; + let d = (hi[0] - lo[0]).max(hi[1] - lo[1]).max(1.0); + let corners = [[c[0] - 20.0 * d, c[1] - 10.0 * d], [c[0] + 20.0 * d, c[1] - 10.0 * d], [c[0], c[1] + 20.0 * d]]; + let mut pts = region.pts.clone(); + let mut p2 = region.p2.clone(); + for k in corners { + let k = [snap_to_domain(k[0]), snap_to_domain(k[1])]; + p2.push(Point2::new(k[0], k[1])?); + pts.push(k); + } + let s = n as u32; + let mut cdt = Cdt { + vt: vec![NONE; pts.len()], + pts, + p2, + tris: vec![Tri { v: [s, s + 1, s + 2], n: [NONE; 3], fixed: [false; 3], inside: false, alive: true }], + free: Vec::new(), + last: 0, + mark: vec![0], + epoch: 0, + max_tris, + }; + for v in [s, s + 1, s + 2] { + if let Some(x) = cdt.vt.get_mut(v as usize) { + *x = 0; + } + } + // Insert in Morton order (short walks), ties by id. + let span = [(hi[0] - lo[0]).max(f64::MIN_POSITIVE), (hi[1] - lo[1]).max(f64::MIN_POSITIVE)]; + let key = |p: &[f64; 2]| { + let q = |x: f64, l: f64, s: f64| ((((x - l) / s) * 65535.0) as u32).min(65535); + let (mut x, mut y) = (q(p[0], lo[0], span[0]), q(p[1], lo[1], span[1])); + let mut k = 0u64; + for bit in 0..16 { + k |= u64::from(x & 1) << (2 * bit) | u64::from(y & 1) << (2 * bit + 1); + x >>= 1; + y >>= 1; + } + k + }; + let mut order: Vec<(u64, u32)> = region.pts.iter().enumerate().map(|(i, p)| (key(p), i as u32)).collect(); + order.sort_unstable(); + for (_, v) in order { + if !cdt.insert(v, None)? { + return Err(err("a boundary point could not be inserted")); + } + } + for l in ®ion.loops { + for (i, &a) in l.iter().enumerate() { + let b = l.get((i + 1) % l.len()).copied().ok_or_else(|| err("loop"))?; + cdt.recover(a, b)?; + } + } + cdt.classify()?; + #[cfg(test)] + assert_eq!(cdt.verify(), Ok(()), "after build"); + Ok(cdt) + } + + fn alloc(&mut self, tri: Tri) -> Result { + if let Some(t) = self.free.pop() { + *self.tm(t)? = tri; + return Ok(t); + } + if self.tris.len() >= self.max_tris { + return Err(Error::Limit(format!("the triangulation exceeds {} triangles", self.max_tris))); + } + self.tris.push(tri); + self.mark.push(0); + Ok(self.tris.len() as u32 - 1) + } + + /// The triangle containing `p` (on its boundary counts), by a visibility walk from `start` + /// with a linear scan as the fallback. + fn locate(&self, p: &Point2, start: u32) -> Result { + let mut t = if self.t(start).is_ok_and(|x| x.alive) { start } else { self.any_alive()? }; + let cap = 64 + 4 * self.tris.len(); + 'walk: for step in 0..cap { + let tri = *self.t(t)?; + for k in 0..3 { + let i = (k + step) % 3; + let (a, b) = (tri.v[(i + 1) % 3], tri.v[(i + 2) % 3]); + if orient2d(self.p(a)?, self.p(b)?, p) == Sign::Negative { + let next = tri.n[i]; + if next == NONE { + return Err(err("point outside the super-triangle")); + } + t = next; + continue 'walk; + } + } + return Ok(t); + } + self.tris.iter().position(|tri| tri.alive && contains(tri, &self.p2, p)).map(|i| i as u32).ok_or_else(|| err("point location failed")) + } + + fn any_alive(&self) -> Result { + self.tris.iter().position(|t| t.alive).map(|i| i as u32).ok_or_else(|| err("no triangles")) + } + + /// Insert vertex `v` (already in `pts`). Boundary points (`steiner == None`) must be inserted. + /// Steiner points (`Some(min_distance)`) may be rejected (outside the region, closer than + /// `min_distance` to a vertex, or no valid cavity), leaving the triangulation unchanged; + /// returns whether the point went in. + fn insert(&mut self, v: u32, steiner: Option) -> Result { + let p = *self.p(v)?; + let t0 = self.locate(&p, self.last)?; + let tri0 = *self.t(t0)?; + for &u in &tri0.v { + if self.p(u)?.coords() == p.coords() { + return match steiner { + Some(_) => Ok(false), + None => Err(err("repeated point")), + }; + } + } + if steiner.is_some() && !tri0.inside { + return Ok(false); + } + let mut excluded: Vec = Vec::new(); + for _ in 0..MAX_CAVITY_ROUNDS { + let (cavity, edges) = self.cavity(t0, &p, &excluded)?; + let mut bad: Vec = Vec::new(); + // A constrained edge inside the cavity: drop the triangle on the far side from p. + for &t in &cavity { + let tri = *self.t(t)?; + for i in 0..3 { + if tri.fixed[i] && tri.n[i] != NONE && self.marked(tri.n[i]) { + let (a, b) = (tri.v[(i + 1) % 3], tri.v[(i + 2) % 3]); + if orient2d(self.p(a)?, self.p(b)?, &p) != Sign::Positive { + bad.push(t); + } + } + } + } + // Every new triangle (a, b, p) must be counterclockwise. + for &(a, b, _, _, owner) in &edges { + if orient2d(self.p(a)?, self.p(b)?, &p) != Sign::Positive { + bad.push(owner); + } + } + if let Some(min_distance) = steiner { + let pc = p.coords(); + for &(a, _, _, _, _) in &edges { + let pa = self.xy(a)?; + if (pa[0] - pc[0]).powi(2) + (pa[1] - pc[1]).powi(2) < min_distance * min_distance { + return Ok(false); + } + } + } + if bad.is_empty() && edges.len() == cavity.len() + 2 { + self.retriangulate(v, &cavity, &edges, tri0.inside)?; + return Ok(true); + } + if bad.contains(&t0) || bad.is_empty() { + break; + } + excluded.extend(bad); + excluded.sort_unstable(); + excluded.dedup(); + } + match steiner { + Some(_) => Ok(false), + None => Err(err("no valid cavity for a boundary point")), + } + } + + fn marked(&self, t: u32) -> bool { + self.mark.get(t as usize).is_some_and(|&m| m == self.epoch) + } + + /// The Bowyer–Watson cavity of `p` from `t0` (never crossing constrained edges, skipping + /// `excluded`), marked with the current epoch, and its boundary edges `(a, b, outside + /// neighbour, constrained, owner)` with the cavity on their left. + #[allow(clippy::type_complexity)] + fn cavity(&mut self, t0: u32, p: &Point2, excluded: &[u32]) -> Result<(Vec, Vec<(u32, u32, u32, bool, u32)>)> { + self.epoch = self.epoch.wrapping_add(1); + if self.epoch == 0 { + self.mark.iter_mut().for_each(|m| *m = 0); + self.epoch = 1; + } + let epoch = self.epoch; + let set = |mark: &mut Vec, t: u32| { + if let Some(m) = mark.get_mut(t as usize) { + *m = epoch; + } + }; + set(&mut self.mark, t0); + let mut cavity = vec![t0]; + let mut k = 0; + while let Some(&t) = cavity.get(k) { + k += 1; + let tri = *self.t(t)?; + for i in 0..3 { + let nb = tri.n[i]; + if nb == NONE || tri.fixed[i] || self.marked(nb) || excluded.binary_search(&nb).is_ok() { + continue; + } + let o = *self.t(nb)?; + if incircle(self.p(o.v[0])?, self.p(o.v[1])?, self.p(o.v[2])?, p) == Sign::Positive { + set(&mut self.mark, nb); + cavity.push(nb); + } + } + } + let mut edges = Vec::new(); + for &t in &cavity { + let tri = *self.t(t)?; + for i in 0..3 { + let nb = tri.n[i]; + // Shared with another cavity triangle (a constrained one is reported by `insert`). + if nb != NONE && self.marked(nb) { + continue; + } + edges.push((tri.v[(i + 1) % 3], tri.v[(i + 2) % 3], nb, tri.fixed[i], t)); + } + } + Ok((cavity, edges)) + } + + fn retriangulate(&mut self, v: u32, cavity: &[u32], edges: &[(u32, u32, u32, bool, u32)], inside: bool) -> Result<()> { + for &t in cavity { + self.tm(t)?.alive = false; + self.free.push(t); + } + let mut made: Vec<(u32, u32, u32)> = Vec::with_capacity(edges.len()); // (a, b, tri) + for &(a, b, nb, fixed, _) in edges { + let t = self.alloc(Tri { v: [a, b, v], n: [NONE, NONE, nb], fixed: [false, false, fixed], inside, alive: true })?; + made.push((a, b, t)); + if nb != NONE { + self.relink(nb, a, b, t)?; + } + } + let mut by_a: Vec<(u32, u32)> = made.iter().map(|&(a, _, t)| (a, t)).collect(); + let mut by_b: Vec<(u32, u32)> = made.iter().map(|&(_, b, t)| (b, t)).collect(); + by_a.sort_unstable(); + by_b.sort_unstable(); + let find = |list: &[(u32, u32)], key: u32| -> Result { + let i = list.binary_search_by_key(&key, |e| e.0).map_err(|_| err("cavity boundary is not a closed loop"))?; + list.get(i).map(|e| e.1).ok_or_else(|| err("cavity lookup")) + }; + for &(a, b, t) in &made { + let next = find(&by_a, b)?; // the triangle (b, c, v): shares edge b–v + let prev = find(&by_b, a)?; // the triangle (z, a, v): shares edge v–a + let tri = self.tm(t)?; + tri.n[0] = next; + tri.n[1] = prev; + for u in [a, b, v] { + if let Some(x) = self.vt.get_mut(u as usize) { + *x = t; + } + } + self.last = t; + } + Ok(()) + } + + /// In triangle `t`, point the edge from `b` to `a` at `to`. + fn relink(&mut self, t: u32, a: u32, b: u32, to: u32) -> Result<()> { + let tri = self.tm(t)?; + let i = edge_index(tri, b, a).ok_or_else(|| err("neighbour does not share the edge"))?; + tri.n[i] = to; + Ok(()) + } + + /// Visit the triangles around vertex `v` counterclockwise, each with `v`'s local index, + /// until `visit` returns a value (or the rotation closes: `None`). + fn around(&self, v: u32, mut visit: impl FnMut(u32, &Tri, usize) -> Result>) -> Result> { + let start = self.vt.get(v as usize).copied().ok_or_else(|| err("vertex"))?; + let mut t = start; + for _ in 0..self.tris.len() + 1 { + let tri = self.t(t)?; + let k = tri.v.iter().position(|&x| x == v).ok_or_else(|| err("vertex map is stale"))?; + if let Some(r) = visit(t, tri, k)? { + return Ok(Some(r)); + } + t = tri.n[(k + 1) % 3]; + if t == NONE || t == start { + return Ok(None); + } + } + Err(err("vertex rotation did not close")) + } + + /// The triangle and local index of the directed edge `a → b` (the edge opposite `v[i]` runs + /// from `v[i+1]` to `v[i+2]`). + fn find_edge(&self, a: u32, b: u32) -> Result> { + self.around(a, |t, tri, k| Ok((tri.v[(k + 1) % 3] == b).then_some((t, (k + 2) % 3)))) + } + + /// Flip the edge opposite `v[i]` of `t`. Returns the new edge's end points. + fn flip(&mut self, t: u32, i: usize) -> Result<(u32, u32)> { + let a = *self.t(t)?; + let (p, q, r) = (a.v[i], a.v[(i + 1) % 3], a.v[(i + 2) % 3]); + let t2 = a.n[i]; + let b = *self.t(t2)?; + let j = edge_index(&b, r, q).ok_or_else(|| err("flip: neighbour mismatch"))?; + let s = b.v[j]; + let (n_pq, f_pq) = (a.n[(i + 2) % 3], a.fixed[(i + 2) % 3]); + let (n_rp, f_rp) = (a.n[(i + 1) % 3], a.fixed[(i + 1) % 3]); + let (n_sr, f_sr) = (b.n[(j + 2) % 3], b.fixed[(j + 2) % 3]); + let (n_qs, f_qs) = (b.n[(j + 1) % 3], b.fixed[(j + 1) % 3]); + *self.tm(t)? = Tri { v: [p, q, s], n: [n_qs, t2, n_pq], fixed: [f_qs, false, f_pq], inside: a.inside, alive: true }; + *self.tm(t2)? = Tri { v: [p, s, r], n: [n_sr, n_rp, t], fixed: [f_sr, f_rp, false], inside: b.inside, alive: true }; + if n_qs != NONE { + self.relink(n_qs, q, s, t)?; + } + if n_rp != NONE { + self.relink(n_rp, r, p, t2)?; + } + for (u, x) in [(p, t), (q, t), (s, t2), (r, t2)] { + if let Some(y) = self.vt.get_mut(u as usize) { + *y = x; + } + } + Ok((p, s)) + } + + /// The edge opposite `v[i]` of `t` with its far vertex, when it is flippable (strictly convex + /// quadrilateral). + fn flippable(&self, t: u32, i: usize) -> Result> { + let a = *self.t(t)?; + let t2 = a.n[i]; + if t2 == NONE || a.fixed[i] { + return Ok(None); + } + let (p, q, r) = (a.v[i], a.v[(i + 1) % 3], a.v[(i + 2) % 3]); + let b = *self.t(t2)?; + let s = b.v[edge_index(&b, r, q).ok_or_else(|| err("neighbour mismatch"))?]; + let (pp, ps) = (self.p(p)?, self.p(s)?); + Ok((orient2d(pp, ps, self.p(q)?) * orient2d(pp, ps, self.p(r)?) == Sign::Negative).then_some(s)) + } + + /// Make `a–b` an edge (Sloan 1993) and mark it constrained. + fn recover(&mut self, a: u32, b: u32) -> Result<()> { + if self.find_edge(a, b)?.is_none() { + let crossing = self.crossing_edges(a, b)?; + let (pa, pb) = (*self.p(a)?, *self.p(b)?); + let mut queue: VecDeque<(u32, u32)> = crossing.into_iter().collect(); + let cap = 1000 + 16 * queue.len() * queue.len(); + let mut fresh = Vec::new(); + let mut steps = 0usize; + while let Some((u, w)) = queue.pop_front() { + steps += 1; + if steps > cap { + return Err(Error::Limit(format!("boundary edge recovery did not finish in {cap} flips"))); + } + let (t, i) = self.find_edge(u, w)?.ok_or_else(|| err("crossing edge vanished"))?; + match self.flippable(t, i)? { + Some(_) => { + let (x, y) = self.flip(t, i)?; + let crosses = + x != a && x != b && y != a && y != b && orient2d(&pa, &pb, self.p(x)?) * orient2d(&pa, &pb, self.p(y)?) == Sign::Negative; + if crosses { queue.push_back((x, y)) } else { fresh.push((x, y)) } + } + None => queue.push_back((u, w)), + } + } + self.mark_fixed(a, b)?; + self.lawson(fresh)?; + return Ok(()); + } + self.mark_fixed(a, b) + } + + fn mark_fixed(&mut self, a: u32, b: u32) -> Result<()> { + let (t, i) = self.find_edge(a, b)?.ok_or_else(|| err("recovered edge is missing"))?; + let nb = self.t(t)?.n[i]; + self.tm(t)?.fixed[i] = true; + if nb != NONE { + let tri = self.tm(nb)?; + if let Some(k) = edge_index(tri, b, a) { + tri.fixed[k] = true; + } + } + Ok(()) + } + + /// The edges crossed by the open segment `a–b`, from `a` to `b`. + fn crossing_edges(&self, a: u32, b: u32) -> Result> { + let (pa, pb) = (*self.p(a)?, *self.p(b)?); + let found = self.around(a, |t, tri, k| { + let (u, w) = (tri.v[(k + 1) % 3], tri.v[(k + 2) % 3]); + let (ou, ow) = (orient2d(&pa, self.p(u)?, &pb), orient2d(&pa, self.p(w)?, &pb)); + if ou.is_zero() || ow.is_zero() { + let x = if ou.is_zero() { u } else { w }; + let px = self.xy(x)?; + let (qa, qb) = (pa.coords(), pb.coords()); + if (px[0] - qa[0]) * (qb[0] - qa[0]) + (px[1] - qa[1]) * (qb[1] - qa[1]) > 0.0 { + return Err(Error::Invalid("a point lies on a boundary edge".into())); + } + } + Ok((ou == Sign::Positive && ow == Sign::Negative).then_some((t, u, w))) + })?; + let (mut t, mut right, mut left) = found.ok_or_else(|| err("segment start not found"))?; + let mut out = vec![(right, left)]; + for _ in 0..self.tris.len() + 1 { + let tri = *self.t(t)?; + let i = tri.v.iter().position(|&x| x != right && x != left).ok_or_else(|| err("walk"))?; + let next = tri.n[i]; + let nt = *self.t(next)?; + let d = *nt.v.iter().find(|&&x| x != right && x != left).ok_or_else(|| err("walk"))?; + if d == b { + return Ok(out); + } + match orient2d(&pa, &pb, self.p(d)?) { + Sign::Positive => left = d, + Sign::Negative => right = d, + Sign::Zero => return Err(Error::Invalid("a point lies on a boundary edge".into())), + } + out.push((right, left)); + t = next; + } + Err(err("segment walk did not reach its end")) + } + + /// Lawson flips from the given edges until every unconstrained edge is locally Delaunay. + fn lawson(&mut self, edges: Vec<(u32, u32)>) -> Result<()> { + let mut stack = edges; + let cap = 1000 + 64 * self.tris.len(); + let mut steps = 0usize; + while let Some((u, w)) = stack.pop() { + steps += 1; + if steps > cap { + return Err(Error::Limit("Delaunay restoration did not finish".into())); + } + let Some((t, i)) = self.find_edge(u, w)? else { continue }; + let tri = *self.t(t)?; + if tri.fixed[i] || tri.n[i] == NONE { + continue; + } + let Some(s) = self.flippable(t, i)? else { continue }; + let [a, b, c] = tri.v; + if incircle(self.p(a)?, self.p(b)?, self.p(c)?, self.p(s)?) != Sign::Positive { + continue; + } + let (p, _) = self.flip(t, i)?; + let (q, r) = (tri.v[(i + 1) % 3], tri.v[(i + 2) % 3]); + stack.extend([(p, q), (q, s), (s, r), (r, p)]); + } + Ok(()) + } + + /// Mark the triangles inside the region: those separated from the super-triangle by an odd + /// number of constrained edges. + fn classify(&mut self) -> Result<()> { + let n = self.tris.len(); + let mut depth = vec![u32::MAX; n]; + let start = self.vt.get(self.pts.len().saturating_sub(1)).copied().ok_or_else(|| err("super vertex"))?; + let mut current = vec![start]; + if let Some(d) = depth.get_mut(start as usize) { + *d = 0; + } + let mut level = 0u32; + while !current.is_empty() { + let mut next = Vec::new(); + let mut k = 0; + while let Some(&t) = current.get(k) { + k += 1; + let tri = *self.t(t)?; + for i in 0..3 { + let nb = tri.n[i]; + if nb == NONE || depth.get(nb as usize).is_none_or(|&d| d != u32::MAX) { + continue; + } + if tri.fixed[i] { + next.push(nb); + } else { + if let Some(d) = depth.get_mut(nb as usize) { + *d = level; + } + current.push(nb); + } + } + } + level += 1; + next.retain(|&t| depth.get(t as usize).is_some_and(|&d| d == u32::MAX)); + next.sort_unstable(); + next.dedup(); + for &t in &next { + if let Some(d) = depth.get_mut(t as usize) { + *d = level; + } + } + current = next; + if level > 4 { + break; + } + } + for (tri, d) in self.tris.iter_mut().zip(depth) { + tri.inside = tri.alive && d == 1; + } + Ok(()) + } + + /// Try to insert Steiner point `p` (located from `hint`), at least `min_distance` from the + /// vertices of its cavity. + fn try_steiner(&mut self, p: Point2, hint: u32, min_distance: f64) -> Result { + let v = self.pts.len() as u32; + self.pts.push(p.coords()); + self.p2.push(p); + self.vt.push(NONE); + self.last = hint; + let ok = self.insert(v, Some(min_distance))?; + if !ok { + self.pts.pop(); + self.p2.pop(); + self.vt.pop(); + } + Ok(ok) + } + + /// Delaunay refinement of the inside triangles to `scale × size field`. A Steiner point `p` + /// for local size `h` is skipped when `allows(p, h)` refuses it (the caller's clearance from + /// the true boundary) or when it comes closer than [`CLEARANCE`]` × h` to a vertex. + pub fn refine(&mut self, size: &dyn SizeField, scale: f64, allows: &dyn Fn(&Point2, f64) -> bool) -> Result<()> { + for _ in 0..MAX_REFINE_PASSES { + let mut inserted = 0usize; + let n = self.tris.len(); + for t in 0..n as u32 { + let tri = *self.t(t)?; + if !tri.alive || !tri.inside { + continue; + } + let [a, b, c] = [self.xy(tri.v[0])?, self.xy(tri.v[1])?, self.xy(tri.v[2])?]; + let centroid = [(a[0] + b[0] + c[0]) / 3.0, (a[1] + b[1] + c[1]) / 3.0]; + let h = scale * checked_size(size, [centroid[0], centroid[1], 0.0])?; + let (bx, by, cx, cy) = (b[0] - a[0], b[1] - a[1], c[0] - a[0], c[1] - a[1]); + let d = 2.0 * (bx * cy - by * cx); + let (b2, c2) = (bx * bx + by * by, cx * cx + cy * cy); + let (ux, uy) = ((cy * b2 - by * c2) / d, (bx * c2 - cx * b2) / d); + let r = (ux * ux + uy * uy).sqrt(); + let shortest = b2.min(c2).min((cx - bx).powi(2) + (cy - by).powi(2)).sqrt(); + if !r.is_finite() || !(d > 0.0) { + continue; + } + let too_big = r > 0.75 * h; + let poor = r > std::f64::consts::SQRT_2 * shortest; + if !too_big && !poor { + continue; + } + let candidates = [Some([a[0] + ux, a[1] + uy]), too_big.then_some(centroid)]; + for c in candidates.into_iter().flatten() { + if let Some(p) = steiner_point(c).filter(|p| allows(p, h)) + && self.try_steiner(p, t, CLEARANCE * h)? + { + inserted += 1; + break; + } + } + } + if inserted == 0 { + break; + } + } + #[cfg(test)] + assert_eq!(self.verify(), Ok(()), "after refinement"); + Ok(()) + } + + /// The inside triangles around vertex `v`. + fn inside_fan(&self, v: u32) -> Result> { + let mut fan = Vec::new(); + self.around(v, |t, tri, _| { + if tri.inside { + fan.push(t); + } + Ok(None::<()>) + })?; + Ok(fan) + } + + /// [`Cdt::ensure_fan`] for every vertex with its count in `need`. + pub fn ensure_fans(&mut self, need: &[usize]) -> Result<()> { + for (v, &k) in need.iter().enumerate() { + if k > 1 { + self.ensure_fan(v as u32, k)?; + } + } + #[cfg(test)] + assert_eq!(self.verify(), Ok(()), "after fans"); + Ok(()) + } + + /// Give boundary vertex `v` at least `count` inside triangles by inserting the centroid of + /// its widest triangle (bounded attempts; the caller's validity gate judges the result). + fn ensure_fan(&mut self, v: u32, count: usize) -> Result<()> { + for _ in 0..count + 2 { + let fan = self.inside_fan(v)?; + if fan.len() >= count { + return Ok(()); + } + let pv = self.xy(v)?; + let mut widest: Option<(f64, u32, [f64; 2])> = None; + for &t in &fan { + let tri = *self.t(t)?; + let k = tri.v.iter().position(|&x| x == v).ok_or_else(|| err("fan"))?; + let (a, b) = (self.xy(tri.v[(k + 1) % 3])?, self.xy(tri.v[(k + 2) % 3])?); + let (u, w) = ([a[0] - pv[0], a[1] - pv[1]], [b[0] - pv[0], b[1] - pv[1]]); + let cos = (u[0] * w[0] + u[1] * w[1]) / ((u[0] * u[0] + u[1] * u[1]).sqrt() * (w[0] * w[0] + w[1] * w[1]).sqrt()); + let c = [(pv[0] + a[0] + b[0]) / 3.0, (pv[1] + a[1] + b[1]) / 3.0]; + if widest.is_none_or(|(m, _, _)| cos < m) { + widest = Some((cos, t, c)); + } + } + let Some((t, p)) = widest.and_then(|(_, t, c)| steiner_point(c).map(|p| (t, p))) else { return Ok(()) }; + if !self.try_steiner(p, t, 0.0)? { + return Ok(()); + } + } + Ok(()) + } + + /// Consistency check for tests: every live triangle counterclockwise (exact), neighbours + /// mutual, constrained flags matching on both sides. + #[cfg(test)] + pub fn verify(&self) -> std::result::Result<(), String> { + for (t, tri) in self.tris.iter().enumerate() { + if !tri.alive { + continue; + } + let [a, b, c] = tri.v; + let o = orient2d(&self.p2[a as usize], &self.p2[b as usize], &self.p2[c as usize]); + if o != Sign::Positive { + return Err(format!("triangle {t} {:?} has orientation {o:?} (inside {})", tri.v, tri.inside)); + } + for i in 0..3 { + let nb = tri.n[i]; + if nb == NONE { + continue; + } + let o = &self.tris[nb as usize]; + let (p, q) = (tri.v[(i + 1) % 3], tri.v[(i + 2) % 3]); + let back = (0..3).find(|&j| o.v[(j + 1) % 3] == q && o.v[(j + 2) % 3] == p); + match back { + Some(j) if o.n[j] == t as u32 && o.alive && o.fixed[j] == tri.fixed[i] => { + let s = o.v[j]; + if !tri.fixed[i] + && incircle(&self.p2[a as usize], &self.p2[b as usize], &self.p2[c as usize], &self.p2[s as usize]) == Sign::Positive + { + return Err(format!("triangle {t} edge {i} is not locally Delaunay (inside {})", tri.inside)); + } + } + _ => return Err(format!("triangle {t} edge {i}: neighbour {nb} does not match")), + } + } + } + Ok(()) + } + + /// The inside triangles. + pub fn into_triangulation(self) -> Triangulation { + let tris = self.tris.iter().filter(|t| t.alive && t.inside).map(|t| t.v).collect(); + Triangulation { pts: self.pts, p2: self.p2, tris } + } +} diff --git a/crates/mesh/src/quad/domain.rs b/crates/mesh/src/quad/domain.rs new file mode 100644 index 0000000..24b9123 --- /dev/null +++ b/crates/mesh/src/quad/domain.rs @@ -0,0 +1,309 @@ +//! Validation of a planar region given as boundary loops: finite, admissible coordinates, no +//! zero-length edges, no repeated points, no crossing or touching edges, holes inside the outer +//! loop and not inside each other. Every geometric decision is an exact `orient2d` or an exact +//! coordinate comparison. + +use std::cmp::Ordering; + +use crate::predicates::{Point2, Sign, orient2d, snap_to_domain}; +use crate::{Error, Result}; + +/// The most boundary points of one region. +pub const MAX_BOUNDARY_POINTS: usize = 1 << 20; +/// The most loops (outer plus holes) of one region. +pub const MAX_LOOPS: usize = 1 << 14; +/// The most segment-cell insertions and pair tests of the intersection check. +const MAX_CHECK_WORK: usize = 1 << 28; + +/// A validated region: the points in input order and the loops as point ids, the outer loop +/// first and counterclockwise, holes clockwise. +#[derive(Debug, Clone)] +pub(crate) struct Region { + pub pts: Vec<[f64; 2]>, + pub p2: Vec, + pub loops: Vec>, +} + +impl Region { + pub fn point(&self, v: u32) -> Result<&Point2> { + self.p2.get(v as usize).ok_or_else(|| Error::Failed(format!("point {v} out of range"))) + } +} + +/// The bounding box `(lo, hi)` of points (`lo > hi` when there are none). +pub(crate) fn bbox(pts: impl IntoIterator) -> ([f64; 2], [f64; 2]) { + pts.into_iter() + .fold(([f64::INFINITY; 2], [f64::NEG_INFINITY; 2]), |(lo, hi), p| ([lo[0].min(p[0]), lo[1].min(p[1])], [hi[0].max(p[0]), hi[1].max(p[1])])) +} + +/// Admissible exact-predicate point from user coordinates (finite, within the model limit; tiny +/// magnitudes rounded to zero). +pub(crate) fn admissible(c: [f64; 2]) -> Result<[f64; 2]> { + modelcraft_model::check_coordinates(&c)?; + Ok([snap_to_domain(c[0]), snap_to_domain(c[1])]) +} + +/// Validate loops given as point lists (the first is the outer boundary, the rest are holes, in +/// either orientation). Point ids number the points in input order, loop by loop. +pub(crate) fn region(input: &[Vec<[f64; 2]>]) -> Result { + if input.is_empty() { + return Err(Error::Invalid("a region needs an outer boundary loop".into())); + } + if input.len() > MAX_LOOPS { + return Err(Error::Limit(format!("{} loops exceed the limit of {MAX_LOOPS}", input.len()))); + } + let total: usize = input.iter().map(Vec::len).sum(); + if total > MAX_BOUNDARY_POINTS { + return Err(Error::Limit(format!("{total} boundary points exceed the limit of {MAX_BOUNDARY_POINTS}"))); + } + let mut pts = Vec::with_capacity(total); + let mut loops = Vec::with_capacity(input.len()); + for (k, l) in input.iter().enumerate() { + if l.len() < 3 { + return Err(Error::Invalid(format!("loop {k} has {} points; a loop needs at least 3", l.len()))); + } + let mut ids = Vec::with_capacity(l.len()); + for (i, c) in l.iter().enumerate() { + let c = admissible(*c).map_err(|e| Error::Invalid(format!("loop {k} point {i}: {e}")))?; + ids.push(pts.len() as u32); + pts.push(c); + } + loops.push(ids); + } + let p2 = pts.iter().map(|c| Point2::new(c[0], c[1])).collect::>>()?; + let mut r = Region { pts, p2, loops }; + check_distinct(&r)?; + check_simple(&r)?; + orient(&mut r)?; + check_nesting(&r)?; + Ok(r) +} + +fn lex(a: &[f64; 2], b: &[f64; 2]) -> Ordering { + a[0].total_cmp(&b[0]).then(a[1].total_cmp(&b[1])) +} + +/// Report a point by its loop and position. +fn locate(r: &Region, v: u32) -> String { + for (k, l) in r.loops.iter().enumerate() { + if let Some(i) = l.iter().position(|&x| x == v) { + return format!("loop {k} point {i}"); + } + } + format!("point {v}") +} + +/// No point appears twice (this also rejects zero-length edges). +fn check_distinct(r: &Region) -> Result<()> { + let mut order: Vec = (0..r.pts.len() as u32).collect(); + order.sort_by(|&a, &b| match (r.pts.get(a as usize), r.pts.get(b as usize)) { + (Some(p), Some(q)) => lex(p, q).then(a.cmp(&b)), + _ => a.cmp(&b), + }); + for w in order.windows(2) { + if let [a, b] = *w + && r.pts.get(a as usize) == r.pts.get(b as usize) + { + let what = if r.loops.iter().any(|l| l.windows(2).any(|e| e == [a, b] || e == [b, a]) || is_closing(l, a, b)) { + "zero-length edge" + } else { + "repeated point" + }; + return Err(Error::Invalid(format!("{what}: {} and {} coincide at {:?}", locate(r, a), locate(r, b), r.pts.get(a as usize)))); + } + } + Ok(()) +} + +fn is_closing(l: &[u32], a: u32, b: u32) -> bool { + matches!((l.first(), l.last()), (Some(&f), Some(&z)) if (f == a && z == b) || (f == b && z == a)) +} + +/// The boundary segments `(a, b)` of every loop, in loop order. +pub(crate) fn segments(loops: &[Vec]) -> Vec<(u32, u32)> { + let mut s = Vec::new(); + for l in loops { + for (i, &a) in l.iter().enumerate() { + if let Some(&b) = l.get((i + 1) % l.len()) { + s.push((a, b)); + } + } + } + s +} + +/// `c` lies within the bounding box of `a`, `b` (for collinear points: on the segment). +fn on_segment(a: [f64; 2], b: [f64; 2], c: [f64; 2]) -> bool { + a[0].min(b[0]) <= c[0] && c[0] <= a[0].max(b[0]) && a[1].min(b[1]) <= c[1] && c[1] <= a[1].max(b[1]) +} + +/// Whether two boundary segments conflict: they cross or touch, or (when they share an end +/// point) they fold back on each other. +fn segments_conflict(r: &Region, s: (u32, u32), t: (u32, u32)) -> Result { + let (a, b, c, d) = (r.point(s.0)?, r.point(s.1)?, r.point(t.0)?, r.point(t.1)?); + let shared = [(s.0 == t.0, b, d, a), (s.0 == t.1, b, c, a), (s.1 == t.0, a, d, b), (s.1 == t.1, a, c, b)]; + if let Some(&(_, x, y, o)) = shared.iter().find(|e| e.0) { + // Segments o–x and o–y: a conflict when collinear and on the same side of o. + if orient2d(o, x, y) != Sign::Zero { + return Ok(false); + } + let (o, x, y) = (o.coords(), x.coords(), y.coords()); + let dot = (x[0] - o[0]) * (y[0] - o[0]) + (x[1] - o[1]) * (y[1] - o[1]); + return Ok(dot > 0.0); + } + let (o1, o2, o3, o4) = (orient2d(a, b, c), orient2d(a, b, d), orient2d(c, d, a), orient2d(c, d, b)); + let proper = (o1 * o2) == Sign::Negative && (o3 * o4) == Sign::Negative; + let (pa, pb, pc, pd) = (a.coords(), b.coords(), c.coords(), d.coords()); + Ok(proper + || (o1.is_zero() && on_segment(pa, pb, pc)) + || (o2.is_zero() && on_segment(pa, pb, pd)) + || (o3.is_zero() && on_segment(pc, pd, pa)) + || (o4.is_zero() && on_segment(pc, pd, pb))) +} + +/// No two boundary segments cross, touch or overlap. Segments are bucketed on a uniform grid +/// over the bounding box; each pair is tested once, in the first cell both bounding boxes cover. +pub(crate) fn check_simple(r: &Region) -> Result<()> { + let segs = segments(&r.loops); + let n = segs.len(); + let (lo, hi) = bbox(r.pts.iter().copied()); + let g = ((n as f64).sqrt().ceil() as usize).clamp(1, 1024); + let cell = |x: f64, axis: usize| -> usize { + let (l, h) = (lo[axis], hi[axis]); + if !(h > l) { + return 0; + } + (((x - l) / (h - l) * g as f64) as usize).min(g - 1) + }; + let mut boxes = Vec::with_capacity(n); + let mut grid: Vec> = vec![Vec::new(); g * g]; + let mut work = 0usize; + for (i, &(a, b)) in segs.iter().enumerate() { + let (pa, pb) = (r.point(a)?.coords(), r.point(b)?.coords()); + let bx = [cell(pa[0].min(pb[0]), 0), cell(pa[1].min(pb[1]), 1), cell(pa[0].max(pb[0]), 0), cell(pa[1].max(pb[1]), 1)]; + for cy in bx[1]..=bx[3] { + for cx in bx[0]..=bx[2] { + work += 1; + if let Some(c) = grid.get_mut(cy * g + cx) { + c.push(i as u32); + } + } + } + if work > MAX_CHECK_WORK { + return Err(Error::Limit("boundary intersection check exceeded its work limit".into())); + } + boxes.push(bx); + } + for (ci, members) in grid.iter().enumerate() { + let (cx, cy) = (ci % g, ci / g); + for (k, &i) in members.iter().enumerate() { + for &j in members.get(k + 1..).unwrap_or(&[]) { + let (Some(bi), Some(bj), Some(&si), Some(&sj)) = + (boxes.get(i as usize), boxes.get(j as usize), segs.get(i as usize), segs.get(j as usize)) + else { + continue; + }; + if bi[0].max(bj[0]) != cx || bi[1].max(bj[1]) != cy { + continue; + } + work += 1; + if work > MAX_CHECK_WORK { + return Err(Error::Limit("boundary intersection check exceeded its work limit".into())); + } + if segments_conflict(r, si, sj)? { + return Err(Error::Invalid(format!( + "boundary edges cross, touch or overlap: {}–{} and {}–{}", + locate(r, si.0), + locate(r, si.1), + locate(r, sj.0), + locate(r, sj.1) + ))); + } + } + } + } + Ok(()) +} + +/// The exact orientation of a simple loop: the sign of the turn at its lexicographically +/// smallest point, a vertex of its convex hull. +pub(crate) fn loop_orientation(r: &Region, l: &[u32]) -> Result { + let (i, _) = l + .iter() + .enumerate() + .filter_map(|(i, &v)| r.pts.get(v as usize).map(|p| (i, p))) + .min_by(|x, y| lex(x.1, y.1)) + .ok_or_else(|| Error::Invalid("empty loop".into()))?; + let n = l.len(); + let at = |k: usize| l.get(k % n).copied().ok_or_else(|| Error::Failed("loop index".into())); + Ok(orient2d(r.point(at(i + n - 1)?)?, r.point(at(i)?)?, r.point(at(i + 1)?)?)) +} + +/// The outer loop is counterclockwise and the holes clockwise (no repair). +pub(crate) fn check_orientation(r: &Region) -> Result<()> { + for (k, l) in r.loops.iter().enumerate() { + let want = if k == 0 { Sign::Positive } else { Sign::Negative }; + if loop_orientation(r, l)? != want { + return Err(Error::Invalid(format!("loop {k} has the wrong orientation"))); + } + } + Ok(()) +} + +/// Outer loop counterclockwise, holes clockwise. +fn orient(r: &mut Region) -> Result<()> { + for k in 0..r.loops.len() { + let l = r.loops.get(k).cloned().unwrap_or_default(); + let want = if k == 0 { Sign::Positive } else { Sign::Negative }; + let s = loop_orientation(r, &l)?; + if s == Sign::Zero { + return Err(Error::Invalid(format!("loop {k} is degenerate (all points collinear)"))); + } + if s != want + && let Some(l) = r.loops.get_mut(k) + { + l.reverse(); + } + } + Ok(()) +} + +/// Winding number of loop `l` around `p` (exact: half-open crossing rule with `orient2d`). +pub(crate) fn winding(r: &Region, l: &[u32], p: &Point2) -> Result { + let y = p.coords()[1]; + let mut w = 0i64; + for (i, &a) in l.iter().enumerate() { + let b = l.get((i + 1) % l.len()).copied().unwrap_or(a); + let (pa, pb) = (r.point(a)?, r.point(b)?); + let (ya, yb) = (pa.coords()[1], pb.coords()[1]); + if ya <= y && yb > y && orient2d(pa, pb, p) == Sign::Positive { + w += 1; + } else if yb <= y && ya > y && orient2d(pa, pb, p) == Sign::Negative { + w -= 1; + } + } + Ok(w) +} + +/// Every hole lies inside the outer loop and outside every other hole. The loops do not +/// intersect, so one point per hole decides. +pub(crate) fn check_nesting(r: &Region) -> Result<()> { + let holes = r.loops.len().saturating_sub(1); + if holes.saturating_mul(r.pts.len()) > MAX_CHECK_WORK { + return Err(Error::Limit(format!("{holes} holes on {} boundary points exceed the nesting-check limit", r.pts.len()))); + } + let outer = r.loops.first().cloned().unwrap_or_default(); + for (k, hole) in r.loops.iter().enumerate().skip(1) { + let first = hole.first().copied().ok_or_else(|| Error::Invalid(format!("loop {k} is empty")))?; + let p = r.point(first)?; + if winding(r, &outer, p)? == 0 { + return Err(Error::Invalid(format!("hole loop {k} lies outside the outer loop"))); + } + for (m, other) in r.loops.iter().enumerate().skip(1) { + if m != k && winding(r, other, p)? != 0 { + return Err(Error::Invalid(format!("hole loop {k} lies inside hole loop {m}"))); + } + } + } + Ok(()) +} diff --git a/crates/mesh/src/quad/free.rs b/crates/mesh/src/quad/free.rs new file mode 100644 index 0000000..fdb3c1e --- /dev/null +++ b/crates/mesh/src/quad/free.rs @@ -0,0 +1,361 @@ +//! The free all-quadrilateral mesher: half-resolution constrained Delaunay triangulation, +//! recombination into quadrilaterals, one uniform split that lands on the prescribed boundary +//! nodes, then constrained smoothing. See the module docs of [`super`] for the method and its +//! sources. + +use std::collections::BTreeMap; + +use super::cdt::{CLEARANCE, Cdt, Triangulation}; +use super::domain::{self, Region}; +use super::guard::Guard; +use super::{QuadMesh, QuadOptions, finish}; +use crate::predicates::{Sign, orient2d}; +use crate::quality::{planar_corner, turns_left}; +use crate::size::SizeField; +use crate::{Error, Result}; + +/// Mesh the region bounded by `loops` (the first is the outer boundary, the rest are holes; +/// either orientation) with quadrilaterals only. See [`super::free_quad_mesh`]. +pub(crate) fn mesh(loops: &[Vec<[f64; 2]>], size: &dyn SizeField, opts: &QuadOptions) -> Result { + let region = domain::region(loops)?; + for (k, l) in region.loops.iter().enumerate() { + if l.len() % 2 == 1 { + return Err(Error::Invalid(format!( + "loop {k} has an odd number ({}) of boundary points; an all-quadrilateral mesh needs an even number on every loop \ + (choose the edge divisions with modelcraft_mesh::size::loop_divisions or even_divisions)", + l.len() + ))); + } + } + if region.loops.len() == 1 && region.loops.first().is_some_and(|l| l.len() == 4) { + let l = region.loops.first().cloned().unwrap_or_default(); + let quad = [l.first(), l.get(1), l.get(2), l.get(3)]; + let quad = match quad { + [Some(&a), Some(&b), Some(&c), Some(&d)] => [a, b, c, d], + _ => return Err(Error::Failed("quad loop".into())), + }; + return finish(region.pts.clone(), vec![quad], region.pts.len(), opts, Default::default()); + } + if let Some((k, l)) = region.loops.iter().enumerate().find(|(_, l)| l.len() < 6) { + return Err(Error::Invalid(format!( + "loop {k} has {} boundary points; the free quad mesher needs at least 6 on every loop of a region with more than one loop \ + (split its sides)", + l.len() + ))); + } + + // Both parity classes, preferred first; the second is tried when the first cannot be made + // valid. + let mut last = None; + for coarse in coarse_regions(®ion)? { + match attempt(®ion, coarse, size, opts) { + Ok(mesh) => return Ok(mesh), + Err(e @ Error::Failed(_)) => last = Some(e), + Err(e) => return Err(e), + } + } + Err(last.unwrap_or_else(|| Error::Failed("no half-resolution boundary".into()))) +} + +/// One run of the pipeline on one half-resolution region. +fn attempt(region: &Region, coarse: Coarse, size: &dyn SizeField, opts: &QuadOptions) -> Result { + let Coarse { region: coarse, to_input, skipped, before } = coarse; + let chords = chords(&coarse, &skipped)?; + // For each coarse boundary edge (coarse ids, smaller first), the input point it skips. + let boundary_mid: BTreeMap<(u32, u32), u32> = chords.iter().map(|&(a, b, m)| ((a.min(b), a.max(b)), m)).collect(); + let guard = + Guard::new(&chords.iter().map(|&(a, b, m)| Ok((*coarse.point(a)?, *coarse.point(b)?, *region.point(m)?))).collect::>>()?)?; + let need = fan_needs(region, &to_input, &before, &skipped)?; + let max_tris = opts.max_elements / 2 + 4 * region.pts.len() + 16; + let mut cdt = Cdt::build(&coarse, max_tris)?; + cdt.refine(size, 2.0, &|p, h| guard.allows(p, CLEARANCE * h))?; + cdt.ensure_fans(&need)?; + let tri = cdt.into_triangulation(); + let (quads, leftovers) = recombine(&tri, &need)?; + let count = 4 * quads.len() + 3 * leftovers.len(); + if count > opts.max_elements { + return Err(Error::Limit(format!("the mesh would have {count} quadrilaterals; the limit is {}", opts.max_elements))); + } + let (nodes, elements) = split(region, &tri, &to_input, &boundary_mid, &quads, &leftovers)?; + let stats = super::QuadStats { recombined: quads.len(), leftover_triangles: leftovers.len(), ..Default::default() }; + finish(nodes, elements, region.pts.len(), opts, stats) +} + +/// `1 − cos θ` of the turn at each loop position (0 on a straight run, 2 on a full reversal). +fn turns(region: &Region, l: &[u32]) -> Result> { + let n = l.len(); + let at = |i: usize| -> Result<[f64; 2]> { + let v = l.get(i % n).copied().ok_or_else(|| Error::Failed("loop index".into()))?; + region.pts.get(v as usize).copied().ok_or_else(|| Error::Failed("point index".into())) + }; + (0..n) + .map(|i| { + let (a, b, c) = (at(i + n - 1)?, at(i)?, at(i + 1)?); + let (u, w) = ([b[0] - a[0], b[1] - a[1]], [c[0] - b[0], c[1] - b[1]]); + let den = (u[0] * u[0] + u[1] * u[1]).sqrt() * (w[0] * w[0] + w[1] * w[1]).sqrt(); + let cos = (u[0] * w[0] + u[1] * w[1]) / den; + Ok(if den > 0.0 && cos.is_finite() { 1.0 - cos } else { 0.0 }) + }) + .collect() +} + +/// The half-resolution regions to try: every other point of each loop, first keeping the parity +/// class with the sharper turns (so loop corners survive when they sit at an even spacing), then +/// the other class of every loop. Choices whose chords cross or flip are dropped. +fn coarse_regions(region: &Region) -> Result> { + let mut preferred = Vec::with_capacity(region.loops.len()); + for l in ®ion.loops { + let t = turns(region, l)?; + let even: f64 = t.iter().step_by(2).sum(); + let odd: f64 = t.iter().skip(1).step_by(2).sum(); + preferred.push(usize::from(odd > even)); + } + let mut out = Vec::new(); + let mut last = None; + for flip in [0, 1] { + let offsets: Vec = preferred.iter().map(|&o| o ^ flip).collect(); + let coarse = pick(region, &offsets)?; + match check_coarse(&coarse.region) { + Ok(()) => out.push(coarse), + Err(e) => last = Some(e), + } + } + if out.is_empty() { + return Err(Error::Failed(format!( + "the boundary is too coarse for the free quad mesher: the chords joining every other boundary point cross or flip ({}); \ + refine the boundary divisions", + last.map(|e| e.to_string()).unwrap_or_default() + ))); + } + Ok(out) +} + +/// The half-resolution region, with each coarse point's input id and the input id of the +/// boundary point after it (the one the coarse edge skips). +struct Coarse { + region: Region, + to_input: Vec, + skipped: Vec, + /// The input point before each coarse point. + before: Vec, +} + +fn pick(region: &Region, offsets: &[usize]) -> Result { + let (mut pts, mut p2, mut loops, mut to_input, mut skipped, mut before) = + (Vec::new(), Vec::new(), Vec::new(), Vec::new(), Vec::new(), Vec::new()); + for (l, &off) in region.loops.iter().zip(offsets) { + let mut ids = Vec::with_capacity(l.len() / 2); + for pos in (off..l.len()).step_by(2) { + let (v, after, prior) = match (l.get(pos), l.get((pos + 1) % l.len()), l.get((pos + l.len() - 1) % l.len())) { + (Some(&v), Some(&a), Some(&b)) => (v, a, b), + _ => return Err(Error::Failed("loop index".into())), + }; + ids.push(pts.len() as u32); + pts.push(region.pts.get(v as usize).copied().ok_or_else(|| Error::Failed("point index".into()))?); + p2.push(*region.point(v)?); + to_input.push(v); + skipped.push(after); + before.push(prior); + } + loops.push(ids); + } + Ok(Coarse { region: Region { pts, p2, loops }, to_input, skipped, before }) +} + +/// The coarse loops must be simple, keep their orientation and nesting. +fn check_coarse(coarse: &Region) -> Result<()> { + domain::check_simple(coarse)?; + domain::check_orientation(coarse)?; + domain::check_nesting(coarse) +} + +/// The coarse boundary edges `(a, b, m)`: coarse ids in loop order and the input point `m` the +/// edge skips. +fn chords(coarse: &Region, skipped: &[u32]) -> Result> { + let mut out = Vec::with_capacity(coarse.pts.len()); + for (a, b) in domain::segments(&coarse.loops) { + out.push((a, b, skipped.get(a as usize).copied().ok_or_else(|| Error::Failed("coarse id".into()))?)); + } + Ok(out) +} + +/// How many coarse elements each coarse boundary vertex needs around it. A coarse vertex `v` +/// gets one final quadrilateral per incident coarse element, and the corner of the first and last +/// of them at `v` opens towards the prescribed neighbours of `v`. A single element would span the +/// whole boundary angle at `v`, so a straight or nearly straight (interior angle above 150°) or +/// reflex vertex needs two, and one above 300° needs three. Decided with the exact `orient2d` and +/// the cosine of the turn. +fn fan_needs(region: &Region, to_input: &[u32], before: &[u32], after: &[u32]) -> Result> { + let mut need = Vec::with_capacity(to_input.len()); + for ((&v, &a), &b) in to_input.iter().zip(before).zip(after) { + let (pa, pv, pb) = (region.point(a)?, region.point(v)?, region.point(b)?); + let o = orient2d(pa, pv, pb); + let (a, v, b) = (pa.coords(), pv.coords(), pb.coords()); + let (u, w) = ([v[0] - a[0], v[1] - a[1]], [b[0] - v[0], b[1] - v[1]]); + let cos = (u[0] * w[0] + u[1] * w[1]) / ((u[0] * u[0] + u[1] * u[1]).sqrt() * (w[0] * w[0] + w[1] * w[1]).sqrt()); + // Interior angle 180° − turn: above 150° when the turn is below 30° (cos > √3/2). + need.push(match o { + Sign::Positive if cos > 0.866_025_403_784_438_6 => 2, + Sign::Positive => 1, + Sign::Zero => 2, + Sign::Negative if cos < -0.5 => 3, + Sign::Negative => 2, + }); + } + Ok(need) +} + +/// The smallest corner scaled Jacobian of a planar quadrilateral (corner `i`: +/// `cross(p[i+1] − p[i], p[i−1] − p[i]) / (|…| |…|)`), in `[-1, 1]`. +/// A degenerate corner (zero-length edge) scores −1, the worst, so optimizers move away from it. +pub(crate) fn min_scaled_jacobian(p: [[f64; 2]; 4]) -> f64 { + (0..4).fold(f64::INFINITY, |m, i| m.min(planar_corner(&p, i).1.unwrap_or(-1.0))) +} + +/// The smallest corner scaled Jacobian of a recombined quadrilateral (sin 15°): a flatter corner +/// would survive the split as a nearly straight angle, which rounding of the mid nodes can turn +/// degenerate. +const MIN_RECOMBINED_QUALITY: f64 = 0.258_819_045_102_520_74; + +/// Greedy recombination of adjacent triangle pairs into strictly convex quadrilaterals, best +/// first by the smallest corner scaled Jacobian (Lee and Lo 1994; Knupp 2001). Returns the +/// quadrilaterals (coarse ids, counterclockwise) and the triangles left over. +/// +/// A merge removes one element around the two shared vertices; it is skipped when a boundary +/// vertex would drop below the elements it needs ([`fan_needs`]). +fn recombine(tri: &Triangulation, need: &[usize]) -> Result<(Vec<[u32; 4]>, Vec<[u32; 3]>)> { + let mut count = vec![0usize; tri.pts.len()]; + for t in &tri.tris { + for &v in t { + if let Some(c) = count.get_mut(v as usize) { + *c += 1; + } + } + } + let mut edges: Vec<(u32, u32, u32, usize)> = Vec::with_capacity(3 * tri.tris.len()); + for (t, v) in tri.tris.iter().enumerate() { + for i in 0..3 { + let (a, b) = (v[(i + 1) % 3], v[(i + 2) % 3]); + edges.push((a.min(b), a.max(b), t as u32, i)); + } + } + edges.sort_unstable(); + let pt = |v: u32| tri.p2.get(v as usize).ok_or_else(|| Error::Failed("point index".into())); + let xy = |v: u32| tri.pts.get(v as usize).copied().ok_or_else(|| Error::Failed("point index".into())); + let mut candidates: Vec<(f64, u32, u32, [u32; 4], (u32, u32))> = Vec::new(); + for w in edges.windows(2) { + let [(a0, b0, t1, i), (a1, b1, t2, j)] = w else { continue }; + if (a0, b0) != (a1, b1) { + continue; + } + let (v1, v2) = match (tri.tris.get(*t1 as usize), tri.tris.get(*t2 as usize)) { + (Some(x), Some(y)) => (x, y), + _ => continue, + }; + let (p, q, r, s) = (v1[*i], v1[(i + 1) % 3], v1[(i + 2) % 3], v2[*j]); + let quad = [p, q, s, r]; + if !turns_left([pt(p)?, pt(q)?, pt(s)?, pt(r)?]) { + continue; + } + let quality = min_scaled_jacobian([xy(p)?, xy(q)?, xy(s)?, xy(r)?]); + if quality < MIN_RECOMBINED_QUALITY { + continue; + } + candidates.push((quality, *t1, *t2, quad, (q, r))); + } + candidates.sort_by(|x, y| y.0.total_cmp(&x.0).then((x.1, x.2).cmp(&(y.1, y.2)))); + let mut used = vec![false; tri.tris.len()]; + let mut quads = Vec::new(); + for (_, t1, t2, quad, (q, r)) in candidates { + let (Some(&u1), Some(&u2)) = (used.get(t1 as usize), used.get(t2 as usize)) else { continue }; + if u1 || u2 { + continue; + } + let keeps = |v: u32| count.get(v as usize).copied().unwrap_or(0) > need.get(v as usize).copied().unwrap_or(0); + if !keeps(q) || !keeps(r) { + continue; + } + for v in [q, r] { + if let Some(c) = count.get_mut(v as usize) { + *c -= 1; + } + } + for t in [t1, t2] { + if let Some(u) = used.get_mut(t as usize) { + *u = true; + } + } + quads.push(quad); + } + let leftovers = tri.tris.iter().zip(&used).filter(|(_, u)| !**u).map(|(t, _)| *t).collect(); + Ok((quads, leftovers)) +} + +/// The uniform split: every coarse edge gets a mid node (on the boundary, the prescribed point +/// it skipped), every element a centre node; a quadrilateral becomes four quadrilaterals, a +/// triangle three. Node order: the input boundary points (input order), the interior coarse +/// vertices, the interior edge mid nodes, the centres. +fn split( + region: &Region, + tri: &Triangulation, + to_input: &[u32], + boundary_mid: &BTreeMap<(u32, u32), u32>, + quads: &[[u32; 4]], + leftovers: &[[u32; 3]], +) -> Result<(Vec<[f64; 2]>, Vec<[u32; 4]>)> { + let mut nodes = region.pts.clone(); + let mut final_id = vec![u32::MAX; tri.pts.len()]; + for (c, &input) in to_input.iter().enumerate() { + if let Some(f) = final_id.get_mut(c) { + *f = input; + } + } + let xy = |v: u32| tri.pts.get(v as usize).copied().ok_or_else(|| Error::Failed("point index".into())); + let elements: Vec<&[u32]> = quads.iter().map(|q| q.as_slice()).chain(leftovers.iter().map(|t| t.as_slice())).collect(); + for e in &elements { + for &v in e.iter() { + let f = final_id.get_mut(v as usize).ok_or_else(|| Error::Failed("vertex id".into()))?; + if *f == u32::MAX { + *f = nodes.len() as u32; + nodes.push(xy(v)?); + } + } + } + let fid = |v: u32| final_id.get(v as usize).copied().filter(|&f| f != u32::MAX).ok_or_else(|| Error::Failed("vertex id".into())); + let mut mids: BTreeMap<(u32, u32), u32> = BTreeMap::new(); + for e in &elements { + for (i, &a) in e.iter().enumerate() { + let b = e.get((i + 1) % e.len()).copied().ok_or_else(|| Error::Failed("element".into()))?; + let key = (a.min(b), a.max(b)); + if let std::collections::btree_map::Entry::Vacant(slot) = mids.entry(key) { + let id = match boundary_mid.get(&key) { + Some(&m) => m, + None => { + let (pa, pb) = (xy(a)?, xy(b)?); + nodes.push([0.5 * (pa[0] + pb[0]), 0.5 * (pa[1] + pb[1])]); + nodes.len() as u32 - 1 + } + }; + slot.insert(id); + } + } + } + let mid = |a: u32, b: u32| mids.get(&(a.min(b), a.max(b))).copied().ok_or_else(|| Error::Failed("mid node".into())); + let mut out = Vec::with_capacity(4 * quads.len() + 3 * leftovers.len()); + for e in &elements { + let k = e.len(); + let mut c = [0.0, 0.0]; + for &v in e.iter() { + let p = xy(v)?; + c = [c[0] + p[0], c[1] + p[1]]; + } + nodes.push([c[0] / k as f64, c[1] / k as f64]); + let centre = nodes.len() as u32 - 1; + for (i, &v) in e.iter().enumerate() { + let next = e.get((i + 1) % k).copied().ok_or_else(|| Error::Failed("element".into()))?; + let prev = e.get((i + k - 1) % k).copied().ok_or_else(|| Error::Failed("element".into()))?; + out.push([fid(v)?, mid(v, next)?, centre, mid(prev, v)?]); + } + } + Ok((nodes, out)) +} diff --git a/crates/mesh/src/quad/guard.rs b/crates/mesh/src/quad/guard.rs new file mode 100644 index 0000000..3d44b2e --- /dev/null +++ b/crates/mesh/src/quad/guard.rs @@ -0,0 +1,121 @@ +//! Where the half-resolution triangulation may place Steiner points. +//! +//! The half-resolution region replaces each pair of boundary edges `a–m–b` by the chord `a–b`; +//! the split later puts `m` back as the chord's mid node. The final mesh therefore covers the true +//! region, not the coarse one, and a Steiner point must respect the true boundary: it may not lie +//! in a *cap* (the triangle `a, b, m` when `m` lies inside the coarse region, which the true region +//! excludes), and it keeps a clearance from the true edges `a–m`, `m–b` of at least half the chord +//! (so the elements along the boundary keep the boundary's spacing and can absorb `m`'s offset +//! from the chord). Lookups go through a uniform grid; the cap test is exact (`orient2d`). + +use crate::predicates::{Point2, Sign, orient2d}; +use crate::{Error, Result}; + +/// Above this many grid cells a clearance query scans every segment instead. +const MAX_QUERY_CELLS: usize = 1024; + +pub(crate) struct Guard { + lo: [f64; 2], + step: [f64; 2], + g: usize, + /// True boundary edges with their own clearance radius. + segs: Vec<([f64; 2], [f64; 2], f64)>, + /// Counterclockwise caps `(a, b, m)`. + caps: Vec<[Point2; 3]>, + seg_cells: Vec>, + cap_cells: Vec>, +} + +fn distance(a: [f64; 2], b: [f64; 2], p: [f64; 2]) -> f64 { + let (d, w) = ([b[0] - a[0], b[1] - a[1]], [p[0] - a[0], p[1] - a[1]]); + let len2 = d[0] * d[0] + d[1] * d[1]; + let t = if len2 > 0.0 { ((w[0] * d[0] + w[1] * d[1]) / len2).clamp(0.0, 1.0) } else { 0.0 }; + ((w[0] - t * d[0]).powi(2) + (w[1] - t * d[1]).powi(2)).sqrt() +} + +impl Guard { + /// The guard of chords `(a, b, m)`: coarse boundary edge `a → b` (region on its left) and the + /// skipped boundary point `m`. + pub fn new(chords: &[(Point2, Point2, Point2)]) -> Result { + let (lo, hi) = super::domain::bbox(chords.iter().flat_map(|(a, b, m)| [a.coords(), b.coords(), m.coords()])); + if chords.is_empty() { + return Err(Error::Failed("guard without a boundary".into())); + } + let g = ((2.0 * chords.len() as f64).sqrt().ceil() as usize).clamp(1, 1024); + let step = [((hi[0] - lo[0]) / g as f64).max(f64::MIN_POSITIVE), ((hi[1] - lo[1]) / g as f64).max(f64::MIN_POSITIVE)]; + let mut guard = + Guard { lo, step, g, segs: Vec::new(), caps: Vec::new(), seg_cells: vec![Vec::new(); g * g], cap_cells: vec![Vec::new(); g * g] }; + for (a, b, m) in chords { + let (pa, pb, pm) = (a.coords(), b.coords(), m.coords()); + let radius = 0.5 * ((pb[0] - pa[0]).powi(2) + (pb[1] - pa[1]).powi(2)).sqrt(); + for (s, t) in [(pa, pm), (pm, pb)] { + let id = guard.segs.len() as u32; + guard.segs.push((s, t, radius)); + let (x0, y0) = guard.cell([s[0].min(t[0]) - radius, s[1].min(t[1]) - radius]); + let (x1, y1) = guard.cell([s[0].max(t[0]) + radius, s[1].max(t[1]) + radius]); + guard.register(true, id, (x0, y0), (x1, y1)); + } + if orient2d(a, b, m) == Sign::Positive { + let id = guard.caps.len() as u32; + guard.caps.push([*a, *b, *m]); + let (x0, y0) = guard.cell([pa[0].min(pb[0]).min(pm[0]), pa[1].min(pb[1]).min(pm[1])]); + let (x1, y1) = guard.cell([pa[0].max(pb[0]).max(pm[0]), pa[1].max(pb[1]).max(pm[1])]); + guard.register(false, id, (x0, y0), (x1, y1)); + } + } + Ok(guard) + } + + fn cell(&self, p: [f64; 2]) -> (usize, usize) { + let f = |x: f64, lo: f64, step: f64| { + let c = ((x - lo) / step).floor(); + if c.is_nan() || c < 0.0 { 0 } else { (c as usize).min(self.g - 1) } + }; + (f(p[0], self.lo[0], self.step[0]), f(p[1], self.lo[1], self.step[1])) + } + + fn register(&mut self, seg: bool, id: u32, from: (usize, usize), to: (usize, usize)) { + let cells = if seg { &mut self.seg_cells } else { &mut self.cap_cells }; + for y in from.1..=to.1 { + for x in from.0..=to.0 { + if let Some(c) = cells.get_mut(y * self.g + x) { + c.push(id); + } + } + } + } + + /// Whether a Steiner point may go at `p`: outside every cap, and farther than both the + /// segment's own radius and `clearance` from every true boundary edge. + pub fn allows(&self, p: &Point2, clearance: f64) -> bool { + let c = p.coords(); + let (cx, cy) = self.cell(c); + if let Some(caps) = self.cap_cells.get(cy * self.g + cx) { + for &k in caps { + if let Some([a, b, m]) = self.caps.get(k as usize) + && orient2d(a, b, p) != Sign::Negative + && orient2d(b, m, p) != Sign::Negative + && orient2d(m, a, p) != Sign::Negative + { + return false; + } + } + } + let too_close = |&(a, b, r): &([f64; 2], [f64; 2], f64)| distance(a, b, c) < r.max(clearance); + let (x0, y0) = self.cell([c[0] - clearance, c[1] - clearance]); + let (x1, y1) = self.cell([c[0] + clearance, c[1] + clearance]); + if (x1 - x0 + 1) * (y1 - y0 + 1) > MAX_QUERY_CELLS { + return !self.segs.iter().any(too_close); + } + for y in y0..=y1 { + for x in x0..=x1 { + for &k in self.seg_cells.get(y * self.g + x).map(Vec::as_slice).unwrap_or(&[]) { + if self.segs.get(k as usize).is_some_and(too_close) { + return false; + } + } + } + } + true + } +} diff --git a/crates/mesh/src/quad/mod.rs b/crates/mesh/src/quad/mod.rs new file mode 100644 index 0000000..313787e --- /dev/null +++ b/crates/mesh/src/quad/mod.rs @@ -0,0 +1,287 @@ +//! Planar all-quadrilateral meshing: mapped (transfinite) meshes of four-sided regions and free +//! meshes of polygons with holes, on prescribed boundary nodes. +//! +//! These are the face meshers under the mapped and swept hexahedral meshers and shell meshing. +//! Both take the boundary nodes as given (shared edges must share nodes), put them first in the +//! output in input order, and return counterclockwise `Quad4` connectivity whose every element is +//! proven valid by the exact corner test ([`crate::quality::planar_quad_validity`]). +//! +//! # Mapped: [`mapped_quad_mesh`], [`transfinite_grid`] +//! Transfinite interpolation (Gordon and Hall 1973), see [`tfi`]. Opposite sides need equal +//! division counts; a region whose TFI folds (strongly concave sides) is an error, not a tangled +//! mesh. +//! +//! # Free: [`free_quad_mesh`] +//! **Method.** An indirect tri-to-quad method in three steps: +//! 1. Keep every other boundary point of each loop (the parity class with the sharper turns, so +//! corners at even spacing survive) and build the constrained Delaunay triangulation of that +//! half-resolution region, refined to twice the target size ([`cdt`]: Bowyer–Watson, Sloan +//! 1993 recovery, Ruppert 1995 refinement without segment splitting). +//! 2. Recombine adjacent triangles into strictly convex quadrilaterals, greedily, best first by +//! the smallest corner scaled Jacobian (Lee and Lo, *A new scheme for the generation of a +//! graded quadrilateral mesh*, Computers & Structures 52(5):847–857, 1994). +//! 3. Split every element once: a quadrilateral into four, a leftover triangle into three, through +//! edge mid nodes and a centre node. On a boundary edge the mid node is the prescribed point the +//! half-resolution loop skipped, so the final boundary nodes are exactly the input points. The +//! split is the Catmull–Clark topology (Catmull and Clark, CAD 10(6):350–355, 1978) and appears +//! as the all-quad completion of recombination in Lee and Lo 1994. +//! +//! Then constrained Laplacian smoothing ([`smooth`]) and the exact validity gate. +//! +//! **Why this method.** The alternatives considered, from their papers: +//! - *Q-Morph* (Owen, Staten, Canann, Saigal, *Q-Morph: an indirect approach to advancing front +//! quad meshing*, IJNME 44(9):1317–1340, 1999) gives better-aligned rows along the boundary, but +//! needs front-edge state classification, side-edge recovery by swaps, seam closing and special +//! handling where fronts collide; the paper reports the cases that need it, and each is a place +//! to fail on concave faces with many holes. An odd front remainder still ends in a triangle +//! to be removed. +//! - *Paving* (Blacker and Stephenson, IJNME 32(4):811–847, 1991) has the same front-collision +//! problem with no background triangulation to fall back on; it is the highest-risk choice. +//! - *Blossom-Quad* (Remacle, Lambrechts, Seny, Marchandise, Johnen, Geuzaine, *Blossom-Quad: a +//! non-uniform quadrilateral mesh generator using a minimum-cost perfect-matching algorithm*, +//! IJNME 89(9):1102–1119, 2012; paper only, Gmsh's code is GPL) gives the best matching, but a +//! perfect matching of the triangle graph need not exist, so it still needs a fallback, and it +//! needs Edmonds' weighted blossom algorithm (cubic in the worst case). +//! +//! The split method is **all-quad by construction**: any triangulation of the half-resolution +//! region becomes quadrilaterals only, whatever the matching found, so the only requirement is an +//! even number of boundary points on every loop. Robustness on concave regions with holes is +//! that of the constrained Delaunay triangulation, which exists for every valid input and is +//! built here with exact predicates. The implementation risk is low: one triangulator (which the +//! surface and tetrahedral work need anyway), a greedy matching and a local split. The cost is +//! quality: the split creates irregular nodes (valence 3 at triangle centres, 6 at some coarse +//! vertices) and triangle-born quads have a 120° corner; smoothing removes most of the +//! distortion. If that proves too weak, the greedy matching can be replaced by Blossom-Quad's +//! minimum-cost matching without touching the rest (fewer leftover triangles, same guarantee). +//! +//! **Requirements and failure modes.** +//! - Every loop has an even number of points, at least 6 when there is more than one loop (a +//! single 4-point loop is meshed as one quadrilateral). An odd loop is an error, never repaired: +//! the points are prescribed by neighbouring faces, so adding one would break conformity. Choose +//! divisions with [`crate::size::loop_divisions`] or [`crate::size::even_divisions`]. +//! - The chords joining every other boundary point must not cross (a boundary too coarse for a +//! thin feature); both parity classes are tried before this is reported. +//! - Every output element is checked with exact predicates; a mesh that smoothing cannot make +//! valid is an error. +//! - Hostile input (NaN, coordinates beyond the model limit, repeated points, zero-length edges, +//! crossing or touching loops, holes outside the outer loop or inside another hole) is an +//! [`Error::Invalid`]; size and work limits are [`Error::Limit`]. + +mod cdt; +mod domain; +mod free; +mod guard; +mod smooth; +pub mod tfi; + +pub use domain::{MAX_BOUNDARY_POINTS, MAX_LOOPS}; +pub use tfi::{MAX_GRID_NODES, transfinite_grid}; + +use crate::predicates::{Point2, Sign, snap_to_domain}; +use crate::quality::{planar_corner, quad_axes, turns_left}; +use crate::size::SizeField; +use crate::{Error, Result}; + +/// Limits and smoothing of the quad meshers. +#[derive(Debug, Clone, PartialEq)] +pub struct QuadOptions { + /// The most quadrilaterals of one mesh ([`Error::Limit`] beyond). + pub max_elements: usize, + /// Smoothing sweeps of the free mesher (0 turns smoothing off). + pub smoothing_sweeps: usize, +} + +impl Default for QuadOptions { + fn default() -> Self { + QuadOptions { max_elements: 4_000_000, smoothing_sweeps: 10 } + } +} + +/// What a quad mesher produced and how good it is. +#[derive(Debug, Clone, PartialEq)] +pub struct QuadStats { + /// Smallest corner scaled Jacobian over the mesh (1 for squares). + pub min_scaled_jacobian: f64, + /// Largest skew (0 for rectangles). + pub max_skew: f64, + /// Largest aspect ratio (1 for squares). + pub max_aspect_ratio: f64, + /// Smoothing sweeps run, and whether smoothing converged before its sweep limit. + pub smoothing_sweeps: usize, + pub smoothing_converged: bool, + /// Free mesher: triangle pairs recombined, and triangles split on their own. + pub recombined: usize, + pub leftover_triangles: usize, +} + +impl Default for QuadStats { + fn default() -> Self { + QuadStats { + min_scaled_jacobian: 1.0, + max_skew: 0.0, + max_aspect_ratio: 1.0, + smoothing_sweeps: 0, + smoothing_converged: true, + recombined: 0, + leftover_triangles: 0, + } + } +} + +/// A planar all-quadrilateral mesh. +#[derive(Debug, Clone, PartialEq)] +pub struct QuadMesh { + /// Node coordinates. `nodes[..boundary_nodes]` are the input boundary points, unchanged and + /// in input order (loop by loop), so node `k` of the mesh is boundary point `k`. + pub nodes: Vec<[f64; 2]>, + /// Counterclockwise corner node ids (the `Quad4` node order). + pub quads: Vec<[u32; 4]>, + pub boundary_nodes: usize, + pub stats: QuadStats, +} + +/// A mapped mesh with its structure: node `(i, j)`, `0 ≤ i ≤ nu`, `0 ≤ j ≤ nv`, is +/// `mesh.nodes[grid[j·(nu + 1) + i]]`, with `i` running from corner 0 to corner 1 and `j` from +/// corner 1 to corner 2 (for a sweep: the source face's structured layout). +#[derive(Debug, Clone, PartialEq)] +pub struct MappedMesh { + pub mesh: QuadMesh, + pub nu: usize, + pub nv: usize, + pub grid: Vec, +} + +/// Free all-quadrilateral mesh of a planar region (see the module docs for the method). +/// +/// `loops[0]` is the outer boundary, the others are holes; each loop lists its boundary nodes +/// once (not closed), in either orientation, with an even count. `size` sets the interior +/// density; the boundary density is the given nodes. +pub fn free_quad_mesh(loops: &[Vec<[f64; 2]>], size: &dyn SizeField, opts: &QuadOptions) -> Result { + free::mesh(loops, size, opts) +} + +/// Mapped mesh of a four-sided region: `boundary` lists the region's boundary nodes once (either +/// orientation), and `corners` are the positions in it of the four corners in boundary order. +/// Opposite sides need the same number of divisions. The grid is transfinite interpolation +/// ([`transfinite_grid`]); a grid that folds is an error. +pub fn mapped_quad_mesh(boundary: &[[f64; 2]], corners: [usize; 4], opts: &QuadOptions) -> Result { + let region = domain::region(&[boundary.to_vec()])?; + let n = boundary.len(); + if let Some(c) = corners.iter().find(|&&c| c >= n) { + return Err(Error::Invalid(format!("corner position {c} is outside the boundary of {n} points"))); + } + let c0 = corners[0]; + let d = corners.map(|c| (c + n - c0) % n); + if !(d[0] == 0 && d[0] < d[1] && d[1] < d[2] && d[2] < d[3]) { + return Err(Error::Invalid("the corners must be four distinct positions in boundary order".into())); + } + let (nu, nv) = (d[1], d[2] - d[1]); + if d[3] - d[2] != nu || n - d[3] != nv { + return Err(Error::Invalid(format!("opposite sides need the same number of divisions: {} and {}, {} and {}", nu, d[3] - d[2], nv, n - d[3]))); + } + let at = |k: usize| region.pts.get((c0 + k) % n).copied().ok_or_else(|| Error::Failed("boundary index".into())); + let side = |from: usize, len: usize, forward: bool| -> Result> { + (0..=len).map(|k| if forward { at(from + k) } else { at(from + n - k) }).collect() + }; + // Boundary positions relative to corner 0: bottom 0..nu, right nu..nu+nv, top runs backwards + // from d[3], left backwards from 0 (= n). + let bottom = side(0, nu, true)?; + let right = side(d[1], nv, true)?; + let top = side(d[3], nu, false)?; + let left = side(n, nv, false)?; + let points = transfinite_grid(&bottom, &right, &top, &left)?; + // Node ids: boundary nodes keep their input positions; interior nodes follow, row by row. + let mut nodes = region.pts.clone(); + let mut grid = Vec::with_capacity(points.len()); + for j in 0..=nv { + for i in 0..=nu { + let pos = if j == 0 { + Some(i) + } else if i == nu { + Some(nu + j) + } else if j == nv { + Some(d[3] - i) + } else if i == 0 { + Some((n - j) % n) + } else { + None + }; + let id = match pos { + Some(p) => ((c0 + p) % n) as u32, + None => { + let x = points.get(j * (nu + 1) + i).copied().ok_or_else(|| Error::Failed("grid index".into()))?; + nodes.push(domain::admissible(x)?); + nodes.len() as u32 - 1 + } + }; + grid.push(id); + } + } + let ccw = domain::loop_orientation(®ion, &(0..n as u32).collect::>())? == Sign::Positive; + let g = |i: usize, j: usize| grid.get(j * (nu + 1) + i).copied().ok_or_else(|| Error::Failed("grid index".into())); + let mut quads = Vec::with_capacity(nu * nv); + for j in 0..nv { + for i in 0..nu { + let (a, b, c, e) = (g(i, j)?, g(i + 1, j)?, g(i + 1, j + 1)?, g(i, j + 1)?); + quads.push(if ccw { [a, b, c, e] } else { [a, e, c, b] }); + } + } + if quads.len() > opts.max_elements { + return Err(Error::Limit(format!("the mesh would have {} quadrilaterals; the limit is {}", quads.len(), opts.max_elements))); + } + let mesh = finish(nodes, quads, n, &QuadOptions { smoothing_sweeps: 0, ..opts.clone() }, QuadStats::default())?; + Ok(MappedMesh { mesh, nu, nv, grid }) +} + +/// Smooth (free nodes only), gate every element with the exact validity test and measure. +pub(crate) fn finish( + mut nodes: Vec<[f64; 2]>, + quads: Vec<[u32; 4]>, + boundary_nodes: usize, + opts: &QuadOptions, + mut stats: QuadStats, +) -> Result { + // Constructed nodes (mid nodes, centres) stay in the predicates' domain. + for n in nodes.iter_mut() { + *n = [snap_to_domain(n[0]), snap_to_domain(n[1])]; + } + if opts.smoothing_sweeps > 0 { + let (sweeps, converged) = smooth::smooth(&mut nodes, &quads, boundary_nodes, opts.smoothing_sweeps)?; + stats.smoothing_sweeps = sweeps; + stats.smoothing_converged = converged; + } + let mut invalid = 0usize; + let mut first_invalid = None; + let (mut sj, mut skew, mut aspect) = (f64::INFINITY, 0.0_f64, 1.0_f64); + for (e, q) in quads.iter().enumerate() { + let mut c = [[0.0; 2]; 4]; + for (k, &v) in q.iter().enumerate() { + c[k] = nodes.get(v as usize).copied().ok_or_else(|| Error::Failed("node index".into()))?; + } + let [a, b, d, f] = c.map(|x| Point2::new(x[0], x[1])); + if !turns_left([&a?, &b?, &d?, &f?]) { + invalid += 1; + first_invalid.get_or_insert(e); + } + sj = sj.min((0..4).fold(f64::INFINITY, |m, i| m.min(planar_corner(&c, i).1.unwrap_or(0.0)))); + let (r, k) = quad_axes(c.map(|x| [x[0], x[1], 0.0])); + aspect = aspect.max(r); + skew = skew.max(k); + } + if invalid > 0 { + let at = first_invalid.and_then(|e| quads.get(e)).and_then(|q| nodes.get(q[0] as usize)).copied().unwrap_or_default(); + return Err(Error::Failed(format!( + "{invalid} of {} quadrilaterals are inverted or degenerate (first near {at:?}); refine the boundary there or use a finer size", + quads.len() + ))); + } + if !quads.is_empty() { + stats.min_scaled_jacobian = sj; + stats.max_skew = skew; + stats.max_aspect_ratio = aspect; + } + Ok(QuadMesh { nodes, quads, boundary_nodes, stats }) +} + +#[cfg(test)] +mod tests; diff --git a/crates/mesh/src/quad/smooth.rs b/crates/mesh/src/quad/smooth.rs new file mode 100644 index 0000000..100b7a9 --- /dev/null +++ b/crates/mesh/src/quad/smooth.rs @@ -0,0 +1,250 @@ +//! Constrained Laplacian smoothing of a planar quadrilateral mesh. +//! +//! Each free node moves towards the average of its edge neighbours (Laplacian smoothing), and +//! the move is kept only when it does not lower the smallest scaled Jacobian of the node's +//! quadrilaterals and keeps them valid by the exact corner test ("smart" Laplacian: Freitag, +//! *On combining Laplacian and optimization-based mesh smoothing techniques*, AMD-Vol. 220, +//! ASME 1997; Canann, Tristano, Staten, *An approach to combined Laplacian and optimization-based +//! smoothing for triangular, quadrilateral, and quad-dominant meshes*, IMR 1998). Half and quarter +//! steps are tried when the full step is refused. Nodes of quadrilaterals that are already +//! invalid accept any step that raises their smallest corner area (see below). +//! +//! Laplacian smoothing cannot untangle a node whose average position is itself bad, and the +//! nodes next to a sharp prescribed boundary point often are. Every node whose quadrilaterals +//! remain inverted or poor (smallest scaled Jacobian below [`POOR`]) then gets optimization-based +//! smoothing (Freitag 1997) by a compass search (eight directions, step halved on failure, +//! bounded). An inverted patch is untangled first by maximizing its smallest corner area, which is +//! linear in the node's position, so its maximum lies in the feasible region when there is one +//! (Freitag and Plassmann, *Local optimization-based simplicial mesh untangling and +//! improvement*, IJNME 49:109–125, 2000); a valid patch then maximizes its smallest scaled +//! Jacobian without leaving the valid set. +//! +//! Gauss–Seidel order by node id; bounded sweeps; deterministic. + +use crate::predicates::{Point2, snap_to_domain}; +use crate::quality::{planar_corner, turns_left}; +use crate::{Error, Result}; + +/// Compressed adjacency: `items[start[v]..start[v + 1]]`. +struct Csr { + start: Vec, + items: Vec, +} + +impl Csr { + fn of(n: usize, pairs: &mut [(u32, u32)]) -> Csr { + pairs.sort_unstable(); + let mut start = vec![0usize; n + 1]; + let mut items = Vec::with_capacity(pairs.len()); + let mut prev = None; + for &(a, b) in pairs.iter() { + if prev == Some((a, b)) { + continue; + } + prev = Some((a, b)); + if let Some(s) = start.get_mut(a as usize + 1) { + *s += 1; + } + items.push(b); + } + for i in 0..n { + let before = start.get(i).copied().unwrap_or(0); + if let Some(s) = start.get_mut(i + 1) { + *s += before; + } + } + Csr { start, items } + } + + fn get(&self, v: usize) -> &[u32] { + match (self.start.get(v), self.start.get(v + 1)) { + (Some(&a), Some(&b)) => self.items.get(a..b).unwrap_or(&[]), + _ => &[], + } + } +} + +fn corners(nodes: &[[f64; 2]], q: &[u32; 4], v: u32, at: [f64; 2]) -> Result<[[f64; 2]; 4]> { + let mut c = [[0.0; 2]; 4]; + for (k, &u) in q.iter().enumerate() { + c[k] = if u == v { at } else { nodes.get(u as usize).copied().ok_or_else(|| Error::Failed("node index".into()))? }; + } + Ok(c) +} + +/// The state of `v`'s quadrilaterals with `v` at some position. +#[derive(Debug, Clone, Copy)] +struct Local { + /// Smallest scaled Jacobian. + sj: f64, + /// Smallest unnormalized corner cross product among the corners that move with `v` (linear + /// in `v`'s position: the untangling objective). + area: f64, + /// All valid by the exact test. + valid: bool, +} + +impl Local { + /// Valid beats invalid; then a valid patch ranks by its worst scaled Jacobian and an invalid + /// one by its worst corner area. + fn key(&self) -> (bool, f64) { + (self.valid, if self.valid { self.sj } else { self.area }) + } + fn better(&self, than: &Local) -> bool { + let (a, b) = (self.key(), than.key()); + (a.0 && !b.0) || (a.0 == b.0 && a.1 > b.1) + } + fn no_worse(&self, than: &Local) -> bool { + !than.better(self) + } +} + +/// `v`'s quadrilaterals with `v` at `at`. +fn local(nodes: &[[f64; 2]], quads: &[[u32; 4]], incident: &[u32], v: u32, at: [f64; 2]) -> Result { + let mut l = Local { sj: f64::INFINITY, area: f64::INFINITY, valid: true }; + for &e in incident { + let quad = quads.get(e as usize).ok_or_else(|| Error::Failed("quad index".into()))?; + let c = corners(nodes, quad, v, at)?; + let k = quad.iter().position(|&u| u == v).unwrap_or(0); + for i in 0..4 { + let (cross, scaled) = planar_corner(&c, i); + l.sj = l.sj.min(scaled.unwrap_or(-1.0)); + // The corner opposite `v` does not move with it. + if i != (k + 2) % 4 { + l.area = l.area.min(cross); + } + } + if l.valid { + let [a, b, c, d] = c.map(|x| Point2::new(x[0], x[1])); + l.valid = turns_left([&a?, &b?, &c?, &d?]); + } + } + Ok(l) +} + +/// Nodes whose quadrilaterals score below this get optimization-based smoothing. +const POOR: f64 = 0.3; +/// Compass-search iterations per node. +const SEARCH_STEPS: usize = 48; + +/// Compass search for the best position of `v` ([`Local::key`]: untangle first, then maximize +/// the smallest scaled Jacobian). Returns the new position when it improves. +fn optimize(nodes: &[[f64; 2]], quads: &[[u32; 4]], incident: &[u32], neighbours: &[u32], v: u32) -> Result> { + let old = nodes.get(v as usize).copied().ok_or_else(|| Error::Failed("node index".into()))?; + let start = local(nodes, quads, incident, v, old)?; + if start.valid && start.sj >= POOR { + return Ok(None); + } + let mut reach = 0.0_f64; + for &u in neighbours { + let p = nodes.get(u as usize).copied().ok_or_else(|| Error::Failed("node index".into()))?; + reach = reach.max(((p[0] - old[0]).powi(2) + (p[1] - old[1]).powi(2)).sqrt()); + } + let (mut at, mut best) = (old, start); + let mut step = 0.25 * reach; + let floor = 1e-4 * step; + const D: f64 = std::f64::consts::FRAC_1_SQRT_2; + const DIRS: [[f64; 2]; 8] = [[1.0, 0.0], [0.0, 1.0], [-1.0, 0.0], [0.0, -1.0], [D, D], [-D, D], [-D, -D], [D, -D]]; + for _ in 0..SEARCH_STEPS { + if !(step > floor) { + break; + } + let mut moved = false; + for d in DIRS { + let to = [snap_to_domain(at[0] + step * d[0]), snap_to_domain(at[1] + step * d[1])]; + let l = local(nodes, quads, incident, v, to)?; + if l.better(&best) { + (at, best) = (to, l); + moved = true; + break; + } + } + if !moved { + step *= 0.5; + } + } + Ok((at != old && best.better(&start)).then_some(at)) +} + +/// Smooth the nodes from `first_free` on for at most `sweeps` sweeps. Returns the sweeps run and +/// whether the largest move of the last sweep fell below `1e-4` of the mean edge length. +pub(crate) fn smooth(nodes: &mut [[f64; 2]], quads: &[[u32; 4]], first_free: usize, sweeps: usize) -> Result<(usize, bool)> { + let n = nodes.len(); + let mut edges = Vec::with_capacity(8 * quads.len()); + let mut inc = Vec::with_capacity(4 * quads.len()); + let mut total = 0.0; + for (e, q) in quads.iter().enumerate() { + for k in 0..4 { + let (a, b) = (q[k], q[(k + 1) % 4]); + edges.push((a, b)); + edges.push((b, a)); + inc.push((a, e as u32)); + if let (Some(pa), Some(pb)) = (nodes.get(a as usize), nodes.get(b as usize)) { + total += ((pa[0] - pb[0]).powi(2) + (pa[1] - pb[1]).powi(2)).sqrt(); + } + } + } + if quads.is_empty() { + return Ok((0, true)); + } + let tolerance = 1e-4 * total / (4 * quads.len()) as f64; + let neighbours = Csr::of(n, &mut edges); + let incident = Csr::of(n, &mut inc); + for sweep in 0..sweeps { + let mut largest = 0.0_f64; + for v in first_free..n { + let nb = neighbours.get(v); + if nb.is_empty() { + continue; + } + let old = nodes.get(v).copied().ok_or_else(|| Error::Failed("node index".into()))?; + let mut avg = [0.0, 0.0]; + for &u in nb { + let p = nodes.get(u as usize).copied().ok_or_else(|| Error::Failed("node index".into()))?; + avg = [avg[0] + p[0], avg[1] + p[1]]; + } + let avg = [avg[0] / nb.len() as f64, avg[1] / nb.len() as f64]; + let quads_of_v = incident.get(v); + let before = local(nodes, quads, quads_of_v, v as u32, old)?; + for step in [1.0, 0.5, 0.25] { + let at = [snap_to_domain(old[0] + step * (avg[0] - old[0])), snap_to_domain(old[1] + step * (avg[1] - old[1]))]; + if at == old { + break; + } + if local(nodes, quads, quads_of_v, v as u32, at)?.no_worse(&before) { + largest = largest.max(((at[0] - old[0]).powi(2) + (at[1] - old[1]).powi(2)).sqrt()); + if let Some(x) = nodes.get_mut(v) { + *x = at; + } + break; + } + } + } + if largest <= tolerance { + let polished = polish(nodes, quads, first_free, &neighbours, &incident, sweeps)?; + return Ok((sweep + 1, polished)); + } + } + polish(nodes, quads, first_free, &neighbours, &incident, sweeps)?; + Ok((sweeps, false)) +} + +/// Optimization-based smoothing of the poor and inverted nodes, at most `sweeps` sweeps. +/// Returns whether a sweep ended with nothing left to improve. +fn polish(nodes: &mut [[f64; 2]], quads: &[[u32; 4]], first_free: usize, neighbours: &Csr, incident: &Csr, sweeps: usize) -> Result { + for _ in 0..sweeps { + let mut changed = false; + for v in first_free..nodes.len() { + if let Some(at) = optimize(nodes, quads, incident.get(v), neighbours.get(v), v as u32)? { + if let Some(x) = nodes.get_mut(v) { + *x = at; + } + changed = true; + } + } + if !changed { + return Ok(true); + } + } + Ok(false) +} diff --git a/crates/mesh/src/quad/tests.rs b/crates/mesh/src/quad/tests.rs new file mode 100644 index 0000000..e9b7942 --- /dev/null +++ b/crates/mesh/src/quad/tests.rs @@ -0,0 +1,539 @@ +use std::collections::BTreeMap; +use std::f64::consts::PI; + +use super::*; +use crate::quality::Validity; +use crate::size::{UniformSize, divide_polyline, even_divisions, ideal_divisions, loop_divisions}; + +fn opts() -> QuadOptions { + QuadOptions::default() +} + +/// A `w × h` rectangle at `(x0, y0)`, counterclockwise, `nx` divisions along x and `ny` along y. +fn rect(x0: f64, y0: f64, w: f64, h: f64, nx: usize, ny: usize) -> Vec<[f64; 2]> { + let mut p = Vec::new(); + for i in 0..nx { + p.push([x0 + w * i as f64 / nx as f64, y0]); + } + for j in 0..ny { + p.push([x0 + w, y0 + h * j as f64 / ny as f64]); + } + for i in 0..nx { + p.push([x0 + w - w * i as f64 / nx as f64, y0 + h]); + } + for j in 0..ny { + p.push([x0, y0 + h - h * j as f64 / ny as f64]); + } + p +} + +fn circle(cx: f64, cy: f64, r: f64, n: usize) -> Vec<[f64; 2]> { + (0..n).map(|k| 2.0 * PI * k as f64 / n as f64).map(|t| [cx + r * t.cos(), cy + r * t.sin()]).collect() +} + +/// A closed polygon with every side divided evenly into `per_side` pieces. +fn subdivide(corners: &[[f64; 2]], per_side: usize) -> Vec<[f64; 2]> { + let mut p = Vec::new(); + for (i, a) in corners.iter().enumerate() { + let b = corners[(i + 1) % corners.len()]; + for k in 0..per_side { + let t = k as f64 / per_side as f64; + p.push([a[0] + (b[0] - a[0]) * t, a[1] + (b[1] - a[1]) * t]); + } + } + p +} + +fn signed_area(p: &[[f64; 2]]) -> f64 { + (0..p.len()).map(|i| p[i][0] * p[(i + 1) % p.len()][1] - p[(i + 1) % p.len()][0] * p[i][1]).sum::() / 2.0 +} + +/// Every mesh invariant: boundary nodes unchanged and first, every element valid by the exact +/// test and counterclockwise, conformity (each interior edge used once in each direction, each +/// boundary edge once), the boundary edges exactly the input loops, and the area. +fn check(mesh: &QuadMesh, loops: &[Vec<[f64; 2]>]) { + let input: Vec<[f64; 2]> = loops.iter().flatten().copied().collect(); + assert_eq!(mesh.boundary_nodes, input.len()); + assert_eq!(&mesh.nodes[..input.len()], input.as_slice()); + let mut area = 0.0; + let mut directed: BTreeMap<(u32, u32), usize> = BTreeMap::new(); + for q in &mesh.quads { + let c = q.map(|v| mesh.nodes[v as usize]); + assert_eq!(crate::quality::planar_quad_validity(&c).unwrap(), Validity::Valid, "{c:?}"); + area += signed_area(&c); + for k in 0..4 { + *directed.entry((q[k], q[(k + 1) % 4])).or_default() += 1; + } + } + assert!(directed.values().all(|&n| n == 1), "an edge is used twice in the same direction"); + let boundary: Vec<(u32, u32)> = directed.keys().filter(|(a, b)| !directed.contains_key(&(*b, *a))).copied().collect(); + let mut expected = Vec::new(); + let mut base = 0u32; + for l in loops { + let n = l.len() as u32; + let ccw = signed_area(l) > 0.0; + let outer = base == 0; + for i in 0..n { + let (a, b) = (base + i, base + (i + 1) % n); + // The outer loop runs counterclockwise around the mesh, holes clockwise. + expected.push(if ccw == outer { (a, b) } else { (b, a) }); + } + base += n; + } + expected.sort_unstable(); + assert_eq!(boundary, expected, "the mesh boundary is not the input loops"); + let want: f64 = loops.iter().enumerate().map(|(k, l)| if k == 0 { signed_area(l).abs() } else { -signed_area(l).abs() }).sum(); + assert!((area - want).abs() <= 1e-9 * want.abs().max(1.0), "area {area} vs {want}"); +} + +#[test] +fn rectangle_free() { + let loops = vec![rect(0.0, 0.0, 10.0, 5.0, 20, 10)]; + let m = free_quad_mesh(&loops, &UniformSize(0.5), &opts()).unwrap(); + check(&m, &loops); + assert!(m.stats.min_scaled_jacobian > 0.3, "{:?}", m.stats); + // About (10 × 5) / 0.5² = 200 quads. + assert!((100..400).contains(&m.quads.len()), "{}", m.quads.len()); +} + +#[test] +fn single_quad_and_a_six_point_loop() { + let quad = vec![vec![[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]]; + let m = free_quad_mesh(&quad, &UniformSize(1.0), &opts()).unwrap(); + assert_eq!(m.quads.len(), 1); + check(&m, &quad); + let six = vec![rect(0.0, 0.0, 2.0, 1.0, 2, 1)]; + let m = free_quad_mesh(&six, &UniformSize(1.0), &opts()).unwrap(); + check(&m, &six); +} + +#[test] +fn clockwise_input_is_accepted() { + let mut l = rect(0.0, 0.0, 4.0, 2.0, 8, 4); + l.reverse(); + let loops = vec![l]; + let m = free_quad_mesh(&loops, &UniformSize(0.5), &opts()).unwrap(); + check(&m, &loops); +} + +#[test] +fn l_shape() { + let corners: [[f64; 2]; 6] = [[0.0, 0.0], [10.0, 0.0], [10.0, 4.0], [4.0, 4.0], [4.0, 10.0], [0.0, 10.0]]; + // Even divisions per side keep the corners at even positions. + let mut l = Vec::new(); + for (i, a) in corners.iter().enumerate() { + let b = corners[(i + 1) % corners.len()]; + let len = ((b[0] - a[0]).powi(2) + (b[1] - a[1]).powi(2)).sqrt(); + let n = even_divisions(len / 0.5).unwrap(); + for k in 0..n { + let t = k as f64 / n as f64; + l.push([a[0] + (b[0] - a[0]) * t, a[1] + (b[1] - a[1]) * t]); + } + } + let loops = vec![l]; + let m = free_quad_mesh(&loops, &UniformSize(0.5), &opts()).unwrap(); + check(&m, &loops); + assert!(m.stats.min_scaled_jacobian > 0.2, "{:?}", m.stats); +} + +#[test] +fn plate_with_a_hole() { + let loops = vec![rect(0.0, 0.0, 20.0, 10.0, 40, 20), circle(10.0, 5.0, 2.5, 32)]; + let m = free_quad_mesh(&loops, &UniformSize(0.5), &opts()).unwrap(); + check(&m, &loops); + assert!(m.stats.min_scaled_jacobian > 0.2, "{:?}", m.stats); +} + +#[test] +fn plate_with_many_holes() { + let mut loops = vec![rect(0.0, 0.0, 30.0, 30.0, 30, 30)]; + for i in 0..3 { + for j in 0..3 { + loops.push(circle(5.0 + 10.0 * i as f64, 5.0 + 10.0 * j as f64, 2.0 + 0.5 * ((i + j) % 2) as f64, 16 + 4 * i)); + } + } + // A square hole too, given clockwise. + loops[5] = rect(13.0, 13.0, 4.0, 4.0, 4, 4).into_iter().rev().collect(); + let m = free_quad_mesh(&loops, &UniformSize(1.0), &opts()).unwrap(); + check(&m, &loops); +} + +#[test] +fn thin_sliver() { + let loops = vec![rect(0.0, 0.0, 100.0, 1.0, 200, 2)]; + let m = free_quad_mesh(&loops, &UniformSize(0.5), &opts()).unwrap(); + check(&m, &loops); + assert_eq!(m.quads.len(), 400); + // A sliver one division across: the split puts a mid node on each short side. + let loops = vec![rect(0.0, 0.0, 50.0, 0.25, 100, 1)]; + let m = free_quad_mesh(&loops, &UniformSize(0.5), &opts()).unwrap(); + check(&m, &loops); +} + +fn star(r_out: f64, r_in: f64) -> Vec<[f64; 2]> { + (0..10).map(|k| (k as f64 * PI / 5.0 + PI / 2.0, if k % 2 == 0 { r_out } else { r_in })).map(|(t, r)| [r * t.cos(), r * t.sin()]).collect() +} + +#[test] +fn concave_star() { + let loops = vec![subdivide(&star(10.0, 4.0), 4)]; + let m = free_quad_mesh(&loops, &UniformSize(1.0), &opts()).unwrap(); + check(&m, &loops); + // A star with a star-shaped hole. + let loops = vec![subdivide(&star(10.0, 4.0), 4), subdivide(&star(2.5, 1.0), 2)]; + let m = free_quad_mesh(&loops, &UniformSize(1.0), &opts()).unwrap(); + check(&m, &loops); +} + +/// A polygon with every side divided evenly at about `h` (even counts, so corners stay at even +/// positions). +fn divided(corners: &[[f64; 2]], h: f64) -> Vec<[f64; 2]> { + let mut l = Vec::new(); + for (i, a) in corners.iter().enumerate() { + let b = corners[(i + 1) % corners.len()]; + let n = even_divisions(((b[0] - a[0]).powi(2) + (b[1] - a[1]).powi(2)).sqrt() / h).unwrap(); + l.extend((0..n).map(|k| k as f64 / n as f64).map(|t| [a[0] + (b[0] - a[0]) * t, a[1] + (b[1] - a[1]) * t])); + } + l +} + +#[test] +fn i_beam_end_face() { + // Flanges 10 × 1, web 1 thick, 12 deep overall: thin walls meeting at reflex corners. + let c = [ + [0.0, 0.0], + [10.0, 0.0], + [10.0, 1.0], + [5.5, 1.0], + [5.5, 11.0], + [10.0, 11.0], + [10.0, 12.0], + [0.0, 12.0], + [0.0, 11.0], + [4.5, 11.0], + [4.5, 1.0], + [0.0, 1.0], + ]; + for h in [0.5, 0.25] { + let loops = vec![divided(&c, h)]; + let m = free_quad_mesh(&loops, &UniformSize(h), &opts()).unwrap(); + check(&m, &loops); + assert!(m.stats.min_scaled_jacobian > 0.2, "{h}: {:?}", m.stats); + } +} + +#[test] +fn perforated_plate() { + // A 40 × 20 plate with a 7 × 3 grid of holes of radius 1 on a 5 mm pitch. + let h = 0.5; + let mut loops = vec![divided(&[[0.0, 0.0], [40.0, 0.0], [40.0, 20.0], [0.0, 20.0]], h)]; + for i in 0..7 { + for j in 0..3 { + loops.push(circle(5.0 + 5.0 * i as f64, 5.0 + 5.0 * j as f64, 1.0, 2 * (PI / h).ceil() as usize)); + } + } + let m = free_quad_mesh(&loops, &UniformSize(h), &opts()).unwrap(); + check(&m, &loops); + assert!(m.stats.min_scaled_jacobian > 0.2, "{:?}", m.stats); +} + +#[test] +fn comb_with_deep_slots() { + // Teeth of width 1 between slots of width 1, 4 deep: strongly concave. + let mut corners = vec![[0.0, 0.0], [9.0, 0.0]]; + for k in (0..5).rev() { + let x = 2.0 * k as f64; + corners.extend([[x + 1.0, 6.0], [x, 6.0]]); + if k > 0 { + corners.extend([[x, 2.0], [x - 1.0, 2.0]]); + } + } + let loops = vec![subdivide(&corners, 4)]; + let m = free_quad_mesh(&loops, &UniformSize(0.5), &opts()).unwrap(); + check(&m, &loops); +} + +#[test] +fn graded_size_field() { + let loops = vec![rect(0.0, 0.0, 20.0, 20.0, 20, 20)]; + let near = |p: [f64; 3]| 0.25 + 0.15 * (p[0] * p[0] + p[1] * p[1]).sqrt(); + let graded = free_quad_mesh(&loops, &near, &opts()).unwrap(); + check(&graded, &loops); + let uniform = free_quad_mesh(&loops, &UniformSize(1.0), &opts()).unwrap(); + // Denser than the uniform mesh near the origin, coarser far from it. + let within = |m: &QuadMesh, r: f64| m.nodes.iter().filter(|p| p[0] * p[0] + p[1] * p[1] < r * r).count(); + assert!(within(&graded, 8.0) > within(&uniform, 8.0), "{} {}", within(&graded, 8.0), within(&uniform, 8.0)); + assert!(graded.nodes.len() - within(&graded, 16.0) < uniform.nodes.len() - within(&uniform, 16.0)); +} + +#[test] +fn boundary_from_the_size_helpers() { + // An arc-and-line boundary whose divisions come from the size field, loop parity fixed. + let field = |p: [f64; 3]| 0.3 + 0.05 * p[0].abs(); + let arc: Vec<[f64; 3]> = (0..=32).map(|k| PI * k as f64 / 32.0).map(|t| [5.0 * t.cos(), 5.0 * t.sin(), 0.0]).collect(); + let base = vec![[-5.0, 0.0, 0.0], [5.0, 0.0, 0.0]]; + let curves = [arc, base]; + let ideal: Vec = curves.iter().map(|c| ideal_divisions(&field, c).unwrap()).collect(); + let counts = loop_divisions(&ideal).unwrap(); + let mut l = Vec::new(); + for (c, n) in curves.iter().zip(counts) { + let pts = divide_polyline(&field, c, n).unwrap(); + l.extend(pts[..pts.len() - 1].iter().map(|p| [p[0], p[1]])); + } + assert_eq!(l.len() % 2, 0); + let loops = vec![l]; + let m = free_quad_mesh(&loops, &field, &opts()).unwrap(); + check(&m, &loops); +} + +#[test] +fn deterministic() { + let loops = vec![rect(0.0, 0.0, 20.0, 10.0, 40, 20), circle(7.0, 5.0, 2.0, 24), circle(14.0, 5.0, 3.0, 30)]; + let a = free_quad_mesh(&loops, &UniformSize(0.6), &opts()).unwrap(); + let b = free_quad_mesh(&loops, &UniformSize(0.6), &opts()).unwrap(); + assert_eq!(a, b); +} + +#[test] +fn no_smoothing_is_still_valid() { + let loops = vec![rect(0.0, 0.0, 20.0, 10.0, 40, 20), circle(10.0, 5.0, 2.5, 32)]; + let m = free_quad_mesh(&loops, &UniformSize(0.5), &QuadOptions { smoothing_sweeps: 0, ..opts() }).unwrap(); + check(&m, &loops); + assert_eq!(m.stats.smoothing_sweeps, 0); +} + +#[test] +fn hostile_free_input() { + let s = UniformSize(1.0); + let ok = rect(0.0, 0.0, 4.0, 4.0, 4, 4); + let invalid = |loops: Vec>| matches!(free_quad_mesh(&loops, &s, &opts()), Err(Error::Invalid(_))); + assert!(invalid(vec![])); + assert!(invalid(vec![vec![[0.0, 0.0], [1.0, 0.0]]])); + // Odd count: an error, never a silent repair. + let mut odd = ok.clone(); + odd.push([0.0, 0.5]); + assert!(invalid(vec![odd])); + let mut nan = ok.clone(); + nan[3] = [f64::NAN, 1.0]; + assert!(invalid(vec![nan])); + let mut inf = ok.clone(); + inf[3] = [f64::INFINITY, 1.0]; + assert!(invalid(vec![inf])); + let mut far = ok.clone(); + far[3] = [1e300, 1.0]; + assert!(invalid(vec![far])); + // Zero-length edge and repeated point. + let mut zero = ok.clone(); + zero[2] = zero[1]; + assert!(invalid(vec![zero])); + let mut repeated = ok.clone(); + repeated[10] = repeated[2]; + assert!(invalid(vec![repeated])); + // Self-intersecting (bow tie) loop. + assert!(invalid(vec![vec![[0.0, 0.0], [2.0, 2.0], [2.0, 1.0], [2.0, 0.0], [0.0, 2.0], [0.0, 1.0]]])); + // A point touching another edge. + assert!(invalid(vec![vec![[0.0, 0.0], [4.0, 0.0], [4.0, 4.0], [2.0, 0.0], [0.0, 4.0], [0.0, 2.0]]])); + // Collinear fold-back. + assert!(invalid(vec![vec![[0.0, 0.0], [4.0, 0.0], [2.0, 0.0], [2.0, 2.0], [1.0, 2.0], [0.0, 2.0]]])); + // All collinear. + assert!(invalid(vec![vec![[0.0, 0.0], [1.0, 0.0], [2.0, 0.0], [3.0, 0.0], [4.0, 0.0], [5.0, 0.0]]])); + // Holes: outside, crossing the outer loop, nested, too few points. + assert!(invalid(vec![ok.clone(), circle(10.0, 10.0, 1.0, 8)])); + assert!(invalid(vec![ok.clone(), circle(4.0, 2.0, 1.0, 8)])); + assert!(invalid(vec![rect(0.0, 0.0, 10.0, 10.0, 10, 10), circle(5.0, 5.0, 3.0, 16), circle(5.0, 5.0, 1.0, 8)])); + assert!(invalid(vec![ok.clone(), rect(1.0, 1.0, 2.0, 2.0, 1, 1)])); + // Two holes crossing each other. + assert!(invalid(vec![rect(0.0, 0.0, 10.0, 10.0, 10, 10), circle(4.0, 5.0, 2.0, 16), circle(6.0, 5.0, 2.0, 16)])); + // Bad size fields. + for h in [0.0, -1.0, f64::NAN, f64::INFINITY] { + assert!(matches!(free_quad_mesh(std::slice::from_ref(&ok), &UniformSize(h), &opts()), Err(Error::Invalid(_))), "{h}"); + } + // Limits. + assert!(matches!(free_quad_mesh(std::slice::from_ref(&ok), &UniformSize(1e-4), &opts()), Err(Error::Limit(_)))); + assert!(matches!(free_quad_mesh(&[ok], &UniformSize(0.1), &QuadOptions { max_elements: 10, ..opts() }), Err(Error::Limit(_)))); +} + +#[test] +fn random_boundaries_never_panic() { + // Pseudo-random jagged loops (a fixed linear congruential sequence), every point a corner: + // every result is a valid mesh or an input error. + let mut seed = 12345u64; + let mut next = || { + seed = seed.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); + ((seed >> 33) as f64) / (1u64 << 31) as f64 + }; + let (mut meshed, mut rejected) = (0, 0); + for case in 0..200 { + let n = 6 + 2 * (case % 10); + let mut loops = vec![ + (0..n) + .map(|k| { + let t = 2.0 * PI * k as f64 / n as f64; + let r = 1.0 + 3.0 * next(); + [r * t.cos(), r * t.sin()] + }) + .collect::>(), + ]; + if case % 3 == 0 { + loops.push(circle(0.0, 0.0, 0.5 * next(), 6)); + } + if case % 7 == 0 { + loops[0] = (0..n).map(|_| [10.0 * next() - 5.0, 10.0 * next() - 5.0]).collect(); + } + match free_quad_mesh(&loops, &UniformSize(0.2 + next()), &opts()) { + Ok(m) => { + check(&m, &loops); + meshed += 1; + } + Err(Error::Invalid(_)) => rejected += 1, + Err(e) => panic!("case {case}: {e}"), + } + } + // Every valid region is meshed; the rest are self-intersecting loops. + assert_eq!(meshed + rejected, 200); + assert!(meshed > 150, "{meshed}"); +} + +#[test] +fn mapped_rectangle_is_the_grid() { + let l = rect(0.0, 0.0, 4.0, 2.0, 4, 2); + let m = mapped_quad_mesh(&l, [0, 4, 6, 10], &opts()).unwrap(); + assert_eq!((m.nu, m.nv), (4, 2)); + assert_eq!(m.mesh.quads.len(), 8); + check(&m.mesh, std::slice::from_ref(&l)); + for j in 0..=2 { + for i in 0..=4 { + assert_eq!(m.mesh.nodes[m.grid[j * 5 + i] as usize], [i as f64, j as f64]); + } + } + assert!((m.mesh.stats.min_scaled_jacobian - 1.0).abs() < 1e-12); + assert!(m.mesh.stats.max_skew < 1e-12); +} + +#[test] +fn mapped_curved_and_graded() { + // A quarter annulus, radii 1 and 2, graded radially. + let (na, nr) = (8, 5); + let r = |j: usize| 1.0 + (j as f64 / nr as f64).powi(2); + let mut l = Vec::new(); + for j in 0..nr { + l.push([r(j), 0.0]); + } + for i in 0..na { + let t = 0.5 * PI * i as f64 / na as f64; + l.push([2.0 * t.cos(), 2.0 * t.sin()]); + } + for j in 0..nr { + l.push([0.0, r(nr - j)]); + } + for i in 0..na { + let t = 0.5 * PI * (na - i) as f64 / na as f64; + l.push([t.cos(), t.sin()]); + } + let corners = [0, nr, nr + na, 2 * nr + na]; + let m = mapped_quad_mesh(&l, corners, &opts()).unwrap(); + check(&m.mesh, std::slice::from_ref(&l)); + // The same region given clockwise. + let mut cw = l.clone(); + cw.reverse(); + let n = cw.len(); + let m = mapped_quad_mesh(&cw, [n - 1 - corners[3], n - 1 - corners[2], n - 1 - corners[1], n - 1 - corners[0]], &opts()).unwrap(); + check(&m.mesh, &[cw]); +} + +#[test] +fn mapped_hostile() { + let l = rect(0.0, 0.0, 4.0, 2.0, 4, 2); + assert!(mapped_quad_mesh(&l, [0, 4, 6, 11], &opts()).is_err()); // unequal opposite sides + assert!(mapped_quad_mesh(&l, [0, 6, 4, 10], &opts()).is_err()); // out of order + assert!(mapped_quad_mesh(&l, [0, 4, 6, 99], &opts()).is_err()); + assert!(mapped_quad_mesh(&l, [0, 0, 6, 10], &opts()).is_err()); + let mut nan = l.clone(); + nan[5] = [f64::NAN, 0.0]; + assert!(mapped_quad_mesh(&nan, [0, 4, 6, 10], &opts()).is_err()); + // A region whose TFI folds: the top side dips almost to the bottom between its corners. + let fold = vec![[0.0, 0.0], [1.0, 0.0], [2.0, 0.0], [3.0, 0.0], [4.0, 0.0], [4.0, 4.0], [3.0, 4.0], [2.2, 0.2], [1.8, 0.2], [1.0, 4.0]]; + let r = mapped_quad_mesh(&fold, [0, 4, 5, 9], &opts()); + assert!(matches!(r, Err(Error::Failed(_)) | Err(Error::Invalid(_))), "{r:?}"); +} + +#[test] +fn transfinite_in_space() { + let bottom = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [2.0, 0.0, 0.0]]; + let top = [[0.0, 1.0, 1.0], [1.0, 1.0, 1.0], [2.0, 1.0, 1.0]]; + let left = [[0.0, 0.0, 0.0], [0.0, 1.0, 1.0]]; + let right = [[2.0, 0.0, 0.0], [2.0, 1.0, 1.0]]; + let g = transfinite_grid(&bottom, &right, &top, &left).unwrap(); + assert_eq!(g.len(), 6); + assert_eq!(g[4], [1.0, 1.0, 1.0]); + let side = |a: [f64; 2], b: [f64; 2]| [a, [(a[0] + b[0]) / 2.0, (a[1] + b[1]) / 2.0], b]; + let g = + transfinite_grid(&side([0.0, 0.0], [2.0, 0.0]), &side([2.0, 0.0], [2.0, 2.0]), &side([0.0, 2.0], [2.0, 2.0]), &side([0.0, 0.0], [0.0, 2.0])) + .unwrap(); + assert_eq!(g[4], [1.0, 1.0]); + // Corners that do not match, a zero-length side, empty sides. + assert!(transfinite_grid(&bottom, &right, &top, &[[0.0, 0.0, 0.1], [0.0, 1.0, 1.0]]).is_err()); + assert!(transfinite_grid(&[[0.0; 2], [0.0; 2]], &[[0.0; 2], [0.0, 1.0]], &[[0.0, 1.0], [0.0, 1.0]], &[[0.0; 2], [0.0, 1.0]]).is_err()); + assert!(transfinite_grid::<2>(&[], &[], &[], &[]).is_err()); +} + +#[test] +fn planar_quad_metrics() { + use crate::quality::{Metric, planar_quad_quality, planar_quad_validity}; + let rect = [[0.0, 0.0], [2.0, 0.0], [2.0, 1.0], [0.0, 1.0]]; + let q = planar_quad_quality(&rect).unwrap(); + assert_eq!(q.validity(), Some(Validity::Valid)); + assert_eq!(q.size, 2.0); + assert_eq!(q.get(Metric::AspectRatio), Some(2.0)); + assert_eq!(q.get(Metric::Skew), Some(0.0)); + assert_eq!(q.get(Metric::ScaledJacobian), Some(1.0)); + assert_eq!(q.get(Metric::JacobianRatio), Some(1.0)); + // A rhombus with 60° corners: its principal axes are its edge directions, so skew cos 60° + // and aspect 1; the scaled Jacobian is sin 60°. + let h = 3.0_f64.sqrt() / 2.0; + let rh = [[0.0, 0.0], [1.0, 0.0], [1.5, h], [0.5, h]]; + let q = planar_quad_quality(&rh).unwrap(); + assert!((q.get(Metric::ScaledJacobian).unwrap() - h).abs() < 1e-12); + assert!((q.get(Metric::Skew).unwrap() - 0.5).abs() < 1e-12); + assert!((q.get(Metric::AspectRatio).unwrap() - 1.0).abs() < 1e-12); + // Clockwise, bow tie, concave (dart), degenerate (a straight corner): invalid by the exact test. + let cw = [rect[0], rect[3], rect[2], rect[1]]; + assert_eq!(planar_quad_validity(&cw).unwrap(), Validity::Invalid); + assert_eq!(planar_quad_validity(&[[0.0, 0.0], [1.0, 1.0], [1.0, 0.0], [0.0, 1.0]]).unwrap(), Validity::Invalid); + assert_eq!(planar_quad_validity(&[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.7, 0.3]]).unwrap(), Validity::Invalid); + assert_eq!(planar_quad_validity(&[[0.0, 0.0], [1.0, 0.0], [2.0, 0.0], [0.0, 1.0]]).unwrap(), Validity::Invalid); + // A corner off by one ulp from straight is decided exactly. + let x = 1.0 + f64::EPSILON; + assert_eq!(planar_quad_validity(&[[0.0, 0.0], [1.0, 0.0], [2.0, x - 1.0], [0.0, 1.0]]).unwrap(), Validity::Valid); + assert!(planar_quad_validity(&[[0.0, 0.0], [1.0, f64::NAN], [1.0, 1.0], [0.0, 1.0]]).is_err()); + let q = planar_quad_quality(&[[0.0, 0.0], [1.0, 0.0], [2.0, 0.0], [0.0, 1.0]]).unwrap(); + assert_eq!(q.get(Metric::JacobianRatio), Some(0.0)); +} + +#[test] +#[ignore = "timing and quality report: cargo test --release -p modelcraft-mesh quad_report -- --ignored --nocapture"] +fn quad_report() { + let cases: Vec<(&str, Vec>, f64)> = vec![ + ("rectangle", vec![rect(0.0, 0.0, 10.0, 5.0, 20, 10)], 0.5), + ("plate with hole", vec![rect(0.0, 0.0, 20.0, 10.0, 40, 20), circle(10.0, 5.0, 2.5, 32)], 0.5), + ("star", vec![subdivide(&star(10.0, 4.0), 4)], 1.0), + ("large plate", vec![rect(0.0, 0.0, 200.0, 100.0, 1000, 500), circle(50.0, 50.0, 20.0, 400), circle(150.0, 50.0, 20.0, 400)], 0.2), + ]; + for (name, loops, h) in cases { + let t = std::time::Instant::now(); + let m = free_quad_mesh(&loops, &UniformSize(h), &opts()).unwrap(); + println!("{name}: {} quads in {:?}, {:?}", m.quads.len(), t.elapsed(), m.stats); + let t = std::time::Instant::now(); + let _ = free_quad_mesh(&loops, &UniformSize(h), &QuadOptions { smoothing_sweeps: 0, ..opts() }).unwrap(); + println!(" without smoothing {:?}", t.elapsed()); + } +} + +#[test] +fn unsmoothed_meshes_are_valid() { + // Regression: recombination once accepted nearly straight corners that the split turned + // degenerate after rounding the mid nodes (this plate had one before the quality floor). + let loops = vec![rect(0.0, 0.0, 60.0, 30.0, 300, 150), circle(15.0, 15.0, 6.0, 240)]; + let m = free_quad_mesh(&loops, &UniformSize(0.2), &QuadOptions { smoothing_sweeps: 0, ..opts() }).unwrap(); + check(&m, &loops); +} diff --git a/crates/mesh/src/quad/tfi.rs b/crates/mesh/src/quad/tfi.rs new file mode 100644 index 0000000..0a828f8 --- /dev/null +++ b/crates/mesh/src/quad/tfi.rs @@ -0,0 +1,108 @@ +//! Mapped meshing of four-sided regions by transfinite interpolation. +//! +//! W. J. Gordon, C. A. Hall, *Construction of curvilinear co-ordinate systems and applications to +//! mesh generation*, IJNME 7(4):461–477, 1973: the bilinearly blended Coons map +//! `X(u, v) = (1−v) B + v T + (1−u) L + u R − [(1−u)(1−v) P00 + u(1−v) P10 + u v P11 + (1−u) v P01]` +//! from the four sides `B` (bottom), `R` (right), `T` (top), `L` (left) and their corners. The +//! discrete form uses the side nodes directly; the blending parameters of node `(i, j)` follow the +//! arc-length distribution of the sides, from the intersection of the straight parameter lines +//! `u = (1−v) u_B(i) + v u_T(i)` and `v = (1−u) v_L(j) + u v_R(j)` (the standard treatment of +//! unequal spacing on opposite sides, e.g. Thompson, Soni, Weatherill, *Handbook of Grid +//! Generation*, CRC 1999, ch. 3). The map is dimension-generic: faces of a block in space use it +//! with `D = 3`. + +use crate::{Error, Result}; + +/// The most nodes of one mapped grid. +pub const MAX_GRID_NODES: usize = 1 << 24; + +fn dist(a: &[f64; D], b: &[f64; D]) -> f64 { + a.iter().zip(b).map(|(x, y)| (x - y) * (x - y)).sum::().sqrt() +} + +/// Normalized cumulative arc length of a side, from 0 to 1. +fn params(side: &[[f64; D]], name: &str) -> Result> { + let mut acc = Vec::with_capacity(side.len()); + let mut s = 0.0; + acc.push(0.0); + for w in side.windows(2) { + if let [a, b] = w { + s += dist(a, b); + } + acc.push(s); + } + if !(s > 0.0) || !s.is_finite() { + return Err(Error::Invalid(format!("the {name} side has zero or non-finite length"))); + } + let mut out: Vec = acc.iter().map(|x| x / s).collect(); + if let Some(l) = out.last_mut() { + *l = 1.0; + } + Ok(out) +} + +/// The transfinite grid of a four-sided region: `(nu + 1) × (nv + 1)` points, node `(i, j)` at +/// index `j·(nu + 1) + i`. `bottom` and `top` run in `+u` (each `nu + 1` points), `left` and +/// `right` in `+v` (each `nv + 1` points); the corners must match exactly (`bottom[0] = +/// left[0]`, `bottom[nu] = right[0]`, `top[0] = left[nv]`, `top[nu] = right[nv]`). Side nodes are +/// copied exactly into the grid. +pub fn transfinite_grid(bottom: &[[f64; D]], right: &[[f64; D]], top: &[[f64; D]], left: &[[f64; D]]) -> Result> { + let (nu, nv) = (bottom.len().saturating_sub(1), left.len().saturating_sub(1)); + if nu == 0 || nv == 0 || top.len() != bottom.len() || right.len() != left.len() { + return Err(Error::Invalid(format!( + "opposite sides need the same number of divisions (at least 1): bottom {nu}, top {}, left {nv}, right {}", + top.len().saturating_sub(1), + right.len().saturating_sub(1) + ))); + } + let count = (nu + 1).checked_mul(nv + 1).filter(|&c| c <= MAX_GRID_NODES); + let count = count.ok_or_else(|| Error::Limit(format!("a {nu} × {nv} grid exceeds {MAX_GRID_NODES} nodes")))?; + for side in [bottom, right, top, left] { + modelcraft_model::check_coordinates(side.as_flattened())?; + } + let corner = |s: &[[f64; D]], last: bool| if last { s.last().copied() } else { s.first().copied() }; + let (p00, p10, p01, p11) = match (corner(bottom, false), corner(bottom, true), corner(top, false), corner(top, true)) { + (Some(a), Some(b), Some(c), Some(d)) => (a, b, c, d), + _ => return Err(Error::Invalid("empty side".into())), + }; + if corner(left, false) != Some(p00) || corner(right, false) != Some(p10) || corner(left, true) != Some(p01) || corner(right, true) != Some(p11) { + return Err(Error::Invalid("the sides of a mapped region must share their corner points".into())); + } + let (ub, ut, vl, vr) = (params(bottom, "bottom")?, params(top, "top")?, params(left, "left")?, params(right, "right")?); + let mut grid = Vec::with_capacity(count); + for j in 0..=nv { + for i in 0..=nu { + let (Some(&b), Some(&t), Some(&l), Some(&r)) = (bottom.get(i), top.get(i), left.get(j), right.get(j)) else { + return Err(Error::Failed("grid index".into())); + }; + let point = if j == 0 { + b + } else if j == nv { + t + } else if i == 0 { + l + } else if i == nu { + r + } else { + let (u0, u1, v0, v1) = ( + ub.get(i).copied().unwrap_or(0.0), + ut.get(i).copied().unwrap_or(0.0), + vl.get(j).copied().unwrap_or(0.0), + vr.get(j).copied().unwrap_or(0.0), + ); + let (du, dv) = (u1 - u0, v1 - v0); + // |du|, |dv| < 1 at interior nodes, so the denominator is positive. + let u = (u0 + v0 * du) / (1.0 - du * dv); + let v = v0 + u * dv; + let mut x = [0.0; D]; + for (k, xk) in x.iter_mut().enumerate() { + *xk = (1.0 - v) * b[k] + v * t[k] + (1.0 - u) * l[k] + u * r[k] + - ((1.0 - u) * (1.0 - v) * p00[k] + u * (1.0 - v) * p10[k] + u * v * p11[k] + (1.0 - u) * v * p01[k]); + } + x + }; + grid.push(point); + } + } + Ok(grid) +} diff --git a/crates/mesh/src/quality/mod.rs b/crates/mesh/src/quality/mod.rs index 35a2431..18ceecf 100644 --- a/crates/mesh/src/quality/mod.rs +++ b/crates/mesh/src/quality/mod.rs @@ -27,8 +27,9 @@ use modelcraft_model::{ElementKind, check_coordinates}; pub use jacobian::{JacobianCheck, Validity}; pub use report::{GateOptions, GateReport, MetricStats, QualityReport, enclosed_volume, quality_report, validity_gate, validity_gate_elements}; +pub(crate) use shape::{planar_corner, quad_axes}; -use crate::predicates::{Point3, Sign, tet_volume_sign}; +use crate::predicates::{Point2, Point3, Sign, orient2d, snap_to_domain, tet_volume_sign}; use crate::{Error, Result}; /// A quality metric. Angles are in degrees. @@ -60,6 +61,9 @@ pub enum Metric { MinAngle, /// Largest corner angle of a triangle or quadrilateral. MaxAngle, + /// Skew of a quadrilateral: `|cos θ|` of the angle between its principal axes, in `[0, 1]`; + /// 0 for a rectangle, 1 when degenerate. + Skew, } impl Metric { @@ -76,19 +80,20 @@ impl Metric { Metric::EdgeRatio => "edge_ratio", Metric::MinAngle => "min_angle", Metric::MaxAngle => "max_angle", + Metric::Skew => "skew", } } /// Whether a larger value is a better element. pub fn higher_is_better(self) -> bool { - !matches!(self, Metric::MaxDihedral | Metric::AspectRatio | Metric::EdgeRatio | Metric::MaxAngle) + !matches!(self, Metric::MaxDihedral | Metric::AspectRatio | Metric::EdgeRatio | Metric::MaxAngle | Metric::Skew) } /// The fixed histogram range; values outside it fall in the first or last bin. pub fn range(self) -> (f64, f64) { match self { Metric::ScaledJacobian | Metric::JacobianRatio | Metric::RadiusRatio | Metric::VolumeLength => (-1.0, 1.0), - Metric::Shape => (0.0, 1.0), + Metric::Shape | Metric::Skew => (0.0, 1.0), Metric::MinDihedral | Metric::MaxDihedral | Metric::MinAngle | Metric::MaxAngle => (0.0, 180.0), Metric::AspectRatio | Metric::EdgeRatio => (1.0, 11.0), } @@ -218,5 +223,52 @@ pub fn element_quality(kind: ElementKind, xyz: &[[f64; 3]]) -> Result 0`, decided by the exact predicate. A zero corner (degenerate) or a negative one +/// (concave, bow-tie or clockwise) is [`Validity::Invalid`]. Coordinates below the predicates' +/// smallest magnitude are rounded to zero first. +pub fn planar_quad_validity(p: &[[f64; 2]; 4]) -> Result { + check_coordinates(p.as_flattened())?; + let [a, b, c, d] = p.map(|c| Point2::new(snap_to_domain(c[0]), snap_to_domain(c[1]))); + Ok(if turns_left([&a?, &b?, &c?, &d?]) { Validity::Valid } else { Validity::Invalid }) +} + +/// The exact corner test of [`planar_quad_validity`] on admissible points: every corner of the +/// quadrilateral turns left. +pub(crate) fn turns_left(p: [&Point2; 4]) -> bool { + (0..4).all(|i| orient2d(p[(i + 3) % 4], p[i], p[(i + 1) % 4]) == Sign::Positive) +} + +/// Every metric of a planar quadrilateral (counterclockwise corners in the xy plane) with its +/// exact validity ([`planar_quad_validity`]): the Jacobian ratio is the smallest corner +/// determinant over the largest magnitude (exact bounds, see there) and the size is the signed +/// area. The shape metrics are those of [`element_quality`] for a `Quad4`. +pub fn planar_quad_quality(p: &[[f64; 2]; 4]) -> Result { + let validity = planar_quad_validity(p)?; + let o = p[0]; + let l = p.map(|c| [c[0] - o[0], c[1] - o[1]]); + let corners = [0, 1, 2, 3].map(|i| shape::planar_corner(&l, i)); + let dets = corners.map(|c| c.0); + let lo = dets.iter().fold(f64::INFINITY, |m, &d| m.min(d)); + let hi = dets.iter().fold(0.0_f64, |m, &d| m.max(d.abs())); + let ratio = if hi > 0.0 { (lo / hi).clamp(-1.0, 1.0) } else { 0.0 }; + let area = 0.5 * ((l[2][0] - l[0][0]) * (l[3][1] - l[1][1]) - (l[2][1] - l[0][1]) * (l[3][0] - l[1][0])); + let lifted = l.map(|c| [c[0], c[1], 0.0]); + let (_, mut metrics) = shape::quad(lifted); + // Corner scaled Jacobians measured against +z (the 3D form measures against the diagonals' + // normal, which flips with a clockwise quad). + let scaled = corners.iter().fold(f64::INFINITY, |m, c| m.min(c.1.unwrap_or(0.0))); + for (m, v) in metrics.iter_mut() { + if *m == Metric::ScaledJacobian { + *v = scaled; + } + } + metrics.push((Metric::JacobianRatio, ratio)); + Ok(ElementQuality { kind: ElementKind::Quad4, size: area, jacobian: Some(JacobianCheck { validity, ratio, volume: area }), metrics }) +} + #[cfg(test)] mod tests; diff --git a/crates/mesh/src/quality/shape.rs b/crates/mesh/src/quality/shape.rs index b0d94c7..5c92cc9 100644 --- a/crates/mesh/src/quality/shape.rs +++ b/crates/mesh/src/quality/shape.rs @@ -6,6 +6,10 @@ //! - Hexahedra, wedges and quadrilaterals: Knupp's algebraic metrics built on the corner Jacobian //! matrices (Knupp, *Algebraic mesh quality metrics*, SIAM J. Sci. Comput. 23(1), 2001): the //! scaled Jacobian `det(A) / (|a1| |a2| |a3|)` and the shape metric `3 det(A)^{2/3} / |A|²_F`. +//! - Quadrilateral aspect ratio and skew from the principal axes `X1 = (p1 − p0) + (p2 − p3)`, +//! `X2 = (p2 − p1) + (p3 − p0)` (Robinson, *CRE method of element testing and the Jacobian +//! shape parameters*, Eng. Comput. 4(2):113–118, 1987): aspect `max(|X1|, |X2|) / min(…)`, skew +//! `|X1 · X2| / (|X1| |X2|)`. //! //! Every metric is defined for degenerate input (no division by zero, no NaN): a degenerate //! element scores 0 on bounded metrics and `+∞` on aspect and edge ratios. @@ -132,6 +136,7 @@ pub(crate) fn quad(p: [P; 4]) -> (f64, Vec<(Metric, f64)>) { angles = (angles.0.min(a), angles.1.max(a)); } let edges = [(0, 1), (1, 2), (2, 3), (3, 0)]; + let (aspect, skew) = quad_axes(p); ( area, vec![ @@ -139,10 +144,36 @@ pub(crate) fn quad(p: [P; 4]) -> (f64, Vec<(Metric, f64)>) { (Metric::MinAngle, angles.0), (Metric::MaxAngle, angles.1), (Metric::EdgeRatio, edge_ratio(&p, &edges)), + (Metric::AspectRatio, aspect), + (Metric::Skew, skew), ], ) } +/// Corner `i` of a planar quadrilateral: the cross product `(p[i+1] − p[i]) × (p[i−1] − p[i])` +/// (twice the corner triangle's signed area, positive for a left turn of a counterclockwise quad) +/// and that cross over the two edge lengths (the corner's scaled Jacobian, `None` when an edge +/// has zero length). +pub(crate) fn planar_corner(p: &[[f64; 2]; 4], i: usize) -> (f64, Option) { + let (a, b, c) = (p[(i + 3) % 4], p[i % 4], p[(i + 1) % 4]); + let (u, w) = ([c[0] - b[0], c[1] - b[1]], [a[0] - b[0], a[1] - b[1]]); + let cross = u[0] * w[1] - u[1] * w[0]; + let den = (u[0] * u[0] + u[1] * u[1]).sqrt() * (w[0] * w[0] + w[1] * w[1]).sqrt(); + let scaled = cross / den; + (cross, (den > 0.0 && scaled.is_finite()).then_some(scaled)) +} + +/// Robinson's aspect ratio and skew of a quadrilateral from its principal axes; `(+∞, 1)` when +/// an axis vanishes. +pub(crate) fn quad_axes(p: [P; 4]) -> (f64, f64) { + let x1 = [p[1][0] - p[0][0] + p[2][0] - p[3][0], p[1][1] - p[0][1] + p[2][1] - p[3][1], p[1][2] - p[0][2] + p[2][2] - p[3][2]]; + let x2 = [p[2][0] - p[1][0] + p[3][0] - p[0][0], p[2][1] - p[1][1] + p[3][1] - p[0][1], p[2][2] - p[1][2] + p[3][2] - p[0][2]]; + let (l1, l2) = (norm(x1), norm(x2)); + let aspect = ratio(l1.max(l2), l1.min(l2), f64::INFINITY).max(1.0); + let skew = ratio(dot(x1, x2).abs(), l1 * l2, 1.0).min(1.0); + (aspect, skew) +} + /// Wedge scaled Jacobian: the smallest determinant of the three unit edge vectors at a corner, /// taken in right-handed order (1 for a right prism on a right-angled corner, sin 60° ≈ 0.866 on /// an equilateral triangle). diff --git a/crates/mesh/src/size.rs b/crates/mesh/src/size.rs new file mode 100644 index 0000000..32fcf37 --- /dev/null +++ b/crates/mesh/src/size.rs @@ -0,0 +1,263 @@ +//! Size fields: the target element edge length at a point, and the edge divisions that follow it. +//! +//! One [`SizeField`] drives every mesher: the quad meshers ([`crate::quad`]) read it for interior +//! density, and the edge-division helpers here turn it into boundary division counts so that +//! faces sharing an edge share its nodes. The tetrahedral refinement reads the same trait. +//! +//! # Divisions +//! The number of divisions of a curve `C` under a size field `h` is the length of `C` measured in +//! units of `h`, `∫_C ds / h(s)` ([`ideal_divisions`]), rounded ([`edge_divisions`]). The nodes +//! of a divided curve split it into pieces of equal measured length ([`divide_polyline`]), so the +//! spacing grades with the field. +//! +//! # Parity +//! An all-quadrilateral mesh of a region needs an even number of boundary edges in total (every +//! quad has four sides, every interior edge is shared by two quads). The free quad mesher +//! ([`crate::quad::free_quad_mesh`]) needs an even count on *every* loop. [`loop_divisions`] +//! rounds the ideal counts of a loop's curves and fixes the parity of the loop with the change +//! that costs the least relative spacing error; [`even_divisions`] makes a single curve even, +//! which also keeps every loop corner at an even position (the free mesher keeps corners exactly +//! when they are there). +//! +//! The integrals use the midpoint rule on [`SAMPLES_PER_SEGMENT`] pieces per polyline segment, +//! with the size taken constant within a piece. +//! Every function is deterministic (no transcendental functions, only `+ − × ÷ √`). + +use crate::{Error, Result}; + +/// The target edge length at a point. Implementations return a finite, positive length; callers +/// treat anything else as invalid input. +pub trait SizeField { + fn size(&self, p: [f64; 3]) -> f64; +} + +impl f64> SizeField for F { + fn size(&self, p: [f64; 3]) -> f64 { + self(p) + } +} + +/// The same size everywhere. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct UniformSize(pub f64); + +impl SizeField for UniformSize { + fn size(&self, _: [f64; 3]) -> f64 { + self.0 + } +} + +/// The largest division count of one curve. +pub const MAX_DIVISIONS: usize = 1 << 20; +/// The most points a polyline handed to the helpers may have. +pub const MAX_POLYLINE_POINTS: usize = 1 << 20; +/// Midpoint-rule samples per polyline segment. +pub const SAMPLES_PER_SEGMENT: usize = 8; + +/// The size at `p`, or an error when the field returns something that is not a finite, positive +/// length. +pub fn checked_size(field: &dyn SizeField, p: [f64; 3]) -> Result { + let h = field.size(p); + if h.is_finite() && h > 0.0 { Ok(h) } else { Err(Error::Invalid(format!("size field returned {h} at {p:?}; sizes must be finite and positive"))) } +} + +fn check_polyline(polyline: &[[f64; 3]]) -> Result<()> { + if polyline.len() < 2 { + return Err(Error::Invalid(format!("a polyline needs at least 2 points, got {}", polyline.len()))); + } + if polyline.len() > MAX_POLYLINE_POINTS { + return Err(Error::Limit(format!("polyline has {} points; the limit is {MAX_POLYLINE_POINTS}", polyline.len()))); + } + modelcraft_model::check_coordinates(polyline.as_flattened())?; + Ok(()) +} + +fn lerp(a: [f64; 3], b: [f64; 3], t: f64) -> [f64; 3] { + [a[0] + (b[0] - a[0]) * t, a[1] + (b[1] - a[1]) * t, a[2] + (b[2] - a[2]) * t] +} + +fn dist(a: [f64; 3], b: [f64; 3]) -> f64 { + let d = [b[0] - a[0], b[1] - a[1], b[2] - a[2]]; + (d[0] * d[0] + d[1] * d[1] + d[2] * d[2]).sqrt() +} + +/// The polyline cut into [`SAMPLES_PER_SEGMENT`] pieces per segment, each with its measured +/// length `∫ ds / h` (midpoint rule): `(start, end, measure)`. +fn pieces(field: &dyn SizeField, polyline: &[[f64; 3]]) -> Result> { + check_polyline(polyline)?; + let k = SAMPLES_PER_SEGMENT as f64; + let mut out = Vec::with_capacity(SAMPLES_PER_SEGMENT * polyline.len()); + for w in polyline.windows(2) { + let [a, b] = w else { continue }; + let piece = dist(*a, *b) / k; + for s in 0..SAMPLES_PER_SEGMENT { + let (p, q) = (lerp(*a, *b, s as f64 / k), lerp(*a, *b, (s + 1) as f64 / k)); + let m = piece / checked_size(field, lerp(*a, *b, (s as f64 + 0.5) / k))?; + if !m.is_finite() { + return Err(Error::Invalid("size field is too small for this polyline".into())); + } + out.push((p, q, m)); + } + } + Ok(out) +} + +/// The ideal (real) number of divisions of a polyline: its length measured in units of the size +/// field. +pub fn ideal_divisions(field: &dyn SizeField, polyline: &[[f64; 3]]) -> Result { + let total: f64 = pieces(field, polyline)?.iter().map(|p| p.2).sum(); + if total.is_finite() { Ok(total) } else { Err(Error::Invalid("size field is too small for this polyline".into())) } +} + +fn rounded(ideal: f64) -> Result { + if !(ideal >= 0.0) || !ideal.is_finite() { + return Err(Error::Invalid(format!("ideal division count {ideal} is not a finite, non-negative number"))); + } + if ideal > MAX_DIVISIONS as f64 { + return Err(Error::Limit(format!("{ideal:.0} divisions exceed the limit of {MAX_DIVISIONS}"))); + } + Ok((ideal.round() as usize).max(1)) +} + +/// The division count of one polyline: its ideal count rounded, at least 1. +pub fn edge_divisions(field: &dyn SizeField, polyline: &[[f64; 3]]) -> Result { + rounded(ideal_divisions(field, polyline)?) +} + +/// The nearest even count to `ideal`, at least 2 (ties go up). +pub fn even_divisions(ideal: f64) -> Result { + let half = rounded(ideal / 2.0)?; + half.checked_mul(2).ok_or_else(|| Error::Limit("division count overflow".into())) +} + +/// The relative spacing error of `n` divisions against `ideal`, as `max(n, x) / min(n, x)`. +fn spacing_error(n: usize, ideal: f64) -> f64 { + let (n, x) = (n as f64, ideal.max(f64::MIN_POSITIVE)); + if n > x { n / x } else { x / n } +} + +/// Division counts for the curves of one loop with an even total: each ideal count is rounded +/// (at least 1); when the sum is odd, the one curve whose count changes by ±1 with the smallest +/// increase of its spacing error is changed (ties: the lower index, then +1 before −1). +pub fn loop_divisions(ideal: &[f64]) -> Result> { + if ideal.is_empty() { + return Err(Error::Invalid("a loop needs at least one curve".into())); + } + let mut counts = ideal.iter().map(|&x| rounded(x)).collect::>>()?; + let total = counts.iter().try_fold(0usize, |s, &n| s.checked_add(n)).ok_or_else(|| Error::Limit("division count overflow".into()))?; + if total % 2 == 0 { + return Ok(counts); + } + let mut best: Option<(f64, usize, bool)> = None; + for (i, (&n, &x)) in counts.iter().zip(ideal).enumerate() { + let now = spacing_error(n, x); + for up in [true, false] { + if !up && n <= 1 { + continue; + } + let m = if up { n + 1 } else { n - 1 }; + let cost = spacing_error(m, x) - now; + if best.is_none_or(|(c, _, _)| cost < c) { + best = Some((cost, i, up)); + } + } + } + if let Some((_, i, up)) = best + && let Some(n) = counts.get_mut(i) + { + *n = if up { *n + 1 } else { *n - 1 }; + } + Ok(counts) +} + +/// `n + 1` points along a polyline (the first and last are its end points, copied exactly) that +/// split it into `n` pieces of equal length measured by the size field. +pub fn divide_polyline(field: &dyn SizeField, polyline: &[[f64; 3]], n: usize) -> Result> { + if n == 0 || n > MAX_DIVISIONS { + return Err(Error::Invalid(format!("division count {n} is outside 1..={MAX_DIVISIONS}"))); + } + let pieces = pieces(field, polyline)?; + let total: f64 = pieces.iter().map(|p| p.2).sum(); + let (first, last) = match (polyline.first(), polyline.last()) { + (Some(a), Some(b)) => (*a, *b), + _ => return Err(Error::Invalid("empty polyline".into())), + }; + if !(total > 0.0) || !total.is_finite() { + return Err(Error::Invalid("polyline has zero length".into())); + } + let mut out = Vec::with_capacity(n + 1); + out.push(first); + let (mut piece, mut before) = (0usize, 0.0); + for k in 1..n { + let target = total * k as f64 / n as f64; + while let Some(&(_, _, m)) = pieces.get(piece) { + if before + m >= target || piece + 1 == pieces.len() { + break; + } + before += m; + piece += 1; + } + let (a, b, m) = pieces.get(piece).copied().ok_or_else(|| Error::Failed("polyline division walked off the end".into()))?; + // The size is taken constant within a piece. + let t = if m > 0.0 { ((target - before) / m).clamp(0.0, 1.0) } else { 0.0 }; + out.push(lerp(a, b, t)); + } + out.push(last); + Ok(out) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn uniform_divisions() { + let line = [[0.0, 0.0, 0.0], [10.0, 0.0, 0.0]]; + assert_eq!(edge_divisions(&UniformSize(1.0), &line).unwrap(), 10); + assert_eq!(edge_divisions(&UniformSize(3.0), &line).unwrap(), 3); + assert_eq!(edge_divisions(&UniformSize(100.0), &line).unwrap(), 1); + let pts = divide_polyline(&UniformSize(1.0), &line, 4).unwrap(); + assert_eq!(pts, vec![[0.0, 0.0, 0.0], [2.5, 0.0, 0.0], [5.0, 0.0, 0.0], [7.5, 0.0, 0.0], [10.0, 0.0, 0.0]]); + } + + #[test] + fn graded_division_is_finer_where_the_size_is_smaller() { + let field = |p: [f64; 3]| 0.5 + 0.1 * p[0]; + let line = [[0.0, 0.0, 0.0], [5.0, 0.0, 0.0], [10.0, 0.0, 0.0]]; + let pts = divide_polyline(&field, &line, 6).unwrap(); + let gaps: Vec = pts.windows(2).map(|w| w[1][0] - w[0][0]).collect(); + assert!(gaps.windows(2).all(|g| g[1] > g[0]), "{gaps:?}"); + assert_eq!(pts.last(), Some(&[10.0, 0.0, 0.0])); + } + + #[test] + fn parity() { + assert_eq!(even_divisions(3.2).unwrap(), 4); + assert_eq!(even_divisions(0.1).unwrap(), 2); + assert_eq!(even_divisions(5.0).unwrap(), 6); + let c = loop_divisions(&[3.0, 2.0, 3.0, 2.0]).unwrap(); + assert_eq!(c, vec![3, 2, 3, 2]); + let c = loop_divisions(&[3.0, 2.0, 3.0, 2.6]).unwrap(); + assert_eq!(c.iter().sum::() % 2, 0); + // The curve whose rounding was the worst absorbs the fix. + assert_eq!(c, vec![3, 2, 3, 2]); + let c = loop_divisions(&[1.0, 1.0, 1.0]).unwrap(); + assert_eq!(c, vec![2, 1, 1]); + } + + #[test] + fn hostile() { + let line = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]]; + assert!(edge_divisions(&UniformSize(0.0), &line).is_err()); + assert!(edge_divisions(&UniformSize(f64::NAN), &line).is_err()); + assert!(edge_divisions(&UniformSize(-1.0), &line).is_err()); + assert!(edge_divisions(&UniformSize(1e-300), &line).is_err()); + assert!(edge_divisions(&UniformSize(1.0), &[[0.0; 3]]).is_err()); + assert!(edge_divisions(&UniformSize(1.0), &[[0.0; 3], [f64::NAN, 0.0, 0.0]]).is_err()); + assert!(divide_polyline(&UniformSize(1.0), &[[0.0; 3], [0.0; 3]], 2).is_err()); + assert!(divide_polyline(&UniformSize(1.0), &line, 0).is_err()); + assert!(loop_divisions(&[]).is_err()); + assert!(loop_divisions(&[f64::INFINITY]).is_err()); + assert!(even_divisions(f64::NAN).is_err()); + } +} diff --git a/docs/architecture.md b/docs/architecture.md index 4af077f..f231cf7 100644 --- a/docs/architecture.md +++ b/docs/architecture.md @@ -68,7 +68,14 @@ Rules beyond the layer order: recombination (Pellerin et al. 2018), then polycube (DualCube labelling) for blocky parts; frame-field/integer-grid-map methods stay research. Validity uses a certified trilinear Jacobian check (Marschner et al. 2020), not corner Jacobians alone. -5. **Quality.** Scaled Jacobian, aspect ratio, dihedral angles; every generated mesh is checked +5. **Structured hex and face meshing.** Faces get all-quad meshes on shared, prescribed edge + nodes (`modelcraft-mesh::quad`): transfinite interpolation for four-sided faces (Gordon and + Hall 1973), and for free faces a half-resolution constrained Delaunay triangulation, + recombination (Lee and Lo 1994) and one uniform split that lands on the boundary nodes, then + smoothing and untangling; every quad is checked with exact predicates. One `SizeField` trait + (`modelcraft-mesh::size`) sets interior sizes and edge divisions, with even counts per loop. + Mapped and swept hex build on these face meshes. +6. **Quality.** Scaled Jacobian, aspect ratio, dihedral angles; every generated mesh is checked for positive Jacobians and conserved volume before it reaches the model. ## Commands diff --git a/docs/meshing/literature.md b/docs/meshing/literature.md index edd7781..f110ab7 100644 --- a/docs/meshing/literature.md +++ b/docs/meshing/literature.md @@ -783,6 +783,73 @@ Delaunay methods with Voronoi/ODT improvement, then (2) all-hex and hex-dominant - **Knupp 2001.** *Algebraic mesh quality metrics.* SIAM J. Sci. Comput. 23(1):193–218. Corner-Jacobian scaled Jacobian and shape metrics for hexes, wedges and quads. + +### B7. Planar quadrilateral meshing (implemented in `modelcraft-mesh::quad`) + +The face meshers under mapped and swept hex meshing and shell meshing. All quads on prescribed +boundary nodes, so faces sharing an edge share its nodes. DOIs below were checked through +`doi.org/api/handles` and Crossref on 2026-10-13; papers without a DOI are cited by venue only. + +**B7.1 Gordon and Hall 1973, transfinite interpolation.** **Verdict: BUILD ON (mapped meshing).** +- **Citation.** W. J. Gordon, C. A. Hall. *Construction of curvilinear co-ordinate systems and applications to mesh + generation.* IJNME 7(4):461–477, 1973. DOI [10.1002/nme.1620070405](https://doi.org/10.1002/nme.1620070405). +- **Use.** The bilinearly blended Coons map of a four-sided region, dimension-generic (`quad::transfinite_grid`), + with arc-length blending parameters for unequal spacing on opposite sides. + +**B7.2 Free all-quad: the methods compared.** +- **Paving.** T. D. Blacker, M. B. Stephenson. *Paving: a new approach to automated quadrilateral mesh generation.* + IJNME 32(4):811–847, 1991. DOI [10.1002/nme.1620320410](https://doi.org/10.1002/nme.1620320410). **Verdict: + reference.** Best boundary rows; front collisions are its failure mode. +- **Q-Morph.** S. J. Owen, M. L. Staten, S. A. Canann, S. Saigal. *Q-Morph: an indirect approach to advancing front + quad meshing.* IJNME 44(9):1317–1340, 1999. DOI + [10.1002/(SICI)1097-0207(19990330)44:9<1317::AID-NME532>3.0.CO;2-N](https://doi.org/10.1002/(SICI)1097-0207(19990330)44:9%3C1317::AID-NME532%3E3.0.CO;2-N). + **Verdict: reference (upgrade path for boundary alignment).** Front states, side-edge recovery by swaps, seam + closing. +- **Blossom-Quad.** J.-F. Remacle, J. Lambrechts, B. Seny, E. Marchandise, A. Johnen, C. Geuzaine. *Blossom-Quad: a + non-uniform quadrilateral mesh generator using a minimum-cost perfect-matching algorithm.* IJNME 89(9):1102–1119, + 2012. DOI [10.1002/nme.3279](https://doi.org/10.1002/nme.3279). **Verdict: reference (upgrade path for the + matching).** Paper only: the implementation lives in Gmsh, which is GPL and forbidden. +- **Recombination and split.** D. T. Lee, S. H. Lo. *A new scheme for the generation of a graded quadrilateral + mesh.* Computers & Structures 52(5):847–857, 1994. DOI + [10.1016/0045-7949(94)90070-1](https://doi.org/10.1016/0045-7949(94)90070-1). **Verdict: BUILD ON.** Triangles + merged into quads, then every element split once, which makes the mesh all-quad whatever the merge found. + The split topology is that of E. Catmull, J. Clark, CAD 10(6):350–355, 1978, DOI + [10.1016/0010-4485(78)90110-0](https://doi.org/10.1016/0010-4485(78)90110-0). +- **Why the split method.** All-quad by construction given an even point count per loop; robustness is that of the + constrained Delaunay triangulation (exists for every valid input; exact predicates); low implementation risk. + Its cost is irregular nodes and 120° corners in triangle-born quads, reduced by smoothing. The reasoning in full + is in the `quad` module docs. + +**B7.3 Planar constrained Delaunay triangulation and refinement.** **Verdict: BUILD ON.** +- A. Bowyer, *Computing Dirichlet tessellations*, Comput. J. 24(2):162–166, 1981, + DOI [10.1093/comjnl/24.2.162](https://doi.org/10.1093/comjnl/24.2.162); D. F. Watson, *Computing the n-dimensional + Delaunay tessellation…*, Comput. J. 24(2):167–172, 1981, DOI + [10.1093/comjnl/24.2.167](https://doi.org/10.1093/comjnl/24.2.167). +- S. W. Sloan. *A fast algorithm for generating constrained Delaunay triangulations.* Computers & Structures + 47(3):441–450, 1993. DOI [10.1016/0045-7949(93)90239-A](https://doi.org/10.1016/0045-7949(93)90239-A). Edge + recovery by flipping. +- J. Ruppert. *A Delaunay refinement algorithm for quality 2-dimensional mesh generation.* J. Algorithms + 18(3):548–585, 1995. DOI [10.1006/jagm.1995.1021](https://doi.org/10.1006/jagm.1995.1021); J. R. Shewchuk, + *Delaunay refinement algorithms for triangular mesh generation*, CGTA 22:21–74, 2002, DOI + [10.1016/S0925-7721(01)00047-5](https://doi.org/10.1016/S0925-7721(01)00047-5). Used without segment splitting + (boundary nodes are prescribed). + +**B7.4 Smoothing and untangling of the quads.** **Verdict: BUILD ON.** +- L. A. Freitag. *On combining Laplacian and optimization-based mesh smoothing techniques.* AMD-Vol. 220, ASME 1997; + S. A. Canann, J. R. Tristano, M. L. Staten, *An approach to combined Laplacian and optimization-based smoothing for + triangular, quadrilateral, and quad-dominant meshes*, IMR 1998. Laplacian moves kept only when the worst element + does not get worse. +- L. A. Freitag, P. Plassmann. *Local optimization-based simplicial mesh untangling and improvement.* IJNME + 49(1–2):109–125, 2000. DOI + [10.1002/1097-0207(20000910/20)49:1/2<109::AID-NME925>3.0.CO;2-U](https://doi.org/10.1002/1097-0207(20000910/20)49:1/2%3C109::AID-NME925%3E3.0.CO;2-U). + Untangling by maximizing the smallest (linear) corner area. + +**B7.5 Quadrilateral quality.** **Verdict: BUILD ON.** +- Knupp 2001 (B6.6) for the scaled Jacobian. J. Robinson, *CRE method of element testing and the Jacobian shape + parameters*, Eng. Comput. 4(2):113–118, 1987, for aspect ratio and skew from the principal axes. Planar validity + is exact: a planar bilinear quad's Jacobian is affine in each reference coordinate, so four exact `orient2d` corner + tests decide it. + --- ## Part C: Datasets and benchmarks