ndmesh/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(
107                                    *cell,
108                                    *degree,
109                                    Continuity::Standard,
110                                    lagrange::Variant::Equispaced,
111                                ),
112                            )
113                        })
114                        .collect::<HashMap<_, _>>()
115                })
116                .collect::<Vec<_>>(),
117            insertion_indices_to_element_indices: s.insertion_indices_to_element_indices,
118            insertion_indices_to_cell_indices: s.insertion_indices_to_cell_indices,
119        }
120    }
121}
122
123impl<T: Scalar, E: MappedFiniteElement> MixedGeometry<T, E> {
124    /// Create single element geometry
125    pub fn new(
126        cell_types_in: &[ReferenceCellType],
127        points: DynArray<T, 2>,
128        cells_input: &[usize],
129        element_families: &[impl ElementFamily<
130            T = T,
131            CellType = ReferenceCellType,
132            FiniteElement = E,
133        >],
134        cell_families: &[usize],
135    ) -> Self {
136        let mut element_info = HashMap::new();
137        let mut elements = vec![];
138        let mut cells = vec![];
139        let mut cell_counts = vec![];
140        let mut points_per_cell = vec![];
141        let mut insertion_indices_to_element_indices = vec![];
142        let mut insertion_indices_to_cell_indices = vec![];
143
144        let mut start = 0;
145
146        for (findex, cell_type) in izip!(cell_families, cell_types_in) {
147            let eindex = *element_info.entry((findex, cell_type)).or_insert_with(|| {
148                let n = elements.len();
149
150                let mut new_e = HashMap::new();
151                for et in reference_cell::entity_types(*cell_type)
152                    .iter()
153                    .skip(1)
154                    .filter(|i| !i.is_empty())
155                {
156                    for e in et {
157                        new_e
158                            .entry(*e)
159                            .or_insert(element_families[*findex].element(*e));
160                    }
161                }
162                points_per_cell.push(new_e[cell_type].dim());
163                elements.push(new_e);
164                cells.push(vec![]);
165                cell_counts.push(0);
166                n
167            });
168
169            insertion_indices_to_element_indices.push(eindex);
170            insertion_indices_to_cell_indices.push(cell_counts[eindex]);
171            for i in 0..points_per_cell[eindex] {
172                cells[eindex].push(cells_input[start + i]);
173            }
174            cell_counts[eindex] += 1;
175            start += points_per_cell[eindex];
176        }
177        let cells = izip!(cell_counts, points_per_cell, cells)
178            .map(|(ncells, npts, c_in)| {
179                let mut c = rlst_dynamic_array!(usize, [npts, ncells]);
180                c.data_mut().unwrap().copy_from_slice(&c_in);
181                c
182            })
183            .collect::<Vec<_>>();
184        Self {
185            points,
186            cells,
187            elements,
188            insertion_indices_to_element_indices,
189            insertion_indices_to_cell_indices,
190        }
191    }
192    /// Geometric dimension
193    pub fn dim(&self) -> usize {
194        self.points().shape()[0]
195    }
196    /// Points
197    pub fn points(&self) -> &DynArray<T, 2> {
198        &self.points
199    }
200    /// Cells
201    pub fn cells(&self, element_index: usize) -> &DynArray<usize, 2> {
202        &self.cells[element_index]
203    }
204    /// Element for a cell
205    pub fn element(&self, element_index: usize) -> &E {
206        for (ct, e) in &self.elements[element_index] {
207            if reference_cell::dim(*ct) == self.dim() {
208                return e;
209            }
210        }
211        panic!("Could not find element");
212    }
213    /// Number of elments
214    pub fn element_count(&self) -> usize {
215        self.elements.len()
216    }
217}
ndmesh/geometry/mixed/geometry.rs

ndmesh/geometry/mixed/
geometry.rs
