diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index 67d25f5..ce518ef 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -27,7 +27,7 @@ jobs: mpi: "mpich" - uses: actions/checkout@v4 - name: Build Rust docs - run: cargo +nightly doc --no-deps -Zunstable-options -Zrustdoc-scrape-examples --all-features + run: RUSTDOCFLAGS="--html-in-header katex-header.html" cargo +nightly doc --no-deps -Zunstable-options -Zrustdoc-scrape-examples --all-features - name: Make index page run: echo "" > target/doc/index.html - name: Set file permissions diff --git a/Cargo.toml b/Cargo.toml index eca6102..d5aa7e8 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -44,10 +44,6 @@ 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.*" diff --git a/src/bindings.rs b/src/bindings.rs index cae544d..8d498ab 100644 --- a/src/bindings.rs +++ b/src/bindings.rs @@ -1,4 +1,4 @@ -//! Binding for C +//! Bindings for C. #![allow(missing_docs)] #![allow(clippy::missing_safety_doc)] diff --git a/src/ciarlet.rs b/src/ciarlet.rs index 958ac20..0a73361 100644 --- a/src/ciarlet.rs +++ b/src/ciarlet.rs @@ -1,4 +1,30 @@ -//! Finite element definitions +//! Ciarlet finite elements. +//! +//! In the Ciarlet definition, a finite element is defined to be a triple +//! $(\mathcal{R}, \mathcal{V}, \mathcal{L})$, where: +//! +//! - $\mathcal{R}\subset\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}}$, with +// each functional associated with a subentity of the reference cell $\mathcal{R}$. +//! +//! 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$, and each functional associated with the relevant vertex. +//! +//! 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) extern crate blas_src; extern crate lapack_src; @@ -74,8 +100,8 @@ where { /// 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. + /// The majority of users will not wish to use this directly, and should insteal 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. diff --git a/src/ciarlet/lagrange.rs b/src/ciarlet/lagrange.rs index 87890cc..af8dfcc 100644 --- a/src/ciarlet/lagrange.rs +++ b/src/ciarlet/lagrange.rs @@ -1,4 +1,4 @@ -//! Lagrange elements +//! Lagrange elements. use super::CiarletElement; use crate::map::IdentityMap; @@ -270,7 +270,8 @@ pub fn create( /// Lagrange element family. /// -/// A factory structure to create new Lagrange elements. +/// A family of Lagrange elements on multiple cell types with appropriate +/// continuity across different cell types. pub struct LagrangeElementFamily { degree: usize, continuity: Continuity, diff --git a/src/ciarlet/nedelec.rs b/src/ciarlet/nedelec.rs index c03e96d..23bd707 100644 --- a/src/ciarlet/nedelec.rs +++ b/src/ciarlet/nedelec.rs @@ -1,4 +1,4 @@ -//! Nedelec elements +//! Nedelec elements. use super::CiarletElement; use crate::map::CovariantPiolaMap; @@ -641,7 +641,8 @@ pub fn create( /// Nedelec (first kind) element family /// -/// A factory structure to create new Nedelec elements. +/// A family of Nedelec elements on multiple cell types with appropriate +/// continuity across different cell types. pub struct NedelecFirstKindElementFamily { degree: usize, continuity: Continuity, diff --git a/src/ciarlet/raviart_thomas.rs b/src/ciarlet/raviart_thomas.rs index c437f63..a56b820 100644 --- a/src/ciarlet/raviart_thomas.rs +++ b/src/ciarlet/raviart_thomas.rs @@ -1,4 +1,4 @@ -//! Raviart-Thomas elements +//! Raviart-Thomas elements. use super::CiarletElement; use crate::map::ContravariantPiolaMap; @@ -484,7 +484,8 @@ pub fn create( /// Raviart-Thomas element family /// -/// A factory structure to create new Raviart-Thomas elements. +/// A family of Raviart-Thomas elements on multiple cell types with appropriate +/// continuity across different cell types. pub struct RaviartThomasElementFamily { degree: usize, continuity: Continuity, diff --git a/src/lib.rs b/src/lib.rs index cf5c53b..8575a52 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -33,7 +33,7 @@ //! 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() +//! basis_values.data() //! ); //! ``` #![cfg_attr(feature = "strict", deny(warnings), deny(unused_crate_dependencies))] diff --git a/src/map.rs b/src/map.rs index f4af76c..12aa282 100644 --- a/src/map.rs +++ b/src/map.rs @@ -1,10 +1,10 @@ -//! Maps from the reference cell to/from physical cells +//! Maps from the reference cell to/from physical cells. use crate::traits::Map; use rlst::{Array, MutableArrayImpl, RlstScalar, ValueArrayImpl}; /// Identity map /// -/// An identity map pushes values from the reference to the physical +/// An identity map pushes values from the reference to the physical /// cell without modifying them. pub struct IdentityMap {} @@ -62,7 +62,7 @@ impl Map for IdentityMap { /// $$ /// J^{-T}\Phi\circ F^{-1} /// $$ -/// The coviariant Piola map preserves tangential continuity. If $J$ +/// The covariant Piola map preserves tangential continuity. If $J$ /// is a rectangular matrix then the pseudo-inverse is used instead of /// the inverse. pub struct CovariantPiolaMap {} @@ -157,9 +157,11 @@ impl Map for CovariantPiolaMap { /// 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} +/// \frac{1}{\det{J}}J\Phi\circ F^{-1} /// $$ -/// The contravariant Piola map preserves normal continuity. +/// The contravariant Piola map preserves normal continuity. If $J$ +/// is a rectangular matrix then $\sqrt{\det{J^TJ}}$ is used instead +/// of $\det{J}$. pub struct ContravariantPiolaMap {} impl Map for ContravariantPiolaMap { diff --git a/src/math.rs b/src/math.rs index 5aa77d2..8614125 100644 --- a/src/math.rs +++ b/src/math.rs @@ -1,4 +1,4 @@ -//! Mathematical functions +//! Mathematical functions. use rlst::{Array, MutableArrayImpl, RlstScalar, ValueArrayImpl}; /// Orthogonalise the rows of a matrix, starting with the row numbered `start` @@ -161,31 +161,6 @@ 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/orientation.rs b/src/orientation.rs index 076c03d..c9a2c99 100644 --- a/src/orientation.rs +++ b/src/orientation.rs @@ -1,4 +1,4 @@ -//! Cell orientation +//! Cell orientation. use crate::{reference_cell, types::ReferenceCellType}; use itertools::izip; diff --git a/src/polynomials.rs b/src/polynomials.rs index 26d49b7..e2d3db7 100644 --- a/src/polynomials.rs +++ b/src/polynomials.rs @@ -1,11 +1,28 @@ -//! Orthonormal polynomials +//! Orthonormal polynomials. +//! +//! The orthonormal polynomials in ndelement span the degree $k$ natural polynomial (or rationomial +//! for a pyramid) space on a cell. As given at , +//! these natural spaces are defined for an interval, triangle, quadrilateral, tetrahedron, +//! hexahedron, triangular prism and square-based pyramid by +//! +//! - $\mathbb{P}^{\text{interval}}_k=\operatorname{span}\left\{x^{p_0}\,\middle|\,p_0\in\mathbb{N},\,p_0\leqslant k\right\},$ +//! - $\mathbb{P}^{\text{triangle}}_k=\operatorname{span}\left\{x^{p_0}y^{p_1}\,\middle|\,p_0,p_1\in\mathbb{N}_0,\,p_0+p_1\leqslant k\right\},$ +//! - $\mathbb{P}^{\text{quadrilateral}}_k=\operatorname{span}\left\{x^{p_0}y^{p_1}\,\middle|\,p_0,p_1\in\mathbb{N}_0,\,p_0\leqslant k,\,p_1\leqslant k\right\},$ +//! - $\mathbb{P}^{\text{tetrahedron}}_k=\operatorname{span}\left\{x^{p_0}y^{p_1}z^{p_2}\,\middle|\,p_0,p_1,p_2\in\mathbb{N}_0,\,p_0+p_1+p_2\leqslant k\right\},$ +//! - $\mathbb{P}^{\text{hexahedron}}_k=\operatorname{span}\left\{x^{p_0}y^{p_1}z^{p_2}\,\middle|\,p_0,p_1,p_2\in\mathbb{N}_0,\,p_0\leqslant k,\,p_1\leqslant k,\,p_2\leqslant k\right\},$ +//! - $\mathbb{P}^{\text{prism}}_k=\operatorname{span}\left\{x^{p_0}y^{p_1}z^{p_2}\,\middle|\,p_0,p_1,p_2\in\mathbb{N}_0,\,p_0+p_1\leqslant k,\,p_2\leqslant k\right\},$ +//! - $\mathbb{P}^{\text{pyramid}}_k=\operatorname{span}\left\{\frac{x^{p_0}y^{p_1}z^{p_2}}{(1-z)^{p_0+p_1}}\,\middle|\,p_0,p_1,p_2\in\mathbb{N}_0,\,p_0\leqslant k,\,p_1\leqslant k,\,p_2\leqslant k\right\}.$ +//! +//! Note that for non-pyramid cells, these coincide with the polynomial space for the degree $k$ +//! Lagrange element on the cell. + mod legendre; pub use legendre::{shape as legendre_shape, tabulate as tabulate_legendre_polynomials}; use crate::reference_cell; use crate::types::ReferenceCellType; -/// Return the number of polynomial terms of the form $x^iy^jz^k$ for a Lagrange type element. +/// Return the number of polynomials in the natural polynomial set for a given cell type and degree. pub fn polynomial_count(cell_type: ReferenceCellType, degree: usize) -> usize { match cell_type { ReferenceCellType::Interval => degree + 1, diff --git a/src/polynomials/legendre.rs b/src/polynomials/legendre.rs index 96e2ead..ae4449b 100644 --- a/src/polynomials/legendre.rs +++ b/src/polynomials/legendre.rs @@ -1,4 +1,4 @@ -//! Orthonormal polynomials +//! Orthonormal polynomials. //! //! Adapted from the C++ code by Chris Richardson and Matthew Scroggs //! diff --git a/src/quadrature.rs b/src/quadrature.rs index 75c304c..d931112 100644 --- a/src/quadrature.rs +++ b/src/quadrature.rs @@ -1,4 +1,4 @@ -//! Quadrature +//! Quadrature. use crate::types::ReferenceCellType; use bempp_quadrature::{ gauss_jacobi, simplex_rules, @@ -50,7 +50,7 @@ pub fn available_simplex_rules(cell: ReferenceCellType) -> Vec { } } -/// Get the points and weights for a Gauss-Jacobi quadrature rule of order `m` on the given cell type. +/// Get the points and weights for a Gauss-Jacobi quadrature rule of degree `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 of order `m` on a given cell type. +/// Get the number of quadrature points for a Gauss-Jacobi rule of degree `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/reference_cell.rs b/src/reference_cell.rs index e42a856..905a719 100644 --- a/src/reference_cell.rs +++ b/src/reference_cell.rs @@ -1,4 +1,13 @@ -//! Cell definitions +//! Cell definitions. +//! +//! The reference cells used by ndelement are defined by: +//! - The reference interval is \(\left\{x\,\middle|\,0\leqslant x\leqslant 1\right\}\). +//! - The reference triangle is \(\left\{(x,y)\,\middle|\,0\leqslant x,\,0\leqslant y,\,x+y\leqslant1\right\}\). +//! - The reference quadrilateral is \(\left\{(x,y)\,\middle|\,0\leqslant x\leqslant1,\,0\leqslant y\leqslant1\right\}\). +//! - The reference tetrahedron is \(\left\{(x,y,z)\,\middle|\,0\leqslant x,\,0\leqslant y,\,0\leqslant z,\,x+y+z\leqslant1\right\}\). +//! - The reference hexahedron is \(\left\{(x,y,z)\,\middle|\,0\leqslant x\leqslant1,\,0\leqslant y\leqslant1,\,0\leqslant z\leqslant1\right\}\). +//! - The reference prism is \(\left\{(x,y,z)\,\middle|\,0\leqslant x,\,0\leqslant y,\,x+y\leqslant1,\,0\leqslant z\leqslant1\right\}\). +//! - The reference pyramid is \(\left\{(x,y,z)\,\middle|\,0\leqslant x,\,0\leqslant y,\,0\leqslant z\leqslant1,\,x+z\leqslant1,\,y+z\leqslant1\right\}\). use crate::types::ReferenceCellType; use rlst::RlstScalar; diff --git a/src/traits.rs b/src/traits.rs index 88be8c2..86f454d 100644 --- a/src/traits.rs +++ b/src/traits.rs @@ -1,33 +1,10 @@ -//! Traits +//! Traits. use crate::types::DofTransformation; 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 { /// The scalar type type T: RlstScalar; @@ -43,7 +20,7 @@ pub trait FiniteElement { /// 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, e.g. one of `Point`, `Interval`, `Triangle`, etc. + /// The reference cell type, eg one of `Point`, `Interval`, `Triangle`, etc. fn cell_type(&self) -> Self::CellType; /// The smallest degree Lagrange space that contains all possible polynomials of the finite element's polynomial space. @@ -62,32 +39,32 @@ pub trait FiniteElement { /// The number of values returned. /// - /// This is the product of the elements of `value_shape`, e.g. if `value_shape = [3, 4]` then `value_size = 3 x 4 = 12`. - /// If `value_shape` returns an empty array (i.e. the shape functions are scalar) the convention is used that the product of the elements of an empty array is `1`. + /// If (for example) `value_shape` is `[3, 4]` then `value_size` is $3\times4 = 12$. + /// If `value_shape` returns an empty array (ie the shape functions are scalar) the + // convention is used that the product of the elements of an empty array is 1. 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 + /// - `data` is the 4-dimensional output array. The first dimension is the total number of partial derivatives, + /// the second dimension is the number of evaluation points, the third dimension is the number of + /// basis functions of the element, and the fourth dimension is the value size of the basis function output. + /// For example, `data[3, 2, 1, 0]` returns the 0th value of the third partial 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. + /// - Triangle: $l = (i + j + 1) * (i + j) / 2 + j$ + /// - Quadrilateral: $l = i * (n + 1) + j$ + /// - Tetrahedron: $l = (i + j + k) * (i + j + k + 1) * (i + j + k + 2) / 6 + (j + k) * (j + k + 1) / 2 + k$ + /// - Hexahedron $l = i * (n + 1) * (n + 1) + j * (n + 1) + k$. /// + /// For the quadrilaterals and hexahedra, the parameter $n$ denotes the Lagrange superdegree. fn tabulate< Array2Impl: ValueArrayImpl<::Real, 2>, Array4MutImpl: MutableArrayImpl, @@ -100,20 +77,20 @@ pub trait FiniteElement { /// 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 = 0` this returns the degrees of freedom (dofs) 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. + /// Note that this does not return dofs on the boundary of an entity, that means that (for example) /// 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. /// - /// 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. + /// This method is similar to [FiniteElement::entity_dofs], but it includes the dofs + /// associated with the boundary of an entity. For an edge (for example) it returns the dofs associated + /// with the vertices at the boundary of the edge (as well as the dofs associated with the edge itself). fn entity_closure_dofs(&self, entity_dim: usize, entity_number: usize) -> Option<&[usize]>; /// Get the required shape for a tabulation array. @@ -121,26 +98,26 @@ pub trait FiniteElement { /// 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. + /// For Lagrange elements, this is an identity map. For many other element types, the function values + /// on the reference cell differ from those on the physical cell: for example Nedlec elements use a covariant + /// Piola transform. This method implements the appropriate transformation for the element. /// - /// - `reference_values`: The values on the reference cell. - /// - `nderivs`: The number (degree) of derivatives. + /// - `reference_values`: The values on the reference cell. The shape of this input is the same as the `data` input to the function + /// [[FiniteElement::tabulate]. + /// - `nderivs`: The number 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. + /// 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. For example, /// 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 + /// - `jacobian_determinants`: The determinant of the jacobian at each point. If the jacobian $J$ is not square, then the + /// determinant is computed using $d=\sqrt{\det(J^TJ)}$. + /// - `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]. + /// inverse Jacobian is the pseudo-inverse of the Jacobian, ie the matrix $J^\dagger$ such that $J^\dagger J = I$. + /// - `physical_values`: The output array of the push operation. This shape of this array is the same as the `reference_values` + /// input, with the [FiniteElement::physical_value_size] used instead of the reference value size. fn push_forward< Array3RealImpl: ValueArrayImpl<::Real, 3>, Array4Impl: ValueArrayImpl, @@ -208,7 +185,8 @@ pub trait FiniteElement { } } -/// A factory that can create elements belonging to a given element family. +/// A family of finite elements defined on various cell types, with appropriate continuity +/// between different cell types. pub trait ElementFamily { /// The scalar type type T: RlstScalar; diff --git a/src/types.rs b/src/types.rs index 50d83f6..e98710e 100644 --- a/src/types.rs +++ b/src/types.rs @@ -1,4 +1,4 @@ -//! Types +//! Types. #[cfg(feature = "mpi")] use mpi::traits::Equivalence; use rlst::{DynArray, RlstScalar};