ogeom_bool/

lib.rs

//! Boolean operations.
//!
//! *Elsewhere:* `BOPDS`, `BOPAlgo`, `BOPTools`, `IntTools` and `BRepAlgoAPI`.
//!
//! The architectural insight worth preserving: **fuse, common, cut, section and
//! split are one algorithm (general fuse) plus a selection predicate.** The
//! pipeline is
//!
//! 1. a data structure of indexed sub-shapes, per-type interference lists
//!    (V/V, V/E, V/F, E/E, E/F, F/F), pave blocks and common blocks;
//! 2. a pave filler running in strictly increasing dimension, ending in face/face
//!    intersection, section-edge construction and pcurve generation;
//! 3. a builder that splits faces in 2D parametric space, unifies same-domain
//!    faces, rebuilds solids from face sets, and repairs tolerances;
//! 4. filters that select from the general-fuse result.
//!
//! # One pipeline, any analytic face
//!
//! The pipeline runs on solids whose faces are any surface §7's intersectors
//! answer for (planes, cylinders, spheres, cones, tori), with the planar
//! case as nothing more than the case where every curve is a line. Sections
//! come from [`intersect_surfaces`](ogeom_intersect::intersect_surfaces), exact
//! with same-parameter pcurves where the projection has a closed form and
//! *marched and fitted to a stated tolerance* where it does not. Paves (the
//! points where sections meet boundary edges and each other) come from the
//! exact curve/curve intersection. Every face is then split in its own
//! parameter space by an arrangement of *strands*: polyline scaffolding, each
//! naming the exact sub-curve it stands for, so the combinatorics run on
//! polylines and the rebuilt result is exact curves, pcurves attached both
//! sides, sewn back into shared topology by the sewing whose flipped-carry
//! bug this crate found and §9 fixed.
//!
//! Pieces are classified against the other solid by the exact ray classifier
//! and the filters select: fuse keeps what is outside, common what is inside,
//! cut flips the tool's contribution. Same-domain contact (a piece lying
//! *on* the other boundary) and tangential touching are refused with an
//! error naming the deferred entry, never silently mishandled.

mod arrange;
mod bins;
mod defeature;
mod half_space;

pub use defeature::remove_faces;

use ogeom_algo::{
    Built, Containment, History, is_shell_closed, make_edge_between, make_face_on, make_vertex,
    make_wire, sew, shape_bounds,
};
use ogeom_core::{OgeomResult, Tolerances, ogeom_bail};
use ogeom_geom::Curve2d as _;
use ogeom_geom::Curve3d as _;
use ogeom_geom::Surface as _;
use ogeom_geom::{Curve, PlanarCurve, SurfaceGeometry, Transformable};
use ogeom_math::{Point, Point2};
use ogeom_topo::{
    EdgeRepr, Filter, Location, Model, NodeData, Shape, ShapeType, explore, explore_unique,
};

use arrange::{
    Strand, Traversal, assemble as arrange_pieces, inside_many, inside_many_slanted, inside_rings,
};

/// Parameter-space chord for the polyline scaffolding.
const SCAFFOLD_CHORD: f64 = 1e-3;

/// Parameter-space node snap for the arrangement.
const PARAM_SNAP: f64 = 1e-6;

// --- gathering ---------------------------------------------------------------

/// An edge occurrence's identity: its node and where it is placed.
///
/// One node placed twice is two edges in space: a prism's far cap is its
/// near cap moved, the same nodes under a displacement. Paves shared by
/// node alone would split the near cap's edges wherever the far cap's were
/// crossed. The placement is folded in as its hash, which keeps the key a
/// plain value.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, PartialOrd, Ord)]
struct EdgeKey {
    node: ogeom_topo::TShapeId,
    placement: u64,
}

impl EdgeKey {
    fn of(edge: &Shape) -> Self {
        use core::hash::{Hash as _, Hasher as _};
        let mut hasher = std::collections::hash_map::DefaultHasher::new();
        edge.location().hash(&mut hasher);
        Self {
            node: edge.node(),
            placement: hasher.finish(),
        }
    }

    /// The node's index, for the debug traces.
    fn index(self) -> u32 {
        self.node.index()
    }
}

/// One boundary strand source: an edge as one face uses it, with the pcurve
/// that side of it.
struct BoundaryEdge {
    /// The edge occurrence's identity, for sharing paves across faces.
    node: EdgeKey,
    /// The curve in world space.
    curve: Curve,
    /// The portion the edge covers, in the curve's parameter.
    crange: (f64, f64),
    /// The pcurve on this face's surface.
    pcurve: PlanarCurve,
    /// The portion of the pcurve, mapped proportionally from `crange`.
    prange: (f64, f64),
    /// The other side of a seam, where this edge is one.
    other_side: Option<(PlanarCurve, (f64, f64))>,
    /// The radius within which the edge honestly lies: a fitted rail can
    /// carry a few dozen microns of construction slop, and every filter that
    /// compares this edge's curve against exact geometry widens by it.
    tolerance: f64,
    /// The looser of the edge's two end vertices' own radii: the doubt an
    /// earlier junction recorded there, which the arrangement's node snap
    /// on this face must reach.
    ends_tolerance: f64,
    /// Each end vertex in space with its own recorded radius: a vertex that
    /// owns a span is a junction of this boolean too, and every strand end
    /// inside it names it.
    ends: [(Point, f64); 2],
    /// A conservative box round the edge: its samples, grown by twice the
    /// sag measured between them and by its own radius.
    bound: ogeom_math::Aabb,
}

/// A pole: an edge that bounds a face in parameter space and collapses to
/// a point in space.
///
/// A sphere's poles and a cone's apex have no curve to split, but they are
/// half the chart's boundary; leave them out and the face's outline does
/// not close, and nothing can be arranged inside it.
struct PoleEdge {
    /// Where the pole sits in space.
    point: Point,
    /// The chart line the pole runs along.
    pcurve: PlanarCurve,
    prange: (f64, f64),
}

/// One face, gathered and vetted.
struct GFace {
    face: Shape,
    /// The surface in world space, placement applied.
    surface: SurfaceGeometry,
    /// A conservative world bound: the boundary edges' sampled extent plus
    /// the surface's own allowance: nothing for a plane, measured sampling
    /// slack for the ruled kinds whose rulings pin them to their boundary's
    /// hull, most of the diagonal for anything that can genuinely bulge, a
    /// dome past its equator. The poles join after. Gates the pair filter
    /// and refusals; `OGEOM_BOOL_AUDIT_BOUNDS` audits its conservatism.
    bound: ogeom_math::Aabb,
    /// The scale the marching chord derives from: deliberately *not* the
    /// filter box's diagonal, so the filter can tighten without silently
    /// tightening the marcher.
    chord_scale: f64,
    edges: Vec<BoundaryEdge>,
    poles: Vec<PoleEdge>,
}

/// An argument solid.
struct GSolid {
    solid: Shape,
    faces: Vec<GFace>,
}

fn gather(model: &Model, solid: &Shape, tol: Tolerances) -> OgeomResult<GSolid> {
    if model.kind_of(solid)? != ShapeType::Solid {
        ogeom_bail!(Construction, "boolean arguments are solids");
    }
    for shell in explore_unique(model, solid, ShapeType::Shell)? {
        if !is_shell_closed(model, &shell)? {
            ogeom_bail!(Construction, "an open shell bounds no volume to operate on");
        }
    }

    let mut faces = Vec::new();
    for face in explore(model, solid, Filter::OfType(ShapeType::Face))? {
        let Some(node) = model.node(&face) else {
            ogeom_bail!(Dangling, "face is not in this model");
        };
        let NodeData::Face(data) = node.data() else {
            ogeom_bail!(Construction, "face node holds no face data");
        };
        let Some(stored) = model.geometry().surface(data.surface) else {
            ogeom_bail!(Dangling, "face refers to a surface not in this model");
        };
        let placement = face.transform(model.datums())?;
        if (placement.scale_factor().abs() - 1.0).abs() > 1e-9 {
            ogeom_bail!(
                NotDone,
                "a scaled placement changes a surface's parameterization out \
                 from under its pcurves; bake the scale before a boolean"
            );
        }
        let surface = stored.transformed(&placement, tol)?;
        let surface_id = data.surface;

        let mut edges = Vec::new();
        let mut poles = Vec::new();
        for edge in explore_unique(model, &face, ShapeType::Edge)? {
            let Some(edge_node) = model.node(&edge) else {
                ogeom_bail!(Dangling, "edge is not in this model");
            };
            let NodeData::Edge(edge_data) = edge_node.data() else {
                ogeom_bail!(Construction, "edge node holds no edge data");
            };
            let Some(EdgeRepr::Curve3d { curve, range, .. }) = edge_data.curve3d() else {
                // A degenerate edge (a sphere's pole, a cone's apex) has
                // no extent to split, but it *does* bound the chart, and a
                // chart whose top is missing bounds nothing.
                if let Some(EdgeRepr::PCurve {
                    curve: pc, range, ..
                }) = edge_data.pcurve_for(surface_id, edge.location())
                {
                    let Some(planar) = model.geometry().pcurve(*pc) else {
                        ogeom_bail!(Dangling, "pcurve is not in this model");
                    };
                    let at = explore_unique(model, &edge, ShapeType::Vertex)?;
                    let Some(point) = at
                        .first()
                        .and_then(|v| model.node(v))
                        .and_then(|n| n.data().as_vertex().map(|d| d.point))
                    else {
                        ogeom_bail!(Construction, "a pole with no vertex is nowhere");
                    };
                    let placed = edge.transform(model.datums())?.apply(point);
                    poles.push(PoleEdge {
                        point: placed,
                        pcurve: planar.clone(),
                        prange: *range,
                    });
                }
                continue;
            };
            let Some(geometry) = model.geometry().curve(*curve) else {
                ogeom_bail!(Dangling, "curve is not in this model");
            };
            let world = geometry.transformed(&edge.transform(model.datums())?, tol)?;
            let (pcurve, prange, other_side) =
                match edge_data.pcurve_for(surface_id, edge.location()) {
                    Some(EdgeRepr::PCurve {
                        curve: pc, range, ..
                    }) => {
                        let Some(planar) = model.geometry().pcurve(*pc) else {
                            ogeom_bail!(Dangling, "pcurve is not in this model");
                        };
                        (planar.clone(), *range, None)
                    }
                    Some(EdgeRepr::Seam {
                        forward,
                        reversed,
                        range,
                        ..
                    }) => {
                        let (Some(f), Some(r)) = (
                            model.geometry().pcurve(*forward),
                            model.geometry().pcurve(*reversed),
                        ) else {
                            ogeom_bail!(Dangling, "seam pcurve is not in this model");
                        };
                        (f.clone(), *range, Some((r.clone(), *range)))
                    }
                    _ => ogeom_bail!(
                        Construction,
                        "an edge with no pcurve on its face cannot be split in \
                         that face's parameter space"
                    ),
                };
            let ends_tolerance = model
                .children_of(&edge)?
                .iter()
                .filter_map(|v| model.node(v).and_then(|n| n.data().as_vertex()))
                .fold(0.0_f64, |acc, d| acc.max(d.tolerance.get()));
            // Each end's own vertex, matched to the curve's ends by position:
            // the edge's vertex order and its curve's direction need not
            // agree once the edge is reversed in its wire.
            let ends = {
                let (a, b) = (world.point_at(range.0, tol)?, world.point_at(range.1, tol)?);
                let mut ends = [(a, tol.confusion()), (b, tol.confusion())];
                for v in model.children_of(&edge)? {
                    if let Some(d) = model.node(&v).and_then(|n| n.data().as_vertex()) {
                        let at = v.transform(model.datums())?.apply(d.point);
                        let k = usize::from(at.distance(b) < at.distance(a));
                        ends[k].1 = ends[k].1.max(d.tolerance.get());
                    }
                }
                ends
            };
            edges.push(BoundaryEdge {
                node: EdgeKey::of(&edge),
                curve: world,
                crange: *range,
                pcurve,
                prange,
                other_side,
                tolerance: edge_data.tolerance.get(),
                ends_tolerance,
                ends,
                bound: ogeom_math::Aabb::EMPTY,
            });
        }
        if edges.is_empty() {
            ogeom_bail!(Construction, "a face with no boundary bounds nothing");
        }

        // A seam edge bounds the chart twice *only when this face wraps*:
        // when its other boundary reaches both columns, as a full drum's
        // rims do. A face that merely sits against the meridian (a sphere
        // octant whose boundary happens to be the seam) uses one column,
        // and feeding the far copy into the arrangement leaves a strand
        // nothing connects to. The decision is made here, once, by chart
        // connectivity, and everything downstream (the arrangement, the
        // trim tests, the rebuild) inherits it.
        let seam_indices: Vec<usize> = edges
            .iter()
            .enumerate()
            .filter(|(_, e)| e.other_side.is_some())
            .map(|(i, _)| i)
            .collect();
        for i in seam_indices {
            let mut pool: Vec<Point2> = Vec::new();
            for (j, e) in edges.iter().enumerate() {
                if j == i {
                    continue;
                }
                for t in [e.prange.0, e.prange.1] {
                    pool.push(e.pcurve.point_at(t, tol)?);
                }
                if let Some((other, orange)) = &e.other_side {
                    for t in [orange.0, orange.1] {
                        pool.push(other.point_at(t, tol)?);
                    }
                }
            }
            for p in &poles {
                for t in [p.prange.0, p.prange.1] {
                    pool.push(p.pcurve.point_at(t, tol)?);
                }
            }
            let connected = |pc: &PlanarCurve, range: (f64, f64)| -> OgeomResult<bool> {
                for t in [range.0, range.1] {
                    let at = pc.point_at(t, tol)?;
                    if !pool.iter().any(|q| q.distance(at) <= PARAM_SNAP) {
                        return Ok(false);
                    }
                }
                Ok(true)
            };
            let e = &edges[i];
            let primary_connects = connected(&e.pcurve, e.prange)?;
            let (other_pc, orange) = e
                .other_side
                .clone()
                .unwrap_or_else(|| (e.pcurve.clone(), e.prange));
            let other_connects = connected(&other_pc, orange)?;
            match (primary_connects, other_connects) {
                // Both columns meet the rest of the boundary: the face
                // wraps, and both belong.
                (true, true) => {}
                // One column is this face's; the far copy would dangle.
                (true, false) => edges[i].other_side = None,
                (false, true) => {
                    let e = &mut edges[i];
                    e.pcurve = other_pc;
                    e.prange = orange;
                    e.other_side = None;
                }
                // Neither connects: leave both, and the arrangement's own
                // refusal names the failure as it always did.
                (false, false) => {}
            }
        }
        let mut bound = ogeom_math::Aabb::EMPTY;
        // For a ruled surface the box will be trusted to the boundary's own
        // hull, so the boundary's sampling slack must be measured: how far
        // the true edge sags from the 16-chord polyline, read at each
        // chord's midpoint and doubled for the sag's asymmetry. Nothing else
        // uses the measurement, and the extra evaluations are priced on
        // spline edges, so nothing else pays for it: an unruled face keeps
        // the exact box it always had, and the exact admit set with it.
        let ruled = matches!(
            &surface,
            SurfaceGeometry::Cylinder(_) | SurfaceGeometry::Cone(_)
        );
        let mut slack = 0.0_f64;
        for e in &mut edges {
            let mut previous: Option<Point> = None;
            let (mut own, mut sag) = (ogeom_math::Aabb::EMPTY, 0.0_f64);
            for i in 0..=16 {
                #[allow(clippy::cast_possible_truncation, clippy::cast_precision_loss)]
                let t = e.crange.0 + (e.crange.1 - e.crange.0) * f64::from(i) / 16.0;
                let p = e.curve.point_at(t, tol)?;
                if let Some(q) = previous {
                    let step = (e.crange.1 - e.crange.0) / 32.0;
                    let mid = e.curve.point_at(t - step, tol)?;
                    let off = mid.distance(Point::midpoint(q, p)) * 2.0;
                    sag = sag.max(off);
                    own = own.with_point(mid);
                    if ruled {
                        slack = slack.max(off);
                        bound = bound.with_point(mid);
                    }
                }
                own = own.with_point(p);
                bound = bound.with_point(p);
                previous = Some(p);
            }
            e.bound = own.expanded(sag + e.tolerance + tol.confusion() * 1e2);
        }
        // A plane never bulges past its boundary. A ruled surface (cylinder,
        // cone) cannot either: every surface point lies on a straight ruling
        // whose ends are on the boundary, so the face sits inside its
        // boundary's hull and only the boundary's own sampling slack is owed.
        // Anything else may genuinely bulge (a dome past its equator) and
        // keeps most of its own diagonal as allowance. The audit behind
        // OGEOM_BOOL_AUDIT_BOUNDS holds every arm to conservatism.
        let bulge = match &surface {
            SurfaceGeometry::Plane(_) => 0.0,
            SurfaceGeometry::Cylinder(_) | SurfaceGeometry::Cone(_) => slack,
            _ => bound.diagonal() * 0.75,
        };
        // The scale the marching chord is derived from, decoupled from the
        // filter box. It reproduces exactly what the heuristic was tuned
        // against (the diagonal as the blanket three-quarter bulge left it)
        // because the chord is a tolerance, not a bound: tightening the
        // filter must not silently tighten the marcher, which is exactly
        // what a single box feeding both does.
        let margin = tol.confusion() * 1e2;
        let chord_scale = match &surface {
            SurfaceGeometry::Plane(_) => bound.expanded(margin).diagonal(),
            _ => bound.expanded(bound.diagonal() * 0.75 + margin).diagonal(),
        };
        // A pole bounds the chart with no edge to sample: a cone drilled to
        // its apex reaches the apex, and a bound that omits it would let the
        // filter drop a pair the apex genuinely meets. It joins *after* the
        // bulge is taken from the boundary's own diagonal: the pole is an
        // exact point, and letting it stretch the diagonal would inflate a
        // dome's allowance by its own height over again.
        let mut bound = bound.expanded(bulge);
        for pole in &poles {
            bound = bound.with_point(pole.point);
        }
        let bound = bound.expanded(tol.confusion() * 1e2);
        faces.push(GFace {
            poles,
            face,
            surface,
            bound,
            chord_scale,
            edges,
        });
    }
    if faces.is_empty() {
        ogeom_bail!(Construction, "a solid with no faces bounds nothing");
    }
    Ok(GSolid {
        solid: solid.clone(),
        faces,
    })
}

/// A parameter brought onto the turn its edge actually covers.
///
/// A periodic curve's own domain and an edge's range over it need not agree
/// on which turn to count from: a sphere's seam runs its half meridian over
/// `[-π/2, π/2]`, while every crossing found on it comes back in `[0, 2π)`.
/// Left alone, a crossing at latitude `-0.96` arrives as `5.32`, sits outside
/// the range, and is discarded, so the seam never splits where a section
/// genuinely meets it, and the arrangement finds the chain hanging.
fn onto_range(t: f64, curve: &Curve, range: (f64, f64), tol: Tolerances) -> f64 {
    if !curve.is_periodic() {
        return t;
    }
    let (a, b) = curve.domain();
    let period = b - a;
    if period <= 0.0 {
        return t;
    }
    for k in [0.0, -1.0, 1.0, -2.0, 2.0] {
        let shifted = period.mul_add(k, t);
        if shifted >= range.0 - tol.parametric() && shifted <= range.1 + tol.parametric() {
            return shifted;
        }
    }
    t
}

/// The part of an overlap that falls within the second curve's own bounded
/// range, stated in the *first* curve's parameter.
///
/// An overlap is between two curves; an edge covers only part of its curve.
/// The correspondence the overlap states is affine, so the edge's range
/// carries across as an interval, and on a periodic curve it carries across
/// up to whole turns, so the shift that meets the overlap is the one meant.
/// `None` where the edge's own stretch and the overlap do not meet.
fn overlap_within(
    overlap: &ogeom_intersect::Overlap,
    range: (f64, f64),
    curve: &Curve,
    tol: Tolerances,
) -> Option<(f64, f64)> {
    let ordered = |r: (f64, f64)| if r.0 <= r.1 { r } else { (r.1, r.0) };
    let span_a = overlap.on_a.1 - overlap.on_a.0;
    let span_b = overlap.on_b.1 - overlap.on_b.0;
    if span_b.abs() <= f64::MIN_POSITIVE {
        return None;
    }
    let to_a = |t: f64| overlap.on_a.0 + span_a * (t - overlap.on_b.0) / span_b;
    let (wlo, whi) = ordered((to_a(range.0), to_a(range.1)));
    let (lo, hi) = ordered(overlap.on_a);
    let period = if curve.is_periodic() {
        let (a, b) = curve.domain();
        // The turn measured in the *first* curve's parameter, which is where
        // the intersection is being taken.
        (b - a) * (span_a / span_b).abs()
    } else {
        0.0
    };
    let mut best: Option<(f64, f64)> = None;
    for k in [0.0, 1.0, -1.0, 2.0, -2.0] {
        let candidate = (
            lo.max(period.mul_add(k, wlo)),
            hi.min(period.mul_add(k, whi)),
        );
        if candidate.1 - candidate.0 > best.map_or(tol.parametric(), |(x, y)| y - x) {
            best = Some(candidate);
        }
        if period == 0.0 {
            break;
        }
    }
    best
}

/// A parameter carried proportionally from one range to another.
fn rescale(t: f64, from: (f64, f64), to: (f64, f64)) -> f64 {
    let span = from.1 - from.0;
    if span.abs() <= f64::MIN_POSITIVE {
        return to.0;
    }
    to.0 + (to.1 - to.0) * (t - from.0) / span
}

/// The pcurve's course over a sub-range of the *curve's* parameters.
/// The stretches of `curve` over `crange` that lie along the edge `e`
/// within the pair's honesty, as overlaps with the edge's own parameters at
/// their ends: the measured twin of the intersector's closed-form
/// coincidence, for a fitted curve on an exact one.
fn measured_overlaps(
    curve: &Curve,
    crange: (f64, f64),
    tolerance: f64,
    e: &BoundaryEdge,
    tol: Tolerances,
) -> OgeomResult<Vec<ogeom_intersect::Overlap>> {
    const SAMPLES: usize = 48;
    let width = (tolerance.max(e.tolerance) * 2.0).max(tol.confusion() * 1e3);
    let mut near: Vec<bool> = Vec::with_capacity(SAMPLES + 1);
    for i in 0..=SAMPLES {
        #[allow(clippy::cast_precision_loss)]
        let t = crange.0 + (crange.1 - crange.0) * (i as f64) / (SAMPLES as f64);
        let p = curve.point_at(t, tol)?;
        near.push(distance_to_edge_curve(&e.curve, e.crange, p, tol)? <= width);
    }
    let mut out = Vec::new();
    let mut i = 0;
    while i <= SAMPLES {
        if !near[i] {
            i += 1;
            continue;
        }
        let start = i;
        while i < SAMPLES && near[i + 1] {
            i += 1;
        }
        let end = i;
        i += 1;
        // Three samples along is a stretch; fewer is a crossing's blur.
        if end - start < 2 {
            continue;
        }
        #[allow(clippy::cast_precision_loss)]
        let at = |k: usize| crange.0 + (crange.1 - crange.0) * (k as f64) / (SAMPLES as f64);
        // Each end refined by bisection into the neighbouring sample gap.
        let refine = |inside: f64, outside: f64| -> OgeomResult<f64> {
            let (mut a, mut b) = (inside, outside);
            for _ in 0..24 {
                let m = f64::midpoint(a, b);
                let p = curve.point_at(m, tol)?;
                if distance_to_edge_curve(&e.curve, e.crange, p, tol)? <= width {
                    a = m;
                } else {
                    b = m;
                }
            }
            Ok(a)
        };
        let lo = if start == 0 {
            at(0)
        } else {
            refine(at(start), at(start - 1))?
        };
        let hi = if end == SAMPLES {
            at(SAMPLES)
        } else {
            refine(at(end), at(end + 1))?
        };
        if hi - lo <= tol.parametric() {
            continue;
        }
        let on_edge = |t: f64| -> OgeomResult<f64> {
            let p = curve.point_at(t, tol)?;
            let foot = ogeom_algo::project_on_curve(&e.curve, p, 64, tol)?;
            Ok(onto_range(foot.parameter, &e.curve, e.crange, tol))
        };
        if *DEBUG_WIRE {
            eprintln!(
                "CONTACT along edge {} measured over ({lo:.6}, {hi:.6}) within {width:.2e}",
                e.node.index()
            );
        }
        out.push(ogeom_intersect::Overlap {
            on_a: (lo, hi),
            on_b: (on_edge(lo)?, on_edge(hi)?),
        });
    }
    Ok(out)
}

fn pcurve_polyline(
    pcurve: &PlanarCurve,
    prange: (f64, f64),
    crange: (f64, f64),
    sub: (f64, f64),
    surface: &SurfaceGeometry,
    tol: Tolerances,
) -> OgeomResult<Vec<Point2>> {
    let lo = rescale(sub.0, crange, prange);
    let hi = rescale(sub.1, crange, prange);
    // Enough samples that the first step approximates the tangent and the
    // scanline interior test has a faithful outline. Straight pcurves get
    // two points; everything else a fixed fine sampling.
    let mut count = match pcurve {
        PlanarCurve::Line(_) => 1,
        _ => {
            let span = (hi - lo).abs().max(SCAFFOLD_CHORD);
            #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
            let n = (span / SCAFFOLD_CHORD).sqrt().ceil() as usize;
            n.clamp(8, 256)
        }
    };
    // And never a step [`unwrap_polyline`] could read as a period jump. A
    // straight pcurve is faithful at two points (a coaxial rim's image on
    // a cylinder is exactly a line across the chart), but one that crosses
    // more than half the turn looks to that reader like a wrap, and the
    // strand comes back running the complementary way: a bore's mouth
    // blended where the mouth circle starts anywhere but the wall's own
    // seam. The chart's span decides the count, not the parameter's.
    let ((ua, ub), (va, vb)) = surface.domain();
    let (first, last) = (pcurve.point_at(lo, tol)?, pcurve.point_at(hi, tol)?);
    for (periodic, period, reach) in [
        (surface.is_periodic_u(), ub - ua, (last.x - first.x).abs()),
        (surface.is_periodic_v(), vb - va, (last.y - first.y).abs()),
    ] {
        if !periodic || period <= 0.0 || !reach.is_finite() {
            continue;
        }
        let steps = (reach / (period / 3.0)).ceil();
        #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
        let steps = steps.clamp(1.0, 4096.0) as usize;
        count = count.max(steps);
    }
    let mut out = Vec::with_capacity(count + 1);
    for i in 0..=count {
        #[allow(clippy::cast_precision_loss)]
        let t = lo + (hi - lo) * i as f64 / count as f64;
        out.push(pcurve.point_at(t, tol)?);
    }
    Ok(out)
}

// --- the filler --------------------------------------------------------------

/// One section curve between a face of `a` and a face of `b`.
struct SectionRec {
    curve: Curve,
    /// The pcurve on the A face's surface, sharing the curve's parameter.
    pc_a: PlanarCurve,
    /// The same on the B face.
    pc_b: PlanarCurve,
    face_a: usize,
    face_b: usize,
    closed: bool,
    /// How far the curve may sit from the true intersection: zero for an
    /// exact section, the trace-plus-fit budget for a marched one. Crossing
    /// filters widen their acceptance by this, or a fitted section would
    /// never register against the edges it genuinely meets.
    tolerance: f64,
}

/// The parameter intervals of a contact curve over which *both* faces
/// actually reach it.
///
/// Two surfaces touch along the whole of their contact; two faces touch
/// along whatever part of it their trims both hold. That part is found by
/// sampling (sixty-four stations along the curve, each asked of both
/// charts) rather than by intersecting the contact with the boundary
/// edges, because a contact meets those boundaries tangentially too and the
/// crossing finder is the wrong instrument for it. The cost of sampling is
/// the usual one: a stretch shorter than a station can be missed, and an
/// endpoint is placed within a station of the truth.
fn contact_intervals(
    fused: &GeneralFused,
    contact: &TangentRec,
    tol: Tolerances,
) -> OgeomResult<Vec<(f64, f64)>> {
    let outline = |face: &GFace| -> OgeomResult<Vec<Vec<Point2>>> {
        let mut lines = Vec::new();
        for e in &face.edges {
            lines.push(pcurve_polyline(
                &e.pcurve,
                e.prange,
                e.crange,
                e.crange,
                &face.surface,
                tol,
            )?);
        }
        Ok(lines)
    };
    let rings_a = outline(&fused.a.faces[contact.face_a])?;
    let rings_b = outline(&fused.b.faces[contact.face_b])?;
    let refs_a: Vec<&[Point2]> = rings_a.iter().map(Vec::as_slice).collect();
    let refs_b: Vec<&[Point2]> = rings_b.iter().map(Vec::as_slice).collect();

    const STATIONS: usize = 64;
    let domain = contact.curve.domain();
    let span = domain.1 - domain.0;
    let mut runs: Vec<(f64, f64)> = Vec::new();
    let mut open: Option<f64> = None;
    for k in 0..=STATIONS {
        #[expect(
            clippy::cast_precision_loss,
            reason = "a station index, far below the mantissa"
        )]
        let t = span.mul_add(k as f64 / STATIONS as f64, domain.0);
        // A periodic chart's parameters run out past its own window; the
        // trim test is only meaningful once they are folded back into it.
        let held = matches!(
            (
                contact.pc_a.point_at(t, tol),
                contact.pc_b.point_at(t, tol)
            ),
            (Ok(pa), Ok(pb))
                if arrange::inside_many(
                    &refs_a,
                    fold_point_into_chart(pa, &fused.a.faces[contact.face_a].surface),
                ) && arrange::inside_many(
                    &refs_b,
                    fold_point_into_chart(pb, &fused.b.faces[contact.face_b].surface),
                )
        );
        match (held, open) {
            (true, None) => open = Some(t),
            (false, Some(from)) => {
                if t - from > tol.parametric() {
                    runs.push((from, t));
                }
                open = None;
            }
            _ => {}
        }
    }
    if let Some(from) = open
        && domain.1 - from > tol.parametric()
    {
        runs.push((from, domain.1));
    }
    Ok(runs)
}

/// One curve along which two faces *touch* without crossing.
///
/// A contact carries no boundary parity (neither face passes through the
/// other), so it takes no part in the arrangement or the classification.
/// It is still a curve that exists on the result, and a section view
/// through a tangency has to show it, so it is carried alongside.
struct TangentRec {
    curve: Curve,
    pc_a: PlanarCurve,
    pc_b: PlanarCurve,
    face_a: usize,
    face_b: usize,
}

/// One kept sub-range of one section.
/// A strand's tolerance as a junction may trust it: a fitted section whose
/// trace failed reports a budget of metres, and a weld that believed it
/// would join every vertex of the model. Nothing this pipeline fits is
/// honestly looser than ten thousand confusions.
fn honest(tolerance: f64, tol: Tolerances) -> f64 {
    tolerance.min(tol.confusion() * 1e4)
}

/// One split an edge is asked for, with how honestly it can be placed.
#[derive(Debug, Clone, Copy)]
struct Pave {
    /// The parameter on the edge's curve.
    t: f64,
    /// How far the strand that asked for it may honestly sit from the exact
    /// junction: its own tolerance, and its tangential doubt where the
    /// crossing was a touch.
    honesty: f64,
}

/// Paves the edge cannot tell apart, as one junction each.
///
/// A tolerant rail meeting a wedge's faces near one corner collects a
/// cluster of crossings inside its own stated radius, and a fitted section
/// and a fitted contact asking for the same corner land a few microns
/// apart. Split at each, the edge shatters into dust no weld downstream can
/// rejoin. So consecutive paves whose gap along the edge in space is within
/// the edge's honesty or either pave's own are one cluster, its first pave
/// speaking for it.
#[derive(Debug, Clone, Copy)]
struct PaveCluster {
    /// The representative parameter: the cluster's first pave.
    t: f64,
    /// Where it sits.
    at: Point,
    /// How far the cluster's members reach from the representative.
    span: f64,
    /// The loosest honesty among the members.
    honesty: f64,
    /// How many paves the cluster holds.
    members: usize,
}

fn cluster_paves(
    curve: &Curve,
    crange: (f64, f64),
    edge_tolerance: f64,
    paves: &[Pave],
    tol: Tolerances,
) -> OgeomResult<Vec<PaveCluster>> {
    let mut ts: Vec<Pave> = paves
        .iter()
        .copied()
        .filter(|p| p.t > crange.0 + tol.parametric() && p.t < crange.1 - tol.parametric())
        .collect();
    ts.sort_by(|x, y| x.t.partial_cmp(&y.t).unwrap_or(core::cmp::Ordering::Equal));
    ts.dedup_by(|x, y| {
        if (x.t - y.t).abs() <= tol.parametric() {
            y.honesty = y.honesty.max(x.honesty);
            true
        } else {
            false
        }
    });
    let floor = edge_tolerance.max(tol.confusion() * 10.0);
    // A pave within the edge's honesty of its own end *is* the end: split
    // there, the sliver between them is dust one face keeps and the
    // face across the edge drops, and the sew finds it used once.
    let ends = [
        curve.point_at(crange.0, tol)?,
        curve.point_at(crange.1, tol)?,
    ];
    let mut clusters: Vec<PaveCluster> = Vec::new();
    let mut prev: Option<(Point, f64)> = None;
    for pave in ts {
        let at = curve.point_at(pave.t, tol)?;
        if ends
            .iter()
            .any(|e| e.distance(at) <= floor.max(pave.honesty))
        {
            continue;
        }
        let joined = prev.is_some_and(|(held, honesty): (Point, f64)| {
            held.distance(at) <= floor.max(honesty).max(pave.honesty)
        });
        if joined && let Some(cluster) = clusters.last_mut() {
            cluster.span = cluster.span.max(cluster.at.distance(at));
            cluster.honesty = cluster.honesty.max(pave.honesty);
            cluster.members += 1;
        } else {
            clusters.push(PaveCluster {
                t: pave.t,
                at,
                span: 0.0,
                honesty: pave.honesty,
                members: 1,
            });
        }
        prev = Some((at, pave.honesty));
    }
    Ok(clusters)
}

#[derive(Clone)]
struct SectionPiece {
    section: usize,
    range: (f64, f64),
    /// Whether this sub-range already *is* boundary on the A face, the B
    /// face, or neither. A piece that is boundary on one side is still the
    /// splitting curve on the other, and only the face it duplicates leaves
    /// it out.
    hugs: [bool; 2],
    /// Set when the piece was admitted by a hug: a split of the other face
    /// along this edge of the hugging face: (edge node, target from A,
    /// target face). Two sections hugging one edge onto one face are the
    /// same split, kept once.
    hug_key: Option<(EdgeKey, bool, usize)>,
}

/// One boundary edge of one argument's face, lying in a face of the other
/// argument's own surface: the splitting curve same-domain contact
/// contributes, since coincident surfaces have no section curve to offer.
struct ContactRec {
    /// The owner's world curve and the sub-range its edge covers.
    curve: Curve,
    crange: (f64, f64),
    /// The curve spoken in the *target* face's chart, over its own window,
    /// mapped proportionally from `crange` exactly as an edge's own stored
    /// pcurve is.
    pcurve: PlanarCurve,
    prange: (f64, f64),
    /// The owner edge's node, whose paves this record shares.
    node: EdgeKey,
    /// The owner edge's own tolerance: how far its curve may honestly sit
    /// from the exact geometry it meets, which is how far a crossing filter
    /// must reach to see a fitted rail cross an exact seam.
    tolerance: f64,
    /// The target: which argument's face list, and which face.
    target_from_a: bool,
    target_face: usize,
}

/// A same-domain contact edge's image in the shared surface's chart, where
/// no closed form exists: fitted by projection into the target's chart, or
/// (where the edge overhangs the target's own window, a blend's leg on a
/// spline host continued past the face) into the owner's, which is the
/// same chart wherever an extension kept its parent's parameters. The
/// owner's image is trusted only where lifting it through the *target*
/// lands back on the edge, sampled over the run the target's window holds.
///
/// `None` for an edge lying wholly outside the target's window, which the
/// target face cannot meet.
///
/// # Errors
///
/// [`OgeomError::NotDone`](ogeom_core::OgeomError::NotDone) if neither
/// chart holds the edge.
fn projected_into_shared_chart(
    curve: &Curve,
    range: (f64, f64),
    owner: &SurfaceGeometry,
    target: &SurfaceGeometry,
    tol: Tolerances,
) -> OgeomResult<Option<PlanarCurve>> {
    let against_target = ogeom_algo::pcurve_fit::fit_projected_pcurve(curve, range, target, tol);
    if let Ok(fitted) = &against_target
        && fitted.2
    {
        return Ok(Some(fitted.0.clone()));
    }
    let Ok(on_owner) = ogeom_algo::pcurve_fit::fit_projected_pcurve(curve, range, owner, tol)
    else {
        ogeom_bail!(
            NotDone,
            "same-domain contact whose edge could not be projected into the \
             shared surface's chart: {}",
            against_target
                .err()
                .map_or_else(String::new, |e| e.to_string())
        );
    };
    let ((ua, ub), (va, vb)) = target.domain();
    let mut checked = 0usize;
    for i in 0..=16 {
        let t = range.0 + (range.1 - range.0) * f64::from(i) / 16.0;
        let uv = on_owner.0.point_at(t, tol)?;
        if uv.x < ua || uv.x > ub || uv.y < va || uv.y > vb {
            continue;
        }
        let lifted = target.point_at(uv.x, uv.y, tol)?;
        if lifted.distance(curve.point_at(t, tol)?) > tol.confusion() * 1e4 {
            ogeom_bail!(
                NotDone,
                "same-domain contact whose edge overhangs the shared surface's \
                 window, on a chart the two surfaces do not share"
            );
        }
        checked += 1;
    }
    if checked == 0 {
        return Ok(None);
    }
    Ok(Some(on_owner.0))
}

/// Whether two surfaces are the *identical chart*: the same
/// parameterization, frame and all, not merely the same point set.
///
/// The five analytics only: a fitted surface never qualifies, because two
/// independent fits of one geometry agree nowhere in parameter space. This
/// is what lets a stored pcurve stand in for a closed-form projection in the
/// same-domain melt: on the identical chart, the owner's pcurve already is
/// the projection.
fn same_chart(a: &SurfaceGeometry, b: &SurfaceGeometry, tol: Tolerances) -> bool {
    use SurfaceGeometry as S;
    let frames = |fa: ogeom_math::Frame, fb: ogeom_math::Frame| -> bool {
        fa.origin().distance(fb.origin()) <= tol.confusion()
            && fa.z().vector().dot(fb.z().vector()) >= 1.0 - tol.angular()
            && fa.x().vector().dot(fb.x().vector()) >= 1.0 - tol.angular()
    };
    match (a, b) {
        (S::Plane(x), S::Plane(y)) => frames(x.plane().frame(), y.plane().frame()),
        (S::Cylinder(x), S::Cylinder(y)) => {
            frames(x.cylinder().frame(), y.cylinder().frame())
                && (x.cylinder().radius() - y.cylinder().radius()).abs() <= tol.confusion()
        }
        (S::Cone(x), S::Cone(y)) => {
            frames(x.cone().frame(), y.cone().frame())
                && (x.cone().reference_radius() - y.cone().reference_radius()).abs()
                    <= tol.confusion()
                && (x.cone().half_angle() - y.cone().half_angle()).abs() <= tol.angular()
        }
        (S::Sphere(x), S::Sphere(y)) => {
            frames(x.sphere().frame(), y.sphere().frame())
                && (x.sphere().radius() - y.sphere().radius()).abs() <= tol.confusion()
        }
        (S::Torus(x), S::Torus(y)) => {
            frames(x.torus().frame(), y.torus().frame())
                && (x.torus().major_radius() - y.torus().major_radius()).abs() <= tol.confusion()
                && (x.torus().minor_radius() - y.torus().minor_radius()).abs() <= tol.confusion()
        }
        // Two patches are one chart only as one patch: the same knots and
        // the same net, to the bit: a blend's leg built on the host's own
        // patch, widened in place and shared. Any other pair of patches
        // that merely coincides is not.
        (S::BSpline(x), S::BSpline(y)) => x == y,
        _ => false,
    }
}

/// What a face's arrangement strand stands for.
#[derive(Clone)]
enum Tag {
    /// A sub-range of boundary edge `edge` (index into the face's edges).
    Boundary { edge: usize, range: (f64, f64) },
    /// A sub-range of a global section.
    Section { section: usize, range: (f64, f64) },
    /// A sub-range of a global contact edge: another face's boundary edge
    /// lying in this face's own surface, splitting it.
    Contact { contact: usize, range: (f64, f64) },
    /// A sub-range of a pole of the face's own chart. The edge is a point in
    /// space whatever the range says; the range is where it runs in the
    /// chart, which is the only place it has length.
    Pole { pole: usize, range: (f64, f64) },
}

/// Where a piece stands relative to the other solid.
///
/// `In` and `Out` are the classifier's words. The two `On` states split the
/// case the classifier cannot decide alone: a piece lying on the other
/// boundary bounds material on one side here and one side there, and whether
/// those sides agree (outward normals aligned) or oppose is what every
/// operation's filter turns on.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum PieceState {
    In,
    Out,
    /// On the other boundary, outward normals aligned: both solids' material
    /// on the same side. One copy of the piece bounds the union and the
    /// intersection alike.
    OnAligned,
    /// On the other boundary, outward normals opposed: material on both
    /// sides. The contact is interior to the union and vanishes from it.
    OnOpposed,
}

/// Fold a param-space polyline into a periodic surface's chart, by one
/// constant offset per axis.
///
/// A section's pcurve is *unwrapped* across a seam: continuous, and allowed
/// to leave the stated domain, because that is what crossing a seam is. The
/// arrangement lives in one chart, and the filler has already split every
/// section at its seam crossings, so each strand spans at most one chart
/// width and a single period shift per axis brings it home. The shift is
/// chosen by the strand's midpoint, so endpoints sitting exactly on the
/// chart's edge stay on whichever side the strand's body is.
/// A point in a polyline's *interior*: its half-way member, except that a
/// straight strand is two points and its half-way member is an endpoint,
/// which may sit exactly on a seam or a trim. The chord midpoint is on the
/// curve for a straight image and strictly inside either way.
fn interior_of(line: &[Point2]) -> Point2 {
    if line.len() == 2 {
        Point2::new(
            f64::midpoint(line[0].x, line[1].x),
            f64::midpoint(line[0].y, line[1].y),
        )
    } else {
        line[line.len() / 2]
    }
}

fn fold_into_chart(line: &mut [Point2], surface: &SurfaceGeometry) {
    let ((ua, ub), (va, vb)) = surface.domain();
    if line.is_empty() {
        return;
    }
    let mid = interior_of(line);
    if surface.is_periodic_u() {
        let span = ub - ua;
        if span > 0.0 {
            let shift = (ua + (mid.x - ua).rem_euclid(span)) - mid.x;
            for p in line.iter_mut() {
                p.x += shift;
            }
        }
    }
    if surface.is_periodic_v() {
        let span = vb - va;
        if span > 0.0 {
            let shift = (va + (mid.y - va).rem_euclid(span)) - mid.y;
            for p in line.iter_mut() {
                p.y += shift;
            }
        }
    }
}

/// The shifts by whole periods worth trying on a chart point after its fold,
/// nearest first: none, then one period along each periodic direction.
fn period_shifts(surface: &SurfaceGeometry) -> Vec<(f64, f64)> {
    let ((ua, ub), (va, vb)) = surface.domain();
    let mut out = vec![(0.0, 0.0)];
    if surface.is_periodic_u() && ub > ua {
        out.push((ub - ua, 0.0));
        out.push((ua - ub, 0.0));
    }
    if surface.is_periodic_v() && vb > va {
        out.push((0.0, vb - va));
        out.push((0.0, va - vb));
    }
    out
}

/// A chart point folded to the side of the face's seam its trim is on.
///
/// The surface's fold puts a point in the period starting at the surface's
/// own zero, which is where a face's trim is only when its seam runs there.
/// A patch whose seam was cut elsewhere (a band opened along the widest gap
/// its mesh left, running diagonally round) covers a period starting at the
/// seam, and a point just short of it belongs a period on. The fold is kept
/// where it lands inside the trim, or where no shift does.
fn fold_inside(p: Point2, surface: &SurfaceGeometry, trim: &[&[Point2]]) -> Point2 {
    let folded = fold_point_into_chart(p, surface);
    if trim.is_empty() {
        return folded;
    }
    period_shifts(surface)
        .into_iter()
        .map(|(du, dv)| Point2::new(folded.x + du, folded.y + dv))
        .find(|q| inside_many(trim, *q))
        .unwrap_or(folded)
}

/// A strand folded into the chart as a whole, to the side of the face's
/// seam its interior lies inside the trim on, as [`fold_inside`] does for
/// a point.
fn fold_line_inside(line: &mut [Point2], surface: &SurfaceGeometry, trim: &[&[Point2]]) {
    fold_into_chart(line, surface);
    if line.is_empty() || trim.is_empty() {
        return;
    }
    let mid = interior_of(line);
    if let Some((du, dv)) = period_shifts(surface)
        .into_iter()
        .find(|(du, dv)| inside_many(trim, Point2::new(mid.x + du, mid.y + dv)))
    {
        for p in line.iter_mut() {
            p.x += du;
            p.y += dv;
        }
    }
}

/// A face's boundary in its chart, each edge's image sampled, both sides of
/// a seam: enough to tell which side of the seam a point is on.
fn face_trim_lines(face: &GFace, tol: Tolerances) -> Vec<Vec<Point2>> {
    let sample = |pcurve: &PlanarCurve, range: (f64, f64)| -> Vec<Point2> {
        (0..=16)
            .filter_map(|k| {
                pcurve
                    .point_at(range.0 + (range.1 - range.0) * f64::from(k) / 16.0, tol)
                    .ok()
            })
            .collect()
    };
    let mut out = Vec::new();
    for e in &face.edges {
        out.push(sample(&e.pcurve, e.prange));
        if let Some((other, range)) = &e.other_side {
            out.push(sample(other, *range));
        }
    }
    out
}

/// Remove period tears from a sampled polyline, axis by axis.
fn unwrap_polyline(line: &mut [Point2], surface: &SurfaceGeometry) {
    let ((ua, ub), (va, vb)) = surface.domain();
    let spans = (
        if surface.is_periodic_u() {
            ub - ua
        } else {
            0.0
        },
        if surface.is_periodic_v() {
            vb - va
        } else {
            0.0
        },
    );
    for i in 1..line.len() {
        if spans.0 > 0.0 {
            while line[i].x - line[i - 1].x > spans.0 * 0.5 {
                line[i].x -= spans.0;
            }
            while line[i].x - line[i - 1].x < -spans.0 * 0.5 {
                line[i].x += spans.0;
            }
        }
        if spans.1 > 0.0 {
            while line[i].y - line[i - 1].y > spans.1 * 0.5 {
                line[i].y -= spans.1;
            }
            while line[i].y - line[i - 1].y < -spans.1 * 0.5 {
                line[i].y += spans.1;
            }
        }
    }
}

/// Fold a point into a periodic surface's chart.
fn fold_point_into_chart(p: Point2, surface: &SurfaceGeometry) -> Point2 {
    let mut one = [p];
    fold_into_chart(&mut one, surface);
    one[0]
}

/// A sub-range of a closed curve, brought into its domain.
///
/// The filler split every wrap interval at the domain end, so a piece fits
/// within one period; the fold of its start may still land the end a hair
/// past the domain, which clamps.
fn folded_range(range: (f64, f64), domain: (f64, f64), closed: bool) -> (f64, f64) {
    if !closed {
        return range;
    }
    let f0 = fold(range.0, domain);
    let f1 = (f0 + (range.1 - range.0)).min(domain.1);
    (f0, f1)
}

/// Fold a parameter into a closed curve's domain.
fn fold(t: f64, domain: (f64, f64)) -> f64 {
    let span = domain.1 - domain.0;
    if span <= 0.0 {
        return domain.0;
    }
    domain.0 + (t - domain.0).rem_euclid(span)
}

/// The parameter to evaluate a section at: folded where the curve is closed
/// and its pieces may run past the domain end, left alone where it is not.
///
/// An open section's own end *is* the domain end, and folding it lands on the
/// domain start instead: the far end of the curve. On a plane's section
/// through a ball's poles, where the marcher hands back two open half circles
/// each running pole to pole, that is the difference between an edge bounded
/// by the pole it reaches and one bounded by the pole on the other side.
fn at_param(t: f64, domain: (f64, f64), closed: bool) -> f64 {
    if closed { fold(t, domain) } else { t }
}

/// The sections between two gathered solids, with the paves they put on
/// boundary edges.
#[allow(clippy::type_complexity)]
fn fill(
    ga: &GSolid,
    gb: &GSolid,
    admit_all: bool,
    tol: Tolerances,
) -> OgeomResult<(
    Vec<SectionRec>,
    Vec<SectionPiece>,
    Vec<ContactRec>,
    Vec<TangentRec>,
    Vec<Vec<(f64, f64)>>,
    std::collections::HashMap<EdgeKey, Vec<Pave>>,
    Vec<Vec<usize>>,
    Vec<Vec<usize>>,
    Vec<Junction>,
)> {
    use ogeom_intersect::{
        CurveCurveOptions, IntersectOptions, SurfaceIntersection, intersect_curves,
        intersect_surfaces,
    };
    let mut sections: Vec<SectionRec> = Vec::new();
    let mut contacts: Vec<ContactRec> = Vec::new();
    let mut tangents: Vec<TangentRec> = Vec::new();
    let mut same_pairs: Vec<(usize, usize)> = Vec::new();
    // Marched sections are fitted; the fit is driven below the confusion
    // tolerance so a fitted curve meets edges, vertices and the mesh welder
    // on the same terms as an exact one. The budget each still carries is
    // recorded per section and widens the crossing filters.
    for (ia, fa) in ga.faces.iter().enumerate() {
        for (ib, fb) in gb.faces.iter().enumerate() {
            let admitted = fa.bound.intersects(&fb.bound);
            if !admitted && !admit_all {
                // The faces cannot meet, whatever their surfaces do.
                continue;
            }
            let scale = fa.chord_scale.min(fb.chord_scale);
            let chord = (scale * 1e-7).max(tol.confusion() * 0.5);
            let options = IntersectOptions {
                // The fit budget scales with the chord: a section against a
                // fitted surface cannot honestly land closer than the
                // geometry it cuts, and whatever it carries is stated on
                // the record and widens every filter downstream.
                tolerance: chord,
                marching: ogeom_intersect::Marching {
                    chord,
                    ..ogeom_intersect::Marching::default()
                },
            };
            // Coincidence is asked before the intersector is, but only where
            // the closed forms have already declined the pair. The marcher
            // documents that it is not the one to answer it (it seeds on
            // sign changes, and a pair that never separates has none), so
            // what it traces over a coincident pair is noise wearing a
            // section's name, and it costs seconds to produce.
            let met = if ogeom_intersect::surface_surface(&fa.surface, &fb.surface, tol).is_err()
                && surfaces_coincide(&fa.surface, &fb.surface, options.tolerance, tol)
            {
                SurfaceIntersection::Same
            } else {
                intersect_surfaces(&fa.surface, &fb.surface, options, tol)?
            };
            match met {
                SurfaceIntersection::Apart => {}
                SurfaceIntersection::Same => {
                    // Coincident surfaces offer no section curve; what splits
                    // each face is the *other* face's boundary. Exact pcurve
                    // projection carries an edge into the other chart, and
                    // planes always have one; a curved same-domain pair whose
                    // edges do not project in closed form is still refused.
                    same_pairs.push((ia, ib));
                    for (owner_from_a, owner, target_from_a, target, target_face) in
                        [(false, fb, true, fa, ia), (true, fa, false, fb, ib)]
                    {
                        let _ = owner_from_a;
                        for e in &owner.edges {
                            let (pcurve, prange) = match ogeom_intersect::exact_pcurve_of(
                                &e.curve,
                                &target.surface,
                                tol,
                            ) {
                                Some(exact) => (exact, e.crange),
                                // A fitted edge has no closed-form projection,
                                // but when the two faces sit on the
                                // *identical chart*, which is exactly the
                                // situation `Same` names for the analytics,
                                // the owner's own stored pcurve already is
                                // the projection, attached at construction,
                                // and it travels with its own window the way
                                // every stored pcurve does. A chart that
                                // merely coincides as a point set is still
                                // refused.
                                None if same_chart(&owner.surface, &target.surface, tol) => {
                                    (e.pcurve.clone(), e.prange)
                                }
                                // Two patches that coincide as point sets
                                // without being one chart (a blend's leg
                                // on a spline host continued past the
                                // face, against the face's own patch) get
                                // the edge fitted by projection into the
                                // target's chart, same-parameter with the
                                // edge, the way every reader derives a
                                // pcurve it was not given.
                                None => {
                                    let Some(pcurve) = projected_into_shared_chart(
                                        &e.curve,
                                        e.crange,
                                        &owner.surface,
                                        &target.surface,
                                        tol,
                                    )?
                                    else {
                                        // Wholly outside the target's window
                                        // is wholly outside the target: a
                                        // blend's run-out past the face it
                                        // melts with splits nothing there.
                                        continue;
                                    };
                                    if *DEBUG_WIRE {
                                        eprintln!(
                                            "SAME owner {} face {} (from_a {owner_from_a}) edge {} of kind {:?}: {:?} .. {:?} against target face {target_face}",
                                            if owner_from_a { "a" } else { "b" },
                                            if owner_from_a { ia } else { ib },
                                            e.node.index(),
                                            core::mem::discriminant(&e.curve),
                                            e.curve.point_at(e.crange.0, tol).ok(),
                                            e.curve.point_at(e.crange.1, tol).ok()
                                        );
                                    }
                                    (pcurve, e.crange)
                                }
                            };
                            contacts.push(ContactRec {
                                curve: e.curve.clone(),
                                crange: e.crange,
                                pcurve,
                                prange,
                                node: e.node,
                                tolerance: e.tolerance,
                                target_from_a,
                                target_face,
                            });
                        }
                    }
                }
                // A touch at a point bounds nothing: no curve to split a
                // face along, no side that is inside on one hand and
                // outside on the other. It contributes to the arrangement
                // exactly what a tangential contact curve does, which is
                // nothing, and the result carries the touch as the
                // non-manifold contact it is.
                SurfaceIntersection::Touching(_) => {}
                SurfaceIntersection::Along(curves) => {
                    for sc in curves {
                        // A tangential curve is contact, not crossing: the
                        // two faces meet along it and neither passes
                        // through the other, so it is set aside from the
                        // arrangement entirely and kept for the consumers
                        // that draw contact rather than classify by it.
                        if sc.tangential {
                            if let (Some(pa), Some(pb)) = (sc.on_a, sc.on_b) {
                                tangents.push(TangentRec {
                                    curve: sc.curve,
                                    pc_a: pa,
                                    pc_b: pb,
                                    face_a: ia,
                                    face_b: ib,
                                });
                            }
                            continue;
                        }
                        // A section can be no longer than a turn round the
                        // faces it cuts: a marched trace that wandered off
                        // beside a chart's pole came back twenty-five times
                        // the circle it stood for, faithfully fitted. That
                        // is no section, and is refused by name.
                        if !sc.exact {
                            use ogeom_geom::Curve3d as _;
                            let (lo, hi) = sc.curve.domain();
                            let mut length = 0.0_f64;
                            let mut last: Option<Point> = None;
                            for k in 0..=64 {
                                let t = if k == 64 {
                                    hi
                                } else {
                                    lo + (hi - lo) * f64::from(k) / 64.0
                                };
                                let p = sc.curve.point_at(t, tol)?;
                                if let Some(q) = last {
                                    length += q.distance(p);
                                }
                                last = Some(p);
                            }
                            let turn = 4.0 * fa.bound.diagonal().max(fb.bound.diagonal());
                            if *DEBUG_WIRE {
                                eprintln!(
                                    "SECTION CHECK faces {ia}/{ib} exact {} closed {} length {length:.4} turn {turn:.4} tol {:.2e}",
                                    sc.exact, sc.closed, sc.tolerance
                                );
                            }
                            if length > turn {
                                ogeom_bail!(
                                    NotDone,
                                    "a marched section of length {length} runs beyond a turn \
                                     round the faces it cuts ({turn}); its trace wandered \
                                     beside a chart's pole; see docs/PARITY.md, bool.booleans"
                                );
                            }
                        }
                        match (sc.on_a.clone(), sc.on_b.clone()) {
                            (Some(pa), Some(pb)) => sections.push(SectionRec {
                                curve: sc.curve,
                                pc_a: pa,
                                pc_b: pb,
                                face_a: ia,
                                face_b: ib,
                                closed: sc.closed,
                                tolerance: sc.tolerance,
                            }),
                            _ => {
                                // A section running through a chart
                                // degeneracy (a plane cutting a ball on its
                                // own axis meets it at both poles) has no
                                // single chart image, because the longitude
                                // jumps half a turn there. Each piece
                                // *between* the poles does have one, and it
                                // is exact. Split first; march only if that
                                // fails.
                                if let Some(split) = split_at_degeneracies(&sc.curve, fa, fb, tol)?
                                {
                                    for (curve, pa, pb) in split {
                                        sections.push(SectionRec {
                                            curve,
                                            pc_a: pa,
                                            pc_b: pb,
                                            face_a: ia,
                                            face_b: ib,
                                            closed: false,
                                            tolerance: 0.0,
                                        });
                                    }
                                    continue;
                                }
                                // An exact curve whose projection has no
                                // closed form keeps its exactness and has
                                // its chart image fitted: the curve's own
                                // points inverted on the surface, read in
                                // the chart, fitted at the curve's own
                                // parameters: a plane's circle passing
                                // beside a sphere chart's pole, whose image
                                // swings fast but is a curve all the same.
                                // Marching the pair, which follows, wandered
                                // beside the pole in both ways there are.
                                let image = |surface: &SurfaceGeometry,
                                             have: Option<&PlanarCurve>|
                                 -> OgeomResult<Option<(PlanarCurve, f64)>> {
                                    if let Some(pc) = have {
                                        return Ok(Some((pc.clone(), 0.0)));
                                    }
                                    fitted_image(&sc.curve, surface, options.tolerance, tol)
                                };
                                if let (Some((pa, ea)), Some((pb, eb))) = (
                                    image(&fa.surface, sc.on_a.as_ref())?,
                                    image(&fb.surface, sc.on_b.as_ref())?,
                                ) {
                                    sections.push(SectionRec {
                                        curve: sc.curve,
                                        pc_a: pa,
                                        pc_b: pb,
                                        face_a: ia,
                                        face_b: ib,
                                        closed: sc.closed,
                                        tolerance: sc.tolerance.max(ea).max(eb),
                                    });
                                    continue;
                                }
                                // No image fits: march the pair instead, so
                                // curve and pcurves are fitted *together*.
                                if *DEBUG_WIRE {
                                    eprintln!(
                                        "MARCH PAIR faces {ia}/{ib} a {:?} b {:?} admitted {admitted}",
                                        core::mem::discriminant(&fa.surface),
                                        core::mem::discriminant(&fb.surface)
                                    );
                                }
                                let shared = if admitted {
                                    fa.bound.intersection(&fb.bound)
                                } else {
                                    // Audit only: disjoint bounds have no
                                    // window, and an empty window would mask
                                    // the very miss being hunted.
                                    fa.bound.union(&fb.bound)
                                };
                                for fitted in march_pair(
                                    &windowed_to(&fa.surface, &shared),
                                    &windowed_to(&fb.surface, &shared),
                                    &options,
                                    tol,
                                )? {
                                    // The same two refusals the intersector's
                                    // own marched sections meet: a fit a
                                    // thousand chords off its trace, and a
                                    // trace longer than a turn round the
                                    // faces. A plane's circle passing beside
                                    // a sphere chart's pole marched both ways
                                    // here (a section six tenths of a
                                    // millimetre off, a section twenty-five
                                    // laps long), and stated as data they
                                    // welded the tool into a point.
                                    let budget =
                                        (options.marching.chord * 1e3).max(options.tolerance * 1e3);
                                    if fitted.fit_error > budget {
                                        ogeom_bail!(
                                            NotDone,
                                            "a marched section's fit misses its trace by {} \
                                             against a chord of {}; a branch passing beside a \
                                             chart's pole fits nothing yet; see docs/PARITY.md, \
                                             boolean.general",
                                            fitted.fit_error,
                                            options.marching.chord
                                        );
                                    }
                                    let curve: Curve = fitted.curve.into();
                                    let length = {
                                        use ogeom_geom::Curve3d as _;
                                        let (lo, hi) = curve.domain();
                                        let mut length = 0.0_f64;
                                        let mut last: Option<Point> = None;
                                        for k in 0..=64 {
                                            let t = if k == 64 {
                                                hi
                                            } else {
                                                lo + (hi - lo) * f64::from(k) / 64.0
                                            };
                                            let p = curve.point_at(t, tol)?;
                                            if let Some(q) = last {
                                                length += q.distance(p);
                                            }
                                            last = Some(p);
                                        }
                                        length
                                    };
                                    let turn = 4.0 * fa.bound.diagonal().max(fb.bound.diagonal());
                                    if length > turn {
                                        ogeom_bail!(
                                            NotDone,
                                            "a marched section of length {length} runs beyond a \
                                             turn round the faces it cuts ({turn}); its trace \
                                             wandered beside a chart's pole; see docs/PARITY.md, \
                                             boolean.general"
                                        );
                                    }
                                    sections.push(SectionRec {
                                        closed: curve.is_closed(tol),
                                        tolerance: options.marching.chord + fitted.fit_error,
                                        curve,
                                        pc_a: fitted.on_a.into(),
                                        pc_b: fitted.on_b.into(),
                                        face_a: ia,
                                        face_b: ib,
                                    });
                                }
                                break;
                            }
                        }
                    }
                }
            }
        }
    }

    // Boundary polylines per face, for the trim tests.
    let outline = |face: &GFace| -> OgeomResult<Vec<Vec<Point2>>> {
        let mut lines = Vec::new();
        for e in &face.edges {
            lines.push(pcurve_polyline(
                &e.pcurve,
                e.prange,
                e.crange,
                e.crange,
                &face.surface,
                tol,
            )?);
            if let Some((other, orange)) = &e.other_side {
                lines.push(pcurve_polyline(
                    other,
                    *orange,
                    e.crange,
                    e.crange,
                    &face.surface,
                    tol,
                )?);
            }
        }
        weld_outline_ends(&mut lines, outline_snap(face, tol));
        Ok(lines)
    };
    let mut outlines_a = Vec::new();
    for f in &ga.faces {
        outlines_a.push(outline(f)?);
    }
    let mut outlines_b = Vec::new();
    for f in &gb.faces {
        outlines_b.push(outline(f)?);
    }

    // Crossings of each section with the boundary edges of both its faces,
    // and with every other section sharing a face.
    let mut paves: std::collections::HashMap<EdgeKey, Vec<Pave>> = std::collections::HashMap::new();
    let mut pieces: Vec<SectionPiece> = Vec::new();
    // Each section's paving depends only on the sections and the two
    // gathered solids, all read-only here, and writes nothing the next
    // section reads. So the measuring runs in parallel and the accumulating
    // runs afterwards in section order: the same split `tessellate` uses,
    // and the same reason: nothing about scheduling can reach the answer.
    type SectionWork = (Vec<(EdgeKey, Pave)>, Vec<SectionPiece>, Vec<Junction>);
    let mut hug_junctions: Vec<Junction> = Vec::new();
    let paved: Vec<OgeomResult<SectionWork>> = ogeom_core::parallel::map_ordered(
        &sections,
        |si, section: &SectionRec| {
            ogeom_core::progress::checkpoint()?;
            let mut paves: Vec<(EdgeKey, Pave)> = Vec::new();
            let mut pieces: Vec<SectionPiece> = Vec::new();
            let mut junctions: Vec<Junction> = Vec::new();
            // A fitted section meets an edge within its own budget, not within
            // rounding.
            let reach = tol.confusion().max(section.tolerance * 2.0);
            let cc = CurveCurveOptions {
                gap: reach.max(CurveCurveOptions::default().gap),
                ..CurveCurveOptions::default()
            };
            let domain = section.curve.domain();
            let mut trim_ts: Vec<f64> = Vec::new();
            // Crossings with boundary edges: side, edge, parameter on the
            // edge, parameter on the section, and how honestly the stop sits.
            let mut hits: Vec<(usize, EdgeKey, f64, f64, f64)> = Vec::new();
            // Spans of the section running *along* a boundary edge. The split
            // such a span would make already exists as boundary (stacked boxes'
            // perpendicular side planes meet exactly at the boxes' own edges),
            // so the span is excluded from that face's strands rather than
            // refused or duplicated.
            //
            // Which face, though, is the whole point. A plane through a bore's
            // axis meets the wall along two rulings, and one of them is the
            // wall's own seam: on the wall that ruling is boundary already, and
            // on the plane it is the curve that separates the two halves the
            // section leaves. Dropped from both, the plane keeps one region where
            // it has two, and the result does not close. So the exclusion is
            // recorded per face.
            let mut along: [Vec<(f64, f64, usize)>; 2] = [Vec::new(), Vec::new()];
            for (side, own) in [
                (0_usize, &ga.faces[section.face_a]),
                (1, &gb.faces[section.face_b]),
            ] {
                // The other face the section must also lie in: a crossing
                // with this edge outside it stops nothing that is kept, and a
                // face bounded by thousands of edges would try every one.
                let across = if side == 0 {
                    &gb.faces[section.face_b]
                } else {
                    &ga.faces[section.face_a]
                };
                for (ei, e) in own.edges.iter().enumerate() {
                    if !admit_all && !e.bound.expanded(reach).intersects(&across.bound) {
                        continue;
                    }
                    // On a plane, an edge within its own radius of the face
                    // (the boundary of a merged group of near-coplanar facets)
                    // stands off the plane the section lies in by as much, and
                    // meets the section only that near; it departs from the
                    // face no other way, so a near miss there is the crossing.
                    // A curved face's fitted rail can pass that near a section
                    // without crossing it, and keeps the section's reach.
                    let edge_reach = if matches!(own.surface, SurfaceGeometry::Plane(_)) {
                        reach + e.tolerance
                    } else {
                        reach
                    };
                    let cc = CurveCurveOptions {
                        gap: cc.gap.max(edge_reach),
                        ..cc
                    };
                    let found = intersect_curves(&section.curve, &e.curve, cc, tol)?;
                    for crossing in &found.crossings {
                        if crossing.gap > edge_reach {
                            continue;
                        }
                        let mut on_b = onto_range(crossing.on_b, &e.curve, e.crange, tol);
                        let mut honesty = honest(section.tolerance, tol);
                        // A crossing at a boundary edge's own end *is* that end's
                        // vertex, exactly. The stop the section keeps must be the
                        // vertex's parameter on the section: the meet of two
                        // curves a fit tolerance apart lands a couple of microns
                        // off, the boundary side keeps its exact vertex, and the
                        // rebuilt wire gapes by the difference.
                        let mut on_a = crossing.on_a;
                        // The window is the fit-slop scale, not the section's
                        // own: the boundary curve may be a fitted intersection
                        // from an earlier boolean carrying a couple of microns
                        // of wobble, and a stop that misses the vertex by that
                        // much gapes the rebuilt wire by the same. A fitted
                        // section's stop misses by its own honesty once more:
                        // the trace is off the true crossing by that, and
                        // the crossing found against the edge by that again.
                        let weld =
                            (reach + honest(section.tolerance, tol)).max(tol.confusion() * 1e2);
                        let mut at_end: Option<(f64, Point)> = None;
                        for end in [e.crange.0, e.crange.1] {
                            let vertex = e.curve.point_at(end, tol)?;
                            if vertex.distance(crossing.point) <= weld + crossing.reach {
                                at_end = Some((end, vertex));
                                let snapped =
                                    ogeom_algo::project_on_curve(&section.curve, vertex, 64, tol)?;
                                if snapped.distance <= weld + crossing.reach {
                                    on_a = snapped.parameter;
                                    // And the edge's own parameter is the
                                    // end's: two rim circles through one
                                    // block corner meet at a shallow angle
                                    // there, the crossing lands a few
                                    // hundredths of a millimetre along the
                                    // rim from the corner, and paved there
                                    // it split the rim into a sliver.
                                    on_b = end;
                                }
                                break;
                            }
                        }
                        // A touch at the section's own end *is* the end. A
                        // band's section through a wall meets the wall's drum
                        // edge tangentially where the band is tangent to the
                        // drum, and the touch of two fitted curves at a
                        // shallow angle wanders along them by far more than
                        // their honesty; but the section stops there because
                        // it leaves its own face, and that stop is the
                        // junction: the crossing takes the end's parameter
                        // and the edge splits under the end itself.
                        let touch = crossing.reach > 0.0
                            || tangential(
                                &section.curve,
                                crossing.on_a,
                                &e.curve,
                                crossing.on_b,
                                tol,
                            )?;
                        if touch {
                            let half = (domain.1 - domain.0) * 0.5;
                            for end in [domain.0, domain.1] {
                                if (crossing.on_a - end).abs() > half {
                                    continue;
                                }
                                let tip = section
                                    .curve
                                    .point_at(at_param(end, domain, section.closed), tol)?;
                                if crossing.reach > 0.0
                                    && tip.distance(crossing.point) > crossing.reach + weld
                                {
                                    continue;
                                }
                                // The tip stands where the marcher left the
                                // other face (on a fitted rail, to the rail's
                                // own honesty), so the edge is asked at the
                                // hug width, not the section's.
                                let foot = ogeom_algo::project_on_curve(&e.curve, tip, 64, tol)?;
                                let width =
                                    (crossing.reach + honest(weld, tol)).max(tol.confusion() * 1e3);
                                if *DEBUG_WIRE {
                                    eprintln!(
                                        "PAVE s{si}: touch on edge {} at {:.6}: tip {:.3e} off the edge, width {width:.2e}",
                                        e.node.index(),
                                        crossing.on_a,
                                        foot.distance
                                    );
                                }
                                if foot.distance <= width {
                                    on_a = end;
                                    on_b = onto_range(foot.parameter, &e.curve, e.crange, tol);
                                    // A touch at the edge's own end vertex
                                    // is that vertex: a section tangent to
                                    // the edge there runs along it for a
                                    // stretch either side, stops somewhere
                                    // on that stretch, and split where it
                                    // stopped the edge keeps a sliver the
                                    // face across the vertex never has.
                                    if let Some((end_b, vertex)) = at_end
                                        && tip.distance(vertex) <= width.min(tol.confusion() * 1e5)
                                    {
                                        on_b = end_b;
                                        let gap = tip.distance(vertex);
                                        if gap > tol.confusion() * 1e2 {
                                            junctions.push(Junction {
                                                at: vertex,
                                                reach: gap + tol.confusion() * 1e2,
                                                onto_vertex: true,
                                            });
                                        }
                                    }
                                    // Along the edge the touch is known only
                                    // as far as the two curves stay together:
                                    // a rail grazing a bore is met there by
                                    // the sections on the faces either side
                                    // of it a hundredth of a millimetre
                                    // apart, and those are one junction.
                                    honesty = honesty.max(foot.distance).max(crossing.reach);
                                }
                                break;
                            }
                        }
                        if *DEBUG_WIRE {
                            eprintln!(
                                "PAVE s{si}: side {side} edge {} crossing at {on_a:.6} (edge {on_b:.6}) gap {:.2e} reach {:.2e} honesty {honesty:.2e} at {:?}",
                                e.node.index(),
                                crossing.gap,
                                crossing.reach,
                                crossing.point
                            );
                        }
                        hits.push((side, e.node, on_b, on_a, honesty));
                    }
                    for overlap in &found.overlaps {
                        // The curves overlap; what is *boundary* is the stretch
                        // the edge actually covers. A sphere's seam and the far
                        // half of the same great circle lie on one curve, and
                        // reading the whole curve as boundary makes the meridian
                        // opposite the seam disappear, which is the octant's own
                        // edge, on a ball cut at its corner.
                        let Some((lo, hi)) = overlap_within(overlap, e.crange, &e.curve, tol)
                        else {
                            continue;
                        };
                        if *DEBUG_WIRE {
                            eprintln!("PAVE s{si}: side {side} overlap ({lo:.6}, {hi:.6})");
                        }
                        // The edge splits where the shared stretch ends, as
                        // the section does: the stretch itself is the edge's
                        // to carry, and the section's next piece must meet
                        // the edge at a vertex the edge actually has. A band
                        // meeting a wall tangentially to the wall's own top
                        // edge hugs it for a few microns from their corner;
                        // without the split at the hug's far end the section
                        // dangles there and is pruned, and the wall never
                        // splits.
                        // The section leaves the edge at the hug's end by
                        // what the hug allowed, which no single strand's
                        // honesty covers: that junction owns the gap.
                        for t in [lo, hi] {
                            let at = section
                                .curve
                                .point_at(at_param(t, domain, section.closed), tol)?;
                            let foot = ogeom_algo::project_on_curve(&e.curve, at, 64, tol)?;
                            if foot.distance <= reach.max(tol.confusion() * 1e3) {
                                let on_b = onto_range(foot.parameter, &e.curve, e.crange, tol);
                                hits.push((
                                    side,
                                    e.node,
                                    on_b,
                                    t,
                                    honest(section.tolerance, tol).max(foot.distance),
                                ));
                                if foot.distance > tol.confusion() * 1e2 {
                                    junctions.push(Junction {
                                        at: foot.point,
                                        reach: foot.distance + tol.confusion() * 1e2,
                                        onto_vertex: false,
                                    });
                                }
                            }
                        }
                        trim_ts.push(lo);
                        trim_ts.push(hi);
                        along[side].push((lo, hi, ei));
                    }
                }
            }
            // Inside a stretch the section runs along one of a face's edges,
            // that face's other edges cannot genuinely cross it: a boundary
            // is a simple loop, and they meet the hugged edge only at its
            // ends, which the stretch's own ends already stop at. What the
            // sampler reports there is the tangential dust of a leg touching
            // the arc it ends on: thirty near-crossings inside a micron,
            // which read as stops would shatter the section and pave the
            // leg at each.
            hits.retain(|(side, _, _, on_a, _)| {
                !along[*side].iter().any(|(lo, hi, _)| {
                    *on_a > lo + tol.parametric() && *on_a < hi - tol.parametric()
                })
            });
            let edge_hits: Vec<(EdgeKey, f64, f64, f64)> = hits
                .iter()
                .map(|(_, node, on_b, on_a, honesty)| (*node, *on_b, *on_a, *honesty))
                .collect();
            trim_ts.extend(hits.iter().map(|(_, _, _, on_a, _)| *on_a));
            let mut cross_ts: Vec<f64> = Vec::new();
            for (sj, other) in sections.iter().enumerate() {
                if sj == si {
                    continue;
                }
                if other.face_a != section.face_a && other.face_b != section.face_b {
                    continue;
                }
                let both = reach.max(tol.confusion().max(other.tolerance * 2.0));
                // Where two sections share one face, a crossing matters only
                // inside the faces they do not share, since a section is kept
                // only where it lies in both its faces: one wall met by a
                // thousand facets carries a thousand sections, and only those
                // whose facets' bounds meet can cross where it counts. The
                // audit's unfiltered pass tries every pair still.
                if !admit_all {
                    let apart = |x: &ogeom_math::Aabb, y: &ogeom_math::Aabb| {
                        !x.expanded(both).intersects(&y.expanded(both))
                    };
                    let unshared_a = other.face_a != section.face_a
                        && apart(
                            &ga.faces[section.face_a].bound,
                            &ga.faces[other.face_a].bound,
                        );
                    let unshared_b = other.face_b != section.face_b
                        && apart(
                            &gb.faces[section.face_b].bound,
                            &gb.faces[other.face_b].bound,
                        );
                    if unshared_a || unshared_b {
                        continue;
                    }
                }
                let cc2 = CurveCurveOptions {
                    gap: both.max(CurveCurveOptions::default().gap),
                    ..CurveCurveOptions::default()
                };
                let found = intersect_curves(&section.curve, &other.curve, cc2, tol)?;
                for crossing in &found.crossings {
                    if crossing.gap > both {
                        continue;
                    }
                    let mut at = crossing.on_a;
                    // A touch at this section's end is the end, as against an
                    // edge: two sections through one wall from a band and
                    // the leg it is tangent to touch where the band's ends.
                    let touch = crossing.reach > 0.0
                        || tangential(
                            &section.curve,
                            crossing.on_a,
                            &other.curve,
                            crossing.on_b,
                            tol,
                        )?;
                    if touch {
                        let half = (domain.1 - domain.0) * 0.5;
                        for end in [domain.0, domain.1] {
                            if (crossing.on_a - end).abs() > half {
                                continue;
                            }
                            let tip = section
                                .curve
                                .point_at(at_param(end, domain, section.closed), tol)?;
                            let foot = ogeom_algo::project_on_curve(&other.curve, tip, 64, tol)?;
                            let width =
                                (crossing.reach + honest(both, tol)).max(tol.confusion() * 1e3);
                            if foot.distance <= width {
                                at = end;
                            }
                            break;
                        }
                    }
                    // Two sections through one face that lie on two different
                    // faces of the other solid cross only on the edge those
                    // two share, where each already stops. Meeting there
                    // tangentially (a wall through a plane and the fillet
                    // tangent to it), their own crossing is ill-conditioned
                    // and lands off the edge, as far off as the two stay
                    // within their doubt of each other; the edge's stop is
                    // the one where they cannot be told apart up to it.
                    let shared = |edges: &[BoundaryEdge], node: EdgeKey| {
                        edges.iter().any(|e| e.node == node)
                    };
                    let apart_on = if other.face_a == section.face_a {
                        Some(&gb.faces[other.face_b].edges)
                    } else if other.face_b == section.face_b {
                        Some(&ga.faces[other.face_a].edges)
                    } else {
                        None
                    };
                    if let Some(edges) = apart_on {
                        let close = both.max(tol.confusion() * 1e2);
                        let mut best: Option<(f64, f64)> = None;
                        for (node, _, on_a, _) in &edge_hits {
                            if !shared(edges, *node)
                                || best.is_some_and(|(bd, _)| (on_a - at).abs() >= bd)
                            {
                                continue;
                            }
                            let mut together = true;
                            for k in 1..8 {
                                let t = (on_a - at).mul_add(f64::from(k) / 8.0, at);
                                let p = section
                                    .curve
                                    .point_at(at_param(t, domain, section.closed), tol)?;
                                if ogeom_algo::project_on_curve(&other.curve, p, 16, tol)?.distance
                                    > close
                                {
                                    together = false;
                                    break;
                                }
                            }
                            if together {
                                best = Some(((on_a - at).abs(), *on_a));
                            }
                        }
                        if let Some((_, on_a)) = best {
                            at = on_a;
                        }
                    }
                    if *DEBUG_WIRE {
                        eprintln!(
                            "PAVE s{si}: cross s{sj} at {at:.6} gap {:.2e}",
                            crossing.gap
                        );
                    }
                    cross_ts.push(at);
                }
            }

            // Candidate intervals between trim crossings, kept where the middle
            // sits inside both faces' trims.
            trim_ts.sort_by(|a, b| a.partial_cmp(b).unwrap_or(core::cmp::Ordering::Equal));
            trim_ts.dedup_by(|a, b| (*a - *b).abs() <= tol.parametric());
            // Stops the section cannot tell apart are one stop: two
            // crossings a few nanometres apart along it would make an
            // interval of dust that one face keeps and the other drops.
            {
                let floor = (honest(section.tolerance, tol) * 3.0).max(tol.confusion() * 10.0);
                let mut kept: Vec<f64> = Vec::with_capacity(trim_ts.len());
                let mut held: Option<Point> = None;
                for t in &trim_ts {
                    let at = section
                        .curve
                        .point_at(at_param(*t, domain, section.closed), tol)?;
                    if held.is_none_or(|h: Point| h.distance(at) > floor) {
                        kept.push(*t);
                        held = Some(at);
                    }
                }
                trim_ts = kept;
            }
            let mut candidates: Vec<(f64, f64)> = Vec::new();
            if section.closed {
                let period = domain.1 - domain.0;
                if trim_ts.is_empty() {
                    candidates.push(domain);
                } else {
                    for i in 0..trim_ts.len() {
                        let lo = trim_ts[i];
                        let hi = if i + 1 < trim_ts.len() {
                            trim_ts[i + 1]
                        } else {
                            trim_ts[0] + period
                        };
                        candidates.push((lo, hi));
                    }
                }
            } else {
                let mut stops = vec![domain.0];
                stops.extend(
                    trim_ts
                        .iter()
                        .copied()
                        .filter(|t| *t > domain.0 && *t < domain.1),
                );
                stops.push(domain.1);
                for pair in stops.windows(2) {
                    candidates.push((pair[0], pair[1]));
                }
            }

            // Whether the section at `t` lies inside each face's trim.
            let inside_each = |t: f64| -> OgeomResult<[bool; 2]> {
                let tf = if section.closed { fold(t, domain) } else { t };
                let la: Vec<&[Point2]> = outlines_a[section.face_a]
                    .iter()
                    .map(Vec::as_slice)
                    .collect();
                let lb: Vec<&[Point2]> = outlines_b[section.face_b]
                    .iter()
                    .map(Vec::as_slice)
                    .collect();
                let ua = fold_inside(
                    section.pc_a.point_at(tf, tol)?,
                    &ga.faces[section.face_a].surface,
                    &la,
                );
                let ub = fold_inside(
                    section.pc_b.point_at(tf, tol)?,
                    &gb.faces[section.face_b].surface,
                    &lb,
                );
                Ok([inside_many(&la, ua), inside_many(&lb, ub)])
            };

            for (lo, hi) in candidates {
                if hi - lo <= tol.parametric() {
                    continue;
                }
                let mid = f64::midpoint(lo, hi);
                let mid_folded = if section.closed {
                    fold(mid, domain)
                } else {
                    mid
                };
                let held = inside_each(mid)?;
                if *DEBUG_WIRE {
                    let tf = if section.closed {
                        fold(mid, domain)
                    } else {
                        mid
                    };
                    eprintln!(
                        "PAVE s{si}: candidate ({lo:.6}, {hi:.6}) inside {held:?} at a {:?} b {:?}",
                        section
                            .pc_a
                            .point_at(tf, tol)
                            .ok()
                            .map(|q| fold_point_into_chart(q, &ga.faces[section.face_a].surface)),
                        section
                            .pc_b
                            .point_at(tf, tol)
                            .ok()
                            .map(|q| fold_point_into_chart(q, &gb.faces[section.face_b].surface))
                    );
                }
                // A section that runs along a boundary edge of a face splits
                // nothing *there*: the split already exists as boundary. The
                // analytic overlap detection above catches the same-support
                // cases; this catches the rest (a fitted section tracing a
                // boundary curve, a surface meeting another exactly at its own
                // trim) by measurement rather than by recognising supports.
                // The hug is asked before the trim: a section on a face's own
                // edge reads inside or outside that face by a hair, and the
                // corner tool's block face meets a band exactly along the arc
                // that bounds it. On the hugging face it is boundary already;
                // what matters is whether the *other* face holds it.
                let mut hugs = [false; 2];
                // The edges each side hugs (a rim already in pieces is
                // several): they must split wherever the section does, and
                // the section wherever they end, or the two faces walk
                // different subdivisions of one curve.
                let mut hugged: [Vec<usize>; 2] = [Vec::new(), Vec::new()];
                for (side, own, side_from_a, side_face) in [
                    (0_usize, &ga.faces[section.face_a], true, section.face_a),
                    (1, &gb.faces[section.face_b], false, section.face_b),
                ] {
                    if let Some((_, _, ei)) = along[side]
                        .iter()
                        .find(|(alo, ahi, _)| mid_folded >= *alo && mid_folded <= *ahi)
                    {
                        hugs[side] = true;
                        hugged[side].push(*ei);
                        continue;
                    }
                    let mut all_near = true;
                    let mut votes: Vec<usize> = vec![0; own.edges.len()];
                    let mut contact_votes: Vec<usize> = vec![0; contacts.len()];
                    for i in 0..=4 {
                        let t = lo + (hi - lo) * f64::from(i) / 4.0;
                        let tf = if section.closed { fold(t, domain) } else { t };
                        let at = section.curve.point_at(tf, tol)?;
                        // Wider than the crossing filters on purpose for a
                        // fitted section: a tangentially-traced curve wobbles
                        // about the boundary it hugs by far more than a fit
                        // budget, and a genuine section keeps a distance of
                        // feature scale, not microns. An exact section keeps
                        // to the edge's own honesty: a plane's ellipse across
                        // a band touches the band's rail tangentially, and
                        // the stretch within a tenth of a millimetre of the
                        // rail is the section, not the rail; read as a hug
                        // it left the cap's crescent unable to close.
                        let floor = if section.tolerance > 0.0 {
                            tol.confusion() * 1e3
                        } else {
                            tol.confusion() * 10.0
                        };
                        let mut near = false;
                        for (ei, e) in own.edges.iter().enumerate() {
                            let width = reach.max(floor).max(e.tolerance * 2.0);
                            let d = distance_to_edge_curve(&e.curve, e.crange, at, tol)?;
                            if d <= width {
                                if *DEBUG_WIRE {
                                    eprintln!(
                                        "PAVE s{si}: side {side} sample {i} hugs edge {} at {d:.2e} (width {width:.2e}) range {:?} point {at:?}",
                                        e.node.index(),
                                        e.crange
                                    );
                                }
                                // Every edge within reach is hugged, not the
                                // first found: where three rims converge on
                                // a pole corner all three come within a
                                // micron of the section, and the one it
                                // actually runs along must not lose its
                                // vote to a neighbour listed before it.
                                near = true;
                                votes[ei] += 1;
                            }
                        }
                        let width = reach.max(floor);
                        // A contact edge is a strand on this face too (the other
                        // solid's boundary, carried into a chart they share), so a
                        // section tracing one would be the same curve twice, and
                        // the arrangement cannot walk a line it meets from both
                        // sides at once. Two boxes side by side put the low one's
                        // lid exactly there.
                        for (ci, c) in contacts.iter().enumerate() {
                            if c.target_from_a != side_from_a || c.target_face != side_face {
                                continue;
                            }
                            if distance_to_edge_curve(&c.curve, c.crange, at, tol)?
                                <= width.max(c.tolerance * 2.0)
                            {
                                near = true;
                                contact_votes[ci] += 1;
                            }
                        }
                        if !near {
                            all_near = false;
                            break;
                        }
                    }
                    // And along one of them all the way: a section winding
                    // round a drum passes the drum's seam once a turn, and a
                    // stretch spanning whole turns puts every sample on the
                    // seam without running along it anywhere between.
                    let one_line = votes.iter().chain(&contact_votes).any(|&v| v == 5);
                    hugs[side] = all_near && one_line;
                    if all_near {
                        hugged[side].extend(
                            votes
                                .iter()
                                .enumerate()
                                .filter(|(_, c)| **c > 0)
                                .map(|(ei, _)| ei),
                        );
                    }
                }
                if *DEBUG_WIRE {
                    eprintln!("PAVE s{si}: candidate ({lo:.6}, {hi:.6}) hugs {hugs:?}");
                }
                // Inside both faces the piece is a section of both. Hugging
                // one face's own edge while that face's trim reads it out
                // (by a hair, on the edge), it is still the other face's
                // split, and is admitted there, unless that edge is already
                // carried onto the other face as a contact: then the split
                // is laid down once already, and a section would lay it
                // twice.
                let mut hug_key: Option<(EdgeKey, bool, usize)> = None;
                let admitted = if held[0] && held[1] {
                    true
                } else if hugs[0] != hugs[1] {
                    let hugging = usize::from(hugs[1]);
                    let other = 1 - hugging;
                    let (own, other_from_a, other_face) = if hugging == 0 {
                        (&ga.faces[section.face_a], false, section.face_b)
                    } else {
                        (&gb.faces[section.face_b], true, section.face_a)
                    };
                    let carried = hugged[hugging].iter().any(|&ei| {
                        contacts.iter().any(|c| {
                            c.node == own.edges[ei].node
                                && c.target_from_a == other_from_a
                                && c.target_face == other_face
                        })
                    });
                    if held[other] && !carried {
                        hug_key = hugged[hugging]
                            .first()
                            .map(|&ei| (own.edges[ei].node, other_from_a, other_face));
                        true
                    } else {
                        false
                    }
                } else {
                    false
                };
                if !admitted {
                    continue;
                }
                if hugs[0] && hugs[1] {
                    // Boundary on both sides: the split exists twice over and
                    // adding it a third time would cancel what it copies.
                    continue;
                }
                // Keep the paves that end a kept interval: those are where edges
                // genuinely split.
                for (node, on_edge, on_section, honesty) in &edge_hits {
                    let s = *on_section;
                    let near = |x: f64| {
                        (s - x).abs() <= tol.parametric()
                            || (section.closed
                                && ((s + (domain.1 - domain.0)) - x).abs() <= tol.parametric())
                    };
                    if near(lo) || near(hi) {
                        paves.push((
                            *node,
                            Pave {
                                t: *on_edge,
                                honesty: *honesty,
                            },
                        ));
                    }
                }
                // Split at section/section crossings inside the kept interval,
                // so every face sees the same subdivision.
                let mut cuts = vec![lo];
                // A wrap interval also splits at the curve's own domain end, so
                // every piece lives within one period and evaluates in-domain
                // after a single fold.
                if section.closed
                    && domain.1 > lo + tol.parametric()
                    && domain.1 < hi - tol.parametric()
                {
                    cuts.push(domain.1);
                }
                for &c in &cross_ts {
                    let c2 = if section.closed && c < lo {
                        c + (domain.1 - domain.0)
                    } else {
                        c
                    };
                    if c2 > lo + tol.parametric() && c2 < hi - tol.parametric() {
                        cuts.push(c2);
                    }
                }
                cuts.push(hi);
                // A hugged edge's own ends are cuts of the section too: the
                // edge may already be several pieces (a sphere's rim split
                // at its chart seam), and the section must walk the same
                // pieces, or the two faces never sew along it.
                for side in 0..2 {
                    if !hugs[side] {
                        continue;
                    }
                    let own = if side == 0 {
                        &ga.faces[section.face_a]
                    } else {
                        &gb.faces[section.face_b]
                    };
                    for &ei in &hugged[side] {
                        let e = &own.edges[ei];
                        for end in [e.crange.0, e.crange.1] {
                            let p = e.curve.point_at(end, tol)?;
                            let foot = ogeom_algo::project_on_curve(&section.curve, p, 64, tol)?;
                            if *DEBUG_WIRE {
                                eprintln!(
                                    "PAVE s{si}: side {side} hugged edge {} end {end:.6} lands at {:.6} off {:.2e} against ({lo:.6}, {hi:.6})",
                                    e.node.index(),
                                    foot.parameter,
                                    foot.distance
                                );
                            }
                            if foot.distance > reach.max(tol.confusion() * 1e3) {
                                continue;
                            }
                            let mut t = foot.parameter;
                            if section.closed {
                                let period = domain.1 - domain.0;
                                while t < lo {
                                    t += period;
                                }
                                while t - period >= lo {
                                    t -= period;
                                }
                            }
                            if t > lo + tol.parametric() && t < hi - tol.parametric() {
                                if *DEBUG_WIRE {
                                    eprintln!(
                                        "PAVE s{si}: side {side} cut at {t:.6} where hugged edge {} ends",
                                        e.node.index()
                                    );
                                }
                                cuts.push(t);
                            }
                        }
                    }
                }
                cuts.sort_by(|a, b| a.partial_cmp(b).unwrap_or(core::cmp::Ordering::Equal));
                // Cuts the section cannot tell apart are one cut: two
                // crossings a few nanometres apart along it would leave a
                // piece of dust one face keeps and the other drops. The
                // interval's own ends stay.
                {
                    let floor = (honest(section.tolerance, tol) * 3.0).max(tol.confusion() * 10.0);
                    let mut kept: Vec<f64> = Vec::with_capacity(cuts.len());
                    let mut held: Option<Point> = None;
                    for (index, c) in cuts.iter().enumerate() {
                        let at = section
                            .curve
                            .point_at(at_param(*c, domain, section.closed), tol)?;
                        let last = index + 1 == cuts.len();
                        if index == 0 || last || held.is_none_or(|h: Point| h.distance(at) > floor)
                        {
                            kept.push(*c);
                            held = Some(at);
                        }
                    }
                    cuts = kept;
                }
                // A hugged edge is boundary where this piece is section: it
                // must split at the piece's interior cuts too, or the face
                // that keeps the edge walks one long strand where its
                // neighbour walks three, and the sew pairs none of them.
                for side in 0..2 {
                    if !hugs[side] {
                        continue;
                    }
                    let own = if side == 0 {
                        &ga.faces[section.face_a]
                    } else {
                        &gb.faces[section.face_b]
                    };
                    for &ei in &hugged[side] {
                        let e = &own.edges[ei];
                        // The piece's own ends as well as its interior cuts: a
                        // stop against another face's edge is a split of the
                        // hugged edge too, or the section's neighbour walks two
                        // strands where the hugged edge's face walks one.
                        for c in &cuts {
                            let at = section
                                .curve
                                .point_at(at_param(*c, domain, section.closed), tol)?;
                            let foot = ogeom_algo::project_on_curve(&e.curve, at, 64, tol)?;
                            // A cut off the edge (past the stretch the
                            // section hugs, on a neighbouring edge) is no
                            // split of it; stated as one with its distance
                            // for honesty it seeds a junction a feature wide.
                            if foot.distance > reach.max(tol.confusion() * 1e3) {
                                continue;
                            }
                            let on_e = onto_range(foot.parameter, &e.curve, e.crange, tol);
                            if on_e <= e.crange.0.min(e.crange.1) + tol.parametric()
                                || on_e >= e.crange.0.max(e.crange.1) - tol.parametric()
                            {
                                continue;
                            }
                            if *DEBUG_WIRE {
                                eprintln!(
                                    "PAVE s{si}: side {side} hugged edge {} paved at {on_e:.6} for the cut at {c:.6}",
                                    e.node.index()
                                );
                            }
                            paves.push((
                                e.node,
                                Pave {
                                    t: on_e,
                                    honesty: honest(section.tolerance, tol).max(foot.distance),
                                },
                            ));
                        }
                    }
                }
                for pair in cuts.windows(2) {
                    let (lo2, hi2) = (pair[0], pair[1]);
                    let from = section
                        .curve
                        .point_at(at_param(lo2, domain, section.closed), tol)?;
                    let to = section
                        .curve
                        .point_at(at_param(hi2, domain, section.closed), tol)?;
                    if from.distance(to) <= tol.confusion() {
                        // A full loop: two arcs, so every strand has two
                        // distinct endpoints.
                        let mid = f64::midpoint(lo2, hi2);
                        pieces.push(SectionPiece {
                            section: si,
                            range: (lo2, mid),
                            hugs,
                            hug_key,
                        });
                        pieces.push(SectionPiece {
                            section: si,
                            range: (mid, hi2),
                            hugs,
                            hug_key,
                        });
                    } else {
                        pieces.push(SectionPiece {
                            section: si,
                            range: (lo2, hi2),
                            hugs,
                            hug_key,
                        });
                    }
                }
            }
            Ok((paves, pieces, junctions))
        },
    );
    for work in paved {
        let (found, made, hugged) = work?;
        hug_junctions.extend(hugged);
        for (node, at) in found {
            paves.entry(node).or_default().push(at);
        }
        pieces.extend(made);
    }
    // A piece running along a face's own edge and that edge are one curve,
    // split twice: the piece where its section's cuts fell, the edge where
    // every section that met it paved it, each worked out apart. A loop
    // hugging a rim is split at its middle, the rim wherever another
    // section crossed it; the two faces then walk different pieces of one
    // circle and never sew. Now that every pave is in, each such piece is
    // cut wherever its edge is split, and the edge split wherever the piece
    // ends, twice over so a split one piece adds reaches the others.
    for _ in 0..2 {
        let mut settled: Vec<SectionPiece> = Vec::with_capacity(pieces.len());
        for piece in pieces.drain(..) {
            if !(piece.hugs[0] || piece.hugs[1]) {
                settled.push(piece);
                continue;
            }
            let section = &sections[piece.section];
            let domain = section.curve.domain();
            let period = domain.1 - domain.0;
            let (lo, hi) = piece.range;
            let at = |t: f64| {
                section
                    .curve
                    .point_at(at_param(t, domain, section.closed), tol)
            };
            let width = (tol.confusion() * 1e3).max(section.tolerance * 3.0);
            let mid = at(f64::midpoint(lo, hi))?;
            let mut hugged: Vec<&BoundaryEdge> = Vec::new();
            for side in 0..2 {
                if !piece.hugs[side] {
                    continue;
                }
                let face = if side == 0 {
                    &ga.faces[section.face_a]
                } else {
                    &gb.faces[section.face_b]
                };
                for e in &face.edges {
                    if distance_to_edge_curve(&e.curve, e.crange, mid, tol)?
                        <= width.max(e.tolerance * 2.0)
                    {
                        hugged.push(e);
                    }
                }
            }
            let mut cuts: Vec<f64> = Vec::new();
            for e in &hugged {
                let mut on_edge: Vec<f64> = vec![e.crange.0, e.crange.1];
                if let Some(list) = paves.get(&e.node) {
                    on_edge.extend(list.iter().map(|pave| pave.t));
                }
                for t in on_edge {
                    let q = e.curve.point_at(t, tol)?;
                    let foot = ogeom_algo::project_on_curve(&section.curve, q, 64, tol)?;
                    if foot.distance > width.max(e.tolerance * 2.0) {
                        continue;
                    }
                    let mut f = foot.parameter;
                    if section.closed && period > 0.0 {
                        while f < lo {
                            f += period;
                        }
                        while f - period >= lo {
                            f -= period;
                        }
                    }
                    if f > lo + tol.parametric() && f < hi - tol.parametric() {
                        cuts.push(f);
                    }
                }
            }
            cuts.push(lo);
            cuts.push(hi);
            cuts.sort_by(f64::total_cmp);
            cuts.dedup_by(|a, b| (*a - *b).abs() <= tol.parametric());
            // The edge split wherever the piece now ends.
            for e in &hugged {
                for &c in &cuts {
                    let q = at(c)?;
                    let foot = ogeom_algo::project_on_curve(&e.curve, q, 64, tol)?;
                    if foot.distance > width.max(e.tolerance * 2.0) {
                        continue;
                    }
                    let on_e = onto_range(foot.parameter, &e.curve, e.crange, tol);
                    if on_e <= e.crange.0.min(e.crange.1) + tol.parametric()
                        || on_e >= e.crange.0.max(e.crange.1) - tol.parametric()
                    {
                        continue;
                    }
                    let list = paves.entry(e.node).or_default();
                    if list
                        .iter()
                        .all(|pave| (pave.t - on_e).abs() > tol.parametric())
                    {
                        list.push(Pave {
                            t: on_e,
                            honesty: honest(section.tolerance, tol).max(foot.distance),
                        });
                    }
                }
            }
            for pair in cuts.windows(2) {
                settled.push(SectionPiece {
                    section: piece.section,
                    range: (pair[0], pair[1]),
                    hugs: piece.hugs,
                    hug_key: piece.hug_key,
                });
            }
        }
        pieces = settled;
    }
    // Two sections hugging onto one face along one line lay the same split
    // down twice (a top face's trace and its band's, both along the rail
    // they share, which need not be one edge node once a shell has been
    // built), and a line met from both sides at once breaks the
    // arrangement. Pieces of different sections admitted by a hug onto the
    // same face, each lying on the other's curve, are one split, and the
    // section covering the most of it speaks for it.
    {
        let width = tol.confusion() * 1e3;
        let mut drop: Vec<usize> = Vec::new();
        let hugged_onto: Vec<usize> = pieces
            .iter()
            .enumerate()
            .filter(|(_, p)| p.hug_key.is_some())
            .map(|(i, _)| i)
            .collect();
        let mid_of = |p: &SectionPiece| -> OgeomResult<Point> {
            let section = &sections[p.section];
            let domain = section.curve.domain();
            section.curve.point_at(
                at_param(f64::midpoint(p.range.0, p.range.1), domain, section.closed),
                tol,
            )
        };
        let length_of = |p: &SectionPiece| -> OgeomResult<f64> {
            let section = &sections[p.section];
            let domain = section.curve.domain();
            let a = section
                .curve
                .point_at(at_param(p.range.0, domain, section.closed), tol)?;
            let b = section
                .curve
                .point_at(at_param(p.range.1, domain, section.closed), tol)?;
            Ok(a.distance(b))
        };
        // A hug-admitted piece against every other piece kept on the same
        // face: a section admitted the ordinary way, inside both faces,
        // already speaks for the split, and the hug-admitted one yields
        // to it; two hug-admitted ones yield to the longer.
        for &x in &hugged_onto {
            let px = &pieces[x];
            let Some((_, target_from_a, target_face)) = px.hug_key else {
                continue;
            };
            for (y, py) in pieces.iter().enumerate() {
                if y == x || py.section == px.section || drop.contains(&x) || drop.contains(&y) {
                    continue;
                }
                let sy = &sections[py.section];
                let on_target = if target_from_a {
                    sy.face_a == target_face && !py.hugs[0]
                } else {
                    sy.face_b == target_face && !py.hugs[1]
                };
                if !on_target {
                    continue;
                }
                let sx = &sections[px.section];
                // A closed section's pieces are ranged past its domain end
                // and folded on use; the distance helper samples the range
                // it is given, so it is given the folded one.
                let stretch = |section: &SectionRec, range: (f64, f64)| -> (f64, f64) {
                    folded_range(range, section.curve.domain(), section.closed)
                };
                let on_y =
                    distance_to_edge_curve(&sy.curve, stretch(sy, py.range), mid_of(px)?, tol)?
                        <= width;
                let on_x =
                    distance_to_edge_curve(&sx.curve, stretch(sx, px.range), mid_of(py)?, tol)?
                        <= width;
                if !(on_x && on_y) {
                    continue;
                }
                let loser = if py.hug_key.is_some() && length_of(px)? > length_of(py)? {
                    y
                } else {
                    x
                };
                if *DEBUG_WIRE {
                    eprintln!(
                        "PAVE: s{} and s{} split one face along one line; s{} speaks for it",
                        px.section,
                        py.section,
                        pieces[if loser == x { y } else { x }].section
                    );
                }
                drop.push(loser);
                if loser == x {
                    break;
                }
            }
        }
        let mut index = 0_usize;
        pieces.retain(|_| {
            let keep = !drop.contains(&index);
            index += 1;
            keep
        });
    }
    // The audit's verdict. A dropped pair whose surfaces intersect is not a
    // filter bug: planes meet along an infinite line the paving then trims
    // to the faces, usually to nothing. A dropped pair whose section
    // *survives paving* is: some of that curve lies inside both faces, so
    // the faces genuinely meet and the filter's boxes failed to. Contacts
    // need no check: a contact is an owner edge lying in the target face,
    // which forces the boxes to overlap where the edge does.
    //
    // This membership check is the audit's first tooth; the second is the
    // strict replay in general_fuse, which runs the *filtered* fill for the
    // production result and diffs this unfiltered one against it, so the
    // non-compositionality of paving (a section's kept intervals see the
    // other sections' paves) is caught rather than stated as a limit.
    if admit_all {
        for piece in &pieces {
            let section = &sections[piece.section];
            let (fa, fb) = (&ga.faces[section.face_a], &gb.faces[section.face_b]);
            assert!(
                fa.bound.intersects(&fb.bound),
                "bound filter audit: faces {}/{} were dropped by the bound \
                 filter, yet their section paved a surviving piece over \
                 {:?}; the filter under-approximates",
                section.face_a,
                section.face_b,
                piece.range,
            );
        }
    }
    // A contact edge splits where it crosses the target face's boundary,
    // and that crossing is a pave on *both* edges, so the owner's own face
    // splits its boundary consistently and the pieces sew back shared.
    let mut contact_along: Vec<Vec<(f64, f64)>> = vec![Vec::new(); contacts.len()];
    for (ci, contact) in contacts.iter().enumerate() {
        let target = if contact.target_from_a {
            &ga.faces[contact.target_face]
        } else {
            &gb.faces[contact.target_face]
        };
        // A fitted contact sits off exact geometry by its own tolerance, and
        // the crossings it genuinely makes gape by the same: the filter and
        // the finder both widen, or a fitted rail never registers against
        // the exact seam it crosses.
        let reach = tol.confusion().max(contact.tolerance * 2.0);
        let cc = CurveCurveOptions {
            gap: reach.max(CurveCurveOptions::default().gap),
            ..CurveCurveOptions::default()
        };
        for e in &target.edges {
            let found = intersect_curves(&contact.curve, &e.curve, cc, tol)?;
            for crossing in &found.crossings {
                if crossing.gap > reach {
                    continue;
                }
                if crossing.on_a < contact.crange.0 + tol.parametric()
                    || crossing.on_a > contact.crange.1 - tol.parametric()
                {
                    continue;
                }
                // Touching is not crossing, at curve level as at surface
                // level: where the two curves run nearly parallel (a blend
                // arc meeting the edge it is tangent to), the "crossing" is
                // numerical noise smeared along the contact, and paving it
                // would plant a vertex a hair off both curves.
                if tangential(&contact.curve, crossing.on_a, &e.curve, crossing.on_b, tol)? {
                    continue;
                }
                let honesty = honest(contact.tolerance, tol).max(crossing.reach);
                paves.entry(contact.node).or_default().push(Pave {
                    t: crossing.on_a,
                    honesty,
                });
                let on_b = onto_range(crossing.on_b, &e.curve, e.crange, tol);
                if on_b > e.crange.0 + tol.parametric() && on_b < e.crange.1 - tol.parametric() {
                    paves
                        .entry(e.node)
                        .or_default()
                        .push(Pave { t: on_b, honesty });
                }
            }
            // A span of the contact running along a target boundary edge
            // splits nothing: it is already boundary on both sides
            // (identically stacked boxes are all such spans), and duplicating
            // it as a strand would cancel the boundary it copies.
            // The intersector answers coincidence in closed form for the
            // analytic pairs only; a fitted curve lying on an exact one
            // (a marched sphere's rim on a wedge cap's arc) comes back as
            // crossings or nothing. That span is measured instead: the
            // stretch of the contact within the pair's honesty of the edge.
            // An analytic overlap survives only where it lands inside the
            // contact's own window; a trimmed circle written a turn up
            // clips to nothing, and is then measured like a fitted one.
            let survives = |overlap: &ogeom_intersect::Overlap| -> bool {
                let (lo, hi) = if overlap.on_a.0 <= overlap.on_a.1 {
                    overlap.on_a
                } else {
                    (overlap.on_a.1, overlap.on_a.0)
                };
                let carried = if contact.curve.is_periodic() {
                    let (dlo, dhi) = contact.curve.domain();
                    let period = dhi - dlo;
                    let turn = ((contact.crange.0 - lo) / period).floor();
                    [turn, turn + 1.0, turn - 1.0]
                        .into_iter()
                        .map(|k| (k.mul_add(period, lo), k.mul_add(period, hi)))
                        .map(|(a, b)| b.min(contact.crange.1) - a.max(contact.crange.0))
                        .fold(f64::NEG_INFINITY, f64::max)
                } else {
                    hi.min(contact.crange.1) - lo.max(contact.crange.0)
                };
                carried > tol.parametric()
            };
            // The closed-form overlap is between the two *curves*; the
            // stretch that is boundary is what the *edge* covers of it. A
            // box cut from an L-bracket flush with the bracket's wall puts
            // the box's wall-side edge on the line of the end face's own
            // edge along the wall, a length below it: read as along that
            // edge over the whole curve, the strip's side was never paved
            // and the end face kept the strip. Clipped through the
            // overlap's own correspondence so the carry below stays affine.
            let clipped: Vec<ogeom_intersect::Overlap> = found
                .overlaps
                .iter()
                .filter_map(|overlap| {
                    let (lo, hi) = overlap_within(overlap, e.crange, &contact.curve, tol)?;
                    let span = overlap.on_a.1 - overlap.on_a.0;
                    if span.abs() <= f64::MIN_POSITIVE {
                        return None;
                    }
                    let to_b = |t: f64| {
                        overlap.on_b.0
                            + (overlap.on_b.1 - overlap.on_b.0) * (t - overlap.on_a.0) / span
                    };
                    Some(ogeom_intersect::Overlap {
                        on_a: (lo, hi),
                        on_b: (to_b(lo), to_b(hi)),
                    })
                })
                .collect();
            // A line and a curved conic share no stretch, however loosely
            // either is held: a straight edge tangent to an arc stays within
            // the measuring width of it for a short run either side of the
            // touch, and read as a span there it splits the arc at a vertex
            // the face across the arc never has.
            let conic = |c: &Curve| {
                matches!(
                    c,
                    Curve::Circle(_) | Curve::Ellipse(_) | Curve::Hyperbola(_) | Curve::Parabola(_)
                )
            };
            let straight_on_conic = (matches!(contact.curve, Curve::Line(_)) && conic(&e.curve))
                || (conic(&contact.curve) && matches!(e.curve, Curve::Line(_)));
            let measured: Vec<ogeom_intersect::Overlap> =
                if clipped.iter().any(survives) || straight_on_conic {
                    Vec::new()
                } else {
                    measured_overlaps(&contact.curve, contact.crange, contact.tolerance, e, tol)?
                };
            for overlap in clipped.iter().chain(measured.iter()) {
                let ordered = |r: (f64, f64)| if r.0 <= r.1 { r } else { (r.1, r.0) };
                let (lo, hi) = ordered(overlap.on_a);
                // The overlap is between the two *curves*; what interferes is
                // the stretch both *edges* actually cover. A hole's arc and
                // the disc that fills it lie on one circle, so the curves
                // overlap over the whole turn while the arc covers three
                // quarters of it), and paving at the turn's ends says nothing,
                // where paving at the arc's ends is exactly the split the
                // other side needs to sew against.
                // A periodic curve's overlap is answered on its base turn;
                // an edge written a turn up (a wedge cap's arc at 2π..2.5π)
                // covers it only once the turn is carried across.
                let (lo, hi) = if contact.curve.is_periodic() {
                    let (dlo, dhi) = contact.curve.domain();
                    let period = dhi - dlo;
                    let turn = ((contact.crange.0 - lo) / period).floor();
                    let mut best = (lo, hi);
                    let mut best_span = f64::NEG_INFINITY;
                    for k in [turn, turn + 1.0, turn - 1.0] {
                        let cand = (k.mul_add(period, lo), k.mul_add(period, hi));
                        let span = cand.1.min(contact.crange.1) - cand.0.max(contact.crange.0);
                        if span > best_span {
                            best_span = span;
                            best = cand;
                        }
                    }
                    best
                } else {
                    (lo, hi)
                };
                let (lo, hi) = (lo.max(contact.crange.0), hi.min(contact.crange.1));
                if hi - lo <= tol.parametric() {
                    continue;
                }
                for t in [lo, hi] {
                    paves.entry(contact.node).or_default().push(Pave {
                        t,
                        honesty: contact.tolerance,
                    });
                }
                if *DEBUG_WIRE {
                    eprintln!(
                        "CONTACT c{ci} along edge {} over ({lo:.6}, {hi:.6}) of {:?}: {:?} .. {:?}",
                        e.node.index(),
                        contact.crange,
                        contact.curve.point_at(lo, tol).ok(),
                        contact.curve.point_at(hi, tol).ok()
                    );
                }
                contact_along[ci].push((lo, hi));
                // The *target* edge splits where the shared stretch ends,
                // exactly as the contact does. Without this, the face across
                // the overlap keeps one long boundary edge where its new
                // neighbours carry two short ones, and sew (which matches
                // edges whole) can pair it with neither. The clamped ends are
                // carried across by the correspondence the overlap itself
                // states, which is affine over the shared stretch.
                let span = overlap.on_a.1 - overlap.on_a.0;
                let carry = |t: f64| -> f64 {
                    if span.abs() <= f64::MIN_POSITIVE {
                        overlap.on_b.0
                    } else {
                        overlap.on_b.0
                            + (overlap.on_b.1 - overlap.on_b.0) * (t - overlap.on_a.0) / span
                    }
                };
                let (mut tlo, mut thi) = (carry(lo), carry(hi));
                if tlo > thi {
                    core::mem::swap(&mut tlo, &mut thi);
                }
                let target_domain = e.curve.domain();
                let periodic = e.curve.is_periodic();
                for t in [tlo, thi] {
                    // A correspondence across a full turn can run the
                    // parameter past the domain; the pave belongs where the
                    // edge actually is.
                    let t = if periodic { fold(t, target_domain) } else { t };
                    if t > e.crange.0 + tol.parametric() && t < e.crange.1 - tol.parametric() {
                        paves.entry(e.node).or_default().push(Pave {
                            t,
                            honesty: contact.tolerance,
                        });
                    }
                }
            }
        }
    }

    let mut same_a: Vec<Vec<usize>> = vec![Vec::new(); ga.faces.len()];
    let mut same_b: Vec<Vec<usize>> = vec![Vec::new(); gb.faces.len()];
    for (ia, ib) in same_pairs {
        same_a[ia].push(ib);
        same_b[ib].push(ia);
    }
    Ok((
        sections,
        pieces,
        contacts,
        tangents,
        contact_along,
        paves,
        same_a,
        same_b,
        hug_junctions,
    ))
}

/// An exact section cut at the chart degeneracies it runs through, each
/// piece carrying exact pcurves on both faces.
///
/// The degeneracies are not guessed from the surfaces: they are the faces'
/// own *pole edges*, the degenerate edges the topology already carries, so a
/// face whose chart has a pole says so and one that has none costs nothing.
///
/// `None` means the split is no use here (no degeneracy on the curve, or a
/// piece whose projection still has no closed form), and the caller falls
/// back to marching, which is the honest answer rather than a fitted pcurve
/// pretending to be exact.
fn split_at_degeneracies(
    curve: &Curve,
    fa: &GFace,
    fb: &GFace,
    tol: Tolerances,
) -> OgeomResult<Option<Vec<(Curve, PlanarCurve, PlanarCurve)>>> {
    let mut stops: Vec<f64> = Vec::new();
    let domain = curve.domain();
    // The split points are the *surfaces'* chart degeneracies, not merely the
    // faces' pole edges: a meridian section runs through both of a sphere's
    // poles, and a face that owns only the north one still cannot chart an
    // arc that wraps through the south. The surface knows where its chart
    // collapses whether or not the face's trim reaches there.
    let mut candidates: Vec<Point> = Vec::new();
    for face in [fa, fb] {
        for pole in &face.poles {
            candidates.push(pole.point);
        }
        match &face.surface {
            SurfaceGeometry::Sphere(s) => {
                let sphere = s.sphere();
                let axis = sphere.frame().z().vector();
                candidates.push(sphere.centre() + axis * sphere.radius());
                candidates.push(sphere.centre() - axis * sphere.radius());
            }
            SurfaceGeometry::Cone(c) => {
                candidates.push(c.cone().apex());
            }
            _ => {}
        }
    }
    for point in candidates {
        let found = ogeom_algo::project_on_curve(curve, point, 256, tol)?;
        if found.distance <= tol.confusion() {
            stops.push(found.parameter);
        }
    }
    if stops.is_empty() {
        return Ok(None);
    }
    stops.sort_by(|a, b| a.partial_cmp(b).unwrap_or(core::cmp::Ordering::Equal));
    stops.dedup_by(|a, b| (*a - *b).abs() <= tol.parametric());

    // A closed curve is cut into the arcs between successive stops, the last
    // wrapping past the domain end; an open one keeps its own ends as stops.
    let closed = curve.is_closed(tol) || curve.is_periodic();
    let mut arcs: Vec<(f64, f64)> = Vec::new();
    if closed {
        if stops.len() < 2 {
            // One stop on a closed curve leaves one arc, from the stop right
            // round to itself, which still has the pole at both ends.
            let period = domain.1 - domain.0;
            arcs.push((stops[0], stops[0] + period));
        } else {
            let period = domain.1 - domain.0;
            for i in 0..stops.len() {
                let lo = stops[i];
                let hi = if i + 1 < stops.len() {
                    stops[i + 1]
                } else {
                    stops[0] + period
                };
                arcs.push((lo, hi));
            }
        }
    } else {
        let mut cuts = vec![domain.0];
        cuts.extend(
            stops
                .iter()
                .copied()
                .filter(|t| *t > domain.0 + tol.parametric() && *t < domain.1 - tol.parametric()),
        );
        cuts.push(domain.1);
        for pair in cuts.windows(2) {
            arcs.push((pair[0], pair[1]));
        }
    }

    let mut out = Vec::with_capacity(arcs.len());
    for (lo, hi) in arcs {
        if hi - lo <= tol.parametric() {
            continue;
        }
        let (Some(pa), Some(pb)) = (
            ogeom_intersect::exact_pcurve_over(curve, (lo, hi), &fa.surface, tol),
            ogeom_intersect::exact_pcurve_over(curve, (lo, hi), &fb.surface, tol),
        ) else {
            return Ok(None);
        };
        let Ok(piece) = ogeom_geom::TrimmedCurve::new(curve.clone(), lo, hi, tol) else {
            return Ok(None);
        };
        out.push((piece.into(), pa, pb));
    }
    if out.is_empty() {
        return Ok(None);
    }
    Ok(Some(out))
}

/// A surface narrowed to the reach of a bound, for the marcher's benefit.
///
/// Seeding samples a surface's *parameter box*, so a plane stored over
/// ±10^9 (which is how an unbounded carrier reaches this code) is sampled
/// at a spacing a hundred million times anything it could be meeting, and
/// no seed ever lands near the curve. Narrowing to where the two faces
/// actually are is what lets the seeding see it. The chart is untouched:
/// only the window the marcher walks changes, so pcurves fitted in it mean
/// the same thing on the surface as stored.
fn windowed_to(surface: &SurfaceGeometry, bound: &ogeom_math::Aabb) -> SurfaceGeometry {
    let corners = bound.corners();
    if corners.is_empty() {
        return surface.clone();
    }
    match surface {
        SurfaceGeometry::Plane(p) => {
            let frame = p.plane().frame();
            let (mut u0, mut u1, mut v0, mut v1) = (
                f64::INFINITY,
                f64::NEG_INFINITY,
                f64::INFINITY,
                f64::NEG_INFINITY,
            );
            for c in &corners {
                let local = frame.to_local(*c);
                u0 = u0.min(local.x);
                u1 = u1.max(local.x);
                v0 = v0.min(local.y);
                v1 = v1.max(local.y);
            }
            let margin = ((u1 - u0) + (v1 - v0)).mul_add(0.25, 1.0);
            let (want_u, want_v) = ((u0 - margin, u1 + margin), (v0 - margin, v1 + margin));
            let (have_u, have_v) = ogeom_geom::Surface::domain(p);
            if have_u.1 - have_u.0 <= want_u.1 - want_u.0
                && have_v.1 - have_v.0 <= want_v.1 - want_v.0
            {
                return surface.clone();
            }
            ogeom_geom::PlaneSurface::over(p.plane(), want_u, want_v)
                .map_or_else(|_| surface.clone(), Into::into)
        }
        SurfaceGeometry::Cylinder(c) => {
            let frame = c.cylinder().frame();
            let (mut h0, mut h1) = (f64::INFINITY, f64::NEG_INFINITY);
            for corner in &corners {
                let h = (*corner - frame.origin()).dot(frame.z().vector());
                h0 = h0.min(h);
                h1 = h1.max(h);
            }
            let margin = (h1 - h0).mul_add(0.25, 1.0);
            let want = (h0 - margin, h1 + margin);
            let have = ogeom_geom::Surface::domain(c).1;
            if have.1 - have.0 <= want.1 - want.0 {
                return surface.clone();
            }
            ogeom_geom::CylinderSurface::new(c.cylinder(), want)
                .map_or_else(|_| surface.clone(), Into::into)
        }
        other => other.clone(),
    }
}

/// March a pair whose exact section has no closed-form pcurve.
fn march_pair(
    a: &SurfaceGeometry,
    b: &SurfaceGeometry,
    options: &ogeom_intersect::IntersectOptions,
    tol: Tolerances,
) -> OgeomResult<Vec<ogeom_intersect::IntersectionCurve>> {
    use ogeom_intersect::{approximate_branch, branches};
    let traced = branches(a, b, options.marching, tol)?;
    let mut out = Vec::new();
    for branch in &traced {
        out.push(approximate_branch(a, b, branch, options.tolerance, tol)?);
    }
    if out.is_empty() {
        // No branch found is two different stories. Surfaces that measurably
        // stand apart over their stated extents (a blend's cylinder and a
        // far-off bevel plane whose ellipse of intersection lies beyond
        // both) simply do not interact, and an empty section is the true
        // answer. Only a pair that comes close and still resolves nothing
        // is beyond the intersector.
        if surfaces_stand_apart(a, b, tol) {
            return Ok(out);
        }
        ogeom_bail!(
            NotDone,
            "an exact section has no closed-form pcurve and marching resolved \
             no branch; the configuration is beyond the intersector's current \
             reach"
        );
    }
    Ok(out)
}

/// Whether two surfaces measurably keep their distance over their stated
/// extents.
///
/// A conservative grid measurement: sample the smaller-extent surface and
/// project each sample onto the other; apart means every sample clears a
/// margin scaled to the extents. Surfaces with unbounded or enormous stated
/// domains (an imported plane's billion units) are never called apart this
/// way, because a grid over them samples nothing.
fn surfaces_stand_apart(a: &SurfaceGeometry, b: &SurfaceGeometry, tol: Tolerances) -> bool {
    use ogeom_geom::Surface as _;
    let extent_of = |s: &SurfaceGeometry| -> f64 {
        let ((ua, ub), (va, vb)) = s.domain();
        (ub - ua).abs().max((vb - va).abs())
    };
    let (sample, against) = if extent_of(a) <= extent_of(b) {
        (a, b)
    } else {
        (b, a)
    };
    let span = extent_of(sample);
    if !span.is_finite() || span > 1e4 {
        return false;
    }
    let ((ua, ub), (va, vb)) = sample.domain();
    const GRID: usize = 9;
    // One seeding grid over the far surface, asked a hundred times: the
    // same seeds and the same Newton the per-call projection would use, so
    // the verdict is bit-identical, at hundreds of evaluations instead of
    // tens of thousands for every genuine near-miss.
    let Ok(seeds) = ogeom_algo::SurfaceSeeds::over(against, 16, tol) else {
        return false;
    };
    let mut clearance = f64::INFINITY;
    for i in 0..=GRID {
        for j in 0..=GRID {
            #[allow(clippy::cast_precision_loss)]
            let u = ua + (ub - ua) * i as f64 / GRID as f64;
            #[allow(clippy::cast_precision_loss)]
            let v = va + (vb - va) * j as f64 / GRID as f64;
            let Ok(p) = sample.point_at(u, v, tol) else {
                return false;
            };
            let Ok(projection) = seeds.project(against, p, tol) else {
                return false;
            };
            clearance = clearance.min(projection.distance);
            if clearance <= tol.confusion() * 1e3 {
                return false;
            }
        }
    }
    clearance > tol.confusion() * 1e3
}

/// Whether two surfaces are one surface wherever they overlap, measured.
///
/// The closed forms answer this for the pairs they know. Where there is no
/// closed form it does not stop being a fair question (two patches restated
/// from one plane are the same surface, and nothing in their control points
/// says so), but it stops being answerable exactly, so it is measured here
/// and only for the pairs the analytic layer has already declined.
///
/// Sampled on the *smaller* window, because the answer is about the region
/// the two share and a stated window is not that region: a plane's own
/// extends for a billion units either way, and a grid over it samples
/// nothing. A sample whose foot lands on the rim of the other is skipped
/// rather than counted against: the other patch simply does not reach that
/// far, and a distance measured to its rim is about the window, not the
/// surface.
///
/// One-sided by construction: a pair that crosses puts interior samples well
/// off the other, so it cannot pass, and a pair this cannot resolve marches
/// exactly as it did before.
fn surfaces_coincide(
    a: &SurfaceGeometry,
    b: &SurfaceGeometry,
    reach: f64,
    tol: Tolerances,
) -> bool {
    use ogeom_geom::Surface as _;
    /// Samples per direction over the window, and how many must land inside
    /// the other before agreement means anything.
    const GRID: usize = 6;
    const EVIDENCE: usize = 4;

    let span = |s: &SurfaceGeometry| -> f64 {
        let ((ua, ub), (va, vb)) = s.domain();
        (ub - ua).abs().max((vb - va).abs())
    };
    let (sampled, against) = if span(a) <= span(b) { (a, b) } else { (b, a) };
    let ((ua, ub), (va, vb)) = sampled.domain();
    if !(ua.is_finite() && ub.is_finite() && va.is_finite() && vb.is_finite()) {
        return false;
    }
    let ((wu0, wu1), (wv0, wv1)) = against.domain();
    // A twentieth of the window in from each rim: enough that a foot the
    // search pinned to the rim is not read as one the surface truly reaches.
    let (mu, mv) = ((wu1 - wu0) * 0.05, (wv1 - wv0) * 0.05);

    let mut evidence = 0_usize;
    for i in 0..=GRID {
        for j in 0..=GRID {
            #[allow(clippy::cast_precision_loss)]
            let u = ua + (ub - ua) * (i as f64 / GRID as f64);
            #[allow(clippy::cast_precision_loss)]
            let v = va + (vb - va) * (j as f64 / GRID as f64);
            let Ok(p) = sampled.point_at(u, v, tol) else {
                return false;
            };
            let Ok(foot) = ogeom_algo::project_on_surface(against, p, 16, tol) else {
                return false;
            };
            let (fu, fv) = foot.parameters;
            if fu <= wu0 + mu || fu >= wu1 - mu || fv <= wv0 + mv || fv >= wv1 - mv {
                continue;
            }
            if foot.distance > reach {
                return false;
            }
            evidence += 1;
        }
    }
    evidence >= EVIDENCE
}

// --- the general fuse --------------------------------------------------------

/// One piece of one argument face, classified against the other argument.
struct FacePiece {
    /// Which argument's face list, and which face.
    from_a: bool,
    face: usize,
    rings: Vec<Vec<Traversal<Tag>>>,
    /// The rings as chart polylines, for the coincidence pairing below.
    outlines: Vec<Vec<Point2>>,
    /// The interior point the state was decided at, in space.
    probe: Point,
    state: PieceState,
    /// Whether the *other* argument already contributes this same patch of
    /// this same surface.
    ///
    /// Two faces lying on one surface with their material on the same side
    /// bound the union and the intersection once between them, not twice,
    /// but "once" is a statement about the *patch*, not about which argument
    /// it came from. Where a tool's cap fills a hole its own bore left, the
    /// part has no piece there at all, and dropping the tool's by argument
    /// identity leaves the result with a hole in it. So the duplicate is
    /// identified by containment, and only a piece genuinely stood in for is
    /// dropped.
    covered: bool,
}

struct GeneralFused {
    a: GSolid,
    b: GSolid,
    sections: Vec<SectionRec>,
    contacts: Vec<ContactRec>,
    /// Curves the two boundaries touch along without crossing. They take no
    /// part in the classification below (that is what tangency means) and
    /// are carried for the consumers that want the contact itself.
    tangents: Vec<TangentRec>,
    pieces: Vec<FacePiece>,
    /// Junctions several paves describe, each resolved once for every
    /// strand that ends in it.
    junctions: Vec<Junction>,
}

/// One junction the paving found several times over.
///
/// A tolerant rail meeting a wedge's faces near one corner collects a
/// cluster of crossings inside its own stated radius: each face's section
/// stops at its own crossing with the rail, and the crossings sit a few
/// tenths of a micron to a few hundred apart along it. The rail's strands
/// split once per cluster, at its first pave; the sections' ends still name
/// their own crossings; and the two descriptions of the junction can sit
/// apart by the cluster's whole span, which no single strand's honesty
/// covers. So the junction is one vertex standing at the first pave, owning
/// the span the paves disagree by, and every strand end inside that span
/// names it, whichever pave its own trim stopped at.
#[derive(Debug, Clone, Copy)]
struct Junction {
    /// Where the cluster's first pave sits.
    at: Point,
    /// How far a strand end may sit from `at` and still be this junction:
    /// the cluster's span plus the rail's own reach.
    reach: f64,
    /// Whether `at` is an edge's own end vertex that a section tangent to
    /// the edge was welded onto: the section's chart image is bent onto it
    /// rather than left where the section stopped.
    onto_vertex: bool,
}

/// Junctions whose balls overlap, merged transitively into one junction
/// each: at the members' centroid, reaching as far as the farthest member's
/// own reach extends from it.
fn merge_junctions(junctions: Vec<Junction>) -> Vec<Junction> {
    let n = junctions.len();
    let mut parent: Vec<usize> = (0..n).collect();
    fn root(parent: &mut [usize], mut i: usize) -> usize {
        while parent[i] != i {
            parent[i] = parent[parent[i]];
            i = parent[i];
        }
        i
    }
    // Swept along x: two balls overlap only if their centres are within
    // both reaches along every axis, so each junction is asked only of those
    // ahead of it by no more than its reach and the widest. The groups are
    // the same as asking every pair; a converted part carries tens of
    // thousands of junctions, and every pair was seconds.
    let widest = junctions.iter().map(|j| j.reach).fold(0.0_f64, f64::max);
    let mut order: Vec<usize> = (0..n).collect();
    order.sort_by(|&x, &y| junctions[x].at.x.total_cmp(&junctions[y].at.x));
    for (k, &i) in order.iter().enumerate() {
        let limit = junctions[i].at.x + junctions[i].reach + widest;
        for &j in &order[k + 1..] {
            if junctions[j].at.x > limit {
                break;
            }
            if junctions[i].at.distance(junctions[j].at) <= junctions[i].reach + junctions[j].reach
            {
                let (ri, rj) = (root(&mut parent, i), root(&mut parent, j));
                if ri != rj {
                    parent[rj] = ri;
                }
            }
        }
    }
    let mut groups: Vec<Vec<usize>> = Vec::new();
    let mut of_root: Vec<Option<usize>> = vec![None; n];
    for i in 0..n {
        let r = root(&mut parent, i);
        match of_root[r] {
            Some(g) => groups[g].push(i),
            None => {
                of_root[r] = Some(groups.len());
                groups.push(vec![i]);
            }
        }
    }
    // A group is one junction only while it stays a junction's size: a
    // chain of overlapping balls along an edge would otherwise merge into
    // one reaching the chain's length, and the rebuild would weld the edge
    // into a point. A group reaching beyond four of its widest member stays
    // as it was.
    let mut out = Vec::with_capacity(groups.len());
    for members in groups {
        if members.len() == 1 {
            out.push(junctions[members[0]]);
            continue;
        }
        let mut sum = ogeom_math::Vector::ZERO;
        for &m in &members {
            sum += junctions[m].at - Point::ORIGIN;
        }
        #[allow(clippy::cast_precision_loss)]
        let at = Point::ORIGIN + sum / (members.len() as f64);
        let reach = members
            .iter()
            .map(|&m| at.distance(junctions[m].at) + junctions[m].reach)
            .fold(0.0_f64, f64::max);
        let widest = members
            .iter()
            .map(|&m| junctions[m].reach)
            .fold(0.0_f64, f64::max);
        if reach > widest * 4.0 {
            out.extend(members.iter().map(|&m| junctions[m]));
        } else {
            out.push(Junction {
                at,
                reach,
                onto_vertex: members.iter().any(|&m| junctions[m].onto_vertex),
            });
        }
    }
    out
}

/// The junctions the paves describe more than once, or more loosely than
/// the edge itself, per edge in face order.
fn pave_junctions(
    ga: &GSolid,
    gb: &GSolid,
    paves: &std::collections::HashMap<EdgeKey, Vec<Pave>>,
    tol: Tolerances,
) -> OgeomResult<Vec<Junction>> {
    let mut seen: std::collections::HashSet<EdgeKey> = std::collections::HashSet::new();
    let mut junctions = Vec::new();
    for e in ga
        .faces
        .iter()
        .chain(gb.faces.iter())
        .flat_map(|f| f.edges.iter())
    {
        if !seen.insert(e.node) {
            continue;
        }
        let Some(ts) = paves.get(&e.node) else {
            continue;
        };
        // A tolerant edge (a merged facet group's boundary, a fitted rail)
        // is met by a section up to its own radius off it, so the section's
        // end and the edge's split point are one junction that far apart,
        // on either side of the edge.
        let tolerant = e.tolerance > tol.confusion() * 1e2;
        let floor = if tolerant {
            e.tolerance * 2.0
        } else {
            e.tolerance
        }
        .max(tol.confusion() * 10.0);
        for cluster in cluster_paves(&e.curve, e.crange, e.tolerance, ts, tol)? {
            if cluster.members > 1 || cluster.honesty > tol.confusion() * 1e2 || tolerant {
                junctions.push(Junction {
                    at: cluster.at,
                    reach: cluster.span + cluster.honesty.max(floor),
                    onto_vertex: false,
                });
            }
        }
    }
    Ok(junctions)
}

/// The face's outward normal at a chart point: the surface's, flipped when
/// the face presents its other side.
fn outward_normal(face: &GFace, at: Point2, tol: Tolerances) -> OgeomResult<ogeom_math::Vector> {
    let n = face.surface.normal_at(at.x, at.y, tol)?.vector();
    Ok(
        if face.face.orientation() == ogeom_topo::Orientation::Reversed {
            -n
        } else {
            n
        },
    )
}

/// The chart point of a world point lying on a planar face, if it lands
/// inside the face's trim.
/// How far apart two of a face's edges may honestly end in its chart: the
/// loosest edge's tolerance, or the hug width where an edge is a fitted
/// section that once ran along another before parting from it.
fn outline_snap(face: &GFace, tol: Tolerances) -> f64 {
    face.edges
        .iter()
        .fold(tol.confusion() * 1e2, |acc, e| acc.max(e.tolerance * 2.0))
        .max(
            if face.edges.iter().any(|e| e.tolerance > tol.confusion()) {
                tol.confusion() * 1e3
            } else {
                0.0
            },
        )
}

/// Close the gaps between a face's outline polylines.
///
/// The edges are polylined one by one, and where two of them meet at a
/// vertex that owns some tolerance (a fitted section's end welded to the
/// edge it hugged, a few dozen microns off), their polylines stop that far
/// apart. A ray cast for containment slips through such a gap, and a probe
/// standing inside the face within a hug's width of the seam reads as
/// outside. Ends within `snap` of another polyline's end are made one
/// point, so the outline is closed exactly as the topology says it is.
fn weld_outline_ends(lines: &mut [Vec<Point2>], snap: f64) {
    if snap <= 0.0 {
        return;
    }
    // Every polyline end, as (line, is its last point, where).
    let ends: Vec<(usize, bool, Point2)> = lines
        .iter()
        .enumerate()
        .flat_map(|(i, line)| {
            let first = line.first().map(|p| (i, false, *p));
            let last = line.last().map(|p| (i, true, *p));
            first.into_iter().chain(last)
        })
        .collect();
    // Each end moves onto the nearest other end within reach that sorts
    // before it, so a matched pair lands on one point rather than trading
    // places.
    let mut moves: Vec<(usize, bool, Point2)> = Vec::new();
    for &(i, end_i, p) in &ends {
        let nearest = ends
            .iter()
            // Ends closer than the parametric snap already meet for every
            // purpose here, and moving one by rounding noise takes a probe
            // that stands exactly on a chart's seam off it.
            .filter(|&&(j, end_j, q)| {
                let d = p.distance(q);
                (j, end_j) < (i, end_i) && j != i && d > PARAM_SNAP && d <= snap
            })
            .min_by(|a, b| {
                p.distance(a.2)
                    .partial_cmp(&p.distance(b.2))
                    .unwrap_or(core::cmp::Ordering::Equal)
            });
        if let Some(&(_, _, q)) = nearest {
            moves.push((i, end_i, q));
        }
    }
    for (i, end_i, q) in moves {
        let slot = if end_i {
            lines[i].last_mut()
        } else {
            lines[i].first_mut()
        };
        if let Some(p) = slot {
            *p = q;
        }
    }
}

fn chart_point_of(face: &GFace, p: Point, tol: Tolerances) -> Option<Point2> {
    // Closed-form inversion for the analytic surfaces: the same-domain
    // resolution asks "where does this probe sit in the partner's chart", and
    // the partner may be any surface a face melts along: a plane against a
    // plane, but equally a wall band against the cylinder it copies. The
    // reach check keeps the answer honest: a point off the surface has no
    // chart position, whatever the inversion returns.
    use ogeom_math::elementary;
    let reach = tol.confusion() * 10.0;
    let raw = match &face.surface {
        SurfaceGeometry::Plane(x) => {
            let local = x.plane().frame().to_local(p);
            if local.z.abs() > reach {
                return None;
            }
            Point2::new(local.x, local.y)
        }
        SurfaceGeometry::Cylinder(x) => {
            let cylinder = x.cylinder();
            if cylinder.distance_to(p) > reach {
                return None;
            }
            let (u, v) = elementary::cylinder_parameters(&cylinder, p, tol).ok()?;
            Point2::new(u, v)
        }
        SurfaceGeometry::Cone(x) => {
            let cone = x.cone();
            if cone.distance_to(p) > reach {
                return None;
            }
            let (u, v) = elementary::cone_parameters(&cone, p, tol).ok()?;
            Point2::new(u, v)
        }
        SurfaceGeometry::Sphere(x) => {
            let sphere = x.sphere();
            if sphere.distance_to(p) > reach {
                return None;
            }
            let (u, v) = elementary::sphere_parameters(&sphere, p, tol).ok()?;
            Point2::new(u, v)
        }
        SurfaceGeometry::Torus(x) => {
            let torus = x.torus();
            if torus.distance_to(p) > reach {
                return None;
            }
            let (u, v) = elementary::torus_parameters(&torus, p, tol).ok()?;
            Point2::new(u, v)
        }
        // No closed form (a fitted patch, a swept or revolved surface)
        // inverts by projection, held to the same reach: a blend's leg on
        // a spline host is a partner like any other.
        other => {
            let foot = ogeom_algo::project_on_surface(other, p, 24, tol).ok()?;
            if foot.distance > reach {
                return None;
            }
            Point2::new(foot.parameters.0, foot.parameters.1)
        }
    };
    let at = fold_point_into_chart(raw, &face.surface);
    let mut lines: Vec<Vec<Point2>> = Vec::new();
    for e in &face.edges {
        lines.push(
            pcurve_polyline(&e.pcurve, e.prange, e.crange, e.crange, &face.surface, tol).ok()?,
        );
        // A seam bounds the chart twice (once per column), and a trim test
        // that sees only one side reads half the band as outside.
        if let Some((other, orange)) = &e.other_side {
            lines.push(
                pcurve_polyline(other, *orange, e.crange, e.crange, &face.surface, tol).ok()?,
            );
        }
    }
    weld_outline_ends(&mut lines, outline_snap(face, tol));
    let borrowed: Vec<&[Point2]> = lines.iter().map(Vec::as_slice).collect();
    // The face's boundary polylines are unwrapped (a winding ring may span
    // any one period's window, not necessarily the chart's canonical one,
    // and a wire chained onto one branch of the chart may sit whole periods
    // away from it), so the probe is carried to the outline's own branch
    // first, and asked at the neighbouring images as well.
    let mut shifts = vec![0.0];
    if face.surface.is_periodic_u() {
        let ((ua, ub), _) = face.surface.domain();
        if ub > ua {
            let period = ub - ua;
            let (sum, count) = lines
                .iter()
                .flatten()
                .fold((0.0_f64, 0_usize), |(s, n), q| (s + q.x, n + 1));
            #[allow(clippy::cast_precision_loss)]
            let centre = if count > 0 { sum / count as f64 } else { at.x };
            let home = ((centre - at.x) / period).round() * period;
            shifts = vec![home, home + period, home - period];
        }
    }
    for shift in shifts {
        let shifted = Point2::new(at.x + shift, at.y);
        if inside_many(&borrowed, shifted) {
            return Some(shifted);
        }
    }
    None
}

/// Whether two curves meet tangentially at a crossing: the parallel-noise
/// gate for pave placement. One degree is far below any deliberate crossing
/// and far above the smear a tangency leaves in the general intersector.
fn tangential(a: &Curve, ta: f64, b: &Curve, tb: f64, tol: Tolerances) -> OgeomResult<bool> {
    let da = a.d1_at(ta, tol)?;
    let db = b.d1_at(tb, tol)?;
    let (ma, mb) = (da.magnitude(), db.magnitude());
    if ma <= tol.confusion() || mb <= tol.confusion() {
        return Ok(true);
    }
    Ok(da.cross(db).magnitude() / (ma * mb) < 2e-2)
}

/// The chart image of an exact curve on a surface, fitted at the curve's
/// own parameters from its points inverted on the surface, with the
/// distance the image honestly sits off the curve, or nothing, when the
/// fit misses its budget by more than a chart's swing beside a pole allows.
fn fitted_image(
    curve: &Curve,
    surface: &SurfaceGeometry,
    budget: f64,
    tol: Tolerances,
) -> OgeomResult<Option<(PlanarCurve, f64)>> {
    use ogeom_geom::Surface as _;
    const SAMPLES: usize = 96;
    let (lo, hi) = curve.domain();
    if !hi.is_finite() || !lo.is_finite() || hi <= lo {
        return Ok(None);
    }
    let ((ua, ub), (va, vb)) = surface.domain();
    let periods = (
        if surface.is_periodic_u() {
            ub - ua
        } else {
            0.0
        },
        if surface.is_periodic_v() {
            vb - va
        } else {
            0.0
        },
    );
    // The image sampled where it moves: beside a pole a circle's image
    // swings half a turn within a hair of arc, and a uniform sampling puts
    // one point in the swing. Between any two samples whose images sit
    // further apart than a small step of the chart a sample is added, until
    // the image is walked at that step or the budget of samples is spent.
    let invert = |t: f64, guess: Option<(f64, f64)>| -> OgeomResult<Option<((f64, f64), f64)>> {
        let p = curve.point_at(t, tol)?;
        // A plane's image is its frame's own coordinates, exact and
        // unbounded: the face's chart window is a trim, not a limit on
        // where the curve may be read, and a projection clamped to it puts
        // the circle's far side off the surface by the window's shortfall.
        if let SurfaceGeometry::Plane(plane) = surface {
            let local = plane.plane().frame().to_local(p);
            let back = plane
                .plane()
                .frame()
                .to_world(ogeom_math::Point::new(local.x, local.y, 0.0));
            return Ok(Some(((local.x, local.y), back.distance(p))));
        }
        // A sphere's image is its frame's longitude and latitude, exact
        // where the projection's Newton cannot refine beside a pole and
        // hands back a coarse scan sample instead; the reading is checked
        // round the trip and left to the projection if the chart's
        // convention is not this one.
        if let SurfaceGeometry::Sphere(ball) = surface {
            let frame = ball.sphere().frame();
            let local = frame.to_local(p);
            let flat = local.x.hypot(local.y);
            // Kept a hair inside the poles: a spline fitted through points
            // on the pole itself overshoots it between them, and the chart
            // refuses a latitude past its end.
            let half = core::f64::consts::FRAC_PI_2 - 1e-6;
            let v = local.z.atan2(flat).clamp(-half, half);
            let u = local.y.atan2(local.x).rem_euclid(core::f64::consts::TAU);
            let back = surface.point_at(u, v, tol)?;
            if back.distance(p) <= budget.max(tol.confusion() * 1e2) {
                return Ok(Some(((u, v), back.distance(p))));
            }
        }
        let mut foot = match guess {
            Some(g) => ogeom_algo::project_on_surface_from(surface, p, g, tol)?,
            None => ogeom_algo::project_on_surface(surface, p, 24, tol)?,
        };
        if foot.distance > budget.max(tol.confusion() * 1e2) && guess.is_some() {
            foot = ogeom_algo::project_on_surface(surface, p, 48, tol)?;
        }
        if foot.distance > budget.max(tol.confusion() * 1e2) {
            if *DEBUG_WIRE {
                eprintln!(
                    "IMAGE FIT: a point sits {:.2e} off the surface; no image",
                    foot.distance
                );
            }
            return Ok(None);
        }
        Ok(Some((foot.parameters, foot.distance)))
    };
    let unwrap = |uv: (f64, f64), before: Option<Point2>| -> Point2 {
        let mut uv = Point2::new(uv.0, uv.1);
        if let Some(prev) = before {
            for (coord, period, before) in [
                (&mut uv.x, periods.0, prev.x),
                (&mut uv.y, periods.1, prev.y),
            ] {
                if period > 0.0 {
                    while *coord - before > period / 2.0 {
                        *coord -= period;
                    }
                    while before - *coord > period / 2.0 {
                        *coord += period;
                    }
                }
            }
        }
        uv
    };
    let mut samples: Vec<(f64, Point2)> = Vec::with_capacity(SAMPLES + 1);
    let mut guess: Option<(f64, f64)> = None;
    for k in 0..=SAMPLES {
        #[allow(clippy::cast_precision_loss)]
        let t = lo + (hi - lo) * (k as f64) / (SAMPLES as f64);
        let Some((uv, _)) = invert(t, guess)? else {
            return Ok(None);
        };
        guess = Some(uv);
        let before = samples.last().map(|(_, p)| *p);
        samples.push((t, unwrap(uv, before)));
    }
    let step = {
        let (su, sv) = (
            if periods.0 > 0.0 { periods.0 } else { ub - ua },
            if periods.1 > 0.0 { periods.1 } else { vb - va },
        );
        su.min(sv) * 0.01
    };
    const CAP: usize = 4096;
    let mut refined = true;
    while refined && samples.len() < CAP {
        refined = false;
        let mut next: Vec<(f64, Point2)> = Vec::with_capacity(samples.len() * 2);
        for pair in samples.windows(2) {
            let ((t0, p0), (t1, p1)) = (pair[0], pair[1]);
            next.push((t0, p0));
            if p0.distance(p1) > step && t1 - t0 > (hi - lo) * 1e-7 && next.len() < CAP {
                let tm = f64::midpoint(t0, t1);
                let Some((uv, _)) = invert(tm, Some((p0.x, p0.y)))? else {
                    return Ok(None);
                };
                next.push((tm, unwrap(uv, Some(p0))));
                refined = true;
            }
        }
        next.push(samples[samples.len() - 1]);
        samples = next;
    }
    let params: Vec<f64> = samples.iter().map(|(t, _)| *t).collect();
    let image: Vec<Point2> = samples.iter().map(|(_, p)| *p).collect();
    let fitted =
        ogeom_geom::fit::fit_points_2d_at(&params, &image, 3, tol.confusion() * 10.0, tol)?;
    let planar: PlanarCurve = fitted.curve.into();
    // The honest error is in space: the surface read through the image
    // against the curve itself, between the samples as well as at them.
    // Checked at the samples and between each pair of them: where the
    // image swings the samples are dense, and a uniform check would step
    // over the swing and read the fit as honest.
    // An image that leaves the chart between its samples (a spline
    // overshooting the pole it was fitted up to) is no image either.
    let mut off = 0.0_f64;
    let mut checks: Vec<f64> = Vec::with_capacity(params.len() * 2);
    for pair in params.windows(2) {
        checks.push(pair[0]);
        checks.push(f64::midpoint(pair[0], pair[1]));
    }
    if let Some(&t) = params.last() {
        checks.push(t);
    }
    for t in checks {
        let uv = planar.point_at(t, tol)?;
        let Ok(on_surface) = surface.point_at(uv.x, uv.y, tol) else {
            if *DEBUG_WIRE {
                eprintln!("IMAGE FIT: the image leaves the chart at {uv:?}; no image");
            }
            return Ok(None);
        };
        off = off.max(on_surface.distance(curve.point_at(t, tol)?));
    }
    if *DEBUG_WIRE {
        eprintln!(
            "IMAGE FIT: image off the curve by {off:.2e} against {:.2e} (fit error {:.2e} met {})",
            budget.max(tol.confusion() * 1e3),
            fitted.error,
            fitted.met
        );
    }
    if off > budget.max(tol.confusion() * 1e3) {
        return Ok(None);
    }
    Ok(Some((planar, off)))
}

/// The distance from a point to a bounded edge curve, through a sampling
/// fine enough for the along-boundary question it answers.
fn distance_to_edge_curve(
    curve: &Curve,
    crange: (f64, f64),
    p: Point,
    tol: Tolerances,
) -> OgeomResult<f64> {
    // Coarse bracket, then two rounds of local refinement: the answer feeds
    // the hug filters, whose widths are fractions of a millimetre, and a
    // long arc's 48-segment polyline sags by more than that on its own.
    let scan = |lo: f64, hi: f64, steps: u32| -> OgeomResult<(f64, f64)> {
        let mut best = f64::INFINITY;
        let mut best_t = lo;
        let mut previous: Option<(f64, Point)> = None;
        for i in 0..=steps {
            let t = lo + (hi - lo) * f64::from(i) / f64::from(steps);
            let at = curve.point_at(t, tol)?;
            if let Some((t0, last)) = previous {
                let d = at - last;
                let len2 = d.dot(d);
                let s = if len2 > 0.0 {
                    ((p - last).dot(d) / len2).clamp(0.0, 1.0)
                } else {
                    0.0
                };
                let dist = p.distance(last + d * s);
                if dist < best {
                    best = dist;
                    best_t = (t - t0).mul_add(s, t0);
                }
            }
            previous = Some((t, at));
        }
        Ok((best, best_t))
    };
    let span = crange.1 - crange.0;
    let (_, t1) = scan(crange.0, crange.1, 48)?;
    let step = span / 48.0;
    let (_, t2) = scan((t1 - step).max(crange.0), (t1 + step).min(crange.1), 16)?;
    let fine = span / (48.0 * 8.0);
    let (best, _) = scan((t2 - fine).max(crange.0), (t2 + fine).min(crange.1), 16)?;
    Ok(best)
}

/// Pair up the pieces two coincident faces contribute for the same patch of
/// one surface, and mark the second argument's copy as stood in for.
///
/// The substitution is by *region*, not by argument: a piece of `b` is a
/// duplicate exactly where a piece of `a`, on the same side and on a surface
/// they share, already covers the point `b`'s piece was classified at. Where
/// `a` has nothing there (a bore refilled by the cylinder that cut it, whose
/// caps fill holes the part no longer has faces for), nothing stands in, and
/// `b`'s piece is the only description of that patch there is.
fn mark_covered_coincidences(ga: &GSolid, pieces: &mut [FacePiece], tol: Tolerances) {
    // Which pieces of the first argument stand on a shared surface.
    let from_a: Vec<(usize, usize)> = pieces
        .iter()
        .enumerate()
        .filter(|(_, p)| p.from_a && p.state == PieceState::OnAligned)
        .map(|(i, p)| (i, p.face))
        .collect();
    if from_a.is_empty() {
        return;
    }
    let mut covered: Vec<usize> = Vec::new();
    for (index, piece) in pieces.iter().enumerate() {
        if piece.from_a || piece.state != PieceState::OnAligned {
            continue;
        }
        for &(other, face_a) in &from_a {
            let host = &ga.faces[face_a];
            // The point `b`'s piece stands at, read in `a`'s face's chart. A
            // point that is not on that surface at all has no chart position,
            // and `chart_point_of` says so.
            let Some(at) = chart_point_of(host, piece.probe, tol) else {
                continue;
            };
            // The host piece's rings, not its boundary strands: a ring is a
            // closed loop whose last point does not repeat its first, so the
            // test has to close it. Asked as though the rings were strands
            // that jointly close, the segment from the ring's end back to its
            // start goes uncounted, and a point the ring plainly encloses
            // comes back outside whenever that missing segment would have
            // been crossed.
            if inside_rings(&pieces[other].outlines, at) {
                covered.push(index);
                break;
            }
        }
    }
    for index in covered {
        pieces[index].covered = true;
    }
}

/// The debug dumps, read once rather than once per face.
///
/// `env::var` takes a process-wide lock and allocates; the strand dump asked
/// it inside the per-face loop, where a large model asks thousands of times
/// to be told no.
/// The face-bound filter's audit: admit every pair, and name any the filter
/// would have dropped that then produces a record. A conservative filter is
/// a correctness precondition (a pair wrongly dropped is absorbed by the
/// empty-result fallback today, and would be a wrong solid if that fallback
/// ever came up empty too), and this is the check that makes the
/// precondition falsifiable. Costs one branch per pair when off.
static AUDIT_BOUNDS: std::sync::LazyLock<bool> =
    std::sync::LazyLock::new(|| std::env::var("OGEOM_BOOL_AUDIT_BOUNDS").is_ok());
static DEBUG_WIRE: std::sync::LazyLock<bool> =
    std::sync::LazyLock::new(|| std::env::var("OGEOM_DEBUG_WIRE").is_ok());
static DEBUG_STRANDS: std::sync::LazyLock<bool> =
    std::sync::LazyLock::new(|| std::env::var("OGEOM_DEBUG_STRANDS").is_ok());
static ARRANGE_DEBUG: std::sync::LazyLock<bool> =
    std::sync::LazyLock::new(|| std::env::var("OGEOM_ARRANGE_DEBUG").is_ok());

/// The strict audit's verdict: the filtered fill and the unfiltered fill
/// keep the same material.
///
/// Pieces are compared as merged parameter intervals per section: the
/// unfiltered run sees more curves, and a dropped pair's dead section can
/// still cross a live one and split its pieces differently, so individual
/// piece boundaries are noise and the kept *union* is the signal. Sections
/// are keyed by their face pair and a sampled point, because indices differ
/// between the runs.
fn audit_fill_equivalence(
    filtered: (&[SectionRec], &[SectionPiece]),
    unfiltered: (&[SectionRec], &[SectionPiece]),
    tol: Tolerances,
) {
    type Key = (usize, usize, [i64; 3]);
    let key_of = |section: &SectionRec, tol: Tolerances| -> Key {
        let (lo, hi) = section.curve.domain();
        let p = section
            .curve
            .point_at(f64::midpoint(lo, hi), tol)
            .unwrap_or(ogeom_math::Point::ORIGIN);
        let grid = (tol.confusion() * 100.0).max(1e-9);
        #[allow(clippy::cast_possible_truncation)]
        let q = |x: f64| (x / grid).round() as i64;
        (section.face_a, section.face_b, [q(p.x), q(p.y), q(p.z)])
    };
    let merged = |sections: &[SectionRec], pieces: &[SectionPiece], tol: Tolerances| {
        let mut kept: std::collections::HashMap<Key, Vec<(f64, f64)>> =
            std::collections::HashMap::new();
        for piece in pieces {
            kept.entry(key_of(&sections[piece.section], tol))
                .or_default()
                .push(piece.range);
        }
        let slop = tol.parametric() * 10.0;
        for ranges in kept.values_mut() {
            ranges.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(core::cmp::Ordering::Equal));
            let mut out: Vec<(f64, f64)> = Vec::with_capacity(ranges.len());
            for &(lo, hi) in ranges.iter() {
                match out.last_mut() {
                    Some(last) if lo <= last.1 + slop => last.1 = last.1.max(hi),
                    _ => out.push((lo, hi)),
                }
            }
            *ranges = out;
        }
        kept
    };
    let a = merged(filtered.0, filtered.1, tol);
    let b = merged(unfiltered.0, unfiltered.1, tol);
    let slop = tol.parametric() * 20.0;
    let matches = |x: &Vec<(f64, f64)>, y: &Vec<(f64, f64)>| {
        x.len() == y.len()
            && x.iter()
                .zip(y)
                .all(|(p, q)| (p.0 - q.0).abs() <= slop && (p.1 - q.1).abs() <= slop)
    };
    for (key, ranges) in &b {
        let held = a.get(key);
        assert!(
            held.is_some_and(|r| matches(r, ranges)),
            "strict bound-filter audit: the unfiltered fill keeps {ranges:?} \
             on section {key:?}, the filtered fill keeps {held:?}; the \
             filter drops material"
        );
    }
    for (key, ranges) in &a {
        assert!(
            b.contains_key(key),
            "strict bound-filter audit: the filtered fill keeps {ranges:?} \
             on section {key:?} the unfiltered fill never made"
        );
    }
}

fn general_fuse(model: &Model, a: &Shape, b: &Shape, tol: Tolerances) -> OgeomResult<GeneralFused> {
    ogeom_core::progress::stage("boolean: gather");
    let ga = gather(model, a, tol)?;
    let gb = gather(model, b, tol)?;
    ogeom_core::progress::stage("boolean: intersect");
    let (sections, section_pieces, contacts, tangents, contact_along, paves, same_a, same_b, hugs) =
        fill(&ga, &gb, false, tol)?;
    let mut junctions = pave_junctions(&ga, &gb, &paves, tol)?;
    if *DEBUG_WIRE {
        for j in &junctions {
            eprintln!("JUNCTION from paves at {:?} reach {:.3e}", j.at, j.reach);
        }
        for j in &hugs {
            eprintln!("JUNCTION from hugs at {:?} reach {:.3e}", j.at, j.reach);
        }
    }
    junctions.extend(hugs);
    for face in ga.faces.iter().chain(gb.faces.iter()) {
        // An input vertex that owns a span (the corner an earlier boolean
        // welded three rims into, each ending a fraction of a micron from
        // the others) is a junction here as well: the rebuild's positional
        // weld reaches only a few honesties, and minted separately those
        // ends stand as three vertices with hairlines between them.
        for e in &face.edges {
            for (at, radius) in e.ends {
                if radius > tol.confusion() * 10.0 {
                    if *DEBUG_WIRE {
                        eprintln!("JUNCTION from vertex at {at:?} radius {radius:.3e}");
                    }
                    junctions.push(Junction {
                        at,
                        reach: radius,
                        onto_vertex: false,
                    });
                }
            }
        }
    }
    // Junctions whose balls overlap are one junction: two pave clusters a
    // fraction of a micron apart at a corner where rims meet, each owning
    // a span that reaches into the other's, would weld a strand end to
    // whichever covers it first, and the corner would stand as two
    // vertices with a hairline between them, kept on one face, dropped
    // on another. Merged, the corner is one vertex whose reach covers
    // every member's span, on every face alike.
    junctions = merge_junctions(junctions);
    if *DEBUG_WIRE {
        for j in &junctions {
            eprintln!("JUNCTION at {:?} reach {:.3e}", j.at, j.reach);
        }
        let kind = |sf: &SurfaceGeometry| match sf {
            SurfaceGeometry::Plane(_) => "plane",
            SurfaceGeometry::Cylinder(_) => "cyl",
            SurfaceGeometry::BSpline(_) => "bspline",
            _ => "other",
        };
        for (si, sec) in sections.iter().enumerate() {
            let fa = &ga.faces[sec.face_a];
            let fb = &gb.faces[sec.face_b];
            let d = sec.curve.domain();
            let mut worst_a: f64 = 0.0;
            let mut worst_b: f64 = 0.0;
            let mut len = 0.0;
            let mut prev: Option<Point> = None;
            for i in 0..=32 {
                let t = if i == 32 {
                    d.1
                } else {
                    d.0 + (d.1 - d.0) * f64::from(i) / 32.0
                };
                let at = sec.curve.point_at(t, tol)?;
                if let Some(p) = prev {
                    len += p.distance(at);
                }
                prev = Some(at);
                let near = |f: &GFace| -> OgeomResult<f64> {
                    let mut best = f64::INFINITY;
                    for e in &f.edges {
                        best = best.min(distance_to_edge_curve(&e.curve, e.crange, at, tol)?);
                    }
                    Ok(best)
                };
                worst_a = worst_a.max(near(fa)?);
                worst_b = worst_b.max(near(fb)?);
            }
            let pieces = section_pieces.iter().filter(|p| p.section == si).count();
            eprintln!(
                "SECTION s{si} a{}({}) x b{}({}) tol {:.2e} closed {} len {:.4} pieces {pieces} hug-dist a {:.3e} b {:.3e}",
                sec.face_a,
                kind(&fa.surface),
                sec.face_b,
                kind(&fb.surface),
                sec.tolerance,
                sec.closed,
                len,
                worst_a,
                worst_b
            );
        }
        eprintln!("TANGENTS {}  CONTACTS {}", tangents.len(), contacts.len());
    }
    // The strict audit: fill again with every pair admitted and demand the
    // same kept material. The production result above is always the
    // filtered run (under audit too), so the audit compares rather than
    // substitutes, and a filtered run *is* replayed exactly. Zero cost with
    // the variable unset.
    if *AUDIT_BOUNDS {
        let (audit_sections, audit_pieces, ..) = fill(&ga, &gb, true, tol)?;
        audit_fill_equivalence(
            (&sections, &section_pieces),
            (&audit_sections, &audit_pieces),
            tol,
        );
    }

    ogeom_core::progress::stage("boolean: split");
    let mut pieces: Vec<FacePiece> = Vec::new();
    // Every face's strands, before any face is arranged: the dust decision
    // (a piece shorter than its face's snap collapses to one node) must
    // be one decision per shared edge piece, and a sub-piece of one edge
    // lies in several charts, each with its own snap and its own metric.
    // A band's rail hugging a wedge's cylinder for a quarter of a micron
    // is dust in the wall's chart and a strand in the cylinder's, and the
    // sew then finds the cylinder's piece used once.
    let strands_of = |from_a: bool,
                      fi: usize,
                      face: &GFace|
     -> OgeomResult<(Vec<Strand<Tag>>, f64)> {
        let mut strands: Vec<Strand<Tag>> = Vec::new();
        // How far apart the paves one junction stands for lie. A section
        // meeting the edge all but tangentially (its fitted end sliding
        // along it) paves it a few microns from where an exact section
        // crosses; the first pave speaks for both, and the other's strand
        // must still reach it.
        let mut spread = 0.0_f64;
        for (ei, e) in face.edges.iter().enumerate() {
            let mut stops = vec![e.crange.0];
            if let Some(ts) = paves.get(&e.node) {
                // Paves the edge itself cannot tell apart are one
                // junction, and the cluster's first pave speaks for it.
                for c in cluster_paves(&e.curve, e.crange, e.tolerance, ts, tol)? {
                    if c.members > 1 {
                        spread = spread.max(c.span + c.honesty);
                    }
                    stops.push(c.t);
                }
            }
            stops.push(e.crange.1);
            // A closed boundary edge (a cap's full circle) needs two
            // distinct endpoints per strand.
            let closed_edge = e
                .curve
                .point_at(e.crange.0, tol)?
                .distance(e.curve.point_at(e.crange.1, tol)?)
                <= tol.confusion();
            if closed_edge && stops.len() == 2 {
                stops.insert(1, f64::midpoint(e.crange.0, e.crange.1));
            }
            for pair in stops.windows(2) {
                let sub = (pair[0], pair[1]);
                if sub.1 - sub.0 <= tol.parametric() {
                    continue;
                }
                strands.push(Strand {
                    polyline: pcurve_polyline(
                        &e.pcurve,
                        e.prange,
                        e.crange,
                        sub,
                        &face.surface,
                        tol,
                    )?,
                    tag: Tag::Boundary {
                        edge: ei,
                        range: sub,
                    },
                    boundary: true,
                });
                if let Some((other_pc, orange)) = &e.other_side {
                    strands.push(Strand {
                        polyline: pcurve_polyline(
                            other_pc,
                            *orange,
                            e.crange,
                            sub,
                            &face.surface,
                            tol,
                        )?,
                        tag: Tag::Boundary {
                            edge: ei,
                            range: sub,
                        },
                        boundary: true,
                    });
                }
            }
        }
        for sp in &section_pieces {
            let section = &sections[sp.section];
            let (belongs, pcurve) = if from_a {
                (section.face_a == fi && !sp.hugs[0], &section.pc_a)
            } else {
                (section.face_b == fi && !sp.hugs[1], &section.pc_b)
            };
            if !belongs {
                continue;
            }
            let domain = section.curve.domain();
            let sub = if section.closed {
                (
                    fold(sp.range.0, domain),
                    sp.range.1 - sp.range.0 + fold(sp.range.0, domain),
                )
            } else {
                sp.range
            };
            // The pcurve shares the curve's parameterization; sampling
            // uses folded parameters for periodic curves.
            let count = 32;
            let mut line = Vec::with_capacity(count + 1);
            for i in 0..=count {
                #[allow(clippy::cast_precision_loss)]
                let t = sub.0 + (sub.1 - sub.0) * i as f64 / count as f64;
                let tf = if section.closed { fold(t, domain) } else { t };
                line.push(pcurve.point_at(tf, tol)?);
            }
            // Folding the parameter can tear the sampled polyline at the
            // period; unwrap it pointwise, then bring the whole strand
            // into the chart with one shift.
            unwrap_polyline(&mut line, &face.surface);
            {
                let trim: Vec<&[Point2]> = strands
                    .iter()
                    .filter(|st| st.boundary)
                    .map(|st| st.polyline.as_slice())
                    .collect();
                fold_line_inside(&mut line, &face.surface, &trim);
            }
            strands.push(Strand {
                polyline: line,
                tag: Tag::Section {
                    section: sp.section,
                    range: sp.range,
                },
                boundary: false,
            });
        }
        // The poles, after the sections, because a section can end *on*
        // a pole (a plane through a ball's axis cuts it exactly there),
        // and the pole has to be cut where that happens or the two meet
        // at no shared node and the arrangement sees a dangling section.
        for (pi, pole) in face.poles.iter().enumerate() {
            let mut stops = vec![pole.prange.0, pole.prange.1];
            if let PlanarCurve::Line(line) = &pole.pcurve {
                let axis = line.axis();
                for strand in &strands {
                    if strand.boundary {
                        continue;
                    }
                    for end in [strand.polyline.first(), strand.polyline.last()]
                        .into_iter()
                        .flatten()
                    {
                        let along = (*end - axis.location).dot(axis.direction.vector());
                        let foot = axis.point_at(along);
                        if foot.distance(*end) <= PARAM_SNAP
                            && along > pole.prange.0 + PARAM_SNAP
                            && along < pole.prange.1 - PARAM_SNAP
                        {
                            stops.push(along);
                        }
                    }
                }
            }
            stops.sort_by(|a, b| a.partial_cmp(b).unwrap_or(core::cmp::Ordering::Equal));
            stops.dedup_by(|a, b| (*a - *b).abs() <= PARAM_SNAP);
            for pair in stops.windows(2) {
                let sub = (pair[0], pair[1]);
                strands.push(Strand {
                    polyline: pcurve_polyline(
                        &pole.pcurve,
                        pole.prange,
                        pole.prange,
                        sub,
                        &face.surface,
                        tol,
                    )?,
                    tag: Tag::Pole {
                        pole: pi,
                        range: sub,
                    },
                    boundary: true,
                });
            }
        }
        for (ci, contact) in contacts.iter().enumerate() {
            if contact.target_from_a != from_a || contact.target_face != fi {
                continue;
            }
            // The owner's paves split its edge; the contact strands split
            // at the same parameters, so the sub-edges rebuilt from both
            // sides are the same edges and sew shared.
            let mut stops = vec![contact.crange.0];
            if let Some(ts) = paves.get(&contact.node) {
                stops.extend(
                    cluster_paves(&contact.curve, contact.crange, contact.tolerance, ts, tol)?
                        .iter()
                        .map(|c| c.t),
                );
            }
            stops.push(contact.crange.1);
            let closed_contact = contact
                .curve
                .point_at(contact.crange.0, tol)?
                .distance(contact.curve.point_at(contact.crange.1, tol)?)
                <= tol.confusion();
            if closed_contact && stops.len() == 2 {
                stops.insert(1, f64::midpoint(contact.crange.0, contact.crange.1));
            }
            for pair in stops.windows(2) {
                let sub = (pair[0], pair[1]);
                if sub.1 - sub.0 <= tol.parametric() {
                    continue;
                }
                let mid_t = f64::midpoint(sub.0, sub.1);
                if contact_along[ci]
                    .iter()
                    .any(|(lo, hi)| mid_t >= *lo && mid_t <= *hi)
                {
                    // Already boundary on both sides.
                    if *DEBUG_STRANDS {
                        eprintln!("CONTACT c{ci} {sub:?} on fi={fi}: along boundary, skipped");
                    }
                    continue;
                }
                // Keep only what lies inside this face's trim; the rest
                // of the owner's boundary splits nothing here.
                let mut line = pcurve_polyline(
                    &contact.pcurve,
                    contact.prange,
                    contact.crange,
                    sub,
                    &face.surface,
                    tol,
                )?;
                unwrap_polyline(&mut line, &face.surface);
                let boundary_lines: Vec<&[Point2]> = strands
                    .iter()
                    .filter(|st| st.boundary)
                    .map(|st| st.polyline.as_slice())
                    .collect();
                fold_line_inside(&mut line, &face.surface, &boundary_lines);
                let mid = interior_of(&line);
                if !inside_many_slanted(&boundary_lines, mid) {
                    if *DEBUG_STRANDS {
                        eprintln!(
                            "CONTACT c{ci} {sub:?} on fi={fi}: outside the trim at {mid:?}, skipped"
                        );
                    }
                    continue;
                }
                strands.push(Strand {
                    polyline: line,
                    tag: Tag::Contact {
                        contact: ci,
                        range: sub,
                    },
                    boundary: false,
                });
            }
        }

        // A tolerant contact's chart image meets the boundary it paved
        // only as closely as its own slop allows; the arrangement's node
        // weld reaches that far on this face, or the strand dangles a
        // few microns from the junction it belongs to.
        // A fitted section that hugs one of this face's edges leaves it
        // at the hug's far end by up to the hug's own width; its strand
        // must still find the edge's node there.
        let doubt_of = |edges: &[BoundaryEdge]| -> f64 {
            edges.iter().fold(0.0_f64, |acc, e| {
                acc.max(e.tolerance * 2.0)
                    .max(e.ends_tolerance + e.tolerance)
            })
        };
        let near = contacts
            .iter()
            .filter(|c| c.target_from_a == from_a && c.target_face == fi)
            .fold(PARAM_SNAP, |acc, c| acc.max(c.tolerance * 2.0))
            // An edge's end lands in the chart within its vertex's own
            // recorded doubt plus its own image's: a projected pcurve is
            // honest to the edge's tolerance, and the vertex it ends at
            // was welded to some earlier gap.
            .max(doubt_of(&face.edges))
            .max(spread)
            .max(
                if sections.iter().any(|s| {
                    s.tolerance > 0.0
                        && if from_a {
                            s.face_a == fi
                        } else {
                            s.face_b == fi
                        }
                }) {
                    tol.confusion() * 1e3
                } else {
                    0.0
                },
            );
        // A section ends where it crosses the other face's boundary,
        // and lands there within that boundary's own doubt: two
        // sections through neighbouring facets of a converted mesh
        // stop on the edge they share a few microns apart, each on its
        // own facet's plane, and meet on this face only if the weld
        // reaches that far.
        let far = |section: usize| -> f64 {
            let s = &sections[section];
            doubt_of(if from_a {
                &gb.faces[s.face_b].edges
            } else {
                &ga.faces[s.face_a].edges
            })
        };
        let face_snap = near.max(
            (0..sections.len())
                .filter(|&k| {
                    if from_a {
                        sections[k].face_a == fi
                    } else {
                        sections[k].face_b == fi
                    }
                })
                .map(far)
                .fold(0.0_f64, f64::max),
        );
        // The doubts above are lengths in space, and the weld runs in the
        // chart. Where the chart stretches every direction (a sphere of a
        // few millimetres runs a radian over its whole radius), a gap in
        // space is that many times narrower in the chart, and a weld taken
        // at the space length would swallow whole edges of a small patch.
        // Where some direction shrinks instead (a pole, a thin cylinder's
        // turn), the length stands.
        let stretch = strands
            .iter()
            .flat_map(|st| st.polyline.iter())
            .step_by(8)
            .filter_map(|p| face.surface.d1_at(p.x, p.y, tol).ok())
            .map(|(du, dv)| du.magnitude().min(dv.magnitude()))
            .fold(f64::INFINITY, f64::min);
        let face_snap = if stretch.is_finite() && stretch > 1.0 {
            (face_snap / stretch).max(PARAM_SNAP)
        } else {
            face_snap
        };
        // Where only one direction stretches (a cylinder's turn against its
        // straight length) no single scale holds, and a strand can be short
        // in the chart yet long in space. And the weld is the loosest doubt
        // of anything on the face: one sloppy edge a section crosses sets it
        // for every other section there. A strand longer in space than its
        // own doubt is no dust, so the weld stays below its chart length,
        // or it collapses and takes its neighbours' ends with it.
        let own_doubt = |tag: &Tag| -> f64 {
            match tag {
                Tag::Section { section, .. } => near.max(far(*section)),
                Tag::Pole { .. } => f64::INFINITY,
                Tag::Boundary { .. } | Tag::Contact { .. } => near,
            }
        };
        let face_snap = strands
            .iter()
            .filter(|st| {
                st.polyline.len() >= 2
                    && st
                        .polyline
                        .windows(2)
                        .map(|w| w[0].distance(w[1]))
                        .sum::<f64>()
                        <= face_snap
            })
            .filter_map(|st| {
                let points: Vec<Point> = st
                    .polyline
                    .iter()
                    .filter_map(|p| face.surface.point_at(p.x, p.y, tol).ok())
                    .collect();
                let space: f64 = points.windows(2).map(|w| w[0].distance(w[1])).sum();
                (space > own_doubt(&st.tag)).then(|| {
                    st.polyline
                        .windows(2)
                        .map(|w| w[0].distance(w[1]))
                        .sum::<f64>()
                        / 2.0
                })
            })
            .fold(face_snap, f64::min)
            .max(PARAM_SNAP);
        // A section that runs along one of this face's edges into its end
        // stops within the hug's width of the node the edge was split at,
        // and two sections meeting a rail that grazes this face stop on it
        // as far apart as the rail stays within tolerance of the face. The
        // junction there says each pair is one point. Where a chart unit is
        // shorter than a millimetre, or the graze long, that gap is wider in
        // the chart than the weld, and the section dangles beside the node it
        // ends at. Its end moves onto a boundary node the junction holds, or
        // onto the first section end already there.
        let mut anchors: Vec<(Option<usize>, Point2, Point)> = strands
            .iter()
            .filter(|st| st.boundary && st.polyline.len() >= 2)
            .flat_map(|st| [st.polyline[0], st.polyline[st.polyline.len() - 1]])
            .filter_map(|p| {
                face.surface
                    .point_at(p.x, p.y, tol)
                    .ok()
                    .map(|q| (None, p, q))
            })
            .collect();
        for st in strands
            .iter_mut()
            .filter(|st| matches!(st.tag, Tag::Section { .. }) && st.polyline.len() >= 2)
        {
            let last = st.polyline.len() - 1;
            for at in [0, last] {
                let p = st.polyline[at];
                if anchors
                    .iter()
                    .any(|(held, n, _)| held.is_none() && n.distance(p) <= face_snap)
                {
                    continue;
                }
                let Ok(q) = face.surface.point_at(p.x, p.y, tol) else {
                    continue;
                };
                let Some((ji, junction)) = junctions
                    .iter()
                    .enumerate()
                    .find(|(_, j)| j.at.distance(q) <= j.reach)
                else {
                    continue;
                };
                // The nearest in space the junction also holds, and reached
                // across the chart without leaving it: the chart's midpoint
                // of the two lies midway in space too, where a node in
                // another period's copy would not.
                let anchor = anchors
                    .iter()
                    .filter(|(held, _, s)| {
                        held.is_none_or(|h| h == ji) && junction.at.distance(*s) <= junction.reach
                    })
                    .filter(|(_, n, s)| {
                        face.surface
                            .point_at((n.x + p.x) / 2.0, (n.y + p.y) / 2.0, tol)
                            .is_ok_and(|m| m.distance(q.midpoint(*s)) <= junction.reach)
                    })
                    .min_by(|a, b| {
                        a.2.distance(q)
                            .partial_cmp(&b.2.distance(q))
                            .unwrap_or(core::cmp::Ordering::Equal)
                    })
                    .map(|(_, n, _)| *n);
                match anchor {
                    Some(n) => st.polyline[at] = n,
                    None => anchors.push((Some(ji), p, q)),
                }
            }
        }
        if *DEBUG_STRANDS {
            eprintln!(
                "FACE-SNAP from_a={from_a} fi={fi}: snap {face_snap:.3e} edges {:?}",
                face.edges
                    .iter()
                    .map(|e| (
                        format!("{:.2e}", e.tolerance),
                        format!("{:.2e}", e.ends_tolerance)
                    ))
                    .collect::<Vec<_>>()
            );
            for (si, st) in strands.iter().enumerate() {
                let tag = match st.tag {
                    Tag::Boundary { edge, range } => format!("Boundary e{edge} {range:?}"),
                    Tag::Contact { contact, range } => format!("Contact c{contact} {range:?}"),
                    Tag::Section { section, range } => format!("Section s{section} {range:?}"),
                    Tag::Pole { pole, range } => format!("Pole p{pole} {range:?}"),
                };
                let (a, b) = (st.polyline[0], st.polyline[st.polyline.len() - 1]);
                eprintln!(
                    "STRAND from_a={from_a} fi={fi} {si}: boundary={} {tag} pts={} {a:?} .. {b:?}",
                    st.boundary,
                    st.polyline.len()
                );
            }
        }
        Ok((strands, face_snap))
    };
    let chart_length =
        |line: &[Point2]| -> f64 { line.windows(2).map(|w| w[0].distance(w[1])).sum::<f64>() };
    // The identity of a piece across faces: a boundary or contact piece by
    // its edge node, a section piece by its section: the same section is
    // split the same way on both faces it cuts, and a piece that is dust
    // in a sphere's chart at the pole must be dust on the plane it also
    // lies in.
    let node_of = |face: &GFace, tag: &Tag| -> Option<(usize, usize, (f64, f64))> {
        match tag {
            Tag::Boundary { edge, range } => {
                Some((0, face.edges[*edge].node.index() as usize, *range))
            }
            Tag::Contact { contact, range } => {
                Some((1, contacts[*contact].node.index() as usize, *range))
            }
            Tag::Section { section, range } => Some((2, *section, *range)),
            Tag::Pole { .. } => None,
        }
    };
    /// A face's strands and its snap, waiting to be arranged.
    type Prepared = Option<(Vec<Strand<Tag>>, f64)>;
    let mut prepared: [Vec<Prepared>; 2] = [Vec::new(), Vec::new()];
    let mut dust: Vec<(usize, usize, (f64, f64))> = Vec::new();
    for (side, solid, from_a) in [(0_usize, &ga, true), (1, &gb, false)] {
        for (fi, face) in solid.faces.iter().enumerate() {
            ogeom_core::progress::checkpoint()?;
            let (strands, snap) = strands_of(from_a, fi, face)?;
            for st in &strands {
                if st.polyline.len() >= 2
                    && chart_length(&st.polyline) <= snap
                    && let Some(key) = node_of(face, &st.tag)
                {
                    dust.push(key);
                }
            }
            prepared[side].push(Some((strands, snap)));
        }
    }
    ogeom_core::progress::stage("boolean: split arrange");
    let same_key = |a: &(usize, usize, (f64, f64)), b: &(usize, usize, (f64, f64))| {
        a.0 == b.0
            && a.1 == b.1
            && (a.2.0 - b.2.0).abs() <= tol.parametric()
            && (a.2.1 - b.2.1).abs() <= tol.parametric()
    };
    for (from_a, own, other) in [(true, &ga, &gb.solid), (false, &gb, &ga.solid)] {
        // The other solid's boundary, prepared once for the whole side. It is
        // asked once per face piece, and what it costs to prepare (every
        // face's trimming rings, polylined) does not depend on the point
        // being asked about. Rebuilt per question it dwarfed the question:
        // 3.5 ms of preparation against 5.6 µs of ray casting.
        let boundary = ogeom_algo::SolidBoundary::of(model, other, tol.confusion() * 1e4, tol)?;
        for (fi, face) in own.faces.iter().enumerate() {
            ogeom_core::progress::checkpoint()?;
            let Some((mut strands, face_snap)) = prepared[usize::from(!from_a)][fi].take() else {
                ogeom_bail!(Construction, "a face was prepared twice");
            };
            // A piece some other chart already collapsed collapses here too:
            // its ends become one node, every neighbour meeting them moves
            // onto it, and one junction owns the span in space so the
            // rebuilt vertices agree on every face.
            // A strand's ends in space, or none for a pole.
            let space_ends = |face: &GFace, tag: &Tag| -> OgeomResult<Option<(Point, Point)>> {
                Ok(Some(match tag {
                    Tag::Boundary { edge, range } => {
                        let e = &face.edges[*edge];
                        (
                            e.curve.point_at(range.0, tol)?,
                            e.curve.point_at(range.1, tol)?,
                        )
                    }
                    Tag::Contact { contact, range } => {
                        let c = &contacts[*contact];
                        (
                            c.curve.point_at(range.0, tol)?,
                            c.curve.point_at(range.1, tol)?,
                        )
                    }
                    Tag::Section { section, range } => {
                        let sec = &sections[*section];
                        let domain = sec.curve.domain();
                        (
                            sec.curve
                                .point_at(at_param(range.0, domain, sec.closed), tol)?,
                            sec.curve
                                .point_at(at_param(range.1, domain, sec.closed), tol)?,
                        )
                    }
                    Tag::Pole { .. } => return Ok(None),
                }))
            };
            // A strand this face's own weld collapses is one point in space
            // as well: its ends become one junction, or the pieces either
            // side of it end on two vertices the chart says are one, and a
            // loop cut at a seam a few microns from its own start closes in
            // the chart and stays open in space.
            for st in &strands {
                if st.polyline.len() >= 2
                    && chart_length(&st.polyline) <= face_snap
                    && let Some((from, to)) = space_ends(face, &st.tag)?
                    && from.distance(to) > tol.confusion() * 1e2
                {
                    junctions.push(Junction {
                        at: from.midpoint(to),
                        reach: from.distance(to) / 2.0 + tol.confusion() * 1e2,
                        onto_vertex: false,
                    });
                }
            }
            let forced: Vec<usize> = strands
                .iter()
                .enumerate()
                .filter(|(_, st)| {
                    st.polyline.len() >= 2
                        && chart_length(&st.polyline) > face_snap
                        && node_of(face, &st.tag)
                            .is_some_and(|k| dust.iter().any(|d| same_key(d, &k)))
                })
                .map(|(i, _)| i)
                .collect();
            if !forced.is_empty() {
                let mut reps: Vec<(Point2, usize)> = Vec::new();
                let canon = |p: Point2, reps: &mut Vec<(Point2, usize)>| -> usize {
                    match reps.iter().position(|(q, _)| q.distance(p) <= face_snap) {
                        Some(i) => i,
                        None => {
                            reps.push((p, reps.len()));
                            reps.len() - 1
                        }
                    }
                };
                fn root(reps: &mut [(Point2, usize)], mut i: usize) -> usize {
                    while reps[i].1 != i {
                        reps[i].1 = reps[reps[i].1].1;
                        i = reps[i].1;
                    }
                    i
                }
                for &i in &forced {
                    let st = &strands[i];
                    let a = canon(st.polyline[0], &mut reps);
                    let b = canon(st.polyline[st.polyline.len() - 1], &mut reps);
                    let (ra, rb) = (root(&mut reps, a), root(&mut reps, b));
                    if ra != rb {
                        reps[rb].1 = ra;
                    }
                    let Some((from, to)) = space_ends(face, &st.tag)? else {
                        continue;
                    };
                    if *DEBUG_STRANDS {
                        eprintln!(
                            "DUST from_a={from_a} fi={fi}: strand {i} collapses with its partners, {from:?} -> {to:?}"
                        );
                    }
                    junctions.push(Junction {
                        at: from.midpoint(to),
                        reach: from.distance(to) / 2.0 + tol.confusion() * 1e2,
                        onto_vertex: false,
                    });
                }
                for (i, st) in strands.iter_mut().enumerate() {
                    if forced.contains(&i) {
                        continue;
                    }
                    let last = st.polyline.len() - 1;
                    for end in [0, last] {
                        if let Some(k) = reps
                            .iter()
                            .position(|(q, _)| q.distance(st.polyline[end]) <= face_snap)
                        {
                            let r = root(&mut reps, k);
                            st.polyline[end] = reps[r].0;
                        }
                    }
                }
                let mut index = 0_usize;
                strands.retain(|_| {
                    let keep = !forced.contains(&index);
                    index += 1;
                    keep
                });
            }
            let split = match arrange_pieces(&strands, face_snap) {
                Ok(split) => split,
                Err(err) => {
                    if *ARRANGE_DEBUG {
                        eprintln!("ARRANGE from_a={from_a} fi={fi} failed: {err}");
                        for (si, st) in strands.iter().enumerate() {
                            let (a, b) = (st.polyline[0], st.polyline[st.polyline.len() - 1]);
                            let describe = |curve: &Curve, range: (f64, f64)| -> String {
                                let mid = curve
                                    .point_at(f64::midpoint(range.0, range.1), tol)
                                    .map(|p| format!("({:.4},{:.4},{:.4})", p.x, p.y, p.z))
                                    .unwrap_or_default();
                                match curve {
                                    Curve::Circle(c) => format!(
                                        "circle centre {:?} r {:.4} mid {mid}",
                                        c.circle().centre(),
                                        c.circle().radius()
                                    ),
                                    other => {
                                        format!("{:?} mid {mid}", core::mem::discriminant(other))
                                    }
                                }
                            };
                            let geometry = match st.tag {
                                Tag::Boundary { edge, range } => {
                                    describe(&face.edges[edge].curve, range)
                                }
                                Tag::Contact { contact, range } => {
                                    format!(
                                        "{} window {:?} along {:?} node {}",
                                        describe(&contacts[contact].curve, range),
                                        contacts[contact].crange,
                                        contact_along[contact],
                                        contacts[contact].node.index()
                                    )
                                }
                                _ => String::new(),
                            };
                            eprintln!("    {geometry}");
                            eprintln!(
                                "  strand {si} boundary={} {:?} {a:?} .. {b:?}",
                                st.boundary,
                                match st.tag {
                                    Tag::Boundary { edge, range } =>
                                        format!("Boundary e{edge} {range:?}"),
                                    Tag::Contact { contact, range } =>
                                        format!("Contact c{contact} {range:?}"),
                                    Tag::Section { section, range } =>
                                        format!("Section s{section} {range:?}"),
                                    Tag::Pole { pole, range } => format!("Pole p{pole} {range:?}"),
                                }
                            );
                        }
                    }
                    return Err(err);
                }
            };
            for piece in split {
                // Where a piece stands is asked at its interior probes in
                // turn. The first is the roomiest, and usually the only one
                // needed; the rest are for the piece that merely *touches*
                // the other solid, whose roomiest probe can land on the
                // contact and read neither in nor out.
                //
                // A partner face is asked before the classifier, not after.
                // Where a piece sits on a surface the other solid also
                // carries, the partner *is* the answer, and getting there
                // through the classifier means every ray grazing the shared
                // face, its whole fan of directions exhausted, before it
                // reports the On the partner list already knew. On a part
                // whose bore is refilled by its own cylinder that is the
                // difference between a tenth of a second and a minute.
                let partners = if from_a { &same_a[fi] } else { &same_b[fi] };
                let mut chosen = None;
                for candidate in &piece.interiors {
                    let at = face.surface.point_at(candidate.x, candidate.y, tol)?;
                    let shared = !partners.is_empty()
                        && partners.iter().any(|&pi| {
                            let partner = if from_a { &gb.faces[pi] } else { &ga.faces[pi] };
                            chart_point_of(partner, at, tol).is_some()
                        });
                    let says = if shared {
                        Containment::On
                    } else {
                        boundary.holds(model, at, tol)?
                    };
                    if chosen.is_none() || !matches!(says, Containment::On) {
                        chosen = Some((*candidate, at, says));
                    }
                    if !matches!(says, Containment::On) {
                        break;
                    }
                }
                let Some((interior, probe, said)) = chosen else {
                    ogeom_bail!(
                        Construction,
                        "a piece of a face has no interior point to classify at"
                    );
                };
                let state = match said {
                    Containment::In => PieceState::In,
                    Containment::Out => PieceState::Out,
                    Containment::On => {
                        // On the other boundary: same-domain contact. The
                        // partner face on the shared surface decides whether
                        // the two materials lie on the same side or oppose.
                        let own_normal = outward_normal(face, interior, tol)?;
                        let mut resolved = None;
                        for &pi in partners {
                            let partner = if from_a { &gb.faces[pi] } else { &ga.faces[pi] };
                            let Some(uv) = chart_point_of(partner, probe, tol) else {
                                continue;
                            };
                            let theirs = outward_normal(partner, uv, tol)?;
                            resolved = Some(if own_normal.dot(theirs) > 0.0 {
                                PieceState::OnAligned
                            } else {
                                PieceState::OnOpposed
                            });
                            break;
                        }
                        let Some(state) = resolved else {
                            // On with no partner containing the probe: the
                            // band's generosity read proximity as
                            // coincidence. Ask again at a width where it
                            // cannot, and only a genuine edge contact
                            // remains refused.
                            match ogeom_algo::classify_in_solid_exact_banded(
                                model,
                                other,
                                probe,
                                tol.confusion() * 10.0,
                                tol,
                            )? {
                                Containment::In => {
                                    pieces.push(FacePiece {
                                        from_a,
                                        face: fi,
                                        rings: piece.rings,
                                        outlines: piece.outlines,
                                        probe,
                                        state: PieceState::In,
                                        covered: false,
                                    });
                                    continue;
                                }
                                Containment::Out => {
                                    pieces.push(FacePiece {
                                        from_a,
                                        face: fi,
                                        rings: piece.rings,
                                        outlines: piece.outlines,
                                        probe,
                                        state: PieceState::Out,
                                        covered: false,
                                    });
                                    continue;
                                }
                                Containment::On => {}
                            }
                            if *DEBUG_WIRE {
                                eprintln!(
                                    "ON-CONTACT piece of {} face {fi} probe {probe:?} partners {partners:?}",
                                    if from_a { "A" } else { "B" }
                                );
                            }
                            ogeom_bail!(
                                NotDone,
                                "a piece lies on the other solid's boundary \
                                 with no coincident partner face to compare \
                                 sides against; edge or vertex contact is \
                                 refused rather than resolved; see the \
                                 remaining work in docs/PLAN.md"
                            );
                        };
                        state
                    }
                };
                pieces.push(FacePiece {
                    from_a,
                    face: fi,
                    rings: piece.rings,
                    outlines: piece.outlines,
                    probe,
                    state,
                    covered: false,
                });
            }
        }
    }
    mark_covered_coincidences(&ga, &mut pieces, tol);
    if *ARRANGE_DEBUG {
        for (fi, partners) in same_a.iter().enumerate() {
            if !partners.is_empty() {
                eprintln!("SAME a{fi} with b{partners:?}");
            }
        }
        for (ci, contact) in contacts.iter().enumerate() {
            eprintln!(
                "CONTACT c{ci} onto {}{} over {:?} prange {:?} node {} pcurve {:?}",
                if contact.target_from_a { "a" } else { "b" },
                contact.target_face,
                contact.crange,
                contact.prange,
                contact.node.index(),
                contact.pcurve
            );
        }
        for (i, p) in pieces.iter().enumerate() {
            let own = if p.from_a {
                &ga.faces[p.face]
            } else {
                &gb.faces[p.face]
            };
            let kind = match &own.surface {
                SurfaceGeometry::Plane(_) => "plane",
                SurfaceGeometry::Cylinder(_) => "cyl",
                SurfaceGeometry::Torus(_) => "torus",
                SurfaceGeometry::BSpline(_) => "bspline",
                _ => "other",
            };
            eprintln!(
                "PIECE {i} from_a={} face={} {kind} state={:?} covered={} probe={:?} rings {:?}",
                p.from_a,
                p.face,
                p.state,
                p.covered,
                p.probe,
                p.rings.iter().map(Vec::len).collect::<Vec<_>>()
            );
        }
    }
    ogeom_core::progress::stage("boolean: classified");
    // Merged again: the arrangement adds junctions of its own (a strand
    // that collapsed with its partners becomes one) after the first
    // merge, and two of those at one triple point, each within the other's
    // reach, welded the ends of one band to two vertices with a hairline
    // between them.
    let junctions = merge_junctions(junctions);
    Ok(GeneralFused {
        a: ga,
        b: gb,
        sections,
        contacts,
        tangents,
        pieces,
        junctions,
    })
}

// --- rebuilding --------------------------------------------------------------

/// Everything the rebuild shares across pieces.
struct Rebuild<'m> {
    model: &'m mut Model,
    /// World surface ids, minted once per source face.
    surfaces_a: Vec<Option<ogeom_topo::SurfaceId>>,
    surfaces_b: Vec<Option<ogeom_topo::SurfaceId>>,
    /// Vertices shared by position: a wire's connectivity is checked by node
    /// identity, so two sub-edges meeting at a point must *name* the same
    /// vertex, not merely coincide there.
    vertices: Vec<(Point, Shape)>,
    /// `vertices` binned by position.
    vertex_bins: bins::Bins,
    /// The widest tolerance any vertex in `vertices` holds: every widening
    /// of one goes through [`Rebuild::vertex`] or [`Rebuild::widen`], so an
    /// end is compared only with the vertices this far and a weld from it.
    widest: f64,
    /// How far two honest descriptions of one junction may sit apart: a
    /// hundred confusions as the floor, widened to three times the loosest
    /// contact edge's or fitted section's own tolerance when one took part:
    /// the strand's own, the crossing's, and the edge's it stopped on.
    weld: f64,
    /// The paving's junctions, each minted as a vertex the first time a
    /// strand end lands inside it.
    junctions: Vec<(Junction, Option<Shape>)>,
    /// `junctions` binned by position.
    junction_bins: bins::Bins,
    /// The widest reach of any junction.
    junction_reach: f64,
    /// Vertices minted from junctions at an edge's own end vertex, which a
    /// tangent section's chart image is bent onto.
    onto_vertex: std::collections::HashSet<ogeom_topo::TShapeId>,
}

impl Rebuild<'_> {
    fn vertex(&mut self, p: Point, tol: Tolerances) -> Shape {
        // A junction several paves described is one vertex owning the span
        // they disagree by, and every strand end inside it names that
        // vertex, whichever pave its own trim stopped at.
        // The rim is as generous as the positional weld below: two strand
        // ends a weld apart must not fall on opposite sides of it.
        let slack = self.weld.max(tol.confusion() * 1e2);
        if *DEBUG_WIRE {
            eprintln!("VERTEX ask {p:?}");
        }
        let reach = self.junction_reach;
        let inside = |slot: &usize| {
            let j = &self.junctions[*slot].0;
            j.at.distance(p) <= j.reach + slack
        };
        let first = match self.junction_bins.near(p, reach + slack) {
            Some(near) => near.into_iter().find(inside),
            None => (0..self.junctions.len()).find(inside),
        };
        if let Some(slot) = first {
            let junction = self.junctions[slot].0;
            let shape = match &self.junctions[slot].1 {
                Some(shape) => shape.clone(),
                None => {
                    let shape = make_vertex(self.model, junction.at).shape;
                    self.junctions[slot].1 = Some(shape.clone());
                    self.remember(junction.at, &shape);
                    shape
                }
            };
            if junction.onto_vertex {
                self.onto_vertex.insert(shape.node());
            }
            // The junction owns its span, and an end welded in from the rim
            // widens it by what it actually sat off by, as any weld does.
            let gap = junction.at.distance(p);
            if let Some(node) = self.model.node_mut(&shape)
                && let ogeom_topo::NodeData::Vertex(data) = node.data_mut()
            {
                // A confusion of slack past the measured gap: an edge built
                // against this vertex measures the same gap through its own
                // rounding, and a tolerance equal to it fails by an ulp.
                data.tolerance = data
                    .tolerance
                    .widen_to(junction.reach.max(gap) + tol.confusion());
            }
            self.note_tolerance(&shape);
            return shape;
        }
        // The weld reach covers what the inputs may honestly disagree by: a
        // boundary curve that is itself a fitted intersection from an
        // earlier boolean carries a couple of microns of slop, and two
        // descriptions of one junction arrive that far apart. A hundred
        // confusions at millimetre tolerances is ten microns, below any
        // feature this pipeline can resolve, and every weld wider than
        // confusion is recorded on the vertex, not papered over.
        // And as far again as the vertex already owns: a vertex widened by
        // the descriptions it has taken in claims its point within that
        // tolerance, an end arriving within that and the weld of it may be
        // the same point, and a boolean cannot show it is not. Three
        // descriptions of one triple point, where a third band met two
        // bands at an apex, arrived a tenth of a micron apart in turn, each
        // within the last's reach and none within the first's, and the
        // wire round the band's end had two vertices where it needed one.
        // Every end taken in is remembered where it arrived, so the next
        // end is measured from the description nearest it.
        let floor = self.weld.max(tol.confusion() * 1e2);
        let near = self
            .vertex_bins
            .near(p, floor + self.widest)
            .unwrap_or_else(|| (0..self.vertices.len()).collect());
        let found = near
            .into_iter()
            .map(|index| &self.vertices[index])
            .filter_map(|(q, shape)| {
                let own = self
                    .model
                    .node(shape)
                    .and_then(|n| n.data().as_vertex())
                    .map_or(0.0, |d| d.tolerance.get());
                let gap = q.distance(p);
                (gap <= floor.max(own + floor)).then_some((gap, shape.clone()))
            })
            .min_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(core::cmp::Ordering::Equal));
        if let Some((gap, shape)) = found {
            // Two descriptions of one junction may disagree by a general
            // crossing's residual; the vertex's tolerance is where that
            // disagreement is recorded, so the sub-edges built against
            // either description still reach it honestly.
            let at = self
                .model
                .node(&shape)
                .and_then(|n| n.data().as_vertex())
                .map_or(p, |d| d.point);
            let off = at.distance(p);
            if off > tol.confusion()
                && let Some(node) = self.model.node_mut(&shape)
                && let ogeom_topo::NodeData::Vertex(data) = node.data_mut()
            {
                data.tolerance = data.tolerance.widen_to(off + tol.confusion());
            }
            self.note_tolerance(&shape);
            if gap > tol.confusion() {
                self.remember(p, &shape);
            }
            return shape;
        }
        let shape = make_vertex(self.model, p).shape;
        self.remember(p, &shape);
        shape
    }

    /// Take in `shape` as described at `p`.
    fn remember(&mut self, p: Point, shape: &Shape) {
        self.vertex_bins.insert(p, self.vertices.len());
        self.vertices.push((p, shape.clone()));
        self.note_tolerance(shape);
    }

    /// Keep `widest` covering `shape`'s tolerance.
    fn note_tolerance(&mut self, shape: &Shape) {
        if let Some(data) = self.model.node(shape).and_then(|n| n.data().as_vertex()) {
            self.widest = self.widest.max(data.tolerance.get());
        }
    }

    /// Widen a vertex this rebuild handed out.
    fn widen(&mut self, vertex: &Shape, to: f64) -> OgeomResult<()> {
        self.model.widen(vertex, ogeom_core::Tolerance::new(to)?)?;
        self.note_tolerance(vertex);
        Ok(())
    }

    fn surface_id(
        &mut self,
        fused: &GeneralFused,
        from_a: bool,
        face: usize,
    ) -> ogeom_topo::SurfaceId {
        let slot = if from_a {
            &mut self.surfaces_a[face]
        } else {
            &mut self.surfaces_b[face]
        };
        if let Some(id) = slot {
            return *id;
        }
        let surface = if from_a {
            fused.a.faces[face].surface.clone()
        } else {
            fused.b.faces[face].surface.clone()
        };
        let id = self.model.geometry_mut().add_surface(surface);
        *slot = Some(id);
        id
    }
}

/// Build one piece as a face, orientation matching its source face's side.
fn build_piece(
    rebuild: &mut Rebuild,
    fused: &GeneralFused,
    piece: &FacePiece,
    tol: Tolerances,
) -> OgeomResult<Shape> {
    let own = if piece.from_a { &fused.a } else { &fused.b };
    let face = &own.faces[piece.face];
    let surface_id = rebuild.surface_id(fused, piece.from_a, piece.face);

    // Sub-edges cached within the piece, so a seam used from both sides is
    // one edge appearing twice.
    let mut cache: Vec<(usize, u8, (f64, f64), Shape)> = Vec::new();
    let mut wires = Vec::new();
    for ring in &piece.rings {
        let mut edges = Vec::with_capacity(ring.len());
        for traversal in ring {
            let (key_edge, key_kind, range) = match &traversal.tag {
                Tag::Boundary { edge, range } => (*edge, 0_u8, *range),
                Tag::Section { section, range } => (*section, 1, *range),
                Tag::Contact { contact, range } => (*contact, 2, *range),
                Tag::Pole { pole, range } => (*pole, 3, *range),
            };
            let near = |a: (f64, f64), b: (f64, f64)| {
                (a.0 - b.0).abs() <= tol.parametric() && (a.1 - b.1).abs() <= tol.parametric()
            };
            if *DEBUG_WIRE && (range.1 - range.0).abs() < 1e-4 {
                let kind = match &traversal.tag {
                    Tag::Boundary { .. } => "boundary",
                    Tag::Section { .. } => "section",
                    Tag::Contact { .. } => "contact",
                    Tag::Pole { .. } => "pole",
                };
                eprintln!(
                    "DUST {kind} {key_edge} range {range:?} in piece from_a={} face={}",
                    piece.from_a, piece.face
                );
            }
            let built = if let Some((.., shape)) = cache
                .iter()
                .find(|(k, s, r, _)| *k == key_edge && *s == key_kind && near(*r, range))
            {
                shape.clone()
            } else {
                let shape = build_sub_edge(
                    rebuild,
                    fused,
                    piece.from_a,
                    face,
                    surface_id,
                    &traversal.tag,
                    tol,
                )?;
                cache.push((key_edge, key_kind, range, shape.clone()));
                shape
            };
            // A piece whose two ends welded to one vertex and whose whole
            // length lies within the vertex's reach is that vertex's dust
            // (a rim's last fraction of a micron before a pole corner) and
            // no edge of the wire; a closed edge on one vertex, a full
            // circle, reaches far from it and stays.
            if !matches!(traversal.tag, Tag::Pole { .. })
                && let Some((v0, v1)) = ogeom_algo::edge_vertices(rebuild.model, &built)?
                && v0.node() == v1.node()
                && let Some(vd) = rebuild
                    .model
                    .node(&v0)
                    .and_then(|n| n.data().as_vertex())
                    .map(|d| (d.point, d.tolerance.get()))
                && let Some(ed) = rebuild.model.node(&built).and_then(|n| n.data().as_edge())
                && let Some(ogeom_topo::EdgeRepr::Curve3d {
                    curve, range: r3, ..
                }) = ed.curve3d()
                && let Some(g) = rebuild.model.geometry().curve(*curve)
                && g.point_at(f64::midpoint(r3.0, r3.1), tol)?.distance(vd.0)
                    <= vd.1.max(tol.confusion() * 1e4)
            {
                if *DEBUG_WIRE {
                    eprintln!(
                        "DUST piece closing on one vertex at {:?} dropped from piece from_a={} face={}",
                        vd.0, piece.from_a, piece.face
                    );
                }
                continue;
            }
            edges.push(if traversal.reversed {
                built.reversed()
            } else {
                built
            });
        }
        let wire = match make_wire(rebuild.model, &edges, tol) {
            Ok(w) => w.shape,
            Err(e) => {
                if *DEBUG_WIRE {
                    eprintln!(
                        "WIRE FAIL piece from_a={} face={}",
                        piece.from_a, piece.face
                    );
                    for built in &edges {
                        if let Some((a, b)) = ogeom_algo::edge_vertices(rebuild.model, built)? {
                            let at = |v: &Shape| {
                                rebuild
                                    .model
                                    .node(v)
                                    .and_then(|n| n.data().as_vertex())
                                    .map(|d| (d.point, d.tolerance.get()))
                            };
                            eprintln!(
                                "   built {:?}{}: v{} {:?} -> v{} {:?}",
                                built.node(),
                                if built.orientation() == ogeom_topo::Orientation::Reversed {
                                    " rev"
                                } else {
                                    ""
                                },
                                a.node().index(),
                                at(&a),
                                b.node().index(),
                                at(&b)
                            );
                        }
                    }
                    for t in ring {
                        let tag = match &t.tag {
                            Tag::Boundary { edge, range } => format!(
                                "Boundary e{edge} {range:?} tol {:.2e}",
                                face.edges[*edge].tolerance
                            ),
                            Tag::Contact { contact, range } => format!(
                                "Contact c{contact} {range:?} tol {:.2e}",
                                fused.contacts[*contact].tolerance
                            ),
                            Tag::Section { section, range } => format!(
                                "Section s{section} {range:?} tol {:.2e}",
                                fused.sections[*section].tolerance
                            ),
                            Tag::Pole { pole, range } => format!("Pole p{pole} {range:?}"),
                        };
                        eprintln!("   {tag} reversed={}", t.reversed);
                    }
                }
                return Err(e);
            }
        };
        wires.push(wire);
    }
    let built = match make_face_on(rebuild.model, surface_id, &wires, tol) {
        Ok(b) => b.shape,
        Err(e) => {
            if *DEBUG_WIRE {
                eprintln!(
                    "FACE FAIL piece from_a={} face={}: {e}",
                    piece.from_a, piece.face
                );
                for (wi, wire) in wires.iter().enumerate() {
                    for edge in rebuild.model.ordered_children_of(wire)? {
                        if let Some((a, b)) = ogeom_algo::edge_vertices(rebuild.model, &edge)? {
                            let at = |v: &Shape| {
                                rebuild
                                    .model
                                    .node(v)
                                    .and_then(|n| n.data().as_vertex())
                                    .map(|d| d.point)
                            };
                            eprintln!(
                                "   wire {wi} edge {:?}{}: v{} {:?} -> v{} {:?}",
                                edge.node(),
                                if edge.orientation() == ogeom_topo::Orientation::Reversed {
                                    " rev"
                                } else {
                                    ""
                                },
                                a.node().index(),
                                at(&a),
                                b.node().index(),
                                at(&b)
                            );
                        }
                    }
                }
            }
            return Err(e);
        }
    };
    Ok(
        if face.face.orientation() == ogeom_topo::Orientation::Reversed {
            built.reversed()
        } else {
            built
        },
    )
}

/// Build the exact sub-edge a tag names, pcurves attached.
fn build_sub_edge(
    rebuild: &mut Rebuild,
    fused: &GeneralFused,
    from_a: bool,
    face: &GFace,
    surface_id: ogeom_topo::SurfaceId,
    tag: &Tag,
    tol: Tolerances,
) -> OgeomResult<Shape> {
    match tag {
        Tag::Pole { pole, range } => {
            // A pole rebuilds as what it was: one vertex, an edge with no
            // curve bounded by it twice, and the chart line that says where
            // it runs in this face's parameters.
            let p = &face.poles[*pole];
            let at = rebuild.vertex(p.point, tol);
            let model = &mut *rebuild.model;
            let mut data = ogeom_topo::EdgeData::new();
            data.degenerate = true;
            let built = model.add_edge(data, &[at.clone(), at])?;
            ogeom_algo::attach_pcurve(
                model,
                &built,
                p.pcurve.clone(),
                surface_id,
                Location::identity(),
                *range,
            )?;
            Ok(built)
        }
        Tag::Boundary { edge, range } => {
            let e = &face.edges[*edge];
            let from = e.curve.point_at(range.0, tol)?;
            let to = e.curve.point_at(range.1, tol)?;
            let v0 = rebuild.vertex(from, tol);
            let v1 = rebuild.vertex(to, tol);
            let model = &mut *rebuild.model;
            let built = make_edge_between(model, e.curve.clone(), *range, &v0, &v1, tol)?.shape;
            // A piece of a tolerant edge is the same curve with the same
            // honest radius: a fitted rail's stated slop must survive the
            // split, or the next boolean over this solid measures the rail
            // against a tolerance it never had. The ends own it too: a
            // junction on a tolerant curve is a junction to the curve's own
            // resolution, and the wires rebuilt through it meet within that.
            if e.tolerance > tol.confusion() {
                if let Some(node) = model.node_mut(&built)
                    && let ogeom_topo::NodeData::Edge(data) = node.data_mut()
                {
                    data.tolerance = data.tolerance.widen_to(e.tolerance);
                }
                for v in [&v0, &v1] {
                    rebuild.widen(v, e.tolerance)?;
                }
            }
            let model = &mut *rebuild.model;
            let sub_p = (
                rescale(range.0, e.crange, e.prange),
                rescale(range.1, e.crange, e.prange),
            );
            match &e.other_side {
                None => ogeom_algo::attach_pcurve(
                    model,
                    &built,
                    e.pcurve.clone(),
                    surface_id,
                    Location::identity(),
                    sub_p,
                )?,
                Some((other, orange)) => {
                    // A seam: both sides attach, and an occurrence picks its
                    // side by its orientation.
                    let sub_o = (
                        rescale(range.0, e.crange, *orange),
                        rescale(range.1, e.crange, *orange),
                    );
                    let _ = sub_o;
                    ogeom_algo::attach_seam(
                        model,
                        &built,
                        e.pcurve.clone(),
                        other.clone(),
                        surface_id,
                        Location::identity(),
                        sub_p,
                    )?;
                }
            }
            Ok(built)
        }
        Tag::Contact { contact, range } => {
            let c = &fused.contacts[*contact];
            let from = c.curve.point_at(range.0, tol)?;
            let to = c.curve.point_at(range.1, tol)?;
            let v0 = rebuild.vertex(from, tol);
            let v1 = rebuild.vertex(to, tol);
            let model = &mut *rebuild.model;
            let built = make_edge_between(model, c.curve.clone(), *range, &v0, &v1, tol)?.shape;
            // A tolerant contact's pieces and ends own its stated slop, as a
            // boundary's and a section's do: three kinds of strand close one
            // ring, and every junction meets at the strands' own honesty.
            if c.tolerance > tol.confusion() {
                if let Some(node) = model.node_mut(&built)
                    && let ogeom_topo::NodeData::Edge(data) = node.data_mut()
                {
                    data.tolerance = data.tolerance.widen_to(c.tolerance);
                }
                for v in [&v0, &v1] {
                    rebuild.widen(v, c.tolerance)?;
                }
            }
            let model = &mut *rebuild.model;
            // The stored image keeps its own window; the attached copy names
            // the sub-window this piece covers under the proportional map.
            let sub_p = (
                rescale(range.0, c.crange, c.prange),
                rescale(range.1, c.crange, c.prange),
            );
            let mid = c.pcurve.point_at(f64::midpoint(sub_p.0, sub_p.1), tol)?;
            let trim = face_trim_lines(face, tol);
            let trim: Vec<&[Point2]> = trim.iter().map(Vec::as_slice).collect();
            let folded = fold_inside(mid, &face.surface, &trim);
            let shifted = c
                .pcurve
                .transformed(&ogeom_math::Transform2::translation(folded - mid), tol)?;
            ogeom_algo::attach_pcurve(
                model,
                &built,
                shifted,
                surface_id,
                Location::identity(),
                sub_p,
            )?;
            Ok(built)
        }
        Tag::Section { section, range } => {
            let s = &fused.sections[*section];
            let domain = s.curve.domain();
            let (f0, f1) = folded_range(*range, domain, s.closed);
            let from = s.curve.point_at(at_param(f0, domain, s.closed), tol)?;
            let to = s.curve.point_at(at_param(f1, domain, s.closed), tol)?;
            let v0 = rebuild.vertex(from, tol);
            let v1 = rebuild.vertex(to, tol);
            let model = &mut *rebuild.model;
            let built = make_edge_between(model, s.curve.clone(), (f0, f1), &v0, &v1, tol)?.shape;
            // A piece welded into a junction some way off its own end (a
            // section tangent to an edge at the edge's vertex, stopped
            // where it began to hug the edge) owns that gap as the vertex
            // does: the face it bounds meets its neighbour there only that
            // closely, in its chart as in space.
            let at = |v: &Shape| {
                model
                    .node(v)
                    .and_then(|n| n.data().as_vertex())
                    .map(|d| d.point)
            };
            let ends = [at(&v0), at(&v1)];
            let welded = [(ends[0], from), (ends[1], to)]
                .iter()
                .filter_map(|(v, p)| v.map(|v| v.distance(*p)))
                .fold(0.0, f64::max);
            let own = s.tolerance.max(welded);
            // A fitted section's pieces and their ends own the section's
            // stated slop, exactly as a tolerant boundary's do: the ring
            // they close alternates between the two, and both sides must
            // meet at the junction's own resolution.
            if own > tol.confusion() {
                if let Some(node) = model.node_mut(&built)
                    && let ogeom_topo::NodeData::Edge(data) = node.data_mut()
                {
                    data.tolerance = data.tolerance.widen_to(own);
                }
                for v in [&v0, &v1] {
                    rebuild.widen(v, own)?;
                }
            }
            let model = &mut *rebuild.model;
            // The section's pcurve is unwrapped across any seam; the face's
            // triangulator lives in one chart, so the attached copy is folded
            // home by the same period shift the arrangement gave this
            // strand's polyline (to the side of the face's seam its trim is
            // on), decided by the sub-range's midpoint, so an
            // endpoint sitting exactly on the chart's edge stays on the side
            // the arc's body is.
            let pcurve = if from_a { &s.pc_a } else { &s.pc_b };
            let mid = pcurve.point_at(f64::midpoint(f0, f1), tol)?;
            let trim = face_trim_lines(face, tol);
            let trim: Vec<&[Point2]> = trim.iter().map(Vec::as_slice).collect();
            let folded = fold_inside(mid, &face.surface, &trim);
            let shifted =
                pcurve.transformed(&ogeom_math::Transform2::translation(folded - mid), tol)?;
            // The chart image ends where the vertex is, as the edge does in
            // space to its tolerance: a welded end left where the section
            // stopped leaves the face's outline open by the weld in its
            // chart, wider there than any later arrangement's snap.
            let loose = s.tolerance.max(tol.confusion() * 1e2);
            let bent = [
                rebuild.onto_vertex.contains(&v0.node()),
                rebuild.onto_vertex.contains(&v1.node()),
            ];
            let targets = [
                ends[0].filter(|v| bent[0] && v.distance(from) > loose),
                ends[1].filter(|v| bent[1] && v.distance(to) > loose),
            ];
            let shifted = if targets.iter().any(Option::is_some) {
                pcurve_onto_ends(&shifted, (f0, f1), &face.surface, targets, tol)?
            } else {
                shifted
            };
            ogeom_algo::attach_pcurve(
                model,
                &built,
                shifted,
                surface_id,
                Location::identity(),
                (f0, f1),
            )?;
            Ok(built)
        }
    }
}

/// A pcurve over `range` bent at its ends onto the chart points of
/// `targets`, where given: each end's correction fades linearly to nothing
/// at the other end, and the result is refitted at the same parameters, so
/// the image stays same-parameter with its edge to the correction's size.
fn pcurve_onto_ends(
    pcurve: &PlanarCurve,
    range: (f64, f64),
    surface: &SurfaceGeometry,
    targets: [Option<Point>; 2],
    tol: Tolerances,
) -> OgeomResult<PlanarCurve> {
    const SAMPLES: u32 = 32;
    let ts: Vec<f64> = (0..=SAMPLES)
        .map(|i| range.0 + (range.1 - range.0) * f64::from(i) / f64::from(SAMPLES))
        .collect();
    let mut points = ts
        .iter()
        .map(|t| pcurve.point_at(*t, tol))
        .collect::<OgeomResult<Vec<Point2>>>()?;
    let last = points.len() - 1;
    for (k, target) in targets.iter().enumerate() {
        let Some(target) = target else {
            continue;
        };
        let end = points[if k == 0 { 0 } else { last }];
        let foot = ogeom_algo::project_on_surface_from(surface, *target, (end.x, end.y), tol)?;
        let raw = Point2::new(foot.parameters.0, foot.parameters.1);
        // The foot may land a turn away from the end on a periodic chart,
        // whether or not the surface's wrapper says it is periodic; a
        // shifted candidate counts where the surface agrees it is the foot.
        let turn = core::f64::consts::TAU;
        let mut shifts = period_shifts(surface);
        for du in [-turn, 0.0, turn] {
            for dv in [-turn, 0.0, turn] {
                shifts.push((du, dv));
            }
        }
        let onto = shifts
            .into_iter()
            .map(|(du, dv)| Point2::new(raw.x + du, raw.y + dv))
            .filter(|c| {
                surface
                    .point_at(c.x, c.y, tol)
                    .is_ok_and(|q| q.distance(foot.point) <= tol.confusion() * 1e2)
            })
            .min_by(|a, b| {
                a.distance(end)
                    .partial_cmp(&b.distance(end))
                    .unwrap_or(core::cmp::Ordering::Equal)
            })
            .unwrap_or(raw);
        let delta = onto - end;
        for (i, p) in points.iter_mut().enumerate() {
            #[allow(clippy::cast_precision_loss)]
            let along = i as f64 / last as f64;
            let weight = if k == 0 { 1.0 - along } else { along };
            *p += delta * weight;
        }
    }
    let fitted = ogeom_geom::fit::fit_points_2d_at(&ts, &points, 3, tol.confusion() * 10.0, tol)?;
    Ok(PlanarCurve::BSpline(fitted.curve))
}

/// Sew kept pieces, demand closure, and nest shells into solids and voids.
fn assemble_result(
    model: &mut Model,
    fused: &GeneralFused,
    kept: &[(usize, bool)],
    a: &Shape,
    b: &Shape,
    tol: Tolerances,
) -> OgeomResult<Built> {
    ogeom_core::progress::stage("boolean: assemble");
    let mut history = History::new();
    let source_face = |piece: &FacePiece| -> Shape {
        if piece.from_a {
            fused.a.faces[piece.face].face.clone()
        } else {
            fused.b.faces[piece.face].face.clone()
        }
    };

    if kept.is_empty() {
        // A legitimate answer: cutting a solid away entirely leaves nothing.
        let empty = model.add_compound(&[])?;
        for piece in &fused.pieces {
            history.delete(&source_face(piece));
        }
        history.modify(a, empty.clone());
        history.modify(b, empty.clone());
        return Ok(Built::new(empty, history));
    }

    let weld = fused
        .contacts
        .iter()
        .map(|c| c.tolerance)
        .chain(fused.sections.iter().map(|s| s.tolerance))
        .fold(0.0_f64, |acc, t| acc.max(honest(t, tol) * 3.0));
    let floor = weld.max(tol.confusion() * 1e2);
    let junction_reach = fused
        .junctions
        .iter()
        .fold(0.0_f64, |acc, j| acc.max(j.reach));
    let mut junction_bins = bins::Bins::new(junction_reach + floor);
    for (index, j) in fused.junctions.iter().enumerate() {
        junction_bins.insert(j.at, index);
    }
    // The vertices' cells are as wide as the loosest tolerance a strand
    // brings to its ends, so an end is compared across a few cells, not
    // every vertex.
    let loosest = fused
        .a
        .faces
        .iter()
        .chain(fused.b.faces.iter())
        .flat_map(|f| f.edges.iter().map(|e| e.tolerance))
        .chain(fused.contacts.iter().map(|c| c.tolerance))
        .chain(fused.sections.iter().map(|s| s.tolerance))
        .fold(junction_reach, f64::max);
    let mut rebuild = Rebuild {
        model,
        surfaces_a: vec![None; fused.a.faces.len()],
        surfaces_b: vec![None; fused.b.faces.len()],
        vertices: Vec::new(),
        vertex_bins: bins::Bins::new(floor + loosest),
        widest: 0.0,
        weld,
        junctions: fused.junctions.iter().map(|j| (*j, None)).collect(),
        junction_bins,
        junction_reach,
        onto_vertex: std::collections::HashSet::new(),
    };
    let mut faces = Vec::new();
    let mut kept_sources: std::collections::HashSet<Shape> = std::collections::HashSet::new();
    for &(index, flip) in kept {
        let piece = &fused.pieces[index];
        let mut built = build_piece(&mut rebuild, fused, piece, tol)?;
        if flip {
            built = built.reversed();
        }
        history.modify(&source_face(piece), built.clone());
        kept_sources.insert(source_face(piece));
        faces.push(built);
    }
    for piece in &fused.pieces {
        let source = source_face(piece);
        if !kept_sources.contains(&source) {
            history.delete(&source);
        }
    }

    let model = rebuild.model;
    let sewn = sew(model, &faces, tol)?;
    for shell in &sewn.shells {
        if !is_shell_closed(model, shell)? {
            // Env-gated forensics: the open shell's unshared edges, the
            // question every failure here starts from.
            if *ARRANGE_DEBUG {
                use ogeom_geom::Curve3d as _;
                for edge in ogeom_topo::explore_unique(model, shell, ShapeType::Edge)? {
                    let users = ogeom_topo::explore(model, shell, Filter::OfType(ShapeType::Face))?
                        .iter()
                        .filter(|f| {
                            ogeom_topo::explore_unique(model, f, ShapeType::Edge)
                                .map(|es| es.iter().any(|e2| e2.node() == edge.node()))
                                .unwrap_or(false)
                        })
                        .count();
                    let mut occurrences = 0_usize;
                    for f in ogeom_topo::explore(model, shell, Filter::OfType(ShapeType::Face))? {
                        for wire in model.children_of(&f)? {
                            for e2 in model.children_of(&wire)? {
                                if e2.node() == edge.node() {
                                    occurrences += 1;
                                }
                            }
                        }
                    }
                    if occurrences % 2 == 1
                        && let Some(data) = model.node(&edge).and_then(|n| n.data().as_edge())
                        && let Some(ogeom_topo::EdgeRepr::Curve3d { curve, range, .. }) =
                            data.curve3d()
                        && let Some(g) = model.geometry().curve(*curve)
                    {
                        let a = g.point_at(range.0, tol)?;
                        let b = g.point_at(range.1, tol)?;
                        eprintln!(
                            "  open edge uses={occurrences} faces={users} {:?} range={range:?}: {a:?} -> {b:?}",
                            core::mem::discriminant(g)
                        );
                        // Twins: other edges of the shell with the same ends,
                        // which the sew should have found to be this edge.
                        for other in ogeom_topo::explore_unique(model, shell, ShapeType::Edge)? {
                            if other.node() == edge.node() {
                                continue;
                            }
                            let Some(od) = model.node(&other).and_then(|n| n.data().as_edge())
                            else {
                                continue;
                            };
                            let Some(ogeom_topo::EdgeRepr::Curve3d {
                                curve: oc,
                                range: orange,
                                ..
                            }) = od.curve3d()
                            else {
                                continue;
                            };
                            let Some(og) = model.geometry().curve(*oc) else {
                                continue;
                            };
                            let oa = og.point_at(orange.0, tol)?;
                            let ob = og.point_at(orange.1, tol)?;
                            let same = (oa.distance(a) <= 1e-3 && ob.distance(b) <= 1e-3)
                                || (oa.distance(b) <= 1e-3 && ob.distance(a) <= 1e-3);
                            if same {
                                let om = og.point_at(f64::midpoint(orange.0, orange.1), tol)?;
                                let m = g.point_at(f64::midpoint(range.0, range.1), tol)?;
                                eprintln!(
                                    "    twin {:?} range={orange:?} mid gap {:.2e} tol {:.2e} vs {:.2e}",
                                    core::mem::discriminant(og),
                                    om.distance(m),
                                    od.tolerance.get(),
                                    data.tolerance.get()
                                );
                            }
                        }
                        for f in ogeom_topo::explore(model, shell, Filter::OfType(ShapeType::Face))?
                        {
                            let uses_it = ogeom_topo::explore_unique(model, &f, ShapeType::Edge)
                                .map(|es| es.iter().any(|e2| e2.node() == edge.node()))
                                .unwrap_or(false);
                            if !uses_it {
                                continue;
                            }
                            if let Some(ogeom_topo::NodeData::Face(fd)) =
                                model.node(&f).map(|n| n.data())
                                && let Some(sg) = model.geometry().surface(fd.surface)
                            {
                                let bound = shape_bounds(model, &f, tol)?;
                                eprintln!(
                                    "    used by face on {:?} bound {:?}",
                                    core::mem::discriminant(sg),
                                    bound
                                );
                            }
                        }
                    }
                }
            }
            if *ARRANGE_DEBUG {
                use ogeom_geom::Curve3d as _;
                for f in ogeom_topo::explore(model, shell, Filter::OfType(ShapeType::Face))? {
                    let Some(ogeom_topo::NodeData::Face(fd)) = model.node(&f).map(|n| n.data())
                    else {
                        continue;
                    };
                    let Some(sg) = model.geometry().surface(fd.surface) else {
                        continue;
                    };
                    let edges = ogeom_topo::explore_unique(model, &f, ShapeType::Edge)?;
                    eprintln!(
                        "  face {:?} with {} edges",
                        core::mem::discriminant(sg),
                        edges.len()
                    );
                    if matches!(sg, SurfaceGeometry::Sphere(_)) {
                        for e in &edges {
                            if let Some(ed) = model.node(e).and_then(|n| n.data().as_edge())
                                && let Some(ogeom_topo::EdgeRepr::Curve3d { curve, range, .. }) =
                                    ed.curve3d()
                                && let Some(g) = model.geometry().curve(*curve)
                            {
                                let a = g.point_at(range.0, tol)?;
                                let b = g.point_at(range.1, tol)?;
                                eprintln!(
                                    "    sphere edge {:?} ({:.4},{:.4},{:.4}) -> ({:.4},{:.4},{:.4})",
                                    core::mem::discriminant(g),
                                    a.x,
                                    a.y,
                                    a.z,
                                    b.x,
                                    b.y,
                                    b.z
                                );
                            }
                        }
                    }
                }
            }
            ogeom_bail!(
                NotDone,
                "the kept pieces did not close into a shell; the configuration \
                 is beyond what the boolean currently resolves"
            );
        }
    }

    // Nest: a shell whose bound sits inside another's is that solid's void.
    let mut bounds = Vec::new();
    for shell in &sewn.shells {
        bounds.push(shape_bounds(model, shell, tol)?);
    }
    let mut solids = Vec::new();
    for (i, shell) in sewn.shells.iter().enumerate() {
        let contained = bounds
            .iter()
            .enumerate()
            .any(|(j, other)| j != i && other.contains_box(&bounds[i]));
        if contained {
            continue;
        }
        let mut group = vec![shell.clone()];
        for (j, candidate) in sewn.shells.iter().enumerate() {
            if j != i && bounds[i].contains_box(&bounds[j]) {
                group.push(candidate.clone());
            }
        }
        solids.push(model.add_solid(&group)?);
    }
    let result = if solids.len() == 1 {
        solids.remove(0)
    } else {
        model.add_compound(&solids)?
    };
    history.modify(a, result.clone());
    history.modify(b, result.clone());
    Ok(Built::new(result, history))
}

/// Solids from an unordered soup of faces: sew, demand closure, nest
/// shells into solids and voids: the pipeline's own final stages offered
/// as a builder.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if a
/// shell fails to close; an open soup encloses no volume, and saying so
/// beats guessing.
pub fn make_volume(model: &mut Model, faces: &[Shape], tol: Tolerances) -> OgeomResult<Built> {
    let sewn = ogeom_algo::sew(model, faces, tol)?;
    for shell in &sewn.shells {
        if !ogeom_algo::is_shell_closed(model, shell)? {
            ogeom_bail!(
                Construction,
                "the faces do not close into shells; an open soup encloses no volume"
            );
        }
    }
    let mut bounds = Vec::new();
    for shell in &sewn.shells {
        bounds.push(ogeom_algo::shape_bounds(model, shell, tol)?);
    }
    let mut solids = Vec::new();
    for (i, shell) in sewn.shells.iter().enumerate() {
        let contained = bounds
            .iter()
            .enumerate()
            .any(|(j, other)| j != i && other.contains_box(&bounds[i]));
        if contained {
            continue;
        }
        let mut group = vec![shell.clone()];
        for (j, candidate) in sewn.shells.iter().enumerate() {
            if j != i && bounds[i].contains_box(&bounds[j]) {
                // A void bounds its solid from inside: material lies outside
                // it, so the sewn outward orientation reverses.
                group.push(candidate.reversed());
            }
        }
        solids.push(model.add_solid(&group)?);
    }
    let mut history = History::new();
    let result = if solids.len() == 1 {
        solids.remove(0)
    } else {
        model.add_compound(&solids)?
    };
    for face in faces {
        history.modify(face, result.clone());
    }
    Ok(Built::new(result, history))
}

/// The three cells two solids cut space into, each one boolean's answer:
/// what is only in `a`, what is only in `b`, and what is in both. Arbitrary
/// set expressions compose by fusing a selection of these.
#[derive(Debug)]
pub struct Cells {
    /// `a` with `b` removed.
    pub a_not_b: Built,
    /// `b` with `a` removed.
    pub b_not_a: Built,
    /// The overlap.
    pub common: Built,
}

/// Split two solids into their three cells.
///
/// # Errors
///
/// As the operations themselves.
pub fn cells(model: &mut Model, a: &Shape, b: &Shape, tol: Tolerances) -> OgeomResult<Cells> {
    Ok(Cells {
        a_not_b: cut(model, a, b, tol)?,
        b_not_a: cut(model, b, a, tol)?,
        common: common(model, a, b, tol)?,
    })
}

/// A tolerance whose confusion *is* the stated fuzz: every gap, pave and
/// weld decision inherits it coherently, which is what a fuzzy boolean
/// means.
fn fuzzed(fuzz: f64, tol: Tolerances) -> OgeomResult<Tolerances> {
    if !fuzz.is_finite() || fuzz <= 0.0 {
        ogeom_bail!(Construction, "a fuzz of {fuzz} is not a distance");
    }
    if fuzz <= tol.confusion() {
        return Ok(tol);
    }
    Tolerances::with_scale(ogeom_core::tolerance::CONFUSION / fuzz)
}

/// [`fuse`] with geometry within `fuzz` of touching counted as touching.
///
/// # Errors
///
/// As [`fuse`], plus a non-positive fuzz.
pub fn fuse_fuzzy(
    model: &mut Model,
    a: &Shape,
    b: &Shape,
    fuzz: f64,
    tol: Tolerances,
) -> OgeomResult<Built> {
    let loosened = fuzzed(fuzz, tol)?;
    fuse(model, a, b, loosened)
}

/// [`cut`] at a stated fuzz.
///
/// # Errors
///
/// As [`fuse_fuzzy`].
pub fn cut_fuzzy(
    model: &mut Model,
    a: &Shape,
    b: &Shape,
    fuzz: f64,
    tol: Tolerances,
) -> OgeomResult<Built> {
    let loosened = fuzzed(fuzz, tol)?;
    cut(model, a, b, loosened)
}

/// A shape repeated `count` times along a direction at a period and fused
/// into one: the periodic pattern as a composition of what exists.
///
/// # Errors
///
/// As [`fuse`], plus an unusable count or period.
pub fn make_periodic(
    model: &mut Model,
    shape: &Shape,
    step: ogeom_math::Vector,
    count: usize,
    tol: Tolerances,
) -> OgeomResult<Built> {
    if count == 0 {
        ogeom_bail!(Construction, "a pattern of zero copies is nothing");
    }
    if step.magnitude() <= tol.confusion() {
        ogeom_bail!(Construction, "a zero period stacks every copy on the first");
    }
    let mut history = History::new();
    let mut result = shape.clone();
    for i in 1..count {
        #[allow(clippy::cast_precision_loss, reason = "pattern counts are small")]
        let offset = step * i as f64;
        let moved =
            ogeom_algo::transformed(model, shape, ogeom_math::Transform::translation(offset))?
                .shape;
        let joined = fuse(model, &result, &moved, tol)?;
        history.modify(&moved, joined.shape.clone());
        result = joined.shape;
    }
    history.generate(shape, result.clone());
    Ok(Built::new(result, history))
}

// --- the operations ----------------------------------------------------------

/// A shape whose placements carry scale, rebuilt with the scale baked
/// into its geometry before the pipeline runs.
///
/// A scale changes a surface's parameterization out from under its
/// pcurves (the old refusal), so the bake is the whole-shape conversion:
/// surfaces restated in world space, edges moved exactly, pcurves
/// re-derived against the new parameterizations. Unscaled shapes pass
/// through untouched.
fn baked_if_scaled(model: &mut Model, shape: &Shape, tol: Tolerances) -> OgeomResult<Shape> {
    let mut restate = false;
    for face in ogeom_topo::explore(model, shape, Filter::OfType(ShapeType::Face))? {
        let placement = face.transform(model.datums())?;
        // A scale changes lengths the melt compares; a reflection flips
        // every chart's natural normal against its face's flag. Either way
        // the operand is restated in world coordinates first, where both
        // effects are already folded in.
        if (placement.scale_factor().abs() - 1.0).abs() > 1e-9 || !placement.preserves_handedness()
        {
            restate = true;
            break;
        }
    }
    if !restate {
        return Ok(shape.clone());
    }
    Ok(ogeom_algo::baked_shape(model, shape, tol)?.shape)
}

/// Whether a shape is a half space: one shell of one face, open by
/// construction or closed inside out.
fn is_half_space(model: &Model, shape: &Shape) -> OgeomResult<bool> {
    Ok(half_space::half_space_face(model, shape, Tolerances::millimetres())?.is_some())
}

/// A half space resolved into the solid the operation can act on.
///
/// A planar boundary becomes a box filling the material side, sized past
/// the other argument's whole reach. The box's plane-side face is
/// *coplanar with the boundary itself*, so the cut the caller sees is the
/// exact plane; its far faces stand outside everything the other shape
/// reaches and never appear in the result. A curved boundary is resolved
/// by [`half_space::resolved`]. A shape that is not a half space passes
/// through untouched.
fn resolved_half_space(
    model: &mut Model,
    shape: &Shape,
    other: &Shape,
    tol: Tolerances,
) -> OgeomResult<Shape> {
    let Some(face) = half_space::half_space_face(model, shape, tol)? else {
        return Ok(shape.clone());
    };
    let planar = model
        .node(&face)
        .and_then(|n| n.data().as_face())
        .and_then(|d| model.geometry().surface(d.surface))
        .is_some_and(|s| matches!(s, SurfaceGeometry::Plane(_)));
    if !planar {
        return half_space::resolved(model, &face, other, tol);
    }
    // The boundary's outward normal points away from the material.
    let (at, outward) = ogeom_algo::face_normal(model, &face, tol)?;
    let bound = ogeom_algo::shape_bounds(model, other, tol)?;
    let Some(centre) = bound.centre() else {
        ogeom_bail!(
            Construction,
            "the other argument has no bound to fill against"
        );
    };
    let reach = bound.diagonal().max(tol.confusion() * 1e3) * 2.0;
    let into = ogeom_math::Direction::new(-outward, tol)?;
    let foot = centre - outward * outward.dot(centre - at);
    let seed = if into.vector().x.abs() < 0.9 {
        ogeom_math::Vector::new(1.0, 0.0, 0.0)
    } else {
        ogeom_math::Vector::new(0.0, 1.0, 0.0)
    };
    let frame_x = ogeom_math::Direction::from_cross(into.vector(), seed, tol)?;
    let oriented = ogeom_math::Frame::new(foot, into, frame_x, tol)?;
    let corner =
        foot - oriented.x().vector() * (reach / 2.0) - oriented.y().vector() * (reach / 2.0);
    let placed = ogeom_math::Frame::new(corner, into, frame_x, tol)?;
    Ok(ogeom_algo::make_box(model, placed, (reach, reach, reach), tol)?.shape)
}

/// The union of two solids.
///
/// # Errors
///
/// [`OgeomError::NotDone`](ogeom_core::OgeomError::NotDone) for configurations the
/// boolean refuses: tangential or same-domain contact, scaled placements;
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) for arguments
/// that are not closed solids.
pub fn fuse(model: &mut Model, a: &Shape, b: &Shape, tol: Tolerances) -> OgeomResult<Built> {
    if is_half_space(model, a)? || is_half_space(model, b)? {
        ogeom_bail!(
            Construction,
            "the union with a half space is unbounded; a half space serves cut, \
             common and section"
        );
    }
    let (a, b) = (
        &baked_if_scaled(model, a, tol)?,
        &baked_if_scaled(model, b, tol)?,
    );
    let fused = general_fuse(model, a, b, tol)?;
    // Outward pieces bound the union; a same-domain pair with aligned
    // material keeps one copy, and one with opposed material is interior to
    // the union and vanishes.
    let kept: Vec<(usize, bool)> = fused
        .pieces
        .iter()
        .enumerate()
        .filter(|(_, p)| {
            p.state == PieceState::Out || (p.state == PieceState::OnAligned && !p.covered)
        })
        .map(|(i, _)| (i, false))
        .collect();
    assemble_result(model, &fused, &kept, a, b, tol)
}

/// The intersection of two solids.
///
/// # Errors
///
/// As [`fuse`].
pub fn common(model: &mut Model, a: &Shape, b: &Shape, tol: Tolerances) -> OgeomResult<Built> {
    let (a, b) = (
        &resolved_half_space(model, a, b, tol)?,
        &resolved_half_space(model, b, a, tol)?,
    );
    let (a, b) = (
        &baked_if_scaled(model, a, tol)?,
        &baked_if_scaled(model, b, tol)?,
    );
    let fused = general_fuse(model, a, b, tol)?;
    // Inward pieces bound the intersection; an aligned same-domain pair
    // bounds it too, once. An opposed pair encloses no volume between them.
    let kept: Vec<(usize, bool)> = fused
        .pieces
        .iter()
        .enumerate()
        .filter(|(_, p)| {
            p.state == PieceState::In || (p.state == PieceState::OnAligned && !p.covered)
        })
        .map(|(i, _)| (i, false))
        .collect();
    assemble_result(model, &fused, &kept, a, b, tol)
}

/// The first solid with the second removed.
///
/// The pieces of `b` that close the cut into `a`'s material bound the removed
/// volume from `b`'s side, so they join the result with their material side
/// flipped.
///
/// # Errors
///
/// As [`fuse`].
pub fn cut(model: &mut Model, a: &Shape, b: &Shape, tol: Tolerances) -> OgeomResult<Built> {
    let (a, b) = (
        &resolved_half_space(model, a, b, tol)?,
        &resolved_half_space(model, b, a, tol)?,
    );
    let (a, b) = (
        &baked_if_scaled(model, a, tol)?,
        &baked_if_scaled(model, b, tol)?,
    );
    let fused = general_fuse(model, a, b, tol)?;
    // The first argument's outward pieces stay; the tool's inward pieces
    // close the cut with their material side flipped. On the shared surface:
    // an opposed pair means the tool's material is entirely on the other
    // side, so the first argument's face survives untouched; an aligned pair
    // means the tool's material backs the same wall, which the cut removes.
    let kept: Vec<(usize, bool)> = fused
        .pieces
        .iter()
        .enumerate()
        .filter_map(|(i, p)| match (p.from_a, p.state) {
            (true, PieceState::Out) => Some((i, false)),
            (true, PieceState::OnOpposed) => Some((i, false)),
            (false, PieceState::In) => Some((i, true)),
            _ => None,
        })
        .collect();
    assemble_result(model, &fused, &kept, a, b, tol)
}

/// The edges where the two solids' boundaries cross.
///
/// # Errors
///
/// As [`fuse`].
pub fn section(model: &mut Model, a: &Shape, b: &Shape, tol: Tolerances) -> OgeomResult<Built> {
    let (a, b) = (
        &resolved_half_space(model, a, b, tol)?,
        &resolved_half_space(model, b, a, tol)?,
    );
    let (a, b) = (
        &baked_if_scaled(model, a, tol)?,
        &baked_if_scaled(model, b, tol)?,
    );
    let fused = general_fuse(model, a, b, tol)?;
    // Every section sub-edge that survived into some piece's ring, built
    // once per distinct sub-range.
    let mut wanted: Vec<(usize, (f64, f64))> = Vec::new();
    for piece in &fused.pieces {
        for ring in &piece.rings {
            for traversal in ring {
                if let Tag::Section { section, range } = &traversal.tag {
                    let near = |a: (f64, f64), b: (f64, f64)| {
                        (a.0 - b.0).abs() <= tol.parametric()
                            && (a.1 - b.1).abs() <= tol.parametric()
                    };
                    if !wanted.iter().any(|(s, r)| s == section && near(*r, *range)) {
                        wanted.push((*section, *range));
                    }
                }
            }
        }
    }
    let mut history = History::new();
    let mut edges = Vec::new();
    for (si, range) in wanted {
        let s = &fused.sections[si];
        let domain = s.curve.domain();
        let (f0, f1) = folded_range(range, domain, s.closed);
        let from = s.curve.point_at(at_param(f0, domain, s.closed), tol)?;
        let to = s.curve.point_at(at_param(f1, domain, s.closed), tol)?;
        let v0 = make_vertex(model, from).shape;
        let v1 = make_vertex(model, to).shape;
        edges.push(make_edge_between(model, s.curve.clone(), (f0, f1), &v0, &v1, tol)?.shape);
    }
    // Contacts are not crossings, so no piece's ring carries them and the
    // loop above cannot see them, but a section through a tangency has a
    // curve in it, and this is where it comes from.
    for contact in &fused.tangents {
        for (lo, hi) in contact_intervals(&fused, contact, tol)? {
            let from = contact.curve.point_at(lo, tol)?;
            let to = contact.curve.point_at(hi, tol)?;
            let v0 = make_vertex(model, from).shape;
            let v1 = if from.distance(to) <= tol.confusion() {
                v0.clone()
            } else {
                make_vertex(model, to).shape
            };
            edges.push(
                make_edge_between(model, contact.curve.clone(), (lo, hi), &v0, &v1, tol)?.shape,
            );
        }
    }
    let result = model.add_compound(&edges)?;
    history.modify(a, result.clone());
    history.modify(b, result.clone());
    Ok(Built::new(result, history))
}

#[cfg(test)]
#[allow(clippy::unwrap_used)]
mod tests {
    use super::*;
    use ogeom_algo::{check, make_box, make_cylinder, volume_properties};
    use ogeom_math::{Direction, Frame};
    use ogeom_mesh::Deflection;

    const T: Tolerances = Tolerances::millimetres();
    const PI: f64 = core::f64::consts::PI;

    #[test]
    fn coincidence_is_measured_over_the_overlap_and_nowhere_else() {
        // Patches restated from planes: the geometry no longer says "plane",
        // which is the whole reason this measurement exists.
        let patch = |plane: ogeom_math::Plane, u: (f64, f64), v: (f64, f64)| {
            let surface: SurfaceGeometry =
                ogeom_geom::PlaneSurface::over(plane, u, v).unwrap().into();
            SurfaceGeometry::from(surface.to_bspline(T).unwrap())
        };
        let reach = T.confusion() * 1e2;

        // Two windows on one plane, overlapping over a quarter of each. They
        // are the same surface exactly where they meet, which is the claim.
        let here = patch(ogeom_math::Plane::XY, (0.0, 10.0), (0.0, 10.0));
        let over = patch(ogeom_math::Plane::XY, (5.0, 15.0), (5.0, 15.0));
        assert!(surfaces_coincide(&here, &over, reach, T));

        // The same plane lifted clear of itself is not the same surface, and
        // a plane square to it crosses rather than coincides: the case that
        // must keep marching, since a crossing has a section to find.
        let above = patch(
            ogeom_math::Plane::new(
                Frame::new(Point::new(0.0, 0.0, 1.0), Direction::Z, Direction::X, T).unwrap(),
            ),
            (0.0, 10.0),
            (0.0, 10.0),
        );
        assert!(!surfaces_coincide(&here, &above, reach, T));
        let across = patch(
            ogeom_math::Plane::new(
                Frame::new(Point::new(5.0, 0.0, 0.0), Direction::X, Direction::Y, T).unwrap(),
            ),
            (0.0, 10.0),
            (0.0, 10.0),
        );
        assert!(!surfaces_coincide(&here, &across, reach, T));
    }

    fn frame_at(origin: Point) -> Frame {
        Frame::new(origin, Direction::Z, Direction::X, T).unwrap()
    }

    fn boxes(model: &mut Model) -> (Shape, Shape) {
        // Overlapping in the corner cube [1,2]^3.
        let a = make_box(model, Frame::WORLD, (2.0, 2.0, 2.0), T).unwrap();
        let b = make_box(
            model,
            frame_at(Point::new(1.0, 1.0, 1.0)),
            (2.0, 2.0, 2.0),
            T,
        )
        .unwrap();
        (a.shape, b.shape)
    }

    fn volume(model: &Model, shape: &Shape) -> f64 {
        let fine = Deflection {
            chord: 1e-3,
            ..Deflection::default()
        };
        volume_properties(model, shape, fine, T).unwrap().mass
    }

    fn assert_valid(model: &Model, shape: &Shape) {
        let diagnosis = check(model, shape, T).unwrap();
        assert!(
            diagnosis.is_valid(),
            "the result fails validity: {:?}",
            diagnosis.problems
        );
    }

    #[test]
    fn fuse_of_overlapping_boxes_has_the_inclusion_exclusion_volume() {
        let mut model = Model::new();
        let (a, b) = boxes(&mut model);
        let fused = fuse(&mut model, &a, &b, T).unwrap();
        assert_valid(&model, &fused.shape);
        assert!((volume(&model, &fused.shape) - 15.0).abs() < 1e-9);
        assert_eq!(
            fused.history.modified(&a),
            std::slice::from_ref(&fused.shape)
        );
    }

    #[test]
    fn common_of_overlapping_boxes_is_the_overlap_cube() {
        let mut model = Model::new();
        let (a, b) = boxes(&mut model);
        let result = common(&mut model, &a, &b, T).unwrap();
        assert_valid(&model, &result.shape);
        assert!((volume(&model, &result.shape) - 1.0).abs() < 1e-9);
    }

    #[test]
    fn cut_removes_the_overlap_from_the_first_argument() {
        let mut model = Model::new();
        let (a, b) = boxes(&mut model);
        let result = cut(&mut model, &a, &b, T).unwrap();
        assert_valid(&model, &result.shape);
        assert!((volume(&model, &result.shape) - 7.0).abs() < 1e-9);
    }

    #[test]
    fn section_of_overlapping_boxes_is_the_six_segment_seam() {
        let mut model = Model::new();
        let (a, b) = boxes(&mut model);
        let result = section(&mut model, &a, &b, T).unwrap();
        let edges = explore(&model, &result.shape, Filter::OfType(ShapeType::Edge)).unwrap();
        assert_eq!(edges.len(), 6);
    }

    #[test]
    fn fuse_of_disjoint_boxes_is_a_compound_of_both() {
        let mut model = Model::new();
        let a = make_box(&mut model, Frame::WORLD, (2.0, 2.0, 2.0), T).unwrap();
        let b = make_box(
            &mut model,
            frame_at(Point::new(5.0, 0.0, 0.0)),
            (2.0, 2.0, 2.0),
            T,
        )
        .unwrap();
        let fused = fuse(&mut model, &a.shape, &b.shape, T).unwrap();
        assert_eq!(model.kind_of(&fused.shape).unwrap(), ShapeType::Compound);
        assert!((volume(&model, &fused.shape) - 16.0).abs() < 1e-9);
    }

    #[test]
    fn cutting_a_through_post_leaves_a_slab_with_a_hole() {
        let mut model = Model::new();
        let slab = make_box(&mut model, Frame::WORLD, (4.0, 4.0, 1.0), T).unwrap();
        let post = make_box(
            &mut model,
            frame_at(Point::new(1.5, 1.5, -1.0)),
            (1.0, 1.0, 3.0),
            T,
        )
        .unwrap();
        let result = cut(&mut model, &slab.shape, &post.shape, T).unwrap();
        assert_valid(&model, &result.shape);
        assert!((volume(&model, &result.shape) - 15.0).abs() < 1e-9);
    }

    #[test]
    fn common_with_a_contained_box_is_that_box() {
        let mut model = Model::new();
        let outer = make_box(&mut model, Frame::WORLD, (6.0, 6.0, 6.0), T).unwrap();
        let inner = make_box(
            &mut model,
            frame_at(Point::new(2.0, 2.0, 2.0)),
            (2.0, 2.0, 2.0),
            T,
        )
        .unwrap();
        let result = common(&mut model, &outer.shape, &inner.shape, T).unwrap();
        assert_valid(&model, &result.shape);
        assert!((volume(&model, &result.shape) - 8.0).abs() < 1e-9);
    }

    #[test]
    fn cutting_everything_away_leaves_an_empty_compound() {
        let mut model = Model::new();
        let outer = make_box(&mut model, Frame::WORLD, (6.0, 6.0, 6.0), T).unwrap();
        let inner = make_box(
            &mut model,
            frame_at(Point::new(2.0, 2.0, 2.0)),
            (2.0, 2.0, 2.0),
            T,
        )
        .unwrap();
        let result = cut(&mut model, &inner.shape, &outer.shape, T).unwrap();
        assert_eq!(model.kind_of(&result.shape).unwrap(), ShapeType::Compound);
        assert!(
            explore(&model, &result.shape, Filter::OfType(ShapeType::Face))
                .unwrap()
                .is_empty()
        );
    }

    #[test]
    fn drilling_a_box_leaves_a_cylindrical_hole() {
        // The curved milestone: box minus a through-post cylinder. The box's
        // top and bottom faces come back with *circular* holes, the hole's
        // wall is the cylinder's own surface with its material side flipped,
        // and the cylinder's seam and both section circles all had to split
        // and sew for the shell to close.
        let mut model = Model::new();
        let block = make_box(&mut model, Frame::WORLD, (4.0, 4.0, 1.0), T).unwrap();
        let drill = make_cylinder(
            &mut model,
            frame_at(Point::new(2.0, 2.0, -1.0)),
            0.5,
            3.0,
            T,
        )
        .unwrap();
        let result = cut(&mut model, &block.shape, &drill.shape, T).unwrap();
        assert_valid(&model, &result.shape);
        let exact = 16.0 - PI * 0.25;
        let got = volume(&model, &result.shape);
        assert!(
            (got - exact).abs() / exact < 2e-3,
            "volume {got} against {exact}"
        );
    }

    #[test]
    fn common_of_a_box_and_a_cylinder_is_the_post_inside_it() {
        let mut model = Model::new();
        let block = make_box(&mut model, Frame::WORLD, (4.0, 4.0, 1.0), T).unwrap();
        let post = make_cylinder(
            &mut model,
            frame_at(Point::new(2.0, 2.0, -1.0)),
            0.5,
            3.0,
            T,
        )
        .unwrap();
        let result = common(&mut model, &block.shape, &post.shape, T).unwrap();
        assert_valid(&model, &result.shape);
        let exact = PI * 0.25;
        let got = volume(&model, &result.shape);
        assert!(
            (got - exact).abs() / exact < 2e-3,
            "volume {got} against {exact}"
        );
    }

    #[test]
    fn fusing_a_post_onto_a_slab_adds_what_stands_proud() {
        let mut model = Model::new();
        let slab = make_box(&mut model, Frame::WORLD, (4.0, 4.0, 1.0), T).unwrap();
        let post = make_cylinder(
            &mut model,
            frame_at(Point::new(2.0, 2.0, -1.0)),
            0.5,
            3.0,
            T,
        )
        .unwrap();
        let result = fuse(&mut model, &slab.shape, &post.shape, T).unwrap();
        assert_valid(&model, &result.shape);
        let exact = PI.mul_add(0.25 * 3.0, 16.0) - PI * 0.25;
        let got = volume(&model, &result.shape);
        assert!(
            (got - exact).abs() / exact < 2e-3,
            "volume {got} against {exact}"
        );
    }

    #[test]
    fn crossed_cylinders_run_through_the_marched_sections() {
        // No closed form exists for cylinder/cylinder: the sections are
        // marched and fitted, and every stage downstream (paves against
        // seams on both charts, splitting, classification, rebuilding with
        // fitted pcurves, sewing, meshing) has to work within the fit's
        // stated budget. Unequal radii keep the tangential branch points
        // away. The volumes have no easy closed form, so the operations are
        // held to each other: fuse = A + B - common and cut = A - common are
        // identities whatever the shapes.
        let mut model = Model::new();
        let upright = make_cylinder(&mut model, Frame::WORLD, 1.0, 4.0, T).unwrap();
        let across_frame =
            Frame::new(Point::new(-2.0, 0.0, 2.0), Direction::X, Direction::Y, T).unwrap();
        let across = make_cylinder(&mut model, across_frame, 0.6, 4.0, T).unwrap();

        let both = fuse(&mut model, &upright.shape, &across.shape, T).unwrap();
        assert_valid(&model, &both.shape);
        let shared = common(&mut model, &upright.shape, &across.shape, T).unwrap();
        assert_valid(&model, &shared.shape);
        let pierced = cut(&mut model, &upright.shape, &across.shape, T).unwrap();
        assert_valid(&model, &pierced.shape);

        let va = volume(&model, &upright.shape);
        let vb = volume(&model, &across.shape);
        let vf = volume(&model, &both.shape);
        let vc = volume(&model, &shared.shape);
        let vx = volume(&model, &pierced.shape);

        assert!(vc > 0.0 && vc < vb, "the overlap is real and partial: {vc}");
        assert!(
            (vf - (va + vb - vc)).abs() / vf < 2e-3,
            "fuse {vf} against A + B - common {}",
            va + vb - vc
        );
        assert!(
            (vx - (va - vc)).abs() / vx < 2e-3,
            "cut {vx} against A - common {}",
            va - vc
        );
    }

    #[test]
    fn stacked_boxes_fuse_into_one_solid_and_the_contact_vanishes() {
        // Same-domain contact with opposed materials: the shared rectangle is
        // interior to the union and no face of the result may carry it.
        let mut model = Model::new();
        let lower = make_box(&mut model, Frame::WORLD, (2.0, 2.0, 1.0), T).unwrap();
        let upper = make_box(
            &mut model,
            frame_at(Point::new(0.0, 0.0, 1.0)),
            (2.0, 2.0, 1.0),
            T,
        )
        .unwrap();
        let fused = fuse(&mut model, &lower.shape, &upper.shape, T).unwrap();
        assert_valid(&model, &fused.shape);
        assert_eq!(model.kind_of(&fused.shape).unwrap(), ShapeType::Solid);
        assert!((volume(&model, &fused.shape) - 8.0).abs() < 1e-9);
    }

    #[test]
    fn a_small_box_on_a_big_one_fuses_with_a_partial_contact() {
        // The big top face splits into the contact rectangle (swallowed)
        // and the surround, which stays and must sew to the small box's
        // walls along the contact's edges.
        let mut model = Model::new();
        let big = make_box(&mut model, Frame::WORLD, (4.0, 4.0, 1.0), T).unwrap();
        let small = make_box(
            &mut model,
            frame_at(Point::new(1.0, 1.0, 1.0)),
            (2.0, 2.0, 1.0),
            T,
        )
        .unwrap();
        let fused = fuse(&mut model, &big.shape, &small.shape, T).unwrap();
        assert_valid(&model, &fused.shape);
        assert!((volume(&model, &fused.shape) - 20.0).abs() < 1e-9);

        // Cutting the same pair removes nothing but the measure-zero contact:
        // the big box survives whole, its top face's contact piece intact.
        let mut model = Model::new();
        let big = make_box(&mut model, Frame::WORLD, (4.0, 4.0, 1.0), T).unwrap();
        let small = make_box(
            &mut model,
            frame_at(Point::new(1.0, 1.0, 1.0)),
            (2.0, 2.0, 1.0),
            T,
        )
        .unwrap();
        let result = cut(&mut model, &big.shape, &small.shape, T).unwrap();
        assert_valid(&model, &result.shape);
        assert!((volume(&model, &result.shape) - 16.0).abs() < 1e-9);
    }

    #[test]
    fn flush_walls_fuse_cut_and_meet_with_aligned_contact() {
        // Overlapping boxes sharing flush walls: same-domain contact with
        // *aligned* materials. A = [0,2]^3, B = [1,3]x[0,2]x[0,2]: the y and
        // z walls of the overlap are coplanar with aligned outward normals.
        let mut model = Model::new();
        let a = make_box(&mut model, Frame::WORLD, (2.0, 2.0, 2.0), T).unwrap();
        let b = make_box(
            &mut model,
            frame_at(Point::new(1.0, 0.0, 0.0)),
            (2.0, 2.0, 2.0),
            T,
        )
        .unwrap();
        let fused = fuse(&mut model, &a.shape, &b.shape, T).unwrap();
        assert_valid(&model, &fused.shape);
        assert!((volume(&model, &fused.shape) - 12.0).abs() < 1e-9);

        let mut model = Model::new();
        let a = make_box(&mut model, Frame::WORLD, (2.0, 2.0, 2.0), T).unwrap();
        let b = make_box(
            &mut model,
            frame_at(Point::new(1.0, 0.0, 0.0)),
            (2.0, 2.0, 2.0),
            T,
        )
        .unwrap();
        let shared = common(&mut model, &a.shape, &b.shape, T).unwrap();
        assert_valid(&model, &shared.shape);
        assert!((volume(&model, &shared.shape) - 4.0).abs() < 1e-9);

        let mut model = Model::new();
        let a = make_box(&mut model, Frame::WORLD, (2.0, 2.0, 2.0), T).unwrap();
        let b = make_box(
            &mut model,
            frame_at(Point::new(1.0, 0.0, 0.0)),
            (2.0, 2.0, 2.0),
            T,
        )
        .unwrap();
        let cut_result = cut(&mut model, &a.shape, &b.shape, T).unwrap();
        assert_valid(&model, &cut_result.shape);
        assert!((volume(&model, &cut_result.shape) - 4.0).abs() < 1e-9);
    }

    #[test]
    fn every_source_face_is_accounted_for_in_the_history() {
        let mut model = Model::new();
        let (a, b) = boxes(&mut model);
        let result = cut(&mut model, &a, &b, T).unwrap();
        for face in explore(&model, &a, Filter::OfType(ShapeType::Face)).unwrap() {
            assert!(
                result.history.is_affected(&face),
                "a face of the first argument vanished from the history"
            );
        }
    }
}