brepkit-blend 3.2.15

//! Shared utilities for fillet and chamfer builders.
//!
//! Functions used by both [`FilletBuilder`](crate::fillet_builder::FilletBuilder)
//! and [`ChamferBuilder`](crate::chamfer_builder::ChamferBuilder) for creating
//! blend faces and sampling contact curves.

use brepkit_math::nurbs::curve::NurbsCurve;
use brepkit_math::traits::ParametricSurface;
use brepkit_math::vec::{Point3, Vec3};
use brepkit_topology::Topology;
use brepkit_topology::edge::{Edge, EdgeCurve};
use brepkit_topology::face::{Face, FaceId, FaceSurface};
use brepkit_topology::vertex::{Vertex, VertexId};
use brepkit_topology::wire::{OrientedEdge, Wire};

use crate::BlendError;
use crate::stripe::Stripe;

/// Sample the start and end points of a NURBS curve.
#[must_use]
pub fn sample_nurbs_endpoints(curve: &NurbsCurve) -> Vec<Point3> {
    let (t0, t1) = curve.domain();
    vec![curve.evaluate(t0), curve.evaluate(t1)]
}

/// Create a blend face from a stripe's surface and contact curves.
///
/// Builds a minimal quadrilateral wire from the four contact-curve endpoints
/// and associates the blend surface with it.
///
/// # Errors
///
/// Returns [`BlendError`] if wire or face construction fails.
/// [`create_blend_face`] that REUSES the trimmers' contact edges when they
/// span the same contacts. Minting fresh edges for curves the trimmed
/// neighbours already carry leaves two edge entities per contact — each used
/// by one face — opening the shell along every blend flank. A trimmer edge
/// is adopted (with its vertices) when its endpoints match the stripe's
/// contact endpoints within the weld band, in either orientation; otherwise
/// that side falls back to a fresh edge.
pub fn create_blend_face_with_contacts(
    topo: &mut Topology,
    stripe: &Stripe,
    contact1_edge: Option<brepkit_topology::edge::EdgeId>,
    contact2_edge: Option<brepkit_topology::edge::EdgeId>,
) -> Result<BlendFaceInfo, BlendError> {
    const WELD: f64 = 1e-5;
    let (t0_1, t1_1) = stripe.contact1.domain();
    let (t0_2, t1_2) = stripe.contact2.domain();

    let p1_start = stripe.contact1.evaluate(t0_1);
    let p1_end = stripe.contact1.evaluate(t1_1);
    let p2_start = stripe.contact2.evaluate(t0_2);
    let p2_end = stripe.contact2.evaluate(t1_2);

    // Adopt a trimmer contact edge when its endpoints match `(want_s, want_e)`
    // in either orientation: returns (edge, forward, start_vid, end_vid) in
    // the WIRE traversal direction.
    let adopt = |topo: &Topology,
                 eid: Option<brepkit_topology::edge::EdgeId>,
                 want_s: Point3,
                 want_e: Point3|
     -> Option<(brepkit_topology::edge::EdgeId, bool, VertexId, VertexId)> {
        let eid = eid?;
        let e = topo.edge(eid).ok()?;
        let (sv, ev) = (e.start(), e.end());
        let sp = topo.vertex(sv).ok()?.point();
        let ep = topo.vertex(ev).ok()?.point();
        if (sp - want_s).length() <= WELD && (ep - want_e).length() <= WELD {
            Some((eid, true, sv, ev))
        } else if (sp - want_e).length() <= WELD && (ep - want_s).length() <= WELD {
            Some((eid, false, ev, sv))
        } else {
            None
        }
    };
    let adopt1 = adopt(topo, contact1_edge, p1_start, p1_end);
    // Contact 2 traverses end -> start in the quad below.
    let adopt2 = adopt(topo, contact2_edge, p2_end, p2_start);

    // A variable-radius stripe can pinch to a point at an end: both
    // contact curves land on the same position. Detect it up front so the
    // pinched end SHARES one vertex entity between the two contact curves
    // (the cross edge is skipped below, and separate entities would leave
    // the wire closed only positionally, not at entity level — see the
    // closure tolerance note: validation treats vertices as coincident at
    // 1e-7, tighter than the 1e-5 weld distance used here).
    let end_degenerate = (p1_end - p2_end).length() < WELD;
    let start_degenerate = (p2_start - p1_start).length() < WELD;

    // Create/reuse vertices (snapshot then allocate).
    let (v1s, v1e) = adopt1.map_or_else(
        || {
            (
                topo.add_vertex(Vertex::new(p1_start, 1e-7)),
                topo.add_vertex(Vertex::new(p1_end, 1e-7)),
            )
        },
        |(_, _, s, e)| (s, e),
    );
    let (v2e, v2s) = adopt2.map_or_else(
        || {
            (
                if end_degenerate {
                    v1e
                } else {
                    topo.add_vertex(Vertex::new(p2_end, 1e-7))
                },
                if start_degenerate {
                    v1s
                } else {
                    topo.add_vertex(Vertex::new(p2_start, 1e-7))
                },
            )
        },
        |(_, _, s, e)| (s, e),
    );

    // Build quad: p1_start -> p1_end -> p2_end -> p2_start -> p1_start.
    // Use actual contact curves for e0 and e2 (the longitudinal edges along
    // the spine direction). Cross edges e1 and e3 are straight lines connecting
    // the two contact curves at the spine endpoints.
    let (e0, e0_fwd) = adopt1.map_or_else(
        || {
            (
                topo.add_edge(Edge::new(
                    v1s,
                    v1e,
                    EdgeCurve::NurbsCurve(stripe.contact1.clone()),
                )),
                true,
            )
        },
        |(eid, fwd, _, _)| (eid, fwd),
    );
    // Cross edges carry the true end cross-section arcs when the stripe has
    // sections: the fillet's end profile is a circular arc, and a straight
    // chord both misrepresents the surface boundary and can never be shared
    // with a notched end cap. The arc's plane normal comes from the two
    // contact endpoints and the section centre.
    let arc_curve =
        |sec: &crate::section::CircSection, a: Point3, b: Point3| -> Option<EdgeCurve> {
            let u = a - sec.center;
            let v = b - sec.center;
            let n = u.cross(v);
            let n = n.normalize().ok()?;
            let circle = brepkit_math::curves::Circle3D::new(sec.center, n, sec.radius).ok()?;
            Some(EdgeCurve::Circle(circle))
        };
    let end_curve = stripe
        .sections
        .last()
        .and_then(|sec| {
            let r = arc_curve(sec, p1_end, p2_end);
            if r.is_none() {
                log::debug!(
                    "cross END line fallback: sec c={:?} r={:.5} a={p1_end:?} b={p2_end:?}",
                    sec.center,
                    sec.radius
                );
            }
            r
        })
        .unwrap_or(EdgeCurve::Line);
    let start_curve = stripe
        .sections
        .first()
        .and_then(|sec| {
            let r = arc_curve(sec, p2_start, p1_start);
            if r.is_none() {
                log::debug!(
                    "cross START line fallback: sec c={:?} r={:.5} a={p2_start:?} b={p1_start:?}",
                    sec.center,
                    sec.radius
                );
            }
            r
        })
        .unwrap_or(EdgeCurve::Line);
    // A pinched end's cross edge would be zero-length: minting it leaves a
    // degenerate use-1 edge no weld can pair; skip it.
    let e1 = if end_degenerate {
        Option::None
    } else {
        Some(topo.add_edge(Edge::new(v1e, v2e, end_curve)))
    };
    let (e2, e2_fwd) = adopt2.map_or_else(
        || {
            (
                topo.add_edge(Edge::new(
                    v2e,
                    v2s,
                    EdgeCurve::NurbsCurve(stripe.contact2.clone()),
                )),
                true,
            )
        },
        |(eid, fwd, _, _)| (eid, fwd),
    );
    let e3 = if start_degenerate {
        Option::None
    } else {
        Some(topo.add_edge(Edge::new(v2s, v1s, start_curve)))
    };

    let mut wire_edges = vec![OrientedEdge::new(e0, e0_fwd)];
    if let Some(e1) = e1 {
        wire_edges.push(OrientedEdge::new(e1, true));
    }
    wire_edges.push(OrientedEdge::new(e2, e2_fwd));
    if let Some(e3) = e3 {
        wire_edges.push(OrientedEdge::new(e3, true));
    }
    let wire = Wire::new(wire_edges, true)?;
    let wire_id = topo.add_wire(wire);

    let face = Face::new(wire_id, Vec::new(), stripe.surface.clone());
    let face_id = topo.add_face(face);

    Ok(BlendFaceInfo {
        face: face_id,
        cross_end: e1.map(|e| (e, v1e, v2e)),
        cross_start: e3.map(|e| (e, v2s, v1s)),
    })
}

/// A created blend face plus its two cross edges (the end cross-section
/// arcs), each with its (from, to) vertices in the blend wire's traversal
/// direction — the handles the end-cap notch surgery needs to SHARE those
/// arcs instead of leaving both sides use-1.
pub struct BlendFaceInfo {
    /// The blend face.
    pub face: FaceId,
    /// Cross edge at the spine end: `(edge, from, to)`.
    pub cross_end: Option<(brepkit_topology::edge::EdgeId, VertexId, VertexId)>,
    /// Cross edge at the spine start: `(edge, from, to)`. `None` when the
    /// stripe pinches to a point at that end and no cross edge exists.
    pub cross_start: Option<(brepkit_topology::edge::EdgeId, VertexId, VertexId)>,
}

/// Replace a face's two-edge corner path `from -> corner -> to` with the
/// single cross-section arc `edge`, notching the fillet's end profile out of
/// an end cap so the cap and the blend share one edge entity. Both replaced
/// edges must be straight (the box corner sides); returns whether a
/// replacement happened.
pub fn notch_face_corner_with_arc(
    topo: &mut Topology,
    face_id: FaceId,
    arc: (brepkit_topology::edge::EdgeId, VertexId, VertexId),
) -> Result<Option<FaceId>, BlendError> {
    let (arc_eid, va, vb) = arc;
    let wire_id = topo.face(face_id)?.outer_wire();
    let oes = topo.wire(wire_id)?.edges().to_vec();
    let n = oes.len();
    if n < 3 {
        return Ok(None);
    }
    let ends = |oe: &OrientedEdge| -> Result<(VertexId, VertexId), BlendError> {
        let e = topo.edge(oe.edge())?;
        Ok((oe.oriented_start(e), oe.oriented_end(e)))
    };
    if std::env::var("BK_NOTCH_TRACE").is_ok() {
        let mut has_a = false;
        let mut has_b = false;
        for oe in &oes {
            let (s, e) = ends(oe)?;
            has_a |= s == va || e == va;
            has_b |= s == vb || e == vb;
        }
        if has_a || has_b {
            log::warn!("NOTCH-TRACE face={face_id:?} has_va={has_a} has_vb={has_b} wire_len={n}");
        }
    }
    for i in 0..n {
        let j = (i + 1) % n;
        let (s0, e0) = ends(&oes[i])?;
        let (s1, e1) = ends(&oes[j])?;
        if e0 != s1 || e0 == va || e0 == vb {
            continue;
        }
        let fwd = s0 == va && e1 == vb;
        let rev = s0 == vb && e1 == va;
        if !(fwd || rev) {
            continue;
        }
        let both_straight = [oes[i].edge(), oes[j].edge()].iter().all(|&eid| {
            topo.edge(eid)
                .is_ok_and(|e| matches!(e.curve(), EdgeCurve::Line))
        });
        if !both_straight {
            continue;
        }
        let mut new_oes: Vec<OrientedEdge> = Vec::with_capacity(n - 1);
        for (k, oe) in oes.iter().enumerate() {
            if k == i {
                new_oes.push(OrientedEdge::new(arc_eid, fwd));
            } else if k != j {
                new_oes.push(*oe);
            }
        }
        let new_wire = topo.add_wire(Wire::new(new_oes, true)?);
        let (surface, reversed, inners) = {
            let f = topo.face(face_id)?;
            (
                f.surface().clone(),
                f.is_reversed(),
                f.inner_wires().to_vec(),
            )
        };
        let new_face = if reversed {
            Face::new_reversed(new_wire, inners, surface)
        } else {
            Face::new(new_wire, inners, surface)
        };
        let nf = topo.add_face(new_face);
        return Ok(Some(nf));
    }
    Ok(None)
}

/// Adapter that provides [`ParametricSurface`] for a `FaceSurface::Plane`.
///
/// Planes store only a normal and signed distance `d`, with no parametric
/// frame.  This adapter builds an orthonormal UV frame from the normal so
/// that the walking engine can evaluate, project, and differentiate the
/// plane surface uniformly.
pub struct PlaneAdapter {
    /// Origin point on the plane (the point closest to the world origin).
    pub origin: Point3,
    /// U-direction tangent (unit vector in the plane).
    pub u_dir: Vec3,
    /// V-direction tangent (unit vector in the plane, orthogonal to `u_dir`).
    pub v_dir: Vec3,
    /// Outward-facing unit normal.
    pub norm: Vec3,
}

impl PlaneAdapter {
    /// Build a `PlaneAdapter` from a plane normal and signed distance.
    ///
    /// The UV frame is constructed by choosing a non-parallel reference vector
    /// and computing the cross products.
    #[must_use]
    pub fn from_normal_and_d(normal: Vec3, d: f64) -> Self {
        let origin = Point3::new(normal.x() * d, normal.y() * d, normal.z() * d);

        // Pick a reference vector that is not parallel to the normal.
        let ref_vec = if normal.x().abs() < 0.9 {
            Vec3::new(1.0, 0.0, 0.0)
        } else {
            Vec3::new(0.0, 1.0, 0.0)
        };

        let u_dir = normal
            .cross(ref_vec)
            .normalize()
            .unwrap_or(Vec3::new(1.0, 0.0, 0.0));
        let v_dir = normal
            .cross(u_dir)
            .normalize()
            .unwrap_or(Vec3::new(0.0, 1.0, 0.0));

        Self {
            origin,
            u_dir,
            v_dir,
            norm: normal,
        }
    }
}

impl ParametricSurface for PlaneAdapter {
    fn evaluate(&self, u: f64, v: f64) -> Point3 {
        self.origin + self.u_dir * u + self.v_dir * v
    }

    fn normal(&self, _u: f64, _v: f64) -> Vec3 {
        self.norm
    }

    fn project_point(&self, point: Point3) -> (f64, f64) {
        let d = point - self.origin;
        (d.dot(self.u_dir), d.dot(self.v_dir))
    }

    fn partial_u(&self, _u: f64, _v: f64) -> Vec3 {
        self.u_dir
    }

    fn partial_v(&self, _u: f64, _v: f64) -> Vec3 {
        self.v_dir
    }
}

/// A [`ParametricSurface`] view that negates the wrapped surface's normal.
///
/// The walking engine's blend constraint places the rolling-ball centre on the
/// `+normal` side of each surface (`centre = p + r·normal`), so the surfaces
/// must present their **inward** (toward-material) normals. `PlaneAdapter`
/// flips a plane via its stored normal, but analytic/NURBS surfaces have an
/// intrinsic outward normal that can't be re-oriented in place — wrapping one
/// here flips it so a fillet against a curved neighbour solves the internal
/// (material-side) branch instead of the external common-tangent one.
pub struct FlippedNormalSurface<'a> {
    inner: &'a dyn ParametricSurface,
}

impl<'a> FlippedNormalSurface<'a> {
    /// Wrap a surface so its normal is negated.
    #[must_use]
    pub const fn new(inner: &'a dyn ParametricSurface) -> Self {
        Self { inner }
    }
}

impl ParametricSurface for FlippedNormalSurface<'_> {
    fn evaluate(&self, u: f64, v: f64) -> Point3 {
        self.inner.evaluate(u, v)
    }

    fn normal(&self, u: f64, v: f64) -> Vec3 {
        -self.inner.normal(u, v)
    }

    fn project_point(&self, point: Point3) -> (f64, f64) {
        self.inner.project_point(point)
    }

    fn partial_u(&self, u: f64, v: f64) -> Vec3 {
        self.inner.partial_u(u, v)
    }

    fn partial_v(&self, u: f64, v: f64) -> Vec3 {
        self.inner.partial_v(u, v)
    }
}

/// Extract a `&dyn ParametricSurface` from a `FaceSurface`, or build a
/// `PlaneAdapter` for plane faces.
///
/// Returns `Ok(adapter)` for planes and `Err(face_id)` for unsupported types.
/// For analytic and NURBS surfaces that already implement `ParametricSurface`,
/// the reference is extracted directly and the adapter is unused.
///
/// # Usage pattern
///
/// ```ignore
/// let mut adapter = None;
/// let surf: &dyn ParametricSurface = surface_ref_or_adapter(&face_surface, &mut adapter);
/// ```
#[must_use]
pub fn surface_ref_or_adapter<'a>(
    surface: &'a FaceSurface,
    adapter_slot: &'a mut Option<PlaneAdapter>,
) -> &'a dyn ParametricSurface {
    // For Plane faces, we need to populate the adapter_slot first,
    // then return a reference to it. For all other variants, we can
    // return a reference directly to the surface inside FaceSurface.
    if let FaceSurface::Plane { normal, d } = surface {
        let adapter = adapter_slot.insert(PlaneAdapter::from_normal_and_d(*normal, *d));
        return adapter as &dyn ParametricSurface;
    }
    match surface {
        FaceSurface::Plane { .. } => {
            // Already handled above; this arm is unreachable.
            adapter_slot.insert(PlaneAdapter::from_normal_and_d(
                Vec3::new(0.0, 0.0, 1.0),
                0.0,
            )) as &dyn ParametricSurface
        }
        FaceSurface::Cylinder(c) => c as &dyn ParametricSurface,
        FaceSurface::Cone(c) => c as &dyn ParametricSurface,
        FaceSurface::Sphere(s) => s as &dyn ParametricSurface,
        FaceSurface::Torus(t) => t as &dyn ParametricSurface,
        FaceSurface::Nurbs(n) => n as &dyn ParametricSurface,
    }
}

/// Weld pairs of free (use-1) edges that trace identical geometry.
///
/// Adjacent blend walls whose terminal sections coincide each mint their own
/// cross edge — same endpoints, same curve, two edge entities each used by
/// one face. Rewrite every wire of the faces to reference one edge per
/// geometric identity. Requires BOTH endpoints and the curve midpoint to
/// match at weld distance, so complementary arcs and genuinely distinct
/// co-endpoint edges are never merged; zero-length edges collapse away
/// entirely when their twin is also zero-length.
/// Split free Line edges whose interior contains another free edge's
/// endpoint (a full corner-edge segment coexisting with its two halves),
/// so the pieces become weldable, then fill CLOSED COPLANAR loops of
/// remaining free edges with an exact plane face — a closed free loop is a
/// hole, and the corner-floor triangles left by pairwise junction patches
/// are planar by construction.
#[allow(
    clippy::redundant_pub_crate,
    clippy::too_many_lines,
    clippy::items_after_statements,
    clippy::type_complexity
)]
pub(crate) fn close_residual_free_loops(
    topo: &mut Topology,
    faces: &mut Vec<FaceId>,
) -> Result<(), BlendError> {
    use brepkit_topology::edge::EdgeId;
    use std::collections::HashMap;

    let free_edges = |topo: &Topology, faces: &[FaceId]| -> Result<Vec<EdgeId>, BlendError> {
        let mut uses: HashMap<EdgeId, usize> = HashMap::new();
        for &fid in faces {
            let face = topo.face(fid)?;
            let mut wires = vec![face.outer_wire()];
            wires.extend_from_slice(face.inner_wires());
            for wid in wires {
                for oe in topo.wire(wid)?.edges() {
                    *uses.entry(oe.edge()).or_insert(0) += 1;
                }
            }
        }
        let mut v: Vec<EdgeId> = uses
            .iter()
            .filter(|&(_, &c)| c == 1)
            .map(|(&e, _)| e)
            .collect();
        v.sort_unstable_by_key(|e| e.index());
        Ok(v)
    };

    // Pass 1: split covering Line edges at interior endpoints of other free
    // edges, then re-weld.
    let frees = free_edges(topo, faces)?;
    let mut endpoints: Vec<(Point3, VertexId)> = Vec::new();
    for &eid in &frees {
        let e = topo.edge(eid)?;
        for v in [e.start(), e.end()] {
            let p = topo.vertex(v)?.point();
            if !endpoints.iter().any(|(q, _)| (*q - p).length() < 1e-6) {
                endpoints.push((p, v));
            }
        }
    }
    for &eid in &frees {
        let e = topo.edge(eid)?;
        if !matches!(e.curve(), EdgeCurve::Line) {
            continue;
        }
        let (sv, ev) = (e.start(), e.end());
        let sp = topo.vertex(sv)?.point();
        let ep = topo.vertex(ev)?.point();
        let dir = ep - sp;
        let len2 = dir.dot(dir);
        if len2 < 1e-18 {
            continue;
        }
        for (p, vid) in endpoints.clone() {
            let t = dir.dot(p - sp) / len2;
            if !(1e-6..=1.0 - 1e-6).contains(&t) {
                continue;
            }
            if (p - (sp + dir * t)).length() > 1e-6 {
                continue;
            }
            let oe = OrientedEdge::new(eid, true);
            let _ = crate::trimmer::split_edge_at(topo, &oe, vid)?;
            break;
        }
    }
    weld_coincident_free_edges(topo, faces)?;

    // Pass 2: fill closed coplanar loops of remaining free edges. Chains
    // connect POSITIONALLY (the loop's edges were minted by different
    // faces and share no vertex ids); the fill face is built from
    // geometry-identical COPIES with its own shared vertices, and the
    // final weld unifies each copy with its free original.
    let frees = free_edges(topo, faces)?;
    let mut used: std::collections::HashSet<EdgeId> = std::collections::HashSet::new();
    let ends_p = |topo: &Topology, eid: EdgeId| -> Result<(Point3, Point3), BlendError> {
        let e = topo.edge(eid)?;
        Ok((
            topo.vertex(e.start())?.point(),
            topo.vertex(e.end())?.point(),
        ))
    };
    let mut filled_any = false;
    for &seed in &frees {
        if used.contains(&seed) {
            continue;
        }
        let (s0, e0) = ends_p(topo, seed)?;
        let mut chain: Vec<(EdgeId, bool)> = vec![(seed, true)];
        let mut cursor = e0;
        let mut guard = 0;
        while (cursor - s0).length() > 1e-6 && guard < 8 {
            guard += 1;
            let mut advanced = false;
            for &c in &frees {
                if used.contains(&c) || chain.iter().any(|(x, _)| *x == c) {
                    continue;
                }
                let Ok((a, b)) = ends_p(topo, c) else {
                    continue;
                };
                if (a - cursor).length() <= 1e-6 {
                    cursor = b;
                    chain.push((c, true));
                    advanced = true;
                    break;
                }
                if (b - cursor).length() <= 1e-6 {
                    cursor = a;
                    chain.push((c, false));
                    advanced = true;
                    break;
                }
            }
            if !advanced {
                break;
            }
        }
        if (cursor - s0).length() > 1e-6 || chain.len() < 2 || chain.len() > 4 {
            continue;
        }
        // Loop corner positions in order, and coplanarity.
        let mut pts: Vec<Point3> = Vec::new();
        for &(eid, fwd) in &chain {
            let (a, b) = ends_p(topo, eid)?;
            pts.push(if fwd { a } else { b });
        }
        // A 2-edge loop (arc + chord lens: a band's straight rail against a
        // rebuilt face's bridge arc) has too few corners to span a plane;
        // take the plane from an arc's own circle instead.
        let nrm = if chain.len() == 2 {
            let mut circle_nrm = Option::None;
            for &(eid, _) in &chain {
                if let EdgeCurve::Circle(c) = topo.edge(eid)?.curve() {
                    let n = c.normal();
                    if circle_nrm.is_some_and(|prev: Vec3| prev.cross(n).length() > 1e-6) {
                        circle_nrm = Option::None;
                        break;
                    }
                    circle_nrm = Some(n);
                }
            }
            let Some(nrm) = circle_nrm else { continue };
            nrm
        } else {
            let n_raw = (pts[1] - pts[0]).cross(pts[2] - pts[0]);
            let Ok(nrm) = n_raw.normalize() else { continue };
            nrm
        };
        if pts.iter().any(|p| ((*p - pts[0]).dot(nrm)).abs() > 1e-6) {
            continue;
        }
        // Mint shared corner vertices and copy edges.
        let vids: Vec<VertexId> = pts
            .iter()
            .map(|&p| topo.add_vertex(Vertex::new(p, 1e-7)))
            .collect();
        let mut oes: Vec<OrientedEdge> = Vec::with_capacity(chain.len());
        let mut ok = true;
        for (k, &(eid, fwd)) in chain.iter().enumerate() {
            let curve = topo.edge(eid)?.curve().clone();
            let (v_from, v_to) = (vids[k], vids[(k + 1) % chain.len()]);
            let new_e = if fwd {
                topo.add_edge(Edge::new(v_from, v_to, curve))
            } else {
                topo.add_edge(Edge::new(v_to, v_from, curve))
            };
            if topo.edge(new_e).is_err() {
                ok = false;
                break;
            }
            oes.push(OrientedEdge::new(new_e, fwd));
        }
        if !ok {
            continue;
        }
        let Ok(wire) = Wire::new(oes, true) else {
            continue;
        };
        let wid = topo.add_wire(wire);
        let d = nrm.dot(Vec3::new(pts[0].x(), pts[0].y(), pts[0].z()));
        let fid = topo.add_face(Face::new(
            wid,
            Vec::new(),
            FaceSurface::Plane { normal: nrm, d },
        ));
        faces.push(fid);
        for &(eid, _) in &chain {
            used.insert(eid);
        }
        filled_any = true;
        log::debug!(
            "residual free loop filled with a plane face ({} edges)",
            chain.len()
        );
    }
    if filled_any {
        weld_coincident_free_edges(topo, faces)?;
    }
    Ok(())
}

#[allow(
    clippy::redundant_pub_crate,
    clippy::items_after_statements,
    clippy::type_complexity
)]
pub(crate) fn weld_coincident_free_edges(
    topo: &mut Topology,
    faces: &[FaceId],
) -> Result<(), BlendError> {
    use brepkit_topology::edge::EdgeId;
    use std::collections::HashMap;

    let mut uses: HashMap<EdgeId, usize> = HashMap::new();
    for &fid in faces {
        let face = topo.face(fid)?;
        let mut wires = vec![face.outer_wire()];
        wires.extend_from_slice(face.inner_wires());
        for wid in wires {
            for oe in topo.wire(wid)?.edges() {
                *uses.entry(oe.edge()).or_insert(0) += 1;
            }
        }
    }

    const WELD: f64 = 1e-6;
    let q = |p: Point3| -> (i64, i64, i64) {
        (
            (p.x() / WELD).round() as i64,
            (p.y() / WELD).round() as i64,
            (p.z() / WELD).round() as i64,
        )
    };

    // Geometry key for every free edge: symmetric endpoint pair + midpoint.
    let mut groups: HashMap<
        ((i64, i64, i64), (i64, i64, i64), (i64, i64, i64)),
        Vec<(EdgeId, VertexId, VertexId)>,
    > = HashMap::new();
    let mut free_edges: Vec<EdgeId> = uses
        .iter()
        .filter(|&(_, &c)| c == 1)
        .map(|(&e, _)| e)
        .collect();
    free_edges.sort_unstable_by_key(|e| e.index());
    for eid in free_edges {
        let e = topo.edge(eid)?;
        let (sv, ev) = (e.start(), e.end());
        let sp = topo.vertex(sv)?.point();
        let ep = topo.vertex(ev)?.point();
        // The geometric identity slot. Stored-curve evaluation is
        // phase-dependent for circles (endpoints do not trim the raw
        // parameterization), so circle edges key on centre + radius + |axis|
        // instead; antipodal endpoint pairs stay unkeyed (minor/major arc
        // ambiguity — the merge-key lesson).
        let mid = match e.curve() {
            EdgeCurve::Circle(c) => {
                let chord_mid = Point3::new(
                    (sp.x() + ep.x()) * 0.5,
                    (sp.y() + ep.y()) * 0.5,
                    (sp.z() + ep.z()) * 0.5,
                );
                if (chord_mid - c.center()).length() < 1e-6 {
                    continue;
                }
                let ax = c.normal();
                c.center() + Vec3::new(ax.x().abs(), ax.y().abs(), ax.z().abs()) * c.radius()
            }
            _ => e.curve().evaluate_with_endpoints(0.5, sp, ep),
        };
        let (ks, ke) = (q(sp), q(ep));
        let key = if ks <= ke {
            (ks, ke, q(mid))
        } else {
            (ke, ks, q(mid))
        };
        groups.entry(key).or_default().push((eid, sv, ev));
    }

    // For each group, rewrite all wires to use the first edge.
    let mut replace: HashMap<EdgeId, (EdgeId, bool)> = HashMap::new();
    for members in groups.values() {
        if members.len() < 2 {
            continue;
        }
        let (keep, keep_sv, _) = members[0];
        let keep_sp = topo.vertex(keep_sv)?.point();
        for &(dup, dup_sv, _) in &members[1..] {
            let dup_sp = topo.vertex(dup_sv)?.point();
            let same_dir = (dup_sp - keep_sp).length() < WELD;
            replace.insert(dup, (keep, same_dir));
        }
    }
    if replace.is_empty() {
        return Ok(());
    }

    for &fid in faces {
        let face = topo.face(fid)?;
        let mut wires = vec![face.outer_wire()];
        wires.extend_from_slice(face.inner_wires());
        for wid in wires {
            let wire = topo.wire(wid)?;
            let mut edges = wire.edges().to_vec();
            let mut changed = false;
            for oe in &mut edges {
                if let Some(&(keep, same_dir)) = replace.get(&oe.edge()) {
                    let fwd = if same_dir {
                        oe.is_forward()
                    } else {
                        !oe.is_forward()
                    };
                    *oe = OrientedEdge::new(keep, fwd);
                    changed = true;
                }
            }
            if changed {
                let closed = wire.is_closed();
                *topo.wire_mut(wid)? = Wire::new(edges, closed)?;
            }
        }
    }
    Ok(())
}

/// Make each new face's effective surface normal agree with the solid's
/// boundary-walk convention.
///
/// Meshers wind triangles by the effective normal (surface normal XOR
/// reversal) while manifold pairing runs on effective wire senses; a face
/// can satisfy sense pairing with a backwards normal, and its mesh then
/// comes out flipped against every neighbour. The interior-side integral
/// sums `(n x t) . (c - p)` along the effective boundary; its sign says
/// which side of the walk the interior lies on. The repair is the
/// sense-preserving triple flip: reverse the wire order, toggle every edge
/// sense, and toggle the reversal flag — effective senses are unchanged
/// (sense XOR reversal is invariant) while the effective normal flips.
#[allow(clippy::redundant_pub_crate, clippy::too_many_lines)]
pub(crate) fn normalize_face_normals(
    topo: &mut Topology,
    faces: &[FaceId],
    seeds: &[FaceId],
) -> Result<(), BlendError> {
    let trace = std::env::var("BK_NORM_TRACE").is_ok();
    // The boundary-walk convention is a property of the INPUT solid, not an
    // absolute: a solid built from a clockwise profile walks its boundaries
    // with the interior on the right, and every face of a valid solid obeys
    // ONE convention. Calibrate the expected sign from the carried-over
    // input faces, then repair only new faces that disagree.
    let seed_set: std::collections::HashSet<FaceId> = seeds.iter().copied().collect();
    let mut convention = 0.0;
    let mut flips: Vec<FaceId> = Vec::new();
    for &fid in seeds.iter().chain(faces.iter()) {
        let face = topo.face(fid)?;
        let rev = face.is_reversed();
        let surface = face.surface().clone();
        let wid = face.outer_wire();
        // The interior-side integral is only meaningful for disk-like
        // boundaries. A closed band (a full-revolution rim fillet's torus:
        // two rim circles joined by a doubled seam edge) or a face with
        // holes is skipped — and must be: the structured two-rim mesher
        // depends on the band's wire layout, which the triple flip would
        // rearrange.
        if !face.inner_wires().is_empty() {
            continue;
        }
        let wire = topo.wire(wid)?;
        {
            let mut seen = std::collections::HashSet::new();
            if wire.edges().iter().any(|oe| !seen.insert(oe.edge())) {
                continue;
            }
        }

        // Sample the outer boundary in EFFECTIVE traversal order: a
        // reversed face's boundary is the wire in reverse order with
        // flipped senses.
        let oes: Vec<_> = if rev {
            wire.edges().iter().rev().copied().collect()
        } else {
            wire.edges().to_vec()
        };
        let mut pts: Vec<Point3> = Vec::new();
        for oe in &oes {
            let e = topo.edge(oe.edge())?;
            let (sp, ep) = (
                topo.vertex(e.start())?.point(),
                topo.vertex(e.end())?.point(),
            );
            let (t0, t1) = e.curve().domain_with_endpoints(sp, ep);
            let n = 8usize;
            let fwd = oe.is_forward() ^ rev;
            for k in 0..n {
                #[allow(clippy::cast_precision_loss)]
                let f = k as f64 / n as f64;
                let t = if fwd {
                    t0 + (t1 - t0) * f
                } else {
                    t1 - (t1 - t0) * f
                };
                pts.push(e.curve().evaluate_with_endpoints(t, sp, ep));
            }
        }
        if pts.len() < 3 {
            continue;
        }
        let inv = 1.0 / {
            #[allow(clippy::cast_precision_loss)]
            let n = pts.len() as f64;
            n
        };
        let mut cx = 0.0;
        let mut cy = 0.0;
        let mut cz = 0.0;
        for p in &pts {
            cx += p.x();
            cy += p.y();
            cz += p.z();
        }
        let c = Point3::new(cx * inv, cy * inv, cz * inv);

        // Interior-left rule: walking the effective boundary with the
        // effective normal up, the face interior lies to the LEFT. The
        // accumulated test integral sums (n x t) . (c - p) over the
        // boundary; a negative total means the effective normal points the
        // wrong way. Unlike a fixed-normal winding (Newell) test, this
        // holds for VALID faces of either wire-winding convention — a
        // trimmed cap whose outwardness is encoded purely in the reversal
        // flag measures positive here, and must not be flipped.
        let normal_at = |p: Point3| -> Option<Vec3> {
            if let FaceSurface::Plane { normal, .. } = &surface {
                Some(*normal)
            } else {
                let (u, v) = surface.project_point(p)?;
                Some(surface.normal(u, v))
            }
        };
        let mut accum = 0.0;
        let mut total_len = 0.0;
        for (k, &a) in pts.iter().enumerate() {
            let b = pts[(k + 1) % pts.len()];
            let seg = b - a;
            let len = seg.length();
            if len < 1e-12 {
                continue;
            }
            let m = Point3::new(
                f64::midpoint(a.x(), b.x()),
                f64::midpoint(a.y(), b.y()),
                f64::midpoint(a.z(), b.z()),
            );
            let Some(mut n) = normal_at(m) else { continue };
            if rev {
                n = -n;
            }
            accum += n.cross(seg).dot(c - m);
            total_len += len;
        }
        if total_len < 1e-12 {
            continue;
        }
        if trace {
            log::debug!(
                "normalize: {fid:?} {} rev={rev} accum={accum:.4} seed={} c=({:.2},{:.2},{:.2})",
                surface.type_tag(),
                seed_set.contains(&fid),
                c.x(),
                c.y(),
                c.z()
            );
        }
        // Ignore slivers whose integral is numerically indecisive.
        if accum.abs() < 1e-9 * total_len * total_len {
            continue;
        }
        if seed_set.contains(&fid) {
            convention += accum.signum();
            continue;
        }
        if convention != 0.0 && accum.signum() != convention.signum() {
            flips.push(fid);
        }
    }

    for fid in &flips {
        let face = topo.face(*fid)?;
        let rev = face.is_reversed();
        let mut wires = vec![face.outer_wire()];
        wires.extend_from_slice(face.inner_wires());
        for wid in wires {
            let wire = topo.wire_mut(wid)?;
            let mut oes: Vec<_> = wire.edges().to_vec();
            oes.reverse();
            for oe in &mut oes {
                *oe = brepkit_topology::wire::OrientedEdge::new(oe.edge(), !oe.is_forward());
            }
            for (slot, oe) in wire.edges_mut().iter_mut().zip(oes) {
                *slot = oe;
            }
        }
        topo.face_mut(*fid)?.set_reversed(!rev);
    }
    if !flips.is_empty() {
        log::debug!(
            "normalize_face_normals: triple-flipped {} faces",
            flips.len()
        );
    }
    Ok(())
}

/// Propagate orientation consistency from `seeds` (faces whose orientation
/// is known-correct, typically untouched input faces) across shared edges.
///
/// A manifold shell's edges must each be traversed once forward and once
/// backward by their two using faces (effective sense = wire sense XOR face
/// reversal). Newly built blend faces (walls, corner patches, bands, fills)
/// and rebuilt originals pick their wire order constructively, so whole
/// faces can come out backwards relative to their neighbours; every such
/// face flips coherently via `set_reversed`, which fixes both its effective
/// wire senses and its effective surface normal (constructions wind their
/// wires consistently with their stored normals). Faces unreachable from
/// any seed (or in sense conflict, which a valid closed shell cannot
/// produce) are left as built.
#[allow(clippy::redundant_pub_crate)]
pub(crate) fn propagate_orientation(
    topo: &mut Topology,
    faces: &[FaceId],
    seeds: &[FaceId],
) -> Result<(), BlendError> {
    use brepkit_topology::edge::EdgeId;
    use std::collections::{HashMap, HashSet, VecDeque};

    // Raw wire senses are immutable during propagation; only the reversal
    // bit changes on a flip. Precompute per-face raw senses and track
    // reversal in a map so the BFS stays O(E) instead of re-walking wires
    // per visited edge.
    let mut raw_senses: HashMap<FaceId, Vec<(EdgeId, bool)>> = HashMap::new();
    let mut revs: HashMap<FaceId, bool> = HashMap::new();
    for &fid in faces {
        let face = topo.face(fid)?;
        revs.insert(fid, face.is_reversed());
        let mut v = Vec::new();
        let mut wires = vec![face.outer_wire()];
        wires.extend_from_slice(face.inner_wires());
        for wid in wires {
            for oe in topo.wire(wid)?.edges() {
                v.push((oe.edge(), oe.is_forward()));
            }
        }
        raw_senses.insert(fid, v);
    }

    let mut edge_users: HashMap<EdgeId, Vec<FaceId>> = HashMap::new();
    for (&fid, senses) in &raw_senses {
        for (eid, _) in senses {
            edge_users.entry(*eid).or_default().push(fid);
        }
    }

    let face_set: HashSet<FaceId> = faces.iter().copied().collect();
    let mut visited: HashSet<FaceId> = seeds
        .iter()
        .copied()
        .filter(|f| face_set.contains(f))
        .collect();
    let mut queue: VecDeque<FaceId> = visited.iter().copied().collect();
    // A shell can have no untouched face (fully consumed input); fall back
    // to the largest-index face as an arbitrary but deterministic seed.
    if queue.is_empty()
        && let Some(&f) = faces.iter().max()
    {
        visited.insert(f);
        queue.push_back(f);
    }

    let mut flipped = 0usize;
    while let Some(fid) = queue.pop_front() {
        let my_rev = revs.get(&fid).copied().unwrap_or(false);
        let senses = raw_senses.get(&fid).cloned().unwrap_or_default();
        for (eid, raw) in senses {
            let my_sense = raw ^ my_rev;
            let Some(users) = edge_users.get(&eid) else {
                continue;
            };
            if users.len() != 2 {
                continue;
            }
            for &other in users {
                if other == fid || visited.contains(&other) {
                    continue;
                }
                let other_rev = revs.get(&other).copied().unwrap_or(false);
                let Some(&(_, other_raw)) = raw_senses
                    .get(&other)
                    .and_then(|v| v.iter().find(|(e, _)| *e == eid))
                else {
                    continue;
                };
                if other_raw ^ other_rev == my_sense {
                    topo.face_mut(other)?.set_reversed(!other_rev);
                    revs.insert(other, !other_rev);
                    flipped += 1;
                }
                visited.insert(other);
                queue.push_back(other);
            }
        }
    }
    if flipped > 0 {
        log::debug!("propagate_orientation: flipped {flipped} faces");
    }
    Ok(())
}