]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Improving and fixing GridTools::compute_point_locations 5468/head
authorGiovanni Alzetta <giovannialzetta@hotmail.it>
Wed, 15 Nov 2017 11:05:59 +0000 (11:05 +0000)
committerGiovanni Alzetta <giovannialzetta@hotmail.it>
Wed, 15 Nov 2017 11:51:24 +0000 (11:51 +0000)
doc/news/changes/minor/20171115GiovanniAlzetta [new file with mode: 0644]
include/deal.II/grid/grid_tools.h
source/grid/grid_tools.cc
source/grid/grid_tools.inst.in

diff --git a/doc/news/changes/minor/20171115GiovanniAlzetta b/doc/news/changes/minor/20171115GiovanniAlzetta
new file mode 100644 (file)
index 0000000..187dd26
--- /dev/null
@@ -0,0 +1,3 @@
+Fixed and improved: The GridTools::compute_point_locations() now always return the correct number of cells. The new algorithm is also significantly faster.
+<br>
+(Giovanni Alzetta, 2017/11/15)
index 98f90d21dc831fa8a820fa69074f1ad13d210751..cf84dd095e488aac66dd137ad7f3b42b4790fa0b 100644 (file)
@@ -622,30 +622,32 @@ 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].
+   * Given a Triangulation's @p cache and a list of @p points create the quadrature rules.
+   *
+   * @param[in] cache The triangulation's GridTools::Cache .
+   * @param[in] points The point's vector.
+   *
+   * @param[out] Tuple containing the following information:
+   *  - Cells, is a vector of a vector cells of the all cells
+   *   containing at least one of the @p points .
+   *  - A vector qpoints of vector of points, containing in @p qpoints[i]
+   *   the reference positions of all points that fall within the cell @P cells[i] .
+   *  - A vector indices of vector of integers, containing the mapping between
+   *   local numbering in qpoints, and global index in points
+   *
+   * If @p points[a] and @p points[b] are the only two points that fall in @p cells[c],
+   * then @p qpoints[c][0] and @p qpoints[c][1] are the reference positions of
+   * @p points[a] and @p points[b] in @p cells[c], and @p indices[c][0] = a,
+   * @p indices[c][1] = b. The function Mapping::tansform_unit_to_real(qpoints[c][0])
+   * returns @p 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.
+   * Pre-sorting points, trying to minimize distances between them, might make the function
+   * extremely 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.
+   * @author Giovanni Alzetta, 2017
    */
   template <int dim, int spacedim>
   std::tuple<
@@ -653,8 +655,7 @@ namespace GridTools
       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);
+                              const std::vector<Point<spacedim> >                               &points);
 
   /**
    * Return a map of index:Point<spacedim>, containing the used vertices of the
index 45317013248a63b9a965391c4890a9203c848dd2..c278b805222dd09f22802c9e26a4d157d34a40d1 100644 (file)
@@ -4921,8 +4921,7 @@ next_cell:
       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)
+                              const std::vector<Point<spacedim> >      &points)
   {
     // How many points are here?
     const unsigned int np = points.size();
@@ -4936,134 +4935,77 @@ next_cell:
     // 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"));
+    // We begin by finding the cell/transform of the first point
+    std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator, Point<dim> >
+    my_pair  = GridTools::find_active_cell_around_point
+               (cache, points[0]);
 
-    // 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<0>(cell_qpoint_map).emplace_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;
+    // Now the second easy case.
+    if (np==1) return cell_qpoint_map;
+    // Computing the cell center and diameter
+    Point<spacedim> cell_center = std::get<0>(cell_qpoint_map)[0]->center();
+    double cell_diameter = std::get<0>(cell_qpoint_map)[0]->diameter()*
+                           (0.5 + std::numeric_limits<double>::epsilon() );
 
-    // 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)
+    // Cycle over all points left
+    for (unsigned int p=1; p< np; ++p)
       {
-        // 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;
-                    }
-                }
+        // Checking if the point is close to the cell center, in which
+        // case calling find active cell with a cell hint
+        if ( cell_center.distance(points[p]) < cell_diameter )
+          my_pair  = GridTools::find_active_cell_around_point
+                     (cache, points[p],std::get<0>(cell_qpoint_map).back());
+        else
+          my_pair  = GridTools::find_active_cell_around_point
+                     (cache, points[p]);
 
-              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)
+        // Assuming the cell is probably the last cell added
+        if ( my_pair.first == std::get<0>(cell_qpoint_map).back() )
+          {
+            // Found in the last cell: adding the data
+            std::get<1>(cell_qpoint_map).back().emplace_back(my_pair.second);
+            std::get<2>(cell_qpoint_map).back().emplace_back(p);
+          }
+        else
           {
-            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;
+            // Check if it is in another cell already found
+            typename std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>::iterator
+            cells_it = std::find(std::get<0>(cell_qpoint_map).begin(),std::get<0>(cell_qpoint_map).end()-1,my_pair.first);
+
+            if ( cells_it == std::get<0>(cell_qpoint_map).end()-1 )
+              {
+                // Cell not found: adding a new cell
+                std::get<0>(cell_qpoint_map).emplace_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, p);
+                // Updating center and radius of the cell
+                cell_center = std::get<0>(cell_qpoint_map).back()->center();
+                cell_diameter = std::get<0>(cell_qpoint_map).back()->diameter()*
+                                (0.5 + std::numeric_limits<double>::epsilon() );
+              }
+            else
+              {
+                unsigned int current_cell = cells_it - std::get<0>(cell_qpoint_map).begin();
+                // Cell found: just adding the point index and qpoint to the list
+                std::get<1>(cell_qpoint_map)[current_cell].emplace_back(my_pair.second);
+                std::get<2>(cell_qpoint_map)[current_cell].emplace_back(p);
+              }
           }
       }
 
-    // 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(std::get<0>(cell_qpoint_map).size() == std::get<2>(cell_qpoint_map).size(),
+           ExcDimensionMismatch(std::get<0>(cell_qpoint_map).size(), 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()));
+    Assert(std::get<0>(cell_qpoint_map).size() == std::get<1>(cell_qpoint_map).size(),
+           ExcDimensionMismatch(std::get<0>(cell_qpoint_map).size(), std::get<1>(cell_qpoint_map).size()));
 
 #ifdef DEBUG
+    unsigned int c = std::get<0>(cell_qpoint_map).size();
     unsigned int qps = 0;
     // The number of points in all
     // the cells must be the same as
index ee1cceff892ca44141965651a410601e7289226c..a695d87a3b4296ee5e37898708609f2d08c323e8 100644 (file)
@@ -143,8 +143,7 @@ for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS
         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 &);
+                                const std::vector< Point< deal_II_space_dimension > > &);
                                                                                                       \}
 
 #endif

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.