]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Remove ExcPointNotFound, and return end() instead.
authorLuca Heltai <luca.heltai@sissa.it>
Mon, 30 Nov 2020 10:00:26 +0000 (11:00 +0100)
committerLuca Heltai <luca.heltai@sissa.it>
Tue, 11 May 2021 16:55:05 +0000 (18:55 +0200)
include/deal.II/grid/grid_tools.h
source/grid/grid_tools.cc
source/grid/grid_tools_dof_handlers.cc
source/particles/particle_handler.cc

index 1140d67fad11f47c2e26280cdfbf24c74903b407..001c55618c2788fe4c76eecabad4e1ddc1176807 100644 (file)
@@ -861,7 +861,10 @@ namespace GridTools
   /*@{*/
 
   /**
-   * Given a Triangulation's @p cache and a list of @p points create the quadrature rules.
+   * 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
+   * to global indices into @p points .
    *
    * @param[in] cache The triangulation's GridTools::Cache .
    * @param[in] points The point's vector.
@@ -883,14 +886,14 @@ namespace GridTools
    * Mapping::transform_unit_to_real(qpoints[c][0])
    * returns @p points[a].
    *
-   * The algorithm assumes it's easier to look for a point in the cell that was
-   * used previously. For this reason random points are, computationally
-   * speaking, the worst case scenario while points grouped by the cell to which
-   * they belong are the best case. Pre-sorting points, trying to minimize
-   * distances between them, might make the function extremely faster.
+   * The algorithm builds an rtree of @p points to sort them spatially, before
+   * attempting to call find_active_cell_around_point().
    *
    * @note If a point is not found inside the mesh, or is lying inside an
-   * artificial cell of a parallel::TriangulationBase, an exception is thrown.
+   * artificial cell of a parallel::TriangulationBase, the point is silently
+   * igored. 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.
    *
    * @note The actual return type of this function, i.e., the type referenced
    * above as @p return_type, is
@@ -927,7 +930,10 @@ namespace GridTools
 
   /**
    * This function is similar to GridTools::compute_point_locations(),
-   * but it tries to find and transform every point of @p points.
+   * but while compute_point_locations() silently ignores all points for which
+   * find_active_cell_around_point() fails, this function also returns a
+   * vector containing the indices of the points for which
+   * find_active_cell_around_point() failed.
    *
    * @return A tuple containing four elements; the first three
    * are documented in GridTools::compute_point_locations().
@@ -1280,11 +1286,6 @@ namespace GridTools
    * adjacent cells.
    * @return A vector of cells that lie adjacent to the given vertex.
    *
-   * @note If the point requested does not lie in any of the cells of the mesh
-   * given, then this function throws an exception of type
-   * GridTools::ExcPointNotFound. You can catch this exception and decide what
-   * to do in that case.
-   *
    * @note It isn't entirely clear at this time whether the function does the
    * right thing with anisotropically refined meshes. It needs to be checked
    * for this case.
@@ -1318,14 +1319,8 @@ namespace GridTools
    * tries to identify the cell that is of highest refinement level.
    *
    * If the point requested does not lie in a locally-owned or ghost cell,
-   * then this function throws an exception of type GridTools::ExcPointNotFound.
-   * You can catch this exception and decide what to do in that case. Hence,
-   * for programs that work with partitioned (parallel) triangulations, this
-   * function should always be called inside a `try`-block unless it is a
-   * priori clear that the point with which it is called must be inside
-   * a locally owned or ghost cell (and not close enough to the boundary
-   * between ghost and artificial cells so that decision which cell it is
-   * on depends on floating point accuracy).
+   * then this function will return the (invalid) MeshType<dim, spacedim>::end()
+   * iterator.
    *
    * @param mapping The mapping used to determine whether the given point is
    *   inside a given cell.
@@ -1447,18 +1442,15 @@ namespace GridTools
    *
    * for(auto p : points)
    * {
-   *   try
-   *   {
-   *     auto cell_and_ref_point = GridTools::find_active_cell_around_point(
-   *       cache, p, cell_hint, marked_vertices, tolerance);
-   *
-   *     // use current cell as hint for the next point
-   *     cell_hint = cell_and_ref_point.first;
-   *   }
-   *   catch(...)
-   *   {
+   *   auto cell_and_ref_point = GridTools::find_active_cell_around_point(
+   *     cache, p, cell_hint, marked_vertices, tolerance);
+   *
+   *   if(cell_and_ref_point.first != triangulation.end()) {
+   *      // use current cell as hint for the next point
+   *      cell_hint = cell_and_ref_point.first;
+   *      // do something with cell_and_ref_point
+   *      ...
    *   }
-   *
    *   ...
    * }
    * @endcode
@@ -1525,9 +1517,9 @@ namespace GridTools
    *
    * This function is used as follows
    * @code
-   *   auto first_cell = GridTools::find_active_cell_around_point(...);
+   *   auto first_pair = GridTools::find_active_cell_around_point(...);
    *   auto all_cells  = GridTools::find_all_active_cells_around_point(
-   *                              mapping, mesh, p, tolerance, first_cell);
+   *                              mapping, mesh, p, tolerance, first_pair);
    * @endcode
    */
   template <int dim, template <int, int> class MeshType, int spacedim>
@@ -3426,24 +3418,6 @@ namespace GridTools
                  double,
                  << "The scaling factor must be positive, but it is " << arg1
                  << ".");
-  /**
-   * Exception
-   */
-  template <int N>
-  DeclException1(ExcPointNotFoundInCoarseGrid,
-                 Point<N>,
-                 << "The point <" << arg1
-                 << "> could not be found inside any of the "
-                 << "coarse grid cells.");
-  /**
-   * Exception
-   */
-  template <int N>
-  DeclException1(ExcPointNotFound,
-                 Point<N>,
-                 << "The point <" << arg1
-                 << "> could not be found inside any of the "
-                 << "subcells of a coarse grid cell.");
 
   /**
    * Exception
index 57769cc554399dc2da55623941b5bc32ea7ec0b4..678eedf6573adb8e8f43ff15c2bc4571d2d577dd 100644 (file)
@@ -5431,16 +5431,9 @@ namespace GridTools
   {
     const auto cqmp = compute_point_locations_try_all(cache, points, cell_hint);
     // Splitting the tuple's components
-    auto &cells          = std::get<0>(cqmp);
-    auto &qpoints        = std::get<1>(cqmp);
-    auto &maps           = std::get<2>(cqmp);
-    auto &missing_points = std::get<3>(cqmp);
-    // If a point was not found, throwing an error, as the old
-    // implementation of compute_point_locations would have done
-    AssertThrow(std::get<3>(cqmp).size() == 0,
-                ExcPointNotFound<spacedim>(points[missing_points[0]]));
-
-    (void)missing_points;
+    auto &cells   = std::get<0>(cqmp);
+    auto &qpoints = std::get<1>(cqmp);
+    auto &maps    = std::get<2>(cqmp);
 
     return std::make_tuple(std::move(cells),
                            std::move(qpoints),
@@ -5465,6 +5458,12 @@ namespace GridTools
     const typename Triangulation<dim, spacedim>::active_cell_iterator
       &cell_hint)
   {
+    // Alias
+    namespace bgi = boost::geometry::index;
+
+    // Get the mapping
+    const auto &mapping = cache.get_mapping();
+
     // How many points are here?
     const unsigned int np = points.size();
 
@@ -5484,284 +5483,102 @@ namespace GridTools
     // For the search we shall use the following tree
     const auto &b_tree = cache.get_cell_bounding_boxes_rtree();
 
-    // We begin by finding the cell/transform of the first point
-    std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
-              Point<dim>>
-      my_pair;
+    // Now make a tree of indices for the points
+    const auto p_tree = pack_rtree_of_indices(points);
 
-    bool         found          = false;
-    unsigned int points_checked = 0;
-
-    // If a hint cell was given, use it
-    if (cell_hint.state() == IteratorState::valid)
-      {
-        try
-          {
-            my_pair = GridTools::find_active_cell_around_point(cache,
-                                                               points[0],
-                                                               cell_hint);
-            found   = true;
-          }
-        catch (const GridTools::ExcPointNotFound<dim> &)
-          {
-            missing_points_out.emplace_back(0);
-          }
-        ++points_checked;
-      }
+    // Keep track of all found points
+    std::vector<bool> found_points(points.size(), false);
 
-    // The tree search returns
-    // - a bounding box covering the cell
-    // - the active cell iterator
-    std::vector<
-      std::pair<BoundingBox<spacedim>,
-                typename Triangulation<dim, spacedim>::active_cell_iterator>>
-      box_cell;
-
-    // This is used as an index for box_cell
-    int cell_candidate_idx = -1;
-    // If any of the cells in box_cell is a ghost cell,
-    // an artificial cell or at the boundary,
-    // we want to use try/catch
-    bool use_try = false;
-
-    while (!found && points_checked < np)
-      {
-        box_cell.clear();
-        b_tree.query(boost::geometry::index::intersects(points[points_checked]),
-                     std::back_inserter(box_cell));
+    // 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];
+    };
 
-        // Checking box_cell result for a suitable candidate
-        cell_candidate_idx = -1;
-        for (unsigned int i = 0; i < box_cell.size(); ++i)
+    // 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) {
+        const auto it = std::find(cells_out.rbegin(), cells_out.rend(), cell);
+        if (it != cells_out.rend())
           {
-            // As a candidate we don't want artificial cells
-            if (!box_cell[i].second->is_artificial())
-              cell_candidate_idx = i;
-
-            // If the cell is not locally owned or at boundary
-            // we check for exceptions
-            if (cell_candidate_idx != -1 &&
-                (!box_cell[i].second->is_locally_owned() ||
-                 box_cell[i].second->at_boundary()))
-              use_try = true;
-
-
-            if (cell_candidate_idx != -1)
-              break;
-          }
-
-        // If a suitable cell was found, use it as hint
-        if (cell_candidate_idx != -1)
-          {
-            if (use_try)
-              {
-                try
-                  {
-                    my_pair = GridTools::find_active_cell_around_point(
-                      cache,
-                      points[points_checked],
-                      box_cell[cell_candidate_idx].second);
-                    found = true;
-                  }
-                catch (const GridTools::ExcPointNotFound<dim> &)
-                  {
-                    missing_points_out.emplace_back(points_checked);
-                  }
-              }
-            else
-              {
-                my_pair = GridTools::find_active_cell_around_point(
-                  cache,
-                  points[points_checked],
-                  box_cell[cell_candidate_idx].second);
-                found = true;
-              }
+            return (cells_out.size() - 1 - (it - cells_out.rbegin()));
           }
         else
           {
-            try
-              {
-                my_pair = GridTools::find_active_cell_around_point(
-                  cache, points[points_checked]);
-                // If we arrive here the cell was not among
-                // the candidates returned by the tree, so we're adding it
-                // by hand
-                found              = true;
-                cell_candidate_idx = box_cell.size();
-                box_cell.push_back(
-                  std::make_pair(my_pair.first->bounding_box(), my_pair.first));
-              }
-            catch (const GridTools::ExcPointNotFound<dim> &)
-              {
-                missing_points_out.emplace_back(points_checked);
-              }
+            cells_out.emplace_back(cell);
+            qpoints_out.resize(cells_out.size());
+            maps_out.resize(cells_out.size());
+            return (cells_out.size() - 1);
           }
+      };
 
-        // Updating the position of the analyzed points
-        ++points_checked;
-      }
-
-    // If the point has been found in a cell, adding it
-    if (found)
-      {
-        cells_out.emplace_back(my_pair.first);
-        qpoints_out.emplace_back(1, my_pair.second);
-        maps_out.emplace_back(1, points_checked - 1);
-      }
-
-    // Now the second easy case.
-    if (np == qpoints_out.size())
-      return std::make_tuple(std::move(cells_out),
-                             std::move(qpoints_out),
-                             std::move(maps_out),
-                             std::move(missing_points_out));
+    // 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;
 
-    // Cycle over all points left
-    for (unsigned int p = points_checked; p < np; ++p)
-      {
-        // We assume the last used cell contains the point: checking it
-        if (cell_candidate_idx != -1)
-          if (!box_cell[cell_candidate_idx].first.point_inside(points[p]))
-            // Point outside candidate cell: we have no candidate
-            cell_candidate_idx = -1;
-
-        // If there's no candidate, run a tree search
-        if (cell_candidate_idx == -1)
+        for (const auto id :
+             p_tree | bgi::adaptors::queried(!bgi::satisfies(already_found) &&
+                                             bgi::intersects(box)))
           {
-            // Using the b_tree to find new candidates
-            box_cell.clear();
-            b_tree.query(boost::geometry::index::intersects(points[p]),
-                         std::back_inserter(box_cell));
-            // Checking the returned bounding boxes/cells
-            use_try            = false;
-            cell_candidate_idx = -1;
-            for (unsigned int i = 0; i < box_cell.size(); ++i)
+            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)
               {
-                // As a candidate we don't want artificial cells
-                if (!box_cell[i].second->is_artificial())
-                  cell_candidate_idx = i;
-
-                // If the cell is not locally owned or at boundary
-                // we check for exceptions
-                if (cell_candidate_idx != -1 &&
-                    (!box_cell[i].second->is_locally_owned() ||
-                     box_cell[i].second->at_boundary()))
-                  use_try = true;
-
-                // If a cell candidate was found we can stop
-                if (cell_candidate_idx != -1)
-                  break;
-              }
-          }
-
-        if (cell_candidate_idx == -1)
-          {
-            // No candidate cell, but the cell might
-            // still be inside the mesh, this is our final check:
-            try
-              {
-                my_pair =
-                  GridTools::find_active_cell_around_point(cache, points[p]);
-                // If we arrive here the cell was not among
-                // the candidates returned by the tree, so we're adding it
-                // by hand
-                cell_candidate_idx = box_cell.size();
-                box_cell.push_back(
-                  std::make_pair(my_pair.first->bounding_box(), my_pair.first));
-              }
-            catch (const GridTools::ExcPointNotFound<dim> &)
-              {
-                missing_points_out.emplace_back(p);
-                continue;
-              }
-          }
-        else
-          {
-            // We have a candidate cell
-            if (use_try)
-              {
-                try
-                  {
-                    my_pair = GridTools::find_active_cell_around_point(
-                      cache, points[p], box_cell[cell_candidate_idx].second);
-                  }
-                catch (const GridTools::ExcPointNotFound<dim> &)
-                  {
-                    missing_points_out.push_back(p);
-                    continue;
-                  }
+                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
               {
-                my_pair = GridTools::find_active_cell_around_point(
-                  cache, points[p], box_cell[cell_candidate_idx].second);
-              }
-
-            // If the point was found in another cell,
-            // updating cell_candidate_idx
-            if (my_pair.first != box_cell[cell_candidate_idx].second)
-              {
-                for (unsigned int i = 0; i < box_cell.size(); ++i)
-                  {
-                    if (my_pair.first == box_cell[i].second)
-                      {
-                        cell_candidate_idx = i;
-                        break;
-                      }
-                  }
-
-                if (my_pair.first != box_cell[cell_candidate_idx].second)
-                  {
-                    // The cell was not among the candidates returned by the
-                    // tree
-                    cell_candidate_idx = box_cell.size();
-                    box_cell.push_back(
-                      std::make_pair(my_pair.first->bounding_box(),
-                                     my_pair.first));
-                  }
+                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)
+      check_all_points_within_box(
+        std::make_pair(mapping.get_bounding_box(cell_hint), cell_hint));
 
-        // Assuming the point is more likely to be in the last
-        // used cell
-        if (my_pair.first == cells_out.back())
-          {
-            // Found in the last cell: adding the data
-            qpoints_out.back().emplace_back(my_pair.second);
-            maps_out.back().emplace_back(p);
-          }
-        else
-          {
-            // Check if it is in another cell already found
-            typename std::vector<typename Triangulation<dim, spacedim>::
-                                   active_cell_iterator>::iterator cells_it =
-              std::find(cells_out.begin(), cells_out.end() - 1, my_pair.first);
-
-            if (cells_it == cells_out.end() - 1)
-              {
-                // Cell not found: adding a new cell
-                cells_out.emplace_back(my_pair.first);
-                qpoints_out.emplace_back(1, my_pair.second);
-                maps_out.emplace_back(1, p);
-              }
-            else
-              {
-                // Cell found: just adding the point index and qpoint to the
-                // list
-                unsigned int current_cell = cells_it - cells_out.begin();
-                qpoints_out[current_cell].emplace_back(my_pair.second);
-                maps_out[current_cell].emplace_back(p);
-              }
-          }
-      }
+    // Now loop over all points that have not been found yet
+    for (unsigned int i = 0; i < np; ++i)
+      if (found_points[i] == false)
+        {
+          // Get the closest cell to this point
+          const auto leaf = b_tree.qbegin(bgi::nearest(points[i], 1));
+          // Now checks all points that fall within this box
+          if (leaf != b_tree.qend())
+            check_all_points_within_box(*leaf);
+          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."));
+            }
+        }
+    // Now make sure we send out the rest of the points that we did not find.
+    for (unsigned int i = 0; i < np; ++i)
+      if (found_points[i] == false)
+        missing_points_out.emplace_back(i);
 
     // Debug Checking
-    Assert(cells_out.size() == maps_out.size(),
-           ExcDimensionMismatch(cells_out.size(), maps_out.size()));
-
-    Assert(cells_out.size() == qpoints_out.size(),
-           ExcDimensionMismatch(cells_out.size(), qpoints_out.size()));
+    AssertDimension(cells_out.size(), maps_out.size());
+    AssertDimension(cells_out.size(), qpoints_out.size());
 
 #ifdef DEBUG
     unsigned int c   = cells_out.size();
@@ -5773,8 +5590,7 @@ namespace GridTools
     // plus the points which were ignored
     for (unsigned int n = 0; n < c; ++n)
       {
-        Assert(qpoints_out[n].size() == maps_out[n].size(),
-               ExcDimensionMismatch(qpoints_out[n].size(), maps_out[n].size()));
+        AssertDimension(qpoints_out[n].size(), maps_out[n].size());
         qps += qpoints_out[n].size();
       }
 
@@ -5786,7 +5602,7 @@ namespace GridTools
                            std::move(qpoints_out),
                            std::move(maps_out),
                            std::move(missing_points_out));
-  }
+  } // namespace GridTools
 
 
 
index b13ea27f35fe42ff92b155313af5df067c329cc1..aa66cd938a0aa7c47104205fe7a3d2a1a6744d86 100644 (file)
@@ -550,9 +550,6 @@ namespace GridTools
           }
       }
 
-    AssertThrow(best_cell.first.state() == IteratorState::valid,
-                ExcPointNotFound<spacedim>(p));
-
     return best_cell;
   }
 
@@ -574,18 +571,14 @@ namespace GridTools
                                      const double                   tolerance,
                                      const std::vector<bool> &marked_vertices)
   {
-    try
-      {
-        const auto cell_and_point = find_active_cell_around_point(
-          mapping, mesh, p, marked_vertices, tolerance);
+    const auto cell_and_point = find_active_cell_around_point(
+      mapping, mesh, p, marked_vertices, tolerance);
 
-        return find_all_active_cells_around_point(
-          mapping, mesh, p, tolerance, cell_and_point);
-      }
-    catch (ExcPointNotFound<spacedim> &)
-      {}
+    if (cell_and_point.first.state() != IteratorState::valid)
+      return {};
 
-    return {};
+    return find_all_active_cells_around_point(
+      mapping, mesh, p, tolerance, cell_and_point);
   }
 
 
@@ -1384,9 +1377,6 @@ namespace GridTools
           }
       }
 
-    AssertThrow(best_cell.first.state() == IteratorState::valid,
-                ExcPointNotFound<spacedim>(p));
-
     return best_cell;
   }
 
index b2f7e4d3c6d21dddeebc1951df73e3c243ced763..a66b824bf68c8f59375c6be1f311437eb4dd0408 100644 (file)
@@ -1088,18 +1088,16 @@ namespace Particles
               // The particle is not in a neighbor of the old cell.
               // Look for the new cell in the whole local domain.
               // This case is rare.
-              try
-                {
-                  const std::pair<const typename Triangulation<dim, spacedim>::
-                                    active_cell_iterator,
-                                  Point<dim>>
-                    current_cell_and_position =
-                      GridTools::find_active_cell_around_point<>(
-                        *mapping, *triangulation, out_particle->get_location());
-                  current_cell               = current_cell_and_position.first;
-                  current_reference_position = current_cell_and_position.second;
-                }
-              catch (GridTools::ExcPointNotFound<spacedim> &)
+              const std::pair<const typename Triangulation<dim, spacedim>::
+                                active_cell_iterator,
+                              Point<dim>>
+                current_cell_and_position =
+                  GridTools::find_active_cell_around_point<>(
+                    *mapping, *triangulation, out_particle->get_location());
+              current_cell               = current_cell_and_position.first;
+              current_reference_position = current_cell_and_position.second;
+
+              if (current_cell.state() != IteratorState::valid)
                 {
                   // We can find no cell for this particle. It has left the
                   // domain due to an integration error or an open boundary.

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.