From bebba3a140d4e10f180a7fb1b9a95a589746f073 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Wed, 18 Oct 2017 19:39:53 +0200 Subject: [PATCH] Find active cell around point with Cache. --- include/deal.II/grid/grid_tools.h | 16 ++++ source/grid/grid_tools.cc | 19 +++++ source/grid/grid_tools.inst.in | 11 ++- .../grid/find_active_cell_around_point_02.cc | 83 +++++++++++++++++++ .../find_active_cell_around_point_02.output | 10 +++ 5 files changed, 138 insertions(+), 1 deletion(-) create mode 100644 tests/grid/find_active_cell_around_point_02.cc create mode 100644 tests/grid/find_active_cell_around_point_02.output diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index 715096c463..06c09e832b 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -92,6 +92,9 @@ namespace internal */ namespace GridTools { + template + class Cache; + /** * @name Information about meshes and cells */ @@ -926,6 +929,19 @@ namespace GridTools const hp::DoFHandler &mesh, const Point &p); + /** + * A version of the previous function that exploits an already existing + * GridTools::Cache object. + * + * @author Luca Heltai, 2017 + */ + template + std::pair::active_cell_iterator, Point > + find_active_cell_around_point (const Cache &cache, + const Point &p, + const typename Triangulation::active_cell_iterator &cell_hint=typename Triangulation::active_cell_iterator(), + const std::vector &marked_vertices = std::vector()); + /** * Return a list of all descendants of the given cell that are active. For * example, if the current cell is once refined but none of its children are diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 95664a9fc8..196abb7329 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -35,6 +35,7 @@ #include #include #include +#include #include #include #include @@ -4946,6 +4947,24 @@ next_cell: return id_and_v->first; } + template + std::pair::active_cell_iterator, Point > + find_active_cell_around_point(const Cache &cache, + const Point &p, + const typename Triangulation::active_cell_iterator &cell_hint, + const std::vector &marked_vertices) + { + const auto &mesh = cache.get_triangulation(); + const auto &mapping = cache.get_mapping(); + const auto &vertex_to_cells = cache.get_vertex_to_cell_map(); + const auto &vertex_to_cell_centers = cache.get_vertex_to_cell_centers_directions(); + + return find_active_cell_around_point(mapping, mesh, p, + vertex_to_cells, + vertex_to_cell_centers, + cell_hint, + marked_vertices); + } } /* namespace GridTools */ diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in index 2e21e8cf52..27020ca365 100644 --- a/source/grid/grid_tools.inst.in +++ b/source/grid/grid_tools.inst.in @@ -131,7 +131,16 @@ for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS extract_used_vertices(const Triangulation&mesh, const Mapping &mapping); - \} + template + std::pair::active_cell_iterator, + Point > + find_active_cell_around_point(const Cache&, + const Point &, + const typename Triangulation::active_cell_iterator &, + const std::vector &); + + + \} #endif } diff --git a/tests/grid/find_active_cell_around_point_02.cc b/tests/grid/find_active_cell_around_point_02.cc new file mode 100644 index 0000000000..cd371b229a --- /dev/null +++ b/tests/grid/find_active_cell_around_point_02.cc @@ -0,0 +1,83 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + + +/* + * Test cached version of find active cell around point + */ +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +template +void 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; + + GridTools::Cache cache(tria, mapping); + + auto cell = tria.begin_active(); + for (auto &p : points) + { + auto c_and_p = GridTools::find_active_cell_around_point(cache, p); + 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_02.output b/tests/grid/find_active_cell_around_point_02.output new file mode 100644 index 0000000000..02007871f7 --- /dev/null +++ b/tests/grid/find_active_cell_around_point_02.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 -- 2.39.5