BREP_kernel 0.2.0

use crate::analytic_surface::{circumcenter, AnalyticSurface};
use crate::topology::{BrepSolid, EdgeRecord, FaceRecord, VertexRecord};
use crate::{make_arc, KernelTolerances, NurbsCurve, NurbsSurface, Vec3};
use rustc_hash::FxHashMap as HashMap;

fn step_string(value: &str) -> String {
    value.replace('\'', "''")
}

fn real(value: f64) -> Result<String, String> {
    if !value.is_finite() {
        return Err(format!("export_step: non-finite number {value}"));
    }
    // Analytic frame axes come from cross products that can round to -0.0;
    // normalize so directions never print a negative zero component.
    if value == 0.0 {
        return Ok("0.".into());
    }
    if value.fract() == 0.0 && value.abs() < 1e15 {
        return Ok(format!("{value:.0}."));
    }
    let mut output = format!("{value:.15}");
    while output.ends_with('0') {
        output.pop();
    }
    if output.ends_with('.') {
        output.push('0');
    }
    if output == "-0.0" {
        output = "0.0".into();
    }
    Ok(output)
}

fn knot_runs(knots: &[f64]) -> (Vec<f64>, Vec<usize>) {
    let mut values = Vec::new();
    let mut multiplicities = Vec::new();
    for &knot in knots {
        if values
            .last()
            .is_some_and(|previous: &f64| (*previous - knot).abs() <= 1e-12)
        {
            *multiplicities.last_mut().unwrap() += 1;
        } else {
            values.push(knot);
            multiplicities.push(1);
        }
    }
    (values, multiplicities)
}

fn edge_subcurve(edge: &EdgeRecord) -> Result<NurbsCurve, String> {
    let [start, end] = edge.curve.domain()?;
    let epsilon = (1e-9 * (end - start)).max(2e-9);
    let mut curve = edge.curve.clone();
    if edge.t0 > start + epsilon && edge.t0 < end - epsilon {
        curve = curve.split(edge.t0)?.1;
    }
    let domain = curve.domain()?;
    if edge.t1 < domain[1] - epsilon && edge.t1 > domain[0] + epsilon {
        curve = curve.split(edge.t1)?.0;
    }
    Ok(curve)
}

#[derive(Default)]
struct StepWriter {
    lines: Vec<String>,
}

impl StepWriter {
    fn add(&mut self, body: impl Into<String>) -> usize {
        let id = self.lines.len() + 1;
        self.lines.push(format!("#{id}={};", body.into()));
        id
    }

    fn data(&self) -> String {
        self.lines.join("\n")
    }
}

fn write_point(writer: &mut StepWriter, point: Vec3) -> Result<usize, String> {
    Ok(writer.add(format!(
        "CARTESIAN_POINT('',({},{},{}))",
        real(point.x)?,
        real(point.y)?,
        real(point.z)?
    )))
}

fn id_list(ids: &[usize]) -> String {
    format!(
        "({})",
        ids.iter()
            .map(|id| format!("#{id}"))
            .collect::<Vec<_>>()
            .join(",")
    )
}

fn write_direction(writer: &mut StepWriter, direction: Vec3) -> Result<usize, String> {
    Ok(writer.add(format!(
        "DIRECTION('',({},{},{}))",
        real(direction.x)?,
        real(direction.y)?,
        real(direction.z)?
    )))
}

fn write_placement(
    writer: &mut StepWriter,
    origin: Vec3,
    axis: Vec3,
    ref_direction: Vec3,
) -> Result<usize, String> {
    let origin = write_point(writer, origin)?;
    let axis = write_direction(writer, axis)?;
    let ref_direction = write_direction(writer, ref_direction)?;
    Ok(writer.add(format!(
        "AXIS2_PLACEMENT_3D('',#{origin},#{axis},#{ref_direction})"
    )))
}

/// Emit the analytic AP214 surface entity (PLANE / CYLINDRICAL_SURFACE /
/// CONICAL_SURFACE / SPHERICAL_SURFACE / TOROIDAL_SURFACE) for a recognized
/// carrier, or `None` when the surface must stay a B-spline. The second value
/// reports whether the STEP-standard orientation of the emitted entity (plane
/// normal along the placement axis; revolution normal outward) is the REVERSE
/// of the stored NURBS orientation, so ADVANCED_FACE can invert `same_sense`
/// and the face normal survives the round trip.
fn write_analytic_surface(
    writer: &mut StepWriter,
    surface: &NurbsSurface,
) -> Result<Option<(usize, bool)>, String> {
    let Some(analytic) = surface.analytic() else {
        return Ok(None);
    };
    match analytic {
        AnalyticSurface::Plane {
            origin,
            u_dir,
            v_dir,
            ..
        } => {
            // STEP planes are unbounded; the importer re-sizes the patch from
            // the face's edges, so only origin/normal/ref matter.
            let (Ok(normal), Ok(x_axis)) = (u_dir.cross(*v_dir).normalized(), u_dir.normalized())
            else {
                return Ok(None);
            };
            let placement = write_placement(writer, *origin, normal, x_axis)?;
            Ok(Some((writer.add(format!("PLANE('',#{placement})")), false)))
        }
        AnalyticSurface::RuledRevolution {
            frame,
            rho0,
            rho1,
            height,
        } => {
            // The recognizer allows height < 0 (descending generatrix), whose
            // normal is the reverse of the standard outward convention the
            // importer reconstructs; report that so the face sense compensates.
            let flipped = *height < 0.0;
            let radius_scale = 1.0 + rho0.abs().max(rho1.abs());
            if (rho1 - rho0).abs() <= 1e-9 * radius_scale {
                if *rho0 <= 0.0 {
                    return Ok(None);
                }
                let placement = write_placement(writer, frame.origin, frame.axis, frame.x_axis)?;
                return Ok(Some((
                    writer.add(format!(
                        "CYLINDRICAL_SURFACE('',#{placement},{})",
                        real(*rho0)?
                    )),
                    flipped,
                )));
            }
            // Cone. The importer re-sizes the carrier from edge samples with a
            // 1e-4-scaled axial margin and clamps a negative extended-end
            // radius to zero — which BENDS the rebuilt slope when the apex sits
            // at an end of the face's axial range. Keep apex-touching cones as
            // exact NURBS instead of exporting a distorted carrier.
            let slope = (rho1 - rho0) / height;
            let apex_margin = 2.0 * (1e-4 * height.abs().max(1.0) + 1e-9) * slope.abs();
            if rho0.min(*rho1) <= apex_margin {
                return Ok(None);
            }
            // Orient the placement axis so the radius grows along +axis: STEP
            // semi-angles are positive. Radius at the placement origin stays
            // rho0 either way because the origin is on-axis at the v = 0 base.
            let axis = if slope >= 0.0 {
                frame.axis
            } else {
                frame.axis.scale(-1.0)
            };
            let placement = write_placement(writer, frame.origin, axis, frame.x_axis)?;
            Ok(Some((
                writer.add(format!(
                    "CONICAL_SURFACE('',#{placement},{},{})",
                    real(*rho0)?,
                    real(slope.abs().atan())?
                )),
                flipped,
            )))
        }
        AnalyticSurface::Sphere { frame, radius } => {
            // Recognition template and importer reconstruction share the same
            // south-to-north meridian construction, so the rebuild is exact.
            let placement = write_placement(writer, frame.origin, frame.axis, frame.x_axis)?;
            Ok(Some((
                writer.add(format!(
                    "SPHERICAL_SURFACE('',#{placement},{})",
                    real(*radius)?
                )),
                false,
            )))
        }
        AnalyticSurface::Torus {
            frame,
            major_radius,
            minor_radius,
        } => {
            let placement = write_placement(writer, frame.origin, frame.axis, frame.x_axis)?;
            Ok(Some((
                writer.add(format!(
                    "TOROIDAL_SURFACE('',#{placement},{},{})",
                    real(*major_radius)?,
                    real(*minor_radius)?
                )),
                false,
            )))
        }
        // No SURFACE_OF_REVOLUTION reader exists yet; general revolutions keep
        // their exact NURBS form.
        AnalyticSurface::Revolution { .. } => Ok(None),
    }
}

/// A curve recognized as an exact `make_arc` product: a circular arc of
/// `radius` about `axis`, starting at angle 0 on `x_axis` and travelling
/// counterclockwise through `sweep` — exactly the CIRCLE parameterization the
/// importer trims between the edge's vertices.
struct CircularArc {
    center: Vec3,
    axis: Vec3,
    x_axis: Vec3,
    radius: f64,
}

fn curve_scale(curve: &NurbsCurve) -> f64 {
    curve
        .control_points
        .iter()
        .map(|p| (p.x.abs() / p.w).max(p.y.abs() / p.w).max(p.z.abs() / p.w))
        .fold(0.0, f64::max)
}

/// Homogeneous-net equality to a scale-relative tolerance (the same
/// reconstruction contract analytic_surface.rs uses): matching nets mean the
/// curves are the SAME exact rational arc, not merely close.
fn curves_match(a: &NurbsCurve, b: &NurbsCurve, scale: f64) -> bool {
    if a.degree != b.degree
        || a.knots.len() != b.knots.len()
        || a.control_points.len() != b.control_points.len()
    {
        return false;
    }
    if a.knots
        .iter()
        .zip(&b.knots)
        .any(|(x, y)| (x - y).abs() > 1e-12)
    {
        return false;
    }
    let tolerance = 1e-9 * scale.max(1.0);
    a.control_points
        .iter()
        .zip(&b.control_points)
        .all(|(p, q)| {
            (p.x - q.x).abs() <= tolerance
                && (p.y - q.y).abs() <= tolerance
                && (p.z - q.z).abs() <= tolerance
                && (p.w - q.w).abs() <= 1e-9
        })
}

/// Recognition by exact reconstruction: extract a candidate circle from three
/// curve points, rebuild it with `make_arc`, and demand the identical net.
/// Split subranges of a circle (whose knots are no longer the pristine
/// make_arc pattern) are rejected and honestly stay NURBS.
fn recognize_circular_arc(curve: &NurbsCurve) -> Option<CircularArc> {
    if curve.degree != 2
        || curve.control_points.len() < 3
        || curve.control_points.len() % 2 == 0
        || (curve.control_points.len() - 1) / 2 > 4
    {
        return None;
    }
    let [t0, t1] = curve.domain().ok()?;
    let at = |fraction: f64| curve.evaluate(t0 + (t1 - t0) * fraction);
    // Three points at < 74% of the sweep apart, so consecutive pairs subtend
    // less than pi and the cross product below gives the travel direction.
    let p0 = at(0.0).ok()?;
    let pa = at(0.35).ok()?;
    let pb = at(0.7).ok()?;
    let center = circumcenter(p0, pa, pb)?;
    let radial = p0.sub(center);
    let radius = radial.length();
    let scale = curve_scale(curve);
    if radius <= 1e-9 * scale.max(1.0) {
        return None;
    }
    let x_axis = radial.scale(1.0 / radius);
    let axis = radial.cross(pa.sub(center)).normalized().ok()?;
    let y_axis = axis.cross(x_axis);
    let p_end = at(1.0).ok()?;
    let sweep = if p_end.sub(p0).length() <= 1e-9 * (1.0 + radius) {
        std::f64::consts::TAU
    } else {
        let closing = p_end.sub(center);
        let mut angle = closing.dot(y_axis).atan2(closing.dot(x_axis));
        if angle < 0.0 {
            angle += std::f64::consts::TAU;
        }
        angle
    };
    let rebuilt = make_arc(center, x_axis, y_axis, radius, 0.0, sweep).ok()?;
    curves_match(curve, &rebuilt, scale).then_some(CircularArc {
        center,
        axis,
        x_axis,
        radius,
    })
}

/// Emit LINE or CIRCLE for a recognized analytic edge curve (already oriented
/// start-to-end by `edge_subcurve`), or `None` for the B-spline fallback.
fn write_analytic_curve(
    writer: &mut StepWriter,
    curve: &NurbsCurve,
) -> Result<Option<usize>, String> {
    if curve.degree == 1
        && curve.control_points.len() == 2
        && curve
            .control_points
            .iter()
            .all(|control| (control.w - 1.0).abs() <= 1e-12)
    {
        let start = curve.control_points[0].point()?;
        let end = curve.control_points[1].point()?;
        let Ok(direction) = end.sub(start).normalized() else {
            return Ok(None);
        };
        let point = write_point(writer, start)?;
        let step_direction = write_direction(writer, direction)?;
        let vector = writer.add(format!(
            "VECTOR('',#{step_direction},{})",
            real(end.sub(start).length())?
        ));
        return Ok(Some(writer.add(format!("LINE('',#{point},#{vector})"))));
    }
    if let Some(arc) = recognize_circular_arc(curve) {
        // ref_direction points at the edge's start vertex and the arc runs
        // counterclockwise about the axis, so the importer's vertex-trimmed
        // CCW rebuild reproduces the same directed curve with sense .T.
        let placement = write_placement(writer, arc.center, arc.axis, arc.x_axis)?;
        return Ok(Some(
            writer.add(format!("CIRCLE('',#{placement},{})", real(arc.radius)?)),
        ));
    }
    Ok(None)
}

fn write_curve(writer: &mut StepWriter, curve: &NurbsCurve) -> Result<usize, String> {
    let points = curve
        .control_points
        .iter()
        .map(|control| write_point(writer, control.point()?))
        .collect::<Result<Vec<_>, _>>()?;
    let (knot_values, multiplicities) = knot_runs(&curve.knots);
    let multiplicities = format!(
        "({})",
        multiplicities
            .iter()
            .map(usize::to_string)
            .collect::<Vec<_>>()
            .join(",")
    );
    let knots = format!(
        "({})",
        knot_values
            .iter()
            .map(|value| real(*value))
            .collect::<Result<Vec<_>, _>>()?
            .join(",")
    );
    let rational = curve
        .control_points
        .iter()
        .any(|control| (control.w - 1.0).abs() > 1e-12);
    if !rational {
        return Ok(writer.add(format!(
            "B_SPLINE_CURVE_WITH_KNOTS('',{},{},.UNSPECIFIED.,.F.,.F.,{multiplicities},{knots},.UNSPECIFIED.)",
            curve.degree,
            id_list(&points),
        )));
    }
    let weights = format!(
        "({})",
        curve
            .control_points
            .iter()
            .map(|control| real(control.w))
            .collect::<Result<Vec<_>, _>>()?
            .join(",")
    );
    Ok(writer.add(format!(
        "(BOUNDED_CURVE()B_SPLINE_CURVE({},{},.UNSPECIFIED.,.F.,.F.)\
         B_SPLINE_CURVE_WITH_KNOTS({multiplicities},{knots},.UNSPECIFIED.)\
         CURVE()GEOMETRIC_REPRESENTATION_ITEM()RATIONAL_B_SPLINE_CURVE({weights})\
         REPRESENTATION_ITEM(''))",
        curve.degree,
        id_list(&points),
    )))
}

fn write_surface(writer: &mut StepWriter, surface: &NurbsSurface) -> Result<usize, String> {
    let rows = surface
        .control_points
        .iter()
        .map(|row| {
            row.iter()
                .map(|control| write_point(writer, control.point()?))
                .collect::<Result<Vec<_>, _>>()
                .map(|ids| id_list(&ids))
        })
        .collect::<Result<Vec<_>, _>>()?;
    let grid = format!("({})", rows.join(","));
    let (u_values, u_multiplicities) = knot_runs(&surface.knots_u);
    let (v_values, v_multiplicities) = knot_runs(&surface.knots_v);
    let multiplicities = |values: &[usize]| {
        format!(
            "({})",
            values
                .iter()
                .map(usize::to_string)
                .collect::<Vec<_>>()
                .join(",")
        )
    };
    let knots = |values: &[f64]| -> Result<String, String> {
        Ok(format!(
            "({})",
            values
                .iter()
                .map(|value| real(*value))
                .collect::<Result<Vec<_>, _>>()?
                .join(",")
        ))
    };
    let u_mults = multiplicities(&u_multiplicities);
    let v_mults = multiplicities(&v_multiplicities);
    let u_knots = knots(&u_values)?;
    let v_knots = knots(&v_values)?;
    let rational = surface
        .control_points
        .iter()
        .flatten()
        .any(|control| (control.w - 1.0).abs() > 1e-12);
    if !rational {
        return Ok(writer.add(format!(
            "B_SPLINE_SURFACE_WITH_KNOTS('',{},{},{grid},.UNSPECIFIED.,.F.,.F.,.F.,\
             {u_mults},{v_mults},{u_knots},{v_knots},.UNSPECIFIED.)",
            surface.degree_u, surface.degree_v,
        )));
    }
    let weights = format!(
        "({})",
        surface
            .control_points
            .iter()
            .map(|row| {
                row.iter()
                    .map(|control| real(control.w))
                    .collect::<Result<Vec<_>, _>>()
                    .map(|values| format!("({})", values.join(",")))
            })
            .collect::<Result<Vec<_>, _>>()?
            .join(",")
    );
    Ok(writer.add(format!(
        "(BOUNDED_SURFACE()B_SPLINE_SURFACE({},{},{grid},.UNSPECIFIED.,.F.,.F.,.F.)\
         B_SPLINE_SURFACE_WITH_KNOTS({u_mults},{v_mults},{u_knots},{v_knots},.UNSPECIFIED.)\
         GEOMETRIC_REPRESENTATION_ITEM()RATIONAL_B_SPLINE_SURFACE({weights})\
         REPRESENTATION_ITEM('')SURFACE())",
        surface.degree_u, surface.degree_v,
    )))
}

fn write_length_unit(writer: &mut StepWriter, unit: &str) -> Result<usize, String> {
    let normalized = unit.to_lowercase();
    if normalized == "meter" || normalized == "metre" {
        return Ok(writer.add("(LENGTH_UNIT()NAMED_UNIT(*)SI_UNIT($,.METRE.))"));
    }
    if normalized == "centimeter" || normalized == "centimetre" {
        return Ok(writer.add("(LENGTH_UNIT()NAMED_UNIT(*)SI_UNIT(.CENTI.,.METRE.))"));
    }
    if matches!(normalized.as_str(), "micron" | "micrometer" | "micrometre") {
        return Ok(writer.add("(LENGTH_UNIT()NAMED_UNIT(*)SI_UNIT(.MICRO.,.METRE.))"));
    }
    if normalized == "inch" || normalized == "foot" {
        let metre = writer.add("(LENGTH_UNIT()NAMED_UNIT(*)SI_UNIT($,.METRE.))");
        let (factor, name) = if normalized == "inch" {
            (0.0254, "INCH")
        } else {
            (0.3048, "FOOT")
        };
        let measure = writer.add(format!(
            "LENGTH_MEASURE_WITH_UNIT(LENGTH_MEASURE({}),#{metre})",
            real(factor)?
        ));
        return Ok(writer.add(format!(
            "(CONVERSION_BASED_UNIT('{name}',#{measure})LENGTH_UNIT()NAMED_UNIT(*))"
        )));
    }
    Ok(writer.add("(LENGTH_UNIT()NAMED_UNIT(*)SI_UNIT(.MILLI.,.METRE.))"))
}

fn vertex_for(solid: &BrepSolid, id: u64) -> Result<&VertexRecord, String> {
    solid
        .vertices
        .iter()
        .find(|vertex| vertex.id == id)
        .ok_or_else(|| format!("export_step: missing vertex {id}"))
}

fn edge_for(solid: &BrepSolid, id: u64) -> Result<&EdgeRecord, String> {
    solid
        .edges
        .iter()
        .find(|edge| edge.id == id)
        .ok_or_else(|| format!("export_step: missing edge {id}"))
}

fn surface_key(face: &FaceRecord) -> usize {
    face as *const FaceRecord as usize
}

/// Serialize exact NURBS BREP topology as an AP214 STEP Part 21 document.
pub fn export_step(
    solids: &[BrepSolid],
    name: &str,
    unit: &str,
    timestamp: &str,
) -> Result<String, String> {
    if solids.is_empty() {
        return Err("export_step: at least one solid is required".into());
    }
    for solid in solids {
        let policy = KernelTolerances::for_solid(solid, 1e-7);
        let issues = solid.validate_with_tolerances(&KernelTolerances {
            pcurve_consistency: policy.export_knit,
            ..policy
        });
        if !issues.is_empty() {
            return Err(format!("export_step: invalid solid: {issues:?}"));
        }
    }
    let mut writer = StepWriter::default();
    let safe_name = step_string(name);
    let application = writer.add("APPLICATION_CONTEXT('automotive design')");
    writer.add(format!(
        "APPLICATION_PROTOCOL_DEFINITION('','automotive_design',2010,#{application})"
    ));
    let product_context = writer.add(format!("PRODUCT_CONTEXT('',#{application},'mechanical')"));
    let product = writer.add(format!(
        "PRODUCT('{safe_name}','{safe_name}','',(#{product_context}))"
    ));
    let formation = writer.add(format!("PRODUCT_DEFINITION_FORMATION('','',#{product})"));
    let definition_context = writer.add(format!(
        "PRODUCT_DEFINITION_CONTEXT('part definition',#{application},'design')"
    ));
    let definition = writer.add(format!(
        "PRODUCT_DEFINITION('design','',#{formation},#{definition_context})"
    ));
    let product_shape = writer.add(format!("PRODUCT_DEFINITION_SHAPE('','',#{definition})"));
    let length_unit = write_length_unit(&mut writer, unit)?;
    let angle_unit = writer.add("(NAMED_UNIT(*)PLANE_ANGLE_UNIT()SI_UNIT($,.RADIAN.))");
    let solid_angle_unit = writer.add("(NAMED_UNIT(*)SI_UNIT($,.STERADIAN.)SOLID_ANGLE_UNIT())");
    let uncertainty = writer.add(format!(
        "UNCERTAINTY_MEASURE_WITH_UNIT(LENGTH_MEASURE(1.E-6),#{length_unit},'distance_accuracy_value','')"
    ));
    let geometry_context = writer.add(format!(
        "(GEOMETRIC_REPRESENTATION_CONTEXT(3)\
         GLOBAL_UNCERTAINTY_ASSIGNED_CONTEXT((#{uncertainty}))\
         GLOBAL_UNIT_ASSIGNED_CONTEXT((#{length_unit},#{angle_unit},#{solid_angle_unit}))\
         REPRESENTATION_CONTEXT('',''))"
    ));
    let origin = write_point(&mut writer, Vec3::default())?;
    let direction_z = writer.add("DIRECTION('',(0.,0.,1.))");
    let direction_x = writer.add("DIRECTION('',(1.,0.,0.))");
    let axis = writer.add(format!(
        "AXIS2_PLACEMENT_3D('',#{origin},#{direction_z},#{direction_x})"
    ));

    let mut solid_ids = Vec::new();
    for solid in solids {
        let mut vertex_ids = HashMap::<u64, usize>::default();
        let mut edge_ids = HashMap::<u64, usize>::default();
        let mut surface_ids = HashMap::<usize, (usize, bool)>::default();
        for shell in &solid.shells {
            let mut face_ids = Vec::new();
            for face in &shell.faces {
                let mut bound_ids = Vec::new();
                for (loop_index, loop_record) in face.loops.iter().enumerate() {
                    let mut oriented_edges = Vec::new();
                    for coedge in &loop_record.coedges {
                        let edge = edge_for(solid, coedge.edge_id)?;
                        if edge.degenerate {
                            continue;
                        }
                        let edge_id = if let Some(id) = edge_ids.get(&edge.id) {
                            *id
                        } else {
                            let subcurve = edge_subcurve(edge)?;
                            let curve = match write_analytic_curve(&mut writer, &subcurve)? {
                                Some(id) => id,
                                None => write_curve(&mut writer, &subcurve)?,
                            };
                            let mut vertex_id =
                                |id: u64, writer: &mut StepWriter| -> Result<usize, String> {
                                    if let Some(step_id) = vertex_ids.get(&id) {
                                        return Ok(*step_id);
                                    }
                                    let point = write_point(writer, vertex_for(solid, id)?.point)?;
                                    let step_id = writer.add(format!("VERTEX_POINT('',#{point})"));
                                    vertex_ids.insert(id, step_id);
                                    Ok(step_id)
                                };
                            let start = vertex_id(edge.start_vertex_id, &mut writer)?;
                            let end = vertex_id(edge.end_vertex_id, &mut writer)?;
                            let step_id =
                                writer.add(format!("EDGE_CURVE('',#{start},#{end},#{curve},.T.)"));
                            edge_ids.insert(edge.id, step_id);
                            step_id
                        };
                        let orientation = if coedge.forward { ".T." } else { ".F." };
                        oriented_edges.push(
                            writer.add(format!("ORIENTED_EDGE('',*,*,#{edge_id},{orientation})")),
                        );
                    }
                    if oriented_edges.is_empty() {
                        continue;
                    }
                    let edge_loop =
                        writer.add(format!("EDGE_LOOP('',{})", id_list(&oriented_edges)));
                    let kind = if loop_index == 0 {
                        "FACE_OUTER_BOUND"
                    } else {
                        "FACE_BOUND"
                    };
                    bound_ids.push(writer.add(format!("{kind}('',#{edge_loop},.T.)")));
                }
                let key = surface_key(face);
                let (surface, flipped) = if let Some(entry) = surface_ids.get(&key) {
                    *entry
                } else {
                    let entry = match write_analytic_surface(&mut writer, &face.surface)? {
                        Some(pair) => pair,
                        None => (write_surface(&mut writer, &face.surface)?, false),
                    };
                    surface_ids.insert(key, entry);
                    entry
                };
                // `flipped` analytic entities are written with the reverse of
                // the stored NURBS orientation, so invert the flag to keep the
                // face normal identical through the round trip.
                let sense = if face.same_sense != flipped {
                    ".T."
                } else {
                    ".F."
                };
                face_ids.push(writer.add(format!(
                    "ADVANCED_FACE('',{},#{surface},{sense})",
                    id_list(&bound_ids)
                )));
            }
            let closed_shell = writer.add(format!("CLOSED_SHELL('',{})", id_list(&face_ids)));
            solid_ids.push(writer.add(format!(
                "MANIFOLD_SOLID_BREP('{safe_name}',#{closed_shell})"
            )));
        }
    }
    let mut items = vec![axis];
    items.extend(&solid_ids);
    let representation = writer.add(format!(
        "ADVANCED_BREP_SHAPE_REPRESENTATION('',{},#{geometry_context})",
        id_list(&items)
    ));
    writer.add(format!(
        "SHAPE_DEFINITION_REPRESENTATION(#{product_shape},#{representation})"
    ));
    let safe_timestamp = step_string(timestamp);
    let output = [
        "ISO-10303-21;".to_string(),
        "HEADER;".to_string(),
        "FILE_DESCRIPTION((''),'2;1');".to_string(),
        format!(
            "FILE_NAME('{safe_name}.step','{safe_timestamp}',(''),(''),'brep-kernel-rs','brep-kernel-rs','');"
        ),
        "FILE_SCHEMA(('AUTOMOTIVE_DESIGN { 1 0 10303 214 1 1 1 1 }'));".to_string(),
        "ENDSEC;".to_string(),
        "DATA;".to_string(),
        writer.data(),
        "ENDSEC;".to_string(),
        "END-ISO-10303-21;".to_string(),
        String::new(),
    ]
    .join("\n");
    let manifold_issues = audit_step_manifold(&output);
    if !manifold_issues.is_empty() {
        return Err(format!(
            "export_step: emitted AP214 manifold audit failed: {}",
            manifold_issues.join("; ")
        ));
    }
    Ok(output)
}

/// Audit the serialized entity graph rather than assuming that valid
/// in-memory topology was necessarily written correctly.  Every EDGE_CURVE
/// in a closed shell must have exactly two ORIENTED_EDGE users with opposite
/// senses.
pub fn audit_step_manifold(step: &str) -> Vec<String> {
    let marker = "ORIENTED_EDGE('',*,*,#";
    let mut uses = HashMap::<u64, Vec<bool>>::default();
    for line in step.lines() {
        let Some(offset) = line.find(marker) else {
            continue;
        };
        let rest = &line[offset + marker.len()..];
        let digits = rest
            .chars()
            .take_while(|character| character.is_ascii_digit())
            .collect::<String>();
        let Ok(edge_id) = digits.parse::<u64>() else {
            continue;
        };
        let suffix = &rest[digits.len()..];
        let sense = suffix.starts_with(",.T.");
        uses.entry(edge_id).or_default().push(sense);
    }
    let mut issues = uses
        .into_iter()
        .filter_map(|(edge, senses)| {
            (senses.len() != 2 || senses[0] == senses[1]).then(|| {
                format!(
                    "EDGE_CURVE #{edge} has {} uses with senses {:?}",
                    senses.len(),
                    senses
                )
            })
        })
        .collect::<Vec<_>>();
    issues.sort();
    issues
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::{
        boolean_operation, import_step, make_box_brep, make_cone_brep, make_cylinder_brep,
        make_cylinder_surface, make_sphere_brep, make_torus_brep, solid_mass_properties,
        BooleanOperation, BooleanOptions,
    };

    #[test]
    fn box_step_contains_exact_manifold_topology() {
        let box_solid = make_box_brep(Vec3::default(), 2.0, 3.0, 4.0).unwrap();
        let step = export_step(&[box_solid], "box", "millimeter", "2026-07-27T00:00:00").unwrap();
        assert!(step.starts_with("ISO-10303-21;\nHEADER;"));
        assert!(audit_step_manifold(&step).is_empty());
        assert!(step.contains("MANIFOLD_SOLID_BREP('box'"));
        assert_eq!(step.matches("ADVANCED_FACE(").count(), 6);
        assert_eq!(step.matches("EDGE_CURVE(").count(), 12);
        assert!(step.ends_with("END-ISO-10303-21;\n"));
    }

    #[test]
    fn step_manifold_audit_rejects_single_and_same_sense_uses() {
        let single = "#1=ORIENTED_EDGE('',*,*,#9,.T.);";
        assert_eq!(audit_step_manifold(single).len(), 1);
        let same = "#1=ORIENTED_EDGE('',*,*,#9,.T.);\n\
                    #2=ORIENTED_EDGE('',*,*,#9,.T.);";
        assert_eq!(audit_step_manifold(same).len(), 1);
        let good = "#1=ORIENTED_EDGE('',*,*,#9,.T.);\n\
                    #2=ORIENTED_EDGE('',*,*,#9,.F.);";
        assert!(audit_step_manifold(good).is_empty());
    }

    #[test]
    fn unrecognized_surfaces_and_curves_still_write_rational_complex_entities() {
        // The all-NURBS fallback writers must stay intact for carriers no
        // analytic entity covers (general revolutions, split arc subranges).
        let mut writer = StepWriter::default();
        let cylinder =
            make_cylinder_surface(Vec3::default(), Vec3::new(0.0, 0.0, 1.0), 2.0, 5.0).unwrap();
        write_surface(&mut writer, &cylinder).unwrap();
        let split_arc = make_arc(
            Vec3::default(),
            Vec3::new(1.0, 0.0, 0.0),
            Vec3::new(0.0, 1.0, 0.0),
            2.0,
            0.0,
            std::f64::consts::TAU,
        )
        .unwrap()
        .split(0.37)
        .unwrap()
        .1;
        assert!(
            recognize_circular_arc(&split_arc).is_none(),
            "a split subrange is not the pristine make_arc net"
        );
        assert!(write_analytic_curve(&mut writer, &split_arc)
            .unwrap()
            .is_none());
        write_curve(&mut writer, &split_arc).unwrap();
        let data = writer.data();
        assert!(data.contains("RATIONAL_B_SPLINE_SURFACE"));
        assert!(data.contains("RATIONAL_B_SPLINE_CURVE"));
    }

    /// Export → assert the analytic entities appear (and, when the solid is
    /// fully analytic, that NO B-spline entity remains) → import → validate,
    /// match volume to 1e-6, and re-recognize every face's carrier.
    fn assert_analytic_round_trip(
        label: &str,
        original: &BrepSolid,
        expected_markers: &[&str],
        forbid_nurbs: bool,
    ) -> BrepSolid {
        let step = export_step(std::slice::from_ref(original), label, "millimeter", "fixed")
            .expect("export");
        for marker in expected_markers {
            assert!(step.contains(marker), "{label}: missing {marker}");
        }
        if forbid_nurbs {
            assert!(
                !step.contains("B_SPLINE"),
                "{label}: expected a fully analytic export"
            );
        }
        assert!(audit_step_manifold(&step).is_empty(), "{label}: audit");
        let imported = import_step(&step).expect("import");
        assert_eq!(imported.len(), 1, "{label}: one solid");
        let solid = imported.into_iter().next().unwrap();
        assert!(
            solid.validate().is_empty(),
            "{label}: imported solid invalid: {:?}",
            solid.validate()
        );
        let original_volume = solid_mass_properties(original).unwrap().volume;
        let volume = solid_mass_properties(&solid).unwrap().volume;
        let relative = ((volume - original_volume) / original_volume).abs();
        assert!(
            relative < 1e-6,
            "{label}: volume {volume} vs {original_volume} (rel {relative:.3e})"
        );
        for shell in &solid.shells {
            for face in &shell.faces {
                assert!(
                    face.surface.analytic().is_some(),
                    "{label}: imported face {} did not re-recognize as analytic",
                    face.id
                );
            }
        }
        solid
    }

    #[test]
    fn box_round_trips_through_plane_and_line_entities() {
        let solid = make_box_brep(Vec3::new(-1.0, 0.5, 2.0), 2.0, 3.0, 4.0).unwrap();
        let step = export_step(std::slice::from_ref(&solid), "box", "millimeter", "fixed").unwrap();
        assert_eq!(step.matches("PLANE(").count(), 6);
        assert_eq!(step.matches("LINE(").count(), 12);
        assert_analytic_round_trip("box", &solid, &["PLANE(", "LINE("], true);
    }

    #[test]
    fn cylinder_round_trips_through_analytic_entities() {
        let solid = make_cylinder_brep(
            Vec3::new(1.0, -2.0, 0.5),
            Vec3::new(0.0, 0.0, 1.0),
            2.0,
            5.0,
        )
        .unwrap();
        assert_analytic_round_trip(
            "cylinder",
            &solid,
            &[
                "CYLINDRICAL_SURFACE(",
                "PLANE(",
                "CIRCLE(",
                "LINE(",
                "VECTOR(",
            ],
            true,
        );
    }

    #[test]
    fn cylinder_export_keeps_unit_conversion_entities() {
        let cylinder =
            make_cylinder_brep(Vec3::default(), Vec3::new(0.0, 0.0, 1.0), 2.0, 5.0).unwrap();
        let step = export_step(&[cylinder], "cylinder", "inch", "fixed").unwrap();
        assert!(step.contains("CYLINDRICAL_SURFACE("));
        assert!(step.contains("CONVERSION_BASED_UNIT('INCH'"));
    }

    #[test]
    fn frustum_round_trips_through_conical_surface() {
        let solid = make_cone_brep(
            Vec3::new(0.5, 0.5, -1.0),
            Vec3::new(0.0, 0.0, 1.0),
            3.0,
            1.5,
            5.0,
        )
        .unwrap();
        assert_analytic_round_trip("frustum", &solid, &["CONICAL_SURFACE(", "PLANE("], true);
    }

    #[test]
    fn pointed_cone_wall_stays_nurbs_but_caps_and_rim_export_analytic() {
        // The importer's axial cover margin clamps a negative apex-end radius
        // to zero, bending the rebuilt slope; an apex-touching CONICAL export
        // would round-trip with ~1e-4 relative volume error, so the wall must
        // honestly stay NURBS while the cap plane and rim circle go analytic.
        let solid =
            make_cone_brep(Vec3::default(), Vec3::new(0.0, 0.0, 1.0), 3.0, 0.0, 5.0).unwrap();
        let step =
            export_step(std::slice::from_ref(&solid), "cone", "millimeter", "fixed").unwrap();
        assert!(!step.contains("CONICAL_SURFACE("));
        assert!(step.contains("B_SPLINE_SURFACE"));
        assert!(step.contains("PLANE("));
        assert!(step.contains("CIRCLE("));
        let imported = import_step(&step).expect("import");
        let volume = solid_mass_properties(&imported[0]).unwrap().volume;
        let expected = solid_mass_properties(&solid).unwrap().volume;
        assert!(((volume - expected) / expected).abs() < 1e-6);
    }

    #[test]
    fn sphere_round_trips_through_spherical_surface() {
        let solid =
            make_sphere_brep(Vec3::new(2.0, 1.0, -1.0), 3.0, Vec3::new(0.0, 0.0, 1.0)).unwrap();
        assert_analytic_round_trip("sphere", &solid, &["SPHERICAL_SURFACE(", "CIRCLE("], true);
    }

    #[test]
    fn torus_round_trips_through_toroidal_surface() {
        let solid =
            make_torus_brep(Vec3::new(0.0, 0.0, 1.0), Vec3::new(0.0, 0.0, 1.0), 5.0, 1.5).unwrap();
        assert_analytic_round_trip("torus", &solid, &["TOROIDAL_SURFACE(", "CIRCLE("], true);
    }

    #[test]
    fn box_minus_cylinder_round_trips_with_analytic_entities() {
        let block = make_box_brep(Vec3::new(-3.0, -3.0, 0.0), 6.0, 6.0, 4.0).unwrap();
        let drill = make_cylinder_brep(
            Vec3::new(0.0, 0.0, -1.0),
            Vec3::new(0.0, 0.0, 1.0),
            1.5,
            6.0,
        )
        .unwrap();
        let cut = boolean_operation(
            &block,
            &drill,
            BooleanOperation::Subtract,
            &BooleanOptions::default(),
        )
        .unwrap();
        // Boolean-produced edges may be split subranges (NURBS fallback), so
        // only the surface entities are required to be analytic here.
        let solid = assert_analytic_round_trip(
            "box_minus_cyl",
            &cut,
            &["CYLINDRICAL_SURFACE(", "PLANE("],
            false,
        );
        assert_eq!(solid.genus, 1, "through-hole genus survives the round trip");
    }
}