From: Johannes Heinz Date: Wed, 22 Jun 2022 12:08:01 +0000 (+0200) Subject: fix find_active_cell_around_point to marked cells or no cell X-Git-Tag: v9.5.0-rc1~1119^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F14059%2Fhead;p=dealii.git fix find_active_cell_around_point to marked cells or no cell Co-authored-by: Marc Fehling Co-authored-by: Peter Munch --- diff --git a/doc/news/changes/minor/20220628Heinz b/doc/news/changes/minor/20220628Heinz new file mode 100644 index 0000000000..9cea2675d9 --- /dev/null +++ b/doc/news/changes/minor/20220628Heinz @@ -0,0 +1,3 @@ +Fixed: The function find_active_cell_around_point() now skips cells without any marked vertices. +
+(Johannes Heinz, 2022/06/28) diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index d8c3279452..81cc38ff33 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -2850,15 +2850,51 @@ namespace GridTools bool found_cell = false; bool approx_cell = false; - unsigned int closest_vertex_index = 0; + unsigned int closest_vertex_index = 0; + // ensure closest vertex index is a marked one, otherwise cell (with vertex + // 0) might be found even though it is not marked. This is only relevant if + // searching with rtree, using find_closest_vertex already can manage not + // finding points + if (marked_vertices.size() && !used_vertices_rtree.empty()) + { + const auto itr = + std::find(marked_vertices.begin(), marked_vertices.end(), true); + Assert(itr != marked_vertices.end(), + dealii::ExcMessage("No vertex has been marked!")); + closest_vertex_index = std::distance(marked_vertices.begin(), itr); + } + Tensor<1, spacedim> vertex_to_point; auto current_cell = cell_hint; + // check whether cell has at least one marked vertex + const auto cell_marked = [&mesh, &marked_vertices](const auto &cell) { + if (marked_vertices.size() == 0) + return true; + + if (cell != mesh.active_cell_iterators().end()) + for (unsigned int i = 0; i < cell->n_vertices(); ++i) + if (marked_vertices[cell->vertex_index(i)]) + return true; + + return false; + }; + + // check whether any cell in collection is marked + const auto any_cell_marked = [&cell_marked](const auto &cells) { + return std::any_of(cells.begin(), + cells.end(), + [&cell_marked](const auto &cell) { + return cell_marked(cell); + }); + }; + 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) + if (current_cell.state() == IteratorState::valid && + cell_marked(cell_hint)) { const auto cell_vertices = mapping.get_vertices(current_cell); const unsigned int closest_vertex = @@ -2888,9 +2924,14 @@ namespace GridTools boost::geometry::index::satisfies(marked), std::back_inserter(res)); - // We should have one and only one result - AssertDimension(res.size(), 1); - closest_vertex_index = res[0].second; + // Searching for a point which is located outside the + // triangulation results in res.size() = 0 + Assert(res.size() < 2, + dealii::ExcMessage("There can not be multiple results")); + + if (res.size() > 0) + if (any_cell_marked(vertex_to_cells[res[0].second])) + closest_vertex_index = res[0].second; } else { @@ -2900,6 +2941,14 @@ namespace GridTools vertex_to_point = p - mesh.get_vertices()[closest_vertex_index]; } +#ifdef DEBUG + { + // Double-check if found index is at marked cell + Assert(any_cell_marked(vertex_to_cells[closest_vertex_index]), + dealii::ExcMessage("Found non-marked vertex")); + } +#endif + const double vertex_point_norm = vertex_to_point.norm(); if (vertex_point_norm > 0) vertex_to_point /= vertex_point_norm; diff --git a/tests/grid/grid_tools_07.cc b/tests/grid/grid_tools_07.cc new file mode 100644 index 0000000000..a40c3c3667 --- /dev/null +++ b/tests/grid/grid_tools_07.cc @@ -0,0 +1,190 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2001 - 2021 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +// +// check find_active_cell_around_point finds only marked cells. +// If there are no marked cell containing the point it should not find +// anything. Latter is especially important if run with MPI. +// + +#include +#include +#include +#include + +#include "../tests.h" + +/* + * Generate a non-matching grid consisting of two elements which + * do not share any face + */ + +void +generate_grid(Triangulation<2> &triangulation) +{ + // +-------------------+ +-------------------+ + // | | | | | + // | 2 | 3 | | | + // | | | | | + // 10|_________|_________|8 9| 4 | + // | | | | | + // | | | | | + // | 0 | 1 | p | | + // | | | | | + // +-------------------+ +-------------------+ + + Triangulation<2> triangulation0; + Triangulation<2> triangulation1; + GridGenerator::hyper_rectangle(triangulation0, {-.5, -.5}, {.5, .5}, true); + for (const auto &face : triangulation0.active_face_iterators()) + if (face->at_boundary() && face->boundary_id() == 1) + face->set_boundary_id(8); + for (const auto &face : triangulation0.active_face_iterators()) + if (face->at_boundary() && face->boundary_id() == 0) + face->set_boundary_id(10); + triangulation0.refine_global(1); + GridGenerator::flatten_triangulation(triangulation0, triangulation1); + + Triangulation<2> triangulation2; + GridGenerator::hyper_rectangle(triangulation2, {.5, -.5}, {1.5, .5}, true); + for (const auto &face : triangulation2.active_face_iterators()) + if (face->at_boundary() && face->boundary_id() == 0) + face->set_boundary_id(9); + + // tolerance 0. to ensure vertices are not merged + GridGenerator::merge_triangulations( + triangulation1, triangulation2, triangulation, 0., false); + + // make sure boundary ids are kept + std::vector boundary_ids; + for (const auto &face : triangulation1.active_face_iterators()) + if (face->at_boundary()) + boundary_ids.emplace_back(face->boundary_id()); + for (const auto &face : triangulation2.active_face_iterators()) + if (face->at_boundary()) + boundary_ids.emplace_back(face->boundary_id()); + + unsigned int i = 0; + for (const auto &face : triangulation.active_face_iterators()) + if (face->at_boundary()) + face->set_boundary_id(boundary_ids[i++]); +} + +std::vector +mark_vertices_at_boundary(const types::boundary_id boundary_id, + const Triangulation<2> & triangulation, + std::vector & marked_vertices) +{ + std::vector marked_cell_idxs; + + // mark_vertices_of_cells_at_boundary + for (const auto &cell : triangulation.active_cell_iterators()) + { + for (unsigned int face = 0; face < cell->n_faces(); ++face) + { + if (cell->face(face)->at_boundary() && + cell->face(face)->boundary_id() == boundary_id) + { + marked_cell_idxs.push_back(cell->index()); + for (unsigned int v = 0; v < cell->face(face)->n_vertices(); ++v) + { + marked_vertices[cell->face(face)->vertex_index(v)] = true; + } + continue; + } + } + } + return marked_cell_idxs; +} + +int +find_cell_at_point(const Point<2> & p, + const Triangulation<2> & triangulation, + const std::vector &marked_vertices, + const Triangulation<2>::active_cell_iterator &cell_hint, + const double tolerance = 1e-6) +{ + const MappingQ<2> mapping{1}; + const GridTools::Cache<2, 2> cache{triangulation, mapping}; + const auto cell_and_pnt = GridTools::find_active_cell_around_point( + cache.get_mapping(), + cache.get_triangulation(), + p, + cache.get_vertex_to_cell_map(), + cache.get_vertex_to_cell_centers_directions(), + cell_hint, + marked_vertices, + cache.get_used_vertices_rtree(), + tolerance, + &cache.get_locally_owned_cell_bounding_boxes_rtree()); + + if (cell_and_pnt.first != triangulation.active_cell_iterators().end()) + return cell_and_pnt.first->index(); + return -1; +} + +int +main() +{ + // Setup + Triangulation<2> triangulation; + generate_grid(triangulation); + + const Point<2> p{.5, -.25}; // point located between cell 4 and 1 + const auto cell_hint = ++triangulation.begin_active(); // hint cell 1 + + initlog(); + deallog << std::setprecision(4); + + // Test1 + deallog << "Test if cell not found:\n\n"; + deallog << "Searching for point with marked cells that do not hold the " + "point:\n"; + { + std::vector marked_vertices(triangulation.n_vertices(), false); + const auto cell_idxs = + mark_vertices_at_boundary(10, triangulation, marked_vertices); + const int found_cell = + find_cell_at_point(p, triangulation, marked_vertices, cell_hint); + + deallog << "Marked cells "; + for (const auto &cell_idx : cell_idxs) + deallog << cell_idx << " "; + deallog << "\nFound cell: " << found_cell << "\n"; + deallog << "\n---------------\n\n"; + } + + // Test2 + deallog << "Test if marked cell is found with cell hint that has the " + "point but is not marked:\n\n"; + { + std::vector marked_vertices(triangulation.n_vertices(), false); + + // mark cell 4 + const auto cell_idxs = + mark_vertices_at_boundary(9, triangulation, marked_vertices); + deallog << "Searching for cells "; + for (const auto &cell_idx : cell_idxs) + deallog << cell_idx << " "; + deallog << "with cell_hint " << cell_hint->index() << ":\n"; + + const int found_cell = + find_cell_at_point(p, triangulation, marked_vertices, cell_hint); + deallog << "Found cell: " << found_cell << "\n"; + } + deallog << std::endl << std::flush; + + return 0; +} diff --git a/tests/grid/grid_tools_07.output b/tests/grid/grid_tools_07.output new file mode 100644 index 0000000000..014804cc42 --- /dev/null +++ b/tests/grid/grid_tools_07.output @@ -0,0 +1,15 @@ + +DEAL::Test if cell not found: + +Searching for point with marked cells that do not hold the point: +Marked cells 0 2 +Found cell: -1 + +--------------- + +Test if marked cell is found with cell hint that has the point but is not marked: + +Searching for cells 4 with cell_hint 1: +Found cell: 4 + +DEAL:: \ No newline at end of file