From 36f9318e08ae72903d5b7b7e315b0229ee793af0 Mon Sep 17 00:00:00 2001
From: Bruno <blais.bruno@gmail.com>
Date: Sun, 25 Oct 2020 14:39:22 -0400
Subject: [PATCH] Added vertex_to_neighbhor_subdomain to cache and use it
 therein the particle_handler

---
 include/deal.II/grid/grid_tools_cache.h       | 16 +++++++++++++
 .../grid/grid_tools_cache_update_flags.h      |  9 +++++--
 source/grid/grid_tools_cache.cc               | 19 +++++++++++++++
 source/particles/particle_handler.cc          | 24 +++++++------------
 4 files changed, 50 insertions(+), 18 deletions(-)

diff --git a/include/deal.II/grid/grid_tools_cache.h b/include/deal.II/grid/grid_tools_cache.h
index 26140d8097..f61b1f3403 100644
--- a/include/deal.II/grid/grid_tools_cache.h
+++ b/include/deal.II/grid/grid_tools_cache.h
@@ -152,6 +152,15 @@ namespace GridTools
                 typename Triangulation<dim, spacedim>::active_cell_iterator>> &
     get_locally_owned_cell_bounding_boxes_rtree() const;
 
+
+    /** Returns the vector of set of integer containing the subdomain id
+     * to which each vertex is connected to. This feature is used extensively
+     * in the particle_handler to detect on which processors ghost particles
+     * must be built
+     */
+    const std::vector<std::set<unsigned int>> &
+    get_vertex_to_neighbor_subdomain() const;
+
     /**
      * Return a reference to the stored triangulation.
      */
@@ -269,6 +278,13 @@ namespace GridTools
                 typename Triangulation<dim, spacedim>::active_cell_iterator>>
       locally_owned_cell_bounding_boxes_rtree;
 
+
+    /**
+     * Store an std::vector of std::set of integer containing the id of all
+     * subdomain to which a vertex is connected to
+     */
+    mutable std::vector<std::set<unsigned int>> vertex_to_neighbor_subdomain;
+
     /**
      * Storage for the status of the triangulation signal.
      */
diff --git a/include/deal.II/grid/grid_tools_cache_update_flags.h b/include/deal.II/grid/grid_tools_cache_update_flags.h
index 2a89e3236c..c8d3e4c05b 100644
--- a/include/deal.II/grid/grid_tools_cache_update_flags.h
+++ b/include/deal.II/grid/grid_tools_cache_update_flags.h
@@ -79,6 +79,11 @@ namespace GridTools
      */
     update_locally_owned_cell_bounding_boxes_rtree = 0x080,
 
+    /**
+     * Update vertex to neighbhor subdomain
+     */
+    update_vertex_to_neighbor_subdomain = 0x100,
+
     /**
      * Update all objects.
      */
@@ -163,8 +168,8 @@ namespace GridTools
    *
    * @ref CacheUpdateFlags
    */
-  inline CacheUpdateFlags operator&(const CacheUpdateFlags f1,
-                                    const CacheUpdateFlags f2)
+  inline CacheUpdateFlags
+  operator&(const CacheUpdateFlags f1, const CacheUpdateFlags f2)
   {
     return static_cast<CacheUpdateFlags>(static_cast<unsigned int>(f1) &
                                          static_cast<unsigned int>(f2));
diff --git a/source/grid/grid_tools_cache.cc b/source/grid/grid_tools_cache.cc
index 9cd4a8704b..46986e3ca4 100644
--- a/source/grid/grid_tools_cache.cc
+++ b/source/grid/grid_tools_cache.cc
@@ -199,6 +199,25 @@ namespace GridTools
     return covering_rtree[level];
   }
 
+  template <int dim, int spacedim>
+  const std::vector<std::set<unsigned int>> &
+  Cache<dim, spacedim>::get_vertex_to_neighbor_subdomain() const
+  {
+    if (update_flags & update_vertex_to_neighbor_subdomain)
+      {
+        vertex_to_neighbor_subdomain.clear();
+        vertex_to_neighbor_subdomain.resize(tria->n_vertices());
+        for (const auto &cell : tria->active_cell_iterators())
+          {
+            if (cell->is_ghost())
+              for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
+                vertex_to_neighbor_subdomain[cell->vertex_index(v)].insert(
+                  cell->subdomain_id());
+          }
+      }
+    return vertex_to_neighbor_subdomain;
+  }
+
 #include "grid_tools_cache.inst"
 
 } // namespace GridTools
diff --git a/source/particles/particle_handler.cc b/source/particles/particle_handler.cc
index 83cfac4a13..3f3503aead 100644
--- a/source/particles/particle_handler.cc
+++ b/source/particles/particle_handler.cc
@@ -1186,16 +1186,8 @@ namespace Particles
         static_cast<typename std::vector<particle_iterator>::size_type>(
           particles.size() * 0.25));
 
-    std::vector<std::set<unsigned int>> vertex_to_neighbor_subdomain(
-      triangulation->n_vertices());
-
-    for (const auto &cell : triangulation->active_cell_iterators())
-      {
-        if (cell->is_ghost())
-          for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
-            vertex_to_neighbor_subdomain[cell->vertex_index(v)].insert(
-              cell->subdomain_id());
-      }
+    const std::vector<std::set<unsigned int>> vertex_to_neighbor_subdomain =
+      triangulation_cache->get_vertex_to_neighbor_subdomain();
 
     for (const auto &cell : triangulation->active_cell_iterators())
       {
@@ -1687,8 +1679,8 @@ namespace Particles
 
     switch (status)
       {
-        case parallel::distributed::Triangulation<dim, spacedim>::CELL_PERSIST:
-          {
+          case parallel::distributed::Triangulation<dim,
+                                                    spacedim>::CELL_PERSIST: {
             auto position_hint = particles.end();
             for (const auto &particle : loaded_particles_on_cell)
               {
@@ -1707,8 +1699,8 @@ namespace Particles
           }
           break;
 
-        case parallel::distributed::Triangulation<dim, spacedim>::CELL_COARSEN:
-          {
+          case parallel::distributed::Triangulation<dim,
+                                                    spacedim>::CELL_COARSEN: {
             typename std::multimap<internal::LevelInd,
                                    Particle<dim, spacedim>>::iterator
               position_hint = particles.end();
@@ -1733,8 +1725,8 @@ namespace Particles
           }
           break;
 
-        case parallel::distributed::Triangulation<dim, spacedim>::CELL_REFINE:
-          {
+          case parallel::distributed::Triangulation<dim,
+                                                    spacedim>::CELL_REFINE: {
             std::vector<
               typename std::multimap<internal::LevelInd,
                                      Particle<dim, spacedim>>::iterator>
-- 
2.39.5