From: goll Date: Sun, 22 Jul 2012 13:21:09 +0000 (+0000) Subject: Fixed bug in find_active_cell_around_point. See the testcases bits/find_active_cell_8... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=92fff7e79d05497c51d23f6e5c68a8de48fa8cae;p=dealii-svn.git Fixed bug in find_active_cell_around_point. See the testcases bits/find_active_cell_8 and _9 git-svn-id: https://svn.dealii.org/trunk@25715 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/source/grid/grid_tools.cc b/deal.II/source/grid/grid_tools.cc index e81c6eb272..eedaf0a5c6 100644 --- a/deal.II/source/grid/grid_tools.cc +++ b/deal.II/source/grid/grid_tools.cc @@ -781,6 +781,54 @@ namespace GridTools + namespace + { + template class Container, int spacedim> + void find_active_cell_around_point_internal(const Container& container, + std::set::active_cell_iterator>& searched_cells, + std::set::active_cell_iterator>& adjacent_cells) + { + typedef typename Container::active_cell_iterator cell_iterator; + + // update the searched cells + searched_cells.insert(adjacent_cells.begin(), adjacent_cells.end()); + // now we to collect all neighbors + // of the cells in adjacent_cells we + // have not yet searched. + std::set adjacent_cells_new; + + typename std::set::const_iterator + cell = adjacent_cells.begin(), + endc = adjacent_cells.end(); + for(; cell != endc; ++cell) + { + std::vector active_neighbors; + get_active_neighbors >(*cell, active_neighbors); + for (unsigned int i=0; i class Container, int spacedim> typename Container::active_cell_iterator find_active_cell_around_point (const Container &container, @@ -811,51 +859,83 @@ namespace GridTools // all adjacent cells unsigned int vertex = find_closest_vertex(container, p); - std::vector adjacent_cells = + std::vector adjacent_cells_tmp = find_cells_adjacent_to_vertex(container, vertex); - - typename std::vector::const_iterator - cell = adjacent_cells.begin(), - endc = adjacent_cells.end(); - - for(; cell != endc; ++cell) - { - try - { - const Point p_cell = mapping.transform_real_to_unit_cell(*cell, p); - - // calculate the infinity norm of - // the distance vector to the unit cell. - const double dist = GeometryInfo::distance_to_unit_cell(p_cell); - - // We compare if the point is inside the - // 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)) - { - best_distance = dist; - best_level = (*cell)->level(); - best_cell = std::make_pair(*cell, p_cell); - } - } - catch (typename MappingQ1::ExcTransformationFailed &) - { - // ok, the transformation - // failed presumably - // because the point we - // are looking for lies - // outside the current - // cell. this means that - // the current cell can't - // be the cell around the - // point, so just ignore - // this cell and move on - // to the next - } - - } + + // 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; + + // Determine the maximal number of cells + // in the grid. + // As long as we have not found + // 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(); + bool found = false; + unsigned int cells_searched = 0; + while (!found && cells_searched < n_cells) + { + typename std::set::const_iterator + cell = adjacent_cells.begin(), + endc = adjacent_cells.end(); + for(; cell != endc; ++cell) + { + try + { + const Point p_cell = mapping.transform_real_to_unit_cell(*cell, p); + + // calculate the infinity norm of + // the distance vector to the unit cell. + const double dist = GeometryInfo::distance_to_unit_cell(p_cell); + + // We compare if the point is inside the + // 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)) + { + found = true; + best_distance = dist; + best_level = (*cell)->level(); + best_cell = std::make_pair(*cell, p_cell); + } + } + catch (typename MappingQ1::ExcTransformationFailed &) + { + // ok, the transformation + // failed presumably + // because the point we + // are looking for lies + // outside the current + // cell. this means that + // the current cell can't + // be the cell around the + // point, so just ignore + // this cell and move on + // to the next + } + } + //udpate 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 + // is for. + if(!found && cells_searched < n_cells) + { + find_active_cell_around_point_internal(container, searched_cells, adjacent_cells); + } + } Assert (best_cell.first.state() == IteratorState::valid, ExcPointNotFound(p)); @@ -901,55 +981,88 @@ namespace GridTools // all adjacent cells unsigned int vertex = find_closest_vertex(container, p); - std::vector adjacent_cells = - find_cells_adjacent_to_vertex(container, vertex); - - typename std::vector::const_iterator - cell = adjacent_cells.begin(), - endc = adjacent_cells.end(); - - for(; cell != endc; ++cell) - { - try - { - const Point p_cell - = mapping[(*cell)->active_fe_index()].transform_real_to_unit_cell(*cell, p); - - // calculate the infinity norm of - // the distance vector to the unit cell. - const double dist = GeometryInfo::distance_to_unit_cell(p_cell); - - // We compare if the point is inside the - // 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)) + std::vector adjacent_cells_tmp = + find_cells_adjacent_to_vertex(container, vertex); + + // 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; + + // Determine the maximal number of cells + // in the grid. + // As long as we have not found + // 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(); + bool found = false; + unsigned int cells_searched = 0; + while (!found && cells_searched < n_cells) + { + typename std::set::const_iterator + cell = adjacent_cells.begin(), + endc = adjacent_cells.end(); + for(; cell != endc; ++cell) { - best_distance = dist; - best_level = (*cell)->level(); - best_cell = std::make_pair(*cell, p_cell); + try + { + const Point p_cell = mapping[(*cell)->active_fe_index()].transform_real_to_unit_cell(*cell, p); + + + // calculate the infinity norm of + // the distance vector to the unit cell. + const double dist = GeometryInfo::distance_to_unit_cell(p_cell); + + // We compare if the point is inside the + // 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)) + { + found = true; + best_distance = dist; + best_level = (*cell)->level(); + best_cell = std::make_pair(*cell, p_cell); + } + } + catch (typename MappingQ1::ExcTransformationFailed &) + { + // ok, the transformation + // failed presumably + // because the point we + // are looking for lies + // outside the current + // cell. this means that + // the current cell can't + // be the cell around the + // point, so just ignore + // this cell and move on + // to the next + } + } + //udpate 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. + if(!found && cells_searched < n_cells) + { + find_active_cell_around_point_internal(container, searched_cells, adjacent_cells); } - } - catch (typename MappingQ1::ExcTransformationFailed &) - { - // ok, the transformation - // failed presumably - // because the point we - // are looking for lies - // outside the current - // cell. this means that - // the current cell can't - // be the cell around the - // point, so just ignore - // this cell and move on - // to the next - } - } - Assert (best_cell.first.state() == IteratorState::valid, - ExcPointNotFound(p)); - } + } + } + + Assert (best_cell.first.state() == IteratorState::valid, + ExcPointNotFound(p)); + return best_cell; }