brepkit-operations 4.6.3

//! Operations on compound entities.
//!
//! Provides utilities for working with compounds of solids:
//! extracting individual solids, fusing all solids in a compound,
//! and computing compound-level measurements.

use brepkit_math::aabb::Aabb3;
use brepkit_topology::Topology;
use brepkit_topology::compound::CompoundId;
use brepkit_topology::solid::SolidId;

/// Extract all solid IDs from a compound.
///
/// # Errors
///
/// Returns an error if the compound ID is invalid.
pub fn explode(
    topo: &Topology,
    compound: CompoundId,
) -> Result<Vec<SolidId>, crate::OperationsError> {
    let comp = topo.compound(compound)?;
    Ok(comp.solids().to_vec())
}

/// Fuse (union) all solids in a compound into a single solid.
///
/// Performs iterative boolean union on all solids. Requires at least
/// one solid in the compound.
///
/// # Errors
///
/// Returns an error if the compound is empty or a boolean operation fails.
pub fn fuse_all(
    topo: &mut Topology,
    compound: CompoundId,
) -> Result<SolidId, crate::OperationsError> {
    let solids = {
        let comp = topo.compound(compound)?;
        comp.solids().to_vec()
    };

    if solids.is_empty() {
        return Err(crate::OperationsError::InvalidInput {
            reason: "compound has no solids to fuse".into(),
        });
    }

    // Partition solids into overlapping groups. Disjoint solids can be merged
    // directly (no boolean needed), while overlapping groups use boolean fuse.
    let bboxes: Vec<Aabb3> = solids
        .iter()
        .map(|&sid| crate::measure::solid_bounding_box(topo, sid))
        .collect::<Result<_, _>>()?;

    // Per-solid polyhedral bounds (plane normals + vertices), or `None` for any
    // solid with a curved face. Lets `partition_touching` prove that two solids
    // whose loose AABBs overlap are actually disjoint (e.g. honeycomb hex prisms
    // packed tighter than their corner-to-corner AABB extent), keeping them off
    // the expensive boolean path.
    let margin = brepkit_math::tolerance::Tolerance::new().linear;
    let poly_bounds: Vec<Option<PolyhedralBounds>> =
        solids.iter().map(|&s| polyhedral_bounds(topo, s)).collect();

    let groups = partition_touching(&bboxes, &poly_bounds, margin);

    let mut group_results: Vec<SolidId> = Vec::new();
    for group in &groups {
        let group_solids: Vec<SolidId> = group.iter().map(|&i| solids[i]).collect();
        if group_solids.len() == 1 {
            group_results.push(group_solids[0]);
            continue;
        }
        // Each group is a connected cluster of interpenetrating/touching solids
        // — fuse it in ONE GFA arrangement (via `fuse_cluster`, N-way with a
        // sequential fallback) instead of a pairwise reduction that re-processes
        // a growing accumulator O(n²).
        //
        // A 2-solid group takes the pairwise contract directly: fuse_cluster's
        // mesh-fallback bail exists to keep a degraded fuse from poisoning the
        // REST of a batch, but a pair has no rest — erroring only makes the
        // caller redo this exact fuse pairwise and pay the fallback twice
        // (brepjs's fuseAllBisect retry did exactly that on every degraded
        // pair fuse).
        group_results.push(match group_solids.as_slice() {
            &[a, b] => crate::boolean::boolean(topo, crate::boolean::BooleanOp::Fuse, a, b)?,
            _ => crate::boolean::fuse_cluster(topo, &group_solids)?,
        });
    }

    if group_results.len() == 1 {
        return Ok(group_results[0]);
    }

    merge_disjoint_solids(topo, &group_results)
}

/// Count the total number of solids in a compound.
///
/// # Errors
///
/// Returns an error if the compound ID is invalid.
pub fn solid_count(topo: &Topology, compound: CompoundId) -> Result<usize, crate::OperationsError> {
    let comp = topo.compound(compound)?;
    Ok(comp.solids().len())
}

/// Compute the combined bounding box of all solids in a compound.
///
/// # Errors
///
/// Returns an error if the compound is empty or measurement fails.
pub fn compound_bounding_box(
    topo: &Topology,
    compound: CompoundId,
) -> Result<brepkit_math::aabb::Aabb3, crate::OperationsError> {
    let comp = topo.compound(compound)?;
    let solids = comp.solids();

    if solids.is_empty() {
        return Err(crate::OperationsError::InvalidInput {
            reason: "compound is empty".into(),
        });
    }

    let mut combined = crate::measure::solid_bounding_box(topo, solids[0])?;
    for &sid in &solids[1..] {
        let bb = crate::measure::solid_bounding_box(topo, sid)?;
        combined = combined.union(bb);
    }

    Ok(combined)
}

/// Union-find path-compressed lookup.
fn uf_find(parent: &mut [usize], mut x: usize) -> usize {
    while parent[x] != x {
        parent[x] = parent[parent[x]];
        x = parent[x];
    }
    x
}

/// Plane normals (candidate separating axes) and boundary vertices of a
/// solid, used to prove disjointness via the separating-axis theorem.
struct PolyhedralBounds {
    normals: Vec<brepkit_math::vec::Vec3>,
    verts: Vec<brepkit_math::vec::Point3>,
}

/// Collect a solid's plane normals and vertices — but only if *every* outer-shell
/// face is planar. A flat-faced solid is contained in the convex hull of its
/// vertices, which makes the vertex-projection separation test (below) sound.
/// A single curved face can bulge past that hull, so any non-`Plane` face makes
/// this return `None` (the caller then falls back to the conservative AABB test).
fn polyhedral_bounds(topo: &Topology, sid: SolidId) -> Option<PolyhedralBounds> {
    use brepkit_topology::face::FaceSurface;

    let solid = topo.solid(sid).ok()?;
    let shell = topo.shell(solid.outer_shell()).ok()?;

    let mut normals = Vec::new();
    let mut vert_ids = std::collections::HashSet::new();
    for &fid in shell.faces() {
        let face = topo.face(fid).ok()?;
        match face.surface() {
            // Normalize: stored plane normals aren't guaranteed unit length (e.g.
            // raw STEP `DIRECTION` data), and `polyhedral_separated` compares
            // projection gaps against a world-space margin, which is only valid
            // for unit axes. Bail the whole solid to the AABB path on a
            // degenerate normal.
            FaceSurface::Plane { normal, .. } => normals.push(normal.normalize().ok()?),
            _ => return None,
        }
        for wid in std::iter::once(face.outer_wire()).chain(face.inner_wires().iter().copied()) {
            let wire = topo.wire(wid).ok()?;
            for oe in wire.edges() {
                let edge = topo.edge(oe.edge()).ok()?;
                vert_ids.insert(edge.start());
                vert_ids.insert(edge.end());
            }
        }
    }

    let mut verts = Vec::with_capacity(vert_ids.len());
    for vid in vert_ids {
        verts.push(topo.vertex(vid).ok()?.point());
    }
    if verts.is_empty() {
        return None;
    }
    Some(PolyhedralBounds { normals, verts })
}

/// Whether two flat-faced solids are provably disjoint: `true` iff some face
/// normal of either separates their vertex projections by a clear `margin`.
///
/// Soundness: each solid lies within the convex hull of its vertices (all faces
/// planar), so a gap between the vertex projections on any axis is a real gap
/// between the solids. Only face-normal axes are tried (not edge-edge cross
/// products), so the test is sound but not complete — an undetected separation
/// just falls through to the boolean, never a false "disjoint" for touching
/// inputs.
fn polyhedral_separated(a: &PolyhedralBounds, b: &PolyhedralBounds, margin: f64) -> bool {
    let project = |verts: &[brepkit_math::vec::Point3], axis: &brepkit_math::vec::Vec3| {
        let mut lo = f64::INFINITY;
        let mut hi = f64::NEG_INFINITY;
        for p in verts {
            let d = p.x() * axis.x() + p.y() * axis.y() + p.z() * axis.z();
            lo = lo.min(d);
            hi = hi.max(d);
        }
        (lo, hi)
    };
    a.normals.iter().chain(b.normals.iter()).any(|axis| {
        let (a_lo, a_hi) = project(&a.verts, axis);
        let (b_lo, b_hi) = project(&b.verts, axis);
        b_lo - a_hi > margin || a_lo - b_hi > margin
    })
}

/// Partition indices into groups that may actually touch (union-find).
///
/// Two solids share a group when their AABBs overlap or come within `margin`,
/// *unless* both are flat-faced and a separating axis proves a real gap
/// between them. This keeps geometrically
/// disjoint pieces whose loose AABBs overlap (honeycomb hex prisms, tightly
/// packed feet) in separate groups, so `fuse_all` merges them via the cheap
/// disjoint-shell path instead of an O(n) chain of boolean unions.
fn partition_touching(
    bboxes: &[Aabb3],
    poly_bounds: &[Option<PolyhedralBounds>],
    margin: f64,
) -> Vec<Vec<usize>> {
    let n = bboxes.len();
    let mut parent: Vec<usize> = (0..n).collect();

    for i in 0..n {
        for j in (i + 1)..n {
            // Boxes within the margin touch: a face two solids share can
            // stand a rounding error apart in each.
            if !bboxes[i].expanded(margin).intersects(bboxes[j]) {
                continue;
            }
            // AABBs overlap. Only keep them apart if we can *prove* a gap.
            if let (Some(pi), Some(pj)) = (&poly_bounds[i], &poly_bounds[j])
                && polyhedral_separated(pi, pj, margin)
            {
                continue;
            }
            let ri = uf_find(&mut parent, i);
            let rj = uf_find(&mut parent, j);
            if ri != rj {
                parent[ri] = rj;
            }
        }
    }

    let mut groups: std::collections::HashMap<usize, Vec<usize>> = std::collections::HashMap::new();
    for i in 0..n {
        groups.entry(uf_find(&mut parent, i)).or_default().push(i);
    }
    groups.into_values().collect()
}

/// Merge disjoint solids into a single solid by combining all faces.
///
/// Note: the resulting outer shell contains disconnected face groups,
/// which technically violates the connected-shell invariant. This is
/// acceptable for volume measurement and tessellation (which iterate
/// faces independently), but algorithms that assume shell connectivity
/// should be aware. A future improvement would return a `Compound`.
///
/// The result references the input solids' existing faces (no deep copy),
/// so callers that need an independent result must pass copies.
pub(crate) fn merge_disjoint_solids(
    topo: &mut Topology,
    solids: &[SolidId],
) -> Result<SolidId, crate::OperationsError> {
    use brepkit_topology::shell::Shell;
    use brepkit_topology::solid::Solid;

    let mut all_faces = Vec::new();
    let mut inner_shell_ids = Vec::new();

    // Snapshot phase: collect all face IDs and inner shell face sets.
    let mut inner_face_sets: Vec<Vec<brepkit_topology::face::FaceId>> = Vec::new();
    for &sid in solids {
        let solid_data = topo.solid(sid)?;
        let outer_shell = topo.shell(solid_data.outer_shell())?;
        all_faces.extend_from_slice(outer_shell.faces());

        let inner_ids: Vec<_> = solid_data.inner_shells().to_vec();
        for inner_id in inner_ids {
            let inner_shell = topo.shell(inner_id)?;
            inner_face_sets.push(inner_shell.faces().to_vec());
        }
    }

    // Allocate phase: create inner shells.
    for faces in inner_face_sets {
        let inner = Shell::new(faces).map_err(crate::OperationsError::Topology)?;
        inner_shell_ids.push(topo.add_shell(inner));
    }

    let outer = Shell::new(all_faces).map_err(crate::OperationsError::Topology)?;
    let outer_id = topo.add_shell(outer);
    Ok(topo.add_solid(Solid::new(outer_id, inner_shell_ids)))
}

/// Split a solid into one solid per separate piece of material.
///
/// A boolean that leaves separate pieces returns them in one solid whose
/// outer shell holds every piece's faces. Each set of faces joined through
/// shared edges bounds a region; a region inside an odd number of closed
/// regions is a cavity of the innermost one around it, and every other
/// region is the outside of a piece. A solid whose outer shell is one
/// connected set of faces comes back as is.
///
/// The pieces share the solid's faces rather than copying them;
/// an empty solid has none.
///
/// # Errors
///
/// Returns an error if the solid holds invalid topology references or a
/// point cannot be classified against one of its regions.
pub fn solid_pieces(
    topo: &mut Topology,
    solid: SolidId,
) -> Result<Vec<SolidId>, crate::OperationsError> {
    use brepkit_check::classify::{ClassifyOptions, PointClassification, SolidClassifier};
    use brepkit_topology::shell::Shell;
    use brepkit_topology::solid::Solid;

    let faces = brepkit_topology::explorer::solid_faces(topo, solid)?;
    if faces.is_empty() {
        return Ok(Vec::new());
    }
    if crate::boolean::assembly::face_components(topo, solid).len() < 2 {
        return Ok(vec![solid]);
    }
    let groups = crate::boolean::assembly::edge_connected_faces(topo, &faces);

    let mut regions = Vec::with_capacity(groups.len());
    for group in groups {
        let shell = topo.add_shell(Shell::new(group).map_err(crate::OperationsError::Topology)?);
        regions.push(topo.add_solid(Solid::new(shell, Vec::new())));
    }

    let mut probes = Vec::with_capacity(regions.len());
    let mut closed = Vec::with_capacity(regions.len());
    for &region in &regions {
        let mut points = Vec::new();
        for vid in brepkit_topology::explorer::solid_vertices(topo, region)? {
            points.push(topo.vertex(vid)?.point());
        }
        probes.push(points);
        closed.push(!crate::boolean::has_free_edges(topo, region)?);
    }

    // `around[i]` lists the closed regions that region `i` lies inside. A
    // probe on the other region's boundary (the two touch there) decides
    // nothing, so the first probe off it does.
    let options = ClassifyOptions {
        tolerance: brepkit_math::tolerance::Tolerance::new().linear,
        ..ClassifyOptions::default()
    };
    // A region inside another lies inside its box, which a NURBS face's
    // control points bound where its sampled extent can fall short; only
    // regions that do are classified.
    let tol = brepkit_math::tolerance::Tolerance::new().linear;
    let mut boxes = Vec::with_capacity(regions.len());
    for &region in &regions {
        let mut bounds = crate::measure::solid_bounding_box(topo, region)?;
        for face in brepkit_topology::explorer::solid_faces(topo, region)? {
            if let brepkit_topology::face::FaceSurface::Nurbs(nurbs) = topo.face(face)?.surface() {
                bounds = bounds.union(nurbs.aabb());
            }
        }
        boxes.push(bounds.expanded(tol));
    }
    let mut around: Vec<Vec<usize>> = vec![Vec::new(); regions.len()];
    for (j, &container) in regions.iter().enumerate() {
        let held: Vec<usize> = (0..regions.len())
            .filter(|&i| {
                closed[j] && i != j && probes[i].iter().all(|&p| boxes[j].contains_point(p))
            })
            .collect();
        if held.is_empty() {
            continue;
        }
        let classifier = SolidClassifier::new(topo, container)?;
        for i in held {
            for &p in &probes[i] {
                match classifier.classify(p, &options)? {
                    PointClassification::Inside => {
                        around[i].push(j);
                        break;
                    }
                    PointClassification::Outside => break,
                    PointClassification::OnBoundary => {}
                }
            }
        }
    }

    let mut pieces = Vec::new();
    for (i, containers) in around.iter().enumerate() {
        if containers.len() % 2 == 0 {
            pieces.push(regions[i]);
            continue;
        }
        let Some(&owner) = containers.iter().max_by_key(|&&j| around[j].len()) else {
            continue;
        };
        let cavity = topo.solid(regions[i])?.outer_shell();
        topo.solid_mut(regions[owner])?.add_inner_shell(cavity);
    }
    Ok(pieces)
}

#[cfg(test)]
mod tests {
    #![allow(clippy::unwrap_used)]

    use brepkit_math::tolerance::Tolerance;
    use brepkit_topology::Topology;
    use brepkit_topology::compound::Compound;

    use super::*;

    #[test]
    fn explode_returns_solids() {
        let mut topo = Topology::new();
        let s1 = crate::primitives::make_box(&mut topo, 1.0, 1.0, 1.0).unwrap();
        let s2 = crate::primitives::make_box(&mut topo, 1.0, 1.0, 1.0).unwrap();
        let cid = topo.add_compound(Compound::new(vec![s1, s2]));

        let solids = explode(&topo, cid).unwrap();
        assert_eq!(solids.len(), 2);
    }

    #[test]
    fn solid_count_works() {
        let mut topo = Topology::new();
        let s1 = crate::primitives::make_box(&mut topo, 1.0, 1.0, 1.0).unwrap();
        let cid = topo.add_compound(Compound::new(vec![s1]));

        assert_eq!(solid_count(&topo, cid).unwrap(), 1);
    }

    #[test]
    fn compound_bbox() {
        let mut topo = Topology::new();
        let s1 = crate::primitives::make_box(&mut topo, 1.0, 1.0, 1.0).unwrap();
        let s2 = crate::primitives::make_box(&mut topo, 1.0, 1.0, 1.0).unwrap();

        crate::transform::transform_solid(
            &mut topo,
            s2,
            &brepkit_math::mat::Mat4::translation(5.0, 0.0, 0.0),
        )
        .unwrap();

        let cid = topo.add_compound(Compound::new(vec![s1, s2]));
        let bb = compound_bounding_box(&topo, cid).unwrap();

        let tol = Tolerance::loose();
        // s1 is [0,1], s2 translated by 5 is [5,6]
        assert!(tol.approx_eq(bb.min.x(), 0.0));
        assert!(tol.approx_eq(bb.max.x(), 6.0));
    }

    /// A sheet-metal bend and the wall it carries share a 50 x 2.5 face, each
    /// placed with the round-off a rigid transform leaves, so the two boxes
    /// stand a rounding error apart. Fused all at once they weld that face
    /// away, as the pair fuse does.
    #[test]
    fn fuse_all_welds_a_face_two_solids_share_after_round_off() {
        use brepkit_math::mat::Mat4;
        use brepkit_math::vec::{Point3, Vec3};

        let mut topo = Topology::new();
        let face = |topo: &mut Topology, pts: &[[f64; 3]]| {
            let pts: Vec<Point3> = pts.iter().map(|p| Point3::new(p[0], p[1], p[2])).collect();
            brepkit_topology::builder::make_planar_face(topo, &pts, 1e-7).unwrap()
        };
        let profile = face(
            &mut topo,
            &[
                [5.0, 55.0, 0.0],
                [4.999_999_999_999_999, 5.0, 0.0],
                [4.999_999_999_999_999, 5.0, 2.5],
                [5.0, 55.0, 2.5],
            ],
        );
        let bend = crate::revolve::revolve(
            &mut topo,
            profile,
            Point3::new(5.0, 55.0, 5.0),
            Vec3::new(1.776_356_839_400_250_5e-17, 1.0, 0.0),
            std::f64::consts::FRAC_PI_2,
        )
        .unwrap();
        let base = face(
            &mut topo,
            &[
                [5.0, 55.732_233_047_033_63, 0.0],
                [-15.0, 55.732_233_047_033_63, 0.0],
                [-15.0, 4.267_766_952_966_369, 0.0],
                [4.999_999_999_999_999, 4.267_766_952_966_369, 0.0],
            ],
        );
        let wall = crate::extrude::extrude(&mut topo, base, Vec3::new(0.0, 0.0, 1.0), 2.5).unwrap();
        let (c, s) = (6.123_233_995_736_766e-17, 1.776_356_839_400_250_2e-17);
        let s2 = 1.776_356_839_400_250_5e-17;
        crate::transform::transform_solid(
            &mut topo,
            wall,
            &Mat4([
                [c, s, 1.0, -8.881_784_197_001_252e-16],
                [s, 1.0, -s2, 0.0],
                [-1.0, s2, c, 10.0],
                [0.0, 0.0, 0.0, 1.0],
            ]),
        )
        .unwrap();

        let area =
            |topo: &Topology, solid| crate::measure::solid_surface_area(topo, solid, 0.01).unwrap();
        let mut pair_topo = topo.clone();
        let pair =
            crate::boolean::boolean(&mut pair_topo, crate::boolean::BooleanOp::Fuse, bend, wall)
                .unwrap();
        let compound = topo.add_compound(Compound::new(vec![bend, wall]));
        let all = fuse_all(&mut topo, compound).unwrap();
        let (pair_area, all_area) = (area(&pair_topo, pair), area(&topo, all));
        assert!(
            (pair_area - 3034.4).abs() < 0.1,
            "pair fuse area {pair_area}"
        );
        assert!(
            (all_area - pair_area).abs() < 1e-6,
            "fused all {all_area}, pair {pair_area}"
        );
    }

    #[test]
    fn fuse_all_two_overlapping_boxes() {
        let mut topo = Topology::new();
        let s1 = crate::primitives::make_box(&mut topo, 1.0, 1.0, 1.0).unwrap();
        let s2 = crate::primitives::make_box(&mut topo, 1.0, 1.0, 1.0).unwrap();

        // Offset s2 slightly — overlapping boxes.
        crate::transform::transform_solid(
            &mut topo,
            s2,
            &brepkit_math::mat::Mat4::translation(0.5, 0.0, 0.0),
        )
        .unwrap();

        let cid = topo.add_compound(Compound::new(vec![s1, s2]));
        let fused = fuse_all(&mut topo, cid).unwrap();

        let vol = crate::measure::solid_volume(&topo, fused, 0.1).unwrap();
        // Two overlapping unit cubes: total should be less than 2.0.
        assert!(
            vol > 1.0 && vol < 2.0,
            "fused volume should be between 1 and 2, got {vol}"
        );
    }

    /// A connected chain of four overlapping unit cubes forms ONE cluster, so
    /// `fuse_all` fuses it via the N-way path. The union is a solid
    /// [0,2.5]×[0,1]×[0,1] bar, so the result must be watertight with volume 2.5.
    #[test]
    fn fuse_all_connected_cluster_is_watertight_bar() {
        use brepkit_math::mat::Mat4;

        let offsets = [0.0, 0.5, 1.0, 1.5];

        // fuse_all (N-way) path.
        let mut topo = Topology::new();
        let boxes: Vec<SolidId> = offsets
            .iter()
            .map(|&dx| {
                let b = crate::primitives::make_box(&mut topo, 1.0, 1.0, 1.0).unwrap();
                crate::transform::transform_solid(&mut topo, b, &Mat4::translation(dx, 0.0, 0.0))
                    .unwrap();
                b
            })
            .collect();
        let cid = topo.add_compound(Compound::new(boxes));
        let fused = fuse_all(&mut topo, cid).unwrap();
        let vol = crate::measure::solid_volume(&topo, fused, 0.01).unwrap();

        // Every edge of a watertight solid is used by exactly two faces.
        let mut uses: std::collections::HashMap<usize, usize> = std::collections::HashMap::new();
        for fid in brepkit_topology::explorer::solid_faces(&topo, fused).unwrap() {
            let face = topo.face(fid).unwrap();
            for wid in std::iter::once(face.outer_wire()).chain(face.inner_wires().iter().copied())
            {
                for oe in topo.wire(wid).unwrap().edges() {
                    *uses.entry(oe.edge().index()).or_default() += 1;
                }
            }
        }
        assert!(
            uses.values().all(|&c| c == 2),
            "fuse_all cluster result must be watertight"
        );
        assert!(
            (vol - 2.5).abs() < 0.01,
            "union of the overlapping row is a [0,2.5] bar (vol 2.5), got {vol}"
        );
    }

    /// Build a hexagonal prism (flat top at z=0..h) centred at the origin via
    /// convex hull — a polyhedral stand-in for a honeycomb pocket.
    fn make_hex_prism(topo: &mut Topology, circumradius: f64, height: f64) -> SolidId {
        use brepkit_math::vec::Point3;
        let mut pts = Vec::with_capacity(12);
        for k in 0..6 {
            let a = std::f64::consts::PI / 3.0 * k as f64;
            let (x, y) = (circumradius * a.cos(), circumradius * a.sin());
            pts.push(Point3::new(x, y, 0.0));
            pts.push(Point3::new(x, y, height));
        }
        crate::primitives::make_convex_hull(topo, &pts).unwrap()
    }

    /// Honeycomb-packed hex prisms with a real gap between every pair, but
    /// corner-to-corner AABBs that overlap. The AABB-only partition collapsed
    /// these into one giant group and unioned them with an O(n) boolean chain;
    /// `partition_touching` proves the gaps with the separating-axis test and
    /// keeps each prism in its own group, so `fuse_all` takes the cheap
    /// disjoint-shell merge.
    #[test]
    fn fuse_all_honeycomb_stays_disjoint() {
        let r = 1.0_f64; // circumradius; across-corners = 2r = 2.0
        let pitch = 2.3_f64; // clear gap on every neighbour, AABBs still overlap
        let height = 4.0_f64;
        let nx = 6;
        let ny = 6;

        let mut topo = Topology::new();
        let mut bboxes = Vec::new();
        let mut solids = Vec::new();
        for j in 0..ny {
            for i in 0..nx {
                let s = make_hex_prism(&mut topo, r, height);
                let x = i as f64 * pitch + (j % 2) as f64 * pitch / 2.0;
                let y = j as f64 * pitch * 0.9;
                crate::transform::transform_solid(
                    &mut topo,
                    s,
                    &brepkit_math::mat::Mat4::translation(x, y, 0.0),
                )
                .unwrap();
                bboxes.push(crate::measure::solid_bounding_box(&topo, s).unwrap());
                solids.push(s);
            }
        }
        let n = solids.len();

        // Every prism is provably disjoint from the others -> one group each.
        let margin = brepkit_math::tolerance::Tolerance::new().linear;
        let pb: Vec<Option<PolyhedralBounds>> = solids
            .iter()
            .map(|&s| polyhedral_bounds(&topo, s))
            .collect();
        let groups = partition_touching(&bboxes, &pb, margin);
        assert_eq!(
            groups.len(),
            n,
            "disjoint hex prisms should each be their own group, got {} groups",
            groups.len()
        );

        // Geometry is still the full disjoint union: volume == n * hex-prism volume.
        let cid = topo.add_compound(Compound::new(solids));
        let fused = fuse_all(&mut topo, cid).unwrap();
        let vol = crate::measure::solid_volume(&topo, fused, 0.05).unwrap();
        let hex_area = 3.0_f64.sqrt() * 1.5 * r * r; // (3*sqrt(3)/2) r^2
        let expected = n as f64 * hex_area * height;
        assert!(
            (vol - expected).abs() < expected * 0.02,
            "fused volume {vol:.2} should match {expected:.2} (n disjoint prisms)"
        );
    }

    fn centred_box(topo: &mut Topology, size: f64, centre: [f64; 3]) -> SolidId {
        let b = crate::primitives::make_box(topo, size, size, size).unwrap();
        let h = size / 2.0;
        crate::transform::transform_solid(
            topo,
            b,
            &brepkit_math::mat::Mat4::translation(centre[0] - h, centre[1] - h, centre[2] - h),
        )
        .unwrap();
        b
    }

    /// Each piece's shell count and volume, smallest piece first.
    fn piece_census(topo: &Topology, pieces: &[SolidId]) -> Vec<(usize, f64)> {
        let mut census: Vec<(usize, f64)> = pieces
            .iter()
            .map(|&p| {
                let shells = 1 + topo.solid(p).unwrap().inner_shells().len();
                (shells, crate::measure::solid_volume(topo, p, 0.01).unwrap())
            })
            .collect();
        census.sort_by(|a, b| a.1.total_cmp(&b.1));
        census
    }

    #[test]
    fn a_slot_through_a_bar_leaves_two_pieces() {
        use crate::boolean::{BooleanOp, boolean};

        let mut topo = Topology::new();
        let bar = crate::primitives::make_box(&mut topo, 80.0, 5.0, 2.5).unwrap();
        let slot = crate::primitives::make_box(&mut topo, 10.0, 20.0, 20.0).unwrap();
        crate::transform::transform_solid(
            &mut topo,
            slot,
            &brepkit_math::mat::Mat4::translation(35.0, -5.0, -5.0),
        )
        .unwrap();
        let cut = boolean(&mut topo, BooleanOp::Cut, bar, slot).unwrap();

        let pieces = solid_pieces(&mut topo, cut).unwrap();
        let census = piece_census(&topo, &pieces);
        assert_eq!(census.len(), 2, "{census:?}");
        for (shells, volume) in census {
            assert_eq!(shells, 1);
            assert!((volume - 437.5).abs() < 1e-6, "{volume}");
        }
    }

    #[test]
    fn a_cube_left_in_a_hollow_boxs_cavity_is_its_own_piece() {
        use crate::boolean::{BooleanOp, boolean};

        let mut topo = Topology::new();
        let outer = centred_box(&mut topo, 20.0, [0.0, 0.0, 10.0]);
        let inner = centred_box(&mut topo, 10.0, [0.0, 0.0, 10.0]);
        let hollow = boolean(&mut topo, BooleanOp::Cut, outer, inner).unwrap();
        assert_eq!(solid_pieces(&mut topo, hollow).unwrap(), vec![hollow]);

        let block = centred_box(&mut topo, 22.0, [0.0, 0.0, 10.0]);
        let left = boolean(&mut topo, BooleanOp::Cut, block, hollow).unwrap();

        let pieces = solid_pieces(&mut topo, left).unwrap();
        let census = piece_census(&topo, &pieces);
        assert_eq!(census.len(), 2, "{census:?}");
        assert_eq!(census[0].0, 1);
        assert!((census[0].1 - 1000.0).abs() < 1e-6, "{census:?}");
        assert_eq!(census[1].0, 2);
        assert!((census[1].1 - 2648.0).abs() < 1e-6, "{census:?}");
    }

    #[test]
    fn a_solid_in_one_piece_comes_back_as_is() {
        let mut topo = Topology::new();
        let b = crate::primitives::make_box(&mut topo, 1.0, 2.0, 3.0).unwrap();
        assert_eq!(solid_pieces(&mut topo, b).unwrap(), vec![b]);
        // A handle the topology does not hold is refused, not handed back.
        assert!(solid_pieces(&mut Topology::new(), b).is_err());
        // The empty solid an intersection of disjoint solids leaves holds no
        // piece.
        let empty = topo.add_empty_solid();
        assert!(solid_pieces(&mut topo, empty).unwrap().is_empty());
    }

    /// A box holding a cavity behind a wall 5e-7 thick, next to a second
    /// box, all in one solid: the cavity is the first box's, by the kernel's
    /// tolerance rather than a coarser one that reads its corners as on the
    /// wall.
    #[test]
    fn a_cavity_behind_a_thin_wall_stays_with_its_piece() {
        use brepkit_topology::shell::Shell;
        use brepkit_topology::solid::Solid;

        let mut topo = Topology::new();
        let wall = 5e-7;
        let outer = centred_box(&mut topo, 10.0, [0.0, 0.0, 0.0]);
        let cavity = centred_box(&mut topo, 2.0f64.mul_add(-wall, 10.0), [0.0, 0.0, 0.0]);
        for f in brepkit_topology::explorer::solid_faces(&topo, cavity).unwrap() {
            let face = topo.face_mut(f).unwrap();
            face.set_reversed(!face.is_reversed());
        }
        let beside = centred_box(&mut topo, 4.0, [20.0, 0.0, 0.0]);
        let mut faces = Vec::new();
        for b in [outer, cavity, beside] {
            faces.extend(brepkit_topology::explorer::solid_faces(&topo, b).unwrap());
        }
        let shell = topo.add_shell(Shell::new(faces).unwrap());
        let all = topo.add_solid(Solid::new(shell, Vec::new()));

        let pieces = solid_pieces(&mut topo, all).unwrap();
        let census = piece_census(&topo, &pieces);
        assert_eq!(census.len(), 2, "{census:?}");
        let hollow = 1000.0 - 2.0f64.mul_add(-wall, 10.0).powi(3);
        assert_eq!(census[0].0, 2);
        assert!((census[0].1 - hollow).abs() < 1e-9, "{census:?}");
        assert_eq!(census[1].0, 1);
        assert!((census[1].1 - 64.0).abs() < 1e-9, "{census:?}");
    }

    /// A ball of radius 10 hollowed to radius 5, fused with a cylinder
    /// standing apart from it: two pieces, each valid and meshing closed,
    /// the ball holding its cavity and the cylinder its own material.
    #[test]
    fn a_hollow_ball_and_a_cylinder_apart_are_two_pieces() {
        use std::f64::consts::PI;

        use brepkit_check::classify::{ClassifyOptions, PointClassification, classify_point};
        use brepkit_math::vec::Point3;

        use crate::boolean::{BooleanOp, boolean};

        let mut topo = Topology::new();
        let outer = crate::primitives::make_sphere(&mut topo, 10.0, 32).unwrap();
        let inner = crate::primitives::make_sphere(&mut topo, 5.0, 32).unwrap();
        let rod = crate::primitives::make_cylinder(&mut topo, 3.0, 8.0).unwrap();
        crate::transform::transform_solid(
            &mut topo,
            rod,
            &brepkit_math::mat::Mat4::translation(30.0, 0.0, -4.0),
        )
        .unwrap();
        // The primitive ball's equator is a polygon of chords, so the ball
        // is measured against its balls as built.
        let [outer_v, inner_v] =
            [outer, inner].map(|b| crate::measure::solid_volume(&topo, b, 0.01).unwrap());
        let hollow = boolean(&mut topo, BooleanOp::Cut, outer, inner).unwrap();
        let both = boolean(&mut topo, BooleanOp::Fuse, hollow, rod).unwrap();

        let pieces = solid_pieces(&mut topo, both).unwrap();
        assert_eq!(pieces.len(), 2);
        let classes = |piece: SolidId| {
            [0.0, 7.5, 30.0].map(|x| {
                classify_point(
                    &topo,
                    piece,
                    Point3::new(x, 0.0, 0.0),
                    &ClassifyOptions::default(),
                )
                .unwrap()
            })
        };
        for &piece in &pieces {
            assert!(
                crate::validate::validate_solid(&topo, piece)
                    .unwrap()
                    .is_valid()
            );
            let mesh = crate::tessellate::tessellate_solid(&topo, piece, 0.01).unwrap();
            assert_eq!(crate::tessellate::boundary_edge_count(&mesh), 0);
            assert_eq!(crate::tessellate::non_manifold_edge_count(&mesh), 0);
            let faces = brepkit_topology::explorer::solid_faces(&topo, piece).unwrap();
            assert!(
                faces
                    .iter()
                    .any(|&f| !topo.face(f).unwrap().surface().is_planar())
            );
            let volume = crate::measure::solid_volume(&topo, piece, 0.01).unwrap();
            let (expected, class) = if topo.solid(piece).unwrap().inner_shells().is_empty() {
                (
                    PI * 9.0 * 8.0,
                    [
                        PointClassification::Outside,
                        PointClassification::Outside,
                        PointClassification::Inside,
                    ],
                )
            } else {
                (
                    outer_v - inner_v,
                    [
                        PointClassification::Outside,
                        PointClassification::Inside,
                        PointClassification::Outside,
                    ],
                )
            };
            assert!(
                (volume - expected).abs() < 1e-3 * expected,
                "volume {volume}, expected {expected}"
            );
            assert_eq!(classes(piece), class);
        }
    }
}