Add function that gives permutation to lexicographic ordering

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]]: ...

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>,

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.

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]]: ...