From d4eca8caaadbed88e8940d9a29e5812c19689355 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Mon, 30 Nov 2020 11:00:26 +0100 Subject: [PATCH] Remove ExcPointNotFound, and return end() instead. --- include/deal.II/grid/grid_tools.h | 78 ++---- source/grid/grid_tools.cc | 356 ++++++------------------- source/grid/grid_tools_dof_handlers.cc | 22 +- source/particles/particle_handler.cc | 22 +- 4 files changed, 128 insertions(+), 350 deletions(-) diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index 1140d67fad..001c55618c 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -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::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 class MeshType, int spacedim> @@ -3426,24 +3418,6 @@ namespace GridTools double, << "The scaling factor must be positive, but it is " << arg1 << "."); - /** - * Exception - */ - template - DeclException1(ExcPointNotFoundInCoarseGrid, - Point, - << "The point <" << arg1 - << "> could not be found inside any of the " - << "coarse grid cells."); - /** - * Exception - */ - template - DeclException1(ExcPointNotFound, - Point, - << "The point <" << arg1 - << "> could not be found inside any of the " - << "subcells of a coarse grid cell."); /** * Exception diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 57769cc554..678eedf657 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -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(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::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::active_cell_iterator, - Point> - 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 &) - { - missing_points_out.emplace_back(0); - } - ++points_checked; - } + // Keep track of all found points + std::vector found_points(points.size(), false); - // The tree search returns - // - a bounding box covering the cell - // - the active cell iterator - std::vector< - std::pair, - typename Triangulation::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::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 &) - { - 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 &) - { - 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 &) - { - 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 &) - { - 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:: - 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 diff --git a/source/grid/grid_tools_dof_handlers.cc b/source/grid/grid_tools_dof_handlers.cc index b13ea27f35..aa66cd938a 100644 --- a/source/grid/grid_tools_dof_handlers.cc +++ b/source/grid/grid_tools_dof_handlers.cc @@ -550,9 +550,6 @@ namespace GridTools } } - AssertThrow(best_cell.first.state() == IteratorState::valid, - ExcPointNotFound(p)); - return best_cell; } @@ -574,18 +571,14 @@ namespace GridTools const double tolerance, const std::vector &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 &) - {} + 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(p)); - return best_cell; } diff --git a/source/particles/particle_handler.cc b/source/particles/particle_handler.cc index b2f7e4d3c6..a66b824bf6 100644 --- a/source/particles/particle_handler.cc +++ b/source/particles/particle_handler.cc @@ -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:: - active_cell_iterator, - Point> - 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 &) + const std::pair:: + active_cell_iterator, + Point> + 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. -- 2.39.5