Updates to the documentation

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"

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());
 }

katex-header.html

@@ -0,0 +1,19 @@
+<link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/katex@0.16.25/dist/katex.min.css" integrity="sha384-WcoG4HRXMzYzfCgiyfrySxx90XSl2rxY5mnVY5TwtWE6KLrArNKn0T/mOgNL0Mmi" crossorigin="anonymous">
+<script defer src="https://cdn.jsdelivr.net/npm/katex@0.16.25/dist/katex.min.js" integrity="sha384-J+9dG2KMoiR9hqcFao0IBLwxt6zpcyN68IgwzsCSkbreXUjmNVRhPFTssqdSGjwQ" crossorigin="anonymous"></script>
+<script defer src="https://cdn.jsdelivr.net/npm/katex@0.16.25/dist/contrib/auto-render.min.js" integrity="sha384-hCXGrW6PitJEwbkoStFjeJxv+fSOOQKOPbJxSfM6G5sWZjAyWhXiTIIAmQqnlLlh" crossorigin="anonymous"></script>
+<script>
+    document.addEventListener("DOMContentLoaded", function() {
+        renderMathInElement(document.body, {
+          // customised options
+          // • auto-render specific keys, e.g.:
+          delimiters: [
+              {left: '$$', right: '$$', display: true},
+              {left: '$', right: '$', display: false},
+              {left: '\\(', right: '\\)', display: false},
+              {left: '\\[', right: '\\]', display: true}
+          ],
+          // • rendering keys, e.g.:
+          throwOnError : false
+        });
+    });
+</script>

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<E: FiniteElement>(element: &E) -> usize {
-        element.embedded_superdegree()
+        element.lagrange_superdegree()
     }
 
     #[concretise_types(

src/ciarlet.rs

@@ -72,7 +72,13 @@ impl<T: RlstScalar, M: Map> CiarletElement<T, M>
 where
     DynArray<T, 2>: Inverse<Output = DynArray<T, 2>>,
 {
-    /// 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<T: RlstScalar, M: Map> FiniteElement for CiarletElement<T, M> {
     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 {

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<T: RlstScalar + Getrf + Getri>(
     cell_type: ReferenceCellType,
     degree: usize,
@@ -268,15 +268,17 @@ pub fn create<T: RlstScalar + Getrf + Getri>(
     )
 }
 
-/// Lagrange element family
+/// Lagrange element family.
+///
+/// A factory structure to create new Lagrange elements.
 pub struct LagrangeElementFamily<T: RlstScalar + Getrf + Getri> {
     degree: usize,
     continuity: Continuity,
     _t: PhantomData<T>,
 }
 
 impl<T: RlstScalar + Getrf + Getri> LagrangeElementFamily<T> {
-    /// Create new family
+    /// Create new family with given `degree` and `continuity`.
     pub fn new(degree: usize, continuity: Continuity) -> Self {
         Self {
             degree,

src/ciarlet/nedelec.rs

@@ -640,14 +640,16 @@ pub fn create<T: RlstScalar + Getrf + Getri>(
 }
 
 /// Nedelec (first kind) element family
+///
+/// A factory structure to create new Nedelec elements.
 pub struct NedelecFirstKindElementFamily<T: RlstScalar + Getrf + Getri> {
     degree: usize,
     continuity: Continuity,
     _t: PhantomData<T>,
 }
 
 impl<T: RlstScalar + Getrf + Getri> NedelecFirstKindElementFamily<T> {
-    /// Create new family
+    /// Create new family with given `degree` and `continuity`.
     pub fn new(degree: usize, continuity: Continuity) -> Self {
         Self {
             degree,

src/ciarlet/raviart_thomas.rs

@@ -483,14 +483,16 @@ pub fn create<T: RlstScalar + Getrf + Getri>(
 }
 
 /// Raviart-Thomas element family
+///
+/// A factory structure to create new Raviart-Thomas elements.
 pub struct RaviartThomasElementFamily<T: RlstScalar + Getrf + Getri> {
     degree: usize,
     continuity: Continuity,
     _t: PhantomData<T>,
 }
 
 impl<T: RlstScalar + Getrf + Getri> RaviartThomasElementFamily<T> {
-    /// Create new family
+    /// Create new family with given `degree` and `continuity`.
     pub fn new(degree: usize, continuity: Continuity) -> Self {
         Self {
             degree,

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::<f64>::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::<f64, 4>::from_shape(element.tabulate_array_shape(0, 1));
+//! // Create array containing the point [1/2, 1/2]
+//!    let mut points = DynArray::<f64, 2>::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)]
 

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 {

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<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;
-            }
-        }
-    }
-}
-
-/// Orthogonalise the rows of a matrix, starting with the row numbered `start`
-pub fn orthogonalise_3<T: RlstScalar, Array3MutImpl: MutableArrayImpl<T, 3>>(
+pub(crate) fn orthogonalise_3<T: RlstScalar, Array3MutImpl: MutableArrayImpl<T, 3>>(
     mat: &mut Array<Array3MutImpl, 3>,
     start: usize,
 ) {
@@ -186,6 +161,31 @@ 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::*;

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

src/quadrature.rs

@@ -50,7 +50,7 @@ pub fn available_simplex_rules(cell: ReferenceCellType) -> Vec<usize> {
     }
 }
 
-/// 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 {