From 592b8103971289bf06db4fcacef05ffa0b8c4e4c Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Tue, 16 Jan 2024 11:46:28 +0000 Subject: [PATCH] Use rlst directly rather than helper functions from tools/arrays (#158) * remove rlst objects and functions from tools/array * parallel grid * types * fmt * update parallel grid example * cargo fmt --- bem/examples/assemble.rs | 8 +- bem/src/assembly.rs | 19 +++-- bem/src/assembly/batched.rs | 78 +++++++++++------ element/src/element.rs | 26 ++++-- element/src/element/lagrange.rs | 118 +++++++++++++++++++------- element/src/element/raviart_thomas.rs | 28 +++--- element/src/polynomials.rs | 34 +++++--- field/src/array.rs | 12 ++- field/src/field.rs | 9 +- grid/examples/curved_cells.rs | 18 ++-- grid/examples/parallel_grid.rs | 6 +- grid/src/grid.rs | 106 +++++++++++++++-------- grid/src/io.rs | 15 +++- grid/src/parallel_grid.rs | 20 +++-- grid/src/shapes.rs | 29 ++++--- kernel/src/laplace_3d.rs | 20 +++-- tools/Cargo.toml | 2 - tools/src/arrays.rs | 33 ------- 18 files changed, 363 insertions(+), 218 deletions(-) diff --git a/bem/examples/assemble.rs b/bem/examples/assemble.rs index 3af85f69..1d23f9ea 100644 --- a/bem/examples/assemble.rs +++ b/bem/examples/assemble.rs @@ -2,11 +2,11 @@ use bempp_bem::assembly::{assemble_batched, BoundaryOperator, PDEType}; use bempp_bem::function_space::SerialFunctionSpace; use bempp_element::element::create_element; use bempp_grid::shapes::regular_sphere; -use bempp_tools::arrays::zero_matrix; use bempp_traits::bem::DofMap; use bempp_traits::bem::FunctionSpace; use bempp_traits::cell::ReferenceCellType; use bempp_traits::element::{Continuity, ElementFamily}; +use rlst_dense::rlst_dynamic_array2; fn main() { println!("Creating grid"); @@ -33,8 +33,10 @@ fn main() { space1.dofmap().global_size(), space0.dofmap().global_size() ); - let mut matrix = - zero_matrix::([space1.dofmap().global_size(), space0.dofmap().global_size()]); + let mut matrix = rlst_dynamic_array2!( + f64, + [space1.dofmap().global_size(), space0.dofmap().global_size()] + ); println!("Assembling dense matrix (complex)"); assemble_batched( diff --git a/bem/src/assembly.rs b/bem/src/assembly.rs index a5e3bae5..bc26bb74 100644 --- a/bem/src/assembly.rs +++ b/bem/src/assembly.rs @@ -1,7 +1,7 @@ pub mod batched; use crate::function_space::SerialFunctionSpace; use bempp_kernel::laplace_3d; -use bempp_tools::arrays::Mat; +use rlst_dense::{array::Array, base_array::BaseArray, data_container::VectorContainer}; #[derive(Debug, PartialEq, Eq, Clone, Copy, Hash)] #[repr(u8)] @@ -24,7 +24,7 @@ pub enum PDEType { /// Assemble an operator into a dense matrix using batched parallelisation pub fn assemble_batched<'a>( // TODO: ouput should be `&mut impl ArrayAccess2D` once such a trait exists - output: &mut Mat, + output: &mut Array, 2>, 2>, operator: BoundaryOperator, pde: PDEType, trial_space: &SerialFunctionSpace<'a>, @@ -61,13 +61,12 @@ mod test { use bempp_element::element::create_element; use bempp_grid::shapes::regular_sphere; use bempp_kernel::laplace_3d::Laplace3dKernel; - use bempp_tools::arrays::zero_matrix; use bempp_traits::bem::DofMap; use bempp_traits::cell::ReferenceCellType; use bempp_traits::element::{Continuity, ElementFamily}; // use num::complex::Complex; use bempp_traits::bem::FunctionSpace; - use rlst_dense::traits::RandomAccessByRef; + use rlst_dense::{rlst_dynamic_array2, traits::RandomAccessByRef}; #[test] fn test_laplace_single_layer() { @@ -87,8 +86,10 @@ mod test { let space0 = SerialFunctionSpace::new(&grid, &element0); let space1 = SerialFunctionSpace::new(&grid, &element1); - let mut matrix = - zero_matrix::([space1.dofmap().global_size(), space0.dofmap().global_size()]); + let mut matrix = rlst_dynamic_array2!( + f64, + [space1.dofmap().global_size(), space0.dofmap().global_size()] + ); batched::assemble( &mut matrix, &Laplace3dKernel::new(), @@ -98,8 +99,10 @@ mod test { &space1, ); - let mut matrix2 = - zero_matrix::([space1.dofmap().global_size(), space0.dofmap().global_size()]); + let mut matrix2 = rlst_dynamic_array2!( + f64, + [space1.dofmap().global_size(), space0.dofmap().global_size()] + ); assemble_batched( &mut matrix2, diff --git a/bem/src/assembly/batched.rs b/bem/src/assembly/batched.rs index f584eb84..f0a4e47e 100644 --- a/bem/src/assembly/batched.rs +++ b/bem/src/assembly/batched.rs @@ -3,7 +3,6 @@ use bempp_quadrature::duffy::quadrilateral::quadrilateral_duffy; use bempp_quadrature::duffy::triangle::triangle_duffy; use bempp_quadrature::simplex_rules::simplex_rule; use bempp_quadrature::types::{CellToCellConnectivity, TestTrialNumericalQuadratureDefinition}; -use bempp_tools::arrays::{transpose_to_matrix, zero_matrix, Array4D, Mat}; use bempp_traits::arrays::AdjacencyListAccess; use bempp_traits::bem::{DofMap, FunctionSpace}; use bempp_traits::cell::ReferenceCellType; @@ -13,9 +12,14 @@ use bempp_traits::kernel::Kernel; use bempp_traits::types::EvalType; use bempp_traits::types::Scalar; use rayon::prelude::*; -use rlst_dense::rlst_dynamic_array4; -use rlst_dense::traits::{ - RandomAccessByRef, RawAccess, RawAccessMut, Shape, UnsafeRandomAccessByRef, +use rlst_dense::{ + array::Array, + base_array::BaseArray, + data_container::VectorContainer, + rlst_dynamic_array2, rlst_dynamic_array4, + traits::{ + RandomAccessByRef, RandomAccessMut, RawAccess, RawAccessMut, Shape, UnsafeRandomAccessByRef, + }, }; fn get_quadrature_rule( @@ -88,11 +92,11 @@ fn assemble_batch_singular<'a>( trial_space: &SerialFunctionSpace<'a>, test_space: &SerialFunctionSpace<'a>, cell_pairs: &[(usize, usize)], - trial_points: &Mat, - test_points: &Mat, + trial_points: &Array, 2>, 2>, + test_points: &Array, 2>, 2>, weights: &[f64], - trial_table: &Array4D, - test_table: &Array4D, + trial_table: &Array, 4>, 4>, + test_table: &Array, 4>, 4>, ) -> usize { let grid = test_space.grid(); let mut k = vec![0.0]; @@ -100,10 +104,10 @@ fn assemble_batch_singular<'a>( // Memory assignment to be moved elsewhere as passed into here mutable? let mut test_jdet = vec![0.0; test_points.shape()[0]]; let mut trial_jdet = vec![0.0; trial_points.shape()[0]]; - let mut test_mapped_pts = zero_matrix([test_points.shape()[0], 3]); - let mut trial_mapped_pts = zero_matrix([trial_points.shape()[0], 3]); - let mut test_normals = zero_matrix([test_points.shape()[0], 3]); - let mut trial_normals = zero_matrix([trial_points.shape()[0], 3]); + let mut test_mapped_pts = rlst_dynamic_array2!(f64, [test_points.shape()[0], 3]); + let mut trial_mapped_pts = rlst_dynamic_array2!(f64, [trial_points.shape()[0], 3]); + let mut test_normals = rlst_dynamic_array2!(f64, [test_points.shape()[0], 3]); + let mut trial_normals = rlst_dynamic_array2!(f64, [trial_points.shape()[0], 3]); for (test_cell, trial_cell) in cell_pairs { let test_cell_tindex = grid.topology().index_map()[*test_cell]; @@ -186,12 +190,12 @@ fn assemble_batch_nonadjacent<'a, const NPTS_TEST: usize, const NPTS_TRIAL: usiz trial_cells: &[usize], test_space: &SerialFunctionSpace<'a>, test_cells: &[usize], - trial_points: &Mat, + trial_points: &Array, 2>, 2>, trial_weights: &[f64], - test_points: &Mat, + test_points: &Array, 2>, 2>, test_weights: &[f64], - trial_table: &Array4D, - test_table: &Array4D, + trial_table: &Array, 4>, 4>, + test_table: &Array, 4>, 4>, ) -> usize { debug_assert!(test_weights.len() == NPTS_TEST); debug_assert!(test_points.shape()[0] == NPTS_TEST); @@ -205,7 +209,7 @@ fn assemble_batch_nonadjacent<'a, const NPTS_TEST: usize, const NPTS_TRIAL: usiz let mut k = vec![0.0; NPTS_TEST * NPTS_TRIAL]; let mut test_jdet = [0.0; NPTS_TEST]; - let mut test_normals = zero_matrix([NPTS_TEST, 3]); + let mut test_normals = rlst_dynamic_array2!(f64, [NPTS_TEST, 3]); let mut test_mapped_pts = rlst_dense::rlst_dynamic_array2![f64, [NPTS_TEST, 3]]; @@ -223,8 +227,8 @@ fn assemble_batch_nonadjacent<'a, const NPTS_TEST: usize, const NPTS_TRIAL: usiz let mut trial_mapped_pts = vec![]; let mut trial_normals = vec![]; for _i in 0..trial_cells.len() { - trial_mapped_pts.push(zero_matrix([NPTS_TRIAL, 3])); - trial_normals.push(zero_matrix([NPTS_TRIAL, 3])); + trial_mapped_pts.push(rlst_dynamic_array2!(f64, [NPTS_TRIAL, 3])); + trial_normals.push(rlst_dynamic_array2!(f64, [NPTS_TRIAL, 3])); } for (trial_cell_i, trial_cell) in trial_cells.iter().enumerate() { @@ -322,7 +326,7 @@ fn assemble_batch_nonadjacent<'a, const NPTS_TEST: usize, const NPTS_TRIAL: usiz } pub fn assemble<'a>( - output: &mut Mat, + output: &mut Array, 2>, 2>, kernel: &impl Kernel, needs_trial_normal: bool, needs_test_normal: bool, @@ -358,7 +362,7 @@ pub fn assemble<'a>( #[allow(clippy::too_many_arguments)] pub fn assemble_nonsingular<'a, const NPTS_TEST: usize, const NPTS_TRIAL: usize>( - output: &mut Mat, + output: &mut Array, 2>, 2>, kernel: &impl Kernel, trial_space: &SerialFunctionSpace<'a>, test_space: &SerialFunctionSpace<'a>, @@ -381,10 +385,20 @@ pub fn assemble_nonsingular<'a, const NPTS_TEST: usize, const NPTS_TRIAL: usize> // TODO: pass cell types into this function let qrule_test = simplex_rule(ReferenceCellType::Triangle, NPTS_TEST).unwrap(); - let qpoints_test = transpose_to_matrix(&qrule_test.points, [NPTS_TEST, 2]); + let mut qpoints_test = rlst_dynamic_array2!(f64, [NPTS_TEST, 2]); + for i in 0..NPTS_TEST { + for j in 0..2 { + *qpoints_test.get_mut([i, j]).unwrap() = qrule_test.points[2 * i + j]; + } + } let qweights_test = qrule_test.weights; let qrule_trial = simplex_rule(ReferenceCellType::Triangle, NPTS_TRIAL).unwrap(); - let qpoints_trial = transpose_to_matrix(&qrule_trial.points, [NPTS_TRIAL, 2]); + let mut qpoints_trial = rlst_dynamic_array2!(f64, [NPTS_TRIAL, 2]); + for i in 0..NPTS_TRIAL { + for j in 0..2 { + *qpoints_trial.get_mut([i, j]).unwrap() = qrule_trial.points[2 * i + j]; + } + } let qweights_trial = qrule_trial.weights; let mut test_table = @@ -460,7 +474,7 @@ pub fn assemble_nonsingular<'a, const NPTS_TEST: usize, const NPTS_TRIAL: usize> #[allow(clippy::too_many_arguments)] pub fn assemble_singular<'a>( - output: &mut Mat, + output: &mut Array, 2>, 2>, kernel: &impl Kernel, needs_trial_normal: bool, needs_test_normal: bool, @@ -534,7 +548,12 @@ pub fn assemble_singular<'a>( npoints, ); - let points = transpose_to_matrix(&qrule.trial_points, [qrule.npoints, 2]); + let mut points = rlst_dynamic_array2!(f64, [qrule.npoints, 2]); + for i in 0..qrule.npoints { + for j in 0..2 { + *points.get_mut([i, j]).unwrap() = qrule.trial_points[2 * i + j]; + } + } let mut table = rlst_dynamic_array4!( f64, trial_space @@ -545,7 +564,12 @@ pub fn assemble_singular<'a>( trial_points.push(points); trial_tables.push(table); - let points = transpose_to_matrix(&qrule.test_points, [qrule.npoints, 2]); + let mut points = rlst_dynamic_array2!(f64, [qrule.npoints, 2]); + for i in 0..qrule.npoints { + for j in 0..2 { + *points.get_mut([i, j]).unwrap() = qrule.test_points[2 * i + j]; + } + } let mut table = rlst_dynamic_array4!( f64, test_space @@ -644,7 +668,7 @@ mod test { let ndofs = space.dofmap().global_size(); - let mut matrix = zero_matrix::([ndofs, ndofs]); + let mut matrix = rlst_dynamic_array2!(f64, [ndofs, ndofs]); assemble( &mut matrix, &Laplace3dKernel::new(), diff --git a/element/src/element.rs b/element/src/element.rs index 6c83f5c5..b841a2ab 100644 --- a/element/src/element.rs +++ b/element/src/element.rs @@ -2,16 +2,24 @@ use crate::cell::create_cell; use crate::polynomials::{legendre_shape, polynomial_count, tabulate_legendre_polynomials}; -use bempp_tools::arrays::{AdjacencyList, Array3D, Mat}; +use bempp_tools::arrays::AdjacencyList; use bempp_traits::arrays::AdjacencyListAccess; use bempp_traits::cell::ReferenceCellType; use bempp_traits::element::{Continuity, ElementFamily, FiniteElement, MapType}; use rlst_dense::linalg::inverse::MatrixInverse; -use rlst_dense::traits::{RandomAccessByRef, RandomAccessMut, Shape, UnsafeRandomAccessMut}; +use rlst_dense::{ + array::Array, + base_array::BaseArray, + data_container::VectorContainer, + traits::{RandomAccessByRef, RandomAccessMut, Shape, UnsafeRandomAccessMut}, +}; use rlst_dense::{rlst_dynamic_array2, rlst_dynamic_array3}; pub mod lagrange; pub mod raviart_thomas; +type EntityPoints = [Vec, 2>, 2>>; 4]; +type EntityWeights = [Vec, 3>, 3>>; 4]; + pub struct CiarletElement { cell_type: ReferenceCellType, degree: usize, @@ -22,10 +30,10 @@ pub struct CiarletElement { family: ElementFamily, continuity: Continuity, dim: usize, - coefficients: Array3D, + coefficients: Array, 3>, 3>, entity_dofs: [AdjacencyList; 4], - // interpolation_points: [Vec>; 4], - // interpolation_weights: [Vec>; 4], + // interpolation_points: EntityPoints, + // interpolation_weights: EntityWeights, } impl CiarletElement { @@ -36,9 +44,9 @@ impl CiarletElement { cell_type: ReferenceCellType, degree: usize, value_shape: Vec, - polynomial_coeffs: Array3D, - interpolation_points: [Vec>; 4], - interpolation_weights: [Vec>; 4], + polynomial_coeffs: Array, 3>, 3>, + interpolation_points: EntityPoints, + interpolation_weights: EntityWeights, map_type: MapType, continuity: Continuity, highest_degree: usize, @@ -69,7 +77,7 @@ impl CiarletElement { } let new_pts = if continuity == Continuity::Discontinuous { - let mut new_pts: [Vec>; 4] = [vec![], vec![], vec![], vec![]]; + let mut new_pts: EntityPoints = [vec![], vec![], vec![], vec![]]; let mut pn = 0; let mut all_pts = rlst_dynamic_array2![f64, [npts, tdim]]; for (i, pts_i) in interpolation_points.iter().take(tdim).enumerate() { diff --git a/element/src/element/lagrange.rs b/element/src/element/lagrange.rs index 3e272307..94248b23 100644 --- a/element/src/element/lagrange.rs +++ b/element/src/element/lagrange.rs @@ -2,11 +2,9 @@ use crate::element::{create_cell, CiarletElement}; use crate::polynomials::polynomial_count; -use bempp_tools::arrays::{to_matrix, zero_matrix}; use bempp_traits::cell::ReferenceCellType; use bempp_traits::element::{Continuity, ElementFamily, MapType}; -use rlst_dense::rlst_dynamic_array3; -use rlst_dense::traits::RandomAccessMut; +use rlst_dense::{rlst_dynamic_array2, rlst_dynamic_array3, traits::RandomAccessMut}; /// Create a Lagrange element pub fn create( @@ -30,18 +28,22 @@ pub fn create( } for d in 0..tdim { for _e in 0..cell.entity_count(d) { - x[d].push(zero_matrix([0, tdim])); + x[d].push(rlst_dynamic_array2!(f64, [0, tdim])); m[d].push(rlst_dynamic_array3!(f64, [0, 1, 0])); } } - x[tdim].push(to_matrix(&cell.midpoint(), [1, tdim])); + let mut midp = rlst_dynamic_array2!(f64, [1, tdim]); + for (i, j) in cell.midpoint().iter().enumerate() { + *midp.get_mut([0, i]).unwrap() = *j; + } + x[tdim].push(midp); let mut mentry = rlst_dynamic_array3!(f64, [1, 1, 1]); *mentry.get_mut([0, 0, 0]).unwrap() = 1.0; m[tdim].push(mentry); } else { // TODO: GLL points for e in 0..cell.entity_count(0) { - let mut pts = zero_matrix([1, tdim]); + let mut pts = rlst_dynamic_array2!(f64, [1, tdim]); for i in 0..tdim { *pts.get_mut([0, i]).unwrap() = cell.vertices()[e * tdim + i]; } @@ -51,7 +53,7 @@ pub fn create( m[0].push(mentry); } for e in 0..cell.entity_count(1) { - let mut pts = zero_matrix([degree - 1, tdim]); + let mut pts = rlst_dynamic_array2!(f64, [degree - 1, tdim]); let vn0 = cell.edges()[2 * e]; let vn1 = cell.edges()[2 * e + 1]; let v0 = &cell.vertices()[vn0 * tdim..(vn0 + 1) * tdim]; @@ -82,7 +84,7 @@ pub fn create( } else { panic!("Unsupported face type"); }; - let mut pts = zero_matrix([npts, tdim]); + let mut pts = rlst_dynamic_array2!(f64, [npts, tdim]); let vn0 = cell.faces()[start]; let vn1 = cell.faces()[start + 1]; @@ -188,7 +190,11 @@ mod test { let e = create(ReferenceCellType::Interval, 0, Continuity::Discontinuous); assert_eq!(e.value_size(), 1); let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 4)); - let points = to_matrix(&[0.0, 0.2, 0.4, 1.0], [4, 1]); + let mut points = rlst_dynamic_array2!(f64, [4, 1]); + *points.get_mut([0, 0]).unwrap() = 0.0; + *points.get_mut([1, 0]).unwrap() = 0.2; + *points.get_mut([2, 0]).unwrap() = 0.4; + *points.get_mut([3, 0]).unwrap() = 1.0; e.tabulate(&points, 0, &mut data); for pt in 0..4 { @@ -202,7 +208,11 @@ mod test { let e = create(ReferenceCellType::Interval, 1, Continuity::Continuous); assert_eq!(e.value_size(), 1); let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 4)); - let points = to_matrix(&[0.0, 0.2, 0.4, 1.0], [4, 1]); + let mut points = rlst_dynamic_array2!(f64, [4, 1]); + *points.get_mut([0, 0]).unwrap() = 0.0; + *points.get_mut([1, 0]).unwrap() = 0.2; + *points.get_mut([2, 0]).unwrap() = 0.4; + *points.get_mut([3, 0]).unwrap() = 1.0; e.tabulate(&points, 0, &mut data); for pt in 0..4 { @@ -223,10 +233,21 @@ mod test { let e = create(ReferenceCellType::Triangle, 0, Continuity::Discontinuous); assert_eq!(e.value_size(), 1); let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 6)); - let points = to_matrix( - &[0.0, 1.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.0, 1.0, 0.0, 0.5, 0.5], - [6, 2], - ); + + let mut points = rlst_dynamic_array2!(f64, [6, 2]); + *points.get_mut([0, 0]).unwrap() = 0.0; + *points.get_mut([0, 1]).unwrap() = 0.0; + *points.get_mut([1, 0]).unwrap() = 1.0; + *points.get_mut([1, 1]).unwrap() = 0.0; + *points.get_mut([2, 0]).unwrap() = 0.0; + *points.get_mut([2, 1]).unwrap() = 1.0; + *points.get_mut([3, 0]).unwrap() = 0.5; + *points.get_mut([3, 1]).unwrap() = 0.0; + *points.get_mut([4, 0]).unwrap() = 0.0; + *points.get_mut([4, 1]).unwrap() = 0.5; + *points.get_mut([5, 0]).unwrap() = 0.5; + *points.get_mut([5, 1]).unwrap() = 0.5; + e.tabulate(&points, 0, &mut data); for pt in 0..6 { @@ -240,10 +261,19 @@ mod test { let e = create(ReferenceCellType::Triangle, 1, Continuity::Continuous); assert_eq!(e.value_size(), 1); let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 6)); - let points = to_matrix( - &[0.0, 1.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.0, 1.0, 0.0, 0.5, 0.5], - [6, 2], - ); + let mut points = rlst_dynamic_array2!(f64, [6, 2]); + *points.get_mut([0, 0]).unwrap() = 0.0; + *points.get_mut([0, 1]).unwrap() = 0.0; + *points.get_mut([1, 0]).unwrap() = 1.0; + *points.get_mut([1, 1]).unwrap() = 0.0; + *points.get_mut([2, 0]).unwrap() = 0.0; + *points.get_mut([2, 1]).unwrap() = 1.0; + *points.get_mut([3, 0]).unwrap() = 0.5; + *points.get_mut([3, 1]).unwrap() = 0.0; + *points.get_mut([4, 0]).unwrap() = 0.0; + *points.get_mut([4, 1]).unwrap() = 0.5; + *points.get_mut([5, 0]).unwrap() = 0.5; + *points.get_mut([5, 1]).unwrap() = 0.5; e.tabulate(&points, 0, &mut data); for pt in 0..6 { @@ -327,10 +357,19 @@ mod test { ); assert_eq!(e.value_size(), 1); let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 6)); - let points = to_matrix( - &[0.0, 1.0, 0.0, 1.0, 0.25, 0.3, 0.0, 0.0, 1.0, 1.0, 0.5, 0.2], - [6, 2], - ); + let mut points = rlst_dynamic_array2!(f64, [6, 2]); + *points.get_mut([0, 0]).unwrap() = 0.0; + *points.get_mut([0, 1]).unwrap() = 0.0; + *points.get_mut([1, 0]).unwrap() = 1.0; + *points.get_mut([1, 1]).unwrap() = 0.0; + *points.get_mut([2, 0]).unwrap() = 0.0; + *points.get_mut([2, 1]).unwrap() = 1.0; + *points.get_mut([3, 0]).unwrap() = 0.5; + *points.get_mut([3, 1]).unwrap() = 0.0; + *points.get_mut([4, 0]).unwrap() = 0.0; + *points.get_mut([4, 1]).unwrap() = 0.5; + *points.get_mut([5, 0]).unwrap() = 0.5; + *points.get_mut([5, 1]).unwrap() = 0.5; e.tabulate(&points, 0, &mut data); for pt in 0..6 { @@ -344,10 +383,20 @@ mod test { let e = create(ReferenceCellType::Quadrilateral, 1, Continuity::Continuous); assert_eq!(e.value_size(), 1); let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 6)); - let points = to_matrix( - &[0.0, 1.0, 0.0, 1.0, 0.25, 0.3, 0.0, 0.0, 1.0, 1.0, 0.5, 0.2], - [6, 2], - ); + let mut points = rlst_dynamic_array2!(f64, [6, 2]); + *points.get_mut([0, 0]).unwrap() = 0.0; + *points.get_mut([0, 1]).unwrap() = 0.0; + *points.get_mut([1, 0]).unwrap() = 1.0; + *points.get_mut([1, 1]).unwrap() = 0.0; + *points.get_mut([2, 0]).unwrap() = 0.0; + *points.get_mut([2, 1]).unwrap() = 1.0; + *points.get_mut([3, 0]).unwrap() = 1.0; + *points.get_mut([3, 1]).unwrap() = 1.0; + *points.get_mut([4, 0]).unwrap() = 0.25; + *points.get_mut([4, 1]).unwrap() = 0.5; + *points.get_mut([5, 0]).unwrap() = 0.3; + *points.get_mut([5, 1]).unwrap() = 0.2; + e.tabulate(&points, 0, &mut data); for pt in 0..6 { @@ -376,10 +425,19 @@ mod test { let e = create(ReferenceCellType::Quadrilateral, 2, Continuity::Continuous); assert_eq!(e.value_size(), 1); let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 6)); - let points = to_matrix( - &[0.0, 1.0, 0.0, 1.0, 0.25, 0.3, 0.0, 0.0, 1.0, 1.0, 0.5, 0.2], - [6, 2], - ); + let mut points = rlst_dynamic_array2!(f64, [6, 2]); + *points.get_mut([0, 0]).unwrap() = 0.0; + *points.get_mut([0, 1]).unwrap() = 0.0; + *points.get_mut([1, 0]).unwrap() = 1.0; + *points.get_mut([1, 1]).unwrap() = 0.0; + *points.get_mut([2, 0]).unwrap() = 0.0; + *points.get_mut([2, 1]).unwrap() = 1.0; + *points.get_mut([3, 0]).unwrap() = 1.0; + *points.get_mut([3, 1]).unwrap() = 1.0; + *points.get_mut([4, 0]).unwrap() = 0.25; + *points.get_mut([4, 1]).unwrap() = 0.5; + *points.get_mut([5, 0]).unwrap() = 0.3; + *points.get_mut([5, 1]).unwrap() = 0.2; e.tabulate(&points, 0, &mut data); for pt in 0..6 { diff --git a/element/src/element/raviart_thomas.rs b/element/src/element/raviart_thomas.rs index 5b7c035e..c80c1559 100644 --- a/element/src/element/raviart_thomas.rs +++ b/element/src/element/raviart_thomas.rs @@ -2,11 +2,9 @@ use crate::element::{create_cell, CiarletElement}; use crate::polynomials::polynomial_count; -use bempp_tools::arrays::zero_matrix; use bempp_traits::cell::ReferenceCellType; use bempp_traits::element::{Continuity, ElementFamily, MapType}; -use rlst_dense::rlst_dynamic_array3; -use rlst_dense::traits::RandomAccessMut; +use rlst_dense::{rlst_dynamic_array2, rlst_dynamic_array3, traits::RandomAccessMut}; /// Create a Raviart-Thomas element pub fn create( @@ -46,12 +44,12 @@ 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(zero_matrix([0, tdim])); + x[0].push(rlst_dynamic_array2!(f64, [0, tdim])); m[0].push(rlst_dynamic_array3!(f64, [0, 2, 0])); } for e in 0..cell.entity_count(1) { - let mut pts = zero_matrix([1, tdim]); + let mut pts = rlst_dynamic_array2!(f64, [1, tdim]); let mut mat = rlst_dynamic_array3!(f64, [1, 2, 1]); let vn0 = cell.edges()[2 * e]; let vn1 = cell.edges()[2 * e + 1]; @@ -67,7 +65,7 @@ pub fn create( } for _e in 0..cell.entity_count(2) { - x[2].push(zero_matrix([0, tdim])); + x[2].push(rlst_dynamic_array2!(f64, [0, tdim])); m[2].push(rlst_dynamic_array3!(f64, [0, 2, 0])) } @@ -90,7 +88,6 @@ mod test { use crate::cell::*; use crate::element::raviart_thomas::*; use approx::*; - use bempp_tools::arrays::to_matrix; use bempp_traits::element::FiniteElement; use rlst_dense::rlst_dynamic_array4; use rlst_dense::traits::RandomAccessByRef; @@ -130,10 +127,19 @@ mod test { let e = create(ReferenceCellType::Triangle, 1, Continuity::Continuous); assert_eq!(e.value_size(), 2); let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 6)); - let points = to_matrix( - &[0.0, 1.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.0, 1.0, 0.0, 0.5, 0.5], - [6, 2], - ); + let mut points = rlst_dynamic_array2!(f64, [6, 2]); + *points.get_mut([0, 0]).unwrap() = 0.0; + *points.get_mut([0, 1]).unwrap() = 0.0; + *points.get_mut([1, 0]).unwrap() = 1.0; + *points.get_mut([1, 1]).unwrap() = 0.0; + *points.get_mut([2, 0]).unwrap() = 0.0; + *points.get_mut([2, 1]).unwrap() = 1.0; + *points.get_mut([3, 0]).unwrap() = 0.5; + *points.get_mut([3, 1]).unwrap() = 0.0; + *points.get_mut([4, 0]).unwrap() = 0.0; + *points.get_mut([4, 1]).unwrap() = 0.5; + *points.get_mut([5, 0]).unwrap() = 0.5; + *points.get_mut([5, 1]).unwrap() = 0.5; e.tabulate(&points, 0, &mut data); for pt in 0..6 { diff --git a/element/src/polynomials.rs b/element/src/polynomials.rs index cfce8498..b4d4ae13 100644 --- a/element/src/polynomials.rs +++ b/element/src/polynomials.rs @@ -444,16 +444,18 @@ mod test { use crate::polynomials::*; use approx::*; use bempp_quadrature::simplex_rules::simplex_rule; - use bempp_tools::arrays::{transpose_to_matrix, zero_matrix}; - use rlst_dense::rlst_dynamic_array3; - use rlst_dense::traits::RandomAccessMut; + use rlst_dense::{rlst_dynamic_array2, rlst_dynamic_array3, traits::RandomAccessMut}; #[test] fn test_legendre_interval() { let degree = 6; let rule = simplex_rule(ReferenceCellType::Interval, degree + 1).unwrap(); - let points = transpose_to_matrix(&rule.points, [rule.npoints, 1]); + + let mut points = rlst_dynamic_array2!(f64, [rule.npoints, 1]); + for (i, j) in rule.points.iter().enumerate() { + *points.get_mut([i, 0]).unwrap() = *j; + } let mut data = rlst_dynamic_array3!( f64, @@ -483,7 +485,12 @@ mod test { let degree = 5; let rule = simplex_rule(ReferenceCellType::Triangle, 79).unwrap(); - let points = transpose_to_matrix(&rule.points, [rule.npoints, 2]); + let mut points = rlst_dynamic_array2!(f64, [rule.npoints, 2]); + for i in 0..rule.npoints { + for j in 0..2 { + *points.get_mut([i, j]).unwrap() = rule.points[i * 2 + j]; + } + } let mut data = rlst_dynamic_array3!( f64, @@ -513,7 +520,12 @@ mod test { let degree = 5; let rule = simplex_rule(ReferenceCellType::Quadrilateral, 85).unwrap(); - let points = transpose_to_matrix(&rule.points, [rule.npoints, 2]); + let mut points = rlst_dynamic_array2!(f64, [rule.npoints, 2]); + for i in 0..rule.npoints { + for j in 0..2 { + *points.get_mut([i, j]).unwrap() = rule.points[i * 2 + j]; + } + } let mut data = rlst_dynamic_array3!( f64, @@ -549,7 +561,7 @@ mod test { let degree = 6; let epsilon = 1e-10; - let mut points = zero_matrix([20, 1]); + let mut points = rlst_dynamic_array2!(f64, [20, 1]); for i in 0..10 { *points.get_mut([2 * i, 0]).unwrap() = i as f64 / 10.0; *points.get_mut([2 * i + 1, 0]).unwrap() = points.get([2 * i, 0]).unwrap() + epsilon; @@ -578,7 +590,7 @@ mod test { let degree = 6; let epsilon = 1e-10; - let mut points = zero_matrix([165, 2]); + let mut points = rlst_dynamic_array2!(f64, [165, 2]); let mut index = 0; for i in 0..10 { for j in 0..10 - i { @@ -623,7 +635,7 @@ mod test { let degree = 6; let epsilon = 1e-10; - let mut points = zero_matrix([300, 2]); + let mut points = rlst_dynamic_array2!(f64, [300, 2]); for i in 0..10 { for j in 0..10 { let index = 10 * i + j; @@ -672,7 +684,7 @@ mod test { fn test_legendre_interval_against_known_polynomials() { let degree = 3; - let mut points = zero_matrix([11, 1]); + let mut points = rlst_dynamic_array2!(f64, [11, 1]); for i in 0..11 { *points.get_mut([i, 0]).unwrap() = i as f64 / 10.0; } @@ -752,7 +764,7 @@ mod test { fn test_legendre_quadrilateral_against_known_polynomials() { let degree = 2; - let mut points = zero_matrix([121, 2]); + let mut points = rlst_dynamic_array2!(f64, [121, 2]); for i in 0..11 { for j in 0..11 { *points.get_mut([11 * i + j, 0]).unwrap() = i as f64 / 10.0; diff --git a/field/src/array.rs b/field/src/array.rs index 00528be4..111bcf03 100644 --- a/field/src/array.rs +++ b/field/src/array.rs @@ -4,9 +4,11 @@ use itertools::Itertools; use num::traits::Num; -use bempp_tools::Array3D; use rlst_common::types::Scalar; use rlst_dense::{ + array::Array, + base_array::BaseArray, + data_container::VectorContainer, rlst_dynamic_array3, traits::{RandomAccessByRef, RandomAccessMut, Shape}, }; @@ -30,10 +32,10 @@ pub fn argsort(arr: &[T]) -> Vec { /// * `pad_size` - The amount of padding to be added along each axis. /// * `pad_index` - The position in the array to start the padding from. pub fn pad3( - arr: &Array3D, + arr: &Array, 3>, 3>, pad_size: (usize, usize, usize), pad_index: (usize, usize, usize), -) -> Array3D +) -> Array, 3>, 3> where T: Clone + Copy + Num, { @@ -62,7 +64,9 @@ where /// /// # Arguments /// * `arr` - An array to be flipped. -pub fn flip3(arr: &Array3D) -> Array3D +pub fn flip3( + arr: &Array, 3>, 3>, +) -> Array, 3>, 3> where T: Clone + Copy + Num, { diff --git a/field/src/field.rs b/field/src/field.rs index a9bf7a0a..f5c4c533 100644 --- a/field/src/field.rs +++ b/field/src/field.rs @@ -17,7 +17,6 @@ use rlst_dense::{ }; use std::collections::HashSet; -use bempp_tools::arrays::Array3D; use bempp_traits::{field::FieldTranslationData, kernel::Kernel, types::EvalType}; use bempp_tree::{ implementations::helpers::find_corners, types::domain::Domain, types::morton::MortonKey, @@ -504,7 +503,7 @@ where order: usize, convolution_grid: &[T], target_pt: [T; 3], - ) -> Array3D { + ) -> Array, 3>, 3> { let n = 2 * order - 1; let npad = n + 1; @@ -538,7 +537,11 @@ where /// # Arguments /// * `order` - The expansion order for the multipole and local expansions. /// * `charges` - A vector of charges. - pub fn compute_signal(&self, order: usize, charges: &[T]) -> Array3D { + pub fn compute_signal( + &self, + order: usize, + charges: &[T], + ) -> Array, 3>, 3> { let n = 2 * order - 1; let npad = n + 1; diff --git a/grid/examples/curved_cells.rs b/grid/examples/curved_cells.rs index 7d7b4e9a..81333dc9 100644 --- a/grid/examples/curved_cells.rs +++ b/grid/examples/curved_cells.rs @@ -1,21 +1,21 @@ use bempp_grid::grid::SerialGrid; use bempp_grid::io::export_as_gmsh; -use bempp_tools::arrays::{to_matrix, AdjacencyList}; +use bempp_tools::arrays::AdjacencyList; use bempp_traits::arrays::AdjacencyListAccess; use bempp_traits::cell::ReferenceCellType; use bempp_traits::grid::{Geometry, Grid, Topology}; +use rlst_dense::{rlst_dynamic_array2, traits::RawAccessMut}; fn main() { // Create a grid of three cells. One cell is a curved triangle, one cell is a flat triangle, the other is a curved quadrilateral + let mut points = rlst_dynamic_array2!(f64, [13, 3]); + points.data_mut().copy_from_slice(&vec![ + 0.0, 0.5, 1.0, 1.5, 0.0, 0.5, 1.0, 2.0, 1.5, 0.0, 0.5, 1.0, 2.0, 0.0, 0.0, 0.0, 0.25, 0.5, + 0.5, 0.5, 0.5, 0.75, 1.0, 1.0, 1.0, -0.5, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5, 0.5, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, + ]); let grid = SerialGrid::new( - to_matrix( - &vec![ - 0.0, 0.5, 1.0, 1.5, 0.0, 0.5, 1.0, 2.0, 1.5, 0.0, 0.5, 1.0, 2.0, 0.0, 0.0, 0.0, - 0.25, 0.5, 0.5, 0.5, 0.5, 0.75, 1.0, 1.0, 1.0, -0.5, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5, - 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - ], - [13, 3], - ), + points, AdjacencyList::from_data( vec![2, 7, 12, 0, 2, 9, 11, 1, 4, 6, 10, 5, 2, 7, 11, 8, 6, 3], vec![0, 3, 12, 18], diff --git a/grid/examples/parallel_grid.rs b/grid/examples/parallel_grid.rs index 06d7456b..fa3000e0 100644 --- a/grid/examples/parallel_grid.rs +++ b/grid/examples/parallel_grid.rs @@ -5,7 +5,7 @@ use approx::*; #[cfg(feature = "mpi")] use bempp_grid::parallel_grid::ParallelGrid; #[cfg(feature = "mpi")] -use bempp_tools::arrays::{zero_matrix, AdjacencyList}; +use bempp_tools::arrays::AdjacencyList; #[cfg(feature = "mpi")] use bempp_traits::arrays::AdjacencyListAccess; #[cfg(feature = "mpi")] @@ -15,7 +15,7 @@ use bempp_traits::grid::{Geometry, Grid, Ownership, Topology}; #[cfg(feature = "mpi")] use mpi::{environment::Universe, request::WaitGuard, topology::Communicator, traits::*}; #[cfg(feature = "mpi")] -use rlst_dense::traits::RandomAccessMut; +use rlst_dense::{rlst_dynamic_array2, traits::RandomAccessMut}; #[cfg(feature = "mpi")] fn test_parallel_grid() { @@ -29,7 +29,7 @@ fn test_parallel_grid() { let n = 10; let grid = if rank == 0 { - let mut pts = zero_matrix([n * n, 3]); + let mut pts = rlst_dynamic_array2!(f64, [n * n, 3]); let mut i = 0; for y in 0..n { for x in 0..n { diff --git a/grid/src/grid.rs b/grid/src/grid.rs index 0a676334..be626336 100644 --- a/grid/src/grid.rs +++ b/grid/src/grid.rs @@ -1,15 +1,20 @@ //! A serial implementation of a grid use bempp_element::cell; use bempp_element::element::{create_element, CiarletElement}; -use bempp_tools::arrays::{zero_matrix, AdjacencyList, Array4D, Mat}; +use bempp_tools::arrays::AdjacencyList; use bempp_traits::arrays::AdjacencyListAccess; use bempp_traits::cell::{ReferenceCell, ReferenceCellType}; use bempp_traits::element::{Continuity, ElementFamily, FiniteElement}; use bempp_traits::grid::{Geometry, GeometryEvaluator, Grid, Ownership, Topology}; use itertools::izip; -use rlst_dense::rlst_dynamic_array4; -use rlst_dense::traits::{ - RandomAccessByRef, RandomAccessMut, Shape, UnsafeRandomAccessByRef, UnsafeRandomAccessMut, +use rlst_dense::{ + array::Array, + base_array::BaseArray, + data_container::VectorContainer, + rlst_dynamic_array2, rlst_dynamic_array4, + traits::{ + RandomAccessByRef, RandomAccessMut, Shape, UnsafeRandomAccessByRef, UnsafeRandomAccessMut, + }, }; use rlst_proc_macro::rlst_static_array; use std::cell::RefCell; @@ -17,24 +22,24 @@ use std::ptr; pub struct EvaluatorTdim2Gdim3<'a> { geometry: &'a SerialGeometry, - points: &'a Mat, - table: Array4D, + points: &'a Array, 2>, 2>, + table: Array, 4>, 4>, npts: usize, - axes: RefCell>, + axes: RefCell, 2>, 2>>, } impl<'a> EvaluatorTdim2Gdim3<'a> { pub fn new( geometry: &'a SerialGeometry, element: &impl FiniteElement, - points: &'a Mat, + points: &'a Array, 2>, 2>, ) -> Self { let npts = points.shape()[0]; assert_eq!(points.shape()[1], 2); assert_eq!(geometry.dim(), 3); let mut table = rlst_dynamic_array4!(f64, element.tabulate_array_shape(1, npts)); element.tabulate(points, 1, &mut table); - let axes = RefCell::new(zero_matrix([2, 3])); + let axes = RefCell::new(rlst_dynamic_array2!(f64, [2, 3])); Self { geometry, @@ -46,12 +51,21 @@ impl<'a> EvaluatorTdim2Gdim3<'a> { } } -impl<'a> GeometryEvaluator, Mat> for EvaluatorTdim2Gdim3<'a> { - fn points(&self) -> &Mat { +impl<'a> + GeometryEvaluator< + Array, 2>, 2>, + Array, 2>, 2>, + > for EvaluatorTdim2Gdim3<'a> +{ + fn points(&self) -> &Array, 2>, 2> { self.points } - fn compute_points(&self, cell_index: usize, points: &mut Mat) { + fn compute_points( + &self, + cell_index: usize, + points: &mut Array, 2>, 2>, + ) { for i in 0..3 { for p in 0..self.npts { unsafe { @@ -76,7 +90,7 @@ impl<'a> GeometryEvaluator, Mat> for EvaluatorTdim2Gdim3<'a> { fn compute_normals_and_jacobian_determinants( &self, cell_index: usize, - normals: &mut Mat, + normals: &mut Array, 2>, 2>, jdets: &mut [f64], ) { let mut axes = self.axes.borrow_mut(); @@ -125,8 +139,8 @@ impl<'a> GeometryEvaluator, Mat> for EvaluatorTdim2Gdim3<'a> { pub struct LinearSimplexEvaluatorTdim2Gdim3<'a> { geometry: &'a SerialGeometry, - points: &'a Mat, - table: Array4D, + points: &'a Array, 2>, 2>, + table: Array, 4>, 4>, npts: usize, js: RefCell<[f64; 6]>, } @@ -135,7 +149,7 @@ impl<'a> LinearSimplexEvaluatorTdim2Gdim3<'a> { pub fn new( geometry: &'a SerialGeometry, element: &impl FiniteElement, - points: &'a Mat, + points: &'a Array, 2>, 2>, ) -> Self { let npts = points.shape()[0]; assert_eq!(points.shape()[1], 2); @@ -172,12 +186,21 @@ impl<'a> LinearSimplexEvaluatorTdim2Gdim3<'a> { } } -impl<'a> GeometryEvaluator, Mat> for LinearSimplexEvaluatorTdim2Gdim3<'a> { - fn points(&self) -> &Mat { +impl<'a> + GeometryEvaluator< + Array, 2>, 2>, + Array, 2>, 2>, + > for LinearSimplexEvaluatorTdim2Gdim3<'a> +{ + fn points(&self) -> &Array, 2>, 2> { self.points } - fn compute_points(&self, cell_index: usize, points: &mut Mat) { + fn compute_points( + &self, + cell_index: usize, + points: &mut Array, 2>, 2>, + ) { for j in 0..3 { for p in 0..self.npts { let mut sum = 0.0; @@ -198,7 +221,7 @@ impl<'a> GeometryEvaluator, Mat> for LinearSimplexEvaluatorTdim2Gd fn compute_normals_and_jacobian_determinants( &self, cell_index: usize, - normals: &mut Mat, + normals: &mut Array, 2>, 2>, jdets: &mut [f64], ) { self.single_jacobian(cell_index); @@ -234,7 +257,7 @@ impl<'a> GeometryEvaluator, Mat> for LinearSimplexEvaluatorTdim2Gd /// Geometry of a serial grid pub struct SerialGeometry { coordinate_elements: Vec, - coordinates: Mat, + coordinates: Array, 2>, 2>, cells: AdjacencyList, element_changes: Vec, index_map: Vec, @@ -257,7 +280,7 @@ fn element_from_npts(cell_type: ReferenceCellType, npts: usize) -> CiarletElemen impl SerialGeometry { pub fn new( - coordinates: Mat, + coordinates: Array, 2>, 2>, cells: &AdjacencyList, cell_types: &[ReferenceCellType], ) -> Self { @@ -317,8 +340,8 @@ impl SerialGeometry { } } impl Geometry for SerialGeometry { - type T = Mat; - type TMut = Mat; + type T = Array, 2>, 2>; + type TMut = Array, 2>, 2>; fn dim(&self) -> usize { self.coordinates.shape()[1] } @@ -501,7 +524,7 @@ impl Geometry for SerialGeometry { if points.shape()[0] != jacobian_determinants.len() { panic!("jacobian_determinants has wrong length."); } - let mut js = zero_matrix([npts, gdim * tdim]); + let mut js = rlst_dynamic_array2!(f64, [npts, gdim * tdim]); self.compute_jacobians(points, cell, &mut js); for (p, jdet) in jacobian_determinants.iter_mut().enumerate() { @@ -584,7 +607,7 @@ impl Geometry for SerialGeometry { && element.degree() == 1 { // Map is affine - let mut js = zero_matrix([npts, gdim * tdim]); + let mut js = rlst_dynamic_array2!(f64, [npts, gdim * tdim]); self.compute_jacobians(points, cell, &mut js); for p in 0..npts { @@ -988,7 +1011,7 @@ pub struct SerialGrid { impl SerialGrid { pub fn new( - coordinates: Mat, + coordinates: Array, 2>, 2>, cells: AdjacencyList, cell_types: Vec, ) -> Self { @@ -1029,7 +1052,16 @@ mod test { use crate::grid::*; use crate::shapes::regular_sphere; use approx::*; - use bempp_tools::arrays::to_matrix; + use rlst_dense::traits::RawAccessMut; + + fn to_matrix( + slice: &[f64], + shape: [usize; 2], + ) -> Array, 2>, 2> { + let mut mat = rlst_dynamic_array2!(f64, shape); + mat.data_mut().copy_from_slice(slice); + mat + } #[test] fn test_connectivity() { @@ -1287,7 +1319,7 @@ mod test { ); // Test compute_points - let mut physical_points = zero_matrix([points.shape()[0], 3]); + let mut physical_points = rlst_dynamic_array2!(f64, [points.shape()[0], 3]); g.geometry() .compute_points(&points, 0, &mut physical_points); assert_relative_eq!( @@ -1414,7 +1446,7 @@ mod test { ); // Test compute_jacobians - let mut jacobians = zero_matrix([points.shape()[0], 6]); + let mut jacobians = rlst_dynamic_array2!(f64, [points.shape()[0], 6]); g.geometry().compute_jacobians(&points, 0, &mut jacobians); for i in 0..3 { assert_relative_eq!(*jacobians.get([i, 0]).unwrap(), 2.0, max_relative = 1e-14); @@ -1448,7 +1480,7 @@ mod test { } // Test compute_jacobian_inverses - let mut jinvs = zero_matrix([points.shape()[0], 6]); + let mut jinvs = rlst_dynamic_array2!(f64, [points.shape()[0], 6]); g.geometry() .compute_jacobian_inverses(&points, 0, &mut jinvs); for i in 0..3 { @@ -1492,7 +1524,7 @@ mod test { let pt = to_matrix(&[1.0 / 3.0, 1.0 / 3.0], [1, 2]); - let mut normal = zero_matrix([1, 3]); + let mut normal = rlst_dynamic_array2!(f64, [1, 3]); g.geometry().compute_normals(&pt, 0, &mut normal); assert_relative_eq!(*normal.get([0, 0]).unwrap(), 0.0); @@ -1524,7 +1556,7 @@ mod test { ); let points = to_matrix(&[0.0, 0.2, 0.5, 0.7, 1.0, 0.0, 0.3, 0.9, 1.0, 0.3], [5, 2]); - let mut normals = zero_matrix([5, 3]); + let mut normals = rlst_dynamic_array2!(f64, [5, 3]); curved_g .geometry() @@ -1595,8 +1627,8 @@ mod test { let pts = to_matrix(&[0.1, 0.2, 0.6, 0.1, 0.4, 0.2], [3, 2]); let e = grid.geometry().get_evaluator(&element, &pts); - let mut points0 = zero_matrix([3, 3]); - let mut points1 = zero_matrix([3, 3]); + let mut points0 = rlst_dynamic_array2!(f64, [3, 3]); + let mut points1 = rlst_dynamic_array2!(f64, [3, 3]); for c in 0..grid.geometry().cell_count() { grid.geometry().compute_points(&pts, c, &mut points0); e.compute_points(c, &mut points1); @@ -1624,8 +1656,8 @@ mod test { let pts = to_matrix(&[0.1, 0.2, 0.6, 0.1, 0.4, 0.2], [3, 2]); let e = grid.geometry().get_evaluator(&element, &pts); - let mut normals0 = zero_matrix([3, 3]); - let mut normals1 = zero_matrix([3, 3]); + let mut normals0 = rlst_dynamic_array2!(f64, [3, 3]); + let mut normals1 = rlst_dynamic_array2!(f64, [3, 3]); let mut jacobian_determinants0 = vec![0.0; 3]; let mut jacobian_determinants1 = vec![0.0; 3]; diff --git a/grid/src/io.rs b/grid/src/io.rs index ab4e3c02..58282ef4 100644 --- a/grid/src/io.rs +++ b/grid/src/io.rs @@ -122,8 +122,21 @@ mod test { use crate::grid::SerialGrid; use crate::io::*; use crate::shapes::regular_sphere; - use bempp_tools::arrays::{to_matrix, AdjacencyList}; + use bempp_tools::arrays::AdjacencyList; use bempp_traits::cell::ReferenceCellType; + use rlst_dense::{ + array::Array, base_array::BaseArray, data_container::VectorContainer, rlst_dynamic_array2, + traits::RawAccessMut, + }; + + fn to_matrix( + slice: &[f64], + shape: [usize; 2], + ) -> Array, 2>, 2> { + let mut mat = rlst_dynamic_array2!(f64, shape); + mat.data_mut().copy_from_slice(slice); + mat + } #[test] fn test_gmsh_output_regular_sphere() { diff --git a/grid/src/parallel_grid.rs b/grid/src/parallel_grid.rs index 6c0e4205..3c17fd60 100644 --- a/grid/src/parallel_grid.rs +++ b/grid/src/parallel_grid.rs @@ -1,13 +1,19 @@ //! A parallel implementation of a grid use crate::grid::{SerialGeometry, SerialTopology}; use bempp_element::element::CiarletElement; -use bempp_tools::arrays::{zero_matrix, AdjacencyList, Mat}; +use bempp_tools::arrays::AdjacencyList; use bempp_traits::arrays::AdjacencyListAccess; use bempp_traits::cell::ReferenceCellType; use bempp_traits::element::FiniteElement; use bempp_traits::grid::{Geometry, GeometryEvaluator, Grid, Ownership, Topology}; use mpi::{request::WaitGuard, topology::Communicator, traits::*}; -use rlst_dense::traits::{RandomAccessByRef, RandomAccessMut, Shape}; +use rlst_dense::{ + array::Array, + base_array::BaseArray, + data_container::VectorContainer, + rlst_dynamic_array2, + traits::{RandomAccessByRef, RandomAccessMut, Shape}, +}; /// Geometry of a parallel grid pub struct ParallelGeometry<'a, C: Communicator> { @@ -18,7 +24,7 @@ pub struct ParallelGeometry<'a, C: Communicator> { impl<'a, C: Communicator> ParallelGeometry<'a, C> { pub fn new( comm: &'a C, - coordinates: Mat, + coordinates: Array, 2>, 2>, cells: &AdjacencyList, cell_types: &Vec, ) -> Self { @@ -42,8 +48,8 @@ impl<'a, C: Communicator> ParallelGeometry<'a, C> { } impl<'a, C: Communicator> Geometry for ParallelGeometry<'a, C> { - type T = Mat; - type TMut = Mat; + type T = Array, 2>, 2>; + type TMut = Array, 2>, 2>; fn dim(&self) -> usize { self.serial_geometry.dim() @@ -203,7 +209,7 @@ pub struct ParallelGrid<'a, C: Communicator> { impl<'a, C: Communicator> ParallelGrid<'a, C> { pub fn new( comm: &'a C, - coordinates: Mat, + coordinates: Array, 2>, 2>, cells: AdjacencyList, cell_types: Vec, cell_owners: Vec, @@ -403,7 +409,7 @@ impl<'a, C: Communicator> ParallelGrid<'a, C> { panic!("Unsupported cell type"); }); } - let mut coordinates = zero_matrix([vertex_indices.len(), gdim]); + let mut coordinates = rlst_dynamic_array2!(f64, [vertex_indices.len(), gdim]); for i in 0..vertex_indices.len() { for j in 0..gdim { *coordinates.get_mut([i, j]).unwrap() = flat_coordinates[i * gdim + j]; diff --git a/grid/src/shapes.rs b/grid/src/shapes.rs index a98d178b..24946fd9 100644 --- a/grid/src/shapes.rs +++ b/grid/src/shapes.rs @@ -2,11 +2,14 @@ use crate::grid::SerialGrid; use bempp_element::cell::Triangle; -use bempp_tools::arrays::{to_matrix, zero_matrix, AdjacencyList}; +use bempp_tools::arrays::AdjacencyList; use bempp_traits::arrays::AdjacencyListAccess; use bempp_traits::cell::{ReferenceCell, ReferenceCellType}; use bempp_traits::grid::{Geometry, Grid, Topology}; -use rlst_dense::traits::{RandomAccessByRef, RandomAccessMut}; +use rlst_dense::{ + rlst_dynamic_array2, + traits::{RandomAccessByRef, RandomAccessMut, RawAccessMut}, +}; /// Create a regular sphere /// @@ -14,14 +17,12 @@ use rlst_dense::traits::{RandomAccessByRef, RandomAccessMut}; /// Each time the grid is refined, each triangle is split into four triangles (by adding lines connecting the midpoints of /// each edge). The new points are then scaled so that they are a distance of 1 from the origin. pub fn regular_sphere(refinement_level: usize) -> SerialGrid { + let mut points = rlst_dynamic_array2!(f64, [6, 3]); + points.data_mut().copy_from_slice(&[ + 0.0, 1.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, -1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, -1.0, + ]); let mut g = SerialGrid::new( - to_matrix( - &[ - 0.0, 1.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, -1.0, 0.0, 1.0, 0.0, 0.0, 0.0, - 0.0, -1.0, - ], - [6, 3], - ), + points, AdjacencyList::from_data( vec![ 0, 1, 2, 0, 2, 3, 0, 3, 4, 0, 4, 1, 5, 2, 1, 5, 3, 2, 5, 4, 3, 5, 1, 4, @@ -35,7 +36,7 @@ pub fn regular_sphere(refinement_level: usize) -> SerialGrid { let nvertices_old = g.topology().entity_count(0); let ncells_old = g.topology().entity_count(2); let nvertices = g.topology().entity_count(0) + g.topology().entity_count(1); - let mut coordinates = zero_matrix([nvertices, 3]); + let mut coordinates = rlst_dynamic_array2!(f64, [nvertices, 3]); let mut cells = AdjacencyList::::new(); for i in 0..ncells_old { @@ -109,10 +110,12 @@ mod test { fn test_normal_is_outward() { for i in 0..3 { let g = regular_sphere(i); - let points = to_matrix(&[1.0 / 3.0, 1.0 / 3.0], [1, 2]); + let mut points = rlst_dynamic_array2!(f64, [1, 2]); + *points.get_mut([0, 0]).unwrap() = 1.0 / 3.0; + *points.get_mut([0, 1]).unwrap() = 1.0 / 3.0; - let mut mapped_pt = zero_matrix([1, 3]); - let mut normal = zero_matrix([1, 3]); + let mut mapped_pt = rlst_dynamic_array2!(f64, [1, 3]); + let mut normal = rlst_dynamic_array2!(f64, [1, 3]); for i in 0..g.geometry().cell_count() { g.geometry().compute_points(&points, i, &mut mapped_pt); diff --git a/kernel/src/laplace_3d.rs b/kernel/src/laplace_3d.rs index 2e1835f0..8147f64e 100644 --- a/kernel/src/laplace_3d.rs +++ b/kernel/src/laplace_3d.rs @@ -365,13 +365,19 @@ mod test { use super::*; use approx::assert_relative_eq; - use bempp_tools::arrays::Mat; use bempp_traits::types::Scalar; use rand::prelude::*; - use rlst_dense::rlst_dynamic_array2; - use rlst_dense::traits::{RandomAccessByRef, RandomAccessMut, RawAccess, RawAccessMut, Shape}; - - fn copy(m_in: &Mat) -> Mat { + use rlst_dense::{ + array::Array, + base_array::BaseArray, + data_container::VectorContainer, + rlst_dynamic_array2, + traits::{RandomAccessByRef, RandomAccessMut, RawAccess, RawAccessMut, Shape}, + }; + + fn copy( + m_in: &Array, 2>, 2>, + ) -> Array, 2>, 2> { let mut m = rlst_dynamic_array2!(f64, m_in.shape()); for i in 0..m_in.shape()[0] { for j in 0..m_in.shape()[1] { @@ -381,7 +387,7 @@ mod test { m } - fn rand_mat(shape: [usize; 2]) -> Mat { + fn rand_mat(shape: [usize; 2]) -> Array, 2>, 2> { let mut m = rlst_dynamic_array2!(f64, shape); let mut rng = rand::thread_rng(); for i in 0..shape[0] { @@ -392,7 +398,7 @@ mod test { m } - fn rand_vec(size: usize) -> Mat { + fn rand_vec(size: usize) -> Array, 2>, 2> { let mut v = rlst_dynamic_array2!(f64, [size, 1]); let mut rng = rand::thread_rng(); for i in 0..size { diff --git a/tools/Cargo.toml b/tools/Cargo.toml index c2c1c0af..0f877963 100644 --- a/tools/Cargo.toml +++ b/tools/Cargo.toml @@ -24,5 +24,3 @@ libc = "0.2" num = "0.4" bempp-traits = { path = "../traits"} rayon = "1.7" -rlst-common = { git = "https://github.com/linalg-rs/rlst.git" } -rlst-dense = { git = "https://github.com/linalg-rs/rlst.git" } diff --git a/tools/src/arrays.rs b/tools/src/arrays.rs index b7917fd2..5f270591 100644 --- a/tools/src/arrays.rs +++ b/tools/src/arrays.rs @@ -1,39 +1,6 @@ //! Containers to store multi-dimensional data use bempp_traits::arrays::AdjacencyListAccess; use num::Num; -use rlst_common::types::Scalar; -use rlst_dense::{ - array::Array, base_array::BaseArray, data_container::VectorContainer, rlst_dynamic_array2, - traits::UnsafeRandomAccessMut, -}; - -pub type Mat = Array, 2>, 2>; -pub type Array3D = Array, 3>, 3>; -pub type Array4D = Array, 4>, 4>; - -pub fn to_matrix(data: &[T], shape: [usize; 2]) -> Mat { - let mut mat = rlst_dynamic_array2![T, shape]; - for (i, d) in data.iter().enumerate() { - unsafe { - *mat.get_unchecked_mut([i % shape[0], i / shape[0]]) = *d; - } - } - mat -} - -pub fn transpose_to_matrix(data: &[T], shape: [usize; 2]) -> Mat { - let mut mat = rlst_dynamic_array2![T, shape]; - for (i, d) in data.iter().enumerate() { - unsafe { - *mat.get_unchecked_mut([i / shape[1], i % shape[1]]) = *d; - } - } - mat -} - -pub fn zero_matrix(shape: [usize; 2]) -> Mat { - rlst_dynamic_array2![T, shape] -} /// An adjacency list ///