ndgrid/grid/local_grid/mixed/
builder.rs

1//! Grid builder
2
3use super::MixedGrid;
4use crate::{
5    geometry::MixedGeometry,
6    topology::mixed::MixedTopology,
7    traits::{Builder, GeometryBuilder, GridBuilder, TopologyBuilder},
8    types::Scalar,
9};
10use itertools::izip;
11use ndelement::{
12    ciarlet::{CiarletElement, LagrangeElementFamily, lagrange},
13    map::IdentityMap,
14    reference_cell,
15    traits::FiniteElement,
16    types::{Continuity, ReferenceCellType},
17};
18use rlst::rlst_dynamic_array;
19use std::collections::{HashMap, HashSet};
20
21/// Grid builder for a grid with a mixture of element types
22#[derive(Debug)]
23pub struct MixedGridBuilder<T: Scalar> {
24    gdim: usize,
25    element_indices: HashMap<(ReferenceCellType, usize), usize>,
26    elements: Vec<CiarletElement<T, IdentityMap, T>>,
27    points_per_cell: Vec<usize>,
28    pub(crate) points: Vec<T>,
29    cells: Vec<usize>,
30    cell_starts: Vec<usize>,
31    cell_families: Vec<usize>,
32    cell_degrees: Vec<usize>,
33    cell_types: Vec<ReferenceCellType>,
34    point_indices_to_ids: Vec<usize>,
35    point_ids_to_indices: HashMap<usize, usize>,
36    cell_indices_to_ids: Vec<usize>,
37    cell_ids_to_indices: HashMap<usize, usize>,
38    point_indices: HashSet<usize>,
39    cell_indices: HashSet<usize>,
40}
41
42impl<T: Scalar> MixedGridBuilder<T> {
43    /// Create a new grid builder
44    pub fn new(gdim: usize) -> Self {
45        Self {
46            gdim,
47            element_indices: HashMap::new(),
48            elements: vec![],
49            points_per_cell: vec![],
50            points: vec![],
51            cells: vec![],
52            cell_starts: vec![],
53            cell_families: vec![],
54            cell_degrees: vec![],
55            cell_types: vec![],
56            point_indices_to_ids: vec![],
57            point_ids_to_indices: HashMap::new(),
58            cell_indices_to_ids: vec![],
59            cell_ids_to_indices: HashMap::new(),
60            point_indices: HashSet::new(),
61            cell_indices: HashSet::new(),
62        }
63    }
64}
65
66impl<T: Scalar> Builder for MixedGridBuilder<T> {
67    type Grid = MixedGrid<T, CiarletElement<T, IdentityMap, T>>;
68    type T = T;
69    type CellData<'a> = (ReferenceCellType, usize, &'a [usize]);
70    type EntityDescriptor = ReferenceCellType;
71
72    fn add_point(&mut self, id: usize, data: &[T]) {
73        if data.len() != self.gdim {
74            panic!("Point has wrong number of coordinates");
75        }
76        if self.point_indices.contains(&id) {
77            panic!("Cannot add point with duplicate id.");
78        }
79        self.point_ids_to_indices
80            .insert(id, self.point_indices_to_ids.len());
81        self.point_indices.insert(id);
82        self.point_indices_to_ids.push(id);
83        self.points.extend_from_slice(data);
84    }
85
86    fn add_cell(&mut self, id: usize, cell_data: (ReferenceCellType, usize, &[usize])) {
87        let element_index = *self
88            .element_indices
89            .entry((cell_data.0, cell_data.1))
90            .or_insert_with(|| {
91                let n = self.cell_indices_to_ids.len();
92                self.elements.push(lagrange::create::<T, T>(
93                    cell_data.0,
94                    cell_data.1,
95                    Continuity::Standard,
96                ));
97                self.points_per_cell.push(self.elements[n].dim());
98                n
99            });
100        self.cell_families.push(element_index);
101        self.cell_starts.push(self.cells.len());
102        self.cell_types.push(cell_data.0);
103        if reference_cell::dim(cell_data.0) != reference_cell::dim(self.cell_types[0]) {
104            panic!("Cannot create grid with cells of different topological dimensions");
105        }
106        self.cell_degrees.push(cell_data.1);
107        if self.cell_indices.contains(&id) {
108            panic!("Cannot add cell with duplicate id.");
109        }
110        assert_eq!(cell_data.2.len(), self.points_per_cell[element_index]);
111        self.cell_ids_to_indices
112            .insert(id, self.cell_indices_to_ids.len());
113        self.cell_indices.insert(id);
114        self.cell_indices_to_ids.push(id);
115        for id in cell_data.2 {
116            self.cells.push(self.point_ids_to_indices[id]);
117        }
118    }
119
120    fn add_cell_from_nodes_and_type(
121        &mut self,
122        id: usize,
123        nodes: &[usize],
124        cell_type: ReferenceCellType,
125        cell_degree: usize,
126    ) {
127        self.add_cell(id, (cell_type, cell_degree, nodes));
128    }
129
130    fn create_grid(&self) -> MixedGrid<T, CiarletElement<T, IdentityMap, T>> {
131        let cell_vertices =
132            self.extract_vertices(&self.cells, &self.cell_types, &self.cell_degrees);
133
134        // Add the vertex ids. But need to make sure that we don't have duplicates.
135        // So first generate a set of vertex ids, then convert it to a vector.
136        // We finally sort the vector to make sure that grid generation is reproducible across runs
137        // as sets do not have a stable ordering.
138        let mut vertex_ids = HashSet::<usize>::new();
139        for v in &cell_vertices {
140            vertex_ids.insert(*v);
141        }
142
143        let vertex_ids = {
144            let mut tmp = Vec::<usize>::with_capacity(vertex_ids.len());
145            tmp.extend(vertex_ids.iter());
146            tmp.sort();
147            tmp
148        };
149
150        let geometry = self.create_geometry(
151            &self.point_indices_to_ids,
152            &self.points,
153            &self.cells,
154            &self.cell_types,
155            &self.cell_degrees,
156        );
157
158        let topology = self.create_topology(
159            vertex_ids,
160            (0..self.cell_count()).collect::<Vec<_>>(),
161            &cell_vertices,
162            &self.cell_types,
163        );
164
165        self.create_grid_from_topology_geometry(topology, geometry)
166    }
167
168    fn point_count(&self) -> usize {
169        self.point_indices_to_ids.len()
170    }
171    fn cell_count(&self) -> usize {
172        self.cell_indices_to_ids.len()
173    }
174    fn point_indices_to_ids(&self) -> &[usize] {
175        &self.point_indices_to_ids
176    }
177    fn cell_indices_to_ids(&self) -> &[usize] {
178        unimplemented!();
179        // &self.cell_indices_to_ids
180    }
181    fn cell_points(&self, index: usize) -> &[usize] {
182        &self.cells[self.cell_starts[index]
183            ..self.cell_starts[index] + self.points_per_cell[self.cell_families[index]]]
184    }
185    fn cell_vertices(&self, index: usize) -> &[usize] {
186        &self.cells[self.cell_starts[index]
187            ..self.cell_starts[index] + reference_cell::entity_counts(self.cell_types[index])[0]]
188    }
189    fn point(&self, index: usize) -> &[T] {
190        &self.points[self.gdim * index..self.gdim * (index + 1)]
191    }
192    fn points(&self) -> &[T] {
193        &self.points
194    }
195    fn cell_type(&self, index: usize) -> ReferenceCellType {
196        self.cell_types[index]
197    }
198    fn cell_degree(&self, index: usize) -> usize {
199        self.cell_degrees[index]
200    }
201    fn gdim(&self) -> usize {
202        self.gdim
203    }
204    fn tdim(&self) -> usize {
205        reference_cell::dim(self.cell_types[0])
206    }
207    fn npts(&self, cell_type: Self::EntityDescriptor, degree: usize) -> usize {
208        self.points_per_cell[self.element_indices[&(cell_type, degree)]]
209    }
210}
211
212impl<T: Scalar> GeometryBuilder for MixedGridBuilder<T> {
213    type GridGeometry = MixedGeometry<T, CiarletElement<T, IdentityMap, T>>;
214    fn create_geometry(
215        &self,
216        point_ids: &[usize],
217        coordinates: &[Self::T],
218        cell_points: &[usize],
219        cell_types: &[ReferenceCellType],
220        cell_degrees: &[usize],
221    ) -> MixedGeometry<T, CiarletElement<T, IdentityMap, T>> {
222        let npts = point_ids.len();
223        let mut points = rlst_dynamic_array!(T, [self.gdim(), npts]);
224        points.data_mut().unwrap().copy_from_slice(coordinates);
225
226        let mut element_families = vec![];
227        let mut ef_indices = HashMap::new();
228        let mut cell_families = vec![];
229        for degree in cell_degrees {
230            cell_families.push(*ef_indices.entry(*degree).or_insert_with(|| {
231                let n = element_families.len();
232                element_families.push(LagrangeElementFamily::<T, T>::new(
233                    *degree,
234                    Continuity::Standard,
235                ));
236                n
237            }))
238        }
239
240        MixedGeometry::<T, CiarletElement<T, IdentityMap, T>>::new(
241            cell_types,
242            points,
243            cell_points,
244            &element_families,
245            &cell_families,
246        )
247    }
248}
249
250impl<T: Scalar> TopologyBuilder for MixedGridBuilder<T> {
251    type GridTopology = MixedTopology;
252    fn create_topology(
253        &self,
254        vertex_ids: Vec<usize>,
255        cell_ids: Vec<usize>,
256        cells: &[usize],
257        _cell_types: &[ReferenceCellType],
258    ) -> MixedTopology {
259        // Create a map from point ids to the corresponding positions in the points array.
260        let vertex_ids_to_pos = {
261            let mut tmp = HashMap::<usize, usize>::new();
262            for (i, id) in vertex_ids.iter().enumerate() {
263                tmp.insert(*id, i);
264            }
265            tmp
266        };
267
268        let cells = {
269            let mut new_cells = Vec::<usize>::with_capacity(cells.len());
270            for id in cells {
271                new_cells.push(vertex_ids_to_pos[id]);
272            }
273            new_cells
274        };
275
276        MixedTopology::new(&cells, &self.cell_types, Some(vertex_ids), Some(cell_ids))
277    }
278
279    fn extract_vertices(
280        &self,
281        cell_points: &[usize],
282        cell_types: &[Self::EntityDescriptor],
283        cell_degrees: &[usize],
284    ) -> Vec<usize> {
285        let mut cell_vertices = vec![];
286        let mut start = 0;
287        for (t, d) in izip!(cell_types, cell_degrees) {
288            let e = &self.elements[self.element_indices[&(*t, *d)]];
289            for v in 0..reference_cell::entity_counts(*t)[0] {
290                for d in e.entity_dofs(0, v).unwrap() {
291                    cell_vertices.push(cell_points[start + *d])
292                }
293            }
294            start += e.dim();
295        }
296        cell_vertices
297    }
298}
299
300impl<T: Scalar> GridBuilder for MixedGridBuilder<T> {
301    fn create_grid_from_topology_geometry(
302        &self,
303        topology: <Self as TopologyBuilder>::GridTopology,
304        geometry: <Self as GeometryBuilder>::GridGeometry,
305    ) -> <Self as Builder>::Grid {
306        MixedGrid::new(topology, geometry)
307    }
308}
309
310#[cfg(test)]
311mod test {
312    use super::*;
313
314    #[test]
315    #[should_panic]
316    fn test_duplicate_point_id() {
317        let mut b = MixedGridBuilder::<f64>::new(3);
318
319        b.add_point(2, &[0.0, 0.0, 0.0]);
320        b.add_point(0, &[1.0, 0.0, 0.0]);
321        b.add_point(1, &[0.0, 1.0, 0.0]);
322        b.add_point(2, &[1.0, 1.0, 0.0]);
323    }
324
325    #[test]
326    #[should_panic]
327    fn test_duplicate_cell_id() {
328        let mut b = MixedGridBuilder::<f64>::new(3);
329
330        b.add_point(0, &[0.0, 0.0, 0.0]);
331        b.add_point(1, &[1.0, 0.0, 0.0]);
332        b.add_point(2, &[0.0, 1.0, 0.0]);
333        b.add_point(3, &[1.0, 1.0, 0.0]);
334
335        b.add_cell(0, (ReferenceCellType::Triangle, 1, &[0, 1, 2]));
336        b.add_cell(0, (ReferenceCellType::Triangle, 1, &[1, 2, 3]));
337    }
338
339    #[test]
340    fn test_non_contiguous_ids() {
341        let mut b = MixedGridBuilder::<f64>::new(3);
342
343        b.add_point(0, &[0.0, 0.0, 0.0]);
344        b.add_point(1, &[1.0, 0.0, 0.0]);
345        b.add_point(2, &[0.0, 1.0, 0.0]);
346        b.add_point(4, &[1.0, 1.0, 0.0]);
347
348        b.add_cell(0, (ReferenceCellType::Triangle, 1, &[0, 1, 2]));
349        b.add_cell(2, (ReferenceCellType::Triangle, 1, &[1, 2, 4]));
350
351        b.create_grid();
352    }
353}
ndgrid/grid/local_grid/mixed/builder.rs

ndgrid/grid/local_grid/mixed/
builder.rs
