ndelement/
orientation.rs

1//! Cell orientation.
2
3use crate::{reference_cell, types::ReferenceCellType};
4use itertools::izip;
5
6/// Compute a 32-bit integer that encodes the orientation differences between the reference cell and a cell with vertices numbered as input
7///
8/// From right to left, the bits of this function's output encode:
9///  - 1 bit for each edge - set to 1 if the edge needs reversing
10///  - 3 bits for each face - two rightmost bits encode number of rotations that need to be applied to face; left bit encodes whether it needs reflecting
11pub fn compute_orientation(entity_type: ReferenceCellType, vertices: &[usize]) -> i32 {
12    let mut orientation = 0;
13    let mut n = 0;
14    let dim = reference_cell::dim(entity_type);
15    if dim > 1 {
16        for v in reference_cell::edges(entity_type) {
17            if vertices[v[0]] > vertices[v[1]] {
18                orientation += 1 << n;
19            }
20            n += 1;
21        }
22    }
23    if dim > 2 {
24        for (t, v) in izip!(
25            &reference_cell::entity_types(entity_type)[2],
26            reference_cell::faces(entity_type)
27        ) {
28            match t {
29                ReferenceCellType::Triangle => {
30                    let m = v
31                        .iter()
32                        .map(|i| vertices[*i])
33                        .enumerate()
34                        .min_by(|(_, a), (_, b)| a.partial_cmp(b).unwrap())
35                        .map(|(index, _)| index)
36                        .unwrap() as i32;
37                    orientation += m << n;
38                    n += 2;
39                    let next = vertices[v[(m as usize + 1) % 3]];
40                    let prev = vertices[v[if m == 0 { 2 } else { m as usize - 1 }]];
41                    if next > prev {
42                        orientation += 1 << n;
43                    }
44                    n += 1;
45                }
46                ReferenceCellType::Quadrilateral => {
47                    let m = v
48                        .iter()
49                        .map(|i| vertices[*i])
50                        .enumerate()
51                        .min_by(|(_, a), (_, b)| a.partial_cmp(b).unwrap())
52                        .map(|(index, _)| index)
53                        .unwrap() as i32;
54                    if m == 2 {
55                        orientation += 3 << n;
56                    } else if m == 3 {
57                        orientation += 2 << n;
58                    } else {
59                        orientation += m << n;
60                    }
61                    n += 2;
62                    let next = vertices[v[[1, 3, 0, 2][m as usize]]];
63                    let prev = vertices[v[[2, 0, 3, 1][m as usize]]];
64                    if next > prev {
65                        orientation += 1 << n;
66                    }
67                    n += 1;
68                }
69                _ => {
70                    panic!("Unsupported face type: {t:?}");
71                }
72            }
73        }
74    }
75    orientation
76}
77
78#[cfg(test)]
79mod test {
80    use super::*;
81
82    #[test]
83    fn test_compute_orientations_interval() {
84        for (v, o) in [([0, 1], 0), ([1, 0], 0)] {
85            assert_eq!(compute_orientation(ReferenceCellType::Interval, &v), o);
86        }
87    }
88
89    #[test]
90    fn test_compute_orientations_triangle() {
91        for (v, o) in [
92            ([0, 1, 2], 0),
93            ([0, 2, 1], 1),
94            ([1, 0, 2], 4),
95            ([1, 2, 0], 3),
96            ([2, 0, 1], 6),
97            ([2, 1, 0], 7),
98        ] {
99            assert_eq!(compute_orientation(ReferenceCellType::Triangle, &v), o);
100        }
101    }
102
103    #[test]
104    fn test_compute_orientations_quadrilateral() {
105        for (v, o) in [
106            ([0, 1, 2, 3], 0),
107            ([0, 1, 3, 2], 8),
108            ([0, 2, 1, 3], 0),
109            ([0, 2, 3, 1], 12),
110            ([0, 3, 1, 2], 4),
111            ([0, 3, 2, 1], 12),
112            ([1, 0, 2, 3], 1),
113            ([1, 0, 3, 2], 9),
114            ([1, 2, 0, 3], 2),
115            ([1, 2, 3, 0], 12),
116            ([1, 3, 0, 2], 6),
117            ([1, 3, 2, 0], 12),
118        ] {
119            assert_eq!(compute_orientation(ReferenceCellType::Quadrilateral, &v), o);
120        }
121    }
122
123    #[test]
124    fn test_compute_orientations_tetrahedron() {
125        for (v, o) in [
126            ([0, 1, 2, 3], 0),
127            ([0, 3, 2, 1], 149895), // 100 100 100 110 000111
128            ([1, 2, 0, 3], 68436),  // 010 000 101 101 010100
129            ([1, 3, 2, 0], 140687), // 100 010 010 110 001111
130        ] {
131            assert_eq!(compute_orientation(ReferenceCellType::Tetrahedron, &v), o);
132        }
133    }
134
135    #[test]
136    fn test_compute_orientations_hexahedron() {
137        for (v, o) in [
138            ([0, 1, 2, 3, 4, 5, 6, 7], 0),
139            ([2, 5, 6, 1, 7, 0, 4, 3], 165899128), // 001 001 111 000 110 110 101101111000
140            ([4, 1, 5, 0, 2, 3, 6, 7], 95215661),  // 000 101 101 011 001 110 000000101101
141            ([6, 5, 1, 7, 2, 4, 3, 0], 306691223), // 010 010 010 001 111 011 110010010111
142        ] {
143            assert_eq!(compute_orientation(ReferenceCellType::Hexahedron, &v), o);
144        }
145    }
146}
ndelement/orientation.rs

ndelement/
orientation.rs
