Introduce dolfinx.fem.create_dofmaps (#3842)

* Wrap create_dofmaps in Python layer

* Type oblivious

* Use python type for reference

* Remove

* Remove empty

* Without checking?

* ruff

* Try again

* Fully qualify

* Adapt demo

* Fix

* Ruff

* Simplify

* One more

* Test case

* Update python/dolfinx/fem/dofmap.py

Co-authored-by: Garth N. Wells <gnw20@cam.ac.uk>

* Remove coordinate elements naming

* Improve docstring

* Update python/dolfinx/fem/dofmap.py

Co-authored-by: Garth N. Wells <gnw20@cam.ac.uk>

* rst

---------

Co-authored-by: Garth N. Wells <gnw20@cam.ac.uk>

Author: schnellerhase

python/demo/demo_mixed-topology.py

@@ -26,14 +26,16 @@
 import ufl
 from dolfinx.cpp.mesh import GhostMode, create_cell_partitioner, create_mesh
 from dolfinx.fem import (
+    FiniteElement,
     FunctionSpace,
     assemble_matrix,
     assemble_vector,
     coordinate_element,
+    create_dofmaps,
     mixed_topology_form,
 )
 from dolfinx.io.utils import cell_perm_vtk
-from dolfinx.mesh import CellType, Mesh
+from dolfinx.mesh import CellType, Mesh, Topology
 
 if MPI.COMM_WORLD.size > 1:
     print("Not yet running in parallel")
@@ -102,12 +104,20 @@
     basix.create_element(basix.ElementFamily.P, basix.CellType.hexahedron, 1),
     basix.create_element(basix.ElementFamily.P, basix.CellType.prism, 1),
 ]
-elements_cpp = [_cpp.fem.FiniteElement_float64(e._e, None, True) for e in elements]
+dolfinx_elements = [
+    FiniteElement(_cpp.fem.FiniteElement_float64(e._e, None, True)) for e in elements
+]
 # NOTE: Both dofmaps have the same IndexMap, but different cell_dofs
-dofmaps = _cpp.fem.create_dofmaps(mesh.comm, mesh.topology, elements_cpp)
+dofmaps = create_dofmaps(
+    mesh.comm,
+    Topology(mesh.topology),
+    dolfinx_elements,
+)
 
 # Create C++ function space
-V_cpp = _cpp.fem.FunctionSpace_float64(mesh, elements_cpp, dofmaps)
+V_cpp = _cpp.fem.FunctionSpace_float64(
+    mesh, [e._cpp_object for e in dolfinx_elements], [dofmap._cpp_object for dofmap in dofmaps]
+)
 
 # Create forms for each cell type.
 # FIXME This hack is required at the moment because UFL does not yet know

python/dolfinx/fem/__init__.py

@@ -5,6 +5,8 @@
 # SPDX-License-Identifier:    LGPL-3.0-or-later
 """Tools for assembling and manipulating finite element forms."""
 
+import typing
+
 import numpy as np
 import numpy.typing as npt
 
@@ -18,7 +20,6 @@
 from dolfinx.cpp.fem import interpolation_matrix as _interpolation_matrix
 from dolfinx.cpp.fem import transpose_dofmap
 from dolfinx.cpp.la import SparsityPattern
-from dolfinx.cpp.mesh import Topology
 from dolfinx.fem.assemble import (
     apply_lifting,
     assemble_matrix,
@@ -37,7 +38,7 @@
     locate_dofs_geometrical,
     locate_dofs_topological,
 )
-from dolfinx.fem.dofmap import DofMap
+from dolfinx.fem.dofmap import DofMap, create_dofmaps
 from dolfinx.fem.element import CoordinateElement, FiniteElement, coordinate_element, finiteelement
 from dolfinx.fem.forms import (
     Form,
@@ -59,6 +60,9 @@
 from dolfinx.geometry import PointOwnershipData as _PointOwnershipData
 from dolfinx.la import MatrixCSR as _MatrixCSR
 
+if typing.TYPE_CHECKING:
+    import dolfinx.mesh
+
 
 def create_sparsity_pattern(a: Form):
     """Create a sparsity pattern from a bilinear form.
@@ -174,7 +178,7 @@ def interpolation_matrix(space0: FunctionSpace, space1: FunctionSpace) -> _Matri
 
 
 def compute_integration_domains(
-    integral_type: IntegralType, topology: Topology, entities: np.ndarray
+    integral_type: IntegralType, topology: "dolfinx.mesh.Topology", entities: np.ndarray
 ):
     """Given an integral type and a set of entities compute integration
     entities.
@@ -231,6 +235,7 @@ def compute_integration_domains(
     "compile_form",
     "compute_integration_domains",
     "coordinate_element",
+    "create_dofmaps",
     "create_form",
     "create_interpolation_data",
     "create_matrix",

python/dolfinx/fem/dofmap.py

@@ -1,10 +1,20 @@
-# Copyright (C) 2018 Michal Habera
+# Copyright (C) 2018-2025 Michal Habera and Paul T. Kühner
 #
 # This file is part of DOLFINx (https://www.fenicsproject.org)
 #
 # SPDX-License-Identifier:    LGPL-3.0-or-later
 
-from dolfinx import cpp as _cpp
+import typing
+from collections.abc import Sequence
+
+from mpi4py.MPI import Comm
+
+from basix.finite_element import FiniteElement
+from dolfinx.cpp.fem import DofMap as _DofMap
+from dolfinx.cpp.fem import create_dofmaps as _create_dofmaps
+
+if typing.TYPE_CHECKING:
+    import dolfinx.mesh
 
 
 class DofMap:
@@ -14,9 +24,9 @@ class DofMap:
     dof map based on a FiniteElement on a specific mesh.
     """
 
-    _cpp_object: _cpp.fem.DofMap
+    _cpp_object: _DofMap
 
-    def __init__(self, dofmap: _cpp.fem.DofMap):
+    def __init__(self, dofmap: _DofMap):
         self._cpp_object = dofmap
 
     def cell_dofs(self, cell_index: int):
@@ -56,3 +66,22 @@ def index_map_bs(self):
     def list(self):
         """Adjacency list with dof indices for each cell."""
         return self._cpp_object.map()
+
+
+def create_dofmaps(
+    comm: Comm, topology: "dolfinx.mesh.Topology", elements: Sequence[FiniteElement]
+) -> list[DofMap]:
+    """Create degree-of-freedom maps on a given topology.
+
+    Args:
+        comm: MPI communicator
+        topology: Mesh topology
+        elements: Sequence of elements
+
+    Returns:
+        List of degree-of-freedom maps where the ``i``-th map is the map
+        for ``elements[i]``.
+    """
+    elements_cpp = [e._cpp_object for e in elements]
+    cpp_dofmaps = _create_dofmaps(comm, topology._cpp_object, elements_cpp)
+    return [DofMap(cpp_object) for cpp_object in cpp_dofmaps]

python/test/unit/fem/test_mixed_mesh_dofmap.py

@@ -3,23 +3,19 @@
 import numpy as np
 
 import basix
-import dolfinx.cpp as _cpp
 from dolfinx.cpp.mesh import Mesh_float64, create_geometry, create_topology
-from dolfinx.fem import coordinate_element
-from dolfinx.fem.dofmap import DofMap
+from dolfinx.fem import coordinate_element, create_dofmaps
+from dolfinx.fem.element import finiteelement
 from dolfinx.log import LogLevel, set_log_level
-from dolfinx.mesh import CellType
+from dolfinx.mesh import CellType, Topology
 
 
 def create_element_dofmap(mesh, cell_types, degree):
-    cpp_elements = []
-    for cell_type in cell_types:
-        ufl_e = basix.ufl.element("P", cell_type, degree, dtype=np.float64)
-        cpp_elements += [_cpp.fem.FiniteElement_float64(ufl_e.basix_element._e, None, False)]
-
-    cpp_dofmaps = _cpp.fem.create_dofmaps(mesh.comm, mesh.topology, cpp_elements)
-
-    return (cpp_elements, cpp_dofmaps)
+    elements = [
+        finiteelement(ct, basix.ufl.element("P", ct, degree), np.float64) for ct in cell_types
+    ]
+    dofmaps = create_dofmaps(mesh.comm, Topology(mesh.topology), elements)
+    return (elements, dofmaps)
 
 
 def test_dofmap_mixed_topology():
@@ -72,8 +68,7 @@ def test_dofmap_mixed_topology():
     assert elements[1].basix_element.cell_type.name == "quadrilateral"
 
     assert len(dofmaps) == 2
-    q0 = DofMap(dofmaps[0])
-    q1 = DofMap(dofmaps[1])
+    q0, q1 = dofmaps
     assert q0.index_map.size_local == q1.index_map.size_local
     # Triangles
     print(q0.list)
@@ -121,10 +116,9 @@ def test_dofmap_prism_mesh():
     mesh = Mesh_float64(MPI.COMM_WORLD, topology, geom)
 
     elements, dofmaps = create_element_dofmap(mesh, [basix.CellType.prism], 2)
-    print()
     assert len(elements) == 1
     assert len(dofmaps) == 1
-    q = DofMap(dofmaps[0])
+    q = dofmaps[0]
     assert q.index_map.size_local == 18
     print(q.list)
     facet_dofs = []