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

ndelement/
reference_cell.rs
