ndgrid/shapes/
screen.rs

1//! Regular sphere grid
2
3#[cfg(feature = "mpi")]
4use crate::{ParallelGridImpl, traits::ParallelBuilder, types::GraphPartitioner};
5use crate::{
6    grid::local_grid::{SingleElementGrid, SingleElementGridBuilder},
7    traits::Builder,
8    types::Scalar,
9};
10#[cfg(feature = "mpi")]
11use mpi::traits::{Communicator, Equivalence};
12use ndelement::{ciarlet::CiarletElement, map::IdentityMap, types::ReferenceCellType};
13
14/// Add points and cells for a square screen to builder
15fn screen_add_points_and_cells<T: Scalar>(
16    b: &mut SingleElementGridBuilder<T>,
17    ncells: usize,
18    cell_type: ReferenceCellType,
19) {
20    let zero = T::from(0.0).unwrap();
21    let n = T::from(ncells).unwrap();
22    for y in 0..ncells + 1 {
23        for x in 0..ncells + 1 {
24            b.add_point(
25                y * (ncells + 1) + x,
26                &[T::from(x).unwrap() / n, T::from(y).unwrap() / n, zero],
27            );
28        }
29    }
30    match cell_type {
31        ReferenceCellType::Triangle => {
32            for y in 0..ncells {
33                for x in 0..ncells {
34                    b.add_cell(
35                        2 * y * ncells + 2 * x,
36                        &[
37                            y * (ncells + 1) + x,
38                            y * (ncells + 1) + x + 1,
39                            y * (ncells + 1) + x + ncells + 2,
40                        ],
41                    );
42                    b.add_cell(
43                        2 * y * ncells + 2 * x + 1,
44                        &[
45                            y * (ncells + 1) + x,
46                            y * (ncells + 1) + x + ncells + 2,
47                            y * (ncells + 1) + x + ncells + 1,
48                        ],
49                    );
50                }
51            }
52        }
53        ReferenceCellType::Quadrilateral => {
54            for y in 0..ncells {
55                for x in 0..ncells {
56                    b.add_cell(
57                        y * ncells + x,
58                        &[
59                            y * (ncells + 1) + x,
60                            y * (ncells + 1) + x + 1,
61                            y * (ncells + 1) + x + ncells + 1,
62                            y * (ncells + 1) + x + ncells + 2,
63                        ],
64                    );
65                }
66            }
67        }
68        _ => {
69            panic!("Unsupported cell type: {cell_type:?}");
70        }
71    }
72}
73
74/// Create a grid of a square screen
75///
76/// Create a grid of the square \[0,1\]^2. The input ncells is the number of cells
77/// along each side of the square.
78pub fn screen<T: Scalar>(
79    ncells: usize,
80    cell_type: ReferenceCellType,
81) -> SingleElementGrid<T, CiarletElement<T, IdentityMap, T>> {
82    let mut b = SingleElementGridBuilder::new_with_capacity(
83        3,
84        (ncells + 1) * (ncells + 1),
85        match cell_type {
86            ReferenceCellType::Quadrilateral => ncells * ncells,
87            ReferenceCellType::Triangle => 2 * ncells * ncells,
88            _ => {
89                panic!("Unsupported cell type: {cell_type:?}");
90            }
91        },
92        (cell_type, 1),
93    );
94    screen_add_points_and_cells(&mut b, ncells, cell_type);
95    b.create_grid()
96}
97
98/// Create a grid of a square screen distributed in parallel
99#[cfg(feature = "mpi")]
100pub fn screen_distributed<T: Scalar + Equivalence, C: Communicator>(
101    comm: &C,
102    partitioner: GraphPartitioner,
103    ncells: usize,
104    cell_type: ReferenceCellType,
105) -> ParallelGridImpl<'_, C, SingleElementGrid<T, CiarletElement<T, IdentityMap, T>>> {
106    let mut b = SingleElementGridBuilder::new(3, (cell_type, 1));
107    if comm.rank() == 0 {
108        screen_add_points_and_cells(&mut b, ncells, cell_type);
109        b.create_parallel_grid_root(comm, partitioner)
110    } else {
111        b.create_parallel_grid(comm, 0)
112    }
113}
114
115#[cfg(test)]
116mod test {
117    use super::*;
118    use crate::traits::{GeometryMap, Grid};
119    use approx::assert_relative_eq;
120    use rlst::rlst_dynamic_array;
121
122    #[test]
123    fn test_screen_triangles() {
124        let _g1 = screen::<f64>(1, ReferenceCellType::Triangle);
125        let _g2 = screen::<f64>(2, ReferenceCellType::Triangle);
126        let _g3 = screen::<f64>(3, ReferenceCellType::Triangle);
127    }
128    #[test]
129    fn test_screen_triangles_normals() {
130        for i in 1..5 {
131            let g = screen::<f64>(i, ReferenceCellType::Triangle);
132            let mut points = rlst_dynamic_array!(f64, [2, 1]);
133            points[[0, 0]] = 1.0 / 3.0;
134            points[[1, 0]] = 1.0 / 3.0;
135            let map = g.geometry_map(ReferenceCellType::Triangle, 1, &points);
136            let mut mapped_pt = rlst_dynamic_array!(f64, [3, 1]);
137            let mut j = rlst_dynamic_array!(f64, [3, 2, 1]);
138            let mut jinv = rlst_dynamic_array!(f64, [2, 3, 1]);
139            let mut jdet = vec![0.0];
140            let mut normal = rlst_dynamic_array!(f64, [3, 1]);
141            for i in 0..g.entity_count(ReferenceCellType::Triangle) {
142                map.physical_points(i, &mut mapped_pt);
143                map.jacobians_inverses_dets_normals(i, &mut j, &mut jinv, &mut jdet, &mut normal);
144                assert!(normal[[2, 0]] > 0.0);
145                assert_relative_eq!(normal[[2, 0]], 1.0);
146            }
147        }
148    }
149
150    #[test]
151    fn test_screen_quadrilaterals() {
152        let _g1 = screen::<f64>(1, ReferenceCellType::Quadrilateral);
153        let _g2 = screen::<f64>(2, ReferenceCellType::Quadrilateral);
154        let _g3 = screen::<f64>(3, ReferenceCellType::Quadrilateral);
155    }
156
157    #[test]
158    fn test_screen_quadrilaterals_normals() {
159        for i in 1..5 {
160            let g = screen::<f64>(i, ReferenceCellType::Quadrilateral);
161            let mut points = rlst_dynamic_array!(f64, [2, 1]);
162            points[[0, 0]] = 1.0 / 3.0;
163            points[[1, 0]] = 1.0 / 3.0;
164            let map = g.geometry_map(ReferenceCellType::Quadrilateral, 1, &points);
165            let mut mapped_pt = rlst_dynamic_array!(f64, [3, 1]);
166            let mut j = rlst_dynamic_array!(f64, [3, 2, 1]);
167            let mut jinv = rlst_dynamic_array!(f64, [2, 3, 1]);
168            let mut jdet = vec![0.0];
169            let mut normal = rlst_dynamic_array!(f64, [3, 1]);
170            for i in 0..g.entity_count(ReferenceCellType::Quadrilateral) {
171                map.physical_points(i, &mut mapped_pt);
172                map.jacobians_inverses_dets_normals(i, &mut j, &mut jinv, &mut jdet, &mut normal);
173                assert!(normal[[2, 0]] > 0.0);
174                assert_relative_eq!(normal[[2, 0]], 1.0);
175            }
176        }
177    }
178}
ndgrid/shapes/screen.rs

ndgrid/shapes/
screen.rs
