ndmesh/mesh/local_mesh/single_element/
mesh.rs

1//! Single element mesh
2#[cfg(feature = "mpi")]
3use crate::ParallelMeshImpl;
4#[cfg(feature = "mpi")]
5use crate::{
6    SingleElementMeshBuilder,
7    traits::{Builder, DistributableMesh, ParallelBuilder},
8    types::GraphPartitioner,
9};
10#[cfg(feature = "serde")]
11use crate::{
12    geometry::single_element::SerializableGeometry, topology::single_type::SerializableTopology,
13    traits::ConvertToSerializable,
14};
15use crate::{
16    geometry::{GeometryMap, SingleElementEntityGeometry, SingleElementGeometry},
17    topology::single_type::{SingleTypeEntityTopology, SingleTypeTopology},
18    traits::{Entity, Mesh},
19    types::{Ownership, Scalar},
20};
21#[cfg(feature = "mpi")]
22use mpi::traits::{Communicator, Equivalence};
23use ndelement::{
24    ciarlet::{CiarletElement, LagrangeElementFamily, LagrangeVariant},
25    map::IdentityMap,
26    reference_cell,
27    traits::{ElementFamily, FiniteElement, MappedFiniteElement},
28    types::{Continuity, ReferenceCellType},
29};
30use rlst::{
31    Array, ValueArrayImpl,
32    dense::{base_array::BaseArray, data_container::VectorContainer},
33    rlst_dynamic_array,
34};
35
36/// Single element mesh entity
37#[derive(Debug)]
38pub struct SingleElementMeshEntity<
39    'a,
40    T: Scalar,
41    E: MappedFiniteElement<CellType = ReferenceCellType, T = T>,
42> {
43    mesh: &'a SingleElementMesh<T, E>,
44    cell_index: usize,
45    entity_dim: usize,
46    entity_index: usize,
47}
48
49impl<'e, T: Scalar, E: MappedFiniteElement<CellType = ReferenceCellType, T = T>>
50    SingleElementMeshEntity<'e, T, E>
51{
52    /// Create new
53    pub fn new(
54        mesh: &'e SingleElementMesh<T, E>,
55        cell_index: usize,
56        entity_dim: usize,
57        entity_index: usize,
58    ) -> Self {
59        Self {
60            mesh,
61            cell_index,
62            entity_dim,
63            entity_index,
64        }
65    }
66}
67impl<T: Scalar, E: MappedFiniteElement<CellType = ReferenceCellType, T = T>> Entity
68    for SingleElementMeshEntity<'_, T, E>
69{
70    type T = T;
71    type EntityDescriptor = ReferenceCellType;
72    type Topology<'a>
73        = SingleTypeEntityTopology<'a>
74    where
75        Self: 'a;
76    type Geometry<'a>
77        = SingleElementEntityGeometry<'a, T, E>
78    where
79        Self: 'a;
80    fn entity_type(&self) -> ReferenceCellType {
81        self.mesh.topology.entity_types()[self.entity_dim]
82    }
83    fn local_index(&self) -> usize {
84        self.mesh
85            .topology
86            .cell_entity_index(self.cell_index, self.entity_dim, self.entity_index)
87    }
88    fn global_index(&self) -> usize {
89        self.local_index()
90    }
91    fn geometry(&self) -> Self::Geometry<'_> {
92        SingleElementEntityGeometry::new(
93            &self.mesh.geometry,
94            self.cell_index,
95            self.entity_dim,
96            self.entity_index,
97        )
98    }
99    fn topology(&self) -> Self::Topology<'_> {
100        SingleTypeEntityTopology::new(&self.mesh.topology, self.entity_type(), self.local_index())
101    }
102    fn ownership(&self) -> Ownership {
103        Ownership::Owned
104    }
105    fn id(&self) -> Option<usize> {
106        self.mesh
107            .topology
108            .entity_id(self.entity_dim, self.local_index())
109    }
110}
111
112/// Single element mesh entity iterator
113#[derive(Debug)]
114pub struct SingleElementMeshEntityIter<
115    'a,
116    T: Scalar,
117    E: MappedFiniteElement<CellType = ReferenceCellType, T = T>,
118> {
119    mesh: &'a SingleElementMesh<T, E>,
120    entity_type: ReferenceCellType,
121    index: usize,
122}
123
124impl<'a, T: Scalar, E: MappedFiniteElement<CellType = ReferenceCellType, T = T>>
125    SingleElementMeshEntityIter<'a, T, E>
126{
127    /// Create new
128    pub fn new(mesh: &'a SingleElementMesh<T, E>, entity_type: ReferenceCellType) -> Self {
129        Self {
130            mesh,
131            entity_type,
132            index: 0,
133        }
134    }
135}
136impl<'a, T: Scalar, E: MappedFiniteElement<CellType = ReferenceCellType, T = T>> Iterator
137    for SingleElementMeshEntityIter<'a, T, E>
138{
139    type Item = SingleElementMeshEntity<'a, T, E>;
140
141    fn next(&mut self) -> Option<SingleElementMeshEntity<'a, T, E>> {
142        self.index += 1;
143        self.mesh.entity(self.entity_type, self.index - 1)
144    }
145}
146
147/// Serial single element mesh
148#[derive(Debug)]
149pub struct SingleElementMesh<T: Scalar, E: MappedFiniteElement<CellType = ReferenceCellType, T = T>>
150{
151    topology: SingleTypeTopology,
152    geometry: SingleElementGeometry<T, E>,
153}
154
155#[cfg(feature = "serde")]
156#[derive(serde::Serialize, Debug, serde::Deserialize)]
157#[serde(bound = "for<'de2> T: serde::Deserialize<'de2>")]
158pub struct SerializableMesh<T: Scalar + serde::Serialize>
159where
160    for<'de2> T: serde::Deserialize<'de2>,
161{
162    topology: SerializableTopology,
163    geometry: SerializableGeometry<T>,
164}
165
166#[cfg(feature = "serde")]
167impl<T: Scalar + serde::Serialize> ConvertToSerializable
168    for SingleElementMesh<T, CiarletElement<T, IdentityMap, T>>
169{
170    type SerializableType = SerializableMesh<T>;
171    fn to_serializable(&self) -> SerializableMesh<T> {
172        SerializableMesh {
173            topology: self.topology.to_serializable(),
174            geometry: self.geometry.to_serializable(),
175        }
176    }
177    fn from_serializable(s: SerializableMesh<T>) -> Self {
178        Self {
179            topology: SingleTypeTopology::from_serializable(s.topology),
180            geometry: SingleElementGeometry::from_serializable(s.geometry),
181        }
182    }
183}
184
185impl<T: Scalar, E: MappedFiniteElement<CellType = ReferenceCellType, T = T>>
186    SingleElementMesh<T, E>
187{
188    /// Create new
189    pub fn new(topology: SingleTypeTopology, geometry: SingleElementGeometry<T, E>) -> Self {
190        Self { topology, geometry }
191    }
192}
193
194impl<T: Scalar> SingleElementMesh<T, CiarletElement<T, IdentityMap, T>> {
195    /// Create new from raw data
196    pub fn new_from_raw_data(
197        coordinates: &[T],
198        gdim: usize,
199        cells: &[usize],
200        cell_type: ReferenceCellType,
201        geometry_degree: usize,
202    ) -> Self {
203        let npts = coordinates.len() / gdim;
204        let mut points = rlst_dynamic_array!(T, [gdim, npts]);
205        points.data_mut().unwrap().copy_from_slice(coordinates);
206
207        let family = LagrangeElementFamily::<T, T>::new(
208            geometry_degree,
209            Continuity::Standard,
210            LagrangeVariant::Equispaced,
211        );
212
213        let geometry = SingleElementGeometry::<T, CiarletElement<T, IdentityMap, T>>::new(
214            cell_type, points, cells, &family,
215        );
216
217        let points_per_cell = family.element(cell_type).dim();
218        let tpoints_per_cell = reference_cell::entity_counts(cell_type)[0];
219        let ncells = cells.len() / points_per_cell;
220
221        let mut tcells = vec![0; tpoints_per_cell * ncells];
222        for c in 0..ncells {
223            tcells[c * tpoints_per_cell..(c + 1) * tpoints_per_cell].copy_from_slice(
224                &cells[c * points_per_cell..c * points_per_cell + tpoints_per_cell],
225            );
226        }
227
228        let topology = SingleTypeTopology::new(&tcells, cell_type, None, None);
229
230        Self { topology, geometry }
231    }
232}
233
234impl<T: Scalar, E: MappedFiniteElement<CellType = ReferenceCellType, T = T>> Mesh
235    for SingleElementMesh<T, E>
236{
237    type T = T;
238    type Entity<'a>
239        = SingleElementMeshEntity<'a, T, E>
240    where
241        Self: 'a;
242    type GeometryMap<'a>
243        = GeometryMap<'a, T, BaseArray<VectorContainer<T>, 2>, BaseArray<VectorContainer<usize>, 2>>
244    where
245        Self: 'a;
246    type EntityDescriptor = ReferenceCellType;
247    type EntityIter<'a>
248        = SingleElementMeshEntityIter<'a, T, E>
249    where
250        Self: 'a;
251
252    fn geometry_dim(&self) -> usize {
253        self.geometry.dim()
254    }
255    fn topology_dim(&self) -> usize {
256        self.topology.dim()
257    }
258
259    fn entity(
260        &self,
261        entity_type: ReferenceCellType,
262        local_index: usize,
263    ) -> Option<Self::Entity<'_>> {
264        let dim = reference_cell::dim(entity_type);
265        if local_index < self.topology.entity_count(entity_type) {
266            if dim == self.topology_dim() {
267                Some(SingleElementMeshEntity::new(self, local_index, dim, 0))
268            } else {
269                let cell = self.topology.upward_connectivity[dim][self.topology_dim() - dim - 1]
270                    [local_index][0];
271                let dc = &self.topology.downward_connectivity[self.topology_dim()][dim];
272                let index = (0..dc.shape()[0])
273                    .position(|i| dc[[i, cell]] == local_index)
274                    .unwrap();
275                Some(SingleElementMeshEntity::new(self, cell, dim, index))
276            }
277        } else {
278            None
279        }
280    }
281
282    fn entity_types(&self, dim: usize) -> &[ReferenceCellType] {
283        &self.topology.entity_types()[dim..dim + 1]
284    }
285
286    fn entity_count(&self, entity_type: ReferenceCellType) -> usize {
287        self.topology.entity_count(entity_type)
288    }
289
290    fn entity_iter(&self, entity_type: ReferenceCellType) -> Self::EntityIter<'_> {
291        SingleElementMeshEntityIter::new(self, entity_type)
292    }
293
294    fn entity_from_id(
295        &self,
296        entity_type: ReferenceCellType,
297        id: usize,
298    ) -> Option<Self::Entity<'_>> {
299        let entity_dim = reference_cell::dim(entity_type);
300        self.topology.ids_to_indices[entity_dim]
301            .get(&id)
302            .map(|i| self.entity(entity_type, *i))?
303    }
304
305    fn geometry_map<Array2Impl: ValueArrayImpl<T, 2>>(
306        &self,
307        entity_type: ReferenceCellType,
308        geometry_degree: usize,
309        points: &Array<Array2Impl, 2>,
310    ) -> GeometryMap<'_, T, BaseArray<VectorContainer<T>, 2>, BaseArray<VectorContainer<usize>, 2>>
311    {
312        debug_assert!(points.shape()[0] == reference_cell::dim(entity_type));
313        debug_assert!(geometry_degree == self.geometry.element().lagrange_superdegree());
314        if entity_type == self.topology.entity_types()[self.topology_dim()] {
315            GeometryMap::new(
316                self.geometry.element(),
317                points,
318                self.geometry.points(),
319                self.geometry.cells(),
320            )
321        } else {
322            unimplemented!();
323        }
324    }
325}
326
327#[cfg(feature = "mpi")]
328impl<T: Scalar + Equivalence, E: MappedFiniteElement<CellType = ReferenceCellType, T = T>>
329    DistributableMesh for SingleElementMesh<T, E>
330{
331    type ParallelMesh<'a, C: Communicator + 'a> =
332        ParallelMeshImpl<'a, C, SingleElementMesh<T, CiarletElement<T, IdentityMap, T>>>;
333
334    fn distribute<'a, C: Communicator>(
335        &self,
336        comm: &'a C,
337        partitioner: GraphPartitioner,
338    ) -> Self::ParallelMesh<'a, C> {
339        let e = self.geometry.element();
340        let pts = self.geometry.points();
341        let cells = self.geometry.cells();
342        let mut b = SingleElementMeshBuilder::<T>::new_with_capacity(
343            self.geometry.dim(),
344            pts.shape()[1],
345            cells.shape()[1],
346            (e.cell_type(), e.lagrange_superdegree()),
347        );
348        for p in 0..pts.shape()[1] {
349            b.add_point(
350                p,
351                &pts.data().unwrap()[p * pts.shape()[0]..(p + 1) * pts.shape()[0]],
352            );
353        }
354        for c in 0..cells.shape()[1] {
355            b.add_cell(
356                c,
357                &cells.data().unwrap()[c * cells.shape()[0]..(c + 1) * cells.shape()[0]],
358            );
359        }
360        b.create_parallel_mesh_root(comm, partitioner)
361    }
362}
363
364#[cfg(test)]
365mod test {
366    use super::*;
367    use crate::traits::Topology;
368    use itertools::izip;
369    use rlst::rlst_dynamic_array;
370
371    fn example_mesh_triangle() -> SingleElementMesh<f64, CiarletElement<f64, IdentityMap, f64>> {
372        let mut points = rlst_dynamic_array!(f64, [3, 4]);
373        *points.get_mut([0, 0]).unwrap() = 0.0;
374        *points.get_mut([1, 0]).unwrap() = 0.0;
375        *points.get_mut([2, 0]).unwrap() = 1.0;
376        *points.get_mut([0, 1]).unwrap() = 1.0;
377        *points.get_mut([1, 1]).unwrap() = 0.0;
378        *points.get_mut([2, 1]).unwrap() = 0.0;
379        *points.get_mut([0, 2]).unwrap() = 0.0;
380        *points.get_mut([1, 2]).unwrap() = 1.0;
381        *points.get_mut([2, 2]).unwrap() = 0.0;
382        *points.get_mut([0, 3]).unwrap() = 2.0;
383        *points.get_mut([1, 3]).unwrap() = 1.0;
384        *points.get_mut([2, 3]).unwrap() = 0.0;
385        let family =
386            LagrangeElementFamily::<f64>::new(1, Continuity::Standard, LagrangeVariant::Equispaced);
387        SingleElementMesh::new(
388            SingleTypeTopology::new(&[0, 1, 2, 2, 1, 3], ReferenceCellType::Triangle, None, None),
389            SingleElementGeometry::<f64, CiarletElement<f64, IdentityMap, f64>>::new(
390                ReferenceCellType::Triangle,
391                points,
392                &[0, 1, 2, 2, 1, 3],
393                &family,
394            ),
395        )
396    }
397
398    #[test]
399    fn test_edges_triangle() {
400        let mesh = example_mesh_triangle();
401        let conn = reference_cell::connectivity(ReferenceCellType::Triangle);
402        for edge in mesh.entity_iter(ReferenceCellType::Interval) {
403            let cell = mesh
404                .entity(ReferenceCellType::Triangle, edge.cell_index)
405                .unwrap();
406            for (i, v) in izip!(
407                &conn[1][edge.entity_index][0],
408                edge.topology().sub_entity_iter(ReferenceCellType::Point)
409            ) {
410                assert_eq!(v, cell.topology().sub_entity(ReferenceCellType::Point, *i));
411            }
412        }
413    }
414}
ndmesh/mesh/local_mesh/single_element/mesh.rs

ndmesh/mesh/local_mesh/single_element/
mesh.rs
