FEniCS
diff --git a/‎cpp/dolfinx/fem/utils.h‎
Lines changed: 71 additions & 0 deletions b/‎cpp/dolfinx/fem/utils.h‎
Lines changed: 71 additions & 0 deletions
diff --git a/‎python/dolfinx/fem/__init__.py‎
Lines changed: 13 additions & 171 deletions b/‎python/dolfinx/fem/__init__.py‎
Lines changed: 13 additions & 171 deletions
@@ -31,6 +31,7 @@
 #include <format>
 #include <functional>
 #include <memory>
+#include <numeric>
 #include <optional>
 #include <ranges>
 #include <set>
@@ -1060,4 +1061,74 @@ Expression<T, U> create_expression(
   return create_expression(e, coeff_map, const_map, argument_space);
 }
 
+/// @brief Take an existing mesh and create a new mesh with its geometry
+/// interpolated into a new coordinate element.
+///
+/// The topology is shared between `mesh` and the returned mesh.
+///
+/// Useful for creating a higher-order mesh from a lower-order one for
+/// computation, or vice-versa, for IO.
+///
+/// @param[in] mesh Input mesh.
+/// @param[in] new_cmap Coordinate element for the new geometry.
+/// @param[in] reorder_fn Optional graph reordering function applied to
+/// the new geometry dofmap.
+/// @return A new mesh sharing the topology of `mesh` and with a
+/// geometry described by `new_cmap`.
+template <std::floating_point T>
+mesh::Mesh<T> interpolate_geometry(
+    std::shared_ptr<mesh::Mesh<T>> mesh, const CoordinateElement<T>& new_cmap,
+    const std::function<std::vector<int>(
+        const graph::AdjacencyList<std::int32_t>&)>& reorder_fn = nullptr)
+{
+  assert(mesh);
+  const CoordinateElement<T>& old_cmap = mesh->geometry().cmap();
+  if (new_cmap.cell_shape() != old_cmap.cell_shape())
+  {
+    throw std::runtime_error(
+        "Cell shape of new coordinate element must match input mesh.");
+  }
+
+  const int gdim = mesh->geometry().dim();
+
+  // Build a vector-valued Lagrange FiniteElement from the coordinate element.
+  basix::FiniteElement<T> b_element = basix::create_element<T>(
+      basix::element::family::P,
+      mesh::cell_type_to_basix_type(new_cmap.cell_shape()), new_cmap.degree(),
+      new_cmap.variant(), basix::element::dpc_variant::unset, false);
+  auto element = std::make_shared<const FiniteElement<T>>(
+      b_element, std::vector<std::size_t>{static_cast<std::size_t>(gdim)});
+
+  FunctionSpace<T> V = create_functionspace(mesh, element, reorder_fn);
+
+  // Tabulate physical coordinates of the new geometry dofs.
+  std::vector<T> x_new = V.tabulate_dof_coordinates(false);
+
+  // Pull the geometry dofmap and index map from V.
+  std::shared_ptr<const DofMap> dm = V.dofmap();
+  assert(dm);
+  std::shared_ptr<const common::IndexMap> new_imap = dm->index_map;
+  assert(new_imap);
+
+  auto map_view = dm->map();
+  std::vector<std::int32_t> dofmap_flat(
+      map_view.data_handle(), map_view.data_handle() + map_view.size());
+
+  // Build input_global_indices as the local-to-global of the new geometry
+  // dofs.
+  const std::int32_t num_nodes
+      = new_imap->size_local() + new_imap->num_ghosts();
+  std::vector<std::int32_t> local(num_nodes);
+  std::iota(local.begin(), local.end(), 0);
+  std::vector<std::int64_t> igi(num_nodes);
+  new_imap->local_to_global(local, igi);
+
+  mesh::Geometry<T> geometry(
+      new_imap, std::vector<std::vector<std::int32_t>>{std::move(dofmap_flat)},
+      std::vector<CoordinateElement<T>>{new_cmap}, std::move(x_new), gdim,
+      std::move(igi));
+
+  return mesh::Mesh<T>(mesh->comm(), mesh->topology(), std::move(geometry));
+}
+
 } // namespace dolfinx::fem
@@ -10,21 +10,8 @@
     dependencies and must be explicitly imported.
 """
 
-import typing
-
-import numpy as np
-import numpy.typing as npt
-
 from dolfinx.cpp.fem import _IntegralType as IntegralType
-from dolfinx.cpp.fem import build_sparsity_pattern as _build_sparsity_pattern
-from dolfinx.cpp.fem import compute_integration_domains as _compute_integration_domains
-from dolfinx.cpp.fem import create_interpolation_data as _create_interpolation_data
-from dolfinx.cpp.fem import create_sparsity_pattern as _create_sparsity_pattern
-from dolfinx.cpp.fem import discrete_curl as _discrete_curl
-from dolfinx.cpp.fem import discrete_gradient as _discrete_gradient
-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.fem.assemble import (
     apply_lifting,
     assemble_matrix,
@@ -62,164 +49,16 @@
     FunctionSpace,
     functionspace,
 )
-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.
-
-    Args:
-        a: Bilinear form to build a sparsity pattern for.
-
-    Returns:
-        Sparsity pattern for the form ``a``.
-
-    Note:
-        The pattern is not finalised, i.e. the caller is responsible for
-        calling ``assemble`` on the sparsity pattern.
-    """
-    return _create_sparsity_pattern(a._cpp_object)
-
-
-def build_sparsity_pattern(pattern: SparsityPattern, a: Form):
-    """Build a sparsity pattern from a bilinear form.
-
-    Args:
-        pattern: The sparsity pattern to add to
-        a: Bilinear form to build a sparsity pattern for.
-
-    Returns:
-        Sparsity pattern for the form ``a``.
-
-    Note:
-        The pattern is not finalised, i.e. the caller is responsible for
-        calling ``assemble`` on the sparsity pattern.
-    """
-    return _build_sparsity_pattern(pattern, a._cpp_object)
-
-
-def create_interpolation_data(
-    V_to: FunctionSpace,
-    V_from: FunctionSpace,
-    cells: npt.NDArray[np.int32],
-    padding: float = 1e-14,
-) -> _PointOwnershipData:
-    """Generate data for interpolating functions on different meshes.
-
-    Args:
-        V_to: Function space to interpolate into.
-        V_from: Function space to interpolate from.
-        cells: Indices of the cells associated with `V_to` on which to
-            interpolate into.
-        padding: Absolute padding of bounding boxes of all entities on
-            mesh_to.
-
-    Returns:
-        Data needed to interpolation functions defined on function
-        spaces on the meshes.
-    """
-    return _PointOwnershipData(
-        _create_interpolation_data(
-            V_to.mesh._cpp_object.geometry,
-            V_to.element._cpp_object,
-            V_from.mesh._cpp_object,
-            cells,
-            padding,
-        )
-    )
-
-
-def discrete_curl(V0: FunctionSpace, V1: FunctionSpace) -> _MatrixCSR:
-    """Assemble a discrete curl operator.
-
-    The discrete curl operator interpolates the curl of H(curl) finite
-    element function into a H(div) space.
-
-    Args:
-        V0: H1(curl) space to interpolate the curl from.
-        V1: H(div) space to interpolate into.
-
-    Returns:
-        Discrete curl operator.
-    """
-    return _MatrixCSR(_discrete_curl(V0._cpp_object, V1._cpp_object))
-
-
-def discrete_gradient(space0: FunctionSpace, space1: FunctionSpace) -> _MatrixCSR:
-    """Assemble a discrete gradient operator.
-
-    The discrete gradient operator interpolates the gradient of a H1
-    finite element function into a H(curl) space. It is assumed that the
-    H1 space uses an identity map and the H(curl) space uses a covariant
-    Piola map.
-
-    Args:
-        space0: H1 space to interpolate the gradient from.
-        space1: H(curl) space to interpolate into.
-
-    Returns:
-        Discrete gradient operator.
-    """
-    return _MatrixCSR(_discrete_gradient(space0._cpp_object, space1._cpp_object))
-
-
-def interpolation_matrix(space0: FunctionSpace, space1: FunctionSpace) -> _MatrixCSR:
-    """Create interpolation matrix between spaces on the same mesh.
-
-    Args:
-        space0: space to interpolate from.
-        space1: space to interpolate into.
-
-    Returns:
-        Interpolation matrix.
-
-    Note:
-        The returned matrix is not finalised, i.e. ghost values are not
-        accumulated.
-    """
-    return _MatrixCSR(_interpolation_matrix(space0._cpp_object, space1._cpp_object))
-
-
-def compute_integration_domains(
-    integral_type: IntegralType, topology: "dolfinx.mesh.Topology", entities: np.ndarray
-):
-    """Determine compute integration entities.
-
-    This function returns a list ``[(id, entities)]``. For cell
-    integrals ``entities`` are the cell indices. For exterior facet
-    integrals, ``entities`` is a list of ``(cell_index,
-    local_facet_index)`` pairs. For interior facet integrals,
-    ``entities`` is a list of ``(cell_index0, local_facet_index0,
-    cell_index1, local_facet_index1)``. ``id`` refers to the subdomain
-    id used in the definition of the integration measures of the
-    variational form.
-
-    Note:
-        Owned mesh entities only are returned. Ghost entities are not
-        included.
-
-    Note:
-        For facet integrals, the topology facet-to-cell and
-        cell-to-facet connectivity must be computed before calling this
-        function.
-
-    Args:
-        integral_type: Integral type.
-        topology: Mesh topology.
-        entities: List of mesh entities. For
-            ``integral_type==IntegralType.cell``, ``entities`` should be
-            cell indices. For other ``IntegralType``s, ``entities``
-            should be facet indices.
-
-    Returns:
-        List of integration entities.
-    """
-    return _compute_integration_domains(integral_type, topology._cpp_object, entities)
-
+from dolfinx.fem.utils import (
+    build_sparsity_pattern,
+    compute_integration_domains,
+    create_interpolation_data,
+    create_sparsity_pattern,
+    discrete_curl,
+    discrete_gradient,
+    interpolate_geometry,
+    interpolation_matrix,
+)
 
 __all__ = [
     "Constant",
@@ -238,6 +77,7 @@ def compute_integration_domains(
     "assemble_scalar",
     "assemble_vector",
     "bcs_by_block",
+    "build_sparsity_pattern",
     "compile_form",
     "compute_integration_domains",
     "coordinate_element",
@@ -255,6 +95,8 @@ def compute_integration_domains(
     "form",
     "form_cpp_class",
     "functionspace",
+    "interpolate_geometry",
+    "interpolation_matrix",
     "locate_dofs_geometrical",
     "locate_dofs_topological",
     "mixed_topology_form",