From fb26adb4d6acde0eb50bbc35396d6de95359a9c4 Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Tue, 25 Nov 2025 10:37:34 +0000 Subject: [PATCH 1/9] remove real from reference_cell --- src/reference_cell.rs | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/reference_cell.rs b/src/reference_cell.rs index 905a719..063e3b0 100644 --- a/src/reference_cell.rs +++ b/src/reference_cell.rs @@ -40,7 +40,7 @@ pub fn is_simplex(cell: ReferenceCellType) -> bool { } /// The vertices of the reference cell -pub fn vertices>(cell: ReferenceCellType) -> Vec> { +pub fn vertices(cell: ReferenceCellType) -> Vec> { let zero = T::from(0.0).unwrap(); let one = T::from(1.0).unwrap(); match cell { @@ -88,7 +88,7 @@ pub fn vertices>(cell: ReferenceCellType) -> Vec> } /// The midpoint of the cell -pub fn midpoint>(cell: ReferenceCellType) -> Vec { +pub fn midpoint(cell: ReferenceCellType) -> Vec { let half = T::from(0.5).unwrap(); let third = T::from(1.0).unwrap() / T::from(3.0).unwrap(); match cell { @@ -258,7 +258,7 @@ pub fn peaks(cell: ReferenceCellType) -> Vec> { } /// The normals to the facets of the reference cell. The length of each normal is the volume of the facet -pub fn facet_normals>(cell: ReferenceCellType) -> Vec> { +pub fn facet_normals(cell: ReferenceCellType) -> Vec> { let zero = T::from(0.0).unwrap(); let one = T::from(1.0).unwrap(); match cell { @@ -303,7 +303,7 @@ pub fn facet_normals>(cell: ReferenceCellType) -> Vec>(cell: ReferenceCellType) -> Vec> { +pub fn facet_unit_normals(cell: ReferenceCellType) -> Vec> { let mut normals = facet_normals::(cell); for normal in normals.iter_mut() { let size = normal.iter().map(|x| x.powi(2)).sum::().sqrt(); From 6ad2c3c5b88c5b3ae7c327e36ba4c133edc66103 Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Tue, 25 Nov 2025 10:44:23 +0000 Subject: [PATCH 2/9] default f64 --- src/ciarlet/lagrange.rs | 2 +- src/ciarlet/nedelec.rs | 2 +- src/ciarlet/raviart_thomas.rs | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/ciarlet/lagrange.rs b/src/ciarlet/lagrange.rs index af8dfcc..0879239 100644 --- a/src/ciarlet/lagrange.rs +++ b/src/ciarlet/lagrange.rs @@ -272,7 +272,7 @@ pub fn create( /// /// A family of Lagrange elements on multiple cell types with appropriate /// continuity across different cell types. -pub struct LagrangeElementFamily { +pub struct LagrangeElementFamily { degree: usize, continuity: Continuity, _t: PhantomData, diff --git a/src/ciarlet/nedelec.rs b/src/ciarlet/nedelec.rs index 23bd707..0aa2567 100644 --- a/src/ciarlet/nedelec.rs +++ b/src/ciarlet/nedelec.rs @@ -643,7 +643,7 @@ pub fn create( /// /// A family of Nedelec elements on multiple cell types with appropriate /// continuity across different cell types. -pub struct NedelecFirstKindElementFamily { +pub struct NedelecFirstKindElementFamily { degree: usize, continuity: Continuity, _t: PhantomData, diff --git a/src/ciarlet/raviart_thomas.rs b/src/ciarlet/raviart_thomas.rs index a56b820..4fb707a 100644 --- a/src/ciarlet/raviart_thomas.rs +++ b/src/ciarlet/raviart_thomas.rs @@ -486,7 +486,7 @@ pub fn create( /// /// A family of Raviart-Thomas elements on multiple cell types with appropriate /// continuity across different cell types. -pub struct RaviartThomasElementFamily { +pub struct RaviartThomasElementFamily { degree: usize, continuity: Continuity, _t: PhantomData, From 94b601030d13a8d274e5e4890fad49570e961231 Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Tue, 25 Nov 2025 11:17:56 +0000 Subject: [PATCH 3/9] replace more real with TGeo scalar type --- src/ciarlet.rs | 130 ++++++++++++++++++------------------ src/map.rs | 54 ++++++++------- src/polynomials/legendre.rs | 18 +++-- src/traits.rs | 39 ++++++----- 4 files changed, 130 insertions(+), 111 deletions(-) diff --git a/src/ciarlet.rs b/src/ciarlet.rs index fcb702f..1259238 100644 --- a/src/ciarlet.rs +++ b/src/ciarlet.rs @@ -35,7 +35,6 @@ use crate::reference_cell; use crate::traits::{FiniteElement, Map, MappedFiniteElement}; use crate::types::{Continuity, DofTransformation, ReferenceCellType, Transformation}; use itertools::izip; -use num::{One, Zero}; use rlst::RlstScalar; use rlst::{Array, DynArray, Inverse, MutableArrayImpl, ValueArrayImpl, rlst_dynamic_array}; use std::collections::HashMap; @@ -66,7 +65,7 @@ fn compute_derivative_count(nderivs: usize, cell_type: ReferenceCellType) -> usi } /// A Ciarlet element -pub struct CiarletElement { +pub struct CiarletElement::Real> { family_name: String, cell_type: ReferenceCellType, degree: usize, @@ -78,13 +77,13 @@ pub struct CiarletElement { coefficients: DynArray, entity_dofs: [Vec>; 4], entity_closure_dofs: [Vec>; 4], - interpolation_points: EntityPoints, + interpolation_points: EntityPoints, interpolation_weights: EntityWeights, map: M, dof_transformations: HashMap<(ReferenceCellType, Transformation), DofTransformation>, } -impl Debug for CiarletElement { +impl Debug for CiarletElement { fn fmt(&self, f: &mut Formatter<'_>) -> Result<(), std::fmt::Error> { f.debug_tuple("CiarletElement") .field(&self.family_name) @@ -94,7 +93,7 @@ impl Debug for CiarletElement { } } -impl CiarletElement +impl CiarletElement where DynArray: Inverse>, { @@ -112,7 +111,7 @@ where degree: usize, value_shape: Vec, polynomial_coeffs: DynArray, - interpolation_points: EntityPoints, + interpolation_points: EntityPoints, interpolation_weights: EntityWeights, continuity: Continuity, embedded_superdegree: usize, @@ -144,11 +143,11 @@ where // Format the interpolation points and weights let new_pts = if continuity == Continuity::Discontinuous { - let mut new_pts: EntityPoints = [vec![], vec![], vec![], vec![]]; - let mut all_pts = rlst_dynamic_array!(T::Real, [tdim, npts]); + let mut new_pts: EntityPoints = [vec![], vec![], vec![], vec![]]; + let mut all_pts = rlst_dynamic_array!(TGeo, [tdim, npts]); for (i, pts_i) in interpolation_points.iter().take(tdim).enumerate() { for _pts in pts_i { - new_pts[i].push(rlst_dynamic_array!(T::Real, [tdim, 0])); + new_pts[i].push(rlst_dynamic_array!(TGeo, [tdim, 0])); } } let mut col = 0; @@ -303,32 +302,32 @@ where ReferenceCellType::Interval, 0, Transformation::Reflection, - (|x| vec![x[1], x[0]]) as fn(&[T::Real]) -> Vec, + (|x| vec![x[1], x[0]]) as fn(&[TGeo]) -> Vec, )], ReferenceCellType::Quadrilateral => vec![( ReferenceCellType::Interval, 0, Transformation::Reflection, - (|x| vec![T::Real::one() - x[0], x[1]]) as fn(&[T::Real]) -> Vec, + (|x| vec![TGeo::one() - x[0], x[1]]) as fn(&[TGeo]) -> Vec, )], ReferenceCellType::Tetrahedron => vec![ ( ReferenceCellType::Interval, 0, Transformation::Reflection, - (|x| vec![x[0], x[2], x[1]]) as fn(&[T::Real]) -> Vec, + (|x| vec![x[0], x[2], x[1]]) as fn(&[TGeo]) -> Vec, ), ( ReferenceCellType::Triangle, 0, Transformation::Rotation, - (|x| vec![x[1], x[2], x[0]]) as fn(&[T::Real]) -> Vec, + (|x| vec![x[1], x[2], x[0]]) as fn(&[TGeo]) -> Vec, ), ( ReferenceCellType::Triangle, 0, Transformation::Reflection, - (|x| vec![x[0], x[2], x[1]]) as fn(&[T::Real]) -> Vec, + (|x| vec![x[0], x[2], x[1]]) as fn(&[TGeo]) -> Vec, ), ], ReferenceCellType::Hexahedron => vec![ @@ -336,19 +335,19 @@ where ReferenceCellType::Interval, 0, Transformation::Reflection, - (|x| vec![T::Real::one() - x[0], x[1], x[2]]) as fn(&[T::Real]) -> Vec, + (|x| vec![TGeo::one() - x[0], x[1], x[2]]) as fn(&[TGeo]) -> Vec, ), ( ReferenceCellType::Quadrilateral, 0, Transformation::Rotation, - (|x| vec![x[1], T::Real::one() - x[0], x[2]]) as fn(&[T::Real]) -> Vec, + (|x| vec![x[1], TGeo::one() - x[0], x[2]]) as fn(&[TGeo]) -> Vec, ), ( ReferenceCellType::Quadrilateral, 0, Transformation::Reflection, - (|x| vec![x[1], x[0], x[2]]) as fn(&[T::Real]) -> Vec, + (|x| vec![x[1], x[0], x[2]]) as fn(&[TGeo]) -> Vec, ), ], ReferenceCellType::Prism => vec![ @@ -356,32 +355,32 @@ where ReferenceCellType::Interval, 0, Transformation::Reflection, - (|x| vec![T::Real::one() - x[0], x[1], x[2]]) as fn(&[T::Real]) -> Vec, + (|x| vec![TGeo::one() - x[0], x[1], x[2]]) as fn(&[TGeo]) -> Vec, ), ( ReferenceCellType::Triangle, 0, Transformation::Rotation, - (|x| vec![x[1], T::Real::one() - x[1] - x[0], x[2]]) - as fn(&[T::Real]) -> Vec, + (|x| vec![x[1], TGeo::one() - x[1] - x[0], x[2]]) + as fn(&[TGeo]) -> Vec, ), ( ReferenceCellType::Triangle, 0, Transformation::Reflection, - (|x| vec![x[1], x[0], x[2]]) as fn(&[T::Real]) -> Vec, + (|x| vec![x[1], x[0], x[2]]) as fn(&[TGeo]) -> Vec, ), ( ReferenceCellType::Quadrilateral, 1, Transformation::Rotation, - (|x| vec![x[2], T::Real::one() - x[1], x[0]]) as fn(&[T::Real]) -> Vec, + (|x| vec![x[2], TGeo::one() - x[1], x[0]]) as fn(&[TGeo]) -> Vec, ), ( ReferenceCellType::Quadrilateral, 1, Transformation::Reflection, - (|x| vec![x[2], x[1], x[0]]) as fn(&[T::Real]) -> Vec, + (|x| vec![x[2], x[1], x[0]]) as fn(&[TGeo]) -> Vec, ), ], ReferenceCellType::Pyramid => vec![ @@ -389,32 +388,32 @@ where ReferenceCellType::Interval, 0, Transformation::Reflection, - (|x| vec![T::Real::one() - x[0], x[1], x[2]]) as fn(&[T::Real]) -> Vec, + (|x| vec![TGeo::one() - x[0], x[1], x[2]]) as fn(&[TGeo]) -> Vec, ), ( ReferenceCellType::Triangle, 1, Transformation::Rotation, - (|x| vec![x[2], x[1], T::Real::one() - x[2] - x[0]]) - as fn(&[T::Real]) -> Vec, + (|x| vec![x[2], x[1], TGeo::one() - x[2] - x[0]]) + as fn(&[TGeo]) -> Vec, ), ( ReferenceCellType::Triangle, 1, Transformation::Reflection, - (|x| vec![x[2], x[1], x[0]]) as fn(&[T::Real]) -> Vec, + (|x| vec![x[2], x[1], x[0]]) as fn(&[TGeo]) -> Vec, ), ( ReferenceCellType::Quadrilateral, 0, Transformation::Rotation, - (|x| vec![x[1], T::Real::one() - x[0], x[2]]) as fn(&[T::Real]) -> Vec, + (|x| vec![x[1], TGeo::one() - x[0], x[2]]) as fn(&[TGeo]) -> Vec, ), ( ReferenceCellType::Quadrilateral, 0, Transformation::Reflection, - (|x| vec![x[1], x[0], x[2]]) as fn(&[T::Real]) -> Vec, + (|x| vec![x[1], x[0], x[2]]) as fn(&[TGeo]) -> Vec, ), ], } { @@ -424,26 +423,26 @@ where let finv = match transform { Transformation::Reflection => { - (|x, f| f(x)) as fn(&[T::Real], fn(&[T::Real]) -> Vec) -> Vec + (|x, f| f(x)) as fn(&[TGeo], fn(&[TGeo]) -> Vec) -> Vec } Transformation::Rotation => match entity { ReferenceCellType::Interval => { (|x, f| f(x)) - as fn(&[T::Real], fn(&[T::Real]) -> Vec) -> Vec + as fn(&[TGeo], fn(&[TGeo]) -> Vec) -> Vec } ReferenceCellType::Triangle => { (|x, f| f(&f(x))) - as fn(&[T::Real], fn(&[T::Real]) -> Vec) -> Vec + as fn(&[TGeo], fn(&[TGeo]) -> Vec) -> Vec } ReferenceCellType::Quadrilateral => { (|x, f| f(&f(&f(x)))) - as fn(&[T::Real], fn(&[T::Real]) -> Vec) -> Vec + as fn(&[TGeo], fn(&[TGeo]) -> Vec) -> Vec } _ => panic!("Unsupported entity: {entity:?}"), }, }; - let mut pts = DynArray::::from_shape(ref_pts.shape()); + let mut pts = DynArray::::from_shape(ref_pts.shape()); for p in 0..npts { for (i, c) in finv( &ref_pts.data().unwrap()[ref_pts.shape()[0] * p..ref_pts.shape()[0] * (p + 1)], @@ -459,13 +458,13 @@ where let wts = &new_wts[edim][entity_id]; let edofs = &entity_dofs[edim][entity_id]; let endofs = edofs.len(); - let mut j = rlst_dynamic_array![T::Real, [npts, tdim, tdim]]; + let mut j = rlst_dynamic_array![TGeo, [npts, tdim, tdim]]; for t_in in 0..tdim { for (t_out, (a, b)) in izip!( - f(&vec![T::Real::zero(); tdim]), + f(&vec![TGeo::zero(); tdim]), f(&{ - let mut axis = vec![T::Real::zero(); tdim]; - axis[t_in] = T::Real::one(); + let mut axis = vec![TGeo::zero(); tdim]; + axis[t_in] = TGeo::one(); axis }) ) @@ -479,19 +478,19 @@ where // f is linear. Use this to compute determinants let jdet = vec![ match transform { - Transformation::Reflection => -T::Real::one(), - Transformation::Rotation => T::Real::one(), + Transformation::Reflection => -TGeo::one(), + Transformation::Rotation => TGeo::one(), }; npts ]; - let mut jinv = rlst_dynamic_array![T::Real, [npts, tdim, tdim]]; + let mut jinv = rlst_dynamic_array![TGeo, [npts, tdim, tdim]]; for t_in in 0..tdim { for (t_out, (a, b)) in izip!( - finv(&vec![T::Real::zero(); tdim], f), + finv(&vec![TGeo::zero(); tdim], f), finv( &{ - let mut axis = vec![T::Real::zero(); tdim]; - axis[t_in] = T::Real::one(); + let mut axis = vec![TGeo::zero(); tdim]; + axis[t_in] = TGeo::one(); axis }, f @@ -500,7 +499,7 @@ where .enumerate() { for p in 0..npts { - *jinv.get_mut([p, t_out, t_in]).unwrap() = (b - a).re(); + *jinv.get_mut([p, t_out, t_in]).unwrap() = TGeo::from(b - a).unwrap(); } } } @@ -585,7 +584,7 @@ where ); } } - CiarletElement:: { + CiarletElement:: { family_name, cell_type, degree, @@ -613,7 +612,7 @@ where self.continuity } /// The interpolation points - pub fn interpolation_points(&self) -> &EntityPoints { + pub fn interpolation_points(&self) -> &EntityPoints { &self.interpolation_points } /// The interpolation weights @@ -622,7 +621,7 @@ where } } -impl FiniteElement for CiarletElement { +impl FiniteElement for CiarletElement { type CellType = ReferenceCellType; type T = T; @@ -643,8 +642,9 @@ impl FiniteElement for CiarletElement { } fn tabulate< - Array2Impl: ValueArrayImpl<::Real, 2>, - Array4MutImpl: MutableArrayImpl, + TGeo: RlstScalar, + Array2Impl: ValueArrayImpl, + Array4MutImpl: MutableArrayImpl, >( &self, points: &Array, @@ -708,7 +708,7 @@ impl FiniteElement for CiarletElement { } } -impl MappedFiniteElement for CiarletElement { +impl MappedFiniteElement for CiarletElement { type TransformationType = Transformation; fn lagrange_superdegree(&self) -> usize { @@ -716,16 +716,17 @@ impl MappedFiniteElement for CiarletElement { } fn push_forward< - Array3RealImpl: ValueArrayImpl<::Real, 3>, - Array4Impl: ValueArrayImpl, - Array4MutImpl: MutableArrayImpl, + TGeo: RlstScalar, + Array3GeoImpl: ValueArrayImpl, + Array4Impl: ValueArrayImpl, + Array4MutImpl: MutableArrayImpl, >( &self, reference_values: &Array, nderivs: usize, - jacobians: &Array, - jacobian_determinants: &[::Real], - inverse_jacobians: &Array, + jacobians: &Array, + jacobian_determinants: &[TGeo], + inverse_jacobians: &Array, physical_values: &mut Array, ) { self.map.push_forward( @@ -739,16 +740,17 @@ impl MappedFiniteElement for CiarletElement { } fn pull_back< - Array3RealImpl: ValueArrayImpl<::Real, 3>, - Array4Impl: ValueArrayImpl, - Array4MutImpl: MutableArrayImpl, + TGeo: RlstScalar, + Array3GeoImpl: ValueArrayImpl, + Array4Impl: ValueArrayImpl, + Array4MutImpl: MutableArrayImpl, >( &self, physical_values: &Array, nderivs: usize, - jacobians: &Array, - jacobian_determinants: &[::Real], - inverse_jacobians: &Array, + jacobians: &Array, + jacobian_determinants: &[TGeo], + inverse_jacobians: &Array, reference_values: &mut Array, ) { self.map.pull_back( @@ -847,7 +849,7 @@ impl MappedFiniteElement for CiarletElement { } } - fn apply_dof_transformations(&self, data: &mut [Self::T], cell_orientation: i32) { + fn apply_dof_transformations(&self, data: &mut [T], cell_orientation: i32) { debug_assert_eq!(data.len() % self.dim, 0); let block_size = data.len() / self.dim; let tdim = reference_cell::dim(self.cell_type); diff --git a/src/map.rs b/src/map.rs index 12aa282..4359d74 100644 --- a/src/map.rs +++ b/src/map.rs @@ -11,16 +11,17 @@ pub struct IdentityMap {} impl Map for IdentityMap { fn push_forward< T: RlstScalar, - Array3RealImpl: ValueArrayImpl, + TGeo: RlstScalar, + Array3GeoImpl: ValueArrayImpl, Array4Impl: ValueArrayImpl, Array4MutImpl: MutableArrayImpl, >( &self, reference_values: &Array, nderivs: usize, - _jacobians: &Array, - _jacobian_determinants: &[T::Real], - _inverse_jacobians: &Array, + _jacobians: &Array, + _jacobian_determinants: &[TGeo], + _inverse_jacobians: &Array, physical_values: &mut Array, ) { if nderivs > 0 { @@ -31,16 +32,17 @@ impl Map for IdentityMap { } fn pull_back< T: RlstScalar, - Array3RealImpl: ValueArrayImpl, + TGeo: RlstScalar, + Array3GeoImpl: ValueArrayImpl, Array4Impl: ValueArrayImpl, Array4MutImpl: MutableArrayImpl, >( &self, physical_values: &Array, nderivs: usize, - _jacobians: &Array, - _jacobian_determinants: &[T::Real], - _inverse_jacobians: &Array, + _jacobians: &Array, + _jacobian_determinants: &[TGeo], + _inverse_jacobians: &Array, reference_values: &mut Array, ) { if nderivs > 0 { @@ -70,16 +72,17 @@ pub struct CovariantPiolaMap {} impl Map for CovariantPiolaMap { fn push_forward< T: RlstScalar, - Array3RealImpl: ValueArrayImpl, + TGeo: RlstScalar, + Array3GeoImpl: ValueArrayImpl, Array4Impl: ValueArrayImpl, Array4MutImpl: MutableArrayImpl, >( &self, reference_values: &Array, nderivs: usize, - _jacobians: &Array, - _jacobian_determinants: &[T::Real], - inverse_jacobians: &Array, + _jacobians: &Array, + _jacobian_determinants: &[TGeo], + inverse_jacobians: &Array, physical_values: &mut Array, ) { let tdim = inverse_jacobians.shape()[1]; @@ -109,16 +112,17 @@ impl Map for CovariantPiolaMap { } fn pull_back< T: RlstScalar, - Array3RealImpl: ValueArrayImpl, + TGeo: RlstScalar, + Array3GeoImpl: ValueArrayImpl, Array4Impl: ValueArrayImpl, Array4MutImpl: MutableArrayImpl, >( &self, physical_values: &Array, nderivs: usize, - jacobians: &Array, - _jacobian_determinants: &[T::Real], - _inverse_jacobians: &Array, + jacobians: &Array, + _jacobian_determinants: &[TGeo], + _inverse_jacobians: &Array, reference_values: &mut Array, ) { let gdim = jacobians.shape()[1]; @@ -167,16 +171,17 @@ pub struct ContravariantPiolaMap {} impl Map for ContravariantPiolaMap { fn push_forward< T: RlstScalar, - Array3RealImpl: ValueArrayImpl, + TGeo: RlstScalar, + Array3GeoImpl: ValueArrayImpl, Array4Impl: ValueArrayImpl, Array4MutImpl: MutableArrayImpl, >( &self, reference_values: &Array, nderivs: usize, - jacobians: &Array, - jacobian_determinants: &[T::Real], - _inverse_jacobians: &Array, + jacobians: &Array, + jacobian_determinants: &[TGeo], + _inverse_jacobians: &Array, physical_values: &mut Array, ) { let gdim = jacobians.shape()[1]; @@ -207,16 +212,17 @@ impl Map for ContravariantPiolaMap { } fn pull_back< T: RlstScalar, - Array3RealImpl: ValueArrayImpl, + TGeo: RlstScalar, + Array3GeoImpl: ValueArrayImpl, Array4Impl: ValueArrayImpl, Array4MutImpl: MutableArrayImpl, >( &self, physical_values: &Array, nderivs: usize, - _jacobians: &Array, - jacobian_determinants: &[T::Real], - inverse_jacobians: &Array, + _jacobians: &Array, + jacobian_determinants: &[TGeo], + inverse_jacobians: &Array, reference_values: &mut Array, ) { let tdim = inverse_jacobians.shape()[1]; diff --git a/src/polynomials/legendre.rs b/src/polynomials/legendre.rs index ae4449b..3865d62 100644 --- a/src/polynomials/legendre.rs +++ b/src/polynomials/legendre.rs @@ -37,7 +37,8 @@ fn jrc(a: usize, n: usize) -> (T, T, T) { /// Tabulate orthonormal polynomials on an interval fn tabulate_interval< T: RlstScalar, - Array2Impl: ValueArrayImpl, + TGeo: RlstScalar, + Array2Impl: ValueArrayImpl, Array3MutImpl: MutableArrayImpl, >( points: &Array, @@ -98,7 +99,8 @@ fn tabulate_interval< /// Tabulate orthonormal polynomials on a quadrilateral fn tabulate_quadrilateral< T: RlstScalar, - Array2Impl: ValueArrayImpl, + TGeo: RlstScalar, + Array2Impl: ValueArrayImpl, Array3MutImpl: MutableArrayImpl, >( points: &Array, @@ -251,7 +253,8 @@ fn tabulate_quadrilateral< /// Tabulate orthonormal polynomials on a triangle fn tabulate_triangle< T: RlstScalar, - Array2Impl: ValueArrayImpl, + TGeo: RlstScalar, + Array2Impl: ValueArrayImpl, Array3MutImpl: MutableArrayImpl, >( points: &Array, @@ -452,7 +455,8 @@ fn tabulate_triangle< /// Tabulate orthonormal polynomials on a tetrahedron fn tabulate_tetrahedron< T: RlstScalar, - Array2Impl: ValueArrayImpl, + TGeo: RlstScalar, + Array2Impl: ValueArrayImpl, Array3MutImpl: MutableArrayImpl, >( points: &Array, @@ -815,7 +819,8 @@ fn tabulate_tetrahedron< /// Tabulate orthonormal polynomials on a hexahedron fn tabulate_hexahedron< T: RlstScalar, - Array2Impl: ValueArrayImpl, + TGeo: RlstScalar, + Array2Impl: ValueArrayImpl, Array3MutImpl: MutableArrayImpl, >( points: &Array, @@ -1049,7 +1054,8 @@ pub fn shape>( /// Tabulate orthonormal polynomials pub fn tabulate< T: RlstScalar, - Array2Impl: ValueArrayImpl, + TGeo: RlstScalar, + Array2Impl: ValueArrayImpl, Array3MutImpl: MutableArrayImpl, >( cell_type: ReferenceCellType, diff --git a/src/traits.rs b/src/traits.rs index 7caae3d..ae21797 100644 --- a/src/traits.rs +++ b/src/traits.rs @@ -54,7 +54,8 @@ pub trait FiniteElement { /// /// For the quadrilaterals and hexahedra, the parameter $n$ denotes the Lagrange superdegree. fn tabulate< - Array2Impl: ValueArrayImpl<::Real, 2>, + TGeo: RlstScalar, + Array2Impl: ValueArrayImpl, Array4MutImpl: MutableArrayImpl, >( &self, @@ -122,16 +123,17 @@ pub trait MappedFiniteElement: FiniteElement { /// - `physical_values`: The output array of the push operation. This shape of this array is the same as the `reference_values` /// input, with the [MappedFiniteElement::physical_value_size] used instead of the reference value size. fn push_forward< - Array3RealImpl: ValueArrayImpl<::Real, 3>, + TGeo: RlstScalar, + Array3GeoImpl: ValueArrayImpl, Array4Impl: ValueArrayImpl, Array4MutImpl: MutableArrayImpl, >( &self, reference_values: &Array, nderivs: usize, - jacobians: &Array, - jacobian_determinants: &[::Real], - inverse_jacobians: &Array, + jacobians: &Array, + jacobian_determinants: &[TGeo], + inverse_jacobians: &Array, physical_values: &mut Array, ); @@ -139,16 +141,17 @@ pub trait MappedFiniteElement: FiniteElement { /// /// This is the inverse operation to [MappedFiniteElement::push_forward]. fn pull_back< - Array3RealImpl: ValueArrayImpl<::Real, 3>, + TGeo: RlstScalar, + Array3GeoImpl: ValueArrayImpl, Array4Impl: ValueArrayImpl, Array4MutImpl: MutableArrayImpl, >( &self, physical_values: &Array, nderivs: usize, - jacobians: &Array, - jacobian_determinants: &[::Real], - inverse_jacobians: &Array, + jacobians: &Array, + jacobian_determinants: &[TGeo], + inverse_jacobians: &Array, reference_values: &mut Array, ); @@ -207,32 +210,34 @@ pub trait Map { /// Push values from the reference cell to physical cells. fn push_forward< T: RlstScalar, - Array3RealImpl: ValueArrayImpl, + TGeo: RlstScalar, + Array3GeoImpl: ValueArrayImpl, Array4Impl: ValueArrayImpl, Array4MutImpl: MutableArrayImpl, >( &self, reference_values: &Array, nderivs: usize, - jacobians: &Array, - jacobian_determinants: &[T::Real], - inverse_jacobians: &Array, + jacobians: &Array, + jacobian_determinants: &[TGeo], + inverse_jacobians: &Array, physical_values: &mut Array, ); /// Pull values back to the reference cell. fn pull_back< T: RlstScalar, - Array3RealImpl: ValueArrayImpl, + TGeo: RlstScalar, + Array3GeoImpl: ValueArrayImpl, Array4Impl: ValueArrayImpl, Array4MutImpl: MutableArrayImpl, >( &self, physical_values: &Array, nderivs: usize, - jacobians: &Array, - jacobian_determinants: &[T::Real], - inverse_jacobians: &Array, + jacobians: &Array, + jacobian_determinants: &[TGeo], + inverse_jacobians: &Array, reference_values: &mut Array, ); From 289b4b3cd58aaacb2bae3856820fb43a30d144b9 Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Tue, 25 Nov 2025 11:20:04 +0000 Subject: [PATCH 4/9] remove another real --- src/map.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/map.rs b/src/map.rs index 4359d74..5ae2ddb 100644 --- a/src/map.rs +++ b/src/map.rs @@ -262,7 +262,7 @@ mod test { use approx::*; use rlst::{Array, rlst_dynamic_array}; - fn fill_jacobians, Array3MutImpl: MutableArrayImpl>( + fn fill_jacobians>( j: &mut Array, jdet: &mut [T], jinv: &mut Array, From 97aded2aadae77e37b6d3e4088eb7a2d97b7d979 Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Tue, 25 Nov 2025 11:50:47 +0000 Subject: [PATCH 5/9] remove real from bindings Also cargo fmt on other files --- python/ndelement/ciarlet.py | 110 ++++++++++----- src/bindings.rs | 261 +++++++++++++++++++++++++++--------- src/ciarlet.rs | 23 ++-- 3 files changed, 288 insertions(+), 106 deletions(-) diff --git a/python/ndelement/ciarlet.py b/python/ndelement/ciarlet.py index 5a492de..e1f585d 100644 --- a/python/ndelement/ciarlet.py +++ b/python/ndelement/ciarlet.py @@ -159,13 +159,26 @@ def tabulate(self, points: npt.NDArray[np.floating], nderivs: int) -> npt.NDArra self._rs_element, nderivs, points.shape[0], _ffi.cast("uintptr_t*", shape.ctypes.data) ) data = np.empty(shape[::-1], dtype=self.dtype) - _lib.ciarlet_element_tabulate( - self._rs_element, - _ffi.cast("void*", points.ctypes.data), - points.shape[0], - nderivs, - _ffi.cast("void*", data.ctypes.data), - ) + + if points.dtype == np.float64: + _lib.ciarlet_element_tabulate_f64( + self._rs_element, + _ffi.cast("void*", points.ctypes.data), + points.shape[0], + nderivs, + _ffi.cast("void*", data.ctypes.data), + ) + elif points.dtype == np.float32: + _lib.ciarlet_element_tabulate_f32( + self._rs_element, + points.ctypes.data, + points.shape[0], + nderivs, + _ffi.cast("void*", data.ctypes.data), + ) + else: + raise TypeError(f"Unsupported dtype: {points.dtype}") + return data def physical_value_size(self, gdim: int) -> int: @@ -205,18 +218,36 @@ def push_forward( shape = (pvs,) + reference_values.shape[1:] data = np.empty(shape, dtype=self.dtype) - _lib.ciarlet_element_push_forward( - self._rs_element, - npts, - nfuncs, - gdim, - _ffi.cast("void*", reference_values.ctypes.data), - nderivs, - _ffi.cast("void*", jacobians.ctypes.data), - _ffi.cast("void*", jacobian_determinants.ctypes.data), - _ffi.cast("void*", inverse_jacobians.ctypes.data), - _ffi.cast("void*", data.ctypes.data), - ) + + if reference_values.dtype == np.float64: + _lib.ciarlet_element_push_forward_f64( + self._rs_element, + npts, + nfuncs, + gdim, + reference_values.ctypes.data, + nderivs, + jacobians.ctypes.data, + jacobian_determinants.ctypes.data, + inverse_jacobians.ctypes.data, + _ffi.cast("void*", data.ctypes.data), + ) + elif reference_values.dtype == np.float32: + _lib.ciarlet_element_push_forward_f32( + self._rs_element, + npts, + nfuncs, + gdim, + reference_values.ctypes.data, + nderivs, + jacobians.ctypes.data, + jacobian_determinants.ctypes.data, + inverse_jacobians.ctypes.data, + _ffi.cast("void*", data.ctypes.data), + ) + else: + raise TypeError(f"Unsupported dtype: {reference_values.dtype}") + return data def pull_back( @@ -244,18 +275,35 @@ def pull_back( shape = (vs,) + physical_values.shape[1:] data = np.empty(shape, dtype=self.dtype) - _lib.ciarlet_element_pull_back( - self._rs_element, - npts, - nfuncs, - gdim, - _ffi.cast("void*", physical_values.ctypes.data), - nderivs, - _ffi.cast("void*", jacobians.ctypes.data), - _ffi.cast("void*", jacobian_determinants.ctypes.data), - _ffi.cast("void*", inverse_jacobians.ctypes.data), - _ffi.cast("void*", data.ctypes.data), - ) + + if physical_values.dtype == np.float64: + _lib.ciarlet_element_push_forward_f64( + self._rs_element, + npts, + nfuncs, + gdim, + physical_values.ctypes.data, + nderivs, + jacobians.ctypes.data, + jacobian_determinants.ctypes.data, + inverse_jacobians.ctypes.data, + _ffi.cast("void*", data.ctypes.data), + ) + elif physical_values.dtype == np.float32: + _lib.ciarlet_element_push_forward_f32( + self._rs_element, + npts, + nfuncs, + gdim, + physical_values.ctypes.data, + nderivs, + jacobians.ctypes.data, + jacobian_determinants.ctypes.data, + inverse_jacobians.ctypes.data, + _ffi.cast("void*", data.ctypes.data), + ) + else: + raise TypeError(f"Unsupported dtype: {physical_values.dtype}") return data def apply_dof_permutations(self, data: npt.NDArray, orientation: np.int32): diff --git a/src/bindings.rs b/src/bindings.rs index fb6110c..cbe44da 100644 --- a/src/bindings.rs +++ b/src/bindings.rs @@ -15,7 +15,7 @@ pub mod reference_cell { pub unsafe extern "C" fn is_simplex(cell: ReferenceCellType) -> bool { reference_cell::is_simplex(cell) } - unsafe fn vertices>(cell: ReferenceCellType, vs: *mut T) { + unsafe fn vertices(cell: ReferenceCellType, vs: *mut T) { let mut i = 0; for v in reference_cell::vertices::(cell) { for c in v { @@ -38,7 +38,7 @@ pub mod reference_cell { vertices(cell, vs); } } - unsafe fn midpoint>(cell: ReferenceCellType, pt: *mut T) { + unsafe fn midpoint(cell: ReferenceCellType, pt: *mut T) { for (i, c) in reference_cell::midpoint(cell).iter().enumerate() { unsafe { *pt.add(i) = *c; @@ -199,16 +199,16 @@ pub mod polynomials { *shape.add(2) = npts; } } - unsafe fn tabulate_legendre_polynomials( + unsafe fn tabulate_legendre_polynomials( cell: ReferenceCellType, - points: *const T::Real, + points: *const TGeo, npts: usize, degree: usize, derivatives: usize, data: *mut T, ) { let tdim = reference_cell::dim(cell); - let points = SliceArray::::from_shape( + let points = SliceArray::::from_shape( unsafe { from_raw_parts(points, npts * tdim) }, [tdim, npts], ); @@ -432,22 +432,19 @@ pub mod ciarlet { } } - #[concretise_types( - gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - 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_tabulate>( + pub fn ciarlet_element_tabulate< + E: FiniteElement, + TGeo: RlstScalar, + >( element: &E, - points: *const c_void, + points: *const TGeo, npoints: usize, nderivs: usize, data: *mut c_void, ) { let tdim = reference_cell::dim(element.cell_type()); - let points = points as *mut ::Real; let data = data as *mut E::T; - let points = SliceArray::<::Real, 2>::from_shape( + let points = SliceArray::::from_shape( unsafe { from_raw_parts(points, npoints * tdim) }, [tdim, npoints], ); @@ -459,6 +456,36 @@ pub mod ciarlet { element.tabulate(&points, nderivs, &mut data); } + #[concretise_types( + gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + 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_tabulate_f32>( + element: &E, + points: *const f32, + npoints: usize, + nderivs: usize, + data: *mut c_void, + ) { + ciarlet_element_tabulate(element, points, npoints, nderivs, data); + } + + #[concretise_types( + gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + 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_tabulate_f64>( + element: &E, + points: *const f64, + npoints: usize, + nderivs: usize, + data: *mut c_void, + ) { + ciarlet_element_tabulate(element, points, npoints, nderivs, data); + } + #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), @@ -531,22 +558,19 @@ pub mod ciarlet { } } - #[allow(clippy::too_many_arguments)] - #[concretise_types( - gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - 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< + E: MappedFiniteElement, + TGeo: RlstScalar, + >( element: &E, npoints: usize, nfunctions: usize, gdim: usize, - reference_values: *const c_void, + reference_values: *const TGeo, nderivs: usize, - j: *const c_void, - jdet: *const c_void, - jinv: *const c_void, + j: *const TGeo, + jdet: *const TGeo, + jinv: *const TGeo, physical_values: *mut c_void, ) { let tdim = reference_cell::dim(element.cell_type()); @@ -562,23 +586,13 @@ pub mod ciarlet { }, [deriv_size, npoints, nfunctions, vs], ); - let j = SliceArray::<::Real, 3>::from_shape( - unsafe { - from_raw_parts( - j as *const ::Real, - npoints * gdim * tdim, - ) - }, + let j = SliceArray::::from_shape( + unsafe { from_raw_parts(j, npoints * gdim * tdim) }, [npoints, gdim, tdim], ); - let jdet = unsafe { from_raw_parts(jdet as *const ::Real, npoints) }; - let jinv = SliceArray::<::Real, 3>::from_shape( - unsafe { - from_raw_parts( - jinv as *const ::Real, - npoints * tdim * gdim, - ) - }, + let jdet = unsafe { from_raw_parts(jdet, npoints) }; + let jinv = SliceArray::::from_shape( + unsafe { from_raw_parts(jinv, npoints * tdim * gdim) }, [npoints, tdim, gdim], ); let mut physical_values = SliceArrayMut::::from_shape( @@ -606,16 +620,81 @@ 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_push_forward_f32< + E: MappedFiniteElement, + >( + element: &E, + npoints: usize, + nfunctions: usize, + gdim: usize, + reference_values: *const f32, + nderivs: usize, + j: *const f32, + jdet: *const f32, + jinv: *const f32, + physical_values: *mut c_void, + ) { + ciarlet_element_push_forward( + element, + npoints, + nfunctions, + gdim, + reference_values, + nderivs, + j, + jdet, + jinv, + physical_values, + ); + } + + #[allow(clippy::too_many_arguments)] + #[concretise_types( + gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + 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_f64< + E: MappedFiniteElement, + >( + element: &E, + npoints: usize, + nfunctions: usize, + gdim: usize, + reference_values: *const f64, + nderivs: usize, + j: *const f64, + jdet: *const f64, + jinv: *const f64, + physical_values: *mut c_void, + ) { + ciarlet_element_push_forward( + element, + npoints, + nfunctions, + gdim, + reference_values, + nderivs, + j, + jdet, + jinv, + physical_values, + ); + } + + pub fn ciarlet_element_pull_back< + E: MappedFiniteElement, + TGeo: RlstScalar, + >( element: &E, npoints: usize, nfunctions: usize, gdim: usize, - physical_values: *const c_void, + physical_values: *const TGeo, nderivs: usize, - j: *const c_void, - jdet: *const c_void, - jinv: *const c_void, + j: *const TGeo, + jdet: *const TGeo, + jinv: *const TGeo, reference_values: *mut c_void, ) { let tdim = reference_cell::dim(element.cell_type()); @@ -631,23 +710,13 @@ pub mod ciarlet { }, [deriv_size, npoints, nfunctions, pvs], ); - let j = SliceArray::<::Real, 3>::from_shape( - unsafe { - from_raw_parts( - j as *const ::Real, - npoints * gdim * tdim, - ) - }, + let j = SliceArray::::from_shape( + unsafe { from_raw_parts(j, npoints * gdim * tdim) }, [npoints, gdim, tdim], ); - let jdet = unsafe { from_raw_parts(jdet as *const ::Real, npoints) }; - let jinv = SliceArray::<::Real, 3>::from_shape( - unsafe { - from_raw_parts( - jinv as *const ::Real, - npoints * tdim * gdim, - ) - }, + let jdet = unsafe { from_raw_parts(jdet, npoints) }; + let jinv = SliceArray::::from_shape( + unsafe { from_raw_parts(jinv, npoints * tdim * gdim) }, [npoints, tdim, gdim], ); let mut reference_values = SliceArrayMut::::from_shape( @@ -669,6 +738,70 @@ pub mod ciarlet { ); } + #[allow(clippy::too_many_arguments)] + #[concretise_types( + gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + 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_f32>( + element: &E, + npoints: usize, + nfunctions: usize, + gdim: usize, + physical_values: *const f32, + nderivs: usize, + j: *const f32, + jdet: *const f32, + jinv: *const f32, + reference_values: *mut c_void, + ) { + ciarlet_element_pull_back( + element, + npoints, + nfunctions, + gdim, + physical_values, + nderivs, + j, + jdet, + jinv, + reference_values, + ); + } + + #[allow(clippy::too_many_arguments)] + #[concretise_types( + gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + 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_f64>( + element: &E, + npoints: usize, + nfunctions: usize, + gdim: usize, + physical_values: *const f64, + nderivs: usize, + j: *const f64, + jdet: *const f64, + jinv: *const f64, + reference_values: *mut c_void, + ) { + ciarlet_element_pull_back( + element, + npoints, + nfunctions, + gdim, + physical_values, + nderivs, + j, + jdet, + jinv, + reference_values, + ); + } + #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), @@ -829,19 +962,21 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]) )] pub fn ciarlet_element_interpolation_points< T: RlstScalar + DTypeIdentifier + Getrf + Getri, + TGeo: RlstScalar + DTypeIdentifier + Getrf + Getri, M: Map, >( - element: &CiarletElement, + element: &CiarletElement, entity_dim: usize, entity_index: usize, points: *mut c_void, ) { - let points = points as *mut T::Real; + let points = points as *mut TGeo; for (i, j) in element.interpolation_points()[entity_dim][entity_index] .data() .unwrap() diff --git a/src/ciarlet.rs b/src/ciarlet.rs index 1259238..94b05cf 100644 --- a/src/ciarlet.rs +++ b/src/ciarlet.rs @@ -361,8 +361,7 @@ where ReferenceCellType::Triangle, 0, Transformation::Rotation, - (|x| vec![x[1], TGeo::one() - x[1] - x[0], x[2]]) - as fn(&[TGeo]) -> Vec, + (|x| vec![x[1], TGeo::one() - x[1] - x[0], x[2]]) as fn(&[TGeo]) -> Vec, ), ( ReferenceCellType::Triangle, @@ -394,8 +393,7 @@ where ReferenceCellType::Triangle, 1, Transformation::Rotation, - (|x| vec![x[2], x[1], TGeo::one() - x[2] - x[0]]) - as fn(&[TGeo]) -> Vec, + (|x| vec![x[2], x[1], TGeo::one() - x[2] - x[0]]) as fn(&[TGeo]) -> Vec, ), ( ReferenceCellType::Triangle, @@ -427,16 +425,13 @@ where } Transformation::Rotation => match entity { ReferenceCellType::Interval => { - (|x, f| f(x)) - as fn(&[TGeo], fn(&[TGeo]) -> Vec) -> Vec + (|x, f| f(x)) as fn(&[TGeo], fn(&[TGeo]) -> Vec) -> Vec } ReferenceCellType::Triangle => { - (|x, f| f(&f(x))) - as fn(&[TGeo], fn(&[TGeo]) -> Vec) -> Vec + (|x, f| f(&f(x))) as fn(&[TGeo], fn(&[TGeo]) -> Vec) -> Vec } ReferenceCellType::Quadrilateral => { - (|x, f| f(&f(&f(x)))) - as fn(&[TGeo], fn(&[TGeo]) -> Vec) -> Vec + (|x, f| f(&f(&f(x)))) as fn(&[TGeo], fn(&[TGeo]) -> Vec) -> Vec } _ => panic!("Unsupported entity: {entity:?}"), }, @@ -621,7 +616,9 @@ where } } -impl FiniteElement for CiarletElement { +impl FiniteElement + for CiarletElement +{ type CellType = ReferenceCellType; type T = T; @@ -708,7 +705,9 @@ impl FiniteElement for CiarletE } } -impl MappedFiniteElement for CiarletElement { +impl MappedFiniteElement + for CiarletElement +{ type TransformationType = Transformation; fn lagrange_superdegree(&self) -> usize { From 7bedf89cbd17476d765458d8b35627ffc114a0f8 Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Tue, 25 Nov 2025 13:33:23 +0000 Subject: [PATCH 6/9] remove more real --- examples/element_family.rs | 2 +- examples/lagrange_element.rs | 3 +- examples/test_high_degree.rs | 6 +- src/ciarlet.rs | 204 +++++++++++++++++++++------------- src/ciarlet/lagrange.rs | 117 +++++++++---------- src/ciarlet/nedelec.rs | 93 +++++++++------- src/ciarlet/raviart_thomas.rs | 87 ++++++++------- 7 files changed, 292 insertions(+), 220 deletions(-) diff --git a/examples/element_family.rs b/examples/element_family.rs index 26874d4..b9b8666 100644 --- a/examples/element_family.rs +++ b/examples/element_family.rs @@ -5,7 +5,7 @@ use ndelement::types::{Continuity, ReferenceCellType}; fn main() { // 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::::new(2, Continuity::Standard); + let family = LagrangeElementFamily::::new(2, Continuity::Standard); // Get the element in the family on a triangle let element = family.element(ReferenceCellType::Triangle); diff --git a/examples/lagrange_element.rs b/examples/lagrange_element.rs index b3adce6..59f3ed4 100644 --- a/examples/lagrange_element.rs +++ b/examples/lagrange_element.rs @@ -7,7 +7,8 @@ use rlst::{DynArray, rlst_dynamic_array}; fn main() { // Create a P2 element on a triangle - let element = lagrange::create::(ReferenceCellType::Triangle, 2, Continuity::Standard); + let element = + lagrange::create::(ReferenceCellType::Triangle, 2, Continuity::Standard); println!("This element has {} basis functions.", element.dim()); diff --git a/examples/test_high_degree.rs b/examples/test_high_degree.rs index 776a4bf..5ebb0cf 100644 --- a/examples/test_high_degree.rs +++ b/examples/test_high_degree.rs @@ -10,7 +10,7 @@ fn main() { paste! { for d in 1..[<$max_degree>] { println!("Constructing Lagrange(degree={d}, cell={:?})", ReferenceCellType::[<$cell>]); - let _e = lagrange::create::(ReferenceCellType::[<$cell>], d, Continuity::Standard); + let _e = lagrange::create::(ReferenceCellType::[<$cell>], d, Continuity::Standard); } } }; @@ -21,7 +21,7 @@ fn main() { paste! { for d in 1..[<$max_degree>] { println!("Constructing RaviartThomas(degree={d}, cell={:?})", ReferenceCellType::[<$cell>]); - let _e = raviart_thomas::create::(ReferenceCellType::[<$cell>], d, Continuity::Standard); + let _e = raviart_thomas::create::(ReferenceCellType::[<$cell>], d, Continuity::Standard); } } }; @@ -32,7 +32,7 @@ fn main() { paste! { for d in 1..[<$max_degree>] { println!("Constructing Nedelec(degree={d}, cell={:?})", ReferenceCellType::[<$cell>]); - let _e = nedelec::create::(ReferenceCellType::[<$cell>], d, Continuity::Standard); + let _e = nedelec::create::(ReferenceCellType::[<$cell>], d, Continuity::Standard); } } }; diff --git a/src/ciarlet.rs b/src/ciarlet.rs index 94b05cf..38e2701 100644 --- a/src/ciarlet.rs +++ b/src/ciarlet.rs @@ -951,13 +951,14 @@ mod test { #[test] fn test_lagrange_1() { - let e = lagrange::create::(ReferenceCellType::Triangle, 1, Continuity::Standard); + let e = lagrange::create::(ReferenceCellType::Triangle, 1, Continuity::Standard); assert_eq!(e.value_size(), 1); } #[test] fn test_lagrange_0_interval() { - let e = lagrange::create::(ReferenceCellType::Interval, 0, Continuity::Discontinuous); + let e = + lagrange::create::(ReferenceCellType::Interval, 0, Continuity::Discontinuous); assert_eq!(e.value_size(), 1); let mut data = DynArray::::from_shape(e.tabulate_array_shape(0, 4)); let mut points = rlst_dynamic_array!(f64, [1, 4]); @@ -975,7 +976,7 @@ mod test { #[test] fn test_lagrange_1_interval() { - let e = lagrange::create::(ReferenceCellType::Interval, 1, Continuity::Standard); + let e = lagrange::create::(ReferenceCellType::Interval, 1, Continuity::Standard); assert_eq!(e.value_size(), 1); let mut data = DynArray::::from_shape(e.tabulate_array_shape(0, 4)); let mut points = rlst_dynamic_array!(f64, [1, 4]); @@ -995,7 +996,8 @@ mod test { #[test] fn test_lagrange_0_triangle() { - let e = lagrange::create::(ReferenceCellType::Triangle, 0, Continuity::Discontinuous); + let e = + lagrange::create::(ReferenceCellType::Triangle, 0, Continuity::Discontinuous); assert_eq!(e.value_size(), 1); let mut data = DynArray::::from_shape(e.tabulate_array_shape(0, 6)); @@ -1023,7 +1025,7 @@ mod test { #[test] fn test_lagrange_1_triangle() { - let e = lagrange::create::(ReferenceCellType::Triangle, 1, Continuity::Standard); + let e = lagrange::create::(ReferenceCellType::Triangle, 1, Continuity::Standard); assert_eq!(e.value_size(), 1); let mut data = DynArray::::from_shape(e.tabulate_array_shape(0, 6)); let mut points = rlst_dynamic_array!(f64, [2, 6]); @@ -1053,7 +1055,7 @@ mod test { #[test] fn test_lagrange_0_quadrilateral() { - let e = lagrange::create::( + let e = lagrange::create::( ReferenceCellType::Quadrilateral, 0, Continuity::Discontinuous, @@ -1083,7 +1085,8 @@ mod test { #[test] fn test_lagrange_1_quadrilateral() { - let e = lagrange::create::(ReferenceCellType::Quadrilateral, 1, Continuity::Standard); + let e = + lagrange::create::(ReferenceCellType::Quadrilateral, 1, Continuity::Standard); assert_eq!(e.value_size(), 1); let mut data = DynArray::::from_shape(e.tabulate_array_shape(0, 6)); let mut points = rlst_dynamic_array!(f64, [2, 6]); @@ -1115,7 +1118,8 @@ mod test { #[test] fn test_lagrange_2_quadrilateral() { - let e = lagrange::create::(ReferenceCellType::Quadrilateral, 2, Continuity::Standard); + let e = + lagrange::create::(ReferenceCellType::Quadrilateral, 2, Continuity::Standard); assert_eq!(e.value_size(), 1); let mut data = DynArray::::from_shape(e.tabulate_array_shape(0, 6)); let mut points = rlst_dynamic_array!(f64, [2, 6]); @@ -1188,8 +1192,11 @@ mod test { #[test] fn test_lagrange_0_tetrahedron() { - let e = - lagrange::create::(ReferenceCellType::Tetrahedron, 0, Continuity::Discontinuous); + let e = lagrange::create::( + ReferenceCellType::Tetrahedron, + 0, + Continuity::Discontinuous, + ); assert_eq!(e.value_size(), 1); let mut data = DynArray::::from_shape(e.tabulate_array_shape(0, 6)); let mut points = rlst_dynamic_array!(f64, [3, 6]); @@ -1221,7 +1228,8 @@ mod test { #[test] fn test_lagrange_1_tetrahedron() { - let e = lagrange::create::(ReferenceCellType::Tetrahedron, 1, Continuity::Standard); + let e = + lagrange::create::(ReferenceCellType::Tetrahedron, 1, Continuity::Standard); assert_eq!(e.value_size(), 1); let mut data = DynArray::::from_shape(e.tabulate_array_shape(0, 6)); let mut points = rlst_dynamic_array!(f64, [3, 6]); @@ -1264,8 +1272,11 @@ mod test { #[test] fn test_lagrange_0_hexahedron() { - let e = - lagrange::create::(ReferenceCellType::Hexahedron, 0, Continuity::Discontinuous); + let e = lagrange::create::( + ReferenceCellType::Hexahedron, + 0, + Continuity::Discontinuous, + ); assert_eq!(e.value_size(), 1); let mut data = DynArray::::from_shape(e.tabulate_array_shape(0, 6)); let mut points = rlst_dynamic_array!(f64, [3, 6]); @@ -1297,7 +1308,8 @@ mod test { #[test] fn test_lagrange_1_hexahedron() { - let e = lagrange::create::(ReferenceCellType::Hexahedron, 1, Continuity::Standard); + let e = + lagrange::create::(ReferenceCellType::Hexahedron, 1, Continuity::Standard); assert_eq!(e.value_size(), 1); let mut data = DynArray::::from_shape(e.tabulate_array_shape(0, 6)); let mut points = rlst_dynamic_array!(f64, [3, 6]); @@ -1372,53 +1384,53 @@ mod test { #[test] fn test_lagrange_higher_degree_triangle() { - lagrange::create::(ReferenceCellType::Triangle, 2, Continuity::Standard); - lagrange::create::(ReferenceCellType::Triangle, 3, Continuity::Standard); - lagrange::create::(ReferenceCellType::Triangle, 4, Continuity::Standard); - lagrange::create::(ReferenceCellType::Triangle, 5, Continuity::Standard); + lagrange::create::(ReferenceCellType::Triangle, 2, Continuity::Standard); + lagrange::create::(ReferenceCellType::Triangle, 3, Continuity::Standard); + lagrange::create::(ReferenceCellType::Triangle, 4, Continuity::Standard); + lagrange::create::(ReferenceCellType::Triangle, 5, Continuity::Standard); - lagrange::create::(ReferenceCellType::Triangle, 2, Continuity::Discontinuous); - lagrange::create::(ReferenceCellType::Triangle, 3, Continuity::Discontinuous); - lagrange::create::(ReferenceCellType::Triangle, 4, Continuity::Discontinuous); - lagrange::create::(ReferenceCellType::Triangle, 5, Continuity::Discontinuous); + lagrange::create::(ReferenceCellType::Triangle, 2, Continuity::Discontinuous); + lagrange::create::(ReferenceCellType::Triangle, 3, Continuity::Discontinuous); + lagrange::create::(ReferenceCellType::Triangle, 4, Continuity::Discontinuous); + lagrange::create::(ReferenceCellType::Triangle, 5, Continuity::Discontinuous); } #[test] fn test_lagrange_higher_degree_interval() { - lagrange::create::(ReferenceCellType::Interval, 2, Continuity::Standard); - lagrange::create::(ReferenceCellType::Interval, 3, Continuity::Standard); - lagrange::create::(ReferenceCellType::Interval, 4, Continuity::Standard); - lagrange::create::(ReferenceCellType::Interval, 5, Continuity::Standard); + lagrange::create::(ReferenceCellType::Interval, 2, Continuity::Standard); + lagrange::create::(ReferenceCellType::Interval, 3, Continuity::Standard); + lagrange::create::(ReferenceCellType::Interval, 4, Continuity::Standard); + lagrange::create::(ReferenceCellType::Interval, 5, Continuity::Standard); - lagrange::create::(ReferenceCellType::Interval, 2, Continuity::Discontinuous); - lagrange::create::(ReferenceCellType::Interval, 3, Continuity::Discontinuous); - lagrange::create::(ReferenceCellType::Interval, 4, Continuity::Discontinuous); - lagrange::create::(ReferenceCellType::Interval, 5, Continuity::Discontinuous); + lagrange::create::(ReferenceCellType::Interval, 2, Continuity::Discontinuous); + lagrange::create::(ReferenceCellType::Interval, 3, Continuity::Discontinuous); + lagrange::create::(ReferenceCellType::Interval, 4, Continuity::Discontinuous); + lagrange::create::(ReferenceCellType::Interval, 5, Continuity::Discontinuous); } #[test] fn test_lagrange_higher_degree_quadrilateral() { - lagrange::create::(ReferenceCellType::Quadrilateral, 2, Continuity::Standard); - lagrange::create::(ReferenceCellType::Quadrilateral, 3, Continuity::Standard); - lagrange::create::(ReferenceCellType::Quadrilateral, 4, Continuity::Standard); - lagrange::create::(ReferenceCellType::Quadrilateral, 5, Continuity::Standard); + lagrange::create::(ReferenceCellType::Quadrilateral, 2, Continuity::Standard); + lagrange::create::(ReferenceCellType::Quadrilateral, 3, Continuity::Standard); + lagrange::create::(ReferenceCellType::Quadrilateral, 4, Continuity::Standard); + lagrange::create::(ReferenceCellType::Quadrilateral, 5, Continuity::Standard); - lagrange::create::( + lagrange::create::( ReferenceCellType::Quadrilateral, 2, Continuity::Discontinuous, ); - lagrange::create::( + lagrange::create::( ReferenceCellType::Quadrilateral, 3, Continuity::Discontinuous, ); - lagrange::create::( + lagrange::create::( ReferenceCellType::Quadrilateral, 4, Continuity::Discontinuous, ); - lagrange::create::( + lagrange::create::( ReferenceCellType::Quadrilateral, 5, Continuity::Discontinuous, @@ -1427,7 +1439,11 @@ mod test { #[test] fn test_raviart_thomas_1_triangle() { - let e = raviart_thomas::create(ReferenceCellType::Triangle, 1, Continuity::Standard); + let e = raviart_thomas::create::( + ReferenceCellType::Triangle, + 1, + Continuity::Standard, + ); assert_eq!(e.value_size(), 2); let mut data = DynArray::::from_shape(e.tabulate_array_shape(0, 6)); let mut points = rlst_dynamic_array!(f64, [2, 6]); @@ -1459,21 +1475,33 @@ mod test { #[test] fn test_raviart_thomas_2_triangle() { - let e = raviart_thomas::create::(ReferenceCellType::Triangle, 2, Continuity::Standard); + let e = raviart_thomas::create::( + ReferenceCellType::Triangle, + 2, + Continuity::Standard, + ); assert_eq!(e.value_size(), 2); check_dofs(e); } #[test] fn test_raviart_thomas_3_triangle() { - let e = raviart_thomas::create::(ReferenceCellType::Triangle, 3, Continuity::Standard); + let e = raviart_thomas::create::( + ReferenceCellType::Triangle, + 3, + Continuity::Standard, + ); assert_eq!(e.value_size(), 2); check_dofs(e); } #[test] fn test_raviart_thomas_1_quadrilateral() { - let e = raviart_thomas::create(ReferenceCellType::Quadrilateral, 1, Continuity::Standard); + let e = raviart_thomas::create::( + ReferenceCellType::Quadrilateral, + 1, + Continuity::Standard, + ); assert_eq!(e.value_size(), 2); let mut data = DynArray::::from_shape(e.tabulate_array_shape(0, 6)); let mut points = rlst_dynamic_array!(f64, [2, 6]); @@ -1508,7 +1536,7 @@ mod test { #[test] fn test_raviart_thomas_2_quadrilateral() { - let e = raviart_thomas::create::( + let e = raviart_thomas::create::( ReferenceCellType::Quadrilateral, 2, Continuity::Standard, @@ -1519,7 +1547,7 @@ mod test { #[test] fn test_raviart_thomas_3_quadrilateral() { - let e = raviart_thomas::create::( + let e = raviart_thomas::create::( ReferenceCellType::Quadrilateral, 3, Continuity::Standard, @@ -1530,31 +1558,44 @@ mod test { #[test] fn test_raviart_thomas_1_tetrahedron() { - let e = - raviart_thomas::create::(ReferenceCellType::Tetrahedron, 1, Continuity::Standard); + let e = raviart_thomas::create::( + ReferenceCellType::Tetrahedron, + 1, + Continuity::Standard, + ); assert_eq!(e.value_size(), 3); check_dofs(e); } #[test] fn test_raviart_thomas_2_tetrahedron() { - let e = - raviart_thomas::create::(ReferenceCellType::Tetrahedron, 2, Continuity::Standard); + let e = raviart_thomas::create::( + ReferenceCellType::Tetrahedron, + 2, + Continuity::Standard, + ); assert_eq!(e.value_size(), 3); check_dofs(e); } #[test] fn test_raviart_thomas_3_tetrahedron() { - let e = - raviart_thomas::create::(ReferenceCellType::Tetrahedron, 3, Continuity::Standard); + let e = raviart_thomas::create::( + ReferenceCellType::Tetrahedron, + 3, + Continuity::Standard, + ); assert_eq!(e.value_size(), 3); check_dofs(e); } #[test] fn test_raviart_thomas_1_hexahedron() { - let e = raviart_thomas::create(ReferenceCellType::Hexahedron, 1, Continuity::Standard); + let e = raviart_thomas::create::( + ReferenceCellType::Hexahedron, + 1, + Continuity::Standard, + ); assert_eq!(e.value_size(), 3); let mut data = DynArray::::from_shape(e.tabulate_array_shape(0, 6)); let mut points = rlst_dynamic_array!(f64, [3, 6]); @@ -1603,15 +1644,18 @@ mod test { #[test] fn test_raviart_thomas_2_hexahedron() { - let e = - raviart_thomas::create::(ReferenceCellType::Hexahedron, 2, Continuity::Standard); + let e = raviart_thomas::create::( + ReferenceCellType::Hexahedron, + 2, + Continuity::Standard, + ); assert_eq!(e.value_size(), 3); check_dofs(e); } #[test] fn test_nedelec_1_triangle() { - let e = nedelec::create(ReferenceCellType::Triangle, 1, Continuity::Standard); + let e = nedelec::create::(ReferenceCellType::Triangle, 1, Continuity::Standard); assert_eq!(e.value_size(), 2); let mut data = DynArray::::from_shape(e.tabulate_array_shape(0, 6)); let mut points = rlst_dynamic_array!(f64, [2, 6]); @@ -1643,21 +1687,22 @@ mod test { #[test] fn test_nedelec_2_triangle() { - let e = nedelec::create::(ReferenceCellType::Triangle, 2, Continuity::Standard); + let e = nedelec::create::(ReferenceCellType::Triangle, 2, Continuity::Standard); assert_eq!(e.value_size(), 2); check_dofs(e); } #[test] fn test_nedelec_3_triangle() { - let e = nedelec::create::(ReferenceCellType::Triangle, 3, Continuity::Standard); + let e = nedelec::create::(ReferenceCellType::Triangle, 3, Continuity::Standard); assert_eq!(e.value_size(), 2); check_dofs(e); } #[test] fn test_nedelec_1_quadrilateral() { - let e = nedelec::create(ReferenceCellType::Quadrilateral, 1, Continuity::Standard); + let e = + nedelec::create::(ReferenceCellType::Quadrilateral, 1, Continuity::Standard); assert_eq!(e.value_size(), 2); let mut data = DynArray::::from_shape(e.tabulate_array_shape(0, 6)); let mut points = rlst_dynamic_array!(f64, [2, 6]); @@ -1692,42 +1737,47 @@ mod test { #[test] fn test_nedelec_2_quadrilateral() { - let e = nedelec::create::(ReferenceCellType::Quadrilateral, 2, Continuity::Standard); + let e = + nedelec::create::(ReferenceCellType::Quadrilateral, 2, Continuity::Standard); assert_eq!(e.value_size(), 2); check_dofs(e); } #[test] fn test_nedelec_3_quadrilateral() { - let e = nedelec::create::(ReferenceCellType::Quadrilateral, 3, Continuity::Standard); + let e = + nedelec::create::(ReferenceCellType::Quadrilateral, 3, Continuity::Standard); assert_eq!(e.value_size(), 2); check_dofs(e); } #[test] fn test_nedelec_1_tetrahedron() { - let e = nedelec::create::(ReferenceCellType::Tetrahedron, 1, Continuity::Standard); + let e = + nedelec::create::(ReferenceCellType::Tetrahedron, 1, Continuity::Standard); assert_eq!(e.value_size(), 3); check_dofs(e); } #[test] fn test_nedelec_2_tetrahedron() { - let e = nedelec::create::(ReferenceCellType::Tetrahedron, 2, Continuity::Standard); + let e = + nedelec::create::(ReferenceCellType::Tetrahedron, 2, Continuity::Standard); assert_eq!(e.value_size(), 3); check_dofs(e); } #[test] fn test_nedelec_3_tetrahedron() { - let e = nedelec::create::(ReferenceCellType::Tetrahedron, 3, Continuity::Standard); + let e = + nedelec::create::(ReferenceCellType::Tetrahedron, 3, Continuity::Standard); assert_eq!(e.value_size(), 3); check_dofs(e); } #[test] fn test_nedelec_1_hexahedron() { - let e = nedelec::create(ReferenceCellType::Hexahedron, 1, Continuity::Standard); + let e = nedelec::create::(ReferenceCellType::Hexahedron, 1, Continuity::Standard); assert_eq!(e.value_size(), 3); let mut data = DynArray::::from_shape(e.tabulate_array_shape(0, 6)); let mut points = rlst_dynamic_array!(f64, [3, 6]); @@ -1782,7 +1832,7 @@ mod test { #[test] fn test_nedelec_2_hexahedron() { - let e = nedelec::create::(ReferenceCellType::Hexahedron, 2, Continuity::Standard); + let e = nedelec::create::(ReferenceCellType::Hexahedron, 2, Continuity::Standard); assert_eq!(e.value_size(), 3); check_dofs(e); } @@ -1792,7 +1842,7 @@ mod test { paste! { #[test] fn []() { - let e = lagrange::create::(ReferenceCellType::[<$cell>], [<$degree>], Continuity::Standard); + let e = lagrange::create::(ReferenceCellType::[<$cell>], [<$degree>], Continuity::Standard); let c = reference_cell::connectivity(ReferenceCellType::[<$cell>]); for (dim, entities) in c.iter().enumerate() { @@ -1841,7 +1891,7 @@ mod test { #[test] fn []() { - let e = [<$element>]::create::(ReferenceCellType::[<$cell>], [<$degree>], Continuity::Standard); + let e = [<$element>]::create::(ReferenceCellType::[<$cell>], [<$degree>], Continuity::Standard); let tdim = reference_cell::dim(ReferenceCellType::[<$cell>]); for edim in 1..tdim { for entity in &reference_cell::entity_types(ReferenceCellType::[<$cell>])[edim] { @@ -1945,7 +1995,7 @@ mod test { #[test] fn test_nedelec1_triangle_dof_transformations() { - let e = nedelec::create::(ReferenceCellType::Triangle, 1, Continuity::Standard); + let e = nedelec::create::(ReferenceCellType::Triangle, 1, Continuity::Standard); if let DofTransformation::Transformation(m, p) = e .dof_transformation(ReferenceCellType::Interval, Transformation::Reflection) .unwrap() @@ -1961,7 +2011,7 @@ mod test { #[test] fn test_nedelec2_triangle_dof_transformations() { - let e = nedelec::create::(ReferenceCellType::Triangle, 2, Continuity::Standard); + let e = nedelec::create::(ReferenceCellType::Triangle, 2, Continuity::Standard); if let DofTransformation::Transformation(m, p) = e .dof_transformation(ReferenceCellType::Interval, Transformation::Reflection) .unwrap() @@ -1977,7 +2027,8 @@ mod test { #[test] fn test_nedelec1_tetrahedron_dof_transformations() { - let e = nedelec::create::(ReferenceCellType::Tetrahedron, 1, Continuity::Standard); + let e = + nedelec::create::(ReferenceCellType::Tetrahedron, 1, Continuity::Standard); if let DofTransformation::Transformation(m, p) = e .dof_transformation(ReferenceCellType::Interval, Transformation::Reflection) .unwrap() @@ -2007,7 +2058,8 @@ mod test { #[test] fn test_nedelec2_tetrahedron_dof_transformations() { - let e = nedelec::create::(ReferenceCellType::Tetrahedron, 2, Continuity::Standard); + let e = + nedelec::create::(ReferenceCellType::Tetrahedron, 2, Continuity::Standard); if let DofTransformation::Transformation(m, p) = e .dof_transformation(ReferenceCellType::Interval, Transformation::Reflection) .unwrap() @@ -2046,7 +2098,7 @@ mod test { #[test] fn test_nedelec2_quadrilateral_dof_transformations() { - let e = nedelec::create::(ReferenceCellType::Hexahedron, 2, Continuity::Standard); + let e = nedelec::create::(ReferenceCellType::Hexahedron, 2, Continuity::Standard); if let DofTransformation::Transformation(m, p) = e .dof_transformation(ReferenceCellType::Interval, Transformation::Reflection) .unwrap() @@ -2062,7 +2114,7 @@ mod test { #[test] fn test_nedelec2_hexahedron_dof_transformations() { - let e = nedelec::create::(ReferenceCellType::Hexahedron, 2, Continuity::Standard); + let e = nedelec::create::(ReferenceCellType::Hexahedron, 2, Continuity::Standard); if let DofTransformation::Transformation(m, p) = e .dof_transformation(ReferenceCellType::Interval, Transformation::Reflection) .unwrap() @@ -2112,7 +2164,8 @@ mod test { #[test] fn test_lagrange4_tetrahedron_dof_transformations() { - let e = lagrange::create::(ReferenceCellType::Tetrahedron, 4, Continuity::Standard); + let e = + lagrange::create::(ReferenceCellType::Tetrahedron, 4, Continuity::Standard); if let DofTransformation::Permutation(p) = e .dof_transformation(ReferenceCellType::Interval, Transformation::Reflection) .unwrap() @@ -2156,7 +2209,7 @@ mod test { #[test] fn test_dof_permuting_triangle() { - let e = lagrange::create::(ReferenceCellType::Triangle, 3, Continuity::Standard); + let e = lagrange::create::(ReferenceCellType::Triangle, 3, Continuity::Standard); let mut n = (0..10).collect::>(); e.apply_dof_permutations(&mut n, 7); @@ -2178,7 +2231,8 @@ mod test { #[test] fn test_dof_permuting_tetrahedron() { - let e = lagrange::create::(ReferenceCellType::Tetrahedron, 3, Continuity::Standard); + let e = + lagrange::create::(ReferenceCellType::Tetrahedron, 3, Continuity::Standard); let mut n = (0..20).collect::>(); e.apply_dof_permutations(&mut n, 63); diff --git a/src/ciarlet/lagrange.rs b/src/ciarlet/lagrange.rs index 0879239..dc9eed1 100644 --- a/src/ciarlet/lagrange.rs +++ b/src/ciarlet/lagrange.rs @@ -11,11 +11,11 @@ use rlst::{RlstScalar, rlst_dynamic_array}; use std::marker::PhantomData; /// Create a Lagrange element. -pub fn create( +pub fn create( cell_type: ReferenceCellType, degree: usize, continuity: Continuity, -) -> CiarletElement { +) -> CiarletElement { let dim = polynomial_count(cell_type, degree); let tdim = reference_cell::dim(cell_type); let mut wcoeffs = rlst_dynamic_array!(T, [dim, 1, dim]); @@ -26,24 +26,24 @@ pub fn create( let mut x = [vec![], vec![], vec![], vec![]]; let mut m = [vec![], vec![], vec![], vec![]]; let entity_counts = reference_cell::entity_counts(cell_type); - let vertices = reference_cell::vertices::(cell_type); + let vertices = reference_cell::vertices::(cell_type); if degree == 0 { if continuity == Continuity::Standard { panic!("Cannot create continuous degree 0 Lagrange element"); } for (d, counts) in entity_counts.iter().enumerate() { for _e in 0..*counts { - x[d].push(rlst_dynamic_array!(T::Real, [tdim, 0])); + x[d].push(rlst_dynamic_array!(TGeo, [tdim, 0])); m[d].push(rlst_dynamic_array!(T, [0, 1, 0])); } } - let mut midp = rlst_dynamic_array!(T::Real, [tdim, 1]); + let mut midp = rlst_dynamic_array!(TGeo, [tdim, 1]); let nvertices = entity_counts[0]; for i in 0..tdim { for vertex in &vertices { - *midp.get_mut([i, 0]).unwrap() += num::cast::<_, T::Real>(vertex[i]).unwrap(); + *midp.get_mut([i, 0]).unwrap() += num::cast::<_, TGeo>(vertex[i]).unwrap(); } - *midp.get_mut([i, 0]).unwrap() /= num::cast::<_, T::Real>(nvertices).unwrap(); + *midp.get_mut([i, 0]).unwrap() /= num::cast::<_, TGeo>(nvertices).unwrap(); } x[tdim].push(midp); let mut mentry = rlst_dynamic_array!(T, [1, 1, 1]); @@ -54,9 +54,9 @@ pub fn create( let faces = reference_cell::faces(cell_type); let volumes = reference_cell::volumes(cell_type); for vertex in &vertices { - let mut pts = rlst_dynamic_array!(T::Real, [tdim, 1]); + let mut pts = rlst_dynamic_array!(TGeo, [tdim, 1]); for (i, v) in vertex.iter().enumerate() { - *pts.get_mut([i, 0]).unwrap() = num::cast::<_, T::Real>(*v).unwrap(); + *pts.get_mut([i, 0]).unwrap() = num::cast::<_, TGeo>(*v).unwrap(); } x[0].push(pts); let mut mentry = rlst_dynamic_array!(T, [1, 1, 1]); @@ -64,7 +64,7 @@ pub fn create( m[0].push(mentry); } for e in &edges { - let mut pts = rlst_dynamic_array!(T::Real, [tdim, degree - 1]); + let mut pts = rlst_dynamic_array!(TGeo, [tdim, degree - 1]); let [vn0, vn1] = e[..] else { panic!(); }; @@ -75,10 +75,9 @@ pub fn create( for i in 1..degree { *ident.get_mut([i - 1, 0, i - 1]).unwrap() = T::from(1.0).unwrap(); for j in 0..tdim { - *pts.get_mut([j, i - 1]).unwrap() = num::cast::<_, T::Real>(v0[j]).unwrap() - + num::cast::<_, T::Real>(i).unwrap() - / num::cast::<_, T::Real>(degree).unwrap() - * num::cast::<_, T::Real>(v1[j] - v0[j]).unwrap(); + *pts.get_mut([j, i - 1]).unwrap() = num::cast::<_, TGeo>(v0[j]).unwrap() + + num::cast::<_, TGeo>(i).unwrap() / num::cast::<_, TGeo>(degree).unwrap() + * num::cast::<_, TGeo>(v1[j] - v0[j]).unwrap(); } } x[1].push(pts); @@ -101,7 +100,7 @@ pub fn create( panic!("Unsupported face type"); } }; - let mut pts = rlst_dynamic_array!(T::Real, [tdim, npts]); + let mut pts = rlst_dynamic_array!(TGeo, [tdim, npts]); let [vn0, vn1, vn2] = faces[e][..3] else { panic!(); @@ -116,14 +115,14 @@ pub fn create( for i0 in 1..degree { for i1 in 1..degree - i0 { for j in 0..tdim { - *pts.get_mut([j, n]).unwrap() = num::cast::<_, T::Real>(v0[j]) + *pts.get_mut([j, n]).unwrap() = num::cast::<_, TGeo>(v0[j]) .unwrap() - + num::cast::<_, T::Real>(i0).unwrap() - / num::cast::<_, T::Real>(degree).unwrap() - * num::cast::<_, T::Real>(v1[j] - v0[j]).unwrap() - + num::cast::<_, T::Real>(i1).unwrap() - / num::cast::<_, T::Real>(degree).unwrap() - * num::cast::<_, T::Real>(v2[j] - v0[j]).unwrap(); + + num::cast::<_, TGeo>(i0).unwrap() + / num::cast::<_, TGeo>(degree).unwrap() + * num::cast::<_, TGeo>(v1[j] - v0[j]).unwrap() + + num::cast::<_, TGeo>(i1).unwrap() + / num::cast::<_, TGeo>(degree).unwrap() + * num::cast::<_, TGeo>(v2[j] - v0[j]).unwrap(); } n += 1; } @@ -134,14 +133,14 @@ pub fn create( for i0 in 1..degree { for i1 in 1..degree { for j in 0..tdim { - *pts.get_mut([j, n]).unwrap() = num::cast::<_, T::Real>(v0[j]) + *pts.get_mut([j, n]).unwrap() = num::cast::<_, TGeo>(v0[j]) .unwrap() - + num::cast::<_, T::Real>(i0).unwrap() - / num::cast::<_, T::Real>(degree).unwrap() - * num::cast::<_, T::Real>(v1[j] - v0[j]).unwrap() - + num::cast::<_, T::Real>(i1).unwrap() - / num::cast::<_, T::Real>(degree).unwrap() - * num::cast::<_, T::Real>(v2[j] - v0[j]).unwrap(); + + num::cast::<_, TGeo>(i0).unwrap() + / num::cast::<_, TGeo>(degree).unwrap() + * num::cast::<_, TGeo>(v1[j] - v0[j]).unwrap() + + num::cast::<_, TGeo>(i1).unwrap() + / num::cast::<_, TGeo>(degree).unwrap() + * num::cast::<_, TGeo>(v2[j] - v0[j]).unwrap(); } n += 1; } @@ -176,7 +175,7 @@ pub fn create( panic!("Unsupported face type"); } }; - let mut pts = rlst_dynamic_array!(T::Real, [tdim, npts]); + let mut pts = rlst_dynamic_array!(TGeo, [tdim, npts]); match volume_type { ReferenceCellType::Tetrahedron => { @@ -193,17 +192,17 @@ pub fn create( for i1 in 1..degree - i0 { for i2 in 1..degree - i0 - i1 { for j in 0..tdim { - *pts.get_mut([j, n]).unwrap() = num::cast::<_, T::Real>(v0[j]) + *pts.get_mut([j, n]).unwrap() = num::cast::<_, TGeo>(v0[j]) .unwrap() - + num::cast::<_, T::Real>(i0).unwrap() - / num::cast::<_, T::Real>(degree).unwrap() - * num::cast::<_, T::Real>(v1[j] - v0[j]).unwrap() - + num::cast::<_, T::Real>(i1).unwrap() - / num::cast::<_, T::Real>(degree).unwrap() - * num::cast::<_, T::Real>(v2[j] - v0[j]).unwrap() - + num::cast::<_, T::Real>(i2).unwrap() - / num::cast::<_, T::Real>(degree).unwrap() - * num::cast::<_, T::Real>(v3[j] - v0[j]).unwrap(); + + num::cast::<_, TGeo>(i0).unwrap() + / num::cast::<_, TGeo>(degree).unwrap() + * num::cast::<_, TGeo>(v1[j] - v0[j]).unwrap() + + num::cast::<_, TGeo>(i1).unwrap() + / num::cast::<_, TGeo>(degree).unwrap() + * num::cast::<_, TGeo>(v2[j] - v0[j]).unwrap() + + num::cast::<_, TGeo>(i2).unwrap() + / num::cast::<_, TGeo>(degree).unwrap() + * num::cast::<_, TGeo>(v3[j] - v0[j]).unwrap(); } n += 1; } @@ -224,17 +223,17 @@ pub fn create( for i1 in 1..degree { for i2 in 1..degree { for j in 0..tdim { - *pts.get_mut([j, n]).unwrap() = num::cast::<_, T::Real>(v0[j]) + *pts.get_mut([j, n]).unwrap() = num::cast::<_, TGeo>(v0[j]) .unwrap() - + num::cast::<_, T::Real>(i0).unwrap() - / num::cast::<_, T::Real>(degree).unwrap() - * num::cast::<_, T::Real>(v1[j] - v0[j]).unwrap() - + num::cast::<_, T::Real>(i1).unwrap() - / num::cast::<_, T::Real>(degree).unwrap() - * num::cast::<_, T::Real>(v2[j] - v0[j]).unwrap() - + num::cast::<_, T::Real>(i2).unwrap() - / num::cast::<_, T::Real>(degree).unwrap() - * num::cast::<_, T::Real>(v3[j] - v0[j]).unwrap(); + + num::cast::<_, TGeo>(i0).unwrap() + / num::cast::<_, TGeo>(degree).unwrap() + * num::cast::<_, TGeo>(v1[j] - v0[j]).unwrap() + + num::cast::<_, TGeo>(i1).unwrap() + / num::cast::<_, TGeo>(degree).unwrap() + * num::cast::<_, TGeo>(v2[j] - v0[j]).unwrap() + + num::cast::<_, TGeo>(i2).unwrap() + / num::cast::<_, TGeo>(degree).unwrap() + * num::cast::<_, TGeo>(v3[j] - v0[j]).unwrap(); } n += 1; } @@ -254,7 +253,7 @@ pub fn create( m[3].push(ident); } } - CiarletElement::::create( + CiarletElement::::create( "Lagrange".to_string(), cell_type, degree, @@ -272,28 +271,32 @@ pub fn create( /// /// A family of Lagrange elements on multiple cell types with appropriate /// continuity across different cell types. -pub struct LagrangeElementFamily { +pub struct LagrangeElementFamily { degree: usize, continuity: Continuity, _t: PhantomData, + _tgeo: PhantomData, } -impl LagrangeElementFamily { +impl LagrangeElementFamily { /// Create new family with given `degree` and `continuity`. pub fn new(degree: usize, continuity: Continuity) -> Self { Self { degree, continuity, _t: PhantomData, + _tgeo: PhantomData, } } } -impl ElementFamily for LagrangeElementFamily { +impl ElementFamily + for LagrangeElementFamily +{ type T = T; - type FiniteElement = CiarletElement; + type FiniteElement = CiarletElement; type CellType = ReferenceCellType; - fn element(&self, cell_type: ReferenceCellType) -> CiarletElement { - create::(cell_type, self.degree, self.continuity) + fn element(&self, cell_type: ReferenceCellType) -> CiarletElement { + create::(cell_type, self.degree, self.continuity) } } diff --git a/src/ciarlet/nedelec.rs b/src/ciarlet/nedelec.rs index 0aa2567..2c31a54 100644 --- a/src/ciarlet/nedelec.rs +++ b/src/ciarlet/nedelec.rs @@ -13,11 +13,11 @@ use rlst::dense::linalg::lapack::interface::{getrf::Getrf, getri::Getri}; use rlst::{DynArray, RlstScalar, SliceArray, rlst_dynamic_array}; use std::marker::PhantomData; -fn create_simplex, T: RlstScalar + Getrf + Getri>( +fn create_simplex( cell_type: ReferenceCellType, degree: usize, continuity: Continuity, -) -> CiarletElement { +) -> CiarletElement { if cell_type != ReferenceCellType::Triangle && cell_type != ReferenceCellType::Tetrahedron { panic!("Invalid cell: {cell_type:?}"); } @@ -31,9 +31,9 @@ fn create_simplex, T: RlstScalar + let pts_t = cell_q .points .iter() - .map(|i| TReal::from(*i).unwrap()) + .map(|i| TGeo::from(*i).unwrap()) .collect::>(); - let pts = SliceArray::::from_shape(&pts_t, [tdim, cell_q.npoints]); + let pts = SliceArray::::from_shape(&pts_t, [tdim, cell_q.npoints]); let mut phi = DynArray::::from_shape(legendre_shape(cell_type, &pts, degree, 0)); tabulate_legendre_polynomials(cell_type, &pts, degree, 0, &mut phi); @@ -178,10 +178,10 @@ fn create_simplex, T: RlstScalar + let mut m = [vec![], vec![], vec![], vec![]]; let entity_counts = reference_cell::entity_counts(cell_type); - let vertices = reference_cell::vertices::(cell_type); + let vertices = reference_cell::vertices::(cell_type); for _ in 0..entity_counts[0] { - x[0].push(rlst_dynamic_array!(T::Real, [tdim, 0])); + x[0].push(rlst_dynamic_array!(TGeo, [tdim, 0])); m[0].push(rlst_dynamic_array!(T, [0, tdim, 0])); } @@ -190,9 +190,9 @@ fn create_simplex, T: RlstScalar + let edge_pts_t = edge_q .points .iter() - .map(|i| TReal::from(*i).unwrap()) + .map(|i| TGeo::from(*i).unwrap()) .collect::>(); - let edge_pts = SliceArray::::from_shape(&edge_pts_t, [1, edge_q.npoints]); + let edge_pts = SliceArray::::from_shape(&edge_pts_t, [1, edge_q.npoints]); let mut edge_phi = DynArray::::from_shape(legendre_shape( ReferenceCellType::Interval, @@ -209,7 +209,7 @@ fn create_simplex, T: RlstScalar + ); for edge in reference_cell::edges(cell_type) { - let mut pts = rlst_dynamic_array!(T::Real, [tdim, edge_q.npoints]); + let mut pts = rlst_dynamic_array!(TGeo, [tdim, edge_q.npoints]); let mut mat = rlst_dynamic_array!(T, [pdim_edge_minus1, tdim, edge_q.npoints]); for (w_i, (pt, wt)) in izip!(&edge_pts_t, &edge_q.weights).enumerate() { @@ -232,7 +232,7 @@ fn create_simplex, T: RlstScalar + // DOFs on faces if degree == 1 { for _ in 0..entity_counts[2] { - x[2].push(rlst_dynamic_array!(T::Real, [tdim, 0])); + x[2].push(rlst_dynamic_array!(TGeo, [tdim, 0])); m[2].push(rlst_dynamic_array!(T, [0, tdim, 0])) } } else { @@ -240,9 +240,9 @@ fn create_simplex, T: RlstScalar + let face_pts_t = face_q .points .iter() - .map(|i| TReal::from(*i).unwrap()) + .map(|i| TGeo::from(*i).unwrap()) .collect::>(); - let face_pts = SliceArray::::from_shape(&face_pts_t, [2, face_q.npoints]); + let face_pts = SliceArray::::from_shape(&face_pts_t, [2, face_q.npoints]); let mut face_phi = DynArray::::from_shape(legendre_shape( ReferenceCellType::Triangle, @@ -259,7 +259,7 @@ fn create_simplex, T: RlstScalar + ); for face in reference_cell::faces(cell_type) { - let mut pts = rlst_dynamic_array!(T::Real, [tdim, face_q.npoints]); + let mut pts = rlst_dynamic_array!(TGeo, [tdim, face_q.npoints]); let mut mat = rlst_dynamic_array!(T, [2 * pdim_face_minus2, tdim, face_q.npoints]); for (w_i, wt) in face_q.weights.iter().enumerate() { @@ -284,14 +284,14 @@ fn create_simplex, T: RlstScalar + } if tdim == 3 { if degree <= 2 { - x[3].push(rlst_dynamic_array!(T::Real, [tdim, 0])); + x[3].push(rlst_dynamic_array!(TGeo, [tdim, 0])); m[3].push(rlst_dynamic_array!(T, [0, tdim, 0])) } else { let internal_q = gauss_jacobi_rule(cell_type, 2 * degree - 3).unwrap(); - let mut pts = rlst_dynamic_array!(T::Real, [tdim, internal_q.npoints]); + let mut pts = rlst_dynamic_array!(TGeo, [tdim, internal_q.npoints]); for p in 0..internal_q.npoints { for d in 0..3 { - pts[[d, p]] = TReal::from(internal_q.points[3 * p + d]).unwrap() + pts[[d, p]] = TGeo::from(internal_q.points[3 * p + d]).unwrap() } } @@ -328,11 +328,11 @@ fn create_simplex, T: RlstScalar + ) } -fn create_tp, T: RlstScalar + Getrf + Getri>( +fn create_tp( cell_type: ReferenceCellType, degree: usize, continuity: Continuity, -) -> CiarletElement { +) -> CiarletElement { if cell_type != ReferenceCellType::Quadrilateral && cell_type != ReferenceCellType::Hexahedron { panic!("Invalid cell: {cell_type:?}"); } @@ -346,9 +346,9 @@ fn create_tp, T: RlstScalar + Getr let pts_t = cell_q .points .iter() - .map(|i| TReal::from(*i).unwrap()) + .map(|i| TGeo::from(*i).unwrap()) .collect::>(); - let pts = SliceArray::::from_shape(&pts_t, [tdim, cell_q.npoints]); + let pts = SliceArray::::from_shape(&pts_t, [tdim, cell_q.npoints]); let mut phi = DynArray::::from_shape(legendre_shape(cell_type, &pts, degree, 0)); tabulate_legendre_polynomials(cell_type, &pts, degree, 0, &mut phi); @@ -425,10 +425,10 @@ fn create_tp, T: RlstScalar + Getr let mut x = [vec![], vec![], vec![], vec![]]; let mut m = [vec![], vec![], vec![], vec![]]; - let vertices = reference_cell::vertices::(cell_type); + let vertices = reference_cell::vertices::(cell_type); for _ in 0..entity_counts[0] { - x[0].push(rlst_dynamic_array!(T::Real, [tdim, 0])); + x[0].push(rlst_dynamic_array!(TGeo, [tdim, 0])); m[0].push(rlst_dynamic_array!(T, [0, tdim, 0])); } @@ -437,9 +437,9 @@ fn create_tp, T: RlstScalar + Getr let edge_pts_t = edge_q .points .iter() - .map(|i| TReal::from(*i).unwrap()) + .map(|i| TGeo::from(*i).unwrap()) .collect::>(); - let edge_pts = SliceArray::::from_shape(&edge_pts_t, [1, edge_q.npoints]); + let edge_pts = SliceArray::::from_shape(&edge_pts_t, [1, edge_q.npoints]); let mut edge_phi = DynArray::::from_shape(legendre_shape( ReferenceCellType::Interval, @@ -456,7 +456,7 @@ fn create_tp, T: RlstScalar + Getr ); for edge in reference_cell::edges(cell_type) { - let mut pts = rlst_dynamic_array!(T::Real, [tdim, edge_q.npoints]); + let mut pts = rlst_dynamic_array!(TGeo, [tdim, edge_q.npoints]); let mut mat = rlst_dynamic_array!(T, [pdim_edge_minus1, tdim, edge_q.npoints]); for (w_i, (pt, wt)) in izip!(&edge_pts_t, &edge_q.weights).enumerate() { @@ -479,7 +479,7 @@ fn create_tp, T: RlstScalar + Getr // DOFs on faces if degree == 1 { for _ in 0..entity_counts[2] { - x[2].push(rlst_dynamic_array!(T::Real, [tdim, 0])); + x[2].push(rlst_dynamic_array!(TGeo, [tdim, 0])); m[2].push(rlst_dynamic_array!(T, [0, tdim, 0])) } } else { @@ -487,9 +487,9 @@ fn create_tp, T: RlstScalar + Getr let face_pts_t = face_q .points .iter() - .map(|i| TReal::from(*i).unwrap()) + .map(|i| TGeo::from(*i).unwrap()) .collect::>(); - let face_pts = SliceArray::::from_shape(&face_pts_t, [2, face_q.npoints]); + let face_pts = SliceArray::::from_shape(&face_pts_t, [2, face_q.npoints]); let mut face_phi = DynArray::::from_shape(legendre_shape( ReferenceCellType::Quadrilateral, @@ -506,7 +506,7 @@ fn create_tp, T: RlstScalar + Getr ); for face in reference_cell::faces(cell_type) { - let mut pts = rlst_dynamic_array!(T::Real, [tdim, face_q.npoints]); + let mut pts = rlst_dynamic_array!(TGeo, [tdim, face_q.npoints]); let mdim = 2 * pdim_edge_minus2 * pdim_edge_minus1; let mut mat = rlst_dynamic_array!(T, [mdim, tdim, face_q.npoints]); @@ -539,7 +539,7 @@ fn create_tp, T: RlstScalar + Getr // DOFs on volume if tdim == 3 { if degree == 1 { - x[3].push(rlst_dynamic_array!(T::Real, [tdim, 0])); + x[3].push(rlst_dynamic_array!(TGeo, [tdim, 0])); m[3].push(rlst_dynamic_array!(T, [0, tdim, 0])) } else { let interior_q = @@ -547,10 +547,10 @@ fn create_tp, T: RlstScalar + Getr let interior_pts_t = interior_q .points .iter() - .map(|i| TReal::from(*i).unwrap()) + .map(|i| TGeo::from(*i).unwrap()) .collect::>(); let interior_pts = - SliceArray::::from_shape(&interior_pts_t, [3, interior_q.npoints]); + SliceArray::::from_shape(&interior_pts_t, [3, interior_q.npoints]); let mut interior_phi = DynArray::::from_shape(legendre_shape( ReferenceCellType::Hexahedron, @@ -566,7 +566,7 @@ fn create_tp, T: RlstScalar + Getr &mut interior_phi, ); - let mut pts = rlst_dynamic_array!(T::Real, [tdim, interior_q.npoints]); + let mut pts = rlst_dynamic_array!(TGeo, [tdim, interior_q.npoints]); let mdim = 3 * pdim_edge_minus2.pow(2) * pdim_edge_minus1; let mut mat = rlst_dynamic_array!(T, [mdim, tdim, interior_q.npoints]); @@ -623,17 +623,17 @@ fn create_tp, T: RlstScalar + Getr } /// Create a Nedelec (first kind) element -pub fn create( +pub fn create( cell_type: ReferenceCellType, degree: usize, continuity: Continuity, -) -> CiarletElement { +) -> CiarletElement { if cell_type == ReferenceCellType::Triangle || cell_type == ReferenceCellType::Tetrahedron { - create_simplex(cell_type, degree, continuity) + create_simplex::(cell_type, degree, continuity) } else if cell_type == ReferenceCellType::Quadrilateral || cell_type == ReferenceCellType::Hexahedron { - create_tp(cell_type, degree, continuity) + create_tp::(cell_type, degree, continuity) } else { panic!("Invalid cell: {cell_type:?}"); } @@ -643,28 +643,35 @@ pub fn create( /// /// A family of Nedelec elements on multiple cell types with appropriate /// continuity across different cell types. -pub struct NedelecFirstKindElementFamily { +pub struct NedelecFirstKindElementFamily< + T: RlstScalar + Getrf + Getri = f64, + TGeo: RlstScalar = f64, +> { degree: usize, continuity: Continuity, _t: PhantomData, + _tgeo: PhantomData, } -impl NedelecFirstKindElementFamily { +impl NedelecFirstKindElementFamily { /// Create new family with given `degree` and `continuity`. pub fn new(degree: usize, continuity: Continuity) -> Self { Self { degree, continuity, _t: PhantomData, + _tgeo: PhantomData, } } } -impl ElementFamily for NedelecFirstKindElementFamily { +impl ElementFamily + for NedelecFirstKindElementFamily +{ type T = T; type CellType = ReferenceCellType; - type FiniteElement = CiarletElement; - fn element(&self, cell_type: ReferenceCellType) -> CiarletElement { - create::(cell_type, self.degree, self.continuity) + type FiniteElement = CiarletElement; + fn element(&self, cell_type: ReferenceCellType) -> CiarletElement { + create::(cell_type, self.degree, self.continuity) } } diff --git a/src/ciarlet/raviart_thomas.rs b/src/ciarlet/raviart_thomas.rs index 4fb707a..b31d8ad 100644 --- a/src/ciarlet/raviart_thomas.rs +++ b/src/ciarlet/raviart_thomas.rs @@ -13,11 +13,11 @@ use rlst::dense::linalg::lapack::interface::{getrf::Getrf, getri::Getri}; use rlst::{DynArray, RlstScalar, SliceArray, rlst_dynamic_array}; use std::marker::PhantomData; -fn create_simplex, T: RlstScalar + Getrf + Getri>( +fn create_simplex( cell_type: ReferenceCellType, degree: usize, continuity: Continuity, -) -> CiarletElement { +) -> CiarletElement { if cell_type != ReferenceCellType::Triangle && cell_type != ReferenceCellType::Tetrahedron { panic!("Invalid cell: {cell_type:?}"); } @@ -44,9 +44,9 @@ fn create_simplex, T: RlstScalar + let pts_t = cell_q .points .iter() - .map(|i| TReal::from(*i).unwrap()) + .map(|i| TGeo::from(*i).unwrap()) .collect::>(); - let pts = SliceArray::::from_shape(&pts_t, [tdim, cell_q.npoints]); + let pts = SliceArray::::from_shape(&pts_t, [tdim, cell_q.npoints]); let mut phi = DynArray::::from_shape(legendre_shape(cell_type, &pts, degree, 0)); tabulate_legendre_polynomials(cell_type, &pts, degree, 0, &mut phi); @@ -87,11 +87,11 @@ fn create_simplex, T: RlstScalar + let mut m = [vec![], vec![], vec![], vec![]]; let entity_counts = reference_cell::entity_counts(cell_type); - let vertices = reference_cell::vertices::(cell_type); + let vertices = reference_cell::vertices::(cell_type); for d in 0..tdim - 1 { for _ in 0..entity_counts[d] { - x[d].push(rlst_dynamic_array!(T::Real, [tdim, 0])); + x[d].push(rlst_dynamic_array!(TGeo, [tdim, 0])); m[d].push(rlst_dynamic_array!(T, [0, tdim, 0])); } } @@ -101,9 +101,9 @@ fn create_simplex, T: RlstScalar + let facet_pts_t = facet_q .points .iter() - .map(|i| TReal::from(*i).unwrap()) + .map(|i| TGeo::from(*i).unwrap()) .collect::>(); - let facet_pts = SliceArray::::from_shape(&facet_pts_t, [tdim - 1, facet_q.npoints]); + let facet_pts = SliceArray::::from_shape(&facet_pts_t, [tdim - 1, facet_q.npoints]); let mut facet_phi = DynArray::::from_shape(legendre_shape(facet_type, &facet_pts, degree - 1, 0)); @@ -115,9 +115,9 @@ fn create_simplex, T: RlstScalar + } else { reference_cell::faces(cell_type) }, - reference_cell::facet_normals::(cell_type), + reference_cell::facet_normals::(cell_type), ) { - let mut pts = rlst_dynamic_array!(T::Real, [tdim, facet_q.npoints]); + let mut pts = rlst_dynamic_array!(TGeo, [tdim, facet_q.npoints]); let mut mat = rlst_dynamic_array!(T, [pdim_facet_minus1, tdim, facet_q.npoints]); for (w_i, wt) in facet_q.weights.iter().enumerate() { @@ -128,7 +128,7 @@ fn create_simplex, T: RlstScalar + &facet_q.points[w_i * (tdim - 1)..(w_i + 1) * (tdim - 1)] ) .map(|(v_i, pt)| { - (vertices[*v_i][i] - vertices[facet[0]][i]) * TReal::from(*pt).unwrap() + (vertices[*v_i][i] - vertices[facet[0]][i]) * TGeo::from(*pt).unwrap() }) .sum(); @@ -146,15 +146,15 @@ fn create_simplex, T: RlstScalar + if degree == 1 { for _ in 0..entity_counts[tdim] { - x[tdim].push(rlst_dynamic_array!(T::Real, [tdim, 0])); + x[tdim].push(rlst_dynamic_array!(TGeo, [tdim, 0])); m[tdim].push(rlst_dynamic_array!(T, [0, tdim, 0])) } } else { let internal_q = gauss_jacobi_rule(cell_type, 2 * degree - 2).unwrap(); - let mut pts = rlst_dynamic_array!(T::Real, [tdim, internal_q.npoints]); + let mut pts = rlst_dynamic_array!(TGeo, [tdim, internal_q.npoints]); for p in 0..internal_q.npoints { for d in 0..tdim { - pts[[d, p]] = TReal::from(internal_q.points[tdim * p + d]).unwrap() + pts[[d, p]] = TGeo::from(internal_q.points[tdim * p + d]).unwrap() } } @@ -190,11 +190,11 @@ fn create_simplex, T: RlstScalar + ) } -fn create_tp, T: RlstScalar + Getrf + Getri>( +fn create_tp( cell_type: ReferenceCellType, degree: usize, continuity: Continuity, -) -> CiarletElement { +) -> CiarletElement { if cell_type != ReferenceCellType::Quadrilateral && cell_type != ReferenceCellType::Hexahedron { panic!("Invalid cell: {cell_type:?}"); } @@ -208,9 +208,9 @@ fn create_tp, T: RlstScalar + Getr let pts_t = cell_q .points .iter() - .map(|i| TReal::from(*i).unwrap()) + .map(|i| TGeo::from(*i).unwrap()) .collect::>(); - let pts = SliceArray::::from_shape(&pts_t, [tdim, cell_q.npoints]); + let pts = SliceArray::::from_shape(&pts_t, [tdim, cell_q.npoints]); let mut phi = DynArray::::from_shape(legendre_shape(cell_type, &pts, degree, 0)); tabulate_legendre_polynomials(cell_type, &pts, degree, 0, &mut phi); @@ -293,11 +293,11 @@ fn create_tp, T: RlstScalar + Getr let mut x = [vec![], vec![], vec![], vec![]]; let mut m = [vec![], vec![], vec![], vec![]]; - let vertices = reference_cell::vertices::(cell_type); + let vertices = reference_cell::vertices::(cell_type); for d in 0..tdim - 1 { for _ in 0..entity_counts[d] { - x[d].push(rlst_dynamic_array!(T::Real, [tdim, 0])); + x[d].push(rlst_dynamic_array!(TGeo, [tdim, 0])); m[d].push(rlst_dynamic_array!(T, [0, tdim, 0])); } } @@ -306,9 +306,9 @@ fn create_tp, T: RlstScalar + Getr let facet_pts_t = facet_q .points .iter() - .map(|i| TReal::from(*i).unwrap()) + .map(|i| TGeo::from(*i).unwrap()) .collect::>(); - let facet_pts = SliceArray::::from_shape(&facet_pts_t, [tdim - 1, facet_q.npoints]); + let facet_pts = SliceArray::::from_shape(&facet_pts_t, [tdim - 1, facet_q.npoints]); let mut facet_phi = DynArray::::from_shape(legendre_shape(facet_type, &facet_pts, degree - 1, 0)); @@ -316,9 +316,9 @@ fn create_tp, T: RlstScalar + Getr for (facet, normal) in izip!( reference_cell::facets(cell_type), - reference_cell::facet_normals::(cell_type), + reference_cell::facet_normals::(cell_type), ) { - let mut pts = rlst_dynamic_array!(T::Real, [tdim, facet_q.npoints]); + let mut pts = rlst_dynamic_array!(TGeo, [tdim, facet_q.npoints]); let mut mat = rlst_dynamic_array!(T, [pdim_facet_minus1, tdim, facet_q.npoints]); for (w_i, wt) in facet_q.weights.iter().enumerate() { @@ -329,7 +329,7 @@ fn create_tp, T: RlstScalar + Getr (vertices[facet[usize::pow(2, x as u32)]][d] - vertices[facet[0]][d]) * facet_pts[[x, w_i]] }) - .sum::(); + .sum::(); for j in 0..facet_phi.shape()[1] { mat[[j, d, w_i]] = T::from(*wt).unwrap() * facet_phi[[0, j, w_i]] @@ -344,16 +344,16 @@ fn create_tp, T: RlstScalar + Getr // DOFs on interior if degree == 1 { - x[tdim].push(rlst_dynamic_array!(T::Real, [tdim, 0])); + x[tdim].push(rlst_dynamic_array!(TGeo, [tdim, 0])); m[tdim].push(rlst_dynamic_array!(T, [0, tdim, 0])) } else if tdim == 2 { let face_q = gauss_jacobi_rule(ReferenceCellType::Quadrilateral, 2 * degree - 1).unwrap(); let face_pts_t = face_q .points .iter() - .map(|i| TReal::from(*i).unwrap()) + .map(|i| TGeo::from(*i).unwrap()) .collect::>(); - let face_pts = SliceArray::::from_shape(&face_pts_t, [2, face_q.npoints]); + let face_pts = SliceArray::::from_shape(&face_pts_t, [2, face_q.npoints]); let mut face_phi = DynArray::::from_shape(legendre_shape( ReferenceCellType::Quadrilateral, @@ -369,7 +369,7 @@ fn create_tp, T: RlstScalar + Getr &mut face_phi, ); - let mut pts = rlst_dynamic_array!(T::Real, [tdim, face_q.npoints]); + let mut pts = rlst_dynamic_array!(TGeo, [tdim, face_q.npoints]); let mdim = 2 * pdim_edge_minus2 * pdim_edge_minus1; let mut mat = rlst_dynamic_array!(T, [mdim, tdim, face_q.npoints]); @@ -394,10 +394,10 @@ fn create_tp, T: RlstScalar + Getr let interior_pts_t = interior_q .points .iter() - .map(|i| TReal::from(*i).unwrap()) + .map(|i| TGeo::from(*i).unwrap()) .collect::>(); let interior_pts = - SliceArray::::from_shape(&interior_pts_t, [3, interior_q.npoints]); + SliceArray::::from_shape(&interior_pts_t, [3, interior_q.npoints]); let mut interior_phi = DynArray::::from_shape(legendre_shape( ReferenceCellType::Hexahedron, @@ -413,7 +413,7 @@ fn create_tp, T: RlstScalar + Getr &mut interior_phi, ); - let mut pts = rlst_dynamic_array!(T::Real, [tdim, interior_q.npoints]); + let mut pts = rlst_dynamic_array!(TGeo, [tdim, interior_q.npoints]); let mdim = 3 * pdim_edge_minus1.pow(2) * pdim_edge_minus2; let mut mat = rlst_dynamic_array!(T, [mdim, tdim, interior_q.npoints]); @@ -466,11 +466,11 @@ fn create_tp, T: RlstScalar + Getr } /// Create a Raviart-Thomas element -pub fn create( +pub fn create( cell_type: ReferenceCellType, degree: usize, continuity: Continuity, -) -> CiarletElement { +) -> CiarletElement { if cell_type == ReferenceCellType::Triangle || cell_type == ReferenceCellType::Tetrahedron { create_simplex(cell_type, degree, continuity) } else if cell_type == ReferenceCellType::Quadrilateral @@ -486,28 +486,35 @@ pub fn create( /// /// A family of Raviart-Thomas elements on multiple cell types with appropriate /// continuity across different cell types. -pub struct RaviartThomasElementFamily { +pub struct RaviartThomasElementFamily { degree: usize, continuity: Continuity, _t: PhantomData, + _tgeo: PhantomData, } -impl RaviartThomasElementFamily { +impl RaviartThomasElementFamily { /// Create new family with given `degree` and `continuity`. pub fn new(degree: usize, continuity: Continuity) -> Self { Self { degree, continuity, _t: PhantomData, + _tgeo: PhantomData, } } } -impl ElementFamily for RaviartThomasElementFamily { +impl ElementFamily + for RaviartThomasElementFamily +{ type T = T; type CellType = ReferenceCellType; - type FiniteElement = CiarletElement; - fn element(&self, cell_type: ReferenceCellType) -> CiarletElement { - create::(cell_type, self.degree, self.continuity) + type FiniteElement = CiarletElement; + fn element( + &self, + cell_type: ReferenceCellType, + ) -> CiarletElement { + create::(cell_type, self.degree, self.continuity) } } From cd9ae7aff65fc68f8c2e1ac03eb024c3ab64df29 Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Tue, 25 Nov 2025 13:39:56 +0000 Subject: [PATCH 7/9] clippy --- src/bindings.rs | 132 +++++++++++++++++++++++++++--------------------- 1 file changed, 74 insertions(+), 58 deletions(-) diff --git a/src/bindings.rs b/src/bindings.rs index cbe44da..d0cc6d2 100644 --- a/src/bindings.rs +++ b/src/bindings.rs @@ -432,7 +432,7 @@ pub mod ciarlet { } } - pub fn ciarlet_element_tabulate< + pub unsafe fn ciarlet_element_tabulate< E: FiniteElement, TGeo: RlstScalar, >( @@ -461,14 +461,16 @@ 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_tabulate_f32>( + pub unsafe fn ciarlet_element_tabulate_f32>( element: &E, points: *const f32, npoints: usize, nderivs: usize, data: *mut c_void, ) { - ciarlet_element_tabulate(element, points, npoints, nderivs, data); + unsafe { + ciarlet_element_tabulate(element, points, npoints, nderivs, data); + } } #[concretise_types( @@ -476,14 +478,16 @@ 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_tabulate_f64>( + pub unsafe fn ciarlet_element_tabulate_f64>( element: &E, points: *const f64, npoints: usize, nderivs: usize, data: *mut c_void, ) { - ciarlet_element_tabulate(element, points, npoints, nderivs, data); + unsafe { + ciarlet_element_tabulate(element, points, npoints, nderivs, data); + } } #[concretise_types( @@ -558,7 +562,8 @@ pub mod ciarlet { } } - pub fn ciarlet_element_push_forward< + #[allow(clippy::too_many_arguments)] + pub unsafe fn ciarlet_element_push_forward< E: MappedFiniteElement, TGeo: RlstScalar, >( @@ -620,7 +625,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_f32< + pub unsafe fn ciarlet_element_push_forward_f32< E: MappedFiniteElement, >( element: &E, @@ -634,18 +639,20 @@ pub mod ciarlet { jinv: *const f32, physical_values: *mut c_void, ) { - ciarlet_element_push_forward( - element, - npoints, - nfunctions, - gdim, - reference_values, - nderivs, - j, - jdet, - jinv, - physical_values, - ); + unsafe { + ciarlet_element_push_forward( + element, + npoints, + nfunctions, + gdim, + reference_values, + nderivs, + j, + jdet, + jinv, + physical_values, + ); + } } #[allow(clippy::too_many_arguments)] @@ -654,7 +661,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_f64< + pub unsafe fn ciarlet_element_push_forward_f64< E: MappedFiniteElement, >( element: &E, @@ -668,21 +675,24 @@ pub mod ciarlet { jinv: *const f64, physical_values: *mut c_void, ) { - ciarlet_element_push_forward( - element, - npoints, - nfunctions, - gdim, - reference_values, - nderivs, - j, - jdet, - jinv, - physical_values, - ); + unsafe { + ciarlet_element_push_forward( + element, + npoints, + nfunctions, + gdim, + reference_values, + nderivs, + j, + jdet, + jinv, + physical_values, + ); + } } - pub fn ciarlet_element_pull_back< + #[allow(clippy::too_many_arguments)] + pub unsafe fn ciarlet_element_pull_back< E: MappedFiniteElement, TGeo: RlstScalar, >( @@ -744,7 +754,9 @@ 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_f32>( + pub unsafe fn ciarlet_element_pull_back_f32< + E: MappedFiniteElement, + >( element: &E, npoints: usize, nfunctions: usize, @@ -756,18 +768,20 @@ pub mod ciarlet { jinv: *const f32, reference_values: *mut c_void, ) { - ciarlet_element_pull_back( - element, - npoints, - nfunctions, - gdim, - physical_values, - nderivs, - j, - jdet, - jinv, - reference_values, - ); + unsafe { + ciarlet_element_pull_back( + element, + npoints, + nfunctions, + gdim, + physical_values, + nderivs, + j, + jdet, + jinv, + reference_values, + ); + } } #[allow(clippy::too_many_arguments)] @@ -788,18 +802,20 @@ pub mod ciarlet { jinv: *const f64, reference_values: *mut c_void, ) { - ciarlet_element_pull_back( - element, - npoints, - nfunctions, - gdim, - physical_values, - nderivs, - j, - jdet, - jinv, - reference_values, - ); + unsafe { + ciarlet_element_pull_back( + element, + npoints, + nfunctions, + gdim, + physical_values, + nderivs, + j, + jdet, + jinv, + reference_values, + ); + } } #[concretise_types( From e5a9abac67ee586579a87b0b762ad04f9e96962c Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Wed, 26 Nov 2025 16:43:22 +0000 Subject: [PATCH 8/9] fix dtypes --- python/ndelement/ciarlet.py | 91 ++++++++++------------ python/test/test_element.py | 6 ++ src/bindings.rs | 151 +++++++++++++++++++++++------------- 3 files changed, 141 insertions(+), 107 deletions(-) diff --git a/python/ndelement/ciarlet.py b/python/ndelement/ciarlet.py index e1f585d..72c709f 100644 --- a/python/ndelement/ciarlet.py +++ b/python/ndelement/ciarlet.py @@ -163,7 +163,7 @@ def tabulate(self, points: npt.NDArray[np.floating], nderivs: int) -> npt.NDArra if points.dtype == np.float64: _lib.ciarlet_element_tabulate_f64( self._rs_element, - _ffi.cast("void*", points.ctypes.data), + _ffi.cast("double*", points.ctypes.data), points.shape[0], nderivs, _ffi.cast("void*", data.ctypes.data), @@ -171,7 +171,7 @@ def tabulate(self, points: npt.NDArray[np.floating], nderivs: int) -> npt.NDArra elif points.dtype == np.float32: _lib.ciarlet_element_tabulate_f32( self._rs_element, - points.ctypes.data, + _ffi.cast("float*", points.ctypes.data), points.shape[0], nderivs, _ffi.cast("void*", data.ctypes.data), @@ -202,7 +202,7 @@ def push_forward( inverse_jacobians: npt.NDArray[np.floating], ) -> npt.NDArray[np.floating]: """Push values forward to a physical cell.""" - if reference_values.dtype != self.dtype(0): + if reference_values.dtype != self.geo_dtype(0): raise TypeError("reference_values has incorrect type") if jacobians.dtype != self.dtype(0).real.dtype: raise TypeError("jacobians has incorrect type") @@ -220,34 +220,27 @@ def push_forward( data = np.empty(shape, dtype=self.dtype) if reference_values.dtype == np.float64: - _lib.ciarlet_element_push_forward_f64( - self._rs_element, - npts, - nfuncs, - gdim, - reference_values.ctypes.data, - nderivs, - jacobians.ctypes.data, - jacobian_determinants.ctypes.data, - inverse_jacobians.ctypes.data, - _ffi.cast("void*", data.ctypes.data), - ) + push_function = _lib.ciarlet_element_push_forward_f64 + geo_type = "double*" elif reference_values.dtype == np.float32: - _lib.ciarlet_element_push_forward_f32( - self._rs_element, - npts, - nfuncs, - gdim, - reference_values.ctypes.data, - nderivs, - jacobians.ctypes.data, - jacobian_determinants.ctypes.data, - inverse_jacobians.ctypes.data, - _ffi.cast("void*", data.ctypes.data), - ) + push_function = _lib.ciarlet_element_push_forward_f32 + geo_type = "float*" else: raise TypeError(f"Unsupported dtype: {reference_values.dtype}") + push_function( + self._rs_element, + npts, + nfuncs, + gdim, + _ffi.cast(geo_type, reference_values.ctypes.data), + nderivs, + _ffi.cast(geo_type, jacobians.ctypes.data), + _ffi.cast(geo_type, jacobian_determinants.ctypes.data), + _ffi.cast(geo_type, inverse_jacobians.ctypes.data), + _ffi.cast("void*", data.ctypes.data), + ) + return data def pull_back( @@ -259,7 +252,7 @@ def pull_back( inverse_jacobians: npt.NDArray[np.floating], ) -> npt.NDArray[np.floating]: """Push values back from a physical cell.""" - if physical_values.dtype != self.dtype(0): + if physical_values.dtype != self.dtype(0).real.dtype: raise TypeError("physical_values has incorrect type") if jacobians.dtype != self.dtype(0).real.dtype: raise TypeError("jacobians has incorrect type") @@ -277,33 +270,27 @@ def pull_back( data = np.empty(shape, dtype=self.dtype) if physical_values.dtype == np.float64: - _lib.ciarlet_element_push_forward_f64( - self._rs_element, - npts, - nfuncs, - gdim, - physical_values.ctypes.data, - nderivs, - jacobians.ctypes.data, - jacobian_determinants.ctypes.data, - inverse_jacobians.ctypes.data, - _ffi.cast("void*", data.ctypes.data), - ) + pull_function = _lib.ciarlet_element_pull_back_f64 + geo_type = "double*" elif physical_values.dtype == np.float32: - _lib.ciarlet_element_push_forward_f32( - self._rs_element, - npts, - nfuncs, - gdim, - physical_values.ctypes.data, - nderivs, - jacobians.ctypes.data, - jacobian_determinants.ctypes.data, - inverse_jacobians.ctypes.data, - _ffi.cast("void*", data.ctypes.data), - ) + pull_function = _lib.ciarlet_element_pull_back_f32 + geo_type = "float*" else: raise TypeError(f"Unsupported dtype: {physical_values.dtype}") + + pull_function( + self._rs_element, + npts, + nfuncs, + gdim, + _ffi.cast(geo_type, physical_values.ctypes.data), + nderivs, + _ffi.cast(geo_type, jacobians.ctypes.data), + _ffi.cast(geo_type, jacobian_determinants.ctypes.data), + _ffi.cast(geo_type, inverse_jacobians.ctypes.data), + _ffi.cast("void*", data.ctypes.data), + ) + return data def apply_dof_permutations(self, data: npt.NDArray, orientation: np.int32): diff --git a/python/test/test_element.py b/python/test/test_element.py index aac526f..745bf76 100644 --- a/python/test/test_element.py +++ b/python/test/test_element.py @@ -38,12 +38,18 @@ def test_value_size(cell, degree): + [(dt0, dt1) for dt0 in dtypes for dt1 in [np.complex64, np.complex128]], ) def test_incompatible_types(dt0, dt1): + print(1) family = create_family(Family.Lagrange, 2, dtype=dt0) + print(2) element = family.element(ReferenceCellType.Triangle) + print(3) points = np.array([[0.0, 0.0], [0.2, 0.1], [0.8, 0.05]], dtype=dt1) + print(4) with pytest.raises(TypeError): + print(5) element.tabulate(points, 1) + print(6) @pytest.mark.parametrize("dtype", dtypes) diff --git a/src/bindings.rs b/src/bindings.rs index d0cc6d2..0018c41 100644 --- a/src/bindings.rs +++ b/src/bindings.rs @@ -279,8 +279,9 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) )] pub fn element_value_size(element: &E) -> usize { element.value_size() @@ -372,7 +373,12 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - field(arg = 0, name = "element_family", wrapper = "ElementFamilyT", replace_with = ["crate::ciarlet::LagrangeElementFamily<{{dtype}}>", "ciarlet::RaviartThomasElementFamily<{{dtype}}>", "ciarlet::NedelecFirstKindElementFamily<{{dtype}}>"]) + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + field(arg = 0, name = "element_family", wrapper = "ElementFamilyT", replace_with = [ + "ciarlet::LagrangeElementFamily<{{dtype}}, {{dtype2}}>", + "ciarlet::RaviartThomasElementFamily<{{dtype}}, {{dtype2}}>", + "ciarlet::NedelecFirstKindElementFamily<{{dtype}}, {{dtype2}}>" + ]) )] pub fn element_family_create_element>( family: &F, @@ -388,7 +394,12 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - field(arg = 0, name = "element_family", wrapper = "ElementFamilyT", replace_with = ["crate::ciarlet::LagrangeElementFamily<{{dtype}}>", "ciarlet::RaviartThomasElementFamily<{{dtype}}>", "ciarlet::NedelecFirstKindElementFamily<{{dtype}}>"]) + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + field(arg = 0, name = "element_family", wrapper = "ElementFamilyT", replace_with = [ + "ciarlet::LagrangeElementFamily<{{dtype}}, {{dtype2}}>", + "ciarlet::RaviartThomasElementFamily<{{dtype}}, {{dtype2}}>", + "ciarlet::NedelecFirstKindElementFamily<{{dtype}}, {{dtype2}}>" + ]) )] pub fn element_family_dtype< T: RlstScalar + DTypeIdentifier, @@ -401,8 +412,9 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) )] pub fn ciarlet_element_dtype>( _elem: &E, @@ -412,8 +424,9 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) )] pub fn ciarlet_element_tabulate_array_shape( element: &E, @@ -458,8 +471,9 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]), + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]), )] pub unsafe fn ciarlet_element_tabulate_f32>( element: &E, @@ -475,8 +489,9 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]), + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]), )] pub unsafe fn ciarlet_element_tabulate_f64>( element: &E, @@ -492,8 +507,9 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) )] pub fn ciarlet_element_value_size(element: &E) -> usize { element.value_size() @@ -501,8 +517,9 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) )] pub fn ciarlet_element_value_rank(element: &E) -> usize { element.value_shape().len() @@ -510,8 +527,9 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) )] pub fn ciarlet_element_value_shape(element: &E, shape: *mut usize) { for (i, j) in element.value_shape().iter().enumerate() { @@ -523,8 +541,9 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) )] pub fn ciarlet_element_physical_value_size( element: &E, @@ -535,8 +554,9 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) )] pub fn ciarlet_element_physical_value_rank( element: &E, @@ -547,8 +567,9 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) )] pub fn ciarlet_element_physical_value_shape( element: &E, @@ -622,8 +643,9 @@ pub mod ciarlet { #[allow(clippy::too_many_arguments)] #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) )] pub unsafe fn ciarlet_element_push_forward_f32< E: MappedFiniteElement, @@ -658,8 +680,9 @@ pub mod ciarlet { #[allow(clippy::too_many_arguments)] #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) )] pub unsafe fn ciarlet_element_push_forward_f64< E: MappedFiniteElement, @@ -751,8 +774,9 @@ pub mod ciarlet { #[allow(clippy::too_many_arguments)] #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) )] pub unsafe fn ciarlet_element_pull_back_f32< E: MappedFiniteElement, @@ -787,8 +811,9 @@ pub mod ciarlet { #[allow(clippy::too_many_arguments)] #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) )] pub fn ciarlet_element_pull_back_f64>( element: &E, @@ -820,19 +845,21 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) )] - pub fn ciarlet_element_degree( - element: &CiarletElement, + pub fn ciarlet_element_degree( + element: &CiarletElement, ) -> usize { element.degree() } #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) )] pub fn ciarlet_element_embedded_superdegree(element: &E) -> usize { element.lagrange_superdegree() @@ -840,8 +867,9 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) )] pub fn ciarlet_element_dim(element: &E) -> usize { element.dim() @@ -849,19 +877,21 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) )] - pub fn ciarlet_element_continuity( - element: &CiarletElement, + pub fn ciarlet_element_continuity( + element: &CiarletElement, ) -> Continuity { element.continuity() } #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) )] pub fn ciarlet_element_cell_type>( element: &E, @@ -871,8 +901,9 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) )] pub fn ciarlet_element_entity_dofs_size( element: &E, @@ -884,8 +915,9 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) )] pub fn ciarlet_element_entity_dofs( element: &E, @@ -907,11 +939,12 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) )] - pub fn ciarlet_element_entity_closure_dofs_size( - element: &CiarletElement, + pub fn ciarlet_element_entity_closure_dofs_size( + element: &CiarletElement, entity_dim: usize, entity_index: usize, ) -> usize { @@ -923,11 +956,12 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) )] - pub fn ciarlet_element_entity_closure_dofs( - element: &CiarletElement, + pub fn ciarlet_element_entity_closure_dofs( + element: &CiarletElement, entity_dim: usize, entity_index: usize, entity_closure_dofs: *mut usize, @@ -946,14 +980,15 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) )] pub fn ciarlet_element_interpolation_npoints< T: RlstScalar + DTypeIdentifier + Getrf + Getri, - M: Map, + M: Map, TGeo: RlstScalar + DTypeIdentifier >( - element: &CiarletElement, + element: &CiarletElement, entity_dim: usize, entity_index: usize, ) -> usize { @@ -962,14 +997,15 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) )] pub fn ciarlet_element_interpolation_ndofs< T: RlstScalar + DTypeIdentifier + Getrf + Getri, - M: Map, + M: Map, TGeo: RlstScalar + DTypeIdentifier >( - element: &CiarletElement, + element: &CiarletElement, entity_dim: usize, entity_index: usize, ) -> usize { @@ -1007,14 +1043,15 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) )] pub fn ciarlet_element_interpolation_weights< T: RlstScalar + DTypeIdentifier + Getrf + Getri, - M: Map, + M: Map, TGeo: RlstScalar + DTypeIdentifier >( - element: &CiarletElement, + element: &CiarletElement, entity_dim: usize, entity_index: usize, weights: *mut c_void, @@ -1034,14 +1071,15 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) )] pub unsafe fn ciarlet_element_apply_dof_permutations_usize< T: RlstScalar + DTypeIdentifier + Getrf + Getri, - M: Map, + M: Map, TGeo: RlstScalar + DTypeIdentifier >( - element: &CiarletElement, + element: &CiarletElement, data: *mut usize, data_size: usize, orientation: i32, @@ -1053,14 +1091,15 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) )] pub unsafe fn ciarlet_element_apply_dof_permutations< T: RlstScalar + DTypeIdentifier + Getrf + Getri, - M: Map, + M: Map, TGeo: RlstScalar + DTypeIdentifier >( - element: &CiarletElement, + element: &CiarletElement, data: *mut c_void, data_size: usize, orientation: i32, @@ -1073,14 +1112,15 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) )] pub unsafe fn ciarlet_element_apply_dof_transformations< T: RlstScalar + DTypeIdentifier + Getrf + Getri, - M: Map, + M: Map, TGeo: RlstScalar + DTypeIdentifier >( - element: &CiarletElement, + element: &CiarletElement, data: *mut c_void, data_size: usize, orientation: i32, @@ -1095,14 +1135,15 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) )] pub unsafe fn ciarlet_element_apply_dof_permutations_and_transformations< T: RlstScalar + DTypeIdentifier + Getrf + Getri, - M: Map, + M: Map, TGeo: RlstScalar + DTypeIdentifier >( - element: &CiarletElement, + element: &CiarletElement, data: *mut c_void, data_size: usize, orientation: i32, From 59f65c1864960df9a99344be9e9ac23483d04492 Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Wed, 26 Nov 2025 16:55:43 +0000 Subject: [PATCH 9/9] geo_dtype --- python/ndelement/ciarlet.py | 27 ++--- python/test/test_element.py | 25 ----- src/bindings.rs | 201 +++++++++++++++++++++--------------- 3 files changed, 134 insertions(+), 119 deletions(-) diff --git a/python/ndelement/ciarlet.py b/python/ndelement/ciarlet.py index 72c709f..8f2306f 100644 --- a/python/ndelement/ciarlet.py +++ b/python/ndelement/ciarlet.py @@ -52,6 +52,11 @@ def dtype(self) -> typing.Type[np.floating]: """Data type.""" return _dtypes[_lib.ciarlet_element_dtype(self._rs_element)] + @property + def geo_dtype(self) -> typing.Type[np.floating]: + """Data type.""" + return _dtypes[_lib.ciarlet_element_geo_dtype(self._rs_element)] + @property def value_size(self) -> int: """Value size of the element.""" @@ -123,7 +128,7 @@ def interpolation_points(self) -> list[list[npt.NDArray]]: points_d = [] for i in range(n): shape = (_lib.ciarlet_element_interpolation_npoints(self._rs_element, d, i), tdim) - points_di = np.empty(shape, dtype=self.dtype(0).real.dtype) + points_di = np.empty(shape, dtype=self.geo_dtype) _lib.ciarlet_element_interpolation_points( self._rs_element, d, i, _ffi.cast("void*", points_di.ctypes.data) ) @@ -152,8 +157,6 @@ def interpolation_weights(self) -> list[list[npt.NDArray]]: def tabulate(self, points: npt.NDArray[np.floating], nderivs: int) -> npt.NDArray: """Tabulate the basis functions at a set of points.""" - if points.dtype != self.dtype(0).real.dtype: - raise TypeError("points has incorrect type") shape = np.empty(4, dtype=np.uintp) _lib.ciarlet_element_tabulate_array_shape( self._rs_element, nderivs, points.shape[0], _ffi.cast("uintptr_t*", shape.ctypes.data) @@ -202,13 +205,12 @@ def push_forward( inverse_jacobians: npt.NDArray[np.floating], ) -> npt.NDArray[np.floating]: """Push values forward to a physical cell.""" - if reference_values.dtype != self.geo_dtype(0): - raise TypeError("reference_values has incorrect type") - if jacobians.dtype != self.dtype(0).real.dtype: + geo_dtype = reference_values.dtype + if jacobians.dtype != geo_dtype: raise TypeError("jacobians has incorrect type") - if jacobian_determinants.dtype != self.dtype(0).real.dtype: + if jacobian_determinants.dtype != geo_dtype: raise TypeError("jacobian_determinants has incorrect type") - if inverse_jacobians.dtype != self.dtype(0).real.dtype: + if inverse_jacobians.dtype != geo_dtype: raise TypeError("inverse_jacobians has incorrect type") gdim = jacobians.shape[1] @@ -252,13 +254,12 @@ def pull_back( inverse_jacobians: npt.NDArray[np.floating], ) -> npt.NDArray[np.floating]: """Push values back from a physical cell.""" - if physical_values.dtype != self.dtype(0).real.dtype: - raise TypeError("physical_values has incorrect type") - if jacobians.dtype != self.dtype(0).real.dtype: + geo_dtype = physical_values.dtype + if jacobians.dtype != geo_dtype: raise TypeError("jacobians has incorrect type") - if jacobian_determinants.dtype != self.dtype(0).real.dtype: + if jacobian_determinants.dtype != geo_dtype: raise TypeError("jacobian_determinants has incorrect type") - if inverse_jacobians.dtype != self.dtype(0).real.dtype: + if inverse_jacobians.dtype != geo_dtype: raise TypeError("inverse_jacobians has incorrect type") gdim = jacobians.shape[1] diff --git a/python/test/test_element.py b/python/test/test_element.py index 745bf76..5f391dd 100644 --- a/python/test/test_element.py +++ b/python/test/test_element.py @@ -27,31 +27,6 @@ def test_value_size(cell, degree): assert element.value_size == 1 -@pytest.mark.parametrize( - "dt0,dt1", - [ - (np.float32, np.float64), - (np.float64, np.float32), - (np.complex64, np.float64), - (np.complex128, np.float32), - ] - + [(dt0, dt1) for dt0 in dtypes for dt1 in [np.complex64, np.complex128]], -) -def test_incompatible_types(dt0, dt1): - print(1) - family = create_family(Family.Lagrange, 2, dtype=dt0) - print(2) - element = family.element(ReferenceCellType.Triangle) - print(3) - points = np.array([[0.0, 0.0], [0.2, 0.1], [0.8, 0.05]], dtype=dt1) - print(4) - - with pytest.raises(TypeError): - print(5) - element.tabulate(points, 1) - print(6) - - @pytest.mark.parametrize("dtype", dtypes) def test_lagrange_2_triangle_tabulate(dtype): family = create_family(Family.Lagrange, 2, dtype=dtype) diff --git a/src/bindings.rs b/src/bindings.rs index 0018c41..68c7c61 100644 --- a/src/bindings.rs +++ b/src/bindings.rs @@ -279,9 +279,9 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]) )] pub fn element_value_size(element: &E) -> usize { element.value_size() @@ -373,11 +373,11 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), field(arg = 0, name = "element_family", wrapper = "ElementFamilyT", replace_with = [ - "ciarlet::LagrangeElementFamily<{{dtype}}, {{dtype2}}>", - "ciarlet::RaviartThomasElementFamily<{{dtype}}, {{dtype2}}>", - "ciarlet::NedelecFirstKindElementFamily<{{dtype}}, {{dtype2}}>" + "ciarlet::LagrangeElementFamily<{{dtype}}, {{geo_dtype}}>", + "ciarlet::RaviartThomasElementFamily<{{dtype}}, {{geo_dtype}}>", + "ciarlet::NedelecFirstKindElementFamily<{{dtype}}, {{geo_dtype}}>" ]) )] pub fn element_family_create_element>( @@ -394,11 +394,11 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), field(arg = 0, name = "element_family", wrapper = "ElementFamilyT", replace_with = [ - "ciarlet::LagrangeElementFamily<{{dtype}}, {{dtype2}}>", - "ciarlet::RaviartThomasElementFamily<{{dtype}}, {{dtype2}}>", - "ciarlet::NedelecFirstKindElementFamily<{{dtype}}, {{dtype2}}>" + "ciarlet::LagrangeElementFamily<{{dtype}}, {{geo_dtype}}>", + "ciarlet::RaviartThomasElementFamily<{{dtype}}, {{geo_dtype}}>", + "ciarlet::NedelecFirstKindElementFamily<{{dtype}}, {{geo_dtype}}>" ]) )] pub fn element_family_dtype< @@ -412,9 +412,9 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]) )] pub fn ciarlet_element_dtype>( _elem: &E, @@ -424,9 +424,25 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]) + )] + pub fn ciarlet_element_geo_dtype< + T: RlstScalar + DTypeIdentifier, + M: Map, + TGeo: RlstScalar + DTypeIdentifier, + >( + _elem: &CiarletElement, + ) -> DType { + ::dtype() + } + + #[concretise_types( + gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]) )] pub fn ciarlet_element_tabulate_array_shape( element: &E, @@ -471,9 +487,9 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]), + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]), )] pub unsafe fn ciarlet_element_tabulate_f32>( element: &E, @@ -489,9 +505,9 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]), + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]), )] pub unsafe fn ciarlet_element_tabulate_f64>( element: &E, @@ -507,9 +523,9 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]) )] pub fn ciarlet_element_value_size(element: &E) -> usize { element.value_size() @@ -517,9 +533,9 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]) )] pub fn ciarlet_element_value_rank(element: &E) -> usize { element.value_shape().len() @@ -527,9 +543,9 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]) )] pub fn ciarlet_element_value_shape(element: &E, shape: *mut usize) { for (i, j) in element.value_shape().iter().enumerate() { @@ -541,9 +557,9 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]) )] pub fn ciarlet_element_physical_value_size( element: &E, @@ -554,9 +570,9 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]) )] pub fn ciarlet_element_physical_value_rank( element: &E, @@ -567,9 +583,9 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]) )] pub fn ciarlet_element_physical_value_shape( element: &E, @@ -643,9 +659,9 @@ pub mod ciarlet { #[allow(clippy::too_many_arguments)] #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]) )] pub unsafe fn ciarlet_element_push_forward_f32< E: MappedFiniteElement, @@ -680,9 +696,9 @@ pub mod ciarlet { #[allow(clippy::too_many_arguments)] #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]) )] pub unsafe fn ciarlet_element_push_forward_f64< E: MappedFiniteElement, @@ -774,9 +790,9 @@ pub mod ciarlet { #[allow(clippy::too_many_arguments)] #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]) )] pub unsafe fn ciarlet_element_pull_back_f32< E: MappedFiniteElement, @@ -811,9 +827,9 @@ pub mod ciarlet { #[allow(clippy::too_many_arguments)] #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]) )] pub fn ciarlet_element_pull_back_f64>( element: &E, @@ -845,11 +861,15 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]) )] - pub fn ciarlet_element_degree( + pub fn ciarlet_element_degree< + T: RlstScalar + DTypeIdentifier + Getrf + Getri, + M: Map, + TGeo: RlstScalar + DTypeIdentifier, + >( element: &CiarletElement, ) -> usize { element.degree() @@ -857,9 +877,9 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]) )] pub fn ciarlet_element_embedded_superdegree(element: &E) -> usize { element.lagrange_superdegree() @@ -867,9 +887,9 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]) )] pub fn ciarlet_element_dim(element: &E) -> usize { element.dim() @@ -877,11 +897,15 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]) )] - pub fn ciarlet_element_continuity( + pub fn ciarlet_element_continuity< + T: RlstScalar + DTypeIdentifier + Getrf + Getri, + M: Map, + TGeo: RlstScalar + DTypeIdentifier, + >( element: &CiarletElement, ) -> Continuity { element.continuity() @@ -889,9 +913,9 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]) )] pub fn ciarlet_element_cell_type>( element: &E, @@ -901,9 +925,9 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]) )] pub fn ciarlet_element_entity_dofs_size( element: &E, @@ -915,9 +939,9 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]) )] pub fn ciarlet_element_entity_dofs( element: &E, @@ -939,11 +963,15 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]) )] - pub fn ciarlet_element_entity_closure_dofs_size( + pub fn ciarlet_element_entity_closure_dofs_size< + T: RlstScalar + DTypeIdentifier, + M: Map, + TGeo: RlstScalar + DTypeIdentifier, + >( element: &CiarletElement, entity_dim: usize, entity_index: usize, @@ -956,11 +984,15 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]) )] - pub fn ciarlet_element_entity_closure_dofs( + pub fn ciarlet_element_entity_closure_dofs< + T: RlstScalar + DTypeIdentifier, + M: Map, + TGeo: RlstScalar + DTypeIdentifier, + >( element: &CiarletElement, entity_dim: usize, entity_index: usize, @@ -980,13 +1012,14 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]) )] pub fn ciarlet_element_interpolation_npoints< T: RlstScalar + DTypeIdentifier + Getrf + Getri, - M: Map, TGeo: RlstScalar + DTypeIdentifier + M: Map, + TGeo: RlstScalar + DTypeIdentifier, >( element: &CiarletElement, entity_dim: usize, @@ -997,13 +1030,14 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]) )] pub fn ciarlet_element_interpolation_ndofs< T: RlstScalar + DTypeIdentifier + Getrf + Getri, - M: Map, TGeo: RlstScalar + DTypeIdentifier + M: Map, + TGeo: RlstScalar + DTypeIdentifier, >( element: &CiarletElement, entity_dim: usize, @@ -1043,13 +1077,14 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]) )] pub fn ciarlet_element_interpolation_weights< T: RlstScalar + DTypeIdentifier + Getrf + Getri, - M: Map, TGeo: RlstScalar + DTypeIdentifier + M: Map, + TGeo: RlstScalar + DTypeIdentifier, >( element: &CiarletElement, entity_dim: usize, @@ -1071,13 +1106,14 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]) )] pub unsafe fn ciarlet_element_apply_dof_permutations_usize< T: RlstScalar + DTypeIdentifier + Getrf + Getri, - M: Map, TGeo: RlstScalar + DTypeIdentifier + M: Map, + TGeo: RlstScalar + DTypeIdentifier, >( element: &CiarletElement, data: *mut usize, @@ -1091,13 +1127,14 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]) )] pub unsafe fn ciarlet_element_apply_dof_permutations< T: RlstScalar + DTypeIdentifier + Getrf + Getri, - M: Map, TGeo: RlstScalar + DTypeIdentifier + M: Map, + TGeo: RlstScalar + DTypeIdentifier, >( element: &CiarletElement, data: *mut c_void, @@ -1112,13 +1149,14 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]) )] pub unsafe fn ciarlet_element_apply_dof_transformations< T: RlstScalar + DTypeIdentifier + Getrf + Getri, - M: Map, TGeo: RlstScalar + DTypeIdentifier + M: Map, + TGeo: RlstScalar + DTypeIdentifier, >( element: &CiarletElement, data: *mut c_void, @@ -1135,13 +1173,14 @@ pub mod ciarlet { #[concretise_types( gen_type(name = "dtype", replace_with = ["f32", "f64", "c32", "c64"]), - gen_type(name = "dtype2", replace_with = ["f32", "f64", "c32", "c64"]), + gen_type(name = "geo_dtype", replace_with = ["f32", "f64", "c32", "c64"]), gen_type(name = "maptype", replace_with = ["IdentityMap", "CovariantPiolaMap", "ContravariantPiolaMap"]), - field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{dtype2}}>"]) + field(arg = 0, name = "element", wrapper = "CiarletElementT", replace_with = ["CiarletElement<{{dtype}}, {{maptype}}, {{geo_dtype}}>"]) )] pub unsafe fn ciarlet_element_apply_dof_permutations_and_transformations< T: RlstScalar + DTypeIdentifier + Getrf + Getri, - M: Map, TGeo: RlstScalar + DTypeIdentifier + M: Map, + TGeo: RlstScalar + DTypeIdentifier, >( element: &CiarletElement, data: *mut c_void,