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
bempp · Jan 16, 2024 · 592b810 · 592b810
1 parent d7d52e2
commit 592b810
Show file tree

Hide file tree

Showing 18 changed files with 363 additions and 218 deletions.
diff --git 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::<f64>([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
@@ -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<f64>,
+    output: &mut Array<f64, BaseArray<f64, VectorContainer<f64>, 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::<f64>([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::<f64>([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
@@ -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,22 +92,22 @@ fn assemble_batch_singular<'a>(
     trial_space: &SerialFunctionSpace<'a>,
     test_space: &SerialFunctionSpace<'a>,
     cell_pairs: &[(usize, usize)],
-    trial_points: &Mat<f64>,
-    test_points: &Mat<f64>,
+    trial_points: &Array<f64, BaseArray<f64, VectorContainer<f64>, 2>, 2>,
+    test_points: &Array<f64, BaseArray<f64, VectorContainer<f64>, 2>, 2>,
     weights: &[f64],
-    trial_table: &Array4D<f64>,
-    test_table: &Array4D<f64>,
+    trial_table: &Array<f64, BaseArray<f64, VectorContainer<f64>, 4>, 4>,
+    test_table: &Array<f64, BaseArray<f64, VectorContainer<f64>, 4>, 4>,
 ) -> usize {
     let grid = test_space.grid();
     let mut k = vec![0.0];
 
     // 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<f64>,
+    trial_points: &Array<f64, BaseArray<f64, VectorContainer<f64>, 2>, 2>,
     trial_weights: &[f64],
-    test_points: &Mat<f64>,
+    test_points: &Array<f64, BaseArray<f64, VectorContainer<f64>, 2>, 2>,
     test_weights: &[f64],
-    trial_table: &Array4D<f64>,
-    test_table: &Array4D<f64>,
+    trial_table: &Array<f64, BaseArray<f64, VectorContainer<f64>, 4>, 4>,
+    test_table: &Array<f64, BaseArray<f64, VectorContainer<f64>, 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<f64>,
+    output: &mut Array<f64, BaseArray<f64, VectorContainer<f64>, 2>, 2>,
     kernel: &impl Kernel<T = f64>,
     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<f64>,
+    output: &mut Array<f64, BaseArray<f64, VectorContainer<f64>, 2>, 2>,
     kernel: &impl Kernel<T = f64>,
     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<f64>,
+    output: &mut Array<f64, BaseArray<f64, VectorContainer<f64>, 2>, 2>,
     kernel: &impl Kernel<T = f64>,
     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::<f64>([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
@@ -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<Array<f64, BaseArray<f64, VectorContainer<f64>, 2>, 2>>; 4];
+type EntityWeights = [Vec<Array<f64, BaseArray<f64, VectorContainer<f64>, 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<f64>,
+    coefficients: Array<f64, BaseArray<f64, VectorContainer<f64>, 3>, 3>,
     entity_dofs: [AdjacencyList<usize>; 4],
-    // interpolation_points: [Vec<Array2D<f64>>; 4],
-    // interpolation_weights: [Vec<Array3D<f64>>; 4],
+    // interpolation_points: EntityPoints,
+    // interpolation_weights: EntityWeights,
 }
 
 impl CiarletElement {
@@ -36,9 +44,9 @@ impl CiarletElement {
         cell_type: ReferenceCellType,
         degree: usize,
         value_shape: Vec<usize>,
-        polynomial_coeffs: Array3D<f64>,
-        interpolation_points: [Vec<Mat<f64>>; 4],
-        interpolation_weights: [Vec<Array3D<f64>>; 4],
+        polynomial_coeffs: Array<f64, BaseArray<f64, VectorContainer<f64>, 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<Mat<f64>>; 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() {