From 18d23d0484d5b6141d0174c1f9462d5581b0e88a Mon Sep 17 00:00:00 2001 From: Giovanni Alzetta Date: Wed, 15 Nov 2017 11:05:59 +0000 Subject: [PATCH] Improving and fixing GridTools::compute_point_locations --- .../changes/minor/20171115GiovanniAlzetta | 3 + include/deal.II/grid/grid_tools.h | 43 ++--- source/grid/grid_tools.cc | 172 ++++++------------ source/grid/grid_tools.inst.in | 3 +- 4 files changed, 83 insertions(+), 138 deletions(-) create mode 100644 doc/news/changes/minor/20171115GiovanniAlzetta diff --git a/doc/news/changes/minor/20171115GiovanniAlzetta b/doc/news/changes/minor/20171115GiovanniAlzetta new file mode 100644 index 0000000000..187dd26c54 --- /dev/null +++ b/doc/news/changes/minor/20171115GiovanniAlzetta @@ -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. +
+(Giovanni Alzetta, 2017/11/15) diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index 98f90d21dc..cf84dd095e 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -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 std::tuple< @@ -653,8 +655,7 @@ namespace GridTools std::vector< std::vector< Point > >, std::vector< std::vector > > compute_point_locations(const Cache &cache, - const std::vector > &points, - const double &distance_factor=0.5); + const std::vector > &points); /** * Return a map of index:Point, containing the used vertices of the diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 4531701324..c278b80522 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -4921,8 +4921,7 @@ next_cell: std::vector< std::vector< Point > >, std::vector< std::vector > > compute_point_locations(const Cache &cache, - const std::vector > &points, - const double &distance_factor) + const std::vector > &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::active_cell_iterator, Point > + my_pair = GridTools::find_active_cell_around_point + (cache, points[0]); - // Keep track of the points we - // found - std::vector point_flags(np, false); - - // Classify the first point: - // find active cell returns a pair (cell,transformed point) - const std::pair::active_cell_iterator, - Point > - 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 cell_center = my_pair.first->center(); - double cell_diameter = my_pair.first->diameter()* - (distance_factor + std::numeric_limits::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 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::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 cell_diameter ) - flag_outside = true; - else - { - //The point is close to the cell: try transforming it - try - { - Point qp =cache.get_mapping().transform_real_to_unit_cell - (std::get<0>(cell_qpoint_map)[c], points[p]); - if (GeometryInfo::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::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::active_cell_iterator, Point > - 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::epsilon() ); - c++; - point_flags[first_outside] = true; + // Check if it is in another cell already found + typename std::vector::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::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 diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in index ee1cceff89..a695d87a3b 100644 --- a/source/grid/grid_tools.inst.in +++ b/source/grid/grid_tools.inst.in @@ -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 -- 2.39.5