ndgrid/shapes/
regular_sphere.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};
13use std::collections::{HashMap, hash_map::Entry::Vacant};
14
15/// Add points and cells for regular sphere to builder
16fn regular_sphere_triangle_add_points_and_cells<T: Scalar>(
17    b: &mut SingleElementGridBuilder<T>,
18    refinement_level: u32,
19) {
20    b.add_point(0, &[T::zero(), T::zero(), T::one()]);
21    b.add_point(1, &[T::one(), T::zero(), T::zero()]);
22    b.add_point(2, &[T::zero(), T::one(), T::zero()]);
23    b.add_point(3, &[-T::one(), T::zero(), T::zero()]);
24    b.add_point(4, &[T::zero(), -T::one(), T::zero()]);
25    b.add_point(5, &[T::zero(), T::zero(), -T::one()]);
26    let mut point_n = 6;
27
28    let mut cells = vec![
29        [0, 1, 2],
30        [0, 2, 3],
31        [0, 3, 4],
32        [0, 4, 1],
33        [5, 2, 1],
34        [5, 3, 2],
35        [5, 4, 3],
36        [5, 1, 4],
37    ];
38    let mut v = [
39        [T::zero(), T::zero(), T::zero()],
40        [T::zero(), T::zero(), T::zero()],
41        [T::zero(), T::zero(), T::zero()],
42    ];
43
44    for _level in 0..refinement_level {
45        let mut edge_points = HashMap::new();
46        let mut new_cells = Vec::with_capacity(4 * cells.len());
47        for c in &cells {
48            for (i, v_i) in v.iter_mut().enumerate() {
49                for (j, v_ij) in v_i.iter_mut().enumerate() {
50                    *v_ij = b.points[3 * c[i] + j];
51                }
52            }
53            let edges = [[1, 2], [0, 2], [0, 1]]
54                .iter()
55                .map(|[i, j]| {
56                    let mut pt_i = c[*i];
57                    let mut pt_j = c[*j];
58                    if pt_i > pt_j {
59                        std::mem::swap(&mut pt_i, &mut pt_j);
60                    }
61                    if let Vacant(e) = edge_points.entry((pt_i, pt_j)) {
62                        let v_i = v[*i];
63                        let v_j = v[*j];
64                        let mut new_pt = [
65                            T::from(0.5).unwrap() * (v_i[0] + v_j[0]),
66                            T::from(0.5).unwrap() * (v_i[1] + v_j[1]),
67                            T::from(0.5).unwrap() * (v_i[2] + v_j[2]),
68                        ];
69                        let size = (new_pt.iter().map(|x| x.powi(2)).sum::<T>()).sqrt();
70                        for i in new_pt.iter_mut() {
71                            *i /= size;
72                        }
73                        b.add_point(point_n, &new_pt);
74                        e.insert(point_n);
75                        point_n += 1;
76                    }
77                    edge_points[&(pt_i, pt_j)]
78                })
79                .collect::<Vec<_>>();
80            new_cells.push([c[0], edges[2], edges[1]]);
81            new_cells.push([c[1], edges[0], edges[2]]);
82            new_cells.push([c[2], edges[1], edges[0]]);
83            new_cells.push([edges[0], edges[1], edges[2]]);
84        }
85        cells = new_cells;
86    }
87    for (i, v) in cells.iter().enumerate() {
88        b.add_cell(i, v);
89    }
90}
91
92/// Add points and cells for regular sphere to builder
93fn regular_sphere_quadrilateral_add_points_and_cells<T: Scalar>(
94    b: &mut SingleElementGridBuilder<T>,
95    refinement_level: u32,
96) {
97    let k = T::from(1.0 / 3.0).unwrap().sqrt();
98    b.add_point(0, &[-k, -k, -k]);
99    b.add_point(1, &[-k, -k, k]);
100    b.add_point(2, &[-k, k, -k]);
101    b.add_point(3, &[-k, k, k]);
102    b.add_point(4, &[k, -k, -k]);
103    b.add_point(5, &[k, -k, k]);
104    b.add_point(6, &[k, k, -k]);
105    b.add_point(7, &[k, k, k]);
106
107    let mut point_n = 8;
108
109    let mut cells = vec![
110        [0, 1, 2, 3],
111        [4, 6, 5, 7],
112        [0, 4, 1, 5],
113        [2, 3, 6, 7],
114        [0, 2, 4, 6],
115        [1, 5, 3, 7],
116    ];
117    let mut v = [
118        [T::zero(), T::zero(), T::zero()],
119        [T::zero(), T::zero(), T::zero()],
120        [T::zero(), T::zero(), T::zero()],
121        [T::zero(), T::zero(), T::zero()],
122    ];
123
124    for _level in 0..refinement_level {
125        let mut edge_points = HashMap::new();
126        let mut new_cells = Vec::with_capacity(4 * cells.len());
127        for c in &cells {
128            for (i, v_i) in v.iter_mut().enumerate() {
129                for (j, v_ij) in v_i.iter_mut().enumerate() {
130                    *v_ij = b.points[3 * c[i] + j];
131                }
132            }
133            let edges = [[0, 1], [0, 2], [1, 3], [2, 3]]
134                .iter()
135                .map(|[i, j]| {
136                    let mut pt_i = c[*i];
137                    let mut pt_j = c[*j];
138                    if pt_i > pt_j {
139                        std::mem::swap(&mut pt_i, &mut pt_j);
140                    }
141                    if let Vacant(e) = edge_points.entry((pt_i, pt_j)) {
142                        let v_i = v[*i];
143                        let v_j = v[*j];
144                        let mut new_pt = [
145                            T::from(0.5).unwrap() * (v_i[0] + v_j[0]),
146                            T::from(0.5).unwrap() * (v_i[1] + v_j[1]),
147                            T::from(0.5).unwrap() * (v_i[2] + v_j[2]),
148                        ];
149                        let size = (new_pt.iter().map(|x| x.powi(2)).sum::<T>()).sqrt();
150                        for i in new_pt.iter_mut() {
151                            *i /= size;
152                        }
153                        b.add_point(point_n, &new_pt);
154                        e.insert(point_n);
155                        point_n += 1;
156                    }
157                    edge_points[&(pt_i, pt_j)]
158                })
159                .collect::<Vec<_>>();
160
161            new_cells.push([c[0], edges[0], edges[1], point_n]);
162            new_cells.push([edges[0], c[1], point_n, edges[2]]);
163            new_cells.push([edges[1], point_n, c[2], edges[3]]);
164            new_cells.push([point_n, edges[2], edges[3], c[3]]);
165            let mut new_pt = [
166                T::from(0.25).unwrap() * v.iter().map(|i| i[0]).sum(),
167                T::from(0.25).unwrap() * v.iter().map(|i| i[1]).sum(),
168                T::from(0.25).unwrap() * v.iter().map(|i| i[2]).sum(),
169            ];
170            let size = (new_pt.iter().map(|x| x.powi(2)).sum::<T>()).sqrt();
171            for i in new_pt.iter_mut() {
172                *i /= size;
173            }
174            b.add_point(point_n, &new_pt);
175            point_n += 1;
176        }
177        cells = new_cells;
178    }
179    for (i, v) in cells.iter().enumerate() {
180        b.add_cell(i, v);
181    }
182}
183/// Create a surface grid of a regular sphere
184///
185/// A regular sphere is created by starting with a regular octahedron. The shape is then refined `refinement_level` times.
186/// Each time the grid is refined, each triangle is split into four triangles (by adding lines connecting the midpoints of
187/// each edge). The new points are then scaled so that they are a distance of 1 from the origin.
188pub fn regular_sphere<T: Scalar>(
189    refinement_level: u32,
190    cell_type: ReferenceCellType,
191) -> SingleElementGrid<T, CiarletElement<T, IdentityMap, T>> {
192    let mut b = SingleElementGridBuilder::new_with_capacity(
193        3,
194        match cell_type {
195            ReferenceCellType::Triangle => 2 + usize::pow(4, refinement_level + 1),
196            ReferenceCellType::Quadrilateral => 2 + 6 * usize::pow(4, refinement_level),
197            _ => {
198                panic!("Unsupported cell type: {cell_type:?}");
199            }
200        },
201        match cell_type {
202            ReferenceCellType::Triangle => 2 * usize::pow(4, refinement_level + 1),
203            ReferenceCellType::Quadrilateral => 6 * usize::pow(4, refinement_level),
204            _ => {
205                panic!("Unsupported cell type: {cell_type:?}");
206            }
207        },
208        (cell_type, 1),
209    );
210    match cell_type {
211        ReferenceCellType::Triangle => {
212            regular_sphere_triangle_add_points_and_cells(&mut b, refinement_level)
213        }
214        ReferenceCellType::Quadrilateral => {
215            regular_sphere_quadrilateral_add_points_and_cells(&mut b, refinement_level)
216        }
217        _ => {
218            panic!("Unsupported cell type: {cell_type:?}");
219        }
220    }
221    b.create_grid()
222}
223
224/// Create a grid of a regular sphere distributed in parallel
225#[cfg(feature = "mpi")]
226pub fn regular_sphere_distributed<T: Scalar + Equivalence, C: Communicator>(
227    comm: &C,
228    partitioner: GraphPartitioner,
229    refinement_level: u32,
230    cell_type: ReferenceCellType,
231) -> ParallelGridImpl<'_, C, SingleElementGrid<T, CiarletElement<T, IdentityMap, T>>> {
232    let mut b = SingleElementGridBuilder::new(3, (cell_type, 1));
233    if comm.rank() == 0 {
234        match cell_type {
235            ReferenceCellType::Triangle => {
236                regular_sphere_triangle_add_points_and_cells(&mut b, refinement_level)
237            }
238            ReferenceCellType::Quadrilateral => {
239                regular_sphere_quadrilateral_add_points_and_cells(&mut b, refinement_level)
240            }
241            _ => {
242                panic!("Unsupported cell type: {cell_type:?}");
243            }
244        }
245        b.create_parallel_grid_root(comm, partitioner)
246    } else {
247        b.create_parallel_grid(comm, 0)
248    }
249}
250
251#[cfg(test)]
252mod test {
253    use super::*;
254    use crate::traits::{GeometryMap, Grid};
255    use approx::assert_relative_eq;
256    use paste::paste;
257    use rlst::rlst_dynamic_array;
258
259    macro_rules! test_regular_sphere {
260        ($cell:ident, $order:expr) => {
261            paste! {
262                #[test]
263                fn [<test_regular_sphere_ $cell:lower _ $order>]() {
264                    let _g = regular_sphere::<f64>([<$order>], ReferenceCellType::[<$cell>]);
265                }
266            }
267        };
268    }
269
270    test_regular_sphere!(Triangle, 0);
271    test_regular_sphere!(Triangle, 1);
272    test_regular_sphere!(Triangle, 2);
273    test_regular_sphere!(Triangle, 3);
274    test_regular_sphere!(Quadrilateral, 0);
275    test_regular_sphere!(Quadrilateral, 1);
276    test_regular_sphere!(Quadrilateral, 2);
277    test_regular_sphere!(Quadrilateral, 3);
278
279    macro_rules! test_normals {
280        ($cell:ident, $order:expr) => {
281            paste! {
282                #[test]
283                fn [<test_normal_is_outward_ $cell:lower _ $order>]() {
284                    let g = regular_sphere::<f64>([<$order>], ReferenceCellType::[<$cell>]);
285                    let mut points = rlst_dynamic_array!(f64, [2, 1]);
286                    points[[0, 0]] = 1.0 / 3.0;
287                    points[[1, 0]] = 1.0 / 3.0;
288                    let map = g.geometry_map(ReferenceCellType::[<$cell>], 1, &points);
289                    let mut mapped_pt = rlst_dynamic_array!(f64, [3, 1]);
290                    let mut j = rlst_dynamic_array!(f64, [3, 2, 1]);
291                    let mut jinv = rlst_dynamic_array!(f64, [2, 3, 1]);
292                    let mut jdet = vec![0.0];
293                    let mut normal = rlst_dynamic_array!(f64, [3, 1]);
294                    for i in 0..g.entity_count(ReferenceCellType::[<$cell>]) {
295                        map.physical_points(i, &mut mapped_pt);
296                        map.jacobians_inverses_dets_normals(i, &mut j, &mut jinv, &mut jdet, &mut normal);
297                        let dot = mapped_pt
298                            .iter_value()
299                            .zip(normal.iter_value())
300                            .map(|(i, j)| i * j)
301                            .sum::<f64>();
302                        assert!(dot > 0.0);
303                    }
304                }
305
306                #[test]
307                fn [<test_normal_is_unit_ $cell:lower _ $order>]() {
308                    let g = regular_sphere::<f64>([<$order>], ReferenceCellType::[<$cell>]);
309                    let mut points = rlst_dynamic_array!(f64, [2, 1]);
310                    points[[0, 0]] = 1.0 / 3.0;
311                    points[[1, 0]] = 1.0 / 3.0;
312                    let map = g.geometry_map(ReferenceCellType::[<$cell>], 1, &points);
313                    let mut mapped_pt = rlst_dynamic_array!(f64, [3, 1]);
314                    let mut j = rlst_dynamic_array!(f64, [3, 2, 1]);
315                    let mut jinv = rlst_dynamic_array!(f64, [2, 3, 1]);
316                    let mut jdet = vec![0.0];
317                    let mut normal = rlst_dynamic_array!(f64, [3, 1]);
318                    for i in 0..g.entity_count(ReferenceCellType::[<$cell>]) {
319                        map.physical_points(i, &mut mapped_pt);
320                        map.jacobians_inverses_dets_normals(i, &mut j, &mut jinv, &mut jdet, &mut normal);
321                        let dot = normal
322                            .iter_value()
323                            .zip(normal.iter_value())
324                            .map(|(i, j)| i * j)
325                            .sum::<f64>();
326                        assert_relative_eq!(dot, 1.0, epsilon = 1e-10);
327                    }
328                }
329            }
330        };
331    }
332
333    test_normals!(Triangle, 0);
334    test_normals!(Triangle, 1);
335    test_normals!(Triangle, 2);
336    test_normals!(Quadrilateral, 0);
337    test_normals!(Quadrilateral, 1);
338    test_normals!(Quadrilateral, 2);
339}
ndgrid/shapes/regular_sphere.rs

ndgrid/shapes/
regular_sphere.rs
