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

ndgrid/grid/local_grid/single_element/
builder.rs
