From 2460704bb6b949adef7ba16bca78ea4227fe90d6 Mon Sep 17 00:00:00 2001 From: Marco Feder Date: Fri, 30 Jun 2023 20:19:35 +0200 Subject: [PATCH] Choose right predicate with std::is_same --- source/grid/grid_tools.cc | 181 ++++++++++++++++++++++++-------------- 1 file changed, 114 insertions(+), 67 deletions(-) diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index d40be81d75..82a58c6e0a 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -5888,90 +5888,136 @@ namespace GridTools std::vector> ranks_and_indices; ranks_and_indices.reserve(entities.size()); -#if !defined(DEAL_II_WITH_ARBORX) +#if defined(DEAL_II_WITH_ARBORX) +# define USE_ARBORX true +#else +# define USE_ARBORX false +#endif - std::vector>> global_bboxes_temp; - auto *global_bboxes_to_be_used = &global_bboxes; - if (global_bboxes.size() == 1) + if constexpr (USE_ARBORX == false) { - global_bboxes_temp = - Utilities::MPI::all_gather(comm, global_bboxes[0]); - global_bboxes_to_be_used = &global_bboxes_temp; - } + std::vector>> global_bboxes_temp; + auto *global_bboxes_to_be_used = &global_bboxes; + if (global_bboxes.size() == 1) + { + global_bboxes_temp = + Utilities::MPI::all_gather(comm, global_bboxes[0]); + global_bboxes_to_be_used = &global_bboxes_temp; + } - // helper function to determine if a bounding box is valid - const auto is_valid = [](const auto &bb) { - for (unsigned int i = 0; i < spacedim; ++i) - if (bb.get_boundary_points().first[i] > - bb.get_boundary_points().second[i]) - return false; + // helper function to determine if a bounding box is valid + const auto is_valid = [](const auto &bb) { + for (unsigned int i = 0; i < spacedim; ++i) + if (bb.get_boundary_points().first[i] > + bb.get_boundary_points().second[i]) + return false; - return true; - }; + return true; + }; - // linearize vector of vectors - std::vector, unsigned int>> - boxes_and_ranks; + // linearize vector of vectors + std::vector, unsigned int>> + boxes_and_ranks; - for (unsigned rank = 0; rank < global_bboxes_to_be_used->size(); ++rank) - for (const auto &box : (*global_bboxes_to_be_used)[rank]) - if (is_valid(box)) - boxes_and_ranks.emplace_back(box, rank); + for (unsigned rank = 0; rank < global_bboxes_to_be_used->size(); + ++rank) + for (const auto &box : (*global_bboxes_to_be_used)[rank]) + if (is_valid(box)) + boxes_and_ranks.emplace_back(box, rank); - // pack boxes into r-tree - const auto tree = pack_rtree(boxes_and_ranks); + // pack boxes into r-tree + const auto tree = pack_rtree(boxes_and_ranks); - // loop over all entities - for (unsigned int i = 0; i < entities.size(); ++i) - { - // create a bounding box with tolerance - const auto bb = - BoundingBox(entities[i]).create_extended(tolerance); + // loop over all entities + for (unsigned int i = 0; i < entities.size(); ++i) + { + // create a bounding box with tolerance + const auto bb = + BoundingBox(entities[i]).create_extended(tolerance); - // determine ranks potentially owning point/bounding box - std::set my_ranks; + // determine ranks potentially owning point/bounding box + std::set my_ranks; - for (const auto &box_and_rank : - tree | boost::geometry::index::adaptors::queried( - boost::geometry::index::intersects(bb))) - my_ranks.insert(box_and_rank.second); + for (const auto &box_and_rank : + tree | boost::geometry::index::adaptors::queried( + boost::geometry::index::intersects(bb))) + my_ranks.insert(box_and_rank.second); - for (const auto rank : my_ranks) - ranks_and_indices.emplace_back(rank, i); + for (const auto rank : my_ranks) + ranks_and_indices.emplace_back(rank, i); + } } - -#else - if (global_bboxes.size() == 1) + else if constexpr (USE_ARBORX == true) { - // Create query bounding boxes out of entities - std::vector> query_bounding_boxes; - for (const auto &entity : entities) - query_bounding_boxes.push_back( - BoundingBox(entity).create_extended(tolerance)); - ArborXWrappers::BoundingBoxIntersectPredicate bb_intersect( - query_bounding_boxes); + auto fill_ranks_and_indices = + [&ranks_and_indices]( + const std::vector & offsets, + const std::vector> &indices_ranks) { + for (unsigned int i = 0; i < offsets.size() - 1; ++i) + for (int j = offsets[i]; j < offsets[i + 1]; ++j) + ranks_and_indices.emplace_back(indices_ranks[j].second, i); + }; - // Create bounding volume hierarchy out of input vector of bounding - ArborXWrappers::DistributedTree distributed_tree(comm, - global_bboxes[0]); + unsigned int idx = (global_bboxes.size() == 1) ? + 0 : + Utilities::MPI::this_mpi_process(comm); + // size of global_bboxes>1 + if constexpr (spacedim == 2) + { + // Create query bounding boxes out of entities + std::vector> query_bounding_boxes; + for (const auto &entity : entities) + query_bounding_boxes.push_back( + BoundingBox<2>(entity).create_extended(tolerance)); + ArborXWrappers::BoundingBoxIntersectPredicate bb_intersect( + query_bounding_boxes); + + ArborXWrappers::DistributedTree distributed_tree( + comm, global_bboxes[idx]); + + auto indices_ranks_offset = distributed_tree.query(bb_intersect); + auto indices_ranks = indices_ranks_offset.first; + auto offsets = indices_ranks_offset.second; + + fill_ranks_and_indices(offsets, indices_ranks); + } + else if constexpr (spacedim == 3) + { + ArborXWrappers::DistributedTree distributed_tree( + comm, global_bboxes[idx]); - auto indices_ranks_offset = distributed_tree.query(bb_intersect); + if constexpr (std::is_same>::value) + { + ArborXWrappers::PointIntersectPredicate bb_intersect( + entities); - auto indices_ranks = indices_ranks_offset.first; - auto offsets = indices_ranks_offset.second; + auto indices_ranks_offset = + distributed_tree.query(bb_intersect); + auto indices_ranks = indices_ranks_offset.first; + auto offsets = indices_ranks_offset.second; - std::set my_ranks; - for (unsigned int i = 0; i < offsets.size() - 1; ++i) - { - for (int j = offsets[i]; j < offsets[i + 1]; ++j) - my_ranks.insert(indices_ranks[j].second); - for (const auto rank : my_ranks) - ranks_and_indices.emplace_back(rank, i); - my_ranks.clear(); + fill_ranks_and_indices(offsets, indices_ranks); + } + else if constexpr (std::is_same>::value) + { + ArborXWrappers::BoundingBoxIntersectPredicate bb_intersect( + entities); + + auto indices_ranks_offset = + distributed_tree.query(bb_intersect); + auto indices_ranks = indices_ranks_offset.first; + auto offsets = indices_ranks_offset.second; + + fill_ranks_and_indices(offsets, indices_ranks); + } } } -#endif + else + { + // We should never get here: + AssertThrow(false, ExcInternalError()); + } // convert to CRS std::sort(ranks_and_indices.begin(), ranks_and_indices.end()); @@ -6912,7 +6958,8 @@ namespace GridTools unsigned int n_bboxes = local_bboxes.size(); // Dimension of the array to be exchanged (number of double) int n_local_data = 2 * spacedim * n_bboxes; - // data array stores each entry of each point describing the bounding boxes + // data array stores each entry of each point describing the bounding + // boxes std::vector loc_data_array(n_local_data); for (unsigned int i = 0; i < n_bboxes; ++i) for (unsigned int d = 0; d < spacedim; ++d) @@ -7013,8 +7060,8 @@ namespace GridTools // Preparing to flatten the vector const unsigned int n_procs = Utilities::MPI::n_mpi_processes(mpi_communicator); - // The i'th element of the following vector contains the index of the first - // local bounding box from the process of rank i + // The i'th element of the following vector contains the index of the + // first local bounding box from the process of rank i std::vector bboxes_position(n_procs); unsigned int tot_bboxes = 0; -- 2.39.5