From d30b1bf412e9ada5ddae81ecb8ecd773898e77b6 Mon Sep 17 00:00:00 2001
From: "daniel.arndt" <daniel.arndt@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Fri, 6 Sep 2013 09:05:25 +0000
Subject: [PATCH] identify_periodic_face_pairs creates periodicity information
 for two given boundary_ids on distributed meshes

git-svn-id: https://svn.dealii.org/trunk@30625 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/include/deal.II/grid/grid_tools.h | 26 ++++++++
 deal.II/source/grid/grid_tools.cc         | 81 +++++++++++++++++++++++
 deal.II/source/grid/grid_tools.inst.in    | 32 +++++++++
 3 files changed, 139 insertions(+)

diff --git a/deal.II/include/deal.II/grid/grid_tools.h b/deal.II/include/deal.II/grid/grid_tools.h
index eabd288760..83e40cf654 100644
--- a/deal.II/include/deal.II/grid/grid_tools.h
+++ b/deal.II/include/deal.II/grid/grid_tools.h
@@ -24,6 +24,7 @@
 #include <deal.II/grid/tria_accessor.h>
 #include <deal.II/grid/tria_iterator.h>
 #include <deal.II/fe/mapping_q1.h>
+#include <deal.II/base/std_cxx1x/tuple.h>
 
 #include <bitset>
 #include <list>
@@ -1153,6 +1154,31 @@ namespace GridTools
                                const int                direction,
                                const dealii::Tensor<1,DH::space_dimension> &offset);
 
+  /**
+   * Add periodicity information to the @p periodicity_vector for the
+   * boundaries with boundary id @p b_id1 and  @p b_id2 in cartesian
+   * direction @p direction.
+   *
+   * This function tries to match all faces belonging to the first
+   * boundary with faces belonging to the second boundary
+   * by comparing the center of the cell faces. To find the correct
+   * corresponding faces, the direction argument indicates in which
+   * cartesian direction periodicity should be set.
+   * ((0,1,2) -> (x,y,z)-direction)
+   *
+   * The output of this function can be used in
+   * parallel::distributed::Triangulation::add_periodicity
+   */
+  template<typename DH>
+  void
+  identify_periodic_face_pairs
+    (const DH &dof_handler,
+     const types::boundary_id b_id1,
+     const types::boundary_id b_id2,
+     const unsigned int direction,
+     std::vector<std_cxx1x::tuple<typename DH::cell_iterator, unsigned int,
+                                  typename DH::cell_iterator, unsigned int> >
+       &periodicity_vector);
 
 
   /**
diff --git a/deal.II/source/grid/grid_tools.cc b/deal.II/source/grid/grid_tools.cc
index a49ede55e2..4dfd9ec119 100644
--- a/deal.II/source/grid/grid_tools.cc
+++ b/deal.II/source/grid/grid_tools.cc
@@ -17,6 +17,7 @@
 #include <deal.II/base/std_cxx1x/array.h>
 #include <deal.II/base/quadrature_lib.h>
 #include <deal.II/grid/tria.h>
+#include <deal.II/distributed/tria.h>
 #include <deal.II/grid/tria_accessor.h>
 #include <deal.II/grid/tria_iterator.h>
 #include <deal.II/grid/tria_boundary.h>
@@ -2540,6 +2541,86 @@ next_cell:
   }
 
 
+  template<typename DH>
+  void
+  identify_periodic_face_pairs
+    (const DH &dof_handler,
+     const types::boundary_id b_id1,
+     const types::boundary_id b_id2,
+     const unsigned int direction,
+     std::vector<std_cxx1x::tuple<typename DH::cell_iterator, unsigned int,
+                                  typename DH::cell_iterator, unsigned int> >
+       &periodicity_vector)
+  {
+    static const unsigned int dim = DH::dimension;
+    static const unsigned int spacedim = DH::space_dimension;
+
+    Assert (0<=direction && direction<spacedim,
+            ExcIndexRange(direction, 0, spacedim));
+
+    Assert (b_id1 != b_id2,
+            ExcMessage ("The boundary indicators b_id1 and b_id2 must be"
+                        "different to denote different boundaries."));
+
+    typename std::map<std::pair<typename DH::cell_iterator,
+                                unsigned int>,
+                      Point<spacedim> > face_locations;
+
+    // Collect faces with boundary_indicator b_id1
+    typename DH::cell_iterator cell = dof_handler.begin();
+    typename DH::cell_iterator endc = dof_handler.end();
+    for (; cell != endc; ++cell)
+      for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+        if(cell->face(f)->boundary_indicator() == b_id1)
+        {
+          Point<spacedim> face_center (cell->face(f)->center());
+          face_center(direction)=0.;
+          const std::pair<typename DH::cell_iterator, unsigned int>
+          cell_face_pair = std::make_pair(cell, f);
+          face_locations[cell_face_pair]=face_center;
+        }
+
+      // Match faces with boundary_indicator b_id2 to the ones in
+      //face_locations
+      cell = dof_handler.begin();
+      for (; cell != endc; ++cell)
+        for(unsigned int f=0;f<GeometryInfo<dim>::faces_per_cell;++f)
+          if(cell->face(f)->boundary_indicator() == b_id2)
+          {
+            typename std::map<std::pair<typename DH::cell_iterator,
+                                        unsigned int>,
+                              Point<spacedim> >::const_iterator p
+              = face_locations.begin();
+
+            Point<spacedim> center2 (cell->face(f)->center());
+            center2(direction)=0.;
+
+            for (; p != face_locations.end(); ++p)
+            {
+              if (center2.distance(p->second) < 1e-4*cell->face(f)->diameter())
+              {
+                const std_cxx1x::tuple<typename DH::cell_iterator, unsigned int,
+                                       typename DH::cell_iterator, unsigned int>
+                  periodic_tuple
+                    = std::make_tuple(p->first.first,
+                                      p->first.second,
+                                      cell,
+                                      f);
+
+                periodicity_vector.push_back(periodic_tuple);
+
+                face_locations.erase(p);
+                break;
+              }
+              Assert (p != face_locations.end(),
+                      ExcMessage ("No corresponding face was found!"));
+            }
+          }
+
+          Assert (face_locations.size() == 0,
+                  ExcMessage ("There are unmatched faces!"));
+  }
+
 
 } /* namespace GridTools */
 
diff --git a/deal.II/source/grid/grid_tools.inst.in b/deal.II/source/grid/grid_tools.inst.in
index 4411897224..9b3bc3b7f6 100644
--- a/deal.II/source/grid/grid_tools.inst.in
+++ b/deal.II/source/grid/grid_tools.inst.in
@@ -271,9 +271,41 @@ for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II
                                   int,
                                   const Tensor<1,deal_II_space_dimension> &);
 
+      template
+      void
+      identify_periodic_face_pairs
+        (const X &,
+         const types::boundary_id,
+         const types::boundary_id,
+         const unsigned int,
+         std::vector<std_cxx1x::tuple<typename X::cell_iterator, unsigned int,
+                                      typename X::cell_iterator, unsigned int> > &);
+
     #endif
 
   \}
 #endif
 }
 
+for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS)
+{
+#if deal_II_dimension <= deal_II_space_dimension
+   #if deal_II_dimension >= 2
+     namespace GridTools \{
+       template
+       void
+       identify_periodic_face_pairs
+         (const parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension> &,
+          const types::boundary_id,
+          const types::boundary_id,
+          const unsigned int,
+          std::vector<std_cxx1x::tuple<typename parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension>::cell_iterator,
+                                       unsigned int,
+                                       typename parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension>::cell_iterator, 
+                                       unsigned int> > &);
+     \}
+   #endif
+#endif
+}
+
+
-- 
2.39.5