]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Make GridTools::delete_duplicate_vertices() faster. 9297/head
authorDavid Wells <drwells@email.unc.edu>
Sun, 12 Jan 2020 04:35:08 +0000 (23:35 -0500)
committerDavid Wells <drwells@email.unc.edu>
Tue, 14 Jan 2020 03:02:48 +0000 (22:02 -0500)
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 [new file with mode: 0644]
include/deal.II/grid/grid_tools.h
source/grid/grid_tools.cc

diff --git a/doc/news/changes/minor/20200112DavidWells b/doc/news/changes/minor/20200112DavidWells
new file mode 100644 (file)
index 0000000..f3025ee
--- /dev/null
@@ -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.
+<br>
+(David Wells, 2020/01/12)
index 495d281137fd9422ba1cc928c63ac21cc7bfed27..7c225bdbbb148795b9940ad24905553f25cc0e05 100644 (file)
@@ -374,10 +374,11 @@ namespace GridTools
    *
    * This function is called by some <tt>GridIn::read_*</tt> 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.
index 4b48ff14134495c4102f91533ca93beae97d081d..1378bcaea0fd1fe9ed39d54e4d8782b28dc34873 100644 (file)
@@ -848,53 +848,111 @@ namespace GridTools
                              std::vector<CellData<dim>> &  cells,
                              SubCellData &                 subcelldata,
                              std::vector<unsigned int> &   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<unsigned int> 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<spacedim> 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<std::pair<unsigned int, Point<spacedim>>> 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<unsigned int, Point<spacedim>> &a,
+                  const std::pair<unsigned int, Point<spacedim>> &b) {
+                return a.second[longest_coordinate_direction] <
+                       b.second[longest_coordinate_direction];
+              });
+
+    auto within_tolerance = [=](const Point<spacedim> &a,
+                                const Point<spacedim> &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<unsigned int, Point<spacedim>> &a,
+                     const std::pair<unsigned int, Point<spacedim>> &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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.