ndgrid/geometry/mixed/
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;
5use itertools::izip;
6#[cfg(feature = "serde")]
7use ndelement::{
8    ciarlet::{CiarletElement, lagrange},
9    map::IdentityMap,
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::collections::HashMap;
19use std::fmt::{Debug, Formatter};
20
21/// Single element geometry
22pub struct MixedGeometry<T: Scalar, E: MappedFiniteElement> {
23    points: DynArray<T, 2>,
24    cells: Vec<DynArray<usize, 2>>,
25    elements: Vec<HashMap<ReferenceCellType, E>>,
26    pub(crate) insertion_indices_to_element_indices: Vec<usize>,
27    pub(crate) insertion_indices_to_cell_indices: Vec<usize>,
28}
29
30impl<T: Scalar, E: MappedFiniteElement> Debug for MixedGeometry<T, E> {
31    fn fmt(&self, f: &mut Formatter<'_>) -> Result<(), std::fmt::Error> {
32        f.debug_struct("MixedGeometry")
33            .field("points", &self.points)
34            .field("cells", &self.cells)
35            .finish()
36    }
37}
38
39#[cfg(feature = "serde")]
40#[derive(serde::Serialize, Debug, serde::Deserialize)]
41#[serde(bound = "for<'de2> T: serde::Deserialize<'de2>")]
42pub struct SerializableGeometry<T: Scalar + serde::Serialize>
43where
44    for<'de2> T: serde::Deserialize<'de2>,
45{
46    points: (Vec<T>, [usize; 2]),
47    cells: Vec<(Vec<usize>, [usize; 2])>,
48    elements: Vec<HashMap<ReferenceCellType, usize>>,
49    insertion_indices_to_element_indices: Vec<usize>,
50    insertion_indices_to_cell_indices: Vec<usize>,
51}
52
53#[cfg(feature = "serde")]
54impl<T: Scalar + serde::Serialize> ConvertToSerializable
55    for MixedGeometry<T, CiarletElement<T, IdentityMap, T>>
56where
57    for<'de2> T: serde::Deserialize<'de2>,
58{
59    type SerializableType = SerializableGeometry<T>;
60    fn to_serializable(&self) -> SerializableGeometry<T> {
61        SerializableGeometry {
62            points: (self.points.data().unwrap().to_vec(), self.points.shape()),
63            cells: self
64                .cells
65                .iter()
66                .map(|c| (c.data().unwrap().to_vec(), c.shape()))
67                .collect::<Vec<_>>(),
68            elements: self
69                .elements
70                .iter()
71                .map(|a| {
72                    a.iter()
73                        .map(|(b, c)| (*b, c.lagrange_superdegree()))
74                        .collect::<HashMap<_, _>>()
75                })
76                .collect::<Vec<_>>(),
77            insertion_indices_to_element_indices: self.insertion_indices_to_element_indices.clone(),
78            insertion_indices_to_cell_indices: self.insertion_indices_to_cell_indices.clone(),
79        }
80    }
81    fn from_serializable(s: SerializableGeometry<T>) -> Self {
82        Self {
83            points: {
84                let (data, shape) = &s.points;
85                let mut p = DynArray::<T, 2>::from_shape(*shape);
86                p.data_mut().unwrap().copy_from_slice(data.as_slice());
87                p
88            },
89            cells: s
90                .cells
91                .iter()
92                .map(|(data, shape)| {
93                    let mut c = DynArray::<usize, 2>::from_shape(*shape);
94                    c.data_mut().unwrap().copy_from_slice(data.as_slice());
95                    c
96                })
97                .collect::<Vec<_>>(),
98            elements: s
99                .elements
100                .iter()
101                .map(|a| {
102                    a.iter()
103                        .map(|(cell, degree)| {
104                            (
105                                *cell,
106                                lagrange::create(*cell, *degree, Continuity::Standard),
107                            )
108                        })
109                        .collect::<HashMap<_, _>>()
110                })
111                .collect::<Vec<_>>(),
112            insertion_indices_to_element_indices: s.insertion_indices_to_element_indices,
113            insertion_indices_to_cell_indices: s.insertion_indices_to_cell_indices,
114        }
115    }
116}
117
118impl<T: Scalar, E: MappedFiniteElement> MixedGeometry<T, E> {
119    /// Create single element geometry
120    pub fn new(
121        cell_types_in: &[ReferenceCellType],
122        points: DynArray<T, 2>,
123        cells_input: &[usize],
124        element_families: &[impl ElementFamily<
125            T = T,
126            CellType = ReferenceCellType,
127            FiniteElement = E,
128        >],
129        cell_families: &[usize],
130    ) -> Self {
131        let mut element_info = HashMap::new();
132        let mut elements = vec![];
133        let mut cells = vec![];
134        let mut cell_counts = vec![];
135        let mut points_per_cell = vec![];
136        let mut insertion_indices_to_element_indices = vec![];
137        let mut insertion_indices_to_cell_indices = vec![];
138
139        let mut start = 0;
140
141        for (findex, cell_type) in izip!(cell_families, cell_types_in) {
142            let eindex = *element_info.entry((findex, cell_type)).or_insert_with(|| {
143                let n = elements.len();
144
145                let mut new_e = HashMap::new();
146                for et in reference_cell::entity_types(*cell_type)
147                    .iter()
148                    .skip(1)
149                    .filter(|i| !i.is_empty())
150                {
151                    for e in et {
152                        new_e
153                            .entry(*e)
154                            .or_insert(element_families[*findex].element(*e));
155                    }
156                }
157                points_per_cell.push(new_e[cell_type].dim());
158                elements.push(new_e);
159                cells.push(vec![]);
160                cell_counts.push(0);
161                n
162            });
163
164            insertion_indices_to_element_indices.push(eindex);
165            insertion_indices_to_cell_indices.push(cell_counts[eindex]);
166            for i in 0..points_per_cell[eindex] {
167                cells[eindex].push(cells_input[start + i]);
168            }
169            cell_counts[eindex] += 1;
170            start += points_per_cell[eindex];
171        }
172        let cells = izip!(cell_counts, points_per_cell, cells)
173            .map(|(ncells, npts, c_in)| {
174                let mut c = rlst_dynamic_array!(usize, [npts, ncells]);
175                c.data_mut().unwrap().copy_from_slice(&c_in);
176                c
177            })
178            .collect::<Vec<_>>();
179        Self {
180            points,
181            cells,
182            elements,
183            insertion_indices_to_element_indices,
184            insertion_indices_to_cell_indices,
185        }
186    }
187    /// Geometric dimension
188    pub fn dim(&self) -> usize {
189        self.points().shape()[0]
190    }
191    /// Points
192    pub fn points(&self) -> &DynArray<T, 2> {
193        &self.points
194    }
195    /// Cells
196    pub fn cells(&self, element_index: usize) -> &DynArray<usize, 2> {
197        &self.cells[element_index]
198    }
199    /// Element for a cell
200    pub fn element(&self, element_index: usize) -> &E {
201        for (ct, e) in &self.elements[element_index] {
202            if reference_cell::dim(*ct) == self.dim() {
203                return e;
204            }
205        }
206        panic!("Could not find element");
207    }
208    /// Number of elments
209    pub fn element_count(&self) -> usize {
210        self.elements.len()
211    }
212}
ndgrid/geometry/mixed/geometry.rs

ndgrid/geometry/mixed/
geometry.rs
