First working version

diff --git a/include/deal.II/grid/grid_tools_cache.h b/include/deal.II/grid/grid_tools_cache.h
index b926b5808bcf8545ab4bf8070be51c780e719d07..4d5ccfcf75dc0feb086fb6463d10ec0abc353e5f 100644 (file)
--- a/include/deal.II/grid/grid_tools_cache.h
+++ b/include/deal.II/grid/grid_tools_cache.h
@@ -75,7 +75,7 @@ namespace GridTools
      * @param mapping The mapping to use when computing cached objects
      */
     Cache(const Triangulation<dim, spacedim> &tria,
-          const Mapping<dim, spacedim> &      mapping =
+          const Mapping<dim, spacedim>       &mapping =
             (ReferenceCells::get_hypercube<dim>()
 #ifndef _MSC_VER
                .template get_default_linear_mapping<dim, spacedim>()
@@ -168,6 +168,15 @@ namespace GridTools
     const std::vector<std::set<unsigned int>> &
     get_vertex_to_neighbor_subdomain() const;
 
+    /**
+     * TODO
+     */
+    const std::map<unsigned int, std::vector<unsigned int>> &
+    get_coinciding_vertex_groups() const;
+
+    const std::map<unsigned int, unsigned int> &
+    get_vertex_to_coinciding_vertex_group() const;
+
     /**
      * Return a reference to the stored triangulation.
      */
@@ -292,6 +301,20 @@ namespace GridTools
      */
     mutable std::vector<std::set<unsigned int>> vertex_to_neighbor_subdomain;
 
+    /**
+     * Store an std::map of std::vector of unsigned integer containing
+     * coinciding vertices labeled by an arbitrary element from them
+     */
+    mutable std::map<unsigned int, std::vector<unsigned int>>
+      coinciding_vertex_groups;
+
+    /**
+     * Store an std::map of unsigned integer containing
+     * the label of a group of coinciding vertices
+     */
+    mutable std::map<unsigned int, unsigned int>
+      vertex_to_coinciding_vertex_group;
+
     /**
      * 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 904073a7540374072669b0e2c0d50fb42a92e82f..381c3543b1e30dce5ea04698cc8083ff210db921 100644 (file)
--- a/include/deal.II/grid/grid_tools_cache_update_flags.h
+++ b/include/deal.II/grid/grid_tools_cache_update_flags.h
@@ -84,6 +84,12 @@ namespace GridTools
      */
     update_vertex_to_neighbor_subdomain = 0x100,
 
+    /**
+     * Update the coinciding vertex groups as well as the
+     * vertex_to_coinciding_vertex_group
+     */
+    update_collection_of_coinciding_vertices = 0x200,
+
     /**
      * Update all objects.
      */

diff --git a/source/grid/grid_tools_cache.cc b/source/grid/grid_tools_cache.cc
index 76bc5507c5a1635e0fec69d928db8c0004260d70..801df23fb10bac9888fc531cb1dcd94a0e5f1bab 100644 (file)
--- a/source/grid/grid_tools_cache.cc
+++ b/source/grid/grid_tools_cache.cc
@@ -26,7 +26,7 @@ namespace GridTools
 {
   template <int dim, int spacedim>
   Cache<dim, spacedim>::Cache(const Triangulation<dim, spacedim> &tria,
-                              const Mapping<dim, spacedim> &      mapping)
+                              const Mapping<dim, spacedim>       &mapping)
     : update_flags(update_all)
     , tria(&tria)
     , mapping(&mapping)
@@ -218,6 +218,36 @@ namespace GridTools
     return vertex_to_neighbor_subdomain;
   }
 
+  template <int dim, int spacedim>
+  const std::map<unsigned int, std::vector<unsigned int>> &
+  Cache<dim, spacedim>::get_coinciding_vertex_groups() const
+  {
+    if (update_flags & update_collection_of_coinciding_vertices)
+      {
+        GridTools::collect_coinciding_vertices(
+          *tria, coinciding_vertex_groups, vertex_to_coinciding_vertex_group);
+
+        update_flags = update_flags & ~update_collection_of_coinciding_vertices;
+      }
+
+    return coinciding_vertex_groups;
+  }
+
+  template <int dim, int spacedim>
+  const std::map<unsigned int, unsigned int> &
+  Cache<dim, spacedim>::get_vertex_to_coinciding_vertex_group() const
+  {
+    if (update_flags & update_collection_of_coinciding_vertices)
+      {
+        GridTools::collect_coinciding_vertices(
+          *tria, coinciding_vertex_groups, vertex_to_coinciding_vertex_group);
+
+        update_flags = update_flags & ~update_collection_of_coinciding_vertices;
+      }
+
+    return vertex_to_coinciding_vertex_group;
+  }
+
 #include "grid_tools_cache.inst"
 
 } // namespace GridTools

diff --git a/source/particles/particle_handler.cc b/source/particles/particle_handler.cc
index 196ba13565b174ed2cb0eacf7a696d87516e701c..fae3bc91bd1ea40d12f9702e56fc5da3dbb2f915 100644 (file)
--- a/source/particles/particle_handler.cc
+++ b/source/particles/particle_handler.cc
@@ -77,7 +77,7 @@ namespace Particles
   template <int dim, int spacedim>
   ParticleHandler<dim, spacedim>::ParticleHandler(
     const Triangulation<dim, spacedim> &triangulation,
-    const Mapping<dim, spacedim> &      mapping,
+    const Mapping<dim, spacedim>       &mapping,
     const unsigned int                  n_properties)
     : triangulation(&triangulation, typeid(*this).name())
     , mapping(&mapping, typeid(*this).name())
@@ -117,7 +117,7 @@ namespace Particles
   void
   ParticleHandler<dim, spacedim>::initialize(
     const Triangulation<dim, spacedim> &new_triangulation,
-    const Mapping<dim, spacedim> &      new_mapping,
+    const Mapping<dim, spacedim>       &new_mapping,
     const unsigned int                  n_properties)
   {
     clear();
@@ -577,7 +577,7 @@ namespace Particles
   template <int dim, int spacedim>
   typename ParticleHandler<dim, spacedim>::particle_iterator
   ParticleHandler<dim, spacedim>::insert_particle(
-    const Particle<dim, spacedim> &                                    particle,
+    const Particle<dim, spacedim>                                     &particle,
     const typename Triangulation<dim, spacedim>::active_cell_iterator &cell)
   {
     return insert_particle(particle.get_location(),
@@ -623,7 +623,7 @@ namespace Particles
   template <int dim, int spacedim>
   typename ParticleHandler<dim, spacedim>::particle_iterator
   ParticleHandler<dim, spacedim>::insert_particle(
-    const void *&                                                      data,
+    const void                                                       *&data,
     const typename Triangulation<dim, spacedim>::active_cell_iterator &cell)
   {
     Assert(triangulation != nullptr, ExcInternalError());
@@ -648,8 +648,8 @@ namespace Particles
   template <int dim, int spacedim>
   typename ParticleHandler<dim, spacedim>::particle_iterator
   ParticleHandler<dim, spacedim>::insert_particle(
-    const Point<spacedim> &     position,
-    const Point<dim> &          reference_position,
+    const Point<spacedim>      &position,
+    const Point<dim>           &reference_position,
     const types::particle_index particle_index,
     const typename Triangulation<dim, spacedim>::active_cell_iterator &cell,
     const ArrayView<const double> &properties)
@@ -764,8 +764,8 @@ namespace Particles
   ParticleHandler<dim, spacedim>::insert_global_particles(
     const std::vector<Point<spacedim>> &positions,
     const std::vector<std::vector<BoundingBox<spacedim>>>
-      &                                       global_bounding_boxes,
-    const std::vector<std::vector<double>> &  properties,
+                                             &global_bounding_boxes,
+    const std::vector<std::vector<double>>   &properties,
     const std::vector<types::particle_index> &ids)
   {
     if (!properties.empty())
@@ -1190,7 +1190,7 @@ namespace Particles
     compare_particle_association(
       const unsigned int                 a,
       const unsigned int                 b,
-      const Tensor<1, dim> &             particle_direction,
+      const Tensor<1, dim>              &particle_direction,
       const std::vector<Tensor<1, dim>> &center_directions)
     {
       const double scalar_product_a = center_directions[a] * particle_direction;
@@ -1542,6 +1542,18 @@ namespace Particles
     ghost_particles_cache.ghost_particles_by_domain.clear();
     ghost_particles_cache.valid = false;
 
+    // In the case of a parallel simulation with periodic boundary conditions
+    // the vertices associated with periodic boundaries are not directly
+    // connected to the ghost cells but they are connected to the ghost cells
+    // through their coinciding vertices. We loop over the coinciding vertices
+    // to identify cells in which ghost particles must be generated
+    const std::map<unsigned int, std::vector<unsigned int>>
+      coinciding_vertex_groups =
+        triangulation_cache->get_coinciding_vertex_groups();
+    const std::map<unsigned int, unsigned int>
+      vertex_to_coinciding_vertex_group =
+        triangulation_cache->get_vertex_to_coinciding_vertex_group();
+
     const std::set<types::subdomain_id> ghost_owners =
       parallel_triangulation->ghost_owners();
     for (const auto ghost_owner : ghost_owners)
@@ -1561,6 +1573,26 @@ namespace Particles
                 cell_to_neighbor_subdomain.insert(
                   vertex_to_neighbor_subdomain[cell->vertex_index(v)].begin(),
                   vertex_to_neighbor_subdomain[cell->vertex_index(v)].end());
+
+                // If the vertex has one or more coinciding vertices, then
+                // all of these vertices, including the current vertex being
+                // looped over are contained inside the
+                // vertex_to_coinciding_vertex_group Consequently, we loop over
+                // the full content of the coinciding vertex group which will
+                // include the current vertex and we skip the current vertex.
+                if (vertex_to_coinciding_vertex_group.find(cell->vertex_index(
+                      v)) != vertex_to_coinciding_vertex_group.end())
+                  {
+                    for (const unsigned int cv : coinciding_vertex_groups.at(
+                           vertex_to_coinciding_vertex_group.at(
+                             cell->vertex_index(v))))
+                      {
+                        if (cv != cell->vertex_index(v))
+                          cell_to_neighbor_subdomain.insert(
+                            vertex_to_neighbor_subdomain[cv].begin(),
+                            vertex_to_neighbor_subdomain[cv].end());
+                      }
+                  }
               }
 
             if (cell_to_neighbor_subdomain.size() > 0)
@@ -1581,6 +1613,8 @@ namespace Particles
           }
       }
 
+
+
     send_recv_particles(
       ghost_particles_cache.ghost_particles_by_domain,
       std::map<
@@ -1634,7 +1668,7 @@ namespace Particles
     const std::map<
       types::subdomain_id,
       std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>>
-      &        send_cells,
+              &send_cells,
     const bool build_cache)
   {
     Assert(triangulation != nullptr, ExcInternalError());
@@ -2038,7 +2072,7 @@ namespace Particles
   template <int dim, int spacedim>
   void
   ParticleHandler<dim, spacedim>::register_additional_store_load_functions(
-    const std::function<std::size_t()> &                            size_callb,
+    const std::function<std::size_t()>                             &size_callb,
     const std::function<void *(const particle_iterator &, void *)> &store_callb,
     const std::function<const void *(const particle_iterator &, const void *)>
       &load_callb)
@@ -2168,7 +2202,7 @@ namespace Particles
     const auto callback_function =
       [this](
         const typename Triangulation<dim, spacedim>::cell_iterator
-          &                                                     cell_iterator,
+                                                               &cell_iterator,
         const typename Triangulation<dim, spacedim>::CellStatus cell_status) {
         return this->pack_callback(cell_iterator, cell_status);
       };
@@ -2325,7 +2359,7 @@ namespace Particles
   template <int dim, int spacedim>
   void
   ParticleHandler<dim, spacedim>::unpack_callback(
-    const typename Triangulation<dim, spacedim>::cell_iterator &    cell,
+    const typename Triangulation<dim, spacedim>::cell_iterator     &cell,
     const typename Triangulation<dim, spacedim>::CellStatus         status,
     const boost::iterator_range<std::vector<char>::const_iterator> &data_range)
   {
@@ -2359,14 +2393,12 @@ namespace Particles
     // now update particle storage location and properties if necessary
     switch (status)
       {
-        case parallel::TriangulationBase<dim, spacedim>::CELL_PERSIST:
-          {
+          case parallel::TriangulationBase<dim, spacedim>::CELL_PERSIST: {
             // all particles are correctly inserted
           }
           break;
 
-        case parallel::TriangulationBase<dim, spacedim>::CELL_COARSEN:
-          {
+          case parallel::TriangulationBase<dim, spacedim>::CELL_COARSEN: {
             // all particles are in correct cell, but their reference location
             // has changed
             for (auto &particle : loaded_particles_on_cell)
@@ -2379,8 +2411,7 @@ namespace Particles
           }
           break;
 
-        case parallel::TriangulationBase<dim, spacedim>::CELL_REFINE:
-          {
+          case parallel::TriangulationBase<dim, spacedim>::CELL_REFINE: {
             // we need to find the correct child to store the particles and
             // their reference location has changed
             typename particle_container::iterator &cache =