From: bangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Mon, 14 Apr 2014 18:51:59 +0000 (+0000)
Subject: Fix a bug in GridTools::find_active_cell_around_point() where we accidentally used... 
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=33ebdd8233df7909b17e9a4509528fe13ed733c9;p=dealii-svn.git

Fix a bug in GridTools::find_active_cell_around_point() where we accidentally used tria.n_cells() instead of tria.n_active_cells().

git-svn-id: https://svn.dealii.org/trunk@32777 0785d39b-7218-0410-832d-ea1e28bc413d
---

diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index 5c2f730fe7..899dd3ce74 100644
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -148,6 +148,13 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li> Fixed: GridTools::find_active_cell_around_point() could get into an infinite
+  loop if the point we are looking for is in fact not within the domain. This is now
+  fixed.
+  <br>
+  (Giorgos Kourakos, Timo Heister, Wolfgang Bangerth, 2014/04/14)
+  </li>
+
   <li> Changed: TableBase<N,T> now uses AlignedVector for storing data
   instead of std::vector, which allows its use for VectorizedArray<Number>
   data fields which require more alignment.
diff --git a/deal.II/source/grid/grid_tools.cc b/deal.II/source/grid/grid_tools.cc
index 34888a24bb..d8ecd6269d 100644
--- a/deal.II/source/grid/grid_tools.cc
+++ b/deal.II/source/grid/grid_tools.cc
@@ -890,7 +890,7 @@ next_cell:
                                  const Container<dim,spacedim> &container,
                                  const Point<spacedim>     &p)
   {
-    typedef typename Container<dim,spacedim>::active_cell_iterator cell_iterator;
+    typedef typename Container<dim,spacedim>::active_cell_iterator active_cell_iterator;
 
     // The best distance is set to the
     // maximum allowable distance from
@@ -898,22 +898,22 @@ next_cell:
     // max. deviation of 1e-10
     double best_distance = 1e-10;
     int    best_level = -1;
-    std::pair<cell_iterator, Point<dim> > best_cell;
+    std::pair<active_cell_iterator, Point<dim> > best_cell;
 
     // Find closest vertex and determine
     // all adjacent cells
-    unsigned int vertex = find_closest_vertex(container, p);
-
-    std::vector<cell_iterator> adjacent_cells_tmp =
-      find_cells_adjacent_to_vertex(container, vertex);
+    std::vector<active_cell_iterator> adjacent_cells_tmp
+    = find_cells_adjacent_to_vertex(container,
+                                    find_closest_vertex(container, p));
 
     // Make sure that we have found
     // at least one cell adjacent to vertex.
     Assert(adjacent_cells_tmp.size()>0, ExcInternalError());
 
     // Copy all the cells into a std::set
-    std::set<cell_iterator> adjacent_cells(adjacent_cells_tmp.begin(), adjacent_cells_tmp.end());
-    std::set<cell_iterator> searched_cells;
+    std::set<active_cell_iterator> adjacent_cells (adjacent_cells_tmp.begin(),
+                                                   adjacent_cells_tmp.end());
+    std::set<active_cell_iterator> searched_cells;
 
     // Determine the maximal number of cells
     // in the grid.
@@ -921,12 +921,12 @@ next_cell:
     // the cell and have not searched
     // every cell in the triangulation,
     // we keep on looking.
-    const unsigned int n_cells =get_tria(container).n_cells();
+    const unsigned int n_active_cells = get_tria(container).n_active_cells();
     bool found = false;
     unsigned int cells_searched = 0;
-    while (!found && cells_searched < n_cells)
+    while (!found && cells_searched < n_active_cells)
       {
-        typename std::set<cell_iterator>::const_iterator
+        typename std::set<active_cell_iterator>::const_iterator
         cell = adjacent_cells.begin(),
         endc = adjacent_cells.end();
         for (; cell != endc; ++cell)
@@ -943,10 +943,13 @@ next_cell:
                 // unit cell (or at least not too far
                 // outside). If it is, it is also checked
                 // that the cell has a more refined state
-                if (dist < best_distance ||
-                    (dist == best_distance && (*cell)->level() > best_level))
+                if ((dist < best_distance)
+                    ||
+                    ((dist == best_distance)
+                     &&
+                     ((*cell)->level() > best_level)))
                   {
-                    found       = true;
+                    found         = true;
                     best_distance = dist;
                     best_level    = (*cell)->level();
                     best_cell     = std::make_pair(*cell, p_cell);
@@ -967,16 +970,18 @@ next_cell:
                 // to the next
               }
           }
-        //udpate the number of cells searched
+
+        // update the number of cells searched
         cells_searched += adjacent_cells.size();
+
         // if we have not found the cell in
         // question and have not yet searched every
         // cell, we expand our search to
         // all the not already searched neighbors of
         // the cells in adjacent_cells. This is
-        // what find_active_cell_around_poin_internal
+        // what find_active_cell_around_point_internal
         // is for.
-        if (!found && cells_searched < n_cells)
+        if (!found && cells_searched < n_active_cells)
           {
             find_active_cell_around_point_internal<dim,Container,spacedim>
             (container, searched_cells, adjacent_cells);
@@ -1974,7 +1979,7 @@ next_cell:
 
 
 
-      void fix_up_faces (const dealii::Triangulation<1,1>::cell_iterator &,
+      void fix_up_faces (const dealii::Triangulation<1,1>::active_cell_iterator &,
                          dealii::internal::int2type<1>,
                          dealii::internal::int2type<1>)
       {
@@ -1988,7 +1993,7 @@ next_cell:
       // a cell by moving around its
       // mid-points
       template <int structdim, int spacedim>
-      void fix_up_faces (const typename dealii::Triangulation<structdim,spacedim>::cell_iterator &cell,
+      void fix_up_faces (const typename dealii::Triangulation<structdim,spacedim>::active_cell_iterator &cell,
                          dealii::internal::int2type<structdim>,
                          dealii::internal::int2type<spacedim>)
       {