p-refinement of mesh geometry

cpp/dolfinx/mesh/utils.cpp

@@ -11,13 +11,22 @@
 #include "cell_types.h"
 #include "graphbuild.h"
 #include <algorithm>
+#include <basix/element-families.h>
+#include <basix/finite-element.h>
 #include <cstdlib>
 #include <dolfinx/common/IndexMap.h>
 #include <dolfinx/common/log.h>
 #include <dolfinx/common/math.h>
+#include <dolfinx/fem/CoordinateElement.h>
+#include <dolfinx/fem/DofMap.h>
 #include <dolfinx/fem/ElementDofLayout.h>
+#include <dolfinx/fem/FiniteElement.h>
+#include <dolfinx/fem/FunctionSpace.h>
+#include <dolfinx/fem/utils.h>
 #include <dolfinx/graph/AdjacencyList.h>
 #include <dolfinx/graph/partition.h>
+#include <memory>
+#include <numeric>
 #include <optional>
 #include <span>
 #include <stdexcept>
@@ -168,3 +177,70 @@
   return entities1;
 }
 //-----------------------------------------------------------------------------
+template <std::floating_point T>
+mesh::Mesh<T> mesh::interpolate_geometry(
+    std::shared_ptr<mesh::Mesh<T>> mesh,
+    const fem::CoordinateElement<T>& new_cmap,
+    const std::function<std::vector<int>(
+        const graph::AdjacencyList<std::int32_t>&)>& reorder_fn)
+{
+  assert(mesh);
+  const fem::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 fem::FiniteElement<T>>(
+      b_element, std::vector<std::size_t>{gdim});
+
+  fem::FunctionSpace<T> V
+      = fem::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 fem::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);
+
+  Geometry<T> geometry(
+      new_imap, std::vector<std::vector<std::int32_t>>{std::move(dofmap_flat)},
+      std::vector<fem::CoordinateElement<T>>{new_cmap}, std::move(x_new), gdim,
+      std::move(igi));
+
+  return Mesh<T>(mesh->comm(), mesh->topology(), std::move(geometry));
+}
+//-----------------------------------------------------------------------------
+template mesh::Mesh<float> mesh::interpolate_geometry(
+    std::shared_ptr<mesh::Mesh<float>>, const fem::CoordinateElement<float>&,
+    const std::function<
+        std::vector<int>(const graph::AdjacencyList<std::int32_t>&)>&);
+template mesh::Mesh<double> mesh::interpolate_geometry(
+    std::shared_ptr<mesh::Mesh<double>>, const fem::CoordinateElement<double>&,
+    const std::function<
+        std::vector<int>(const graph::AdjacencyList<std::int32_t>&)>&);
+//-----------------------------------------------------------------------------

cpp/dolfinx/mesh/utils.h

@@ -1410,6 +1410,25 @@ create_subgeometry(const Mesh<T>& mesh, int dim,
           std::move(subx_to_x_dofmap)};
 }
 
+/// @brief Make a copy of `mesh` whose geometry is interpolated into a
+/// new coordinate element.
+///
+/// The topology is shared (pointer-level) between the new and old mesh.
+///
+/// @param[in] mesh Input mesh.
+/// @param[in] new_cmap Target 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<T>
+interpolate_geometry(std::shared_ptr<Mesh<T>> mesh,
+                     const fem::CoordinateElement<T>& new_cmap,
+                     const std::function<std::vector<int>(
+                         const graph::AdjacencyList<std::int32_t>&)>& reorder_fn
+                     = nullptr);
+
 /// @brief Create a new mesh consisting of a subset of entities in a
 /// mesh.
 /// @param[in] mesh The mesh.

python/dolfinx/mesh.py

@@ -73,6 +73,7 @@
     "create_unit_square",
     "entities_to_geometry",
     "exterior_facet_indices",
+    "interpolate_geometry",
     "locate_entities",
     "locate_entities_boundary",
     "mark_maximum",
@@ -873,6 +874,39 @@ def create_submesh(
     return (Mesh(submsh, submsh_domain), EntityMap(entity_map), EntityMap(vertex_map), geom_map)
 
 
+def interpolate_geometry(msh: Mesh, cmap: _CoordinateElement) -> Mesh:
+    """Create a copy of a mesh with geometry interpolated into cmap.
+
+    Useful for lifting first-order meshes from built-in mesh generation
+    functions to higher-order geometry.
+
+    Note:
+        The new topology is shared (no copy) with the input mesh. The new
+        geometry is created by interpolating the existing geometry at the
+        reference interpolation points of ``cmap``.
+
+    Args:
+        msh: Input mesh.
+        cmap: Target coordinate element for the new geometry.
+
+    Returns:
+        A new mesh sharing the topology of ``msh`` and with geometry in
+        ``cmap``.
+    """
+    new_msh = _cpp.mesh.interpolate_geometry(msh._cpp_object, cmap._cpp_object)
+    domain = ufl.Mesh(
+        basix.ufl.element(
+            "Lagrange",
+            to_string(new_msh.topology.cell_type),
+            new_msh.geometry.cmap().degree,
+            basix.LagrangeVariant(new_msh.geometry.cmap().variant),
+            shape=(new_msh.geometry.dim,),
+            dtype=new_msh.geometry.x.dtype,
+        )
+    )
+    return Mesh(new_msh, domain)
+
+
 def meshtags(
     msh: Mesh,
     dim: int,

python/test/unit/mesh/test_interpolate_geometry.py

@@ -0,0 +1,71 @@
+# Copyright (C) 2026 Jack S. Hale
+#
+# This file is part of DOLFINx (https://www.fenicsproject.org)
+#
+# SPDX-License-Identifier:    LGPL-3.0-or-later
+
+from mpi4py import MPI
+
+import numpy as np
+import pytest
+
+from dolfinx.fem import coordinate_element
+from dolfinx.mesh import CellType, create_unit_square, interpolate_geometry
+
+
+def _assert_close_up_to_row_permutation(a, b, atol):
+    """Assert rows of a and b are equal, up to a row permutation"""
+    a = np.asarray(a)
+    b = np.asarray(b)
+    assert a.shape == b.shape
+    used = [False] * len(b)
+    for ra in a:
+        for j, rb in enumerate(b):
+            if not used[j] and np.allclose(ra, rb, atol=atol, rtol=0.0):
+                used[j] = True
+                break
+        else:
+            raise AssertionError(f"No match for {ra} in {b}")
+
+
+@pytest.mark.parametrize("dtype", [np.float32, np.float64])
+def test_interpolate_geometry_p1_to_p2(dtype):
+    msh = create_unit_square(MPI.COMM_WORLD, 4, 4, dtype=dtype)
+    cmap = coordinate_element(CellType.triangle, 2, dtype=dtype)
+
+    new_msh = interpolate_geometry(msh, cmap)
+    assert new_msh.geometry.cmap().degree == 2
+
+    # Topology should be the same cpp object.
+    assert new_msh.topology._cpp_object is msh.topology._cpp_object
+
+    x_old, x_new = msh.geometry.x, new_msh.geometry.x
+    dm_old, dm_new = msh.geometry.dofmap, new_msh.geometry.dofmap
+    assert dm_old.shape[0] == dm_new.shape[0]
+    assert dm_new.shape[1] == 6
+
+    atol = 10 * np.finfo(dtype).eps
+    for c in range(dm_new.shape[0]):
+        v = x_old[dm_old[c]]
+        # P2 geometry has original vertices + midpoints.
+        expected = np.vstack([v, 0.5 * (v[0] + v[1]), 0.5 * (v[0] + v[2]), 0.5 * (v[1] + v[2])])
+        _assert_close_up_to_row_permutation(x_new[dm_new[c]], expected, atol=atol)
+
+
+@pytest.mark.parametrize("dtype", [np.float32, np.float64])
+def test_interpolate_geometry_p1_roundtrip(dtype):
+    msh = create_unit_square(MPI.COMM_WORLD, 4, 4, dtype=dtype)
+    cmap = coordinate_element(CellType.triangle, 1, dtype=dtype)
+
+    new_msh = interpolate_geometry(msh, cmap)
+
+    assert new_msh.geometry.cmap().degree == 1
+    assert new_msh.topology._cpp_object is msh.topology._cpp_object
+
+    x_old, x_new = msh.geometry.x, new_msh.geometry.x
+    dm_old, dm_new = msh.geometry.dofmap, new_msh.geometry.dofmap
+    assert dm_old.shape == dm_new.shape
+
+    atol = 10 * np.finfo(dtype).eps
+    for c in range(dm_old.shape[0]):
+        np.testing.assert_allclose(x_new[dm_new[c]], x_old[dm_old[c]], atol=atol, rtol=0.0)