diff --git a/demo/python/demo_custom_element.py b/demo/python/demo_custom_element.py
index af869a0c0..b4855d14f 100644
--- a/demo/python/demo_custom_element.py
+++ b/demo/python/demo_custom_element.py
@@ -67,7 +67,7 @@
 # We compute these integrals using a degree 4 quadrature rule (this is the largest degree
 # that the integrand will be, so these integrals will be exact).
 
-pts, wts = basix.quadrature.make_quadrature(CellType.quadrilateral, 4)
+pts, wts = basix.make_quadrature(CellType.quadrilateral, 4)
 poly = basix.tabulate_polynomials(PolynomialType.legendre, CellType.quadrilateral, 2, pts)
 x = pts[:, 0]
 y = pts[:, 1]
@@ -190,7 +190,7 @@
 wcoeffs[0, 0] = 1
 wcoeffs[1, 3] = 1
 
-pts, wts = basix.quadrature.make_quadrature(CellType.triangle, 2)
+pts, wts = basix.make_quadrature(CellType.triangle, 2)
 poly = basix.tabulate_polynomials(PolynomialType.legendre, CellType.triangle, 1, pts)
 x = pts[:, 0]
 y = pts[:, 1]
@@ -205,7 +205,7 @@
 # the element are integrals. We begin by defining a degree 1 quadrature rule on an interval.
 # This quadrature rule will be used to integrate on the edges of the triangle.
 
-pts, wts = basix.quadrature.make_quadrature(CellType.interval, 1)
+pts, wts = basix.make_quadrature(CellType.interval, 1)
 
 # The points associated with each edge are calculated by mapping the quadrature points to each edge.
 
diff --git a/demo/python/demo_custom_element_conforming_cr.py b/demo/python/demo_custom_element_conforming_cr.py
index 335dbec66..f75cf896d 100644
--- a/demo/python/demo_custom_element_conforming_cr.py
+++ b/demo/python/demo_custom_element_conforming_cr.py
@@ -63,7 +63,7 @@ def create_ccr_triangle(degree):
         wcoeffs[dof_n, dof_n] = 1
         dof_n += 1
 
-    pts, wts = basix.quadrature.make_quadrature(CellType.triangle, 2 * (degree + 1))
+    pts, wts = basix.make_quadrature(CellType.triangle, 2 * (degree + 1))
     poly = basix.tabulate_polynomials(PolynomialType.legendre, CellType.triangle, degree + 1, pts)
     for i in range(1, degree):
         x = pts[:, 0]
diff --git a/demo/python/demo_facet_integral.py b/demo/python/demo_facet_integral.py
index 70a279f09..883626317 100644
--- a/demo/python/demo_facet_integral.py
+++ b/demo/python/demo_facet_integral.py
@@ -26,7 +26,7 @@
 # rule on a triangle. We use an order 3 rule so that we can integrate the
 # basis functions in our space exactly.
 
-points, weights = basix.quadrature.make_quadrature(CellType.triangle, 3)
+points, weights = basix.make_quadrature(CellType.triangle, 3)
 
 # Next, we must map the quadrature points to our facet. We use the function
 # `geometry` to get the coordinates of the vertices of the tetrahedron, and
diff --git a/demo/python/demo_quadrature.py b/demo/python/demo_quadrature.py
index a24fd4d35..df63e1e72 100644
--- a/demo/python/demo_quadrature.py
+++ b/demo/python/demo_quadrature.py
@@ -28,13 +28,13 @@
 # `make_quadrature` returns two values: the points and the weights of the
 # quadrature rule.
 
-points, weights = basix.quadrature.make_quadrature(CellType.triangle, 4)
+points, weights = basix.make_quadrature(CellType.triangle, 4)
 
 # If we want to control the type of quadrature used, we can pass in three
 # inputs to `make_quadrautre`. For example, the following code would force basix
 # to use a Gauss-Jacobi quadrature rule:
 
-points, weights = basix.quadrature.make_quadrature(
+points, weights = basix.make_quadrature(
     CellType.triangle, 4, rule=basix.QuadratureType.gauss_jacobi)
 
 # We now use this quadrature rule to integrate the functions :math:`f(x,y)=x^3y`
diff --git a/python/basix/__init__.py b/python/basix/__init__.py
index 05176dc17..daa0d9f3b 100644
--- a/python/basix/__init__.py
+++ b/python/basix/__init__.py
@@ -14,3 +14,4 @@
 
 from ._basixcpp import superset as polyset_superset
 from ._basixcpp import restriction as polyset_restriction
+from .quadrature import make_quadrature
diff --git a/test/test_custom_element.py b/test/test_custom_element.py
index 19980a4f2..28baddbe2 100644
--- a/test/test_custom_element.py
+++ b/test/test_custom_element.py
@@ -111,13 +111,13 @@ def test_raviart_thomas_triangle_degree1():
     wcoeffs[0, 0] = 1
     wcoeffs[1, 3] = 1
 
-    pts, wts = basix.quadrature.make_quadrature(basix.CellType.triangle, 2)
+    pts, wts = basix.make_quadrature(basix.CellType.triangle, 2)
     poly = basix.tabulate_polynomials(basix.PolynomialType.legendre, basix.CellType.triangle, 1, pts)
     for i in range(3):
         wcoeffs[2, i] = sum(pts[:, 0] * poly[i, :] * wts)
         wcoeffs[2, 3 + i] = sum(pts[:, 1] * poly[i, :] * wts)
 
-    pts, wts = basix.quadrature.make_quadrature(basix.CellType.interval, 2)
+    pts, wts = basix.make_quadrature(basix.CellType.interval, 2)
 
     x = [[], [], [], []]
     for _ in range(3):
diff --git a/test/test_lagrange.py b/test/test_lagrange.py
index 28557b601..2ddbc53f2 100644
--- a/test/test_lagrange.py
+++ b/test/test_lagrange.py
@@ -583,7 +583,7 @@ def test_legendre_lagrange_variant(celltype, degree):
         basix.ElementFamily.P, celltype, degree, basix.LagrangeVariant.legendre, discontinuous=True)
 
     # Test that the basis functions are orthogonal
-    pts, wts = basix.quadrature.make_quadrature(celltype, degree * 2)
+    pts, wts = basix.make_quadrature(celltype, degree * 2)
     values = e.tabulate(0, pts)[0, :, :, 0].T
     for i, row_i in enumerate(values):
         for j, row_j in enumerate(values):
@@ -655,7 +655,7 @@ def test_legendre_dpc_variant(celltype, degree):
         basix.ElementFamily.DPC, celltype, degree, dpc_variant=basix.DPCVariant.legendre, discontinuous=True)
 
     # Test that the basis functions are orthogonal
-    pts, wts = basix.quadrature.make_quadrature(celltype, degree * 2)
+    pts, wts = basix.make_quadrature(celltype, degree * 2)
     values = e.tabulate(0, pts)[0, :, :, 0].T
     for i, row_i in enumerate(values):
         for j, row_j in enumerate(values):
diff --git a/test/test_orthonormal_basis.py b/test/test_orthonormal_basis.py
index 24e532047..548c3ffc7 100644
--- a/test/test_orthonormal_basis.py
+++ b/test/test_orthonormal_basis.py
@@ -9,7 +9,7 @@
 
 @pytest.mark.parametrize("order", [1, 2, 3])
 def test_quad(order):
-    Lpts, Lwts = basix.quadrature.make_quadrature(basix.CellType.interval, 2*order + 1)
+    Lpts, Lwts = basix.make_quadrature(basix.CellType.interval, 2*order + 1)
     Qwts = []
     Qpts = []
     for p, u in zip(Lpts, Lwts):
@@ -30,7 +30,7 @@ def test_quad(order):
 
 @pytest.mark.parametrize("order", [1, 2, 3, 4])
 def test_pyramid(order):
-    Lpts, Lwts = basix.quadrature.make_quadrature(basix.CellType.interval, 4*order + 2)
+    Lpts, Lwts = basix.make_quadrature(basix.CellType.interval, 4*order + 2)
     Qwts = []
     Qpts = []
     for p, u in zip(Lpts, Lwts):
@@ -53,7 +53,7 @@ def test_pyramid(order):
 
 @pytest.mark.parametrize("order", [1, 2, 3])
 def test_hex(order):
-    Lpts, Lwts = basix.quadrature.make_quadrature(basix.CellType.interval, 2*order + 1)
+    Lpts, Lwts = basix.make_quadrature(basix.CellType.interval, 2*order + 1)
     Qwts = []
     Qpts = []
     for p, u in zip(Lpts, Lwts):
@@ -75,8 +75,8 @@ def test_hex(order):
 
 @pytest.mark.parametrize("order", [1, 2, 3])
 def test_prism(order):
-    Tpts, Twts = basix.quadrature.make_quadrature(basix.CellType.triangle, 2*order + 1)
-    Lpts, Lwts = basix.quadrature.make_quadrature(basix.CellType.interval, 2*order + 1)
+    Tpts, Twts = basix.make_quadrature(basix.CellType.triangle, 2*order + 1)
+    Lpts, Lwts = basix.make_quadrature(basix.CellType.interval, 2*order + 1)
     Qwts = []
     Qpts = []
     for p, u in zip(Tpts, Twts):
@@ -105,7 +105,7 @@ def test_prism(order):
 ])
 @pytest.mark.parametrize("order", [0, 1, 2, 3, 4])
 def test_standard(cell_type, order):
-    Qpts, Qwts = basix.quadrature.make_quadrature(cell_type, 2*order + 1)
+    Qpts, Qwts = basix.make_quadrature(cell_type, 2*order + 1)
     basis = basix._basixcpp.tabulate_polynomial_set(
         cell_type, basix.PolysetType.standard, order, 0, Qpts)[0]
 
@@ -123,7 +123,7 @@ def test_standard(cell_type, order):
 ])
 @pytest.mark.parametrize("order", [0, 1, 2, 3, 4])
 def test_macroedge(cell_type, order):
-    Qpts, Qwts = basix.quadrature.make_quadrature(cell_type, 2*order + 1, polyset_type=basix.PolysetType.macroedge)
+    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]
 
diff --git a/test/test_polynomials.py b/test/test_polynomials.py
index 5db88ed77..0952c72b1 100644
--- a/test/test_polynomials.py
+++ b/test/test_polynomials.py
@@ -19,7 +19,7 @@
     basix.CellType.tetrahedron, basix.CellType.hexahedron,
 ])
 def test_legendre(cell_type, degree):
-    points, weights = basix.quadrature.make_quadrature(cell_type, 2 * degree)
+    points, weights = basix.make_quadrature(cell_type, 2 * degree)
 
     polys = basix.tabulate_polynomials(basix.PolynomialType.legendre, cell_type, degree, points)
 
@@ -63,7 +63,7 @@ def evaluate(function, pt):
                             x**2, x**2 * z, x**2 * z**2], 2],
 ])
 def test_order(cell_type, functions, degree):
-    points, weights = basix.quadrature.make_quadrature(cell_type, 2 * degree)
+    points, weights = basix.make_quadrature(cell_type, 2 * degree)
     polys = basix.tabulate_polynomials(basix.PolynomialType.legendre, cell_type, degree, points)
 
     assert len(functions) == polys.shape[0]
diff --git a/test/test_quadrature.py b/test/test_quadrature.py
index 7f93118fa..f8f0f3186 100644
--- a/test/test_quadrature.py
+++ b/test/test_quadrature.py
@@ -16,14 +16,14 @@
                                       (basix.CellType.tetrahedron, 1.0/6.0)])
 @pytest.mark.parametrize("order", range(9))
 def test_cell_quadrature(celltype, order):
-    Qpts, Qwts = basix.quadrature.make_quadrature(celltype[0], order)
+    Qpts, Qwts = basix.make_quadrature(celltype[0], order)
     assert np.isclose(sum(Qwts), celltype[1])
 
 
 @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.quadrature.make_quadrature(basix.CellType.interval, m, rule=scheme)
+    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.quadrature.make_quadrature(basix.CellType.triangle, m, rule=scheme)
+    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.quadrature.make_quadrature(basix.CellType.triangle, m, rule=scheme)
+    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.quadrature.make_quadrature(basix.CellType.tetrahedron, m, rule=scheme)
+    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.quadrature.make_quadrature(basix.CellType.tetrahedron, m, rule=scheme)
+    Qpts, Qwts = basix.make_quadrature(basix.CellType.tetrahedron, m, rule=scheme)
     x = sympy.Symbol('x')
     y = sympy.Symbol('y')
     z = sympy.Symbol('z')
@@ -92,7 +92,7 @@ def test_qorder_tet(m, scheme):
 
 
 def test_quadrature_function():
-    Qpts, Qwts = basix.quadrature.make_quadrature(basix.CellType.interval, 3)
+    Qpts, Qwts = basix.make_quadrature(basix.CellType.interval, 3)
     # Scale to interval [0.0, 2.0]
     Qpts *= 2.0
     Qwts *= 2.0
@@ -109,7 +109,7 @@ def test_gll():
     m = 5
 
     # 1D interval
-    pts, wts = basix.quadrature.make_quadrature(
+    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)])
@@ -120,7 +120,7 @@ def test_gll():
     assert np.isclose(sum(wts), 2)
 
     # 2D quad
-    pts, wts = basix.quadrature.make_quadrature(
+    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])
@@ -131,7 +131,7 @@ def test_gll():
     assert np.isclose(sum(wts), 4)
 
     # 3D hex
-    pts, wts = basix.quadrature.make_quadrature(
+    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])
diff --git a/test/test_tensor_products.py b/test/test_tensor_products.py
index f4496e800..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.quadrature.make_quadrature(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.quadrature.make_quadrature(cell1d, Q)
+    points, _ = basix.make_quadrature(cell1d, Q)
     data = element0.tabulate(1, points)
     phi0 = data[0, :, :, 0]
     dphi0 = data[1, :, :, 0]
@@ -129,7 +129,7 @@ def test_tensor_product_factorisation_hexahedron(degree):
 
     # Quadrature degree
     Q = 2 * P + 2
-    points, _ = basix.quadrature.make_quadrature(
+    points, _ = basix.make_quadrature(
         basix.CellType.hexahedron, Q)
 
     # FIXME: This test assumes all factors formed by a single element
@@ -149,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.quadrature.make_quadrature(cell1d, Q)
+    points, w = basix.make_quadrature(cell1d, Q)
     data = element0.tabulate(1, points)
     phi0 = data[0, :, :, 0]
     dphi0 = data[1, :, :, 0]