Added two cached variables to the particle handler that stores the required informati...

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.

diff --git a/include/deal.II/particles/particle_handler.h b/include/deal.II/particles/particle_handler.h
index 551f39133b2c15ac4f8335e0fa0457e05e45099c..5fec6138200054830d6c37ba02a8c020bc68c0a3 100644 (file)
--- 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 c0ed663d28416d5ab885704acc77422c8bc12678..90b8729b504c01754789ec21297c6ebfbae82a27 100644 (file)
--- 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"