diff --git a/README.md b/README.md
index 8fb87c3cd..502974ce3 100644
--- a/README.md
+++ b/README.md
@@ -151,6 +151,7 @@ The following elements are supported on a hexahedron:
   - [Bubble](https://defelement.com/elements/bubble.html)
   - [DPC](https://defelement.com/elements/dpc.html)
   - [Serendipity](https://defelement.com/elements/serendipity.html)
+  - [iso](https://defelement.com/elements/p1-iso-p2.html)
 
 
 ### Prism
diff --git a/cpp/basix/e-lagrange.cpp b/cpp/basix/e-lagrange.cpp
index c479a65ba..0d9dd5e1a 100644
--- a/cpp/basix/e-lagrange.cpp
+++ b/cpp/basix/e-lagrange.cpp
@@ -851,7 +851,7 @@ create_d_iso(cell::type celltype, int degree, element::lagrange_variant variant,
     _M(i, 0, i, 0) = 1.0;
 
   return FiniteElement(
-      element::family::P, celltype, polyset::type::macroedge, degree, {},
+      element::family::iso, celltype, polyset::type::macroedge, degree, {},
       impl::mdspan_t<const T, 2>(math::eye<T>(ndofs).data(), ndofs, ndofs),
       impl::to_mdspan(x), impl::to_mdspan(M), 0, maps::type::identity,
       sobolev::space::L2, true, degree, degree, variant,
@@ -1434,10 +1434,11 @@ FiniteElement<T> basix::element::create_iso(cell::type celltype, int degree,
                                             lagrange_variant variant,
                                             bool discontinuous)
 {
-  if (celltype != cell::type::interval && celltype != cell::type::quadrilateral)
+  if (celltype != cell::type::interval && celltype != cell::type::quadrilateral
+      && celltype != cell::type::hexahedron)
   {
     throw std::runtime_error("Can currently only create iso elements on "
-                             "intervals and quadrilaterals");
+                             "intervals, quadrilaterals, and hexahedra");
   }
 
   if (variant == lagrange_variant::unset)
@@ -1575,7 +1576,7 @@ FiniteElement<T> basix::element::create_iso(cell::type celltype, int degree,
   auto tensor_factors
       = create_tensor_product_factors<T>(celltype, degree, variant);
   return FiniteElement<T>(
-      family::P, celltype, polyset::type::macroedge, degree, {},
+      family::iso, celltype, polyset::type::macroedge, degree, {},
       impl::mdspan_t<T, 2>(math::eye<T>(ndofs).data(), ndofs, ndofs), xview,
       Mview, 0, maps::type::identity, space, discontinuous, degree, degree,
       variant, dpc_variant::unset, tensor_factors);
diff --git a/cpp/basix/polyset.cpp b/cpp/basix/polyset.cpp
index 726e94845..147d06cd1 100644
--- a/cpp/basix/polyset.cpp
+++ b/cpp/basix/polyset.cpp
@@ -227,8 +227,8 @@ void tabulate_polyset_line_macroedge_derivs(
 //-----------------------------------------------------------------------------
 
 /// Compute the complete set of derivatives from 0 to nderiv, for all
-/// the piecewise polynomials up to order n on a line segment split into two
-/// parts.
+/// the piecewise polynomials up to order n on a quadrilateral split into 4 by
+/// splitting each edge into so parts
 template <typename T>
 void tabulate_polyset_quadrilateral_macroedge_derivs(
     stdex::mdspan<T, stdex::dextents<std::size_t, 3>> P, std::size_t n,
@@ -282,7 +282,7 @@ void tabulate_polyset_quadrilateral_macroedge_derivs(
           value += factorials[k] * x_term;
         }
       }
-      for (std::size_t dy = 0; dy <= nderiv; ++dy)
+      for (std::size_t dy = 0; dy <= nderiv - dx; ++dy)
         for (std::size_t jy = 0; jy < 2 * n + 1; ++jy)
           P(idx(dx, dy), quad_idx(0, jy), p) = value;
     }
@@ -314,7 +314,7 @@ void tabulate_polyset_quadrilateral_macroedge_derivs(
           value *= pow(0.5 - x0[p], n - j - dx);
           for (std::size_t i = n - j; i > n - j - dx; --i)
             value *= -i;
-          for (std::size_t dy = 0; dy <= nderiv; ++dy)
+          for (std::size_t dy = 0; dy <= nderiv - dx; ++dy)
             for (std::size_t jy = 0; jy < 2 * n + 1; ++jy)
               P(idx(dx, dy), quad_idx(j + 1, jy), p) = value;
         }
@@ -331,7 +331,7 @@ void tabulate_polyset_quadrilateral_macroedge_derivs(
           value *= pow(x0[p] - 0.5, n - j - dx);
           for (std::size_t i = n - j; i > n - j - dx; --i)
             value *= i;
-          for (std::size_t dy = 0; dy <= nderiv; ++dy)
+          for (std::size_t dy = 0; dy <= nderiv - dx; ++dy)
             for (std::size_t jy = 0; jy < 2 * n + 1; ++jy)
               P(idx(dx, dy), quad_idx(j + n + 1, jy), p) = value;
         }
@@ -370,7 +370,7 @@ void tabulate_polyset_quadrilateral_macroedge_derivs(
           value += factorials[k] * y_term;
         }
       }
-      for (std::size_t dx = 0; dx <= nderiv; ++dx)
+      for (std::size_t dx = 0; dx <= nderiv - dy; ++dx)
         for (std::size_t jx = 0; jx < 2 * n + 1; ++jx)
           P(idx(dx, dy), quad_idx(jx, 0), p) *= value;
     }
@@ -402,12 +402,14 @@ void tabulate_polyset_quadrilateral_macroedge_derivs(
           value *= pow(0.5 - x1[p], n - j - dy);
           for (std::size_t i = n - j; i > n - j - dy; --i)
             value *= -i;
-          for (std::size_t dx = 0; dx <= nderiv; ++dx)
+          for (std::size_t dx = 0; dx <= nderiv - dy; ++dx)
+          {
             for (std::size_t jx = 0; jx < 2 * n + 1; ++jx)
             {
               P(idx(dx, dy), quad_idx(jx, j + 1), p) *= value;
               P(idx(dx, dy), quad_idx(jx, j + n + 1), p) *= 0.0;
             }
+          }
         }
         else
         {
@@ -422,7 +424,7 @@ void tabulate_polyset_quadrilateral_macroedge_derivs(
           value *= pow(x1[p] - 0.5, n - j - dy);
           for (std::size_t i = n - j; i > n - j - dy; --i)
             value *= i;
-          for (std::size_t dx = 0; dx <= nderiv; ++dx)
+          for (std::size_t dx = 0; dx <= nderiv - dy; ++dx)
           {
             for (std::size_t jx = 0; jx < 2 * n + 1; ++jx)
             {
@@ -436,6 +438,340 @@ void tabulate_polyset_quadrilateral_macroedge_derivs(
   }
 }
 //-----------------------------------------------------------------------------
+/// Compute the complete set of derivatives from 0 to nderiv, for all
+/// the piecewise polynomials up to order n on a hexahedron split into 4 by
+/// splitting each edge into so parts
+template <typename T>
+void tabulate_polyset_hexahedron_macroedge_derivs(
+    stdex::mdspan<T, stdex::dextents<std::size_t, 3>> P, std::size_t n,
+    std::size_t nderiv,
+    stdex::mdspan<const T, stdex::dextents<std::size_t, 2>> x)
+{
+  assert(x.extent(0) > 0);
+  assert(P.extent(0) == nderiv + 1);
+  assert(P.extent(1) == (2 * n + 1) * (2 * n + 1) * (2 * n + 1));
+  assert(P.extent(2) == x.extent(0));
+
+  // Indexing for hexahedral basis functions
+  auto hex_idx
+      = [n](std::size_t px, std::size_t py, std::size_t pz) -> std::size_t
+  { return (2 * n + 1) * (2 * n + 1) * px + (2 * n + 1) * py + pz; };
+
+  auto x0 = stdex::submdspan(x, stdex::full_extent, 0);
+  auto x1 = stdex::submdspan(x, stdex::full_extent, 1);
+  auto x2 = stdex::submdspan(x, stdex::full_extent, 2);
+
+  std::fill(P.data_handle(), P.data_handle() + P.size(), 0.0);
+
+  std::vector<T> factorials(n + 1, 0.0);
+
+  // Fill with values of polynomials in x
+  for (std::size_t k = 0; k <= n; ++k)
+  {
+    factorials[k] = (k % 2 == 0 ? 1 : -1) * single_choose(2 * n + 1 - k, n - k)
+                    * single_choose(n, k) * pow(2, n - k);
+  }
+  for (std::size_t dx = 0; dx <= nderiv; ++dx)
+  {
+    for (std::size_t p = 0; p < P.extent(2); ++p)
+    {
+      T value = 0.0;
+      if (x0[p] <= 0.5)
+      {
+        for (std::size_t k = 0; k + dx <= n; ++k)
+        {
+          T x_term = pow(x0[p], n - k - dx);
+          for (std::size_t i = n - k; i > n - k - dx; --i)
+            x_term *= i;
+          value += factorials[k] * x_term;
+        }
+      }
+      else
+      {
+        for (std::size_t k = 0; k + dx <= n; ++k)
+        {
+          T x_term = pow(1.0 - x0[p], n - k - dx);
+          for (std::size_t i = n - k; i > n - k - dx; --i)
+            x_term *= -i;
+          value += factorials[k] * x_term;
+        }
+      }
+      for (std::size_t dy = 0; dy <= nderiv - dx; ++dy)
+        for (std::size_t dz = 0; dz <= nderiv - dy - dx; ++dz)
+          for (std::size_t jy = 0; jy < 2 * n + 1; ++jy)
+            for (std::size_t jz = 0; jz < 2 * n + 1; ++jz)
+              P(idx(dx, dy, dz), hex_idx(0, jy, jz), p) = value;
+    }
+  }
+
+  for (std::size_t j = 0; j < n; ++j)
+  {
+    for (std::size_t k = 0; k <= j; ++k)
+    {
+      factorials[k] = (k % 2 == 0 ? 1 : -1)
+                      * single_choose(2 * n + 1 - k, j - k)
+                      * single_choose(j, k) * pow(2, j - k) * pow(2, n - j)
+                      * sqrt(4 * (n - j) + 2);
+    }
+    for (std::size_t dx = 0; dx <= nderiv; ++dx)
+    {
+      for (std::size_t p = 0; p < P.extent(2); ++p)
+      {
+        if (x0[p] <= 0.5)
+        {
+          T value = 0.0;
+          for (std::size_t k = 0; k + dx <= j; ++k)
+          {
+            T x_term = pow(x0[p], j - k - dx);
+            for (std::size_t i = j - k; i > j - k - dx; --i)
+              x_term *= i;
+            value += factorials[k] * x_term;
+          }
+          value *= pow(0.5 - x0[p], n - j - dx);
+          for (std::size_t i = n - j; i > n - j - dx; --i)
+            value *= -i;
+          for (std::size_t dy = 0; dy <= nderiv - dx; ++dy)
+            for (std::size_t dz = 0; dz <= nderiv - dy - dx; ++dz)
+              for (std::size_t jy = 0; jy < 2 * n + 1; ++jy)
+                for (std::size_t jz = 0; jz < 2 * n + 1; ++jz)
+                  P(idx(dx, dy, dz), hex_idx(j + 1, jy, jz), p) = value;
+        }
+        else
+        {
+          T value = 0.0;
+          for (std::size_t k = 0; k + dx <= j; ++k)
+          {
+            T x_term = pow(1.0 - x0[p], j - k - dx);
+            for (std::size_t i = j - k; i > j - k - dx; --i)
+              x_term *= -i;
+            value += factorials[k] * x_term;
+          }
+          value *= pow(x0[p] - 0.5, n - j - dx);
+          for (std::size_t i = n - j; i > n - j - dx; --i)
+            value *= i;
+          for (std::size_t dy = 0; dy <= nderiv - dx; ++dy)
+            for (std::size_t dz = 0; dz <= nderiv - dy - dx; ++dz)
+              for (std::size_t jy = 0; jy < 2 * n + 1; ++jy)
+                for (std::size_t jz = 0; jz < 2 * n + 1; ++jz)
+                  P(idx(dx, dy, dz), hex_idx(j + n + 1, jy, jz), p) = value;
+        }
+      }
+    }
+  }
+
+  // Multiply by values of polynomials in y
+  for (std::size_t k = 0; k <= n; ++k)
+  {
+    factorials[k] = (k % 2 == 0 ? 1 : -1) * single_choose(2 * n + 1 - k, n - k)
+                    * single_choose(n, k) * pow(2, n - k);
+  }
+  for (std::size_t dy = 0; dy <= nderiv; ++dy)
+  {
+    for (std::size_t p = 0; p < P.extent(2); ++p)
+    {
+      T value = 0.0;
+      if (x1[p] <= 0.5)
+      {
+        for (std::size_t k = 0; k + dy <= n; ++k)
+        {
+          T y_term = pow(x1[p], n - k - dy);
+          for (std::size_t i = n - k; i > n - k - dy; --i)
+            y_term *= i;
+          value += factorials[k] * y_term;
+        }
+      }
+      else
+      {
+        for (std::size_t k = 0; k + dy <= n; ++k)
+        {
+          T y_term = pow(1.0 - x1[p], n - k - dy);
+          for (std::size_t i = n - k; i > n - k - dy; --i)
+            y_term *= -i;
+          value += factorials[k] * y_term;
+        }
+      }
+      for (std::size_t dx = 0; dx <= nderiv - dy; ++dx)
+        for (std::size_t dz = 0; dz <= nderiv - dx - dy; ++dz)
+          for (std::size_t jx = 0; jx < 2 * n + 1; ++jx)
+            for (std::size_t jz = 0; jz < 2 * n + 1; ++jz)
+              P(idx(dx, dy, dz), hex_idx(jx, 0, jz), p) *= value;
+    }
+  }
+
+  for (std::size_t j = 0; j < n; ++j)
+  {
+    for (std::size_t k = 0; k <= j; ++k)
+    {
+      factorials[k] = (k % 2 == 0 ? 1 : -1)
+                      * single_choose(2 * n + 1 - k, j - k)
+                      * single_choose(j, k) * pow(2, j - k) * pow(2, n - j)
+                      * sqrt(4 * (n - j) + 2);
+    }
+    for (std::size_t dy = 0; dy <= nderiv; ++dy)
+    {
+      for (std::size_t p = 0; p < P.extent(2); ++p)
+      {
+        if (x1[p] <= 0.5)
+        {
+          T value = 0.0;
+          for (std::size_t k = 0; k + dy <= j; ++k)
+          {
+            T y_term = pow(x1[p], j - k - dy);
+            for (std::size_t i = j - k; i > j - k - dy; --i)
+              y_term *= i;
+            value += factorials[k] * y_term;
+          }
+          value *= pow(0.5 - x1[p], n - j - dy);
+          for (std::size_t i = n - j; i > n - j - dy; --i)
+            value *= -i;
+          for (std::size_t dx = 0; dx <= nderiv - dy; ++dx)
+            for (std::size_t dz = 0; dz <= nderiv - dx - dy; ++dz)
+              for (std::size_t jx = 0; jx < 2 * n + 1; ++jx)
+                for (std::size_t jz = 0; jz < 2 * n + 1; ++jz)
+                {
+                  P(idx(dx, dy, dz), hex_idx(jx, j + 1, jz), p) *= value;
+                  P(idx(dx, dy, dz), hex_idx(jx, j + n + 1, jz), p) *= 0.0;
+                }
+        }
+        else
+        {
+          T value = 0.0;
+          for (std::size_t k = 0; k + dy <= j; ++k)
+          {
+            T y_term = pow(1.0 - x1[p], j - k - dy);
+            for (std::size_t i = j - k; i > j - k - dy; --i)
+              y_term *= -i;
+            value += factorials[k] * y_term;
+          }
+          value *= pow(x1[p] - 0.5, n - j - dy);
+          for (std::size_t i = n - j; i > n - j - dy; --i)
+            value *= i;
+          for (std::size_t dx = 0; dx <= nderiv - dy; ++dx)
+          {
+            for (std::size_t dz = 0; dz <= nderiv - dx - dy; ++dz)
+            {
+              for (std::size_t jx = 0; jx < 2 * n + 1; ++jx)
+              {
+                for (std::size_t jz = 0; jz < 2 * n + 1; ++jz)
+                {
+                  P(idx(dx, dy, dz), hex_idx(jx, j + 1, jz), p) *= 0.0;
+                  P(idx(dx, dy, dz), hex_idx(jx, j + n + 1, jz), p) *= value;
+                }
+              }
+            }
+          }
+        }
+      }
+    }
+  }
+
+  // Multiply by values of polynomials in z
+  for (std::size_t k = 0; k <= n; ++k)
+  {
+    factorials[k] = (k % 2 == 0 ? 1 : -1) * single_choose(2 * n + 1 - k, n - k)
+                    * single_choose(n, k) * pow(2, n - k);
+  }
+  for (std::size_t dz = 0; dz <= nderiv; ++dz)
+  {
+    for (std::size_t p = 0; p < P.extent(2); ++p)
+    {
+      T value = 0.0;
+      if (x2[p] <= 0.5)
+      {
+        for (std::size_t k = 0; k + dz <= n; ++k)
+        {
+          T z_term = pow(x2[p], n - k - dz);
+          for (std::size_t i = n - k; i > n - k - dz; --i)
+            z_term *= i;
+          value += factorials[k] * z_term;
+        }
+      }
+      else
+      {
+        for (std::size_t k = 0; k + dz <= n; ++k)
+        {
+          T z_term = pow(1.0 - x2[p], n - k - dz);
+          for (std::size_t i = n - k; i > n - k - dz; --i)
+            z_term *= -i;
+          value += factorials[k] * z_term;
+        }
+      }
+      for (std::size_t dx = 0; dx <= nderiv - dz; ++dx)
+        for (std::size_t dy = 0; dy <= nderiv - dx - dz; ++dy)
+          for (std::size_t jx = 0; jx < 2 * n + 1; ++jx)
+            for (std::size_t jy = 0; jy < 2 * n + 1; ++jy)
+              P(idx(dx, dy, dz), hex_idx(jx, jy, 0), p) *= value;
+    }
+  }
+
+  for (std::size_t j = 0; j < n; ++j)
+  {
+    for (std::size_t k = 0; k <= j; ++k)
+    {
+      factorials[k] = (k % 2 == 0 ? 1 : -1)
+                      * single_choose(2 * n + 1 - k, j - k)
+                      * single_choose(j, k) * pow(2, j - k) * pow(2, n - j)
+                      * sqrt(4 * (n - j) + 2);
+    }
+    for (std::size_t dz = 0; dz <= nderiv; ++dz)
+    {
+      for (std::size_t p = 0; p < P.extent(2); ++p)
+      {
+        if (x2[p] <= 0.5)
+        {
+          T value = 0.0;
+          for (std::size_t k = 0; k + dz <= j; ++k)
+          {
+            T z_term = pow(x2[p], j - k - dz);
+            for (std::size_t i = j - k; i > j - k - dz; --i)
+              z_term *= i;
+            value += factorials[k] * z_term;
+          }
+          value *= pow(0.5 - x2[p], n - j - dz);
+          for (std::size_t i = n - j; i > n - j - dz; --i)
+            value *= -i;
+          for (std::size_t dx = 0; dx <= nderiv - dz; ++dx)
+            for (std::size_t dy = 0; dy <= nderiv - dx - dz; ++dy)
+              for (std::size_t jx = 0; jx < 2 * n + 1; ++jx)
+                for (std::size_t jy = 0; jy < 2 * n + 1; ++jy)
+                {
+                  P(idx(dx, dy, dz), hex_idx(jx, jy, j + 1), p) *= value;
+                  P(idx(dx, dy, dz), hex_idx(jx, jy, j + n + 1), p) *= 0.0;
+                }
+        }
+        else
+        {
+          T value = 0.0;
+          for (std::size_t k = 0; k + dz <= j; ++k)
+          {
+            T z_term = pow(1.0 - x2[p], j - k - dz);
+            for (std::size_t i = j - k; i > j - k - dz; --i)
+              z_term *= -i;
+            value += factorials[k] * z_term;
+          }
+          value *= pow(x2[p] - 0.5, n - j - dz);
+          for (std::size_t i = n - j; i > n - j - dz; --i)
+            value *= i;
+          for (std::size_t dx = 0; dx <= nderiv - dz; ++dx)
+          {
+            for (std::size_t dy = 0; dy <= nderiv - dx - dz; ++dy)
+            {
+              for (std::size_t jx = 0; jx < 2 * n + 1; ++jx)
+              {
+                for (std::size_t jy = 0; jy < 2 * n + 1; ++jy)
+                {
+                  P(idx(dx, dy, dz), hex_idx(jx, jy, j + 1), p) *= 0.0;
+                  P(idx(dx, dy, dz), hex_idx(jx, jy, j + n + 1), p) *= value;
+                }
+              }
+            }
+          }
+        }
+      }
+    }
+  }
+}
+//-----------------------------------------------------------------------------
 
 /// Compute the complete set of derivatives from 0 to nderiv, for all
 /// the polynomials up to order n on a triangle in [0, 1][0, 1]. The
@@ -1854,6 +2190,9 @@ void polyset::tabulate(
     case cell::type::quadrilateral:
       tabulate_polyset_quadrilateral_macroedge_derivs(P, d, n, x);
       return;
+    case cell::type::hexahedron:
+      tabulate_polyset_hexahedron_macroedge_derivs(P, d, n, x);
+      return;
     default:
       throw std::runtime_error("Polynomial set: unsupported cell type");
     }
@@ -1926,6 +2265,8 @@ int polyset::dim(cell::type celltype, polyset::type ptype, int d)
       return 2 * d + 1;
     case cell::type::quadrilateral:
       return (2 * d + 1) * (2 * d + 1);
+    case cell::type::hexahedron:
+      return (2 * d + 1) * (2 * d + 1) * (2 * d + 1);
     default:
       return 1;
     }
diff --git a/cpp/basix/quadrature.cpp b/cpp/basix/quadrature.cpp
index 5bbb6fb28..2cf4d4ee8 100644
--- a/cpp/basix/quadrature.cpp
+++ b/cpp/basix/quadrature.cpp
@@ -4760,6 +4760,52 @@ make_macroedge_quadrature(quadrature::type rule, cell::type celltype, int m)
     }
     return {std::move(x), std::move(w)};
   }
+  case cell::type::hexahedron:
+  {
+    if (m == 0)
+    {
+      return standard_q;
+    }
+    const std::size_t npts = standard_q[0].size() / 3;
+    std::vector<T> x(npts * 24);
+    std::vector<T> w(npts * 8);
+    for (std::size_t i = 0; i < npts; ++i)
+    {
+      x[3 * i] = 0.5 * standard_q[0][3 * i];
+      x[3 * i + 1] = 0.5 * standard_q[0][3 * i + 1];
+      x[3 * i + 2] = 0.5 * standard_q[0][3 * i + 2];
+      x[3 * (npts + i)] = 0.5 + 0.5 * standard_q[0][3 * i];
+      x[3 * (npts + i) + 1] = 0.5 * standard_q[0][3 * i + 1];
+      x[3 * (npts + i) + 2] = 0.5 * standard_q[0][3 * i + 2];
+      x[3 * (2 * npts + i)] = 0.5 * standard_q[0][3 * i];
+      x[3 * (2 * npts + i) + 1] = 0.5 + 0.5 * standard_q[0][3 * i + 1];
+      x[3 * (2 * npts + i) + 2] = 0.5 * standard_q[0][3 * i + 2];
+      x[3 * (3 * npts + i)] = 0.5 + 0.5 * standard_q[0][3 * i];
+      x[3 * (3 * npts + i) + 1] = 0.5 + 0.5 * standard_q[0][3 * i + 1];
+      x[3 * (3 * npts + i) + 2] = 0.5 * standard_q[0][3 * i + 2];
+      x[3 * (4 * npts + i)] = 0.5 * standard_q[0][3 * i];
+      x[3 * (4 * npts + i) + 1] = 0.5 * standard_q[0][3 * i + 1];
+      x[3 * (4 * npts + i) + 2] = 0.5 + 0.5 * standard_q[0][3 * i + 2];
+      x[3 * (5 * npts + i)] = 0.5 + 0.5 * standard_q[0][3 * i];
+      x[3 * (5 * npts + i) + 1] = 0.5 * standard_q[0][3 * i + 1];
+      x[3 * (5 * npts + i) + 2] = 0.5 + 0.5 * standard_q[0][3 * i + 2];
+      x[3 * (6 * npts + i)] = 0.5 * standard_q[0][3 * i];
+      x[3 * (6 * npts + i) + 1] = 0.5 + 0.5 * standard_q[0][3 * i + 1];
+      x[3 * (6 * npts + i) + 2] = 0.5 + 0.5 * standard_q[0][3 * i + 2];
+      x[3 * (7 * npts + i)] = 0.5 + 0.5 * standard_q[0][3 * i];
+      x[3 * (7 * npts + i) + 1] = 0.5 + 0.5 * standard_q[0][3 * i + 1];
+      x[3 * (7 * npts + i) + 2] = 0.5 + 0.5 * standard_q[0][3 * i + 2];
+      w[i] = 0.125 * standard_q[1][i];
+      w[npts + i] = 0.125 * standard_q[1][i];
+      w[2 * npts + i] = 0.125 * standard_q[1][i];
+      w[3 * npts + i] = 0.125 * standard_q[1][i];
+      w[4 * npts + i] = 0.125 * standard_q[1][i];
+      w[5 * npts + i] = 0.125 * standard_q[1][i];
+      w[6 * npts + i] = 0.125 * standard_q[1][i];
+      w[7 * npts + i] = 0.125 * standard_q[1][i];
+    }
+    return {std::move(x), std::move(w)};
+  }
   default:
     throw std::runtime_error("Macro quadrature not supported on this cell.");
   }
diff --git a/test/test_iso.py b/test/test_iso.py
index fd00d75a0..b97c72925 100644
--- a/test/test_iso.py
+++ b/test/test_iso.py
@@ -4,7 +4,9 @@
 
 
 @pytest.mark.parametrize("degree", range(1, 5))
-@pytest.mark.parametrize("cell", [basix.CellType.interval, basix.CellType.quadrilateral])
+@pytest.mark.parametrize("cell", [
+    basix.CellType.interval, basix.CellType.quadrilateral,
+    basix.CellType.hexahedron])
 def test_iso_element(degree, cell):
     e = basix.create_element(basix.ElementFamily.iso, cell, degree, basix.LagrangeVariant.gll_warped)
     e2 = basix.create_element(basix.ElementFamily.P, cell, 2 * degree, basix.LagrangeVariant.gll_warped)
diff --git a/test/test_orthonormal_basis.py b/test/test_orthonormal_basis.py
index 4cc2e898d..81466cf04 100644
--- a/test/test_orthonormal_basis.py
+++ b/test/test_orthonormal_basis.py
@@ -121,6 +121,7 @@ def test_standard(cell_type, order):
 @pytest.mark.parametrize("cell_type", [
     basix.CellType.interval,
     basix.CellType.quadrilateral,
+    basix.CellType.hexahedron,
 ])
 @pytest.mark.parametrize("order", [0, 1, 2, 3, 4])
 def test_macroedge(cell_type, order):
@@ -133,5 +134,4 @@ def test_macroedge(cell_type, order):
     for i in range(ndofs):
         for j in range(ndofs):
             mat[i, j] = sum(basis[i, :] * basis[j, :] * Qwts)
-
     assert np.allclose(mat, np.eye(ndofs))