]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Added tolerance and test for GridTools::find_active_cell_around_point and test 5512/head
authorGiovanni Alzetta <giovannialzetta@hotmail.it>
Tue, 21 Nov 2017 14:56:57 +0000 (15:56 +0100)
committerGiovanni Alzetta <giovannialzetta@hotmail.it>
Sat, 25 Nov 2017 10:19:58 +0000 (11:19 +0100)
doc/news/changes/minor/20171121GiovanniAlzetta [new file with mode: 0644]
source/grid/grid_tools.cc
tests/grid/find_active_cell_around_point_03.cc [new file with mode: 0644]
tests/grid/find_active_cell_around_point_03.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20171121GiovanniAlzetta b/doc/news/changes/minor/20171121GiovanniAlzetta
new file mode 100644 (file)
index 0000000..7d32c09
--- /dev/null
@@ -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.
+<br>
+(Giovanni Alzetta, 2017/11/21)
index cc6c40ff50a90b4aac4c7a460bb71d0ff7de29f8..99f2e415e1ca7648f87bd488331bb4d2a4217a90 100644 (file)
@@ -1601,8 +1601,11 @@ next_cell:
                                  const std::vector<bool>                                                         &marked_vertices)
   {
     std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim> > cell_and_position;
+    // To handle points at the border we keep track of points which are close to the unit cell:
+    std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim> > 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<dim>::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<dim>::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 (file)
index 0000000..789d0c5
--- /dev/null
@@ -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 <deal.II/dofs/dof_handler.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/grid/grid_tools_cache.h>
+
+template <int dim, int spacedim>
+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<dim,spacedim> triangulation;
+  GridGenerator::hyper_cube(triangulation, a, b);
+  triangulation.refine_global(dim == 2 ? 2 : 1);
+
+  MappingQ<dim> mapping(1);
+
+  DoFHandler<spacedim> dof_handler(triangulation);
+
+  GridTools::Cache<dim,spacedim> cache(triangulation,mapping);
+  Point<spacedim> p1;
+  Point<spacedim> p2;
+  Point<spacedim> p3;
+  for (unsigned int d=0; d<spacedim; ++d)
+    {
+      p1[d] = a;
+      if (d==1)
+        p2[d] = b;
+      else
+        p2[d] = (a+b)*0.5;
+      p3[d] = b;
+    }
+
+  std::pair< typename Triangulation<dim,spacedim>::active_cell_iterator,
+      Point< dim > >
+      res_1  = GridTools::find_active_cell_around_point
+               (mapping, dof_handler, p1);
+  std::pair< typename Triangulation<dim,spacedim>::active_cell_iterator,
+      Point< dim > >
+      res_2  = GridTools::find_active_cell_around_point
+               (mapping, dof_handler, p2);
+  std::pair< typename Triangulation<dim,spacedim>::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 (file)
index 0000000..15b5566
--- /dev/null
@@ -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 

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.