From 27d1f4f2e99a445e827e592091bccd14351ad10b Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Sat, 22 Jun 2024 14:45:59 +0100 Subject: [PATCH] Copy elements from bempp-rs (#2) * Add FiniteElement, ElementFamily * Copy elements from bempp-rs * use quadrature from bempp in test --- Cargo.toml | 3 + src/ciarlet.rs | 713 ++++++++++++++++++++++++ src/ciarlet/lagrange.rs | 214 ++++++++ src/ciarlet/raviart_thomas.rs | 135 +++++ src/lib.rs | 12 + src/polynomials.rs | 981 ++++++++++++++++++++++++++++++++++ src/reference_cell.rs | 522 ++++++++++++++++++ src/traits.rs | 5 + src/traits/element.rs | 59 ++ src/types.rs | 112 ++++ 10 files changed, 2756 insertions(+) create mode 100644 src/ciarlet.rs create mode 100644 src/ciarlet/lagrange.rs create mode 100644 src/ciarlet/raviart_thomas.rs create mode 100644 src/polynomials.rs create mode 100644 src/reference_cell.rs create mode 100644 src/traits.rs create mode 100644 src/traits/element.rs create mode 100644 src/types.rs diff --git a/Cargo.toml b/Cargo.toml index 3e9a4ba..fdfb332 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -34,6 +34,9 @@ rand = "0.8.5" rlst = "0.1.0" thiserror="1.*" +[dev-dependencies] +bempp = "0.1.0" + [package.metadata.docs.rs] cargo-args = ["-Zunstable-options", "-Zrustdoc-scrape-examples"] diff --git a/src/ciarlet.rs b/src/ciarlet.rs new file mode 100644 index 0000000..1b9a5b2 --- /dev/null +++ b/src/ciarlet.rs @@ -0,0 +1,713 @@ +//! Finite element definitions + +use crate::polynomials::{legendre_shape, polynomial_count, tabulate_legendre_polynomials}; +use crate::reference_cell; +use crate::traits::FiniteElement; +use crate::types::{Continuity, MapType, ReferenceCellType}; +use rlst::{ + dense::array::views::ArrayViewMut, rlst_dynamic_array2, rlst_dynamic_array3, Array, BaseArray, + MatrixInverse, RandomAccessByRef, RandomAccessMut, RlstScalar, Shape, VectorContainer, +}; + +pub mod lagrange; +pub mod raviart_thomas; +pub use lagrange::LagrangeElementFamily; +pub use raviart_thomas::RaviartThomasElementFamily; + +type EntityPoints = [Vec, 2>, 2>>; 4]; +type EntityWeights = [Vec, 3>, 3>>; 4]; + +/// Compute the number of derivatives for a cell +fn compute_derivative_count(nderivs: usize, cell_type: ReferenceCellType) -> usize { + match cell_type { + ReferenceCellType::Point => 0, + ReferenceCellType::Interval => nderivs + 1, + ReferenceCellType::Triangle => (nderivs + 1) * (nderivs + 2) / 2, + ReferenceCellType::Quadrilateral => (nderivs + 1) * (nderivs + 2) / 2, + ReferenceCellType::Tetrahedron => (nderivs + 1) * (nderivs + 2) * (nderivs + 3) / 6, + ReferenceCellType::Hexahedron => (nderivs + 1) * (nderivs + 2) * (nderivs + 3) / 6, + ReferenceCellType::Prism => (nderivs + 1) * (nderivs + 2) * (nderivs + 3) / 6, + ReferenceCellType::Pyramid => (nderivs + 1) * (nderivs + 2) * (nderivs + 3) / 6, + } +} + +/// A Ciarlet element +pub struct CiarletElement { + cell_type: ReferenceCellType, + degree: usize, + embedded_superdegree: usize, + map_type: MapType, + value_shape: Vec, + value_size: usize, + continuity: Continuity, + dim: usize, + coefficients: Array, 3>, 3>, + entity_dofs: [Vec>; 4], + interpolation_points: EntityPoints, + interpolation_weights: EntityWeights, +} + +impl CiarletElement +where + for<'a> Array, 2>, 2>, 2>: MatrixInverse, +{ + /// Create a Ciarlet element + #[allow(clippy::too_many_arguments)] + pub fn create( + cell_type: ReferenceCellType, + degree: usize, + value_shape: Vec, + polynomial_coeffs: Array, 3>, 3>, + interpolation_points: EntityPoints, + interpolation_weights: EntityWeights, + map_type: MapType, + continuity: Continuity, + embedded_superdegree: usize, + ) -> Self { + let mut dim = 0; + let mut npts = 0; + + for emats in &interpolation_weights { + for mat in emats { + dim += mat.shape()[0]; + npts += mat.shape()[2]; + } + } + let tdim = reference_cell::dim(cell_type); + + let mut value_size = 1; + for i in &value_shape { + value_size *= *i; + } + + for matrices in &interpolation_weights { + for mat in matrices { + if mat.shape()[1] != value_size { + panic!("Incompatible value size"); + } + } + } + + let new_pts = if continuity == Continuity::Discontinuous { + let mut new_pts: EntityPoints = [vec![], vec![], vec![], vec![]]; + let mut all_pts = rlst_dynamic_array2![T::Real, [npts, tdim]]; + for (i, pts_i) in interpolation_points.iter().take(tdim).enumerate() { + for _pts in pts_i { + new_pts[i].push(rlst_dynamic_array2![T::Real, [0, tdim]]); + } + } + let mut row = 0; + for pts_i in interpolation_points.iter() { + for pts in pts_i { + let nrows = pts.shape()[0]; + all_pts + .view_mut() + .into_subview([row, 0], [nrows, tdim]) + .fill_from(pts.view()); + row += nrows; + } + } + new_pts[tdim].push(all_pts); + new_pts + } else { + interpolation_points + }; + let new_wts = if continuity == Continuity::Discontinuous { + let mut new_wts = [vec![], vec![], vec![], vec![]]; + let mut pn = 0; + let mut dn = 0; + let mut all_mat = rlst_dynamic_array3!(T, [dim, value_size, npts]); + for (i, mi) in interpolation_weights.iter().take(tdim).enumerate() { + for _mat in mi { + new_wts[i].push(rlst_dynamic_array3!(T, [0, value_size, 0])); + } + } + for mi in interpolation_weights.iter() { + for mat in mi { + let dim0 = mat.shape()[0]; + let dim2 = mat.shape()[2]; + all_mat + .view_mut() + .into_subview([dn, 0, pn], [dim0, value_size, dim2]) + .fill_from(mat.view()); + dn += dim0; + pn += dim2; + } + } + new_wts[tdim].push(all_mat); + new_wts + } else { + interpolation_weights + }; + + // Compute the dual matrix + let pdim = polynomial_count(cell_type, embedded_superdegree); + let mut d_matrix = rlst_dynamic_array3!(T, [value_size, pdim, dim]); + + let mut dof = 0; + for d in 0..4 { + for (e, pts) in new_pts[d].iter().enumerate() { + if pts.shape()[0] > 0 { + let mut table = rlst_dynamic_array3!(T, [1, pdim, pts.shape()[0]]); + tabulate_legendre_polynomials( + cell_type, + pts, + embedded_superdegree, + 0, + &mut table, + ); + let mat = &new_wts[d][e]; + for i in 0..mat.shape()[0] { + for l in 0..pdim { + for j in 0..value_size { + // d_matrix[j, l, dof + i] = inner(mat[i, j, :], table[0, l, :]) + *d_matrix.get_mut([j, l, dof + i]).unwrap() = mat + .view() + .slice(0, i) + .slice(0, j) + .inner(table.view().slice(0, 0).slice(0, l)); + } + } + } + dof += mat.shape()[0]; + } + } + } + + let mut inverse = rlst::rlst_dynamic_array2!(T, [dim, dim]); + + for i in 0..dim { + for j in 0..dim { + let entry = inverse.get_mut([i, j]).unwrap(); + *entry = T::from(0.0).unwrap(); + for k in 0..value_size { + for l in 0..pdim { + *entry += *polynomial_coeffs.get([i, k, l]).unwrap() + * *d_matrix.get([k, l, j]).unwrap(); + } + } + } + } + + inverse.view_mut().into_inverse_alloc().unwrap(); + + let mut coefficients = rlst_dynamic_array3!(T, [dim, value_size, pdim]); + for i in 0..dim { + for j in 0..value_size { + for k in 0..pdim { + // coefficients[i, j, k] = inner(inverse[i, :], polynomial_coeffs[:, j, k]) + *coefficients.get_mut([i, j, k]).unwrap() = inverse + .view() + .slice(0, i) + .inner(polynomial_coeffs.view().slice(1, j).slice(1, k)); + } + } + } + + let mut entity_dofs = [vec![], vec![], vec![], vec![]]; + let mut dof = 0; + for i in 0..4 { + for pts in &new_pts[i] { + let dofs = (dof..dof + pts.shape()[0]).collect::>(); + entity_dofs[i].push(dofs); + dof += pts.shape()[0]; + } + } + CiarletElement:: { + cell_type, + degree, + embedded_superdegree, + map_type, + value_shape, + value_size, + continuity, + dim, + coefficients, + entity_dofs, + interpolation_points: new_pts, + interpolation_weights: new_wts, + } + } + + /// The polynomial degree + pub fn degree(&self) -> usize { + self.degree + } + /// The continuity of the element + pub fn continuity(&self) -> Continuity { + self.continuity + } + /// The interpolation points + pub fn interpolation_points(&self) -> &EntityPoints { + &self.interpolation_points + } + /// The interpolation weights + pub fn interpolation_weights(&self) -> &EntityWeights { + &self.interpolation_weights + } +} +impl FiniteElement for CiarletElement { + type CellType = ReferenceCellType; + type MapType = MapType; + type T = T; + fn value_shape(&self) -> &[usize] { + &self.value_shape + } + fn value_size(&self) -> usize { + self.value_size + } + fn map_type(&self) -> MapType { + self.map_type + } + + fn cell_type(&self) -> ReferenceCellType { + self.cell_type + } + fn embedded_superdegree(&self) -> usize { + self.embedded_superdegree + } + fn dim(&self) -> usize { + self.dim + } + fn tabulate + Shape<2>>( + &self, + points: &Array2, + nderivs: usize, + data: &mut impl RandomAccessMut<4, Item = T>, + ) { + let mut table = rlst_dynamic_array3!( + T, + legendre_shape(self.cell_type, points, self.embedded_superdegree, nderivs,) + ); + tabulate_legendre_polynomials( + self.cell_type, + points, + self.embedded_superdegree, + nderivs, + &mut table, + ); + + for d in 0..table.shape()[0] { + for p in 0..points.shape()[0] { + for j in 0..self.value_size { + for b in 0..self.dim { + // data[d, p, b, j] = inner(self.coefficients[b, j, :], table[d, :, p]) + *data.get_mut([d, p, b, j]).unwrap() = self + .coefficients + .view() + .slice(0, b) + .slice(0, j) + .inner(table.view().slice(0, d).slice(1, p)); + } + } + } + } + } + fn entity_dofs(&self, entity_dim: usize, entity_number: usize) -> Option<&[usize]> { + if entity_dim < 4 && entity_number < self.entity_dofs[entity_dim].len() { + Some(&self.entity_dofs[entity_dim][entity_number]) + } else { + None + } + } + fn tabulate_array_shape(&self, nderivs: usize, npoints: usize) -> [usize; 4] { + let deriv_count = compute_derivative_count(nderivs, self.cell_type()); + let point_count = npoints; + let basis_count = self.dim(); + let value_size = self.value_size(); + [deriv_count, point_count, basis_count, value_size] + } +} + +#[cfg(test)] +mod test { + use super::*; + use approx::*; + use rlst::rlst_dynamic_array4; + + fn check_dofs(e: impl FiniteElement) { + let mut ndofs = 0; + for (dim, entity_count) in match e.cell_type() { + ReferenceCellType::Point => vec![1], + ReferenceCellType::Interval => vec![2, 1], + ReferenceCellType::Triangle => vec![3, 3, 1], + ReferenceCellType::Quadrilateral => vec![4, 4, 1], + ReferenceCellType::Tetrahedron => vec![4, 6, 4, 1], + ReferenceCellType::Hexahedron => vec![8, 12, 6, 1], + ReferenceCellType::Prism => vec![6, 9, 5, 1], + ReferenceCellType::Pyramid => vec![5, 8, 5, 1], + } + .iter() + .enumerate() + { + for entity in 0..*entity_count { + ndofs += e.entity_dofs(dim, entity).unwrap().len(); + } + } + assert_eq!(ndofs, e.dim()); + } + + #[test] + fn test_lagrange_1() { + let e = lagrange::create::(ReferenceCellType::Triangle, 1, Continuity::Continuous); + assert_eq!(e.value_size(), 1); + } + + #[test] + fn test_lagrange_0_interval() { + let e = lagrange::create::(ReferenceCellType::Interval, 0, Continuity::Discontinuous); + assert_eq!(e.value_size(), 1); + let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 4)); + let mut points = rlst_dynamic_array2!(f64, [4, 1]); + *points.get_mut([0, 0]).unwrap() = 0.0; + *points.get_mut([1, 0]).unwrap() = 0.2; + *points.get_mut([2, 0]).unwrap() = 0.4; + *points.get_mut([3, 0]).unwrap() = 1.0; + e.tabulate(&points, 0, &mut data); + + for pt in 0..4 { + assert_relative_eq!(*data.get([0, pt, 0, 0]).unwrap(), 1.0); + } + check_dofs(e); + } + + #[test] + fn test_lagrange_1_interval() { + let e = lagrange::create::(ReferenceCellType::Interval, 1, Continuity::Continuous); + assert_eq!(e.value_size(), 1); + let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 4)); + let mut points = rlst_dynamic_array2!(f64, [4, 1]); + *points.get_mut([0, 0]).unwrap() = 0.0; + *points.get_mut([1, 0]).unwrap() = 0.2; + *points.get_mut([2, 0]).unwrap() = 0.4; + *points.get_mut([3, 0]).unwrap() = 1.0; + e.tabulate(&points, 0, &mut data); + + for pt in 0..4 { + assert_relative_eq!( + *data.get([0, pt, 0, 0]).unwrap(), + 1.0 - *points.get([pt, 0]).unwrap() + ); + assert_relative_eq!( + *data.get([0, pt, 1, 0]).unwrap(), + *points.get([pt, 0]).unwrap() + ); + } + check_dofs(e); + } + + #[test] + fn test_lagrange_0_triangle() { + let e = lagrange::create::(ReferenceCellType::Triangle, 0, Continuity::Discontinuous); + assert_eq!(e.value_size(), 1); + let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 6)); + + let mut points = rlst_dynamic_array2!(f64, [6, 2]); + *points.get_mut([0, 0]).unwrap() = 0.0; + *points.get_mut([0, 1]).unwrap() = 0.0; + *points.get_mut([1, 0]).unwrap() = 1.0; + *points.get_mut([1, 1]).unwrap() = 0.0; + *points.get_mut([2, 0]).unwrap() = 0.0; + *points.get_mut([2, 1]).unwrap() = 1.0; + *points.get_mut([3, 0]).unwrap() = 0.5; + *points.get_mut([3, 1]).unwrap() = 0.0; + *points.get_mut([4, 0]).unwrap() = 0.0; + *points.get_mut([4, 1]).unwrap() = 0.5; + *points.get_mut([5, 0]).unwrap() = 0.5; + *points.get_mut([5, 1]).unwrap() = 0.5; + + e.tabulate(&points, 0, &mut data); + + for pt in 0..6 { + assert_relative_eq!(*data.get([0, pt, 0, 0]).unwrap(), 1.0); + } + check_dofs(e); + } + + #[test] + fn test_lagrange_1_triangle() { + let e = lagrange::create::(ReferenceCellType::Triangle, 1, Continuity::Continuous); + assert_eq!(e.value_size(), 1); + let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 6)); + let mut points = rlst_dynamic_array2!(f64, [6, 2]); + *points.get_mut([0, 0]).unwrap() = 0.0; + *points.get_mut([0, 1]).unwrap() = 0.0; + *points.get_mut([1, 0]).unwrap() = 1.0; + *points.get_mut([1, 1]).unwrap() = 0.0; + *points.get_mut([2, 0]).unwrap() = 0.0; + *points.get_mut([2, 1]).unwrap() = 1.0; + *points.get_mut([3, 0]).unwrap() = 0.5; + *points.get_mut([3, 1]).unwrap() = 0.0; + *points.get_mut([4, 0]).unwrap() = 0.0; + *points.get_mut([4, 1]).unwrap() = 0.5; + *points.get_mut([5, 0]).unwrap() = 0.5; + *points.get_mut([5, 1]).unwrap() = 0.5; + e.tabulate(&points, 0, &mut data); + + for pt in 0..6 { + assert_relative_eq!( + *data.get([0, pt, 0, 0]).unwrap(), + 1.0 - *points.get([pt, 0]).unwrap() - *points.get([pt, 1]).unwrap() + ); + assert_relative_eq!( + *data.get([0, pt, 1, 0]).unwrap(), + *points.get([pt, 0]).unwrap() + ); + assert_relative_eq!( + *data.get([0, pt, 2, 0]).unwrap(), + *points.get([pt, 1]).unwrap() + ); + } + check_dofs(e); + } + + #[test] + fn test_lagrange_higher_degree_triangle() { + lagrange::create::(ReferenceCellType::Triangle, 2, Continuity::Continuous); + lagrange::create::(ReferenceCellType::Triangle, 3, Continuity::Continuous); + lagrange::create::(ReferenceCellType::Triangle, 4, Continuity::Continuous); + lagrange::create::(ReferenceCellType::Triangle, 5, Continuity::Continuous); + + lagrange::create::(ReferenceCellType::Triangle, 2, Continuity::Discontinuous); + lagrange::create::(ReferenceCellType::Triangle, 3, Continuity::Discontinuous); + lagrange::create::(ReferenceCellType::Triangle, 4, Continuity::Discontinuous); + lagrange::create::(ReferenceCellType::Triangle, 5, Continuity::Discontinuous); + } + + #[test] + fn test_lagrange_higher_degree_interval() { + lagrange::create::(ReferenceCellType::Interval, 2, Continuity::Continuous); + lagrange::create::(ReferenceCellType::Interval, 3, Continuity::Continuous); + lagrange::create::(ReferenceCellType::Interval, 4, Continuity::Continuous); + lagrange::create::(ReferenceCellType::Interval, 5, Continuity::Continuous); + + lagrange::create::(ReferenceCellType::Interval, 2, Continuity::Discontinuous); + lagrange::create::(ReferenceCellType::Interval, 3, Continuity::Discontinuous); + lagrange::create::(ReferenceCellType::Interval, 4, Continuity::Discontinuous); + lagrange::create::(ReferenceCellType::Interval, 5, Continuity::Discontinuous); + } + + #[test] + fn test_lagrange_higher_degree_quadrilateral() { + lagrange::create::(ReferenceCellType::Quadrilateral, 2, Continuity::Continuous); + lagrange::create::(ReferenceCellType::Quadrilateral, 3, Continuity::Continuous); + lagrange::create::(ReferenceCellType::Quadrilateral, 4, Continuity::Continuous); + lagrange::create::(ReferenceCellType::Quadrilateral, 5, Continuity::Continuous); + + lagrange::create::( + ReferenceCellType::Quadrilateral, + 2, + Continuity::Discontinuous, + ); + lagrange::create::( + ReferenceCellType::Quadrilateral, + 3, + Continuity::Discontinuous, + ); + lagrange::create::( + ReferenceCellType::Quadrilateral, + 4, + Continuity::Discontinuous, + ); + lagrange::create::( + ReferenceCellType::Quadrilateral, + 5, + Continuity::Discontinuous, + ); + } + + #[test] + fn test_lagrange_0_quadrilateral() { + let e = lagrange::create::( + ReferenceCellType::Quadrilateral, + 0, + Continuity::Discontinuous, + ); + assert_eq!(e.value_size(), 1); + let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 6)); + let mut points = rlst_dynamic_array2!(f64, [6, 2]); + *points.get_mut([0, 0]).unwrap() = 0.0; + *points.get_mut([0, 1]).unwrap() = 0.0; + *points.get_mut([1, 0]).unwrap() = 1.0; + *points.get_mut([1, 1]).unwrap() = 0.0; + *points.get_mut([2, 0]).unwrap() = 0.0; + *points.get_mut([2, 1]).unwrap() = 1.0; + *points.get_mut([3, 0]).unwrap() = 0.5; + *points.get_mut([3, 1]).unwrap() = 0.0; + *points.get_mut([4, 0]).unwrap() = 0.0; + *points.get_mut([4, 1]).unwrap() = 0.5; + *points.get_mut([5, 0]).unwrap() = 0.5; + *points.get_mut([5, 1]).unwrap() = 0.5; + e.tabulate(&points, 0, &mut data); + + for pt in 0..6 { + assert_relative_eq!(*data.get([0, pt, 0, 0]).unwrap(), 1.0); + } + check_dofs(e); + } + + #[test] + fn test_lagrange_1_quadrilateral() { + let e = + lagrange::create::(ReferenceCellType::Quadrilateral, 1, Continuity::Continuous); + assert_eq!(e.value_size(), 1); + let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 6)); + let mut points = rlst_dynamic_array2!(f64, [6, 2]); + *points.get_mut([0, 0]).unwrap() = 0.0; + *points.get_mut([0, 1]).unwrap() = 0.0; + *points.get_mut([1, 0]).unwrap() = 1.0; + *points.get_mut([1, 1]).unwrap() = 0.0; + *points.get_mut([2, 0]).unwrap() = 0.0; + *points.get_mut([2, 1]).unwrap() = 1.0; + *points.get_mut([3, 0]).unwrap() = 1.0; + *points.get_mut([3, 1]).unwrap() = 1.0; + *points.get_mut([4, 0]).unwrap() = 0.25; + *points.get_mut([4, 1]).unwrap() = 0.5; + *points.get_mut([5, 0]).unwrap() = 0.3; + *points.get_mut([5, 1]).unwrap() = 0.2; + + e.tabulate(&points, 0, &mut data); + + for pt in 0..6 { + assert_relative_eq!( + *data.get([0, pt, 0, 0]).unwrap(), + (1.0 - *points.get([pt, 0]).unwrap()) * (1.0 - *points.get([pt, 1]).unwrap()) + ); + assert_relative_eq!( + *data.get([0, pt, 1, 0]).unwrap(), + *points.get([pt, 0]).unwrap() * (1.0 - *points.get([pt, 1]).unwrap()) + ); + assert_relative_eq!( + *data.get([0, pt, 2, 0]).unwrap(), + (1.0 - *points.get([pt, 0]).unwrap()) * *points.get([pt, 1]).unwrap() + ); + assert_relative_eq!( + *data.get([0, pt, 3, 0]).unwrap(), + *points.get([pt, 0]).unwrap() * *points.get([pt, 1]).unwrap() + ); + } + check_dofs(e); + } + + #[test] + fn test_lagrange_2_quadrilateral() { + let e = + lagrange::create::(ReferenceCellType::Quadrilateral, 2, Continuity::Continuous); + assert_eq!(e.value_size(), 1); + let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 6)); + let mut points = rlst_dynamic_array2!(f64, [6, 2]); + *points.get_mut([0, 0]).unwrap() = 0.0; + *points.get_mut([0, 1]).unwrap() = 0.0; + *points.get_mut([1, 0]).unwrap() = 1.0; + *points.get_mut([1, 1]).unwrap() = 0.0; + *points.get_mut([2, 0]).unwrap() = 0.0; + *points.get_mut([2, 1]).unwrap() = 1.0; + *points.get_mut([3, 0]).unwrap() = 1.0; + *points.get_mut([3, 1]).unwrap() = 1.0; + *points.get_mut([4, 0]).unwrap() = 0.25; + *points.get_mut([4, 1]).unwrap() = 0.5; + *points.get_mut([5, 0]).unwrap() = 0.3; + *points.get_mut([5, 1]).unwrap() = 0.2; + e.tabulate(&points, 0, &mut data); + + for pt in 0..6 { + let x = *points.get([pt, 0]).unwrap(); + let y = *points.get([pt, 1]).unwrap(); + assert_relative_eq!( + *data.get([0, pt, 0, 0]).unwrap(), + (1.0 - x) * (1.0 - 2.0 * x) * (1.0 - y) * (1.0 - 2.0 * y), + epsilon = 1e-14 + ); + assert_relative_eq!( + *data.get([0, pt, 1, 0]).unwrap(), + x * (2.0 * x - 1.0) * (1.0 - y) * (1.0 - 2.0 * y), + epsilon = 1e-14 + ); + assert_relative_eq!( + *data.get([0, pt, 2, 0]).unwrap(), + (1.0 - x) * (1.0 - 2.0 * x) * y * (2.0 * y - 1.0), + epsilon = 1e-14 + ); + assert_relative_eq!( + *data.get([0, pt, 3, 0]).unwrap(), + x * (2.0 * x - 1.0) * y * (2.0 * y - 1.0), + epsilon = 1e-14 + ); + assert_relative_eq!( + *data.get([0, pt, 4, 0]).unwrap(), + 4.0 * x * (1.0 - x) * (1.0 - y) * (1.0 - 2.0 * y), + epsilon = 1e-14 + ); + assert_relative_eq!( + *data.get([0, pt, 5, 0]).unwrap(), + (1.0 - x) * (1.0 - 2.0 * x) * 4.0 * y * (1.0 - y), + epsilon = 1e-14 + ); + assert_relative_eq!( + *data.get([0, pt, 6, 0]).unwrap(), + x * (2.0 * x - 1.0) * 4.0 * y * (1.0 - y), + epsilon = 1e-14 + ); + assert_relative_eq!( + *data.get([0, pt, 7, 0]).unwrap(), + 4.0 * x * (1.0 - x) * y * (2.0 * y - 1.0), + epsilon = 1e-14 + ); + assert_relative_eq!( + *data.get([0, pt, 8, 0]).unwrap(), + 4.0 * x * (1.0 - x) * 4.0 * y * (1.0 - y), + epsilon = 1e-14 + ); + } + check_dofs(e); + } + + #[test] + fn test_raviart_thomas_1_triangle() { + let e = raviart_thomas::create(ReferenceCellType::Triangle, 1, Continuity::Continuous); + assert_eq!(e.value_size(), 2); + let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 6)); + let mut points = rlst_dynamic_array2!(f64, [6, 2]); + *points.get_mut([0, 0]).unwrap() = 0.0; + *points.get_mut([0, 1]).unwrap() = 0.0; + *points.get_mut([1, 0]).unwrap() = 1.0; + *points.get_mut([1, 1]).unwrap() = 0.0; + *points.get_mut([2, 0]).unwrap() = 0.0; + *points.get_mut([2, 1]).unwrap() = 1.0; + *points.get_mut([3, 0]).unwrap() = 0.5; + *points.get_mut([3, 1]).unwrap() = 0.0; + *points.get_mut([4, 0]).unwrap() = 0.0; + *points.get_mut([4, 1]).unwrap() = 0.5; + *points.get_mut([5, 0]).unwrap() = 0.5; + *points.get_mut([5, 1]).unwrap() = 0.5; + e.tabulate(&points, 0, &mut data); + + for pt in 0..6 { + assert_relative_eq!( + *data.get([0, pt, 0, 0]).unwrap(), + -*points.get([pt, 0]).unwrap() + ); + assert_relative_eq!( + *data.get([0, pt, 0, 1]).unwrap(), + -*points.get([pt, 1]).unwrap() + ); + assert_relative_eq!( + *data.get([0, pt, 1, 0]).unwrap(), + *points.get([pt, 0]).unwrap() - 1.0 + ); + assert_relative_eq!( + *data.get([0, pt, 1, 1]).unwrap(), + *points.get([pt, 1]).unwrap() + ); + assert_relative_eq!( + *data.get([0, pt, 2, 0]).unwrap(), + -*points.get([pt, 0]).unwrap() + ); + assert_relative_eq!( + *data.get([0, pt, 2, 1]).unwrap(), + 1.0 - *points.get([pt, 1]).unwrap() + ); + } + check_dofs(e); + } +} diff --git a/src/ciarlet/lagrange.rs b/src/ciarlet/lagrange.rs new file mode 100644 index 0000000..b0e7b78 --- /dev/null +++ b/src/ciarlet/lagrange.rs @@ -0,0 +1,214 @@ +//! Lagrange elements + +use super::CiarletElement; +use crate::polynomials::polynomial_count; +use crate::reference_cell; +use crate::traits::ElementFamily; +use crate::types::{Continuity, MapType, ReferenceCellType}; +use rlst::{ + dense::array::views::ArrayViewMut, rlst_dynamic_array2, rlst_dynamic_array3, Array, BaseArray, + MatrixInverse, RandomAccessMut, RlstScalar, VectorContainer, +}; +use std::marker::PhantomData; + +/// Create a Lagrange element +pub fn create( + cell_type: ReferenceCellType, + degree: usize, + continuity: Continuity, +) -> CiarletElement +where + for<'a> Array, 2>, 2>, 2>: MatrixInverse, +{ + let dim = polynomial_count(cell_type, degree); + let tdim = reference_cell::dim(cell_type); + let mut wcoeffs = rlst_dynamic_array3!(T, [dim, 1, dim]); + for i in 0..dim { + *wcoeffs.get_mut([i, 0, i]).unwrap() = T::from(1.0).unwrap(); + } + + let mut x = [vec![], vec![], vec![], vec![]]; + let mut m = [vec![], vec![], vec![], vec![]]; + let entity_counts = reference_cell::entity_counts(cell_type); + let vertices = reference_cell::vertices::(cell_type); + if degree == 0 { + if continuity == Continuity::Continuous { + panic!("Cannot create continuous degree 0 Lagrange element"); + } + for (d, counts) in entity_counts.iter().enumerate() { + for _e in 0..*counts { + x[d].push(rlst_dynamic_array2!(T::Real, [0, tdim])); + m[d].push(rlst_dynamic_array3!(T, [0, 1, 0])); + } + } + let mut midp = rlst_dynamic_array2!(T::Real, [1, tdim]); + let nvertices = entity_counts[0]; + for i in 0..tdim { + for vertex in &vertices { + *midp.get_mut([0, i]).unwrap() += num::cast::<_, T::Real>(vertex[i]).unwrap(); + } + *midp.get_mut([0, i]).unwrap() /= num::cast::<_, T::Real>(nvertices).unwrap(); + } + x[tdim].push(midp); + let mut mentry = rlst_dynamic_array3!(T, [1, 1, 1]); + *mentry.get_mut([0, 0, 0]).unwrap() = T::from(1.0).unwrap(); + m[tdim].push(mentry); + } else { + let edges = reference_cell::edges(cell_type); + let faces = reference_cell::faces(cell_type); + // TODO: GLL points + for vertex in &vertices { + let mut pts = rlst_dynamic_array2!(T::Real, [1, tdim]); + for (i, v) in vertex.iter().enumerate() { + *pts.get_mut([0, i]).unwrap() = num::cast::<_, T::Real>(*v).unwrap(); + } + x[0].push(pts); + let mut mentry = rlst_dynamic_array3!(T, [1, 1, 1]); + *mentry.get_mut([0, 0, 0]).unwrap() = T::from(1.0).unwrap(); + m[0].push(mentry); + } + for e in &edges { + let mut pts = rlst_dynamic_array2!(T::Real, [degree - 1, tdim]); + let [vn0, vn1] = e[..] else { + panic!(); + }; + let v0 = &vertices[vn0]; + let v1 = &vertices[vn1]; + let mut ident = rlst_dynamic_array3!(T, [degree - 1, 1, degree - 1]); + + for i in 1..degree { + *ident.get_mut([i - 1, 0, i - 1]).unwrap() = T::from(1.0).unwrap(); + for j in 0..tdim { + *pts.get_mut([i - 1, j]).unwrap() = num::cast::<_, T::Real>(v0[j]).unwrap() + + num::cast::<_, T::Real>(i).unwrap() + / num::cast::<_, T::Real>(degree).unwrap() + * num::cast::<_, T::Real>(v1[j] - v0[j]).unwrap(); + } + } + x[1].push(pts); + m[1].push(ident); + } + for (e, face_type) in reference_cell::entity_types(cell_type)[2] + .iter() + .enumerate() + { + let npts = match face_type { + ReferenceCellType::Triangle => { + if degree > 2 { + (degree - 1) * (degree - 2) / 2 + } else { + 0 + } + } + ReferenceCellType::Quadrilateral => (degree - 1).pow(2), + _ => { + panic!("Unsupported face type"); + } + }; + let mut pts = rlst_dynamic_array2!(T::Real, [npts, tdim]); + + let [vn0, vn1, vn2] = faces[e][..3] else { + panic!(); + }; + let v0 = &vertices[vn0]; + let v1 = &vertices[vn1]; + let v2 = &vertices[vn2]; + + match face_type { + ReferenceCellType::Triangle => { + let mut n = 0; + for i0 in 1..degree { + for i1 in 1..degree - i0 { + for j in 0..tdim { + *pts.get_mut([n, j]).unwrap() = num::cast::<_, T::Real>(v0[j]) + .unwrap() + + num::cast::<_, T::Real>(i0).unwrap() + / num::cast::<_, T::Real>(degree).unwrap() + * num::cast::<_, T::Real>(v1[j] - v0[j]).unwrap() + + num::cast::<_, T::Real>(i1).unwrap() + / num::cast::<_, T::Real>(degree).unwrap() + * num::cast::<_, T::Real>(v2[j] - v0[j]).unwrap(); + } + n += 1; + } + } + } + ReferenceCellType::Quadrilateral => { + let mut n = 0; + for i0 in 1..degree { + for i1 in 1..degree { + for j in 0..tdim { + *pts.get_mut([n, j]).unwrap() = num::cast::<_, T::Real>(v0[j]) + .unwrap() + + num::cast::<_, T::Real>(i0).unwrap() + / num::cast::<_, T::Real>(degree).unwrap() + * num::cast::<_, T::Real>(v1[j] - v0[j]).unwrap() + + num::cast::<_, T::Real>(i1).unwrap() + / num::cast::<_, T::Real>(degree).unwrap() + * num::cast::<_, T::Real>(v2[j] - v0[j]).unwrap(); + } + n += 1; + } + } + } + _ => { + panic!("Unsupported face type."); + } + }; + + let mut ident = rlst_dynamic_array3!(T, [npts, 1, npts]); + for i in 0..npts { + *ident.get_mut([i, 0, i]).unwrap() = T::from(1.0).unwrap(); + } + x[2].push(pts); + m[2].push(ident); + } + } + CiarletElement::::create( + cell_type, + degree, + vec![], + wcoeffs, + x, + m, + MapType::Identity, + continuity, + degree, + ) +} + +/// Lagrange element family +pub struct LagrangeElementFamily +where + for<'a> Array, 2>, 2>, 2>: MatrixInverse, +{ + degree: usize, + continuity: Continuity, + _t: PhantomData, +} + +impl LagrangeElementFamily +where + for<'a> Array, 2>, 2>, 2>: MatrixInverse, +{ + /// Create new family + pub fn new(degree: usize, continuity: Continuity) -> Self { + Self { + degree, + continuity, + _t: PhantomData, + } + } +} + +impl ElementFamily for LagrangeElementFamily +where + for<'a> Array, 2>, 2>, 2>: MatrixInverse, +{ + type T = T; + type FiniteElement = CiarletElement; + type CellType = ReferenceCellType; + fn element(&self, cell_type: ReferenceCellType) -> CiarletElement { + create::(cell_type, self.degree, self.continuity) + } +} diff --git a/src/ciarlet/raviart_thomas.rs b/src/ciarlet/raviart_thomas.rs new file mode 100644 index 0000000..eb391b0 --- /dev/null +++ b/src/ciarlet/raviart_thomas.rs @@ -0,0 +1,135 @@ +//! Raviart-Thomas elements + +use super::CiarletElement; +use crate::polynomials::polynomial_count; +use crate::reference_cell; +use crate::traits::ElementFamily; +use crate::types::{Continuity, MapType, ReferenceCellType}; +use rlst::MatrixInverse; +use rlst::RlstScalar; +use rlst::{ + dense::array::views::ArrayViewMut, rlst_dynamic_array2, rlst_dynamic_array3, Array, BaseArray, + RandomAccessMut, VectorContainer, +}; +use std::marker::PhantomData; + +/// Create a Raviart-Thomas element +pub fn create( + cell_type: ReferenceCellType, + degree: usize, + continuity: Continuity, +) -> CiarletElement +where + for<'a> Array, 2>, 2>, 2>: MatrixInverse, +{ + if cell_type != ReferenceCellType::Triangle && cell_type != ReferenceCellType::Quadrilateral { + panic!("Unsupported cell type"); + } + + if cell_type != ReferenceCellType::Triangle { + panic!("RT elements on quadrilaterals not implemented yet"); + } + if degree != 1 { + panic!("Degree > 1 RT elements not implemented yet"); + } + + let pdim = polynomial_count(cell_type, degree); + let tdim = reference_cell::dim(cell_type); + let edim = tdim * polynomial_count(cell_type, degree - 1) + degree; + + let mut wcoeffs = rlst_dynamic_array3!(T, [edim, tdim, pdim]); + + // [sqrt(2), 6*y - 2, 4*sqrt(3)*(x + y/2 - 1/2)] + + // norm(x**2 + y**2) + // sqrt(70)/30 + + *wcoeffs.get_mut([0, 0, 0]).unwrap() = T::from(1.0).unwrap(); + *wcoeffs.get_mut([1, 1, 0]).unwrap() = T::from(1.0).unwrap(); + *wcoeffs.get_mut([2, 0, 1]).unwrap() = T::from(-0.5).unwrap() / T::sqrt(T::from(2.0).unwrap()); + *wcoeffs.get_mut([2, 0, 2]).unwrap() = T::from(0.5).unwrap() * T::sqrt(T::from(1.5).unwrap()); + *wcoeffs.get_mut([2, 1, 1]).unwrap() = T::from(1.0).unwrap() / T::sqrt(T::from(2.0).unwrap()); + + let mut x = [vec![], vec![], vec![], vec![]]; + let mut m = [vec![], vec![], vec![], vec![]]; + + let entity_counts = reference_cell::entity_counts(cell_type); + let vertices = reference_cell::vertices::(cell_type); + let edges = reference_cell::edges(cell_type); + + for _e in 0..entity_counts[0] { + x[0].push(rlst_dynamic_array2!(T::Real, [0, tdim])); + m[0].push(rlst_dynamic_array3!(T, [0, 2, 0])); + } + + for e in &edges { + let mut pts = rlst_dynamic_array2!(T::Real, [1, tdim]); + let mut mat = rlst_dynamic_array3!(T, [1, 2, 1]); + let [vn0, vn1] = e[..] else { + panic!(); + }; + let v0 = &vertices[vn0]; + let v1 = &vertices[vn1]; + for i in 0..tdim { + *pts.get_mut([0, i]).unwrap() = num::cast::<_, T::Real>(v0[i] + v1[i]).unwrap() + / num::cast::<_, T::Real>(2.0).unwrap(); + } + *mat.get_mut([0, 0, 0]).unwrap() = T::from(v0[1] - v1[1]).unwrap(); + *mat.get_mut([0, 1, 0]).unwrap() = T::from(v1[0] - v0[0]).unwrap(); + x[1].push(pts); + m[1].push(mat); + } + + for _e in 0..entity_counts[2] { + x[2].push(rlst_dynamic_array2!(T::Real, [0, tdim])); + m[2].push(rlst_dynamic_array3!(T, [0, 2, 0])) + } + + CiarletElement::create( + cell_type, + degree, + vec![2], + wcoeffs, + x, + m, + MapType::ContravariantPiola, + continuity, + degree, + ) +} + +/// Raviart-Thomas element family +pub struct RaviartThomasElementFamily +where + for<'a> Array, 2>, 2>, 2>: MatrixInverse, +{ + degree: usize, + continuity: Continuity, + _t: PhantomData, +} + +impl RaviartThomasElementFamily +where + for<'a> Array, 2>, 2>, 2>: MatrixInverse, +{ + /// Create new family + pub fn new(degree: usize, continuity: Continuity) -> Self { + Self { + degree, + continuity, + _t: PhantomData, + } + } +} + +impl ElementFamily for RaviartThomasElementFamily +where + for<'a> Array, 2>, 2>, 2>: MatrixInverse, +{ + type T = T; + type CellType = ReferenceCellType; + type FiniteElement = CiarletElement; + fn element(&self, cell_type: ReferenceCellType) -> CiarletElement { + create::(cell_type, self.degree, self.continuity) + } +} diff --git a/src/lib.rs b/src/lib.rs index 04bebcf..0b206db 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,3 +1,15 @@ //! n-dimensional grid #![cfg_attr(feature = "strict", deny(warnings))] #![warn(missing_docs)] + +pub mod ciarlet; +pub mod polynomials; +pub mod reference_cell; +pub mod traits; +pub mod types; + +#[cfg(test)] +mod test { + extern crate blas_src; + extern crate lapack_src; +} diff --git a/src/polynomials.rs b/src/polynomials.rs new file mode 100644 index 0000000..b22abc1 --- /dev/null +++ b/src/polynomials.rs @@ -0,0 +1,981 @@ +//! Orthonormal polynomials + +use crate::types::ReferenceCellType; +use rlst::RlstScalar; +use rlst::{RandomAccessByRef, RandomAccessMut, Shape}; +/// Tabulate orthonormal polynomials on a interval +fn tabulate_legendre_polynomials_interval< + T: RlstScalar, + Array2: RandomAccessByRef<2, Item = T::Real> + Shape<2>, + Array3Mut: RandomAccessMut<3, Item = T> + RandomAccessByRef<3, Item = T> + Shape<3>, +>( + points: &Array2, + degree: usize, + derivatives: usize, + data: &mut Array3Mut, +) { + assert_eq!(data.shape()[0], derivatives + 1); + assert_eq!(data.shape()[1], degree + 1); + assert_eq!(data.shape()[2], points.shape()[0]); + assert_eq!(points.shape()[1], 1); + + for i in 0..data.shape()[2] { + *data.get_mut([0, 0, i]).unwrap() = T::from(1.0).unwrap(); + } + for k in 1..data.shape()[0] { + for i in 0..data.shape()[2] { + *data.get_mut([k, 0, i]).unwrap() = T::from(0.0).unwrap(); + } + } + + for k in 0..derivatives + 1 { + for p in 1..degree + 1 { + let a = T::from(1.0).unwrap() - T::from(1.0).unwrap() / T::from(p).unwrap(); + let b = (a + T::from(1.0).unwrap()) + * ((T::from(2.0).unwrap() * T::from(p).unwrap() + T::from(1.0).unwrap()) + / (T::from(2.0).unwrap() * T::from(p).unwrap() - T::from(1.0).unwrap())) + .sqrt(); + for i in 0..data.shape()[2] { + let d = *data.get([k, p - 1, i]).unwrap(); + *data.get_mut([k, p, i]).unwrap() = + (T::from(*points.get([i, 0]).unwrap()).unwrap() * T::from(2.0).unwrap() + - T::from(1.0).unwrap()) + * d + * b; + } + if p > 1 { + let c = a + * ((T::from(2.0).unwrap() * T::from(p).unwrap() + T::from(1.0).unwrap()) + / (T::from(2.0).unwrap() * T::from(p).unwrap() - T::from(3.0).unwrap())) + .sqrt(); + for i in 0..data.shape()[2] { + let d = *data.get([k, p - 2, i]).unwrap(); + *data.get_mut([k, p, i]).unwrap() -= d * c; + } + } + if k > 0 { + for i in 0..data.shape()[2] { + let d = *data.get([k - 1, p - 1, i]).unwrap(); + *data.get_mut([k, p, i]).unwrap() += + T::from(2.0).unwrap() * T::from(k).unwrap() * d * b; + } + } + } + } +} + +fn tri_index(i: usize, j: usize) -> usize { + (i + j + 1) * (i + j) / 2 + j +} + +fn quad_index(i: usize, j: usize, n: usize) -> usize { + j * (n + 1) + i +} + +/// Tabulate orthonormal polynomials on a quadrilateral +fn tabulate_legendre_polynomials_quadrilateral< + T: RlstScalar, + Array2: RandomAccessByRef<2, Item = T::Real> + Shape<2>, + Array3Mut: RandomAccessMut<3, Item = T> + RandomAccessByRef<3, Item = T> + Shape<3>, +>( + points: &Array2, + degree: usize, + derivatives: usize, + data: &mut Array3Mut, +) { + assert_eq!(data.shape()[0], (derivatives + 1) * (derivatives + 2) / 2); + assert_eq!(data.shape()[1], (degree + 1) * (degree + 1)); + assert_eq!(data.shape()[2], points.shape()[0]); + assert_eq!(points.shape()[1], 2); + + for i in 0..data.shape()[2] { + *data + .get_mut([tri_index(0, 0), quad_index(0, 0, degree), i]) + .unwrap() = T::from(1.0).unwrap(); + } + + // Tabulate polynomials in x + for k in 1..derivatives + 1 { + for i in 0..data.shape()[2] { + *data + .get_mut([tri_index(k, 0), quad_index(0, 0, degree), i]) + .unwrap() = T::from(0.0).unwrap(); + } + } + + for k in 0..derivatives + 1 { + for p in 1..degree + 1 { + let a = T::from(1.0).unwrap() - T::from(1.0).unwrap() / T::from(p).unwrap(); + let b = (a + T::from(1.0).unwrap()) + * ((T::from(2.0).unwrap() * T::from(p).unwrap() + T::from(1.0).unwrap()) + / (T::from(2.0).unwrap() * T::from(p).unwrap() - T::from(1.0).unwrap())) + .sqrt(); + for i in 0..data.shape()[2] { + let d = *data + .get([tri_index(k, 0), quad_index(p - 1, 0, degree), i]) + .unwrap(); + *data + .get_mut([tri_index(k, 0), quad_index(p, 0, degree), i]) + .unwrap() = (T::from(*points.get([i, 0]).unwrap()).unwrap() + * T::from(2.0).unwrap() + - T::from(1.0).unwrap()) + * d + * b; + } + if p > 1 { + let c = a + * ((T::from(2.0).unwrap() * T::from(p).unwrap() + T::from(1.0).unwrap()) + / (T::from(2.0).unwrap() * T::from(p).unwrap() - T::from(3.0).unwrap())) + .sqrt(); + for i in 0..data.shape()[2] { + let d = *data + .get([tri_index(k, 0), quad_index(p - 2, 0, degree), i]) + .unwrap(); + *data + .get_mut([tri_index(k, 0), quad_index(p, 0, degree), i]) + .unwrap() -= d * c; + } + } + if k > 0 { + for i in 0..data.shape()[2] { + let d = *data + .get([tri_index(k - 1, 0), quad_index(p - 1, 0, degree), i]) + .unwrap(); + *data + .get_mut([tri_index(k, 0), quad_index(p, 0, degree), i]) + .unwrap() += T::from(2.0).unwrap() * T::from(k).unwrap() * d * b; + } + } + } + } + + // Tabulate polynomials in y + for k in 1..derivatives + 1 { + for i in 0..data.shape()[2] { + *data + .get_mut([tri_index(0, k), quad_index(0, 0, degree), i]) + .unwrap() = T::from(0.0).unwrap(); + } + } + + for k in 0..derivatives + 1 { + for p in 1..degree + 1 { + let a = T::from(1.0).unwrap() - T::from(1.0).unwrap() / T::from(p).unwrap(); + let b = (a + T::from(1.0).unwrap()) + * ((T::from(2.0).unwrap() * T::from(p).unwrap() + T::from(1.0).unwrap()) + / (T::from(2.0).unwrap() * T::from(p).unwrap() - T::from(1.0).unwrap())) + .sqrt(); + for i in 0..data.shape()[2] { + let d = *data + .get([tri_index(0, k), quad_index(0, p - 1, degree), i]) + .unwrap(); + *data + .get_mut([tri_index(0, k), quad_index(0, p, degree), i]) + .unwrap() = (T::from(*points.get([i, 1]).unwrap()).unwrap() + * T::from(2.0).unwrap() + - T::from(1.0).unwrap()) + * d + * b; + } + if p > 1 { + let c = a + * ((T::from(2.0).unwrap() * T::from(p).unwrap() + T::from(1.0).unwrap()) + / (T::from(2.0).unwrap() * T::from(p).unwrap() - T::from(3.0).unwrap())) + .sqrt(); + for i in 0..data.shape()[2] { + let d = *data + .get([tri_index(0, k), quad_index(0, p - 2, degree), i]) + .unwrap(); + *data + .get_mut([tri_index(0, k), quad_index(0, p, degree), i]) + .unwrap() -= d * c; + } + } + if k > 0 { + for i in 0..data.shape()[2] { + let d = *data + .get([tri_index(0, k - 1), quad_index(0, p - 1, degree), i]) + .unwrap(); + *data + .get_mut([tri_index(0, k), quad_index(0, p, degree), i]) + .unwrap() += T::from(2.0).unwrap() * T::from(k).unwrap() * d * b; + } + } + } + } + + // Fill in the rest of the values as products + for kx in 0..derivatives + 1 { + for ky in 0..derivatives + 1 - kx { + for px in 1..degree + 1 { + for py in 1..degree + 1 { + for i in 0..data.shape()[2] { + let d = *data + .get([tri_index(0, ky), quad_index(0, py, degree), i]) + .unwrap(); + *data + .get_mut([tri_index(kx, ky), quad_index(px, py, degree), i]) + .unwrap() = *data + .get([tri_index(kx, 0), quad_index(px, 0, degree), i]) + .unwrap() + * d; + } + } + } + } + } +} +/// Tabulate orthonormal polynomials on a triangle +fn tabulate_legendre_polynomials_triangle< + T: RlstScalar, + Array2: RandomAccessByRef<2, Item = T::Real> + Shape<2>, + Array3Mut: RandomAccessMut<3, Item = T> + RandomAccessByRef<3, Item = T> + Shape<3>, +>( + points: &Array2, + degree: usize, + derivatives: usize, + data: &mut Array3Mut, +) { + assert_eq!(data.shape()[0], (derivatives + 1) * (derivatives + 2) / 2); + assert_eq!(data.shape()[1], (degree + 1) * (degree + 2) / 2); + assert_eq!(data.shape()[2], points.shape()[0]); + assert_eq!(points.shape()[1], 2); + + for i in 0..data.shape()[2] { + *data.get_mut([tri_index(0, 0), tri_index(0, 0), i]).unwrap() = + T::sqrt(T::from(2.0).unwrap()); + } + + for k in 1..data.shape()[0] { + for i in 0..data.shape()[2] { + *data.get_mut([k, tri_index(0, 0), i]).unwrap() = T::from(0.0).unwrap(); + } + } + + for kx in 0..derivatives + 1 { + for ky in 0..derivatives + 1 - kx { + for p in 1..degree + 1 { + let a = T::from(2.0).unwrap() - T::from(1.0).unwrap() / T::from(p).unwrap(); + let scale1 = T::sqrt( + (T::from(p).unwrap() + T::from(0.5).unwrap()) + * (T::from(p).unwrap() + T::from(1.0).unwrap()) + / ((T::from(p).unwrap() - T::from(0.5).unwrap()) * T::from(p).unwrap()), + ); + for i in 0..data.shape()[2] { + let d = *data + .get([tri_index(kx, ky), tri_index(0, p - 1), i]) + .unwrap(); + *data + .get_mut([tri_index(kx, ky), tri_index(0, p), i]) + .unwrap() = (T::from(*points.get([i, 0]).unwrap()).unwrap() + * T::from(2.0).unwrap() + + T::from(*points.get([i, 1]).unwrap()).unwrap() + - T::from(1.0).unwrap()) + * d + * a + * scale1; + } + if kx > 0 { + for i in 0..data.shape()[2] { + let d = *data + .get([tri_index(kx - 1, ky), tri_index(0, p - 1), i]) + .unwrap(); + *data + .get_mut([tri_index(kx, ky), tri_index(0, p), i]) + .unwrap() += + T::from(2.0).unwrap() * T::from(kx).unwrap() * a * d * scale1; + } + } + if ky > 0 { + for i in 0..data.shape()[2] { + let d = *data + .get([tri_index(kx, ky - 1), tri_index(0, p - 1), i]) + .unwrap(); + *data + .get_mut([tri_index(kx, ky), tri_index(0, p), i]) + .unwrap() += T::from(ky).unwrap() * a * d * scale1; + } + } + if p > 1 { + let scale2 = T::sqrt( + (T::from(p).unwrap() + T::from(0.5).unwrap()) + * (T::from(p).unwrap() + T::from(1.0).unwrap()), + ) / T::sqrt( + (T::from(p).unwrap() - T::from(1.5).unwrap()) + * (T::from(p).unwrap() - T::from(1.0).unwrap()), + ); + + for i in 0..data.shape()[2] { + let b = + T::from(1.0).unwrap() - T::from(*points.get([i, 1]).unwrap()).unwrap(); + let d = *data + .get([tri_index(kx, ky), tri_index(0, p - 2), i]) + .unwrap(); + *data + .get_mut([tri_index(kx, ky), tri_index(0, p), i]) + .unwrap() -= b * b * d * (a - T::from(1.0).unwrap()) * scale2; + } + if ky > 0 { + for i in 0..data.shape()[2] { + let d = *data + .get([tri_index(kx, ky - 1), tri_index(0, p - 2), i]) + .unwrap(); + *data + .get_mut([tri_index(kx, ky), tri_index(0, p), i]) + .unwrap() -= T::from(2.0).unwrap() + * T::from(ky).unwrap() + * (T::from(*points.get([i, 1]).unwrap()).unwrap() + - T::from(1.0).unwrap()) + * d + * scale2 + * (a - T::from(1.0).unwrap()); + } + } + if ky > 1 { + for i in 0..data.shape()[2] { + let d = *data + .get([tri_index(kx, ky - 2), tri_index(0, p - 2), i]) + .unwrap(); + *data + .get_mut([tri_index(kx, ky), tri_index(0, p), i]) + .unwrap() -= T::from(ky).unwrap() + * (T::from(ky).unwrap() - T::from(1.0).unwrap()) + * d + * scale2 + * (a - T::from(1.0).unwrap()); + } + } + } + } + for p in 0..degree { + let scale3 = T::sqrt( + (T::from(p).unwrap() + T::from(2.0).unwrap()) + / (T::from(p).unwrap() + T::from(1.0).unwrap()), + ); + for i in 0..data.shape()[2] { + *data + .get_mut([tri_index(kx, ky), tri_index(1, p), i]) + .unwrap() = *data.get([tri_index(kx, ky), tri_index(0, p), i]).unwrap() + * scale3 + * ((T::from(*points.get([i, 1]).unwrap()).unwrap() + * T::from(2.0).unwrap() + - T::from(1.0).unwrap()) + * (T::from(1.5).unwrap() + T::from(p).unwrap()) + + T::from(0.5).unwrap() + + T::from(p).unwrap()); + } + if ky > 0 { + for i in 0..data.shape()[2] { + let d = *data + .get([tri_index(kx, ky - 1), tri_index(0, p), i]) + .unwrap(); + *data + .get_mut([tri_index(kx, ky), tri_index(1, p), i]) + .unwrap() += T::from(2.0).unwrap() + * T::from(ky).unwrap() + * (T::from(1.5).unwrap() + T::from(p).unwrap()) + * d + * scale3; + } + } + for q in 1..degree - p { + let scale4 = T::sqrt( + (T::from(p).unwrap() + T::from(q).unwrap() + T::from(2.0).unwrap()) + / (T::from(p).unwrap() + T::from(q).unwrap() + T::from(1.0).unwrap()), + ); + let scale5 = T::sqrt( + (T::from(p).unwrap() + T::from(q).unwrap() + T::from(2.0).unwrap()) + / (T::from(p).unwrap() + T::from(q).unwrap()), + ); + let a1 = T::from((p + q + 1) * (2 * p + 2 * q + 3)).unwrap() + / T::from((q + 1) * (2 * p + q + 2)).unwrap(); + let a2 = T::from((2 * p + 1) * (2 * p + 1) * (p + q + 1)).unwrap() + / T::from((q + 1) * (2 * p + q + 2) * (2 * p + 2 * q + 1)).unwrap(); + let a3 = T::from(q * (2 * p + q + 1) * (2 * p + 2 * q + 3)).unwrap() + / T::from((q + 1) * (2 * p + q + 2) * (2 * p + 2 * q + 1)).unwrap(); + + for i in 0..data.shape()[2] { + let d = *data.get([tri_index(kx, ky), tri_index(q, p), i]).unwrap(); + *data + .get_mut([tri_index(kx, ky), tri_index(q + 1, p), i]) + .unwrap() = d + * scale4 + * ((T::from(*points.get([i, 1]).unwrap()).unwrap() + * T::from(T::from(2.0).unwrap()).unwrap() + - T::from(T::from(1.0).unwrap()).unwrap()) + * a1 + + a2) + - *data + .get([tri_index(kx, ky), tri_index(q - 1, p), i]) + .unwrap() + * scale5 + * a3; + } + if ky > 0 { + for i in 0..data.shape()[2] { + let d = *data + .get([tri_index(kx, ky - 1), tri_index(q, p), i]) + .unwrap(); + *data + .get_mut([tri_index(kx, ky), tri_index(q + 1, p), i]) + .unwrap() += T::from(T::from(2.0).unwrap() * T::from(ky).unwrap()) + .unwrap() + * a1 + * d + * scale4; + } + } + } + } + } + } +} + +/// The number of polynomials +pub fn polynomial_count(cell_type: ReferenceCellType, degree: usize) -> usize { + match cell_type { + ReferenceCellType::Interval => degree + 1, + ReferenceCellType::Triangle => (degree + 1) * (degree + 2) / 2, + ReferenceCellType::Quadrilateral => (degree + 1) * (degree + 1), + _ => { + panic!("Unsupported cell type"); + } + } +} + +/// The total number of partial derivatives up to a give degree +pub fn derivative_count(cell_type: ReferenceCellType, derivatives: usize) -> usize { + match cell_type { + ReferenceCellType::Interval => derivatives + 1, + ReferenceCellType::Triangle => (derivatives + 1) * (derivatives + 2) / 2, + ReferenceCellType::Quadrilateral => (derivatives + 1) * (derivatives + 2) / 2, + _ => { + panic!("Unsupported cell type"); + } + } +} + +/// The shape of a table containing the values of Legendre polynomials +pub fn legendre_shape + Shape<2>>( + cell_type: ReferenceCellType, + points: &Array2, + degree: usize, + derivatives: usize, +) -> [usize; 3] { + [ + derivative_count(cell_type, derivatives), + polynomial_count(cell_type, degree), + points.shape()[0], + ] +} + +/// Tabulate orthonormal polynomials +pub fn tabulate_legendre_polynomials< + T: RlstScalar, + Array2: RandomAccessByRef<2, Item = T::Real> + Shape<2>, + Array3Mut: RandomAccessMut<3, Item = T> + RandomAccessByRef<3, Item = T> + Shape<3>, +>( + cell_type: ReferenceCellType, + points: &Array2, + degree: usize, + derivatives: usize, + data: &mut Array3Mut, +) { + match cell_type { + ReferenceCellType::Interval => { + tabulate_legendre_polynomials_interval(points, degree, derivatives, data) + } + ReferenceCellType::Triangle => { + tabulate_legendre_polynomials_triangle(points, degree, derivatives, data) + } + ReferenceCellType::Quadrilateral => { + tabulate_legendre_polynomials_quadrilateral(points, degree, derivatives, data) + } + _ => { + panic!("Unsupported cell type"); + } + }; +} + +#[cfg(test)] +mod test { + use super::*; + 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, [rule.npoints, 1]); + for (i, j) in rule.points.iter().enumerate() { + *points.get_mut([i, 0]).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); + } + } + } + } + + #[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, [rule.npoints, 2]); + for i in 0..rule.npoints { + for j in 0..2 { + *points.get_mut([i, j]).unwrap() = rule.points[i * 2 + j]; + } + } + + 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); + } + } + } + } + + #[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, [rule.npoints, 2]); + for i in 0..rule.npoints { + for j in 0..2 { + *points.get_mut([i, j]).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] + fn test_legendre_interval_derivative() { + let degree = 6; + + let epsilon = 1e-10; + let mut points = rlst_dynamic_array2!(f64, [20, 1]); + for i in 0..10 { + *points.get_mut([2 * i, 0]).unwrap() = i as f64 / 10.0; + *points.get_mut([2 * i + 1, 0]).unwrap() = points.get([2 * i, 0]).unwrap() + epsilon; + } + + let mut data = rlst_dynamic_array3!( + f64, + legendre_shape(ReferenceCellType::Interval, &points, degree, 1,) + ); + tabulate_legendre_polynomials(ReferenceCellType::Interval, &points, degree, 1, &mut data); + + for i in 0..degree + 1 { + for k in 0..points.shape()[0] / 2 { + assert_relative_eq!( + *data.get([1, i, 2 * k]).unwrap(), + (data.get([0, i, 2 * k + 1]).unwrap() - data.get([0, i, 2 * k]).unwrap()) + / epsilon, + epsilon = 1e-4 + ); + } + } + } + + #[test] + fn test_legendre_triangle_derivative() { + let degree = 6; + + let epsilon = 1e-10; + let mut points = rlst_dynamic_array2!(f64, [165, 2]); + let mut index = 0; + for i in 0..10 { + for j in 0..10 - i { + *points.get_mut([3 * index, 0]).unwrap() = i as f64 / 10.0; + *points.get_mut([3 * index, 1]).unwrap() = j as f64 / 10.0; + *points.get_mut([3 * index + 1, 0]).unwrap() = + *points.get([3 * index, 0]).unwrap() + epsilon; + *points.get_mut([3 * index + 1, 1]).unwrap() = *points.get([3 * index, 1]).unwrap(); + *points.get_mut([3 * index + 2, 0]).unwrap() = *points.get([3 * index, 0]).unwrap(); + *points.get_mut([3 * index + 2, 1]).unwrap() = + *points.get([3 * index, 1]).unwrap() + epsilon; + index += 1; + } + } + + let mut data = rlst_dynamic_array3!( + f64, + legendre_shape(ReferenceCellType::Triangle, &points, degree, 1,) + ); + tabulate_legendre_polynomials(ReferenceCellType::Triangle, &points, degree, 1, &mut data); + + for i in 0..degree + 1 { + for k in 0..points.shape()[0] / 3 { + assert_relative_eq!( + *data.get([1, i, 3 * k]).unwrap(), + (data.get([0, i, 3 * k + 1]).unwrap() - data.get([0, i, 3 * k]).unwrap()) + / epsilon, + epsilon = 1e-4 + ); + assert_relative_eq!( + *data.get([2, i, 3 * k]).unwrap(), + (data.get([0, i, 3 * k + 2]).unwrap() - data.get([0, i, 3 * k]).unwrap()) + / epsilon, + epsilon = 1e-4 + ); + } + } + } + + #[test] + fn test_legendre_quadrilateral_derivative() { + let degree = 6; + + let epsilon = 1e-10; + let mut points = rlst_dynamic_array2!(f64, [300, 2]); + for i in 0..10 { + for j in 0..10 { + let index = 10 * i + j; + *points.get_mut([3 * index, 0]).unwrap() = i as f64 / 10.0; + *points.get_mut([3 * index, 1]).unwrap() = j as f64 / 10.0; + *points.get_mut([3 * index + 1, 0]).unwrap() = + *points.get([3 * index, 0]).unwrap() + epsilon; + *points.get_mut([3 * index + 1, 1]).unwrap() = *points.get([3 * index, 1]).unwrap(); + *points.get_mut([3 * index + 2, 0]).unwrap() = *points.get([3 * index, 0]).unwrap(); + *points.get_mut([3 * index + 2, 1]).unwrap() = + *points.get([3 * index, 1]).unwrap() + epsilon; + } + } + + let mut data = rlst_dynamic_array3!( + f64, + legendre_shape(ReferenceCellType::Quadrilateral, &points, degree, 1,) + ); + tabulate_legendre_polynomials( + ReferenceCellType::Quadrilateral, + &points, + degree, + 1, + &mut data, + ); + + for i in 0..degree + 1 { + for k in 0..points.shape()[0] / 3 { + assert_relative_eq!( + *data.get([1, i, 3 * k]).unwrap(), + (data.get([0, i, 3 * k + 1]).unwrap() - data.get([0, i, 3 * k]).unwrap()) + / epsilon, + epsilon = 1e-4 + ); + assert_relative_eq!( + *data.get([2, i, 3 * k]).unwrap(), + (data.get([0, i, 3 * k + 2]).unwrap() - data.get([0, i, 3 * k]).unwrap()) + / epsilon, + epsilon = 1e-4 + ); + } + } + } + + #[test] + fn test_legendre_interval_against_known_polynomials() { + let degree = 3; + + let mut points = rlst_dynamic_array2!(f64, [11, 1]); + for i in 0..11 { + *points.get_mut([i, 0]).unwrap() = i as f64 / 10.0; + } + + let mut data = rlst_dynamic_array3!( + f64, + legendre_shape(ReferenceCellType::Interval, &points, degree, 3,) + ); + tabulate_legendre_polynomials(ReferenceCellType::Interval, &points, degree, 3, &mut data); + + for k in 0..points.shape()[0] { + let x = *points.get([k, 0]).unwrap(); + + // 0 => 1 + assert_relative_eq!(*data.get([0, 0, k]).unwrap(), 1.0, epsilon = 1e-12); + assert_relative_eq!(*data.get([1, 0, k]).unwrap(), 0.0, epsilon = 1e-12); + assert_relative_eq!(*data.get([2, 0, k]).unwrap(), 0.0, epsilon = 1e-12); + assert_relative_eq!(*data.get([3, 0, k]).unwrap(), 0.0, epsilon = 1e-12); + + // 1 => sqrt(3)*(2x - 1) + assert_relative_eq!( + *data.get([0, 1, k]).unwrap(), + f64::sqrt(3.0) * (2.0 * x - 1.0), + epsilon = 1e-12 + ); + assert_relative_eq!( + *data.get([1, 1, k]).unwrap(), + 2.0 * f64::sqrt(3.0), + epsilon = 1e-12 + ); + assert_relative_eq!(*data.get([2, 1, k]).unwrap(), 0.0, epsilon = 1e-12); + assert_relative_eq!(*data.get([3, 1, k]).unwrap(), 0.0, epsilon = 1e-12); + + // 2 => sqrt(5)*(6x^2 - 6x + 1) + assert_relative_eq!( + *data.get([0, 2, k]).unwrap(), + f64::sqrt(5.0) * (6.0 * x * x - 6.0 * x + 1.0), + epsilon = 1e-12 + ); + assert_relative_eq!( + *data.get([1, 2, k]).unwrap(), + f64::sqrt(5.0) * (12.0 * x - 6.0), + epsilon = 1e-12 + ); + assert_relative_eq!( + *data.get([2, 2, k]).unwrap(), + f64::sqrt(5.0) * 12.0, + epsilon = 1e-12 + ); + assert_relative_eq!(*data.get([3, 2, k]).unwrap(), 0.0, epsilon = 1e-12); + + // 3 => sqrt(7)*(20x^3 - 30x^2 + 12x - 1) + assert_relative_eq!( + *data.get([0, 3, k]).unwrap(), + f64::sqrt(7.0) * (20.0 * x * x * x - 30.0 * x * x + 12.0 * x - 1.0), + epsilon = 1e-12 + ); + assert_relative_eq!( + *data.get([1, 3, k]).unwrap(), + f64::sqrt(7.0) * (60.0 * x * x - 60.0 * x + 12.0), + epsilon = 1e-12 + ); + assert_relative_eq!( + *data.get([2, 3, k]).unwrap(), + f64::sqrt(7.0) * (120.0 * x - 60.0), + epsilon = 1e-12 + ); + assert_relative_eq!( + *data.get([3, 3, k]).unwrap(), + f64::sqrt(7.0) * 120.0, + epsilon = 1e-12 + ); + } + } + + #[test] + fn test_legendre_quadrilateral_against_known_polynomials() { + let degree = 2; + + let mut points = rlst_dynamic_array2!(f64, [121, 2]); + for i in 0..11 { + for j in 0..11 { + *points.get_mut([11 * i + j, 0]).unwrap() = i as f64 / 10.0; + *points.get_mut([11 * i + j, 1]).unwrap() = j as f64 / 10.0; + } + } + + let mut data = rlst_dynamic_array3!( + f64, + legendre_shape(ReferenceCellType::Quadrilateral, &points, degree, 1,) + ); + tabulate_legendre_polynomials( + ReferenceCellType::Quadrilateral, + &points, + degree, + 1, + &mut data, + ); + + for k in 0..points.shape()[0] { + let x = *points.get([k, 0]).unwrap(); + let y = *points.get([k, 1]).unwrap(); + + // 0 => 1 + assert_relative_eq!(*data.get([0, 0, k]).unwrap(), 1.0, epsilon = 1e-12); + assert_relative_eq!(*data.get([1, 0, k]).unwrap(), 0.0, epsilon = 1e-12); + assert_relative_eq!(*data.get([2, 0, k]).unwrap(), 0.0, epsilon = 1e-12); + + // 1 => sqrt(3)*(2x - 1) + assert_relative_eq!( + *data.get([0, 1, k]).unwrap(), + f64::sqrt(3.0) * (2.0 * x - 1.0), + epsilon = 1e-12 + ); + assert_relative_eq!( + *data.get([1, 1, k]).unwrap(), + 2.0 * f64::sqrt(3.0), + epsilon = 1e-12 + ); + assert_relative_eq!(*data.get([2, 1, k]).unwrap(), 0.0, epsilon = 1e-12); + + // 2 => sqrt(5)*(6x^2 - 6x + 1) + assert_relative_eq!( + *data.get([0, 2, k]).unwrap(), + f64::sqrt(5.0) * (6.0 * x * x - 6.0 * x + 1.0), + epsilon = 1e-12 + ); + assert_relative_eq!( + *data.get([1, 2, k]).unwrap(), + f64::sqrt(5.0) * (12.0 * x - 6.0), + epsilon = 1e-12 + ); + assert_relative_eq!(*data.get([2, 2, k]).unwrap(), 0.0, epsilon = 1e-12); + + // 3 => sqrt(3)*(2y - 1) + assert_relative_eq!( + *data.get([0, 3, k]).unwrap(), + f64::sqrt(3.0) * (2.0 * y - 1.0), + epsilon = 1e-12 + ); + + assert_relative_eq!(*data.get([1, 3, k]).unwrap(), 0.0, epsilon = 1e-12); + assert_relative_eq!( + *data.get([2, 3, k]).unwrap(), + 2.0 * f64::sqrt(3.0), + epsilon = 1e-12 + ); + + // 4 => 3*(2x - 1)*(2y - 1) + assert_relative_eq!( + *data.get([0, 4, k]).unwrap(), + 3.0 * (2.0 * x - 1.0) * (2.0 * y - 1.0), + epsilon = 1e-12 + ); + assert_relative_eq!( + *data.get([1, 4, k]).unwrap(), + 6.0 * (2.0 * y - 1.0), + epsilon = 1e-12 + ); + assert_relative_eq!( + *data.get([2, 4, k]).unwrap(), + 6.0 * (2.0 * x - 1.0), + epsilon = 1e-12 + ); + + // 5 => sqrt(15)*(6x^2 - 6x + 1)*(2y - 1) + assert_relative_eq!( + *data.get([0, 5, k]).unwrap(), + f64::sqrt(15.0) * (6.0 * x * x - 6.0 * x + 1.0) * (2.0 * y - 1.0), + epsilon = 1e-12 + ); + assert_relative_eq!( + *data.get([1, 5, k]).unwrap(), + f64::sqrt(15.0) * (12.0 * x - 6.0) * (2.0 * y - 1.0), + epsilon = 1e-12 + ); + assert_relative_eq!( + *data.get([2, 5, k]).unwrap(), + 2.0 * f64::sqrt(15.0) * (6.0 * x * x - 6.0 * x + 1.0), + epsilon = 1e-12 + ); + + // 6 => sqrt(5)*(6y^2 - 6y + 1) + assert_relative_eq!( + *data.get([0, 6, k]).unwrap(), + f64::sqrt(5.0) * (6.0 * y * y - 6.0 * y + 1.0), + epsilon = 1e-12 + ); + assert_relative_eq!(*data.get([1, 6, k]).unwrap(), 0.0, epsilon = 1e-12); + assert_relative_eq!( + *data.get([2, 6, k]).unwrap(), + f64::sqrt(5.0) * (12.0 * y - 6.0), + epsilon = 1e-12 + ); + + // 7 => sqrt(15)*(2x - 1)*(6y^2 - 6y + 1) + assert_relative_eq!( + *data.get([0, 7, k]).unwrap(), + f64::sqrt(15.0) * (2.0 * x - 1.0) * (6.0 * y * y - 6.0 * y + 1.0), + epsilon = 1e-12 + ); + assert_relative_eq!( + *data.get([1, 7, k]).unwrap(), + 2.0 * f64::sqrt(15.0) * (6.0 * y * y - 6.0 * y + 1.0), + epsilon = 1e-12 + ); + assert_relative_eq!( + *data.get([2, 7, k]).unwrap(), + f64::sqrt(15.0) * (2.0 * x - 1.0) * (12.0 * y - 6.0), + epsilon = 1e-12 + ); + + // 8 => 5*(6x^2 - 6x + 1)*(6y^2 - 6y + 1) + assert_relative_eq!( + *data.get([0, 8, k]).unwrap(), + 5.0 * (6.0 * x * x - 6.0 * x + 1.0) * (6.0 * y * y - 6.0 * y + 1.0), + epsilon = 1e-12 + ); + assert_relative_eq!( + *data.get([1, 8, k]).unwrap(), + 5.0 * (12.0 * x - 6.0) * (6.0 * y * y - 6.0 * y + 1.0), + epsilon = 1e-12 + ); + assert_relative_eq!( + *data.get([2, 8, k]).unwrap(), + 5.0 * (12.0 * y - 6.0) * (6.0 * x * x - 6.0 * x + 1.0), + epsilon = 1e-12 + ); + } + } +} diff --git a/src/reference_cell.rs b/src/reference_cell.rs new file mode 100644 index 0000000..5854b88 --- /dev/null +++ b/src/reference_cell.rs @@ -0,0 +1,522 @@ +//! Cell definitions + +use crate::types::ReferenceCellType; +use rlst::RlstScalar; + +/// The topological dimension of the cell +pub fn dim(cell: ReferenceCellType) -> usize { + match cell { + ReferenceCellType::Point => 0, + ReferenceCellType::Interval => 1, + ReferenceCellType::Triangle => 2, + ReferenceCellType::Quadrilateral => 2, + ReferenceCellType::Tetrahedron => 3, + ReferenceCellType::Hexahedron => 3, + ReferenceCellType::Prism => 3, + ReferenceCellType::Pyramid => 3, + } +} +/// Is the cell a simplex? +pub fn is_simplex(cell: ReferenceCellType) -> bool { + match cell { + ReferenceCellType::Point => true, + ReferenceCellType::Interval => true, + ReferenceCellType::Triangle => true, + ReferenceCellType::Quadrilateral => false, + ReferenceCellType::Tetrahedron => true, + ReferenceCellType::Hexahedron => false, + ReferenceCellType::Prism => false, + ReferenceCellType::Pyramid => false, + } +} + +/// The vertices of the reference cell +pub fn vertices>(cell: ReferenceCellType) -> Vec> { + let zero = T::from(0.0).unwrap(); + let one = T::from(1.0).unwrap(); + match cell { + ReferenceCellType::Point => vec![], + ReferenceCellType::Interval => vec![vec![zero], vec![one]], + ReferenceCellType::Triangle => vec![vec![zero, zero], vec![one, zero], vec![zero, one]], + ReferenceCellType::Quadrilateral => vec![ + vec![zero, zero], + vec![one, zero], + vec![zero, one], + vec![one, one], + ], + ReferenceCellType::Tetrahedron => vec![ + vec![zero, zero, zero], + vec![one, zero, zero], + vec![zero, one, zero], + vec![zero, zero, one], + ], + ReferenceCellType::Hexahedron => vec![ + vec![zero, zero, zero], + vec![one, zero, zero], + vec![zero, one, zero], + vec![one, one, zero], + vec![zero, zero, one], + vec![one, zero, one], + vec![zero, one, one], + vec![one, one, one], + ], + ReferenceCellType::Prism => vec![ + vec![zero, zero, zero], + vec![one, zero, zero], + vec![zero, one, zero], + vec![zero, zero, one], + vec![one, zero, one], + vec![zero, one, one], + ], + ReferenceCellType::Pyramid => vec![ + vec![zero, zero, zero], + vec![one, zero, zero], + vec![zero, one, zero], + vec![one, one, zero], + vec![zero, zero, one], + ], + } +} + +/// The midpoint of the cell +pub fn midpoint>(cell: ReferenceCellType) -> Vec { + let half = T::from(0.5).unwrap(); + let third = T::from(1.0).unwrap() / T::from(3.0).unwrap(); + match cell { + ReferenceCellType::Point => vec![], + ReferenceCellType::Interval => vec![half], + ReferenceCellType::Triangle => vec![third; 2], + ReferenceCellType::Quadrilateral => vec![half; 2], + ReferenceCellType::Tetrahedron => vec![T::from(1.0).unwrap() / T::from(6.0).unwrap(); 3], + ReferenceCellType::Hexahedron => vec![half; 3], + ReferenceCellType::Prism => vec![third, third, half], + ReferenceCellType::Pyramid => vec![ + T::from(0.4).unwrap(), + T::from(0.4).unwrap(), + T::from(0.2).unwrap(), + ], + } +} + +/// The edges of the reference cell +pub fn edges(cell: ReferenceCellType) -> Vec> { + match cell { + ReferenceCellType::Point => vec![], + ReferenceCellType::Interval => vec![vec![0, 1]], + ReferenceCellType::Triangle => vec![vec![1, 2], vec![0, 2], vec![0, 1]], + ReferenceCellType::Quadrilateral => vec![vec![0, 1], vec![0, 2], vec![1, 3], vec![2, 3]], + ReferenceCellType::Tetrahedron => vec![ + vec![2, 3], + vec![1, 3], + vec![1, 2], + vec![0, 3], + vec![0, 2], + vec![0, 1], + ], + ReferenceCellType::Hexahedron => vec![ + vec![0, 1], + vec![0, 2], + vec![0, 4], + vec![1, 3], + vec![1, 5], + vec![2, 3], + vec![2, 6], + vec![3, 7], + vec![4, 5], + vec![4, 6], + vec![5, 7], + vec![6, 7], + ], + ReferenceCellType::Prism => vec![ + vec![0, 1], + vec![0, 2], + vec![0, 3], + vec![1, 2], + vec![1, 4], + vec![2, 5], + vec![3, 4], + vec![3, 5], + vec![4, 5], + ], + ReferenceCellType::Pyramid => vec![ + vec![0, 1], + vec![0, 2], + vec![0, 4], + vec![1, 3], + vec![1, 4], + vec![2, 3], + vec![2, 4], + vec![3, 4], + ], + } +} + +/// The faces of the reference cell +pub fn faces(cell: ReferenceCellType) -> Vec> { + match cell { + ReferenceCellType::Point => vec![], + ReferenceCellType::Interval => vec![], + ReferenceCellType::Triangle => vec![vec![0, 1, 2]], + ReferenceCellType::Quadrilateral => vec![vec![0, 1, 2, 3]], + ReferenceCellType::Tetrahedron => { + vec![vec![1, 2, 3], vec![0, 2, 3], vec![0, 1, 3], vec![0, 1, 2]] + } + ReferenceCellType::Hexahedron => vec![ + vec![0, 1, 2, 3], + vec![0, 1, 4, 5], + vec![0, 2, 4, 6], + vec![1, 3, 5, 7], + vec![2, 3, 6, 7], + vec![4, 5, 6, 7], + ], + ReferenceCellType::Prism => vec![ + vec![0, 1, 2], + vec![0, 1, 3, 4], + vec![0, 2, 3, 5], + vec![1, 2, 4, 5], + vec![3, 4, 5], + ], + ReferenceCellType::Pyramid => vec![ + vec![0, 1, 2, 3], + vec![0, 1, 4], + vec![0, 2, 4], + vec![1, 3, 4], + vec![2, 3, 4], + ], + } +} + +/// The types of the subentities of the reference cell +pub fn entity_types(cell: ReferenceCellType) -> Vec> { + match cell { + ReferenceCellType::Point => vec![vec![ReferenceCellType::Point], vec![], vec![], vec![]], + ReferenceCellType::Interval => vec![ + vec![ReferenceCellType::Point; 2], + vec![ReferenceCellType::Interval], + vec![], + vec![], + ], + ReferenceCellType::Triangle => vec![ + vec![ReferenceCellType::Point; 3], + vec![ReferenceCellType::Interval; 3], + vec![ReferenceCellType::Triangle], + vec![], + ], + ReferenceCellType::Quadrilateral => vec![ + vec![ReferenceCellType::Point; 4], + vec![ReferenceCellType::Interval; 4], + vec![ReferenceCellType::Quadrilateral], + vec![], + ], + ReferenceCellType::Tetrahedron => vec![ + vec![ReferenceCellType::Point; 4], + vec![ReferenceCellType::Interval; 6], + vec![ReferenceCellType::Triangle; 4], + vec![ReferenceCellType::Tetrahedron], + ], + ReferenceCellType::Hexahedron => vec![ + vec![ReferenceCellType::Point; 8], + vec![ReferenceCellType::Interval; 12], + vec![ReferenceCellType::Quadrilateral; 6], + vec![ReferenceCellType::Hexahedron], + ], + ReferenceCellType::Prism => vec![ + vec![ReferenceCellType::Point; 6], + vec![ReferenceCellType::Interval; 9], + vec![ + ReferenceCellType::Triangle, + ReferenceCellType::Quadrilateral, + ReferenceCellType::Quadrilateral, + ReferenceCellType::Quadrilateral, + ReferenceCellType::Triangle, + ], + vec![ReferenceCellType::Prism], + ], + ReferenceCellType::Pyramid => vec![ + vec![ReferenceCellType::Point; 5], + vec![ReferenceCellType::Interval; 8], + vec![ + ReferenceCellType::Quadrilateral, + ReferenceCellType::Triangle, + ReferenceCellType::Triangle, + ReferenceCellType::Triangle, + ReferenceCellType::Triangle, + ], + vec![ReferenceCellType::Pyramid], + ], + } +} + +/// The number of subentities of each dimension +pub fn entity_counts(cell: ReferenceCellType) -> Vec { + match cell { + ReferenceCellType::Point => vec![1, 0, 0, 0], + ReferenceCellType::Interval => vec![2, 1, 0, 0], + ReferenceCellType::Triangle => vec![3, 3, 1, 0], + ReferenceCellType::Quadrilateral => vec![4, 4, 1, 0], + ReferenceCellType::Tetrahedron => vec![4, 6, 4, 1], + ReferenceCellType::Hexahedron => vec![8, 12, 6, 1], + ReferenceCellType::Prism => vec![6, 9, 5, 1], + ReferenceCellType::Pyramid => vec![5, 8, 5, 1], + } +} + +/// The connectivity of the reference cell +/// +/// The indices of the result are \[i\]\[j\]\[k\]\[l\] +pub fn connectivity(cell: ReferenceCellType) -> Vec>>> { + match cell { + ReferenceCellType::Point => vec![vec![vec![vec![0]]]], + ReferenceCellType::Interval => vec![ + vec![vec![vec![0], vec![0]], vec![vec![1], vec![0]]], + vec![vec![vec![0, 1], vec![0]]], + ], + ReferenceCellType::Triangle => vec![ + vec![ + vec![vec![0], vec![1, 2], vec![0]], + vec![vec![1], vec![0, 2], vec![0]], + vec![vec![2], vec![0, 1], vec![0]], + ], + vec![ + vec![vec![1, 2], vec![0], vec![0]], + vec![vec![0, 2], vec![1], vec![0]], + vec![vec![0, 1], vec![2], vec![0]], + ], + vec![vec![vec![0, 1, 2], vec![0, 1, 2], vec![0]]], + ], + ReferenceCellType::Quadrilateral => vec![ + vec![ + vec![vec![0], vec![0, 1], vec![0]], + vec![vec![1], vec![0, 2], vec![0]], + vec![vec![2], vec![1, 3], vec![0]], + vec![vec![3], vec![2, 3], vec![0]], + ], + vec![ + vec![vec![0, 1], vec![0], vec![0]], + vec![vec![0, 2], vec![1], vec![0]], + vec![vec![1, 3], vec![2], vec![0]], + vec![vec![2, 3], vec![3], vec![0]], + ], + vec![vec![vec![0, 1, 2, 3], vec![0, 1, 2, 3], vec![0]]], + ], + ReferenceCellType::Tetrahedron => vec![ + vec![ + vec![vec![0], vec![3, 4, 5], vec![1, 2, 3], vec![0]], + vec![vec![1], vec![1, 2, 5], vec![0, 2, 3], vec![0]], + vec![vec![2], vec![0, 2, 4], vec![0, 1, 3], vec![0]], + vec![vec![3], vec![0, 1, 3], vec![0, 1, 2], vec![0]], + ], + vec![ + vec![vec![2, 3], vec![0], vec![0, 1], vec![0]], + vec![vec![1, 3], vec![1], vec![0, 2], vec![0]], + vec![vec![1, 2], vec![2], vec![0, 3], vec![0]], + vec![vec![0, 3], vec![3], vec![1, 2], vec![0]], + vec![vec![0, 2], vec![4], vec![1, 3], vec![0]], + vec![vec![0, 1], vec![5], vec![2, 3], vec![0]], + ], + vec![ + vec![vec![1, 2, 3], vec![0, 1, 2], vec![0], vec![0]], + vec![vec![0, 2, 3], vec![0, 3, 4], vec![1], vec![0]], + vec![vec![0, 1, 3], vec![1, 3, 5], vec![2], vec![0]], + vec![vec![0, 1, 2], vec![2, 4, 5], vec![3], vec![0]], + ], + vec![vec![ + vec![0, 1, 2, 3], + vec![0, 1, 2, 3, 4, 5], + vec![0, 1, 2, 3], + vec![0], + ]], + ], + ReferenceCellType::Hexahedron => vec![ + vec![ + vec![vec![0], vec![0, 1, 2], vec![0, 1, 2], vec![0]], + vec![vec![1], vec![0, 3, 4], vec![0, 1, 3], vec![0]], + vec![vec![2], vec![1, 5, 6], vec![0, 2, 4], vec![0]], + vec![vec![3], vec![3, 5, 7], vec![0, 3, 4], vec![0]], + vec![vec![4], vec![2, 8, 9], vec![1, 2, 5], vec![0]], + vec![vec![5], vec![4, 8, 10], vec![1, 3, 5], vec![0]], + vec![vec![6], vec![6, 9, 11], vec![2, 4, 5], vec![0]], + vec![vec![7], vec![7, 10, 11], vec![3, 4, 5], vec![0]], + ], + vec![ + vec![vec![0, 1], vec![0], vec![0, 1], vec![0]], + vec![vec![0, 2], vec![1], vec![0, 2], vec![0]], + vec![vec![0, 4], vec![2], vec![1, 2], vec![0]], + vec![vec![1, 3], vec![3], vec![0, 3], vec![0]], + vec![vec![1, 5], vec![4], vec![1, 3], vec![0]], + vec![vec![2, 3], vec![5], vec![0, 4], vec![0]], + vec![vec![2, 6], vec![6], vec![2, 4], vec![0]], + vec![vec![3, 7], vec![7], vec![3, 4], vec![0]], + vec![vec![4, 5], vec![8], vec![1, 5], vec![0]], + vec![vec![4, 6], vec![9], vec![2, 5], vec![0]], + vec![vec![5, 7], vec![10], vec![3, 5], vec![0]], + vec![vec![6, 7], vec![11], vec![4, 5], vec![0]], + ], + vec![ + vec![vec![0, 1, 2, 3], vec![0, 1, 3, 5], vec![0], vec![0]], + vec![vec![0, 1, 4, 5], vec![0, 2, 4, 8], vec![1], vec![0]], + vec![vec![0, 2, 4, 6], vec![1, 2, 6, 9], vec![2], vec![0]], + vec![vec![1, 3, 5, 7], vec![3, 4, 7, 10], vec![3], vec![0]], + vec![vec![2, 3, 6, 7], vec![5, 6, 7, 11], vec![4], vec![0]], + vec![vec![4, 5, 6, 7], vec![8, 9, 10, 11], vec![5], vec![0]], + ], + vec![vec![ + vec![0, 1, 2, 3, 4, 5, 6, 7], + vec![0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11], + vec![0, 1, 2, 3, 4, 5], + vec![0], + ]], + ], + ReferenceCellType::Prism => vec![ + vec![ + vec![vec![0], vec![0, 1, 2], vec![0, 1, 2], vec![0]], + vec![vec![1], vec![0, 3, 4], vec![0, 1, 3], vec![0]], + vec![vec![2], vec![1, 3, 5], vec![0, 2, 3], vec![0]], + vec![vec![3], vec![2, 6, 7], vec![1, 2, 4], vec![0]], + vec![vec![4], vec![4, 6, 8], vec![1, 3, 4], vec![0]], + vec![vec![5], vec![5, 7, 8], vec![2, 3, 4], vec![0]], + ], + vec![ + vec![vec![0, 1], vec![0], vec![0, 1], vec![0]], + vec![vec![0, 2], vec![1], vec![0, 2], vec![0]], + vec![vec![0, 3], vec![2], vec![1, 2], vec![0]], + vec![vec![1, 2], vec![3], vec![0, 3], vec![0]], + vec![vec![1, 4], vec![4], vec![1, 3], vec![0]], + vec![vec![2, 5], vec![5], vec![2, 3], vec![0]], + vec![vec![3, 4], vec![6], vec![1, 4], vec![0]], + vec![vec![3, 5], vec![7], vec![2, 4], vec![0]], + vec![vec![4, 5], vec![8], vec![3, 4], vec![0]], + ], + vec![ + vec![vec![0, 1, 2], vec![0, 1, 3], vec![0], vec![0]], + vec![vec![0, 1, 3, 4], vec![0, 2, 4, 6], vec![1], vec![0]], + vec![vec![0, 2, 3, 5], vec![1, 2, 5, 7], vec![2], vec![0]], + vec![vec![1, 2, 4, 5], vec![3, 4, 5, 8], vec![3], vec![0]], + vec![vec![3, 4, 5], vec![6, 7, 8], vec![4], vec![0]], + ], + vec![vec![ + vec![0, 1, 2, 3, 4, 5], + vec![0, 1, 2, 3, 4, 5, 6, 7, 8], + vec![0, 1, 2, 3, 4], + vec![0], + ]], + ], + ReferenceCellType::Pyramid => vec![ + vec![ + vec![vec![0], vec![0, 1, 2], vec![0, 1, 2], vec![0]], + vec![vec![1], vec![0, 3, 4], vec![0, 1, 3], vec![0]], + vec![vec![2], vec![1, 5, 6], vec![0, 2, 4], vec![0]], + vec![vec![3], vec![3, 5, 7], vec![0, 3, 4], vec![0]], + vec![vec![4], vec![2, 4, 6, 7], vec![1, 2, 3, 4], vec![0]], + ], + vec![ + vec![vec![0, 1], vec![0], vec![0, 1], vec![0]], + vec![vec![0, 2], vec![1], vec![0, 2], vec![0]], + vec![vec![0, 4], vec![2], vec![1, 2], vec![0]], + vec![vec![1, 3], vec![3], vec![0, 3], vec![0]], + vec![vec![1, 4], vec![4], vec![1, 3], vec![0]], + vec![vec![2, 3], vec![5], vec![0, 4], vec![0]], + vec![vec![2, 4], vec![6], vec![2, 4], vec![0]], + vec![vec![3, 4], vec![7], vec![3, 4], vec![0]], + ], + vec![ + vec![vec![0, 1, 2, 3], vec![0, 1, 3, 5], vec![0], vec![0]], + vec![vec![0, 1, 4], vec![0, 2, 4], vec![1], vec![0]], + vec![vec![0, 2, 4], vec![1, 2, 6], vec![2], vec![0]], + vec![vec![1, 3, 4], vec![3, 4, 7], vec![3], vec![0]], + vec![vec![2, 3, 4], vec![5, 6, 7], vec![4], vec![0]], + ], + vec![vec![ + vec![0, 1, 2, 3, 4], + vec![0, 1, 2, 3, 4, 5, 6, 7], + vec![0, 1, 2, 3, 4], + vec![0], + ]], + ], + } +} + +#[cfg(test)] +mod test { + use super::*; + use paste::paste; + + macro_rules! test_cell { + + ($($cell:ident),+) => { + + $( + paste! { + + #[test] + fn []() { + let v = vertices::(ReferenceCellType::[<$cell>]); + let d = dim(ReferenceCellType::[<$cell>]); + let ec = entity_counts(ReferenceCellType::[<$cell>]); + let et = entity_types(ReferenceCellType::[<$cell>]); + let conn = connectivity(ReferenceCellType::[<$cell>]); + for i in 0..d+1 { + assert_eq!(ec[i], et[i].len()); + assert_eq!(ec[i], conn[i].len()); + } + assert_eq!(ec[0], v.len()); + for i in v { + assert_eq!(i.len(), d); + } + + for v_n in 0..ec[0] { + let v = &conn[0][v_n][0]; + assert_eq!(v, &[v_n]); + } + for e_n in 0..ec[1] { + let vs = &conn[1][e_n][0]; + assert_eq!(vs, &edges(ReferenceCellType::[<$cell>])[e_n]); + } + for f_n in 0..ec[2] { + let vs = &conn[2][f_n][0]; + assert_eq!(vs, &faces(ReferenceCellType::[<$cell>])[f_n]); + } + + for e_dim in 0..d { + for e_n in 0..ec[e_dim] { + let e_vertices = &conn[e_dim][e_n][0]; + for c_dim in 0..d + 1 { + let connectivity = &conn[e_dim][e_n][c_dim]; + if e_dim == c_dim { + assert_eq!(connectivity.len(), 1); + assert_eq!(connectivity[0], e_n) + } else { + for c_n in connectivity { + let c_vertices = &conn[c_dim][*c_n][0]; + if e_dim < c_dim { + for i in e_vertices { + assert!(c_vertices.contains(&i)); + } + } else { + for i in c_vertices { + assert!(e_vertices.contains(&i)); + } + } + assert!(connectivity.contains(&c_n)); + } + } + } + } + } + } + + } + )* + }; + } + + test_cell!( + Interval, + Triangle, + Quadrilateral, + Tetrahedron, + Hexahedron, + Prism, + Pyramid + ); +} diff --git a/src/traits.rs b/src/traits.rs new file mode 100644 index 0000000..ee3c79e --- /dev/null +++ b/src/traits.rs @@ -0,0 +1,5 @@ +//! Traits + +mod element; + +pub use element::{ElementFamily, FiniteElement}; diff --git a/src/traits/element.rs b/src/traits/element.rs new file mode 100644 index 0000000..bc9f7fd --- /dev/null +++ b/src/traits/element.rs @@ -0,0 +1,59 @@ +//! Finite element trait +use rlst::{RandomAccessByRef, RandomAccessMut, RlstScalar, Shape}; +use std::fmt::Debug; +use std::hash::Hash; + +pub trait FiniteElement { + //! A finite element defined on a reference cell + /// The scalar type + type T: RlstScalar; + /// Cell type + type CellType: Debug + PartialEq + Eq + Clone + Copy + Hash; + /// Map type + type MapType: Debug + PartialEq + Eq + Clone + Copy + Hash; + + /// The reference cell type + fn cell_type(&self) -> Self::CellType; + + /// The highest degree polynomial in the element's polynomial set + fn embedded_superdegree(&self) -> usize; + + /// The number of basis functions + fn dim(&self) -> usize; + + /// The value shape + fn value_shape(&self) -> &[usize]; + + /// The value size + fn value_size(&self) -> usize; + + /// Tabulate the values of the basis functions and their derivatives at a set of points + fn tabulate::Real> + Shape<2>>( + &self, + points: &Array2, + nderivs: usize, + data: &mut impl RandomAccessMut<4, Item = Self::T>, + ); + + /// The DOFs that are associated with a subentity of the reference cell + fn entity_dofs(&self, entity_dim: usize, entity_number: usize) -> Option<&[usize]>; + + /// The push forward / pull back map to use for this element + fn map_type(&self) -> Self::MapType; + + /// Get the required shape for a tabulation array + fn tabulate_array_shape(&self, nderivs: usize, npoints: usize) -> [usize; 4]; +} + +pub trait ElementFamily { + //! A family of finite elements + /// The scalar type + type T: RlstScalar; + /// Cell type + type CellType: Debug + PartialEq + Eq + Clone + Copy + Hash; + /// The finite element type + type FiniteElement: FiniteElement; + + /// Get an elenent for a cell type + fn element(&self, cell_type: Self::CellType) -> Self::FiniteElement; +} diff --git a/src/types.rs b/src/types.rs new file mode 100644 index 0000000..7bb0ed4 --- /dev/null +++ b/src/types.rs @@ -0,0 +1,112 @@ +//! Types + +/// Continuity type +#[derive(Debug, PartialEq, Eq, Clone, Copy, Hash)] +#[repr(u8)] +pub enum Continuity { + /// The element has standard continuity between cells + /// + /// For some element, this option does not indicate that the values are fully continuous. + /// For example, for Raviart-Thomas elements it only indicates that the normal components + /// are continuous across edges + Continuous = 0, + /// The element is discontinuous betweeen cells + Discontinuous = 1, +} + +/// The map type used by an element +#[derive(Debug, PartialEq, Eq, Clone, Copy, Hash)] +#[repr(u8)] +pub enum MapType { + /// Identity map + Identity = 0, + /// Covariant Piola map + /// + /// This map is used by H(curl) elements + CovariantPiola = 1, + /// Contravariant Piola map + /// + /// This map is used by H(div) elements + ContravariantPiola = 2, + /// L2 Piola map + L2Piola = 3, +} + +/// The type of a reference cell +#[derive(Debug, PartialEq, Eq, Clone, Copy, Hash)] +#[repr(u8)] +pub enum ReferenceCellType { + /// A point + Point = 0, + /// A line interval + Interval = 1, + /// A triangle + Triangle = 2, + /// A quadrilateral + Quadrilateral = 3, + /// A tetrahedron (whose faces are all triangles) + Tetrahedron = 4, + /// A hexahedron (whose faces are all quadrilaterals) + Hexahedron = 5, + /// A triangular prism + Prism = 6, + /// A square-based pyramid + Pyramid = 7, +} + +impl ReferenceCellType { + /// Create a reference cell type from a u8 + pub fn from(i: u8) -> Option { + match i { + 0 => Some(ReferenceCellType::Point), + 1 => Some(ReferenceCellType::Interval), + 2 => Some(ReferenceCellType::Triangle), + 3 => Some(ReferenceCellType::Quadrilateral), + 4 => Some(ReferenceCellType::Tetrahedron), + 5 => Some(ReferenceCellType::Hexahedron), + 6 => Some(ReferenceCellType::Prism), + 7 => Some(ReferenceCellType::Pyramid), + _ => None, + } + } +} + +#[cfg(test)] +mod test { + use super::*; + #[test] + fn test_reference_cell_type() { + assert_eq!( + ReferenceCellType::Point, + ReferenceCellType::from(ReferenceCellType::Point as u8).unwrap() + ); + assert_eq!( + ReferenceCellType::Interval, + ReferenceCellType::from(ReferenceCellType::Interval as u8).unwrap() + ); + assert_eq!( + ReferenceCellType::Triangle, + ReferenceCellType::from(ReferenceCellType::Triangle as u8).unwrap() + ); + assert_eq!( + ReferenceCellType::Quadrilateral, + ReferenceCellType::from(ReferenceCellType::Quadrilateral as u8).unwrap() + ); + assert_eq!( + ReferenceCellType::Tetrahedron, + ReferenceCellType::from(ReferenceCellType::Tetrahedron as u8).unwrap() + ); + assert_eq!( + ReferenceCellType::Hexahedron, + ReferenceCellType::from(ReferenceCellType::Hexahedron as u8).unwrap() + ); + assert_eq!( + ReferenceCellType::Prism, + ReferenceCellType::from(ReferenceCellType::Prism as u8).unwrap() + ); + assert_eq!( + ReferenceCellType::Pyramid, + ReferenceCellType::from(ReferenceCellType::Pyramid as u8).unwrap() + ); + } +}