geop_lib_booleans/

boolean.rs

//! Boolean operations on two solids, built on top of [`remesh`].
//!
//! The algorithm is deliberately simple, and it is simple *because* remesh
//! has already done the hard part. After remeshing, every face of either
//! solid lies wholly inside the other, wholly outside it, or exactly on its
//! boundary — no face straddles, because every intersection curve has been
//! imprinted and every face it crossed has been split. That turns a boolean
//! into a per-face classification followed by a keep/drop table.
//!
//! 1. [`remesh`] the two solids against each other.
//! 2. Classify each face with [`classify_face`], by taking a point strictly
//!    inside its trimmed region and asking where it sits relative to the
//!    other solid.
//! 3. Keep the faces the operator wants (see `keeps`), reversing
//!    them where the operator needs the material on the other side.
//! 4. Assemble the survivors into a new solid and discard everything else.
//!
//! [`boolean_up_to_next`] is the same with one more step between 3 and 4:
//! of the second solid, a tool, it keeps only the first piece its start
//! reaches on the side of the first solid's boundary the operator wants.

use std::collections::{HashMap, HashSet};

use geop_core_geometry::{intersection::curve_surface_intersect, nurb_surface::NurbSurface3D};
use geop_core_math::{
    geop_error::{GeopError, GeopResult, WithContext},
    scalars::Scalar,
    union_find::UnionFind,
    vector::Vector3,
    with_context,
};
use geop_core_topology::{
    EdgeId, FaceId, Model, SolidId, VertexId,
    contains::{
        face::face_interior_point_where,
        shell::{Classified, PointClassification as ShellPoint, solid_contains, solid_face_at},
    },
};
use geop_ops_model::{BooleanOp, NamedModel, Namer};

use crate::naming::{name_lone_pieces, origins_by_name};
use crate::remesh::{Origins, RemeshParams, remesh};

/// Where one solid's face sits relative to the *other* solid.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub(crate) enum FaceClassification {
    /// Strictly inside the other solid.
    Inside,
    /// Strictly outside it.
    Outside,
    /// On its boundary, with the two surfaces' normals pointing the same way
    /// — the two solids touch and lie on the same side of the shared patch.
    OnSameNormal,
    /// On its boundary, with the normals opposed — the solids meet along the
    /// patch from opposite sides.
    OnOppositeNormal,
}

/// What to do with a classified face.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
enum Keep {
    /// Keep it as it is.
    AsIs,
    /// Keep it with its normal flipped — the operator wants the material on
    /// the other side of this patch.
    Reversed,
    /// Drop it.
    Drop,
}

/// Whether to keep a face classified as `class`, and with which
/// orientation. `from_a` says which solid the face came from, because the
/// two are not symmetric for [`BooleanOp::Difference`], and because a
/// coincident patch must be kept exactly once rather than from both.
///
/// The coincident rules are the only subtle ones. A patch shared by both
/// solids with **matching** normals bounds the same material on the same
/// side, so union and intersection each keep exactly one copy (`a`'s, by
/// convention) and difference deletes it — the material behind it is
/// removed from both sides at once. A patch shared with **opposing**
/// normals is the reverse: it separates the two solids, so union and
/// intersection drop it (it is interior to the result) while difference
/// keeps `a`'s copy, since that is exactly the surface where `a` is left
/// open by removing `b`.
fn keeps(op: BooleanOp, class: FaceClassification, from_a: bool) -> Keep {
    use FaceClassification::*;
    match (op, class) {
        // Union: the result is bounded by whatever is outside the other.
        (BooleanOp::Union, Outside) => Keep::AsIs,
        (BooleanOp::Union, OnSameNormal) if from_a => Keep::AsIs,
        (BooleanOp::Union, Inside | OnSameNormal | OnOppositeNormal) => Keep::Drop,

        // Intersection: bounded by whatever is inside the other.
        (BooleanOp::Intersection, Inside) => Keep::AsIs,
        (BooleanOp::Intersection, OnSameNormal) if from_a => Keep::AsIs,
        (BooleanOp::Intersection, Outside | OnSameNormal | OnOppositeNormal) => Keep::Drop,

        // Difference (a - b) = a intersected with the complement of b, so
        // `a`'s faces behave as for intersection-with-the-outside, and
        // `b`'s surviving faces are the ones inside `a`, turned around to
        // face into the cavity they now bound.
        (BooleanOp::Difference, Outside | OnOppositeNormal) if from_a => Keep::AsIs,
        (BooleanOp::Difference, Inside) if !from_a => Keep::Reversed,
        (BooleanOp::Difference, Inside | Outside | OnSameNormal | OnOppositeNormal) => Keep::Drop,

        // Symmetric difference: (a - b) and (b - a), each bounded by its
        // own faces outside the other and the other's faces inside it,
        // turned around. Where the two lie along each other, the material
        // is on one side of both (in both, so in neither piece) or on
        // either side (in one each, so inside the result): no boundary.
        (BooleanOp::SymmetricDifference, Outside) => Keep::AsIs,
        (BooleanOp::SymmetricDifference, Inside) => Keep::Reversed,
        (BooleanOp::SymmetricDifference, OnSameNormal | OnOppositeNormal) => Keep::Drop,
    }
}

/// Fixed seed for the ray casting behind every containment query here. Both
/// `face_contains` and `shell_contains` retry until they find a ray grazing
/// nothing, so their answers are seed-independent; a constant keeps a boolean
/// reproducible run to run.
const SEED: u64 = 0xB001_EA47_0000_0001;

/// `solid_a` combined with `solid_b` under `op`, as a new solid in `model`.
///
/// `Ok(None)` means the result is **empty**, which is an answer rather than a
/// failure: intersecting two solids that do not overlap, or subtracting a
/// solid that wholly contains the first, legitimately leaves nothing. Callers
/// that treat "no solid" as an error would reject the majority of scenes in
/// this crate's own test set.
///
/// Both input solids are consumed either way: their faces are transferred to
/// the result or deleted, so neither id is valid afterwards.
///
/// The result is named `namer`'s root name, and everything the boolean
/// creates on the way is named as `naming` describes; every face,
/// edge and vertex it keeps keeps its name.
pub fn boolean(
    model: &mut NamedModel,
    namer: &Namer,
    solid_a: SolidId,
    solid_b: SolidId,
    op: BooleanOp,
    params: RemeshParams,
) -> GeopResult<Option<SolidId>> {
    combine(model, namer, solid_a, solid_b, op, None, params)?.stopped()
}

/// A tool going up to the next face — the second operand of a boolean — by
/// the names its start and end faces had before it (see `naming`),
/// and whether only the piece of it taken is kept, `alone`, as a solid of
/// its own, rather than combined with the first operand.
#[derive(Clone, Copy, Debug)]
struct ToolEnds<'a> {
    start: &'a str,
    end: &'a str,
    alone: bool,
}

/// What combining a tool up to the next face came to: stopped all round,
/// and what that left — or not, the target meeting part of the profile, or
/// none, and the rest going on past it to the tool's end. A caller may then
/// go another way (see `combine_sides` of the extrude and revolve); the
/// error says where the tool ran on, for one that does not.
#[derive(Debug)]
pub enum Reach<T> {
    Stopped(T),
    RunsOn(GeopError),
}

impl<T> Reach<T> {
    /// What it left, if stopped all round, else why not as an error.
    pub fn stopped(self) -> GeopResult<T> {
        match self {
            Reach::Stopped(t) => Ok(t),
            Reach::RunsOn(e) => Err(e),
        }
    }
}

/// Like [`boolean`] — `target` combined with `tool` — but "up to next":
/// only the first piece of `tool` its face `start` reaches is combined.
///
/// For a [`BooleanOp::Union`] that is the first piece of the tool outside
/// the target — growing from `start` until it runs into the target; for a
/// [`BooleanOp::Difference`] or [`BooleanOp::Intersection`], the first
/// piece inside it — cutting, or keeping, from where the tool enters the
/// target until it comes out of it again. The tool has to be long enough to
/// reach past whatever stops it: if the piece also reaches its face `end`,
/// nothing stopped it, and it [`Reach::RunsOn`]. `start` and `end` are the
/// names those faces had before the boolean (see `naming`).
pub fn boolean_up_to_next(
    model: &mut NamedModel,
    namer: &Namer,
    target: SolidId,
    tool: SolidId,
    op: BooleanOp,
    (start, end): (&str, &str),
    params: RemeshParams,
) -> GeopResult<Reach<Option<SolidId>>> {
    let up_to_next = ToolEnds {
        start,
        end,
        alone: false,
    };
    combine(model, namer, target, tool, op, Some(up_to_next), params)
}

/// The first piece of `tool` its face `start` reaches outside `stops` — up
/// to the next face of any of them — as a solid of its own, named
/// `namer`'s root, `stops` left as they are: what a new body going up to the
/// next face is. [`Reach::RunsOn`], as [`boolean_up_to_next`] does, if
/// nothing stops the tool before its face `end`.
///
/// The tool is cut by a copy of the solids it stops at, so that only the
/// copy is imprinted and consumed.
pub fn piece_up_to_next(
    model: &mut NamedModel,
    namer: &Namer,
    tool: SolidId,
    stops: &[SolidId],
    (start, end): (&str, &str),
    params: RemeshParams,
) -> GeopResult<Reach<SolidId>> {
    let ctx = with_context!(
        "piece_up_to_next(name={}, tool={tool}, stops={stops:?})",
        namer.root()
    );
    let copy = copy_of(model, namer, stops).with_context(&ctx)?;
    // As a join of the tool to the copy, picking only the piece — the copy
    // first, as in every other combination with a tool.
    let up_to_next = ToolEnds {
        start,
        end,
        alone: true,
    };
    let reach = combine(
        model,
        namer,
        copy,
        tool,
        BooleanOp::Union,
        Some(up_to_next),
        params,
    )
    .with_context(&ctx)?;
    Ok(match reach {
        Reach::Stopped(piece) => Reach::Stopped(
            piece.ok_or_else(|| ctx(GeopError::new("up to next: nothing of the tool is left")))?,
        ),
        Reach::RunsOn(e) => Reach::RunsOn(ctx(e)),
    })
}

/// A copy of the solids `solids` as one solid that shares nothing with
/// them, named by `namer` as a copy: what an operation that leaves them as
/// they are imprints and consumes instead.
pub(crate) fn copy_of(
    model: &mut NamedModel,
    namer: &Namer,
    solids: &[SolidId],
) -> GeopResult<SolidId> {
    let mut faces = Vec::new();
    for &solid in solids {
        faces.extend(model.topology().solid_faces(solid)?);
    }
    let copy = model.copy_faces(&faces, Some(namer.name(&["copy"])), |name| {
        namer.name(&["copy", name])
    })?;
    Ok(copy.solid.expect("built as a solid"))
}

/// [`boolean`], or with `up_to_next` [`boolean_up_to_next`] — or, with the
/// tool kept alone, [`piece_up_to_next`].
fn combine(
    model: &mut NamedModel,
    namer: &Namer,
    solid_a: SolidId,
    solid_b: SolidId,
    op: BooleanOp,
    up_to_next: Option<ToolEnds>,
    params: RemeshParams,
) -> GeopResult<Reach<Option<SolidId>>> {
    let ctx = with_context!(
        "boolean(name={}, solid_a={solid_a}, solid_b={solid_b}, op={op:?}, up_to_next={up_to_next:?})",
        namer.root()
    );
    let (parted, between_two) = given_topology(model, &[solid_a, solid_b]).with_context(&ctx)?;

    // A remesh error names entities by id, and ids say nothing of what the
    // entities are: name them, as they are where it failed.
    let origins = remesh(model, namer, solid_a, solid_b, params)
        .map_err(|e| {
            let names = model.names_mentioned(&e.to_string());
            e.with_context(format!("where it names {}", names.join(", ")))
        })
        .with_context(&ctx)?;
    let topology = model.topology();

    let faces_a = topology.solid_faces(solid_a).with_context(&ctx)?;
    let faces_b = topology.solid_faces(solid_b).with_context(&ctx)?;

    let mut decisions: HashMap<FaceId, (FaceClassification, Keep)> = HashMap::new();
    let mut keep: Vec<FaceId> = Vec::new();
    let mut reverse: Vec<FaceId> = Vec::new();
    let mut classified = Classified::default();
    // Every face that could meet the other solid's boundary is classified
    // at its own interior points, each checked against the others: a face
    // a missed intersection curve leaves straddling is named. A face whose
    // trimmed patch is apart from every one of the other solid's cannot
    // straddle it, and lies on the side of any face it shares an edge with
    // that is inside or outside, all of the edge being apart from the
    // other solid's boundary too. Those take their side from such a
    // neighbour, and only a set of them no classified face borders is
    // classified, at one face: what a cut leaves untouched — all but a few
    // faces of a gear once it has many teeth — costs no ray.
    let (patches_a, patches_b) = (
        trimmed_patches(topology, &faces_a).with_context(&ctx)?,
        trimmed_patches(topology, &faces_b).with_context(&ctx)?,
    );
    for (patches, from_a, other, others) in [
        (&patches_a, true, solid_b, &patches_b),
        (&patches_b, false, solid_a, &patches_a),
    ] {
        let classify = |face_id: FaceId, classified: &mut Classified| {
            classify_face(topology, face_id, other, params, classified)
                .map_err(|e| {
                    let names = model.names_mentioned(&e.to_string());
                    e.with_context(format!("where it names {}", names.join(", ")))
                })
                .with_context(&ctx)
                .with_context(with_context!(
                    "classifying face {face_id} ({:?}, of {})",
                    model.name_of(face_id),
                    if from_a {
                        "the first solid"
                    } else {
                        "the second solid"
                    }
                ))
        };
        let mut apart = std::collections::BTreeSet::new();
        let mut sided = std::collections::VecDeque::new();
        for (face_id, patch) in patches {
            if others.iter().any(|(_, o)| patch.could_overlap(o)) {
                let class = classify(*face_id, &mut classified)?;
                decisions.insert(*face_id, (class, keeps(op, class, from_a)));
                sided.push_back((*face_id, class));
            } else {
                apart.insert(*face_id);
            }
        }
        loop {
            while let Some((face_id, class)) = sided.pop_front() {
                if !matches!(
                    class,
                    FaceClassification::Inside | FaceClassification::Outside
                ) {
                    continue;
                }
                for edge in topology.face_edges(face_id) {
                    for coedge in topology.coedges_of_edge(edge) {
                        let next = topology.get_coedge(coedge).with_context(&ctx)?.face;
                        if apart.remove(&next) {
                            decisions.insert(next, (class, keeps(op, class, from_a)));
                            sided.push_back((next, class));
                        }
                    }
                }
            }
            let Some(face_id) = apart.pop_first() else {
                break;
            };
            let class = classify(face_id, &mut classified)?;
            decisions.insert(face_id, (class, keeps(op, class, from_a)));
            sided.push_back((face_id, class));
        }
    }
    if let Some(ToolEnds { start, end, alone }) = up_to_next {
        let tool: HashSet<FaceId> = faces_b.iter().copied().collect();
        // The piece of the kind the operator takes from the target, and
        // what of the target the tool leaves untouched — but for an
        // intersection, which keeps only what both hold, or a piece kept
        // alone.
        let (wanted_inside, keep_untouched) = match op {
            BooleanOp::Union => (false, !alone),
            BooleanOp::Difference => (true, true),
            BooleanOp::Intersection => (true, false),
            BooleanOp::SymmetricDifference => {
                return Err(ctx(GeopError::new(
                    "up to next: a symmetric difference keeps every piece of the tool, so it reaches no next face; use a blind extent",
                )));
            }
        };
        let reached = match reach_up_to_next(
            topology,
            &tool,
            (wanted_inside, keep_untouched),
            &origins.faces,
            (start, end),
            &mut decisions,
        )
        .with_context(&ctx)?
        {
            Reach::Stopped(reached) => reached,
            Reach::RunsOn(e) => return Ok(Reach::RunsOn(ctx(e))),
        };
        if alone {
            keep_alone(&reached, &mut decisions);
        }
    }
    // In the order the faces were classified, so the result's shell lists
    // its faces the same way every run.
    for &face_id in faces_a.iter().chain(&faces_b) {
        match decisions[&face_id].1 {
            Keep::AsIs => keep.push(face_id),
            Keep::Reversed => {
                keep.push(face_id);
                reverse.push(face_id);
            }
            Keep::Drop => {}
        }
    }
    check_closed(model, &decisions).with_context(&ctx)?;

    for &face_id in &reverse {
        model.reverse_face(face_id).with_context(&ctx)?;
    }

    let Some(solid) = model
        .assemble_solid(&[solid_a.into(), solid_b.into()], &keep, namer.root())
        .with_context(&ctx)?
    else {
        return Ok(Reach::Stopped(None));
    };
    // What the boolean left in two is made one where it lies on one surface:
    // pieces of one face side by side, split along an edge of the other
    // solid that is no edge of the result, and a face of each solid running
    // on into the other's. What the operands already had stays as they were
    // given: an edge that parted the same two faces there, a cylinder's
    // quarters or a face a step split on purpose, and a vertex that was
    // already only between two edges.
    let by_name = origins_by_name(model, &origins, solid).with_context(&ctx)?;
    let stays = edges_parting_as_given(model, solid, &origins, &parted).with_context(&ctx)?;
    let merged = model
        .merge_faces_on_one_surface(
            solid.into(),
            |e| !stays.contains(&e),
            |v| !between_two.contains(&v),
        )
        .with_context(&ctx)?;
    match merged {
        geop_core_topology::Body::Solid(solid) => {
            name_lone_pieces(model, &by_name).with_context(&ctx)?;
            Ok(Reach::Stopped(Some(solid)))
        }
        geop_core_topology::Body::Sheet(_) => Err(ctx(GeopError::new(
            "merging the faces of the result on one surface left no solid",
        ))),
    }
}

/// See [`given_topology`].
type GivenTopology = (HashMap<String, [String; 2]>, HashSet<VertexId>);

/// What of `solids` a boolean leaves as it was given: per edge, by name, the
/// names of the two faces it parts, sorted; and the vertices only two edges
/// meet at.
fn given_topology(model: &NamedModel, solids: &[SolidId]) -> GeopResult<GivenTopology> {
    let topology = model.topology();
    let mut parted: HashMap<String, Vec<String>> = HashMap::new();
    let mut edges_at: HashMap<VertexId, HashSet<EdgeId>> = HashMap::new();
    for &solid in solids {
        for face in topology.solid_faces(solid)? {
            let name = model.name(face)?;
            for edge in topology.face_edges(face) {
                parted
                    .entry(model.name(edge)?.to_string())
                    .or_default()
                    .push(name.to_string());
                let e = topology.get_edge(edge)?;
                for v in [e.start_vertex, e.end_vertex] {
                    edges_at.entry(v).or_default().insert(edge);
                }
            }
        }
    }
    let parted = (parted.into_iter())
        .filter_map(|(edge, mut faces)| {
            faces.sort();
            let [f, g] = <[String; 2]>::try_from(faces).ok()?;
            Some((edge, [f, g]))
        })
        .collect();
    let between_two = (edges_at.into_iter())
        .filter(|(_, edges)| edges.len() == 2)
        .map(|(v, _)| v)
        .collect();
    Ok((parted, between_two))
}

/// The edges of `solid`, a boolean's result, that part pieces of the two
/// faces the edge they are a piece of parted in the operands (`parted`, see
/// [`given_topology`]): an edge the operands had between those faces, which
/// the boolean did not make. An edge of one operand that now parts a face
/// of each is the boolean's seam, whatever curve it runs along.
fn edges_parting_as_given(
    model: &NamedModel,
    solid: SolidId,
    origins: &Origins,
    parted: &HashMap<String, [String; 2]>,
) -> GeopResult<HashSet<EdgeId>> {
    let topology = model.topology();
    let mut faces_of: HashMap<EdgeId, Vec<&String>> = HashMap::new();
    for face in topology.solid_faces(solid)? {
        if let Some(origin) = origins.faces.get(&face) {
            for edge in topology.face_edges(face) {
                faces_of.entry(edge).or_default().push(origin);
            }
        }
    }
    Ok((faces_of.into_iter())
        .filter_map(|(edge, mut faces)| {
            let given = parted.get(origins.edges.get(&edge)?)?;
            faces.sort();
            (faces.len() == 2 && *faces[0] == given[0] && *faces[1] == given[1]).then_some(edge)
        })
        .collect())
}

/// Narrows `decisions` down to the piece of the tool (the faces `tool`) that
/// its start reaches first, see [`boolean_up_to_next`].
///
/// The target's boundary cuts the tool's material into pieces, alternately
/// outside and inside the target. Each piece is bounded by tool faces —
/// all on one side of the target, the side its classification says — and
/// by target faces inside the tool, each of which bounds two pieces: the
/// inside one on its material's side, the outside one on the other. (Remesh
/// split every face along every crossing, so no face straddles two pieces.)
/// Faces sharing an edge bound the same piece if their material lies on the
/// same side of the target, which makes the pieces connected sets of face
/// *sides*, and two pieces neighbours when a target face lies between them.
///
/// From the piece bounded by the start face, the piece wanted is the first
/// of the kind wanted — outside the target, growing until it runs into the
/// target, or inside it, cutting or keeping from where it enters the target
/// until it leaves it. That is the start's own piece, if it is of that
/// kind, or else the pieces next to it. The tool's other pieces are
/// dropped, and the target's faces between them kept if `keep_untouched`,
/// as if the tool had not been there, or else dropped too.
fn reach_up_to_next(
    model: &Model,
    tool: &HashSet<FaceId>,
    (wanted_inside, keep_untouched): (bool, bool),
    origins: &HashMap<FaceId, String>,
    (start, end): (&str, &str),
    decisions: &mut HashMap<FaceId, (FaceClassification, Keep)>,
) -> GeopResult<Reach<Reached>> {
    use FaceClassification::*;
    // Which side of the target the material behind a tool face lies on.
    let inside = |class: FaceClassification| matches!(class, Inside | OnSameNormal);
    let mut tool_faces: Vec<FaceId> = tool.iter().copied().collect();
    tool_faces.sort_by_key(|f| f.0);
    let mut between: Vec<FaceId> = decisions
        .iter()
        .filter(|(face, (class, _))| !tool.contains(face) && *class == Inside)
        .map(|(&face, _)| face)
        .collect();
    between.sort_by_key(|f| f.0);
    // Nodes: every tool face, then each target face between pieces by its
    // inner side, then by its outer side.
    let (n, m) = (tool_faces.len(), between.len());
    let (inner, outer) = (|j: usize| n + j, |j: usize| n + m + j);
    let tool_index: HashMap<FaceId, usize> = tool_faces
        .iter()
        .enumerate()
        .map(|(i, &f)| (f, i))
        .collect();
    let between_index: HashMap<FaceId, usize> =
        between.iter().enumerate().map(|(j, &f)| (f, j)).collect();
    let class = |face: &FaceId| decisions[face].0;

    let mut users: HashMap<EdgeId, Vec<FaceId>> = HashMap::new();
    for &face in tool_faces.iter().chain(&between) {
        for edge in model.face_edges(face) {
            users.entry(edge).or_default().push(face);
        }
    }
    let mut pieces = UnionFind::new(n + 2 * m);
    for faces in users.values() {
        let (mut in_node, mut out_node) = (None, None);
        let join = |slot: &mut Option<usize>, node: usize, pieces: &mut UnionFind| match *slot {
            Some(other) => pieces.union(other, node),
            None => *slot = Some(node),
        };
        for face in faces {
            if let Some(&i) = tool_index.get(face) {
                let slot = if inside(class(face)) {
                    &mut in_node
                } else {
                    &mut out_node
                };
                join(slot, i, &mut pieces);
            } else {
                let j = between_index[face];
                join(&mut in_node, inner(j), &mut pieces);
                join(&mut out_node, outer(j), &mut pieces);
            }
        }
    }
    // A tool face's piece is on its material's side; a target face's inner
    // side is in an inside piece, its outer side in an outside one.
    let piece_inside = |node: usize| match node {
        i if i < n => inside(class(&tool_faces[i])),
        i => i < n + m,
    };
    let mut kind: HashMap<usize, bool> = HashMap::new();
    for node in 0..n + 2 * m {
        kind.insert(pieces.find(node), piece_inside(node));
    }
    let mut neighbours: HashMap<usize, HashSet<usize>> = HashMap::new();
    for j in 0..m {
        let (a, b) = (pieces.find(inner(j)), pieces.find(outer(j)));
        neighbours.entry(a).or_default().insert(b);
        neighbours.entry(b).or_default().insert(a);
    }

    let origin_is = |face: &FaceId, name: &str| origins.get(face).is_some_and(|o| o == name);
    // The start may lie partly on either side of the target — a rib's
    // profile running on into the walls it stands between. Its pieces of
    // the kind wanted are the first ones it reaches; only where it has none
    // are the pieces next to it the first. Taking the pieces next to its
    // other fragments too would take a second piece, reached through the
    // target: the rest of the tool on the walls' far side.
    let start_pieces: HashSet<usize> = tool_faces
        .iter()
        .enumerate()
        .filter(|(_, face)| origin_is(face, start))
        .map(|(i, _)| pieces.find(i))
        .collect();
    let mut reached: HashSet<usize> = start_pieces
        .iter()
        .copied()
        .filter(|p| kind[p] == wanted_inside)
        .collect();
    if reached.is_empty() {
        for piece in &start_pieces {
            let next = neighbours.get(piece).into_iter().flatten();
            reached.extend(next.filter(|p| kind[p] == wanted_inside));
        }
    }
    if reached.is_empty() {
        let starts: Vec<String> = tool_faces
            .iter()
            .enumerate()
            .filter(|(_, face)| origin_is(face, start))
            .map(|(i, face)| {
                let piece = pieces.find(i);
                format!(
                    "{face} ({:?}, in an {} piece with {} neighbour(s))",
                    class(face),
                    if kind[&piece] { "inside" } else { "outside" },
                    neighbours.get(&piece).map_or(0, HashSet::len)
                )
            })
            .collect();
        let mut tool_origins: Vec<&str> = tool_faces
            .iter()
            .filter_map(|f| origins.get(f).map(String::as_str))
            .collect();
        tool_origins.sort_unstable();
        tool_origins.dedup();
        return Err(GeopError::new(format!(
            "up to next: from its start {start:?}, the profile meets nothing of the target to {} — does it go the right way? (the start's faces: [{}]; {} piece(s), {m} target face(s) between them; the tool's faces come from {tool_origins:?})",
            if wanted_inside {
                "cut into"
            } else {
                "grow up to"
            },
            starts.join(", "),
            kind.len(),
        )));
    }
    for (i, face) in tool_faces.iter().enumerate() {
        let in_reach = reached.contains(&pieces.find(i));
        if in_reach && origin_is(face, end) {
            return Ok(Reach::RunsOn(GeopError::new(format!(
                "up to next: nothing stops the profile: it reaches the end {end:?} without meeting a face of the target all around"
            ))));
        }
        if !in_reach {
            decisions.get_mut(face).expect("classified").1 = Keep::Drop;
        }
    }
    let mut bounding = Reached {
        tool: HashSet::new(),
        target: HashSet::new(),
    };
    for (i, face) in tool_faces.iter().enumerate() {
        if reached.contains(&pieces.find(i)) {
            bounding.tool.insert(*face);
        }
    }
    for (j, face) in between.iter().enumerate() {
        let facing = if wanted_inside { inner(j) } else { outer(j) };
        if reached.contains(&pieces.find(facing)) {
            bounding.target.insert(*face);
        } else {
            // Untouched by the piece taking part: the target's boundary
            // there stays — unless only what the tool reaches is kept.
            decisions.get_mut(face).expect("classified").1 = if keep_untouched {
                Keep::AsIs
            } else {
                Keep::Drop
            };
        }
    }
    Ok(Reach::Stopped(bounding))
}

/// What [`reach_up_to_next`] found the tool's taken piece bounded by: the
/// tool's faces, and the target's faces between it and the rest.
struct Reached {
    tool: HashSet<FaceId>,
    target: HashSet<FaceId>,
}

/// Turns the decisions of a union of a target with the piece of a tool it
/// `reached` into those of the piece alone: its own faces as they are —
/// where it lies on the target too, its start drawn on a face of the target,
/// included — and the target's faces bounding it turned around to face out
/// of it. Nothing else of either.
fn keep_alone(reached: &Reached, decisions: &mut HashMap<FaceId, (FaceClassification, Keep)>) {
    use FaceClassification::*;
    for (face, (class, keep)) in decisions.iter_mut() {
        *keep = if reached.tool.contains(face) && matches!(class, Outside | OnOppositeNormal) {
            Keep::AsIs
        } else if reached.target.contains(face) {
            Keep::Reversed
        } else {
            Keep::Drop
        };
    }
}

/// Fails unless the kept faces close up: every edge they use must be used an
/// even number of times by them — twice where two faces meet, four times
/// where the result touches itself along the edge, which a solid may (see
/// `validation::manifold_check`). An odd count is a hole in
/// the result, so a face was kept or dropped wrongly, or remesh left a face
/// straddling the other solid; the error lists every face along that edge
/// with its name and classification, which tells the two apart.
fn check_closed(
    model: &NamedModel,
    decisions: &HashMap<FaceId, (FaceClassification, Keep)>,
) -> GeopResult<()> {
    let topology = model.topology();
    let mut uses: HashMap<EdgeId, usize> = HashMap::new();
    for (&face_id, &(_, decision)) in decisions {
        if decision == Keep::Drop {
            continue;
        }
        for edge in topology.face_edges(face_id) {
            *uses.entry(edge).or_default() += 1;
        }
    }
    let Some((&edge, &count)) = uses
        .iter()
        .filter(|(_, n)| *n % 2 != 0)
        .min_by_key(|(e, _)| e.0)
    else {
        return Ok(());
    };
    let e = topology.get_edge(edge)?;
    let faces: Vec<String> = topology
        .coedges_of_edge(edge)
        .into_iter()
        .map(|coedge| {
            let face = topology.get_coedge(coedge)?.face;
            let name = model.name_of(face).unwrap_or("?");
            let decision = match decisions.get(&face) {
                Some((class, decision)) => format!("{class:?} -> {decision:?}"),
                None => "not in either solid".to_string(),
            };
            Ok(format!(
                "{face} {name:?}: {decision}, bounded by {:?}",
                topology.face_corners(face)?
            ))
        })
        .collect::<GeopResult<_>>()?;
    Err(GeopError::new(format!(
        "boolean: the kept faces use edge {edge} (from {:?} to {:?}) {count} time(s), leaving the result open there; faces along it: [{}]",
        topology.get_vertex(e.start_vertex)?.point,
        topology.get_vertex(e.end_vertex)?.point,
        faces.join(", ")
    )))
}

/// Each of `faces` with the part of its surface over its trim's parameter
/// box: the box around its pcurves' control points, which holds the loops
/// and so all the face (a face bounded by a bare vertex is its whole
/// surface). Its control points' box holds the face, and is apart from
/// another's where the faces are, however large the surface they were cut
/// from — a gear's cap cut by every tooth is one surface.
fn trimmed_patches(model: &Model, faces: &[FaceId]) -> GeopResult<Vec<(FaceId, NurbSurface3D)>> {
    faces
        .iter()
        .map(|&face_id| {
            let surface = &model.get_face(face_id)?.surface;
            let mut bounds: Option<[Scalar; 2]> = None;
            for coedge in model.iterate_face_coedges(face_id) {
                for p in &model.get_coedge(coedge)?.pcurve.control_points {
                    let uv = [p[0].div(p[2])?, p[1].div(p[2])?];
                    bounds = Some(match bounds {
                        None => uv,
                        Some(b) => [b[0].union(uv[0]), b[1].union(uv[1])],
                    });
                }
            }
            let patch = match bounds {
                Some([u, v]) => {
                    surface.sub_surface((u.lower(), u.upper()), (v.lower(), v.upper()))?
                }
                None => surface.clone(),
            };
            Ok((face_id, patch))
        })
        .collect()
}

thread_local! {
    static FACES_CLASSIFIED: std::cell::Cell<usize> = const { std::cell::Cell::new(0) };
}

/// How many faces this thread's booleans have classified at points of
/// their own, against the other solid, ever: what classifying costs grows
/// with this, so a test can bound the work of a cut by the faces it
/// touches rather than time it.
pub fn faces_classified() -> usize {
    FACES_CLASSIFIED.get()
}

/// Where `face_id` sits relative to `other_solid`, decided at points strictly
/// inside the face's trimmed region.
///
/// One point would be enough *because remesh ran first*: every curve along
/// which the other solid's boundary crosses this face has been imprinted, and
/// the face split along it, so the face no longer straddles anything. Without
/// that guarantee this would be unsound, which is why it lives here rather
/// than as a general-purpose query. Every interior point offered (one or two
/// per boundary coedge, see `face_interior_point_where`) is classified anyway,
/// and points that disagree are an error: a missed intersection curve then
/// names the face it left straddling, instead of surfacing as a result that
/// is open somewhere else — or, by luck of which point came first, not at all.
///
/// A point that lands *on* the other solid's boundary doesn't decide it by
/// itself: the face may lie along that boundary over an area (coincident
/// patches, told apart by their normals), or merely touch it at a point or
/// along a curve — a cube face resting on a sphere's pole, say, where the
/// sphere has no normal at all. So such a point is set aside and the next
/// interior point tried (one per boundary coedge, see
/// `face_interior_point_where`); any point off the other solid's boundary
/// classifies the whole face. Only if every one lies on it are the patches
/// coincident, and their normals are compared.
///
/// Each point is classified from the nearest point `classified` already
/// holds — of this face or one classified before — along the segment
/// between them (see `shell_contains_from`), and added to it: one ray
/// across the whole solid for the first point of a boolean, and short ones
/// after it, instead of one across the solid for every point.
pub(crate) fn classify_face(
    model: &Model,
    face_id: FaceId,
    other_solid: SolidId,
    params: RemeshParams,
    classified: &mut Classified,
) -> GeopResult<FaceClassification> {
    FACES_CLASSIFIED.set(FACES_CLASSIFIED.get() + 1);
    let face = model.get_face(face_id)?;
    let mut decided: Vec<(FaceClassification, Vector3)> = Vec::new();
    let mut on_boundary = Vec::new();
    let (max_nodes, epsilon) = (params.max_nodes, params.curve_curve_min_subdivision_size);
    face_interior_point_where(model, face_id, max_nodes, epsilon, SEED, |u, v| {
        let point = face.surface.evaluate(u, v)?;
        match solid_contains(
            model,
            other_solid,
            point,
            max_nodes,
            epsilon,
            SEED,
            classified,
        )? {
            ShellPoint::Inside => decided.push((FaceClassification::Inside, point)),
            ShellPoint::Outside => decided.push((FaceClassification::Outside, point)),
            ShellPoint::OnFace | ShellPoint::OnEdge | ShellPoint::OnVertex => {
                on_boundary.push((u, v, point));
            }
        }
        Ok(false)
    })?;
    if let Some(&(classification, point)) = decided.first() {
        if let Some((other, other_point)) = decided.iter().find(|(c, _)| *c != classification) {
            // The face's boundary: which curves were imprinted into it, and
            // so — between the two points — which one was not.
            let loops = model.face_corners(face_id)?;
            // Where the missing curve may run from: the other solid's edges
            // at the face's corners that do not bound it, and what the
            // edge x face search makes of each against the face's surface —
            // an edge lying on the face that was not imprinted shows here.
            let other_edges: HashSet<EdgeId> = model.iter_body_edges(other_solid)?.collect();
            let bounding: HashSet<EdgeId> = model.face_edges(face_id).collect();
            let mut corners: Vec<_> = loops.iter().flatten().map(|&(v, _)| v).collect();
            corners.sort();
            corners.dedup();
            let mut leaving = Vec::new();
            for vertex in corners {
                for edge in model.edges_at_vertex(vertex) {
                    if !other_edges.contains(&edge) || bounding.contains(&edge) {
                        continue;
                    }
                    let found = curve_surface_intersect(
                        &model.get_edge(edge)?.curve,
                        &face.surface,
                        params.max_edge_intersections,
                        max_nodes,
                        epsilon,
                    );
                    leaving.push(format!("{edge} at {vertex}: {found:?}"));
                }
            }
            return Err(GeopError::new(format!(
                "classify_face: face {face_id} straddles solid {other_solid}'s boundary — {point:?} is {classification:?}, {other_point:?} is {other:?} — so remesh left an intersection curve unimprinted; the face is bounded by (corners, outer loop first) {loops:?}; the other solid's edges at its corners, against its surface: {leaving:?}"
            )));
        }
        return Ok(classification);
    }

    // Every interior point lies on the other solid's boundary: the patches
    // coincide, and which way the shared patch faces is what separates "these
    // solids touch" from "these solids overlap along this patch". Any point
    // of the shared patch shows that equally well, so the first where both
    // normals exist is used — a revolve's cap collapses to its pole, where
    // it has none.
    let mut undefined = None;
    for &(u, v, point) in &on_boundary {
        let normals = face
            .surface
            .normal(u, v)
            .and_then(|n| Ok((n, solid_normal_at(model, other_solid, &point, params)?)));
        let (this_normal, other_normal) = match normals {
            Ok(normals) => normals,
            Err(e) => {
                undefined = Some(e);
                continue;
            }
        };
        let alignment = this_normal.prod_dot(&other_normal);
        return if alignment.definitely_greater(Scalar::ZERO) {
            Ok(FaceClassification::OnSameNormal)
        } else if alignment.definitely_less(Scalar::ZERO) {
            Ok(FaceClassification::OnOppositeNormal)
        } else {
            Err(GeopError::new(format!(
                "boolean: face {face_id} lies on solid {other_solid}'s boundary at {point:?} (its interior point uv=({u:?}, {v:?})), but the two normals ({this_normal:?} and {other_normal:?}) are too close to perpendicular to tell which side is which"
            )))
        };
    }
    Err(match undefined {
        Some(e) => e.with_context(format!(
            "classify_face: face {face_id} lies on solid {other_solid}'s boundary at each of its {} interior points tried, and no normal comparison could be made at any of them",
            on_boundary.len()
        )),
        None => GeopError::new(format!(
            "classify_face: face {face_id} yielded interior points, yet none was classified or set aside"
        )),
    })
}

/// The normal of whichever face of `solid_id` contains `point`.
fn solid_normal_at(
    model: &Model,
    solid_id: SolidId,
    point: &Vector3,
    params: RemeshParams,
) -> GeopResult<Vector3> {
    let at = solid_face_at(
        model,
        solid_id,
        point,
        params.max_nodes,
        params.curve_curve_min_subdivision_size,
        SEED,
    )?;
    match at {
        Some((face_id, u, v)) => model.get_face(face_id)?.surface.normal(u, v),
        None => Err(GeopError::new(format!(
            "boolean: no face of solid {solid_id} contains {point:?}, although the solid reported the point on its boundary"
        ))),
    }
}

#[cfg(test)]
mod tests {
    use super::{BooleanOp, Classified, FaceClassification, boolean, classify_face};
    use crate::{
        remesh::RemeshParams,
        scenes::{scene, validation},
    };
    use geop_core_math::{scalars::Scalar, vector::Vector3};
    use geop_core_topology::{contains::rng::Rng, validation::validate_fast};
    use geop_lib_sweep::shapes::{Axis, block, v3};
    use geop_ops_model::Namer;

    /// A fresh operation id, so that every solid and boolean a test builds
    /// gets names of its own.
    fn fresh_id() -> String {
        static NEXT: std::sync::atomic::AtomicUsize = std::sync::atomic::AtomicUsize::new(0);
        format!(
            "op{}",
            NEXT.fetch_add(1, std::sync::atomic::Ordering::Relaxed)
        )
    }

    fn namer() -> Namer {
        Namer::new("boolean", &fresh_id()).unwrap()
    }

    /// Every face of an un-remeshed cube is outside a cylinder that only
    /// overlaps part of it — the classification itself, before any operator.
    #[test]
    fn classify_face_separates_inside_from_outside() {
        let mut s = scene("box_cylinder_drilled_hole_through");
        let params = RemeshParams::default();
        crate::remesh::remesh(&mut s.model, &namer(), s.solid_a, s.solid_b, params).unwrap();

        let mut inside = 0;
        let mut outside = 0;
        for face_id in s.model.topology().solid_faces(s.solid_a).unwrap() {
            let mut classified = Classified::default();
            match classify_face(
                s.model.topology(),
                face_id,
                s.solid_b,
                params,
                &mut classified,
            )
            .unwrap()
            {
                FaceClassification::Inside => inside += 1,
                FaceClassification::Outside => outside += 1,
                _ => {}
            }
        }
        assert!(
            inside > 0 && outside > 0,
            "a cylinder drilled through a cube must leave cube faces on both sides: {inside} inside, {outside} outside"
        );
    }

    /// Every face must yield an interior point, or classification is
    /// meaningless — this is the part of the algorithm with no fallback.
    #[test]
    fn every_remeshed_face_has_an_interior_point() {
        let mut s = scene("box_cylinder_drilled_hole_through");
        let params = RemeshParams::default();
        crate::remesh::remesh(&mut s.model, &namer(), s.solid_a, s.solid_b, params).unwrap();

        for solid in [s.solid_a, s.solid_b] {
            for face_id in s.model.topology().solid_faces(solid).unwrap() {
                let (u, v) = geop_core_topology::contains::face::face_interior_point(
                    s.model.topology(),
                    face_id,
                    params.max_nodes,
                    params.curve_curve_min_subdivision_size,
                    1234,
                )
                .unwrap_or_else(|e| panic!("face {face_id}: {e}"));
                let _ = s
                    .model
                    .topology()
                    .get_face(face_id)
                    .unwrap()
                    .surface
                    .evaluate(u, v)
                    .unwrap();
            }
        }
    }

    fn check_op(name: &str, op: BooleanOp) {
        let mut s = scene(name);
        let params = RemeshParams::default();
        let result = boolean(&mut s.model, &namer(), s.solid_a, s.solid_b, op, params)
            .unwrap_or_else(|e| panic!("{name} {op:?}: {e}"))
            .unwrap_or_else(|| panic!("{name} {op:?}: result is empty"));

        assert!(
            !s.model.topology().solid_faces(result).unwrap().is_empty(),
            "{name} {op:?}: result has no faces"
        );
        if let Err(errors) = validate_fast(&validation(), s.model.topology()) {
            panic!(
                "{name} {op:?}: {} validate_fast error(s): {}",
                errors.len(),
                errors[0]
            );
        }
    }

    /// Where a probe point lands relative to the result solid.
    fn contains(
        model: &geop_core_topology::Model,
        solid: geop_core_topology::SolidId,
        p: (f64, f64, f64),
    ) -> bool {
        let params = RemeshParams::default();
        let point = geop_core_math::vector::Vector3::from_array([
            Scalar::from_f64(p.0),
            Scalar::from_f64(p.1),
            Scalar::from_f64(p.2),
        ]);
        let shell = model.get_solid(solid).unwrap().shells[0];
        matches!(
            geop_core_topology::contains::shell::shell_contains(
                model,
                shell,
                point,
                params.max_nodes,
                params.curve_curve_min_subdivision_size,
                0xA5A5_1234,
            )
            .unwrap(),
            geop_core_topology::contains::shell::PointClassification::Inside
        )
    }

    /// The cube is `[-0.5, 0.5]^3`; the cylinder has radius `0.2` on the z
    /// axis and spans `z in [-1, 1]`. So each probe below is unambiguously in
    /// one region, and the three operators must disagree about them in
    /// exactly the way their definitions say. Structural validity says
    /// nothing about *which* faces were kept — this is what does.
    const IN_CUBE_ONLY: (f64, f64, f64) = (0.4, 0.4, 0.0);
    const IN_CYLINDER_ONLY: (f64, f64, f64) = (0.0, 0.0, 0.8);
    const IN_BOTH: (f64, f64, f64) = (0.0, 0.0, 0.0);

    fn run(op: BooleanOp) -> (geop_core_topology::Model, geop_core_topology::SolidId) {
        let mut s = scene("box_cylinder_drilled_hole_through");
        let params = RemeshParams::default();
        let result = boolean(&mut s.model, &namer(), s.solid_a, s.solid_b, op, params)
            .unwrap()
            .expect("this scene's operands overlap, so no operator is empty");
        (s.model.topology().clone(), result)
    }

    /// Regression: two cubes meeting along a shared edge, where the traced
    /// intersection curve *is* that edge — already part of both faces'
    /// boundaries. Splicing it carved a zero-area sliver instead of splitting
    /// anything. The direction check was passing because the corrector
    /// returned a sharpened `(u, v)` sitting 6e-17 off the trim boundary,
    /// which `face_contains` read as `Inside`; see
    /// `predictor_corrector_step`.
    #[test]
    fn box_grid_n1p00_n0p50_n1p00_difference_succeeds() {
        let mut s = scene("box_grid_n1p00_n0p50_n1p00");
        let params = RemeshParams::default();
        boolean(
            &mut s.model,
            &namer(),
            s.solid_a,
            s.solid_b,
            BooleanOp::Difference,
            params,
        )
        .unwrap_or_else(|e| panic!("{e}"));
    }

    /// The defining property of a boolean, checked by sampling space rather
    /// than by inspecting topology: for a point `p`, membership in the result
    /// is a pure function of membership in the two operands.
    ///
    /// `union` keeps `p` iff either operand held it, `intersection` iff both
    /// did, `difference` iff `a` held it and `b` did not. Nothing about faces,
    /// coedges or orientation enters into it — which is exactly why it is
    /// worth checking. A result can be structurally perfect and still enclose
    /// the wrong region, if `BooleanOp::keeps` drops a face it should have
    /// kept or leaves one facing inward, and no structural check can see that.
    ///
    /// Points landing *on* a boundary are skipped rather than asserted about:
    /// membership there is genuinely ambiguous, and both the operands and the
    /// result may legitimately disagree about a surface they share.
    fn check_boolean_matches_point_membership(name: &str, op: BooleanOp, seed: u64) {
        let mut s = scene(name);
        let params = RemeshParams::default();
        let (solid_a, solid_b) = (s.solid_a, s.solid_b);

        // Sampled against the operands *before* the boolean consumes them.
        let mut rng = Rng::new(seed);
        let mut samples = Vec::new();
        while samples.len() < SAMPLE_COUNT {
            let p = Vector3::from_array([
                Scalar::from_f64(rng.next_f64() * 4.0 - 2.0),
                Scalar::from_f64(rng.next_f64() * 4.0 - 2.0),
                Scalar::from_f64(rng.next_f64() * 4.0 - 2.0),
            ]);
            let (Some(in_a), Some(in_b)) = (
                strictly_inside(s.model.topology(), solid_a, p),
                strictly_inside(s.model.topology(), solid_b, p),
            ) else {
                continue;
            };
            samples.push((p, in_a, in_b));
        }

        let result = boolean(&mut s.model, &namer(), solid_a, solid_b, op, params)
            .unwrap_or_else(|e| panic!("{name} {op:?}: {e}"));

        for (p, in_a, in_b) in samples {
            let expected = match op {
                BooleanOp::Union => in_a || in_b,
                BooleanOp::Intersection => in_a && in_b,
                BooleanOp::Difference => in_a && !in_b,
                BooleanOp::SymmetricDifference => in_a != in_b,
            };
            let actual = match result {
                Some(solid) => match strictly_inside(s.model.topology(), solid, p) {
                    Some(inside) => inside,
                    // On the result's own boundary — ambiguous, no claim.
                    None => continue,
                },
                // An empty result contains nothing.
                None => false,
            };
            assert_eq!(
                actual,
                expected,
                "{name} {op:?}: point {p:?} is {} solid A and {} solid B, so the result should {} contain it",
                if in_a { "inside" } else { "outside" },
                if in_b { "inside" } else { "outside" },
                if expected { "" } else { "not" }
            );
        }
    }

    /// How many points each membership check samples.
    const SAMPLE_COUNT: usize = 20;

    /// `Some(true)`/`Some(false)` for a point strictly inside/outside every
    /// shell of `solid`; `None` if it lands on a boundary, where membership is
    /// not a yes/no question.
    fn strictly_inside(
        model: &geop_core_topology::Model,
        solid: geop_core_topology::SolidId,
        p: Vector3,
    ) -> Option<bool> {
        let params = RemeshParams::default();
        let mut inside = false;
        for shell in model.get_solid(solid).ok()?.shells.clone() {
            match geop_core_topology::contains::shell::shell_contains(
                model,
                shell,
                p,
                params.max_nodes,
                params.curve_curve_min_subdivision_size,
                0x5A3D_1234,
            ) {
                Ok(geop_core_topology::contains::shell::PointClassification::Inside) => {
                    inside = true
                }
                Ok(geop_core_topology::contains::shell::PointClassification::Outside) => {}
                _ => return None,
            }
        }
        Some(inside)
    }

    /// Scenes with genuinely overlapping operands, so every operator has
    /// something to do and the samples land on both sides of each boundary.
    const MEMBERSHIP_SCENES: &[&str] = &[
        "box_cylinder_drilled_hole_through",
        "box_cylinder_blind_hole",
        "figure8_cylinder_through_neck",
    ];

    #[test]
    fn union_matches_point_membership() {
        for (i, name) in MEMBERSHIP_SCENES.iter().enumerate() {
            check_boolean_matches_point_membership(name, BooleanOp::Union, 0xB001 + i as u64);
        }
    }

    /// What is in one solid but not the other: the drilled box's material
    /// outside the cylinder and the cylinder's outside the box, as one solid
    /// whose two pieces meet along the rims.
    #[test]
    fn symmetric_difference_matches_point_membership() {
        for (i, name) in MEMBERSHIP_SCENES.iter().enumerate() {
            check_boolean_matches_point_membership(
                name,
                BooleanOp::SymmetricDifference,
                0xB201 + i as u64,
            );
        }
    }

    #[test]
    fn intersection_matches_point_membership() {
        for (i, name) in MEMBERSHIP_SCENES.iter().enumerate() {
            check_boolean_matches_point_membership(
                name,
                BooleanOp::Intersection,
                0xB101 + i as u64,
            );
        }
    }

    #[test]
    fn difference_matches_point_membership() {
        for (i, name) in MEMBERSHIP_SCENES.iter().enumerate() {
            check_boolean_matches_point_membership(name, BooleanOp::Difference, 0xB201 + i as u64);
        }
    }

    #[test]
    fn union_contains_either_operand() {
        let (model, r) = run(BooleanOp::Union);
        assert!(contains(&model, r, IN_CUBE_ONLY), "cube-only point");
        assert!(contains(&model, r, IN_CYLINDER_ONLY), "cylinder-only point");
        assert!(contains(&model, r, IN_BOTH), "shared point");
    }

    #[test]
    fn intersection_contains_only_the_overlap() {
        let (model, r) = run(BooleanOp::Intersection);
        assert!(
            !contains(&model, r, IN_CUBE_ONLY),
            "cube-only point must be out"
        );
        assert!(
            !contains(&model, r, IN_CYLINDER_ONLY),
            "cylinder-only point must be out"
        );
        assert!(contains(&model, r, IN_BOTH), "shared point must be in");
    }

    #[test]
    fn difference_removes_the_second_operand() {
        let (model, r) = run(BooleanOp::Difference);
        assert!(
            contains(&model, r, IN_CUBE_ONLY),
            "cube-only point must remain"
        );
        assert!(
            !contains(&model, r, IN_CYLINDER_ONLY),
            "cylinder-only point must be out"
        );
        assert!(
            !contains(&model, r, IN_BOTH),
            "the drilled-out region must be gone"
        );
    }

    #[test]
    fn union_of_box_and_cylinder_is_valid() {
        check_op("box_cylinder_drilled_hole_through", BooleanOp::Union);
    }

    #[test]
    fn intersection_of_box_and_cylinder_is_valid() {
        check_op("box_cylinder_drilled_hole_through", BooleanOp::Intersection);
    }

    #[test]
    fn difference_of_box_and_cylinder_is_valid() {
        check_op("box_cylinder_drilled_hole_through", BooleanOp::Difference);
    }

    /// A bore through a block leaves rims where its side meets the block's
    /// top and bottom, whether the bore's caps lie in those faces (the rims
    /// are the caps' own edges, imprinted) or stand out of them (the rims
    /// are traced where the side crosses the faces, and fitted like any
    /// other traced curve). Either way four rim edges lie on each face, and
    /// the volume is the block's less the bore's.
    #[test]
    fn rims_of_a_bore_whether_flush_or_standing_out() {
        use geop_core_topology::Body;
        let radius = 0.3;
        for (base, height) in [(0.0, 1.0), (-0.5, 2.0)] {
            let mut model = M::new();
            let block = cube(&mut model, [0., 0., 0.], [1., 1., 1.]);
            let bore = cylinder(&mut model, [0.5, 0.5, base], radius, height, Axis::Z);
            let result = op(&mut model, block, bore, BooleanOp::Difference);
            let topology = model.topology();
            let curved = topology
                .iter_body_edges(Body::Solid(result))
                .unwrap()
                .filter(|&e| {
                    topology
                        .get_edge(e)
                        .unwrap()
                        .curve
                        .as_line()
                        .unwrap()
                        .is_none()
                })
                .count();
            // Four quarters on each rim.
            assert_eq!(curved, 8, "base {base}");
            let r = Scalar::from_f64(radius);
            let expected = Scalar::ONE.sub(Scalar::PI.mul(r).mul(r));
            let volume = topology.mass_properties(result).unwrap().volume;
            assert!(volume.could_be_equal(expected), "base {base}: {volume:?}");
        }
    }

    // ── Chained booleans on the basic shapes, as the CAD front end builds them ──
    //
    // Each solid is created just before the boolean that uses it: every step
    // validates the *whole* model, and a solid not yet combined with anything
    // legitimately overlaps the others in space — which the validation would
    // report as edges crossing faces.

    type M = geop_ops_model::NamedModel;

    /// `CreateCube(offset, dims)`: the axis-aligned box from `offset` to
    /// `offset + dims`.
    fn cube(model: &mut M, offset: [f64; 3], dims: [f64; 3]) -> geop_core_topology::SolidId {
        let max = [0, 1, 2].map(|k| offset[k] + dims[k]);
        block(model, &fresh_id(), offset, max)
    }

    /// `CreateSphere(offset, r)`.
    fn sphere(model: &mut M, center: [f64; 3], r: f64) -> geop_core_topology::SolidId {
        let [x, y, z] = center;
        let r = Scalar::from_f64(r);
        geop_lib_sweep::shapes::sphere_solid(model, &fresh_id(), v3(x, y, z), r).unwrap()
    }

    /// `CreateCylinder(offset, r, h, axis)`: `offset` is the bottom cap's centre.
    fn cylinder(
        model: &mut M,
        base: [f64; 3],
        r: f64,
        h: f64,
        axis: Axis,
    ) -> geop_core_topology::SolidId {
        let [x, y, z] = base;
        geop_lib_sweep::shapes::revolved_cylinder_along_axis(
            model,
            &fresh_id(),
            v3(x, y, z),
            Scalar::from_f64(r),
            Scalar::from_f64(h),
            axis,
        )
        .unwrap()
    }

    /// Run one boolean, requiring it to succeed with a non-empty solid that
    /// passes `validate_fast` and the full `validate`.
    fn op(
        model: &mut M,
        a: geop_core_topology::SolidId,
        b: geop_core_topology::SolidId,
        op: BooleanOp,
    ) -> geop_core_topology::SolidId {
        let result = boolean(model, &namer(), a, b, op, RemeshParams::default())
            .unwrap_or_else(|e| panic!("{op:?} failed: {e:?}"))
            .unwrap_or_else(|| panic!("{op:?} produced an empty solid"));
        model.check_names().unwrap();
        let topology = model.topology();
        if let Err(errors) = validate_fast(&validation(), topology) {
            panic!(
                "{op:?}: {} validate_fast error(s): {}",
                errors.len(),
                errors[0]
            );
        }
        // The full validation too, as `assert_valid` holds a solid to it:
        // `validate_fast` checks each face on its own, and passes an edge
        // running through the middle of another face, or two edges crossing
        // away from any vertex. A result touching itself along a line,
        // its edge shared by four faces, is a valid solid.
        if let Err(errors) =
            geop_core_topology::validation::validate_manifold(&validation(), topology)
        {
            panic!("{op:?}: {} validate error(s): {}", errors.len(), errors[0]);
        }
        // Closed, which neither check covers: every edge of the result is used
        // by an even number of coedges — two where two faces meet, four where
        // the solid touches itself along the edge. An odd count is a hole.
        for face in topology.solid_faces(result).unwrap() {
            for coedge in topology.iterate_face_coedges(face) {
                // A bare-vertex boundary has no edge to share.
                let Ok(edge) = topology.get_coedge(coedge).unwrap().edge() else {
                    continue;
                };
                let uses = topology.coedges_of_edge(edge).len();
                if !uses.is_multiple_of(2) {
                    let e = topology.get_edge(edge).unwrap();
                    let at = |v| {
                        let p = topology.get_vertex(v).unwrap().point;
                        [p[0].to_f64(), p[1].to_f64(), p[2].to_f64()]
                    };
                    panic!(
                        "{op:?}: edge {edge} of face {face}, from {:?} to {:?}, is used by {uses} coedge(s) — the solid is open there",
                        at(e.start_vertex),
                        at(e.end_vertex)
                    );
                }
            }
        }
        result
    }

    /// Reported from the web front end: a cube with a bore through it, the
    /// sphere inscribed in the cube added back, then the `y > 0` half cut
    /// away. The last difference failed in `remesh_edges_x_edges` with "could
    /// not locate vertex on coedge's pcurve" at the cube corner
    /// `(-0.5, 0, -0.5)`, where the cutting cube's face meets an edge of the
    /// first result.
    #[test]
    #[ignore = "slow: 6 s — run with `cargo test -- --ignored`"]
    fn bored_cube_plus_inscribed_sphere_minus_half() {
        let mut model = M::new();
        let block = cube(&mut model, [-0.5, -0.5, -0.5], [1.0, 1.0, 1.0]);
        let bore = cylinder(&mut model, [0.0, 0.0, -0.875], 0.5, 1.75, Axis::Z);
        let bored = op(&mut model, block, bore, BooleanOp::Difference);
        let ball = sphere(&mut model, [0.0, 0.0, 0.0], 0.5);
        let filled = op(&mut model, bored, ball, BooleanOp::Union);
        let half = cube(&mut model, [-0.5, 0.0, -0.5], [1.0, 1.0, 1.0]);
        op(&mut model, filled, half, BooleanOp::Difference);
    }

    /// Reported from the web front end: a cube with a flush inscribed bore
    /// (tangent to the four sides, caps coplanar with top and bottom), the
    /// inscribed sphere added back, then everything with `x < 0.05` cut away
    /// like a section view. The result looked open — two faces missing —
    /// with spurious triangles in the cutting plane. Not a manifold even when
    /// right: the bore touches the cube's sides along lines, and the one at
    /// `x = 0.5` survives the cut.
    #[test]
    #[ignore = "slow: 5 s — run with `cargo test -- --ignored`"]
    fn flush_bored_cube_plus_inscribed_sphere_section() {
        let mut model = M::new();
        let block = cube(&mut model, CORNER, UNIT);
        let bore = cylinder(&mut model, [0.0, 0.0, -0.5], 0.5, 1.0, Axis::Z);
        let bored = op(&mut model, block, bore, BooleanOp::Difference);
        let ball = sphere(&mut model, [0.0, 0.0, 0.0], 0.5);
        let filled = op(&mut model, bored, ball, BooleanOp::Union);
        let cutter = cube(&mut model, [-1.45, -1.05, -0.82], [1.5, 2.1, 1.65]);
        op(&mut model, filled, cutter, BooleanOp::Difference);
    }

    /// Like [`op`], for a boolean whose correct result is empty.
    fn op_empty(
        model: &mut M,
        a: geop_core_topology::SolidId,
        b: geop_core_topology::SolidId,
        op: BooleanOp,
    ) {
        let result = boolean(model, &namer(), a, b, op, RemeshParams::default())
            .unwrap_or_else(|e| panic!("{op:?} failed: {e:?}"));
        assert!(result.is_none(), "{op:?} should be empty");
    }

    const UNIT: [f64; 3] = [1.0, 1.0, 1.0];
    const CORNER: [f64; 3] = [-0.5, -0.5, -0.5];

    // Tangent contact. Where two surfaces touch without crossing, the
    // intersection is a curve or point along which a tiny error in one
    // surface moves the contact a long way (a perpendicular error `d` shifts a
    // tangency by `~sqrt(2 R d)`), which is what broke
    // `bored_cube_plus_inscribed_sphere_minus_half`.

    /// The cylinder touches all four side faces along vertical lines.
    #[test]
    fn cube_minus_inscribed_cylinder() {
        let mut model = M::new();
        let block = cube(&mut model, CORNER, UNIT);
        let bore = cylinder(&mut model, [0.0, 0.0, -0.875], 0.5, 1.75, Axis::Z);
        op(&mut model, block, bore, BooleanOp::Difference);
    }

    /// The cylinder is tangent to the four sides *and* its caps are flush
    /// with the top and bottom: each cap's rim lies in the cube's face and
    /// touches that face's four edges at their midpoints. The first step of
    /// `flush_bored_cube_plus_inscribed_sphere_section`.
    #[test]
    fn cube_minus_flush_inscribed_cylinder() {
        let mut model = M::new();
        let block = cube(&mut model, CORNER, UNIT);
        let bore = cylinder(&mut model, [0.0, 0.0, -0.5], 0.5, 1.0, Axis::Z);
        op(&mut model, block, bore, BooleanOp::Difference);
    }

    /// The sphere touches each face of the cube at its centre.
    #[test]
    fn cube_minus_inscribed_sphere() {
        let mut model = M::new();
        let block = cube(&mut model, CORNER, UNIT);
        let ball = sphere(&mut model, [0.0, 0.0, 0.0], 0.5);
        op(&mut model, block, ball, BooleanOp::Difference);
    }

    /// The sphere passes through all eight corners of the cube.
    #[test]
    fn cube_intersect_circumscribed_sphere() {
        let mut model = M::new();
        let block = cube(&mut model, CORNER, UNIT);
        let ball = sphere(&mut model, [0.0, 0.0, 0.0], 3f64.sqrt() / 2.0);
        op(&mut model, block, ball, BooleanOp::Intersection);
    }

    /// The cylinder's surface contains the cube's four vertical edges.
    #[test]
    fn cube_intersect_cylinder_through_its_edges() {
        let mut model = M::new();
        let block = cube(&mut model, CORNER, UNIT);
        let tube = cylinder(&mut model, [0.0, 0.0, -1.0], 0.5f64.sqrt(), 2.0, Axis::Z);
        op(&mut model, block, tube, BooleanOp::Intersection);
    }

    /// A coaxial cylinder of the sphere's radius touches it along the
    /// equator, a whole circle of tangency.
    #[test]
    fn sphere_union_tangent_cylinder() {
        let mut model = M::new();
        let ball = sphere(&mut model, [0.0, 0.0, 0.0], 0.5);
        let tube = cylinder(&mut model, [0.0, 0.0, -1.0], 0.5, 2.0, Axis::Z);
        op(&mut model, ball, tube, BooleanOp::Union);
    }

    /// A cube resting on the sphere's pole, where the sphere's patches all
    /// collapse to one point: tangent contact at a singular point.
    #[test]
    fn sphere_union_cube_touching_its_pole() {
        let mut model = M::new();
        let ball = sphere(&mut model, [0.0, 0.0, 0.0], 0.5);
        let block = cube(&mut model, [-0.5, -0.5, 0.5], UNIT);
        op(&mut model, ball, block, BooleanOp::Union);
    }

    /// A cylinder along x lying on the cube's top face: tangent along a line.
    #[test]
    fn cube_union_cylinder_lying_on_top() {
        let mut model = M::new();
        let block = cube(&mut model, CORNER, UNIT);
        let log = cylinder(&mut model, [-1.0, 0.0, 0.75], 0.25, 2.0, Axis::X);
        op(&mut model, block, log, BooleanOp::Union);
    }

    // Coincident features: shared faces, edges and corners, coplanar caps, and
    // identical solids.

    #[test]
    fn cubes_sharing_a_face_union() {
        let mut model = M::new();
        let a = cube(&mut model, CORNER, UNIT);
        let b = cube(&mut model, [0.5, -0.5, -0.5], UNIT);
        op(&mut model, a, b, BooleanOp::Union);
    }

    #[test]
    fn cubes_sharing_an_edge_union() {
        let mut model = M::new();
        let a = cube(&mut model, CORNER, UNIT);
        let b = cube(&mut model, [0.5, 0.5, -0.5], UNIT);
        op(&mut model, a, b, BooleanOp::Union);
    }

    #[test]
    fn cubes_sharing_a_corner_union() {
        let mut model = M::new();
        let a = cube(&mut model, CORNER, UNIT);
        let b = cube(&mut model, [0.5, 0.5, 0.5], UNIT);
        op(&mut model, a, b, BooleanOp::Union);
    }

    /// A half-overlapping cube: two of its faces are coplanar with the first
    /// cube's, over part of their area.
    #[test]
    fn cube_minus_offset_cube_with_coplanar_faces() {
        let mut model = M::new();
        let a = cube(&mut model, CORNER, UNIT);
        let b = cube(&mut model, [0.0, -0.5, -0.5], UNIT);
        op(&mut model, a, b, BooleanOp::Difference);
    }

    /// A hole whose caps are flush with the cube's top and bottom.
    #[test]
    fn cube_minus_flush_cylinder() {
        let mut model = M::new();
        let block = cube(&mut model, CORNER, UNIT);
        let bore = cylinder(&mut model, [0.0, 0.0, -0.5], 0.3, 1.0, Axis::Z);
        op(&mut model, block, bore, BooleanOp::Difference);
    }

    /// A plate 6 by 3 and a blind hole to drill into it from its top: a
    /// cylinder of radius 0.33 around `(3, 1.5)`, its start cap flush with
    /// the top.
    fn plate_and_blind_hole(
        model: &mut M,
    ) -> (geop_core_topology::SolidId, geop_core_topology::SolidId) {
        let plate = cube(model, [0.0, 0.0, 0.0], [6.0, 3.0, 0.5]);
        let hole = cylinder(model, [3.0, 1.5, 0.2], 0.33, 0.3, Axis::Z);
        (plate, hole)
    }

    /// The plate and hole of [`plate_and_blind_hole`] remeshed, and the disc
    /// the remesh cuts from the plate's top, inside the hole's rim.
    fn plate_and_blind_hole_remeshed(
        model: &mut M,
    ) -> (geop_core_topology::SolidId, geop_core_topology::FaceId) {
        let (plate, hole) = plate_and_blind_hole(model);
        crate::remesh::remesh(model, &namer(), plate, hole, RemeshParams::default()).unwrap();
        let topology = model.topology();
        if let Err(errors) = validate_fast(&validation(), topology) {
            panic!("{} validate_fast error(s): {}", errors.len(), errors[0]);
        }
        let in_rim = |p: &[f64; 3]| (p[0] - 3.0).hypot(p[1] - 1.5) < 0.34 && p[2] > 0.49;
        let disc = topology
            .solid_faces(plate)
            .unwrap()
            .into_iter()
            .find(|&face| {
                let corners = topology.face_corners(face).unwrap();
                corners.len() == 1 && corners[0].iter().all(|(_, p)| in_rim(p))
            })
            .expect("a face of the plate inside the hole's rim");
        (hole, disc)
    }

    /// Reported as the parametric plate made 6 wide with an M6 hole: the
    /// result was left open along the hole's rim, as the disc cut from the
    /// top was classified outside the hole rather than on its cap.
    #[test]
    fn plate_minus_blind_hole_from_its_top() {
        let mut model = M::new();
        let (plate, hole) = plate_and_blind_hole(&mut model);
        op(&mut model, plate, hole, BooleanOp::Difference);
    }

    /// A bolt hole patterned round a revolved flange: the sixth copy's centre
    /// landed `3.6e-4` beside the flange cap's quarter seam, the angle of the
    /// seed given to a thousandth of a degree. The hole's cap is flush with
    /// the flange's, and its own quarter seam, a radial line from its centre
    /// pole, crosses the flange's seam `3.6e-4` from that pole: a genuine
    /// transversal crossing of two straight edges, which the remesh refused
    /// as a near-duplicate of the pole and panicked on. Here on a plain
    /// cylinder: the tool's `+x` seam from `(-3.6e-4, 28)` crosses the
    /// cylinder's `+y` seam at `(0, 28)`, and the piece of the tool's seam
    /// between the two is `3.6e-4` long.
    #[test]
    fn blind_hole_whose_cap_seam_crosses_the_caps_seam_beside_its_centre() {
        let mut model = M::new();
        let disc = cylinder(&mut model, [0.0, 0.0, 0.0], 40.0, 6.0, Axis::Z);
        let hole = cylinder(
            &mut model,
            [-0.000_355_652_98, 27.999_383_99, 3.0],
            3.3,
            3.0,
            Axis::Z,
        );
        op(&mut model, disc, hole, BooleanOp::Difference);
    }

    /// The disc cut from the plate's top lies on the hole's start cap, so
    /// every point a boolean classifies it by is on the hole's boundary.
    #[test]
    fn blind_hole_cap_covers_the_disc_cut_from_the_top() {
        use geop_core_topology::contains::{
            face::face_interior_point_where,
            shell::{PointClassification, shell_contains},
        };
        let mut model = M::new();
        let (hole, disc) = plate_and_blind_hole_remeshed(&mut model);
        let topology = model.topology();
        let params = RemeshParams::default();
        let shell = topology.get_solid(hole).unwrap().shells[0];
        let mut checked = 0;
        face_interior_point_where(
            topology,
            disc,
            params.max_nodes,
            params.curve_curve_min_subdivision_size,
            super::SEED,
            |u, v| {
                let point = topology.get_face(disc)?.surface.evaluate(u, v)?;
                let class = shell_contains(
                    topology,
                    shell,
                    point,
                    params.max_nodes,
                    params.curve_curve_min_subdivision_size,
                    super::SEED,
                )?;
                assert!(
                    !matches!(
                        class,
                        PointClassification::Inside | PointClassification::Outside
                    ),
                    "{point:?} of the disc, at uv=({u:?}, {v:?}), is {class:?}"
                );
                checked += 1;
                Ok(false)
            },
        )
        .unwrap();
        assert!(checked > 0);
    }

    /// The point the disc cut from the plate's top was once classified by:
    /// beside it — 0.44 from the hole's axis, the rim 0.33 — and outside it,
    /// as a point and as the box around it `face_interior_point_where` asks
    /// about. With the boolean's seed, the first ray from the box grazes the
    /// rim; cast from the box itself, it counted that as one crossing.
    #[test]
    fn point_beside_the_disc_cut_from_the_top_is_outside_it() {
        use geop_core_topology::contains::face::{PointClassification, face_contains};
        let mut model = M::new();
        let (_, disc) = plate_and_blind_hole_remeshed(&mut model);
        let topology = model.topology();
        let (u, v) = (
            Scalar::from_f64(0.42698015887332674),
            Scalar::from_f64(0.5220523370528238),
        );
        let point = topology
            .get_face(disc)
            .unwrap()
            .surface
            .evaluate(u, v)
            .unwrap();
        assert!(
            (point[0].to_f64() - 3.0).hypot(point[1].to_f64() - 1.5) > 0.4,
            "{point:?}"
        );
        let params = RemeshParams::default();
        let epsilon = params.curve_curve_min_subdivision_size;
        let neighbourhood = |t: Scalar| t.sub(epsilon).union(t.add(epsilon));
        let mut wrong = Vec::new();
        for (u, v) in [(u, v), (neighbourhood(u), neighbourhood(v))] {
            for seed in (0..64).chain([super::SEED]) {
                let class =
                    face_contains(topology, disc, u, v, params.max_nodes, epsilon, seed).unwrap();
                if class != PointClassification::Outside {
                    wrong.push(format!("{class:?} at ({u:?}, {v:?}) with seed {seed}"));
                }
            }
        }
        assert!(wrong.is_empty(), "{}", wrong.join("\n"));
    }

    #[test]
    fn identical_spheres_union() {
        let mut model = M::new();
        let a = sphere(&mut model, [0.0, 0.0, 0.0], 0.5);
        let b = sphere(&mut model, [0.0, 0.0, 0.0], 0.5);
        op(&mut model, a, b, BooleanOp::Union);
    }

    #[test]
    fn identical_cubes_difference_is_empty() {
        let mut model = M::new();
        let a = cube(&mut model, CORNER, UNIT);
        let b = cube(&mut model, CORNER, UNIT);
        op_empty(&mut model, a, b, BooleanOp::Difference);
    }

    /// Every operator on two copies of one shape — same geometry, separate
    /// topology with ids of their own — so that each face of either copy
    /// coincides with a face of the other over its whole area, and every edge
    /// and vertex lies exactly on one of the other copy's. Union and
    /// intersection must give the shape back, difference must be empty.
    fn check_duplicated_shape(make: impl Fn(&mut M) -> geop_core_topology::SolidId) {
        let probes = [(0.0, 0.0, 0.0), (0.2, -0.1, 0.15), (0.9, 0.0, 0.0)];
        let inside_original = {
            let mut model = M::new();
            let shape = make(&mut model);
            probes.map(|p| contains(model.topology(), shape, p))
        };
        for operator in [BooleanOp::Union, BooleanOp::Intersection] {
            let mut model = M::new();
            let a = make(&mut model);
            let b = make(&mut model);
            let result = op(&mut model, a, b, operator);
            for (p, inside) in probes.iter().zip(inside_original) {
                assert_eq!(
                    contains(model.topology(), result, *p),
                    inside,
                    "{operator:?} of two copies must contain {p:?} exactly when one copy does"
                );
            }
        }
        let mut model = M::new();
        let a = make(&mut model);
        let b = make(&mut model);
        op_empty(&mut model, a, b, BooleanOp::Difference);
    }

    #[test]
    fn duplicated_cube_all_operators() {
        check_duplicated_shape(|model| cube(model, CORNER, UNIT));
    }

    #[test]
    #[ignore = "slow: 6 s — run with `cargo test -- --ignored`"]
    fn duplicated_sphere_all_operators() {
        check_duplicated_shape(|model| sphere(model, [0.0, 0.0, 0.0], 0.5));
    }

    #[test]
    #[ignore = "slow: 8 s — run with `cargo test -- --ignored`"]
    fn duplicated_cylinder_all_operators() {
        check_duplicated_shape(|model| cylinder(model, [0.0, 0.0, -0.5], 0.5, 1.0, Axis::Z));
    }

    /// A sphere centred on the cube's corner: its three coordinate-plane
    /// seams lie exactly in three faces of the cube.
    #[test]
    fn cube_minus_sphere_on_its_corner() {
        let mut model = M::new();
        let block = cube(&mut model, CORNER, UNIT);
        let ball = sphere(&mut model, [0.5, 0.5, 0.5], 0.5);
        op(&mut model, block, ball, BooleanOp::Difference);
    }

    // Curved against curved.

    /// Two equal cylinders crossing at right angles: their intersection
    /// curves (two ellipses) cross each other at two singular points.
    #[test]
    fn steinmetz_cylinders_union() {
        let mut model = M::new();
        let a = cylinder(&mut model, [0.0, 0.0, -1.0], 0.5, 2.0, Axis::Z);
        let b = cylinder(&mut model, [-1.0, 0.0, 0.0], 0.5, 2.0, Axis::X);
        op(&mut model, a, b, BooleanOp::Union);
    }

    #[test]
    fn steinmetz_cylinders_intersection() {
        let mut model = M::new();
        let a = cylinder(&mut model, [0.0, 0.0, -1.0], 0.5, 2.0, Axis::Z);
        let b = cylinder(&mut model, [-1.0, 0.0, 0.0], 0.5, 2.0, Axis::X);
        op(&mut model, a, b, BooleanOp::Intersection);
    }

    /// A bore along the sphere's axis, through both poles.
    #[test]
    fn sphere_minus_cylinder_through_its_poles() {
        let mut model = M::new();
        let ball = sphere(&mut model, [0.0, 0.0, 0.0], 0.5);
        let bore = cylinder(&mut model, [0.0, 0.0, -1.0], 0.2, 2.0, Axis::Z);
        op(&mut model, ball, bore, BooleanOp::Difference);
    }

    /// A sphere centred on a cylinder's cap: the cap cuts it along its
    /// equator, which is also where its patches meet.
    #[test]
    fn cylinder_union_sphere_on_its_cap() {
        let mut model = M::new();
        let tube = cylinder(&mut model, [0.0, 0.0, -1.0], 0.5, 1.0, Axis::Z);
        let ball = sphere(&mut model, [0.0, 0.0, 0.0], 0.3);
        op(&mut model, tube, ball, BooleanOp::Union);
    }

    /// The same with the sphere's radius equal to the cylinder's: the
    /// sphere's equator coincides with the cap's rim, *and* the sphere is
    /// tangent to the cylinder's side along that same circle.
    #[test]
    fn cylinder_union_equal_sphere_on_its_cap() {
        let mut model = M::new();
        let tube = cylinder(&mut model, [0.0, 0.0, -1.0], 0.5, 1.0, Axis::Z);
        let ball = sphere(&mut model, [0.0, 0.0, 0.0], 0.5);
        op(&mut model, tube, ball, BooleanOp::Union);
    }

    /// External tangency at a point of the sphere's equator (not a pole).
    #[test]
    fn cube_union_sphere_touching_a_face() {
        let mut model = M::new();
        let block = cube(&mut model, CORNER, UNIT);
        let ball = sphere(&mut model, [1.0, 0.0, 0.0], 0.5);
        op(&mut model, block, ball, BooleanOp::Union);
    }

    /// A sphere half sunk into the cube through the centre of a face.
    #[test]
    fn cube_union_sphere_on_a_face() {
        let mut model = M::new();
        let block = cube(&mut model, CORNER, UNIT);
        let ball = sphere(&mut model, [0.5, 0.0, 0.0], 0.3);
        op(&mut model, block, ball, BooleanOp::Union);
    }

    /// Two inscribed bores at right angles: the first two steps of
    /// `cube_minus_three_inscribed_bores`. The second bore's wall meets the
    /// first's in two Steinmetz ellipses, which end exactly where both bores
    /// touch the cube's faces — at the midpoints of its edges.
    #[test]
    #[ignore = "slow: 4 s — run with `cargo test -- --ignored`"]
    fn cube_minus_two_inscribed_bores() {
        let mut model = M::new();
        let block = cube(&mut model, CORNER, UNIT);
        let z = cylinder(&mut model, [0.0, 0.0, -0.875], 0.5, 1.75, Axis::Z);
        let a = op(&mut model, block, z, BooleanOp::Difference);
        let x = cylinder(&mut model, [-0.875, 0.0, 0.0], 0.5, 1.75, Axis::X);
        op(&mut model, a, x, BooleanOp::Difference);
    }

    /// Remesh `a` against `b`, then require an edge shared by faces of both
    /// between each pair of points in `arcs`, and the full validation to
    /// pass.
    ///
    /// The full validation cannot see a missing intersection arc that ends
    /// at vertices both solids share: faces sharing a vertex are not
    /// searched for crossings. A boolean only notices one as a face it cannot
    /// classify, or a result that is open, far from the cause. So the arcs a
    /// scene must have are named, and on failure the edges it does have are
    /// listed.
    fn check_remesh_imprints(
        model: &mut M,
        a: geop_core_topology::SolidId,
        b: geop_core_topology::SolidId,
        arcs: &[[[f64; 3]; 2]],
    ) {
        crate::remesh::remesh(model, &namer(), a, b, RemeshParams::default())
            .unwrap_or_else(|e| panic!("remesh failed: {e}"));
        let topology = model.topology();
        let faces_a = topology.solid_faces(a).unwrap();
        let faces_b = topology.solid_faces(b).unwrap();
        let mut shared = Vec::new();
        for (&edge_id, edge) in &topology.edges {
            let faces: Vec<_> = topology
                .coedges_of_edge(edge_id)
                .iter()
                .map(|&c| topology.get_coedge(c).unwrap().face)
                .collect();
            if faces.iter().any(|f| faces_a.contains(f))
                && faces.iter().any(|f| faces_b.contains(f))
            {
                shared.push([edge.start_vertex, edge.end_vertex].map(|v| {
                    let p = topology.get_vertex(v).unwrap().point;
                    [p[0].to_f64(), p[1].to_f64(), p[2].to_f64()]
                }));
            }
        }
        // Telling which known point a vertex is, not a kernel comparison:
        // the expected points are well apart, so any coarse tolerance names
        // them unambiguously.
        let near = |p: [f64; 3], q: [f64; 3]| (0..3).all(|i| (p[i] - q[i]).abs() < 1e-4);
        let missing: Vec<_> = arcs
            .iter()
            .filter(|[p, q]| {
                !shared
                    .iter()
                    .any(|[s, e]| (near(*s, *p) && near(*e, *q)) || (near(*s, *q) && near(*e, *p)))
            })
            .collect();
        assert!(
            missing.is_empty(),
            "no edge shared by both solids between {missing:?}; the shared edges are {shared:?}"
        );

        if let Err(errors) = geop_core_topology::validation::validate(&validation(), topology) {
            let all: Vec<String> = errors.iter().map(|e| format!("{e}")).collect();
            panic!(
                "{} validate error(s):\n{}",
                errors.len(),
                all.join("\n---\n")
            );
        }
    }

    /// The remesh behind `cube_minus_two_inscribed_bores`' second step: each
    /// Steinmetz ellipse `x = ±z` of the two walls runs through midpoints
    /// `(±0.5, 0, ±0.5)` of the cube's edges and through `(0, ±0.5, 0)`,
    /// where the two ellipses cross with the walls tangent — so each of its
    /// eight arcs can only be traced from the midpoint end.
    #[test]
    #[ignore = "slow: 5 s — run with `cargo test -- --ignored`"]
    fn two_inscribed_bores_remesh_imprints_every_arc() {
        let mut model = M::new();
        let block = cube(&mut model, CORNER, UNIT);
        let z = cylinder(&mut model, [0.0, 0.0, -0.875], 0.5, 1.75, Axis::Z);
        let a = op(&mut model, block, z, BooleanOp::Difference);
        let x = cylinder(&mut model, [-0.875, 0.0, 0.0], 0.5, 1.75, Axis::X);
        let mut arcs = Vec::new();
        for (mx, mz) in [(0.5, 0.5), (0.5, -0.5), (-0.5, 0.5), (-0.5, -0.5)] {
            for cy in [0.5, -0.5] {
                arcs.push([[mx, 0.0, mz], [0.0, cy, 0.0]]);
            }
        }
        check_remesh_imprints(&mut model, a, x, &arcs);
    }

    /// The remesh behind `cube_minus_three_inscribed_bores`' third step. In
    /// each octant the three walls meet at one point, `(±1, ±1, ±1) / sqrt 8`,
    /// where the arc the third wall cuts from the first meets the arc it cuts
    /// from the second: one runs from there to `(0, ±0.5, ±0.5)`, the other
    /// to `(±0.5, ±0.5, 0)`. That point is interior to the third wall, so both
    /// arcs leave it into the same face.
    #[test]
    #[ignore = "slow: 9 s — run with `cargo test -- --ignored`"]
    fn three_inscribed_bores_remesh_imprints_every_arc() {
        let mut model = M::new();
        let block = cube(&mut model, CORNER, UNIT);
        let z = cylinder(&mut model, [0.0, 0.0, -0.875], 0.5, 1.75, Axis::Z);
        let a = op(&mut model, block, z, BooleanOp::Difference);
        let x = cylinder(&mut model, [-0.875, 0.0, 0.0], 0.5, 1.75, Axis::X);
        let b = op(&mut model, a, x, BooleanOp::Difference);
        let y = cylinder(&mut model, [0.0, -0.875, 0.0], 0.5, 1.75, Axis::Y);
        let t = 0.125f64.sqrt();
        let mut arcs = Vec::new();
        for sx in [1.0, -1.0] {
            for sy in [1.0, -1.0] {
                for sz in [1.0, -1.0] {
                    let triple = [sx * t, sy * t, sz * t];
                    arcs.push([triple, [0.0, sy * 0.5, sz * 0.5]]);
                    arcs.push([triple, [sx * 0.5, sy * 0.5, 0.0]]);
                }
            }
        }
        check_remesh_imprints(&mut model, b, y, &arcs);
    }

    /// Inscribed bores along all three axes, one after another (a "jack").
    /// Each bore is tangent to four faces, and each later bore crosses the
    /// earlier ones at Steinmetz points.
    #[test]
    #[ignore = "slow: 8 s — run with `cargo test -- --ignored`"]
    fn cube_minus_three_inscribed_bores() {
        let mut model = M::new();
        let block = cube(&mut model, CORNER, UNIT);
        let z = cylinder(&mut model, [0.0, 0.0, -0.875], 0.5, 1.75, Axis::Z);
        let a = op(&mut model, block, z, BooleanOp::Difference);
        let x = cylinder(&mut model, [-0.875, 0.0, 0.0], 0.5, 1.75, Axis::X);
        let b = op(&mut model, a, x, BooleanOp::Difference);
        let y = cylinder(&mut model, [0.0, -0.875, 0.0], 0.5, 1.75, Axis::Y);
        op(&mut model, b, y, BooleanOp::Difference);
    }

    /// Thinner bores, so the cube keeps its faces: two crossing bores.
    #[test]
    fn cube_minus_two_crossing_bores() {
        let mut model = M::new();
        let block = cube(&mut model, CORNER, UNIT);
        let z = cylinder(&mut model, [0.0, 0.0, -0.875], 0.3, 1.75, Axis::Z);
        let a = op(&mut model, block, z, BooleanOp::Difference);
        let x = cylinder(&mut model, [-0.875, 0.0, 0.0], 0.3, 1.75, Axis::X);
        op(&mut model, a, x, BooleanOp::Difference);
    }

    // More tangencies, and a longer chain.

    /// Two spheres touching externally at a single point.
    #[test]
    fn touching_spheres_union() {
        let mut model = M::new();
        let a = sphere(&mut model, [0.0, 0.0, 0.0], 0.5);
        let b = sphere(&mut model, [1.0, 0.0, 0.0], 0.5);
        op(&mut model, a, b, BooleanOp::Union);
    }

    /// A sphere inside a cylinder of the same radius, touching it along the
    /// equator from inside.
    #[test]
    fn cylinder_minus_inscribed_sphere() {
        let mut model = M::new();
        let tube = cylinder(&mut model, [0.0, 0.0, -1.0], 0.5, 2.0, Axis::Z);
        let ball = sphere(&mut model, [0.0, 0.0, 0.0], 0.5);
        op(&mut model, tube, ball, BooleanOp::Difference);
    }

    /// Two parallel cylinders touching along a line.
    #[test]
    fn touching_parallel_cylinders_union() {
        let mut model = M::new();
        let a = cylinder(&mut model, [0.0, 0.0, -0.5], 0.5, 1.0, Axis::Z);
        let b = cylinder(&mut model, [1.0, 0.0, -0.5], 0.5, 1.0, Axis::Z);
        op(&mut model, a, b, BooleanOp::Union);
    }

    /// A sphere touching one of the cube's edges from outside.
    #[test]
    fn cube_union_sphere_touching_an_edge() {
        let mut model = M::new();
        let block = cube(&mut model, CORNER, UNIT);
        let d = 0.5 + 0.5 / 2f64.sqrt();
        let ball = sphere(&mut model, [d, d, 0.0], 0.5);
        op(&mut model, block, ball, BooleanOp::Union);
    }

    /// A small part as a user might build it: a plate with a boss and a
    /// through hole, then a rounded cap on the boss.
    #[test]
    fn plate_with_boss_hole_and_cap() {
        let mut model = M::new();
        let plate = cube(&mut model, [-1.0, -1.0, -0.25], [2.0, 2.0, 0.5]);
        let boss = cylinder(&mut model, [0.0, 0.0, 0.25], 0.5, 0.5, Axis::Z);
        let a = op(&mut model, plate, boss, BooleanOp::Union);
        let cap = sphere(&mut model, [0.0, 0.0, 0.75], 0.5);
        let b = op(&mut model, a, cap, BooleanOp::Union);
        let hole = cylinder(&mut model, [0.0, 0.0, -0.5], 0.25, 1.5, Axis::Z);
        op(&mut model, b, hole, BooleanOp::Difference);
    }

    /// The second step of `chained_differences_block_with_two_slots_and_a_sphere`,
    /// stopped after remesh: a point of the second slot's side face, above the
    /// block, must classify as outside the slotted block. It came back
    /// `Inside` — a ray-parity miscount that the boolean only stopped hiding
    /// once `classify_face` checked every interior point it sampled.
    #[test]
    fn chained_slots_point_above_the_block_is_outside() {
        let mut model = M::new();
        let a = cube(&mut model, [-0.50, -0.50, -0.50], [1.00, 1.00, 1.00]);
        let b = cube(&mut model, [-1.13, -0.30, -0.33], [2.25, 0.60, 0.65]);
        let c = op(&mut model, a, b, BooleanOp::Difference);
        let d = cube(&mut model, [-0.33, -0.28, -1.15], [0.65, 0.55, 2.30]);
        crate::remesh::remesh(&mut model, &namer(), c, d, RemeshParams::default()).unwrap();

        let params = RemeshParams::default();
        let shell = model.topology().get_solid(c).unwrap().shells[0];
        let classification = geop_core_topology::contains::shell::shell_contains(
            model.topology(),
            shell,
            v3(-0.005, -0.28, 0.825),
            params.max_nodes,
            params.curve_curve_min_subdivision_size,
            super::SEED,
        )
        .unwrap();
        assert_eq!(
            classification,
            geop_core_topology::contains::shell::PointClassification::Outside
        );
    }

    #[test]
    fn chained_differences_block_with_two_slots_and_a_sphere() {
        let f = Scalar::from_f64;
        let mut model = M::new();
        let a = cube(&mut model, [-0.50, -0.50, -0.50], [1.00, 1.00, 1.00]);
        let b = cube(&mut model, [-1.13, -0.30, -0.33], [2.25, 0.60, 0.65]);
        let params = RemeshParams::default();
        let c = boolean(&mut model, &namer(), a, b, BooleanOp::Difference, params)
            .unwrap()
            .expect("block minus the first slot must be non-empty");

        let d = cube(&mut model, [-0.33, -0.28, -1.15], [0.65, 0.55, 2.30]);
        let e = boolean(&mut model, &namer(), c, d, BooleanOp::Difference, params)
            .unwrap()
            .expect("minus the second slot must be non-empty");

        let sphere =
            geop_lib_sweep::shapes::sphere_solid(&mut model, "s", Vector3::zero(), f(0.45))
                .unwrap();
        let result = boolean(
            &mut model,
            &namer(),
            e,
            sphere,
            BooleanOp::Difference,
            params,
        )
        .unwrap()
        .expect("minus the sphere must be non-empty");
        model.check_names().unwrap();

        let topology = model.topology();
        assert!(!topology.solid_faces(result).unwrap().is_empty());
        if let Err(errors) = validate_fast(&validation(), topology) {
            panic!("{} validate_fast error(s): {}", errors.len(), errors[0]);
        }

        let scene = geop_core_rasterize::rasterize(topology, 8)
            .unwrap()
            .scene(|_| geop_core_rasterize::debug::Color10::Blue);
        std::fs::create_dir_all("outputs").unwrap();
        scene
            .save_to_file("outputs/chained_differences_block_with_two_slots_and_a_sphere.html")
            .unwrap();
    }

    /// Regression test, captured from a browser session: two cubes
    /// differenced, then a sphere differenced out of that, then a thin slab
    /// cube differenced out of *that*. The second cube and the sphere were
    /// anchored to corners of the first cube, picked by kernel id
    /// (`VertexId(13)` and `VertexId(9)`, its corners `p3,start` and
    /// `p2,start`); their coordinates as the session read them are given
    /// here, to the last bit.
    ///
    /// Used to fail inside the final boolean's `classify_face`:
    /// `shell_contains` reported a point on the slab's boundary that no face
    /// of that shell's `face_contains` agreed contained. Root cause:
    /// `shell_contains`'s face coincidence pre-check only asked
    /// `surface_could_contain` — proximity to a face's *untrimmed* surface —
    /// without checking the point fell within the face's trim. Proximity is
    /// not membership (see `AGENTS.md`).
    #[test]
    fn chained_differences_with_anchored_shapes_and_thin_slab_cutter_succeeds() {
        let mut model = M::new();
        let block = cube(&mut model, [-0.5, -0.5, -0.5], [1.0, 1.0, 1.0]);
        let [x, y, z] = [-0.5000000000000001, 0.5000000000000002, 0.5000000000000001];
        let second = cube(&mut model, [x - 0.5, y - 0.5, z - 0.5], [1.0, 1.0, 1.0]);
        let center = [0.5000000000000003, 0.5000000000000002, 0.5000000000000001];
        let ball = sphere(&mut model, center, 0.5);
        // Only success is asserted, as when this was captured: the full
        // `validate` that `op` runs finds an edge of the second result
        // crossing a face unsplit, which this regression never covered.
        let difference = |model: &mut M, a, b| {
            boolean(
                model,
                &namer(),
                a,
                b,
                BooleanOp::Difference,
                RemeshParams::default(),
            )
            .unwrap()
            .expect("the result is not empty")
        };
        let cut = difference(&mut model, block, second);
        let cut = difference(&mut model, cut, ball);
        let slab = cube(&mut model, [-1.13, -0.10, -0.15], [2.25, 0.20, 0.30]);
        difference(&mut model, cut, slab);
        model.check_names().unwrap();
    }

    /// Regression test: a cylinder drilled through a cube, its height
    /// chosen so *both* caps sit exactly flush with the cube's own
    /// opposite faces (`cylinder z in [-0.5, 0.5]` exactly matching the
    /// cube's), rather than the usual "drilled hole" test scenes (e.g.
    /// `box_cylinder_drilled_hole_through`) where the cylinder
    /// deliberately extends *past* both faces.
    ///
    /// That exact-coplanar-cap case used to (1) take ~30s — over 1000x the
    /// sub-second time every other boolean test in this file takes — and
    /// (2) produce a *wrong* result, cutting the hole through only one of
    /// the two coplanar-cap faces and leaving the other fully solid.
    ///
    /// Root cause: `revolve_at_oriented` builds any flat cap (a
    /// doubly-curved pole, like a sphere's, genuinely needs a fan of
    /// wedges meeting at a center vertex — a flat cap doesn't, but got the
    /// same treatment) as 4 wedges with radial "spoke" edges from rim to
    /// center. `find_coincident_pair` (`remesh_edges_x_faces.rs`) couldn't
    /// tell those spokes apart from the genuine circular rim boundary, so
    /// it tried imprinting all of them — the slowdown, and (via
    /// interleaved imprints across the two caps corrupting
    /// `Model::splice_edge_into_face`'s face-splitting bookkeeping) the
    /// wrong result. Fixed by `seam_plane`/`faces_are_coplanar`
    /// there: an edge whose two neighboring faces (within its own solid)
    /// are already coplanar with each other is never the only carrier of a
    /// genuine coincidence, so it's skipped — cutting straight to the rim
    /// arcs.
    #[test]
    fn cube_minus_z_cylinder_with_coplanar_cap_is_fast_and_correct() {
        let mut model = M::new();
        let a = cube(&mut model, [-0.5, -0.5, -0.5], [1.0, 1.0, 1.0]);
        let b = cylinder(&mut model, [0.0, 0.0, -0.5], 0.3, 1.0, Axis::Z);
        let params = RemeshParams::default();

        let t0 = std::time::Instant::now();
        let result = boolean(&mut model, &namer(), a, b, BooleanOp::Difference, params)
            .unwrap()
            .unwrap();
        let elapsed = t0.elapsed();
        let topology = model.topology();
        // Ballpark-matches `difference_of_box_and_cylinder_is_valid`'s own
        // (non-coincident-cap) ~2s in isolation; running inside the full
        // suite under CPU contention from other parallel tests has been
        // observed up to ~12s. The threshold stays generous — the point is
        // catching a regression back toward the old ~32s (a *further*
        // 1000x-ish blowup on top of ordinary parallel-run noise), not
        // pinning down exact timing.
        assert!(
            elapsed.as_secs() < 20,
            "boolean took {elapsed:?}, expected well under 20s"
        );

        // Every point on both cap planes' hole boundary must be *outside*
        // the result (the cylinder's full bore is open at both ends) —
        // catches "hole only cut on one end" directly, unlike
        // `validate_fast`.
        let shell = topology.get_solid(result).unwrap().shells[0];
        for z in [-0.45, 0.45] {
            let p = v3(0.0, 0.0, z);
            let outside = matches!(
                geop_core_topology::contains::shell::shell_contains(
                    topology,
                    shell,
                    p,
                    params.max_nodes,
                    params.curve_curve_min_subdivision_size,
                    0xC7D1
                )
                .unwrap(),
                geop_core_topology::contains::shell::PointClassification::Outside
            );
            assert!(
                outside,
                "point {p:?} (near a cap's bore) should be outside the drilled result"
            );
        }

        let scene = geop_core_rasterize::rasterize(topology, 8)
            .unwrap()
            .scene(|_| geop_core_rasterize::debug::Color10::Blue);
        std::fs::create_dir_all("outputs").unwrap();
        scene
            .save_to_file("outputs/cube_minus_coplanar_cap_cylinder.html")
            .unwrap();
    }

    /// A bore through a block must actually be empty in the rendered mesh:
    /// no triangle may cover the hole. End-to-end cover for the trimming
    /// that `geop_core_rasterize::clip` does per grid cell — where a
    /// duplicated vertex in a clipped hole outline, or a concave fragment
    /// classified as a whole, used to leave a flap of surface hanging
    /// across an opening (isolated in that module's own tests).
    #[test]
    fn bore_renders_as_a_hole() {
        let mut model = M::new();
        let block = cube(&mut model, [-1.0, -1.0, 0.0], [2.0, 2.0, 0.5]);
        // Radius 0.5 of the 2x2 footprint puts the bore's outline exactly
        // through grid cell corners at this resolution, which is what it
        // takes to produce the duplicated vertex.
        let bore = cylinder(&mut model, [0.0, 0.0, -0.5], 0.5, 1.5, Axis::Z);
        op(&mut model, block, bore, BooleanOp::Difference);

        let rasterized = geop_core_rasterize::rasterize(model.topology(), 24).unwrap();
        for triangle in rasterized.faces.values().flatten() {
            for p in [triangle.a, triangle.b, triangle.c] {
                let (x, y, z) = (p[0].to_f64(), p[1].to_f64(), p[2].to_f64());
                let r = x.hypot(y);
                // Strictly inside the bore, and within the block's height:
                // nothing may be drawn there. The margin keeps the bore
                // wall's own triangles (at r = 0.6, faceted slightly inwards)
                // out of the test.
                assert!(
                    r > 0.45 || !(0.01..0.49).contains(&z),
                    "a triangle corner sits inside the bore at ({x}, {y}, {z})"
                );
            }
        }
    }

    /// A box cut into the top of a cube, its corners known only to within
    /// 1e-12 — as a tool computed from a solid's own (interval) vertices
    /// is: a lip's groove offset from a shelled enclosure's rim. The same
    /// box with sharp corners was always cut cleanly.
    ///
    /// The march fed each corrector step the previous step's honest
    /// `(u, v)` as its first seed, unsharpened, so the width compounded
    /// step after step: from the corners' 1e-12 to a whole patch within a
    /// dozen steps, until the corrector no longer converged (see
    /// `predictor_corrector_step`).
    #[test]
    fn cube_minus_box_with_wide_corners() {
        let wide = |x: f64| Scalar::new(x - 1e-12, x + 1e-12);
        let corner = |p: [f64; 3]| Vector3::from_array(p.map(wide));
        let mut model = M::new();
        let block = cube(&mut model, [0.0, 0.0, 0.0], [2.0, 2.0, 1.0]);
        let tool = geop_lib_sweep::shapes::cube_solid(
            &mut model,
            &fresh_id(),
            corner([1.72, 0.1, 0.88]),
            corner([1.9, 1.9, 1.08]),
        )
        .unwrap();
        op(&mut model, block, tool, BooleanOp::Difference);
    }
}