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

ndgrid/grid/local_grid/mixed/
grid.rs
