pub fn regular_sphere<T: Scalar>(
refinement_level: u32,
) -> SingleElementGrid<T, CiarletElement<T, IdentityMap, T>>Expand description
Create a surface grid of a regular sphere
A regular sphere is created by starting with a regular octahedron. The shape is then refined refinement_level times.
Each time the grid is refined, each triangle is split into four triangles (by adding lines connecting the midpoints of
each edge). The new points are then scaled so that they are a distance of 1 from the origin.
Examples found in repository?
ndfunctionspace/examples/test_mass_matrix.rs (line 19)
18fn test_lagrange_mass_matrix() {
19 let grid = regular_sphere(0);
20
21 let family = LagrangeElementFamily::<f64>::new(1, Continuity::Standard);
22 let space = FunctionSpaceImpl::new(&grid, &family);
23
24 let mut mass_matrix = rlst_dynamic_array!(f64, [space.local_size(), space.local_size()]);
25
26 let element = &space.elements()[0];
27
28 let (p, w) = single_integral_quadrature(
29 QuadratureRule::XiaoGimbutas,
30 Domain::Triangle,
31 2 * element.lagrange_superdegree(),
32 )
33 .unwrap();
34 let npts = w.len();
35 let mut pts = rlst_dynamic_array!(f64, [2, npts]);
36 for i in 0..w.len() {
37 for j in 0..2 {
38 *pts.get_mut([j, i]).unwrap() = p[3 * i + j];
39 }
40 }
41 let wts = w.iter().map(|i| *i / 2.0).collect::<Vec<_>>();
42
43 let mut table = DynArray::<f64, 4>::from_shape(element.tabulate_array_shape(0, npts));
44 element.tabulate(&pts, 0, &mut table);
45
46 let gmap = grid.geometry_map(ReferenceCellType::Triangle, 1, &pts);
47 let mut jacobians = rlst_dynamic_array!(f64, [grid.geometry_dim(), grid.topology_dim(), npts]);
48 let mut jinv = rlst_dynamic_array!(f64, [grid.topology_dim(), grid.geometry_dim(), npts]);
49 let mut jdets = vec![0.0; npts];
50
51 for cell in grid.entity_iter(ReferenceCellType::Triangle) {
52 let dofs = space
53 .entity_closure_dofs(ReferenceCellType::Triangle, cell.local_index())
54 .unwrap();
55 gmap.jacobians_inverses_dets(cell.local_index(), &mut jacobians, &mut jinv, &mut jdets);
56 for (test_i, test_dof) in dofs.iter().enumerate() {
57 for (trial_i, trial_dof) in dofs.iter().enumerate() {
58 *mass_matrix.get_mut([*test_dof, *trial_dof]).unwrap() += wts
59 .iter()
60 .enumerate()
61 .map(|(i, w)| {
62 jdets[i]
63 * *w
64 * *table.get([0, i, test_i, 0]).unwrap()
65 * *table.get([0, i, trial_i, 0]).unwrap()
66 })
67 .sum::<f64>();
68 }
69 }
70 }
71
72 // Compare matrix entries to values from Bempp
73 for i in 0..6 {
74 assert_relative_eq!(mass_matrix[[i, i]], 0.5773502691896255, epsilon = 1e-10);
75 }
76 for i in 0..6 {
77 for j in 0..6 {
78 if i != j && mass_matrix[[i, j]].abs() > 0.001 {
79 assert_relative_eq!(mass_matrix[[i, j]], 0.1443375672974061, epsilon = 1e-10);
80 }
81 }
82 }
83}
84
85/// Test values in Raviart-Thomas mass matrix
86fn test_rt_mass_matrix() {
87 let grid = regular_sphere(0);
88
89 let family = RaviartThomasElementFamily::<f64>::new(1, Continuity::Standard);
90 let space = FunctionSpaceImpl::new(&grid, &family);
91
92 let mut mass_matrix = rlst_dynamic_array!(f64, [space.local_size(), space.local_size()]);
93
94 let element = &space.elements()[0];
95
96 let (p, w) = single_integral_quadrature(
97 QuadratureRule::XiaoGimbutas,
98 Domain::Triangle,
99 2 * element.lagrange_superdegree(),
100 )
101 .unwrap();
102 let npts = w.len();
103 let mut pts = rlst_dynamic_array!(f64, [2, npts]);
104 for i in 0..w.len() {
105 for j in 0..2 {
106 *pts.get_mut([j, i]).unwrap() = p[3 * i + j];
107 }
108 }
109 let wts = w.iter().map(|i| *i / 2.0).collect::<Vec<_>>();
110
111 let mut table = DynArray::<f64, 4>::from_shape(element.tabulate_array_shape(0, npts));
112 element.tabulate(&pts, 0, &mut table);
113 let mut pushed_table = rlst_dynamic_array!(
114 f64,
115 [table.shape()[0], table.shape()[1], table.shape()[2], 3]
116 );
117
118 let gmap = grid.geometry_map(ReferenceCellType::Triangle, 1, &pts);
119 let mut jacobians = rlst_dynamic_array!(f64, [grid.geometry_dim(), grid.topology_dim(), npts]);
120 let mut jinv = rlst_dynamic_array!(f64, [grid.topology_dim(), grid.geometry_dim(), npts]);
121 let mut jdets = vec![0.0; npts];
122
123 for cell in grid.entity_iter(ReferenceCellType::Triangle) {
124 let dofs = space
125 .entity_closure_dofs(ReferenceCellType::Triangle, cell.local_index())
126 .unwrap();
127 gmap.jacobians_inverses_dets(cell.local_index(), &mut jacobians, &mut jinv, &mut jdets);
128 element.push_forward(&table, 0, &jacobians, &jdets, &jinv, &mut pushed_table);
129
130 for (test_i, test_dof) in dofs.iter().enumerate() {
131 for (trial_i, trial_dof) in dofs.iter().enumerate() {
132 *mass_matrix.get_mut([*test_dof, *trial_dof]).unwrap() += wts
133 .iter()
134 .enumerate()
135 .map(|(i, w)| {
136 jdets[i]
137 * *w
138 * (0..3)
139 .map(|j| {
140 *pushed_table.get([0, i, test_i, j]).unwrap()
141 * *pushed_table.get([0, i, trial_i, j]).unwrap()
142 })
143 .sum::<f64>()
144 })
145 .sum::<f64>();
146 }
147 }
148 }
149
150 // Compare matrix entries to FEniCS
151 for i in 0..12 {
152 assert_relative_eq!(mass_matrix[[i, i]], 0.4811252243246884, epsilon = 1e-10);
153 }
154 for i in 0..12 {
155 for j in 0..12 {
156 if i != j && mass_matrix[[i, j]].abs() > 0.001 {
157 assert_relative_eq!(
158 mass_matrix[[i, j]].abs(),
159 0.0481125224324689,
160 epsilon = 1e-10
161 );
162 }
163 }
164 }
165}