From 855892954fc7add876d05989ebb0c13c1ed77ef4 Mon Sep 17 00:00:00 2001
From: Magdalena Schreter <schreter.magdalena@gmail.com>
Date: Thu, 29 Sep 2022 21:20:59 +0200
Subject: [PATCH] update

---
 include/deal.II/grid/grid_tools.h | 18 ++++++-
 source/grid/grid_tools.cc         | 82 ++++++++++++-------------------
 2 files changed, 48 insertions(+), 52 deletions(-)

diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h
index 25aa2220ce..989afb5076 100644
--- a/include/deal.II/grid/grid_tools.h
+++ b/include/deal.II/grid/grid_tools.h
@@ -449,6 +449,18 @@ namespace GridTools
                              std::vector<unsigned int> &   considered_vertices,
                              const double                  tol = 1e-12);
 
+  /**
+   * Remove vertices that are duplicated.
+   *
+   * Two vertices are considered equal if their difference in each coordinate
+   * direction is less than @p tol. This implies that nothing happens if
+   * the tolerance is set to zero.
+   */
+  template <int dim>
+  void
+  delete_duplicated_vertices(std::vector<Point<dim>> &vertices,
+                             const double             tol = 1e-12);
+
   /**
    * Grids generated by grid generators may have an orientation of cells which
    * is the inverse of the orientation required by deal.II.
@@ -3441,9 +3453,11 @@ namespace GridTools
      * Process all locally-owned cells and fill @p vertices and @p cells for all
      * cells that are cut.
      *
-     * @note: This function is only implemented for dim>1. Use
+     * @note This function is only implemented for dim>1. Use
      * process(background_dof_handler, ls_vector, iso_level, vertices) for
      * dim==1.
+     *
+     * @note Duplicate vertices are not deleted.
      */
     void
     process(const DoFHandler<dim> &  background_dof_handler,
@@ -3468,7 +3482,7 @@ namespace GridTools
      *
      * @note The resulting vectors are empty if the cell is not cut.
      *
-     * @note: This function is only implemented for dim>1. Use
+     * @note This function is only implemented for dim>1. Use
      * process_cell(cell, ls_vector, iso_level, vertices) for dim==1.
      */
     void
diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc
index 6eb7200e6d..ea1cc381ca 100644
--- a/source/grid/grid_tools.cc
+++ b/source/grid/grid_tools.cc
@@ -868,6 +868,37 @@ namespace GridTools
 
 
 
+  template <int dim>
+  void
+  delete_duplicated_vertices(std::vector<Point<dim>> &vertices,
+                             const double             tol)
+  {
+    if (vertices.size() == 0)
+      return;
+
+    // 1) map point to local vertex index
+    std::map<Point<dim>, unsigned int, FloatingPointComparator<double>>
+      map_point_to_local_vertex_index{FloatingPointComparator<double>(tol)};
+
+    // 2) initialize map with existing points uniquely
+    for (unsigned int i = 0; i < vertices.size(); ++i)
+      map_point_to_local_vertex_index[vertices[i]] = i;
+
+    // no duplicate points are found
+    if (map_point_to_local_vertex_index.size() == vertices.size())
+      return;
+
+    // 3) remove duplicate entries from vertices
+    vertices.resize(map_point_to_local_vertex_index.size());
+    {
+      unsigned int j = 0;
+      for (const auto &p : map_point_to_local_vertex_index)
+        vertices[j++] = p.first;
+    }
+  }
+
+
+
   template <int dim, int spacedim>
   std::size_t
   invert_cells_with_negative_measure(
@@ -6744,51 +6775,6 @@ namespace GridTools
             }
         }
     }
-
-
-
-    template <int dim>
-    void
-    merge_duplicate_points(std::vector<Point<dim>> &vertices,
-                           std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells)
-    {
-      if (vertices.size() == 0)
-        return;
-
-      // 1) map point to local vertex index
-      std::map<Point<dim>, unsigned int, FloatingPointComparator<double>>
-        map_point_to_local_vertex_index(FloatingPointComparator<double>(1e-10));
-
-      // 2) initialize map with existing points uniquely
-      for (unsigned int i = 0; i < vertices.size(); ++i)
-        map_point_to_local_vertex_index[vertices[i]] = i;
-
-      // no duplicate points are found
-      if (map_point_to_local_vertex_index.size() == vertices.size())
-        return;
-
-      // 3) remove duplicate entries from vertices
-      vertices.resize(map_point_to_local_vertex_index.size());
-      {
-        unsigned int j = 0;
-        for (const auto &p : map_point_to_local_vertex_index)
-          vertices[j++] = p.first;
-      }
-
-      // 4) renumber vertices in CellData object
-      if (cells.size() > 0)
-        {
-          map_point_to_local_vertex_index.clear();
-          // initialize map with unique points
-          for (unsigned int i = 0; i < vertices.size(); ++i)
-            map_point_to_local_vertex_index[vertices[i]] = i;
-
-          // overwrite vertex indices in CellData object
-          for (auto &cell : cells)
-            for (auto &v : cell.vertices)
-              v = map_point_to_local_vertex_index[vertices[v]];
-        }
-    }
   } // namespace internal
 
 
@@ -6862,8 +6848,6 @@ namespace GridTools
     for (const auto &cell : background_dof_handler.active_cell_iterators() |
                               IteratorFilters::LocallyOwnedCell())
       process_cell(cell, ls_vector, iso_level, vertices, cells);
-
-    internal::merge_duplicate_points(vertices, cells);
   }
 
   template <int dim, typename VectorType>
@@ -6878,9 +6862,7 @@ namespace GridTools
                               IteratorFilters::LocallyOwnedCell())
       process_cell(cell, ls_vector, iso_level, vertices);
 
-    // This vector is just a placeholder to reuse the process_cell function.
-    std::vector<CellData<dim == 1 ? 1 : dim - 1>> dummy_cells;
-    internal::merge_duplicate_points(vertices, dummy_cells);
+    delete_duplicated_vertices(vertices, 1e-10 /*tol*/);
   }
 
 
-- 
2.39.5