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 local indices of the subentities of the reference cell, with each entity type indexed separately
379pub fn entity_indices(cell: ReferenceCellType) -> Vec<Vec<usize>> {
380    match cell {
381        ReferenceCellType::Point => vec![vec![0], vec![], vec![], vec![]],
382        ReferenceCellType::Interval => vec![vec![0, 1], vec![0], vec![], vec![]],
383        ReferenceCellType::Triangle => vec![vec![0, 1, 2], vec![0, 1, 2], vec![0], vec![]],
384        ReferenceCellType::Quadrilateral => {
385            vec![vec![0, 1, 2, 3], vec![0, 1, 2, 3], vec![0], vec![]]
386        }
387        ReferenceCellType::Tetrahedron => vec![
388            vec![0, 1, 2, 3],
389            vec![0, 1, 2, 3, 4, 5],
390            vec![0, 1, 2, 3],
391            vec![0],
392        ],
393        ReferenceCellType::Hexahedron => vec![
394            vec![0, 1, 2, 3, 4, 5, 6, 7],
395            vec![0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11],
396            vec![0, 1, 2, 3, 4, 5],
397            vec![0],
398        ],
399        ReferenceCellType::Prism => vec![
400            vec![0, 1, 2, 3, 4, 5],
401            vec![0, 1, 2, 3, 4, 5, 6, 7, 8],
402            vec![0, 0, 1, 2, 1],
403            vec![0],
404        ],
405        ReferenceCellType::Pyramid => vec![
406            vec![0, 1, 2, 3, 4],
407            vec![0, 1, 2, 3, 4, 5, 6, 7],
408            vec![0, 0, 1, 2, 3],
409            vec![0],
410        ],
411    }
412}
413
414/// The number of subentities of each dimension
415pub fn entity_counts(cell: ReferenceCellType) -> Vec<usize> {
416    match cell {
417        ReferenceCellType::Point => vec![1, 0, 0, 0],
418        ReferenceCellType::Interval => vec![2, 1, 0, 0],
419        ReferenceCellType::Triangle => vec![3, 3, 1, 0],
420        ReferenceCellType::Quadrilateral => vec![4, 4, 1, 0],
421        ReferenceCellType::Tetrahedron => vec![4, 6, 4, 1],
422        ReferenceCellType::Hexahedron => vec![8, 12, 6, 1],
423        ReferenceCellType::Prism => vec![6, 9, 5, 1],
424        ReferenceCellType::Pyramid => vec![5, 8, 5, 1],
425    }
426}
427
428/// The connectivity of the reference cell
429///
430/// The indices of the result are \[dim0\]\[index0\]\[dim1\]\[\]
431pub fn connectivity(cell: ReferenceCellType) -> Vec<Vec<Vec<Vec<usize>>>> {
432    match cell {
433        ReferenceCellType::Point => vec![vec![vec![vec![0]]]],
434        ReferenceCellType::Interval => vec![
435            vec![vec![vec![0], vec![0]], vec![vec![1], vec![0]]],
436            vec![vec![vec![0, 1], vec![0]]],
437        ],
438        ReferenceCellType::Triangle => vec![
439            vec![
440                vec![vec![0], vec![1, 2], vec![0]],
441                vec![vec![1], vec![0, 2], vec![0]],
442                vec![vec![2], vec![0, 1], vec![0]],
443            ],
444            vec![
445                vec![vec![1, 2], vec![0], vec![0]],
446                vec![vec![0, 2], vec![1], vec![0]],
447                vec![vec![0, 1], vec![2], vec![0]],
448            ],
449            vec![vec![vec![0, 1, 2], vec![0, 1, 2], vec![0]]],
450        ],
451        ReferenceCellType::Quadrilateral => vec![
452            vec![
453                vec![vec![0], vec![0, 1], vec![0]],
454                vec![vec![1], vec![0, 2], vec![0]],
455                vec![vec![2], vec![1, 3], vec![0]],
456                vec![vec![3], vec![2, 3], vec![0]],
457            ],
458            vec![
459                vec![vec![0, 1], vec![0], vec![0]],
460                vec![vec![0, 2], vec![1], vec![0]],
461                vec![vec![1, 3], vec![2], vec![0]],
462                vec![vec![2, 3], vec![3], vec![0]],
463            ],
464            vec![vec![vec![0, 1, 2, 3], vec![0, 1, 2, 3], vec![0]]],
465        ],
466        ReferenceCellType::Tetrahedron => vec![
467            vec![
468                vec![vec![0], vec![3, 4, 5], vec![1, 2, 3], vec![0]],
469                vec![vec![1], vec![1, 2, 5], vec![0, 2, 3], vec![0]],
470                vec![vec![2], vec![0, 2, 4], vec![0, 1, 3], vec![0]],
471                vec![vec![3], vec![0, 1, 3], vec![0, 1, 2], vec![0]],
472            ],
473            vec![
474                vec![vec![2, 3], vec![0], vec![0, 1], vec![0]],
475                vec![vec![1, 3], vec![1], vec![0, 2], vec![0]],
476                vec![vec![1, 2], vec![2], vec![0, 3], vec![0]],
477                vec![vec![0, 3], vec![3], vec![1, 2], vec![0]],
478                vec![vec![0, 2], vec![4], vec![1, 3], vec![0]],
479                vec![vec![0, 1], vec![5], vec![2, 3], vec![0]],
480            ],
481            vec![
482                vec![vec![1, 2, 3], vec![0, 1, 2], vec![0], vec![0]],
483                vec![vec![0, 2, 3], vec![0, 3, 4], vec![1], vec![0]],
484                vec![vec![0, 1, 3], vec![1, 3, 5], vec![2], vec![0]],
485                vec![vec![0, 1, 2], vec![2, 4, 5], vec![3], vec![0]],
486            ],
487            vec![vec![
488                vec![0, 1, 2, 3],
489                vec![0, 1, 2, 3, 4, 5],
490                vec![0, 1, 2, 3],
491                vec![0],
492            ]],
493        ],
494        ReferenceCellType::Hexahedron => vec![
495            vec![
496                vec![vec![0], vec![0, 1, 2], vec![0, 1, 2], vec![0]],
497                vec![vec![1], vec![0, 3, 4], vec![0, 1, 3], vec![0]],
498                vec![vec![2], vec![1, 5, 6], vec![0, 2, 4], vec![0]],
499                vec![vec![3], vec![3, 5, 7], vec![0, 3, 4], vec![0]],
500                vec![vec![4], vec![2, 8, 9], vec![1, 2, 5], vec![0]],
501                vec![vec![5], vec![4, 8, 10], vec![1, 3, 5], vec![0]],
502                vec![vec![6], vec![6, 9, 11], vec![2, 4, 5], vec![0]],
503                vec![vec![7], vec![7, 10, 11], vec![3, 4, 5], vec![0]],
504            ],
505            vec![
506                vec![vec![0, 1], vec![0], vec![0, 1], vec![0]],
507                vec![vec![0, 2], vec![1], vec![0, 2], vec![0]],
508                vec![vec![0, 4], vec![2], vec![1, 2], vec![0]],
509                vec![vec![1, 3], vec![3], vec![0, 3], vec![0]],
510                vec![vec![1, 5], vec![4], vec![1, 3], vec![0]],
511                vec![vec![2, 3], vec![5], vec![0, 4], vec![0]],
512                vec![vec![2, 6], vec![6], vec![2, 4], vec![0]],
513                vec![vec![3, 7], vec![7], vec![3, 4], vec![0]],
514                vec![vec![4, 5], vec![8], vec![1, 5], vec![0]],
515                vec![vec![4, 6], vec![9], vec![2, 5], vec![0]],
516                vec![vec![5, 7], vec![10], vec![3, 5], vec![0]],
517                vec![vec![6, 7], vec![11], vec![4, 5], vec![0]],
518            ],
519            vec![
520                vec![vec![0, 1, 2, 3], vec![0, 1, 3, 5], vec![0], vec![0]],
521                vec![vec![0, 1, 4, 5], vec![0, 2, 4, 8], vec![1], vec![0]],
522                vec![vec![0, 2, 4, 6], vec![1, 2, 6, 9], vec![2], vec![0]],
523                vec![vec![1, 3, 5, 7], vec![3, 4, 7, 10], vec![3], vec![0]],
524                vec![vec![2, 3, 6, 7], vec![5, 6, 7, 11], vec![4], vec![0]],
525                vec![vec![4, 5, 6, 7], vec![8, 9, 10, 11], vec![5], vec![0]],
526            ],
527            vec![vec![
528                vec![0, 1, 2, 3, 4, 5, 6, 7],
529                vec![0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11],
530                vec![0, 1, 2, 3, 4, 5],
531                vec![0],
532            ]],
533        ],
534        ReferenceCellType::Prism => vec![
535            vec![
536                vec![vec![0], vec![0, 1, 2], vec![0, 1, 2], vec![0]],
537                vec![vec![1], vec![0, 3, 4], vec![0, 1, 3], vec![0]],
538                vec![vec![2], vec![1, 3, 5], vec![0, 2, 3], vec![0]],
539                vec![vec![3], vec![2, 6, 7], vec![1, 2, 4], vec![0]],
540                vec![vec![4], vec![4, 6, 8], vec![1, 3, 4], vec![0]],
541                vec![vec![5], vec![5, 7, 8], vec![2, 3, 4], vec![0]],
542            ],
543            vec![
544                vec![vec![0, 1], vec![0], vec![0, 1], vec![0]],
545                vec![vec![0, 2], vec![1], vec![0, 2], vec![0]],
546                vec![vec![0, 3], vec![2], vec![1, 2], vec![0]],
547                vec![vec![1, 2], vec![3], vec![0, 3], vec![0]],
548                vec![vec![1, 4], vec![4], vec![1, 3], vec![0]],
549                vec![vec![2, 5], vec![5], vec![2, 3], vec![0]],
550                vec![vec![3, 4], vec![6], vec![1, 4], vec![0]],
551                vec![vec![3, 5], vec![7], vec![2, 4], vec![0]],
552                vec![vec![4, 5], vec![8], vec![3, 4], vec![0]],
553            ],
554            vec![
555                vec![vec![0, 1, 2], vec![0, 1, 3], vec![0], vec![0]],
556                vec![vec![0, 1, 3, 4], vec![0, 2, 4, 6], vec![1], vec![0]],
557                vec![vec![0, 2, 3, 5], vec![1, 2, 5, 7], vec![2], vec![0]],
558                vec![vec![1, 2, 4, 5], vec![3, 4, 5, 8], vec![3], vec![0]],
559                vec![vec![3, 4, 5], vec![6, 7, 8], vec![4], vec![0]],
560            ],
561            vec![vec![
562                vec![0, 1, 2, 3, 4, 5],
563                vec![0, 1, 2, 3, 4, 5, 6, 7, 8],
564                vec![0, 1, 2, 3, 4],
565                vec![0],
566            ]],
567        ],
568        ReferenceCellType::Pyramid => vec![
569            vec![
570                vec![vec![0], vec![0, 1, 2], vec![0, 1, 2], vec![0]],
571                vec![vec![1], vec![0, 3, 4], vec![0, 1, 3], vec![0]],
572                vec![vec![2], vec![1, 5, 6], vec![0, 2, 4], vec![0]],
573                vec![vec![3], vec![3, 5, 7], vec![0, 3, 4], vec![0]],
574                vec![vec![4], vec![2, 4, 6, 7], vec![1, 2, 3, 4], vec![0]],
575            ],
576            vec![
577                vec![vec![0, 1], vec![0], vec![0, 1], vec![0]],
578                vec![vec![0, 2], vec![1], vec![0, 2], vec![0]],
579                vec![vec![0, 4], vec![2], vec![1, 2], vec![0]],
580                vec![vec![1, 3], vec![3], vec![0, 3], vec![0]],
581                vec![vec![1, 4], vec![4], vec![1, 3], vec![0]],
582                vec![vec![2, 3], vec![5], vec![0, 4], vec![0]],
583                vec![vec![2, 4], vec![6], vec![2, 4], vec![0]],
584                vec![vec![3, 4], vec![7], vec![3, 4], vec![0]],
585            ],
586            vec![
587                vec![vec![0, 1, 2, 3], vec![0, 1, 3, 5], vec![0], vec![0]],
588                vec![vec![0, 1, 4], vec![0, 2, 4], vec![1], vec![0]],
589                vec![vec![0, 2, 4], vec![1, 2, 6], vec![2], vec![0]],
590                vec![vec![1, 3, 4], vec![3, 4, 7], vec![3], vec![0]],
591                vec![vec![2, 3, 4], vec![5, 6, 7], vec![4], vec![0]],
592            ],
593            vec![vec![
594                vec![0, 1, 2, 3, 4],
595                vec![0, 1, 2, 3, 4, 5, 6, 7],
596                vec![0, 1, 2, 3, 4],
597                vec![0],
598            ]],
599        ],
600    }
601}
602
603#[cfg(test)]
604mod test {
605    use super::*;
606    use approx::*;
607    use itertools::izip;
608    use paste::paste;
609
610    macro_rules! test_cell {
611        ($cell:ident) => {
612            paste! {
613
614                #[test]
615                fn [<test_ $cell:lower>]() {
616                    let v = vertices::<f64>(ReferenceCellType::[<$cell>]);
617                    let d = dim(ReferenceCellType::[<$cell>]);
618                    let ec = entity_counts(ReferenceCellType::[<$cell>]);
619                    let et = entity_types(ReferenceCellType::[<$cell>]);
620                    let conn = connectivity(ReferenceCellType::[<$cell>]);
621                    for i in 0..d+1 {
622                        assert_eq!(ec[i], et[i].len());
623                        assert_eq!(ec[i], conn[i].len());
624                    }
625                    assert_eq!(ec[0], v.len());
626                    for i in v {
627                        assert_eq!(i.len(), d);
628                    }
629
630                    for v_n in 0..ec[0] {
631                        let v = &conn[0][v_n][0];
632                        assert_eq!(v, &[v_n]);
633                    }
634                    for e_n in 0..ec[1] {
635                        let vs = &conn[1][e_n][0];
636                        assert_eq!(vs, &edges(ReferenceCellType::[<$cell>])[e_n]);
637                    }
638                    for f_n in 0..ec[2] {
639                        let vs = &conn[2][f_n][0];
640                        assert_eq!(vs, &faces(ReferenceCellType::[<$cell>])[f_n]);
641                    }
642
643                    for e_dim in 0..d {
644                        for e_n in 0..ec[e_dim] {
645                            let e_vertices = &conn[e_dim][e_n][0];
646                            for c_dim in 0..d + 1 {
647                                let connectivity = &conn[e_dim][e_n][c_dim];
648                                if e_dim == c_dim {
649                                    assert_eq!(connectivity.len(), 1);
650                                    assert_eq!(connectivity[0], e_n)
651                                } else {
652                                    for c_n in connectivity {
653                                        let c_vertices = &conn[c_dim][*c_n][0];
654                                        if e_dim < c_dim {
655                                            for i in e_vertices {
656                                                assert!(c_vertices.contains(&i));
657                                            }
658                                        } else {
659                                            for i in c_vertices {
660                                                assert!(e_vertices.contains(&i));
661                                            }
662                                        }
663                                        assert!(connectivity.contains(&c_n));
664                                    }
665                                }
666                            }
667                        }
668                    }
669                }
670            }
671        };
672    }
673
674    test_cell!(Interval);
675    test_cell!(Triangle);
676    test_cell!(Quadrilateral);
677    test_cell!(Tetrahedron);
678    test_cell!(Hexahedron);
679    test_cell!(Prism);
680    test_cell!(Pyramid);
681
682    macro_rules! test_facet_normals_2d {
683        ($cell:ident) => {
684            paste! {
685
686                #[test]
687                fn [<test_ $cell:lower _facet_normals>]() {
688                    let v = vertices::<f64>(ReferenceCellType::[<$cell>]);
689                    let ns = facet_normals::<f64>(ReferenceCellType::[<$cell>]);
690                    let es = edges(ReferenceCellType::[<$cell>]);
691
692                    for (e, n) in izip!(es, ns) {
693                        let tangent = (0..2).map(|i| v[e[1]][i] - v[e[0]][i]).collect::<Vec<_>>();
694                        assert_relative_eq!(n[0], -tangent[1]);
695                        assert_relative_eq!(n[1], tangent[0]);
696                    }
697                }
698            }
699        };
700    }
701
702    fn cross<T: RlstScalar>(a: &[T], b: &[T]) -> Vec<T> {
703        assert_eq!(a.len(), 3);
704        assert_eq!(b.len(), 3);
705        vec![
706            a[1] * b[2] - a[2] * b[1],
707            a[2] * b[0] - a[0] * b[2],
708            a[0] * b[1] - a[1] * b[0],
709        ]
710    }
711
712    macro_rules! test_facet_normals_3d {
713        ($cell:ident) => {
714            paste! {
715
716                #[test]
717                fn [<test_ $cell:lower _facet_normals>]() {
718                    let v = vertices::<f64>(ReferenceCellType::[<$cell>]);
719                    let ns = facet_normals::<f64>(ReferenceCellType::[<$cell>]);
720                    let fs = faces(ReferenceCellType::[<$cell>]);
721
722                    for (f, n) in izip!(fs, ns) {
723                        let tangent0 = (0..3).map(|i| v[f[1]][i] - v[f[0]][i]).collect::<Vec<_>>();
724                        let tangent1 = (0..3).map(|i| v[f[2]][i] - v[f[0]][i]).collect::<Vec<_>>();
725                        for (i, j) in izip!(n, cross(&tangent0, &tangent1)) {
726                            assert_relative_eq!(i, j);
727                        }
728                    }
729                }
730            }
731        };
732    }
733
734    test_facet_normals_2d!(Triangle);
735    test_facet_normals_2d!(Quadrilateral);
736    test_facet_normals_3d!(Tetrahedron);
737    test_facet_normals_3d!(Hexahedron);
738    test_facet_normals_3d!(Prism);
739    test_facet_normals_3d!(Pyramid);
740}
ndelement/reference_cell.rs

ndelement/
reference_cell.rs
