Skip to main content

dyntri_core/triangulation/
simplices.rs

1use serde::{Deserialize, Serialize};
2use serde_with::serde_as;
3
4/// Euclidean (N-1)-simplex, built from `N` vertices.
5///
6/// The intended `neighbours` and `vertices` ordering is such that the vertex corresponding
7/// to a neighbouring simplex is the vertex opposite the simplex. Additionally the overall
8/// ordering is supposed to be consistent, such that the simplices are *oriented*.
9#[serde_as]
10#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Serialize, Deserialize)]
11pub struct Simplex<const N: usize> {
12    #[serde_as(as = "[_; N]")]
13    pub neighbours: [usize; N],
14    #[serde_as(as = "[_; N]")]
15    pub vertices: [usize; N],
16}
17
18impl<const N: usize> Default for Simplex<N> {
19    fn default() -> Self {
20        Self {
21            neighbours: [0; N],
22            vertices: [0; N],
23        }
24    }
25}
26
27pub type Triangle = Simplex<3>;
28pub type Simplex2 = Simplex<3>;
29pub type Tetrahedron = Simplex<4>;
30pub type Simplex3 = Simplex<4>;
31pub type Simplex4 = Simplex<5>;
32
33/// Methods for working with [`Simplex<N>`] when `neighbours` are encoded by `half-faces`,
34/// see [`Triangulation`](super::triangulation::Triangulation)
35impl<const N: usize> Simplex<N> {
36    /// Returns the label of the neighbouring (N-1)-simplex on the `index` side.
37    ///
38    /// # Panic
39    /// This method will panic `index >= N`
40    #[inline]
41    pub fn get_neighbour(&self, index: usize) -> usize {
42        let nbr = self.neighbours[index];
43        nbr / N
44    }
45
46    /// Returns the label and backindex of neighbouring (N-1)-simplex on the `index` side.
47    ///
48    /// # Panic
49    /// This method will panic `index >= N`
50    #[inline]
51    pub fn get_neighbour_backindex(&self, index: usize) -> (usize, usize) {
52        let nbr = self.neighbours[index];
53        (nbr / N, nbr % N)
54    }
55
56    /// Returns the neighbouring (N-1)-simplex labels of the [`Simplex`]
57    #[inline]
58    pub fn get_neighbours(&self) -> [usize; N] {
59        let adjacent_edges = self.neighbours;
60        adjacent_edges.map(|nbr| nbr / N)
61    }
62
63    /// Returns the neighbouring (N-1)-simplex labels and back indices in tuple pairs
64    #[inline]
65    pub fn get_neighbours_backindices(&self) -> [(usize, usize); N] {
66        let adjacent_edges = self.neighbours;
67        adjacent_edges.map(|nbr| (nbr / N, nbr % N))
68    }
69}
70
71/// Lorentian/causal (N-1)-simplex, built from N vertices.
72///
73/// The intended `neighbours` and `vertices` ordering is such that the vertex corresponding
74/// to a neighbouring simplex is the vertex opposite the simplex. Additionally the overall
75/// ordering is supposed to be consistent, such that the simplices are *oriented*.
76#[serde_as]
77#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Serialize, Deserialize)]
78pub struct CausalSimplex<K, const N: usize> {
79    pub t: u16,
80    pub kind: K,
81    #[serde_as(as = "[_; N]")]
82    pub neighbours: [usize; N],
83    #[serde_as(as = "[_; N]")]
84    pub vertices: [usize; N],
85}
86
87impl<K: Default, const N: usize> Default for CausalSimplex<K, N> {
88    fn default() -> Self {
89        Self {
90            t: u16::default(),
91            kind: K::default(),
92            neighbours: [0; N],
93            vertices: [0; N],
94        }
95    }
96}
97
98/// Methods for working with [`CausalSimplex<N>`] when `neighbours` are encoded
99/// by `half-faces`, see [`Triangulation`](super::triangulation::Triangulation)
100impl<K, const N: usize> CausalSimplex<K, N> {
101    /// Returns the label of the neighbouring (N-1)-simplex on the `index` side.
102    ///
103    /// # Panic
104    /// This method will panic `index >= N`
105    #[inline]
106    pub fn get_neighbour(&self, index: usize) -> usize {
107        let nbr = self.neighbours[index];
108        nbr / N
109    }
110
111    /// Returns the label and backindex of neighbouring (N-1)-simplex on the `index` side.
112    ///
113    /// # Panic
114    /// This method will panic `index >= N`
115    #[inline]
116    pub fn get_neighbour_backindex(&self, index: usize) -> (usize, usize) {
117        let nbr = self.neighbours[index];
118        (nbr / N, nbr % N)
119    }
120
121    /// Returns the neighbouring (N-1)-simplex labels of the [`Simplex`]
122    #[inline]
123    pub fn get_neighbours(&self) -> [usize; N] {
124        let adjacent_edges = self.neighbours;
125        adjacent_edges.map(|nbr| nbr / N)
126    }
127
128    /// Returns the neighbouring (N-1)-simplex labels and back indices in tuple pairs
129    #[inline]
130    pub fn get_neighbours_backindices(&self) -> [(usize, usize); N] {
131        let adjacent_edges = self.neighbours;
132        adjacent_edges.map(|nbr| (nbr / N, nbr % N))
133    }
134}
135
136/// Kind/type/variant specification of a triangle/2-simplex.
137///
138/// Variant `S21` refers to a triangle with 2 vertices at the previous time-slice
139/// to the other one, where variant `S12` is the opposite.
140#[derive(
141    Debug, Clone, Copy, PartialEq, Eq, Hash, Default, PartialOrd, Ord, Serialize, Deserialize,
142)]
143pub enum CausalSimplex2Kind {
144    #[default]
145    S21 = 0,
146    S12 = 1,
147}
148pub type CausalTriangle = CausalSimplex<CausalSimplex2Kind, 3>;
149pub type CausalSimplex2 = CausalSimplex<CausalSimplex2Kind, 3>;
150
151/// Kind/type/variant specification of a tetrahedron/3-simplex.
152///
153/// Variant `S31` refers to a tetrahedron with 3 vertices at the previous time-slice
154/// to 1 vertex in the next, and variant `S13` is the opposite.
155/// Variant `S22` refers to a tetrahedron with 2 vertices at both time-slices.
156#[derive(
157    Debug, Clone, Copy, PartialEq, Eq, Hash, Default, PartialOrd, Ord, Serialize, Deserialize,
158)]
159pub enum CausalSimplex3Kind {
160    #[default]
161    S31 = 0,
162    S22 = 1,
163    S13 = 2,
164}
165pub type CausalTetrahedron = CausalSimplex<CausalSimplex3Kind, 4>;
166pub type CausalSimplex3 = CausalSimplex<CausalSimplex3Kind, 4>;
167
168/// Kind/type/variant specification of a 4-simplex.
169///
170/// Variants `S41` and `S14` refer to a 4-simplex with 4 vertices at the previous time-slice
171/// to 1 vertex in the next and the reverse respectivally, i.e tertahedron to point.
172/// Variants `S32` and `S23` refer to a 4-simplex with 3 vertices at the previous time-slice
173/// to 2 vertices in the next and the reverse respectivally, i.e. triangle to edge.
174#[derive(
175    Debug, Clone, Copy, PartialEq, Eq, Hash, Default, PartialOrd, Ord, Serialize, Deserialize,
176)]
177pub enum CausalSimplex4Kind {
178    #[default]
179    S41 = 0,
180    S32 = 1,
181    S23 = 2,
182    S14 = 3,
183}
184pub type CausalSimplex4 = CausalSimplex<CausalSimplex4Kind, 5>;
185
186impl<K, const N: usize> From<CausalSimplex<K, N>> for Simplex<N> {
187    /// Convert [`CausalSimplex`] to [`Simplex`] by simply dropping `time` and `kind` information.
188    fn from(value: CausalSimplex<K, N>) -> Self {
189        Self {
190            neighbours: value.neighbours,
191            vertices: value.vertices,
192        }
193    }
194}
195
196#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Default, PartialOrd, Ord)]
197pub enum CausalOrientation {
198    Past = 0,
199    #[default]
200    Space = 1,
201    Future = 2,
202}
203
204pub trait CausalSimplexKind {
205    /// Return spacetime orientation representing the orientation of the simplex
206    ///
207    /// Consider a path starting inside the simplex moving through the simplex faces:
208    /// The orientation is `Past` if there is a face moving through which goes towards the past;
209    /// it is `Future` if there is a face moving through which goes towards the future;
210    /// and it is `Space` if moving through all faces moves in spacelike direction.
211    /// Note: it is indeed only possible for CDT simplices to be only one of these 3.
212    fn orientation(&self) -> CausalOrientation;
213    fn signature(&self) -> (u8, u8);
214}
215
216impl CausalSimplexKind for CausalSimplex2Kind {
217    fn orientation(&self) -> CausalOrientation {
218        match self {
219            CausalSimplex2Kind::S21 => CausalOrientation::Past,
220            CausalSimplex2Kind::S12 => CausalOrientation::Future,
221        }
222    }
223
224    fn signature(&self) -> (u8, u8) {
225        match self {
226            CausalSimplex2Kind::S21 => (2, 1),
227            CausalSimplex2Kind::S12 => (1, 2),
228        }
229    }
230}
231
232impl CausalSimplexKind for CausalSimplex3Kind {
233    fn orientation(&self) -> CausalOrientation {
234        match self {
235            CausalSimplex3Kind::S31 => CausalOrientation::Past,
236            CausalSimplex3Kind::S22 => CausalOrientation::Space,
237            CausalSimplex3Kind::S13 => CausalOrientation::Future,
238        }
239    }
240
241    fn signature(&self) -> (u8, u8) {
242        match self {
243            CausalSimplex3Kind::S31 => (3, 1),
244            CausalSimplex3Kind::S22 => (2, 2),
245            CausalSimplex3Kind::S13 => (1, 3),
246        }
247    }
248}
249
250impl CausalSimplexKind for CausalSimplex4Kind {
251    fn orientation(&self) -> CausalOrientation {
252        match self {
253            CausalSimplex4Kind::S41 => CausalOrientation::Past,
254            CausalSimplex4Kind::S32 => CausalOrientation::Space,
255            CausalSimplex4Kind::S23 => CausalOrientation::Space,
256            CausalSimplex4Kind::S14 => CausalOrientation::Future,
257        }
258    }
259    fn signature(&self) -> (u8, u8) {
260        match self {
261            CausalSimplex4Kind::S41 => (4, 1),
262            CausalSimplex4Kind::S32 => (3, 2),
263            CausalSimplex4Kind::S23 => (2, 3),
264            CausalSimplex4Kind::S14 => (1, 4),
265        }
266    }
267}
268
269/// Orders the vertices of the simplex from small to big to give a canonical ordering for comparison
270pub(crate) fn canonical_simplex<T: Ord, const N: usize>(mut simplex: [T; N]) -> [T; N] {
271    simplex.sort();
272    simplex
273}
274
275/// Orders the vertices of the simplex from small to big to give a canonical ordering for comparison,
276/// additionally giving the parity of the permulation to have ordering (true-even, false-odd)
277pub(crate) fn canonical_oriented_simplex<T: Ord, const N: usize>(
278    mut simplex: [T; N],
279) -> (bool, [T; N]) {
280    // Use bubble sort
281    let mut parity = true;
282    loop {
283        let mut unchanged = true;
284        for i in 0..(N - 1) {
285            if simplex[i] > simplex[i + 1] {
286                parity = !parity;
287                simplex.swap(i, i + 1);
288                unchanged = false;
289            }
290        }
291        if unchanged {
292            break;
293        }
294    }
295    (parity, simplex)
296}
297
298#[cfg(test)]
299mod tests {
300    use crate::triangulation::simplices::canonical_oriented_simplex;
301
302    #[test]
303    fn test_canonical_oriented_simplex() {
304        assert_eq!(
305            canonical_oriented_simplex([0, 1, 2, 3, 4]),
306            (true, [0, 1, 2, 3, 4])
307        );
308        assert_eq!(
309            canonical_oriented_simplex([1, 0, 2, 3, 4]),
310            (false, [0, 1, 2, 3, 4])
311        );
312        assert_eq!(
313            canonical_oriented_simplex([3, 1, 2, 0, 4]),
314            (false, [0, 1, 2, 3, 4])
315        );
316        assert_eq!(
317            canonical_oriented_simplex([4, 1, 2, 3, 0]),
318            (false, [0, 1, 2, 3, 4])
319        );
320        assert_eq!(
321            canonical_oriented_simplex([1, 3, 0, 4, 2]),
322            (true, [0, 1, 2, 3, 4])
323        );
324    }
325}