]> https://gitweb.dealii.org/ - dealii.git/commitdiff
fix find_active_cell_around_point to marked cells or no cell 14059/head
authorJohannes Heinz <johannes.heinz@tuwien.ac.at>
Wed, 22 Jun 2022 12:08:01 +0000 (14:08 +0200)
committerJohannes Heinz <johannes.heinz@tuwien.ac.at>
Mon, 4 Jul 2022 10:26:00 +0000 (12:26 +0200)
Co-authored-by: Marc Fehling <mafehling.git@gmail.com>
Co-authored-by: Peter Munch <peterrmuench@gmail.com>
doc/news/changes/minor/20220628Heinz [new file with mode: 0644]
source/grid/grid_tools.cc
tests/grid/grid_tools_07.cc [new file with mode: 0644]
tests/grid/grid_tools_07.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20220628Heinz b/doc/news/changes/minor/20220628Heinz
new file mode 100644 (file)
index 0000000..9cea267
--- /dev/null
@@ -0,0 +1,3 @@
+Fixed: The function find_active_cell_around_point() now skips cells without any marked vertices.
+<br>
+(Johannes Heinz, 2022/06/28)
index d8c32794526653257a72360fb9acc6c31f8545c7..81cc38ff33b0ad8b0c7c77e9812abc22cf79e8ef 100644 (file)
@@ -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 (file)
index 0000000..a40c3c3
--- /dev/null
@@ -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 <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_tools_cache.h>
+#include <deal.II/grid/tria.h>
+
+#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<types::boundary_id> 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<unsigned int>
+mark_vertices_at_boundary(const types::boundary_id boundary_id,
+                          const Triangulation<2> & triangulation,
+                          std::vector<bool> &      marked_vertices)
+{
+  std::vector<unsigned int> 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<bool> &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<bool> 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<bool> 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 (file)
index 0000000..014804c
--- /dev/null
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.