From: Niklas Fehn Date: Sun, 7 Jun 2020 20:33:41 +0000 (+0200) Subject: GridTools::find_active_cell_around_point(): introduce tolerance as parameter X-Git-Tag: v9.3.0-rc1~1456^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e982fc9f31ca568f3b30723065825df45517505f;p=dealii.git GridTools::find_active_cell_around_point(): introduce tolerance as parameter --- diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index 1a19ca6c1b..b46b273cbc 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -1112,17 +1112,18 @@ namespace GridTools # ifndef _MSC_VER typename MeshType::active_cell_iterator # else - typename dealii::internal::ActiveCellIterator>::type + typename dealii::internal:: + ActiveCellIterator>::type # endif find_active_cell_around_point(const MeshType &mesh, const Point & p, - const std::vector &marked_vertices = {}); + const std::vector &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 @@ -1150,6 +1151,14 @@ namespace GridTools * 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 @@ -1185,27 +1194,15 @@ namespace GridTools * @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 exactly 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 class MeshType, int spacedim> # ifndef _MSC_VER @@ -1218,7 +1215,8 @@ namespace GridTools find_active_cell_around_point(const Mapping & mapping, const MeshType &mesh, const Point & p, - const std::vector &marked_vertices = {}); + const std::vector &marked_vertices = {}, + const double tolerance = 1.e-10); /** * A version of the previous function that exploits an already existing @@ -1249,7 +1247,8 @@ namespace GridTools typename MeshType::active_cell_iterator(), const std::vector & marked_vertices = {}, const RTree, unsigned int>> &used_vertices_rtree = - RTree, unsigned int>>{}); + RTree, unsigned int>>{}, + const double tolerance = 1.e-10); /** * A version of the previous function where we use that mapping on a given @@ -1263,7 +1262,8 @@ namespace GridTools find_active_cell_around_point( const hp::MappingCollection &mapping, const hp::DoFHandler & mesh, - const Point & p); + const Point & p, + const double tolerance = 1.e-10); /** * A version of the previous function that exploits an already existing @@ -1277,7 +1277,8 @@ namespace GridTools const Point & p, const typename Triangulation::active_cell_iterator & cell_hint = typename Triangulation::active_cell_iterator(), - const std::vector &marked_vertices = {}); + const std::vector &marked_vertices = {}, + const double tolerance = 1.e-10); /** * As compared to the functions above, this function identifies all cells diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 9f6beb99d2..535e1d4392 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -1580,7 +1580,8 @@ namespace GridTools const std::vector>> &vertex_to_cell_centers, const typename MeshType::active_cell_iterator &cell_hint, const std::vector & marked_vertices, - const RTree, unsigned int>> &used_vertices_rtree) + const RTree, unsigned int>> &used_vertices_rtree, + const double tolerance) { std::pair::active_cell_iterator, Point> @@ -1670,9 +1671,8 @@ namespace GridTools comp); // It is possible the vertex is close // to an edge, thus we add a tolerance - // setting it initially to 1e-10 // to keep also the "best" cell - double best_distance = 1e-10; + double best_distance = tolerance; // Search all of the cells adjacent to the closest vertex of the cell // hint Most likely we will find the point in them. @@ -1684,7 +1684,7 @@ namespace GridTools std::advance(cell, neighbor_permutation[i]); const Point p_unit = mapping.transform_real_to_unit_cell(*cell, p); - if (GeometryInfo::is_inside_unit_cell(p_unit)) + if (GeometryInfo::is_inside_unit_cell(p_unit, tolerance)) { cell_and_position.first = *cell; cell_and_position.second = p_unit; @@ -1727,10 +1727,8 @@ namespace GridTools // domain, or that we've had problems with the algorithm above. Try as a // last resort the other (simpler) algorithm. if (current_cell.state() != IteratorState::valid) - return find_active_cell_around_point(mapping, - mesh, - p, - marked_vertices); + return find_active_cell_around_point( + mapping, mesh, p, marked_vertices, tolerance); current_cell = typename MeshType::active_cell_iterator(); } @@ -5225,7 +5223,8 @@ namespace GridTools const Point & p, const typename Triangulation::active_cell_iterator & cell_hint, - const std::vector &marked_vertices) + const std::vector &marked_vertices, + const double tolerance) { const auto &mesh = cache.get_triangulation(); const auto &mapping = cache.get_mapping(); @@ -5241,7 +5240,8 @@ namespace GridTools vertex_to_cell_centers, cell_hint, marked_vertices, - used_vertices_rtree); + used_vertices_rtree, + tolerance); } template diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in index e08f457d1c..fb09c045f6 100644 --- a/source/grid/grid_tools.inst.in +++ b/source/grid/grid_tools.inst.in @@ -37,7 +37,8 @@ for (X : TRIANGULATIONS; deal_II_dimension : DIMENSIONS; deal_II_space_dimension, X>::type &, const std::vector &, - const RTree, unsigned int>> &); + const RTree, unsigned int>> &, + const double); template std::vector> compute_mesh_predicate_bounding_box( @@ -92,7 +93,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) const typename Triangulation< deal_II_dimension, deal_II_space_dimension>::active_cell_iterator &, - const std::vector &); + const std::vector &, + const double); template std::tuple class MeshType, int spacedim> -#ifndef _MSC_VER - std::pair::active_cell_iterator, - Point> -#else - std::pair< - typename dealii::internal:: - ActiveCellIterator>::type, - Point> -#endif - find_active_cell_around_point_tolerance( - const Mapping & mapping, - const MeshType &mesh, - const Point & p, - const std::vector & marked_vertices, - const double tolerance) - { - using active_cell_iterator = typename dealii::internal:: - ActiveCellIterator>::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> best_cell; - - // Find closest vertex and determine - // all adjacent cells - std::vector 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 adjacent_cells(adjacent_cells_tmp.begin(), - adjacent_cells_tmp.end()); - std::set 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::const_iterator - cell = adjacent_cells.begin(), - endc = adjacent_cells.end(); - for (; cell != endc; ++cell) - { - if ((*cell)->is_artificial() == false) - { - try - { - const Point 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::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:: - 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( - mesh, searched_cells, adjacent_cells); - } - } - - AssertThrow(best_cell.first.state() == IteratorState::valid, - ExcPointNotFound(p)); - - return best_cell; - } } // namespace @@ -550,10 +415,15 @@ namespace GridTools #endif find_active_cell_around_point(const MeshType &mesh, const Point & p, - const std::vector & marked_vertices) + const std::vector & marked_vertices, + const double tolerance) { return find_active_cell_around_point( - StaticMappingQ1::mapping, mesh, p, marked_vertices) + StaticMappingQ1::mapping, + mesh, + p, + marked_vertices, + tolerance) .first; } @@ -570,10 +440,124 @@ namespace GridTools find_active_cell_around_point(const Mapping & mapping, const MeshType &mesh, const Point & p, - const std::vector & marked_vertices) + const std::vector & 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>::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> best_cell; + + // Find closest vertex and determine + // all adjacent cells + std::vector 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 adjacent_cells(adjacent_cells_tmp.begin(), + adjacent_cells_tmp.end()); + std::set 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::const_iterator + cell = adjacent_cells.begin(), + endc = adjacent_cells.end(); + for (; cell != endc; ++cell) + { + if ((*cell)->is_artificial() == false) + { + try + { + const Point 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::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::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( + mesh, searched_cells, adjacent_cells); + } + } + + AssertThrow(best_cell.first.state() == IteratorState::valid, + ExcPointNotFound(p)); + + return best_cell; } @@ -596,7 +580,7 @@ namespace GridTools { 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( @@ -1285,7 +1269,8 @@ namespace GridTools find_active_cell_around_point( const hp::MappingCollection &mapping, const hp::DoFHandler & mesh, - const Point & p) + const Point & p, + const double tolerance) { Assert((mapping.size() == 1) || (mapping.size() == mesh.get_fe_collection().size()), @@ -1301,14 +1286,17 @@ namespace GridTools // 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 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; diff --git a/source/grid/grid_tools_dof_handlers.inst.in b/source/grid/grid_tools_dof_handlers.inst.in index 30ca321b91..eadd159d05 100644 --- a/source/grid/grid_tools_dof_handlers.inst.in +++ b/source/grid/grid_tools_dof_handlers.inst.in @@ -43,7 +43,8 @@ for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS; ActiveCellIterator::type find_active_cell_around_point(const X &, const Point &, - const std::vector &); + const std::vector &, + const double); template std::pair< dealii::internal::ActiveCellIterator &, const X &, const Point &, - const std::vector &); + const std::vector &, + const double); template std::vector< std::pair &, const hp::DoFHandler &, - const Point &); + const Point &, + const double); \} #endif