ndmesh/mesh/local_mesh/mixed/
builder.rs

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

ndmesh/mesh/local_mesh/mixed/
builder.rs
