From 042dce82a5f9de0e30cc722de84744ddcabe692b Mon Sep 17 00:00:00 2001 From: Johannes Heinz Date: Fri, 28 Apr 2023 11:02:36 +0200 Subject: [PATCH] make it possible to fill RPE with GT::DistributedComputePointLocationsInternal --- .../base/mpi_remote_point_evaluation.h | 24 ++ include/deal.II/grid/grid_tools.h | 16 ++ source/base/mpi_remote_point_evaluation.cc | 25 +- source/grid/grid_tools.cc | 219 ++++++++++-------- source/grid/grid_tools.inst.in | 4 + 5 files changed, 186 insertions(+), 102 deletions(-) diff --git a/include/deal.II/base/mpi_remote_point_evaluation.h b/include/deal.II/base/mpi_remote_point_evaluation.h index 7dfbf220f9..b357d28f6d 100644 --- a/include/deal.II/base/mpi_remote_point_evaluation.h +++ b/include/deal.II/base/mpi_remote_point_evaluation.h @@ -25,6 +25,15 @@ DEAL_II_NAMESPACE_OPEN +// Forward declarations +namespace GridTools +{ + namespace internal + { + template + struct DistributedComputePointLocationsInternal; + } +} // namespace GridTools namespace Utilities { @@ -86,6 +95,21 @@ namespace Utilities const Triangulation &tria, const Mapping & mapping); + /** + * Set up internal data structures and communication pattern based on + * GridTools::internal::DistributedComputePointLocationsInternal. + * + * This function is called internally by the reinit() function above. + * Having it as a seperate function makes it possible to setup the class + * if it is known in which cells corresponding reference points are + * located (e.g. if intersections of cells are known). + */ + void + reinit(const GridTools::internal:: + DistributedComputePointLocationsInternal &data, + const Triangulation & tria, + const Mapping &mapping); + /** * Data of points positioned in a cell. */ diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index 7bf25de6f9..d48e80213d 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -1213,6 +1213,22 @@ namespace GridTools template struct DistributedComputePointLocationsInternal { + DistributedComputePointLocationsInternal(); + + /** + * Function which sets up @p send_ranks, @p send_ptrs, @p recv_ranks, + * and @p recv_ptrs from @p send_components, @p recv_components, + * and @p n_searched_points. Internally @p send_components and @p recv_components + * are sorted and enumerated. + */ + void + finalize_setup(); + + /** + * Number of searched point locations. + */ + unsigned int n_searched_points; + /** * Information of each point on sending/evaluation side. The elements of * the tuple are as follows: 0) cell level and index, 1) rank of the diff --git a/source/base/mpi_remote_point_evaluation.cc b/source/base/mpi_remote_point_evaluation.cc index 656be3e2a6..43565aa536 100644 --- a/source/base/mpi_remote_point_evaluation.cc +++ b/source/base/mpi_remote_point_evaluation.cc @@ -78,9 +78,6 @@ namespace Utilities tria_signal = tria.signals.any_change.connect([&]() { this->ready_flag = false; }); - this->tria = &tria; - this->mapping = &mapping; - std::vector> local_boxes; for (const auto &cell : tria.active_cell_iterators() | IteratorFilters::LocallyOwnedCell()) @@ -105,6 +102,23 @@ namespace Utilities true, enforce_unique_mapping); + this->reinit(data, tria, mapping); +#endif + } + + + + template + void + RemotePointEvaluation::reinit( + const GridTools::internal:: + DistributedComputePointLocationsInternal &data, + const Triangulation & tria, + const Mapping & mapping) + { + this->tria = &tria; + this->mapping = &mapping; + this->recv_ranks = data.recv_ranks; this->recv_ptrs = data.recv_ptrs; @@ -113,7 +127,7 @@ namespace Utilities this->recv_permutation = {}; this->recv_permutation.resize(data.recv_components.size()); - this->point_ptrs.assign(points.size() + 1, 0); + this->point_ptrs.assign(data.n_searched_points + 1, 0); for (unsigned int i = 0; i < data.recv_components.size(); ++i) { AssertIndexRange(std::get<2>(data.recv_components[i]), @@ -130,7 +144,7 @@ namespace Utilities std::tuple n_owning_processes_local = n_owning_processes_default; - for (unsigned int i = 0; i < points.size(); ++i) + for (unsigned int i = 0; i < data.n_searched_points; ++i) { std::get<0>(n_owning_processes_local) = std::min(std::get<0>(n_owning_processes_local), @@ -196,7 +210,6 @@ namespace Utilities cell_data.reference_point_values.size()); this->ready_flag = true; -#endif } diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index dfce19ff5d..8afe92502e 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -6168,6 +6168,127 @@ namespace GridTools return locally_owned_active_cells_around_point; } + template + DistributedComputePointLocationsInternal:: + DistributedComputePointLocationsInternal() + : n_searched_points(numbers::invalid_unsigned_int) + {} + + template + void + DistributedComputePointLocationsInternal::finalize_setup() + { + // before reshuffeling the data check if data.recv_components and + // n_searched_points are in a valid state. + Assert(n_searched_points != numbers::invalid_unsigned_int, + ExcInternalError()); + Assert(recv_components.empty() || + std::get<1>(*std::max_element(recv_components.begin(), + recv_components.end(), + [](const auto &a, const auto &b) { + return std::get<1>(a) < + std::get<1>(b); + })) < n_searched_points, + ExcInternalError()); + + send_ranks.clear(); + recv_ranks.clear(); + send_ptrs.clear(); + recv_ptrs.clear(); + + if (true) + { + // sort according to rank (and point index and cell) -> make + // deterministic + std::sort(send_components.begin(), + send_components.end(), + [&](const auto &a, const auto &b) { + if (std::get<1>(a) != std::get<1>(b)) // rank + return std::get<1>(a) < std::get<1>(b); + + if (std::get<2>(a) != std::get<2>(b)) // point index + return std::get<2>(a) < std::get<2>(b); + + return std::get<0>(a) < std::get<0>(b); // cell + }); + + // perform enumeration and extract rank information + for (unsigned int i = 0, dummy = numbers::invalid_unsigned_int; + i < send_components.size(); + ++i) + { + std::get<5>(send_components[i]) = i; + + if (dummy != std::get<1>(send_components[i])) + { + dummy = std::get<1>(send_components[i]); + send_ranks.push_back(dummy); + send_ptrs.push_back(i); + } + } + send_ptrs.push_back(send_components.size()); + + // sort according to cell, rank, point index (while keeping + // partial ordering) + std::sort(send_components.begin(), + send_components.end(), + [&](const auto &a, const auto &b) { + if (std::get<0>(a) != std::get<0>(b)) + return std::get<0>(a) < std::get<0>(b); // cell + + if (std::get<1>(a) != std::get<1>(b)) + return std::get<1>(a) < std::get<1>(b); // rank + + if (std::get<2>(a) != std::get<2>(b)) + return std::get<2>(a) < std::get<2>(b); // point index + + return std::get<5>(a) < std::get<5>(b); // enumeration + }); + } + + if (recv_components.size() > 0) + { + // sort according to rank (and point index) -> make deterministic + std::sort(recv_components.begin(), + recv_components.end(), + [&](const auto &a, const auto &b) { + if (std::get<0>(a) != std::get<0>(b)) + return std::get<0>(a) < std::get<0>(b); // rank + + return std::get<1>(a) < std::get<1>(b); // point index + }); + + // perform enumeration and extract rank information + for (unsigned int i = 0, dummy = numbers::invalid_unsigned_int; + i < recv_components.size(); + ++i) + { + std::get<2>(recv_components[i]) = i; + + if (dummy != std::get<0>(recv_components[i])) + { + dummy = std::get<0>(recv_components[i]); + recv_ranks.push_back(dummy); + recv_ptrs.push_back(i); + } + } + recv_ptrs.push_back(recv_components.size()); + + // sort according to point index and rank (while keeping partial + // ordering) + std::sort(recv_components.begin(), + recv_components.end(), + [&](const auto &a, const auto &b) { + if (std::get<1>(a) != std::get<1>(b)) + return std::get<1>(a) < std::get<1>(b); // point index + + if (std::get<0>(a) != std::get<0>(b)) + return std::get<0>(a) < std::get<0>(b); // rank + + return std::get<2>(a) < std::get<2>(b); // enumeration + }); + } + } template @@ -6182,13 +6303,10 @@ namespace GridTools const bool enforce_unique_mapping) { DistributedComputePointLocationsInternal result; + result.n_searched_points = points.size(); auto &send_components = result.send_components; - auto &send_ranks = result.send_ranks; - auto &send_ptrs = result.send_ptrs; auto &recv_components = result.recv_components; - auto &recv_ranks = result.recv_ranks; - auto &recv_ptrs = result.recv_ptrs; const auto comm = cache.get_triangulation().get_communicator(); @@ -6330,98 +6448,7 @@ namespace GridTools process_answer, comm); - if (true) - { - // sort according to rank (and point index and cell) -> make - // deterministic - std::sort(send_components.begin(), - send_components.end(), - [&](const auto &a, const auto &b) { - if (std::get<1>(a) != std::get<1>(b)) // rank - return std::get<1>(a) < std::get<1>(b); - - if (std::get<2>(a) != std::get<2>(b)) // point index - return std::get<2>(a) < std::get<2>(b); - - return std::get<0>(a) < std::get<0>(b); // cell - }); - - // perform enumeration and extract rank information - for (unsigned int i = 0, dummy = numbers::invalid_unsigned_int; - i < send_components.size(); - ++i) - { - std::get<5>(send_components[i]) = i; - - if (dummy != std::get<1>(send_components[i])) - { - dummy = std::get<1>(send_components[i]); - send_ranks.push_back(dummy); - send_ptrs.push_back(i); - } - } - send_ptrs.push_back(send_components.size()); - - // sort according to cell, rank, point index (while keeping - // partial ordering) - std::sort(send_components.begin(), - send_components.end(), - [&](const auto &a, const auto &b) { - if (std::get<0>(a) != std::get<0>(b)) - return std::get<0>(a) < std::get<0>(b); // cell - - if (std::get<1>(a) != std::get<1>(b)) - return std::get<1>(a) < std::get<1>(b); // rank - - if (std::get<2>(a) != std::get<2>(b)) - return std::get<2>(a) < std::get<2>(b); // point index - - return std::get<5>(a) < std::get<5>(b); // enumeration - }); - } - - if (perform_handshake) - { - // sort according to rank (and point index) -> make deterministic - std::sort(recv_components.begin(), - recv_components.end(), - [&](const auto &a, const auto &b) { - if (std::get<0>(a) != std::get<0>(b)) - return std::get<0>(a) < std::get<0>(b); // rank - - return std::get<1>(a) < std::get<1>(b); // point index - }); - - // perform enumeration and extract rank information - for (unsigned int i = 0, dummy = numbers::invalid_unsigned_int; - i < recv_components.size(); - ++i) - { - std::get<2>(recv_components[i]) = i; - - if (dummy != std::get<0>(recv_components[i])) - { - dummy = std::get<0>(recv_components[i]); - recv_ranks.push_back(dummy); - recv_ptrs.push_back(i); - } - } - recv_ptrs.push_back(recv_components.size()); - - // sort according to point index and rank (while keeping partial - // ordering) - std::sort(recv_components.begin(), - recv_components.end(), - [&](const auto &a, const auto &b) { - if (std::get<1>(a) != std::get<1>(b)) - return std::get<1>(a) < std::get<1>(b); // point index - - if (std::get<0>(a) != std::get<0>(b)) - return std::get<0>(a) < std::get<0>(b); // rank - - return std::get<2>(a) < std::get<2>(b); // enumeration - }); - } + result.finalize_setup(); return result; } diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in index ec06814639..9e5d8fe338 100644 --- a/source/grid/grid_tools.inst.in +++ b/source/grid/grid_tools.inst.in @@ -141,6 +141,10 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) const std::vector &marked_vertices, const bool enforce_unique_mapping); + template struct internal::DistributedComputePointLocationsInternal< + deal_II_dimension, + deal_II_space_dimension>; + template internal::DistributedComputePointLocationsInternal< deal_II_dimension, deal_II_space_dimension> -- 2.39.5