ndmesh/mesh/local_mesh/mixed/
mesh.rs

1//! Mixed mesh
2#[cfg(feature = "mpi")]
3use crate::ParallelMeshImpl;
4#[cfg(feature = "mpi")]
5use crate::{
6    MixedMeshBuilder,
7    traits::{Builder, DistributableMesh, 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, Mesh},
19    types::{Ownership, Scalar},
20};
21use itertools::izip;
22#[cfg(feature = "mpi")]
23use mpi::traits::{Communicator, Equivalence};
24use ndelement::{
25    ciarlet::{CiarletElement, LagrangeElementFamily, LagrangeVariant},
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 mesh entity
39#[derive(Debug)]
40pub struct MixedMeshEntity<
41    'a,
42    T: Scalar,
43    E: MappedFiniteElement<CellType = ReferenceCellType, T = T>,
44> {
45    mesh: &'a MixedMesh<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    MixedMeshEntity<'e, T, E>
56{
57    /// Create new
58    pub fn new(
59        mesh: &'e MixedMesh<T, E>,
60        cell_type: ReferenceCellType,
61        cell_index: usize,
62        entity_type: ReferenceCellType,
63        entity_index: usize,
64    ) -> Self {
65        Self {
66            mesh,
67            cell_type,
68            cell_index,
69            entity_type,
70            entity_index,
71            geometry_element_index: mesh.geometry.insertion_indices_to_element_indices
72                [mesh.topology.insertion_indices[&cell_type][cell_index]],
73            geometry_cell_index: mesh.geometry.insertion_indices_to_cell_indices
74                [mesh.topology.insertion_indices[&cell_type][cell_index]],
75        }
76    }
77}
78impl<T: Scalar, E: MappedFiniteElement<CellType = ReferenceCellType, T = T>> Entity
79    for MixedMeshEntity<'_, 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.mesh.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.mesh.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.mesh.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.mesh
122            .topology
123            .entity_id(self.entity_type, self.local_index())
124    }
125}
126
127/// Mixed mesh entity iterator
128#[derive(Debug)]
129pub struct MixedMeshEntityIter<
130    'a,
131    T: Scalar,
132    E: MappedFiniteElement<CellType = ReferenceCellType, T = T>,
133> {
134    mesh: &'a MixedMesh<T, E>,
135    entity_type: ReferenceCellType,
136    index: usize,
137}
138
139impl<'a, T: Scalar, E: MappedFiniteElement<CellType = ReferenceCellType, T = T>>
140    MixedMeshEntityIter<'a, T, E>
141{
142    /// Create new
143    pub fn new(mesh: &'a MixedMesh<T, E>, entity_type: ReferenceCellType) -> Self {
144        Self {
145            mesh,
146            entity_type,
147            index: 0,
148        }
149    }
150}
151impl<'a, T: Scalar, E: MappedFiniteElement<CellType = ReferenceCellType, T = T>> Iterator
152    for MixedMeshEntityIter<'a, T, E>
153{
154    type Item = MixedMeshEntity<'a, T, E>;
155
156    fn next(&mut self) -> Option<MixedMeshEntity<'a, T, E>> {
157        self.index += 1;
158        self.mesh.entity(self.entity_type, self.index - 1)
159    }
160}
161
162/// Serial mixed element mesh
163#[derive(Debug)]
164pub struct MixedMesh<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 SerializableMesh<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 MixedMesh<T, CiarletElement<T, IdentityMap, T>>
183{
184    type SerializableType = SerializableMesh<T>;
185    fn to_serializable(&self) -> SerializableMesh<T> {
186        SerializableMesh {
187            topology: self.topology.to_serializable(),
188            geometry: self.geometry.to_serializable(),
189        }
190    }
191    fn from_serializable(s: SerializableMesh<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>> MixedMesh<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> MixedMesh<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.push(LagrangeElementFamily::<T, T>::new(
228                        *d,
229                        Continuity::Standard,
230                        LagrangeVariant::Equispaced,
231                    ));
232                    index
233                })
234            })
235            .collect::<Vec<_>>();
236
237        let geometry = MixedGeometry::<T, CiarletElement<T, IdentityMap, T>>::new(
238            cell_types,
239            points,
240            cells,
241            &element_families,
242            &cell_families,
243        );
244
245        let mut start = 0;
246        let mut tcells = vec![];
247        for (t, f) in izip!(cell_types, &cell_families) {
248            let tpoints_per_cell = reference_cell::entity_counts(*t)[0];
249            for i in 0..tpoints_per_cell {
250                tcells.push(cells[start + i]);
251            }
252            start += element_families[*f].element(*t).dim();
253        }
254
255        let topology = MixedTopology::new(&tcells, cell_types, None, None);
256
257        Self { topology, geometry }
258    }
259}
260
261impl<T: Scalar, E: MappedFiniteElement<CellType = ReferenceCellType, T = T>> Mesh
262    for MixedMesh<T, E>
263{
264    type T = T;
265    type Entity<'a>
266        = MixedMeshEntity<'a, T, E>
267    where
268        Self: 'a;
269    type GeometryMap<'a>
270        = GeometryMap<'a, T, BaseArray<VectorContainer<T>, 2>, BaseArray<VectorContainer<usize>, 2>>
271    where
272        Self: 'a;
273    type EntityDescriptor = ReferenceCellType;
274    type EntityIter<'a>
275        = MixedMeshEntityIter<'a, T, E>
276    where
277        Self: 'a;
278
279    fn geometry_dim(&self) -> usize {
280        self.geometry.dim()
281    }
282    fn topology_dim(&self) -> usize {
283        self.topology.dim()
284    }
285
286    fn entity(
287        &self,
288        entity_type: ReferenceCellType,
289        local_index: usize,
290    ) -> Option<Self::Entity<'_>> {
291        let dim = reference_cell::dim(entity_type);
292        if local_index < self.topology.entity_count(entity_type) {
293            if dim == self.topology_dim() {
294                Some(MixedMeshEntity::new(
295                    self,
296                    entity_type,
297                    local_index,
298                    entity_type,
299                    0,
300                ))
301            } else {
302                for t in &self.topology.entity_types()[self.topology_dim()] {
303                    if let Some(cell) =
304                        self.topology.upward_connectivity[&entity_type][t][local_index].first()
305                    {
306                        let dc = &self.topology.downward_connectivity[t][&entity_type];
307                        if let Some(index) =
308                            (0..dc.shape()[0]).position(|i| dc[[i, *cell]] == local_index)
309                        {
310                            return Some(MixedMeshEntity::new(self, *t, *cell, entity_type, index));
311                        }
312                    }
313                }
314                None
315            }
316        } else {
317            None
318        }
319    }
320
321    fn entity_types(&self, dim: usize) -> &[ReferenceCellType] {
322        &self.topology.entity_types()[dim]
323    }
324
325    fn entity_count(&self, entity_type: ReferenceCellType) -> usize {
326        self.topology.entity_count(entity_type)
327    }
328
329    fn entity_iter(&self, entity_type: ReferenceCellType) -> Self::EntityIter<'_> {
330        MixedMeshEntityIter::new(self, entity_type)
331    }
332
333    fn entity_from_id(
334        &self,
335        entity_type: ReferenceCellType,
336        id: usize,
337    ) -> Option<Self::Entity<'_>> {
338        self.topology.ids_to_indices[&entity_type]
339            .get(&id)
340            .map(|i| self.entity(entity_type, *i))?
341    }
342
343    fn geometry_map<Array2Impl: ValueArrayImpl<T, 2>>(
344        &self,
345        entity_type: ReferenceCellType,
346        geometry_degree: usize,
347        points: &Array<Array2Impl, 2>,
348    ) -> GeometryMap<'_, T, BaseArray<VectorContainer<T>, 2>, BaseArray<VectorContainer<usize>, 2>>
349    {
350        debug_assert!(points.shape()[0] == reference_cell::dim(entity_type));
351
352        for i in 0..self.geometry.element_count() {
353            let e = self.geometry.element(i);
354            if e.cell_type() == entity_type && e.lagrange_superdegree() == geometry_degree {
355                return GeometryMap::new(e, points, self.geometry.points(), self.geometry.cells(i));
356            }
357        }
358        unimplemented!();
359    }
360}
361
362#[cfg(feature = "mpi")]
363impl<T: Scalar + Equivalence, E: MappedFiniteElement<CellType = ReferenceCellType, T = T>>
364    DistributableMesh for MixedMesh<T, E>
365{
366    type ParallelMesh<'a, C: Communicator + 'a> =
367        ParallelMeshImpl<'a, C, MixedMesh<T, CiarletElement<T, IdentityMap, T>>>;
368
369    fn distribute<'a, C: Communicator>(
370        &self,
371        comm: &'a C,
372        partitioner: GraphPartitioner,
373    ) -> Self::ParallelMesh<'a, C> {
374        let mut b = MixedMeshBuilder::<T>::new(self.geometry.dim());
375        let pts = self.geometry.points();
376        for p in 0..pts.shape()[1] {
377            b.add_point(
378                p,
379                &pts.data().unwrap()[p * pts.shape()[0]..(p + 1) * pts.shape()[0]],
380            );
381        }
382
383        for (c, (element_i, cell_i)) in izip!(
384            &self.geometry.insertion_indices_to_element_indices,
385            &self.geometry.insertion_indices_to_cell_indices
386        )
387        .enumerate()
388        {
389            let e = self.geometry.element(*element_i);
390            let cells = self.geometry.cells(*element_i);
391            b.add_cell(
392                c,
393                (
394                    e.cell_type(),
395                    e.lagrange_superdegree(),
396                    &cells.data().unwrap()
397                        [cell_i * cells.shape()[0]..(cell_i + 1) * cells.shape()[0]],
398                ),
399            );
400        }
401        b.create_parallel_mesh_root(comm, partitioner)
402    }
403}
404
405#[cfg(test)]
406mod test {
407    use super::*;
408    use crate::traits::{GeometryMap, Topology};
409    use approx::*;
410    use itertools::izip;
411    use rlst::rlst_dynamic_array;
412
413    fn example_mesh_triangle() -> MixedMesh<f64, CiarletElement<f64, IdentityMap, f64>> {
414        let mut points = rlst_dynamic_array!(f64, [3, 4]);
415        *points.get_mut([0, 0]).unwrap() = 0.0;
416        *points.get_mut([1, 0]).unwrap() = 0.0;
417        *points.get_mut([2, 0]).unwrap() = 1.0;
418        *points.get_mut([0, 1]).unwrap() = 1.0;
419        *points.get_mut([1, 1]).unwrap() = 0.0;
420        *points.get_mut([2, 1]).unwrap() = 0.0;
421        *points.get_mut([0, 2]).unwrap() = 0.0;
422        *points.get_mut([1, 2]).unwrap() = 1.0;
423        *points.get_mut([2, 2]).unwrap() = 0.0;
424        *points.get_mut([0, 3]).unwrap() = 2.0;
425        *points.get_mut([1, 3]).unwrap() = 1.0;
426        *points.get_mut([2, 3]).unwrap() = 0.0;
427        let family =
428            LagrangeElementFamily::<f64>::new(1, Continuity::Standard, LagrangeVariant::Equispaced);
429        MixedMesh::new(
430            MixedTopology::new(
431                &[0, 1, 2, 2, 1, 3],
432                &[ReferenceCellType::Triangle, ReferenceCellType::Triangle],
433                None,
434                None,
435            ),
436            MixedGeometry::<f64, CiarletElement<f64, IdentityMap, f64>>::new(
437                &[ReferenceCellType::Triangle, ReferenceCellType::Triangle],
438                points,
439                &[0, 1, 2, 2, 1, 3],
440                &[family],
441                &[0, 0],
442            ),
443        )
444    }
445
446    fn example_mesh_mixed() -> MixedMesh<f64, CiarletElement<f64, IdentityMap, f64>> {
447        let mut points = rlst_dynamic_array!(f64, [2, 8]);
448        *points.get_mut([0, 0]).unwrap() = 0.0;
449        *points.get_mut([1, 0]).unwrap() = 0.0;
450        *points.get_mut([0, 1]).unwrap() = 1.0;
451        *points.get_mut([1, 1]).unwrap() = 0.0;
452        *points.get_mut([0, 2]).unwrap() = 2.0;
453        *points.get_mut([1, 2]).unwrap() = 0.0;
454        *points.get_mut([0, 3]).unwrap() = 4.0;
455        *points.get_mut([1, 3]).unwrap() = 0.0;
456        *points.get_mut([0, 4]).unwrap() = 0.0;
457        *points.get_mut([1, 4]).unwrap() = 1.0;
458        *points.get_mut([0, 5]).unwrap() = 1.0;
459        *points.get_mut([1, 5]).unwrap() = 1.0;
460        *points.get_mut([0, 6]).unwrap() = 2.0;
461        *points.get_mut([1, 6]).unwrap() = 1.0;
462        *points.get_mut([0, 7]).unwrap() = 4.0;
463        *points.get_mut([1, 7]).unwrap() = 1.0;
464        let family =
465            LagrangeElementFamily::<f64>::new(1, Continuity::Standard, LagrangeVariant::Equispaced);
466        MixedMesh::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 mesh = example_mesh_triangle();
496        let conn = reference_cell::connectivity(ReferenceCellType::Triangle);
497        for edge in mesh.entity_iter(ReferenceCellType::Interval) {
498            let cell = mesh
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 mesh = example_mesh_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 = mesh.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 = mesh.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}
ndmesh/mesh/local_mesh/mixed/mesh.rs

ndmesh/mesh/local_mesh/mixed/
mesh.rs
