Add ridge permutations

cpp/dolfinx/mesh/Topology.cpp

@@ -886,6 +886,20 @@ const std::vector<std::uint8_t>& Topology::get_facet_permutations() const
   return _facet_permutations;
 }
 //-----------------------------------------------------------------------------
+const std::vector<std::uint8_t>& Topology::get_ridge_permutations() const
+{
+  if (auto i_map = this->index_map(this->dim() - 2);
+      !i_map
+      or (_ridge_permutations.empty()
+          and i_map->size_local() + i_map->num_ghosts() > 0))
+  {
+    throw std::runtime_error(
+        "create_entity_permutations must be called before using this data.");
+  }
+
+  return _ridge_permutations;
+}
+//-----------------------------------------------------------------------------
 const std::vector<std::int32_t>& Topology::interprocess_facets(int index) const
 {
   if (_interprocess_facets.empty())
@@ -1009,9 +1023,10 @@ void Topology::create_entity_permutations()
   for (int d = 0; d < tdim; ++d)
     create_entities(d);
 
-  auto [facet_permutations, cell_permutations]
+  auto [ridge_permutations, facet_permutations, cell_permutations]
       = compute_entity_permutations(*this);
   _facet_permutations = std::move(facet_permutations);
+  _ridge_permutations = std::move(ridge_permutations);
   _cell_permutations = std::move(cell_permutations);
 }
 //-----------------------------------------------------------------------------

cpp/dolfinx/mesh/Topology.h

@@ -159,6 +159,21 @@ class Topology
   /// computed
   const std::vector<std::uint8_t>& get_facet_permutations() const;
 
+  /// @brief Get the numbers that encode the number of permutations to
+  /// apply to ridges.
+  ///
+  /// The permutations are encoded so that:
+  ///
+  /// @todo Matthew needs to document this.
+  ///
+  /// The data is stored in a flattened 2D array, so that `data[cell_index *
+  /// ridges_per_cell + ridge_index]` contains the ridge with index
+  /// `ridge_index` of the cell with index `cell_index`.
+  /// @return The encoded permutation info
+  /// @note An exception is raised if the permutations have not been
+  /// computed
+  const std::vector<std::uint8_t>& get_ridge_permutations() const;
+
   /// @brief Get the types of cells in the topology
   /// @return The cell types
   std::vector<CellType> cell_types() const;
@@ -236,6 +251,11 @@ class Topology
   // celln_0, celln_1, ,celln_2,]
   std::vector<std::uint8_t> _facet_permutations;
 
+  // The ridge permutations (local ridge, cell)
+  // [cell0_0, cell0_1, ,cell0_2, cell1_0, cell1_1, ,cell1_2, ...,
+  // celln_0, celln_1, ,celln_2,]
+  std::vector<std::uint8_t> _ridge_permutations;
+
   // Cell permutation info. See the documentation for
   // get_cell_permutation_info for documentation of how this is encoded.
   std::vector<std::uint32_t> _cell_permutations;

cpp/dolfinx/mesh/permutationcomputation.cpp

@@ -284,17 +284,22 @@ compute_face_permutations(const mesh::Topology& topology)
 } // namespace
 
 //-----------------------------------------------------------------------------
-std::pair<std::vector<std::uint8_t>, std::vector<std::uint32_t>>
+std::tuple<std::vector<std::uint8_t>, std::vector<std::uint8_t>,
+           std::vector<std::uint32_t>>
 mesh::compute_entity_permutations(const mesh::Topology& topology)
 {
+
   common::Timer t_perm("Compute entity permutations");
   const int tdim = topology.dim();
   CellType cell_type = topology.cell_type();
   const std::int32_t num_cells = topology.connectivity(tdim, 0)->num_nodes();
   const int facets_per_cell = cell_num_entities(cell_type, tdim - 1);
+  const int ridges_per_cell = cell_num_entities(cell_type, tdim - 2);
 
   std::vector<std::uint32_t> cell_permutation_info(num_cells, 0);
   std::vector<std::uint8_t> facet_permutations(num_cells * facets_per_cell);
+  std::vector<std::uint8_t> ridge_permutations(num_cells * ridges_per_cell);
+
   std::int32_t used_bits = 0;
   if (tdim > 2)
   {
@@ -335,9 +340,18 @@ mesh::compute_entity_permutations(const mesh::Topology& topology)
           facet_permutations[c * facets_per_cell + i] = edge_perm[c][i];
       }
     }
+    if (tdim == 3)
+    {
+      for (int c = 0; c < num_cells; ++c)
+      {
+        for (int i = 0; i < ridges_per_cell; ++i)
+          ridge_permutations[c * ridges_per_cell + i] = edge_perm[c][i];
+      }
+    }
   }
   assert(used_bits < _BITSETSIZE);
 
-  return {std::move(facet_permutations), std::move(cell_permutation_info)};
+  return {std::move(ridge_permutations), std::move(facet_permutations),
+          std::move(cell_permutation_info)};
 }
 //-----------------------------------------------------------------------------

cpp/dolfinx/mesh/permutationcomputation.h

@@ -66,8 +66,9 @@ class Topology;
 ///    This data is used to correct the direction of vector function
 ///    on permuted facets.
 ///
-/// @return Facet permutation and cells permutations
-std::pair<std::vector<std::uint8_t>, std::vector<std::uint32_t>>
+/// @param[in] topology The mesh topology.
+std::tuple<std::vector<std::uint8_t>, std::vector<std::uint8_t>,
+           std::vector<std::uint32_t>>
 compute_entity_permutations(const Topology& topology);
 
 } // namespace dolfinx::mesh

python/dolfinx/mesh.py

@@ -173,6 +173,21 @@ def get_facet_permutations(self) -> npt.NDArray[np.uint8]:
         """
         return self._cpp_object.get_facet_permutations()
 
+    def get_ridge_permutations(self) -> npt.NDArray[np.uint8]:
+        """Get the permutation integer to apply to ridges.
+
+        The bits of each integer describes the number of reflections
+        that has to be applied to a ridge to map between a
+        ridge in the mesh (relative to a cell) and the corresponding
+        ridge on the reference element. The data has the shape
+        ``(num_cells, num_ridges_per_cell)``, flattened row-wise. The
+        number of cells include potential ghost cells.
+
+        Note:
+            The data can be unpacked with ``numpy.unpack_bits``.
+        """
+        return self._cpp_object.get_ridge_permutations()
+
     def index_map(self, dim: int) -> _cpp.common.IndexMap:
         """Get the IndexMap that describes the parallel distribution of the
         mesh entities.

python/dolfinx/wrappers/mesh.cpp

@@ -668,6 +668,15 @@ void mesh(nb::module_& m)
                                                               {p.size()});
           },
           nb::rv_policy::reference_internal)
+      .def(
+          "get_ridge_permutations",
+          [](const dolfinx::mesh::Topology& self)
+          {
+            const std::vector<std::uint8_t>& p = self.get_ridge_permutations();
+            return nb::ndarray<const std::uint8_t, nb::numpy>(p.data(),
+                                                              {p.size()});
+          },
+          nb::rv_policy::reference_internal)
       .def(
           "get_cell_permutation_info",
           [](const dolfinx::mesh::Topology& self)