From 1d914a5f19c540a67e05ee4b923798d3ef626794 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Thu, 22 Apr 2021 17:47:39 +0200 Subject: [PATCH] Apply suggestions from code review Co-authored-by: Wolfgang Bangerth --- .../incompatibilities/20201130LucaHeltai | 3 +- include/deal.II/grid/grid_tools.h | 11 +- source/grid/grid_tools.cc | 100 +++++++++--------- 3 files changed, 60 insertions(+), 54 deletions(-) diff --git a/doc/news/changes/incompatibilities/20201130LucaHeltai b/doc/news/changes/incompatibilities/20201130LucaHeltai index 613a9af8f7..1c9399aa76 100644 --- a/doc/news/changes/incompatibilities/20201130LucaHeltai +++ b/doc/news/changes/incompatibilities/20201130LucaHeltai @@ -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.
(Luca Heltai, 2020/11/30) diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index 001c55618c..25e94b1bd2 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -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 diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 678eedf657..34d5b7b6de 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -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, 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 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::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::active_cell_iterator &cell, + const Point & 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>({ref_point})); + maps_out.emplace_back(std::vector({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 + } -- 2.39.5