mesh-sieve 4.0.1

//! Reference refinement templates and mesh refinement helpers.
//!
//! # Expected invariants
//! - The coarse topology is a cell complex where each cell points directly to its vertices
//!   or to intermediate edges/faces that can be reduced to vertices.
//! - `cell_types` contains **exactly one** [`CellType`] entry for every cell to refine.
//! - Each supported cell type has the expected number of vertices in its closure.
//! - Vertex ordering in each cone is **consistent** (e.g., counter-clockwise) so the
//!   reference templates preserve orientation. When intermediate edges are present,
//!   we try to recover this ordering from the edge cycle; otherwise vertices are
//!   canonicalized by PointId order.
//! - Shared edges/faces must appear with the same vertex sets so mid-edge and face-center
//!   points can be reused.
//!
//! The refinement map returned by [`refine_mesh`] is suitable for
//! [`crate::data::refine::SievedArray::try_refine_with_sifter`], enabling transfer of
//! cell-wise data from the coarse mesh to the refined mesh.

use crate::algs::interpolate::{InterpolationResult, interpolate_edges_faces};
use crate::data::atlas::Atlas;
use crate::data::coordinates::{Coordinates, HighOrderCoordinates};
use crate::data::hanging_node_constraints::{
    HangingNodeConstraints, constraints_from_topological_anchors,
};
use crate::data::section::Section;
use crate::data::storage::{Storage, VecStorage};
use crate::geometry::quality::validate_cell_geometry;
use crate::mesh_error::MeshSieveError;
use crate::topology::anchors::{AnchorKind, TopologicalAnchors};
use crate::topology::arrow::Polarity;
use crate::topology::cell_type::CellType;
use crate::topology::ownership::PointOwnership;
use crate::topology::point::PointId;
use crate::topology::sieve::{MeshSieve, OrientedSieve, Sieve};
use crate::topology::validation::debug_validate_overlap_ownership_topology;
use std::collections::{HashMap, HashSet};

/// Mapping from coarse cells to refined cells, with orientation hints.
pub type RefinementMap = Vec<(PointId, Vec<(PointId, Polarity)>)>;

/// Output of [`refine_mesh`], including the refined topology and transfer map.
#[derive(Clone, Debug)]
pub struct RefinedMesh {
    /// Refined topology (cells → vertices).
    pub sieve: MeshSieve,
    /// Mapping from coarse cell to refined cells (all `Polarity::Forward`).
    pub cell_refinement: RefinementMap,
    /// Optional coordinates for the refined mesh vertices.
    pub coordinates: Option<Coordinates<f64, VecStorage<f64>>>,
    /// Parent/constrained-point anchors produced by refinement.
    pub anchors: TopologicalAnchors,
    /// Automatically generated scalar hanging-node constraints from anchors.
    pub hanging_constraints: HangingNodeConstraints<f64>,
}

/// Output of [`refine_mesh_with_ownership`], including ownership metadata updates.
#[derive(Clone, Debug)]
pub struct RefinedMeshWithOwnership {
    /// Refined topology (cells → vertices).
    pub sieve: MeshSieve,
    /// Mapping from coarse cell to refined cells (all `Polarity::Forward`).
    pub cell_refinement: RefinementMap,
    /// Updated ownership metadata including new points.
    pub ownership: PointOwnership,
    /// Optional coordinates for the refined mesh vertices.
    pub coordinates: Option<Coordinates<f64, VecStorage<f64>>>,
    /// Parent/constrained-point anchors produced by refinement.
    pub anchors: TopologicalAnchors,
    /// Automatically generated scalar hanging-node constraints from anchors.
    pub hanging_constraints: HangingNodeConstraints<f64>,
}

/// Output of [`refine_mesh_full_topology`], including a full topology and types.
#[derive(Clone, Debug)]
pub struct RefinedMeshWithTopology {
    /// Refined topology (cells → faces → edges → vertices).
    pub sieve: MeshSieve,
    /// Mapping from coarse cell to refined cells (all `Polarity::Forward`).
    pub cell_refinement: RefinementMap,
    /// Cell types for the refined topology (cells, edges, faces, vertices).
    pub cell_types: Section<CellType, VecStorage<CellType>>,
    /// Optional coordinates for the refined mesh vertices.
    pub coordinates: Option<Coordinates<f64, VecStorage<f64>>>,
    /// Parent/constrained-point anchors produced by refinement.
    pub anchors: TopologicalAnchors,
    /// Automatically generated scalar hanging-node constraints from anchors.
    pub hanging_constraints: HangingNodeConstraints<f64>,
}

/// Optional anisotropic split hints for refinement.
#[derive(Clone, Debug, Default)]
pub struct AnisotropicSplitHints {
    /// Per-cell edge splits represented by vertex pairs.
    pub edge_splits: HashMap<PointId, Vec<[PointId; 2]>>,
    /// Per-cell face splits represented by ordered vertex loops.
    pub face_splits: HashMap<PointId, Vec<Vec<PointId>>>,
}

impl AnisotropicSplitHints {
    /// Returns true if there are no split hints recorded.
    pub fn is_empty(&self) -> bool {
        self.edge_splits.is_empty() && self.face_splits.is_empty()
    }
}

/// Optional settings for refinement.
#[derive(Clone, Debug, Default)]
pub struct RefineOptions {
    /// When enabled, validate cell geometry and reject inverted or degenerate elements.
    pub check_geometry: bool,
    /// Optional anisotropic split hints for edges/faces.
    pub anisotropic_splits: Option<AnisotropicSplitHints>,
}

/// Pre-interpolate a cell→vertex topology to a full cell→face→edge→vertex mesh.
pub fn pre_interpolate_topology<S, CtSt>(
    sieve: &mut S,
    cell_types: &mut Section<CellType, CtSt>,
) -> Result<InterpolationResult, MeshSieveError>
where
    S: crate::topology::sieve::MutableSieve<Point = PointId> + OrientedSieve<Orient = i32>,
    S::Payload: Default,
    CtSt: Storage<CellType> + Clone,
{
    interpolate_edges_faces(sieve, cell_types)
}

/// Collapse a full topology to cell→vertex connectivity for refinement.
pub fn collapse_to_cell_vertices<CtSt>(
    sieve: &mut impl Sieve<Point = PointId>,
    cell_types: &Section<CellType, CtSt>,
) -> Result<MeshSieve, MeshSieveError>
where
    CtSt: Storage<CellType>,
{
    let mesh_dimension = crate::topology::utils::dimension(sieve)?;
    let cells = collect_top_cells(cell_types, Some(mesh_dimension))?;
    let mut collapsed = MeshSieve::default();

    for (cell, cell_type) in cells {
        let vertices = match cell_type {
            CellType::Polygon(n) => cell_vertices(sieve, cell, n as usize)?,
            CellType::Polyhedron => cell_vertices_from_closure(sieve, cell)?,
            _ => {
                let expected = expected_vertex_count(cell_type)
                    .ok_or(MeshSieveError::UnsupportedRefinementCellType { cell, cell_type })?;
                cell_vertices(sieve, cell, expected)?
            }
        };
        for v in vertices {
            collapsed.add_arrow(cell, v, ());
        }
    }

    Ok(collapsed)
}

/// Reference 1→4 subdivision of a triangle given vertices and mid-edge points.
pub fn triangle_subdivision(vertices: [PointId; 3], midpoints: [PointId; 3]) -> [[PointId; 3]; 4] {
    let [v0, v1, v2] = vertices;
    let [m01, m12, m20] = midpoints;
    [
        [v0, m01, m20],
        [v1, m12, m01],
        [v2, m20, m12],
        [m01, m12, m20],
    ]
}

/// Reference 1→4 subdivision of a quad given vertices, mid-edge points, and a center point.
pub fn quadrilateral_subdivision(
    vertices: [PointId; 4],
    midpoints: [PointId; 4],
    center: PointId,
) -> [[PointId; 4]; 4] {
    let [v0, v1, v2, v3] = vertices;
    let [m01, m12, m23, m30] = midpoints;
    [
        [v0, m01, center, m30],
        [m01, v1, m12, center],
        [center, m12, v2, m23],
        [m30, center, m23, v3],
    ]
}

/// Reference 1→8 subdivision of a tetrahedron given vertices and mid-edge points.
pub fn tetrahedron_subdivision(
    vertices: [PointId; 4],
    midpoints: [PointId; 6],
) -> [[PointId; 4]; 8] {
    let [v0, v1, v2, v3] = vertices;
    let [m01, m12, m20, m03, m13, m23] = midpoints;
    [
        [v0, m01, m20, m03],
        [v1, m12, m01, m13],
        [v2, m20, m12, m23],
        [v3, m03, m13, m23],
        [m01, m12, m13, m23],
        [m01, m13, m03, m23],
        [m01, m03, m20, m23],
        [m01, m20, m12, m23],
    ]
}

/// Reference 1→8 subdivision of a hex using edge midpoints, face centers, and cell center.
pub fn hexahedron_subdivision(
    vertices: [PointId; 8],
    edge_midpoints: [PointId; 12],
    face_centers: [PointId; 6],
    center: PointId,
) -> [[PointId; 8]; 8] {
    let [v0, v1, v2, v3, v4, v5, v6, v7] = vertices;
    let [m01, m12, m23, m30, m45, m56, m67, m74, m04, m15, m26, m37] = edge_midpoints;
    let [f0123, f4567, f0154, f1265, f2376, f3047] = face_centers;

    let mut grid: HashMap<(u8, u8, u8), PointId> = HashMap::new();
    let mut insert = |x: u8, y: u8, z: u8, p: PointId| {
        grid.insert((x, y, z), p);
    };

    insert(0, 0, 0, v0);
    insert(2, 0, 0, v1);
    insert(2, 2, 0, v2);
    insert(0, 2, 0, v3);
    insert(0, 0, 2, v4);
    insert(2, 0, 2, v5);
    insert(2, 2, 2, v6);
    insert(0, 2, 2, v7);

    insert(1, 0, 0, m01);
    insert(2, 1, 0, m12);
    insert(1, 2, 0, m23);
    insert(0, 1, 0, m30);
    insert(1, 0, 2, m45);
    insert(2, 1, 2, m56);
    insert(1, 2, 2, m67);
    insert(0, 1, 2, m74);
    insert(0, 0, 1, m04);
    insert(2, 0, 1, m15);
    insert(2, 2, 1, m26);
    insert(0, 2, 1, m37);

    insert(1, 1, 0, f0123);
    insert(1, 1, 2, f4567);
    insert(1, 0, 1, f0154);
    insert(2, 1, 1, f1265);
    insert(1, 2, 1, f2376);
    insert(0, 1, 1, f3047);
    insert(1, 1, 1, center);

    let mut sub_hexes = Vec::with_capacity(8);
    for (ix, iy, iz) in [
        (0u8, 0u8, 0u8),
        (1, 0, 0),
        (1, 1, 0),
        (0, 1, 0),
        (0, 0, 1),
        (1, 0, 1),
        (1, 1, 1),
        (0, 1, 1),
    ] {
        let (x0, x1) = if ix == 0 { (0, 1) } else { (1, 2) };
        let (y0, y1) = if iy == 0 { (0, 1) } else { (1, 2) };
        let (z0, z1) = if iz == 0 { (0, 1) } else { (1, 2) };
        let verts = [
            grid[&(x0, y0, z0)],
            grid[&(x1, y0, z0)],
            grid[&(x1, y1, z0)],
            grid[&(x0, y1, z0)],
            grid[&(x0, y0, z1)],
            grid[&(x1, y0, z1)],
            grid[&(x1, y1, z1)],
            grid[&(x0, y1, z1)],
        ];
        sub_hexes.push(verts);
    }
    sub_hexes
        .try_into()
        .expect("hexahedron subdivision should yield 8 sub-hexes")
}

/// Reference 1→4 subdivision of a prism via triangle subdivisions on bottom/top.
pub fn prism_subdivision(
    vertices: [PointId; 6],
    bottom_midpoints: [PointId; 3],
    top_midpoints: [PointId; 3],
) -> [[PointId; 6]; 4] {
    let [v0, v1, v2, v3, v4, v5] = vertices;
    let bottom = triangle_subdivision([v0, v1, v2], bottom_midpoints);
    let top = triangle_subdivision([v3, v4, v5], top_midpoints);
    let mut prisms = Vec::with_capacity(4);
    for i in 0..4 {
        prisms.push([
            bottom[i][0],
            bottom[i][1],
            bottom[i][2],
            top[i][0],
            top[i][1],
            top[i][2],
        ]);
    }
    prisms
        .try_into()
        .expect("prism subdivision should yield 4 sub-prisms")
}

/// Reference 1→4 subdivision of a pyramid via quad subdivision of the base.
pub fn pyramid_subdivision(
    vertices: [PointId; 5],
    base_midpoints: [PointId; 4],
    base_center: PointId,
) -> [[PointId; 5]; 4] {
    let [v0, v1, v2, v3, apex] = vertices;
    let base = quadrilateral_subdivision([v0, v1, v2, v3], base_midpoints, base_center);
    let mut pyramids = Vec::with_capacity(4);
    for quad in base {
        pyramids.push([quad[0], quad[1], quad[2], quad[3], apex]);
    }
    pyramids
        .try_into()
        .expect("pyramid subdivision should yield 4 sub-pyramids")
}

/// Refine a 2D/3D mesh using reference subdivision templates.
///
/// Returns a refined [`InMemorySieve`] plus a refinement map for carrying
/// cell-wise data via [`crate::data::refine::SievedArray::try_refine_with_sifter`].
///
/// # Errors
/// - [`MeshSieveError::UnsupportedRefinementCellType`] for unsupported cell types.
/// - [`MeshSieveError::RefinementTopologyMismatch`] if vertex counts do not match templates.
/// - [`MeshSieveError::SliceLengthMismatch`] if the cell-type section is not length-1.
pub fn refine_mesh<S>(
    coarse: &mut impl Sieve<Point = PointId>,
    cell_types: &Section<CellType, S>,
) -> Result<RefinedMesh, MeshSieveError>
where
    S: Storage<CellType>,
{
    refine_mesh_with_options::<S, crate::data::storage::VecStorage<f64>>(
        coarse,
        cell_types,
        None,
        RefineOptions::default(),
    )
}

/// Refine a 2D/3D mesh using reference subdivision templates with optional checks.
///
/// Returns a refined [`InMemorySieve`] plus a refinement map for carrying
/// cell-wise data via [`crate::data::refine::SievedArray::try_refine_with_sifter`].
///
/// # Errors
/// - [`MeshSieveError::InvalidGeometry`] when geometry checks are requested and fail.
/// - [`MeshSieveError::UnsupportedRefinementCellType`] for unsupported cell types.
/// - [`MeshSieveError::RefinementTopologyMismatch`] if vertex counts do not match templates.
/// - [`MeshSieveError::SliceLengthMismatch`] if the cell-type section is not length-1.
pub fn refine_mesh_with_options<S, Cs>(
    coarse: &mut impl Sieve<Point = PointId>,
    cell_types: &Section<CellType, S>,
    coordinates: Option<&Coordinates<f64, Cs>>,
    options: RefineOptions,
) -> Result<RefinedMesh, MeshSieveError>
where
    S: Storage<CellType>,
    Cs: Storage<f64>,
{
    if let Some(hints) = &options.anisotropic_splits
        && !hints.is_empty()
    {
        let cell_ids: HashSet<_> = cell_types.iter().map(|(cell, _)| cell).collect();
        for cell in hints.edge_splits.keys().chain(hints.face_splits.keys()) {
            if !cell_ids.contains(cell) {
                return Err(MeshSieveError::UnknownPoint(format!(
                    "anisotropic split hint for missing cell {cell:?}"
                )));
            }
        }
    }
    let mut dimension: Option<u8> = None;
    for (cell, cell_slice) in cell_types.iter() {
        if cell_slice.len() != 1 {
            return Err(MeshSieveError::SliceLengthMismatch {
                point: cell,
                expected: 1,
                found: cell_slice.len(),
            });
        }
        let cell_dim = cell_slice[0].dimension();
        if let Some(existing) = dimension {
            if existing != cell_dim {
                return Err(MeshSieveError::InvalidGeometry(format!(
                    "refinement requires a uniform cell dimension, found {existing}D and {cell_dim}D"
                )));
            }
        } else {
            dimension = Some(cell_dim);
        }
    }
    if options.check_geometry {
        let coords = coordinates.ok_or_else(|| {
            MeshSieveError::InvalidGeometry("geometry checks requested without coordinates".into())
        })?;
        for (cell, cell_slice) in cell_types.iter() {
            if cell_slice.len() != 1 {
                return Err(MeshSieveError::SliceLengthMismatch {
                    point: cell,
                    expected: 1,
                    found: cell_slice.len(),
                });
            }
            let cell_type = cell_slice[0];
            let vertices = match cell_type {
                CellType::Polygon(_) | CellType::Polyhedron => {
                    return Err(MeshSieveError::InvalidGeometry(format!(
                        "geometry checks are not supported for {cell_type:?} cells"
                    )));
                }
                _ => {
                    let expected = expected_vertex_count(cell_type)
                        .ok_or(MeshSieveError::UnsupportedRefinementCellType { cell, cell_type })?;
                    cell_vertices(coarse, cell, expected)?
                }
            };
            if let Err(err) = validate_cell_geometry(cell_type, &vertices, coords) {
                return Err(MeshSieveError::InvalidGeometry(format!(
                    "cell {cell:?}: {err}"
                )));
            }
        }
    }

    let mut max_id = 0u64;
    for p in coarse.chart_points()? {
        max_id = max_id.max(p.get());
    }
    let mut next_id = max_id
        .checked_add(1)
        .ok_or(MeshSieveError::InvalidPointId)?;

    let mut refined = MeshSieve::default();
    let mut refinement_map: RefinementMap = Vec::new();
    let mut edge_midpoints: HashMap<(PointId, PointId), PointId> = HashMap::new();
    let mut face_centers: HashMap<Vec<PointId>, PointId> = HashMap::new();
    let mut point_sources: HashMap<PointId, Vec<PointId>> = HashMap::new();
    let hinted_edge_keys = collect_hinted_edge_keys(options.anisotropic_splits.as_ref());

    for (cell, cell_slice) in cell_types.iter() {
        if cell_slice.len() != 1 {
            return Err(MeshSieveError::SliceLengthMismatch {
                point: cell,
                expected: 1,
                found: cell_slice.len(),
            });
        }
        let cell_type = cell_slice[0];
        let vertices = match cell_type {
            CellType::Polygon(n) => cell_vertices(coarse, cell, n as usize)?,
            CellType::Polyhedron => cell_vertices_from_closure(coarse, cell)?,
            _ => {
                let expected = expected_vertex_count(cell_type)
                    .ok_or(MeshSieveError::UnsupportedRefinementCellType { cell, cell_type })?;
                cell_vertices(coarse, cell, expected)?
            }
        };
        let mut fine_cells = Vec::new();

        match cell_type {
            CellType::Triangle => {
                let vertices = [vertices[0], vertices[1], vertices[2]];
                let midpoints = [
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[0],
                        vertices[1],
                    )?,
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[1],
                        vertices[2],
                    )?,
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[2],
                        vertices[0],
                    )?,
                ];
                let split_indices =
                    triangle_split_indices(cell, vertices, options.anisotropic_splits.as_ref());
                let templates =
                    anisotropic_triangle_subdivision(vertices, midpoints, &split_indices);
                for verts in templates {
                    let fine_cell = alloc_point(&mut next_id)?;
                    for v in verts {
                        refined.add_arrow(fine_cell, v, ());
                    }
                    fine_cells.push((fine_cell, Polarity::Forward));
                }
            }
            CellType::Quadrilateral => {
                let vertices = [vertices[0], vertices[1], vertices[2], vertices[3]];
                let midpoints = [
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[0],
                        vertices[1],
                    )?,
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[1],
                        vertices[2],
                    )?,
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[2],
                        vertices[3],
                    )?,
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[3],
                        vertices[0],
                    )?,
                ];
                let center = alloc_point(&mut next_id)?;
                point_sources.insert(center, vertices.to_vec());
                let split_indices = quadrilateral_split_indices(
                    cell,
                    vertices,
                    options.anisotropic_splits.as_ref(),
                );
                let templates = anisotropic_quadrilateral_subdivision(
                    vertices,
                    midpoints,
                    center,
                    &split_indices,
                );
                for verts in templates {
                    let fine_cell = alloc_point(&mut next_id)?;
                    for v in verts {
                        refined.add_arrow(fine_cell, v, ());
                    }
                    fine_cells.push((fine_cell, Polarity::Forward));
                }
            }
            CellType::Tetrahedron => {
                let vertices = [vertices[0], vertices[1], vertices[2], vertices[3]];
                let midpoints = [
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[0],
                        vertices[1],
                    )?,
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[1],
                        vertices[2],
                    )?,
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[2],
                        vertices[0],
                    )?,
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[0],
                        vertices[3],
                    )?,
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[1],
                        vertices[3],
                    )?,
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[2],
                        vertices[3],
                    )?,
                ];
                for verts in tetrahedron_subdivision(vertices, midpoints) {
                    let fine_cell = alloc_point(&mut next_id)?;
                    for v in verts {
                        refined.add_arrow(fine_cell, v, ());
                    }
                    fine_cells.push((fine_cell, Polarity::Forward));
                }
            }
            CellType::Hexahedron => {
                let vertices = [
                    vertices[0],
                    vertices[1],
                    vertices[2],
                    vertices[3],
                    vertices[4],
                    vertices[5],
                    vertices[6],
                    vertices[7],
                ];
                let edge_midpoints = [
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[0],
                        vertices[1],
                    )?,
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[1],
                        vertices[2],
                    )?,
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[2],
                        vertices[3],
                    )?,
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[3],
                        vertices[0],
                    )?,
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[4],
                        vertices[5],
                    )?,
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[5],
                        vertices[6],
                    )?,
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[6],
                        vertices[7],
                    )?,
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[7],
                        vertices[4],
                    )?,
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[0],
                        vertices[4],
                    )?,
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[1],
                        vertices[5],
                    )?,
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[2],
                        vertices[6],
                    )?,
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[3],
                        vertices[7],
                    )?,
                ];
                let face_centers = [
                    face_center(
                        &mut face_centers,
                        &mut point_sources,
                        &mut next_id,
                        &[vertices[0], vertices[1], vertices[2], vertices[3]],
                    )?,
                    face_center(
                        &mut face_centers,
                        &mut point_sources,
                        &mut next_id,
                        &[vertices[4], vertices[5], vertices[6], vertices[7]],
                    )?,
                    face_center(
                        &mut face_centers,
                        &mut point_sources,
                        &mut next_id,
                        &[vertices[0], vertices[1], vertices[5], vertices[4]],
                    )?,
                    face_center(
                        &mut face_centers,
                        &mut point_sources,
                        &mut next_id,
                        &[vertices[1], vertices[2], vertices[6], vertices[5]],
                    )?,
                    face_center(
                        &mut face_centers,
                        &mut point_sources,
                        &mut next_id,
                        &[vertices[2], vertices[3], vertices[7], vertices[6]],
                    )?,
                    face_center(
                        &mut face_centers,
                        &mut point_sources,
                        &mut next_id,
                        &[vertices[3], vertices[0], vertices[4], vertices[7]],
                    )?,
                ];
                let center = alloc_point(&mut next_id)?;
                point_sources.insert(center, vertices.to_vec());
                for verts in hexahedron_subdivision(vertices, edge_midpoints, face_centers, center)
                {
                    let fine_cell = alloc_point(&mut next_id)?;
                    for v in verts {
                        refined.add_arrow(fine_cell, v, ());
                    }
                    fine_cells.push((fine_cell, Polarity::Forward));
                }
            }
            CellType::Prism => {
                let vertices = [
                    vertices[0],
                    vertices[1],
                    vertices[2],
                    vertices[3],
                    vertices[4],
                    vertices[5],
                ];
                let bottom_midpoints = [
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[0],
                        vertices[1],
                    )?,
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[1],
                        vertices[2],
                    )?,
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[2],
                        vertices[0],
                    )?,
                ];
                let top_midpoints = [
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[3],
                        vertices[4],
                    )?,
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[4],
                        vertices[5],
                    )?,
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[5],
                        vertices[3],
                    )?,
                ];
                for verts in prism_subdivision(vertices, bottom_midpoints, top_midpoints) {
                    let fine_cell = alloc_point(&mut next_id)?;
                    for v in verts {
                        refined.add_arrow(fine_cell, v, ());
                    }
                    fine_cells.push((fine_cell, Polarity::Forward));
                }
            }
            CellType::Pyramid => {
                let vertices = [
                    vertices[0],
                    vertices[1],
                    vertices[2],
                    vertices[3],
                    vertices[4],
                ];
                let base_midpoints = [
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[0],
                        vertices[1],
                    )?,
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[1],
                        vertices[2],
                    )?,
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[2],
                        vertices[3],
                    )?,
                    midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[3],
                        vertices[0],
                    )?,
                ];
                let base_center = face_center(
                    &mut face_centers,
                    &mut point_sources,
                    &mut next_id,
                    &[vertices[0], vertices[1], vertices[2], vertices[3]],
                )?;
                for verts in pyramid_subdivision(vertices, base_midpoints, base_center) {
                    let fine_cell = alloc_point(&mut next_id)?;
                    for v in verts {
                        refined.add_arrow(fine_cell, v, ());
                    }
                    fine_cells.push((fine_cell, Polarity::Forward));
                }
            }
            CellType::Polygon(_) => {
                let n = vertices.len();
                if n < 3 {
                    return Err(MeshSieveError::RefinementTopologyMismatch {
                        cell,
                        template: "polygon",
                        expected: 3,
                        found: n,
                    });
                }
                let mut edge_midpoints_vec = Vec::with_capacity(n);
                for i in 0..n {
                    edge_midpoints_vec.push(midpoint(
                        &mut edge_midpoints,
                        &mut point_sources,
                        &mut next_id,
                        vertices[i],
                        vertices[(i + 1) % n],
                    )?);
                }
                let center = face_center(
                    &mut face_centers,
                    &mut point_sources,
                    &mut next_id,
                    &vertices,
                )?;
                for i in 0..n {
                    let verts = [
                        vertices[i],
                        edge_midpoints_vec[i],
                        center,
                        edge_midpoints_vec[(i + n - 1) % n],
                    ];
                    let fine_cell = alloc_point(&mut next_id)?;
                    for v in verts {
                        refined.add_arrow(fine_cell, v, ());
                    }
                    fine_cells.push((fine_cell, Polarity::Forward));
                }
            }
            CellType::Polyhedron => {
                let faces = polyhedron_faces(coarse, cell)?;
                let center = alloc_point(&mut next_id)?;
                point_sources.insert(center, vertices.clone());
                for face_vertices in faces {
                    let n = face_vertices.len();
                    if n < 3 {
                        return Err(MeshSieveError::RefinementTopologyMismatch {
                            cell,
                            template: "polyhedron face",
                            expected: 3,
                            found: n,
                        });
                    }
                    let face_center = face_center(
                        &mut face_centers,
                        &mut point_sources,
                        &mut next_id,
                        &face_vertices,
                    )?;
                    let mut edge_midpoints_vec = Vec::with_capacity(n);
                    for i in 0..n {
                        edge_midpoints_vec.push(midpoint(
                            &mut edge_midpoints,
                            &mut point_sources,
                            &mut next_id,
                            face_vertices[i],
                            face_vertices[(i + 1) % n],
                        )?);
                    }
                    for i in 0..n {
                        let v0 = face_vertices[i];
                        let v1 = face_vertices[(i + 1) % n];
                        let m = edge_midpoints_vec[i];
                        for tri in [[v0, m, face_center], [m, v1, face_center]] {
                            let fine_cell = alloc_point(&mut next_id)?;
                            for v in [center, tri[0], tri[1], tri[2]] {
                                refined.add_arrow(fine_cell, v, ());
                            }
                            fine_cells.push((fine_cell, Polarity::Forward));
                        }
                    }
                }
            }
            _ => return Err(MeshSieveError::UnsupportedRefinementCellType { cell, cell_type }),
        }

        refinement_map.push((cell, fine_cells));
    }

    let refined_coordinates = if let Some(coords) = coordinates {
        Some(build_refined_coordinates(
            &refined,
            coords,
            &point_sources,
            &refinement_map,
        )?)
    } else {
        None
    };

    let mut anchors = TopologicalAnchors::default();
    for (point, sources) in &point_sources {
        let kind =
            if sources.len() == 2 && hinted_edge_keys.contains(&edge_key(sources[0], sources[1])) {
                AnchorKind::Hanging
            } else {
                AnchorKind::Refined
            };
        anchors.insert(*point, sources.iter().copied(), kind);
    }
    let hanging_constraints = constraints_from_topological_anchors(&anchors, 1);

    Ok(RefinedMesh {
        sieve: refined,
        cell_refinement: refinement_map,
        coordinates: refined_coordinates,
        anchors,
        hanging_constraints,
    })
}

/// Refine a mesh and return a full cell→face→edge→vertex topology.
///
/// This path is suitable for mixed-element meshes that already carry edges and faces.
/// The refined mesh includes interpolated edges/faces that are kept consistent across
/// shared boundaries via canonical vertex keys.
pub fn refine_mesh_full_topology<S>(
    coarse: &mut impl Sieve<Point = PointId>,
    cell_types: &Section<CellType, S>,
) -> Result<RefinedMeshWithTopology, MeshSieveError>
where
    S: Storage<CellType>,
{
    refine_mesh_full_topology_with_options::<S, crate::data::storage::VecStorage<f64>>(
        coarse,
        cell_types,
        None,
        RefineOptions::default(),
    )
}

/// Refine a mesh and return a full cell→face→edge→vertex topology with options.
pub fn refine_mesh_full_topology_with_options<S, Cs>(
    coarse: &mut impl Sieve<Point = PointId>,
    cell_types: &Section<CellType, S>,
    coordinates: Option<&Coordinates<f64, Cs>>,
    options: RefineOptions,
) -> Result<RefinedMeshWithTopology, MeshSieveError>
where
    S: Storage<CellType>,
    Cs: Storage<f64>,
{
    let cell_only = cell_only_section(coarse, cell_types)?;
    let mut refined = refine_mesh_with_options(coarse, &cell_only, coordinates, options)?;
    let mut refined_cell_types = build_refined_cell_types(&refined, &cell_only)?;
    interpolate_edges_faces(&mut refined.sieve, &mut refined_cell_types)?;

    Ok(RefinedMeshWithTopology {
        sieve: refined.sieve,
        cell_refinement: refined.cell_refinement,
        cell_types: refined_cell_types,
        coordinates: refined.coordinates,
        anchors: refined.anchors,
        hanging_constraints: refined.hanging_constraints,
    })
}

/// Refine a mesh and update point ownership metadata consistently.
///
/// Ownership for new points is inherited from the owning rank of the coarse
/// cell that spawned them. When a new point is shared by multiple refined cells,
/// the smallest owner rank encountered is retained to guarantee deterministic
/// ownership assignment.
pub fn refine_mesh_with_ownership<S>(
    coarse: &mut impl Sieve<Point = PointId>,
    cell_types: &Section<CellType, S>,
    ownership: &PointOwnership,
    my_rank: usize,
) -> Result<RefinedMeshWithOwnership, MeshSieveError>
where
    S: Storage<CellType>,
{
    let refined = refine_mesh_with_options::<S, crate::data::storage::VecStorage<f64>>(
        coarse,
        cell_types,
        None,
        RefineOptions::default(),
    )?;
    let mut refined_ownership = ownership.clone();

    for (cell, fine_cells) in &refined.cell_refinement {
        let owner = refined_ownership.owner_or_err(*cell)?;
        for (fine_cell, _) in fine_cells {
            refined_ownership.set_owner_min(*fine_cell, owner, my_rank)?;
            for p in refined.sieve.cone_points(*fine_cell) {
                if refined_ownership.entry(p).is_none() {
                    refined_ownership.set_owner_min(p, owner, my_rank)?;
                }
            }
        }
    }

    let refined_ownership = refined_ownership.filtered_to_points(refined.sieve.points())?;
    debug_validate_overlap_ownership_topology(&refined.sieve, &refined_ownership, None, my_rank)?;

    Ok(RefinedMeshWithOwnership {
        sieve: refined.sieve,
        cell_refinement: refined.cell_refinement,
        ownership: refined_ownership,
        coordinates: refined.coordinates,
        anchors: refined.anchors,
        hanging_constraints: refined.hanging_constraints,
    })
}

/// Refine a mesh and update point ownership metadata with optional geometry checks.
pub fn refine_mesh_with_ownership_and_options<S, Cs>(
    coarse: &mut impl Sieve<Point = PointId>,
    cell_types: &Section<CellType, S>,
    ownership: &PointOwnership,
    my_rank: usize,
    coordinates: Option<&Coordinates<f64, Cs>>,
    options: RefineOptions,
) -> Result<RefinedMeshWithOwnership, MeshSieveError>
where
    S: Storage<CellType>,
    Cs: Storage<f64>,
{
    let refined = refine_mesh_with_options(coarse, cell_types, coordinates, options)?;
    let mut refined_ownership = ownership.clone();

    for (cell, fine_cells) in &refined.cell_refinement {
        let owner = refined_ownership.owner_or_err(*cell)?;
        for (fine_cell, _) in fine_cells {
            refined_ownership.set_owner_min(*fine_cell, owner, my_rank)?;
            for p in refined.sieve.cone_points(*fine_cell) {
                if refined_ownership.entry(p).is_none() {
                    refined_ownership.set_owner_min(p, owner, my_rank)?;
                }
            }
        }
    }

    let refined_ownership = refined_ownership.filtered_to_points(refined.sieve.points())?;
    debug_validate_overlap_ownership_topology(&refined.sieve, &refined_ownership, None, my_rank)?;

    Ok(RefinedMeshWithOwnership {
        sieve: refined.sieve,
        cell_refinement: refined.cell_refinement,
        ownership: refined_ownership,
        coordinates: refined.coordinates,
        anchors: refined.anchors,
        hanging_constraints: refined.hanging_constraints,
    })
}

fn edge_key(a: PointId, b: PointId) -> (PointId, PointId) {
    if a < b { (a, b) } else { (b, a) }
}

fn collect_hinted_edge_keys(hints: Option<&AnisotropicSplitHints>) -> HashSet<(PointId, PointId)> {
    let mut keys = HashSet::new();
    if let Some(hints) = hints {
        for edges in hints.edge_splits.values() {
            for [a, b] in edges {
                keys.insert(edge_key(*a, *b));
            }
        }
    }
    keys
}

fn triangle_split_indices(
    cell: PointId,
    vertices: [PointId; 3],
    hints: Option<&AnisotropicSplitHints>,
) -> Vec<usize> {
    let Some(edges) = hints.and_then(|h| h.edge_splits.get(&cell)) else {
        return vec![0, 1, 2];
    };
    let canonical = [
        edge_key(vertices[0], vertices[1]),
        edge_key(vertices[1], vertices[2]),
        edge_key(vertices[2], vertices[0]),
    ];
    let requested: HashSet<_> = edges.iter().map(|[a, b]| edge_key(*a, *b)).collect();
    let mut out = Vec::new();
    for (idx, key) in canonical.iter().enumerate() {
        if requested.contains(key) {
            out.push(idx);
        }
    }
    if out.is_empty() { vec![0, 1, 2] } else { out }
}

fn quadrilateral_split_indices(
    cell: PointId,
    vertices: [PointId; 4],
    hints: Option<&AnisotropicSplitHints>,
) -> Vec<usize> {
    let Some(edges) = hints.and_then(|h| h.edge_splits.get(&cell)) else {
        return vec![0, 1, 2, 3];
    };
    let canonical = [
        edge_key(vertices[0], vertices[1]),
        edge_key(vertices[1], vertices[2]),
        edge_key(vertices[2], vertices[3]),
        edge_key(vertices[3], vertices[0]),
    ];
    let requested: HashSet<_> = edges.iter().map(|[a, b]| edge_key(*a, *b)).collect();
    let mut out = Vec::new();
    for (idx, key) in canonical.iter().enumerate() {
        if requested.contains(key) {
            out.push(idx);
        }
    }
    if out.is_empty() {
        vec![0, 1, 2, 3]
    } else {
        out
    }
}

fn anisotropic_triangle_subdivision(
    vertices: [PointId; 3],
    midpoints: [PointId; 3],
    split_indices: &[usize],
) -> Vec<[PointId; 3]> {
    let [v0, v1, v2] = vertices;
    let [m01, m12, m20] = midpoints;
    if split_indices.len() == 1 {
        match split_indices[0] {
            0 => return vec![[v0, m01, v2], [m01, v1, v2]],
            1 => return vec![[v1, m12, v0], [m12, v2, v0]],
            2 => return vec![[v2, m20, v1], [m20, v0, v1]],
            _ => {}
        }
    }
    triangle_subdivision(vertices, midpoints).to_vec()
}

fn anisotropic_quadrilateral_subdivision(
    vertices: [PointId; 4],
    midpoints: [PointId; 4],
    center: PointId,
    split_indices: &[usize],
) -> Vec<[PointId; 4]> {
    let [v0, v1, v2, v3] = vertices;
    let [m01, m12, m23, m30] = midpoints;
    let splits: HashSet<_> = split_indices.iter().copied().collect();
    if splits.len() == 2 && splits.contains(&0) && splits.contains(&2) {
        return vec![[v0, m01, m23, v3], [m01, v1, v2, m23]];
    }
    if splits.len() == 2 && splits.contains(&1) && splits.contains(&3) {
        return vec![[v0, v1, m12, m30], [m30, m12, v2, v3]];
    }
    quadrilateral_subdivision(vertices, midpoints, center).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 midpoint(
    edge_midpoints: &mut HashMap<(PointId, PointId), PointId>,
    point_sources: &mut HashMap<PointId, Vec<PointId>>,
    next_id: &mut u64,
    a: PointId,
    b: PointId,
) -> Result<PointId, MeshSieveError> {
    let key = if a < b { (a, b) } else { (b, a) };
    if let Some(p) = edge_midpoints.get(&key) {
        return Ok(*p);
    }
    let p = alloc_point(next_id)?;
    edge_midpoints.insert(key, p);
    point_sources.insert(p, vec![a, b]);
    Ok(p)
}

fn face_center(
    face_centers: &mut HashMap<Vec<PointId>, PointId>,
    point_sources: &mut HashMap<PointId, Vec<PointId>>,
    next_id: &mut u64,
    vertices: &[PointId],
) -> Result<PointId, MeshSieveError> {
    let mut key = vertices.to_vec();
    key.sort();
    if let Some(p) = face_centers.get(&key) {
        return Ok(*p);
    }
    let p = alloc_point(next_id)?;
    face_centers.insert(key, p);
    point_sources.insert(p, vertices.to_vec());
    Ok(p)
}

fn expected_vertex_count(cell_type: CellType) -> Option<usize> {
    match cell_type {
        CellType::Triangle => Some(3),
        CellType::Quadrilateral => Some(4),
        CellType::Tetrahedron => Some(4),
        CellType::Hexahedron => Some(8),
        CellType::Prism => Some(6),
        CellType::Pyramid => Some(5),
        CellType::Polygon(n) => Some(n as usize),
        _ => None,
    }
}

fn cell_vertices_from_closure(
    coarse: &mut impl Sieve<Point = PointId>,
    cell: PointId,
) -> Result<Vec<PointId>, MeshSieveError> {
    let mut vertices = Vec::new();
    for p in coarse.closure(std::iter::once(cell)) {
        if coarse.cone_points(p).next().is_none() {
            vertices.push(p);
        }
    }
    vertices.sort();
    vertices.dedup();
    if vertices.len() < 3 {
        return Err(MeshSieveError::RefinementTopologyMismatch {
            cell,
            template: "cell",
            expected: 3,
            found: vertices.len(),
        });
    }
    Ok(vertices)
}

fn polyhedron_faces(
    coarse: &mut impl Sieve<Point = PointId>,
    cell: PointId,
) -> Result<Vec<Vec<PointId>>, MeshSieveError> {
    let faces: Vec<PointId> = coarse.cone_points(cell).collect();
    if faces.is_empty() {
        return Err(MeshSieveError::RefinementTopologyMismatch {
            cell,
            template: "polyhedron",
            expected: 1,
            found: 0,
        });
    }
    if faces
        .iter()
        .all(|face| coarse.cone_points(*face).next().is_none())
    {
        return Err(MeshSieveError::InvalidGeometry(
            "polyhedron refinement requires explicit face entities".into(),
        ));
    }
    let mut out = Vec::with_capacity(faces.len());
    for face in faces {
        let vertices = face_vertices(coarse, face)?;
        out.push(vertices);
    }
    Ok(out)
}

fn face_vertices(
    coarse: &mut impl Sieve<Point = PointId>,
    face: PointId,
) -> Result<Vec<PointId>, MeshSieveError> {
    let cone: Vec<PointId> = coarse.cone_points(face).collect();
    if cone.is_empty() {
        return Err(MeshSieveError::RefinementTopologyMismatch {
            cell: face,
            template: "face",
            expected: 3,
            found: 0,
        });
    }
    if cone.iter().all(|p| coarse.cone_points(*p).next().is_none()) {
        return Ok(cone);
    }
    if cone.iter().all(|p| coarse.cone_points(*p).count() == 2)
        && let Some(ordered) = ordered_vertices_from_edges(coarse, &cone)
    {
        return Ok(ordered);
    }
    let mut vertices = Vec::new();
    for p in coarse.closure(std::iter::once(face)) {
        if coarse.cone_points(p).next().is_none() {
            vertices.push(p);
        }
    }
    vertices.sort();
    vertices.dedup();
    if vertices.len() < 3 {
        return Err(MeshSieveError::RefinementTopologyMismatch {
            cell: face,
            template: "face",
            expected: 3,
            found: vertices.len(),
        });
    }
    Ok(vertices)
}

fn collect_top_cells<CtSt>(
    cell_types: &Section<CellType, CtSt>,
    mesh_dimension: Option<u32>,
) -> Result<Vec<(PointId, CellType)>, MeshSieveError>
where
    CtSt: Storage<CellType>,
{
    let mut cells = Vec::new();
    let mut max_dim: Option<u8> = None;
    let mesh_dim = mesh_dimension.and_then(|dim| u8::try_from(dim).ok());
    let mut has_mesh_dim = false;

    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];
        let dim = cell_type.dimension();
        if mesh_dim == Some(dim) {
            has_mesh_dim = true;
        }
        max_dim = Some(max_dim.map_or(dim, |current| current.max(dim)));
        cells.push((point, cell_type));
    }

    let target_dim = if has_mesh_dim {
        mesh_dim.unwrap_or(0)
    } else {
        max_dim.unwrap_or(0)
    };
    cells.retain(|(_, cell_type)| cell_type.dimension() == target_dim);
    Ok(cells)
}

fn cell_only_section<CtSt>(
    sieve: &mut impl Sieve<Point = PointId>,
    cell_types: &Section<CellType, CtSt>,
) -> Result<Section<CellType, VecStorage<CellType>>, MeshSieveError>
where
    CtSt: Storage<CellType>,
{
    let mesh_dimension = crate::topology::utils::dimension(sieve)?;
    let cells = collect_top_cells(cell_types, Some(mesh_dimension))?;
    let mut atlas = Atlas::default();
    for (cell, _) in &cells {
        atlas
            .try_insert(*cell, 1)
            .map_err(|e| MeshSieveError::AtlasInsertionFailed(*cell, Box::new(e)))?;
    }
    let mut filtered = Section::<CellType, VecStorage<CellType>>::new(atlas);
    for (cell, cell_type) in cells {
        filtered.try_set(cell, &[cell_type])?;
    }
    Ok(filtered)
}

fn build_refined_cell_types<CtSt>(
    refined: &RefinedMesh,
    coarse_cells: &Section<CellType, CtSt>,
) -> Result<Section<CellType, VecStorage<CellType>>, MeshSieveError>
where
    CtSt: Storage<CellType>,
{
    let mut cell_map = HashMap::new();
    for (cell, cell_slice) in coarse_cells.iter() {
        if cell_slice.len() != 1 {
            return Err(MeshSieveError::SliceLengthMismatch {
                point: cell,
                expected: 1,
                found: cell_slice.len(),
            });
        }
        cell_map.insert(cell, cell_slice[0]);
    }

    let mut atlas = Atlas::default();
    let mut points = HashSet::new();
    for p in refined.sieve.points() {
        if points.insert(p) {
            atlas
                .try_insert(p, 1)
                .map_err(|e| MeshSieveError::AtlasInsertionFailed(p, Box::new(e)))?;
        }
    }
    let mut cell_types = Section::<CellType, VecStorage<CellType>>::new(atlas);

    for p in refined.sieve.points() {
        if refined.sieve.cone_points(p).next().is_none() {
            cell_types.try_set(p, &[CellType::Vertex])?;
        }
    }

    for (cell, fine_cells) in &refined.cell_refinement {
        let cell_type = *cell_map
            .get(cell)
            .ok_or(MeshSieveError::PointNotInAtlas(*cell))?;
        for (fine_cell, _) in fine_cells {
            let refined_type = match cell_type {
                CellType::Polygon(_) => {
                    let count = refined.sieve.cone_points(*fine_cell).count();
                    let count = u8::try_from(count).map_err(|_| {
                        MeshSieveError::InvalidGeometry(format!(
                            "refined polygon has too many vertices: {count}"
                        ))
                    })?;
                    CellType::Polygon(count)
                }
                _ => cell_type,
            };
            cell_types.try_set(*fine_cell, &[refined_type])?;
        }
    }

    Ok(cell_types)
}

fn build_refined_coordinates<Cs>(
    refined: &MeshSieve,
    coarse_coords: &Coordinates<f64, Cs>,
    point_sources: &HashMap<PointId, Vec<PointId>>,
    refinement_map: &RefinementMap,
) -> Result<Coordinates<f64, VecStorage<f64>>, MeshSieveError>
where
    Cs: Storage<f64>,
{
    let dimension = coarse_coords.embedding_dimension();
    let mut atlas = Atlas::default();
    let mut vertices = Vec::new();
    for p in refined.points() {
        if refined.cone_points(p).next().is_none() {
            atlas.try_insert(p, dimension)?;
            vertices.push(p);
        }
    }

    let mut section = Section::<f64, VecStorage<f64>>::new(atlas);
    for vertex in vertices {
        if coarse_coords.section().atlas().contains(vertex) {
            let values = coarse_coords.try_restrict(vertex)?;
            section.try_set(vertex, values)?;
            continue;
        }
        let sources = point_sources.get(&vertex).ok_or_else(|| {
            MeshSieveError::InvalidGeometry(format!(
                "missing coordinate sources for refined vertex {vertex:?}"
            ))
        })?;
        if sources.is_empty() {
            return Err(MeshSieveError::InvalidGeometry(format!(
                "missing coordinate sources for refined vertex {vertex:?}"
            )));
        }
        let mut accum = vec![0.0; dimension];
        for source in sources {
            let values = coarse_coords.try_restrict(*source)?;
            for (idx, value) in values.iter().enumerate().take(dimension) {
                accum[idx] += value;
            }
        }
        let denom = sources.len() as f64;
        for value in &mut accum {
            *value /= denom;
        }
        section.try_set(vertex, &accum)?;
    }

    let mut out =
        Coordinates::from_section(coarse_coords.topological_dimension(), dimension, section)?;
    if let Some(high_order) = coarse_coords.high_order() {
        let mut ho_entries = Vec::new();
        for (cell, fine_cells) in refinement_map {
            if !high_order.section().atlas().contains(*cell) {
                continue;
            }
            let values = high_order.section().try_restrict(*cell)?;
            for (fine_cell, _) in fine_cells {
                ho_entries.push((*fine_cell, values.to_vec()));
            }
        }
        if !ho_entries.is_empty() {
            let mut ho_atlas = Atlas::default();
            for (cell, values) in &ho_entries {
                ho_atlas.try_insert(*cell, values.len())?;
            }
            let mut ho_section = Section::<f64, VecStorage<f64>>::new(ho_atlas);
            for (cell, values) in ho_entries {
                ho_section.try_set(cell, &values)?;
            }
            let high_order = HighOrderCoordinates::from_section(dimension, ho_section)?;
            out.set_high_order(high_order)?;
        }
    }
    Ok(out)
}

fn cell_vertices(
    coarse: &mut impl Sieve<Point = PointId>,
    cell: PointId,
    expected: usize,
) -> Result<Vec<PointId>, MeshSieveError> {
    let cone: Vec<PointId> = coarse.cone_points(cell).collect();
    if cone.len() == expected && cone.iter().all(|p| coarse.cone_points(*p).next().is_none()) {
        return Ok(cone);
    }

    if (expected == 3 || expected == 4)
        && !cone.is_empty()
        && cone.iter().all(|p| coarse.cone_points(*p).count() == 2)
        && let Some(ordered) = ordered_vertices_from_edges(coarse, &cone)
        && ordered.len() == expected
    {
        return Ok(ordered);
    }

    let mut vertices = Vec::new();
    for p in coarse.closure(std::iter::once(cell)) {
        if coarse.cone_points(p).next().is_none() {
            vertices.push(p);
        }
    }
    vertices.sort();
    vertices.dedup();
    if vertices.len() != expected {
        return Err(MeshSieveError::RefinementTopologyMismatch {
            cell,
            template: "cell",
            expected,
            found: vertices.len(),
        });
    }
    Ok(vertices)
}

fn ordered_vertices_from_edges(
    coarse: &mut impl Sieve<Point = PointId>,
    edges: &[PointId],
) -> Option<Vec<PointId>> {
    let mut edge_vertices = Vec::with_capacity(edges.len());
    for edge in edges {
        let verts: Vec<PointId> = coarse.cone_points(*edge).collect();
        if verts.len() != 2 {
            return None;
        }
        edge_vertices.push([verts[0], verts[1]]);
    }
    let mut ordered = Vec::with_capacity(edges.len());
    let first = edge_vertices.first()?.to_owned();
    ordered.push(first[0]);
    ordered.push(first[1]);
    for verts in edge_vertices.iter().skip(1) {
        let last = *ordered.last()?;
        if verts[0] == last {
            ordered.push(verts[1]);
        } else if verts[1] == last {
            ordered.push(verts[0]);
        } else {
            return None;
        }
    }
    if ordered.len() > 1 && ordered.first()? != ordered.last()? {
        return None;
    }
    ordered.pop();
    Some(ordered)
}