From 7878df99e97c19cd1e6a0e9c785b489e70981cb3 Mon Sep 17 00:00:00 2001
From: Daniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Date: Tue, 22 Mar 2016 17:44:45 +0100
Subject: [PATCH] Use periodic_face_map also for the multigrid structures in
 distributed/tria.h, compare

---
 include/deal.II/distributed/tria.h |   4 +-
 include/deal.II/grid/tria.h        |  33 +++--
 source/distributed/tria.cc         | 226 ++++++++++++++++++++++++++---
 source/grid/tria.cc                |  17 ++-
 tests/mpi/periodicity_04.cc        |   7 +-
 5 files changed, 244 insertions(+), 43 deletions(-)

diff --git a/include/deal.II/distributed/tria.h b/include/deal.II/distributed/tria.h
index 4ed9440966..3a12ebf34f 100644
--- a/include/deal.II/distributed/tria.h
+++ b/include/deal.II/distributed/tria.h
@@ -946,7 +946,7 @@ namespace parallel
        */
       void
       fill_level_vertices_with_ghost_neighbors
-      (const unsigned int level,
+      (const int level,
        std::map<unsigned int, std::set<dealii::types::subdomain_id> >
        &vertices_with_ghost_neighbors);
 
@@ -960,7 +960,7 @@ namespace parallel
        * Used by DoFHandler::Policy::ParallelDistributed.
        */
       std::vector<bool>
-      mark_locally_active_vertices_on_level(const unsigned int level) const;
+      mark_locally_active_vertices_on_level(const int level) const;
 
       template <int, int> friend class dealii::internal::DoFHandler::Policy::ParallelDistributed;
     };
diff --git a/include/deal.II/grid/tria.h b/include/deal.II/grid/tria.h
index a20531642d..7ecbe2acf9 100644
--- a/include/deal.II/grid/tria.h
+++ b/include/deal.II/grid/tria.h
@@ -1471,19 +1471,6 @@ public:
     distorted_cells;
   };
 
-  /**
-   * If add_periodicity() is called, this variable stores the given
-   * periodic face pairs on level 0 for later access during the
-   * identification of ghost cells for the multigrid hierarchy and for
-   * setting up the periodic_face_map.
-   */
-  std::vector<GridTools::PeriodicFacePair<cell_iterator> > periodic_face_pairs_level_0;
-
-  /**
-   * If add_periodicity() is called, this variable stores the active periodic face pairs.
-   */
-  std::map<std::pair<cell_iterator, unsigned int>, std::pair<std::pair<cell_iterator, unsigned int>, std::bitset<3> > > periodic_face_map;
-
   /**
    * Make the dimension available in function templates.
    */
@@ -2984,7 +2971,6 @@ public:
     * Declare the (coarse) face pairs given in the argument of this function
     * as periodic. This way it it possible to obtain neighbors across periodic
     * boundaries.
-    * This function initializes periodic_face_map that stores the active periodic faces.
     *
     * The vector can be filled by the function
     * GridTools::collect_periodic_faces.
@@ -3000,6 +2986,12 @@ public:
   add_periodicity
   (const std::vector<GridTools::PeriodicFacePair<cell_iterator> > &);
 
+  /**
+    * Return the periodic_face_map.
+    */
+  const std::map<std::pair<cell_iterator, unsigned int>,std::pair<std::pair<cell_iterator, unsigned int>, std::bitset<3> > > &
+  get_periodic_face_map() const;
+
 
   BOOST_SERIALIZATION_SPLIT_MEMBER()
 
@@ -3113,6 +3105,19 @@ protected:
 
 
 private:
+  /**
+    * If add_periodicity() is called, this variable stores the given
+    * periodic face pairs on level 0 for later access during the
+    * identification of ghost cells for the multigrid hierarchy and for
+    * setting up the periodic_face_map.
+    */
+  std::vector<GridTools::PeriodicFacePair<cell_iterator> > periodic_face_pairs_level_0;
+
+  /**
+   * If add_periodicity() is called, this variable stores the active periodic face pairs.
+   */
+  std::map<std::pair<cell_iterator, unsigned int>, std::pair<std::pair<cell_iterator, unsigned int>, std::bitset<3> > > periodic_face_map;
+
   /**
    * @name Cell iterator functions for internal use
    * @{
diff --git a/source/distributed/tria.cc b/source/distributed/tria.cc
index 0bdab3a7b8..2e052d49e2 100644
--- a/source/distributed/tria.cc
+++ b/source/distributed/tria.cc
@@ -2244,7 +2244,7 @@ namespace parallel
         }
 
       this->update_number_cache ();
-      Triangulation<dim, spacedim>::update_periodic_face_map();
+      this->update_periodic_face_map();
     }
 
 
@@ -2773,7 +2773,7 @@ namespace parallel
         }
 
       this->update_number_cache ();
-      Triangulation<dim, spacedim>::update_periodic_face_map();
+      this->update_periodic_face_map();
     }
 
 
@@ -3197,6 +3197,13 @@ namespace parallel
             topological_vertex_numbering[periodic.cell[0]->face(periodic.face_idx[0])->vertex_index(vface0)]
               = topological_vertex_numbering[periodic.cell[1]->face(periodic.face_idx[1])->vertex_index(v)]
                 = min_index;
+
+            /*  std::cout << "Old identified vertices " << periodic.cell[0]->face(periodic.face_idx[0])->vertex_index(vface0)
+                        << " at "                 << periodic.cell[0]->face(periodic.face_idx[0])->vertex(vface0)
+                        << " and "                << periodic.cell[1]->face(periodic.face_idx[1])->vertex_index(v)
+                        << " at "                 << periodic.cell[1]->face(periodic.face_idx[1])->vertex(v)
+                        << std::endl;*/
+
           }
       }
 
@@ -3206,10 +3213,9 @@ namespace parallel
       // artificial cell layer (the active cells are taken care of by p4est)
       template <int dim, int spacedim>
       bool enforce_mesh_balance_over_periodic_boundaries
-      (Triangulation<dim,spacedim> &tria,
-       const std::vector<GridTools::PeriodicFacePair<typename dealii::Triangulation<dim,spacedim>::cell_iterator> > &periodic_face_pairs_level_0)
+      (Triangulation<dim,spacedim> &tria)
       {
-        if (periodic_face_pairs_level_0.empty())
+        if (tria.get_periodic_face_map().size()==0)
           return false;
 
         std::vector<bool> flags_before[2];
@@ -3217,14 +3223,100 @@ namespace parallel
         tria.save_refine_flags (flags_before[1]);
 
         std::vector<unsigned int> topological_vertex_numbering(tria.n_vertices());
+        std::vector<unsigned int> topological_vertex_numbering_cmp(tria.n_vertices());
         for (unsigned int i=0; i<topological_vertex_numbering.size(); ++i)
-          topological_vertex_numbering[i] = i;
+          {
+            topological_vertex_numbering[i] = i;
+            topological_vertex_numbering_cmp[i] = i;
+          }
+        //std::cout << std::endl;
         for (unsigned int i=0; i<periodic_face_pairs_level_0.size(); ++i)
           {
             identify_periodic_vertices_recursively(periodic_face_pairs_level_0[i],
-                                                   topological_vertex_numbering);
+                                                   topological_vertex_numbering_cmp);
+          }
+        // combine vertices that have different locations (and thus, different
+        // vertex_index) but represent the same topological entity over periodic
+        // boundaries. The vector topological_vertex_numbering contains a linear
+        // map from 0 to n_vertices at input and at output relates periodic
+        // vertices with only one vertex index. The output is used to always
+        // identify the same vertex according to the periodicity, e.g. when
+        // finding the maximum cell level around a vertex.
+        //
+        // Example: On a 3D cell with vertices numbered from 0 to 7 and periodic
+        // boundary conditions in x direction, the vector
+        // topological_vertex_numbering will contain the numbers
+        // {0,0,2,2,4,4,6,6} (because the vertex pairs {0,1}, {2,3}, {4,5},
+        // {6,7} belong together, respectively). If periodicity is set in x and
+        // z direction, the output is {0,0,2,2,0,0,2,2}, and if periodicity is
+        // in all directions, the output is simply {0,0,0,0,0,0,0,0}.
+        typedef typename Triangulation<dim, spacedim>::cell_iterator cell_iterator;
+        typename std::map<std::pair<cell_iterator, unsigned int>,
+                 std::pair<std::pair<cell_iterator,unsigned int>, std::bitset<3> > >::const_iterator it;
+        for (it = tria.get_periodic_face_map().begin(); it!= tria.get_periodic_face_map().end(); ++it)
+          {
+            const cell_iterator &cell_1 = it->first.first;
+            const unsigned int face_no_1 = it->first.second;
+            const cell_iterator &cell_2 = it->second.first.first;
+            const unsigned int face_no_2 = it->second.first.second;
+            const std::bitset<3> face_orientation = it->second.second;
+
+            if (cell_1->level() == cell_2->level())
+              {
+                for (unsigned int v=0; v<GeometryInfo<dim-1>::vertices_per_cell; ++v)
+                  {
+                    // take possible non-standard orientation of face on cell[0] into
+                    // account
+                    const unsigned int vface0 =
+                      GeometryInfo<dim>::standard_to_real_face_vertex(v,face_orientation[0],
+                                                                      face_orientation[1],
+                                                                      face_orientation[2]);
+                    const unsigned int vi0 = topological_vertex_numbering[cell_1->face(face_no_1)->vertex_index(vface0)];
+                    const unsigned int vi1 = topological_vertex_numbering[cell_2->face(face_no_2)->vertex_index(v)];
+                    const unsigned int min_index = std::min(vi0, vi1);
+                    topological_vertex_numbering[cell_1->face(face_no_1)->vertex_index(vface0)]
+                      = topological_vertex_numbering[cell_2->face(face_no_2)->vertex_index(v)]
+                        = min_index;
+                    /*  std::cout << "Newly identified vertices " << cell_1->face(face_no_1)->vertex_index(vface0)
+                              << " at "                 << cell_1->face(face_no_1)->vertex(vface0)
+                              << " and "                << cell_2->face(face_no_2)->vertex_index(v)
+                              << " at "                 << cell_2->face(face_no_2)->vertex(v)
+                              << std::endl;*/
+                  }
+              }
+          }
+
+        // There must not be any chains!
+        for (unsigned int i=0; i<topological_vertex_numbering.size(); ++i)
+          {
+            const unsigned int j = topological_vertex_numbering[i];
+            if (j != i)
+              Assert(topological_vertex_numbering[j] == j,
+                     ExcInternalError());
           }
 
+        for (unsigned int i=0; i<topological_vertex_numbering_cmp.size(); ++i)
+          {
+            const unsigned int j = topological_vertex_numbering_cmp[i];
+            if (j != i)
+              Assert(topological_vertex_numbering_cmp[j] == j,
+                     ExcInternalError());
+          }
+
+        if (topological_vertex_numbering != topological_vertex_numbering_cmp)
+          {
+            for (unsigned int i=0; i<topological_vertex_numbering.size(); ++i)
+              {
+                if (topological_vertex_numbering[i] != topological_vertex_numbering_cmp[i])
+                  std::cout << i
+                            << " old: " << topological_vertex_numbering_cmp[i]
+                            << " new: " << topological_vertex_numbering[i]
+                            << std::endl;
+              }
+            Assert( false, ExcInternalError());
+          }
+
+
         // this code is replicated from grid/tria.cc but using an indirection
         // for periodic boundary conditions
         bool continue_iterating = true;
@@ -3356,12 +3448,11 @@ namespace parallel
       do
         {
           this->dealii::Triangulation<dim,spacedim>::prepare_coarsening_and_refinement();
-
+          this->update_periodic_face_map();
           // enforce 2:1 mesh balance over periodic boundaries
           if (this->smooth_grid &
               dealii::Triangulation<dim,spacedim>::limit_level_difference_at_vertices)
-            mesh_changed = enforce_mesh_balance_over_periodic_boundaries(*this,
-                           this->periodic_face_pairs_level_0);
+            mesh_changed = enforce_mesh_balance_over_periodic_boundaries(*this);
         }
       while (mesh_changed);
 
@@ -3895,7 +3986,7 @@ namespace parallel
 
       refinement_in_progress = false;
       this->update_number_cache ();
-      Triangulation<dim, spacedim>::update_periodic_face_map();
+      this->update_periodic_face_map();
     }
 
     template <int dim, int spacedim>
@@ -3965,7 +4056,7 @@ namespace parallel
 
       // update how many cells, edges, etc, we store locally
       this->update_number_cache ();
-      Triangulation<dim, spacedim>::update_periodic_face_map();
+      this->update_periodic_face_map();
     }
 
 
@@ -4678,7 +4769,7 @@ namespace parallel
     void
     Triangulation<dim,spacedim>::
     fill_level_vertices_with_ghost_neighbors
-    (const unsigned int level,
+    (const int level,
      std::map<unsigned int, std::set<dealii::types::subdomain_id> >
      &vertices_with_ghost_neighbors)
     {
@@ -4694,9 +4785,72 @@ namespace parallel
               vertices_with_ghost_neighbors[cell->vertex_index(v)]
               .insert (cell->level_subdomain_id());
 
+      std::map<unsigned int, std::set<dealii::types::subdomain_id> >
+      vertices_with_ghost_neighbors_cmp = vertices_with_ghost_neighbors;
+
+
       for (unsigned int i=0; i<this->periodic_face_pairs_level_0.size(); ++i)
         set_periodic_ghost_neighbors_recursively(this->periodic_face_pairs_level_0[i],
-                                                 level, vertices_with_ghost_neighbors);
+                                                 level, vertices_with_ghost_neighbors_cmp);
+
+      //now for the vertices on periodic faces
+      typename std::map<std::pair<cell_iterator, unsigned int>,
+               std::pair<std::pair<cell_iterator,unsigned int>, std::bitset<3> > >::const_iterator it;
+
+      for (it = this->get_periodic_face_map().begin(); it!= this->get_periodic_face_map().end(); ++it)
+        {
+          const cell_iterator &cell_1 = it->first.first;
+          const unsigned int face_no_1 = it->first.second;
+          const cell_iterator &cell_2 = it->second.first.first;
+          const unsigned int face_no_2 = it->second.first.second;
+          const std::bitset<3> face_orientation = it->second.second;
+
+          if (cell_1->level() == level &&
+              cell_2->level() == level)
+            {
+              for (unsigned int v=0; v<GeometryInfo<dim-1>::vertices_per_cell; ++v)
+                {
+                  // take possible non-standard orientation of faces into account
+                  const unsigned int vface0 =
+                    GeometryInfo<dim>::standard_to_real_face_vertex(v,face_orientation[0],
+                                                                    face_orientation[1],
+                                                                    face_orientation[2]);
+                  const unsigned int idx0 = cell_1->face(face_no_1)->vertex_index(vface0);
+                  const unsigned int idx1 = cell_2->face(face_no_2)->vertex_index(v);
+                  if (vertices_with_ghost_neighbors.find(idx0) != vertices_with_ghost_neighbors.end())
+                    vertices_with_ghost_neighbors[idx1].insert(vertices_with_ghost_neighbors[idx0].begin(),
+                                                               vertices_with_ghost_neighbors[idx0].end());
+                  if (vertices_with_ghost_neighbors.find(idx1) != vertices_with_ghost_neighbors.end())
+                    vertices_with_ghost_neighbors[idx0].insert(vertices_with_ghost_neighbors[idx1].begin(),
+                                                               vertices_with_ghost_neighbors[idx1].end());
+                }
+            }
+        }
+
+      if (vertices_with_ghost_neighbors != vertices_with_ghost_neighbors_cmp)
+        {
+          std::cout << "Print vertices_with_ghost_neighbors: " << std::endl;
+          std::map<unsigned int, std::set<dealii::types::subdomain_id> >::const_iterator it;
+          for (it = vertices_with_ghost_neighbors.begin(); it!=vertices_with_ghost_neighbors.end(); ++it)
+            {
+              std::cout << "Vertex " << it->first << " is on subdomains ";
+              std::set<dealii::types::subdomain_id>::const_iterator it2;
+              for (it2 = it->second.begin(); it2!=it->second.end(); ++it2)
+                std::cout << *it2 << " ";
+              std::cout << std::endl;
+            }
+
+          std::cout << "Print vertices_with_ghost_neighbors_cmp: " << std::endl;
+          for (it = vertices_with_ghost_neighbors_cmp.begin(); it!=vertices_with_ghost_neighbors_cmp.end(); ++it)
+            {
+              std::cout << "Vertex " << it->first << " is on subdomains ";
+              std::set<dealii::types::subdomain_id>::const_iterator it2;
+              for (it2 = it->second.begin(); it2!=it->second.end(); ++it2)
+                std::cout << *it2 << " ";
+              std::cout << std::endl;
+            }
+          Assert(false, ExcInternalError());
+        }
     }
 
 
@@ -4704,7 +4858,7 @@ namespace parallel
     template<int dim, int spacedim>
     std::vector<bool>
     Triangulation<dim,spacedim>
-    ::mark_locally_active_vertices_on_level (const unsigned int level) const
+    ::mark_locally_active_vertices_on_level (const int level) const
     {
       Assert (dim>1, ExcNotImplemented());
 
@@ -4716,10 +4870,48 @@ namespace parallel
           for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
             marked_vertices[cell->vertex_index(v)] = true;
 
+      std::vector<bool> marked_vertices_cmp = marked_vertices;
+
       for (unsigned int i=0; i<this->periodic_face_pairs_level_0.size(); ++i)
         mark_periodic_vertices_recursively(this->periodic_face_pairs_level_0[i],
-                                           level, marked_vertices);
+                                           level, marked_vertices_cmp);
+
+      /**
+       * ensure that if one of the two vertices on a periodic face is marked
+       * as active (i.e., belonging to an owned level cell), also the other
+       * one is active
+       */
+      typename std::map<std::pair<cell_iterator, unsigned int>,
+               std::pair<std::pair<cell_iterator,unsigned int>, std::bitset<3> > >::const_iterator it;
+
+      for (it = this->get_periodic_face_map().begin(); it!= this->get_periodic_face_map().end(); ++it)
+        {
+          const cell_iterator &cell_1 = it->first.first;
+          const unsigned int face_no_1 = it->first.second;
+          const cell_iterator &cell_2 = it->second.first.first;
+          const unsigned int face_no_2 = it->second.first.second;
+          const std::bitset<3> &face_orientation = it->second.second;
+
+          if (cell_1->level() == level &&
+              cell_2->level() == level)
+            {
+              for (unsigned int v=0; v<GeometryInfo<dim-1>::vertices_per_cell; ++v)
+                {
+                  // take possible non-standard orientation of faces into account
+                  const unsigned int vface0 =
+                    GeometryInfo<dim>::standard_to_real_face_vertex(v,face_orientation[0],
+                                                                    face_orientation[1],
+                                                                    face_orientation[2]);
+                  if (marked_vertices[cell_1->face(face_no_1)->vertex_index(vface0)] ||
+                      marked_vertices[cell_2->face(face_no_2)->vertex_index(v)])
+                    marked_vertices[cell_1->face(face_no_1)->vertex_index(vface0)]
+                      = marked_vertices[cell_2->face(face_no_2)->vertex_index(v)]
+                        = true;
+                }
+            }
+        }
 
+      Assert(marked_vertices == marked_vertices_cmp, ExcInternalError());
       return marked_vertices;
     }
 
@@ -4965,7 +5157,7 @@ namespace parallel
         }
 
       this->update_number_cache ();
-      Triangulation<dim, spacedim>::update_periodic_face_map();
+      this->update_periodic_face_map();
     }
 
 
diff --git a/source/grid/tria.cc b/source/grid/tria.cc
index 4e06ab94fe..99a3e47e25 100644
--- a/source/grid/tria.cc
+++ b/source/grid/tria.cc
@@ -11744,6 +11744,14 @@ Triangulation<dim, spacedim>::add_periodicity
   update_periodic_face_map();
 }
 
+template <int dim, int spacedim>
+const typename std::map<std::pair<typename Triangulation<dim, spacedim>::cell_iterator, unsigned int>,
+      std::pair<std::pair<typename Triangulation<dim, spacedim>::cell_iterator, unsigned int>, std::bitset<3> > > &
+      Triangulation<dim, spacedim>::get_periodic_face_map() const
+{
+  return periodic_face_map;
+}
+
 
 template <int dim, int spacedim>
 void
@@ -11784,12 +11792,7 @@ Triangulation<dim, spacedim>::execute_coarsening_and_refinement ()
   AssertThrow (cells_with_distorted_children.distorted_cells.size() == 0,
                cells_with_distorted_children);
 
-  // For parallel::distributed::Triangulations we update
-  // periodic_face_map later.
-  const parallel::Triangulation< dim, spacedim > *distributed_triangulation
-    = dynamic_cast<const parallel::Triangulation< dim, spacedim > *> (this);
-  if (!distributed_triangulation)
-    update_periodic_face_map();
+  update_periodic_face_map();
 }
 
 
@@ -11819,7 +11822,7 @@ Triangulation<dim,spacedim>::update_periodic_face_map ()
   //first empty the currently stored objects
   periodic_face_map.clear();
 
-  typename std::vector<GridTools::PeriodicFacePair<typename Triangulation<dim,spacedim>::cell_iterator> >::const_iterator it;
+  typename std::vector<GridTools::PeriodicFacePair<cell_iterator> >::const_iterator it;
   for (it=periodic_face_pairs_level_0.begin(); it!=periodic_face_pairs_level_0.end(); ++it)
     {
       update_periodic_face_map_recursively<dim, spacedim>
diff --git a/tests/mpi/periodicity_04.cc b/tests/mpi/periodicity_04.cc
index 7211af7135..844355036c 100644
--- a/tests/mpi/periodicity_04.cc
+++ b/tests/mpi/periodicity_04.cc
@@ -241,7 +241,7 @@ void check
   triangulation.execute_coarsening_and_refinement();
 
   typedef std::pair<typename Triangulation<dim>::cell_iterator, unsigned int> CellFace;
-  const typename std::map<CellFace, std::pair<CellFace, std::bitset<3> > > &face_map = triangulation.periodic_face_map;
+  const typename std::map<CellFace, std::pair<CellFace, std::bitset<3> > > &face_map = triangulation.get_periodic_face_map();
   typename std::map<CellFace, std::pair<CellFace, std::bitset<3> > >::const_iterator it;
   int sum_of_pairs_local = face_map.size();
   int sum_of_pairs_global;
@@ -263,8 +263,9 @@ void check
             {
               std::cout << "face_center_1: " << face_center_1 << std::endl;
               std::cout << "face_center_2: " << face_center_2 << std::endl;
-              for (auto it = triangulation.periodic_face_map.begin();
-                   it !=  triangulation.periodic_face_map.end(); ++it)
+              typename std::map<CellFace, std::pair<CellFace, std::bitset<3> > >::const_iterator it;
+              for (it = triangulation.get_periodic_face_map().begin();
+                   it !=  triangulation.get_periodic_face_map().end(); ++it)
                 {
                   std::cout << "The cell with center " << it->first.first->center()
                             << " has on face " << it->first.second
-- 
2.39.5