From aa7bddda8ccd0423693e742a5b7a27fd59d474ef Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Sun, 1 Mar 2026 09:50:10 +0000 Subject: [PATCH 01/23] Make partitioner options explicit --- cpp/demo/biharmonic/main.cpp | 2 +- cpp/demo/codim_0_assembly/main.cpp | 2 +- cpp/demo/hyperelasticity/main.cpp | 2 +- cpp/demo/interpolation-io/main.cpp | 4 ++-- cpp/demo/poisson/main.cpp | 2 +- cpp/demo/poisson_matrix_free/main.cpp | 2 +- cpp/dolfinx/mesh/utils.h | 3 +-- cpp/test/fem/form.cpp | 2 +- cpp/test/mesh/distributed_mesh.cpp | 2 +- cpp/test/mesh/refinement/interval.cpp | 2 +- cpp/test/mesh/refinement/rectangle.cpp | 3 ++- 11 files changed, 13 insertions(+), 13 deletions(-) diff --git a/cpp/demo/biharmonic/main.cpp b/cpp/demo/biharmonic/main.cpp index 77277910836..88644089380 100644 --- a/cpp/demo/biharmonic/main.cpp +++ b/cpp/demo/biharmonic/main.cpp @@ -157,7 +157,7 @@ int main(int argc, char* argv[]) PetscInitialize(&argc, &argv, nullptr, nullptr); { // Create mesh - auto part = mesh::create_cell_partitioner(mesh::GhostMode::shared_facet); + auto part = mesh::create_cell_partitioner(mesh::GhostMode::shared_facet, 2); auto mesh = std::make_shared>( mesh::create_rectangle(MPI_COMM_WORLD, {{{0.0, 0.0}, {1.0, 1.0}}}, {32, 32}, mesh::CellType::triangle, part)); diff --git a/cpp/demo/codim_0_assembly/main.cpp b/cpp/demo/codim_0_assembly/main.cpp index 5be7086e5e2..caeb81afbf0 100644 --- a/cpp/demo/codim_0_assembly/main.cpp +++ b/cpp/demo/codim_0_assembly/main.cpp @@ -35,7 +35,7 @@ int main(int argc, char* argv[]) auto mesh = std::make_shared>(mesh::create_rectangle( MPI_COMM_WORLD, {{{0.0, 0.0}, {2.0, 1.0}}}, {1, 4}, mesh::CellType::quadrilateral, - mesh::create_cell_partitioner(mesh::GhostMode::shared_facet))); + mesh::create_cell_partitioner(mesh::GhostMode::shared_facet, 2))); basix::FiniteElement element = basix::create_element( basix::element::family::P, basix::cell::type::quadrilateral, 1, diff --git a/cpp/demo/hyperelasticity/main.cpp b/cpp/demo/hyperelasticity/main.cpp index df95ae8b51b..2978e7dc58c 100644 --- a/cpp/demo/hyperelasticity/main.cpp +++ b/cpp/demo/hyperelasticity/main.cpp @@ -159,7 +159,7 @@ int main(int argc, char* argv[]) auto mesh = std::make_shared>(mesh::create_box( MPI_COMM_WORLD, {{{0.0, 0.0, 0.0}, {1.0, 1.0, 1.0}}}, {10, 10, 10}, mesh::CellType::tetrahedron, - mesh::create_cell_partitioner(mesh::GhostMode::none))); + mesh::create_cell_partitioner(mesh::GhostMode::none, 2))); auto element = basix::create_element( basix::element::family::P, basix::cell::type::tetrahedron, 1, diff --git a/cpp/demo/interpolation-io/main.cpp b/cpp/demo/interpolation-io/main.cpp index 2fe2625e0c4..50666a56c7e 100644 --- a/cpp/demo/interpolation-io/main.cpp +++ b/cpp/demo/interpolation-io/main.cpp @@ -209,7 +209,7 @@ int main(int argc, char* argv[]) = std::make_shared>(mesh::create_rectangle( MPI_COMM_WORLD, {{{0.0, 0.0}, {1.0, 1.0}}}, {32, 4}, mesh::CellType::triangle, - mesh::create_cell_partitioner(mesh::GhostMode::none))); + mesh::create_cell_partitioner(mesh::GhostMode::none, 2))); // Create mesh using same topology as mesh0, but with different // scalar type for geometry @@ -217,7 +217,7 @@ int main(int argc, char* argv[]) = std::make_shared>(mesh::create_rectangle( MPI_COMM_WORLD, {{{0.0, 0.0}, {1.0, 1.0}}}, {32, 4}, mesh::CellType::triangle, - mesh::create_cell_partitioner(mesh::GhostMode::none))); + mesh::create_cell_partitioner(mesh::GhostMode::none, 2))); // Interpolate a function in a scalar Lagrange space and output the // result to file for visualisation using different types diff --git a/cpp/demo/poisson/main.cpp b/cpp/demo/poisson/main.cpp index 2684403557f..e33a3aba712 100644 --- a/cpp/demo/poisson/main.cpp +++ b/cpp/demo/poisson/main.cpp @@ -111,7 +111,7 @@ int main(int argc, char* argv[]) { // Create mesh and function space - auto part = mesh::create_cell_partitioner(mesh::GhostMode::shared_facet); + auto part = mesh::create_cell_partitioner(mesh::GhostMode::shared_facet, 2); auto mesh = std::make_shared>( mesh::create_rectangle(MPI_COMM_WORLD, {{{0.0, 0.0}, {2.0, 1.0}}}, {32, 16}, mesh::CellType::triangle, part)); diff --git a/cpp/demo/poisson_matrix_free/main.cpp b/cpp/demo/poisson_matrix_free/main.cpp index 7e199f32c00..2181d956bcd 100644 --- a/cpp/demo/poisson_matrix_free/main.cpp +++ b/cpp/demo/poisson_matrix_free/main.cpp @@ -138,7 +138,7 @@ void solver(MPI_Comm comm) // Create mesh and function space auto mesh = std::make_shared>(mesh::create_rectangle( comm, {{{0.0, 0.0}, {1.0, 1.0}}}, {10, 10}, mesh::CellType::triangle, - mesh::create_cell_partitioner(mesh::GhostMode::none))); + mesh::create_cell_partitioner(mesh::GhostMode::none, 2))); auto element = basix::create_element( basix::element::family::P, basix::cell::type::triangle, 2, basix::element::lagrange_variant::unset, diff --git a/cpp/dolfinx/mesh/utils.h b/cpp/dolfinx/mesh/utils.h index b67c40aa756..513a6c69531 100644 --- a/cpp/dolfinx/mesh/utils.h +++ b/cpp/dolfinx/mesh/utils.h @@ -982,8 +982,7 @@ create_cell_partitioner(mesh::GhostMode ghost_mode, /// @return Function that computes the destination ranks for each cell. CellPartitionFunction create_cell_partitioner(mesh::GhostMode ghost_mode, - std::optional max_facet_to_cell_links - = 2); + std::optional max_facet_to_cell_links); /// @brief Compute incident entities. /// @param[in] topology The topology. diff --git a/cpp/test/fem/form.cpp b/cpp/test/fem/form.cpp index 348ae15c5eb..c39aa1f1186 100644 --- a/cpp/test/fem/form.cpp +++ b/cpp/test/fem/form.cpp @@ -54,7 +54,7 @@ TEST_CASE("Create Expression/Form (mismatch of mesh geometry)", auto mesh = std::make_shared>(mesh::create_box( MPI_COMM_WORLD, {{{0.0, 0.0, 0.0}, {0.96, 4.5, 2.0}}}, {2, 4, 5}, mesh::CellType::hexahedron, - mesh::create_cell_partitioner(mesh::GhostMode::none))); + mesh::create_cell_partitioner(mesh::GhostMode::none, 2))); auto element = basix::create_element( basix::element::family::P, basix::cell::type::hexahedron, 1, basix::element::lagrange_variant::unset, diff --git a/cpp/test/mesh/distributed_mesh.cpp b/cpp/test/mesh/distributed_mesh.cpp index d51b32b929d..a19788f9d1c 100644 --- a/cpp/test/mesh/distributed_mesh.cpp +++ b/cpp/test/mesh/distributed_mesh.cpp @@ -27,7 +27,7 @@ constexpr int N = 8; [[maybe_unused]] void create_mesh_file(MPI_Comm comm) { // Create mesh using all processes and save xdmf - auto part = mesh::create_cell_partitioner(mesh::GhostMode::shared_facet); + auto part = mesh::create_cell_partitioner(mesh::GhostMode::shared_facet, 2); auto mesh = std::make_shared>( mesh::create_rectangle(comm, {{{0.0, 0.0}, {1.0, 1.0}}}, {N, N}, mesh::CellType::triangle, part)); diff --git a/cpp/test/mesh/refinement/interval.cpp b/cpp/test/mesh/refinement/interval.cpp index a3297d49de8..d5f180a1895 100644 --- a/cpp/test/mesh/refinement/interval.cpp +++ b/cpp/test/mesh/refinement/interval.cpp @@ -115,7 +115,7 @@ TEMPLATE_TEST_CASE("Interval adaptive refinement", // TODO: parent_facet auto [refined_mesh, parent_edge, parent_facet] = refinement::refine( mesh, std::span(edges), - mesh::create_cell_partitioner(mesh::GhostMode::shared_facet), + mesh::create_cell_partitioner(mesh::GhostMode::shared_facet, 2), refinement::Option::parent_cell); std::vector expected_x = { diff --git a/cpp/test/mesh/refinement/rectangle.cpp b/cpp/test/mesh/refinement/rectangle.cpp index 8c0534fa951..5686ffceb45 100644 --- a/cpp/test/mesh/refinement/rectangle.cpp +++ b/cpp/test/mesh/refinement/rectangle.cpp @@ -98,7 +98,8 @@ plotter.show() mesh.topology()->create_entities(1, std::thread::hardware_concurrency()); auto [mesh_fine, parent_cell, parent_facet] = refinement::refine( - mesh, std::nullopt, mesh::create_cell_partitioner(mesh::GhostMode::none), + mesh, std::nullopt, + mesh::create_cell_partitioner(mesh::GhostMode::none, 2), refinement::Option::parent_cell_and_facet); // vertex layout: From af1ce94c6d144b116a8f2c06d5d28f8e85ffa8c8 Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Sun, 1 Mar 2026 10:16:17 +0000 Subject: [PATCH 02/23] Graph build updates --- cpp/dolfinx/mesh/graphbuild.cpp | 185 ++++++++++++++++++++++++-------- cpp/dolfinx/mesh/graphbuild.h | 10 +- 2 files changed, 151 insertions(+), 44 deletions(-) diff --git a/cpp/dolfinx/mesh/graphbuild.cpp b/cpp/dolfinx/mesh/graphbuild.cpp index e0f3e98ccaa..527811a22d8 100644 --- a/cpp/dolfinx/mesh/graphbuild.cpp +++ b/cpp/dolfinx/mesh/graphbuild.cpp @@ -7,6 +7,7 @@ #include "graphbuild.h" #include "cell_types.h" #include +#include #include #include #include @@ -17,6 +18,7 @@ #include #include #include +#include #include #include @@ -183,8 +185,8 @@ graph::AdjacencyList compute_nonlocal_dual_graph( } std::ranges::sort(dest_to_index); - // Build list of dest ranks and count number of items (facets+cell) to - // send to each dest post office (by neighbourhood rank) + // Build list of dest ranks and count number of items (facets+cell) + // to send to each dest post office (by neighbourhood rank) for (auto it = dest_to_index.begin(); it != dest_to_index.end();) { const int neigh_rank = dest.size(); @@ -222,8 +224,7 @@ graph::AdjacencyList compute_nonlocal_dual_graph( static_cast(dest.size()) / comm_size, static_cast(src.size()) / comm_size); - // Create neighbourhood communicator for sending data to - // post offices + // Create neighbourhood communicator for sending data to post offices MPI_Comm comm_po_post; MPI_Dist_graph_create_adjacent(comm, src.size(), src.data(), MPI_UNWEIGHTED, dest.size(), dest.data(), MPI_UNWEIGHTED, @@ -286,8 +287,8 @@ graph::AdjacencyList compute_nonlocal_dual_graph( // Search for consecutive facets (-> dual graph edge between cells) // and pack into send buffer. We store for every cell the number of - // matches, the offsets of each cell and the continuous data. - // Note: 'deges' is short for dual edges. + // matches, the offsets of each cell and the continuous data. Note: + // 'deges' is short for dual edges. std::vector dedge_send_count(recv_disp.back()); std::vector dedge_send_displs(dedge_send_count.size() + 1, 0); std::vector dedge_send_data; @@ -498,8 +499,9 @@ graph::AdjacencyList compute_nonlocal_dual_graph( offset += dedge_recv_count[i]; } - // local connections are possibly introduced again by remote -> remove - // duplicates + + // Local connections are possibly introduced again by remote -> + // remove duplicates std::size_t duplicates_count = 0; for (std::size_t node = 0; node < offsets.size() - 1; node++) { @@ -530,10 +532,11 @@ std::tuple, std::vector, int, mesh::build_local_dual_graph( std::span celltypes, const std::vector>& cells, - std::optional max_facet_to_cell_links) + std::optional max_facet_to_cell_links, int num_threads) { - spdlog::info("Build local part of mesh dual graph (mixed)"); - common::Timer timer("Compute local part of mesh dual graph (mixed)"); + num_threads = 0; + spdlog::info("Build local part of mesh dual graph"); + common::Timer timer("Compute local part of mesh dual graph"); if (std::size_t ncells_local = std::accumulate(cells.begin(), cells.end(), 0, @@ -558,6 +561,8 @@ mesh::build_local_dual_graph( // number of vertices per facet -> size computations for later on // used data structures + common::Timer timer0("Compute local part of mesh dual graph: 0"); + // TODO: cell_offsets can be removed? std::vector cell_offsets{0}; cell_offsets.reserve(cells.size() + 1); @@ -585,6 +590,9 @@ mesh::build_local_dual_graph( { max = std::max(max, cell_facets.num_links(node)); }); } + timer0.stop(); + timer0.flush(); + // 2) Build a list of (all) facets, defined by sorted vertices, with // the connected cell index after the vertices. For v_ij the j-th // vertex of the i-th facet. The last index is the cell index (non @@ -594,47 +602,133 @@ mesh::build_local_dual_graph( // ⋮ ⋮ ⋮ ⋮ ⋱ ⋮ ⋮ // v_n1, v_n2, -1, -1, ..., -1, n] - const int shape1 = max_vertices_per_facet + 1; - std::vector facets; - facets.reserve(facet_count * shape1); - constexpr std::int32_t padding_value = -1; + common::Timer timer1("Compute local part of mesh dual graph: 1"); - for (std::size_t j = 0; j < cells.size(); ++j) + auto build_facets_fn + = [](int shape1, int num_cell_vertices, std::size_t cell_offset, + const graph::AdjacencyList& cell_facets, + std::span cells, std::span facets) { - const CellType& cell_type = celltypes[j]; - std::span _cells = cells[j]; - - int num_cell_vertices = mesh::cell_num_entities(cell_type, 0); - std::int32_t num_cells = _cells.size() / num_cell_vertices; - graph::AdjacencyList cell_facets - = mesh::get_entity_vertices(cell_type, tdim - 1); - - for (std::int32_t c = 0; c < num_cells; ++c) + constexpr std::int32_t padding_value = -1; + for (std::size_t c = 0; c < cells.size() / num_cell_vertices; ++c) { // Loop over cell facets - std::span v = _cells.subspan(num_cell_vertices * c, num_cell_vertices); + auto v = cells.subspan(num_cell_vertices * c, num_cell_vertices); for (int f = 0; f < cell_facets.num_nodes(); ++f) { std::span facet_vertices = cell_facets.links(f); - std::ranges::transform(facet_vertices, std::back_inserter(facets), + auto facet_c = facets.subspan( + (c * cell_facets.num_nodes() + f) * shape1, shape1); + std::ranges::transform(facet_vertices, facet_c.begin(), [v](auto idx) { return v[idx]; }); - std::sort(std::prev(facets.end(), facet_vertices.size()), facets.end()); - facets.insert(facets.end(), - max_vertices_per_facet - facet_vertices.size(), - padding_value); - facets.push_back(c + cell_offsets[j]); + + auto it = std::next(facet_c.begin(), facet_vertices.size()); + std::sort(facet_c.begin(), it); + std::fill(it, facet_c.end(), padding_value); + facet_c.back() = c + cell_offset; } } + }; + + const int shape1 = max_vertices_per_facet + 1; + std::vector facets(facet_count * shape1); + for (std::size_t j = 0; j < cells.size(); ++j) + { + CellType cell_type = celltypes[j]; + int num_cell_vertices = mesh::cell_num_entities(cell_type, 0); + graph::AdjacencyList cell_facets + = mesh::get_entity_vertices(cell_type, tdim - 1); + std::span _cells = cells[j]; + if (num_threads > 0) + { + std::vector threads; + for (int i = 0; i < num_threads; ++i) + { + auto [c0, c1] = dolfinx::MPI::local_range( + i, _cells.size() / num_cell_vertices, num_threads); + threads.emplace_back( + build_facets_fn, shape1, num_cell_vertices, cell_offsets[j] + c0, + std::cref(cell_facets), + _cells.subspan(c0 * num_cell_vertices, + (c1 - c0) * num_cell_vertices), + std::span(facets.data() + (c0 * cell_facets.num_nodes()) * shape1, + ((c1 - c0) * cell_facets.num_nodes()) * shape1)); + } + } + else + { + build_facets_fn(shape1, num_cell_vertices, cell_offsets[j], + std::cref(cell_facets), _cells, facets); + } } + timer1.stop(); + timer1.flush(); + + // const int shape1 = max_vertices_per_facet + 1; + // std::vector facets; + // facets.reserve(facet_count * shape1); + // constexpr std::int32_t padding_value = -1; + + // for (std::size_t j = 0; j < cells.size(); ++j) + // { + // const CellType& cell_type = celltypes[j]; + // std::span _cells = cells[j]; + + // int num_cell_vertices = mesh::cell_num_entities(cell_type, 0); + // std::int32_t num_cells = _cells.size() / num_cell_vertices; + // graph::AdjacencyList cell_facets + // = mesh::get_entity_vertices(cell_type, tdim - 1); + + // for (std::int32_t c = 0; c < num_cells; ++c) + // { + // // Loop over cell facets + // std::span v = _cells.subspan(num_cell_vertices * c, num_cell_vertices); + // for (int f = 0; f < cell_facets.num_nodes(); ++f) + // { + // std::span facet_vertices = cell_facets.links(f); + // std::ranges::transform(facet_vertices, std::back_inserter(facets), + // [v](auto idx) { return v[idx]; }); + // std::sort(std::prev(facets.end(), facet_vertices.size()), + // facets.end()); facets.insert(facets.end(), + // max_vertices_per_facet - facet_vertices.size(), + // padding_value); + // facets.push_back(c + cell_offsets[j]); + // } + // } + // } + // 3) Sort facets by vertex key - std::vector perm - = dolfinx::sort_by_perm(std::span(facets), shape1); + common::Timer timer3("Compute local part of mesh dual graph: 3"); + + std::vector perm; + if (num_threads > 0) + { + perm.resize(facets.size() / shape1, 0); + std::iota(perm.begin(), perm.end(), 0); + boost::sort::block_indirect_sort( + perm.begin(), perm.end(), + [facets = std::cref(facets), shape1](auto f0, auto f1) + { + auto it0 = std::next(facets.get().begin(), f0 * shape1); + auto it1 = std::next(facets.get().begin(), f1 * shape1); + return std::lexicographical_compare(it0, std::next(it0, shape1), it1, + std::next(it1, shape1)); + }, + num_threads); + } + else + perm = dolfinx::sort_by_perm(std::span(facets), shape1); + + timer3.stop(); + timer3.flush(); // 4) Iterate over sorted list of facets. Facets shared by more than // one cell lead to a graph edge to be added. Facets that are not // shared are stored as these might be shared by a cell on another // process. + common::Timer timer4("Compute local part of mesh dual graph: 4"); + std::vector unmatched_facets; std::vector local_cells; std::vector> edges; @@ -659,8 +753,7 @@ mesh::build_local_dual_graph( std::int32_t cell_count = matching_facets.size(); assert(cell_count >= 1); - if (!max_facet_to_cell_links.has_value() - or (cell_count < *max_facet_to_cell_links)) + if (!max_facet_to_cell_links or cell_count < *max_facet_to_cell_links) { // Store unmatched facets and the attached cell for (std::int32_t i = 0; i < cell_count; i++) @@ -676,12 +769,12 @@ mesh::build_local_dual_graph( // added later). In the range [it, it_next_facet), all // combinations are added. for (auto facet_a_it = it; facet_a_it != matching_facets.end(); - facet_a_it++) + ++facet_a_it) { std::span facet_a(facets.data() + *facet_a_it * shape1, shape1); std::int32_t cell_a = facet_a.back(); for (auto facet_b_it = std::next(facet_a_it); - facet_b_it != matching_facets.end(); facet_b_it++) + facet_b_it != matching_facets.end(); ++facet_b_it) { std::span facet_b(facets.data() + *facet_b_it * shape1, shape1); std::int32_t cell_b = facet_b.back(); @@ -694,13 +787,17 @@ mesh::build_local_dual_graph( } } + timer4.stop(); + timer4.flush(); + // 5) Build adjacency list data. Prepare data structure and assemble // into. Important: we have only computed one direction of the dual // edges, we add both forward and backward to the final data // structure. - std::vector num_links(cell_offsets.back(), 0); + common::Timer timer5("Compute local part of mesh dual graph: 5"); + std::vector num_links(cell_offsets.back(), 0); for (auto [a, b] : edges) { ++num_links[a]; @@ -718,6 +815,9 @@ mesh::build_local_dual_graph( data[pos[e[1]]++] = e[0]; }); + timer5.stop(); + timer5.flush(); + return {graph::AdjacencyList(std::move(data), std::move(offsets)), std::move(unmatched_facets), max_vertices_per_facet, std::move(local_cells)}; @@ -726,14 +826,15 @@ mesh::build_local_dual_graph( graph::AdjacencyList mesh::build_dual_graph(MPI_Comm comm, std::span celltypes, const std::vector>& cells, - std::optional max_facet_to_cell_links) + std::optional max_facet_to_cell_links, + int num_threads) { spdlog::info("Building mesh dual graph"); // Compute local part of dual graph (cells are graph nodes, and edges // are connections by facet) - auto [local_graph, facets, shape1, fcells] - = mesh::build_local_dual_graph(celltypes, cells, max_facet_to_cell_links); + auto [local_graph, facets, shape1, fcells] = mesh::build_local_dual_graph( + celltypes, cells, max_facet_to_cell_links, num_threads); // Extend with nonlocal edges and convert to global indices graph::AdjacencyList graph diff --git a/cpp/dolfinx/mesh/graphbuild.h b/cpp/dolfinx/mesh/graphbuild.h index 7137ca99d67..53f6716aacd 100644 --- a/cpp/dolfinx/mesh/graphbuild.h +++ b/cpp/dolfinx/mesh/graphbuild.h @@ -33,6 +33,8 @@ enum class CellType : std::int8_t; /// Passing std::nullopt (no upper bound) corresponds. /// to `max_facet_to_cell_links`=∞, i.e. every facet is considered /// unmatched. +/// @param[in] num_threads Number of threads to use. Use 0 to not launch +/// threads. /// /// @return /// 1. Local dual graph @@ -60,7 +62,8 @@ std::tuple, std::vector, int, std::vector> build_local_dual_graph(std::span celltypes, const std::vector>& cells, - std::optional max_facet_to_cell_links); + std::optional max_facet_to_cell_links, + int num_threads = 0); /// @brief Build distributed mesh dual graph (cell-cell connections via /// facets) from minimal mesh data. @@ -83,6 +86,8 @@ build_local_dual_graph(std::span celltypes, /// manifold meshes. Passing std::nullopt (no upper bound) corresponds /// to `max_facet_to_cell_links`=∞, i.e. every facet is considered /// unmatched. +/// @param[in] num_threads Number of threads to use. Use 0 to not launch +/// threads. /// /// @return The dual graph. /// @@ -96,6 +101,7 @@ build_local_dual_graph(std::span celltypes, graph::AdjacencyList build_dual_graph(MPI_Comm comm, std::span celltypes, const std::vector>& cells, - std::optional max_facet_to_cell_links); + std::optional max_facet_to_cell_links, + int num_threads = 0); } // namespace dolfinx::mesh From 89822dd02a638e1d26d163dbf0cb228fc90df06b Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Sun, 1 Mar 2026 10:22:12 +0000 Subject: [PATCH 03/23] Topology updates --- cpp/dolfinx/mesh/Topology.cpp | 122 ++++++++++++++++++++++++++-------- 1 file changed, 95 insertions(+), 27 deletions(-) diff --git a/cpp/dolfinx/mesh/Topology.cpp b/cpp/dolfinx/mesh/Topology.cpp index f6f0a31624d..69003724465 100644 --- a/cpp/dolfinx/mesh/Topology.cpp +++ b/cpp/dolfinx/mesh/Topology.cpp @@ -10,6 +10,8 @@ #include "topologycomputation.h" #include "utils.h" #include +#include +#include #include #include #include @@ -19,6 +21,7 @@ #include #include #include +#include using namespace dolfinx; using namespace dolfinx::mesh; @@ -40,7 +43,8 @@ namespace /// @return Map from global index to sharing ranks for each index in /// indices. The owner rank is the first as the first in the of ranks. graph::AdjacencyList -determine_sharing_ranks(MPI_Comm comm, std::span indices) +determine_sharing_ranks(MPI_Comm comm, std::span indices, + int num_threads) { common::Timer timer("Topology: determine shared index ownership"); @@ -63,7 +67,13 @@ determine_sharing_ranks(MPI_Comm comm, std::span indices) int dest = dolfinx::MPI::index_owner(size, idx, global_range); dest_to_index.push_back({dest, static_cast(dest_to_index.size())}); } - std::ranges::sort(dest_to_index); + if (num_threads > 0) + { + boost::sort::block_indirect_sort(dest_to_index.begin(), + dest_to_index.end(), num_threads); + } + else + std::ranges::sort(dest_to_index); } // Build list of neighbour dest ranks and count number of indices to @@ -135,10 +145,13 @@ determine_sharing_ranks(MPI_Comm comm, std::span indices) // Build {global index, pos, src} list std::vector> indices_list; - for (std::size_t p = 0; p < recv_disp0.size() - 1; ++p) - for (std::int32_t i = recv_disp0[p]; i < recv_disp0[p + 1]; ++i) - indices_list.push_back({recv_buffer0[i], i, int(p)}); - std::ranges::sort(indices_list); + { + common::Timer timer("Topology: XXXXX"); + for (std::size_t p = 0; p < recv_disp0.size() - 1; ++p) + for (std::int32_t i = recv_disp0[p]; i < recv_disp0[p + 1]; ++i) + indices_list.push_back({recv_buffer0[i], i, int(p)}); + std::ranges::sort(indices_list); + } // Find which ranks have each index std::vector num_items_per_dest1(recv_disp0.size() - 1, 0); @@ -286,7 +299,7 @@ determine_sharing_ranks(MPI_Comm comm, std::span indices) std::array, 2> vertex_ownership_groups( const std::vector>& cells_owned, const std::vector>& cells_ghost, - std::span boundary_vertices) + std::span boundary_vertices, int num_threads) { common::Timer timer("Topology: determine vertex ownership groups (owned, " "undetermined, unowned)"); @@ -300,10 +313,17 @@ std::array, 2> vertex_ownership_groups( local_vertex_set.insert(local_vertex_set.end(), c.begin(), c.end()); { - dolfinx::radix_sort(local_vertex_set); + if (num_threads > 0) + { + boost::sort::block_indirect_sort(local_vertex_set.begin(), + local_vertex_set.end(), 10); + } + else + dolfinx::radix_sort(local_vertex_set); auto [unique_end, range_end] = std::ranges::unique(local_vertex_set); local_vertex_set.erase(unique_end, range_end); } + // Build set of ghost cell vertices (attached to a ghost cell) std::vector ghost_vertex_set; ghost_vertex_set.reserve( @@ -313,10 +333,19 @@ std::array, 2> vertex_ownership_groups( ghost_vertex_set.insert(ghost_vertex_set.end(), c.begin(), c.end()); { - dolfinx::radix_sort(ghost_vertex_set); + if (num_threads > 0) + { + + boost::sort::block_indirect_sort(ghost_vertex_set.begin(), + ghost_vertex_set.end(), 10); + } + else + dolfinx::radix_sort(ghost_vertex_set); + auto [unique_end, range_end] = std::ranges::unique(ghost_vertex_set); ghost_vertex_set.erase(unique_end, range_end); } + // Build difference 1: Vertices attached only to owned cells, and // therefore owned by this rank std::vector owned_vertices; @@ -329,7 +358,7 @@ std::array, 2> vertex_ownership_groups( std::ranges::set_difference(ghost_vertex_set, local_vertex_set, std::back_inserter(unowned_vertices)); - // TODO Check this in debug mode only? +#ifndef NDEBUG // Sanity check // No vertices in unowned should also be in boundary... std::vector unowned_vertices_in_error; @@ -341,6 +370,7 @@ std::array, 2> vertex_ownership_groups( throw std::runtime_error( "Adding boundary vertices in ghost cells not allowed."); } +#endif return {std::move(owned_vertices), std::move(unowned_vertices)}; } @@ -711,12 +741,12 @@ std::vector convert_to_local_indexing( std::vector data(g.size()); if (num_threads > 0) { - std::vector threads(num_threads); + std::vector threads; for (int i = 0; i < num_threads; ++i) { auto [c0, c1] = dolfinx::MPI::local_range(i, g.size(), num_threads); - threads[i] = std::jthread(transform, std::span(data.data() + c0, c1 - c0), - g.subspan(c0, c1 - c0), global_to_local); + threads.emplace_back(transform, std::span(data.data() + c0, c1 - c0), + g.subspan(c0, c1 - c0), global_to_local); } } else @@ -1091,17 +1121,17 @@ Topology mesh::create_topology( // Create sets of owned and unowned vertices from the cell ownership // and the list of boundary vertices - auto [owned_vertices, unowned_vertices] - = vertex_ownership_groups(owned_cells, ghost_cells, boundary_vertices); + auto [owned_vertices, unowned_vertices] = vertex_ownership_groups( + owned_cells, ghost_cells, boundary_vertices, num_threads); + + timer1.stop(); + timer1.flush(); // For each vertex whose ownership needs determining, find the sharing // ranks. The first index in the list of ranks for a vertex is the // owner (as determined by determine_sharing_ranks). const graph::AdjacencyList global_vertex_to_ranks - = determine_sharing_ranks(comm, boundary_vertices); - - timer1.stop(); - timer1.flush(); + = determine_sharing_ranks(comm, boundary_vertices, num_threads); // Iterate over vertices that have 'unknown' ownership, and if flagged // as owned by determine_sharing_ranks update ownership status @@ -1134,20 +1164,58 @@ Topology mesh::create_topology( // NOTE: This block is relatively expensive // Number all owned vertices, iterating over vertices cell-wise + // std::vector local_vertex_indices(owned_vertices.size(), -1); + // { + // common::Timer timer3("Topology: 3"); + // std::int32_t v = 0; + // for (std::span cells_t : cells) + // { + // for (auto vtx : cells_t) + // { + // if (auto it = std::ranges::lower_bound(owned_vertices, vtx); + // it != owned_vertices.end() and *it == vtx) + // { + // std::size_t pos = std::distance(owned_vertices.begin(), it); + // if (local_vertex_indices[pos] < 0) + // local_vertex_indices[pos] = v++; + // } + // } + // } + // } + + // NEW: Number all owned vertices + // TODO: check if this re-ordering affects geometry retrieval + // performance std::vector local_vertex_indices(owned_vertices.size(), -1); { - common::Timer timer3("Topology: 3 "); + common::Timer timer3("Topology: 3"); + + common::Timer timer("Topology: number owned vertices"); std::int32_t v = 0; - for (std::span cells_t : cells) + for (auto _cells : cells) { - for (auto vtx : cells_t) + std::vector vertex_data; + for (auto vtx : _cells) + vertex_data.push_back(vtx); + if (num_threads == 0) + dolfinx::radix_sort(vertex_data); + else + { + boost::sort::block_indirect_sort(vertex_data.begin(), vertex_data.end(), + num_threads); + } + + auto ret = std::ranges::unique(vertex_data); + auto it0 = owned_vertices.cbegin(); + for (auto it = vertex_data.begin(); it != ret.begin(); ++it) { - if (auto it = std::ranges::lower_bound(owned_vertices, vtx); - it != owned_vertices.end() and *it == vtx) + if (*it == *it0) { - std::size_t pos = std::distance(owned_vertices.begin(), it); - if (local_vertex_indices[pos] < 0) - local_vertex_indices[pos] = v++; + std::size_t d = std::distance(owned_vertices.cbegin(), it0); + assert(d < owned_vertices.size()); + assert(local_vertex_indices[d] == -1); + local_vertex_indices[d] = v++; + ++it0; } } } From 607b3f64b645473d34bf952b6633fb700df4a80d Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Sun, 1 Mar 2026 10:44:01 +0000 Subject: [PATCH 04/23] Updates --- cpp/dolfinx/mesh/graphbuild.h | 2 +- cpp/dolfinx/mesh/utils.h | 2 +- cpp/test/mesh/branching_manifold.cpp | 10 +++++----- python/dolfinx/wrappers/mesh.cpp | 16 ++++++++++------ python/test/unit/mesh/test_dual_graph.py | 2 +- 5 files changed, 18 insertions(+), 14 deletions(-) diff --git a/cpp/dolfinx/mesh/graphbuild.h b/cpp/dolfinx/mesh/graphbuild.h index 53f6716aacd..5c45500660b 100644 --- a/cpp/dolfinx/mesh/graphbuild.h +++ b/cpp/dolfinx/mesh/graphbuild.h @@ -63,7 +63,7 @@ std::tuple, std::vector, int, build_local_dual_graph(std::span celltypes, const std::vector>& cells, std::optional max_facet_to_cell_links, - int num_threads = 0); + int num_threads); /// @brief Build distributed mesh dual graph (cell-cell connections via /// facets) from minimal mesh data. diff --git a/cpp/dolfinx/mesh/utils.h b/cpp/dolfinx/mesh/utils.h index 513a6c69531..be40cbc5b4d 100644 --- a/cpp/dolfinx/mesh/utils.h +++ b/cpp/dolfinx/mesh/utils.h @@ -276,7 +276,7 @@ create_boundary_vertices_fn(const CellReorderFunction& reorder_fn, auto [graph, unmatched_facets, max_v, _facet_attached_cells] = build_local_dual_graph(std::vector{celltypes[i]}, std::vector{cells1_v_local.back()}, - max_facet_to_cell_links); + max_facet_to_cell_links, 0); // Store unmatched_facets for current cell type facets.emplace_back(std::move(unmatched_facets), max_v); diff --git a/cpp/test/mesh/branching_manifold.cpp b/cpp/test/mesh/branching_manifold.cpp index 6ed5fab9bef..8d52fd47631 100644 --- a/cpp/test/mesh/branching_manifold.cpp +++ b/cpp/test/mesh/branching_manifold.cpp @@ -48,7 +48,7 @@ TEST_CASE("dual_graph_branching") { // default auto [dual_graph, unmatched_facets, max_vertices_per_facet, cell_data] - = mesh::build_local_dual_graph(celltypes, {cells}, 2); + = mesh::build_local_dual_graph(celltypes, {cells}, 2, 0); CHECK(dual_graph.num_nodes() == 4); @@ -83,7 +83,7 @@ TEST_CASE("dual_graph_branching") // max_facet_to_cell_links = 3 // Note: additionally facet (2) is now considered unmatched auto [dual_graph, unmatched_facets, max_vertices_per_facet, cell_data] - = mesh::build_local_dual_graph(celltypes, {cells}, 3); + = mesh::build_local_dual_graph(celltypes, {cells}, 3, 0); CHECK(dual_graph.num_nodes() == 4); @@ -125,7 +125,7 @@ TEST_CASE("dual_graph_branching") auto [dual_graph, unmatched_facets, max_vertices_per_facet, cell_data] = mesh::build_local_dual_graph(celltypes, {cells}, - max_facet_to_cell_links); + max_facet_to_cell_links, 0); CHECK(dual_graph.num_nodes() == 4); @@ -184,7 +184,7 @@ TEST_CASE("dual_graph_self_dual") { auto [dual_graph, unmatched_facets, max_vertices_per_facet, cell_data] = mesh::build_local_dual_graph(celltypes, {cells}, - max_facet_to_cell_links); + max_facet_to_cell_links, 0); CHECK(max_vertices_per_facet == 1); CHECK(dual_graph.num_nodes() == 3); @@ -254,7 +254,7 @@ TEST_CASE("dual_graph_branching_parallel") // Check local dual graphs. auto [dual_graph, unmatched_facets, max_vertices_per_facet, cell_data] - = mesh::build_local_dual_graph(celltypes, {cells}, 3); + = mesh::build_local_dual_graph(celltypes, {cells}, 3, 0); CHECK(max_vertices_per_facet == 1); diff --git a/python/dolfinx/wrappers/mesh.cpp b/python/dolfinx/wrappers/mesh.cpp index d86035de9bd..c6628c9f3a0 100644 --- a/python/dolfinx/wrappers/mesh.cpp +++ b/python/dolfinx/wrappers/mesh.cpp @@ -105,14 +105,16 @@ void mesh(nb::module_& m) "build_dual_graph", [](const MPICommWrapper comm, dolfinx::mesh::CellType cell_type, const dolfinx::graph::AdjacencyList& cells, - std::optional max_facet_to_cell_links) + std::optional max_facet_to_cell_links, int num_threads) { std::vector c = {cell_type}; return dolfinx::mesh::build_dual_graph( - comm.get(), std::span{c}, {cells.array()}, max_facet_to_cell_links); + comm.get(), std::span{c}, {cells.array()}, max_facet_to_cell_links, + num_threads); }, nb::arg("comm"), nb::arg("cell_type"), nb::arg("cells"), - nb::arg("max_facet_to_cell_links").none(), "Build dual graph for cells"); + nb::arg("max_facet_to_cell_links").none(), nb::arg("num_threads"), + "Build dual graph for cells"); m.def( "build_dual_graph", @@ -120,7 +122,7 @@ void mesh(nb::module_& m) std::vector& cell_types, const std::vector< nb::ndarray, nb::c_contig>>& cells, - std::optional max_facet_to_cell_links) + std::optional max_facet_to_cell_links, int num_threads) { std::vector> cell_span(cells.size()); for (std::size_t i = 0; i < cells.size(); ++i) @@ -129,10 +131,12 @@ void mesh(nb::module_& m) = std::span(cells[i].data(), cells[i].size()); } return dolfinx::mesh::build_dual_graph( - comm.get(), cell_types, cell_span, max_facet_to_cell_links); + comm.get(), cell_types, cell_span, max_facet_to_cell_links, + num_threads); }, nb::arg("comm"), nb::arg("cell_types"), nb::arg("cells"), - nb::arg("max_facet_to_cell_links").none(), "Build dual graph for cells"); + nb::arg("max_facet_to_cell_links").none(), nb::arg("num_threads"), + "Build dual graph for cells"); // dolfinx::mesh::GhostMode enums nb::enum_(m, "GhostMode") diff --git a/python/test/unit/mesh/test_dual_graph.py b/python/test/unit/mesh/test_dual_graph.py index a4ce550d949..42f10013245 100644 --- a/python/test/unit/mesh/test_dual_graph.py +++ b/python/test/unit/mesh/test_dual_graph.py @@ -26,7 +26,7 @@ def test_dgrsph_1d(): # Circular chain of interval cells cells = [[n0, n0 + 1], [n0 + 1, n0 + 2], [n0 + 2, x]] w = mesh.build_dual_graph( - MPI.COMM_WORLD, mesh.CellType.interval, to_adj(cells, np.int64)._cpp_object, 2 + MPI.COMM_WORLD, mesh.CellType.interval, to_adj(cells, np.int64)._cpp_object, 2, 0 ) assert w.num_nodes == 3 for i in range(w.num_nodes): From eeee2eee9c98406af743e83c483e1e46557e12d9 Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Sat, 18 Jul 2026 14:38:50 +0100 Subject: [PATCH 05/23] Build fix --- cpp/dolfinx/mesh/Topology.cpp | 3 ++- cpp/dolfinx/mesh/graphbuild.cpp | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/cpp/dolfinx/mesh/Topology.cpp b/cpp/dolfinx/mesh/Topology.cpp index 22f3458a079..6170d3be35b 100644 --- a/cpp/dolfinx/mesh/Topology.cpp +++ b/cpp/dolfinx/mesh/Topology.cpp @@ -13,6 +13,7 @@ #include #include #include +#include #include #include #include @@ -744,7 +745,7 @@ std::vector convert_to_local_indexing( std::vector threads; for (int i = 0; i < num_threads; ++i) { - auto [c0, c1] = dolfinx::MPI::local_range(i, g.size(), num_threads); + auto [c0, c1] = dolfinx::common::local_range(i, g.size(), num_threads); threads.emplace_back(transform, std::span(data.data() + c0, c1 - c0), g.subspan(c0, c1 - c0), global_to_local); } diff --git a/cpp/dolfinx/mesh/graphbuild.cpp b/cpp/dolfinx/mesh/graphbuild.cpp index 527811a22d8..c38dd57db80 100644 --- a/cpp/dolfinx/mesh/graphbuild.cpp +++ b/cpp/dolfinx/mesh/graphbuild.cpp @@ -10,6 +10,7 @@ #include #include #include +#include #include #include #include @@ -644,7 +645,7 @@ mesh::build_local_dual_graph( std::vector threads; for (int i = 0; i < num_threads; ++i) { - auto [c0, c1] = dolfinx::MPI::local_range( + auto [c0, c1] = dolfinx::common::local_range( i, _cells.size() / num_cell_vertices, num_threads); threads.emplace_back( build_facets_fn, shape1, num_cell_vertices, cell_offsets[j] + c0, From 420c159c09bf41a0291194e542ab2a6842a6d93a Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Sat, 18 Jul 2026 17:52:32 +0100 Subject: [PATCH 06/23] Fix broken vertex numbering in create_topology The experimental replacement for the cell-wise owned-vertex numbering loop broke correctness: it reset the merge iterator into the globally sorted owned_vertices array back to cbegin() for every cell, and never advanced it when a mismatch occurred, so vertices were mismatched to the wrong position. This caused CHECK_adjacency_list_equal failures in "Interval adaptive refinement" (interval.cpp) with cyclically-shifted connectivity. Restore the previously working, correctness-verified implementation (the old code was left commented out above the broken block). Co-Authored-By: Claude Sonnet 5 --- cpp/dolfinx/mesh/Topology.cpp | 54 ++++++----------------------------- 1 file changed, 8 insertions(+), 46 deletions(-) diff --git a/cpp/dolfinx/mesh/Topology.cpp b/cpp/dolfinx/mesh/Topology.cpp index 6170d3be35b..3b1aa8ea735 100644 --- a/cpp/dolfinx/mesh/Topology.cpp +++ b/cpp/dolfinx/mesh/Topology.cpp @@ -1165,58 +1165,20 @@ Topology mesh::create_topology( // NOTE: This block is relatively expensive // Number all owned vertices, iterating over vertices cell-wise - // std::vector local_vertex_indices(owned_vertices.size(), -1); - // { - // common::Timer timer3("Topology: 3"); - // std::int32_t v = 0; - // for (std::span cells_t : cells) - // { - // for (auto vtx : cells_t) - // { - // if (auto it = std::ranges::lower_bound(owned_vertices, vtx); - // it != owned_vertices.end() and *it == vtx) - // { - // std::size_t pos = std::distance(owned_vertices.begin(), it); - // if (local_vertex_indices[pos] < 0) - // local_vertex_indices[pos] = v++; - // } - // } - // } - // } - - // NEW: Number all owned vertices - // TODO: check if this re-ordering affects geometry retrieval - // performance std::vector local_vertex_indices(owned_vertices.size(), -1); { - common::Timer timer3("Topology: 3"); - - common::Timer timer("Topology: number owned vertices"); + common::Timer timer3("Topology: number owned vertices"); std::int32_t v = 0; - for (auto _cells : cells) + for (std::span cells_t : cells) { - std::vector vertex_data; - for (auto vtx : _cells) - vertex_data.push_back(vtx); - if (num_threads == 0) - dolfinx::radix_sort(vertex_data); - else - { - boost::sort::block_indirect_sort(vertex_data.begin(), vertex_data.end(), - num_threads); - } - - auto ret = std::ranges::unique(vertex_data); - auto it0 = owned_vertices.cbegin(); - for (auto it = vertex_data.begin(); it != ret.begin(); ++it) + for (auto vtx : cells_t) { - if (*it == *it0) + if (auto it = std::ranges::lower_bound(owned_vertices, vtx); + it != owned_vertices.end() and *it == vtx) { - std::size_t d = std::distance(owned_vertices.cbegin(), it0); - assert(d < owned_vertices.size()); - assert(local_vertex_indices[d] == -1); - local_vertex_indices[d] = v++; - ++it0; + std::size_t pos = std::distance(owned_vertices.begin(), it); + if (local_vertex_indices[pos] < 0) + local_vertex_indices[pos] = v++; } } } From 274d6520a86054ddc32f08b4422e0f3c0d8b6b28 Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Sun, 19 Jul 2026 10:55:44 +0100 Subject: [PATCH 07/23] Pass num_threads --- cpp/dolfinx/mesh/Topology.cpp | 3 +-- cpp/dolfinx/mesh/graphbuild.cpp | 36 ++------------------------------- cpp/dolfinx/mesh/utils.h | 8 ++++---- 3 files changed, 7 insertions(+), 40 deletions(-) diff --git a/cpp/dolfinx/mesh/Topology.cpp b/cpp/dolfinx/mesh/Topology.cpp index 3b1aa8ea735..1406cce5305 100644 --- a/cpp/dolfinx/mesh/Topology.cpp +++ b/cpp/dolfinx/mesh/Topology.cpp @@ -336,9 +336,8 @@ std::array, 2> vertex_ownership_groups( { if (num_threads > 0) { - boost::sort::block_indirect_sort(ghost_vertex_set.begin(), - ghost_vertex_set.end(), 10); + ghost_vertex_set.end(), num_threads); } else dolfinx::radix_sort(ghost_vertex_set); diff --git a/cpp/dolfinx/mesh/graphbuild.cpp b/cpp/dolfinx/mesh/graphbuild.cpp index c38dd57db80..339be7e5fa4 100644 --- a/cpp/dolfinx/mesh/graphbuild.cpp +++ b/cpp/dolfinx/mesh/graphbuild.cpp @@ -535,7 +535,6 @@ mesh::build_local_dual_graph( const std::vector>& cells, std::optional max_facet_to_cell_links, int num_threads) { - num_threads = 0; spdlog::info("Build local part of mesh dual graph"); common::Timer timer("Compute local part of mesh dual graph"); @@ -666,39 +665,6 @@ mesh::build_local_dual_graph( timer1.stop(); timer1.flush(); - // const int shape1 = max_vertices_per_facet + 1; - // std::vector facets; - // facets.reserve(facet_count * shape1); - // constexpr std::int32_t padding_value = -1; - - // for (std::size_t j = 0; j < cells.size(); ++j) - // { - // const CellType& cell_type = celltypes[j]; - // std::span _cells = cells[j]; - - // int num_cell_vertices = mesh::cell_num_entities(cell_type, 0); - // std::int32_t num_cells = _cells.size() / num_cell_vertices; - // graph::AdjacencyList cell_facets - // = mesh::get_entity_vertices(cell_type, tdim - 1); - - // for (std::int32_t c = 0; c < num_cells; ++c) - // { - // // Loop over cell facets - // std::span v = _cells.subspan(num_cell_vertices * c, num_cell_vertices); - // for (int f = 0; f < cell_facets.num_nodes(); ++f) - // { - // std::span facet_vertices = cell_facets.links(f); - // std::ranges::transform(facet_vertices, std::back_inserter(facets), - // [v](auto idx) { return v[idx]; }); - // std::sort(std::prev(facets.end(), facet_vertices.size()), - // facets.end()); facets.insert(facets.end(), - // max_vertices_per_facet - facet_vertices.size(), - // padding_value); - // facets.push_back(c + cell_offsets[j]); - // } - // } - // } - // 3) Sort facets by vertex key common::Timer timer3("Compute local part of mesh dual graph: 3"); @@ -719,7 +685,9 @@ mesh::build_local_dual_graph( num_threads); } else + { perm = dolfinx::sort_by_perm(std::span(facets), shape1); + } timer3.stop(); timer3.flush(); diff --git a/cpp/dolfinx/mesh/utils.h b/cpp/dolfinx/mesh/utils.h index 90bf2397329..b3b876eb47b 100644 --- a/cpp/dolfinx/mesh/utils.h +++ b/cpp/dolfinx/mesh/utils.h @@ -253,8 +253,8 @@ create_boundary_vertices_fn(const CellReorderFunction& reorder_fn, const std::vector>& ghost_owners, std::vector>& cells, std::vector>& cells_v, - std::vector>& original_idx) - -> std::vector + std::vector>& original_idx, + int num_threads) -> std::vector { // Build local dual graph for owned cells to (i) get list of vertices // on the process boundary and (ii) apply re-ordering to cells for @@ -279,7 +279,7 @@ create_boundary_vertices_fn(const CellReorderFunction& reorder_fn, auto [graph, unmatched_facets, max_v, _facet_attached_cells] = build_local_dual_graph(std::vector{celltypes[i]}, std::vector{cells1_v_local.back()}, - max_facet_to_cell_links, 0); + max_facet_to_cell_links, num_threads); // Store unmatched_facets for current cell type facets.emplace_back(std::move(unmatched_facets), max_v); @@ -1158,7 +1158,7 @@ Mesh> create_mesh( auto boundary_v_fn = create_boundary_vertices_fn(reorder_fn, max_facet_to_cell_links); const std::vector boundary_v = boundary_v_fn( - celltypes, doflayouts, ghost_owners, cells1, cells1_v, original_idx1); + celltypes, doflayouts, ghost_owners, cells1, cells1_v, original_idx1, 0); spdlog::debug("Got {} boundary vertices", boundary_v.size()); From 62d7cd245d22f26ca59999175dede6d534b843d1 Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Sun, 19 Jul 2026 11:54:04 +0100 Subject: [PATCH 08/23] Pass num threads --- cpp/dolfinx/io/VTKHDF.h | 2 +- cpp/dolfinx/mesh/generation.h | 24 +++++++++---------- cpp/dolfinx/mesh/utils.h | 17 +++++++------ cpp/dolfinx/refinement/refine.h | 2 +- cpp/dolfinx/refinement/uniform.cpp | 2 +- .../dolfinx/wrappers/dolfinx_wrappers/mesh.h | 14 +++++------ 6 files changed, 32 insertions(+), 29 deletions(-) diff --git a/cpp/dolfinx/io/VTKHDF.h b/cpp/dolfinx/io/VTKHDF.h index 77184e10ae2..488751b6ee2 100644 --- a/cpp/dolfinx/io/VTKHDF.h +++ b/cpp/dolfinx/io/VTKHDF.h @@ -465,6 +465,6 @@ mesh::Mesh read_mesh(MPI_Comm comm, const std::filesystem::path& filename, cells_local.end()); return mesh::create_mesh(comm, comm, cells_span, coordinate_elements, comm, points_pruned, {(std::size_t)x_shape[0], gdim}, part, - max_facet_to_cell_links); + max_facet_to_cell_links, 0); } } // namespace dolfinx::io::VTKHDF diff --git a/cpp/dolfinx/mesh/generation.h b/cpp/dolfinx/mesh/generation.h index 5d0d8a68bfa..49ccd23b59b 100644 --- a/cpp/dolfinx/mesh/generation.h +++ b/cpp/dolfinx/mesh/generation.h @@ -285,20 +285,20 @@ create_interval(MPI_Comm comm, std::int64_t n, std::array p, if (gdim == 1) { return create_mesh(comm, MPI_COMM_SELF, cells, element, MPI_COMM_SELF, - x1d, {npts, 1}, partitioner, 2, reorder_fn); + x1d, {npts, 1}, partitioner, 2, 0, reorder_fn); } std::vector x(npts * gdim, T(0)); for (std::size_t i = 0; i < npts; i++) x[i * gdim] = x1d[i]; return create_mesh(comm, MPI_COMM_SELF, cells, element, MPI_COMM_SELF, x, {npts, static_cast(gdim)}, partitioner, 2, - reorder_fn); + 0, reorder_fn); } else { return create_mesh(comm, MPI_COMM_NULL, {}, element, MPI_COMM_NULL, std::vector{}, {0, static_cast(gdim)}, - partitioner, 2, reorder_fn); + partitioner, 2, 0, reorder_fn); } } @@ -422,7 +422,7 @@ Mesh build_tet(MPI_Comm comm, MPI_Comm subcomm, } return create_mesh(comm, subcomm, cells, element, subcomm, x, - {x.size() / 3, 3}, partitioner, 2, reorder_fn); + {x.size() / 3, 3}, partitioner, 2, 0, reorder_fn); } template @@ -467,7 +467,7 @@ build_hex(MPI_Comm comm, MPI_Comm subcomm, std::array, 2> p, } return create_mesh(comm, subcomm, cells, element, subcomm, x, - {x.size() / 3, 3}, partitioner, 2, reorder_fn); + {x.size() / 3, 3}, partitioner, 2, 0, reorder_fn); } template @@ -515,7 +515,7 @@ Mesh build_prism(MPI_Comm comm, MPI_Comm subcomm, } return create_mesh(comm, subcomm, cells, element, subcomm, x, - {x.size() / 3, 3}, partitioner, 2, reorder_fn); + {x.size() / 3, 3}, partitioner, 2, 0, reorder_fn); } template @@ -672,7 +672,7 @@ Mesh build_tri(MPI_Comm comm, std::array, 2> p, if (gdim == 2) { return create_mesh(comm, MPI_COMM_SELF, cells, element, MPI_COMM_SELF, x, - {npts, 2}, partitioner, 2, reorder_fn); + {npts, 2}, partitioner, 2, 0, reorder_fn); } std::vector xg(npts * gdim, T(0)); for (std::size_t i = 0; i < npts; i++) @@ -682,13 +682,13 @@ Mesh build_tri(MPI_Comm comm, std::array, 2> p, } return create_mesh(comm, MPI_COMM_SELF, cells, element, MPI_COMM_SELF, xg, {npts, static_cast(gdim)}, partitioner, 2, - reorder_fn); + 0, reorder_fn); } else { return create_mesh(comm, MPI_COMM_NULL, {}, element, MPI_COMM_NULL, std::vector{}, {0, static_cast(gdim)}, - partitioner, 2, reorder_fn); + partitioner, 2, 0, reorder_fn); } } @@ -738,7 +738,7 @@ Mesh build_quad(MPI_Comm comm, std::array, 2> p, if (gdim == 2) { return create_mesh(comm, MPI_COMM_SELF, cells, element, MPI_COMM_SELF, x, - {npts, 2}, partitioner, 2, reorder_fn); + {npts, 2}, partitioner, 2, 0, reorder_fn); } std::vector xg(npts * gdim, T(0)); for (std::size_t i = 0; i < npts; i++) @@ -748,13 +748,13 @@ Mesh build_quad(MPI_Comm comm, std::array, 2> p, } return create_mesh(comm, MPI_COMM_SELF, cells, element, MPI_COMM_SELF, xg, {npts, static_cast(gdim)}, partitioner, 2, - reorder_fn); + 0, reorder_fn); } else { return create_mesh(comm, MPI_COMM_NULL, {}, element, MPI_COMM_NULL, std::vector{}, {0, static_cast(gdim)}, - partitioner, 2, reorder_fn); + partitioner, 2, 0, reorder_fn); } } } // namespace impl diff --git a/cpp/dolfinx/mesh/utils.h b/cpp/dolfinx/mesh/utils.h index b3b876eb47b..9ca2974a595 100644 --- a/cpp/dolfinx/mesh/utils.h +++ b/cpp/dolfinx/mesh/utils.h @@ -1039,6 +1039,8 @@ compute_incident_entities(const Topology& topology, /// redistributed. /// @param[in] max_facet_to_cell_links Bound on the number of cells a /// facet can be connected to. +/// @param[in] num_threads Number threads to use in mesh construction. +/// It 0, no threads are launched. /// @param[in] reorder_fn Function that reorders (locally) cells that /// are owned by this process. /// @return A mesh distributed on the communicator `comm`. @@ -1050,7 +1052,7 @@ Mesh> create_mesh( typename std::remove_reference_t>>& elements, MPI_Comm commg, const U& x, std::array xshape, const CellPartitionFunction& partitioner, - std::optional max_facet_to_cell_links, + std::optional max_facet_to_cell_links, int num_threads, const CellReorderFunction& reorder_fn = graph::reorder_rcm) { assert(cells.size() == elements.size()); @@ -1157,8 +1159,9 @@ Mesh> create_mesh( auto boundary_v_fn = create_boundary_vertices_fn(reorder_fn, max_facet_to_cell_links); - const std::vector boundary_v = boundary_v_fn( - celltypes, doflayouts, ghost_owners, cells1, cells1_v, original_idx1, 0); + const std::vector boundary_v + = boundary_v_fn(celltypes, doflayouts, ghost_owners, cells1, cells1_v, + original_idx1, num_threads); spdlog::debug("Got {} boundary vertices", boundary_v.size()); @@ -1264,12 +1267,12 @@ Mesh> create_mesh( typename std::remove_reference_t>& element, MPI_Comm commg, const U& x, std::array xshape, const CellPartitionFunction& partitioner, - std::optional max_facet_to_cell_links, + std::optional max_facet_to_cell_links, int num_threads, const CellReorderFunction& reorder_fn = graph::reorder_rcm) { return create_mesh(comm, commt, std::vector{cells}, std::vector{element}, commg, x, xshape, partitioner, max_facet_to_cell_links, - reorder_fn); + num_threads, reorder_fn); } /// @brief Create a distributed mesh from mesh data using the default @@ -1303,14 +1306,14 @@ create_mesh(MPI_Comm comm, std::span cells, if (dolfinx::MPI::size(comm) == 1) { return create_mesh(comm, comm, std::vector{cells}, std::vector{elements}, - comm, x, xshape, nullptr, max_facet_to_cell_links); + comm, x, xshape, nullptr, max_facet_to_cell_links, 0); } else { return create_mesh( comm, comm, std::vector{cells}, std::vector{elements}, comm, x, xshape, create_cell_partitioner(ghost_mode, max_facet_to_cell_links), - max_facet_to_cell_links); + max_facet_to_cell_links, 0); } } diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index cce82d36ee8..dd9a3941e10 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -152,7 +152,7 @@ refine(const mesh::Mesh& mesh, mesh::Mesh mesh1 = mesh::create_mesh( mesh.comm(), mesh.comm(), cell_adj.array(), mesh.geometry().cmaps().front(), mesh.comm(), new_vertex_coords, xshape, - std::get(partitioner), 2); + std::get(partitioner), 2, 0); // Report the number of refined cells const int D = topology->dim(); diff --git a/cpp/dolfinx/refinement/uniform.cpp b/cpp/dolfinx/refinement/uniform.cpp index 819393d30c8..1deacaf8886 100644 --- a/cpp/dolfinx/refinement/uniform.cpp +++ b/cpp/dolfinx/refinement/uniform.cpp @@ -301,7 +301,7 @@ refinement::uniform_refine(const mesh::Mesh& mesh, geometry_cmaps.push_back(cm); mesh::Mesh new_mesh = mesh::create_mesh( mesh.comm(), mesh.comm(), topo_span, geometry_cmaps, mesh.comm(), new_x, - {new_x.size() / 3, 3}, partitioner, 2); + {new_x.size() / 3, 3}, partitioner, 2, 0); return new_mesh; } diff --git a/python/dolfinx/wrappers/dolfinx_wrappers/mesh.h b/python/dolfinx/wrappers/dolfinx_wrappers/mesh.h index ad1b315cc1f..79d5350972b 100644 --- a/python/dolfinx/wrappers/dolfinx_wrappers/mesh.h +++ b/python/dolfinx/wrappers/dolfinx_wrappers/mesh.h @@ -314,7 +314,7 @@ void declare_mesh(nb::module_& m, std::string type) const std::vector>& elements, nb::ndarray x, const part::impl::PythonCellPartitionFunction& p, - std::optional max_facet_to_cell_links) + std::optional max_facet_to_cell_links, int num_threads) { std::size_t shape1 = x.ndim() == 1 ? 1 : x.shape(1); @@ -326,12 +326,12 @@ void declare_mesh(nb::module_& m, std::string type) return dolfinx::mesh::create_mesh( comm.get(), comm.get(), cells, elements, comm.get(), std::span(x.data(), x.size()), {x.shape(0), shape1}, - part::impl::create_cell_partitioner_cpp(p), - max_facet_to_cell_links); + part::impl::create_cell_partitioner_cpp(p), max_facet_to_cell_links, + num_threads); }, nb::arg("comm"), nb::arg("cells"), nb::arg("elements"), nb::arg("x").noconvert(), nb::arg("partitioner").none(), - nb::arg("max_facet_to_cell_links").none(), + nb::arg("max_facet_to_cell_links").none(), nb::arg("num_threads"), "Helper function for creating a mixed topology mesh."); m.def( @@ -341,18 +341,18 @@ void declare_mesh(nb::module_& m, std::string type) const dolfinx::fem::CoordinateElement& element, nb::ndarray x, const part::impl::PythonCellPartitionFunction& p, - std::optional max_facet_to_cell_links) + std::optional max_facet_to_cell_links, int num_threads) { std::size_t shape1 = x.ndim() == 1 ? 1 : x.shape(1); return dolfinx::mesh::create_mesh( comm.get(), comm.get(), std::span(cells.data(), cells.size()), element, comm.get(), std::span(x.data(), x.size()), {x.shape(0), shape1}, part::impl::create_cell_partitioner_cpp(p), - max_facet_to_cell_links); + max_facet_to_cell_links, num_threads); }, nb::arg("comm"), nb::arg("cells"), nb::arg("element"), nb::arg("x").noconvert(), nb::arg("partitioner").none(), - nb::arg("max_facet_to_cell_links").none(), + nb::arg("max_facet_to_cell_links").none(), nb::arg("num_threads"), "Helper function for creating meshes."); m.def( "create_submesh", From 7bb3412107056c91399561efef53699c8985835b Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Sun, 19 Jul 2026 12:32:58 +0100 Subject: [PATCH 09/23] Fixes --- cpp/test/mesh/distributed_mesh.cpp | 2 +- cpp/test/mesh/refinement/interval.cpp | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/cpp/test/mesh/distributed_mesh.cpp b/cpp/test/mesh/distributed_mesh.cpp index a19788f9d1c..6eeb1659228 100644 --- a/cpp/test/mesh/distributed_mesh.cpp +++ b/cpp/test/mesh/distributed_mesh.cpp @@ -139,7 +139,7 @@ void test_distributed_mesh(const mesh::CellPartitionFunction& partitioner) // Build mesh mesh::Mesh mesh = mesh::create_mesh(comm, subset_comm, cells, cmap, comm, x, - xshape, partitioner, 2); + xshape, partitioner, 2, 0); auto t = mesh.topology(); int tdim = t->dim(); CHECK(t->index_map(tdim)->size_global() == 2 * N * N); diff --git a/cpp/test/mesh/refinement/interval.cpp b/cpp/test/mesh/refinement/interval.cpp index 4650fc7c59b..59c49d87024 100644 --- a/cpp/test/mesh/refinement/interval.cpp +++ b/cpp/test/mesh/refinement/interval.cpp @@ -61,7 +61,7 @@ mesh::Mesh create_3_vertex_interval_mesh() return mesh::create_mesh( MPI_COMM_SELF, MPI_COMM_SELF, cells, element, MPI_COMM_SELF, x, {x.size() / 3, 3}, - mesh::create_cell_partitioner(mesh::GhostMode::none, 2), 2); + mesh::create_cell_partitioner(mesh::GhostMode::none, 2), 2, 0); } TEMPLATE_TEST_CASE("Interval uniform refinement", @@ -152,7 +152,7 @@ TEMPLATE_TEST_CASE("Interval Refinement (parallel)", MPI_Comm commt = rank == 0 ? MPI_COMM_SELF : MPI_COMM_NULL; return mesh::create_mesh(MPI_COMM_WORLD, commt, cells, element, commt, x, - {x.size() / 3, 3}, partitioner, 2); + {x.size() / 3, 3}, partitioner, 2, 0); }; mesh::Mesh mesh = create_mesh(); From 7044ed0666e2c505ebe6e30412c03f7469a1323f Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Sun, 19 Jul 2026 12:47:36 +0100 Subject: [PATCH 10/23] Fixes --- python/dolfinx/io/gmsh.py | 1 + python/dolfinx/io/utils.py | 1 + python/dolfinx/mesh.py | 5 ++++- python/test/unit/io/test_vtkhdf.py | 1 + python/test/unit/mesh/test_mixed_topology.py | 2 +- 5 files changed, 8 insertions(+), 2 deletions(-) diff --git a/python/dolfinx/io/gmsh.py b/python/dolfinx/io/gmsh.py index d6074c07218..75c2ac55dc6 100644 --- a/python/dolfinx/io/gmsh.py +++ b/python/dolfinx/io/gmsh.py @@ -441,6 +441,7 @@ def model_to_mesh( x[:, :gdim].astype(dtype).copy(), partitioner, max_facet_to_cell_links, + 0, ) mesh = Mesh(cpp_mesh, None) diff --git a/python/dolfinx/io/utils.py b/python/dolfinx/io/utils.py index a6c4276d295..82381feb9d7 100644 --- a/python/dolfinx/io/utils.py +++ b/python/dolfinx/io/utils.py @@ -273,6 +273,7 @@ def read_mesh( x, create_cell_partitioner(ghost_mode, max_facet_to_cell_links), # type: ignore max_facet_to_cell_links, + 0, ) msh.name = name domain = ufl.Mesh(basix_el) diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index 8afffd43e28..1402f1b6bb6 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -789,6 +789,7 @@ def create_mesh( x: npt.NDArray[np.floating], partitioner: Callable | None = None, max_facet_to_cell_links: int = 2, + num_threads: int = 0, ) -> Mesh: """Create a mesh from topology and geometry arrays. @@ -804,6 +805,8 @@ def create_mesh( cells across MPI ranks. max_facet_to_cell_links: Maximum number of cells a facet can be connected to. + num_threads: Number of threads to use tp build mesh. Set to 0 to + not spawn threads. Note: If required, the coordinates ``x`` will be cast to the same @@ -856,7 +859,7 @@ def create_mesh( x = np.asarray(x, dtype=dtype, order="C") cells = np.asarray(cells, dtype=np.int64, order="C") msh: _cpp.mesh.Mesh_float32 | _cpp.mesh.Mesh_float64 = _cpp.mesh.create_mesh( - comm, cells, cmap._cpp_object, x, partitioner, max_facet_to_cell_links + comm, cells, cmap._cpp_object, x, partitioner, max_facet_to_cell_links, num_threads ) return Mesh(msh, domain) # type: ignore diff --git a/python/test/unit/io/test_vtkhdf.py b/python/test/unit/io/test_vtkhdf.py index 91f655de2e6..8143a581d54 100644 --- a/python/test/unit/io/test_vtkhdf.py +++ b/python/test/unit/io/test_vtkhdf.py @@ -91,6 +91,7 @@ def test_read_write_higher_order(): geom, part, max_cells_per_facet, + 0, ) py_mesh = Mesh(mesh, None) diff --git a/python/test/unit/mesh/test_mixed_topology.py b/python/test/unit/mesh/test_mixed_topology.py index 96e9d90889b..0316a8aef48 100644 --- a/python/test/unit/mesh/test_mixed_topology.py +++ b/python/test/unit/mesh/test_mixed_topology.py @@ -344,7 +344,7 @@ def test_locate_entities(dtype): comm = MPI.COMM_WORLD max_cells_per_facet = 2 mesh = create_mesh( - comm, cells, [hexahedron._cpp_object, prism._cpp_object], geom, part, max_cells_per_facet + comm, cells, [hexahedron._cpp_object, prism._cpp_object], geom, part, max_cells_per_facet, 0 ) fdim = mesh.topology.dim - 1 From e2a1da134a02ac826317d23cfaaf40d25db251fe Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Sun, 19 Jul 2026 13:37:07 +0100 Subject: [PATCH 11/23] Fix --- python/test/unit/conftest.py | 1 + python/test/unit/mesh/test_create_mixed_mesh.py | 1 + 2 files changed, 2 insertions(+) diff --git a/python/test/unit/conftest.py b/python/test/unit/conftest.py index 93f4758f75f..30510e3472c 100644 --- a/python/test/unit/conftest.py +++ b/python/test/unit/conftest.py @@ -120,4 +120,5 @@ def mixed_topology_mesh(): geomx, part, max_cells_per_facet, + 0, ) diff --git a/python/test/unit/mesh/test_create_mixed_mesh.py b/python/test/unit/mesh/test_create_mixed_mesh.py index 7ea83c967a8..b2038f5ad3a 100644 --- a/python/test/unit/mesh/test_create_mixed_mesh.py +++ b/python/test/unit/mesh/test_create_mixed_mesh.py @@ -91,6 +91,7 @@ def test_create_mixed_mesh(dtype): geomx, part, max_cells_per_facet, + 0, ) entity_types = mesh.topology.entity_types[3] From 766fec19869310dd46a16ee8746d143c1e6f25d7 Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Sun, 19 Jul 2026 13:38:02 +0100 Subject: [PATCH 12/23] Demo fix --- python/demo/demo_mixed-topology.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/demo/demo_mixed-topology.py b/python/demo/demo_mixed-topology.py index f725b90465d..2b6a3efc652 100644 --- a/python/demo/demo_mixed-topology.py +++ b/python/demo/demo_mixed-topology.py @@ -112,7 +112,7 @@ part = create_cell_partitioner(GhostMode.none, 2) # type: ignore mesh = create_mesh( - MPI.COMM_WORLD, cells_np, [hexahedron._cpp_object, prism._cpp_object], geomx, part, 2 + MPI.COMM_WORLD, cells_np, [hexahedron._cpp_object, prism._cpp_object], geomx, part, 2, 0 ) # - From 3f82db30ab371ce8db9d094ef1f23ecdcd08eb58 Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Sun, 19 Jul 2026 17:02:22 +0100 Subject: [PATCH 13/23] Doc fix --- cpp/dolfinx/mesh/utils.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/cpp/dolfinx/mesh/utils.h b/cpp/dolfinx/mesh/utils.h index 9ca2974a595..c0b64ddb1c7 100644 --- a/cpp/dolfinx/mesh/utils.h +++ b/cpp/dolfinx/mesh/utils.h @@ -1257,6 +1257,8 @@ Mesh> create_mesh( /// rank for each cell. If not callable, cells are not redistributed. /// @param[in] max_facet_to_cell_links Bound on the number of cells a /// facet can be connected to. +/// @param[in] num_threads Number threads to use in mesh construction. If 0, no +/// threads are launched. /// @param[in] reorder_fn Function that reorders (locally) cells that /// are owned by this process. /// @return A mesh distributed on the communicator `comm`. From aae46c429a5b6e43d842257f767bf9bb7433efc9 Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Sun, 19 Jul 2026 17:05:28 +0100 Subject: [PATCH 14/23] Update --- cpp/dolfinx/mesh/Topology.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/cpp/dolfinx/mesh/Topology.cpp b/cpp/dolfinx/mesh/Topology.cpp index 1406cce5305..ebec9a58dd4 100644 --- a/cpp/dolfinx/mesh/Topology.cpp +++ b/cpp/dolfinx/mesh/Topology.cpp @@ -317,10 +317,11 @@ std::array, 2> vertex_ownership_groups( if (num_threads > 0) { boost::sort::block_indirect_sort(local_vertex_set.begin(), - local_vertex_set.end(), 10); + local_vertex_set.end(), num_threads); } else dolfinx::radix_sort(local_vertex_set); + auto [unique_end, range_end] = std::ranges::unique(local_vertex_set); local_vertex_set.erase(unique_end, range_end); } From 02f2214442a85d35f08c5c8c991031e8e2caf029 Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Mon, 20 Jul 2026 10:25:45 +0200 Subject: [PATCH 15/23] Updates --- cpp/dolfinx/io/VTKHDF.h | 2 +- cpp/dolfinx/mesh/Topology.cpp | 6 ++--- cpp/dolfinx/mesh/generation.h | 24 +++++++++--------- cpp/dolfinx/mesh/graphbuild.cpp | 38 +++++++++++++--------------- cpp/dolfinx/mesh/utils.h | 4 +-- cpp/test/mesh/branching_manifold.cpp | 4 +-- 6 files changed, 38 insertions(+), 40 deletions(-) diff --git a/cpp/dolfinx/io/VTKHDF.h b/cpp/dolfinx/io/VTKHDF.h index ad2f6c9ef4f..b1948dd2f79 100644 --- a/cpp/dolfinx/io/VTKHDF.h +++ b/cpp/dolfinx/io/VTKHDF.h @@ -468,6 +468,6 @@ mesh::Mesh read_mesh(MPI_Comm comm, const std::filesystem::path& filename, cells_local.end()); return mesh::create_mesh(comm, comm, cells_span, coordinate_elements, comm, points_pruned, {(std::size_t)x_shape[0], gdim}, part, - max_facet_to_cell_links, 0); + max_facet_to_cell_links, 1); } } // namespace dolfinx::io::VTKHDF diff --git a/cpp/dolfinx/mesh/Topology.cpp b/cpp/dolfinx/mesh/Topology.cpp index a16dbcc13a7..01cd6c671db 100644 --- a/cpp/dolfinx/mesh/Topology.cpp +++ b/cpp/dolfinx/mesh/Topology.cpp @@ -68,7 +68,7 @@ determine_sharing_ranks(MPI_Comm comm, std::span indices, int dest = dolfinx::MPI::index_owner(size, idx, global_range); dest_to_index.push_back({dest, static_cast(dest_to_index.size())}); } - if (num_threads > 0) + if (num_threads > 1) { boost::sort::block_indirect_sort(dest_to_index.begin(), dest_to_index.end(), num_threads); @@ -314,7 +314,7 @@ std::array, 2> vertex_ownership_groups( local_vertex_set.insert(local_vertex_set.end(), c.begin(), c.end()); { - if (num_threads > 0) + if (num_threads > 1) { boost::sort::block_indirect_sort(local_vertex_set.begin(), local_vertex_set.end(), num_threads); @@ -335,7 +335,7 @@ std::array, 2> vertex_ownership_groups( ghost_vertex_set.insert(ghost_vertex_set.end(), c.begin(), c.end()); { - if (num_threads > 0) + if (num_threads > 1) { boost::sort::block_indirect_sort(ghost_vertex_set.begin(), ghost_vertex_set.end(), num_threads); diff --git a/cpp/dolfinx/mesh/generation.h b/cpp/dolfinx/mesh/generation.h index f76328a3e9b..36ae80cba70 100644 --- a/cpp/dolfinx/mesh/generation.h +++ b/cpp/dolfinx/mesh/generation.h @@ -285,20 +285,20 @@ create_interval(MPI_Comm comm, std::int64_t n, std::array p, if (gdim == 1) { return create_mesh(comm, MPI_COMM_SELF, cells, element, MPI_COMM_SELF, - x1d, {npts, 1}, partitioner, 2, 0, reorder_fn); + x1d, {npts, 1}, partitioner, 2, 1, reorder_fn); } std::vector x(npts * gdim, T(0)); for (std::size_t i = 0; i < npts; i++) x[i * gdim] = x1d[i]; return create_mesh(comm, MPI_COMM_SELF, cells, element, MPI_COMM_SELF, x, {npts, static_cast(gdim)}, partitioner, 2, - 0, reorder_fn); + 1, reorder_fn); } else { return create_mesh(comm, MPI_COMM_NULL, {}, element, MPI_COMM_NULL, std::vector{}, {0, static_cast(gdim)}, - partitioner, 2, 0, reorder_fn); + partitioner, 2, 1, reorder_fn); } } @@ -422,7 +422,7 @@ Mesh build_tet(MPI_Comm comm, MPI_Comm subcomm, } return create_mesh(comm, subcomm, cells, element, subcomm, x, - {x.size() / 3, 3}, partitioner, 2, 0, reorder_fn); + {x.size() / 3, 3}, partitioner, 2, 1, reorder_fn); } template @@ -467,7 +467,7 @@ build_hex(MPI_Comm comm, MPI_Comm subcomm, std::array, 2> p, } return create_mesh(comm, subcomm, cells, element, subcomm, x, - {x.size() / 3, 3}, partitioner, 2, 0, reorder_fn); + {x.size() / 3, 3}, partitioner, 2, 1, reorder_fn); } template @@ -515,7 +515,7 @@ Mesh build_prism(MPI_Comm comm, MPI_Comm subcomm, } return create_mesh(comm, subcomm, cells, element, subcomm, x, - {x.size() / 3, 3}, partitioner, 2, 0, reorder_fn); + {x.size() / 3, 3}, partitioner, 2, 1, reorder_fn); } template @@ -672,7 +672,7 @@ Mesh build_tri(MPI_Comm comm, std::array, 2> p, if (gdim == 2) { return create_mesh(comm, MPI_COMM_SELF, cells, element, MPI_COMM_SELF, x, - {npts, 2}, partitioner, 2, 0, reorder_fn); + {npts, 2}, partitioner, 2, 1, reorder_fn); } std::vector xg(npts * gdim, T(0)); for (std::size_t i = 0; i < npts; i++) @@ -682,13 +682,13 @@ Mesh build_tri(MPI_Comm comm, std::array, 2> p, } return create_mesh(comm, MPI_COMM_SELF, cells, element, MPI_COMM_SELF, xg, {npts, static_cast(gdim)}, partitioner, 2, - 0, reorder_fn); + 1, reorder_fn); } else { return create_mesh(comm, MPI_COMM_NULL, {}, element, MPI_COMM_NULL, std::vector{}, {0, static_cast(gdim)}, - partitioner, 2, 0, reorder_fn); + partitioner, 2, 1, reorder_fn); } } @@ -738,7 +738,7 @@ Mesh build_quad(MPI_Comm comm, std::array, 2> p, if (gdim == 2) { return create_mesh(comm, MPI_COMM_SELF, cells, element, MPI_COMM_SELF, x, - {npts, 2}, partitioner, 2, 0, reorder_fn); + {npts, 2}, partitioner, 2, 1, reorder_fn); } std::vector xg(npts * gdim, T(0)); for (std::size_t i = 0; i < npts; i++) @@ -748,13 +748,13 @@ Mesh build_quad(MPI_Comm comm, std::array, 2> p, } return create_mesh(comm, MPI_COMM_SELF, cells, element, MPI_COMM_SELF, xg, {npts, static_cast(gdim)}, partitioner, 2, - 0, reorder_fn); + 1, reorder_fn); } else { return create_mesh(comm, MPI_COMM_NULL, {}, element, MPI_COMM_NULL, std::vector{}, {0, static_cast(gdim)}, - partitioner, 2, 0, reorder_fn); + partitioner, 2, 1, reorder_fn); } } } // namespace impl diff --git a/cpp/dolfinx/mesh/graphbuild.cpp b/cpp/dolfinx/mesh/graphbuild.cpp index 339be7e5fa4..f86aa61b354 100644 --- a/cpp/dolfinx/mesh/graphbuild.cpp +++ b/cpp/dolfinx/mesh/graphbuild.cpp @@ -639,27 +639,25 @@ mesh::build_local_dual_graph( graph::AdjacencyList cell_facets = mesh::get_entity_vertices(cell_type, tdim - 1); std::span _cells = cells[j]; - if (num_threads > 0) + std::vector threads; + for (int i = 1; i < num_threads; ++i) { - std::vector threads; - for (int i = 0; i < num_threads; ++i) - { - auto [c0, c1] = dolfinx::common::local_range( - i, _cells.size() / num_cell_vertices, num_threads); - threads.emplace_back( - build_facets_fn, shape1, num_cell_vertices, cell_offsets[j] + c0, - std::cref(cell_facets), - _cells.subspan(c0 * num_cell_vertices, - (c1 - c0) * num_cell_vertices), - std::span(facets.data() + (c0 * cell_facets.num_nodes()) * shape1, - ((c1 - c0) * cell_facets.num_nodes()) * shape1)); - } - } - else - { - build_facets_fn(shape1, num_cell_vertices, cell_offsets[j], - std::cref(cell_facets), _cells, facets); + auto [c0, c1] = common::local_range(i, _cells.size() / num_cell_vertices, + num_threads); + threads.emplace_back( + build_facets_fn, shape1, num_cell_vertices, cell_offsets[j] + c0, + std::cref(cell_facets), + _cells.subspan(c0 * num_cell_vertices, (c1 - c0) * num_cell_vertices), + std::span(facets.data() + (c0 * cell_facets.num_nodes()) * shape1, + ((c1 - c0) * cell_facets.num_nodes()) * shape1)); } + auto [c0, c1] = common::local_range(0, _cells.size() / num_cell_vertices, + num_threads); + build_facets_fn( + shape1, num_cell_vertices, cell_offsets[j] + c0, std::cref(cell_facets), + _cells.subspan(c0 * num_cell_vertices, (c1 - c0) * num_cell_vertices), + std::span(facets.data() + (c0 * cell_facets.num_nodes()) * shape1, + ((c1 - c0) * cell_facets.num_nodes()) * shape1)); } timer1.stop(); @@ -669,7 +667,7 @@ mesh::build_local_dual_graph( common::Timer timer3("Compute local part of mesh dual graph: 3"); std::vector perm; - if (num_threads > 0) + if (num_threads > 1) { perm.resize(facets.size() / shape1, 0); std::iota(perm.begin(), perm.end(), 0); diff --git a/cpp/dolfinx/mesh/utils.h b/cpp/dolfinx/mesh/utils.h index b72077531b7..4c3d1f8393a 100644 --- a/cpp/dolfinx/mesh/utils.h +++ b/cpp/dolfinx/mesh/utils.h @@ -1310,14 +1310,14 @@ create_mesh(MPI_Comm comm, std::span cells, if (dolfinx::MPI::size(comm) == 1) { return create_mesh(comm, comm, std::vector{cells}, std::vector{elements}, - comm, x, xshape, nullptr, max_facet_to_cell_links, 0); + comm, x, xshape, nullptr, max_facet_to_cell_links, 1); } else { return create_mesh( comm, comm, std::vector{cells}, std::vector{elements}, comm, x, xshape, create_cell_partitioner(ghost_mode, max_facet_to_cell_links), - max_facet_to_cell_links, 0); + max_facet_to_cell_links, 1); } } diff --git a/cpp/test/mesh/branching_manifold.cpp b/cpp/test/mesh/branching_manifold.cpp index 8d52fd47631..e30255760e7 100644 --- a/cpp/test/mesh/branching_manifold.cpp +++ b/cpp/test/mesh/branching_manifold.cpp @@ -125,7 +125,7 @@ TEST_CASE("dual_graph_branching") auto [dual_graph, unmatched_facets, max_vertices_per_facet, cell_data] = mesh::build_local_dual_graph(celltypes, {cells}, - max_facet_to_cell_links, 0); + max_facet_to_cell_links, 1); CHECK(dual_graph.num_nodes() == 4); @@ -184,7 +184,7 @@ TEST_CASE("dual_graph_self_dual") { auto [dual_graph, unmatched_facets, max_vertices_per_facet, cell_data] = mesh::build_local_dual_graph(celltypes, {cells}, - max_facet_to_cell_links, 0); + max_facet_to_cell_links, 1); CHECK(max_vertices_per_facet == 1); CHECK(dual_graph.num_nodes() == 3); From 810b27c6c67426fe0a44e9336b0322c2793165c5 Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Mon, 20 Jul 2026 11:08:38 +0200 Subject: [PATCH 16/23] Update --- cpp/dolfinx/mesh/graphbuild.h | 2 +- cpp/dolfinx/refinement/refine.h | 2 +- cpp/dolfinx/refinement/uniform.cpp | 2 +- cpp/test/mesh/branching_manifold.cpp | 8 +++----- cpp/test/mesh/distributed_mesh.cpp | 2 +- cpp/test/mesh/refinement/interval.cpp | 4 ++-- python/demo/demo_mixed-topology.py | 2 +- python/doc/source/release_notes.md | 19 +++++++++---------- python/dolfinx/io/gmsh.py | 2 +- python/dolfinx/io/utils.py | 2 +- python/test/unit/mesh/test_dual_graph.py | 2 +- python/test/unit/mesh/test_mixed_topology.py | 2 +- 12 files changed, 23 insertions(+), 26 deletions(-) diff --git a/cpp/dolfinx/mesh/graphbuild.h b/cpp/dolfinx/mesh/graphbuild.h index 2d75bc2d33b..77bf175275f 100644 --- a/cpp/dolfinx/mesh/graphbuild.h +++ b/cpp/dolfinx/mesh/graphbuild.h @@ -102,6 +102,6 @@ graph::AdjacencyList build_dual_graph(MPI_Comm comm, std::span celltypes, const std::vector>& cells, std::optional max_facet_to_cell_links, - int num_threads = 0); + int num_threads = 1); } // namespace dolfinx::mesh diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index dd9a3941e10..a606ce32eb7 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -152,7 +152,7 @@ refine(const mesh::Mesh& mesh, mesh::Mesh mesh1 = mesh::create_mesh( mesh.comm(), mesh.comm(), cell_adj.array(), mesh.geometry().cmaps().front(), mesh.comm(), new_vertex_coords, xshape, - std::get(partitioner), 2, 0); + std::get(partitioner), 2, 1s); // Report the number of refined cells const int D = topology->dim(); diff --git a/cpp/dolfinx/refinement/uniform.cpp b/cpp/dolfinx/refinement/uniform.cpp index 85f9ecd7f7b..073968d2edf 100644 --- a/cpp/dolfinx/refinement/uniform.cpp +++ b/cpp/dolfinx/refinement/uniform.cpp @@ -304,7 +304,7 @@ refinement::uniform_refine(const mesh::Mesh& mesh, geometry_cmaps.push_back(cm); mesh::Mesh new_mesh = mesh::create_mesh( mesh.comm(), mesh.comm(), topo_span, geometry_cmaps, mesh.comm(), new_x, - {new_x.size() / 3, 3}, partitioner, 2, 0); + {new_x.size() / 3, 3}, partitioner, 2, 1); return new_mesh; } diff --git a/cpp/test/mesh/branching_manifold.cpp b/cpp/test/mesh/branching_manifold.cpp index e30255760e7..e36d54b08e6 100644 --- a/cpp/test/mesh/branching_manifold.cpp +++ b/cpp/test/mesh/branching_manifold.cpp @@ -48,7 +48,7 @@ TEST_CASE("dual_graph_branching") { // default auto [dual_graph, unmatched_facets, max_vertices_per_facet, cell_data] - = mesh::build_local_dual_graph(celltypes, {cells}, 2, 0); + = mesh::build_local_dual_graph(celltypes, {cells}, 2, 1); CHECK(dual_graph.num_nodes() == 4); @@ -83,7 +83,7 @@ TEST_CASE("dual_graph_branching") // max_facet_to_cell_links = 3 // Note: additionally facet (2) is now considered unmatched auto [dual_graph, unmatched_facets, max_vertices_per_facet, cell_data] - = mesh::build_local_dual_graph(celltypes, {cells}, 3, 0); + = mesh::build_local_dual_graph(celltypes, {cells}, 3, 1); CHECK(dual_graph.num_nodes() == 4); @@ -254,12 +254,10 @@ TEST_CASE("dual_graph_branching_parallel") // Check local dual graphs. auto [dual_graph, unmatched_facets, max_vertices_per_facet, cell_data] - = mesh::build_local_dual_graph(celltypes, {cells}, 3, 0); + = mesh::build_local_dual_graph(celltypes, {cells}, 3, 1); CHECK(max_vertices_per_facet == 1); - CHECK(dual_graph.num_nodes() == 2); - CHECK(dual_graph.num_links(0) == 1); CHECK_THAT(dual_graph.links(0), Catch::Matchers::RangeEquals(std::array{1})); diff --git a/cpp/test/mesh/distributed_mesh.cpp b/cpp/test/mesh/distributed_mesh.cpp index 6eeb1659228..0a6734f89ca 100644 --- a/cpp/test/mesh/distributed_mesh.cpp +++ b/cpp/test/mesh/distributed_mesh.cpp @@ -139,7 +139,7 @@ void test_distributed_mesh(const mesh::CellPartitionFunction& partitioner) // Build mesh mesh::Mesh mesh = mesh::create_mesh(comm, subset_comm, cells, cmap, comm, x, - xshape, partitioner, 2, 0); + xshape, partitioner, 2, 1); auto t = mesh.topology(); int tdim = t->dim(); CHECK(t->index_map(tdim)->size_global() == 2 * N * N); diff --git a/cpp/test/mesh/refinement/interval.cpp b/cpp/test/mesh/refinement/interval.cpp index 59c49d87024..4e4b022055f 100644 --- a/cpp/test/mesh/refinement/interval.cpp +++ b/cpp/test/mesh/refinement/interval.cpp @@ -61,7 +61,7 @@ mesh::Mesh create_3_vertex_interval_mesh() return mesh::create_mesh( MPI_COMM_SELF, MPI_COMM_SELF, cells, element, MPI_COMM_SELF, x, {x.size() / 3, 3}, - mesh::create_cell_partitioner(mesh::GhostMode::none, 2), 2, 0); + mesh::create_cell_partitioner(mesh::GhostMode::none, 2), 2, 1); } TEMPLATE_TEST_CASE("Interval uniform refinement", @@ -152,7 +152,7 @@ TEMPLATE_TEST_CASE("Interval Refinement (parallel)", MPI_Comm commt = rank == 0 ? MPI_COMM_SELF : MPI_COMM_NULL; return mesh::create_mesh(MPI_COMM_WORLD, commt, cells, element, commt, x, - {x.size() / 3, 3}, partitioner, 2, 0); + {x.size() / 3, 3}, partitioner, 2, 1); }; mesh::Mesh mesh = create_mesh(); diff --git a/python/demo/demo_mixed-topology.py b/python/demo/demo_mixed-topology.py index 2b6a3efc652..09aa79053aa 100644 --- a/python/demo/demo_mixed-topology.py +++ b/python/demo/demo_mixed-topology.py @@ -112,7 +112,7 @@ part = create_cell_partitioner(GhostMode.none, 2) # type: ignore mesh = create_mesh( - MPI.COMM_WORLD, cells_np, [hexahedron._cpp_object, prism._cpp_object], geomx, part, 2, 0 + MPI.COMM_WORLD, cells_np, [hexahedron._cpp_object, prism._cpp_object], geomx, part, 2, 1 ) # - diff --git a/python/doc/source/release_notes.md b/python/doc/source/release_notes.md index acd47f0b5d1..385341943e0 100644 --- a/python/doc/source/release_notes.md +++ b/python/doc/source/release_notes.md @@ -28,7 +28,7 @@ The following constructors and classes have been added: solver options, and {py:class}`entity_maps`. ### Built-in matrix support - + **Authors**: [Chris Richardson](https://github.com/chrisrichardson) - Adds {py:meth}`A.transpose()`, {py:meth}`A.mult(x, y, transpose=True)` @@ -69,13 +69,12 @@ and [Garth N. Wells](https://github.com/garth-wells) For a long time, DOLFINx has been exclusively using MPI for the distribution of computational load. However, with the computational landscape evolving to more and more heterogeneous systems, -the need for additional parallelisation methods are required. +the need for additional parallelisation methods are required. In this release, we introduce initial threading support using [std::jthread`](https://en.cppreference.com/cpp/thread/jthread) in the following methods: - {py:meth}`dolfinx.mesh.Topology.create_entities` - {py:func}`dolfinx.geometry.compute_distances_gjk` which both take an optional argument `num_threads` which specifies how many CPU threads should be used. -If set to 0, threads are not spawned. ### New and improved demos @@ -153,8 +152,8 @@ reading {py:class}`entity tags` on mixed-topology grids t **Authors**: [Jørgen S. Dokken](https://github.com/jorgensd) A set of users have had issues with non-affine geometries, in particular higher order grids, -and the tolerance and maximum number of iterations in the -{py:meth}`dolfinx.fem.CoordinateElement.pull_back`, {py:meth}`dolfinx.fem.Function.interpolate_nonmatching` +and the tolerance and maximum number of iterations in the +{py:meth}`dolfinx.fem.CoordinateElement.pull_back`, {py:meth}`dolfinx.fem.Function.interpolate_nonmatching` and {py:meth}`dolfinx.fem.Function.eval` yielding errors such as: ```bash RuntimeError: Newton method failed to converge for non-affine geometry @@ -175,7 +174,7 @@ Furthermore, {py:func}`dolfinx.fem.apply_lifting` and {py:func}`dolfinx.fem.asse ### Meshes -**Authors**: [Jørgen S. Dokken](https://github.com/jorgensd) and [Jack Hale](https://github.com/jhale) +**Authors**: [Jørgen S. Dokken](https://github.com/jorgensd) and [Jack Hale](https://github.com/jhale) - New function {py:func}`dolfinx.mesh.create_point_mesh` to create a point cloud mesh with no points shared between the different processes. Useful for reading in point measures or outputting data. @@ -219,7 +218,7 @@ In addition to the changes below, the ever-lasting quest of improving performanc [Jack Hale](https://github.com/jhale) and [Garth N. Wells](https://github.com/garth-wells) Mapping data between {py:class}`PETSc.Vec` and {py:class}`dolfinx.fem.Function`s is now -trivial for blocked problems by using {py:func}`dolfinx.fem.petsc.assign`. +trivial for blocked problems by using {py:func}`dolfinx.fem.petsc.assign`. Both solvers and assembly routines interfacing with PETSc has received a drastic make-over to improve usability and maintenance, both for developers and end-users. @@ -387,7 +386,7 @@ The GMSH interface to DOLFINx has received a major upgrade. #### VTKHDF5 -**Authors**: [Chris Richardson](https://github.com/chrisrichardson) and [Jørgen S. Dokken](https://github.com/jorgensd) +**Authors**: [Chris Richardson](https://github.com/chrisrichardson) and [Jørgen S. Dokken](https://github.com/jorgensd) As Kitware has stated that [VTKHDF](https://www.kitware.com/vtk-hdf-reader/) is the future format they want to support, we have started the transition to this format. @@ -409,12 +408,12 @@ from another (unique) function space. **Author**: [Massimiliano Leoni](https://github.com/mleoni-pf) and [Paul T. Kühner](https://github.com/schnellerhase) - When using {py:meth}`dolfinx.io.XDMFFIle.read_meshtags` one can now specify the attribute name, if the grid has -multiple tags assigned to it. +multiple tags assigned to it. - Flushing data to file is now possible with {py:meth}`dolfinx.io.XDMFFile.flush`. This is useful when wanting to visualize long-running jobs in Paraview. #### Remove Fides backend -As we unfortunately haven't seen an expanding set of features for the +As we unfortunately haven't seen an expanding set of features for the [Fides Reader](https://fides.readthedocs.io/en/latest/paraview/paraview.html) in Paraview, we have decided to remove it from DOLFINx. diff --git a/python/dolfinx/io/gmsh.py b/python/dolfinx/io/gmsh.py index 75c2ac55dc6..b946ae1bee9 100644 --- a/python/dolfinx/io/gmsh.py +++ b/python/dolfinx/io/gmsh.py @@ -441,7 +441,7 @@ def model_to_mesh( x[:, :gdim].astype(dtype).copy(), partitioner, max_facet_to_cell_links, - 0, + 1, ) mesh = Mesh(cpp_mesh, None) diff --git a/python/dolfinx/io/utils.py b/python/dolfinx/io/utils.py index 82381feb9d7..8476d317611 100644 --- a/python/dolfinx/io/utils.py +++ b/python/dolfinx/io/utils.py @@ -273,7 +273,7 @@ def read_mesh( x, create_cell_partitioner(ghost_mode, max_facet_to_cell_links), # type: ignore max_facet_to_cell_links, - 0, + 1, ) msh.name = name domain = ufl.Mesh(basix_el) diff --git a/python/test/unit/mesh/test_dual_graph.py b/python/test/unit/mesh/test_dual_graph.py index 42f10013245..646d4dd3fd5 100644 --- a/python/test/unit/mesh/test_dual_graph.py +++ b/python/test/unit/mesh/test_dual_graph.py @@ -26,7 +26,7 @@ def test_dgrsph_1d(): # Circular chain of interval cells cells = [[n0, n0 + 1], [n0 + 1, n0 + 2], [n0 + 2, x]] w = mesh.build_dual_graph( - MPI.COMM_WORLD, mesh.CellType.interval, to_adj(cells, np.int64)._cpp_object, 2, 0 + MPI.COMM_WORLD, mesh.CellType.interval, to_adj(cells, np.int64)._cpp_object, 2, 1 ) assert w.num_nodes == 3 for i in range(w.num_nodes): diff --git a/python/test/unit/mesh/test_mixed_topology.py b/python/test/unit/mesh/test_mixed_topology.py index 6d07d178932..488e9000ec2 100644 --- a/python/test/unit/mesh/test_mixed_topology.py +++ b/python/test/unit/mesh/test_mixed_topology.py @@ -344,7 +344,7 @@ def test_locate_entities(dtype): comm = MPI.COMM_WORLD max_cells_per_facet = 2 mesh = create_mesh( - comm, cells, [hexahedron._cpp_object, prism._cpp_object], geom, part, max_cells_per_facet, 0 + comm, cells, [hexahedron._cpp_object, prism._cpp_object], geom, part, max_cells_per_facet, 1 ) fdim = mesh.topology.dim - 1 From 74feb1611f14c8ad650a5c7753988d8494106c01 Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Mon, 20 Jul 2026 11:20:09 +0200 Subject: [PATCH 17/23] Fix --- 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 a606ce32eb7..f41e27614a1 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -152,7 +152,7 @@ refine(const mesh::Mesh& mesh, mesh::Mesh mesh1 = mesh::create_mesh( mesh.comm(), mesh.comm(), cell_adj.array(), mesh.geometry().cmaps().front(), mesh.comm(), new_vertex_coords, xshape, - std::get(partitioner), 2, 1s); + std::get(partitioner), 2, 1); // Report the number of refined cells const int D = topology->dim(); From 8059b9853b07d8431fcf31de4c3ade95df8b2dc5 Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Mon, 20 Jul 2026 11:22:14 +0200 Subject: [PATCH 18/23] Doc fixes --- cpp/dolfinx/mesh/graphbuild.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/cpp/dolfinx/mesh/graphbuild.h b/cpp/dolfinx/mesh/graphbuild.h index 77bf175275f..3ebbb7db1a3 100644 --- a/cpp/dolfinx/mesh/graphbuild.h +++ b/cpp/dolfinx/mesh/graphbuild.h @@ -33,8 +33,8 @@ enum class CellType : std::int8_t; /// Passing std::nullopt (no upper bound) corresponds. /// to `max_facet_to_cell_links`=∞, i.e. every facet is considered /// unmatched. -/// @param[in] num_threads Number of threads to use. Use 0 to not launch -/// threads. +/// @param[in] num_threads Number of threads to use. Must be greater +/// than 0. /// /// @return /// 1. Local dual graph @@ -86,8 +86,8 @@ build_local_dual_graph(std::span celltypes, /// manifold meshes. Passing std::nullopt (no upper bound) corresponds /// to `max_facet_to_cell_links`=∞, i.e. every facet is considered /// unmatched. -/// @param[in] num_threads Number of threads to use. Use 0 to not launch -/// threads. +/// @param[in] num_threads Number of threads to use. Must be greater +/// than 0. /// /// @return The dual graph. /// From 75f51dd9f1b574114fe00ed0add7085e637ec3c7 Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Mon, 20 Jul 2026 11:23:45 +0200 Subject: [PATCH 19/23] Python fix --- python/dolfinx/mesh.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index 945e63d3ac7..d3d44477497 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -792,7 +792,7 @@ def create_mesh( x: npt.NDArray[np.floating], partitioner: Callable | None = None, max_facet_to_cell_links: int = 2, - num_threads: int = 0, + num_threads: int = 1, ) -> Mesh: """Create a mesh from topology and geometry arrays. @@ -808,8 +808,8 @@ def create_mesh( cells across MPI ranks. max_facet_to_cell_links: Maximum number of cells a facet can be connected to. - num_threads: Number of threads to use tp build mesh. Set to 0 to - not spawn threads. + num_threads: Number of threads to use tp build mesh. Must be + greater than 0. Note: If required, the coordinates ``x`` will be cast to the same From 1b8875afc64daa8440cb2a991af67217eebf3a71 Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Mon, 20 Jul 2026 11:29:28 +0200 Subject: [PATCH 20/23] Fixes --- python/test/unit/conftest.py | 2 +- python/test/unit/io/test_vtkhdf.py | 2 +- python/test/unit/mesh/test_create_mixed_mesh.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/python/test/unit/conftest.py b/python/test/unit/conftest.py index 30510e3472c..ed1c2cbd21f 100644 --- a/python/test/unit/conftest.py +++ b/python/test/unit/conftest.py @@ -120,5 +120,5 @@ def mixed_topology_mesh(): geomx, part, max_cells_per_facet, - 0, + num_threads=1, ) diff --git a/python/test/unit/io/test_vtkhdf.py b/python/test/unit/io/test_vtkhdf.py index 8143a581d54..d21d4190d12 100644 --- a/python/test/unit/io/test_vtkhdf.py +++ b/python/test/unit/io/test_vtkhdf.py @@ -91,7 +91,7 @@ def test_read_write_higher_order(): geom, part, max_cells_per_facet, - 0, + num_threads=1, ) py_mesh = Mesh(mesh, None) diff --git a/python/test/unit/mesh/test_create_mixed_mesh.py b/python/test/unit/mesh/test_create_mixed_mesh.py index b2038f5ad3a..e8a018bdb00 100644 --- a/python/test/unit/mesh/test_create_mixed_mesh.py +++ b/python/test/unit/mesh/test_create_mixed_mesh.py @@ -91,7 +91,7 @@ def test_create_mixed_mesh(dtype): geomx, part, max_cells_per_facet, - 0, + num_threads=0, ) entity_types = mesh.topology.entity_types[3] From 89fd9abf20d98c6d6212de1066e3c840ac3c6d17 Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Mon, 20 Jul 2026 11:35:27 +0200 Subject: [PATCH 21/23] Another fix --- python/test/unit/mesh/test_create_mixed_mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/test/unit/mesh/test_create_mixed_mesh.py b/python/test/unit/mesh/test_create_mixed_mesh.py index e8a018bdb00..1fe66d7cac6 100644 --- a/python/test/unit/mesh/test_create_mixed_mesh.py +++ b/python/test/unit/mesh/test_create_mixed_mesh.py @@ -91,7 +91,7 @@ def test_create_mixed_mesh(dtype): geomx, part, max_cells_per_facet, - num_threads=0, + num_threads=1, ) entity_types = mesh.topology.entity_types[3] From 4218a1f9cc6530887fe06f505a85586dd041d7a8 Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Mon, 20 Jul 2026 12:28:31 +0200 Subject: [PATCH 22/23] Fixes --- cpp/dolfinx/mesh/Topology.cpp | 6 +++--- cpp/dolfinx/mesh/graphbuild.cpp | 21 +++++++++++++++------ cpp/dolfinx/mesh/utils.h | 6 +++--- 3 files changed, 21 insertions(+), 12 deletions(-) diff --git a/cpp/dolfinx/mesh/Topology.cpp b/cpp/dolfinx/mesh/Topology.cpp index 01cd6c671db..faa5fc84c5f 100644 --- a/cpp/dolfinx/mesh/Topology.cpp +++ b/cpp/dolfinx/mesh/Topology.cpp @@ -147,7 +147,7 @@ determine_sharing_ranks(MPI_Comm comm, std::span indices, // Build {global index, pos, src} list std::vector> indices_list; { - common::Timer timer("Topology: XXXXX"); + common::Timer timer("Topology: build and sort transposed index list"); for (std::size_t p = 0; p < recv_disp0.size() - 1; ++p) for (std::int32_t i = recv_disp0[p]; i < recv_disp0[p + 1]; ++i) indices_list.push_back({recv_buffer0[i], i, int(p)}); @@ -360,8 +360,8 @@ std::array, 2> vertex_ownership_groups( std::back_inserter(unowned_vertices)); #ifndef NDEBUG - // Sanity check - // No vertices in unowned should also be in boundary... + // Sanity check: no vertices in unowned should also be in boundary. + // Test in DEBUG mode only because of cost. std::vector unowned_vertices_in_error; std::ranges::set_intersection(unowned_vertices, boundary_vertices, std::back_inserter(unowned_vertices_in_error)); diff --git a/cpp/dolfinx/mesh/graphbuild.cpp b/cpp/dolfinx/mesh/graphbuild.cpp index f86aa61b354..c9ed2964aa5 100644 --- a/cpp/dolfinx/mesh/graphbuild.cpp +++ b/cpp/dolfinx/mesh/graphbuild.cpp @@ -19,6 +19,7 @@ #include #include #include +#include #include #include #include @@ -535,6 +536,9 @@ mesh::build_local_dual_graph( const std::vector>& cells, std::optional max_facet_to_cell_links, int num_threads) { + if (num_threads < 1) + throw std::runtime_error("num_threads must be >= 1."); + spdlog::info("Build local part of mesh dual graph"); common::Timer timer("Compute local part of mesh dual graph"); @@ -632,6 +636,7 @@ mesh::build_local_dual_graph( const int shape1 = max_vertices_per_facet + 1; std::vector facets(facet_count * shape1); + std::size_t facet_offset = 0; for (std::size_t j = 0; j < cells.size(); ++j) { CellType cell_type = celltypes[j]; @@ -639,25 +644,29 @@ mesh::build_local_dual_graph( graph::AdjacencyList cell_facets = mesh::get_entity_vertices(cell_type, tdim - 1); std::span _cells = cells[j]; + std::size_t num_cells_j = _cells.size() / num_cell_vertices; std::vector threads; for (int i = 1; i < num_threads; ++i) { - auto [c0, c1] = common::local_range(i, _cells.size() / num_cell_vertices, - num_threads); + auto [c0, c1] = common::local_range(i, num_cells_j, num_threads); threads.emplace_back( build_facets_fn, shape1, num_cell_vertices, cell_offsets[j] + c0, std::cref(cell_facets), _cells.subspan(c0 * num_cell_vertices, (c1 - c0) * num_cell_vertices), - std::span(facets.data() + (c0 * cell_facets.num_nodes()) * shape1, + std::span(facets.data() + + (facet_offset + c0 * cell_facets.num_nodes()) + * shape1, ((c1 - c0) * cell_facets.num_nodes()) * shape1)); } - auto [c0, c1] = common::local_range(0, _cells.size() / num_cell_vertices, - num_threads); + auto [c0, c1] = common::local_range(0, num_cells_j, num_threads); build_facets_fn( shape1, num_cell_vertices, cell_offsets[j] + c0, std::cref(cell_facets), _cells.subspan(c0 * num_cell_vertices, (c1 - c0) * num_cell_vertices), - std::span(facets.data() + (c0 * cell_facets.num_nodes()) * shape1, + std::span(facets.data() + + (facet_offset + c0 * cell_facets.num_nodes()) * shape1, ((c1 - c0) * cell_facets.num_nodes()) * shape1)); + + facet_offset += num_cells_j * cell_facets.num_nodes(); } timer1.stop(); diff --git a/cpp/dolfinx/mesh/utils.h b/cpp/dolfinx/mesh/utils.h index 4c3d1f8393a..1fe7abcca58 100644 --- a/cpp/dolfinx/mesh/utils.h +++ b/cpp/dolfinx/mesh/utils.h @@ -1042,7 +1042,7 @@ compute_incident_entities(const Topology& topology, /// @param[in] max_facet_to_cell_links Bound on the number of cells a /// facet can be connected to. /// @param[in] num_threads Number threads to use in mesh construction. -/// It 0, no threads are launched. +/// Must be >= 1. /// @param[in] reorder_fn Function that reorders (locally) cells that /// are owned by this process. /// @return A mesh distributed on the communicator `comm`. @@ -1259,8 +1259,8 @@ Mesh> create_mesh( /// rank for each cell. If not callable, cells are not redistributed. /// @param[in] max_facet_to_cell_links Bound on the number of cells a /// facet can be connected to. -/// @param[in] num_threads Number threads to use in mesh construction. If 0, no -/// threads are launched. +/// @param[in] num_threads Number threads to use in mesh construction. +/// Must be >= 1. /// @param[in] reorder_fn Function that reorders (locally) cells that /// are owned by this process. /// @return A mesh distributed on the communicator `comm`. From 36a613d959a38277bf5cd335249864cb974eaa0f Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Tue, 21 Jul 2026 08:52:49 +0200 Subject: [PATCH 23/23] Update python/dolfinx/mesh.py Co-authored-by: Chris Richardson --- 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 d3d44477497..2c7ba98e9f3 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -808,7 +808,7 @@ def create_mesh( cells across MPI ranks. max_facet_to_cell_links: Maximum number of cells a facet can be connected to. - num_threads: Number of threads to use tp build mesh. Must be + num_threads: Number of threads to use to build mesh. Must be greater than 0. Note: