From 22f747152393081fb3e6a44b4914ac579225c221 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Tue, 11 Mar 2025 13:34:09 +0100 Subject: [PATCH 01/36] Sketch out strategy and copy maintain_coarse_partitioner (for non ghosted case) --- cpp/dolfinx/refinement/refine.h | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index 430bcc495bc..618204022b9 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -84,4 +84,34 @@ refine(const mesh::Mesh& mesh, return {std::move(mesh1), std::move(parent_cell), std::move(parent_facet)}; } + +template +mesh::CellPartitionFunction +create_identity_partitioner(const mesh::Mesh& parent_mesh, + std::vector parent_cell) +{ + // for every non ghost index maintain_coarse_partitioner below should be doing + // the correct thing. ghost cells need to given the rank, which previously the + // COARSE cell was assigned. These ranks can be computed for parent_mesh with + // compute_destination_ranks and the parent_cell[i] is the cell idx in + // parent_mesh to access it. + + return + [&](MPI_Comm comm, int /*nparts*/, std::vector cell_types, + std::vector> cells) + -> graph::AdjacencyList + { + int num_cell_vertices = mesh::num_cell_vertices(cell_types.front()); + std::int32_t num_cells = cells.front().size() / num_cell_vertices; + std::vector destinations(num_cells, dolfinx::MPI::rank(comm)); + // TODO: iterate over previous ghost cells here and replace destination + // ranks with previous coarse ranks, i.e. do not change ownership of ghosts + // cells to local process + std::vector 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 From 6160a96a61231f32acb1d9922d6740c531237ce5 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Tue, 11 Mar 2025 14:53:34 +0100 Subject: [PATCH 02/36] A first version --- cpp/dolfinx/refinement/refine.h | 229 +++++++++++++++++++++++++++++--- 1 file changed, 214 insertions(+), 15 deletions(-) diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index 618204022b9..bab24157164 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -6,12 +6,14 @@ #pragma once +#include "common/MPI.h" #include "dolfinx/graph/AdjacencyList.h" #include "dolfinx/mesh/Mesh.h" #include "dolfinx/mesh/Topology.h" #include "dolfinx/mesh/cell_types.h" #include "dolfinx/mesh/utils.h" #include "interval.h" +#include "mesh/graphbuild.h" #include "plaza.h" #include #include @@ -19,6 +21,176 @@ #include #include +// TODO: Remove once works +namespace +{ +/// @todo Is it un-documented that the owning rank must come first in +/// reach list of edges? +/// +/// @param[in] comm The communicator +/// @param[in] graph Graph, using global indices for graph edges +/// @param[in] node_disp The distribution of graph nodes across MPI +/// ranks. The global index `gidx` of local index `lidx` is `lidx + +/// node_disp[my_rank]`. +/// @param[in] part The destination rank for owned nodes, i.e. `dest[i]` +/// is the destination of the node with local index `i`. +/// @return Destination ranks for each local node. +template +graph::AdjacencyList compute_destination_ranks( + MPI_Comm comm, const graph::AdjacencyList& graph, + const std::vector& node_disp, const std::vector& part) +{ + common::Timer timer("Extend graph destination ranks for halo"); + + const int rank = dolfinx::MPI::rank(comm); + const std::int64_t range0 = node_disp[rank]; + const std::int64_t range1 = node_disp[rank + 1]; + assert(static_cast(range1 - range0) == graph.num_nodes()); + + // Wherever an owned 'node' goes, so must the nodes connected to it by + // an edge ('node1'). Task is to let the owner of node1 know the extra + // ranks that it needs to send node1 to. + std::vector> node_to_dest; + for (int node0 = 0; node0 < graph.num_nodes(); ++node0) + { + // Wherever 'node' goes to, so must the attached 'node1' + for (auto node1 : graph.links(node0)) + { + if (node1 < range0 or node1 >= range1) + { + auto it = std::ranges::upper_bound(node_disp, node1); + int remote_rank = std::distance(node_disp.begin(), it) - 1; + node_to_dest.push_back( + {remote_rank, node1, static_cast(part[node0])}); + } + else + node_to_dest.push_back( + {rank, node1, static_cast(part[node0])}); + } + } + + std::ranges::sort(node_to_dest); + auto [unique_end, range_end] = std::ranges::unique(node_to_dest); + node_to_dest.erase(unique_end, range_end); + + // Build send data and buffer + std::vector dest, send_sizes; + std::vector send_buffer; + { + auto it = node_to_dest.begin(); + while (it != node_to_dest.end()) + { + // Current destination rank + dest.push_back(it->front()); + + // Find iterator to next destination rank and pack send data + auto it1 + = std::find_if(it, node_to_dest.end(), [r0 = dest.back()](auto& idx) + { return idx[0] != r0; }); + send_sizes.push_back(2 * std::distance(it, it1)); + for (auto itx = it; itx != it1; ++itx) + { + send_buffer.push_back(itx->at(1)); + send_buffer.push_back(itx->at(2)); + } + + it = it1; + } + } + + // Prepare send displacements + std::vector send_disp(send_sizes.size() + 1, 0); + std::partial_sum(send_sizes.begin(), send_sizes.end(), + std::next(send_disp.begin())); + + // Discover src ranks. ParMETIS/KaHIP are not scalable (holding an + // array of size equal to the comm size), so no extra harm in using + // non-scalable neighbourhood detection (which might be faster for + // small rank counts). + const std::vector src + = dolfinx::MPI::compute_graph_edges_pcx(comm, dest); + + // Create neighbourhood communicator + MPI_Comm neigh_comm; + MPI_Dist_graph_create_adjacent(comm, src.size(), src.data(), MPI_UNWEIGHTED, + dest.size(), dest.data(), MPI_UNWEIGHTED, + MPI_INFO_NULL, false, &neigh_comm); + + // Determine receives sizes + std::vector recv_sizes(dest.size()); + send_sizes.reserve(1); + recv_sizes.reserve(1); + MPI_Neighbor_alltoall(send_sizes.data(), 1, MPI_INT, recv_sizes.data(), 1, + MPI_INT, neigh_comm); + + // Prepare receive displacements + std::vector recv_disp(recv_sizes.size() + 1, 0); + std::partial_sum(recv_sizes.begin(), recv_sizes.end(), + std::next(recv_disp.begin())); + + // Send/receive data + std::vector recv_buffer(recv_disp.back()); + MPI_Neighbor_alltoallv(send_buffer.data(), send_sizes.data(), + send_disp.data(), MPI_INT64_T, recv_buffer.data(), + recv_sizes.data(), recv_disp.data(), MPI_INT64_T, + neigh_comm); + MPI_Comm_free(&neigh_comm); + + // Prepare (local node index, destination rank) array. Add local data, + // then add the received data, and the make unique. + std::vector> local_node_to_dest; + for (auto d : part) + { + local_node_to_dest.push_back( + {static_cast(local_node_to_dest.size()), static_cast(d)}); + } + for (std::size_t i = 0; i < recv_buffer.size(); i += 2) + { + std::int64_t idx = recv_buffer[i]; + int d = recv_buffer[i + 1]; + assert(idx >= range0 and idx < range1); + std::int32_t idx_local = idx - range0; + local_node_to_dest.push_back({idx_local, d}); + } + + { + std::ranges::sort(local_node_to_dest); + auto [unique_end, range_end] = std::ranges::unique(local_node_to_dest); + local_node_to_dest.erase(unique_end, range_end); + } + // Compute offsets + std::vector offsets(graph.num_nodes() + 1, 0); + { + std::vector num_dests(graph.num_nodes(), 0); + for (auto x : local_node_to_dest) + ++num_dests[x[0]]; + std::partial_sum(num_dests.begin(), num_dests.end(), + std::next(offsets.begin())); + } + + // Fill data array + std::vector data(offsets.back()); + { + std::vector pos = offsets; + for (auto [x0, x1] : local_node_to_dest) + data[pos[x0]++] = x1; + } + + graph::AdjacencyList g(std::move(data), std::move(offsets)); + + // Make sure the owning rank comes first for each node + for (std::int32_t i = 0; i < g.num_nodes(); ++i) + { + auto d = g.links(i); + auto it = std::find(d.begin(), d.end(), part[i]); + assert(it != d.end()); + std::iter_swap(d.begin(), it); + } + + return g; +} +} // namespace + namespace dolfinx::refinement { /// @brief Refine a mesh with markers. @@ -56,8 +228,7 @@ std::tuple, std::optional>, std::optional>> refine(const mesh::Mesh& mesh, std::optional> edges, - const mesh::CellPartitionFunction& partitioner - = mesh::create_cell_partitioner(mesh::GhostMode::none), + const mesh::CellPartitionFunction& partitioner = nullptr, Option option = Option::none) { auto topology = mesh.topology(); @@ -70,6 +241,11 @@ refine(const mesh::Mesh& mesh, ? interval::compute_refinement_data(mesh, edges, option) : plaza::compute_refinement_data(mesh, edges, option); + if (!partitioner) + { + partitioner = create_identity_partitioner(mesh, parent_cell); + } + mesh::Mesh mesh1 = mesh::create_mesh( mesh.comm(), mesh.comm(), cell_adj.array(), mesh.geometry().cmap(), mesh.comm(), new_vertex_coords, xshape, partitioner); @@ -88,29 +264,52 @@ refine(const mesh::Mesh& mesh, template mesh::CellPartitionFunction create_identity_partitioner(const mesh::Mesh& parent_mesh, - std::vector parent_cell) + std::span parent_cell) { - // for every non ghost index maintain_coarse_partitioner below should be doing - // the correct thing. ghost cells need to given the rank, which previously the - // COARSE cell was assigned. These ranks can be computed for parent_mesh with - // compute_destination_ranks and the parent_cell[i] is the cell idx in - // parent_mesh to access it. + // TODO: optimize for non ghosted mesh? return [&](MPI_Comm comm, int /*nparts*/, std::vector cell_types, std::vector> cells) -> graph::AdjacencyList { + auto parent_top = parent_mesh.topology(); + auto parent_cell_im = parent_top.index_map(parent_top.dim()); + int num_cell_vertices = mesh::num_cell_vertices(cell_types.front()); std::int32_t num_cells = cells.front().size() / num_cell_vertices; - std::vector destinations(num_cells, dolfinx::MPI::rank(comm)); - // TODO: iterate over previous ghost cells here and replace destination - // ranks with previous coarse ranks, i.e. do not change ownership of ghosts - // cells to local process + std::vector destinations(num_cells); + std::vector dest_offsets(num_cells + 1); - std::iota(dest_offsets.begin(), dest_offsets.end(), 0); - return graph::AdjacencyList(std::move(destinations), - std::move(dest_offsets)); + int rank = dolfinx::MPI::rank(comm); + for (std::int32_t i = 0; i < destinations.size(); i++) + { + bool ghost_parent_cell = parent_cell[i] > parent_cell_im.size_local(); + if (ghost_parent_cell) + { + destinations[i] + = parent_cell_im + .owners()[parent_cell[i] - parent_cell_im.size_local()]; + } + else + { + destinations[i] = rank; + } + } + + auto dual_graph = mesh::build_dual_graph(comm, cell_types, cells); + std::vector node_disp; + node_disp = std::vector(MPI::size(comm) + 1, 0); + MPI_Allgather(&cells.front().size(), 1, dolfinx::MPI::mpi_t, + node_disp.data() + 1, 1, dolfinx::MPI::mpi_t, + comm); + std::partial_sum(node_disp.begin(), node_disp.end(), node_disp.begin()); + + return compute_destination_ranks(comm, dual_graph, node_disp, destinations); + + // std::iota(dest_offsets.begin(), dest_offsets.end(), 0); + // return graph::AdjacencyList(std::move(destinations), + // std::move(dest_offsets)); }; } From c7f6e321f05490014c4977df64631db813e601c1 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Tue, 11 Mar 2025 16:49:00 +0100 Subject: [PATCH 03/36] Compiling... --- cpp/dolfinx/refinement/refine.h | 116 +++++++++++++++++--------------- 1 file changed, 60 insertions(+), 56 deletions(-) diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index bab24157164..f982c76726b 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -6,14 +6,14 @@ #pragma once -#include "common/MPI.h" +#include "dolfinx/common/MPI.h" #include "dolfinx/graph/AdjacencyList.h" #include "dolfinx/mesh/Mesh.h" #include "dolfinx/mesh/Topology.h" #include "dolfinx/mesh/cell_types.h" +#include "dolfinx/mesh/graphbuild.h" #include "dolfinx/mesh/utils.h" #include "interval.h" -#include "mesh/graphbuild.h" #include "plaza.h" #include #include @@ -22,6 +22,7 @@ #include // TODO: Remove once works +namespace dolfinx { namespace { /// @todo Is it un-documented that the owning rank must come first in @@ -190,9 +191,63 @@ graph::AdjacencyList compute_destination_ranks( return g; } } // namespace +} namespace dolfinx::refinement { + +template +mesh::CellPartitionFunction +create_identity_partitioner(const mesh::Mesh& parent_mesh, + std::span parent_cell) +{ + // TODO: optimize for non ghosted mesh? + + return + [&](MPI_Comm comm, int /*nparts*/, std::vector cell_types, + std::vector> cells) + -> graph::AdjacencyList + { + auto parent_top = parent_mesh.topology(); + auto parent_cell_im = parent_top->index_map(parent_top->dim()); + + int num_cell_vertices = mesh::num_cell_vertices(cell_types.front()); + std::int32_t num_cells = cells.front().size() / num_cell_vertices; + std::vector destinations(num_cells); + + std::vector dest_offsets(num_cells + 1); + int rank = dolfinx::MPI::rank(comm); + for (std::int32_t i = 0; i < destinations.size(); i++) + { + bool ghost_parent_cell = parent_cell[i] > parent_cell_im->size_local(); + if (ghost_parent_cell) + { + destinations[i] + = parent_cell_im->owners()[parent_cell[i] - parent_cell_im->size_local()]; + } + else + { + destinations[i] = rank; + } + } + + auto dual_graph = mesh::build_dual_graph(comm, cell_types, cells); + std::vector node_disp; + node_disp = std::vector(MPI::size(comm) + 1, 0); + std::int32_t local_size = cells.front().size(); + MPI_Allgather(&local_size, 1, dolfinx::MPI::mpi_t, + node_disp.data() + 1, 1, dolfinx::MPI::mpi_t, + comm); + std::partial_sum(node_disp.begin(), node_disp.end(), node_disp.begin()); + + return compute_destination_ranks(comm, dual_graph, node_disp, destinations); + + // std::iota(dest_offsets.begin(), dest_offsets.end(), 0); + // return graph::AdjacencyList(std::move(destinations), + // std::move(dest_offsets)); + }; +} + /// @brief Refine a mesh with markers. /// /// The refined mesh can be optionally re-partitioned across processes. @@ -228,7 +283,7 @@ std::tuple, std::optional>, std::optional>> refine(const mesh::Mesh& mesh, std::optional> edges, - const mesh::CellPartitionFunction& partitioner = nullptr, + mesh::CellPartitionFunction partitioner = nullptr, Option option = Option::none) { auto topology = mesh.topology(); @@ -243,7 +298,8 @@ refine(const mesh::Mesh& mesh, if (!partitioner) { - partitioner = create_identity_partitioner(mesh, parent_cell); + assert(parent_cell); + partitioner = create_identity_partitioner(mesh, parent_cell.value()); } mesh::Mesh mesh1 = mesh::create_mesh( @@ -261,56 +317,4 @@ refine(const mesh::Mesh& mesh, return {std::move(mesh1), std::move(parent_cell), std::move(parent_facet)}; } -template -mesh::CellPartitionFunction -create_identity_partitioner(const mesh::Mesh& parent_mesh, - std::span parent_cell) -{ - // TODO: optimize for non ghosted mesh? - - return - [&](MPI_Comm comm, int /*nparts*/, std::vector cell_types, - std::vector> cells) - -> graph::AdjacencyList - { - auto parent_top = parent_mesh.topology(); - auto parent_cell_im = parent_top.index_map(parent_top.dim()); - - int num_cell_vertices = mesh::num_cell_vertices(cell_types.front()); - std::int32_t num_cells = cells.front().size() / num_cell_vertices; - std::vector destinations(num_cells); - - std::vector dest_offsets(num_cells + 1); - int rank = dolfinx::MPI::rank(comm); - for (std::int32_t i = 0; i < destinations.size(); i++) - { - bool ghost_parent_cell = parent_cell[i] > parent_cell_im.size_local(); - if (ghost_parent_cell) - { - destinations[i] - = parent_cell_im - .owners()[parent_cell[i] - parent_cell_im.size_local()]; - } - else - { - destinations[i] = rank; - } - } - - auto dual_graph = mesh::build_dual_graph(comm, cell_types, cells); - std::vector node_disp; - node_disp = std::vector(MPI::size(comm) + 1, 0); - MPI_Allgather(&cells.front().size(), 1, dolfinx::MPI::mpi_t, - node_disp.data() + 1, 1, dolfinx::MPI::mpi_t, - comm); - std::partial_sum(node_disp.begin(), node_disp.end(), node_disp.begin()); - - return compute_destination_ranks(comm, dual_graph, node_disp, destinations); - - // std::iota(dest_offsets.begin(), dest_offsets.end(), 0); - // return graph::AdjacencyList(std::move(destinations), - // std::move(dest_offsets)); - }; -} - } // namespace dolfinx::refinement From c4afbcb059c3adb8349dcf0562037f7f5bbef6da Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Tue, 11 Mar 2025 21:58:53 +0100 Subject: [PATCH 04/36] Fix: python layer default value of partitioner does not align with cpp layer --- python/dolfinx/mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index 46385f0d345..cc11abd7db5 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -518,7 +518,7 @@ def transfer_meshtag( def refine( msh: Mesh, edges: typing.Optional[np.ndarray] = None, - partitioner: typing.Optional[typing.Callable] = create_cell_partitioner(GhostMode.none), + partitioner: typing.Optional[typing.Callable] = None, option: RefinementOption = RefinementOption.none, ) -> tuple[Mesh, npt.NDArray[np.int32], npt.NDArray[np.int8]]: """Refine a mesh. From d8650569577f49de782e933824bc2943b425fe01 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Tue, 11 Mar 2025 22:10:26 +0100 Subject: [PATCH 05/36] Debug --- cpp/dolfinx/refinement/refine.h | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index f982c76726b..22e9ed49c5d 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -22,7 +22,8 @@ #include // TODO: Remove once works -namespace dolfinx { +namespace dolfinx +{ namespace { /// @todo Is it un-documented that the owning rank must come first in @@ -191,7 +192,7 @@ graph::AdjacencyList compute_destination_ranks( return g; } } // namespace -} +} // namespace dolfinx namespace dolfinx::refinement { @@ -203,10 +204,11 @@ create_identity_partitioner(const mesh::Mesh& parent_mesh, { // TODO: optimize for non ghosted mesh? - return - [&](MPI_Comm comm, int /*nparts*/, std::vector cell_types, - std::vector> cells) - -> graph::AdjacencyList + return [&parent_mesh, + parent_cell](MPI_Comm comm, int /*nparts*/, + std::vector cell_types, + std::vector> cells) + -> graph::AdjacencyList { auto parent_top = parent_mesh.topology(); auto parent_cell_im = parent_top->index_map(parent_top->dim()); @@ -215,7 +217,6 @@ create_identity_partitioner(const mesh::Mesh& parent_mesh, std::int32_t num_cells = cells.front().size() / num_cell_vertices; std::vector destinations(num_cells); - std::vector dest_offsets(num_cells + 1); int rank = dolfinx::MPI::rank(comm); for (std::int32_t i = 0; i < destinations.size(); i++) { @@ -223,7 +224,8 @@ create_identity_partitioner(const mesh::Mesh& parent_mesh, if (ghost_parent_cell) { destinations[i] - = parent_cell_im->owners()[parent_cell[i] - parent_cell_im->size_local()]; + = parent_cell_im + ->owners()[parent_cell[i] - parent_cell_im->size_local()]; } else { @@ -234,7 +236,7 @@ create_identity_partitioner(const mesh::Mesh& parent_mesh, auto dual_graph = mesh::build_dual_graph(comm, cell_types, cells); std::vector node_disp; node_disp = std::vector(MPI::size(comm) + 1, 0); - std::int32_t local_size = cells.front().size(); + std::int32_t local_size = dual_graph.num_nodes(); MPI_Allgather(&local_size, 1, dolfinx::MPI::mpi_t, node_disp.data() + 1, 1, dolfinx::MPI::mpi_t, comm); From 9b8314b34317b5d9ed26abc111127643e2b7c7ea Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Tue, 11 Mar 2025 22:28:42 +0100 Subject: [PATCH 06/36] Move compute_destination_ranks out of anonymous namespace --- cpp/dolfinx/graph/partitioners.cpp | 10 +- cpp/dolfinx/graph/partitioners.h | 6 + cpp/dolfinx/refinement/refine.h | 174 +---------------------------- 3 files changed, 11 insertions(+), 179 deletions(-) diff --git a/cpp/dolfinx/graph/partitioners.cpp b/cpp/dolfinx/graph/partitioners.cpp index fade5e160cf..e86b6ad6f2b 100644 --- a/cpp/dolfinx/graph/partitioners.cpp +++ b/cpp/dolfinx/graph/partitioners.cpp @@ -35,8 +35,6 @@ extern "C" using namespace dolfinx; -namespace -{ /// @todo Is it un-documented that the owning rank must come first in /// reach list of edges? /// @@ -49,7 +47,7 @@ namespace /// is the destination of the node with local index `i`. /// @return Destination ranks for each local node. template -graph::AdjacencyList compute_destination_ranks( +graph::AdjacencyList dolfinx::graph::compute_destination_ranks( MPI_Comm comm, const graph::AdjacencyList& graph, const std::vector& node_disp, const std::vector& part) { @@ -202,6 +200,7 @@ graph::AdjacencyList compute_destination_ranks( return g; } + //----------------------------------------------------------------------------- #ifdef HAS_PARMETIS template @@ -304,7 +303,6 @@ std::vector refine(MPI_Comm comm, const graph::AdjacencyList& adj_graph) //----------------------------------------------------------------------------- } #endif -} // namespace //----------------------------------------------------------------------------- #ifdef HAS_PTSCOTCH @@ -587,7 +585,7 @@ graph::partition_fn graph::parmetis::partitioner(double imbalance, { // FIXME: Is it implicit that the first entry is the owner? graph::AdjacencyList dest - = compute_destination_ranks(pcomm, graph, node_disp, part); + = graph::compute_destination_ranks(pcomm, graph, node_disp, part); if (split_comm) MPI_Comm_free(&pcomm); return dest; @@ -653,7 +651,7 @@ graph::partition_fn graph::kahip::partitioner(int mode, int seed, timer2.stop(); if (ghosting) - return compute_destination_ranks(comm, graph, node_disp, part); + return graph::compute_destination_ranks(comm, graph, node_disp, part); else { return regular_adjacency_list(std::vector(part.begin(), part.end()), diff --git a/cpp/dolfinx/graph/partitioners.h b/cpp/dolfinx/graph/partitioners.h index aba51a8430a..29aa62af97f 100644 --- a/cpp/dolfinx/graph/partitioners.h +++ b/cpp/dolfinx/graph/partitioners.h @@ -11,6 +11,12 @@ namespace dolfinx::graph { + +template +graph::AdjacencyList compute_destination_ranks( + MPI_Comm comm, const graph::AdjacencyList& graph, + const std::vector& node_disp, const std::vector& part); + namespace scotch { #ifdef HAS_PTSCOTCH diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index 22e9ed49c5d..39cd82045a1 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -8,6 +8,7 @@ #include "dolfinx/common/MPI.h" #include "dolfinx/graph/AdjacencyList.h" +#include "dolfinx/graph/partitioners.h" #include "dolfinx/mesh/Mesh.h" #include "dolfinx/mesh/Topology.h" #include "dolfinx/mesh/cell_types.h" @@ -21,179 +22,6 @@ #include #include -// TODO: Remove once works -namespace dolfinx -{ -namespace -{ -/// @todo Is it un-documented that the owning rank must come first in -/// reach list of edges? -/// -/// @param[in] comm The communicator -/// @param[in] graph Graph, using global indices for graph edges -/// @param[in] node_disp The distribution of graph nodes across MPI -/// ranks. The global index `gidx` of local index `lidx` is `lidx + -/// node_disp[my_rank]`. -/// @param[in] part The destination rank for owned nodes, i.e. `dest[i]` -/// is the destination of the node with local index `i`. -/// @return Destination ranks for each local node. -template -graph::AdjacencyList compute_destination_ranks( - MPI_Comm comm, const graph::AdjacencyList& graph, - const std::vector& node_disp, const std::vector& part) -{ - common::Timer timer("Extend graph destination ranks for halo"); - - const int rank = dolfinx::MPI::rank(comm); - const std::int64_t range0 = node_disp[rank]; - const std::int64_t range1 = node_disp[rank + 1]; - assert(static_cast(range1 - range0) == graph.num_nodes()); - - // Wherever an owned 'node' goes, so must the nodes connected to it by - // an edge ('node1'). Task is to let the owner of node1 know the extra - // ranks that it needs to send node1 to. - std::vector> node_to_dest; - for (int node0 = 0; node0 < graph.num_nodes(); ++node0) - { - // Wherever 'node' goes to, so must the attached 'node1' - for (auto node1 : graph.links(node0)) - { - if (node1 < range0 or node1 >= range1) - { - auto it = std::ranges::upper_bound(node_disp, node1); - int remote_rank = std::distance(node_disp.begin(), it) - 1; - node_to_dest.push_back( - {remote_rank, node1, static_cast(part[node0])}); - } - else - node_to_dest.push_back( - {rank, node1, static_cast(part[node0])}); - } - } - - std::ranges::sort(node_to_dest); - auto [unique_end, range_end] = std::ranges::unique(node_to_dest); - node_to_dest.erase(unique_end, range_end); - - // Build send data and buffer - std::vector dest, send_sizes; - std::vector send_buffer; - { - auto it = node_to_dest.begin(); - while (it != node_to_dest.end()) - { - // Current destination rank - dest.push_back(it->front()); - - // Find iterator to next destination rank and pack send data - auto it1 - = std::find_if(it, node_to_dest.end(), [r0 = dest.back()](auto& idx) - { return idx[0] != r0; }); - send_sizes.push_back(2 * std::distance(it, it1)); - for (auto itx = it; itx != it1; ++itx) - { - send_buffer.push_back(itx->at(1)); - send_buffer.push_back(itx->at(2)); - } - - it = it1; - } - } - - // Prepare send displacements - std::vector send_disp(send_sizes.size() + 1, 0); - std::partial_sum(send_sizes.begin(), send_sizes.end(), - std::next(send_disp.begin())); - - // Discover src ranks. ParMETIS/KaHIP are not scalable (holding an - // array of size equal to the comm size), so no extra harm in using - // non-scalable neighbourhood detection (which might be faster for - // small rank counts). - const std::vector src - = dolfinx::MPI::compute_graph_edges_pcx(comm, dest); - - // Create neighbourhood communicator - MPI_Comm neigh_comm; - MPI_Dist_graph_create_adjacent(comm, src.size(), src.data(), MPI_UNWEIGHTED, - dest.size(), dest.data(), MPI_UNWEIGHTED, - MPI_INFO_NULL, false, &neigh_comm); - - // Determine receives sizes - std::vector recv_sizes(dest.size()); - send_sizes.reserve(1); - recv_sizes.reserve(1); - MPI_Neighbor_alltoall(send_sizes.data(), 1, MPI_INT, recv_sizes.data(), 1, - MPI_INT, neigh_comm); - - // Prepare receive displacements - std::vector recv_disp(recv_sizes.size() + 1, 0); - std::partial_sum(recv_sizes.begin(), recv_sizes.end(), - std::next(recv_disp.begin())); - - // Send/receive data - std::vector recv_buffer(recv_disp.back()); - MPI_Neighbor_alltoallv(send_buffer.data(), send_sizes.data(), - send_disp.data(), MPI_INT64_T, recv_buffer.data(), - recv_sizes.data(), recv_disp.data(), MPI_INT64_T, - neigh_comm); - MPI_Comm_free(&neigh_comm); - - // Prepare (local node index, destination rank) array. Add local data, - // then add the received data, and the make unique. - std::vector> local_node_to_dest; - for (auto d : part) - { - local_node_to_dest.push_back( - {static_cast(local_node_to_dest.size()), static_cast(d)}); - } - for (std::size_t i = 0; i < recv_buffer.size(); i += 2) - { - std::int64_t idx = recv_buffer[i]; - int d = recv_buffer[i + 1]; - assert(idx >= range0 and idx < range1); - std::int32_t idx_local = idx - range0; - local_node_to_dest.push_back({idx_local, d}); - } - - { - std::ranges::sort(local_node_to_dest); - auto [unique_end, range_end] = std::ranges::unique(local_node_to_dest); - local_node_to_dest.erase(unique_end, range_end); - } - // Compute offsets - std::vector offsets(graph.num_nodes() + 1, 0); - { - std::vector num_dests(graph.num_nodes(), 0); - for (auto x : local_node_to_dest) - ++num_dests[x[0]]; - std::partial_sum(num_dests.begin(), num_dests.end(), - std::next(offsets.begin())); - } - - // Fill data array - std::vector data(offsets.back()); - { - std::vector pos = offsets; - for (auto [x0, x1] : local_node_to_dest) - data[pos[x0]++] = x1; - } - - graph::AdjacencyList g(std::move(data), std::move(offsets)); - - // Make sure the owning rank comes first for each node - for (std::int32_t i = 0; i < g.num_nodes(); ++i) - { - auto d = g.links(i); - auto it = std::find(d.begin(), d.end(), part[i]); - assert(it != d.end()); - std::iter_swap(d.begin(), it); - } - - return g; -} -} // namespace -} // namespace dolfinx - namespace dolfinx::refinement { From c4766dcad705d0e12a97b8f56e44c87262ee5f00 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Wed, 12 Mar 2025 09:08:16 +0100 Subject: [PATCH 07/36] Add docstring --- cpp/dolfinx/refinement/refine.h | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index 39cd82045a1..077d5bd5a1c 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -25,6 +25,17 @@ namespace dolfinx::refinement { +/** + * @brief Create a cell partitioner which maintains the partition of a coarse + * mesh. + * + * @tparam T floating point type of mesh geometry. + * @param parent_mesh mesh indicating the partition scheme to use, i.e. the + * coarse mesh. + * @param parent_cell list of cell indices mapping cells of the new refined mesh + * into the coarse mesh, usually output of `refinement::refine`. + * @return The created cell partition function. + */ template mesh::CellPartitionFunction create_identity_partitioner(const mesh::Mesh& parent_mesh, @@ -71,10 +82,6 @@ create_identity_partitioner(const mesh::Mesh& parent_mesh, std::partial_sum(node_disp.begin(), node_disp.end(), node_disp.begin()); return compute_destination_ranks(comm, dual_graph, node_disp, destinations); - - // std::iota(dest_offsets.begin(), dest_offsets.end(), 0); - // return graph::AdjacencyList(std::move(destinations), - // std::move(dest_offsets)); }; } From bb79adacd9a5f923e578546792672691de12b0e8 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Wed, 12 Mar 2025 09:23:40 +0100 Subject: [PATCH 08/36] Improve sequential code path --- cpp/dolfinx/refinement/refine.h | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index 077d5bd5a1c..eae18da556f 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -18,6 +18,7 @@ #include "plaza.h" #include #include +#include #include #include #include @@ -49,22 +50,20 @@ create_identity_partitioner(const mesh::Mesh& parent_mesh, std::vector> cells) -> graph::AdjacencyList { - auto parent_top = parent_mesh.topology(); - auto parent_cell_im = parent_top->index_map(parent_top->dim()); + auto cell_im + = parent_mesh.topology()->index_map(parent_mesh.topology()->dim()); - int num_cell_vertices = mesh::num_cell_vertices(cell_types.front()); - std::int32_t num_cells = cells.front().size() / num_cell_vertices; + std::int32_t num_cells = cells.front().size() / mesh::num_cell_vertices(cell_types.front()); std::vector destinations(num_cells); int rank = dolfinx::MPI::rank(comm); for (std::int32_t i = 0; i < destinations.size(); i++) { - bool ghost_parent_cell = parent_cell[i] > parent_cell_im->size_local(); + bool ghost_parent_cell = parent_cell[i] > cell_im->size_local(); if (ghost_parent_cell) { destinations[i] - = parent_cell_im - ->owners()[parent_cell[i] - parent_cell_im->size_local()]; + = cell_im->owners()[parent_cell[i] - cell_im->size_local()]; } else { @@ -72,15 +71,18 @@ create_identity_partitioner(const mesh::Mesh& parent_mesh, } } + if (comm == MPI_COMM_NULL) + { + return graph::regular_adjacency_list(std::move(destinations), 1); + } + auto dual_graph = mesh::build_dual_graph(comm, cell_types, cells); - std::vector node_disp; - node_disp = std::vector(MPI::size(comm) + 1, 0); + std::vector node_disp(MPI::size(comm) + 1, 0); std::int32_t local_size = dual_graph.num_nodes(); MPI_Allgather(&local_size, 1, dolfinx::MPI::mpi_t, node_disp.data() + 1, 1, dolfinx::MPI::mpi_t, comm); std::partial_sum(node_disp.begin(), node_disp.end(), node_disp.begin()); - return compute_destination_ranks(comm, dual_graph, node_disp, destinations); }; } From d2ac3336561e80dffcd3f29be4e0800ff8b9dc00 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Wed, 12 Mar 2025 09:28:59 +0100 Subject: [PATCH 09/36] Add explicit instantiations --- cpp/dolfinx/graph/partitioners.cpp | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/cpp/dolfinx/graph/partitioners.cpp b/cpp/dolfinx/graph/partitioners.cpp index e86b6ad6f2b..a013068afaf 100644 --- a/cpp/dolfinx/graph/partitioners.cpp +++ b/cpp/dolfinx/graph/partitioners.cpp @@ -201,6 +201,15 @@ graph::AdjacencyList dolfinx::graph::compute_destination_ranks( return g; } +template graph::AdjacencyList dolfinx::graph::compute_destination_ranks( + MPI_Comm comm, const graph::AdjacencyList& graph, + const std::vector& node_disp, const std::vector& part); + +template graph::AdjacencyList dolfinx::graph::compute_destination_ranks( + MPI_Comm comm, const graph::AdjacencyList& graph, + const std::vector& node_disp, + const std::vector& part); + //----------------------------------------------------------------------------- #ifdef HAS_PARMETIS template From b59a1675e253b86ea08e8a40430c5221c9606eb8 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Wed, 12 Mar 2025 10:01:27 +0100 Subject: [PATCH 10/36] Doxygen fix explicit instantiation --- cpp/dolfinx/graph/partitioners.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/cpp/dolfinx/graph/partitioners.cpp b/cpp/dolfinx/graph/partitioners.cpp index a013068afaf..6d568fc5102 100644 --- a/cpp/dolfinx/graph/partitioners.cpp +++ b/cpp/dolfinx/graph/partitioners.cpp @@ -201,6 +201,7 @@ graph::AdjacencyList dolfinx::graph::compute_destination_ranks( return g; } +/// @cond template graph::AdjacencyList dolfinx::graph::compute_destination_ranks( MPI_Comm comm, const graph::AdjacencyList& graph, const std::vector& node_disp, const std::vector& part); @@ -209,6 +210,7 @@ template graph::AdjacencyList dolfinx::graph::compute_destination_ranks( MPI_Comm comm, const graph::AdjacencyList& graph, const std::vector& node_disp, const std::vector& part); +/// @endcond //----------------------------------------------------------------------------- #ifdef HAS_PARMETIS From 5419dee23c89e4db3cfc1a69ed8c3c5043c54659 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Wed, 12 Mar 2025 10:02:09 +0100 Subject: [PATCH 11/36] Move docstring to header --- cpp/dolfinx/graph/partitioners.h | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/cpp/dolfinx/graph/partitioners.h b/cpp/dolfinx/graph/partitioners.h index 29aa62af97f..26a78309790 100644 --- a/cpp/dolfinx/graph/partitioners.h +++ b/cpp/dolfinx/graph/partitioners.h @@ -12,6 +12,17 @@ namespace dolfinx::graph { +/// @todo Is it un-documented that the owning rank must come first in +/// reach list of edges? +/// +/// @param[in] comm The communicator +/// @param[in] graph Graph, using global indices for graph edges +/// @param[in] node_disp The distribution of graph nodes across MPI +/// ranks. The global index `gidx` of local index `lidx` is `lidx + +/// node_disp[my_rank]`. +/// @param[in] part The destination rank for owned nodes, i.e. `dest[i]` +/// is the destination of the node with local index `i`. +/// @return Destination ranks for each local node. template graph::AdjacencyList compute_destination_ranks( MPI_Comm comm, const graph::AdjacencyList& graph, From 6f1fcbcfb9694af5f07aaee6c16cb846f3285653 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Wed, 12 Mar 2025 10:07:48 +0100 Subject: [PATCH 12/36] Remove cpp docstring --- cpp/dolfinx/graph/partitioners.cpp | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/cpp/dolfinx/graph/partitioners.cpp b/cpp/dolfinx/graph/partitioners.cpp index 6d568fc5102..425ca4e2910 100644 --- a/cpp/dolfinx/graph/partitioners.cpp +++ b/cpp/dolfinx/graph/partitioners.cpp @@ -35,17 +35,6 @@ extern "C" using namespace dolfinx; -/// @todo Is it un-documented that the owning rank must come first in -/// reach list of edges? -/// -/// @param[in] comm The communicator -/// @param[in] graph Graph, using global indices for graph edges -/// @param[in] node_disp The distribution of graph nodes across MPI -/// ranks. The global index `gidx` of local index `lidx` is `lidx + -/// node_disp[my_rank]`. -/// @param[in] part The destination rank for owned nodes, i.e. `dest[i]` -/// is the destination of the node with local index `i`. -/// @return Destination ranks for each local node. template graph::AdjacencyList dolfinx::graph::compute_destination_ranks( MPI_Comm comm, const graph::AdjacencyList& graph, From 701e39906b90610adf22ab3032a214bc7677acd0 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Wed, 12 Mar 2025 15:40:36 +0100 Subject: [PATCH 13/36] Change defaults and add special case of config --- cpp/dolfinx/refinement/refine.h | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index eae18da556f..1b7f44fa91a 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -21,6 +21,7 @@ #include #include #include +#include #include namespace dolfinx::refinement @@ -53,7 +54,8 @@ create_identity_partitioner(const mesh::Mesh& parent_mesh, auto cell_im = parent_mesh.topology()->index_map(parent_mesh.topology()->dim()); - std::int32_t num_cells = cells.front().size() / mesh::num_cell_vertices(cell_types.front()); + std::int32_t num_cells + = cells.front().size() / mesh::num_cell_vertices(cell_types.front()); std::vector destinations(num_cells); int rank = dolfinx::MPI::rank(comm); @@ -123,7 +125,7 @@ std::tuple, std::optional>, refine(const mesh::Mesh& mesh, std::optional> edges, mesh::CellPartitionFunction partitioner = nullptr, - Option option = Option::none) + Option option = Option::parent_cell) { auto topology = mesh.topology(); assert(topology); @@ -137,6 +139,11 @@ refine(const mesh::Mesh& mesh, if (!partitioner) { + if (!parent_cell) + { + throw std::runtime_error( + "Identity partitioner relies on parent cell computation"); + } assert(parent_cell); partitioner = create_identity_partitioner(mesh, parent_cell.value()); } From 1d19a2128f0012fe3c35172b083d4609d4d9381f Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Wed, 12 Mar 2025 15:51:37 +0100 Subject: [PATCH 14/36] Switch to optional to handle cases correctly --- cpp/dolfinx/refinement/refine.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index 1b7f44fa91a..8fa4e61432b 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -124,7 +124,7 @@ std::tuple, std::optional>, std::optional>> refine(const mesh::Mesh& mesh, std::optional> edges, - mesh::CellPartitionFunction partitioner = nullptr, + std::optional partitioner = std::nullopt, Option option = Option::parent_cell) { auto topology = mesh.topology(); @@ -137,7 +137,7 @@ refine(const mesh::Mesh& mesh, ? interval::compute_refinement_data(mesh, edges, option) : plaza::compute_refinement_data(mesh, edges, option); - if (!partitioner) + if (!partitioner.has_value()) { if (!parent_cell) { @@ -150,7 +150,7 @@ refine(const mesh::Mesh& mesh, mesh::Mesh mesh1 = mesh::create_mesh( mesh.comm(), mesh.comm(), cell_adj.array(), mesh.geometry().cmap(), - mesh.comm(), new_vertex_coords, xshape, partitioner); + mesh.comm(), new_vertex_coords, xshape, *partitioner); // Report the number of refined cells const int D = topology->dim(); From 5629a0e5e82495f612971fe36cc7deb5806e6276 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Wed, 12 Mar 2025 15:54:48 +0100 Subject: [PATCH 15/36] Update docstring --- cpp/dolfinx/refinement/refine.h | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index 8fa4e61432b..8675f16876b 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -103,20 +103,17 @@ create_identity_partitioner(const mesh::Mesh& parent_mesh, /// refined mesh will **not** include ghosts cells (cells connected by /// facet to an owned cell) even if the parent mesh is ghosted. /// -/// @warning Passing `nullptr` for `partitioner`, the refined mesh will -/// **not** have ghosts cells even if the parent mesh is ghosted. The -/// possibility to not re-partition the refined mesh and include ghost -/// cells in the refined mesh will be added in a future release. /// /// @param[in] mesh Input mesh to be refined. /// @param[in] edges Indices of the edges that should be split in the /// refinement. If not provided (`std::nullopt`), uniform refinement is /// performed. -/// @param[in] partitioner Partitioner to be used to distribute the -/// refined mesh. If not callable, refined cells will be on the same +/// @param[in] partitioner (Optional) partitioner to be used to distribute the +/// refined mesh. If not provided (`std::nullopt`) the partition of the coarse +/// mesh will be maintained. If not callable, refined cells will be on the same /// process as the parent cell. /// @param[in] option Control the computation of parent facets, parent -/// cells. If an option is not selected, an empty list is returned. +/// cells. /// @return New mesh, and optional parent cell indices and parent facet /// indices. template From c386c6a94e0c6ca61094cdb5d8388820131f93d6 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Tue, 18 Mar 2025 19:36:09 +0100 Subject: [PATCH 16/36] Fix Python export for optional --- python/dolfinx/wrappers/refinement.cpp | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/python/dolfinx/wrappers/refinement.cpp b/python/dolfinx/wrappers/refinement.cpp index 1e7c73a8e96..3be42f40aae 100644 --- a/python/dolfinx/wrappers/refinement.cpp +++ b/python/dolfinx/wrappers/refinement.cpp @@ -51,11 +51,14 @@ void export_refinement(nb::module_& m) std::span(edges.value().data(), edges.value().size())); } - dolfinx_wrappers::part::impl::CppCellPartitionFunction cpp_partitioner - = partitioner.has_value() - ? dolfinx_wrappers::part::impl::create_cell_partitioner_cpp( - partitioner.value()) - : nullptr; + std::optional + cpp_partitioner(std::nullopt); + if (partitioner.has_value()) + { + cpp_partitioner + = dolfinx_wrappers::part::impl::create_cell_partitioner_cpp( + partitioner.value()); + } auto [mesh1, cell, facet] = dolfinx::refinement::refine( mesh, cpp_edges, cpp_partitioner, option); From 6a79e104053919e1f7f060daefae05985af6b6ee Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Tue, 18 Mar 2025 20:10:37 +0100 Subject: [PATCH 17/36] Update python default value for refinement optio --- python/dolfinx/mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index cc11abd7db5..d50d4ee4ce0 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -519,7 +519,7 @@ def refine( msh: Mesh, edges: typing.Optional[np.ndarray] = None, partitioner: typing.Optional[typing.Callable] = None, - option: RefinementOption = RefinementOption.none, + option: RefinementOption = RefinementOption.parent_cell, ) -> tuple[Mesh, npt.NDArray[np.int32], npt.NDArray[np.int8]]: """Refine a mesh. From dba96cce97734b044d2b05c3108ad50eb9422f5c Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Tue, 18 Mar 2025 20:10:48 +0100 Subject: [PATCH 18/36] Fix interval refinement test --- python/test/unit/refinement/test_interval.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/python/test/unit/refinement/test_interval.py b/python/test/unit/refinement/test_interval.py index 27040fb0b90..9496df8c8fd 100644 --- a/python/test/unit/refinement/test_interval.py +++ b/python/test/unit/refinement/test_interval.py @@ -20,7 +20,8 @@ "ghost_mode_refined", [mesh.GhostMode.none, mesh.GhostMode.shared_vertex, mesh.GhostMode.shared_facet], ) -@pytest.mark.parametrize("option", [mesh.RefinementOption.none, mesh.RefinementOption.parent_cell]) +@pytest.mark.parametrize("option", [mesh.RefinementOption.parent_cell]) +# TODO: reactivate mesh.mesh.RefinementOption.none def test_refine_interval(n, ghost_mode, ghost_mode_refined, option): msh = mesh.create_interval(MPI.COMM_WORLD, n, [0, 1], ghost_mode=ghost_mode) msh_refined, edges, vertices = mesh.refine( From c2917ddcbd1a2f9cbaf753b4e31be35231ced5c2 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Mon, 24 Mar 2025 15:45:44 +0100 Subject: [PATCH 19/36] Partial fix --- python/dolfinx/mesh.py | 3 +- python/dolfinx/wrappers/refinement.cpp | 8 + .../test/unit/refinement/test_refinement.py | 197 +++++++++--------- 3 files changed, 109 insertions(+), 99 deletions(-) diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index d50d4ee4ce0..fa3301cdeb6 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -30,7 +30,7 @@ to_string, to_type, ) -from dolfinx.cpp.refinement import RefinementOption +from dolfinx.cpp.refinement import RefinementOption, empty_partitioner from dolfinx.fem import CoordinateElement as _CoordinateElement from dolfinx.fem import coordinate_element as _coordinate_element from dolfinx.graph import AdjacencyList @@ -56,6 +56,7 @@ "create_unit_cube", "create_unit_interval", "create_unit_square", + "empty_partitioner", "entities_to_geometry", "exterior_facet_indices", "locate_entities", diff --git a/python/dolfinx/wrappers/refinement.cpp b/python/dolfinx/wrappers/refinement.cpp index 3be42f40aae..02c00e3c5b2 100644 --- a/python/dolfinx/wrappers/refinement.cpp +++ b/python/dolfinx/wrappers/refinement.cpp @@ -94,6 +94,14 @@ void refinement(nb::module_& m) export_refinement(m); export_refinement(m); + m.def("empty_partitioner", + [&] + { + return std::make_optional< + dolfinx_wrappers::part::impl::PythonCellPartitionFunction>( + nullptr); + }); + nb::enum_(m, "RefinementOption") .value("none", dolfinx::refinement::Option::none) .value("parent_facet", dolfinx::refinement::Option::parent_facet) diff --git a/python/test/unit/refinement/test_refinement.py b/python/test/unit/refinement/test_refinement.py index 0a8ee67b79a..5db12147581 100644 --- a/python/test/unit/refinement/test_refinement.py +++ b/python/test/unit/refinement/test_refinement.py @@ -21,6 +21,7 @@ compute_incident_entities, create_unit_cube, create_unit_square, + empty_partitioner, locate_entities, locate_entities_boundary, meshtags, @@ -33,7 +34,7 @@ def test_refine_create_unit_square(): """Refine mesh of unit square.""" mesh = create_unit_square(MPI.COMM_WORLD, 5, 7, ghost_mode=GhostMode.none) mesh.topology.create_entities(1) - mesh_refined, _, _ = refine(mesh) + mesh_refined, _, _ = refine(mesh, partitioner=empty_partitioner()) assert mesh_refined.topology.index_map(0).size_global == 165 assert mesh_refined.topology.index_map(2).size_global == 280 @@ -59,7 +60,7 @@ def test_refine_create_form(): """Check that forms can be assembled on refined mesh""" mesh = create_unit_cube(MPI.COMM_WORLD, 3, 3, 3) mesh.topology.create_entities(1) - mesh, _, _ = refine(mesh) + mesh, _, _ = refine(mesh, partitioner=empty_partitioner()) V = functionspace(mesh, ("Lagrange", 1)) # Define variational problem @@ -83,7 +84,7 @@ def left_corner_edge(x, tol=1e-7): if MPI.COMM_WORLD.size == 1: assert edges == 1 - msh1, _, _ = refine(msh, edges) + msh1, _, _ = refine(msh, edges, partitioner=empty_partitioner()) assert msh.topology.index_map(2).size_global + 3 == msh1.topology.index_map(2).size_global @@ -105,110 +106,110 @@ def left_side(x, tol=1e-16): edges = compute_incident_entities(msh.topology, cells, 2, 1) if MPI.COMM_WORLD.size == 1: assert len(edges) == Nx // 2 * (2 * Ny + 1) + (Nx // 2 + 1) * Ny - mesh2, _, _ = refine(msh, edges) + mesh2, _, _ = refine(msh, edges, partitioner=empty_partitioner()) num_cells_global = mesh2.topology.index_map(2).size_global actual_cells = 3 * (Nx * Ny) + 3 * Ny + 2 * Nx * Ny assert num_cells_global == actual_cells -@pytest.mark.parametrize("tdim", [2, 3]) -@pytest.mark.parametrize( - "refine_plaza_wrapper", - [ - lambda msh: refine(msh, partitioner=None, option=RefinementOption.parent_cell_and_facet), - lambda msh: refine( - msh, - edges=np.arange(msh.topology.index_map(1).size_local), - partitioner=None, - option=RefinementOption.parent_cell_and_facet, - ), - ], -) -def test_refine_facet_meshtag(tdim, refine_plaza_wrapper): - if tdim == 3: - msh = create_unit_cube( - MPI.COMM_WORLD, 2, 3, 5, CellType.tetrahedron, ghost_mode=GhostMode.none - ) - else: - msh = create_unit_square(MPI.COMM_WORLD, 2, 5, CellType.triangle, ghost_mode=GhostMode.none) - msh.topology.create_entities(tdim - 1) - msh.topology.create_connectivity(tdim - 1, tdim) - msh.topology.create_entities(1) - f_to_c = msh.topology.connectivity(tdim - 1, tdim) - facet_indices = [] - for f in range(msh.topology.index_map(tdim - 1).size_local): - if len(f_to_c.links(f)) == 1: - facet_indices += [f] - meshtag = meshtags( - msh, - tdim - 1, - np.array(facet_indices, dtype=np.int32), - np.arange(len(facet_indices), dtype=np.int32), - ) - - msh1, parent_cell, parent_facet = refine_plaza_wrapper(msh) - - msh1.topology.create_entities(tdim - 1) - new_meshtag = transfer_meshtag(meshtag, msh1, parent_cell, parent_facet) - assert len(new_meshtag.indices) == (tdim * 2 - 2) * len(meshtag.indices) - - # New tags should be on facets with one cell (i.e. exterior) - msh1.topology.create_connectivity(tdim - 1, tdim) - new_f_to_c = msh1.topology.connectivity(tdim - 1, tdim) - for f in new_meshtag.indices: - assert len(new_f_to_c.links(f)) == 1 - - # Now mark all facets (including internal) - facet_indices = np.arange(msh.topology.index_map(tdim - 1).size_local) - meshtag = meshtags( - msh, - tdim - 1, - np.array(facet_indices, dtype=np.int32), - np.arange(len(facet_indices), dtype=np.int32), - ) - new_meshtag = transfer_meshtag(meshtag, msh1, parent_cell, parent_facet) - assert len(new_meshtag.indices) == (tdim * 2 - 2) * len(meshtag.indices) - - -@pytest.mark.parametrize("tdim", [2, 3]) -@pytest.mark.parametrize( - "refine_plaza_wrapper", - [ - lambda msh: refine(msh, partitioner=None, option=RefinementOption.parent_cell_and_facet), - lambda msh: refine( - msh, - np.arange(msh.topology.index_map(1).size_local), - partitioner=None, - option=RefinementOption.parent_cell_and_facet, - ), - ], -) -def test_refine_cell_meshtag(tdim, refine_plaza_wrapper): - if tdim == 3: - msh = create_unit_cube( - MPI.COMM_WORLD, 2, 3, 5, CellType.tetrahedron, ghost_mode=GhostMode.none - ) - else: - msh = create_unit_square(MPI.COMM_WORLD, 2, 5, CellType.triangle, ghost_mode=GhostMode.none) - - msh.topology.create_entities(1) - cell_indices = np.arange(msh.topology.index_map(tdim).size_local) - meshtag = meshtags( - msh, - tdim, - np.array(cell_indices, dtype=np.int32), - np.arange(len(cell_indices), dtype=np.int32), - ) - - msh1, parent_cell, _ = refine_plaza_wrapper(msh) - new_meshtag = transfer_meshtag(meshtag, msh1, parent_cell) - assert sum(new_meshtag.values) == (tdim * 4 - 4) * sum(meshtag.values) - assert len(new_meshtag.indices) == (tdim * 4 - 4) * len(meshtag.indices) +# @pytest.mark.parametrize("tdim", [2, 3]) +# @pytest.mark.parametrize( +# "refine_plaza_wrapper", +# [ +# lambda msh: refine(msh, partitioner=empty_partitioner(), option=RefinementOption.parent_cell_and_facet), +# lambda msh: refine( +# msh, +# edges=np.arange(msh.topology.index_map(1).size_local), +# partitioner=empty_partitioner(), +# option=RefinementOption.parent_cell_and_facet, +# ), +# ], +# ) +# def test_refine_facet_meshtag(tdim, refine_plaza_wrapper): +# if tdim == 3: +# msh = create_unit_cube( +# MPI.COMM_WORLD, 2, 3, 5, CellType.tetrahedron, ghost_mode=GhostMode.none +# ) +# else: +# msh = create_unit_square(MPI.COMM_WORLD, 2, 5, CellType.triangle, ghost_mode=GhostMode.none) +# msh.topology.create_entities(tdim - 1) +# msh.topology.create_connectivity(tdim - 1, tdim) +# msh.topology.create_entities(1) +# f_to_c = msh.topology.connectivity(tdim - 1, tdim) +# facet_indices = [] +# for f in range(msh.topology.index_map(tdim - 1).size_local): +# if len(f_to_c.links(f)) == 1: +# facet_indices += [f] +# meshtag = meshtags( +# msh, +# tdim - 1, +# np.array(facet_indices, dtype=np.int32), +# np.arange(len(facet_indices), dtype=np.int32), +# ) + +# msh1, parent_cell, parent_facet = refine_plaza_wrapper(msh) + +# msh1.topology.create_entities(tdim - 1) +# new_meshtag = transfer_meshtag(meshtag, msh1, parent_cell, parent_facet) +# assert len(new_meshtag.indices) == (tdim * 2 - 2) * len(meshtag.indices) + +# # New tags should be on facets with one cell (i.e. exterior) +# msh1.topology.create_connectivity(tdim - 1, tdim) +# new_f_to_c = msh1.topology.connectivity(tdim - 1, tdim) +# for f in new_meshtag.indices: +# assert len(new_f_to_c.links(f)) == 1 + +# # Now mark all facets (including internal) +# facet_indices = np.arange(msh.topology.index_map(tdim - 1).size_local) +# meshtag = meshtags( +# msh, +# tdim - 1, +# np.array(facet_indices, dtype=np.int32), +# np.arange(len(facet_indices), dtype=np.int32), +# ) +# new_meshtag = transfer_meshtag(meshtag, msh1, parent_cell, parent_facet) +# assert len(new_meshtag.indices) == (tdim * 2 - 2) * len(meshtag.indices) + + +# @pytest.mark.parametrize("tdim", [2, 3]) +# @pytest.mark.parametrize( +# "refine_plaza_wrapper", +# [ +# lambda msh: refine(msh, partitioner=empty_partitioner(), option=RefinementOption.parent_cell_and_facet), +# lambda msh: refine( +# msh, +# np.arange(msh.topology.index_map(1).size_local), +# partitioner=empty_partitioner(), +# option=RefinementOption.parent_cell_and_facet, +# ), +# ], +# ) +# def test_refine_cell_meshtag(tdim, refine_plaza_wrapper): +# if tdim == 3: +# msh = create_unit_cube( +# MPI.COMM_WORLD, 2, 3, 5, CellType.tetrahedron, ghost_mode=GhostMode.none +# ) +# else: +# msh = create_unit_square(MPI.COMM_WORLD, 2, 5, CellType.triangle, ghost_mode=GhostMode.none) + +# msh.topology.create_entities(1) +# cell_indices = np.arange(msh.topology.index_map(tdim).size_local) +# meshtag = meshtags( +# msh, +# tdim, +# np.array(cell_indices, dtype=np.int32), +# np.arange(len(cell_indices), dtype=np.int32), +# ) + +# msh1, parent_cell, _ = refine_plaza_wrapper(msh) +# new_meshtag = transfer_meshtag(meshtag, msh1, parent_cell) +# assert sum(new_meshtag.values) == (tdim * 4 - 4) * sum(meshtag.values) +# assert len(new_meshtag.indices) == (tdim * 4 - 4) * len(meshtag.indices) def test_refine_ufl_cargo(): msh = create_unit_cube(MPI.COMM_WORLD, 4, 3, 3) msh.topology.create_entities(1) - msh1, _, _ = refine(msh) + msh1, _, _ = refine(msh, partitioner=empty_partitioner()) assert msh1.ufl_domain().ufl_cargo() != msh.ufl_domain().ufl_cargo() From 8973edc68338fa2e9106702187b10b4f9e5e3fa1 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Tue, 25 Mar 2025 11:01:07 +0100 Subject: [PATCH 20/36] Track down further --- .../test/unit/refinement/test_refinement.py | 196 +++++++++--------- 1 file changed, 101 insertions(+), 95 deletions(-) diff --git a/python/test/unit/refinement/test_refinement.py b/python/test/unit/refinement/test_refinement.py index 5db12147581..ddefb10953b 100644 --- a/python/test/unit/refinement/test_refinement.py +++ b/python/test/unit/refinement/test_refinement.py @@ -24,9 +24,9 @@ empty_partitioner, locate_entities, locate_entities_boundary, - meshtags, + # meshtags, refine, - transfer_meshtag, + # transfer_meshtag, ) @@ -113,99 +113,105 @@ def left_side(x, tol=1e-16): assert num_cells_global == actual_cells -# @pytest.mark.parametrize("tdim", [2, 3]) -# @pytest.mark.parametrize( -# "refine_plaza_wrapper", -# [ -# lambda msh: refine(msh, partitioner=empty_partitioner(), option=RefinementOption.parent_cell_and_facet), -# lambda msh: refine( -# msh, -# edges=np.arange(msh.topology.index_map(1).size_local), -# partitioner=empty_partitioner(), -# option=RefinementOption.parent_cell_and_facet, -# ), -# ], -# ) -# def test_refine_facet_meshtag(tdim, refine_plaza_wrapper): -# if tdim == 3: -# msh = create_unit_cube( -# MPI.COMM_WORLD, 2, 3, 5, CellType.tetrahedron, ghost_mode=GhostMode.none -# ) -# else: -# msh = create_unit_square(MPI.COMM_WORLD, 2, 5, CellType.triangle, ghost_mode=GhostMode.none) -# msh.topology.create_entities(tdim - 1) -# msh.topology.create_connectivity(tdim - 1, tdim) -# msh.topology.create_entities(1) -# f_to_c = msh.topology.connectivity(tdim - 1, tdim) -# facet_indices = [] -# for f in range(msh.topology.index_map(tdim - 1).size_local): -# if len(f_to_c.links(f)) == 1: -# facet_indices += [f] -# meshtag = meshtags( -# msh, -# tdim - 1, -# np.array(facet_indices, dtype=np.int32), -# np.arange(len(facet_indices), dtype=np.int32), -# ) - -# msh1, parent_cell, parent_facet = refine_plaza_wrapper(msh) - -# msh1.topology.create_entities(tdim - 1) -# new_meshtag = transfer_meshtag(meshtag, msh1, parent_cell, parent_facet) -# assert len(new_meshtag.indices) == (tdim * 2 - 2) * len(meshtag.indices) - -# # New tags should be on facets with one cell (i.e. exterior) -# msh1.topology.create_connectivity(tdim - 1, tdim) -# new_f_to_c = msh1.topology.connectivity(tdim - 1, tdim) -# for f in new_meshtag.indices: -# assert len(new_f_to_c.links(f)) == 1 - -# # Now mark all facets (including internal) -# facet_indices = np.arange(msh.topology.index_map(tdim - 1).size_local) -# meshtag = meshtags( -# msh, -# tdim - 1, -# np.array(facet_indices, dtype=np.int32), -# np.arange(len(facet_indices), dtype=np.int32), -# ) -# new_meshtag = transfer_meshtag(meshtag, msh1, parent_cell, parent_facet) -# assert len(new_meshtag.indices) == (tdim * 2 - 2) * len(meshtag.indices) - - -# @pytest.mark.parametrize("tdim", [2, 3]) -# @pytest.mark.parametrize( -# "refine_plaza_wrapper", -# [ -# lambda msh: refine(msh, partitioner=empty_partitioner(), option=RefinementOption.parent_cell_and_facet), -# lambda msh: refine( -# msh, -# np.arange(msh.topology.index_map(1).size_local), -# partitioner=empty_partitioner(), -# option=RefinementOption.parent_cell_and_facet, -# ), -# ], -# ) -# def test_refine_cell_meshtag(tdim, refine_plaza_wrapper): -# if tdim == 3: -# msh = create_unit_cube( -# MPI.COMM_WORLD, 2, 3, 5, CellType.tetrahedron, ghost_mode=GhostMode.none -# ) -# else: -# msh = create_unit_square(MPI.COMM_WORLD, 2, 5, CellType.triangle, ghost_mode=GhostMode.none) - -# msh.topology.create_entities(1) -# cell_indices = np.arange(msh.topology.index_map(tdim).size_local) -# meshtag = meshtags( -# msh, -# tdim, -# np.array(cell_indices, dtype=np.int32), -# np.arange(len(cell_indices), dtype=np.int32), -# ) - -# msh1, parent_cell, _ = refine_plaza_wrapper(msh) -# new_meshtag = transfer_meshtag(meshtag, msh1, parent_cell) -# assert sum(new_meshtag.values) == (tdim * 4 - 4) * sum(meshtag.values) -# assert len(new_meshtag.indices) == (tdim * 4 - 4) * len(meshtag.indices) +@pytest.mark.parametrize("tdim", [2, 3]) +@pytest.mark.parametrize( + "refine_plaza_wrapper", + [ + lambda msh: refine( + msh, partitioner=empty_partitioner(), option=RefinementOption.parent_cell_and_facet + ), + lambda msh: refine( + msh, + edges=np.arange(msh.topology.index_map(1).size_local), + partitioner=empty_partitioner(), + option=RefinementOption.parent_cell_and_facet, + ), + ], +) +def test_refine_facet_meshtag(tdim, refine_plaza_wrapper): + if tdim == 3: + msh = create_unit_cube( + MPI.COMM_WORLD, 2, 3, 5, CellType.tetrahedron, ghost_mode=GhostMode.none + ) + else: + msh = create_unit_square(MPI.COMM_WORLD, 2, 5, CellType.triangle, ghost_mode=GhostMode.none) + msh.topology.create_entities(tdim - 1) + msh.topology.create_connectivity(tdim - 1, tdim) + msh.topology.create_entities(1) + f_to_c = msh.topology.connectivity(tdim - 1, tdim) + facet_indices = [] + for f in range(msh.topology.index_map(tdim - 1).size_local): + if len(f_to_c.links(f)) == 1: + facet_indices += [f] + # meshtag = meshtags( + # msh, + # tdim - 1, + # np.array(facet_indices, dtype=np.int32), + # np.arange(len(facet_indices), dtype=np.int32), + # ) + + msh1, parent_cell, parent_facet = refine_plaza_wrapper(msh) + + msh1.topology.create_entities(tdim - 1) + # TODO: crashes below + # new_meshtag = transfer_meshtag(meshtag, msh1, parent_cell, parent_facet) + # assert len(new_meshtag.indices) == (tdim * 2 - 2) * len(meshtag.indices) + + # # New tags should be on facets with one cell (i.e. exterior) + # msh1.topology.create_connectivity(tdim - 1, tdim) + # new_f_to_c = msh1.topology.connectivity(tdim - 1, tdim) + # for f in new_meshtag.indices: + # assert len(new_f_to_c.links(f)) == 1 + + # # Now mark all facets (including internal) + # facet_indices = np.arange(msh.topology.index_map(tdim - 1).size_local) + # meshtag = meshtags( + # msh, + # tdim - 1, + # np.array(facet_indices, dtype=np.int32), + # np.arange(len(facet_indices), dtype=np.int32), + # ) + # new_meshtag = transfer_meshtag(meshtag, msh1, parent_cell, parent_facet) + # assert len(new_meshtag.indices) == (tdim * 2 - 2) * len(meshtag.indices) + + +@pytest.mark.parametrize("tdim", [2, 3]) +@pytest.mark.parametrize( + "refine_plaza_wrapper", + [ + lambda msh: refine( + msh, partitioner=empty_partitioner(), option=RefinementOption.parent_cell_and_facet + ), + lambda msh: refine( + msh, + np.arange(msh.topology.index_map(1).size_local), + partitioner=empty_partitioner(), + option=RefinementOption.parent_cell_and_facet, + ), + ], +) +def test_refine_cell_meshtag(tdim, refine_plaza_wrapper): + if tdim == 3: + msh = create_unit_cube( + MPI.COMM_WORLD, 2, 3, 5, CellType.tetrahedron, ghost_mode=GhostMode.none + ) + else: + msh = create_unit_square(MPI.COMM_WORLD, 2, 5, CellType.triangle, ghost_mode=GhostMode.none) + + msh.topology.create_entities(1) + # cell_indices = np.arange(msh.topology.index_map(tdim).size_local) + # meshtag = meshtags( + # msh, + # tdim, + # np.array(cell_indices, dtype=np.int32), + # np.arange(len(cell_indices), dtype=np.int32), + # ) + + msh1, parent_cell, _ = refine_plaza_wrapper(msh) + # TODO: crashes below + # new_meshtag = transfer_meshtag(meshtag, msh1, parent_cell) + # assert sum(new_meshtag.values) == (tdim * 4 - 4) * sum(meshtag.values) + # assert len(new_meshtag.indices) == (tdim * 4 - 4) * len(meshtag.indices) def test_refine_ufl_cargo(): From b769eeb333298754c28c1d5a09f56b56073087fe Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Tue, 25 Mar 2025 11:34:59 +0100 Subject: [PATCH 21/36] Showcase None export problem --- python/dolfinx/wrappers/refinement.cpp | 12 ++++++++++-- python/test/unit/refinement/test_refinement.py | 4 ++++ 2 files changed, 14 insertions(+), 2 deletions(-) diff --git a/python/dolfinx/wrappers/refinement.cpp b/python/dolfinx/wrappers/refinement.cpp index 02c00e3c5b2..e6c20edb81f 100644 --- a/python/dolfinx/wrappers/refinement.cpp +++ b/python/dolfinx/wrappers/refinement.cpp @@ -25,6 +25,7 @@ #include #include #include +#include namespace nb = nanobind; @@ -95,11 +96,18 @@ void refinement(nb::module_& m) export_refinement(m); m.def("empty_partitioner", - [&] + []() -> std::optional< + dolfinx_wrappers::part::impl::PythonCellPartitionFunction> { - return std::make_optional< + auto empty_partitioner = std::make_optional< dolfinx_wrappers::part::impl::PythonCellPartitionFunction>( nullptr); + assert(empty_partitioner.has_value()); + if (!empty_partitioner.has_value()) + { + throw std::runtime_error("empty_partitioner has no value!"); + } + return empty_partitioner; }); nb::enum_(m, "RefinementOption") diff --git a/python/test/unit/refinement/test_refinement.py b/python/test/unit/refinement/test_refinement.py index ddefb10953b..2e9c3109539 100644 --- a/python/test/unit/refinement/test_refinement.py +++ b/python/test/unit/refinement/test_refinement.py @@ -30,6 +30,10 @@ ) +def test_empty_partitioner(): + assert empty_partitioner() is not None + + def test_refine_create_unit_square(): """Refine mesh of unit square.""" mesh = create_unit_square(MPI.COMM_WORLD, 5, 7, ghost_mode=GhostMode.none) From fe6de20f53d7251f11571a630cc15c8c251eb0ad Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Mon, 31 Mar 2025 16:08:02 +0200 Subject: [PATCH 22/36] Remove empty_partitioner and introduce helper variable as sentinel --- python/dolfinx/mesh.py | 5 +- python/dolfinx/wrappers/refinement.cpp | 35 +++--- .../test/unit/refinement/test_refinement.py | 103 ++++++++---------- 3 files changed, 61 insertions(+), 82 deletions(-) diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index fa3301cdeb6..f2b2b24994b 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -30,7 +30,7 @@ to_string, to_type, ) -from dolfinx.cpp.refinement import RefinementOption, empty_partitioner +from dolfinx.cpp.refinement import RefinementOption from dolfinx.fem import CoordinateElement as _CoordinateElement from dolfinx.fem import coordinate_element as _coordinate_element from dolfinx.graph import AdjacencyList @@ -519,6 +519,7 @@ def transfer_meshtag( def refine( msh: Mesh, edges: typing.Optional[np.ndarray] = None, + redistribute: bool = False, partitioner: typing.Optional[typing.Callable] = None, option: RefinementOption = RefinementOption.parent_cell, ) -> tuple[Mesh, npt.NDArray[np.int32], npt.NDArray[np.int8]]: @@ -553,7 +554,7 @@ def refine( Refined mesh, (optional) parent cells, (optional) parent facets """ mesh1, parent_cell, parent_facet = _cpp.refinement.refine( - msh._cpp_object, edges, partitioner, option + msh._cpp_object, edges, redistribute, partitioner, option ) # Create new ufl domain as it will carry a reference to the C++ mesh # in the ufl_cargo diff --git a/python/dolfinx/wrappers/refinement.cpp b/python/dolfinx/wrappers/refinement.cpp index e6c20edb81f..28d2cbe7954 100644 --- a/python/dolfinx/wrappers/refinement.cpp +++ b/python/dolfinx/wrappers/refinement.cpp @@ -40,11 +40,16 @@ void export_refinement(nb::module_& m) std::optional< nb::ndarray, nb::c_contig>> edges, - std::optional< - dolfinx_wrappers::part::impl::PythonCellPartitionFunction> - partitioner, + bool redistribute, + dolfinx_wrappers::part::impl::PythonCellPartitionFunction partitioner, dolfinx::refinement::Option option) { + if (!redistribute && partitioner) + { + throw std::runtime_error("Providing a partitioner and passing " + "redistribute=False is bugprone."); + } + std::optional> cpp_edges(std::nullopt); if (edges.has_value()) { @@ -54,12 +59,13 @@ void export_refinement(nb::module_& m) std::optional cpp_partitioner(std::nullopt); - if (partitioner.has_value()) + if (redistribute) { cpp_partitioner = dolfinx_wrappers::part::impl::create_cell_partitioner_cpp( - partitioner.value()); + partitioner); } + auto [mesh1, cell, facet] = dolfinx::refinement::refine( mesh, cpp_edges, cpp_partitioner, option); @@ -82,8 +88,8 @@ void export_refinement(nb::module_& m) return std::tuple{std::move(mesh1), std::move(python_cell), std::move(python_facet)}; }, - nb::arg("mesh"), nb::arg("edges").none(), nb::arg("partitioner").none(), - nb::arg("option")); + nb::arg("mesh"), nb::arg("edges").none(), nb::arg("redistribute"), + nb::arg("partitioner").none(), nb::arg("option")); } } // namespace @@ -95,21 +101,6 @@ void refinement(nb::module_& m) export_refinement(m); export_refinement(m); - m.def("empty_partitioner", - []() -> std::optional< - dolfinx_wrappers::part::impl::PythonCellPartitionFunction> - { - auto empty_partitioner = std::make_optional< - dolfinx_wrappers::part::impl::PythonCellPartitionFunction>( - nullptr); - assert(empty_partitioner.has_value()); - if (!empty_partitioner.has_value()) - { - throw std::runtime_error("empty_partitioner has no value!"); - } - return empty_partitioner; - }); - nb::enum_(m, "RefinementOption") .value("none", dolfinx::refinement::Option::none) .value("parent_facet", dolfinx::refinement::Option::parent_facet) diff --git a/python/test/unit/refinement/test_refinement.py b/python/test/unit/refinement/test_refinement.py index 2e9c3109539..06f852014f5 100644 --- a/python/test/unit/refinement/test_refinement.py +++ b/python/test/unit/refinement/test_refinement.py @@ -21,24 +21,19 @@ compute_incident_entities, create_unit_cube, create_unit_square, - empty_partitioner, locate_entities, locate_entities_boundary, - # meshtags, + meshtags, refine, - # transfer_meshtag, + transfer_meshtag, ) -def test_empty_partitioner(): - assert empty_partitioner() is not None - - def test_refine_create_unit_square(): """Refine mesh of unit square.""" mesh = create_unit_square(MPI.COMM_WORLD, 5, 7, ghost_mode=GhostMode.none) mesh.topology.create_entities(1) - mesh_refined, _, _ = refine(mesh, partitioner=empty_partitioner()) + mesh_refined, _, _ = refine(mesh) assert mesh_refined.topology.index_map(0).size_global == 165 assert mesh_refined.topology.index_map(2).size_global == 280 @@ -52,7 +47,7 @@ def test_refine_create_unit_cube(ghost_mode, redistribute): """Refine mesh of unit cube.""" mesh = create_unit_cube(MPI.COMM_WORLD, 5, 7, 9, ghost_mode=ghost_mode) mesh.topology.create_entities(1) - mesh, _, _ = refine(mesh, partitioner=create_cell_partitioner(ghost_mode)) + mesh, _, _ = refine(mesh, partitioner=create_cell_partitioner(ghost_mode), redistribute=True) assert mesh.topology.index_map(0).size_global == 3135 assert mesh.topology.index_map(3).size_global == 15120 @@ -64,7 +59,7 @@ def test_refine_create_form(): """Check that forms can be assembled on refined mesh""" mesh = create_unit_cube(MPI.COMM_WORLD, 3, 3, 3) mesh.topology.create_entities(1) - mesh, _, _ = refine(mesh, partitioner=empty_partitioner()) + mesh, _, _ = refine(mesh) V = functionspace(mesh, ("Lagrange", 1)) # Define variational problem @@ -88,7 +83,7 @@ def left_corner_edge(x, tol=1e-7): if MPI.COMM_WORLD.size == 1: assert edges == 1 - msh1, _, _ = refine(msh, edges, partitioner=empty_partitioner()) + msh1, _, _ = refine(msh, edges) assert msh.topology.index_map(2).size_global + 3 == msh1.topology.index_map(2).size_global @@ -110,7 +105,7 @@ def left_side(x, tol=1e-16): edges = compute_incident_entities(msh.topology, cells, 2, 1) if MPI.COMM_WORLD.size == 1: assert len(edges) == Nx // 2 * (2 * Ny + 1) + (Nx // 2 + 1) * Ny - mesh2, _, _ = refine(msh, edges, partitioner=empty_partitioner()) + mesh2, _, _ = refine(msh, edges) num_cells_global = mesh2.topology.index_map(2).size_global actual_cells = 3 * (Nx * Ny) + 3 * Ny + 2 * Nx * Ny @@ -121,13 +116,10 @@ def left_side(x, tol=1e-16): @pytest.mark.parametrize( "refine_plaza_wrapper", [ - lambda msh: refine( - msh, partitioner=empty_partitioner(), option=RefinementOption.parent_cell_and_facet - ), + lambda msh: refine(msh, option=RefinementOption.parent_cell_and_facet), lambda msh: refine( msh, edges=np.arange(msh.topology.index_map(1).size_local), - partitioner=empty_partitioner(), option=RefinementOption.parent_cell_and_facet, ), ], @@ -147,49 +139,45 @@ def test_refine_facet_meshtag(tdim, refine_plaza_wrapper): for f in range(msh.topology.index_map(tdim - 1).size_local): if len(f_to_c.links(f)) == 1: facet_indices += [f] - # meshtag = meshtags( - # msh, - # tdim - 1, - # np.array(facet_indices, dtype=np.int32), - # np.arange(len(facet_indices), dtype=np.int32), - # ) + meshtag = meshtags( + msh, + tdim - 1, + np.array(facet_indices, dtype=np.int32), + np.arange(len(facet_indices), dtype=np.int32), + ) msh1, parent_cell, parent_facet = refine_plaza_wrapper(msh) msh1.topology.create_entities(tdim - 1) - # TODO: crashes below - # new_meshtag = transfer_meshtag(meshtag, msh1, parent_cell, parent_facet) - # assert len(new_meshtag.indices) == (tdim * 2 - 2) * len(meshtag.indices) - - # # New tags should be on facets with one cell (i.e. exterior) - # msh1.topology.create_connectivity(tdim - 1, tdim) - # new_f_to_c = msh1.topology.connectivity(tdim - 1, tdim) - # for f in new_meshtag.indices: - # assert len(new_f_to_c.links(f)) == 1 - - # # Now mark all facets (including internal) - # facet_indices = np.arange(msh.topology.index_map(tdim - 1).size_local) - # meshtag = meshtags( - # msh, - # tdim - 1, - # np.array(facet_indices, dtype=np.int32), - # np.arange(len(facet_indices), dtype=np.int32), - # ) - # new_meshtag = transfer_meshtag(meshtag, msh1, parent_cell, parent_facet) - # assert len(new_meshtag.indices) == (tdim * 2 - 2) * len(meshtag.indices) + new_meshtag = transfer_meshtag(meshtag, msh1, parent_cell, parent_facet) + assert len(new_meshtag.indices) == (tdim * 2 - 2) * len(meshtag.indices) + + # New tags should be on facets with one cell (i.e. exterior) + msh1.topology.create_connectivity(tdim - 1, tdim) + new_f_to_c = msh1.topology.connectivity(tdim - 1, tdim) + for f in new_meshtag.indices: + assert len(new_f_to_c.links(f)) == 1 + + # Now mark all facets (including internal) + facet_indices = np.arange(msh.topology.index_map(tdim - 1).size_local) + meshtag = meshtags( + msh, + tdim - 1, + np.array(facet_indices, dtype=np.int32), + np.arange(len(facet_indices), dtype=np.int32), + ) + new_meshtag = transfer_meshtag(meshtag, msh1, parent_cell, parent_facet) + assert len(new_meshtag.indices) == (tdim * 2 - 2) * len(meshtag.indices) @pytest.mark.parametrize("tdim", [2, 3]) @pytest.mark.parametrize( "refine_plaza_wrapper", [ - lambda msh: refine( - msh, partitioner=empty_partitioner(), option=RefinementOption.parent_cell_and_facet - ), + lambda msh: refine(msh, option=RefinementOption.parent_cell_and_facet), lambda msh: refine( msh, np.arange(msh.topology.index_map(1).size_local), - partitioner=empty_partitioner(), option=RefinementOption.parent_cell_and_facet, ), ], @@ -203,23 +191,22 @@ def test_refine_cell_meshtag(tdim, refine_plaza_wrapper): msh = create_unit_square(MPI.COMM_WORLD, 2, 5, CellType.triangle, ghost_mode=GhostMode.none) msh.topology.create_entities(1) - # cell_indices = np.arange(msh.topology.index_map(tdim).size_local) - # meshtag = meshtags( - # msh, - # tdim, - # np.array(cell_indices, dtype=np.int32), - # np.arange(len(cell_indices), dtype=np.int32), - # ) + cell_indices = np.arange(msh.topology.index_map(tdim).size_local) + meshtag = meshtags( + msh, + tdim, + np.array(cell_indices, dtype=np.int32), + np.arange(len(cell_indices), dtype=np.int32), + ) msh1, parent_cell, _ = refine_plaza_wrapper(msh) - # TODO: crashes below - # new_meshtag = transfer_meshtag(meshtag, msh1, parent_cell) - # assert sum(new_meshtag.values) == (tdim * 4 - 4) * sum(meshtag.values) - # assert len(new_meshtag.indices) == (tdim * 4 - 4) * len(meshtag.indices) + new_meshtag = transfer_meshtag(meshtag, msh1, parent_cell) + assert sum(new_meshtag.values) == (tdim * 4 - 4) * sum(meshtag.values) + assert len(new_meshtag.indices) == (tdim * 4 - 4) * len(meshtag.indices) def test_refine_ufl_cargo(): msh = create_unit_cube(MPI.COMM_WORLD, 4, 3, 3) msh.topology.create_entities(1) - msh1, _, _ = refine(msh, partitioner=empty_partitioner()) + msh1, _, _ = refine(msh) assert msh1.ufl_domain().ufl_cargo() != msh.ufl_domain().ufl_cargo() From 91d02f0e9ab7767215d2fe40fa1f8a2ec6e0d388 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Mon, 31 Mar 2025 16:09:30 +0200 Subject: [PATCH 23/36] Ruff --- python/dolfinx/mesh.py | 1 - 1 file changed, 1 deletion(-) diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index f2b2b24994b..1272f55d9b9 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -56,7 +56,6 @@ "create_unit_cube", "create_unit_interval", "create_unit_square", - "empty_partitioner", "entities_to_geometry", "exterior_facet_indices", "locate_entities", From 6004efbb532207a7819d11a72b56ced1132580ac Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Mon, 31 Mar 2025 16:41:40 +0200 Subject: [PATCH 24/36] Tidy docstring --- python/dolfinx/mesh.py | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index 1272f55d9b9..b57a2d8feef 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -532,17 +532,14 @@ def refine( mesh will **not** include ghosts cells (cells connected by facet to an owned cells) even if the parent mesh is ghosted. - Warning: - Passing ``None`` for ``partitioner``, the refined mesh will - **not** have ghosts cells even if the parent mesh has ghost - cells. The possibility to not re-partition the refined mesh and - include ghost cells in the refined mesh will be added in a - future release. - Args: msh: Mesh from which to create the refined mesh. edges: Indices of edges to split during refinement. If ``None``, mesh refinement is uniform. + redistribute: Controls whether the refined mesh should be redistributed during creation. If + ``False`` (default) the refined mesh will exchange ghosts cells, but distribution of the + coarse mesh is maintained, otherwise the passed ``partitioner`` is used to redistribute + the refined mesh. partitioner: Partitioner to distribute the refined mesh. If ``None`` no redistribution is performed, i.e. refined cells remain on the same process as the parent cell. From fd71eeb5cc1ed61826f9dc0e9a35bb18190a0473 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Mon, 31 Mar 2025 16:41:54 +0200 Subject: [PATCH 25/36] Return to previous test state --- python/test/unit/refinement/test_refinement.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/python/test/unit/refinement/test_refinement.py b/python/test/unit/refinement/test_refinement.py index 06f852014f5..f5484ab7427 100644 --- a/python/test/unit/refinement/test_refinement.py +++ b/python/test/unit/refinement/test_refinement.py @@ -116,11 +116,12 @@ def left_side(x, tol=1e-16): @pytest.mark.parametrize( "refine_plaza_wrapper", [ - lambda msh: refine(msh, option=RefinementOption.parent_cell_and_facet), + lambda msh: refine(msh, option=RefinementOption.parent_cell_and_facet, redistribute=True), lambda msh: refine( msh, edges=np.arange(msh.topology.index_map(1).size_local), option=RefinementOption.parent_cell_and_facet, + redistribute=True, ), ], ) @@ -174,11 +175,12 @@ def test_refine_facet_meshtag(tdim, refine_plaza_wrapper): @pytest.mark.parametrize( "refine_plaza_wrapper", [ - lambda msh: refine(msh, option=RefinementOption.parent_cell_and_facet), + lambda msh: refine(msh, option=RefinementOption.parent_cell_and_facet, redistribute=True), lambda msh: refine( msh, np.arange(msh.topology.index_map(1).size_local), option=RefinementOption.parent_cell_and_facet, + redistribute=True, ), ], ) From fca16712a9cff7950dc99232a6b202c2cb18be68 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Wed, 2 Apr 2025 19:56:44 +0200 Subject: [PATCH 26/36] Add unit test for identity_partitioner --- .../test/unit/refinement/test_refinement.py | 32 ++++++++++++++++++- 1 file changed, 31 insertions(+), 1 deletion(-) diff --git a/python/test/unit/refinement/test_refinement.py b/python/test/unit/refinement/test_refinement.py index f5484ab7427..2f8fe45f99f 100644 --- a/python/test/unit/refinement/test_refinement.py +++ b/python/test/unit/refinement/test_refinement.py @@ -1,4 +1,4 @@ -# Copyright (C) 2018-2024 Chris N Richardson and Jørgen S. Dokken +# Copyright (C) 2018-2025 Chris N Richardson, Jørgen S. Dokken and Paul T. Kühner # This file is part of DOLFINx (https://www.fenicsproject.org) # @@ -20,6 +20,7 @@ RefinementOption, compute_incident_entities, create_unit_cube, + create_unit_interval, create_unit_square, locate_entities, locate_entities_boundary, @@ -212,3 +213,32 @@ def test_refine_ufl_cargo(): msh.topology.create_entities(1) msh1, _, _ = refine(msh) assert msh1.ufl_domain().ufl_cargo() != msh.ufl_domain().ufl_cargo() + + +@pytest.mark.parametrize("tdim", [1, 2, 3]) +@pytest.mark.parametrize( + "ghost_mode", [GhostMode.none, GhostMode.shared_vertex, GhostMode.shared_facet] +) +def test_identity_partitioner(tdim, ghost_mode): + n = 2 + if tdim == 1: + mesh = create_unit_interval(MPI.COMM_WORLD, n) + elif tdim == 2: + mesh = create_unit_square(MPI.COMM_WORLD, n, n) + else: + mesh = create_unit_cube(MPI.COMM_WORLD, n, n, n) + + cells = mesh.topology.index_map(mesh.topology.dim) + owning_process = np.full(cells.size_local + cells.num_ghosts, MPI.COMM_WORLD.rank) + owning_process[cells.size_local :] = cells.owners + + mesh.topology.create_entities(1) + [mesh_fine, parent_cells, _] = refine(mesh, redistribute=False, partitioner=None) + + cells_fine = mesh_fine.topology.index_map(mesh.topology.dim) + owning_process_fine = np.full( + cells_fine.size_local + cells_fine.num_ghosts, MPI.COMM_WORLD.rank + ) + owning_process_fine[cells_fine.size_local :] = cells_fine.owners + + assert np.all(owning_process_fine[: cells_fine.size_local] == owning_process[parent_cells]) From cd8c56f826baf20343102ded23258588d32a44e4 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Thu, 29 May 2025 22:43:21 +0200 Subject: [PATCH 27/36] Fix sign compare --- cpp/dolfinx/refinement/refine.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index 8675f16876b..66ca5f035b8 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -59,7 +59,7 @@ create_identity_partitioner(const mesh::Mesh& parent_mesh, std::vector destinations(num_cells); int rank = dolfinx::MPI::rank(comm); - for (std::int32_t i = 0; i < destinations.size(); i++) + for (std::size_t i = 0; i < destinations.size(); i++) { bool ghost_parent_cell = parent_cell[i] > cell_im->size_local(); if (ghost_parent_cell) From c754ba897f7755ac1f443ffebc558f0a2dde5748 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Thu, 29 May 2025 22:48:52 +0200 Subject: [PATCH 28/36] Tidy logic and consistent naming --- cpp/dolfinx/refinement/refine.h | 17 ++++++----------- 1 file changed, 6 insertions(+), 11 deletions(-) diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index 66ca5f035b8..444bafc07a0 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -51,8 +51,10 @@ create_identity_partitioner(const mesh::Mesh& parent_mesh, std::vector> cells) -> graph::AdjacencyList { - auto cell_im + const auto parent_cell_im = parent_mesh.topology()->index_map(parent_mesh.topology()->dim()); + std::int32_t parent_num_cells = parent_cell_im->size_local(); + std::span parent_cell_owners = parent_cell_im->owners(); std::int32_t num_cells = cells.front().size() / mesh::num_cell_vertices(cell_types.front()); @@ -61,22 +63,15 @@ create_identity_partitioner(const mesh::Mesh& parent_mesh, int rank = dolfinx::MPI::rank(comm); for (std::size_t i = 0; i < destinations.size(); i++) { - bool ghost_parent_cell = parent_cell[i] > cell_im->size_local(); - if (ghost_parent_cell) - { - destinations[i] - = cell_im->owners()[parent_cell[i] - cell_im->size_local()]; - } + bool parent_is_ghost_cell = parent_cell[i] > parent_num_cells; + if (parent_is_ghost_cell) + destinations[i] = parent_cell_owners[parent_cell[i] - parent_num_cells]; else - { destinations[i] = rank; - } } if (comm == MPI_COMM_NULL) - { return graph::regular_adjacency_list(std::move(destinations), 1); - } auto dual_graph = mesh::build_dual_graph(comm, cell_types, cells); std::vector node_disp(MPI::size(comm) + 1, 0); From 6aa0e1cef9e9b5ca3fd3cac5fac06c27e9b6cb6c Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Fri, 20 Jun 2025 13:51:04 +0200 Subject: [PATCH 29/36] Adapt to docstrng length --- python/dolfinx/mesh.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index b2e9a52c1b2..980034796da 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -546,10 +546,11 @@ def refine( msh: Mesh from which to create the refined mesh. edges: Indices of edges to split during refinement. If ``None``, mesh refinement is uniform. - redistribute: Controls whether the refined mesh should be redistributed during creation. If - ``False`` (default) the refined mesh will exchange ghosts cells, but distribution of the - coarse mesh is maintained, otherwise the passed ``partitioner`` is used to redistribute - the refined mesh. + redistribute: Controls whether the refined mesh should be + redistributed during creation. If ``False`` (default) the + refined mesh will exchange ghosts cells, but distribution of + the coarse mesh is maintained, otherwise the passed + ``partitioner`` is used to redistribute the refined mesh. partitioner: Partitioner to distribute the refined mesh. If ``None`` no redistribution is performed, i.e. refined cells remain on the same process as the parent cell. From 5bd2649afe844314075cd4731b2c8d61f7b77ee6 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Mon, 7 Jul 2025 15:25:17 +0200 Subject: [PATCH 30/36] Change to placehodler + variant --- cpp/dolfinx/refinement/refine.h | 17 +++++++++--- python/dolfinx/wrappers/refinement.cpp | 37 +++++++++++++++----------- 2 files changed, 36 insertions(+), 18 deletions(-) diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index 444bafc07a0..6dd5d00961c 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -23,6 +23,7 @@ #include #include #include +#include namespace dolfinx::refinement { @@ -84,6 +85,11 @@ create_identity_partitioner(const mesh::Mesh& parent_mesh, }; } +/// @brief Placeholder for the creation of an identity partitioner in refine. +struct IdentityPartitionerPlaceholder +{ +}; + /// @brief Refine a mesh with markers. /// /// The refined mesh can be optionally re-partitioned across processes. @@ -116,7 +122,9 @@ std::tuple, std::optional>, std::optional>> refine(const mesh::Mesh& mesh, std::optional> edges, - std::optional partitioner = std::nullopt, + std::variant + partitioner + = IdentityPartitionerPlaceholder(), Option option = Option::parent_cell) { auto topology = mesh.topology(); @@ -129,7 +137,7 @@ refine(const mesh::Mesh& mesh, ? interval::compute_refinement_data(mesh, edges, option) : plaza::compute_refinement_data(mesh, edges, option); - if (!partitioner.has_value()) + if (std::holds_alternative(partitioner)) { if (!parent_cell) { @@ -140,9 +148,12 @@ refine(const mesh::Mesh& mesh, partitioner = create_identity_partitioner(mesh, parent_cell.value()); } + assert(std::holds_alternative(partitioner)); + mesh::Mesh mesh1 = mesh::create_mesh( mesh.comm(), mesh.comm(), cell_adj.array(), mesh.geometry().cmap(), - mesh.comm(), new_vertex_coords, xshape, *partitioner); + mesh.comm(), new_vertex_coords, xshape, + std::get(partitioner)); // Report the number of refined cells const int D = topology->dim(); diff --git a/python/dolfinx/wrappers/refinement.cpp b/python/dolfinx/wrappers/refinement.cpp index 9e256d5a120..f24d38fdcf3 100644 --- a/python/dolfinx/wrappers/refinement.cpp +++ b/python/dolfinx/wrappers/refinement.cpp @@ -1,4 +1,5 @@ -// Copyright (C) 2018-2024 Chris N. Richardson and Garth N. Wells +// Copyright (C) 2018-2024 Chris N. Richardson, Garth N. Wells and Paul T. +// Kühner // // This file is part of DOLFINx (https://www.fenicsproject.org) // @@ -23,10 +24,12 @@ #include #include #include +#include #include #include #include #include +#include namespace nb = nanobind; @@ -40,22 +43,21 @@ void export_refinement(nb::module_& m) { return dolfinx::refinement::uniform_refine(mesh); }, nb::arg("mesh")); + nb::class_( + m, "IdentityPartitionerPlaceholder") + .def(nb::init<>()); + m.def( "refine", [](const dolfinx::mesh::Mesh& mesh, std::optional< nb::ndarray, nb::c_contig>> edges, - bool redistribute, - dolfinx_wrappers::part::impl::PythonCellPartitionFunction partitioner, + std::variant + partitioner, dolfinx::refinement::Option option) { - if (!redistribute && partitioner) - { - throw std::runtime_error("Providing a partitioner and passing " - "redistribute=False is bugprone."); - } - std::optional> cpp_edges(std::nullopt); if (edges.has_value()) { @@ -63,13 +65,18 @@ void export_refinement(nb::module_& m) std::span(edges.value().data(), edges.value().size())); } - std::optional - cpp_partitioner(std::nullopt); - if (redistribute) + std::variant + cpp_partitioner + = dolfinx::refinement::IdentityPartitionerPlaceholder(); + if (std::holds_alternative< + dolfinx_wrappers::part::impl::PythonCellPartitionFunction>( + partitioner)) { cpp_partitioner = dolfinx_wrappers::part::impl::create_cell_partitioner_cpp( - partitioner); + std::get(partitioner)); } auto [mesh1, cell, facet] = dolfinx::refinement::refine( @@ -94,8 +101,8 @@ void export_refinement(nb::module_& m) return std::tuple{std::move(mesh1), std::move(python_cell), std::move(python_facet)}; }, - nb::arg("mesh"), nb::arg("edges").none(), nb::arg("redistribute"), - nb::arg("partitioner").none(), nb::arg("option")); + nb::arg("mesh"), nb::arg("edges").none(), nb::arg("partitioner"), + nb::arg("option")); } } // namespace From 0bc2f9b30f89dbe1c173702a5c497f52ec95ab2f Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Mon, 7 Jul 2025 15:42:21 +0200 Subject: [PATCH 31/36] Adapt interface --- python/dolfinx/mesh.py | 9 +++++---- python/test/unit/refinement/test_refinement.py | 13 +++++-------- 2 files changed, 10 insertions(+), 12 deletions(-) diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index 261e5119210..7edb07da0cf 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -31,7 +31,7 @@ to_string, to_type, ) -from dolfinx.cpp.refinement import RefinementOption +from dolfinx.cpp.refinement import IdentityPartitionerPlaceholder, RefinementOption from dolfinx.fem import CoordinateElement as _CoordinateElement from dolfinx.fem import coordinate_element as _coordinate_element from dolfinx.graph import AdjacencyList @@ -550,8 +550,9 @@ def transfer_meshtag( def refine( msh: Mesh, edges: typing.Optional[np.ndarray] = None, - redistribute: bool = False, - partitioner: typing.Optional[typing.Callable] = None, + partitioner: typing.Union[ + typing.Callable, IdentityPartitionerPlaceholder + ] = IdentityPartitionerPlaceholder(), option: RefinementOption = RefinementOption.parent_cell, ) -> tuple[Mesh, npt.NDArray[np.int32], npt.NDArray[np.int8]]: """Refine a mesh. @@ -583,7 +584,7 @@ def refine( Refined mesh, (optional) parent cells, (optional) parent facets """ mesh1, parent_cell, parent_facet = _cpp.refinement.refine( - msh._cpp_object, edges, redistribute, partitioner, option + msh._cpp_object, edges, partitioner, option ) # Create new ufl domain as it will carry a reference to the C++ mesh # in the ufl_cargo diff --git a/python/test/unit/refinement/test_refinement.py b/python/test/unit/refinement/test_refinement.py index 2f8fe45f99f..cdecb556947 100644 --- a/python/test/unit/refinement/test_refinement.py +++ b/python/test/unit/refinement/test_refinement.py @@ -43,12 +43,11 @@ def test_refine_create_unit_square(): @pytest.mark.parametrize("ghost_mode", [GhostMode.none, GhostMode.shared_facet]) -@pytest.mark.parametrize("redistribute", [True, False]) -def test_refine_create_unit_cube(ghost_mode, redistribute): +def test_refine_create_unit_cube(ghost_mode): """Refine mesh of unit cube.""" mesh = create_unit_cube(MPI.COMM_WORLD, 5, 7, 9, ghost_mode=ghost_mode) mesh.topology.create_entities(1) - mesh, _, _ = refine(mesh, partitioner=create_cell_partitioner(ghost_mode), redistribute=True) + mesh, _, _ = refine(mesh, partitioner=create_cell_partitioner(ghost_mode)) assert mesh.topology.index_map(0).size_global == 3135 assert mesh.topology.index_map(3).size_global == 15120 @@ -117,12 +116,11 @@ def left_side(x, tol=1e-16): @pytest.mark.parametrize( "refine_plaza_wrapper", [ - lambda msh: refine(msh, option=RefinementOption.parent_cell_and_facet, redistribute=True), + lambda msh: refine(msh, option=RefinementOption.parent_cell_and_facet), lambda msh: refine( msh, edges=np.arange(msh.topology.index_map(1).size_local), option=RefinementOption.parent_cell_and_facet, - redistribute=True, ), ], ) @@ -176,12 +174,11 @@ def test_refine_facet_meshtag(tdim, refine_plaza_wrapper): @pytest.mark.parametrize( "refine_plaza_wrapper", [ - lambda msh: refine(msh, option=RefinementOption.parent_cell_and_facet, redistribute=True), + lambda msh: refine(msh, option=RefinementOption.parent_cell_and_facet), lambda msh: refine( msh, np.arange(msh.topology.index_map(1).size_local), option=RefinementOption.parent_cell_and_facet, - redistribute=True, ), ], ) @@ -233,7 +230,7 @@ def test_identity_partitioner(tdim, ghost_mode): owning_process[cells.size_local :] = cells.owners mesh.topology.create_entities(1) - [mesh_fine, parent_cells, _] = refine(mesh, redistribute=False, partitioner=None) + [mesh_fine, parent_cells, _] = refine(mesh) cells_fine = mesh_fine.topology.index_map(mesh.topology.dim) owning_process_fine = np.full( From 1f3fd0a8ed6030861612fc32182e403f60d00b42 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Mon, 7 Jul 2025 15:46:07 +0200 Subject: [PATCH 32/36] Update docstring --- python/dolfinx/mesh.py | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index 7edb07da0cf..a48ca5276f3 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -569,14 +569,13 @@ def refine( msh: Mesh from which to create the refined mesh. edges: Indices of edges to split during refinement. If ``None``, mesh refinement is uniform. - redistribute: Controls whether the refined mesh should be - redistributed during creation. If ``False`` (default) the - refined mesh will exchange ghosts cells, but distribution of - the coarse mesh is maintained, otherwise the passed - ``partitioner`` is used to redistribute the refined mesh. - partitioner: Partitioner to distribute the refined mesh. If - ``None`` no redistribution is performed, i.e. refined cells - remain on the same process as the parent cell. + partitioner: Partitioner to distribute the refined mesh. If a + ``IdentityPartitionerPlaceholder`` is passed (default) no + redistribution is performed, i.e. refined cells remain on the + same process as the parent cell, but the ghost layer is + updated. If a custom partitioner is passed, it will be used for + distributing the refined mesh. If ``None`` is passed no + redistribution will happen. option: Controls whether parent cells and/or parent facets are computed. From a3346dbfae5fc4c4712c8ff2aee6e636b5d0f22e Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Mon, 7 Jul 2025 19:50:39 +0200 Subject: [PATCH 33/36] Only once --- python/dolfinx/wrappers/refinement.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/python/dolfinx/wrappers/refinement.cpp b/python/dolfinx/wrappers/refinement.cpp index f24d38fdcf3..2bb2ac63101 100644 --- a/python/dolfinx/wrappers/refinement.cpp +++ b/python/dolfinx/wrappers/refinement.cpp @@ -43,10 +43,6 @@ void export_refinement(nb::module_& m) { return dolfinx::refinement::uniform_refine(mesh); }, nb::arg("mesh")); - nb::class_( - m, "IdentityPartitionerPlaceholder") - .def(nb::init<>()); - m.def( "refine", [](const dolfinx::mesh::Mesh& mesh, @@ -114,6 +110,10 @@ void refinement(nb::module_& m) export_refinement(m); export_refinement(m); + nb::class_( + m, "IdentityPartitionerPlaceholder") + .def(nb::init<>()); + nb::enum_(m, "RefinementOption") .value("none", dolfinx::refinement::Option::none) .value("parent_facet", dolfinx::refinement::Option::parent_facet) From d5f3a25de3a498cd80aebd0340bc06be802cd7da Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Tue, 8 Jul 2025 09:11:24 +0200 Subject: [PATCH 34/36] Remvoe not necessary optional and fix tests --- python/dolfinx/wrappers/refinement.cpp | 26 +++++++++++++------ .../test/unit/refinement/test_refinement.py | 6 +++-- 2 files changed, 22 insertions(+), 10 deletions(-) diff --git a/python/dolfinx/wrappers/refinement.cpp b/python/dolfinx/wrappers/refinement.cpp index 2bb2ac63101..bddc79a24c9 100644 --- a/python/dolfinx/wrappers/refinement.cpp +++ b/python/dolfinx/wrappers/refinement.cpp @@ -50,7 +50,8 @@ void export_refinement(nb::module_& m) nb::ndarray, nb::c_contig>> edges, std::variant + std::optional> partitioner, dolfinx::refinement::Option option) { @@ -65,14 +66,23 @@ void export_refinement(nb::module_& m) dolfinx_wrappers::part::impl::CppCellPartitionFunction> cpp_partitioner = dolfinx::refinement::IdentityPartitionerPlaceholder(); - if (std::holds_alternative< - dolfinx_wrappers::part::impl::PythonCellPartitionFunction>( + if (std::holds_alternative>( partitioner)) { - cpp_partitioner - = dolfinx_wrappers::part::impl::create_cell_partitioner_cpp( - std::get(partitioner)); + auto optional = std::get>( + partitioner); + if (!optional.has_value()) + cpp_partitioner + = dolfinx_wrappers::part::impl::CppCellPartitionFunction( + nullptr); + else + { + cpp_partitioner + = dolfinx_wrappers::part::impl::create_cell_partitioner_cpp( + optional.value()); + } } auto [mesh1, cell, facet] = dolfinx::refinement::refine( @@ -97,7 +107,7 @@ void export_refinement(nb::module_& m) return std::tuple{std::move(mesh1), std::move(python_cell), std::move(python_facet)}; }, - nb::arg("mesh"), nb::arg("edges").none(), nb::arg("partitioner"), + nb::arg("mesh"), nb::arg("edges").none(), nb::arg("partitioner").none(), nb::arg("option")); } } // namespace diff --git a/python/test/unit/refinement/test_refinement.py b/python/test/unit/refinement/test_refinement.py index cdecb556947..f23ed562d23 100644 --- a/python/test/unit/refinement/test_refinement.py +++ b/python/test/unit/refinement/test_refinement.py @@ -116,10 +116,11 @@ def left_side(x, tol=1e-16): @pytest.mark.parametrize( "refine_plaza_wrapper", [ - lambda msh: refine(msh, option=RefinementOption.parent_cell_and_facet), + lambda msh: refine(msh, partitioner=None, option=RefinementOption.parent_cell_and_facet), lambda msh: refine( msh, edges=np.arange(msh.topology.index_map(1).size_local), + partitioner=None, option=RefinementOption.parent_cell_and_facet, ), ], @@ -174,10 +175,11 @@ def test_refine_facet_meshtag(tdim, refine_plaza_wrapper): @pytest.mark.parametrize( "refine_plaza_wrapper", [ - lambda msh: refine(msh, option=RefinementOption.parent_cell_and_facet), + lambda msh: refine(msh, partitioner=None, option=RefinementOption.parent_cell_and_facet), lambda msh: refine( msh, np.arange(msh.topology.index_map(1).size_local), + partitioner=None, option=RefinementOption.parent_cell_and_facet, ), ], From 79f25c3eac6722f6ef4c3d586199fe7b21dced78 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Tue, 8 Jul 2025 09:15:25 +0200 Subject: [PATCH 35/36] Reactivate interval refinement test --- python/test/unit/refinement/test_interval.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/python/test/unit/refinement/test_interval.py b/python/test/unit/refinement/test_interval.py index 9496df8c8fd..c8267af60e8 100644 --- a/python/test/unit/refinement/test_interval.py +++ b/python/test/unit/refinement/test_interval.py @@ -20,13 +20,11 @@ "ghost_mode_refined", [mesh.GhostMode.none, mesh.GhostMode.shared_vertex, mesh.GhostMode.shared_facet], ) -@pytest.mark.parametrize("option", [mesh.RefinementOption.parent_cell]) -# TODO: reactivate mesh.mesh.RefinementOption.none +@pytest.mark.parametrize("option", [mesh.RefinementOption.none, mesh.RefinementOption.parent_cell]) def test_refine_interval(n, ghost_mode, ghost_mode_refined, option): msh = mesh.create_interval(MPI.COMM_WORLD, n, [0, 1], ghost_mode=ghost_mode) msh_refined, edges, vertices = mesh.refine( - msh, - option=option, + msh, option=option, partitioner=mesh.create_cell_partitioner(ghost_mode_refined) ) # TODO: add create_cell_partitioner(ghost_mode) when works From 001524b3df6a37e58fbcc4826727fe144dad62d0 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Wed, 16 Jul 2025 18:40:11 +0200 Subject: [PATCH 36/36] Remove ghostmode shared_vertex --- python/test/unit/refinement/test_refinement.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/python/test/unit/refinement/test_refinement.py b/python/test/unit/refinement/test_refinement.py index f23ed562d23..9b35e9870eb 100644 --- a/python/test/unit/refinement/test_refinement.py +++ b/python/test/unit/refinement/test_refinement.py @@ -215,9 +215,7 @@ def test_refine_ufl_cargo(): @pytest.mark.parametrize("tdim", [1, 2, 3]) -@pytest.mark.parametrize( - "ghost_mode", [GhostMode.none, GhostMode.shared_vertex, GhostMode.shared_facet] -) +@pytest.mark.parametrize("ghost_mode", [GhostMode.none, GhostMode.shared_facet]) def test_identity_partitioner(tdim, ghost_mode): n = 2 if tdim == 1: