diff --git a/src/bindings.rs b/src/bindings.rs index 1dd670e..fb6110c 100644 --- a/src/bindings.rs +++ b/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( + pub fn ciarlet_element_physical_value_size( 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( + pub fn ciarlet_element_physical_value_rank( 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( + pub fn ciarlet_element_physical_value_shape( 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>( + pub fn ciarlet_element_push_forward>( 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>( + pub fn ciarlet_element_pull_back>( 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(element: &E) -> usize { + pub fn ciarlet_element_embedded_superdegree(element: &E) -> usize { element.lagrange_superdegree() } diff --git a/src/ciarlet.rs b/src/ciarlet.rs index 3f269f4..fcb702f 100644 --- a/src/ciarlet.rs +++ b/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,13 +621,15 @@ where &self.interpolation_weights } } + impl FiniteElement for CiarletElement { 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 } @@ -635,12 +637,11 @@ impl FiniteElement for CiarletElement { 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<::Real, 2>, Array4MutImpl: MutableArrayImpl, @@ -681,6 +682,15 @@ impl FiniteElement for CiarletElement { } } } + + 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]) @@ -688,6 +698,7 @@ impl FiniteElement for CiarletElement { 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]) @@ -695,13 +706,15 @@ impl FiniteElement for CiarletElement { 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 MappedFiniteElement for CiarletElement { + type TransformationType = Transformation; + + fn lagrange_superdegree(&self) -> usize { + self.embedded_superdegree } + fn push_forward< Array3RealImpl: ValueArrayImpl<::Real, 3>, Array4Impl: ValueArrayImpl, @@ -724,6 +737,7 @@ impl FiniteElement for CiarletElement { physical_values, ) } + fn pull_back< Array3RealImpl: ValueArrayImpl<::Real, 3>, Array4Impl: ValueArrayImpl, @@ -746,9 +760,11 @@ impl FiniteElement for CiarletElement { reference_values, ) } + fn physical_value_shape(&self, gdim: usize) -> Vec { self.map.physical_value_shape(gdim) } + fn dof_transformation( &self, entity: ReferenceCellType, diff --git a/src/traits.rs b/src/traits.rs index 86f454d..7caae3d 100644 --- a/src/traits.rs +++ b/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, ); + /// 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<::Real, 3>, Array4Impl: ValueArrayImpl, @@ -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<::Real, 3>, Array4Impl: ValueArrayImpl,