Introduction of transfer operator assembly for geometry multigrid operations

cpp/dolfinx/CMakeLists.txt

@@ -18,6 +18,7 @@ set(DOLFINX_DIRS
     mesh
     nls
     refinement
+    transfer
 )
 
 # Add source to dolfinx target, and get sets of header files

cpp/dolfinx/refinement/CMakeLists.txt

@@ -11,4 +11,5 @@ set(HEADERS_refinement
 target_sources(
   dolfinx PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/plaza.cpp
                   ${CMAKE_CURRENT_SOURCE_DIR}/utils.cpp
+                  ${CMAKE_CURRENT_SOURCE_DIR}/refine.cpp
 )

cpp/dolfinx/refinement/refine.cpp

@@ -0,0 +1,25 @@
+// Copyright (C) 2024 Paul T. Kühner
+//
+// This file is part of DOLFINx (https://www.fenicsproject.org)
+//
+// SPDX-License-Identifier:    LGPL-3.0-or-later
+
+#include "refine.h"
+
+namespace dolfinx::refinement
+{
+
+graph::AdjacencyList<std::int32_t> maintain_coarse_partitioner(
+    MPI_Comm comm, int, const std::vector<mesh::CellType>& cell_types,
+    const std::vector<std::span<const std::int64_t>>& cell_topology)
+{
+  int mpi_rank = MPI::rank(comm);
+  int num_cell_vertices = mesh::num_cell_vertices(cell_types.front());
+  std::int32_t num_cells = cell_topology.front().size() / num_cell_vertices;
+  std::vector<std::int32_t> destinations(num_cells, mpi_rank);
+  std::vector<std::int32_t> dest_offsets(num_cells + 1);
+  std::iota(dest_offsets.begin(), dest_offsets.end(), 0);
+  return graph::AdjacencyList(std::move(destinations), std::move(dest_offsets));
+}
+
+} // namespace dolfinx::refinement

cpp/dolfinx/refinement/refine.h

@@ -1,11 +1,12 @@
-// Copyright (C) 2010-2023 Garth N. Wells
+// Copyright (C) 2010-2024 Garth N. Wells and Paul T. Kühner
 //
 // This file is part of DOLFINx (https://www.fenicsproject.org)
 //
 // SPDX-License-Identifier:    LGPL-3.0-or-later
 
 #pragma once
 
+#include "dolfinx/graph/AdjacencyList.h"
 #include "dolfinx/mesh/Mesh.h"
 #include "dolfinx/mesh/Topology.h"
 #include "dolfinx/mesh/cell_types.h"
@@ -19,44 +20,28 @@
 
 namespace dolfinx::refinement
 {
-namespace impl
-{
-/// @brief create the refined mesh by optionally redistributing it
-template <std::floating_point T>
-mesh::Mesh<T>
-create_refined_mesh(const mesh::Mesh<T>& mesh,
-                    const graph::AdjacencyList<std::int64_t>& cell_adj,
-                    std::span<const T> new_vertex_coords,
-                    std::array<std::size_t, 2> xshape, bool redistribute,
-                    mesh::GhostMode ghost_mode)
-{
-  if (dolfinx::MPI::size(mesh.comm()) == 1)
-  {
-    // No parallel construction necessary, i.e. redistribute also has no
-    // effect
-    return mesh::create_mesh(mesh.comm(), cell_adj.array(),
-                             mesh.geometry().cmap(), new_vertex_coords, xshape,
-                             ghost_mode);
-  }
-  else
-  {
-    // Let partition handle the parallel construction of the mesh
-    return partition<T>(mesh, cell_adj, new_vertex_coords, xshape, redistribute,
-                        ghost_mode);
-  }
-}
-} // namespace impl
+
+/// @brief Partitioner that maintains the regional destribution for a refined
+/// mesh. Meaning, process local data is just maintained and no redistribution
+/// happens.
+///
+/// @param[in] comm MPI Communicator
+/// @param[in] nparts Number of partitions (input has no effect)
+/// @param[in] cell_types Cell types in the mesh
+/// @param[in] cells Lists of cells of each cell type.
+/// @return Destination ranks for each cell on this process
+graph::AdjacencyList<std::int32_t> maintain_coarse_partitioner(
+    MPI_Comm comm, int nparts, const std::vector<mesh::CellType>& cell_types,
+    const std::vector<std::span<const std::int64_t>>& cell_topology);
 
 /// @brief Refine with markers, optionally redistributing, and
 /// optionally calculating the parent-child relationships.
 ///
-/// @param[in] mesh Input mesh to be refined
+/// @param[in] mesh Input mesh to be refined.
 /// @param[in] edges Optional indices of the edges that should be split
 /// by this refinement, if optional is not set, a uniform refinement is
 /// performed, same behavior as passing a list of all indices.
-/// @param[in] redistribute Flag to call the Mesh Partitioner to
-/// redistribute after refinement.
-/// @param[in] ghost_mode Ghost mode of the refined mesh.
+/// @param[in] partitioner partitioner to be used for the refined mesh
 /// @param[in] option Control the computation of parent facets, parent
 /// cells. If an option is unselected, an empty list is returned.
 /// @return New Mesh and optional parent cell index, parent facet
@@ -65,26 +50,26 @@ template <std::floating_point T>
 std::tuple<mesh::Mesh<T>, std::optional<std::vector<std::int32_t>>,
            std::optional<std::vector<std::int8_t>>>
 refine(const mesh::Mesh<T>& mesh,
-       std::optional<std::span<const std::int32_t>> edges, bool redistribute,
-       mesh::GhostMode ghost_mode = mesh::GhostMode::shared_facet,
+       std::optional<std::span<const std::int32_t>> edges,
+       mesh::CellPartitionFunction partitioner = maintain_coarse_partitioner,
        Option option = Option::none)
 {
   auto topology = mesh.topology();
   assert(topology);
-
-  mesh::CellType cell_t = topology->cell_type();
-  if (!mesh::is_simplex(cell_t))
+  if (!mesh::is_simplex(topology->cell_type()))
     throw std::runtime_error("Refinement only defined for simplices");
+
   bool oned = topology->cell_type() == mesh::CellType::interval;
   auto [cell_adj, new_vertex_coords, xshape, parent_cell, parent_facet]
       = oned ? interval::compute_refinement_data(mesh, edges, option)
              : plaza::compute_refinement_data(mesh, edges, option);
 
-  mesh::Mesh<T> refined_mesh = impl::create_refined_mesh(
-      mesh, cell_adj, std::span<const T>(new_vertex_coords), xshape,
-      redistribute, ghost_mode);
+  mesh::Mesh<T> refined_mesh
+      = mesh::create_mesh(mesh.comm(), mesh.comm(), std::move(cell_adj).array(),
+                          mesh.geometry().cmap(), mesh.comm(),
+                          std::move(new_vertex_coords), xshape, partitioner);
 
-  // Report the number of refined cellse
+  // Report the number of refined cells
   const int D = topology->dim();
   const std::int64_t n0 = topology->index_map(D)->size_global();
   const std::int64_t n1 = refined_mesh.topology()->index_map(D)->size_global();

cpp/dolfinx/refinement/utils.h

@@ -264,51 +264,6 @@ create_new_vertices(MPI_Comm comm,
           xshape};
 }
 
-/// Use vertex and topology data to partition new mesh across
-/// processes
-/// @param[in] old_mesh
-/// @param[in] cell_topology Topology of cells, (vertex indices)
-/// @param[in] new_coords New coordinates, row-major storage
-/// @param[in] xshape The shape of `new_coords`
-/// @param[in] redistribute Call graph partitioner if true
-/// @param[in] ghost_mode None or shared_facet
-/// @return New mesh
-template <std::floating_point T>
-mesh::Mesh<T> partition(const mesh::Mesh<T>& old_mesh,
-                        const graph::AdjacencyList<std::int64_t>& cell_topology,
-                        std::span<const T> new_coords,
-                        std::array<std::size_t, 2> xshape, bool redistribute,
-                        mesh::GhostMode ghost_mode)
-{
-  if (redistribute)
-  {
-    return mesh::create_mesh(old_mesh.comm(), cell_topology.array(),
-                             old_mesh.geometry().cmap(), new_coords, xshape,
-                             ghost_mode);
-  }
-  else
-  {
-    auto partitioner
-        = [](MPI_Comm comm, int, const std::vector<mesh::CellType>& cell_types,
-             const std::vector<std::span<const std::int64_t>>& cell_topology)
-    {
-      std::int32_t mpi_rank = MPI::rank(comm);
-      std::int32_t num_cell_vertices
-          = mesh::num_cell_vertices(cell_types.front());
-      std::int32_t num_cells = cell_topology.front().size() / num_cell_vertices;
-      std::vector<std::int32_t> destinations(num_cells, mpi_rank);
-      std::vector<std::int32_t> dest_offsets(num_cells + 1);
-      std::iota(dest_offsets.begin(), dest_offsets.end(), 0);
-      return graph::AdjacencyList(std::move(destinations),
-                                  std::move(dest_offsets));
-    };
-
-    return mesh::create_mesh(old_mesh.comm(), old_mesh.comm(),
-                             cell_topology.array(), old_mesh.geometry().cmap(),
-                             old_mesh.comm(), new_coords, xshape, partitioner);
-  }
-}
-
 /// @brief Given an index map, add "n" extra indices at the end of local range
 ///
 /// @note The returned global indices (local and ghosts) are adjust for the

cpp/dolfinx/transfer/CMakeLists.txt

@@ -0,0 +1,4 @@
+set(HEADERS_transfer
+    ${CMAKE_CURRENT_SOURCE_DIR}/transfer_matrix.h
+    PARENT_SCOPE
+)