mesh-sieve 4.0.2

//! Construct intermediate entities (edges/faces) from cell→vertex meshes.
//!
//! # Expected cell types
//! This module assumes cells reference their vertices in a standard ordering:
//!
//! - [`CellType::Triangle`]: `(0,1,2)`
//! - [`CellType::Quadrilateral`]: `(0,1,2,3)`
//! - [`CellType::Tetrahedron`]: `(0,1,2,3)`
//! - [`CellType::Hexahedron`]: `(0,1,2,3,4,5,6,7)` with
//!   `0..=3` the bottom face and `4..=7` the top face.
//! - [`CellType::Polygon(n)`]: `(0..n)` in cyclic order.
//! - [`CellType::Polyhedron`]: requires explicit face connectivity metadata.
//!
//! The interpolation routine uses these orderings to derive edges and faces,
//! then inserts the appropriate arrows so the sieve has
//! `cell → face → edge → vertex` (for 3D cells) or `cell → edge → vertex`
//! (for 2D cells) connectivity.
//!
//! # Example
//! ```rust
//! # fn try_main() -> Result<(), mesh_sieve::mesh_error::MeshSieveError> {
//! use mesh_sieve::algs::interpolate::interpolate_edges_faces;
//! use mesh_sieve::data::atlas::Atlas;
//! use mesh_sieve::data::section::Section;
//! use mesh_sieve::data::storage::VecStorage;
//! use mesh_sieve::topology::cell_type::CellType;
//! use mesh_sieve::topology::point::PointId;
//! use mesh_sieve::topology::sieve::{InMemorySieve, MutableSieve, Sieve};
//!
//! let mut sieve = InMemorySieve::<PointId, ()>::default();
//! let cell = PointId::new(1)?;
//! let v = |i| PointId::new(i).unwrap();
//! for p in [cell, v(2), v(3), v(4), v(5)] {
//!     MutableSieve::add_point(&mut sieve, p);
//! }
//! for vert in [v(2), v(3), v(4), v(5)] {
//!     sieve.add_arrow(cell, vert, ());
//! }
//!
//! let mut atlas = Atlas::default();
//! for p in [cell, v(2), v(3), v(4), v(5)] {
//!     atlas.try_insert(p, 1)?;
//! }
//! let mut cell_types = Section::<CellType, VecStorage<CellType>>::new(atlas);
//! cell_types.try_set(cell, &[CellType::Quadrilateral])?;
//! for p in [v(2), v(3), v(4), v(5)] {
//!     cell_types.try_set(p, &[CellType::Vertex])?;
//! }
//!
//! let result = interpolate_edges_faces(&mut sieve, &mut cell_types)?;
//! assert_eq!(result.edge_points.len(), 4);
//! # Ok(())
//! # }
//! ```

use crate::data::coordinates::Coordinates;
use crate::data::section::Section;
use crate::data::storage::Storage;
use crate::geometry::cell_vertex_ordering::topology_cell_vertices;
use crate::mesh_error::MeshSieveError;
use crate::topology::cell_type::CellType;
use crate::topology::point::PointId;
use crate::topology::sieve::{MutableSieve, OrientedSieve};
use std::collections::BTreeMap;

/// Edge/face bookkeeping produced during interpolation.
#[derive(Debug, Default)]
pub struct InterpolationResult {
    /// Canonical edge key `(min,max)` → edge point.
    pub edge_points: BTreeMap<(PointId, PointId), PointId>,
    /// Canonical face key (sorted vertices) → face point.
    pub face_points: BTreeMap<Vec<PointId>, PointId>,
    /// First-seen oriented loop for each shared face.
    face_directions: BTreeMap<Vec<PointId>, Vec<PointId>>,
}

/// Optional per-cell ordering/connectivity metadata for interpolation.
#[derive(Debug, Default)]
pub struct InterpolationOrdering {
    /// Per-cell permutation of the cell's cone vertices.
    pub vertex_permutations: BTreeMap<PointId, Vec<usize>>,
    /// Per-cell face connectivity for polyhedra, expressed as vertex indices.
    pub polyhedron_faces: BTreeMap<PointId, Vec<Vec<usize>>>,
}

/// Insert intermediate edges/faces into a cell→vertex mesh.
pub fn interpolate_edges_faces<S, CtSt>(
    sieve: &mut S,
    cell_types: &mut Section<CellType, CtSt>,
) -> Result<InterpolationResult, MeshSieveError>
where
    S: MutableSieve<Point = PointId> + OrientedSieve<Orient = i32>,
    S::Payload: Default,
    CtSt: Storage<CellType> + Clone,
{
    interpolate_edges_faces_with_ordering(sieve, cell_types, None)
}

/// Insert intermediate edges/faces into a cell→vertex mesh with explicit ordering metadata.
pub fn interpolate_edges_faces_with_ordering<S, CtSt>(
    sieve: &mut S,
    cell_types: &mut Section<CellType, CtSt>,
    ordering: Option<&InterpolationOrdering>,
) -> Result<InterpolationResult, MeshSieveError>
where
    S: MutableSieve<Point = PointId> + OrientedSieve<Orient = i32>,
    S::Payload: Default,
    CtSt: Storage<CellType> + Clone,
{
    let mut max_id = 0u64;
    for p in sieve.points() {
        max_id = max_id.max(p.get());
    }
    let mut next_id = max_id
        .checked_add(1)
        .ok_or(MeshSieveError::InvalidPointId)?;

    let cells = collect_cells(cell_types)?;
    let mut result = InterpolationResult::default();
    let preserve_orientation = sieve.supports_nontrivial_orientation();

    for (cell, cell_type) in cells {
        let permutation = ordering.and_then(|meta| {
            meta.vertex_permutations
                .get(&cell)
                .map(|perm| perm.as_slice())
        });
        let topology_vertices = topology_cell_vertices(sieve, cell, cell_type)?;
        let vertices = if let Some(permutation) = permutation {
            apply_permutation(cell, &topology_vertices, permutation)?
        } else {
            topology_vertices
        };

        let poly_faces = ordering.and_then(|meta| {
            meta.polyhedron_faces
                .get(&cell)
                .map(|faces| faces.as_slice())
        });
        let edges = cell_edges(cell_type, &vertices, poly_faces)?;
        for [a, b] in edges {
            let edge = edge_point(
                &mut result.edge_points,
                &mut next_id,
                sieve,
                cell_types,
                a,
                b,
            )?;
            sieve.add_arrow_o(
                cell,
                edge,
                S::Payload::default(),
                if preserve_orientation {
                    edge_orientation(&vertices, a, b)
                } else {
                    0
                },
            )?;
        }

        for (face_vertices, face_type) in cell_faces(cell_type, &vertices, poly_faces)? {
            let (face, canonical_face_vertices) = face_point(
                &mut result.face_points,
                &mut result.face_directions,
                &mut next_id,
                sieve,
                cell_types,
                &face_vertices,
                face_type,
            )?;
            for [a, b] in face_edges(&canonical_face_vertices)? {
                let edge = edge_point(
                    &mut result.edge_points,
                    &mut next_id,
                    sieve,
                    cell_types,
                    a,
                    b,
                )?;
                sieve.add_arrow_o(
                    face,
                    edge,
                    S::Payload::default(),
                    if preserve_orientation {
                        edge_orientation(&canonical_face_vertices, a, b)
                    } else {
                        0
                    },
                )?;
            }
            sieve.add_arrow_o(
                cell,
                face,
                S::Payload::default(),
                if preserve_orientation {
                    face_orientation(&face_vertices, &canonical_face_vertices)
                } else {
                    0
                },
            )?;
        }
    }

    Ok(result)
}

fn edge_orientation(vertices: &[PointId], a: PointId, b: PointId) -> i32 {
    let ia = vertices.iter().position(|p| *p == a);
    let ib = vertices.iter().position(|p| *p == b);
    match (ia, ib) {
        (Some(a), Some(b)) if a < b => 0,
        (Some(_), Some(_)) => -1,
        _ => 0,
    }
}

/// Interpolate a mesh using the shared topology-and-coordinate canonical
/// ordering.  This variant is suitable for geometry-sensitive pipelines;
/// [`interpolate_edges_faces`] remains available for topology-only callers.
pub fn interpolate_edges_faces_with_coordinates<S, CtSt, CoordSt>(
    sieve: &mut S,
    cell_types: &mut Section<CellType, CtSt>,
    coordinates: &Coordinates<f64, CoordSt>,
) -> Result<InterpolationResult, MeshSieveError>
where
    S: MutableSieve<Point = PointId> + OrientedSieve<Orient = i32>,
    S::Payload: Default,
    CtSt: Storage<CellType> + Clone,
    CoordSt: Storage<f64>,
{
    let max_id = sieve.points().map(|p| p.get()).max().unwrap_or(0);
    let mut next_id = max_id
        .checked_add(1)
        .ok_or(MeshSieveError::InvalidPointId)?;
    let cells = collect_cells(cell_types)?;
    let mut result = InterpolationResult::default();
    let preserve_orientation = sieve.supports_nontrivial_orientation();
    for (cell, cell_type) in cells {
        let canonical = crate::geometry::cell_vertex_ordering::canonical_cell_vertices(
            sieve,
            cell,
            cell_type,
            coordinates,
        )?;
        let vertices = canonical.vertices;
        let edges = cell_edges(cell_type, &vertices, None)?;
        for [a, b] in edges {
            let edge = edge_point(
                &mut result.edge_points,
                &mut next_id,
                sieve,
                cell_types,
                a,
                b,
            )?;
            sieve.add_arrow_o(
                cell,
                edge,
                S::Payload::default(),
                if preserve_orientation {
                    edge_orientation(&vertices, a, b)
                } else {
                    0
                },
            )?;
        }
        for (face_vertices, face_type) in cell_faces(cell_type, &vertices, None)? {
            let local_face_vertices = face_vertices;
            let face_vertices = canonical_face_loop(&local_face_vertices, face_type, coordinates)?;
            let (face, canonical_face_vertices) = face_point(
                &mut result.face_points,
                &mut result.face_directions,
                &mut next_id,
                sieve,
                cell_types,
                &face_vertices,
                face_type,
            )?;
            for [a, b] in face_edges(&canonical_face_vertices)? {
                let edge = edge_point(
                    &mut result.edge_points,
                    &mut next_id,
                    sieve,
                    cell_types,
                    a,
                    b,
                )?;
                sieve.add_arrow_o(
                    face,
                    edge,
                    S::Payload::default(),
                    if preserve_orientation {
                        edge_orientation(&canonical_face_vertices, a, b)
                    } else {
                        0
                    },
                )?;
            }
            sieve.add_arrow_o(
                cell,
                face,
                S::Payload::default(),
                if preserve_orientation {
                    face_orientation(&local_face_vertices, &canonical_face_vertices)
                } else {
                    0
                },
            )?;
        }
    }
    Ok(result)
}

fn face_orientation(local: &[PointId], canonical: &[PointId]) -> i32 {
    if same_cyclic_orientation(local, canonical) {
        0
    } else if local.len() == canonical.len() {
        -1
    } else {
        0
    }
}

fn same_cyclic_orientation(a: &[PointId], b: &[PointId]) -> bool {
    if a.len() != b.len() || a.is_empty() {
        return false;
    }
    let Some(start) = a.iter().position(|p| *p == b[0]) else {
        return false;
    };
    (0..a.len()).all(|i| a[(start + i) % a.len()] == b[i])
}

fn canonical_face_loop<St: Storage<f64>>(
    vertices: &[PointId],
    face_type: CellType,
    coordinates: &Coordinates<f64, St>,
) -> Result<Vec<PointId>, MeshSieveError> {
    if !matches!(face_type, CellType::Triangle | CellType::Quadrilateral) || vertices.len() < 3 {
        return Ok(vertices.to_vec());
    }
    let points: Vec<[f64; 3]> = vertices
        .iter()
        .map(|&point| {
            let values = coordinates.try_restrict(point)?;
            if values.len() != 2 && values.len() != 3 {
                return Err(MeshSieveError::InvalidGeometry(
                    "expected 2D/3D coordinates".into(),
                ));
            }
            Ok([values[0], values[1], values.get(2).copied().unwrap_or(0.0)])
        })
        .collect::<Result<_, _>>()?;
    let centroid = points.iter().fold([0.0; 3], |mut sum, point| {
        for axis in 0..3 {
            sum[axis] += point[axis] / points.len() as f64;
        }
        sum
    });
    let mut normal = [0.0; 3];
    'normal: for i in 1..points.len() {
        for j in (i + 1)..points.len() {
            normal = cross3(sub3(points[i], points[0]), sub3(points[j], points[0]));
            if normal.iter().any(|value| value.abs() > 1e-12) {
                break 'normal;
            }
        }
    }
    let axis = (0..3)
        .max_by(|&a, &b| normal[a].abs().total_cmp(&normal[b].abs()))
        .unwrap_or(2);
    let (u, v) = match axis {
        0 => (1, 2),
        1 => (0, 2),
        _ => (0, 1),
    };
    let mut order: Vec<usize> = (0..points.len()).collect();
    order.sort_by(|&a, &b| {
        (points[a][v] - centroid[v])
            .atan2(points[a][u] - centroid[u])
            .total_cmp(&(points[b][v] - centroid[v]).atan2(points[b][u] - centroid[u]))
    });
    let mut best: Option<(Vec<usize>, Vec<f64>)> = None;
    for reversed in [false, true] {
        let mut candidate = order.clone();
        if reversed {
            candidate.reverse();
        }
        for start in 0..candidate.len() {
            candidate.rotate_left(start);
            let key: Vec<f64> = candidate.iter().flat_map(|&index| points[index]).collect();
            if best.as_ref().is_none_or(|(_, old)| lex_less(&key, old)) {
                best = Some((candidate.clone(), key));
            }
        }
    }
    Ok(best
        .map(|(order, _)| order.into_iter().map(|index| vertices[index]).collect())
        .unwrap_or_else(|| vertices.to_vec()))
}

fn sub3(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {
    [a[0] - b[0], a[1] - b[1], a[2] - b[2]]
}

fn cross3(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {
    [
        a[1] * b[2] - a[2] * b[1],
        a[2] * b[0] - a[0] * b[2],
        a[0] * b[1] - a[1] * b[0],
    ]
}

fn lex_less(a: &[f64], b: &[f64]) -> bool {
    a.iter().zip(b).find_map(|(left, right)| {
        let ordering = left.total_cmp(right);
        (ordering != std::cmp::Ordering::Equal).then_some(ordering == std::cmp::Ordering::Less)
    }) == Some(true)
}

fn collect_cells<CtSt>(
    cell_types: &Section<CellType, CtSt>,
) -> Result<Vec<(PointId, CellType)>, MeshSieveError>
where
    CtSt: Storage<CellType>,
{
    let mut cells = Vec::new();
    for (point, cell_slice) in cell_types.iter() {
        if cell_slice.len() != 1 {
            return Err(MeshSieveError::SliceLengthMismatch {
                point,
                expected: 1,
                found: cell_slice.len(),
            });
        }
        let cell_type = cell_slice[0];
        if matches!(
            cell_type,
            CellType::Triangle
                | CellType::Quadrilateral
                | CellType::Tetrahedron
                | CellType::Hexahedron
                | CellType::Prism
                | CellType::Pyramid
                | CellType::Polygon(_)
                | CellType::Polyhedron
                | CellType::Simplex(2)
                | CellType::Simplex(3)
        ) {
            cells.push((point, cell_type));
        } else if cell_type.dimension() >= 2 {
            return Err(MeshSieveError::InvalidGeometry(format!(
                "unsupported cell type: {cell_type:?}"
            )));
        }
    }
    Ok(cells)
}

fn cell_edges(
    cell_type: CellType,
    vertices: &[PointId],
    polyhedron_faces: Option<&[Vec<usize>]>,
) -> Result<Vec<[PointId; 2]>, MeshSieveError> {
    let edges = match cell_type {
        CellType::Triangle | CellType::Simplex(2) => vec![
            [vertices[0], vertices[1]],
            [vertices[1], vertices[2]],
            [vertices[2], vertices[0]],
        ],
        CellType::Quadrilateral => vec![
            [vertices[0], vertices[1]],
            [vertices[1], vertices[2]],
            [vertices[2], vertices[3]],
            [vertices[3], vertices[0]],
        ],
        CellType::Polygon(_) => polygon_edges(vertices),
        CellType::Tetrahedron | CellType::Simplex(3) => vec![
            [vertices[0], vertices[1]],
            [vertices[1], vertices[2]],
            [vertices[2], vertices[0]],
            [vertices[0], vertices[3]],
            [vertices[1], vertices[3]],
            [vertices[2], vertices[3]],
        ],
        CellType::Hexahedron => vec![
            [vertices[0], vertices[1]],
            [vertices[1], vertices[2]],
            [vertices[2], vertices[3]],
            [vertices[3], vertices[0]],
            [vertices[4], vertices[5]],
            [vertices[5], vertices[6]],
            [vertices[6], vertices[7]],
            [vertices[7], vertices[4]],
            [vertices[0], vertices[4]],
            [vertices[1], vertices[5]],
            [vertices[2], vertices[6]],
            [vertices[3], vertices[7]],
        ],
        CellType::Prism => vec![
            [vertices[0], vertices[1]],
            [vertices[1], vertices[2]],
            [vertices[2], vertices[0]],
            [vertices[3], vertices[4]],
            [vertices[4], vertices[5]],
            [vertices[5], vertices[3]],
            [vertices[0], vertices[3]],
            [vertices[1], vertices[4]],
            [vertices[2], vertices[5]],
        ],
        CellType::Pyramid => vec![
            [vertices[0], vertices[1]],
            [vertices[1], vertices[2]],
            [vertices[2], vertices[3]],
            [vertices[3], vertices[0]],
            [vertices[0], vertices[4]],
            [vertices[1], vertices[4]],
            [vertices[2], vertices[4]],
            [vertices[3], vertices[4]],
        ],
        CellType::Polyhedron => polyhedron_edges(vertices, polyhedron_faces)?,
        _ => {
            return Err(MeshSieveError::InvalidGeometry(format!(
                "unsupported cell type: {cell_type:?}"
            )));
        }
    };
    Ok(edges)
}

fn cell_faces(
    cell_type: CellType,
    vertices: &[PointId],
    polyhedron_faces: Option<&[Vec<usize>]>,
) -> Result<Vec<(Vec<PointId>, CellType)>, MeshSieveError> {
    let faces = match cell_type {
        CellType::Triangle
        | CellType::Simplex(2)
        | CellType::Quadrilateral
        | CellType::Polygon(_) => Vec::new(),
        CellType::Tetrahedron | CellType::Simplex(3) => vec![
            (
                vec![vertices[0], vertices[1], vertices[2]],
                CellType::Triangle,
            ),
            (
                vec![vertices[0], vertices[1], vertices[3]],
                CellType::Triangle,
            ),
            (
                vec![vertices[1], vertices[2], vertices[3]],
                CellType::Triangle,
            ),
            (
                vec![vertices[0], vertices[2], vertices[3]],
                CellType::Triangle,
            ),
        ],
        CellType::Hexahedron => vec![
            (
                vec![vertices[0], vertices[1], vertices[2], vertices[3]],
                CellType::Quadrilateral,
            ),
            (
                vec![vertices[4], vertices[5], vertices[6], vertices[7]],
                CellType::Quadrilateral,
            ),
            (
                vec![vertices[0], vertices[1], vertices[5], vertices[4]],
                CellType::Quadrilateral,
            ),
            (
                vec![vertices[1], vertices[2], vertices[6], vertices[5]],
                CellType::Quadrilateral,
            ),
            (
                vec![vertices[2], vertices[3], vertices[7], vertices[6]],
                CellType::Quadrilateral,
            ),
            (
                vec![vertices[3], vertices[0], vertices[4], vertices[7]],
                CellType::Quadrilateral,
            ),
        ],
        CellType::Prism => vec![
            (
                vec![vertices[0], vertices[1], vertices[2]],
                CellType::Triangle,
            ),
            (
                vec![vertices[3], vertices[4], vertices[5]],
                CellType::Triangle,
            ),
            (
                vec![vertices[0], vertices[1], vertices[4], vertices[3]],
                CellType::Quadrilateral,
            ),
            (
                vec![vertices[1], vertices[2], vertices[5], vertices[4]],
                CellType::Quadrilateral,
            ),
            (
                vec![vertices[2], vertices[0], vertices[3], vertices[5]],
                CellType::Quadrilateral,
            ),
        ],
        CellType::Pyramid => vec![
            (
                vec![vertices[0], vertices[1], vertices[2], vertices[3]],
                CellType::Quadrilateral,
            ),
            (
                vec![vertices[0], vertices[1], vertices[4]],
                CellType::Triangle,
            ),
            (
                vec![vertices[1], vertices[2], vertices[4]],
                CellType::Triangle,
            ),
            (
                vec![vertices[2], vertices[3], vertices[4]],
                CellType::Triangle,
            ),
            (
                vec![vertices[3], vertices[0], vertices[4]],
                CellType::Triangle,
            ),
        ],
        CellType::Polyhedron => polyhedron_faces_to_vertices(vertices, polyhedron_faces)?,
        _ => {
            return Err(MeshSieveError::InvalidGeometry(format!(
                "unsupported cell type: {cell_type:?}"
            )));
        }
    };
    Ok(faces)
}

fn face_edges(vertices: &[PointId]) -> Result<Vec<[PointId; 2]>, MeshSieveError> {
    let n = vertices.len();
    if n < 3 {
        return Err(MeshSieveError::InvalidGeometry(format!(
            "unsupported face vertex count: {n}"
        )));
    }
    let mut edges = Vec::with_capacity(n);
    for i in 0..n {
        edges.push([vertices[i], vertices[(i + 1) % n]]);
    }
    Ok(edges)
}

fn edge_point<S, CtSt>(
    edges: &mut BTreeMap<(PointId, PointId), PointId>,
    next_id: &mut u64,
    sieve: &mut S,
    cell_types: &mut Section<CellType, CtSt>,
    a: PointId,
    b: PointId,
) -> Result<PointId, MeshSieveError>
where
    S: MutableSieve<Point = PointId> + OrientedSieve<Orient = i32>,
    S::Payload: Default,
    CtSt: Storage<CellType> + Clone,
{
    let key = if a < b { (a, b) } else { (b, a) };
    if let Some(p) = edges.get(&key) {
        return Ok(*p);
    }
    let edge = alloc_point(next_id)?;
    edges.insert(key, edge);
    MutableSieve::add_point(sieve, edge);
    sieve.add_arrow_o(edge, a, S::Payload::default(), 0)?;
    sieve.add_arrow_o(edge, b, S::Payload::default(), 0)?;
    cell_types.try_add_point(edge, 1)?;
    cell_types.try_set(edge, &[CellType::Segment])?;
    Ok(edge)
}

fn face_point<S, CtSt>(
    faces: &mut BTreeMap<Vec<PointId>, PointId>,
    directions: &mut BTreeMap<Vec<PointId>, Vec<PointId>>,
    next_id: &mut u64,
    sieve: &mut S,
    cell_types: &mut Section<CellType, CtSt>,
    vertices: &[PointId],
    face_type: CellType,
) -> Result<(PointId, Vec<PointId>), MeshSieveError>
where
    S: MutableSieve<Point = PointId> + OrientedSieve<Orient = i32>,
    S::Payload: Default,
    CtSt: Storage<CellType> + Clone,
{
    let mut key = vertices.to_vec();
    key.sort();
    if let Some(p) = faces.get(&key) {
        let direction = directions
            .get(&key)
            .cloned()
            .unwrap_or_else(|| vertices.to_vec());
        return Ok((*p, direction));
    }
    let face = alloc_point(next_id)?;
    directions.insert(key.clone(), vertices.to_vec());
    faces.insert(key, face);
    MutableSieve::add_point(sieve, face);
    cell_types.try_add_point(face, 1)?;
    cell_types.try_set(face, &[face_type])?;
    Ok((face, vertices.to_vec()))
}

fn alloc_point(next_id: &mut u64) -> Result<PointId, MeshSieveError> {
    let id = PointId::new(*next_id)?;
    *next_id = next_id
        .checked_add(1)
        .ok_or(MeshSieveError::InvalidPointId)?;
    Ok(id)
}

fn apply_permutation(
    cell: PointId,
    vertices: &[PointId],
    permutation: &[usize],
) -> Result<Vec<PointId>, MeshSieveError> {
    if permutation.len() != vertices.len() {
        return Err(MeshSieveError::InvalidGeometry(format!(
            "cell {cell:?} permutation length {} does not match vertex count {}",
            permutation.len(),
            vertices.len()
        )));
    }
    let mut seen = vec![false; vertices.len()];
    let mut ordered = Vec::with_capacity(vertices.len());
    for &idx in permutation {
        if idx >= vertices.len() {
            return Err(MeshSieveError::InvalidGeometry(format!(
                "cell {cell:?} permutation index {idx} out of bounds"
            )));
        }
        if seen[idx] {
            return Err(MeshSieveError::InvalidGeometry(format!(
                "cell {cell:?} permutation index {idx} appears multiple times"
            )));
        }
        seen[idx] = true;
        ordered.push(vertices[idx]);
    }
    Ok(ordered)
}

fn polygon_edges(vertices: &[PointId]) -> Vec<[PointId; 2]> {
    let n = vertices.len();
    let mut edges = Vec::with_capacity(n);
    for i in 0..n {
        edges.push([vertices[i], vertices[(i + 1) % n]]);
    }
    edges
}

fn polyhedron_edges(
    vertices: &[PointId],
    polyhedron_faces: Option<&[Vec<usize>]>,
) -> Result<Vec<[PointId; 2]>, MeshSieveError> {
    let faces = polyhedron_faces.ok_or_else(|| {
        MeshSieveError::InvalidGeometry("polyhedron faces metadata required".to_string())
    })?;
    let mut edges = BTreeMap::new();
    for face in faces {
        let face_vertices = face_indices_to_vertices(vertices, face)?;
        for [a, b] in face_edges(&face_vertices)? {
            let key = if a < b { (a, b) } else { (b, a) };
            edges.entry(key).or_insert([a, b]);
        }
    }
    Ok(edges.values().copied().collect())
}

fn polyhedron_faces_to_vertices(
    vertices: &[PointId],
    polyhedron_faces: Option<&[Vec<usize>]>,
) -> Result<Vec<(Vec<PointId>, CellType)>, MeshSieveError> {
    let faces = polyhedron_faces.ok_or_else(|| {
        MeshSieveError::InvalidGeometry("polyhedron faces metadata required".to_string())
    })?;
    let mut result = Vec::with_capacity(faces.len());
    for face in faces {
        let face_vertices = face_indices_to_vertices(vertices, face)?;
        let face_type = face_type_for_vertices(face_vertices.len())?;
        result.push((face_vertices, face_type));
    }
    Ok(result)
}

fn face_indices_to_vertices(
    vertices: &[PointId],
    face: &[usize],
) -> Result<Vec<PointId>, MeshSieveError> {
    if face.len() < 3 {
        return Err(MeshSieveError::InvalidGeometry(format!(
            "unsupported face vertex count: {}",
            face.len()
        )));
    }
    let mut face_vertices = Vec::with_capacity(face.len());
    for &idx in face {
        let v = vertices.get(idx).ok_or_else(|| {
            MeshSieveError::InvalidGeometry(format!(
                "face index {idx} out of bounds for {} vertices",
                vertices.len()
            ))
        })?;
        face_vertices.push(*v);
    }
    Ok(face_vertices)
}

fn face_type_for_vertices(count: usize) -> Result<CellType, MeshSieveError> {
    match count {
        3 => Ok(CellType::Triangle),
        4 => Ok(CellType::Quadrilateral),
        n => {
            let n_u8 = u8::try_from(n).map_err(|_| {
                MeshSieveError::InvalidGeometry(format!("polygon face vertex count {n} exceeds u8"))
            })?;
            Ok(CellType::Polygon(n_u8))
        }
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::data::atlas::Atlas;
    use crate::data::storage::VecStorage;
    use crate::topology::sieve::{InMemorySieve, Sieve};

    fn v(raw: u64) -> PointId {
        PointId::new(raw).unwrap()
    }

    #[test]
    fn interpolates_quad_edges() {
        let mut sieve = InMemorySieve::<PointId, ()>::default();
        let cell = v(1);
        let vertices = [v(2), v(3), v(4), v(5)];
        MutableSieve::add_point(&mut sieve, cell);
        for p in vertices {
            MutableSieve::add_point(&mut sieve, p);
            sieve.add_arrow(cell, p, ()).unwrap();
        }

        let mut atlas = Atlas::default();
        atlas.try_insert(cell, 1).unwrap();
        for p in vertices {
            atlas.try_insert(p, 1).unwrap();
        }
        let mut cell_types = Section::<CellType, VecStorage<CellType>>::new(atlas);
        cell_types
            .try_set(cell, &[CellType::Quadrilateral])
            .unwrap();
        for p in vertices {
            cell_types.try_set(p, &[CellType::Vertex]).unwrap();
        }

        let result = interpolate_edges_faces(&mut sieve, &mut cell_types).unwrap();
        assert_eq!(result.edge_points.len(), 4);
        assert!(result.face_points.is_empty());
        for ((a, b), edge) in &result.edge_points {
            assert!(sieve.has_arrow(*edge, *a));
            assert!(sieve.has_arrow(*edge, *b));
            assert!(sieve.has_arrow(cell, *edge));
        }
    }

    #[test]
    fn interpolates_mixed_elements_with_ordering() {
        let mut sieve = InMemorySieve::<PointId, ()>::default();
        let cell_quad = v(1);
        let quad_vertices = [v(2), v(3), v(4), v(5)];
        let cell_polygon = v(10);
        let polygon_vertices = [v(11), v(12), v(13), v(14), v(15)];
        let cell_polyhedron = v(20);
        let poly_vertices = [v(21), v(22), v(23), v(24)];

        let mut all_points = vec![cell_quad, cell_polygon, cell_polyhedron];
        all_points.extend(quad_vertices);
        all_points.extend(polygon_vertices);
        all_points.extend(poly_vertices);
        for p in &all_points {
            MutableSieve::add_point(&mut sieve, *p);
        }
        for p in quad_vertices {
            sieve.add_arrow(cell_quad, p, ()).unwrap();
        }
        for p in polygon_vertices {
            sieve.add_arrow(cell_polygon, p, ()).unwrap();
        }
        for p in poly_vertices {
            sieve.add_arrow(cell_polyhedron, p, ()).unwrap();
        }

        let mut atlas = Atlas::default();
        for p in &all_points {
            atlas.try_insert(*p, 1).unwrap();
        }
        let mut cell_types = Section::<CellType, VecStorage<CellType>>::new(atlas);
        cell_types
            .try_set(cell_quad, &[CellType::Quadrilateral])
            .unwrap();
        cell_types
            .try_set(cell_polygon, &[CellType::Polygon(5)])
            .unwrap();
        cell_types
            .try_set(cell_polyhedron, &[CellType::Polyhedron])
            .unwrap();
        for p in quad_vertices
            .into_iter()
            .chain(polygon_vertices)
            .chain(poly_vertices)
        {
            cell_types.try_set(p, &[CellType::Vertex]).unwrap();
        }

        let mut ordering = InterpolationOrdering::default();
        ordering
            .vertex_permutations
            .insert(cell_polygon, vec![0, 2, 3, 4, 1]);
        ordering.polyhedron_faces.insert(
            cell_polyhedron,
            vec![vec![0, 1, 2], vec![0, 1, 3], vec![1, 2, 3], vec![0, 2, 3]],
        );

        let result =
            interpolate_edges_faces_with_ordering(&mut sieve, &mut cell_types, Some(&ordering))
                .unwrap();

        assert_eq!(result.edge_points.len(), 15);
        assert_eq!(result.face_points.len(), 4);

        let key = if polygon_vertices[0] < polygon_vertices[2] {
            (polygon_vertices[0], polygon_vertices[2])
        } else {
            (polygon_vertices[2], polygon_vertices[0])
        };
        assert!(result.edge_points.contains_key(&key));

        for face in result.face_points.values() {
            assert!(sieve.has_arrow(cell_polyhedron, *face));
            let face_type = cell_types.try_restrict(*face).unwrap()[0];
            assert!(matches!(face_type, CellType::Triangle));
            assert!(sieve.cone_points(*face).next().is_some());
        }
    }
}