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_add_points_and_cells<T: Scalar>(
17    b: &mut SingleElementGridBuilder<T>,
18    refinement_level: u32,
19) {
20    let zero = T::from(0.0).unwrap();
21    let one = T::from(1.0).unwrap();
22    let half = T::from(0.5).unwrap();
23    b.add_point(0, &[zero, zero, one]);
24    b.add_point(1, &[one, zero, zero]);
25    b.add_point(2, &[zero, one, zero]);
26    b.add_point(3, &[-one, zero, zero]);
27    b.add_point(4, &[zero, -one, zero]);
28    b.add_point(5, &[zero, zero, -one]);
29    let mut point_n = 6;
30
31    let mut cells = vec![
32        [0, 1, 2],
33        [0, 2, 3],
34        [0, 3, 4],
35        [0, 4, 1],
36        [5, 2, 1],
37        [5, 3, 2],
38        [5, 4, 3],
39        [5, 1, 4],
40    ];
41    let mut v = [[zero, zero, zero], [zero, zero, zero], [zero, zero, zero]];
42
43    for level in 0..refinement_level {
44        let mut edge_points = HashMap::new();
45        let mut new_cells = Vec::with_capacity(8 * usize::pow(6, level));
46        for c in &cells {
47            for (i, v_i) in v.iter_mut().enumerate() {
48                for (j, v_ij) in v_i.iter_mut().enumerate() {
49                    *v_ij = b.points[3 * c[i] + j];
50                }
51            }
52            let edges = [[1, 2], [0, 2], [0, 1]]
53                .iter()
54                .map(|[i, j]| {
55                    let mut pt_i = c[*i];
56                    let mut pt_j = c[*j];
57                    if pt_i > pt_j {
58                        std::mem::swap(&mut pt_i, &mut pt_j);
59                    }
60                    if let Vacant(e) = edge_points.entry((pt_i, pt_j)) {
61                        let v_i = v[*i];
62                        let v_j = v[*j];
63                        let mut new_pt = [
64                            half * (v_i[0] + v_j[0]),
65                            half * (v_i[1] + v_j[1]),
66                            half * (v_i[2] + v_j[2]),
67                        ];
68                        let size = (new_pt.iter().map(|x| x.powi(2)).sum::<T>()).sqrt();
69                        for i in new_pt.iter_mut() {
70                            *i /= size;
71                        }
72                        b.add_point(point_n, &new_pt);
73                        e.insert(point_n);
74                        point_n += 1;
75                    }
76                    edge_points[&(pt_i, pt_j)]
77                })
78                .collect::<Vec<_>>();
79            new_cells.push([c[0], edges[2], edges[1]]);
80            new_cells.push([c[1], edges[0], edges[2]]);
81            new_cells.push([c[2], edges[1], edges[0]]);
82            new_cells.push([edges[0], edges[1], edges[2]]);
83        }
84        cells = new_cells;
85    }
86    for (i, v) in cells.iter().enumerate() {
87        b.add_cell(i, v);
88    }
89}
90
91/// Create a surface grid of a regular sphere
92///
93/// A regular sphere is created by starting with a regular octahedron. The shape is then refined `refinement_level` times.
94/// Each time the grid is refined, each triangle is split into four triangles (by adding lines connecting the midpoints of
95/// each edge). The new points are then scaled so that they are a distance of 1 from the origin.
96pub fn regular_sphere<T: Scalar>(
97    refinement_level: u32,
98) -> SingleElementGrid<T, CiarletElement<T, IdentityMap, T>> {
99    let mut b = SingleElementGridBuilder::new_with_capacity(
100        3,
101        2 + usize::pow(4, refinement_level + 1),
102        8 * usize::pow(4, refinement_level),
103        (ReferenceCellType::Triangle, 1),
104    );
105    regular_sphere_add_points_and_cells(&mut b, refinement_level);
106    b.create_grid()
107}
108
109/// Create a grid of a regular sphere distributed in parallel
110#[cfg(feature = "mpi")]
111pub fn regular_sphere_distributed<T: Scalar + Equivalence, C: Communicator>(
112    comm: &C,
113    partitioner: GraphPartitioner,
114    refinement_level: u32,
115) -> ParallelGridImpl<'_, C, SingleElementGrid<T, CiarletElement<T, IdentityMap, T>>> {
116    let mut b = SingleElementGridBuilder::new(3, (ReferenceCellType::Triangle, 1));
117    if comm.rank() == 0 {
118        regular_sphere_add_points_and_cells(&mut b, refinement_level);
119        b.create_parallel_grid_root(comm, partitioner)
120    } else {
121        b.create_parallel_grid(comm, 0)
122    }
123}
124
125#[cfg(test)]
126mod test {
127    use super::*;
128    use crate::traits::{GeometryMap, Grid};
129    use approx::assert_relative_eq;
130    use rlst::rlst_dynamic_array;
131
132    #[test]
133    fn test_regular_sphere_0() {
134        let _g = regular_sphere::<f64>(0);
135    }
136
137    #[test]
138    fn test_regular_spheres() {
139        let _g1 = regular_sphere::<f64>(1);
140        let _g2 = regular_sphere::<f64>(2);
141        let _g3 = regular_sphere::<f64>(3);
142    }
143
144    #[test]
145    fn test_normal_is_outward() {
146        for i in 0..3 {
147            let g = regular_sphere::<f64>(i);
148            let mut points = rlst_dynamic_array!(f64, [2, 1]);
149            points[[0, 0]] = 1.0 / 3.0;
150            points[[1, 0]] = 1.0 / 3.0;
151            let map = g.geometry_map(ReferenceCellType::Triangle, 1, &points);
152            let mut mapped_pt = rlst_dynamic_array!(f64, [3, 1]);
153            let mut j = rlst_dynamic_array!(f64, [3, 2, 1]);
154            let mut jinv = rlst_dynamic_array!(f64, [2, 3, 1]);
155            let mut jdet = vec![0.0];
156            let mut normal = rlst_dynamic_array!(f64, [3, 1]);
157            for i in 0..g.entity_count(ReferenceCellType::Triangle) {
158                map.physical_points(i, &mut mapped_pt);
159                map.jacobians_inverses_dets_normals(i, &mut j, &mut jinv, &mut jdet, &mut normal);
160                let dot = mapped_pt
161                    .iter_value()
162                    .zip(normal.iter_value())
163                    .map(|(i, j)| i * j)
164                    .sum::<f64>();
165                assert!(dot > 0.0);
166            }
167        }
168    }
169
170    #[test]
171    fn test_normal_is_unit() {
172        for i in 0..3 {
173            let g = regular_sphere::<f64>(i);
174            let mut points = rlst_dynamic_array!(f64, [2, 1]);
175            points[[0, 0]] = 1.0 / 3.0;
176            points[[1, 0]] = 1.0 / 3.0;
177            let map = g.geometry_map(ReferenceCellType::Triangle, 1, &points);
178            let mut mapped_pt = rlst_dynamic_array!(f64, [3, 1]);
179            let mut j = rlst_dynamic_array!(f64, [3, 2, 1]);
180            let mut jinv = rlst_dynamic_array!(f64, [2, 3, 1]);
181            let mut jdet = vec![0.0];
182            let mut normal = rlst_dynamic_array!(f64, [3, 1]);
183            for i in 0..g.entity_count(ReferenceCellType::Triangle) {
184                map.physical_points(i, &mut mapped_pt);
185                map.jacobians_inverses_dets_normals(i, &mut j, &mut jinv, &mut jdet, &mut normal);
186                let dot = normal
187                    .iter_value()
188                    .zip(normal.iter_value())
189                    .map(|(i, j)| i * j)
190                    .sum::<f64>();
191                assert_relative_eq!(dot, 1.0, epsilon = 1e-10);
192            }
193        }
194    }
195}
ndgrid/shapes/regular_sphere.rs

ndgrid/shapes/
regular_sphere.rs
