From e1ed42ad55951be975cdf5481ba86f0ee0289e61 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Wed, 21 Aug 2024 13:13:07 +0200 Subject: [PATCH 01/12] Remove refine option `distribute` and `ghost_mode` in favor of an explicit partitioner --- cpp/dolfinx/refinement/refine.h | 67 ++-- cpp/dolfinx/refinement/utils.h | 45 --- cpp/test/mesh/refinement/interval.cpp | 7 +- cpp/test/mesh/refinement/rectangle.cpp | 6 +- cpp/test/transfer/transfer_matrix.cpp | 315 ++++++++++++++++++ python/dolfinx/mesh.py | 12 +- python/dolfinx/wrappers/refinement.cpp | 78 ++++- python/test/unit/refinement/test_interval.py | 12 +- .../test/unit/refinement/test_refinement.py | 20 +- 9 files changed, 434 insertions(+), 128 deletions(-) create mode 100644 cpp/test/transfer/transfer_matrix.cpp diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index 32f18d89827..50d2181207c 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -17,46 +17,40 @@ #include #include +#include "dolfinx/graph/AdjacencyList.h" +#include "dolfinx/mesh/Mesh.h" +#include "dolfinx/mesh/Topology.h" +#include "dolfinx/mesh/cell_types.h" +#include "dolfinx/mesh/utils.h" + +#include "interval.h" +#include "plaza.h" + namespace dolfinx::refinement { -namespace impl -{ -/// @brief create the refined mesh by optionally redistributing it -template -mesh::Mesh -create_refined_mesh(const mesh::Mesh& mesh, - const graph::AdjacencyList& cell_adj, - std::span new_vertex_coords, - std::array xshape, bool redistribute, - mesh::GhostMode ghost_mode) + +// TODO: move to cpp? +inline graph::AdjacencyList maintain_coarse_partitioner( + MPI_Comm comm, int, const std::vector& cell_types, + const std::vector>& cell_topology) { - if (dolfinx::MPI::size(mesh.comm()) == 1) - { - // No parallel construction necessary, i.e. redistribute also has no - // effect - return mesh::create_mesh(mesh.comm(), cell_adj.array(), - mesh.geometry().cmap(), new_vertex_coords, xshape, - ghost_mode); - } - else - { - // Let partition handle the parallel construction of the mesh - return partition(mesh, cell_adj, new_vertex_coords, xshape, redistribute, - ghost_mode); - } + std::int32_t mpi_rank = MPI::rank(comm); + std::int32_t num_cell_vertices = mesh::num_cell_vertices(cell_types.front()); + std::int32_t num_cells = cell_topology.front().size() / num_cell_vertices; + std::vector destinations(num_cells, mpi_rank); + std::vector dest_offsets(num_cells + 1); + std::iota(dest_offsets.begin(), dest_offsets.end(), 0); + return graph::AdjacencyList(std::move(destinations), std::move(dest_offsets)); } -} // namespace impl /// @brief Refine with markers, optionally redistributing, and /// optionally calculating the parent-child relationships. /// /// @param[in] mesh Input mesh to be refined -/// @param[in] edges Optional indices of the edges that should be split -/// by this refinement, if optional is not set, a uniform refinement is -/// performed, same behavior as passing a list of all indices. -/// @param[in] redistribute Flag to call the Mesh Partitioner to -/// redistribute after refinement. -/// @param[in] ghost_mode Ghost mode of the refined mesh. +/// @param[in] edges Optional indices of the edges that should be split by this +/// refinement, if optional is not set, a uniform refinement is performend, same +/// behavior as passing a list of all indices. +/// @param[in] partitioner partitioner to be used for the refined mesh /// @param[in] option Control the computation of parent facets, parent /// cells. If an option is unselected, an empty list is returned. /// @return New Mesh and optional parent cell index, parent facet @@ -65,8 +59,8 @@ template std::tuple, std::optional>, std::optional>> refine(const mesh::Mesh& mesh, - std::optional> edges, bool redistribute, - mesh::GhostMode ghost_mode = mesh::GhostMode::shared_facet, + std::optional> edges, + mesh::CellPartitionFunction partitioner = maintain_coarse_partitioner, Option option = Option::none) { auto topology = mesh.topology(); @@ -80,9 +74,10 @@ refine(const mesh::Mesh& mesh, = oned ? interval::compute_refinement_data(mesh, edges, option) : plaza::compute_refinement_data(mesh, edges, option); - mesh::Mesh refined_mesh = impl::create_refined_mesh( - mesh, cell_adj, std::span(new_vertex_coords), xshape, - redistribute, ghost_mode); + mesh::Mesh refined_mesh + = mesh::create_mesh(mesh.comm(), mesh.comm(), std::move(cell_adj).array(), + mesh.geometry().cmap(), mesh.comm(), + std::move(new_vertex_coords), xshape, partitioner); // Report the number of refined cellse const int D = topology->dim(); diff --git a/cpp/dolfinx/refinement/utils.h b/cpp/dolfinx/refinement/utils.h index 275b40d7cc2..1c6af077551 100644 --- a/cpp/dolfinx/refinement/utils.h +++ b/cpp/dolfinx/refinement/utils.h @@ -264,51 +264,6 @@ create_new_vertices(MPI_Comm comm, xshape}; } -/// Use vertex and topology data to partition new mesh across -/// processes -/// @param[in] old_mesh -/// @param[in] cell_topology Topology of cells, (vertex indices) -/// @param[in] new_coords New coordinates, row-major storage -/// @param[in] xshape The shape of `new_coords` -/// @param[in] redistribute Call graph partitioner if true -/// @param[in] ghost_mode None or shared_facet -/// @return New mesh -template -mesh::Mesh partition(const mesh::Mesh& old_mesh, - const graph::AdjacencyList& cell_topology, - std::span new_coords, - std::array xshape, bool redistribute, - mesh::GhostMode ghost_mode) -{ - if (redistribute) - { - return mesh::create_mesh(old_mesh.comm(), cell_topology.array(), - old_mesh.geometry().cmap(), new_coords, xshape, - ghost_mode); - } - else - { - auto partitioner - = [](MPI_Comm comm, int, const std::vector& cell_types, - const std::vector>& cell_topology) - { - std::int32_t mpi_rank = MPI::rank(comm); - std::int32_t num_cell_vertices - = mesh::num_cell_vertices(cell_types.front()); - std::int32_t num_cells = cell_topology.front().size() / num_cell_vertices; - std::vector destinations(num_cells, mpi_rank); - std::vector dest_offsets(num_cells + 1); - std::iota(dest_offsets.begin(), dest_offsets.end(), 0); - return graph::AdjacencyList(std::move(destinations), - std::move(dest_offsets)); - }; - - return mesh::create_mesh(old_mesh.comm(), old_mesh.comm(), - cell_topology.array(), old_mesh.geometry().cmap(), - old_mesh.comm(), new_coords, xshape, partitioner); - } -} - /// @brief Given an index map, add "n" extra indices at the end of local range /// /// @note The returned global indices (local and ghosts) are adjust for the diff --git a/cpp/test/mesh/refinement/interval.cpp b/cpp/test/mesh/refinement/interval.cpp index 763418ac4b3..74b3f434627 100644 --- a/cpp/test/mesh/refinement/interval.cpp +++ b/cpp/test/mesh/refinement/interval.cpp @@ -73,7 +73,7 @@ TEMPLATE_TEST_CASE("Interval uniform refinement", // TODO: parent_facet auto [refined_mesh, parent_edge, parent_facet] = refinement::refine( - mesh, std::nullopt, false, mesh::GhostMode::shared_facet, + mesh, std::nullopt, &refinement::maintain_coarse_partitioner, refinement::Option::parent_cell); std::vector expected_x = { @@ -114,7 +114,8 @@ TEMPLATE_TEST_CASE("Interval adaptive refinement", std::vector edges{1}; // TODO: parent_facet auto [refined_mesh, parent_edge, parent_facet] = refinement::refine( - mesh, std::span(edges), false, mesh::GhostMode::shared_facet, + mesh, std::span(edges), + mesh::create_cell_partitioner(mesh::GhostMode::shared_facet), refinement::Option::parent_cell); std::vector expected_x = { @@ -192,7 +193,7 @@ TEMPLATE_TEST_CASE("Interval Refinement (parallel)", // TODO: parent_facet auto [refined_mesh, parent_edges, parent_facet] = refinement::refine( - mesh, std::nullopt, false, mesh::GhostMode::shared_facet, + mesh, std::nullopt, &refinement::maintain_coarse_partitioner, refinement::Option::parent_cell); T rank_d = static_cast(rank); diff --git a/cpp/test/mesh/refinement/rectangle.cpp b/cpp/test/mesh/refinement/rectangle.cpp index 954879abb1d..ff2a67979df 100644 --- a/cpp/test/mesh/refinement/rectangle.cpp +++ b/cpp/test/mesh/refinement/rectangle.cpp @@ -99,9 +99,9 @@ plotter.show() // plaza requires the edges to be pre initialized! mesh.topology()->create_entities(1); - auto [mesh_fine, parent_cell, parent_facet] - = refinement::refine(mesh, std::nullopt, false, mesh::GhostMode::none, - refinement::Option::parent_cell_and_facet); + auto [mesh_fine, parent_cell, parent_facet] = refinement::refine( + mesh, std::nullopt, mesh::create_cell_partitioner(mesh::GhostMode::none), + refinement::Option::parent_cell_and_facet); // vertex layout: // 8---7---5 diff --git a/cpp/test/transfer/transfer_matrix.cpp b/cpp/test/transfer/transfer_matrix.cpp new file mode 100644 index 00000000000..9747e776fcb --- /dev/null +++ b/cpp/test/transfer/transfer_matrix.cpp @@ -0,0 +1,315 @@ +// Copyright (C) 2024 Paul Kühner +// +// This file is part of DOLFINX (https://www.fenicsproject.org) +// +// SPDX-License-Identifier: LGPL-3.0-or-later + +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../mesh/util.h" + +using namespace dolfinx; +using namespace Catch::Matchers; + +template +constexpr auto weight = [](std::int32_t d) -> T { return d == 0 ? 1. : .5; }; + +TEMPLATE_TEST_CASE("Transfer Matrix 1D", "[transfer_matrix]", + double) // TODO: float +{ + using T = TestType; + + auto mesh_coarse + = std::make_shared>(dolfinx::mesh::create_interval( + MPI_COMM_SELF, 2, {0.0, 1.0}, mesh::GhostMode::none)); + + auto [mesh_fine, parent_cell, parent_facet] = refinement::refine( + *mesh_coarse, std::nullopt, &refinement::maintain_coarse_partitioner, mesh::GhostMode::none, + refinement::Option::parent_cell); + + auto element = basix::create_element( + basix::element::family::P, basix::cell::type::interval, 1, + basix::element::lagrange_variant::unset, + basix::element::dpc_variant::unset, false); + + auto V_coarse = std::make_shared>( + fem::create_functionspace(mesh_coarse, element, {})); + auto V_fine + = std::make_shared>(fem::create_functionspace( + std::make_shared>(mesh_fine), element, {})); + + mesh_fine.topology()->create_connectivity(1, 0); + mesh_fine.topology()->create_connectivity(0, 1); + + std::vector from_to_map{0, 2, 4}; // TODO: general computation! + + la::MatrixCSR transfer_matrix = transfer::create_transfer_matrix( + *V_coarse, *V_fine, from_to_map, weight); + + std::vector expected{1.0, .5, 0, 0, 0, 0, .5, 1, .5, 0, 0, 0, 0, .5, 1}; + CHECK_THAT(transfer_matrix.to_dense(), + RangeEquals(expected, [](auto a, auto b) + { return std::abs(a - b) <= EPS; })); +} + +TEMPLATE_TEST_CASE("Transfer Matrix 1D (parallel)", "[transfer_matrix]", + double) // TODO: float +{ + using T = TestType; + + int comm_size = dolfinx::MPI::size(MPI_COMM_WORLD); + int rank = dolfinx::MPI::rank(MPI_COMM_WORLD); + + // TODO: see https://github.com/FEniCS/dolfinx/issues/3358 + // auto part = mesh::create_cell_partitioner(mesh::GhostMode::shared_facet); + mesh::CellPartitionFunction part + = [&](MPI_Comm /* comm */, int /* nparts */, + const std::vector& /* cell_types */, + const std::vector>& /* cells */) + { + std::vector> data; + if (comm_size == 1) + data = {{0}, {0}, {0}, {0}}; + else if (comm_size == 2) + data = {{0}, {0}, {0}, {0}, {0, 1}, {1, 0}, {1}, {1}, {1}}; + else if (comm_size == 3) + data = {{1}, {1}, {1}, {1}, {1, 2}, {2, 1}, {2}, + {2}, {2, 0}, {0, 2}, {0}, {0}, {0}, {0}}; + else + FAIL("Test only supports <= 3 processes"); + return graph::AdjacencyList(std::move(data)); + }; + + auto mesh_coarse = std::make_shared>( + mesh::create_interval(MPI_COMM_WORLD, 5 * comm_size - 1, {0., 1.}, + mesh::GhostMode::shared_facet, part)); + + auto [mesh_fine, parent_cell, parent_facet] = refinement::refine( + *mesh_coarse, std::nullopt, part, mesh::GhostMode::shared_facet, + refinement::Option::none); + + auto element = basix::create_element( + basix::element::family::P, basix::cell::type::interval, 1, + basix::element::lagrange_variant::unset, + basix::element::dpc_variant::unset, false); + + auto V_coarse = std::make_shared>( + fem::create_functionspace(mesh_coarse, element, {})); + auto V_fine + = std::make_shared>(fem::create_functionspace( + std::make_shared>(mesh_fine), element, {})); + + mesh_fine.topology()->create_connectivity(1, 0); + mesh_fine.topology()->create_connectivity(0, 1); + + std::vector from_to_map; + switch (comm_size) + { + case 1: + from_to_map = {0, 2, 4, 6, 8}; + break; + case 2: + from_to_map = {0, 2, 4, 6, 8, 10, 12, 14, 16, 18}; + break; + case 3: + from_to_map = {0, 2, 4, 6, 8, 9, 11, 13, 15, 17, 19, 28, 26, 24, 22}; + break; + } + + la::MatrixCSR transfer_matrix = transfer::create_transfer_matrix( + *V_coarse, *V_fine, from_to_map, weight); + + std::array, 3> expected; + if (comm_size == 1) + { + // clang-format off + expected[0] = { + /* row_0 */ 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + /* row_1 */ 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, + /* row_2 */ 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, + /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, + /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0 + }; + // clang-format on + } + else if (comm_size == 2) + { + // clang-format off + expected[0] = { + /* row_0 */ 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + /* row_1 */ 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + /* row_2 */ 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + /* row_6 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + }; + expected[1] = { + /* row_0 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, + /* row_1 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, + /* row_2 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, + /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, + /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + }; + // clang-format on + } + else if (comm_size == 3) + { + // clang-format off + expected[0] = { + /* row_0 */ 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + /* row_1 */ 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + /* row_2 */ 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + /* row_6 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + }; + expected[1] = { + /* row_0 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + /* row_1 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + /* row_2 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, + /* row_6 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + }; + expected[2] = { + /* row_0 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, + /* row_1 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, + /* row_2 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, + /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, + /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + /* row_6 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + }; + // clang-format on + } + else + { + CHECK(false); + } + + CHECK_THAT(transfer_matrix.to_dense(), + RangeEquals(expected[rank], [](auto a, auto b) + { return std::abs(a - b) <= EPS; })); +} + +TEMPLATE_TEST_CASE("Transfer Matrix 2D", "[transfer_matrix]", double) +{ + using T = TestType; + + auto mesh_coarse + = std::make_shared>(dolfinx::mesh::create_rectangle( + MPI_COMM_SELF, {{{0, 0}, {1, 1}}}, {1, 1}, mesh::CellType::triangle)); + mesh_coarse->topology()->create_entities(1); + + auto [mesh_fine, parent_cell, parent_facet] = refinement::refine( + *mesh_coarse, std::nullopt, &refinement::maintain_coarse_partitioner, mesh::GhostMode::none, + refinement::Option::parent_cell); + + auto element = basix::create_element( + basix::element::family::P, basix::cell::type::triangle, 1, + basix::element::lagrange_variant::unset, + basix::element::dpc_variant::unset, false); + + auto V_coarse = std::make_shared>( + fem::create_functionspace(mesh_coarse, element, {})); + auto V_fine + = std::make_shared>(fem::create_functionspace( + std::make_shared>(mesh_fine), element, {})); + + mesh_fine.topology()->create_connectivity(1, 0); + mesh_fine.topology()->create_connectivity(0, 1); + + std::vector from_to_map{4, 1, 5, + 8}; // TODO: general computation! + + la::MatrixCSR transfer_matrix = transfer::create_transfer_matrix( + *V_coarse, *V_fine, from_to_map, weight); + + std::vector expected{0.5, 0.0, 0.5, 0.0, 1.0, 0.0, 0.5, 0.0, 0.0, + 0.5, 1.0, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.5, 0.0, 0.0, 0.5, 0.0, 1.0, 0.0, 0.5, 0.0, + 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5, 1.0}; + CHECK_THAT(transfer_matrix.to_dense(), + RangeEquals(expected, [](auto a, auto b) + { return std::abs(a - b) <= EPS; })); +} + +TEMPLATE_TEST_CASE("Transfer Matrix 3D", "[transfer_matrix]", double) +{ + using T = TestType; + + auto mesh_coarse = std::make_shared>( + dolfinx::mesh::create_box(MPI_COMM_SELF, {{{0, 0, 0}, {1, 1, 1}}}, + {1, 1, 1}, mesh::CellType::tetrahedron)); + mesh_coarse->topology()->create_entities(1); + + auto [mesh_fine, parent_cell, parent_facet] = refinement::refine( + *mesh_coarse, std::nullopt, &refinement::maintain_coarse_partitioner, mesh::GhostMode::none, + refinement::Option::parent_cell); + + auto element = basix::create_element( + basix::element::family::P, basix::cell::type::tetrahedron, 1, + basix::element::lagrange_variant::unset, + basix::element::dpc_variant::unset, false); + + auto V_coarse = std::make_shared>( + fem::create_functionspace(mesh_coarse, element, {})); + auto V_fine + = std::make_shared>(fem::create_functionspace( + std::make_shared>(mesh_fine), element, {})); + + mesh_fine.topology()->create_connectivity(1, 0); + mesh_fine.topology()->create_connectivity(0, 1); + + std::vector from_to_map{ + 0, 6, 15, 25, 17, 9, 11, 22}; // TODO: general computation! + + la::MatrixCSR transfer_matrix = transfer::create_transfer_matrix( + *V_coarse, *V_fine, from_to_map, weight); + + std::vector expected{ + 1.0, 0.5, 0.5, 0.5, 0.5, 0.5, 0.0, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5, + 0.5, 0.5, 0.0, 1.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.5, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.5, 0.5, 1.0, 0.0, 0.0, 0.0, 0.5, 0.0, + 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.5, + 0.5, 1.0, 0.5, 0.0, 0.5, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.5, 0.0, 0.5, 0.0, 0.5, 1.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.5, 0.0, 0.0, + 0.0, 0.5, 0.5, 0.0, 0.0, 0.5, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0, 0.0, 0.5, 0.0, + 0.0, 0.0, 0.0, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, + 0.0, 0.5, 0.0, 0.0, 0.5, 0.5, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, + 0.0, 0.0, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5, 0.5, + 0.5, 1.0, 0.0, 0.0, 0.0, 0.5}; + CHECK_THAT(transfer_matrix.to_dense(), + RangeEquals(expected, [](auto a, auto b) + { return std::abs(a - b) <= EPS; })); +} diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index b9cefce09b3..46fec0bed34 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -370,8 +370,7 @@ def transfer_meshtag( def refine( mesh: Mesh, edges: typing.Optional[np.ndarray] = None, - redistribute: bool = True, - ghost_mode: GhostMode = GhostMode.shared_facet, + partitioner: typing.Optional[typing.Callable] = None, option: RefinementOption = RefinementOption.none, ) -> tuple[Mesh, npt.NDArray[np.int32], npt.NDArray[np.int8]]: """Refine a mesh. @@ -380,17 +379,16 @@ def refine( mesh: Mesh from which to create the refined mesh. edges: Indices of edges to split during refinement. If ``None``, mesh refinement is uniform. - redistribute: - Refined mesh is re-partitioned if ``True``. - ghost_mode: Ghost mode to use for the refined mesh. + partitioner: + partitioner to use for the refined mesh, If ``None`` no redistribution is performed, + i.e. previous local mesh is equally parallelized now with new vertices.s option: Controls whether parent cells and/or parent facets are computed. - Returns: Refined mesh, (optional) parent cells, (optional) parent facets """ mesh1, parent_cell, parent_facet = _cpp.refinement.refine( - mesh._cpp_object, edges, redistribute, ghost_mode, option + mesh._cpp_object, edges, partitioner, option ) return Mesh(mesh1, mesh._ufl_domain), parent_cell, parent_facet diff --git a/python/dolfinx/wrappers/refinement.cpp b/python/dolfinx/wrappers/refinement.cpp index 37ca3937992..7921ffb22e5 100644 --- a/python/dolfinx/wrappers/refinement.cpp +++ b/python/dolfinx/wrappers/refinement.cpp @@ -4,9 +4,20 @@ // // SPDX-License-Identifier: LGPL-3.0-or-later -#include "array.h" #include #include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + #include #include #include @@ -14,30 +25,67 @@ #include #include #include -#include -#include -#include -#include -#include -#include -#include -#include + +#include "MPICommWrapper.h" +#include "array.h" namespace nb = nanobind; namespace dolfinx_wrappers { +namespace +{ + +using PythonCellPartitionFunction + = std::function( + dolfinx_wrappers::MPICommWrapper, int, + const std::vector&, + std::vector>)>; + +using CppCellPartitionFunction + = std::function( + MPI_Comm, int, const std::vector& q, + const std::vector>&)>; + +/// Wrap a Python graph partitioning function as a C++ function +CppCellPartitionFunction +create_cell_partitioner_cpp(const PythonCellPartitionFunction& p) +{ + if (p) + { + return [p](MPI_Comm comm, int n, + const std::vector& cell_types, + const std::vector>& cells) + { + std::vector> cells_nb; + std::ranges::transform( + cells, std::back_inserter(cells_nb), + [](auto c) + { + return nb::ndarray( + c.data(), {c.size()}, nb::handle()); + }); + + return p(dolfinx_wrappers::MPICommWrapper(comm), n, cell_types, cells_nb); + }; + } + else + return nullptr; +} +} // namespace + template void export_refinement_with_variable_mesh_type(nb::module_& m) { + m.def( "refine", [](const dolfinx::mesh::Mesh& mesh, std::optional< nb::ndarray, nb::c_contig>> edges, - bool redistribute, dolfinx::mesh::GhostMode ghost_mode, + std::optional partitioner, dolfinx::refinement::Option option) { std::optional> cpp_edges(std::nullopt); @@ -47,8 +95,12 @@ void export_refinement_with_variable_mesh_type(nb::module_& m) std::span(edges.value().data(), edges.value().size())); } + CppCellPartitionFunction cpp_partitioner + = partitioner.has_value() + ? create_cell_partitioner_cpp(partitioner.value()) + : dolfinx::refinement::maintain_coarse_partitioner; auto [mesh1, cell, facet] = dolfinx::refinement::refine( - mesh, cpp_edges, redistribute, ghost_mode, option); + mesh, cpp_edges, cpp_partitioner, option); std::optional> python_cell( std::nullopt); @@ -63,8 +115,8 @@ void export_refinement_with_variable_mesh_type(nb::module_& m) return std::tuple{std::move(mesh1), std::move(python_cell), std::move(python_facet)}; }, - nb::arg("mesh"), nb::arg("edges") = nb::none(), nb::arg("redistribute"), - nb::arg("ghost_mode"), nb::arg("option")); + nb::arg("mesh"), nb::arg("edges") = nb::none(), + nb::arg("partitioner") = nb::none(), nb::arg("option")); } void refinement(nb::module_& m) diff --git a/python/test/unit/refinement/test_interval.py b/python/test/unit/refinement/test_interval.py index 28779644840..27040fb0b90 100644 --- a/python/test/unit/refinement/test_interval.py +++ b/python/test/unit/refinement/test_interval.py @@ -20,16 +20,14 @@ "ghost_mode_refined", [mesh.GhostMode.none, mesh.GhostMode.shared_vertex, mesh.GhostMode.shared_facet], ) -@pytest.mark.parametrize("redistribute", [True, False]) @pytest.mark.parametrize("option", [mesh.RefinementOption.none, mesh.RefinementOption.parent_cell]) -def test_refine_interval(n, ghost_mode, redistribute, ghost_mode_refined, option): +def test_refine_interval(n, ghost_mode, ghost_mode_refined, option): msh = mesh.create_interval(MPI.COMM_WORLD, n, [0, 1], ghost_mode=ghost_mode) msh_refined, edges, vertices = mesh.refine( msh, - redistribute=redistribute, - ghost_mode=ghost_mode_refined, option=option, ) + # TODO: add create_cell_partitioner(ghost_mode) when works # vertex count assert msh_refined.topology.index_map(0).size_global == 2 * n + 1 @@ -52,16 +50,14 @@ def test_refine_interval(n, ghost_mode, redistribute, ghost_mode_refined, option "ghost_mode_refined", [mesh.GhostMode.none, mesh.GhostMode.shared_vertex, mesh.GhostMode.shared_facet], ) -@pytest.mark.parametrize("redistribute", [True, False]) -def test_refine_interval_adaptive(n, ghost_mode, redistribute, ghost_mode_refined): +def test_refine_interval_adaptive(n, ghost_mode, ghost_mode_refined): msh = mesh.create_interval(MPI.COMM_WORLD, n, [0, 1], ghost_mode=ghost_mode) msh_refined, edges, vertices = mesh.refine( msh, np.arange(10, dtype=np.int32), - redistribute=redistribute, - ghost_mode=ghost_mode_refined, option=mesh.RefinementOption.parent_cell, ) + # TODO: add create_cell_partitioner(ghost_mode) when works # vertex count assert msh_refined.topology.index_map(0).size_global == n + 1 + 10 * MPI.COMM_WORLD.size diff --git a/python/test/unit/refinement/test_refinement.py b/python/test/unit/refinement/test_refinement.py index 440057ce265..eb5188d5bf1 100644 --- a/python/test/unit/refinement/test_refinement.py +++ b/python/test/unit/refinement/test_refinement.py @@ -32,7 +32,7 @@ def test_refine_create_unit_square(): """Refine mesh of unit square.""" mesh = create_unit_square(MPI.COMM_WORLD, 5, 7, ghost_mode=GhostMode.none) mesh.topology.create_entities(1) - mesh_refined, _, _ = refine(mesh, redistribute=False) + mesh_refined, _, _ = refine(mesh) assert mesh_refined.topology.index_map(0).size_global == 165 assert mesh_refined.topology.index_map(2).size_global == 280 @@ -46,7 +46,7 @@ def test_refine_create_unit_cube(ghost_mode, redistribute): """Refine mesh of unit cube.""" mesh = create_unit_cube(MPI.COMM_WORLD, 5, 7, 9, ghost_mode=ghost_mode) mesh.topology.create_entities(1) - mesh, _, _ = refine(mesh, redistribute=redistribute) + mesh, _, _ = refine(mesh) # , partitioner=create_cell_partitioner(ghost_mode)) assert mesh.topology.index_map(0).size_global == 3135 assert mesh.topology.index_map(3).size_global == 15120 @@ -58,7 +58,7 @@ def test_refine_create_form(): """Check that forms can be assembled on refined mesh""" mesh = create_unit_cube(MPI.COMM_WORLD, 3, 3, 3) mesh.topology.create_entities(1) - mesh, _, _ = refine(mesh, redistribute=True) + mesh, _, _ = refine(mesh) # TODO: recover redistribute=True behavior V = functionspace(mesh, ("Lagrange", 1)) @@ -83,7 +83,7 @@ def left_corner_edge(x, tol=1e-7): if MPI.COMM_WORLD.size == 0: assert edges == 1 - mesh2, _, _ = refine(mesh, edges, redistribute=False) + mesh2, _, _ = refine(mesh, edges) assert mesh.topology.index_map(2).size_global + 3 == mesh2.topology.index_map(2).size_global @@ -105,7 +105,7 @@ def left_side(x, tol=1e-16): edges = compute_incident_entities(msh.topology, cells, 2, 1) if MPI.COMM_WORLD.size == 0: assert edges.__len__() == Nx // 2 * (2 * Ny + 1) + (Nx // 2 + 1) * Ny - mesh2, _, _ = refine(msh, edges, redistribute=True) + mesh2, _, _ = refine(msh, edges) # TODO: redistribute=True num_cells_global = mesh2.topology.index_map(2).size_global actual_cells = 3 * (Nx * Ny) + 3 * Ny + 2 * Nx * Ny @@ -116,13 +116,10 @@ def left_side(x, tol=1e-16): @pytest.mark.parametrize( "refine_plaza_wrapper", [ - lambda mesh: refine( - mesh, redistribute=False, option=RefinementOption.parent_cell_and_facet - ), + lambda mesh: refine(mesh, option=RefinementOption.parent_cell_and_facet), lambda mesh: refine( mesh, np.arange(mesh.topology.index_map(1).size_local), - redistribute=False, option=RefinementOption.parent_cell_and_facet, ), ], @@ -178,13 +175,10 @@ def test_refine_facet_meshtag(tdim, refine_plaza_wrapper): @pytest.mark.parametrize( "refine_plaza_wrapper", [ - lambda mesh: refine( - mesh, redistribute=False, option=RefinementOption.parent_cell_and_facet - ), + lambda mesh: refine(mesh, option=RefinementOption.parent_cell_and_facet), lambda mesh: refine( mesh, np.arange(mesh.topology.index_map(1).size_local), - redistribute=False, option=RefinementOption.parent_cell_and_facet, ), ], From 32fca7cb63bbf64732a24541066f61ca254f20df Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Tue, 17 Sep 2024 22:43:48 +0200 Subject: [PATCH 02/12] Remove wrongly added transfer_matrix -> later PR --- cpp/test/transfer/transfer_matrix.cpp | 315 -------------------------- 1 file changed, 315 deletions(-) delete mode 100644 cpp/test/transfer/transfer_matrix.cpp diff --git a/cpp/test/transfer/transfer_matrix.cpp b/cpp/test/transfer/transfer_matrix.cpp deleted file mode 100644 index 9747e776fcb..00000000000 --- a/cpp/test/transfer/transfer_matrix.cpp +++ /dev/null @@ -1,315 +0,0 @@ -// Copyright (C) 2024 Paul Kühner -// -// This file is part of DOLFINX (https://www.fenicsproject.org) -// -// SPDX-License-Identifier: LGPL-3.0-or-later - -#include -#include -#include -#include - -#include - -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include "../mesh/util.h" - -using namespace dolfinx; -using namespace Catch::Matchers; - -template -constexpr auto weight = [](std::int32_t d) -> T { return d == 0 ? 1. : .5; }; - -TEMPLATE_TEST_CASE("Transfer Matrix 1D", "[transfer_matrix]", - double) // TODO: float -{ - using T = TestType; - - auto mesh_coarse - = std::make_shared>(dolfinx::mesh::create_interval( - MPI_COMM_SELF, 2, {0.0, 1.0}, mesh::GhostMode::none)); - - auto [mesh_fine, parent_cell, parent_facet] = refinement::refine( - *mesh_coarse, std::nullopt, &refinement::maintain_coarse_partitioner, mesh::GhostMode::none, - refinement::Option::parent_cell); - - auto element = basix::create_element( - basix::element::family::P, basix::cell::type::interval, 1, - basix::element::lagrange_variant::unset, - basix::element::dpc_variant::unset, false); - - auto V_coarse = std::make_shared>( - fem::create_functionspace(mesh_coarse, element, {})); - auto V_fine - = std::make_shared>(fem::create_functionspace( - std::make_shared>(mesh_fine), element, {})); - - mesh_fine.topology()->create_connectivity(1, 0); - mesh_fine.topology()->create_connectivity(0, 1); - - std::vector from_to_map{0, 2, 4}; // TODO: general computation! - - la::MatrixCSR transfer_matrix = transfer::create_transfer_matrix( - *V_coarse, *V_fine, from_to_map, weight); - - std::vector expected{1.0, .5, 0, 0, 0, 0, .5, 1, .5, 0, 0, 0, 0, .5, 1}; - CHECK_THAT(transfer_matrix.to_dense(), - RangeEquals(expected, [](auto a, auto b) - { return std::abs(a - b) <= EPS; })); -} - -TEMPLATE_TEST_CASE("Transfer Matrix 1D (parallel)", "[transfer_matrix]", - double) // TODO: float -{ - using T = TestType; - - int comm_size = dolfinx::MPI::size(MPI_COMM_WORLD); - int rank = dolfinx::MPI::rank(MPI_COMM_WORLD); - - // TODO: see https://github.com/FEniCS/dolfinx/issues/3358 - // auto part = mesh::create_cell_partitioner(mesh::GhostMode::shared_facet); - mesh::CellPartitionFunction part - = [&](MPI_Comm /* comm */, int /* nparts */, - const std::vector& /* cell_types */, - const std::vector>& /* cells */) - { - std::vector> data; - if (comm_size == 1) - data = {{0}, {0}, {0}, {0}}; - else if (comm_size == 2) - data = {{0}, {0}, {0}, {0}, {0, 1}, {1, 0}, {1}, {1}, {1}}; - else if (comm_size == 3) - data = {{1}, {1}, {1}, {1}, {1, 2}, {2, 1}, {2}, - {2}, {2, 0}, {0, 2}, {0}, {0}, {0}, {0}}; - else - FAIL("Test only supports <= 3 processes"); - return graph::AdjacencyList(std::move(data)); - }; - - auto mesh_coarse = std::make_shared>( - mesh::create_interval(MPI_COMM_WORLD, 5 * comm_size - 1, {0., 1.}, - mesh::GhostMode::shared_facet, part)); - - auto [mesh_fine, parent_cell, parent_facet] = refinement::refine( - *mesh_coarse, std::nullopt, part, mesh::GhostMode::shared_facet, - refinement::Option::none); - - auto element = basix::create_element( - basix::element::family::P, basix::cell::type::interval, 1, - basix::element::lagrange_variant::unset, - basix::element::dpc_variant::unset, false); - - auto V_coarse = std::make_shared>( - fem::create_functionspace(mesh_coarse, element, {})); - auto V_fine - = std::make_shared>(fem::create_functionspace( - std::make_shared>(mesh_fine), element, {})); - - mesh_fine.topology()->create_connectivity(1, 0); - mesh_fine.topology()->create_connectivity(0, 1); - - std::vector from_to_map; - switch (comm_size) - { - case 1: - from_to_map = {0, 2, 4, 6, 8}; - break; - case 2: - from_to_map = {0, 2, 4, 6, 8, 10, 12, 14, 16, 18}; - break; - case 3: - from_to_map = {0, 2, 4, 6, 8, 9, 11, 13, 15, 17, 19, 28, 26, 24, 22}; - break; - } - - la::MatrixCSR transfer_matrix = transfer::create_transfer_matrix( - *V_coarse, *V_fine, from_to_map, weight); - - std::array, 3> expected; - if (comm_size == 1) - { - // clang-format off - expected[0] = { - /* row_0 */ 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - /* row_1 */ 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, - /* row_2 */ 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, - /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, - /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0 - }; - // clang-format on - } - else if (comm_size == 2) - { - // clang-format off - expected[0] = { - /* row_0 */ 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - /* row_1 */ 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - /* row_2 */ 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - /* row_6 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 - }; - expected[1] = { - /* row_0 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, - /* row_1 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, - /* row_2 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, - /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, - /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 - }; - // clang-format on - } - else if (comm_size == 3) - { - // clang-format off - expected[0] = { - /* row_0 */ 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - /* row_1 */ 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - /* row_2 */ 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - /* row_6 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - }; - expected[1] = { - /* row_0 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - /* row_1 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - /* row_2 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, - /* row_6 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - }; - expected[2] = { - /* row_0 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, - /* row_1 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, - /* row_2 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, - /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, - /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - /* row_6 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - }; - // clang-format on - } - else - { - CHECK(false); - } - - CHECK_THAT(transfer_matrix.to_dense(), - RangeEquals(expected[rank], [](auto a, auto b) - { return std::abs(a - b) <= EPS; })); -} - -TEMPLATE_TEST_CASE("Transfer Matrix 2D", "[transfer_matrix]", double) -{ - using T = TestType; - - auto mesh_coarse - = std::make_shared>(dolfinx::mesh::create_rectangle( - MPI_COMM_SELF, {{{0, 0}, {1, 1}}}, {1, 1}, mesh::CellType::triangle)); - mesh_coarse->topology()->create_entities(1); - - auto [mesh_fine, parent_cell, parent_facet] = refinement::refine( - *mesh_coarse, std::nullopt, &refinement::maintain_coarse_partitioner, mesh::GhostMode::none, - refinement::Option::parent_cell); - - auto element = basix::create_element( - basix::element::family::P, basix::cell::type::triangle, 1, - basix::element::lagrange_variant::unset, - basix::element::dpc_variant::unset, false); - - auto V_coarse = std::make_shared>( - fem::create_functionspace(mesh_coarse, element, {})); - auto V_fine - = std::make_shared>(fem::create_functionspace( - std::make_shared>(mesh_fine), element, {})); - - mesh_fine.topology()->create_connectivity(1, 0); - mesh_fine.topology()->create_connectivity(0, 1); - - std::vector from_to_map{4, 1, 5, - 8}; // TODO: general computation! - - la::MatrixCSR transfer_matrix = transfer::create_transfer_matrix( - *V_coarse, *V_fine, from_to_map, weight); - - std::vector expected{0.5, 0.0, 0.5, 0.0, 1.0, 0.0, 0.5, 0.0, 0.0, - 0.5, 1.0, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, - 0.5, 0.0, 0.0, 0.5, 0.0, 1.0, 0.0, 0.5, 0.0, - 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5, 1.0}; - CHECK_THAT(transfer_matrix.to_dense(), - RangeEquals(expected, [](auto a, auto b) - { return std::abs(a - b) <= EPS; })); -} - -TEMPLATE_TEST_CASE("Transfer Matrix 3D", "[transfer_matrix]", double) -{ - using T = TestType; - - auto mesh_coarse = std::make_shared>( - dolfinx::mesh::create_box(MPI_COMM_SELF, {{{0, 0, 0}, {1, 1, 1}}}, - {1, 1, 1}, mesh::CellType::tetrahedron)); - mesh_coarse->topology()->create_entities(1); - - auto [mesh_fine, parent_cell, parent_facet] = refinement::refine( - *mesh_coarse, std::nullopt, &refinement::maintain_coarse_partitioner, mesh::GhostMode::none, - refinement::Option::parent_cell); - - auto element = basix::create_element( - basix::element::family::P, basix::cell::type::tetrahedron, 1, - basix::element::lagrange_variant::unset, - basix::element::dpc_variant::unset, false); - - auto V_coarse = std::make_shared>( - fem::create_functionspace(mesh_coarse, element, {})); - auto V_fine - = std::make_shared>(fem::create_functionspace( - std::make_shared>(mesh_fine), element, {})); - - mesh_fine.topology()->create_connectivity(1, 0); - mesh_fine.topology()->create_connectivity(0, 1); - - std::vector from_to_map{ - 0, 6, 15, 25, 17, 9, 11, 22}; // TODO: general computation! - - la::MatrixCSR transfer_matrix = transfer::create_transfer_matrix( - *V_coarse, *V_fine, from_to_map, weight); - - std::vector expected{ - 1.0, 0.5, 0.5, 0.5, 0.5, 0.5, 0.0, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5, - 0.5, 0.5, 0.0, 1.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.5, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.5, 0.5, 1.0, 0.0, 0.0, 0.0, 0.5, 0.0, - 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.5, - 0.5, 1.0, 0.5, 0.0, 0.5, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - 0.5, 0.0, 0.5, 0.0, 0.5, 1.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.5, 0.0, 0.0, - 0.0, 0.5, 0.5, 0.0, 0.0, 0.5, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0, 0.0, 0.5, 0.0, - 0.0, 0.0, 0.0, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, - 0.0, 0.5, 0.0, 0.0, 0.5, 0.5, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, - 0.0, 0.0, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5, 0.5, - 0.5, 1.0, 0.0, 0.0, 0.0, 0.5}; - CHECK_THAT(transfer_matrix.to_dense(), - RangeEquals(expected, [](auto a, auto b) - { return std::abs(a - b) <= EPS; })); -} From 5bd8ba5387c894bd993c2de87747baca345394cf Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Tue, 17 Sep 2024 22:46:02 +0200 Subject: [PATCH 03/12] Add author name --- 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 50d2181207c..12d5ccce196 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -1,4 +1,4 @@ -// Copyright (C) 2010-2023 Garth N. Wells +// Copyright (C) 2010-2024 Garth N. Wells and Paul T. Kühner // // This file is part of DOLFINx (https://www.fenicsproject.org) // From 124630aea480a1823c43e1dcc51371bba1c5efec Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Wed, 18 Sep 2024 11:49:28 +0100 Subject: [PATCH 04/12] Doc fix --- cpp/dolfinx/refinement/refine.h | 29 ++++++++++------------------- 1 file changed, 10 insertions(+), 19 deletions(-) diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index 12d5ccce196..4e72c35837c 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -6,6 +6,7 @@ #pragma once +#include "dolfinx/graph/AdjacencyList.h" #include "dolfinx/mesh/Mesh.h" #include "dolfinx/mesh/Topology.h" #include "dolfinx/mesh/cell_types.h" @@ -17,25 +18,16 @@ #include #include -#include "dolfinx/graph/AdjacencyList.h" -#include "dolfinx/mesh/Mesh.h" -#include "dolfinx/mesh/Topology.h" -#include "dolfinx/mesh/cell_types.h" -#include "dolfinx/mesh/utils.h" - -#include "interval.h" -#include "plaza.h" - namespace dolfinx::refinement { -// TODO: move to cpp? +/// TODO: Document function inline graph::AdjacencyList maintain_coarse_partitioner( MPI_Comm comm, int, const std::vector& cell_types, const std::vector>& cell_topology) { - std::int32_t mpi_rank = MPI::rank(comm); - std::int32_t num_cell_vertices = mesh::num_cell_vertices(cell_types.front()); + int mpi_rank = MPI::rank(comm); + int num_cell_vertices = mesh::num_cell_vertices(cell_types.front()); std::int32_t num_cells = cell_topology.front().size() / num_cell_vertices; std::vector destinations(num_cells, mpi_rank); std::vector dest_offsets(num_cells + 1); @@ -46,10 +38,10 @@ inline graph::AdjacencyList maintain_coarse_partitioner( /// @brief Refine with markers, optionally redistributing, and /// optionally calculating the parent-child relationships. /// -/// @param[in] mesh Input mesh to be refined -/// @param[in] edges Optional indices of the edges that should be split by this -/// refinement, if optional is not set, a uniform refinement is performend, same -/// behavior as passing a list of all indices. +/// @param[in] mesh Input mesh to be refined. +/// @param[in] edges Optional indices of the edges that should be split +/// by this refinement, if optional is not set, a uniform refinement is +/// performed, same behavior as passing a list of all indices. /// @param[in] partitioner partitioner to be used for the refined mesh /// @param[in] option Control the computation of parent facets, parent /// cells. If an option is unselected, an empty list is returned. @@ -65,10 +57,9 @@ refine(const mesh::Mesh& mesh, { auto topology = mesh.topology(); assert(topology); - - mesh::CellType cell_t = topology->cell_type(); - if (!mesh::is_simplex(cell_t)) + if (!mesh::is_simplex(topology->cell_type())) throw std::runtime_error("Refinement only defined for simplices"); + bool oned = topology->cell_type() == mesh::CellType::interval; auto [cell_adj, new_vertex_coords, xshape, parent_cell, parent_facet] = oned ? interval::compute_refinement_data(mesh, edges, option) From 8a706621c4c064a6b2d66b74e1bab206ac8c8764 Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Wed, 18 Sep 2024 11:49:58 +0100 Subject: [PATCH 05/12] Fix typi --- 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 4e72c35837c..1473a4db4e0 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -70,7 +70,7 @@ refine(const mesh::Mesh& mesh, mesh.geometry().cmap(), mesh.comm(), std::move(new_vertex_coords), xshape, partitioner); - // Report the number of refined cellse + // Report the number of refined cells const int D = topology->dim(); const std::int64_t n0 = topology->index_map(D)->size_global(); const std::int64_t n1 = refined_mesh.topology()->index_map(D)->size_global(); From 973f16dca8705a87c23cbbbc5d4e0401520ad139 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Thu, 19 Sep 2024 07:44:53 +0200 Subject: [PATCH 06/12] Move maintain_coarse_partitioner to refine.cpp --- cpp/dolfinx/refinement/CMakeLists.txt | 1 + cpp/dolfinx/refinement/refine.h | 13 ++----------- 2 files changed, 3 insertions(+), 11 deletions(-) diff --git a/cpp/dolfinx/refinement/CMakeLists.txt b/cpp/dolfinx/refinement/CMakeLists.txt index 9560e45c823..b11546361db 100644 --- a/cpp/dolfinx/refinement/CMakeLists.txt +++ b/cpp/dolfinx/refinement/CMakeLists.txt @@ -11,4 +11,5 @@ set(HEADERS_refinement target_sources( dolfinx PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/plaza.cpp ${CMAKE_CURRENT_SOURCE_DIR}/utils.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/refine.cpp ) diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index 1473a4db4e0..e55245d4750 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -22,18 +22,9 @@ namespace dolfinx::refinement { /// TODO: Document function -inline graph::AdjacencyList maintain_coarse_partitioner( +graph::AdjacencyList maintain_coarse_partitioner( MPI_Comm comm, int, const std::vector& cell_types, - const std::vector>& cell_topology) -{ - int mpi_rank = MPI::rank(comm); - int num_cell_vertices = mesh::num_cell_vertices(cell_types.front()); - std::int32_t num_cells = cell_topology.front().size() / num_cell_vertices; - std::vector destinations(num_cells, mpi_rank); - std::vector dest_offsets(num_cells + 1); - std::iota(dest_offsets.begin(), dest_offsets.end(), 0); - return graph::AdjacencyList(std::move(destinations), std::move(dest_offsets)); -} + const std::vector>& cell_topology); /// @brief Refine with markers, optionally redistributing, and /// optionally calculating the parent-child relationships. From face987aba59853b04a65bf65459da9057a1f76c Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Thu, 19 Sep 2024 07:47:41 +0200 Subject: [PATCH 07/12] Works better with the actual file --- cpp/dolfinx/refinement/refine.cpp | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) create mode 100644 cpp/dolfinx/refinement/refine.cpp diff --git a/cpp/dolfinx/refinement/refine.cpp b/cpp/dolfinx/refinement/refine.cpp new file mode 100644 index 00000000000..961f5cb98fc --- /dev/null +++ b/cpp/dolfinx/refinement/refine.cpp @@ -0,0 +1,25 @@ +// Copyright (C) 2024 Paul T. Kühner +// +// This file is part of DOLFINx (https://www.fenicsproject.org) +// +// SPDX-License-Identifier: LGPL-3.0-or-later + +#include "refine.h" + +namespace dolfinx::refinement +{ + +graph::AdjacencyList maintain_coarse_partitioner( + MPI_Comm comm, int, const std::vector& cell_types, + const std::vector>& cell_topology) +{ + int mpi_rank = MPI::rank(comm); + int num_cell_vertices = mesh::num_cell_vertices(cell_types.front()); + std::int32_t num_cells = cell_topology.front().size() / num_cell_vertices; + std::vector destinations(num_cells, mpi_rank); + std::vector dest_offsets(num_cells + 1); + std::iota(dest_offsets.begin(), dest_offsets.end(), 0); + return graph::AdjacencyList(std::move(destinations), std::move(dest_offsets)); +} + +} // namespace dolfinx::refinement From 965c28bdcb59166da33a36a2cb51efffe044c39c Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Thu, 19 Sep 2024 07:53:48 +0200 Subject: [PATCH 08/12] Add doc string to partitioner --- cpp/dolfinx/refinement/refine.h | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index e55245d4750..4d1420e7460 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -21,7 +21,15 @@ namespace dolfinx::refinement { -/// TODO: Document function +/// @brief Partitioner that maintains the regional destribution for a refined +/// mesh. Meaning, process local data is just maintained and no redistribution +/// happens. +/// +/// @param[in] comm MPI Communicator +/// @param[in] nparts Number of partitions (input has no effect) +/// @param[in] cell_types Cell types in the mesh +/// @param[in] cells Lists of cells of each cell type. +/// @return Destination ranks for each cell on this process graph::AdjacencyList maintain_coarse_partitioner( MPI_Comm comm, int, const std::vector& cell_types, const std::vector>& cell_topology); From ba75f6afd86dd897bcfcd6a7aa595d0be789ee0b Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Thu, 19 Sep 2024 08:01:37 +0200 Subject: [PATCH 09/12] Make doxygen happy --- 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 4d1420e7460..dddfa908104 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -31,7 +31,7 @@ namespace dolfinx::refinement /// @param[in] cells Lists of cells of each cell type. /// @return Destination ranks for each cell on this process graph::AdjacencyList maintain_coarse_partitioner( - MPI_Comm comm, int, const std::vector& cell_types, + MPI_Comm comm, int nparts, const std::vector& cell_types, const std::vector>& cell_topology); /// @brief Refine with markers, optionally redistributing, and From 1f001cae4e81949fe6400a0e97957d8b542ace4f Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Thu, 19 Sep 2024 08:13:11 +0200 Subject: [PATCH 10/12] Switch to general partitioner naming --- cpp/dolfinx/refinement/refine.cpp | 4 ++-- cpp/dolfinx/refinement/refine.h | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/cpp/dolfinx/refinement/refine.cpp b/cpp/dolfinx/refinement/refine.cpp index 961f5cb98fc..2b5ce85d3c3 100644 --- a/cpp/dolfinx/refinement/refine.cpp +++ b/cpp/dolfinx/refinement/refine.cpp @@ -11,11 +11,11 @@ namespace dolfinx::refinement graph::AdjacencyList maintain_coarse_partitioner( MPI_Comm comm, int, const std::vector& cell_types, - const std::vector>& cell_topology) + const std::vector>& cells) { int mpi_rank = MPI::rank(comm); int num_cell_vertices = mesh::num_cell_vertices(cell_types.front()); - std::int32_t num_cells = cell_topology.front().size() / num_cell_vertices; + std::int32_t num_cells = cells.front().size() / num_cell_vertices; std::vector destinations(num_cells, mpi_rank); std::vector dest_offsets(num_cells + 1); std::iota(dest_offsets.begin(), dest_offsets.end(), 0); diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index dddfa908104..c23bc0a3c30 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -32,7 +32,7 @@ namespace dolfinx::refinement /// @return Destination ranks for each cell on this process graph::AdjacencyList maintain_coarse_partitioner( MPI_Comm comm, int nparts, const std::vector& cell_types, - const std::vector>& cell_topology); + const std::vector>& cells); /// @brief Refine with markers, optionally redistributing, and /// optionally calculating the parent-child relationships. From 22127f7248d495fd71ed485c5490e5f60b665d9b Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Mon, 23 Sep 2024 20:03:05 +0100 Subject: [PATCH 11/12] Small updates --- cpp/dolfinx/refinement/refine.cpp | 15 ++++++-------- cpp/dolfinx/refinement/refine.h | 14 ++++++------- python/dolfinx/wrappers/refinement.cpp | 27 +++++++++++--------------- 3 files changed, 24 insertions(+), 32 deletions(-) diff --git a/cpp/dolfinx/refinement/refine.cpp b/cpp/dolfinx/refinement/refine.cpp index 2b5ce85d3c3..3fba0e2410e 100644 --- a/cpp/dolfinx/refinement/refine.cpp +++ b/cpp/dolfinx/refinement/refine.cpp @@ -5,21 +5,18 @@ // SPDX-License-Identifier: LGPL-3.0-or-later #include "refine.h" +#include -namespace dolfinx::refinement -{ +using namespace dolfinx; -graph::AdjacencyList maintain_coarse_partitioner( - MPI_Comm comm, int, const std::vector& cell_types, - const std::vector>& cells) +graph::AdjacencyList refinement::maintain_coarse_partitioner( + MPI_Comm comm, int /*nparts*/, std::vector cell_types, + std::vector> cells) { - int mpi_rank = MPI::rank(comm); int num_cell_vertices = mesh::num_cell_vertices(cell_types.front()); std::int32_t num_cells = cells.front().size() / num_cell_vertices; - std::vector destinations(num_cells, mpi_rank); + std::vector destinations(num_cells, dolfinx::MPI::rank(comm)); std::vector dest_offsets(num_cells + 1); std::iota(dest_offsets.begin(), dest_offsets.end(), 0); return graph::AdjacencyList(std::move(destinations), std::move(dest_offsets)); } - -} // namespace dolfinx::refinement diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index c23bc0a3c30..07a5b411a86 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -20,19 +20,19 @@ namespace dolfinx::refinement { - -/// @brief Partitioner that maintains the regional destribution for a refined -/// mesh. Meaning, process local data is just maintained and no redistribution -/// happens. +/// @brief Partitioner that maintains the regional distribution for a +/// refined mesh. Meaning, process local data is just maintained and no +/// redistribution happens. /// /// @param[in] comm MPI Communicator /// @param[in] nparts Number of partitions (input has no effect) /// @param[in] cell_types Cell types in the mesh /// @param[in] cells Lists of cells of each cell type. /// @return Destination ranks for each cell on this process -graph::AdjacencyList maintain_coarse_partitioner( - MPI_Comm comm, int nparts, const std::vector& cell_types, - const std::vector>& cells); +graph::AdjacencyList +maintain_coarse_partitioner(MPI_Comm comm, int nparts, + std::vector cell_types, + std::vector> cells); /// @brief Refine with markers, optionally redistributing, and /// optionally calculating the parent-child relationships. diff --git a/python/dolfinx/wrappers/refinement.cpp b/python/dolfinx/wrappers/refinement.cpp index 7921ffb22e5..8bb89b7260d 100644 --- a/python/dolfinx/wrappers/refinement.cpp +++ b/python/dolfinx/wrappers/refinement.cpp @@ -4,12 +4,18 @@ // // SPDX-License-Identifier: LGPL-3.0-or-later +#include "MPICommWrapper.h" +#include "array.h" #include #include +#include +#include +#include +#include +#include +#include +#include #include -#include -#include - #include #include #include @@ -17,26 +23,15 @@ #include #include #include - -#include -#include -#include -#include -#include -#include -#include - -#include "MPICommWrapper.h" -#include "array.h" +#include +#include namespace nb = nanobind; namespace dolfinx_wrappers { - namespace { - using PythonCellPartitionFunction = std::function( dolfinx_wrappers::MPICommWrapper, int, From 3b74e494b7fd1bf2bbd6e4f1ca9d01a8da4aed2b Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Fri, 27 Sep 2024 09:38:33 +0100 Subject: [PATCH 12/12] Fix test failure --- python/test/unit/refinement/test_refinement.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/test/unit/refinement/test_refinement.py b/python/test/unit/refinement/test_refinement.py index 25b1d460bbc..202fed23393 100644 --- a/python/test/unit/refinement/test_refinement.py +++ b/python/test/unit/refinement/test_refinement.py @@ -211,5 +211,5 @@ def test_refine_cell_meshtag(tdim, refine_plaza_wrapper): def test_refine_ufl_cargo(): mesh = create_unit_cube(MPI.COMM_WORLD, 4, 3, 3) mesh.topology.create_entities(1) - refined_mesh, _, _ = refine(mesh, redistribute=False) + refined_mesh, _, _ = refine(mesh) assert refined_mesh.ufl_domain().ufl_cargo() != mesh.ufl_domain().ufl_cargo()