Add MappedFiniteElement trait (#109)

* add MappedFiniteElement trait

* doc

Author: mscroggs

src/bindings.rs

@@ -254,7 +254,7 @@ pub mod ciarlet {
         ciarlet::CiarletElement,
         map::{ContravariantPiolaMap, CovariantPiolaMap, IdentityMap},
         reference_cell,
-        traits::{ElementFamily, FiniteElement, Map},
+        traits::{ElementFamily, FiniteElement, Map, MappedFiniteElement},
         types::{Continuity, ReferenceCellType},
     };
     use c_api_tools::{DType, DTypeIdentifier, cfuncs, concretise_types};
@@ -495,7 +495,7 @@ pub mod ciarlet {
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
         field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
     )]
-    pub fn ciarlet_element_physical_value_size<E: FiniteElement>(
+    pub fn ciarlet_element_physical_value_size<E: MappedFiniteElement>(
         element: &E,
         gdim: usize,
     ) -> usize {
@@ -507,7 +507,7 @@ pub mod ciarlet {
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
         field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
     )]
-    pub fn ciarlet_element_physical_value_rank<E: FiniteElement>(
+    pub fn ciarlet_element_physical_value_rank<E: MappedFiniteElement>(
         element: &E,
         gdim: usize,
     ) -> usize {
@@ -519,7 +519,7 @@ pub mod ciarlet {
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
         field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
     )]
-    pub fn ciarlet_element_physical_value_shape<E: FiniteElement>(
+    pub fn ciarlet_element_physical_value_shape<E: MappedFiniteElement>(
         element: &E,
         gdim: usize,
         shape: *mut usize,
@@ -537,7 +537,7 @@ pub mod ciarlet {
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
         field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
     )]
-    pub fn ciarlet_element_push_forward<E: FiniteElement<CellType = ReferenceCellType>>(
+    pub fn ciarlet_element_push_forward<E: MappedFiniteElement<CellType = ReferenceCellType>>(
         element: &E,
         npoints: usize,
         nfunctions: usize,
@@ -606,7 +606,7 @@ pub mod ciarlet {
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
         field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
     )]
-    pub fn ciarlet_element_pull_back<E: FiniteElement<CellType = ReferenceCellType>>(
+    pub fn ciarlet_element_pull_back<E: MappedFiniteElement<CellType = ReferenceCellType>>(
         element: &E,
         npoints: usize,
         nfunctions: usize,
@@ -685,7 +685,7 @@ pub mod ciarlet {
         gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]),
         field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"])
     )]
-    pub fn ciarlet_element_embedded_superdegree<E: FiniteElement>(element: &E) -> usize {
+    pub fn ciarlet_element_embedded_superdegree<E: MappedFiniteElement>(element: &E) -> usize {
         element.lagrange_superdegree()
     }
 

src/ciarlet.rs

@@ -32,7 +32,7 @@ extern crate lapack_src;
 use crate::math;
 use crate::polynomials::{legendre_shape, polynomial_count, tabulate_legendre_polynomials};
 use crate::reference_cell;
-use crate::traits::{FiniteElement, Map};
+use crate::traits::{FiniteElement, Map, MappedFiniteElement};
 use crate::types::{Continuity, DofTransformation, ReferenceCellType, Transformation};
 use itertools::izip;
 use num::{One, Zero};
@@ -621,26 +621,27 @@ where
         &self.interpolation_weights
     }
 }
+
 impl<T: RlstScalar, M: Map> FiniteElement for CiarletElement<T, M> {
     type CellType = ReferenceCellType;
-    type TransformationType = Transformation;
     type T = T;
+
     fn value_shape(&self) -> &[usize] {
         &self.value_shape
     }
+
     fn value_size(&self) -> usize {
         self.value_size
     }
 
     fn cell_type(&self) -> ReferenceCellType {
         self.cell_type
     }
-    fn lagrange_superdegree(&self) -> usize {
-        self.embedded_superdegree
-    }
+
     fn dim(&self) -> usize {
         self.dim
     }
+
     fn tabulate<
         Array2Impl: ValueArrayImpl<<Self::T as RlstScalar>::Real, 2>,
         Array4MutImpl: MutableArrayImpl<Self::T, 4>,
@@ -681,27 +682,39 @@ impl<T: RlstScalar, M: Map> FiniteElement for CiarletElement<T, M> {
             }
         }
     }
+
+    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]
+    }
+
     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 entity_closure_dofs(&self, entity_dim: usize, entity_number: usize) -> Option<&[usize]> {
         if entity_dim < 4 && entity_number < self.entity_closure_dofs[entity_dim].len() {
             Some(&self.entity_closure_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]
+}
+
+impl<T: RlstScalar, M: Map> MappedFiniteElement for CiarletElement<T, M> {
+    type TransformationType = Transformation;
+
+    fn lagrange_superdegree(&self) -> usize {
+        self.embedded_superdegree
     }
+
     fn push_forward<
         Array3RealImpl: ValueArrayImpl<<Self::T as RlstScalar>::Real, 3>,
         Array4Impl: ValueArrayImpl<Self::T, 4>,
@@ -724,6 +737,7 @@ impl<T: RlstScalar, M: Map> FiniteElement for CiarletElement<T, M> {
             physical_values,
         )
     }
+
     fn pull_back<
         Array3RealImpl: ValueArrayImpl<<Self::T as RlstScalar>::Real, 3>,
         Array4Impl: ValueArrayImpl<Self::T, 4>,
@@ -746,9 +760,11 @@ impl<T: RlstScalar, M: Map> FiniteElement for CiarletElement<T, M> {
             reference_values,
         )
     }
+
     fn physical_value_shape(&self, gdim: usize) -> Vec<usize> {
         self.map.physical_value_shape(gdim)
     }
+
     fn dof_transformation(
         &self,
         entity: ReferenceCellType,

src/traits.rs

@@ -4,7 +4,7 @@ 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.
 pub trait FiniteElement {
     /// The scalar type
     type T: RlstScalar;
@@ -14,21 +14,9 @@ pub trait FiniteElement {
     /// 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, 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.
-    ///
-    /// 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.
     fn dim(&self) -> usize;
 
@@ -75,6 +63,9 @@ pub trait FiniteElement {
         data: &mut Array<Array4MutImpl, 4>,
     );
 
+    /// Get the required shape for a tabulation array.
+    fn tabulate_array_shape(&self, nderivs: usize, npoints: usize) -> [usize; 4];
+
     /// 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 degrees of freedom (dofs) associated with the corresponding point.
@@ -92,9 +83,21 @@ pub trait FiniteElement {
     /// 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.
-    fn tabulate_array_shape(&self, nderivs: usize, npoints: usize) -> [usize; 4];
+/// A finite element that is mapped from a reference cell.
+pub trait MappedFiniteElement: FiniteElement {
+    /// 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 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;
 
     /// Push function values forward to a physical cell.
     ///
@@ -117,7 +120,7 @@ pub trait FiniteElement {
     ///   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, 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.
+    ///   input, with the [MappedFiniteElement::physical_value_size] used instead of the reference value size.
     fn push_forward<
         Array3RealImpl: ValueArrayImpl<<Self::T as RlstScalar>::Real, 3>,
         Array4Impl: ValueArrayImpl<Self::T, 4>,
@@ -134,7 +137,7 @@ pub trait FiniteElement {
 
     /// Pull function values back to the reference cell.
     ///
-    /// This is the inverse operation to [FiniteElement::push_forward].
+    /// This is the inverse operation to [MappedFiniteElement::push_forward].
     fn pull_back<
         Array3RealImpl: ValueArrayImpl<<Self::T as RlstScalar>::Real, 3>,
         Array4Impl: ValueArrayImpl<Self::T, 4>,