diff --git a/cpp/demo/interpolation_different_meshes/main.cpp b/cpp/demo/interpolation_different_meshes/main.cpp
index aaf927d1b5f..30133693218 100644
--- a/cpp/demo/interpolation_different_meshes/main.cpp
+++ b/cpp/demo/interpolation_different_meshes/main.cpp
@@ -66,7 +66,9 @@ int main(int argc, char* argv[])
     // Interpolate from u_tet to u_hex
     auto nmm_interpolation_data
         = fem::create_nonmatching_meshes_interpolation_data(
-            *u_hex->function_space(), *u_tet->function_space());
+            *u_hex->function_space()->mesh(),
+            *u_hex->function_space()->element(),
+            *u_tet->function_space()->mesh());
     u_hex->interpolate(*u_tet, nmm_interpolation_data);
 
 #ifdef HAS_ADIOS2
diff --git a/cpp/dolfinx/fem/Function.h b/cpp/dolfinx/fem/Function.h
index 8891dd96138..e6f0c96ab5b 100644
--- a/cpp/dolfinx/fem/Function.h
+++ b/cpp/dolfinx/fem/Function.h
@@ -189,8 +189,9 @@ class Function
     assert(_function_space);
     assert(_function_space->element());
     assert(_function_space->mesh());
-    const std::vector<double> x = fem::interpolation_coords(
-        *_function_space->element(), *_function_space->mesh(), cells);
+    const std::vector<double> x
+        = fem::interpolation_coords(*_function_space->element(),
+                                    _function_space->mesh()->geometry(), cells);
     namespace stdex = std::experimental;
     stdex::mdspan<const double,
                   stdex::extents<std::size_t, 3, stdex::dynamic_extent>>
diff --git a/cpp/dolfinx/fem/interpolate.cpp b/cpp/dolfinx/fem/interpolate.cpp
index a4fcd970154..7e8c069c8d0 100644
--- a/cpp/dolfinx/fem/interpolate.cpp
+++ b/cpp/dolfinx/fem/interpolate.cpp
@@ -14,15 +14,15 @@ using namespace dolfinx;
 
 //-----------------------------------------------------------------------------
 std::vector<double>
-fem::interpolation_coords(const FiniteElement& element, const mesh::Mesh& mesh,
+fem::interpolation_coords(const FiniteElement& element,
+                          const mesh::Geometry<double>& geometry,
                           std::span<const std::int32_t> cells)
 {
-  // Get mesh geometry data and the element coordinate map
-  const std::size_t gdim = mesh.geometry().dim();
-  const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
-
-  std::span<const double> x_g = mesh.geometry().x();
-  const CoordinateElement& cmap = mesh.geometry().cmap();
+  // Get geometry data and the element coordinate map
+  const std::size_t gdim = geometry.dim();
+  const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
+  std::span<const double> x_g = geometry.x();
+  const CoordinateElement& cmap = geometry.cmap();
   const std::size_t num_dofs_g = cmap.dim();
 
   // Get the interpolation points on the reference cells
@@ -71,49 +71,37 @@ fem::interpolation_coords(const FiniteElement& element, const mesh::Mesh& mesh,
 }
 //-----------------------------------------------------------------------------
 fem::nmm_interpolation_data_t fem::create_nonmatching_meshes_interpolation_data(
-    const fem::FunctionSpace& Vu, const fem::FunctionSpace& Vv,
-    std::span<const std::int32_t> cells)
+    const mesh::Geometry<double>& geometry0, const FiniteElement& element0,
+    const mesh::Mesh& mesh1, std::span<const std::int32_t> cells)
 {
-
   // Collect all the points at which values are needed to define the
   // interpolating function
-  auto element_u = Vu.element();
-  assert(element_u);
-  auto mesh = Vu.mesh();
-  assert(mesh);
-  const std::vector<double> coords_b
-      = interpolation_coords(*element_u, *mesh, cells);
-
-  namespace stdex = std::experimental;
-  using cmdspan2_t
-      = stdex::mdspan<const double, stdex::dextents<std::size_t, 2>>;
-  using mdspan2_t = stdex::mdspan<double, stdex::dextents<std::size_t, 2>>;
-  cmdspan2_t coords(coords_b.data(), 3, coords_b.size() / 3);
+  const std::vector<double> coords
+      = interpolation_coords(element0, geometry0, cells);
 
   // Transpose interpolation coords
   std::vector<double> x(coords.size());
-  mdspan2_t _x(x.data(), coords_b.size() / 3, 3);
-  for (std::size_t j = 0; j < coords.extent(1); ++j)
-    for (std::size_t i = 0; i < 3; ++i)
-      _x(j, i) = coords(i, j);
+  std::size_t num_points = coords.size() / 3;
+  for (std::size_t i = 0; i < num_points; ++i)
+    for (std::size_t j = 0; j < 3; ++j)
+      x[3 * i + j] = coords[i + j * num_points];
 
   // Determine ownership of each point
-  auto mesh_v = Vv.mesh();
-  assert(mesh_v);
-  return geometry::determine_point_ownership(*mesh_v, x);
+  return geometry::determine_point_ownership(mesh1, x);
 }
 //-----------------------------------------------------------------------------
 fem::nmm_interpolation_data_t
-fem::create_nonmatching_meshes_interpolation_data(const FunctionSpace& Vu,
-                                                  const FunctionSpace& Vv)
+fem::create_nonmatching_meshes_interpolation_data(const mesh::Mesh& mesh0,
+                                                  const FiniteElement& element0,
+                                                  const mesh::Mesh& mesh1)
 {
-  assert(Vu.mesh());
-  int tdim = Vu.mesh()->topology().dim();
-  auto cell_map = Vu.mesh()->topology().index_map(tdim);
+  int tdim = mesh0.topology().dim();
+  auto cell_map = mesh0.topology().index_map(tdim);
   assert(cell_map);
   std::int32_t num_cells = cell_map->size_local() + cell_map->num_ghosts();
   std::vector<std::int32_t> cells(num_cells, 0);
   std::iota(cells.begin(), cells.end(), 0);
-  return create_nonmatching_meshes_interpolation_data(Vu, Vv, cells);
+  return create_nonmatching_meshes_interpolation_data(mesh0.geometry(),
+                                                      element0, mesh1, cells);
 }
 //-----------------------------------------------------------------------------
diff --git a/cpp/dolfinx/fem/interpolate.h b/cpp/dolfinx/fem/interpolate.h
index dd532b0780b..c32e427fd45 100644
--- a/cpp/dolfinx/fem/interpolate.h
+++ b/cpp/dolfinx/fem/interpolate.h
@@ -25,24 +25,24 @@ namespace dolfinx::fem
 template <typename T>
 class Function;
 
-/// Compute the evaluation points in the physical space at which an
-/// expression should be computed to interpolate it in a finite element
-/// space.
+/// @brief Compute the evaluation points in the physical space at which
+/// an expression should be computed to interpolate it in a finite
+/// element space.
 ///
 /// @param[in] element The element to be interpolated into
-/// @param[in] mesh The domain
+/// @param[in] geometry Mesh geometry
 /// @param[in] cells Indices of the cells in the mesh to compute
 /// interpolation coordinates for
 /// @return The coordinates in the physical space at which to evaluate
 /// an expression. The shape is (3, num_points) and storage is row-major.
 std::vector<double> interpolation_coords(const fem::FiniteElement& element,
-                                         const mesh::Mesh& mesh,
+                                         const mesh::Geometry<double>& geometry,
                                          std::span<const std::int32_t> cells);
 
 /// Helper type for the data that can be cached to speed up repeated
 /// interpolation of discrete functions on nonmatching meshes
-using nmm_interpolation_data_t = decltype(std::function{
-    dolfinx::geometry::determine_point_ownership})::result_type;
+using nmm_interpolation_data_t
+    = decltype(std::function{geometry::determine_point_ownership})::result_type;
 
 /// Forward declaration
 template <typename T>
@@ -52,18 +52,18 @@ void interpolate(Function<T>& u, std::span<const T> f,
 
 namespace impl
 {
-/// @brief Scatter data into non-contiguous memory
+/// @brief Scatter data into non-contiguous memory.
 ///
-/// Scatter blocked data `send_values` to its
-/// corresponding src_rank and insert the data into `recv_values`.
-/// The insert location in `recv_values` is determined by `dest_ranks`.
-/// If the j-th dest rank is -1, then
-/// `recv_values[j*block_size:(j+1)*block_size]) = 0.
+/// Scatter blocked data `send_values` to its corresponding src_rank and
+/// insert the data into `recv_values`. The insert location in
+/// `recv_values` is determined by `dest_ranks`. If the j-th dest rank
+/// is -1, then `recv_values[j*block_size:(j+1)*block_size]) = 0.
 ///
 /// @param[in] comm The mpi communicator
-/// @param[in] src_ranks The rank owning the values of each row in send_values
-/// @param[in] dest_ranks List of ranks receiving data. Size of array is how
-/// many values we are receiving (not unrolled for blcok_size).
+/// @param[in] src_ranks The rank owning the values of each row in
+/// send_values
+/// @param[in] dest_ranks List of ranks receiving data. Size of array is
+/// how many values we are receiving (not unrolled for blcok_size).
 /// @param[in] send_values The values to send back to owner. Shape
 /// (src_ranks.size(), block_size). Storage is row-major.
 /// @param[in] s_shape Shape of send_values
@@ -73,11 +73,10 @@ namespace impl
 /// @note dest_ranks can contain repeated entries
 /// @note dest_ranks might contain -1 (no process owns the point)
 template <typename T>
-void scatter_values(const MPI_Comm& comm,
-                    std::span<const std::int32_t> src_ranks,
+void scatter_values(MPI_Comm comm, std::span<const std::int32_t> src_ranks,
                     std::span<const std::int32_t> dest_ranks,
                     std::span<const T> send_values,
-                    const std::array<std::size_t, 2>& s_shape,
+                    std::array<std::size_t, 2> s_shape,
                     std::span<T> recv_values)
 {
   const std::size_t block_size = s_shape[1];
@@ -916,29 +915,29 @@ void interpolate(Function<T>& u, std::span<const T> f,
   }
 }
 
-/// Generate data needed to interpolate discrete functions across
+/// @brief Generate data needed to interpolate discrete functions across
 /// different meshes.
-///
-/// @param[out] Vu The function space of the function to interpolate
-/// into
-/// @param[in] Vv The function space of the function to interpolate from
+/// @param[in] geometry0 Mesh geometry of the space to interpolate into
+/// @param[in] element0 Element of the space to interpolate into
+/// @param[in] mesh1 Mesh of the function to interpolate from
 /// @param[in] cells Indices of the cells in the destination mesh on
 /// which to interpolate. Should be the same as the list used when
 /// calling fem::interpolation_coords.
 nmm_interpolation_data_t create_nonmatching_meshes_interpolation_data(
-    const FunctionSpace& Vu, const FunctionSpace& Vv,
-    std::span<const std::int32_t> cells);
-
-/// Generate data needed to interpolate discrete functions defined on
-/// different meshes. Interpolate on all cells in the mesh.
-///
-/// @param[out] Vu The function space of the function to interpolate into
-/// @param[in] Vv The function space of the function to interpolate from
+    const mesh::Geometry<double>& geometry0, const FiniteElement& element0,
+    const mesh::Mesh& mesh1, std::span<const std::int32_t> cells);
+
+/// @brief Generate data needed to interpolate discrete functions
+/// defined on different meshes. Interpolate on all cells in the mesh.
+/// @param[in] mesh0 Mesh of the space to interpolate into
+/// @param[in] element0 Element of the space to interpolate into
+/// @param[in] mesh1 Mesh of the function to interpolate from
 nmm_interpolation_data_t
-create_nonmatching_meshes_interpolation_data(const FunctionSpace& Vu,
-                                             const FunctionSpace& Vv);
+create_nonmatching_meshes_interpolation_data(const mesh::Mesh& mesh0,
+                                             const FiniteElement& element0,
+                                             const mesh::Mesh& mesh1);
 
-/// Interpolate from one finite element Function to another one
+/// @brief Interpolate from one finite element Function to another one.
 /// @param[out] u The function to interpolate into
 /// @param[in] v The function to be interpolated
 /// @param[in] cells List of cell indices to interpolate on
diff --git a/cpp/dolfinx/geometry/BoundingBoxTree.cpp b/cpp/dolfinx/geometry/BoundingBoxTree.cpp
index 66b70677dc5..1bcbaf41fbd 100644
--- a/cpp/dolfinx/geometry/BoundingBoxTree.cpp
+++ b/cpp/dolfinx/geometry/BoundingBoxTree.cpp
@@ -20,12 +20,10 @@ using namespace dolfinx::geometry;
 namespace
 {
 //-----------------------------------------------------------------------------
-std::vector<std::int32_t> range(const mesh::Mesh& mesh, int tdim)
+std::vector<std::int32_t> range(mesh::Topology& topology, int tdim)
 {
-  // Initialize entities of given dimension if they don't exist
-  mesh.topology_mutable().create_entities(tdim);
-
-  auto map = mesh.topology().index_map(tdim);
+  topology.create_entities(tdim);
+  auto map = topology.index_map(tdim);
   assert(map);
   const std::int32_t num_entities = map->size_local() + map->num_ghosts();
   std::vector<std::int32_t> r(num_entities);
@@ -207,7 +205,8 @@ std::int32_t _build_from_point(
 //-----------------------------------------------------------------------------
 BoundingBoxTree::BoundingBoxTree(const mesh::Mesh& mesh, int tdim,
                                  double padding)
-    : BoundingBoxTree::BoundingBoxTree(mesh, tdim, range(mesh, tdim), padding)
+    : BoundingBoxTree::BoundingBoxTree(
+        mesh, tdim, range(mesh.topology_mutable(), tdim), padding)
 {
   // Do nothing
 }
diff --git a/cpp/dolfinx/geometry/utils.cpp b/cpp/dolfinx/geometry/utils.cpp
index 475fd9506dc..6314176e561 100644
--- a/cpp/dolfinx/geometry/utils.cpp
+++ b/cpp/dolfinx/geometry/utils.cpp
@@ -532,10 +532,10 @@ int geometry::compute_first_colliding_cell(
       if (norm < eps2)
         return cell;
     }
+
+    return -1;
   }
-  return -1;
 }
-
 //-------------------------------------------------------------------------------
 std::tuple<std::vector<std::int32_t>, std::vector<std::int32_t>,
            std::vector<double>, std::vector<std::int32_t>>
@@ -649,8 +649,9 @@ geometry::determine_point_ownership(const mesh::Mesh& mesh,
     cell_indicator[p / 3] = (colliding_cell >= 0) ? rank : -1;
     colliding_cells[p / 3] = colliding_cell;
   }
-  // Create neighborhood communicator in the reverse direction: send back col to
-  // requesting processes
+
+  // Create neighborhood communicator in the reverse direction: send
+  // back col to requesting processes
   MPI_Comm reverse_comm;
   MPI_Dist_graph_create_adjacent(
       comm, out_ranks.size(), out_ranks.data(), MPI_UNWEIGHTED, in_ranks.size(),
@@ -689,6 +690,7 @@ geometry::determine_point_ownership(const mesh::Mesh& mesh,
     if ((recv_ranks[i] >= 0) && (point_owners[pos] == -1))
       point_owners[pos] = recv_ranks[i];
   }
+
   // Communication is reversed again to send dest ranks to all processes
   std::swap(send_sizes, recv_sizes);
   std::swap(send_offsets, recv_offsets);
@@ -739,4 +741,5 @@ geometry::determine_point_ownership(const mesh::Mesh& mesh,
 
   return std::make_tuple(point_owners, owned_recv_ranks, owned_recv_points,
                          owned_recv_cells);
-};
\ No newline at end of file
+};
+//-------------------------------------------------------------------------------
diff --git a/cpp/dolfinx/io/ADIOS2Writers.cpp b/cpp/dolfinx/io/ADIOS2Writers.cpp
index deea238cf27..a0e85ddb548 100644
--- a/cpp/dolfinx/io/ADIOS2Writers.cpp
+++ b/cpp/dolfinx/io/ADIOS2Writers.cpp
@@ -161,7 +161,8 @@ void vtx_write_mesh(adios2::IO& io, adios2::Engine& engine,
       io, "NumberOfNodes", {adios2::LocalValueDim});
   engine.Put<std::uint32_t>(vertices, num_vertices);
 
-  const auto [vtkcells, shape] = io::extract_vtk_connectivity(mesh);
+  const auto [vtkcells, shape] = io::extract_vtk_connectivity(
+      mesh.geometry(), mesh.topology().cell_type());
 
   // Add cell metadata
   const int tdim = topology.dim();
@@ -552,7 +553,8 @@ void fides_write_mesh(adios2::IO& io, adios2::Engine& engine,
   const int tdim = topology.dim();
   const std::int32_t num_cells = topology.index_map(tdim)->size_local();
   const int num_nodes = geometry.cmap().dim();
-  const auto [cells, shape] = io::extract_vtk_connectivity(mesh);
+  const auto [cells, shape] = io::extract_vtk_connectivity(
+      mesh.geometry(), mesh.topology().cell_type());
 
   // "Put" topology data in the result in the ADIOS2 file
   adios2::Variable<std::int64_t> local_topology = define_variable<std::int64_t>(
diff --git a/cpp/dolfinx/io/VTKFile.cpp b/cpp/dolfinx/io/VTKFile.cpp
index d8d666aebef..46e126c4a11 100644
--- a/cpp/dolfinx/io/VTKFile.cpp
+++ b/cpp/dolfinx/io/VTKFile.cpp
@@ -419,7 +419,8 @@ void write_function(
   if (is_cellwise(*V0))
   {
     std::vector<std::int64_t> tmp;
-    std::tie(tmp, cshape) = io::extract_vtk_connectivity(*mesh0);
+    std::tie(tmp, cshape) = io::extract_vtk_connectivity(
+        mesh0->geometry(), mesh0->topology().cell_type());
     cells.assign(tmp.begin(), tmp.end());
     const mesh::Geometry<double>& geometry = mesh0->geometry();
     x.assign(geometry.x().begin(), geometry.x().end());
@@ -770,7 +771,8 @@ void io::VTKFile::write(const mesh::Mesh& mesh, double time)
   piece_node.append_attribute("NumberOfCells") = num_cells;
 
   // Add mesh data to "Piece" node
-  const auto [cells, cshape] = extract_vtk_connectivity(mesh);
+  const auto [cells, cshape]
+      = extract_vtk_connectivity(mesh.geometry(), mesh.topology().cell_type());
   std::array<std::size_t, 2> xshape = {geometry.x().size() / 3, 3};
   std::vector<std::uint8_t> x_ghost(xshape[0], 0);
   std::fill(std::next(x_ghost.begin(), xmap->size_local()), x_ghost.end(), 1);
diff --git a/cpp/dolfinx/io/vtk_utils.cpp b/cpp/dolfinx/io/vtk_utils.cpp
index 40b349c1dbb..024973a10e8 100644
--- a/cpp/dolfinx/io/vtk_utils.cpp
+++ b/cpp/dolfinx/io/vtk_utils.cpp
@@ -188,20 +188,18 @@ io::vtk_mesh_from_space(const fem::FunctionSpace& V)
 }
 //-----------------------------------------------------------------------------
 std::pair<std::vector<std::int64_t>, std::array<std::size_t, 2>>
-io::extract_vtk_connectivity(const mesh::Mesh& mesh)
+io::extract_vtk_connectivity(const mesh::Geometry<double>& geometry,
+                             mesh::CellType cell_type)
 {
   // Get DOLFINx to VTK permutation
   // FIXME: Use better way to get number of nodes
-  const graph::AdjacencyList<std::int32_t>& dofmap_x = mesh.geometry().dofmap();
-  const std::size_t num_nodes = mesh.geometry().cmap().dim();
-  mesh::CellType cell_type = mesh.topology().cell_type();
+  const graph::AdjacencyList<std::int32_t>& dofmap_x = geometry.dofmap();
+  const std::size_t num_nodes = geometry.cmap().dim();
   std::vector vtkmap
       = io::cells::transpose(io::cells::perm_vtk(cell_type, num_nodes));
 
   // Extract mesh 'nodes'
-  const int tdim = mesh.topology().dim();
-  const std::size_t num_cells = mesh.topology().index_map(tdim)->size_local()
-                                + mesh.topology().index_map(tdim)->num_ghosts();
+  const std::size_t num_cells = dofmap_x.num_nodes();
 
   // Build mesh connectivity
 
diff --git a/cpp/dolfinx/io/vtk_utils.h b/cpp/dolfinx/io/vtk_utils.h
index bb0fa86b289..ef7d359f964 100644
--- a/cpp/dolfinx/io/vtk_utils.h
+++ b/cpp/dolfinx/io/vtk_utils.h
@@ -20,6 +20,9 @@ class FunctionSpace;
 
 namespace mesh
 {
+enum class CellType;
+template <typename T>
+class Geometry;
 class Mesh;
 } // namespace mesh
 
@@ -44,11 +47,14 @@ std::tuple<std::vector<double>, std::array<std::size_t, 2>,
            std::vector<std::int64_t>, std::array<std::size_t, 2>>
 vtk_mesh_from_space(const fem::FunctionSpace& V);
 
-/// Extract the cell topology (connectivity) in VTK ordering for all
-/// cells the mesh. The 'topology' includes higher-order 'nodes'. The
-/// index of a 'node' corresponds to the index of DOLFINx geometry
+/// @brief Extract the cell topology (connectivity) in VTK ordering for
+/// all cells the mesh. The 'topology' includes higher-order 'nodes'.
+///
+/// The index of a 'node' corresponds to the index of DOLFINx geometry
 /// 'nodes'.
-/// @param [in] mesh The mesh
+///
+/// @param[in] geometry The mesh geometry
+/// @param[in] cell_type The cell type
 /// @return The cell topology in VTK ordering and in term of the DOLFINx
 /// geometry 'nodes'
 /// @note The indices in the return array correspond to the point
@@ -56,7 +62,7 @@ vtk_mesh_from_space(const fem::FunctionSpace& V);
 /// @note Even if the indices are local (int32), both Fides and VTX
 /// require int64 as local input
 std::pair<std::vector<std::int64_t>, std::array<std::size_t, 2>>
-extract_vtk_connectivity(const mesh::Mesh& mesh);
-
+extract_vtk_connectivity(const mesh::Geometry<double>& geometry,
+                         mesh::CellType cell_type);
 } // namespace io
 } // namespace dolfinx
diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py
index 7accac8156e..4f922b7c25e 100644
--- a/python/dolfinx/fem/function.py
+++ b/python/dolfinx/fem/function.py
@@ -375,7 +375,7 @@ def _(expr: Expression, cells: typing.Optional[np.ndarray] = None):
         except TypeError:
             # u is callable
             assert callable(u)
-            x = _cpp.fem.interpolation_coords(self._V.element, self._V.mesh._cpp_object, cells)
+            x = _cpp.fem.interpolation_coords(self._V.element, self._V.mesh.geometry, cells)
             self._cpp_object.interpolate(np.asarray(u(x), dtype=self.dtype), cells)
 
     def copy(self) -> Function:
diff --git a/python/dolfinx/wrappers/fem.cpp b/python/dolfinx/wrappers/fem.cpp
index f186fdf00f0..ae6324cfdcc 100644
--- a/python/dolfinx/wrappers/fem.cpp
+++ b/python/dolfinx/wrappers/fem.cpp
@@ -199,7 +199,7 @@ void declare_objects(py::module& m, const std::string& type)
 
             assert(self.function_space()->mesh());
             const std::vector<double> x = dolfinx::fem::interpolation_coords(
-                *element, *self.function_space()->mesh(),
+                *element, self.function_space()->mesh()->geometry(),
                 std::span(cells.data(), cells.size()));
 
             std::array<std::size_t, 2> shape = {value_size, x.size() / 3};
@@ -627,32 +627,34 @@ void fem(py::module& m)
       py::arg("integral_type"), py::arg("meshtags"));
   m.def(
       "create_nonmatching_meshes_interpolation_data",
-      [](const dolfinx::fem::FunctionSpace& Vu,
-         const dolfinx::fem::FunctionSpace& Vv)
+      [](const dolfinx::mesh::Mesh& mesh0,
+         const dolfinx::fem::FiniteElement& element0,
+         const dolfinx::mesh::Mesh& mesh1)
       {
-        assert(Vu.mesh());
-        int tdim = Vu.mesh()->topology().dim();
-        auto cell_map = Vu.mesh()->topology().index_map(tdim);
+        int tdim = mesh0.topology().dim();
+        auto cell_map = mesh0.topology().index_map(tdim);
         assert(cell_map);
         std::int32_t num_cells
             = cell_map->size_local() + cell_map->num_ghosts();
         std::vector<std::int32_t> cells(num_cells, 0);
         std::iota(cells.begin(), cells.end(), 0);
-
         return dolfinx::fem::create_nonmatching_meshes_interpolation_data(
-            Vu, Vv, std::span(cells.data(), cells.size()));
+            mesh0.geometry(), element0, mesh1,
+            std::span(cells.data(), cells.size()));
       },
-      py::arg("Vu"), py::arg("Vv"));
+      py::arg("mesh0"), py::arg("element0"), py::arg("mesh1"));
   m.def(
       "create_nonmatching_meshes_interpolation_data",
-      [](const dolfinx::fem::FunctionSpace& Vu,
-         const dolfinx::fem::FunctionSpace& Vv,
+      [](const dolfinx::mesh::Geometry<double>& geometry0,
+         const dolfinx::fem::FiniteElement& element0,
+         const dolfinx::mesh::Mesh& mesh1,
          const py::array_t<std::int32_t, py::array::c_style>& cells)
       {
         return dolfinx::fem::create_nonmatching_meshes_interpolation_data(
-            Vu, Vv, std::span(cells.data(), cells.size()));
+            geometry0, element0, mesh1, std::span(cells.data(), cells.size()));
       },
-      py::arg("Vu"), py::arg("Vv"), py::arg("cells"));
+      py::arg("geometry0"), py::arg("element0"), py::arg("mesh1"),
+      py::arg("cells"));
 
   // dolfinx::fem::FiniteElement
   py::class_<dolfinx::fem::FiniteElement,
@@ -959,11 +961,12 @@ void fem(py::module& m)
 
   m.def(
       "interpolation_coords",
-      [](const dolfinx::fem::FiniteElement& e, const dolfinx::mesh::Mesh& mesh,
+      [](const dolfinx::fem::FiniteElement& e,
+         const dolfinx::mesh::Geometry<double>& geometry,
          py::array_t<std::int32_t, py::array::c_style> cells)
       {
         std::vector<double> x = dolfinx::fem::interpolation_coords(
-            e, mesh, std::span(cells.data(), cells.size()));
+            e, geometry, std::span(cells.data(), cells.size()));
         return as_pyarray(std::move(x),
                           std::array<std::size_t, 2>{3, x.size() / 3});
       },
diff --git a/python/dolfinx/wrappers/io.cpp b/python/dolfinx/wrappers/io.cpp
index 86ea8bc429c..b5fb4531666 100644
--- a/python/dolfinx/wrappers/io.cpp
+++ b/python/dolfinx/wrappers/io.cpp
@@ -40,13 +40,14 @@ void io(py::module& m)
 
   m.def(
       "extract_vtk_connectivity",
-      [](const dolfinx::mesh::Mesh& mesh)
+      [](const dolfinx::mesh::Geometry<double>& x, dolfinx::mesh::CellType cell)
       {
-        auto [cells, shape] = dolfinx::io::extract_vtk_connectivity(mesh);
+        auto [cells, shape] = dolfinx::io::extract_vtk_connectivity(x, cell);
         return as_pyarray(std::move(cells), shape);
       },
-      py::arg("mesh"),
-      "Extract the mesh topology with VTK ordering using geometry indices");
+      py::arg("x"), py::arg("celltype"),
+      "Extract the mesh topology with VTK ordering using "
+      "geometry indices");
 
   // dolfinx::io::cell permutation functions
   m.def("perm_vtk", &dolfinx::io::cells::perm_vtk, py::arg("type"),
diff --git a/python/test/unit/fem/test_function.py b/python/test/unit/fem/test_function.py
index 646ea6acd4f..43ba3becae7 100644
--- a/python/test/unit/fem/test_function.py
+++ b/python/test/unit/fem/test_function.py
@@ -197,8 +197,9 @@ def f(x):
     # Interpolate 3D->2D
     u1 = Function(V1)
     u1.interpolate(u0, nmm_interpolation_data=create_nonmatching_meshes_interpolation_data(
-        u1.function_space._cpp_object,
-        u0.function_space._cpp_object))
+        u1.function_space.mesh._cpp_object,
+        u1.function_space.element,
+        u0.function_space.mesh._cpp_object))
     u1.x.scatter_forward()
 
     # Exact interpolation on 2D mesh
@@ -211,8 +212,9 @@ def f(x):
     # Interpolate 2D->3D
     u0_2 = Function(V0)
     u0_2.interpolate(u1, nmm_interpolation_data=create_nonmatching_meshes_interpolation_data(
-        u0_2.function_space._cpp_object,
-        u1.function_space._cpp_object))
+        u0_2.function_space.mesh._cpp_object,
+        u0_2.function_space.element,
+        u1.function_space.mesh._cpp_object))
 
     # Check that function values over facets of 3D mesh of the twice interpolated property is preserved
     def locate_bottom_facets(x):