From eb0dd9604bcc9f74833b8485aa11dccbf235d4be Mon Sep 17 00:00:00 2001
From: Luca Heltai <luca.heltai@sissa.it>
Date: Sat, 12 Aug 2017 15:52:49 -0600
Subject: [PATCH] Get closest vertex to cell.

---
 include/deal.II/grid/grid_tools.h | 14 ++++++++++++++
 source/grid/grid_tools.cc         | 22 ++++++++++++++++++++++
 2 files changed, 36 insertions(+)

diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h
index 1a2532d306..32d3bc5245 100644
--- a/include/deal.II/grid/grid_tools.h
+++ b/include/deal.II/grid/grid_tools.h
@@ -1023,6 +1023,20 @@ namespace GridTools
   std::vector<std::set<typename Triangulation<dim,spacedim>::active_cell_iterator> >
   vertex_to_cell_map(const Triangulation<dim,spacedim> &triangulation);
 
+  /**
+
+
+  /**
+   * Returns the local vertex index of cell @p cell that is closest to
+   * the given location @p position.
+   *
+   * @author Rene Gassmoeller, Luca Heltai, 2017.
+   */
+  template <int dim, int spacedim>
+  unsigned int
+  get_closest_vertex_of_cell(const typename Triangulation<dim,spacedim>::active_cell_iterator &cell,
+                             const Point<spacedim> &position);
+
   /**
    * Compute a globally unique index for each vertex and hanging node
    * associated with a locally owned active cell. The vertices of a ghost cell
diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc
index 85a9194788..6c9e9d5a38 100644
--- a/source/grid/grid_tools.cc
+++ b/source/grid/grid_tools.cc
@@ -1460,6 +1460,28 @@ next_cell:
 
 
 
+  template <int dim, int spacedim>
+  unsigned int
+  get_closest_vertex_of_cell(const typename Triangulation<dim,spacedim>::active_cell_iterator &cell,
+                             const Point<spacedim> &position)
+  {
+    double minimum_distance = std::numeric_limits<double>::max();
+    unsigned int closest_vertex = numbers::invalid_unsigned_int;
+
+    for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+      {
+        const double vertex_distance = position.distance(cell->vertex(v));
+        if (vertex_distance < minimum_distance)
+          {
+            closest_vertex = v;
+            minimum_distance = vertex_distance;
+          }
+      }
+    return closest_vertex;
+  }
+
+
+
   template <int dim, int spacedim>
   std::pair<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator, Point<dim> >
   find_active_cell_around_point (const hp::MappingCollection<dim,spacedim>   &mapping,
-- 
2.39.5