Create topology Python wrapper

python/dolfinx/fem/__init__.py

@@ -155,7 +155,7 @@ def compute_integration_domains(
     Returns:
         List of integration entities
     """
-    return _compute_integration_domains(integral_type, topology, entities, dim)
+    return _compute_integration_domains(integral_type, topology._cpp_object, entities, dim)
 
 
 __all__ = [

python/dolfinx/fem/function.py

@@ -647,7 +647,7 @@ def functionspace(
 
     cpp_element = _create_dolfinx_element(mesh.comm, mesh.topology.cell_type, ufl_e, dtype)
 
-    cpp_dofmap = _cpp.fem.create_dofmap(mesh.comm, mesh.topology, cpp_element)
+    cpp_dofmap = _cpp.fem.create_dofmap(mesh.comm, mesh.topology._cpp_object, cpp_element)
 
     assert np.issubdtype(
         mesh.geometry.x.dtype, cpp_element.dtype

python/dolfinx/mesh.py

@@ -27,7 +27,6 @@
     build_dual_graph,
     cell_dim,
     create_cell_partitioner,
-    exterior_facet_indices,
     to_string,
     to_type,
 )
@@ -68,6 +67,154 @@
 ]
 
 
+class Topology:
+    """Topology for a :class:`dolfinx.mesh.Mesh`"""
+
+    _cpp_object: _cpp.mesh.Topology
+
+    def __init__(self, topology: _cpp.mesh.Topology):
+        """Initialize a topology from a C++ topology.
+        Args: The C++ topology object
+        Note:
+            Topology objects should usually be constructed with the
+            :func:`dolfinx.cpp.mesh.create_topology` and not this class initializer.
+        """
+        self._cpp_object = topology
+
+    def cell_name(self) -> str:
+        """String representation of the cell-type of the topology"""
+        return to_string(self._cpp_object.cell_type)
+
+    def connectivity(self, d0: int, d1: int) -> _cpp.graph.AdjacencyList_int32:
+        """Return connectivity from entities of dimension ``d0`` to entities of dimension ``d1``.
+
+        Note:
+            Assumes only one entity type per dimension.
+
+        Args:
+            d0: Dimension of entity one is mapping from
+            d1: Dimension of entity one is mapping to
+        """
+        if (conn := self._cpp_object.connectivity(d0, d1)) is not None:
+            return conn
+        else:
+            raise RuntimeError(
+                f"Connectivity between dimension {d0} and {d1} has not been computed.",
+                f"Please call `dolfinx.mesh.Topology.create_connectivity({d0}, {d1})` first.",
+            )
+
+    @property
+    def comm(self):
+        return self._cpp_object.comm
+
+    def create_connectivity(self, d0: int, d1: int):
+        """Create connectivity between given pair of dimensions, ``d0`` and  ``d1``.
+
+        Args:
+            d0: Dimension of entities one is mapping from
+            d1: Dimension of entities one is mapping to
+        """
+        self._cpp_object.create_connectivity(d0, d1)
+
+    def create_entities(self, dim: int) -> int:
+        """Create entities of given topological dimension.
+
+        Args:
+            dim: Topological dimension
+
+        Returns:
+            Number of newly created entities, returns -1 if entities already existed
+        """
+        return self._cpp_object.create_entities(dim)
+
+    def create_entity_permutations(self):
+        """Compute entity permutations and reflections."""
+        self._cpp_object.create_entity_permutations()
+
+    @property
+    def dim(self) -> int:
+        """Return the topological dimension of the mesh."""
+        return self._cpp_object.dim
+
+    @property
+    def entity_types(self) -> list[list[CellType]]:
+        """Get the entity types in the topology for all topological dimensions."""
+        return self._cpp_object.entity_types
+
+    def get_cell_permutation_info(self) -> npt.NDArray[np.uint32]:
+        """Returns the permutation information
+
+        The returned data is used for packing coefficients and assembling of tensors.
+        The bits of each integer encodes the number of reflections and permutations
+        for each sub-entity of the cell to be able to map it to the reference element.
+        """
+        return self._cpp_object.get_cell_permutation_info()
+
+    def get_facet_permutations(self) -> npt.NDArray[np.uint8]:
+        """Get the permutation number to apply to facets.
+
+        The bits of each integer describes the number of reflections and rotations that has to
+        be applied to a facet to map between a facet in the mesh (relative to a cell)
+        and the corresponding facet on the reference element. The data has
+        the shape ``(num_cells, num_facets_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_facet_permutations()
+
+    def index_map(self, dim: int) -> _cpp.common.IndexMap:
+        """Get the IndexMap that described the parallel distribution of the mesh entities.
+
+        Args:
+            dim: Topological dimension
+
+        Returns:
+            Index map for the entities of dimension ``dim``.
+        """
+        if (imap := self._cpp_object.index_map(dim)) is not None:
+            return imap
+        else:
+            raise RuntimeError(
+                f"Entities of dimension {dim} has not been computed."
+                f"Call `dolfinx.mesh.Topology.create_entities({dim}) first."
+            )
+
+    def interprocess_facets(self) -> npt.NDArray[np.int32]:
+        """List of inter-process facets, if facet topology has been computed."""
+        return self._cpp_object.interprocess_facets()
+
+    @property
+    def original_cell_index(self) -> npt.NDArray[np.int64]:
+        """Get the original cell index"""
+        return self._cpp_object.original_cell_index
+
+    def set_connectivity(self, graph: _cpp.graph.AdjacencyList_int32, d0: int, d1: int):
+        """Set connectivity for given pair of topological dimensions.
+
+        Args:
+            graph: Connectivity graph
+            d0: Topological dimension mapping from
+            d1: Topological dimension mapping to
+        """
+        self._cpp_object.set_connectivity(graph, d0, d1)
+
+    def set_index_map(self, dim: int, index_map: _cpp.common.IndexMap):
+        """Set the IndexMap for dimension ``dim``.
+
+        Args:
+            dim: Topological dimension of entity
+            index_map: IndexMap
+        """
+        return self._cpp_object.set_index_map(dim, index_map)
+
+    @property
+    def cell_type(self) -> CellType:
+        """Get the cell type of the topology"""
+        return self._cpp_object.cell_type
+
+
 class Geometry:
     """The geometry of a :class:`dolfinx.mesh.Mesh`"""
 
@@ -123,6 +270,7 @@ class Mesh:
     """A mesh."""
 
     _mesh: typing.Union[_cpp.mesh.Mesh_float32, _cpp.mesh.Mesh_float64]
+    _topology: Topology
     _geometry: Geometry
     _ufl_domain: typing.Optional[ufl.Mesh]
 
@@ -139,6 +287,7 @@ def __init__(self, mesh, domain: typing.Optional[ufl.Mesh]):
             that depend on the scalar type used in the Mesh.
         """
         self._cpp_object = mesh
+        self._topology = Topology(self._cpp_object.topology)
         self._geometry = Geometry(self._cpp_object.geometry)
         self._ufl_domain = domain
         if self._ufl_domain is not None:
@@ -193,9 +342,9 @@ def h(self, dim: int, entities: npt.NDArray[np.int32]) -> npt.NDArray[np.float64
         return _cpp.mesh.h(self._cpp_object, dim, entities)
 
     @property
-    def topology(self):
+    def topology(self) -> Topology:
         "Mesh topology."
-        return self._cpp_object.topology
+        return self._topology
 
     @property
     def geometry(self) -> Geometry:
@@ -269,8 +418,19 @@ def find(self, value) -> npt.NDArray[np.int32]:
         return self._cpp_object.find(value)
 
 
-def compute_incident_entities(topology, entities: npt.NDArray[np.int32], d0: int, d1: int):
-    return _cpp.mesh.compute_incident_entities(topology, entities, d0, d1)
+def compute_incident_entities(
+    topology: Topology, entities: npt.NDArray[np.int32], d0: int, d1: int
+) -> npt.NDArray[np.int32]:
+    """Compute all entities of ``d1`` connected to ``entities`` of dimension ``d0``.
+
+    Args:
+        topology: The topology
+        entities: List of entities fo dimension ``d0``.
+        d0: Topological dimension
+        d1: Topological dimension to map to
+
+    """
+    return _cpp.mesh.compute_incident_entities(topology._cpp_object, entities, d0, d1)
 
 
 def compute_midpoints(mesh: Mesh, dim: int, entities: npt.NDArray[np.int32]):
@@ -321,19 +481,6 @@ def locate_entities_boundary(mesh: Mesh, dim: int, marker: typing.Callable) -> n
     return _cpp.mesh.locate_entities_boundary(mesh._cpp_object, dim, marker)
 
 
-_uflcell_to_dolfinxcell = {
-    "interval": CellType.interval,
-    "interval2D": CellType.interval,
-    "interval3D": CellType.interval,
-    "triangle": CellType.triangle,
-    "triangle3D": CellType.triangle,
-    "quadrilateral": CellType.quadrilateral,
-    "quadrilateral3D": CellType.quadrilateral,
-    "tetrahedron": CellType.tetrahedron,
-    "hexahedron": CellType.hexahedron,
-}
-
-
 def transfer_meshtag(
     meshtag: MeshTags,
     mesh1: Mesh,
@@ -355,12 +502,14 @@ def transfer_meshtag(
         Mesh tags on the refined mesh.
     """
     if meshtag.dim == meshtag.topology.dim:
-        mt = _cpp.refinement.transfer_cell_meshtag(meshtag._cpp_object, mesh1.topology, parent_cell)
+        mt = _cpp.refinement.transfer_cell_meshtag(
+            meshtag._cpp_object, mesh1.topology._cpp_object, parent_cell
+        )
         return MeshTags(mt)
     elif meshtag.dim == meshtag.topology.dim - 1:
         assert parent_facet is not None
         mt = _cpp.refinement.transfer_facet_meshtag(
-            meshtag._cpp_object, mesh1.topology, parent_cell, parent_facet
+            meshtag._cpp_object, mesh1.topology._cpp_object, parent_cell, parent_facet
         )
         return MeshTags(mt)
     else:
@@ -489,7 +638,7 @@ def create_submesh(
     submsh, entity_map, vertex_map, geom_map = _cpp.mesh.create_submesh(
         msh._cpp_object, dim, entities
     )
-    submsh_ufl_cell = ufl.Cell(submsh.topology.cell_name())
+    submsh_ufl_cell = ufl.Cell(to_string(submsh.topology.cell_type))
     submsh_domain = ufl.Mesh(
         basix.ufl.element(
             "Lagrange",
@@ -546,7 +695,9 @@ def meshtags(
     else:
         raise NotImplementedError(f"Type {values.dtype} not supported.")
 
-    return MeshTags(ftype(mesh.topology, dim, np.asarray(entities, dtype=np.int32), values))
+    return MeshTags(
+        ftype(mesh.topology._cpp_object, dim, np.asarray(entities, dtype=np.int32), values)
+    )
 
 
 def meshtags_from_entities(
@@ -577,7 +728,7 @@ def meshtags_from_entities(
     elif isinstance(values, float):
         values = np.full(entities.num_nodes, values, dtype=np.double)
     values = np.asarray(values)
-    return MeshTags(_cpp.mesh.create_meshtags(mesh.topology, dim, entities, values))
+    return MeshTags(_cpp.mesh.create_meshtags(mesh.topology._cpp_object, dim, entities, values))
 
 
 def create_interval(
@@ -821,6 +972,24 @@ def entities_to_geometry(
     return _cpp.mesh.entities_to_geometry(mesh._cpp_object, dim, entities, permute)
 
 
+def exterior_facet_indices(topology: Topology) -> npt.NDArray[np.int32]:
+    """Compute the indices of all exterior facets that are owned by the caller.
+
+    An exterior facet (co-dimension 1) is one that is connected globally to
+    only one cell of co-dimension 0).
+
+    Note:
+        This is a collective operation that should be called on all processes.
+
+    Args:
+        topology: The topology
+
+    Returns:
+        Sorted list of owned facet indices that are exterior facets of the mesh.
+    """
+    return _cpp.mesh.exterior_facet_indices(topology._cpp_object)
+
+
 def create_geometry(
     index_map: _IndexMap,
     dofmap: npt.NDArray[np.int32],

python/dolfinx/wrappers/mesh.cpp

@@ -685,8 +685,6 @@ void mesh(nb::module_& m)
               entity_types.push_back(self.entity_types(i));
             return entity_types;
           })
-      .def("cell_name", [](const dolfinx::mesh::Topology& self)
-           { return dolfinx::mesh::to_string(self.cell_type()); })
       .def(
           "interprocess_facets",
           [](const dolfinx::mesh::Topology& self)

python/test/unit/fem/test_assemble_submesh.py

@@ -13,7 +13,7 @@
 
 import ufl
 from dolfinx import default_scalar_type, fem, la
-from dolfinx.cpp.fem import compute_integration_domains
+from dolfinx.fem import compute_integration_domains
 from dolfinx.mesh import (
     GhostMode,
     compute_incident_entities,

python/test/unit/mesh/test_mixed_topology.py

@@ -242,7 +242,7 @@ def test_create_entities():
 
     quad_v = mesh.topology.connectivity((2, qi), (0, 0))
     tri_v = mesh.topology.connectivity((2, ti), (0, 0))
-    ims = mesh.topology.index_maps(2)
+    ims = mesh.topology._cpp_object.index_maps(2)
     assert len(ims) == 2
     assert ims[qi].size_global == 32
     assert len(quad_v.links(0)) == 4