diff --git a/src/ciarlet.rs b/src/ciarlet.rs
index bc1332e..5fe82a2 100644
--- a/src/ciarlet.rs
+++ b/src/ciarlet.rs
@@ -427,14 +427,9 @@ mod test {
         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([0, pt]).unwrap()
-            );
-            assert_relative_eq!(
-                *data.get([0, pt, 1, 0]).unwrap(),
-                *points.get([0, pt]).unwrap()
-            );
+            let x = *points.get([0, pt]).unwrap();
+            assert_relative_eq!(*data.get([0, pt, 0, 0]).unwrap(), 1.0 - x);
+            assert_relative_eq!(*data.get([0, pt, 1, 0]).unwrap(), x);
         }
         check_dofs(e);
     }
@@ -488,18 +483,11 @@ mod test {
         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([0, pt]).unwrap() - *points.get([1, pt]).unwrap()
-            );
-            assert_relative_eq!(
-                *data.get([0, pt, 1, 0]).unwrap(),
-                *points.get([0, pt]).unwrap()
-            );
-            assert_relative_eq!(
-                *data.get([0, pt, 2, 0]).unwrap(),
-                *points.get([1, pt]).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 - y);
+            assert_relative_eq!(*data.get([0, pt, 1, 0]).unwrap(), x);
+            assert_relative_eq!(*data.get([0, pt, 2, 0]).unwrap(), y);
         }
         check_dofs(e);
     }
@@ -556,22 +544,12 @@ mod test {
         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([0, pt]).unwrap()) * (1.0 - *points.get([1, pt]).unwrap())
-            );
-            assert_relative_eq!(
-                *data.get([0, pt, 1, 0]).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([0, pt]).unwrap()) * *points.get([1, pt]).unwrap()
-            );
-            assert_relative_eq!(
-                *data.get([0, pt, 3, 0]).unwrap(),
-                *points.get([0, pt]).unwrap() * *points.get([1, pt]).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 - y));
+            assert_relative_eq!(*data.get([0, pt, 1, 0]).unwrap(), x * (1.0 - y));
+            assert_relative_eq!(*data.get([0, pt, 2, 0]).unwrap(), (1.0 - x) * y);
+            assert_relative_eq!(*data.get([0, pt, 3, 0]).unwrap(), x * y);
         }
         check_dofs(e);
     }
@@ -710,28 +688,17 @@ mod test {
         e.tabulate(&points, 0, &mut data);
 
         for pt in 0..6 {
+            let x = *points.get([0, pt]).unwrap();
+            let y = *points.get([1, pt]).unwrap();
+            let z = *points.get([2, pt]).unwrap();
             assert_relative_eq!(
                 *data.get([0, pt, 0, 0]).unwrap(),
-                1.0 - *points.get([0, pt]).unwrap()
-                    - *points.get([1, pt]).unwrap()
-                    - *points.get([2, pt]).unwrap(),
-                epsilon = 1e-14
-            );
-            assert_relative_eq!(
-                *data.get([0, pt, 1, 0]).unwrap(),
-                *points.get([0, pt]).unwrap(),
-                epsilon = 1e-14
-            );
-            assert_relative_eq!(
-                *data.get([0, pt, 2, 0]).unwrap(),
-                *points.get([1, pt]).unwrap(),
-                epsilon = 1e-14
-            );
-            assert_relative_eq!(
-                *data.get([0, pt, 3, 0]).unwrap(),
-                *points.get([2, pt]).unwrap(),
+                1.0 - x - y - z,
                 epsilon = 1e-14
             );
+            assert_relative_eq!(*data.get([0, pt, 1, 0]).unwrap(), x, epsilon = 1e-14);
+            assert_relative_eq!(*data.get([0, pt, 2, 0]).unwrap(), y, epsilon = 1e-14);
+            assert_relative_eq!(*data.get([0, pt, 3, 0]).unwrap(), z, epsilon = 1e-14);
         }
         check_dofs(e);
     }
@@ -797,60 +764,47 @@ mod test {
         e.tabulate(&points, 0, &mut data);
 
         for pt in 0..6 {
+            let x = *points.get([0, pt]).unwrap();
+            let y = *points.get([1, pt]).unwrap();
+            let z = *points.get([2, pt]).unwrap();
             assert_relative_eq!(
                 *data.get([0, pt, 0, 0]).unwrap(),
-                (1.0 - *points.get([0, pt]).unwrap())
-                    * (1.0 - *points.get([1, pt]).unwrap())
-                    * (1.0 - *points.get([2, pt]).unwrap()),
+                (1.0 - x) * (1.0 - y) * (1.0 - z),
                 epsilon = 1e-14
             );
             assert_relative_eq!(
                 *data.get([0, pt, 1, 0]).unwrap(),
-                *points.get([0, pt]).unwrap()
-                    * (1.0 - *points.get([1, pt]).unwrap())
-                    * (1.0 - *points.get([2, pt]).unwrap()),
+                x * (1.0 - y) * (1.0 - z),
                 epsilon = 1e-14
             );
             assert_relative_eq!(
                 *data.get([0, pt, 2, 0]).unwrap(),
-                (1.0 - *points.get([0, pt]).unwrap())
-                    * *points.get([1, pt]).unwrap()
-                    * (1.0 - *points.get([2, pt]).unwrap()),
+                (1.0 - x) * y * (1.0 - z),
                 epsilon = 1e-14
             );
             assert_relative_eq!(
                 *data.get([0, pt, 3, 0]).unwrap(),
-                *points.get([0, pt]).unwrap()
-                    * *points.get([1, pt]).unwrap()
-                    * (1.0 - *points.get([2, pt]).unwrap()),
+                x * y * (1.0 - z),
                 epsilon = 1e-14
             );
             assert_relative_eq!(
                 *data.get([0, pt, 4, 0]).unwrap(),
-                (1.0 - *points.get([0, pt]).unwrap())
-                    * (1.0 - *points.get([1, pt]).unwrap())
-                    * *points.get([2, pt]).unwrap(),
+                (1.0 - x) * (1.0 - y) * z,
                 epsilon = 1e-14
             );
             assert_relative_eq!(
                 *data.get([0, pt, 5, 0]).unwrap(),
-                *points.get([0, pt]).unwrap()
-                    * (1.0 - *points.get([1, pt]).unwrap())
-                    * *points.get([2, pt]).unwrap(),
+                x * (1.0 - y) * z,
                 epsilon = 1e-14
             );
             assert_relative_eq!(
                 *data.get([0, pt, 6, 0]).unwrap(),
-                (1.0 - *points.get([0, pt]).unwrap())
-                    * *points.get([1, pt]).unwrap()
-                    * *points.get([2, pt]).unwrap(),
+                (1.0 - x) * y * z,
                 epsilon = 1e-14
             );
             assert_relative_eq!(
                 *data.get([0, pt, 7, 0]).unwrap(),
-                *points.get([0, pt]).unwrap()
-                    * *points.get([1, pt]).unwrap()
-                    * *points.get([2, pt]).unwrap(),
+                x * y * z,
                 epsilon = 1e-14
             );
         }
@@ -933,36 +887,13 @@ mod test {
         e.tabulate(&points, 0, &mut data);
 
         for pt in 0..6 {
-            assert_relative_eq!(
-                *data.get([0, pt, 0, 0]).unwrap(),
-                -*points.get([0, pt]).unwrap(),
-                epsilon = 1e-14
-            );
-            assert_relative_eq!(
-                *data.get([0, pt, 0, 1]).unwrap(),
-                -*points.get([1, pt]).unwrap(),
-                epsilon = 1e-14
-            );
-            assert_relative_eq!(
-                *data.get([0, pt, 1, 0]).unwrap(),
-                *points.get([0, pt]).unwrap() - 1.0,
-                epsilon = 1e-14
-            );
-            assert_relative_eq!(
-                *data.get([0, pt, 1, 1]).unwrap(),
-                *points.get([1, pt]).unwrap(),
-                epsilon = 1e-14
-            );
-            assert_relative_eq!(
-                *data.get([0, pt, 2, 0]).unwrap(),
-                -*points.get([0, pt]).unwrap(),
-                epsilon = 1e-14
-            );
-            assert_relative_eq!(
-                *data.get([0, pt, 2, 1]).unwrap(),
-                1.0 - *points.get([1, pt]).unwrap(),
-                epsilon = 1e-14
-            );
+            let x = *points.get([0, pt]).unwrap();
+            let y = *points.get([1, pt]).unwrap();
+            for (i, basis_f) in [[-x, -y], [x - 1.0, y], [-x, 1.0 - y]].iter().enumerate() {
+                for (d, value) in basis_f.iter().enumerate() {
+                    assert_relative_eq!(*data.get([0, pt, i, d]).unwrap(), value, epsilon = 1e-14);
+                }
+            }
         }
         check_dofs(e);
     }
@@ -981,6 +912,63 @@ mod test {
         check_dofs(e);
     }
 
+    #[test]
+    fn test_raviart_thomas_1_quadrilateral() {
+        let e = raviart_thomas::create(ReferenceCellType::Quadrilateral, 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, [2, 6]);
+        *points.get_mut([0, 0]).unwrap() = 0.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([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() = 1.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 {
+            let x = *points.get([0, pt]).unwrap();
+            let y = *points.get([1, pt]).unwrap();
+            for (i, basis_f) in [[0.0, 1.0 - y], [x - 1.0, 0.0], [-x, 0.0], [0.0, y]]
+                .iter()
+                .enumerate()
+            {
+                for (d, value) in basis_f.iter().enumerate() {
+                    assert_relative_eq!(*data.get([0, pt, i, d]).unwrap(), value, epsilon = 1e-14);
+                }
+            }
+        }
+        check_dofs(e);
+    }
+
+    #[test]
+    fn test_raviart_thomas_2_quadrilateral() {
+        let e = raviart_thomas::create::<f64>(
+            ReferenceCellType::Quadrilateral,
+            2,
+            Continuity::Standard,
+        );
+        assert_eq!(e.value_size(), 2);
+        check_dofs(e);
+    }
+
+    #[test]
+    fn test_raviart_thomas_3_quadrilateral() {
+        let e = raviart_thomas::create::<f64>(
+            ReferenceCellType::Quadrilateral,
+            3,
+            Continuity::Standard,
+        );
+        assert_eq!(e.value_size(), 2);
+        check_dofs(e);
+    }
+
     #[test]
     fn test_raviart_thomas_1_tetrahedron() {
         let e =
@@ -1005,6 +993,63 @@ mod test {
         check_dofs(e);
     }
 
+    #[test]
+    fn test_raviart_thomas_1_hexahedron() {
+        let e = raviart_thomas::create(ReferenceCellType::Hexahedron, 1, Continuity::Standard);
+        assert_eq!(e.value_size(), 3);
+        let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 6));
+        let mut points = rlst_dynamic_array2!(f64, [3, 6]);
+        *points.get_mut([0, 0]).unwrap() = 0.0;
+        *points.get_mut([1, 0]).unwrap() = 0.0;
+        *points.get_mut([2, 0]).unwrap() = 0.0;
+        *points.get_mut([0, 1]).unwrap() = 1.0;
+        *points.get_mut([1, 1]).unwrap() = 0.0;
+        *points.get_mut([2, 1]).unwrap() = 0.8;
+        *points.get_mut([0, 2]).unwrap() = 0.0;
+        *points.get_mut([1, 2]).unwrap() = 1.0;
+        *points.get_mut([2, 2]).unwrap() = 1.0;
+        *points.get_mut([0, 3]).unwrap() = 0.5;
+        *points.get_mut([1, 3]).unwrap() = 0.0;
+        *points.get_mut([2, 3]).unwrap() = 0.1;
+        *points.get_mut([0, 4]).unwrap() = 1.0;
+        *points.get_mut([1, 4]).unwrap() = 0.5;
+        *points.get_mut([2, 4]).unwrap() = 0.5;
+        *points.get_mut([0, 5]).unwrap() = 0.5;
+        *points.get_mut([1, 5]).unwrap() = 0.5;
+        *points.get_mut([2, 5]).unwrap() = 1.0;
+        e.tabulate(&points, 0, &mut data);
+
+        for pt in 0..6 {
+            let x = *points.get([0, pt]).unwrap();
+            let y = *points.get([1, pt]).unwrap();
+            let z = *points.get([2, pt]).unwrap();
+            for (i, basis_f) in [
+                [0.0, 0.0, 1.0 - z],
+                [0.0, y - 1.0, 0.0],
+                [1.0 - x, 0.0, 0.0],
+                [x, 0.0, 0.0],
+                [0.0, -y, 0.0],
+                [0.0, 0.0, z],
+            ]
+            .iter()
+            .enumerate()
+            {
+                for (d, value) in basis_f.iter().enumerate() {
+                    assert_relative_eq!(*data.get([0, pt, i, d]).unwrap(), value, epsilon = 1e-14);
+                }
+            }
+        }
+        check_dofs(e);
+    }
+
+    #[test]
+    fn test_raviart_thomas_2_hexahedron() {
+        let e =
+            raviart_thomas::create::<f64>(ReferenceCellType::Hexahedron, 2, Continuity::Standard);
+        assert_eq!(e.value_size(), 3);
+        check_dofs(e);
+    }
+
     #[test]
     fn test_nedelec_1_triangle() {
         let e = nedelec::create(ReferenceCellType::Triangle, 1, Continuity::Standard);
@@ -1026,36 +1071,13 @@ mod test {
         e.tabulate(&points, 0, &mut data);
 
         for pt in 0..6 {
-            assert_relative_eq!(
-                *data.get([0, pt, 0, 0]).unwrap(),
-                -*points.get([1, pt]).unwrap(),
-                epsilon = 1e-14
-            );
-            assert_relative_eq!(
-                *data.get([0, pt, 0, 1]).unwrap(),
-                *points.get([0, pt]).unwrap(),
-                epsilon = 1e-14
-            );
-            assert_relative_eq!(
-                *data.get([0, pt, 1, 0]).unwrap(),
-                *points.get([1, pt]).unwrap(),
-                epsilon = 1e-14
-            );
-            assert_relative_eq!(
-                *data.get([0, pt, 1, 1]).unwrap(),
-                1.0 - *points.get([0, pt]).unwrap(),
-                epsilon = 1e-14
-            );
-            assert_relative_eq!(
-                *data.get([0, pt, 2, 0]).unwrap(),
-                1.0 - *points.get([1, pt]).unwrap(),
-                epsilon = 1e-14
-            );
-            assert_relative_eq!(
-                *data.get([0, pt, 2, 1]).unwrap(),
-                *points.get([0, pt]).unwrap(),
-                epsilon = 1e-14
-            );
+            let x = *points.get([0, pt]).unwrap();
+            let y = *points.get([1, pt]).unwrap();
+            for (i, basis_f) in [[-y, x], [y, 1.0 - x], [1.0 - y, x]].iter().enumerate() {
+                for (d, value) in basis_f.iter().enumerate() {
+                    assert_relative_eq!(*data.get([0, pt, i, d]).unwrap(), value, epsilon = 1e-14);
+                }
+            }
         }
         check_dofs(e);
     }
@@ -1074,6 +1096,55 @@ mod test {
         check_dofs(e);
     }
 
+    #[test]
+    fn test_nedelec_1_quadrilateral() {
+        let e = nedelec::create(ReferenceCellType::Quadrilateral, 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, [2, 6]);
+        *points.get_mut([0, 0]).unwrap() = 0.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([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() = 1.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 {
+            let x = *points.get([0, pt]).unwrap();
+            let y = *points.get([1, pt]).unwrap();
+            for (i, basis_f) in [[1.0 - y, 0.0], [0.0, 1.0 - x], [0.0, x], [y, 0.0]]
+                .iter()
+                .enumerate()
+            {
+                for (d, value) in basis_f.iter().enumerate() {
+                    assert_relative_eq!(*data.get([0, pt, i, d]).unwrap(), value, epsilon = 1e-14);
+                }
+            }
+        }
+        check_dofs(e);
+    }
+
+    #[test]
+    fn test_nedelec_2_quadrilateral() {
+        let e = nedelec::create::<f64>(ReferenceCellType::Quadrilateral, 2, Continuity::Standard);
+        assert_eq!(e.value_size(), 2);
+        check_dofs(e);
+    }
+
+    #[test]
+    fn test_nedelec_3_quadrilateral() {
+        let e = nedelec::create::<f64>(ReferenceCellType::Quadrilateral, 3, Continuity::Standard);
+        assert_eq!(e.value_size(), 2);
+        check_dofs(e);
+    }
+
     #[test]
     fn test_nedelec_1_tetrahedron() {
         let e = nedelec::create::<f64>(ReferenceCellType::Tetrahedron, 1, Continuity::Standard);
@@ -1095,6 +1166,68 @@ mod test {
         check_dofs(e);
     }
 
+    #[test]
+    fn test_nedelec_1_hexahedron() {
+        let e = nedelec::create(ReferenceCellType::Hexahedron, 1, Continuity::Standard);
+        assert_eq!(e.value_size(), 3);
+        let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 6));
+        let mut points = rlst_dynamic_array2!(f64, [3, 6]);
+        *points.get_mut([0, 0]).unwrap() = 0.0;
+        *points.get_mut([1, 0]).unwrap() = 0.0;
+        *points.get_mut([2, 0]).unwrap() = 0.0;
+        *points.get_mut([0, 1]).unwrap() = 1.0;
+        *points.get_mut([1, 1]).unwrap() = 0.0;
+        *points.get_mut([2, 1]).unwrap() = 0.8;
+        *points.get_mut([0, 2]).unwrap() = 0.0;
+        *points.get_mut([1, 2]).unwrap() = 1.0;
+        *points.get_mut([2, 2]).unwrap() = 1.0;
+        *points.get_mut([0, 3]).unwrap() = 0.5;
+        *points.get_mut([1, 3]).unwrap() = 0.0;
+        *points.get_mut([2, 3]).unwrap() = 0.1;
+        *points.get_mut([0, 4]).unwrap() = 1.0;
+        *points.get_mut([1, 4]).unwrap() = 0.5;
+        *points.get_mut([2, 4]).unwrap() = 0.5;
+        *points.get_mut([0, 5]).unwrap() = 0.5;
+        *points.get_mut([1, 5]).unwrap() = 0.5;
+        *points.get_mut([2, 5]).unwrap() = 1.0;
+        e.tabulate(&points, 0, &mut data);
+
+        for pt in 0..6 {
+            let x = *points.get([0, pt]).unwrap();
+            let y = *points.get([1, pt]).unwrap();
+            let z = *points.get([2, pt]).unwrap();
+            for (i, basis_f) in [
+                [(1.0 - y) * (1.0 - z), 0.0, 0.0],
+                [0.0, (1.0 - x) * (1.0 - z), 0.0],
+                [0.0, 0.0, (1.0 - x) * (1.0 - y)],
+                [0.0, x * (1.0 - z), 0.0],
+                [0.0, 0.0, x * (1.0 - y)],
+                [y * (1.0 - z), 0.0, 0.0],
+                [0.0, 0.0, (1.0 - x) * y],
+                [0.0, 0.0, x * y],
+                [(1.0 - y) * z, 0.0, 0.0],
+                [0.0, (1.0 - x) * z, 0.0],
+                [0.0, x * z, 0.0],
+                [y * z, 0.0, 0.0],
+            ]
+            .iter()
+            .enumerate()
+            {
+                for (d, value) in basis_f.iter().enumerate() {
+                    assert_relative_eq!(*data.get([0, pt, i, d]).unwrap(), value, epsilon = 1e-14);
+                }
+            }
+        }
+        check_dofs(e);
+    }
+
+    #[test]
+    fn test_nedelec_2_hexahedron() {
+        let e = nedelec::create::<f64>(ReferenceCellType::Hexahedron, 2, Continuity::Standard);
+        assert_eq!(e.value_size(), 3);
+        check_dofs(e);
+    }
+
     macro_rules! test_entity_closure_dofs_lagrange {
         ($cell:ident, $degree:expr) => {
             paste! {
diff --git a/src/ciarlet/nedelec.rs b/src/ciarlet/nedelec.rs
index bb9a996..ed32b72 100644
--- a/src/ciarlet/nedelec.rs
+++ b/src/ciarlet/nedelec.rs
@@ -265,13 +265,13 @@ fn create_simplex<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> +
             let mut pts = rlst_dynamic_array2!(T::Real, [tdim, face_q.npoints]);
             let mut mat = rlst_dynamic_array3!(T, [2 * pdim_face_minus2, tdim, face_q.npoints]);
 
-            for tangent in 0..2 {
-                for (w_i, wt) in face_q.weights.iter().enumerate() {
-                    for i in 0..tdim {
-                        pts[[i, w_i]] = vertices[face[0]][i]
-                            + (vertices[face[1]][i] - vertices[face[0]][i]) * face_pts[[0, w_i]]
-                            + (vertices[face[2]][i] - vertices[face[0]][i]) * face_pts[[1, w_i]];
+            for (w_i, wt) in face_q.weights.iter().enumerate() {
+                for i in 0..tdim {
+                    pts[[i, w_i]] = vertices[face[0]][i]
+                        + (vertices[face[1]][i] - vertices[face[0]][i]) * face_pts[[0, w_i]]
+                        + (vertices[face[2]][i] - vertices[face[0]][i]) * face_pts[[1, w_i]];
 
+                    for tangent in 0..2 {
                         for j in 0..pdim_face_minus2 {
                             mat[[tangent * pdim_face_minus2 + j, i, w_i]] = T::from(*wt).unwrap()
                                 * face_phi[[0, j, w_i]]
@@ -331,6 +331,311 @@ fn create_simplex<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> +
     )
 }
 
+fn create_tp<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> + MatrixInverse>(
+    cell_type: ReferenceCellType,
+    degree: usize,
+    continuity: Continuity,
+) -> CiarletElement<T> {
+    if cell_type != ReferenceCellType::Quadrilateral && cell_type != ReferenceCellType::Hexahedron {
+        panic!("Invalid cell: {cell_type:?}");
+    }
+
+    if degree < 1 {
+        panic!("Degree must be at least 1");
+    }
+
+    let tdim = reference_cell::dim(cell_type);
+    let cell_q = gauss_jacobi_rule(cell_type, 2 * degree).unwrap();
+    let pts_t = cell_q
+        .points
+        .iter()
+        .map(|i| TReal::from(*i).unwrap())
+        .collect::<Vec<_>>();
+    let pts = rlst_array_from_slice2!(&pts_t, [tdim, cell_q.npoints]);
+
+    let mut phi = rlst_dynamic_array3![T, legendre_shape(cell_type, &pts, degree, 0)];
+    tabulate_legendre_polynomials(cell_type, &pts, degree, 0, &mut phi);
+
+    let pdim = phi.shape()[1];
+
+    let pdim_edge = polynomial_count(ReferenceCellType::Interval, degree);
+    let pdim_edge_minus1 = polynomial_count(ReferenceCellType::Interval, degree - 1);
+    let pdim_edge_minus2 = if degree < 2 {
+        0
+    } else {
+        polynomial_count(ReferenceCellType::Interval, degree - 2)
+    };
+
+    let entity_counts = reference_cell::entity_counts(cell_type);
+
+    let mut wcoeffs = rlst_dynamic_array3!(
+        T,
+        [
+            entity_counts[1] * pdim_edge_minus1
+                + entity_counts[2] * 2 * pdim_edge_minus1 * pdim_edge_minus2
+                + entity_counts[3] * 3 * pdim_edge_minus1 * pdim_edge_minus2 * pdim_edge_minus2,
+            tdim,
+            pdim
+        ]
+    );
+
+    // vector polynomials of degree <= n-1
+    if tdim == 2 {
+        for i in 0..pdim_edge_minus1 {
+            for j in 0..pdim_edge {
+                *wcoeffs
+                    .get_mut([i * pdim_edge + j, 0, i * pdim_edge + j])
+                    .unwrap() = T::from(1.0).unwrap();
+                *wcoeffs
+                    .get_mut([
+                        pdim_edge_minus1 * pdim_edge + j * pdim_edge_minus1 + i,
+                        1,
+                        j * pdim_edge + i,
+                    ])
+                    .unwrap() = T::from(1.0).unwrap();
+            }
+        }
+    } else {
+        for i in 0..pdim_edge_minus1 {
+            for j in 0..pdim_edge {
+                for k in 0..pdim_edge {
+                    *wcoeffs
+                        .get_mut([
+                            i * pdim_edge.pow(2) + j * pdim_edge + k,
+                            0,
+                            i * pdim_edge.pow(2) + j * pdim_edge + k,
+                        ])
+                        .unwrap() = T::from(1.0).unwrap();
+                    *wcoeffs
+                        .get_mut([
+                            pdim_edge.pow(2) * pdim_edge_minus1
+                                + k * pdim_edge * pdim_edge_minus1
+                                + i * pdim_edge
+                                + j,
+                            1,
+                            k * pdim_edge.pow(2) + i * pdim_edge + j,
+                        ])
+                        .unwrap() = T::from(1.0).unwrap();
+                    *wcoeffs
+                        .get_mut([
+                            pdim_edge.pow(2) * pdim_edge_minus1 * 2
+                                + j * pdim_edge * pdim_edge_minus1
+                                + k * pdim_edge_minus1
+                                + i,
+                            2,
+                            j * pdim_edge.pow(2) + k * pdim_edge + i,
+                        ])
+                        .unwrap() = T::from(1.0).unwrap();
+                }
+            }
+        }
+    }
+
+    let mut x = [vec![], vec![], vec![], vec![]];
+    let mut m = [vec![], vec![], vec![], vec![]];
+
+    let vertices = reference_cell::vertices::<T::Real>(cell_type);
+
+    for _ in 0..entity_counts[0] {
+        x[0].push(rlst_dynamic_array2!(T::Real, [tdim, 0]));
+        m[0].push(rlst_dynamic_array3!(T, [0, tdim, 0]));
+    }
+
+    // DOFs on edges
+    let edge_q = gauss_jacobi_rule(ReferenceCellType::Interval, 2 * degree - 1).unwrap();
+    let edge_pts_t = edge_q
+        .points
+        .iter()
+        .map(|i| TReal::from(*i).unwrap())
+        .collect::<Vec<_>>();
+    let edge_pts = rlst_array_from_slice2!(&edge_pts_t, [1, edge_q.npoints]);
+
+    let mut edge_phi = rlst_dynamic_array3![
+        T,
+        legendre_shape(ReferenceCellType::Interval, &edge_pts, degree - 1, 0)
+    ];
+    tabulate_legendre_polynomials(
+        ReferenceCellType::Interval,
+        &edge_pts,
+        degree - 1,
+        0,
+        &mut edge_phi,
+    );
+
+    for edge in reference_cell::edges(cell_type) {
+        let mut pts = rlst_dynamic_array2!(T::Real, [tdim, edge_q.npoints]);
+        let mut mat = rlst_dynamic_array3!(T, [pdim_edge_minus1, tdim, edge_q.npoints]);
+
+        for (w_i, (pt, wt)) in izip!(&edge_pts_t, &edge_q.weights).enumerate() {
+            for i in 0..tdim {
+                pts[[i, w_i]] =
+                    vertices[edge[0]][i] + (vertices[edge[1]][i] - vertices[edge[0]][i]) * *pt;
+
+                for j in 0..pdim_edge_minus1 {
+                    mat[[j, i, w_i]] = T::from(*wt).unwrap()
+                        * edge_phi[[0, j, w_i]]
+                        * T::from(vertices[edge[1]][i] - vertices[edge[0]][i]).unwrap();
+                }
+            }
+        }
+
+        x[1].push(pts);
+        m[1].push(mat);
+    }
+
+    // DOFs on faces
+    if degree == 1 {
+        for _ in 0..entity_counts[2] {
+            x[2].push(rlst_dynamic_array2!(T::Real, [tdim, 0]));
+            m[2].push(rlst_dynamic_array3!(T, [0, tdim, 0]))
+        }
+    } else {
+        let face_q = gauss_jacobi_rule(ReferenceCellType::Quadrilateral, 2 * degree - 1).unwrap();
+        let face_pts_t = face_q
+            .points
+            .iter()
+            .map(|i| TReal::from(*i).unwrap())
+            .collect::<Vec<_>>();
+        let face_pts = rlst_array_from_slice2!(&face_pts_t, [2, face_q.npoints]);
+
+        let mut face_phi = rlst_dynamic_array3![
+            T,
+            legendre_shape(ReferenceCellType::Quadrilateral, &face_pts, degree - 1, 0)
+        ];
+        tabulate_legendre_polynomials(
+            ReferenceCellType::Quadrilateral,
+            &face_pts,
+            degree - 1,
+            0,
+            &mut face_phi,
+        );
+
+        for face in reference_cell::faces(cell_type) {
+            let mut pts = rlst_dynamic_array2!(T::Real, [tdim, face_q.npoints]);
+            let mut mat = rlst_dynamic_array3!(
+                T,
+                [
+                    2 * pdim_edge_minus2 * pdim_edge_minus1,
+                    tdim,
+                    face_q.npoints
+                ]
+            );
+
+            for (w_i, wt) in face_q.weights.iter().enumerate() {
+                for i in 0..tdim {
+                    pts[[i, w_i]] = vertices[face[0]][i]
+                        + (vertices[face[1]][i] - vertices[face[0]][i]) * face_pts[[0, w_i]]
+                        + (vertices[face[2]][i] - vertices[face[0]][i]) * face_pts[[1, w_i]];
+                }
+                for i in 0..pdim_edge_minus2 {
+                    for j in 0..pdim_edge_minus1 {
+                        let index = 2 * (i * pdim_edge_minus1 + j);
+                        let entry0 =
+                            T::from(*wt).unwrap() * face_phi[[0, j * pdim_edge_minus1 + i, w_i]];
+                        let entry1 =
+                            T::from(*wt).unwrap() * face_phi[[0, i * pdim_edge_minus1 + j, w_i]];
+                        for d in 0..tdim {
+                            mat[[index, d, w_i]] = entry0
+                                * T::from(vertices[face[1]][d] - vertices[face[0]][d]).unwrap();
+                            mat[[index + 1, d, w_i]] = entry1
+                                * T::from(vertices[face[2]][d] - vertices[face[0]][d]).unwrap();
+                        }
+                    }
+                }
+            }
+            x[2].push(pts);
+            m[2].push(mat);
+        }
+    }
+    // DOFs on volume
+    if tdim == 3 {
+        if degree == 1 {
+            x[3].push(rlst_dynamic_array2!(T::Real, [tdim, 0]));
+            m[3].push(rlst_dynamic_array3!(T, [0, tdim, 0]))
+        } else {
+            let interior_q =
+                gauss_jacobi_rule(ReferenceCellType::Hexahedron, 2 * degree - 1).unwrap();
+            let interior_pts_t = interior_q
+                .points
+                .iter()
+                .map(|i| TReal::from(*i).unwrap())
+                .collect::<Vec<_>>();
+            let interior_pts = rlst_array_from_slice2!(&interior_pts_t, [3, interior_q.npoints]);
+
+            let mut interior_phi = rlst_dynamic_array3![
+                T,
+                legendre_shape(ReferenceCellType::Hexahedron, &interior_pts, degree - 1, 0)
+            ];
+            tabulate_legendre_polynomials(
+                ReferenceCellType::Hexahedron,
+                &interior_pts,
+                degree - 1,
+                0,
+                &mut interior_phi,
+            );
+
+            let mut pts = rlst_dynamic_array2!(T::Real, [tdim, interior_q.npoints]);
+            let mut mat = rlst_dynamic_array3!(
+                T,
+                [
+                    3 * pdim_edge_minus2.pow(2) * pdim_edge_minus1,
+                    tdim,
+                    interior_q.npoints
+                ]
+            );
+
+            for (w_i, wt) in interior_q.weights.iter().enumerate() {
+                for i in 0..tdim {
+                    pts[[i, w_i]] = interior_pts[[i, w_i]];
+                }
+                for i in 0..pdim_edge_minus2 {
+                    for j in 0..pdim_edge_minus2 {
+                        for k in 0..pdim_edge_minus1 {
+                            let index = 3
+                                * (i * pdim_edge_minus1 * pdim_edge_minus2
+                                    + j * pdim_edge_minus1
+                                    + k);
+                            mat[[index, 0, w_i]] = T::from(*wt).unwrap()
+                                * interior_phi[[
+                                    0,
+                                    k * pdim_edge_minus1.pow(2) + j * pdim_edge_minus1 + i,
+                                    w_i,
+                                ]];
+                            mat[[index + 1, 1, w_i]] = T::from(*wt).unwrap()
+                                * interior_phi[[
+                                    0,
+                                    i * pdim_edge_minus1.pow(2) + k * pdim_edge_minus1 + j,
+                                    w_i,
+                                ]];
+                            mat[[index + 2, 2, w_i]] = T::from(*wt).unwrap()
+                                * interior_phi[[
+                                    0,
+                                    j * pdim_edge_minus1.pow(2) + i * pdim_edge_minus1 + k,
+                                    w_i,
+                                ]];
+                        }
+                    }
+                }
+            }
+            x[3].push(pts);
+            m[3].push(mat);
+        }
+    }
+
+    CiarletElement::create(
+        "Nedelec (first kind)".to_string(),
+        cell_type,
+        degree,
+        vec![tdim],
+        wcoeffs,
+        x,
+        m,
+        MapType::CovariantPiola,
+        continuity,
+        degree,
+    )
+}
+
 /// Create a Nedelec (first kind) element
 pub fn create<T: RlstScalar + MatrixInverse>(
     cell_type: ReferenceCellType,
@@ -339,6 +644,10 @@ pub fn create<T: RlstScalar + MatrixInverse>(
 ) -> CiarletElement<T> {
     if cell_type == ReferenceCellType::Triangle || cell_type == ReferenceCellType::Tetrahedron {
         create_simplex(cell_type, degree, continuity)
+    } else if cell_type == ReferenceCellType::Quadrilateral
+        || cell_type == ReferenceCellType::Hexahedron
+    {
+        create_tp(cell_type, degree, continuity)
     } else {
         panic!("Invalid cell: {cell_type:?}");
     }
diff --git a/src/ciarlet/raviart_thomas.rs b/src/ciarlet/raviart_thomas.rs
index e1662fb..e2b87c2 100644
--- a/src/ciarlet/raviart_thomas.rs
+++ b/src/ciarlet/raviart_thomas.rs
@@ -191,6 +191,294 @@ fn create_simplex<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> +
     )
 }
 
+fn create_tp<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> + MatrixInverse>(
+    cell_type: ReferenceCellType,
+    degree: usize,
+    continuity: Continuity,
+) -> CiarletElement<T> {
+    if cell_type != ReferenceCellType::Quadrilateral && cell_type != ReferenceCellType::Hexahedron {
+        panic!("Invalid cell: {cell_type:?}");
+    }
+
+    if degree < 1 {
+        panic!("Degree must be at least 1");
+    }
+
+    let tdim = reference_cell::dim(cell_type);
+    let cell_q = gauss_jacobi_rule(cell_type, 2 * degree).unwrap();
+    let pts_t = cell_q
+        .points
+        .iter()
+        .map(|i| TReal::from(*i).unwrap())
+        .collect::<Vec<_>>();
+    let pts = rlst_array_from_slice2!(&pts_t, [tdim, cell_q.npoints]);
+
+    let mut phi = rlst_dynamic_array3![T, legendre_shape(cell_type, &pts, degree, 0)];
+    tabulate_legendre_polynomials(cell_type, &pts, degree, 0, &mut phi);
+
+    let pdim = phi.shape()[1];
+
+    let pdim_edge = polynomial_count(ReferenceCellType::Interval, degree);
+    let pdim_edge_minus1 = polynomial_count(ReferenceCellType::Interval, degree - 1);
+    let pdim_edge_minus2 = if degree < 2 {
+        0
+    } else {
+        polynomial_count(ReferenceCellType::Interval, degree - 2)
+    };
+
+    let facet_type = if tdim == 2 {
+        ReferenceCellType::Interval
+    } else {
+        ReferenceCellType::Quadrilateral
+    };
+    let pdim_facet_minus1 = polynomial_count(facet_type, degree - 1);
+    let pdim_internal = tdim * pdim_edge_minus2 * pdim_edge_minus1.pow(tdim as u32 - 1);
+
+    let entity_counts = reference_cell::entity_counts(cell_type);
+
+    let mut wcoeffs = rlst_dynamic_array3!(
+        T,
+        [
+            entity_counts[tdim - 1] * pdim_facet_minus1 + entity_counts[tdim] * pdim_internal,
+            tdim,
+            pdim,
+        ]
+    );
+
+    // vector polynomials of degree <= n-1
+    if tdim == 2 {
+        for i in 0..pdim_edge_minus1 {
+            for j in 0..pdim_edge {
+                *wcoeffs
+                    .get_mut([j * pdim_edge_minus1 + i, 0, j * pdim_edge + i])
+                    .unwrap() = T::from(1.0).unwrap();
+                *wcoeffs
+                    .get_mut([
+                        pdim_edge_minus1 * pdim_edge + i * pdim_edge + j,
+                        1,
+                        i * pdim_edge + j,
+                    ])
+                    .unwrap() = T::from(1.0).unwrap();
+            }
+        }
+    } else {
+        for i in 0..pdim_edge_minus1 {
+            for j in 0..pdim_edge_minus1 {
+                for k in 0..pdim_edge {
+                    *wcoeffs
+                        .get_mut([
+                            k * pdim_edge_minus1.pow(2) + j * pdim_edge_minus1 + i,
+                            0,
+                            k * pdim_edge.pow(2) + j * pdim_edge + i,
+                        ])
+                        .unwrap() = T::from(1.0).unwrap();
+                    *wcoeffs
+                        .get_mut([
+                            pdim_edge * pdim_edge_minus1.pow(2)
+                                + i * pdim_edge * pdim_edge_minus1
+                                + k * pdim_edge_minus1
+                                + j,
+                            1,
+                            i * pdim_edge.pow(2) + k * pdim_edge + j,
+                        ])
+                        .unwrap() = T::from(1.0).unwrap();
+                    *wcoeffs
+                        .get_mut([
+                            pdim_edge * pdim_edge_minus1.pow(2) * 2
+                                + j * pdim_edge * pdim_edge_minus1
+                                + i * pdim_edge
+                                + k,
+                            2,
+                            j * pdim_edge.pow(2) + i * pdim_edge + k,
+                        ])
+                        .unwrap() = T::from(1.0).unwrap();
+                }
+            }
+        }
+    }
+
+    let mut x = [vec![], vec![], vec![], vec![]];
+    let mut m = [vec![], vec![], vec![], vec![]];
+
+    let vertices = reference_cell::vertices::<T::Real>(cell_type);
+
+    for d in 0..tdim - 1 {
+        for _ in 0..entity_counts[d] {
+            x[d].push(rlst_dynamic_array2!(T::Real, [tdim, 0]));
+            m[d].push(rlst_dynamic_array3!(T, [0, tdim, 0]));
+        }
+    }
+    // DOFs on facets
+    let facet_q = gauss_jacobi_rule(facet_type, 2 * degree - 1).unwrap();
+    let facet_pts_t = facet_q
+        .points
+        .iter()
+        .map(|i| TReal::from(*i).unwrap())
+        .collect::<Vec<_>>();
+    let facet_pts = rlst_array_from_slice2!(&facet_pts_t, [tdim - 1, facet_q.npoints]);
+
+    let mut facet_phi =
+        rlst_dynamic_array3![T, legendre_shape(facet_type, &facet_pts, degree - 1, 0)];
+    tabulate_legendre_polynomials(facet_type, &facet_pts, degree - 1, 0, &mut facet_phi);
+
+    for (facet, normal) in izip!(
+        reference_cell::facets(cell_type),
+        reference_cell::facet_normals::<T::Real>(cell_type),
+    ) {
+        let mut pts = rlst_dynamic_array2!(T::Real, [tdim, facet_q.npoints]);
+        let mut mat = rlst_dynamic_array3!(T, [pdim_facet_minus1, tdim, facet_q.npoints]);
+
+        for (w_i, wt) in facet_q.weights.iter().enumerate() {
+            for d in 0..tdim {
+                pts[[d, w_i]] = vertices[facet[0]][d]
+                    + (0..tdim - 1)
+                        .map(|x| {
+                            (vertices[facet[usize::pow(2, x as u32)]][d] - vertices[facet[0]][d])
+                                * facet_pts[[x, w_i]]
+                        })
+                        .sum::<T::Real>();
+                for j in 0..facet_phi.shape()[1] {
+                    mat[[j, d, w_i]] = T::from(*wt).unwrap()
+                        * facet_phi[[0, j, w_i]]
+                        * T::from(normal[d]).unwrap();
+                }
+            }
+        }
+
+        x[tdim - 1].push(pts);
+        m[tdim - 1].push(mat);
+    }
+
+    // DOFs on interior
+    if degree == 1 {
+        x[tdim].push(rlst_dynamic_array2!(T::Real, [tdim, 0]));
+        m[tdim].push(rlst_dynamic_array3!(T, [0, tdim, 0]))
+    } else if tdim == 2 {
+        let face_q = gauss_jacobi_rule(ReferenceCellType::Quadrilateral, 2 * degree - 1).unwrap();
+        let face_pts_t = face_q
+            .points
+            .iter()
+            .map(|i| TReal::from(*i).unwrap())
+            .collect::<Vec<_>>();
+        let face_pts = rlst_array_from_slice2!(&face_pts_t, [2, face_q.npoints]);
+
+        let mut face_phi = rlst_dynamic_array3![
+            T,
+            legendre_shape(ReferenceCellType::Quadrilateral, &face_pts, degree - 1, 0)
+        ];
+        tabulate_legendre_polynomials(
+            ReferenceCellType::Quadrilateral,
+            &face_pts,
+            degree - 1,
+            0,
+            &mut face_phi,
+        );
+
+        let mut pts = rlst_dynamic_array2!(T::Real, [tdim, face_q.npoints]);
+        let mut mat = rlst_dynamic_array3!(
+            T,
+            [
+                2 * pdim_edge_minus2 * pdim_edge_minus1,
+                tdim,
+                face_q.npoints
+            ]
+        );
+
+        for (w_i, wt) in face_q.weights.iter().enumerate() {
+            for i in 0..tdim {
+                pts[[i, w_i]] = face_pts[[i, w_i]];
+            }
+            for i in 0..pdim_edge_minus2 {
+                for j in 0..pdim_edge_minus1 {
+                    let index = 2 * (i * pdim_edge_minus1 + j);
+                    mat[[index, 0, w_i]] =
+                        T::from(*wt).unwrap() * face_phi[[0, i * pdim_edge_minus1 + j, w_i]];
+                    mat[[index + 1, 1, w_i]] =
+                        T::from(*wt).unwrap() * face_phi[[0, j * pdim_edge_minus1 + i, w_i]];
+                }
+            }
+        }
+        x[2].push(pts);
+        m[2].push(mat);
+    } else {
+        let interior_q = gauss_jacobi_rule(ReferenceCellType::Hexahedron, 2 * degree - 1).unwrap();
+        let interior_pts_t = interior_q
+            .points
+            .iter()
+            .map(|i| TReal::from(*i).unwrap())
+            .collect::<Vec<_>>();
+        let interior_pts = rlst_array_from_slice2!(&interior_pts_t, [3, interior_q.npoints]);
+
+        let mut interior_phi = rlst_dynamic_array3![
+            T,
+            legendre_shape(ReferenceCellType::Hexahedron, &interior_pts, degree - 1, 0)
+        ];
+        tabulate_legendre_polynomials(
+            ReferenceCellType::Hexahedron,
+            &interior_pts,
+            degree - 1,
+            0,
+            &mut interior_phi,
+        );
+
+        let mut pts = rlst_dynamic_array2!(T::Real, [tdim, interior_q.npoints]);
+        let mut mat = rlst_dynamic_array3!(
+            T,
+            [
+                3 * pdim_edge_minus1.pow(2) * pdim_edge_minus2,
+                tdim,
+                interior_q.npoints
+            ]
+        );
+
+        for (w_i, wt) in interior_q.weights.iter().enumerate() {
+            for i in 0..tdim {
+                pts[[i, w_i]] = interior_pts[[i, w_i]];
+            }
+            for i in 0..pdim_edge_minus2 {
+                for j in 0..pdim_edge_minus1 {
+                    for k in 0..pdim_edge_minus1 {
+                        let index = 3 * (i * pdim_edge_minus1.pow(2) + j * pdim_edge_minus1 + k);
+                        mat[[index, 0, w_i]] = T::from(*wt).unwrap()
+                            * interior_phi[[
+                                0,
+                                i * pdim_edge_minus1.pow(2) + j * pdim_edge_minus1 + k,
+                                w_i,
+                            ]];
+                        mat[[index + 1, 1, w_i]] = T::from(*wt).unwrap()
+                            * interior_phi[[
+                                0,
+                                k * pdim_edge_minus1.pow(2) + i * pdim_edge_minus1 + j,
+                                w_i,
+                            ]];
+                        mat[[index + 2, 2, w_i]] = T::from(*wt).unwrap()
+                            * interior_phi[[
+                                0,
+                                j * pdim_edge_minus1.pow(2) + k * pdim_edge_minus1 + i,
+                                w_i,
+                            ]];
+                    }
+                }
+            }
+        }
+        x[3].push(pts);
+        m[3].push(mat);
+    }
+
+    CiarletElement::create(
+        "Raviart-Thomas".to_string(),
+        cell_type,
+        degree,
+        vec![tdim],
+        wcoeffs,
+        x,
+        m,
+        MapType::ContravariantPiola,
+        continuity,
+        degree,
+    )
+}
+
 /// Create a Raviart-Thomas element
 pub fn create<T: RlstScalar + MatrixInverse>(
     cell_type: ReferenceCellType,
@@ -199,6 +487,10 @@ pub fn create<T: RlstScalar + MatrixInverse>(
 ) -> CiarletElement<T> {
     if cell_type == ReferenceCellType::Triangle || cell_type == ReferenceCellType::Tetrahedron {
         create_simplex(cell_type, degree, continuity)
+    } else if cell_type == ReferenceCellType::Quadrilateral
+        || cell_type == ReferenceCellType::Hexahedron
+    {
+        create_tp(cell_type, degree, continuity)
     } else {
         panic!("Invalid cell: {cell_type:?}");
     }
diff --git a/src/polynomials.rs b/src/polynomials.rs
index fadb1ef..df379a3 100644
--- a/src/polynomials.rs
+++ b/src/polynomials.rs
@@ -341,7 +341,7 @@ mod test {
 
         let mut data = rlst_dynamic_array3!(
             f64,
-            legendre_shape(ReferenceCellType::Quadrilateral, &points, degree, 1,)
+            legendre_shape(ReferenceCellType::Quadrilateral, &points, degree, 1)
         );
         tabulate_legendre_polynomials(
             ReferenceCellType::Quadrilateral,
@@ -360,42 +360,42 @@ mod test {
             assert_relative_eq!(*data.get([1, 0, k]).unwrap(), 0.0, epsilon = 1e-12);
             assert_relative_eq!(*data.get([2, 0, k]).unwrap(), 0.0, epsilon = 1e-12);
 
-            // 1 => sqrt(3)*(2x - 1)
+            // 3 => sqrt(3)*(2x - 1)
             assert_relative_eq!(
-                *data.get([0, 1, k]).unwrap(),
+                *data.get([0, 3, k]).unwrap(),
                 f64::sqrt(3.0) * (2.0 * x - 1.0),
                 epsilon = 1e-12
             );
             assert_relative_eq!(
-                *data.get([1, 1, k]).unwrap(),
+                *data.get([1, 3, k]).unwrap(),
                 2.0 * f64::sqrt(3.0),
                 epsilon = 1e-12
             );
-            assert_relative_eq!(*data.get([2, 1, k]).unwrap(), 0.0, epsilon = 1e-12);
+            assert_relative_eq!(*data.get([2, 3, k]).unwrap(), 0.0, epsilon = 1e-12);
 
-            // 2 => sqrt(5)*(6x^2 - 6x + 1)
+            // 6 => sqrt(5)*(6x^2 - 6x + 1)
             assert_relative_eq!(
-                *data.get([0, 2, k]).unwrap(),
+                *data.get([0, 6, k]).unwrap(),
                 f64::sqrt(5.0) * (6.0 * x * x - 6.0 * x + 1.0),
                 epsilon = 1e-12
             );
             assert_relative_eq!(
-                *data.get([1, 2, k]).unwrap(),
+                *data.get([1, 6, k]).unwrap(),
                 f64::sqrt(5.0) * (12.0 * x - 6.0),
                 epsilon = 1e-12
             );
-            assert_relative_eq!(*data.get([2, 2, k]).unwrap(), 0.0, epsilon = 1e-12);
+            assert_relative_eq!(*data.get([2, 6, k]).unwrap(), 0.0, epsilon = 1e-12);
 
-            // 3 => sqrt(3)*(2y - 1)
+            // 1 => sqrt(3)*(2y - 1)
             assert_relative_eq!(
-                *data.get([0, 3, k]).unwrap(),
+                *data.get([0, 1, k]).unwrap(),
                 f64::sqrt(3.0) * (2.0 * y - 1.0),
                 epsilon = 1e-12
             );
 
-            assert_relative_eq!(*data.get([1, 3, k]).unwrap(), 0.0, epsilon = 1e-12);
+            assert_relative_eq!(*data.get([1, 1, k]).unwrap(), 0.0, epsilon = 1e-12);
             assert_relative_eq!(
-                *data.get([2, 3, k]).unwrap(),
+                *data.get([2, 1, k]).unwrap(),
                 2.0 * f64::sqrt(3.0),
                 epsilon = 1e-12
             );
@@ -417,49 +417,49 @@ mod test {
                 epsilon = 1e-12
             );
 
-            // 5 => sqrt(15)*(6x^2 - 6x + 1)*(2y - 1)
+            // 7 => sqrt(15)*(6x^2 - 6x + 1)*(2y - 1)
             assert_relative_eq!(
-                *data.get([0, 5, k]).unwrap(),
+                *data.get([0, 7, k]).unwrap(),
                 f64::sqrt(15.0) * (6.0 * x * x - 6.0 * x + 1.0) * (2.0 * y - 1.0),
                 epsilon = 1e-12
             );
             assert_relative_eq!(
-                *data.get([1, 5, k]).unwrap(),
+                *data.get([1, 7, k]).unwrap(),
                 f64::sqrt(15.0) * (12.0 * x - 6.0) * (2.0 * y - 1.0),
                 epsilon = 1e-12
             );
             assert_relative_eq!(
-                *data.get([2, 5, k]).unwrap(),
+                *data.get([2, 7, k]).unwrap(),
                 2.0 * f64::sqrt(15.0) * (6.0 * x * x - 6.0 * x + 1.0),
                 epsilon = 1e-12
             );
 
-            // 6 => sqrt(5)*(6y^2 - 6y + 1)
+            // 2 => sqrt(5)*(6y^2 - 6y + 1)
             assert_relative_eq!(
-                *data.get([0, 6, k]).unwrap(),
+                *data.get([0, 2, k]).unwrap(),
                 f64::sqrt(5.0) * (6.0 * y * y - 6.0 * y + 1.0),
                 epsilon = 1e-12
             );
-            assert_relative_eq!(*data.get([1, 6, k]).unwrap(), 0.0, epsilon = 1e-12);
+            assert_relative_eq!(*data.get([1, 2, k]).unwrap(), 0.0, epsilon = 1e-12);
             assert_relative_eq!(
-                *data.get([2, 6, k]).unwrap(),
+                *data.get([2, 2, k]).unwrap(),
                 f64::sqrt(5.0) * (12.0 * y - 6.0),
                 epsilon = 1e-12
             );
 
-            // 7 => sqrt(15)*(2x - 1)*(6y^2 - 6y + 1)
+            // 5 => sqrt(15)*(2x - 1)*(6y^2 - 6y + 1)
             assert_relative_eq!(
-                *data.get([0, 7, k]).unwrap(),
+                *data.get([0, 5, k]).unwrap(),
                 f64::sqrt(15.0) * (2.0 * x - 1.0) * (6.0 * y * y - 6.0 * y + 1.0),
                 epsilon = 1e-12
             );
             assert_relative_eq!(
-                *data.get([1, 7, k]).unwrap(),
+                *data.get([1, 5, k]).unwrap(),
                 2.0 * f64::sqrt(15.0) * (6.0 * y * y - 6.0 * y + 1.0),
                 epsilon = 1e-12
             );
             assert_relative_eq!(
-                *data.get([2, 7, k]).unwrap(),
+                *data.get([2, 5, k]).unwrap(),
                 f64::sqrt(15.0) * (2.0 * x - 1.0) * (12.0 * y - 6.0),
                 epsilon = 1e-12
             );
diff --git a/src/polynomials/legendre.rs b/src/polynomials/legendre.rs
index e388bd1..e5307e2 100644
--- a/src/polynomials/legendre.rs
+++ b/src/polynomials/legendre.rs
@@ -12,7 +12,7 @@ fn tri_index(i: usize, j: usize) -> usize {
 }
 
 fn quad_index(i: usize, j: usize, n: usize) -> usize {
-    j * (n + 1) + i
+    i * (n + 1) + j
 }
 
 fn tet_index(i: usize, j: usize, k: usize) -> usize {
@@ -20,7 +20,7 @@ fn tet_index(i: usize, j: usize, k: usize) -> usize {
 }
 
 fn hex_index(i: usize, j: usize, k: usize, n: usize) -> usize {
-    k * (n + 1) * (n + 1) + j * (n + 1) + i
+    i * (n + 1) * (n + 1) + j * (n + 1) + k
 }
 
 /// The coefficients in the Jacobi Polynomial recurrence relation
diff --git a/src/reference_cell.rs b/src/reference_cell.rs
index b590629..4d77d20 100644
--- a/src/reference_cell.rs
+++ b/src/reference_cell.rs
@@ -98,7 +98,7 @@ pub fn midpoint<T: RlstScalar<Real = T>>(cell: ReferenceCellType) -> Vec<T> {
     }
 }
 
-/// The edges of the reference cell
+/// The edges (dimension 1 entities) of the reference cell
 pub fn edges(cell: ReferenceCellType) -> Vec<Vec<usize>> {
     match cell {
         ReferenceCellType::Point => vec![],
@@ -151,7 +151,7 @@ pub fn edges(cell: ReferenceCellType) -> Vec<Vec<usize>> {
     }
 }
 
-/// The faces of the reference cell
+/// The faces (dimension 2 entities) of the reference cell
 pub fn faces(cell: ReferenceCellType) -> Vec<Vec<usize>> {
     match cell {
         ReferenceCellType::Point => vec![],
@@ -186,6 +186,68 @@ pub fn faces(cell: ReferenceCellType) -> Vec<Vec<usize>> {
     }
 }
 
+/// The volumes (dimension 3 entities) of the reference cell
+pub fn volumes(cell: ReferenceCellType) -> Vec<Vec<usize>> {
+    match cell {
+        ReferenceCellType::Point => vec![],
+        ReferenceCellType::Interval => vec![],
+        ReferenceCellType::Triangle => vec![],
+        ReferenceCellType::Quadrilateral => vec![],
+        ReferenceCellType::Tetrahedron => vec![vec![0, 1, 2, 3]],
+        ReferenceCellType::Hexahedron => vec![vec![0, 1, 2, 3, 4, 5, 6, 7]],
+        ReferenceCellType::Prism => vec![vec![0, 1, 2, 3, 4, 5]],
+        ReferenceCellType::Pyramid => vec![vec![0, 1, 2, 3, 4]],
+    }
+}
+
+/// The facets (dimension dim-1 entities) of the reference cell
+pub fn facets(cell: ReferenceCellType) -> Vec<Vec<usize>> {
+    match dim(cell) {
+        0 => vec![],
+        1 => volumes(cell)[0]
+            .iter()
+            .map(|x| vec![*x])
+            .collect::<Vec<_>>(),
+        2 => edges(cell),
+        3 => faces(cell),
+        _ => {
+            panic!("Invalid dimension");
+        }
+    }
+}
+
+/// The ridges (dimension dim-2 entities) of the reference cell
+pub fn ridges(cell: ReferenceCellType) -> Vec<Vec<usize>> {
+    match dim(cell) {
+        0 => vec![],
+        1 => vec![],
+        2 => volumes(cell)[0]
+            .iter()
+            .map(|x| vec![*x])
+            .collect::<Vec<_>>(),
+        3 => edges(cell),
+        _ => {
+            panic!("Invalid dimension");
+        }
+    }
+}
+
+/// The peaks (dimension dim-3 entities) of the reference cell
+pub fn peaks(cell: ReferenceCellType) -> Vec<Vec<usize>> {
+    match dim(cell) {
+        0 => vec![],
+        1 => vec![],
+        2 => vec![],
+        3 => volumes(cell)[0]
+            .iter()
+            .map(|x| vec![*x])
+            .collect::<Vec<_>>(),
+        _ => {
+            panic!("Invalid dimension");
+        }
+    }
+}
+
 /// The normals to the facets of the reference cell. The length of each normal is the volume of the facet
 pub fn facet_normals<T: RlstScalar<Real = T>>(cell: ReferenceCellType) -> Vec<Vec<T>> {
     let zero = T::from(0.0).unwrap();
@@ -244,20 +306,6 @@ pub fn facet_unit_normals<T: RlstScalar<Real = T>>(cell: ReferenceCellType) -> V
     normals
 }
 
-/// The faces of the reference cell
-pub fn volumes(cell: ReferenceCellType) -> Vec<Vec<usize>> {
-    match cell {
-        ReferenceCellType::Point => vec![],
-        ReferenceCellType::Interval => vec![],
-        ReferenceCellType::Triangle => vec![],
-        ReferenceCellType::Quadrilateral => vec![],
-        ReferenceCellType::Tetrahedron => vec![vec![0, 1, 2, 3]],
-        ReferenceCellType::Hexahedron => vec![vec![0, 1, 2, 3, 4, 5, 6, 7]],
-        ReferenceCellType::Prism => vec![vec![0, 1, 2, 3, 4, 5]],
-        ReferenceCellType::Pyramid => vec![vec![0, 1, 2, 3, 4]],
-    }
-}
-
 /// The types of the subentities of the reference cell
 pub fn entity_types(cell: ReferenceCellType) -> Vec<Vec<ReferenceCellType>> {
     match cell {