From: Luca Heltai <luca.heltai@sissa.it>
Date: Wed, 18 Oct 2017 17:39:53 +0000 (+0200)
Subject: Find active cell around point with Cache.
X-Git-Tag: v9.0.0-rc1~932^2~1
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=bebba3a140d4e10f180a7fb1b9a95a589746f073;p=dealii.git

Find active cell around point with Cache.
---

diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h
index 715096c463..06c09e832b 100644
--- a/include/deal.II/grid/grid_tools.h
+++ b/include/deal.II/grid/grid_tools.h
@@ -92,6 +92,9 @@ namespace internal
  */
 namespace GridTools
 {
+  template <int dim, int spacedim>
+  class Cache;
+
   /**
    * @name Information about meshes and cells
    */
@@ -926,6 +929,19 @@ namespace GridTools
                                  const hp::DoFHandler<dim,spacedim>        &mesh,
                                  const Point<spacedim>                     &p);
 
+  /**
+   * A version of the previous function that exploits an already existing
+   * GridTools::Cache<dim,spacedim> object.
+   *
+   * @author Luca Heltai, 2017
+   */
+  template <int dim, int spacedim>
+  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator, Point<dim> >
+  find_active_cell_around_point (const Cache<dim,spacedim>                                         &cache,
+                                 const Point<spacedim>                                             &p,
+                                 const typename Triangulation<dim, spacedim>::active_cell_iterator &cell_hint=typename Triangulation<dim, spacedim>::active_cell_iterator(),
+                                 const std::vector<bool>                                           &marked_vertices = std::vector<bool>());
+
   /**
    * Return a list of all descendants of the given cell that are active. For
    * example, if the current cell is once refined but none of its children are
diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc
index 95664a9fc8..196abb7329 100644
--- a/source/grid/grid_tools.cc
+++ b/source/grid/grid_tools.cc
@@ -35,6 +35,7 @@
 #include <deal.II/grid/tria_boundary.h>
 #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/grid_reordering.h>
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/dofs/dof_accessor.h>
@@ -4946,6 +4947,24 @@ next_cell:
     return id_and_v->first;
   }
 
+  template<int dim, int spacedim>
+  std::pair<typename Triangulation<dim,spacedim>::active_cell_iterator, Point<dim> >
+  find_active_cell_around_point(const Cache<dim,spacedim> &cache,
+                                const Point<spacedim> &p,
+                                const typename Triangulation<dim,spacedim>::active_cell_iterator &cell_hint,
+                                const std::vector<bool> &marked_vertices)
+  {
+    const auto &mesh = cache.get_triangulation();
+    const auto &mapping = cache.get_mapping();
+    const auto &vertex_to_cells = cache.get_vertex_to_cell_map();
+    const auto &vertex_to_cell_centers = cache.get_vertex_to_cell_centers_directions();
+
+    return find_active_cell_around_point(mapping, mesh, p,
+                                         vertex_to_cells,
+                                         vertex_to_cell_centers,
+                                         cell_hint,
+                                         marked_vertices);
+  }
 } /* namespace GridTools */
 
 
diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in
index 2e21e8cf52..27020ca365 100644
--- a/source/grid/grid_tools.inst.in
+++ b/source/grid/grid_tools.inst.in
@@ -131,7 +131,16 @@ for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS
         extract_used_vertices(const Triangulation<deal_II_dimension,deal_II_space_dimension>&mesh,
                               const Mapping<deal_II_dimension,deal_II_space_dimension> &mapping);
 
-                                    \}
+        template
+        std::pair<typename Triangulation<deal_II_dimension,deal_II_space_dimension>::active_cell_iterator,
+                  Point<deal_II_dimension> >
+        find_active_cell_around_point(const Cache<deal_II_dimension,deal_II_space_dimension>&,
+                                      const Point<deal_II_space_dimension> &,
+                                      const typename Triangulation<deal_II_dimension,deal_II_space_dimension>::active_cell_iterator &,
+                                      const std::vector<bool>  &);
+
+
+                       \}
 
 #endif
 }
diff --git a/tests/grid/find_active_cell_around_point_02.cc b/tests/grid/find_active_cell_around_point_02.cc
new file mode 100644
index 0000000000..cd371b229a
--- /dev/null
+++ b/tests/grid/find_active_cell_around_point_02.cc
@@ -0,0 +1,83 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+/*
+ * Test cached version of find active cell around point
+ */
+#include "../tests.h"
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_tools_cache.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/fe_values.h>
+#include <deal.II/hp/mapping_collection.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/fe_field_function.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/numerics/vector_tools.h>
+
+
+template <int dim, int spacedim>
+void test(unsigned int n_ref, unsigned int n_points)
+{
+  deallog << "Testing " << dim << ", " << spacedim << std::endl;
+
+  Triangulation<dim,spacedim> tria;
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(n_ref);
+  DoFHandler<dim,spacedim> dof_handler(tria);
+
+  std::vector<Point<spacedim>> points;
+
+  deallog << "Points in study: " << n_points << std::endl;
+  for (size_t i=0; i<n_points; ++i)
+    {
+      //We need points in the square [0,1]x[0,1]: this is achieved normalizing with RAND_MAX
+      //RAND_MAX for the used algorithm is 2147483647
+      Point<spacedim> p;
+      for (unsigned int d=0; d<spacedim; ++d)
+        p[d] = double(Testing::rand())/RAND_MAX; //Normalizing the value
+      points.push_back(p);
+    }
+
+  auto &mapping = StaticMappingQ1<dim,spacedim>::mapping;
+
+  GridTools::Cache<dim,spacedim> cache(tria, mapping);
+
+  auto cell = tria.begin_active();
+  for (auto &p : points)
+    {
+      auto c_and_p = GridTools::find_active_cell_around_point(cache, p);
+      auto p2 = mapping.transform_unit_to_real_cell(c_and_p.first, c_and_p.second);
+      if (p2.distance(p)>1e-10)
+        deallog << "NOT OK!" << p << ", " << p2
+                << ", " << c_and_p.first << std::endl;
+      cell = c_and_p.first;
+    }
+  deallog << "OK" << std::endl;
+}
+
+int main ()
+{
+  initlog();
+
+  test<1,1> (4,10);
+  test<2,2> (3,20);
+  test<3,3> (2,30);
+}
diff --git a/tests/grid/find_active_cell_around_point_02.output b/tests/grid/find_active_cell_around_point_02.output
new file mode 100644
index 0000000000..02007871f7
--- /dev/null
+++ b/tests/grid/find_active_cell_around_point_02.output
@@ -0,0 +1,10 @@
+
+DEAL::Testing 1, 1
+DEAL::Points in study: 10
+DEAL::OK
+DEAL::Testing 2, 2
+DEAL::Points in study: 20
+DEAL::OK
+DEAL::Testing 3, 3
+DEAL::Points in study: 30
+DEAL::OK