ndgrid/geometry/single_element/
geometry.rs

1//! Geometry where each entity of a given dimension is represented by the same element
2#[cfg(feature = "serde")]
3use crate::traits::ConvertToSerializable;
4use crate::types::Scalar;
5#[cfg(feature = "serde")]
6use ndelement::{
7    ciarlet::{CiarletElement, lagrange},
8    map::IdentityMap,
9    traits::FiniteElement,
10    types::Continuity,
11};
12use ndelement::{
13    reference_cell,
14    traits::{ElementFamily, MappedFiniteElement},
15    types::ReferenceCellType,
16};
17use rlst::{DynArray, rlst_dynamic_array};
18use std::fmt::{Debug, Formatter};
19
20/// Single element geometry
21pub struct SingleElementGeometry<T: Scalar, E: MappedFiniteElement> {
22    points: DynArray<T, 2>,
23    cells: DynArray<usize, 2>,
24    elements: Vec<E>,
25}
26
27impl<T: Scalar, E: MappedFiniteElement> Debug for SingleElementGeometry<T, E> {
28    fn fmt(&self, f: &mut Formatter<'_>) -> Result<(), std::fmt::Error> {
29        f.debug_struct("SingleElementGeometry")
30            .field("points", &self.points)
31            .field("cells", &self.cells)
32            .finish()
33    }
34}
35
36#[cfg(feature = "serde")]
37#[derive(serde::Serialize, Debug, serde::Deserialize)]
38#[serde(bound = "for<'de2> T: serde::Deserialize<'de2>")]
39pub struct SerializableGeometry<T: Scalar + serde::Serialize>
40where
41    for<'de2> T: serde::Deserialize<'de2>,
42{
43    points: (Vec<T>, [usize; 2]),
44    cells: (Vec<usize>, [usize; 2]),
45    elements: Vec<(ReferenceCellType, usize)>,
46}
47
48#[cfg(feature = "serde")]
49impl<T: Scalar + serde::Serialize> ConvertToSerializable
50    for SingleElementGeometry<T, CiarletElement<T, IdentityMap, T>>
51where
52    for<'de2> T: serde::Deserialize<'de2>,
53{
54    type SerializableType = SerializableGeometry<T>;
55    fn to_serializable(&self) -> SerializableGeometry<T> {
56        SerializableGeometry {
57            points: (self.points.data().unwrap().to_vec(), self.points.shape()),
58            cells: (self.cells.data().unwrap().to_vec(), self.cells.shape()),
59            elements: self
60                .elements
61                .iter()
62                .map(|e| (e.cell_type(), e.lagrange_superdegree()))
63                .collect::<Vec<_>>(),
64        }
65    }
66    fn from_serializable(s: SerializableGeometry<T>) -> Self {
67        Self {
68            points: {
69                let (data, shape) = &s.points;
70                let mut p = DynArray::<T, 2>::from_shape(*shape);
71                p.data_mut().unwrap().copy_from_slice(data.as_slice());
72                p
73            },
74            cells: {
75                let (data, shape) = &s.cells;
76                let mut c = DynArray::<usize, 2>::from_shape(*shape);
77                c.data_mut().unwrap().copy_from_slice(data.as_slice());
78                c
79            },
80            elements: s
81                .elements
82                .iter()
83                .map(|(cell, degree)| lagrange::create(*cell, *degree, Continuity::Standard))
84                .collect::<Vec<_>>(),
85        }
86    }
87}
88
89impl<T: Scalar, E: MappedFiniteElement> SingleElementGeometry<T, E> {
90    /// Create single element geometry
91    pub fn new(
92        cell_type: ReferenceCellType,
93        points: DynArray<T, 2>,
94        cells_input: &[usize],
95        element_family: &impl ElementFamily<T = T, CellType = ReferenceCellType, FiniteElement = E>,
96    ) -> Self {
97        let mut elements = vec![];
98        for et in reference_cell::entity_types(cell_type)
99            .iter()
100            .skip(1)
101            .filter(|i| !i.is_empty())
102        {
103            for c in et.iter().skip(1) {
104                if *c != et[0] {
105                    panic!("Unsupported cell type for SingleElementGeometry: {cell_type:?}");
106                }
107            }
108            elements.push(element_family.element(et[0]));
109        }
110        let points_per_cell = elements[elements.len() - 1].dim();
111        let mut cells = rlst_dynamic_array!(
112            usize,
113            [points_per_cell, cells_input.len() / points_per_cell]
114        );
115        cells.data_mut().unwrap().copy_from_slice(cells_input);
116        Self {
117            points,
118            cells,
119            elements,
120        }
121    }
122    /// Geometric dimension
123    pub fn dim(&self) -> usize {
124        self.points().shape()[0]
125    }
126    /// Points
127    pub fn points(&self) -> &DynArray<T, 2> {
128        &self.points
129    }
130    /// Cells
131    pub fn cells(&self) -> &DynArray<usize, 2> {
132        &self.cells
133    }
134    /// Element for a sub-entity
135    pub fn entity_element(&self, tdim: usize) -> &E {
136        &self.elements[tdim - 1]
137    }
138    /// Element for a cell
139    pub fn element(&self) -> &E {
140        &self.elements[self.elements.len() - 1]
141    }
142}