From 33ebdd8233df7909b17e9a4509528fe13ed733c9 Mon Sep 17 00:00:00 2001 From: bangerth Date: Mon, 14 Apr 2014 18:51:59 +0000 Subject: [PATCH] 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 --- deal.II/doc/news/changes.h | 7 +++++ deal.II/source/grid/grid_tools.cc | 43 +++++++++++++++++-------------- 2 files changed, 31 insertions(+), 19 deletions(-) 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.

Specific improvements

    +
  1. 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. +
    + (Giorgos Kourakos, Timo Heister, Wolfgang Bangerth, 2014/04/14) +
    2. +
  2. Changed: TableBase now uses AlignedVector for storing data instead of std::vector, which allows its use for VectorizedArray 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 &container, const Point &p) { - typedef typename Container::active_cell_iterator cell_iterator; + typedef typename Container::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 > best_cell; + std::pair > best_cell; // Find closest vertex and determine // all adjacent cells - unsigned int vertex = find_closest_vertex(container, p); - - std::vector adjacent_cells_tmp = - find_cells_adjacent_to_vertex(container, vertex); + std::vector 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 adjacent_cells(adjacent_cells_tmp.begin(), adjacent_cells_tmp.end()); - std::set searched_cells; + std::set adjacent_cells (adjacent_cells_tmp.begin(), + adjacent_cells_tmp.end()); + std::set 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::const_iterator + typename std::set::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 (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 - void fix_up_faces (const typename dealii::Triangulation::cell_iterator &cell, + void fix_up_faces (const typename dealii::Triangulation::active_cell_iterator &cell, dealii::internal::int2type, dealii::internal::int2type) { -- 2.39.5