ndgrid/grid/local_grid/mixed/
grid.rs

1//! Mixed grid
2#[cfg(feature = "mpi")]
3use crate::ParallelGridImpl;
4#[cfg(feature = "mpi")]
5use crate::{
6    MixedGridBuilder,
7    traits::{Builder, DistributableGrid, ParallelBuilder},
8    types::GraphPartitioner,
9};
10#[cfg(feature = "serde")]
11use crate::{
12    geometry::mixed::SerializableGeometry, topology::mixed::SerializableTopology,
13    traits::ConvertToSerializable,
14};
15use crate::{
16    geometry::{GeometryMap, MixedEntityGeometry, MixedGeometry},
17    topology::mixed::{MixedEntityTopology, MixedTopology},
18    traits::{Entity, Grid},
19    types::{Ownership, Scalar},
20};
21use itertools::izip;
22#[cfg(feature = "mpi")]
23use mpi::traits::{Communicator, Equivalence};
24use ndelement::{
25    ciarlet::{CiarletElement, LagrangeElementFamily},
26    map::IdentityMap,
27    reference_cell,
28    traits::{ElementFamily, FiniteElement, MappedFiniteElement},
29    types::{Continuity, ReferenceCellType},
30};
31use rlst::{
32    Array, ValueArrayImpl,
33    dense::{base_array::BaseArray, data_container::VectorContainer},
34    rlst_dynamic_array,
35};
36use std::collections::HashMap;
37
38/// Mixed grid entity
39#[derive(Debug)]
40pub struct MixedGridEntity<
41    'a,
42    T: Scalar,
43    E: MappedFiniteElement<CellType = ReferenceCellType, T = T>,
44> {
45    grid: &'a MixedGrid<T, E>,
46    cell_type: ReferenceCellType,
47    cell_index: usize,
48    entity_type: ReferenceCellType,
49    entity_index: usize,
50    geometry_element_index: usize,
51    geometry_cell_index: usize,
52}
53
54impl<'e, T: Scalar, E: MappedFiniteElement<CellType = ReferenceCellType, T = T>>
55    MixedGridEntity<'e, T, E>
56{
57    /// Create new
58    pub fn new(
59        grid: &'e MixedGrid<T, E>,
60        cell_type: ReferenceCellType,
61        cell_index: usize,
62        entity_type: ReferenceCellType,
63        entity_index: usize,
64    ) -> Self {
65        Self {
66            grid,
67            cell_type,
68            cell_index,
69            entity_type,
70            entity_index,
71            geometry_element_index: grid.geometry.insertion_indices_to_element_indices
72                [grid.topology.insertion_indices[&cell_type][cell_index]],
73            geometry_cell_index: grid.geometry.insertion_indices_to_cell_indices
74                [grid.topology.insertion_indices[&cell_type][cell_index]],
75        }
76    }
77}
78impl<T: Scalar, E: MappedFiniteElement<CellType = ReferenceCellType, T = T>> Entity
79    for MixedGridEntity<'_, T, E>
80{
81    type T = T;
82    type EntityDescriptor = ReferenceCellType;
83    type Topology<'a>
84        = MixedEntityTopology<'a>
85    where
86        Self: 'a;
87    type Geometry<'a>
88        = MixedEntityGeometry<'a, T, E>
89    where
90        Self: 'a;
91    fn entity_type(&self) -> ReferenceCellType {
92        self.entity_type
93    }
94    fn local_index(&self) -> usize {
95        self.grid.topology.cell_entity_index(
96            self.cell_type,
97            self.cell_index,
98            self.entity_type,
99            self.entity_index,
100        )
101    }
102    fn global_index(&self) -> usize {
103        self.local_index()
104    }
105    fn geometry(&self) -> Self::Geometry<'_> {
106        MixedEntityGeometry::new(
107            &self.grid.geometry,
108            self.geometry_element_index,
109            self.geometry_cell_index,
110            reference_cell::dim(self.entity_type),
111            self.entity_index,
112        )
113    }
114    fn topology(&self) -> Self::Topology<'_> {
115        MixedEntityTopology::new(&self.grid.topology, self.entity_type(), self.local_index())
116    }
117    fn ownership(&self) -> Ownership {
118        Ownership::Owned
119    }
120    fn id(&self) -> Option<usize> {
121        self.grid
122            .topology
123            .entity_id(self.entity_type, self.local_index())
124    }
125}
126
127/// Mixed grid entity iterator
128#[derive(Debug)]
129pub struct MixedGridEntityIter<
130    'a,
131    T: Scalar,
132    E: MappedFiniteElement<CellType = ReferenceCellType, T = T>,
133> {
134    grid: &'a MixedGrid<T, E>,
135    entity_type: ReferenceCellType,
136    index: usize,
137}
138
139impl<'a, T: Scalar, E: MappedFiniteElement<CellType = ReferenceCellType, T = T>>
140    MixedGridEntityIter<'a, T, E>
141{
142    /// Create new
143    pub fn new(grid: &'a MixedGrid<T, E>, entity_type: ReferenceCellType) -> Self {
144        Self {
145            grid,
146            entity_type,
147            index: 0,
148        }
149    }
150}
151impl<'a, T: Scalar, E: MappedFiniteElement<CellType = ReferenceCellType, T = T>> Iterator
152    for MixedGridEntityIter<'a, T, E>
153{
154    type Item = MixedGridEntity<'a, T, E>;
155
156    fn next(&mut self) -> Option<MixedGridEntity<'a, T, E>> {
157        self.index += 1;
158        self.grid.entity(self.entity_type, self.index - 1)
159    }
160}
161
162/// Serial mixed element grid
163#[derive(Debug)]
164pub struct MixedGrid<T: Scalar, E: MappedFiniteElement<CellType = ReferenceCellType, T = T>> {
165    topology: MixedTopology,
166    geometry: MixedGeometry<T, E>,
167}
168
169#[cfg(feature = "serde")]
170#[derive(serde::Serialize, Debug, serde::Deserialize)]
171#[serde(bound = "for<'de2> T: serde::Deserialize<'de2>")]
172pub struct SerializableGrid<T: Scalar + serde::Serialize>
173where
174    for<'de2> T: serde::Deserialize<'de2>,
175{
176    topology: SerializableTopology,
177    geometry: SerializableGeometry<T>,
178}
179
180#[cfg(feature = "serde")]
181impl<T: Scalar + serde::Serialize> ConvertToSerializable
182    for MixedGrid<T, CiarletElement<T, IdentityMap, T>>
183{
184    type SerializableType = SerializableGrid<T>;
185    fn to_serializable(&self) -> SerializableGrid<T> {
186        SerializableGrid {
187            topology: self.topology.to_serializable(),
188            geometry: self.geometry.to_serializable(),
189        }
190    }
191    fn from_serializable(s: SerializableGrid<T>) -> Self {
192        Self {
193            topology: MixedTopology::from_serializable(s.topology),
194            geometry: MixedGeometry::from_serializable(s.geometry),
195        }
196    }
197}
198
199impl<T: Scalar, E: MappedFiniteElement<CellType = ReferenceCellType, T = T>> MixedGrid<T, E> {
200    /// Create new
201    pub fn new(topology: MixedTopology, geometry: MixedGeometry<T, E>) -> Self {
202        Self { topology, geometry }
203    }
204}
205
206impl<T: Scalar> MixedGrid<T, CiarletElement<T, IdentityMap, T>> {
207    /// Create new from raw data
208    pub fn new_from_raw_data(
209        coordinates: &[T],
210        gdim: usize,
211        cells: &[usize],
212        cell_types: &[ReferenceCellType],
213        cell_degrees: &[usize],
214    ) -> Self {
215        let npts = coordinates.len() / gdim;
216        let mut points = rlst_dynamic_array!(T, [gdim, npts]);
217        points.data_mut().unwrap().copy_from_slice(coordinates);
218
219        let mut element_families = vec![];
220        let mut element_family_indices = HashMap::new();
221
222        let cell_families = cell_degrees
223            .iter()
224            .map(|d| {
225                *element_family_indices.entry(*d).or_insert_with(|| {
226                    let index = element_families.len();
227                    element_families
228                        .push(LagrangeElementFamily::<T, T>::new(*d, Continuity::Standard));
229                    index
230                })
231            })
232            .collect::<Vec<_>>();
233
234        let geometry = MixedGeometry::<T, CiarletElement<T, IdentityMap, T>>::new(
235            cell_types,
236            points,
237            cells,
238            &element_families,
239            &cell_families,
240        );
241
242        let mut start = 0;
243        let mut tcells = vec![];
244        for (t, f) in izip!(cell_types, &cell_families) {
245            let tpoints_per_cell = reference_cell::entity_counts(*t)[0];
246            for i in 0..tpoints_per_cell {
247                tcells.push(cells[start + i]);
248            }
249            start += element_families[*f].element(*t).dim();
250        }
251
252        let topology = MixedTopology::new(&tcells, cell_types, None, None);
253
254        Self { topology, geometry }
255    }
256}
257
258impl<T: Scalar, E: MappedFiniteElement<CellType = ReferenceCellType, T = T>> Grid
259    for MixedGrid<T, E>
260{
261    type T = T;
262    type Entity<'a>
263        = MixedGridEntity<'a, T, E>
264    where
265        Self: 'a;
266    type GeometryMap<'a>
267        = GeometryMap<'a, T, BaseArray<VectorContainer<T>, 2>, BaseArray<VectorContainer<usize>, 2>>
268    where
269        Self: 'a;
270    type EntityDescriptor = ReferenceCellType;
271    type EntityIter<'a>
272        = MixedGridEntityIter<'a, T, E>
273    where
274        Self: 'a;
275
276    fn geometry_dim(&self) -> usize {
277        self.geometry.dim()
278    }
279    fn topology_dim(&self) -> usize {
280        self.topology.dim()
281    }
282
283    fn entity(
284        &self,
285        entity_type: ReferenceCellType,
286        local_index: usize,
287    ) -> Option<Self::Entity<'_>> {
288        let dim = reference_cell::dim(entity_type);
289        if local_index < self.topology.entity_count(entity_type) {
290            if dim == self.topology_dim() {
291                Some(MixedGridEntity::new(
292                    self,
293                    entity_type,
294                    local_index,
295                    entity_type,
296                    0,
297                ))
298            } else {
299                for t in &self.topology.entity_types()[self.topology_dim()] {
300                    if let Some(cell) =
301                        self.topology.upward_connectivity[&entity_type][t][local_index].first()
302                    {
303                        let dc = &self.topology.downward_connectivity[t][&entity_type];
304                        if let Some(index) =
305                            (0..dc.shape()[0]).position(|i| dc[[i, *cell]] == local_index)
306                        {
307                            return Some(MixedGridEntity::new(self, *t, *cell, entity_type, index));
308                        }
309                    }
310                }
311                None
312            }
313        } else {
314            None
315        }
316    }
317
318    fn entity_types(&self, dim: usize) -> &[ReferenceCellType] {
319        &self.topology.entity_types()[dim]
320    }
321
322    fn entity_count(&self, entity_type: ReferenceCellType) -> usize {
323        self.topology.entity_count(entity_type)
324    }
325
326    fn entity_iter(&self, entity_type: ReferenceCellType) -> Self::EntityIter<'_> {
327        MixedGridEntityIter::new(self, entity_type)
328    }
329
330    fn entity_from_id(
331        &self,
332        entity_type: ReferenceCellType,
333        id: usize,
334    ) -> Option<Self::Entity<'_>> {
335        self.topology.ids_to_indices[&entity_type]
336            .get(&id)
337            .map(|i| self.entity(entity_type, *i))?
338    }
339
340    fn geometry_map<Array2Impl: ValueArrayImpl<T, 2>>(
341        &self,
342        entity_type: ReferenceCellType,
343        geometry_degree: usize,
344        points: &Array<Array2Impl, 2>,
345    ) -> GeometryMap<'_, T, BaseArray<VectorContainer<T>, 2>, BaseArray<VectorContainer<usize>, 2>>
346    {
347        debug_assert!(points.shape()[0] == reference_cell::dim(entity_type));
348
349        for i in 0..self.geometry.element_count() {
350            let e = self.geometry.element(i);
351            if e.cell_type() == entity_type && e.lagrange_superdegree() == geometry_degree {
352                return GeometryMap::new(e, points, self.geometry.points(), self.geometry.cells(i));
353            }
354        }
355        unimplemented!();
356    }
357}
358
359#[cfg(feature = "mpi")]
360impl<T: Scalar + Equivalence, E: MappedFiniteElement<CellType = ReferenceCellType, T = T>>
361    DistributableGrid for MixedGrid<T, E>
362{
363    type ParallelGrid<'a, C: Communicator + 'a> =
364        ParallelGridImpl<'a, C, MixedGrid<T, CiarletElement<T, IdentityMap, T>>>;
365
366    fn distribute<'a, C: Communicator>(
367        &self,
368        comm: &'a C,
369        partitioner: GraphPartitioner,
370    ) -> Self::ParallelGrid<'a, C> {
371        let mut b = MixedGridBuilder::<T>::new(self.geometry.dim());
372        let pts = self.geometry.points();
373        for p in 0..pts.shape()[1] {
374            b.add_point(
375                p,
376                &pts.data().unwrap()[p * pts.shape()[0]..(p + 1) * pts.shape()[0]],
377            );
378        }
379
380        for (c, (element_i, cell_i)) in izip!(
381            &self.geometry.insertion_indices_to_element_indices,
382            &self.geometry.insertion_indices_to_cell_indices
383        )
384        .enumerate()
385        {
386            let e = self.geometry.element(*element_i);
387            let cells = self.geometry.cells(*element_i);
388            b.add_cell(
389                c,
390                (
391                    e.cell_type(),
392                    e.lagrange_superdegree(),
393                    &cells.data().unwrap()
394                        [cell_i * cells.shape()[0]..(cell_i + 1) * cells.shape()[0]],
395                ),
396            );
397        }
398        b.create_parallel_grid_root(comm, partitioner)
399    }
400}
401
402#[cfg(test)]
403mod test {
404    use super::*;
405    use crate::traits::{GeometryMap, Topology};
406    use approx::*;
407    use itertools::izip;
408    use ndelement::{
409        ciarlet::{CiarletElement, LagrangeElementFamily},
410        reference_cell,
411        types::Continuity,
412    };
413    use rlst::rlst_dynamic_array;
414
415    fn example_grid_triangle() -> MixedGrid<f64, CiarletElement<f64, IdentityMap, f64>> {
416        let mut points = rlst_dynamic_array!(f64, [3, 4]);
417        *points.get_mut([0, 0]).unwrap() = 0.0;
418        *points.get_mut([1, 0]).unwrap() = 0.0;
419        *points.get_mut([2, 0]).unwrap() = 1.0;
420        *points.get_mut([0, 1]).unwrap() = 1.0;
421        *points.get_mut([1, 1]).unwrap() = 0.0;
422        *points.get_mut([2, 1]).unwrap() = 0.0;
423        *points.get_mut([0, 2]).unwrap() = 0.0;
424        *points.get_mut([1, 2]).unwrap() = 1.0;
425        *points.get_mut([2, 2]).unwrap() = 0.0;
426        *points.get_mut([0, 3]).unwrap() = 2.0;
427        *points.get_mut([1, 3]).unwrap() = 1.0;
428        *points.get_mut([2, 3]).unwrap() = 0.0;
429        let family = LagrangeElementFamily::<f64>::new(1, Continuity::Standard);
430        MixedGrid::new(
431            MixedTopology::new(
432                &[0, 1, 2, 2, 1, 3],
433                &[ReferenceCellType::Triangle, ReferenceCellType::Triangle],
434                None,
435                None,
436            ),
437            MixedGeometry::<f64, CiarletElement<f64, IdentityMap, f64>>::new(
438                &[ReferenceCellType::Triangle, ReferenceCellType::Triangle],
439                points,
440                &[0, 1, 2, 2, 1, 3],
441                &[family],
442                &[0, 0],
443            ),
444        )
445    }
446
447    fn example_grid_mixed() -> MixedGrid<f64, CiarletElement<f64, IdentityMap, f64>> {
448        let mut points = rlst_dynamic_array!(f64, [2, 8]);
449        *points.get_mut([0, 0]).unwrap() = 0.0;
450        *points.get_mut([1, 0]).unwrap() = 0.0;
451        *points.get_mut([0, 1]).unwrap() = 1.0;
452        *points.get_mut([1, 1]).unwrap() = 0.0;
453        *points.get_mut([0, 2]).unwrap() = 2.0;
454        *points.get_mut([1, 2]).unwrap() = 0.0;
455        *points.get_mut([0, 3]).unwrap() = 4.0;
456        *points.get_mut([1, 3]).unwrap() = 0.0;
457        *points.get_mut([0, 4]).unwrap() = 0.0;
458        *points.get_mut([1, 4]).unwrap() = 1.0;
459        *points.get_mut([0, 5]).unwrap() = 1.0;
460        *points.get_mut([1, 5]).unwrap() = 1.0;
461        *points.get_mut([0, 6]).unwrap() = 2.0;
462        *points.get_mut([1, 6]).unwrap() = 1.0;
463        *points.get_mut([0, 7]).unwrap() = 4.0;
464        *points.get_mut([1, 7]).unwrap() = 1.0;
465        let family = LagrangeElementFamily::<f64>::new(1, Continuity::Standard);
466        MixedGrid::new(
467            MixedTopology::new(
468                &[0, 5, 4, 0, 1, 5, 1, 2, 5, 6, 2, 3, 6, 7],
469                &[
470                    ReferenceCellType::Triangle,
471                    ReferenceCellType::Triangle,
472                    ReferenceCellType::Quadrilateral,
473                    ReferenceCellType::Quadrilateral,
474                ],
475                None,
476                None,
477            ),
478            MixedGeometry::<f64, CiarletElement<f64, IdentityMap, f64>>::new(
479                &[
480                    ReferenceCellType::Triangle,
481                    ReferenceCellType::Triangle,
482                    ReferenceCellType::Quadrilateral,
483                    ReferenceCellType::Quadrilateral,
484                ],
485                points,
486                &[0, 5, 4, 0, 1, 5, 1, 2, 5, 6, 2, 3, 6, 7],
487                &[family],
488                &[0, 0, 0, 0],
489            ),
490        )
491    }
492
493    #[test]
494    fn test_edges_triangle() {
495        let grid = example_grid_triangle();
496        let conn = reference_cell::connectivity(ReferenceCellType::Triangle);
497        for edge in grid.entity_iter(ReferenceCellType::Interval) {
498            let cell = grid
499                .entity(ReferenceCellType::Triangle, edge.cell_index)
500                .unwrap();
501            for (i, v) in izip!(
502                &conn[1][edge.entity_index][0],
503                edge.topology().sub_entity_iter(ReferenceCellType::Point)
504            ) {
505                assert_eq!(v, cell.topology().sub_entity(ReferenceCellType::Point, *i));
506            }
507        }
508    }
509
510    #[test]
511    fn test_geometry_map() {
512        let grid = example_grid_mixed();
513        let mut mapped_pts = rlst_dynamic_array!(f64, [2, 2]);
514
515        let mut pts = rlst_dynamic_array!(f64, [2, 2]);
516        pts[[0, 0]] = 0.0;
517        pts[[1, 0]] = 0.0;
518        pts[[0, 1]] = 0.5;
519        pts[[1, 1]] = 0.5;
520
521        let map = grid.geometry_map(ReferenceCellType::Triangle, 1, &pts);
522
523        map.physical_points(0, &mut mapped_pts);
524        assert_relative_eq!(mapped_pts[[0, 0]], 0.0, epsilon = 1e-10);
525        assert_relative_eq!(mapped_pts[[1, 0]], 0.0, epsilon = 1e-10);
526        assert_relative_eq!(mapped_pts[[0, 1]], 0.5, epsilon = 1e-10);
527        assert_relative_eq!(mapped_pts[[1, 1]], 1.0, epsilon = 1e-10);
528
529        map.physical_points(1, &mut mapped_pts);
530        assert_relative_eq!(mapped_pts[[0, 0]], 0.0, epsilon = 1e-10);
531        assert_relative_eq!(mapped_pts[[1, 0]], 0.0, epsilon = 1e-10);
532        assert_relative_eq!(mapped_pts[[0, 1]], 1.0, epsilon = 1e-10);
533        assert_relative_eq!(mapped_pts[[1, 1]], 0.5, epsilon = 1e-10);
534
535        let mut pts = rlst_dynamic_array!(f64, [2, 2]);
536        pts[[0, 0]] = 0.5;
537        pts[[1, 0]] = 0.0;
538        pts[[0, 1]] = 1.0;
539        pts[[1, 1]] = 1.0;
540        let map = grid.geometry_map(ReferenceCellType::Quadrilateral, 1, &pts);
541
542        map.physical_points(0, &mut mapped_pts);
543        assert_relative_eq!(mapped_pts[[0, 0]], 1.5, epsilon = 1e-10);
544        assert_relative_eq!(mapped_pts[[1, 0]], 0.0, epsilon = 1e-10);
545        assert_relative_eq!(mapped_pts[[0, 1]], 2.0, epsilon = 1e-10);
546        assert_relative_eq!(mapped_pts[[1, 1]], 1.0, epsilon = 1e-10);
547
548        map.physical_points(1, &mut mapped_pts);
549        assert_relative_eq!(mapped_pts[[0, 0]], 3.0, epsilon = 1e-10);
550        assert_relative_eq!(mapped_pts[[1, 0]], 0.0, epsilon = 1e-10);
551        assert_relative_eq!(mapped_pts[[0, 1]], 4.0, epsilon = 1e-10);
552        assert_relative_eq!(mapped_pts[[1, 1]], 1.0, epsilon = 1e-10);
553    }
554}
ndgrid/grid/local_grid/mixed/grid.rs

ndgrid/grid/local_grid/mixed/
grid.rs
