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::{
31    Array, ValueArrayImpl,
32    dense::{base_array::BaseArray, data_container::VectorContainer},
33    rlst_dynamic_array,
34};
35
36/// Single element grid entity
37#[derive(Debug)]
38pub struct SingleElementGridEntity<
39    'a,
40    T: Scalar,
41    E: MappedFiniteElement<CellType = ReferenceCellType, T = T>,
42> {
43    grid: &'a SingleElementGrid<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    SingleElementGridEntity<'e, T, E>
51{
52    /// Create new
53    pub fn new(
54        grid: &'e SingleElementGrid<T, E>,
55        cell_index: usize,
56        entity_dim: usize,
57        entity_index: usize,
58    ) -> Self {
59        Self {
60            grid,
61            cell_index,
62            entity_dim,
63            entity_index,
64        }
65    }
66}
67impl<T: Scalar, E: MappedFiniteElement<CellType = ReferenceCellType, T = T>> Entity
68    for SingleElementGridEntity<'_, 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.grid.topology.entity_types()[self.entity_dim]
82    }
83    fn local_index(&self) -> usize {
84        self.grid
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.grid.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.grid.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.grid
107            .topology
108            .entity_id(self.entity_dim, self.local_index())
109    }
110}
111
112/// Single element grid entity iterator
113#[derive(Debug)]
114pub struct SingleElementGridEntityIter<
115    'a,
116    T: Scalar,
117    E: MappedFiniteElement<CellType = ReferenceCellType, T = T>,
118> {
119    grid: &'a SingleElementGrid<T, E>,
120    entity_type: ReferenceCellType,
121    index: usize,
122}
123
124impl<'a, T: Scalar, E: MappedFiniteElement<CellType = ReferenceCellType, T = T>>
125    SingleElementGridEntityIter<'a, T, E>
126{
127    /// Create new
128    pub fn new(grid: &'a SingleElementGrid<T, E>, entity_type: ReferenceCellType) -> Self {
129        Self {
130            grid,
131            entity_type,
132            index: 0,
133        }
134    }
135}
136impl<'a, T: Scalar, E: MappedFiniteElement<CellType = ReferenceCellType, T = T>> Iterator
137    for SingleElementGridEntityIter<'a, T, E>
138{
139    type Item = SingleElementGridEntity<'a, T, E>;
140
141    fn next(&mut self) -> Option<SingleElementGridEntity<'a, T, E>> {
142        self.index += 1;
143        self.grid.entity(self.entity_type, self.index - 1)
144    }
145}
146
147/// Serial single element grid
148#[derive(Debug)]
149pub struct SingleElementGrid<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 SerializableGrid<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 SingleElementGrid<T, CiarletElement<T, IdentityMap, T>>
169{
170    type SerializableType = SerializableGrid<T>;
171    fn to_serializable(&self) -> SerializableGrid<T> {
172        SerializableGrid {
173            topology: self.topology.to_serializable(),
174            geometry: self.geometry.to_serializable(),
175        }
176    }
177    fn from_serializable(s: SerializableGrid<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    SingleElementGrid<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> SingleElementGrid<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(geometry_degree, Continuity::Standard);
208
209        let geometry = SingleElementGeometry::<T, CiarletElement<T, IdentityMap, T>>::new(
210            cell_type, points, cells, &family,
211        );
212
213        let points_per_cell = family.element(cell_type).dim();
214        let tpoints_per_cell = reference_cell::entity_counts(cell_type)[0];
215        let ncells = cells.len() / points_per_cell;
216
217        let mut tcells = vec![0; tpoints_per_cell * ncells];
218        for c in 0..ncells {
219            tcells[c * tpoints_per_cell..(c + 1) * tpoints_per_cell].copy_from_slice(
220                &cells[c * points_per_cell..c * points_per_cell + tpoints_per_cell],
221            );
222        }
223
224        let topology = SingleTypeTopology::new(&tcells, cell_type, None, None);
225
226        Self { topology, geometry }
227    }
228}
229
230impl<T: Scalar, E: MappedFiniteElement<CellType = ReferenceCellType, T = T>> Grid
231    for SingleElementGrid<T, E>
232{
233    type T = T;
234    type Entity<'a>
235        = SingleElementGridEntity<'a, T, E>
236    where
237        Self: 'a;
238    type GeometryMap<'a>
239        = GeometryMap<'a, T, BaseArray<VectorContainer<T>, 2>, BaseArray<VectorContainer<usize>, 2>>
240    where
241        Self: 'a;
242    type EntityDescriptor = ReferenceCellType;
243    type EntityIter<'a>
244        = SingleElementGridEntityIter<'a, T, E>
245    where
246        Self: 'a;
247
248    fn geometry_dim(&self) -> usize {
249        self.geometry.dim()
250    }
251    fn topology_dim(&self) -> usize {
252        self.topology.dim()
253    }
254
255    fn entity(
256        &self,
257        entity_type: ReferenceCellType,
258        local_index: usize,
259    ) -> Option<Self::Entity<'_>> {
260        let dim = reference_cell::dim(entity_type);
261        if local_index < self.topology.entity_count(entity_type) {
262            if dim == self.topology_dim() {
263                Some(SingleElementGridEntity::new(self, local_index, dim, 0))
264            } else {
265                let cell = self.topology.upward_connectivity[dim][self.topology_dim() - dim - 1]
266                    [local_index][0];
267                let dc = &self.topology.downward_connectivity[self.topology_dim()][dim];
268                let index = (0..dc.shape()[0])
269                    .position(|i| dc[[i, cell]] == local_index)
270                    .unwrap();
271                Some(SingleElementGridEntity::new(self, cell, dim, index))
272            }
273        } else {
274            None
275        }
276    }
277
278    fn entity_types(&self, dim: usize) -> &[ReferenceCellType] {
279        &self.topology.entity_types()[dim..dim + 1]
280    }
281
282    fn entity_count(&self, entity_type: ReferenceCellType) -> usize {
283        self.topology.entity_count(entity_type)
284    }
285
286    fn entity_iter(&self, entity_type: ReferenceCellType) -> Self::EntityIter<'_> {
287        SingleElementGridEntityIter::new(self, entity_type)
288    }
289
290    fn entity_from_id(
291        &self,
292        entity_type: ReferenceCellType,
293        id: usize,
294    ) -> Option<Self::Entity<'_>> {
295        let entity_dim = reference_cell::dim(entity_type);
296        self.topology.ids_to_indices[entity_dim]
297            .get(&id)
298            .map(|i| self.entity(entity_type, *i))?
299    }
300
301    fn geometry_map<Array2Impl: ValueArrayImpl<T, 2>>(
302        &self,
303        entity_type: ReferenceCellType,
304        geometry_degree: usize,
305        points: &Array<Array2Impl, 2>,
306    ) -> GeometryMap<'_, T, BaseArray<VectorContainer<T>, 2>, BaseArray<VectorContainer<usize>, 2>>
307    {
308        debug_assert!(points.shape()[0] == reference_cell::dim(entity_type));
309        debug_assert!(geometry_degree == self.geometry.element().lagrange_superdegree());
310        if entity_type == self.topology.entity_types()[self.topology_dim()] {
311            GeometryMap::new(
312                self.geometry.element(),
313                points,
314                self.geometry.points(),
315                self.geometry.cells(),
316            )
317        } else {
318            unimplemented!();
319        }
320    }
321}
322
323#[cfg(feature = "mpi")]
324impl<T: Scalar + Equivalence, E: MappedFiniteElement<CellType = ReferenceCellType, T = T>>
325    DistributableGrid for SingleElementGrid<T, E>
326{
327    type ParallelGrid<'a, C: Communicator + 'a> =
328        ParallelGridImpl<'a, C, SingleElementGrid<T, CiarletElement<T, IdentityMap, T>>>;
329
330    fn distribute<'a, C: Communicator>(
331        &self,
332        comm: &'a C,
333        partitioner: GraphPartitioner,
334    ) -> Self::ParallelGrid<'a, C> {
335        let e = self.geometry.element();
336        let pts = self.geometry.points();
337        let cells = self.geometry.cells();
338        let mut b = SingleElementGridBuilder::<T>::new_with_capacity(
339            self.geometry.dim(),
340            pts.shape()[1],
341            cells.shape()[1],
342            (e.cell_type(), e.lagrange_superdegree()),
343        );
344        for p in 0..pts.shape()[1] {
345            b.add_point(
346                p,
347                &pts.data().unwrap()[p * pts.shape()[0]..(p + 1) * pts.shape()[0]],
348            );
349        }
350        for c in 0..cells.shape()[1] {
351            b.add_cell(
352                c,
353                &cells.data().unwrap()[c * cells.shape()[0]..(c + 1) * cells.shape()[0]],
354            );
355        }
356        b.create_parallel_grid_root(comm, partitioner)
357    }
358}
359
360#[cfg(test)]
361mod test {
362    use super::*;
363    use crate::traits::Topology;
364    use itertools::izip;
365    use ndelement::{
366        ciarlet::{CiarletElement, LagrangeElementFamily},
367        reference_cell,
368        types::Continuity,
369    };
370    use rlst::rlst_dynamic_array;
371
372    fn example_grid_triangle() -> SingleElementGrid<f64, CiarletElement<f64, IdentityMap, f64>> {
373        let mut points = rlst_dynamic_array!(f64, [3, 4]);
374        *points.get_mut([0, 0]).unwrap() = 0.0;
375        *points.get_mut([1, 0]).unwrap() = 0.0;
376        *points.get_mut([2, 0]).unwrap() = 1.0;
377        *points.get_mut([0, 1]).unwrap() = 1.0;
378        *points.get_mut([1, 1]).unwrap() = 0.0;
379        *points.get_mut([2, 1]).unwrap() = 0.0;
380        *points.get_mut([0, 2]).unwrap() = 0.0;
381        *points.get_mut([1, 2]).unwrap() = 1.0;
382        *points.get_mut([2, 2]).unwrap() = 0.0;
383        *points.get_mut([0, 3]).unwrap() = 2.0;
384        *points.get_mut([1, 3]).unwrap() = 1.0;
385        *points.get_mut([2, 3]).unwrap() = 0.0;
386        let family = LagrangeElementFamily::<f64>::new(1, Continuity::Standard);
387        SingleElementGrid::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 grid = example_grid_triangle();
401        let conn = reference_cell::connectivity(ReferenceCellType::Triangle);
402        for edge in grid.entity_iter(ReferenceCellType::Interval) {
403            let cell = grid
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}
ndgrid/grid/local_grid/single_element/grid.rs

ndgrid/grid/local_grid/single_element/
grid.rs
