]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Changed the data type for periodic face pairs
authordaniel.arndt <daniel.arndt@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 16 Sep 2013 17:02:05 +0000 (17:02 +0000)
committerdaniel.arndt <daniel.arndt@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 16 Sep 2013 17:02:05 +0000 (17:02 +0000)
git-svn-id: https://svn.dealii.org/trunk@30732 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/distributed/tria.h
deal.II/include/deal.II/dofs/dof_tools.h
deal.II/include/deal.II/grid/grid_tools.h
deal.II/source/distributed/tria.cc
deal.II/source/dofs/dof_tools_constraints.cc
deal.II/source/dofs/dof_tools_constraints.inst.in
deal.II/source/grid/grid_tools.cc
deal.II/source/grid/grid_tools.inst.in

index 95c34afc1d6f467d17712411155702cbb7e629f0..5d647258f1003d70f162268c4994f98fdaecbfa6 100644 (file)
@@ -158,7 +158,8 @@ namespace internal
   }
 }
 
-
+//forward declaration of the data type for periodic face pairs
+namespace GridTools {template <typename CellIterator> struct PeriodicFacePair;}
 
 namespace parallel
 {
@@ -700,15 +701,36 @@ namespace parallel
       get_p4est_tree_to_coarse_cell_permutation() const;
 
 
+
       /**
        * Join faces in the p4est forest due to periodic boundary conditions.
        *
+       * The vector can be filled by the function
+       * GridTools::collect_periodic_faces.
+       *
+       * @todo At the moment just default orientation is implemented.
+       *
+       * @note Before this function can be used the triangulation has to be
+       * initialized and must not be refined.
+       * Calling this function more than once is possible, but not recommended:
+       * The function destroys and rebuilds the p4est forest each time it is called.
+       */
+      void
+      add_periodicity
+      (const std::vector<GridTools::PeriodicFacePair<cell_iterator> >&);
+
+      /**
+       * Same as the function above, but takes a different argument.
+       *
        * The entries in the std::vector should have the form
        * std_cxx1x::tuple<cell1, face_no1, cell2, face_no2>.
        *
        * The vector can be filled by the function
-       * DoFTools::identify_periodic_face_pairs.
-       *
+       * GridTools::identify_periodic_face_pairs.
+       * 
+       * @note This function can only be used if the faces are in
+       * default orientation.
+       * 
        * @note Before this function can be used the triangulation has to be
        * initialized and must not be refined.
        * Calling this function more than once is possible, but not recommended:
@@ -720,6 +742,7 @@ namespace parallel
                                             cell_iterator, unsigned int> >&);
 
 
+
     private:
       /**
        * MPI communicator to be
index e4d237faf3f11a4eee83eaa1f9e2b495c23268d8..9b391971049b0a8f283b4c85802e73f4060a83b2 100644 (file)
@@ -53,6 +53,7 @@ class ConstraintMatrix;
 template <class GridClass> class InterGridMap;
 template <int dim, int spacedim> class Mapping;
 
+namespace GridTools {template <typename CellIterator> struct PeriodicFacePair;}
 
 //TODO: map_support_points_to_dofs should generate a multimap, rather than just a map, since several dofs may be located at the same support point
 
@@ -1087,13 +1088,15 @@ namespace DoFTools
    */
   template<typename FaceIterator>
   void
-  make_periodicity_constraints (const FaceIterator                          &face_1,
-                                const typename identity<FaceIterator>::type &face_2,
-                                dealii::ConstraintMatrix                    &constraint_matrix,
-                                const ComponentMask                         &component_mask = ComponentMask(),
-                                const bool face_orientation = true,
-                                const bool face_flip = false,
-                                const bool face_rotation = false);
+  make_periodicity_constraints
+  (const FaceIterator                          &face_1,
+   const typename identity<FaceIterator>::type &face_2,
+   dealii::ConstraintMatrix                    &constraint_matrix,
+   const ComponentMask                         &component_mask = ComponentMask(),
+   const bool                                  face_orientation = true,
+   const bool                                  face_flip = false,
+   const bool                                  face_rotation = false);
+
 
 
   /**
@@ -1143,12 +1146,13 @@ namespace DoFTools
    */
   template<typename DH>
   void
-  make_periodicity_constraints (const DH                 &dof_handler,
-                                const types::boundary_id b_id1,
-                                const types::boundary_id b_id2,
-                                const int                direction,
-                                dealii::ConstraintMatrix &constraint_matrix,
-                                const ComponentMask      &component_mask = ComponentMask());
+  make_periodicity_constraints
+  (const DH                 &dof_handler,
+   const types::boundary_id b_id1,
+   const types::boundary_id b_id2,
+   const int                direction,
+   dealii::ConstraintMatrix &constraint_matrix,
+   const ComponentMask      &component_mask = ComponentMask());
 
 
   /**
@@ -1167,17 +1171,18 @@ namespace DoFTools
    *
    * @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
    *
-   * @author Matthias Maier, 2012
+   * @author Daniel Arndt, 2013, Matthias Maier, 2012
    */
   template<typename DH>
   void
-  make_periodicity_constraints (const DH                 &dof_handler,
-                                const types::boundary_id b_id1,
-                                const types::boundary_id b_id2,
-                                const int                direction,
-                                dealii::Tensor<1,DH::space_dimension> &offset,
-                                dealii::ConstraintMatrix &constraint_matrix,
-                                const ComponentMask      &component_mask = ComponentMask());
+  make_periodicity_constraints
+  (const DH                              &dof_handler,
+   const types::boundary_id              b_id1,
+   const types::boundary_id              b_id2,
+   const int                             direction,
+   dealii::Tensor<1,DH::space_dimension> &offset,
+   dealii::ConstraintMatrix              &constraint_matrix,
+   const ComponentMask                   &component_mask = ComponentMask());
 
 
   /**
@@ -1204,11 +1209,12 @@ namespace DoFTools
    */
   template<typename DH>
   void
-  make_periodicity_constraints (const DH                 &dof_handler,
-                                const types::boundary_id b_id,
-                                const int                direction,
-                                dealii::ConstraintMatrix &constraint_matrix,
-                                const ComponentMask      &component_mask = ComponentMask());
+  make_periodicity_constraints
+  (const DH                 &dof_handler,
+   const types::boundary_id b_id,
+   const int                direction,
+   dealii::ConstraintMatrix &constraint_matrix,
+   const ComponentMask      &component_mask = ComponentMask());
 
 
   /**
@@ -1233,12 +1239,31 @@ namespace DoFTools
    */
   template<typename DH>
   void
-  make_periodicity_constraints (const DH                 &dof_handler,
-                                const types::boundary_id b_id,
-                                const int                direction,
-                                dealii::Tensor<1,DH::space_dimension> &offset,
-                                dealii::ConstraintMatrix &constraint_matrix,
-                                const ComponentMask      &component_mask = ComponentMask());
+  make_periodicity_constraints
+  (const DH                              &dof_handler,
+   const types::boundary_id              b_id,
+   const int                             direction,
+   dealii::Tensor<1,DH::space_dimension> &offset,
+   dealii::ConstraintMatrix              &constraint_matrix,
+   const ComponentMask                   &component_mask = ComponentMask());
+
+  /**
+   * Same as above but the periodicity information is given by
+   * @p periodic_faces. This std::vector can be created by
+   * GridTools::collect_periodic_faces.
+   *
+   * @note For DoFHandler objects that are built on a
+   * parallel::distributed::Triangulation object
+   * parallel::distributed::Triangulation::add_periodicity has to be called
+   * before.
+   */
+  template<typename DH>
+  void
+  make_periodicity_constraints
+  (const std::vector<GridTools::PeriodicFacePair<typename DH::cell_iterator> > 
+     &periodic_faces,
+   dealii::ConstraintMatrix &constraint_matrix,
+   const ComponentMask      &component_mask = ComponentMask());
 
 
   //@}
index 83e40cf654f02171e39c098de3afc821857c08ff..fc43a3cd464ba500789d920a7ce3fe0c502991f8 100644 (file)
@@ -990,7 +990,18 @@ namespace GridTools
                              const std::set<types::boundary_id> &boundary_ids
                              = std::set<types::boundary_id>());
 
-
+  /**
+   * Data type that provides all the information that is needed
+   * to create periodicity constraints and a periodic p4est forest
+   * with respect to two periodic cell faces
+   */
+  template<typename CellIterator>
+  struct PeriodicFacePair
+  {
+    CellIterator cell[2];
+    unsigned int face_idx[2];
+    std::bitset<3> orientation;
+  };
 
   /**
    * An orthogonal equality test for faces.
@@ -1076,8 +1087,8 @@ namespace GridTools
 
 
   /**
-   * This function loops over all faces from @p begin to @p end and
-   * collects a set of periodic face pairs for @p direction:
+   * This function will collect periodic face pairs on the highest (i.e.
+   * coarsest) mesh level.
    *
    * Define a 'first' boundary as all boundary faces having boundary_id
    * @p b_id1 and a 'second' boundary consisting of all faces belonging
@@ -1087,45 +1098,33 @@ namespace GridTools
    * boundary with faces belonging to the second boundary with the help
    * of orthogonal_equality().
    *
-   * The bitset that is returned together with the second face encodes the
+   * The bitset that is returned inside of PeriodicFacePair encodes the
    * _relative_ orientation of the first face with respect to the second
    * face, see the documentation of orthogonal_equality for further details.
    *
+   * The @p direction refers to the space direction in which periodicity
+   * is enforced.
+   *
    * The @p offset is a vector tangential to the faces that is added to the
    * location of vertices of the 'first' boundary when attempting to match
    * them to the corresponding vertices of the 'second' boundary. This can
    * be used to implement conditions such as $u(0,y)=u(1,y+1)$.
    *
-   * @author Matthias Maier, 2012
-   */
-  template<typename FaceIterator>
-  std::map<FaceIterator, std::pair<FaceIterator, std::bitset<3> > >
-  collect_periodic_face_pairs (const FaceIterator                          &begin,
-                               const typename identity<FaceIterator>::type &end,
-                               const types::boundary_id                    b_id1,
-                               const types::boundary_id                    b_id2,
-                               const int                                   direction,
-                               const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset);
-
-
-  /**
-   * Same function as above, but accepts a Triangulation or DoFHandler
-   * object @p container (a container is a collection of objects, here a
-   * collection of cells) instead of an explicit face iterator range.
-   *
-   * This function will collect periodic face pairs on the highest (i.e.
-   * coarsest) mesh level.
+   * @note The created std::vector can be used in
+   * DoFTools::make_periodicity_constraints and in
+   * parallel::distributed::Triangulation::add_periodicity to enforce
+   * periodicity algebraically.
    *
-   * @author Matthias Maier, 2012
+   * @author Daniel Arndt, 2013
    */
   template<typename DH>
-  std::map<typename DH::face_iterator, std::pair<typename DH::face_iterator, std::bitset<3> > >
-  collect_periodic_face_pairs (const DH                 &container,
-                               const types::boundary_id b_id1,
-                               const types::boundary_id b_id2,
-                               const int                direction,
-                               const dealii::Tensor<1,DH::space_dimension> &offset);
-
+  std::vector<PeriodicFacePair<typename DH::cell_iterator> >
+  collect_periodic_faces
+  (const DH                 &dof_handler,
+   const types::boundary_id b_id1,
+   const types::boundary_id b_id2,
+   const unsigned int       direction,
+   const dealii::Tensor<1,DH::space_dimension> &offset);
 
   /**
    * This compatibility version of collect_periodic_face_pairs only works
@@ -1145,40 +1144,101 @@ namespace GridTools
    * meshes with cells not in @ref GlossFaceOrientation
    * "standard orientation".
    *
+   * @author Daniel Arndt, 2013
+   */
+  template<typename DH>
+  std::vector<PeriodicFacePair<typename DH::cell_iterator> >
+  collect_periodic_faces
+  (const DH                 &dof_handler,
+   const types::boundary_id b_id,
+   const unsigned int       direction,
+   const dealii::Tensor<1,DH::space_dimension> &offset);
+
+  /**
+   * This function does the same as collect_periodic_faces but returns a
+   * different data type.
+   *
+   * @author Matthias Maier, 2012
+   *
+   * @note The returned data type is not compatible with
+   * DoFTools::make_periodicity_constraints
+   * 
+   * @deprecated
+   */
+  template<typename DH>
+  std::map<typename DH::face_iterator, std::pair<typename DH::face_iterator, std::bitset<3> > >
+  collect_periodic_face_pairs
+  (const DH                 &container,
+   const types::boundary_id b_id1,
+   const types::boundary_id b_id2,
+   int                      direction,
+   const dealii::Tensor<1,DH::space_dimension> &offset) DEAL_II_DEPRECATED;
+
+  /**
+   * This compatibility version of collect_periodic_face_pairs only works
+   * on grids with cells in @ref GlossFaceOrientation "standard orientation".
+   *
    * @author Matthias Maier, 2012
+   *
+   * @note The returned data type is not compatible with
+   * DoFTools::make_periodicity_constraints
+   *
+   * @deprecated
    */
   template<typename DH>
   std::map<typename DH::face_iterator, typename DH::face_iterator>
-  collect_periodic_face_pairs (const DH                 &dof_handler, /*TODO: Name*/
-                               const types::boundary_id b_id,
-                               const int                direction,
-                               const dealii::Tensor<1,DH::space_dimension> &offset);
+  collect_periodic_face_pairs
+  (const DH                 &dof_handler, /*TODO: Name*/
+   const types::boundary_id b_id,
+   int                      direction,
+   const dealii::Tensor<1,DH::space_dimension> &offset) DEAL_II_DEPRECATED;
 
   /**
-   * 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 version loops over all faces from @p begin to @p end
+   * instead of accepting a DoFHandler or a Triangulation.
    *
-   * 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)
+   * @author Matthias Maier, 2012
+   *
+   * @note This function cannot produce the return as the other
+   * collect_periodic_faces functions.
+   *
+   * @deprecated
+   */
+  template<typename FaceIterator>
+  std::map<FaceIterator, std::pair<FaceIterator, std::bitset<3> > >
+  collect_periodic_face_pairs
+  (const FaceIterator                          &begin,
+   const typename identity<FaceIterator>::type &end,
+   const types::boundary_id                    b_id1,
+   const types::boundary_id                    b_id2,
+   const int                                   direction,
+   const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset)
+  DEAL_II_DEPRECATED;
+
+  /**
+   * This function does the same as collect_periodic_faces but returns a
+   * different data type.
+   *
+   * @author Daniel Arndt, 2013
    *
-   * The output of this function can be used in
+   * @note The returned data type is not compatible with
+   * DoFTools::make_periodicity_constraints, but with a version of
    * parallel::distributed::Triangulation::add_periodicity
+   *
+   * @note Use collect_periodic_faces instead.
+   * 
+   * @deprecated 
    */
   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);
+  (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) DEAL_II_DEPRECATED;
 
 
   /**
index ce4cbc6e975a34da9ccb9eea5cf90e6282a1a77c..b036739d6f70776a417b7bb43d69c3da80b788c3 100644 (file)
@@ -23,8 +23,8 @@
 #include <deal.II/grid/tria.h>
 #include <deal.II/grid/tria_accessor.h>
 #include <deal.II/grid/tria_iterator.h>
-#include <deal.II/distributed/tria.h>
 #include <deal.II/grid/grid_tools.h>
+#include <deal.II/distributed/tria.h>
 
 #ifdef DEAL_II_WITH_P4EST
 #  include <p4est_bits.h>
@@ -3436,39 +3436,30 @@ namespace parallel
       return mpi_communicator;
     }
 
-
-    template <int dim, int spacedim>
+    template<int dim, int spacedim>
     void
     Triangulation<dim,spacedim>::add_periodicity
-       (const std::vector<std_cxx1x::tuple<cell_iterator, unsigned int,
-                                           cell_iterator, unsigned int> >&
-          periodicity_vector)
+      (const std::vector<GridTools::PeriodicFacePair<cell_iterator> >&
+        periodicity_vector)
     {
-#if DEAL_II_P4EST_VERSION_GTE(0,3,4,1)
+      #if DEAL_II_P4EST_VERSION_GTE(0,3,4,1)
       Assert (triangulation_has_content == true,
               ExcMessage ("The triangulation is empty!"));
       Assert (this->n_levels() == 1,
               ExcMessage ("The triangulation is refined!"));
-
-      typename std::vector
-        <typename std_cxx1x::tuple
-          <cell_iterator, unsigned int,
-           cell_iterator, unsigned int> >::const_iterator periodic_tuple;
-      periodic_tuple = periodicity_vector.begin();
-
-      typename std::vector
-        <typename std_cxx1x::tuple
-          <cell_iterator, unsigned int,
-           cell_iterator, unsigned int> >::const_iterator periodic_end;
-      periodic_end = periodicity_vector.end();
-
-      for (; periodic_tuple<periodic_end; ++periodic_tuple)
+      
+      typedef std::vector<GridTools::PeriodicFacePair<cell_iterator> >
+        FaceVector;
+      typename FaceVector::const_iterator it, periodic_end;
+      it = periodicity_vector.begin();
+      
+      for (; it<periodic_end; ++it)
       {
-        const cell_iterator first_cell=std_cxx1x::get<0>(*periodic_tuple);
-        const cell_iterator second_cell=std_cxx1x::get<2>(*periodic_tuple);
-        const unsigned int face_right=std_cxx1x::get<3>(*periodic_tuple);
-        const unsigned int face_left=std_cxx1x::get<1>(*periodic_tuple);
-
+        const cell_iterator first_cell = it->cell[0];
+        const cell_iterator second_cell = it->cell[1];
+        const unsigned int face_left = it->face_idx[0];
+        const unsigned int face_right = it->face_idx[1];
+        
         //respective cells of the matching faces in p4est
         const unsigned int tree_left
           = coarse_cell_to_p4est_tree_permutation[std::distance(this->begin(),
@@ -3477,6 +3468,12 @@ namespace parallel
           = coarse_cell_to_p4est_tree_permutation[std::distance(this->begin(),
                                                                 second_cell)];
 
+        //TODO Add support for non default orientation.
+        Assert(it->orientation == 1,
+               ExcMessage("Found a face match with non standard orientation. "
+                          "This function is only suitable for meshes with "
+                          "cells in default orientation"));
+
         dealii::internal::p4est::functions<dim>::
           connectivity_join_faces (connectivity,
                                    tree_left,
@@ -3484,6 +3481,7 @@ namespace parallel
                                    face_left,
                                    face_right,
                                    /* orientation */ 0);
+
         /* The orientation parameter above describes the difference between
          * the cell on the left and the cell on the right would number of the
          * corners of the face.  In the periodic domains most users will want,
@@ -3495,25 +3493,60 @@ namespace parallel
          * date.
          */
       }
-
-
+      
+      
       Assert(dealii::internal::p4est::functions<dim>::connectivity_is_valid
-               (connectivity) == 1,
-             ExcInternalError());
-
-        // now create a forest out of the connectivity data structure
+               (connectivity) == 1, ExcInternalError());
+      
+      // now create a forest out of the connectivity data structure
       dealii::internal::p4est::functions<dim>::destroy (parallel_forest);
       parallel_forest
         = dealii::internal::p4est::functions<dim>::
             new_forest (mpi_communicator,
                         connectivity,
-                        /* minimum initial number of quadrants per tree */ 0,
-                        /* minimum level of upfront refinement */ 0,
-                        /* use uniform upfront refinement */ 1,
-                        /* user_data_size = */ 0,
-                        /* user_data_constructor = */ NULL,
-                        /* user_pointer */ this);
+                       /* minimum initial number of quadrants per tree */ 0,
+                       /* minimum level of upfront refinement */ 0,
+                       /* use uniform upfront refinement */ 1,
+                       /* user_data_size = */ 0,
+                       /* user_data_constructor = */ NULL,
+                       /* user_pointer */ this);
+      
+      #else
+      Assert(false, ExcMessage ("Need p4est version >= 0.3.4.1!"));
+      #endif
+    }
+
+    
+    template <int dim, int spacedim>
+    void
+    Triangulation<dim,spacedim>::add_periodicity
+       (const std::vector<std_cxx1x::tuple<cell_iterator, unsigned int,
+                                           cell_iterator, unsigned int> >&
+          periodicity_vector)
+    {
+#if DEAL_II_P4EST_VERSION_GTE(0,3,4,1)
+      typedef std::vector<std_cxx1x::tuple <cell_iterator, unsigned int,
+                                            cell_iterator, unsigned int> >
+        FaceVector;
+
+      typename FaceVector::const_iterator it, end_periodic;
+      it = periodicity_vector.begin();
+      end_periodic = periodicity_vector.end();
+
+      std::vector<GridTools::PeriodicFacePair<cell_iterator> > periodic_faces;
+
+      for(; it!=end_periodic; ++it)
+      {
+        const cell_iterator cell1 = std::get<0> (*it);
+        const cell_iterator cell2 = std::get<2> (*it);
+        const unsigned int face_idx1 = std::get<1> (*it);
+        const unsigned int face_idx2 = std::get<3> (*it);
+        const GridTools::PeriodicFacePair<cell_iterator> matched_face
+          = {{cell1, cell2},{face_idx1, face_idx2}, 1};
+        periodic_faces.push_back(matched_face);
+      }
 
+      add_periodicity(periodic_faces);
 #else
       Assert(false, ExcMessage ("Need p4est version >= 0.3.4.1!"));
 #endif
index 402557bee3c6c35d079709870a4f4945c263c834..adaf93b7f0084543316153f47dd9b5b4aa133e8f 100644 (file)
@@ -1970,7 +1970,7 @@ namespace DoFTools
   }
 
 
-
+  
   template<typename DH>
   void
   make_periodicity_constraints (const DH                       &dof_handler,
@@ -2010,24 +2010,23 @@ namespace DoFTools
             ExcMessage ("The boundary indicators b_id1 and b_id2 must be"
                         "different to denote different boundaries."));
 
-    typedef typename DH::face_iterator FaceIterator;
-    typedef std::map<FaceIterator, std::pair<FaceIterator, std::bitset<3> > > FaceMap;
+    typedef std::vector<GridTools::PeriodicFacePair
+                          <typename DH::cell_iterator> > FaceVector;
 
     // Collect matching periodic cells on the coarsest level:
-    FaceMap matched_cells =
-      GridTools::collect_periodic_face_pairs(dof_handler,
-                                             b_id1, b_id2,
-                                             direction, offset);
+    FaceVector matched_cells =
+      GridTools::collect_periodic_faces(dof_handler, b_id1, b_id2,
+                                        direction, offset);
 
     // And apply the low level make_periodicity_constraints function to
     // every matching pair:
-    for (typename FaceMap::iterator it = matched_cells.begin();
+    for (typename FaceVector::iterator it = matched_cells.begin();
          it != matched_cells.end(); ++it)
       {
         typedef typename DH::face_iterator FaceIterator;
-        const FaceIterator &face_1 = it->first;
-        const FaceIterator &face_2 = it->second.first;
-        const std::bitset<3> &orientation = it->second.second;
+        const FaceIterator &face_1 = it->cell[0]->face(it->face_idx[0]);
+        const FaceIterator &face_2 = it->cell[1]->face(it->face_idx[1]);
+        const std::bitset<3> &orientation = it->orientation;
 
         Assert(face_1->at_boundary() && face_2->at_boundary(),
                ExcInternalError());
@@ -2088,23 +2087,22 @@ namespace DoFTools
     Assert(dim == space_dim,
            ExcNotImplemented());
 
-    typedef typename DH::face_iterator FaceIterator;
-    typedef std::map<FaceIterator, FaceIterator> FaceMap;
+    typedef std::vector<GridTools::PeriodicFacePair
+                          <typename DH::cell_iterator> > FaceVector;
 
     // Collect matching periodic cells on the coarsest level:
-    FaceMap matched_cells =
-      GridTools::collect_periodic_face_pairs(dof_handler,
-                                             b_id,
-                                             direction, offset);
+    const FaceVector matched_cells =
+      GridTools::collect_periodic_faces(dof_handler, b_id,
+                                        direction, offset);
 
     // And apply the low level make_periodicity_constraints function to
     // every matching pair:
-    for (typename FaceMap::iterator it = matched_cells.begin();
+    for (typename FaceVector::const_iterator it = matched_cells.begin();
          it != matched_cells.end(); ++it)
       {
         typedef typename DH::face_iterator FaceIterator;
-        const FaceIterator &face_1 = it->first;
-        const FaceIterator &face_2 = it->second;
+        const FaceIterator &face_1 = it->cell[0]->face(it->face_idx[0]);
+        const FaceIterator &face_2 = it->cell[1]->face(it->face_idx[1]);
 
         Assert(face_1->at_boundary() && face_2->at_boundary(),
                ExcInternalError());
@@ -2126,6 +2124,38 @@ namespace DoFTools
 
 
 
+  template<typename DH>
+  void
+  make_periodicity_constraints
+  (const std::vector<GridTools::PeriodicFacePair<typename DH::cell_iterator> >
+    &periodic_faces,
+   dealii::ConstraintMatrix &constraint_matrix,
+   const ComponentMask      &component_mask)
+  {
+    typedef std::vector<GridTools::PeriodicFacePair<typename DH::cell_iterator> >
+    FaceVector;
+    typename FaceVector::const_iterator it, end_periodic;
+    it = periodic_faces.begin();
+    end_periodic = periodic_faces.end();
+    
+    for(; it!=end_periodic; ++it)
+    {
+      typedef typename DH::face_iterator FaceIterator;
+      const FaceIterator face_1 = it->cell[0]->face(it->face_idx[0]);
+      const FaceIterator face_2 = it->cell[1]->face(it->face_idx[1]);
+      
+      make_periodicity_constraints(face_1,
+                                   face_2,
+                                   constraint_matrix,
+                                   component_mask,
+                                   it->orientation[0],
+                                   it->orientation[1],
+                                   it->orientation[2]);
+    }
+  }
+
+
+
   namespace internal
   {
     namespace
index 246f53aa02a43e677781f2242d5d8c7747a8f7dc..5c9bdd7c889c7d335ff8324f30168843920e68bb 100644 (file)
@@ -68,6 +68,14 @@ for (DH : DOFHANDLERS; deal_II_dimension : DIMENSIONS)
                                          dealii::Tensor<1,DH::space_dimension> &,
                                          dealii::ConstraintMatrix &,
                                          const ComponentMask &);
+
+  template
+  void
+  DoFTools::make_periodicity_constraints<DH>
+  (const std::vector<GridTools::PeriodicFacePair<typename DH::cell_iterator> > &,
+   dealii::ConstraintMatrix &,
+   const ComponentMask &);
+  
 #endif
 }
 
index eeeb0dc6adecd9868dcda3cdc473f868f9715a3d..292b113d624b43120faa3507ad2121cfc883a7b1 100644 (file)
@@ -2217,7 +2217,6 @@ next_cell:
   }
 
 
-
   /*
    * Internally used in orthogonal_equality
    *
@@ -2226,7 +2225,7 @@ next_cell:
    *
    * See the comment on the next function as well as the detailed
    * documentation of make_periodicity_constraints and
-   * collect_periodic_face_pairs for details
+   * collect_periodic_faces for details
    */
   template<int dim> struct OrientationLookupTable {};
 
@@ -2355,14 +2354,82 @@ next_cell:
 
 
   /*
-   * Internally used in collect_periodic_face_pairs
+   * Internally used in collect_periodic_faces
+   */
+  template<typename CellIterator>
+  std::vector<PeriodicFacePair<CellIterator> >
+  match_periodic_face_pairs 
+  (std::set<std::pair<CellIterator, unsigned int> >
+       &pairs1,
+   std::set<std::pair<typename identity<CellIterator>::type, unsigned int> >
+       &pairs2,
+   int direction,
+   const dealii::Tensor<1,CellIterator::AccessorType::space_dimension>
+       &offset)
+  {
+    static const int space_dim = CellIterator::AccessorType::space_dimension;
+    Assert (0<=direction && direction<space_dim,
+            ExcIndexRange (direction, 0, space_dim));
+
+    Assert (pairs1.size() == pairs2.size(),
+            ExcMessage ("Unmatched faces on periodic boundaries"));
+
+    typename std::vector<PeriodicFacePair<CellIterator> > matched_faces;
+
+    // Match with a complexity of O(n^2). This could be improved...
+    std::bitset<3> orientation;
+    typedef typename std::set
+      <std::pair<CellIterator, unsigned int> >::const_iterator PairIterator;
+    for (PairIterator it1 = pairs1.begin(); it1 != pairs1.end(); ++it1)
+      {
+        for (PairIterator it2 = pairs2.begin(); it2 != pairs2.end(); ++it2)
+          {
+            const CellIterator cell1 = it1->first;
+            const CellIterator cell2 = it2->first;
+            const unsigned int face_idx1 = it1->second;
+            const unsigned int face_idx2 = it2->second;
+            if (GridTools::orthogonal_equality(orientation,
+                                               cell1->face(face_idx1),
+                                               cell2->face(face_idx2),
+                                               direction, offset))
+              {
+                // We have a match, so insert the matching pairs and
+                // remove the matched cell in pairs2 to speed up the
+                // matching:
+                const PeriodicFacePair<CellIterator> matched_face
+                  = {{cell1, cell2},{face_idx1, face_idx2}, orientation};
+                matched_faces.push_back(matched_face);
+                pairs2.erase(it2);
+                break;
+              }
+          }
+      }
+
+    AssertThrow (matched_faces.size() == pairs1.size() && pairs2.size() == 0,
+                 ExcMessage ("Unmatched faces on periodic boundaries"));
+
+    return matched_faces;
+  }
+
+  /* Deprecated version of the function above with different return value.
+   * It is used the deprecated collect_periodic_face_pairs.
    */
   template<typename FaceIterator>
   std::map<FaceIterator, std::pair<FaceIterator, std::bitset<3> > >
-  match_periodic_face_pairs (std::set<FaceIterator> &faces1, /* not const! */
-                             std::set<typename identity<FaceIterator>::type> &faces2, /* not const! */
-                             int                         direction,
-                             const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset)
+  match_periodic_face_pairs
+  (std::set<FaceIterator>                          &faces1, /* not const! */
+   std::set<typename identity<FaceIterator>::type> &faces2, /* not const! */
+   int                                             direction,
+   const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset)
+  DEAL_II_DEPRECATED;
+
+  template<typename FaceIterator>
+  std::map<FaceIterator, std::pair<FaceIterator, std::bitset<3> > >
+  match_periodic_face_pairs
+  (std::set<FaceIterator> &faces1, /* not const! */
+   std::set<typename identity<FaceIterator>::type> &faces2, /* not const! */
+   int                    direction,
+   const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset)
   {
     static const int space_dim = FaceIterator::AccessorType::space_dimension;
     Assert (0<=direction && direction<space_dim,
@@ -2401,6 +2468,119 @@ next_cell:
   }
 
 
+  template<typename DH>
+  std::vector<PeriodicFacePair<typename DH::cell_iterator> >
+  collect_periodic_faces
+    (const DH                 &dof_handler,
+     const types::boundary_id b_id1,
+     const types::boundary_id b_id2,
+     const unsigned int       direction,
+     const dealii::Tensor<1,DH::space_dimension> &offset)
+  {
+    static const unsigned int dim = DH::dimension;
+    static const unsigned int space_dim = DH::space_dimension;
+    Assert (0<=direction && direction<space_dim,
+            ExcIndexRange (direction, 0, space_dim));
+
+    // Loop over all cells on the highest level and collect all boundary
+    // faces belonging to b_id1 and b_id2:
+
+    std::set<std::pair<typename DH::cell_iterator, unsigned int> > pairs1;
+    std::set<std::pair<typename DH::cell_iterator, unsigned int> > pairs2;
+
+    for (typename DH::cell_iterator cell = dof_handler.begin(0);
+         cell != dof_handler.end(0); ++cell)
+      {
+        for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
+          {
+            const typename DH::face_iterator face = cell->face(i);
+            if (face->at_boundary() && face->boundary_indicator() == b_id1)
+            {
+              const std::pair<typename DH::cell_iterator, unsigned int> pair1
+                = std::make_pair(cell, i);
+              pairs1.insert(pair1);
+            } 
+
+            if (face->at_boundary() && face->boundary_indicator() == b_id2)
+            {
+              const std::pair<typename DH::cell_iterator, unsigned int> pair2
+                = std::make_pair(cell, i);
+              pairs2.insert(pair2);
+            }
+          }
+      }
+
+    Assert (pairs1.size() == pairs2.size(),
+            ExcMessage ("Unmatched faces on periodic boundaries"));
+
+    // and call match_periodic_face_pairs that does the actual matching:
+    return match_periodic_face_pairs(pairs1, pairs2, direction, offset);
+  }
+
+  template<typename DH>
+  std::vector<PeriodicFacePair<typename DH::cell_iterator> >
+  collect_periodic_faces (const DH                 &dof_handler,
+                          const types::boundary_id b_id,
+                          const unsigned int       direction,
+                          const dealii::Tensor<1,DH::space_dimension> &offset)
+  {
+    static const unsigned int dim = DH::dimension;
+    static const unsigned int space_dim = DH::space_dimension;
+    Assert (0<=direction && direction<space_dim,
+            ExcIndexRange (direction, 0, space_dim));
+    
+    Assert(dim == space_dim,
+           ExcNotImplemented());
+    
+    // Loop over all cells on the highest level and collect all boundary
+    // faces 2*direction and 2*direction*1:
+    
+    std::set<std::pair<typename DH::cell_iterator, unsigned int> > pairs1;
+    std::set<std::pair<typename DH::cell_iterator, unsigned int> > pairs2;
+    
+    for (typename DH::cell_iterator cell = dof_handler.begin(0);
+         cell != dof_handler.end(0); ++cell)
+         {
+           const typename DH::face_iterator face_1 = cell->face(2*direction);
+           const typename DH::face_iterator face_2 = cell->face(2*direction+1);
+           
+           if (face_1->at_boundary() && face_1->boundary_indicator() == b_id)
+           {
+             const std::pair<typename DH::cell_iterator, unsigned int> pair1
+             = std::make_pair(cell, 2*direction);
+             pairs1.insert(pair1);
+           }
+           
+           if (face_2->at_boundary() && face_2->boundary_indicator() == b_id)
+           {
+             const std::pair<typename DH::cell_iterator, unsigned int> pair1
+             = std::make_pair(cell, 2*direction+1);
+             pairs1.insert(pair1);
+           }
+         }
+         
+    Assert (pairs1.size() == pairs2.size(),
+            ExcMessage ("Unmatched faces on periodic boundaries"));
+                 
+    // and call match_periodic_face_pairs that does the actual matching:
+
+    typedef std::vector<PeriodicFacePair<typename DH::cell_iterator> >
+      FaceVector;
+
+    FaceVector matching = match_periodic_face_pairs(pairs1, pairs2,
+                                                    direction, offset);
+
+    for (typename FaceVector::iterator pairing = matching.begin();
+         pairing != matching.end(); ++pairing)
+    {
+      Assert(pairing->orientation == 1,
+      ExcMessage("Found a face match with non standard orientation. "
+                 "This function is only suitable for meshes with cells "
+                 "in default orientation"));
+    }
+                      
+    return matching;
+  }
 
   template<typename FaceIterator>
   std::map<FaceIterator, std::pair<FaceIterator, std::bitset<3> > >
@@ -2436,8 +2616,6 @@ next_cell:
     return match_periodic_face_pairs(faces1, faces2, direction, offset);
   }
 
-
-
   template<typename DH>
   std::map<typename DH::face_iterator, std::pair<typename DH::face_iterator, std::bitset<3> > >
   collect_periodic_face_pairs (const DH                 &dof_handler,
@@ -2446,41 +2624,26 @@ next_cell:
                                int                      direction,
                                const dealii::Tensor<1,DH::space_dimension> &offset)
   {
-    static const int dim = DH::dimension;
-    static const int space_dim = DH::space_dimension;
-    Assert (0<=direction && direction<space_dim,
-            ExcIndexRange (direction, 0, space_dim));
-
-    // Loop over all cells on the highest level and collect all boundary
-    // faces belonging to b_id1 and b_id2:
+    typedef std::vector<PeriodicFacePair<typename DH::cell_iterator> > FaceVector;
 
-    std::set<typename DH::face_iterator> faces1;
-    std::set<typename DH::face_iterator> faces2;
+    const FaceVector face_vector
+      = collect_periodic_faces (dof_handler, b_id1, b_id2, direction, offset);
 
-    for (typename DH::cell_iterator cell = dof_handler.begin(0);
-         cell != dof_handler.end(0); ++cell)
-      {
-        for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
-          {
-            const typename DH::face_iterator face = cell->face(i);
-
-            if (face->at_boundary() && face->boundary_indicator() == b_id1)
-              faces1.insert(face);
-
-            if (face->at_boundary() && face->boundary_indicator() == b_id2)
-              faces2.insert(face);
-          }
-      }
-
-    Assert (faces1.size() == faces2.size(),
-            ExcMessage ("Unmatched faces on periodic boundaries"));
+    std::map<typename DH::face_iterator,
+             std::pair<typename DH::face_iterator, std::bitset<3> > >
+      return_value;
+    for(typename FaceVector::const_iterator it = face_vector.begin();
+        it != face_vector.end(); ++it)
+    {
+      const typename DH::face_iterator face1 = it->cell[0]->face(it->face_idx[0]);
+      const typename DH::face_iterator face2 = it->cell[1]->face(it->face_idx[1]);
+      return_value[face1] = std::make_pair(face2, it->orientation);
+    }
 
-    // and call match_periodic_face_pairs that does the actual matching:
-    return match_periodic_face_pairs(faces1, faces2, direction, offset);
+    return return_value;
   }
 
-
-
+  
   template<typename DH>
   std::map<typename DH::face_iterator, typename DH::face_iterator>
   collect_periodic_face_pairs (const DH                 &dof_handler,
@@ -2488,59 +2651,24 @@ next_cell:
                                int                      direction,
                                const dealii::Tensor<1,DH::space_dimension> &offset)
   {
-    static const int dim = DH::dimension;
-    static const int space_dim = DH::space_dimension;
-    Assert (0<=direction && direction<space_dim,
-            ExcIndexRange (direction, 0, space_dim));
-
-    Assert(dim == space_dim,
-           ExcNotImplemented());
-
-    // Loop over all cells on the highest level and collect all boundary
-    // faces 2*direction and 2*direction*1:
-
-    std::set<typename DH::face_iterator> faces1;
-    std::set<typename DH::face_iterator> faces2;
-
-    for (typename DH::cell_iterator cell = dof_handler.begin(0);
-         cell != dof_handler.end(0); ++cell)
-      {
-        const typename DH::face_iterator face_1 = cell->face(2*direction);
-        const typename DH::face_iterator face_2 = cell->face(2*direction+1);
-
-        if (face_1->at_boundary() && face_1->boundary_indicator() == b_id)
-          faces1.insert(face_1);
-
-        if (face_2->at_boundary() && face_2->boundary_indicator() == b_id)
-          faces2.insert(face_2);
-      }
-
-    Assert (faces1.size() == faces2.size(),
-            ExcMessage ("Unmatched faces on periodic boundaries"));
-
-    // and call match_periodic_face_pairs that does the actual matching:
-
-    typedef std::map<typename DH::face_iterator, std::pair<typename DH::face_iterator, std::bitset<3> > > FaceMap;
-    FaceMap matching = match_periodic_face_pairs(faces1, faces2, direction, offset);
-
-    std::map<typename DH::face_iterator, typename DH::face_iterator>
-    return_value;
-
-    for (typename FaceMap::iterator pairing = matching.begin();
-         pairing != matching.end(); ++pairing)
-      {
-        Assert(pairing->second.second == 1,
-               ExcMessage("Found a face match with non standard orientation. "
-                          "This function is only suitable for meshes with cells "
-                          "in default orientation"));
-
-        return_value[pairing->first] = pairing->second.first;
-      }
-
+    typedef std::vector<PeriodicFacePair<typename DH::cell_iterator> > FaceVector;
+    
+    const FaceVector face_vector
+      = collect_periodic_faces (dof_handler, b_id, direction, offset);
+    
+    std::map<typename DH::face_iterator, typename DH::face_iterator> return_value;
+    for(typename FaceVector::const_iterator it = face_vector.begin();
+        it != face_vector.end(); ++it)
+    {
+      const typename DH::face_iterator face1 = it->cell[0]->face(it->face_idx[0]);
+      const typename DH::face_iterator face2 = it->cell[1]->face(it->face_idx[1]);
+      return_value[face1] = face2;
+    }
+    
     return return_value;
   }
 
-
   template<typename DH>
   void
   identify_periodic_face_pairs
@@ -2552,72 +2680,26 @@ next_cell:
                                   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."));
-
-    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> >::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 (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!"));
-            }
-          }
+    typedef std::vector<PeriodicFacePair<typename DH::cell_iterator> >
+      FaceVector;
+    const FaceVector periodic_faces
+      = collect_periodic_faces(dof_handler, b_id1, b_id2, direction,
+                               dealii::Tensor<1,DH::space_dimension> ());
+
+    typename FaceVector::const_iterator it, end_faces;
+    it = periodic_faces.begin();
+    end_faces = periodic_faces.end();
+    
+    for(; it!=end_faces; ++it)
+    {
 
-          Assert (face_locations.size() == 0,
-                  ExcMessage ("There are unmatched faces!"));
+      const std_cxx1x::tuple<typename DH::cell_iterator, unsigned int,
+                             typename DH::cell_iterator, unsigned int>
+        periodic_tuple (it->cell[0], it->face_idx[0],
+                        it->cell[1], it->face_idx[1]);
+      
+      periodicity_vector.push_back(periodic_tuple);
+    }
   }
 
 
index fdaecb93c5f189abf133c1bd1cd79d5d79153aed..465931c953538800ef5f21bcdef15bdc655f4299 100644 (file)
@@ -234,6 +234,23 @@ for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II
 
     #if deal_II_dimension >= 2
 
+      template
+      std::vector<PeriodicFacePair<X::cell_iterator> >
+      collect_periodic_faces
+                                 (const X &,
+                                  const types::boundary_id,
+                                  const types::boundary_id,
+                                  unsigned int,
+                                  const Tensor<1,deal_II_space_dimension> &);
+    
+      template<typename DH>
+      std::vector<PeriodicFacePair<X::cell_iterator> >
+      collect_periodic_faces
+                                 (const X &,
+                                  const types::boundary_id,
+                                  unsigned int,
+                                  const Tensor<1,deal_II_space_dimension> &);
+
       template
       std::map<X::face_iterator, std::pair<X::face_iterator, std::bitset<3> > >
       collect_periodic_face_pairs
@@ -243,7 +260,6 @@ for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II
                                   const types::boundary_id,
                                   int,
                                   const Tensor<1,deal_II_space_dimension> &);
-
       template
       std::map<X::active_face_iterator, std::pair<X::active_face_iterator, std::bitset<3> > >
       collect_periodic_face_pairs
@@ -270,7 +286,7 @@ for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II
                                   const types::boundary_id,
                                   int,
                                   const Tensor<1,deal_II_space_dimension> &);
-
       template
       void
       identify_periodic_face_pairs

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.