brepkit-operations 4.3.17

//! Parametric primitive shape builders.
//!
//! Provides factory functions for creating standard CAD primitives as
//! proper B-Rep solids: boxes, cylinders, cones, spheres, and tori.
//!
//! Each function returns a [`SolidId`] for a newly created solid in
//! the given [`Topology`]. The solids are built with correct manifold
//! topology and outward-pointing face normals.

use std::collections::{HashMap, HashSet};
use std::f64::consts::{FRAC_PI_2, PI, TAU};

use brepkit_math::tolerance::Tolerance;
use brepkit_math::vec::{Point3, Vec3};
use brepkit_topology::Topology;
use brepkit_topology::edge::{Edge, EdgeCurve};
use brepkit_topology::face::{Face, FaceSurface};
use brepkit_topology::shell::Shell;
use brepkit_topology::solid::{Solid, SolidId};
use brepkit_topology::vertex::Vertex;
use brepkit_topology::wire::{OrientedEdge, Wire};

/// Create a box solid with one corner at the origin.
///
/// The box extends from `(0, 0, 0)` to `(dx, dy, dz)`.
/// One corner sits at the origin; the opposite at `(dx, dy, dz)`.
///
/// # Errors
///
/// Returns an error if any dimension is zero or negative.
pub fn make_box(
    topo: &mut Topology,
    dx: f64,
    dy: f64,
    dz: f64,
) -> Result<SolidId, crate::OperationsError> {
    let tol = Tolerance::new();

    if dx <= tol.linear || dy <= tol.linear || dz <= tol.linear {
        return Err(crate::OperationsError::InvalidInput {
            reason: format!("box dimensions must be positive, got ({dx}, {dy}, {dz})"),
        });
    }

    let v = [
        topo.add_vertex(Vertex::new(Point3::new(0.0, 0.0, 0.0), tol.linear)),
        topo.add_vertex(Vertex::new(Point3::new(dx, 0.0, 0.0), tol.linear)),
        topo.add_vertex(Vertex::new(Point3::new(dx, dy, 0.0), tol.linear)),
        topo.add_vertex(Vertex::new(Point3::new(0.0, dy, 0.0), tol.linear)),
        topo.add_vertex(Vertex::new(Point3::new(0.0, 0.0, dz), tol.linear)),
        topo.add_vertex(Vertex::new(Point3::new(dx, 0.0, dz), tol.linear)),
        topo.add_vertex(Vertex::new(Point3::new(dx, dy, dz), tol.linear)),
        topo.add_vertex(Vertex::new(Point3::new(0.0, dy, dz), tol.linear)),
    ];

    // Bottom ring (z=0)
    let eb0 = topo.add_edge(Edge::new(v[0], v[1], EdgeCurve::Line));
    let eb1 = topo.add_edge(Edge::new(v[1], v[2], EdgeCurve::Line));
    let eb2 = topo.add_edge(Edge::new(v[2], v[3], EdgeCurve::Line));
    let eb3 = topo.add_edge(Edge::new(v[3], v[0], EdgeCurve::Line));
    // Top ring (z=dz)
    let et0 = topo.add_edge(Edge::new(v[4], v[5], EdgeCurve::Line));
    let et1 = topo.add_edge(Edge::new(v[5], v[6], EdgeCurve::Line));
    let et2 = topo.add_edge(Edge::new(v[6], v[7], EdgeCurve::Line));
    let et3 = topo.add_edge(Edge::new(v[7], v[4], EdgeCurve::Line));
    // Verticals
    let ev0 = topo.add_edge(Edge::new(v[0], v[4], EdgeCurve::Line));
    let ev1 = topo.add_edge(Edge::new(v[1], v[5], EdgeCurve::Line));
    let ev2 = topo.add_edge(Edge::new(v[2], v[6], EdgeCurve::Line));
    let ev3 = topo.add_edge(Edge::new(v[3], v[7], EdgeCurve::Line));

    let mk_face = |topo: &mut Topology,
                   edges: [(brepkit_topology::edge::EdgeId, bool); 4],
                   normal: Vec3,
                   d: f64|
     -> Result<brepkit_topology::face::FaceId, crate::OperationsError> {
        let wire = Wire::new(
            edges
                .iter()
                .map(|&(eid, fwd)| OrientedEdge::new(eid, fwd))
                .collect(),
            true,
        )
        .map_err(crate::OperationsError::Topology)?;
        let wid = topo.add_wire(wire);
        Ok(topo.add_face(Face::new(wid, vec![], FaceSurface::Plane { normal, d })))
    };

    // Plane 'd' values: signed distance from origin along normal
    let bottom = mk_face(
        topo,
        [(eb0, false), (eb3, false), (eb2, false), (eb1, false)],
        Vec3::new(0.0, 0.0, -1.0),
        0.0,
    )?;
    let top = mk_face(
        topo,
        [(et0, true), (et1, true), (et2, true), (et3, true)],
        Vec3::new(0.0, 0.0, 1.0),
        dz,
    )?;
    let front = mk_face(
        topo,
        [(eb0, true), (ev1, true), (et0, false), (ev0, false)],
        Vec3::new(0.0, -1.0, 0.0),
        0.0,
    )?;
    let back = mk_face(
        topo,
        [(eb2, true), (ev3, true), (et2, false), (ev2, false)],
        Vec3::new(0.0, 1.0, 0.0),
        dy,
    )?;
    let left = mk_face(
        topo,
        [(eb3, true), (ev0, true), (et3, false), (ev3, false)],
        Vec3::new(-1.0, 0.0, 0.0),
        0.0,
    )?;
    let right = mk_face(
        topo,
        [(eb1, true), (ev2, true), (et1, false), (ev1, false)],
        Vec3::new(1.0, 0.0, 0.0),
        dx,
    )?;

    let shell = Shell::new(vec![bottom, top, front, back, left, right])
        .map_err(crate::OperationsError::Topology)?;
    let shell_id = topo.add_shell(shell);
    Ok(topo.add_solid(Solid::new(shell_id, vec![])))
}

/// Create a cylinder solid with its axis along +Z, base at the origin.
///
/// The cylinder extends from `z = 0` to `z = height`.
/// Built with one `CylindricalSurface` lateral face and two planar cap faces.
///
/// # Errors
///
/// Returns an error if `radius` or `height` is zero or negative.
pub fn make_cylinder(
    topo: &mut Topology,
    radius: f64,
    height: f64,
) -> Result<SolidId, crate::OperationsError> {
    let tol = Tolerance::new();

    if radius <= tol.linear {
        return Err(crate::OperationsError::InvalidInput {
            reason: format!("cylinder radius must be positive, got {radius}"),
        });
    }
    if height <= tol.linear {
        return Err(crate::OperationsError::InvalidInput {
            reason: format!("cylinder height must be positive, got {height}"),
        });
    }

    // Cylinder base at z=0, top at z=height.
    // (brepjs drill and placement code assumes this convention.)

    let cyl_surface = brepkit_math::surfaces::CylindricalSurface::new(
        Point3::new(0.0, 0.0, 0.0),
        Vec3::new(0.0, 0.0, 1.0),
        radius,
    )
    .map_err(crate::OperationsError::Math)?;

    // --- Lateral face: single face with degenerate seam wire ---
    let v_bot = topo.add_vertex(Vertex::new(Point3::new(radius, 0.0, 0.0), tol.linear));
    let v_top = topo.add_vertex(Vertex::new(Point3::new(radius, 0.0, height), tol.linear));

    let bot_circle = brepkit_math::curves::Circle3D::new(
        Point3::new(0.0, 0.0, 0.0),
        Vec3::new(0.0, 0.0, 1.0),
        radius,
    )
    .map_err(crate::OperationsError::Math)?;
    let top_circle = brepkit_math::curves::Circle3D::new(
        Point3::new(0.0, 0.0, height),
        Vec3::new(0.0, 0.0, 1.0),
        radius,
    )
    .map_err(crate::OperationsError::Math)?;

    let e_bot_circle = topo.add_edge(Edge::new(v_bot, v_bot, EdgeCurve::Circle(bot_circle)));
    let e_top_circle = topo.add_edge(Edge::new(v_top, v_top, EdgeCurve::Circle(top_circle)));
    let e_seam = topo.add_edge(Edge::new(v_bot, v_top, EdgeCurve::Line));

    let lateral_wire = Wire::new(
        vec![
            OrientedEdge::new(e_bot_circle, true),
            OrientedEdge::new(e_seam, true),
            OrientedEdge::new(e_top_circle, false),
            OrientedEdge::new(e_seam, false),
        ],
        true,
    )
    .map_err(crate::OperationsError::Topology)?;
    let lateral_wid = topo.add_wire(lateral_wire);
    let lateral_face = topo.add_face(Face::new(
        lateral_wid,
        vec![],
        FaceSurface::Cylinder(cyl_surface),
    ));

    // --- Bottom cap (z = 0, normal pointing down) ---
    // Reuse the same circle edge as the lateral face for watertight topology.
    // Reversed orientation: CW from +z corresponds to outward normal -z.
    let bot_cap_wire = Wire::new(vec![OrientedEdge::new(e_bot_circle, false)], true)
        .map_err(crate::OperationsError::Topology)?;
    let bot_wid = topo.add_wire(bot_cap_wire);
    let bot_face = topo.add_face(Face::new(
        bot_wid,
        vec![],
        FaceSurface::Plane {
            normal: Vec3::new(0.0, 0.0, -1.0),
            d: 0.0,
        },
    ));

    // --- Top cap (z = height, normal pointing up) ---
    // Reuse the same circle edge; forward orientation gives outward normal +z.
    let top_cap_wire = Wire::new(vec![OrientedEdge::new(e_top_circle, true)], true)
        .map_err(crate::OperationsError::Topology)?;
    let top_wid = topo.add_wire(top_cap_wire);
    let top_face = topo.add_face(Face::new(
        top_wid,
        vec![],
        FaceSurface::Plane {
            normal: Vec3::new(0.0, 0.0, 1.0),
            d: height,
        },
    ));

    let shell = Shell::new(vec![lateral_face, bot_face, top_face])
        .map_err(crate::OperationsError::Topology)?;
    let shell_id = topo.add_shell(shell);
    Ok(topo.add_solid(Solid::new(shell_id, vec![])))
}

/// Create a cone solid with its base at the origin, axis along +Z.
///
/// The cone has `bottom_radius` at `z = 0` and `top_radius` at
/// `z = height`. Setting `top_radius = 0` creates a pointed cone;
/// setting it to a positive value creates a truncated cone (frustum).
/// The base is at the origin.
///
/// # Errors
///
/// Returns an error if `height` is non-positive, both radii are zero,
/// or any radius is negative.
#[allow(clippy::too_many_lines)]
pub fn make_cone(
    topo: &mut Topology,
    bottom_radius: f64,
    top_radius: f64,
    height: f64,
) -> Result<SolidId, crate::OperationsError> {
    let tol = Tolerance::new();

    if height <= tol.linear {
        return Err(crate::OperationsError::InvalidInput {
            reason: format!("cone height must be positive, got {height}"),
        });
    }
    if bottom_radius < 0.0 || top_radius < 0.0 {
        return Err(crate::OperationsError::InvalidInput {
            reason: format!(
                "cone radii must be non-negative, got bottom={bottom_radius}, top={top_radius}"
            ),
        });
    }
    if bottom_radius <= tol.linear && top_radius <= tol.linear {
        return Err(crate::OperationsError::InvalidInput {
            reason: "cone must have at least one non-zero radius".into(),
        });
    }

    // Determine which end is larger and compute virtual apex + half-angle
    // Base at z=0, top at z=height.
    let (r_big, r_small, big_z, small_z, axis_sign) = if bottom_radius >= top_radius {
        (bottom_radius, top_radius, 0.0, height, -1.0_f64)
    } else {
        (top_radius, bottom_radius, height, 0.0, 1.0_f64)
    };

    // The ConicalSurface formula is:
    //   P(u,v) = apex + v*(cos(a)*radial(u) + sin(a)*axis)
    // where `a` is the angle from the radial plane to the surface generator.
    // For a cone with base radius R at axial distance H from the apex:
    //   tan(a) = H / R  →  a = atan2(H, R)
    let half_angle = if r_small <= tol.linear {
        // Pointed cone: apex at the small end
        height.atan2(r_big)
    } else {
        // Frustum: virtual apex beyond the small end
        let axial_to_apex = r_small * height / (r_big - r_small);
        (axial_to_apex + height).atan2(r_big)
    };

    // half_angle must be in (0, π/2) for ConicalSurface
    if half_angle <= tol.angular || half_angle >= FRAC_PI_2 {
        // Degenerate case — fall back to revolve approach
        let face_id = make_trapezoid_xz_face(topo, bottom_radius, top_radius, 0.0, height)?;
        return crate::revolve::revolve(
            topo,
            face_id,
            Point3::new(0.0, 0.0, 0.0),
            Vec3::new(0.0, 0.0, 1.0),
            2.0 * PI,
        );
    }

    let apex_pos = if r_small <= tol.linear {
        Point3::new(0.0, 0.0, small_z)
    } else {
        let axial_to_apex = r_small * height / (r_big - r_small);
        // Apex is beyond the small end, away from the big end.
        // axis_sign points big→small, so we negate it to go small→apex.
        Point3::new(0.0, 0.0, small_z - axis_sign * axial_to_apex)
    };

    // Axis points from apex toward the base (big end), so that the
    // surface generator v>0 sweeps from apex outward to the base.
    let axis_dir = Vec3::new(0.0, 0.0, axis_sign);
    let cone_surface = brepkit_math::surfaces::ConicalSurface::new(apex_pos, axis_dir, half_angle)
        .map_err(crate::OperationsError::Math)?;

    let mut faces = Vec::new();

    // --- Lateral conical face ---
    if r_small <= tol.linear {
        // Pointed cone: degenerate wire from base circle to apex
        let v_apex = topo.add_vertex(Vertex::new(apex_pos, tol.linear));
        let v_base = topo.add_vertex(Vertex::new(Point3::new(r_big, 0.0, big_z), tol.linear));

        let base_circle = brepkit_math::curves::Circle3D::new(
            Point3::new(0.0, 0.0, big_z),
            Vec3::new(0.0, 0.0, 1.0),
            r_big,
        )
        .map_err(crate::OperationsError::Math)?;
        let e_circle = topo.add_edge(Edge::new(v_base, v_base, EdgeCurve::Circle(base_circle)));
        let e_seam = topo.add_edge(Edge::new(v_base, v_apex, EdgeCurve::Line));

        let lateral_wire = Wire::new(
            vec![
                OrientedEdge::new(e_circle, true),
                OrientedEdge::new(e_seam, true),
                OrientedEdge::new(e_seam, false),
            ],
            true,
        )
        .map_err(crate::OperationsError::Topology)?;
        let lateral_wid = topo.add_wire(lateral_wire);
        faces.push(topo.add_face(Face::new(
            lateral_wid,
            vec![],
            FaceSurface::Cone(cone_surface),
        )));

        // Base cap: reuse the same circle edge for watertight topology.
        // axis_sign < 0 means bottom is bigger (cap at z=0, outward normal -z, reversed edge).
        // axis_sign > 0 means top is bigger (cap at z=height, outward normal +z, forward edge).
        let cap_forward = axis_sign > 0.0;
        let cap_wire = Wire::new(vec![OrientedEdge::new(e_circle, cap_forward)], true)
            .map_err(crate::OperationsError::Topology)?;
        let cap_wid = topo.add_wire(cap_wire);
        let cap_normal = Vec3::new(0.0, 0.0, axis_sign);
        faces.push(topo.add_face(Face::new(
            cap_wid,
            vec![],
            FaceSurface::Plane {
                normal: cap_normal,
                d: big_z,
            },
        )));
    } else {
        // Frustum: two circles
        let v_bot = topo.add_vertex(Vertex::new(
            Point3::new(bottom_radius, 0.0, 0.0),
            tol.linear,
        ));
        let v_top = topo.add_vertex(Vertex::new(
            Point3::new(top_radius, 0.0, height),
            tol.linear,
        ));

        let bot_circle = brepkit_math::curves::Circle3D::new(
            Point3::new(0.0, 0.0, 0.0),
            Vec3::new(0.0, 0.0, 1.0),
            bottom_radius,
        )
        .map_err(crate::OperationsError::Math)?;
        let top_circle = brepkit_math::curves::Circle3D::new(
            Point3::new(0.0, 0.0, height),
            Vec3::new(0.0, 0.0, 1.0),
            top_radius,
        )
        .map_err(crate::OperationsError::Math)?;
        let e_bot = topo.add_edge(Edge::new(v_bot, v_bot, EdgeCurve::Circle(bot_circle)));
        let e_top = topo.add_edge(Edge::new(v_top, v_top, EdgeCurve::Circle(top_circle)));
        let e_seam = topo.add_edge(Edge::new(v_bot, v_top, EdgeCurve::Line));

        let lateral_wire = Wire::new(
            vec![
                OrientedEdge::new(e_bot, true),
                OrientedEdge::new(e_seam, true),
                OrientedEdge::new(e_top, false),
                OrientedEdge::new(e_seam, false),
            ],
            true,
        )
        .map_err(crate::OperationsError::Topology)?;
        let lateral_wid = topo.add_wire(lateral_wire);
        faces.push(topo.add_face(Face::new(
            lateral_wid,
            vec![],
            FaceSurface::Cone(cone_surface),
        )));

        // Bottom cap (z=0): reuse the same circle edge for watertight topology.
        let bot_cap_wire = Wire::new(vec![OrientedEdge::new(e_bot, false)], true)
            .map_err(crate::OperationsError::Topology)?;
        let bot_wid = topo.add_wire(bot_cap_wire);
        faces.push(topo.add_face(Face::new(
            bot_wid,
            vec![],
            FaceSurface::Plane {
                normal: Vec3::new(0.0, 0.0, -1.0),
                d: 0.0,
            },
        )));

        // Top cap (z=height): reuse the same circle edge for watertight topology.
        let top_cap_wire = Wire::new(vec![OrientedEdge::new(e_top, true)], true)
            .map_err(crate::OperationsError::Topology)?;
        let top_wid = topo.add_wire(top_cap_wire);
        faces.push(topo.add_face(Face::new(
            top_wid,
            vec![],
            FaceSurface::Plane {
                normal: Vec3::new(0.0, 0.0, 1.0),
                d: height,
            },
        )));
    }

    let shell = Shell::new(faces).map_err(crate::OperationsError::Topology)?;
    let shell_id = topo.add_shell(shell);
    Ok(topo.add_solid(Solid::new(shell_id, vec![])))
}

/// Create a sphere solid centered at the origin.
///
/// Built as a single spherical face with a degenerate boundary wire,
/// boundary wire, yielding exact `SphericalSurface` geometry with proper
/// topology for boolean operations.
///
/// The `segments` parameter controls the number of vertices on the
/// equatorial boundary polygon (minimum 4). Higher values improve
/// boolean accuracy for curved intersections. Tessellation density
/// is separately controlled by the `deflection` parameter at
/// tessellation time.
///
/// # Errors
///
/// Returns an error if `radius` is zero or negative, or `segments < 4`.
#[allow(clippy::too_many_lines)]
pub fn make_sphere(
    topo: &mut Topology,
    radius: f64,
    segments: usize,
) -> Result<SolidId, crate::OperationsError> {
    let tol = Tolerance::new();

    if radius <= tol.linear {
        return Err(crate::OperationsError::InvalidInput {
            reason: format!("sphere radius must be positive, got {radius}"),
        });
    }
    if segments < 4 {
        return Err(crate::OperationsError::InvalidInput {
            reason: format!("sphere needs at least 4 segments, got {segments}"),
        });
    }

    let surface_n =
        brepkit_math::surfaces::SphericalSurface::new(Point3::new(0.0, 0.0, 0.0), radius)
            .map_err(crate::OperationsError::Math)?;
    let surface_s =
        brepkit_math::surfaces::SphericalSurface::new(Point3::new(0.0, 0.0, 0.0), radius)
            .map_err(crate::OperationsError::Math)?;

    // Equatorial polygon: `segments` vertices evenly spaced on the circle
    // at z = 0, forming the shared boundary between north and south hemispheres.
    let eq_verts: Vec<_> = (0..segments)
        .map(|i| {
            let theta = TAU * i as f64 / segments as f64;
            let pt = Point3::new(radius * theta.cos(), radius * theta.sin(), 0.0);
            topo.add_vertex(Vertex::new(pt, tol.linear))
        })
        .collect();

    // Line edges connecting consecutive equatorial vertices (closed polygon).
    let eq_edges: Vec<_> = (0..segments)
        .map(|i| {
            let j = (i + 1) % segments;
            topo.add_edge(Edge::new(eq_verts[i], eq_verts[j], EdgeCurve::Line))
        })
        .collect();

    // North hemisphere: equatorial wire traversed forward (CCW from +Z).
    // Outward normal points upward (+Z hemisphere).
    let north_wire = Wire::new(
        eq_edges
            .iter()
            .map(|&eid| OrientedEdge::new(eid, true))
            .collect(),
        true,
    )
    .map_err(crate::OperationsError::Topology)?;
    let north_wid = topo.add_wire(north_wire);
    let north_face = topo.add_face(Face::new(north_wid, vec![], FaceSurface::Sphere(surface_n)));

    // South hemisphere: equatorial wire traversed backward (CW from +Z).
    // Outward normal points downward (-Z hemisphere).
    let south_wire = Wire::new(
        eq_edges
            .iter()
            .rev()
            .map(|&eid| OrientedEdge::new(eid, false))
            .collect(),
        true,
    )
    .map_err(crate::OperationsError::Topology)?;
    let south_wid = topo.add_wire(south_wire);
    let south_face = topo.add_face(Face::new(south_wid, vec![], FaceSurface::Sphere(surface_s)));

    let shell =
        Shell::new(vec![north_face, south_face]).map_err(crate::OperationsError::Topology)?;
    let shell_id = topo.add_shell(shell);
    Ok(topo.add_solid(Solid::new(shell_id, vec![])))
}

/// Create a torus solid centered at the origin in the XY plane.
///
/// Built as a single `ToroidalSurface` face with exact analytic geometry,
/// yielding fast tessellation.
///
/// The `segments` parameter is accepted for API compatibility but ignored —
/// tessellation density is controlled by the `deflection` parameter.
///
/// # Errors
///
/// Returns an error if either radius is non-positive, or if the minor
/// radius is greater than the major radius (self-intersecting torus).
pub fn make_torus(
    topo: &mut Topology,
    major_radius: f64,
    minor_radius: f64,
    segments: usize,
) -> Result<SolidId, crate::OperationsError> {
    let tol = Tolerance::new();

    if major_radius <= tol.linear {
        return Err(crate::OperationsError::InvalidInput {
            reason: format!("torus major radius must be positive, got {major_radius}"),
        });
    }
    if minor_radius <= tol.linear {
        return Err(crate::OperationsError::InvalidInput {
            reason: format!("torus minor radius must be positive, got {minor_radius}"),
        });
    }
    if minor_radius >= major_radius {
        return Err(crate::OperationsError::InvalidInput {
            reason: format!(
                "torus minor radius ({minor_radius}) must be less than major radius ({major_radius})"
            ),
        });
    }
    if segments < 4 {
        return Err(crate::OperationsError::InvalidInput {
            reason: format!("torus needs at least 4 segments, got {segments}"),
        });
    }

    let surface = brepkit_math::surfaces::ToroidalSurface::new(
        Point3::new(0.0, 0.0, 0.0),
        major_radius,
        minor_radius,
    )
    .map_err(crate::OperationsError::Math)?;

    // The torus surface is doubly periodic (two seam curves).
    // Minimal CW complex: 1 vertex, 2 edges, 1 face → V-E+F = 0 (genus 1).
    // The boundary wire follows the fundamental polygon: a → b → a⁻¹ → b⁻¹.
    let v0 = topo.add_vertex(Vertex::new(
        Point3::new(major_radius + minor_radius, 0.0, 0.0),
        tol.linear,
    ));
    // Seam edge a (longitudinal — around the tube)
    let ea = topo.add_edge(Edge::new(v0, v0, EdgeCurve::Line));
    // Seam edge b (meridional — around the ring)
    let eb = topo.add_edge(Edge::new(v0, v0, EdgeCurve::Line));

    let wire = Wire::new(
        vec![
            OrientedEdge::new(ea, true),
            OrientedEdge::new(eb, true),
            OrientedEdge::new(ea, false),
            OrientedEdge::new(eb, false),
        ],
        true,
    )
    .map_err(crate::OperationsError::Topology)?;
    let wid = topo.add_wire(wire);

    let face_id = topo.add_face(Face::new(wid, vec![], FaceSurface::Torus(surface)));

    let shell = Shell::new(vec![face_id]).map_err(crate::OperationsError::Topology)?;
    let shell_id = topo.add_shell(shell);
    Ok(topo.add_solid(Solid::new(shell_id, vec![])))
}

/// Build a trapezoid face in the XZ plane (y=0) for cone profiles.
///
/// Vertices: `(0,-hz)`, `(bottom_r,-hz)`, `(top_r,+hz)`, `(0,+hz)`
/// CCW winding when viewed from -Y.
fn make_trapezoid_xz_face(
    topo: &mut Topology,
    bottom_radius: f64,
    top_radius: f64,
    z_bottom: f64,
    z_top: f64,
) -> Result<brepkit_topology::face::FaceId, crate::OperationsError> {
    let tol = Tolerance::new();

    let v0 = topo.add_vertex(Vertex::new(Point3::new(0.0, 0.0, z_bottom), tol.linear));
    let v1 = topo.add_vertex(Vertex::new(
        Point3::new(bottom_radius, 0.0, z_bottom),
        tol.linear,
    ));
    let v2 = topo.add_vertex(Vertex::new(Point3::new(top_radius, 0.0, z_top), tol.linear));
    let v3 = topo.add_vertex(Vertex::new(Point3::new(0.0, 0.0, z_top), tol.linear));

    let e0 = topo.add_edge(Edge::new(v0, v1, EdgeCurve::Line));
    let e1 = topo.add_edge(Edge::new(v1, v2, EdgeCurve::Line));
    let e2 = topo.add_edge(Edge::new(v2, v3, EdgeCurve::Line));
    let e3 = topo.add_edge(Edge::new(v3, v0, EdgeCurve::Line));

    let wire = Wire::new(
        vec![
            OrientedEdge::new(e0, true),
            OrientedEdge::new(e1, true),
            OrientedEdge::new(e2, true),
            OrientedEdge::new(e3, true),
        ],
        true,
    )
    .map_err(crate::OperationsError::Topology)?;
    let wid = topo.add_wire(wire);

    let normal = Vec3::new(0.0, -1.0, 0.0);
    Ok(topo.add_face(Face::new(
        wid,
        vec![],
        FaceSurface::Plane { normal, d: 0.0 },
    )))
}

#[allow(dead_code)]
/// Build a rectangular face in the XZ plane (y=0) from corners.
fn make_rect_xz_face(
    topo: &mut Topology,
    x0: f64,
    z0: f64,
    x1: f64,
    z1: f64,
) -> Result<brepkit_topology::face::FaceId, crate::OperationsError> {
    let tol = Tolerance::new();

    let v0 = topo.add_vertex(Vertex::new(Point3::new(x0, 0.0, z0), tol.linear));
    let v1 = topo.add_vertex(Vertex::new(Point3::new(x1, 0.0, z0), tol.linear));
    let v2 = topo.add_vertex(Vertex::new(Point3::new(x1, 0.0, z1), tol.linear));
    let v3 = topo.add_vertex(Vertex::new(Point3::new(x0, 0.0, z1), tol.linear));

    let e0 = topo.add_edge(Edge::new(v0, v1, EdgeCurve::Line));
    let e1 = topo.add_edge(Edge::new(v1, v2, EdgeCurve::Line));
    let e2 = topo.add_edge(Edge::new(v2, v3, EdgeCurve::Line));
    let e3 = topo.add_edge(Edge::new(v3, v0, EdgeCurve::Line));

    let wire = Wire::new(
        vec![
            OrientedEdge::new(e0, true),
            OrientedEdge::new(e1, true),
            OrientedEdge::new(e2, true),
            OrientedEdge::new(e3, true),
        ],
        true,
    )
    .map_err(crate::OperationsError::Topology)?;
    let wid = topo.add_wire(wire);

    let normal = Vec3::new(0.0, -1.0, 0.0);
    let d = 0.0;
    Ok(topo.add_face(Face::new(wid, vec![], FaceSurface::Plane { normal, d })))
}

/// Create a solid from the 3D convex hull of a point cloud.
///
/// Uses the Quickhull algorithm to compute the convex hull, then
/// converts the resulting triangulated surface to a B-Rep solid with one
/// planar face per facet, coplanar triangles merged into one polygon.
///
/// # Errors
///
/// Returns an error if fewer than 4 non-coplanar points are provided.
pub fn make_convex_hull(
    topo: &mut Topology,
    points: &[Point3],
) -> Result<SolidId, crate::OperationsError> {
    fn root(facet: &mut [usize], mut i: usize) -> usize {
        while facet[i] != i {
            facet[i] = facet[facet[i]];
            i = facet[i];
        }
        i
    }
    let hull = brepkit_math::convex_hull::convex_hull_3d(points).ok_or_else(|| {
        crate::OperationsError::InvalidInput {
            reason: "points are coplanar or degenerate — cannot form a 3D convex hull".into(),
        }
    })?;

    let tol = Tolerance::new();

    // Each triangle's plane and twice its area, and the facets: triangles
    // joined across a shared edge where each one's far corner lies on the
    // other's plane.
    let plane_of = |&[a, b, c]: &[usize; 3]| {
        let (pa, pb, pc) = (hull.vertices[a], hull.vertices[b], hull.vertices[c]);
        let cross = (pb - pa).cross(pc - pa);
        let normal = cross.normalize().unwrap_or(Vec3::new(0.0, 0.0, 1.0));
        let d = normal
            .x()
            .mul_add(pa.x(), normal.y().mul_add(pa.y(), normal.z() * pa.z()));
        (normal, d, cross.length())
    };
    let planes: Vec<(Vec3, f64, f64)> = hull.faces.iter().map(plane_of).collect();
    let scale = hull
        .vertices
        .iter()
        .map(|p| p.x().abs().max(p.y().abs()).max(p.z().abs()))
        .fold(1.0_f64, f64::max);
    let on_plane = |t: usize, v: usize| {
        let ((n, d, _), p) = (planes[t], hull.vertices[v]);
        (n.x().mul_add(p.x(), n.y().mul_add(p.y(), n.z() * p.z())) - d).abs() <= 1e-9 * scale
    };
    let far_corner = |t: usize, a: usize, b: usize| {
        hull.faces[t]
            .into_iter()
            .find(|&v| v != a && v != b)
            .unwrap_or(a)
    };
    let mut facet: Vec<usize> = (0..hull.faces.len()).collect();
    let mut by_edge: HashMap<(usize, usize), usize> = HashMap::new();
    for (t, tri) in hull.faces.iter().enumerate() {
        for k in 0..3 {
            let (a, b) = (tri[k], tri[(k + 1) % 3]);
            if let Some(&other) = by_edge.get(&(a.min(b), a.max(b))) {
                if planes[t].0.dot(planes[other].0) > 0.0
                    && on_plane(t, far_corner(other, a, b))
                    && on_plane(other, far_corner(t, a, b))
                {
                    let (r1, r2) = (root(&mut facet, t), root(&mut facet, other));
                    facet[r1] = r2;
                }
            } else {
                by_edge.insert((a.min(b), a.max(b)), t);
            }
        }
    }
    // Each facet's boundary: its triangles' edges less those two of them
    // share, chained into one loop (a convex polygon), on the plane of its
    // largest triangle. A facet whose corners stray off that plane, merged
    // through a chain of neighbours each nearly level with the next, or
    // whose edges do not chain into one loop, keeps its triangles.
    let mut members: Vec<Vec<usize>> = vec![Vec::new(); hull.faces.len()];
    for t in 0..hull.faces.len() {
        members[root(&mut facet, t)].push(t);
    }
    let mut loops: Vec<(Vec<usize>, usize)> = Vec::new();
    for tris in members.iter().filter(|m| !m.is_empty()) {
        let largest = tris
            .iter()
            .copied()
            .max_by(|&s, &t| planes[s].2.total_cmp(&planes[t].2))
            .unwrap_or(tris[0]);
        let directed: HashSet<(usize, usize)> = tris
            .iter()
            .flat_map(|&t| {
                let [a, b, c] = hull.faces[t];
                [(a, b), (b, c), (c, a)]
            })
            .collect();
        let rim: HashMap<usize, usize> = directed
            .iter()
            .filter(|&&(a, b)| !directed.contains(&(b, a)))
            .copied()
            .collect();
        let first = rim.keys().copied().min();
        let mut ring = Vec::with_capacity(rim.len());
        let mut at = first;
        while let Some(v) = at {
            ring.push(v);
            at = rim.get(&v).copied().filter(|&n| Some(n) != first);
            if ring.len() > rim.len() {
                break;
            }
        }
        let level = tris
            .iter()
            .all(|&t| hull.faces[t].into_iter().all(|v| on_plane(largest, v)));
        if level && ring.len() == rim.len() && rim.len() >= 3 {
            loops.push((ring, largest));
        } else {
            loops.extend(tris.iter().map(|&t| (hull.faces[t].to_vec(), t)));
        }
    }

    // A hull point on an edge between two facets (a vertex sum landing midway
    // along a Minkowski sum's edge) leaves both facets' rims straight through
    // it; it is no corner of the solid, and both rims drop it.
    let mut facets_at: HashMap<usize, usize> = HashMap::new();
    for (ring, _) in &loops {
        for &v in ring {
            *facets_at.entry(v).or_default() += 1;
        }
    }
    let straight = |a: usize, v: usize, b: usize| {
        let (pa, pv, pb) = (hull.vertices[a], hull.vertices[v], hull.vertices[b]);
        (pa - pv).cross(pb - pv).length() <= 1e-12 * scale * scale
    };
    for (ring, _) in &mut loops {
        let mut k = 0;
        while ring.len() > 3 && k < ring.len() {
            let n = ring.len();
            let (a, v, b) = (ring[(k + n - 1) % n], ring[k], ring[(k + 1) % n]);
            if facets_at.get(&v) == Some(&2) && straight(a, v, b) {
                ring.remove(k);
            } else {
                k += 1;
            }
        }
    }

    let mut vertex_ids: HashMap<usize, brepkit_topology::vertex::VertexId> = HashMap::new();
    for (ring, _) in &loops {
        for &v in ring {
            vertex_ids
                .entry(v)
                .or_insert_with(|| topo.add_vertex(Vertex::new(hull.vertices[v], tol.linear)));
        }
    }

    // Each undirected edge (min_vertex, max_vertex) is created once; the second
    // face that references it uses reverse orientation.
    let mut edge_map: HashMap<(usize, usize), brepkit_topology::edge::EdgeId> = HashMap::new();
    let mut face_ids = Vec::with_capacity(loops.len());
    for (ring, tri) in &loops {
        let mut oriented_edges = Vec::with_capacity(ring.len());
        for (k, &idx_a) in ring.iter().enumerate() {
            let idx_b = ring[(k + 1) % ring.len()];
            let key = (idx_a.min(idx_b), idx_a.max(idx_b));
            let (eid, forward) = if let Some(&existing) = edge_map.get(&key) {
                (existing, false)
            } else {
                let eid = topo.add_edge(Edge::new(
                    vertex_ids[&idx_a],
                    vertex_ids[&idx_b],
                    EdgeCurve::Line,
                ));
                edge_map.insert(key, eid);
                (eid, true)
            };
            oriented_edges.push(OrientedEdge::new(eid, forward));
        }

        let wire = Wire::new(oriented_edges, true).map_err(crate::OperationsError::Topology)?;
        let wid = topo.add_wire(wire);
        let (normal, d, _) = planes[*tri];
        let fid = topo.add_face(Face::new(wid, vec![], FaceSurface::Plane { normal, d }));
        face_ids.push(fid);
    }

    let shell = Shell::new(face_ids).map_err(crate::OperationsError::Topology)?;
    let shell_id = topo.add_shell(shell);
    let solid = Solid::new(shell_id, vec![]);
    Ok(topo.add_solid(solid))
}

/// Convex Minkowski sum of two solids: `A ⊕ B = { a + b | a ∈ A, b ∈ B }`.
///
/// Computed as the convex hull of all pairwise vertex sums. This is exact when
/// both inputs are convex **polytopes** (e.g. boxes, or a shape with a
/// tessellated sphere as the rolling tool); for a smooth convex solid the
/// vertices under-sample the surface, and for non-convex inputs it is a convex
/// over-approximation.
///
/// The vertex sum is `O(|A|·|B|)` points — large meshes can be expensive.
///
/// # Errors
///
/// Returns [`crate::OperationsError::InvalidInput`] if either solid has no
/// vertices, and propagates [`make_convex_hull`] errors (e.g. when the summed
/// points are degenerate / coplanar and no hull can be built).
pub fn make_minkowski_sum(
    topo: &mut Topology,
    a: SolidId,
    b: SolidId,
) -> Result<SolidId, crate::OperationsError> {
    use brepkit_topology::explorer::solid_vertices;

    let point_of = |topo: &Topology, ids: &[brepkit_topology::vertex::VertexId]| {
        ids.iter()
            .map(|&v| topo.vertex(v).map(brepkit_topology::vertex::Vertex::point))
            .collect::<Result<Vec<Point3>, _>>()
    };
    let a_pts = point_of(topo, &solid_vertices(topo, a)?)?;
    let b_pts = point_of(topo, &solid_vertices(topo, b)?)?;

    if a_pts.is_empty() || b_pts.is_empty() {
        return Err(crate::OperationsError::InvalidInput {
            reason: "minkowski sum requires two non-empty solids".into(),
        });
    }

    // A genuinely overflowing vertex product is an impossibly large input —
    // return a clean error rather than panicking in `Vec::with_capacity`.
    let capacity = a_pts.len().checked_mul(b_pts.len()).ok_or_else(|| {
        crate::OperationsError::InvalidInput {
            reason: "minkowski sum input too large: vertex-count product overflows".into(),
        }
    })?;
    let mut sums = Vec::with_capacity(capacity);
    for &p in &a_pts {
        for &q in &b_pts {
            sums.push(Point3::new(p.x() + q.x(), p.y() + q.y(), p.z() + q.z()));
        }
    }
    make_convex_hull(topo, &sums)
}

#[cfg(test)]
mod tests;