]> https://gitweb.dealii.org/ - dealii.git/commitdiff
New find active cell around point. 4794/head
authorLuca Heltai <luca.heltai@sissa.it>
Sat, 12 Aug 2017 21:53:32 +0000 (15:53 -0600)
committerLuca Heltai <luca.heltai@sissa.it>
Mon, 14 Aug 2017 19:09:12 +0000 (13:09 -0600)
cmake/config/template-arguments.in
doc/news/changes/minor/20170813LucaHeltai [new file with mode: 0644]
include/deal.II/grid/grid_tools.h
source/grid/grid_tools.cc
source/grid/grid_tools.inst.in
tests/grid/find_active_cell_around_point_01.cc [new file with mode: 0644]
tests/grid/find_active_cell_around_point_01.output [new file with mode: 0644]

index f6a6523a486e16da76855a8ba06aeeba5ec769ef..a07bc908c2099f61e07bd1c0cddc5a2d94629c11 100644 (file)
@@ -156,6 +156,12 @@ TRIANGULATION_AND_DOFHANDLERS := { Triangulation<deal_II_dimension, deal_II_spac
                                    DoFHandler<deal_II_dimension, deal_II_space_dimension>;
                                    hp::DoFHandler<deal_II_dimension, deal_II_space_dimension> }
 
+// concrete Triangulation types (all types you can iterate over cells from)
+// with hard-coded <deal_II_dimension, deal_II_space_dimension>
+TRIANGULATIONS := { Triangulation<deal_II_dimension, deal_II_space_dimension>;
+                    parallel::shared::Triangulation<deal_II_dimension, deal_II_space_dimension>;
+                    parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension>;}
+
 // serial and parallel sparsity pattern types
 SPARSITY_PATTERNS := { SparsityPattern;
                        DynamicSparsityPattern;
diff --git a/doc/news/changes/minor/20170813LucaHeltai b/doc/news/changes/minor/20170813LucaHeltai
new file mode 100644 (file)
index 0000000..4485d82
--- /dev/null
@@ -0,0 +1,5 @@
+New: A new version of GridTools::find_active_cell_around_point()
+has been added that exploits local maps constructed 
+using standard GridTools utilities.
+<br>
+(Rene Gassmoeller, Luca Heltai, 2017/08/13)
index 4e16025de79205bbe01a930c186267d1d53d1927..77cc9e5153125938416ea693a1ac54d50e9b9e5e 100644 (file)
@@ -747,6 +747,28 @@ namespace GridTools
                                  const Point<spacedim>        &p,
                                  const std::vector<bool>      &marked_vertices = std::vector<bool>());
 
+  /**
+   * A version of the previous function that exploits an already existing
+   * map between vertices and cells, constructed using the function
+   * GridTools::vertex_to_cell_map, a map of vertex_to_cell_centers, obtained
+   * through GridTools::vertex_to_cells_center_directions, and a guess `cell_hint`.
+   *
+   * @author Luca Heltai, Rene Gassmoeller, 2017
+   */
+  template <int dim, template <int, int> class MeshType, int spacedim>
+#ifndef _MSC_VER
+  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim> >
+#else
+  std::pair<typename dealii::internal::ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim> >::type, Point<dim> >
+#endif
+  find_active_cell_around_point (const Mapping<dim,spacedim>                                                          &mapping,
+                                 const MeshType<dim,spacedim>                                                         &mesh,
+                                 const Point<spacedim>                                                                &p,
+                                 const std::vector<std::set<typename MeshType<dim,spacedim>::active_cell_iterator > > &vertex_to_cell_map,
+                                 const std::vector<std::vector<Tensor<1,spacedim> > >                                 &vertex_to_cell_centers,
+                                 const typename MeshType<dim, spacedim>::active_cell_iterator                         &cell_hint=typename MeshType<dim, spacedim>::active_cell_iterator(),
+                                 const std::vector<bool>                                                              &marked_vertices = std::vector<bool>());
+
   /**
    * A version of the previous function where we use that mapping on a given
    * cell that corresponds to the active finite element index of that cell.
@@ -1024,10 +1046,16 @@ namespace GridTools
   vertex_to_cell_map(const Triangulation<dim,spacedim> &triangulation);
 
   /**
-   * Returns a vector that contains a tensor for every vertex-cell combination
-   * of the output of GridTools::vertex_to_cell_map() (which is expected as
-   * input parameter for this function). Each tensor represents a geometric
-   * vector from the vertex to the respective cell center.
+   * Returns a vector of normalized tensors for each vertex-cell combination of
+   * the output of GridTools::vertex_to_cell_map() (which is expected as input
+   * parameter for this function). Each tensor represents a geometric vector
+   * from the vertex to the respective cell center.
+   *
+   * An assertion will be thrown if the size of the input vector is not equal to
+   * the number of vertices of the triangulation.
+   *
+   * result[v][c] is a unit Tensor for vertex index v, indicating the direction of
+   * the center of the c-th cell with respect to the vertex v.
    *
    * @author Rene Gassmoeller, Luca Heltai, 2017.
    */
@@ -1045,8 +1073,8 @@ namespace GridTools
    */
   template <int dim, int spacedim>
   unsigned int
-  get_closest_vertex_of_cell(const typename Triangulation<dim,spacedim>::active_cell_iterator &cell,
-                             const Point<spacedim> &position);
+  find_closest_vertex_of_cell(const typename Triangulation<dim,spacedim>::active_cell_iterator &cell,
+                              const Point<spacedim> &position);
 
   /**
    * Compute a globally unique index for each vertex and hanging node
index 920fd0f78303d72bae77cc6fdf7430254ad56db0..91728bb39b3fbdd47c41633b108f02da601738da 100644 (file)
@@ -1466,6 +1466,9 @@ next_cell:
     const std::vector<Point<spacedim> > &vertices = mesh.get_vertices();
     const unsigned int n_vertices = vertex_to_cells.size();
 
+    AssertDimension(vertices.size(), n_vertices);
+
+
     std::vector<std::vector<Tensor<1,spacedim> > > vertex_to_cell_centers(n_vertices);
     for (unsigned int vertex=0; vertex<n_vertices; ++vertex)
       if (mesh.vertex_used(vertex))
@@ -1484,18 +1487,137 @@ next_cell:
   }
 
 
+  namespace
+  {
+    template <int spacedim>
+    bool
+    compare_point_association(const unsigned int a,
+                              const unsigned int b,
+                              const Tensor<1,spacedim> &point_direction,
+                              const std::vector<Tensor<1,spacedim> > &center_directions)
+    {
+      const double scalar_product_a = center_directions[a] * point_direction;
+      const double scalar_product_b = center_directions[b] * point_direction;
+
+      // The function is supposed to return if a is before b. We are looking
+      // for the alignment of point direction and center direction, therefore
+      // return if the scalar product of a is larger.
+      return (scalar_product_a > scalar_product_b);
+    }
+  }
+
+  template <int dim, template <int, int> class MeshType, int spacedim>
+#ifndef _MSC_VER
+  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim> >
+#else
+  std::pair<typename dealii::internal::ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim> >::type, Point<dim> >
+#endif
+  find_active_cell_around_point (const Mapping<dim,spacedim>                                                    &mapping,
+                                 const MeshType<dim,spacedim>                                                   &mesh,
+                                 const Point<spacedim>                                                          &p,
+                                 const std::vector<std::set<typename MeshType<dim,spacedim>::active_cell_iterator > > &vertex_to_cells,
+                                 const std::vector<std::vector<Tensor<1,spacedim> > >                            &vertex_to_cell_centers,
+                                 const typename MeshType<dim, spacedim>::active_cell_iterator                    &cell_hint ,
+                                 const std::vector<bool>                                                         &marked_vertices)
+  {
+    std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim> > cell_and_position;
+
+    bool found_cell = false;
+
+    unsigned int closest_vertex_index = 0;
+    Tensor<1,spacedim> vertex_to_point;
+    auto current_cell = cell_hint;
+
+    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)
+          {
+            const unsigned int closest_vertex = find_closest_vertex_of_cell<dim,spacedim>(current_cell , p);
+            vertex_to_point = p - current_cell ->vertex(closest_vertex);
+            closest_vertex_index = current_cell ->vertex_index(closest_vertex);
+          }
+        else
+          {
+            closest_vertex_index = GridTools::find_closest_vertex(mesh,p,marked_vertices);
+            vertex_to_point = p - mesh.get_vertices()[closest_vertex_index];
+          }
+
+        vertex_to_point /= vertex_to_point.norm();
+        const unsigned int n_neighbor_cells = vertex_to_cells[closest_vertex_index].size();
+
+        // Create a corresponding map of vectors from vertex to cell center
+        std::vector<unsigned int> neighbor_permutation(n_neighbor_cells);
+
+        for (unsigned int i=0; i<n_neighbor_cells; ++i)
+          neighbor_permutation[i] = i;
+
+        auto comp = [&](const unsigned int a, const unsigned int b) -> bool
+        {
+          return compare_point_association<spacedim>(a,b,vertex_to_point,vertex_to_cell_centers[closest_vertex_index]);
+        };
+
+        std::sort(neighbor_permutation.begin(),
+                  neighbor_permutation.end(),
+                  comp);
+
+        // Search all of the cells adjacent to the closest vertex of the cell hint
+        // Most likely we will find the point in them.
+        for (unsigned int i=0; i<n_neighbor_cells; ++i)
+          {
+            try
+              {
+                auto cell = vertex_to_cells[closest_vertex_index].begin();
+                std::advance(cell,neighbor_permutation[i]);
+                const Point<dim> p_unit = mapping.transform_real_to_unit_cell(*cell, p);
+                if (GeometryInfo<dim>::is_inside_unit_cell(p_unit))
+                  {
+                    cell_and_position.first = *cell;
+                    cell_and_position.second = p_unit;
+                    found_cell = true;
+                    break;
+                  }
+              }
+            catch (typename Mapping<dim>::ExcTransformationFailed &)
+              {}
+          }
+
+        if (found_cell == true)
+          return cell_and_position;
+
+        // 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
+        // for a global vertex close to the point. If that does not work, we are in
+        // trouble, and just throw an exception.
+        //
+        // If we got here, then we did not find the point. If the
+        // current_cell.state() here is not IteratorState::valid, it means that
+        // the user did not provide a hint_cell, and at the beginning of the
+        // while loop we performed an actual global search on the mesh
+        // vertices. Not finding the point then means the point is outside the
+        // domain.
+        AssertThrow(current_cell.state() == IteratorState::valid,
+                    ExcPointNotFound<spacedim>(p));
+
+        current_cell = typename MeshType<dim,spacedim>::active_cell_iterator();
+      }
+    return cell_and_position;
+  }
+
+
 
   template <int dim, int spacedim>
   unsigned int
-  get_closest_vertex_of_cell(const typename Triangulation<dim,spacedim>::active_cell_iterator &cell,
-                             const Point<spacedim> &position)
+  find_closest_vertex_of_cell(const typename Triangulation<dim,spacedim>::active_cell_iterator &cell,
+                              const Point<spacedim> &position)
   {
-    double minimum_distance = std::numeric_limits<double>::max();
-    unsigned int closest_vertex = numbers::invalid_unsigned_int;
+    double minimum_distance = position.distance_square(cell->vertex(0));
+    unsigned int closest_vertex = 0;
 
-    for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+    for (unsigned int v=1; v<GeometryInfo<dim>::vertices_per_cell; ++v)
       {
-        const double vertex_distance = position.distance(cell->vertex(v));
+        const double vertex_distance = position.distance_square(cell->vertex(v));
         if (vertex_distance < minimum_distance)
           {
             closest_vertex = v;
index 5f97b7adb165e4f86270c94f356740df4dd62a1f..2b3f29b8e1d8ecdea841ae7bd871531a30e6ab59 100644 (file)
 //
 // ---------------------------------------------------------------------
 
+for (X : TRIANGULATIONS; deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS)
+{
+
+#if deal_II_dimension <= deal_II_space_dimension
+    namespace GridTools \{
+    template
+    std::pair<dealii::internal::ActiveCellIterator<deal_II_dimension, deal_II_space_dimension, X>::type, Point<deal_II_dimension> >
+    find_active_cell_around_point (const Mapping<deal_II_dimension, deal_II_space_dimension> &,
+                                   const X &,
+                                   const Point<deal_II_space_dimension> &,
+                                   const std::vector<std::set<typename dealii::internal::ActiveCellIterator<deal_II_dimension, deal_II_space_dimension, X>::type> > &,
+                                   const std::vector<std::vector<Tensor<1,deal_II_space_dimension> > > &,
+                                   const dealii::internal::ActiveCellIterator<deal_II_dimension, deal_II_space_dimension, X>::type &,
+                                   const std::vector<bool> &);
+    \}
+
+#endif
+}
 
 
 for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS)
@@ -101,6 +119,8 @@ for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS
         template
         std::map<unsigned int,types::global_vertex_index>
         compute_local_to_global_vertex_index_map(const parallel::distributed::Triangulation<deal_II_dimension,deal_II_space_dimension> &triangulation);
+
+
                 \}
 
 #endif
@@ -241,6 +261,12 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
     std::vector<std::set<Triangulation<deal_II_dimension,deal_II_space_dimension>::active_cell_iterator> >
     vertex_to_cell_map(const Triangulation<deal_II_dimension,deal_II_space_dimension> &triangulation);
 
+    template
+    std::vector<std::vector<Tensor<1,deal_II_space_dimension> > >
+    vertex_to_cell_centers_directions(const Triangulation<deal_II_dimension,deal_II_space_dimension> &mesh,
+                                      const std::vector<std::set<typename Triangulation<deal_II_dimension,
+                                      deal_II_space_dimension>::active_cell_iterator> > &vertex_to_cells);
+
 #if deal_II_dimension == deal_II_space_dimension
 #  if deal_II_dimension > 1
     template
diff --git a/tests/grid/find_active_cell_around_point_01.cc b/tests/grid/find_active_cell_around_point_01.cc
new file mode 100644 (file)
index 0000000..e8cec95
--- /dev/null
@@ -0,0 +1,85 @@
+// ---------------------------------------------------------------------
+//
+// 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/base/utilities.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/grid/grid_tools.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>
+double 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 v_to_c = GridTools::vertex_to_cell_map (tria);
+  auto v_to_c_d = GridTools::vertex_to_cell_centers_directions (tria, v_to_c);
+
+  auto &mapping = StaticMappingQ1<dim,spacedim>::mapping;
+  auto cell = tria.begin_active();
+  for (auto &p : points)
+    {
+      auto c_and_p = GridTools::find_active_cell_around_point(mapping, tria, p,
+                                                              v_to_c, v_to_c_d);
+      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_01.output b/tests/grid/find_active_cell_around_point_01.output
new file mode 100644 (file)
index 0000000..0200787
--- /dev/null
@@ -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

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.