diff --git a/examples/lagrange_element.rs b/examples/lagrange_element.rs
index 8901146..fed061d 100644
--- a/examples/lagrange_element.rs
+++ b/examples/lagrange_element.rs
@@ -17,9 +17,9 @@ fn main() {
     // Create an array to store the basis function values
     let mut basis_values = rlst_dynamic_array4!(f64, element.tabulate_array_shape(0, 1));
     // Create array containing the point [1/3, 1/3]
-    let mut points = rlst_dynamic_array2!(f64, [1, 2]);
+    let mut points = rlst_dynamic_array2!(f64, [2, 1]);
     points[[0, 0]] = 1.0 / 3.0;
-    points[[0, 1]] = 1.0 / 3.0;
+    points[[1, 0]] = 1.0 / 3.0;
     // Tabulate the element's basis functions at the point
     element.tabulate(&points, 0, &mut basis_values);
     println!(
@@ -29,7 +29,7 @@ fn main() {
 
     // Set point to [1, 0]
     points[[0, 0]] = 1.0;
-    points[[0, 1]] = 0.0;
+    points[[1, 0]] = 0.0;
     // Tabulate the element's basis functions at the point
     element.tabulate(&points, 0, &mut basis_values);
     println!(
diff --git a/src/ciarlet.rs b/src/ciarlet.rs
index bc3e481..066ec52 100644
--- a/src/ciarlet.rs
+++ b/src/ciarlet.rs
@@ -88,22 +88,23 @@ impl<T: RlstScalar + MatrixInverse> CiarletElement<T> {
         }
 
         let new_pts = if continuity == Continuity::Discontinuous {
+            println!("A");
             let mut new_pts: EntityPoints<T::Real> = [vec![], vec![], vec![], vec![]];
-            let mut all_pts = rlst_dynamic_array2![T::Real, [npts, tdim]];
+            let mut all_pts = rlst_dynamic_array2![T::Real, [tdim, npts]];
             for (i, pts_i) in interpolation_points.iter().take(tdim).enumerate() {
                 for _pts in pts_i {
-                    new_pts[i].push(rlst_dynamic_array2![T::Real, [0, tdim]]);
+                    new_pts[i].push(rlst_dynamic_array2![T::Real, [tdim, 0]]);
                 }
             }
-            let mut row = 0;
+            let mut col = 0;
             for pts_i in interpolation_points.iter() {
                 for pts in pts_i {
-                    let nrows = pts.shape()[0];
+                    let ncols = pts.shape()[1];
                     all_pts
                         .view_mut()
-                        .into_subview([row, 0], [nrows, tdim])
+                        .into_subview([0, col], [tdim, ncols])
                         .fill_from(pts.view());
-                    row += nrows;
+                    col += ncols;
                 }
             }
             new_pts[tdim].push(all_pts);
@@ -146,8 +147,8 @@ impl<T: RlstScalar + MatrixInverse> CiarletElement<T> {
         let mut dof = 0;
         for d in 0..4 {
             for (e, pts) in new_pts[d].iter().enumerate() {
-                if pts.shape()[0] > 0 {
-                    let mut table = rlst_dynamic_array3!(T, [1, pdim, pts.shape()[0]]);
+                if pts.shape()[1] > 0 {
+                    let mut table = rlst_dynamic_array3!(T, [1, pdim, pts.shape()[1]]);
                     tabulate_legendre_polynomials(
                         cell_type,
                         pts,
@@ -207,9 +208,9 @@ impl<T: RlstScalar + MatrixInverse> CiarletElement<T> {
         let mut dof = 0;
         for i in 0..4 {
             for pts in &new_pts[i] {
-                let dofs = (dof..dof + pts.shape()[0]).collect::<Vec<_>>();
+                let dofs = (dof..dof + pts.shape()[1]).collect::<Vec<_>>();
                 entity_dofs[i].push(dofs);
-                dof += pts.shape()[0];
+                dof += pts.shape()[1];
             }
         }
         let connectivity = reference_cell::connectivity(cell_type);
@@ -307,7 +308,7 @@ impl<T: RlstScalar + MatrixInverse> FiniteElement for CiarletElement<T> {
         );
 
         for d in 0..table.shape()[0] {
-            for p in 0..points.shape()[0] {
+            for p in 0..points.shape()[1] {
                 for j in 0..self.value_size {
                     for b in 0..self.dim {
                         // data[d, p, b, j] = inner(self.coefficients[b, j, :], table[d, :, p])
@@ -384,11 +385,11 @@ mod test {
         let e = lagrange::create::<f64>(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 mut points = rlst_dynamic_array2!(f64, [4, 1]);
+        let mut points = rlst_dynamic_array2!(f64, [1, 4]);
         *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;
+        *points.get_mut([0, 1]).unwrap() = 0.2;
+        *points.get_mut([0, 2]).unwrap() = 0.4;
+        *points.get_mut([0, 3]).unwrap() = 1.0;
         e.tabulate(&points, 0, &mut data);
 
         for pt in 0..4 {
@@ -402,21 +403,21 @@ mod test {
         let e = lagrange::create::<f64>(ReferenceCellType::Interval, 1, Continuity::Standard);
         assert_eq!(e.value_size(), 1);
         let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 4));
-        let mut points = rlst_dynamic_array2!(f64, [4, 1]);
+        let mut points = rlst_dynamic_array2!(f64, [1, 4]);
         *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;
+        *points.get_mut([0, 1]).unwrap() = 0.2;
+        *points.get_mut([0, 2]).unwrap() = 0.4;
+        *points.get_mut([0, 3]).unwrap() = 1.0;
         e.tabulate(&points, 0, &mut data);
 
         for pt in 0..4 {
             assert_relative_eq!(
                 *data.get([0, pt, 0, 0]).unwrap(),
-                1.0 - *points.get([pt, 0]).unwrap()
+                1.0 - *points.get([0, pt]).unwrap()
             );
             assert_relative_eq!(
                 *data.get([0, pt, 1, 0]).unwrap(),
-                *points.get([pt, 0]).unwrap()
+                *points.get([0, pt]).unwrap()
             );
         }
         check_dofs(e);
@@ -428,19 +429,19 @@ mod test {
         assert_eq!(e.value_size(), 1);
         let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 6));
 
-        let mut points = rlst_dynamic_array2!(f64, [6, 2]);
+        let mut points = rlst_dynamic_array2!(f64, [2, 6]);
         *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, 0]).unwrap() = 0.0;
+        *points.get_mut([0, 1]).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;
+        *points.get_mut([0, 2]).unwrap() = 0.0;
+        *points.get_mut([1, 2]).unwrap() = 1.0;
+        *points.get_mut([0, 3]).unwrap() = 0.5;
+        *points.get_mut([1, 3]).unwrap() = 0.0;
+        *points.get_mut([0, 4]).unwrap() = 0.0;
+        *points.get_mut([1, 4]).unwrap() = 0.5;
+        *points.get_mut([0, 5]).unwrap() = 0.5;
+        *points.get_mut([1, 5]).unwrap() = 0.5;
 
         e.tabulate(&points, 0, &mut data);
 
@@ -455,33 +456,33 @@ mod test {
         let e = lagrange::create::<f64>(ReferenceCellType::Triangle, 1, Continuity::Standard);
         assert_eq!(e.value_size(), 1);
         let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 6));
-        let mut points = rlst_dynamic_array2!(f64, [6, 2]);
+        let mut points = rlst_dynamic_array2!(f64, [2, 6]);
         *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, 0]).unwrap() = 0.0;
+        *points.get_mut([0, 1]).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;
+        *points.get_mut([0, 2]).unwrap() = 0.0;
+        *points.get_mut([1, 2]).unwrap() = 1.0;
+        *points.get_mut([0, 3]).unwrap() = 0.5;
+        *points.get_mut([1, 3]).unwrap() = 0.0;
+        *points.get_mut([0, 4]).unwrap() = 0.0;
+        *points.get_mut([1, 4]).unwrap() = 0.5;
+        *points.get_mut([0, 5]).unwrap() = 0.5;
+        *points.get_mut([1, 5]).unwrap() = 0.5;
         e.tabulate(&points, 0, &mut data);
 
         for pt in 0..6 {
             assert_relative_eq!(
                 *data.get([0, pt, 0, 0]).unwrap(),
-                1.0 - *points.get([pt, 0]).unwrap() - *points.get([pt, 1]).unwrap()
+                1.0 - *points.get([0, pt]).unwrap() - *points.get([1, pt]).unwrap()
             );
             assert_relative_eq!(
                 *data.get([0, pt, 1, 0]).unwrap(),
-                *points.get([pt, 0]).unwrap()
+                *points.get([0, pt]).unwrap()
             );
             assert_relative_eq!(
                 *data.get([0, pt, 2, 0]).unwrap(),
-                *points.get([pt, 1]).unwrap()
+                *points.get([1, pt]).unwrap()
             );
         }
         check_dofs(e);
@@ -551,19 +552,19 @@ mod test {
         );
         assert_eq!(e.value_size(), 1);
         let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 6));
-        let mut points = rlst_dynamic_array2!(f64, [6, 2]);
+        let mut points = rlst_dynamic_array2!(f64, [2, 6]);
         *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, 0]).unwrap() = 0.0;
+        *points.get_mut([0, 1]).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;
+        *points.get_mut([0, 2]).unwrap() = 0.0;
+        *points.get_mut([1, 2]).unwrap() = 1.0;
+        *points.get_mut([0, 3]).unwrap() = 0.5;
+        *points.get_mut([1, 3]).unwrap() = 0.0;
+        *points.get_mut([0, 4]).unwrap() = 0.0;
+        *points.get_mut([1, 4]).unwrap() = 0.5;
+        *points.get_mut([0, 5]).unwrap() = 0.5;
+        *points.get_mut([1, 5]).unwrap() = 0.5;
         e.tabulate(&points, 0, &mut data);
 
         for pt in 0..6 {
@@ -577,38 +578,38 @@ mod test {
         let e = lagrange::create::<f64>(ReferenceCellType::Quadrilateral, 1, Continuity::Standard);
         assert_eq!(e.value_size(), 1);
         let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 6));
-        let mut points = rlst_dynamic_array2!(f64, [6, 2]);
+        let mut points = rlst_dynamic_array2!(f64, [2, 6]);
         *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, 0]).unwrap() = 0.0;
+        *points.get_mut([0, 1]).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;
+        *points.get_mut([0, 2]).unwrap() = 0.0;
+        *points.get_mut([1, 2]).unwrap() = 1.0;
+        *points.get_mut([0, 3]).unwrap() = 1.0;
+        *points.get_mut([1, 3]).unwrap() = 1.0;
+        *points.get_mut([0, 4]).unwrap() = 0.25;
+        *points.get_mut([1, 4]).unwrap() = 0.5;
+        *points.get_mut([0, 5]).unwrap() = 0.3;
+        *points.get_mut([1, 5]).unwrap() = 0.2;
 
         e.tabulate(&points, 0, &mut data);
 
         for pt in 0..6 {
             assert_relative_eq!(
                 *data.get([0, pt, 0, 0]).unwrap(),
-                (1.0 - *points.get([pt, 0]).unwrap()) * (1.0 - *points.get([pt, 1]).unwrap())
+                (1.0 - *points.get([0, pt]).unwrap()) * (1.0 - *points.get([1, pt]).unwrap())
             );
             assert_relative_eq!(
                 *data.get([0, pt, 1, 0]).unwrap(),
-                *points.get([pt, 0]).unwrap() * (1.0 - *points.get([pt, 1]).unwrap())
+                *points.get([0, pt]).unwrap() * (1.0 - *points.get([1, pt]).unwrap())
             );
             assert_relative_eq!(
                 *data.get([0, pt, 2, 0]).unwrap(),
-                (1.0 - *points.get([pt, 0]).unwrap()) * *points.get([pt, 1]).unwrap()
+                (1.0 - *points.get([0, pt]).unwrap()) * *points.get([1, pt]).unwrap()
             );
             assert_relative_eq!(
                 *data.get([0, pt, 3, 0]).unwrap(),
-                *points.get([pt, 0]).unwrap() * *points.get([pt, 1]).unwrap()
+                *points.get([0, pt]).unwrap() * *points.get([1, pt]).unwrap()
             );
         }
         check_dofs(e);
@@ -619,24 +620,25 @@ mod test {
         let e = lagrange::create::<f64>(ReferenceCellType::Quadrilateral, 2, Continuity::Standard);
         assert_eq!(e.value_size(), 1);
         let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 6));
-        let mut points = rlst_dynamic_array2!(f64, [6, 2]);
+        let mut points = rlst_dynamic_array2!(f64, [2, 6]);
         *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, 0]).unwrap() = 0.0;
+        *points.get_mut([0, 1]).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;
+        *points.get_mut([0, 2]).unwrap() = 0.0;
+        *points.get_mut([1, 2]).unwrap() = 1.0;
+        *points.get_mut([0, 3]).unwrap() = 1.0;
+        *points.get_mut([1, 3]).unwrap() = 1.0;
+        *points.get_mut([0, 4]).unwrap() = 0.25;
+        *points.get_mut([1, 4]).unwrap() = 0.5;
+        *points.get_mut([0, 5]).unwrap() = 0.3;
+        *points.get_mut([1, 5]).unwrap() = 0.2;
+
         e.tabulate(&points, 0, &mut data);
 
         for pt in 0..6 {
-            let x = *points.get([pt, 0]).unwrap();
-            let y = *points.get([pt, 1]).unwrap();
+            let x = *points.get([0, pt]).unwrap();
+            let y = *points.get([1, pt]).unwrap();
             assert_relative_eq!(
                 *data.get([0, pt, 0, 0]).unwrap(),
                 (1.0 - x) * (1.0 - 2.0 * x) * (1.0 - y) * (1.0 - 2.0 * y),
@@ -691,45 +693,45 @@ mod test {
         let e = raviart_thomas::create(ReferenceCellType::Triangle, 1, Continuity::Standard);
         assert_eq!(e.value_size(), 2);
         let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 6));
-        let mut points = rlst_dynamic_array2!(f64, [6, 2]);
+        let mut points = rlst_dynamic_array2!(f64, [2, 6]);
         *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, 0]).unwrap() = 0.0;
+        *points.get_mut([0, 1]).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;
+        *points.get_mut([0, 2]).unwrap() = 0.0;
+        *points.get_mut([1, 2]).unwrap() = 1.0;
+        *points.get_mut([0, 3]).unwrap() = 0.5;
+        *points.get_mut([1, 3]).unwrap() = 0.0;
+        *points.get_mut([0, 4]).unwrap() = 0.0;
+        *points.get_mut([1, 4]).unwrap() = 0.5;
+        *points.get_mut([0, 5]).unwrap() = 0.5;
+        *points.get_mut([1, 5]).unwrap() = 0.5;
         e.tabulate(&points, 0, &mut data);
 
         for pt in 0..6 {
             assert_relative_eq!(
                 *data.get([0, pt, 0, 0]).unwrap(),
-                -*points.get([pt, 0]).unwrap()
+                -*points.get([0, pt]).unwrap()
             );
             assert_relative_eq!(
                 *data.get([0, pt, 0, 1]).unwrap(),
-                -*points.get([pt, 1]).unwrap()
+                -*points.get([1, pt]).unwrap()
             );
             assert_relative_eq!(
                 *data.get([0, pt, 1, 0]).unwrap(),
-                *points.get([pt, 0]).unwrap() - 1.0
+                *points.get([0, pt]).unwrap() - 1.0
             );
             assert_relative_eq!(
                 *data.get([0, pt, 1, 1]).unwrap(),
-                *points.get([pt, 1]).unwrap()
+                *points.get([1, pt]).unwrap()
             );
             assert_relative_eq!(
                 *data.get([0, pt, 2, 0]).unwrap(),
-                -*points.get([pt, 0]).unwrap()
+                -*points.get([0, pt]).unwrap()
             );
             assert_relative_eq!(
                 *data.get([0, pt, 2, 1]).unwrap(),
-                1.0 - *points.get([pt, 1]).unwrap()
+                1.0 - *points.get([1, pt]).unwrap()
             );
         }
         check_dofs(e);
diff --git a/src/ciarlet/lagrange.rs b/src/ciarlet/lagrange.rs
index 7869612..6fa3a42 100644
--- a/src/ciarlet/lagrange.rs
+++ b/src/ciarlet/lagrange.rs
@@ -31,17 +31,17 @@ pub fn create<T: RlstScalar + MatrixInverse>(
         }
         for (d, counts) in entity_counts.iter().enumerate() {
             for _e in 0..*counts {
-                x[d].push(rlst_dynamic_array2!(T::Real, [0, tdim]));
+                x[d].push(rlst_dynamic_array2!(T::Real, [tdim, 0]));
                 m[d].push(rlst_dynamic_array3!(T, [0, 1, 0]));
             }
         }
-        let mut midp = rlst_dynamic_array2!(T::Real, [1, tdim]);
+        let mut midp = rlst_dynamic_array2!(T::Real, [tdim, 1]);
         let nvertices = entity_counts[0];
         for i in 0..tdim {
             for vertex in &vertices {
-                *midp.get_mut([0, i]).unwrap() += num::cast::<_, T::Real>(vertex[i]).unwrap();
+                *midp.get_mut([i, 0]).unwrap() += num::cast::<_, T::Real>(vertex[i]).unwrap();
             }
-            *midp.get_mut([0, i]).unwrap() /= num::cast::<_, T::Real>(nvertices).unwrap();
+            *midp.get_mut([i, 0]).unwrap() /= num::cast::<_, T::Real>(nvertices).unwrap();
         }
         x[tdim].push(midp);
         let mut mentry = rlst_dynamic_array3!(T, [1, 1, 1]);
@@ -52,9 +52,9 @@ pub fn create<T: RlstScalar + MatrixInverse>(
         let faces = reference_cell::faces(cell_type);
         // TODO: GLL points
         for vertex in &vertices {
-            let mut pts = rlst_dynamic_array2!(T::Real, [1, tdim]);
+            let mut pts = rlst_dynamic_array2!(T::Real, [tdim, 1]);
             for (i, v) in vertex.iter().enumerate() {
-                *pts.get_mut([0, i]).unwrap() = num::cast::<_, T::Real>(*v).unwrap();
+                *pts.get_mut([i, 0]).unwrap() = num::cast::<_, T::Real>(*v).unwrap();
             }
             x[0].push(pts);
             let mut mentry = rlst_dynamic_array3!(T, [1, 1, 1]);
@@ -62,7 +62,7 @@ pub fn create<T: RlstScalar + MatrixInverse>(
             m[0].push(mentry);
         }
         for e in &edges {
-            let mut pts = rlst_dynamic_array2!(T::Real, [degree - 1, tdim]);
+            let mut pts = rlst_dynamic_array2!(T::Real, [tdim, degree - 1]);
             let [vn0, vn1] = e[..] else {
                 panic!();
             };
@@ -73,7 +73,7 @@ pub fn create<T: RlstScalar + MatrixInverse>(
             for i in 1..degree {
                 *ident.get_mut([i - 1, 0, i - 1]).unwrap() = T::from(1.0).unwrap();
                 for j in 0..tdim {
-                    *pts.get_mut([i - 1, j]).unwrap() = num::cast::<_, T::Real>(v0[j]).unwrap()
+                    *pts.get_mut([j, i - 1]).unwrap() = num::cast::<_, T::Real>(v0[j]).unwrap()
                         + num::cast::<_, T::Real>(i).unwrap()
                             / num::cast::<_, T::Real>(degree).unwrap()
                             * num::cast::<_, T::Real>(v1[j] - v0[j]).unwrap();
@@ -99,7 +99,7 @@ pub fn create<T: RlstScalar + MatrixInverse>(
                     panic!("Unsupported face type");
                 }
             };
-            let mut pts = rlst_dynamic_array2!(T::Real, [npts, tdim]);
+            let mut pts = rlst_dynamic_array2!(T::Real, [tdim, npts]);
 
             let [vn0, vn1, vn2] = faces[e][..3] else {
                 panic!();
@@ -114,7 +114,7 @@ pub fn create<T: RlstScalar + MatrixInverse>(
                     for i0 in 1..degree {
                         for i1 in 1..degree - i0 {
                             for j in 0..tdim {
-                                *pts.get_mut([n, j]).unwrap() = num::cast::<_, T::Real>(v0[j])
+                                *pts.get_mut([j, n]).unwrap() = num::cast::<_, T::Real>(v0[j])
                                     .unwrap()
                                     + num::cast::<_, T::Real>(i0).unwrap()
                                         / num::cast::<_, T::Real>(degree).unwrap()
@@ -132,7 +132,7 @@ pub fn create<T: RlstScalar + MatrixInverse>(
                     for i0 in 1..degree {
                         for i1 in 1..degree {
                             for j in 0..tdim {
-                                *pts.get_mut([n, j]).unwrap() = num::cast::<_, T::Real>(v0[j])
+                                *pts.get_mut([j, n]).unwrap() = num::cast::<_, T::Real>(v0[j])
                                     .unwrap()
                                     + num::cast::<_, T::Real>(i0).unwrap()
                                         / num::cast::<_, T::Real>(degree).unwrap()
diff --git a/src/ciarlet/raviart_thomas.rs b/src/ciarlet/raviart_thomas.rs
index c7c9bed..c3d0a5a 100644
--- a/src/ciarlet/raviart_thomas.rs
+++ b/src/ciarlet/raviart_thomas.rs
@@ -52,12 +52,12 @@ pub fn create<T: RlstScalar + MatrixInverse>(
     let edges = reference_cell::edges(cell_type);
 
     for _e in 0..entity_counts[0] {
-        x[0].push(rlst_dynamic_array2!(T::Real, [0, tdim]));
+        x[0].push(rlst_dynamic_array2!(T::Real, [tdim, 0]));
         m[0].push(rlst_dynamic_array3!(T, [0, 2, 0]));
     }
 
     for e in &edges {
-        let mut pts = rlst_dynamic_array2!(T::Real, [1, tdim]);
+        let mut pts = rlst_dynamic_array2!(T::Real, [tdim, 1]);
         let mut mat = rlst_dynamic_array3!(T, [1, 2, 1]);
         let [vn0, vn1] = e[..] else {
             panic!();
@@ -65,7 +65,7 @@ pub fn create<T: RlstScalar + MatrixInverse>(
         let v0 = &vertices[vn0];
         let v1 = &vertices[vn1];
         for i in 0..tdim {
-            *pts.get_mut([0, i]).unwrap() = num::cast::<_, T::Real>(v0[i] + v1[i]).unwrap()
+            *pts.get_mut([i, 0]).unwrap() = num::cast::<_, T::Real>(v0[i] + v1[i]).unwrap()
                 / num::cast::<_, T::Real>(2.0).unwrap();
         }
         *mat.get_mut([0, 0, 0]).unwrap() = T::from(v0[1] - v1[1]).unwrap();
@@ -75,7 +75,7 @@ pub fn create<T: RlstScalar + MatrixInverse>(
     }
 
     for _e in 0..entity_counts[2] {
-        x[2].push(rlst_dynamic_array2!(T::Real, [0, tdim]));
+        x[2].push(rlst_dynamic_array2!(T::Real, [tdim, 0]));
         m[2].push(rlst_dynamic_array3!(T, [0, 2, 0]))
     }
 
diff --git a/src/polynomials.rs b/src/polynomials.rs
index d872709..0cd8df6 100644
--- a/src/polynomials.rs
+++ b/src/polynomials.rs
@@ -16,8 +16,8 @@ fn tabulate_legendre_polynomials_interval<
 ) {
     assert_eq!(data.shape()[0], derivatives + 1);
     assert_eq!(data.shape()[1], degree + 1);
-    assert_eq!(data.shape()[2], points.shape()[0]);
-    assert_eq!(points.shape()[1], 1);
+    assert_eq!(data.shape()[2], points.shape()[1]);
+    assert_eq!(points.shape()[0], 1);
 
     for i in 0..data.shape()[2] {
         *data.get_mut([0, 0, i]).unwrap() = T::from(1.0).unwrap();
@@ -38,7 +38,7 @@ fn tabulate_legendre_polynomials_interval<
             for i in 0..data.shape()[2] {
                 let d = *data.get([k, p - 1, i]).unwrap();
                 *data.get_mut([k, p, i]).unwrap() =
-                    (T::from(*points.get([i, 0]).unwrap()).unwrap() * T::from(2.0).unwrap()
+                    (T::from(*points.get([0, i]).unwrap()).unwrap() * T::from(2.0).unwrap()
                         - T::from(1.0).unwrap())
                         * d
                         * b;
@@ -85,8 +85,8 @@ fn tabulate_legendre_polynomials_quadrilateral<
 ) {
     assert_eq!(data.shape()[0], (derivatives + 1) * (derivatives + 2) / 2);
     assert_eq!(data.shape()[1], (degree + 1) * (degree + 1));
-    assert_eq!(data.shape()[2], points.shape()[0]);
-    assert_eq!(points.shape()[1], 2);
+    assert_eq!(data.shape()[2], points.shape()[1]);
+    assert_eq!(points.shape()[0], 2);
 
     for i in 0..data.shape()[2] {
         *data
@@ -116,7 +116,7 @@ fn tabulate_legendre_polynomials_quadrilateral<
                     .unwrap();
                 *data
                     .get_mut([tri_index(k, 0), quad_index(p, 0, degree), i])
-                    .unwrap() = (T::from(*points.get([i, 0]).unwrap()).unwrap()
+                    .unwrap() = (T::from(*points.get([0, i]).unwrap()).unwrap()
                     * T::from(2.0).unwrap()
                     - T::from(1.0).unwrap())
                     * d
@@ -171,7 +171,7 @@ fn tabulate_legendre_polynomials_quadrilateral<
                     .unwrap();
                 *data
                     .get_mut([tri_index(0, k), quad_index(0, p, degree), i])
-                    .unwrap() = (T::from(*points.get([i, 1]).unwrap()).unwrap()
+                    .unwrap() = (T::from(*points.get([1, i]).unwrap()).unwrap()
                     * T::from(2.0).unwrap()
                     - T::from(1.0).unwrap())
                     * d
@@ -238,8 +238,8 @@ fn tabulate_legendre_polynomials_triangle<
 ) {
     assert_eq!(data.shape()[0], (derivatives + 1) * (derivatives + 2) / 2);
     assert_eq!(data.shape()[1], (degree + 1) * (degree + 2) / 2);
-    assert_eq!(data.shape()[2], points.shape()[0]);
-    assert_eq!(points.shape()[1], 2);
+    assert_eq!(data.shape()[2], points.shape()[1]);
+    assert_eq!(points.shape()[0], 2);
 
     for i in 0..data.shape()[2] {
         *data.get_mut([tri_index(0, 0), tri_index(0, 0), i]).unwrap() =
@@ -267,9 +267,9 @@ fn tabulate_legendre_polynomials_triangle<
                         .unwrap();
                     *data
                         .get_mut([tri_index(kx, ky), tri_index(0, p), i])
-                        .unwrap() = (T::from(*points.get([i, 0]).unwrap()).unwrap()
+                        .unwrap() = (T::from(*points.get([0, i]).unwrap()).unwrap()
                         * T::from(2.0).unwrap()
-                        + T::from(*points.get([i, 1]).unwrap()).unwrap()
+                        + T::from(*points.get([1, i]).unwrap()).unwrap()
                         - T::from(1.0).unwrap())
                         * d
                         * a
@@ -307,7 +307,7 @@ fn tabulate_legendre_polynomials_triangle<
 
                     for i in 0..data.shape()[2] {
                         let b =
-                            T::from(1.0).unwrap() - T::from(*points.get([i, 1]).unwrap()).unwrap();
+                            T::from(1.0).unwrap() - T::from(*points.get([1, i]).unwrap()).unwrap();
                         let d = *data
                             .get([tri_index(kx, ky), tri_index(0, p - 2), i])
                             .unwrap();
@@ -324,7 +324,7 @@ fn tabulate_legendre_polynomials_triangle<
                                 .get_mut([tri_index(kx, ky), tri_index(0, p), i])
                                 .unwrap() -= T::from(2.0).unwrap()
                                 * T::from(ky).unwrap()
-                                * (T::from(*points.get([i, 1]).unwrap()).unwrap()
+                                * (T::from(*points.get([1, i]).unwrap()).unwrap()
                                     - T::from(1.0).unwrap())
                                 * d
                                 * scale2
@@ -357,7 +357,7 @@ fn tabulate_legendre_polynomials_triangle<
                         .get_mut([tri_index(kx, ky), tri_index(1, p), i])
                         .unwrap() = *data.get([tri_index(kx, ky), tri_index(0, p), i]).unwrap()
                         * scale3
-                        * ((T::from(*points.get([i, 1]).unwrap()).unwrap()
+                        * ((T::from(*points.get([1, i]).unwrap()).unwrap()
                             * T::from(2.0).unwrap()
                             - T::from(1.0).unwrap())
                             * (T::from(1.5).unwrap() + T::from(p).unwrap())
@@ -400,7 +400,7 @@ fn tabulate_legendre_polynomials_triangle<
                             .get_mut([tri_index(kx, ky), tri_index(q + 1, p), i])
                             .unwrap() = d
                             * scale4
-                            * ((T::from(*points.get([i, 1]).unwrap()).unwrap()
+                            * ((T::from(*points.get([1, i]).unwrap()).unwrap()
                                 * T::from(T::from(2.0).unwrap()).unwrap()
                                 - T::from(T::from(1.0).unwrap()).unwrap())
                                 * a1
@@ -465,7 +465,7 @@ pub fn legendre_shape<T, Array2: RandomAccessByRef<2, Item = T> + Shape<2>>(
     [
         derivative_count(cell_type, derivatives),
         polynomial_count(cell_type, degree),
-        points.shape()[0],
+        points.shape()[1],
     ]
 }
 
@@ -514,9 +514,9 @@ mod test {
         )
         .unwrap();
 
-        let mut points = rlst_dynamic_array2!(f64, [rule.npoints, 1]);
+        let mut points = rlst_dynamic_array2!(f64, [1, rule.npoints]);
         for (i, j) in rule.points.iter().enumerate() {
-            *points.get_mut([i, 0]).unwrap() = *j;
+            *points.get_mut([0, i]).unwrap() = *j;
         }
 
         let mut data = rlst_dynamic_array3!(
@@ -547,10 +547,10 @@ mod test {
         let degree = 5;
 
         let rule = simplex_rule(bempp::traits::types::ReferenceCellType::Triangle, 79).unwrap();
-        let mut points = rlst_dynamic_array2!(f64, [rule.npoints, 2]);
+        let mut points = rlst_dynamic_array2!(f64, [2, rule.npoints]);
         for i in 0..rule.npoints {
             for j in 0..2 {
-                *points.get_mut([i, j]).unwrap() = rule.points[i * 2 + j];
+                *points.get_mut([j, i]).unwrap() = rule.points[i * 2 + j];
             }
         }
 
@@ -583,10 +583,10 @@ mod test {
 
         let rule =
             simplex_rule(bempp::traits::types::ReferenceCellType::Quadrilateral, 85).unwrap();
-        let mut points = rlst_dynamic_array2!(f64, [rule.npoints, 2]);
+        let mut points = rlst_dynamic_array2!(f64, [2, rule.npoints]);
         for i in 0..rule.npoints {
             for j in 0..2 {
-                *points.get_mut([i, j]).unwrap() = rule.points[i * 2 + j];
+                *points.get_mut([j, i]).unwrap() = rule.points[i * 2 + j];
             }
         }
 
@@ -624,10 +624,10 @@ mod test {
         let degree = 6;
 
         let epsilon = 1e-10;
-        let mut points = rlst_dynamic_array2!(f64, [20, 1]);
+        let mut points = rlst_dynamic_array2!(f64, [1, 20]);
         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;
+            *points.get_mut([0, 2 * i]).unwrap() = i as f64 / 10.0;
+            *points.get_mut([0, 2 * i + 1]).unwrap() = points.get([0, 2 * i]).unwrap() + epsilon;
         }
 
         let mut data = rlst_dynamic_array3!(
@@ -637,7 +637,7 @@ mod test {
         tabulate_legendre_polynomials(ReferenceCellType::Interval, &points, degree, 1, &mut data);
 
         for i in 0..degree + 1 {
-            for k in 0..points.shape()[0] / 2 {
+            for k in 0..points.shape()[1] / 2 {
                 assert_relative_eq!(
                     *data.get([1, i, 2 * k]).unwrap(),
                     (data.get([0, i, 2 * k + 1]).unwrap() - data.get([0, i, 2 * k]).unwrap())
@@ -653,18 +653,18 @@ mod test {
         let degree = 6;
 
         let epsilon = 1e-10;
-        let mut points = rlst_dynamic_array2!(f64, [165, 2]);
+        let mut points = rlst_dynamic_array2!(f64, [2, 165]);
         let mut index = 0;
         for i in 0..10 {
             for j in 0..10 - i {
-                *points.get_mut([3 * index, 0]).unwrap() = i as f64 / 10.0;
-                *points.get_mut([3 * index, 1]).unwrap() = j as f64 / 10.0;
-                *points.get_mut([3 * index + 1, 0]).unwrap() =
-                    *points.get([3 * index, 0]).unwrap() + epsilon;
-                *points.get_mut([3 * index + 1, 1]).unwrap() = *points.get([3 * index, 1]).unwrap();
-                *points.get_mut([3 * index + 2, 0]).unwrap() = *points.get([3 * index, 0]).unwrap();
-                *points.get_mut([3 * index + 2, 1]).unwrap() =
-                    *points.get([3 * index, 1]).unwrap() + epsilon;
+                *points.get_mut([0, 3 * index]).unwrap() = i as f64 / 10.0;
+                *points.get_mut([1, 3 * index]).unwrap() = j as f64 / 10.0;
+                *points.get_mut([0, 3 * index + 1]).unwrap() =
+                    *points.get([0, 3 * index]).unwrap() + epsilon;
+                *points.get_mut([1, 3 * index + 1]).unwrap() = *points.get([1, 3 * index]).unwrap();
+                *points.get_mut([0, 3 * index + 2]).unwrap() = *points.get([0, 3 * index]).unwrap();
+                *points.get_mut([1, 3 * index + 2]).unwrap() =
+                    *points.get([1, 3 * index]).unwrap() + epsilon;
                 index += 1;
             }
         }
@@ -676,7 +676,7 @@ mod test {
         tabulate_legendre_polynomials(ReferenceCellType::Triangle, &points, degree, 1, &mut data);
 
         for i in 0..degree + 1 {
-            for k in 0..points.shape()[0] / 3 {
+            for k in 0..points.shape()[1] / 3 {
                 assert_relative_eq!(
                     *data.get([1, i, 3 * k]).unwrap(),
                     (data.get([0, i, 3 * k + 1]).unwrap() - data.get([0, i, 3 * k]).unwrap())
@@ -698,18 +698,18 @@ mod test {
         let degree = 6;
 
         let epsilon = 1e-10;
-        let mut points = rlst_dynamic_array2!(f64, [300, 2]);
+        let mut points = rlst_dynamic_array2!(f64, [2, 300]);
         for i in 0..10 {
             for j in 0..10 {
                 let index = 10 * i + j;
-                *points.get_mut([3 * index, 0]).unwrap() = i as f64 / 10.0;
-                *points.get_mut([3 * index, 1]).unwrap() = j as f64 / 10.0;
-                *points.get_mut([3 * index + 1, 0]).unwrap() =
-                    *points.get([3 * index, 0]).unwrap() + epsilon;
-                *points.get_mut([3 * index + 1, 1]).unwrap() = *points.get([3 * index, 1]).unwrap();
-                *points.get_mut([3 * index + 2, 0]).unwrap() = *points.get([3 * index, 0]).unwrap();
-                *points.get_mut([3 * index + 2, 1]).unwrap() =
-                    *points.get([3 * index, 1]).unwrap() + epsilon;
+                *points.get_mut([0, 3 * index]).unwrap() = i as f64 / 10.0;
+                *points.get_mut([1, 3 * index]).unwrap() = j as f64 / 10.0;
+                *points.get_mut([0, 3 * index + 1]).unwrap() =
+                    *points.get([0, 3 * index]).unwrap() + epsilon;
+                *points.get_mut([1, 3 * index + 1]).unwrap() = *points.get([1, 3 * index]).unwrap();
+                *points.get_mut([0, 3 * index + 2]).unwrap() = *points.get([0, 3 * index]).unwrap();
+                *points.get_mut([1, 3 * index + 2]).unwrap() =
+                    *points.get([1, 3 * index]).unwrap() + epsilon;
             }
         }
 
@@ -726,7 +726,7 @@ mod test {
         );
 
         for i in 0..degree + 1 {
-            for k in 0..points.shape()[0] / 3 {
+            for k in 0..points.shape()[1] / 3 {
                 assert_relative_eq!(
                     *data.get([1, i, 3 * k]).unwrap(),
                     (data.get([0, i, 3 * k + 1]).unwrap() - data.get([0, i, 3 * k]).unwrap())
@@ -747,9 +747,9 @@ mod test {
     fn test_legendre_interval_against_known_polynomials() {
         let degree = 3;
 
-        let mut points = rlst_dynamic_array2!(f64, [11, 1]);
+        let mut points = rlst_dynamic_array2!(f64, [1, 11]);
         for i in 0..11 {
-            *points.get_mut([i, 0]).unwrap() = i as f64 / 10.0;
+            *points.get_mut([0, i]).unwrap() = i as f64 / 10.0;
         }
 
         let mut data = rlst_dynamic_array3!(
@@ -759,7 +759,7 @@ mod test {
         tabulate_legendre_polynomials(ReferenceCellType::Interval, &points, degree, 3, &mut data);
 
         for k in 0..points.shape()[0] {
-            let x = *points.get([k, 0]).unwrap();
+            let x = *points.get([0, k]).unwrap();
 
             // 0 => 1
             assert_relative_eq!(*data.get([0, 0, k]).unwrap(), 1.0, epsilon = 1e-12);
@@ -827,11 +827,11 @@ mod test {
     fn test_legendre_quadrilateral_against_known_polynomials() {
         let degree = 2;
 
-        let mut points = rlst_dynamic_array2!(f64, [121, 2]);
+        let mut points = rlst_dynamic_array2!(f64, [2, 121]);
         for i in 0..11 {
             for j in 0..11 {
-                *points.get_mut([11 * i + j, 0]).unwrap() = i as f64 / 10.0;
-                *points.get_mut([11 * i + j, 1]).unwrap() = j as f64 / 10.0;
+                *points.get_mut([0, 11 * i + j]).unwrap() = i as f64 / 10.0;
+                *points.get_mut([1, 11 * i + j]).unwrap() = j as f64 / 10.0;
             }
         }
 
@@ -847,9 +847,9 @@ mod test {
             &mut data,
         );
 
-        for k in 0..points.shape()[0] {
-            let x = *points.get([k, 0]).unwrap();
-            let y = *points.get([k, 1]).unwrap();
+        for k in 0..points.shape()[1] {
+            let x = *points.get([0, k]).unwrap();
+            let y = *points.get([1, k]).unwrap();
 
             // 0 => 1
             assert_relative_eq!(*data.get([0, 0, k]).unwrap(), 1.0, epsilon = 1e-12);