mesh-sieve 4.0.2

//! Canonical geometric vertex ordering for supported reference cells.
//!
//! The order returned by this module is independent of point identifiers and
//! of adjacency insertion order.  Topology is used to constrain admissible
//! permutations (when explicit facets are present); coordinates then select a
//! deterministic, positively oriented representative.

use std::collections::HashSet;

use crate::data::coordinates::Coordinates;
use crate::data::storage::Storage;
use crate::mesh_error::MeshSieveError;
use crate::topology::cell_type::CellType;
use crate::topology::point::PointId;
use crate::topology::sieve::OrientedSieve;

/// Canonical vertices and their coordinates in the same order.
#[derive(Clone, Debug, PartialEq)]
pub struct CanonicalCellVertices {
    pub vertices: Vec<PointId>,
    pub coordinates: Vec<[f64; 3]>,
}

/// Options for canonical ordering.  The default tolerance is scale-aware and
/// is only used for degeneracy tests, never for choosing a point-id order.
#[derive(Clone, Copy, Debug)]
pub struct CanonicalOrderingOptions {
    pub tolerance: f64,
}

impl Default for CanonicalOrderingOptions {
    fn default() -> Self {
        Self { tolerance: 1.0e-12 }
    }
}

/// Return the canonical geometric order for one supported cell.
pub fn canonical_cell_vertices<S, St>(
    sieve: &S,
    cell: PointId,
    cell_type: CellType,
    coordinates: &Coordinates<f64, St>,
) -> Result<CanonicalCellVertices, MeshSieveError>
where
    S: OrientedSieve<Point = PointId>,
    St: Storage<f64>,
{
    canonical_cell_vertices_with_options(
        sieve,
        cell,
        cell_type,
        coordinates,
        CanonicalOrderingOptions::default(),
    )
}

/// Canonical ordering with an explicit degeneracy tolerance.
pub fn canonical_cell_vertices_with_options<S, St>(
    sieve: &S,
    cell: PointId,
    cell_type: CellType,
    coordinates: &Coordinates<f64, St>,
    options: CanonicalOrderingOptions,
) -> Result<CanonicalCellVertices, MeshSieveError>
where
    S: OrientedSieve<Point = PointId>,
    St: Storage<f64>,
{
    if !options.tolerance.is_finite() || options.tolerance < 0.0 {
        return Err(MeshSieveError::InvalidGeometry(
            "canonical ordering tolerance must be finite and non-negative".into(),
        ));
    }
    let expected = expected_vertices(cell_type)
        .ok_or(MeshSieveError::UnsupportedCanonicalCellType { cell, cell_type })?;

    // closure_o is orientation-aware.  Leaf points are the actual geometric
    // vertices, regardless of whether the cell is direct or interpolated.
    let mut vertices = Vec::new();
    for (point, _) in sieve.closure_o([cell]) {
        if sieve.cone_o(point).next().is_none() && point != cell {
            vertices.push(point);
        }
    }
    if vertices.len() != expected {
        return Err(MeshSieveError::MalformedCellTopology {
            cell,
            cell_type,
            reason: format!(
                "expected {expected} leaf vertices, found {}",
                vertices.len()
            ),
        });
    }
    let unique: HashSet<_> = vertices.iter().copied().collect();
    if unique.len() != expected {
        return Err(MeshSieveError::MalformedCellTopology {
            cell,
            cell_type,
            reason: "duplicate vertex in closure".into(),
        });
    }

    let mut points = Vec::with_capacity(expected);
    let mut coords = Vec::with_capacity(expected);
    for &point in &vertices {
        let values =
            coordinates
                .try_restrict(point)
                .map_err(|_| MeshSieveError::MalformedCellTopology {
                    cell,
                    cell_type,
                    reason: format!("missing coordinates for vertex {point:?}"),
                })?;
        if values.len() != 2 && values.len() != 3 {
            return Err(MeshSieveError::MalformedCellTopology {
                cell,
                cell_type,
                reason: format!(
                    "coordinate dimension must be 2 or 3, found {}",
                    values.len()
                ),
            });
        }
        let xyz = [values[0], values[1], values.get(2).copied().unwrap_or(0.0)];
        if !xyz.iter().all(|v| v.is_finite()) {
            return Err(MeshSieveError::DegenerateCellGeometry {
                cell,
                cell_type,
                reason: "non-finite vertex coordinate".into(),
            });
        }
        points.push(point);
        coords.push(xyz);
    }

    // Coincident geometric vertices cannot define a unique cell ordering.  Do
    // this check before any permutation search so a tie can never fall back to
    // the point-id order returned by the topology traversal.
    if !coordinates_are_distinct(&coords, options.tolerance) {
        return Err(MeshSieveError::DegenerateCellGeometry {
            cell,
            cell_type,
            reason: "coincident vertex coordinates".into(),
        });
    }

    let facets = topology_facets(sieve, cell, cell_type, &vertices)?;
    if facets.is_empty() {
        match geometric_template(cell_type, &coords, options.tolerance) {
            Some(perm) => {
                let candidate: Vec<[f64; 3]> = perm.iter().map(|&i| coords[i]).collect();
                if positive_and_nondegenerate(cell_type, &candidate, options.tolerance) {
                    return Ok(CanonicalCellVertices {
                        vertices: perm.iter().map(|&i| points[i]).collect(),
                        coordinates: candidate,
                    });
                }
            }
            None if matches!(
                cell_type,
                CellType::Hexahedron | CellType::Prism | CellType::Pyramid
            ) && coordinates_are_distinct(&coords, options.tolerance)
                && affine_rank(&coords, options.tolerance) < 3 =>
            {
                return Err(MeshSieveError::DegenerateCellGeometry {
                    cell,
                    cell_type,
                    reason: "cell vertices do not span a three-dimensional volume".into(),
                });
            }
            None if matches!(
                cell_type,
                CellType::Hexahedron | CellType::Prism | CellType::Pyramid
            ) && coordinates_are_distinct(&coords, options.tolerance)
                && coordinate_span(&coords) > options.tolerance =>
            {
                return Err(MeshSieveError::AmbiguousCellOrdering {
                    cell,
                    cell_type,
                    reason: "direct cell coordinates do not determine a unique reference frame"
                        .into(),
                });
            }
            None => {}
        }
    }
    let mut best: Option<(Vec<usize>, Vec<f64>)> = None;
    let mut permutation = (0..expected).collect::<Vec<_>>();
    enumerate_permutations(&mut permutation, 0, &mut |perm| {
        if !facets_match(cell_type, perm, &vertices, &facets) {
            return;
        }
        let candidate: Vec<[f64; 3]> = perm.iter().map(|&i| coords[i]).collect();
        if !positive_and_nondegenerate(cell_type, &candidate, options.tolerance) {
            return;
        }
        let key: Vec<f64> = candidate.iter().flat_map(|p| p.iter().copied()).collect();
        if best.as_ref().is_none_or(|(_, old)| lex_less(&key, old)) {
            best = Some((perm.to_vec(), key));
        }
    });

    let Some((perm, _)) = best else {
        return Err(MeshSieveError::DegenerateCellGeometry {
            cell,
            cell_type,
            reason: "no positively oriented non-degenerate topology-compatible ordering".into(),
        });
    };
    let ordered_vertices = perm.iter().map(|&i| points[i]).collect();
    let ordered_coordinates = perm.iter().map(|&i| coords[i]).collect();
    Ok(CanonicalCellVertices {
        vertices: ordered_vertices,
        coordinates: ordered_coordinates,
    })
}

/// Return the vertex order implied by oriented connectivity when coordinates
/// are intentionally unavailable.
///
/// This is the compatibility path for topology-only interpolation and mesh
/// construction APIs.  It never compares or sorts `PointId` values: the
/// order is the first-arrival order of the oriented cone traversal.  Geometry
/// consumers should use [`canonical_cell_vertices`] instead, because only the
/// coordinate-aware API can make the result invariant under arbitrary cone
/// insertion order.
pub fn topology_cell_vertices<S>(
    sieve: &S,
    cell: PointId,
    cell_type: CellType,
) -> Result<Vec<PointId>, MeshSieveError>
where
    S: OrientedSieve<Point = PointId>,
{
    let expected = match cell_type {
        CellType::Polyhedron => None,
        _ => Some(
            expected_vertices(cell_type)
                .ok_or(MeshSieveError::UnsupportedCanonicalCellType { cell, cell_type })?,
        ),
    };
    let mut vertices = Vec::with_capacity(expected.unwrap_or(0));
    let mut seen = HashSet::new();
    fn visit<S: OrientedSieve<Point = PointId>>(
        sieve: &S,
        root: PointId,
        point: PointId,
        seen: &mut HashSet<PointId>,
        vertices: &mut Vec<PointId>,
    ) {
        if !seen.insert(point) {
            return;
        }
        let children: Vec<_> = sieve.cone_o(point).map(|(child, _)| child).collect();
        if children.is_empty() {
            if point != root {
                vertices.push(point);
            }
            return;
        }
        for child in children {
            visit(sieve, root, child, seen, vertices);
        }
    }
    visit(sieve, cell, cell, &mut seen, &mut vertices);
    if let Some(expected) = expected
        && vertices.len() != expected
    {
        return Err(MeshSieveError::MalformedCellTopology {
            cell,
            cell_type,
            reason: format!(
                "expected {expected} leaf vertices, found {}",
                vertices.len()
            ),
        });
    }
    Ok(vertices)
}

fn coordinates_are_distinct(coords: &[[f64; 3]], tol: f64) -> bool {
    let tol2 = tol * tol;
    coords
        .iter()
        .enumerate()
        .all(|(i, a)| coords.iter().skip(i + 1).all(|b| distance2(a, b) > tol2))
}

fn coordinate_span(coords: &[[f64; 3]]) -> f64 {
    let mut lo = [f64::INFINITY; 3];
    let mut hi = [f64::NEG_INFINITY; 3];
    for p in coords {
        for d in 0..3 {
            lo[d] = lo[d].min(p[d]);
            hi[d] = hi[d].max(p[d]);
        }
    }
    (0..3).map(|d| hi[d] - lo[d]).fold(0.0, f64::max)
}

fn affine_rank(coords: &[[f64; 3]], tol: f64) -> usize {
    let Some(&origin) = coords.first() else {
        return 0;
    };
    let vectors: Vec<_> = coords.iter().skip(1).map(|p| sub3(*p, origin)).collect();
    let scale = coordinate_span(coords).max(1.0);
    let eps = tol * scale;
    let Some(first) = vectors
        .iter()
        .copied()
        .find(|v| distance2(v, &[0.0; 3]) > eps * eps)
    else {
        return 0;
    };
    let Some(second) = vectors
        .iter()
        .copied()
        .map(|v| cross3(first, v))
        .find(|v| distance2(v, &[0.0; 3]) > eps * eps)
    else {
        return 1;
    };
    if vectors.iter().any(|v| dot3(*v, second).abs() > eps * eps) {
        3
    } else {
        2
    }
}

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

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

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

fn topology_facets<S: OrientedSieve<Point = PointId>>(
    sieve: &S,
    cell: PointId,
    cell_type: CellType,
    vertices: &[PointId],
) -> Result<Vec<HashSet<PointId>>, MeshSieveError> {
    let expected_facets = match cell_type {
        CellType::Tetrahedron | CellType::Simplex(3) => 4,
        CellType::Hexahedron => 6,
        CellType::Prism => 5,
        CellType::Pyramid => 5,
        _ => 0,
    };
    if expected_facets == 0 {
        return Ok(Vec::new());
    }
    let children: Vec<_> = sieve.cone_o(cell).map(|(p, _)| p).collect();
    if children.is_empty() {
        return Ok(Vec::new());
    }
    // A direct cell stores its vertices as the immediate cone.  There are no
    // explicit facets in that representation, so leave topology unconstrained
    // and let the coordinates choose the canonical frame.
    if children
        .iter()
        .all(|child| sieve.cone_o(*child).next().is_none())
    {
        return Ok(Vec::new());
    }
    let allowed: HashSet<_> = vertices.iter().copied().collect();
    let mut out = Vec::with_capacity(expected_facets);
    let valid_sizes = match cell_type {
        CellType::Tetrahedron | CellType::Simplex(3) => &[3][..],
        CellType::Hexahedron => &[4][..],
        CellType::Prism | CellType::Pyramid => &[3, 4][..],
        _ => &[][..],
    };
    for child in children {
        if sieve.cone_o(child).next().is_none() {
            continue;
        }
        let mut set = HashSet::new();
        for (p, _) in sieve.closure_o([child]) {
            if sieve.cone_o(p).next().is_none() && p != child {
                set.insert(p);
            }
        }
        if !valid_sizes.contains(&set.len()) {
            continue;
        }
        if !set.is_subset(&allowed) {
            return Err(MeshSieveError::MalformedCellTopology {
                cell,
                cell_type,
                reason: "facet does not have a valid vertex set".into(),
            });
        }
        out.push(set);
    }
    if out.len() != expected_facets {
        return Err(MeshSieveError::MalformedCellTopology {
            cell,
            cell_type,
            reason: format!("expected {expected_facets} facets, found {}", out.len()),
        });
    }
    Ok(out)
}

fn facets_match(
    cell_type: CellType,
    perm: &[usize],
    vertices: &[PointId],
    actual: &[HashSet<PointId>],
) -> bool {
    if actual.is_empty() {
        return true;
    }
    let patterns: &[&[usize]] = match cell_type {
        CellType::Tetrahedron | CellType::Simplex(3) => {
            &[&[0, 1, 2], &[0, 1, 3], &[1, 2, 3], &[0, 2, 3]]
        }
        CellType::Hexahedron => &[
            &[0, 1, 2, 3],
            &[4, 5, 6, 7],
            &[0, 1, 5, 4],
            &[1, 2, 6, 5],
            &[2, 3, 7, 6],
            &[3, 0, 4, 7],
        ],
        CellType::Prism => &[
            &[0, 1, 2],
            &[3, 4, 5],
            &[0, 1, 4, 3],
            &[1, 2, 5, 4],
            &[2, 0, 3, 5],
        ],
        CellType::Pyramid => &[
            &[0, 1, 2, 3],
            &[0, 1, 4],
            &[1, 2, 4],
            &[2, 3, 4],
            &[3, 0, 4],
        ],
        _ => return true,
    };
    let expected: Vec<HashSet<PointId>> = patterns
        .iter()
        .map(|pattern| pattern.iter().map(|&i| vertices[perm[i]]).collect())
        .collect();
    expected
        .iter()
        .all(|set| actual.iter().any(|got| got == set))
}

/// Infer the standard connectivity for direct tensor/product cells from the
/// coordinate frame.  A direct cell has no facet entities, so this is the
/// only geometric information available beyond its vertex set.
fn geometric_template(cell_type: CellType, coords: &[[f64; 3]], tol: f64) -> Option<Vec<usize>> {
    match cell_type {
        CellType::Segment | CellType::Simplex(1) => {
            let mut out = vec![0, 1];
            out.sort_by(|&a, &b| coord_cmp(coords[a], coords[b]));
            Some(out)
        }
        CellType::Triangle | CellType::Simplex(2) | CellType::Quadrilateral => {
            let centroid = coords.iter().fold([0.0; 3], |mut acc, p| {
                for d in 0..3 {
                    acc[d] += p[d];
                }
                acc
            });
            let inv = 1.0 / coords.len() as f64;
            let mut normal = [0.0; 3];
            'outer: for i in 1..coords.len() {
                for j in (i + 1)..coords.len() {
                    normal = cross3(sub3(coords[i], coords[0]), sub3(coords[j], coords[0]));
                    if normal.iter().any(|v| v.abs() > tol) {
                        break 'outer;
                    }
                }
            }
            let axis = (0..3).max_by(|&a, &b| normal[a].abs().total_cmp(&normal[b].abs()))?;
            let (u, v) = match axis {
                0 => (1, 2),
                1 => (0, 2),
                _ => (0, 1),
            };
            let center = [centroid[0] * inv, centroid[1] * inv, centroid[2] * inv];
            let mut out: Vec<usize> = (0..coords.len()).collect();
            out.sort_by(|&a, &b| {
                (coords[a][v] - center[v])
                    .atan2(coords[a][u] - center[u])
                    .total_cmp(&(coords[b][v] - center[v]).atan2(coords[b][u] - center[u]))
            });
            let ordered: Vec<_> = out.iter().map(|&i| coords[i]).collect();
            if !positive_and_nondegenerate(cell_type, &ordered, tol) {
                out.reverse();
            }
            let start = out
                .iter()
                .enumerate()
                .min_by(|a, b| coord_cmp(coords[*a.1], coords[*b.1]))
                .map(|(i, _)| i)
                .unwrap_or(0);
            out.rotate_left(start);
            Some(out)
        }
        CellType::Hexahedron => {
            let mut lo = [f64::INFINITY; 3];
            let mut hi = [f64::NEG_INFINITY; 3];
            for p in coords {
                for d in 0..3 {
                    lo[d] = lo[d].min(p[d]);
                    hi[d] = hi[d].max(p[d]);
                }
            }
            let mut slots = [usize::MAX; 8];
            for (i, p) in coords.iter().enumerate() {
                let mut bits = 0usize;
                for d in 0..3 {
                    let span = hi[d] - lo[d];
                    if span <= tol {
                        return None;
                    }
                    if (p[d] - lo[d]) > 0.5 * span {
                        bits |= 1 << d;
                    }
                }
                let slot = match bits {
                    0 => 0,
                    1 => 1,
                    3 => 2,
                    2 => 3,
                    4 => 4,
                    5 => 5,
                    7 => 6,
                    6 => 7,
                    _ => return None,
                };
                if slots[slot] != usize::MAX {
                    return None;
                }
                slots[slot] = i;
            }
            slots
                .iter()
                .all(|i| *i != usize::MAX)
                .then_some(slots.to_vec())
        }
        CellType::Prism => {
            let zmin = coords.iter().map(|p| p[2]).fold(f64::INFINITY, f64::min);
            let zmax = coords
                .iter()
                .map(|p| p[2])
                .fold(f64::NEG_INFINITY, f64::max);
            let span = zmax - zmin;
            if span <= tol {
                return None;
            }
            let base: Vec<_> = coords
                .iter()
                .enumerate()
                .filter(|(_, p)| p[2] <= zmin + 0.5 * span)
                .map(|(i, _)| i)
                .collect();
            let top: Vec<_> = coords
                .iter()
                .enumerate()
                .filter(|(_, p)| p[2] > zmin + 0.5 * span)
                .map(|(i, _)| i)
                .collect();
            if base.len() != 3 || top.len() != 3 {
                return None;
            }
            let mut best = None;
            let mut perm = base.clone();
            enumerate_permutations(&mut perm, 0, &mut |candidate| {
                let p = candidate.iter().map(|&i| coords[i]).collect::<Vec<_>>();
                if !positive_and_nondegenerate(CellType::Triangle, &p, tol) {
                    return;
                }
                let mut matched = Vec::with_capacity(3);
                for &i in candidate {
                    let Some(j) = top.iter().copied().min_by(|&a, &b| {
                        distance2_xy(coords[i], coords[a])
                            .total_cmp(&distance2_xy(coords[i], coords[b]))
                    }) else {
                        return;
                    };
                    matched.push(j);
                }
                if matched.iter().collect::<HashSet<_>>().len() == 3 {
                    let key = p.iter().flat_map(|x| x.iter().copied()).collect::<Vec<_>>();
                    let replace = best.as_ref().is_none_or(|old: &Vec<usize>| {
                        let old_key = old
                            .iter()
                            .map(|&i| coords[i])
                            .flat_map(|x| x.into_iter())
                            .collect::<Vec<_>>();
                        lex_less(&key, &old_key)
                    });
                    if replace {
                        best = Some(candidate.to_vec());
                    }
                }
            });
            let base = best?;
            let top_order: Vec<_> = base
                .iter()
                .map(|&i| {
                    top.iter()
                        .copied()
                        .min_by(|&a, &b| {
                            distance2_xy(coords[i], coords[a])
                                .total_cmp(&distance2_xy(coords[i], coords[b]))
                        })
                        .unwrap()
                })
                .collect();
            Some(base.into_iter().chain(top_order).collect())
        }
        CellType::Pyramid => {
            let zmin = coords.iter().map(|p| p[2]).fold(f64::INFINITY, f64::min);
            let zmax = coords
                .iter()
                .map(|p| p[2])
                .fold(f64::NEG_INFINITY, f64::max);
            let span = zmax - zmin;
            if span <= tol {
                return None;
            }
            let base: Vec<_> = coords
                .iter()
                .enumerate()
                .filter(|(_, p)| p[2] <= zmin + 0.25 * span)
                .map(|(i, _)| i)
                .collect();
            let apex = coords
                .iter()
                .enumerate()
                .max_by(|(_, a), (_, b)| a[2].total_cmp(&b[2]))
                .map(|(i, _)| i)?;
            if base.len() != 4 {
                return None;
            }
            let mut best = None;
            let mut perm = base;
            enumerate_permutations(&mut perm, 0, &mut |candidate| {
                let mut p = candidate.iter().map(|&i| coords[i]).collect::<Vec<_>>();
                p.push(coords[apex]);
                if positive_and_nondegenerate(CellType::Pyramid, &p, tol) {
                    let key = p.iter().flat_map(|x| x.iter().copied()).collect::<Vec<_>>();
                    let replace = best.as_ref().is_none_or(|old: &Vec<usize>| {
                        let old_key = old
                            .iter()
                            .map(|&i| coords[i])
                            .chain([coords[apex]])
                            .flat_map(|x| x.into_iter())
                            .collect::<Vec<_>>();
                        lex_less(&key, &old_key)
                    });
                    if replace {
                        best = Some(candidate.to_vec());
                    }
                }
            });
            best.map(|base| base.into_iter().chain([apex]).collect())
        }
        _ => None,
    }
}

fn distance2_xy(a: [f64; 3], b: [f64; 3]) -> f64 {
    (a[0] - b[0]).mul_add(a[0] - b[0], (a[1] - b[1]) * (a[1] - b[1]))
}

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 coord_cmp(a: [f64; 3], b: [f64; 3]) -> std::cmp::Ordering {
    a[0].total_cmp(&b[0])
        .then_with(|| a[1].total_cmp(&b[1]))
        .then_with(|| a[2].total_cmp(&b[2]))
}

fn enumerate_permutations<F: FnMut(&[usize])>(values: &mut [usize], at: usize, f: &mut F) {
    if at == values.len() {
        f(values);
        return;
    }
    for i in at..values.len() {
        values.swap(at, i);
        enumerate_permutations(values, at + 1, f);
        values.swap(at, i);
    }
}

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

fn positive_and_nondegenerate(cell_type: CellType, p: &[[f64; 3]], tol: f64) -> bool {
    let scale = p
        .iter()
        .flat_map(|x| x.iter())
        .fold(1.0_f64, |m, x| m.max(x.abs()));
    let eps = tol * scale.max(1.0).powi(3);
    let cross = |a: [f64; 3], b: [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],
        ]
    };
    let sub = |a: [f64; 3], b: [f64; 3]| [a[0] - b[0], a[1] - b[1], a[2] - b[2]];
    let dot = |a: [f64; 3], b: [f64; 3]| a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
    match cell_type {
        CellType::Segment | CellType::Simplex(1) => {
            let delta = sub(p[1], p[0]);
            dot(delta, delta).sqrt() > tol * scale.max(1.0)
        }
        CellType::Triangle | CellType::Simplex(2) | CellType::Quadrilateral => {
            let n = cross(sub(p[1], p[0]), sub(p[2], p[0]));
            let area = dot(n, n).sqrt();
            if area <= tol * scale.max(1.0).powi(2) {
                return false;
            }
            let axis = (0..3)
                .max_by(|&i, &j| n[i].abs().total_cmp(&n[j].abs()))
                .unwrap();
            n[axis] > tol * scale.max(1.0).powi(2)
        }
        CellType::Tetrahedron | CellType::Simplex(3) => {
            dot(cross(sub(p[1], p[0]), sub(p[2], p[0])), sub(p[3], p[0])) > eps
        }
        CellType::Hexahedron => dot(cross(sub(p[1], p[0]), sub(p[3], p[0])), sub(p[4], p[0])) > eps,
        CellType::Prism => dot(cross(sub(p[1], p[0]), sub(p[2], p[0])), sub(p[3], p[0])) > eps,
        CellType::Pyramid => {
            // The diagonal fan must have positive total volume.  Checking
            // only one base diagonal admits self-crossing permutations whose
            // two tetrahedra cancel to zero volume.
            let first = dot(cross(sub(p[1], p[0]), sub(p[2], p[0])), sub(p[4], p[0])) / 6.0;
            let second = dot(cross(sub(p[2], p[0]), sub(p[3], p[0])), sub(p[4], p[0])) / 6.0;
            first + second > eps
        }
        _ => false,
    }
}