From 8dd583930e4c6c092387b4de82245bc71edd8316 Mon Sep 17 00:00:00 2001 From: Irakli Kourtis Date: Tue, 13 Oct 2026 14:45:54 -0400 Subject: [PATCH] Display tessellation plans (2/4): size rules and the facet angle (#803) Per-face target side lengths, the caller's facet angle, refined grid lines, and fixes for underflowing coordinates and straight edges ending exactly on their vertex. Nothing drawn changes yet. Refs #303 Co-Authored-By: Claude Opus 5.5 --- Cargo.lock | 1 + Cargo.toml | 1 + crates/kernel/Cargo.toml | 3 + crates/kernel/src/body.rs | 15 ++- crates/kernel/src/par.rs | 62 +++------ crates/kernel/src/tess/cdt.rs | 207 +++++++++++++++++++++++----- crates/kernel/src/tess/coons.rs | 57 ++++++-- crates/kernel/src/tess/edges.rs | 139 +++++++++++++++---- crates/kernel/src/tess/mod.rs | 173 ++++++++++++++++++------ crates/kernel/src/tess/size.rs | 231 ++++++++++++++++++++++++++++++++ crates/kernel/src/tess/tests.rs | 40 +++++- crates/kernel/src/tess/uv.rs | 130 ++++++++++++++---- 12 files changed, 870 insertions(+), 189 deletions(-) create mode 100644 crates/kernel/src/tess/size.rs diff --git a/Cargo.lock b/Cargo.lock index 1e74cafe..026c8bef 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -3365,6 +3365,7 @@ name = "solidcraft-kernel" version = "0.8.0" dependencies = [ "log", + "rayon", "serde", "serde_json", "solidcraft-geom", diff --git a/Cargo.toml b/Cargo.toml index bc4f8618..0269567f 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -40,6 +40,7 @@ truck-topology = "0.6" # constrained Delaunay + Ruppert refinement for trimmed faces (MIT OR Apache-2.0); already in the # tree through truck-meshalgo. spade = { version = "2.15", default-features = false, features = ["std"] } +rayon = "1.12" # packages (3MF) zip = { version = "3", default-features = false, features = ["deflate"] } diff --git a/crates/kernel/Cargo.toml b/crates/kernel/Cargo.toml index 18712825..a728a365 100644 --- a/crates/kernel/Cargo.toml +++ b/crates/kernel/Cargo.toml @@ -20,3 +20,6 @@ truck-meshalgo = { workspace = true } truck-stepio = { workspace = true } truck-topology = { workspace = true } spade = { workspace = true } + +[target.'cfg(not(target_arch = "wasm32"))'.dependencies] +rayon = { workspace = true } diff --git a/crates/kernel/src/body.rs b/crates/kernel/src/body.rs index 31d8a0bd..0929eed8 100644 --- a/crates/kernel/src/body.rs +++ b/crates/kernel/src/body.rs @@ -1051,12 +1051,25 @@ impl Body { /// triangulation in its own parameters, with each edge sampled once (chord tolerance /// `tol`, mm) so faces meet exactly. pub fn display_plan(&self, tol: f64) -> Result { + self.plan_within(tol, None) + } + + /// Like [`Body::display_plan`], and no facet turns more than `angle` radians (see + /// [`Body::display_mesh_angular`]). + pub fn display_plan_angular(&self, tol: f64, angle: f64) -> Result { + self.plan_within(tol, Some(angle)) + } + + fn plan_within(&self, tol: f64, angle: Option) -> Result { if self.is_mesh() { return Err(KernelError::Invalid("a mesh body has no surfaces to tessellate".into())); } let tol = if tol.is_finite() && tol > 0.0 { tol } else { 0.05 }; let tol = tol.max(self.size() * 1e-6); - guard("tessellate", || crate::tess::plan(&self.solid, tol).map(|(plan, _)| plan)) + // The same angle and floor as the reference mesher's angular display mesh. + let turn = + angle.filter(|a| a.is_finite()).map(|a| crate::tess::size::Turn { angle: a.max(MIN_DISPLAY_ANGLE), floor: tol * DISPLAY_ANGLE_FLOOR }); + guard("tessellate", || crate::tess::plan(&self.solid, tol, turn).map(|(plan, _)| plan)) } /// Faces with no triangles in a coarse tessellation (the mesher failed or gave up). diff --git a/crates/kernel/src/par.rs b/crates/kernel/src/par.rs index 930a9868..b0f281c8 100644 --- a/crates/kernel/src/par.rs +++ b/crates/kernel/src/par.rs @@ -1,53 +1,22 @@ -//! Ordered, bounded parallel map over independent items (desktop only; the browser build runs -//! serially). Each item runs under [`crate::guard`]; a worker that cannot be spawned runs its -//! chunk on the calling thread, and a worker that dies reports an internal error for its items. +//! Ordered parallel map over independent items, on the process-wide worker pool (desktop +//! only; the browser build runs serially). Each item runs under [`crate::guard`], so a panic in +//! one item is that item's error and never reaches the pool. -use crate::{KernelError, Result, guard}; +use crate::{Result, guard}; -/// Most worker threads. -const MAX_WORKERS: usize = 8; - -pub(crate) fn workers() -> usize { - #[cfg(not(target_arch = "wasm32"))] - { - std::thread::available_parallelism().map_or(1, |n| n.get().min(MAX_WORKERS)) - } - #[cfg(target_arch = "wasm32")] - { - 1 - } -} - -/// `f` over `items`, results in item order. Runs serially below `min_items` items. -pub(crate) fn map(what: &'static str, items: &[T], min_items: usize, f: impl Fn(&T) -> Result + Sync) -> Vec> { +/// `f` over `items`, results in item order: on the pool when `parallel`, else on the calling +/// thread. Waking the pool's workers costs tens of microseconds, so callers ask for it only +/// when the items carry more work than that (decided from the work, not the item count). +pub(crate) fn map(what: &'static str, items: &[T], parallel: bool, f: impl Fn(&T) -> Result + Sync) -> Vec> { let one = |item: &T| guard(what, || f(item)); + // The browser build has no worker pool: everything runs on the calling thread. + #[cfg(target_arch = "wasm32")] + let _ = parallel; #[cfg(not(target_arch = "wasm32"))] { - let workers = workers(); - if items.len() >= min_items.max(2) && workers > 1 { - return std::thread::scope(|scope| { - let size = items.len().div_ceil(workers); - let one = &one; - let pending: Vec<_> = items - .chunks(size) - .map(|chunk| { - let handle = - std::thread::Builder::new().name(what.into()).spawn_scoped(scope, move || chunk.iter().map(one).collect::>()).ok(); - (chunk, handle) - }) - .collect(); - let mut out = Vec::with_capacity(items.len()); - for (chunk, handle) in pending { - match handle { - Some(h) => match h.join() { - Ok(results) => out.extend(results), - Err(_) => out.extend(chunk.iter().map(|_| Err(KernelError::Internal(format!("{what}: worker stopped"))))), - }, - None => out.extend(chunk.iter().map(one)), - } - } - out - }); + use rayon::prelude::*; + if parallel && items.len() > 1 && rayon::current_num_threads() > 1 { + return items.par_iter().map(one).collect(); } } items.iter().map(one).collect() @@ -56,11 +25,12 @@ pub(crate) fn map(what: &'static str, items: &[T], min_items: #[cfg(test)] mod tests { use super::*; + use crate::KernelError; #[test] fn keeps_order_and_errors() { let items: Vec = (0..1000).collect(); - let out = map("test", &items, 1, |x| if *x == 500 { Err(KernelError::Invalid("five hundred".into())) } else { Ok(x * 2) }); + let out = map("test", &items, true, |x| if *x == 500 { Err(KernelError::Invalid("five hundred".into())) } else { Ok(x * 2) }); assert_eq!(out.len(), 1000); for (i, r) in out.iter().enumerate() { match r { diff --git a/crates/kernel/src/tess/cdt.rs b/crates/kernel/src/tess/cdt.rs index 1c317ced..1dfed9c0 100644 --- a/crates/kernel/src/tess/cdt.rs +++ b/crates/kernel/src/tess/cdt.rs @@ -1,28 +1,41 @@ -//! Faces the structured grid cannot take (trimmed faces, holes): a constrained Delaunay -//! triangulation of the boundary loops in (u, v) scaled by the surface's mean metric, with -//! interior seed points, refinement wherever a triangle sags more than the tolerance from the -//! surface, then Ruppert refinement for shape. Constraint edges are never split, so the +//! Faces the structured grid cannot take (planar faces, trimmed faces, holes): a constrained +//! Delaunay triangulation of the boundary loops in (u, v) scaled by the surface's mean metric, +//! seeded with an equilateral lattice, refined wherever a triangle sags more than the +//! tolerance from the surface, has a side longer than the target or bends more than +//! the facet angle, then Ruppert refinement for shape. Constraint edges are never split, so the //! boundary stays exactly the shared edge nodes. use std::collections::HashMap; use solidcraft_geom::Vec3; +use super::size::Turn; use spade::handles::FixedVertexHandle; use spade::{AngleLimit, ConstrainedDelaunayTriangulation, Point2, RefinementParameters, Triangulation}; type Cdt = ConstrainedDelaunayTriangulation>; /// Smallest angle Ruppert refinement aims for. -const MIN_ANGLE_DEG: f64 = 20.0; +const MIN_ANGLE_DEG: f64 = 25.0; +/// A side is split once it is this many times the target length (the lattice's own sides, +/// exactly the target, and the slightly longer ones where it meets the boundary, stay whole). +const LONG: f64 = 1.5; +/// Seeds closer than this fraction of the side length to the boundary are left out (the +/// boundary nodes are already there; a seed beside them makes a sliver). +const CLEARANCE: f64 = 0.5; /// Rounds of sag refinement. const SAG_ROUNDS: usize = 16; /// The surface, for refining where the triangles leave it. pub(crate) struct Sag<'a> { - pub eval: &'a dyn Fn(f64, f64) -> Option, + /// Point and unit normal (when not degenerate) at (u, v). + pub eval: &'a dyn Fn(f64, f64) -> Option<(Vec3, Option)>, pub tol: f64, + /// Longest triangle side (mm). + pub h: f64, + /// The facet angle, if any: corner normals turning more get a point in the middle. + pub turn: Option, /// Scale of u and v in the triangulation (the surface's mean |S_u|, |S_v|). pub scale: [f64; 2], } @@ -33,6 +46,10 @@ pub(crate) struct Fill<'a> { pub seeds: &'a [[f64; 2]], /// Refinement, adding at most this many points in all. pub refine: Option<(Sag<'a>, usize)>, + /// Ruppert refinement for shape after it. Shape matters for shading only where the normal + /// varies: a plane is shaded the same by any triangles, so planar faces skip it (around many + /// small holes it grades out from each rim and multiplies the triangles). + pub shape: bool, } /// Triangulate the region bounded by `loops` (ranges of `boundary`, each closed). Boundary @@ -46,7 +63,9 @@ pub(crate) fn triangulate(boundary: &[[f64; 2]], loops: &[std::ops::Range Some((sag, _)) if sag.scale.iter().all(|k| k.is_finite() && *k > 0.0) => sag.scale, _ => [1.0, 1.0], }; - let to = |p: &[f64; 2]| Point2::new(p[0] * scale[0], p[1] * scale[1]); + // spade refuses coordinates below about 1e-43 in magnitude (other than zero), which + // projection leaves where a parameter is zero in exact arithmetic. + let to = |p: &[f64; 2]| spade::mitigate_underflow(Point2::new(p[0] * scale[0], p[1] * scale[1])); let from = |p: Point2| [p.x / scale[0], p.y / scale[1]]; let mut cdt = Cdt::new(); let mut index: HashMap = HashMap::new(); @@ -77,12 +96,14 @@ pub(crate) fn triangulate(boundary: &[[f64; 2]], loops: &[std::ops::Range if let Some((sag, budget)) = &fill.refine { let constraints = cdt.num_constraints(); let added = refine_sag(&mut cdt, sag, *budget, &from); - let params = RefinementParameters::::new() - .keep_constraint_edges() - .exclude_outer_faces(true) - .with_angle_limit(AngleLimit::from_deg(MIN_ANGLE_DEG)) - .with_max_additional_vertices(budget.saturating_sub(added)); - cdt.refine(params); + if fill.shape { + let params = RefinementParameters::::new() + .keep_constraint_edges() + .exclude_outer_faces(true) + .with_angle_limit(AngleLimit::from_deg(MIN_ANGLE_DEG)) + .with_max_additional_vertices(budget.saturating_sub(added)); + cdt.refine(params); + } // A split constraint would put a point on the boundary the neighbour does not have. if cdt.num_constraints() != constraints { return None; @@ -114,12 +135,18 @@ pub(crate) fn triangulate(boundary: &[[f64; 2]], loops: &[std::ops::Range Some((extra, tris)) } -/// Insert the parameters of the surface point under every triangle centre and every -/// non-constraint edge middle that lies more than the tolerance from the triangle, until none -/// does or `budget` points were added. Returns the points added. +/// Insert the parameters of the surface point under every triangle centre that lies more than +/// the tolerance from the triangle or whose corner normals turn more than the facet angle, and +/// under the middle of every non-constraint side that sags more than the tolerance or is +/// longer than [`LONG`] times `h`, until none does or `budget` points were added. Returns the points added. fn refine_sag(cdt: &mut Cdt, sag: &Sag, budget: usize, from: &dyn Fn(Point2) -> [f64; 2]) -> usize { let mut added = 0; - let mut cache: HashMap> = HashMap::new(); + let mut cache: HashMap)>> = HashMap::new(); + // The widest turn between two corner normals. + let spread = |n: [Option; 3]| match n { + [Some(a), Some(b), Some(c)] => [(a, b), (b, c), (c, a)].iter().map(|(x, y)| x.dot(*y).clamp(-1.0, 1.0).acos()).fold(0.0, f64::max), + _ => 0.0, + }; for _ in 0..SAG_ROUNDS { let mut todo: Vec> = Vec::new(); let mut at = |h: FixedVertexHandle, p: Point2| { @@ -136,12 +163,17 @@ fn refine_sag(cdt: &mut Cdt, sag: &Sag, budget: usize, from: &dyn Fn(Point2 let vs = face.vertices(); let ps = vs.map(|v| v.position()); let c = Point2::new((ps[0].x + ps[1].x + ps[2].x) / 3.0, (ps[0].y + ps[1].y + ps[2].y) / 3.0); - let (Some(a), Some(b), Some(d)) = (at(vs[0].fix(), ps[0]), at(vs[1].fix(), ps[1]), at(vs[2].fix(), ps[2])) else { continue }; - let off = |q: Point2, chord: Vec3| { - let uv = from(q); - (sag.eval)(uv[0], uv[1]).is_some_and(|s| s.dist(chord) > sag.tol) + let (Some((a, na)), Some((b, nb)), Some((d, nd))) = (at(vs[0].fix(), ps[0]), at(vs[1].fix(), ps[1]), at(vs[2].fix(), ps[2])) else { + continue; }; - if off(c, (a + b + d) * (1.0 / 3.0)) { + // How far the surface lies from the triangle at `q`. + let gap = |q: Point2, chord: Vec3| { + let uv = from(q); + (sag.eval)(uv[0], uv[1]).map_or(0.0, |(s, _)| s.dist(chord)) + }; + let off = |q: Point2, chord: Vec3| gap(q, chord) > sag.tol; + let middle = gap(c, (a + b + d) * (1.0 / 3.0)); + if middle > sag.tol || sag.turn.is_some_and(|t| t.exceeds(spread([na, nb, nd]), middle)) { todo.push(c); continue; } @@ -151,7 +183,7 @@ fn refine_sag(cdt: &mut Cdt, sag: &Sag, budget: usize, from: &dyn Fn(Point2 } let [p, q] = e.positions(); let m = Point2::new(0.5 * (p.x + q.x), 0.5 * (p.y + q.y)); - if off(m, (pa + pb) * 0.5) { + if pa.dist(pb) > LONG * sag.h || off(m, (pa + pb) * 0.5) { todo.push(m); } } @@ -165,7 +197,7 @@ fn refine_sag(cdt: &mut Cdt, sag: &Sag, budget: usize, from: &dyn Fn(Point2 if added >= budget { return added; } - if cdt.insert(p).is_ok() { + if cdt.insert(spade::mitigate_underflow(p)).is_ok() { added += 1; } } @@ -202,12 +234,24 @@ fn region(cdt: &Cdt) -> Vec { parity.into_iter().map(|p| p.unwrap_or(false)).collect() } -/// The grid points `us` × `vs` inside the region bounded by the loops (even–odd, so holes are -/// out): per grid line, the loops' crossings sorted once, then each point by binary search. -pub(crate) fn inside_grid(us: &[f64], vs: &[f64], polys: &[&[[f64; 2]]]) -> Vec<[f64; 2]> { +/// Seeds for a face: an equilateral (hexagonal) lattice over the parameter box `bbox` with +/// side `step` in each parameter (rows `step[1]·√3/2` apart, every other row shifted half a +/// step), keeping the points inside the region bounded by the loops (even–odd, so holes are +/// out) and at least [`CLEARANCE`] of a side from the boundary in the metric `scale`. +pub(crate) fn lattice(bbox: [(f64, f64); 2], step: [f64; 2], polys: &[&[[f64; 2]]], scale: [f64; 2]) -> Vec<[f64; 2]> { let mut out = Vec::new(); + if !step.iter().all(|s| s.is_finite() && *s > 0.0) || !bbox.iter().all(|(a, b)| a.is_finite() && b.is_finite() && b > a) { + return out; + } + let row = step[1] * 3f64.sqrt() * 0.5; + let (rows, cols) = (((bbox[1].1 - bbox[1].0) / row).ceil(), ((bbox[0].1 - bbox[0].0) / step[0]).ceil()); + if !(rows * cols <= MAX_SEEDS as f64) { + return out; + } + let near = Clearance::new(polys, scale, CLEARANCE * step[0] * scale[0]); let mut xs: Vec = Vec::new(); - for &v in vs { + for j in 1..rows as usize { + let v = bbox[1].0 + row * j as f64; xs.clear(); for poly in polys { let n = poly.len(); @@ -219,8 +263,10 @@ pub(crate) fn inside_grid(us: &[f64], vs: &[f64], polys: &[&[[f64; 2]]]) -> Vec< } } xs.sort_by(f64::total_cmp); - for &u in us { - if xs.partition_point(|x| *x < u) % 2 == 1 { + let shift = if j % 2 == 1 { 0.5 * step[0] } else { 0.0 }; + for i in 0..=cols as usize { + let u = bbox[0].0 + shift + step[0] * i as f64; + if u < bbox[0].1 && xs.partition_point(|x| *x < u) % 2 == 1 && !near.within([u, v]) { out.push([u, v]); } } @@ -228,6 +274,71 @@ pub(crate) fn inside_grid(us: &[f64], vs: &[f64], polys: &[&[[f64; 2]]]) -> Vec< out } +/// Most lattice seeds of one face. +const MAX_SEEDS: usize = 200_000; + +/// Boundary segments bucketed on a square grid (cell = the clearance), for "is this point +/// within the clearance of the boundary" in constant time per point. +struct Clearance { + cells: HashMap<(i64, i64), Vec<([f64; 2], [f64; 2])>>, + scale: [f64; 2], + r: f64, +} + +impl Clearance { + fn new(polys: &[&[[f64; 2]]], scale: [f64; 2], r: f64) -> Clearance { + let mut cells: HashMap<(i64, i64), Vec<([f64; 2], [f64; 2])>> = HashMap::new(); + if !(r.is_finite() && r > 0.0) { + return Clearance { cells, scale, r: 0.0 }; + } + for poly in polys { + let n = poly.len(); + for k in 0..n { + let (Some(a), Some(b)) = (poly.get(k), poly.get((k + 1) % n)) else { continue }; + let (a, b) = ([a[0] * scale[0], a[1] * scale[1]], [b[0] * scale[0], b[1] * scale[1]]); + let cell = |x: f64| (x / r).floor() as i64; + let (x0, x1) = (cell(a[0].min(b[0])), cell(a[0].max(b[0]))); + let (y0, y1) = (cell(a[1].min(b[1])), cell(a[1].max(b[1]))); + // A segment much longer than the clearance (a straight side of a large face) + // is checked from every cell of its box; cap that box. + if (x1 - x0 + 1).saturating_mul(y1 - y0 + 1) > 4096 { + continue; + } + for x in x0..=x1 { + for y in y0..=y1 { + cells.entry((x, y)).or_default().push((a, b)); + } + } + } + } + Clearance { cells, scale, r } + } + + fn within(&self, p: [f64; 2]) -> bool { + if self.r == 0.0 { + return false; + } + let p = [p[0] * self.scale[0], p[1] * self.scale[1]]; + let (cx, cy) = ((p[0] / self.r).floor() as i64, (p[1] / self.r).floor() as i64); + for x in cx - 1..=cx + 1 { + for y in cy - 1..=cy + 1 { + let Some(segs) = self.cells.get(&(x, y)) else { continue }; + if segs.iter().any(|(a, b)| dist_to_segment(p, *a, *b) < self.r) { + return true; + } + } + } + false + } +} + +fn dist_to_segment(p: [f64; 2], a: [f64; 2], b: [f64; 2]) -> f64 { + let (dx, dy) = (b[0] - a[0], b[1] - a[1]); + let l2 = dx * dx + dy * dy; + let t = if l2 > 0.0 { (((p[0] - a[0]) * dx + (p[1] - a[1]) * dy) / l2).clamp(0.0, 1.0) } else { 0.0 }; + (p[0] - a[0] - t * dx).hypot(p[1] - a[1] - t * dy) +} + #[cfg(test)] mod tests { use super::*; @@ -246,11 +357,13 @@ mod tests { let mut b = vec![[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]; b.extend([[0.4, 0.4], [0.4, 0.6], [0.6, 0.6], [0.6, 0.4]]); let loops = [0..4, 4..8]; - let flat = |u: f64, v: f64| Some(Vec3::new(u, v, 0.0)); + let flat = |u: f64, v: f64| Some((Vec3::new(u, v, 0.0), Some(Vec3::Z))); for refine in [false, true] { - let refine = refine.then_some((Sag { eval: &flat, tol: 1e-3, scale: [2.0, 1.0] }, 200)); + let refine = refine.then_some((Sag { eval: &flat, tol: 1e-3, h: 0.1, turn: None, scale: [2.0, 1.0] }, 2000)); let label = refine.is_some(); - let Some((extra, tris)) = triangulate(&b, &loops, Fill { seeds: &[[0.2, 0.2], [0.5, 0.5]], refine }) else { panic!("meshes") }; + let Some((extra, tris)) = triangulate(&b, &loops, Fill { seeds: &[[0.2, 0.2], [0.5, 0.5]], refine, shape: true }) else { + panic!("meshes") + }; let mut uv = b.clone(); uv.extend(extra); assert!((area(&uv, &tris) - 0.96).abs() < 1e-12, "refined: {label}"); @@ -271,6 +384,32 @@ mod tests { fn crossing_loops_are_refused() { let b = [[0.0, 0.0], [1.0, 1.0], [1.0, 0.0], [0.0, 1.0]]; let whole = std::slice::from_ref(&(0..4)); - assert!(triangulate(&b, whole, Fill { seeds: &[], refine: None }).is_none()); + assert!(triangulate(&b, whole, Fill { seeds: &[], refine: None, shape: false }).is_none()); + } + + #[test] + fn lattice_seeds_are_equilateral_and_clear_of_the_boundary() { + let outer = [[0.0, 0.0], [10.0, 0.0], [10.0, 10.0], [0.0, 10.0]]; + let hole = [[4.0, 4.0], [4.0, 6.0], [6.0, 6.0], [6.0, 4.0]]; + let seeds = lattice([(0.0, 10.0), (0.0, 10.0)], [1.0, 1.0], &[&outer, &hole], [1.0, 1.0]); + assert!(!seeds.is_empty()); + for p in &seeds { + let clear = |poly: &[[f64; 2]]| (0..4).all(|k| dist_to_segment(*p, poly[k], poly[(k + 1) % 4]) >= 0.5); + assert!(clear(&outer) && clear(&hole), "{p:?}"); + } + // Nearest neighbours one side apart, six of them in the open. + let near = |p: [f64; 2]| seeds.iter().filter(|q| ((q[0] - p[0]).hypot(q[1] - p[1]) - 1.0).abs() < 1e-9).count(); + assert!(seeds.iter().any(|p| near(*p) == 6)); + assert!(seeds.iter().all(|p| seeds.iter().all(|q| p == q || (q[0] - p[0]).hypot(q[1] - p[1]) > 1.0 - 1e-9))); + } + + #[test] + fn coordinates_that_underflow_are_taken_as_zero() { + // A parameter that is zero in exact arithmetic, left at 1e-300 by projection. + let b = [[1e-300, 0.0], [1.0, -1e-310], [1.0, 1.0], [0.0, 1.0]]; + let whole = std::slice::from_ref(&(0..4)); + let Some((extra, tris)) = triangulate(&b, whole, Fill { seeds: &[], refine: None, shape: false }) else { panic!("meshes") }; + assert!(extra.is_empty()); + assert_eq!(tris.len(), 2); } } diff --git a/crates/kernel/src/tess/coons.rs b/crates/kernel/src/tess/coons.rs index a7056600..6f79d518 100644 --- a/crates/kernel/src/tess/coons.rs +++ b/crates/kernel/src/tess/coons.rs @@ -88,7 +88,15 @@ fn lines(a: &Side, b: &Side) -> Vec { /// Triangulate the loop `uv` (with its fixed 3D nodes) as a Coons patch: boundary vertices /// keep their indices, the interior grid follows. Triangles counter-clockwise in (u, v). -pub(crate) fn coons(uv: &[[f64; 2]], fixed: &[Vec3], scale: [f64; 2], max_points: usize) -> Option<(Vec<[f64; 2]>, Vec<[u32; 3]>)> { +/// `fine(a, b)`: whether the surface between two parameter points needs no point between them +/// (the grid lines, taken from the sides, are split inside the face where it does not). +pub(crate) fn coons( + uv: &[[f64; 2]], + fixed: &[Vec3], + scale: [f64; 2], + max_points: usize, + fine: &dyn Fn([f64; 2], [f64; 2]) -> bool, +) -> Option<(Vec<[f64; 2]>, Vec<[u32; 3]>)> { let n = uv.len(); let [c0, c1, c2, c3] = corners(uv, scale)?; // Sides in loop order; the bottom (0) and right (1) run with the loop, the top (2) and @@ -97,22 +105,27 @@ pub(crate) fn coons(uv: &[[f64; 2]], fixed: &[Vec3], scale: [f64; 2], max_points let right = side(fixed, n, c1, c2)?; let top = reversed(side(fixed, n, c2, c3)?); let left = reversed(side(fixed, n, c3, c0)?); - let (ss, ts) = (lines(&bottom, &top), lines(&left, &right)); + let corner = |i: usize| uv.get(i).copied(); + let (p00, p10, p11, p01) = (corner(c0)?, corner(c1)?, corner(c2)?, corner(c3)?); + // The blend at (s, t): the sides' points, less the corners counted twice. + let at = |s: f64, t: f64| -> Option<[f64; 2]> { + let (b, tp, l, r) = (along(&bottom, uv, s)?, along(&top, uv, s)?, along(&left, uv, t)?, along(&right, uv, t)?); + let blend = |d: usize| { + (1.0 - t) * b[d] + t * tp[d] + (1.0 - s) * l[d] + s * r[d] + - ((1.0 - s) * (1.0 - t) * p00[d] + s * (1.0 - t) * p10[d] + (1.0 - s) * t * p01[d] + s * t * p11[d]) + }; + Some([blend(0), blend(1)]) + }; + let ss = refine(lines(&bottom, &top), max_points, &|a, b, x| Some(fine(at(a, x)?, at(b, x)?))); + let ts = refine(lines(&left, &right), max_points, &|a, b, x| Some(fine(at(x, a)?, at(x, b)?))); if ss.len().checked_mul(ts.len())? > max_points { return None; } - let corner = |i: usize| uv.get(i).copied(); - let (p00, p10, p11, p01) = (corner(c0)?, corner(c1)?, corner(c2)?, corner(c3)?); let mut out: Vec<[f64; 2]> = uv.to_vec(); let base = u32::try_from(n).ok()?; for &t in &ts { for &s in &ss { - let (b, tp, l, r) = (along(&bottom, uv, s)?, along(&top, uv, s)?, along(&left, uv, t)?, along(&right, uv, t)?); - let blend = |d: usize| { - (1.0 - t) * b[d] + t * tp[d] + (1.0 - s) * l[d] + s * r[d] - - ((1.0 - s) * (1.0 - t) * p00[d] + s * (1.0 - t) * p10[d] + (1.0 - s) * t * p01[d] + s * t * p11[d]) - }; - out.push([blend(0), blend(1)]); + out.push(at(s, t)?); } } let (ns, nt) = (ss.len(), ts.len()); @@ -154,6 +167,26 @@ pub(crate) fn coons(uv: &[[f64; 2]], fixed: &[Vec3], scale: [f64; 2], max_points Some((extra, tris)) } +/// Interior lines (strictly between 0 and 1) with every interval, ends included, split while +/// `fine(a, b, x)` fails at some position `x` across (or the lines reach `max`). +fn refine(lines: Vec, max: usize, fine: &dyn Fn(f64, f64, f64) -> Option) -> Vec { + const ACROSS: [f64; 5] = [0.0, 0.25, 0.5, 0.75, 1.0]; + let all: Vec = std::iter::once(0.0).chain(lines).chain(std::iter::once(1.0)).collect(); + let mut out = vec![0.0]; + let mut stack: Vec<(f64, f64, u32)> = all.windows(2).rev().map(|w| (w[0], w[1], 0)).collect(); + while let Some((a, b, depth)) = stack.pop() { + let ok = depth >= 16 || out.len() + stack.len() >= max || ACROSS.iter().all(|x| fine(a, b, *x).unwrap_or(true)); + if ok { + out.push(b); + } else { + let m = 0.5 * (a + b); + stack.push((m, b, depth + 1)); + stack.push((a, m, depth + 1)); + } + } + out.into_iter().filter(|x| *x > 0.0 && *x < 1.0).collect() +} + fn reversed(s: Side) -> Side { Side { nodes: s.nodes.into_iter().rev().collect(), at: s.at.into_iter().rev().map(|x| 1.0 - x).collect() } } @@ -218,7 +251,7 @@ mod tests { b.reverse(); } let fixed: Vec = b.iter().map(|p| Vec3::new(p[0], p[1], 0.0)).collect(); - let Some((extra, tris)) = coons(&b, &fixed, [1.0, 1.0], 10_000) else { panic!("meshes (k {k})") }; + let Some((extra, tris)) = coons(&b, &fixed, [1.0, 1.0], 10_000, &|_, _| true) else { panic!("meshes (k {k})") }; let mut uv = b.clone(); uv.extend(extra); assert!((area(&uv, &tris) - 16.0).abs() < 1e-9, "k {k} reverse {reverse}: {}", area(&uv, &tris)); @@ -230,6 +263,6 @@ mod tests { fn five_corners_is_not_a_coons_patch() { let b = [[0.0, 0.0], [1.0, 0.0], [1.0, 0.5], [0.5, 1.0], [0.0, 1.0]]; let fixed: Vec = b.iter().map(|p| Vec3::new(p[0], p[1], 0.0)).collect(); - assert!(coons(&b, &fixed, [1.0, 1.0], 100).is_none()); + assert!(coons(&b, &fixed, [1.0, 1.0], 100, &|_, _| true).is_none()); } } diff --git a/crates/kernel/src/tess/edges.rs b/crates/kernel/src/tess/edges.rs index 6dc20ce5..1955c7b6 100644 --- a/crates/kernel/src/tess/edges.rs +++ b/crates/kernel/src/tess/edges.rs @@ -5,6 +5,7 @@ use solidcraft_geom::Vec3; use truck_meshalgo::prelude::*; use truck_modeling as mt; +use super::size::Rules; use crate::body::from_p3; /// Most nodes on one edge. @@ -12,15 +13,26 @@ pub(crate) const MAX_EDGE_NODES: usize = 20_000; /// Deepest bisection of one edge segment. const MAX_DEPTH: u32 = 24; -/// Nodes of `edge` from its absolute front to its absolute back, every chord within `tol` of -/// the curve (checked at the quarter points, so an S-bend in one segment is not missed). The -/// end nodes are the vertices' own points, so edges meeting at a vertex meet exactly. -pub(crate) fn sample_edge(edge: &mt::Edge, tol: f64) -> Vec { +/// Nodes of `edge` from its absolute front to its absolute back: every chord within the +/// tolerance of the curve (checked at the quarter points, so an S-bend in one segment is not +/// missed), at most `h` long, and within the rules' facet angle. The end nodes are the +/// vertices' own points, so edges meeting at a vertex meet exactly. +pub(crate) fn sample_edge(edge: &mt::Edge, rules: &Rules, h: f64) -> Vec { + let tol = rules.tol; let front = from_p3(edge.absolute_front().point()); let back = from_p3(edge.absolute_back().point()); let curve = edge.curve(); if matches!(curve, mt::Curve::Line(_)) { - return vec![front, back]; + // Evenly, as few segments as keep each at most `h`. The ends are the vertices' own + // points: `front + (back − front) · 1` is not `back` to the last bit, and the face that + // runs the edge backwards starts its loop with that node, where the face running it + // forwards takes the next edge's first node (the vertex) — a crack. + let n = if h > 0.0 { (front.dist(back) / h).ceil().clamp(1.0, (MAX_EDGE_NODES - 1) as f64) as usize } else { 1 }; + let mut pts: Vec = (0..=n).map(|i| front + (back - front) * (i as f64 / n as f64)).collect(); + if let Some(p) = pts.last_mut() { + *p = back; + } + return pts; } let (t0, t1) = curve.range_tuple(); if !(t0.is_finite() && t1.is_finite() && t1 > t0) { @@ -36,36 +48,77 @@ pub(crate) fn sample_edge(edge: &mt::Edge, tol: f64) -> Vec { } .filter(|p| p.is_finite()) }; - // A closed curve's chord between its ends is empty: start from a few segments. - let initial: u32 = 4; + // A closed curve's chord between its ends is empty: start it from a few segments. An open + // one starts whole, so an arc the tolerance leaves flat stays one segment. + let initial: u32 = if edge.absolute_front().id() == edge.absolute_back().id() { 4 } else { 1 }; let t_at = |i: u32| t0 + (t1 - t0) * f64::from(i) / f64::from(initial); - let mut out_t = vec![t0]; - // Depth-first bisection, segment by segment, keeping the output ordered. - let mut stack: Vec<(f64, f64, u32)> = (0..initial).rev().map(|i| (t_at(i), t_at(i + 1), 0)).collect(); - while let Some((a, b, depth)) = stack.pop() { - let flat = depth >= MAX_DEPTH || out_t.len() + stack.len() >= MAX_EDGE_NODES || { - let (Some(pa), Some(pb)) = (at(a), at(b)) else { return vec![front, back] }; - [0.25, 0.5, 0.75].iter().all(|f| at(a + (b - a) * f).is_none_or(|q| q.dist_to_segment(pa, pb) <= tol)) + let mut ends: Vec = Vec::with_capacity(initial as usize + 1); + for i in 0..=initial { + let Some(p) = at(t_at(i)) else { return vec![front, back] }; + ends.push(p); + } + // Depth-first bisection, segment by segment, keeping the output ordered. Each segment + // carries its end points and middle, so every parameter is evaluated once (an + // intersection curve's point is a Newton solve). + struct Seg { + a: f64, + b: f64, + pa: Vec3, + pb: Vec3, + mid: Option, + depth: u32, + } + let mut out: Vec = vec![front]; + let mut stack: Vec = Vec::with_capacity(64); + for i in (0..initial).rev() { + let (Some(pa), Some(pb)) = (<[Vec3]>::get(&ends, i as usize).copied(), <[Vec3]>::get(&ends, i as usize + 1).copied()) else { + return vec![front, back]; }; - if flat { - out_t.push(b); - } else { - let m = 0.5 * (a + b); - stack.push((m, b, depth + 1)); - stack.push((a, m, depth + 1)); + stack.push(Seg { a: t_at(i), b: t_at(i + 1), pa, pb, mid: None, depth: 0 }); + } + while let Some(sg) = stack.pop() { + let (a, b, pa, pb) = (sg.a, sg.b, sg.pa, sg.pb); + let m = 0.5 * (a + b); + let pm = match sg.mid { + Some(p) => Some(p), + None => at(m), + }; + let (q1, q3) = (at(a + 0.25 * (b - a)), at(a + 0.75 * (b - a))); + let done = sg.depth >= MAX_DEPTH || out.len() + stack.len() >= MAX_EDGE_NODES; + let fine = || { + let sag = [q1, pm, q3].iter().flatten().map(|q| q.dist_to_segment(pa, pb)).fold(0.0, f64::max); + let short = !(pa.dist(pb) > h * (1.0 + 1e-9)); + // The two half chords of an arc turning θ differ by θ/2. + let straight = match (rules.turn, pm) { + (Some(t), Some(q)) => !t.exceeds(2.0 * turn(q - pa, pb - q), sag), + _ => true, + }; + short && straight && sag <= tol + }; + match pm { + Some(pm) if !done && !fine() => { + stack.push(Seg { a: m, b, pa: pm, pb, mid: q3, depth: sg.depth + 1 }); + stack.push(Seg { a, b: m, pa, pb: pm, mid: q1, depth: sg.depth + 1 }); + } + _ => out.push(pb), } } - let mut pts: Vec = out_t.iter().filter_map(|&t| at(t)).collect(); - if pts.len() < 2 { + if out.len() < 2 { return vec![front, back]; } - if let Some(p) = pts.first_mut() { - *p = front; - } - if let Some(p) = pts.last_mut() { + if let Some(p) = out.last_mut() { *p = back; } - pts + out +} + +/// Angle (radians) between two directions; zero when either is empty. +fn turn(a: Vec3, b: Vec3) -> f64 { + let (la, lb) = (a.len(), b.len()); + if !(la > 0.0 && lb > 0.0) { + return 0.0; + } + (a.dot(b) / (la * lb)).clamp(-1.0, 1.0).acos() } #[cfg(test)] @@ -77,7 +130,7 @@ mod tests { let Ok(cyl) = crate::cylinder(Vec3::ZERO, Vec3::Z, 10.0, 5.0) else { panic!("cylinder") }; for tol in [0.1, 0.01, 0.001] { for e in cyl.solid.edge_iter() { - let pts = sample_edge(&e, tol); + let pts = sample_edge(&e, &Rules { tol, max_len: f64::INFINITY, turn: None }, f64::INFINITY); assert!(pts.len() >= 2); assert_eq!(pts.first().copied(), Some(from_p3(e.absolute_front().point()))); assert_eq!(pts.last().copied(), Some(from_p3(e.absolute_back().point()))); @@ -94,4 +147,34 @@ mod tests { } } } + + #[test] + fn edges_are_split_to_the_side_length_and_turn() { + let Ok(cyl) = crate::cylinder(Vec3::ZERO, Vec3::Z, 10.0, 5.0) else { panic!("cylinder") }; + for e in cyl.solid.edge_iter() { + let pts = sample_edge(&e, &Rules { tol: 1.0, max_len: f64::INFINITY, turn: None }, 1.2); + for w in pts.windows(2) { + assert!(w[0].dist(w[1]) <= 1.2 + 1e-9, "side {}", w[0].dist(w[1])); + } + if matches!(e.curve(), mt::Curve::Line(_)) { + // 5 mm in sides of at most 1.2: five equal sides. + assert_eq!(pts.len(), 6); + } + // The ends are the vertices' own points, to the last bit. + assert_eq!(pts.first().copied(), Some(from_p3(e.absolute_front().point()))); + assert_eq!(pts.last().copied(), Some(from_p3(e.absolute_back().point()))); + } + // A coarse tolerance alone would leave a few long chords; a 30° facet angle keeps every + // side within 30° of arc: at least 12 around each of the two circles. + let thirty = + Rules { tol: 5.0, max_len: f64::INFINITY, turn: Some(crate::tess::size::Turn { angle: std::f64::consts::PI / 6.0, floor: 0.0 }) }; + let mut seen = std::collections::HashSet::new(); + let sides: usize = cyl + .solid + .edge_iter() + .filter(|e| seen.insert(e.id()) && !matches!(e.curve(), mt::Curve::Line(_))) + .map(|e| sample_edge(&e, &thirty, f64::INFINITY).len() - 1) + .sum(); + assert!(sides >= 24, "{sides}"); + } } diff --git a/crates/kernel/src/tess/mod.rs b/crates/kernel/src/tess/mod.rs index 06509a1f..6516d4ee 100644 --- a/crates/kernel/src/tess/mod.rs +++ b/crates/kernel/src/tess/mod.rs @@ -20,6 +20,7 @@ pub(crate) mod cdt; pub(crate) mod coons; pub(crate) mod edges; pub(crate) mod grid; +pub(crate) mod size; pub(crate) mod surface; pub(crate) mod uv; @@ -40,6 +41,16 @@ pub struct PlanStats { pub delaunay: usize, } +/// Curved faces from which their sizes are worked out on the worker pool. +const PARALLEL_SIZES: usize = 8; +/// Curved faces from which faces are planned on the worker pool (measured on the oracle +/// parts: two curved faces plan faster serially, four break even, six run twice as fast). +const PARALLEL_FACES: usize = 4; +/// Intersection curves (a Newton solve per point), or edges of any kind, from which edges are +/// sampled on the worker pool. +const PARALLEL_CURVES: usize = 2; +const PARALLEL_EDGES: usize = 256; + /// Most lines a curved face's grid gets in one direction. const MAX_LINES: usize = 1024; /// Most points Ruppert refinement may add to one face, beyond its seeds. @@ -49,16 +60,34 @@ fn unsupported(what: impl std::fmt::Display) -> KernelError { KernelError::Failed(format!("not supported yet: tessellation: {what}")) } -/// Plan the display tessellation of `solid` with chord tolerance `tol`. -pub(crate) fn plan(solid: &Solid, tol: f64) -> Result<(TessPlan, PlanStats)> { +/// Plan the display tessellation of `solid` with chord tolerance `tol` and, when given, a +/// facet angle no facet may turn more than (where it sags above the angle's floor). +pub(crate) fn plan(solid: &Solid, tol: f64, turn: Option) -> Result<(TessPlan, PlanStats)> { let faces: Vec<&mt::Face> = solid.face_iter().collect(); let (edges, edge_faces) = edge_order(faces.iter().copied()); - // Every edge once, in parallel; faces then read their nodes. - let sampled = crate::par::map("edge nodes", &edges, 64, |e| Ok(edges::sample_edge(e, tol))); + let rules = size::Rules::new(tol, size::extent(solid.vertex_iter().map(|v| crate::body::from_p3(v.point()))), turn); + // Each face's surface and target side length, then every edge once at the smaller size of + // its faces; faces then read their nodes. A planar face costs microseconds and a curved + // one about a hundred, so the pool is used only where curved faces (or intersection + // curves) make up the work. + let evals: Vec> = faces.iter().map(|f| surface::eval_surface(&f.surface())).collect(); + let curved = evals.iter().filter(|e| !matches!(e, Some(EvalSurface::Plane { .. }))).count(); + let sizes = + crate::par::map("face sizes", &evals, curved >= PARALLEL_SIZES, |eval| Ok(eval.as_ref().map(|s| size::face_size(s, &uv::domain(s), &rules)))); + let sized: Vec<(Option, Option<(f64, [f64; 2])>)> = + evals.into_iter().zip(sizes).map(|(e, s)| s.map(|s| (e, s))).collect::>()?; + let face_h = |f: usize| sized.get(f).and_then(|(_, s)| s.map(|(h, _)| h)).unwrap_or(rules.max_len); + let edge_h: Vec = edge_faces.iter().map(|fs| fs.iter().map(|&f| face_h(f)).fold(rules.max_len, f64::min)).collect(); + let jobs: Vec<(&mt::Edge, f64)> = edges.iter().zip(edge_h).collect(); + let solved = edges.iter().filter(|e| matches!(e.curve(), mt::Curve::IntersectionCurve(_))).count(); + let sampled = crate::par::map("edge nodes", &jobs, solved >= PARALLEL_CURVES || edges.len() >= PARALLEL_EDGES, |(e, h)| { + Ok(edges::sample_edge(e, &rules, *h)) + }); let mut nodes: HashMap<_, Vec> = HashMap::with_capacity(edges.len()); for (e, pts) in edges.iter().zip(sampled) { nodes.insert(e.id(), pts?); } + relax_crowded_faces(&faces, &jobs, &rules, &mut nodes); let oriented = |e: &mt::Edge| -> Option> { let mut pts = nodes.get(&e.id())?.clone(); if !e.orientation() { @@ -67,10 +96,12 @@ pub(crate) fn plan(solid: &Solid, tol: f64) -> Result<(TessPlan, PlanStats)> { Some(pts) }; let indexed: Vec<(usize, &mt::Face)> = faces.iter().copied().enumerate().collect(); - let planned = crate::par::map("face tessellation", &indexed, 8, |(fi, face)| { + let planned = crate::par::map("face tessellation", &indexed, curved >= PARALLEL_FACES, |(fi, face)| { let fi = *fi; - let surface = face.surface(); - let eval = surface::eval_surface(&surface).ok_or_else(|| unsupported(format!("face {fi}: surface not evaluable")))?; + let (eval, size) = match sized.get(fi) { + Some((Some(eval), Some(size))) => (eval.clone(), *size), + _ => return Err(unsupported(format!("face {fi}: surface not evaluable"))), + }; // Boundary nodes, loop by loop (each edge's last node is the next edge's first). let mut loops3: Vec> = Vec::new(); for wire in face.absolute_boundaries() { @@ -85,7 +116,7 @@ pub(crate) fn plan(solid: &Solid, tol: f64) -> Result<(TessPlan, PlanStats)> { } } let mut stats = PlanStats::default(); - let fp = plan_face(fi, &eval, &loops3, tol, &mut stats)?; + let fp = plan_face(fi, &eval, &loops3, &rules, size, &mut stats)?; let fp = FacePlan { flip: !face.orientation(), triangles: if face.orientation() { fp.triangles } else { fp.triangles.into_iter().map(|[a, b, c]| [a, c, b]).collect() }, @@ -109,8 +140,56 @@ pub(crate) fn plan(solid: &Solid, tol: f64) -> Result<(TessPlan, PlanStats)> { Ok((plan, stats)) } +/// Bounded work for the facet angle, as the reference mesher bounds it: a face whose boundary +/// would get more than `MAX_ANGLE_BOUNDARY` segments has its curved edges sampled again at twice +/// the angle, as often as it takes, down to the tolerance alone (a plate with a thousand small +/// holes is not worth a hundred times the work to draw each hole round). An edge keeps one set +/// of nodes for both its faces, at the coarser of their levels, so every rim of a crowded face +/// gets the same share whatever the order of the faces. +fn relax_crowded_faces(faces: &[&mt::Face], jobs: &[(&mt::Edge, f64)], rules: &size::Rules, nodes: &mut HashMap>) { + use truck_meshalgo::tessellation::Deflection; + let Some(turn) = rules.turn else { return }; + let h_of: HashMap = jobs.iter().map(|(e, h)| (e.id(), (*e, *h))).collect(); + let boundaries: Vec> = + faces.iter().map(|f| f.absolute_boundaries().iter().flat_map(|w| w.edge_iter().map(|e| e.id())).collect()).collect(); + let mut level: HashMap = HashMap::new(); + for l in 1..=Deflection::MAX_RELAX { + let segments = |id: &mt::EdgeID| nodes.get(id).map_or(0, |n| n.len().saturating_sub(1)); + let mut crowded: Vec = boundaries + .iter() + .filter(|ids| ids.iter().map(segments).sum::() > Deflection::MAX_ANGLE_BOUNDARY) + .flatten() + .filter(|id| level.get(*id).copied().unwrap_or(0) < l) + .filter(|id| h_of.get(*id).is_some_and(|(e, _)| !matches!(e.curve(), mt::Curve::Line(_)))) + .copied() + .collect(); + crowded.sort_unstable_by_key(|id| format!("{id:?}")); + crowded.dedup(); + if crowded.is_empty() { + break; + } + let relaxed = size::Rules { turn: turn.relaxed(l), ..*rules }; + for id in crowded { + if let Some((e, h)) = h_of.get(&id) { + nodes.insert(id, edges::sample_edge(e, &relaxed, *h)); + level.insert(id, l); + } + } + } +} + /// One face: boundary loops in 3D (shared edge nodes) to parameters and triangles. -fn plan_face(fi: usize, s: &EvalSurface, loops3: &[Vec], tol: f64, stats: &mut PlanStats) -> Result { +/// `size`: the face's target side length (mm) and the surface's mean speed along u and v. +fn plan_face( + fi: usize, + s: &EvalSurface, + loops3: &[Vec], + rules: &size::Rules, + size: (f64, [f64; 2]), + stats: &mut PlanStats, +) -> Result { + let tol = rules.tol; + let (h, speed) = size; let dom = uv::domain(s); let accept = (tol * 10.0).max(1e-6); let fail = |why: &str| unsupported(format!("face {fi}: {why}")); @@ -172,26 +251,43 @@ fn plan_face(fi: usize, s: &EvalSurface, loops3: &[Vec], tol: f64, stats: if uvs.len() < 3 { return Err(fail("degenerate boundary")); } - let planar = matches!(s, EvalSurface::Plane { .. }); - let (extra, triangles) = if planar { - stats.planar += 1; - cdt::triangulate(&uvs, &ranges, cdt::Fill { seeds: &[], refine: None }).ok_or_else(|| fail("boundary crosses itself"))? - } else { - let mut bbox = [(f64::MAX, f64::MIN); 2]; - for p in &uvs { - for d in 0..2 { - bbox[d] = (bbox[d].0.min(p[d]), bbox[d].1.max(p[d])); - } + let mut bbox = [(f64::MAX, f64::MIN); 2]; + for p in &uvs { + for d in 0..2 { + bbox[d] = (bbox[d].0.min(p[d]), bbox[d].1.max(p[d])); } - let mut us = uv::divide(s, 0, bbox[0], bbox[1], tol, MAX_LINES); - let mut vs = uv::divide(s, 1, bbox[1], bbox[0], tol, MAX_LINES); + } + // Constrained Delaunay in (u, v) scaled to millimetres: equilateral seeds of side `h`, + // refined where triangles leave the surface, are too long or bend too much, then shaped by + // Ruppert refinement. + let delaunay = || { + let polys: Vec<&[[f64; 2]]> = ranges.iter().filter_map(|r| uvs.get(r.clone())).collect(); + let seeds = cdt::lattice(bbox, [h / speed[0], h / speed[1]], &polys, speed); + let eval = |u: f64, v: f64| s.derivs(u, v).map(|(p, su, sv)| (p, su.cross(sv).normalized())); + let sag = cdt::Sag { eval: &eval, tol, h, turn: rules.turn, scale: speed }; + let fill = cdt::Fill { seeds: &seeds, refine: Some((sag, MAX_REFINE)), shape: !matches!(s, EvalSurface::Plane { .. }) }; + cdt::triangulate(&uvs, &ranges, fill).ok_or_else(|| fail("boundary crosses itself")) + }; + // Planar and curved faces alike: a structured grid where the boundary runs along the + // parameter box (a rectangle, a cylinder half), a Coons grid on four corners, else Delaunay. + let planar = matches!(s, EvalSurface::Plane { .. }); + let (extra, triangles) = { + let mut us = uv::divide(s, 0, bbox[0], bbox[1], rules, h, MAX_LINES); + let mut vs = uv::divide(s, 1, bbox[1], bbox[0], rules, h, MAX_LINES); thin(&mut us, &mut vs, grid::MAX_FACE_GRID); let structured = match (ranges.len(), grid::quad(&uvs)) { (1, Some(q)) => { - // Interior lines: the boundary's own where opposite sides agree, else the - // tolerance division. + // Interior lines: the boundary's own where opposite sides agree (refined where the + // surface turns more between them inside the face), else the division. let lines = |d: usize, div: &Vec| { - grid::aligned_lines(&uvs, &q, d).unwrap_or_else(|| div.iter().copied().filter(|x| *x > bbox[d].0 && *x < bbox[d].1).collect()) + let inner = |l: Vec| l.into_iter().filter(|x| *x > bbox[d].0 && *x < bbox[d].1).collect::>(); + match grid::aligned_lines(&uvs, &q, d) { + Some(aligned) => { + let ends: Vec = std::iter::once(bbox[d].0).chain(aligned).chain(std::iter::once(bbox[d].1)).collect(); + inner(uv::refine(s, d, &ends, bbox[1 - d], rules, h, MAX_LINES)) + } + None => inner(div.clone()), + } }; let (lu, lv) = (lines(0, &us), lines(1, &vs)); let strip = match (lu.is_empty(), lv.is_empty()) { @@ -203,28 +299,29 @@ fn plan_face(fi: usize, s: &EvalSurface, loops3: &[Vec], tol: f64, stats: } _ => None, }; - let cells = |lines: &Vec, (a, b): (f64, f64)| (lines.len().saturating_sub(1).max(1) as f64) / (b - a).max(f64::MIN_POSITIVE); - let scale = [cells(&us, bbox[0]), cells(&vs, bbox[1])]; let structured = match structured { - None if ranges.len() == 1 => coons::coons(&uvs, &fixed, scale, grid::MAX_FACE_GRID), + None if ranges.len() == 1 => { + let fine = |a: [f64; 2], b: [f64; 2]| uv::fine_between(s, rules, h, a, b); + coons::coons(&uvs, &fixed, speed, grid::MAX_FACE_GRID, &fine) + } other => other, }; match structured { Some(r) => { - stats.structured += 1; + if planar { + stats.planar += 1; + } else { + stats.structured += 1; + } r } None => { - stats.delaunay += 1; - let polys: Vec<&[[f64; 2]]> = ranges.iter().filter_map(|r| uvs.get(r.clone())).collect(); - let inner = |lines: &Vec, (a, b): (f64, f64)| lines.iter().copied().filter(|x| *x > a && *x < b).collect::>(); - let seeds = cdt::inside_grid(&inner(&us, bbox[0]), &inner(&vs, bbox[1]), &polys); - let eval = |u: f64, v: f64| s.derivs(u, v).map(|(p, _, _)| p); - // One unit per tolerance cell in each direction: a straight direction (one - // cell) keeps long triangles, Ruppert shapes them in that metric. - let sag = cdt::Sag { eval: &eval, tol, scale }; - let fill = cdt::Fill { seeds: &seeds, refine: Some((sag, MAX_REFINE)) }; - cdt::triangulate(&uvs, &ranges, fill).ok_or_else(|| fail("boundary crosses itself"))? + if planar { + stats.planar += 1; + } else { + stats.delaunay += 1; + } + delaunay()? } } }; diff --git a/crates/kernel/src/tess/size.rs b/crates/kernel/src/tess/size.rs new file mode 100644 index 00000000..37129878 --- /dev/null +++ b/crates/kernel/src/tess/size.rs @@ -0,0 +1,231 @@ +//! Target triangle size: one side length per face, so triangles come out about equilateral +//! (well-shaped triangles shade evenly) and as small as the surface needs. +//! +//! A face's size is the smallest of the body's longest side, the length over which the surface +//! normal turns the caller's facet angle (when one is given, see [`Turn`]), and the chord that +//! sags the tolerance; an edge takes the smaller size of its two faces, so a curved face and +//! its flat neighbour meet at the same spacing. + +use solidcraft_geom::Vec3; +use solidcraft_geom::tess::EvalSurface; + +use super::uv::Domain; + +/// Longest side, as a fraction of the body's extent. +const MAX_LEN_FRACTION: f64 = 1.0 / 5.0; +/// Longest side relative to the side the curvature asks for: along a straight direction (a +/// cylinder's axis) triangles may be this many times longer than across it. With a 30° facet +/// angle and [`MAX_LEN_FRACTION`] this puts the oracle parts' display meshes at about the +/// reference mesher's triangle count (1.2×; 5× with equilateral triangles and 15°). +const ASPECT: f64 = 4.0; +/// Shortest side the aspect limit may ask for, as a fraction of the longest side. A bend whose +/// radius is tiny next to the body (a hem) would otherwise tie the whole length of its +/// straight edges to that radius; along a straight direction the normal does not turn, so +/// long triangles there shade exactly, and the turn and sag limits still divide the bend +/// across. +const MIN_LEN_FRACTION: f64 = 1.0 / 4.0; +/// Curvature samples per parameter direction. +const SAMPLES: usize = 9; + +/// The most one facet may turn (radians): an edge segment's tangent or a cell's surface normal. +/// The angle divides only where a facet would sag more than `floor` (mm), so a radius far below +/// the tolerance keeps its chord facets. +#[derive(Clone, Copy, Debug)] +pub(crate) struct Turn { + pub angle: f64, + pub floor: f64, +} + +impl Turn { + /// Whether a facet that turns `angle` and sags `sag` from the surface must be divided. + pub(crate) fn exceeds(&self, angle: f64, sag: f64) -> bool { + sag > self.floor && angle > self.angle + } + + /// The angle doubled `level` times, as the reference mesher relaxes it for a crowded face; + /// from its last level, or past half a turn, no angle at all (the tolerance alone). + pub(crate) fn relaxed(&self, level: u8) -> Option { + use truck_meshalgo::tessellation::Deflection; + let angle = self.angle * f64::from(1u32 << u32::from(level.min(Deflection::MAX_RELAX))); + (level == 0 || (level < Deflection::MAX_RELAX && angle < std::f64::consts::PI)).then_some(Turn { angle, ..*self }) + } +} + +/// The chord tolerance, the longest side and the facet angle of one body. +#[derive(Clone, Copy, Debug)] +pub(crate) struct Rules { + pub tol: f64, + pub max_len: f64, + pub turn: Option, +} + +impl Rules { + /// Rules for a body of extent `extent` (see [`extent`]). + pub(crate) fn new(tol: f64, extent: f64, turn: Option) -> Rules { + let max_len = if extent.is_finite() && extent > 0.0 { extent * MAX_LEN_FRACTION } else { f64::INFINITY }; + Rules { tol, max_len: max_len.max(tol), turn } + } + + /// The side length for a curvature `k` (1/mm): the chord sags at most the tolerance across + /// it, and the normal turns at most the facet angle where that turn sags above its floor. + fn for_curvature(&self, k: f64) -> f64 { + if !(k.is_finite() && k > 0.0) { + return self.max_len; + } + let r = 1.0 / k; + let mut across = (8.0 * r * self.tol).sqrt(); + if let Some(t) = self.turn + && r * (1.0 - (0.5 * t.angle).cos()) > t.floor + { + across = across.min(r * t.angle); + } + let aspect = ASPECT * across; + self.max_len.min(aspect.max(self.max_len * MIN_LEN_FRACTION)).max(self.tol) + } +} + +/// The body's extent: twice the largest distance of a vertex from the vertices' centroid. Unlike +/// a bounding box diagonal it does not change when the body is turned, so a placed copy meshes +/// like the original. +pub(crate) fn extent(points: impl Iterator) -> f64 { + let points: Vec = points.filter(|p| p.is_finite()).collect(); + if points.is_empty() { + return 0.0; + } + let c = points.iter().fold(Vec3::ZERO, |a, p| a + *p) * (1.0 / points.len() as f64); + 2.0 * points.iter().map(|p| p.dist(c)).fold(0.0, f64::max) +} + +/// A face's target side length and the surface's mean speed along u and v (mm per parameter +/// unit), from samples over the surface's parameter box (the whole box of a trimmed surface, +/// so the size is never larger than the face needs). +pub(crate) fn face_size(s: &EvalSurface, dom: &Domain, rules: &Rules) -> (f64, [f64; 2]) { + if let EvalSurface::Plane { du, dv, .. } = s { + return (rules.max_len, [du.len(), dv.len()]); + } + let (Some(ru), Some(rv)) = (dom.range[0], dom.range[1]) else { return (rules.max_len, [1.0, 1.0]) }; + // A flat spline (sheet-metal walls, chamfers): no curvature to look for, only the speeds. + let samples = if planar(s) { 3 } else { SAMPLES }; + let normal = |u: f64, v: f64| -> Option<(Vec3, Vec3, [f64; 2])> { + let (p, su, sv) = s.derivs(u, v)?; + Some((p, su.cross(sv).normalized()?, [su.len(), sv.len()])) + }; + let step = [(ru.1 - ru.0) * 1e-3, (rv.1 - rv.0) * 1e-3]; + let mut k_max: f64 = 0.0; + let (mut speed, mut count) = ([0.0; 2], 0usize); + for i in 0..samples { + for j in 0..samples { + // Cell centres, clear of the box edges (poles, seams). + let u = ru.0 + (ru.1 - ru.0) * (i as f64 + 0.5) / samples as f64; + let v = rv.0 + (rv.1 - rv.0) * (j as f64 + 0.5) / samples as f64; + let Some((p, n, sp)) = normal(u, v) else { continue }; + speed = [speed[0] + sp[0], speed[1] + sp[1]]; + count += 1; + if samples < SAMPLES { + continue; + } + // How fast the normal turns per mm, along u and along v (twist included). + for (du, dv) in [(step[0], 0.0), (0.0, step[1])] { + let Some((q, m, _)) = normal(u + du, v + dv) else { continue }; + let ds = p.dist(q); + if ds > 0.0 { + k_max = k_max.max((m - n).len() / ds); + } + } + } + } + let speed = if count > 0 { [speed[0] / count as f64, speed[1] / count as f64] } else { [1.0, 1.0] }; + (rules.for_curvature(k_max), speed.map(|x| if x.is_finite() && x > 0.0 { x } else { 1.0 })) +} + +/// A spline surface whose control points lie in one plane, with positive weights: the +/// surface lies in that plane (it is inside the control points' convex hull). +fn planar(s: &EvalSurface) -> bool { + let EvalSurface::Nurbs { ctrl, .. } = s else { return false }; + if ctrl.iter().any(|h| !(h[3] > 0.0)) { + return false; + } + let pts: Vec = ctrl.iter().map(|h| Vec3::new(h[0] / h[3], h[1] / h[3], h[2] / h[3])).collect(); + let Some(&a) = pts.first() else { return false }; + let size = pts.iter().map(|p| p.dist(a)).fold(0.0, f64::max); + // The two points that span the widest triangle with the first give the plane. + let Some(&b) = pts.iter().max_by(|p, q| p.dist(a).total_cmp(&q.dist(a))) else { return false }; + let Some(n) = pts.iter().map(|p| (b - a).cross(*p - a)).max_by(|x, y| x.len().total_cmp(&y.len())).and_then(Vec3::normalized) else { + return false; + }; + size > 0.0 && pts.iter().all(|p| (*p - a).dot(n).abs() <= 1e-9 * size) +} + +#[cfg(test)] +mod tests { + use super::*; + use solidcraft_geom::tess::EvalCurve; + + const THIRTY: Option = Some(Turn { angle: std::f64::consts::PI / 6.0, floor: 0.0 }); + + #[test] + fn a_cylinder_is_sized_by_its_radius() { + let s = EvalSurface::Revolved { + profile: EvalCurve { degree: 1, knots: vec![0.0, 0.0, 1.0, 1.0], ctrl: vec![[3.0, 0.0, 0.0, 1.0], [3.0, 0.0, 40.0, 1.0]] }, + origin: Vec3::ZERO, + axis: Vec3::Z, + transform: [[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0]], + swap: false, + }; + // A 20 mm body: sides from 0.5 mm (the floor) to 4 mm. + let rules = Rules::new(0.01, 20.0, THIRTY); + let (h, speed) = face_size(&s, &super::super::uv::domain(&s), &rules); + // r = 3: 30° of turn is 1.57 mm; the 0.01 mm sag chord is √(8·3·0.01) = 0.49 mm. + assert!((h - ASPECT * (8.0 * 3.0 * 0.01f64).sqrt()).abs() < 0.01, "{h}"); + assert!((speed[0] - 40.0).abs() < 1e-6 && (speed[1] - 3.0).abs() < 1e-6, "{speed:?}"); + // Coarser: the turn limit governs. + let (h, _) = face_size(&s, &super::super::uv::domain(&s), &Rules::new(0.5, 100.0, THIRTY)); + assert!((h - ASPECT * 3.0 * std::f64::consts::PI / 6.0).abs() < 0.01, "{h}"); + // Without a facet angle the sag alone: √(8·3·0.5) = 3.46 mm, four times along. + let (h, _) = face_size(&s, &super::super::uv::domain(&s), &Rules::new(0.5, 100.0, None)); + assert!((h - ASPECT * (8.0 * 3.0 * 0.5f64).sqrt()).abs() < 0.01, "{h}"); + // A turn that sags below its floor does not divide. + let floored = Some(Turn { angle: std::f64::consts::PI / 6.0, floor: 1.0 }); + let (h, _) = face_size(&s, &super::super::uv::domain(&s), &Rules::new(0.5, 100.0, floored)); + assert!((h - ASPECT * (8.0 * 3.0 * 0.5f64).sqrt()).abs() < 0.01, "{h}"); + } + + #[test] + fn a_flat_spline_is_planar_and_a_bent_one_is_not() { + let patch = |z: f64| EvalSurface::Nurbs { + degree: [1, 1], + knots: [vec![0.0, 0.0, 1.0, 1.0], vec![0.0, 0.0, 1.0, 1.0]], + count: [2, 2], + ctrl: vec![[0.0, 0.0, 0.0, 1.0], [4.0, 0.0, 0.0, 1.0], [0.0, 3.0, 0.0, 1.0], [4.0, 3.0, z, 1.0]], + }; + assert!(planar(&patch(0.0))); + assert!(!planar(&patch(0.5))); + let rules = Rules::new(0.01, 100.0, THIRTY); + let (h, speed) = face_size(&patch(0.0), &super::super::uv::domain(&patch(0.0)), &rules); + assert_eq!(h, rules.max_len); + assert!((speed[0] - 4.0).abs() < 1e-9 && (speed[1] - 3.0).abs() < 1e-9, "{speed:?}"); + // The twisted patch (a hyperbolic paraboloid) has curvature: a smaller side. + assert!(face_size(&patch(0.5), &super::super::uv::domain(&patch(0.5)), &rules).0 < rules.max_len); + } + + #[test] + fn a_tiny_radius_does_not_shrink_the_side_below_the_floor() { + let s = EvalSurface::Revolved { + profile: EvalCurve { degree: 1, knots: vec![0.0, 0.0, 1.0, 1.0], ctrl: vec![[0.05, 0.0, 0.0, 1.0], [0.05, 0.0, 80.0, 1.0]] }, + origin: Vec3::ZERO, + axis: Vec3::Z, + transform: [[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0]], + swap: false, + }; + let rules = Rules::new(0.1, 100.0, THIRTY); + let (h, _) = face_size(&s, &super::super::uv::domain(&s), &rules); + assert!((h - rules.max_len * MIN_LEN_FRACTION).abs() < 1e-9, "{h}"); + } + + #[test] + fn a_plane_takes_the_longest_side() { + let s = EvalSurface::Plane { origin: Vec3::ZERO, du: Vec3::X, dv: Vec3::Y }; + let rules = Rules::new(0.01, 100.0, THIRTY); + assert_eq!(face_size(&s, &super::super::uv::domain(&s), &rules).0, 20.0); + } +} diff --git a/crates/kernel/src/tess/tests.rs b/crates/kernel/src/tess/tests.rs index 8b756cfb..6580299c 100644 --- a/crates/kernel/src/tess/tests.rs +++ b/crates/kernel/src/tess/tests.rs @@ -27,7 +27,7 @@ fn assert_watertight(m: &Mesh, name: &str) { } fn plan_body(b: &Body, tol: f64) -> (Mesh, PlanStats) { - match plan(&b.solid, tol) { + match plan(&b.solid, tol, Some(crate::tess::size::Turn { angle: crate::DISPLAY_ANGLE, floor: tol * 1e-3 })) { Ok((p, stats)) => (p.evaluate(), stats), Err(e) => panic!("plan: {e}"), } @@ -63,6 +63,33 @@ fn cut(a: &Body, b: &Body) -> Body { const PI: f64 = std::f64::consts::PI; +/// Smallest corner angle (degrees) over the mesh. +fn min_angle(m: &Mesh) -> f64 { + let mut least = 180.0f64; + for t in &m.triangles { + let [a, b, c] = t.map(|i| m.positions[i as usize]); + for (p, q, r) in [(a, b, c), (b, c, a), (c, a, b)] { + let (u, v) = (q - p, r - p); + if let (Some(u), Some(v)) = (u.normalized(), v.normalized()) { + least = least.min(u.dot(v).clamp(-1.0, 1.0).acos().to_degrees()); + } + } + } + least +} + +#[test] +fn a_box_is_a_regular_grid_on_every_face() { + // 40 × 30 × 20: sides at most 1/5 of the diagonal (10.77 mm), so 4 × 3 × 2 cells. + let b = ok(crate::box_solid(Vec3::ZERO, Vec3::new(40.0, 30.0, 20.0))); + let (m, s) = plan_body(&b, 0.05); + assert_eq!(s.planar, 6); + assert_watertight(&m, "box"); + let cells = 2 * (4 * 3 + 4 * 2 + 3 * 2); + assert_eq!(m.triangles.len(), 2 * cells, "two triangles per grid cell"); + assert!(min_angle(&m) > 40.0, "{}", min_angle(&m)); +} + #[test] fn primitives_are_watertight_with_the_right_volume() { for tol in [0.05, 0.005] { @@ -148,7 +175,7 @@ fn render_wireframes() { let Ok(old) = body.tessellate(tol) else { panic!("reference") }; let t_old = t0.elapsed().as_secs_f64() * 1e3; let t0 = std::time::Instant::now(); - let new = match plan(&body.solid, tol) { + let new = match plan(&body.solid, tol, Some(crate::tess::size::Turn { angle: crate::DISPLAY_ANGLE, floor: tol * 1e-3 })) { Ok((p, st)) => Some((p.evaluate(), st)), Err(e) => { println!("{name}: no plan ({e})"); @@ -243,3 +270,12 @@ fn a_drill_point_apex_leaves_no_zero_area_triangles() { } } } + +#[test] +fn a_placed_body_meshes_like_the_original() { + let b = ok(crate::box_solid(Vec3::ZERO, Vec3::new(2.0, 3.0, 4.0))); + let placed = ok(crate::transform(&b, Vec3::new(15.0, -2.0, 8.0), Vec3::ZERO, Vec3::X, 0.7)); + let (m, _) = plan_body(&b, 0.005); + let (p, _) = plan_body(&placed, 0.005); + assert_eq!(m.tri_face, p.tri_face); +} diff --git a/crates/kernel/src/tess/uv.rs b/crates/kernel/src/tess/uv.rs index 2b2457fb..b495aedd 100644 --- a/crates/kernel/src/tess/uv.rs +++ b/crates/kernel/src/tess/uv.rs @@ -5,6 +5,8 @@ use solidcraft_geom::Vec3; use solidcraft_geom::tess::EvalSurface; +use super::size::Rules; + /// The parameter box of a surface and its periods. #[derive(Clone, Copy, Debug)] pub(crate) struct Domain { @@ -144,36 +146,43 @@ pub(crate) fn degenerate(s: &EvalSurface, uv: [f64; 2], d: usize) -> bool { this.len() <= 1e-9 * other.len().max(1e-12) } -/// Division of `range` along direction `d` so that the surface between neighbouring lines -/// sags at most `tol` from its chord, checked along iso-lines across `across`; at most -/// `max_lines` lines, ends included. -pub(crate) fn divide(s: &EvalSurface, d: usize, range: (f64, f64), across: (f64, f64), tol: f64, max_lines: usize) -> Vec { +/// Division of `range` along direction `d` so that between neighbouring lines the surface sags +/// at most the tolerance from its chord, the chord is at most `h` long and the normal turns +/// within the rules' facet angle, checked along iso-lines across `across`; at most `max_lines` +/// lines, ends included. +pub(crate) fn divide(s: &EvalSurface, d: usize, range: (f64, f64), across: (f64, f64), rules: &Rules, h: f64, max_lines: usize) -> Vec { let (a0, a1) = range; - if matches!(s, EvalSurface::Plane { .. }) || !(a1 > a0) { + if !(a1 > a0) { return vec![a0, a1]; } - let at = |x: f64, y: f64| { - let (u, v) = if d == 0 { (x, y) } else { (y, x) }; - s.derivs(u, v).map(|(p, _, _)| p) - }; - const ACROSS: usize = 5; - let ys: Vec = (0..ACROSS).map(|k| across.0 + (across.1 - across.0) * k as f64 / (ACROSS - 1) as f64).collect(); - let flat = |a: f64, b: f64| { - ys.iter().all(|&y| { - let (Some(pa), Some(pb)) = (at(a, y), at(b, y)) else { return true }; - [0.25, 0.5, 0.75].iter().all(|f| at(a + (b - a) * f, y).is_none_or(|q| q.dist_to_segment(pa, pb) <= tol)) - }) - }; - // Start from quarter turns at most, so a closed iso-line is not taken for a flat one. - let initial = ((a1 - a0) / std::f64::consts::FRAC_PI_2).ceil().clamp(1.0, 64.0) as u32; + let flat = limits(s, d, across, rules, h); + // Start from quarter turns at most, so a closed iso-line is not taken for a flat one (a + // plane's parameters are lengths and its lines are never closed). + let initial = if matches!(s, EvalSurface::Plane { .. }) { 1 } else { ((a1 - a0) / std::f64::consts::FRAC_PI_2).ceil().clamp(1.0, 64.0) as u32 }; // Both ends exactly (a0 + (a1 − a0) · 1 may round away from a1): the end lines are the // face's own parameter box, which the boundary nodes sit on. - let x_at = |i: u32| { - let t = f64::from(i) / f64::from(initial); - a0 * (1.0 - t) + a1 * t - }; - let mut out = vec![a0]; - let mut stack: Vec<(f64, f64, u32)> = (0..initial).rev().map(|i| (x_at(i), x_at(i + 1), 0)).collect(); + let starts: Vec = (0..=initial) + .map(|i| { + let t = f64::from(i) / f64::from(initial); + a0 * (1.0 - t) + a1 * t + }) + .collect(); + split(&starts, &flat, max_lines) +} + +/// `lines` (increasing, ends included) with every interval that fails the same limits as +/// [`divide`] split until it passes: lines taken from a face's edge nodes meet the edges node for +/// node, but nothing else checks how far the surface turns between them inside the face. +pub(crate) fn refine(s: &EvalSurface, d: usize, lines: &[f64], across: (f64, f64), rules: &Rules, h: f64, max_lines: usize) -> Vec { + let flat = limits(s, d, across, rules, h); + split(lines, &flat, max_lines) +} + +/// Depth-first bisection of the intervals between `starts` until `flat` holds for each. +fn split(starts: &[f64], flat: &dyn Fn(f64, f64) -> bool, max_lines: usize) -> Vec { + let Some(&first) = starts.first() else { return Vec::new() }; + let mut out = vec![first]; + let mut stack: Vec<(f64, f64, u32)> = starts.windows(2).rev().map(|w| (w[0], w[1], 0)).collect(); while let Some((a, b, depth)) = stack.pop() { if depth >= 20 || out.len() + stack.len() >= max_lines.max(2) || flat(a, b) { out.push(b); @@ -186,6 +195,66 @@ pub(crate) fn divide(s: &EvalSurface, d: usize, range: (f64, f64), across: (f64, out } +/// Whether the surface between the parameter points `a` and `b` is within the tolerance (at +/// their middle), the side length `h` and the facet angle: no point is needed between them. +pub(crate) fn fine_between(s: &EvalSurface, rules: &Rules, h: f64, a: [f64; 2], b: [f64; 2]) -> bool { + let at = |q: [f64; 2]| s.derivs(q[0], q[1]); + let (Some((pa, _, _)), Some((pb, _, _))) = (at(a), at(b)) else { return true }; + if pa.dist(pb) > h * (1.0 + 1e-9) { + return false; + } + let mid = [0.5 * (a[0] + b[0]), 0.5 * (a[1] + b[1])]; + let sag = at(mid).map_or(0.0, |(m, _, _)| m.dist((pa + pb) * 0.5)); + if sag > rules.tol { + return false; + } + // At a pole the normal is undefined: take it from just inside the segment. + let normal = |q: [f64; 2], toward: [f64; 2]| { + let n = |q: [f64; 2]| at(q).and_then(|(_, su, sv)| su.cross(sv).normalized()); + n(q).or_else(|| n([q[0] + 1e-4 * (toward[0] - q[0]), q[1] + 1e-4 * (toward[1] - q[1])])) + }; + match (rules.turn, normal(a, b), normal(b, a)) { + (Some(t), Some(na), Some(nb)) => !t.exceeds(na.dot(nb).clamp(-1.0, 1.0).acos(), sag), + _ => true, + } +} + +/// Whether the surface between the iso-lines at `a` and `b` (direction `d`, sampled along +/// iso-lines across `across`) is within the tolerance, the side length `h` and the facet angle. +fn limits<'a>(s: &'a EvalSurface, d: usize, across: (f64, f64), rules: &'a Rules, h: f64) -> impl Fn(f64, f64) -> bool + 'a { + let tol = rules.tol; + let uv = move |x: f64, y: f64| if d == 0 { (x, y) } else { (y, x) }; + let at = move |x: f64, y: f64| { + let (u, v) = uv(x, y); + s.derivs(u, v).map(|(p, _, _)| p) + }; + let normal = move |x: f64, y: f64| { + let (u, v) = uv(x, y); + s.derivs(u, v).and_then(|(_, su, sv)| su.cross(sv).normalized()) + }; + const ACROSS: usize = 5; + let ys: Vec = (0..ACROSS).map(|k| across.0 + (across.1 - across.0) * k as f64 / (ACROSS - 1) as f64).collect(); + move |a: f64, b: f64| { + ys.iter().all(|&y| { + let (Some(pa), Some(pb)) = (at(a, y), at(b, y)) else { return true }; + if pa.dist(pb) > h * (1.0 + 1e-9) { + return false; + } + let sag = [0.25, 0.5, 0.75].iter().filter_map(|f| at(a + (b - a) * f, y)).map(|q| q.dist_to_segment(pa, pb)).fold(0.0, f64::max); + if sag > tol { + return false; + } + // At a pole the normal is undefined: take it from just inside the segment. + let na = normal(a, y).or_else(|| normal(a + 1e-4 * (b - a), y)); + let nb = normal(b, y).or_else(|| normal(b - 1e-4 * (b - a), y)); + match (rules.turn, na, nb) { + (Some(t), Some(na), Some(nb)) => !t.exceeds(na.dot(nb).clamp(-1.0, 1.0).acos(), sag), + _ => true, + } + }) + } +} + #[cfg(test)] mod tests { use super::*; @@ -219,11 +288,16 @@ mod tests { #[test] fn divides_the_turn_but_not_the_straight_direction() { let s = cylinder(); - let turn = divide(&s, 1, (0.0, std::f64::consts::TAU), (0.0, 1.0), 0.01, 10_000); - let along = divide(&s, 0, (0.0, 1.0), (0.0, std::f64::consts::TAU), 0.01, 10_000); + let turn = + divide(&s, 1, (0.0, std::f64::consts::TAU), (0.0, 1.0), &Rules { tol: 0.01, max_len: f64::INFINITY, turn: None }, f64::INFINITY, 10_000); + let along = + divide(&s, 0, (0.0, 1.0), (0.0, std::f64::consts::TAU), &Rules { tol: 0.01, max_len: f64::INFINITY, turn: None }, f64::INFINITY, 10_000); assert_eq!(along, vec![0.0, 1.0]); // Chord sag r(1 − cos(θ/2)) ≤ 0.01 with r = 3 needs θ ≤ 0.163: at least 39 segments. assert!(turn.len() >= 40 && turn.len() < 80, "{}", turn.len()); + // With a longest side the straight direction is divided too: 4 mm in sides ≤ 0.5. + let along = divide(&s, 0, (0.0, 1.0), (0.0, std::f64::consts::TAU), &Rules { tol: 0.01, max_len: f64::INFINITY, turn: None }, 0.5, 10_000); + assert_eq!(along.len(), 9); } /// The first and last lines are the range's ends to the last bit: the boundary nodes sit on @@ -232,7 +306,7 @@ mod tests { #[test] fn divides_with_the_ends_exactly() { let (a0, a1) = (-1.010485317689877, 1.083909784703316); - let lines = divide(&cylinder(), 1, (a0, a1), (0.0, 1.0), 0.01, 10_000); + let lines = divide(&cylinder(), 1, (a0, a1), (0.0, 1.0), &Rules { tol: 0.01, max_len: f64::INFINITY, turn: None }, f64::INFINITY, 10_000); assert_eq!(lines.first().copied(), Some(a0)); assert_eq!(lines.last().copied(), Some(a1)); }