From b1d9bc8e69302918aea77b72c03b521847c9609a Mon Sep 17 00:00:00 2001 From: David Wells Date: Sat, 11 Jan 2020 23:35:08 -0500 Subject: [PATCH] Make GridTools::delete_duplicate_vertices() faster. This change makes this function, for cube-like geometries, O(n^(3/2)) in 2D and O(n^(5/3)) in 3D rather than O(n^2) by presorting vertices by value in a single component and then only comparing vertices which are nearby in that component. --- doc/news/changes/minor/20200112DavidWells | 5 + include/deal.II/grid/grid_tools.h | 9 +- source/grid/grid_tools.cc | 130 ++++++++++++++++------ 3 files changed, 104 insertions(+), 40 deletions(-) create mode 100644 doc/news/changes/minor/20200112DavidWells diff --git a/doc/news/changes/minor/20200112DavidWells b/doc/news/changes/minor/20200112DavidWells new file mode 100644 index 0000000000..f3025ee255 --- /dev/null +++ b/doc/news/changes/minor/20200112DavidWells @@ -0,0 +1,5 @@ +Improved: GridTools::delete_duplicated_vertices() now runs, for cubelike +geometries, in $O(n^{3/2})$ time in 2D and $O(n^(5/3))$ time in 3D instead +of $O(n^2)$ time. +
+(David Wells, 2020/01/12) diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index 495d281137..7c225bdbbb 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -374,10 +374,11 @@ namespace GridTools * * This function is called by some GridIn::read_* functions. Only * the vertices with indices in @p considered_vertices are tested for - * equality. This speeds up the algorithm, which is quadratic and thus quite - * slow to begin with. However, if you wish to consider all vertices, simply - * pass an empty vector. In that case, the function fills - * @p considered_vertices with all vertices. + * equality. This speeds up the algorithm, which is, for worst-case hyper + * cube geometries $O(N^{3/2})$ in 2D and $O(N^{5/3})$ in 3D: quite slow. + * However, if you wish to consider all vertices, simply pass an empty + * vector. In that case, the function fills @p considered_vertices with all + * vertices. * * Two vertices are considered equal if their difference in each coordinate * direction is less than @p tol. diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 4b48ff1413..1378bcaea0 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -848,53 +848,111 @@ namespace GridTools std::vector> & cells, SubCellData & subcelldata, std::vector & considered_vertices, - double tol) + const double tol) { - // create a vector of vertex - // indices. initialize it to the identity, - // later on change that if necessary. + AssertIndexRange(2, vertices.size()); + // create a vector of vertex indices. initialize it to the identity, later + // on change that if necessary. std::vector new_vertex_numbers(vertices.size()); - for (unsigned int i = 0; i < vertices.size(); ++i) - new_vertex_numbers[i] = i; + std::iota(new_vertex_numbers.begin(), new_vertex_numbers.end(), 0); - // if the considered_vertices vector is - // empty, consider all vertices + // if the considered_vertices vector is empty, consider all vertices if (considered_vertices.size() == 0) considered_vertices = new_vertex_numbers; - Assert(considered_vertices.size() <= vertices.size(), ExcInternalError()); + // The algorithm below improves upon the naive O(n^2) algorithm by first + // sorting vertices by their value in one component and then only + // comparing vertices for equality which are nearly equal in that + // component. For example, if @p vertices form a cube, then we will only + // compare points that have the same x coordinate when we try to find + // duplicated vertices. + + // Start by finding the longest coordinate direction. This minimizes the + // number of points that need to be compared against each-other in a + // single set for typical geometries. + const BoundingBox bbox(vertices); + const auto & min = bbox.get_boundary_points().first; + const auto & max = bbox.get_boundary_points().second; + + unsigned int longest_coordinate_direction = 0; + double longest_coordinate_length = max[0] - min[0]; + for (unsigned int d = 1; d < spacedim; ++d) + { + const double coordinate_length = max[d] - min[d]; + if (longest_coordinate_length < coordinate_length) + { + longest_coordinate_length = coordinate_length; + longest_coordinate_direction = d; + } + } - // now loop over all vertices to be - // considered and try to find an identical - // one - for (unsigned int i = 0; i < considered_vertices.size(); ++i) + // Sort vertices (while preserving their vertex numbers) along that + // coordinate direction: + std::vector>> sorted_vertices; + sorted_vertices.reserve(vertices.size()); + for (const unsigned int vertex_n : considered_vertices) { - Assert(considered_vertices[i] < vertices.size(), ExcInternalError()); - if (new_vertex_numbers[considered_vertices[i]] != - considered_vertices[i]) - // this vertex has been identified with - // another one already, skip it in the - // test - continue; - // this vertex is not identified with - // another one so far. search in the list - // of remaining vertices. if a duplicate - // vertex is found, set the new vertex - // index for that vertex to this vertex' - // index. - for (unsigned int j = i + 1; j < considered_vertices.size(); ++j) + AssertIndexRange(vertex_n, vertices.size()); + sorted_vertices.emplace_back(vertex_n, vertices[vertex_n]); + } + std::sort(sorted_vertices.begin(), + sorted_vertices.end(), + [&](const std::pair> &a, + const std::pair> &b) { + return a.second[longest_coordinate_direction] < + b.second[longest_coordinate_direction]; + }); + + auto within_tolerance = [=](const Point &a, + const Point &b) { + for (unsigned int d = 0; d < spacedim; ++d) + if (std::abs(a[d] - b[d]) > tol) + return false; + return true; + }; + + // Find a range of numbers that have the same component in the longest + // coordinate direction: + auto range_start = sorted_vertices.begin(); + while (range_start != sorted_vertices.end()) + { + auto range_end = range_start + 1; + while (range_end != sorted_vertices.end() && + std::abs(range_end->second[longest_coordinate_direction] - + range_start->second[longest_coordinate_direction]) < + tol) + ++range_end; + + // preserve behavior with older versions of this function by replacing + // higher vertex numbers by lower vertex numbers + std::sort(range_start, + range_end, + [](const std::pair> &a, + const std::pair> &b) { + return a.first < b.first; + }); + + // Now de-duplicate [range_start, range_end) + // + // We have identified all points that are within a strip of width 'tol' + // in one coordinate direction. Now we need to figure out which of these + // are also close in other coordinate directions. If two are close, we + // can mark the second one for deletion. + for (auto reference = range_start; reference != range_end; ++reference) { - bool equal = true; - for (unsigned int d = 0; d < spacedim; ++d) - equal &= (std::abs(vertices[considered_vertices[j]](d) - - vertices[considered_vertices[i]](d)) < tol); - if (equal) - { - new_vertex_numbers[considered_vertices[j]] = - considered_vertices[i]; - } + if (reference->first != numbers::invalid_unsigned_int) + for (auto it = reference + 1; it != range_end; ++it) + { + if (within_tolerance(reference->second, it->second)) + { + new_vertex_numbers[it->first] = reference->first; + // skip the replaced vertex in the future + it->first = numbers::invalid_unsigned_int; + } + } } + range_start = range_end; } // now we got a renumbering list. simply renumber all vertices -- 2.39.5