M2: exact planar boolean fallback for coincident faces (BSP polygons rebuilt into B-rep)

Co-Authored-By: Claude Opus 5.5 <[email protected]>
Claude-Session: https://claude.ai/code/session_01AS5Tn3xieewUSg3uxrMVdL

.gitignore

@@ -1,4 +1,4 @@
-/target
+/target*/
 /plan/
 /dist/
 *~

crates/doc/src/eval.rs

@@ -616,7 +616,9 @@ fn feature_tools(vals: &BTreeMap<String, Value>, 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] = &regions[..] 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<String, Value>, 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<Vec<k
     // Orient the chain: the first segment starts at the end nearest `start`, each next one
     // starts where the previous ended.
     let w = |p: Vec2| ps.plane.to_world(p);
-    if let (Some(f), true) = (segs.first().copied(), segs.len() > 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<Vec<k
         .iter()
         .map(|s| match *s {
             solvecraft_geom::Seg2::Line { a, b } => 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())
 }

crates/engine/src/cmd/features.rs

@@ -597,7 +597,18 @@ fn sweep(s: &mut Session, p: &Value) -> Result<Value> {
     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<Value> {

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<mt::Face>, w0: &[mt::Wire], w1: &[mt::Wire]) -> Result<Body> {
-    let bottom = builder::try_attach_plane(&w0.iter().map(|w| w.inverse()).collect::<Vec<_>>()).map_err(|e| KernelError::Failed(format!("cap: {e}")))?;
+    let bottom =
+        builder::try_attach_plane(&w0.iter().map(|w| w.inverse()).collect::<Vec<_>>()).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<Body> {
 /// 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<Body>
                     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<f64> {
     }
     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<mt::BSplineCurve<mt::Point3>> {
             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 {

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

crates/kernel/src/ops.rs

@@ -136,6 +136,19 @@ pub fn boolean(a: &Body, b: &Body, op: BoolOp) -> Result<Option<Body>> {
             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}")))
 }
 

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<Vec3>,
+    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<Poly>, cb: &mut Vec<Poly>, f: &mut Vec<Poly>, b: &mut Vec<Poly>) {
+    let mut ptype = 0u8;
+    let types: Vec<u8> =
+        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<PlaneEq>,
+    front: Option<Box<Node>>,
+    back: Option<Box<Node>>,
+    polys: Vec<Poly>,
+}
+
+impl Node {
+    fn build(&mut self, polys: Vec<Poly>, 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<Poly>, eps: f64) -> Vec<Poly> {
+        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<Poly>) {
+        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<Vec<Poly>> {
+    let tol = (b.size() * 1e-3).max(1e-3);
+    let faces = b.faces(tol).ok()?;
+    let mesh = b.tessellate(tol).ok()?;
+    let planes: Vec<PlaneEq> = faces.iter().map(|f| f.plane_normal.map(|n| PlaneEq { n, w: n.dot(f.centroid) })).collect::<Option<_>>()?;
+    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<Option<Body>> {
+    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<Poly>, size: f64) -> Result<Body> {
+    let q = size * 1e-7;
+    // Global vertex table.
+    let mut index: HashMap<(i64, i64, i64), usize> = HashMap::new();
+    let mut verts: Vec<Vec3> = Vec::new();
+    let mut id = |p: Vec3, verts: &mut Vec<Vec3>| -> usize {
+        let k = key(p, q);
+        *index.entry(k).or_insert_with(|| {
+            verts.push(p);
+            verts.len() - 1
+        })
+    };
+    let mut loops: Vec<(PlaneEq, Vec<usize>)> = Vec::new();
+    for p in &polys {
+        let mut l: Vec<usize> = 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<PlaneEq> = 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<mt::Vertex> = verts.iter().map(|v| builder::vertex(p3(*v))).collect();
+        let mut edges: HashMap<(usize, usize), mt::Edge> = HashMap::new();
+        let mut faces: Vec<mt::Face> = 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<mt::Wire> = Vec::new();
+            let mut loops_g: Vec<Vec<usize>> = 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<usize>| -> 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<usize>, Vec<usize>) = (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<usize>| {
+                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)
+    })
+}

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, &notch, BoolOp::Cut).unwrap().unwrap();
+    assert!(rel(measure(&c).unwrap().volume, 1000.0 - 250.0) < 1e-9);
+    let i = boolean(&a, &notch, BoolOp::Intersect).unwrap().unwrap();
+    assert!(rel(measure(&i).unwrap().volume, 250.0) < 1e-9);
+}
+
 #[test]
 #[ignore]
 fn debug_convergence() {