]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Choose right predicate with std::is_same
authorMarco Feder <marco.feder@sissa.it>
Fri, 30 Jun 2023 18:19:35 +0000 (20:19 +0200)
committerMarco Feder <marco.feder@sissa.it>
Fri, 14 Jul 2023 07:58:27 +0000 (09:58 +0200)
source/grid/grid_tools.cc

index d40be81d755b8705ab0c4805eb0fdeafc29af8a7..82a58c6e0aa087c161970d1d2db0e15eb5a4debe 100644 (file)
@@ -5888,90 +5888,136 @@ namespace GridTools
       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());
@@ -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<double> 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<unsigned int> bboxes_position(n_procs);
 
     unsigned int tot_bboxes = 0;

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.