From: vishalkenchan Date: Wed, 26 Apr 2017 08:31:58 +0000 (+0200) Subject: extended find_closest_vertex() and find_active_cell_around_point() X-Git-Tag: v9.0.0-rc1~1647^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=22ca754e915afdfb0763f0c80a0b147956e53367;p=dealii.git extended find_closest_vertex() and find_active_cell_around_point() to take an optional custom mask for vertices. --- diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index 2282ad977f..469fb214d1 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -468,20 +468,30 @@ namespace GridTools /*@{*/ /** - * Find and return the number of the used vertex in a given mesh that is - * located closest to a given point. + * Find and return the number of the used vertex (or marked vertex) in a + * given mesh that is located closest to a given point. * * @param mesh A variable of a type that satisfies the requirements of the * @ref ConceptMeshType "MeshType concept". * @param p The point for which we want to find the closest vertex. + * @param marked_vertices An array of bools indicating which + * vertices of @p mesh will be considered within the search + * as the potentially closest vertex. On receiving a non-empty + * @p marked_vertices, the function will + * only search among @p marked_vertices for the closest vertex. + * The size of this array should be equal to the value returned by + * Triangulation::n_vertices() for the triangulation underlying the given mesh + * (as opposed to the value returned by Triangulation::n_used_vertices()). * @return The index of the closest vertex found. * + * * @author Ralf B. Schulz, 2006 */ template class MeshType, int spacedim> unsigned int find_closest_vertex (const MeshType &mesh, - const Point &p); + const Point &p, + const std::vector &marked_vertices = std::vector()); /** * Find and return a vector of iterators to active cells that surround a @@ -532,6 +542,14 @@ namespace GridTools * @param mesh A variable of a type that satisfies the requirements of the * @ref ConceptMeshType "MeshType concept". * @param p The point for which we want to find the surrounding cell. + * @param marked_vertices An array of bools indicating whether an + * entry in the vertex array should be considered + * (and the others must be ignored) as the potentially + * closest vertex to the specified point. On specifying a non-default + * @p marked_vertices, find_closest_vertex() would + * 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() ). * @return An iterator into the mesh that points to the surrounding cell. * * @note If the point requested does not lie in any of the cells of the mesh @@ -556,7 +574,8 @@ namespace GridTools typename dealii::internal::ActiveCellIterator >::type #endif find_active_cell_around_point (const MeshType &mesh, - const Point &p); + const Point &p, + const std::vector &marked_vertices = std::vector()); /** * Find and return an iterator to the active cell that surrounds a given @@ -580,6 +599,14 @@ namespace GridTools * @param mesh A variable of a type that satisfies the requirements of the * @ref ConceptMeshType "MeshType concept". * @param p The point for which we want to find the surrounding cell. + * @param marked_vertices An array of bools indicating whether an + * entry in the vertex array should be considered + * (and the others must be ignored) as the potentially + * closest vertex to the specified point. On specifying a non-default + * @p marked_vertices, find_closest_vertex() would + * 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() ). * @return An pair of an iterators into the mesh that points to the * surrounding cell, and of the coordinates of that point inside the cell in * the reference coordinates of that cell. This local position might be @@ -588,6 +615,19 @@ namespace GridTools * the unit cell, using GeometryInfo::project_to_unit_cell(). This is not * automatically performed by the algorithm. * + * @note When @p marked_vertices is specified the function should always be + * called inside a try block to catch the exception that the function might + * throw in the case it couldn't find an active cell surrounding the point. + * The motivation of using @p marked_vertices is to cut down the search space + * of vertices if one has a priori knowledge of a collection of vertices that + * the point of interest may be close to. For instance, in the case when a + * parallel::shared::Triangulation is employed and we are looking for a point + * that we know is inside the locally owned part of the mesh, then it would + * make sense to pass an array for @p marked_vertices that flags only the + * vertices of all locally owned active cells. If, however, the function + * throws an exception, then that would imply that the point lies outside + * locally owned active cells. + * * @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 @@ -611,7 +651,8 @@ namespace GridTools #endif find_active_cell_around_point (const Mapping &mapping, const MeshType &mesh, - const Point &p); + const Point &p, + const std::vector &marked_vertices = std::vector()); /** * A version of the previous function where we use that mapping on a given diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 6c0467c0c5..20b4021d3d 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -1047,11 +1047,11 @@ namespace GridTools } - template class MeshType, int spacedim> unsigned int find_closest_vertex (const MeshType &mesh, - const Point &p) + const Point &p, + const std::vector &marked_vertices) { // first get the underlying // triangulation from the @@ -1060,7 +1060,37 @@ namespace GridTools const Triangulation &tria = mesh.get_triangulation(); const std::vector< Point > &vertices = tria.get_vertices(); - const std::vector< bool > &used = tria.get_used_vertices(); + + Assert ( tria.get_vertices().size() == marked_vertices.size() || marked_vertices.size() ==0, + ExcDimensionMismatch(tria.get_vertices().size(), marked_vertices.size())); + + // If p is an element of marked_vertices, + // and q is that of used_Vertices, + // the vector marked_vertices does NOT + // contain unused vertices if p implies q. + // I.e., if p is true q must be true + // (if p is false, q could be false or true). + // p implies q logic is encapsulated in ~p|q. + Assert( marked_vertices.size()==0 + || + std::equal( marked_vertices.begin(), + marked_vertices.end(), + tria.get_used_vertices().begin(), + [](bool p, bool q) + { + return !p || q; + }), + ExcMessage("marked_vertices should be a subset of used vertices in the triangulation " + "but marked_vertices contains one or more vertices that are not used vertices!") ); + + // In addition, if a vector bools + // is specified (marked_vertices) + // marking all the vertices which + // could be the potentially closest + // vertex to the point, use it instead + // of used vertices + const std::vector &used = + (marked_vertices.size()==0) ? tria.get_used_vertices() : marked_vertices; // At the beginning, the first // used vertex is the closest one @@ -1296,12 +1326,13 @@ next_cell: typename dealii::internal::ActiveCellIterator >::type #endif find_active_cell_around_point (const MeshType &mesh, - const Point &p) + const Point &p, + const std::vector &marked_vertices) { return find_active_cell_around_point (StaticMappingQ1::mapping, - mesh, p).first; + mesh, p, marked_vertices).first; } @@ -1313,7 +1344,8 @@ next_cell: #endif find_active_cell_around_point (const Mapping &mapping, const MeshType &mesh, - const Point &p) + const Point &p, + const std::vector &marked_vertices) { typedef typename dealii::internal::ActiveCellIterator >::type active_cell_iterator; @@ -1329,7 +1361,7 @@ next_cell: // all adjacent cells std::vector adjacent_cells_tmp = find_cells_adjacent_to_vertex(mesh, - find_closest_vertex(mesh, p)); + find_closest_vertex(mesh, p, marked_vertices)); // Make sure that we have found // at least one cell adjacent to vertex. @@ -1399,6 +1431,12 @@ next_cell: // 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 diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in index 9d557cd056..8b30e95e2b 100644 --- a/source/grid/grid_tools.inst.in +++ b/source/grid/grid_tools.inst.in @@ -24,7 +24,8 @@ for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II template unsigned int find_closest_vertex (const X &, - const Point &); + const Point &, + const std::vector &); template std::vector::type> @@ -32,13 +33,14 @@ for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II template dealii::internal::ActiveCellIterator::type - find_active_cell_around_point (const X &, const Point &p); + find_active_cell_around_point (const X &, const Point &, const std::vector &); template std::pair::type, Point > find_active_cell_around_point (const Mapping &, const X &, - const Point &); + const Point &, + const std::vector &); template std::vector::type>