From: Luca Heltai Date: Tue, 19 Sep 2017 14:37:50 +0000 (+0200) Subject: GridTools::find_closest_vertex with mapping argument. X-Git-Tag: v9.0.0-rc1~1018^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=52e6dd9564fbf9aa70dd18ec46638d4d84ca502e;p=dealii.git GridTools::find_closest_vertex with mapping argument. --- diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index f6fdc81b45..9586e715c8 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -601,6 +601,8 @@ namespace GridTools * extract_used_vertices() * @param p The target point * @return The index of the vertex that is closest to the target point `p` + * + * @author Luca Heltai, 2017. */ template unsigned int @@ -633,6 +635,35 @@ namespace GridTools const Point &p, const std::vector &marked_vertices = std::vector()); + /** + * Find and return the number of the used vertex (or marked vertex) in a + * given mesh that is located closest to a given point. Use the given + * mapping to compute the actual location of the vertices. + * + * @param mapping A mapping used to compute the vertex locations + * @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 Luca Heltai, 2017 + */ + template class MeshType, int spacedim> + unsigned int + find_closest_vertex (const Mapping &mapping, + const MeshType &mesh, + const Point &p, + const std::vector &marked_vertices = std::vector()); + + /** * Find and return a vector of iterators to active cells that surround a * given vertex with index @p vertex_index. diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index cac93c628a..612d039126 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -1048,6 +1048,7 @@ namespace GridTools } + template class MeshType, int spacedim> unsigned int find_closest_vertex (const MeshType &mesh, @@ -1122,6 +1123,67 @@ namespace GridTools } + + template class MeshType, int spacedim> + unsigned int + find_closest_vertex (const Mapping &mapping, + const MeshType &mesh, + const Point &p, + const std::vector &marked_vertices) + { + // Take a shortcut in the simple case. + if (mapping.preserves_vertex_locations() == true) + return find_closest_vertex(mesh, p, marked_vertices); + + // first get the underlying + // triangulation from the + // mesh and determine vertices + // and used vertices + const Triangulation &tria = mesh.get_triangulation(); + + auto vertices = extract_used_vertices(tria, mapping); + + 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!") ); + + // Remove from the map unwanted elements. + if (marked_vertices.size()) + for (auto it = vertices.begin(); it != vertices.end(); ) + { + if (marked_vertices[it->first] == false) + { + vertices.erase(it++); + } + else + { + ++it; + } + } + + return find_closest_vertex(vertices, p); + } + + + template class MeshType, int spacedim> #ifndef _MSC_VER std::vector::active_cell_iterator> @@ -1362,7 +1424,7 @@ next_cell: // all adjacent cells std::vector adjacent_cells_tmp = find_cells_adjacent_to_vertex(mesh, - find_closest_vertex(mesh, p, marked_vertices)); + find_closest_vertex(mapping, mesh, p, marked_vertices)); // Make sure that we have found // at least one cell adjacent to vertex. diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in index 3938a71118..7c47ceaf4e 100644 --- a/source/grid/grid_tools.inst.in +++ b/source/grid/grid_tools.inst.in @@ -45,6 +45,13 @@ for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II const Point &, const std::vector &); + template + unsigned int + find_closest_vertex (const Mapping&, + const X &, + const Point &, + const std::vector &); + template std::vector::type> find_cells_adjacent_to_vertex(const X &, const unsigned int); diff --git a/tests/grid/grid_tools_find_closest_vertex_01.cc b/tests/grid/grid_tools_find_closest_vertex_01.cc new file mode 100644 index 0000000000..01fcac645f --- /dev/null +++ b/tests/grid/grid_tools_find_closest_vertex_01.cc @@ -0,0 +1,76 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2001 - 2015 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// Check extract used_vertices and find_closest_vertex when a mapping +// is present + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include + + +template +void test (const Point &p) +{ + deallog << "dim: " << dim << ", spacedim: " + << spacedim << std::endl; + + Triangulation tria; + GridGenerator::hyper_cube(tria); + tria.refine_global(1); + + // Vector fe + FESystem fe(FE_Q(1)^spacedim); + DoFHandler dh(tria); + dh.distribute_dofs(fe); + + // Displacement vector. + Vector disp(dh.n_dofs()); + for (unsigned int i=0; i, spacedim> mapping(1, dh, disp); + + auto m = GridTools::extract_used_vertices(tria, mapping); + + std::vector selected_vertices(tria.n_vertices(), false); + selected_vertices[0] = true; + selected_vertices[5] = true; + + for (auto &e: m) + deallog << "Vertex: " << e.first << ": " << e.second << std::endl; + + auto i = GridTools::find_closest_vertex(mapping,tria,p,selected_vertices); + deallog << "Closest vertex to " << p + << ", v["<< i << "] :" + << m[i] << std::endl; +}; + + +int main () +{ + initlog(); + test<2,2> (Point<2>(.2,.2)); + test<2,2> (Point<2>(.6,.9)); + return 0; +} + diff --git a/tests/grid/grid_tools_find_closest_vertex_01.output b/tests/grid/grid_tools_find_closest_vertex_01.output new file mode 100644 index 0000000000..0331bd7a3a --- /dev/null +++ b/tests/grid/grid_tools_find_closest_vertex_01.output @@ -0,0 +1,23 @@ + +DEAL::dim: 2, spacedim: 2 +DEAL::Vertex: 0: 0.500000 0.500000 +DEAL::Vertex: 1: 1.50000 0.500000 +DEAL::Vertex: 2: 0.500000 1.50000 +DEAL::Vertex: 3: 1.50000 1.50000 +DEAL::Vertex: 4: 1.00000 0.500000 +DEAL::Vertex: 5: 0.500000 1.00000 +DEAL::Vertex: 6: 1.50000 1.00000 +DEAL::Vertex: 7: 1.00000 1.50000 +DEAL::Vertex: 8: 1.00000 1.00000 +DEAL::Closest vertex to 0.200000 0.200000, v[0] :0.500000 0.500000 +DEAL::dim: 2, spacedim: 2 +DEAL::Vertex: 0: 0.500000 0.500000 +DEAL::Vertex: 1: 1.50000 0.500000 +DEAL::Vertex: 2: 0.500000 1.50000 +DEAL::Vertex: 3: 1.50000 1.50000 +DEAL::Vertex: 4: 1.00000 0.500000 +DEAL::Vertex: 5: 0.500000 1.00000 +DEAL::Vertex: 6: 1.50000 1.00000 +DEAL::Vertex: 7: 1.00000 1.50000 +DEAL::Vertex: 8: 1.00000 1.00000 +DEAL::Closest vertex to 0.600000 0.900000, v[5] :0.500000 1.00000