First WIP version. Issue with exchange ghosts

diff --git a/include/deal.II/particles/particle_accessor.h b/include/deal.II/particles/particle_accessor.h
index d2c7455e1f3ce9d499b66549d92be11f2b6afe14..41c9e125ee5672e3d1ac0ea0f2bbe6b940112faa 100644 (file)
--- a/include/deal.II/particles/particle_accessor.h
+++ b/include/deal.II/particles/particle_accessor.h
@@ -713,8 +713,11 @@ namespace Particles
 
     ++particle_index_within_cell;
 
+
     if (particle_index_within_cell > (*particles)[active_cell_index].size() - 1)
       {
+        const bool initial_cell_is_owned = cell->is_locally_owned();
+
         particle_index_within_cell = 0;
 
         do
@@ -723,7 +726,8 @@ namespace Particles
             ++active_cell_index;
           }
         while (cell.state() == IteratorState::valid &&
-               (*particles)[active_cell_index].size() == 0);
+               ((*particles)[active_cell_index].size() == 0 ||
+                cell->is_locally_owned() != initial_cell_is_owned));
       }
   }
 
@@ -739,6 +743,8 @@ namespace Particles
       --particle_index_within_cell;
     else
       {
+        const bool initial_cell_is_owned = cell->is_locally_owned();
+
         do
           {
             --cell;
@@ -750,14 +756,16 @@ namespace Particles
                 break;
               }
 
-            if ((*particles)[active_cell_index].size() > 0)
+            if ((*particles)[active_cell_index].size() > 0 &&
+                cell->is_locally_owned() == initial_cell_is_owned)
               {
                 particle_index_within_cell =
                   (*particles)[active_cell_index].size() - 1;
                 break;
               }
           }
-        while ((*particles)[active_cell_index].size() == 0);
+        while ((*particles)[active_cell_index].size() == 0 &&
+               cell->is_locally_owned() != initial_cell_is_owned);
       }
   }
 

diff --git a/include/deal.II/particles/particle_handler.h b/include/deal.II/particles/particle_handler.h
index 4cc8088af27b59d130ca73fbf71a648e13849a24..b86a57e9dcfcda29ff1aec0a1a29795fefbe5a99 100644 (file)
--- a/include/deal.II/particles/particle_handler.h
+++ b/include/deal.II/particles/particle_handler.h
@@ -208,7 +208,9 @@ namespace Particles
     /**
      * Return the number of particles that live on the given cell.
      *
-     * @note While this function is used in step-19, it is not an efficient
+     * @note
+     *  //TODO this not true now no?
+     *  While this function is used in step-19, it is not an efficient
      *   function to use if the number of particles is large. That is because
      *   to find the particles that are located in one cell costs
      *   ${\cal O}(\log N)$ where $N$ is the number of overall particles. Since
@@ -847,17 +849,11 @@ namespace Particles
     std::unique_ptr<PropertyPool<dim, spacedim>> property_pool;
 
     /**
-     * Set of particles currently living in the local domain, organized by
-     * the level/index of the cell they are in.
+     * Set of particles currently living in the local domain including ghost
+     * cells , organized by the level/index of the cell they are in.
      */
     particle_container particles;
 
-    /**
-     * Set of particles currently living in the local domain, organized by
-     * the level/index of the cell they are in.
-     */
-    particle_container ghost_particles;
-
     /**
      * This variable stores how many particles are stored globally. It is
      * calculated by update_cached_numbers().
@@ -1113,8 +1109,8 @@ namespace Particles
 
     for (const auto &cell : triangulation->active_cell_iterators())
       if (cell->is_locally_owned() == false &&
-          ghost_particles[cell->active_cell_index()].size() != 0)
-        return particle_iterator(ghost_particles, *property_pool, cell, 0);
+          particles[cell->active_cell_index()].size() != 0)
+        return particle_iterator(particles, *property_pool, cell, 0);
 
     return end_ghost();
   }
@@ -1134,7 +1130,7 @@ namespace Particles
   inline typename ParticleHandler<dim, spacedim>::particle_iterator
   ParticleHandler<dim, spacedim>::end_ghost()
   {
-    return particle_iterator(ghost_particles,
+    return particle_iterator(particles,
                              *property_pool,
                              triangulation->end(),
                              0);

diff --git a/source/particles/particle_handler.cc b/source/particles/particle_handler.cc
index ebab1fa0809fff8efa703d74de04892531786ba9..c820e18ef866f2db7d9f7299ae433c4e0700cb4f 100644 (file)
--- a/source/particles/particle_handler.cc
+++ b/source/particles/particle_handler.cc
@@ -94,7 +94,6 @@ namespace Particles
     , mapping()
     , property_pool(std::make_unique<PropertyPool<dim, spacedim>>(0))
     , particles()
-    , ghost_particles()
     , global_number_of_particles(0)
     , local_number_of_particles(0)
     , global_max_particles_per_cell(0)
@@ -116,7 +115,6 @@ namespace Particles
     , mapping(&mapping, typeid(*this).name())
     , property_pool(std::make_unique<PropertyPool<dim, spacedim>>(n_properties))
     , particles(triangulation.n_active_cells())
-    , ghost_particles(triangulation.n_active_cells())
     , global_number_of_particles(0)
     , local_number_of_particles(0)
     , global_max_particles_per_cell(0)
@@ -160,7 +158,6 @@ namespace Particles
                                                         new_mapping);
 
     particles.resize(triangulation->n_active_cells());
-    ghost_particles.resize(triangulation->n_active_cells());
   }
 
 
@@ -189,7 +186,6 @@ namespace Particles
       particle_handler.global_max_particles_per_cell;
     next_free_particle_index = particle_handler.next_free_particle_index;
     particles                = particle_handler.particles;
-    ghost_particles          = particle_handler.ghost_particles;
 
     ghost_particles_cache.ghost_particles_by_domain =
       particle_handler.ghost_particles_cache.ghost_particles_by_domain;
@@ -224,15 +220,6 @@ namespace Particles
     if (triangulation != nullptr)
       particles.resize(triangulation->n_active_cells());
 
-    for (auto &particles_in_cell : ghost_particles)
-      for (auto &particle : particles_in_cell)
-        if (particle != PropertyPool<dim, spacedim>::invalid_handle)
-          property_pool->deregister_particle(particle);
-
-    ghost_particles.clear();
-    if (triangulation != nullptr)
-      ghost_particles.resize(triangulation->n_active_cells());
-
     // the particle properties have already been deleted by their destructor,
     // but the memory is still allocated. Return the memory as well.
     property_pool->clear();
@@ -314,10 +301,7 @@ namespace Particles
 
     if (cell->is_artificial() == false)
       {
-        if (cell->is_locally_owned() == true)
-          return particles[cell->active_cell_index()].size();
-        else
-          return ghost_particles[cell->active_cell_index()].size();
+        return particles[cell->active_cell_index()].size();
       }
     else
       AssertThrow(false,
@@ -351,8 +335,7 @@ namespace Particles
 
     if (cell->is_artificial() == false)
       {
-        particle_container &container =
-          (cell->is_locally_owned() == true) ? particles : ghost_particles;
+        particle_container &container = particles;
 
         if (container[active_cell_index].size() == 0)
           {
@@ -1328,15 +1311,6 @@ namespace Particles
 #ifndef DEAL_II_WITH_MPI
     (void)enable_cache;
 #else
-    // First clear the current ghost_particle information
-    for (auto &particles_in_cell : ghost_particles)
-      for (auto &particle : particles_in_cell)
-        if (particle != PropertyPool<dim, spacedim>::invalid_handle)
-          property_pool->deregister_particle(particle);
-    ghost_particles.clear();
-    ghost_particles.resize(triangulation->n_active_cells());
-
-    // Clear ghost particles data structures and invalidate cache
     ghost_particles_cache.ghost_particles_by_domain.clear();
     ghost_particles_cache.valid = false;
 
@@ -1381,9 +1355,11 @@ namespace Particles
           }
       }
 
+
+    // TODO - To fix
     send_recv_particles(
       ghost_particles_cache.ghost_particles_by_domain,
-      ghost_particles,
+      particles,
       std::map<
         types::subdomain_id,
         std::vector<
@@ -1418,7 +1394,7 @@ namespace Particles
 
 
     send_recv_particles_properties_and_location(
-      ghost_particles_cache.ghost_particles_by_domain, ghost_particles);
+      ghost_particles_cache.ghost_particles_by_domain, particles);
 #endif
   }