diff --git a/Cargo.toml b/Cargo.toml index 2af9ed3..2b7ee4c 100644 --- a/Cargo.toml +++ b/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" diff --git a/src/lib.rs b/src/lib.rs index be9bbbf..305aee9 100644 --- a/src/lib.rs +++ b/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; diff --git a/src/polynomials.rs b/src/polynomials.rs index 0cd8df6..ca17eb3 100644 --- a/src/polynomials.rs +++ b/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 []() { + 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() { diff --git a/src/quadrature.rs b/src/quadrature.rs new file mode 100644 index 0000000..00bf00a --- /dev/null +++ b/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> { + points: Vec, + npoints: usize, + dim: usize, + weights: Vec, +} + +impl> QuadratureRule { + /// Create new + pub fn new(points: Vec, weights: Vec) -> Self { + let npoints = weights.len(); + debug_assert!(points.len() % npoints == 0); + let dim = points.len() / npoints; + Self { + points, + npoints, + dim, + weights, + } + } +} +impl> QuadratureRuleTrait for QuadratureRule { + 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 + } +} diff --git a/src/quadrature/gauss_jacobi.rs b/src/quadrature/gauss_jacobi.rs new file mode 100644 index 0000000..621daea --- /dev/null +++ b/src/quadrature/gauss_jacobi.rs @@ -0,0 +1,259 @@ +//! Gauss-Jacobi quadrature +//! +//! Adapted from the C++ code by Chris Richardson +//! +use super::QuadratureRule; +use crate::traits::QuadratureRule as QuadratureRuleTrait; +use crate::types::{Array2D, ReferenceCellType}; +use itertools::izip; +use num::traits::FloatConst; +use rlst::{rlst_dynamic_array2, rlst_dynamic_array3, DefaultIteratorMut, RlstScalar}; +use std::cmp::PartialOrd; + +/// Evaluate the nth Jacobi polynomial and derivatives with weight +/// parameters (a, 0) at points x +fn compute_jacobi_deriv>( + a: T, + n: usize, + nderiv: usize, + x: &[T], +) -> Array2D { + let mut j_all = rlst_dynamic_array3!(T, [nderiv + 1, n + 1, x.len()]); + + let one = T::from(1.0).unwrap(); + let two = T::from(2.0).unwrap(); + + for i in 0..=nderiv { + if i == 0 { + for j in 0..x.len() { + j_all[[i, 0, j]] = one; + } + } + + if n > 0 { + if i == 0 { + for (j, x_j) in x.iter().enumerate() { + j_all[[i, 1, j]] = (*x_j * (a + two) + a) / two; + } + } else if i == 1 { + for j in 0..x.len() { + j_all[[i, 1, j]] = a / two + one; + } + } + } + for j in 2..=n { + let j_t = T::from(j).unwrap(); + let a1 = two * j_t * (j_t + a) * (two * j_t + a - two); + let a2 = (two * j_t + a - one) * (a * a) / a1; + let a3 = (two * j_t + a - one) * (two * j_t + a) / (two * j_t * (j_t + a)); + let a4 = two * (j_t + a - one) * (j_t - one) * (two * j_t + a) / a1; + for (k, x_k) in x.iter().enumerate() { + j_all[[i, j, k]] = + j_all[[i, j - 1, k]] * (*x_k * a3 + a2) - j_all[[i, j - 2, k]] * a4; + } + if i > 0 { + for (k, _) in x.iter().enumerate() { + let add = T::from(i).unwrap() * a3 * j_all[[i - 1, j - 1, k]]; + j_all[[i, j, k]] += add; + } + } + } + } + + let mut j = rlst_dynamic_array2!(T, [nderiv + 1, x.len()]); + for (i1, mut col) in j.col_iter_mut().enumerate() { + for (i0, entry) in col.iter_mut().enumerate() { + *entry = j_all[[i0, n, i1]]; + } + } + + j +} + +/// Computes the m roots of \f$P_{m}^{a,0}\f$ on [-1,1] by Newton's +/// method. The initial guesses are the Chebyshev points. Algorithm +/// implemented from the pseudocode given by Karniadakis and Sherwin. +fn compute_gauss_jacobi_points + FloatConst + PartialOrd>( + a: T, + m: usize, +) -> Vec { + let eps = T::from(1.0e-8).unwrap(); + let max_iter = 100; + let two = T::from(2.0).unwrap(); + let one = T::from(1.0).unwrap(); + + let mut x = vec![T::zero(); m]; + + for k in 0..m { + // Initial guess + x[k] = -T::cos(T::from(2 * k + 1).unwrap() * T::PI() / T::from(2 * m).unwrap()); + if k > 0 { + x[k] = (x[k] + x[k - 1]) / two; + } + + for _ in 0..max_iter { + let s = x.iter().take(k).map(|i| one / (x[k] - *i)).sum::(); + let f = compute_jacobi_deriv(a, m, 1, &x[k..k + 1]); + let delta = f[[0, 0]] / (f[[1, 0]] - f[[0, 0]] * s); + x[k] -= delta; + if delta.abs() < eps { + break; + } + } + } + x +} + +/// Note: computes on [-1, 1] +fn compute_gauss_jacobi_rule + FloatConst + PartialOrd>( + a: T, + m: usize, +) -> (Vec, Vec) { + let one = T::from(1.0).unwrap(); + let two = T::from(2.0).unwrap(); + + let pts = compute_gauss_jacobi_points(a, m); + let j_d = compute_jacobi_deriv(a, m, 1, &pts); + let a1 = T::pow(two, a + one); + let wts = pts + .iter() + .enumerate() + .map(|(i, x)| a1 / (one - x.powi(2)) / j_d[[1, i]].powi(2)) + .collect::>(); + + (pts, wts) +} + +fn make_quadrature_line + FloatConst + PartialOrd>( + m: usize, +) -> QuadratureRule { + let (mut pts, mut wts) = compute_gauss_jacobi_rule(T::zero(), m); + + let half = T::from(0.5).unwrap(); + let one = T::from(1.0).unwrap(); + for p in pts.iter_mut() { + *p = half * (*p + one); + } + for w in wts.iter_mut() { + *w *= half; + } + QuadratureRule::new(pts, wts) +} + +fn make_quadrature_triangle_collapsed + FloatConst + PartialOrd>( + m: usize, +) -> QuadratureRule { + let one = T::from(1.0).unwrap(); + + let (ptx, wtx) = compute_gauss_jacobi_rule(T::zero(), m); + let (pty, wty) = compute_gauss_jacobi_rule(T::from(1).unwrap(), m); + + let mut pts = vec![T::zero(); m.pow(2) * 2]; + let mut wts = vec![T::zero(); m.pow(2)]; + + for (i, (px, wx)) in izip!(&ptx, &wtx).enumerate() { + for (j, (py, wy)) in izip!(&pty, &wty).enumerate() { + let index = i * wty.len() + j; + pts[2 * index] = T::from(0.25).unwrap() * (one + *px) * (one - *py); + pts[2 * index + 1] = T::from(0.5).unwrap() * (one + *py); + wts[index] = *wx * *wy * T::from(0.125).unwrap(); + } + } + QuadratureRule::new(pts, wts) +} + +fn make_quadrature_tetrahedron_collapsed + FloatConst + PartialOrd>( + m: usize, +) -> QuadratureRule { + let one = T::from(1.0).unwrap(); + + let (ptx, wtx) = compute_gauss_jacobi_rule(T::zero(), m); + let (pty, wty) = compute_gauss_jacobi_rule(T::from(1).unwrap(), m); + let (ptz, wtz) = compute_gauss_jacobi_rule(T::from(2).unwrap(), m); + + let mut pts = vec![T::zero(); m.pow(3) * 3]; + let mut wts = vec![T::zero(); m.pow(3)]; + + for (i, (px, wx)) in izip!(&ptx, &wtx).enumerate() { + for (j, (py, wy)) in izip!(&pty, &wty).enumerate() { + for (k, (pz, wz)) in izip!(&ptz, &wtz).enumerate() { + let index = i * wty.len() * wtz.len() + j * wtz.len() + k; + pts[3 * index] = T::from(0.125).unwrap() * (one + *px) * (one - *py) * (one - *pz); + pts[3 * index + 1] = T::from(0.25).unwrap() * (one + *py) * (one - *pz); + pts[3 * index + 1] = T::from(0.5).unwrap() * (one + *pz); + wts[index] = *wx * *wy * *wz * T::from(0.015625).unwrap(); + } + } + } + QuadratureRule::new(pts, wts) +} + +/// Get the points and weights for a Gauss-Jacobi quadrature rule +pub fn make_gauss_jacobi_quadrature + FloatConst + PartialOrd>( + celltype: ReferenceCellType, + m: usize, +) -> QuadratureRule { + let np = (m + 2) / 2; + match celltype { + ReferenceCellType::Interval => make_quadrature_line::(np), + ReferenceCellType::Quadrilateral => { + let rule = make_quadrature_line::(np); + let mut pts = vec![T::zero(); np.pow(2) * 2]; + let mut wts = vec![T::zero(); np.pow(2)]; + + for (i, (pi, wi)) in izip!(rule.points(), rule.weights()).enumerate() { + for (j, (pj, wj)) in izip!(rule.points(), rule.weights()).enumerate() { + let index = i * np + j; + pts[2 * index] = *pi; + pts[2 * index + 1] = *pj; + wts[index] = *wi * *wj; + } + } + QuadratureRule::new(pts, wts) + } + ReferenceCellType::Hexahedron => { + let rule = make_quadrature_line::(np); + let mut pts = vec![T::zero(); np.pow(3) * 3]; + let mut wts = vec![T::zero(); np.pow(3)]; + + for (i, (pi, wi)) in izip!(rule.points(), rule.weights()).enumerate() { + for (j, (pj, wj)) in izip!(rule.points(), rule.weights()).enumerate() { + for (k, (pk, wk)) in izip!(rule.points(), rule.weights()).enumerate() { + let index = i * np.pow(2) + j * np + k; + pts[3 * index] = *pi; + pts[3 * index + 1] = *pj; + pts[3 * index + 2] = *pk; + wts[index] = *wi * *wj * *wk; + } + } + } + QuadratureRule::new(pts, wts) + } + ReferenceCellType::Triangle => make_quadrature_triangle_collapsed::(np), + ReferenceCellType::Tetrahedron => make_quadrature_tetrahedron_collapsed::(np), + _ => { + panic!("Unsupported cell type"); + } + } +} + +#[cfg(test)] +mod test { + use super::*; + use approx::*; + + #[test] + fn test_interval_3() { + let rule = make_quadrature_line::(3); + + for (p, q) in izip!(rule.points(), [0.1127016653792583, 0.5, 0.8872983346207417]) { + assert_relative_eq!(*p, q); + } + for (w, v) in izip!( + rule.weights(), + [0.2777777777777777, 0.4444444444444444, 0.2777777777777777] + ) { + assert_relative_eq!(*w, v); + } + } +} diff --git a/src/traits.rs b/src/traits.rs index 01e8d37..ee7eadb 100644 --- a/src/traits.rs +++ b/src/traits.rs @@ -60,3 +60,17 @@ pub trait ElementFamily { /// Get an elenent for a cell type fn element(&self, cell_type: Self::CellType) -> Self::FiniteElement; } + +pub trait QuadratureRule { + //! A quadrature rule + /// The scalar type + type T: RlstScalar; + /// Quadrature points + fn points(&self) -> &[Self::T]; + /// Quadrature weights + fn weights(&self) -> &[Self::T]; + /// Number of quadrature points + fn npoints(&self) -> usize; + /// Topological dimension of cell (ie number of components of each point) + fn dim(&self) -> usize; +} diff --git a/src/types.rs b/src/types.rs index cd6d349..2731014 100644 --- a/src/types.rs +++ b/src/types.rs @@ -1,6 +1,12 @@ //! Types #[cfg(feature = "mpi")] use mpi::traits::Equivalence; +use rlst::{Array, BaseArray, VectorContainer}; + +/// An N-dimensional array +pub type ArrayND = Array, N>, N>; +/// A 2-dimensional array +pub type Array2D = ArrayND<2, T>; /// Continuity type #[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]