ndgrid/grid/local_grid/single_element/
grid.rs

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

ndgrid/grid/local_grid/single_element/
grid.rs
