ndgrid/topology/
single_type.rs

1//! Topology for grids where entities of each tdim are the same type
2
3#[cfg(feature = "serde")]
4use crate::traits::ConvertToSerializable;
5use crate::traits::Topology;
6use itertools::izip;
7use ndelement::{orientation::compute_orientation, reference_cell, types::ReferenceCellType};
8use rlst::{DynArray, rlst_dynamic_array};
9use std::collections::HashMap;
10use std::iter::Copied;
11
12/// Topology of a single element grid
13#[derive(Debug)]
14pub struct SingleTypeTopology {
15    dim: usize,
16    pub(crate) ids: Vec<Option<Vec<usize>>>,
17    pub(crate) ids_to_indices: Vec<HashMap<usize, usize>>,
18    entity_types: Vec<ReferenceCellType>,
19    entity_counts: Vec<usize>,
20    pub(crate) downward_connectivity: Vec<Vec<DynArray<usize, 2>>>,
21    pub(crate) upward_connectivity: Vec<Vec<Vec<Vec<usize>>>>,
22    pub(crate) orientation: Vec<Vec<i32>>,
23}
24
25#[cfg(feature = "serde")]
26#[derive(serde::Serialize, Debug, serde::Deserialize)]
27/// Serde serializable topology
28pub struct SerializableTopology {
29    dim: usize,
30    ids: Vec<Option<Vec<usize>>>,
31    ids_to_indices: Vec<HashMap<usize, usize>>,
32    entity_types: Vec<ReferenceCellType>,
33    entity_counts: Vec<usize>,
34    downward_connectivity: Vec<Vec<(Vec<usize>, [usize; 2])>>,
35    upward_connectivity: Vec<Vec<Vec<Vec<usize>>>>,
36    orientation: Vec<Vec<i32>>,
37}
38
39#[cfg(feature = "serde")]
40impl ConvertToSerializable for SingleTypeTopology {
41    type SerializableType = SerializableTopology;
42    fn to_serializable(&self) -> SerializableTopology {
43        SerializableTopology {
44            dim: self.dim,
45            ids: self.ids.clone(),
46            ids_to_indices: self.ids_to_indices.clone(),
47            entity_types: self.entity_types.clone(),
48            entity_counts: self.entity_counts.clone(),
49            downward_connectivity: self
50                .downward_connectivity
51                .iter()
52                .map(|a| {
53                    a.iter()
54                        .map(|b| (b.data().unwrap().to_vec(), b.shape()))
55                        .collect::<Vec<_>>()
56                })
57                .collect::<Vec<_>>(),
58            upward_connectivity: self.upward_connectivity.clone(),
59            orientation: self.orientation.clone(),
60        }
61    }
62    fn from_serializable(s: SerializableTopology) -> Self {
63        Self {
64            dim: s.dim,
65            ids: s.ids,
66            ids_to_indices: s.ids_to_indices,
67            entity_types: s.entity_types,
68            entity_counts: s.entity_counts,
69            downward_connectivity: s
70                .downward_connectivity
71                .iter()
72                .map(|a| {
73                    a.iter()
74                        .map(|(data, shape)| {
75                            let mut c = DynArray::<usize, 2>::from_shape(*shape);
76                            c.data_mut().unwrap().copy_from_slice(data);
77                            c
78                        })
79                        .collect::<Vec<_>>()
80                })
81                .collect::<Vec<_>>(),
82            upward_connectivity: s.upward_connectivity,
83            orientation: s.orientation,
84        }
85    }
86}
87
88impl SingleTypeTopology {
89    /// Create a topology
90    pub fn new(
91        cells: &[usize],
92        cell_type: ReferenceCellType,
93        vertex_ids: Option<Vec<usize>>,
94        cell_ids: Option<Vec<usize>>,
95    ) -> Self {
96        let size = reference_cell::entity_counts(cell_type)[0];
97        let ncells = cells.len() / size;
98        let dim = reference_cell::dim(cell_type);
99
100        // The reference connectivity is a 4-dimensional array.
101        // The first dimension is the topological dimension of the entities for which we want the connectivity.
102        // The second dimension is the index of the entity for which we want to enquire connectivity.
103        // The third dimension is the dimension of the sub-entities for which we want to enquire connectivity.
104        // Hence, for a triangle ref_connectivity[0][1][1] gives us the indices of the two lines that point one
105        // is connected two.
106        let ref_conn = reference_cell::connectivity(cell_type);
107
108        // This gives us all the subentity-types of the cell type. etypes[0] is a list of
109        // all the entity types of dimension 0 (i.e. vertices), etypes[1] is a list of all
110        // the entity types of dimension 1, etc.
111        let etypes = reference_cell::entity_types(cell_type);
112
113        // The following iterates over all the entity types of each dimension and makes
114        // sure that within a dimension all entity types are identical.
115        // Cells where faces are mixture of triangles and quads not supported
116        for t in &etypes {
117            for i in t {
118                if *i != t[0] {
119                    panic!("Unsupported cell type for SingleTypeTopology: {cell_type:?}");
120                }
121            }
122        }
123
124        // This simply collects all the different entity types that are possible within an element.
125        let entity_types = etypes
126            .iter()
127            .filter(|i| !i.is_empty())
128            .map(|i| i[0])
129            .collect::<Vec<_>>();
130
131        // List of entities by dimension
132        // Entity are stored by dimension. We are not storing vertices. Hence, dim - 1.
133        // entities[0] are a HashMap of all zero dimensional entities to their respective entity index.
134        // The key of this HashMap is the ordered set of associated vertex indices for this entity.
135        let mut entities = Vec::<HashMap<Vec<usize>, usize>>::new();
136        if dim > 0 {
137            for _ in 0..dim - 1 {
138                entities.push(HashMap::new());
139            }
140        }
141        // We setup entity counters for each dimension.
142
143        let mut entity_counter = vec![0; entities.len()];
144
145        // Now iterate through all the cells
146        for cell_index in 0..ncells {
147            // Get all the vertices of the cell.
148            let cell = &cells[cell_index * size..(cell_index + 1) * size];
149            // Iterate over topological dimension.
150            // The variable e_i is a mutable reference to the list of entities of dimension i.
151            // The variable rc_i is a reference to the connectivity of the cell type for entities of dimension i.
152            // Hence, rc_i[0] is the connectivity information about the entity with reference index 0 and
153            // dimension i.
154            for (entity_dim, (e_i, rc_i)) in izip!(
155                entities.iter_mut(),
156                // Skip the vertices in the connectivity
157                ref_conn.iter().take(dim).skip(1),
158            )
159            .enumerate()
160            {
161                // We iterate through the concrete entities of dimension i and the corresponding
162                // entity types. c_ij is a reference to the connectivity information of the jth entity
163                // with dimension i.
164                for c_ij in rc_i {
165                    // c_ij[0] is the list of reference cell vertex indices that are associated with the jth entity of dimension i.
166                    // cell[*i] below maps the local reference cell vertex index to the actual id of the vertex.
167                    // Hence, the following command gives us all the vertex indicies of the entity.
168                    let mut entity = c_ij[0].iter().map(|i| cell[*i]).collect::<Vec<_>>();
169                    // We sort entities so that entities with same vertices but different order of vertices
170                    // are treated as the same entity.
171                    entity.sort();
172                    // An entity only gets added if it has not been added before.
173                    e_i.entry(entity).or_insert_with(|| {
174                        let old = entity_counter[entity_dim];
175                        entity_counter[entity_dim] += 1;
176                        old
177                    });
178                }
179            }
180        }
181        // Number of entities by dimension
182        let mut entity_counts = vec![0; dim + 1];
183        // Vertices are indexed consecutively starting from 0.
184        // So to get the number of vertices take the highest index
185        // and add 1.
186        entity_counts[0] = cells.iter().max().unwrap() + 1;
187        entity_counts[dim] = ncells;
188        for d in 1..dim {
189            // In entities the dimension 0 is skipped.
190            // So entity_counts[d] is the number of entities in
191            // entities[d - 1].
192            entity_counts[d] = entities[d - 1].len();
193        }
194
195        // Downward connectivity: The entities of dimension dim1 that are subentities of
196        // entities  of dimension dim0 (with dim0>=dim1) (eg edges of a triangle, vertices
197        // of a tetrahedron, etc)
198        // downward_connectivity[dim0][dim1][[.., dim0_entity_index]] = [dim1_entity_index]
199        // Hence, downward_connectivity[2][1][0, 5] for a triangle grid gives the index of the
200        // first edge of the 6th triangle.
201        // The following loop only creates space for the downward connectivity. It does not fill
202        // the actual entries.
203        let mut downward_connectivity = entity_counts
204            .iter()
205            .enumerate()
206            // i is the topological dimension of the entities for which we want the connectivity.
207            // j is the reference cell index of the entity of dimension j.
208            .map(|(i, j)| {
209                // ref_conn[i][0] gives us the first entity of dimension i. We don't need to look
210                // at other entities of dimension j since the information about number of
211                // subentities is the same.
212                ref_conn[i][0]
213                    // We iterate through all the dimensions it can be connected to.
214                    .iter()
215                    // We take i + 1 dimensions of the connection.
216                    // If i = 2 (triangle) then we take vertices, edges, and triangles.
217                    // That's why take(i + 1).
218                    .take(i + 1)
219                    // We now iterate through the connected dimensions, e.g. in the case
220                    // of i = 2, the dimensions 0, 1, 2. These are in r.
221                    // r.len() is the number of these entities connected to our i-dim entity.
222                    // and j is the actual number of i-dimensional entities.
223                    .map(|r| rlst_dynamic_array!(usize, [r.len(), *j]))
224                    .collect::<Vec<_>>()
225            })
226            .collect::<Vec<_>>();
227
228        // Upward connectivity: The entities of dimension dim1 that are superentities of
229        // entities of dimension dim0 (with dim0<dim1) (eg triangles connected to an edge,
230        // tetrahedra connected to a vertex, etc)
231        // upward_connectivity[dim0][dim1 - dim0 - 1][dim0_entity_index][..] = [dim1_entity_index]
232        let mut upward_connectivity = entity_counts
233            .iter()
234            .take(dim)
235            .enumerate()
236            .map(|(i, j)| vec![vec![vec![]; *j]; dim - i])
237            .collect::<Vec<Vec<Vec<Vec<usize>>>>>();
238
239        // downward_connectivity[d][d][i] = [i] (ie each entity is a sub-entity of itself)
240        for (d, dc) in downward_connectivity.iter_mut().enumerate() {
241            for (i, mut j) in dc[d].col_iter_mut().enumerate() {
242                // Each entity has exactly one sub-entity of the same dimension, namely itself.
243                // Hence, the matrix is a 1x1 matrix that refers back to itself.
244                j[[0]] = i;
245            }
246        }
247
248        // downward_connectivity[dim][0] = vertices of each cell
249        // This is just the vertices that each cell is made of.
250        for (i, mut dc_d0i) in downward_connectivity[dim][0].col_iter_mut().enumerate() {
251            for (dc_d0ij, c_j) in izip!(dc_d0i.iter_mut(), &cells[i * size..(i + 1) * size]) {
252                *dc_d0ij = *c_j;
253                // We can also fill out the upward connectivity if dim > 0.
254                if dim > 0 {
255                    upward_connectivity[0][dim - 1][*c_j].push(i);
256                }
257            }
258        }
259
260        // downward_connectivity[i][0] = vertices of entity
261        for (i, (es_i, dc_i)) in izip!(
262            entities.iter(),
263            downward_connectivity.iter_mut().take(dim).skip(1)
264        )
265        .enumerate()
266        {
267            for es_ij in es_i.iter() {
268                let entity_vertices = es_ij.0;
269                let entity_index = es_ij.1;
270                let mut dc_i0j = dc_i[0].r_mut().slice::<1>(1, *entity_index);
271                for (vertex, dc_i0jk) in izip!(entity_vertices.iter(), dc_i0j.iter_mut()) {
272                    *dc_i0jk = *vertex;
273                    if !upward_connectivity[0][i][*vertex].contains(entity_index) {
274                        upward_connectivity[0][i][*vertex].push(*entity_index);
275                    }
276                }
277            }
278        }
279
280        // We now have to fill the upward and downward connectivity for the other cases. This can only happen if dim > 0.
281        if dim > 0 {
282            // We collect the cell entities in an array. We do not need the connectivity of the vertices as that
283            // has already been dealt with. Hence, the skip(1). The index i specifies the subentity of the cell.
284            // and i.len() is the number of subentities of the given dimension.
285            // Hence, cell_entities[0] is a vector of length of all subentities of dimension 1 in the reference cell.
286            let mut cell_entities = ref_conn
287                .iter()
288                .skip(1)
289                .map(|i| vec![0; i.len()])
290                .collect::<Vec<_>>();
291            // Now we iterate through each individual cell.
292            for cell_index in 0..ncells {
293                // Collect indices of each subentity of the cell.
294                // The subentities of dimension the cell itself are just the cell indicies. Remember that a dimension index of
295                // dim - 1 denotes actually the dim dimensional subentity as we skipped vertices.
296                cell_entities[dim - 1][0] = cell_index;
297                // We now get all the vertices of the current cell.
298                let cell = &cells[cell_index * size..(cell_index + 1) * size];
299                // We now iterate through the connectivity dimensions. e_i is all entities of dimension i + 1.
300                // ce_i is the list of all cell entities of dimension i + 1. rc_i is the connectivity information
301                // of dimension i + 1 (because we do skip(1))
302                for (e_i, ce_i, rc_i) in izip!(
303                    entities.iter(),
304                    cell_entities.iter_mut(),
305                    ref_conn.iter().skip(1),
306                ) {
307                    // This iterates over the actual sub entities.
308                    // - ce_ij is the index of the jth subentity of dimension i + 1 (that's what we want to get in this loop).
309                    // - rc_ij is the reference connectivity information of the ith subentity of dimension i + 1.
310                    for (ce_ij, rc_ij) in izip!(ce_i.iter_mut(), rc_i) {
311                        // We get the actual entity by mapping reference cell vertices to actual vertices
312                        let mut entity = rc_ij[0].iter().map(|i| cell[*i]).collect::<Vec<_>>();
313                        entity.sort();
314                        *ce_ij = *e_i.get(&entity).unwrap();
315                    }
316                }
317                // The following is a double loop that enters all k-dimensional entities connected to all
318                // i-dimensional entities into the downward connectivities.
319                // Copy these indices into connectivity for each dim
320                // Loop over dim0
321                // We now fill out the downward connectivity.
322                // dc_i is downward connectivity for dim i.
323                // rc_i is the reference connectivity for dim i.
324                // ce_i are the cell entitiies for dim i.
325                // The loop ignores dim 0 and dim 1 entities as everything for dim 1 entities
326                // has been processed at this point (entities with themselves have been done, and entities with vertices).
327                // The index i can only be 0 or one. We have 4 dimensions in total (0, 1, 2, 3). If we skip 0 and 1 then
328                // the possible dimensions are 2 and 3. Since i is a counter i=0 is associated with dimension 2 and i=1 is
329                // associated with dimension 3.
330                for (i, (dc_i, rc_i, ce_i)) in izip!(
331                    downward_connectivity.iter_mut().skip(2),
332                    ref_conn.iter().skip(2),
333                    cell_entities.iter().skip(1)
334                )
335                .enumerate()
336                {
337                    // ce_ij is the jth subentity of dimension i.
338                    // rc_ij is the corresponding reference connectivity of the jth subentity of dimension i.
339                    for (ce_ij, rc_ij) in izip!(ce_i, rc_i) {
340                        // k loops over the dimensions of the cell entities.
341                        // rc_ijk is the kth dimensional reference connectivity of the jth subentity of dimension i.
342                        // dc_ik is the kth dimensional downward connectivity of the subentities of dimension i.
343                        // If i = 0 we consider 2-dimensional entities (since we have skipped vertices and edges). Hence,
344                        // dc_i.iter_mut().take(i+2) takes the first two dimensions of the downward connectivity,
345                        // that is vertices and edges. It then skips over the vertices. So k = 0 is only option.
346                        // If i = 1 we consider 3-dimensional entities. dc_i.iter_mut().take(i+2) takes the first 3
347                        // dimensions, that is vertices, edges, faces. It then skips over the vertices. So k can be 0 and 1.
348                        for (k, (dc_ik, rc_ijk, ce_k)) in izip!(
349                            dc_i.iter_mut().take(i + 2).skip(1),
350                            rc_ij.iter().take(i + 2).skip(1),
351                            &cell_entities
352                        )
353                        .enumerate()
354                        {
355                            // l is the concrete index of the kth dimensional connectivity of the jth subentity of dimension i.
356                            // We copy this into the right position of the downward connectivity.
357                            // Remember that ce_ij is the jth subentity of dimension i.
358                            for (l, rc_ijkl) in rc_ijk.iter().enumerate() {
359                                dc_ik[[l, *ce_ij]] = ce_k[*rc_ijkl];
360                                // Now will in reverse the corresponding upward connectivity.
361                                // If i = 0 then we consider 2-dimensional entities (faces). k=0 corresponds to edges.
362                                // The upward connectivity between the two is upward_connectivity[1][0].
363                                // The general setup is upward_connectivity[dim0][dim1 - dim0 - 1][dim0_entity_index][..] = [dim1_entity_index].
364                                // We hae dim0 = k + 1 and dim1 = i + 2. Hence, dim1 - dim0 - 1 = i - k.
365                                if !upward_connectivity[k + 1][i - k][ce_k[*rc_ijkl]]
366                                    .contains(ce_ij)
367                                {
368                                    upward_connectivity[k + 1][i - k][ce_k[*rc_ijkl]].push(*ce_ij);
369                                }
370                            }
371                        }
372                    }
373                }
374            }
375        }
376
377        let mut ids = vec![vertex_ids];
378        for _ in 1..dim {
379            ids.push(None);
380        }
381        ids.push(cell_ids);
382
383        let ids_to_indices = ids
384            .iter()
385            .map(|ids_option| {
386                let mut ie = HashMap::new();
387                if let Some(ids) = ids_option {
388                    for (i, j) in ids.iter().enumerate() {
389                        ie.insert(*j, i);
390                    }
391                }
392                ie
393            })
394            .collect::<Vec<_>>();
395
396        let mut orientation = vec![];
397
398        for d in 1..dim + 1 {
399            let dc = &downward_connectivity[d][0];
400            orientation.push(
401                (0..dc.shape()[1])
402                    .map(|i| {
403                        compute_orientation(
404                            entity_types[d],
405                            &dc.data().unwrap()[i * dc.shape()[0]..(i + 1) * dc.shape()[0]],
406                        )
407                    })
408                    .collect::<Vec<_>>(),
409            );
410        }
411
412        Self {
413            dim,
414            ids,
415            ids_to_indices,
416            entity_types,
417            entity_counts,
418            downward_connectivity,
419            upward_connectivity,
420            orientation,
421        }
422    }
423    /// Topological dimension
424    pub fn dim(&self) -> usize {
425        self.dim
426    }
427    /// Entity types
428    pub fn entity_types(&self) -> &[ReferenceCellType] {
429        &self.entity_types
430    }
431    /// Entity count
432    pub fn entity_counts(&self) -> &[usize] {
433        &self.entity_counts
434    }
435    /// Entity count
436    pub fn entity_count(&self, entity_type: ReferenceCellType) -> usize {
437        if !self.entity_types.contains(&entity_type) {
438            0
439        } else {
440            self.entity_counts[reference_cell::dim(entity_type)]
441        }
442    }
443    /// Downward connectivity
444    pub fn downward_connectivity(&self) -> &[Vec<DynArray<usize, 2>>] {
445        &self.downward_connectivity
446    }
447    /// Upward connectivity
448    pub fn upward_connectivity(&self) -> &[Vec<Vec<Vec<usize>>>] {
449        &self.upward_connectivity
450    }
451    /// Cell sub-entity index
452    pub fn cell_entity_index(
453        &self,
454        cell_index: usize,
455        entity_dim: usize,
456        entity_index: usize,
457    ) -> usize {
458        self.downward_connectivity[self.dim][entity_dim][[entity_index, cell_index]]
459    }
460    /// Entity id
461    pub fn entity_id(&self, entity_dim: usize, entity_index: usize) -> Option<usize> {
462        self.ids[entity_dim].as_ref().map(|a| a[entity_index])
463    }
464}
465
466/// Topology of a cell
467#[derive(Debug)]
468pub struct SingleTypeEntityTopology<'a> {
469    topology: &'a SingleTypeTopology,
470    entity_index: usize,
471    dim: usize,
472}
473
474impl<'t> SingleTypeEntityTopology<'t> {
475    /// Create new
476    pub fn new(
477        topology: &'t SingleTypeTopology,
478        entity_type: ReferenceCellType,
479        entity_index: usize,
480    ) -> Self {
481        Self {
482            topology,
483            entity_index,
484            dim: reference_cell::dim(entity_type),
485        }
486    }
487}
488impl Topology for SingleTypeEntityTopology<'_> {
489    type EntityDescriptor = ReferenceCellType;
490    type EntityIndexIter<'a>
491        = Copied<std::slice::Iter<'a, usize>>
492    where
493        Self: 'a;
494
495    type ConnectedEntityIndexIter<'a>
496        = Copied<std::slice::Iter<'a, usize>>
497    where
498        Self: 'a;
499
500    fn connected_entity_iter(
501        &self,
502        entity_type: ReferenceCellType,
503    ) -> Copied<std::slice::Iter<'_, usize>> {
504        let dim = reference_cell::dim(entity_type);
505        if entity_type == self.topology.entity_types()[dim] {
506            self.topology.upward_connectivity[self.dim][dim - self.dim - 1][self.entity_index]
507                .iter()
508                .copied()
509        } else {
510            [].iter().copied()
511        }
512    }
513
514    fn sub_entity_iter(
515        &self,
516        entity_type: ReferenceCellType,
517    ) -> Copied<std::slice::Iter<'_, usize>> {
518        let dim = reference_cell::dim(entity_type);
519        if entity_type == self.topology.entity_types()[dim] {
520            let rows = self.topology.downward_connectivity[self.dim][dim].shape()[0];
521            self.topology.downward_connectivity[self.dim][dim]
522                .data()
523                .unwrap()[rows * self.entity_index..rows * (self.entity_index + 1)]
524                .iter()
525                .copied()
526        } else {
527            [].iter().copied()
528        }
529    }
530
531    fn sub_entity(&self, entity_type: ReferenceCellType, index: usize) -> usize {
532        let dim = reference_cell::dim(entity_type);
533        if entity_type != self.topology.entity_types()[dim] {
534            panic!("Invalid entity type");
535        }
536        self.topology.downward_connectivity[self.dim][dim][[index, self.entity_index]]
537    }
538
539    fn orientation(&self) -> i32 {
540        if self.dim == 0 {
541            0
542        } else {
543            self.topology.orientation[self.dim - 1][self.entity_index]
544        }
545    }
546}
547
548#[cfg(test)]
549mod test {
550    use super::*;
551
552    fn example_topology_point() -> SingleTypeTopology {
553        //! An example topology
554        SingleTypeTopology::new(&[0, 1], ReferenceCellType::Point, None, None)
555    }
556
557    fn example_topology_interval() -> SingleTypeTopology {
558        //! An example topology
559        SingleTypeTopology::new(&[0, 1, 1, 2], ReferenceCellType::Interval, None, None)
560    }
561
562    fn example_topology_triangle() -> SingleTypeTopology {
563        //! An example topology
564        SingleTypeTopology::new(&[0, 1, 2, 2, 1, 3], ReferenceCellType::Triangle, None, None)
565    }
566
567    fn example_topology_tetrahedron() -> SingleTypeTopology {
568        //! An example topology
569        SingleTypeTopology::new(
570            &[0, 1, 2, 3, 4, 0, 2, 3],
571            ReferenceCellType::Tetrahedron,
572            None,
573            None,
574        )
575    }
576
577    #[test]
578    #[should_panic]
579    fn test_prism() {
580        let _ = SingleTypeTopology::new(&[0, 1, 2, 3, 4, 5], ReferenceCellType::Prism, None, None);
581    }
582    #[test]
583    #[should_panic]
584    fn test_pyramid() {
585        let _ = SingleTypeTopology::new(&[0, 1, 2, 3, 4], ReferenceCellType::Pyramid, None, None);
586    }
587
588    #[test]
589    fn test_sub_entities_point() {
590        //! Test sub-entities of a point
591        let t = example_topology_point();
592        let cell0 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Point, 0);
593        assert_eq!(cell0.sub_entity(ReferenceCellType::Point, 0), 0);
594        let cell1 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Point, 1);
595        assert_eq!(cell1.sub_entity(ReferenceCellType::Point, 0), 1);
596    }
597
598    #[test]
599    fn test_sub_entities_interval() {
600        //! Test sub-entities of an interval
601        let t = example_topology_interval();
602        let cell0 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Interval, 0);
603        assert_eq!(cell0.sub_entity(ReferenceCellType::Point, 0), 0);
604        assert_eq!(cell0.sub_entity(ReferenceCellType::Point, 1), 1);
605        assert_eq!(cell0.sub_entity(ReferenceCellType::Interval, 0), 0);
606        let cell1 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Interval, 1);
607        assert_eq!(cell1.sub_entity(ReferenceCellType::Point, 0), 1);
608        assert_eq!(cell1.sub_entity(ReferenceCellType::Point, 1), 2);
609        assert_eq!(cell1.sub_entity(ReferenceCellType::Interval, 0), 1);
610
611        let vertex0 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Point, 0);
612        assert_eq!(vertex0.sub_entity(ReferenceCellType::Point, 0), 0);
613        let vertex1 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Point, 1);
614        assert_eq!(vertex1.sub_entity(ReferenceCellType::Point, 0), 1);
615        let vertex2 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Point, 2);
616        assert_eq!(vertex2.sub_entity(ReferenceCellType::Point, 0), 2);
617    }
618
619    #[test]
620    fn test_sub_entities_triangle() {
621        //! Test sub-entities of a triangle
622        let t = example_topology_triangle();
623        let cell0 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Triangle, 0);
624        assert_eq!(cell0.sub_entity(ReferenceCellType::Point, 0), 0);
625        assert_eq!(cell0.sub_entity(ReferenceCellType::Point, 1), 1);
626        assert_eq!(cell0.sub_entity(ReferenceCellType::Point, 2), 2);
627        assert_eq!(cell0.sub_entity(ReferenceCellType::Interval, 0), 0);
628        assert_eq!(cell0.sub_entity(ReferenceCellType::Interval, 1), 1);
629        assert_eq!(cell0.sub_entity(ReferenceCellType::Interval, 2), 2);
630        assert_eq!(cell0.sub_entity(ReferenceCellType::Triangle, 0), 0);
631        let cell1 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Triangle, 1);
632        assert_eq!(cell1.sub_entity(ReferenceCellType::Point, 0), 2);
633        assert_eq!(cell1.sub_entity(ReferenceCellType::Point, 1), 1);
634        assert_eq!(cell1.sub_entity(ReferenceCellType::Point, 2), 3);
635        assert_eq!(cell1.sub_entity(ReferenceCellType::Interval, 0), 3);
636        assert_eq!(cell1.sub_entity(ReferenceCellType::Interval, 1), 4);
637        assert_eq!(cell1.sub_entity(ReferenceCellType::Interval, 2), 0);
638        assert_eq!(cell1.sub_entity(ReferenceCellType::Triangle, 0), 1);
639
640        let edge0 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Interval, 0);
641        assert_eq!(edge0.sub_entity(ReferenceCellType::Point, 0), 1);
642        assert_eq!(edge0.sub_entity(ReferenceCellType::Point, 1), 2);
643        assert_eq!(edge0.sub_entity(ReferenceCellType::Interval, 0), 0);
644        let edge1 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Interval, 1);
645        assert_eq!(edge1.sub_entity(ReferenceCellType::Point, 0), 0);
646        assert_eq!(edge1.sub_entity(ReferenceCellType::Point, 1), 2);
647        assert_eq!(edge1.sub_entity(ReferenceCellType::Interval, 0), 1);
648        let edge2 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Interval, 2);
649        assert_eq!(edge2.sub_entity(ReferenceCellType::Point, 0), 0);
650        assert_eq!(edge2.sub_entity(ReferenceCellType::Point, 1), 1);
651        assert_eq!(edge2.sub_entity(ReferenceCellType::Interval, 0), 2);
652        let edge3 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Interval, 3);
653        assert_eq!(edge3.sub_entity(ReferenceCellType::Point, 0), 1);
654        assert_eq!(edge3.sub_entity(ReferenceCellType::Point, 1), 3);
655        assert_eq!(edge3.sub_entity(ReferenceCellType::Interval, 0), 3);
656        let edge4 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Interval, 4);
657        assert_eq!(edge4.sub_entity(ReferenceCellType::Point, 0), 2);
658        assert_eq!(edge4.sub_entity(ReferenceCellType::Point, 1), 3);
659        assert_eq!(edge4.sub_entity(ReferenceCellType::Interval, 0), 4);
660
661        let vertex0 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Point, 0);
662        assert_eq!(vertex0.sub_entity(ReferenceCellType::Point, 0), 0);
663        let vertex1 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Point, 1);
664        assert_eq!(vertex1.sub_entity(ReferenceCellType::Point, 0), 1);
665        let vertex2 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Point, 2);
666        assert_eq!(vertex2.sub_entity(ReferenceCellType::Point, 0), 2);
667        let vertex3 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Point, 3);
668        assert_eq!(vertex3.sub_entity(ReferenceCellType::Point, 0), 3);
669    }
670
671    #[test]
672    fn test_sub_entities_tetrahedron() {
673        //! Test sub-entities of a tetrahedron
674        let t = example_topology_tetrahedron();
675        let cell0 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Tetrahedron, 0);
676        assert_eq!(cell0.sub_entity(ReferenceCellType::Point, 0), 0);
677        assert_eq!(cell0.sub_entity(ReferenceCellType::Point, 1), 1);
678        assert_eq!(cell0.sub_entity(ReferenceCellType::Point, 2), 2);
679        assert_eq!(cell0.sub_entity(ReferenceCellType::Point, 3), 3);
680        assert_eq!(cell0.sub_entity(ReferenceCellType::Interval, 0), 0);
681        assert_eq!(cell0.sub_entity(ReferenceCellType::Interval, 1), 1);
682        assert_eq!(cell0.sub_entity(ReferenceCellType::Interval, 2), 2);
683        assert_eq!(cell0.sub_entity(ReferenceCellType::Interval, 3), 3);
684        assert_eq!(cell0.sub_entity(ReferenceCellType::Interval, 4), 4);
685        assert_eq!(cell0.sub_entity(ReferenceCellType::Interval, 5), 5);
686        assert_eq!(cell0.sub_entity(ReferenceCellType::Triangle, 0), 0);
687        assert_eq!(cell0.sub_entity(ReferenceCellType::Triangle, 1), 1);
688        assert_eq!(cell0.sub_entity(ReferenceCellType::Triangle, 2), 2);
689        assert_eq!(cell0.sub_entity(ReferenceCellType::Triangle, 3), 3);
690        assert_eq!(cell0.sub_entity(ReferenceCellType::Tetrahedron, 0), 0);
691        let cell1 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Tetrahedron, 1);
692        assert_eq!(cell1.sub_entity(ReferenceCellType::Point, 0), 4);
693        assert_eq!(cell1.sub_entity(ReferenceCellType::Point, 1), 0);
694        assert_eq!(cell1.sub_entity(ReferenceCellType::Point, 2), 2);
695        assert_eq!(cell1.sub_entity(ReferenceCellType::Point, 3), 3);
696        assert_eq!(cell1.sub_entity(ReferenceCellType::Interval, 0), 0);
697        assert_eq!(cell1.sub_entity(ReferenceCellType::Interval, 1), 3);
698        assert_eq!(cell1.sub_entity(ReferenceCellType::Interval, 2), 4);
699        assert_eq!(cell1.sub_entity(ReferenceCellType::Interval, 3), 6);
700        assert_eq!(cell1.sub_entity(ReferenceCellType::Interval, 4), 7);
701        assert_eq!(cell1.sub_entity(ReferenceCellType::Interval, 5), 8);
702        assert_eq!(cell1.sub_entity(ReferenceCellType::Triangle, 0), 1);
703        assert_eq!(cell1.sub_entity(ReferenceCellType::Triangle, 1), 4);
704        assert_eq!(cell1.sub_entity(ReferenceCellType::Triangle, 2), 5);
705        assert_eq!(cell1.sub_entity(ReferenceCellType::Triangle, 3), 6);
706        assert_eq!(cell1.sub_entity(ReferenceCellType::Tetrahedron, 0), 1);
707
708        let face0 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Triangle, 0);
709        assert_eq!(face0.sub_entity(ReferenceCellType::Point, 0), 1);
710        assert_eq!(face0.sub_entity(ReferenceCellType::Point, 1), 2);
711        assert_eq!(face0.sub_entity(ReferenceCellType::Point, 2), 3);
712        assert_eq!(face0.sub_entity(ReferenceCellType::Interval, 0), 0);
713        assert_eq!(face0.sub_entity(ReferenceCellType::Interval, 1), 1);
714        assert_eq!(face0.sub_entity(ReferenceCellType::Interval, 2), 2);
715        assert_eq!(face0.sub_entity(ReferenceCellType::Triangle, 0), 0);
716        let face1 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Triangle, 1);
717        assert_eq!(face1.sub_entity(ReferenceCellType::Point, 0), 0);
718        assert_eq!(face1.sub_entity(ReferenceCellType::Point, 1), 2);
719        assert_eq!(face1.sub_entity(ReferenceCellType::Point, 2), 3);
720        assert_eq!(face1.sub_entity(ReferenceCellType::Interval, 0), 0);
721        assert_eq!(face1.sub_entity(ReferenceCellType::Interval, 1), 3);
722        assert_eq!(face1.sub_entity(ReferenceCellType::Interval, 2), 4);
723        assert_eq!(face1.sub_entity(ReferenceCellType::Triangle, 0), 1);
724        let face2 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Triangle, 2);
725        assert_eq!(face2.sub_entity(ReferenceCellType::Point, 0), 0);
726        assert_eq!(face2.sub_entity(ReferenceCellType::Point, 1), 1);
727        assert_eq!(face2.sub_entity(ReferenceCellType::Point, 2), 3);
728        assert_eq!(face2.sub_entity(ReferenceCellType::Interval, 0), 1);
729        assert_eq!(face2.sub_entity(ReferenceCellType::Interval, 1), 3);
730        assert_eq!(face2.sub_entity(ReferenceCellType::Interval, 2), 5);
731        assert_eq!(face2.sub_entity(ReferenceCellType::Triangle, 0), 2);
732        let face3 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Triangle, 3);
733        assert_eq!(face3.sub_entity(ReferenceCellType::Point, 0), 0);
734        assert_eq!(face3.sub_entity(ReferenceCellType::Point, 1), 1);
735        assert_eq!(face3.sub_entity(ReferenceCellType::Point, 2), 2);
736        assert_eq!(face3.sub_entity(ReferenceCellType::Interval, 0), 2);
737        assert_eq!(face3.sub_entity(ReferenceCellType::Interval, 1), 4);
738        assert_eq!(face3.sub_entity(ReferenceCellType::Interval, 2), 5);
739        assert_eq!(face3.sub_entity(ReferenceCellType::Triangle, 0), 3);
740        let face4 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Triangle, 4);
741        assert_eq!(face4.sub_entity(ReferenceCellType::Point, 0), 2);
742        assert_eq!(face4.sub_entity(ReferenceCellType::Point, 1), 3);
743        assert_eq!(face4.sub_entity(ReferenceCellType::Point, 2), 4);
744        assert_eq!(face4.sub_entity(ReferenceCellType::Interval, 0), 0);
745        assert_eq!(face4.sub_entity(ReferenceCellType::Interval, 1), 6);
746        assert_eq!(face4.sub_entity(ReferenceCellType::Interval, 2), 7);
747        assert_eq!(face4.sub_entity(ReferenceCellType::Triangle, 0), 4);
748        let face5 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Triangle, 5);
749        assert_eq!(face5.sub_entity(ReferenceCellType::Point, 0), 0);
750        assert_eq!(face5.sub_entity(ReferenceCellType::Point, 1), 3);
751        assert_eq!(face5.sub_entity(ReferenceCellType::Point, 2), 4);
752        assert_eq!(face5.sub_entity(ReferenceCellType::Interval, 0), 3);
753        assert_eq!(face5.sub_entity(ReferenceCellType::Interval, 1), 6);
754        assert_eq!(face5.sub_entity(ReferenceCellType::Interval, 2), 8);
755        assert_eq!(face5.sub_entity(ReferenceCellType::Triangle, 0), 5);
756        let face6 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Triangle, 6);
757        assert_eq!(face6.sub_entity(ReferenceCellType::Point, 0), 0);
758        assert_eq!(face6.sub_entity(ReferenceCellType::Point, 1), 2);
759        assert_eq!(face6.sub_entity(ReferenceCellType::Point, 2), 4);
760        assert_eq!(face6.sub_entity(ReferenceCellType::Interval, 0), 4);
761        assert_eq!(face6.sub_entity(ReferenceCellType::Interval, 1), 7);
762        assert_eq!(face6.sub_entity(ReferenceCellType::Interval, 2), 8);
763        assert_eq!(face6.sub_entity(ReferenceCellType::Triangle, 0), 6);
764
765        let edge0 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Interval, 0);
766        assert_eq!(edge0.sub_entity(ReferenceCellType::Point, 0), 2);
767        assert_eq!(edge0.sub_entity(ReferenceCellType::Point, 1), 3);
768        assert_eq!(edge0.sub_entity(ReferenceCellType::Interval, 0), 0);
769        let edge1 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Interval, 1);
770        assert_eq!(edge1.sub_entity(ReferenceCellType::Point, 0), 1);
771        assert_eq!(edge1.sub_entity(ReferenceCellType::Point, 1), 3);
772        assert_eq!(edge1.sub_entity(ReferenceCellType::Interval, 0), 1);
773        let edge2 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Interval, 2);
774        assert_eq!(edge2.sub_entity(ReferenceCellType::Point, 0), 1);
775        assert_eq!(edge2.sub_entity(ReferenceCellType::Point, 1), 2);
776        assert_eq!(edge2.sub_entity(ReferenceCellType::Interval, 0), 2);
777        let edge3 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Interval, 3);
778        assert_eq!(edge3.sub_entity(ReferenceCellType::Point, 0), 0);
779        assert_eq!(edge3.sub_entity(ReferenceCellType::Point, 1), 3);
780        assert_eq!(edge3.sub_entity(ReferenceCellType::Interval, 0), 3);
781        let edge4 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Interval, 4);
782        assert_eq!(edge4.sub_entity(ReferenceCellType::Point, 0), 0);
783        assert_eq!(edge4.sub_entity(ReferenceCellType::Point, 1), 2);
784        assert_eq!(edge4.sub_entity(ReferenceCellType::Interval, 0), 4);
785        let edge5 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Interval, 5);
786        assert_eq!(edge5.sub_entity(ReferenceCellType::Point, 0), 0);
787        assert_eq!(edge5.sub_entity(ReferenceCellType::Point, 1), 1);
788        assert_eq!(edge5.sub_entity(ReferenceCellType::Interval, 0), 5);
789        let edge6 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Interval, 6);
790        assert_eq!(edge6.sub_entity(ReferenceCellType::Point, 0), 3);
791        assert_eq!(edge6.sub_entity(ReferenceCellType::Point, 1), 4);
792        assert_eq!(edge6.sub_entity(ReferenceCellType::Interval, 0), 6);
793        let edge7 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Interval, 7);
794        assert_eq!(edge7.sub_entity(ReferenceCellType::Point, 0), 2);
795        assert_eq!(edge7.sub_entity(ReferenceCellType::Point, 1), 4);
796        assert_eq!(edge7.sub_entity(ReferenceCellType::Interval, 0), 7);
797        let edge8 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Interval, 8);
798        assert_eq!(edge8.sub_entity(ReferenceCellType::Point, 0), 0);
799        assert_eq!(edge8.sub_entity(ReferenceCellType::Point, 1), 4);
800        assert_eq!(edge8.sub_entity(ReferenceCellType::Interval, 0), 8);
801
802        let vertex0 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Point, 0);
803        assert_eq!(vertex0.sub_entity(ReferenceCellType::Point, 0), 0);
804        let vertex1 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Point, 1);
805        assert_eq!(vertex1.sub_entity(ReferenceCellType::Point, 0), 1);
806        let vertex2 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Point, 2);
807        assert_eq!(vertex2.sub_entity(ReferenceCellType::Point, 0), 2);
808        let vertex3 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Point, 3);
809        assert_eq!(vertex3.sub_entity(ReferenceCellType::Point, 0), 3);
810        let vertex4 = SingleTypeEntityTopology::new(&t, ReferenceCellType::Point, 4);
811        assert_eq!(vertex4.sub_entity(ReferenceCellType::Point, 0), 4);
812    }
813
814    macro_rules! make_tests {
815        ($cellname:ident) => {
816            paste::item! {
817                #[test]
818                fn [< test_up_and_down_ $cellname >]() {
819                    //! Test that upward and downward connectivities agree
820                    let t = [< example_topology_ $cellname >]();
821
822                    for (i, dc_i) in t.downward_connectivity.iter().enumerate() {
823                        for (j, dc_ij) in dc_i.iter().enumerate() {
824                            if i != j {
825                                let uc_ji = &t.upward_connectivity[j][i - j - 1];
826                                for (c, col) in dc_ij.col_iter().enumerate() {
827                                    for value in col.iter_ref() {
828                                        assert!(uc_ji[*value].contains(&c));
829                                    }
830                                }
831                                for (k, uc_jik) in uc_ji.iter().enumerate() {
832                                    for value in uc_jik {
833                                        assert!(dc_ij.r().slice::<1>(1, *value).iter_ref()
834                                            .position(|&i| i == k).is_some());
835                                    }
836                                }
837                            }
838                        }
839                    }
840                }
841                #[test]
842                fn [< test_sub_entity_iter_ $cellname >]() {
843                    //! Test sub-entity iterators
844                    let t = [< example_topology_ $cellname >]();
845                    for cell_type in t.entity_types() {
846                        for index in 0..t.entity_count(*cell_type) {
847                            let cell = SingleTypeEntityTopology::new(&t, *cell_type, index);
848                            for dim in 0..reference_cell::dim(*cell_type) + 1 {
849                                let sub_entity_type = t.entity_types()[dim];
850                                for (i, j) in cell.sub_entity_iter(sub_entity_type).enumerate() {
851                                    assert_eq!(j, cell.sub_entity(sub_entity_type, i));
852                                }
853                            }
854                        }
855                    }
856                }
857            }
858        };
859    }
860
861    make_tests!(point);
862    make_tests!(interval);
863    make_tests!(triangle);
864    make_tests!(tetrahedron);
865
866    #[test]
867    fn test_orientation_triangle() {
868        let t =
869            SingleTypeTopology::new(&[0, 1, 2, 0, 2, 1], ReferenceCellType::Triangle, None, None);
870        assert_ne!(t.orientation[1][0], t.orientation[1][1]);
871    }
872
873    #[test]
874    fn test_orientation_quadrilateral() {
875        let t = SingleTypeTopology::new(
876            &[0, 1, 2, 3, 0, 3, 1, 2],
877            ReferenceCellType::Quadrilateral,
878            None,
879            None,
880        );
881        assert_ne!(t.orientation[1][0], t.orientation[1][1]);
882    }
883}
ndgrid/topology/single_type.rs

ndgrid/topology/
single_type.rs
