BREP_kernel 0.2.0

use crate::fit::solve_dense;
use crate::topology::{BrepSolid, CoedgeRecord, EdgeRecord, FaceRecord, LoopRecord, VertexRecord};
use crate::{interpolate_curve, KnotVector, NurbsCurve, NurbsSurface, Vec3, Vec4};
use rustc_hash::FxHashMap as HashMap;
use serde::Serialize;

fn domains(surface: &NurbsSurface) -> Result<([f64; 2], [f64; 2]), String> {
    Ok((
        KnotVector::new(surface.knots_u.clone(), surface.degree_u)?.domain(),
        KnotVector::new(surface.knots_v.clone(), surface.degree_v)?.domain(),
    ))
}

fn stable_face_normal(face: &FaceRecord, u: f64, v: f64) -> Result<Vec3, String> {
    let normal_at = |u, v| face.surface.normal(u, v).ok();
    let mut normal = normal_at(u, v);
    let ([u0, u1], [v0, v1]) = domains(&face.surface)?;
    if normal.is_none() {
        let du = (u1 - u0) * 1e-5;
        let dv = (v1 - v0) * 1e-5;
        for (candidate_u, candidate_v) in [
            ((u + du).clamp(u0, u1), v),
            ((u - du).clamp(u0, u1), v),
            (u, (v + dv).clamp(v0, v1)),
            (u, (v - dv).clamp(v0, v1)),
        ] {
            normal = normal_at(candidate_u, candidate_v);
            if normal.is_some() {
                break;
            }
        }
    }
    // SINGULAR-ROW override: at a surface singularity where a whole
    // parameter row collapses to one point (a cone apex), the at-point /
    // nudged normal is the cross of a vanishing partial with noise — the
    // AXIS direction instead of the ruling normal, which offsets the apex
    // row straight down the axis and bends the fitted surface by exactly
    // d·cos(half-angle). Detect the collapse by local point spread, walk
    // DEEP inward for the true per-ruling limit, and replace the at-point
    // value only when the two genuinely DISAGREE. A sphere/dome pole also
    // reads as collapsed, but there the at-point normal (the axis) IS the
    // limit — agreement keeps the exact baseline value.
    let singular_here = {
        let du = (u1 - u0) * 1e-4;
        let dv = (v1 - v0) * 1e-4;
        let here = face.surface.evaluate(u, v)?;
        let along_u = face
            .surface
            .evaluate((u + du).clamp(u0, u1), v)?
            .sub(here)
            .length()
            .max(
                face.surface
                    .evaluate((u - du).clamp(u0, u1), v)?
                    .sub(here)
                    .length(),
            );
        let along_v = face
            .surface
            .evaluate(u, (v + dv).clamp(v0, v1))?
            .sub(here)
            .length()
            .max(
                face.surface
                    .evaluate(u, (v - dv).clamp(v0, v1))?
                    .sub(here)
                    .length(),
            );
        let scale = along_u.max(along_v);
        scale > 0.0 && along_u.min(along_v) < scale * 1e-6
    };
    if singular_here {
        let v_mid = (v0 + v1) * 0.5;
        let u_mid = (u0 + u1) * 0.5;
        let mut interior = None;
        for fraction in [1e-3, 1e-2, 5e-2, 0.25] {
            let candidate_v = v + (v_mid - v) * fraction;
            let candidate_u = u + (u_mid - u) * fraction;
            for (cu, cv) in [(u, candidate_v), (candidate_u, v), (candidate_u, candidate_v)] {
                if let Some(candidate) = normal_at(cu, cv) {
                    interior = Some(candidate);
                    break;
                }
            }
            if interior.is_some() {
                break;
            }
        }
        normal = match (normal, interior) {
            (Some(at_point), Some(interior)) if at_point.dot(interior) > 1.0 - 1e-6 => {
                Some(at_point)
            }
            (_, Some(interior)) => Some(interior),
            (at_point, None) => at_point,
        };
    }
    let normal =
        normal.ok_or_else(|| "offset_surface: cannot determine surface normal".to_string())?;
    Ok(if face.same_sense {
        normal
    } else {
        normal.scale(-1.0)
    })
}

fn greville_parameters(knots: &KnotVector) -> Vec<f64> {
    let mut parameters = (0..knots.control_point_count())
        .map(|index| {
            knots.knots[index + 1..=index + knots.degree]
                .iter()
                .sum::<f64>()
                / knots.degree as f64
        })
        .collect::<Vec<_>>();
    let domain = knots.domain();
    parameters[0] = domain[0];
    *parameters.last_mut().unwrap() = domain[1];
    parameters
}

/// Rational collocation matrix: rows are the rational basis functions
/// R_i(t) = N_i(t)·w_i / Σ_k N_k(t)·w_k evaluated at each parameter. With
/// uniform weights this reduces to the ordinary B-spline collocation matrix.
fn collocation_matrix(knots: &KnotVector, parameters: &[f64], weights: &[f64]) -> Vec<Vec<f64>> {
    parameters
        .iter()
        .map(|parameter| {
            let mut row = vec![0.0; knots.control_point_count()];
            let span = knots.find_span(*parameter);
            for (offset, value) in knots
                .basis_functions(span, *parameter)
                .into_iter()
                .enumerate()
            {
                let index = span - knots.degree + offset;
                row[index] = value * weights[index];
            }
            let denominator: f64 = row.iter().sum();
            if denominator.abs() > 0.0 {
                for value in &mut row {
                    *value /= denominator;
                }
            }
            row
        })
        .collect()
}

/// Split the weight grid into per-direction factors when it is separable
/// (w_ij = a_i·b_j), which covers every tensor surface built from rational
/// profile/rail curves (cylinders, cones, spheres, tori, revolves).
fn separable_weights(weights: &[Vec<f64>]) -> Option<(Vec<f64>, Vec<f64>)> {
    let first_row = weights.first()?;
    let anchor = *first_row.first()?;
    if anchor.abs() <= 1e-12 {
        return None;
    }
    let a: Vec<f64> = weights.iter().map(|row| row[0]).collect();
    let b: Vec<f64> = first_row.iter().map(|w| w / anchor).collect();
    for (i, row) in weights.iter().enumerate() {
        for (j, &w) in row.iter().enumerate() {
            if (w - a[i] * b[j]).abs() > 1e-10 * (1.0 + w.abs()) {
                return None;
            }
        }
    }
    Some((a, b))
}

/// Interpolate the sample grid in the SOURCE surface's rational basis (same
/// knots and weights). When the true offset is representable in that basis —
/// planes, cylinders, cones, spheres, tori — collocation at the Greville grid
/// recovers it EXACTLY, so offset carriers stay real analytic surfaces
/// instead of non-rational approximations with span-scale wobble.
fn interpolate_tensor(
    knot_u: &KnotVector,
    knot_v: &KnotVector,
    parameters_u: &[f64],
    parameters_v: &[f64],
    samples: &[Vec<Vec3>],
    weights: &[Vec<f64>],
) -> Result<Vec<Vec<Vec4>>, String> {
    let count_u = parameters_u.len();
    let count_v = parameters_v.len();
    if let Some((weights_u, weights_v)) = separable_weights(weights) {
        let matrix_u = collocation_matrix(knot_u, parameters_u, &weights_u);
        let matrix_v = collocation_matrix(knot_v, parameters_v, &weights_v);
        let mut intermediate = vec![vec![Vec3::default(); count_v]; count_u];
        for column in 0..count_v {
            let solve_axis = |axis: fn(Vec3) -> f64| {
                solve_dense(
                    matrix_u.clone(),
                    samples.iter().map(|row| axis(row[column])).collect(),
                )
            };
            let x = solve_axis(|point| point.x)?;
            let y = solve_axis(|point| point.y)?;
            let z = solve_axis(|point| point.z)?;
            for row in 0..count_u {
                intermediate[row][column] = Vec3::new(x[row], y[row], z[row]);
            }
        }
        let mut controls = vec![vec![Vec4::from_point(Vec3::default(), 1.0); count_v]; count_u];
        for row in 0..count_u {
            let solve_axis = |axis: fn(Vec3) -> f64| {
                solve_dense(
                    matrix_v.clone(),
                    intermediate[row].iter().copied().map(axis).collect(),
                )
            };
            let x = solve_axis(|point| point.x)?;
            let y = solve_axis(|point| point.y)?;
            let z = solve_axis(|point| point.z)?;
            for column in 0..count_v {
                controls[row][column] = Vec4::from_point(
                    Vec3::new(x[column], y[column], z[column]),
                    weights[row][column],
                );
            }
        }
        return Ok(controls);
    }

    // Non-separable weights: solve the full tensor collocation system with
    // the exact 2D rational basis. Nets are small in practice.
    let unknowns = count_u * count_v;
    let mut matrix = vec![vec![0.0; unknowns]; unknowns];
    for (k, &u) in parameters_u.iter().enumerate() {
        let span_u = knot_u.find_span(u);
        let basis_u = knot_u.basis_functions(span_u, u);
        for (l, &v) in parameters_v.iter().enumerate() {
            let span_v = knot_v.find_span(v);
            let basis_v = knot_v.basis_functions(span_v, v);
            let row = &mut matrix[k * count_v + l];
            let mut denominator = 0.0;
            for (du, value_u) in basis_u.iter().enumerate() {
                let i = span_u - knot_u.degree + du;
                for (dv, value_v) in basis_v.iter().enumerate() {
                    let j = span_v - knot_v.degree + dv;
                    let entry = value_u * value_v * weights[i][j];
                    row[i * count_v + j] = entry;
                    denominator += entry;
                }
            }
            if denominator.abs() > 0.0 {
                for value in row.iter_mut() {
                    *value /= denominator;
                }
            }
        }
    }
    let solve_axis = |axis: fn(Vec3) -> f64| {
        solve_dense(
            matrix.clone(),
            samples
                .iter()
                .flat_map(|row| row.iter().copied().map(axis))
                .collect(),
        )
    };
    let x = solve_axis(|point| point.x)?;
    let y = solve_axis(|point| point.y)?;
    let z = solve_axis(|point| point.z)?;
    let mut controls = vec![vec![Vec4::from_point(Vec3::default(), 1.0); count_v]; count_u];
    for row in 0..count_u {
        for column in 0..count_v {
            let index = row * count_v + column;
            controls[row][column] = Vec4::from_point(
                Vec3::new(x[index], y[index], z[index]),
                weights[row][column],
            );
        }
    }
    Ok(controls)
}

/// Construct the same fitted offset carrier surface as the reference shell
/// implementation. Positive distance follows its convention and moves
/// opposite the face's outward normal.
pub fn offset_surface(
    face: &FaceRecord,
    distance: f64,
    planar_extension: f64,
) -> Result<NurbsSurface, String> {
    let source = &face.surface;
    if source.is_affine()? {
        let ([u0, u1], [v0, v1]) = domains(source)?;
        let normal = stable_face_normal(face, (u0 + u1) / 2.0, (v0 + v1) / 2.0)?;
        let shift = normal.scale(-distance);
        let mut points = source
            .control_points
            .iter()
            .map(|row| {
                row.iter()
                    .map(|control| Ok(control.point()?.add(shift)))
                    .collect::<Result<Vec<_>, String>>()
            })
            .collect::<Result<Vec<_>, String>>()?;
        if planar_extension > 0.0 {
            let p00 = points[0][0];
            let p01 = points[0][1];
            let p10 = points[1][0];
            let direction_u = p10.sub(p00).normalized()?;
            let direction_v = p01.sub(p00).normalized()?;
            points[0][0] = p00
                .sub(direction_u.scale(planar_extension))
                .sub(direction_v.scale(planar_extension));
            points[0][1] = p01
                .sub(direction_u.scale(planar_extension))
                .add(direction_v.scale(planar_extension));
            points[1][0] = p10
                .add(direction_u.scale(planar_extension))
                .sub(direction_v.scale(planar_extension));
            points[1][1] = points[1][1]
                .add(direction_u.scale(planar_extension))
                .add(direction_v.scale(planar_extension));
        }
        let controls = points
            .into_iter()
            .enumerate()
            .map(|(row, points)| {
                points
                    .into_iter()
                    .enumerate()
                    .map(|(column, point)| {
                        Vec4::from_point(point, source.control_points[row][column].w)
                    })
                    .collect()
            })
            .collect();
        return NurbsSurface::new(
            source.degree_u,
            source.degree_v,
            source.knots_u.clone(),
            source.knots_v.clone(),
            controls,
        );
    }

    let knot_u = KnotVector::new(source.knots_u.clone(), source.degree_u)?;
    let knot_v = KnotVector::new(source.knots_v.clone(), source.degree_v)?;
    let parameters_u = greville_parameters(&knot_u);
    let parameters_v = greville_parameters(&knot_v);
    let mut samples = Vec::new();
    for &u in &parameters_u {
        let mut row = Vec::new();
        for &v in &parameters_v {
            row.push(
                source
                    .evaluate(u, v)?
                    .add(stable_face_normal(face, u, v)?.scale(-distance)),
            );
        }
        samples.push(row);
    }
    // APEX-CONE PINCH RETRIM: offsetting an apex cone INWARD moves each
    // ruling past the axis — the sampled far row becomes a ring on the far
    // side (radius d·cos half-angle, mirrored through the axis) and the
    // offset surface self-pinches inside the v-domain. The genuine cavity
    // ends AT the pinch (the offset cone's own apex). For a linear-v net
    // (two sample rows — every made/booleaned cone) the pinch lies on each
    // ruling at the fraction where the radial vector vanishes: detect the
    // inversion (far-row radials anti-parallel to near-row radials about the
    // row centroids) and pull the far row back to the pinch point, so the
    // fitted surface ends in a proper degenerate apex row instead of a
    // parasitic inverted tip ending in an unweldable ring.
    if parameters_v.len() == 2 && parameters_u.len() >= 3 {
        let centroid = |column: usize| {
            let mut sum = Vec3::default();
            for row in &samples {
                sum = sum.add(row[column]);
            }
            sum.scale(1.0 / samples.len() as f64)
        };
        let near_centroid = centroid(0);
        let far_centroid = centroid(1);
        let mut inverted = true;
        let mut pinch_fraction = 0.0f64;
        let mut near_mean = 0.0f64;
        let mut far_mean = 0.0f64;
        for row in &samples {
            let near_radial = row[0].sub(near_centroid);
            let far_radial = row[1].sub(far_centroid);
            let near_len = near_radial.length();
            let far_len = far_radial.length();
            if near_len <= 1e-9 || far_len <= 1e-9 {
                inverted = false;
                break;
            }
            if near_radial.dot(far_radial) >= 0.0 {
                inverted = false;
                break;
            }
            pinch_fraction += near_len / (near_len + far_len) / samples.len() as f64;
            near_mean += near_len / samples.len() as f64;
            far_mean += far_len / samples.len() as f64;
        }
        if inverted {
            // Pull the crossed (smaller-ring, past-the-pinch) end back to the
            // pinch point on each ruling.
            let retrim_far = far_mean <= near_mean;
            for row in &mut samples {
                let near = row[0];
                let far = row[1];
                let pinch = near.add(far.sub(near).scale(pinch_fraction));
                if retrim_far {
                    row[1] = pinch;
                } else {
                    row[0] = pinch;
                }
            }
        }
        // RULED EXTENSION: `planar_extension` is a no-op for curved carriers
        // above, but a cone/cylinder lateral joined at a reflex edge needs
        // its offset skin to GROW past the source rim exactly like a plane
        // (a cylinder piercing a cone: the two offsets only meet past both
        // cloned rims). A linear-v net is ruled — stretching each sampled
        // ruling beyond both ends stays ON the same surface, so the fitted
        // carrier keeps its parameterization (knots/pcurves untouched) while
        // its world image (and with it the cloned trim's image) inflates.
        if planar_extension > 0.0 && !inverted {
            let mut min_ruling = f64::MAX;
            let mut back_allowance = f64::MAX;
            let mut forward_allowance = f64::MAX;
            let mut extendable = true;
            for row in &samples {
                let ruling = row[1].sub(row[0]);
                let length = ruling.length();
                min_ruling = min_ruling.min(length);
                // Radii about the row centroids expose a converging (conic)
                // ruling sheaf; the extension must stop short of its apex or
                // the sheet folds through it.
                let near_radial = row[0].sub(near_centroid).length();
                let far_radial = row[1].sub(far_centroid).length();
                if (far_radial - near_radial).abs() > 1e-9 {
                    let apex_at = near_radial / (near_radial - far_radial);
                    if (-1e-9..=1.0 + 1e-9).contains(&apex_at) {
                        // Apex inside the span: degenerate sheet, do not touch.
                        extendable = false;
                        break;
                    }
                    if apex_at < 0.0 {
                        back_allowance = back_allowance.min(0.9 * -apex_at);
                    } else {
                        forward_allowance = forward_allowance.min(0.9 * (apex_at - 1.0));
                    }
                }
            }
            if extendable && min_ruling > 1e-9 {
                let stretch = planar_extension / min_ruling;
                let back = stretch.min(back_allowance);
                let forward = stretch.min(forward_allowance);
                for row in &mut samples {
                    let ruling = row[1].sub(row[0]);
                    row[0] = row[0].sub(ruling.scale(back));
                    row[1] = row[1].add(ruling.scale(forward));
                }
            }
        }
    }
    let weights = source
        .control_points
        .iter()
        .map(|row| row.iter().map(|point| point.w).collect::<Vec<_>>())
        .collect::<Vec<_>>();
    NurbsSurface::new(
        source.degree_u,
        source.degree_v,
        source.knots_u.clone(),
        source.knots_v.clone(),
        interpolate_tensor(
            &knot_u,
            &knot_v,
            &parameters_u,
            &parameters_v,
            &samples,
            &weights,
        )?,
    )
}

fn mapped_pcurve_polyline(
    surface: &NurbsSurface,
    pcurve: &NurbsCurve,
    degenerate: bool,
) -> Result<(Vec<Vec3>, Vec<f64>), String> {
    let [start, end] = pcurve.domain()?;
    let evaluate = |fraction: f64| {
        let uv = pcurve.evaluate(start + (end - start) * fraction)?;
        surface.evaluate(uv.x, uv.y)
    };
    let first = evaluate(0.0)?;
    let last = evaluate(1.0)?;
    if degenerate {
        return Ok((vec![first, last], vec![0.0, 1.0]));
    }
    fn append(
        evaluate: &impl Fn(f64) -> Result<Vec3, String>,
        a_fraction: f64,
        a: Vec3,
        b_fraction: f64,
        b: Vec3,
        depth: usize,
        parameters: &mut Vec<f64>,
        points: &mut Vec<Vec3>,
    ) -> Result<(), String> {
        let fractions =
            [0.25, 0.5, 0.75].map(|local| a_fraction + (b_fraction - a_fraction) * local);
        let samples = fractions
            .map(evaluate)
            .into_iter()
            .collect::<Result<Vec<_>, String>>()?;
        let deviation = samples
            .iter()
            .enumerate()
            .map(|(index, point)| {
                point
                    .sub(a.add(b.sub(a).scale((index + 1) as f64 * 0.25)))
                    .length()
            })
            .fold(0.0, f64::max);
        if deviation <= 5e-4 || depth >= 10 {
            parameters.push(b_fraction);
            points.push(b);
            return Ok(());
        }
        append(
            evaluate,
            a_fraction,
            a,
            fractions[1],
            samples[1],
            depth + 1,
            parameters,
            points,
        )?;
        append(
            evaluate,
            fractions[1],
            samples[1],
            b_fraction,
            b,
            depth + 1,
            parameters,
            points,
        )
    }
    let mut parameters = vec![0.0];
    let mut points = vec![first];
    append(
        &evaluate,
        0.0,
        first,
        1.0,
        last,
        0,
        &mut parameters,
        &mut points,
    )?;
    Ok((points, parameters))
}

#[derive(Clone, Debug, Serialize)]
pub struct OffsetFaceCarrier {
    pub vertices: Vec<VertexRecord>,
    pub edges: Vec<EdgeRecord>,
    pub face: FaceRecord,
}

fn claim_vertex_image(
    source_id: u64,
    point: Vec3,
    vertex_images: &mut HashMap<u64, u64>,
    vertices: &mut Vec<VertexRecord>,
    next_id: &mut u64,
) -> u64 {
    if let Some(id) = vertex_images.get(&source_id) {
        return *id;
    }
    let id = *next_id;
    *next_id += 1;
    vertices.push(VertexRecord { id, point });
    vertex_images.insert(source_id, id);
    id
}

pub fn offset_face_carrier(
    solid: &BrepSolid,
    face_id: u64,
    distance: f64,
    planar_extension: f64,
) -> Result<OffsetFaceCarrier, String> {
    let source = solid
        .shells
        .iter()
        .flat_map(|shell| &shell.faces)
        .find(|face| face.id == face_id)
        .ok_or_else(|| format!("offset_face_carrier: missing face {face_id}"))?;
    let surface = offset_surface(source, distance, planar_extension)?;
    let source_edges = solid
        .edges
        .iter()
        .map(|edge| (edge.id, edge))
        .collect::<HashMap<_, _>>();
    let source_vertices = solid
        .vertices
        .iter()
        .map(|vertex| (vertex.id, vertex))
        .collect::<HashMap<_, _>>();
    let mut vertices = Vec::new();
    let mut vertex_images = HashMap::default();
    let mut edges = Vec::new();
    let mut edge_images = HashMap::default();
    let mut loops = Vec::new();
    let mut next_id = 1u64;

    for source_loop in &source.loops {
        let mut coedges = Vec::new();
        for source_coedge in &source_loop.coedges {
            let source_edge = source_edges
                .get(&source_coedge.edge_id)
                .ok_or_else(|| "offset_face_carrier: missing source edge".to_string())?;
            let (source_start, source_end) = if source_coedge.forward {
                (source_edge.start_vertex_id, source_edge.end_vertex_id)
            } else {
                (source_edge.end_vertex_id, source_edge.start_vertex_id)
            };
            if !source_vertices.contains_key(&source_start)
                || !source_vertices.contains_key(&source_end)
            {
                return Err("offset_face_carrier: missing source vertex".into());
            }
            // Map even DEGENERATE source edges through the full polyline: a
            // cone apex's image on the offset surface is a genuine CIRCLE
            // (radius d·cos half-angle), not a point — shortcutting to the
            // two endpoints would collapse the ring and leave the carrier's
            // topology inconsistent with its surface. Edges whose image truly
            // collapses (sphere poles, planar corners) still interpolate to a
            // point-sized curve and keep their degenerate flag below.
            // Map even DEGENERATE source edges through the full polyline: an
            // EXTERIOR cone offset turns the apex point into a genuine RING
            // (radius d·cos half-angle) — shortcutting to the endpoints would
            // collapse it and leave the carrier topology inconsistent with
            // its surface (and the ring imprint would be dropped as
            // boundary-coincident with a "degenerate" edge). Images that
            // truly collapse (sphere poles; interior apexes after the pinch
            // retrim) stay degenerate below. The threshold scales with the
            // offset distance: a real ring measures ~d·cos α, while fitted
            // pole rows wobble ~1e-4 absolute.
            let (points, parameters) =
                mapped_pcurve_polyline(&surface, &source_coedge.pcurve, false)?;
            let collapse_tolerance = 1e-6f64.max(distance.abs() * 1e-2);
            let image_collapsed = points
                .iter()
                .all(|point| point.sub(points[0]).length() <= collapse_tolerance);
            let (edge_id, forward) =
                if let Some((edge_id, edge_start_vertex_id)) = edge_images.get(&source_edge.id) {
                    (
                        *edge_id,
                        vertex_images.get(&source_start) == Some(edge_start_vertex_id),
                    )
                } else {
                    let curve = if image_collapsed {
                        NurbsCurve::new(
                            1,
                            vec![0.0, 0.0, 1.0, 1.0],
                            vec![
                                Vec4::from_point(points[0], 1.0),
                                Vec4::from_point(points[0], 1.0),
                            ],
                        )?
                    } else {
                        interpolate_curve(&points, 1, &parameters)?
                    };
                    let start_vertex_id = claim_vertex_image(
                        source_start,
                        points[0],
                        &mut vertex_images,
                        &mut vertices,
                        &mut next_id,
                    );
                    let end_vertex_id = claim_vertex_image(
                        source_end,
                        points[points.len() - 1],
                        &mut vertex_images,
                        &mut vertices,
                        &mut next_id,
                    );
                    let id = next_id;
                    next_id += 1;
                    let domain = curve.domain()?;
                    edges.push(EdgeRecord {
                        id,
                        curve,
                        t0: domain[0],
                        t1: domain[1],
                        start_vertex_id,
                        end_vertex_id,
                        // Degenerate only if the IMAGE collapsed too — a cone
                        // apex maps to a real ring on the offset surface and
                        // must carry a real closed edge.
                        degenerate: source_edge.degenerate && image_collapsed,
                        // Image of a named source edge on the offset carrier;
                        // suffixed so it cannot collide with the source edge
                        // when both faces survive into one solid.
                        name: source_edge
                            .name
                            .as_ref()
                            .map(|name| format!("{name}_Offset")),
                    });
                    edge_images.insert(source_edge.id, (id, start_vertex_id));
                    (id, true)
                };
            let id = next_id;
            next_id += 1;
            coedges.push(CoedgeRecord {
                id,
                edge_id,
                forward,
                pcurve: source_coedge.pcurve.clone(),
            });
        }
        let id = next_id;
        next_id += 1;
        loops.push(LoopRecord { id, coedges });
    }
    Ok(OffsetFaceCarrier {
        vertices,
        edges,
        face: FaceRecord {
            id: next_id,
            surface,
            same_sense: source.same_sense,
            loops,
            name: source.name.as_ref().map(|name| format!("{name}_Offset")),
        },
    })
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::{make_box_brep, make_cylinder_brep};

    #[test]
    fn affine_offset_is_exact_and_preserves_weights() {
        let solid = make_box_brep(Vec3::default(), 4.0, 4.0, 4.0).unwrap();
        let face = &solid.shells[0].faces[0];
        let offset = offset_surface(face, 0.75, 0.0).unwrap();
        let domain_u = KnotVector::new(face.surface.knots_u.clone(), 1)
            .unwrap()
            .domain();
        let domain_v = KnotVector::new(face.surface.knots_v.clone(), 1)
            .unwrap()
            .domain();
        let u = (domain_u[0] + domain_u[1]) / 2.0;
        let v = (domain_v[0] + domain_v[1]) / 2.0;
        let displacement = offset
            .evaluate(u, v)
            .unwrap()
            .sub(face.surface.evaluate(u, v).unwrap());
        assert!((displacement.length() - 0.75).abs() < 1e-12);
    }

    #[test]
    fn curved_offset_carrier_maps_every_trim_to_new_surface() {
        let solid =
            make_cylinder_brep(Vec3::default(), Vec3::new(0.0, 0.0, 1.0), 2.0, 4.0).unwrap();
        let side = &solid.shells[0].faces[0];
        let carrier = offset_face_carrier(&solid, side.id, 0.5, 0.0).unwrap();
        for coedge in carrier
            .face
            .loops
            .iter()
            .flat_map(|loop_record| &loop_record.coedges)
        {
            let edge = carrier
                .edges
                .iter()
                .find(|edge| edge.id == coedge.edge_id)
                .unwrap();
            for fraction in [0.0, 0.3, 0.8, 1.0] {
                let uv = coedge.pcurve.evaluate(fraction).unwrap();
                let on_surface = carrier.face.surface.evaluate(uv.x, uv.y).unwrap();
                let parameter = if coedge.forward {
                    edge.t0 + (edge.t1 - edge.t0) * fraction
                } else {
                    edge.t1 - (edge.t1 - edge.t0) * fraction
                };
                assert!(
                    on_surface
                        .sub(edge.curve.evaluate(parameter).unwrap())
                        .length()
                        < 7e-4
                );
            }
        }
    }
}