Add Gauss-Jacobi quadrature (#25)

* start adding quadrature

* Make Gauss-Jacobi quadrature work

* remove dependency

* hyperlink

* add tetrahedron

Cargo.toml

@@ -28,7 +28,6 @@ rlst = "0.2.0"
 serde = { version = "1", features = ["derive"], optional = true }
 
 [dev-dependencies]
-bempp = "0.1.0"
 paste = "1.*"
 approx = "0.5"
 

src/lib.rs

@@ -4,6 +4,7 @@
 
 pub mod ciarlet;
 pub mod polynomials;
+pub mod quadrature;
 pub mod reference_cell;
 pub mod traits;
 pub mod types;

src/polynomials.rs

@@ -500,124 +500,63 @@ pub fn tabulate_legendre_polynomials<
 #[cfg(test)]
 mod test {
     use super::*;
+    use crate::{quadrature::make_gauss_jacobi_quadrature, traits::QuadratureRule};
     use approx::*;
-    use bempp::quadrature::simplex_rules::simplex_rule;
-    use rlst::{rlst_dynamic_array2, rlst_dynamic_array3};
-
-    #[test]
-    fn test_legendre_interval() {
-        let degree = 6;
-
-        let rule = simplex_rule(
-            bempp::traits::types::ReferenceCellType::Interval,
-            degree + 1,
-        )
-        .unwrap();
-
-        let mut points = rlst_dynamic_array2!(f64, [1, rule.npoints]);
-        for (i, j) in rule.points.iter().enumerate() {
-            *points.get_mut([0, i]).unwrap() = *j;
-        }
-
-        let mut data = rlst_dynamic_array3!(
-            f64,
-            legendre_shape(ReferenceCellType::Interval, &points, degree, 0,)
-        );
-        tabulate_legendre_polynomials(ReferenceCellType::Interval, &points, degree, 0, &mut data);
-
-        for i in 0..degree + 1 {
-            for j in 0..degree + 1 {
-                let mut product = 0.0;
-                for k in 0..rule.npoints {
-                    product += data.get([0, i, k]).unwrap()
-                        * data.get([0, j, k]).unwrap()
-                        * rule.weights[k];
-                }
-                if i == j {
-                    assert_relative_eq!(product, 1.0, epsilon = 1e-12);
-                } else {
-                    assert_relative_eq!(product, 0.0, epsilon = 1e-12);
-                }
-            }
-        }
-    }
+    use paste::paste;
+    use rlst::{rlst_array_from_slice2, rlst_dynamic_array2, rlst_dynamic_array3};
+    macro_rules! test_orthogonal {
+        ($cell:ident, $degree:expr) => {
+            paste! {
+                #[test]
+                fn [<test_orthogonal_ $cell:lower _ $degree>]() {
+                    let rule = make_gauss_jacobi_quadrature(
+                        ReferenceCellType::[<$cell>],
+                        2 * [<$degree>],
+                    );
 
-    #[test]
-    fn test_legendre_triangle() {
-        let degree = 5;
-
-        let rule = simplex_rule(bempp::traits::types::ReferenceCellType::Triangle, 79).unwrap();
-        let mut points = rlst_dynamic_array2!(f64, [2, rule.npoints]);
-        for i in 0..rule.npoints {
-            for j in 0..2 {
-                *points.get_mut([j, i]).unwrap() = rule.points[i * 2 + j];
-            }
-        }
+                    let points = rlst_array_from_slice2!(rule.points(), [rule.dim(), rule.npoints()]);
 
-        let mut data = rlst_dynamic_array3!(
-            f64,
-            legendre_shape(ReferenceCellType::Triangle, &points, degree, 0,)
-        );
-        tabulate_legendre_polynomials(ReferenceCellType::Triangle, &points, degree, 0, &mut data);
-
-        for i in 0..data.shape()[1] {
-            for j in 0..data.shape()[1] {
-                let mut product = 0.0;
-                for k in 0..rule.npoints {
-                    product += data.get([0, i, k]).unwrap()
-                        * data.get([0, j, k]).unwrap()
-                        * rule.weights[k];
-                }
-                if i == j {
-                    assert_relative_eq!(product, 1.0, epsilon = 1e-12);
-                } else {
-                    assert_relative_eq!(product, 0.0, epsilon = 1e-12);
+                    let mut data = rlst_dynamic_array3!(
+                        f64,
+                        legendre_shape(ReferenceCellType::[<$cell>], &points, [<$degree>], 0,)
+                    );
+                    tabulate_legendre_polynomials(ReferenceCellType::[<$cell>], &points, [<$degree>], 0, &mut data);
+
+                    for i in 0..[<$degree>] + 1 {
+                        for j in 0..[<$degree>] + 1 {
+                            let mut product = 0.0;
+                            for k in 0..rule.npoints() {
+                                product += data.get([0, i, k]).unwrap()
+                                    * data.get([0, j, k]).unwrap()
+                                    * rule.weights()[k];
+                            }
+                            if i == j {
+                                assert_relative_eq!(product, 1.0, epsilon = 1e-12);
+                            } else {
+                                assert_relative_eq!(product, 0.0, epsilon = 1e-12);
+                            }
+                        }
+                    }
                 }
             }
-        }
+        };
     }
 
-    #[test]
-    fn test_legendre_quadrilateral() {
-        let degree = 5;
-
-        let rule =
-            simplex_rule(bempp::traits::types::ReferenceCellType::Quadrilateral, 85).unwrap();
-        let mut points = rlst_dynamic_array2!(f64, [2, rule.npoints]);
-        for i in 0..rule.npoints {
-            for j in 0..2 {
-                *points.get_mut([j, i]).unwrap() = rule.points[i * 2 + j];
-            }
-        }
-
-        let mut data = rlst_dynamic_array3!(
-            f64,
-            legendre_shape(ReferenceCellType::Quadrilateral, &points, degree, 0,)
-        );
-        tabulate_legendre_polynomials(
-            ReferenceCellType::Quadrilateral,
-            &points,
-            degree,
-            0,
-            &mut data,
-        );
-
-        for i in 0..data.shape()[1] {
-            for j in 0..data.shape()[1] {
-                let mut product = 0.0;
-                for k in 0..rule.npoints {
-                    product += data.get([0, i, k]).unwrap()
-                        * data.get([0, j, k]).unwrap()
-                        * rule.weights[k];
-                }
-                if i == j {
-                    assert_relative_eq!(product, 1.0, epsilon = 1e-12);
-                } else {
-                    assert_relative_eq!(product, 0.0, epsilon = 1e-12);
-                }
-            }
-        }
-    }
+    test_orthogonal!(Interval, 2);
+    test_orthogonal!(Interval, 3);
+    test_orthogonal!(Interval, 4);
+    test_orthogonal!(Interval, 5);
+    test_orthogonal!(Interval, 6);
+    test_orthogonal!(Triangle, 2);
+    test_orthogonal!(Triangle, 3);
+    test_orthogonal!(Triangle, 4);
+    test_orthogonal!(Triangle, 5);
+    test_orthogonal!(Triangle, 6);
+    test_orthogonal!(Quadrilateral, 2);
+    test_orthogonal!(Quadrilateral, 3);
+    test_orthogonal!(Quadrilateral, 4);
+    test_orthogonal!(Quadrilateral, 5);
+    test_orthogonal!(Quadrilateral, 6);
 
     #[test]
     fn test_legendre_interval_derivative() {

src/quadrature.rs

@@ -0,0 +1,44 @@
+//! Quadrature
+mod gauss_jacobi;
+pub use gauss_jacobi::make_gauss_jacobi_quadrature;
+
+use crate::traits::QuadratureRule as QuadratureRuleTrait;
+use rlst::RlstScalar;
+
+/// Quadrature rule
+pub struct QuadratureRule<T: RlstScalar<Real = T>> {
+    points: Vec<T>,
+    npoints: usize,
+    dim: usize,
+    weights: Vec<T>,
+}
+
+impl<T: RlstScalar<Real = T>> QuadratureRule<T> {
+    /// Create new
+    pub fn new(points: Vec<T>, weights: Vec<T>) -> Self {
+        let npoints = weights.len();
+        debug_assert!(points.len() % npoints == 0);
+        let dim = points.len() / npoints;
+        Self {
+            points,
+            npoints,
+            dim,
+            weights,
+        }
+    }
+}
+impl<T: RlstScalar<Real = T>> QuadratureRuleTrait for QuadratureRule<T> {
+    type T = T;
+    fn points(&self) -> &[T] {
+        &self.points
+    }
+    fn weights(&self) -> &[T] {
+        &self.weights
+    }
+    fn npoints(&self) -> usize {
+        self.npoints
+    }
+    fn dim(&self) -> usize {
+        self.dim
+    }
+}