Add Python-interface for geometry constructor and geometry class (FEn…

…iCS#3403)

* Expose geometry constructor

* Accept 2D arrays for dofmap and points + simplify wrapper

* Revert weird linebreak

* Send in ((num_points, gdim) nodes and pad in nanobind

* Add Python-interface for geometry constructor and geometry class with documentation

* Fix geometry access

* Fix docstring

* Add python-interface of all of CoordinateElement.

* Doc fix

* Ruff formatting

* Fix class doc

* Wrap geometry once

* Fix geometry wrapper

* Apply suggestions from code review

* Apply suggestions from code review

* Apply suggestions from code review

* Modifying docstrings

* Apply suggestions from code review

* Add full stop

* Another full stop

* Stops all around

* Final set of full stops?

* Add type hints for push_forward/pull_back

* Use old typing to work on redhat

* Apply suggestions from @garth-wells

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

* Adhere to suggestions

* Add hsape

* Ruff formatting

* Fix indentation

* Small doc updates

* Doc formatting fix

---------

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

cpp/dolfinx/mesh/Geometry.h

@@ -113,7 +113,7 @@ class Geometry
   /// Move Assignment
   Geometry& operator=(Geometry&&) = default;
 
-  /// Return Euclidean dimension of coordinate system
+  /// Return dimension of the Euclidean coordinate system
   int dim() const { return _dim; }
 
   /// @brief  DofMap for the geometry

python/dolfinx/fem/element.py

@@ -16,7 +16,7 @@
 
 
 class CoordinateElement:
-    """Coordinate element describing the geometry map for mesh cells"""
+    """Coordinate element describing the geometry map for mesh cells."""
 
     _cpp_object: typing.Union[
         _cpp.fem.CoordinateElement_float32, _cpp.fem.CoordinateElement_float64
@@ -46,6 +46,75 @@ def dtype(self) -> np.dtype:
         """Scalar type for the coordinate element."""
         return np.dtype(self._cpp_object.dtype)
 
+    @property
+    def dim(self) -> int:
+        """Dimension of the coordinate element space.
+
+        This is number of basis functions that span the coordinate
+        space, e.g., for a linear triangle cell the dimension will be 3.
+        """
+        return self._cpp_object.dim
+
+    def create_dof_layout(self) -> _cpp.fem.ElementDofLayout:
+        """Compute and return the dof layout"""
+        return self._cpp_object.create_dof_layout()
+
+    def push_forward(
+        self,
+        X: typing.Union[npt.NDArray[np.float32], npt.NDArray[np.float64]],
+        cell_geometry: typing.Union[npt.NDArray[np.float32], npt.NDArray[np.float64]],
+    ) -> typing.Union[npt.NDArray[np.float32], npt.NDArray[np.float64]]:
+        """Push points on the reference cell forward to the physical cell.
+
+        Args:
+            X: Coordinates of points on the reference cell,
+                ``shape=(num_points, topological_dimension)``.
+            cell_geometry: Coordinate 'degrees-of-freedom' (nodes) of
+                the cell, ``shape=(num_geometry_basis_functions,
+                geometrical_dimension)``. Can be created by accessing
+                ``geometry.x[geometry.dofmap.cell_dofs(i)]``,
+
+        Returns:
+            Physical coordinates of the points reference points ``X``.
+        """
+        return self._cpp_object.push_forward(X, cell_geometry)
+
+    def pull_back(
+        self,
+        x: typing.Union[npt.NDArray[np.float32], npt.NDArray[np.float64]],
+        cell_geometry: typing.Union[npt.NDArray[np.float32], npt.NDArray[np.float64]],
+    ) -> typing.Union[npt.NDArray[np.float32], npt.NDArray[np.float64]]:
+        """Pull points on the physical cell back to the reference cell.
+
+        For non-affine cells, the pull-back is a nonlinear operation.
+
+        Args:
+            x: Physical coordinates to pull back to the reference cells,
+                ``shape=(num_points, geometrical_dimension)``.
+            cell_geometry: Physical coordinates describing the cell,
+                shape ``(num_of_geometry_basis_functions, geometrical_dimension)``
+                They can be created by accessing `geometry.x[geometry.dofmap.cell_dofs(i)]`,
+
+        Returns:
+            Reference coordinates of the physical points ``x``.
+        """
+        return self._cpp_object.pull_back(x, cell_geometry)
+
+    @property
+    def variant(self) -> int:
+        """Lagrange variant of the coordinate element.
+
+        Note:
+            The return type is an integer. A Basix enum can be created using
+            ``basix.LagrangeVariant(value)``.
+        """
+        return self._cpp_object.variant
+
+    @property
+    def degree(self) -> int:
+        """Polynomial degree of the coordinate element."""
+        return self._cpp_object.degree
+
 
 @singledispatch
 def coordinate_element(

python/dolfinx/fem/function.py

@@ -461,7 +461,9 @@ def _(e0: Expression):
         except TypeError:
             # u0 is callable
             assert callable(u0)
-            x = _cpp.fem.interpolation_coords(self._V.element, self._V.mesh.geometry, cells0)
+            x = _cpp.fem.interpolation_coords(
+                self._V.element, self._V.mesh.geometry._cpp_object, cells0
+            )
             self._cpp_object.interpolate(np.asarray(u0(x), dtype=self.dtype), cells0)  # type: ignore
 
     def copy(self) -> Function:

python/dolfinx/io/utils.py

@@ -21,7 +21,7 @@
 from dolfinx.cpp.io import perm_gmsh as cell_perm_gmsh
 from dolfinx.cpp.io import perm_vtk as cell_perm_vtk
 from dolfinx.fem import Function
-from dolfinx.mesh import GhostMode, Mesh, MeshTags
+from dolfinx.mesh import Geometry, GhostMode, Mesh, MeshTags
 
 __all__ = ["VTKFile", "XDMFFile", "cell_perm_gmsh", "cell_perm_vtk", "distribute_entity_data"]
 
@@ -225,12 +225,12 @@ def write_mesh(self, mesh: Mesh, xpath: str = "/Xdmf/Domain") -> None:
     def write_meshtags(
         self,
         tags: MeshTags,
-        x: typing.Union[_cpp.mesh.Geometry_float32, _cpp.mesh.Geometry_float64],
+        x: Geometry,
         geometry_xpath: str = "/Xdmf/Domain/Grid/Geometry",
         xpath: str = "/Xdmf/Domain",
     ) -> None:
         """Write mesh tags to file"""
-        super().write_meshtags(tags._cpp_object, x, geometry_xpath, xpath)
+        super().write_meshtags(tags._cpp_object, x._cpp_object, geometry_xpath, xpath)
 
     def write_function(
         self, u: Function, t: float = 0.0, mesh_xpath="/Xdmf/Domain/Grid[@GridType='Uniform'][1]"

python/dolfinx/mesh.py

@@ -1,4 +1,4 @@
-# Copyright (C) 2017-2021 Chris N. Richardson and Garth N. Wells
+# Copyright (C) 2017-2024 Chris N. Richardson, Garth N. Wells and Jørgen S. Dokken
 #
 # This file is part of DOLFINx (https://www.fenicsproject.org)
 #
@@ -19,6 +19,7 @@
 import ufl
 from dolfinx import cpp as _cpp
 from dolfinx import default_real_type
+from dolfinx.common import IndexMap as _IndexMap
 from dolfinx.cpp.mesh import (
     CellType,
     DiagonalType,
@@ -62,13 +63,67 @@
     "to_string",
     "transfer_meshtag",
     "entities_to_geometry",
+    "create_geometry",
+    "Geometry",
 ]
 
 
+class Geometry:
+    """The geometry of a :class:`dolfinx.mesh.Mesh`"""
+
+    _cpp_object: typing.Union[_cpp.mesh.Geometry_float32, _cpp.mesh.Geometry_float64]
+
+    def __init__(
+        self, geometry: typing.Union[_cpp.mesh.Geometry_float32, _cpp.mesh.Geometry_float64]
+    ):
+        """Initialize a geometry from a C++ geometry.
+
+        Args:
+            geometry: The C++ geometry object.
+
+        Note:
+            Geometry objects should usually be constructed with the :func:`create_geometry`
+            and not using this class initializer. This class is combined with different base classes
+            that depend on the scalar type used in the Geometry.
+        """
+
+        self._cpp_object = geometry
+
+    @property
+    def cmap(self) -> _CoordinateElement:
+        """Element that describes the geometry map."""
+        return _CoordinateElement(self._cpp_object.cmap)
+
+    @property
+    def dim(self):
+        """Dimension of the Euclidean coordinate system."""
+        return self._cpp_object.dim
+
+    @property
+    def dofmap(self) -> npt.NDArray[np.int32]:
+        """Dofmap for the geometry, shape ``(num_cells, dofs_per_cell)``."""
+        return self._cpp_object.dofmap
+
+    def index_map(self) -> _IndexMap:
+        """Index map describing the layout of the geometry points (nodes)."""
+        return self._cpp_object.index_map()
+
+    @property
+    def input_global_indices(self) -> npt.NDArray[np.int64]:
+        """Global input indices of the geometry nodes."""
+        return self._cpp_object.input_global_indices
+
+    @property
+    def x(self) -> typing.Union[npt.NDArray[np.float32], npt.NDArray[np.float64]]:
+        """Geometry coordinate points,  ``shape=(num_points, 3)``."""
+        return self._cpp_object.x
+
+
 class Mesh:
     """A mesh."""
 
     _mesh: typing.Union[_cpp.mesh.Mesh_float32, _cpp.mesh.Mesh_float64]
+    _geometry: Geometry
     _ufl_domain: typing.Optional[ufl.Mesh]
 
     def __init__(self, mesh, domain: typing.Optional[ufl.Mesh]):
@@ -79,10 +134,12 @@ def __init__(self, mesh, domain: typing.Optional[ufl.Mesh]):
             domain: A UFL domain.
 
         Note:
-            Mesh objects should not usually be created using this
-            initializer directly.
+            Mesh objects should usually be constructed using :func:`create_mesh`
+            and not using this class initializer. This class is combined with different base classes
+            that depend on the scalar type used in the Mesh.
         """
         self._cpp_object = mesh
+        self._geometry = Geometry(self._cpp_object.geometry)
         self._ufl_domain = domain
         if self._ufl_domain is not None:
             self._ufl_domain._ufl_cargo = self._cpp_object  # type: ignore
@@ -141,9 +198,9 @@ def topology(self):
         return self._cpp_object.topology
 
     @property
-    def geometry(self):
+    def geometry(self) -> Geometry:
         "Mesh geometry."
-        return self._cpp_object.geometry
+        return self._geometry
 
 
 class MeshTags:
@@ -762,3 +819,34 @@ def entities_to_geometry(
         The geometric DOFs associated with the closure of the entities in `entities`.
     """
     return _cpp.mesh.entities_to_geometry(mesh._cpp_object, dim, entities, permute)
+
+
+def create_geometry(
+    index_map: _IndexMap,
+    dofmap: npt.NDArray[np.int32],
+    element: _CoordinateElement,
+    x: np.ndarray,
+    input_global_indices: npt.NDArray[np.int64],
+) -> Geometry:
+    """Create a Geometry object.
+
+    Args:
+        index_map: Index map describing the layout of the geometry points (nodes).
+        dofmap: The geometry (point) dofmap. For a cell, it gives the row
+            in the point coordinates ``x`` of each local geometry node.
+            ``shape=(num_cells, num_dofs_per_cell)``.
+        element: Element that describes the cell geometry map.
+        x: The point coordinates. The shape is ``(num_points, geometric_dimension).``
+        input_global_indices: The 'global' input index of each point, commonly
+            from a mesh input file.
+    """
+    if x.dtype == np.float64:
+        ftype = _cpp.mesh.Geometry_float64
+    elif x.dtype == np.float32:
+        ftype = _cpp.mesh.Geometry_float64
+    else:
+        raise ValueError("Unknown floating type for geometry, got: {x.dtype}")
+
+    if (dtype := np.dtype(element.dtype)) != x.dtype:
+        raise ValueError(f"Mismatch in x dtype ({x.dtype}) and coordinate element ({dtype})")
+    return Geometry(ftype(index_map, dofmap, element, x, input_global_indices))

python/dolfinx/wrappers/mesh.cpp

@@ -239,6 +239,9 @@ void declare_mesh(nb::module_& m, std::string type)
       .def_prop_ro(
           "cmap", [](dolfinx::mesh::Geometry<T>& self) { return self.cmap(); },
           "The coordinate map")
+      .def(
+          "cmaps", [](dolfinx::mesh::Geometry<T>& self, int i)
+          { return self.cmap(i); }, "The ith coordinate map")
       .def_prop_ro(
           "input_global_indices",
           [](const dolfinx::mesh::Geometry<T>& self)