]> https://gitweb.dealii.org/ - dealii.git/commitdiff
refactoring of GridTools::find_all_active_cells_around_point() 10477/head
authorNiklas Fehn <fehn@lnm.mw.tum.de>
Sun, 7 Jun 2020 17:22:33 +0000 (19:22 +0200)
committerNiklas Fehn <fehn@lnm.mw.tum.de>
Tue, 9 Jun 2020 08:34:38 +0000 (10:34 +0200)
include/deal.II/grid/grid_tools.h
source/grid/grid_tools_dof_handlers.cc
source/grid/grid_tools_dof_handlers.inst.in
tests/grid/find_all_active_cells_around_point_01_b.cc [new file with mode: 0644]
tests/grid/find_all_active_cells_around_point_01_b.output [new file with mode: 0644]

index a24b8c9a965aec85c339f4fd28b1bb4d9a68e625..1a19ca6c1bf064a4ad12e5aaf7b155a5a92bf7b5 100644 (file)
@@ -1280,19 +1280,48 @@ namespace GridTools
     const std::vector<bool> &marked_vertices = {});
 
   /**
-   * A variant of the previous find_active_cell_around_point() function that,
-   * instead of returning only the first matching cell, identifies all cells
+   * As compared to the functions above, this function identifies all cells
    * around a point for a given tolerance level `tolerance` in terms of unit
-   * coordinates. More precisely, whenever the point returned by
-   * find_active_cell_around_point() is within the given tolerance from the
-   * surface of the unit cell, all corresponding neighbors are also
-   * identified, including the location of the point in unit coordinates on
-   * any of these cells.
+   * coordinates. Given a first cell with reference coordinates as parameter
+   * @p first_cell, e.g. obtained by one of the functions above, all
+   * corresponding neighboring cells with points in unit coordinates are also
+   * identified.
    *
    * This function is useful e.g. for discontinuous function spaces where, for
-   * the case the given point `p` coincides with a vertex or an edge, several
+   * the case the given point `p` lies on a vertex, edge or face, several
    * cells might hold independent values of the solution that get combined in
    * some way in a user code.
+   *
+   * This function is used as follows
+   * @code
+   *   auto first_cell = GridTools::find_active_cell_around_point(...);
+   *   auto all_cells = GridTools::find_all_active_cells_around_point(mapping,
+   * mesh, p, tolerance, first_cell);
+   * @endcode
+   */
+  template <int dim, template <int, int> class MeshType, int spacedim>
+#  ifndef _MSC_VER
+  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
+                        Point<dim>>>
+#  else
+  std::vector<std::pair<
+    typename dealii::internal::
+      ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
+    Point<dim>>>
+#  endif
+  find_all_active_cells_around_point(
+    const Mapping<dim, spacedim> & mapping,
+    const MeshType<dim, spacedim> &mesh,
+    const Point<spacedim> &        p,
+    const double                   tolerance,
+    const std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
+                    Point<dim>> &  first_cell);
+
+  /**
+   * A variant of the previous function that internally calls one of the
+   * functions find_active_cell_around_point() to obtain a first cell, and
+   * subsequently adds all other cells by calling the function
+   * find_all_active_cells_around_point() above.
    */
   template <int dim, template <int, int> class MeshType, int spacedim>
 #  ifndef _MSC_VER
index 059bf89b50c50a9cdacff79cec6de31804146cbf..c2bdd9cadf9a48834000196d3491ae5a197162bd 100644 (file)
@@ -592,133 +592,156 @@ namespace GridTools
                                      const double                   tolerance,
                                      const std::vector<bool> &marked_vertices)
   {
-    // first use the result of the single point function as a guess. In order
-    // not to make the other find_all_active_cells_around_point more expensive
-    // and avoid some additional logic there, we first start with one cell as
-    // given by that other function (that possibly goes through a larger set
-    // of cells) and later add a list of more cells as appropriate.
-    std::vector<
-      std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
-                Point<dim>>>
-      cells_and_points;
     try
       {
-        cells_and_points.push_back(find_active_cell_around_point_tolerance(
-          mapping, mesh, p, marked_vertices, tolerance));
+        const auto cell_and_point = find_active_cell_around_point_tolerance(
+          mapping, mesh, p, marked_vertices, tolerance);
+
+        return find_all_active_cells_around_point(
+          mapping, mesh, p, tolerance, cell_and_point);
       }
     catch (ExcPointNotFound<spacedim> &)
       {}
 
-    if (!cells_and_points.empty())
+    return {};
+  }
+
+
+
+  template <int dim, template <int, int> class MeshType, int spacedim>
+#ifndef _MSC_VER
+  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
+                        Point<dim>>>
+#else
+  std::vector<std::pair<
+    typename dealii::internal::
+      ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
+    Point<dim>>>
+#endif
+  find_all_active_cells_around_point(
+    const Mapping<dim, spacedim> & mapping,
+    const MeshType<dim, spacedim> &mesh,
+    const Point<spacedim> &        p,
+    const double                   tolerance,
+    const std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
+                    Point<dim>> &  first_cell)
+  {
+    std::vector<
+      std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
+                Point<dim>>>
+      cells_and_points;
+
+    // insert the fist cell and point into the vector
+    cells_and_points.push_back(first_cell);
+
+    // check if the given point is on the surface of the unit cell. if yes,
+    // need to find all neighbors
+    const Point<dim> unit_point = cells_and_points.front().second;
+    const auto       my_cell    = cells_and_points.front().first;
+    Tensor<1, dim>   distance_to_center;
+    unsigned int     n_dirs_at_threshold     = 0;
+    unsigned int     last_point_at_threshold = numbers::invalid_unsigned_int;
+    for (unsigned int d = 0; d < dim; ++d)
       {
-        // check if the given point is on the surface of the unit cell. if yes,
-        // need to find all neighbors
-        const Point<dim> unit_point = cells_and_points.front().second;
-        const auto       my_cell    = cells_and_points.front().first;
-        Tensor<1, dim>   distance_to_center;
-        unsigned int     n_dirs_at_threshold = 0;
-        unsigned int last_point_at_threshold = numbers::invalid_unsigned_int;
-        for (unsigned int d = 0; d < dim; ++d)
+        distance_to_center[d] = std::abs(unit_point[d] - 0.5);
+        if (distance_to_center[d] > 0.5 - tolerance)
           {
-            distance_to_center[d] = std::abs(unit_point[d] - 0.5);
-            if (distance_to_center[d] > 0.5 - tolerance)
-              {
-                ++n_dirs_at_threshold;
-                last_point_at_threshold = d;
-              }
+            ++n_dirs_at_threshold;
+            last_point_at_threshold = d;
           }
+      }
 
+    std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
+      cells_to_add;
+    // point is within face -> only need neighbor
+    if (n_dirs_at_threshold == 1)
+      {
+        unsigned int neighbor_index =
+          2 * last_point_at_threshold +
+          (unit_point[last_point_at_threshold] > 0.5 ? 1 : 0);
+        if (!my_cell->at_boundary(neighbor_index))
+          cells_to_add.push_back(my_cell->neighbor(neighbor_index));
+      }
+    // corner point -> use all neighbors
+    else if (n_dirs_at_threshold == dim)
+      {
+        unsigned int local_vertex_index = 0;
+        for (unsigned int d = 0; d < dim; ++d)
+          local_vertex_index += (unit_point[d] > 0.5 ? 1 : 0) << d;
         std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
-          cells_to_add;
-        // point is within face -> only need neighbor
-        if (n_dirs_at_threshold == 1)
-          {
-            unsigned int neighbor_index =
-              2 * last_point_at_threshold +
-              (unit_point[last_point_at_threshold] > 0.5 ? 1 : 0);
-            if (!my_cell->at_boundary(neighbor_index))
-              cells_to_add.push_back(my_cell->neighbor(neighbor_index));
-          }
-        // corner point -> use all neighbors
-        else if (n_dirs_at_threshold == dim)
-          {
-            unsigned int local_vertex_index = 0;
-            for (unsigned int d = 0; d < dim; ++d)
-              local_vertex_index += (unit_point[d] > 0.5 ? 1 : 0) << d;
-            std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
-              cells = find_cells_adjacent_to_vertex(
-                mesh, my_cell->vertex_index(local_vertex_index));
-            for (const auto &cell : cells)
-              if (cell != my_cell)
-                cells_to_add.push_back(cell);
-          }
-        // point on line in 3D: We cannot simply take the intersection between
-        // the two vertices of cels because of hanging nodes. So instead we
-        // list the vertices around both points and then select the
-        // appropriate cells according to the result of read_to_unit_cell
-        // below.
-        else if (n_dirs_at_threshold == 2)
+          cells = find_cells_adjacent_to_vertex(
+            mesh, my_cell->vertex_index(local_vertex_index));
+        for (const auto &cell : cells)
+          if (cell != my_cell)
+            cells_to_add.push_back(cell);
+      }
+    // point on line in 3D: We cannot simply take the intersection between
+    // the two vertices of cells because of hanging nodes. So instead we
+    // list the vertices around both points and then select the
+    // appropriate cells according to the result of read_to_unit_cell
+    // below.
+    else if (n_dirs_at_threshold == 2)
+      {
+        std::pair<unsigned int, unsigned int> vertex_indices[3];
+        unsigned int                          count_vertex_indices = 0;
+        unsigned int free_direction = numbers::invalid_unsigned_int;
+        for (unsigned int d = 0; d < dim; ++d)
           {
-            std::pair<unsigned int, unsigned int> vertex_indices[3];
-            unsigned int                          count_vertex_indices = 0;
-            unsigned int free_direction = numbers::invalid_unsigned_int;
-            for (unsigned int d = 0; d < dim; ++d)
+            if (distance_to_center[d] > 0.5 - tolerance)
               {
-                if (distance_to_center[d] > 0.5 - tolerance)
-                  {
-                    vertex_indices[count_vertex_indices].first = d;
-                    vertex_indices[count_vertex_indices].second =
-                      unit_point[d] > 0.5 ? 1 : 0;
-                    count_vertex_indices++;
-                  }
-                else
-                  free_direction = d;
+                vertex_indices[count_vertex_indices].first = d;
+                vertex_indices[count_vertex_indices].second =
+                  unit_point[d] > 0.5 ? 1 : 0;
+                count_vertex_indices++;
               }
+            else
+              free_direction = d;
+          }
 
-            AssertDimension(count_vertex_indices, 2);
-            Assert(free_direction != numbers::invalid_unsigned_int,
-                   ExcInternalError());
+        AssertDimension(count_vertex_indices, 2);
+        Assert(free_direction != numbers::invalid_unsigned_int,
+               ExcInternalError());
 
-            const unsigned int first_vertex =
-              (vertex_indices[0].second << vertex_indices[0].first) +
-              (vertex_indices[1].second << vertex_indices[1].first);
-            for (unsigned int d = 0; d < 2; ++d)
+        const unsigned int first_vertex =
+          (vertex_indices[0].second << vertex_indices[0].first) +
+          (vertex_indices[1].second << vertex_indices[1].first);
+        for (unsigned int d = 0; d < 2; ++d)
+          {
+            auto tentative_cells = find_cells_adjacent_to_vertex(
+              mesh,
+              my_cell->vertex_index(first_vertex + (d << free_direction)));
+            for (const auto &cell : tentative_cells)
               {
-                auto tentative_cells = find_cells_adjacent_to_vertex(
-                  mesh,
-                  my_cell->vertex_index(first_vertex + (d << free_direction)));
-                for (const auto &cell : tentative_cells)
-                  {
-                    bool cell_not_yet_present = true;
-                    for (const auto &other_cell : cells_to_add)
-                      if (cell == other_cell)
-                        {
-                          cell_not_yet_present = false;
-                          break;
-                        }
-                    if (cell_not_yet_present)
-                      cells_to_add.push_back(cell);
-                  }
+                bool cell_not_yet_present = true;
+                for (const auto &other_cell : cells_to_add)
+                  if (cell == other_cell)
+                    {
+                      cell_not_yet_present = false;
+                      break;
+                    }
+                if (cell_not_yet_present)
+                  cells_to_add.push_back(cell);
               }
           }
+      }
 
-        const double original_distance_to_unit_cell =
-          GeometryInfo<dim>::distance_to_unit_cell(unit_point);
-        for (const auto &cell : cells_to_add)
-          {
-            if (cell != my_cell)
-              try
-                {
-                  const Point<dim> p_unit =
-                    mapping.transform_real_to_unit_cell(cell, p);
-                  if (GeometryInfo<dim>::distance_to_unit_cell(p_unit) <
-                      original_distance_to_unit_cell + tolerance)
-                    cells_and_points.emplace_back(cell, p_unit);
-                }
-              catch (typename Mapping<dim>::ExcTransformationFailed &)
-                {}
-          }
+    const double original_distance_to_unit_cell =
+      GeometryInfo<dim>::distance_to_unit_cell(unit_point);
+    for (const auto &cell : cells_to_add)
+      {
+        if (cell != my_cell)
+          try
+            {
+              const Point<dim> p_unit =
+                mapping.transform_real_to_unit_cell(cell, p);
+              if (GeometryInfo<dim>::distance_to_unit_cell(p_unit) <
+                  original_distance_to_unit_cell + tolerance)
+                cells_and_points.emplace_back(cell, p_unit);
+            }
+          catch (typename Mapping<dim>::ExcTransformationFailed &)
+            {}
       }
+
     std::sort(
       cells_and_points.begin(),
       cells_and_points.end(),
index 400f0c2984773e57faeb565e1ddfd82cd2304269..30ca321b913cf9a13390237580211bf63f5d1244 100644 (file)
@@ -56,6 +56,19 @@ for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS;
         const Point<deal_II_space_dimension> &,
         const std::vector<bool> &);
 
+      template std::vector<
+        std::pair<dealii::internal::ActiveCellIterator<deal_II_dimension,
+                                                       deal_II_space_dimension,
+                                                       X>::type,
+                  Point<deal_II_dimension>>>
+      find_all_active_cells_around_point(
+        const Mapping<deal_II_dimension, deal_II_space_dimension> &,
+        const X &,
+        const Point<deal_II_space_dimension> &,
+        const double,
+        const std::pair<typename X::active_cell_iterator,
+                        Point<deal_II_dimension>> &);
+
       template std::vector<
         std::pair<dealii::internal::ActiveCellIterator<deal_II_dimension,
                                                        deal_II_space_dimension,
diff --git a/tests/grid/find_all_active_cells_around_point_01_b.cc b/tests/grid/find_all_active_cells_around_point_01_b.cc
new file mode 100644 (file)
index 0000000..0538662
--- /dev/null
@@ -0,0 +1,77 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 - 2018 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.
+//
+// ---------------------------------------------------------------------
+
+
+// find_all_active_cells_around_point for a simple cube mesh
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/fe/mapping_q_generic.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include "../tests.h"
+
+
+
+template <int dim, int spacedim>
+void
+print_result(const Mapping<dim, spacedim> &      mapping,
+             const Triangulation<dim, spacedim> &tria,
+             const Point<dim>                    p,
+             const double                        tolerance = 1.e-10)
+{
+  deallog << "Testing " << dim << "D with point " << p << " tolerance "
+          << tolerance << std::endl;
+  auto first_cell = GridTools::find_active_cell_around_point(mapping, tria, p);
+  auto c_p        = GridTools::find_all_active_cells_around_point(
+    mapping, tria, p, tolerance, first_cell);
+  for (auto i : c_p)
+    deallog << "Cell: " << i.first->id() << " unit point " << i.second
+            << std::endl;
+  deallog << std::endl;
+}
+
+
+
+template <int dim, int spacedim>
+void
+test(unsigned int n_ref)
+{
+  Triangulation<dim, spacedim> tria;
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(n_ref);
+  MappingQGeneric<dim, spacedim> mapping(3);
+
+  Point<dim> p;
+  {
+    for (unsigned int d = 0; d < dim; ++d)
+      p[d] = 0.5;
+    print_result(mapping, tria, p);
+  }
+}
+
+
+
+int
+main()
+{
+  initlog();
+  deallog << std::setprecision(10);
+
+  test<3, 3>(1);
+}
diff --git a/tests/grid/find_all_active_cells_around_point_01_b.output b/tests/grid/find_all_active_cells_around_point_01_b.output
new file mode 100644 (file)
index 0000000..70cfe14
--- /dev/null
@@ -0,0 +1,11 @@
+
+DEAL::Testing 3D with point 0.5000000000 0.5000000000 0.5000000000 tolerance 1.000000000e-10
+DEAL::Cell: 0_1:0 unit point 1.000000000 1.000000000 1.000000000
+DEAL::Cell: 0_1:1 unit point 0.000000000 1.000000000 1.000000000
+DEAL::Cell: 0_1:2 unit point 1.000000000 0.000000000 1.000000000
+DEAL::Cell: 0_1:3 unit point 0.000000000 0.000000000 1.000000000
+DEAL::Cell: 0_1:4 unit point 1.000000000 1.000000000 0.000000000
+DEAL::Cell: 0_1:5 unit point 0.000000000 1.000000000 0.000000000
+DEAL::Cell: 0_1:6 unit point 1.000000000 0.000000000 0.000000000
+DEAL::Cell: 0_1:7 unit point 0.000000000 0.000000000 0.000000000
+DEAL::

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.