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, LagrangeVariant, RaviartThomasElementFamily},
6    traits::{FiniteElement, MappedFiniteElement},
7    types::{Continuity, ReferenceCellType},
8};
9use ndfunctionspace::{FunctionSpaceImpl, traits::FunctionSpace};
10use ndmesh::{
11    shapes::regular_sphere,
12    traits::{Entity, GeometryMap, Mesh},
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 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}
167
168/// Run tests
169fn main() {
170    println!("Testing Lagrange mass matrix");
171    test_lagrange_mass_matrix();
172
173    println!("Testing Raviart-Thomas mass matrix");
174    test_rt_mass_matrix();
175}
test_mass_matrix/test_mass_matrix.rs

test_mass_matrix/
test_mass_matrix.rs
