diff --git a/Cargo.toml b/Cargo.toml index 9531684..eca6102 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -44,11 +44,16 @@ blas-src = { version = "0.14", features = ["blis"]} paste = "1.*" approx = "0.5" +[build] +rustdocflags = [ "--html-in-header", "katex-header.html" ] + + [build-dependencies] cbindgen = "0.29.*" [package.metadata.docs.rs] cargo-args = ["-Zunstable-options", "-Zrustdoc-scrape-examples"] +rustdoc-args = [ "--html-in-header", "katex-header.html" ] [lints.clippy] wildcard_imports = "forbid" diff --git a/examples/element_family.rs b/examples/element_family.rs index 60c29b8..26874d4 100644 --- a/examples/element_family.rs +++ b/examples/element_family.rs @@ -11,7 +11,7 @@ fn main() { let element = family.element(ReferenceCellType::Triangle); println!("Cell: {:?}", element.cell_type()); - // Get the element in the family on a triangle + // Get the element in the family on a quadrilateral let element = family.element(ReferenceCellType::Quadrilateral); println!("Cell: {:?}", element.cell_type()); } diff --git a/katex-header.html b/katex-header.html new file mode 100644 index 0000000..da412ac --- /dev/null +++ b/katex-header.html @@ -0,0 +1,19 @@ + + + + diff --git a/src/bindings.rs b/src/bindings.rs index 9e6bf53..cae544d 100644 --- a/src/bindings.rs +++ b/src/bindings.rs @@ -686,7 +686,7 @@ pub mod ciarlet { field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]) )] pub fn ciarlet_element_embedded_superdegree(element: &E) -> usize { - element.embedded_superdegree() + element.lagrange_superdegree() } #[concretise_types( diff --git a/src/ciarlet.rs b/src/ciarlet.rs index 22197b0..958ac20 100644 --- a/src/ciarlet.rs +++ b/src/ciarlet.rs @@ -72,7 +72,13 @@ impl CiarletElement where DynArray: Inverse>, { - /// Create a Ciarlet element + /// Create a Ciarlet element. + /// + /// This should not be used directly. Instead users should call the `create` + /// function for one of the following special cases of a general Ciarlet element. + /// - [crate::ciarlet::lagrange::create]: Create a new Lagrange element. + /// - [crate::ciarlet::nedelec::create]: Create a new Nedelec element. + /// - [crate::ciarlet::raviart_thomas::create]: Create a Raviart-Thomas element. #[allow(clippy::too_many_arguments)] pub fn create( family_name: String, @@ -603,7 +609,7 @@ impl FiniteElement for CiarletElement { fn cell_type(&self) -> ReferenceCellType { self.cell_type } - fn embedded_superdegree(&self) -> usize { + fn lagrange_superdegree(&self) -> usize { self.embedded_superdegree } fn dim(&self) -> usize { diff --git a/src/ciarlet/lagrange.rs b/src/ciarlet/lagrange.rs index d4ea274..87890cc 100644 --- a/src/ciarlet/lagrange.rs +++ b/src/ciarlet/lagrange.rs @@ -10,7 +10,7 @@ use rlst::dense::linalg::lapack::interface::{getrf::Getrf, getri::Getri}; use rlst::{RlstScalar, rlst_dynamic_array}; use std::marker::PhantomData; -/// Create a Lagrange element +/// Create a Lagrange element. pub fn create( cell_type: ReferenceCellType, degree: usize, @@ -268,7 +268,9 @@ pub fn create( ) } -/// Lagrange element family +/// Lagrange element family. +/// +/// A factory structure to create new Lagrange elements. pub struct LagrangeElementFamily { degree: usize, continuity: Continuity, @@ -276,7 +278,7 @@ pub struct LagrangeElementFamily { } impl LagrangeElementFamily { - /// Create new family + /// Create new family with given `degree` and `continuity`. pub fn new(degree: usize, continuity: Continuity) -> Self { Self { degree, diff --git a/src/ciarlet/nedelec.rs b/src/ciarlet/nedelec.rs index 67ecde2..c03e96d 100644 --- a/src/ciarlet/nedelec.rs +++ b/src/ciarlet/nedelec.rs @@ -640,6 +640,8 @@ pub fn create( } /// Nedelec (first kind) element family +/// +/// A factory structure to create new Nedelec elements. pub struct NedelecFirstKindElementFamily { degree: usize, continuity: Continuity, @@ -647,7 +649,7 @@ pub struct NedelecFirstKindElementFamily { } impl NedelecFirstKindElementFamily { - /// Create new family + /// Create new family with given `degree` and `continuity`. pub fn new(degree: usize, continuity: Continuity) -> Self { Self { degree, diff --git a/src/ciarlet/raviart_thomas.rs b/src/ciarlet/raviart_thomas.rs index c91ad6e..c437f63 100644 --- a/src/ciarlet/raviart_thomas.rs +++ b/src/ciarlet/raviart_thomas.rs @@ -483,6 +483,8 @@ pub fn create( } /// Raviart-Thomas element family +/// +/// A factory structure to create new Raviart-Thomas elements. pub struct RaviartThomasElementFamily { degree: usize, continuity: Continuity, @@ -490,7 +492,7 @@ pub struct RaviartThomasElementFamily { } impl RaviartThomasElementFamily { - /// Create new family + /// Create new family with given `degree` and `continuity`. pub fn new(degree: usize, continuity: Continuity) -> Self { Self { degree, diff --git a/src/lib.rs b/src/lib.rs index e460bbe..cf5c53b 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,4 +1,41 @@ -//! n-dimensional grid +//! A library for the definition and tabulation of finite elements in Rust. +//! +//! `ndelement` provides definition of many frequently used low and high order finite elements +//! and provides routines for the tabulation of their values. +//! +//! The following presents a simple example of how to use `ndelement`. +//! +//! ``` +//! use ndelement::ciarlet::LagrangeElementFamily; +//! use ndelement::traits::{ElementFamily, FiniteElement}; +//! use ndelement::types::{Continuity, ReferenceCellType}; +//! use rlst::DynArray; +//! +//! // Create the degree 2 Lagrange element family. A family is a set of finite elements with the +//! // same family type, degree, and continuity across a set of cells +//! let family = LagrangeElementFamily::::new(2, Continuity::Standard); +//! +//! // Get the element in the family on a triangle +//! let element = family.element(ReferenceCellType::Triangle); +//! println!("Cell: {:?}", element.cell_type()); +//! +//! // Get the element in the family on a quadrilateral +//! let element = family.element(ReferenceCellType::Quadrilateral); +//! println!("Cell: {:?}", element.cell_type()); +//! +//! // Create an array to store the basis function values +//! let mut basis_values = DynArray::::from_shape(element.tabulate_array_shape(0, 1)); +//! // Create array containing the point [1/2, 1/2] +//! let mut points = DynArray::::from_shape([2, 1]); +//! points[[0, 0]] = 1.0 / 2.0; +//! points[[1, 0]] = 1.0 / 2.0; +//! // Tabulate the element's basis functions at the point +//! element.tabulate(&points, 0, &mut basis_values); +//! println!( +//! "The values of the basis functions at the point (1/2, 1/2) are: {:?}", +//! basis_values.data() +//! ); +//! ``` #![cfg_attr(feature = "strict", deny(warnings), deny(unused_crate_dependencies))] #![warn(missing_docs)] diff --git a/src/map.rs b/src/map.rs index 1ab6866..f4af76c 100644 --- a/src/map.rs +++ b/src/map.rs @@ -3,6 +3,9 @@ use crate::traits::Map; use rlst::{Array, MutableArrayImpl, RlstScalar, ValueArrayImpl}; /// Identity map +/// +/// An identity map pushes values from the reference to the physical +/// cell without modifying them. pub struct IdentityMap {} impl Map for IdentityMap { @@ -51,7 +54,17 @@ impl Map for IdentityMap { } } -/// CovariantPiola map +/// Covariant Piola map. +/// +/// Let $F$ be the map from the reference cell to the physical cell +/// and let $J$ be its Jacobian. Let $\Phi$ be the function values +/// on the reference cell. The covariant Piola map is defined by +/// $$ +/// J^{-T}\Phi\circ F^{-1} +/// $$ +/// The coviariant Piola map preserves tangential continuity. If $J$ +/// is a rectangular matrix then the pseudo-inverse is used instead of +/// the inverse. pub struct CovariantPiolaMap {} impl Map for CovariantPiolaMap { @@ -138,7 +151,15 @@ impl Map for CovariantPiolaMap { } } -/// ContravariantPiola map +/// Contravariant Piola map. +/// +/// Let $F$ be the map from the reference cell to the physical cell +/// and let $J$ be its Jacobian. Let $\Phi$ be the function values +/// on the reference cell. The contravariant Piola map is defined by +/// $$ +/// \frac{1}{\sqrt{\det{J^TJ}}}J\Phi\circ F^{-1} +/// $$ +/// The contravariant Piola map preserves normal continuity. pub struct ContravariantPiolaMap {} impl Map for ContravariantPiolaMap { diff --git a/src/math.rs b/src/math.rs index aad6236..5aa77d2 100644 --- a/src/math.rs +++ b/src/math.rs @@ -2,32 +2,7 @@ use rlst::{Array, MutableArrayImpl, RlstScalar, ValueArrayImpl}; /// Orthogonalise the rows of a matrix, starting with the row numbered `start` -pub fn orthogonalise>( - mat: &mut Array, - start: usize, -) { - for row in start..mat.shape()[0] { - let norm = (0..mat.shape()[1]) - .map(|i| mat.get([row, i]).unwrap().powi(2)) - .sum::() - .sqrt(); - for i in 0..mat.shape()[1] { - *mat.get_mut([row, i]).unwrap() /= norm; - } - for r in row + 1..mat.shape()[0] { - let dot = (0..mat.shape()[1]) - .map(|i| *mat.get([row, i]).unwrap() * *mat.get([r, i]).unwrap()) - .sum::(); - for i in 0..mat.shape()[1] { - let sub = dot * *mat.get([row, i]).unwrap(); - *mat.get_mut([r, i]).unwrap() -= sub; - } - } - } -} - -/// Orthogonalise the rows of a matrix, starting with the row numbered `start` -pub fn orthogonalise_3>( +pub(crate) fn orthogonalise_3>( mat: &mut Array, start: usize, ) { @@ -186,6 +161,31 @@ pub fn apply_matrix>( } } +// /// Orthogonalise the rows of a matrix, starting with the row numbered `start` +// pub(crate) fn orthogonalise>( +// mat: &mut Array, +// start: usize, +// ) { +// for row in start..mat.shape()[0] { +// let norm = (0..mat.shape()[1]) +// .map(|i| mat.get([row, i]).unwrap().powi(2)) +// .sum::() +// .sqrt(); +// for i in 0..mat.shape()[1] { +// *mat.get_mut([row, i]).unwrap() /= norm; +// } +// for r in row + 1..mat.shape()[0] { +// let dot = (0..mat.shape()[1]) +// .map(|i| *mat.get([row, i]).unwrap() * *mat.get([r, i]).unwrap()) +// .sum::(); +// for i in 0..mat.shape()[1] { +// let sub = dot * *mat.get([row, i]).unwrap(); +// *mat.get_mut([r, i]).unwrap() -= sub; +// } +// } +// } +// } + #[cfg(test)] mod test { use super::*; diff --git a/src/polynomials.rs b/src/polynomials.rs index a56686d..26d49b7 100644 --- a/src/polynomials.rs +++ b/src/polynomials.rs @@ -5,7 +5,7 @@ pub use legendre::{shape as legendre_shape, tabulate as tabulate_legendre_polyno use crate::reference_cell; use crate::types::ReferenceCellType; -/// The number of polynomials +/// Return the number of polynomial terms of the form $x^iy^jz^k$ for a Lagrange type element. pub fn polynomial_count(cell_type: ReferenceCellType, degree: usize) -> usize { match cell_type { ReferenceCellType::Interval => degree + 1, @@ -19,12 +19,12 @@ pub fn polynomial_count(cell_type: ReferenceCellType, degree: usize) -> usize { } } -/// The total number of partial derivatives up to a give degree -pub fn derivative_count(cell_type: ReferenceCellType, derivatives: usize) -> usize { +/// Return the total number of partial derivatives up to a given degree. +pub fn derivative_count(cell_type: ReferenceCellType, degree: usize) -> usize { let mut num = 1; let mut denom = 1; for i in 0..reference_cell::dim(cell_type) { - num *= derivatives + i + 1; + num *= degree + i + 1; denom *= i + 1; } num / denom diff --git a/src/quadrature.rs b/src/quadrature.rs index 26cf9e2..75c304c 100644 --- a/src/quadrature.rs +++ b/src/quadrature.rs @@ -50,7 +50,7 @@ pub fn available_simplex_rules(cell: ReferenceCellType) -> Vec { } } -/// Get the points and weights for a Gauss-Jacobi quadrature rule +/// Get the points and weights for a Gauss-Jacobi quadrature rule of order `m` on the given cell type. pub fn gauss_jacobi_rule( celltype: ReferenceCellType, m: usize, @@ -66,7 +66,7 @@ pub fn gauss_jacobi_rule( } } -/// Get the number of quadrature points for a Gauss-Jacobi rule +/// Get the number of quadrature points for a Gauss-Jacobi rule of order `m` on a given cell type. pub fn gauss_jacobi_npoints(celltype: ReferenceCellType, m: usize) -> usize { let np = (m + 2) / 2; match celltype { diff --git a/src/traits.rs b/src/traits.rs index 76c9bd3..9512659 100644 --- a/src/traits.rs +++ b/src/traits.rs @@ -4,31 +4,87 @@ use rlst::{Array, MutableArrayImpl, RlstScalar, ValueArrayImpl}; use std::fmt::Debug; use std::hash::Hash; +/// This trait provides the definition of a finite element. +/// +/// A finite element in the Ciarlet definition can be described as a triplet +/// $(\mathcal{R}, \mathcal{V}, \mathcal{L})$. Here, $\mathcal{R}\in\mathbb{R}^d$ is the reference cell, +/// $\mathcal{V}$ is a finite dimensional function space on $\mathcal{R}$ of dimension $n$, usually polynomials, +/// $\mathcal{L} = \{\ell_0,\dots, \ell_{n-1}\}$ is a basis of the dual space $\mathcal{V}^* = \set{f:\mathcal{V} -> \mathbb{R}: f\text{ is linear}}$. +/// The basis functions $\phi_0,\dots, \phi_{n-1}$ of the finite element space are defined by +/// $$\ell_i(\phi_j) = \begin{cases}1 &\text{if }~i = j \newline +/// 0 &\text{otherwise}\end{cases}. +/// $$ +/// +/// ## Example +/// +/// The order 1 [Lagrange space](https://defelement.org/elements/lagrange.html) on a triangle is +/// defined by: +/// - $\mathcal{R}$ is a triangle with vertices $(0, 0)$, $(1, 0)$, $(0, 1)$. +/// - $\mathcal{V} = \text{span}\set{1, x, y}$. +/// - $\mathcal{L} = \set{\ell_0, \ell_0, \ell_1}$ with $\ell_j$ the pointwise evaluation at vertex $v_j$. +/// +/// This gives the basis functions $\phi_0(x, y) = 1 - x - y$, $\phi_1(x, y) = x$, $\phi_2(x, y) = y$. +/// +/// ## References +/// +/// - [https://defelement.org/ciarlet.html](https://defelement.org/ciarlet.html) pub trait FiniteElement { - //! A finite element defined on a reference cell /// The scalar type type T: RlstScalar; + /// Cell type + /// + /// The cell type is typically defined through [ReferenceCellType](crate::types::ReferenceCellType). type CellType: Debug + PartialEq + Eq + Clone + Copy + Hash; + /// Transformation type + /// + /// The Transformation type specifies possible transformations of the dofs on the reference element. + /// In most cases these will be rotations and reflections as defined in [Transformation](crate::types::Transformation). type TransformationType: Debug + PartialEq + Eq + Clone + Copy + Hash; - /// The reference cell type + /// The reference cell type, e.g. one of `Point`, `Interval`, `Triangle`, etc. fn cell_type(&self) -> Self::CellType; - /// The highest degree polynomial in the element's polynomial set - fn embedded_superdegree(&self) -> usize; + /// The smallest degree Lagrange space that contains all possible polynomials of the finite element's polynomial space. + /// + /// Details on the definition of the degree of Lagrange spaces of finite elements are + /// given [here](https://defelement.org/ciarlet.html#The+degree+of+a+finite+element). + fn lagrange_superdegree(&self) -> usize; - /// The number of basis functions + /// The number of basis functions. fn dim(&self) -> usize; - /// The value shape + /// The shape of the values returned by functions in $\mathcal{V}$. fn value_shape(&self) -> &[usize]; - /// The value size + /// The number of values returned. + /// + /// If e.g. `value_shape = [3, 4]` then `value_size = 3 x 4 = 12`. fn value_size(&self) -> usize; /// Tabulate the values of the basis functions and their derivatives at a set of points + /// + /// - `points` is a two-dimensional array where each column contains the coordinates of one point. + /// - `nderivs` is the desired number of derivatives (0 for no derivatives, 1 for all first order derivatives, etc.). + /// - `data` is the 4-dimensional output array. The first dimension is the number of derivatives, + /// the second dimension is the number of evaluation points, the third dimension is the number `n` of + /// basis functions on the element and the last dimension is the value size of the basis function output. + /// For example, `data[3, 2, 1, 0]` returns the 0th value of the third derivative on the point with index 2 for the + /// basis function with index 1. + /// + /// ## Remark + /// + /// Let $d^{i + k} = dx^{i}dy^{j}$ be a derivative with respect to $x$, $y$ in two dimensions and + /// $d^{i + k + j} = dx^{i}dy^{j}dz^{k}$ be a derivative with respect to $x$, $y$, and $z$ in three dimensions. + /// Then the corresponding index $\ell$ in the first dimension of the `data` array is computed as follows. + /// - Triangle: $\ell = (i + j + 1) * (i + j) / 2 + j$ + /// - Quadrilateral: $\ell = i * (n + 1) + j$ + /// - Tetrahedron: $\ell = (i + j + k) * (i + j + k + 1) * (i + j + k + 2) / 6 + (j + k) * (j + k + 1) / 2 + k$ + /// - Hexahedron $\ell = i * (n + 1) * (n + 1) + j * (n + 1) + k$. + /// + /// For the quadrilateral and hexahedron the parameter $n$ denotes the degree of the Lagrange space. + /// fn tabulate< Array2Impl: ValueArrayImpl<::Real, 2>, Array4MutImpl: MutableArrayImpl, @@ -39,16 +95,49 @@ pub trait FiniteElement { data: &mut Array, ); - /// The DOFs that are associated with a subentity of the reference cell + /// Return the dof indices that are associated with the subentity with index `entity_number` and dimension `entity_dim`. + /// + /// - For `entity_dim = 0` this returns the dof associated with the corresponding point. + /// - For `entity_dim = 1` this returns the dofs associated with the corresponding edge. + /// - For `entity_dim = 2` this returns the dofs associated with the corresponding face. + /// + /// Note that this does not return dofs on the boundary of an entity, that means e.g. + /// for an edge the dofs associated with the two vertices at the boundary of the edge are not returned. + /// To return also the boundary dofs use [FiniteElement::entity_closure_dofs]. fn entity_dofs(&self, entity_dim: usize, entity_number: usize) -> Option<&[usize]>; - /// The DOFs that are associated with a closure of a subentity of the reference cell + /// The DOFs that are associated with a closure of a subentity of the reference cell. + /// + /// This method is similar to [FiniteElement::entity_dofs]. But it returns additionally the dofs + /// associated with the boundary of an entity, e.g. for an edge it returns also the dofs associated + /// with the boundary vertices of they exist. fn entity_closure_dofs(&self, entity_dim: usize, entity_number: usize) -> Option<&[usize]>; - /// Get the required shape for a tabulation array + /// Get the required shape for a tabulation array. fn tabulate_array_shape(&self, nderivs: usize, npoints: usize) -> [usize; 4]; - /// Push function values forward to a physical cell + /// Push function values forward to a physical cell. + /// + /// Usually this is just an identity map. But for certain element types function values + /// on the reference cell differ from those on the physical cell, e.g. in the case of a + /// Piola transform. This method implements the corresponding transformation or an identity + /// map if no transformation is required. + /// + /// - `reference_values`: The values on the reference cell. + /// - `nderivs`: The number (degree) of derivatives. + /// - `jacobians:` A three-dimensional array of jacobians of the map from reference to physical cell. + /// The first dimension is the reference point. The second dimension is the geometric dimension of the physical space and + /// the third dimension is the topological dimension of the reference element, e.g. + /// for the map of 5 points from the reference triangle to a physical surface triangle embedded in 3d space the dimension + /// of `jacobians` is `[5, 3, 2]`. + /// - `jacobian_determinants`: Let $J$ be the Jacobian from the map of the reference to the physical element at a given point. + /// The corresponding Jacobian determinant is given as $d = \sqrt{\det(J^TJ)}$. `jacobian_determinants[j]` stores the Jacobian + /// determinant at position `j`. + /// - `inverse_jacobians`: A three dimensional array that stores the inverse Jacobian for the point with index j at position + /// `inverse_jacobians[j, :, :]`. The first dimension of `inverse_jacobians` is the point index. The second dimension + /// is the topological dimension, and the third dimension is the geometric dimension. If the Jacobian is rectangular then the + /// inverse Jacobian is the pseudo-inverse of the Jacobian such that $J^\dagger J = I$. + /// - `physical_values`: The output array of the push operation. Its required shape can be queried with [FiniteElement::physical_value_shape]. fn push_forward< Array3RealImpl: ValueArrayImpl<::Real, 3>, Array4Impl: ValueArrayImpl, @@ -63,7 +152,9 @@ pub trait FiniteElement { physical_values: &mut Array, ); - /// Pull function values back to the reference cell + /// Pull function values back to the reference cell. + /// + /// This is the inverse operation to [FiniteElement::push_forward]. fn pull_back< Array3RealImpl: ValueArrayImpl<::Real, 3>, Array4Impl: ValueArrayImpl, @@ -90,20 +181,20 @@ pub trait FiniteElement { vs } - /// The DOF transformation for a sub-entity + /// The DOF transformation for a sub-entity. fn dof_transformation( &self, entity: Self::CellType, transformation: Self::TransformationType, ) -> Option<&DofTransformation>; - /// Apply permutation parts of DOF transformations + /// Apply permutation parts of DOF transformations. fn apply_dof_permutations(&self, data: &mut [T], cell_orientation: i32); - /// Apply non-permutation parts of DOF transformations + /// Apply non-permutation parts of DOF transformations. fn apply_dof_transformations(&self, data: &mut [Self::T], cell_orientation: i32); - /// Apply DOF transformations + /// Apply DOF transformations. fn apply_dof_permutations_and_transformations( &self, data: &mut [Self::T], @@ -114,8 +205,8 @@ pub trait FiniteElement { } } +/// A factory that can create elements belonging to a given element family. pub trait ElementFamily { - //! A family of finite elements /// The scalar type type T: RlstScalar; /// Cell type @@ -123,28 +214,13 @@ pub trait ElementFamily { /// The finite element type type FiniteElement: FiniteElement + 'static; - /// Get an elenent for a cell type + /// Create an element for the given 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; -} - +/// A map from the reference cell to physical cells. pub trait Map { - //! A map from the reference cell to physical cells - - /// Push function values forward to a physical cell + /// Push values from the reference cell to physical cells. fn push_forward< T: RlstScalar, Array3RealImpl: ValueArrayImpl, @@ -160,7 +236,7 @@ pub trait Map { physical_values: &mut Array, ); - /// Pull function values back to the reference cell + /// Pull values back to the reference cell. fn pull_back< T: RlstScalar, Array3RealImpl: ValueArrayImpl, @@ -176,6 +252,6 @@ pub trait Map { reference_values: &mut Array, ); - /// The value shape on a physical cell + /// The value shape on a physical cell. fn physical_value_shape(&self, gdim: usize) -> Vec; }