/*@{*/
/**
- * 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.
* 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
/**
* 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().
* 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.
* 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.
*
* 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
*
* 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>
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
{
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),
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();
// 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();
// 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();
}
std::move(qpoints_out),
std::move(maps_out),
std::move(missing_points_out));
- }
+ } // namespace GridTools
projects / dealii.git / commitdiff