From b249254406c91dcf2cb4c503a2e6766e94ad9077 Mon Sep 17 00:00:00 2001 From: Timo Betcke Date: Thu, 6 Nov 2025 17:55:41 +0000 Subject: [PATCH 1/4] Cleaned up RLST traits --- Cargo.toml | 2 + src/ciarlet.rs | 49 ++++++------ src/map.rs | 154 ++++++++++++++---------------------- src/math.rs | 15 ++-- src/polynomials/legendre.rs | 110 +++++++++++++------------- src/traits.rs | 69 ++++++++-------- 6 files changed, 183 insertions(+), 216 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 8b2c82a..dac1328 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -33,10 +33,12 @@ strum_macros = "0.27" c-api-tools = { version = "0.1.0" } [target.'cfg(target_os = "macos")'.dependencies] +blas-src = { version = "0.14", features = ["accelerate"]} lapack-src = { version = "0.13", features = ["accelerate"]} [target.'cfg(target_os = "linux")'.dependencies] lapack-src = { version = "0.13", features = ["netlib"]} +blas-src = { version = "0.14", features = ["blis"]} [dev-dependencies] paste = "1.*" diff --git a/src/ciarlet.rs b/src/ciarlet.rs index cfeb183..e95c945 100644 --- a/src/ciarlet.rs +++ b/src/ciarlet.rs @@ -9,12 +9,8 @@ use crate::traits::{FiniteElement, Map}; use crate::types::{Continuity, DofTransformation, ReferenceCellType, Transformation}; use itertools::izip; use num::{One, Zero}; -use rlst::{ - Array, DynArray, Inverse, RandomAccessByRef, RandomAccessMut, RlstScalar, Shape, - UnsafeRandomAccessByRef, UnsafeRandomAccessMut, - dense::linalg::lapack::interface::{getrf::Getrf, getri::Getri}, - rlst_dynamic_array, -}; +use rlst::RlstScalar; +use rlst::{Array, DynArray, Inverse, MutableArrayImpl, ValueArrayImpl, rlst_dynamic_array}; use std::collections::HashMap; use std::fmt::{Debug, Formatter}; @@ -71,7 +67,10 @@ impl Debug for CiarletElement { } } -impl CiarletElement { +impl CiarletElement +where + DynArray: Inverse>, +{ /// Create a Ciarlet element #[allow(clippy::too_many_arguments)] pub fn create( @@ -610,8 +609,8 @@ impl FiniteElement for CiarletElement { self.dim } fn tabulate< - Array2: RandomAccessByRef<2, Item = T::Real> + Shape<2>, - Array4Mut: UnsafeRandomAccessMut<4, Item = T> + RandomAccessMut<4, Item = T>, + Array2: ValueArrayImpl<::Real, 2>, + Array4Mut: MutableArrayImpl, >( &self, points: &Array, @@ -671,17 +670,17 @@ impl FiniteElement for CiarletElement { [deriv_count, point_count, basis_count, value_size] } fn push_forward< - Array3Real: UnsafeRandomAccessByRef<3, Item = ::Real> + Shape<3>, - Array4: UnsafeRandomAccessByRef<4, Item = Self::T> + Shape<4>, - Array4Mut: UnsafeRandomAccessMut<4, Item = Self::T> + Shape<4>, + Array3RealImpl: ValueArrayImpl<::Real, 3>, + Array4Impl: ValueArrayImpl, + Array4MutImpl: MutableArrayImpl, >( &self, - reference_values: &Array, + reference_values: &Array, nderivs: usize, - jacobians: &Array, + jacobians: &Array, jacobian_determinants: &[::Real], - inverse_jacobians: &Array, - physical_values: &mut Array, + inverse_jacobians: &Array, + physical_values: &mut Array, ) { self.map.push_forward( reference_values, @@ -693,17 +692,17 @@ impl FiniteElement for CiarletElement { ) } fn pull_back< - Array3Real: UnsafeRandomAccessByRef<3, Item = ::Real> + Shape<3>, - Array4: UnsafeRandomAccessByRef<4, Item = Self::T> + Shape<4>, - Array4Mut: UnsafeRandomAccessMut<4, Item = Self::T> + Shape<4>, + Array3RealImpl: ValueArrayImpl<::Real, 3>, + Array4Impl: ValueArrayImpl, + Array4MutImpl: MutableArrayImpl, >( &self, - physical_values: &Array, + physical_values: &Array, nderivs: usize, - jacobians: &Array, + jacobians: &Array, jacobian_determinants: &[::Real], - inverse_jacobians: &Array, - reference_values: &mut Array, + inverse_jacobians: &Array, + reference_values: &mut Array, ) { self.map.pull_back( physical_values, @@ -1875,9 +1874,9 @@ mod test { test_dof_transformations!(Hexahedron, raviart_thomas, 1); test_dof_transformations!(Hexahedron, raviart_thomas, 2); - fn assert_dof_transformation_equal + Shape<2>>( + fn assert_dof_transformation_equal>( p: &[usize], - m: &Array, + m: &Array, expected: Vec>, ) { let ndofs = p.len(); diff --git a/src/map.rs b/src/map.rs index d1ffcc3..648e21d 100644 --- a/src/map.rs +++ b/src/map.rs @@ -1,6 +1,6 @@ //! Maps from the reference cell to/from physical cells use crate::traits::Map; -use rlst::{Array, RlstScalar, Shape, UnsafeRandomAccessByRef, UnsafeRandomAccessMut}; +use rlst::{Array, MutableArrayImpl, RlstScalar, ValueArrayImpl}; /// Identity map pub struct IdentityMap {} @@ -8,67 +8,43 @@ pub struct IdentityMap {} impl Map for IdentityMap { fn push_forward< T: RlstScalar, - Array3Real: UnsafeRandomAccessByRef<3, Item = T::Real> + Shape<3>, - Array4: UnsafeRandomAccessByRef<4, Item = T> + Shape<4>, - Array4Mut: UnsafeRandomAccessMut<4, Item = T> + Shape<4>, + Array3RealImpl: ValueArrayImpl, + Array4Impl: ValueArrayImpl, + Array4MutImpl: MutableArrayImpl, >( &self, - reference_values: &Array, + reference_values: &Array, nderivs: usize, - _jacobians: &Array, + _jacobians: &Array, _jacobian_determinants: &[T::Real], - _inverse_jacobians: &Array, - physical_values: &mut Array, + _inverse_jacobians: &Array, + physical_values: &mut Array, ) { - assert_eq!(reference_values.shape()[0], physical_values.shape()[0]); - assert_eq!(reference_values.shape()[1], physical_values.shape()[1]); - assert_eq!(reference_values.shape()[2], physical_values.shape()[2]); - assert_eq!(reference_values.shape()[3], physical_values.shape()[3]); if nderivs > 0 { unimplemented!(); } - for i0 in 0..reference_values.shape()[0] { - for i1 in 0..reference_values.shape()[1] { - for i2 in 0..reference_values.shape()[2] { - for i3 in 0..reference_values.shape()[3] { - *physical_values.get_mut([i0, i1, i2, i3]).unwrap() = - *reference_values.get([i0, i1, i2, i3]).unwrap(); - } - } - } - } + + physical_values.fill_from(reference_values); } fn pull_back< T: RlstScalar, - Array3Real: UnsafeRandomAccessByRef<3, Item = T::Real> + Shape<3>, - Array4: UnsafeRandomAccessByRef<4, Item = T> + Shape<4>, - Array4Mut: UnsafeRandomAccessMut<4, Item = T> + Shape<4>, + Array3RealImpl: ValueArrayImpl, + Array4Impl: ValueArrayImpl, + Array4MutImpl: MutableArrayImpl, >( &self, - physical_values: &Array, + physical_values: &Array, nderivs: usize, - _jacobians: &Array, + _jacobians: &Array, _jacobian_determinants: &[T::Real], - _inverse_jacobians: &Array, - reference_values: &mut Array, + _inverse_jacobians: &Array, + reference_values: &mut Array, ) { - assert_eq!(reference_values.shape()[0], physical_values.shape()[0]); - assert_eq!(reference_values.shape()[1], physical_values.shape()[1]); - assert_eq!(reference_values.shape()[2], physical_values.shape()[2]); - assert_eq!(reference_values.shape()[3], physical_values.shape()[3]); if nderivs > 0 { unimplemented!(); } - for i0 in 0..physical_values.shape()[0] { - for i1 in 0..physical_values.shape()[1] { - for i2 in 0..physical_values.shape()[2] { - for i3 in 0..physical_values.shape()[3] { - *reference_values.get_mut([i0, i1, i2, i3]).unwrap() = - *physical_values.get([i0, i1, i2, i3]).unwrap(); - } - } - } - } + + reference_values.fill_from(physical_values); } fn physical_value_shape(&self, _gdim: usize) -> Vec { vec![1] @@ -81,17 +57,17 @@ pub struct CovariantPiolaMap {} impl Map for CovariantPiolaMap { fn push_forward< T: RlstScalar, - Array3Real: UnsafeRandomAccessByRef<3, Item = T::Real> + Shape<3>, - Array4: UnsafeRandomAccessByRef<4, Item = T> + Shape<4>, - Array4Mut: UnsafeRandomAccessMut<4, Item = T> + Shape<4>, + Array3RealImpl: ValueArrayImpl, + Array4Impl: ValueArrayImpl, + Array4MutImpl: MutableArrayImpl, >( &self, - reference_values: &Array, + reference_values: &Array, nderivs: usize, - _jacobians: &Array, + _jacobians: &Array, _jacobian_determinants: &[T::Real], - inverse_jacobians: &Array, - physical_values: &mut Array, + inverse_jacobians: &Array, + physical_values: &mut Array, ) { let tdim = inverse_jacobians.shape()[1]; let gdim = inverse_jacobians.shape()[2]; @@ -109,8 +85,8 @@ impl Map for CovariantPiolaMap { unsafe { *physical_values.get_unchecked_mut([0, p, b, gd]) = (0..tdim) .map(|td| { - T::from(*inverse_jacobians.get_unchecked([p, td, gd])).unwrap() - * *reference_values.get_unchecked([0, p, b, td]) + T::from(inverse_jacobians.get_value_unchecked([p, td, gd])).unwrap() + * reference_values.get_value_unchecked([0, p, b, td]) }) .sum::(); } @@ -120,17 +96,17 @@ impl Map for CovariantPiolaMap { } fn pull_back< T: RlstScalar, - Array3Real: UnsafeRandomAccessByRef<3, Item = T::Real> + Shape<3>, - Array4: UnsafeRandomAccessByRef<4, Item = T> + Shape<4>, - Array4Mut: UnsafeRandomAccessMut<4, Item = T> + Shape<4>, + Array3RealImpl: ValueArrayImpl, + Array4Impl: ValueArrayImpl, + Array4MutImpl: MutableArrayImpl, >( &self, - physical_values: &Array, + physical_values: &Array, nderivs: usize, - jacobians: &Array, + jacobians: &Array, _jacobian_determinants: &[T::Real], - _inverse_jacobians: &Array, - reference_values: &mut Array, + _inverse_jacobians: &Array, + reference_values: &mut Array, ) { let gdim = jacobians.shape()[1]; let tdim = jacobians.shape()[2]; @@ -148,8 +124,8 @@ impl Map for CovariantPiolaMap { unsafe { *reference_values.get_unchecked_mut([0, p, b, td]) = (0..gdim) .map(|gd| { - T::from(*jacobians.get_unchecked([p, gd, td])).unwrap() - * *physical_values.get_unchecked([0, p, b, gd]) + T::from(jacobians.get_value_unchecked([p, gd, td])).unwrap() + * physical_values.get_value_unchecked([0, p, b, gd]) }) .sum::(); } @@ -168,17 +144,17 @@ pub struct ContravariantPiolaMap {} impl Map for ContravariantPiolaMap { fn push_forward< T: RlstScalar, - Array3Real: UnsafeRandomAccessByRef<3, Item = T::Real> + Shape<3>, - Array4: UnsafeRandomAccessByRef<4, Item = T> + Shape<4>, - Array4Mut: UnsafeRandomAccessMut<4, Item = T> + Shape<4>, + Array3RealImpl: ValueArrayImpl, + Array4Impl: ValueArrayImpl, + Array4MutImpl: MutableArrayImpl, >( &self, - reference_values: &Array, + reference_values: &Array, nderivs: usize, - jacobians: &Array, + jacobians: &Array, jacobian_determinants: &[T::Real], - _inverse_jacobians: &Array, - physical_values: &mut Array, + _inverse_jacobians: &Array, + physical_values: &mut Array, ) { let gdim = jacobians.shape()[1]; let tdim = jacobians.shape()[2]; @@ -196,8 +172,8 @@ impl Map for ContravariantPiolaMap { unsafe { *physical_values.get_unchecked_mut([0, p, b, gd]) = (0..tdim) .map(|td| { - T::from(*jacobians.get_unchecked([p, gd, td])).unwrap() - * *reference_values.get_unchecked([0, p, b, td]) + T::from(jacobians.get_value_unchecked([p, gd, td])).unwrap() + * reference_values.get_value_unchecked([0, p, b, td]) }) .sum::() / T::from(*jacobian_determinants.get_unchecked(p)).unwrap(); @@ -208,17 +184,17 @@ impl Map for ContravariantPiolaMap { } fn pull_back< T: RlstScalar, - Array3Real: UnsafeRandomAccessByRef<3, Item = T::Real> + Shape<3>, - Array4: UnsafeRandomAccessByRef<4, Item = T> + Shape<4>, - Array4Mut: UnsafeRandomAccessMut<4, Item = T> + Shape<4>, + Array3RealImpl: ValueArrayImpl, + Array4Impl: ValueArrayImpl, + Array4MutImpl: MutableArrayImpl, >( &self, - physical_values: &Array, + physical_values: &Array, nderivs: usize, - _jacobians: &Array, + _jacobians: &Array, jacobian_determinants: &[T::Real], - inverse_jacobians: &Array, - reference_values: &mut Array, + inverse_jacobians: &Array, + reference_values: &mut Array, ) { let tdim = inverse_jacobians.shape()[1]; let gdim = inverse_jacobians.shape()[2]; @@ -236,8 +212,8 @@ impl Map for ContravariantPiolaMap { unsafe { *reference_values.get_unchecked_mut([0, p, b, td]) = (0..gdim) .map(|gd| { - T::from(*inverse_jacobians.get_unchecked([p, td, gd])).unwrap() - * *physical_values.get_unchecked([0, p, b, gd]) + T::from(inverse_jacobians.get_value_unchecked([p, td, gd])).unwrap() + * physical_values.get_value_unchecked([0, p, b, gd]) }) .sum::() * T::from(*jacobian_determinants.get_unchecked(p)).unwrap(); @@ -257,19 +233,6 @@ mod test { use approx::*; use rlst::{Array, RandomAccessMut, rlst_dynamic_array}; - fn set_to_zero + Shape<4>>( - data: &mut Array, - ) { - for i0 in 0..data.shape()[0] { - for i1 in 0..data.shape()[1] { - for i2 in 0..data.shape()[2] { - for i3 in 0..data.shape()[3] { - *data.get_mut([i0, i1, i2, i3]).unwrap() = T::from(0.0).unwrap(); - } - } - } - } - } fn fill_jacobians< T: RlstScalar, Array3: RandomAccessMut<3, Item = T>, @@ -338,7 +301,8 @@ mod test { epsilon = 1e-14 ); - set_to_zero(&mut values); + values.set_zero(); + map.pull_back(&physical_values, 0, &j, &jdet, &jinv, &mut values); assert_relative_eq!(*values.get([0, 0, 0, 0]).unwrap(), 1.0, epsilon = 1e-14); @@ -410,7 +374,7 @@ mod test { epsilon = 1e-14 ); - set_to_zero(&mut values); + values.set_zero(); map.pull_back(&physical_values, 0, &j, &jdet, &jinv, &mut values); assert_relative_eq!(*values.get([0, 0, 0, 0]).unwrap(), 1.0, epsilon = 1e-14); @@ -486,7 +450,7 @@ mod test { epsilon = 1e-14 ); - set_to_zero(&mut values); + values.set_zero(); map.pull_back(&physical_values, 0, &j, &jdet, &jinv, &mut values); assert_relative_eq!(*values.get([0, 0, 0, 0]).unwrap(), 1.0, epsilon = 1e-14); diff --git a/src/math.rs b/src/math.rs index 39c8ca2..344d34e 100644 --- a/src/math.rs +++ b/src/math.rs @@ -1,7 +1,7 @@ //! Mathematical functions use rlst::{ Array, RandomAccessByRef, RandomAccessMut, RlstScalar, Shape, UnsafeRandomAccessByRef, - UnsafeRandomAccessMut, + UnsafeRandomAccessMut, ValueArrayImpl, }; /// Orthogonalise the rows of a matrix, starting with the row numbered `start` @@ -167,8 +167,8 @@ pub fn prepare_matrix< } /// Apply a permutation and a matrix to some data -pub fn apply_perm_and_matrix + Shape<2>>( - mat: &Array, +pub fn apply_perm_and_matrix>( + mat: &Array, perm: &[usize], data: &mut [T], ) { @@ -177,7 +177,7 @@ pub fn apply_perm_and_matrix + Shape<2>>( +pub fn apply_matrix>( mat: &Array, data: &mut [T], ) { @@ -187,18 +187,19 @@ pub fn apply_matrix + Shap for i in 0..dim { for j in i + 1..dim { for k in 0..block_size { - data[i * block_size + k] += *mat.get([i, j]).unwrap() * data[j * block_size + k]; + data[i * block_size + k] += + mat.get_value([i, j]).unwrap() * data[j * block_size + k]; } } } for i in 1..=dim { for k in 0..block_size { - data[(dim - i) * block_size + k] *= *mat.get([dim - i, dim - i]).unwrap(); + data[(dim - i) * block_size + k] *= mat.get_value([dim - i, dim - i]).unwrap(); } for j in 0..dim - i { for k in 0..block_size { data[(dim - i) * block_size + k] += - *mat.get([dim - i, j]).unwrap() * data[j * block_size + k]; + mat.get_value([dim - i, j]).unwrap() * data[j * block_size + k]; } } } diff --git a/src/polynomials/legendre.rs b/src/polynomials/legendre.rs index f5731cd..f6267af 100644 --- a/src/polynomials/legendre.rs +++ b/src/polynomials/legendre.rs @@ -4,7 +4,7 @@ //! use super::{derivative_count, polynomial_count}; use crate::types::ReferenceCellType; -use rlst::{Array, RandomAccessByRef, RandomAccessMut, RlstScalar, Shape}; +use rlst::{Array, MutableArrayImpl, RlstScalar, ValueArrayImpl}; fn tri_index(i: usize, j: usize) -> usize { (i + j + 1) * (i + j) / 2 + j @@ -37,13 +37,13 @@ fn jrc(a: usize, n: usize) -> (T, T, T) { /// Tabulate orthonormal polynomials on an interval fn tabulate_interval< T: RlstScalar, - Array2: RandomAccessByRef<2, Item = T::Real> + Shape<2>, - Array3Mut: RandomAccessMut<3, Item = T> + RandomAccessByRef<3, Item = T> + Shape<3>, + Array2Impl: ValueArrayImpl, + Array3MutImpl: MutableArrayImpl, >( - points: &Array, + points: &Array, degree: usize, derivatives: usize, - data: &mut Array, + data: &mut Array, ) { debug_assert!(data.shape()[0] == derivatives + 1); debug_assert!(data.shape()[1] == degree + 1); @@ -69,7 +69,7 @@ fn tabulate_interval< for i in 0..data.shape()[2] { let d = *data.get([k, p - 1, i]).unwrap(); *data.get_mut([k, p, i]).unwrap() = - (T::from(*points.get([0, i]).unwrap()).unwrap() * T::from(2.0).unwrap() + (T::from(points.get_value([0, i]).unwrap()).unwrap() * T::from(2.0).unwrap() - T::from(1.0).unwrap()) * d * b; @@ -98,13 +98,13 @@ fn tabulate_interval< /// Tabulate orthonormal polynomials on a quadrilateral fn tabulate_quadrilateral< T: RlstScalar, - Array2: RandomAccessByRef<2, Item = T::Real> + Shape<2>, - Array3Mut: RandomAccessMut<3, Item = T> + RandomAccessByRef<3, Item = T> + Shape<3>, + Array2Impl: ValueArrayImpl, + Array3MutImpl: MutableArrayImpl, >( - points: &Array, + points: &Array, degree: usize, derivatives: usize, - data: &mut Array, + data: &mut Array, ) { debug_assert!(data.shape()[0] == (derivatives + 1) * (derivatives + 2) / 2); debug_assert!(data.shape()[1] == (degree + 1) * (degree + 1)); @@ -139,7 +139,7 @@ fn tabulate_quadrilateral< .unwrap(); *data .get_mut([tri_index(k, 0), quad_index(p, 0, degree), i]) - .unwrap() = (T::from(*points.get([0, i]).unwrap()).unwrap() + .unwrap() = (T::from(points.get_value([0, i]).unwrap()).unwrap() * T::from(2.0).unwrap() - T::from(1.0).unwrap()) * d @@ -194,7 +194,7 @@ fn tabulate_quadrilateral< .unwrap(); *data .get_mut([tri_index(0, k), quad_index(0, p, degree), i]) - .unwrap() = (T::from(*points.get([1, i]).unwrap()).unwrap() + .unwrap() = (T::from(points.get_value([1, i]).unwrap()).unwrap() * T::from(2.0).unwrap() - T::from(1.0).unwrap()) * d @@ -251,13 +251,13 @@ fn tabulate_quadrilateral< /// Tabulate orthonormal polynomials on a triangle fn tabulate_triangle< T: RlstScalar, - Array2: RandomAccessByRef<2, Item = T::Real> + Shape<2>, - Array3Mut: RandomAccessMut<3, Item = T> + RandomAccessByRef<3, Item = T> + Shape<3>, + Array2Impl: ValueArrayImpl, + Array3MutImpl: MutableArrayImpl, >( - points: &Array, + points: &Array, degree: usize, derivatives: usize, - data: &mut Array, + data: &mut Array, ) { debug_assert!(data.shape()[0] == (derivatives + 1) * (derivatives + 2) / 2); debug_assert!(data.shape()[1] == (degree + 1) * (degree + 2) / 2); @@ -290,9 +290,9 @@ fn tabulate_triangle< .unwrap(); *data .get_mut([tri_index(kx, ky), tri_index(0, p), i]) - .unwrap() = (T::from(*points.get([0, i]).unwrap()).unwrap() + .unwrap() = (T::from(points.get_value([0, i]).unwrap()).unwrap() * T::from(2.0).unwrap() - + T::from(*points.get([1, i]).unwrap()).unwrap() + + T::from(points.get_value([1, i]).unwrap()).unwrap() - T::from(1.0).unwrap()) * d * a @@ -329,8 +329,8 @@ fn tabulate_triangle< ); for i in 0..data.shape()[2] { - let b = - T::from(1.0).unwrap() - T::from(*points.get([1, i]).unwrap()).unwrap(); + let b = T::from(1.0).unwrap() + - T::from(points.get_value([1, i]).unwrap()).unwrap(); let d = *data .get([tri_index(kx, ky), tri_index(0, p - 2), i]) .unwrap(); @@ -347,7 +347,7 @@ fn tabulate_triangle< .get_mut([tri_index(kx, ky), tri_index(0, p), i]) .unwrap() -= T::from(2.0).unwrap() * T::from(ky).unwrap() - * (T::from(*points.get([1, i]).unwrap()).unwrap() + * (T::from(points.get_value([1, i]).unwrap()).unwrap() - T::from(1.0).unwrap()) * d * scale2 @@ -380,7 +380,7 @@ fn tabulate_triangle< .get_mut([tri_index(kx, ky), tri_index(1, p), i]) .unwrap() = *data.get([tri_index(kx, ky), tri_index(0, p), i]).unwrap() * scale3 - * ((T::from(*points.get([1, i]).unwrap()).unwrap() + * ((T::from(points.get_value([1, i]).unwrap()).unwrap() * T::from(2.0).unwrap() - T::from(1.0).unwrap()) * (T::from(1.5).unwrap() + T::from(p).unwrap()) @@ -418,7 +418,7 @@ fn tabulate_triangle< .get_mut([tri_index(kx, ky), tri_index(q + 1, p), i]) .unwrap() = d * scale4 - * ((T::from(*points.get([1, i]).unwrap()).unwrap() + * ((T::from(points.get_value([1, i]).unwrap()).unwrap() * T::from(T::from(2.0).unwrap()).unwrap() - T::from(T::from(1.0).unwrap()).unwrap()) * a1 @@ -452,13 +452,13 @@ fn tabulate_triangle< /// Tabulate orthonormal polynomials on a tetrahedron fn tabulate_tetrahedron< T: RlstScalar, - Array2: RandomAccessByRef<2, Item = T::Real> + Shape<2>, - Array3Mut: RandomAccessMut<3, Item = T> + RandomAccessByRef<3, Item = T> + Shape<3>, + Array2Impl: ValueArrayImpl, + Array3MutImpl: MutableArrayImpl, >( - points: &Array, + points: &Array, degree: usize, derivatives: usize, - data: &mut Array, + data: &mut Array, ) { debug_assert!(data.shape()[0] == (derivatives + 1) * (derivatives + 2) * (derivatives + 3) / 6); debug_assert!(data.shape()[1] == (degree + 1) * (degree + 2) * (degree + 3) / 6); @@ -488,10 +488,10 @@ fn tabulate_tetrahedron< .unwrap(); *data .get_mut([tet_index(kx, ky, kz), tet_index(0, 0, p), i]) - .unwrap() = (T::from(*points.get([0, i]).unwrap()).unwrap() + .unwrap() = (T::from(points.get_value([0, i]).unwrap()).unwrap() * T::from(2.0).unwrap() - + T::from(*points.get([1, i]).unwrap()).unwrap() - + T::from(*points.get([2, i]).unwrap()).unwrap() + + T::from(points.get_value([1, i]).unwrap()).unwrap() + + T::from(points.get_value([2, i]).unwrap()).unwrap() - T::from(1.0).unwrap()) * a * d; @@ -534,7 +534,8 @@ fn tabulate_tetrahedron< *data .get_mut([tet_index(kx, ky, kz), tet_index(0, 0, p), i]) .unwrap() -= (T::from( - *points.get([1, i]).unwrap() + *points.get([2, i]).unwrap(), + points.get_value([1, i]).unwrap() + + points.get_value([2, i]).unwrap(), ) .unwrap() - T::from(1.0).unwrap()) @@ -551,7 +552,8 @@ fn tabulate_tetrahedron< .get_mut([tet_index(kx, ky, kz), tet_index(0, 0, p), i]) .unwrap() -= T::from(ky * 2).unwrap() * (T::from( - *points.get([1, i]).unwrap() + *points.get([2, i]).unwrap(), + points.get_value([1, i]).unwrap() + + points.get_value([2, i]).unwrap(), ) .unwrap() - T::from(1.0).unwrap()) @@ -580,7 +582,8 @@ fn tabulate_tetrahedron< .get_mut([tet_index(kx, ky, kz), tet_index(0, 0, p), i]) .unwrap() -= T::from(kz * 2).unwrap() * (T::from( - *points.get([1, i]).unwrap() + *points.get([2, i]).unwrap(), + points.get_value([1, i]).unwrap() + + points.get_value([2, i]).unwrap(), ) .unwrap() - T::from(1.0).unwrap()) @@ -621,9 +624,9 @@ fn tabulate_tetrahedron< *data .get_mut([tet_index(kx, ky, kz), tet_index(0, 1, p), i]) .unwrap() = d - * (T::from(*points.get([1, i]).unwrap()).unwrap() + * (T::from(points.get_value([1, i]).unwrap()).unwrap() * T::from(2 * p + 3).unwrap() - + T::from(*points.get([2, i]).unwrap()).unwrap() + + T::from(points.get_value([2, i]).unwrap()).unwrap() - T::from(1).unwrap()); } if ky > 0 { @@ -663,15 +666,15 @@ fn tabulate_tetrahedron< .get_mut([tet_index(kx, ky, kz), tet_index(0, q + 1, p), i]) .unwrap() = d * (aq - * (T::from(*points.get([1, i]).unwrap()).unwrap() + * (T::from(points.get_value([1, i]).unwrap()).unwrap() * T::from(2.0).unwrap() - T::from(1.0).unwrap() - + T::from(*points.get([2, i]).unwrap()).unwrap()) + + T::from(points.get_value([2, i]).unwrap()).unwrap()) + bq * (T::from(1.0).unwrap() - - T::from(*points.get([2, i]).unwrap()).unwrap())) + - T::from(points.get_value([2, i]).unwrap()).unwrap())) - d2 * cq * (T::from(1.0).unwrap() - - T::from(*points.get([2, i]).unwrap()).unwrap()) + - T::from(points.get_value([2, i]).unwrap()).unwrap()) .powi(2); } @@ -698,7 +701,7 @@ fn tabulate_tetrahedron< .unwrap() += T::from(kz).unwrap() * d * (aq - bq) + T::from(2 * kz).unwrap() * (T::from(1.0).unwrap() - - T::from(*points.get([2, i]).unwrap()).unwrap()) + - T::from(points.get_value([2, i]).unwrap()).unwrap()) * d2 * cq; } @@ -725,7 +728,7 @@ fn tabulate_tetrahedron< *data .get_mut([tet_index(kx, ky, kz), tet_index(1, q, p), i]) .unwrap() = d - * (T::from(*points.get([2, i]).unwrap()).unwrap() + * (T::from(points.get_value([2, i]).unwrap()).unwrap() * T::from(2 + p + q).unwrap() * T::from(2.0).unwrap() - T::from(1.0).unwrap()); @@ -761,7 +764,8 @@ fn tabulate_tetrahedron< .unwrap() = d * (ar * (T::from(2.0).unwrap() - * T::from(*points.get([2, i]).unwrap()).unwrap() + * T::from(points.get_value([2, i]).unwrap()) + .unwrap() - T::from(1.0).unwrap()) + br) - d2 * cr; @@ -811,13 +815,13 @@ fn tabulate_tetrahedron< /// Tabulate orthonormal polynomials on a hexahedron fn tabulate_hexahedron< T: RlstScalar, - Array2: RandomAccessByRef<2, Item = T::Real> + Shape<2>, - Array3Mut: RandomAccessMut<3, Item = T> + RandomAccessByRef<3, Item = T> + Shape<3>, + Array2Impl: ValueArrayImpl, + Array3MutImpl: MutableArrayImpl, >( - points: &Array, + points: &Array, degree: usize, derivatives: usize, - data: &mut Array, + data: &mut Array, ) { debug_assert!(data.shape()[0] == (derivatives + 1) * (derivatives + 2) * (derivatives + 3) / 6); debug_assert!(data.shape()[1] == (degree + 1) * (degree + 1) * (degree + 1)); @@ -852,7 +856,7 @@ fn tabulate_hexahedron< .unwrap(); *data .get_mut([tet_index(k, 0, 0), hex_index(p, 0, 0, degree), i]) - .unwrap() = (T::from(*points.get([0, i]).unwrap()).unwrap() + .unwrap() = (T::from(points.get_value([0, i]).unwrap()).unwrap() * T::from(2.0).unwrap() - T::from(1.0).unwrap()) * d @@ -907,7 +911,7 @@ fn tabulate_hexahedron< .unwrap(); *data .get_mut([tet_index(0, k, 0), hex_index(0, p, 0, degree), i]) - .unwrap() = (T::from(*points.get([1, i]).unwrap()).unwrap() + .unwrap() = (T::from(points.get_value([1, i]).unwrap()).unwrap() * T::from(2.0).unwrap() - T::from(1.0).unwrap()) * d @@ -962,7 +966,7 @@ fn tabulate_hexahedron< .unwrap(); *data .get_mut([tet_index(0, 0, k), hex_index(0, 0, p, degree), i]) - .unwrap() = (T::from(*points.get([2, i]).unwrap()).unwrap() + .unwrap() = (T::from(points.get_value([2, i]).unwrap()).unwrap() * T::from(2.0).unwrap() - T::from(1.0).unwrap()) * d @@ -1029,9 +1033,9 @@ fn tabulate_hexahedron< } /// The shape of a table containing the values of Legendre polynomials -pub fn shape + Shape<2>>( +pub fn shape>( cell_type: ReferenceCellType, - points: &Array, + points: &Array, degree: usize, derivatives: usize, ) -> [usize; 3] { @@ -1045,8 +1049,8 @@ pub fn shape + Shape<2>>( /// Tabulate orthonormal polynomials pub fn tabulate< T: RlstScalar, - Array2: RandomAccessByRef<2, Item = T::Real> + Shape<2>, - Array3Mut: RandomAccessMut<3, Item = T> + RandomAccessByRef<3, Item = T> + Shape<3>, + Array2: ValueArrayImpl, + Array3Mut: MutableArrayImpl, >( cell_type: ReferenceCellType, points: &Array, diff --git a/src/traits.rs b/src/traits.rs index 90f706e..76c9bd3 100644 --- a/src/traits.rs +++ b/src/traits.rs @@ -1,9 +1,6 @@ //! Traits use crate::types::DofTransformation; -use rlst::{ - Array, RandomAccessByRef, RandomAccessMut, RlstScalar, Shape, UnsafeRandomAccessByRef, - UnsafeRandomAccessMut, -}; +use rlst::{Array, MutableArrayImpl, RlstScalar, ValueArrayImpl}; use std::fmt::Debug; use std::hash::Hash; @@ -33,13 +30,13 @@ pub trait FiniteElement { /// Tabulate the values of the basis functions and their derivatives at a set of points fn tabulate< - Array2: RandomAccessByRef<2, Item = ::Real> + Shape<2>, - Array4Mut: RandomAccessMut<4, Item = Self::T> + UnsafeRandomAccessMut<4, Item = Self::T>, + Array2Impl: ValueArrayImpl<::Real, 2>, + Array4MutImpl: MutableArrayImpl, >( &self, - points: &Array, + points: &Array, nderivs: usize, - data: &mut Array, + data: &mut Array, ); /// The DOFs that are associated with a subentity of the reference cell @@ -53,32 +50,32 @@ pub trait FiniteElement { /// Push function values forward to a physical cell fn push_forward< - Array3Real: UnsafeRandomAccessByRef<3, Item = ::Real> + Shape<3>, - Array4: UnsafeRandomAccessByRef<4, Item = Self::T> + Shape<4>, - Array4Mut: UnsafeRandomAccessMut<4, Item = Self::T> + Shape<4>, + Array3RealImpl: ValueArrayImpl<::Real, 3>, + Array4Impl: ValueArrayImpl, + Array4MutImpl: MutableArrayImpl, >( &self, - reference_values: &Array, + reference_values: &Array, nderivs: usize, - jacobians: &Array, + jacobians: &Array, jacobian_determinants: &[::Real], - inverse_jacobians: &Array, - physical_values: &mut Array, + inverse_jacobians: &Array, + physical_values: &mut Array, ); /// Pull function values back to the reference cell fn pull_back< - Array3Real: UnsafeRandomAccessByRef<3, Item = ::Real> + Shape<3>, - Array4: UnsafeRandomAccessByRef<4, Item = Self::T> + Shape<4>, - Array4Mut: UnsafeRandomAccessMut<4, Item = Self::T> + Shape<4>, + Array3RealImpl: ValueArrayImpl<::Real, 3>, + Array4Impl: ValueArrayImpl, + Array4MutImpl: MutableArrayImpl, >( &self, - physical_values: &Array, + physical_values: &Array, nderivs: usize, - jacobians: &Array, + jacobians: &Array, jacobian_determinants: &[::Real], - inverse_jacobians: &Array, - reference_values: &mut Array, + inverse_jacobians: &Array, + reference_values: &mut Array, ); /// The value shape on a physical cell @@ -150,33 +147,33 @@ pub trait Map { /// Push function values forward to a physical cell fn push_forward< T: RlstScalar, - Array3Real: UnsafeRandomAccessByRef<3, Item = T::Real> + Shape<3>, - Array4: UnsafeRandomAccessByRef<4, Item = T> + Shape<4>, - Array4Mut: UnsafeRandomAccessMut<4, Item = T> + Shape<4>, + Array3RealImpl: ValueArrayImpl, + Array4Impl: ValueArrayImpl, + Array4MutImpl: MutableArrayImpl, >( &self, - reference_values: &Array, + reference_values: &Array, nderivs: usize, - jacobians: &Array, + jacobians: &Array, jacobian_determinants: &[T::Real], - inverse_jacobians: &Array, - physical_values: &mut Array, + inverse_jacobians: &Array, + physical_values: &mut Array, ); /// Pull function values back to the reference cell fn pull_back< T: RlstScalar, - Array3Real: UnsafeRandomAccessByRef<3, Item = T::Real> + Shape<3>, - Array4: UnsafeRandomAccessByRef<4, Item = T> + Shape<4>, - Array4Mut: UnsafeRandomAccessMut<4, Item = T> + Shape<4>, + Array3RealImpl: ValueArrayImpl, + Array4Impl: ValueArrayImpl, + Array4MutImpl: MutableArrayImpl, >( &self, - physical_values: &Array, + physical_values: &Array, nderivs: usize, - jacobians: &Array, + jacobians: &Array, jacobian_determinants: &[T::Real], - inverse_jacobians: &Array, - reference_values: &mut Array, + inverse_jacobians: &Array, + reference_values: &mut Array, ); /// The value shape on a physical cell From 69784f91a77dda18797b89c1a8e8379a98a241dd Mon Sep 17 00:00:00 2001 From: Timo Betcke Date: Thu, 6 Nov 2025 18:15:28 +0000 Subject: [PATCH 2/4] More traits clean up --- src/math.rs | 37 +++++++++---------------------------- 1 file changed, 9 insertions(+), 28 deletions(-) diff --git a/src/math.rs b/src/math.rs index 344d34e..7de6c7f 100644 --- a/src/math.rs +++ b/src/math.rs @@ -1,15 +1,9 @@ //! Mathematical functions -use rlst::{ - Array, RandomAccessByRef, RandomAccessMut, RlstScalar, Shape, UnsafeRandomAccessByRef, - UnsafeRandomAccessMut, ValueArrayImpl, -}; +use rlst::{Array, MutableArrayImpl, RlstScalar, ValueArrayImpl}; /// Orthogonalise the rows of a matrix, starting with the row numbered `start` -pub fn orthogonalise< - T: RlstScalar, - Array2Mut: RandomAccessByRef<2, Item = T> + RandomAccessMut<2, Item = T> + Shape<2>, ->( - mat: &mut Array, +pub fn orthogonalise>( + mat: &mut Array, start: usize, ) { for row in start..mat.shape()[0] { @@ -33,11 +27,8 @@ pub fn orthogonalise< } /// Orthogonalise the rows of a matrix, starting with the row numbered `start` -pub fn orthogonalise_3< - T: RlstScalar, - Array3Mut: RandomAccessByRef<3, Item = T> + RandomAccessMut<3, Item = T> + Shape<3>, ->( - mat: &mut Array, +pub fn orthogonalise_3>( + mat: &mut Array, start: usize, ) { for row in start..mat.shape()[0] { @@ -73,11 +64,7 @@ pub fn orthogonalise_3< } /// Swap two entries in a matrix -unsafe fn entry_swap< - const N: usize, - T: RlstScalar, - ArrayMut: UnsafeRandomAccessMut + UnsafeRandomAccessByRef + Shape, ->( +unsafe fn entry_swap>( mat: &mut Array, mindex0: [usize; N], mindex1: [usize; N], @@ -90,11 +77,8 @@ unsafe fn entry_swap< } /// Compute the LU decomposition of the transpose of a square matrix -pub fn lu_transpose< - T: RlstScalar, - Array2Mut: UnsafeRandomAccessMut<2, Item = T> + UnsafeRandomAccessByRef<2, Item = T> + Shape<2>, ->( - mat: &mut Array, +pub fn lu_transpose>( + mat: &mut Array, ) -> Vec { let dim = mat.shape()[0]; assert_eq!(mat.shape()[1], dim); @@ -155,10 +139,7 @@ pub fn apply_permutation(perm: &[usize], data: &mut [T]) { } /// Convert a linear transformation info the format expected by `apply_matrix` and return the premutation to pass into `apply_matrix` -pub fn prepare_matrix< - T: RlstScalar, - Array2Mut: UnsafeRandomAccessMut<2, Item = T> + UnsafeRandomAccessByRef<2, Item = T> + Shape<2>, ->( +pub fn prepare_matrix>( mat: &mut Array, ) -> Vec { let mut perm = lu_transpose(mat); From 7f6030742b32b7ba13b4d5b2bdd0553f396dca9e Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Fri, 7 Nov 2025 07:06:15 +0000 Subject: [PATCH 3/4] external crate blas_src, consistent template naming --- src/ciarlet.rs | 9 +++++---- src/map.rs | 8 ++++---- src/math.rs | 4 ++-- src/polynomials/legendre.rs | 8 ++++---- 4 files changed, 15 insertions(+), 14 deletions(-) diff --git a/src/ciarlet.rs b/src/ciarlet.rs index e95c945..22197b0 100644 --- a/src/ciarlet.rs +++ b/src/ciarlet.rs @@ -1,5 +1,6 @@ //! Finite element definitions +extern crate blas_src; extern crate lapack_src; use crate::math; @@ -609,13 +610,13 @@ impl FiniteElement for CiarletElement { self.dim } fn tabulate< - Array2: ValueArrayImpl<::Real, 2>, - Array4Mut: MutableArrayImpl, + Array2Impl: ValueArrayImpl<::Real, 2>, + Array4MutImpl: MutableArrayImpl, >( &self, - points: &Array, + points: &Array, nderivs: usize, - data: &mut Array, + data: &mut Array, ) { let mut table = DynArray::::from_shape(legendre_shape( self.cell_type, diff --git a/src/map.rs b/src/map.rs index 648e21d..f8567c3 100644 --- a/src/map.rs +++ b/src/map.rs @@ -235,12 +235,12 @@ mod test { fn fill_jacobians< T: RlstScalar, - Array3: RandomAccessMut<3, Item = T>, - Array3Mut: RandomAccessMut<3, Item = T>, + Array3Impl: ValueArrayImpl<3, Item = T>, + Array3MutImpl: MutableArrayImpl<3, Item = T>, >( - j: &mut Array, + j: &mut Array, jdet: &mut [T], - jinv: &mut Array, + jinv: &mut Array, ) { *j.get_mut([0, 0, 0]).unwrap() = T::from(1.0).unwrap(); *j.get_mut([0, 0, 1]).unwrap() = T::from(1.0).unwrap(); diff --git a/src/math.rs b/src/math.rs index 7de6c7f..aad6236 100644 --- a/src/math.rs +++ b/src/math.rs @@ -158,8 +158,8 @@ pub fn apply_perm_and_matrix>( } /// Apply a matrix to some data -pub fn apply_matrix>( - mat: &Array, +pub fn apply_matrix>( + mat: &Array, data: &mut [T], ) { let dim = mat.shape()[0]; diff --git a/src/polynomials/legendre.rs b/src/polynomials/legendre.rs index f6267af..96e2ead 100644 --- a/src/polynomials/legendre.rs +++ b/src/polynomials/legendre.rs @@ -1049,14 +1049,14 @@ pub fn shape>( /// Tabulate orthonormal polynomials pub fn tabulate< T: RlstScalar, - Array2: ValueArrayImpl, - Array3Mut: MutableArrayImpl, + Array2Impl: ValueArrayImpl, + Array3MutImpl: MutableArrayImpl, >( cell_type: ReferenceCellType, - points: &Array, + points: &Array, degree: usize, derivatives: usize, - data: &mut Array, + data: &mut Array, ) { match cell_type { ReferenceCellType::Interval => tabulate_interval(points, degree, derivatives, data), From b684b980e6c3901df84c7fa23859bf4ad1fd2404 Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Fri, 7 Nov 2025 07:19:02 +0000 Subject: [PATCH 4/4] correct --- src/map.rs | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/src/map.rs b/src/map.rs index f8567c3..1ab6866 100644 --- a/src/map.rs +++ b/src/map.rs @@ -231,14 +231,10 @@ impl Map for ContravariantPiolaMap { mod test { use super::*; use approx::*; - use rlst::{Array, RandomAccessMut, rlst_dynamic_array}; + use rlst::{Array, rlst_dynamic_array}; - fn fill_jacobians< - T: RlstScalar, - Array3Impl: ValueArrayImpl<3, Item = T>, - Array3MutImpl: MutableArrayImpl<3, Item = T>, - >( - j: &mut Array, + fn fill_jacobians, Array3MutImpl: MutableArrayImpl>( + j: &mut Array, jdet: &mut [T], jinv: &mut Array, ) {