ogeom_algo/

build.rs

//! Building topology from geometry.
//!
//! One level above [`Model`]'s raw builders: these take geometry, derive the
//! topology it implies, check the invariants that need geometry to check, and
//! report history.
//!
//! # What "checked" means here
//!
//! [`Model::add_wire`] verifies that its children are edges. That is all it can
//! do, because it has no geometry. [`make_wire`] additionally verifies that
//! consecutive edges actually *meet*, which is the property that makes a wire
//! a connected path rather than a bag of edges, and which every algorithm
//! downstream assumes without checking. A wire whose edges do not join produces
//! a face with a gap in its boundary, and the first thing to notice is usually
//! a boolean, several operations later.

use ogeom_core::FastMap;

use ogeom_core::{OgeomResult, Tolerances, ogeom_bail};
use ogeom_geom::{Curve, Curve3d, SurfaceGeometry};
use ogeom_math::{Point, Point2};
use ogeom_topo::{EdgeData, EdgeRepr, FaceData, Location, Model, Shape, ShapeType, VertexData};

use crate::history::{Built, History};

/// Roles a builder assigns, so a rebuild can match entities up.
pub mod roles {
    use ogeom_core::Role;

    /// The vertex an edge starts at.
    pub const EDGE_START: Role = Role::op_defined(0);
    /// The vertex an edge ends at.
    pub const EDGE_END: Role = Role::op_defined(1);
}

/// Add a vertex at `point`.
pub fn make_vertex(model: &mut Model, point: Point) -> Built {
    Built::from_nothing(model.add_vertex(VertexData::new(point)))
}

/// Build an edge on `curve`, bounded by its own endpoints.
///
/// Creates the two bounding vertices from the curve's ends, so the edge's
/// topology and its geometry cannot disagree about where it starts and stops. A
/// closed curve gets one vertex named twice, which is what keeps "walk to the
/// end" meaningful for a full circle.
///
/// # Errors
///
/// [`OgeomError::Domain`](ogeom_core::OgeomError::Domain) if `range` leaves the curve's
/// domain, and [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if it
/// is empty.
pub fn make_edge(
    model: &mut Model,
    curve: Curve,
    range: (f64, f64),
    tol: Tolerances,
) -> OgeomResult<Built> {
    let (lo, hi) = range;
    if !lo.is_finite() || !hi.is_finite() || hi <= lo + tol.parametric() {
        ogeom_bail!(Construction, "edge range [{lo}, {hi}] is empty");
    }
    let start = curve.point_at(lo, tol)?;
    let end = curve.point_at(hi, tol)?;

    let start_vertex = model.add_vertex(VertexData::new(start));
    // A closed edge names one vertex twice rather than two coincident ones:
    // two would leave the wire looking open at the join.
    let end_vertex = if start.is_equal(end, tol) {
        start_vertex.clone()
    } else {
        model.add_vertex(VertexData::new(end))
    };

    let id = model.geometry_mut().add_curve(curve);
    let data = EdgeData::on_curve(id, Location::identity(), range);
    let edge = model.add_edge(data, &[start_vertex.clone(), end_vertex.clone()])?;

    model.set_derived(
        &start_vertex,
        std::slice::from_ref(&edge),
        roles::EDGE_START,
    )?;
    if !end_vertex.is_same(&start_vertex) {
        model.set_derived(&end_vertex, std::slice::from_ref(&edge), roles::EDGE_END)?;
    }

    let mut history = History::new();
    history.generate(&edge, start_vertex);
    if !end_vertex.is_same(&edge) {
        history.generate(&edge, end_vertex);
    }
    Ok(Built::new(edge, history))
}

/// Build an edge between two existing vertices, on `curve`.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if either shape is
/// not a vertex, or the curve's ends do not reach them within tolerance.
pub fn make_edge_between(
    model: &mut Model,
    curve: Curve,
    range: (f64, f64),
    from: &Shape,
    to: &Shape,
    tol: Tolerances,
) -> OgeomResult<Built> {
    for v in [from, to] {
        if model.kind_of(v)? != ShapeType::Vertex {
            ogeom_bail!(Construction, "an edge is bounded by vertices");
        }
    }
    // The geometry has to actually reach the vertices it claims to join. An
    // edge whose curve stops short leaves a gap that only shows up later, in
    // whatever first tries to walk the boundary.
    let ends = [(range.0, from), (range.1, to)];
    for (parameter, vertex) in ends {
        let Some(node) = model.node(vertex) else {
            ogeom_bail!(Dangling, "vertex is not in this model");
        };
        let Some(data) = node.data().as_vertex() else {
            ogeom_bail!(Construction, "vertex node holds no point");
        };
        // Through the vertex's own placement, not against its stored point. A
        // vertex is a triple like any other shape, and the same node appears at
        // different places: the two ends of a prism are one vertex twice.
        // Comparing against the stored point would put both ends at the origin
        // of the placement and reject every edge that joins them.
        let placed = vertex.transform(model.datums())?.apply(data.point);
        let on_curve = curve.point_at(parameter, tol)?;
        let reach = data.tolerance.get().max(tol.confusion());
        if !on_curve.is_within(placed, reach) {
            if std::env::var_os("OGEOM_DEBUG_EDGE").is_some() {
                eprintln!(
                    "EDGE MISS range {range:?} at {parameter}: curve {on_curve:?} vertex {placed:?} kind {}",
                    match &curve {
                        Curve::Line(_) => "line",
                        Curve::Circle(_) => "circle",
                        Curve::BSpline(b) =>
                            if b.degree() == 1 {
                                "bspline-1"
                            } else {
                                "bspline"
                            },
                        _ => "other",
                    }
                );
            }
            ogeom_bail!(
                Construction,
                "curve at {parameter} is {} from the vertex it should meet, \
                 outside its tolerance of {reach}",
                on_curve.distance(placed)
            );
        }
    }

    let id = model.geometry_mut().add_curve(curve);
    let data = EdgeData::on_curve(id, Location::identity(), range);
    let edge = model.add_edge(data, &[from.clone(), to.clone()])?;
    Ok(Built::from_nothing(edge))
}

/// The vertices an edge runs between, in the direction it is traversed.
///
/// A reversed edge runs the other way, so its start is its stored second
/// bound. Every caller that walks a boundary needs this, and getting it from
/// the raw children instead is how a wire ends up appearing disconnected.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if `edge` is not
/// an edge.
pub fn edge_vertices(model: &Model, edge: &Shape) -> OgeomResult<Option<(Shape, Shape)>> {
    if model.kind_of(edge)? != ShapeType::Edge {
        ogeom_bail!(Construction, "expected an edge");
    }
    let bounds = model.children_of(edge)?;
    let (first, last) = match bounds.len() {
        0 => return Ok(None),
        1 => (bounds[0].clone(), bounds[0].clone()),
        _ => (bounds[0].clone(), bounds[bounds.len() - 1].clone()),
    };
    Ok(Some(
        if edge.orientation() == ogeom_topo::Orientation::Reversed {
            (last, first)
        } else {
            (first, last)
        },
    ))
}

/// An edge occurrence's curve where the occurrence stands, and the range of
/// it the edge covers; `None` for an edge with no curve in space.
///
/// The curve is moved by the curve's own location and then by the
/// occurrence's, so a reader building geometry from it builds it where the
/// edge is. A placement with a scale rescales a length-parameterized
/// curve's domain, and the range is carried with it. An unplaced edge's
/// curve comes back as stored.
///
/// # Errors
///
/// [`OgeomError::Dangling`](ogeom_core::OgeomError::Dangling) if the edge,
/// its curve or a datum is not in `model`; as
/// [`Transformable::transformed`](ogeom_geom::Transformable::transformed).
pub fn placed_edge_curve(
    model: &Model,
    edge: &Shape,
    tol: Tolerances,
) -> OgeomResult<Option<(Curve, (f64, f64))>> {
    use ogeom_geom::Transformable as _;
    let Some(data) = model.node(edge).and_then(|n| n.data().as_edge()) else {
        ogeom_bail!(Dangling, "edge is not in this model");
    };
    let Some(EdgeRepr::Curve3d {
        curve,
        location,
        range,
    }) = data.curve3d()
    else {
        return Ok(None);
    };
    let Some(stored) = model.geometry().curve(*curve) else {
        ogeom_bail!(Dangling, "curve is not in this model");
    };
    if location.is_identity() && edge.location().is_identity() {
        return Ok(Some((stored.clone(), *range)));
    }
    let placement = edge.transform(model.datums())? * location.composed(model.datums())?;
    let placed = stored.transformed(&placement, tol)?;
    let (from, to) = (stored.domain(), placed.domain());
    let range = (
        ogeom_geom::carried(range.0, from, to),
        ogeom_geom::carried(range.1, from, to),
    );
    Ok(Some((placed, range)))
}

/// An edge's curve where its occurrence places it in the face it bounds,
/// for a pcurve on that face's surface: `edge` is the occurrence as the
/// face's wire holds it. The pcurve runs on the edge's own parameter, so an
/// occurrence placed with a scale, which carries the curve to another
/// parameter, is refused.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the
/// edge has no curve in space or its placement rescales its parameter; as
/// [`placed_edge_curve`].
pub(crate) fn charted_edge_curve(
    model: &Model,
    edge: &Shape,
    tol: Tolerances,
) -> OgeomResult<(Curve, (f64, f64))> {
    let Some((curve, range)) = placed_edge_curve(model, edge, tol)? else {
        ogeom_bail!(Construction, "a face edge has no 3D curve");
    };
    let stored = model
        .node(edge)
        .and_then(|n| n.data().as_edge())
        .and_then(EdgeData::curve3d)
        .and_then(EdgeRepr::range);
    if stored != Some(range) {
        ogeom_bail!(
            Construction,
            "an edge placed with a scale runs on another parameter than its \
             curve; bake its placement before charting it on a face"
        );
    }
    Ok((curve, range))
}

/// Build a wire of straight segments through a sequence of points.
///
/// `closed` adds a final segment back to the first point, and does it by
/// naming the *first vertex again* rather than making a coincident second one,
/// which is what keeps the wire closed under [`is_wire_closed`] rather than
/// merely looking closed.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if there are fewer
/// than two points, fewer than three for a closed polygon, if consecutive
/// points coincide within tolerance (a zero-length edge is not an edge), or if
/// a closed polygon's ends do not meet.
pub fn make_polygon(
    model: &mut Model,
    points: &[Point],
    closed: bool,
    tol: Tolerances,
) -> OgeomResult<Built> {
    let least = if closed { 3 } else { 2 };
    if points.len() < least {
        ogeom_bail!(
            Construction,
            "a {} polygon needs at least {least} points, got {}",
            if closed { "closed" } else { "open" },
            points.len()
        );
    }
    for i in 1..points.len() {
        if points[i].is_equal(points[i - 1], tol) {
            ogeom_bail!(
                Construction,
                "points {} and {i} coincide, so the edge between them has no \
                 length",
                i - 1
            );
        }
    }
    // Whether or not `closed` was asked for. A caller that repeated the first
    // point at the end wants a loop, and building it as written would give the
    // wire two vertices in the same place, which every later boundary walk
    // treats as a gap that happens to be zero wide. `closed` produces the loop
    // by naming the first vertex again, which is the thing that actually
    // closes.
    if points[0].is_equal(points[points.len() - 1], tol) {
        ogeom_bail!(
            Construction,
            "the first and last points coincide; pass `closed` and drop the \
             repeat, or the wire gets two vertices in one place rather than a \
             closed loop"
        );
    }

    let mut vertices: Vec<Shape> = Vec::with_capacity(points.len());
    for point in points {
        vertices.push(model.add_vertex(VertexData::new(*point)));
    }

    let mut edges = Vec::with_capacity(points.len());
    let segments = if closed {
        points.len()
    } else {
        points.len() - 1
    };
    for i in 0..segments {
        let (from, to) = (points[i], points[(i + 1) % points.len()]);
        let curve: Curve = ogeom_geom::LineCurve::segment(from, to, tol)?.into();
        edges.push(
            make_edge_between(
                model,
                curve,
                (0.0, from.distance(to)),
                &vertices[i],
                // The closing segment names the first vertex again rather than
                // a second one at the same place.
                &vertices[(i + 1) % points.len()],
                tol,
            )?
            .shape,
        );
    }
    make_wire(model, &edges, tol)
}

/// The plane a shape lies in, if it lies in one.
///
/// Fitted to the shape's geometry and then *checked*: the fit always produces a
/// plane, and the answer is only useful if every sample is within tolerance of
/// it. `None` means the shape is not planar, which is a fact about the shape
/// rather than a failure.
///
/// Curved edges are sampled along their length, not only at their ends. A
/// circular arc's endpoints lie in a great many planes that the arc does not.
///
/// # Errors
///
/// [`OgeomError::Dangling`](ogeom_core::OgeomError::Dangling) if a handle fails to
/// resolve.
pub fn find_plane(
    model: &Model,
    shape: &Shape,
    tol: Tolerances,
) -> OgeomResult<Option<ogeom_math::Plane>> {
    let points = sample_shape(model, shape, tol)?;
    if points.len() < 3 {
        return Ok(None);
    }
    let Some((centroid, normal)) = crate::measure::least_squares_plane(&points, tol) else {
        return Ok(None);
    };
    // Fitting always yields a plane. Whether the shape is *in* it is the
    // question, and it is answered by measuring, not by having fitted.
    let reach = tol.confusion().max(
        points
            .iter()
            .map(|p| normal.dot_vector(*p - centroid).abs())
            .fold(0.0_f64, f64::max),
    );
    if reach > tol.confusion() {
        return Ok(None);
    }
    Ok(Some(ogeom_math::Plane::new(ogeom_math::Frame::about(
        centroid, normal,
    ))))
}

/// Points along a shape's geometry: every vertex, and every curved edge sampled
/// along its length.
fn sample_shape(model: &Model, shape: &Shape, tol: Tolerances) -> OgeomResult<Vec<Point>> {
    /// Enough to catch a curve leaving a candidate plane, without the cost of a
    /// real discretization; this is a yes-or-no question, not a mesh.
    const ALONG_EDGE: usize = 8;

    let mut points = Vec::new();
    for vertex in ogeom_topo::explore_unique(model, shape, ShapeType::Vertex)? {
        if let Some(data) = model.node(&vertex).and_then(|n| n.data().as_vertex()) {
            points.push(vertex.transform(model.datums())?.apply(data.point));
        }
    }
    for edge in ogeom_topo::explore_unique(model, shape, ShapeType::Edge)? {
        let Some(data) = model.node(&edge).and_then(|n| n.data().as_edge()) else {
            continue;
        };
        let Some(EdgeRepr::Curve3d { curve, range, .. }) = data.curve3d() else {
            continue;
        };
        let Some(geometry) = model.geometry().curve(*curve) else {
            continue;
        };
        if geometry.kind() == ogeom_geom::CurveKind::Line {
            continue;
        }
        let placement = edge.transform(model.datums())?;
        for i in 1..ALONG_EDGE {
            #[allow(clippy::cast_precision_loss)]
            let t = i as f64 / ALONG_EDGE as f64;
            let at = range.0 + (range.1 - range.0) * t;
            points.push(placement.apply(geometry.point_at(at, tol)?));
        }
    }
    Ok(points)
}

/// Build a wire from edges that meet end to end.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if a child is not
/// an edge, the list is empty, or consecutive edges do not share a vertex.
pub fn make_wire(model: &mut Model, edges: &[Shape], tol: Tolerances) -> OgeomResult<Built> {
    if edges.is_empty() {
        ogeom_bail!(Construction, "a wire needs at least one edge");
    }
    check_connected(model, edges, tol)?;

    let wire = model.add_wire(edges)?;
    let mut history = History::new();
    for edge in edges {
        // The edges are not consumed: they are shared, and an edge between two
        // faces belongs to both wires. Reporting them modified into the wire
        // would claim the edge ceased to exist.
        history.generate(edge, wire.clone());
    }
    Ok(Built::new(wire, history))
}

/// Verify that consecutive edges share a vertex.
fn check_connected(model: &Model, edges: &[Shape], tol: Tolerances) -> OgeomResult<()> {
    let mut ends: Vec<Option<(Shape, Shape)>> = Vec::with_capacity(edges.len());
    for edge in edges {
        ends.push(edge_vertices(model, edge)?);
    }
    for i in 0..edges.len().saturating_sub(1) {
        let (Some((_, end)), Some((next_start, _))) = (&ends[i], &ends[i + 1]) else {
            // An unbounded edge cannot be shown to join anything, and pretending
            // otherwise is worse than saying so.
            ogeom_bail!(
                Construction,
                "edge {i} or {} has no bounding vertices, so the wire cannot be \
                 shown to connect",
                i + 1
            );
        };
        if !end.is_same(next_start)
            && !model.same_position(end, next_start, tol)?
            && !one_point(model, end, next_start, tol)?
        {
            if std::env::var("OGEOM_DEBUG_WIRE").is_ok()
                && let (Some(a), Some(b)) = (
                    model.node(end).and_then(|n| n.data().as_vertex().cloned()),
                    model
                        .node(next_start)
                        .and_then(|n| n.data().as_vertex().cloned()),
                )
            {
                eprintln!(
                    "WIRE GAP: {:?} (tol {:.2e}) vs {:?} (tol {:.2e}), gap {:.3e}",
                    a.point,
                    a.tolerance.get(),
                    b.point,
                    b.tolerance.get(),
                    a.point.distance(b.point)
                );
            }
            ogeom_bail!(
                Construction,
                "edge {i} ends where edge {} does not begin; a wire whose edges \
                 do not meet leaves a gap in every face built on it",
                i + 1
            );
        }
    }
    Ok(())
}

/// Whether a wire's last edge returns to its first edge's start.
///
/// # Errors
///
/// As [`edge_vertices`].
pub fn is_wire_closed(model: &Model, wire: &Shape, tol: Tolerances) -> OgeomResult<bool> {
    let edges = model.children_of(wire)?;
    let (Some(first), Some(last)) = (edges.first(), edges.last()) else {
        return Ok(false);
    };
    let (Some((start, _)), Some((_, end))) =
        (edge_vertices(model, first)?, edge_vertices(model, last)?)
    else {
        return Ok(false);
    };
    Ok(start.is_same(&end)
        || model.same_position(&start, &end, tol)?
        || one_point(model, &start, &end, tol)?)
}

/// Whether two occurrences of one vertex node land on one point: a vertex
/// on a revolution's axis, turned, is itself under another placement.
pub(crate) fn one_point(model: &Model, a: &Shape, b: &Shape, tol: Tolerances) -> OgeomResult<bool> {
    if a.node() != b.node() {
        return Ok(false);
    }
    let Some(data) = model.node(a).and_then(|n| n.data().as_vertex()) else {
        return Ok(false);
    };
    let pa = a.transform(model.datums())?.apply(data.point);
    let pb = b.transform(model.datums())?.apply(data.point);
    Ok(pa.distance(pb) <= data.tolerance.get().max(tol.confusion()))
}

/// A fresh wire walking `wire`'s edges the other way round.
pub(crate) fn walked_back(model: &mut Model, wire: &Shape, tol: Tolerances) -> OgeomResult<Shape> {
    let edges = model.ordered_children_of(&wire.reversed())?;
    Ok(make_wire(model, &edges, tol)?.shape)
}

/// How a wire turns about `axis`: twice the area it encloses projected
/// square to the axis, signed by the right-hand rule, from points sampled
/// along its edges in traversal order.
pub(crate) fn wire_turn(
    model: &Model,
    wire: &Shape,
    axis: ogeom_math::Vector,
    tol: Tolerances,
) -> OgeomResult<f64> {
    edges_turn(model, &model.ordered_children_of(wire)?, axis, tol)
}

/// How a closed walk of edge occurrences turns about `axis`: twice the
/// area it encloses projected square to the axis, signed by the right-hand
/// rule, from points sampled along the edges in traversal order. Positive
/// for a ring walked counter-clockwise about `axis`.
///
/// # Errors
///
/// As the model's placements and the curves' evaluation.
pub fn edges_turn(
    model: &Model,
    edges: &[Shape],
    axis: ogeom_math::Vector,
    tol: Tolerances,
) -> OgeomResult<f64> {
    let mut points: Vec<ogeom_math::Point> = Vec::new();
    for edge in edges {
        let Some(EdgeRepr::Curve3d { curve, range, .. }) = model
            .node(edge)
            .and_then(|n| n.data().as_edge())
            .and_then(|d| d.curve3d())
        else {
            continue;
        };
        let Some(geometry) = model.geometry().curve(*curve) else {
            continue;
        };
        let placement = edge.transform(model.datums())?;
        const SAMPLES: u32 = 16;
        for k in 0..SAMPLES {
            let f = f64::from(k) / f64::from(SAMPLES);
            let t = if edge.orientation() == ogeom_topo::Orientation::Reversed {
                range.1 + (range.0 - range.1) * f
            } else {
                range.0 + (range.1 - range.0) * f
            };
            points.push(placement.apply(geometry.point_at(t, tol)?));
        }
    }
    // About the walk's own first point: the sum is the same about any
    // point, and about one far off it is a difference of huge products.
    let Some(&anchor) = points.first() else {
        return Ok(0.0);
    };
    let mut newell = ogeom_math::Vector::ZERO;
    for i in 0..points.len() {
        let (a, b) = (points[i], points[(i + 1) % points.len()]);
        newell += (a - anchor).cross(b - anchor);
    }
    Ok(newell.dot(axis))
}

/// Build a face on `surface`, bounded by `wires`.
///
/// The first wire is the outer boundary; any others are holes. Every wire must
/// be closed, since an open boundary encloses nothing. On a plane the face
/// faces the plane's normal and walks its outer boundary counter-clockwise
/// about it and its holes clockwise, whichever way the wires were walked.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if a wire is open
/// or is not a wire.
pub fn make_face(
    model: &mut Model,
    surface: SurfaceGeometry,
    wires: &[Shape],
    tol: Tolerances,
) -> OgeomResult<Built> {
    for (i, wire) in wires.iter().enumerate() {
        if model.kind_of(wire)? != ShapeType::Wire {
            ogeom_bail!(Construction, "a face is bounded by wires");
        }
        if !is_wire_closed(model, wire, tol)? {
            ogeom_bail!(
                Construction,
                "wire {i} is open; an open boundary encloses no area"
            );
        }
    }

    let wound = wound_to_plane(model, &surface, wires, tol)?;
    let id = model.geometry_mut().add_surface(surface);
    let mut built = make_face_on(model, id, &wound, tol)?;
    for (given, used) in wires.iter().zip(&wound) {
        if !given.is_same(used) {
            built.history.generate(given, built.shape.clone());
        }
    }
    Ok(built)
}

/// `wires` as a face on `surface` walks them: on a plane, the outer ring
/// counter-clockwise about the plane's normal and each hole clockwise, so
/// the material lies on the left of every ring seen from the side the
/// face faces. A ring wound the other way is replaced by a fresh wire
/// walking its edges back; the plane, and so the face's normal, is kept.
/// Any other surface, and a ring enclosing no area, is left as given.
fn wound_to_plane(
    model: &mut Model,
    surface: &SurfaceGeometry,
    wires: &[Shape],
    tol: Tolerances,
) -> OgeomResult<Vec<Shape>> {
    use ogeom_geom::Surface as _;
    let SurfaceGeometry::Plane(plane) = surface else {
        return Ok(wires.to_vec());
    };
    let normal = plane.normal_at(0.0, 0.0, tol)?.vector();
    let mut wound = Vec::with_capacity(wires.len());
    for (i, wire) in wires.iter().enumerate() {
        let turn = wire_turn(model, wire, normal, tol)?;
        let want = if i == 0 { 1.0 } else { -1.0 };
        // Twice the enclosed area: a ring closing on itself within the
        // confusion distance encloses nothing to wind about.
        let against = turn * want < -tol.confusion() * tol.confusion();
        wound.push(if against {
            walked_back(model, wire, tol)?
        } else {
            wire.clone()
        });
    }
    Ok(wound)
}

/// Build a face on a surface the model already holds.
///
/// The distinction from [`make_face`] matters more than it looks: a pcurve
/// names the surface it is drawn on by id, so an edge's pcurve and the face
/// bounded by that edge have to name the *same* id. Registering the surface
/// twice gives two ids for one surface, and every lookup of "the pcurve on this
/// face" comes back empty even though the pcurve is right there.
///
/// # Errors
///
/// As [`make_face`].
pub fn make_face_on(
    model: &mut Model,
    surface: ogeom_topo::SurfaceId,
    wires: &[Shape],
    tol: Tolerances,
) -> OgeomResult<Built> {
    for (i, wire) in wires.iter().enumerate() {
        if model.kind_of(wire)? != ShapeType::Wire {
            ogeom_bail!(Construction, "a face is bounded by wires");
        }
        if !is_wire_closed(model, wire, tol)? {
            ogeom_bail!(
                Construction,
                "wire {i} is open; an open boundary encloses no area"
            );
        }
    }

    let data = FaceData::new(surface, Location::identity());
    let face = model.add_face(data, wires)?;

    let mut history = History::new();
    for wire in wires {
        history.generate(wire, face.clone());
    }
    Ok(Built::new(face, history))
}

/// Build a face on `surface` from per-wire edge lists, attaching an exact
/// same-parameter pcurve to every edge.
///
/// The construction path for faces whose curves were *chosen* to have
/// closed-form charts: blend wedges, offset rebuilds. Every edge's curve
/// must lie on the surface in a configuration
/// [`ogeom_intersect::exact_pcurve_of`] recognises; a fitted pcurve here would
/// manufacture disagreement where none exists, so an edge with no closed
/// form is refused instead.
///
/// Each ring is walked as its edges are listed, on a plane too: a rebuilt
/// face whose ring winds against its surface was driven through itself,
/// and the face keeps that for the checks downstream to find.
///
/// # Errors
///
/// As [`make_wire`] and [`make_face`], and
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if an edge has
/// no 3D curve or no closed-form pcurve on `surface`.
pub fn make_face_with_pcurves(
    model: &mut Model,
    surface: SurfaceGeometry,
    wires: &[Vec<Shape>],
    tol: Tolerances,
) -> OgeomResult<Built> {
    let mut rings: Vec<Shape> = Vec::with_capacity(wires.len());
    for edges in wires {
        rings.push(make_wire(model, edges, tol)?.shape);
    }
    let id = model.geometry_mut().add_surface(surface.clone());
    let built = make_face_on(model, id, &rings, tol)?;
    let surface_id = {
        let Some(node) = model.node(&built.shape) else {
            ogeom_bail!(Dangling, "the face just built is not in this model");
        };
        let ogeom_topo::NodeData::Face(data) = node.data() else {
            ogeom_bail!(Construction, "the face holds no face data");
        };
        data.surface
    };
    // Walked in each wire's own traversal order, so that on a periodic
    // chart every pcurve lands on the *same branch* as its neighbour: an
    // exact inversion answers with whatever phase its construction likes,
    // and two images a full period apart describe the same points while
    // tearing the boundary the arrangement walks. Each pcurve after the
    // first is shifted by whole periods until its start meets the previous
    // traversal's end.
    use ogeom_geom::Curve2d as _;
    use ogeom_geom::Surface as _;
    let ((ua, ub), (va, vb)) = surface.domain();
    let u_period = surface.is_periodic_u().then_some(ub - ua);
    let v_period = surface.is_periodic_v().then_some(vb - va);
    let mut done: Vec<Shape> = Vec::new();
    for edges in wires {
        let mut prev_end: Option<ogeom_math::Point2> = None;
        for edge in edges {
            let (curve, prange) = charted_edge_curve(model, edge, tol)?;
            let mut pcurve = match ogeom_intersect::exact_pcurve_over(&curve, prange, &surface, tol)
            {
                Some(exact) => exact,
                // No closed form: a fitted surface, mostly. The edge lies on
                // the surface by construction here (a rebuilt boundary, a
                // recovered intersection), so the projected fit speaks it: the
                // same machinery the exchange readers trust, at the same cap.
                // The edge is widened to the offset its samples sat at and to
                // how far the fitted trim, lifted, stands from the curve
                // anywhere along it.
                None => {
                    let (fitted, _, _, worst_off, _) =
                        crate::pcurve_fit::fit_projected_pcurve(&curve, prange, &surface, tol)?;
                    let gap = crate::pcurve_fit::lifted_gap(
                        (&curve, prange),
                        (&fitted, prange),
                        &surface,
                        false,
                        tol,
                    )?;
                    let worst_off = worst_off.max(gap);
                    if worst_off > tol.confusion() {
                        // The edge owns the offset, and so must the vertices
                        // that bound it: the containment rule.
                        let widened = ogeom_core::Tolerance::new(worst_off + tol.confusion())?;
                        model.widen(edge, widened)?;
                        if let Some((a, b)) = edge_vertices(model, edge)? {
                            model.widen(&a, widened)?;
                            model.widen(&b, widened)?;
                        }
                    }
                    fitted
                }
            };
            let (t_start, t_end) = if edge.orientation() == ogeom_topo::Orientation::Reversed {
                (prange.1, prange.0)
            } else {
                (prange.0, prange.1)
            };
            if u_period.is_some() || v_period.is_some() {
                let start = pcurve.point_at(t_start, tol)?;
                if let Some(prev) = prev_end {
                    let shift = ogeom_math::Vector2::new(
                        u_period.map_or(0.0, |p| ((prev.x - start.x) / p).round() * p),
                        v_period.map_or(0.0, |p| ((prev.y - start.y) / p).round() * p),
                    );
                    if shift.x != 0.0 || shift.y != 0.0 {
                        pcurve =
                            pcurve.transformed(&ogeom_math::Transform2::translation(shift), tol)?;
                    }
                }
                prev_end = Some(pcurve.point_at(t_end, tol)?);
            }
            if done.iter().any(|e| e.is_same(edge)) {
                continue;
            }
            done.push(edge.clone());
            attach_pcurve(
                model,
                edge,
                pcurve,
                surface_id,
                edge.location().clone(),
                prange,
            )?;
        }
    }
    Ok(built)
}

/// Put a face's wires on one branch of a periodic chart.
///
/// An exchange file's pcurves, or a fit of them, answer with whatever
/// phase the inversion likes: a hole loop that straddles a drum's seam
/// comes back half on each branch (a loop that never closes in the chart,
/// which no mesher can cut), and a slitted wall's outer wire hops branches
/// at every slit. Walked in each wire's own order, an image whose start
/// misses the previous traversal's end by close to a whole period is
/// shifted by that period; a miss of half a period is a pole, where two
/// meridians meet at one point with no edge between them, and is left
/// alone. A seam representation carries both branches by design: the
/// walk passes through it on the side the occurrence uses and moves
/// nothing. Every wire after the first is then carried whole onto the
/// first wire's branch, the middle of the extent its images run over
/// brought within half a period of the first's, so rims, slits, holes and
/// rings closed by one edge all read in one chart.
/// The images are rewritten in place (`GeometryStore::pcurve_mut`), once
/// each; a slit uses one image twice.
///
/// # Errors
///
/// [`OgeomError::Dangling`](ogeom_core::OgeomError::Dangling) if a wire or
/// edge is not in this model.
pub fn chain_wire_branches(
    model: &mut Model,
    surface: ogeom_topo::SurfaceId,
    wires: &[Shape],
    tol: Tolerances,
) -> OgeomResult<()> {
    use ogeom_geom::Curve2d as _;
    use ogeom_geom::Surface as _;
    let Some(geometry) = model.geometry().surface(surface) else {
        ogeom_bail!(Dangling, "surface is not in this model");
    };
    let ((ua, ub), (va, vb)) = geometry.domain();
    let u_period = geometry.is_periodic_u().then_some(ub - ua);
    let v_period = geometry.is_periodic_v().then_some(vb - va);
    if u_period.is_none() && v_period.is_none() {
        return Ok(());
    }
    // A whole number of periods, when the gap is within a quarter period
    // of one; nothing for a pole's half-period hop.
    let whole = |gap: f64, period: Option<f64>| -> f64 {
        period.map_or(0.0, |p| {
            let k = (gap / p).round();
            if k != 0.0 && (gap - k * p).abs() <= p * 0.25 {
                k * p
            } else {
                0.0
            }
        })
    };
    // A wire's images in traversal order: the pcurve to move (none for a
    // seam), its start and end in the edge's own direction, and points
    // along it, which place a loop closed by one edge whose ends say
    // nothing of where it runs.
    type Image = (
        Option<ogeom_topo::PCurveId>,
        ogeom_math::Point2,
        ogeom_math::Point2,
        Vec<ogeom_math::Point2>,
    );
    /// How many steps each image is sampled at for the wire's extent.
    const STEPS: u32 = 8;
    let images = |model: &Model, wire: &Shape| -> OgeomResult<Vec<Image>> {
        let mut out = Vec::new();
        for edge in model.ordered_children_of(wire)? {
            let Some(data) = model.node(&edge).and_then(|n| n.data().as_edge()) else {
                ogeom_bail!(Dangling, "edge is not in this model");
            };
            let reversed = edge.orientation() == ogeom_topo::Orientation::Reversed;
            let (id, movable, range) = match data.pcurve_for(surface, edge.location()) {
                Some(EdgeRepr::PCurve { curve, range, .. }) => (*curve, true, *range),
                Some(EdgeRepr::Seam {
                    forward,
                    reversed: back,
                    range,
                    ..
                }) => (if reversed { *back } else { *forward }, false, *range),
                _ => continue,
            };
            let Some(planar) = model.geometry().pcurve(id) else {
                ogeom_bail!(Dangling, "pcurve is not in this model");
            };
            let (t_start, t_end) = if reversed {
                (range.1, range.0)
            } else {
                (range.0, range.1)
            };
            let mut along = Vec::with_capacity(STEPS as usize - 1);
            for k in 1..STEPS {
                let t = range.0 + (range.1 - range.0) * f64::from(k) / f64::from(STEPS);
                along.push(planar.point_at(t, tol)?);
            }
            out.push((
                movable.then_some(id),
                planar.point_at(t_start, tol)?,
                planar.point_at(t_end, tol)?,
                along,
            ));
        }
        Ok(out)
    };
    let mut first_centre: Option<ogeom_math::Point2> = None;
    for wire in wires {
        let wire_images = images(model, wire)?;
        if wire_images.is_empty() {
            continue;
        }
        // The chain, dry: each image's shift so its start meets the
        // previous end, an image already placed taking its earlier shift.
        let mut shifts: Vec<(ogeom_topo::PCurveId, ogeom_math::Vector2)> = Vec::new();
        let mut prev_end: Option<ogeom_math::Point2> = None;
        let mut lo = ogeom_math::Point2::new(f64::INFINITY, f64::INFINITY);
        let mut hi = ogeom_math::Point2::new(f64::NEG_INFINITY, f64::NEG_INFINITY);
        for (id, start, end, along) in &wire_images {
            let shift = match id {
                None => ogeom_math::Vector2::new(0.0, 0.0),
                Some(id) => match shifts.iter().find(|(seen, _)| seen == id) {
                    Some((_, s)) => *s,
                    None => {
                        let s = prev_end.map_or(ogeom_math::Vector2::new(0.0, 0.0), |prev| {
                            ogeom_math::Vector2::new(
                                whole(prev.x - start.x, u_period),
                                whole(prev.y - start.y, v_period),
                            )
                        });
                        shifts.push((*id, s));
                        s
                    }
                },
            };
            let (s, e) = (*start + shift, *end + shift);
            for p in along.iter().map(|p| *p + shift).chain([s, e]) {
                lo = ogeom_math::Point2::new(lo.x.min(p.x), lo.y.min(p.y));
                hi = ogeom_math::Point2::new(hi.x.max(p.x), hi.y.max(p.y));
            }
            prev_end = Some(e);
        }
        // Then the whole wire onto the first wire's branch.
        let centre = ogeom_math::Point2::new(f64::midpoint(lo.x, hi.x), f64::midpoint(lo.y, hi.y));
        let carry = match first_centre {
            None => {
                first_centre = Some(centre);
                ogeom_math::Vector2::new(0.0, 0.0)
            }
            Some(first) => ogeom_math::Vector2::new(
                u_period.map_or(0.0, |p| ((first.x - centre.x) / p).round() * p),
                v_period.map_or(0.0, |p| ((first.y - centre.y) / p).round() * p),
            ),
        };
        for (id, shift) in shifts {
            let total = shift + carry;
            if total.x == 0.0 && total.y == 0.0 {
                continue;
            }
            let Some(planar) = model.geometry().pcurve(id) else {
                ogeom_bail!(Dangling, "pcurve is not in this model");
            };
            let shifted = planar.transformed(&ogeom_math::Transform2::translation(total), tol)?;
            if let Some(slot) = model.geometry_mut().pcurve_mut(id) {
                *slot = shifted;
            }
        }
    }
    Ok(())
}

/// A wire's edge uses with the chart closed where the surface collapses.
///
/// A face whose boundary passes through a pole or an apex meets it in space
/// at one vertex, but in the chart the use arriving there and the use
/// leaving it sit apart along the row (or column) that collapses to that
/// point: a sphere bounded by its seam walked both ways arrives at the north
/// pole on one side of the chart and leaves it on the other. An exchange
/// file has no edge for that stretch, and without one the boundary does not
/// close in the chart. Between consecutive uses that share a vertex and
/// whose images sit apart on a stretch the surface maps to that one point, a
/// degenerate edge is put in, its pcurve the straight run across. A wire
/// whose images break anywhere else (an image a period off its neighbours,
/// which the mesher folds back) has no run to trust, and comes back as it
/// is, as does one with a use that has no image on `surface`.
///
/// # Errors
///
/// [`OgeomError::Dangling`](ogeom_core::OgeomError::Dangling) if the surface,
/// an edge or a pcurve is not in this model.
pub fn closed_at_poles(
    model: &mut Model,
    surface: ogeom_topo::SurfaceId,
    edges: &[Shape],
    tol: Tolerances,
) -> OgeomResult<Vec<Shape>> {
    use ogeom_geom::Curve2d as _;
    use ogeom_geom::Surface as _;
    let Some(geometry) = model.geometry().surface(surface).cloned() else {
        ogeom_bail!(Dangling, "surface is not in this model");
    };
    let ((ua, ub), (va, vb)) = geometry.domain();
    let span = (ub - ua).abs().max((vb - va).abs());
    // Each use's chart start and end, in the direction it is walked.
    let mut ends = Vec::with_capacity(edges.len());
    for edge in edges {
        let Some(data) = model.node(edge).and_then(|n| n.data().as_edge()) else {
            ogeom_bail!(Dangling, "edge is not in this model");
        };
        let reversed = edge.orientation() == ogeom_topo::Orientation::Reversed;
        let (id, range) = match data.pcurve_for(surface, edge.location()) {
            Some(EdgeRepr::PCurve { curve, range, .. }) => (*curve, *range),
            Some(EdgeRepr::Seam {
                forward,
                reversed: back,
                range,
                ..
            }) => (if reversed { *back } else { *forward }, *range),
            _ => return Ok(edges.to_vec()),
        };
        let Some(planar) = model.geometry().pcurve(id) else {
            ogeom_bail!(Dangling, "pcurve is not in this model");
        };
        let (t_start, t_end) = if reversed {
            (range.1, range.0)
        } else {
            (range.0, range.1)
        };
        ends.push((planar.point_at(t_start, tol)?, planar.point_at(t_end, tol)?));
    }
    let collapses = |from: Point2, to: Point2| -> bool {
        let Ok(at) = geometry.point_at(from.x, from.y, tol) else {
            return false;
        };
        (1..=4).all(|k| {
            let s = f64::from(k) / 4.0;
            let p = Point2::new(from.x + (to.x - from.x) * s, from.y + (to.y - from.y) * s);
            geometry
                .point_at(p.x, p.y, tol)
                .is_ok_and(|q| q.distance(at) <= tol.confusion())
        })
    };
    // Each junction: whether the images meet there, or which pole run
    // closes it.
    let mut runs: Vec<Option<(Shape, Point2, Point2)>> = Vec::with_capacity(edges.len());
    for (i, edge) in edges.iter().enumerate() {
        let j = (i + 1) % edges.len();
        // A miss that small is fitting slop between images that meet, and
        // collapses to a point on any surface. So is one too short in the
        // chart to be an edge of its own whose ends stand at one point in
        // space: on a chart spanning a fortieth of a unit, a millionth of
        // the span is shorter than the confusion distance.
        let (from, to) = (ends[i].1, ends[j].0);
        let gap = from.distance(to);
        let meets_in_space = || match (
            geometry.point_at(from.x, from.y, tol),
            geometry.point_at(to.x, to.y, tol),
        ) {
            (Ok(p), Ok(q)) => p.distance(q) <= tol.confusion(),
            _ => false,
        };
        if gap <= span * 1e-6 || (gap <= tol.confusion() && meets_in_space()) {
            runs.push(None);
            continue;
        }
        let (Some((_, arriving)), Some((leaving, _))) = (
            edge_vertices(model, edge)?,
            edge_vertices(model, &edges[j])?,
        ) else {
            return Ok(edges.to_vec());
        };
        if arriving.node() != leaving.node() || !collapses(from, to) {
            return Ok(edges.to_vec());
        }
        runs.push(Some((arriving, from, to)));
    }
    let mut out = Vec::with_capacity(edges.len() + 2);
    for (edge, run) in edges.iter().zip(runs) {
        out.push(edge.clone());
        let Some((at, from, to)) = run else {
            continue;
        };
        let mut data = EdgeData::new();
        data.degenerate = true;
        let pole = model.add_edge(data, &[at.clone(), at])?;
        attach_pcurve(
            model,
            &pole,
            ogeom_geom::Line2d::segment(from, to, tol)?.into(),
            surface,
            Location::identity(),
            (0.0, from.distance(to)),
        )?;
        out.push(pole);
    }
    Ok(out)
}

/// The one wire a face on a B-spline closed round in `u` takes when its
/// file bounds it by rings alone and writes no seam: each ring runs the
/// whole period along the chart, from one end column to the other, and
/// without a seam no ring meets another, so the boundary encloses nothing
/// in the chart. The seam is the surface's own column at the period's
/// join, exact (the iso-curve the control columns trace), from where one
/// ring's walk ends to where the next begins; walked up one end column and
/// back down the other, it closes the two rings into one loop. A single
/// ring opposite a row the surface collapses to a point (a dome's pole)
/// takes the seam to that point and back, and the pole's run across the
/// chart is laid by [`closed_at_poles`].
///
/// `None` where the face is not that shape: another surface, one closed
/// the other way or periodic, a wire of more than one edge, or a ring
/// whose image does not run end column to end column along one row.
///
/// # Errors
///
/// [`OgeomError::Dangling`](ogeom_core::OgeomError::Dangling) if the
/// surface, an edge or a pcurve is not in this model.
pub fn belted_wire(
    model: &mut Model,
    surface: ogeom_topo::SurfaceId,
    wires: &[Shape],
    tol: Tolerances,
) -> OgeomResult<Option<Shape>> {
    use ogeom_geom::Curve2d as _;
    use ogeom_geom::Surface as _;
    let Some(SurfaceGeometry::BSpline(spline)) = model.geometry().surface(surface).cloned() else {
        return Ok(None);
    };
    if !spline.is_closed_u(tol) || spline.is_periodic_u() || spline.is_closed_v(tol) {
        return Ok(None);
    }
    if wires.is_empty() || wires.len() > 2 {
        return Ok(None);
    }
    let ((u0, u1), (v0, v1)) = spline.domain();
    let slack = (u1 - u0).abs().max((v1 - v0).abs()) * 1e-6;
    // Each ring: its use, the vertex it starts and ends at, and its walk's
    // row, start column and end column.
    struct Belt {
        edge: Shape,
        vertex: Shape,
        row: f64,
        from: f64,
        to: f64,
    }
    let mut belts = Vec::with_capacity(2);
    for wire in wires {
        let edges = model.ordered_children_of(wire)?;
        let [edge] = edges.as_slice() else {
            return Ok(None);
        };
        let Some((a, b)) = edge_vertices(model, edge)? else {
            return Ok(None);
        };
        if !a.is_same(&b) {
            return Ok(None);
        }
        let Some(data) = model.node(edge).and_then(|n| n.data().as_edge()) else {
            ogeom_bail!(Dangling, "edge is not in this model");
        };
        if data.degenerate {
            return Ok(None);
        }
        let Some(EdgeRepr::PCurve { curve, range, .. }) = data.pcurve_for(surface, edge.location())
        else {
            return Ok(None);
        };
        let Some(planar) = model.geometry().pcurve(*curve) else {
            ogeom_bail!(Dangling, "pcurve is not in this model");
        };
        let reversed = edge.orientation() == ogeom_topo::Orientation::Reversed;
        let (s, e) = (
            planar.point_at(range.0, tol)?,
            planar.point_at(range.1, tol)?,
        );
        let (s, e) = if reversed { (e, s) } else { (s, e) };
        let column = |x: f64| -> Option<f64> {
            if (x - u0).abs() <= slack {
                Some(u0)
            } else if (x - u1).abs() <= slack {
                Some(u1)
            } else {
                None
            }
        };
        let (Some(from), Some(to)) = (column(s.x), column(e.x)) else {
            return Ok(None);
        };
        if from == to || (s.y - e.y).abs() > slack {
            return Ok(None);
        }
        belts.push(Belt {
            edge: edge.clone(),
            vertex: a,
            row: f64::midpoint(s.y, e.y),
            from,
            to,
        });
    }
    // The far end of the seam: the second ring's vertex, or the point the
    // row opposite a lone ring collapses to.
    let (end_vertex, end_row) = if let [_, second] = belts.as_slice() {
        if second.from != belts[0].to {
            return Ok(None);
        }
        (second.vertex.clone(), second.row)
    } else {
        let row = if (belts[0].row - v0).abs() <= (belts[0].row - v1).abs() {
            v1
        } else {
            v0
        };
        let across = spline.iso_v_curve(row, tol)?;
        let points: Vec<Point> = (0..=8)
            .map(|k| {
                let t = u0 + (u1 - u0) * f64::from(k) / 8.0;
                ogeom_geom::Curve3d::point_at(&across, t, tol)
            })
            .collect::<OgeomResult<_>>()?;
        if points
            .iter()
            .any(|p| p.distance(points[0]) > tol.confusion())
        {
            return Ok(None);
        }
        (model.add_vertex(VertexData::new(points[0])), row)
    };
    let first = &belts[0];
    let up = first.row < end_row;
    let (lo, hi) = if up {
        (first.vertex.clone(), end_vertex)
    } else {
        (end_vertex, first.vertex.clone())
    };
    let column_curve: Curve = spline.iso_u_curve(first.to, tol)?.into();
    let range = (first.row.min(end_row), first.row.max(end_row));
    let seam = make_edge_between(model, column_curve, range, &lo, &hi, tol)?.shape;
    let column = |u: f64| -> OgeomResult<ogeom_geom::PlanarCurve> {
        Ok(ogeom_geom::Line2d::over(
            ogeom_math::Axis2::new(
                Point2::new(u, 0.0),
                ogeom_math::Direction2::new(ogeom_math::Vector2::new(0.0, 1.0), tol)?,
            ),
            range.0,
            range.1,
        )?
        .into())
    };
    // The seam's Forward use climbs; walked up the first ring's end column
    // where the walk climbs, and down its start column where it falls.
    let (forward_col, reversed_col) = if up {
        (first.to, first.from)
    } else {
        (first.from, first.to)
    };
    attach_seam(
        model,
        &seam,
        column(forward_col)?,
        column(reversed_col)?,
        surface,
        Location::identity(),
        range,
    )?;
    let climb = if up { seam.clone() } else { seam.reversed() };
    let mut walk = vec![first.edge.clone(), climb.clone()];
    if let [_, second] = belts.as_slice() {
        walk.push(second.edge.clone());
    }
    walk.push(climb.reversed());
    Ok(Some(make_wire(model, &walk, tol)?.shape))
}

/// Build a face covering the whole of `surface`, with no trimming.
///
/// # Errors
///
/// Never fails for a well-formed surface; the signature matches [`make_face`]
/// so the two are interchangeable at a call site.
pub fn make_natural_face(model: &mut Model, surface: SurfaceGeometry) -> OgeomResult<Built> {
    let id = model.geometry_mut().add_surface(surface);
    let face = model.add_face(FaceData::natural(id, Location::identity()), &[])?;
    Ok(Built::from_nothing(face))
}

/// Build a shell from faces.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if a child is not
/// a face, or the list is empty.
/// Build a compsolid: solids gluing along shared faces.
///
/// A compsolid is more than a bag of solids: its members must actually
/// glue. Every pair connected through the whole is connected through shared
/// *face nodes*: the same face entity bounding two solids, once from each
/// side, which is the sharing a boolean or a sew produces. A set of solids
/// that merely touch, faces coincident but distinct, is a compound, not a
/// compsolid, and is refused as one.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if
/// fewer than two solids are given, a member is not a solid, or the members
/// do not glue into one connected whole through shared faces.
pub fn make_compsolid(model: &mut Model, solids: &[Shape]) -> OgeomResult<Built> {
    if solids.len() < 2 {
        ogeom_bail!(
            Construction,
            "a compsolid glues at least two solids; one solid is just a solid"
        );
    }
    // Which face nodes bound each solid.
    let mut faces_of: Vec<ogeom_core::FastSet<ogeom_topo::TShapeId>> = Vec::new();
    for solid in solids {
        if model.kind_of(solid)? != ShapeType::Solid {
            ogeom_bail!(Construction, "a compsolid's members are solids");
        }
        let set = ogeom_topo::explore(model, solid, ogeom_topo::Filter::OfType(ShapeType::Face))?
            .iter()
            .map(Shape::node)
            .collect();
        faces_of.push(set);
    }
    // Connectivity over shared face nodes: a flood fill from the first.
    let mut joined = vec![false; solids.len()];
    joined[0] = true;
    let mut frontier = vec![0_usize];
    while let Some(current) = frontier.pop() {
        for (other, other_faces) in faces_of.iter().enumerate() {
            if joined[other] {
                continue;
            }
            if !faces_of[current].is_disjoint(other_faces) {
                joined[other] = true;
                frontier.push(other);
            }
        }
    }
    if joined.iter().any(|j| !j) {
        ogeom_bail!(
            Construction,
            "the solids do not glue into one connected whole: at least one              shares no face with the rest. Coincident-but-distinct faces are              a compound's arrangement, not a compsolid's; sew or fuse first."
        );
    }

    let compsolid = model.add_compsolid(solids)?;
    let mut history = History::new();
    for solid in solids {
        history.generate(solid, compsolid.clone());
    }
    Ok(Built::new(compsolid, history))
}

/// Build a compound from any shapes, with history.
///
/// The grouping container: a compound holds anything, orders nothing, and
/// claims nothing about closure. What this adds over the raw model call is
/// the same thing every builder adds: a history that says what went in.
///
/// # Errors
///
/// As [`Model::add_compound`].
pub fn make_compound(model: &mut Model, shapes: &[Shape]) -> OgeomResult<Built> {
    let compound = model.add_compound(shapes)?;
    let mut history = History::new();
    for shape in shapes {
        history.generate(shape, compound.clone());
    }
    Ok(Built::new(compound, history))
}

/// Build a shell from faces.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if a child is
/// not a face, or the list is empty.
pub fn make_shell(model: &mut Model, faces: &[Shape]) -> OgeomResult<Built> {
    let shell = model.add_shell(faces)?;
    let mut history = History::new();
    for face in faces {
        history.generate(face, shell.clone());
    }
    Ok(Built::new(shell, history))
}

/// Build a solid from shells.
///
/// The first shell is the outer boundary; any others are voids inside it.
///
/// Every face of a solid faces away from its material: the outer shell
/// out of the volume it bounds, a void's shell into its cavity. A closed
/// shell whose faces agree with each other (every edge walked once each
/// way) but face the other way is used turned over: the outer shell where
/// its meshed volume comes out negative, a void's where it comes out
/// positive. A shell that is open, whose faces disagree, or whose mesh
/// does not close is used as given; so is anything put together with
/// [`Model::add_solid`], for a builder whose shell faces in on purpose
/// (the material outside a closed face).
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if a child is not
/// a shell, or the list is empty.
pub fn make_solid(model: &mut Model, shells: &[Shape]) -> OgeomResult<Built> {
    let mut used = Vec::with_capacity(shells.len());
    for (i, shell) in shells.iter().enumerate() {
        let encloses = encloses_material(model, shell)?;
        let inside_out = encloses.is_some_and(|out| if i == 0 { !out } else { out });
        used.push(if inside_out {
            shell.reversed()
        } else {
            shell.clone()
        });
    }
    let solid = model.add_solid(&used)?;
    let mut history = History::new();
    for shell in shells {
        history.generate(shell, solid.clone());
    }
    Ok(Built::new(solid, history))
}

/// The signed volume a closed shell's mesh encloses, positive where its
/// faces face out; `None` for anything else: an open shell, one with an
/// edge both its faces walk the same way (no whole turn makes it right),
/// or one whose mesh does not close.
fn encloses_material(model: &Model, shell: &Shape) -> OgeomResult<Option<bool>> {
    /// How many of the shell's faces are probed before it is meshed.
    const PROBED: usize = 4;
    if model.node(shell).map(ogeom_topo::TShape::kind) != Some(ShapeType::Shell)
        || !is_shell_closed(model, shell)?
        || !crate::check::walked_one_way(model, shell)?.is_empty()
    {
        return Ok(None);
    }
    // Its faces agree with each other, so any one of them says which way
    // the whole shell faces: one probed off its two sides against the
    // region the shell encloses faces out of it or into it. Probing a face
    // costs that face's mesh and a few rays, where the shell's volume costs
    // every face's mesh. But a shell that passes through itself reads out
    // at some places and in at others, so no one place decides: the probes
    // settle it only where several spread over the faces agree, and where
    // they disagree or too few read, the volume the whole shell encloses
    // does.
    let tol = Tolerances::default();
    let mut read = crate::check::Readings::default();
    if let Ok(boundary) = crate::check::probe_boundary(model, shell, tol) {
        for face in ogeom_topo::explore_unique(model, shell, ShapeType::Face)?
            .iter()
            .take(PROBED)
        {
            let wanted = crate::check::AGREEING - read.count();
            if let Ok(more) = crate::check::face_readings(model, face, &boundary, wanted, tol) {
                read.out += more.out;
                read.inward += more.inward;
            }
            if read.count() >= crate::check::AGREEING {
                break;
            }
        }
    }
    if read.count() >= crate::check::AGREEING
        && let Some(inward) = read.verdict()
    {
        return Ok(Some(!inward));
    }
    Ok(enclosed_volume(model, shell)
        .map(|v| v > 0.0)
        .or_else(|| read.verdict().map(|inward| !inward)))
}

/// The volume a closed shell whose faces agree encloses, from its mesh:
/// positive where they face out of it.
fn enclosed_volume(model: &Model, shell: &Shape) -> Option<f64> {
    ogeom_mesh::triangulate(
        model,
        shell,
        ogeom_mesh::Deflection::default(),
        Tolerances::default(),
    )
    .ok()
    .filter(ogeom_topo::Triangulation::is_closed)
    .map(|mesh| mesh.volume())
}

/// Whether every edge in a shell is shared by exactly two faces.
///
/// The defining property of a closed shell, and the one that decides whether it
/// can bound a solid. An edge used once is a free boundary: the shell has a
/// hole. An edge used three or more times is non-manifold, which is legitimate
/// topology but not something that encloses a volume.
///
/// # Errors
///
/// As [`ogeom_topo::explore`].
pub fn is_shell_closed(model: &Model, shell: &Shape) -> OgeomResult<bool> {
    // Counted by *use*, not by how many distinct faces an edge belongs to. A
    // seam edge bounds one face twice (up one side of its parameter rectangle
    // and down the other), so counting faces would call every cylinder, sphere
    // and torus open, which is precisely backwards.
    let mut uses: FastMap<ogeom_topo::TShapeId, usize> = FastMap::default();
    // Read from the nodes as they are stored: what is counted is which edge
    // node each use names, which no placement changes.
    let children = |shape: &Shape| -> OgeomResult<&[Shape]> {
        match model.node(shape) {
            Some(node) => Ok(node.children()),
            None => ogeom_bail!(Dangling, "shape refers to a node not in this model"),
        }
    };
    for face in ogeom_topo::explore(model, shell, ogeom_topo::Filter::OfType(ShapeType::Face))? {
        for wire in children(&face)? {
            for edge in children(wire)? {
                *uses.entry(edge.node()).or_default() += 1;
            }
        }
    }
    // A degenerate edge (a sphere's pole, a cone's apex) has no length, so
    // there is no gap along it for a second face to close. Counting it would
    // call every sphere and every true cone open, and the thing that is
    // actually open, isn't. Only an edge used an odd number of times is
    // asked.
    let degenerate = |id: ogeom_topo::TShapeId| {
        model
            .node_by_id(id)
            .and_then(|n| n.data().as_edge())
            .is_some_and(|d| d.degenerate)
    };
    let counted = uses.keys().any(|id| !degenerate(*id));
    Ok(counted && uses.iter().all(|(id, n)| n % 2 == 0 || degenerate(*id)))
}

/// Whether a face covers the whole of its surface: it has a wire, and
/// every edge its wires walk has no length or is a seam of that surface
/// walked once each way.
///
/// Such a face has no boundary on the surface: a seam's two walks run
/// along one line in opposite senses and cancel, and an edge of no length
/// bounds nothing. On a connected closed surface (a sphere, a torus) a
/// face without boundary is the whole surface. Any other edge, or a seam
/// walked otherwise, and the answer is `false`.
///
/// # Errors
///
/// [`OgeomError::Dangling`](ogeom_core::OgeomError::Dangling) if the face,
/// a wire or an edge is not in this model.
pub fn covers_whole_surface(model: &Model, face: &Shape) -> OgeomResult<bool> {
    use ogeom_topo::Orientation;
    let Some(data) = model.node(face).and_then(|n| n.data().as_face()) else {
        ogeom_bail!(Dangling, "face is not in this model");
    };
    let wires = model.children_of(face)?;
    if wires.is_empty() {
        return Ok(false);
    }
    // Each seam's walks, forward and reversed.
    let mut walks: FastMap<ogeom_topo::TShapeId, (u32, u32)> = FastMap::default();
    for wire in &wires {
        for edge in model.children_of(wire)? {
            let Some(edge_data) = model.node(&edge).and_then(|n| n.data().as_edge()) else {
                ogeom_bail!(Dangling, "edge is not in this model");
            };
            if edge_data.degenerate {
                continue;
            }
            if !matches!(
                edge_data.pcurve_for(data.surface, edge.location()),
                Some(EdgeRepr::Seam { .. })
            ) {
                return Ok(false);
            }
            let walked = walks.entry(edge.node()).or_default();
            match edge.orientation() {
                Orientation::Forward => walked.0 += 1,
                Orientation::Reversed => walked.1 += 1,
                Orientation::Internal | Orientation::External => return Ok(false),
            }
        }
    }
    Ok(walks.values().all(|&w| w == (1, 1)))
}

/// Give every edge of `face` that has no trim on the face's surface its
/// exact one, where a closed form exists: a profile built from bare
/// curves becomes a face a sweep can keep as an end.
///
/// # Errors
///
/// [`OgeomError::Dangling`](ogeom_core::OgeomError::Dangling) if the face
/// or an edge is not in this model.
pub fn trimmed_where_bare(model: &mut Model, face: &Shape, tol: Tolerances) -> OgeomResult<()> {
    trim_bare_edges(model, face, false, tol)
}

/// As [`trimmed_where_bare`], with every edge left without a closed form
/// given a pcurve fitted by projection at the curve's own parameters: a
/// face kept whole in a solid (a sweep's end) is only valid with a trim on
/// every edge.
///
/// The fit's miss is measured, the lifted pcurve against the curve at the
/// same parameter, and the edge and its vertices widened to it. An edge
/// sitting off the surface by more than its own tolerance is not the
/// face's boundary, and refuses.
///
/// # Errors
///
/// As [`trimmed_where_bare`], and
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if an
/// edge with no closed-form trim sits off the face's surface by more than
/// its tolerance, or cannot be projected onto it.
pub fn trimmed_or_fitted_where_bare(
    model: &mut Model,
    face: &Shape,
    tol: Tolerances,
) -> OgeomResult<()> {
    trim_bare_edges(model, face, true, tol)
}

fn trim_bare_edges(model: &mut Model, face: &Shape, fit: bool, tol: Tolerances) -> OgeomResult<()> {
    let Some(data) = model.node(face).and_then(|n| n.data().as_face().cloned()) else {
        ogeom_bail!(Dangling, "face is not in this model");
    };
    let Some(surface) = model.geometry().surface(data.surface).cloned() else {
        ogeom_bail!(Dangling, "face refers to a surface not in this model");
    };
    // The edges are reached through the face, carrying its placement; the
    // surface is stated unplaced, so each curve is taken into the surface's
    // own frame before it is read in the chart.
    let unplace = face.transform(model.datums())?.inverse()?;
    for edge in ogeom_topo::explore_unique(model, face, ShapeType::Edge)? {
        let Some(edge_data) = model.node(&edge).and_then(|n| n.data().as_edge().cloned()) else {
            continue;
        };
        if edge_data
            .pcurve_for(data.surface, edge.location())
            .is_some()
        {
            continue;
        }
        let Some(EdgeRepr::Curve3d { curve, range, .. }) = edge_data.curve3d() else {
            continue;
        };
        let Some(geometry) = model.geometry().curve(*curve).cloned() else {
            continue;
        };
        use ogeom_geom::Transformable as _;
        let placed = geometry.transformed(&(unplace * edge.transform(model.datums())?), tol)?;
        // A plane's closed forms project the curve rather than test that it
        // lies there, so where a fit could take over, the projection stands
        // only if its lift meets the curve within the edge's tolerance at
        // the same parameter, the measure `check` holds it to.
        let exact = ogeom_intersect::exact_pcurve_of(&placed, &surface, tol).filter(|exact| {
            !fit || crate::pcurve_fit::lifted_gap(
                (&placed, *range),
                (exact, *range),
                &surface,
                true,
                tol,
            )
            .is_ok_and(|gap| gap <= edge_data.tolerance.get().max(tol.confusion()))
        });
        let pcurve = match exact {
            Some(exact) => exact,
            None if fit => {
                fitted_where_bare(model, &edge, &edge_data, (&placed, *range), &surface, tol)?
            }
            None => continue,
        };
        attach_pcurve(
            model,
            &edge,
            pcurve,
            data.surface,
            edge.location().clone(),
            *range,
        )?;
    }
    Ok(())
}

/// The pcurve of a bare edge with no closed form, fitted at the curve's
/// own parameters, the edge and its vertices widened to the fit's
/// measured miss.
fn fitted_where_bare(
    model: &mut Model,
    edge: &Shape,
    edge_data: &EdgeData,
    (curve, range): (&Curve, (f64, f64)),
    surface: &SurfaceGeometry,
    tol: Tolerances,
) -> OgeomResult<ogeom_geom::PlanarCurve> {
    let reach = edge_data.tolerance.get().max(tol.confusion());
    let fitted =
        match crate::pcurve_fit::fit_projected_pcurve_within(curve, range, surface, reach, tol)? {
            crate::pcurve_fit::CappedFit::Fitted((fitted, ..)) => fitted,
            crate::pcurve_fit::CappedFit::TooFar(off) => ogeom_bail!(
                Construction,
                "an edge of the face sits {off:.2e} off its surface, outside the edge's \
                 tolerance of {reach:.2e}, so it has no trim there"
            ),
        };
    // Held to the same parameter, the stricter of the two measures `check`
    // takes, so the width holds whether or not the edge claims it.
    let gap = crate::pcurve_fit::lifted_gap((curve, range), (&fitted, range), surface, true, tol)?;
    if gap > tol.confusion() {
        // The edge owns the miss, and so must the vertices that bound it:
        // the containment rule.
        let widened = ogeom_core::Tolerance::new(gap + tol.confusion())?;
        model.widen(edge, widened)?;
        if let Some((a, b)) = edge_vertices(model, edge)? {
            model.widen(&a, widened)?;
            model.widen(&b, widened)?;
        }
    }
    Ok(fitted)
}

/// Attach a pcurve to an edge, describing it in a surface's parameter space.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if `edge` is not
/// an edge; [`OgeomError::Dangling`](ogeom_core::OgeomError::Dangling) if it does not
/// resolve.
pub fn attach_pcurve(
    model: &mut Model,
    edge: &Shape,
    pcurve: ogeom_geom::PlanarCurve,
    surface: ogeom_topo::SurfaceId,
    location: Location,
    range: (f64, f64),
) -> OgeomResult<()> {
    if model.kind_of(edge)? != ShapeType::Edge {
        ogeom_bail!(Construction, "pcurves attach to edges");
    }
    let id = model.geometry_mut().add_pcurve(pcurve);
    let Some(node) = model.node_mut(edge) else {
        ogeom_bail!(Dangling, "edge is not in this model");
    };
    let ogeom_topo::NodeData::Edge(data) = node.data_mut() else {
        ogeom_bail!(Construction, "edge node holds no edge data");
    };
    data.add(EdgeRepr::PCurve {
        curve: id,
        surface,
        location,
        range,
    });
    Ok(())
}

/// Attach a seam representation to an edge: one pcurve per side of a closed
/// surface's join.
///
/// The counterpart of [`attach_pcurve`] for the edge that bounds one face
/// twice: up one side of the parameter rectangle and down the other. Which
/// pcurve applies to an occurrence is decided by that occurrence's
/// orientation, which is the only thing distinguishing the two.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if `edge` is not
/// an edge; [`OgeomError::Dangling`](ogeom_core::OgeomError::Dangling) if it is not in
/// this model.
pub fn attach_seam(
    model: &mut Model,
    edge: &Shape,
    forward: ogeom_geom::PlanarCurve,
    reversed: ogeom_geom::PlanarCurve,
    surface: ogeom_topo::SurfaceId,
    location: Location,
    range: (f64, f64),
) -> OgeomResult<()> {
    if model.kind_of(edge)? != ShapeType::Edge {
        ogeom_bail!(Construction, "seams attach to edges");
    }
    let forward = model.geometry_mut().add_pcurve(forward);
    let reversed = model.geometry_mut().add_pcurve(reversed);
    let Some(node) = model.node_mut(edge) else {
        ogeom_bail!(Dangling, "edge is not in this model");
    };
    let ogeom_topo::NodeData::Edge(data) = node.data_mut() else {
        ogeom_bail!(Construction, "edge node holds no edge data");
    };
    data.add(EdgeRepr::Seam {
        forward,
        reversed,
        surface,
        location,
        range,
    });
    Ok(())
}

/// Whether a circle is a *parallel* of a revolved surface (its axis the
/// revolution axis, its centre on it) rather than a circle that merely
/// lies on the surface.
fn circle_is_parallel_of(
    circle: ogeom_math::Circle,
    surface: &ogeom_geom::SurfaceGeometry,
    axis_z: ogeom_math::Vector,
    tol: Tolerances,
) -> bool {
    if circle.frame().z().vector().cross(axis_z).magnitude() > tol.angular() {
        return false;
    }
    match surface_axis_origin(surface) {
        Some(origin) => {
            let off = circle.centre() - origin;
            let radial = off - axis_z * off.dot(axis_z);
            radial.magnitude() <= tol.confusion() * 10.0
        }
        None => true,
    }
}

/// Whether two closed ring edges are parallels of `surface`, so that
/// [`make_revolution_band`] can build a band between them. A reader asks
/// this first: a periodic face bounded by two closed circles that are *not*
/// parallels (a button head's rims, square to the screw on a sphere whose
/// chart runs along z) is a legitimate face on its own bounds, not a band
/// missing its seam, and deserves no warning.
///
/// # Errors
///
/// [`OgeomError::Dangling`](ogeom_core::OgeomError::Dangling) if an edge is
/// not in the model.
pub fn rings_are_parallels(
    model: &Model,
    surface: &ogeom_geom::SurfaceGeometry,
    rings: &[&Shape],
    tol: Tolerances,
) -> OgeomResult<bool> {
    let Some(axis_z) = surface_iso_axis(surface) else {
        return Ok(false);
    };
    for ring in rings {
        let Some(node) = model.node(ring) else {
            ogeom_bail!(Dangling, "edge is not in this model");
        };
        let Some(data) = node.data().as_edge() else {
            return Ok(false);
        };
        if data.degenerate {
            continue;
        }
        let Some(EdgeRepr::Curve3d { curve, .. }) = data.curve3d() else {
            return Ok(false);
        };
        let Some(ogeom_geom::Curve::Circle(c)) = model.geometry().curve(*curve) else {
            return Ok(false);
        };
        if !circle_is_parallel_of(c.circle(), surface, axis_z, tol) {
            return Ok(false);
        }
    }
    Ok(true)
}

/// A point on the revolution axis of an analytic surface, where one is
/// stated: the frame origin of a cylinder, cone or torus, a sphere's centre.
fn surface_axis_origin(surface: &ogeom_geom::SurfaceGeometry) -> Option<Point> {
    use ogeom_geom::SurfaceGeometry as S;
    match surface {
        S::Cylinder(c) => Some(c.cylinder().frame().origin()),
        S::Cone(c) => Some(c.cone().frame().origin()),
        S::Sphere(s) => Some(s.sphere().centre()),
        S::Torus(t) => Some(t.torus().frame().origin()),
        _ => None,
    }
}

/// Build the face of a revolution band: two closed rings joined by a
/// synthesised seam, pcurves attached window-coherently.
///
/// The one construction of a band's seam, for every caller. The chart walk
/// decides everything: the bottom ring is
/// traversed forward, and where its chart line *ends* is where the first
/// seam column stands; the top ring's occurrence direction is chosen so its
/// walk starts there (rings winding the same way traverse opposite, rings
/// winding opposite traverse alike), and the seam's two pcurves are assigned
/// to match which occurrence the triangulator will hand them to. The face is
/// built on a fresh copy of the surface, so no annotation left by an earlier
/// use of the same rings can apply.
///
/// One ring may be *degenerate* (an edge with no curve, both ends the same
/// vertex, flagged as such), standing for an apex or a pole: a rim of no
/// length that still bounds the face in parameter space. It takes the row
/// the collapsed point sits on and traverses the chart opposite the real
/// ring, exactly as native cones and spheres bound their tips.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if a ring is
/// neither a circle nor degenerate, both rings are degenerate, the rings are
/// not rings of this surface, their vertices stand on different columns of
/// it, or the surface's iso-curve has no closed form. A refused band leaves
/// the rings' edges as they were.
pub fn make_revolution_band(
    model: &mut Model,
    surface: &ogeom_geom::SurfaceGeometry,
    ring_lo: &Shape,
    ring_hi: &Shape,
    tol: Tolerances,
) -> OgeomResult<Shape> {
    revolution_band(model, surface, ring_lo, ring_hi, None, tol)
}

/// [`make_revolution_band`] on a curve revolved a whole turn about an
/// axis ([`RevolutionSurface`](ogeom_geom::RevolutionSurface)), its rings
/// the parallels through the profile's two ends.
///
/// The seam is the profile turned to the column of the first ring's
/// vertex, on the profile's own parameter, and each ring's row is the
/// profile's end it is the parallel of. The rings must be the parallels
/// through the profile's ends, each a circle about the axis.
///
/// # Errors
///
/// As [`make_revolution_band`], and
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the
/// surface is not a whole turn of a revolution, or a ring is not the
/// parallel through an end of its profile.
pub fn make_profile_revolution_band(
    model: &mut Model,
    surface: &ogeom_geom::SurfaceGeometry,
    ring_lo: &Shape,
    ring_hi: &Shape,
    tol: Tolerances,
) -> OgeomResult<Shape> {
    use ogeom_geom::Surface as _;
    let SurfaceGeometry::Revolution(revolved) = surface else {
        ogeom_bail!(
            Construction,
            "a profile's band asked of a surface no profile sweeps"
        );
    };
    let ((ua, ub), _) = surface.domain();
    if !surface.is_periodic_u() || ((ub - ua) - core::f64::consts::TAU).abs() > tol.angular() {
        ogeom_bail!(
            Construction,
            "a profile's band asked of less than a whole turn"
        );
    }
    revolution_band(model, surface, ring_lo, ring_hi, Some(revolved), tol)
}

/// Where a point on a whole turn of a revolved profile stands in its
/// chart: its angle about the axis from the profile's own half-plane, and
/// the profile's parameter at the end of it the point turns from. `None`
/// for a point on the parallel of neither end.
fn profile_chart_of(
    revolved: &ogeom_geom::RevolutionSurface,
    p: Point,
    tol: Tolerances,
) -> Option<ogeom_math::Point2> {
    use ogeom_geom::Curve3d as _;
    let axis = revolved.axis();
    let curve = revolved.curve();
    let (lo, hi) = curve.domain();
    // A point's height along the axis and its radius off it.
    let place = |q: Point| {
        let off = q - axis.location;
        let height = off.dot(axis.direction.vector());
        let radial = off - axis.direction.vector() * height;
        (height, radial)
    };
    let (height, radial) = place(p);
    let mut found = None;
    for end in [lo, hi] {
        let q = curve.point_at(end, tol).ok()?;
        let (h, r) = place(q);
        if (h - height).abs() <= tol.confusion()
            && (r.magnitude() - radial.magnitude()).abs() <= tol.confusion()
        {
            let z = axis.direction.vector();
            let x = r / r.magnitude();
            let y = z.cross(x);
            let angle = radial.dot(y).atan2(radial.dot(x));
            found = Some(ogeom_math::Point2::new(
                angle.rem_euclid(core::f64::consts::TAU),
                end,
            ));
            break;
        }
    }
    found
}

/// [`make_revolution_band`], or [`make_profile_revolution_band`] where
/// `profile` is the surface's own revolution.
fn revolution_band(
    model: &mut Model,
    surface: &ogeom_geom::SurfaceGeometry,
    ring_lo: &Shape,
    ring_hi: &Shape,
    profile: Option<&ogeom_geom::RevolutionSurface>,
    tol: Tolerances,
) -> OgeomResult<Shape> {
    use ogeom_geom::Surface as _;

    let ((ua_dom, ub_dom), _) = surface.domain();
    let span = ub_dom - ua_dom;
    let axis_z = surface_iso_axis(surface)
        .or_else(|| profile.map(|r| r.axis().direction.vector()))
        .ok_or_else(|| {
            ogeom_core::ogeom_err!(Construction, "the surface has no revolution axis")
        })?;
    // The chart's closed-form inversion: the analytic kinds', or the
    // profile's at its ends.
    let chart_of = |p: Point| -> Option<ogeom_math::Point2> {
        match profile {
            Some(revolved) => profile_chart_of(revolved, p, tol),
            None => analytic_chart_of(surface, p),
        }
    };

    // Each ring's curve, range, winding, vertex, and chart row.
    struct Ring {
        edge: Shape,
        vertex: Shape,
        crange: (f64, f64),
        winding: f64,
        row: f64,
        degenerate: bool,
    }
    let mut rings = Vec::new();
    for used in [ring_lo, ring_hi] {
        // The caller's edges may carry orientation from their old wire uses;
        // the band is built on the curves' own directions, and the walk
        // chooses each occurrence's orientation itself.
        let edge = &if used.orientation() == ogeom_topo::Orientation::Reversed {
            used.reversed()
        } else {
            used.clone()
        };
        let Some((vertex, other)) = edge_vertices(model, edge)? else {
            ogeom_bail!(Construction, "a band ring has no vertex");
        };
        if !vertex.is_same(&other) {
            ogeom_bail!(Construction, "a band ring is not closed");
        }
        let at = {
            let Some(node) = model.node(&vertex) else {
                ogeom_bail!(Dangling, "vertex is not in this model");
            };
            let Some(data) = node.data().as_vertex() else {
                ogeom_bail!(Construction, "vertex node holds no vertex data");
            };
            data.point
        };
        let degenerate = {
            let Some(node) = model.node(edge) else {
                ogeom_bail!(Dangling, "edge is not in this model");
            };
            let Some(data) = node.data().as_edge() else {
                ogeom_bail!(Construction, "edge node holds no edge data");
            };
            data.degenerate
        };
        if degenerate {
            // An apex or a pole: no curve to read, no winding of its own.
            // Its row comes from the collapsed point, which sits on the
            // axis, where iterative projection has no nearest angle, so only
            // the closed-form inversion can place it.
            let Some(uv) = chart_of(at) else {
                ogeom_bail!(
                    Construction,
                    "a degenerate ring's row cannot be found on this surface"
                );
            };
            rings.push(Ring {
                edge: edge.clone(),
                vertex,
                crange: (0.0, span),
                winding: 0.0,
                row: uv.y,
                degenerate: true,
            });
            continue;
        }
        let (curve, crange) = {
            let Some(node) = model.node(edge) else {
                ogeom_bail!(Dangling, "edge is not in this model");
            };
            let Some(data) = node.data().as_edge() else {
                ogeom_bail!(Construction, "edge node holds no edge data");
            };
            let Some(EdgeRepr::Curve3d { curve, range, .. }) = data.curve3d() else {
                ogeom_bail!(Construction, "a band ring has no curve");
            };
            let Some(geometry) = model.geometry().curve(*curve) else {
                ogeom_bail!(Dangling, "curve is not in this model");
            };
            (geometry.clone(), *range)
        };
        let ogeom_geom::Curve::Circle(c) = &curve else {
            ogeom_bail!(Construction, "a band ring is not a circle");
        };
        // A ring of *this* surface is a parallel: its axis the revolution
        // axis, its centre on it. A circle merely lying on the surface
        // (a button head's rim, cut square to the screw while the sphere's
        // chart runs along z) is a closed loop in the chart, not a row of
        // it, and building a band on it would hand every rim a latitude
        // line it never follows. Refused here, so a reader falls through to
        // the face's own bounds, which triangulate as the nested loops they
        // are.
        let off_axis = profile.is_some_and(|r| {
            let off = c.circle().centre() - r.axis().location;
            (off - axis_z * off.dot(axis_z)).magnitude() > tol.confusion() * 10.0
        });
        if off_axis || !circle_is_parallel_of(c.circle(), surface, axis_z, tol) {
            ogeom_bail!(
                Construction,
                "a band ring is a circle on the surface but not a parallel \
                 of it; its axis or centre is off the revolution axis"
            );
        }
        // A circle turns about its frame's spin, and back about it where
        // the curve is reversed.
        let turned = if c.is_reversed() { -1.0 } else { 1.0 };
        let winding = turned * spin(&c.circle().frame()).dot(axis_z).signum();
        let row = match chart_of(at) {
            Some(uv) => uv.y,
            None => {
                crate::measure::project_on_surface(surface, at, 32, tol)?
                    .parameters
                    .1
            }
        };
        rings.push(Ring {
            edge: edge.clone(),
            vertex,
            crange,
            winding,
            row,
            degenerate: false,
        });
    }
    // A degenerate ring bounds nothing by itself, and the chart anchor must
    // come from a rim that has an angle: the real ring goes first, and the
    // degenerate one traverses opposite it.
    if rings[0].degenerate && rings[1].degenerate {
        ogeom_bail!(Construction, "both band rings are degenerate");
    }
    if rings[0].degenerate {
        rings.swap(0, 1);
    }
    if rings[1].degenerate {
        rings[1].winding = -rings[0].winding;
    }
    // A row on the seam of a chart periodic in `v` (a torus's outer
    // equator) reads at either end of it, and which one rounding picked
    // decides whether the band runs a quarter of the tube or the other
    // three. Such a row is taken at the end nearer the other ring's row.
    if surface.is_periodic_v() {
        let (_, (va, vb)) = surface.domain();
        let on_seam =
            |row: f64| (row - va).abs() <= tol.angular() || (vb - row).abs() <= tol.angular();
        for (this, other) in [(0, 1), (1, 0)] {
            if on_seam(rings[this].row) && !on_seam(rings[other].row) {
                let toward = rings[other].row;
                rings[this].row = if (toward - va).abs() <= (vb - toward).abs() {
                    va
                } else {
                    vb
                };
            }
        }
    }
    let anchor = {
        let Some(node) = model.node(&rings[0].vertex) else {
            ogeom_bail!(Dangling, "vertex is not in this model");
        };
        let Some(data) = node.data().as_vertex() else {
            ogeom_bail!(Construction, "vertex node holds no vertex data");
        };
        data.point
    };
    let ua = match chart_of(anchor) {
        Some(uv) => uv.x,
        None => {
            crate::measure::project_on_surface(surface, anchor, 32, tol)?
                .parameters
                .0
        }
    };
    if !rings[1].degenerate && circle_crossing_of(model, &rings[1].edge, anchor, tol)?.is_some() {
        ogeom_bail!(
            Construction,
            "the band rings' vertices stand on different columns, which no \
             seam joins; split_ring_on_column gives one rim a vertex on the \
             other's column"
        );
    }

    // The seam runs along the surface's own iso-curve at the anchor angle,
    // parameterized by `v`, built along increasing `v`.
    let (va, vb) = (rings[0].row, rings[1].row);
    let turned = match profile {
        Some(revolved) => {
            use ogeom_geom::Transformable as _;
            let turn = ogeom_math::Transform::rotation(revolved.axis(), ua);
            Some(revolved.curve().clone().transformed(&turn, tol)?)
        }
        None => seam_iso_u_curve(surface, ua, tol),
    };
    let Some(seam_curve) = turned else {
        ogeom_bail!(
            Construction,
            "the surface's iso-curve has no closed form; no seam can be built"
        );
    };
    let (range, from, to, downward) = if va <= vb {
        (
            (va, vb),
            rings[0].vertex.clone(),
            rings[1].vertex.clone(),
            false,
        )
    } else {
        (
            (vb, va),
            rings[1].vertex.clone(),
            rings[0].vertex.clone(),
            true,
        )
    };
    // The edge lives on the curve's own parameterization; the chart rows map
    // onto it linearly, which is exactly the rescale a pcurve range states.
    let curve_range = (
        iso_curve_parameter_at(surface, range.0),
        iso_curve_parameter_at(surface, range.1),
    );
    let seam = make_edge_between(model, seam_curve, curve_range, &from, &to, tol)?.shape;

    // Nothing above touches the rings, so a refusal (ring vertices on
    // different columns, which no seam joins) leaves them as they came,
    // with no pcurve on a surface no face uses.
    let surface_id = model.geometry_mut().add_surface(surface.clone());

    // Window-coherent ring pcurves: u(t) spans [ua, ua + span] whichever way
    // each ring winds.
    for ring in &rings {
        let u_start = if ring.winding > 0.0 { ua } else { ua + span };
        let origin =
            ogeom_math::Point2::new(ring.winding.mul_add(-ring.crange.0, u_start), ring.row);
        let pcurve: ogeom_geom::PlanarCurve = ogeom_geom::Line2d::over(
            ogeom_math::Axis2::new(
                origin,
                ogeom_math::Direction2::new(ogeom_math::Vector2::new(ring.winding, 0.0), tol)?,
            ),
            ring.crange.0,
            ring.crange.1,
        )?
        .into();
        attach_pcurve(
            model,
            &ring.edge,
            pcurve,
            surface_id,
            Location::identity(),
            ring.crange,
        )?;
    }

    // The walk closes only if the top ring's traversal starts where the
    // bottom's ends, and the seam sides sit at the columns the walk visits.
    let bottom_end = if rings[0].winding > 0.0 {
        ua + span
    } else {
        ua
    };
    let other_col = if rings[0].winding > 0.0 {
        ua
    } else {
        ua + span
    };
    let hi_reversed = (rings[1].winding - rings[0].winding).abs() < 0.5;
    let column = |u: f64| -> OgeomResult<ogeom_geom::PlanarCurve> {
        Ok(ogeom_geom::Line2d::over(
            ogeom_math::Axis2::new(
                ogeom_math::Point2::new(u, 0.0),
                ogeom_math::Direction2::new(ogeom_math::Vector2::new(0.0, 1.0), tol)?,
            ),
            range.0 - 1.0,
            range.1 + 1.0,
        )?
        .into())
    };
    // The walk runs counter-clockwise in the chart where the first ring
    // turns with `u` along the lower row, or against it along the upper; the
    // other way round, the ring is walked back, so the face keeps its
    // material on the left about the surface's normal.
    let counter_clockwise = (rings[0].winding > 0.0) == (va < vb);
    // The first seam occurrence in the wire is `up`; the triangulator hands
    // a Forward occurrence the `forward` pcurve. `up` is Forward exactly
    // when the seam was built upward, and it stands where the first ring's
    // walk ends, or where it starts when the ring is walked back.
    let (forward_col, reversed_col) = if downward == counter_clockwise {
        (other_col, bottom_end)
    } else {
        (bottom_end, other_col)
    };
    attach_seam(
        model,
        &seam,
        column(forward_col)?,
        column(reversed_col)?,
        surface_id,
        Location::identity(),
        range,
    )?;

    let up = if downward {
        seam.reversed()
    } else {
        seam.clone()
    };
    let top = if hi_reversed {
        rings[1].edge.reversed()
    } else {
        rings[1].edge.clone()
    };
    let ring = if counter_clockwise {
        vec![rings[0].edge.clone(), up.clone(), top, up.reversed()]
    } else {
        vec![
            rings[0].edge.reversed(),
            up.clone(),
            top.reversed(),
            up.reversed(),
        ]
    };
    let wire = make_wire(model, &ring, tol)?.shape;
    Ok(make_face_on(model, surface_id, &[wire], tol)?.shape)
}

/// A circle edge's curve and range, the parameter where it crosses a
/// column, and the point there.
type RingCrossing = (Curve, (f64, f64), f64, Point);

/// Where a closed circle edge crosses the column of `anchor`: the circle's
/// parameter there, inside the edge's range, and the point. `None` where
/// the edge's own vertex already stands on that column.
fn circle_crossing_of(
    model: &Model,
    ring: &Shape,
    anchor: Point,
    tol: Tolerances,
) -> OgeomResult<Option<RingCrossing>> {
    let Some(data) = model.node(ring).and_then(|n| n.data().as_edge()) else {
        ogeom_bail!(Dangling, "edge is not in this model");
    };
    let Some(EdgeRepr::Curve3d { curve, range, .. }) = data.curve3d() else {
        ogeom_bail!(Construction, "a band ring has no curve");
    };
    let range = *range;
    let Some(geometry) = model.geometry().curve(*curve) else {
        ogeom_bail!(Dangling, "curve is not in this model");
    };
    let Curve::Circle(c) = geometry else {
        ogeom_bail!(Construction, "a band ring is not a circle");
    };
    let circle = c.circle();
    let local = circle.frame().to_local(anchor);
    if local.x.hypot(local.y) <= tol.confusion() {
        ogeom_bail!(
            Construction,
            "the anchor stands on the ring's axis, where every column meets"
        );
    }
    let angle = local.y.atan2(local.x);
    let along = if c.is_reversed() { -angle } else { angle };
    let period = core::f64::consts::TAU;
    let t = range.0 + (along - range.0).rem_euclid(period);
    let point = geometry.point_at(t, tol)?;
    let start = geometry.point_at(range.0, tol)?;
    let end = geometry.point_at(range.1, tol)?;
    if point.distance(start) <= tol.confusion() || point.distance(end) <= tol.confusion() {
        return Ok(None);
    }
    Ok(Some((geometry.clone(), range, t, point)))
}

/// Split a closed circle rim of a revolution band where it crosses the
/// column of `anchor`, the vertex of the band's other rim, so both rims
/// have a vertex on one column and a seam can join them.
///
/// The rim must be a parallel of `surface` (its axis the revolution axis,
/// its centre on it), as [`make_revolution_band`] requires. It becomes two
/// edges on its own curve: the first from its vertex to a new vertex on the
/// anchor's column, the second back. The result is a closed wire of the two
/// in the curve's direction starting on the anchor's column, the second
/// piece first, ready to be a ring of [`make_band_of_rings`]. The history
/// records the rim as modified into both pieces.
///
/// The pieces carry no pcurves, and every face holding the rim must take
/// both pieces in its place, or the shell opens along the rim: a caller
/// splitting a rim that bounds a face already built rebuilds that face.
///
/// Returns `None` when the rim's vertex already stands on the anchor's
/// column.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the
/// rim is not a closed circle that is a parallel of `surface`, or the
/// anchor stands on its axis.
pub fn split_ring_on_column(
    model: &mut Model,
    surface: &SurfaceGeometry,
    ring: &Shape,
    anchor: &Shape,
    tol: Tolerances,
) -> OgeomResult<Option<Built>> {
    let rim = ring.oriented(ogeom_topo::Orientation::Forward);
    let Some((vertex, other)) = edge_vertices(model, &rim)? else {
        ogeom_bail!(Construction, "a band ring has no vertex");
    };
    if !vertex.is_same(&other) {
        ogeom_bail!(Construction, "a band ring is not closed");
    }
    if !rings_are_parallels(model, surface, &[&rim], tol)? {
        ogeom_bail!(
            Construction,
            "a band ring is not a parallel of the surface it bounds"
        );
    }
    let anchor_point = {
        let Some(data) = model.node(anchor).and_then(|n| n.data().as_vertex()) else {
            ogeom_bail!(Construction, "the anchor is not a vertex of this model");
        };
        anchor.transform(model.datums())?.apply(data.point)
    };
    let Some((curve, range, t, point)) = circle_crossing_of(model, &rim, anchor_point, tol)? else {
        return Ok(None);
    };
    let stated = model
        .node(&rim)
        .and_then(|n| n.data().as_edge())
        .map(|d| d.tolerance);
    let column = make_vertex(model, point).shape;
    let head = make_edge_between(model, curve.clone(), (range.0, t), &vertex, &column, tol)?.shape;
    let tail = make_edge_between(model, curve, (t, range.1), &column, &vertex, tol)?.shape;
    if let Some(stated) = stated {
        model.widen(&head, stated)?;
        model.widen(&tail, stated)?;
    }
    let wire = make_wire(model, &[tail.clone(), head.clone()], tol)?.shape;
    let mut history = History::new();
    history.modify(&rim, head);
    history.modify(&rim, tail);
    Ok(Some(Built::new(wire, history)))
}

/// Build the face of a band between two closed rings of edges on a periodic
/// analytic surface, joined by a synthesised seam.
///
/// [`make_revolution_band`] for rings of any number of edges, each edge a
/// curve with a closed-form pcurve on the surface: a notched rim of arcs and
/// rulings, not only a whole circle. Each ring is the edge occurrences of a
/// closed wire in order, the first starting on the seam's column; a ring of
/// one degenerate edge stands for a pole or an apex. The first ring is
/// walked as given and fixes the seam's column and which way round the band
/// runs; the second is walked back the other way, reversed if it was given
/// winding the same way, and the seam joins the two start vertices along the
/// surface's iso-curve. The face is built on a fresh copy of the surface.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the
/// surface is not a periodic analytic surface, a ring does not go once round
/// it, the rings' start vertices do not share a column, an edge has no
/// closed-form pcurve, or the iso-curve has no closed form.
pub fn make_band_of_rings(
    model: &mut Model,
    surface: &SurfaceGeometry,
    ring_a: &[Shape],
    ring_b: &[Shape],
    tol: Tolerances,
) -> OgeomResult<Shape> {
    use ogeom_geom::Curve2d as _;
    use ogeom_geom::Surface as _;
    use ogeom_math::{Axis2, Direction2, Point2, Transform2, Vector2};

    let ((ua_dom, ub_dom), _) = surface.domain();
    let span = ub_dom - ua_dom;
    if !surface.is_periodic_u() || span <= 0.0 || ring_a.is_empty() || ring_b.is_empty() {
        ogeom_bail!(
            Construction,
            "a band needs two rings on a surface closing in u"
        );
    }
    let point_of = |model: &Model, vertex: &Shape| -> OgeomResult<Point> {
        let Some(data) = model.node(vertex).and_then(|n| n.data().as_vertex()) else {
            ogeom_bail!(Dangling, "vertex is not in this model");
        };
        Ok(vertex.transform(model.datums())?.apply(data.point))
    };
    let start_of = |model: &Model, ring: &[Shape]| -> OgeomResult<Shape> {
        let Some((start, _)) = edge_vertices(model, &ring[0])? else {
            ogeom_bail!(Construction, "a band ring's edge has no vertex");
        };
        Ok(start)
    };
    let (start_a, start_b) = (start_of(model, ring_a)?, start_of(model, ring_b)?);
    // The face's copy of the surface is turned about its axis so the seam
    // stands at the start of its window: the band then lies in the
    // surface's own period, as a primitive's does.
    let surface = {
        let Some(at) = analytic_chart_of(surface, point_of(model, &start_a)?) else {
            ogeom_bail!(Construction, "the band's surface has no closed-form chart");
        };
        let (Some(origin), Some(axis)) = (surface_axis_origin(surface), surface_iso_axis(surface))
        else {
            ogeom_bail!(Construction, "the band's surface has no axis");
        };
        let turn = ogeom_math::Transform::rotation(
            ogeom_math::Axis::new(origin, ogeom_math::Direction::new(axis, tol)?),
            at.x - ua_dom,
        );
        ogeom_geom::Transformable::transformed(surface, &turn, tol)?
    };
    let surface = &surface;
    let (Some(chart_a), Some(chart_b)) = (
        analytic_chart_of(surface, point_of(model, &start_a)?),
        analytic_chart_of(surface, point_of(model, &start_b)?),
    ) else {
        ogeom_bail!(Construction, "the band's surface has no closed-form chart");
    };
    // At the window's start, give or take a rounding either side of it.
    let column = if (chart_a.x - ub_dom).abs() <= 1e-9 * span {
        ua_dom
    } else {
        chart_a.x
    };
    let degenerate_b = ring_b.len() == 1
        && model
            .node(&ring_b[0])
            .and_then(|n| n.data().as_edge())
            .is_some_and(|d| d.degenerate);
    // A pole's column is any column; a rim's start must stand on the seam.
    let apart = (chart_b.x - column).rem_euclid(span);
    if !degenerate_b && apart.min(span - apart) > 1e-6 * span {
        ogeom_bail!(
            Construction,
            "the band rings' start vertices are not on one column"
        );
    }

    // A ring's pcurves, each moved by whole periods to continue the walk
    // from `from`, and where the walk ends.
    type Walk = (Vec<(Shape, ogeom_geom::PlanarCurve, (f64, f64))>, Point2);
    let walk = |model: &Model, ring: &[Shape], from: Point2| -> OgeomResult<Walk> {
        let mut out = Vec::with_capacity(ring.len());
        let mut at = from;
        for occurrence in ring {
            let Some(data) = model.node(occurrence).and_then(|n| n.data().as_edge()) else {
                ogeom_bail!(Dangling, "edge is not in this model");
            };
            let Some(EdgeRepr::Curve3d { curve, range, .. }) = data.curve3d() else {
                ogeom_bail!(Construction, "a band ring's edge has no curve");
            };
            let Some(geometry) = model.geometry().curve(*curve) else {
                ogeom_bail!(Dangling, "curve is not in this model");
            };
            let geometry = ogeom_geom::Transformable::transformed(
                geometry,
                &occurrence.transform(model.datums())?,
                tol,
            )?;
            let Some(pcurve) = ogeom_intersect::exact_pcurve_over(&geometry, *range, surface, tol)
            else {
                ogeom_bail!(Construction, "a band ring's edge has no closed-form pcurve");
            };
            // A line pcurve's stated domain must cover the edge's range,
            // which for a circle anchored anywhere runs past one period.
            let pcurve = match &pcurve {
                ogeom_geom::PlanarCurve::Line(l) => {
                    let (lo, hi) = (l.domain().0.min(range.0), l.domain().1.max(range.1));
                    ogeom_geom::Line2d::over(l.axis(), lo, hi)?.into()
                }
                _ => pcurve,
            };
            let reversed = occurrence.orientation() == ogeom_topo::Orientation::Reversed;
            let (t0, t1) = if reversed {
                (range.1, range.0)
            } else {
                (range.0, range.1)
            };
            let start = pcurve.point_at(t0, tol)?;
            let turns = ((at.x - start.x) / span).round();
            let pcurve = pcurve.transformed(
                &Transform2::translation(Vector2::new(turns * span, 0.0)),
                tol,
            )?;
            at = pcurve.point_at(t1, tol)?;
            out.push((occurrence.clone(), pcurve, *range));
        }
        Ok((out, at))
    };
    let once_round = |from: f64, to: f64| ((to - from).abs() - span).abs() <= 1e-6 * span;

    // The band stands in the period from the column up, the surface's own
    // window's way round: a ring winding backward starts a period on.
    let (mut walked_a, mut end_a) = walk(model, ring_a, Point2::new(column, chart_a.y))?;
    let column = if end_a.x < column {
        (walked_a, end_a) = walk(model, ring_a, Point2::new(column + span, chart_a.y))?;
        column + span
    } else {
        column
    };
    if !once_round(column, end_a.x) {
        ogeom_bail!(Construction, "the first band ring does not go once round");
    }
    let end_column = end_a.x;

    // The second ring runs back from where the first ended.
    let mut ring_b: Vec<Shape> = ring_b.to_vec();
    let walked_b = if degenerate_b {
        let toward = (column - end_column).signum();
        let pole: ogeom_geom::PlanarCurve = ogeom_geom::Line2d::over(
            Axis2::new(
                Point2::new(end_column, chart_b.y),
                Direction2::new(Vector2::new(toward, 0.0), tol)?,
            ),
            0.0,
            span,
        )?
        .into();
        vec![(ring_b[0].clone(), pole, (0.0, span))]
    } else {
        let mut walked = walk(model, &ring_b, Point2::new(end_column, chart_b.y))?;
        if !once_round(end_column, walked.1.x) {
            ogeom_bail!(Construction, "the second band ring does not go once round");
        }
        if (walked.1.x - column).abs() > 1e-6 * span {
            ring_b = ring_b.iter().rev().map(Shape::reversed).collect();
            walked = walk(model, &ring_b, Point2::new(end_column, chart_b.y))?;
        }
        walked.0
    };

    let surface_id = model.geometry_mut().add_surface(surface.clone());
    for (edge, pcurve, range) in walked_a.into_iter().chain(walked_b) {
        attach_pcurve(
            model,
            &edge,
            pcurve,
            surface_id,
            Location::identity(),
            range,
        )?;
    }

    // The seam along the iso-curve at the column, built upward in `v`; the
    // occurrence climbing from the first ring to the second stands on the
    // column the first ring's walk ended at.
    let (va, vb) = (chart_a.y, chart_b.y);
    let Some(seam_curve) = seam_iso_u_curve(surface, column, tol) else {
        ogeom_bail!(
            Construction,
            "the surface's iso-curve has no closed form; no seam can be built"
        );
    };
    let (rows, from, to, downward) = if va <= vb {
        ((va, vb), start_a.clone(), start_b.clone(), false)
    } else {
        ((vb, va), start_b.clone(), start_a.clone(), true)
    };
    let curve_range = (
        iso_curve_parameter_at(surface, rows.0),
        iso_curve_parameter_at(surface, rows.1),
    );
    let seam = make_edge_between(model, seam_curve, curve_range, &from, &to, tol)?.shape;
    let side = |u: f64| -> OgeomResult<ogeom_geom::PlanarCurve> {
        Ok(ogeom_geom::Line2d::over(
            Axis2::new(
                Point2::new(u, 0.0),
                Direction2::new(Vector2::new(0.0, 1.0), tol)?,
            ),
            rows.0 - 1.0,
            rows.1 + 1.0,
        )?
        .into())
    };
    // The rings run as given, so the face's neighbours still meet them
    // the other way; which way the loop then goes round the chart says
    // which side the face faces. Clockwise, it faces against the surface,
    // and stores the loop walked back so it keeps the face on its left in
    // the chart: the climbing occurrence then stands where the walk
    // starts, on the first ring's column.
    let clockwise = (end_column - column) * (vb - va) < 0.0;
    let (forward_side, reversed_side) = if downward != clockwise {
        (column, end_column)
    } else {
        (end_column, column)
    };
    attach_seam(
        model,
        &seam,
        side(forward_side)?,
        side(reversed_side)?,
        surface_id,
        Location::identity(),
        rows,
    )?;
    let up = if downward {
        seam.reversed()
    } else {
        seam.clone()
    };
    let mut edges: Vec<Shape> = ring_a.to_vec();
    edges.push(up.clone());
    edges.extend(ring_b);
    edges.push(up.reversed());
    if clockwise {
        edges = edges.iter().rev().map(Shape::reversed).collect();
    }
    let wire = make_wire(model, &edges, tol)?.shape;
    let face = make_face_on(model, surface_id, &[wire], tol)?.shape;
    Ok(if clockwise { face.reversed() } else { face })
}

/// Build the face of a band between two closed rings that need not be
/// circles, each with a caller-supplied chart image.
///
/// [`make_revolution_band`]'s general sibling: where the band derives each
/// circle's row and pcurve itself, this takes the rings as they come (a
/// fitted tangency curve, a wavy trim), with the pcurve the caller already
/// knows for each, same-parameter over the edge's own range. The seam stays
/// an exact iso-column of the surface, which is what routes around the
/// closed-form refusal in [`make_face_with_pcurves`].
///
/// The caller's contract: both rings are closed, both chart images run the
/// full period *forward* in `u`, and both start on the same column; the
/// seam runs there, between the two start vertices. One ring may be a pole:
/// a degenerate edge whose image is the pole's row over `(0, period)`. The
/// face walks the lower ring forward and the upper one back, whichever
/// order they come in, so it keeps itself on the left of its walk.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if a
/// ring is not closed, the images disagree about the start column or do not
/// run the period forward, or the surface's iso-curve has no closed form.
pub fn make_band_between(
    model: &mut Model,
    surface: &SurfaceGeometry,
    rings: [(&Shape, ogeom_geom::PlanarCurve); 2],
    tol: Tolerances,
) -> OgeomResult<Shape> {
    use ogeom_geom::Curve2d as _;
    use ogeom_geom::Surface as _;

    let surface_id = model.geometry_mut().add_surface(surface.clone());
    let ((ua_dom, ub_dom), _) = surface.domain();
    let span = ub_dom - ua_dom;

    struct Ring {
        edge: Shape,
        vertex: Shape,
        crange: (f64, f64),
        pcurve: ogeom_geom::PlanarCurve,
        start: ogeom_math::Point2,
    }
    let mut prepared: Vec<Ring> = Vec::with_capacity(2);
    for (used, pcurve) in rings {
        let edge = if used.orientation() == ogeom_topo::Orientation::Reversed {
            used.reversed()
        } else {
            used.clone()
        };
        let Some((vertex, other)) = edge_vertices(model, &edge)? else {
            ogeom_bail!(Construction, "a band ring has no vertex");
        };
        if !vertex.is_same(&other) {
            ogeom_bail!(Construction, "a band ring is not closed");
        }
        let crange = {
            let Some(data) = model.node(&edge).and_then(|n| n.data().as_edge()) else {
                ogeom_bail!(Construction, "a band ring holds no edge data");
            };
            // A pole: no curve, and its chart image the pole's row run over
            // one period.
            if data.degenerate {
                (0.0, span)
            } else {
                let Some(EdgeRepr::Curve3d { range, .. }) = data.curve3d() else {
                    ogeom_bail!(Construction, "a band ring has no curve");
                };
                *range
            }
        };
        let start = pcurve.point_at(crange.0, tol)?;
        let end = pcurve.point_at(crange.1, tol)?;
        if (end.x - start.x - span).abs() > 1e-3 {
            ogeom_bail!(
                Construction,
                "a band ring's chart image must run the period forward in u"
            );
        }
        prepared.push(Ring {
            edge,
            vertex,
            crange,
            pcurve,
            start,
        });
    }
    // The band keeps itself on the left of its walk in the chart, which
    // runs the lower ring forward and the upper one back: the rings do not
    // cross, so their start rows say which is lower.
    if prepared[1].start.y < prepared[0].start.y {
        prepared.swap(0, 1);
    }
    // The connector between the rings' starts. On one column it is the
    // surface's own iso; on different columns it is the straight chart
    // segment: a ruling, a parallel arc, or a helix, all exact on a
    // cylinder, and refused by name elsewhere.
    // The rings close round the period, so a start a whole period from the
    // other's is the same column: the second ring's image is moved into the
    // first's period, or the connector would wind a full turn round to it.
    let turns = ((prepared[1].start.x - prepared[0].start.x) / span).round();
    if turns != 0.0 {
        let shift =
            ogeom_math::Transform2::translation(ogeom_math::Vector2::new(-turns * span, 0.0));
        prepared[1].pcurve = prepared[1].pcurve.transformed(&shift, tol)?;
        prepared[1].start = shift.apply(prepared[1].start);
    }
    let (start0, start1) = (prepared[0].start, prepared[1].start);
    let dcol = start1.x - start0.x;
    // The iso is for rings sharing a column to within the confusion
    // distance: the offset times the chart's stretch there, since an
    // off-column iso misses a start vertex by exactly that. Starts further
    // apart take the chart segment, which meets both exactly; a segment
    // spanning less than confusion would be an edge with no length.
    let stretch = {
        let (du, _) = surface.d1_at(start0.x, start0.y, tol)?;
        du.magnitude().max(1.0)
    };
    let (seam, up_is_forward, chart_from, chart_to) = if dcol.abs() * stretch <= tol.confusion() {
        let column = start0.x;
        let (va, vb) = (start0.y, start1.y);
        let Some(seam_curve) = seam_iso_u_curve(surface, column, tol) else {
            ogeom_bail!(
                Construction,
                "the surface's iso-curve has no closed form; no seam can be built"
            );
        };
        let (range, from, to, downward) = if va <= vb {
            (
                (va, vb),
                prepared[0].vertex.clone(),
                prepared[1].vertex.clone(),
                false,
            )
        } else {
            (
                (vb, va),
                prepared[1].vertex.clone(),
                prepared[0].vertex.clone(),
                true,
            )
        };
        let curve_range = (
            iso_curve_parameter_at(surface, range.0),
            iso_curve_parameter_at(surface, range.1),
        );
        // A ring fitted through its stations starts a fit error off the
        // exact column; its start vertex owns that slop.
        for (vertex, at) in [(&from, curve_range.0), (&to, curve_range.1)] {
            let Some(data) = model.node(vertex).and_then(|n| n.data().as_vertex()) else {
                ogeom_bail!(Construction, "a band ring's start vertex holds no point");
            };
            let placed = vertex.transform(model.datums())?.apply(data.point);
            let miss = seam_curve.point_at(at, tol)?.distance(placed);
            if miss > tol.confusion() {
                model.widen(vertex, ogeom_core::Tolerance::new(miss * 2.0)?)?;
            }
        }
        let seam = make_edge_between(model, seam_curve, curve_range, &from, &to, tol)?.shape;
        let (a, b) = (
            ogeom_math::Point2::new(column, range.0),
            ogeom_math::Point2::new(column, range.1),
        );
        (seam, !downward, a, b)
    } else if !matches!(surface, SurfaceGeometry::Cylinder(_)) {
        // On any other surface the chart segment lifts to a curve with no
        // closed form (a loxodrome on a sphere, a skew run on a torus)
        // and is fitted through it, at the chart's own arc length so the
        // curve and its images share one parameter. The fit's error is the
        // vertices' to carry.
        let (a, b, from, to, forward) = if dcol > 0.0 {
            (
                start0,
                start1,
                prepared[0].vertex.clone(),
                prepared[1].vertex.clone(),
                true,
            )
        } else {
            (
                start1,
                start0,
                prepared[1].vertex.clone(),
                prepared[0].vertex.clone(),
                false,
            )
        };
        let length = (b.x - a.x).hypot(b.y - a.y);
        const SAMPLES: usize = 64;
        // Held to the target between the samples as well as at them; the
        // error is the worst measured, and the edge states it below.
        let target = tol.confusion() * 10.0;
        let stations = crate::traced::stations_between((0.0, length), &[], 3, SAMPLES);
        let fitted = crate::traced::fit_traced(
            |t| {
                let f = t / length;
                surface.point_at(a.x + (b.x - a.x) * f, a.y + (b.y - a.y) * f, tol)
            },
            &stations,
            3,
            target,
            tol,
        )?;
        if fitted.error > tol.confusion() * 1e3 {
            ogeom_bail!(
                NotDone,
                "a band's connector reached {} against a target of {target}",
                fitted.error
            );
        }
        let fit_error = fitted.error;
        let connector: ogeom_geom::Curve = ogeom_geom::Curve::BSpline(fitted.curve);
        // Against the vertices' own points, not the lifted chart points: a
        // ring's start vertex stands off the host by the ring's own slop
        // (a fitted seam a few microns off the surface it trims), and the
        // connector, which lies on the host, misses it by that much.
        for (vertex, at) in [(&from, 0.0), (&to, length)] {
            let Some(data) = model.node(vertex).and_then(|n| n.data().as_vertex()) else {
                ogeom_bail!(Construction, "a band ring's start vertex holds no point");
            };
            let placed = vertex.transform(model.datums())?.apply(data.point);
            let miss = connector.point_at(at, tol)?.distance(placed);
            if miss > tol.confusion() {
                model.widen(vertex, ogeom_core::Tolerance::new(miss * 2.0)?)?;
            }
        }
        let seam = make_edge_between(model, connector, (0.0, length), &from, &to, tol)?.shape;
        if fit_error > tol.confusion() {
            model.widen(
                &seam,
                ogeom_core::Tolerance::new(fit_error + tol.confusion())?,
            )?;
        }
        (seam, forward, a, b)
    } else {
        let SurfaceGeometry::Cylinder(c) = surface else {
            ogeom_bail!(
                Construction,
                "a band whose rings start on different columns needs a chart \
                 connector; the surface above was not a cylinder"
            );
        };
        let cylinder = c.cylinder();
        let frame = cylinder.frame();
        // The chart segment, parameterized by the angle so the pcurve's
        // linear map onto it is exact.
        let (a, b, from, to, forward) = if dcol > 0.0 {
            (
                start0,
                start1,
                prepared[0].vertex.clone(),
                prepared[1].vertex.clone(),
                true,
            )
        } else {
            (
                start1,
                start0,
                prepared[1].vertex.clone(),
                prepared[0].vertex.clone(),
                false,
            )
        };
        let drow = b.y - a.y;
        let connector: ogeom_geom::Curve = if drow.abs() <= 1e-9 {
            // A parallel: an arc at one height.
            let lifted = ogeom_math::Frame::new(
                frame.origin() + frame.z().vector() * a.y,
                frame.z(),
                frame.x(),
                tol,
            )?;
            ogeom_geom::CircleCurve::new(ogeom_math::Circle::new(lifted, cylinder.radius(), tol)?)
                .into()
        } else {
            // A helix through both chart points: the pitch is the slope, and
            // the frame rides down so height zero lands where it must.
            let pitch = core::f64::consts::TAU * drow / (b.x - a.x);
            let dropped = ogeom_math::Frame::new(
                frame.origin()
                    + frame.z().vector() * (pitch / core::f64::consts::TAU).mul_add(-a.x, a.y),
                frame.z(),
                frame.x(),
                tol,
            )?;
            ogeom_geom::HelixCurve::over(dropped, cylinder.radius(), pitch, a.x, b.x)?.into()
        };
        let seam = make_edge_between(model, connector, (a.x, b.x), &from, &to, tol)?.shape;
        (seam, forward, a, b)
    };

    // Both rings run forward, so the walk closes as [lo, up, hi rev, down]:
    // the bottom ends at the far column, and the seam's two chart images sit
    // one period apart, assigned to whichever occurrence the walk meets
    // first.
    let segment_line = |shift: f64| -> OgeomResult<ogeom_geom::PlanarCurve> {
        let a = ogeom_math::Point2::new(chart_from.x + shift, chart_from.y);
        let b = ogeom_math::Point2::new(chart_to.x + shift, chart_to.y);
        let d = ogeom_math::Vector2::new(b.x - a.x, b.y - a.y);
        let m = d.x.hypot(d.y);
        if m <= 1e-12 {
            ogeom_bail!(Construction, "the seam's chart image has no length");
        }
        Ok(ogeom_geom::Line2d::over(
            ogeom_math::Axis2::new(a, ogeom_math::Direction2::new(d, tol)?),
            0.0,
            m,
        )?
        .into())
    };
    let (forward_shift, reversed_shift) = if up_is_forward {
        (span, 0.0)
    } else {
        (0.0, span)
    };
    // The pcurves run over their own arc length, start-for-start with the
    // seam's curve range, so that is the window the attachment states.
    let seam_length = (chart_to.x - chart_from.x).hypot(chart_to.y - chart_from.y);
    attach_seam(
        model,
        &seam,
        segment_line(forward_shift)?,
        segment_line(reversed_shift)?,
        surface_id,
        Location::identity(),
        (0.0, seam_length),
    )?;
    crate::state_pcurve_gaps_of(model, core::slice::from_ref(&seam), tol)?;
    for ring in &prepared {
        attach_pcurve(
            model,
            &ring.edge,
            ring.pcurve.clone(),
            surface_id,
            Location::identity(),
            ring.crange,
        )?;
    }

    let up = if up_is_forward {
        seam.clone()
    } else {
        seam.reversed()
    };
    let walk = vec![
        prepared[0].edge.clone(),
        up.clone(),
        prepared[1].edge.reversed(),
        up.reversed(),
    ];
    let wire = make_wire(model, &walk, tol)?.shape;
    Ok(make_face_on(model, surface_id, &[wire], tol)?.shape)
}

/// Build the face of a revolution cap: one closed ring belting a cone, with
/// the apex the file never wrote synthesised as the degenerate ring the band
/// needs.
///
/// A cone face bounded by a single circle has exactly one other boundary the
/// geometry permits, the apex, because the region away from the apex is
/// unbounded. The apex becomes a vertex and a degenerate edge, and the rest
/// is [`make_revolution_band`], one authority for the seam either way.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the surface
/// is not a cone, or the ring resists the band construction.
pub fn make_apex_band(
    model: &mut Model,
    surface: &SurfaceGeometry,
    ring: &Shape,
    tol: Tolerances,
) -> OgeomResult<Shape> {
    let SurfaceGeometry::Cone(c) = surface else {
        ogeom_bail!(
            Construction,
            "only a cone has the single apex a one-ring cap implies"
        );
    };
    let apex = model.add_vertex(VertexData::new(c.cone().apex()));
    let mut data = EdgeData::new();
    data.degenerate = true;
    let edge = model.add_edge(data, &[apex.clone(), apex])?;
    make_revolution_band(model, surface, ring, &edge, tol)
}

/// The chart coordinates of a point on a periodic analytic surface, by
/// closed-form inversion, folded into the surface's own stated window.
///
/// The band synthesis needs the *row* each ring stands on and the *column*
/// the seam anchors to; iterative projection over an imported surface's
/// enormous stated extents can converge to a clamped boundary and place a
/// seam whole units from the vertex it must meet, so the analytic kinds
/// invert exactly instead.
fn analytic_chart_of(
    surface: &SurfaceGeometry,
    p: ogeom_math::Point,
) -> Option<ogeom_math::Point2> {
    use ogeom_geom::Surface as _;
    let tau = core::f64::consts::TAU;
    let raw = match surface {
        SurfaceGeometry::Cylinder(s) => {
            let l = s.cylinder().frame().to_local(p);
            ogeom_math::Point2::new(l.y.atan2(l.x), l.z)
        }
        SurfaceGeometry::Cone(s) => {
            let l = s.cone().frame().to_local(p);
            ogeom_math::Point2::new(l.y.atan2(l.x), l.z)
        }
        SurfaceGeometry::Sphere(s) => {
            let sphere = s.sphere();
            let l = sphere.frame().to_local(p);
            let lat = (l.z / sphere.radius()).clamp(-1.0, 1.0).asin();
            ogeom_math::Point2::new(l.y.atan2(l.x), lat)
        }
        SurfaceGeometry::Torus(s) => {
            let torus = s.torus();
            let l = torus.frame().to_local(p);
            let radial = l.x.hypot(l.y) - torus.major_radius();
            ogeom_math::Point2::new(l.y.atan2(l.x), l.z.atan2(radial))
        }
        _ => return None,
    };
    let ((ua, _), (va, vb)) = surface.domain();
    let u = ua + (raw.x - ua).rem_euclid(tau);
    let v = if surface.is_periodic_v() {
        va + (raw.y - va).rem_euclid(vb - va)
    } else {
        raw.y
    };
    Some(ogeom_math::Point2::new(u, v))
}

/// The revolution axis direction of a periodic analytic surface, the way
/// its angle `u` turns about: [`spin`] of its frame.
fn surface_iso_axis(surface: &ogeom_geom::SurfaceGeometry) -> Option<ogeom_math::Vector> {
    match surface {
        ogeom_geom::SurfaceGeometry::Cylinder(c) => Some(spin(&c.cylinder().frame())),
        ogeom_geom::SurfaceGeometry::Cone(c) => Some(spin(&c.cone().frame())),
        ogeom_geom::SurfaceGeometry::Torus(t) => Some(spin(&t.torus().frame())),
        ogeom_geom::SurfaceGeometry::Sphere(s) => Some(spin(&s.sphere().frame())),
        _ => None,
    }
}

/// The axis a frame's angle turns about: `x × y`, which is `z` on a
/// right-handed frame and `-z` on a left-handed one.
fn spin(frame: &ogeom_math::Frame) -> ogeom_math::Vector {
    frame.z().vector() * frame.handedness().sign()
}

/// The surface's `u = at` iso-curve, parameterized by `v` exactly.
///
/// A ruling on a cylinder, a tube circle on a torus: the curve a seam runs
/// along. `None` where no closed form exists. Shared by the STEP reader's
/// seam synthesis and the healer's ring re-anchoring, so the two cannot
/// disagree about what a seam is.
pub fn surface_iso_u_curve(
    surface: &SurfaceGeometry,
    at: f64,
    tol: Tolerances,
) -> Option<ogeom_geom::Curve> {
    match surface {
        ogeom_geom::SurfaceGeometry::Cylinder(c) => {
            let cylinder = c.cylinder();
            let frame = cylinder.frame();
            let radial = frame.x().vector() * at.cos() + frame.y().vector() * at.sin();
            let location = frame.origin() + radial * cylinder.radius();
            let axis = ogeom_math::Axis {
                location,
                direction: frame.z(),
            };
            Some(ogeom_geom::LineCurve::new(axis).into())
        }
        ogeom_geom::SurfaceGeometry::Torus(t) => {
            let torus = t.torus();
            let frame = torus.frame();
            let radial = frame.x().vector() * at.cos() + frame.y().vector() * at.sin();
            let centre = frame.origin() + radial * torus.major_radius();
            // The tube circle framed so its own angle *is* the surface's v:
            // x toward the outer equator, y along the axis, which makes
            // z = x cross y the tangential direction.
            let circle_frame = ogeom_math::Frame::new(
                centre,
                ogeom_math::Direction::new(radial.cross(frame.z().vector()), tol).ok()?,
                ogeom_math::Direction::new(radial, tol).ok()?,
                tol,
            )
            .ok()?;
            let circle = ogeom_math::Circle::new(circle_frame, torus.minor_radius(), tol).ok()?;
            Some(ogeom_geom::CircleCurve::new(circle).into())
        }
        ogeom_geom::SurfaceGeometry::Cone(c) => {
            let cone = c.cone();
            let frame = cone.frame();
            let radial = frame.x().vector() * at.cos() + frame.y().vector() * at.sin();
            // The ruling through v = 0, arc-length parameterized: the chart's
            // v maps onto it linearly, by t = v / cos(half angle), which the
            // seam construction recovers through `iso_curve_parameter_at`.
            let location = frame.origin() + radial * cone.radius_at(0.0);
            let direction = ogeom_math::Direction::new(
                frame.z().vector() + radial * cone.half_angle().tan(),
                tol,
            )
            .ok()?;
            Some(
                ogeom_geom::LineCurve::new(ogeom_math::Axis {
                    location,
                    direction,
                })
                .into(),
            )
        }
        ogeom_geom::SurfaceGeometry::BSpline(b) => {
            Some(ogeom_geom::Curve::BSpline(b.iso_u_curve(at, tol).ok()?))
        }
        ogeom_geom::SurfaceGeometry::Sphere(sp) => {
            let sphere = sp.sphere();
            let frame = sphere.frame();
            let radial = frame.x().vector() * at.cos() + frame.y().vector() * at.sin();
            // The meridian framed so its own angle is the latitude exactly:
            // x one radius out along the parallel, y toward the north pole.
            let circle_frame = ogeom_math::Frame::new(
                frame.origin(),
                ogeom_math::Direction::new(radial.cross(frame.z().vector()), tol).ok()?,
                ogeom_math::Direction::new(radial, tol).ok()?,
                tol,
            )
            .ok()?;
            let circle = ogeom_math::Circle::new(circle_frame, sphere.radius(), tol).ok()?;
            Some(ogeom_geom::CircleCurve::new(circle).into())
        }
        _ => None,
    }
}

/// The iso-curve a band's seam runs along at column `at`: the closed forms
/// [`surface_iso_u_curve`] names, and an extrusion's generator, the line
/// through its profile's point at `at` along the sweep, on which the
/// chart's `v` is the distance swept.
fn seam_iso_u_curve(
    surface: &SurfaceGeometry,
    at: f64,
    tol: Tolerances,
) -> Option<ogeom_geom::Curve> {
    surface_iso_u_curve(surface, at, tol).or_else(|| match surface {
        ogeom_geom::SurfaceGeometry::Extrusion(swept) => {
            use ogeom_geom::Curve3d as _;
            let location = swept.curve().point_at(at, tol).ok()?;
            Some(
                ogeom_geom::LineCurve::new(ogeom_math::Axis {
                    location,
                    direction: swept.direction(),
                })
                .into(),
            )
        }
        _ => None,
    })
}

/// The parameter on a surface's iso-curve that lands at chart row `v`.
///
/// Identity for the kinds whose iso-curve is parameterized by `v` itself
/// (a cylinder ruling, a torus tube circle, a sphere meridian, an
/// extrusion's generator) and the slant
/// rescale for a cone, whose ruling is arc-length parameterized while the
/// chart's `v` is the height.
fn iso_curve_parameter_at(surface: &SurfaceGeometry, v: f64) -> f64 {
    match surface {
        ogeom_geom::SurfaceGeometry::Cone(c) => v / c.cone().half_angle().cos(),
        _ => v,
    }
}

#[cfg(test)]
#[allow(clippy::unwrap_used, clippy::expect_used)]
mod tests {
    use super::*;

    /// A ring's turn is twice its area wherever it stands: a square 0.125 on
    /// a side turns by 2 × 0.125² about its normal at the origin and 1e8 away
    /// alike. The side is small against the offset, so a sum whose rounding
    /// grows with the distance from the origin would swamp the area.
    #[test]
    fn a_small_ring_far_from_the_origin_turns_by_its_area() {
        let tol = ogeom_core::Tolerances::millimetres();
        let side = 0.125;
        for offset in [0.0, 1.0e6, 1.0e7, 1.0e8] {
            let mut model = Model::new();
            let corners = [(0.0, 0.0), (side, 0.0), (side, side), (0.0, side)]
                .map(|(x, y)| Point::new(offset + x, -offset + y, 0.0));
            let ring = make_polygon(&mut model, &corners, true, tol).unwrap().shape;
            let z = ogeom_math::Vector::new(0.0, 0.0, 1.0);
            // Corners 1e8 out are exact (0.125 is a power of two), but the
            // points sampled between them are rounded to 1.5e-8 there, a
            // relative 1e-7 of the side.
            approx::assert_relative_eq!(
                wire_turn(&model, &ring, z, tol).unwrap(),
                2.0 * side * side,
                max_relative = 1e-6
            );
        }
    }

    /// A plane keeps its normal and the face walks its rings about it:
    /// an outline wound clockwise and a hole wound counter-clockwise are
    /// both walked back.
    #[test]
    fn a_planar_face_walks_its_rings_about_its_plane() {
        let tol = ogeom_core::Tolerances::millimetres();
        let mut model = Model::new();
        let ring = |model: &mut Model, corners: &[(f64, f64)]| {
            let points: Vec<Point> = corners
                .iter()
                .map(|&(x, y)| Point::new(x, y, 0.0))
                .collect();
            make_polygon(model, &points, true, tol).unwrap().shape
        };
        let outline = ring(
            &mut model,
            &[(0.0, 0.0), (0.0, 4.0), (4.0, 4.0), (4.0, 0.0)],
        );
        let hole = ring(
            &mut model,
            &[(1.0, 1.0), (3.0, 1.0), (3.0, 3.0), (1.0, 3.0)],
        );
        let plane = ogeom_math::Plane::new(ogeom_math::Frame::WORLD);
        let face = make_face(
            &mut model,
            ogeom_geom::PlaneSurface::new(plane).into(),
            &[outline.clone(), hole.clone()],
            tol,
        )
        .unwrap();
        let z = ogeom_math::Vector::new(0.0, 0.0, 1.0);
        let wires = model.children_of(&face.shape).unwrap();
        assert!(wire_turn(&model, &wires[0], z, tol).unwrap() > 0.0);
        assert!(wire_turn(&model, &wires[1], z, tol).unwrap() < 0.0);
        assert_eq!(
            face.history.generated(&outline),
            std::slice::from_ref(&face.shape)
        );
        assert_eq!(
            face.history.generated(&hole),
            std::slice::from_ref(&face.shape)
        );
        use ogeom_geom::Surface as _;
        let data = model.node(&face.shape).unwrap().data().as_face().unwrap();
        let normal = model
            .geometry()
            .surface(data.surface)
            .unwrap()
            .normal_at(0.0, 0.0, tol)
            .unwrap();
        assert!(normal.vector().dot(z) > 0.0);
    }

    /// A face on a fitted surface builds with fitted trims: the closed-form
    /// refusal falls back to the projected fit, which is what lets a
    /// defeaturing rebuild stand a face on a spline neighbour.
    #[test]
    fn a_face_on_a_fitted_surface_gains_projected_pcurves() {
        let tol = ogeom_core::Tolerances::millimetres();
        let mut model = Model::new();
        // A gently wavy fitted patch, nothing analytic recognises.
        let rows: Vec<Vec<ogeom_math::Point>> = (0..5)
            .map(|j| {
                (0..5)
                    .map(|i| {
                        let (x, y) = (f64::from(i) * 2.5, f64::from(j) * 2.5);
                        ogeom_math::Point::new(x, y, 0.4 * (x * 0.7).sin() * (y * 0.5).cos())
                    })
                    .collect()
            })
            .collect();
        let surface: ogeom_geom::SurfaceGeometry =
            ogeom_geom::fit::fit_surface_grid(&rows, 3, 1e-6, tol)
                .unwrap()
                .curve
                .into();
        // The border iso-curves, straight off the fitted chart.
        use ogeom_geom::Surface as _;
        let ((ua, ub), (va, vb)) = surface.domain();
        let iso = |fixed_u: Option<f64>, fixed_v: Option<f64>| -> Curve {
            const N: usize = 33;
            let pts: Vec<ogeom_math::Point> = (0..N)
                .map(|k| {
                    let t = f64::from(u32::try_from(k).unwrap())
                        / f64::from(u32::try_from(N - 1).unwrap());
                    let (u, v) = match (fixed_u, fixed_v) {
                        (Some(u), None) => (u, va + (vb - va) * t),
                        (None, Some(v)) => (ua + (ub - ua) * t, v),
                        _ => unreachable!(),
                    };
                    surface.point_at(u, v, tol).unwrap()
                })
                .collect();
            Curve::BSpline(
                ogeom_geom::fit::fit_points(&pts, 3, 1e-9, tol)
                    .unwrap()
                    .curve,
            )
        };
        // Shared corner vertices, widened to absorb the border fits' own
        // slack, so the ring connects by node rather than by luck.
        let corners: Vec<Shape> = [(ua, va), (ub, va), (ub, vb), (ua, vb)]
            .into_iter()
            .map(|(u, v)| {
                let p = surface.point_at(u, v, tol).unwrap();
                let vertex = model.add_vertex(ogeom_topo::VertexData::new(p));
                if let Some(node) = model.node_mut(&vertex)
                    && let ogeom_topo::NodeData::Vertex(data) = node.data_mut()
                {
                    data.tolerance = data.tolerance.widen_to(1e-4);
                }
                vertex
            })
            .collect();
        let ring = [
            (iso(None, Some(va)), 0usize, 1usize),
            (iso(Some(ub), None), 1, 2),
            (iso(None, Some(vb)), 3, 2),
            (iso(Some(ua), None), 0, 3),
        ];
        let edges: Vec<Shape> = ring
            .into_iter()
            .enumerate()
            .map(|(k, (c, from, to))| {
                use ogeom_geom::Curve3d as _;
                let domain = c.domain();
                let edge =
                    make_edge_between(&mut model, c, domain, &corners[from], &corners[to], tol)
                        .unwrap()
                        .shape;
                if k >= 2 { edge.reversed() } else { edge }
            })
            .collect();
        let built =
            make_face_with_pcurves(&mut model, surface, std::slice::from_ref(&edges), tol).unwrap();
        // Every edge carries a pcurve for the face's surface.
        let data = model
            .node(&built.shape)
            .unwrap()
            .data()
            .as_face()
            .unwrap()
            .clone();
        for edge in &edges {
            let e = model.node(edge).unwrap().data().as_edge().unwrap();
            assert!(
                e.pcurve_for(data.surface, edge.location()).is_some(),
                "a fitted trim was attached"
            );
        }
    }

    #[test]
    fn a_band_whose_rings_start_apart_gets_a_helical_connector() {
        // Two circles at different heights whose parameter origins differ:
        // the seam has no single column to run down, so it runs the straight
        // chart segment: a helix on the cylinder, exactly.
        let mut model = Model::new();
        let radius = 2.0;
        let cylinder = ogeom_geom::CylinderSurface::new(
            ogeom_math::Cylinder::new(Frame::WORLD, radius, T).unwrap(),
            (-1.0, 4.0),
        )
        .unwrap();
        let surface: SurfaceGeometry = cylinder.into();

        let ring_at = |model: &mut Model, height: f64, phase: f64| {
            let x = ogeom_math::Direction::new(
                ogeom_math::Vector::new(phase.cos(), phase.sin(), 0.0),
                T,
            )
            .unwrap();
            let frame =
                Frame::new(Point::new(0.0, 0.0, height), ogeom_math::Direction::Z, x, T).unwrap();
            let curve: Curve = CircleCurve::new(Circle::new(frame, radius, T).unwrap()).into();
            let domain = curve.domain();
            let edge = make_edge(model, curve, domain, T).unwrap().shape;
            let pcurve: ogeom_geom::PlanarCurve = ogeom_geom::Line2d::over(
                ogeom_math::Axis2::new(
                    ogeom_math::Point2::new(phase, height),
                    ogeom_math::Direction2::new(ogeom_math::Vector2::new(1.0, 0.0), T).unwrap(),
                ),
                0.0,
                core::f64::consts::TAU,
            )
            .unwrap()
            .into();
            (edge, pcurve)
        };
        let (lo, lo_pcurve) = ring_at(&mut model, 0.0, 0.0);
        let (hi, hi_pcurve) = ring_at(&mut model, 2.0, 0.4);

        let band = make_band_between(
            &mut model,
            &surface,
            [(&lo, lo_pcurve), (&hi, hi_pcurve)],
            T,
        )
        .unwrap();

        assert_eq!(
            explore_unique(&model, &band, ShapeType::Edge)
                .unwrap()
                .len(),
            3
        );
        // The connector's curve is a genuine helix through both starts.
        let helix = explore_unique(&model, &band, ShapeType::Edge)
            .unwrap()
            .into_iter()
            .find_map(|e| {
                let data = model.node(&e)?.data().as_edge()?;
                let ogeom_topo::EdgeRepr::Curve3d { curve, .. } = data.curve3d()? else {
                    return None;
                };
                match model.geometry().curve(*curve)? {
                    Curve::Helix(h) => Some(*h),
                    _ => None,
                }
            })
            .expect("the connector is a helix");
        use ogeom_geom::Curve3d as _;
        let at_start = Curve::Helix(helix).point_at(0.0, T).unwrap();
        let at_end = Curve::Helix(helix).point_at(0.4, T).unwrap();
        assert!(at_start.distance(Point::new(2.0, 0.0, 0.0)) < 1e-9);
        assert!(at_end.distance(Point::new(2.0 * 0.4_f64.cos(), 2.0 * 0.4_f64.sin(), 2.0)) < 1e-9);
    }

    /// A boundary with more detail than the projected fit's samples, on a
    /// patch no closed form speaks: each fitted trim, lifted, stands within
    /// what its edge states.
    #[test]
    fn a_fitted_trim_states_its_gap_along_the_whole_edge() {
        let mut model = Model::new();
        let patch = ogeom_geom::BSplineSurface::new(
            ogeom_math::KnotVector::clamped_uniform(1, 2).unwrap(),
            ogeom_math::KnotVector::clamped_uniform(1, 2).unwrap(),
            &ogeom_math::ControlGrid::new(
                vec![
                    Point::new(-1.0, -1.0, 0.0),
                    Point::new(-1.0, 11.0, 0.0),
                    Point::new(11.0, -1.0, 0.0),
                    Point::new(11.0, 11.0, 0.0),
                ],
                2,
                2,
            )
            .unwrap(),
            T,
        )
        .unwrap();
        const SPANS: usize = 300;
        let control: Vec<Point> = (0..SPANS + 3)
            .map(|i| {
                #[allow(clippy::cast_precision_loss)]
                let x = 10.0 * i as f64 / (SPANS + 2) as f64;
                let y = if i == 0 || i == SPANS + 2 {
                    0.0
                } else if i.is_multiple_of(2) {
                    0.03
                } else {
                    -0.03
                };
                Point::new(x, y, 0.0)
            })
            .collect();
        let wiggle: Curve = ogeom_geom::BSplineCurve::new(
            ogeom_math::KnotVector::clamped_uniform(3, control.len()).unwrap(),
            control,
            T,
        )
        .unwrap()
        .into();
        let corners = [
            Point::new(0.0, 0.0, 0.0),
            Point::new(10.0, 0.0, 0.0),
            Point::new(10.0, 10.0, 0.0),
            Point::new(0.0, 10.0, 0.0),
        ];
        let v: Vec<Shape> = corners
            .iter()
            .map(|p| make_vertex(&mut model, *p).shape)
            .collect();
        let mut edges = vec![
            make_edge_between(&mut model, wiggle, (0.0, 1.0), &v[0], &v[1], T)
                .unwrap()
                .shape,
        ];
        for i in 1..4 {
            let (a, b) = (corners[i], corners[(i + 1) % 4]);
            let line: Curve = LineCurve::segment(a, b, T).unwrap().into();
            edges.push(
                make_edge_between(
                    &mut model,
                    line,
                    (0.0, a.distance(b)),
                    &v[i],
                    &v[(i + 1) % 4],
                    T,
                )
                .unwrap()
                .shape,
            );
        }
        let face = make_face_with_pcurves(&mut model, SurfaceGeometry::BSpline(patch), &[edges], T)
            .unwrap()
            .shape;
        for edge in explore_unique(&model, &face, ShapeType::Edge).unwrap() {
            let stated = model
                .node(&edge)
                .unwrap()
                .data()
                .as_edge()
                .unwrap()
                .tolerance
                .get();
            let gap = crate::edge_pcurve_gap(&model, &edge, T).unwrap().unwrap();
            assert!(gap <= stated, "a trim {gap} off an edge stating {stated}");
        }
        let diagnosis = crate::check(&model, &face, T).unwrap();
        assert!(diagnosis.is_valid(), "{:?}", diagnosis.problems);
    }

    /// On a sphere the chart segment between two starts lifts to a curve
    /// with no closed form, and the connector is fitted through it: its
    /// charts stand within what the edge states, everywhere along it.
    #[test]
    fn a_fitted_connector_states_its_charts_gap() {
        let mut model = Model::new();
        let radius = 5.0;
        let surface: SurfaceGeometry = ogeom_geom::SphereSurface::new(
            ogeom_math::Sphere::new(Frame::WORLD, radius, T).unwrap(),
        )
        .into();
        let ring_at = |model: &mut Model, lat: f64, phase: f64| {
            let x = ogeom_math::Direction::new(
                ogeom_math::Vector::new(phase.cos(), phase.sin(), 0.0),
                T,
            )
            .unwrap();
            let frame = Frame::new(
                Point::new(0.0, 0.0, radius * lat.sin()),
                ogeom_math::Direction::Z,
                x,
                T,
            )
            .unwrap();
            let curve: Curve =
                CircleCurve::new(Circle::new(frame, radius * lat.cos(), T).unwrap()).into();
            let domain = curve.domain();
            let edge = make_edge(model, curve, domain, T).unwrap().shape;
            let pcurve: ogeom_geom::PlanarCurve = ogeom_geom::Line2d::over(
                ogeom_math::Axis2::new(
                    ogeom_math::Point2::new(phase, lat),
                    ogeom_math::Direction2::new(ogeom_math::Vector2::new(1.0, 0.0), T).unwrap(),
                ),
                0.0,
                core::f64::consts::TAU,
            )
            .unwrap()
            .into();
            (edge, pcurve)
        };
        let (lo, lo_pcurve) = ring_at(&mut model, -1.2, 0.0);
        let (hi, hi_pcurve) = ring_at(&mut model, 1.2, 3.0);
        let band = make_band_between(
            &mut model,
            &surface,
            [(&lo, lo_pcurve), (&hi, hi_pcurve)],
            T,
        )
        .unwrap();
        for edge in explore_unique(&model, &band, ShapeType::Edge).unwrap() {
            let stated = model
                .node(&edge)
                .unwrap()
                .data()
                .as_edge()
                .unwrap()
                .tolerance
                .get();
            let gap = crate::edge_pcurve_gap(&model, &edge, T).unwrap().unwrap();
            assert!(gap <= stated, "a chart {gap} off an edge stating {stated}");
        }
        let diagnosis = crate::check(&model, &band, T).unwrap();
        assert!(diagnosis.is_valid(), "{:?}", diagnosis.problems);
    }

    #[test]
    fn a_band_between_a_circle_and_a_fitted_wavy_ring_closes_with_its_seam() {
        use ogeom_geom::Curve2d as _;
        use ogeom_geom::Curve3d as _;
        use ogeom_geom::Surface as _;

        let mut model = Model::new();
        let radius = 2.0;
        let cylinder = ogeom_geom::CylinderSurface::new(
            ogeom_math::Cylinder::new(Frame::WORLD, radius, T).unwrap(),
            (-1.0, 4.0),
        )
        .unwrap();
        let surface: SurfaceGeometry = cylinder.into();

        // The lower ring: the exact circle at z = 0, its chart image the row.
        let lo_curve: Curve =
            CircleCurve::new(Circle::new(Frame::WORLD, radius, T).unwrap()).into();
        let lo_domain = lo_curve.domain();
        let lo = make_edge(&mut model, lo_curve, lo_domain, T).unwrap().shape;
        let lo_pcurve: ogeom_geom::PlanarCurve = ogeom_geom::Line2d::over(
            ogeom_math::Axis2::new(
                ogeom_math::Point2::new(0.0, 0.0),
                ogeom_math::Direction2::new(ogeom_math::Vector2::new(1.0, 0.0), T).unwrap(),
            ),
            lo_domain.0,
            lo_domain.1,
        )
        .unwrap()
        .into();

        // The upper ring: a wavy loop on the cylinder, fitted jointly with
        // its own chart image so curve and pcurve stay same-parameter.
        let tau = core::f64::consts::TAU;
        let samples = 96;
        let mut points = Vec::with_capacity(samples + 1);
        let mut chart = Vec::with_capacity(samples + 1);
        for i in 0..=samples {
            #[allow(clippy::cast_precision_loss)]
            let u = tau * (i as f64) / (samples as f64);
            let v = (3.0 * u).sin().mul_add(0.2, 1.5);
            points.push(Point::new(radius * u.cos(), radius * u.sin(), v));
            chart.push(ogeom_math::Point2::new(u, v));
        }
        let (fitted, hi_chart, _) =
            ogeom_geom::fit::fit_points_joint_closed(&points, &chart, &chart, 3, 1e-4, T).unwrap();
        assert!(fitted.met, "the wavy ring fit should meet its tolerance");
        let hi_curve: Curve = fitted.curve.into();
        let hi_domain = hi_curve.domain();
        let hi = make_edge(&mut model, hi_curve.clone(), hi_domain, T)
            .unwrap()
            .shape;

        let band = make_band_between(
            &mut model,
            &surface,
            [(&lo, lo_pcurve), (&hi, hi_chart.into())],
            T,
        )
        .unwrap();

        // Four traversals, three distinct edges: two rings and the seam,
        // used once per side.
        let occurrences =
            ogeom_topo::explore(&model, &band, ogeom_topo::Filter::OfType(ShapeType::Edge))
                .unwrap();
        assert_eq!(occurrences.len(), 4);
        assert_eq!(
            explore_unique(&model, &band, ShapeType::Edge)
                .unwrap()
                .len(),
            3
        );

        // The wavy ring's chart image lands on the curve itself: the fit is
        // same-parameter, and the surface evaluates the image back onto it.
        let hi_pcurve = {
            let data = model.node(&hi).unwrap().data().as_edge().unwrap();
            let repr = data
                .representations
                .iter()
                .find_map(|r| match r {
                    ogeom_topo::EdgeRepr::PCurve { curve, .. } => Some(*curve),
                    _ => None,
                })
                .expect("the ring carries its pcurve");
            model.geometry().pcurve(repr).unwrap().clone()
        };
        for i in 0..8 {
            #[allow(clippy::cast_precision_loss)]
            let t = hi_domain.0 + (hi_domain.1 - hi_domain.0) * (i as f64) / 8.0;
            let uv = hi_pcurve.point_at(t, T).unwrap();
            let through_chart = surface.point_at(uv.x, uv.y, T).unwrap();
            let direct = hi_curve.point_at(t, T).unwrap();
            assert!(
                through_chart.distance(direct) < 1e-4,
                "same-parameter drift at t = {t}"
            );
        }
    }

    use ogeom_geom::{CircleCurve, LineCurve, PlaneSurface};
    use ogeom_math::{Circle, Frame, Plane};
    use ogeom_topo::explore_unique;

    const T: Tolerances = Tolerances::millimetres();

    fn segment(a: Point, b: Point) -> Curve {
        LineCurve::segment(a, b, T).unwrap().into()
    }

    /// Four edges forming a closed square in the xy plane.
    fn square_edges(model: &mut Model) -> Vec<Shape> {
        let corners = [
            Point::new(0.0, 0.0, 0.0),
            Point::new(1.0, 0.0, 0.0),
            Point::new(1.0, 1.0, 0.0),
            Point::new(0.0, 1.0, 0.0),
        ];
        let mut vertices: Vec<Shape> = corners
            .iter()
            .map(|p| model.add_vertex(VertexData::new(*p)))
            .collect();
        vertices.push(vertices[0].clone());

        (0..4)
            .map(|i| {
                let curve = segment(corners[i], corners[(i + 1) % 4]);
                let length = corners[i].distance(corners[(i + 1) % 4]);
                make_edge_between(
                    model,
                    curve,
                    (0.0, length),
                    &vertices[i],
                    &vertices[i + 1],
                    T,
                )
                .unwrap()
                .shape
            })
            .collect()
    }

    /// A face on a patch whose chart spans a fortieth of a unit, its
    /// pcurves meeting at one corner 6e-8 apart in the chart and less than
    /// the confusion distance apart in space: the miss is fitting slop, as
    /// a miss of a millionth of the chart's span is. It is too short in the
    /// chart to be an edge of its own, so no pole is laid there, and the
    /// wire comes back as it was.
    #[test]
    fn a_small_chart_miss_whose_ends_meet_in_space_lays_no_pole() {
        use ogeom_geom::{BSplineSurface, Line2d, SurfaceGeometry};
        use ogeom_math::{ControlGrid, KnotVector, Point2};
        let side = 0.025;
        let mut model = Model::new();
        // The patch p(u, v) = (u, v, 0) over [0, side] each way.
        let knots = || KnotVector::new(vec![0.0, 0.0, side, side], 1).unwrap();
        let grid = ControlGrid::new(
            vec![
                Point::new(0.0, 0.0, 0.0),
                Point::new(0.0, side, 0.0),
                Point::new(side, 0.0, 0.0),
                Point::new(side, side, 0.0),
            ],
            2,
            2,
        )
        .unwrap();
        let surface: SurfaceGeometry = BSplineSurface::new(knots(), knots(), &grid, T)
            .unwrap()
            .into();
        let surface_id = model.geometry_mut().add_surface(surface);
        let corners = [(0.0, 0.0), (side, 0.0), (side, side), (0.0, side)];
        let vertices: Vec<Shape> = corners
            .iter()
            .map(|&(x, y)| model.add_vertex(VertexData::new(Point::new(x, y, 0.0))))
            .collect();
        let miss = 6e-8;
        let edges: Vec<Shape> = (0..4)
            .map(|i| {
                let (a, b) = (corners[i], corners[(i + 1) % 4]);
                let (pa, pb) = (Point::new(a.0, a.1, 0.0), Point::new(b.0, b.1, 0.0));
                let edge = make_edge_between(
                    &mut model,
                    segment(pa, pb),
                    (0.0, side),
                    &vertices[i],
                    &vertices[(i + 1) % 4],
                    T,
                )
                .unwrap()
                .shape;
                // The second edge's pcurve starts a hair up the chart from
                // where the first one's ends.
                let start = if i == 1 {
                    Point2::new(a.0, a.1 + miss)
                } else {
                    Point2::new(a.0, a.1)
                };
                let end = Point2::new(b.0, b.1);
                attach_pcurve(
                    &mut model,
                    &edge,
                    Line2d::segment(start, end, T).unwrap().into(),
                    surface_id,
                    Location::identity(),
                    (0.0, start.distance(end)),
                )
                .unwrap();
                edge
            })
            .collect();
        let closed = closed_at_poles(&mut model, surface_id, &edges, T).unwrap();
        assert_eq!(closed.len(), 4, "no pole edge is laid in the miss");
        for (got, given) in closed.iter().zip(&edges) {
            assert!(got.is_same(given));
        }
    }

    #[test]
    fn an_edge_takes_its_vertices_from_its_own_geometry() {
        // Deriving them means the topology and the geometry cannot disagree
        // about where the edge starts and stops.
        let mut model = Model::new();
        let built = make_edge(
            &mut model,
            segment(Point::ORIGIN, Point::new(3.0, 4.0, 0.0)),
            (0.0, 5.0),
            T,
        )
        .unwrap();

        let (start, end) = edge_vertices(&model, &built.shape).unwrap().unwrap();
        let point_of = |v: &Shape| model.node(v).unwrap().data().as_vertex().unwrap().point;
        assert!(point_of(&start).is_equal(Point::ORIGIN, T));
        assert!(point_of(&end).is_equal(Point::new(3.0, 4.0, 0.0), T));
        assert_eq!(built.history.generated(&built.shape).len(), 2);
    }

    #[test]
    fn a_closed_edge_names_one_vertex_twice() {
        // Two coincident vertices would leave the wire looking open at the join,
        // and the gap only surfaces when something tries to walk the boundary.
        let mut model = Model::new();
        let circle: Curve = CircleCurve::new(Circle::new(Frame::WORLD, 2.0, T).unwrap()).into();
        let built = make_edge(&mut model, circle, (0.0, core::f64::consts::TAU), T).unwrap();

        let (start, end) = edge_vertices(&model, &built.shape).unwrap().unwrap();
        assert!(start.is_same(&end), "a closed edge starts where it ends");
        assert_eq!(
            explore_unique(&model, &built.shape, ShapeType::Vertex)
                .unwrap()
                .len(),
            1
        );
    }

    #[test]
    fn an_edge_whose_curve_does_not_reach_its_vertices_is_refused() {
        // The gap would show up much later, in whatever first walks the
        // boundary, with nothing left to say where it came from.
        let mut model = Model::new();
        let a = model.add_vertex(VertexData::new(Point::ORIGIN));
        let b = model.add_vertex(VertexData::new(Point::new(10.0, 0.0, 0.0)));
        let curve = segment(Point::ORIGIN, Point::new(5.0, 0.0, 0.0));

        assert!(
            make_edge_between(&mut model, curve.clone(), (0.0, 5.0), &a, &b, T).is_err(),
            "the curve stops half way"
        );
        let c = model.add_vertex(VertexData::new(Point::new(5.0, 0.0, 0.0)));
        assert!(make_edge_between(&mut model, curve, (0.0, 5.0), &a, &c, T).is_ok());
    }

    #[test]
    fn a_vertex_with_a_loose_tolerance_admits_a_curve_that_stops_short() {
        // The reach is the vertex's own tolerance, not a fixed epsilon: a
        // vertex that has been widened by an earlier repair genuinely does
        // occupy that much space.
        let mut model = Model::new();
        let a = model.add_vertex(VertexData::new(Point::ORIGIN));
        let b = model
            .add_vertex(VertexData::with_tolerance(Point::new(5.001, 0.0, 0.0), 1e-2).unwrap());
        let curve = segment(Point::ORIGIN, Point::new(5.0, 0.0, 0.0));
        assert!(make_edge_between(&mut model, curve, (0.0, 5.0), &a, &b, T).is_ok());
    }

    #[test]
    fn edge_vertices_follow_the_direction_of_travel() {
        let mut model = Model::new();
        let built = make_edge(
            &mut model,
            segment(Point::ORIGIN, Point::new(1.0, 0.0, 0.0)),
            (0.0, 1.0),
            T,
        )
        .unwrap();

        let (forward_start, forward_end) = edge_vertices(&model, &built.shape).unwrap().unwrap();
        let (back_start, back_end) = edge_vertices(&model, &built.shape.reversed())
            .unwrap()
            .unwrap();

        assert!(back_start.is_same(&forward_end));
        assert!(back_end.is_same(&forward_start));
    }

    #[test]
    fn a_wire_whose_edges_do_not_meet_is_refused() {
        // Model::add_wire cannot catch this: it has no geometry. This is the
        // check that keeps a face from being built on a boundary with a gap.
        let mut model = Model::new();
        let joined = make_edge(
            &mut model,
            segment(Point::ORIGIN, Point::new(1.0, 0.0, 0.0)),
            (0.0, 1.0),
            T,
        )
        .unwrap()
        .shape;
        let apart = make_edge(
            &mut model,
            segment(Point::new(5.0, 0.0, 0.0), Point::new(6.0, 0.0, 0.0)),
            (0.0, 1.0),
            T,
        )
        .unwrap()
        .shape;

        let err = make_wire(&mut model, &[joined.clone(), apart], T).unwrap_err();
        assert!(err.to_string().contains("gap"), "unexpected message: {err}");
        assert!(make_wire(&mut model, &[joined], T).is_ok());
    }

    #[test]
    fn a_wire_of_edges_that_meet_is_accepted_and_reported_closed() {
        let mut model = Model::new();
        let edges = square_edges(&mut model);
        let wire = make_wire(&mut model, &edges, T).unwrap();
        assert!(is_wire_closed(&model, &wire.shape, T).unwrap());

        // The edges are generated into the wire, not consumed by it: an edge
        // between two faces belongs to both wires.
        for edge in &edges {
            assert!(!wire.history.is_deleted(edge));
            assert_eq!(wire.history.generated(edge).len(), 1);
        }
    }

    #[test]
    fn an_open_wire_is_reported_open() {
        let mut model = Model::new();
        let edges = square_edges(&mut model);
        let open = make_wire(&mut model, &edges[..3], T).unwrap();
        assert!(!is_wire_closed(&model, &open.shape, T).unwrap());
    }

    #[test]
    fn a_face_refuses_an_open_boundary() {
        let mut model = Model::new();
        let edges = square_edges(&mut model);
        let open = make_wire(&mut model, &edges[..3], T).unwrap().shape;
        let plane: SurfaceGeometry = PlaneSurface::new(Plane::new(Frame::WORLD)).into();

        let err = make_face(&mut model, plane.clone(), &[open], T).unwrap_err();
        assert!(
            err.to_string().contains("open"),
            "unexpected message: {err}"
        );

        let closed = make_wire(&mut model, &edges, T).unwrap().shape;
        assert!(make_face(&mut model, plane, &[closed], T).is_ok());
    }

    #[test]
    fn a_natural_face_needs_no_wires_at_all() {
        let mut model = Model::new();
        let built = make_natural_face(
            &mut model,
            PlaneSurface::new(Plane::new(Frame::WORLD)).into(),
        )
        .unwrap();
        let node = model.node(&built.shape).unwrap();
        assert!(node.data().as_face().unwrap().natural_restriction);
        assert!(built.history.is_empty());
    }

    #[test]
    fn a_shell_of_one_face_is_not_closed() {
        // A single face has free edges all round, so it bounds nothing.
        let mut model = Model::new();
        let edges = square_edges(&mut model);
        let wire = make_wire(&mut model, &edges, T).unwrap().shape;
        let face = make_face(
            &mut model,
            PlaneSurface::new(Plane::new(Frame::WORLD)).into(),
            &[wire],
            T,
        )
        .unwrap()
        .shape;
        let shell = make_shell(&mut model, &[face]).unwrap();
        assert!(!is_shell_closed(&model, &shell.shape).unwrap());
    }

    #[test]
    fn two_faces_sharing_every_edge_form_a_closed_shell() {
        // The degenerate closed shell: the same square from both sides. Every
        // edge is used exactly twice, which is what closure means.
        let mut model = Model::new();
        let edges = square_edges(&mut model);
        let plane: SurfaceGeometry = PlaneSurface::new(Plane::new(Frame::WORLD)).into();

        let front_wire = make_wire(&mut model, &edges, T).unwrap().shape;
        let reversed: Vec<Shape> = edges.iter().rev().map(Shape::reversed).collect();
        let back_wire = make_wire(&mut model, &reversed, T).unwrap().shape;

        let front = make_face(&mut model, plane.clone(), &[front_wire], T)
            .unwrap()
            .shape;
        let back = make_face(&mut model, plane, &[back_wire], T).unwrap().shape;
        let shell = make_shell(&mut model, &[front, back]).unwrap().shape;

        assert!(is_shell_closed(&model, &shell).unwrap());
        assert_eq!(
            explore_unique(&model, &shell, ShapeType::Edge)
                .unwrap()
                .len(),
            4
        );
        assert!(make_solid(&mut model, &[shell]).is_ok());
    }

    #[test]
    fn a_pcurve_can_be_attached_after_the_edge_exists() {
        // An edge learns about a face's parameter space when it joins that
        // face, not when it is created; it may join several.
        let mut model = Model::new();
        let edge = make_edge(
            &mut model,
            segment(Point::ORIGIN, Point::new(1.0, 0.0, 0.0)),
            (0.0, 1.0),
            T,
        )
        .unwrap()
        .shape;

        let surface = model
            .geometry_mut()
            .add_surface(PlaneSurface::new(Plane::new(Frame::WORLD)).into());
        let pcurve = ogeom_geom::Line2d::segment(
            ogeom_math::Point2::ORIGIN,
            ogeom_math::Point2::new(1.0, 0.0),
            T,
        )
        .unwrap()
        .into();
        attach_pcurve(
            &mut model,
            &edge,
            pcurve,
            surface,
            Location::identity(),
            (0.0, 1.0),
        )
        .unwrap();

        let data = model.node(&edge).unwrap().data().as_edge().unwrap();
        assert_eq!(data.representations.len(), 2, "the curve and now a pcurve");
        assert!(data.pcurve_on(surface).is_some());
        assert!(
            !data.same_parameter(),
            "the new representation has not been shown to agree with the curve"
        );
    }

    #[test]
    fn degenerate_edge_ranges_are_refused() {
        let mut model = Model::new();
        let curve = segment(Point::ORIGIN, Point::new(1.0, 0.0, 0.0));
        assert!(make_edge(&mut model, curve.clone(), (1.0, 0.0), T).is_err());
        assert!(make_edge(&mut model, curve.clone(), (0.5, 0.5), T).is_err());
        assert!(make_edge(&mut model, curve.clone(), (0.0, f64::NAN), T).is_err());
        assert!(make_edge(&mut model, curve, (0.0, 1.0), T).is_ok());
    }

    #[test]
    fn builders_reject_children_of_the_wrong_type() {
        let mut model = Model::new();
        let vertex = model.add_vertex(VertexData::new(Point::ORIGIN));
        assert!(make_wire(&mut model, std::slice::from_ref(&vertex), T).is_err());
        assert!(make_shell(&mut model, std::slice::from_ref(&vertex)).is_err());
        assert!(make_solid(&mut model, &[vertex]).is_err());
        assert!(make_wire(&mut model, &[], T).is_err());
    }

    #[test]
    fn an_edges_vertices_record_where_they_came_from() {
        // Provenance is what a rebuild matches against; a vertex that does not
        // say which edge produced it cannot be found again.
        let mut model = Model::new();
        model.begin_operation();
        let built = make_edge(
            &mut model,
            segment(Point::ORIGIN, Point::new(1.0, 0.0, 0.0)),
            (0.0, 1.0),
            T,
        )
        .unwrap();

        let (start, end) = edge_vertices(&model, &built.shape).unwrap().unwrap();
        for (vertex, role) in [(start, roles::EDGE_START), (end, roles::EDGE_END)] {
            let provenance = model.provenance_of(&vertex).unwrap();
            assert!(
                matches!(provenance, ogeom_core::Provenance::Derived { role: r, .. } if *r == role),
                "expected a derived vertex with role {role:?}, got {provenance:?}"
            );
        }
    }
}

#[cfg(test)]
#[allow(clippy::unwrap_used, clippy::expect_used)]
mod polygon_tests {
    use super::*;
    use ogeom_geom::{CircleCurve, PlaneSurface};
    use ogeom_math::{Circle, Frame, Plane};
    use ogeom_topo::explore_unique;

    const T: Tolerances = Tolerances::millimetres();

    fn square() -> Vec<Point> {
        vec![
            Point::new(0.0, 0.0, 0.0),
            Point::new(2.0, 0.0, 0.0),
            Point::new(2.0, 2.0, 0.0),
            Point::new(0.0, 2.0, 0.0),
        ]
    }

    #[test]
    fn a_closed_polygon_names_its_first_vertex_again_rather_than_a_second_one() {
        // Two coincident vertices would leave the wire looking open at the
        // join, and the gap only surfaces when something walks the boundary.
        let mut model = Model::new();
        let built = make_polygon(&mut model, &square(), true, T).unwrap();

        assert_eq!(
            explore_unique(&model, &built.shape, ShapeType::Edge)
                .unwrap()
                .len(),
            4
        );
        assert_eq!(
            explore_unique(&model, &built.shape, ShapeType::Vertex)
                .unwrap()
                .len(),
            4,
            "four corners, not five"
        );
        assert!(is_wire_closed(&model, &built.shape, T).unwrap());
    }

    #[test]
    fn an_open_polygon_is_one_edge_short_and_says_it_is_open() {
        let mut model = Model::new();
        let built = make_polygon(&mut model, &square(), false, T).unwrap();
        assert_eq!(
            explore_unique(&model, &built.shape, ShapeType::Edge)
                .unwrap()
                .len(),
            3
        );
        assert!(!is_wire_closed(&model, &built.shape, T).unwrap());
    }

    #[test]
    fn a_polygon_that_describes_nothing_is_refused() {
        let mut model = Model::new();
        assert!(make_polygon(&mut model, &[], false, T).is_err());
        assert!(make_polygon(&mut model, &square()[..1], false, T).is_err());
        assert!(
            make_polygon(&mut model, &square()[..2], true, T).is_err(),
            "two points do not enclose anything"
        );

        // A repeated point is a zero-length edge, which is not an edge.
        let mut doubled = square();
        doubled.insert(2, doubled[1]);
        assert!(make_polygon(&mut model, &doubled, true, T).is_err());

        // Repeating the first point at the end is refused either way. Closed,
        // it would add a zero-length segment; open, it would leave the wire
        // with two vertices in one place, which every later boundary walk
        // reads as a gap that happens to be zero wide.
        let mut wrapped = square();
        wrapped.push(wrapped[0]);
        for closed in [true, false] {
            let err = make_polygon(&mut model, &wrapped, closed, T).unwrap_err();
            assert!(
                err.to_string().contains("first and last points coincide"),
                "unexpected message: {err}"
            );
        }
    }

    #[test]
    fn a_planar_shape_reports_the_plane_it_lies_in() {
        let mut model = Model::new();
        let wire = make_polygon(&mut model, &square(), true, T).unwrap().shape;
        let plane = find_plane(&model, &wire, T)
            .unwrap()
            .expect("a flat square");
        assert!(
            plane.normal().is_parallel(ogeom_math::Direction::Z, T),
            "got {:?}",
            plane.normal()
        );
        for p in square() {
            assert!(plane.distance_to(p) < 1e-9);
        }
    }

    #[test]
    fn a_shape_that_is_not_planar_says_so_rather_than_fitting_one_anyway() {
        // Every set of three or more points has a best-fit plane, including a
        // set nowhere near one. Returning it would be a confident wrong answer.
        let mut model = Model::new();
        let mut skew = square();
        skew[2] = Point::new(2.0, 2.0, 1.0);
        let wire = make_polygon(&mut model, &skew, true, T).unwrap().shape;
        assert!(find_plane(&model, &wire, T).unwrap().is_none());
    }

    #[test]
    fn a_curved_edge_is_sampled_along_its_length_not_only_at_its_ends() {
        // An arc's endpoints lie in a great many planes the arc itself does
        // not. Checking only the vertices would call this shape planar.
        let mut model = Model::new();
        let circle = Circle::new(
            Frame::new(
                Point::ORIGIN,
                ogeom_math::Direction::Z,
                ogeom_math::Direction::X,
                T,
            )
            .unwrap(),
            2.0,
            T,
        )
        .unwrap();
        let arc = make_edge(
            &mut model,
            CircleCurve::new(circle).into(),
            (0.0, std::f64::consts::PI),
            T,
        )
        .unwrap()
        .shape;
        // In its own plane, it is planar.
        assert!(find_plane(&model, &arc, T).unwrap().is_some());

        // The two endpoints alone would admit the plane through them and the
        // z axis; the arc does not lie in it, and sampling catches that.
        let plane = find_plane(&model, &arc, T).unwrap().unwrap();
        assert!(plane.normal().is_parallel(ogeom_math::Direction::Z, T));
    }

    #[test]
    fn a_face_is_planar_when_its_surface_is() {
        let mut model = Model::new();
        let wire = make_polygon(&mut model, &square(), true, T).unwrap().shape;
        let face = make_face(
            &mut model,
            PlaneSurface::new(Plane::new(Frame::WORLD)).into(),
            std::slice::from_ref(&wire),
            T,
        )
        .unwrap()
        .shape;
        assert!(find_plane(&model, &face, T).unwrap().is_some());

        let solid = crate::make_box(&mut model, Frame::WORLD, (1.0, 1.0, 1.0), T)
            .unwrap()
            .shape;
        assert!(
            find_plane(&model, &solid, T).unwrap().is_none(),
            "a box is not planar"
        );
    }
}
#[cfg(test)]
#[allow(clippy::unwrap_used, clippy::expect_used)]
mod compsolid_tests {
    use super::*;
    use ogeom_geom::PlaneSurface;
    use ogeom_math::{Frame, Plane, Point};

    const T: Tolerances = Tolerances::millimetres();

    /// Two unit cubes glued along one shared face node at `x = 1`.
    fn glued_cubes(model: &mut Model) -> (Shape, Shape) {
        let a = crate::make_box(model, Frame::WORLD, (1.0, 1.0, 1.0), T)
            .unwrap()
            .shape;
        // The face of `a` at x = 1: the one whose every vertex has x = 1.
        let shared = ogeom_topo::explore(model, &a, ogeom_topo::Filter::OfType(ShapeType::Face))
            .unwrap()
            .into_iter()
            .find(|f| {
                ogeom_topo::explore(model, f, ogeom_topo::Filter::OfType(ShapeType::Vertex))
                    .unwrap()
                    .iter()
                    .all(|v| {
                        let p = model.node(v).unwrap().data().as_vertex().unwrap().point;
                        let world = v.transform(model.datums()).unwrap().apply(p);
                        (world.x - 1.0).abs() < 1e-9
                    })
            })
            .expect("the box has a face at x = 1");

        // The second cube: the shared face seen from the other side, plus
        // five new faces over the shared boundary vertices.
        let corner = |x: f64, y: f64, z: f64| Point::new(x, y, z);
        let quad = |model: &mut Model, points: [Point; 4]| -> Shape {
            let wire = crate::make_polygon(model, &points, true, T).unwrap().shape;
            let plane = crate::find_plane(model, &wire, T).unwrap().unwrap();
            make_face(
                model,
                PlaneSurface::over(Plane::new(plane.frame()), (-4.0, 4.0), (-4.0, 4.0))
                    .unwrap()
                    .into(),
                std::slice::from_ref(&wire),
                T,
            )
            .unwrap()
            .shape
        };
        let faces = [
            quad(
                model,
                [
                    corner(2.0, 0.0, 0.0),
                    corner(2.0, 1.0, 0.0),
                    corner(2.0, 1.0, 1.0),
                    corner(2.0, 0.0, 1.0),
                ],
            ),
            quad(
                model,
                [
                    corner(1.0, 0.0, 0.0),
                    corner(2.0, 0.0, 0.0),
                    corner(2.0, 0.0, 1.0),
                    corner(1.0, 0.0, 1.0),
                ],
            ),
            quad(
                model,
                [
                    corner(1.0, 1.0, 0.0),
                    corner(1.0, 1.0, 1.0),
                    corner(2.0, 1.0, 1.0),
                    corner(2.0, 1.0, 0.0),
                ],
            ),
            quad(
                model,
                [
                    corner(1.0, 0.0, 0.0),
                    corner(1.0, 1.0, 0.0),
                    corner(2.0, 1.0, 0.0),
                    corner(2.0, 0.0, 0.0),
                ],
            ),
            quad(
                model,
                [
                    corner(1.0, 0.0, 1.0),
                    corner(2.0, 0.0, 1.0),
                    corner(2.0, 1.0, 1.0),
                    corner(1.0, 1.0, 1.0),
                ],
            ),
        ];
        let mut shell_faces = vec![shared.reversed()];
        shell_faces.extend(faces);
        let shell = make_shell(model, &shell_faces).unwrap().shape;
        let b = make_solid(model, std::slice::from_ref(&shell))
            .unwrap()
            .shape;
        (a, b)
    }

    #[test]
    fn glued_solids_build_a_compsolid_and_loose_ones_are_refused() {
        let mut model = Model::new();
        let (a, b) = glued_cubes(&mut model);
        let built = make_compsolid(&mut model, &[a.clone(), b.clone()]).unwrap();
        assert_eq!(model.kind_of(&built.shape).unwrap(), ShapeType::CompSolid);
        assert_eq!(
            built.history.generated(&a),
            std::slice::from_ref(&built.shape)
        );

        // Two boxes merely sitting apart share nothing and are refused.
        let mut model = Model::new();
        let a = crate::make_box(&mut model, Frame::WORLD, (1.0, 1.0, 1.0), T)
            .unwrap()
            .shape;
        let far = Frame::new(
            Point::new(5.0, 0.0, 0.0),
            ogeom_math::Direction::Z,
            ogeom_math::Direction::X,
            T,
        )
        .unwrap();
        let b = crate::make_box(&mut model, far, (1.0, 1.0, 1.0), T)
            .unwrap()
            .shape;
        assert!(make_compsolid(&mut model, &[a.clone(), b]).is_err());
        assert!(make_compsolid(&mut model, std::slice::from_ref(&a)).is_err());
    }
}

#[cfg(test)]
#[allow(clippy::unwrap_used, clippy::expect_used)]
mod band_tests {
    use super::*;
    use approx::assert_relative_eq;
    use ogeom_core::Tolerances;
    use ogeom_geom::CircleCurve;
    use ogeom_math::{Circle, Cone, Frame, Sphere};

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

    fn mesh_area(model: &Model, face: &Shape) -> f64 {
        let fine = ogeom_mesh::Deflection {
            chord: 1e-3,
            ..ogeom_mesh::Deflection::default()
        };
        let mesh = ogeom_mesh::triangulate(model, face, fine, T).unwrap();
        mesh.triangles
            .iter()
            .map(|t| {
                let [a, b, c] = t.map(|i| mesh.positions[i as usize]);
                (b - a).cross(c - a).magnitude() / 2.0
            })
            .sum()
    }

    /// A full circle edge at `frame`'s origin, radius `r`, single vertex.
    fn ring(model: &mut Model, frame: Frame, r: f64) -> Shape {
        let circle = Circle::new(frame, r, T).unwrap();
        let vertex = make_vertex(model, frame.origin() + frame.x() * r).shape;
        make_edge_between(
            model,
            CircleCurve::new(circle).into(),
            (0.0, TAU),
            &vertex,
            &vertex,
            T,
        )
        .unwrap()
        .shape
    }

    /// The area a face's stored rings enclose in its surface's chart,
    /// walked as stored: positive where they keep the face on their left.
    fn chart_area(model: &Model, face: &Shape) -> f64 {
        use ogeom_geom::Curve2d as _;
        let surface = model.node(face).unwrap().data().as_face().unwrap().surface;
        let mut points = Vec::new();
        for wire in model
            .ordered_children_of(&face.oriented(ogeom_topo::Orientation::Forward))
            .unwrap()
        {
            for edge in model.ordered_children_of(&wire).unwrap() {
                let data = model.node(&edge).unwrap().data().as_edge().unwrap();
                let forward = edge.orientation() == ogeom_topo::Orientation::Forward;
                let (curve, range) = match data.pcurve_for(surface, edge.location()).unwrap() {
                    EdgeRepr::PCurve { curve, range, .. } => (*curve, *range),
                    EdgeRepr::Seam {
                        forward: f,
                        reversed: r,
                        range,
                        ..
                    } => (if forward { *f } else { *r }, *range),
                    _ => unreachable!("a pcurve was asked for"),
                };
                let pcurve = model.geometry().pcurve(curve).unwrap();
                for k in 0..32 {
                    let f = f64::from(k) / 32.0;
                    let t = if forward {
                        range.0 + (range.1 - range.0) * f
                    } else {
                        range.1 + (range.0 - range.1) * f
                    };
                    points.push(pcurve.point_at(t, T).unwrap());
                }
            }
        }
        (0..points.len())
            .map(|i| {
                let (a, b) = (points[i], points[(i + 1) % points.len()]);
                a.x * b.y - b.x * a.y
            })
            .sum::<f64>()
            / 2.0
    }

    /// A band keeps its face on the left of its ring in the surface's
    /// chart, whichever ring comes first and whichever way each circle
    /// turns: a drum's wall between two rims, a quarter turn of chart
    /// height two wide.
    #[test]
    fn a_band_keeps_its_face_on_the_left_whichever_way_its_rings_come() {
        let mut model = Model::new();
        let cylinder = ogeom_math::Cylinder::new(Frame::WORLD, 2.0, T).unwrap();
        let surface: SurfaceGeometry = ogeom_geom::CylinderSurface::new(cylinder, (-1.0, 3.0))
            .unwrap()
            .into();
        let at = |z: f64, up: bool| {
            let axis = if up {
                ogeom_math::Direction::Z
            } else {
                -ogeom_math::Direction::Z
            };
            Frame::new(Point::new(0.0, 0.0, z), axis, ogeom_math::Direction::X, T).unwrap()
        };
        for (first_up, second_up) in [(true, true), (false, true), (true, false)] {
            let low = ring(&mut model, at(0.0, first_up), 2.0);
            let high = ring(&mut model, at(2.0, second_up), 2.0);
            for (a, b) in [(&low, &high), (&high, &low)] {
                let face = make_revolution_band(&mut model, &surface, a, b, T).unwrap();
                let area = chart_area(&model, &face);
                assert_relative_eq!(area, TAU * 2.0, max_relative = 1e-9);
            }
        }
    }

    /// A band on a cylinder of either handedness, between rims spelled on
    /// frames of either handedness and either axis direction, each circle
    /// forward or reversed: each rim's pcurve lands on its circle at every
    /// parameter, and the face keeps its chart area `2π · 2` on its left.
    /// On a left-handed frame the angle turns about `-z`, and a reversed
    /// circle turns back, so a winding read from `z` alone runs such a
    /// rim's pcurve the wrong way round.
    #[test]
    fn a_band_follows_its_rims_on_frames_of_either_handedness() {
        use ogeom_geom::{Curve2d as _, Curve3d as _, Surface as _};
        let spellings = |f: Frame| {
            let (o, x, y, z) = (f.origin(), f.x(), f.y(), f.z());
            [
                f,
                f.mirrored(),
                Frame::from_axes(o, x, y, -z, T).unwrap(),
                Frame::from_axes(o, x, -y, -z, T).unwrap(),
            ]
        };
        for wall in spellings(Frame::WORLD) {
            let cylinder = ogeom_math::Cylinder::new(wall, 2.0, T).unwrap();
            let heights = if wall.z().vector().z > 0.0 {
                (-1.0, 3.0)
            } else {
                (-3.0, 1.0)
            };
            let surface: SurfaceGeometry = ogeom_geom::CylinderSurface::new(cylinder, heights)
                .unwrap()
                .into();
            let lows = spellings(Frame::WORLD);
            let highs = spellings(Frame::WORLD.with_origin(Point::new(0.0, 0.0, 2.0)));
            let pairs = lows.into_iter().zip(highs.into_iter().rev());
            for ((low, high), reversed) in pairs.flat_map(|p| [(p, false), (p, true)]) {
                let mut model = Model::new();
                let mut rim = |frame: Frame| {
                    use ogeom_geom::Reversible as _;
                    let circle: ogeom_geom::Curve =
                        CircleCurve::new(Circle::new(frame, 2.0, T).unwrap()).into();
                    let circle = if reversed { circle.reversed() } else { circle };
                    let vertex = make_vertex(&mut model, frame.origin() + frame.x() * 2.0).shape;
                    make_edge_between(&mut model, circle, (0.0, TAU), &vertex, &vertex, T)
                        .unwrap()
                        .shape
                };
                let (low, high) = (rim(low), rim(high));
                let face = make_revolution_band(&mut model, &surface, &low, &high, T).unwrap();
                let id = model.node(&face).unwrap().data().as_face().unwrap().surface;
                for rim in [&low, &high] {
                    let data = model.node(rim).unwrap().data().as_edge().unwrap();
                    let Some(EdgeRepr::Curve3d { curve, .. }) = data.curve3d() else {
                        unreachable!("a rim has a curve");
                    };
                    let curve = model.geometry().curve(*curve).unwrap();
                    let Some(EdgeRepr::PCurve { curve: p, .. }) =
                        data.pcurve_for(id, &Location::identity())
                    else {
                        unreachable!("a rim has a pcurve on the band");
                    };
                    let pcurve = model.geometry().pcurve(*p).unwrap();
                    for k in 0..=16 {
                        let t = TAU * f64::from(k) / 16.0;
                        let uv = pcurve.point_at(t, T).unwrap();
                        let lifted = surface.point_at(uv.x, uv.y, T).unwrap();
                        let on = curve.point_at(t, T).unwrap();
                        assert!(
                            lifted.distance(on) < 1e-9,
                            "wall {wall:?}: a rim's pcurve misses its circle by {}",
                            lifted.distance(on)
                        );
                    }
                }
                assert_relative_eq!(chart_area(&model, &face), TAU * 2.0, max_relative = 1e-9);
            }
        }
    }

    /// Two rims of a drum whose vertices stand half a turn apart: no seam
    /// along one column joins them, so the band is refused, and the rims
    /// come back with no pcurve on the surface. A pcurve left behind would
    /// start at the wrong vertex column, half a turn off its circle at
    /// every matched parameter, on a surface no face uses.
    #[test]
    fn a_refused_band_leaves_its_rings_without_pcurves() {
        let mut model = Model::new();
        let cylinder = ogeom_math::Cylinder::new(Frame::WORLD, 2.0, T).unwrap();
        let surface: SurfaceGeometry = ogeom_geom::CylinderSurface::new(cylinder, (-1.0, 3.0))
            .unwrap()
            .into();
        let at = |z: f64, x: ogeom_math::Direction| {
            Frame::new(Point::new(0.0, 0.0, z), ogeom_math::Direction::Z, x, T).unwrap()
        };
        let low = ring(&mut model, at(0.0, ogeom_math::Direction::X), 2.0);
        let high = ring(&mut model, at(2.0, -ogeom_math::Direction::X), 2.0);
        let refused = make_revolution_band(&mut model, &surface, &low, &high, T).unwrap_err();
        assert!(
            refused.to_string().contains("different columns"),
            "the refusal names why: {refused}"
        );
        for rim in [&low, &high] {
            let data = model.node(rim).unwrap().data().as_edge().unwrap();
            assert!(
                data.representations
                    .iter()
                    .all(|r| !matches!(r, EdgeRepr::PCurve { .. } | EdgeRepr::Seam { .. })),
                "a refused band left a pcurve on a rim"
            );
        }
    }

    /// A band of rings on a cylinder of either handedness, its rims'
    /// vertices a quarter turn round from the frame's `x`: the face's copy
    /// of the surface is turned so the seam stands at the start of its
    /// window, and every rim pcurve lies in the window and lands on its
    /// circle. On a left-handed frame the angle turns about `-z`, so a turn
    /// about `z` puts the seam half a turn off the window's start.
    #[test]
    fn a_band_of_rings_seams_at_its_window_start_on_either_handedness() {
        use ogeom_geom::{Curve2d as _, Curve3d as _, Surface as _};
        for wall in [Frame::WORLD, Frame::WORLD.mirrored()] {
            let mut model = Model::new();
            let cylinder = ogeom_math::Cylinder::new(wall, 2.0, T).unwrap();
            let surface: SurfaceGeometry = ogeom_geom::CylinderSurface::new(cylinder, (-1.0, 3.0))
                .unwrap()
                .into();
            let at = |z: f64| {
                Frame::new(
                    Point::new(0.0, 0.0, z),
                    ogeom_math::Direction::Z,
                    ogeom_math::Direction::Y,
                    T,
                )
                .unwrap()
            };
            let low = ring(&mut model, at(0.0), 2.0);
            let high = ring(&mut model, at(2.0), 2.0);
            let face = make_band_of_rings(
                &mut model,
                &surface,
                core::slice::from_ref(&low),
                core::slice::from_ref(&high),
                T,
            )
            .unwrap();
            let id = model.node(&face).unwrap().data().as_face().unwrap().surface;
            let band = model.geometry().surface(id).unwrap().clone();
            let ((ua, ub), _) = band.domain();
            for rim in [&low, &high] {
                let data = model.node(rim).unwrap().data().as_edge().unwrap();
                let Some(EdgeRepr::Curve3d { curve, .. }) = data.curve3d() else {
                    unreachable!("a rim has a curve");
                };
                let curve = model.geometry().curve(*curve).unwrap();
                let Some(EdgeRepr::PCurve { curve: p, .. }) =
                    data.pcurve_for(id, &Location::identity())
                else {
                    unreachable!("a rim has a pcurve on the band");
                };
                let pcurve = model.geometry().pcurve(*p).unwrap();
                for k in 0..=16 {
                    let t = TAU * f64::from(k) / 16.0;
                    let uv = pcurve.point_at(t, T).unwrap();
                    assert!(
                        uv.x >= ua - 1e-9 && uv.x <= ub + 1e-9,
                        "wall {wall:?}: u {} off the window",
                        uv.x
                    );
                    let lifted = band.point_at(uv.x, uv.y, T).unwrap();
                    assert!(lifted.distance(curve.point_at(t, T).unwrap()) < 1e-9);
                }
            }
            assert_relative_eq!(chart_area(&model, &face), TAU * 2.0, max_relative = 1e-9);
        }
    }

    /// A drum whose rims' vertices stand half a turn apart: the top rim is
    /// split on the bottom vertex's column, the band seams there, and the
    /// top cap takes both pieces, so the solid closes, checks valid and
    /// holds pi r^2 h. The history traces the top rim to its two pieces.
    #[test]
    fn a_rim_split_on_the_other_rims_column_closes_the_drum() {
        let mut model = Model::new();
        let (r, h) = (2.0, 3.0);
        let cylinder = ogeom_math::Cylinder::new(Frame::WORLD, r, T).unwrap();
        let surface: SurfaceGeometry = ogeom_geom::CylinderSurface::new(cylinder, (-1.0, 4.0))
            .unwrap()
            .into();
        let at = |z: f64, x: ogeom_math::Direction| {
            Frame::new(Point::new(0.0, 0.0, z), ogeom_math::Direction::Z, x, T).unwrap()
        };
        let low = ring(&mut model, at(0.0, ogeom_math::Direction::X), r);
        let high = ring(&mut model, at(h, -ogeom_math::Direction::X), r);
        let (anchor, _) = edge_vertices(&model, &low).unwrap().unwrap();

        let split = split_ring_on_column(&mut model, &surface, &high, &anchor, T)
            .unwrap()
            .expect("the rims stand on different columns");
        let pieces = model.ordered_children_of(&split.shape).unwrap();
        assert_eq!(pieces.len(), 2);
        assert_eq!(split.history.modified(&high).len(), 2);
        // The wire starts on the anchor's column: x = +r.
        let (start, _) = edge_vertices(&model, &pieces[0]).unwrap().unwrap();
        let start = model
            .node(&start)
            .unwrap()
            .data()
            .as_vertex()
            .unwrap()
            .point;
        assert!(start.distance(Point::new(r, 0.0, h)) <= T.confusion());
        // Aligned now: a second split has nothing to do.
        assert!(
            split_ring_on_column(&mut model, &surface, &low, &anchor, T)
                .unwrap()
                .is_none()
        );

        let band = make_band_of_rings(
            &mut model,
            &surface,
            core::slice::from_ref(&low),
            &pieces,
            T,
        )
        .unwrap();
        let down = Frame::new(
            Point::ORIGIN,
            -ogeom_math::Direction::Z,
            ogeom_math::Direction::X,
            T,
        )
        .unwrap();
        let bottom = make_face_with_pcurves(
            &mut model,
            ogeom_geom::PlaneSurface::new(ogeom_math::Plane::new(down)).into(),
            &[vec![low.reversed()]],
            T,
        )
        .unwrap()
        .shape;
        let top = make_face_with_pcurves(
            &mut model,
            ogeom_geom::PlaneSurface::new(ogeom_math::Plane::new(at(h, ogeom_math::Direction::X)))
                .into(),
            core::slice::from_ref(&pieces),
            T,
        )
        .unwrap()
        .shape;
        let shell = make_shell(&mut model, &[bottom, band, top]).unwrap().shape;
        assert!(is_shell_closed(&model, &shell).unwrap());
        let solid = make_solid(&mut model, &[shell]).unwrap().shape;
        let diagnosis = crate::check(&model, &solid, T).unwrap();
        assert!(diagnosis.is_valid(), "{:?}", diagnosis.problems);
        let fine = ogeom_mesh::Deflection {
            chord: 1e-3,
            ..ogeom_mesh::Deflection::default()
        };
        let volume = crate::volume_properties(&model, &solid, fine, T)
            .unwrap()
            .mass;
        assert_relative_eq!(
            volume,
            core::f64::consts::PI * r * r * h,
            max_relative = 1e-3
        );
    }

    #[test]
    fn a_cone_cap_takes_its_apex_as_a_degenerate_ring() {
        // Half angle 45 degrees, reference radius 2: apex at z = -2, slant
        // length 2*sqrt(2), lateral area pi * r * slant.
        let mut model = Model::new();
        let cone = Cone::new(Frame::WORLD, 2.0, core::f64::consts::FRAC_PI_4, T).unwrap();
        let surface: SurfaceGeometry = ogeom_geom::ConeSurface::new(cone, (-3.0, 3.0))
            .unwrap()
            .into();
        let rim = ring(&mut model, Frame::WORLD, 2.0);
        let face = crate::make_apex_band(&mut model, &surface, &rim, T).unwrap();
        assert_relative_eq!(
            mesh_area(&model, &face),
            core::f64::consts::PI * 2.0 * 2.0 * core::f64::consts::SQRT_2,
            max_relative = 1e-2
        );
    }

    #[test]
    fn a_ring_square_to_the_axis_is_not_a_band_ring() {
        // A button head: a sphere whose chart runs along z, bounded by two
        // circles cut square to the screw along x. Each lies on the sphere
        // and is closed, but neither is a parallel; building a band on
        // them hands every rim a latitude line it never follows. Refused,
        // and the reader's question answers the same.
        let mut model = Model::new();
        let sphere = ogeom_math::Sphere::centred(Point::ORIGIN, 5.0, T).unwrap();
        let surface: SurfaceGeometry = ogeom_geom::SphereSurface::new(sphere).into();
        // Rims in planes x = 3 and x = 4: radii 4 and 3, axis along x.
        let rim_at = |model: &mut Model, x: f64| {
            let frame = Frame::new(
                Point::new(x, 0.0, 0.0),
                ogeom_math::Direction::X,
                ogeom_math::Direction::Z,
                T,
            )
            .unwrap();
            ring(model, frame, (25.0 - x * x).sqrt())
        };
        let lo = rim_at(&mut model, 3.0);
        let hi = rim_at(&mut model, 4.0);
        assert!(!crate::rings_are_parallels(&model, &surface, &[&lo, &hi], T).unwrap());
        assert!(make_revolution_band(&mut model, &surface, &lo, &hi, T).is_err());
        // The genuine article still passes: parallels of the same sphere.
        let p_lo = ring(
            &mut model,
            Frame::new(
                Point::new(0.0, 0.0, 3.0),
                ogeom_math::Direction::Z,
                ogeom_math::Direction::X,
                T,
            )
            .unwrap(),
            4.0,
        );
        let p_hi = ring(
            &mut model,
            Frame::new(
                Point::new(0.0, 0.0, 4.0),
                ogeom_math::Direction::Z,
                ogeom_math::Direction::X,
                T,
            )
            .unwrap(),
            3.0,
        );
        assert!(crate::rings_are_parallels(&model, &surface, &[&p_lo, &p_hi], T).unwrap());
    }

    #[test]
    fn a_degenerate_ring_is_accepted_whichever_side_it_is_passed_on() {
        let mut model = Model::new();
        let cone = Cone::new(Frame::WORLD, 2.0, core::f64::consts::FRAC_PI_4, T).unwrap();
        let surface: SurfaceGeometry = ogeom_geom::ConeSurface::new(cone, (-3.0, 3.0))
            .unwrap()
            .into();
        let rim = ring(&mut model, Frame::WORLD, 2.0);
        let apex = make_vertex(&mut model, cone.apex()).shape;
        let mut data = EdgeData::new();
        data.degenerate = true;
        let degenerate = model.add_edge(data, &[apex.clone(), apex]).unwrap();
        // Degenerate first: the band swaps it into place rather than asking
        // the caller to know which ring is real.
        let face = make_revolution_band(&mut model, &surface, &degenerate, &rim, T).unwrap();
        assert_relative_eq!(
            mesh_area(&model, &face),
            core::f64::consts::PI * 2.0 * 2.0 * core::f64::consts::SQRT_2,
            max_relative = 1e-2
        );
    }

    #[test]
    fn a_sphere_cap_closes_against_its_pole() {
        // The equator ring and the north pole: a hemisphere, area 2 pi r^2.
        let mut model = Model::new();
        let sphere = Sphere::new(Frame::WORLD, 3.0, T).unwrap();
        let surface: SurfaceGeometry = ogeom_geom::SphereSurface::new(sphere).into();
        let rim = ring(&mut model, Frame::WORLD, 3.0);
        let pole = make_vertex(&mut model, Point::new(0.0, 0.0, 3.0)).shape;
        let mut data = EdgeData::new();
        data.degenerate = true;
        let degenerate = model.add_edge(data, &[pole.clone(), pole]).unwrap();
        let face = make_revolution_band(&mut model, &surface, &rim, &degenerate, T).unwrap();
        assert_relative_eq!(
            mesh_area(&model, &face),
            2.0 * core::f64::consts::PI * 9.0,
            max_relative = 1e-2
        );
    }

    #[test]
    fn two_degenerate_rings_bound_nothing_and_are_refused() {
        let mut model = Model::new();
        let cone = Cone::new(Frame::WORLD, 2.0, core::f64::consts::FRAC_PI_4, T).unwrap();
        let surface: SurfaceGeometry = ogeom_geom::ConeSurface::new(cone, (-3.0, 3.0))
            .unwrap()
            .into();
        let apex = make_vertex(&mut model, cone.apex()).shape;
        let mut data = EdgeData::new();
        data.degenerate = true;
        let a = model
            .add_edge(data.clone(), &[apex.clone(), apex.clone()])
            .unwrap();
        let b = model.add_edge(data, &[apex.clone(), apex]).unwrap();
        assert!(make_revolution_band(&mut model, &surface, &a, &b, T).is_err());
    }

    #[test]
    fn an_apex_band_needs_a_cone() {
        let mut model = Model::new();
        let sphere = Sphere::new(Frame::WORLD, 3.0, T).unwrap();
        let surface: SurfaceGeometry = ogeom_geom::SphereSurface::new(sphere).into();
        let rim = ring(&mut model, Frame::WORLD, 3.0);
        assert!(crate::make_apex_band(&mut model, &surface, &rim, T).is_err());
    }
}

#[cfg(test)]
#[allow(clippy::unwrap_used, clippy::expect_used)]
mod solid_tests {
    use super::*;
    use ogeom_math::Frame;

    const T: Tolerances = Tolerances::millimetres();

    /// The one shell of a box `size` on a side with its corner at `corner`.
    fn box_shell(model: &mut Model, corner: f64, size: f64) -> Shape {
        let frame = Frame::new(
            Point::new(corner, corner, corner),
            ogeom_math::Direction::Z,
            ogeom_math::Direction::X,
            T,
        )
        .unwrap();
        let solid = crate::make_box(model, frame, (size, size, size), T)
            .unwrap()
            .shape;
        model.children_of(&solid).unwrap()[0].clone()
    }

    /// A 10 mm cube with a 2 mm cubic void, the void given facing into
    /// its cavity or out of it and the outer shell facing out or in: the
    /// solid comes back with its outer shell facing out and its void
    /// facing in, valid, 992 mm3.
    #[test]
    fn a_solid_faces_out_and_its_void_faces_in_however_given() {
        for (outer_in, void_out) in [(false, false), (false, true), (true, false), (true, true)] {
            let mut model = Model::new();
            let outer = box_shell(&mut model, 0.0, 10.0);
            let void = box_shell(&mut model, 4.0, 2.0);
            let outer = if outer_in { outer.reversed() } else { outer };
            let void = if void_out { void } else { void.reversed() };
            let solid = make_solid(&mut model, &[outer, void]).unwrap().shape;
            let d = crate::check(&model, &solid, T).unwrap();
            assert!(d.is_valid(), "{outer_in} {void_out}: {d}");
            let shells = model.children_of(&solid).unwrap();
            assert!(enclosed_volume(&model, &shells[0]).unwrap() > 0.0);
            assert!(enclosed_volume(&model, &shells[1]).unwrap() < 0.0);
            let v = crate::volume_properties(&model, &solid, ogeom_mesh::Deflection::default(), T)
                .unwrap()
                .mass;
            assert!((v - 992.0).abs() < 1e-9, "{outer_in} {void_out}: {v}");
        }
    }
}