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

ndgrid/grid/local_grid/mixed/
builder.rs
