# ifndef _MSC_VER
typename MeshType<dim, spacedim>::active_cell_iterator
# else
- typename dealii::internal::ActiveCellIterator<dim,
- spacedim,
- MeshType<dim, spacedim>>::type
+ typename dealii::internal::
+ ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type
# endif
find_active_cell_around_point(const MeshType<dim, spacedim> &mesh,
const Point<spacedim> & p,
- const std::vector<bool> &marked_vertices = {});
+ const std::vector<bool> &marked_vertices = {},
+ const double tolerance = 1.e-10);
/**
- * Find and return an iterator to the active cell that surrounds a given
- * point @p p.
+ * Find an active cell that surrounds a given point @p p. The return type
+ * is a pair of an iterator to the active cell along with the unit cell
+ * coordinates of the point.
*
* The algorithm used in this function proceeds by first looking for the
* vertex located closest to the given point, see
* only search among @p marked_vertices for the closest vertex.
* The size of this array should be equal to n_vertices() of the
* triangulation (as opposed to n_used_vertices() ).
+ * @param tolerance Tolerance in terms of unit cell coordinates. Depending
+ * on the problem, it might be necessary to adjust the tolerance in order
+ * to be able to identify a cell. Floating
+ * point arithmetic implies that a point will, in general, not lie exactly
+ * on a vertex, edge, or face. In either case, it is not predictable which
+ * of the cells adjacent to a vertex or an edge/face this function returns.
+ * Consequently, algorithms that call this function need to take into
+ * account that the returned cell will only contain the point approximately.
*
* @return A pair of an iterators into the mesh that points to the
* surrounding cell, and of the coordinates of that point inside the cell in
* @ref GlossGhostCell).
* If so, many of the operations one may want to do on this cell (e.g.,
* evaluating the solution) may not be possible and you will have to decide
- * what to do in that case.
- *
- * @note Floating point arithmetic implies that a point will, in general,
- * never lie <i>exactly</i> on an edge or a face. It may, however, lie
- * on a vertex of a cell. In either case, it is not predictable which
- * of the cells adjacent to a vertex or an edge/face this function returns
- * when given a point that lies on a vertex or within floating point
- * precision of an edge or face. Consequently, algorithms that call
- * this function need to take into account that the returned cell
- * will only contain the point approximately (to within round-off error)
- * and that these cells may also be ghost cells or artificial cells
- * if the triangulation is a parallel one. The latter may even be true
- * if the given point is in fact a vertex of a locally owned cell: the
- * returned cell may still be a ghost cell that happens to share this
- * vertex with a locally owned one. The reason for this behavior is that
- * it is the only way to guarantee that all processors that participate
- * in a parallel triangulation will agree which cell contains a point.
- * In other words, two processors that own two cells that come together
- * at one vertex will return the same cell when called with this vertex.
- * One of them will then return a locally owned cell and the other one
- * a ghost cell.
+ * what to do in that case. This might even be the case if the given point is
+ * a vertex of a locally owned cell: the returned cell may still be a ghost
+ * cell that happens to share this vertex with a locally owned one. The
+ * reason for this behavior is that it is the only way to guarantee that all
+ * processors that participate in a parallel triangulation will agree which
+ * cell contains a point. In other words, two processors that own two cells
+ * that come together at one vertex will return the same cell when called
+ * with this vertex. One of them will then return a locally owned cell and
+ * the other one a ghost cell.
*/
template <int dim, template <int, int> class MeshType, int spacedim>
# ifndef _MSC_VER
find_active_cell_around_point(const Mapping<dim, spacedim> & mapping,
const MeshType<dim, spacedim> &mesh,
const Point<spacedim> & p,
- const std::vector<bool> &marked_vertices = {});
+ const std::vector<bool> &marked_vertices = {},
+ const double tolerance = 1.e-10);
/**
* A version of the previous function that exploits an already existing
typename MeshType<dim, spacedim>::active_cell_iterator(),
const std::vector<bool> & marked_vertices = {},
const RTree<std::pair<Point<spacedim>, unsigned int>> &used_vertices_rtree =
- RTree<std::pair<Point<spacedim>, unsigned int>>{});
+ RTree<std::pair<Point<spacedim>, unsigned int>>{},
+ const double tolerance = 1.e-10);
/**
* A version of the previous function where we use that mapping on a given
find_active_cell_around_point(
const hp::MappingCollection<dim, spacedim> &mapping,
const hp::DoFHandler<dim, spacedim> & mesh,
- const Point<spacedim> & p);
+ const Point<spacedim> & p,
+ const double tolerance = 1.e-10);
/**
* A version of the previous function that exploits an already existing
const Point<spacedim> & p,
const typename Triangulation<dim, spacedim>::active_cell_iterator &
cell_hint = typename Triangulation<dim, spacedim>::active_cell_iterator(),
- const std::vector<bool> &marked_vertices = {});
+ const std::vector<bool> &marked_vertices = {},
+ const double tolerance = 1.e-10);
/**
* As compared to the functions above, this function identifies all cells
}
}
}
-
-
-
- template <int dim, template <int, int> class MeshType, int spacedim>
-#ifndef _MSC_VER
- std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
- Point<dim>>
-#else
- std::pair<
- typename dealii::internal::
- ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
- Point<dim>>
-#endif
- find_active_cell_around_point_tolerance(
- const Mapping<dim, spacedim> & mapping,
- const MeshType<dim, spacedim> &mesh,
- const Point<spacedim> & p,
- const std::vector<bool> & marked_vertices,
- const double tolerance)
- {
- using active_cell_iterator = typename dealii::internal::
- ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type;
-
- // The best distance is set to the
- // maximum allowable distance from
- // the unit cell; we assume a
- // max. deviation of the given tolerance
- double best_distance = tolerance;
- int best_level = -1;
- std::pair<active_cell_iterator, Point<dim>> best_cell;
-
- // Find closest vertex and determine
- // all adjacent cells
- std::vector<active_cell_iterator> adjacent_cells_tmp =
- find_cells_adjacent_to_vertex(
- mesh, find_closest_vertex(mapping, mesh, p, marked_vertices));
-
- // Make sure that we have found
- // at least one cell adjacent to vertex.
- Assert(adjacent_cells_tmp.size() > 0, ExcInternalError());
-
- // Copy all the cells into a std::set
- std::set<active_cell_iterator> adjacent_cells(adjacent_cells_tmp.begin(),
- adjacent_cells_tmp.end());
- std::set<active_cell_iterator> searched_cells;
-
- // Determine the maximal number of cells
- // in the grid.
- // As long as we have not found
- // the cell and have not searched
- // every cell in the triangulation,
- // we keep on looking.
- const unsigned int n_active_cells =
- mesh.get_triangulation().n_active_cells();
- bool found = false;
- unsigned int cells_searched = 0;
- while (!found && cells_searched < n_active_cells)
- {
- typename std::set<active_cell_iterator>::const_iterator
- cell = adjacent_cells.begin(),
- endc = adjacent_cells.end();
- for (; cell != endc; ++cell)
- {
- if ((*cell)->is_artificial() == false)
- {
- try
- {
- const Point<dim> p_cell =
- mapping.transform_real_to_unit_cell(*cell, p);
-
- // calculate the infinity norm of
- // the distance vector to the unit cell.
- const double dist =
- GeometryInfo<dim>::distance_to_unit_cell(p_cell);
-
- // We compare if the point is inside the
- // unit cell (or at least not too far
- // outside). If it is, it is also checked
- // that the cell has a more refined state
- if ((dist < best_distance) ||
- ((dist == best_distance) &&
- ((*cell)->level() > best_level)))
- {
- found = true;
- best_distance = dist;
- best_level = (*cell)->level();
- best_cell = std::make_pair(*cell, p_cell);
- }
- }
- catch (typename MappingQGeneric<dim, spacedim>::
- ExcTransformationFailed &)
- {
- // ok, the transformation
- // failed presumably
- // because the point we
- // are looking for lies
- // outside the current
- // cell. this means that
- // the current cell can't
- // be the cell around the
- // point, so just ignore
- // this cell and move on
- // to the next
- }
- }
- }
-
- // update the number of cells searched
- cells_searched += adjacent_cells.size();
-
- // if the user provided a custom mask for vertices,
- // terminate the search without trying to expand the search
- // to all cells of the triangulation, as done below.
- if (marked_vertices.size() > 0)
- cells_searched = n_active_cells;
-
- // if we have not found the cell in
- // question and have not yet searched every
- // cell, we expand our search to
- // all the not already searched neighbors of
- // the cells in adjacent_cells. This is
- // what find_active_cell_around_point_internal
- // is for.
- if (!found && cells_searched < n_active_cells)
- {
- find_active_cell_around_point_internal<dim, MeshType, spacedim>(
- mesh, searched_cells, adjacent_cells);
- }
- }
-
- AssertThrow(best_cell.first.state() == IteratorState::valid,
- ExcPointNotFound<spacedim>(p));
-
- return best_cell;
- }
} // namespace
#endif
find_active_cell_around_point(const MeshType<dim, spacedim> &mesh,
const Point<spacedim> & p,
- const std::vector<bool> & marked_vertices)
+ const std::vector<bool> & marked_vertices,
+ const double tolerance)
{
return find_active_cell_around_point<dim, MeshType, spacedim>(
- StaticMappingQ1<dim, spacedim>::mapping, mesh, p, marked_vertices)
+ StaticMappingQ1<dim, spacedim>::mapping,
+ mesh,
+ p,
+ marked_vertices,
+ tolerance)
.first;
}
find_active_cell_around_point(const Mapping<dim, spacedim> & mapping,
const MeshType<dim, spacedim> &mesh,
const Point<spacedim> & p,
- const std::vector<bool> & marked_vertices)
+ const std::vector<bool> & marked_vertices,
+ const double tolerance)
{
- return find_active_cell_around_point_tolerance(
- mapping, mesh, p, marked_vertices, 1e-10);
+ using active_cell_iterator = typename dealii::internal::
+ ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type;
+
+ // The best distance is set to the
+ // maximum allowable distance from
+ // the unit cell; we assume a
+ // max. deviation of the given tolerance
+ double best_distance = tolerance;
+ int best_level = -1;
+ std::pair<active_cell_iterator, Point<dim>> best_cell;
+
+ // Find closest vertex and determine
+ // all adjacent cells
+ std::vector<active_cell_iterator> adjacent_cells_tmp =
+ find_cells_adjacent_to_vertex(
+ mesh, find_closest_vertex(mapping, mesh, p, marked_vertices));
+
+ // Make sure that we have found
+ // at least one cell adjacent to vertex.
+ Assert(adjacent_cells_tmp.size() > 0, ExcInternalError());
+
+ // Copy all the cells into a std::set
+ std::set<active_cell_iterator> adjacent_cells(adjacent_cells_tmp.begin(),
+ adjacent_cells_tmp.end());
+ std::set<active_cell_iterator> searched_cells;
+
+ // Determine the maximal number of cells
+ // in the grid.
+ // As long as we have not found
+ // the cell and have not searched
+ // every cell in the triangulation,
+ // we keep on looking.
+ const unsigned int n_active_cells =
+ mesh.get_triangulation().n_active_cells();
+ bool found = false;
+ unsigned int cells_searched = 0;
+ while (!found && cells_searched < n_active_cells)
+ {
+ typename std::set<active_cell_iterator>::const_iterator
+ cell = adjacent_cells.begin(),
+ endc = adjacent_cells.end();
+ for (; cell != endc; ++cell)
+ {
+ if ((*cell)->is_artificial() == false)
+ {
+ try
+ {
+ const Point<dim> p_cell =
+ mapping.transform_real_to_unit_cell(*cell, p);
+
+ // calculate the infinity norm of
+ // the distance vector to the unit cell.
+ const double dist =
+ GeometryInfo<dim>::distance_to_unit_cell(p_cell);
+
+ // We compare if the point is inside the
+ // unit cell (or at least not too far
+ // outside). If it is, it is also checked
+ // that the cell has a more refined state
+ if ((dist < best_distance) ||
+ ((dist == best_distance) &&
+ ((*cell)->level() > best_level)))
+ {
+ found = true;
+ best_distance = dist;
+ best_level = (*cell)->level();
+ best_cell = std::make_pair(*cell, p_cell);
+ }
+ }
+ catch (
+ typename MappingQGeneric<dim,
+ spacedim>::ExcTransformationFailed &)
+ {
+ // ok, the transformation
+ // failed presumably
+ // because the point we
+ // are looking for lies
+ // outside the current
+ // cell. this means that
+ // the current cell can't
+ // be the cell around the
+ // point, so just ignore
+ // this cell and move on
+ // to the next
+ }
+ }
+ }
+
+ // update the number of cells searched
+ cells_searched += adjacent_cells.size();
+
+ // if the user provided a custom mask for vertices,
+ // terminate the search without trying to expand the search
+ // to all cells of the triangulation, as done below.
+ if (marked_vertices.size() > 0)
+ cells_searched = n_active_cells;
+
+ // if we have not found the cell in
+ // question and have not yet searched every
+ // cell, we expand our search to
+ // all the not already searched neighbors of
+ // the cells in adjacent_cells. This is
+ // what find_active_cell_around_point_internal
+ // is for.
+ if (!found && cells_searched < n_active_cells)
+ {
+ find_active_cell_around_point_internal<dim, MeshType, spacedim>(
+ mesh, searched_cells, adjacent_cells);
+ }
+ }
+
+ AssertThrow(best_cell.first.state() == IteratorState::valid,
+ ExcPointNotFound<spacedim>(p));
+
+ return best_cell;
}
{
try
{
- const auto cell_and_point = find_active_cell_around_point_tolerance(
+ const auto cell_and_point = find_active_cell_around_point(
mapping, mesh, p, marked_vertices, tolerance);
return find_all_active_cells_around_point(
find_active_cell_around_point(
const hp::MappingCollection<dim, spacedim> &mapping,
const hp::DoFHandler<dim, spacedim> & mesh,
- const Point<spacedim> & p)
+ const Point<spacedim> & p,
+ const double tolerance)
{
Assert((mapping.size() == 1) ||
(mapping.size() == mesh.get_fe_collection().size()),
// we use find_active_cell_around_point using only one
// mapping.
if (mapping.size() == 1)
- best_cell = find_active_cell_around_point(mapping[0], mesh, p);
+ {
+ const std::vector<bool> marked_vertices = {};
+ best_cell = find_active_cell_around_point(
+ mapping[0], mesh, p, marked_vertices, tolerance);
+ }
else
{
// The best distance is set to the
// maximum allowable distance from
- // the unit cell; we assume a
- // max. deviation of 1e-10
- double best_distance = 1e-10;
+ // the unit cell
+ double best_distance = tolerance;
int best_level = -1;
projects / dealii.git / commitdiff