Mixed topology: find dS cells (#3632)

* compute mixed cell pairs

* Add test

* Update docs

cpp/dolfinx/mesh/Topology.cpp

@@ -1411,3 +1411,85 @@ mesh::entities_to_index(const Topology& topology, int dim,
   return indices;
 }
 //-----------------------------------------------------------------------------
+std::vector<std::vector<std::int32_t>>
+mesh::compute_mixed_cell_pairs(const Topology& topology,
+                               mesh::CellType facet_type)
+{
+  int tdim = topology.dim();
+  std::vector<mesh::CellType> cell_types = topology.entity_types(tdim);
+  std::vector<mesh::CellType> facet_types = topology.entity_types(tdim - 1);
+
+  int facet_index = -1;
+  for (std::size_t i = 0; i < facet_types.size(); ++i)
+  {
+    if (facet_types[i] == facet_type)
+    {
+      facet_index = i;
+      break;
+    }
+  }
+  if (facet_index == -1)
+    throw std::runtime_error("Cannot find facet type in topology");
+
+  std::vector<std::vector<std::int32_t>> facet_pair_lists;
+  for (std::size_t i = 0; i < cell_types.size(); ++i)
+    for (std::size_t j = 0; j < cell_types.size(); ++j)
+    {
+      std::vector<std::int32_t> facet_pairs_ij;
+      auto fci = topology.connectivity({tdim - 1, facet_index},
+                                       {tdim, static_cast<int>(i)});
+      auto cfi = topology.connectivity({tdim, static_cast<int>(i)},
+                                       {tdim - 1, facet_index});
+
+      auto local_facet = [](auto cf, std::int32_t c, std::int32_t f)
+      {
+        auto it = std::find(cf->links(c).begin(), cf->links(c).end(), f);
+        if (it == cf->links(c).end())
+          throw std::runtime_error("Bad connectivity");
+        return std::distance(cf->links(c).begin(), it);
+      };
+
+      if (i == j)
+      {
+        if (fci)
+        {
+          for (std::int32_t k = 0; k < fci->num_nodes(); ++k)
+          {
+            if (fci->num_links(k) == 2)
+            {
+              std::int32_t c0 = fci->links(k)[0], c1 = fci->links(k)[1];
+              facet_pairs_ij.push_back(c0);
+              facet_pairs_ij.push_back(local_facet(cfi, c0, k));
+              facet_pairs_ij.push_back(c1);
+              facet_pairs_ij.push_back(local_facet(cfi, c1, k));
+            }
+          }
+        }
+      }
+      else
+      {
+        auto fcj = topology.connectivity({tdim - 1, facet_index},
+                                         {tdim, static_cast<int>(j)});
+        auto cfj = topology.connectivity({tdim, static_cast<int>(j)},
+                                         {tdim - 1, facet_index});
+        if (fci and fcj)
+        {
+          assert(fci->num_nodes() == fcj->num_nodes());
+          for (std::int32_t k = 0; k < fci->num_nodes(); ++k)
+          {
+            if (fci->num_links(k) == 1 and fcj->num_links(k) == 1)
+            {
+              std::int32_t ci = fci->links(k)[0], cj = fcj->links(k)[0];
+              facet_pairs_ij.push_back(ci);
+              facet_pairs_ij.push_back(local_facet(cfi, ci, k));
+              facet_pairs_ij.push_back(cj);
+              facet_pairs_ij.push_back(local_facet(cfj, cj, k));
+            }
+          }
+        }
+      }
+      facet_pair_lists.push_back(facet_pairs_ij);
+    }
+
+  return facet_pair_lists;
+}

cpp/dolfinx/mesh/Topology.h

@@ -341,4 +341,25 @@ create_subtopology(const Topology& topology, int dim,
 std::vector<std::int32_t>
 entities_to_index(const Topology& topology, int dim,
                   std::span<const std::int32_t> entities);
+
+/// @brief Compute a list of cell-cell connections for each possible combination
+/// in the topology which have the same connecting facet type.
+/// @param topology A mesh topology
+/// @param facet_type Type of facet connection between cells
+/// @return A list for each possible cell-cell connection, arranged
+/// in the ordering given by `cell_types()` in the topology. i.e. if the cell
+/// types are tet, prism, hex, and the facet_type is quadrilateral, then the
+/// lists returned are: tet-tet (empty), tet-prism (empty), tet-hex (empty),
+/// prism-tet (empty), prism-prism, prism-hex, hex-tet (empty), hex-prism,
+/// hex-hex. Note there are empty lists for the invalid cases, and also that
+/// there is currently redundant data, since the transpose is also computed,
+/// i.e. prism-hex as well as hex-prism.
+/// Each list contains a flattened array with data in the order (cell0, local
+/// facet0, cell1, local facet1) for each connection between cells.
+/// @note All facet-cell connectivity and cell-facet connectivity must be
+/// computed beforehand in the topology.
+/// @todo Remove redundant data.
+std::vector<std::vector<std::int32_t>>
+compute_mixed_cell_pairs(const Topology& topology, mesh::CellType facet_type);
+
 } // namespace dolfinx::mesh

python/dolfinx/wrappers/mesh.cpp

@@ -696,6 +696,8 @@ void mesh(nb::module_& m)
               ghost_owners_span, boundary_vertices_span);
         });
 
+  m.def("compute_mixed_cell_pairs", &dolfinx::mesh::compute_mixed_cell_pairs);
+
   declare_meshtags<std::int8_t>(m, "int8");
   declare_meshtags<std::int32_t>(m, "int32");
   declare_meshtags<std::int64_t>(m, "int64");

python/test/unit/mesh/test_mixed_topology.py

@@ -3,10 +3,15 @@
 import numpy as np
 
 from dolfinx.cpp.log import set_thread_name
-from dolfinx.cpp.mesh import Mesh_float64, create_geometry, create_topology
+from dolfinx.cpp.mesh import (
+    Mesh_float64,
+    compute_mixed_cell_pairs,
+    create_geometry,
+    create_topology,
+)
 from dolfinx.fem import coordinate_element
 from dolfinx.log import LogLevel, set_log_level
-from dolfinx.mesh import CellType, GhostMode, create_unit_cube
+from dolfinx.mesh import CellType, GhostMode, Mesh, create_unit_cube
 
 
 def test_mixed_topology_mesh():
@@ -271,3 +276,27 @@ def test_create_entities():
 
     # Triangle and quad to prism (facet->cell)
     mesh.topology.create_connectivity(2, 3)
+
+
+def test_mixed_cell_pairs(mixed_topology_mesh):
+    mesh = Mesh(mixed_topology_mesh, None)
+    mesh.topology.create_entities(2)
+    mesh.topology.create_connectivity(2, 3)
+    cell_types = mesh.topology.entity_types[3]
+    facet_types = mesh.topology.entity_types[2]
+    print(cell_types, facet_types)
+
+    # For each facet type
+    for f, ft in enumerate(facet_types):
+        cell_pairs = compute_mixed_cell_pairs(mesh.topology._cpp_object, ft)
+        for i, cti in enumerate(cell_types):
+            for j, ctj in enumerate(cell_types):
+                idx = i * len(cell_types) + j
+                num_conns = len(cell_pairs[idx]) // 4
+                print(f"Connectivity ({ft}) from {cti} to {ctj} : {num_conns}")
+                if len(cell_pairs[idx]) > 0:
+                    connection = np.array(cell_pairs[idx]).reshape((num_conns, -1))
+                    f0 = mesh.topology.connectivity((3, i), (2, f))
+                    f1 = mesh.topology.connectivity((3, j), (2, f))
+                    for row in connection:
+                        assert f0.links(row[0])[row[1]] == f1.links(row[2])[row[3]]