From 92d9cc7dfbd9f72dd1518ba279273f864c2ad22e Mon Sep 17 00:00:00 2001
From: Matthew Scroggs <matthew.w.scroggs@gmail.com>
Date: Fri, 19 Apr 2024 08:11:09 +0100
Subject: [PATCH 01/11] start on new function

---
 cpp/basix/finite-element.h | 20 ++++++++++++++++++++
 1 file changed, 20 insertions(+)

diff --git a/cpp/basix/finite-element.h b/cpp/basix/finite-element.h
index be269858f..d290e3c0e 100644
--- a/cpp/basix/finite-element.h
+++ b/cpp/basix/finite-element.h
@@ -805,6 +805,26 @@ class FiniteElement
     permute_data<std::int32_t, true>(dofs, 1, cell_info, _eperm_rev);
   }
 
+  /// Permute the dof numbering on the closure of a subentity of a cell
+  ///
+  /// @note This function is designed to be called at runtime, so its
+  /// performance is critical.
+  ///
+  /// @param[in,out] dofs The dof numbering for the closure of the subentity
+  /// @param cell_info The permutation info for the subentity
+  /// @param tdim The topological dimension of the subentity
+  void permute_closure_dofs(std::span<std::int32_t> dofs,
+                            std::uint32_t cell_info, int tdim) const
+  {
+    if (!_dof_transformations_are_permutations)
+    {
+      throw std::runtime_error(
+          "The DOF transformations for this element are not permutations");
+    }
+
+    if (_dof_transformations_are_identity)
+      return;
+  }
   /// Multiply data by DOF transformation matrix from the left
   ///
   /// @note This function is designed to be called at runtime, so its

From 7fdf05abd9f1396057f3237a507edca4503531bb Mon Sep 17 00:00:00 2001
From: Matthew Scroggs <matthew.w.scroggs@gmail.com>
Date: Fri, 19 Apr 2024 17:00:43 +0100
Subject: [PATCH 02/11] Working on subentity closure permutations

---
 cpp/basix/finite-element.h       | 117 ++++++++++++++++++++++++++++---
 test/test_dof_transformations.py |  76 ++++++++++++++++----
 2 files changed, 168 insertions(+), 25 deletions(-)

diff --git a/cpp/basix/finite-element.h b/cpp/basix/finite-element.h
index cc334c476..fb7fedf09 100644
--- a/cpp/basix/finite-element.h
+++ b/cpp/basix/finite-element.h
@@ -863,21 +863,118 @@ class FiniteElement
           "The DOF transformations for this element are not permutations");
     }
 
-    auto dofs = entity_closure_dofs()[entity_dim][entity_n];
-    std::cout << "dofs = { ";
-    for (auto i: dofs) std::cout << i << " ";
-    std::cout << "}\n";
-
+    // This is for a triangle
+
+    int dof_n = 0;
+    const auto ed = entity_dofs();
+    const auto conn
+        = cell::sub_entity_connectivity(_cell_type)[entity_dim][entity_n];
+    std::vector<std::vector<std::vector<int>>> dofs;
+    dofs.resize(entity_dim + 1);
     for (int dim = 0; dim <= entity_dim; ++dim)
     {
+      dofs[dim].resize(conn[dim].size());
+      for (std::size_t i = 0; i < conn[dim].size(); ++i)
+      {
+        const int e = conn[dim][i];
+        for (int j : ed[dim][e])
+        {
+          std::ignore = j;
+          dofs[dim][i].push_back(dof_n++);
+        }
+      }
+      std::cout << "<" << dim << " " << dof_n << ">\n";
+    }
+
+    std::vector<std::size_t> rot;
+    // Vertices
+    for (int i : dofs[0][1])
+      rot.push_back(i);
+    for (int i : dofs[0][2])
+      rot.push_back(i);
+    for (int i : dofs[0][0])
+      rot.push_back(i);
+    // Edges
+    for (int i : dofs[1][1])
+      rot.push_back(i);
+    for (int i : dofs[1][2])
+      rot.push_back(i);
+    for (int i : dofs[1][0])
+      rot.push_back(i);
+    // Face
+    for (int i : dofs[2][0])
+      rot.push_back(i);
+
+    std::vector<std::size_t> ref;
+    for (int i : dofs[0][0])
+      ref.push_back(i);
+    for (int i : dofs[0][2])
+      ref.push_back(i);
+    for (int i : dofs[0][1])
+      ref.push_back(i);
+    // Edges
+    for (int i : dofs[1][0])
+      ref.push_back(i);
+    for (int i : dofs[1][2])
+      ref.push_back(i);
+    for (int i : dofs[1][1])
+      ref.push_back(i);
+    // Face
+    for (int i : dofs[2][0])
+      ref.push_back(i);
+
+    if (!_dof_transformations_are_identity)
+    {
+      auto& trans1 = _eperm.at(cell::type::interval)[0];
+      auto& trans2 = _eperm.at(cell::type::triangle);
+
+      precompute::apply_permutation(trans1, std::span(rot), dofs[1][1][0]);
+      precompute::apply_permutation(trans1, std::span(rot), dofs[1][0][0]);
+      precompute::apply_permutation(trans2[0], std::span(rot), dofs[2][0][0]);
+
+      precompute::apply_permutation(trans1, std::span(ref), dofs[1][0][0]);
+      precompute::apply_permutation(trans2[1], std::span(ref), dofs[2][0][0]);
     }
 
-    d[0] = 13;
-    d[1] = 13;
-    std::cout << cell_info << entity_dim << entity_n << "\n";
+    precompute::prepare_permutation(rot);
+    precompute::prepare_permutation(ref);
+    // if (_dof_transformations_are_identity)
+    //   return;
 
-    if (_dof_transformations_are_identity)
-      return;
+    if (entity_dim == 1)
+    {
+      throw std::runtime_error("TODO");
+      // auto& trans = eperm.at(cell::type::interval)[0];
+      // for (std::size_t e = 0; e < _edofs[1].size(); ++e)
+      //{
+      //  Reverse an edge
+      // if (cell_info >> (face_start + e) & 1)
+      //{
+      //   precompute::apply_permutation_mapped(trans, data, _edofs[1][e],
+      //                                        block_size);
+      // }
+      //}
+    }
+    else if (entity_dim == 2)
+    {
+      // Rotate a face
+      for (std::uint32_t r = 0; r < (cell_info >> 1 & 3); ++r)
+      {
+        std::cout << r << "\n";
+        precompute::apply_permutation(rot, d);
+      }
+
+      // Reflect a face (post rotate)
+      if (cell_info & 1)
+      {
+        precompute::apply_permutation(ref, d);
+      }
+    }
+    else if (entity_dim != 0)
+    {
+      throw std::runtime_error(
+          "Invalid dimension for permute_subentity_closure");
+    }
   }
 
   /// @brief Transform basis functions from the reference element
diff --git a/test/test_dof_transformations.py b/test/test_dof_transformations.py
index 1cb4dda07..799d2f9cb 100644
--- a/test/test_dof_transformations.py
+++ b/test/test_dof_transformations.py
@@ -495,23 +495,69 @@ def test_transformation_of_tabulated_data_pyramid(element_type, degree, element_
                     (bt[9].dot(i_slice))[start : start + ndofs], j_slice[start : start + ndofs]
                 )
 
+
 @pytest.mark.parametrize(
-    "family, cell_Type, degree, args, cell_info, entity_dim, entity_n, result", [
-        (basix.ElementFamily.P, basix.CellType.tetrahedron, 1, (), 0, 2, 0, [0, 1, 2]),
-        (basix.ElementFamily.P, basix.CellType.tetrahedron, 1, (), 1, 2, 0, [1, 2, 0]),
-        (basix.ElementFamily.P, basix.CellType.tetrahedron, 1, (), 2, 2, 0, [2, 0, 1]),
-        (basix.ElementFamily.P, basix.CellType.tetrahedron, 1, (), 4, 2, 0, [0, 2, 1]),
-        (basix.ElementFamily.P, basix.CellType.tetrahedron, 1, (), 5, 2, 0, [1, 0, 2]),
-        (basix.ElementFamily.P, basix.CellType.tetrahedron, 1, (), 6, 2, 0, [2, 1, 0]),
-    ]
+    "family, cell_Type, degree, args, entity_dim, entity_n, results",
+    [
+        (
+            basix.ElementFamily.P,
+            basix.CellType.tetrahedron,
+            1,
+            (),
+            2,
+            0,
+            {
+                0: [0, 1, 2],
+                2: [1, 2, 0],
+                4: [2, 0, 1],
+                1: [0, 2, 1],
+                3: [1, 0, 2],
+                5: [2, 1, 0],
+            },
+        ),
+        (
+            basix.ElementFamily.P,
+            basix.CellType.tetrahedron,
+            2,
+            (),
+            2,
+            0,
+            {
+                0: [0, 1, 2, 3, 4, 5],
+                2: [1, 2, 0, 4, 5, 3],
+                4: [2, 0, 1, 5, 3, 4],
+                1: [0, 2, 1, 3, 5, 4],
+                3: [1, 0, 2, 4, 3, 5],
+                5: [2, 1, 0, 5, 4, 3],
+            },
+        ),
+        (
+            basix.ElementFamily.P,
+            basix.CellType.tetrahedron,
+            3,
+            (basix.LagrangeVariant.equispaced,),
+            2,
+            0,
+            {
+                0: [0, 1, 2, 3, 4, 5, 6, 7, 8],
+                2: [1, 2, 0, 6, 5, 8, 7, 3, 4],
+                4: [2, 0, 1, 7, 8, 4, 3, 6, 5],
+                1: [0, 2, 1, 4, 3, 7, 8, 5, 6],
+                3: [1, 0, 2, 5, 6, 3, 4, 8, 7],
+                5: [2, 1, 0, 8, 7, 6, 5, 4, 3],
+            },
+        ),
+    ],
 )
-def test_permute_subentity_closure(
-    family, cell_Type, degree, args, cell_info, entity_dim, entity_n, result
-):
+def test_permute_subentity_closure(family, cell_Type, degree, args, entity_dim, entity_n, results):
     e = basix.create_element(family, cell_Type, degree, *args)
 
-    data = np.array(list(range(len(result))))
-    data = e._e.permute_subentity_closure(data, cell_info, entity_dim, entity_n)
+    for cell_info, result in results.items():
+        data = np.array(list(range(len(result))))
+        data = e._e.permute_subentity_closure(data, cell_info, entity_dim, entity_n)
 
-    for i, j in zip(data, result):
-        assert i == j
+        print(cell_info)
+        print(data)
+        print(result)
+        for i, j in zip(data, result):
+            assert i == j

From 083b31d38e191fc1181665d996b395b94c4946cf Mon Sep 17 00:00:00 2001
From: Matthew Scroggs <matthew.w.scroggs@gmail.com>
Date: Fri, 19 Apr 2024 17:50:32 +0100
Subject: [PATCH 03/11] subentity_closure_perm for all elements

---
 cpp/basix/finite-element.cpp     | 220 +++++++++++++++++++++++++++++++
 cpp/basix/finite-element.h       | 115 +++-------------
 python/wrapper.cpp               |   6 +-
 test/test_dof_transformations.py |  56 +++++---
 4 files changed, 282 insertions(+), 115 deletions(-)

diff --git a/cpp/basix/finite-element.cpp b/cpp/basix/finite-element.cpp
index 8df56b43a..6285438bb 100644
--- a/cpp/basix/finite-element.cpp
+++ b/cpp/basix/finite-element.cpp
@@ -1133,6 +1133,226 @@ FiniteElement<F>::FiniteElement(
     }
   }
 
+  // If DOF transformations are permutations, compute the subentity closure
+  // permutations
+  if (_dof_transformations_are_permutations)
+  {
+    if (_cell_tdim > 1)
+    {
+      // interval
+      int dof_n = 0;
+      const auto conn = cell::sub_entity_connectivity(_cell_type)[1][0];
+      std::vector<std::vector<std::vector<int>>> dofs;
+      dofs.resize(2);
+      for (int dim = 0; dim <= 1; ++dim)
+      {
+        dofs[dim].resize(conn[dim].size());
+        for (std::size_t i = 0; i < conn[dim].size(); ++i)
+        {
+          const int e = conn[dim][i];
+          for (int j : _edofs[dim][e])
+          {
+            std::ignore = j;
+            dofs[dim][i].push_back(dof_n++);
+          }
+        }
+      }
+
+      std::vector<std::size_t> ref;
+      for (int i : dofs[0][1])
+        ref.push_back(i);
+      for (int i : dofs[0][0])
+        ref.push_back(i);
+      // Edges
+      for (int i : dofs[1][0])
+        ref.push_back(i);
+
+      if (!_dof_transformations_are_identity)
+      {
+        auto& trans1 = _eperm.at(cell::type::interval)[0];
+        precompute::apply_permutation(trans1, std::span(ref), dofs[1][0][0]);
+      }
+
+      precompute::prepare_permutation(ref);
+
+      auto& secp = _subentity_closure_perm.try_emplace(cell::type::interval)
+                       .first->second;
+      secp.push_back(ref);
+    }
+    if (_cell_type == cell::type::tetrahedron || cell_type == cell::type::prism
+        || cell_type == cell::type::pyramid)
+    {
+      // triangle
+      const int face_n = cell_type == cell::type::pyramid ? 1 : 0;
+
+      int dof_n = 0;
+      const auto conn = cell::sub_entity_connectivity(_cell_type)[2][face_n];
+      std::vector<std::vector<std::vector<int>>> dofs;
+      dofs.resize(3);
+      for (int dim = 0; dim <= 2; ++dim)
+      {
+        dofs[dim].resize(conn[dim].size());
+        for (std::size_t i = 0; i < conn[dim].size(); ++i)
+        {
+          const int e = conn[dim][i];
+          for (int j : _edofs[dim][e])
+          {
+            std::ignore = j;
+            dofs[dim][i].push_back(dof_n++);
+          }
+        }
+      }
+
+      std::vector<std::size_t> rot;
+      // Vertices
+      for (int i : dofs[0][1])
+        rot.push_back(i);
+      for (int i : dofs[0][2])
+        rot.push_back(i);
+      for (int i : dofs[0][0])
+        rot.push_back(i);
+      // Edges
+      for (int i : dofs[1][1])
+        rot.push_back(i);
+      for (int i : dofs[1][2])
+        rot.push_back(i);
+      for (int i : dofs[1][0])
+        rot.push_back(i);
+      // Face
+      for (int i : dofs[2][0])
+        rot.push_back(i);
+
+      std::vector<std::size_t> ref;
+      for (int i : dofs[0][0])
+        ref.push_back(i);
+      for (int i : dofs[0][2])
+        ref.push_back(i);
+      for (int i : dofs[0][1])
+        ref.push_back(i);
+      // Edges
+      for (int i : dofs[1][0])
+        ref.push_back(i);
+      for (int i : dofs[1][2])
+        ref.push_back(i);
+      for (int i : dofs[1][1])
+        ref.push_back(i);
+      // Face
+      for (int i : dofs[2][0])
+        ref.push_back(i);
+
+      if (!_dof_transformations_are_identity)
+      {
+        auto& trans1 = _eperm.at(cell::type::interval)[0];
+        auto& trans2 = _eperm.at(cell::type::triangle);
+
+        precompute::apply_permutation(trans1, std::span(rot), dofs[1][1][0]);
+        precompute::apply_permutation(trans1, std::span(rot), dofs[1][0][0]);
+        precompute::apply_permutation(trans2[0], std::span(rot), dofs[2][0][0]);
+
+        precompute::apply_permutation(trans1, std::span(ref), dofs[1][0][0]);
+        precompute::apply_permutation(trans2[1], std::span(ref), dofs[2][0][0]);
+      }
+
+      precompute::prepare_permutation(rot);
+      precompute::prepare_permutation(ref);
+
+      auto& secp = _subentity_closure_perm.try_emplace(cell::type::triangle)
+                       .first->second;
+      secp.push_back(rot);
+      secp.push_back(ref);
+    }
+    if (_cell_type == cell::type::hexahedron || cell_type == cell::type::prism
+        || cell_type == cell::type::pyramid)
+    {
+      // quadrilateral
+      const int face_n = cell_type == cell::type::prism ? 1 : 0;
+
+      int dof_n = 0;
+      const auto conn = cell::sub_entity_connectivity(_cell_type)[2][face_n];
+      std::vector<std::vector<std::vector<int>>> dofs;
+      dofs.resize(3);
+      for (int dim = 0; dim <= 2; ++dim)
+      {
+        dofs[dim].resize(conn[dim].size());
+        for (std::size_t i = 0; i < conn[dim].size(); ++i)
+        {
+          const int e = conn[dim][i];
+          for (int j : _edofs[dim][e])
+          {
+            std::ignore = j;
+            dofs[dim][i].push_back(dof_n++);
+          }
+        }
+      }
+
+      std::vector<std::size_t> rot;
+      // Vertices
+      for (int i : dofs[0][1])
+        rot.push_back(i);
+      for (int i : dofs[0][3])
+        rot.push_back(i);
+      for (int i : dofs[0][0])
+        rot.push_back(i);
+      for (int i : dofs[0][2])
+        rot.push_back(i);
+      // Edges
+      for (int i : dofs[1][2])
+        rot.push_back(i);
+      for (int i : dofs[1][0])
+        rot.push_back(i);
+      for (int i : dofs[1][3])
+        rot.push_back(i);
+      for (int i : dofs[1][1])
+        rot.push_back(i);
+      // Face
+      for (int i : dofs[2][0])
+        rot.push_back(i);
+
+      std::vector<std::size_t> ref;
+      for (int i : dofs[0][0])
+        ref.push_back(i);
+      for (int i : dofs[0][2])
+        ref.push_back(i);
+      for (int i : dofs[0][1])
+        ref.push_back(i);
+      for (int i : dofs[0][3])
+        ref.push_back(i);
+      // Edges
+      for (int i : dofs[1][1])
+        ref.push_back(i);
+      for (int i : dofs[1][0])
+        ref.push_back(i);
+      for (int i : dofs[1][3])
+        ref.push_back(i);
+      for (int i : dofs[1][2])
+        ref.push_back(i);
+      // Face
+      for (int i : dofs[2][0])
+        ref.push_back(i);
+
+      if (!_dof_transformations_are_identity)
+      {
+        auto& trans1 = _eperm.at(cell::type::interval)[0];
+        auto& trans2 = _eperm.at(cell::type::quadrilateral);
+
+        precompute::apply_permutation(trans1, std::span(rot), dofs[1][1][0]);
+        precompute::apply_permutation(trans1, std::span(rot), dofs[1][2][0]);
+        precompute::apply_permutation(trans2[0], std::span(rot), dofs[2][0][0]);
+
+        precompute::apply_permutation(trans2[1], std::span(ref), dofs[2][0][0]);
+      }
+
+      precompute::prepare_permutation(rot);
+      precompute::prepare_permutation(ref);
+
+      auto& secp
+          = _subentity_closure_perm.try_emplace(cell::type::quadrilateral)
+                .first->second;
+      secp.push_back(rot);
+      secp.push_back(ref);
+    }
+  }
+
   // Check if interpolation matrix is the identity
   mdspan_t<const F, 2> matM(_matM.first.data(), _matM.second);
   _interpolation_is_identity = matM.extent(0) == matM.extent(1);
diff --git a/cpp/basix/finite-element.h b/cpp/basix/finite-element.h
index fb7fedf09..15f8adf75 100644
--- a/cpp/basix/finite-element.h
+++ b/cpp/basix/finite-element.h
@@ -23,8 +23,6 @@
 #include <utility>
 #include <vector>
 
-#include <iostream>
-
 /// Basix: FEniCS runtime basis evaluation library
 namespace basix
 {
@@ -852,10 +850,10 @@ class FiniteElement
   /// degree-of-freedom (in). Indices associated with each physical
   /// element degree-of-freedom (out).
   /// @param cell_info Permutation info for the cell
-  /// @param entity_dim The dimension of the sub-entity
-  /// @param entity_n The number of the sub-entity
+  /// @param entity_type The cell type of the sub-entity
   void permute_subentity_closure(std::span<std::int32_t> d,
-                                 std::uint32_t cell_info, int entity_dim, int entity_n) const
+                                 std::uint32_t cell_info,
+                                 cell::type entity_type) const
   {
     if (!_dof_transformations_are_permutations)
     {
@@ -863,114 +861,34 @@ class FiniteElement
           "The DOF transformations for this element are not permutations");
     }
 
-    // This is for a triangle
-
-    int dof_n = 0;
-    const auto ed = entity_dofs();
-    const auto conn
-        = cell::sub_entity_connectivity(_cell_type)[entity_dim][entity_n];
-    std::vector<std::vector<std::vector<int>>> dofs;
-    dofs.resize(entity_dim + 1);
-    for (int dim = 0; dim <= entity_dim; ++dim)
-    {
-      dofs[dim].resize(conn[dim].size());
-      for (std::size_t i = 0; i < conn[dim].size(); ++i)
-      {
-        const int e = conn[dim][i];
-        for (int j : ed[dim][e])
-        {
-          std::ignore = j;
-          dofs[dim][i].push_back(dof_n++);
-        }
-      }
-      std::cout << "<" << dim << " " << dof_n << ">\n";
-    }
-
-    std::vector<std::size_t> rot;
-    // Vertices
-    for (int i : dofs[0][1])
-      rot.push_back(i);
-    for (int i : dofs[0][2])
-      rot.push_back(i);
-    for (int i : dofs[0][0])
-      rot.push_back(i);
-    // Edges
-    for (int i : dofs[1][1])
-      rot.push_back(i);
-    for (int i : dofs[1][2])
-      rot.push_back(i);
-    for (int i : dofs[1][0])
-      rot.push_back(i);
-    // Face
-    for (int i : dofs[2][0])
-      rot.push_back(i);
-
-    std::vector<std::size_t> ref;
-    for (int i : dofs[0][0])
-      ref.push_back(i);
-    for (int i : dofs[0][2])
-      ref.push_back(i);
-    for (int i : dofs[0][1])
-      ref.push_back(i);
-    // Edges
-    for (int i : dofs[1][0])
-      ref.push_back(i);
-    for (int i : dofs[1][2])
-      ref.push_back(i);
-    for (int i : dofs[1][1])
-      ref.push_back(i);
-    // Face
-    for (int i : dofs[2][0])
-      ref.push_back(i);
-
-    if (!_dof_transformations_are_identity)
-    {
-      auto& trans1 = _eperm.at(cell::type::interval)[0];
-      auto& trans2 = _eperm.at(cell::type::triangle);
+    const int entity_dim = cell::topological_dimension(entity_type);
 
-      precompute::apply_permutation(trans1, std::span(rot), dofs[1][1][0]);
-      precompute::apply_permutation(trans1, std::span(rot), dofs[1][0][0]);
-      precompute::apply_permutation(trans2[0], std::span(rot), dofs[2][0][0]);
-
-      precompute::apply_permutation(trans1, std::span(ref), dofs[1][0][0]);
-      precompute::apply_permutation(trans2[1], std::span(ref), dofs[2][0][0]);
-    }
-
-    precompute::prepare_permutation(rot);
-    precompute::prepare_permutation(ref);
-    // if (_dof_transformations_are_identity)
-    //   return;
+    if (entity_dim == 0)
+      return;
 
+    auto& perm = _subentity_closure_perm.at(entity_type);
     if (entity_dim == 1)
     {
-      throw std::runtime_error("TODO");
-      // auto& trans = eperm.at(cell::type::interval)[0];
-      // for (std::size_t e = 0; e < _edofs[1].size(); ++e)
-      //{
-      //  Reverse an edge
-      // if (cell_info >> (face_start + e) & 1)
-      //{
-      //   precompute::apply_permutation_mapped(trans, data, _edofs[1][e],
-      //                                        block_size);
-      // }
-      //}
+      if (cell_info & 1)
+      {
+        precompute::apply_permutation(perm[0], d);
+      }
     }
     else if (entity_dim == 2)
     {
       // Rotate a face
       for (std::uint32_t r = 0; r < (cell_info >> 1 & 3); ++r)
       {
-        std::cout << r << "\n";
-        precompute::apply_permutation(rot, d);
+        precompute::apply_permutation(perm[0], d);
       }
 
       // Reflect a face (post rotate)
       if (cell_info & 1)
       {
-        precompute::apply_permutation(ref, d);
+        precompute::apply_permutation(perm[1], d);
       }
     }
-    else if (entity_dim != 0)
+    else
     {
       throw std::runtime_error(
           "Invalid dimension for permute_subentity_closure");
@@ -1493,6 +1411,11 @@ class FiniteElement
   // The inverse transpose entity transformations in precomputed form
   std::map<cell::type, trans_data_t> _etrans_invT;
 
+  // The subentity closure permutations (factorised). This will only be set if
+  // _dof_transformations_are_permutations is True
+  std::map<cell::type, std::vector<std::vector<std::size_t>>>
+      _subentity_closure_perm;
+
   // Indicates whether or not this is the discontinuous version of the
   // element
   bool _discontinuous;
diff --git a/python/wrapper.cpp b/python/wrapper.cpp
index 5e54c9a84..45a93a2f7 100644
--- a/python/wrapper.cpp
+++ b/python/wrapper.cpp
@@ -101,12 +101,10 @@ void declare_float(nb::module_& m, std::string type)
       .def("permute_subentity_closure",
            [](const FiniteElement<T>& self,
               nb::ndarray<std::int32_t, nb::ndim<1>, nb::c_contig> d,
-              std::uint32_t cell_info,
-              int entity_dim,
-              int entity_n)
+              std::uint32_t cell_info, cell::type entity_type)
            {
              std::span<std::int32_t> _d(d.data(), d.shape(0));
-             self.permute_subentity_closure(_d, cell_info, entity_dim, entity_n);
+             self.permute_subentity_closure(_d, cell_info, entity_type);
              return as_nbarray(std::move(_d));
            })
       .def("push_forward",
diff --git a/test/test_dof_transformations.py b/test/test_dof_transformations.py
index 799d2f9cb..008bce0f1 100644
--- a/test/test_dof_transformations.py
+++ b/test/test_dof_transformations.py
@@ -497,15 +497,14 @@ def test_transformation_of_tabulated_data_pyramid(element_type, degree, element_
 
 
 @pytest.mark.parametrize(
-    "family, cell_Type, degree, args, entity_dim, entity_n, results",
+    "family, cell_type, degree, args, subentity, results",
     [
         (
             basix.ElementFamily.P,
             basix.CellType.tetrahedron,
             1,
             (),
-            2,
-            0,
+            basix.CellType.triangle,
             {
                 0: [0, 1, 2],
                 2: [1, 2, 0],
@@ -520,8 +519,7 @@ def test_transformation_of_tabulated_data_pyramid(element_type, degree, element_
             basix.CellType.tetrahedron,
             2,
             (),
-            2,
-            0,
+            basix.CellType.triangle,
             {
                 0: [0, 1, 2, 3, 4, 5],
                 2: [1, 2, 0, 4, 5, 3],
@@ -536,8 +534,7 @@ def test_transformation_of_tabulated_data_pyramid(element_type, degree, element_
             basix.CellType.tetrahedron,
             3,
             (basix.LagrangeVariant.equispaced,),
-            2,
-            0,
+            basix.CellType.triangle,
             {
                 0: [0, 1, 2, 3, 4, 5, 6, 7, 8],
                 2: [1, 2, 0, 6, 5, 8, 7, 3, 4],
@@ -547,17 +544,46 @@ def test_transformation_of_tabulated_data_pyramid(element_type, degree, element_
                 5: [2, 1, 0, 8, 7, 6, 5, 4, 3],
             },
         ),
+        (
+            basix.ElementFamily.P,
+            basix.CellType.tetrahedron,
+            1,
+            (basix.LagrangeVariant.equispaced,),
+            basix.CellType.interval,
+            {
+                0: [0, 1],
+                1: [1, 0],
+            },
+        ),
+        (
+            basix.ElementFamily.P,
+            basix.CellType.tetrahedron,
+            2,
+            (basix.LagrangeVariant.equispaced,),
+            basix.CellType.interval,
+            {
+                0: [0, 1, 2],
+                1: [1, 0, 2],
+            },
+        ),
+        (
+            basix.ElementFamily.P,
+            basix.CellType.tetrahedron,
+            3,
+            (basix.LagrangeVariant.equispaced,),
+            basix.CellType.interval,
+            {
+                0: [0, 1, 2, 3],
+                1: [1, 0, 3, 2],
+            },
+        ),
     ],
 )
-def test_permute_subentity_closure(family, cell_Type, degree, args, entity_dim, entity_n, results):
-    e = basix.create_element(family, cell_Type, degree, *args)
+def test_permute_subentity_closure(family, cell_type, degree, args, subentity, results):
+    e = basix.create_element(family, cell_type, degree, *args)
 
     for cell_info, result in results.items():
         data = np.array(list(range(len(result))))
-        data = e._e.permute_subentity_closure(data, cell_info, entity_dim, entity_n)
+        data = list(e._e.permute_subentity_closure(data, cell_info, subentity.value))
 
-        print(cell_info)
-        print(data)
-        print(result)
-        for i, j in zip(data, result):
-            assert i == j
+        assert result == data

From fe4e0ec4a46571733b6d86b5e6a268091486c606 Mon Sep 17 00:00:00 2001
From: Matthew Scroggs <matthew.w.scroggs@gmail.com>
Date: Mon, 22 Apr 2024 10:18:40 +0100
Subject: [PATCH 04/11] don't try to permute when no dofs are on the entity

---
 cpp/basix/finite-element.cpp | 38 ++++++++++++++++++++++++------------
 1 file changed, 26 insertions(+), 12 deletions(-)

diff --git a/cpp/basix/finite-element.cpp b/cpp/basix/finite-element.cpp
index 6285438bb..fbca4e22a 100644
--- a/cpp/basix/finite-element.cpp
+++ b/cpp/basix/finite-element.cpp
@@ -1170,7 +1170,8 @@ FiniteElement<F>::FiniteElement(
       if (!_dof_transformations_are_identity)
       {
         auto& trans1 = _eperm.at(cell::type::interval)[0];
-        precompute::apply_permutation(trans1, std::span(ref), dofs[1][0][0]);
+        if (dofs[1][0].size() > 0)
+          precompute::apply_permutation(trans1, std::span(ref), dofs[1][0][0]);
       }
 
       precompute::prepare_permutation(ref);
@@ -1245,12 +1246,19 @@ FiniteElement<F>::FiniteElement(
         auto& trans1 = _eperm.at(cell::type::interval)[0];
         auto& trans2 = _eperm.at(cell::type::triangle);
 
-        precompute::apply_permutation(trans1, std::span(rot), dofs[1][1][0]);
-        precompute::apply_permutation(trans1, std::span(rot), dofs[1][0][0]);
-        precompute::apply_permutation(trans2[0], std::span(rot), dofs[2][0][0]);
-
-        precompute::apply_permutation(trans1, std::span(ref), dofs[1][0][0]);
-        precompute::apply_permutation(trans2[1], std::span(ref), dofs[2][0][0]);
+        if (dofs[1][0].size() > 0)
+          precompute::apply_permutation(trans1, std::span(rot), dofs[1][0][0]);
+        if (dofs[1][1].size() > 0)
+          precompute::apply_permutation(trans1, std::span(rot), dofs[1][1][0]);
+        if (dofs[2][0].size() > 0)
+          precompute::apply_permutation(trans2[0], std::span(rot),
+                                        dofs[2][0][0]);
+
+        if (dofs[1][0].size() > 0)
+          precompute::apply_permutation(trans1, std::span(ref), dofs[1][0][0]);
+        if (dofs[2][0].size() > 0)
+          precompute::apply_permutation(trans2[1], std::span(ref),
+                                        dofs[2][0][0]);
       }
 
       precompute::prepare_permutation(rot);
@@ -1335,11 +1343,17 @@ FiniteElement<F>::FiniteElement(
         auto& trans1 = _eperm.at(cell::type::interval)[0];
         auto& trans2 = _eperm.at(cell::type::quadrilateral);
 
-        precompute::apply_permutation(trans1, std::span(rot), dofs[1][1][0]);
-        precompute::apply_permutation(trans1, std::span(rot), dofs[1][2][0]);
-        precompute::apply_permutation(trans2[0], std::span(rot), dofs[2][0][0]);
-
-        precompute::apply_permutation(trans2[1], std::span(ref), dofs[2][0][0]);
+        if (dofs[1][1].size() > 0)
+          precompute::apply_permutation(trans1, std::span(rot), dofs[1][1][0]);
+        if (dofs[1][2].size() > 0)
+          precompute::apply_permutation(trans1, std::span(rot), dofs[1][2][0]);
+        if (dofs[2][0].size() > 0)
+          precompute::apply_permutation(trans2[0], std::span(rot),
+                                        dofs[2][0][0]);
+
+        if (dofs[2][0].size() > 0)
+          precompute::apply_permutation(trans2[1], std::span(ref),
+                                        dofs[2][0][0]);
       }
 
       precompute::prepare_permutation(rot);

From fbd189066faa664126366a141596f7f8ef79090b Mon Sep 17 00:00:00 2001
From: Matthew Scroggs <matthew.w.scroggs@gmail.com>
Date: Mon, 22 Apr 2024 10:19:33 +0100
Subject: [PATCH 05/11] {}

---
 cpp/basix/finite-element.cpp | 21 ++++++++++++++++++++-
 1 file changed, 20 insertions(+), 1 deletion(-)

diff --git a/cpp/basix/finite-element.cpp b/cpp/basix/finite-element.cpp
index fbca4e22a..a3e93f238 100644
--- a/cpp/basix/finite-element.cpp
+++ b/cpp/basix/finite-element.cpp
@@ -1171,7 +1171,9 @@ FiniteElement<F>::FiniteElement(
       {
         auto& trans1 = _eperm.at(cell::type::interval)[0];
         if (dofs[1][0].size() > 0)
+        {
           precompute::apply_permutation(trans1, std::span(ref), dofs[1][0][0]);
+        }
       }
 
       precompute::prepare_permutation(ref);
@@ -1247,18 +1249,27 @@ FiniteElement<F>::FiniteElement(
         auto& trans2 = _eperm.at(cell::type::triangle);
 
         if (dofs[1][0].size() > 0)
+        {
           precompute::apply_permutation(trans1, std::span(rot), dofs[1][0][0]);
+        }
         if (dofs[1][1].size() > 0)
+        {
           precompute::apply_permutation(trans1, std::span(rot), dofs[1][1][0]);
+        }
         if (dofs[2][0].size() > 0)
+        {
           precompute::apply_permutation(trans2[0], std::span(rot),
                                         dofs[2][0][0]);
-
+        }
         if (dofs[1][0].size() > 0)
+        {
           precompute::apply_permutation(trans1, std::span(ref), dofs[1][0][0]);
+        }
         if (dofs[2][0].size() > 0)
+        {
           precompute::apply_permutation(trans2[1], std::span(ref),
                                         dofs[2][0][0]);
+        }
       }
 
       precompute::prepare_permutation(rot);
@@ -1344,16 +1355,24 @@ FiniteElement<F>::FiniteElement(
         auto& trans2 = _eperm.at(cell::type::quadrilateral);
 
         if (dofs[1][1].size() > 0)
+        {
           precompute::apply_permutation(trans1, std::span(rot), dofs[1][1][0]);
+        }
         if (dofs[1][2].size() > 0)
+        {
           precompute::apply_permutation(trans1, std::span(rot), dofs[1][2][0]);
+        }
         if (dofs[2][0].size() > 0)
+        {
           precompute::apply_permutation(trans2[0], std::span(rot),
                                         dofs[2][0][0]);
+        }
 
         if (dofs[2][0].size() > 0)
+        {
           precompute::apply_permutation(trans2[1], std::span(ref),
                                         dofs[2][0][0]);
+        }
       }
 
       precompute::prepare_permutation(rot);

From af54479d035e67c658786bc7cc14dfcd7904d65a Mon Sep 17 00:00:00 2001
From: Matthew Scroggs <matthew.w.scroggs@gmail.com>
Date: Thu, 25 Apr 2024 11:02:38 +0100
Subject: [PATCH 06/11] Update cpp/basix/finite-element.cpp

Co-authored-by: Garth N. Wells <gnw20@cam.ac.uk>
---
 cpp/basix/finite-element.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/cpp/basix/finite-element.cpp b/cpp/basix/finite-element.cpp
index a3e93f238..93f997f92 100644
--- a/cpp/basix/finite-element.cpp
+++ b/cpp/basix/finite-element.cpp
@@ -1170,7 +1170,7 @@ FiniteElement<F>::FiniteElement(
       if (!_dof_transformations_are_identity)
       {
         auto& trans1 = _eperm.at(cell::type::interval)[0];
-        if (dofs[1][0].size() > 0)
+        if (!dofs[1][0].empty())
         {
           precompute::apply_permutation(trans1, std::span(ref), dofs[1][0][0]);
         }

From 3af8ea3955cbe4081c3b95246d5ecebe8ae2f4f5 Mon Sep 17 00:00:00 2001
From: Matthew Scroggs <matthew.w.scroggs@gmail.com>
Date: Thu, 25 Apr 2024 11:04:06 +0100
Subject: [PATCH 07/11] more !*.empty()

---
 cpp/basix/finite-element.cpp | 22 +++++++++++-----------
 1 file changed, 11 insertions(+), 11 deletions(-)

diff --git a/cpp/basix/finite-element.cpp b/cpp/basix/finite-element.cpp
index 93f997f92..3e22e0730 100644
--- a/cpp/basix/finite-element.cpp
+++ b/cpp/basix/finite-element.cpp
@@ -230,7 +230,7 @@ basix::create_element(element::family family, cell::type cell, int degree,
       throw std::runtime_error("Cannot pass a DPC variant to this element.");
   }
 
-  if (dof_ordering.size() > 0 and family != element::family::P)
+  if (!dof_ordering.empty() and family != element::family::P)
   {
     throw std::runtime_error("DOF ordering only supported for Lagrange");
   }
@@ -353,7 +353,7 @@ basix::tp_factors(element::family family, cell::type cell, int degree,
 {
   std::vector<int> tp_dofs = tp_dof_ordering(family, cell, degree, lvariant,
                                              dvariant, discontinuous);
-  if (tp_dofs.size() > 0 && tp_dofs == dof_ordering)
+  if (!tp_dofs.empty() && tp_dofs == dof_ordering)
   {
     switch (family)
     {
@@ -1248,24 +1248,24 @@ FiniteElement<F>::FiniteElement(
         auto& trans1 = _eperm.at(cell::type::interval)[0];
         auto& trans2 = _eperm.at(cell::type::triangle);
 
-        if (dofs[1][0].size() > 0)
+        if (!dofs[1][0].empty())
         {
           precompute::apply_permutation(trans1, std::span(rot), dofs[1][0][0]);
         }
-        if (dofs[1][1].size() > 0)
+        if (!dofs[1][1].empty())
         {
           precompute::apply_permutation(trans1, std::span(rot), dofs[1][1][0]);
         }
-        if (dofs[2][0].size() > 0)
+        if (!dofs[2][0].empty())
         {
           precompute::apply_permutation(trans2[0], std::span(rot),
                                         dofs[2][0][0]);
         }
-        if (dofs[1][0].size() > 0)
+        if (!dofs[1][0].empty())
         {
           precompute::apply_permutation(trans1, std::span(ref), dofs[1][0][0]);
         }
-        if (dofs[2][0].size() > 0)
+        if (!dofs[2][0].empty())
         {
           precompute::apply_permutation(trans2[1], std::span(ref),
                                         dofs[2][0][0]);
@@ -1354,21 +1354,21 @@ FiniteElement<F>::FiniteElement(
         auto& trans1 = _eperm.at(cell::type::interval)[0];
         auto& trans2 = _eperm.at(cell::type::quadrilateral);
 
-        if (dofs[1][1].size() > 0)
+        if (!dofs[1][1].empty())
         {
           precompute::apply_permutation(trans1, std::span(rot), dofs[1][1][0]);
         }
-        if (dofs[1][2].size() > 0)
+        if (!dofs[1][2].empty())
         {
           precompute::apply_permutation(trans1, std::span(rot), dofs[1][2][0]);
         }
-        if (dofs[2][0].size() > 0)
+        if (!dofs[2][0].empty())
         {
           precompute::apply_permutation(trans2[0], std::span(rot),
                                         dofs[2][0][0]);
         }
 
-        if (dofs[2][0].size() > 0)
+        if (!dofs[2][0].empty())
         {
           precompute::apply_permutation(trans2[1], std::span(ref),
                                         dofs[2][0][0]);

From 28864faec7f8df99ca919b15577ba571e028858f Mon Sep 17 00:00:00 2001
From: Matthew Scroggs <matthew.w.scroggs@gmail.com>
Date: Wed, 1 May 2024 15:15:23 +0100
Subject: [PATCH 08/11] Code simplifications

---
 cpp/basix/finite-element.cpp | 109 ++++++++++++-----------------------
 1 file changed, 37 insertions(+), 72 deletions(-)

diff --git a/cpp/basix/finite-element.cpp b/cpp/basix/finite-element.cpp
index 3e22e0730..05b39eec3 100644
--- a/cpp/basix/finite-element.cpp
+++ b/cpp/basix/finite-element.cpp
@@ -1150,22 +1150,18 @@ FiniteElement<F>::FiniteElement(
         for (std::size_t i = 0; i < conn[dim].size(); ++i)
         {
           const int e = conn[dim][i];
-          for (int j : _edofs[dim][e])
+          for (std::size_t j = 0; j < _edofs[dim][e].size(); ++j)
           {
-            std::ignore = j;
             dofs[dim][i].push_back(dof_n++);
           }
         }
       }
 
       std::vector<std::size_t> ref;
-      for (int i : dofs[0][1])
-        ref.push_back(i);
-      for (int i : dofs[0][0])
-        ref.push_back(i);
+      ref.insert(ref.end(), dofs[0][1].begin(), dofs[0][1].end());
+      ref.insert(ref.end(), dofs[0][0].begin(), dofs[0][0].end());
       // Edges
-      for (int i : dofs[1][0])
-        ref.push_back(i);
+      ref.insert(ref.end(), dofs[1][0].begin(), dofs[1][0].end());
 
       if (!_dof_transformations_are_identity)
       {
@@ -1208,40 +1204,27 @@ FiniteElement<F>::FiniteElement(
 
       std::vector<std::size_t> rot;
       // Vertices
-      for (int i : dofs[0][1])
-        rot.push_back(i);
-      for (int i : dofs[0][2])
-        rot.push_back(i);
-      for (int i : dofs[0][0])
-        rot.push_back(i);
+      rot.insert(rot.end(), dofs[0][1].begin(), dofs[0][1].end());
+      rot.insert(rot.end(), dofs[0][2].begin(), dofs[0][2].end());
+      rot.insert(rot.end(), dofs[0][0].begin(), dofs[0][0].end());
       // Edges
-      for (int i : dofs[1][1])
-        rot.push_back(i);
-      for (int i : dofs[1][2])
-        rot.push_back(i);
-      for (int i : dofs[1][0])
-        rot.push_back(i);
+      rot.insert(rot.end(), dofs[1][1].begin(), dofs[1][1].end());
+      rot.insert(rot.end(), dofs[1][2].begin(), dofs[1][2].end());
+      rot.insert(rot.end(), dofs[1][0].begin(), dofs[1][0].end());
       // Face
-      for (int i : dofs[2][0])
-        rot.push_back(i);
+      rot.insert(rot.end(), dofs[2][0].begin(), dofs[2][0].end());
 
       std::vector<std::size_t> ref;
-      for (int i : dofs[0][0])
-        ref.push_back(i);
-      for (int i : dofs[0][2])
-        ref.push_back(i);
-      for (int i : dofs[0][1])
-        ref.push_back(i);
+      // Vertices
+      ref.insert(ref.end(), dofs[0][0].begin(), dofs[0][0].end());
+      ref.insert(ref.end(), dofs[0][2].begin(), dofs[0][2].end());
+      ref.insert(ref.end(), dofs[0][1].begin(), dofs[0][1].end());
       // Edges
-      for (int i : dofs[1][0])
-        ref.push_back(i);
-      for (int i : dofs[1][2])
-        ref.push_back(i);
-      for (int i : dofs[1][1])
-        ref.push_back(i);
+      ref.insert(ref.end(), dofs[1][0].begin(), dofs[1][0].end());
+      ref.insert(ref.end(), dofs[1][2].begin(), dofs[1][2].end());
+      ref.insert(ref.end(), dofs[1][1].begin(), dofs[1][1].end());
       // Face
-      for (int i : dofs[2][0])
-        ref.push_back(i);
+      ref.insert(ref.end(), dofs[2][0].begin(), dofs[2][0].end());
 
       if (!_dof_transformations_are_identity)
       {
@@ -1306,48 +1289,30 @@ FiniteElement<F>::FiniteElement(
 
       std::vector<std::size_t> rot;
       // Vertices
-      for (int i : dofs[0][1])
-        rot.push_back(i);
-      for (int i : dofs[0][3])
-        rot.push_back(i);
-      for (int i : dofs[0][0])
-        rot.push_back(i);
-      for (int i : dofs[0][2])
-        rot.push_back(i);
+      rot.insert(rot.end(), dofs[0][1].begin(), dofs[0][1].end());
+      rot.insert(rot.end(), dofs[0][3].begin(), dofs[0][3].end());
+      rot.insert(rot.end(), dofs[0][0].begin(), dofs[0][0].end());
+      rot.insert(rot.end(), dofs[0][2].begin(), dofs[0][2].end());
       // Edges
-      for (int i : dofs[1][2])
-        rot.push_back(i);
-      for (int i : dofs[1][0])
-        rot.push_back(i);
-      for (int i : dofs[1][3])
-        rot.push_back(i);
-      for (int i : dofs[1][1])
-        rot.push_back(i);
+      rot.insert(rot.end(), dofs[1][2].begin(), dofs[1][2].end());
+      rot.insert(rot.end(), dofs[1][0].begin(), dofs[1][0].end());
+      rot.insert(rot.end(), dofs[1][3].begin(), dofs[1][3].end());
+      rot.insert(rot.end(), dofs[1][1].begin(), dofs[1][1].end());
       // Face
-      for (int i : dofs[2][0])
-        rot.push_back(i);
+      rot.insert(rot.end(), dofs[2][0].begin(), dofs[2][0].end());
 
       std::vector<std::size_t> ref;
-      for (int i : dofs[0][0])
-        ref.push_back(i);
-      for (int i : dofs[0][2])
-        ref.push_back(i);
-      for (int i : dofs[0][1])
-        ref.push_back(i);
-      for (int i : dofs[0][3])
-        ref.push_back(i);
+      ref.insert(ref.end(), dofs[0][0].begin(), dofs[0][0].end());
+      ref.insert(ref.end(), dofs[0][2].begin(), dofs[0][2].end());
+      ref.insert(ref.end(), dofs[0][1].begin(), dofs[0][1].end());
+      ref.insert(ref.end(), dofs[0][3].begin(), dofs[0][3].end());
       // Edges
-      for (int i : dofs[1][1])
-        ref.push_back(i);
-      for (int i : dofs[1][0])
-        ref.push_back(i);
-      for (int i : dofs[1][3])
-        ref.push_back(i);
-      for (int i : dofs[1][2])
-        ref.push_back(i);
+      ref.insert(ref.end(), dofs[1][1].begin(), dofs[1][1].end());
+      ref.insert(ref.end(), dofs[1][0].begin(), dofs[1][0].end());
+      ref.insert(ref.end(), dofs[1][3].begin(), dofs[1][3].end());
+      ref.insert(ref.end(), dofs[1][2].begin(), dofs[1][2].end());
       // Face
-      for (int i : dofs[2][0])
-        ref.push_back(i);
+      ref.insert(ref.end(), dofs[2][0].begin(), dofs[2][0].end());
 
       if (!_dof_transformations_are_identity)
       {

From d3495d9732d72dba37d9a2ce5f42728d52e5d1bf Mon Sep 17 00:00:00 2001
From: Matthew Scroggs <matthew.w.scroggs@gmail.com>
Date: Fri, 10 May 2024 13:31:02 +0100
Subject: [PATCH 09/11] add permute_subentity_closure_inv

---
 cpp/basix/finite-element.cpp     | 89 +++++++++++++++++++++++++++++---
 cpp/basix/finite-element.h       | 82 +++++++++++++++++++++++++----
 python/wrapper.cpp               | 24 +++++----
 test/test_dof_transformations.py | 11 ++--
 4 files changed, 175 insertions(+), 31 deletions(-)

diff --git a/cpp/basix/finite-element.cpp b/cpp/basix/finite-element.cpp
index 05b39eec3..604cdc9c6 100644
--- a/cpp/basix/finite-element.cpp
+++ b/cpp/basix/finite-element.cpp
@@ -992,7 +992,7 @@ FiniteElement<F>::FiniteElement(
         for (std::size_t i = 0; i < trans.extent(0); ++i)
         {
           std::vector<std::size_t> perm(trans.extent(1));
-          std::vector<std::size_t> rev_perm(trans.extent(1));
+          std::vector<std::size_t> inv_perm(trans.extent(1));
           for (std::size_t row = 0; row < trans.extent(1); ++row)
           {
             for (std::size_t col = 0; col < trans.extent(1); ++col)
@@ -1000,7 +1000,7 @@ FiniteElement<F>::FiniteElement(
               if (trans(i, row, col) > 0.5)
               {
                 perm[row] = col;
-                rev_perm[col] = row;
+                inv_perm[col] = row;
                 break;
               }
             }
@@ -1008,13 +1008,13 @@ FiniteElement<F>::FiniteElement(
 
           // Factorise the permutations
           precompute::prepare_permutation(perm);
-          precompute::prepare_permutation(rev_perm);
+          precompute::prepare_permutation(inv_perm);
 
           // Store the permutations
           auto& eperm = _eperm.try_emplace(ctype).first->second;
-          auto& eperm_rev = _eperm_rev.try_emplace(ctype).first->second;
+          auto& eperm_inv = _eperm_inv.try_emplace(ctype).first->second;
           eperm.push_back(perm);
-          eperm_rev.push_back(rev_perm);
+          eperm_inv.push_back(inv_perm);
 
           // Generate the entity transformations from the permutations
           std::pair<std::vector<F>, std::array<std::size_t, 2>> identity
@@ -1030,8 +1030,8 @@ FiniteElement<F>::FiniteElement(
           auto& etrans_inv = _etrans_inv.try_emplace(ctype).first->second;
           etrans.push_back({perm, identity});
           etrans_invT.push_back({perm, identity});
-          etransT.push_back({rev_perm, identity});
-          etrans_inv.push_back({rev_perm, identity});
+          etransT.push_back({inv_perm, identity});
+          etrans_inv.push_back({inv_perm, identity});
         }
       }
     }
@@ -1177,6 +1177,10 @@ FiniteElement<F>::FiniteElement(
       auto& secp = _subentity_closure_perm.try_emplace(cell::type::interval)
                        .first->second;
       secp.push_back(ref);
+      auto& secpi
+          = _subentity_closure_perm_inv.try_emplace(cell::type::interval)
+                .first->second;
+      secpi.push_back(ref);
     }
     if (_cell_type == cell::type::tetrahedron || cell_type == cell::type::prism
         || cell_type == cell::type::pyramid)
@@ -1214,6 +1218,18 @@ FiniteElement<F>::FiniteElement(
       // Face
       rot.insert(rot.end(), dofs[2][0].begin(), dofs[2][0].end());
 
+      std::vector<std::size_t> rot_inv;
+      // Vertices
+      rot_inv.insert(rot_inv.end(), dofs[0][2].begin(), dofs[0][2].end());
+      rot_inv.insert(rot_inv.end(), dofs[0][0].begin(), dofs[0][0].end());
+      rot_inv.insert(rot_inv.end(), dofs[0][1].begin(), dofs[0][1].end());
+      // Edges
+      rot_inv.insert(rot_inv.end(), dofs[1][2].begin(), dofs[1][2].end());
+      rot_inv.insert(rot_inv.end(), dofs[1][0].begin(), dofs[1][0].end());
+      rot_inv.insert(rot_inv.end(), dofs[1][1].begin(), dofs[1][1].end());
+      // Face
+      rot_inv.insert(rot_inv.end(), dofs[2][0].begin(), dofs[2][0].end());
+
       std::vector<std::size_t> ref;
       // Vertices
       ref.insert(ref.end(), dofs[0][0].begin(), dofs[0][0].end());
@@ -1230,6 +1246,7 @@ FiniteElement<F>::FiniteElement(
       {
         auto& trans1 = _eperm.at(cell::type::interval)[0];
         auto& trans2 = _eperm.at(cell::type::triangle);
+        auto& trans3 = _eperm_inv.at(cell::type::triangle);
 
         if (!dofs[1][0].empty())
         {
@@ -1253,15 +1270,36 @@ FiniteElement<F>::FiniteElement(
           precompute::apply_permutation(trans2[1], std::span(ref),
                                         dofs[2][0][0]);
         }
+        if (!dofs[1][1].empty())
+        {
+          precompute::apply_permutation(trans1, std::span(rot_inv),
+                                        dofs[1][1][0]);
+        }
+        if (!dofs[1][2].empty())
+        {
+          precompute::apply_permutation(trans1, std::span(rot_inv),
+                                        dofs[1][2][0]);
+        }
+        if (!dofs[2][0].empty())
+        {
+          precompute::apply_permutation(trans3[0], std::span(rot_inv),
+                                        dofs[2][0][0]);
+        }
       }
 
       precompute::prepare_permutation(rot);
+      precompute::prepare_permutation(rot_inv);
       precompute::prepare_permutation(ref);
 
       auto& secp = _subentity_closure_perm.try_emplace(cell::type::triangle)
                        .first->second;
       secp.push_back(rot);
       secp.push_back(ref);
+      auto& secpi
+          = _subentity_closure_perm_inv.try_emplace(cell::type::triangle)
+                .first->second;
+      secpi.push_back(rot_inv);
+      secpi.push_back(ref);
     }
     if (_cell_type == cell::type::hexahedron || cell_type == cell::type::prism
         || cell_type == cell::type::pyramid)
@@ -1301,6 +1339,20 @@ FiniteElement<F>::FiniteElement(
       // Face
       rot.insert(rot.end(), dofs[2][0].begin(), dofs[2][0].end());
 
+      std::vector<std::size_t> rot_inv;
+      // Vertices
+      rot_inv.insert(rot_inv.end(), dofs[0][2].begin(), dofs[0][2].end());
+      rot_inv.insert(rot_inv.end(), dofs[0][0].begin(), dofs[0][0].end());
+      rot_inv.insert(rot_inv.end(), dofs[0][3].begin(), dofs[0][3].end());
+      rot_inv.insert(rot_inv.end(), dofs[0][1].begin(), dofs[0][1].end());
+      // Edges
+      rot_inv.insert(rot_inv.end(), dofs[1][1].begin(), dofs[1][1].end());
+      rot_inv.insert(rot_inv.end(), dofs[1][3].begin(), dofs[1][3].end());
+      rot_inv.insert(rot_inv.end(), dofs[1][0].begin(), dofs[1][0].end());
+      rot_inv.insert(rot_inv.end(), dofs[1][2].begin(), dofs[1][2].end());
+      // Face
+      rot_inv.insert(rot_inv.end(), dofs[2][0].begin(), dofs[2][0].end());
+
       std::vector<std::size_t> ref;
       ref.insert(ref.end(), dofs[0][0].begin(), dofs[0][0].end());
       ref.insert(ref.end(), dofs[0][2].begin(), dofs[0][2].end());
@@ -1318,6 +1370,7 @@ FiniteElement<F>::FiniteElement(
       {
         auto& trans1 = _eperm.at(cell::type::interval)[0];
         auto& trans2 = _eperm.at(cell::type::quadrilateral);
+        auto& trans3 = _eperm_inv.at(cell::type::quadrilateral);
 
         if (!dofs[1][1].empty())
         {
@@ -1338,9 +1391,26 @@ FiniteElement<F>::FiniteElement(
           precompute::apply_permutation(trans2[1], std::span(ref),
                                         dofs[2][0][0]);
         }
+
+        if (!dofs[1][1].empty())
+        {
+          precompute::apply_permutation(trans1, std::span(rot_inv),
+                                        dofs[1][1][0]);
+        }
+        if (!dofs[1][2].empty())
+        {
+          precompute::apply_permutation(trans1, std::span(rot_inv),
+                                        dofs[1][2][0]);
+        }
+        if (!dofs[2][0].empty())
+        {
+          precompute::apply_permutation(trans3[0], std::span(rot_inv),
+                                        dofs[2][0][0]);
+        }
       }
 
       precompute::prepare_permutation(rot);
+      precompute::prepare_permutation(rot_inv);
       precompute::prepare_permutation(ref);
 
       auto& secp
@@ -1348,6 +1418,11 @@ FiniteElement<F>::FiniteElement(
                 .first->second;
       secp.push_back(rot);
       secp.push_back(ref);
+      auto& secpi
+          = _subentity_closure_perm_inv.try_emplace(cell::type::quadrilateral)
+                .first->second;
+      secpi.push_back(rot_inv);
+      secpi.push_back(ref);
     }
   }
 
diff --git a/cpp/basix/finite-element.h b/cpp/basix/finite-element.h
index a8a69958f..dc849ebc5 100644
--- a/cpp/basix/finite-element.h
+++ b/cpp/basix/finite-element.h
@@ -834,7 +834,7 @@ class FiniteElement
     if (_dof_transformations_are_identity)
       return;
     else
-      permute_data<std::int32_t, true>(d, 1, cell_info, _eperm_rev);
+      permute_data<std::int32_t, true>(d, 1, cell_info, _eperm_inv);
   }
 
   /// @brief Permute indices associated with degree-of-freedoms on the
@@ -849,10 +849,10 @@ class FiniteElement
   /// @param[in,out] d Indices associated with each reference element
   /// degree-of-freedom (in). Indices associated with each physical
   /// element degree-of-freedom (out).
-  /// @param cell_info Permutation info for the cell
+  /// @param entity_info Permutation info for the cell
   /// @param entity_type The cell type of the sub-entity
   void permute_subentity_closure(std::span<std::int32_t> d,
-                                 std::uint32_t cell_info,
+                                 std::uint32_t entity_info,
                                  cell::type entity_type) const
   {
     if (!_dof_transformations_are_permutations)
@@ -869,7 +869,7 @@ class FiniteElement
     auto& perm = _subentity_closure_perm.at(entity_type);
     if (entity_dim == 1)
     {
-      if (cell_info & 1)
+      if (entity_info & 1)
       {
         precompute::apply_permutation(perm[0], d);
       }
@@ -877,13 +877,13 @@ class FiniteElement
     else if (entity_dim == 2)
     {
       // Rotate a face
-      for (std::uint32_t r = 0; r < (cell_info >> 1 & 3); ++r)
+      for (std::uint32_t r = 0; r < (entity_info >> 1 & 3); ++r)
       {
         precompute::apply_permutation(perm[0], d);
       }
 
       // Reflect a face (post rotate)
-      if (cell_info & 1)
+      if (entity_info & 1)
       {
         precompute::apply_permutation(perm[1], d);
       }
@@ -895,6 +895,61 @@ class FiniteElement
     }
   }
 
+  /// @brief Perform the inverse of the operation applied by
+  /// permute_subentity_closure().
+  ///
+  /// @note This function is designed to be called at runtime, so its
+  /// performance is critical.
+  ///
+  /// @param[in,out] d Indices associated with each reference element
+  /// degree-of-freedom (in). Indices associated with each physical
+  /// element degree-of-freedom (out).
+  /// @param entity_info Permutation info for the cell
+  /// @param entity_type The cell type of the sub-entity
+  void permute_subentity_closure_inv(std::span<std::int32_t> d,
+                                     std::uint32_t entity_info,
+                                     cell::type entity_type) const
+  {
+    if (!_dof_transformations_are_permutations)
+    {
+      throw std::runtime_error(
+          "The DOF transformations for this element are not permutations");
+    }
+
+    const int entity_dim = cell::topological_dimension(entity_type);
+
+    if (entity_dim == 0)
+      return;
+
+    auto& perm = _subentity_closure_perm_inv.at(entity_type);
+    if (entity_dim == 1)
+    {
+      if (entity_info & 1)
+      {
+        precompute::apply_permutation(perm[0], d);
+      }
+    }
+    else if (entity_dim == 2)
+    {
+      // Reflect a face (pre rotate)
+      if (entity_info & 1)
+      {
+        precompute::apply_permutation(perm[1], d);
+      }
+
+      // Rotate a face
+      for (std::uint32_t r = 0; r < (entity_info >> 1 & 3); ++r)
+      {
+        precompute::apply_permutation(perm[0], d);
+      }
+    }
+    else
+    {
+      throw std::runtime_error(
+          "Invalid dimension for permute_subentity_closure");
+    }
+  }
+
   /// @brief Transform basis functions from the reference element
   /// ordering and orientation to the globally consistent physical
   /// element ordering and orientation.
@@ -1367,7 +1422,7 @@ class FiniteElement
   // The reverse entity permutations (factorised). This will only be set
   // if _dof_transformations_are_permutations is True and
   // _dof_transformations_are_identity is False
-  std::map<cell::type, std::vector<std::vector<std::size_t>>> _eperm_rev;
+  std::map<cell::type, std::vector<std::vector<std::size_t>>> _eperm_inv;
 
   // The entity transformations in precomputed form
   std::map<cell::type, trans_data_t> _etrans;
@@ -1386,6 +1441,11 @@ class FiniteElement
   std::map<cell::type, std::vector<std::vector<std::size_t>>>
       _subentity_closure_perm;
 
+  // The inverse subentity closure permutations (factorised). This will only be
+  // set if _dof_transformations_are_permutations is True
+  std::map<cell::type, std::vector<std::vector<std::size_t>>>
+      _subentity_closure_perm_inv;
+
   // Indicates whether or not this is the discontinuous version of the
   // element
   bool _discontinuous;
@@ -1691,7 +1751,7 @@ void FiniteElement<F>::Tt_apply(std::span<T> u, int n,
   if (_dof_transformations_are_identity)
     return;
   else if (_dof_transformations_are_permutations)
-    permute_data<T, true>(u, n, cell_info, _eperm_rev);
+    permute_data<T, true>(u, n, cell_info, _eperm_inv);
   else
   {
     transform_data<T, true>(u, n, cell_info, _etransT,
@@ -1723,7 +1783,7 @@ void FiniteElement<F>::Tinv_apply(std::span<T> u, int n,
   if (_dof_transformations_are_identity)
     return;
   else if (_dof_transformations_are_permutations)
-    permute_data<T, true>(u, n, cell_info, _eperm_rev);
+    permute_data<T, true>(u, n, cell_info, _eperm_inv);
   else
   {
     transform_data<T, true>(u, n, cell_info, _etrans_inv,
@@ -1793,7 +1853,7 @@ void FiniteElement<F>::T_apply_right(std::span<T> u, int n,
     for (int i = 0; i < n; ++i)
     {
       std::span<T> dblock(u.data() + i * step, step);
-      permute_data<T, true>(dblock, 1, cell_info, _eperm_rev);
+      permute_data<T, true>(dblock, 1, cell_info, _eperm_inv);
     }
   }
   else
@@ -1817,7 +1877,7 @@ void FiniteElement<F>::Tt_inv_apply_right(std::span<T> u, int n,
     for (int i = 0; i < n; ++i)
     {
       std::span<T> dblock(u.data() + i * step, step);
-      permute_data<T, true>(dblock, 1, cell_info, _eperm_rev);
+      permute_data<T, true>(dblock, 1, cell_info, _eperm_inv);
     }
   }
   else
diff --git a/python/wrapper.cpp b/python/wrapper.cpp
index 118be40f9..9f7f50687 100644
--- a/python/wrapper.cpp
+++ b/python/wrapper.cpp
@@ -101,10 +101,19 @@ void declare_float(nb::module_& m, std::string type)
       .def("permute_subentity_closure",
            [](const FiniteElement<T>& self,
               nb::ndarray<std::int32_t, nb::ndim<1>, nb::c_contig> d,
-              std::uint32_t cell_info, cell::type entity_type)
+              std::uint32_t entity_info, cell::type entity_type)
            {
              std::span<std::int32_t> _d(d.data(), d.shape(0));
-             self.permute_subentity_closure(_d, cell_info, entity_type);
+             self.permute_subentity_closure(_d, entity_info, entity_type);
+             return as_nbarray(std::move(_d));
+           })
+      .def("permute_subentity_closure_inv",
+           [](const FiniteElement<T>& self,
+              nb::ndarray<std::int32_t, nb::ndim<1>, nb::c_contig> d,
+              std::uint32_t entity_info, cell::type entity_type)
+           {
+             std::span<std::int32_t> _d(d.data(), d.shape(0));
+             self.permute_subentity_closure_inv(_d, entity_info, entity_type);
              return as_nbarray(std::move(_d));
            })
       .def("push_forward",
@@ -152,13 +161,10 @@ void declare_float(nb::module_& m, std::string type)
              self.Tt_apply_right(std::span(u.data(), u.size()), n,
                                 cell_info);
            })
-      .def("Tt_inv_apply",
-           [](const FiniteElement<T>& self,
-              nb::ndarray<T, nb::ndim<1>, nb::c_contig> u, int n,
-              std::uint32_t cell_info) {
-             self.Tt_inv_apply(std::span(u.data(), u.size()), n,
-                               cell_info);
-           })
+      .def("Tt_inv_apply", [](const FiniteElement<T>& self,
+                              nb::ndarray<T, nb::ndim<1>, nb::c_contig> u,
+                              int n, std::uint32_t cell_info)
+           { self.Tt_inv_apply(std::span(u.data(), u.size()), n, cell_info); })
       .def("base_transformations", [](const FiniteElement<T>& self)
            { return as_nbarrayp(self.base_transformations()); })
       .def("entity_transformations",
diff --git a/test/test_dof_transformations.py b/test/test_dof_transformations.py
index 1d722544b..677bd53a2 100644
--- a/test/test_dof_transformations.py
+++ b/test/test_dof_transformations.py
@@ -582,8 +582,11 @@ def test_transformation_of_tabulated_data_pyramid(element_type, degree, element_
 def test_permute_subentity_closure(family, cell_type, degree, args, subentity, results):
     e = basix.create_element(family, cell_type, degree, *args)
 
-    for cell_info, result in results.items():
-        data = np.array(list(range(len(result))))
-        data = list(e._e.permute_subentity_closure(data, cell_info, subentity.value))
-
+    for entity_info, result in results.items():
+        data = np.arange(len(result))
+        data = list(e._e.permute_subentity_closure(data, entity_info, subentity.value))
         assert result == data
+
+        data = np.array(result)
+        data = list(e._e.permute_subentity_closure_inv(data, entity_info, subentity.value))
+        assert data == list(range(len(result)))

From bc48978293ceee09b9bda50c9c9b4f117a6bfa37 Mon Sep 17 00:00:00 2001
From: Matthew Scroggs <matthew.w.scroggs@gmail.com>
Date: Fri, 10 May 2024 13:43:46 +0100
Subject: [PATCH 10/11] overloaded versions that take cell info and entity
 index

---
 cpp/basix/finite-element.h | 83 +++++++++++++++++++++++++++++++++++++-
 1 file changed, 81 insertions(+), 2 deletions(-)

diff --git a/cpp/basix/finite-element.h b/cpp/basix/finite-element.h
index dc849ebc5..7992ce3e6 100644
--- a/cpp/basix/finite-element.h
+++ b/cpp/basix/finite-element.h
@@ -849,8 +849,87 @@ class FiniteElement
   /// @param[in,out] d Indices associated with each reference element
   /// degree-of-freedom (in). Indices associated with each physical
   /// element degree-of-freedom (out).
-  /// @param entity_info Permutation info for the cell
+  /// @param cell_info Permutation info for the cell
+  /// @param entity_type The cell type of the sub-entity
+  /// @param entity_index The index of the entity
+  void permute_subentity_closure(std::span<std::int32_t> d,
+                                 std::uint32_t cell_info,
+                                 cell::type entity_type, int entity_index) const
+  {
+    const int entity_dim = cell::topological_dimension(entity_type);
+
+    int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
+
+    std::uint32_t entity_info = 0;
+    switch (entity_dim)
+    {
+    case 0:
+      return;
+    case 1:
+      entity_info = cell_info >> (face_start + entity_index) & 1;
+      break;
+    case 2:
+      entity_info = cell_info >> (3 * entity_index) & 7;
+      break;
+    default:
+      throw std::runtime_error("Unsupported cell dimension");
+    }
+    permute_subentity_closure(d, entity_info, entity_type);
+  }
+
+  /// @brief Perform the inverse of the operation applied by
+  /// permute_subentity_closure().
+  ///
+  /// @note This function is designed to be called at runtime, so its
+  /// performance is critical.
+  ///
+  /// @param[in,out] d Indices associated with each reference element
+  /// degree-of-freedom (in). Indices associated with each physical
+  /// element degree-of-freedom (out).
+  /// @param cell_info Permutation info for the cell
+  /// @param entity_type The cell type of the sub-entity
+  /// @param entity_index The index of the entity
+  void permute_subentity_closure_inv(std::span<std::int32_t> d,
+                                     std::uint32_t cell_info,
+                                     cell::type entity_type,
+                                     int entity_index) const
+  {
+    const int entity_dim = cell::topological_dimension(entity_type);
+
+    int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
+
+    std::uint32_t entity_info;
+    switch (entity_dim)
+    {
+    case 0:
+      return;
+    case 1:
+      entity_info = cell_info >> (face_start + entity_index) & 1;
+      break;
+    case 2:
+      entity_info = cell_info >> (3 * entity_index) & 7;
+      break;
+    default:
+      throw std::runtime_error("Unsupported cell dimension");
+    }
+    permute_subentity_closure(d, entity_info, entity_type);
+  }
+
+  /// @brief Permute indices associated with degree-of-freedoms on the
+  /// closure of a sub-entity of the reference element
+  ///
+  /// This function performs a similar permutation to permute() but
+  /// additionally permutes the positions of vertices and edges
+  ///
+  /// @note This function is designed to be called at runtime, so its
+  /// performance is critical.
+  ///
+  /// @param[in,out] d Indices associated with each reference element
+  /// degree-of-freedom (in). Indices associated with each physical
+  /// element degree-of-freedom (out).
+  /// @param entity_info Permutation info for the entity
   /// @param entity_type The cell type of the sub-entity
+
   void permute_subentity_closure(std::span<std::int32_t> d,
                                  std::uint32_t entity_info,
                                  cell::type entity_type) const
@@ -904,7 +983,7 @@ class FiniteElement
   /// @param[in,out] d Indices associated with each reference element
   /// degree-of-freedom (in). Indices associated with each physical
   /// element degree-of-freedom (out).
-  /// @param entity_info Permutation info for the cell
+  /// @param entity_info Permutation info for the entity
   /// @param entity_type The cell type of the sub-entity
   void permute_subentity_closure_inv(std::span<std::int32_t> d,
                                      std::uint32_t entity_info,

From bdb83a20a14a078d3e87c6c5f7640d336752cc95 Mon Sep 17 00:00:00 2001
From: Matthew Scroggs <matthew.w.scroggs@gmail.com>
Date: Fri, 10 May 2024 14:38:55 +0100
Subject: [PATCH 11/11] _inv

---
 cpp/basix/finite-element.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/cpp/basix/finite-element.h b/cpp/basix/finite-element.h
index 7992ce3e6..fc7188198 100644
--- a/cpp/basix/finite-element.h
+++ b/cpp/basix/finite-element.h
@@ -912,7 +912,7 @@ class FiniteElement
     default:
       throw std::runtime_error("Unsupported cell dimension");
     }
-    permute_subentity_closure(d, entity_info, entity_type);
+    permute_subentity_closure_inv(d, entity_info, entity_type);
   }
 
   /// @brief Permute indices associated with degree-of-freedoms on the