Implement updating periodic_face_map

diff --git a/include/deal.II/distributed/tria.h b/include/deal.II/distributed/tria.h
index 95e5859e657439d4322e2d85bac7832c307aa150..4ed944096617b9a3009ddc8fb0579233914dab6d 100644 (file)
--- a/include/deal.II/distributed/tria.h
+++ b/include/deal.II/distributed/tria.h
@@ -878,13 +878,6 @@ namespace parallel
       std::vector<types::global_dof_index> coarse_cell_to_p4est_tree_permutation;
       std::vector<types::global_dof_index> p4est_tree_to_coarse_cell_permutation;
 
-      /**
-       * 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.
-       */
-      std::vector<GridTools::PeriodicFacePair<cell_iterator> > periodic_face_pairs_level_0;
-
       /**
        * Return a pointer to the p4est tree that belongs to the given
        * dealii_coarse_cell_index()

diff --git a/include/deal.II/grid/tria.h b/include/deal.II/grid/tria.h
index 0678a6abfb3961b6dd1ea7df3d4111a1081bfb89..bd6ece4cf14e38ae1dc7c5a9cd10ab5edbd19c46 100644 (file)
--- a/include/deal.II/grid/tria.h
+++ b/include/deal.II/grid/tria.h
@@ -46,6 +46,11 @@ template <int dim, int spacedim> class Boundary;
 template <int dim, int spacedim> class StraightBoundary;
 template <int dim, int spacedim> class Manifold;
 
+namespace GridTools
+{
+  template<typename CellIterator>  struct PeriodicFacePair;
+}
+
 template <int, int, int> class TriaAccessor;
 template <int spacedim> class TriaAccessor<0,1,spacedim>;
 template <int, int, int> class TriaAccessorBase;
@@ -1468,15 +1473,15 @@ public:
   /**
    * 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.
+   * 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.
-   * The variable is updated on refinement.
    */
-  std::map<GridTools::PeriodicFacePair<cell_iterator>, GridTools::PeriodicFacePair<cell_iterator> > periodic_face_map;
+  std::map<std::pair<cell_iterator, unsigned int>, std::pair<cell_iterator, unsigned int> > periodic_face_map;
 
   /**
    * Make the dimension available in function templates.
@@ -2978,7 +2983,7 @@ 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 ActivePeriodicFacePairs and on refinement this
+    * This function initializes periodic_face_map that stores the active periodic faces.
     *
     * The vector can be filled by the function
     * GridTools::collect_periodic_faces.
@@ -3100,6 +3105,13 @@ protected:
                                 const unsigned int       magic_number2,
                                 std::istream            &in);
 
+  /**
+   * Recreate information about periodic neighbors from periodic_face_pairs_level_0.
+   * Only faces that belong to locally owned or ghosted cells are stored.
+   */
+  void update_periodic_face_map ();
+
+
 private:
   /**
    * @name Cell iterator functions for internal use

diff --git a/source/distributed/tria.cc b/source/distributed/tria.cc
index 3423854a295c524f75275509190d61c18ab8ed1e..0bdab3a7b846d52a391b668cc6e414157f079ab2 100644 (file)
--- a/source/distributed/tria.cc
+++ b/source/distributed/tria.cc
@@ -2244,6 +2244,7 @@ namespace parallel
         }
 
       this->update_number_cache ();
+      Triangulation<dim, spacedim>::update_periodic_face_map();
     }
 
 
@@ -2582,8 +2583,6 @@ namespace parallel
       coarse_cell_to_p4est_tree_permutation.resize (0);
       p4est_tree_to_coarse_cell_permutation.resize (0);
 
-      periodic_face_pairs_level_0.clear();
-
       dealii::Triangulation<dim,spacedim>::clear ();
 
       this->update_number_cache ();
@@ -2774,6 +2773,7 @@ namespace parallel
         }
 
       this->update_number_cache ();
+      Triangulation<dim, spacedim>::update_periodic_face_map();
     }
 
 
@@ -3207,7 +3207,7 @@ namespace parallel
       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)
+       const std::vector<GridTools::PeriodicFacePair<typename dealii::Triangulation<dim,spacedim>::cell_iterator> > &periodic_face_pairs_level_0)
       {
         if (periodic_face_pairs_level_0.empty())
           return false;
@@ -3361,7 +3361,7 @@ namespace parallel
           if (this->smooth_grid &
               dealii::Triangulation<dim,spacedim>::limit_level_difference_at_vertices)
             mesh_changed = enforce_mesh_balance_over_periodic_boundaries(*this,
-                           periodic_face_pairs_level_0);
+                           this->periodic_face_pairs_level_0);
         }
       while (mesh_changed);
 
@@ -3894,8 +3894,8 @@ namespace parallel
 
 
       refinement_in_progress = false;
-
       this->update_number_cache ();
+      Triangulation<dim, spacedim>::update_periodic_face_map();
     }
 
     template <int dim, int spacedim>
@@ -3965,6 +3965,7 @@ namespace parallel
 
       // update how many cells, edges, etc, we store locally
       this->update_number_cache ();
+      Triangulation<dim, spacedim>::update_periodic_face_map();
     }
 
 
@@ -4693,8 +4694,8 @@ namespace parallel
               vertices_with_ghost_neighbors[cell->vertex_index(v)]
               .insert (cell->level_subdomain_id());
 
-      for (unsigned int i=0; i<periodic_face_pairs_level_0.size(); ++i)
-        set_periodic_ghost_neighbors_recursively(periodic_face_pairs_level_0[i],
+      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);
     }
 
@@ -4715,8 +4716,8 @@ namespace parallel
           for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
             marked_vertices[cell->vertex_index(v)] = true;
 
-      for (unsigned int i=0; i<periodic_face_pairs_level_0.size(); ++i)
-        mark_periodic_vertices_recursively(periodic_face_pairs_level_0[i],
+      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);
 
       return marked_vertices;
@@ -4862,7 +4863,7 @@ namespace parallel
         }
 
       //finally call the base class for storing the periodicity information
-      Triangulation<dim, spacedim>::add_periodicity(periodicity_vector);
+      dealii::Triangulation<dim, spacedim>::add_periodicity(periodicity_vector);
 #else
       Assert(false, ExcMessage ("Need p4est version >= 0.3.4.1!"));
 #endif
@@ -4964,6 +4965,7 @@ namespace parallel
         }
 
       this->update_number_cache ();
+      Triangulation<dim, spacedim>::update_periodic_face_map();
     }
 
 

diff --git a/source/grid/tria.cc b/source/grid/tria.cc
index 78c065079c5ca2036d05a02e5b0067e04522eecc..69852420dc142b17544734463dcb9fc8ab481fa2 100644 (file)
--- a/source/grid/tria.cc
+++ b/source/grid/tria.cc
@@ -1,4 +1,5 @@
 // ---------------------------------------------------------------------
+
 //
 // Copyright (C) 1998 - 2016 by the deal.II authors
 //
@@ -907,6 +908,158 @@ namespace
         }
   }
 
+
+  template <int dim, int spacedim>
+  void update_periodic_face_map_recursively
+  (const typename Triangulation<dim,spacedim>::cell_iterator &cell_1,
+   const typename Triangulation<dim,spacedim>::cell_iterator &cell_2,
+   unsigned int n_face_1, unsigned int n_face_2,
+   const bool face_orientation, const bool face_flip, const bool face_rotation,
+   typename std::map<std::pair<typename Triangulation<dim, spacedim>::cell_iterator, unsigned int>,
+   std::pair<typename Triangulation<dim, spacedim>::cell_iterator, unsigned int> > &periodic_face_map)
+  {
+    typedef typename Triangulation<dim, spacedim>::face_iterator FaceIterator;
+    const FaceIterator face_1 = cell_1->face(n_face_1);
+    const FaceIterator face_2 = cell_2->face(n_face_2);
+
+    Assert((dim != 1) || (face_orientation == true && face_flip == false && face_rotation == false),
+           ExcMessage ("The supplied orientation " "(face_orientation, face_flip, face_rotation) " "is invalid for 1D"));
+
+    Assert((dim != 2) || (face_orientation == true && face_rotation == false),
+           ExcMessage ("The supplied orientation " "(face_orientation, face_flip, face_rotation) " "is invalid for 2D"));
+
+    Assert(face_1 != face_2,
+           ExcMessage ("face_1 and face_2 are equal!"));
+
+    Assert(face_1->at_boundary()  &&face_2->at_boundary(),
+           ExcMessage ("Periodic faces must be on the boundary"));
+
+    // A lookup table on how to go through the child cells depending on the
+    // orientation:
+    // see Documentation of GeometryInfo for details
+
+    static const int lookup_table_2d[2][2] =
+      //               flip:
+    {
+      { 0, 1 },   // false
+      { 1, 0 }    // true
+    };
+
+    static const int lookup_table_3d[2][2][2][4] =
+      //                           orientation flip  rotation
+    {
+      { { { 0, 2, 1, 3 },       // false       false false
+          { 2, 3, 0, 1 }        // false       false true
+        },
+        { { 3, 1, 2, 0 },       // false       true  false
+          { 1, 0, 3, 2 }        // false       true  true
+        }
+      },
+      { { { 0, 1, 2, 3 },       // true        false false
+          { 1, 3, 0, 2 }        // true        false true
+        },
+        { { 3, 2, 1, 0 },       // true        true  false
+          { 2, 0, 3, 1 }        // true        true  true
+        }
+      }
+    };
+
+    if (cell_1->has_children())
+      {
+        if (cell_2->has_children())
+          {
+            // In the case that both faces have children, we loop over all
+            // children and apply update_periodic_face_map_recursively recursively:
+
+            Assert(face_1->n_children() == GeometryInfo<dim>::max_children_per_face
+                   && face_2->n_children() == GeometryInfo<dim>::max_children_per_face,
+                   ExcNotImplemented());
+
+            for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
+              {
+                // Lookup the index for the second face
+                unsigned int j = 0;
+                switch (dim)
+                  {
+                  case 2:
+                    j = lookup_table_2d[face_flip][i];
+                    break;
+                  case 3:
+                    j = lookup_table_3d[face_orientation][face_flip][face_rotation][i];
+                    break;
+                  default:
+                    AssertThrow(false, ExcNotImplemented());
+                  }
+
+                // find subcell ids that belong to the subface indices
+                unsigned int child_cell_1
+                  = GeometryInfo<dim>::child_cell_on_face(cell_1->refinement_case(), n_face_1, i, face_orientation,
+                                                          face_flip, face_rotation, face_1->refinement_case());
+                unsigned int child_cell_2
+                  = GeometryInfo<dim>::child_cell_on_face(cell_2->refinement_case(), n_face_2, j, face_orientation,
+                                                          face_flip, face_rotation, face_2->refinement_case());
+
+                // precondition: subcell has the same orientation as cell (so that the face numbers coincide)
+                // recursive call
+                update_periodic_face_map_recursively<dim, spacedim>
+                (cell_1->child(child_cell_1), cell_2->child(child_cell_2),
+                 n_face_1, n_face_2, face_orientation, face_flip, face_rotation,
+                 periodic_face_map);
+              }
+          }
+        else //only face_1 has children
+          {
+            for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
+              {
+                // find subcell ids that belong to the subface indices
+                unsigned int child_cell_1
+                  = GeometryInfo<dim>::child_cell_on_face(cell_1->refinement_case(), n_face_1, i, face_orientation,
+                                                          face_flip, face_rotation, face_1->refinement_case());
+
+                // recursive call
+                update_periodic_face_map_recursively<dim, spacedim>
+                (cell_1->child(child_cell_1), cell_2,
+                 n_face_1, n_face_2, face_orientation, face_flip, face_rotation,
+                 periodic_face_map);
+              }
+          }
+
+      }
+    else
+      {
+        // Otherwise at least one of the two faces is active and
+        // we need to do some work and enter the periodic face pairs!
+
+        //periodic boundary can only be on the boundary of the domain
+        Assert(cell_1->at_boundary(n_face_1), ExcInternalError());
+        Assert(cell_2->at_boundary(n_face_2), ExcInternalError());
+
+        //Ignore artificial cells, cell_2 may not be active
+        if (!cell_1->is_artificial())
+          {
+            // if cell_2 is active it must not be artificial
+            if (!cell_2->has_children() && cell_2->is_artificial())
+              {
+                std::cout<<"The artificial neighbor of the cell with center " <<cell_1->center()
+                         <<" is the cell with center "<<cell_2->center()<<std::endl;
+                //Assert(!cell_2->is_artificial(), ExcInternalError());
+              }
+            if (cell_2->has_children() || !cell_2->is_artificial())
+              {
+                // insert periodic face pair for both cells
+                typedef std::pair<typename Triangulation<dim,spacedim>::cell_iterator, unsigned int> CellFace;
+                const CellFace cell_face_1 (cell_1, n_face_1);
+                const CellFace cell_face_2 (cell_2, n_face_2);
+                const std::pair<CellFace, CellFace> periodic_faces (cell_face_1, cell_face_2);
+
+                // Only one periodic neighbor is allowed
+                Assert(periodic_face_map.count(cell_face_1) == 0, ExcInternalError());
+                periodic_face_map.insert(periodic_faces);
+              }
+          }
+      }
+  }
+
 }// end of anonymous namespace
 
 
@@ -8823,6 +8976,8 @@ void Triangulation<dim, spacedim>::clear ()
 {
   clear_despite_subscriptions();
   signals.clear();
+  periodic_face_pairs_level_0.clear();
+  periodic_face_map.clear();
 }
 
 
@@ -11588,14 +11743,36 @@ Triangulation<dim,spacedim>::get_triangulation () const
 
 
 template <int dim, int spacedim>
-virtual void
-Triangulation<dim, spacedim>::add_periodicity()
+void
+Triangulation<dim, spacedim>::add_periodicity
+(const std::vector<GridTools::PeriodicFacePair<cell_iterator> > &
+ periodicity_vector)
 {
+  typedef std::pair<cell_iterator, unsigned int> CellFace;
   periodic_face_pairs_level_0.insert(periodic_face_pairs_level_0.end(),
                                      periodicity_vector.begin(),
                                      periodicity_vector.end());
+
+  //Now initialize periodic_face_map
+  typename std::vector<GridTools::PeriodicFacePair<cell_iterator> >::const_iterator it;
+  for (it=periodicity_vector.begin(); it!=periodicity_vector.end(); ++it)
+    {
+      const CellFace cell_face_1 (it->cell[0], it->face_idx[0]);
+      const CellFace cell_face_2 (it->cell[1], it->face_idx[1]);
+
+      //Only one periodic neighbor is allowed
+      Assert(periodic_face_map.count(cell_face_1)==0, ExcInternalError());
+      Assert(periodic_face_map.count(cell_face_2)==0, ExcInternalError());
+
+      const std::pair<CellFace, CellFace> periodic_faces_1 (cell_face_1, cell_face_2);
+      const std::pair<CellFace, CellFace> periodic_faces_2 (cell_face_2, cell_face_1);
+
+      periodic_face_map.insert(periodic_faces_1);
+      periodic_face_map.insert(periodic_faces_2);
+    }
 }
 
+
 template <int dim, int spacedim>
 void
 Triangulation<dim, spacedim>::execute_coarsening_and_refinement ()
@@ -11634,6 +11811,14 @@ 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. At this point cells might be
+  // artificial which prevents us from storing the correct faces
+  const parallel::Triangulation< dim, spacedim > *distributed_triangulation
+    = dynamic_cast<const parallel::Triangulation< dim, spacedim > *> (this);
+  if (!distributed_triangulation)
+    update_periodic_face_map();
 }
 
 
@@ -11656,6 +11841,55 @@ Triangulation<dim,spacedim>::reset_active_cell_indices ()
 }
 
 
+template <int dim, int spacedim>
+void
+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;
+  for (it=periodic_face_pairs_level_0.begin(); it!=periodic_face_pairs_level_0.end(); ++it)
+    {
+      update_periodic_face_map_recursively<dim, spacedim>
+      (it->cell[0], it->cell[1], it->face_idx[0], it->face_idx[1],
+       it->orientation[0], it->orientation[1], it->orientation[2],
+       periodic_face_map);
+      //for the other way, we need to invert the orientation
+      bool orientation, flip, rotation;
+      orientation = it->orientation[0];
+      rotation = it->orientation[2];
+      if (it->orientation[0] && it->orientation [2])
+        flip = !it->orientation[1];
+      else
+        flip = it->orientation[1];
+      update_periodic_face_map_recursively<dim, spacedim>
+      (it->cell[1], it->cell[0], it->face_idx[1], it->face_idx[0],
+       orientation, flip, rotation,
+       periodic_face_map);
+    }
+
+  //check consistency
+  std::cout << "check consistency" << std::endl;
+  //cell->neighbor(neighbor)->neighbor_child_on_subface(face_no, subface_no)==cell.
+  typename std::map<std::pair<cell_iterator, unsigned int>,
+           std::pair<cell_iterator, unsigned int> >::const_iterator it_test;
+  for (it_test=periodic_face_map.begin(); it_test!=periodic_face_map.end(); ++it_test)
+    {
+      if (!it_test->second.first->has_children())
+        {
+          if (!periodic_face_map[it_test->second].first->has_children())
+            {
+              // if the neighbor is active as well, the same pair
+              // order swapped has to be in the map
+              Assert(periodic_face_map[it_test->second] == it_test->first,
+                     ExcInternalError());
+            }
+        }
+    }
+}
+
+
 
 template<int dim, int spacedim>
 void