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<f64, rlst_dense::base_matrix::BaseMatrix<f64, rlst_dense::VectorContainer<f64>, 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<f64>,
+            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::<NPTS_TEST, NPTS_TRIAL>(
@@ -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<usize>,
         polynomial_coeffs: Array3D<f64>,
-        interpolation_points: [Vec<Array2D<f64>>; 4],
+        interpolation_points: [Vec<Mat<f64>>; 4],
         interpolation_weights: [Vec<Array3D<f64>>; 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<Mat<f64>>; 4] = [vec![], vec![], vec![], vec![]];
             let mut pn = 0;
-            let mut all_pts = Array2D::<f64>::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::<f64>::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<T: RandomAccessByRef<Item = f64> + Shape>(
         &self,
-        points: &impl Array2DAccess<'a, f64>,
+        points: &T,
         nderivs: usize,
         data: &mut impl Array4DAccess<f64>,
     ) {
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::<f64>::new((0, tdim)));
+                x[d].push(to_matrix(&[], (0, tdim)));
                 m[d].push(Array3D::<f64>::new((0, 1, 0)));
             }
         }
-        x[tdim].push(Array2D::<f64>::from_data(cell.midpoint(), (1, tdim)));
+        x[tdim].push(to_matrix(&cell.midpoint(), (1, tdim)));
         m[tdim].push(Array3D::<f64>::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::<f64>::from_data(pts, (1, tdim)));
+            x[0].push(to_matrix(&pts, (1, tdim)));
             m[0].push(Array3D::<f64>::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::<f64>::from_data(pts, (degree - 1, tdim)));
+            x[1].push(to_matrix(&pts, (degree - 1, tdim)));
             m[1].push(Array3D::<f64>::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::<f64>::from_data(pts, (npts, tdim)));
+            x[2].push(to_matrix(&pts, (npts, tdim)));
             m[2].push(Array3D::<f64>::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::<f64>::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::<f64>::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::<f64>::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::<f64>::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::<f64>::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::<f64>::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::<f64>::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::<f64>::new((0, tdim)));
+        x[0].push(to_matrix(&[], (0, tdim)));
         m[0].push(Array3D::<f64>::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::<f64>::from_data(pts, (1, tdim)));
+        x[1].push(to_matrix(&pts, (1, tdim)));
         m[1].push(Array3D::<f64>::from_data(mat, (1, 2, 1)));
     }
 
     for _e in 0..cell.entity_count(2) {
-        x[2].push(Array2D::<f64>::new((0, tdim)));
+        x[2].push(to_matrix(&[], (0, tdim)));
         m[2].push(Array3D::<f64>::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::<f64>::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<T: RandomAccessByRef<Item = f64> + Shape>(
+    points: &T,
     degree: usize,
     derivatives: usize,
     data: &mut impl Array3DAccess<f64>,
@@ -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<T: RandomAccessByRef<Item = f64> + Shape>(
+    points: &T,
     degree: usize,
     derivatives: usize,
     data: &mut impl Array3DAccess<f64>,
@@ -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<T: RandomAccessByRef<Item = f64> + Shape>(
+    points: &T,
     degree: usize,
     derivatives: usize,
     data: &mut impl Array3DAccess<f64>,
@@ -368,9 +369,9 @@ pub fn derivative_count(cell_type: ReferenceCellType, derivatives: usize) -> usi
     }
 }
 
-pub fn legendre_shape<'a>(
+pub fn legendre_shape<T: RandomAccessByRef<Item = f64> + 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<T: RandomAccessByRef<Item = f64> + Shape>(
     cell_type: ReferenceCellType,
-    points: &impl Array2DAccess<'a, f64>,
+    points: &T,
     degree: usize,
     derivatives: usize,
     data: &mut impl Array3DAccess<f64>,
@@ -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::<f64>::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::<f64>::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::<f64>::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::<f64>::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::<f64>::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::<f64>::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::<f64>::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::<f64>::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<T: RandomAccessMut<Item = f64> + Shape>(
+        &self,
+        element: &impl FiniteElement,
+        points: &T,
+    ) -> Box<dyn Fn(usize, &mut T)> {
+        let npts = points.shape().0;
+        let mut table = Array4D::<f64>::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::<usize>::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::<usize>()
-        );
-    }
-}
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<T: Num> {
-    /// The data in the array, in row-major order
-    pub data: Vec<T>,
-    /// The shape of the array
-    shape: (usize, usize),
-}
-
-impl<T: Num + Clone> Array2D<T> {
-    /// 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<T> = rlst_dense::Matrix<
+    T,
+    rlst_dense::base_matrix::BaseMatrix<T, rlst_dense::VectorContainer<T>, rlst_dense::Dynamic>,
+    rlst_dense::Dynamic,
+>;
+
+pub fn to_matrix<T: Scalar>(data: &[T], shape: (usize, usize)) -> Mat<T> {
+    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<T>, 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<T> {
-    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::<T> {
-            array: self,
-            index: 0,
-        }
-    }
-}
-
-pub struct Array2DRowIterator<'a, T: Num> {
-    array: &'a Array2D<T>,
-    index: usize,
-}
-
-impl<'a, T: Num> Iterator for Array2DRowIterator<'a, T> {
-    type Item = &'a [T];
-    fn next(&mut self) -> Option<Self::Item> {
-        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::<usize>::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::<f32>::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<T: Num> {
     /// 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<T: RandomAccessByRef<Item = f64> + Shape>(
         &self,
-        points: &impl Array2DAccess<'a, f64>,
+        points: &T,
         nderivs: usize,
         data: &mut impl Array4DAccess<f64>,
     );
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<Item = f64> + Shape,
+        TMut: RandomAccessMut<Item = f64> + 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<Item = f64> + Shape>(
-        &self,
-        points: &impl Array2DAccess<'a, f64>,
-        cell: usize,
-        physical_points: &mut T,
-    );
+        element: &impl FiniteElement,
+        points: &T,
+    ) -> Box<dyn Fn(usize, &mut TMut)>;
 
     /// Compute the physical coordinates of a set of points in a given cell
-    fn compute_points_rlst<'a, T: RandomAccessMut<Item = f64> + Shape>(
+    fn compute_points<
+        'a,
+        T: RandomAccessByRef<Item = f64> + Shape,
+        TMut: RandomAccessMut<Item = f64> + 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<Item = f64> + Shape,
+        TMut: RandomAccessMut<Item = f64> + 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<Item = f64> + Shape,
+        TMut: RandomAccessMut<Item = f64> + 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<T: RandomAccessByRef<Item = f64> + 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<Item = f64> + Shape,
+        TMut: RandomAccessMut<Item = f64> + Shape,
+    >(
         &self,
-        points: &impl Array2DAccess<'a, f64>,
+        points: &T,
         cell: usize,
-        jacobian_inverses: &mut impl Array2DAccess<'a, f64>,
+        jacobian_inverses: &mut TMut,
     );
 }