]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Apply suggestions from code review
authorLuca Heltai <luca.heltai@sissa.it>
Thu, 22 Apr 2021 15:47:39 +0000 (17:47 +0200)
committerLuca Heltai <luca.heltai@sissa.it>
Tue, 11 May 2021 16:55:39 +0000 (18:55 +0200)
Co-authored-by: Wolfgang Bangerth <bangerth@colostate.edu>
doc/news/changes/incompatibilities/20201130LucaHeltai
include/deal.II/grid/grid_tools.h
source/grid/grid_tools.cc

index 613a9af8f713259095aed46ca01933e1e1ba43c4..1c9399aa76f25c972ef701df60d9f3065d4b604d 100644 (file)
@@ -1,4 +1,5 @@
 Changed: GridTools::find_active_cell_around_point() no longer throws an exception, but returns an invalid iterator.
-User codes should check that the returned cell is valid instead of relying on exceptions. This solves issues with MPI, and HDF5.
+User codes should now check that the returned cell is valid instead of relying on exceptions, instead of trying
+to catch exceptions that this function may have thrown.
 <br>
 (Luca Heltai, 2020/11/30)
index 001c55618c2788fe4c76eecabad4e1ddc1176807..25e94b1bd252bf26ca3ca988adf5fdaeefa84df2 100644 (file)
@@ -862,8 +862,8 @@ namespace GridTools
 
   /**
    * Given a Triangulation's @p cache and a list of @p points, call
-   * find_active_cell_around_point() on each element of @p points , and return
-   * @p cells , referece positions @p qpoints , and a mapping @p maps from local
+   * find_active_cell_around_point() on each element of @p points, and return
+   * @p cells, reference positions @p qpoints, and a mapping @p maps from local
    * to global indices into @p points .
    *
    * @param[in] cache The triangulation's GridTools::Cache .
@@ -891,7 +891,7 @@ namespace GridTools
    *
    * @note If a point is not found inside the mesh, or is lying inside an
    * artificial cell of a parallel::TriangulationBase, the point is silently
-   * igored. If you want to infer for which points the search failed, use the
+   * ignored. If you want to infer for which points the search failed, use the
    * function compute_point_locations_try_all() that also returns a vector of
    * indices indicating the points for which the search failed.
    *
@@ -1451,6 +1451,11 @@ namespace GridTools
    *      // do something with cell_and_ref_point
    *      ...
    *   }
+   *  else
+   *    {
+   *       // The function did not find a locally owned or ghost cell in which
+   *       // the point is located. We ought to handle this somehow here.
+   *    }
    *   ...
    * }
    * @endcode
index 678eedf6573adb8e8f43ff15c2bc4571d2d577dd..34d5b7b6deb741b50e3f52d04a099491270a3fe9 100644 (file)
@@ -5484,69 +5484,72 @@ namespace GridTools
     const auto &b_tree = cache.get_cell_bounding_boxes_rtree();
 
     // Now make a tree of indices for the points
-    const auto p_tree = pack_rtree_of_indices(points);
+    // [TODO] This would work better with pack_rtree_of_indices, but
+    // windows does not like it. Build a tree with pairs of point and id
+    std::vector<std::pair<Point<spacedim>, unsigned int>> points_and_ids(np);
+    for (unsigned int i = 0; i < np; ++i)
+      points_and_ids[i] = std::make_pair(points[i], i);
+    const auto p_tree = pack_rtree(points_and_ids);
 
     // Keep track of all found points
     std::vector<bool> found_points(points.size(), false);
 
     // Check if a point was found
-    const auto already_found = [&found_points](const unsigned int &id) {
-      AssertIndexRange(id, found_points.size());
-      return found_points[id];
+    const auto already_found = [&found_points](const auto &id) {
+      AssertIndexRange(id.second, found_points.size());
+      return found_points[id.second];
     };
 
-    // check if the cell was already in the vector before. If so, returns its
-    // index, otherwise adds it to the list of cells, return the new index,
-    // and resize also the other two vectors. Start from last entry.
-    const auto get_local_cell_index =
-      [&](const typename Triangulation<dim, spacedim>::active_cell_iterator
-            &cell) {
+    // check if the given cell was already in the vector of cells before. If so,
+    // insert in the corresponding vectors the reference point and the id.
+    // Otherwise append a new entry to all vectors.
+    const auto store_cell_point_and_id =
+      [&](
+        const typename Triangulation<dim, spacedim>::active_cell_iterator &cell,
+        const Point<dim> &  ref_point,
+        const unsigned int &id) {
         const auto it = std::find(cells_out.rbegin(), cells_out.rend(), cell);
         if (it != cells_out.rend())
           {
-            return (cells_out.size() - 1 - (it - cells_out.rbegin()));
+            const auto cell_id =
+              (cells_out.size() - 1 - (it - cells_out.rbegin()));
+            qpoints_out[cell_id].emplace_back(ref_point);
+            maps_out[cell_id].emplace_back(id);
           }
         else
           {
             cells_out.emplace_back(cell);
-            qpoints_out.resize(cells_out.size());
-            maps_out.resize(cells_out.size());
-            return (cells_out.size() - 1);
+            qpoints_out.emplace_back(std::vector<Point<dim>>({ref_point}));
+            maps_out.emplace_back(std::vector<unsigned int>({id}));
           }
       };
 
     // Check all points within a given pair of box and cell
-    const auto check_all_points_within_box =
-      [&](const decltype(*b_tree.begin()) &leaf) {
-        const auto &box       = leaf.first;
-        const auto &cell_hint = leaf.second;
-
-        for (const auto id :
-             p_tree | bgi::adaptors::queried(!bgi::satisfies(already_found) &&
-                                             bgi::intersects(box)))
-          {
-            const auto cell_and_ref =
-              GridTools::find_active_cell_around_point(cache,
-                                                       points[id],
-                                                       cell_hint);
-            const auto &cell      = cell_and_ref.first;
-            const auto &ref_point = cell_and_ref.second;
-
-            if (cell.state() == IteratorState::valid)
-              {
-                const auto cell_id = get_local_cell_index(cell);
-                qpoints_out[cell_id].emplace_back(ref_point);
-                maps_out[cell_id].emplace_back(id);
-                found_points[id] = true;
-              }
-            else
-              {
-                missing_points_out.emplace_back(id);
-                // Don't look anymore for this point
-                found_points[id] = true;
-              }
-          }
-      };
+    const auto check_all_points_within_box = [&](const auto &leaf) {
+      const auto &box       = leaf.first;
+      const auto &cell_hint = leaf.second;
+
+      for (const auto &point_and_id :
+           p_tree | bgi::adaptors::queried(!bgi::satisfies(already_found) &&
+                                           bgi::intersects(box)))
+        {
+          const auto id = point_and_id.second;
+          const auto cell_and_ref =
+            GridTools::find_active_cell_around_point(cache,
+                                                     points[id],
+                                                     cell_hint);
+          const auto &cell      = cell_and_ref.first;
+          const auto &ref_point = cell_and_ref.second;
+
+          if (cell.state() == IteratorState::valid)
+            store_cell_point_and_id(cell, ref_point, id);
+          else
+            missing_points_out.emplace_back(id);
+
+          // Don't look anymore for this point
+          found_points[id] = true;
+        }
+    };
 
     // If a hint cell was given, use it
     if (cell_hint.state() == IteratorState::valid)
@@ -5565,10 +5568,7 @@ namespace GridTools
           else
             {
               // We should not get here. Throw an error.
-              Assert(false,
-                     ExcInternalError(
-                       "I cannot find the closest cell to a point."
-                       "Something is rotten inside."));
+              Assert(false, ExcInternalError());
             }
         }
     // Now make sure we send out the rest of the points that we did not find.
@@ -5602,7 +5602,7 @@ namespace GridTools
                            std::move(qpoints_out),
                            std::move(maps_out),
                            std::move(missing_points_out));
-  } // namespace GridTools
+  }
 
 
 

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.