From: Luca Heltai Date: Sat, 12 Aug 2017 21:53:32 +0000 (-0600) Subject: New find active cell around point. X-Git-Tag: v9.0.0-rc1~1252^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=bb0a926f4e7e3db180d080a78c3cad3409b44f5e;p=dealii.git New find active cell around point. --- diff --git a/cmake/config/template-arguments.in b/cmake/config/template-arguments.in index f6a6523a48..a07bc908c2 100644 --- a/cmake/config/template-arguments.in +++ b/cmake/config/template-arguments.in @@ -156,6 +156,12 @@ TRIANGULATION_AND_DOFHANDLERS := { Triangulation; hp::DoFHandler } +// concrete Triangulation types (all types you can iterate over cells from) +// with hard-coded +TRIANGULATIONS := { Triangulation; + parallel::shared::Triangulation; + parallel::distributed::Triangulation;} + // serial and parallel sparsity pattern types SPARSITY_PATTERNS := { SparsityPattern; DynamicSparsityPattern; diff --git a/doc/news/changes/minor/20170813LucaHeltai b/doc/news/changes/minor/20170813LucaHeltai new file mode 100644 index 0000000000..4485d8253d --- /dev/null +++ b/doc/news/changes/minor/20170813LucaHeltai @@ -0,0 +1,5 @@ +New: A new version of GridTools::find_active_cell_around_point() +has been added that exploits local maps constructed +using standard GridTools utilities. +
+(Rene Gassmoeller, Luca Heltai, 2017/08/13) diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index 4e16025de7..77cc9e5153 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -747,6 +747,28 @@ namespace GridTools const Point &p, const std::vector &marked_vertices = std::vector()); + /** + * A version of the previous function that exploits an already existing + * map between vertices and cells, constructed using the function + * GridTools::vertex_to_cell_map, a map of vertex_to_cell_centers, obtained + * through GridTools::vertex_to_cells_center_directions, and a guess `cell_hint`. + * + * @author Luca Heltai, Rene Gassmoeller, 2017 + */ + template class MeshType, int spacedim> +#ifndef _MSC_VER + std::pair::active_cell_iterator, Point > +#else + std::pair >::type, Point > +#endif + find_active_cell_around_point (const Mapping &mapping, + const MeshType &mesh, + const Point &p, + const std::vector::active_cell_iterator > > &vertex_to_cell_map, + const std::vector > > &vertex_to_cell_centers, + const typename MeshType::active_cell_iterator &cell_hint=typename MeshType::active_cell_iterator(), + const std::vector &marked_vertices = std::vector()); + /** * A version of the previous function where we use that mapping on a given * cell that corresponds to the active finite element index of that cell. @@ -1024,10 +1046,16 @@ namespace GridTools vertex_to_cell_map(const Triangulation &triangulation); /** - * Returns a vector that contains a tensor for every vertex-cell combination - * of the output of GridTools::vertex_to_cell_map() (which is expected as - * input parameter for this function). Each tensor represents a geometric - * vector from the vertex to the respective cell center. + * Returns a vector of normalized tensors for each vertex-cell combination of + * the output of GridTools::vertex_to_cell_map() (which is expected as input + * parameter for this function). Each tensor represents a geometric vector + * from the vertex to the respective cell center. + * + * An assertion will be thrown if the size of the input vector is not equal to + * the number of vertices of the triangulation. + * + * result[v][c] is a unit Tensor for vertex index v, indicating the direction of + * the center of the c-th cell with respect to the vertex v. * * @author Rene Gassmoeller, Luca Heltai, 2017. */ @@ -1045,8 +1073,8 @@ namespace GridTools */ template unsigned int - get_closest_vertex_of_cell(const typename Triangulation::active_cell_iterator &cell, - const Point &position); + find_closest_vertex_of_cell(const typename Triangulation::active_cell_iterator &cell, + const Point &position); /** * Compute a globally unique index for each vertex and hanging node diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 920fd0f783..91728bb39b 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -1466,6 +1466,9 @@ next_cell: const std::vector > &vertices = mesh.get_vertices(); const unsigned int n_vertices = vertex_to_cells.size(); + AssertDimension(vertices.size(), n_vertices); + + std::vector > > vertex_to_cell_centers(n_vertices); for (unsigned int vertex=0; vertex + bool + compare_point_association(const unsigned int a, + const unsigned int b, + const Tensor<1,spacedim> &point_direction, + const std::vector > ¢er_directions) + { + const double scalar_product_a = center_directions[a] * point_direction; + const double scalar_product_b = center_directions[b] * point_direction; + + // The function is supposed to return if a is before b. We are looking + // for the alignment of point direction and center direction, therefore + // return if the scalar product of a is larger. + return (scalar_product_a > scalar_product_b); + } + } + + template class MeshType, int spacedim> +#ifndef _MSC_VER + std::pair::active_cell_iterator, Point > +#else + std::pair >::type, Point > +#endif + find_active_cell_around_point (const Mapping &mapping, + const MeshType &mesh, + const Point &p, + const std::vector::active_cell_iterator > > &vertex_to_cells, + const std::vector > > &vertex_to_cell_centers, + const typename MeshType::active_cell_iterator &cell_hint , + const std::vector &marked_vertices) + { + std::pair::active_cell_iterator, Point > cell_and_position; + + bool found_cell = false; + + unsigned int closest_vertex_index = 0; + Tensor<1,spacedim> vertex_to_point; + auto current_cell = cell_hint; + + while (found_cell == false) + { + // First look at the vertices of the cell cell_hint. If it's an + // invalid cell, then query for the closest global vertex + if (current_cell.state() == IteratorState::valid) + { + const unsigned int closest_vertex = find_closest_vertex_of_cell(current_cell , p); + vertex_to_point = p - current_cell ->vertex(closest_vertex); + closest_vertex_index = current_cell ->vertex_index(closest_vertex); + } + else + { + closest_vertex_index = GridTools::find_closest_vertex(mesh,p,marked_vertices); + vertex_to_point = p - mesh.get_vertices()[closest_vertex_index]; + } + + vertex_to_point /= vertex_to_point.norm(); + const unsigned int n_neighbor_cells = vertex_to_cells[closest_vertex_index].size(); + + // Create a corresponding map of vectors from vertex to cell center + std::vector neighbor_permutation(n_neighbor_cells); + + for (unsigned int i=0; i bool + { + return compare_point_association(a,b,vertex_to_point,vertex_to_cell_centers[closest_vertex_index]); + }; + + std::sort(neighbor_permutation.begin(), + neighbor_permutation.end(), + comp); + + // Search all of the cells adjacent to the closest vertex of the cell hint + // Most likely we will find the point in them. + for (unsigned int i=0; i p_unit = mapping.transform_real_to_unit_cell(*cell, p); + if (GeometryInfo::is_inside_unit_cell(p_unit)) + { + cell_and_position.first = *cell; + cell_and_position.second = p_unit; + found_cell = true; + break; + } + } + catch (typename Mapping::ExcTransformationFailed &) + {} + } + + if (found_cell == true) + return cell_and_position; + + // The first time around, we check for vertices in the hint_cell. If that + // does not work, we set the cell iterator to an invalid one, and look + // for a global vertex close to the point. If that does not work, we are in + // trouble, and just throw an exception. + // + // If we got here, then we did not find the point. If the + // current_cell.state() here is not IteratorState::valid, it means that + // the user did not provide a hint_cell, and at the beginning of the + // while loop we performed an actual global search on the mesh + // vertices. Not finding the point then means the point is outside the + // domain. + AssertThrow(current_cell.state() == IteratorState::valid, + ExcPointNotFound(p)); + + current_cell = typename MeshType::active_cell_iterator(); + } + return cell_and_position; + } + + template unsigned int - get_closest_vertex_of_cell(const typename Triangulation::active_cell_iterator &cell, - const Point &position) + find_closest_vertex_of_cell(const typename Triangulation::active_cell_iterator &cell, + const Point &position) { - double minimum_distance = std::numeric_limits::max(); - unsigned int closest_vertex = numbers::invalid_unsigned_int; + double minimum_distance = position.distance_square(cell->vertex(0)); + unsigned int closest_vertex = 0; - for (unsigned int v=0; v::vertices_per_cell; ++v) + for (unsigned int v=1; v::vertices_per_cell; ++v) { - const double vertex_distance = position.distance(cell->vertex(v)); + const double vertex_distance = position.distance_square(cell->vertex(v)); if (vertex_distance < minimum_distance) { closest_vertex = v; diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in index 5f97b7adb1..2b3f29b8e1 100644 --- a/source/grid/grid_tools.inst.in +++ b/source/grid/grid_tools.inst.in @@ -13,6 +13,24 @@ // // --------------------------------------------------------------------- +for (X : TRIANGULATIONS; deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS) +{ + +#if deal_II_dimension <= deal_II_space_dimension + namespace GridTools \{ + template + std::pair::type, Point > + find_active_cell_around_point (const Mapping &, + const X &, + const Point &, + const std::vector::type> > &, + const std::vector > > &, + const dealii::internal::ActiveCellIterator::type &, + const std::vector &); + \} + +#endif +} for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS) @@ -101,6 +119,8 @@ for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS template std::map compute_local_to_global_vertex_index_map(const parallel::distributed::Triangulation &triangulation); + + \} #endif @@ -241,6 +261,12 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS std::vector::active_cell_iterator> > vertex_to_cell_map(const Triangulation &triangulation); + template + std::vector > > + vertex_to_cell_centers_directions(const Triangulation &mesh, + const std::vector::active_cell_iterator> > &vertex_to_cells); + #if deal_II_dimension == deal_II_space_dimension # if deal_II_dimension > 1 template diff --git a/tests/grid/find_active_cell_around_point_01.cc b/tests/grid/find_active_cell_around_point_01.cc new file mode 100644 index 0000000000..e8cec9574f --- /dev/null +++ b/tests/grid/find_active_cell_around_point_01.cc @@ -0,0 +1,85 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 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. +// +// --------------------------------------------------------------------- + + +/* + * Given the number of refinements and the number of random points + * it benchmarks the time needed to run the function FCT + * which can be point_locator_D2 (or point_locator when it shall be written) + */ +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +template +double test(unsigned int n_ref, unsigned int n_points) +{ + deallog << "Testing " << dim << ", " << spacedim << std::endl; + + Triangulation tria; + GridGenerator::hyper_cube(tria); + tria.refine_global(n_ref); + DoFHandler dof_handler(tria); + + std::vector> points; + + deallog << "Points in study: " << n_points << std::endl; + for (size_t i=0; i p; + for (unsigned int d=0; d::mapping; + auto cell = tria.begin_active(); + for (auto &p : points) + { + auto c_and_p = GridTools::find_active_cell_around_point(mapping, tria, p, + v_to_c, v_to_c_d); + auto p2 = mapping.transform_unit_to_real_cell(c_and_p.first, c_and_p.second); + if (p2.distance(p)>1e-10) + deallog << "NOT OK!" << p << ", " << p2 + << ", " << c_and_p.first << std::endl; + cell = c_and_p.first; + } + deallog << "OK" << std::endl; +} + +int main () +{ + initlog(); + + test<1,1> (4,10); + test<2,2> (3,20); + test<3,3> (2,30); +} diff --git a/tests/grid/find_active_cell_around_point_01.output b/tests/grid/find_active_cell_around_point_01.output new file mode 100644 index 0000000000..02007871f7 --- /dev/null +++ b/tests/grid/find_active_cell_around_point_01.output @@ -0,0 +1,10 @@ + +DEAL::Testing 1, 1 +DEAL::Points in study: 10 +DEAL::OK +DEAL::Testing 2, 2 +DEAL::Points in study: 20 +DEAL::OK +DEAL::Testing 3, 3 +DEAL::Points in study: 30 +DEAL::OK