std::vector<std::pair<unsigned int, unsigned int>> 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<std::vector<BoundingBox<spacedim>>> 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<std::vector<BoundingBox<spacedim>>> 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<std::pair<BoundingBox<spacedim>, unsigned int>>
- boxes_and_ranks;
+ // linearize vector of vectors
+ std::vector<std::pair<BoundingBox<spacedim>, 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<spacedim>(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<spacedim>(entities[i]).create_extended(tolerance);
- // determine ranks potentially owning point/bounding box
- std::set<unsigned int> my_ranks;
+ // determine ranks potentially owning point/bounding box
+ std::set<unsigned int> 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<BoundingBox<spacedim>> query_bounding_boxes;
- for (const auto &entity : entities)
- query_bounding_boxes.push_back(
- BoundingBox<spacedim>(entity).create_extended(tolerance));
- ArborXWrappers::BoundingBoxIntersectPredicate bb_intersect(
- query_bounding_boxes);
+ auto fill_ranks_and_indices =
+ [&ranks_and_indices](
+ const std::vector<int> & offsets,
+ const std::vector<std::pair<int, int>> &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<BoundingBox<2>> 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<T, Point<3>>::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<unsigned int> 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<T, BoundingBox<3>>::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());
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<double> loc_data_array(n_local_data);
for (unsigned int i = 0; i < n_bboxes; ++i)
for (unsigned int d = 0; d < spacedim; ++d)
// 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<unsigned int> bboxes_position(n_procs);
unsigned int tot_bboxes = 0;