diff --git a/crates/kernel/src/body.rs b/crates/kernel/src/body.rs index 775841b0..d38f03a9 100644 --- a/crates/kernel/src/body.rs +++ b/crates/kernel/src/body.rs @@ -1043,6 +1043,17 @@ impl Body { if let Some(m) = &self.mesh { return Ok(m.to_mesh()); } + // A face trimmed by a helix (a modelled thread from another CAD system) turns more than + // once around its cylinder; truck meshes it with triangles across the solid (#1289). + // Such a body is meshed from its display plan, which cuts the face along a seam. + if guard("tessellate", || Ok(crate::tess::seam::turns_more_than_once(&self.solid))).unwrap_or(false) { + let turn = (angle.is_finite() && angle > 0.0) + .then(|| crate::tess::size::Turn { angle: angle.max(MIN_DISPLAY_ANGLE), floor: tol * DISPLAY_ANGLE_FLOOR }); + match guard("tessellate", || crate::tess::plan(&self.solid, tol, turn)) { + Ok((plan, _)) => return Ok(plan.evaluate()), + Err(e) => log::debug!("helical faces: no display plan ({e}); meshing the reference way"), + } + } let deflection = if angle.is_finite() && angle > 0.0 { Deflection { chord: tol, angle: angle.max(MIN_DISPLAY_ANGLE), floor: tol * DISPLAY_ANGLE_FLOOR, rulings: true } } else { diff --git a/crates/kernel/src/lib.rs b/crates/kernel/src/lib.rs index 9a3c4e29..89bc7bbe 100644 --- a/crates/kernel/src/lib.rs +++ b/crates/kernel/src/lib.rs @@ -103,6 +103,8 @@ mod sweep_path; mod tangent_cylinders; mod tess; #[cfg(test)] +mod thread_display_tests; +#[cfg(test)] mod thread_tests; mod topo; mod tspline; diff --git a/crates/kernel/src/tess/cdt.rs b/crates/kernel/src/tess/cdt.rs index ec0df6d4..17a547fe 100644 --- a/crates/kernel/src/tess/cdt.rs +++ b/crates/kernel/src/tess/cdt.rs @@ -102,7 +102,14 @@ pub(crate) fn triangulate(boundary: &[[f64; 2]], loops: &[std::ops::Range .exclude_outer_faces(true) .with_angle_limit(AngleLimit::from_deg(MIN_ANGLE_DEG)) .with_max_additional_vertices(budget.saturating_sub(added)); + let before = cdt.num_vertices(); cdt.refine(params); + let shaped = cdt.num_vertices().saturating_sub(before); + // Shaping flips sides and adds points of its own; the triangles it leaves may + // sag again (a long strip round a cylinder), so the sag is checked once more. + if shaped > 0 { + refine_sag(&mut cdt, sag, budget.saturating_sub(added + shaped), &from); + } } // A split constraint would put a point on the boundary the neighbour does not have. if cdt.num_constraints() != constraints { diff --git a/crates/kernel/src/tess/mod.rs b/crates/kernel/src/tess/mod.rs index a5421b54..8c1cfc73 100644 --- a/crates/kernel/src/tess/mod.rs +++ b/crates/kernel/src/tess/mod.rs @@ -24,6 +24,7 @@ pub(crate) mod edges; pub(crate) mod grid; pub(crate) mod motorcycle; pub(crate) mod pieces; +pub(crate) mod seam; pub(crate) mod size; pub(crate) mod ssi; pub(crate) mod surface; @@ -237,10 +238,8 @@ fn plan_face( let dom = uv::domain(s); let accept = (tol * 10.0).max(1e-6); let fail = |why: &str| unsupported(format!("face {fi}: {why}")); - // Project every loop, unwrapping periodic parameters along it and splitting poles. - let mut fixed: Vec = Vec::new(); - let mut uvs: Vec<[f64; 2]> = Vec::new(); - let mut ranges: Vec> = Vec::new(); + // Project every loop, unwrapping periodic parameters along it. + let mut projected: Vec = Vec::with_capacity(loops3.len()); for (li, pts) in loops3.iter().enumerate() { let mut luv: Vec<[f64; 2]> = Vec::with_capacity(pts.len()); let mut prev: Option<[f64; 2]> = None; @@ -260,16 +259,28 @@ fn plan_face( luv.push(q); prev = Some(q); } - // A loop that winds once around a periodic direction does not close in (u, v). - if let (Some(first), Some(last)) = (luv.first(), luv.last()) { - for d in 0..2 { - if let Some(per) = dom.period[d] - && (uv::unwrap(first[d], last[d], Some(per)) - first[d]).abs() > per * 0.25 - { - return Err(fail("boundary winds around the surface")); - } - } + projected.push(seam::Loop { uv: luv, pts: pts.clone() }); + } + // A loop that winds around a periodic direction does not close in (u, v): the face is cut + // along a parameter line into one period (a seamless cylinder face, a helical trim). + let wound: Vec = (0..2).filter(|&d| dom.period[d].is_some_and(|per| seam::winds(&projected, d, per))).collect(); + // A face cut along a seam is a long winding strip whose sides are the cut lines and the + // boundary in turn: Delaunay with sag refinement meshes it (the patch layouts below assume + // boundary all round, and leave fans across the strip's corners off the surface). + let seamed = !wound.is_empty(); + match wound.as_slice() { + [] => {} + [d] => { + let per = dom.period[*d].unwrap_or(std::f64::consts::TAU); + projected = seam::cut(s, &projected, *d, per, rules, h).ok_or_else(|| fail("boundary winds around the surface"))?; } + _ => return Err(fail("boundary winds around the surface both ways")), + } + // Split poles. + let mut fixed: Vec = Vec::new(); + let mut uvs: Vec<[f64; 2]> = Vec::new(); + let mut ranges: Vec> = Vec::new(); + for seam::Loop { uv: luv, pts } in &projected { let start = uvs.len(); let n = luv.len(); for (k, (q, p)) in luv.iter().zip(pts).enumerate() { @@ -325,7 +336,7 @@ fn plan_face( 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)) { + let structured = match (ranges.len(), if seamed { None } else { grid::quad(&uvs) }) { (1, Some(q)) => { // Interior lines: the boundary's own where opposite sides agree (refined where the // surface turns more between them inside the face), else the division. @@ -350,7 +361,7 @@ fn plan_face( _ => None, }; let structured = match structured { - None if ranges.len() == 1 => { + None if ranges.len() == 1 && !seamed => { 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, h, Some(&fine)) } @@ -371,7 +382,9 @@ fn plan_face( // else Delaunay. A plane: the butterfly of a round face, else Delaunay (a plane // is shaded the same by any triangles, and its Delaunay needs no refinement). let face = Region { s, uvs: &uvs, fixed: &fixed, ranges: &ranges, rules, h, speed, us: &us, vs: &vs }; - let trimmed = if planar { + let trimmed = if seamed { + None + } else if planar { if ranges.len() == 1 { plan_pieces(&face) } else { None } } else { plan_pieces(&face).or_else(|| { diff --git a/crates/kernel/src/tess/seam.rs b/crates/kernel/src/tess/seam.rs new file mode 100644 index 00000000..21cb6ef1 --- /dev/null +++ b/crates/kernel/src/tess/seam.rs @@ -0,0 +1,405 @@ +//! Faces on a closed surface whose boundary winds around it: a cylinder face bounded by two +//! circles and no seam edge, or one trimmed by a helix that turns several times (a modelled +//! thread, as other CAD systems write them to STEP). In the surface's unwrapped parameters such +//! a face repeats every period, and its loops do not close. It is cut along a parameter line +//! (the seam it never had) into one period: every loop, and every copy of it a whole number of +//! periods along, is clipped to the strip between the line and its next copy, and the parts of +//! the two lines inside the face join the pieces into closed loops. Both lines are one curve +//! on the surface, so their nodes are made once and shared: the face stays closed across it. + +use solidcraft_geom::Vec3; +use solidcraft_geom::tess::EvalSurface; + +use super::size::Rules; + +/// One loop of a face: parameters (unwrapped along the loop) and the 3D node each stands for. +#[derive(Clone, Debug, Default)] +pub(crate) struct Loop { + pub uv: Vec<[f64; 2]>, + pub pts: Vec, +} + +/// Most loop copies one face may be clipped into (a helix of this many turns). +const MAX_COPIES: i64 = 2048; +/// Most nodes a cut line may get between two crossings. +const MAX_SEAM_NODES: usize = 4096; + +/// How far `l` drifts along `d` going once around it (a whole number of periods for a closed +/// loop: zero unless it winds around the surface). +pub(crate) fn drift(l: &Loop, d: usize, per: f64) -> f64 { + let (Some(first), Some(last)) = (l.uv.first(), l.uv.last()) else { return 0.0 }; + let close = super::uv::unwrap(first[d], last[d], Some(per)) - last[d]; + last[d] - first[d] + close +} + +/// Whether any loop winds around the closed direction `d`. +pub(crate) fn winds(loops: &[Loop], d: usize, per: f64) -> bool { + loops.iter().any(|l| drift(l, d, per).abs() > per * 0.5) +} + +/// A point where the boundary crosses a cut line: its place along the line, its 3D point, and +/// whether the face lies above it (towards larger parameters across the line) or below. +#[derive(Clone, Copy, Debug)] +struct Crossing { + y: f64, + p: Vec3, + above: bool, +} + +/// A boundary run across the strip: it enters through one cut line and leaves through one. +struct Piece { + uv: Vec<[f64; 2]>, + pts: Vec, +} + +/// Cut a face whose loops wind around direction `d` (period `per`) into one period. The face +/// lies left of its loops in (u, v). `None` when the loops cannot be cut cleanly (the face would +/// be unbounded along the line, or the two lines disagree); the caller then fails as before. +pub(crate) fn cut(s: &EvalSurface, loops: &[Loop], d: usize, per: f64, rules: &Rules, h: f64) -> Option> { + if !(per > 0.0 && per.is_finite()) || loops.is_empty() { + return None; + } + // Work in (x, y) with x along the closed direction; swapping (u, v) mirrors the plane, so the + // loops are reversed to keep the face on their left. + let swap = d == 1; + let local: Vec = loops + .iter() + .map(|l| { + let mut uv: Vec<[f64; 2]> = l.uv.iter().map(|q| if swap { [q[1], q[0]] } else { *q }).collect(); + let mut pts = l.pts.clone(); + if swap { + uv.reverse(); + pts.reverse(); + } + Loop { uv, pts } + }) + .collect(); + let x0 = cut_position(&local, per)?; + let x1 = x0 + per; + let mut closed: Vec = Vec::new(); + let mut pieces: Vec = Vec::new(); + for l in &local { + let n = l.uv.len(); + if n < 2 || l.pts.len() != n { + return None; + } + let w = drift(l, 0, per); + let (lo, hi) = l.uv.iter().fold((f64::MAX, f64::MIN), |(a, b), q| (a.min(q[0]), b.max(q[0]))); + if !(lo.is_finite() && hi.is_finite()) { + return None; + } + if w.abs() <= per * 0.5 { + // A closed loop: each copy along x that reaches into the strip. + let k0 = ((x0 - hi) / per).floor() as i64; + let k1 = ((x1 - lo) / per).ceil() as i64; + if k1 - k0 > MAX_COPIES { + return None; + } + for k in k0..=k1 { + let shift = per * k as f64; + let uv: Vec<[f64; 2]> = l.uv.iter().map(|q| [q[0] + shift, q[1]]).collect(); + if uv.iter().all(|q| q[0] > x0 && q[0] < x1) { + closed.push(Loop { uv, pts: l.pts.clone() }); + } else if uv.iter().any(|q| q[0] > x0 && q[0] < x1) || crosses(&uv, x0, x1) { + // Start at a node outside the strip, so every run inside is whole. + let start = uv.iter().position(|q| !(q[0] > x0 && q[0] < x1))?; + let idx: Vec = (0..=n).map(|i| (start + i) % n).collect(); + let path: Vec<[f64; 2]> = idx.iter().filter_map(|&i| uv.get(i).copied()).collect(); + let pts: Vec = idx.iter().filter_map(|&i| l.pts.get(i).copied()).collect(); + clip(&path, &pts, x0, x1, &mut pieces); + } + } + } else { + // A loop winding around: the endless path it repeats along, over enough turns to + // start and end outside the strip. + let turns = (w / per).round(); + if turns.abs() != 1.0 { + return None; + } + let reach = (hi - lo) / per; + let span = (reach.ceil() as i64) + 2; + if span > MAX_COPIES { + return None; + } + // Copies j with lo + j·w ≤ x1 and hi + j·w ≥ x0, plus one each side. + let (ja, jb) = if w > 0.0 { + (((x0 - hi) / w).floor() as i64 - 1, ((x1 - lo) / w).ceil() as i64 + 1) + } else { + (((x1 - lo) / w).floor() as i64 - 1, ((x0 - hi) / w).ceil() as i64 + 1) + }; + if jb - ja > MAX_COPIES { + return None; + } + let mut path = Vec::new(); + let mut pts = Vec::new(); + for j in ja..=jb { + let shift = w * j as f64; + path.extend(l.uv.iter().map(|q| [q[0] + shift, q[1]])); + pts.extend(l.pts.iter().copied()); + } + let (Some(a), Some(b)) = (path.first(), path.last()) else { return None }; + if (a[0] > x0 && a[0] < x1) || (b[0] > x0 && b[0] < x1) { + return None; + } + clip(&path, &pts, x0, x1, &mut pieces); + } + } + // The crossings on each line, sorted along it; both lines must see the same ones. + let mut left: Vec<(Crossing, usize, bool)> = Vec::new(); + let mut right: Vec<(Crossing, usize, bool)> = Vec::new(); + for (i, p) in pieces.iter().enumerate() { + let (Some(a), Some(b), Some(pa), Some(pb)) = (p.uv.first(), p.uv.last(), p.pts.first(), p.pts.last()) else { return None }; + // Entering through the left line (moving +x): the face is above; through the right + // line (moving −x): below. Leaving through the left (−x): below; the right (+x): above. + let enter = Crossing { y: a[1], p: *pa, above: a[0] == x0 }; + let leave = Crossing { y: b[1], p: *pb, above: b[0] == x1 }; + if a[0] == x0 { + left.push((enter, i, true)) + } else { + right.push((enter, i, true)) + } + if b[0] == x0 { left.push((leave, i, false)) } else { right.push((leave, i, false)) } + } + left.sort_by(|a, b| a.0.y.total_cmp(&b.0.y)); + right.sort_by(|a, b| a.0.y.total_cmp(&b.0.y)); + if left.len() != right.len() || !left.len().is_multiple_of(2) { + return None; + } + let scale = loops.iter().flat_map(|l| l.uv.iter()).fold(1.0f64, |m, q| m.max(q[0].abs()).max(q[1].abs())); + for (a, b) in left.iter().zip(&right) { + if a.0.above != b.0.above || (a.0.y - b.0.y).abs() > 1e-7 * scale { + return None; + } + } + // The same crossing on both lines is one point. + for (a, b) in left.clone().iter().zip(right.iter_mut()) { + b.0.y = a.0.y; + b.0.p = a.0.p; + } + for (c, i, start) in left.iter().chain(&right) { + let p = pieces.get_mut(*i)?; + let k = if *start { 0 } else { p.uv.len().checked_sub(1)? }; + if let (Some(q), Some(pt)) = (p.uv.get_mut(k), p.pts.get_mut(k)) { + q[1] = c.y; + *pt = c.p; + } + } + // Inside the face along each line: from a crossing with the face above to the next one + // up (which must have it below). + let mut links: Vec<(usize, Vec<[f64; 2]>, Vec, usize)> = Vec::new(); + for pair in left.chunks(2).zip(right.chunks(2)) { + let (lp, rp) = pair; + let (Some(la), Some(lb), Some(ra), Some(rb)) = (lp.first(), lp.get(1), rp.first(), rp.get(1)) else { return None }; + if !(la.0.above && !lb.0.above) { + return None; + } + // Nodes strictly between, made once for both lines. + let inner = seam_nodes(s, swap, x0, la.0.y, lb.0.y, rules, h)?; + // Left line, downwards: from the piece leaving at the top to the one entering at the bottom. + let (top_leave, bottom_enter) = (lb.1, la.1); + let down_uv: Vec<[f64; 2]> = inner.iter().rev().map(|(y, _)| [x0, *y]).collect(); + let down_pts: Vec = inner.iter().rev().map(|(_, p)| *p).collect(); + if lb.2 || !la.2 { + return None; + } + links.push((top_leave, down_uv, down_pts, bottom_enter)); + // Right line, upwards: from the piece leaving at the bottom to the one entering at the top. + if ra.2 || !rb.2 { + return None; + } + let up_uv: Vec<[f64; 2]> = inner.iter().map(|(y, _)| [x1, *y]).collect(); + let up_pts: Vec = inner.iter().map(|(_, p)| *p).collect(); + links.push((ra.1, up_uv, up_pts, rb.1)); + } + // Chain pieces through the links into closed loops. + let mut next: Vec> = vec![None; pieces.len()]; + let mut via: Vec> = vec![None; pieces.len()]; + for (li, (from, _, _, to)) in links.iter().enumerate() { + if next.get(*from).copied().flatten().is_some() { + return None; + } + *next.get_mut(*from)? = Some(*to); + *via.get_mut(*from)? = Some(li); + } + let mut used = vec![false; pieces.len()]; + for start in 0..pieces.len() { + if used.get(start).copied().unwrap_or(true) { + continue; + } + let mut out = Loop::default(); + let mut at = start; + for _ in 0..=pieces.len() { + if *used.get(at)? { + break; + } + *used.get_mut(at)? = true; + let p = pieces.get(at)?; + out.uv.extend(&p.uv); + out.pts.extend(&p.pts); + let (_, luv, lpts, _) = links.get((*via.get(at)?)?)?; + out.uv.extend(luv); + out.pts.extend(lpts); + at = (*next.get(at)?)?; + } + if at != start { + return None; + } + closed.push(out); + } + // Back to (u, v). + let out: Vec = closed + .into_iter() + .map(|mut l| { + if swap { + l.uv = l.uv.iter().rev().map(|q| [q[1], q[0]]).collect(); + l.pts.reverse(); + } + l + }) + .collect(); + // The face must come out with positive area (left of its loops). + let area: f64 = out + .iter() + .map(|l| { + let n = l.uv.len(); + (0..n) + .map(|i| { + let (a, b) = (l.uv[i], l.uv[(i + 1) % n]); + a[0] * b[1] - b[0] * a[1] + }) + .sum::() + }) + .sum(); + (area > 0.0).then_some(out) +} + +/// Where to cut: the middle of the widest gap between the loops' nodes around the period, as +/// far from any node as can be. +fn cut_position(loops: &[Loop], per: f64) -> Option { + let mut xs: Vec = loops.iter().flat_map(|l| l.uv.iter().map(|q| q[0].rem_euclid(per))).collect(); + if xs.is_empty() { + return None; + } + xs.sort_by(f64::total_cmp); + let mut best = (0.0, f64::NEG_INFINITY); + for (i, x) in xs.iter().enumerate() { + let next = xs.get(i + 1).copied().unwrap_or(xs[0] + per); + if next - x > best.1 { + best = (0.5 * (x + next), next - x); + } + } + (best.1 > 0.0).then_some(best.0) +} + +/// Whether a path jumps over the strip in one segment. +fn crosses(uv: &[[f64; 2]], x0: f64, x1: f64) -> bool { + uv.windows(2).any(|w| (w[0][0] <= x0 && w[1][0] >= x1) || (w[1][0] <= x0 && w[0][0] >= x1)) +} + +/// The runs of an open path (starting and ending outside the strip) inside the strip, each +/// from the point where it enters to the one where it leaves. The crossing points lie on the +/// boundary's own straight 3D segments, so the neighbouring face's edge passes through them. +fn clip(path: &[[f64; 2]], pts: &[Vec3], x0: f64, x1: f64, out: &mut Vec) { + let inside = |q: &[f64; 2]| q[0] > x0 && q[0] < x1; + let mut cur: Option = None; + for k in 0..path.len().saturating_sub(1) { + let (Some(a), Some(b), Some(pa), Some(pb)) = (path.get(k), path.get(k + 1), pts.get(k), pts.get(k + 1)) else { break }; + let at = |x: f64| { + let t = ((x - a[0]) / (b[0] - a[0])).clamp(0.0, 1.0); + ([x, a[1] + (b[1] - a[1]) * t], *pa + (*pb - *pa) * t) + }; + // Lines crossed by this segment, in order along it. + let mut hits: Vec = [x0, x1].into_iter().filter(|x| (a[0] - x) * (b[0] - x) < 0.0).collect(); + if b[0] < a[0] { + hits.reverse(); + } + for x in hits { + let (q, p) = at(x); + match cur.take() { + Some(mut piece) => { + piece.uv.push(q); + piece.pts.push(p); + out.push(piece); + } + None => cur = Some(Piece { uv: vec![q], pts: vec![p] }), + } + } + if inside(b) + && let Some(piece) = cur.as_mut() + { + piece.uv.push(*b); + piece.pts.push(*pb); + } + } +} + +/// Nodes along the cut line between two crossings (`y0` < `y1`, ends excluded): within the +/// tolerance, the side length and the facet angle of the surface along the line. +fn seam_nodes(s: &EvalSurface, swap: bool, x: f64, y0: f64, y1: f64, rules: &Rules, h: f64) -> Option> { + let uv = |y: f64| if swap { [y, x] } else { [x, y] }; + let mut ys = vec![y0, y1]; + let mut i = 0; + while i + 1 < ys.len() { + let (a, b) = (ys[i], ys[i + 1]); + if ys.len() < MAX_SEAM_NODES && !super::uv::fine_between(s, rules, h, uv(a), uv(b)) && b - a > 1e-12 * (1.0 + a.abs()) { + ys.insert(i + 1, 0.5 * (a + b)); + } else { + i += 1; + } + } + let inner = ys.get(1..ys.len() - 1)?; + Some(inner.iter().map(|&y| (y, s.eval(uv(y)[0], uv(y)[1]).0)).collect()) +} + +/// Whether a face of `solid` lies on a surface of revolution and has a boundary loop turning +/// more than once around its axis (a helical trim). Truck's mesher cannot mesh such a face (its +/// triangles cut across the solid), so the reference mesh is made from the display plan instead. +pub(crate) fn turns_more_than_once(solid: &crate::body::Solid) -> bool { + use truck_meshalgo::prelude::*; + use truck_modeling as mt; + let tau = std::f64::consts::TAU; + for face in solid.face_iter() { + let mt::Surface::RevolutedCurve(proc) = face.surface() else { continue }; + let rev = proc.entity(); + let m = *proc.transform(); + let o = crate::body::from_p3(m.transform_point(rev.origin())); + let a = m.transform_vector(rev.axis()); + let Some(k) = Vec3::new(a.x, a.y, a.z).normalized() else { continue }; + let e1 = k.cross(if k.x.abs() < 0.9 { Vec3::X } else { Vec3::Y }).normalized().unwrap_or(Vec3::X); + let e2 = k.cross(e1); + let angle = |p: mt::Point3| { + let d = crate::body::from_p3(p) - o; + d.dot(e2).atan2(d.dot(e1)) + }; + for wire in face.boundaries() { + let (mut total, mut lo, mut hi) = (0.0f64, 0.0f64, 0.0f64); + let mut prev: Option = None; + for e in wire.edge_iter() { + let c = e.oriented_curve(); + let (t0, t1) = c.range_tuple(); + // Samples dense enough that no step turns half way round. + let mut n = 16usize; + let steps = loop { + let xs: Vec = (0..=n).map(|i| angle(c.subs(t0 + (t1 - t0) * i as f64 / n as f64))).collect(); + let ok = xs.windows(2).all(|w| ((w[1] - w[0] + std::f64::consts::PI).rem_euclid(tau) - std::f64::consts::PI).abs() < 1.0); + if ok || n >= 4096 { + break xs; + } + n *= 4; + }; + for x in steps { + if let Some(p) = prev { + total += (x - p + std::f64::consts::PI).rem_euclid(tau) - std::f64::consts::PI; + lo = lo.min(total); + hi = hi.max(total); + } + prev = Some(x); + } + } + if hi - lo > tau * 1.01 { + return true; + } + } + } + false +} diff --git a/crates/kernel/src/thread_display_tests.rs b/crates/kernel/src/thread_display_tests.rs new file mode 100644 index 00000000..2b696b2f --- /dev/null +++ b/crates/kernel/src/thread_display_tests.rs @@ -0,0 +1,258 @@ +//! Issue 1289: the display mesh of a modelled screw thread (helical B-spline flanks). + +use solidcraft_geom::tess::TessPlan; +use solidcraft_geom::{Mesh, Vec3}; + +use std::f64::consts::{PI, TAU}; +use std::sync::Arc; + +use crate::{Body, DISPLAY_ANGLE, EdgeSpec, FaceSpec, SurfSpec, sew}; + +/// A helix about Z of radius `r` from angle `a0` to `a1`, rising `pitch` per turn from `z0` at +/// angle 0, as cubic Bézier pieces of at most `TAU / 16`. +fn helix(r: f64, z0: f64, pitch: f64, a0: f64, a1: f64) -> EdgeSpec { + let n = (((a1 - a0).abs() / (TAU / 16.0)).ceil() as usize).max(1); + let h = (a1 - a0) / n as f64; + let k = 4.0 / 3.0 * (h / 4.0).tan(); + let at = |a: f64| Vec3::new(r * a.cos(), r * a.sin(), z0 + pitch * a / TAU); + let rise = Vec3::new(0.0, 0.0, pitch * h / TAU / 3.0); + let turn = |a: f64| Vec3::new(-a.sin(), a.cos(), 0.0) * (r * k); + let mut ctrl = vec![at(a0)]; + for i in 0..n { + let (s, e) = (a0 + h * i as f64, a0 + h * (i + 1) as f64); + ctrl.extend([at(s) + turn(s) + rise, at(e) - turn(e) - rise, at(e)]); + } + EdgeSpec::Chain { ctrl: ctrl.into() } +} + +fn rev(e: &EdgeSpec) -> EdgeSpec { + match e { + EdgeSpec::Chain { ctrl } => EdgeSpec::Chain { ctrl: ctrl.iter().rev().copied().collect::>() }, + EdgeSpec::Line { a, b } => EdgeSpec::Line { a: *b, b: *a }, + EdgeSpec::Arc { a, b, mid } => EdgeSpec::Arc { a: *b, b: *a, mid: *mid }, + } +} + +fn line(a: Vec3, b: Vec3) -> EdgeSpec { + EdgeSpec::Line { a, b } +} + +fn ends(e: &EdgeSpec) -> (Vec3, Vec3) { + match e { + EdgeSpec::Chain { ctrl } => (ctrl[0], ctrl[ctrl.len() - 1]), + EdgeSpec::Line { a, b } | EdgeSpec::Arc { a, b, .. } => (*a, *b), + } +} + +/// A full circle about Z at height `z`, as two arcs. +fn circle(r: f64, z: f64) -> Vec { + let p = |a: f64| Vec3::new(r * a.cos(), r * a.sin(), z); + vec![EdgeSpec::Arc { a: p(0.3), b: p(0.3 + PI), mid: p(0.3 + PI / 2.0) }, EdgeSpec::Arc { a: p(0.3 + PI), b: p(0.3), mid: p(0.3 + 1.5 * PI) }] +} + +/// A threaded shaft the way a modelled thread usually comes out of a CAD system: an ISO metric +/// ridge (60° flanks, P/8 crest) of `turns` turns wound on a core of the minor diameter, the +/// flanks ruled B-spline surfaces between helices, the crest and the core cylinders trimmed by +/// those helices. +fn ridge_thread(major: f64, pitch: f64, turns: f64) -> Body { + let h = pitch * 3f64.sqrt() / 2.0; + let (rc, rm) = (major / 2.0 - 5.0 * h / 8.0, major / 2.0); + let (foot, crest) = ((pitch - pitch / 4.0) / 2.0, pitch / 16.0); + let z0 = pitch; + let len = z0 + pitch * turns + pitch; + let end = TAU * turns; + let fl = helix(rc, z0 - foot, pitch, 0.0, end); + let fu = helix(rc, z0 + foot, pitch, 0.0, end); + let cl = helix(rm, z0 - crest, pitch, 0.0, end); + let cu = helix(rm, z0 + crest, pitch, 0.0, end); + let (fl0, fl1) = ends(&fl); + let (fu0, fu1) = ends(&fu); + let (cl0, cl1) = ends(&cl); + let (cu0, cu1) = ends(&cu); + let at = |r: f64, a: f64, z: f64| Vec3::new(r * a.cos(), r * a.sin(), z); + let mid = z0 + pitch / 2.0; + let out = |a: f64| Vec3::new(a.cos(), a.sin(), 0.0); + let flank = 0.5f64; + let faces = vec![ + FaceSpec { loops: vec![circle(rc, 0.0).iter().rev().map(rev).collect()], surface: SurfSpec::Plane, probe: None }, + FaceSpec { loops: vec![circle(rc, len)], surface: SurfSpec::Plane, probe: None }, + FaceSpec { + loops: vec![circle(rc, 0.0), circle(rc, len).iter().rev().map(rev).collect(), vec![line(fl0, fu0), fu.clone(), line(fu1, fl1), rev(&fl)]], + surface: SurfSpec::Revolved { + a: Vec3::new(rc, 0.0, -1.0), + b: Vec3::new(rc, 0.0, len + 1.0), + axis_origin: Vec3::ZERO, + axis: Vec3::Z, + angle: TAU, + }, + probe: Some((at(rc, PI, z0 * 0.5), out(PI))), + }, + FaceSpec { + loops: vec![vec![fl.clone(), line(fl1, cl1), rev(&cl), line(cl0, fl0)]], + surface: SurfSpec::Ruled { a: fl.clone(), b: cl.clone() }, + probe: Some((at((rc + rm) / 2.0, PI, mid - (foot + crest) / 2.0), out(PI) * flank - Vec3::Z * 0.866)), + }, + FaceSpec { + loops: vec![vec![cu.clone(), line(cu1, fu1), rev(&fu), line(fu0, cu0)]], + surface: SurfSpec::Ruled { a: fu.clone(), b: cu.clone() }, + probe: Some((at((rc + rm) / 2.0, PI, mid + (foot + crest) / 2.0), out(PI) * flank + Vec3::Z * 0.866)), + }, + FaceSpec { + loops: vec![vec![cl.clone(), line(cl1, cu1), rev(&cu), line(cu0, cl0)]], + surface: SurfSpec::Revolved { + a: Vec3::new(rm, 0.0, -1.0), + b: Vec3::new(rm, 0.0, len + 1.0), + axis_origin: Vec3::ZERO, + axis: Vec3::Z, + angle: TAU, + }, + probe: Some((at(rm, PI, mid), out(PI))), + }, + FaceSpec { loops: vec![vec![line(fl0, cl0), line(cl0, cu0), line(cu0, fu0), line(fu0, fl0)]], surface: SurfSpec::Plane, probe: None }, + FaceSpec { loops: vec![vec![line(fl1, fu1), line(fu1, cu1), line(cu1, cl1), line(cl1, fl1)]], surface: SurfSpec::Plane, probe: None }, + ]; + sew(&faces).unwrap() +} + +fn tol(b: &Body) -> f64 { + (b.size() * 1e-3).max(1e-3) +} + +struct Quality { + tris: usize, + slivers: usize, + flipped: usize, + worst_sag: f64, + min_angle: f64, +} + +/// Triangles thinner than 1e-3 of their longest side; triangles whose winding disagrees with +/// their vertex normals; the largest gap between a triangle's centre and the surface at the +/// centre's parameters. +fn quality(plan: &TessPlan, mesh: &Mesh) -> Quality { + let mut q = Quality { tris: mesh.triangles.len(), slivers: 0, flipped: 0, worst_sag: 0.0, min_angle: 180.0 }; + for t in &mesh.triangles { + let [a, b, c] = t.map(|i| mesh.positions[i as usize]); + let n = (b - a).cross(c - a); + let longest = (b - a).len().max((c - b).len()).max((a - c).len()); + let ang = |p: Vec3, q: Vec3, r: Vec3| { + ((q - p).normalized().unwrap_or(Vec3::X)).dot((r - p).normalized().unwrap_or(Vec3::X)).clamp(-1.0, 1.0).acos() + }; + let mina = ang(a, b, c).min(ang(b, c, a)).min(ang(c, a, b)).to_degrees(); + if mina < q.min_angle { + q.min_angle = mina; + } + if n.len() < 1e-3 * longest * longest { + q.slivers += 1; + continue; + } + let vn = t.iter().fold(Vec3::ZERO, |s, &i| s + mesh.normals[i as usize]); + if n.dot(vn) < 0.0 { + q.flipped += 1; + } + } + for f in &plan.faces { + let Some(s) = f.surface.and_then(|i| plan.surfaces.get(i as usize)) else { continue }; + let pos = |k: u32| f.fixed.get(k as usize).copied().unwrap_or_else(|| s.eval(f.uv[k as usize][0], f.uv[k as usize][1]).0); + for t in &f.triangles { + let uv = t.iter().fold([0.0, 0.0], |m, &k| [m[0] + f.uv[k as usize][0] / 3.0, m[1] + f.uv[k as usize][1] / 3.0]); + let mid = (pos(t[0]) + pos(t[1]) + pos(t[2])) * (1.0 / 3.0); + let (on, n) = s.eval(uv[0], uv[1]); + let sag = (on - mid).dot(n).abs(); + q.worst_sag = q.worst_sag.max(sag); + } + } + q +} + +/// Volume of the ridge thread: the core, and per turn of the ridge 2π ∫ r · width(r) dr (a +/// screw motion of a section in an axial plane). +fn ridge_volume(major: f64, pitch: f64, turns: f64) -> f64 { + let h = pitch * 3f64.sqrt() / 2.0; + let (rc, rm) = (major / 2.0 - 5.0 * h / 8.0, major / 2.0); + let (foot, crest) = ((pitch - pitch / 4.0) / 2.0, pitch / 16.0); + let len = pitch * turns + 2.0 * pitch; + let n = 1000; + let ridge: f64 = (0..n) + .map(|i| rc + (rm - rc) * (i as f64 + 0.5) / n as f64) + .map(|r| TAU * r * 2.0 * (foot + (crest - foot) * (r - rc) / (rm - rc)) * (rm - rc) / n as f64) + .sum(); + PI * rc * rc * len + ridge * turns +} + +/// The thread, as made and after a STEP round trip (which is how such a thread reaches us). +fn threads() -> Vec<(String, Body, f64, f64, f64)> { + let mut out = Vec::new(); + for (d, p, turns) in [(10.0, 1.5, 6.0), (12.0, 1.75, 8.0)] { + let b = ridge_thread(d, p, turns); + let step = crate::step_import(&crate::step_export(&[&b], "thread").unwrap()).unwrap().bodies.remove(0).body; + out.push((format!("M{d}x{p}"), b, d, p, turns)); + out.push((format!("M{d}x{p} from STEP"), step, d, p, turns)); + } + out +} + +/// Issue 1289: a modelled thread's faces turn several times around the shaft. The core and the +/// crest cylinders, trimmed by helices (and the core bounded by two circles without a seam +/// edge), had no display plan, and the reference mesher drew them with triangles across the +/// solid. They are now cut along a seam: the display mesh stays on the surfaces. +#[test] +fn thread_display_mesh_stays_on_its_surfaces() { + for (name, b, d, p, turns) in threads() { + let tol = tol(&b); + let plan = b.display_plan_angular(tol, DISPLAY_ANGLE).unwrap_or_else(|e| panic!("{name}: no display plan: {e}")); + let mesh = plan.evaluate(); + let q = quality(&plan, &mesh); + assert!(q.worst_sag <= 1.5 * tol, "{name}: a triangle lies {} mm off its surface (tolerance {tol})", q.worst_sag); + // Every vertex of the cylinders on its cylinder, every triangle centre within the + // tolerance inside it. + let h = p * 3f64.sqrt() / 2.0; + let (rc, rm) = (d / 2.0 - 5.0 * h / 8.0, d / 2.0); + for (ti, t) in mesh.triangles.iter().enumerate() { + let r = |v: Vec3| (v.x * v.x + v.y * v.y).sqrt(); + let [a, b2, c] = t.map(|i| mesh.positions[i as usize]); + let (ra, rb, rcn) = (r(a), r(b2), r(c)); + let centre = r((a + b2 + c) * (1.0 / 3.0)); + for want in [rc, rm] { + if [ra, rb, rcn].iter().all(|x| (x - want).abs() < 1e-4) { + assert!(centre > want - 2.0 * tol, "{name}: triangle {ti} on the r = {want} cylinder cuts in to {centre}"); + } + } + let end_disc = [a, b2, c].iter().all(|v| v.z.abs() < 1e-9) || [a, b2, c].iter().all(|v| (v.z - (turns + 2.0) * p).abs() < 1e-9); + assert!(end_disc || (centre > rc - 2.0 * tol && centre < rm + 1e-6), "{name}: triangle {ti} at radius {centre}, outside {rc}…{rm}"); + } + // A mesh within the tolerance of the surfaces is within area × tolerance of the volume. + let m = mesh.measure(); + let want = ridge_volume(d, p, turns); + assert!((m.volume - want).abs() < m.area * tol, "{name}: display mesh volume {}, the thread's {want}", m.volume); + // The reference mesh (picking, measuring, the viewport without GPU tessellation) is the same. + let reference = b.reference_display_mesh_angular(tol, DISPLAY_ANGLE).unwrap(); + assert_eq!(reference.triangles.len(), mesh.triangles.len(), "{name}"); + let measured = b.tessellate(tol * 0.1).unwrap().measure().volume; + assert!((measured - want).abs() < want * 1e-3, "{name}: measured volume {measured}, the thread's {want}"); + } +} + +/// Issue 1289: no triangle of the thread's display mesh is degenerate, a sliver or turned +/// against its surface's normal. +#[test] +fn thread_display_mesh_has_no_degenerate_or_flipped_triangles() { + for (name, b, ..) in threads() { + let plan = b.display_plan_angular(tol(&b), DISPLAY_ANGLE).unwrap_or_else(|e| panic!("{name}: no display plan: {e}")); + let mesh = plan.evaluate(); + let q = quality(&plan, &mesh); + assert!(q.tris > 1000, "{name}: {} triangles", q.tris); + assert_eq!(q.slivers, 0, "{name}: degenerate triangles"); + assert_eq!(q.flipped, 0, "{name}: triangles facing against their normals"); + assert!(q.min_angle > 1.0, "{name}: a triangle with a {:.3} degree corner", q.min_angle); + } +} + +/// Writes the M10 thread as STEP to `$THREAD_STEP` (for looking at it in the app). +#[test] +#[ignore] +fn write_thread_step() { + if let Ok(path) = std::env::var("THREAD_STEP") { + std::fs::write(path, crate::step_export(&[&ridge_thread(10.0, 1.5, 6.0)], "M10x1.5 thread").unwrap()).unwrap(); + } +}