diff --git a/basix/_basixcpp.pyi b/basix/_basixcpp.pyi
index 445c2da51..04e8e023e 100644
--- a/basix/_basixcpp.pyi
+++ b/basix/_basixcpp.pyi
@@ -520,4 +520,6 @@ def topology(arg: CellType, /) -> list[list[list[int]]]: ...
 
 def tp_dof_ordering(arg0: ElementFamily, arg1: CellType, arg2: int, arg3: LagrangeVariant, arg4: DPCVariant, arg5: bool, /) -> list[int]: ...
 
+def lex_dof_ordering(arg0: ElementFamily, arg1: CellType, arg2: int, arg3: LagrangeVariant, arg4: DPCVariant, arg5: bool, /) -> list[int]: ...
+
 def tp_factors(arg0: ElementFamily, arg1: CellType, arg2: int, arg3: LagrangeVariant, arg4: DPCVariant, arg5: bool, arg6: Sequence[int], arg7: str, /) -> list[list[FiniteElement_float32]] | list[list[FiniteElement_float64]]: ...
diff --git a/cpp/basix/finite-element.cpp b/cpp/basix/finite-element.cpp
index 96ab41370..d97d9dbe1 100644
--- a/cpp/basix/finite-element.cpp
+++ b/cpp/basix/finite-element.cpp
@@ -509,6 +509,308 @@ std::vector<int> basix::tp_dof_ordering(element::family family, cell::type cell,
   return dof_ordering;
 }
 //-----------------------------------------------------------------------------
+std::vector<int> basix::lex_dof_ordering(element::family family, cell::type cell,
+                                         int degree, element::lagrange_variant,
+                                         element::dpc_variant, bool)
+{
+  std::vector<int> dof_ordering;
+  std::vector<int> perm;
+
+  switch (family)
+  {
+  case element::family::P:
+  {
+    switch (cell)
+    {
+    case cell::type::interval:
+    {
+      perm.push_back(0);
+      if (degree > 0)
+      {
+        for (int i = 2; i <= degree; ++i)
+          perm.push_back(i);
+        perm.push_back(1);
+      }
+      break;
+    }
+    case cell::type::quadrilateral:
+    {
+      perm.push_back(0);
+      if (degree > 0)
+      {
+        int n = degree - 1;
+        for (int i = 0; i < n; ++i)
+          perm.push_back(4 + i);
+        perm.push_back(1);
+        for (int i = 0; i < n; ++i)
+        {
+          perm.push_back(4 + n + i);
+          for (int j = 0; j < n; ++j)
+            perm.push_back(4 + j + (4 + i) * n);
+          perm.push_back(4 + 2 * n + i);
+        }
+        perm.push_back(2);
+        for (int i = 0; i < n; ++i)
+          perm.push_back(4 + 3 * n + i);
+        perm.push_back(3);
+      }
+      assert((int)perm.size() == (degree + 1) * (degree + 1));
+      break;
+    }
+    case cell::type::triangle:
+    {
+      perm.push_back(0);
+      if (degree > 0)
+      {
+        int n = degree - 1;
+        for (int i = 0; i < n; ++i)
+          perm.push_back(3 + 2 * n + i);
+        perm.push_back(1);
+        int dof = 3 + 3 * n;
+        for (int i = 0; i < n; ++i)
+        {
+          perm.push_back(3 + n + i);
+          for (int j = 0; j < n - 1 - i; ++j)
+            perm.push_back(dof++);
+          perm.push_back(3 + i);
+        }
+        perm.push_back(2);
+      }
+
+      assert((int)perm.size() == (degree + 1) * (degree + 2) / 2);
+      break;
+    }
+    case cell::type::hexahedron:
+    {
+      perm.push_back(0);
+      if (degree > 0)
+      {
+        int n = degree - 1;
+        for (int i = 0; i < n; ++i)
+          perm.push_back(8 + i);
+        perm.push_back(1);
+        for (int i = 0; i < n; ++i)
+        {
+          perm.push_back(8 + n + i);
+          for (int j = 0; j < n; ++j)
+            perm.push_back(8 + 12 * n + n * i + j);
+          perm.push_back(8 + 3 * n + i);
+        }
+        perm.push_back(2);
+        for (int i = 0; i < n; ++i)
+          perm.push_back(8 + 5 * n + i);
+        perm.push_back(3);
+
+        for (int i = 0; i < n; ++i)
+        {
+          perm.push_back(8 + 2 * n + i);
+          for (int j = 0; j < n; ++j)
+            perm.push_back(8 + 12 * n + n * n + n * i + j);
+          perm.push_back(8 + 4 * n + i);
+          for (int j = 0; j < n; ++j)
+          {
+            perm.push_back(8 + 12 * n + 2 * n * n + n * i + j);
+            for (int k = 0; k < n; ++k)
+              perm.push_back(8 + 12 * n + 6 * n * n + i * n * n + j * n + k);
+            perm.push_back(8 + 12 * n + 3 * n * n + n * i + j);
+          }
+          perm.push_back(8 + 6 * n + i);
+          for (int j = 0; j < n; ++j)
+            perm.push_back(8 + 12 * n + 4 * n * n + n * i + j);
+          perm.push_back(8 + 7 * n + i);
+        }
+        perm.push_back(4);
+        for (int i = 0; i < n; ++i)
+          perm.push_back(8 + 8 * n + i);
+        perm.push_back(5);
+        for (int i = 0; i < n; ++i)
+        {
+          perm.push_back(8 + 9 * n + i);
+          for (int j = 0; j < n; ++j)
+            perm.push_back(8 + 12 * n + 5 * n * n + n * i + j);
+          perm.push_back(8 + 10 * n + i);
+        }
+        perm.push_back(6);
+        for (int i = 0; i < n; ++i)
+          perm.push_back(8 + 11 * n + i);
+        perm.push_back(7);
+      }
+
+      assert((int)perm.size() == (degree + 1) * (degree + 1) * (degree + 1));
+      break;
+    }
+    case cell::type::tetrahedron:
+    {
+      perm.push_back(0);
+      if (degree > 0)
+      {
+        int n = degree - 1;
+        int face0 = 4 + 6 * n;
+        int face1 = 4 + 6 * n + n * (n - 1) / 2;
+        int face2 = 4 + 6 * n + n * (n - 1);
+        int face3 = 4 + 6 * n + n * (n - 1) * 3 / 2;
+        int interior = 4 + 6 * n + n * (n - 1) * 2;
+        for (int i = 0; i < n; ++i)
+          perm.push_back(4 + 5 * n + i);
+        perm.push_back(1);
+        for (int i = 0; i < n; ++i)
+        {
+          perm.push_back(4 + 4 * n + i);
+          for (int j = 0; j < n - 1 - i; ++j)
+            perm.push_back(face3++);
+          perm.push_back(4 + 2 * n + i);
+        }
+        perm.push_back(2);
+        for (int i = 0; i < n; ++i)
+        {
+          perm.push_back(4 + 3 * n + i);
+          for (int j = 0; j < n - 1 - i; ++j)
+            perm.push_back(face2++);
+          perm.push_back(4 + n + i);
+          for (int j = 0; j < n - 1 - i; ++j)
+          {
+            perm.push_back(face1++);
+            for (int k = 0; k < n - 2 - i - j; ++k)
+              perm.push_back(interior++);
+            perm.push_back(face0++);
+          }
+          perm.push_back(4 + i);
+        }
+        perm.push_back(3);
+      }
+
+      assert((int)perm.size() == (degree + 1) * (degree + 2) * (degree + 3) / 6);
+      break;
+    }
+    case cell::type::prism:
+    {
+      perm.push_back(0);
+      if (degree > 0)
+      {
+        int n = degree - 1;
+        int face0 = 6 + 9 * n;
+        int face1 = 6 + 9 * n + n * (n - 1) / 2;
+        int face2 = 6 + 9 * n + n * (n - 1) / 2 + n * n;
+        int face3 = 6 + 9 * n + n * (n - 1) / 2 + 2 * n * n;
+        int face4 = 6 + 9 * n + n * (n - 1) / 2 + 3 * n * n;
+        int interior = 6 + 9 * n + n * (n - 1) + 3 * n * n;
+        for (int i = 0; i < n; ++i)
+          perm.push_back(6 + i);
+        perm.push_back(1);
+        for (int i = 0; i < n; ++i)
+        {
+          perm.push_back(6 + n + i);
+          for (int j = 0; j < n - 1 - i; ++j)
+            perm.push_back(face0++);
+          perm.push_back(6 + 3 * n + i);
+        }
+        perm.push_back(2);
+        for (int i = 0; i < n; ++i)
+        {
+          perm.push_back(6 + 2 * n + i);
+          for (int j = 0; j < n; ++j)
+            perm.push_back(face1++);
+          perm.push_back(6 + 4 * n + i);
+          for (int j = 0; j < n; ++j)
+          {
+            perm.push_back(face2++);
+            for (int k = 0; k < n - 1 - j; ++k)
+              perm.push_back(interior++);
+            perm.push_back(face3++);
+          }
+          perm.push_back(6 + 5 * n + i);
+        }
+        perm.push_back(3);
+        for (int i = 0; i < n; ++i)
+          perm.push_back(6 + 6 * n + i);
+        perm.push_back(4);
+        for (int i = 0; i < n; ++i)
+        {
+          perm.push_back(6 + 7 * n + i);
+          for (int j = 0; j < n - 1 - i; ++j)
+            perm.push_back(face4++);
+          perm.push_back(6 + 8 * n + i);
+        }
+        perm.push_back(5);
+      }
+
+      assert((int)perm.size()
+             == (degree + 1) * (degree + 1) * (degree + 2) / 2);
+      break;
+    }
+
+    case cell::type::pyramid:
+    {
+      perm.push_back(0);
+      if (degree > 0)
+      {
+        int n = degree - 1;
+        int face0 = 5 + 8 * n;
+        int face1 = 5 + 8 * n + n * n;
+        int face2 = 5 + 8 * n + n * n + n * (n - 1) / 2;
+        int face3 = 5 + 8 * n + n * n + n * (n - 1);
+        int face4 = 5 + 8 * n + n * n + n * (n - 1) * 3 / 2;
+        int interior = 5 + 8 * n + n * n + n * (n - 1) * 2;
+        for (int i = 0; i < n; ++i)
+          perm.push_back(5 + i);
+        perm.push_back(1);
+        for (int i = 0; i < n; ++i)
+        {
+          perm.push_back(5 + n + i);
+          for (int j = 0; j < n; ++j)
+            perm.push_back(face0++);
+          perm.push_back(5 + 3 * n + i);
+        }
+        perm.push_back(2);
+        for (int i = 0; i < n; ++i)
+          perm.push_back(5 + 5 * n + i);
+        perm.push_back(3);
+        for (int i = 0; i < n; ++i)
+        {
+          perm.push_back(5 + 2 * n + i);
+          for (int j = 0; j < n - 1 - i; ++j)
+            perm.push_back(face1++);
+          perm.push_back(5 + 4 * n + i);
+          for (int j = 0; j < n - 1 - i; ++j)
+          {
+            perm.push_back(face2++);
+            for (int k = 0; k < n - 1 - i; ++k)
+              perm.push_back(interior++);
+            perm.push_back(face3++);
+          }
+          perm.push_back(5 + 6 * n + i);
+          for (int j = 0; j < n - 1 - i; ++j)
+            perm.push_back(face4++);
+          perm.push_back(5 + 7 * n + i);
+        }
+        perm.push_back(4);
+      }
+      assert((int)perm.size()
+             == (degree + 1) * (degree + 2) * (2 * degree + 3) / 6);
+      break;
+    }
+    default:
+    {
+    }
+    }
+    break;
+  }
+  default:
+  {
+  }
+  }
+
+  if (perm.size() == 0)
+  {
+    throw std::runtime_error(
+        "Lexicographic ordering not implemented for this element.");
+  }
+  dof_ordering.resize(perm.size());
+  for (std::size_t i = 0; i < perm.size(); ++i)
+    dof_ordering[perm[i]] = i;
+  return dof_ordering;
+}
+//-----------------------------------------------------------------------------
 template <std::floating_point T>
 std::tuple<std::array<std::vector<std::vector<T>>, 4>,
            std::array<std::vector<std::array<std::size_t, 2>>, 4>,
diff --git a/cpp/basix/finite-element.h b/cpp/basix/finite-element.h
index 58a255d9e..b56b46231 100644
--- a/cpp/basix/finite-element.h
+++ b/cpp/basix/finite-element.h
@@ -1629,6 +1629,21 @@ std::vector<int> tp_dof_ordering(element::family family, cell::type cell,
                                  element::dpc_variant dvariant,
                                  bool discontinuous);
 
+/// Get the lexicographic DOF ordering for an element
+/// @param[in] family The element family
+/// @param[in] cell The reference cell type that the element is defined on.
+/// @param[in] degree The degree of the element
+/// @param[in] lvariant The variant of Lagrange to use
+/// @param[in] dvariant The variant of DPC to use
+/// @param[in] discontinuous Indicates whether the element is discontinuous
+/// between cells points of the element. The discontinuous element will have the
+/// same DOFs, but they will all be associated with the interior of the cell.
+/// @return A vector containing the dof ordering
+std::vector<int> lex_dof_ordering(element::family family, cell::type cell,
+                                  int degree, element::lagrange_variant lvariant,
+                                  element::dpc_variant dvariant,
+                                  bool discontinuous);
+
 /// Get the tensor factors of an element
 /// @param[in] family The element family
 /// @param[in] cell The reference cell type that the element is defined on.
diff --git a/python/basix/_basixcpp.pyi b/python/basix/_basixcpp.pyi
index 3f6cbeb70..ce18bbf2a 100644
--- a/python/basix/_basixcpp.pyi
+++ b/python/basix/_basixcpp.pyi
@@ -520,4 +520,6 @@ def topology(arg: CellType, /) -> list[list[list[int]]]: ...
 
 def tp_dof_ordering(arg0: ElementFamily, arg1: CellType, arg2: int, arg3: LagrangeVariant, arg4: DPCVariant, arg5: bool, /) -> list[int]: ...
 
+def lex_dof_ordering(arg0: ElementFamily, arg1: CellType, arg2: int, arg3: LagrangeVariant, arg4: DPCVariant, arg5: bool, /) -> list[int]: ...
+
 def tp_factors(arg0: ElementFamily, arg1: CellType, arg2: int, arg3: LagrangeVariant, arg4: DPCVariant, arg5: bool, arg6: Sequence[int], arg7: str, /) -> list[list[FiniteElement_float32]] | list[list[FiniteElement_float64]]: ...
diff --git a/python/basix/finite_element.py b/python/basix/finite_element.py
index fbbb3fd22..093034808 100644
--- a/python/basix/finite_element.py
+++ b/python/basix/finite_element.py
@@ -23,6 +23,7 @@
 from basix._basixcpp import create_element as _create_element
 from basix._basixcpp import create_tp_element as _create_tp_element
 from basix._basixcpp import tp_dof_ordering as _tp_dof_ordering
+from basix._basixcpp import lex_dof_ordering as _lex_dof_ordering
 from basix._basixcpp import tp_factors as _tp_factors
 from basix.cell import CellType, geometry, topology, facet_outward_normals
 from basix import MapType
@@ -34,6 +35,7 @@
     "create_element",
     "create_custom_element",
     "create_tp_element",
+    "lex_dof_ordering",
     "string_to_family",
     "tp_factors",
     "tp_dof_ordering",
@@ -779,7 +781,7 @@ def tp_factors(
 ) -> list[list[FiniteElement]]:
     """Elements in the tensor product factorisation of an element.
 
-    If the element has no factorisation, an empty list is returned.
+    If the element has no factorisation, raises a RuntimeError.
 
     Args:
         family: Finite element family.
@@ -825,8 +827,7 @@ def tp_dof_ordering(
     This DOF ordering can be passed into create_element to create the
     element with DOFs ordered in a tensor product order.
 
-    If the element has no tensor product factorisation, an empty list is
-    returned.
+    If the element has no factorisation, raises a RuntimeError.
 
     Args:
         family: Finite element family.
@@ -852,6 +853,46 @@ def tp_dof_ordering(
     )
 
 
+def lex_dof_ordering(
+    family: ElementFamily,
+    celltype: CellType,
+    degree: int,
+    lagrange_variant: LagrangeVariant = LagrangeVariant.unset,
+    dpc_variant: DPCVariant = DPCVariant.unset,
+    discontinuous: bool = False,
+) -> list[int]:
+    """Lexicographic DOF ordering for an element.
+
+    This DOF ordering can be passed into create_element to create the
+    element with DOFs ordered in a lexicographic order.
+
+    If the element contains DOFs that are not point evaluations, raises a
+    RuntimeError.
+
+    Args:
+        family: Finite element family.
+        celltype: Reference cell type that the element is defined on
+        degree: Polynomial degree of the element.
+        lagrange_variant: Lagrange variant type.
+        dpc_variant: DPC variant type.
+        discontinuous: If `True` element is discontinuous. The
+            discontinuous element will have the same DOFs as a
+            continuous element, but the DOFs will all be associated with
+            the interior of the cell.
+
+    Returns:
+        The DOF ordering.
+    """
+    return _lex_dof_ordering(
+        family,
+        celltype,
+        degree,
+        lagrange_variant,
+        dpc_variant,
+        discontinuous,
+    )
+
+
 def string_to_family(family: str, cell: str) -> ElementFamily:
     """Basix ElementFamily enum representing the family type on the given cell.
 
diff --git a/python/wrapper.cpp b/python/wrapper.cpp
index dcd71b3ab..a42fcef6d 100644
--- a/python/wrapper.cpp
+++ b/python/wrapper.cpp
@@ -670,6 +670,17 @@ NB_MODULE(_basixcpp, m)
                                         discontinuous);
         });
 
+  m.def("lex_dof_ordering",
+        [](element::family family_name, cell::type cell, int degree,
+           element::lagrange_variant lagrange_variant,
+           element::dpc_variant dpc_variant,
+           bool discontinuous) -> std::vector<int>
+        {
+          return basix::lex_dof_ordering(family_name, cell, degree,
+                                         lagrange_variant, dpc_variant,
+                                         discontinuous);
+        });
+
   nb::enum_<polyset::type>(m, "PolysetType", nb::is_arithmetic(),
                            "Polyset type.")
       .value("standard", polyset::type::standard)
diff --git a/test/test_dof_ordering.py b/test/test_dof_ordering.py
index ed1024307..a84b7db34 100644
--- a/test/test_dof_ordering.py
+++ b/test/test_dof_ordering.py
@@ -56,4 +56,3 @@ def test_ordering(lagrange_variant):
     presult[:, :, order, :] = result2
 
     assert np.allclose(result1, presult)
-    print(result1)
diff --git a/test/test_lexicographic.py b/test/test_lexicographic.py
new file mode 100644
index 000000000..4ceecf23e
--- /dev/null
+++ b/test/test_lexicographic.py
@@ -0,0 +1,49 @@
+# Copyright (c) 2025 Matthew Scroggs
+# FEniCS Project
+# SPDX-License-Identifier: MIT
+
+from functools import cmp_to_key
+
+import numpy as np
+import pytest
+
+import basix
+
+
+@pytest.mark.parametrize(
+    "cell_type",
+    [
+        basix.CellType.interval,
+        basix.CellType.triangle,
+        basix.CellType.quadrilateral,
+        basix.CellType.tetrahedron,
+        basix.CellType.hexahedron,
+        basix.CellType.prism,
+        basix.CellType.pyramid,
+    ],
+)
+@pytest.mark.parametrize("degree", range(1, 10))
+def test_dof_ordering(cell_type, degree):
+    element = basix.finite_element.create_element(
+        basix.ElementFamily.P, cell_type, degree, basix.LagrangeVariant.equispaced
+    )
+
+    def cmp(x, y):
+        for i, j in zip(x[1][::-1], y[1][::-1]):
+            if not np.isclose(i, j):
+                if i > j:
+                    return 1
+                else:
+                    return -1
+        return 0
+
+    dof_points = list(enumerate(element.points))
+    dof_points.sort(key=cmp_to_key(cmp))
+
+    lex_order = [i[0] for i in dof_points]
+
+    ordering = basix.finite_element.lex_dof_ordering(
+        basix.ElementFamily.P, cell_type, degree, basix.LagrangeVariant.equispaced
+    )
+    for i, p in enumerate(lex_order):
+        assert ordering[p] == i