ndelement/
reference_cell.rs

1//! Cell definitions.
2//!
3//! The reference cells used by ndelement are defined by:
4//! - The reference interval is \(\left\{x\,\middle|\,0\leqslant x\leqslant 1\right\}\).
5//! - The reference triangle is \(\left\{(x,y)\,\middle|\,0\leqslant x,\,0\leqslant y,\,x+y\leqslant1\right\}\).
6//! - The reference quadrilateral is \(\left\{(x,y)\,\middle|\,0\leqslant x\leqslant1,\,0\leqslant y\leqslant1\right\}\).
7//! - The reference tetrahedron is \(\left\{(x,y,z)\,\middle|\,0\leqslant x,\,0\leqslant y,\,0\leqslant z,\,x+y+z\leqslant1\right\}\).
8//! - The reference hexahedron is \(\left\{(x,y,z)\,\middle|\,0\leqslant x\leqslant1,\,0\leqslant y\leqslant1,\,0\leqslant z\leqslant1\right\}\).
9//! - The reference prism is \(\left\{(x,y,z)\,\middle|\,0\leqslant x,\,0\leqslant y,\,x+y\leqslant1,\,0\leqslant z\leqslant1\right\}\).
10//! - The reference pyramid is \(\left\{(x,y,z)\,\middle|\,0\leqslant x,\,0\leqslant y,\,0\leqslant z\leqslant1,\,x+z\leqslant1,\,y+z\leqslant1\right\}\).
11
12use crate::types::ReferenceCellType;
13use rlst::RlstScalar;
14
15/// The topological dimension of the cell
16pub fn dim(cell: ReferenceCellType) -> usize {
17    match cell {
18        ReferenceCellType::Point => 0,
19        ReferenceCellType::Interval => 1,
20        ReferenceCellType::Triangle => 2,
21        ReferenceCellType::Quadrilateral => 2,
22        ReferenceCellType::Tetrahedron => 3,
23        ReferenceCellType::Hexahedron => 3,
24        ReferenceCellType::Prism => 3,
25        ReferenceCellType::Pyramid => 3,
26    }
27}
28/// Is the cell a simplex?
29pub fn is_simplex(cell: ReferenceCellType) -> bool {
30    match cell {
31        ReferenceCellType::Point => true,
32        ReferenceCellType::Interval => true,
33        ReferenceCellType::Triangle => true,
34        ReferenceCellType::Quadrilateral => false,
35        ReferenceCellType::Tetrahedron => true,
36        ReferenceCellType::Hexahedron => false,
37        ReferenceCellType::Prism => false,
38        ReferenceCellType::Pyramid => false,
39    }
40}
41
42/// The vertices of the reference cell
43pub fn vertices<T: RlstScalar>(cell: ReferenceCellType) -> Vec<Vec<T>> {
44    let zero = T::from(0.0).unwrap();
45    let one = T::from(1.0).unwrap();
46    match cell {
47        ReferenceCellType::Point => vec![],
48        ReferenceCellType::Interval => vec![vec![zero], vec![one]],
49        ReferenceCellType::Triangle => vec![vec![zero, zero], vec![one, zero], vec![zero, one]],
50        ReferenceCellType::Quadrilateral => vec![
51            vec![zero, zero],
52            vec![one, zero],
53            vec![zero, one],
54            vec![one, one],
55        ],
56        ReferenceCellType::Tetrahedron => vec![
57            vec![zero, zero, zero],
58            vec![one, zero, zero],
59            vec![zero, one, zero],
60            vec![zero, zero, one],
61        ],
62        ReferenceCellType::Hexahedron => vec![
63            vec![zero, zero, zero],
64            vec![one, zero, zero],
65            vec![zero, one, zero],
66            vec![one, one, zero],
67            vec![zero, zero, one],
68            vec![one, zero, one],
69            vec![zero, one, one],
70            vec![one, one, one],
71        ],
72        ReferenceCellType::Prism => vec![
73            vec![zero, zero, zero],
74            vec![one, zero, zero],
75            vec![zero, one, zero],
76            vec![zero, zero, one],
77            vec![one, zero, one],
78            vec![zero, one, one],
79        ],
80        ReferenceCellType::Pyramid => vec![
81            vec![zero, zero, zero],
82            vec![one, zero, zero],
83            vec![zero, one, zero],
84            vec![one, one, zero],
85            vec![zero, zero, one],
86        ],
87    }
88}
89
90/// The midpoint of the cell
91pub fn midpoint<T: RlstScalar>(cell: ReferenceCellType) -> Vec<T> {
92    let half = T::from(0.5).unwrap();
93    let third = T::from(1.0).unwrap() / T::from(3.0).unwrap();
94    match cell {
95        ReferenceCellType::Point => vec![],
96        ReferenceCellType::Interval => vec![half],
97        ReferenceCellType::Triangle => vec![third; 2],
98        ReferenceCellType::Quadrilateral => vec![half; 2],
99        ReferenceCellType::Tetrahedron => vec![T::from(1.0).unwrap() / T::from(4.0).unwrap(); 3],
100        ReferenceCellType::Hexahedron => vec![half; 3],
101        ReferenceCellType::Prism => vec![third, third, half],
102        ReferenceCellType::Pyramid => vec![
103            T::from(0.4).unwrap(),
104            T::from(0.4).unwrap(),
105            T::from(0.2).unwrap(),
106        ],
107    }
108}
109
110/// The edges (dimension 1 entities) of the reference cell
111pub fn edges(cell: ReferenceCellType) -> Vec<Vec<usize>> {
112    match cell {
113        ReferenceCellType::Point => vec![],
114        ReferenceCellType::Interval => vec![vec![0, 1]],
115        ReferenceCellType::Triangle => vec![vec![1, 2], vec![0, 2], vec![0, 1]],
116        ReferenceCellType::Quadrilateral => vec![vec![0, 1], vec![0, 2], vec![1, 3], vec![2, 3]],
117        ReferenceCellType::Tetrahedron => vec![
118            vec![2, 3],
119            vec![1, 3],
120            vec![1, 2],
121            vec![0, 3],
122            vec![0, 2],
123            vec![0, 1],
124        ],
125        ReferenceCellType::Hexahedron => vec![
126            vec![0, 1],
127            vec![0, 2],
128            vec![0, 4],
129            vec![1, 3],
130            vec![1, 5],
131            vec![2, 3],
132            vec![2, 6],
133            vec![3, 7],
134            vec![4, 5],
135            vec![4, 6],
136            vec![5, 7],
137            vec![6, 7],
138        ],
139        ReferenceCellType::Prism => vec![
140            vec![0, 1],
141            vec![0, 2],
142            vec![0, 3],
143            vec![1, 2],
144            vec![1, 4],
145            vec![2, 5],
146            vec![3, 4],
147            vec![3, 5],
148            vec![4, 5],
149        ],
150        ReferenceCellType::Pyramid => vec![
151            vec![0, 1],
152            vec![0, 2],
153            vec![0, 4],
154            vec![1, 3],
155            vec![1, 4],
156            vec![2, 3],
157            vec![2, 4],
158            vec![3, 4],
159        ],
160    }
161}
162
163/// The faces (dimension 2 entities) of the reference cell
164pub fn faces(cell: ReferenceCellType) -> Vec<Vec<usize>> {
165    match cell {
166        ReferenceCellType::Point => vec![],
167        ReferenceCellType::Interval => vec![],
168        ReferenceCellType::Triangle => vec![vec![0, 1, 2]],
169        ReferenceCellType::Quadrilateral => vec![vec![0, 1, 2, 3]],
170        ReferenceCellType::Tetrahedron => {
171            vec![vec![1, 2, 3], vec![0, 2, 3], vec![0, 1, 3], vec![0, 1, 2]]
172        }
173        ReferenceCellType::Hexahedron => vec![
174            vec![0, 1, 2, 3],
175            vec![0, 1, 4, 5],
176            vec![0, 2, 4, 6],
177            vec![1, 3, 5, 7],
178            vec![2, 3, 6, 7],
179            vec![4, 5, 6, 7],
180        ],
181        ReferenceCellType::Prism => vec![
182            vec![0, 1, 2],
183            vec![0, 1, 3, 4],
184            vec![0, 2, 3, 5],
185            vec![1, 2, 4, 5],
186            vec![3, 4, 5],
187        ],
188        ReferenceCellType::Pyramid => vec![
189            vec![0, 1, 2, 3],
190            vec![0, 1, 4],
191            vec![0, 2, 4],
192            vec![1, 3, 4],
193            vec![2, 3, 4],
194        ],
195    }
196}
197
198/// The volumes (dimension 3 entities) of the reference cell
199pub fn volumes(cell: ReferenceCellType) -> Vec<Vec<usize>> {
200    match cell {
201        ReferenceCellType::Point => vec![],
202        ReferenceCellType::Interval => vec![],
203        ReferenceCellType::Triangle => vec![],
204        ReferenceCellType::Quadrilateral => vec![],
205        ReferenceCellType::Tetrahedron => vec![vec![0, 1, 2, 3]],
206        ReferenceCellType::Hexahedron => vec![vec![0, 1, 2, 3, 4, 5, 6, 7]],
207        ReferenceCellType::Prism => vec![vec![0, 1, 2, 3, 4, 5]],
208        ReferenceCellType::Pyramid => vec![vec![0, 1, 2, 3, 4]],
209    }
210}
211
212/// The facets (dimension dim-1 entities) of the reference cell
213pub fn facets(cell: ReferenceCellType) -> Vec<Vec<usize>> {
214    match dim(cell) {
215        0 => vec![],
216        1 => volumes(cell)[0]
217            .iter()
218            .map(|x| vec![*x])
219            .collect::<Vec<_>>(),
220        2 => edges(cell),
221        3 => faces(cell),
222        _ => {
223            panic!("Invalid dimension");
224        }
225    }
226}
227
228/// The ridges (dimension dim-2 entities) of the reference cell
229pub fn ridges(cell: ReferenceCellType) -> Vec<Vec<usize>> {
230    match dim(cell) {
231        0 => vec![],
232        1 => vec![],
233        2 => volumes(cell)[0]
234            .iter()
235            .map(|x| vec![*x])
236            .collect::<Vec<_>>(),
237        3 => edges(cell),
238        _ => {
239            panic!("Invalid dimension");
240        }
241    }
242}
243
244/// The peaks (dimension dim-3 entities) of the reference cell
245pub fn peaks(cell: ReferenceCellType) -> Vec<Vec<usize>> {
246    match dim(cell) {
247        0 => vec![],
248        1 => vec![],
249        2 => vec![],
250        3 => volumes(cell)[0]
251            .iter()
252            .map(|x| vec![*x])
253            .collect::<Vec<_>>(),
254        _ => {
255            panic!("Invalid dimension");
256        }
257    }
258}
259
260/// The normals to the facets of the reference cell. The length of each normal is the volume of the facet
261pub fn facet_normals<T: RlstScalar>(cell: ReferenceCellType) -> Vec<Vec<T>> {
262    let zero = T::from(0.0).unwrap();
263    let one = T::from(1.0).unwrap();
264    match cell {
265        ReferenceCellType::Point => vec![],
266        ReferenceCellType::Interval => vec![vec![one], vec![one]],
267        ReferenceCellType::Triangle => vec![vec![-one, -one], vec![-one, zero], vec![zero, one]],
268        ReferenceCellType::Quadrilateral => vec![
269            vec![zero, one],
270            vec![-one, zero],
271            vec![-one, zero],
272            vec![zero, one],
273        ],
274        ReferenceCellType::Tetrahedron => vec![
275            vec![one, one, one],
276            vec![one, zero, zero],
277            vec![zero, -one, zero],
278            vec![zero, zero, one],
279        ],
280        ReferenceCellType::Hexahedron => vec![
281            vec![zero, zero, one],
282            vec![zero, -one, zero],
283            vec![one, zero, zero],
284            vec![one, zero, zero],
285            vec![zero, -one, zero],
286            vec![zero, zero, one],
287        ],
288        ReferenceCellType::Prism => vec![
289            vec![zero, zero, one],
290            vec![zero, -one, zero],
291            vec![one, zero, zero],
292            vec![one, one, zero],
293            vec![zero, zero, one],
294        ],
295        ReferenceCellType::Pyramid => vec![
296            vec![zero, zero, one],
297            vec![zero, -one, zero],
298            vec![one, zero, zero],
299            vec![one, zero, one],
300            vec![zero, -one, -one],
301        ],
302    }
303}
304
305/// The unit normals to the facets of the reference cell
306pub fn facet_unit_normals<T: RlstScalar>(cell: ReferenceCellType) -> Vec<Vec<T>> {
307    let mut normals = facet_normals::<T>(cell);
308    for normal in normals.iter_mut() {
309        let size = normal.iter().map(|x| x.powi(2)).sum::<T>().sqrt();
310        for value in normal.iter_mut() {
311            *value /= size;
312        }
313    }
314    normals
315}
316
317/// The types of the subentities of the reference cell
318pub fn entity_types(cell: ReferenceCellType) -> Vec<Vec<ReferenceCellType>> {
319    match cell {
320        ReferenceCellType::Point => vec![vec![ReferenceCellType::Point], vec![], vec![], vec![]],
321        ReferenceCellType::Interval => vec![
322            vec![ReferenceCellType::Point; 2],
323            vec![ReferenceCellType::Interval],
324            vec![],
325            vec![],
326        ],
327        ReferenceCellType::Triangle => vec![
328            vec![ReferenceCellType::Point; 3],
329            vec![ReferenceCellType::Interval; 3],
330            vec![ReferenceCellType::Triangle],
331            vec![],
332        ],
333        ReferenceCellType::Quadrilateral => vec![
334            vec![ReferenceCellType::Point; 4],
335            vec![ReferenceCellType::Interval; 4],
336            vec![ReferenceCellType::Quadrilateral],
337            vec![],
338        ],
339        ReferenceCellType::Tetrahedron => vec![
340            vec![ReferenceCellType::Point; 4],
341            vec![ReferenceCellType::Interval; 6],
342            vec![ReferenceCellType::Triangle; 4],
343            vec![ReferenceCellType::Tetrahedron],
344        ],
345        ReferenceCellType::Hexahedron => vec![
346            vec![ReferenceCellType::Point; 8],
347            vec![ReferenceCellType::Interval; 12],
348            vec![ReferenceCellType::Quadrilateral; 6],
349            vec![ReferenceCellType::Hexahedron],
350        ],
351        ReferenceCellType::Prism => vec![
352            vec![ReferenceCellType::Point; 6],
353            vec![ReferenceCellType::Interval; 9],
354            vec![
355                ReferenceCellType::Triangle,
356                ReferenceCellType::Quadrilateral,
357                ReferenceCellType::Quadrilateral,
358                ReferenceCellType::Quadrilateral,
359                ReferenceCellType::Triangle,
360            ],
361            vec![ReferenceCellType::Prism],
362        ],
363        ReferenceCellType::Pyramid => vec![
364            vec![ReferenceCellType::Point; 5],
365            vec![ReferenceCellType::Interval; 8],
366            vec![
367                ReferenceCellType::Quadrilateral,
368                ReferenceCellType::Triangle,
369                ReferenceCellType::Triangle,
370                ReferenceCellType::Triangle,
371                ReferenceCellType::Triangle,
372            ],
373            vec![ReferenceCellType::Pyramid],
374        ],
375    }
376}
377
378/// The number of subentities of each dimension
379pub fn entity_counts(cell: ReferenceCellType) -> Vec<usize> {
380    match cell {
381        ReferenceCellType::Point => vec![1, 0, 0, 0],
382        ReferenceCellType::Interval => vec![2, 1, 0, 0],
383        ReferenceCellType::Triangle => vec![3, 3, 1, 0],
384        ReferenceCellType::Quadrilateral => vec![4, 4, 1, 0],
385        ReferenceCellType::Tetrahedron => vec![4, 6, 4, 1],
386        ReferenceCellType::Hexahedron => vec![8, 12, 6, 1],
387        ReferenceCellType::Prism => vec![6, 9, 5, 1],
388        ReferenceCellType::Pyramid => vec![5, 8, 5, 1],
389    }
390}
391
392/// The connectivity of the reference cell
393///
394/// The indices of the result are \[dim0\]\[index0\]\[dim1\]\[\]
395pub fn connectivity(cell: ReferenceCellType) -> Vec<Vec<Vec<Vec<usize>>>> {
396    match cell {
397        ReferenceCellType::Point => vec![vec![vec![vec![0]]]],
398        ReferenceCellType::Interval => vec![
399            vec![vec![vec![0], vec![0]], vec![vec![1], vec![0]]],
400            vec![vec![vec![0, 1], vec![0]]],
401        ],
402        ReferenceCellType::Triangle => vec![
403            vec![
404                vec![vec![0], vec![1, 2], vec![0]],
405                vec![vec![1], vec![0, 2], vec![0]],
406                vec![vec![2], vec![0, 1], vec![0]],
407            ],
408            vec![
409                vec![vec![1, 2], vec![0], vec![0]],
410                vec![vec![0, 2], vec![1], vec![0]],
411                vec![vec![0, 1], vec![2], vec![0]],
412            ],
413            vec![vec![vec![0, 1, 2], vec![0, 1, 2], vec![0]]],
414        ],
415        ReferenceCellType::Quadrilateral => vec![
416            vec![
417                vec![vec![0], vec![0, 1], vec![0]],
418                vec![vec![1], vec![0, 2], vec![0]],
419                vec![vec![2], vec![1, 3], vec![0]],
420                vec![vec![3], vec![2, 3], vec![0]],
421            ],
422            vec![
423                vec![vec![0, 1], vec![0], vec![0]],
424                vec![vec![0, 2], vec![1], vec![0]],
425                vec![vec![1, 3], vec![2], vec![0]],
426                vec![vec![2, 3], vec![3], vec![0]],
427            ],
428            vec![vec![vec![0, 1, 2, 3], vec![0, 1, 2, 3], vec![0]]],
429        ],
430        ReferenceCellType::Tetrahedron => vec![
431            vec![
432                vec![vec![0], vec![3, 4, 5], vec![1, 2, 3], vec![0]],
433                vec![vec![1], vec![1, 2, 5], vec![0, 2, 3], vec![0]],
434                vec![vec![2], vec![0, 2, 4], vec![0, 1, 3], vec![0]],
435                vec![vec![3], vec![0, 1, 3], vec![0, 1, 2], vec![0]],
436            ],
437            vec![
438                vec![vec![2, 3], vec![0], vec![0, 1], vec![0]],
439                vec![vec![1, 3], vec![1], vec![0, 2], vec![0]],
440                vec![vec![1, 2], vec![2], vec![0, 3], vec![0]],
441                vec![vec![0, 3], vec![3], vec![1, 2], vec![0]],
442                vec![vec![0, 2], vec![4], vec![1, 3], vec![0]],
443                vec![vec![0, 1], vec![5], vec![2, 3], vec![0]],
444            ],
445            vec![
446                vec![vec![1, 2, 3], vec![0, 1, 2], vec![0], vec![0]],
447                vec![vec![0, 2, 3], vec![0, 3, 4], vec![1], vec![0]],
448                vec![vec![0, 1, 3], vec![1, 3, 5], vec![2], vec![0]],
449                vec![vec![0, 1, 2], vec![2, 4, 5], vec![3], vec![0]],
450            ],
451            vec![vec![
452                vec![0, 1, 2, 3],
453                vec![0, 1, 2, 3, 4, 5],
454                vec![0, 1, 2, 3],
455                vec![0],
456            ]],
457        ],
458        ReferenceCellType::Hexahedron => vec![
459            vec![
460                vec![vec![0], vec![0, 1, 2], vec![0, 1, 2], vec![0]],
461                vec![vec![1], vec![0, 3, 4], vec![0, 1, 3], vec![0]],
462                vec![vec![2], vec![1, 5, 6], vec![0, 2, 4], vec![0]],
463                vec![vec![3], vec![3, 5, 7], vec![0, 3, 4], vec![0]],
464                vec![vec![4], vec![2, 8, 9], vec![1, 2, 5], vec![0]],
465                vec![vec![5], vec![4, 8, 10], vec![1, 3, 5], vec![0]],
466                vec![vec![6], vec![6, 9, 11], vec![2, 4, 5], vec![0]],
467                vec![vec![7], vec![7, 10, 11], vec![3, 4, 5], vec![0]],
468            ],
469            vec![
470                vec![vec![0, 1], vec![0], vec![0, 1], vec![0]],
471                vec![vec![0, 2], vec![1], vec![0, 2], vec![0]],
472                vec![vec![0, 4], vec![2], vec![1, 2], vec![0]],
473                vec![vec![1, 3], vec![3], vec![0, 3], vec![0]],
474                vec![vec![1, 5], vec![4], vec![1, 3], vec![0]],
475                vec![vec![2, 3], vec![5], vec![0, 4], vec![0]],
476                vec![vec![2, 6], vec![6], vec![2, 4], vec![0]],
477                vec![vec![3, 7], vec![7], vec![3, 4], vec![0]],
478                vec![vec![4, 5], vec![8], vec![1, 5], vec![0]],
479                vec![vec![4, 6], vec![9], vec![2, 5], vec![0]],
480                vec![vec![5, 7], vec![10], vec![3, 5], vec![0]],
481                vec![vec![6, 7], vec![11], vec![4, 5], vec![0]],
482            ],
483            vec![
484                vec![vec![0, 1, 2, 3], vec![0, 1, 3, 5], vec![0], vec![0]],
485                vec![vec![0, 1, 4, 5], vec![0, 2, 4, 8], vec![1], vec![0]],
486                vec![vec![0, 2, 4, 6], vec![1, 2, 6, 9], vec![2], vec![0]],
487                vec![vec![1, 3, 5, 7], vec![3, 4, 7, 10], vec![3], vec![0]],
488                vec![vec![2, 3, 6, 7], vec![5, 6, 7, 11], vec![4], vec![0]],
489                vec![vec![4, 5, 6, 7], vec![8, 9, 10, 11], vec![5], vec![0]],
490            ],
491            vec![vec![
492                vec![0, 1, 2, 3, 4, 5, 6, 7],
493                vec![0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11],
494                vec![0, 1, 2, 3, 4, 5],
495                vec![0],
496            ]],
497        ],
498        ReferenceCellType::Prism => vec![
499            vec![
500                vec![vec![0], vec![0, 1, 2], vec![0, 1, 2], vec![0]],
501                vec![vec![1], vec![0, 3, 4], vec![0, 1, 3], vec![0]],
502                vec![vec![2], vec![1, 3, 5], vec![0, 2, 3], vec![0]],
503                vec![vec![3], vec![2, 6, 7], vec![1, 2, 4], vec![0]],
504                vec![vec![4], vec![4, 6, 8], vec![1, 3, 4], vec![0]],
505                vec![vec![5], vec![5, 7, 8], vec![2, 3, 4], vec![0]],
506            ],
507            vec![
508                vec![vec![0, 1], vec![0], vec![0, 1], vec![0]],
509                vec![vec![0, 2], vec![1], vec![0, 2], vec![0]],
510                vec![vec![0, 3], vec![2], vec![1, 2], vec![0]],
511                vec![vec![1, 2], vec![3], vec![0, 3], vec![0]],
512                vec![vec![1, 4], vec![4], vec![1, 3], vec![0]],
513                vec![vec![2, 5], vec![5], vec![2, 3], vec![0]],
514                vec![vec![3, 4], vec![6], vec![1, 4], vec![0]],
515                vec![vec![3, 5], vec![7], vec![2, 4], vec![0]],
516                vec![vec![4, 5], vec![8], vec![3, 4], vec![0]],
517            ],
518            vec![
519                vec![vec![0, 1, 2], vec![0, 1, 3], vec![0], vec![0]],
520                vec![vec![0, 1, 3, 4], vec![0, 2, 4, 6], vec![1], vec![0]],
521                vec![vec![0, 2, 3, 5], vec![1, 2, 5, 7], vec![2], vec![0]],
522                vec![vec![1, 2, 4, 5], vec![3, 4, 5, 8], vec![3], vec![0]],
523                vec![vec![3, 4, 5], vec![6, 7, 8], vec![4], vec![0]],
524            ],
525            vec![vec![
526                vec![0, 1, 2, 3, 4, 5],
527                vec![0, 1, 2, 3, 4, 5, 6, 7, 8],
528                vec![0, 1, 2, 3, 4],
529                vec![0],
530            ]],
531        ],
532        ReferenceCellType::Pyramid => vec![
533            vec![
534                vec![vec![0], vec![0, 1, 2], vec![0, 1, 2], vec![0]],
535                vec![vec![1], vec![0, 3, 4], vec![0, 1, 3], vec![0]],
536                vec![vec![2], vec![1, 5, 6], vec![0, 2, 4], vec![0]],
537                vec![vec![3], vec![3, 5, 7], vec![0, 3, 4], vec![0]],
538                vec![vec![4], vec![2, 4, 6, 7], vec![1, 2, 3, 4], vec![0]],
539            ],
540            vec![
541                vec![vec![0, 1], vec![0], vec![0, 1], vec![0]],
542                vec![vec![0, 2], vec![1], vec![0, 2], vec![0]],
543                vec![vec![0, 4], vec![2], vec![1, 2], vec![0]],
544                vec![vec![1, 3], vec![3], vec![0, 3], vec![0]],
545                vec![vec![1, 4], vec![4], vec![1, 3], vec![0]],
546                vec![vec![2, 3], vec![5], vec![0, 4], vec![0]],
547                vec![vec![2, 4], vec![6], vec![2, 4], vec![0]],
548                vec![vec![3, 4], vec![7], vec![3, 4], vec![0]],
549            ],
550            vec![
551                vec![vec![0, 1, 2, 3], vec![0, 1, 3, 5], vec![0], vec![0]],
552                vec![vec![0, 1, 4], vec![0, 2, 4], vec![1], vec![0]],
553                vec![vec![0, 2, 4], vec![1, 2, 6], vec![2], vec![0]],
554                vec![vec![1, 3, 4], vec![3, 4, 7], vec![3], vec![0]],
555                vec![vec![2, 3, 4], vec![5, 6, 7], vec![4], vec![0]],
556            ],
557            vec![vec![
558                vec![0, 1, 2, 3, 4],
559                vec![0, 1, 2, 3, 4, 5, 6, 7],
560                vec![0, 1, 2, 3, 4],
561                vec![0],
562            ]],
563        ],
564    }
565}
566
567#[cfg(test)]
568mod test {
569    use super::*;
570    use approx::*;
571    use itertools::izip;
572    use paste::paste;
573
574    macro_rules! test_cell {
575        ($cell:ident) => {
576            paste! {
577
578                #[test]
579                fn [<test_ $cell:lower>]() {
580                    let v = vertices::<f64>(ReferenceCellType::[<$cell>]);
581                    let d = dim(ReferenceCellType::[<$cell>]);
582                    let ec = entity_counts(ReferenceCellType::[<$cell>]);
583                    let et = entity_types(ReferenceCellType::[<$cell>]);
584                    let conn = connectivity(ReferenceCellType::[<$cell>]);
585                    for i in 0..d+1 {
586                        assert_eq!(ec[i], et[i].len());
587                        assert_eq!(ec[i], conn[i].len());
588                    }
589                    assert_eq!(ec[0], v.len());
590                    for i in v {
591                        assert_eq!(i.len(), d);
592                    }
593
594                    for v_n in 0..ec[0] {
595                        let v = &conn[0][v_n][0];
596                        assert_eq!(v, &[v_n]);
597                    }
598                    for e_n in 0..ec[1] {
599                        let vs = &conn[1][e_n][0];
600                        assert_eq!(vs, &edges(ReferenceCellType::[<$cell>])[e_n]);
601                    }
602                    for f_n in 0..ec[2] {
603                        let vs = &conn[2][f_n][0];
604                        assert_eq!(vs, &faces(ReferenceCellType::[<$cell>])[f_n]);
605                    }
606
607                    for e_dim in 0..d {
608                        for e_n in 0..ec[e_dim] {
609                            let e_vertices = &conn[e_dim][e_n][0];
610                            for c_dim in 0..d + 1 {
611                                let connectivity = &conn[e_dim][e_n][c_dim];
612                                if e_dim == c_dim {
613                                    assert_eq!(connectivity.len(), 1);
614                                    assert_eq!(connectivity[0], e_n)
615                                } else {
616                                    for c_n in connectivity {
617                                        let c_vertices = &conn[c_dim][*c_n][0];
618                                        if e_dim < c_dim {
619                                            for i in e_vertices {
620                                                assert!(c_vertices.contains(&i));
621                                            }
622                                        } else {
623                                            for i in c_vertices {
624                                                assert!(e_vertices.contains(&i));
625                                            }
626                                        }
627                                        assert!(connectivity.contains(&c_n));
628                                    }
629                                }
630                            }
631                        }
632                    }
633                }
634            }
635        };
636    }
637
638    test_cell!(Interval);
639    test_cell!(Triangle);
640    test_cell!(Quadrilateral);
641    test_cell!(Tetrahedron);
642    test_cell!(Hexahedron);
643    test_cell!(Prism);
644    test_cell!(Pyramid);
645
646    macro_rules! test_facet_normals_2d {
647        ($cell:ident) => {
648            paste! {
649
650                #[test]
651                fn [<test_ $cell:lower _facet_normals>]() {
652                    let v = vertices::<f64>(ReferenceCellType::[<$cell>]);
653                    let ns = facet_normals::<f64>(ReferenceCellType::[<$cell>]);
654                    let es = edges(ReferenceCellType::[<$cell>]);
655
656                    for (e, n) in izip!(es, ns) {
657                        let tangent = (0..2).map(|i| v[e[1]][i] - v[e[0]][i]).collect::<Vec<_>>();
658                        assert_relative_eq!(n[0], -tangent[1]);
659                        assert_relative_eq!(n[1], tangent[0]);
660                    }
661                }
662            }
663        };
664    }
665
666    fn cross<T: RlstScalar>(a: &[T], b: &[T]) -> Vec<T> {
667        assert_eq!(a.len(), 3);
668        assert_eq!(b.len(), 3);
669        vec![
670            a[1] * b[2] - a[2] * b[1],
671            a[2] * b[0] - a[0] * b[2],
672            a[0] * b[1] - a[1] * b[0],
673        ]
674    }
675
676    macro_rules! test_facet_normals_3d {
677        ($cell:ident) => {
678            paste! {
679
680                #[test]
681                fn [<test_ $cell:lower _facet_normals>]() {
682                    let v = vertices::<f64>(ReferenceCellType::[<$cell>]);
683                    let ns = facet_normals::<f64>(ReferenceCellType::[<$cell>]);
684                    let fs = faces(ReferenceCellType::[<$cell>]);
685
686                    for (f, n) in izip!(fs, ns) {
687                        let tangent0 = (0..3).map(|i| v[f[1]][i] - v[f[0]][i]).collect::<Vec<_>>();
688                        let tangent1 = (0..3).map(|i| v[f[2]][i] - v[f[0]][i]).collect::<Vec<_>>();
689                        for (i, j) in izip!(n, cross(&tangent0, &tangent1)) {
690                            assert_relative_eq!(i, j);
691                        }
692                    }
693                }
694            }
695        };
696    }
697
698    test_facet_normals_2d!(Triangle);
699    test_facet_normals_2d!(Quadrilateral);
700    test_facet_normals_3d!(Tetrahedron);
701    test_facet_normals_3d!(Hexahedron);
702    test_facet_normals_3d!(Prism);
703    test_facet_normals_3d!(Pyramid);
704}
ndelement/reference_cell.rs

ndelement/
reference_cell.rs
