From 4b5affac367c9f47315319567c1f9157e3299adb Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Tue, 17 Oct 2023 14:13:36 +0100 Subject: [PATCH] get element working without Array2D --- bem/src/assembly/batched.rs | 23 ++-- element/src/element.rs | 18 +-- element/src/element/lagrange.rs | 41 +++---- element/src/element/raviart_thomas.rs | 17 +-- element/src/polynomials.rs | 42 +++---- grid/src/grid.rs | 34 ++++++ tools/Cargo.toml | 1 + tools/examples/array2d.rs | 52 --------- tools/src/arrays.rs | 151 +++----------------------- traits/src/arrays.rs | 65 ----------- traits/src/element.rs | 7 +- traits/src/grid.rs | 76 +++++++------ 12 files changed, 166 insertions(+), 361 deletions(-) delete mode 100644 tools/examples/array2d.rs diff --git a/bem/src/assembly/batched.rs b/bem/src/assembly/batched.rs index 7392d9f6..bd53246e 100644 --- a/bem/src/assembly/batched.rs +++ b/bem/src/assembly/batched.rs @@ -215,7 +215,16 @@ fn assemble_batch_nonadjacent<'a, const NPTS_TEST: usize, const NPTS_TRIAL: usiz let gdim = test_grid.geometry().dim(); // TODO: move this to grid.get_compute_points_function(test_points) - let test_compute_points = |cell: usize, pts: &mut rlst_dense::Matrix, rlst_dense::Dynamic>, rlst_dense::Dynamic>| { + let test_compute_points = |cell: usize, + pts: &mut rlst_dense::Matrix< + f64, + rlst_dense::base_matrix::BaseMatrix< + f64, + rlst_dense::VectorContainer, + rlst_dense::Dynamic, + >, + rlst_dense::Dynamic, + >| { for p in 0..NPTS_TEST { for i in 0..gdim { unsafe { @@ -454,8 +463,8 @@ pub fn assemble_nonsingular<'a, const NPTS_TEST: usize, const NPTS_TRIAL: usize> test_start = test_end } - let numthreads = test_cells.len(); - let r: usize = (0..numthreads) + let numtasks = test_cells.len(); + let r: usize = (0..numtasks) .into_par_iter() .map(&|t| { assemble_batch_nonadjacent::( @@ -476,7 +485,7 @@ pub fn assemble_nonsingular<'a, const NPTS_TEST: usize, const NPTS_TRIAL: usize> ) }) .sum(); - assert_eq!(r, numthreads); + assert_eq!(r, numtasks); } } } @@ -616,8 +625,8 @@ pub fn assemble_singular<'a>( start = end; } - let numthreads = cell_blocks.len(); - let r: usize = (0..numthreads) + let numtasks = cell_blocks.len(); + let r: usize = (0..numtasks) .into_par_iter() .map(&|t| { assemble_batch_singular( @@ -636,7 +645,7 @@ pub fn assemble_singular<'a>( ) }) .sum(); - assert_eq!(r, numthreads); + assert_eq!(r, numtasks); } } } diff --git a/element/src/element.rs b/element/src/element.rs index 503c6205..32269d33 100644 --- a/element/src/element.rs +++ b/element/src/element.rs @@ -2,13 +2,13 @@ use crate::cell::create_cell; use crate::polynomials::{legendre_shape, polynomial_count, tabulate_legendre_polynomials}; -use bempp_tools::arrays::{AdjacencyList, Array2D, Array3D}; -use bempp_traits::arrays::{AdjacencyListAccess, Array2DAccess, Array3DAccess, Array4DAccess}; +use bempp_tools::arrays::{AdjacencyList, Array3D, Mat}; +use bempp_traits::arrays::{AdjacencyListAccess, Array3DAccess, Array4DAccess}; use bempp_traits::cell::ReferenceCellType; use bempp_traits::element::{Continuity, ElementFamily, FiniteElement, MapType}; use rlst_algorithms::linalg::LinAlg; use rlst_algorithms::traits::inverse::Inverse; -use rlst_dense::{RandomAccessByRef, RandomAccessMut}; +use rlst_dense::{rlst_dynamic_mat, RandomAccessByRef, RandomAccessMut, Shape}; pub mod lagrange; pub mod raviart_thomas; @@ -37,7 +37,7 @@ impl CiarletElement { degree: usize, value_shape: Vec, polynomial_coeffs: Array3D, - interpolation_points: [Vec>; 4], + interpolation_points: [Vec>; 4], interpolation_weights: [Vec>; 4], map_type: MapType, continuity: Continuity, @@ -69,12 +69,12 @@ impl CiarletElement { } let new_pts = if continuity == Continuity::Discontinuous { - let mut new_pts = [vec![], vec![], vec![], vec![]]; + let mut new_pts: [Vec>; 4] = [vec![], vec![], vec![], vec![]]; let mut pn = 0; - let mut all_pts = Array2D::::new((npts, tdim)); + let mut all_pts = rlst_dynamic_mat![f64, (npts, tdim)]; for (i, pts_i) in interpolation_points.iter().take(tdim).enumerate() { for _pts in pts_i { - new_pts[i].push(Array2D::::new((0, tdim))); + new_pts[i].push(rlst_dynamic_mat![f64, (0, tdim)]); } } for pts_i in interpolation_points.iter() { @@ -240,9 +240,9 @@ impl FiniteElement for CiarletElement { fn dim(&self) -> usize { self.dim } - fn tabulate<'a>( + fn tabulate + Shape>( &self, - points: &impl Array2DAccess<'a, f64>, + points: &T, nderivs: usize, data: &mut impl Array4DAccess, ) { diff --git a/element/src/element/lagrange.rs b/element/src/element/lagrange.rs index d41d0e86..29178e54 100644 --- a/element/src/element/lagrange.rs +++ b/element/src/element/lagrange.rs @@ -2,7 +2,7 @@ use crate::element::{create_cell, CiarletElement}; use crate::polynomials::polynomial_count; -use bempp_tools::arrays::{Array2D, Array3D}; +use bempp_tools::arrays::{to_matrix, Array3D}; use bempp_traits::arrays::Array3DAccess; use bempp_traits::cell::ReferenceCellType; use bempp_traits::element::{Continuity, ElementFamily, MapType}; @@ -29,18 +29,18 @@ pub fn create( } for d in 0..tdim { for _e in 0..cell.entity_count(d) { - x[d].push(Array2D::::new((0, tdim))); + x[d].push(to_matrix(&[], (0, tdim))); m[d].push(Array3D::::new((0, 1, 0))); } } - x[tdim].push(Array2D::::from_data(cell.midpoint(), (1, tdim))); + x[tdim].push(to_matrix(&cell.midpoint(), (1, tdim))); m[tdim].push(Array3D::::from_data(vec![1.0], (1, 1, 1))); } else { // TODO: GLL points for e in 0..cell.entity_count(0) { let mut pts = vec![0.0; tdim]; pts.copy_from_slice(&cell.vertices()[e * tdim..(e + 1) * tdim]); - x[0].push(Array2D::::from_data(pts, (1, tdim))); + x[0].push(to_matrix(&pts, (1, tdim))); m[0].push(Array3D::::from_data(vec![1.0], (1, 1, 1))); } for e in 0..cell.entity_count(1) { @@ -56,7 +56,7 @@ pub fn create( pts[(i - 1) * tdim + j] = v0[j] + i as f64 / degree as f64 * (v1[j] - v0[j]); } } - x[1].push(Array2D::::from_data(pts, (degree - 1, tdim))); + x[1].push(to_matrix(&pts, (degree - 1, tdim))); m[1].push(Array3D::::from_data( ident, (degree - 1, 1, degree - 1), @@ -119,7 +119,7 @@ pub fn create( for i in 0..npts { ident[i * npts + i] = 1.0; } - x[2].push(Array2D::::from_data(pts, (npts, tdim))); + x[2].push(to_matrix(&pts, (npts, tdim))); m[2].push(Array3D::::from_data(ident, (npts, 1, npts))); start += nvertices; } @@ -143,9 +143,10 @@ mod test { use crate::cell::*; use crate::element::lagrange::*; use approx::*; - use bempp_tools::arrays::{Array2D, Array4D}; - use bempp_traits::arrays::{Array2DAccess, Array4DAccess}; + use bempp_tools::arrays::Array4D; + use bempp_traits::arrays::Array4DAccess; use bempp_traits::element::FiniteElement; + use rlst_dense::RandomAccessByRef; fn check_dofs(e: impl FiniteElement) { let cell_dim = match e.cell_type() { @@ -182,7 +183,7 @@ mod test { let e = create(ReferenceCellType::Interval, 0, Continuity::Discontinuous); assert_eq!(e.value_size(), 1); let mut data = Array4D::::new(e.tabulate_array_shape(0, 4)); - let points = Array2D::from_data(vec![0.0, 0.2, 0.4, 1.0], (4, 1)); + let points = to_matrix(&vec![0.0, 0.2, 0.4, 1.0], (4, 1)); e.tabulate(&points, 0, &mut data); for pt in 0..4 { @@ -196,7 +197,7 @@ mod test { let e = create(ReferenceCellType::Interval, 1, Continuity::Continuous); assert_eq!(e.value_size(), 1); let mut data = Array4D::::new(e.tabulate_array_shape(0, 4)); - let points = Array2D::from_data(vec![0.0, 0.2, 0.4, 1.0], (4, 1)); + let points = to_matrix(&vec![0.0, 0.2, 0.4, 1.0], (4, 1)); e.tabulate(&points, 0, &mut data); for pt in 0..4 { @@ -214,8 +215,8 @@ mod test { let e = create(ReferenceCellType::Triangle, 0, Continuity::Discontinuous); assert_eq!(e.value_size(), 1); let mut data = Array4D::::new(e.tabulate_array_shape(0, 6)); - let points = Array2D::from_data( - vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.0, 0.5, 0.5, 0.5], + let points = to_matrix( + &vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.0, 0.5, 0.5, 0.5], (6, 2), ); e.tabulate(&points, 0, &mut data); @@ -231,8 +232,8 @@ mod test { let e = create(ReferenceCellType::Triangle, 1, Continuity::Continuous); assert_eq!(e.value_size(), 1); let mut data = Array4D::::new(e.tabulate_array_shape(0, 6)); - let points = Array2D::from_data( - vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.0, 0.5, 0.5, 0.5], + let points = to_matrix( + &vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.0, 0.5, 0.5, 0.5], (6, 2), ); e.tabulate(&points, 0, &mut data); @@ -312,8 +313,8 @@ mod test { ); assert_eq!(e.value_size(), 1); let mut data = Array4D::::new(e.tabulate_array_shape(0, 6)); - let points = Array2D::from_data( - vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.25, 0.5, 0.3, 0.2], + let points = to_matrix( + &vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.25, 0.5, 0.3, 0.2], (6, 2), ); e.tabulate(&points, 0, &mut data); @@ -329,8 +330,8 @@ mod test { let e = create(ReferenceCellType::Quadrilateral, 1, Continuity::Continuous); assert_eq!(e.value_size(), 1); let mut data = Array4D::::new(e.tabulate_array_shape(0, 6)); - let points = Array2D::from_data( - vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.25, 0.5, 0.3, 0.2], + let points = to_matrix( + &vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.25, 0.5, 0.3, 0.2], (6, 2), ); e.tabulate(&points, 0, &mut data); @@ -361,8 +362,8 @@ mod test { let e = create(ReferenceCellType::Quadrilateral, 2, Continuity::Continuous); assert_eq!(e.value_size(), 1); let mut data = Array4D::::new(e.tabulate_array_shape(0, 6)); - let points = Array2D::from_data( - vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.25, 0.5, 0.3, 0.2], + let points = to_matrix( + &vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.25, 0.5, 0.3, 0.2], (6, 2), ); e.tabulate(&points, 0, &mut data); diff --git a/element/src/element/raviart_thomas.rs b/element/src/element/raviart_thomas.rs index 2964bc7a..b2331bd4 100644 --- a/element/src/element/raviart_thomas.rs +++ b/element/src/element/raviart_thomas.rs @@ -2,7 +2,7 @@ use crate::element::{create_cell, CiarletElement}; use crate::polynomials::polynomial_count; -use bempp_tools::arrays::{Array2D, Array3D}; +use bempp_tools::arrays::{to_matrix, Array3D}; use bempp_traits::arrays::Array3DAccess; use bempp_traits::cell::ReferenceCellType; use bempp_traits::element::{Continuity, ElementFamily, MapType}; @@ -45,7 +45,7 @@ pub fn create( let mut x = [vec![], vec![], vec![], vec![]]; let mut m = [vec![], vec![], vec![], vec![]]; for _e in 0..cell.entity_count(0) { - x[0].push(Array2D::::new((0, tdim))); + x[0].push(to_matrix(&[], (0, tdim))); m[0].push(Array3D::::new((0, 2, 0))); } @@ -61,12 +61,12 @@ pub fn create( } mat[0] = v0[1] - v1[1]; mat[1] = v1[0] - v0[0]; - x[1].push(Array2D::::from_data(pts, (1, tdim))); + x[1].push(to_matrix(&pts, (1, tdim))); m[1].push(Array3D::::from_data(mat, (1, 2, 1))); } for _e in 0..cell.entity_count(2) { - x[2].push(Array2D::::new((0, tdim))); + x[2].push(to_matrix(&[], (0, tdim))); m[2].push(Array3D::::new((0, 2, 0))); } @@ -89,9 +89,10 @@ mod test { use crate::cell::*; use crate::element::raviart_thomas::*; use approx::*; - use bempp_tools::arrays::{Array2D, Array4D}; - use bempp_traits::arrays::{Array2DAccess, Array4DAccess}; + use bempp_tools::arrays::Array4D; + use bempp_traits::arrays::Array4DAccess; use bempp_traits::element::FiniteElement; + use rlst_dense::RandomAccessByRef; fn check_dofs(e: impl FiniteElement) { let cell_dim = match e.cell_type() { @@ -128,8 +129,8 @@ mod test { let e = create(ReferenceCellType::Triangle, 1, Continuity::Continuous); assert_eq!(e.value_size(), 2); let mut data = Array4D::::new(e.tabulate_array_shape(0, 6)); - let points = Array2D::from_data( - vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.0, 0.5, 0.5, 0.5], + let points = to_matrix( + &vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.0, 0.5, 0.5, 0.5], (6, 2), ); e.tabulate(&points, 0, &mut data); diff --git a/element/src/polynomials.rs b/element/src/polynomials.rs index 2c9c60c8..b4d4144e 100644 --- a/element/src/polynomials.rs +++ b/element/src/polynomials.rs @@ -1,11 +1,12 @@ //! Orthonormal polynomials -use bempp_traits::arrays::{Array2DAccess, Array3DAccess}; +use bempp_traits::arrays::Array3DAccess; use bempp_traits::cell::ReferenceCellType; +use rlst_dense::{RandomAccessByRef, Shape}; /// Tabulate orthonormal polynomials on a interval -fn tabulate_legendre_polynomials_interval<'a>( - points: &impl Array2DAccess<'a, f64>, +fn tabulate_legendre_polynomials_interval + Shape>( + points: &T, degree: usize, derivatives: usize, data: &mut impl Array3DAccess, @@ -57,8 +58,8 @@ fn quad_index(i: usize, j: usize, n: usize) -> usize { } /// Tabulate orthonormal polynomials on a quadrilateral -fn tabulate_legendre_polynomials_quadrilateral<'a>( - points: &impl Array2DAccess<'a, f64>, +fn tabulate_legendre_polynomials_quadrilateral + Shape>( + points: &T, degree: usize, derivatives: usize, data: &mut impl Array3DAccess, @@ -191,8 +192,8 @@ fn tabulate_legendre_polynomials_quadrilateral<'a>( } } /// Tabulate orthonormal polynomials on a triangle -fn tabulate_legendre_polynomials_triangle<'a>( - points: &impl Array2DAccess<'a, f64>, +fn tabulate_legendre_polynomials_triangle + Shape>( + points: &T, degree: usize, derivatives: usize, data: &mut impl Array3DAccess, @@ -368,9 +369,9 @@ pub fn derivative_count(cell_type: ReferenceCellType, derivatives: usize) -> usi } } -pub fn legendre_shape<'a>( +pub fn legendre_shape + Shape>( cell_type: ReferenceCellType, - points: &impl Array2DAccess<'a, f64>, + points: &T, degree: usize, derivatives: usize, ) -> (usize, usize, usize) { @@ -382,9 +383,9 @@ pub fn legendre_shape<'a>( } /// Tabulate orthonormal polynomials -pub fn tabulate_legendre_polynomials<'a>( +pub fn tabulate_legendre_polynomials + Shape>( cell_type: ReferenceCellType, - points: &impl Array2DAccess<'a, f64>, + points: &T, degree: usize, derivatives: usize, data: &mut impl Array3DAccess, @@ -410,14 +411,13 @@ mod test { use crate::polynomials::*; use approx::*; use bempp_quadrature::simplex_rules::simplex_rule; - use bempp_tools::arrays::{Array2D, Array3D}; - + use bempp_tools::arrays::{to_matrix, Array3D}; #[test] fn test_legendre_interval() { let degree = 6; let rule = simplex_rule(ReferenceCellType::Interval, degree + 1).unwrap(); - let points = Array2D::from_data(rule.points, (rule.npoints, 1)); + let points = to_matrix(&rule.points, (rule.npoints, 1)); let mut data = Array3D::::new(legendre_shape( ReferenceCellType::Interval, @@ -448,7 +448,7 @@ mod test { let degree = 5; let rule = simplex_rule(ReferenceCellType::Triangle, 79).unwrap(); - let points = Array2D::from_data(rule.points, (rule.npoints, 2)); + let points = to_matrix(&rule.points, (rule.npoints, 2)); let mut data = Array3D::::new(legendre_shape( ReferenceCellType::Triangle, @@ -479,7 +479,7 @@ mod test { let degree = 5; let rule = simplex_rule(ReferenceCellType::Quadrilateral, 85).unwrap(); - let points = Array2D::from_data(rule.points, (rule.npoints, 2)); + let points = to_matrix(&rule.points, (rule.npoints, 2)); let mut data = Array3D::::new(legendre_shape( ReferenceCellType::Quadrilateral, @@ -521,7 +521,7 @@ mod test { p[2 * i] = i as f64 / 10.0; p[2 * i + 1] = p[2 * i] + epsilon; } - let points = Array2D::from_data(p, (20, 1)); + let points = to_matrix(&p, (20, 1)); let mut data = Array3D::::new(legendre_shape( ReferenceCellType::Interval, @@ -560,7 +560,7 @@ mod test { index += 1; } } - let points = Array2D::from_data(p, (165, 2)); + let points = to_matrix(&p, (165, 2)); let mut data = Array3D::::new(legendre_shape( ReferenceCellType::Triangle, @@ -603,7 +603,7 @@ mod test { p[6 * index + 5] = p[6 * index + 1] + epsilon; } } - let points = Array2D::from_data(p, (300, 2)); + let points = to_matrix(&p, (300, 2)); let mut data = Array3D::::new(legendre_shape( ReferenceCellType::Quadrilateral, @@ -643,7 +643,7 @@ mod test { for (i, pi) in p.iter_mut().enumerate() { *pi = i as f64 / 10.0; } - let points = Array2D::from_data(p, (11, 1)); + let points = to_matrix(&p, (11, 1)); let mut data = Array3D::::new(legendre_shape( ReferenceCellType::Interval, @@ -729,7 +729,7 @@ mod test { p[2 * (11 * i + j) + 1] = j as f64 / 10.0; } } - let points = Array2D::from_data(p, (121, 2)); + let points = to_matrix(&p, (121, 2)); let mut data = Array3D::::new(legendre_shape( ReferenceCellType::Quadrilateral, diff --git a/grid/src/grid.rs b/grid/src/grid.rs index dbe2a089..1f193a9a 100644 --- a/grid/src/grid.rs +++ b/grid/src/grid.rs @@ -112,6 +112,40 @@ impl Geometry for SerialGeometry { fn index_map(&self) -> &[usize] { &self.index_map } + + fn get_compute_points_function + Shape>( + &self, + element: &impl FiniteElement, + points: &T, + ) -> Box { + let npts = points.shape().0; + let mut table = Array4D::::new(element.tabulate_array_shape(0, npts)); + element.tabulate_rlst(points, 0, &mut table); + let gdim = self.dim(); + + Box::new(|cell: usize, pts: &mut T| { + for p in 0..npts { + for i in 0..gdim { + unsafe { + *pts.get_unchecked_mut(p, i) = 0.0; + } + } + } + let vertices = self.cell_vertices(cell).unwrap(); + for (i, n) in vertices.iter().enumerate() { + let pt = self.point(*n).unwrap(); + for p in 0..points.shape().0 { + for (j, pt_j) in pt.iter().enumerate() { + unsafe { + *pts.get_unchecked_mut(p, j) += + *pt_j * *table.get_unchecked(0, p, i, 0); + } + } + } + } + }) + } + fn compute_points<'a>( &self, points: &impl Array2DAccess<'a, f64>, diff --git a/tools/Cargo.toml b/tools/Cargo.toml index 0f877963..ebc5b177 100644 --- a/tools/Cargo.toml +++ b/tools/Cargo.toml @@ -24,3 +24,4 @@ libc = "0.2" num = "0.4" bempp-traits = { path = "../traits"} rayon = "1.7" +rlst-dense = { git = "https://github.com/linalg-rs/rlst.git" } diff --git a/tools/examples/array2d.rs b/tools/examples/array2d.rs deleted file mode 100644 index d594c9e0..00000000 --- a/tools/examples/array2d.rs +++ /dev/null @@ -1,52 +0,0 @@ -use bempp_tools::arrays::Array2D; -use bempp_traits::arrays::Array2DAccess; - -fn main() { - // Create an 2D array of zero integers of a given size - // A 2D array is a two-dimensional storage format that allows items to be accessed using two indices (a row and a column) - let mut a = Array2D::::new((5, 3)); - - // Set values in the first row of the array - *a.get_mut(0, 0).unwrap() = 3; - *a.get_mut(0, 1).unwrap() = 1; - *a.get_mut(0, 2).unwrap() = 4; - - // Print values from the first column of the array - println!( - "The first column of a is {} {} {} {} {}", - *a.get(0, 0).unwrap(), - *a.get(1, 0).unwrap(), - *a.get(2, 0).unwrap(), - *a.get(3, 0).unwrap(), - *a.get(4, 0).unwrap() - ); - - // Items can also be accessed without bound checking using the unchecked methods - println!( - "The entry in the first row and first column of a is {}", - unsafe { *a.get_unchecked(0, 0) } - ); - - // Print values from the first row of the array - // This can be done in two ways - println!( - "The first row of a is {} {} {}", - *a.get(0, 0).unwrap(), - *a.get(0, 1).unwrap(), - *a.get(0, 2).unwrap() - ); - let row0 = a.row(0).unwrap(); - println!("The first row of a is {} {} {}", row0[0], row0[1], row0[2]); - - // Print the shape of the array - println!("a has {} rows and {} columns", a.shape().0, a.shape().1); - - // Iterate through the rows of a - for (i, row) in a.iter_rows().enumerate() { - println!( - "The sum of the values in row {} of a is {}", - i, - row.iter().sum::() - ); - } -} diff --git a/tools/src/arrays.rs b/tools/src/arrays.rs index 008e12ef..b6ed547e 100644 --- a/tools/src/arrays.rs +++ b/tools/src/arrays.rs @@ -1,95 +1,23 @@ //! Containers to store multi-dimensional data -use bempp_traits::arrays::{AdjacencyListAccess, Array2DAccess, Array3DAccess, Array4DAccess}; +use bempp_traits::arrays::{AdjacencyListAccess, Array3DAccess, Array4DAccess}; use num::Num; +use rlst_dense::{operations::transpose::Scalar, rlst_dynamic_mat, UnsafeRandomAccessMut}; use std::clone::Clone; -/// A two-dimensional rectangular array -#[derive(Clone)] -pub struct Array2D { - /// The data in the array, in row-major order - pub data: Vec, - /// The shape of the array - shape: (usize, usize), -} - -impl Array2D { - /// Create an array from a data vector - pub fn new(shape: (usize, usize)) -> Self { - Self { - data: vec![T::zero(); shape.0 * shape.1], - shape, +pub type Mat = rlst_dense::Matrix< + T, + rlst_dense::base_matrix::BaseMatrix, rlst_dense::Dynamic>, + rlst_dense::Dynamic, +>; + +pub fn to_matrix(data: &[T], shape: (usize, usize)) -> Mat { + let mut mat = rlst_dynamic_mat![T, shape]; + for i in 0..shape.0 * shape.1 { + unsafe { + *mat.get_unchecked_mut(i / shape.1, i % shape.1) = data[i]; } } - - /// Create an array from a data vector - pub fn from_data(data: Vec, shape: (usize, usize)) -> Self { - assert_eq!(data.len(), shape.0 * shape.1); - Self { data, shape } - } -} - -impl<'a, T: Num + 'a> Array2DAccess<'a, T> for Array2D { - type I = Array2DRowIterator<'a, T>; - - fn get(&self, index0: usize, index1: usize) -> Option<&T> { - if index0 >= self.shape.0 || index1 >= self.shape.1 { - None - } else { - unsafe { Some(self.get_unchecked(index0, index1)) } - } - } - /// Get a mutable item from the array - fn get_mut(&mut self, index0: usize, index1: usize) -> Option<&mut T> { - if index0 >= self.shape.0 || index1 >= self.shape.1 { - None - } else { - unsafe { Some(self.get_unchecked_mut(index0, index1)) } - } - } - /// Get a row of the array - fn row(&self, index: usize) -> Option<&[T]> { - if index >= self.shape.0 { - None - } else { - unsafe { Some(self.row_unchecked(index)) } - } - } - /// Get an item from the array without checking bounds - unsafe fn get_unchecked(&self, index0: usize, index1: usize) -> &T { - self.data.get_unchecked(index0 * self.shape.1 + index1) - } - /// Get a mutable item from the array without checking bounds - unsafe fn get_unchecked_mut(&mut self, index0: usize, index1: usize) -> &mut T { - self.data.get_unchecked_mut(index0 * self.shape.1 + index1) - } - /// Get a row of the array without checking bounds - unsafe fn row_unchecked(&self, index: usize) -> &[T] { - &self.data[index * self.shape.1..(index + 1) * self.shape.1] - } - /// Get the shape of the array - fn shape(&self) -> &(usize, usize) { - &self.shape - } - /// Iterate through the rows - fn iter_rows(&'a self) -> Self::I { - Array2DRowIterator:: { - array: self, - index: 0, - } - } -} - -pub struct Array2DRowIterator<'a, T: Num> { - array: &'a Array2D, - index: usize, -} - -impl<'a, T: Num> Iterator for Array2DRowIterator<'a, T> { - type Item = &'a [T]; - fn next(&mut self) -> Option { - self.index += 1; - self.array.row(self.index - 1) - } + mat } /// A three-dimensional rectangular array @@ -338,57 +266,6 @@ impl<'a, T: Num> Iterator for AdjacencyListRowIterator<'a, T> { mod test { use crate::arrays::*; - #[test] - fn test_array_2d() { - let mut arr = Array2D::from_data(vec![1, 2, 3, 4, 5, 6], (2, 3)); - assert_eq!(*arr.get(0, 0).unwrap(), 1); - assert_eq!(*arr.get(0, 1).unwrap(), 2); - assert_eq!(*arr.get(0, 2).unwrap(), 3); - assert_eq!(*arr.get(1, 0).unwrap(), 4); - assert_eq!(*arr.get(1, 1).unwrap(), 5); - assert_eq!(*arr.get(1, 2).unwrap(), 6); - - let row1 = arr.row(1).unwrap(); - assert_eq!(row1.len(), 3); - assert_eq!(row1[0], 4); - assert_eq!(row1[1], 5); - assert_eq!(row1[2], 6); - - *arr.get_mut(1, 2).unwrap() = 7; - assert_eq!(*arr.get(1, 2).unwrap(), 7); - - for (index, row) in arr.iter_rows().enumerate() { - assert_eq!(*arr.get(index, 0).unwrap(), row[0]); - } - - let mut arr2 = Array2D::from_data(vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0], (2, 3)); - assert_eq!(*arr2.get(0, 0).unwrap(), 1.0); - assert_eq!(*arr2.get(0, 1).unwrap(), 2.0); - assert_eq!(*arr2.get(0, 2).unwrap(), 3.0); - assert_eq!(*arr2.get(1, 0).unwrap(), 4.0); - assert_eq!(*arr2.get(1, 1).unwrap(), 5.0); - assert_eq!(*arr2.get(1, 2).unwrap(), 6.0); - - let row1 = arr2.row(1).unwrap(); - assert_eq!(row1.len(), 3); - assert_eq!(row1[0], 4.0); - assert_eq!(row1[1], 5.0); - assert_eq!(row1[2], 6.0); - - *arr2.get_mut(1, 2).unwrap() = 7.; - assert_eq!(*arr2.get(1, 2).unwrap(), 7.0); - - let mut arr3 = Array2D::::new((4, 5)); - assert_eq!(*arr3.get(1, 2).unwrap(), 0); - *arr3.get_mut(1, 2).unwrap() = 5; - assert_eq!(*arr3.get(1, 2).unwrap(), 5); - - let mut arr4 = Array2D::::new((4, 5)); - assert_eq!(*arr4.get(1, 2).unwrap(), 0.0); - *arr4.get_mut(1, 2).unwrap() = 5.0; - assert_eq!(*arr4.get(1, 2).unwrap(), 5.0); - } - #[test] fn test_adjacency_list() { let mut arr = AdjacencyList::from_data(vec![1, 2, 3, 4, 5, 6], vec![0, 2, 3, 6]); diff --git a/traits/src/arrays.rs b/traits/src/arrays.rs index a9625cbb..4af2a01c 100644 --- a/traits/src/arrays.rs +++ b/traits/src/arrays.rs @@ -1,71 +1,6 @@ //! Containers to store multi-dimensional data use num::Num; -pub trait Array1DAccess<'a, T: Num> { - type I: Iterator; - - /// Get an item from the array - fn get(&self, index: usize) -> Option<&T>; - - /// Get a mutable item from the array - fn get_mut(&mut self, index: usize) -> Option<&mut T>; - - /// Get an item from the array without checking bounds - /// - /// # Safety - /// This function does not perform bound checks - unsafe fn get_unchecked(&self, index: usize) -> &T; - - /// Get a mutable item from the array without checking bounds - /// - /// # Safety - /// This function does not perform bound checks - unsafe fn get_unchecked_mut(&mut self, index: usize) -> &mut T; - - /// Get the shape of the array - fn shape(&self) -> &usize; - - /// Iterate through the values - fn iter(&'a self) -> Self::I; -} - -pub trait Array2DAccess<'a, T: Num> { - type I: Iterator; - - /// Get an item from the array - fn get(&self, index0: usize, index1: usize) -> Option<&T>; - - /// Get a mutable item from the array - fn get_mut(&mut self, index0: usize, index1: usize) -> Option<&mut T>; - - /// Get a row of the array - fn row(&self, index: usize) -> Option<&[T]>; - - /// Get an item from the array without checking bounds - /// - /// # Safety - /// This function does not perform bound checks - unsafe fn get_unchecked(&self, index0: usize, index1: usize) -> &T; - - /// Get a mutable item from the array without checking bounds - /// - /// # Safety - /// This function does not perform bound checks - unsafe fn get_unchecked_mut(&mut self, index0: usize, index1: usize) -> &mut T; - - /// Get a row of the array without checking bounds - /// - /// # Safety - /// This function does not perform bound checks - unsafe fn row_unchecked(&self, index: usize) -> &[T]; - - /// Get the shape of the array - fn shape(&self) -> &(usize, usize); - - /// Iterate through the rows - fn iter_rows(&'a self) -> Self::I; -} - pub trait Array3DAccess { /// Get an item from the array fn get(&self, index0: usize, index1: usize, index2: usize) -> Option<&T>; diff --git a/traits/src/element.rs b/traits/src/element.rs index 24c27fd9..91bd4600 100644 --- a/traits/src/element.rs +++ b/traits/src/element.rs @@ -1,7 +1,8 @@ //! Finite element definitions -use crate::arrays::{Array2DAccess, Array4DAccess}; +use crate::arrays::Array4DAccess; use crate::cell::ReferenceCellType; +use rlst_common::traits::{RandomAccessByRef, Shape}; /// The family of an element #[derive(Debug, PartialEq, Eq, Clone, Copy, Hash)] @@ -70,9 +71,9 @@ pub trait FiniteElement { fn value_size(&self) -> usize; /// Tabulate the values of the basis functions and their derivatives at a set of points - fn tabulate<'a>( + fn tabulate + Shape>( &self, - points: &impl Array2DAccess<'a, f64>, + points: &T, nderivs: usize, data: &mut impl Array4DAccess, ); diff --git a/traits/src/grid.rs b/traits/src/grid.rs index e01587f1..4a27cfa6 100644 --- a/traits/src/grid.rs +++ b/traits/src/grid.rs @@ -1,8 +1,9 @@ //! Geometry and topology definitions -use crate::arrays::{AdjacencyListAccess, Array2DAccess}; +use crate::arrays::AdjacencyListAccess; use crate::cell::ReferenceCellType; -use rlst_common::traits::{RandomAccessMut, Shape}; +use crate::element::FiniteElement; +use rlst_common::traits::{RandomAccessByRef, RandomAccessMut, Shape}; /// The ownership of a mesh entity #[derive(Debug, PartialEq, Eq, Clone, Copy, Hash)] @@ -35,63 +36,57 @@ pub trait Geometry { fn index_map(&self) -> &[usize]; /// Compute the physical coordinates of a set of points in a given cell - fn compute_points<'a>( + fn get_compute_points_function< + T: RandomAccessByRef + Shape, + TMut: RandomAccessMut + Shape, + >( &self, - points: &impl Array2DAccess<'a, f64>, - cell: usize, - physical_points: &mut impl Array2DAccess<'a, f64>, - ); - - /// Compute the physical coordinates of a set of points in a given cell, - /// with return in transposed arrangement - fn compute_points_transpose<'a>( - &self, - points: &impl Array2DAccess<'a, f64>, - cell: usize, - physical_points: &mut impl Array2DAccess<'a, f64>, - ); - - /// Compute the physical coordinates of a set of points in a given cell, - /// with return in transposed arrangement - fn compute_points_transpose_rlst<'a, T: RandomAccessMut + Shape>( - &self, - points: &impl Array2DAccess<'a, f64>, - cell: usize, - physical_points: &mut T, - ); + element: &impl FiniteElement, + points: &T, + ) -> Box; /// Compute the physical coordinates of a set of points in a given cell - fn compute_points_rlst<'a, T: RandomAccessMut + Shape>( + fn compute_points< + 'a, + T: RandomAccessByRef + Shape, + TMut: RandomAccessMut + Shape, + >( &self, - points: &impl Array2DAccess<'a, f64>, + points: &T, cell: usize, - physical_points: &mut T, + physical_points: &mut TMut, ); /// Compute the normals to a set of points in a given cell - fn compute_normals<'a>( + fn compute_normals< + T: RandomAccessByRef + Shape, + TMut: RandomAccessMut + Shape, + >( &self, - points: &impl Array2DAccess<'a, f64>, + points: &T, cell: usize, - normals: &mut impl Array2DAccess<'a, f64>, + normals: &mut TMut, ); /// Evaluate the jacobian at a set of points in a given cell /// /// The input points should be given using coordinates on the reference element - fn compute_jacobians<'a>( + fn compute_jacobians< + T: RandomAccessByRef + Shape, + TMut: RandomAccessMut + Shape, + >( &self, - points: &impl Array2DAccess<'a, f64>, + points: &T, cell: usize, - jacobians: &mut impl Array2DAccess<'a, f64>, + jacobians: &mut TMut, ); /// Evaluate the determinand of the jacobian at a set of points in a given cell /// /// The input points should be given using coordinates on the reference element - fn compute_jacobian_determinants<'a>( + fn compute_jacobian_determinants + Shape>( &self, - points: &impl Array2DAccess<'a, f64>, + points: &T, cell: usize, jacobian_determinants: &mut [f64], ); @@ -99,11 +94,14 @@ pub trait Geometry { /// Evaluate the jacobian inverse at a set of points in a given cell /// /// The input points should be given using coordinates on the reference element - fn compute_jacobian_inverses<'a>( + fn compute_jacobian_inverses< + T: RandomAccessByRef + Shape, + TMut: RandomAccessMut + Shape, + >( &self, - points: &impl Array2DAccess<'a, f64>, + points: &T, cell: usize, - jacobian_inverses: &mut impl Array2DAccess<'a, f64>, + jacobian_inverses: &mut TMut, ); }