]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Bugfix: Make DoFTools::make_periodicity_constraints orientation aware
authorMatthias Maier <tamiko@kyomu.43-1.org>
Sat, 1 Dec 2012 17:17:25 +0000 (17:17 +0000)
committerMatthias Maier <tamiko@kyomu.43-1.org>
Sat, 1 Dec 2012 17:17:25 +0000 (17:17 +0000)
With this it is now possible to prescribe an orientation for the low level
interface, the high level interface will detect the necessary orientation
automatically, *yay*.

git-svn-id: https://svn.dealii.org/trunk@27717 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/dofs/dof_tools.h
deal.II/source/dofs/dof_tools.cc
deal.II/source/dofs/dof_tools.inst.in

index c39b590c478a225d541deef902dc228527c90ca1..5a9b032536ba87bbc4263ecee2007b1695c6c4ab 100644 (file)
@@ -936,197 +936,233 @@ namespace DoFTools
    */
 
   /**
-   * Insert the (algebraic) constraints due
-   * to periodic boundary conditions into
-   * a ConstraintMatrix @p
-   * constraint_matrix.
-   *
-   * Given a pair of not necessarily
-   * active faces @p face_1 and @p face_2,
-   * this functions constrains all DoFs
-   * associated with the boundary
-   * described by @p face_1 to the
-   * respective DoFs of the boundary
-   * described by @p face_2. More
-   * precisely:
-   *
-   * If @p face_1 and @p face_2 are both
-   * active faces it adds the DoFs of
-   * @p face_1 to the list of constrained
-   * DoFs in @p constraint_matrix and adds
-   * lines to constrain them to the
-   * corresponding values of the DoFs on
-   * @p face_2.
-   * Otherwise it loops recursively over
-   * the children of @p face_1 and @p face_2.
-   *
-   * For this to work @p face_1 and @p face_2
-   * must have the same refinement
-   * history, i.e. either @p face_1 and
-   * @p face_2 must be active faces or
-   * must be isotropically refined and
-   * have the same number of child faces
-   * that recursively obey this rule.
-   * (The anisotropic case is not yet
-   * implemented.)
-   *
-   * This routine only constrains DoFs that
-   * are not already constrained.
-   * If this routine encounters a DoF that
-   * already is constrained (for instance
-   * by Dirichlet boundary conditions),
-   * the old setting of the constraint
-   * (dofs the entry is constrained to,
-   * inhomogeneities) is kept and nothing
-   * happens.
-   *
-   * Furthermore, no DoFs belonging to (or
-   * belonging to any descendant of) @p
-   * face_2 get constrained or get marked
-   * as being constrained.
-   *
-   * The flags in the last parameter,
-   * @p component_mask (see @ref GlossComponentMask)
-   * denote which
-   * components of the finite element space
-   * shall be constrained with periodic
-   * boundary conditions. If it is left as
-   * specified by the default value all components are
-   * constrained. If it is different from
-   * the default value, it is assumed that
-   * the number of entries equals the number
-   * of components in the boundary functions
-   * and the finite element, and those
-   * components in the given boundary
-   * function will be used for which the
-   * respective flag was set in the component
-   * mask.
+   * Insert the (algebraic) constraints due to periodic boundary
+   * conditions into a ConstraintMatrix @p constraint_matrix.
    *
-   * @note This function will not work
-   * for DoFHandler objects that are
-   * built on a
-   * parallel::distributed::Triangulation
-   * object.
+   * Given a pair of not necessarily active faces @p face_1 and @p
+   * face_2, this functions constrains all DoFs associated with the boundary
+   * described by @p face_1 to the respective DoFs of the boundary described
+   * by @p face_2. More precisely:
+   *
+   * If @p face_1 and @p face_2 are both active faces it adds the DoFs
+   * of @p face_1 to the list of constrained DoFs in @p constraint_matrix
+   * and adds lines to constrain them to the corresponding values of the
+   * DoFs on @p face_2. This happens on a purely algebraic level, meaning,
+   * the global DoF with (local face) index <tt>i</tt> on @p face_1 gets
+   * constraint to the DoF with (local face) index <tt>i</tt> on @p face_2
+   * (possibly corrected for orientation, see below).
+   *
+   * Otherwise, if @p face_1 and @p face_2 are not active faces, this
+   * function loops recursively over the children of @p face_1 and @p face_2.
+   *
+   * For this to work @p face_1 and @p face_2 must have the same
+   * refinement history, i.e. either @p face_1 and @p face_2 must be active
+   * faces or must be isotropically refined and have the same number of
+   * child faces that recursively obey this rule. (The anisotropic case
+   * is not yet implemented.)
+   *
+   * This routine only constrains DoFs that are not already constrained.
+   * If this routine encounters a DoF that already is constrained (for
+   * instance by Dirichlet boundary conditions), the old setting of the
+   * constraint (dofs the entry is constrained to, inhomogeneities) is
+   * kept and nothing happens.
+   *
+   * Furthermore, no DoFs belonging to (or belonging to any descendant
+   * of) @p face_2 get constrained or get marked as being constrained.
+   *
+   * The flags in the @p component_mask (see @ref GlossComponentMask)
+   * denote which components of the finite element space shall be
+   * constrained with periodic boundary conditions. If it is left as
+   * specified by the default value all components are constrained. If it
+   * is different from the default value, it is assumed that the number
+   * of entries equals the number of components in the boundary functions
+   * and the finite element, and those components in the given boundary
+   * function will be used for which the respective flag was set in the
+   * component mask.
+   *
+   * @p face_orientation, @p face_flip and @p face_rotation describe an
+   * orientation that should be applied to @p face_1 prior to matching and
+   * constraining DoFs. More precisely, this matches local face DoF indices
+   * in the following manner:
+   *
+   * In 2d: <tt>face_orientation</tt> must always be <tt>true</tt>,
+   * <tt>face_rotation</tt> is always <tt>false</tt>, and face_flip has the
+   * meaning of <tt>line_flip</tt>; this implies e.g. for <tt>Q1</tt>:
+   *
+   * @code
+   *
+   * face_orientation = true, face_flip = false, face_rotation = false:
+   *
+   *     face1:           face2:
+   *
+   *     1                1
+   *     |        <-->    |
+   *     0                0
+   *
+   *     Resulting constraints: 0 <-> 0, 1 <-> 1
+   *
+   *     (Numbers denote local face DoF indices.)
+   *
+   *
+   * face_orientation = true, face_flip = true, face_rotation = false:
+   *
+   *     face1:           face2:
+   *
+   *     0                1
+   *     |        <-->    |
+   *     1                0
+   *
+   *     Resulting constraints: 1 <-> 0, 0 <-> 1
+   * @endcode
+   *
+   * And simliarly for the case of Q1 in 3d:
+   *
+   * @code
+   *
+   * face_orientation = true, face_flip = false, face_rotation = false:
+   *
+   *     face1:           face2:
+   *
+   *     2 - 3            2 - 3
+   *     |   |    <-->    |   |
+   *     0 - 1            0 - 1
+   *
+   *     Resulting constraints: 0 <-> 0, 1 <-> 1, 2 <-> 2, 3 <-> 3
+   *
+   *     (Numbers denote local face DoF indices.)
+   *
+   *
+   * face_orientation = false, face_flip = false, face_rotation = false:
+   *
+   *     face1:           face2:
+   *
+   *     1 - 3            2 - 3
+   *     |   |    <-->    |   |
+   *     0 - 2            0 - 1
+   *
+   *     Resulting constraints: 0 <-> 0, 2 <-> 1, 1 <-> 2, 3 <-> 3
+   *
+   *
+   * face_orientation = true, face_flip = true, face_rotation = false:
+   *
+   *     face1:           face2:
+   *
+   *     1 - 0            2 - 3
+   *     |   |    <-->    |   |
+   *     3 - 2            0 - 1
+   *
+   *     Resulting constraints: 3 <-> 0, 2 <-> 1, 1 <-> 2, 0 <-> 3
+   *
+   *
+   * face_orientation = true, face_flip = false, face_rotation = true
+   *
+   *     face1:           face2:
+   *
+   *     0 - 2            2 - 3
+   *     |   |    <-->    |   |
+   *     1 - 3            0 - 1
+   *
+   *     Resulting constraints: 1 <-> 0, 3 <-> 1, 0 <-> 2, 2 <-> 3
+   *
+   * and any combination of that...
+   * @endcode
+   *
+   * More information on the topic can be found in the
+   * @ref GlossFaceOrientation "glossary" article.
+   *
+   * @note This function will not work for DoFHandler objects that are
+   * built on a parallel::distributed::Triangulation object.
+   *
+   * @author Matthias Maier, 2012
    */
   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());
-
-
-  /**
-   * Insert the (algebraic) constraints due
-   * to periodic boundary conditions into
-   * a ConstraintMatrix @p
-   * constraint_matrix.
-   *
-   * This function serves as a high level
-   * interface for the
-   * make_periodicity_constraints function
-   * that takes face_iterator arguments.
-   *
-   * Given a @p direction,
-   * define a 'left' boundary as all
-   * boundary faces belonging to
-   * @p boundary_component with corresponding
-   * unit normal (of the @ref
-   * GlossReferenceCell "reference cell") in
-   * negative @p direction, i.e. all boundary
-   * faces with
-   * <tt>face(2*direction)->at_boundary()==true</tt>.
-   * Analogously, a 'right' boundary
-   * consisting of all faces of @p
-   * boundary_component with unit normal
-   * in positive @p direction,
-   * <tt>face(2*direction+1)->at_boundary()==true</tt>.
-   *
-   * This function tries to match
-   * all faces belonging to the 'left' and
-   * 'right' boundary with the help of an
-   * orthogonal equality relation:
-   * Two faces do match if their vertices
-   * can be transformed into each other by
-   * parallel translation in @p direction.
-   *
-   * If this matching is successfull it
-   * constrains all DoFs associated with the
-   * 'left' boundary to the respective DoFs
-   * of the 'right' boundary.
-   *
-   * This routine only constrains DoFs that
-   * are not already constrained.
-   * If this routine encounters a DoF that
-   * already is constrained (for instance
-   * by Dirichlet boundary conditions),
-   * the old setting of the constraint
-   * (dofs the entry is constrained to,
-   * inhomogeneities) is kept and nothing
-   * happens.
-   *
-   * Furthermore, no DoFs belonging to the
-   * 'right' boundary get constrained or get
-   * marked as being constrained.
-   *
-   * The flags in the last parameter,
-   * @p component_mask (see @ref GlossComponentMask)
-   * denote which
-   * components of the finite element space
-   * shall be constrained with periodic
-   * boundary conditions. If it is left as
-   * specified by the default value all components are
-   * constrained. If it is different from
-   * the default value, it is assumed that
-   * the number of entries equals the number
-   * of components in the boundary functions
-   * and the finite element, and those
-   * components in the given boundary
-   * function will be used for which the
-   * respective flag was set in the component
-   * mask.
+                                const ComponentMask                         &component_mask = ComponentMask(),
+                                bool face_orientation = true,
+                                bool face_flip = false,
+                                bool face_rotation = false);
+
+
+  /**
+   * Insert the (algebraic) constraints due to periodic boundary
+   * conditions into a ConstraintMatrix @p constraint_matrix.
    *
-   * @note This function will not work
-   * for DoFHandler objects that are
-   * built on a
-   * parallel::distributed::Triangulation
-   * object.
+   * This function serves as a high level interface for the
+   * make_periodicity_constraints function that takes face_iterator
+   * arguments.
+   *
+   * Define a 'first' boundary as all boundary faces having boundary_id
+   * @p b_id1 and a 'second' boundary consisting of all faces belonging
+   * to @p b_id2.
+   *
+   * This function tries to match all faces belonging to the first
+   * boundary with faces belonging to the second boundary with the help
+   * of @p orthogonal_equality.
+   *
+   * If this matching is successfull it constrains all DoFs associated
+   * with the 'first' boundary to the respective DoFs of the 'second'
+   * boundary respecting the relative orientation of the two faces.
+   *
+   * This routine only constrains DoFs that are not already constrained.
+   * If this routine encounters a DoF that already is constrained (for
+   * instance by Dirichlet boundary conditions), the old setting of the
+   * constraint (dofs the entry is constrained to, inhomogeneities) is
+   * kept and nothing happens.
+   *
+   * Furthermore, no DoFs belonging to the 'second' boundary get
+   * constrained or get marked as being constrained.
+   *
+   * The flags in the last parameter, @p component_mask (see @ref
+   * GlossComponentMask) denote which components of the finite element space
+   * shall be constrained with periodic boundary conditions. If it is left
+   * as specified by the default value all components are constrained. If
+   * it is different from the default value, it is assumed that the number
+   * of entries equals the number of components in the boundary functions
+   * and the finite element, and those components in the given boundary
+   * function will be used for which the respective flag was set in the
+   * component mask.
+   *
+   * @note This function will not work for DoFHandler objects that are
+   * built on a parallel::distributed::Triangulation object.
+   *
+   * @author Matthias Maier, 2012
    */
   template<typename DH>
   void
-  make_periodicity_constraints (const DH                   &dof_handler,
-                                const types::boundary_id boundary_component,
-                                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());
 
 
   /**
-   * Same as above but with an optional
-   * argument @p offset.
-   * The @p offset is a vector tangential to
-   * the faces that is added to the location
-   * of vertices of the 'left' boundary when
-   * attempting to match them to the
-   * corresponding vertices of the 'right'
-   * boundary. This can be used to implement
-   * conditions such as $u(0,y)=u(1,y+1)$.
+   * Same as above but with an optional argument @p offset.
    *
-   * @note This function will not work
-   * for DoFHandler objects that are
-   * built on a
-   * parallel::distributed::Triangulation
-   * object.
+   * 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 via
+   * @p orthogonal_equality. This can be used to implement conditions such
+   * as $u(0,y)=u(1,y+1)$.
+   *
+   * @note This function will not work for DoFHandler objects that are
+   * built on a parallel::distributed::Triangulation object.
+   *
+   * @author Matthias Maier, 2012
    */
   template<typename DH>
   void
-  make_periodicity_constraints (const DH                       &dof_handler,
-                                const types::boundary_id     boundary_component,
-                                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());
   //@}
 
 
index b9da2481007d90c65891c99435cf61d9361d003a..4e28cc84a7da8043deece2c5b17576ae31899a6f 100644 (file)
@@ -3290,23 +3290,56 @@ 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)
+                                const ComponentMask                         &component_mask,
+                                bool                                        face_orientation,
+                                bool                                        face_flip,
+                                bool                                        face_rotation)
   {
     static const int dim = FaceIterator::AccessorType::dimension;
 
+    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->at_boundary() && face_2->at_boundary(),
            ExcMessage ("Faces for periodicity constraints must be on the boundary"));
 
 
-    // In the case that both faces have
-    // children, we loop over all children
-    // and apply make_periodicty_constrains
-    // recursively:
+    // A lookup table on how to go through the child faces depending on the
+    // orientation:
+
+    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
+      };
+
+    // In the case that both faces have children, we loop over all
+    // children and apply make_periodicty_constrains recursively:
     if (face_1->has_children() && face_2->has_children())
       {
         Assert(face_1->n_children() == GeometryInfo<dim>::max_children_per_face &&
@@ -3315,29 +3348,44 @@ namespace DoFTools
 
         for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
           {
+            // Lookup the index for the second face
+            unsigned int j;
+            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());
+              }
+
             make_periodicity_constraints (face_1->child(i),
-                                          face_2->child(i),
+                                          face_2->child(j),
                                           constraint_matrix,
-                                          component_mask);
+                                          component_mask,
+                                          face_orientation,
+                                          face_flip,
+                                          face_rotation);
           }
         return;
       }
-    // .. otherwise we should be in the case
-    // were both faces are active and have
-    // no children ..
+
+    // .. otherwise we should be in the case were both faces are active
+    // and have no children ..
     Assert (!face_1->has_children() && !face_2->has_children(),
             ExcNotImplemented());
 
     Assert (face_1->n_active_fe_indices() == 1 && face_2->n_active_fe_indices() == 1,
             ExcInternalError());
 
-    // .. then we match the
-    // corresponding DoFs of both faces ..
+    // .. then we match the corresponding DoFs of both faces ..
     const unsigned int face_1_index = face_1->nth_active_fe_index(0);
     const unsigned int face_2_index = face_2->nth_active_fe_index(0);
-    Assert (   face_1->get_fe(face_1_index)
-               == face_2->get_fe(face_1_index),
-               ExcMessage ("Matching periodic cells need to use the same finite element"));
+    Assert (face_1->get_fe(face_1_index) == face_2->get_fe(face_1_index),
+            ExcMessage ("Matching periodic cells need to use the same finite element"));
 
     const dealii::FiniteElement<dim> &fe = face_1->get_fe(face_1_index);
 
@@ -3354,10 +3402,50 @@ namespace DoFTools
     face_1->get_dof_indices(dofs_1, face_1_index);
     face_2->get_dof_indices(dofs_2, face_2_index);
 
-    // .. and constrain them (respecting
-    // component_mask):
+
+    // Well, this is a hack:
+    //
+    // There is no
+    //   face_to_face_index(face_index,
+    //                      face_orientation,
+    //                      face_flip,
+    //                      face_rotation)
+    // function in FiniteElementData, so we have to use
+    //   face_to_cell_index(face_index, face
+    //                      face_orientation,
+    //                      face_flip,
+    //                      face_rotation)
+    // But this will give us an index on a cell - something we cannot work
+    // with. But luckily we can match them back :-]
+
+    std::map<unsigned int, unsigned int> cell_to_rotated_face_index;
+
+    // Build up a cell to face index for face_2:
+    for (unsigned int i = 0; i < dofs_per_face; ++i)
+      {
+        const unsigned int cell_index = fe.
+          face_to_cell_index(i,
+                             0, // It doesn't really matter, just assume
+                                // we're on the first face...
+                             true, false, false // default orientation
+                             );
+        cell_to_rotated_face_index[cell_index] = i;
+      }
+
     for (unsigned int i = 0; i < dofs_per_face; ++i)
       {
+        // Query the correct face_index on face_2 respecting the given
+        // orientation:
+        const unsigned int j =
+          cell_to_rotated_face_index[fe.
+                face_to_cell_index(i,
+                                   0, // It doesn't really matter, just assume
+                                      // we're on the first face...
+                                   face_orientation,
+                                   face_flip,
+                                   face_rotation)];
+
+        // And finally constrain the two DoFs respecting component_mask:
         if ((component_mask.n_selected_components(fe.n_components()) == fe.n_components())
             ||
             (component_mask[fe.face_system_to_component_index(i).first] == true))
@@ -3365,24 +3453,27 @@ namespace DoFTools
             if (!constraint_matrix.is_constrained(dofs_1[i]))
               {
                 constraint_matrix.add_line(dofs_1[i]);
-                constraint_matrix.add_entry(dofs_1[i], dofs_2[i], 1.0);
+                constraint_matrix.add_entry(dofs_1[i], dofs_2[j], 1.0);
               }
           }
       }
   }
 
 
+
   template<typename DH>
   void
   make_periodicity_constraints (const DH                       &dof_handler,
-                                const types::boundary_id     boundary_component,
+                                const types::boundary_id       b_id1,
+                                const types::boundary_id       b_id2,
                                 const int                      direction,
                                 dealii::ConstraintMatrix       &constraint_matrix,
                                 const ComponentMask            &component_mask)
   {
     Tensor<1,DH::space_dimension> dummy;
     make_periodicity_constraints (dof_handler,
-                                  boundary_component,
+                                  b_id1,
+                                  b_id2,
                                   direction,
                                   dummy,
                                   constraint_matrix,
@@ -3390,15 +3481,16 @@ namespace DoFTools
   }
 
 
+
   template<typename DH>
   void
-  make_periodicity_constraints (const DH                       &dof_handler,
-                                const types::boundary_id     boundary_component,
-                                const int                      direction,
-                                dealii::Tensor<1,DH::space_dimension>
-                                &offset,
-                                dealii::ConstraintMatrix       &constraint_matrix,
-                                const ComponentMask            &component_mask)
+  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)
   {
     static const int space_dim = DH::space_dimension;
     Assert (0<=direction && direction<space_dim,
@@ -3410,39 +3502,40 @@ namespace DoFTools
             ExcMessage ("This function can not be used with distributed triangulations."
                         "See the documentation for more information."));
 
-    typedef typename DH::cell_iterator CellIterator;
-
-    // We collect matching periodic cells on
-    // the coarsest level:
-    std::map<CellIterator, CellIterator>
-    matched_cells =
-      GridTools::collect_periodic_cell_pairs(dof_handler.begin(0),
-                                             dof_handler.end(0),
-                                             boundary_component,
-                                             direction,
-                                             offset);
-
-    // And apply the low level
-    // make_periodicity_constraints function
-    // to every matching pair:
-    for (typename std::map<CellIterator, CellIterator>::iterator it = matched_cells.begin();
+    typedef typename DH::face_iterator FaceIterator;
+    typedef std::map<FaceIterator, std::pair<FaceIterator, std::bitset<3> > > FaceMap;
+
+    // Collect matching periodic cells on the coarsest level:
+    FaceMap matched_cells = GridTools::collect_periodic_face_pairs(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();
          it != matched_cells.end(); ++it)
       {
         typedef typename DH::face_iterator FaceIterator;
-        FaceIterator face_1 = it->first->face(2*direction);
-        FaceIterator face_2 = it->second->face(2*direction+1);
+        const FaceIterator &face_1 = it->first;
+        const FaceIterator &face_2 = it->second.first;
+        const std::bitset<3> &orientation = it->second.second;
 
         Assert(face_1->at_boundary() && face_2->at_boundary(),
                ExcInternalError());
 
-        Assert (face_1->boundary_indicator() == boundary_component &&
-                face_2->boundary_indicator() == boundary_component,
+        Assert (face_1->boundary_indicator() == b_id1 &&
+                face_2->boundary_indicator() == b_id2,
                 ExcInternalError());
 
         make_periodicity_constraints(face_1,
                                      face_2,
                                      constraint_matrix,
-                                     component_mask);
+                                     component_mask,
+                                     orientation[0],
+                                     orientation[1],
+                                     orientation[2]);
       }
   }
 
index 1b9a51004342865a6cefb87513b8d1f4cfce6b01..7aa2cb456a8d38c237236d2b71f21b24ebbfacc8 100644 (file)
@@ -326,11 +326,13 @@ for (DH : DOFHANDLERS; deal_II_dimension : DIMENSIONS)
   DoFTools::make_periodicity_constraints (const DH::face_iterator &,
                                           const DH::face_iterator &,
                                           dealii::ConstraintMatrix &,
-                                          const ComponentMask &);
+                                          const ComponentMask &,
+                                          bool, bool, bool);
 
   template
   void
   DoFTools::make_periodicity_constraints(const DH &,
+                                         const types::boundary_id,
                                          const types::boundary_id,
                                          const int,
                                          dealii::ConstraintMatrix &,
@@ -339,6 +341,7 @@ for (DH : DOFHANDLERS; deal_II_dimension : DIMENSIONS)
   template
   void
   DoFTools::make_periodicity_constraints(const DH &,
+                                         const types::boundary_id,
                                          const types::boundary_id,
                                          const int,
                                          dealii::Tensor<1,DH::space_dimension> &,

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.