test_mass_matrix/
test_mass_matrix.rs

1//! Test values in mass matrices against values computing using Bempp--cl and FEniCS
2
3use approx::assert_relative_eq;
4use ndelement::{
5    ciarlet::{LagrangeElementFamily, RaviartThomasElementFamily},
6    traits::{FiniteElement, MappedFiniteElement},
7    types::{Continuity, ReferenceCellType},
8};
9use ndfunctionspace::{FunctionSpaceImpl, traits::FunctionSpace};
10use ndgrid::{
11    shapes::regular_sphere,
12    traits::{Entity, GeometryMap, Grid},
13};
14use quadraturerules::{Domain, QuadratureRule, single_integral_quadrature};
15use rlst::{DynArray, rlst_dynamic_array};
16
17/// Test values in Lagrange mass matrix
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}
166
167/// Run tests
168fn main() {
169    println!("Testing Lagrange mass matrix");
170    test_lagrange_mass_matrix();
171
172    println!("Testing Raviart-Thomas mass matrix");
173    test_rt_mass_matrix();
174}
test_mass_matrix/test_mass_matrix.rs

test_mass_matrix/
test_mass_matrix.rs
