dyntri_core/triangulation/

checks.rs

use crate::triangulation::boundary_map::get_face;
use crate::triangulation::simplices::{
    CausalSimplex, CausalSimplexKind, canonical_oriented_simplex, canonical_simplex,
};
use crate::triangulation::{
    CausalTriangulation, CausalTriangulation2D, CausalTriangulation3D, CausalTriangulation4D,
};
use crate::triangulation::{Triangulation, Triangulation2D, Triangulation3D, Triangulation4D};

#[derive(Debug, thiserror::Error)]
#[error("Neighbour {to} not correctly connected back to {from}.")]
pub struct SimplexConnectivityError {
    from: usize,
    to: usize,
}

/// Checks to check the validity of the [`Triangulation`] structure
impl<const N: usize> Triangulation<N> {
    /// Check if the [`Triangulation`] is properly biconnective, i.e. that each simplex
    /// neighbour points back to the 4-simplex at the correct index.
    pub fn check_biconnectivity(&self) -> Result<(), SimplexConnectivityError> {
        for (label, simplex) in self.simplices.iter().enumerate() {
            for (nbr_label, backindex) in simplex.get_neighbours_backindices() {
                if self.get_neighbour(nbr_label, backindex) != label {
                    return Err(SimplexConnectivityError {
                        from: label,
                        to: nbr_label,
                    });
                }
            }
        }
        Ok(())
    }
}

/// Checks to check the validity of the [`CausalTriangulation`] structure
impl<K, const N: usize> CausalTriangulation<K, N> {
    /// Check if the [`CausalTriangulation`] is properly biconnective, i.e. that each simplex
    /// neighbour points back to the 4-simplex at the correct index.
    pub fn check_biconnectivity(&self) -> Result<(), SimplexConnectivityError> {
        for (label, simplex) in self.simplices.iter().enumerate() {
            for (nbr_label, backindex) in simplex.get_neighbours_backindices() {
                if self.get_neighbour(nbr_label, backindex) != label {
                    return Err(SimplexConnectivityError {
                        from: label,
                        to: nbr_label,
                    });
                }
            }
        }
        Ok(())
    }
}

#[derive(Debug, thiserror::Error)]
pub enum SimplexCheckError {
    #[error("Vertex count mismatch: expected {expected}, found {found}")]
    VertexCountMismatch { expected: usize, found: usize },

    #[error("Face count mismatch: expected {expected}, found {found} (2*N_(D-1) = (D+1)*N_(D))")]
    FaceCountMismatch { expected: usize, found: usize },

    #[error("Dehn-Sommerville relation violated")]
    DehnSommervilleViolation,
}

impl Triangulation2D {
    /// Checks that the number of simplices are all correct
    pub fn check_simplex_counts(&self) -> Result<(), SimplexCheckError> {
        let n0 = self.count_vertices();
        let n1 = self.count_n1();

        if n0 != self.num_vertices() {
            return Err(SimplexCheckError::VertexCountMismatch {
                expected: n0,
                found: self.num_vertices(),
            });
        }
        if n1 != self.num_faces() {
            return Err(SimplexCheckError::FaceCountMismatch {
                expected: n1,
                found: self.num_faces(),
            });
        }

        Ok(())
    }
}

impl Triangulation3D {
    /// Checks that the number of simplices are all correct
    pub fn check_simplex_counts(&self) -> Result<(), SimplexCheckError> {
        let n0 = self.count_vertices();
        let n1 = self.count_n1();
        let n2 = self.count_n2();

        if n0 != self.num_vertices() {
            return Err(SimplexCheckError::VertexCountMismatch {
                expected: n0,
                found: self.num_vertices(),
            });
        }
        if n2 != self.num_faces() {
            return Err(SimplexCheckError::FaceCountMismatch {
                expected: n2,
                found: self.num_faces(),
            });
        }
        if 2 * n0 + n2 != 2 * n1 {
            return Err(SimplexCheckError::DehnSommervilleViolation);
        }

        Ok(())
    }
}

impl Triangulation4D {
    /// Checks that the number of simplices are all correct
    pub fn check_simplex_counts(&self) -> Result<(), SimplexCheckError> {
        let n0 = self.count_vertices();
        let n1 = self.count_n1();
        let n2 = self.count_n2();
        let n3 = self.count_n3();

        if n0 != self.num_vertices() {
            return Err(SimplexCheckError::VertexCountMismatch {
                expected: n0,
                found: self.num_vertices(),
            });
        }
        if n3 != self.num_faces() {
            return Err(SimplexCheckError::FaceCountMismatch {
                expected: n3,
                found: self.num_faces(),
            });
        }
        if 2 * n1 + 2 * n3 != 3 * n2 {
            return Err(SimplexCheckError::DehnSommervilleViolation);
        }

        Ok(())
    }
}

impl CausalTriangulation2D {
    /// Checks that the number of simplices are all correct
    pub fn check_simplex_counts(&self) -> Result<(), SimplexCheckError> {
        let n0 = self.count_vertices();
        let n1 = self.count_n1();

        if n0 != self.num_vertices() {
            return Err(SimplexCheckError::VertexCountMismatch {
                expected: n0,
                found: self.num_vertices(),
            });
        }
        if n1 != self.num_faces() {
            return Err(SimplexCheckError::FaceCountMismatch {
                expected: n1,
                found: self.num_faces(),
            });
        }

        Ok(())
    }
}

impl CausalTriangulation3D {
    /// Checks that the number of simplices are all correct
    pub fn check_simplex_counts(&self) -> Result<(), SimplexCheckError> {
        let n0 = self.count_vertices();
        let n1 = self.count_n1();
        let n2 = self.count_n2();

        if n0 != self.num_vertices() {
            return Err(SimplexCheckError::VertexCountMismatch {
                expected: n0,
                found: self.num_vertices(),
            });
        }
        if n2 != self.num_faces() {
            return Err(SimplexCheckError::FaceCountMismatch {
                expected: n2,
                found: self.num_faces(),
            });
        }
        if 2 * n0 + n2 != 2 * n1 {
            return Err(SimplexCheckError::DehnSommervilleViolation);
        }

        Ok(())
    }
}

impl CausalTriangulation4D {
    /// Checks that the number of simplices are all correct
    pub fn check_simplex_counts(&self) -> Result<(), SimplexCheckError> {
        let n0 = self.count_vertices();
        let n1 = self.count_n1();
        let n2 = self.count_n2();
        let n3 = self.count_n3();

        if n0 != self.num_vertices() {
            return Err(SimplexCheckError::VertexCountMismatch {
                expected: n0,
                found: self.num_vertices(),
            });
        }
        if n3 != self.num_faces() {
            return Err(SimplexCheckError::FaceCountMismatch {
                expected: n3,
                found: self.num_faces(),
            });
        }
        if 2 * n1 + 2 * n3 != 3 * n2 {
            return Err(SimplexCheckError::DehnSommervilleViolation);
        }

        Ok(())
    }
}

#[derive(Debug, thiserror::Error)]
pub enum VertexConsistencyError<const N: usize> {
    #[error("Trianglation has incorrectly matched vertices: {face:?} {nbr_face:?}")]
    VertexMatching {
        face: [usize; N],
        nbr_face: [usize; N],
    },
    #[error("Trianglation has incorrectly oriented vertices: {face:?} {nbr_face:?}")]
    Orientation {
        face: [usize; N],
        nbr_face: [usize; N],
    },
}

macro_rules! impl_vertex_consistency_check {
    ($triangulation:ty, $dim:expr) => {
        impl $triangulation {
            /// Checks whether the vertices shared between neighbouring simplices are the same
            pub fn check_vertex_consistency(&self) -> Result<(), VertexConsistencyError<$dim>> {
                for simplex in self.simplices.iter() {
                    for (index, (nbr_label, backindex)) in
                        simplex.get_neighbours_backindices().into_iter().enumerate()
                    {
                        let vertices = simplex.vertices;
                        let nbr_vertices = self.get_simplex_vertices(nbr_label);

                        let face = get_face(vertices, index);
                        let nbr_face = get_face(nbr_vertices, backindex);
                        let face_canon = canonical_simplex(face);
                        let nbr_face_canon = canonical_simplex(nbr_face);

                        if face_canon != nbr_face_canon {
                            return Err(VertexConsistencyError::VertexMatching { face, nbr_face });
                        }
                    }
                }

                Ok(())
            }

            /// Checks whether the vertices shared between neighbouring simplices are the same
            /// and they are consistently ordered.
            pub fn check_vertex_consistency_ordered(
                &self,
            ) -> Result<(), VertexConsistencyError<$dim>> {
                for simplex in self.simplices.iter() {
                    for (index, (nbr_label, backindex)) in
                        simplex.get_neighbours_backindices().into_iter().enumerate()
                    {
                        let vertices = simplex.vertices;
                        let nbr_vertices = self.get_simplex_vertices(nbr_label);

                        let face = get_face(vertices, index);
                        let nbr_face = get_face(nbr_vertices, backindex);
                        let (mut parity, face_canon) = canonical_oriented_simplex(face);
                        let (mut nbr_parity, nbr_face_canon) = canonical_oriented_simplex(nbr_face);
                        parity = !(parity ^ index.is_multiple_of(2));
                        nbr_parity = !(nbr_parity ^ backindex.is_multiple_of(2));
                        if face_canon != nbr_face_canon {
                            return Err(VertexConsistencyError::VertexMatching { face, nbr_face });
                        }
                        // Check parity is opposite
                        if nbr_parity == parity {
                            return Err(VertexConsistencyError::Orientation { face, nbr_face });
                        }
                    }
                }

                Ok(())
            }
        }
    };
}

impl_vertex_consistency_check!(Triangulation2D, 2);
impl_vertex_consistency_check!(Triangulation3D, 3);
impl_vertex_consistency_check!(Triangulation4D, 4);

impl_vertex_consistency_check!(CausalTriangulation2D, 2);
impl_vertex_consistency_check!(CausalTriangulation3D, 3);
impl_vertex_consistency_check!(CausalTriangulation4D, 4);

//=============================
// Checks of absolute ordering
//=============================

#[derive(Debug, Clone, thiserror::Error)]
#[error(
    "Simplex ({label}: {simplex:?}) has incorrect absolute order of neighbour \
     ({label_nbr}: {simplex_nbr:?}) or the incorrect SimplexKind"
)]
pub struct NeighbourAbsoluteOrderingError<K, const N: usize> {
    label: usize,
    simplex: CausalSimplex<K, N>,
    label_nbr: usize,
    simplex_nbr: CausalSimplex<K, N>,
}

#[derive(Debug, Clone, thiserror::Error)]
pub enum SimplexKindError<K> {
    #[error("Simplex vertex have time {tfound} which is inconsistent with {tsimplex}")]
    IncorrectTime { tsimplex: u16, tfound: u16 },
    #[error("Simplex kind {kind:?} inconsistent with found signature {signature:?}")]
    InconsistentKind { kind: K, signature: (u8, u8) },
}

impl<K: Copy + CausalSimplexKind, const N: usize> CausalTriangulation<K, N> {
    /// Checks if the triangulation obeys the absolute ordering of simplex neighbours such that
    /// the first (index 0) neighbour is the time-like neighbour and time values are correctly set.
    ///
    /// See also notes on fields of [`CausalSimplex`](crate::triangulation::simplices::CausalSimplex).
    pub fn check_neighbour_absolute_ordering(
        &self,
    ) -> Result<(), NeighbourAbsoluteOrderingError<K, N>> {
        for label in 0..self.num_simplices() {
            let t = self.get_simplex_time(label);
            // Compare to neighbour 0
            let label_nbr = self.get_neighbour(label, 0);
            let tnbr = self.get_simplex_time(label_nbr);

            match self.get_simplex_kind(label).orientation() {
                super::simplices::CausalOrientation::Past => {
                    if (tnbr + 1) % self.num_timeslices() != t {
                        return Err(NeighbourAbsoluteOrderingError {
                            label,
                            simplex: *self.get_simplex(label),
                            label_nbr,
                            simplex_nbr: *self.get_simplex(label_nbr),
                        });
                    }
                }
                super::simplices::CausalOrientation::Space => {
                    if t != tnbr {
                        return Err(NeighbourAbsoluteOrderingError {
                            label,
                            simplex: *self.get_simplex(label),
                            label_nbr,
                            simplex_nbr: *self.get_simplex(label_nbr),
                        });
                    }
                }
                super::simplices::CausalOrientation::Future => {
                    if tnbr != (t + 1) % self.num_timeslices() {
                        return Err(NeighbourAbsoluteOrderingError {
                            label,
                            simplex: *self.get_simplex(label),
                            label_nbr,
                            simplex_nbr: *self.get_simplex(label_nbr),
                        });
                    }
                }
            }
        }
        Ok(())
    }

    /// Returns whether the triangulation obeys the absolute ordering of triangle neighbours
    /// such that the first (index 0) neighbour is the time-like neighbour.
    ///
    /// See also notes on fields of [`CausalTriangle`](super::simplices::CausalTriangle).
    pub fn check_simplex_kinds(&self) -> Result<(), SimplexKindError<K>> {
        for label in 0..self.num_simplices() {
            // Get lower vertex times
            let t_lower = self.get_simplex_time(label);
            // Get upper vertex times
            let t_upper = (t_lower + 1) % self.num_timeslices();

            let mut lower_count = 0u8;
            let mut upper_count = 0u8;
            for vtime in self
                .get_simplex_vertices(label)
                .map(|vlabel| self.get_vertex_time(vlabel))
            {
                if vtime == t_lower {
                    lower_count += 1
                } else if vtime == t_upper {
                    upper_count += 1
                } else {
                    return Err(SimplexKindError::IncorrectTime {
                        tsimplex: t_lower,
                        tfound: vtime,
                    });
                }
            }

            if self.get_simplex_kind(label).signature() != (lower_count, upper_count) {
                return Err(SimplexKindError::InconsistentKind {
                    kind: self.get_simplex_kind(label),
                    signature: (lower_count, upper_count),
                });
            }
        }
        Ok(())
    }
}