]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Added new function compute point locations to GridTools 5314/head
authorGiovanni Alzetta <giovannialzetta@hotmail.it>
Tue, 24 Oct 2017 18:14:57 +0000 (18:14 +0000)
committerGiovanni Alzetta <giovannialzetta@hotmail.it>
Fri, 27 Oct 2017 10:33:12 +0000 (10:33 +0000)
doc/news/changes/minor/20171024GiovanniAlzetta [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/compute_point_locations_01.cc [new file with mode: 0644]
tests/grid/compute_point_locations_01.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20171024GiovanniAlzetta b/doc/news/changes/minor/20171024GiovanniAlzetta
new file mode 100644 (file)
index 0000000..6718a08
--- /dev/null
@@ -0,0 +1,3 @@
+Added: Function GridTools::compute_point_locations and a test for it. The function uses the recently added Cache to improve performance.
+<br>
+(Giovanni Alzetta, 2017/10/24)
index 8b27c9f9b4ffad59c24910e76e8d617fcc43d5ad..ca8f434d9d545d0c05d13bc589ff7e0ae1dc1c49 100644 (file)
@@ -605,6 +605,41 @@ namespace GridTools
    */
   /*@{*/
 
+  /**
+   * Given a @p cache and a list of @p points create the quadrature rules. This
+   * function returns a tuple containing the following elements:
+   * - The first ( get<0>), call it cells, is a vector of a vector cells of the all cells
+   *  that contain at least one of the @p points
+   * - The second a vector qpoints of vector of points, containing in qpoints[i]
+   *  the reference positions of all points that fall within the cell cells[i]
+   * - The third a vector indices of vector of integers, containing the mapping between
+   *  local numbering in qpoints, and global index in points
+   *
+   * If points[a] and points[b] are the only two points that fall in cells[c], then
+   * qpoints[c][0] and qpoints[c][1] are the reference positions of points[a] and points[b]
+   * in cells[c], and indices[c][0] = a, indices[c][1] = b. The function
+   * Mapping::tansform_unit_to_real(qpoints[c][0]) returns points[a].
+   *
+   * The algorithm assumes it's easier to look for a point in the cell that was used previously.
+   * For this reason random points are, computationally speaking, the worst case scenario while
+   * points grouped by the cell to which they belong are the best case.
+   * Pre-sorting points, trying to minimize distances, can make the algorithm much faster.
+   *
+   * Notice: given the center of a cell, points at distance greater then
+   * @p distance_factor * cell.diameter() are considered outside the cell. In some cases, such
+   * as curved meshes, this leads to cell's repetitions in the output tuple.
+   * To avoid this either enlarge @p distance_factor (depending on the regularity of the
+   * triangulation) and/or merge the repetitions contained in the output.
+   */
+  template <int dim, int spacedim>
+  std::tuple<
+  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator >,
+      std::vector< std::vector< Point<dim> > >,
+      std::vector< std::vector<unsigned int> > >
+      compute_point_locations(const Cache<dim,spacedim>                                         &cache,
+                              const std::vector<Point<spacedim> >                               &points,
+                              const double                                                      &distance_factor=0.5);
+
   /**
    * Return a map of index:Point<spacedim>, containing the used vertices of the
    * given `container`. The key of the returned map is the global index in the
index 196abb73298fea5e103529c437f985b86d34ac79..6dd0b32123b5ccfb197e41b9b70ea8af0f4a693e 100644 (file)
@@ -4915,6 +4915,177 @@ next_cell:
 
 
 
+  template <int dim, int spacedim>
+  std::tuple<
+  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator >,
+      std::vector< std::vector< Point<dim> > >,
+      std::vector< std::vector<unsigned int> > >
+      compute_point_locations(const Cache<dim,spacedim>                &cache,
+                              const std::vector<Point<spacedim> >      &points,
+                              const double                             &distance_factor)
+  {
+    // How many points are here?
+    const unsigned int np = points.size();
+
+    std::tuple<
+    std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator >,
+        std::vector< std::vector< Point<dim> > >,
+        std::vector< std::vector<unsigned int> > >
+        cell_qpoint_map;
+
+    // Now the easy case.
+    if (np==0) return cell_qpoint_map;
+
+    // If distance_factor is too small we might avoid points which are inside the cell
+    Assert(distance_factor >= 0.5,
+           ExcMessage("distance_factor value must be >= 0.5"));
+
+    // Keep track of the points we
+    // found
+    std::vector<bool> point_flags(np, false);
+
+    // Classify the first point:
+    // find active cell returns a pair (cell,transformed point)
+    const std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
+          Point<dim> >
+          my_pair  = GridTools::find_active_cell_around_point
+                     (cache, points[0]);
+
+    // Add data about the first point
+    std::get<0>(cell_qpoint_map).push_back(my_pair.first);
+    std::get<1>(cell_qpoint_map).emplace_back(1, my_pair.second);
+    std::get<2>(cell_qpoint_map).emplace_back(1, 0);
+
+    // Case only one point is now done
+    if (np==1)
+      return cell_qpoint_map;
+
+    // Compute the information about the current cell
+    Point<spacedim> cell_center = my_pair.first->center();
+    double cell_diameter = my_pair.first->diameter()*
+                           (distance_factor + std::numeric_limits<double>::epsilon() );
+
+    // Set this to true until all
+    // points have been classified
+    bool left_over = true;
+
+    // Flag to signal that a point is not in the current cell
+    bool flag_outside = false;
+
+    // This is the first index of a non processed point
+    unsigned int first_outside = 1;
+
+    // And this is the index of the current cell
+    unsigned int c = 0;
+
+    while (left_over == true)
+      {
+        // Assume this is the last one
+        left_over = false;
+        Assert(first_outside < np,
+               ExcIndexRange(first_outside, 0, np));
+
+        // If we found one in this cell, keep looking in the same cell
+        for (unsigned int p=first_outside; p<np; ++p)
+          if (point_flags[p] == false)
+            {
+              // Assume the point is in the cell
+              flag_outside = false;
+
+              // If the point is too far from the cell center it can't be inside it
+              if ( cell_center.distance(points[p]) > cell_diameter )
+                flag_outside = true;
+              else
+                {
+                  //The point is close to the cell: try transforming it
+                  try
+                    {
+                      Point<dim> qp =cache.get_mapping().transform_real_to_unit_cell
+                                     (std::get<0>(cell_qpoint_map)[c], points[p]);
+                      if (GeometryInfo<dim>::is_inside_unit_cell(qp))
+                        {
+                          point_flags[p] = true;
+                          std::get<1>(cell_qpoint_map)[c].push_back(qp);
+                          std::get<2>(cell_qpoint_map)[c].push_back(p);
+                        }
+                      else
+                        // The point is outside the cell
+                        flag_outside = true;
+                    }
+                  catch (const typename Mapping<dim>::ExcTransformationFailed &)
+                    {
+                      // transformation failed: assume the point is outside
+                      flag_outside = true;
+                    }
+                }
+
+              if (flag_outside)
+                {
+                  // Set things up for next round
+                  if (left_over == false)
+                    first_outside = p;
+                  left_over = true;
+                }
+            }
+        // If we got here and there is no left over, we are
+        // done. Else we need to find the next cell
+        if (left_over == true)
+          {
+            const std::pair<
+            typename Triangulation<dim, spacedim>::active_cell_iterator, Point<dim> >
+            my_pair
+              = GridTools::find_active_cell_around_point (cache, points[first_outside]);
+
+            std::get<0>(cell_qpoint_map).push_back(my_pair.first);
+            std::get<1>(cell_qpoint_map).emplace_back(1, my_pair.second);
+            std::get<2>(cell_qpoint_map).emplace_back(1, first_outside);
+            // Check if we can exit the loop now because only the last point was missing
+            if (first_outside == np-1)
+              left_over = false;
+
+            // Update the information about the current cell
+            cell_center = my_pair.first->center();
+            cell_diameter = my_pair.first->diameter()*
+                            (distance_factor + std::numeric_limits<double>::epsilon() );
+            c++;
+            point_flags[first_outside] = true;
+          }
+      }
+
+    // Augment of one the number of cells
+    ++c;
+    // Debug Checking
+    Assert(c == std::get<0>(cell_qpoint_map).size(), ExcInternalError());
+
+    Assert(c == std::get<2>(cell_qpoint_map).size(),
+           ExcDimensionMismatch(c, std::get<2>(cell_qpoint_map).size()));
+
+    Assert(c == std::get<1>(cell_qpoint_map).size(),
+           ExcDimensionMismatch(c, std::get<1>(cell_qpoint_map).size()));
+
+#ifdef DEBUG
+    unsigned int qps = 0;
+    // The number of points in all
+    // the cells must be the same as
+    // the number of points we
+    // started off from.
+    for (unsigned int n=0; n<c; ++n)
+      {
+        Assert(std::get<1>(cell_qpoint_map)[n].size() ==
+               std::get<2>(cell_qpoint_map)[n].size(),
+               ExcDimensionMismatch(std::get<1>(cell_qpoint_map)[n].size(),
+                                    std::get<2>(cell_qpoint_map)[n].size()));
+        qps += std::get<1>(cell_qpoint_map)[n].size();
+      }
+    Assert(qps == np,
+           ExcDimensionMismatch(qps, np));
+#endif
+
+    return cell_qpoint_map;
+  }
+
+
+
   template<int dim, int spacedim>
   std::map<unsigned int, Point<spacedim> >
   extract_used_vertices(const Triangulation<dim, spacedim> &container,
index 27020ca365057f8041ad93e93bd149f13e394253..c4de24d47308639897b204c5733ffc1b0d7329f5 100644 (file)
@@ -139,8 +139,13 @@ for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS
                                       const typename Triangulation<deal_II_dimension,deal_II_space_dimension>::active_cell_iterator &,
                                       const std::vector<bool>  &);
 
-
-                       \}
+        template
+        std::tuple< std::vector< typename Triangulation< deal_II_dimension, deal_II_space_dimension>::active_cell_iterator >,
+                    std::vector< std::vector< Point< deal_II_dimension > > >, std::vector< std::vector< unsigned int > > >
+        compute_point_locations(const Cache< deal_II_dimension, deal_II_space_dimension > &,
+                                const std::vector< Point< deal_II_space_dimension > > &,
+                                const double &);
+                                                                                                      \}
 
 #endif
 }
diff --git a/tests/grid/compute_point_locations_01.cc b/tests/grid/compute_point_locations_01.cc
new file mode 100644 (file)
index 0000000..334a19a
--- /dev/null
@@ -0,0 +1,112 @@
+// ---------------------------------------------------------------------
+//
+// 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 Functions::FEFieldFunction::compute_point_locations
+
+#include "../tests.h"
+#include <iostream>
+
+#include <deal.II/grid/tria.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/grid_tools_cache.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/numerics/fe_field_function.h>
+
+#include <deal.II/fe/mapping_q.h>
+
+using namespace dealii;
+
+template <int dim>
+void test_compute_pt_loc(unsigned int n_points)
+{
+  deallog << "Testing for dim = " << dim << std::endl;
+  deallog << "Testing on: " << n_points << " points." << std::endl;
+
+  // Creating a grid in the square [0,1]x[0,1]
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(std::max(6-dim,2));
+
+  //Creating the finite elements needed:
+  FE_Q<dim> fe(1);
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+
+  // Creating the random points
+  std::vector<Point<dim>> points;
+
+  for (size_t i=0; i<n_points; ++i)
+    {
+      Point<dim> p;
+      for (unsigned int d=0; d<dim; ++d)
+        p[d] = double(Testing::rand())/RAND_MAX; //Normalizing the value
+      points.push_back(p);
+    }
+
+  // Initializing the cache
+  GridTools::Cache<dim,dim> cache(tria);
+
+  std::tuple< std::vector<typename Triangulation<dim, dim>::active_cell_iterator >, std::vector< std::vector< Point<dim> > >, std::vector<std::vector<unsigned int> > >
+  cell_qpoint_map = GridTools::compute_point_locations(cache,points);
+  size_t n_cells = std::get<0>(cell_qpoint_map).size();
+
+  deallog << "Points found in " << n_cells << " cells" << std::endl;
+
+  // testing if the transformation is correct:
+  // For each cell check if the quadrature points in the i-th FE
+  // are the same as std::get<2>(cell_qpoint_map)[i]
+  for (unsigned int i=0; i<std::get<0>(cell_qpoint_map).size(); ++i)
+    {
+      auto &cell = std::get<0>(cell_qpoint_map)[i];
+      auto &quad = std::get<1>(cell_qpoint_map)[i];
+      auto &local_map = std::get<2>(cell_qpoint_map)[i];
+
+      // Given the std::get<1>(cell_qpoint_map) of the current cell, compute the real points
+      FEValues<dim> fev(fe, quad, update_quadrature_points);
+      fev.reinit(cell);
+      const auto &real_quad = fev.get_quadrature_points();
+
+      for (unsigned int q=0; q<real_quad.size(); ++q)
+        {
+          // Check if points are the same as real points
+          if (real_quad[q].distance(points[local_map[q]]) > 1e-10)
+            deallog << "Error on cell : " << cell
+                    << " at local point " << i
+                    << ", corresponding to real point "
+                    << points[local_map[q]]
+                    << ", that got transformed to "
+                    << real_quad[q]
+                    << " instead." << std::endl;
+        }
+    }
+  deallog << "Test finished" << std::endl;
+}
+
+int main()
+{
+  initlog();
+
+  deallog << "Deal.II compute_pt_loc:" << std::endl;
+  test_compute_pt_loc<2>(100);
+  test_compute_pt_loc<3>(200);
+}
diff --git a/tests/grid/compute_point_locations_01.output b/tests/grid/compute_point_locations_01.output
new file mode 100644 (file)
index 0000000..0a6e0e1
--- /dev/null
@@ -0,0 +1,10 @@
+
+DEAL::Deal.II compute_pt_loc:
+DEAL::Testing for dim = 2
+DEAL::Testing on: 100 points.
+DEAL::Points found in 81 cells
+DEAL::Test finished
+DEAL::Testing for dim = 3
+DEAL::Testing on: 200 points.
+DEAL::Points found in 170 cells
+DEAL::Test finished

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.