ndmesh/mesh/local_mesh/single_element/
builder.rs

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

ndmesh/mesh/local_mesh/single_element/
builder.rs
