diff --git a/crates/engine/src/tests.rs b/crates/engine/src/tests.rs index e9bef257..000b8b58 100644 --- a/crates/engine/src/tests.rs +++ b/crates/engine/src/tests.rs @@ -307,6 +307,7 @@ mod blend_fit; mod boolean_diagnostics; mod component_frames; mod curved_seam_round; +mod flush_join_loft; mod gltf_import; mod import_component; mod large_step_import; diff --git a/crates/engine/src/tests/flush_join_loft.rs b/crates/engine/src/tests/flush_join_loft.rs new file mode 100644 index 00000000..18598e5d --- /dev/null +++ b/crates/engine/src/tests/flush_join_loft.rs @@ -0,0 +1,99 @@ +//! A flange joined flush with the end of a bent, tapered lofted half tube (#1391). +//! +//! The tube is a half annulus with a 6 mm wall, lofted through three sections along a bend of radius 300 mm and 20 degrees, its radius +//! growing from 60 to 75 mm. The flange is a thicker half annulus (3 mm into the wall, 25 mm beyond it), 8 mm thick, extruded back from +//! the tube's end plane and joined. Its front face lies in that plane, coplanar with the tube's end face, and the neighbours of that +//! face are faces of the loft, which curve and lean: the faces cannot be pushed apart, and the join failed with "Union: no result". +//! (With a 40 degree bend the retry with the tool displaced joins it; with this bend that retry fails too.) + +use std::f64::consts::PI; + +use super::*; + +/// Radius of the bend of the tube's axis, and the turn of it (radians). +const BEND: f64 = 300.0; +const TURN: f64 = 20.0 * PI / 180.0; +/// Inner radius of the tube at its start and at its end, and the wall. +const R0: f64 = 60.0; +const R1: f64 = 75.0; +const WALL: f64 = 6.0; +/// How far the flange reaches into the tube's wall, how far beyond its outer face, and how thick it is. +const LAND_IN: f64 = 3.0; +const LAND_OUT: f64 = 25.0; +const THICK: f64 = 8.0; + +/// Sketch coordinates are given to nine decimals, as a design's own would be. +fn r9(x: f64) -> f64 { + (x * 1e9).round() / 1e9 +} + +/// The plane of the section at the bend angle `a`, moved `over` along the axis (the sketch's x is up, its y across the tube). +fn plane(a: f64, over: f64) -> Value { + let (tx, ty) = (a.cos(), a.sin()); + json!({ + "origin": [r9(BEND * a.sin() + over * tx), r9(BEND * (1.0 - a.cos()) + over * ty), 0], + "x_dir": [0, 0, 1], + "y_dir": [r9(ty), r9(-tx), 0], + }) +} + +/// A half annulus (the half above z = 0) on `plane`, between the radii `r_in` and `r_out`. +fn half_annulus(s: &mut Session, name: &str, plane: Value, r_in: f64, r_out: f64) { + let (r_in, r_out) = (r9(r_in), r9(r_out)); + run(s, "sketch.create", json!({"plane": plane, "name": name})); + run(s, "sketch.arc.three_point", json!({"start": [0, -r_out], "end": [0, r_out], "through": [r_out, 0]})); + run(s, "sketch.line", json!({"points": [[0, r_out], [0, r_in]]})); + run(s, "sketch.arc.three_point", json!({"start": [0, r_in], "end": [0, -r_in], "through": [r_in, 0]})); + run(s, "sketch.line", json!({"points": [[0, -r_in], [0, -r_out]]})); + run(s, "sketch.finish", json!({})); +} + +/// The tube: a loft through three sections. +fn tube(s: &mut Session) { + for i in 0..3 { + let t = i as f64 / 2.0; + let r_in = R0 + (R1 - R0) * t; + half_annulus(s, &format!("S{i}"), plane(TURN * t, 0.0), r_in, r_in + WALL); + } + let sections: Vec = (0..3).map(|i| json!({"sketch": format!("S{i}")})).collect(); + run(s, "solid.loft", json!({"sections": sections, "operation": "new", "name": "TubeLoft", "body_name": "Tube"})); +} + +/// The flange, joined to the tube. Its front face is `over` mm past the tube's end plane (0: in it). +fn flange(s: &mut Session, over: f64) { + half_annulus(s, "Flange", plane(TURN, over), R1 + LAND_IN, R1 + WALL + LAND_OUT); + let join = json!({"sketch": "Flange", "distance": THICK + over, "direction": "negative", "operation": "join", "targets": ["Tube"], "name": "JoinFlange"}); + run(s, "solid.extrude", join); +} + +#[test] +fn a_flange_flush_with_the_end_of_a_bent_tapered_tube_joins() { + // The flange's face, the half annulus, has this area. + let (r_in, r_out) = (R1 + LAND_IN, R1 + WALL + LAND_OUT); + let footprint = PI / 2.0 * (r_out * r_out - r_in * r_in); + + // The reference: with the flange 2 mm past the end plane nothing is coplanar with the tube's end face (that join is made by another + // route), and what it adds beyond that plane is the footprint, 2 mm deep. + let mut s = Session::default(); + tube(&mut s); + flange(&mut s, 2.0); + let reference = volume(&mut s) - 2.0 * footprint; + + let mut s = Session::default(); + tube(&mut s); + flange(&mut s, 0.0); + let v = volume(&mut s); + assert!(rel(v, reference) < 3e-4, "{v} against {reference}"); + + // One valid body in one piece, and nothing of it beyond the end plane (to a micrometre: the curved faces end on it within the + // tolerance of their edges): the face the flange and the tube share is one flat face. + let state = s.model.state(); + assert_eq!(state.bodies.len(), 1); + let body = &state.bodies[0].body; + assert!(body.validity().is_empty(), "{:?}", body.validity()); + assert_eq!(body.lumps().unwrap().len(), 1); + let (origin, normal) = (Vec3::new(BEND * TURN.sin(), BEND * (1.0 - TURN.cos()), 0.0), Vec3::new(TURN.cos(), TURN.sin(), 0.0)); + let mesh = body.tessellate(0.05).unwrap(); + let beyond = mesh.positions.iter().map(|p| (*p - origin).dot(normal)).fold(f64::MIN, f64::max); + assert!(beyond < 1e-3, "the join reaches {beyond} mm past the end plane"); +} diff --git a/crates/kernel/src/clearance.rs b/crates/kernel/src/clearance.rs index a2843ff2..c33c9c96 100644 --- a/crates/kernel/src/clearance.rs +++ b/crates/kernel/src/clearance.rs @@ -110,6 +110,8 @@ pub(crate) struct Clearance<'a> { /// How far a curved face may stand off its triangles. sag: f64, tol: f64, + /// Whether curved faces are judged by their exact bounds (see [`Clearance::with_exact_bounds`]). + exact: bool, } /// Pairs of triangles one question may compare; beyond it there is no answer. @@ -119,7 +121,18 @@ impl<'a> Clearance<'a> { /// `chord`: the chord error of the mesh the body is asked on. pub(crate) fn new(body: &'a Body, chord: f64) -> Option> { let mesh = body.tessellate_with(chord, true).ok()?; - (!mesh.triangles.is_empty()).then(|| Clearance { body, mesh, sag: 2.0 * chord, tol: (body.size() * 1e-7).max(1e-9) }) + (!mesh.triangles.is_empty()).then(|| Clearance { body, mesh, sag: 2.0 * chord, tol: (body.size() * 1e-7).max(1e-9), exact: false }) + } + + /// Judge curved faces by their exact bounds where `clear_of_slab` is asked, not by the chord error each way. A B-spline or NURBS + /// face lies in the hull of its control net, so a face whose net ends before a region begins (a loft's side running out in the + /// plane the region starts from) can't reach into it, however its triangles stand off it. A face of a sphere or a cylinder + /// stands off its triangles away from the centre or the axis only, so a triangle that lies clear of the region and leans away + /// from it is clear (a dome standing on the plane the region starts from). Without this such a face counts as reaching in by + /// the chord error. + pub(crate) fn with_exact_bounds(mut self, on: bool) -> Self { + self.exact = on; + self } /// How the body's material sits in the space that `faces` move through when they go along @@ -194,18 +207,44 @@ impl<'a> Clearance<'a> { /// a unit vector), anywhere? pub(crate) fn clear_of_slab(&self, n: Vec3, d: f64, lo: f64, hi: f64) -> bool { let mut reach: HashMap = HashMap::new(); + let mut round: HashMap> = HashMap::new(); self.mesh.triangles.iter().zip(&self.mesh.tri_face).all(|(t, face)| { let Some(tri) = self.mesh.tri(t) else { return false }; - let pad = match *reach.entry(*face).or_insert_with(|| self.reach(*face, n, d + lo, d + hi)) { + // How far the face may stand off its triangle towards larger and towards smaller heights. + let (up, down) = match *reach.entry(*face).or_insert_with(|| self.reach(*face, n, d + lo, d + hi)) { Reach::Never => return true, - Reach::Around | Reach::Beyond { .. } => self.sag, - Reach::Exact | Reach::Sideways => 0.0, + Reach::Around => (self.sag, self.sag), + Reach::Beyond { .. } => match self.exact.then(|| round.entry(*face).or_insert_with(|| self.round_of(*face))).and_then(|s| s.as_ref()) + { + Some(s) => self.bulge(s, (tri[0] + tri[1] + tri[2]) / 3.0, n), + None => (self.sag, self.sag), + }, + Reach::Exact | Reach::Sideways => (0.0, 0.0), }; let heights = tri.map(|p| n.dot(p) - d); - heights.iter().all(|s| s + pad < lo) || heights.iter().all(|s| s - pad > hi) + heights.iter().all(|s| s + up < lo) || heights.iter().all(|s| s - down > hi) }) } + /// The round surface (sphere, cylinder, cone, torus) a face lies on, if it is one. + fn round_of(&self, face: u32) -> Option { + let f = self.body.solid.face_iter().nth(face as usize)?; + surf_of(f, self.tol * 100.0) + } + + /// How far the surface may stand off a triangle with its centroid at `p`, towards larger and towards smaller heights along `n`: a + /// sphere or a cylinder lies radially outward from its triangles, within the chord error. + fn bulge(&self, s: &Surf, p: Vec3, n: Vec3) -> (f64, f64) { + let radial = match *s { + Surf::Sphere { c, .. } => p - c, + Surf::Cylinder { o, a, .. } => (p - o) - a * (p - o).dot(a), + _ => return (self.sag, self.sag), + }; + let Some(out) = radial.normalized() else { return (self.sag, self.sag) }; + let lean = out.dot(n); + (self.sag * lean.max(0.0), self.sag * (-lean).max(0.0)) + } + /// How the face's triangles stand for it in a region between the heights `lo` and `hi` /// along `along` (measured as `along · p`). fn reach(&self, face: u32, along: Vec3, lo: f64, hi: f64) -> Reach { @@ -218,6 +257,12 @@ impl<'a> Clearance<'a> { if matches!(surface, mt::Surface::Plane(_)) { return Ok(Reach::Exact); } + if self.exact + && let Some((least, greatest)) = crate::planar_faces::control_net_span(&surface, along) + && (greatest <= lo || least >= hi) + { + return Ok(Reach::Never); + } let round = surf_of(f, self.tol * 100.0); let straight = matches!(round, Some(Surf::Cylinder { a, .. }) if a.dot(along).abs() > 1.0 - 1e-12) || crate::prism::wall_along(&surface, along, self.tol * 10.0, &mut 262_144); diff --git a/crates/kernel/src/coplanar.rs b/crates/kernel/src/coplanar.rs index b834963d..530cc4a4 100644 --- a/crates/kernel/src/coplanar.rs +++ b/crates/kernel/src/coplanar.rs @@ -240,6 +240,483 @@ pub(crate) fn push_faces_exact(b: &Body, n: Vec3, d: f64, delta: f64) -> Result< }) } +/// Move the planar faces of `b` that lie in the plane (n, d), n their outward normal, out by `delta`, whatever their +/// neighbours are. [`push_faces`] stretches the neighbours, which must stand square to the face along straight edges (a +/// box, an extruded outline); the neighbours of a lofted tube's end cap curve and lean. Here nothing is stretched: the +/// faces move, and each edge one shares with a face that stays gets a new face between them, the edge swept along `n` +/// (a strip). The neighbours keep their surfaces and meet the strips at a crease along the old edge, so the body is +/// the old one with a short prism of the faces' outline on each of them. Exact, with no boolean. The result has no +/// history. +pub(crate) fn grow_faces_exact(b: &Body, n: Vec3, d: f64, delta: f64) -> Result { + use mt::{BoundedCurve, ParametricCurve, ParametricSurface3D, SearchNearestParameter}; + let size = b.size(); + let tol = (size * 1e-7).max(1e-9); + let fail = |m: &str| KernelError::Failed(format!("growing a face: {m}")); + guard("grow face", || { + let solid = b.deep_copy(); + let shift = n * delta; + let mat = mt::Matrix4::from_translation(v3(shift)); + let mut shells = Vec::new(); + for sh in solid.boundaries() { + let faces: Vec = sh.face_iter().cloned().collect(); + let on_plane = |f: &mt::Face| match f.oriented_surface() { + mt::Surface::Plane(p) => { + let m = p.normal(); + let m = Vec3::new(m.x, m.y, m.z); + m.dot(n) > 1.0 - 1e-9 && (n.dot(from_p3(p.origin())) - d).abs() < tol * 10.0 + } + _ => false, + }; + let moving: Vec = (0..faces.len()).filter(|i| faces.get(*i).is_some_and(on_plane)).collect(); + if moving.is_empty() { + shells.push(sh.clone()); + continue; + } + let stays = |i: &usize| !moving.contains(i); + // The faces each edge bounds. + let mut bounds: HashMap> = HashMap::new(); + for (i, f) in faces.iter().enumerate() { + for e in f.edge_iter() { + bounds.entry(e.id()).or_default().push(i); + } + } + // The moved copies of the moving faces' vertices and edges. + let mut vmap: HashMap = HashMap::new(); + let mut emap: HashMap = HashMap::new(); + for i in &moving { + let f = faces.get(*i).ok_or_else(|| fail("face"))?; + for v in f.vertex_iter() { + vmap.entry(v.id()).or_insert_with(|| builder::vertex(p3(from_p3(v.point()) + shift))); + } + } + for i in &moving { + let f = faces.get(*i).ok_or_else(|| fail("face"))?; + for e in f.edge_iter() { + if emap.contains_key(&e.id()) { + continue; + } + let (a, bb) = ( + vmap.get(&e.absolute_front().id()).ok_or_else(|| fail("vertex"))?, + vmap.get(&e.absolute_back().id()).ok_or_else(|| fail("vertex"))?, + ); + emap.insert(e.id(), mt::Edge::new(a, bb, mt::Transformed::transformed(&e.curve(), mat))); + } + } + // The moved copy of an edge, running the way the edge does. + let moved = |e: &mt::Edge| -> Result { + let ne = emap.get(&e.id()).ok_or_else(|| fail("edge"))?; + Ok(if e.front() == e.absolute_front() { ne.clone() } else { ne.inverse() }) + }; + // A strip on every edge between a moving face and one that stays; the edges from a vertex to its moved copy are + // shared by the two strips that meet there. + let mut lateral: HashMap = HashMap::new(); + let mut strips: Vec = Vec::new(); + for i in &moving { + let f = faces.get(*i).ok_or_else(|| fail("face"))?; + for w in f.boundaries() { + for e in w.edge_iter() { + let Some(&beside) = bounds.get(&e.id()).and_then(|fs| fs.iter().find(|j| stays(j))) else { continue }; + // The moving face runs along the edge from p to q; the strip runs along it the same way. + let (p, q) = (e.front().clone(), e.back().clone()); + let mut lat = |v: &mt::Vertex| -> Result { + let to = vmap.get(&v.id()).ok_or_else(|| fail("vertex"))?; + Ok(lateral.entry(v.id()).or_insert_with(|| builder::line(v, to)).clone()) + }; + let (lp, lq) = (lat(&p)?, lat(&q)?); + let wire: mt::Wire = vec![e.clone(), lq, moved(e)?.inverse(), lp.inverse()].into(); + // Its surface is the edge's curve swept along `n`, facing the way the face beside it does; a straight edge + // sweeps to a plane (truck makes that a B-spline patch, which does not mesh). + let curve = e.curve(); + let (t0, t1) = curve.range_tuple(); + let mid = from_p3(curve.subs(0.5 * (t0 + t1))); + let old = faces.get(beside).ok_or_else(|| fail("face"))?.oriented_surface(); + let Some((u0, v0)) = old.search_nearest_parameter(p3(mid), None, 100) else { + return Err(fail("which way a strip faces")); + }; + let outward = old.normal(u0, v0); + let outward = Vec3::new(outward.x, outward.y, outward.z); + // (a loft's straight section edge is a B-spline of degree one, so ask the points, not the type) + let straight = [0.2, 0.4, 0.6, 0.8].iter().all(|t| { + let x = from_p3(curve.subs(t0 + t * (t1 - t0))); + let chord = from_p3(q.point()) - from_p3(p.point()); + let off = (x - from_p3(p.point())).cross(chord).len(); + off <= tol * 10.0 * chord.len() + }); + let s = if straight { + let dir = (from_p3(q.point()) - from_p3(p.point())).normalized().ok_or_else(|| fail("a zero edge"))?; + let mut m = dir.cross(n).normalized().ok_or_else(|| fail("an edge along the normal"))?; + if m.dot(outward) < 0.0 { + m = -m; + } + let o = from_p3(p.point()); + mt::Surface::Plane(mt::Plane::new(p3(o), p3(o + dir), p3(o + m.cross(dir)))) + } else { + let swept: mt::Face = builder::tsweep(e, v3(shift)); + let mut s = swept.surface(); + let Some((u1, v1)) = s.search_nearest_parameter(p3(mid + shift * 0.5), None, 100) else { + return Err(fail("which way a strip faces")); + }; + let now = s.normal(u1, v1); + if outward.dot(Vec3::new(now.x, now.y, now.z)) < 0.0 { + s = mt::Invertible::inverse(&s); + } + s + }; + strips.push(mt::Face::try_new(vec![wire], s).map_err(|e| fail(&e.to_string()))?); + } + } + } + let mut out: Vec = Vec::new(); + for (i, f) in faces.iter().enumerate() { + if !moving.contains(&i) { + out.push(f.clone()); + continue; + } + let wires = f + .absolute_boundaries() + .iter() + .map(|w| w.edge_iter().map(&moved).collect::>>().map(mt::Wire::from)) + .collect::>>()?; + let mut nf = mt::Face::try_new(wires, mt::Transformed::transformed(&f.surface(), mat)).map_err(|e| fail(&e.to_string()))?; + if !f.orientation() { + nf.invert(); + } + out.push(nf); + } + out.extend(strips); + shells.push(out.into()); + } + let s = Solid::try_new(shells).map_err(|e| fail(&e.to_string()))?; + let r = Body::new(s)?; + let problems = r.validity(); + if !problems.is_empty() { + return Err(fail(&format!("an unsound result ({})", problems.iter().take(2).cloned().collect::>().join("; ")))); + } + Ok(r) + }) +} + +/// An edge's points from its start to its end along its curve, `k + 1` of them. +fn edge_samples(e: &mt::Edge, k: u32) -> Vec { + use mt::{BoundedCurve, ParametricCurve}; + let c = e.curve(); + let (t0, t1) = c.range_tuple(); + (0..=k).map(|i| from_p3(c.subs(t0 + (t1 - t0) * f64::from(i) / f64::from(k)))).collect() +} + +/// Where the edges of `b` cross the plane (n, d), one point for every time (edges are cut there before a body is trimmed). +fn plane_crossings(b: &Body, n: Vec3, d: f64, tol: f64) -> Vec { + use mt::{BoundedCurve, ParametricCurve}; + const STEPS: u32 = 64; + let mut out = Vec::new(); + for e in b.solid.edge_iter() { + let c = e.curve(); + let (t0, t1) = c.range_tuple(); + let at = |t: f64| from_p3(c.subs(t)); + let h = |t: f64| n.dot(at(t)) - d; + let ts: Vec = (0..=STEPS).map(|i| t0 + (t1 - t0) * f64::from(i) / f64::from(STEPS)).collect(); + // The samples clear of the plane, and between each two that lie on different sides the crossing. + let clear: Vec<(f64, f64)> = ts.iter().map(|t| (*t, h(*t))).filter(|(_, x)| x.abs() > tol).collect(); + for pair in clear.windows(2) { + let [(ta, ha), (tb, hb)] = pair else { continue }; + if (*ha < 0.0) == (*hb < 0.0) { + continue; + } + let (mut lo, mut hi, mut hlo) = (*ta, *tb, *ha); + for _ in 0..80 { + let mid = 0.5 * (lo + hi); + let hm = h(mid); + if (hm < 0.0) == (hlo < 0.0) { + lo = mid; + hlo = hm; + } else { + hi = mid; + } + } + out.push(at(0.5 * (lo + hi))); + } + } + out +} + +/// Which side of the plane (n, d) an edge lies on, once edges that cross it are cut: -1 near, +1 far, 0 in it. +fn edge_side(e: &mt::Edge, n: Vec3, d: f64, tol: f64) -> Result { + let hs: Vec = edge_samples(e, 8).iter().map(|p| n.dot(*p) - d).collect(); + let (lo, hi) = hs.iter().fold((f64::MAX, f64::MIN), |(l, h), x| (l.min(*x), h.max(*x))); + if hi <= tol && lo >= -tol { + Ok(0) + } else if hi <= tol { + Ok(-1) + } else if lo >= -tol { + Ok(1) + } else { + Err(KernelError::Failed("trimming at a plane: an edge still crosses it".into())) + } +} + +/// Whether the point lies on the surface, to within `tol`. +fn on_surface(s: &mt::Surface, p: Vec3, tol: f64) -> bool { + use mt::{ParametricSurface, SearchNearestParameter}; + if let mt::Surface::Plane(pl) = s { + let m = pl.normal(); + return (Vec3::new(m.x, m.y, m.z).dot(p - from_p3(pl.origin()))).abs() <= tol; + } + s.search_nearest_parameter(p3(p), None, 100).is_some_and(|(u, v)| from_p3(s.subs(u, v)).dist(p) <= tol) +} + +/// The edge that closes the boundary of the flat face `face` from the vertex `from` to the vertex `to` (both in the plane) where a run of +/// its edges beyond the plane went: the segment between them, which lies in the face. (A curved face that crosses the plane has been +/// divided along it already.) Nothing when they are the same vertex. +fn bridge(face: &mt::Face, from: &mt::Vertex, to: &mt::Vertex, tol: f64) -> Result> { + if from.id() == to.id() { + return Ok(Vec::new()); + } + let (p, q) = (from_p3(from.point()), from_p3(to.point())); + let surface = face.oriented_surface(); + if [0.25, 0.5, 0.75].iter().all(|t| on_surface(&surface, p + (q - p) * *t, tol * 10.0)) { + Ok(vec![builder::line(from, to)]) + } else { + Err(KernelError::Failed("trimming at a plane: a curved face crosses the plane where it could not be divided".into())) + } +} + +/// One boundary of the flat face `face`, which crosses the plane (n, d), clipped to the near side: the edges beyond the plane go and the +/// segment in the plane from where the boundary leaves it to where it comes back takes their place. `None` when nothing of the boundary +/// is on the near side. +fn clip_wire(face: &mt::Face, w: &mt::Wire, n: Vec3, d: f64, tol: f64) -> Result>> { + let fail = |m: &str| KernelError::Failed(format!("trimming at a plane: {m}")); + let mut edges: Vec<(mt::Edge, i8)> = Vec::new(); + for e in w.edge_iter() { + let side = edge_side(e, n, d, tol)?; + edges.push((e.clone(), side)); + } + if !edges.iter().any(|(_, side)| *side < 0) { + return Ok(None); + } + let m = edges.len(); + let gone = |i: usize| edges.get(i % m).is_some_and(|(_, side)| *side > 0); + // Begin where a run of edges beyond the plane ends, so that no run is split between the end and the start. + let Some(start) = (0..m).find(|i| !gone(*i) && gone(*i + m - 1)) else { + return Ok(Some(edges.into_iter().map(|(e, _)| e).collect())); + }; + let mut out: Vec = Vec::new(); + let mut skipping = false; + for k in 0..m { + let Some((e, side)) = edges.get((start + k) % m) else { continue }; + if *side > 0 { + skipping = true; + continue; + } + if skipping { + let from = out.last().ok_or_else(|| fail("a run beyond the plane first"))?.back().clone(); + out.extend(bridge(face, &from, e.front(), tol)?); + skipping = false; + } + out.push(e.clone()); + } + if skipping { + let from = out.last().ok_or_else(|| fail("nothing kept"))?.back().clone(); + let to = out.first().ok_or_else(|| fail("nothing kept"))?.front().clone(); + out.extend(bridge(face, &from, &to, tol)?); + } + Ok(Some(out)) +} + +/// `b` with everything beyond the plane (n, d) cut off, for a body that was grown or pushed past the plane (the faces that +/// [`grow_faces_exact`] and [`push_faces`] moved) and joined there. A cut with a boolean would run along the edges where the grown +/// faces meet their neighbours, which the mesh intersection cannot take. Here the edges that cross the plane are cut at it, the +/// faces beyond it go, each face that crosses it gets one edge in the plane where its far part was, and the rim they leave (edges +/// lying in the plane) is closed with flat faces. The result has no history. +pub(crate) fn trim_exact(b: &Body, n: Vec3, d: f64) -> Result { + let fail = |m: &str| KernelError::Failed(format!("trimming at a plane: {m}")); + let tol = (b.size() * 1e-6).max(1e-9); + // A curved face that crosses the plane (a sphere or a cylinder stretched past it) is divided along the plane where they meet, exactly + // ([`crate::splitface::split_faces`]); the flat ones are clipped below. + let curved: Vec = b + .solid + .face_iter() + .enumerate() + .filter(|(_, f)| !matches!(f.oriented_surface(), mt::Surface::Plane(_))) + .filter(|(_, f)| { + let (mut near, mut far) = (false, false); + for e in f.edge_iter() { + for p in edge_samples(&e, 8) { + near |= n.dot(p) - d < -tol; + far |= n.dot(p) - d > tol; + } + } + near && far + }) + .map(|(i, _)| i) + .collect(); + let divided = if curved.is_empty() { + b.clone() + } else { + let u = n.any_perp(); + let plane = solidcraft_geom::Plane::new(n * d, u, n.cross(u)).ok_or_else(|| fail("the plane"))?; + crate::splitface::split_faces(b, &curved, &crate::SplitTool::Plane(plane))? + }; + let cut = crate::splitface::split_edges_at(÷d, &plane_crossings(÷d, n, d, tol))?; + guard("trim at a plane", || { + let mut kept: Vec = Vec::new(); + for f in cut.solid.face_iter() { + let mut sides: Vec = Vec::new(); + for e in f.edge_iter() { + sides.push(edge_side(&e, n, d, tol)?); + } + match (sides.contains(&-1), sides.contains(&1)) { + (true, false) => kept.push(f.clone()), + (false, true) => {} + // A face lying in the plane (a sliver the join left there) goes with the others: the rim it leaves is closed below. + (false, false) => {} + (true, true) => { + let mut wires: Vec = Vec::new(); + for (i, w) in f.absolute_boundaries().iter().enumerate() { + match clip_wire(f, w, n, d, tol)? { + Some(edges) => wires.push(edges.into()), + None if i == 0 => return Err(fail("a face with no edge left on the near side")), + None => {} + } + } + let mut nf = mt::Face::try_new(wires, f.surface()).map_err(|x| fail(&x.to_string()))?; + if !f.orientation() { + nf.invert(); + } + kept.push(nf); + } + } + } + close_rim(kept, n, d, tol) + }) +} + +/// The body whose faces are `faces`, which leave an opening in the plane (n, d): the edges only one face is on (all lying in the plane) +/// are closed with flat faces. Counterclockwise loops about the normal are outlines, the others holes in the outline around them. +fn close_rim(mut faces: Vec, n: Vec3, d: f64, tol: f64) -> Result { + let fail = |m: &str| KernelError::Failed(format!("closing the cut: {m}")); + let mut uses: HashMap = HashMap::new(); + for f in &faces { + for e in f.edge_iter() { + *uses.entry(e.id()).or_default() += 1; + } + } + if uses.values().any(|c| *c > 2) { + return Err(fail("an edge on more than two faces")); + } + let mut open: Vec = Vec::new(); + for f in &faces { + for w in f.boundaries() { + for e in w.edge_iter() { + if uses.get(&e.id()) == Some(&1) { + if edge_samples(e, 8).iter().any(|p| (n.dot(*p) - d).abs() > tol) { + return Err(fail("an edge left open off the plane")); + } + // The cap runs along it the other way. + open.push(e.inverse()); + } + } + } + } + if open.is_empty() { + return Err(fail("nothing was cut off")); + } + let mut loops: Vec> = Vec::new(); + while let Some(first) = open.pop() { + let start = first.front().id(); + let mut at = first.back().id(); + let mut chain = vec![first]; + while at != start { + let next: Vec = (0..open.len()).filter(|i| open.get(*i).is_some_and(|e| e.front().id() == at)).collect(); + let [i] = next[..] else { return Err(fail("the rim does not close up, or pinches")) }; + let e = open.remove(i); + at = e.back().id(); + chain.push(e); + } + loops.push(chain); + } + let u = n.any_perp(); + let v = n.cross(u); + let origin = n * d; + let polys: Vec> = loops + .iter() + .map(|chain| { + chain + .iter() + .flat_map(|e| { + let mut pts = edge_samples(e, 8); + if !e.orientation() { + pts.reverse(); + } + pts.pop(); + pts + }) + .map(|p| ((p - origin).dot(u), (p - origin).dot(v))) + .collect() + }) + .collect(); + let areas: Vec = polys + .iter() + .map(|poly| { + let m = poly.len(); + (0..m) + .map(|i| match (poly.get(i), poly.get((i + 1) % m)) { + (Some(a), Some(b)) => a.0 * b.1 - b.0 * a.1, + _ => 0.0, + }) + .sum::() + * 0.5 + }) + .collect(); + let inside = |poly: &[(f64, f64)], q: (f64, f64)| -> bool { + let m = poly.len(); + (0..m) + .filter(|i| match (poly.get(*i), poly.get((*i + 1) % m)) { + (Some(a), Some(b)) => (a.1 > q.1) != (b.1 > q.1) && q.0 < (b.0 - a.0) * (q.1 - a.1) / (b.1 - a.1) + a.0, + _ => false, + }) + .count() + % 2 + == 1 + }; + let outers: Vec = (0..loops.len()).filter(|i| areas.get(*i).is_some_and(|a| *a > 0.0)).collect(); + let mut holes_of: HashMap> = HashMap::new(); + for i in (0..loops.len()).filter(|i| !outers.contains(i)) { + let q = polys.get(i).and_then(|p| p.first()).copied().ok_or_else(|| fail("an empty loop"))?; + let around = outers + .iter() + .filter(|o| polys.get(**o).is_some_and(|p| inside(p, q))) + .min_by(|x, y| areas.get(**x).partial_cmp(&areas.get(**y)).unwrap_or(std::cmp::Ordering::Equal)) + .ok_or_else(|| fail("a hole outside every outline"))?; + holes_of.entry(*around).or_default().push(i); + } + let surface = mt::Surface::Plane(mt::Plane::new(p3(origin), p3(origin + u), p3(origin + v))); + for o in &outers { + let mut wires: Vec = vec![loops.get(*o).cloned().ok_or_else(|| fail("a loop"))?.into()]; + for h in holes_of.get(o).into_iter().flatten() { + wires.push(loops.get(*h).cloned().ok_or_else(|| fail("a loop"))?.into()); + } + faces.push(mt::Face::try_new(wires, surface.clone()).map_err(|x| fail(&x.to_string()))?); + } + let shell: mt::Shell = faces.into(); + let solid = Solid::try_new(crate::polybool::components_in_order(&shell)).map_err(|x| fail(&x.to_string()))?; + let r = Body::new(solid)?; + let problems = r.validity(); + if !problems.is_empty() { + return Err(fail(&format!("an unsound result ({})", problems.iter().take(2).cloned().collect::>().join("; ")))); + } + Ok(r) +} + +/// Move the planar faces of `x` in the plane (n, d) out by `delta`: stretched ([`push_faces`]) where their neighbours +/// allow it, else, when `grow`, grown ([`grow_faces_exact`]). Whether it grew them (the body then has no history). +fn move_out(x: &Body, n: Vec3, d: f64, delta: f64, grow: bool) -> Result<(Body, bool)> { + match push_faces(x, n, d, delta) { + Ok(r) => Ok((r, false)), + Err(e) if grow => grow_faces_exact(x, n, d, delta).map(|r| (r, true)).map_err(|_| e), + Err(e) => Err(e), + } +} + /// A box covering everything beyond the plane (n, d) by up to `depth`. fn slab_beyond(n: Vec3, d: f64, depth: f64, extent: f64) -> Result { let u = n.any_perp(); @@ -326,6 +803,23 @@ pub(crate) const NO_ROOM: &str = "coincident faces were left in place: pushing t /// coincident faces or the faces can't be pushed; the error [`NO_ROOM`] when a face would /// move through something that changes the result. pub(crate) fn boolean_apart(a: &Body, b: &Body, op: BoolOp) -> Option>> { + apart_with(a, b, op, false) +} + +/// [`boolean_apart`] for a join whose coincident faces cannot be pushed, because the neighbours of one of them curve +/// (a flange flush with the end of a lofted tube): such a face is grown instead ([`grow_faces_exact`]) and the join's excess beyond +/// the plane is cut off again by editing faces ([`trim_exact`]), where a boolean would fail on the edges the grown faces meet their +/// neighbours along. The planar faces that sit on B-spline surfaces (a loft's flat sides) are put on planes first so that they +/// are seen. Asked after the other routes have failed, so that what they take stays as it was. `None` when it does not apply. +pub(crate) fn boolean_grown(a: &Body, b: &Body, op: BoolOp) -> Option>> { + if op != BoolOp::Union { + return None; + } + let (pa, pb) = (crate::planar_faces::planarized(a), crate::planar_faces::planarized(b)); + apart_with(pa.as_ref().unwrap_or(a), pb.as_ref().unwrap_or(b), op, true) +} + +fn apart_with(a: &Body, b: &Body, op: BoolOp, grow: bool) -> Option>> { let size = a.size().max(b.size()); let tol = (size * 1e-7).max(1e-9); let coincident = coincident_planes(a, b); @@ -341,7 +835,7 @@ pub(crate) fn boolean_apart(a: &Body, b: &Body, op: BoolOp) -> Option Result> { // The target's faces, meshed when the first face asks what it would move through. let faces_of_a = std::cell::OnceCell::new(); - let target = || faces_of_a.get_or_init(|| Clearance::new(a, chord)).as_ref().ok_or_else(no_room); + let target = || faces_of_a.get_or_init(|| Clearance::new(a, chord).map(|c| c.with_exact_bounds(grow))).as_ref().ok_or_else(no_room); let (mut a2, mut b2) = (a.clone(), b.clone()); // What the pushed bodies are made of: their histories over [a] and [b]. let (mut ha, mut hb) = (crate::boolean_record::kept(a, 0, 1), crate::boolean_record::kept(b, 0, 1)); @@ -351,7 +845,7 @@ pub(crate) fn boolean_apart(a: &Body, b: &Body, op: BoolOp) -> Option Option Option Body { + let ball = sphere(Vec3::ZERO, r).unwrap(); + let below = box_solid(Vec3::new(-2.0 * r, -2.0 * r, -2.0 * r), Vec3::new(2.0 * r, 2.0 * r, 0.0)).unwrap(); + boolean(&ball, &below, BoolOp::Cut).unwrap().unwrap() +} + +/// A hemisphere (r 20) and a block (x 5 to 35, y -30 to 30, z 0 to 25) that both stand on z = 0, the block's foot over most of the +/// hemisphere's flat face. The neighbour of that face is the sphere, which cannot be stretched, so its face is grown instead (a +/// strip along its rim) and what the join leaves below the plane is cut off by editing faces. The block takes the part of the dome +/// beyond x = 5, half of a cap 15 high: the union is the dome and the block less that part. +#[test] +fn a_block_flush_with_the_flat_face_of_a_hemisphere_joins_exactly() { + let (r, x0) = (20.0, 5.0); + let block = box_solid(Vec3::new(x0, -30.0, 0.0), Vec3::new(35.0, 30.0, 25.0)).unwrap(); + let dome = hemisphere(r); + let h = r - x0; + let expected = 2.0 / 3.0 * std::f64::consts::PI * r.powi(3) - std::f64::consts::PI * h * h * (3.0 * r - h) / 6.0 + 30.0 * 60.0 * 25.0; + for (a, b) in [(&dome, &block), (&block, &dome)] { + let joined = crate::coplanar::boolean_grown(a, b, BoolOp::Union).expect("the route applies").unwrap().unwrap(); + let v = measure(&joined).unwrap().volume; + assert!((v - expected).abs() < expected * 1e-4, "{v} against {expected}"); + assert_eq!(joined.lumps().unwrap().len(), 1); + assert!(joined.validity().is_empty(), "{:?}", joined.validity()); + // The flat faces are one plane again: nothing of the join stays below z = 0. + let mesh = joined.tessellate(0.05).unwrap(); + let lowest = mesh.positions.iter().map(|p| p.z).fold(f64::MAX, f64::min); + assert!(lowest > -1e-6, "the join reaches down to z = {lowest}"); + } +} + +/// Growing the flat face of a dome past its plane and cutting what that adds off again gives the dome back: the same volume, the same +/// faces, nothing below the plane. +#[test] +fn what_growing_a_face_adds_is_cut_off_again_exactly() { + let dome = hemisphere(20.0); + let measured = measure(&dome).unwrap().volume; + let grown = crate::coplanar::grow_faces_exact(&dome, -Vec3::Z, 0.0, 3.0).unwrap(); + assert!(measure(&grown).unwrap().volume > measured + 3.0 * 300.0, "the stub adds volume"); + let back = crate::coplanar::trim_exact(&grown, -Vec3::Z, 0.0).unwrap(); + let v = measure(&back).unwrap().volume; + assert!((v - measured).abs() < measured * 1e-6, "{v} against {measured}"); + assert_eq!(back.face_count(), dome.face_count()); + assert!(back.validity().is_empty(), "{:?}", back.validity()); + let lowest = back.tessellate(0.05).unwrap().positions.iter().map(|p| p.z).fold(f64::MAX, f64::min); + assert!(lowest > -1e-6, "the body reaches down to z = {lowest}"); +} + +/// A plane the body stands on has nothing beyond it, and one clear of the body is refused (nothing is cut off). +#[test] +fn trimming_at_a_plane_with_nothing_beyond_it() { + let dome = hemisphere(20.0); + let measured = measure(&dome).unwrap().volume; + let same = crate::coplanar::trim_exact(&dome, -Vec3::Z, 0.0).unwrap(); + assert!((measure(&same).unwrap().volume - measured).abs() < measured * 1e-9); + assert!(crate::coplanar::trim_exact(&dome, -Vec3::Z, 5.0).is_err()); +} diff --git a/crates/kernel/src/ops.rs b/crates/kernel/src/ops.rs index 0afa9ee2..491890e5 100644 --- a/crates/kernel/src/ops.rs +++ b/crates/kernel/src/ops.rs @@ -977,6 +977,30 @@ fn boolean_whole(a: &Body, b: &Body, op: BoolOp) -> Result> { Err(e) => last = format!("{last}; nudged: {e}"), } } + // A join whose coincident planar faces cannot be pushed apart because the neighbours of one of them curve (a flange + // flush with the end of a lofted tube): the faces are grown instead, and what the growth leaves beyond the plane is cut + // off by editing the faces. After every route above, the nudge included, has failed, so that a join any of them takes + // stays as it was. + if op == BoolOp::Union && !APART.with(|c| c.get()) { + APART.with(|c| c.set(true)); + let (r, diagnostics) = capture_boolean_diagnostics(|| crate::coplanar::boolean_grown(a, b, op)); + APART.with(|c| c.set(false)); + match r { + Some(Ok(Some(body))) => { + let (ok, v) = check(&body); + if ok { + merge(diagnostics); + let mut event = BooleanFallback::new(BooleanStrategy::CoplanarSeparation); + event.face_shift_mm = size * 0.0137; + record(event); + return Ok(Some(body)); + } + last = format!("{last}; faces grown apart gave an implausible volume {v:.4}"); + } + Some(Err(e)) => last = format!("{last}; {e}"), + Some(Ok(None)) | None => {} + } + } // The tries at the meeting region's scale come last, after every strategy that answered // before them, so a result any of those gave is kept as it was. if let Some(body) = fine(&mut last) { diff --git a/crates/kernel/src/planar_faces.rs b/crates/kernel/src/planar_faces.rs index d4df8fa6..17e1b313 100644 --- a/crates/kernel/src/planar_faces.rs +++ b/crates/kernel/src/planar_faces.rs @@ -36,6 +36,13 @@ fn control_points(s: &mt::Surface) -> Option> { } } +/// The least and the greatest `along · p` over the control points of a B-spline or NURBS surface; the surface lies +/// between them (it is in the convex hull of its control net). `None` for any other surface. +pub(crate) fn control_net_span(s: &mt::Surface, along: Vec3) -> Option<(f64, f64)> { + let pts = control_points(s)?; + pts.iter().map(|p| along.dot(*p)).fold(None, |acc, h| Some(acc.map_or((h, h), |(lo, hi): (f64, f64)| (lo.min(h), hi.max(h))))) +} + /// A point and the unit normal of the plane every point lies within `tol` of, if there is one. fn common_plane(pts: &[Vec3], tol: f64) -> Option<(Vec3, Vec3)> { let p0 = *pts.first()?;