diff --git a/.github/workflows/dolfin-tests.yml b/.github/workflows/dolfin-tests.yml
index d2896b841..9bdb5bb23 100644
--- a/.github/workflows/dolfin-tests.yml
+++ b/.github/workflows/dolfin-tests.yml
@@ -51,7 +51,7 @@ jobs:
         if: github.event_name != 'workflow_dispatch'
         run: |
           python3 -m pip install git+https://github.com/FEniCS/ufl.git
-          python3 -m pip install git+https://github.com/FEniCS/ffcx.git
+          python3 -m pip install git+https://github.com/FEniCS/ffcx.git@mscroggs/make_quadrature
       - name: Install FEniCS Python components
         if: github.event_name == 'workflow_dispatch'
         run: |
@@ -64,7 +64,7 @@ jobs:
         with:
           path: ./dolfinx
           repository: FEniCS/dolfinx
-          ref: main
+          ref: mscroggs/make_quadrature
       - name: Get DOLFINx
         if: github.event_name == 'workflow_dispatch'
         uses: actions/checkout@v3
diff --git a/.github/workflows/ffcx-tests.yml b/.github/workflows/ffcx-tests.yml
index a5b0d3444..d5bed8d93 100644
--- a/.github/workflows/ffcx-tests.yml
+++ b/.github/workflows/ffcx-tests.yml
@@ -54,7 +54,7 @@ jobs:
         with:
           path: ./ffcx
           repository: FEniCS/ffcx
-          ref: main
+          ref: mscroggs/make_quadrature
       - name: Get FFCx source (specified branch)
         if: github.event_name == 'workflow_dispatch'
         uses: actions/checkout@v3
diff --git a/README.md b/README.md
index f4a302b8f..060d9c6a2 100644
--- a/README.md
+++ b/README.md
@@ -71,7 +71,7 @@ The following elements are supported on an interval:
   - [Bubble](https://defelement.com/elements/bubble.html)
   - [Serendipity](https://defelement.com/elements/serendipity.html)
   - [Hermite](https://defelement.com/elements/hermite.html)
-
+  - [iso](https://defelement.com/elements/p1-iso-p2.html)
 
 ### Triangle
 
diff --git a/cpp/basix/docs.h b/cpp/basix/docs.h
index c9b7e8d1c..0ffe007b8 100644
--- a/cpp/basix/docs.h
+++ b/cpp/basix/docs.h
@@ -54,7 +54,7 @@ Sub-entity of a cell, given by topological dimension and index
 )";
 
 const std::string create_lattice__celltype_n_type_exterior = R"(
-@brief Create a lattice of points on a reference cell optionally
+Create a lattice of points on a reference cell optionally
 including the outer surface points.
 
 For a given `celltype`, this creates a set of points on a regular
@@ -80,7 +80,7 @@ type::equispaced when `n < 3`.
 )";
 
 const std::string create_lattice__celltype_n_type_exterior_method = R"(
-@brief Create a lattice of points on a reference cell optionally
+Create a lattice of points on a reference cell optionally
 including the outer surface points.
 
 For a given `celltype`, this creates a set of points on a regular
@@ -169,7 +169,7 @@ Get the jacobians of the facets of a reference cell
 )";
 
 const std::string FiniteElement__tabulate = R"(
-@brief Compute basis values and derivatives at set of points.
+Compute basis values and derivatives at set of points.
 
 NOTE: The version of `FiniteElement::tabulate` with the basis data
 as an out argument should be preferred for repeated call where
@@ -270,7 +270,7 @@ performance is critical.
 )";
 
 const std::string FiniteElement__base_transformations = R"(
-@brief Get the base transformations.
+Get the base transformations.
 
 The base transformations represent the effect of rotating or reflecting
 a subentity of the cell on the numbering and orientation of the DOFs.
@@ -409,6 +409,7 @@ Create a custom finite element
     discontinuous: Indicates whether or not this is the discontinuous version of the element
     highest_complete_degree: The highest degree n such that a Lagrange (or vector Lagrange) element of degree n is a subspace of this element
     highest_degree: The degree of a polynomial in this element's polyset
+    poly_type: The type of polyset to use for this element
 
 Returns:
     A custom finite element
@@ -434,30 +435,30 @@ Create an element using a given Lagrange variant and a given DPC variant
 )";
 
 const std::string compute_interpolation_operator = R"(
-Computes a matrix that represents the interpolation between two
-elements.
+Compute a matrix that represents the interpolation between
+two elements.
 
 If the two elements have the same value size, this function returns
 the interpolation between them.
 
-If element_from has value size 1 and element_to has value size > 1, then
-this function returns a matrix to interpolate from a blocked element_from
-(ie multiple copies of element_from) into element_to.
+If element_from has value size 1 and element_to has value size > 1,
+then this function returns a matrix to interpolate from a blocked
+element_from (ie multiple copies of element_from) into element_to.
 
-If element_to has value size 1 and element_from has value size > 1, then
-this function returns a matrix that interpolates the components of
-element_from into copies of element_to.
+If element_to has value size 1 and element_from has value size > 1,
+then this function returns a matrix that interpolates the components
+of element_from into copies of element_to.
 
 NOTE: If the elements have different value sizes and both are
 greater than 1, this function throws a runtime error
 
-In order to interpolate functions between finite element spaces on arbitrary
-cells, the functions must be pulled back to the reference element (this pull
-back includes applying DOF transformations). The matrix that this function
-returns can then be applied, then the result pushed forward to the cell. If
-element_from and element_to have the same map type, then only the DOF
-transformations need to be applied, as the pull back and push forward cancel
-each other out.
+In order to interpolate functions between finite element spaces on
+arbitrary cells, the functions must be pulled back to the reference
+element (this pull back includes applying DOF transformations). The
+matrix that this function returns can then be applied, then the
+result pushed forward to the cell. If element_from and element_to
+have the same map type, then only the DOF transformations need to be
+applied, as the pull back and push forward cancel each other out.
 
 Args:
     element_from: The element to interpolate from
@@ -470,7 +471,7 @@ each other out.
 )";
 
 const std::string tabulate_polynomial_set = R"(
-@brief Tabulate the orthonormal polynomial basis, and derivatives,
+Tabulate the orthonormal polynomial basis, and derivatives,
 at points on the reference cell.
 
 All derivatives up to the given order are computed. If derivatives
@@ -486,6 +487,7 @@ there are `(nderiv + 1)(nderiv + 2)(nderiv + 3)/6`. The ordering is
 
 Args:
     celltype: Cell type
+    ptype: The polynomial type
     d: Polynomial degree
     n: Maximum derivative order. Use n = 0 for the basis only.
     x: Points at which to evaluate the basis. The shape is (number of points, geometric dimension).
@@ -512,7 +514,7 @@ there are `(nderiv + 1)(nderiv + 2)(nderiv + 3)/6`. The ordering is
 )";
 
 const std::string tabulate_polynomials = R"(
-@brief Tabulate a set of polynomials.
+Tabulate a set of polynomials.
 
 Args:
     polytype: Polynomial type
@@ -526,7 +528,7 @@ const std::string tabulate_polynomials = R"(
 )";
 
 const std::string polynomials_dim = R"(
-@brief Dimension of a polynomial space.
+Dimension of a polynomial space.
 
 Args:
     polytype: The polynomial type
@@ -538,24 +540,13 @@ const std::string polynomials_dim = R"(
     polynomial degree @p d
 )";
 
-const std::string make_quadrature__rule_celltype_m = R"(
+const std::string make_quadrature__rule_celltype_polytype_m = R"(
 Make a quadrature rule on a reference cell
 
 Args:
     rule: Type of quadrature rule (or use quadrature::Default)
     celltype: The cell type
-    m: Maximum degree of polynomial that this quadrature rule will integrate exactly
-
-Returns:
-    List of points and list of weights. The number of points
-    arrays has shape (num points, gdim)
-)";
-
-const std::string make_quadrature__celltype_m = R"(
-Make a default quadrature rule on reference cell
-
-Args:
-    celltype: The cell type
+    polytype: The polyset type
     m: Maximum degree of polynomial that this quadrature rule will integrate exactly
 
 Returns:
@@ -574,9 +565,9 @@ Compute trivial indexing in a 1D array (for completeness)
 )";
 
 const std::string index__p_q = R"(
-Compute indexing in a 2D triangular array compressed into a 1D array.
-This can be used to find the index of a derivative returned by
-`FiniteElement::tabulate`. For instance to find d2N/dx2, use
+Compute indexing in a 2D triangular array compressed into a 1D
+array. This can be used to find the index of a derivative returned
+by `FiniteElement::tabulate`. For instance to find d2N/dx2, use
 `FiniteElement::tabulate(2, points)[idx(2, 0)];`
 
 Args:
@@ -588,7 +579,8 @@ This can be used to find the index of a derivative returned by
 )";
 
 const std::string index__p_q_r = R"(
-Compute indexing in a 3D tetrahedral array compressed into a 1D array
+Compute indexing in a 3D tetrahedral array compressed into a 1D
+array
 
 Args:
     p: Index in x
@@ -610,4 +602,32 @@ Get the intersection of two Sobolev spaces.
     Intersection of the spaces
 )";
 
+const std::string superset = R"(
+Get the polyset types that is a superset of two types on the given
+cell
+
+Args:
+    cell: The cell type
+    type1: The first polyset type
+    type2: The second polyset type
+
+Returns::
+    The superset type
+)";
+
+const std::string restriction = R"(
+Get the polyset type that represents the restrictions of a type on a
+subentity
+
+Args:
+    ptype: The polyset type
+    cell: The cell type
+    restriction_cell: The cell type of the subentity
+
+Returns::
+    The restricted polyset type
+)";
+
+
+
 } // namespace basix::docstring
diff --git a/cpp/basix/dof-transformations.cpp b/cpp/basix/dof-transformations.cpp
index 8c9af1fc1..385fd6d70 100644
--- a/cpp/basix/dof-transformations.cpp
+++ b/cpp/basix/dof-transformations.cpp
@@ -327,7 +327,7 @@ std::pair<std::vector<T>, std::array<std::size_t, 2>> compute_transformation(
 
   const std::size_t ndofs = imat.extent(0);
   const std::size_t npts = pts.extent(0);
-  const int psize = polyset::dim(cell_type, degree);
+  const int psize = polyset::dim(cell_type, polyset::type::standard, degree);
 
   std::size_t dofstart = 0;
   for (int d = 0; d < tdim; ++d)
@@ -352,7 +352,7 @@ std::pair<std::vector<T>, std::array<std::size_t, 2>> compute_transformation(
   }
 
   auto [polyset_vals_b, polyset_shape] = polyset::tabulate(
-      cell_type, degree, 0,
+      cell_type, polyset::type::standard, degree, 0,
       mdspan_t<const T, 2>(mapped_pts.data(), mapped_pts.extents()));
   assert(polyset_shape[0] == 1);
   mdspan_t<const T, 2> polyset_vals(polyset_vals_b.data(), polyset_shape[1],
diff --git a/cpp/basix/e-brezzi-douglas-marini.cpp b/cpp/basix/e-brezzi-douglas-marini.cpp
index d16c966ed..2272ebb9e 100644
--- a/cpp/basix/e-brezzi-douglas-marini.cpp
+++ b/cpp/basix/e-brezzi-douglas-marini.cpp
@@ -41,7 +41,8 @@ FiniteElement<T> element::create_bdm(cell::type celltype, int degree,
       = create_lagrange<T>(facettype, degree, lvariant, true);
   {
     auto [_x, xshape, _M, Mshape] = moments::make_normal_integral_moments<T>(
-        facet_moment_space, celltype, tdim, degree * 2);
+        facet_moment_space, celltype, polyset::type::standard, tdim,
+        degree * 2);
     assert(_x.size() == _M.size());
     for (std::size_t i = 0; i < _x.size(); ++i)
     {
@@ -56,8 +57,8 @@ FiniteElement<T> element::create_bdm(cell::type celltype, int degree,
   {
     // Interior integral moment
     auto [_x, xshape, _M, Mshape] = moments::make_dot_integral_moments<T>(
-        create_nedelec<T>(celltype, degree - 1, lvariant, true), celltype, tdim,
-        2 * degree - 1);
+        create_nedelec<T>(celltype, degree - 1, lvariant, true), celltype,
+        polyset::type::standard, tdim, 2 * degree - 1);
     assert(_x.size() == _M.size());
     for (std::size_t i = 0; i < _x.size(); ++i)
     {
@@ -90,11 +91,12 @@ FiniteElement<T> element::create_bdm(cell::type celltype, int degree,
   }
 
   // The number of order (degree) scalar polynomials
-  const std::size_t ndofs = tdim * polyset::dim(celltype, degree);
+  const std::size_t ndofs
+      = tdim * polyset::dim(celltype, polyset::type::standard, degree);
   sobolev::space space
       = discontinuous ? sobolev::space::L2 : sobolev::space::HDiv;
   return FiniteElement<T>(
-      family::BDM, celltype, degree, {tdim},
+      family::BDM, celltype, polyset::type::standard, degree, {tdim},
       impl::mdspan_t<T, 2>(math::eye<T>(ndofs).data(), ndofs, ndofs), xview,
       Mview, 0, maps::type::contravariantPiola, space, discontinuous, degree,
       degree, lvariant, element::dpc_variant::unset);
diff --git a/cpp/basix/e-bubble.cpp b/cpp/basix/e-bubble.cpp
index 474e0dcc3..450ff2cb2 100644
--- a/cpp/basix/e-bubble.cpp
+++ b/cpp/basix/e-bubble.cpp
@@ -67,10 +67,11 @@ FiniteElement<T> basix::element::create_bubble(cell::type celltype, int degree,
 
   // Evaluate the expansion polynomials at the quadrature points
   const auto [_pts, wts] = quadrature::make_quadrature<T>(
-      quadrature::type::Default, celltype, 2 * degree);
+      quadrature::type::Default, celltype, polyset::type::standard, 2 * degree);
   impl::mdspan_t<const T, 2> pts(_pts.data(), wts.size(),
                                  _pts.size() / wts.size());
-  const auto [_phi, shape] = polyset::tabulate(celltype, degree, 0, pts);
+  const auto [_phi, shape]
+      = polyset::tabulate(celltype, polyset::type::standard, degree, 0, pts);
   impl::mdspan_t<const T, 3> phi(_phi.data(), shape);
 
   // The number of order (degree) polynomials
@@ -103,7 +104,8 @@ FiniteElement<T> basix::element::create_bubble(cell::type celltype, int degree,
   {
   case cell::type::interval:
   {
-    const auto [_phi1, shape] = polyset::tabulate(celltype, degree - 2, 0, pts);
+    const auto [_phi1, shape] = polyset::tabulate(
+        celltype, polyset::type::standard, degree - 2, 0, pts);
     impl::mdspan_t<const T, 3> p1(_phi1.data(), shape);
     phi1 = create_phi1(p1, phi1_buffer);
     for (std::size_t i = 0; i < pts.extent(0); ++i)
@@ -115,7 +117,8 @@ FiniteElement<T> basix::element::create_bubble(cell::type celltype, int degree,
   }
   case cell::type::triangle:
   {
-    const auto [_phi1, shape] = polyset::tabulate(celltype, degree - 3, 0, pts);
+    const auto [_phi1, shape] = polyset::tabulate(
+        celltype, polyset::type::standard, degree - 3, 0, pts);
     impl::mdspan_t<const T, 3> p1(_phi1.data(), shape);
     phi1 = create_phi1(p1, phi1_buffer);
     for (std::size_t i = 0; i < pts.extent(0); ++i)
@@ -128,7 +131,8 @@ FiniteElement<T> basix::element::create_bubble(cell::type celltype, int degree,
   }
   case cell::type::tetrahedron:
   {
-    const auto [_phi1, shape] = polyset::tabulate(celltype, degree - 4, 0, pts);
+    const auto [_phi1, shape] = polyset::tabulate(
+        celltype, polyset::type::standard, degree - 4, 0, pts);
     impl::mdspan_t<const T, 3> p1(_phi1.data(), shape);
     phi1 = create_phi1(p1, phi1_buffer);
     for (std::size_t i = 0; i < pts.extent(0); ++i)
@@ -142,7 +146,8 @@ FiniteElement<T> basix::element::create_bubble(cell::type celltype, int degree,
   }
   case cell::type::quadrilateral:
   {
-    const auto [_phi1, shape] = polyset::tabulate(celltype, degree - 2, 0, pts);
+    const auto [_phi1, shape] = polyset::tabulate(
+        celltype, polyset::type::standard, degree - 2, 0, pts);
     impl::mdspan_t<const T, 3> p1(_phi1.data(), shape);
     phi1 = create_phi1(p1, phi1_buffer);
     for (std::size_t i = 0; i < pts.extent(0); ++i)
@@ -155,7 +160,8 @@ FiniteElement<T> basix::element::create_bubble(cell::type celltype, int degree,
   }
   case cell::type::hexahedron:
   {
-    const auto [_phi1, shape] = polyset::tabulate(celltype, degree - 2, 0, pts);
+    const auto [_phi1, shape] = polyset::tabulate(
+        celltype, polyset::type::standard, degree - 2, 0, pts);
     impl::mdspan_t<const T, 3> p1(_phi1.data(), shape);
     phi1 = create_phi1(p1, phi1_buffer);
     for (std::size_t i = 0; i < pts.extent(0); ++i)
@@ -185,9 +191,10 @@ FiniteElement<T> basix::element::create_bubble(cell::type celltype, int degree,
   sobolev::space space
       = discontinuous ? sobolev::space::L2 : sobolev::space::H1;
   return FiniteElement<T>(
-      element::family::bubble, celltype, degree, {}, wview, impl::to_mdspan(x),
-      impl::to_mdspan(M), 0, maps::type::identity, space, discontinuous, -1,
-      degree, element::lagrange_variant::unset, element::dpc_variant::unset);
+      element::family::bubble, celltype, polyset::type::standard, degree, {},
+      wview, impl::to_mdspan(x), impl::to_mdspan(M), 0, maps::type::identity,
+      space, discontinuous, -1, degree, element::lagrange_variant::unset,
+      element::dpc_variant::unset);
 }
 //-----------------------------------------------------------------------------
 template FiniteElement<float> element::create_bubble(cell::type, int, bool);
diff --git a/cpp/basix/e-crouzeix-raviart.cpp b/cpp/basix/e-crouzeix-raviart.cpp
index a795ab50c..8a534fdb7 100644
--- a/cpp/basix/e-crouzeix-raviart.cpp
+++ b/cpp/basix/e-crouzeix-raviart.cpp
@@ -88,7 +88,7 @@ FiniteElement<T> basix::element::create_cr(cell::type celltype, int degree,
   }
 
   return FiniteElement<T>(
-      element::family::CR, celltype, 1, {},
+      element::family::CR, celltype, polyset::type::standard, 1, {},
       impl::mdspan_t<T, 2>(math::eye<T>(ndofs).data(), ndofs, ndofs), xview,
       Mview, 0, maps::type::identity, sobolev::space::L2, discontinuous, degree,
       degree, element::lagrange_variant::unset, element::dpc_variant::unset);
diff --git a/cpp/basix/e-hermite.cpp b/cpp/basix/e-hermite.cpp
index 7c923ca96..5934187b5 100644
--- a/cpp/basix/e-hermite.cpp
+++ b/cpp/basix/e-hermite.cpp
@@ -31,7 +31,8 @@ FiniteElement<T> basix::element::create_hermite(cell::type celltype, int degree,
     throw std::runtime_error("Hermite element must have degree 3");
 
   const std::size_t tdim = cell::topological_dimension(celltype);
-  const std::size_t ndofs = polyset::dim(celltype, degree);
+  const std::size_t ndofs
+      = polyset::dim(celltype, polyset::type::standard, degree);
   const std::vector<std::vector<std::vector<int>>> topology
       = cell::topology(celltype);
 
@@ -99,7 +100,7 @@ FiniteElement<T> basix::element::create_hermite(cell::type celltype, int degree,
   sobolev::space space
       = discontinuous ? sobolev::space::L2 : sobolev::space::H2;
   return FiniteElement<T>(
-      element::family::Hermite, celltype, degree, {},
+      element::family::Hermite, celltype, polyset::type::standard, degree, {},
       impl::mdspan_t<T, 2>(math::eye<T>(ndofs).data(), ndofs, ndofs), xview,
       Mview, 1, maps::type::identity, space, discontinuous, -1, degree,
       element::lagrange_variant::unset, element::dpc_variant::unset);
diff --git a/cpp/basix/e-hhj.cpp b/cpp/basix/e-hhj.cpp
index 3dae5d50d..16ec6d3e7 100644
--- a/cpp/basix/e-hhj.cpp
+++ b/cpp/basix/e-hhj.cpp
@@ -24,7 +24,8 @@ FiniteElement<T> basix::element::create_hhj(cell::type celltype, int degree,
   const std::size_t tdim = cell::topological_dimension(celltype);
 
   const int nc = tdim * (tdim + 1) / 2;
-  const int basis_size = polyset::dim(celltype, degree);
+  const int basis_size
+      = polyset::dim(celltype, polyset::type::standard, degree);
   const std::size_t ndofs = basis_size * nc;
   const std::size_t psize = basis_size * tdim * tdim;
 
@@ -84,9 +85,10 @@ FiniteElement<T> basix::element::create_hhj(cell::type celltype, int degree,
         // Tabulate points in lattice
         cell::type ct = cell::sub_entity_type(celltype, d, e);
 
-        const std::size_t ndofs = polyset::dim(ct, degree + 1 - d);
-        const auto [ptsbuffer, wts]
-            = quadrature::make_quadrature<T>(ct, degree + (degree + 1 - d));
+        const std::size_t ndofs
+            = polyset::dim(ct, polyset::type::standard, degree + 1 - d);
+        const auto [ptsbuffer, wts] = quadrature::make_quadrature<T>(
+            quadrature::type::Default, ct, polyset::type::standard, degree + (degree + 1 - d));
         impl::mdspan_t<const T, 2> pts(ptsbuffer.data(), wts.size(),
                                        ptsbuffer.size() / wts.size());
 
@@ -177,10 +179,11 @@ FiniteElement<T> basix::element::create_hhj(cell::type celltype, int degree,
   sobolev::space space
       = discontinuous ? sobolev::space::L2 : sobolev::space::HDivDiv;
   return FiniteElement<T>(
-      element::family::HHJ, celltype, degree, {tdim, tdim},
-      impl::mdspan_t<T, 2>(wcoeffs.data(), wcoeffs.extents()), xview, Mview, 0,
-      maps::type::doubleContravariantPiola, space, discontinuous, -1, degree,
-      element::lagrange_variant::unset, element::dpc_variant::unset);
+      element::family::HHJ, celltype, polyset::type::standard, degree,
+      {tdim, tdim}, impl::mdspan_t<T, 2>(wcoeffs.data(), wcoeffs.extents()),
+      xview, Mview, 0, maps::type::doubleContravariantPiola, space,
+      discontinuous, -1, degree, element::lagrange_variant::unset,
+      element::dpc_variant::unset);
 }
 //-----------------------------------------------------------------------------
 template FiniteElement<float> element::create_hhj(cell::type, int, bool);
diff --git a/cpp/basix/e-lagrange.cpp b/cpp/basix/e-lagrange.cpp
index e74f7a95e..f74b746d0 100644
--- a/cpp/basix/e-lagrange.cpp
+++ b/cpp/basix/e-lagrange.cpp
@@ -27,7 +27,8 @@ impl::mdarray_t<T, 2> vtk_triangle_points(std::size_t degree)
   if (degree == 0)
     return basix::impl::mdarray_t<T, 2>({d, d}, 1, 2);
 
-  const std::size_t npoints = polyset::dim(cell::type::triangle, degree);
+  const std::size_t npoints
+      = polyset::dim(cell::type::triangle, polyset::type::standard, degree);
   impl::mdarray_t<T, 2> out(npoints, 2);
 
   out(0, 0) = d;
@@ -84,7 +85,8 @@ vtk_tetrahedron_points(std::size_t degree)
                                                                   2);
   }
 
-  const std::size_t npoints = polyset::dim(cell::type::tetrahedron, degree);
+  const std::size_t npoints
+      = polyset::dim(cell::type::tetrahedron, polyset::type::standard, degree);
   stdex::mdarray<T, stdex::extents<std::size_t, stdex::dynamic_extent, 3>> out(
       npoints, 3);
 
@@ -765,7 +767,8 @@ FiniteElement<T> create_d_lagrange(cell::type celltype, int degree,
   }
 
   const std::size_t tdim = cell::topological_dimension(celltype);
-  const std::size_t ndofs = polyset::dim(celltype, degree);
+  const std::size_t ndofs
+      = polyset::dim(celltype, polyset::type::standard, degree);
   const std::vector<std::vector<std::vector<int>>> topology
       = cell::topology(celltype);
 
@@ -796,7 +799,59 @@ FiniteElement<T> create_d_lagrange(cell::type celltype, int degree,
     _M(i, 0, i, 0) = 1.0;
 
   return FiniteElement(
-      element::family::P, celltype, degree, {},
+      element::family::P, celltype, polyset::type::standard, 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,
+      element::dpc_variant::unset);
+}
+//----------------------------------------------------------------------------
+template <std::floating_point T>
+FiniteElement<T>
+create_d_iso(cell::type celltype, int degree, element::lagrange_variant variant,
+             lattice::type lattice_type, lattice::simplex_method simplex_method)
+{
+  if (celltype == cell::type::prism or celltype == cell::type::pyramid)
+  {
+    throw std::runtime_error(
+        "This variant is not yet supported on prisms and pyramids.");
+  }
+
+  const std::size_t tdim = cell::topological_dimension(celltype);
+  const std::size_t ndofs
+      = polyset::dim(celltype, polyset::type::macroedge, degree);
+  const std::vector<std::vector<std::vector<int>>> topology
+      = cell::topology(celltype);
+
+  std::array<std::vector<impl::mdarray_t<T, 2>>, 4> x;
+  std::array<std::vector<impl::mdarray_t<T, 4>>, 4> M;
+  for (std::size_t i = 0; i < tdim; ++i)
+  {
+    std::size_t num_ent = cell::num_sub_entities(celltype, i);
+    x[i] = std::vector<impl::mdarray_t<T, 2>>(num_ent,
+                                              impl::mdarray_t<T, 2>(0, tdim));
+    M[i] = std::vector<impl::mdarray_t<T, 4>>(
+        num_ent, impl::mdarray_t<T, 4>(0, 1, 0, 1));
+  }
+
+  const int lattice_degree
+      = celltype == cell::type::triangle
+            ? 2 * degree + 3
+            : (celltype == cell::type::tetrahedron ? 2 * degree + 4
+                                                   : 2 * degree + 2);
+
+  // Create points in interior
+  const auto [pt, shape] = lattice::create<T>(
+      celltype, lattice_degree, lattice_type, false, simplex_method);
+  x[tdim].emplace_back(pt, shape);
+
+  const std::size_t num_dofs = shape[0];
+  auto& _M = M[tdim].emplace_back(num_dofs, 1, num_dofs, 1);
+  for (std::size_t i = 0; i < _M.extent(0); ++i)
+    _M(i, 0, i, 0) = 1.0;
+
+  return FiniteElement(
+      element::family::P, 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,
@@ -945,13 +1000,14 @@ FiniteElement<T> create_vtk_element(cell::type celltype, std::size_t degree,
   }
 
   const std::size_t tdim = cell::topological_dimension(celltype);
-  const std::size_t ndofs = polyset::dim(celltype, degree);
+  const std::size_t ndofs
+      = polyset::dim(celltype, polyset::type::standard, degree);
   if (discontinuous)
   {
     auto [_x, _xshape, _M, _Mshape] = element::make_discontinuous(
         impl::to_mdspan(x), impl::to_mdspan(M), tdim, 1);
     return FiniteElement(
-        element::family::P, celltype, degree, {},
+        element::family::P, celltype, polyset::type::standard, degree, {},
         impl::mdspan_t<const T, 2>(math::eye<T>(ndofs).data(), ndofs, ndofs),
         impl::to_mdspan(_x, _xshape), impl::to_mdspan(_M, _Mshape), 0,
         maps::type::identity, sobolev::space::L2, discontinuous, degree, degree,
@@ -960,7 +1016,7 @@ FiniteElement<T> create_vtk_element(cell::type celltype, std::size_t degree,
   else
   {
     return FiniteElement(
-        element::family::P, celltype, degree, {},
+        element::family::P, celltype, polyset::type::standard, 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::H1, discontinuous, degree, degree,
@@ -976,7 +1032,8 @@ FiniteElement<T> create_legendre(cell::type celltype, int degree,
     throw std::runtime_error("Legendre variant must be discontinuous");
 
   const std::size_t tdim = cell::topological_dimension(celltype);
-  const std::size_t ndofs = polyset::dim(celltype, degree);
+  const std::size_t ndofs
+      = polyset::dim(celltype, polyset::type::standard, degree);
   const std::vector<std::vector<std::vector<int>>> topology
       = cell::topology(celltype);
 
@@ -985,7 +1042,7 @@ FiniteElement<T> create_legendre(cell::type celltype, int degree,
 
   // Evaluate moment space at quadrature points
   const auto [_pts, wts] = quadrature::make_quadrature<T>(
-      quadrature::type::Default, celltype, degree * 2);
+      quadrature::type::Default, celltype, polyset::type::standard, degree * 2);
   assert(!wts.empty());
   impl::mdspan_t<const T, 2> pts(_pts.data(), wts.size(),
                                  _pts.size() / wts.size());
@@ -1011,7 +1068,7 @@ FiniteElement<T> create_legendre(cell::type celltype, int degree,
   sobolev::space space
       = discontinuous ? sobolev::space::L2 : sobolev::space::H1;
   return FiniteElement<T>(
-      element::family::P, celltype, degree, {},
+      element::family::P, celltype, polyset::type::standard, degree, {},
       impl::mdspan_t<T, 2>(math::eye<T>(ndofs).data(), ndofs, ndofs),
       impl::to_mdspan(x), impl::to_mdspan(M), 0, maps::type::identity, space,
       discontinuous, degree, degree, element::lagrange_variant::legendre,
@@ -1114,8 +1171,9 @@ FiniteElement<T> create_bernstein(cell::type celltype, int degree,
     }
     else
     {
-      const auto [_pts, wts] = quadrature::make_quadrature<T>(
-          quadrature::type::Default, ct[d], degree * 2);
+      const auto [_pts, wts]
+          = quadrature::make_quadrature<T>(quadrature::type::Default, ct[d],
+                                           polyset::type::standard, degree * 2);
       assert(!wts.empty());
       impl::mdspan_t<const T, 2> pts(_pts.data(), wts.size(),
                                      _pts.size() / wts.size());
@@ -1180,9 +1238,10 @@ FiniteElement<T> create_bernstein(cell::type celltype, int degree,
 
   sobolev::space space
       = discontinuous ? sobolev::space::L2 : sobolev::space::H1;
-  const std::size_t ndofs = polyset::dim(celltype, degree);
+  const std::size_t ndofs
+      = polyset::dim(celltype, polyset::type::standard, degree);
   return FiniteElement<T>(
-      element::family::P, celltype, degree, {},
+      element::family::P, celltype, polyset::type::standard, degree, {},
       impl::mdspan_t<T, 2>(math::eye<T>(ndofs).data(), ndofs, ndofs),
       impl::to_mdspan(x), impl::to_mdspan(M), 0, maps::type::identity, space,
       discontinuous, degree, degree, element::lagrange_variant::bernstein,
@@ -1208,7 +1267,7 @@ basix::element::create_lagrange(cell::type celltype, int degree,
     x[0].emplace_back(1, 0);
     M[0].emplace_back(std::vector<T>{1.0}, 1, 1, 1, 1);
     return FiniteElement<T>(
-        family::P, cell::type::point, 0, {},
+        family::P, cell::type::point, polyset::type::standard, 0, {},
         impl::mdspan_t<T, 2>(math::eye<T>(1).data(), 1, 1), impl::to_mdspan(x),
         impl::to_mdspan(M), 0, maps::type::identity, sobolev::space::H1,
         discontinuous, degree, degree, element::lagrange_variant::unset,
@@ -1257,7 +1316,8 @@ basix::element::create_lagrange(cell::type celltype, int degree,
   }
 
   const std::size_t tdim = cell::topological_dimension(celltype);
-  const std::size_t ndofs = polyset::dim(celltype, degree);
+  const std::size_t ndofs
+      = polyset::dim(celltype, polyset::type::standard, degree);
   const std::vector<std::vector<std::vector<int>>> topology
       = cell::topology(celltype);
 
@@ -1363,16 +1423,172 @@ basix::element::create_lagrange(cell::type celltype, int degree,
   auto tensor_factors
       = create_tensor_product_factors<T>(celltype, degree, variant);
   return FiniteElement<T>(
-      family::P, celltype, degree, {},
+      family::P, celltype, polyset::type::standard, 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, dof_ordering);
 }
 //-----------------------------------------------------------------------------
+template <std::floating_point T>
+FiniteElement<T> basix::element::create_iso(cell::type celltype, int degree,
+                                            lagrange_variant variant,
+                                            bool discontinuous)
+{
+  if (celltype != cell::type::interval)
+  {
+    throw std::runtime_error(
+        "Can currently only create iso elements on an interval");
+  }
+
+  if (variant == lagrange_variant::unset)
+  {
+    if (degree < 3)
+      variant = element::lagrange_variant::gll_warped;
+    else
+    {
+      throw std::runtime_error(
+          "Lagrange elements of degree > 2 need to be given a variant.");
+    }
+  }
+
+  auto [lattice_type, simplex_method, exterior]
+      = variant_to_lattice(celltype, variant);
+
+  if (!exterior)
+  {
+    // Points used to define this variant are all interior to the cell,
+    // so this variant requires that the element is discontinuous
+    if (!discontinuous)
+    {
+      throw std::runtime_error("This variant of Lagrange is only supported for "
+                               "discontinuous elements");
+    }
+    return create_d_iso<T>(celltype, degree, variant, lattice_type,
+                           simplex_method);
+  }
+
+  const std::size_t tdim = cell::topological_dimension(celltype);
+  const std::size_t ndofs
+      = polyset::dim(celltype, polyset::type::macroedge, degree);
+  const std::vector<std::vector<std::vector<int>>> topology
+      = cell::topology(celltype);
+
+  std::array<std::vector<impl::mdarray_t<T, 2>>, 4> x;
+  std::array<std::vector<impl::mdarray_t<T, 4>>, 4> M;
+  if (degree == 0)
+  {
+    if (!discontinuous)
+    {
+      throw std::runtime_error(
+          "Cannot create a continuous order 0 Lagrange basis function");
+    }
+
+    for (std::size_t i = 0; i < tdim; ++i)
+    {
+      std::size_t num_entities = cell::num_sub_entities(celltype, i);
+      x[i] = std::vector(num_entities, impl::mdarray_t<T, 2>(0, tdim));
+      M[i] = std::vector(num_entities, impl::mdarray_t<T, 4>(0, 1, 0, 1));
+    }
+
+    const auto [pt, shape]
+        = lattice::create<T>(celltype, 0, lattice_type, true, simplex_method);
+    x[tdim].emplace_back(pt, shape[0], shape[1]);
+    auto& _M = M[tdim].emplace_back(shape[0], 1, shape[0], 1);
+    std::fill(_M.data(), _M.data() + _M.size(), 0);
+    for (std::size_t i = 0; i < shape[0]; ++i)
+      _M(i, 0, i, 0) = 1;
+  }
+  else
+  {
+    // Create points at nodes, ordered by topology (vertices first)
+    for (std::size_t dim = 0; dim <= tdim; ++dim)
+    {
+      // Loop over entities of dimension 'dim'
+      for (std::size_t e = 0; e < topology[dim].size(); ++e)
+      {
+        const auto [entity_x, entity_x_shape]
+            = cell::sub_entity_geometry<T>(celltype, dim, e);
+        if (dim == 0)
+        {
+          x[dim].emplace_back(entity_x, entity_x_shape[0], entity_x_shape[1]);
+          auto& _M
+              = M[dim].emplace_back(entity_x_shape[0], 1, entity_x_shape[0], 1);
+          std::fill(_M.data(), _M.data() + _M.size(), 0);
+          for (std::size_t i = 0; i < entity_x_shape[0]; ++i)
+            _M(i, 0, i, 0) = 1;
+        }
+        else if (dim == tdim)
+        {
+          const auto [pt, shape] = lattice::create<T>(
+              celltype, 2 * degree, lattice_type, false, simplex_method);
+          x[dim].emplace_back(pt, shape[0], shape[1]);
+          auto& _M = M[dim].emplace_back(shape[0], 1, shape[0], 1);
+          std::fill(_M.data(), _M.data() + _M.size(), 0);
+          for (std::size_t i = 0; i < shape[0]; ++i)
+            _M(i, 0, i, 0) = 1;
+        }
+        else
+        {
+          cell::type ct = cell::sub_entity_type(celltype, dim, e);
+          const auto [pt, shape] = lattice::create<T>(
+              ct, 2 * degree, lattice_type, false, simplex_method);
+          impl::mdspan_t<const T, 2> lattice(pt.data(), shape);
+          std::span<const T> x0(entity_x.data(), entity_x_shape[1]);
+          impl::mdspan_t<const T, 2> entity_x_view(entity_x.data(),
+                                                   entity_x_shape);
+
+          auto& _x = x[dim].emplace_back(shape[0], entity_x_shape[1]);
+          for (std::size_t i = 0; i < shape[0]; ++i)
+            for (std::size_t j = 0; j < entity_x_shape[1]; ++j)
+              _x(i, j) = x0[j];
+
+          for (std::size_t j = 0; j < shape[0]; ++j)
+            for (std::size_t k = 0; k < shape[1]; ++k)
+              for (std::size_t q = 0; q < tdim; ++q)
+                _x(j, q) += (entity_x_view(k + 1, q) - x0[q]) * lattice(j, k);
+
+          auto& _M = M[dim].emplace_back(shape[0], 1, shape[0], 1);
+          std::fill(_M.data(), _M.data() + _M.size(), 0);
+          for (std::size_t i = 0; i < shape[0]; ++i)
+            _M(i, 0, i, 0) = 1;
+        }
+      }
+    }
+  }
+
+  std::array<std::vector<mdspan_t<const T, 2>>, 4> xview = impl::to_mdspan(x);
+  std::array<std::vector<mdspan_t<const T, 4>>, 4> Mview = impl::to_mdspan(M);
+  std::array<std::vector<std::vector<T>>, 4> xbuffer;
+  std::array<std::vector<std::vector<T>>, 4> Mbuffer;
+  if (discontinuous)
+  {
+    std::array<std::vector<std::array<std::size_t, 2>>, 4> xshape;
+    std::array<std::vector<std::array<std::size_t, 4>>, 4> Mshape;
+    std::tie(xbuffer, xshape, Mbuffer, Mshape)
+        = make_discontinuous(xview, Mview, tdim, 1);
+    xview = impl::to_mdspan(xbuffer, xshape);
+    Mview = impl::to_mdspan(Mbuffer, Mshape);
+  }
+
+  sobolev::space space
+      = discontinuous ? sobolev::space::L2 : sobolev::space::H1;
+  auto tensor_factors
+      = create_tensor_product_factors<T>(celltype, degree, variant);
+  return FiniteElement<T>(
+      family::P, 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);
+}
+//-----------------------------------------------------------------------------
 template FiniteElement<float> element::create_lagrange(cell::type, int,
                                                        lagrange_variant, bool,
                                                        std::vector<int>);
 template FiniteElement<double> element::create_lagrange(cell::type, int,
                                                         lagrange_variant, bool,
                                                         std::vector<int>);
+template FiniteElement<float> element::create_iso(cell::type, int,
+                                                  lagrange_variant, bool);
+template FiniteElement<double> element::create_iso(cell::type, int,
+                                                   lagrange_variant, bool);
 //-----------------------------------------------------------------------------
diff --git a/cpp/basix/e-lagrange.h b/cpp/basix/e-lagrange.h
index 42d7de203..28a3083b9 100644
--- a/cpp/basix/e-lagrange.h
+++ b/cpp/basix/e-lagrange.h
@@ -22,4 +22,14 @@ template <std::floating_point T>
 FiniteElement<T> create_lagrange(cell::type celltype, int degree,
                                  lagrange_variant variant, bool discontinuous,
                                  std::vector<int> dof_ordering = {});
+
+/// @brief Create an iso macro element on cell with given degree
+/// @param[in] celltype The element cell type
+/// @param[in] degree The degree of the element
+/// @param[in] variant The variant of the element to be created
+/// @param[in] discontinuous True if the is discontinuous
+/// @return A finite element
+template <std::floating_point T>
+FiniteElement<T> create_iso(cell::type celltype, int degree,
+                            lagrange_variant variant, bool discontinuous);
 } // namespace basix::element
diff --git a/cpp/basix/e-nce-rtc.cpp b/cpp/basix/e-nce-rtc.cpp
index 47b4d59e5..3c23eee1a 100644
--- a/cpp/basix/e-nce-rtc.cpp
+++ b/cpp/basix/e-nce-rtc.cpp
@@ -37,25 +37,29 @@ FiniteElement<T> basix::element::create_rtc(cell::type celltype, int degree,
 
   // Evaluate the expansion polynomials at the quadrature points
   const auto [_pts, qwts] = quadrature::make_quadrature<T>(
-      quadrature::type::Default, celltype, 2 * degree);
+      quadrature::type::Default, celltype, polyset::type::standard, 2 * degree);
   impl::mdspan_t<const T, 2> pts(_pts.data(), qwts.size(),
                                  _pts.size() / qwts.size());
-  const auto [_phi, shape] = polyset::tabulate(celltype, degree, 0, pts);
+  const auto [_phi, shape]
+      = polyset::tabulate(celltype, polyset::type::standard, degree, 0, pts);
   impl::mdspan_t<const T, 3> phi(_phi.data(), shape);
 
   // The number of order (degree) polynomials
   const std::size_t psize = phi.extent(1);
 
   const int facet_count = tdim == 2 ? 4 : 6;
-  const int facet_dofs = polyset::dim(facettype, degree - 1);
+  const int facet_dofs
+      = polyset::dim(facettype, polyset::type::standard, degree - 1);
   const int internal_dofs = tdim == 2 ? 2 * degree * (degree - 1)
                                       : 3 * degree * degree * (degree - 1);
   const std::size_t ndofs = facet_count * facet_dofs + internal_dofs;
 
   // Create coefficients for order (degree-1) vector polynomials
   impl::mdarray_t<T, 2> wcoeffs(ndofs, psize * tdim);
-  const int nv_interval = polyset::dim(cell::type::interval, degree);
-  const int ns_interval = polyset::dim(cell::type::interval, degree - 1);
+  const int nv_interval
+      = polyset::dim(cell::type::interval, polyset::type::standard, degree);
+  const int ns_interval
+      = polyset::dim(cell::type::interval, polyset::type::standard, degree - 1);
   int dof = 0;
   if (tdim == 2)
   {
@@ -129,7 +133,7 @@ FiniteElement<T> basix::element::create_rtc(cell::type celltype, int degree,
     FiniteElement<T> moment_space
         = element::create_lagrange<T>(facettype, degree - 1, lvariant, true);
     auto [_x, xshape, _M, Mshape] = moments::make_normal_integral_moments<T>(
-        moment_space, celltype, tdim, 2 * degree - 1);
+        moment_space, celltype, polyset::type::standard, tdim, 2 * degree - 1);
     assert(_x.size() == _M.size());
     for (std::size_t i = 0; i < _x.size(); ++i)
     {
@@ -144,7 +148,7 @@ FiniteElement<T> basix::element::create_rtc(cell::type celltype, int degree,
   {
     auto [_x, xshape, _M, Mshape] = moments::make_dot_integral_moments<T>(
         element::create_nce<T>(celltype, degree - 1, lvariant, true), celltype,
-        tdim, 2 * degree - 1);
+        polyset::type::standard, tdim, 2 * degree - 1);
     assert(_x.size() == _M.size());
     for (std::size_t i = 0; i < _x.size(); ++i)
     {
@@ -176,7 +180,7 @@ FiniteElement<T> basix::element::create_rtc(cell::type celltype, int degree,
   sobolev::space space
       = discontinuous ? sobolev::space::L2 : sobolev::space::HCurl;
   return FiniteElement<T>(
-      element::family::RT, celltype, degree, {tdim},
+      element::family::RT, celltype, polyset::type::standard, degree, {tdim},
       impl::mdspan_t<T, 2>(wcoeffs.data(), wcoeffs.extents()), xview, Mview, 0,
       maps::type::contravariantPiola, space, discontinuous, degree - 1, degree,
       lvariant, element::dpc_variant::unset);
@@ -200,17 +204,19 @@ FiniteElement<T> basix::element::create_nce(cell::type celltype, int degree,
 
   // Evaluate the expansion polynomials at the quadrature points
   const auto [_pts, wts] = quadrature::make_quadrature<T>(
-      quadrature::type::Default, celltype, 2 * degree);
+      quadrature::type::Default, celltype, polyset::type::standard, 2 * degree);
   impl::mdspan_t<const T, 2> pts(_pts.data(), wts.size(),
                                  _pts.size() / wts.size());
-  const auto [_phi, shape] = polyset::tabulate(celltype, degree, 0, pts);
+  const auto [_phi, shape]
+      = polyset::tabulate(celltype, polyset::type::standard, degree, 0, pts);
   impl::mdspan_t<const T, 3> phi(_phi.data(), shape);
 
   // The number of order (degree) polynomials
   const int psize = phi.extent(1);
 
   const int edge_count = tdim == 2 ? 4 : 12;
-  const int edge_dofs = polyset::dim(cell::type::interval, degree - 1);
+  const int edge_dofs
+      = polyset::dim(cell::type::interval, polyset::type::standard, degree - 1);
   const int face_count = tdim == 2 ? 1 : 6;
   const int face_dofs = 2 * degree * (degree - 1);
   const int volume_count = tdim == 2 ? 0 : 1;
@@ -221,8 +227,10 @@ FiniteElement<T> basix::element::create_nce(cell::type celltype, int degree,
   // Create coefficients for order (degree-1) vector polynomials
   impl::mdarray_t<T, 2> wcoeffs(ndofs, psize * tdim);
 
-  const int nv_interval = polyset::dim(cell::type::interval, degree);
-  const int ns_interval = polyset::dim(cell::type::interval, degree - 1);
+  const int nv_interval
+      = polyset::dim(cell::type::interval, polyset::type::standard, degree);
+  const int ns_interval
+      = polyset::dim(cell::type::interval, polyset::type::standard, degree - 1);
   int dof = 0;
   if (tdim == 2)
   {
@@ -332,7 +340,8 @@ FiniteElement<T> basix::element::create_nce(cell::type celltype, int degree,
     FiniteElement<T> edge_moment_space = element::create_lagrange<T>(
         cell::type::interval, degree - 1, lvariant, true);
     auto [_x, xshape, _M, Mshape] = moments::make_tangent_integral_moments<T>(
-        edge_moment_space, celltype, tdim, 2 * degree - 1);
+        edge_moment_space, celltype, polyset::type::standard, tdim,
+        2 * degree - 1);
     assert(_x.size() == _M.size());
     for (std::size_t i = 0; i < _x.size(); ++i)
     {
@@ -348,7 +357,7 @@ FiniteElement<T> basix::element::create_nce(cell::type celltype, int degree,
     FiniteElement<T> moment_space = element::create_rtc<T>(
         cell::type::quadrilateral, degree - 1, lvariant, true);
     auto [_x, xshape, _M, Mshape] = moments::make_dot_integral_moments<T>(
-        moment_space, celltype, tdim, 2 * degree - 1);
+        moment_space, celltype, polyset::type::standard, tdim, 2 * degree - 1);
     assert(_x.size() == _M.size());
     for (std::size_t i = 0; i < _x.size(); ++i)
     {
@@ -369,7 +378,8 @@ FiniteElement<T> basix::element::create_nce(cell::type celltype, int degree,
       FiniteElement<T> moment_space = element::create_rtc<T>(
           cell::type::hexahedron, degree - 1, lvariant, true);
       auto [_x, xshape, _M, Mshape] = moments::make_dot_integral_moments<T>(
-          moment_space, celltype, tdim, 2 * degree - 1);
+          moment_space, celltype, polyset::type::standard, tdim,
+          2 * degree - 1);
       assert(_x.size() == _M.size());
       for (std::size_t i = 0; i < _x.size(); ++i)
       {
@@ -402,7 +412,7 @@ FiniteElement<T> basix::element::create_nce(cell::type celltype, int degree,
   sobolev::space space
       = discontinuous ? sobolev::space::L2 : sobolev::space::HCurl;
   return FiniteElement<T>(
-      element::family::N1E, celltype, degree, {tdim},
+      element::family::N1E, celltype, polyset::type::standard, degree, {tdim},
       impl::mdspan_t<T, 2>(wcoeffs.data(), wcoeffs.extents()), xview, Mview, 0,
       maps::type::covariantPiola, space, discontinuous, degree - 1, degree,
       lvariant, element::dpc_variant::unset);
diff --git a/cpp/basix/e-nedelec.cpp b/cpp/basix/e-nedelec.cpp
index 89199c4f2..064e4482d 100644
--- a/cpp/basix/e-nedelec.cpp
+++ b/cpp/basix/e-nedelec.cpp
@@ -36,11 +36,12 @@ impl::mdarray_t<T, 2> create_nedelec_2d_space(int degree)
 
   // Tabulate polynomial set at quadrature points
   const auto [_pts, wts] = quadrature::make_quadrature<T>(
-      quadrature::type::Default, cell::type::triangle, 2 * degree);
+      quadrature::type::Default, cell::type::triangle, polyset::type::standard,
+      2 * degree);
   impl::mdspan_t<const T, 2> pts(_pts.data(), wts.size(),
                                  _pts.size() / wts.size());
-  const auto [_phi, shape]
-      = polyset::tabulate(cell::type::triangle, degree, 0, pts);
+  const auto [_phi, shape] = polyset::tabulate(
+      cell::type::triangle, polyset::type::standard, degree, 0, pts);
   impl::mdspan_t<const T, 3> phi(_phi.data(), shape);
 
   const std::size_t psize = phi.extent(1);
@@ -97,11 +98,12 @@ impl::mdarray_t<T, 2> create_nedelec_3d_space(int degree)
 
   // Tabulate polynomial basis at quadrature points
   const auto [_pts, wts] = quadrature::make_quadrature<T>(
-      quadrature::type::Default, cell::type::tetrahedron, 2 * degree);
+      quadrature::type::Default, cell::type::tetrahedron,
+      polyset::type::standard, 2 * degree);
   impl::mdspan_t<const T, 2> pts(_pts.data(), wts.size(),
                                  _pts.size() / wts.size());
-  const auto [_phi, shape]
-      = polyset::tabulate(cell::type::tetrahedron, degree, 0, pts);
+  const auto [_phi, shape] = polyset::tabulate(
+      cell::type::tetrahedron, polyset::type::standard, degree, 0, pts);
   impl::mdspan_t<const T, 3> phi(_phi.data(), shape);
   const std::size_t psize = phi.extent(1);
 
@@ -209,7 +211,7 @@ FiniteElement<T> element::create_nedelec(cell::type celltype, int degree,
     FiniteElement edge_space = element::create_lagrange<T>(
         cell::type::interval, degree - 1, lvariant, true);
     auto [_x, xshape, _M, Mshape] = moments::make_tangent_integral_moments<T>(
-        edge_space, celltype, tdim, 2 * degree - 1);
+        edge_space, celltype, polyset::type::standard, tdim, 2 * degree - 1);
     assert(_x.size() == _M.size());
     for (std::size_t i = 0; i < _x.size(); ++i)
     {
@@ -224,7 +226,7 @@ FiniteElement<T> element::create_nedelec(cell::type celltype, int degree,
     FiniteElement face_space = element::create_lagrange<T>(
         cell::type::triangle, degree - 2, lvariant, true);
     auto [_x, xshape, _M, Mshape] = moments::make_integral_moments<T>(
-        face_space, celltype, tdim, 2 * degree - 2);
+        face_space, celltype, polyset::type::standard, tdim, 2 * degree - 2);
     assert(_x.size() == _M.size());
     for (std::size_t i = 0; i < _x.size(); ++i)
     {
@@ -247,7 +249,7 @@ FiniteElement<T> element::create_nedelec(cell::type celltype, int degree,
       auto [_x, xshape, _M, Mshape] = moments::make_integral_moments<T>(
           element::create_lagrange<T>(cell::type::tetrahedron, degree - 3,
                                       lvariant, true),
-          cell::type::tetrahedron, 3, 2 * degree - 3);
+          cell::type::tetrahedron, polyset::type::standard, 3, 2 * degree - 3);
       assert(_x.size() == _M.size());
       for (std::size_t i = 0; i < _x.size(); ++i)
       {
@@ -280,11 +282,11 @@ FiniteElement<T> element::create_nedelec(cell::type celltype, int degree,
 
   sobolev::space space
       = discontinuous ? sobolev::space::L2 : sobolev::space::HCurl;
-  return FiniteElement<T>(element::family::N1E, celltype, degree, {tdim},
-                          impl::mdspan_t<T, 2>(wcoeffs.data(), wshape), xview,
-                          Mview, 0, maps::type::covariantPiola, space,
-                          discontinuous, degree - 1, degree, lvariant,
-                          element::dpc_variant::unset);
+  return FiniteElement<T>(
+      element::family::N1E, celltype, polyset::type::standard, degree, {tdim},
+      impl::mdspan_t<T, 2>(wcoeffs.data(), wshape), xview, Mview, 0,
+      maps::type::covariantPiola, space, discontinuous, degree - 1, degree,
+      lvariant, element::dpc_variant::unset);
 }
 //-----------------------------------------------------------------------------
 template <std::floating_point T>
@@ -314,7 +316,7 @@ FiniteElement<T> element::create_nedelec2(cell::type celltype, int degree,
     FiniteElement edge_space = element::create_lagrange<T>(
         cell::type::interval, degree, lvariant, true);
     auto [_x, xshape, _M, Mshape] = moments::make_tangent_integral_moments<T>(
-        edge_space, celltype, tdim, 2 * degree);
+        edge_space, celltype, polyset::type::standard, tdim, 2 * degree);
     assert(_x.size() == _M.size());
     for (std::size_t i = 0; i < _x.size(); ++i)
     {
@@ -329,7 +331,7 @@ FiniteElement<T> element::create_nedelec2(cell::type celltype, int degree,
     FiniteElement face_space = element::create_rt<T>(
         cell::type::triangle, degree - 1, lvariant, true);
     auto [_x, xshape, _M, Mshape] = moments::make_dot_integral_moments<T>(
-        face_space, celltype, tdim, 2 * degree - 1);
+        face_space, celltype, polyset::type::standard, tdim, 2 * degree - 1);
     assert(_x.size() == _M.size());
     for (std::size_t i = 0; i < _x.size(); ++i)
     {
@@ -351,7 +353,7 @@ FiniteElement<T> element::create_nedelec2(cell::type celltype, int degree,
       auto [_x, xshape, _M, Mshape] = moments::make_dot_integral_moments<T>(
           element::create_rt<T>(cell::type::tetrahedron, degree - 2, lvariant,
                                 true),
-          celltype, tdim, 2 * degree - 2);
+          celltype, polyset::type::standard, tdim, 2 * degree - 2);
       assert(_x.size() == _M.size());
       for (std::size_t i = 0; i < _x.size(); ++i)
       {
@@ -381,9 +383,10 @@ FiniteElement<T> element::create_nedelec2(cell::type celltype, int degree,
     Mview = impl::to_mdspan(Mbuffer, Mshape);
   }
 
-  const std::size_t psize = polyset::dim(celltype, degree);
+  const std::size_t psize
+      = polyset::dim(celltype, polyset::type::standard, degree);
   return FiniteElement<T>(
-      element::family::N2E, celltype, degree, {tdim},
+      element::family::N2E, celltype, polyset::type::standard, degree, {tdim},
       impl::mdspan_t<T, 2>(math::eye<T>(tdim * psize).data(), tdim * psize,
                            tdim * psize),
       xview, Mview, 0, maps::type::covariantPiola, sobolev::space::HCurl,
diff --git a/cpp/basix/e-raviart-thomas.cpp b/cpp/basix/e-raviart-thomas.cpp
index d96e88eea..b46f51b3a 100644
--- a/cpp/basix/e-raviart-thomas.cpp
+++ b/cpp/basix/e-raviart-thomas.cpp
@@ -31,21 +31,25 @@ FiniteElement<T> basix::element::create_rt(cell::type celltype, int degree,
       = (tdim == 2) ? cell::type::interval : cell::type::triangle;
 
   // The number of order (degree-1) scalar polynomials
-  const std::size_t nv = polyset::dim(celltype, degree - 1);
+  const std::size_t nv
+      = polyset::dim(celltype, polyset::type::standard, degree - 1);
 
   // The number of order (degree-2) scalar polynomials
-  const std::size_t ns0 = polyset::dim(celltype, degree - 2);
+  const std::size_t ns0
+      = polyset::dim(celltype, polyset::type::standard, degree - 2);
 
   // The number of additional polynomials in the polynomial basis for
   // Raviart-Thomas
-  const std::size_t ns = polyset::dim(facettype, degree - 1);
+  const std::size_t ns
+      = polyset::dim(facettype, polyset::type::standard, degree - 1);
 
   // Evaluate the expansion polynomials at the quadrature points
   const auto [_pts, wts] = quadrature::make_quadrature<T>(
-      quadrature::type::Default, celltype, 2 * degree);
+      quadrature::type::Default, celltype, polyset::type::standard, 2 * degree);
   impl::mdspan_t<const T, 2> pts(_pts.data(), wts.size(),
                                  _pts.size() / wts.size());
-  const auto [_phi, shape] = polyset::tabulate(celltype, degree, 0, pts);
+  const auto [_phi, shape]
+      = polyset::tabulate(celltype, polyset::type::standard, degree, 0, pts);
   impl::mdspan_t<const T, 3> phi(_phi.data(), shape);
 
   // The number of order (degree) polynomials
@@ -89,7 +93,8 @@ FiniteElement<T> basix::element::create_rt(cell::type celltype, int degree,
     const FiniteElement facet_moment_space
         = element::create_lagrange<T>(facettype, degree - 1, lvariant, true);
     auto [_x, xshape, _M, Mshape] = moments::make_normal_integral_moments<T>(
-        facet_moment_space, celltype, tdim, 2 * degree - 1);
+        facet_moment_space, celltype, polyset::type::standard, tdim,
+        2 * degree - 1);
     assert(_x.size() == _M.size());
     for (std::size_t i = 0; i < _x.size(); ++i)
     {
@@ -105,7 +110,7 @@ FiniteElement<T> basix::element::create_rt(cell::type celltype, int degree,
     // Interior integral moment
     auto [_x, xshape, _M, Mshape] = moments::make_integral_moments<T>(
         element::create_lagrange<T>(celltype, degree - 2, lvariant, true),
-        celltype, tdim, 2 * degree - 2);
+        celltype, polyset::type::standard, tdim, 2 * degree - 2);
     assert(_x.size() == _M.size());
     for (std::size_t i = 0; i < _x.size(); ++i)
     {
@@ -136,11 +141,11 @@ FiniteElement<T> basix::element::create_rt(cell::type celltype, int degree,
 
   sobolev::space space
       = discontinuous ? sobolev::space::L2 : sobolev::space::HDiv;
-  return FiniteElement<T>(element::family::RT, celltype, degree, {tdim},
-                          impl::mdspan_t<T, 2>(B.data(), B.extents()), xview,
-                          Mview, 0, maps::type::contravariantPiola, space,
-                          discontinuous, degree - 1, degree, lvariant,
-                          element::dpc_variant::unset);
+  return FiniteElement<T>(
+      element::family::RT, celltype, polyset::type::standard, degree, {tdim},
+      impl::mdspan_t<T, 2>(B.data(), B.extents()), xview, Mview, 0,
+      maps::type::contravariantPiola, space, discontinuous, degree - 1, degree,
+      lvariant, element::dpc_variant::unset);
 }
 //-----------------------------------------------------------------------------
 template FiniteElement<float>
diff --git a/cpp/basix/e-regge.cpp b/cpp/basix/e-regge.cpp
index 87ffe41ad..a6dd960f0 100644
--- a/cpp/basix/e-regge.cpp
+++ b/cpp/basix/e-regge.cpp
@@ -25,7 +25,8 @@ FiniteElement<T> element::create_regge(cell::type celltype, int degree,
   const std::size_t tdim = cell::topological_dimension(celltype);
 
   const int nc = tdim * (tdim + 1) / 2;
-  const int basis_size = polyset::dim(celltype, degree);
+  const int basis_size
+      = polyset::dim(celltype, polyset::type::standard, degree);
   const std::size_t ndofs = basis_size * nc;
   const std::size_t psize = basis_size * tdim * tdim;
 
@@ -84,9 +85,10 @@ FiniteElement<T> element::create_regge(cell::type celltype, int degree,
         // Tabulate points in lattice
         cell::type ct = cell::sub_entity_type(celltype, d, e);
 
-        const std::size_t ndofs = polyset::dim(ct, degree + 1 - d);
-        const auto [_pts, wts]
-            = quadrature::make_quadrature<T>(ct, degree + (degree + 1 - d));
+        const std::size_t ndofs
+            = polyset::dim(ct, polyset::type::standard, degree + 1 - d);
+        const auto [_pts, wts] = quadrature::make_quadrature<T>(
+            quadrature::type::Default, ct, polyset::type::standard, degree + (degree + 1 - d));
         impl::mdspan_t<const T, 2> pts(_pts.data(), wts.size(),
                                        _pts.size() / wts.size());
 
@@ -178,10 +180,11 @@ FiniteElement<T> element::create_regge(cell::type celltype, int degree,
   sobolev::space space
       = discontinuous ? sobolev::space::L2 : sobolev::space::HEin;
   return FiniteElement<T>(
-      element::family::Regge, celltype, degree, {tdim, tdim},
-      impl::mdspan_t<T, 2>(wcoeffs.data(), wcoeffs.extents()), xview, Mview, 0,
-      maps::type::doubleCovariantPiola, space, discontinuous, -1, degree,
-      element::lagrange_variant::unset, element::dpc_variant::unset);
+      element::family::Regge, celltype, polyset::type::standard, degree,
+      {tdim, tdim}, impl::mdspan_t<T, 2>(wcoeffs.data(), wcoeffs.extents()),
+      xview, Mview, 0, maps::type::doubleCovariantPiola, space, discontinuous,
+      -1, degree, element::lagrange_variant::unset,
+      element::dpc_variant::unset);
 }
 //-----------------------------------------------------------------------------
 template FiniteElement<float> element::create_regge(cell::type, int, bool);
diff --git a/cpp/basix/e-serendipity.cpp b/cpp/basix/e-serendipity.cpp
index a5b1647e3..bdde4f5a5 100644
--- a/cpp/basix/e-serendipity.cpp
+++ b/cpp/basix/e-serendipity.cpp
@@ -28,11 +28,12 @@ impl::mdarray_t<T, 2> make_serendipity_space_2d(int degree)
 
   // Evaluate the expansion polynomials at the quadrature points
   const auto [_pts, wts] = quadrature::make_quadrature<T>(
-      quadrature::type::Default, cell::type::quadrilateral, 2 * degree);
+      quadrature::type::Default, cell::type::quadrilateral,
+      polyset::type::standard, 2 * degree);
   impl::mdspan_t<const T, 2> pts(_pts.data(), wts.size(),
                                  _pts.size() / wts.size());
-  const auto [_Pq, shape]
-      = polyset::tabulate(cell::type::quadrilateral, degree, 0, pts);
+  const auto [_Pq, shape] = polyset::tabulate(
+      cell::type::quadrilateral, polyset::type::standard, degree, 0, pts);
   impl::mdspan_t<const T, 3> Pq(_Pq.data(), shape);
 
   const std::size_t psize = Pq.extent(1);
@@ -127,12 +128,13 @@ impl::mdarray_t<T, 2> make_serendipity_space_3d(int degree)
 
   // Evaluate the expansion polynomials at the quadrature points
   const auto [_pts, wts] = quadrature::make_quadrature<T>(
-      quadrature::type::Default, cell::type::hexahedron, 2 * degree);
+      quadrature::type::Default, cell::type::hexahedron,
+      polyset::type::standard, 2 * degree);
   impl::mdspan_t<const T, 2> pts(_pts.data(), wts.size(),
                                  _pts.size() / wts.size());
 
-  const auto [_Ph, shape]
-      = polyset::tabulate(cell::type::hexahedron, degree, 0, pts);
+  const auto [_Ph, shape] = polyset::tabulate(
+      cell::type::hexahedron, polyset::type::standard, degree, 0, pts);
   impl::mdspan_t<const T, 3> Ph(_Ph.data(), shape);
   const std::size_t psize = Ph.extent(1);
 
@@ -182,15 +184,17 @@ impl::mdarray_t<T, 2> make_serendipity_div_space_2d(int degree)
 
   // Evaluate the expansion polynomials at the quadrature points
   auto [_pts, wts] = quadrature::make_quadrature<T>(
-      quadrature::type::Default, cell::type::quadrilateral, 2 * degree + 2);
+      quadrature::type::Default, cell::type::quadrilateral,
+      polyset::type::standard, 2 * degree + 2);
   impl::mdspan_t<const T, 2> pts(_pts.data(), wts.size(),
                                  _pts.size() / wts.size());
-  const auto [_Pq, shape]
-      = polyset::tabulate(cell::type::quadrilateral, degree + 1, 0, pts);
+  const auto [_Pq, shape] = polyset::tabulate(
+      cell::type::quadrilateral, polyset::type::standard, degree + 1, 0, pts);
   impl::mdspan_t<const T, 3> Pq(_Pq.data(), shape);
 
   const std::size_t psize = Pq.extent(1);
-  const std::size_t nv = polyset::dim(cell::type::triangle, degree);
+  const std::size_t nv
+      = polyset::dim(cell::type::triangle, polyset::type::standard, degree);
 
   // Create coefficients for order (degree) vector polynomials
   impl::mdarray_t<T, 2> wcoeffs(ndofs, psize * 2);
@@ -251,16 +255,18 @@ impl::mdarray_t<T, 2> make_serendipity_div_space_3d(int degree)
 
   // Evaluate the expansion polynomials at the quadrature points
   const auto [_pts, wts] = quadrature::make_quadrature<T>(
-      quadrature::type::Default, cell::type::hexahedron, 2 * degree + 2);
+      quadrature::type::Default, cell::type::hexahedron,
+      polyset::type::standard, 2 * degree + 2);
   impl::mdspan_t<const T, 2> pts(_pts.data(), wts.size(),
                                  _pts.size() / wts.size());
 
-  const auto [_Pq, shape]
-      = polyset::tabulate(cell::type::hexahedron, degree + 1, 0, pts);
+  const auto [_Pq, shape] = polyset::tabulate(
+      cell::type::hexahedron, polyset::type::standard, degree + 1, 0, pts);
   impl::mdspan_t<const T, 3> Pq(_Pq.data(), shape);
 
   const std::size_t psize = Pq.extent(1);
-  const std::size_t nv = polyset::dim(cell::type::tetrahedron, degree);
+  const std::size_t nv
+      = polyset::dim(cell::type::tetrahedron, polyset::type::standard, degree);
 
   // Create coefficients for order (degree) vector polynomials
   impl::mdarray_t<T, 2> wcoeffs(ndofs, psize * 3);
@@ -398,15 +404,17 @@ impl::mdarray_t<T, 2> make_serendipity_curl_space_2d(int degree)
 
   // Evaluate the expansion polynomials at the quadrature points
   const auto [_pts, wts] = quadrature::make_quadrature<T>(
-      quadrature::type::Default, cell::type::quadrilateral, 2 * degree + 2);
+      quadrature::type::Default, cell::type::quadrilateral,
+      polyset::type::standard, 2 * degree + 2);
   impl::mdspan_t<const T, 2> pts(_pts.data(), wts.size(),
                                  _pts.size() / wts.size());
-  const auto [_Pq, shape]
-      = polyset::tabulate(cell::type::quadrilateral, degree + 1, 0, pts);
+  const auto [_Pq, shape] = polyset::tabulate(
+      cell::type::quadrilateral, polyset::type::standard, degree + 1, 0, pts);
   impl::mdspan_t<const T, 3> Pq(_Pq.data(), shape);
 
   const std::size_t psize = Pq.extent(1);
-  const std::size_t nv = polyset::dim(cell::type::triangle, degree);
+  const std::size_t nv
+      = polyset::dim(cell::type::triangle, polyset::type::standard, degree);
 
   // Create coefficients for order (degree) vector polynomials
   impl::mdarray_t<T, 2> wcoeffs(ndofs, psize * 2);
@@ -472,15 +480,17 @@ impl::mdarray_t<T, 2> make_serendipity_curl_space_3d(int degree)
 
   // Evaluate the expansion polynomials at the quadrature points
   const auto [_pts, wts] = quadrature::make_quadrature<T>(
-      quadrature::type::Default, cell::type::hexahedron, 2 * degree + 2);
+      quadrature::type::Default, cell::type::hexahedron,
+      polyset::type::standard, 2 * degree + 2);
   impl::mdspan_t<const T, 2> pts(_pts.data(), wts.size(),
                                  _pts.size() / wts.size());
-  const auto [_Pq, shape]
-      = polyset::tabulate(cell::type::hexahedron, degree + 1, 0, pts);
+  const auto [_Pq, shape] = polyset::tabulate(
+      cell::type::hexahedron, polyset::type::standard, degree + 1, 0, pts);
   impl::mdspan_t<const T, 3> Pq(_Pq.data(), shape);
 
   const std::size_t psize = Pq.extent(1);
-  const std::size_t nv = polyset::dim(cell::type::tetrahedron, degree);
+  const std::size_t nv
+      = polyset::dim(cell::type::tetrahedron, polyset::type::standard, degree);
 
   // Create coefficients for order (degree) vector polynomials
   impl::mdarray_t<T, 2> wcoeffs(ndofs, psize * 3);
@@ -673,11 +683,13 @@ FiniteElement<T> create_legendre_dpc(cell::type celltype, int degree,
   }
 
   const std::size_t tdim = cell::topological_dimension(celltype);
-  const std::size_t psize = polyset::dim(celltype, degree);
-  const std::size_t ndofs = polyset::dim(simplex_type, degree);
+  const std::size_t psize
+      = polyset::dim(celltype, polyset::type::standard, degree);
+  const std::size_t ndofs
+      = polyset::dim(simplex_type, polyset::type::standard, degree);
 
   const auto [_pts, wts] = quadrature::make_quadrature<T>(
-      quadrature::type::Default, celltype, degree * 2);
+      quadrature::type::Default, celltype, polyset::type::standard, degree * 2);
   impl::mdspan_t<const T, 2> pts(_pts.data(), wts.size(),
                                  _pts.size() / wts.size());
 
@@ -738,7 +750,7 @@ FiniteElement<T> create_legendre_dpc(cell::type celltype, int degree,
   }
 
   return FiniteElement<T>(
-      element::family::DPC, celltype, degree, {},
+      element::family::DPC, celltype, polyset::type::standard, degree, {},
       impl::mdspan_t<T, 2>(wcoeffs.data(), wcoeffs.extents()),
       impl::to_mdspan(x), impl::to_mdspan(M), 0, maps::type::identity,
       sobolev::space::L2, discontinuous, degree, degree,
@@ -1010,7 +1022,7 @@ FiniteElement<T> element::create_serendipity(cell::type celltype, int degree,
     FiniteElement moment_space = element::create_lagrange<T>(
         cell::type::interval, degree - 2, lvariant, true);
     auto [_x, xshape, _M, Mshape] = moments::make_integral_moments<T>(
-        moment_space, celltype, 1, 2 * degree - 2);
+        moment_space, celltype, polyset::type::standard, 1, 2 * degree - 2);
     assert(_x.size() == _M.size());
     for (std::size_t i = 0; i < _x.size(); ++i)
     {
@@ -1032,7 +1044,7 @@ FiniteElement<T> element::create_serendipity(cell::type celltype, int degree,
       FiniteElement moment_space = element::create_dpc<T>(
           cell::type::quadrilateral, degree - 4, dvariant, true);
       auto [_x, xshape, _M, Mshape] = moments::make_integral_moments<T>(
-          moment_space, celltype, 1, 2 * degree - 4);
+          moment_space, celltype, polyset::type::standard, 1, 2 * degree - 4);
       assert(_x.size() == _M.size());
       for (std::size_t i = 0; i < _x.size(); ++i)
       {
@@ -1055,7 +1067,7 @@ FiniteElement<T> element::create_serendipity(cell::type celltype, int degree,
       auto [_x, xshape, _M, Mshape] = moments::make_integral_moments<T>(
           element::create_dpc<T>(cell::type::hexahedron, degree - 6, dvariant,
                                  true),
-          celltype, 1, 2 * degree - 6);
+          celltype, polyset::type::standard, 1, 2 * degree - 6);
       assert(_x.size() == _M.size());
       for (std::size_t i = 0; i < _x.size(); ++i)
       {
@@ -1112,12 +1124,12 @@ FiniteElement<T> element::create_serendipity(cell::type celltype, int degree,
 
   sobolev::space space
       = discontinuous ? sobolev::space::L2 : sobolev::space::H1;
-  return FiniteElement<T>(element::family::serendipity, celltype, degree, {},
-                          impl::mdspan_t<const T, 2>(wbuffer.data(), wshape),
-                          xview, Mview, 0, maps::type::identity, space,
-                          discontinuous,
-                          degree < static_cast<int>(tdim) ? 1 : degree / tdim,
-                          degree, lvariant, dvariant);
+  return FiniteElement<T>(
+      element::family::serendipity, celltype, polyset::type::standard, degree,
+      {}, impl::mdspan_t<const T, 2>(wbuffer.data(), wshape), xview, Mview, 0,
+      maps::type::identity, space, discontinuous,
+      degree < static_cast<int>(tdim) ? 1 : degree / tdim, degree, lvariant,
+      dvariant);
 }
 //----------------------------------------------------------------------------
 template <std::floating_point T>
@@ -1157,8 +1169,10 @@ FiniteElement<T> element::create_dpc(cell::type celltype, int degree,
   if (variant == element::dpc_variant::legendre)
     return create_legendre_dpc<T>(celltype, degree, discontinuous);
 
-  const std::size_t ndofs = polyset::dim(simplex_type, degree);
-  const std::size_t psize = polyset::dim(celltype, degree);
+  const std::size_t ndofs
+      = polyset::dim(simplex_type, polyset::type::standard, degree);
+  const std::size_t psize
+      = polyset::dim(celltype, polyset::type::standard, degree);
   impl::mdarray_t<T, 2> wcoeffs(ndofs, psize);
   if (celltype == cell::type::quadrilateral)
   {
@@ -1202,7 +1216,7 @@ FiniteElement<T> element::create_dpc(cell::type celltype, int degree,
   impl::mdarray_t<T, 2> pt = make_dpc_points<T>(celltype, degree, variant);
   x[tdim].push_back(pt);
   return FiniteElement<T>(
-      element::family::DPC, celltype, degree, {},
+      element::family::DPC, celltype, polyset::type::standard, degree, {},
       impl::mdspan_t<T, 2>(wcoeffs.data(), wcoeffs.extents()),
       impl::to_mdspan(x), impl::to_mdspan(M), 0, maps::type::identity,
       sobolev::space::L2, discontinuous, degree, degree,
@@ -1243,7 +1257,8 @@ FiniteElement<T> element::create_serendipity_div(
               ? element::create_lagrange<T>(facettype, degree, lvariant, true)
               : element::create_dpc<T>(facettype, degree, dvariant, true);
     auto [_x, xshape, _M, Mshape] = moments::make_normal_integral_moments<T>(
-        facet_moment_space, celltype, tdim, 2 * degree + 1);
+        facet_moment_space, celltype, polyset::type::standard, tdim,
+        2 * degree + 1);
     assert(_x.size() == _M.size());
     for (std::size_t i = 0; i < _x.size(); ++i)
     {
@@ -1258,7 +1273,8 @@ FiniteElement<T> element::create_serendipity_div(
     FiniteElement cell_moment_space
         = element::create_dpc<T>(celltype, degree - 2, dvariant, true);
     auto [_x, xshape, _M, Mshape] = moments::make_integral_moments<T>(
-        cell_moment_space, celltype, tdim, 2 * degree - 1);
+        cell_moment_space, celltype, polyset::type::standard, tdim,
+        2 * degree - 1);
     assert(_x.size() == _M.size());
     for (std::size_t i = 0; i < _x.size(); ++i)
     {
@@ -1306,11 +1322,11 @@ FiniteElement<T> element::create_serendipity_div(
 
   sobolev::space space
       = discontinuous ? sobolev::space::L2 : sobolev::space::HDiv;
-  return FiniteElement<T>(element::family::BDM, celltype, degree, {tdim},
-                          impl::mdspan_t<const T, 2>(wbuffer.data(), wshape),
-                          xview, Mview, 0, maps::type::contravariantPiola,
-                          space, discontinuous, degree / tdim, degree + 1,
-                          lvariant, dvariant);
+  return FiniteElement<T>(
+      element::family::BDM, celltype, polyset::type::standard, degree, {tdim},
+      impl::mdspan_t<const T, 2>(wbuffer.data(), wshape), xview, Mview, 0,
+      maps::type::contravariantPiola, space, discontinuous, degree / tdim,
+      degree + 1, lvariant, dvariant);
 }
 //-----------------------------------------------------------------------------
 template <std::floating_point T>
@@ -1330,8 +1346,9 @@ FiniteElement<T> element::create_serendipity_curl(
   const std::size_t tdim = cell::topological_dimension(celltype);
 
   // Evaluate the expansion polynomials at the quadrature points
-  const auto [_Qpts, wts] = quadrature::make_quadrature<T>(
-      quadrature::type::Default, celltype, 2 * degree + 1);
+  const auto [_Qpts, wts]
+      = quadrature::make_quadrature<T>(quadrature::type::Default, celltype,
+                                       polyset::type::standard, 2 * degree + 1);
   impl::mdspan_t<const T, 2> Qpts(_Qpts.data(), wts.size(),
                                   _Qpts.size() / wts.size());
 
@@ -1366,7 +1383,8 @@ FiniteElement<T> element::create_serendipity_curl(
     FiniteElement edge_moment_space = element::create_lagrange<T>(
         cell::type::interval, degree, lvariant, true);
     auto [_x, xshape, _M, Mshape] = moments::make_tangent_integral_moments<T>(
-        edge_moment_space, celltype, tdim, 2 * degree + 1);
+        edge_moment_space, celltype, polyset::type::standard, tdim,
+        2 * degree + 1);
     assert(_x.size() == _M.size());
     for (std::size_t i = 0; i < _x.size(); ++i)
     {
@@ -1381,7 +1399,7 @@ FiniteElement<T> element::create_serendipity_curl(
     FiniteElement moment_space = element::create_dpc<T>(
         cell::type::quadrilateral, degree - 2, dvariant, true);
     auto [_x, xshape, _M, Mshape] = moments::make_integral_moments<T>(
-        moment_space, celltype, tdim, 2 * degree - 1);
+        moment_space, celltype, polyset::type::standard, tdim, 2 * degree - 1);
     assert(_x.size() == _M.size());
     for (std::size_t i = 0; i < _x.size(); ++i)
     {
@@ -1404,7 +1422,7 @@ FiniteElement<T> element::create_serendipity_curl(
       auto [_x, xshape, _M, Mshape] = moments::make_integral_moments<T>(
           element::create_dpc<T>(cell::type::hexahedron, degree - 4, dvariant,
                                  true),
-          celltype, tdim, 2 * degree - 3);
+          celltype, polyset::type::standard, tdim, 2 * degree - 3);
       assert(_x.size() == _M.size());
       for (std::size_t i = 0; i < _x.size(); ++i)
       {
@@ -1436,12 +1454,12 @@ FiniteElement<T> element::create_serendipity_curl(
 
   sobolev::space space
       = discontinuous ? sobolev::space::L2 : sobolev::space::HCurl;
-  return FiniteElement<T>(element::family::N2E, celltype, degree, {tdim},
-                          impl::mdspan_t<const T, 2>(wbuffer.data(), wshape),
-                          xview, Mview, 0, maps::type::covariantPiola, space,
-                          discontinuous,
-                          (degree == 2 && tdim == 3) ? 1 : degree / tdim,
-                          degree + 1, lvariant, dvariant);
+  return FiniteElement<T>(
+      element::family::N2E, celltype, polyset::type::standard, degree, {tdim},
+      impl::mdspan_t<const T, 2>(wbuffer.data(), wshape), xview, Mview, 0,
+      maps::type::covariantPiola, space, discontinuous,
+      (degree == 2 && tdim == 3) ? 1 : degree / tdim, degree + 1, lvariant,
+      dvariant);
 }
 //-----------------------------------------------------------------------------
 template FiniteElement<float>
diff --git a/cpp/basix/element-families.h b/cpp/basix/element-families.h
index 1a92d35e9..c73516aa9 100644
--- a/cpp/basix/element-families.h
+++ b/cpp/basix/element-families.h
@@ -56,6 +56,7 @@ enum class family
   bubble = 9,
   serendipity = 10,
   HHJ = 11,
-  Hermite = 12
+  Hermite = 12,
+  iso = 13,
 };
 } // namespace basix::element
diff --git a/cpp/basix/finite-element.cpp b/cpp/basix/finite-element.cpp
index e1a22318c..f942e88c6 100644
--- a/cpp/basix/finite-element.cpp
+++ b/cpp/basix/finite-element.cpp
@@ -107,7 +107,7 @@ constexpr int num_transformations(cell::type cell_type)
 //-----------------------------------------------------------------------------
 template <std::floating_point T>
 std::pair<std::vector<T>, std::array<std::size_t, 2>> compute_dual_matrix(
-    cell::type cell_type, mdspan_t<const T, 2> B,
+    cell::type cell_type, polyset::type poly_type, mdspan_t<const T, 2> B,
     const std::array<std::vector<impl::mdspan_t<const T, 2>>, 4>& x,
     const std::array<std::vector<impl::mdspan_t<const T, 4>>, 4>& M, int degree,
     int nderivs)
@@ -125,7 +125,7 @@ std::pair<std::vector<T>, std::array<std::size_t, 2>> compute_dual_matrix(
     }
   }
 
-  std::size_t pdim = polyset::dim(cell_type, degree);
+  std::size_t pdim = polyset::dim(cell_type, poly_type, degree);
   mdarray_t<T, 3> D(vs, pdim, num_dofs);
   std::fill(D.data(), D.data() + D.size(), 0);
   std::vector<T> Pb;
@@ -144,7 +144,7 @@ std::pair<std::vector<T>, std::array<std::size_t, 2>> compute_dual_matrix(
       {
         std::array<std::size_t, 3> shape;
         std::tie(Pb, shape)
-            = polyset::tabulate(cell_type, degree, nderivs, x_e);
+            = polyset::tabulate(cell_type, poly_type, degree, nderivs, x_e);
         P = mdspan_t<const T, 3>(Pb.data(), shape);
       }
 
@@ -233,7 +233,8 @@ basix::create_element(element::family family, cell::type cell, int degree,
          {element::family::HHJ, {false, false}},
          {element::family::CR, {false, false}},
          {element::family::bubble, {false, false}},
-         {element::family::Hermite, {false, false}}};
+         {element::family::Hermite, {false, false}},
+         {element::family::iso, {true, false}}};
   if (auto it = has_variant.find(family); it != has_variant.end())
   {
     if (it->second[0] == false and lvariant != element::lagrange_variant::unset)
@@ -322,6 +323,8 @@ basix::create_element(element::family family, cell::type cell, int degree,
     return element::create_cr<T>(cell, degree, discontinuous);
   case element::family::bubble:
     return element::create_bubble<T>(cell, degree, discontinuous);
+  case element::family::iso:
+    return element::create_iso<T>(cell, degree, lvariant, discontinuous);
   case element::family::Hermite:
     return element::create_hermite<T>(cell, degree, discontinuous);
   default:
@@ -435,10 +438,10 @@ FiniteElement<T> basix::create_custom_element(
     const std::array<std::vector<impl::mdspan_t<const T, 4>>, 4>& M,
     int interpolation_nderivs, maps::type map_type,
     sobolev::space sobolev_space, bool discontinuous,
-    int highest_complete_degree, int highest_degree)
+    int highest_complete_degree, int highest_degree, polyset::type poly_type)
 {
   // Check that inputs are valid
-  const std::size_t psize = polyset::dim(cell_type, highest_degree);
+  const std::size_t psize = polyset::dim(cell_type, poly_type, highest_degree);
   const std::size_t value_size = std::reduce(
       value_shape.begin(), value_shape.end(), 1, std::multiplies{});
   const std::size_t deriv_count
@@ -501,8 +504,9 @@ FiniteElement<T> basix::create_custom_element(
     }
   }
 
-  auto [dualmatrix, dualshape] = compute_dual_matrix(
-      cell_type, wcoeffs, x, M, highest_degree, interpolation_nderivs);
+  auto [dualmatrix, dualshape]
+      = compute_dual_matrix(cell_type, poly_type, wcoeffs, x, M, highest_degree,
+                            interpolation_nderivs);
   if (math::is_singular(mdspan_t<const T, 2>(dualmatrix.data(), dualshape)))
   {
     throw std::runtime_error(
@@ -510,10 +514,10 @@ FiniteElement<T> basix::create_custom_element(
   }
 
   return basix::FiniteElement<T>(
-      element::family::custom, cell_type, highest_degree, value_shape, wcoeffs,
-      x, M, interpolation_nderivs, map_type, sobolev_space, discontinuous,
-      highest_complete_degree, highest_degree, element::lagrange_variant::unset,
-      element::dpc_variant::unset);
+      element::family::custom, cell_type, poly_type, highest_degree,
+      value_shape, wcoeffs, x, M, interpolation_nderivs, map_type,
+      sobolev_space, discontinuous, highest_complete_degree, highest_degree,
+      element::lagrange_variant::unset, element::dpc_variant::unset);
 }
 //-----------------------------------------------------------------------------
 /// @cond
@@ -522,20 +526,21 @@ template FiniteElement<float> basix::create_custom_element(
     impl::mdspan_t<const float, 2> wcoeffs,
     const std::array<std::vector<impl::mdspan_t<const float, 2>>, 4>&,
     const std::array<std::vector<impl::mdspan_t<const float, 4>>, 4>&, int,
-    maps::type, sobolev::space sobolev_space, bool, int, int);
+    maps::type, sobolev::space sobolev_space, bool, int, int, polyset::type);
 template FiniteElement<double> basix::create_custom_element(
     cell::type, const std::vector<std::size_t>&,
     impl::mdspan_t<const double, 2> wcoeffs,
     const std::array<std::vector<impl::mdspan_t<const double, 2>>, 4>&,
     const std::array<std::vector<impl::mdspan_t<const double, 4>>, 4>&, int,
-    maps::type, sobolev::space sobolev_space, bool, int, int);
+    maps::type, sobolev::space sobolev_space, bool, int, int, polyset::type);
 /// @endcond
 //-----------------------------------------------------------------------------
 /// @cond
 template <std::floating_point F>
 FiniteElement<F>::FiniteElement(
-    element::family family, cell::type cell_type, int degree,
-    const std::vector<std::size_t>& value_shape, mdspan_t<const F, 2> wcoeffs,
+    element::family family, cell::type cell_type, polyset::type poly_type,
+    int degree, const std::vector<std::size_t>& value_shape,
+    mdspan_t<const F, 2> wcoeffs,
     const std::array<std::vector<mdspan_t<const F, 2>>, 4>& x,
     const std::array<std::vector<mdspan_t<const F, 4>>, 4>& M,
     int interpolation_nderivs, maps::type map_type,
@@ -545,7 +550,8 @@ FiniteElement<F>::FiniteElement(
     std::vector<std::tuple<std::vector<FiniteElement<F>>, std::vector<int>>>
         tensor_factors,
     std::vector<int> dof_ordering)
-    : _cell_type(cell_type), _cell_tdim(cell::topological_dimension(cell_type)),
+    : _cell_type(cell_type), _poly_type(poly_type),
+      _cell_tdim(cell::topological_dimension(cell_type)),
       _cell_subentity_types(cell::subentity_types(cell_type)), _family(family),
       _lagrange_variant(lvariant), _dpc_variant(dvariant), _degree(degree),
       _interpolation_nderivs(interpolation_nderivs),
@@ -589,8 +595,9 @@ FiniteElement<F>::FiniteElement(
             wcoeffs_ortho_b.begin());
   if (family != element::family::P)
     orthogonalise(wcoeffs_ortho);
-  _dual_matrix = compute_dual_matrix<F>(cell_type, wcoeffs_ortho, x, M,
-                                        highest_degree, interpolation_nderivs);
+  _dual_matrix
+      = compute_dual_matrix<F>(cell_type, poly_type, wcoeffs_ortho, x, M,
+                               highest_degree, interpolation_nderivs);
 
   _wcoeffs
       = {wcoeffs_ortho_b, {wcoeffs_ortho.extent(0), wcoeffs_ortho.extent(1)}};
@@ -1029,12 +1036,13 @@ void FiniteElement<F>::tabulate(int nd, impl::mdspan_t<const F, 2> x,
                              + std::to_string(_cell_tdim) + ").");
   }
 
-  const std::size_t psize = polyset::dim(_cell_type, _highest_degree);
+  const std::size_t psize
+      = polyset::dim(_cell_type, _poly_type, _highest_degree);
   const std::array<std::size_t, 3> bsize
       = {(std::size_t)polyset::nderivs(_cell_type, nd), psize, x.extent(0)};
   std::vector<F> basis_b(bsize[0] * bsize[1] * bsize[2]);
   mdspan_t<F, 3> basis(basis_b.data(), bsize);
-  polyset::tabulate(basis, _cell_type, _highest_degree, nd, x);
+  polyset::tabulate(basis, _cell_type, _poly_type, _highest_degree, nd, x);
   const int vs = std::accumulate(_value_shape.begin(), _value_shape.end(), 1,
                                  std::multiplies{});
 
diff --git a/cpp/basix/finite-element.h b/cpp/basix/finite-element.h
index 64d74f102..15ad4905b 100644
--- a/cpp/basix/finite-element.h
+++ b/cpp/basix/finite-element.h
@@ -8,6 +8,7 @@
 #include "element-families.h"
 #include "maps.h"
 #include "mdspan.hpp"
+#include "polyset.h"
 #include "precompute.h"
 #include "sobolev-spaces.h"
 #include <array>
@@ -284,6 +285,7 @@ class FiniteElement
   ///
   /// @param[in] family The element family
   /// @param[in] cell_type The cell type
+  /// @param[in] poly_type The polyset type
   /// @param[in] degree The degree of the element
   /// @param[in] interpolation_nderivs The number of derivatives that
   /// need to be used during interpolation
@@ -315,8 +317,9 @@ class FiniteElement
   /// @param[in] dof_ordering DOF reordering: a mapping from the
   /// reference order to a new permuted order
   FiniteElement(
-      element::family family, cell::type cell_type, int degree,
-      const std::vector<std::size_t>& value_shape, mdspan_t<const F, 2> wcoeffs,
+      element::family family, cell::type cell_type, polyset::type poly_type,
+      int degree, const std::vector<std::size_t>& value_shape,
+      mdspan_t<const F, 2> wcoeffs,
       const std::array<std::vector<mdspan_t<const F, 2>>, 4>& x,
       const std::array<std::vector<mdspan_t<const F, 4>>, 4>& M,
       int interpolation_nderivs, maps::type map_type,
@@ -483,6 +486,10 @@ class FiniteElement
   /// @return The cell type
   cell::type cell_type() const { return _cell_type; }
 
+  /// Get the element polyset type
+  /// @return The polyset
+  polyset::type polyset_type() const { return _poly_type; }
+
   /// Get the element polynomial degree
   /// @return Polynomial degree
   int degree() const { return _degree; }
@@ -1152,6 +1159,9 @@ class FiniteElement
   // Cell type
   cell::type _cell_type;
 
+  // Polyset type
+  polyset::type _poly_type;
+
   // Topological dimension of the cell
   std::size_t _cell_tdim;
 
@@ -1311,6 +1321,7 @@ class FiniteElement
 /// element
 /// @param[in] highest_degree The degree of a polynomial in this element's
 /// polyset
+/// @param[in] poly_type The type of polyset to use for this element
 /// @return A custom finite element
 template <std::floating_point T>
 FiniteElement<T> create_custom_element(
@@ -1320,7 +1331,7 @@ FiniteElement<T> create_custom_element(
     const std::array<std::vector<impl::mdspan_t<const T, 4>>, 4>& M,
     int interpolation_nderivs, maps::type map_type,
     sobolev::space sobolev_space, bool discontinuous,
-    int highest_complete_degree, int highest_degree);
+    int highest_complete_degree, int highest_degree, polyset::type poly_type);
 
 /// Create an element using a given Lagrange variant and a given DPC variant
 /// @param[in] family The element family
diff --git a/cpp/basix/lattice.cpp b/cpp/basix/lattice.cpp
index 05daa820f..da4a4919e 100644
--- a/cpp/basix/lattice.cpp
+++ b/cpp/basix/lattice.cpp
@@ -173,8 +173,8 @@ tabulate_dlagrange(std::size_t n, std::span<const T> x)
   stdex::mdspan<const T, stdex::dextents<std::size_t, 2>> equi_pts_v(
       equi_pts.data(), n + 1, 1);
 
-  const auto [dual_values_b, dshape]
-      = polyset::tabulate<T>(cell::type::interval, n, 0, equi_pts_v);
+  const auto [dual_values_b, dshape] = polyset::tabulate<T>(
+      cell::type::interval, polyset::type::standard, n, 0, equi_pts_v);
   stdex::mdspan<const T, stdex::dextents<std::size_t, 3>> dual_values(
       dual_values_b.data(), dshape);
 
@@ -188,8 +188,9 @@ tabulate_dlagrange(std::size_t n, std::span<const T> x)
   using cmdspan2_t
       = stdex::mdspan<const T,
                       stdex::extents<std::size_t, stdex::dynamic_extent, 1>>;
-  const auto [tabulated_values_b, tshape] = polyset::tabulate<T>(
-      cell::type::interval, n, 0, cmdspan2_t(x.data(), x.size(), 1));
+  const auto [tabulated_values_b, tshape]
+      = polyset::tabulate<T>(cell::type::interval, polyset::type::standard, n,
+                             0, cmdspan2_t(x.data(), x.size(), 1));
   stdex::mdspan<const T, stdex::dextents<std::size_t, 3>> tabulated_values(
       tabulated_values_b.data(), tshape);
 
diff --git a/cpp/basix/moments.cpp b/cpp/basix/moments.cpp
index 56d23419b..5739dc360 100644
--- a/cpp/basix/moments.cpp
+++ b/cpp/basix/moments.cpp
@@ -100,7 +100,8 @@ template <std::floating_point T>
 std::tuple<std::vector<std::vector<T>>, std::array<std::size_t, 2>,
            std::vector<std::vector<T>>, std::array<std::size_t, 4>>
 moments::make_integral_moments(const FiniteElement<T>& V, cell::type celltype,
-                               std::size_t value_size, int q_deg)
+                               polyset::type ptype, std::size_t value_size,
+                               int q_deg)
 {
   const cell::type sub_celltype = V.cell_type();
   const std::size_t entity_dim = cell::topological_dimension(sub_celltype);
@@ -110,7 +111,10 @@ moments::make_integral_moments(const FiniteElement<T>& V, cell::type celltype,
 
   // Get the quadrature points and weights
   const auto [_pts, wts] = quadrature::make_quadrature<T>(
-      quadrature::type::Default, sub_celltype, q_deg);
+      quadrature::type::Default, sub_celltype,
+      polyset::superset(sub_celltype, V.polyset_type(),
+                        polyset::restriction(ptype, celltype, sub_celltype)),
+      q_deg);
   mdspan_t<const T, 2> pts(_pts.data(), wts.size(), _pts.size() / wts.size());
 
   // Evaluate moment space at quadrature points
@@ -188,15 +192,18 @@ template <std::floating_point T>
 std::tuple<std::vector<std::vector<T>>, std::array<std::size_t, 2>,
            std::vector<std::vector<T>>, std::array<std::size_t, 4>>
 moments::make_dot_integral_moments(const FiniteElement<T>& V,
-                                   cell::type celltype, std::size_t value_size,
-                                   int q_deg)
+                                   cell::type celltype, polyset::type ptype,
+                                   std::size_t value_size, int q_deg)
 {
   const cell::type sub_celltype = V.cell_type();
   const std::size_t entity_dim = cell::topological_dimension(sub_celltype);
   const std::size_t num_entities = cell::num_sub_entities(celltype, entity_dim);
 
   const auto [_pts, wts] = quadrature::make_quadrature<T>(
-      quadrature::type::Default, sub_celltype, q_deg);
+      quadrature::type::Default, sub_celltype,
+      polyset::superset(sub_celltype, V.polyset_type(),
+                        polyset::restriction(ptype, celltype, sub_celltype)),
+      q_deg);
   mdspan_t<const T, 2> pts(_pts.data(), wts.size(), _pts.size() / wts.size());
 
   // If this is always true, value_size input can be removed
@@ -262,7 +269,7 @@ template <std::floating_point T>
 std::tuple<std::vector<std::vector<T>>, std::array<std::size_t, 2>,
            std::vector<std::vector<T>>, std::array<std::size_t, 4>>
 moments::make_tangent_integral_moments(const FiniteElement<T>& V,
-                                       cell::type celltype,
+                                       cell::type celltype, polyset::type ptype,
                                        std::size_t value_size, int q_deg)
 {
   const cell::type sub_celltype = V.cell_type();
@@ -277,7 +284,10 @@ moments::make_tangent_integral_moments(const FiniteElement<T>& V,
     throw std::runtime_error("Tangent is only well-defined on an edge.");
 
   const auto [_pts, wts] = quadrature::make_quadrature<T>(
-      quadrature::type::Default, cell::type::interval, q_deg);
+      quadrature::type::Default, cell::type::interval,
+      polyset::superset(sub_celltype, V.polyset_type(),
+                        polyset::restriction(ptype, celltype, sub_celltype)),
+      q_deg);
   mdspan_t<const T, 2> pts(_pts.data(), wts.size(), _pts.size() / wts.size());
 
   // Evaluate moment space at quadrature points
@@ -334,7 +344,7 @@ template <std::floating_point T>
 std::tuple<std::vector<std::vector<T>>, std::array<std::size_t, 2>,
            std::vector<std::vector<T>>, std::array<std::size_t, 4>>
 moments::make_normal_integral_moments(const FiniteElement<T>& V,
-                                      cell::type celltype,
+                                      cell::type celltype, polyset::type ptype,
                                       std::size_t value_size, int q_deg)
 {
   const std::size_t tdim = cell::topological_dimension(celltype);
@@ -348,7 +358,10 @@ moments::make_normal_integral_moments(const FiniteElement<T>& V,
 
   // Compute quadrature points for evaluating integral
   const auto [_pts, wts] = quadrature::make_quadrature<T>(
-      quadrature::type::Default, sub_celltype, q_deg);
+      quadrature::type::Default, sub_celltype,
+      polyset::superset(sub_celltype, V.polyset_type(),
+                        polyset::restriction(ptype, celltype, sub_celltype)),
+      q_deg);
   mdspan_t<const T, 2> pts(_pts.data(), wts.size(), _pts.size() / wts.size());
 
   // Evaluate moment space at quadrature points
@@ -429,41 +442,41 @@ moments::make_normal_integral_moments(const FiniteElement<T>& V,
 template std::tuple<std::vector<std::vector<float>>, std::array<std::size_t, 2>,
                     std::vector<std::vector<float>>, std::array<std::size_t, 4>>
 moments::make_integral_moments(const FiniteElement<float>&, cell::type,
-                               std::size_t, int);
+                               polyset::type, std::size_t, int);
 template std::tuple<
     std::vector<std::vector<double>>, std::array<std::size_t, 2>,
     std::vector<std::vector<double>>, std::array<std::size_t, 4>>
 moments::make_integral_moments(const FiniteElement<double>&, cell::type,
-                               std::size_t, int);
+                               polyset::type, std::size_t, int);
 
 template std::tuple<std::vector<std::vector<float>>, std::array<std::size_t, 2>,
                     std::vector<std::vector<float>>, std::array<std::size_t, 4>>
 moments::make_dot_integral_moments(const FiniteElement<float>&, cell::type,
-                                   std::size_t, int);
+                                   polyset::type, std::size_t, int);
 template std::tuple<
     std::vector<std::vector<double>>, std::array<std::size_t, 2>,
     std::vector<std::vector<double>>, std::array<std::size_t, 4>>
 moments::make_dot_integral_moments(const FiniteElement<double>&, cell::type,
-                                   std::size_t, int);
+                                   polyset::type, std::size_t, int);
 
 template std::tuple<std::vector<std::vector<float>>, std::array<std::size_t, 2>,
                     std::vector<std::vector<float>>, std::array<std::size_t, 4>>
 moments::make_tangent_integral_moments(const FiniteElement<float>&, cell::type,
-                                       std::size_t, int);
+                                       polyset::type, std::size_t, int);
 template std::tuple<
     std::vector<std::vector<double>>, std::array<std::size_t, 2>,
     std::vector<std::vector<double>>, std::array<std::size_t, 4>>
 moments::make_tangent_integral_moments(const FiniteElement<double>&, cell::type,
-                                       std::size_t, int);
+                                       polyset::type, std::size_t, int);
 
 template std::tuple<std::vector<std::vector<float>>, std::array<std::size_t, 2>,
                     std::vector<std::vector<float>>, std::array<std::size_t, 4>>
 moments::make_normal_integral_moments(const FiniteElement<float>&, cell::type,
-                                      std::size_t, int);
+                                      polyset::type, std::size_t, int);
 template std::tuple<
     std::vector<std::vector<double>>, std::array<std::size_t, 2>,
     std::vector<std::vector<double>>, std::array<std::size_t, 4>>
 moments::make_normal_integral_moments(const FiniteElement<double>&, cell::type,
-                                      std::size_t, int);
+                                      polyset::type, std::size_t, int);
 /// @endcond
 //----------------------------------------------------------------------------
diff --git a/cpp/basix/moments.h b/cpp/basix/moments.h
index b1caa8460..98d206980 100644
--- a/cpp/basix/moments.h
+++ b/cpp/basix/moments.h
@@ -5,6 +5,7 @@
 #pragma once
 
 #include "cell.h"
+#include "polyset.h"
 #include <array>
 #include <concepts>
 #include <tuple>
@@ -33,6 +34,8 @@ namespace moments
 /// against
 /// @param celltype The cell type of the cell on which the space is
 /// being defined
+/// @param ptype The polyset type of the element this moment is being used to
+/// define
 /// @param value_size The value size of the space being defined
 /// @param q_deg The quadrature degree used for the integrals
 /// @return (interpolation points, interpolation matrix). The indices of
@@ -43,7 +46,7 @@ template <std::floating_point T>
 std::tuple<std::vector<std::vector<T>>, std::array<std::size_t, 2>,
            std::vector<std::vector<T>>, std::array<std::size_t, 4>>
 make_integral_moments(const FiniteElement<T>& moment_space, cell::type celltype,
-                      std::size_t value_size, int q_deg);
+                      polyset::type ptype, std::size_t value_size, int q_deg);
 
 /// @brief Make interpolation points and weights for dot product
 /// integral moments.
@@ -60,6 +63,8 @@ make_integral_moments(const FiniteElement<T>& moment_space, cell::type celltype,
 /// @param V The space to compute the integral moments against
 /// @param celltype The cell type of the cell on which the space is
 /// being defined
+/// @param ptype The polyset type of the element this moment is being used to
+/// define
 /// @param value_size The value size of the space being defined
 /// @param q_deg The quadrature degree used for the integrals
 /// @return (interpolation points, interpolation shape,  interpolation
@@ -71,7 +76,8 @@ template <std::floating_point T>
 std::tuple<std::vector<std::vector<T>>, std::array<std::size_t, 2>,
            std::vector<std::vector<T>>, std::array<std::size_t, 4>>
 make_dot_integral_moments(const FiniteElement<T>& V, cell::type celltype,
-                          std::size_t value_size, int q_deg);
+                          polyset::type ptype, std::size_t value_size,
+                          int q_deg);
 
 /// @brief Make interpolation points and weights for tangent integral
 /// moments.
@@ -82,6 +88,8 @@ make_dot_integral_moments(const FiniteElement<T>& V, cell::type celltype,
 /// @param V The space to compute the integral moments against
 /// @param celltype The cell type of the cell on which the space is
 /// being defined
+/// @param ptype The polyset type of the element this moment is being used to
+/// define
 /// @param value_size The value size of the space being defined the
 /// space
 /// @param q_deg The quadrature degree used for the integrals
@@ -94,7 +102,8 @@ template <std::floating_point T>
 std::tuple<std::vector<std::vector<T>>, std::array<std::size_t, 2>,
            std::vector<std::vector<T>>, std::array<std::size_t, 4>>
 make_tangent_integral_moments(const FiniteElement<T>& V, cell::type celltype,
-                              std::size_t value_size, int q_deg);
+                              polyset::type ptype, std::size_t value_size,
+                              int q_deg);
 
 ///  @brief Compute interpolation points and weights for normal integral
 ///  moments.
@@ -105,6 +114,8 @@ make_tangent_integral_moments(const FiniteElement<T>& V, cell::type celltype,
 /// @param[in] V The space to compute the integral moments against
 /// @param[in] celltype The cell type of the cell on which the space is
 /// being defined
+/// @param ptype The polyset type of the element this moment is being used to
+/// define
 /// @param[in] value_size The value size of the space being defined
 /// @param[in] q_deg The quadrature degree used for the integrals
 /// @return (interpolation points, interpolation shape,  interpolation
@@ -116,7 +127,8 @@ template <std::floating_point T>
 std::tuple<std::vector<std::vector<T>>, std::array<std::size_t, 2>,
            std::vector<std::vector<T>>, std::array<std::size_t, 4>>
 make_normal_integral_moments(const FiniteElement<T>& V, cell::type celltype,
-                             std::size_t value_size, int q_deg);
+                             polyset::type ptype, std::size_t value_size,
+                             int q_deg);
 
 } // namespace moments
 } // namespace basix
diff --git a/cpp/basix/polynomials.cpp b/cpp/basix/polynomials.cpp
index 4d1c8c34e..31c42c0c9 100644
--- a/cpp/basix/polynomials.cpp
+++ b/cpp/basix/polynomials.cpp
@@ -23,9 +23,9 @@ using mdspan_t = stdex::mdspan<T, stdex::dextents<std::size_t, d>>;
 constexpr int single_choose(int n, int k)
 {
   int out = 1;
-  for (int i = k + 1; i <= n; ++i)
+  for (int i = n + 1 - k; i <= n; ++i)
     out *= i;
-  for (int i = 1; i <= n - k; ++i)
+  for (int i = 1; i <= k; ++i)
     out /= i;
   return out;
 }
@@ -119,7 +119,8 @@ std::pair<std::vector<T>, std::array<std::size_t, 2>> polynomials::tabulate(
   {
   case polynomials::type::legendre:
   {
-    auto [values, shape] = polyset::tabulate(celltype, d, 0, x);
+    auto [values, shape]
+        = polyset::tabulate(celltype, polyset::type::standard, d, 0, x);
     assert(shape[0] == 1);
     return {std::move(values), {shape[1], shape[2]}};
   }
@@ -135,7 +136,7 @@ std::pair<std::vector<T>, std::array<std::size_t, 2>> polynomials::tabulate(
 //-----------------------------------------------------------------------------
 int polynomials::dim(polynomials::type, cell::type cell, int d)
 {
-  return polyset::dim(cell, d);
+  return polyset::dim(cell, polyset::type::standard, d);
 }
 //-----------------------------------------------------------------------------
 /// @cond
diff --git a/cpp/basix/polyset.cpp b/cpp/basix/polyset.cpp
index eda7d846e..f1159b1e0 100644
--- a/cpp/basix/polyset.cpp
+++ b/cpp/basix/polyset.cpp
@@ -17,6 +17,17 @@ namespace stdex = std::experimental;
 
 namespace
 {
+//-----------------------------------------------------------------------------
+constexpr int single_choose(int n, int k)
+{
+  int out = 1;
+  for (int i = n + 1 - k; i <= n; ++i)
+    out *= i;
+  for (int i = 1; i <= k; ++i)
+    out /= i;
+  return out;
+}
+//-----------------------------------------------------------------------------
 /// Compute coefficients in the Jacobi Polynomial recurrence relation
 template <typename T>
 constexpr std::array<T, 3> jrc(int a, int n)
@@ -112,6 +123,105 @@ void tabulate_polyset_line_derivs(
 }
 //-----------------------------------------------------------------------------
 
+/// Compute the complete set of derivatives from 0 to nderiv, for all
+/// the polynomials up to order n on a line segment. The polynomials
+/// used are Legendre Polynomials, with the recurrence relation given by
+/// n P(n) = (2n - 1) x P_{n-1} - (n - 1) P_{n-2} in the interval [-1,
+/// 1]. The range is rescaled here to [0, 1].
+template <typename T>
+void tabulate_polyset_line_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);
+  assert(P.extent(2) == x.extent(0));
+
+  auto x0 = stdex::submdspan(x, stdex::full_extent, 0);
+
+  std::fill(P.data_handle(), P.data_handle() + P.size(), 0.0);
+
+  std::vector<T> factorials(n + 1, 0.0);
+
+  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 d = 0; d <= nderiv; ++d)
+  {
+    for (std::size_t p = 0; p < P.extent(2); ++p)
+    {
+      if (x0[p] <= 0.5)
+      {
+        for (std::size_t k = 0; k + d <= n; ++k)
+        {
+          T x_term = pow(x0[p], n - k - d);
+          for (std::size_t i = n - k; i > n - k - d; --i)
+            x_term *= i;
+          P(d, 0, p) += factorials[k] * x_term;
+        }
+      }
+      else
+      {
+        for (std::size_t k = 0; k + d <= n; ++k)
+        {
+          T x_term = pow(1.0 - x0[p], n - k - d);
+          for (std::size_t i = n - k; i > n - k - d; --i)
+            x_term *= -i;
+          P(d, 0, p) += factorials[k] * x_term;
+        }
+      }
+    }
+  }
+
+  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 d = 0; d <= nderiv; ++d)
+    {
+      for (std::size_t p = 0; p < P.extent(2); ++p)
+      {
+        if (x0[p] <= 0.5)
+        {
+          for (std::size_t k = 0; k + d <= j; ++k)
+          {
+            T x_term = pow(x0[p], j - k - d);
+            for (std::size_t i = j - k; i > j - k - d; --i)
+              x_term *= i;
+            P(d, j + 1, p) += factorials[k] * x_term;
+          }
+          P(d, j + 1, p) *= pow(0.5 - x0[p], n - j - d);
+          for (std::size_t i = n - j; i > n - j - d; --i)
+            P(d, j + 1, p) *= -i;
+        }
+        else
+        {
+          for (std::size_t k = 0; k + d <= j; ++k)
+          {
+            T x_term = pow(1.0 - x0[p], j - k - d);
+            for (std::size_t i = j - k; i > j - k - d; --i)
+              x_term *= -i;
+            P(d, j + n + 1, p) += factorials[k] * x_term;
+          }
+          P(d, j + n + 1, p) *= pow(x0[p] - 0.5, n - j - d);
+          for (std::size_t i = n - j; i > n - j - d; --i)
+            P(d, j + n + 1, p) *= i;
+        }
+      }
+    }
+  }
+}
+//-----------------------------------------------------------------------------
+
 /// 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
 /// polynomials P_{pq} are built up in sequence, firstly along q = 0,
@@ -1480,91 +1590,125 @@ void tabulate_polyset_prism_derivs(
 template <std::floating_point T>
 void polyset::tabulate(
     std::experimental::mdspan<T, std::experimental::dextents<std::size_t, 3>> P,
-    cell::type celltype, int d, int n,
+    cell::type celltype, polyset::type ptype, int d, int n,
     std::experimental::mdspan<const T,
                               std::experimental::dextents<std::size_t, 2>>
         x)
 {
-  switch (celltype)
+  switch (ptype)
   {
-  case cell::type::point:
-    tabulate_polyset_point_derivs(P, d, n, x);
-    return;
-  case cell::type::interval:
-    tabulate_polyset_line_derivs(P, d, n, x);
-    return;
-  case cell::type::triangle:
-    tabulate_polyset_triangle_derivs(P, d, n, x);
-    return;
-  case cell::type::tetrahedron:
-    tabulate_polyset_tetrahedron_derivs(P, d, n, x);
-    return;
-  case cell::type::quadrilateral:
-    tabulate_polyset_quad_derivs(P, d, n, x);
-    return;
-  case cell::type::prism:
-    tabulate_polyset_prism_derivs(P, d, n, x);
-    return;
-  case cell::type::pyramid:
-    tabulate_polyset_pyramid_derivs(P, d, n, x);
-    return;
-  case cell::type::hexahedron:
-    tabulate_polyset_hex_derivs(P, d, n, x);
-    return;
+  case polyset::type::standard:
+    switch (celltype)
+    {
+    case cell::type::point:
+      tabulate_polyset_point_derivs(P, d, n, x);
+      return;
+    case cell::type::interval:
+      tabulate_polyset_line_derivs(P, d, n, x);
+      return;
+    case cell::type::triangle:
+      tabulate_polyset_triangle_derivs(P, d, n, x);
+      return;
+    case cell::type::tetrahedron:
+      tabulate_polyset_tetrahedron_derivs(P, d, n, x);
+      return;
+    case cell::type::quadrilateral:
+      tabulate_polyset_quad_derivs(P, d, n, x);
+      return;
+    case cell::type::prism:
+      tabulate_polyset_prism_derivs(P, d, n, x);
+      return;
+    case cell::type::pyramid:
+      tabulate_polyset_pyramid_derivs(P, d, n, x);
+      return;
+    case cell::type::hexahedron:
+      tabulate_polyset_hex_derivs(P, d, n, x);
+      return;
+    default:
+      throw std::runtime_error("Polynomial set: unsupported cell type");
+    }
+  case polyset::type::macroedge:
+    switch (celltype)
+    {
+    case cell::type::point:
+      tabulate_polyset_point_derivs(P, d, n, x);
+      return;
+    case cell::type::interval:
+      tabulate_polyset_line_macroedge_derivs(P, d, n, x);
+      return;
+    default:
+      throw std::runtime_error("Polynomial set: unsupported cell type");
+    }
   default:
-    throw std::runtime_error("Polynomial set: unsupported cell type");
+    throw std::runtime_error("Polynomial set: unsupported polynomial type.");
   }
 }
 //-----------------------------------------------------------------------------
 template <std::floating_point T>
 std::pair<std::vector<T>, std::array<std::size_t, 3>> polyset::tabulate(
-    cell::type celltype, int d, int n,
+    cell::type celltype, polyset::type ptype, int d, int n,
     std::experimental::mdspan<const T,
                               std::experimental::dextents<std::size_t, 2>>
         x)
 {
   std::array<std::size_t, 3> shape
       = {(std::size_t)polyset::nderivs(celltype, n),
-         (std::size_t)polyset::dim(celltype, d), x.extent(0)};
+         (std::size_t)polyset::dim(celltype, ptype, d), x.extent(0)};
   std::vector<T> P(shape[0] * shape[1] * shape[2]);
   stdex::mdspan<T, stdex::dextents<std::size_t, 3>> _P(P.data(), shape);
-  polyset::tabulate(_P, celltype, d, n, x);
+  polyset::tabulate(_P, celltype, ptype, d, n, x);
   return {std::move(P), std::move(shape)};
 }
 //-----------------------------------------------------------------------------
 /// @cond
 template std::pair<std::vector<float>, std::array<std::size_t, 3>>
 polyset::tabulate(
-    cell::type, int, int,
+    cell::type, polyset::type, int, int,
     std::experimental::mdspan<const float,
                               std::experimental::dextents<std::size_t, 2>>);
 template std::pair<std::vector<double>, std::array<std::size_t, 3>>
 polyset::tabulate(
-    cell::type, int, int,
+    cell::type, polyset::type, int, int,
     std::experimental::mdspan<const double,
                               std::experimental::dextents<std::size_t, 2>>);
 /// @endcond
 //-----------------------------------------------------------------------------
-int polyset::dim(cell::type celltype, int d)
+int polyset::dim(cell::type celltype, polyset::type ptype, int d)
 {
-  switch (celltype)
+  switch (ptype)
   {
-  case cell::type::point:
-    return 1;
-  case cell::type::triangle:
-    return (d + 1) * (d + 2) / 2;
-  case cell::type::tetrahedron:
-    return (d + 1) * (d + 2) * (d + 3) / 6;
-  case cell::type::prism:
-    return (d + 1) * (d + 1) * (d + 2) / 2;
-  case cell::type::pyramid:
-    return (d + 1) * (d + 2) * (2 * d + 3) / 6;
-  case cell::type::interval:
-    return (d + 1);
-  case cell::type::quadrilateral:
-    return (d + 1) * (d + 1);
-  case cell::type::hexahedron:
-    return (d + 1) * (d + 1) * (d + 1);
+  case polyset::type::standard:
+    switch (celltype)
+    {
+    case cell::type::point:
+      return 1;
+    case cell::type::triangle:
+      return (d + 1) * (d + 2) / 2;
+    case cell::type::tetrahedron:
+      return (d + 1) * (d + 2) * (d + 3) / 6;
+    case cell::type::prism:
+      return (d + 1) * (d + 1) * (d + 2) / 2;
+    case cell::type::pyramid:
+      return (d + 1) * (d + 2) * (2 * d + 3) / 6;
+    case cell::type::interval:
+      return (d + 1);
+    case cell::type::quadrilateral:
+      return (d + 1) * (d + 1);
+    case cell::type::hexahedron:
+      return (d + 1) * (d + 1) * (d + 1);
+    default:
+      return 1;
+    }
+  case polyset::type::macroedge:
+    switch (celltype)
+    {
+    case cell::type::point:
+      return 1;
+    case cell::type::interval:
+      return 2 * d + 1;
+    default:
+      return 1;
+    }
   default:
     return 1;
   }
@@ -1595,3 +1739,25 @@ int polyset::nderivs(cell::type celltype, int n)
   }
 }
 //-----------------------------------------------------------------------------
+polyset::type polyset::superset(cell::type, polyset::type type1,
+                                polyset::type type2)
+{
+  if (type1 == type2)
+    return type1;
+  if (type1 == polyset::type::standard)
+    return type2;
+  if (type2 == polyset::type::standard)
+    return type1;
+  throw std::runtime_error("Unsupported superset of polynomial sets.");
+}
+//-----------------------------------------------------------------------------
+polyset::type polyset::restriction(polyset::type ptype, cell::type cell,
+                                   cell::type restriction_cell)
+{
+  if (ptype == polyset::type::standard)
+    return polyset::type::standard;
+  if (cell == restriction_cell)
+    return ptype;
+  throw std::runtime_error("Unsupported restriction of polynomial sets.");
+}
+//-----------------------------------------------------------------------------
diff --git a/cpp/basix/polyset.h b/cpp/basix/polyset.h
index 704f7f55a..ba6731588 100644
--- a/cpp/basix/polyset.h
+++ b/cpp/basix/polyset.h
@@ -130,6 +130,14 @@
 /// quadratic terms, and cross-terms in two and three dimensions.
 namespace basix::polyset
 {
+
+/// Cell type
+enum class type
+{
+  standard = 0,
+  macroedge = 1,
+};
+
 /// @brief Tabulate the orthonormal polynomial basis, and derivatives,
 /// at points on the reference cell.
 ///
@@ -145,6 +153,7 @@ namespace basix::polyset
 /// 'triangular' with the lower derivatives appearing first.
 ///
 /// @param[in] celltype Cell type
+/// @param[in] ptype The polynomial type
 /// @param[in] d Polynomial degree
 /// @param[in] n Maximum derivative order. Use n = 0 for the basis only.
 /// @param[in] x Points at which to evaluate the basis. The shape is
@@ -169,7 +178,7 @@ namespace basix::polyset
 /// @todo Does the order for the third index need to be documented?
 template <std::floating_point T>
 std::pair<std::vector<T>, std::array<std::size_t, 3>>
-tabulate(cell::type celltype, int d, int n,
+tabulate(cell::type celltype, polyset::type ptype, int d, int n,
          std::experimental::mdspan<const T,
                                    std::experimental::dextents<std::size_t, 2>>
              x);
@@ -210,6 +219,7 @@ tabulate(cell::type celltype, int d, int n,
 /// - The third index is the basis function index.
 /// @todo Does the order for the third index need to be documented?
 /// @param[in] celltype Cell type
+/// @param[in] ptype The polynomial type
 /// @param[in] d Polynomial degree
 /// @param[in] n Maximum derivative order. Use n = 0 for the basis only.
 /// @param[in] x Points at which to evaluate the basis. The shape is
@@ -217,17 +227,18 @@ tabulate(cell::type celltype, int d, int n,
 template <std::floating_point T>
 void tabulate(
     std::experimental::mdspan<T, std::experimental::dextents<std::size_t, 3>> P,
-    cell::type celltype, int d, int n,
+    cell::type celltype, polyset::type ptype, int d, int n,
     std::experimental::mdspan<const T,
                               std::experimental::dextents<std::size_t, 2>>
         x);
 
 /// @brief Dimension of a polynomial space
 /// @param[in] cell The cell type
+/// @param[in] ptype The polynomial type
 /// @param[in] d The polynomial degree
 /// @return The number of terms in the basis spanning a space of
 /// polynomial degree @p d
-int dim(cell::type cell, int d);
+int dim(cell::type cell, polyset::type ptype, int d);
 
 /// @brief Number of derivatives that the orthonormal basis will have on
 /// the given cell.
@@ -237,4 +248,22 @@ int dim(cell::type cell, int d);
 /// polynomial degree @p d
 int nderivs(cell::type cell, int d);
 
+/// @brief Get the polyset types that is a superset of two types on the given
+/// cell
+/// @param[in] cell The cell type
+/// @param[in] type1 The first polyset type
+/// @param[in] type2 The second polyset type
+/// @return The superset type
+polyset::type superset(cell::type cell, polyset::type type1,
+                       polyset::type type2);
+
+/// @brief Get the polyset type that represents the restrictions of a type on a
+/// subentity
+/// @param[in] ptype The polyset type
+/// @param[in] cell The cell type
+/// @param[in] restriction_cell The cell type of the subentity
+/// @return The restricted polyset type
+polyset::type restriction(polyset::type ptype, cell::type cell,
+                          cell::type restriction_cell);
+
 } // namespace basix::polyset
diff --git a/cpp/basix/quadrature.cpp b/cpp/basix/quadrature.cpp
index a9e0f29a1..39c40643e 100644
--- a/cpp/basix/quadrature.cpp
+++ b/cpp/basix/quadrature.cpp
@@ -4706,7 +4706,38 @@ std::array<std::vector<T>, 2> make_xiao_gimbutas_quadrature(cell::type celltype,
     throw std::runtime_error(
         "Xiao-Gimbutas is only implemented for triangles.");
   }
-} // namespace
+}
+//-----------------------------------------------------------------------------
+template <std::floating_point T>
+std::array<std::vector<T>, 2>
+make_macroedge_quadrature(quadrature::type rule, cell::type celltype, int m)
+{
+  auto standard_q = quadrature::make_quadrature<T>(rule, celltype,
+                                                   polyset::type::standard, m);
+  switch (celltype)
+  {
+  case cell::type::interval:
+  {
+    if (m == 0)
+    {
+      return standard_q;
+    }
+    const std::size_t npts = standard_q[0].size();
+    std::vector<T> x(npts * 2);
+    std::vector<T> w(npts * 2);
+    for (std::size_t i = 0; i < npts; ++i)
+    {
+      x[i] = 0.5 * standard_q[0][i];
+      x[npts + i] = 0.5 + 0.5 * standard_q[0][i];
+      w[i] = 0.5 * standard_q[1][i];
+      w[npts + i] = 0.5 * standard_q[1][i];
+    }
+    return {std::move(x), std::move(w)};
+  }
+  default:
+    throw std::runtime_error("Macro quadrature not supported on this cell.");
+  }
+}
 //-----------------------------------------------------------------------------
 } // namespace
 //-----------------------------------------------------------------------------
@@ -4740,37 +4771,40 @@ quadrature::type quadrature::get_default_rule(cell::type celltype, int m)
 //-----------------------------------------------------------------------------
 template <std::floating_point T>
 std::array<std::vector<T>, 2>
-quadrature::make_quadrature(quadrature::type rule, cell::type celltype, int m)
+quadrature::make_quadrature(quadrature::type rule, cell::type celltype,
+                            polyset::type polytype, int m)
 {
-  switch (rule)
+  switch (polytype)
   {
-  case quadrature::type::Default:
-    return make_quadrature<T>(get_default_rule(celltype, m), celltype, m);
-  case quadrature::type::gauss_jacobi:
-    return make_gauss_jacobi_quadrature<T>(celltype, m);
-  case quadrature::type::gll:
-    return make_gll_quadrature<T>(celltype, m);
-  case quadrature::type::xiao_gimbutas:
-    return make_xiao_gimbutas_quadrature<T>(celltype, m);
-  case quadrature::type::zienkiewicz_taylor:
-    return make_zienkiewicz_taylor_quadrature<T>(celltype, m);
-  case quadrature::type::keast:
-    return make_keast_quadrature<T>(celltype, m);
-  case quadrature::type::strang_fix:
-    return make_strang_fix_quadrature<T>(celltype, m);
+  case polyset::type::standard:
+    switch (rule)
+    {
+    case quadrature::type::Default:
+      return make_quadrature<T>(get_default_rule(celltype, m), celltype,
+                                polytype, m);
+    case quadrature::type::gauss_jacobi:
+      return make_gauss_jacobi_quadrature<T>(celltype, m);
+    case quadrature::type::gll:
+      return make_gll_quadrature<T>(celltype, m);
+    case quadrature::type::xiao_gimbutas:
+      return make_xiao_gimbutas_quadrature<T>(celltype, m);
+    case quadrature::type::zienkiewicz_taylor:
+      return make_zienkiewicz_taylor_quadrature<T>(celltype, m);
+    case quadrature::type::keast:
+      return make_keast_quadrature<T>(celltype, m);
+    case quadrature::type::strang_fix:
+      return make_strang_fix_quadrature<T>(celltype, m);
+    default:
+      throw std::runtime_error("Unknown quadrature rule");
+    }
+  case polyset::type::macroedge:
+    return make_macroedge_quadrature<T>(rule, celltype, m);
   default:
-    throw std::runtime_error("Unknown quadrature rule");
+    throw std::runtime_error("Unsupported polyset type");
   }
 }
 //-----------------------------------------------------------------------------
 template <std::floating_point T>
-std::array<std::vector<T>, 2> quadrature::make_quadrature(cell::type celltype,
-                                                          int m)
-{
-  return make_quadrature<T>(quadrature::type::Default, celltype, m);
-}
-//-----------------------------------------------------------------------------
-template <std::floating_point T>
 std::vector<T> quadrature::get_gl_points(int m)
 {
   std::vector<T> pts = compute_gauss_jacobi_points<T>(0, m);
@@ -4787,14 +4821,9 @@ std::vector<T> quadrature::get_gll_points(int m)
 }
 //-----------------------------------------------------------------------------
 template std::array<std::vector<float>, 2>
-quadrature::make_quadrature(quadrature::type, cell::type, int);
-template std::array<std::vector<double>, 2>
-quadrature::make_quadrature(quadrature::type, cell::type, int);
-
-template std::array<std::vector<float>, 2>
-quadrature::make_quadrature(cell::type, int);
+quadrature::make_quadrature(quadrature::type, cell::type, polyset::type, int);
 template std::array<std::vector<double>, 2>
-quadrature::make_quadrature(cell::type, int);
+quadrature::make_quadrature(quadrature::type, cell::type, polyset::type, int);
 
 template std::vector<float> quadrature::get_gl_points(int);
 template std::vector<double> quadrature::get_gl_points(int);
diff --git a/cpp/basix/quadrature.h b/cpp/basix/quadrature.h
index 041bc5509..cdac684db 100644
--- a/cpp/basix/quadrature.h
+++ b/cpp/basix/quadrature.h
@@ -5,6 +5,7 @@
 #pragma once
 
 #include "cell.h"
+#include "polyset.h"
 #include <array>
 #include <concepts>
 #include <vector>
@@ -28,22 +29,15 @@ enum class type
 /// Make a quadrature rule on a reference cell
 /// @param[in] rule Type of quadrature rule (or use quadrature::Default)
 /// @param[in] celltype The cell type
+/// @param[in] polytype The polyset type
 /// @param[in] m Maximum degree of polynomial that this quadrature rule
 /// will integrate exactly
 /// @return List of points and list of weights. The number of points
 /// arrays has shape (num points, gdim)
 template <std::floating_point T>
 std::array<std::vector<T>, 2> make_quadrature(const quadrature::type rule,
-                                              cell::type celltype, int m);
-
-/// Make a default quadrature rule on reference cell
-/// @param[in] celltype The cell type
-/// @param[in] m Maximum degree of polynomial that this quadrature rule
-/// will integrate exactly
-/// @return List of points and list of weights. The number of points
-/// arrays has shape (num points, gdim)
-template <std::floating_point T>
-std::array<std::vector<T>, 2> make_quadrature(cell::type celltype, int m);
+                                              cell::type celltype,
+                                              polyset::type polytype, int m);
 
 /// Get the default quadrature type for the given cell and order
 /// @param[in] celltype The cell type
diff --git a/cpp/docs.template b/cpp/docs.template
index 3b2495c91..b110f118f 100644
--- a/cpp/docs.template
+++ b/cpp/docs.template
@@ -228,23 +228,24 @@ Returns:
 )";
 
 {{DOCTYPE}} create_custom_element = R"(
-{{finite-element.h > create_custom_element(cell_type, value_shape, wcoeffs, x, M, interpolation_nderivs, map_type, sobolev_space, discontinuous, highest_complete_degree, highest_degree) > doc}}
+{{finite-element.h > create_custom_element(cell_type, value_shape, wcoeffs, x, M, interpolation_nderivs, map_type, sobolev_space, discontinuous, highest_complete_degree, highest_degree, poly_type) > doc}}
 
 Args:
-    {{finite-element.h > create_custom_element(cell_type, value_shape, wcoeffs, x, M, interpolation_nderivs, map_type, sobolev_space, discontinuous, highest_complete_degree, highest_degree) > param > cell_type}}
-    {{finite-element.h > create_custom_element(cell_type, value_shape, wcoeffs, x, M, interpolation_nderivs, map_type, sobolev_space, discontinuous, highest_complete_degree, highest_degree) > param > value_shape}}
-    {{finite-element.h > create_custom_element(cell_type, value_shape, wcoeffs, x, M, interpolation_nderivs, map_type, sobolev_space, discontinuous, highest_complete_degree, highest_degree) > param > wcoeffs}}
-    {{finite-element.h > create_custom_element(cell_type, value_shape, wcoeffs, x, M, interpolation_nderivs, map_type, sobolev_space, discontinuous, highest_complete_degree, highest_degree) > param > x}}
-    {{finite-element.h > create_custom_element(cell_type, value_shape, wcoeffs, x, M, interpolation_nderivs, map_type, sobolev_space, discontinuous, highest_complete_degree, highest_degree) > param > M}}
-    {{finite-element.h > create_custom_element(cell_type, value_shape, wcoeffs, x, M, interpolation_nderivs, map_type, sobolev_space, discontinuous, highest_complete_degree, highest_degree) > param > interpolation_nderivs}}
-    {{finite-element.h > create_custom_element(cell_type, value_shape, wcoeffs, x, M, interpolation_nderivs, map_type, sobolev_space, discontinuous, highest_complete_degree, highest_degree) > param > map_type}}
-    {{finite-element.h > create_custom_element(cell_type, value_shape, wcoeffs, x, M, interpolation_nderivs, map_type, sobolev_space, discontinuous, highest_complete_degree, highest_degree) > param > sobolev_space}}
-    {{finite-element.h > create_custom_element(cell_type, value_shape, wcoeffs, x, M, interpolation_nderivs, map_type, sobolev_space, discontinuous, highest_complete_degree, highest_degree) > param > discontinuous}}
-    {{finite-element.h > create_custom_element(cell_type, value_shape, wcoeffs, x, M, interpolation_nderivs, map_type, sobolev_space, discontinuous, highest_complete_degree, highest_degree) > param > highest_complete_degree}}
-    {{finite-element.h > create_custom_element(cell_type, value_shape, wcoeffs, x, M, interpolation_nderivs, map_type, sobolev_space, discontinuous, highest_complete_degree, highest_degree) > param > highest_degree}}
+    {{finite-element.h > create_custom_element(cell_type, value_shape, wcoeffs, x, M, interpolation_nderivs, map_type, sobolev_space, discontinuous, highest_complete_degree, highest_degree, poly_type) > param > cell_type}}
+    {{finite-element.h > create_custom_element(cell_type, value_shape, wcoeffs, x, M, interpolation_nderivs, map_type, sobolev_space, discontinuous, highest_complete_degree, highest_degree, poly_type) > param > value_shape}}
+    {{finite-element.h > create_custom_element(cell_type, value_shape, wcoeffs, x, M, interpolation_nderivs, map_type, sobolev_space, discontinuous, highest_complete_degree, highest_degree, poly_type) > param > wcoeffs}}
+    {{finite-element.h > create_custom_element(cell_type, value_shape, wcoeffs, x, M, interpolation_nderivs, map_type, sobolev_space, discontinuous, highest_complete_degree, highest_degree, poly_type) > param > x}}
+    {{finite-element.h > create_custom_element(cell_type, value_shape, wcoeffs, x, M, interpolation_nderivs, map_type, sobolev_space, discontinuous, highest_complete_degree, highest_degree, poly_type) > param > M}}
+    {{finite-element.h > create_custom_element(cell_type, value_shape, wcoeffs, x, M, interpolation_nderivs, map_type, sobolev_space, discontinuous, highest_complete_degree, highest_degree, poly_type) > param > interpolation_nderivs}}
+    {{finite-element.h > create_custom_element(cell_type, value_shape, wcoeffs, x, M, interpolation_nderivs, map_type, sobolev_space, discontinuous, highest_complete_degree, highest_degree, poly_type) > param > map_type}}
+    {{finite-element.h > create_custom_element(cell_type, value_shape, wcoeffs, x, M, interpolation_nderivs, map_type, sobolev_space, discontinuous, highest_complete_degree, highest_degree, poly_type) > param > sobolev_space}}
+    {{finite-element.h > create_custom_element(cell_type, value_shape, wcoeffs, x, M, interpolation_nderivs, map_type, sobolev_space, discontinuous, highest_complete_degree, highest_degree, poly_type) > param > discontinuous}}
+    {{finite-element.h > create_custom_element(cell_type, value_shape, wcoeffs, x, M, interpolation_nderivs, map_type, sobolev_space, discontinuous, highest_complete_degree, highest_degree, poly_type) > param > highest_complete_degree}}
+    {{finite-element.h > create_custom_element(cell_type, value_shape, wcoeffs, x, M, interpolation_nderivs, map_type, sobolev_space, discontinuous, highest_complete_degree, highest_degree, poly_type) > param > highest_degree}}
+    {{finite-element.h > create_custom_element(cell_type, value_shape, wcoeffs, x, M, interpolation_nderivs, map_type, sobolev_space, discontinuous, highest_complete_degree, highest_degree, poly_type) > param > poly_type}}
 
 Returns:
-    {{finite-element.h > create_custom_element(cell_type, value_shape, wcoeffs, x, M, interpolation_nderivs, map_type, sobolev_space, discontinuous, highest_complete_degree, highest_degree) > return}}
+    {{finite-element.h > create_custom_element(cell_type, value_shape, wcoeffs, x, M, interpolation_nderivs, map_type, sobolev_space, discontinuous, highest_complete_degree, highest_degree, poly_type) > return}}
 )";
 
 
@@ -282,6 +283,7 @@ Returns:
 
 Args:
     {{polyset.h > tabulate > param > celltype}}
+    {{polyset.h > tabulate > param > ptype}}
     {{polyset.h > tabulate > param > d}}
     {{polyset.h > tabulate > param > n}}
     {{polyset.h > tabulate > param > x}}
@@ -315,27 +317,17 @@ Returns:
     {{polynomials.h > dim > return}}
 )";
 
-{{DOCTYPE}} make_quadrature__rule_celltype_m = R"(
-{{quadrature.h > make_quadrature(rule, celltype, m) > doc}}
+{{DOCTYPE}} make_quadrature__rule_celltype_polytype_m = R"(
+{{quadrature.h > make_quadrature(rule, celltype, polytype, m) > doc}}
 
 Args:
-    {{quadrature.h > make_quadrature(rule, celltype, m) > param > rule}}
-    {{quadrature.h > make_quadrature(rule, celltype, m) > param > celltype}}
-    {{quadrature.h > make_quadrature(rule, celltype, m) > param > m}}
+    {{quadrature.h > make_quadrature(rule, celltype, polytype, m) > param > rule}}
+    {{quadrature.h > make_quadrature(rule, celltype, polytype, m) > param > celltype}}
+    {{quadrature.h > make_quadrature(rule, celltype, polytype, m) > param > polytype}}
+    {{quadrature.h > make_quadrature(rule, celltype, polytype, m) > param > m}}
 
 Returns:
-    {{quadrature.h > make_quadrature(rule, celltype, m) > return}}
-)";
-
-{{DOCTYPE}} make_quadrature__celltype_m = R"(
-{{quadrature.h > make_quadrature(celltype, m) > doc}}
-
-Args:
-    {{quadrature.h > make_quadrature(celltype, m) > param > celltype}}
-    {{quadrature.h > make_quadrature(celltype, m) > param > m}}
-
-Returns:
-    {{quadrature.h > make_quadrature(rule, celltype, m) > return}}
+    {{quadrature.h > make_quadrature(rule, celltype, polytype, m) > return}}
 )";
 
 {{DOCTYPE}} index__p = R"(
@@ -382,4 +374,30 @@ Returns:
     {{sobolev-spaces.h > space_intersection(space1, space2) > return}}
 )";
 
+{{DOCTYPE}} superset = R"(
+{{polyset.h > superset > doc}}
+
+Args:
+    {{polyset.h > superset > param > cell}}
+    {{polyset.h > superset > param > type1}}
+    {{polyset.h > superset > param > type2}}
+
+Returns::
+    {{polyset.h > superset > return}}
+)";
+
+{{DOCTYPE}} restriction = R"(
+{{polyset.h > restriction > doc}}
+
+Args:
+    {{polyset.h > restriction > param > ptype}}
+    {{polyset.h > restriction > param > cell}}
+    {{polyset.h > restriction > param > restriction_cell}}
+
+Returns::
+    {{polyset.h > restriction > return}}
+)";
+
+
+
 } // namespace basix::docstring
diff --git a/cpp/generate_docs.py b/cpp/generate_docs.py
index 099f9827c..aa9c8cd4b 100644
--- a/cpp/generate_docs.py
+++ b/cpp/generate_docs.py
@@ -89,6 +89,7 @@ def get_docstring(matches):
     doc = doc.replace("@f$", "$")
     doc = doc.replace("@f[", "\\[")
     doc = doc.replace("@f]", "\\]")
+    doc = doc.replace("@brief", "")
     doc = doc.replace("@note", "NOTE:")
     doc = doc.replace("@todo", "TODO:")
     doc = doc.replace("@warning", "WARNING:")
diff --git a/demo/python/demo_custom_element.py b/demo/python/demo_custom_element.py
index fa7b4b98e..b4855d14f 100644
--- a/demo/python/demo_custom_element.py
+++ b/demo/python/demo_custom_element.py
@@ -9,7 +9,7 @@
 
 import basix
 import numpy as np
-from basix import CellType, MapType, PolynomialType, LatticeType, SobolevSpace
+from basix import CellType, MapType, PolynomialType, LatticeType, SobolevSpace, PolysetType
 
 # Lagrange element with bubble
 # ============================
@@ -140,9 +140,11 @@
 # - The highest degree of a polynomial in the element. In this example, this is 2. It is
 #   important that this value is correct, as it will be used to determine the number of
 #   polynomials to use when creating and tabulating the element.
+# - The type of polynomial set to use. In this example, we use standard polynomials.
 
 element = basix.create_custom_element(
-    CellType.quadrilateral, [], wcoeffs, x, M, 0, MapType.identity, SobolevSpace.H1, False, 1, 2)
+    CellType.quadrilateral, [], wcoeffs, x, M, 0, MapType.identity, SobolevSpace.H1, False, 1, 2,
+    PolysetType.standard)
 
 # We can now use this element in the same way we can use a built-in element. For example, we
 # can tabulate the element at a set of points. If the points we use are the same as the points
@@ -235,7 +237,7 @@
 
 element = basix.create_custom_element(
     CellType.triangle, [2], wcoeffs, x, M, 0, MapType.contravariantPiola, SobolevSpace.HDiv,
-    False, 0, 1)
+    False, 0, 1, PolysetType.standard)
 
 # To confirm that we have defined this element correctly, we compare it to the built-in
 # Raviart--Thomas element.
diff --git a/demo/python/demo_custom_element_conforming_cr.py b/demo/python/demo_custom_element_conforming_cr.py
index 4f310c3b1..f75cf896d 100644
--- a/demo/python/demo_custom_element_conforming_cr.py
+++ b/demo/python/demo_custom_element_conforming_cr.py
@@ -9,7 +9,7 @@
 
 import basix
 import numpy as np
-from basix import CellType, MapType, PolynomialType, LatticeType, SobolevSpace
+from basix import CellType, MapType, PolynomialType, LatticeType, SobolevSpace, PolysetType
 
 from mpl_toolkits import mplot3d  # noqa: F401
 import matplotlib.pyplot as plt
@@ -51,7 +51,8 @@ def create_ccr_triangle(degree):
         M[2].append(np.zeros((0, 1, 0, 1)))
 
         return basix.create_custom_element(
-            CellType.triangle, [], wcoeffs, x, M, 0, MapType.identity, SobolevSpace.L2, False, 1, 1)
+            CellType.triangle, [], wcoeffs, x, M, 0, MapType.identity, SobolevSpace.L2, False, 1, 1,
+            PolysetType.standard)
 
     npoly = (degree + 2) * (degree + 3) // 2
     ndofs = degree * (degree + 5) // 2
@@ -99,7 +100,7 @@ def create_ccr_triangle(degree):
 
     return basix.create_custom_element(
         CellType.triangle, [], wcoeffs, x, M, 0, MapType.identity, SobolevSpace.L2, False,
-        degree, degree + 1)
+        degree, degree + 1, PolysetType.standard)
 
 
 # We can then create a degree 2 conforming CR element.
diff --git a/demo/python/demo_quadrature.py b/demo/python/demo_quadrature.py
index de616251e..df63e1e72 100644
--- a/demo/python/demo_quadrature.py
+++ b/demo/python/demo_quadrature.py
@@ -35,7 +35,7 @@
 # to use a Gauss-Jacobi quadrature rule:
 
 points, weights = basix.make_quadrature(
-    basix.QuadratureType.gauss_jacobi, CellType.triangle, 4)
+    CellType.triangle, 4, rule=basix.QuadratureType.gauss_jacobi)
 
 # We now use this quadrature rule to integrate the functions :math:`f(x,y)=x^3y`
 # and :math:`g(x,y)=x^3y^2` over the triangle. The exact values of these integrals
diff --git a/python/basix/__init__.py b/python/basix/__init__.py
index 6874ce3dc..daa0d9f3b 100644
--- a/python/basix/__init__.py
+++ b/python/basix/__init__.py
@@ -7,6 +7,11 @@
 from ._basixcpp import __version__
 from . import cell, finite_element, lattice, polynomials, quadrature, sobolev_spaces, variants
 from ._basixcpp import (CellType, LatticeType, LatticeSimplexMethod, ElementFamily, LagrangeVariant,
-                        DPCVariant, QuadratureType, PolynomialType, MapType, SobolevSpace)
+                        DPCVariant, QuadratureType, PolynomialType, MapType, SobolevSpace,
+                        PolysetType)
 from ._basixcpp import (create_lattice, create_element, compute_interpolation_operator, topology,
-                        geometry, make_quadrature, index, tabulate_polynomials, create_custom_element)
+                        geometry, index, tabulate_polynomials, create_custom_element)
+
+from ._basixcpp import superset as polyset_superset
+from ._basixcpp import restriction as polyset_restriction
+from .quadrature import make_quadrature
diff --git a/python/basix/_basixcpp.pyi b/python/basix/_basixcpp.pyi
index b19f08f95..192cb1eea 100644
--- a/python/basix/_basixcpp.pyi
+++ b/python/basix/_basixcpp.pyi
@@ -62,6 +62,7 @@ class ElementFamily:
     DPC: ClassVar[ElementFamily] = ...
     HHJ: ClassVar[ElementFamily] = ...
     Hermite: ClassVar[ElementFamily] = ...
+    iso: ClassVar[ElementFamily] = ...
     N1E: ClassVar[ElementFamily] = ...
     N2E: ClassVar[ElementFamily] = ...
     P: ClassVar[ElementFamily] = ...
@@ -84,6 +85,24 @@ class ElementFamily:
     @property
     def value(self) -> int: ...
 
+class PolysetType:
+    __members__: ClassVar[dict] = ...  # read-only
+    __entries: ClassVar[dict] = ...
+    standard: ClassVar[PolysetType] = ...
+    macroedge: ClassVar[PolysetType] = ...
+    def __init__(self, value: int) -> None: ...
+    def __eq__(self, other: object) -> bool: ...
+    def __getstate__(self) -> int: ...
+    def __hash__(self) -> int: ...
+    def __index__(self) -> int: ...
+    def __int__(self) -> int: ...
+    def __ne__(self, other: object) -> bool: ...
+    def __setstate__(self, state: int) -> None: ...
+    @property
+    def name(self) -> str: ...
+    @property
+    def value(self) -> int: ...
+
 class FiniteElement:
     def __hash__(self, object, int): ClassVar[None] = ...
     def __init__(self, *args, **kwargs) -> None: ...
@@ -102,6 +121,8 @@ class FiniteElement:
     @property
     def cell_type(self) -> CellType: ...
     @property
+    def polyset_type(self) -> PolysetType: ...
+    @property
     def coefficient_matrix(self) -> npt.NDArray[numpy.float64]: ...
     @property
     def degree(self) -> int: ...
@@ -317,7 +338,7 @@ def cell_facet_reference_volumes(arg0: CellType) -> npt.NDArray[numpy.float64]:
 def cell_volume(arg0: CellType) -> float: ...
 def sobolev_space_intersection(arg0: SobolevSpace, arg1: SobolevSpace) -> SobolevSpace: ...
 def compute_interpolation_operator(arg0: FiniteElement, arg1: FiniteElement) -> npt.NDArray[numpy.float64]: ...
-def create_custom_element(arg0: CellType, arg1: List[int], arg2: npt.NDArray[numpy.float64], arg3: List[List[npt.NDArray[numpy.float64]]], arg4: List[List[npt.NDArray[numpy.float64]]], arg5: int, arg6: MapType, arg7: SobolevSpace, arg8: bool, arg9: int, arg10: int) -> FiniteElement: ...
+def create_custom_element(arg0: CellType, arg1: List[int], arg2: npt.NDArray[numpy.float64], arg3: List[List[npt.NDArray[numpy.float64]]], arg4: List[List[npt.NDArray[numpy.float64]]], arg5: int, arg6: MapType, arg7: SobolevSpace, arg8: bool, arg9: int, arg10: int, arg11: PolysetType) -> FiniteElement: ...
 @overload
 def create_element(arg0: ElementFamily, arg1: CellType, arg2: int, arg3: bool) -> FiniteElement: ...
 @overload
@@ -345,10 +366,7 @@ def index(arg0: int) -> int: ...
 def index(arg0: int, arg1: int) -> int: ...
 @overload
 def index(arg0: int, arg1: int, arg2: int) -> int: ...
-@overload
-def make_quadrature(arg0: QuadratureType, arg1: CellType, arg2: int) -> Tuple[npt.NDArray[numpy.float64],npt.NDArray[numpy.float64]]: ...
-@overload
-def make_quadrature(arg0: CellType, arg1: int) -> Tuple[npt.NDArray[numpy.float64],npt.NDArray[numpy.float64]]: ...
+def make_quadrature(arg0: QuadratureType, arg1: CellType, arg2: PolysetType, arg3: int) -> Tuple[npt.NDArray[numpy.float64],npt.NDArray[numpy.float64]]: ...
 def sub_entity_connectivity(arg0) -> List[List[List[List[int]]]]: ...
 def sub_entity_geometry(arg0, arg1: int, arg2: int) -> npt.NDArray[numpy.float64]: ...
 def tabulate_polynomial_set(arg0: CellType, arg1: int, arg2: int, arg3: npt.NDArray[numpy.float64]) -> npt.NDArray[numpy.float64]: ...
@@ -358,3 +376,5 @@ def polynomials_dim(arg0: PolynomialType, arg1: CellType, arg2: int) -> int: ...
 def topology(arg0) -> List[List[List[int]]]: ...
 @overload
 def topology(vertexindices) -> Any: ...
+def superset(arg0: CellType, arg1: PolysetType, arg2: PolysetType) -> PolysetType: ...
+def restriction(arg0: PolysetType, arg1: CellType, arg2: CellType) -> PolysetType: ...
diff --git a/python/basix/finite_element.py b/python/basix/finite_element.py
index 067fa0eb9..5eeedb64b 100644
--- a/python/basix/finite_element.py
+++ b/python/basix/finite_element.py
@@ -45,6 +45,7 @@ def string_to_family(family: str, cell: str) -> _EF:
     if cell == "interval":
         families.update({
             "DPC": _EF.P,
+            "iso": _EF.iso,
         })
     # Family names that are valid for tensor product cells
     if cell in ["interval", "quadrilateral", "hexahedron"]:
diff --git a/python/basix/quadrature.py b/python/basix/quadrature.py
index 0cf442c58..0dcab02b8 100644
--- a/python/basix/quadrature.py
+++ b/python/basix/quadrature.py
@@ -1,6 +1,13 @@
 """Functions to manipulate quadrature types."""
 
+import typing as _typing
+import numpy as _np
+import numpy.typing as _npt
+
 from ._basixcpp import QuadratureType as _QT
+from ._basixcpp import CellType as _CT
+from ._basixcpp import PolysetType as _PT
+from ._basixcpp import make_quadrature as _mq
 
 
 def string_to_type(rule: str) -> _QT:
@@ -37,3 +44,20 @@ def type_to_string(quadraturetype: _QT) -> str:
         The quadrature rule as a string.
     """
     return quadraturetype.name
+
+
+def make_quadrature(
+    cell: _CT, degree: int, rule: _QT = _QT.Default, polyset_type: _PT = _PT.standard
+) -> _typing.Tuple[_npt.NDArray[_np.float64], _npt.NDArray[_np.float64]]:
+    """Create a quadrature rule.
+
+    Args:
+        cell: The cell type
+        degree: The maximum polynomial degree that will be integrated exactly
+        rule: The quadrature rule
+        polyset_type: The type of polynomial that will be integrated exactly
+
+    Returns:
+        The quadrature points and weights
+    """
+    return _mq(rule, cell, polyset_type, degree)
diff --git a/python/basix/ufl.py b/python/basix/ufl.py
index 1f75fa710..2b0a0bd7e 100644
--- a/python/basix/ufl.py
+++ b/python/basix/ufl.py
@@ -270,6 +270,11 @@ def highest_degree(self) -> int:
         """The highest degree of the element."""
         raise NotImplementedError()
 
+    @property
+    @_abstractmethod
+    def polyset_type(self) -> _basix.PolysetType:
+        """The polyset type of the element."""
+
     @property
     def _wcoeffs(self) -> _npt.NDArray[_np.float64]:
         """The coefficients used to define the polynomial set."""
@@ -491,6 +496,10 @@ def highest_degree(self) -> int:
         """The highest degree of the element."""
         return self.element.highest_degree
 
+    @property
+    def polyset_type(self) -> _basix.PolysetType:
+        return self.element.polyset_type
+
     @property
     def _wcoeffs(self) -> _npt.NDArray[_np.float64]:
         """The coefficients used to define the polynomial set."""
@@ -663,6 +672,10 @@ def cell_type(self) -> _basix.CellType:
         """Basix cell type used to initialise the element."""
         return self.element.cell_type
 
+    @property
+    def polyset_type(self) -> _basix.PolysetType:
+        return self.element.polyset_type
+
     @property
     def discontinuous(self) -> bool:
         """True if the discontinuous version of the element is used."""
@@ -886,6 +899,13 @@ def interpolation_nderivs(self) -> int:
         """Number of derivatives needed when interpolating."""
         return max([e.interpolation_nderivs for e in self._sub_elements])
 
+    @property
+    def polyset_type(self) -> _basix.PolysetType:
+        pt = _basix.PolysetType.standard
+        for e in self._sub_elements:
+            pt = _basix.polyset_superset(self.cell_type, pt, e.polyset_type)
+        return pt
+
 
 class _BlockedElement(_ElementBase):
     """Element with a block size that contains multiple copies of a sub element.
@@ -1112,6 +1132,10 @@ def highest_degree(self) -> int:
         """The highest degree of the element."""
         return self.sub_element.highest_degree
 
+    @property
+    def polyset_type(self) -> _basix.PolysetType:
+        return self.sub_element.polyset_type
+
     @property
     def _wcoeffs(self) -> _npt.NDArray[_np.float64]:
         """Coefficients used to define the polynomial set."""
@@ -1332,6 +1356,7 @@ def enriched_element(elements: _typing.List[_ElementBase],
             set equal to the topological dimension of the cell.
     """
     ct = elements[0].cell_type
+    ptype = elements[0].polyset_type
     vshape = elements[0].value_shape()
     vsize = elements[0].value_size
     if map_type is None:
@@ -1349,6 +1374,8 @@ def enriched_element(elements: _typing.List[_ElementBase],
             discontinuous = False
         if e.cell_type != ct:
             raise ValueError("Enriched elements on different cell types not supported.")
+        if e.polyset_type != ptype:
+            raise ValueError("Enriched elements on different polyset types not supported.")
         if e.value_shape() != vshape or e.value_size != vsize:
             raise ValueError("Enriched elements on different value shapes not supported.")
     nderivs = max(e.interpolation_nderivs for e in elements)
@@ -1382,7 +1409,7 @@ def enriched_element(elements: _typing.List[_ElementBase],
         row += e.dim
 
     return custom_element(ct, list(vshape), wcoeffs, x, M, nderivs,
-                          map_type, ss, discontinuous, hcd, hd, gdim=gdim)
+                          map_type, ss, discontinuous, hcd, hd, ptype, gdim=gdim)
 
 
 def custom_element(cell_type: _basix.CellType, value_shape: _typing.Union[_typing.List[int], _typing.Tuple[int, ...]],
@@ -1390,6 +1417,7 @@ def custom_element(cell_type: _basix.CellType, value_shape: _typing.Union[_typin
                    M: _typing.List[_typing.List[_npt.NDArray[_np.float64]]], interpolation_nderivs: int,
                    map_type: _basix.MapType, sobolev_space: _basix.SobolevSpace, discontinuous: bool,
                    highest_complete_degree: int, highest_degree: int,
+                   polyset_type: _basix.PolysetType = _basix.PolysetType.standard,
                    gdim: _typing.Optional[int] = None) -> _ElementBase:
     """Create a UFL compatible custom Basix element.
 
@@ -1408,12 +1436,13 @@ def custom_element(cell_type: _basix.CellType, value_shape: _typing.Union[_typin
         highest_complete_degree: The highest degree n such that a Lagrange (or vector Lagrange)
         element of degree n is a subspace of this element
         highest_degree: The degree of a polynomial in this element's polyset
+        polyset_type: The polyset type for the element
         gdim: Geometric dimension. If not set the geometric dimension is
             set equal to the topological dimension of the cell.
     """
     return _BasixElement(_basix.create_custom_element(
         cell_type, list(value_shape), wcoeffs, x, M, interpolation_nderivs,
-        map_type, sobolev_space, discontinuous, highest_complete_degree, highest_degree), gdim=gdim)
+        map_type, sobolev_space, discontinuous, highest_complete_degree, highest_degree, polyset_type), gdim=gdim)
 
 
 def mixed_element(elements: _typing.List[_ElementBase], gdim: _typing.Optional[int] = None) -> _ElementBase:
diff --git a/python/wrapper.cpp b/python/wrapper.cpp
index eb204f0bc..de4bbd921 100644
--- a/python/wrapper.cpp
+++ b/python/wrapper.cpp
@@ -228,7 +228,8 @@ Interface to the Basix C++ library.
       .value("serendipity", element::family::serendipity)
       .value("DPC", element::family::DPC)
       .value("CR", element::family::CR)
-      .value("Hermite", element::family::Hermite);
+      .value("Hermite", element::family::Hermite)
+      .value("iso", element::family::iso);
 
   py::class_<FiniteElement<double>>(m, "FiniteElement", "Finite Element")
       .def(
@@ -351,6 +352,8 @@ Interface to the Basix C++ library.
       .def_property_readonly("highest_complete_degree",
                              &FiniteElement<double>::highest_complete_degree)
       .def_property_readonly("cell_type", &FiniteElement<double>::cell_type)
+      .def_property_readonly("polyset_type",
+                             &FiniteElement<double>::polyset_type)
       .def_property_readonly("dim", &FiniteElement<double>::dim)
       .def_property_readonly("num_entity_dofs",
                              [](const FiniteElement<double>& self)
@@ -535,8 +538,8 @@ Interface to the Basix C++ library.
              std::vector<py::array_t<double, py::array::c_style>>>& M,
          int interpolation_nderivs, maps::type map_type,
          sobolev::space sobolev_space, bool discontinuous,
-         int highest_complete_degree,
-         int highest_degree) -> FiniteElement<double>
+         int highest_complete_degree, int highest_degree,
+         polyset::type poly_type) -> FiniteElement<double>
       {
         if (x.size() != 4)
           throw std::runtime_error("x has the wrong size");
@@ -577,7 +580,7 @@ Interface to the Basix C++ library.
             mdspan_t<const double, 2>(wcoeffs.data(), wcoeffs.shape(0),
                                       wcoeffs.shape(1)),
             _x, _M, interpolation_nderivs, map_type, sobolev_space,
-            discontinuous, highest_complete_degree, highest_degree);
+            discontinuous, highest_complete_degree, highest_degree, poly_type);
       },
       basix::docstring::create_custom_element.c_str());
 
@@ -612,41 +615,52 @@ Interface to the Basix C++ library.
       },
       basix::docstring::compute_interpolation_operator.c_str());
 
+  py::enum_<polyset::type>(m, "PolysetType")
+      .value("standard", polyset::type::standard)
+      .value("macroedge", polyset::type::macroedge);
+
+  m.def(
+      "superset",
+      [](cell::type cell, polyset::type type1, polyset::type type2)
+      {
+        return polyset::superset(cell, type1, type2);
+      },
+      basix::docstring::superset.c_str());
+
+  m.def(
+      "restriction",
+      [](polyset::type ptype, cell::type cell, cell::type restriction_cell)
+      {
+        return polyset::restriction(ptype, cell, restriction_cell);
+      },
+      basix::docstring::restriction.c_str());
+
   m.def(
       "tabulate_polynomial_set",
-      [](cell::type celltype, int d, int n,
+      [](cell::type celltype, polyset::type polytype, int d, int n,
          const py::array_t<double, py::array::c_style>& x)
       {
         if (x.ndim() != 2)
           throw std::runtime_error("x has the wrong number of dimensions");
         stdex::mdspan<const double, stdex::dextents<std::size_t, 2>> _x(
             x.data(), x.shape(0), x.shape(1));
-        auto [p, shape] = polyset::tabulate(celltype, d, n, _x);
+        auto [p, shape] = polyset::tabulate(celltype, polytype, d, n, _x);
         return py::array_t<double>(shape, p.data());
       },
       basix::docstring::tabulate_polynomial_set.c_str());
 
   m.def(
       "make_quadrature",
-      [](quadrature::type rule, cell::type celltype, int m)
-      {
-        auto [pts, w] = quadrature::make_quadrature<double>(rule, celltype, m);
-        std::array<std::size_t, 2> shape = {w.size(), pts.size() / w.size()};
-        return std::pair(py::array_t<double>(shape, pts.data()),
-                         py::array_t<double>(w.size(), w.data()));
-      },
-      basix::docstring::make_quadrature__rule_celltype_m.c_str());
-
-  m.def(
-      "make_quadrature",
-      [](cell::type celltype, int m)
+      [](quadrature::type rule, cell::type celltype, polyset::type polytype,
+         int m)
       {
-        auto [pts, w] = quadrature::make_quadrature<double>(celltype, m);
+        auto [pts, w]
+            = quadrature::make_quadrature<double>(rule, celltype, polytype, m);
         std::array<std::size_t, 2> shape = {w.size(), pts.size() / w.size()};
         return std::pair(py::array_t<double>(shape, pts.data()),
                          py::array_t<double>(w.size(), w.data()));
       },
-      basix::docstring::make_quadrature__celltype_m.c_str());
+      basix::docstring::make_quadrature__rule_celltype_polytype_m.c_str());
 
   m.def("index", py::overload_cast<int>(&basix::indexing::idx),
         basix::docstring::index__p.c_str())
diff --git a/test/test_create.py b/test/test_create.py
index d4a3a58a4..9fdb2e744 100644
--- a/test/test_create.py
+++ b/test/test_create.py
@@ -29,6 +29,7 @@
     basix.ElementFamily.DPC,
     basix.ElementFamily.CR,
     basix.ElementFamily.Hermite,
+    basix.ElementFamily.iso,
     basix.ElementFamily.custom,
 ]
 
diff --git a/test/test_custom_element.py b/test/test_custom_element.py
index bad7e52f7..28baddbe2 100644
--- a/test/test_custom_element.py
+++ b/test/test_custom_element.py
@@ -23,7 +23,8 @@ def test_lagrange_custom_triangle_degree1():
 
     element = basix.create_custom_element(
         basix.CellType.triangle, [], wcoeffs,
-        x, M, 0, basix.MapType.identity, basix.SobolevSpace.H1, False, 1, 1)
+        x, M, 0, basix.MapType.identity, basix.SobolevSpace.H1, False, 1, 1,
+        basix.PolysetType.standard)
 
     points = basix.create_lattice(basix.CellType.triangle, 5, basix.LatticeType.equispaced, True)
     assert np.allclose(lagrange.tabulate(1, points), element.tabulate(1, points))
@@ -46,7 +47,8 @@ def test_lagrange_custom_triangle_degree1_l2piola():
 
     element = basix.create_custom_element(
         basix.CellType.triangle, [], wcoeffs,
-        x, M, 0, basix.MapType.L2Piola, basix.SobolevSpace.L2, False, 1, 1)
+        x, M, 0, basix.MapType.L2Piola, basix.SobolevSpace.L2, False, 1, 1,
+        basix.PolysetType.standard)
 
     points = basix.create_lattice(basix.CellType.triangle, 5, basix.LatticeType.equispaced, True)
     assert np.allclose(lagrange.tabulate(1, points), element.tabulate(1, points))
@@ -70,7 +72,8 @@ def test_lagrange_custom_triangle_degree4():
 
     element = basix.create_custom_element(
         basix.CellType.triangle, [], wcoeffs,
-        x, M, 0, basix.MapType.identity, basix.SobolevSpace.H1, False, 4, 4)
+        x, M, 0, basix.MapType.identity, basix.SobolevSpace.H1, False, 4, 4,
+        basix.PolysetType.standard)
 
     points = basix.create_lattice(basix.CellType.triangle, 5, basix.LatticeType.equispaced, True)
     assert np.allclose(lagrange.tabulate(1, points), element.tabulate(1, points))
@@ -93,7 +96,8 @@ def test_lagrange_custom_quadrilateral_degree1():
 
     element = basix.create_custom_element(
         basix.CellType.quadrilateral, [], wcoeffs,
-        x, M, 0, basix.MapType.identity, basix.SobolevSpace.H1, False, 1, 1)
+        x, M, 0, basix.MapType.identity, basix.SobolevSpace.H1, False, 1, 1,
+        basix.PolysetType.standard)
 
     points = basix.create_lattice(basix.CellType.quadrilateral, 5, basix.LatticeType.equispaced, True)
     assert np.allclose(lagrange.tabulate(1, points), element.tabulate(1, points))
@@ -135,7 +139,7 @@ def test_raviart_thomas_triangle_degree1():
 
     element = basix.create_custom_element(
         basix.CellType.triangle, [2], wcoeffs, x, M, 0, basix.MapType.contravariantPiola,
-        basix.SobolevSpace.HDiv, False, 0, 1)
+        basix.SobolevSpace.HDiv, False, 0, 1, basix.PolysetType.standard)
 
     rt = basix.create_element(
         basix.ElementFamily.RT, basix.CellType.triangle, 1)
@@ -162,7 +166,7 @@ def create_lagrange1_quad(cell_type=basix.CellType.quadrilateral, degree=1, wcoe
         value_shape = []
     basix.create_custom_element(
         cell_type, value_shape, wcoeffs, x, M, interpolation_nderivs, basix.MapType.identity,
-        basix.SobolevSpace.H1, discontinuous, 1, degree)
+        basix.SobolevSpace.H1, discontinuous, 1, degree, basix.PolysetType.standard)
 
 
 def test_create_lagrange1_quad():
diff --git a/test/test_equality.py b/test/test_equality.py
index 046c9b932..90c84c5f4 100644
--- a/test/test_equality.py
+++ b/test/test_equality.py
@@ -35,10 +35,10 @@ def test_custom_element_equality():
 
     p1_custom = basix.create_custom_element(
         basix.CellType.triangle, [], wcoeffs,
-        x, M, 0, basix.MapType.identity, basix.SobolevSpace.H1, False, 1, 1)
+        x, M, 0, basix.MapType.identity, basix.SobolevSpace.H1, False, 1, 1, basix.PolysetType.standard)
     p1_custom_again = basix.create_custom_element(
         basix.CellType.triangle, [], wcoeffs,
-        x, M, 0, basix.MapType.identity, basix.SobolevSpace.H1, False, 1, 1)
+        x, M, 0, basix.MapType.identity, basix.SobolevSpace.H1, False, 1, 1, basix.PolysetType.standard)
 
     wcoeffs = np.eye(3)
     z = np.zeros((0, 2))
@@ -48,7 +48,7 @@ def test_custom_element_equality():
 
     cr_custom = basix.create_custom_element(
         basix.CellType.triangle, [], wcoeffs,
-        x, M, 0, basix.MapType.identity, basix.SobolevSpace.L2, False, 1, 1)
+        x, M, 0, basix.MapType.identity, basix.SobolevSpace.L2, False, 1, 1, basix.PolysetType.standard)
 
     assert p1_custom == p1_custom_again
     assert p1_custom != cr_custom
diff --git a/test/test_iso.py b/test/test_iso.py
new file mode 100644
index 000000000..a3fea7114
--- /dev/null
+++ b/test/test_iso.py
@@ -0,0 +1,30 @@
+import numpy as np
+import basix
+import pytest
+
+
+@pytest.mark.parametrize("degree", range(1, 5))
+@pytest.mark.parametrize("cell", [basix.CellType.interval])
+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)
+
+    assert e.dim == e2.dim
+
+
+def test_iso_interval_1():
+    e = basix.create_element(basix.ElementFamily.iso, basix.CellType.interval, 1)
+    pts = np.array([[i/20] for i in range(21)])
+
+    values = e.tabulate(0, pts)
+    print(values.shape)
+
+    for n, p in enumerate(pts):
+        if p[0] <= 0.5:
+            assert np.isclose(values[0, n, 0, 0], 1 - 2 * p[0])
+            assert np.isclose(values[0, n, 1, 0], 0.0)
+            assert np.isclose(values[0, n, 2, 0], 2 * p[0])
+        else:
+            assert np.isclose(values[0, n, 0, 0], 0.0)
+            assert np.isclose(values[0, n, 1, 0], 2 * p[0] - 1)
+            assert np.isclose(values[0, n, 2, 0], 2 - 2 * p[0])
diff --git a/test/test_orthonormal_basis.py b/test/test_orthonormal_basis.py
index d19a920e9..548c3ffc7 100644
--- a/test/test_orthonormal_basis.py
+++ b/test/test_orthonormal_basis.py
@@ -16,8 +16,8 @@ def test_quad(order):
         for q, v in zip(Lpts, Lwts):
             Qpts.append([p[0], q[0]])
             Qwts.append(u * v)
-    basis = basix._basixcpp.tabulate_polynomial_set(basix.CellType.quadrilateral,
-                                                    order, 0, Qpts)[0]
+    basis = basix._basixcpp.tabulate_polynomial_set(
+        basix.CellType.quadrilateral, basix.PolysetType.standard, order, 0, Qpts)[0]
     ndofs = basis.shape[0]
 
     mat = np.zeros((ndofs, ndofs))
@@ -39,8 +39,8 @@ def test_pyramid(order):
                 sc = (1.0 - r[0])
                 Qpts.append([p[0] * sc, q[0] * sc, r[0]])
                 Qwts.append(u * v * sc * sc * w)
-    basis = basix._basixcpp.tabulate_polynomial_set(basix.CellType.pyramid,
-                                                    order, 0, Qpts)[0]
+    basis = basix._basixcpp.tabulate_polynomial_set(
+        basix.CellType.pyramid, basix.PolysetType.standard,  order, 0, Qpts)[0]
     ndofs = basis.shape[0]
 
     mat = np.zeros((ndofs, ndofs))
@@ -61,8 +61,8 @@ def test_hex(order):
             for r, w in zip(Lpts, Lwts):
                 Qpts.append([p[0], q[0], r[0]])
                 Qwts.append(u * v * w)
-    basis = basix._basixcpp.tabulate_polynomial_set(basix.CellType.hexahedron,
-                                                    order, 0, Qpts)[0]
+    basis = basix._basixcpp.tabulate_polynomial_set(
+        basix.CellType.hexahedron, basix.PolysetType.standard, order, 0, Qpts)[0]
     ndofs = basis.shape[0]
 
     mat = np.zeros((ndofs, ndofs))
@@ -83,8 +83,8 @@ def test_prism(order):
         for q, v in zip(Lpts, Lwts):
             Qpts.append([p[0], p[1], q[0]])
             Qwts.append(u * v)
-    basis = basix._basixcpp.tabulate_polynomial_set(basix.CellType.prism,
-                                                    order, 0, Qpts)[0]
+    basis = basix._basixcpp.tabulate_polynomial_set(
+        basix.CellType.prism, basix.PolysetType.standard, order, 0, Qpts)[0]
     ndofs = basis.shape[0]
 
     mat = np.zeros((ndofs, ndofs))
@@ -104,9 +104,28 @@ def test_prism(order):
     basix.CellType.prism,
 ])
 @pytest.mark.parametrize("order", [0, 1, 2, 3, 4])
-def test_cell(cell_type, order):
+def test_standard(cell_type, order):
     Qpts, Qwts = basix.make_quadrature(cell_type, 2*order + 1)
-    basis = basix._basixcpp.tabulate_polynomial_set(cell_type, order, 0, Qpts)[0]
+    basis = basix._basixcpp.tabulate_polynomial_set(
+        cell_type, basix.PolysetType.standard, order, 0, Qpts)[0]
+
+    ndofs = basis.shape[0]
+    mat = np.zeros((ndofs, ndofs))
+    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))
+
+
+@pytest.mark.parametrize("cell_type", [
+    basix.CellType.interval,
+])
+@pytest.mark.parametrize("order", [0, 1, 2, 3, 4])
+def test_macroedge(cell_type, order):
+    Qpts, Qwts = basix.make_quadrature(cell_type, 2*order + 1, polyset_type=basix.PolysetType.macroedge)
+    basis = basix._basixcpp.tabulate_polynomial_set(
+        cell_type, basix.PolysetType.macroedge, order, 0, Qpts)[0]
 
     ndofs = basis.shape[0]
     mat = np.zeros((ndofs, ndofs))
diff --git a/test/test_orthonormal_basis_symbolic.py b/test/test_orthonormal_basis_symbolic.py
index f5e58c94b..60b97bd1a 100644
--- a/test/test_orthonormal_basis_symbolic.py
+++ b/test/test_orthonormal_basis_symbolic.py
@@ -27,7 +27,8 @@ def test_symbolic_interval():
 
     cell = basix.CellType.interval
     pts0 = basix.create_lattice(cell, 10, basix.LatticeType.equispaced, True)
-    wtab = basix._basixcpp.tabulate_polynomial_set(cell, n, nderiv, pts0)
+    wtab = basix._basixcpp.tabulate_polynomial_set(
+        cell, basix.PolysetType.standard, n, nderiv, pts0)
 
     for k in range(nderiv + 1):
         wsym = np.zeros_like(wtab[k])
@@ -51,7 +52,8 @@ def test_symbolic_quad():
     m = (n + 1)**2
     cell = basix.CellType.quadrilateral
     pts0 = basix.create_lattice(cell, 2, basix.LatticeType.equispaced, True)
-    wtab = basix._basixcpp.tabulate_polynomial_set(cell, n, nderiv, pts0)
+    wtab = basix._basixcpp.tabulate_polynomial_set(
+        cell, basix.PolysetType.standard, n, nderiv, pts0)
 
     for kx in range(nderiv):
         for ky in range(0, nderiv - kx):
diff --git a/test/test_quadrature.py b/test/test_quadrature.py
index 3cb942ea5..f8f0f3186 100644
--- a/test/test_quadrature.py
+++ b/test/test_quadrature.py
@@ -23,7 +23,7 @@ def test_cell_quadrature(celltype, order):
 @pytest.mark.parametrize("m", range(7))
 @pytest.mark.parametrize("scheme", [basix.QuadratureType.Default, basix.QuadratureType.gll])
 def test_qorder_line(m, scheme):
-    Qpts, Qwts = basix.make_quadrature(scheme, basix.CellType.interval, m)
+    Qpts, Qwts = basix.make_quadrature(basix.CellType.interval, m, rule=scheme)
     x = sympy.Symbol('x')
     f = x**m
     q = sympy.integrate(f, (x, 0, (1)))
@@ -36,7 +36,7 @@ def test_qorder_line(m, scheme):
 @pytest.mark.parametrize("m", range(6))
 @pytest.mark.parametrize("scheme", [basix.QuadratureType.Default, basix.QuadratureType.gauss_jacobi])
 def test_qorder_tri(m, scheme):
-    Qpts, Qwts = basix.make_quadrature(scheme, basix.CellType.triangle, m)
+    Qpts, Qwts = basix.make_quadrature(basix.CellType.triangle, m, rule=scheme)
     x = sympy.Symbol('x')
     y = sympy.Symbol('y')
     f = x**m + y**m
@@ -50,7 +50,7 @@ def test_qorder_tri(m, scheme):
 @pytest.mark.parametrize("m", range(1, 20))
 @pytest.mark.parametrize("scheme", [basix.QuadratureType.xiao_gimbutas])
 def test_xiao_gimbutas_tri(m, scheme):
-    Qpts, Qwts = basix.make_quadrature(scheme, basix.CellType.triangle, m)
+    Qpts, Qwts = basix.make_quadrature(basix.CellType.triangle, m, rule=scheme)
     x = sympy.Symbol('x')
     y = sympy.Symbol('y')
     f = x**m + 2 * y**m
@@ -64,7 +64,7 @@ def test_xiao_gimbutas_tri(m, scheme):
 @pytest.mark.parametrize("m", range(1, 16))
 @pytest.mark.parametrize("scheme", [basix.QuadratureType.xiao_gimbutas])
 def test_xiao_gimbutas_tet(m, scheme):
-    Qpts, Qwts = basix.make_quadrature(scheme, basix.CellType.tetrahedron, m)
+    Qpts, Qwts = basix.make_quadrature(basix.CellType.tetrahedron, m, rule=scheme)
     x = sympy.Symbol('x')
     y = sympy.Symbol('y')
     z = sympy.Symbol('z')
@@ -79,7 +79,7 @@ def test_xiao_gimbutas_tet(m, scheme):
 @pytest.mark.parametrize("m", range(9))
 @pytest.mark.parametrize("scheme", [basix.QuadratureType.Default, basix.QuadratureType.gauss_jacobi])
 def test_qorder_tet(m, scheme):
-    Qpts, Qwts = basix.make_quadrature(scheme, basix.CellType.tetrahedron, m)
+    Qpts, Qwts = basix.make_quadrature(basix.CellType.tetrahedron, m, rule=scheme)
     x = sympy.Symbol('x')
     y = sympy.Symbol('y')
     z = sympy.Symbol('z')
@@ -109,7 +109,8 @@ def test_gll():
     m = 5
 
     # 1D interval
-    pts, wts = basix.make_quadrature(basix.QuadratureType.gll, basix.CellType.interval, m+1)
+    pts, wts = basix.make_quadrature(
+        basix.CellType.interval, m+1, rule=basix.QuadratureType.gll)
     pts, wts = 2*pts.flatten()-1, 2*wts.flatten()
     ref_pts = np.array([-1., 1., -np.sqrt(3/7), 0.0, np.sqrt(3/7)])
     assert (np.allclose(pts.flatten(), ref_pts))
@@ -119,7 +120,8 @@ def test_gll():
     assert np.isclose(sum(wts), 2)
 
     # 2D quad
-    pts, wts = basix.make_quadrature(basix.QuadratureType.gll, basix.CellType.quadrilateral, m+1)
+    pts, wts = basix.make_quadrature(
+        basix.CellType.quadrilateral, m+1, rule=basix.QuadratureType.gll)
     pts, wts = 2*pts-1, 4*wts
     ref_pts2 = np.array([[x, y] for x in ref_pts for y in ref_pts])
     assert (np.allclose(pts, ref_pts2))
@@ -129,7 +131,8 @@ def test_gll():
     assert np.isclose(sum(wts), 4)
 
     # 3D hex
-    pts, wts = basix.make_quadrature(basix.QuadratureType.gll, basix.CellType.hexahedron, m+1)
+    pts, wts = basix.make_quadrature(
+        basix.CellType.hexahedron, m+1, rule=basix.QuadratureType.gll)
     pts, wts = 2*pts-1, 8*wts
     ref_pts3 = np.array([[x, y, z] for x in ref_pts for y in ref_pts for z in ref_pts])
     assert (np.allclose(pts, ref_pts3))
diff --git a/test/test_tensor_products.py b/test/test_tensor_products.py
index 6bccd1945..915ef2229 100644
--- a/test/test_tensor_products.py
+++ b/test/test_tensor_products.py
@@ -70,7 +70,7 @@ def test_tensor_product_factorisation_quadrilateral(degree):
 
     # Quadrature degree
     Q = 2 * P + 2
-    points, w = basix.make_quadrature(basix.QuadratureType.Default, cell_type, Q)
+    points, w = basix.make_quadrature(cell_type, Q)
 
     # FIXME: This test assumes all factors formed by a single element
     perm = factors[1]
@@ -88,7 +88,7 @@ def test_tensor_product_factorisation_quadrilateral(degree):
     assert points.shape[0] == (P+2) * (P+2)
 
     cell1d = element0.cell_type
-    points, _ = basix.make_quadrature(basix.QuadratureType.Default, cell1d, Q)
+    points, _ = basix.make_quadrature(cell1d, Q)
     data = element0.tabulate(1, points)
     phi0 = data[0, :, :, 0]
     dphi0 = data[1, :, :, 0]
@@ -129,7 +129,8 @@ def test_tensor_product_factorisation_hexahedron(degree):
 
     # Quadrature degree
     Q = 2 * P + 2
-    points, _ = basix.make_quadrature(basix.QuadratureType.Default, basix.CellType.hexahedron, Q)
+    points, _ = basix.make_quadrature(
+        basix.CellType.hexahedron, Q)
 
     # FIXME: This test assumes all factors formed by a single element
     perm = factors[1]
@@ -148,7 +149,7 @@ def test_tensor_product_factorisation_hexahedron(degree):
     assert points.shape[0] == (P+2) * (P+2) * (P+2)
 
     cell1d = element0.cell_type
-    points, w = basix.make_quadrature(basix.QuadratureType.Default, cell1d, Q)
+    points, w = basix.make_quadrature(cell1d, Q)
     data = element0.tabulate(1, points)
     phi0 = data[0, :, :, 0]
     dphi0 = data[1, :, :, 0]