Store orientation as well

diff --git a/include/deal.II/grid/tria.h b/include/deal.II/grid/tria.h
index 83665a8c2cafc10fa76808b97aa7d41ffd41e5d3..a20531642d2a224a07c97f4559ab11a388e73c9e 100644 (file)
--- a/include/deal.II/grid/tria.h
+++ b/include/deal.II/grid/tria.h
@@ -38,6 +38,7 @@ DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
 #include <list>
 #include <map>
 #include <numeric>
+#include <bitset>
 
 
 DEAL_II_NAMESPACE_OPEN
@@ -1481,7 +1482,7 @@ public:
   /**
    * If add_periodicity() is called, this variable stores the active periodic face pairs.
    */
-  std::map<std::pair<cell_iterator, unsigned int>, std::pair<cell_iterator, unsigned int> > periodic_face_map;
+  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.

diff --git a/source/grid/tria.cc b/source/grid/tria.cc
index c7760c1c8f5fb1e4002abd4852087ebf354e6246..4e06ab94fe5bc22ca069a3cfdc73f882bfd143f4 100644 (file)
--- a/source/grid/tria.cc
+++ b/source/grid/tria.cc
@@ -914,14 +914,18 @@ namespace
   (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,
+   const std::bitset<3> &orientation,
    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)
+   std::pair<std::pair<typename Triangulation<dim, spacedim>::cell_iterator, unsigned int>, std::bitset<3> > > &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);
 
+    const bool face_orientation = orientation[0];
+    const bool face_flip = orientation[1];
+    const bool face_rotation = orientation[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"));
 
@@ -938,7 +942,9 @@ namespace
     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);
+    const std::pair<CellFace, std::bitset<3> > cell_face_orientation_2 (cell_face_2, orientation);
+
+    const std::pair<CellFace, std::pair<CellFace, std::bitset<3> > > periodic_faces (cell_face_1, cell_face_orientation_2);
 
     // Only one periodic neighbor is allowed
     Assert(periodic_face_map.count(cell_face_1) == 0, ExcInternalError());
@@ -1003,18 +1009,22 @@ namespace
 
                 // 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());
+                  = GeometryInfo<dim>::child_cell_on_face
+                    (cell_1->refinement_case(), n_face_1, i, cell_1->face_orientation(n_face_1),
+                     cell_1->face_flip(n_face_1), cell_1->face_rotation(n_face_1), 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());
+                  = GeometryInfo<dim>::child_cell_on_face
+                    (cell_2->refinement_case(), n_face_2, j, cell_2->face_orientation(n_face_2),
+                     cell_2->face_flip(n_face_2), cell_2->face_rotation(n_face_2), face_2->refinement_case());
+
+                Assert(cell_1->child(child_cell_1)->face(n_face_1) == face_1->child(i), ExcInternalError());
+                Assert(cell_2->child(child_cell_2)->face(n_face_2) == face_2->child(j), ExcInternalError());
 
                 // 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);
+                 n_face_1, n_face_2, orientation, periodic_face_map);
               }
           }
         else //only face_1 has children
@@ -1023,14 +1033,13 @@ namespace
               {
                 // 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());
+                  = GeometryInfo<dim>::child_cell_on_face(cell_1->refinement_case(), n_face_1, i, cell_1->face_orientation(n_face_1),
+                                                          cell_1->face_flip(n_face_1), cell_1->face_rotation(n_face_1), 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);
+                 n_face_1, n_face_2, orientation, periodic_face_map);
               }
           }
 
@@ -11732,22 +11741,7 @@ Triangulation<dim, spacedim>::add_periodicity
                                      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);
-    }
+  update_periodic_face_map();
 }
 
 
@@ -11830,38 +11824,37 @@ Triangulation<dim,spacedim>::update_periodic_face_map ()
     {
       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);
+       it->orientation, 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];
+      std::bitset<3> inverted_orientation;
+      {
+        bool orientation, flip, rotation;
+        orientation = it->orientation[0];
+        rotation = it->orientation[2];
+        flip = orientation ? rotation ^ it->orientation[1] : it->orientation[1];
+        inverted_orientation[0] = orientation;
+        inverted_orientation[1] = flip;
+        inverted_orientation[2] = rotation;
+      }
       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);
+       inverted_orientation, 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;
+           std::pair<std::pair<cell_iterator, unsigned int>, std::bitset<3> >  >::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())
+      const Triangulation<dim, spacedim>::cell_iterator cell_1 = it_test->first.first;
+      const Triangulation<dim, spacedim>::cell_iterator cell_2 = it_test->second.first.first;
+      if (cell_1->level() == cell_2->level())
         {
-          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());
-            }
+          // if both cells have the same neighbor, then the same pair
+          // order swapped has to be in the map
+          Assert (periodic_face_map[it_test->second.first].first == it_test->first,
+                  ExcInternalError());
         }
     }
 }