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