ndgrid/grid/local_grid/single_element/
builder.rs

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

ndgrid/grid/local_grid/single_element/
builder.rs
