diff --git a/crates/doc/src/eval.rs b/crates/doc/src/eval.rs index 8542351e..fad29d84 100644 --- a/crates/doc/src/eval.rs +++ b/crates/doc/src/eval.rs @@ -62,12 +62,14 @@ impl ModelBody { pub fn definition_name(&self) -> &str { self.source_name.as_deref().unwrap_or(&self.name) } - /// Display mesh (chord tolerance 1/1000 of the body size), computed once. + /// Display mesh, computed once: chord tolerance 1/1000 of the body size, and no facet of a + /// curved face or edge turning more than [`kernel::DISPLAY_ANGLE`], so that a bend or a + /// hole that is small beside the body is still drawn round. pub fn mesh(&self) -> Arc { self.mesh .get_or_init(|| { let tol = (self.body.size() * 1e-3).max(1e-3); - Arc::new(self.body.display_mesh(tol).unwrap_or_default()) + Arc::new(self.body.display_mesh_angular(tol, kernel::DISPLAY_ANGLE).unwrap_or_default()) }) .clone() } diff --git a/crates/engine/src/display_mesh_tests.rs b/crates/engine/src/display_mesh_tests.rs new file mode 100644 index 00000000..d0902eee --- /dev/null +++ b/crates/engine/src/display_mesh_tests.rs @@ -0,0 +1,153 @@ +//! What the application draws of a large sheet with small bends and holes, and that what +//! reads the same mesh (bounds, picks, inside tests) answers as before. + +use std::f64::consts::{FRAC_PI_2, TAU}; + +use serde_json::{Value, json}; +use solidcraft_geom::{Mesh, Vec3}; +use solidcraft_kernel::{DISPLAY_ANGLE, FaceSurface, face_surface}; + +use crate::Session; + +fn run(s: &mut Session, id: &str, p: Value) -> Value { + match s.execute(id, &p) { + Ok(v) => v, + Err(e) => panic!("{id} {p}: {e}"), + } +} + +const HOLES: [[f64; 2]; 4] = [[30.0, 20.0], [1500.0, 20.0], [2970.0, 980.0], [30.0, 980.0]]; + +/// A 3000 x 1000 x 2 sheet with bend radius 1, four holes of diameter 4 and a 40 mm flange on +/// each long side. +fn large_sheet() -> Session { + let mut s = Session::default(); + run(&mut s, "sheet.manage_rules", json!({"name": "AL2", "thickness": "2 mm", "k_factor": "0.5", "bend_radius": "1 mm", "active": true})); + run(&mut s, "sketch.create", json!({"plane": "XY", "name": "Base"})); + run(&mut s, "sketch.rectangle.two_point", json!({"p0": [0, 0], "p1": [3000, 1000]})); + run(&mut s, "sketch.finish", json!({})); + run(&mut s, "sheet.flange", json!({"sketch": "Base"})); + run(&mut s, "sketch.create", json!({"plane": "XY", "name": "Holes"})); + for c in HOLES { + run(&mut s, "sketch.circle.center", json!({"center": c, "diameter": 4})); + } + run(&mut s, "sketch.finish", json!({})); + run(&mut s, "solid.extrude", json!({"sketch": "Holes", "distance": 10, "through_all": true, "operation": "cut"})); + run(&mut s, "sheet.flange", json!({"edges": [[1500, 0, 2], [1500, 1000, 2]], "height": 40})); + s +} + +/// The widest angle between neighbouring lines of vertices around a cylinder face, and how +/// many such steps the face has. +fn steps_around(m: &Mesh, face: usize, origin: Vec3, axis: Vec3) -> (usize, f64) { + let u = axis.any_perp(); + let v = axis.cross(u); + let mut a: Vec = m + .triangles + .iter() + .zip(&m.tri_face) + .filter(|(_, f)| **f as usize == face) + .flat_map(|(t, _)| t.iter()) + .map(|k| m.positions[*k as usize] - origin) + .map(|d| d.dot(v).atan2(d.dot(u))) + .collect(); + a.sort_by(f64::total_cmp); + a.dedup_by(|x, y| (*x - *y).abs() < 1e-7); + // Leave out the widest gap: the part of the turn the face does not cover. + let n = a.len(); + let mut gaps: Vec = (0..n).map(|i| (a[(i + 1) % n] - a[i]).rem_euclid(TAU)).collect(); + gaps.sort_by(f64::total_cmp); + gaps.pop(); + (gaps.len(), gaps.last().copied().unwrap_or(0.0)) +} + +#[test] +fn a_large_sheet_shows_its_small_bends_and_holes_round() { + let s = large_sheet(); + let st = s.world_state(); + let b = &st.bodies[0]; + let m = b.mesh(); + assert!(b.body.size() * 1e-3 > 3.0, "the chord tolerance is more than any radius here"); + let (mut bends, mut walls) = (0, 0); + for face in 0..b.body.face_count() { + let Some(FaceSurface::Cylinder { origin, axis, radius }) = face_surface(&b.body, face) else { continue }; + let (steps, widest) = steps_around(&m, face, origin, axis); + assert!(widest <= DISPLAY_ANGLE * 1.001, "radius {radius}: a facet spans {} degrees", widest.to_degrees()); + if (radius - 2.0).abs() < 1e-6 { + // Half of a hole's wall. + walls += 1; + assert!(steps >= (TAU / 2.0 / DISPLAY_ANGLE).ceil() as usize, "{steps} facets around half a hole"); + } else { + // A bend's inner (1 mm) or outer (3 mm) side: a quarter turn. + bends += 1; + assert!((radius - 1.0).abs() < 1e-6 || (radius - 3.0).abs() < 1e-6, "radius {radius}"); + assert!(steps >= (FRAC_PI_2 / DISPLAY_ANGLE).ceil() as usize, "radius {radius}: {steps} facets across the bend"); + } + } + assert_eq!((bends, walls), (4, 8)); + // The edges the viewport draws: every hole's rim is round, on both sides of the sheet. + for [x, y] in HOLES { + for z in [0.0, 2.0] { + let on = |p: &Vec3| (p.z - z).abs() < 1e-6 && ((p.x - x).hypot(p.y - y) - 2.0).abs() < 1e-6; + let rim: Vec<&Vec> = m.edges.iter().filter(|e| e.iter().all(on)).collect(); + let segments: usize = rim.iter().map(|e| e.len() - 1).sum(); + assert!(segments >= (TAU / DISPLAY_ANGLE).ceil() as usize, "{segments} segments around the rim at {x}, {y}, {z}"); + let widest = rim.iter().flat_map(|e| e.windows(2)).map(|w| 2.0 * (w[0].dist(w[1]) / 4.0).asin()).fold(0.0, f64::max); + assert!(widest <= DISPLAY_ANGLE * 1.001, "a rim segment spans {} degrees", widest.to_degrees()); + } + } + // And the scene the renderers take draws those same polylines. + let scene = crate::view::scene(&s, &crate::view::home_camera(&s)); + let drawn: usize = + scene.lines.iter().filter(|l| l.points.iter().all(|p| ((p.x - 30.0).hypot(p.y - 20.0) - 2.0).abs() < 1e-6)).map(|l| l.points.len() - 1).sum(); + assert!(drawn >= 2 * (TAU / DISPLAY_ANGLE).ceil() as usize, "{drawn} rim segments in the scene"); +} + +#[test] +fn what_reads_the_display_mesh_of_a_large_sheet_answers_as_before() { + let s = large_sheet(); + let st = s.world_state(); + let b = &st.bodies[0]; + let m = b.mesh(); + let coarse = b.body.display_mesh((b.body.size() * 1e-3).max(1e-3)).unwrap(); + // Bounds: the sheet's box is set by flat faces, so the rounder bends move nothing. + let (now, before) = (m.bounds(), coarse.bounds()); + assert!(now.min.dist(before.min) < 1e-9 && now.max.dist(before.max) < 1e-9, "{now:?} {before:?}"); + let exact = solidcraft_kernel::measure(&b.body).unwrap().bbox; + assert!(now.min.dist(exact.min) < 1e-6 && now.max.dist(exact.max) < 1e-6, "{now:?} {exact:?}"); + // The cost: the sheet takes 2212 triangles against the chord mesh's 92 when measured, a factor of + // 24; what grows is the number of round faces and the facets across each, not the size of the body. + assert!( + m.triangles.len() <= 30 * coarse.triangles.len() && m.triangles.len() < 4000, + "{} triangles, {} at the chord tolerance", + m.triangles.len(), + coarse.triangles.len() + ); + // The same faces and edges, in the same order, between the same faces, with the same ends. + assert_eq!(m.edge_faces, coarse.edge_faces); + assert_eq!(m.tri_face.iter().max(), coarse.tri_face.iter().max()); + assert_eq!(m.edges.len(), coarse.edges.len()); + for (fine, rough) in m.edges.iter().zip(&coarse.edges) { + assert!(fine[0].dist(rough[0]) < 1e-9 && fine[fine.len() - 1].dist(rough[rough.len() - 1]) < 1e-9); + } + // Picks: a point on a round face (on the true cylinder, which a flat strip across a bend + // misses by up to 0.9 mm) is found on that face, as is one in the middle of a flat face. + let mut round = 0; + for face in 0..b.body.face_count() { + let Some(FaceSurface::Cylinder { origin, axis, radius }) = face_surface(&b.body, face) else { continue }; + let corners: Vec = + m.triangles.iter().zip(&m.tri_face).filter(|(_, f)| **f as usize == face).flat_map(|(t, _)| t.map(|k| m.positions[k as usize])).collect(); + let middle = corners.iter().fold(Vec3::ZERO, |a, p| a + *p) / corners.len() as f64; + // The middle of the face's corners, pushed out onto the cylinder. + let along = origin + axis * (middle - origin).dot(axis); + let p = along + (middle - along).normalized().unwrap() * radius; + assert_eq!(solidcraft_doc::appearance::face_index_at(b, p), Some(face), "radius {radius} at {p:?}"); + round += 1; + } + assert_eq!(round, 12); + let top = solidcraft_doc::appearance::face_index_at(b, Vec3::new(1500.0, 500.0, 2.0)).unwrap(); + assert!(matches!(face_surface(&b.body, top), Some(FaceSurface::Plane { normal, .. }) if normal.dist(Vec3::Z) < 1e-9)); + // Inside and outside: in the sheet, in a flange's wall, in a hole and above the sheet. + assert!(m.contains(Vec3::new(1500.0, 500.0, 1.0)) && m.contains(Vec3::new(1500.0, now.min.y + 1.0, now.max.z - 10.0))); + assert!(!m.contains(Vec3::new(30.0, 20.0, 1.0)) && !m.contains(Vec3::new(1500.0, 500.0, 3.0))); +} diff --git a/crates/engine/src/lib.rs b/crates/engine/src/lib.rs index 8e6f73ff..a2a7a669 100644 --- a/crates/engine/src/lib.rs +++ b/crates/engine/src/lib.rs @@ -793,6 +793,8 @@ mod construct_tests; #[cfg(test)] mod delete_face_tests; #[cfg(test)] +mod display_mesh_tests; +#[cfg(test)] mod extent_tests; #[cfg(test)] mod hole_face_tests; diff --git a/crates/kernel/src/body.rs b/crates/kernel/src/body.rs index 38861a51..f9bb71fa 100644 --- a/crates/kernel/src/body.rs +++ b/crates/kernel/src/body.rs @@ -148,17 +148,38 @@ pub struct FaceInfo { type MeshedShell = truck_topology::Shell, Option>; +/// The most one facet of the display mesh turns (radians): a curved face or edge is divided +/// until no facet spans more of it than this, however small its radius is beside the body. +/// 20 degrees: a quarter turn is drawn with at least five facets, a full circle with at least 18. +pub const DISPLAY_ANGLE: f64 = std::f64::consts::PI / 9.0; + +/// The sag, as a fraction of the chord tolerance, below which the angle no longer divides a +/// facet: at [`DISPLAY_ANGLE`] and the application's tolerance of 1/1000 of the body, a radius +/// under about 1/15000 of the body gets fewer facets than the angle asks for, and one under +/// a millionth of it only those of the chord tolerance. +const DISPLAY_ANGLE_FLOOR: f64 = 1e-3; + +/// The finest angle a display mesh is divided by (one degree). +const MIN_DISPLAY_ANGLE: f64 = std::f64::consts::PI / 180.0; + /// Triangulate each shell (not the solid: truck's solid constructor panics on the open or /// non-manifold shells an import can produce). Boundaries that lie on their surface only within /// tolerance (imports, boolean results) are projected onto it. fn mesh_shells(solid: &Solid, tol: f64) -> Vec { + mesh_shells_within(solid, Deflection::chord(tol)) +} + +/// [`mesh_shells`] divided by a chord tolerance and, when it has one, a bound on the angle a +/// facet spans. Each edge is divided once for both its faces, so the faces still share its +/// points. +fn mesh_shells_within(solid: &Solid, deflection: Deflection) -> Vec { solid .boundaries() .iter() // The robust search tries the exact one first, so exact geometry costs nothing extra; a // face can also come back as an empty mesh rather than none (boolean results with // approximate intersection curves), so always use it. - .map(|sh| sh.robust_triangulation(tol)) + .map(|sh| sh.robust_triangulation_within(deflection)) .collect() } @@ -838,13 +859,24 @@ impl Body { /// Like [`Body::tessellate`]; with `refine`, triangles that cut through curved faces are /// split until within `tol` of the surface (for measuring). pub fn tessellate_with(&self, tol: f64, refine: bool) -> Result { + self.tessellate_within(tol, f64::INFINITY, refine) + } + + /// [`Body::tessellate_with`] in which no facet turns more than `angle` radians either + /// (infinite: the chord tolerance alone). + fn tessellate_within(&self, tol: f64, angle: f64, refine: bool) -> Result { let tol = if tol.is_finite() && tol > 0.0 { tol } else { 0.05 }; let tol = tol.max(self.size() * 1e-6); if let Some(m) = &self.mesh { return Ok(m.to_mesh()); } + let deflection = if angle.is_finite() && angle > 0.0 { + Deflection { chord: tol, angle: angle.max(MIN_DISPLAY_ANGLE), floor: tol * DISPLAY_ANGLE_FLOOR } + } else { + Deflection::chord(tol) + }; guard("tessellate", || { - let meshed = mesh_shells(&self.solid, tol); + let meshed = mesh_shells_within(&self.solid, deflection); let mut out = Mesh::default(); let solid_faces: Vec<&mt::Face> = self.solid.face_iter().collect(); for (fi, face) in meshed.iter().flat_map(|s| s.face_iter()).enumerate() { @@ -941,7 +973,20 @@ impl Body { /// Display mesh: like [`Body::tessellate`], with seam edges flagged. pub fn display_mesh(&self, tol: f64) -> Result { - let mut m = self.tessellate(tol)?; + self.flagged(self.tessellate(tol)?) + } + + /// Display mesh whose curved faces and edges look curved whatever the body's size: like + /// [`Body::display_mesh`], and no facet turns more than `angle` radians (see + /// [`DISPLAY_ANGLE`]). A chord tolerance that follows the body's size is larger than the + /// sagitta of a small radius on a large body (a sheet-metal bend, a small hole), which the + /// tolerance alone leaves as one flat facet. + pub fn display_mesh_angular(&self, tol: f64, angle: f64) -> Result { + self.flagged(self.tessellate_within(tol, angle, false)?) + } + + /// A tessellation with its edges flagged for display (seams, tangent joins, protected creases). + fn flagged(&self, mut m: Mesh) -> Result { if self.is_mesh() { return Ok(m); } @@ -1221,6 +1266,10 @@ mod lump_tests; #[path = "refine_cache_tests.rs"] mod refine_cache_tests; +#[cfg(test)] +#[path = "display_mesh_tests.rs"] +mod display_mesh_tests; + #[cfg(test)] #[path = "regularization_protection_tests.rs"] mod regularization_protection_tests; diff --git a/crates/kernel/src/display_mesh_tests.rs b/crates/kernel/src/display_mesh_tests.rs new file mode 100644 index 00000000..86e15341 --- /dev/null +++ b/crates/kernel/src/display_mesh_tests.rs @@ -0,0 +1,350 @@ +//! The display mesh of a large body with small radii: bends and holes look round, the mesh +//! stays closed, and the cost stays bounded. + +use std::collections::HashMap; +use std::f64::consts::{FRAC_PI_2, PI, TAU}; + +use solidcraft_geom::{Loop2, Mesh, Plane, Region2, Seg2, Vec2, Vec3}; + +use super::{Body, DISPLAY_ANGLE}; +use crate::{FaceSurface, extrude, face_surface}; + +/// The chord tolerance the application draws a body with. +fn display_tol(b: &Body) -> f64 { + (b.size() * 1e-3).max(1e-3) +} + +fn rect(w: f64, h: f64) -> Loop2 { + Loop2::polygon(&[Vec2::new(0.0, 0.0), Vec2::new(w, 0.0), Vec2::new(w, h), Vec2::new(0.0, h)]) +} + +/// A `length` long angle of 2 mm sheet: a 100 mm leg along X, a 90 degree bend of inner +/// radius 1 (outer 3) about the line x = y = 3, and a 40 mm leg along Y, extruded along Z. +fn bent_sheet(length: f64) -> Body { + let line = |a: (f64, f64), b: (f64, f64)| Seg2::Line { a: Vec2::new(a.0, a.1), b: Vec2::new(b.0, b.1) }; + let c = Vec2::new(3.0, 3.0); + let outline = Loop2 { + segs: vec![ + line((100.0, 0.0), (100.0, 2.0)), + line((100.0, 2.0), (3.0, 2.0)), + Seg2::Arc { center: c, radius: 1.0, start: 1.5 * PI, sweep: -FRAC_PI_2 }, + line((2.0, 3.0), (2.0, 40.0)), + line((2.0, 40.0), (0.0, 40.0)), + line((0.0, 40.0), (0.0, 3.0)), + Seg2::Arc { center: c, radius: 3.0, start: PI, sweep: FRAC_PI_2 }, + line((3.0, 0.0), (100.0, 0.0)), + ], + }; + assert!(outline.signed_area() > 0.0); + extrude(&Plane::XY, &[Region2 { outer: outline, holes: vec![] }], 0.0, length).unwrap().pop().unwrap() +} + +/// A 3000 x 1000 x 2 plate with round holes of `radius` at `centres`. +fn plate_with_holes(centres: &[Vec2], radius: f64) -> Body { + let holes = centres.iter().map(|c| Loop2::circle(*c, radius).reversed()).collect(); + extrude(&Plane::XY, &[Region2 { outer: rect(3000.0, 1000.0), holes }], 0.0, 2.0).unwrap().pop().unwrap() +} + +/// The faces of a body on cylinders of this radius. +fn cylinder_faces(b: &Body, radius: f64) -> Vec<(usize, Vec3, Vec3)> { + (0..b.face_count()) + .filter_map(|i| match face_surface(b, i)? { + FaceSurface::Cylinder { origin, axis, radius: r } if (r - radius).abs() < 1e-6 => Some((i, origin, axis)), + _ => None, + }) + .collect() +} + +/// Distance of a point from an axis. +fn radial(p: Vec3, origin: Vec3, axis: Vec3) -> f64 { + let d = p - origin; + (d - axis * d.dot(axis)).len() +} + +/// How a face on a cylinder is divided around its axis: the angles (radians, sorted, from the +/// first) at which it has vertices. +fn angles_around(m: &Mesh, face: usize, origin: Vec3, axis: Vec3) -> Vec { + let u = axis.any_perp(); + let v = axis.cross(u); + let mut a: Vec = m + .triangles + .iter() + .zip(&m.tri_face) + .filter(|(_, f)| **f as usize == face) + .flat_map(|(t, _)| t.iter()) + .filter_map(|k| m.positions.get(*k as usize)) + .map(|p| (*p - origin).dot(v).atan2((*p - origin).dot(u))) + .collect(); + a.sort_by(f64::total_cmp); + // Start at the widest gap, so an arc that crosses ±π stays in one piece. + let n = a.len(); + let gap = |i: usize| (a[(i + 1) % n] - a[i]).rem_euclid(TAU); + let start = (0..n).max_by(|i, j| gap(*i).total_cmp(&gap(*j))).map_or(0, |i| (i + 1) % n); + let first = a[start]; + let mut out: Vec = (0..n).map(|i| (a[(start + i) % n] - first).rem_euclid(TAU)).collect(); + out.dedup_by(|x, y| (*x - *y).abs() < 1e-7); + out +} + +/// The most a triangle of `face` lies inside the cylinder it is on (its corners are on it). +fn sag(m: &Mesh, face: usize, origin: Vec3, axis: Vec3, radius: f64) -> f64 { + let mut worst = 0.0f64; + for (t, _) in m.triangles.iter().zip(&m.tri_face).filter(|(_, f)| **f as usize == face) { + let [a, b, c] = m.tri(t).unwrap(); + for p in [a, b, c] { + assert!((radial(p, origin, axis) - radius).abs() < 1e-6, "a corner off the cylinder"); + } + for p in [(a + b) * 0.5, (b + c) * 0.5, (c + a) * 0.5, (a + b + c) / 3.0] { + worst = worst.max((radius - radial(p, origin, axis)).abs()); + } + } + worst +} + +/// Triangle sides that are not shared by exactly two triangles once corners closer than `eps` +/// are joined (a closed mesh has none). +fn open_sides(m: &Mesh, eps: f64) -> usize { + let cell = |p: &Vec3| ((p.x / eps).floor() as i64, (p.y / eps).floor() as i64, (p.z / eps).floor() as i64); + let mut grid: HashMap<(i64, i64, i64), Vec> = HashMap::new(); + let mut points: Vec = Vec::new(); + let ids: Vec = m + .positions + .iter() + .map(|p| { + let (cx, cy, cz) = cell(p); + for (dx, dy, dz) in (-1..=1).flat_map(|x| (-1..=1).flat_map(move |y| (-1..=1).map(move |z| (x, y, z)))) { + if let Some(i) = grid.get(&(cx + dx, cy + dy, cz + dz)).and_then(|v| v.iter().find(|i| points[**i].dist(*p) < eps)) { + return *i; + } + } + points.push(*p); + grid.entry((cx, cy, cz)).or_default().push(points.len() - 1); + points.len() - 1 + }) + .collect(); + let mut sides: HashMap<(usize, usize), i32> = HashMap::new(); + for t in &m.triangles { + let [a, b, c] = t.map(|k| ids[k as usize]); + // A triangle at a pole has two corners in one point: no area, no sides. + if a == b || b == c || c == a { + continue; + } + for (p, q) in [(a, b), (b, c), (c, a)] { + *sides.entry((p.min(q), p.max(q))).or_default() += 1; + } + } + sides.values().filter(|n| **n != 2).count() +} + +/// The polyline segments around a circle of `radius` about `centre` in the plane z = `z`, +/// and the widest angle one of them spans. +fn rim(m: &Mesh, centre: Vec2, radius: f64, z: f64) -> (usize, f64) { + let on = |p: &Vec3| (p.z - z).abs() < 1e-6 && (Vec2::new(p.x, p.y).dist(centre) - radius).abs() < 1e-6; + let (mut segments, mut widest) = (0, 0.0f64); + for e in m.edges.iter().filter(|e| e.len() >= 2 && e.iter().all(on)) { + for w in e.windows(2) { + segments += 1; + widest = widest.max(2.0 * (w[0].dist(w[1]) / (2.0 * radius)).clamp(0.0, 1.0).asin()); + } + } + (segments, widest) +} + +#[test] +fn a_small_bend_on_a_large_sheet_is_meshed_round() { + let b = bent_sheet(3000.0); + let tol = display_tol(&b); + // The tolerance is larger than either side's sagitta: 1 - cos 45° of the radius. + assert!(tol > 3.0 * (1.0 - (PI / 4.0).cos())); + let coarse = b.display_mesh(tol).unwrap(); + let round = b.display_mesh_angular(tol, DISPLAY_ANGLE).unwrap(); + // A quarter turn in facets of at most DISPLAY_ANGLE. + let facets = (FRAC_PI_2 / DISPLAY_ANGLE).ceil() as usize; + assert_eq!(facets, 5); + for radius in [1.0, 3.0] { + let faces = cylinder_faces(&b, radius); + assert_eq!(faces.len(), 1, "one bend face of radius {radius}"); + let (face, origin, axis) = faces[0]; + // The chord tolerance alone: one flat strip from tangent line to tangent line. + assert_eq!(angles_around(&coarse, face, origin, axis).len(), 2, "radius {radius}"); + assert!(sag(&coarse, face, origin, axis, radius) > 0.29 * radius); + let a = angles_around(&round, face, origin, axis); + assert!(a.len() > facets, "radius {radius}: {} facets across the bend", a.len() - 1); + assert!((a[a.len() - 1] - FRAC_PI_2).abs() < 1e-6, "the bend is a quarter turn: {a:?}"); + let widest = a.windows(2).map(|w| w[1] - w[0]).fold(0.0, f64::max); + assert!(widest <= DISPLAY_ANGLE * 1.001, "radius {radius}: a facet spans {} degrees", widest.to_degrees()); + // Within the chord tolerance of the cylinder, and by much: a facet that turns no more + // than the angle sags no more than 1 - cos(angle / 2) of the radius. + let s = sag(&round, face, origin, axis, radius); + assert!(s <= radius * (1.0 - (DISPLAY_ANGLE / 2.0).cos()) * 1.001 && s < tol, "radius {radius}: sag {s}"); + } + // The flat faces along the sheet meet the bends in straight lines and need nothing more: + // only the bends, and the two ends where their arcs now have more points, gain triangles. + let count = |m: &Mesh, face: usize| m.tri_face.iter().filter(|f| **f as usize == face).count(); + let mut along = 0; + for face in 0..b.face_count() { + if let Some(FaceSurface::Plane { normal, .. }) = face_surface(&b, face) + && normal.dot(Vec3::Z).abs() < 0.5 + { + along += 1; + assert_eq!(count(&round, face), count(&coarse, face), "face {face}"); + } + } + assert_eq!(along, 6); +} + +#[test] +fn the_round_display_mesh_stays_closed_and_keeps_its_edges_and_bounds() { + let b = bent_sheet(3000.0); + let tol = display_tol(&b); + let coarse = b.display_mesh(tol).unwrap(); + let round = b.display_mesh_angular(tol, DISPLAY_ANGLE).unwrap(); + // Neighbouring faces share the points of the edge between them. + assert_eq!(open_sides(&round, 1e-6), 0); + // The same edges in the same order, between the same faces, with the same ends and flags. + assert_eq!(round.edges.len(), coarse.edges.len()); + assert_eq!(round.edge_faces, coarse.edge_faces); + assert_eq!((&round.seams, &round.smooth_edges, &round.sharp_edges), (&coarse.seams, &coarse.smooth_edges, &coarse.sharp_edges)); + for (fine, rough) in round.edges.iter().zip(&coarse.edges) { + assert!(fine.len() >= rough.len()); + assert!(fine[0].dist(rough[0]) < 1e-9 && fine[fine.len() - 1].dist(rough[rough.len() - 1]) < 1e-9); + } + // The body's box: flat faces and tangent lines bound it, so nothing moves. + let (r, c) = (round.bounds(), coarse.bounds()); + assert!(r.min.dist(c.min) < 1e-9 && r.max.dist(c.max) < 1e-9); + // Closer to the solid: the flat strips cut the outer bend's corner off and fill the inner one. + let exact = crate::measure(&b).unwrap().volume; + assert!((round.measure().volume - exact).abs() < (coarse.measure().volume - exact).abs()); +} + +#[test] +fn a_small_hole_on_a_large_plate_has_a_round_rim() { + let centre = Vec2::new(30.0, 20.0); + let b = plate_with_holes(&[centre], 2.0); + let tol = display_tol(&b); + assert!(tol > 2.0, "the hole's radius is below the chord tolerance"); + let coarse = b.display_mesh(tol).unwrap(); + let round = b.display_mesh_angular(tol, DISPLAY_ANGLE).unwrap(); + let segments = (TAU / DISPLAY_ANGLE).ceil() as usize; + assert_eq!(segments, 18); + for z in [0.0, 2.0] { + // The chord tolerance alone: each half circle is one segment, the rim a line. + assert_eq!(rim(&coarse, centre, 2.0, z).0, 2); + let (n, widest) = rim(&round, centre, 2.0, z); + assert!(n >= segments, "{n} segments around the rim at z = {z}"); + assert!(widest <= DISPLAY_ANGLE * 1.001, "a segment spans {} degrees", widest.to_degrees()); + } + // The hole's wall follows its rims. + for (face, origin, axis) in cylinder_faces(&b, 2.0) { + let a = angles_around(&round, face, origin, axis); + assert!(a.windows(2).map(|w| w[1] - w[0]).fold(0.0, f64::max) <= DISPLAY_ANGLE * 1.001); + assert!(sag(&round, face, origin, axis, 2.0) <= 2.0 * (1.0 - (DISPLAY_ANGLE / 2.0).cos()) * 1.001); + } + assert_eq!(open_sides(&round, 1e-6), 0); +} + +#[test] +fn an_arc_the_chord_tolerance_already_divides_finely_is_left_alone() { + // A 40 mm cylinder at a tolerance of 0.01 mm: 64 facets around, far inside the angle. + let b = crate::cylinder(Vec3::ZERO, Vec3::Z, 20.0, 30.0).unwrap(); + let coarse = b.display_mesh(0.01).unwrap(); + let round = b.display_mesh_angular(0.01, DISPLAY_ANGLE).unwrap(); + assert_eq!(round.positions, coarse.positions); + assert_eq!(round.triangles, coarse.triangles); + assert_eq!(round.edges, coarse.edges); + // And a body with nothing curved at any tolerance. + let b = crate::box_solid(Vec3::ZERO, Vec3::new(3000.0, 1000.0, 2.0)).unwrap(); + let tol = display_tol(&b); + assert_eq!(b.display_mesh_angular(tol, DISPLAY_ANGLE).unwrap().positions, b.display_mesh(tol).unwrap().positions); +} + +#[test] +fn a_radius_far_below_the_tolerance_keeps_its_chord_facets() { + // A half circle of radius 0.002 in a 3 m plate sags less than a thousandth of the chord + // tolerance: below that floor nothing is divided for the angle, however far the body + // outgrows its smallest radius. + let centre = Vec2::new(30.0, 20.0); + let b = plate_with_holes(&[centre], 0.002); + let tol = display_tol(&b); + let coarse = b.display_mesh(tol).unwrap(); + let round = b.display_mesh_angular(tol, DISPLAY_ANGLE).unwrap(); + assert_eq!(round.triangles.len(), coarse.triangles.len()); + assert_eq!(round.edges, coarse.edges); + // Just above the floor the rim is divided, but coarsely: the angle stops where the sag + // falls below the floor, not where the facets reach the angle. + let b = plate_with_holes(&[centre], 0.01); + let round = b.display_mesh_angular(display_tol(&b), DISPLAY_ANGLE).unwrap(); + let (n, _) = rim(&round, centre, 0.01, 2.0); + assert!((3..18).contains(&n), "{n} segments"); +} + +#[test] +fn a_face_crowded_with_holes_gets_a_coarser_angle_and_bounded_work() { + // 1000 holes of diameter 4 in a 3 m plate. At the full angle the top and bottom faces + // would each be bounded by more than 20 000 segments, and the work of triangulating a + // face grows faster than its boundary. + let centres: Vec = (0..50).flat_map(|i| (0..20).map(move |j| Vec2::new(30.0 + 60.0 * i as f64, 25.0 + 50.0 * j as f64))).collect(); + let b = plate_with_holes(¢res, 2.0); + let tol = display_tol(&b); + let coarse = b.display_mesh(tol).unwrap(); + let round = b.display_mesh_angular(tol, DISPLAY_ANGLE).unwrap(); + let boundary = |m: &Mesh, z: f64| m.edges.iter().filter(|e| e.iter().all(|p| (p.z - z).abs() < 1e-6)).map(|e| e.len() - 1).sum::(); + let budget = truck_meshalgo::tessellation::Deflection::MAX_ANGLE_BOUNDARY; + for z in [0.0, 2.0] { + assert_eq!(boundary(&coarse, z), 4 + 2 * centres.len()); + let n = boundary(&round, z); + assert!(n > boundary(&coarse, z), "the rims are still rounder than the chord tolerance leaves them: {n}"); + assert!(n <= budget, "{n} segments bound the face at z = {z}"); + } + // Every rim got the same share, and the walls follow their rims. + let per_rim: Vec = centres.iter().map(|c| rim(&round, *c, 2.0, 2.0).0).collect(); + assert!(per_rim.iter().all(|n| *n == per_rim[0] && *n >= 3), "{:?}", &per_rim[..4]); + assert!( + round.triangles.len() < 4 * coarse.triangles.len(), + "{} triangles, {} at the chord tolerance", + round.triangles.len(), + coarse.triangles.len() + ); + assert_eq!(open_sides(&round, 1e-6), 0); + // A plate with few holes keeps the full angle. + let few = plate_with_holes(¢res[..40], 2.0); + let round = few.display_mesh_angular(display_tol(&few), DISPLAY_ANGLE).unwrap(); + assert!(rim(&round, centres[0], 2.0, 2.0).0 >= 18); +} + +#[test] +fn round_faces_of_every_kind_are_divided_by_the_angle() { + // A sphere, a torus and a cone a few millimetres across, drawn at the 3 mm tolerance of a + // body 3 m long: more than any of their sagittas. + let disc = Region2 { outer: Loop2::circle(Vec2::ZERO, 4.0), holes: vec![] }; + let bodies = [ + ("sphere", crate::sphere(Vec3::ZERO, 4.0).unwrap()), + ("torus", crate::torus(Vec3::ZERO, 6.0, 2.0).unwrap()), + ("cone", crate::extrude_tapered(&Plane::XY, &disc, 6.0, 1.0, -0.3).unwrap()), + ]; + // The most the surface normal turns along one side of a triangle of a curved face. + let turn = |b: &Body, m: &Mesh| { + let flat: Vec = (0..b.face_count()).map(|i| matches!(face_surface(b, i), Some(FaceSurface::Plane { .. }))).collect(); + m.triangles + .iter() + .zip(&m.tri_face) + .filter(|(_, f)| !flat[**f as usize]) + .flat_map(|(t, _)| [(t[0], t[1]), (t[1], t[2]), (t[2], t[0])]) + .filter(|(a, b)| m.positions[*a as usize].dist(m.positions[*b as usize]) > 1e-9) + .map(|(a, b)| m.normals[a as usize].dot(m.normals[b as usize]).clamp(-1.0, 1.0).acos()) + .fold(0.0, f64::max) + }; + let check = |name: &str, b: &Body| { + let coarse = b.display_mesh(3.0).unwrap(); + let round = b.display_mesh_angular(3.0, DISPLAY_ANGLE).unwrap(); + assert!(turn(b, &coarse) > 2.0 * DISPLAY_ANGLE, "{name}: the chord tolerance leaves facets of {} degrees", turn(b, &coarse).to_degrees()); + // A side crosses at most one facet in each direction. + assert!(turn(b, &round) <= 1.5 * DISPLAY_ANGLE, "{name}: a side turns {} degrees", turn(b, &round).to_degrees()); + assert_eq!(open_sides(&round, 1e-6), 0, "{name}"); + // Bounded: no more than a full turn of facets each way on every face. + assert!(round.triangles.len() <= b.face_count() * 2 * 36 * 36, "{name}: {} triangles", round.triangles.len()); + }; + for (name, b) in &bodies { + check(name, b); + } +} diff --git a/crates/kernel/src/lib.rs b/crates/kernel/src/lib.rs index 367bc6f3..f7b0f74e 100644 --- a/crates/kernel/src/lib.rs +++ b/crates/kernel/src/lib.rs @@ -73,7 +73,7 @@ mod tspline; pub use blend::{ChamferSide, chamfer, chamfer_sides, fillet, fillet_max}; pub use blend_corner::fillet_radii; pub use blend_variable::{fillet_chord, fillet_variable, fillet_variable_points}; -pub use body::{Body, CapSignature, EdgeInfo, FaceInfo, FacePaint, FaceSource, Paint}; +pub use body::{Body, CapSignature, DISPLAY_ANGLE, EdgeInfo, FaceInfo, FacePaint, FaceSource, Paint}; pub use boolean_diagnostics::{BooleanDiagnostics, BooleanFallback, BooleanStrategy, capture_boolean_diagnostics}; pub use build::{PathSeg, box_solid, canonical_revolve_axis, cylinder, extrude, extrude_tapered, loft, revolve, sphere, sweep, torus}; pub use delete_face::delete_faces; diff --git a/vendor/README.md b/vendor/README.md index 9839a2e2..c3c2c70d 100644 --- a/vendor/README.md +++ b/vendor/README.md @@ -62,6 +62,16 @@ Changes: point that does not evaluate counts as flat instead of recursing). Surface grids refine only the direction that bends, and a ruling left undivided gets cells comparable to the other direction (long triangles from a trim loop to the grid's far ends cut through cylinders). +- `tessellation/mod.rs`, `tessellation/triangulation.rs`, for display: a `Deflection` is the + chord tolerance plus an optional bound on the angle one facet spans, and + `RoundMeshableShape::robust_triangulation_within` divides by it: an edge or a grid direction + within the chord tolerance is still halved while a piece turns more than the angle, until its + sag is under a floor (a radius that is small beside the tolerance was one flat facet). + Bounded work: a face that would be bounded by more than 4096 segments has its curved edges + divided at twice the angle, and again, down to the chord tolerance alone, and the angle + stops refining a face's grid at 16 384 points. Each edge is still divided once for both its + faces. `triangulation` and `robust_triangulation` divide by the chord tolerance alone, as + before. - `lib.rs`: compiler warnings no longer fail the build. Example, test and bench targets are dropped (their sources are not vendored). diff --git a/vendor/truck-meshalgo/src/tessellation/mod.rs b/vendor/truck-meshalgo/src/tessellation/mod.rs index 2e70d1a0..6a32259a 100644 --- a/vendor/truck-meshalgo/src/tessellation/mod.rs +++ b/vendor/truck-meshalgo/src/tessellation/mod.rs @@ -44,6 +44,71 @@ impl> RobustMesh type PolylineCurve = truck_polymesh::PolylineCurve; +/// SolidCraft: how finely curves and surfaces are divided into facets. `chord` is the usual +/// tolerance: the most a facet may sag from the true shape. `angle`, when finite, also bounds +/// what one facet spans: a facet that turns more than `angle` is split whatever its size (a +/// radius that is small beside the tolerance would otherwise stay one flat facet), until its +/// sag is below `floor`. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct Deflection { + /// Most a facet may sag from the curve or surface. + pub chord: f64, + /// Most a facet may turn (radians: the arc it spans on a circle); infinite for no bound. + pub angle: f64, + /// Sag below which `angle` no longer splits a facet. + pub floor: f64, +} + +impl Deflection { + /// Most boundary segments a face may get from the angle bound. Triangulating a face tests + /// each triangle against the boundary loops around it, so the work grows faster than the + /// boundary (with its square for one long loop): the curved edges of a face that would + /// get more are divided at twice the angle, and again, down to the chord tolerance alone. + pub const MAX_ANGLE_BOUNDARY: usize = 4096; + /// Most grid points a face may get from the angle bound (beyond them only `chord` refines). + pub const MAX_ANGLE_GRID: usize = 16_384; + /// How often the angle is doubled for a crowded face before only `chord` is left. + pub const MAX_RELAX: u8 = 4; + + /// The chord tolerance alone: what `triangulation(tol)` divides by. + pub const fn chord(tol: f64) -> Self { + Self { + chord: tol, + angle: f64::INFINITY, + floor: 0.0, + } + } + + /// Whether `angle` bounds anything. + pub fn angular(&self) -> bool { self.angle.is_finite() && self.angle > 0.0 } + + /// The same with the angle doubled `level` times; from [`Deflection::MAX_RELAX`] times, or + /// past half a turn, the chord tolerance alone. + pub fn relaxed(&self, level: u8) -> Self { + let angle = self.angle * f64::from(1u32 << level.min(Self::MAX_RELAX)); + match level == 0 || (level < Self::MAX_RELAX && angle < std::f64::consts::PI) { + true => Self { angle, ..*self }, + false => Self::chord(self.chord), + } + } + + /// The squared sag over squared length above which a facet turns more than `angle`: a + /// chord spanning `a` of a circle sags `tan(a / 4) / 2` of its length. Infinite without + /// an angle bound. + pub(super) fn turn2(&self) -> f64 { + match self.angular() && self.angle < std::f64::consts::PI { + true => (f64::tan(self.angle / 4.0) / 2.0).powi(2), + false => f64::INFINITY, + } + } + + /// Whether the angle bound splits a facet of squared sag `sag2` and squared length `len2` + /// (`turn2` from [`Deflection::turn2`]). + pub(super) fn turns(&self, turn2: f64, sag2: f64, len2: f64) -> bool { + turn2.is_finite() && sag2 > self.floor * self.floor && sag2 > turn2 * len2 + } +} + /// Trait for converting tessellated shape into polygon. pub trait MeshedShape { /// Converts tessellated shape into polygon. @@ -258,11 +323,15 @@ impl MeshableShape for Shell Self::MeshedShape { nonpositive_tolerance!(tol); #[cfg(not(target_arch = "wasm32"))] - let res = triangulation::shell_tessellation(self, tol, triangulation::by_search_parameter); + let res = triangulation::shell_tessellation( + self, + Deflection::chord(tol), + triangulation::by_search_parameter, + ); #[cfg(target_arch = "wasm32")] let res = triangulation::shell_tessellation_single_thread( self, - tol, + Deflection::chord(tol), triangulation::by_search_parameter, ); res @@ -272,17 +341,39 @@ impl MeshableShape for Shell RobustMeshableShape for Shell { type MeshedShape = Shell>; fn robust_triangulation(&self, tol: f64) -> Self::MeshedShape { - nonpositive_tolerance!(tol); + self.robust_triangulation_within(Deflection::chord(tol)) + } +} + +/// SolidCraft: tessellation with a bound on the angle one facet spans as well as on its sag. +pub trait RoundMeshableShape { + /// Shape whose edges are made polylines and faces polygon surface. + type MeshedShape: MeshedShape; + /// [`RobustMeshableShape::robust_triangulation`] divided by `deflection`: with + /// [`Deflection::chord`] the very same shape; with an angle, curved edges and faces whose + /// radius is small beside the chord tolerance are divided as well. Every edge is still + /// divided once, for both faces it bounds, so neighbouring faces keep sharing its points. + /// + /// # Panics + /// + /// `deflection.chord` must be more than `TOLERANCE`. + fn robust_triangulation_within(&self, deflection: Deflection) -> Self::MeshedShape; +} + +impl RoundMeshableShape for Shell { + type MeshedShape = Shell>; + fn robust_triangulation_within(&self, deflection: Deflection) -> Self::MeshedShape { + nonpositive_tolerance!(deflection.chord); #[cfg(not(target_arch = "wasm32"))] let res = triangulation::shell_tessellation( self, - tol, + deflection, triangulation::by_search_nearest_parameter, ); #[cfg(target_arch = "wasm32")] let res = triangulation::shell_tessellation_single_thread( self, - tol, + deflection, triangulation::by_search_nearest_parameter, ); res diff --git a/vendor/truck-meshalgo/src/tessellation/triangulation.rs b/vendor/truck-meshalgo/src/tessellation/triangulation.rs index 8b7d0533..8523718d 100644 --- a/vendor/truck-meshalgo/src/tessellation/triangulation.rs +++ b/vendor/truck-meshalgo/src/tessellation/triangulation.rs @@ -91,7 +91,7 @@ where #[cfg(not(target_arch = "wasm32"))] pub(super) fn shell_tessellation<'a, C, S>( shell: &Shell, - tol: f64, + tol: Deflection, sp: impl SP, ) -> MeshedShell where @@ -104,18 +104,28 @@ where .map(|v| (v.id(), v.mapped(Point3::clone))) .collect(); let eset: HashMap<_, _> = shell.edge_par_iter().map(move |e| (e.id(), e)).collect(); + // SolidCraft: the polylines first, so that the edges of a face crowded by the angle bound + // can be divided again before any face is meshed. + let mut polys: HashMap<_, _> = eset + .par_iter() + .map(|(id, edge)| { + let curve = edge.curve(); + (*id, (bounded_polyline(&curve, curve.range_tuple(), tol), 0)) + }) + .collect(); + relax_crowded_faces(shell, &mut polys, tol); let edge_map: HashMap<_, _> = eset .into_par_iter() .map(move |(id, edge)| { let v0 = vmap.get(&edge.absolute_front().id()).unwrap(); let v1 = vmap.get(&edge.absolute_back().id()).unwrap(); - let curve = edge.curve(); - let poly = bounded_polyline(&curve, curve.range_tuple(), tol); - (id, Edge::new_unchecked(v0, v1, poly)) + let (poly, level) = polys.get(&id).unwrap().clone(); + (id, (Edge::new_unchecked(v0, v1, poly), level)) }) .collect(); + let edge_map = &edge_map; let create_edge = |edge: &Edge| -> Edge<_, _> { - let new_edge = edge_map.get(&edge.id()).unwrap(); + let (new_edge, _) = edge_map.get(&edge.id()).unwrap(); match edge.orientation() { true => new_edge.clone(), false => new_edge.inverse(), @@ -129,16 +139,83 @@ where .iter() .map(create_boundary) .collect(); - shell_create_polygon(&face.surface(), wires, face.orientation(), tol, &sp, cap) + // SolidCraft: a face is divided like the coarsest of its edges. + let level = face + .absolute_boundaries() + .iter() + .flat_map(Wire::edge_iter) + .filter_map(|edge| edge_map.get(&edge.id()).map(|(_, level)| *level)) + .max() + .unwrap_or(0); + shell_create_polygon( + &face.surface(), + wires, + face.orientation(), + tol.relaxed(level), + &sp, + cap, + ) }; shell.face_par_iter().map(create_face).collect() } +/// SolidCraft: bounded work for the angle bound. A face bounded by more than +/// `Deflection::MAX_ANGLE_BOUNDARY` segments has its curved edges divided again at twice the +/// angle, as often as it takes, down to the chord tolerance alone (a plate with a thousand +/// small holes is not worth a hundred times the work to draw each hole round). `polys` holds +/// each edge's polyline and how often its angle was doubled. Every edge still has one +/// polyline for both its faces, and the outcome does not depend on the order of the faces. +fn relax_crowded_faces( + shell: &Shell, + polys: &mut HashMap, (PolylineCurve, u8)>, + tol: Deflection, +) where + C: PolylineableCurve, +{ + if !tol.angular() { + return; + } + let faces: Vec>> = shell + .face_iter() + .map(|face| { + let wires = face.absolute_boundaries(); + wires.iter().flat_map(Wire::edge_iter).cloned().collect() + }) + .collect(); + for level in 1..=Deflection::MAX_RELAX { + let segments = |edge: &Edge| match polys.get(&edge.id()) { + // Straight at the full angle: nothing to relax, and no say in how its faces are divided. + Some((poly, 0)) if poly.len() <= 2 => (1, false), + Some((poly, done)) => (poly.len().saturating_sub(1), *done < level), + None => (0, false), + }; + let crowded: Vec<&Edge> = faces + .iter() + .filter(|edges| { + edges.iter().map(|edge| segments(edge).0).sum::() + > Deflection::MAX_ANGLE_BOUNDARY + }) + .flatten() + .filter(|edge| segments(edge).1) + .collect(); + if crowded.is_empty() { + break; + } + let relaxed = tol.relaxed(level); + for edge in crowded { + if let Some(slot) = polys.get_mut(&edge.id()).filter(|(_, done)| *done < level) { + let curve = edge.curve(); + *slot = (bounded_polyline(&curve, curve.range_tuple(), relaxed), level); + } + } + } +} + /// Tessellates faces #[cfg(any(target_arch = "wasm32", test))] pub(super) fn shell_tessellation_single_thread<'a, C, S>( shell: &'a Shell, - tol: f64, + tol: Deflection, sp: impl SP, ) -> MeshedShell where @@ -148,6 +225,16 @@ where use truck_base::entry_map::FxEntryMap as EntryMap; use truck_topology::Vertex as TVertex; let cap = face_cap(shell.len()); + // SolidCraft: the polylines first (see `shell_tessellation`). + let mut polys: HashMap<_, _> = shell + .edge_iter() + .map(|edge| { + let curve = edge.curve(); + (edge.id(), (bounded_polyline(&curve, curve.range_tuple(), tol), 0)) + }) + .collect(); + relax_crowded_faces(shell, &mut polys, tol); + let polys = &polys; let mut vmap = EntryMap::new( move |v: &TVertex| v.id(), move |v| v.mapped(Point3::clone), @@ -159,8 +246,7 @@ where let v0 = vmap.entry_or_insert(vf).clone(); let vb = edge.absolute_back(); let v1 = vmap.entry_or_insert(vb).clone(); - let curve = edge.curve(); - let poly = bounded_polyline(&curve, curve.range_tuple(), tol); + let (poly, _) = polys.get(&edge.id()).unwrap().clone(); Edge::new_unchecked(&v0, &v1, poly) }, ); @@ -180,7 +266,21 @@ where .iter() .map(&mut create_boundary) .collect(); - shell_create_polygon(&face.surface(), wires, face.orientation(), tol, &sp, cap) + let level = face + .absolute_boundaries() + .iter() + .flat_map(Wire::edge_iter) + .filter_map(|edge| polys.get(&edge.id()).map(|(_, level)| *level)) + .max() + .unwrap_or(0); + shell_create_polygon( + &face.surface(), + wires, + face.orientation(), + tol.relaxed(level), + &sp, + cap, + ) }; shell.face_iter().map(create_face).collect() } @@ -197,6 +297,7 @@ where { let vertices = shell.vertices.clone(); let cap = face_cap(shell.faces.len()); + let tol = Deflection::chord(tol); let tessellate_edge = |edge: &CompressedEdge| { let curve = &edge.curve; CompressedEdge { @@ -245,7 +346,7 @@ fn shell_create_polygon( surface: &S, wires: Vec>, orientation: bool, - tol: f64, + tol: Deflection, sp: impl SP, cap: usize, ) -> Face> { @@ -382,7 +483,11 @@ fn loop_orientation(curve: &[Point2]) -> bool { } impl PolyBoundary { - fn new(pieces: Vec, surface: &impl PreMeshableSurface, tol: f64) -> Self { + fn new( + pieces: Vec, + surface: &impl PreMeshableSurface, + tol: Deflection, + ) -> Self { let (mut closed, mut open) = (Vec::new(), Vec::new()); pieces.into_iter().for_each(|PolyBoundaryPiece(mut vec)| { match vec[0].distance(vec[vec.len() - 1]) < 1.0e-3 { @@ -577,8 +682,15 @@ fn spade_round(x: f64) -> f64 { } /// Tessellates one surface trimmed by polyline. -fn trimming_tessellation(surface: &S, polyboundary: &PolyBoundary, tol: f64, cap: usize) -> PolygonMesh -where S: PreMeshableSurface { +fn trimming_tessellation( + surface: &S, + polyboundary: &PolyBoundary, + tol: Deflection, + cap: usize, +) -> PolygonMesh +where + S: PreMeshableSurface, +{ let mut triangulation = Cdt::new(); polyboundary.insert_to(&mut triangulation); insert_surface(&mut triangulation, surface, polyboundary, tol, cap); @@ -597,7 +709,7 @@ fn insert_surface( triangulation: &mut Cdt, surface: impl PreMeshableSurface, polyline: &PolyBoundary, - tol: f64, + tol: Deflection, cap: usize, ) { let bdb: BoundingBox = polyline.0.iter().flatten().collect(); @@ -656,8 +768,10 @@ pub(super) fn face_cap(faces: usize) -> usize { /// SolidCraft: truck's curve division, bounded: at most `MAX_EDGE_POINTS` points and a /// bisection depth of 24; a midpoint that does not evaluate (NaN) counts as flat (truck's -/// version recursed on it 100 levels deep, which never finishes). -fn bounded_polyline(curve: &C, range: (f64, f64), tol: f64) -> PolylineCurve { +/// version recursed on it 100 levels deep, which never finishes). A piece within the chord +/// tolerance is still halved while it turns more than the deflection's angle. +fn bounded_polyline(curve: &C, range: (f64, f64), d: Deflection) -> PolylineCurve { + let (tol, turn2) = (d.chord, d.turn2()); let mut out = vec![curve.subs(range.0)]; let mut stack: Vec<(f64, f64, usize)> = vec![(range.0, range.1, 0)]; while let Some((a, b, depth)) = stack.pop() { @@ -668,6 +782,7 @@ fn bounded_polyline(curve: &C, range: (f64, f64), tol: f64 let tm = 0.5 * (a + b); let dm = curve.subs(tm).distance2(pa.midpoint(pb)); let flat = !(d2 >= tol * tol || dm >= tol * tol); + let flat = flat && !d.turns(turn2, d2.max(dm), pa.distance2(pb)); if flat || depth >= 24 || out.len() + stack.len() >= MAX_EDGE_POINTS { out.push(pb); } else { @@ -680,12 +795,17 @@ fn bounded_polyline(curve: &C, range: (f64, f64), tol: f64 } /// SolidCraft: truck's surface division, bounded: stops refining when the grid would exceed -/// `cap` points or after 24 rounds; cells that do not evaluate (NaN) are not refined. -fn bounded_surface_division(surface: &S, (ur, vr): ((f64, f64), (f64, f64)), tol: f64, cap: usize) -> (Vec, Vec) { +/// `cap` points or after 24 rounds; cells that do not evaluate (NaN) are not refined. A +/// direction within the chord tolerance is still halved where a cell turns more than the +/// deflection's angle along it, while the grid stays within `Deflection::MAX_ANGLE_GRID`. +fn bounded_surface_division(surface: &S, (ur, vr): ((f64, f64), (f64, f64)), d: Deflection, cap: usize) -> (Vec, Vec) { + let (tol, turn2) = (d.chord, d.turn2()); let (mut udiv, mut vdiv) = (vec![ur.0, ur.1], vec![vr.0, vr.1]); for _ in 0..24 { let mut fu = vec![false; udiv.len() - 1]; let mut fv = vec![false; vdiv.len() - 1]; + // What the angle alone asks for. + let (mut au, mut av) = (fu.clone(), fv.clone()); for (i, u) in udiv.windows(2).enumerate() { for (j, v) in vdiv.windows(2).enumerate() { if fu[i] && fv[j] { @@ -702,10 +822,12 @@ fn bounded_surface_division(surface: &S, (ur, vr): ((f64, let far = p0.distance2(pt) > tol * tol; // Which direction bends: refine only that one (a cylinder needs no division // along its rulings). - let alu = surface.subs(u[0], v0).to_vec() * (1.0 - p) + surface.subs(u[1], v0).to_vec() * p; - let alv = surface.subs(u0, v[0]).to_vec() * (1.0 - q) + surface.subs(u0, v[1]).to_vec() * q; - let bend_u = p0.distance2(Point3::from_vec(alu)) > tol * tol; - let bend_v = p0.distance2(Point3::from_vec(alv)) > tol * tol; + let (su, sv) = ([surface.subs(u[0], v0), surface.subs(u[1], v0)], [surface.subs(u0, v[0]), surface.subs(u0, v[1])]); + let alu = su[0].to_vec() * (1.0 - p) + su[1].to_vec() * p; + let alv = sv[0].to_vec() * (1.0 - q) + sv[1].to_vec() * q; + let (du, dv) = (p0.distance2(Point3::from_vec(alu)), p0.distance2(Point3::from_vec(alv))); + let bend_u = du > tol * tol; + let bend_v = dv > tol * tol; if bend_u || bend_v { fu[i] |= bend_u; fv[j] |= bend_v; @@ -713,8 +835,20 @@ fn bounded_surface_division(surface: &S, (ur, vr): ((f64, fu[i] = true; fv[j] = true; } + au[i] |= d.turns(turn2, du, su[0].distance2(su[1])); + av[j] |= d.turns(turn2, dv, sv[0].distance2(sv[1])); } } + let count = |f: &[bool], a: &[bool]| f.iter().zip(a).filter(|(f, a)| **f || **a).count(); + let (nu, nv) = (udiv.len() + count(&fu, &au), vdiv.len() + count(&fv, &av)); + // The angle refines only while the grid stays small: past that the chord tolerance alone. + let turning = turn2.is_finite() && nu.saturating_mul(nv) <= cap.min(Deflection::MAX_ANGLE_GRID); + if !turning { + au.fill(false); + av.fill(false); + } + fu.iter_mut().zip(&au).for_each(|(f, a)| *f |= *a); + fv.iter_mut().zip(&av).for_each(|(f, a)| *f |= *a); let nu = udiv.len() + fu.iter().filter(|x| **x).count(); let nv = vdiv.len() + fv.iter().filter(|x| **x).count(); if nu == udiv.len() && nv == vdiv.len() { @@ -848,15 +982,42 @@ fn polyline_on_surface( surface: impl PreMeshableSurface, p: Point2, q: Point2, - tol: f64, + tol: Deflection, ) -> Vec { use truck_geometry::prelude::*; let line = Line(p, q); let pcurve = PCurve::new(line, surface); - let (vec, _) = pcurve.parameter_division(pcurve.range_tuple(), tol); + let (vec, _) = pcurve.parameter_division(pcurve.range_tuple(), tol.chord); + let vec = halved_while_turning(|t| pcurve.subs(t), vec, tol); vec.into_iter().map(|t| line.subs(t)).collect() } +/// SolidCraft: `params` with every interval halved while the curve turns more over it than the +/// deflection's angle (at most 12 times each, and `MAX_EDGE_POINTS` in all). Unchanged without +/// an angle. +fn halved_while_turning(curve: impl Fn(f64) -> Point3, params: Vec, d: Deflection) -> Vec { + let turn2 = d.turn2(); + if !turn2.is_finite() || params.len() < 2 { + return params; + } + let mut out = vec![params[0]]; + for w in params.windows(2) { + let mut stack: Vec<(f64, f64, usize)> = vec![(w[0], w[1], 0)]; + while let Some((a, b, depth)) = stack.pop() { + let (pa, pb, tm) = (curve(a), curve(b), 0.5 * (a + b)); + let sag2 = curve(tm).distance2(pa.midpoint(pb)); + let split = d.turns(turn2, sag2, pa.distance2(pb)); + if split && depth < 12 && out.len() + stack.len() < MAX_EDGE_POINTS { + stack.push((tm, b, depth + 1)); + stack.push((a, tm, depth + 1)); + } else { + out.push(b); + } + } + } + out +} + #[test] #[ignore] #[cfg(not(target_arch = "wasm32"))] @@ -872,13 +1033,13 @@ fn par_bench() { let instant = Instant::now(); (0..100).for_each(|_| { - let _shell = shell_tessellation(&shell, 0.01, by_search_parameter); + let _shell = shell_tessellation(&shell, Deflection::chord(0.01), by_search_parameter); }); println!("{}ms", instant.elapsed().as_millis()); let instant = Instant::now(); (0..100).for_each(|_| { - let _shell = shell_tessellation_single_thread(&shell, 0.01, by_search_parameter); + let _shell = shell_tessellation_single_thread(&shell, Deflection::chord(0.01), by_search_parameter); }); println!("{}ms", instant.elapsed().as_millis()); }