From 15ef31fa7342cb4528a447925efaa01a643ab299 Mon Sep 17 00:00:00 2001
From: Bruno Blais <blais.bruno@gmail.com>
Date: Tue, 10 Mar 2020 17:25:47 -0400
Subject: [PATCH] Added two cached variables to the particle handler that
 stores the required information of the triangulation that is used to identify
 in which cells or subdomain the particle lies. This function is directly
 connected to the triangulation which enforces a reconstruction of the cached
 variable if the triangulation is changed in any way This change greatly
 boosts execution speed in the case where the number of cells and the number
 of particles is comparable.

---
 include/deal.II/particles/particle_handler.h | 40 ++++++++++++++++-
 source/particles/particle_handler.cc         | 45 ++++++++++++++------
 2 files changed, 72 insertions(+), 13 deletions(-)

diff --git a/include/deal.II/particles/particle_handler.h b/include/deal.II/particles/particle_handler.h
index 551f39133b..5fec613820 100644
--- a/include/deal.II/particles/particle_handler.h
+++ b/include/deal.II/particles/particle_handler.h
@@ -86,7 +86,10 @@ namespace Particles
     /**
      * Destructor.
      */
-    virtual ~ParticleHandler() override = default;
+    virtual ~ParticleHandler() override
+    {
+      connection.disconnect();
+    }
 
     /**
      * Initialize the particle handler. This function does not clear the
@@ -595,6 +598,16 @@ namespace Particles
     void
     serialize(Archive &ar, const unsigned int version);
 
+    /**
+     * Updates the structures that were generated to search in which cells the
+     * particles belong to. The call to this function updates the
+     * vertex_to_cells and the vertex_to_cell_centers data structure. This
+     * function is called every time the triangulation is altered.
+     *
+     */
+    void
+    update_triangulation_cache();
+
   private:
     /**
      * Address of the triangulation to work on.
@@ -696,6 +709,31 @@ namespace Particles
      */
     unsigned int handle;
 
+    /**
+     * A connection to the corresponding any_changes signal of the Triangulation
+     * which is attached to the particle_handler
+     */
+    boost::signals2::connection connection;
+
+    /** This variable stores a set from vertices to adjacent cells. It is used
+     * to locate the cells in which the particle reside. This structured is kept
+     * in memory to prevent its recalculation every time we
+     * sort_into_subdomain_and_cells(). This structure is updated everytime
+     * we call update_triangulation_cache.
+     */
+    std::vector<
+      std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
+      vertex_to_cells;
+
+
+    /** This variable stores a vector of vectors to cell centers. It is used
+     * to locate the cells in which the particle reside. This structured is kept
+     * in memory to prevent its recalculation every time we
+     * sort_into_subdomain_and_cells(). This structure is updated everytime
+     * we call update_triangulation_cache.
+     */
+    std::vector<std::vector<Tensor<1, spacedim>>> vertex_to_cell_centers;
+
 #ifdef DEAL_II_WITH_MPI
     /**
      * Transfer particles that have crossed subdomain boundaries to other
diff --git a/source/particles/particle_handler.cc b/source/particles/particle_handler.cc
index c0ed663d28..90b8729b50 100644
--- a/source/particles/particle_handler.cc
+++ b/source/particles/particle_handler.cc
@@ -119,7 +119,15 @@ namespace Particles
     , store_callback()
     , load_callback()
     , handle(numbers::invalid_unsigned_int)
-  {}
+  {
+    // Update the triangulation cache to ensure that you can sort the particles
+    // correctly
+    update_triangulation_cache();
+
+    connection = triangulation.signals.any_change.connect(
+      std::bind(&ParticleHandler<dim, spacedim>::update_triangulation_cache,
+                std::ref(*this)));
+  }
 
 
 
@@ -135,6 +143,14 @@ namespace Particles
 
     // Create the memory pool that will store all particle properties
     property_pool = std_cxx14::make_unique<PropertyPool>(n_properties);
+
+    // Update the triangulation cache to ensure that you can sort the particles
+    // correctly
+    update_triangulation_cache();
+
+    connection = triangulation->signals.any_change.connect(
+      std::bind(&ParticleHandler<dim, spacedim>::update_triangulation_cache,
+                std::ref(*this)));
   }
 
 
@@ -910,17 +926,7 @@ namespace Particles
         static_cast<vector_size>(particles_out_of_cell.size() * 0.25));
 
     {
-      // Create a map from vertices to adjacent cells
-      const std::vector<
-        std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
-        vertex_to_cells(GridTools::vertex_to_cell_map(*triangulation));
-
-      // Create a corresponding map of vectors from vertex to cell center
-      const std::vector<std::vector<Tensor<1, spacedim>>>
-        vertex_to_cell_centers(
-          GridTools::vertex_to_cell_centers_directions(*triangulation,
-                                                       vertex_to_cells));
-
+      // update_triangulation_cache();
       std::vector<unsigned int> neighbor_permutation;
 
       // Find the cells that the particles moved to.
@@ -1745,6 +1751,21 @@ namespace Particles
           break;
       }
   }
+
+  template <int dim, int spacedim>
+  void
+  ParticleHandler<dim, spacedim>::update_triangulation_cache()
+  {
+    // Create a map from vertices to adjacent cells
+    vertex_to_cells = GridTools::vertex_to_cell_map(*triangulation);
+
+    // Create a corresponding map of vectors from vertex to cell center
+    vertex_to_cell_centers =
+      GridTools::vertex_to_cell_centers_directions(*triangulation,
+                                                   vertex_to_cells);
+  }
+
+
 } // namespace Particles
 
 #include "particle_handler.inst"
-- 
2.39.5