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 <noreply@anthropic.com>

Cargo.lock

@@ -3365,6 +3365,7 @@ name = "solidcraft-kernel"
 version = "0.8.0"
 dependencies = [
  "log",
+ "rayon",
  "serde",
  "serde_json",
  "solidcraft-geom",

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"] }

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 }

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<solidcraft_geom::tess::TessPlan> {
+        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<solidcraft_geom::tess::TessPlan> {
+        self.plan_within(tol, Some(angle))
+    }
+
+    fn plan_within(&self, tol: f64, angle: Option<f64>) -> Result<solidcraft_geom::tess::TessPlan> {
         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).

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<T: Sync, R: Send>(what: &'static str, items: &[T], min_items: usize, f: impl Fn(&T) -> Result<R> + Sync) -> Vec<Result<R>> {
+/// `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<T: Sync, R: Send>(what: &'static str, items: &[T], parallel: bool, f: impl Fn(&T) -> Result<R> + Sync) -> Vec<Result<R>> {
     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::<Vec<_>>()).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<T: Sync, R: Send>(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<u32> = (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 {

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<Point2<f64>>;
 
 /// 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<Vec3>,
+    /// Point and unit normal (when not degenerate) at (u, v).
+    pub eval: &'a dyn Fn(f64, f64) -> Option<(Vec3, Option<Vec3>)>,
     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<Turn>,
     /// 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<usize>
         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<f64>| [p.x / scale[0], p.y / scale[1]];
     let mut cdt = Cdt::new();
     let mut index: HashMap<FixedVertexHandle, u32> = HashMap::new();
@@ -77,12 +96,14 @@ pub(crate) fn triangulate(boundary: &[[f64; 2]], loops: &[std::ops::Range<usize>
     if let Some((sag, budget)) = &fill.refine {
         let constraints = cdt.num_constraints();
         let added = refine_sag(&mut cdt, sag, *budget, &from);
+        if fill.shape {
             let params = RefinementParameters::<f64>::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<usize>
     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>) -> [f64; 2]) -> usize {
     let mut added = 0;
-    let mut cache: HashMap<FixedVertexHandle, Option<Vec3>> = HashMap::new();
+    let mut cache: HashMap<FixedVertexHandle, Option<(Vec3, Option<Vec3>)>> = HashMap::new();
+    // The widest turn between two corner normals.
+    let spread = |n: [Option<Vec3>; 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<Point2<f64>> = Vec::new();
         let mut at = |h: FixedVertexHandle, p: Point2<f64>| {
@@ -136,12 +163,17 @@ fn refine_sag(cdt: &mut Cdt, sag: &Sag, budget: usize, from: &dyn Fn(Point2<f64>
             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<f64>, 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<f64>, chord: Vec3| {
+                let uv = from(q);
+                (sag.eval)(uv[0], uv[1]).map_or(0.0, |(s, _)| s.dist(chord))
+            };
+            let off = |q: Point2<f64>, 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<f64>
                 }
                 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<f64>
             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<bool> {
     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<f64> = 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);
     }
 }

crates/kernel/src/tess/coons.rs

@@ -88,7 +88,15 @@ fn lines(a: &Side, b: &Side) -> Vec<f64> {
 
 /// 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));
-    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 {
+    // 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])
         };
-            out.push([blend(0), blend(1)]);
+        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 mut out: Vec<[f64; 2]> = uv.to_vec();
+    let base = u32::try_from(n).ok()?;
+    for &t in &ts {
+        for &s in &ss {
+            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<f64>, max: usize, fine: &dyn Fn(f64, f64, f64) -> Option<bool>) -> Vec<f64> {
+    const ACROSS: [f64; 5] = [0.0, 0.25, 0.5, 0.75, 1.0];
+    let all: Vec<f64> = 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<Vec3> = 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<Vec3> = 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());
     }
 }

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<Vec3> {
+/// 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<Vec3> {
+    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<Vec3> = (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<Vec3> {
         }
         .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<Vec3> = 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<Vec3>,
+        depth: u32,
+    }
+    let mut out: Vec<Vec3> = vec![front];
+    let mut stack: Vec<Seg> = 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 {
+        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);
-            stack.push((m, b, depth + 1));
-            stack.push((a, m, depth + 1));
+        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<Vec3> = 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}");
+    }
 }

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<size::Turn>) -> 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<Option<EvalSurface>> = 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<EvalSurface>, Option<(f64, [f64; 2])>)> =
+        evals.into_iter().zip(sizes).map(|(e, s)| s.map(|s| (e, s))).collect::<Result<_>>()?;
+    let face_h = |f: usize| sized.get(f).and_then(|(_, s)| s.map(|(h, _)| h)).unwrap_or(rules.max_len);
+    let edge_h: Vec<f64> = 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<Vec3>> = 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<Vec<Vec3>> {
         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<Vec3>> = 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<mt::EdgeID, Vec<Vec3>>) {
+    use truck_meshalgo::tessellation::Deflection;
+    let Some(turn) = rules.turn else { return };
+    let h_of: HashMap<mt::EdgeID, (&mt::Edge, f64)> = jobs.iter().map(|(e, h)| (e.id(), (*e, *h))).collect();
+    let boundaries: Vec<Vec<mt::EdgeID>> =
+        faces.iter().map(|f| f.absolute_boundaries().iter().flat_map(|w| w.edge_iter().map(|e| e.id())).collect()).collect();
+    let mut level: HashMap<mt::EdgeID, u8> = 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<mt::EdgeID> = boundaries
+            .iter()
+            .filter(|ids| ids.iter().map(segments).sum::<usize>() > 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<Vec3>], tol: f64, stats: &mut PlanStats) -> Result<FacePlan> {
+/// `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<Vec3>],
+    rules: &size::Rules,
+    size: (f64, [f64; 2]),
+    stats: &mut PlanStats,
+) -> Result<FacePlan> {
+    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<Vec3>], 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 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<f64>| {
-                    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<f64>| l.into_iter().filter(|x| *x > bbox[d].0 && *x < bbox[d].1).collect::<Vec<f64>>();
+                    match grid::aligned_lines(&uvs, &q, d) {
+                        Some(aligned) => {
+                            let ends: Vec<f64> = 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<Vec3>], tol: f64, stats:
             }
             _ => None,
         };
-        let cells = |lines: &Vec<f64>, (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) => {
+                if planar {
+                    stats.planar += 1;
+                } else {
                     stats.structured += 1;
+                }
                 r
             }
             None => {
+                if planar {
+                    stats.planar += 1;
+                } else {
                     stats.delaunay += 1;
-                let polys: Vec<&[[f64; 2]]> = ranges.iter().filter_map(|r| uvs.get(r.clone())).collect();
-                let inner = |lines: &Vec<f64>, (a, b): (f64, f64)| lines.iter().copied().filter(|x| *x > a && *x < b).collect::<Vec<f64>>();
-                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"))?
+                }
+                delaunay()?
             }
         }
     };

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<Turn> {
+        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<Turn>,
+}
+
+impl Rules {
+    /// Rules for a body of extent `extent` (see [`extent`]).
+    pub(crate) fn new(tol: f64, extent: f64, turn: Option<Turn>) -> 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<Item = Vec3>) -> f64 {
+    let points: Vec<Vec3> = 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<Vec3> = 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<Turn> = 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);
+    }
+}

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);
+}

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<f64> {
+/// 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<f64> {
     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<f64> = (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 starts: Vec<f64> = (0..=initial)
+        .map(|i| {
             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();
+        })
+        .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<f64> {
+    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<f64> {
+    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<f64> = (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));
     }