bempp · mscroggs · Nov 18, 2025 · Nov 18, 2025 · Nov 18, 2025 · Nov 18, 2025
diff --git 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 "<html><head><meta http-equiv='refresh' content='0; URL=ndelement'></head></html>" > target/doc/index.html
       - name: Set file permissions

diff --git 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
@@ -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
@@ -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
@@ -1,4 +1,4 @@
-//! Lagrange elements
+//! Lagrange elements.
 
 use super::CiarletElement;
 use crate::map::IdentityMap;
@@ -270,7 +270,8 @@ pub fn create<T: RlstScalar + Getrf + Getri>(
 
 /// 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<T: RlstScalar + Getrf + Getri> {
     degree: usize,
     continuity: Continuity,

diff --git 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<T: RlstScalar + Getrf + Getri>(
 
 /// 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<T: RlstScalar + Getrf + Getri> {
     degree: usize,
     continuity: Continuity,

diff --git 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<T: RlstScalar + Getrf + Getri>(
 
 /// 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<T: RlstScalar + Getrf + Getri> {
     degree: usize,
     continuity: Continuity,

diff --git 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
@@ -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
@@ -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<T: RlstScalar, Array2Impl: ValueArrayImpl<T, 2>>(
     }
 }
 
-// /// Orthogonalise the rows of a matrix, starting with the row numbered `start`
-// pub(crate) fn orthogonalise<T: RlstScalar, Array2MutImpl: MutableArrayImpl<T, 2>>(
-//     mat: &mut Array<Array2MutImpl, 2>,
-//     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::<T>()
-//             .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::<T>();
-//             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
@@ -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
@@ -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 <https://defelement.org/ciarlet.html#The+degree+of+a+finite+element>,
+//! 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
@@ -1,4 +1,4 @@
-//! Orthonormal polynomials
+//! Orthonormal polynomials.
 //!
 //! Adapted from the C++ code by Chris Richardson and Matthew Scroggs
 //! <https://github.com/FEniCS/basix/blob/main/cpp/basix/polyset.cpp>

diff --git 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<usize> {
     }
 }
 
-/// 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
@@ -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;