From b31381847baaee0eebff10d7e789ce959c119f1a Mon Sep 17 00:00:00 2001 From: Giovanni Alzetta Date: Tue, 21 Nov 2017 15:56:57 +0100 Subject: [PATCH] Added tolerance and test for GridTools::find_active_cell_around_point and test --- .../changes/minor/20171121GiovanniAlzetta | 3 + source/grid/grid_tools.cc | 22 +++++ .../grid/find_active_cell_around_point_03.cc | 99 +++++++++++++++++++ .../find_active_cell_around_point_03.output | 13 +++ 4 files changed, 137 insertions(+) create mode 100644 doc/news/changes/minor/20171121GiovanniAlzetta create mode 100644 tests/grid/find_active_cell_around_point_03.cc create mode 100644 tests/grid/find_active_cell_around_point_03.output diff --git a/doc/news/changes/minor/20171121GiovanniAlzetta b/doc/news/changes/minor/20171121GiovanniAlzetta new file mode 100644 index 0000000000..7d32c09792 --- /dev/null +++ b/doc/news/changes/minor/20171121GiovanniAlzetta @@ -0,0 +1,3 @@ +Fixed: Added a tolerance inside GridTools::find_active_cell_around_point to errors with boundary points and a test with such points. +
+(Giovanni Alzetta, 2017/11/21) diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index cc6c40ff50..99f2e415e1 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -1601,8 +1601,11 @@ next_cell: const std::vector &marked_vertices) { std::pair::active_cell_iterator, Point > cell_and_position; + // To handle points at the border we keep track of points which are close to the unit cell: + std::pair::active_cell_iterator, Point > cell_and_position_approx; bool found_cell = false; + bool approx_cell = false; unsigned int closest_vertex_index = 0; Tensor<1,spacedim> vertex_to_point; @@ -1644,6 +1647,11 @@ next_cell: std::sort(neighbor_permutation.begin(), neighbor_permutation.end(), 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; // Search all of the cells adjacent to the closest vertex of the cell hint // Most likely we will find the point in them. @@ -1659,8 +1667,20 @@ next_cell: cell_and_position.first = *cell; cell_and_position.second = p_unit; found_cell = true; + approx_cell = false; break; } + // The point is not inside this cell: checking how far outside it is + // and whether we want to use this cell as a backup if we can't find + // a cell within which the point lies. + const double dist = GeometryInfo::distance_to_unit_cell(p_unit); + if (dist < best_distance) + { + best_distance = dist; + cell_and_position_approx.first = *cell; + cell_and_position_approx.second = p_unit; + approx_cell = true; + } } catch (typename Mapping::ExcTransformationFailed &) {} @@ -1668,6 +1688,8 @@ next_cell: if (found_cell == true) return cell_and_position; + else if (approx_cell == true) + return cell_and_position_approx; // 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 diff --git a/tests/grid/find_active_cell_around_point_03.cc b/tests/grid/find_active_cell_around_point_03.cc new file mode 100644 index 0000000000..789d0c55da --- /dev/null +++ b/tests/grid/find_active_cell_around_point_03.cc @@ -0,0 +1,99 @@ +// --------------------------------------------------------------------- +// +// 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 + +template +void test() +{ + deallog << "Testing " << dim << ", " << spacedim << std::endl; + // The hypercube is [a,b]^spacedim + const double a = -0.3; + const double b = 0.7; + + + Triangulation triangulation; + GridGenerator::hyper_cube(triangulation, a, b); + triangulation.refine_global(dim == 2 ? 2 : 1); + + MappingQ mapping(1); + + DoFHandler dof_handler(triangulation); + + GridTools::Cache cache(triangulation,mapping); + Point p1; + Point p2; + Point p3; + for (unsigned int d=0; d::active_cell_iterator, + Point< dim > > + res_1 = GridTools::find_active_cell_around_point + (mapping, dof_handler, p1); + std::pair< typename Triangulation::active_cell_iterator, + Point< dim > > + res_2 = GridTools::find_active_cell_around_point + (mapping, dof_handler, p2); + std::pair< typename Triangulation::active_cell_iterator, + Point< dim > > + res_3 = GridTools::find_active_cell_around_point + (mapping, dof_handler, p3); + + deallog << "Standard find active cell around point: OK " << std::endl; + + // Failing with cache + auto res_c_1 = GridTools::find_active_cell_around_point (cache, p1); + auto res_c_2 = GridTools::find_active_cell_around_point (cache, p2); + auto res_c_3 = GridTools::find_active_cell_around_point (cache, p3); + + deallog << "Cache find active cell around point: OK " << std::endl; + + if ( res_1.first != res_c_1.first) + deallog << "Different cells were found for p1" << std::endl; + else if ( res_2.first != res_c_2.first) + deallog << "Different cells were found for p2" << std::endl; + else if ( res_3.first != res_c_3.first ) + deallog << "Different cells were found for p3" << std::endl; + + deallog << "Test: OK " << std::endl; +} + +int main () +{ + initlog(); + test<1,1> (); + test<2,2> (); + test<3,3> (); +} diff --git a/tests/grid/find_active_cell_around_point_03.output b/tests/grid/find_active_cell_around_point_03.output new file mode 100644 index 0000000000..15b55669ab --- /dev/null +++ b/tests/grid/find_active_cell_around_point_03.output @@ -0,0 +1,13 @@ + +DEAL::Testing 1, 1 +DEAL::Standard find active cell around point: OK +DEAL::Cache find active cell around point: OK +DEAL::Test: OK +DEAL::Testing 2, 2 +DEAL::Standard find active cell around point: OK +DEAL::Cache find active cell around point: OK +DEAL::Test: OK +DEAL::Testing 3, 3 +DEAL::Standard find active cell around point: OK +DEAL::Cache find active cell around point: OK +DEAL::Test: OK -- 2.39.5