From d585afd006a3cccca3b363a190c2c2f884242d65 Mon Sep 17 00:00:00 2001 From: bherbruck Date: Wed, 7 Oct 2026 19:22:26 +0000 Subject: [PATCH] M2: exact planar boolean fallback for coincident faces (BSP polygons rebuilt into B-rep) Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01AS5Tn3xieewUSg3uxrMVdL --- .gitignore | 2 +- crates/doc/src/eval.rs | 13 +- crates/engine/src/cmd/features.rs | 13 +- crates/kernel/src/build.rs | 38 ++- crates/kernel/src/lib.rs | 2 + crates/kernel/src/ops.rs | 13 + crates/kernel/src/polybool.rs | 457 ++++++++++++++++++++++++++++++ crates/kernel/src/tests.rs | 32 ++- 8 files changed, 552 insertions(+), 18 deletions(-) create mode 100644 crates/kernel/src/polybool.rs diff --git a/.gitignore b/.gitignore index 6ddd1c8b..af430a81 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,4 @@ -/target +/target*/ /plan/ /dist/ *~ diff --git a/crates/doc/src/eval.rs b/crates/doc/src/eval.rs index edfe6b55..3df68cf4 100644 --- a/crates/doc/src/eval.rs +++ b/crates/doc/src/eval.rs @@ -616,7 +616,9 @@ fn feature_tools(vals: &BTreeMap, f: &Feature, st: &ModelState) - FeatureKind::Loft { sections, .. } => { let mut secs = Vec::new(); for sec in sections { - let ss = st.sketch(sec.sketch).ok_or_else(|| DocError::Unknown(format!("sketch {} (it must come earlier in the timeline)", sec.sketch)))?; + let ss = st + .sketch(sec.sketch) + .ok_or_else(|| DocError::Unknown(format!("sketch {} (it must come earlier in the timeline)", sec.sketch)))?; let regions = solvecraft_sketch::merge_regions(&select_profiles(ss, &sec.profiles)?); let [r] = ®ions[..] else { return Err(DocError::Invalid("each loft section must be one profile".into())) }; secs.push((ss.plane, r.outer.clone())); @@ -625,7 +627,8 @@ fn feature_tools(vals: &BTreeMap, f: &Feature, st: &ModelState) - } FeatureKind::Sweep { sketch, profiles, path_sketch, path, .. } => { let ss = st.sketch(*sketch).ok_or_else(|| DocError::Unknown(format!("sketch {sketch} (it must come earlier in the timeline)")))?; - let ps = st.sketch(*path_sketch).ok_or_else(|| DocError::Unknown(format!("sketch {path_sketch} (it must come earlier in the timeline)")))?; + let ps = + st.sketch(*path_sketch).ok_or_else(|| DocError::Unknown(format!("sketch {path_sketch} (it must come earlier in the timeline)")))?; let regions = solvecraft_sketch::merge_regions(&select_profiles(ss, profiles)?); let start = regions.first().map(|r| ss.plane.to_world(r.centroid())).unwrap_or_default(); let segs = path_segments(ps, path, start)?; @@ -729,7 +732,7 @@ fn path_segments(ps: &SolvedSketch, ids: &[String], start: Vec3) -> Result 0) + if let (Some(f), true) = (segs.first().copied(), !segs.is_empty()) && w(f.end()).dist(start) < w(f.start()).dist(start) && (segs.len() == 1 || segs.get(1).is_some_and(|n| n.start().dist(f.start()) < 1e-6 || n.end().dist(f.start()) < 1e-6)) { @@ -749,7 +752,9 @@ fn path_segments(ps: &SolvedSketch, ids: &[String], start: Vec3) -> Result kernel::PathSeg::Line { a: w(a), b: w(b) }, - solvecraft_geom::Seg2::Arc { center, sweep, .. } => kernel::PathSeg::Arc { a: w(s.start()), center: w(center), axis: n * sweep.signum(), angle: sweep.abs() }, + solvecraft_geom::Seg2::Arc { center, sweep, .. } => { + kernel::PathSeg::Arc { a: w(s.start()), center: w(center), axis: n * sweep.signum(), angle: sweep.abs() } + } }) .collect()) } diff --git a/crates/engine/src/cmd/features.rs b/crates/engine/src/cmd/features.rs index fe586c2b..76d31716 100644 --- a/crates/engine/src/cmd/features.rs +++ b/crates/engine/src/cmd/features.rs @@ -597,7 +597,18 @@ fn sweep(s: &mut Session, p: &Value) -> Result { if path.is_empty() { return Err(bad(cmd, "`path` must list curve ids")); } - add_feature(s, p, FeatureKind::Sweep { sketch, profiles: profiles(p, cmd)?, path_sketch, path, operation: operation(p, cmd)?, targets: string_list(p, "targets") }) + add_feature( + s, + p, + FeatureKind::Sweep { + sketch, + profiles: profiles(p, cmd)?, + path_sketch, + path, + operation: operation(p, cmd)?, + targets: string_list(p, "targets"), + }, + ) } fn loft(s: &mut Session, p: &Value) -> Result { diff --git a/crates/kernel/src/build.rs b/crates/kernel/src/build.rs index 0a44a2b1..224956fe 100644 --- a/crates/kernel/src/build.rs +++ b/crates/kernel/src/build.rs @@ -398,7 +398,8 @@ fn revolve_touching_axis(plane: &Plane, lp: &Loop2, axis_origin: Vec2, d: Vec2) /// Close a ruled shell (side faces between matching wires) with planar caps into a body. fn cap_and_close(mut faces: Vec, w0: &[mt::Wire], w1: &[mt::Wire]) -> Result { - let bottom = builder::try_attach_plane(&w0.iter().map(|w| w.inverse()).collect::>()).map_err(|e| KernelError::Failed(format!("cap: {e}")))?; + let bottom = + builder::try_attach_plane(&w0.iter().map(|w| w.inverse()).collect::>()).map_err(|e| KernelError::Failed(format!("cap: {e}")))?; let top = builder::try_attach_plane(w1).map_err(|e| KernelError::Failed(format!("cap: {e}")))?; let mut last = String::new(); for (b, t) in [(bottom.clone(), top.clone()), (bottom.inverse(), top.inverse()), (bottom.clone(), top.inverse()), (bottom.inverse(), top)] { @@ -450,9 +451,10 @@ fn resample(lp: &Loop2, n: usize, angles: &[f64]) -> Loop2 { let m = a.lerp(b, 0.5); (Seg2::Line { a, b: m }, Seg2::Line { a: m, b }) } - Seg2::Arc { center, radius, start, sweep } => { - (Seg2::Arc { center, radius, start, sweep: sweep / 2.0 }, Seg2::Arc { center, radius, start: start + sweep / 2.0, sweep: sweep / 2.0 }) - } + Seg2::Arc { center, radius, start, sweep } => ( + Seg2::Arc { center, radius, start, sweep: sweep / 2.0 }, + Seg2::Arc { center, radius, start: start + sweep / 2.0, sweep: sweep / 2.0 }, + ), }; segs.insert(k, q); segs.insert(k, p); @@ -519,9 +521,17 @@ pub fn loft(sections: &[(Plane, Loop2)]) -> Result { /// A 3D path segment for sweeps. #[derive(Clone, Copy, Debug)] pub enum PathSeg { - Line { a: Vec3, b: Vec3 }, + Line { + a: Vec3, + b: Vec3, + }, /// Arc around `center`, rotating by `angle` about `axis` (right hand) from `a`. - Arc { a: Vec3, center: Vec3, axis: Vec3, angle: f64 }, + Arc { + a: Vec3, + center: Vec3, + axis: Vec3, + angle: f64, + }, } /// Sweep a planar region (outer loop) along a chain of path segments, keeping the profile @@ -543,7 +553,11 @@ pub fn sweep(plane: &Plane, region: &Region2, path: &[PathSeg]) -> Result builder::rsweep(&cur, p3(center), v3(ax), mt::Rad(angle)) } }; - let next = part.extract_boundaries().into_iter().find(|w| w.edge_iter().all(|e| !cur.edge_iter().any(|c| c.id() == e.id()))).ok_or_else(|| KernelError::Failed("sweep: lost the profile".into()))?; + let next = part + .extract_boundaries() + .into_iter() + .find(|w| w.edge_iter().all(|e| !cur.edge_iter().any(|c| c.id() == e.id()))) + .ok_or_else(|| KernelError::Failed("sweep: lost the profile".into()))?; shell.append(&mut part); cur = next.inverse(); } @@ -567,7 +581,9 @@ fn cubic_basis(knots: &[f64], n_ctrl: usize, u: f64) -> Vec { } for k in 1..=p { for j in 0..m.saturating_sub(1 + k) { - let (Some(&tj), Some(&tjk), Some(&tj1), Some(&tjk1)) = (knots.get(j), knots.get(j + k), knots.get(j + 1), knots.get(j + k + 1)) else { continue }; + let (Some(&tj), Some(&tjk), Some(&tj1), Some(&tjk1)) = (knots.get(j), knots.get(j + k), knots.get(j + 1), knots.get(j + k + 1)) else { + continue; + }; let left = if tjk - tj > 0.0 { (u - tj) / (tjk - tj) * n.get(j).copied().unwrap_or(0.0) } else { 0.0 }; let right = if tjk1 - tj1 > 0.0 { (tjk1 - u) / (tjk1 - tj1) * n.get(j + 1).copied().unwrap_or(0.0) } else { 0.0 }; if let Some(v) = n.get_mut(j) { @@ -609,9 +625,9 @@ fn interpolate_cubic(pts: &[Vec3]) -> Option> { if f == 0.0 { continue; } - for c in col..m { - let v = a[col][c]; - a[r][c] -= f * v; + let prow = a[col].clone(); + for (x, v) in a[r].iter_mut().zip(&prow).skip(col) { + *x -= f * v; } let qc = q[col]; for k in 0..3 { diff --git a/crates/kernel/src/lib.rs b/crates/kernel/src/lib.rs index e86b21a2..ba76df67 100644 --- a/crates/kernel/src/lib.rs +++ b/crates/kernel/src/lib.rs @@ -17,6 +17,7 @@ mod body; mod build; mod measure; mod ops; +mod polybool; mod polyhedron; mod step; mod topo; @@ -26,6 +27,7 @@ pub use body::{Body, EdgeInfo, FaceInfo}; pub use build::{PathSeg, box_solid, cylinder, extrude, extrude_tapered, loft, revolve, sphere, sweep, torus}; pub use measure::{BodyMeasure, measure}; pub use ops::{BoolOp, boolean, split_by_plane, transform, transform_matrix}; +pub use polybool::planar_boolean; pub use polyhedron::{HalfSpace, convex_polyhedron, draft, shell}; pub use step::step_export; pub use topo::{TopoCounts, merged_topology, seam_flags}; diff --git a/crates/kernel/src/ops.rs b/crates/kernel/src/ops.rs index fa2a698c..37704a05 100644 --- a/crates/kernel/src/ops.rs +++ b/crates/kernel/src/ops.rs @@ -136,6 +136,19 @@ pub fn boolean(a: &Body, b: &Body, op: BoolOp) -> Result> { return Ok(None); } } + // Coincident faces defeat the B-rep boolean; bodies with only planar faces have an exact + // polygon fallback. + match crate::polybool::planar_boolean(a, b, op) { + Ok(Some(body)) => { + let v = volume(&body); + if v > 0.0 && plausible(v, &body) { + return Ok(Some(body)); + } + last = format!("planar fallback gave an implausible volume {v:.4}"); + } + Ok(None) => return Ok(None), + Err(e) => last = format!("{last}; {e}"), + } Err(KernelError::Failed(format!("boolean {op:?}: {last}"))) } diff --git a/crates/kernel/src/polybool.rs b/crates/kernel/src/polybool.rs new file mode 100644 index 00000000..1f539b51 --- /dev/null +++ b/crates/kernel/src/polybool.rs @@ -0,0 +1,457 @@ +//! Exact-plane booleans for bodies whose faces are all planar (BSP-tree CSG on polygons, in +//! the classic style), used when the B-rep boolean fails — typically on coincident faces, such +//! as a wall flush with the side of its base. The polygons are rebuilt into a B-rep: coplanar +//! pieces are merged into faces, T-junctions are split so neighbouring faces share edges. + +use std::collections::HashMap; + +use solvecraft_geom::Vec3; +use truck_modeling::{self as mt, builder}; + +use crate::body::{Body, Solid, p3}; +use crate::{BoolOp, KernelError, Result, guard}; + +const MAX_DEPTH: usize = 4000; +const MAX_POLYS: usize = 200_000; + +#[derive(Clone, Copy, Debug)] +struct PlaneEq { + n: Vec3, + w: f64, +} + +impl PlaneEq { + fn flipped(self) -> PlaneEq { + PlaneEq { n: -self.n, w: -self.w } + } +} + +#[derive(Clone, Debug)] +struct Poly { + v: Vec, + plane: PlaneEq, +} + +impl Poly { + fn flip(&mut self) { + self.v.reverse(); + self.plane = self.plane.flipped(); + } +} + +const COPLANAR: u8 = 0; +const FRONT: u8 = 1; +const BACK: u8 = 2; +const SPANNING: u8 = 3; + +/// Split `p` by `pl` into the four buckets. +fn split(pl: PlaneEq, p: &Poly, eps: f64, cf: &mut Vec, cb: &mut Vec, f: &mut Vec, b: &mut Vec) { + let mut ptype = 0u8; + let types: Vec = + p.v.iter() + .map(|v| { + let t = pl.n.dot(*v) - pl.w; + let ty = if t < -eps { + BACK + } else if t > eps { + FRONT + } else { + COPLANAR + }; + ptype |= ty; + ty + }) + .collect(); + match ptype { + COPLANAR => { + if pl.n.dot(p.plane.n) > 0.0 { + cf.push(p.clone()); + } else { + cb.push(p.clone()); + } + } + FRONT => f.push(p.clone()), + BACK => b.push(p.clone()), + _ => { + let (mut fv, mut bv) = (Vec::new(), Vec::new()); + let n = p.v.len(); + for i in 0..n { + let j = (i + 1) % n; + let (Some(&ti), Some(&tj), Some(&vi), Some(&vj)) = (types.get(i), types.get(j), p.v.get(i), p.v.get(j)) else { continue }; + if ti != BACK { + fv.push(vi); + } + if ti != FRONT { + bv.push(vi); + } + if (ti | tj) == SPANNING { + let den = pl.n.dot(vj - vi); + if den.abs() > 1e-300 { + let t = (pl.w - pl.n.dot(vi)) / den; + let x = vi.lerp(vj, t); + fv.push(x); + bv.push(x); + } + } + } + if fv.len() >= 3 { + f.push(Poly { v: fv, plane: p.plane }); + } + if bv.len() >= 3 { + b.push(Poly { v: bv, plane: p.plane }); + } + } + } +} + +#[derive(Default)] +struct Node { + plane: Option, + front: Option>, + back: Option>, + polys: Vec, +} + +impl Node { + fn build(&mut self, polys: Vec, eps: f64, depth: usize) -> Result<()> { + if polys.is_empty() { + return Ok(()); + } + if depth > MAX_DEPTH { + return Err(KernelError::Failed("planar boolean: tree too deep".into())); + } + let pl = match self.plane { + Some(p) => p, + None => { + let p = polys.first().map(|p| p.plane).ok_or_else(|| KernelError::Failed("bsp".into()))?; + self.plane = Some(p); + p + } + }; + let (mut f, mut b) = (Vec::new(), Vec::new()); + let mut cf = Vec::new(); + let mut cb = Vec::new(); + for p in &polys { + split(pl, p, eps, &mut cf, &mut cb, &mut f, &mut b); + } + self.polys.extend(cf); + self.polys.extend(cb); + if !f.is_empty() { + self.front.get_or_insert_with(Default::default).build(f, eps, depth + 1)?; + } + if !b.is_empty() { + self.back.get_or_insert_with(Default::default).build(b, eps, depth + 1)?; + } + Ok(()) + } + fn invert(&mut self) { + for p in &mut self.polys { + p.flip(); + } + self.plane = self.plane.map(PlaneEq::flipped); + if let Some(f) = &mut self.front { + f.invert(); + } + if let Some(b) = &mut self.back { + b.invert(); + } + std::mem::swap(&mut self.front, &mut self.back); + } + fn clip_polys(&self, polys: Vec, eps: f64) -> Vec { + let Some(pl) = self.plane else { return polys }; + let (mut f, mut b) = (Vec::new(), Vec::new()); + let (mut cf, mut cb) = (Vec::new(), Vec::new()); + for p in &polys { + split(pl, p, eps, &mut cf, &mut cb, &mut f, &mut b); + } + f.extend(cf); + b.extend(cb); + let f = match &self.front { + Some(n) => n.clip_polys(f, eps), + None => f, + }; + let b = match &self.back { + Some(n) => n.clip_polys(b, eps), + None => Vec::new(), + }; + let mut out = f; + out.extend(b); + out + } + fn clip_to(&mut self, other: &Node, eps: f64) { + let polys = std::mem::take(&mut self.polys); + self.polys = other.clip_polys(polys, eps); + if let Some(f) = &mut self.front { + f.clip_to(other, eps); + } + if let Some(b) = &mut self.back { + b.clip_to(other, eps); + } + } + fn all(&self, out: &mut Vec) { + out.extend(self.polys.iter().cloned()); + if let Some(f) = &self.front { + f.all(out); + } + if let Some(b) = &self.back { + b.all(out); + } + } +} + +/// Polygons of a body whose faces are all planar (None otherwise). +fn planar_polys(b: &Body) -> Option> { + let tol = (b.size() * 1e-3).max(1e-3); + let faces = b.faces(tol).ok()?; + let mesh = b.tessellate(tol).ok()?; + let planes: Vec = faces.iter().map(|f| f.plane_normal.map(|n| PlaneEq { n, w: n.dot(f.centroid) })).collect::>()?; + let mut out = Vec::new(); + for (t, fi) in mesh.triangles.iter().zip(&mesh.tri_face) { + let [a, bb, c] = mesh.tri(t)?; + let plane = *planes.get(*fi as usize)?; + // Snap the vertices onto the exact face plane. + let snap = |p: Vec3| p - plane.n * (plane.n.dot(p) - plane.w); + let mut v = vec![snap(a), snap(bb), snap(c)]; + if (v[1] - v[0]).cross(v[2] - v[0]).dot(plane.n) < 0.0 { + v.reverse(); + } + out.push(Poly { v, plane }); + } + (out.len() <= MAX_POLYS).then_some(out) +} + +/// Boolean of two all-planar bodies. +pub fn planar_boolean(a: &Body, b: &Body, op: BoolOp) -> Result> { + let (Some(pa), Some(pb)) = (planar_polys(a), planar_polys(b)) else { + return Err(KernelError::Failed("not supported yet: coincident faces on curved bodies".into())); + }; + let size = a.size().max(b.size()); + let eps = size * 1e-9; + let (mut na, mut nb) = (Node::default(), Node::default()); + na.build(pa, eps, 0)?; + nb.build(pb, eps, 0)?; + match op { + BoolOp::Union => { + na.clip_to(&nb, eps); + nb.clip_to(&na, eps); + nb.invert(); + nb.clip_to(&na, eps); + nb.invert(); + let mut extra = Vec::new(); + nb.all(&mut extra); + na.build(extra, eps, 0)?; + } + BoolOp::Cut => { + na.invert(); + na.clip_to(&nb, eps); + nb.clip_to(&na, eps); + nb.invert(); + nb.clip_to(&na, eps); + nb.invert(); + let mut extra = Vec::new(); + nb.all(&mut extra); + na.build(extra, eps, 0)?; + na.invert(); + } + BoolOp::Intersect => { + na.invert(); + nb.clip_to(&na, eps); + nb.invert(); + na.clip_to(&nb, eps); + nb.clip_to(&na, eps); + let mut extra = Vec::new(); + nb.all(&mut extra); + na.build(extra, eps, 0)?; + na.invert(); + } + } + let mut polys = Vec::new(); + na.all(&mut polys); + if polys.is_empty() { + return Ok(None); + } + rebuild(polys, size).map(Some) +} + +/// Position key for vertex merging. +fn key(p: Vec3, q: f64) -> (i64, i64, i64) { + ((p.x / q).round() as i64, (p.y / q).round() as i64, (p.z / q).round() as i64) +} + +/// Turn polygons into a B-rep solid: merge by plane, split T-junctions, chain boundary loops. +fn rebuild(polys: Vec, size: f64) -> Result { + let q = size * 1e-7; + // Global vertex table. + let mut index: HashMap<(i64, i64, i64), usize> = HashMap::new(); + let mut verts: Vec = Vec::new(); + let mut id = |p: Vec3, verts: &mut Vec| -> usize { + let k = key(p, q); + *index.entry(k).or_insert_with(|| { + verts.push(p); + verts.len() - 1 + }) + }; + let mut loops: Vec<(PlaneEq, Vec)> = Vec::new(); + for p in &polys { + let mut l: Vec = p.v.iter().map(|v| id(*v, &mut verts)).collect(); + l.dedup(); + if l.first() == l.last() && l.len() > 1 { + l.pop(); + } + if l.len() >= 3 { + loops.push((p.plane, l)); + } + } + // Split every polygon edge at vertices lying on it (T-junctions). + let tol = size * 1e-7; + let mut directed: Vec<(usize, usize, usize)> = Vec::new(); // (plane group, from, to) + let mut groups: Vec = Vec::new(); + for (pl, l) in &loops { + let g = match groups.iter().position(|g| g.n.dot(pl.n) > 1.0 - 1e-9 && (g.w - pl.w).abs() < tol * 10.0) { + Some(g) => g, + None => { + groups.push(*pl); + groups.len() - 1 + } + }; + let n = l.len(); + for i in 0..n { + let (Some(&a), Some(&b)) = (l.get(i), l.get((i + 1) % n)) else { continue }; + let (Some(&pa), Some(&pb)) = (verts.get(a), verts.get(b)) else { continue }; + let d = pb - pa; + let len2 = d.len2(); + let mut on: Vec<(f64, usize)> = verts + .iter() + .enumerate() + .filter(|(k, _)| *k != a && *k != b) + .filter_map(|(k, v)| { + let t = (*v - pa).dot(d) / len2.max(1e-300); + (t > 1e-9 && t < 1.0 - 1e-9 && v.dist_to_segment(pa, pb) < tol).then_some((t, k)) + }) + .collect(); + on.sort_by(|x, y| x.0.total_cmp(&y.0)); + let mut prev = a; + for (_, k) in on { + directed.push((g, prev, k)); + prev = k; + } + directed.push((g, prev, b)); + } + } + // Interior edges cancel against their reverse within the same plane group. + let mut count: HashMap<(usize, usize, usize), i64> = HashMap::new(); + for (g, a, b) in &directed { + *count.entry((*g, *a, *b)).or_insert(0) += 1; + } + let mut boundary: Vec<(usize, usize, usize)> = Vec::new(); + for (&(g, a, b), &c) in &count { + let rev = count.get(&(g, b, a)).copied().unwrap_or(0); + for _ in 0..(c - rev).max(0) { + boundary.push((g, a, b)); + } + } + boundary.sort_unstable(); + guard("planar boolean rebuild", || { + let tv: Vec = verts.iter().map(|v| builder::vertex(p3(*v))).collect(); + let mut edges: HashMap<(usize, usize), mt::Edge> = HashMap::new(); + let mut faces: Vec = Vec::new(); + for (g, plane) in groups.iter().enumerate() { + let mut segs: Vec<(usize, usize)> = boundary.iter().filter(|e| e.0 == g).map(|e| (e.1, e.2)).collect(); + let mut wires: Vec = Vec::new(); + let mut loops_g: Vec> = Vec::new(); + while let Some((s, t)) = segs.pop() { + let mut lp = vec![s, t]; + let mut cur = t; + let mut guard_n = 0; + while cur != s && guard_n < 100_000 { + guard_n += 1; + let Some(k) = segs.iter().position(|(a, _)| *a == cur) else { break }; + let (_, nb) = segs.remove(k); + lp.push(nb); + cur = nb; + } + if cur == s { + lp.pop(); + loops_g.push(lp); + } + } + for lp in &loops_g { + let n = lp.len(); + let mut es = Vec::with_capacity(n); + for i in 0..n { + let (Some(&a), Some(&b)) = (lp.get(i), lp.get((i + 1) % n)) else { continue }; + let k = (a.min(b), a.max(b)); + let e = match edges.get(&k) { + Some(e) => e.clone(), + None => { + let (Some(va), Some(vb)) = (tv.get(k.0), tv.get(k.1)) else { continue }; + let e = builder::line(va, vb); + edges.insert(k, e.clone()); + e + } + }; + es.push(if a < b { e } else { e.inverse() }); + } + wires.push(es.into()); + } + if wires.is_empty() { + continue; + } + // Outer loops are counter-clockwise about the plane normal; holes go with the + // outer loop that contains them. + let area = |lp: &Vec| -> f64 { + let n = lp.len(); + let mut s = Vec3::ZERO; + for i in 0..n { + if let (Some(a), Some(b)) = (lp.get(i).and_then(|k| verts.get(*k)), lp.get((i + 1) % n).and_then(|k| verts.get(*k))) { + s += a.cross(*b); + } + } + s.dot(plane.n) * 0.5 + }; + let (outer, holes): (Vec, Vec) = (0..loops_g.len()).partition(|i| loops_g.get(*i).is_some_and(|l| area(l) > 0.0)); + let u = plane.n.any_perp(); + let w = plane.n.cross(u); + let to2 = |p: Vec3| (p.dot(u), p.dot(w)); + let inside = |pt: Vec3, lp: &Vec| { + let (px, py) = to2(pt); + let n = lp.len(); + let mut c = false; + for i in 0..n { + let (Some(a), Some(b)) = (lp.get(i).and_then(|k| verts.get(*k)), lp.get((i + 1) % n).and_then(|k| verts.get(*k))) else { + continue; + }; + let ((ax, ay), (bx, by)) = (to2(*a), to2(*b)); + if (ay > py) != (by > py) && px < ax + (py - ay) / (by - ay) * (bx - ax) { + c = !c; + } + } + c + }; + for o in &outer { + let Some(ol) = loops_g.get(*o) else { continue }; + let mut ws = vec![wires.get(*o).cloned().ok_or_else(|| KernelError::Failed("wire".into()))?]; + for h in &holes { + let Some(hl) = loops_g.get(*h) else { continue }; + let probe = hl.first().and_then(|k| verts.get(*k)).copied().unwrap_or_default(); + // The innermost outer loop containing the hole. + let owner = outer.iter().filter(|x| loops_g.get(**x).is_some_and(|l| inside(probe, l))).min_by(|x, y| { + let ax = loops_g.get(**x).map(area).unwrap_or(0.0); + let ay = loops_g.get(**y).map(area).unwrap_or(0.0); + ax.total_cmp(&ay) + }); + if owner == Some(o) { + ws.push(wires.get(*h).cloned().ok_or_else(|| KernelError::Failed("wire".into()))?); + } + } + let _ = ol; + let face = builder::try_attach_plane(&ws).map_err(|e| KernelError::Failed(format!("planar boolean face: {e}")))?; + faces.push(face); + } + } + let shell: mt::Shell = faces.into(); + let parts = shell.connected_components(); + let solid = Solid::try_new(parts).map_err(|e| KernelError::Failed(format!("planar boolean: {e}")))?; + Body::new(solid) + }) +} diff --git a/crates/kernel/src/tests.rs b/crates/kernel/src/tests.rs index f2b1efd5..9a09a279 100644 --- a/crates/kernel/src/tests.rs +++ b/crates/kernel/src/tests.rs @@ -204,7 +204,12 @@ fn loft_and_sweep() { let disc = Region2 { outer: Loop2::circle(Vec2::ZERO, 3.0), holes: vec![] }; let path = [ PathSeg::Line { a: Vec3::ZERO, b: Vec3::new(0.0, 0.0, 30.0) }, - PathSeg::Arc { a: Vec3::new(0.0, 0.0, 30.0), center: Vec3::new(20.0, 0.0, 30.0), axis: Vec3::new(0.0, 1.0, 0.0), angle: std::f64::consts::FRAC_PI_2 }, + PathSeg::Arc { + a: Vec3::new(0.0, 0.0, 30.0), + center: Vec3::new(20.0, 0.0, 30.0), + axis: Vec3::new(0.0, 1.0, 0.0), + angle: std::f64::consts::FRAC_PI_2, + }, PathSeg::Line { a: Vec3::new(20.0, 0.0, 50.0), b: Vec3::new(50.0, 0.0, 50.0) }, ]; let t = sweep(&Plane::XY, &disc, &path).unwrap(); @@ -216,6 +221,31 @@ fn loft_and_sweep() { assert!(rel(measure(&t).unwrap().volume, expect) < 1e-3, "{} vs {expect}", measure(&t).unwrap().volume); } +#[test] +fn planar_booleans_with_coincident_faces() { + // L-bracket: a wall flush with three sides of its base. + let base = box_solid(Vec3::ZERO, Vec3::new(80.0, 50.0, 8.0)).unwrap(); + let wall = box_solid(Vec3::new(0.0, 42.0, 8.0), Vec3::new(80.0, 50.0, 50.0)).unwrap(); + let u = planar_boolean(&base, &wall, BoolOp::Union).unwrap().unwrap(); + let m = measure(&u).unwrap(); + assert!(rel(m.volume, 80.0 * 50.0 * 8.0 + 80.0 * 8.0 * 42.0) < 1e-9, "{}", m.volume); + assert_eq!(m.merged.faces, 8, "{:?}", m.merged); + // Through the public boolean (B-rep first, polygon fallback). + let u2 = boolean(&base, &wall, BoolOp::Union).unwrap().unwrap(); + assert!(rel(measure(&u2).unwrap().volume, m.volume) < 1e-9); + // Two boxes sharing a face, and a cut flush with a side. + let a = box_solid(Vec3::ZERO, Vec3::new(10.0, 10.0, 10.0)).unwrap(); + let b = box_solid(Vec3::new(10.0, 0.0, 0.0), Vec3::new(20.0, 10.0, 10.0)).unwrap(); + let ab = boolean(&a, &b, BoolOp::Union).unwrap().unwrap(); + assert!(rel(measure(&ab).unwrap().volume, 2000.0) < 1e-9); + assert_eq!(measure(&ab).unwrap().merged.faces, 6); + let notch = box_solid(Vec3::new(0.0, 0.0, 5.0), Vec3::new(5.0, 10.0, 10.0)).unwrap(); + let c = boolean(&a, ¬ch, BoolOp::Cut).unwrap().unwrap(); + assert!(rel(measure(&c).unwrap().volume, 1000.0 - 250.0) < 1e-9); + let i = boolean(&a, ¬ch, BoolOp::Intersect).unwrap().unwrap(); + assert!(rel(measure(&i).unwrap().volume, 250.0) < 1e-9); +} + #[test] #[ignore] fn debug_convergence() {