ndgrid/shapes/
regular_sphere.rs

1//! Regular sphere grid
2
3use crate::{
4    grid::local_grid::{SingleElementGrid, SingleElementGridBuilder},
5    traits::Builder,
6    types::Scalar,
7};
8use ndelement::{ciarlet::CiarletElement, map::IdentityMap, types::ReferenceCellType};
9use std::collections::{HashMap, hash_map::Entry::Vacant};
10
11/// Create a surface grid of a regular sphere
12///
13/// A regular sphere is created by starting with a regular octahedron. The shape is then refined `refinement_level` times.
14/// Each time the grid is refined, each triangle is split into four triangles (by adding lines connecting the midpoints of
15/// each edge). The new points are then scaled so that they are a distance of 1 from the origin.
16pub fn regular_sphere<T: Scalar>(
17    refinement_level: u32,
18) -> SingleElementGrid<T, CiarletElement<T, IdentityMap, T>> {
19    let mut b = SingleElementGridBuilder::new_with_capacity(
20        3,
21        2 + usize::pow(4, refinement_level + 1),
22        8 * usize::pow(4, refinement_level),
23        (ReferenceCellType::Triangle, 1),
24    );
25    let zero = T::from(0.0).unwrap();
26    let one = T::from(1.0).unwrap();
27    let half = T::from(0.5).unwrap();
28    b.add_point(0, &[zero, zero, one]);
29    b.add_point(1, &[one, zero, zero]);
30    b.add_point(2, &[zero, one, zero]);
31    b.add_point(3, &[-one, zero, zero]);
32    b.add_point(4, &[zero, -one, zero]);
33    b.add_point(5, &[zero, zero, -one]);
34    let mut point_n = 6;
35
36    let mut cells = vec![
37        [0, 1, 2],
38        [0, 2, 3],
39        [0, 3, 4],
40        [0, 4, 1],
41        [5, 2, 1],
42        [5, 3, 2],
43        [5, 4, 3],
44        [5, 1, 4],
45    ];
46    let mut v = [[zero, zero, zero], [zero, zero, zero], [zero, zero, zero]];
47
48    for level in 0..refinement_level {
49        let mut edge_points = HashMap::new();
50        let mut new_cells = Vec::with_capacity(8 * usize::pow(6, level));
51        for c in &cells {
52            for (i, v_i) in v.iter_mut().enumerate() {
53                for (j, v_ij) in v_i.iter_mut().enumerate() {
54                    *v_ij = b.points[3 * c[i] + j];
55                }
56            }
57            let edges = [[1, 2], [0, 2], [0, 1]]
58                .iter()
59                .map(|[i, j]| {
60                    let mut pt_i = c[*i];
61                    let mut pt_j = c[*j];
62                    if pt_i > pt_j {
63                        std::mem::swap(&mut pt_i, &mut pt_j);
64                    }
65                    if let Vacant(e) = edge_points.entry((pt_i, pt_j)) {
66                        let v_i = v[*i];
67                        let v_j = v[*j];
68                        let mut new_pt = [
69                            half * (v_i[0] + v_j[0]),
70                            half * (v_i[1] + v_j[1]),
71                            half * (v_i[2] + v_j[2]),
72                        ];
73                        let size = (new_pt.iter().map(|x| x.powi(2)).sum::<T>()).sqrt();
74                        for i in new_pt.iter_mut() {
75                            *i /= size;
76                        }
77                        b.add_point(point_n, &new_pt);
78                        e.insert(point_n);
79                        point_n += 1;
80                    }
81                    edge_points[&(pt_i, pt_j)]
82                })
83                .collect::<Vec<_>>();
84            new_cells.push([c[0], edges[2], edges[1]]);
85            new_cells.push([c[1], edges[0], edges[2]]);
86            new_cells.push([c[2], edges[1], edges[0]]);
87            new_cells.push([edges[0], edges[1], edges[2]]);
88        }
89        cells = new_cells;
90    }
91    for (i, v) in cells.iter().enumerate() {
92        b.add_cell(i, v);
93    }
94
95    b.create_grid()
96}
97
98#[cfg(test)]
99mod test {
100    use super::*;
101    use crate::traits::{GeometryMap, Grid};
102    use approx::assert_relative_eq;
103
104    #[test]
105    fn test_regular_sphere_0() {
106        let _g = regular_sphere::<f64>(0);
107    }
108
109    #[test]
110    fn test_regular_spheres() {
111        let _g1 = regular_sphere::<f64>(1);
112        let _g2 = regular_sphere::<f64>(2);
113        let _g3 = regular_sphere::<f64>(3);
114    }
115
116    #[test]
117    fn test_normal_is_outward() {
118        for i in 0..3 {
119            let g = regular_sphere::<f64>(i);
120            let points = vec![1.0 / 3.0, 1.0 / 3.0];
121            let map = g.geometry_map(ReferenceCellType::Triangle, 1, &points);
122            let mut mapped_pt = vec![0.0; 3];
123            let mut j = vec![0.0; 6];
124            let mut jdet = vec![0.0];
125            let mut normal = vec![0.0; 3];
126            for i in 0..g.entity_count(ReferenceCellType::Triangle) {
127                map.physical_points(i, &mut mapped_pt);
128                map.jacobians_dets_normals(i, &mut j, &mut jdet, &mut normal);
129                let dot = mapped_pt
130                    .iter()
131                    .zip(&normal)
132                    .map(|(i, j)| i * j)
133                    .sum::<f64>();
134                assert!(dot > 0.0);
135            }
136        }
137    }
138
139    #[test]
140    fn test_normal_is_unit() {
141        for i in 0..3 {
142            let g = regular_sphere::<f64>(i);
143            let points = vec![1.0 / 3.0, 1.0 / 3.0];
144            let map = g.geometry_map(ReferenceCellType::Triangle, 1, &points);
145            let mut j = vec![0.0; 6];
146            let mut jdet = vec![0.0];
147            let mut mapped_pt = vec![0.0; 3];
148            let mut normal = vec![0.0; 3];
149            for i in 0..g.entity_count(ReferenceCellType::Triangle) {
150                map.physical_points(i, &mut mapped_pt);
151                map.jacobians_dets_normals(i, &mut j, &mut jdet, &mut normal);
152                let dot = normal.iter().zip(&normal).map(|(i, j)| i * j).sum::<f64>();
153                assert_relative_eq!(dot, 1.0, epsilon = 1e-10);
154            }
155        }
156    }
157}
ndgrid/shapes/regular_sphere.rs

ndgrid/shapes/
regular_sphere.rs
