]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement rotation for vector valued dofs on periodic bc
authorDaniel Arndt <d.arndt@math.uni-goettingen.de>
Wed, 19 Nov 2014 20:14:12 +0000 (21:14 +0100)
committerMatthias Maier <tamiko@kyomu.43-1.org>
Wed, 19 Nov 2014 20:58:47 +0000 (21:58 +0100)
This commit completes the interface for periodic boundary conditions. It
exposes the interpolation matrix (that is internally used to satisfy
hanging node constraints on periodic boundaries) and provides a possibility
to generate it from a rotation matrix and an index for the vector component
in question (in a similar vein as we already do this for different vector
valued finite elements).

Signed-off-by: Matthias Maier <tamiko@kyomu.43-1.org>
include/deal.II/dofs/dof_tools.h
include/deal.II/grid/grid_tools.h
source/dofs/dof_tools_constraints.cc
source/dofs/dof_tools_constraints.inst.in
source/grid/grid_tools.cc
source/grid/grid_tools.inst.in

index 85bdd1a1ff7c53af5c748e6e5956f98f03125976..e7a696146316df5f63db4c7d4049ae7793b6e303 100644 (file)
@@ -848,10 +848,10 @@ namespace DoFTools
    * 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 the finite element. This
-   * can be used to enforce periodicity in only one variable in a system
-   * of equations.
+   * is different from the default value, it is assumed that the number of
+   * entries equals the number of components of the finite element. This
+   * can be used to enforce periodicity in only one variable in a system of
+   * equations.
    *
    * @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
@@ -945,15 +945,21 @@ namespace DoFTools
    * and any combination of that...
    * @endcode
    *
-   * More information on the topic can be found in the
-   * @ref GlossFaceOrientation "glossary" article.
+   * Optionally a matrix @p matrix along with an std::vector @p
+   * first_vector_components can be specified that describes how DoFs on @p
+   * face_1 should be modified prior to constraining to the DoFs of @p
+   * face_2. If the std::vector first_vector_components is non empty the
+   * matrix is interpreted as a rotation matrix that is applied to all
+   * vector valued blocks listed in first_vector_components of the
+   * FESystem.
    *
-   * @note For DoFHandler objects that are built on a
-   * parallel::distributed::Triangulation object
-   * parallel::distributed::Triangulation::add_periodicity has to be called
-   * before.
+   * Detailed information can be found in the @see @ref
+   * GlossPeriodicConstraints "Glossary entry on periodic boundary
+   * conditions".
    *
-   * @author Matthias Maier, 2012
+   * @todo: Reference to soon be written example step and glossary article.
+   *
+   * @author Matthias Maier, 2012, 2014
    */
   template<typename FaceIterator>
   void
@@ -964,7 +970,39 @@ namespace DoFTools
    const ComponentMask                         &component_mask = ComponentMask(),
    const bool                                  face_orientation = true,
    const bool                                  face_flip = false,
-   const bool                                  face_rotation = false);
+   const bool                                  face_rotation = false,
+   const FullMatrix<double>                    &matrix = FullMatrix<double>(),
+   const std::vector<unsigned int>             &first_vector_components = std::vector<unsigned int>());
+
+
+
+  /**
+   * Insert the (algebraic) constraints due to periodic boundary
+   * conditions into a ConstraintMatrix @p constraint_matrix.
+   *
+   * This is the main high level interface for above low level variant of
+   * make_periodicity_constraints. It takes an std::vector @p
+   * periodic_faces as argument and applies above
+   * make_periodicity_constraints on each entry. The std::vector @p
+   * periodic_faces 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.
+   *
+   * @see @ref GlossPeriodicConstraints "Glossary entry on periodic
+   * boundary conditions" for further information.
+   *
+   * @author Daniel Arndt, Matthias Maier, 2013, 2014
+   */
+  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());
 
 
 
@@ -973,8 +1011,7 @@ namespace DoFTools
    * 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.
+   * make_periodicity_constraints function.
    *
    * Define a 'first' boundary as all boundary faces having boundary_id
    * @p b_id1 and a 'second' boundary consisting of all faces belonging
@@ -1004,12 +1041,14 @@ namespace DoFTools
    * function will be used for which the respective flag was set in the
    * component mask.
    *
-   * @note For DoFHandler objects that are built on a
-   * parallel::distributed::Triangulation object
-   * parallel::distributed::Triangulation::add_periodicity has to be called
-   * before.
+   * @note: This function is a convenience wrapper. It internally calls
+   * GridTools::collect_periodic_faces with the supplied paramaters and
+   * feeds the output to above make_periodicity_constraints variant. If you
+   * need more functionality use GridTools::collect_periodic_faces
+   * directly.
    *
-   * @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
+   * @see @ref GlossPeriodicConstraints "Glossary entry on periodic
+   * boundary conditions" for further information.
    *
    * @author Matthias Maier, 2012
    */
@@ -1024,35 +1063,6 @@ namespace DoFTools
    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 '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 For DoFHandler objects that are built on a
-   * parallel::distributed::Triangulation object
-   * parallel::distributed::Triangulation::add_periodicity has to be called
-   * before.
-   *
-   * @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
-   *
-   * @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());
-
 
   /**
    * This compatibility version of make_periodicity_constraints only works
@@ -1069,12 +1079,14 @@ namespace DoFTools
    * meshes with cells not in @ref GlossFaceOrientation
    * "standard orientation".
    *
-   * @note For DoFHandler objects that are built on a
-   * parallel::distributed::Triangulation object
-   * parallel::distributed::Triangulation::add_periodicity has to be called
-   * before.
+   * @note: This function is a convenience wrapper. It internally calls
+   * GridTools::collect_periodic_faces with the supplied paramaters and
+   * feeds the output to above make_periodicity_constraints variant. If you
+   * need more functionality use GridTools::collect_periodic_faces
+   * directly.
    *
-   * @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
+   * @see @ref GlossPeriodicConstraints "Glossary entry on periodic
+   * boundary conditions" for further information.
    */
   template<typename DH>
   void
@@ -1086,53 +1098,54 @@ namespace DoFTools
    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 '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 version of make_periodicity_constraints  will not work on
-   * meshes with cells not in @ref GlossFaceOrientation
-   * "standard orientation".
-   *
-   * @note For DoFHandler objects that are built on a
-   * parallel::distributed::Triangulation object
-   * parallel::distributed::Triangulation::add_periodicity has to be called
-   * before.
-   *
-   * @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
+   * @deprecated This function is deprecated. Use
+   * GridTools::collect_periodic_faces in conjunction with
+   * make_periodicity_constraints instead.
    */
   template<typename DH>
   void
   make_periodicity_constraints
   (const DH                              &dof_handler,
-   const types::boundary_id              b_id,
+   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());
+   const ComponentMask                   &component_mask = ComponentMask()) DEAL_II_DEPRECATED;
+
 
   /**
-   * Same as above but the periodicity information is given by
-   * @p periodic_faces. This std::vector can be created by
-   * GridTools::collect_periodic_faces.
+   * 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 For DoFHandler objects that are built on a
-   * parallel::distributed::Triangulation object
-   * parallel::distributed::Triangulation::add_periodicity has to be called
-   * before.
+   * @note This version of make_periodicity_constraints  will not work on
+   * meshes with cells not in @ref GlossFaceOrientation
+   * "standard orientation".
+   *
+   * @deprecated This function is deprecated. Use
+   * GridTools::collect_periodic_faces in conjunction with
+   * make_periodicity_constraints instead.
    */
   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());
+  (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()) DEAL_II_DEPRECATED;
 
 
   //@}
index 460c9886180de5640e8b91702ee1d3de9dcdc8fc..fc1b9507de3abad522a6829c11227db33f762ff2 100644 (file)
@@ -1052,18 +1052,51 @@ namespace GridTools
   /*@{*/
 
   /**
-   * 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
+   * Data type that provides all information necessary to create
+   * periodicity constraints and a periodic p4est forest with respect to
+   * two 'periodic' cell faces.
    */
   template<typename CellIterator>
   struct PeriodicFacePair
   {
+    /**
+     * The cells associated with the two 'periodic' faces.
+     */
     CellIterator cell[2];
+
+    /**
+     * The local face indices (with respect to the specified cells) of the
+     * two 'periodic' faces.
+     */
     unsigned int face_idx[2];
+
+    /**
+     * The relative orientation of the first face with respect to the
+     * second face as described in orthogonal_equality and
+     * make_periodicity_constraints (and stored as a bitset).
+     */
     std::bitset<3> orientation;
+
+    /**
+     * A matrix that describes how vector valued DoFs of the first face
+     * should be modified prior to constraining to the DoFs of the second
+     * face. If the std::vector first_vector_components is non empty the
+     * matrix is interpreted as a rotation matrix that is applied to all
+     * vector valued blocks listed in first_vector_components of the
+     * FESystem. For more details see make_periodicity_constraints and the
+     * glossary @ref GlossPeriodicConstraints "glossary entry on periodic
+     * boundary conditions".
+     */
+    FullMatrix<double> matrix;
+
+    /**
+     * A vector of unsigned ints pointing to the first components of all
+     * vector valued blocks that should be rotated by matrix.
+     */
+    std::vector<unsigned int> first_vector_components;
   };
 
+
   /**
    * An orthogonal equality test for faces.
    *
@@ -1173,6 +1206,16 @@ namespace GridTools
    * 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)$.
    *
+   * Optionally a rotation matrix @p matrix along with a vector @p
+   * first_vector_components can be specified that describes how vector
+   * valued DoFs of the first face should be modified prior to constraining
+   * to the DoFs of the second face. If @p first_vector_components is non
+   * empty the matrix is interpreted as a rotation matrix that is applied
+   * to all vector valued blocks listet in @p first_vector_components of
+   * the FESystem. For more details see make_periodicity_constraints and
+   * the glossary @ref GlossPeriodicConstraints "glossary entry on periodic
+   * boundary conditions".
+   *
    * @tparam Container A type that satisfies the
    *   requirements of a mesh container (see @ref GlossMeshAsAContainer).
    *
@@ -1186,17 +1229,19 @@ namespace GridTools
    * times with different boundary ids to generate a vector with all periodic
    * pairs.
    *
-   * @author Daniel Arndt, Matthias Maier, 2013
+   * @author Daniel Arndt, Matthias Maier, 2013, 2014
    */
-  template <typename Container>
+  template <typename CONTAINER>
   void
   collect_periodic_faces
-  (const Container &container,
+  (const CONTAINER          &container,
    const types::boundary_id b_id1,
    const types::boundary_id b_id2,
    const int                direction,
-   std::vector<PeriodicFacePair<typename Container::cell_iterator> > &matched_pairs,
-   const dealii::Tensor<1,Container::space_dimension> &offset = dealii::Tensor<1,Container::space_dimension>());
+   std::vector<PeriodicFacePair<typename CONTAINER::cell_iterator> > &matched_pairs,
+   const Tensor<1,CONTAINER::space_dimension> &offset = dealii::Tensor<1,CONTAINER::space_dimension>(),
+   const FullMatrix<double> &matrix = FullMatrix<double>(IdentityMatrix(CONTAINER::space_dimension)),
+   const std::vector<unsigned int> &first_vector_components = std::vector<unsigned int>());
 
 
   /**
@@ -1213,6 +1258,16 @@ namespace GridTools
    * This function will collect periodic face pairs on the coarsest mesh level
    * and add them to @p matched_pairs leaving the original contents intact.
    *
+   * Optionally a rotation matrix @p matrix along with a vector @p
+   * first_vector_components can be specified that describes how vector
+   * valued DoFs of the first face should be modified prior to constraining
+   * to the DoFs of the second face. If @p first_vector_components is non
+   * empty the matrix is interpreted as a rotation matrix that is applied
+   * to all vector valued blocks listet in @p first_vector_components of
+   * the FESystem. For more details see make_periodicity_constraints and
+   * the glossary @ref GlossPeriodicConstraints "glossary entry on periodic
+   * boundary conditions".
+   *
    * @tparam Container A type that satisfies the
    *   requirements of a mesh container (see @ref GlossMeshAsAContainer).
    *
@@ -1220,16 +1275,19 @@ namespace GridTools
    * meshes with cells not in @ref GlossFaceOrientation
    * "standard orientation".
    *
-   * @author Daniel Arndt, Matthias Maier, 2013
+   * @author Daniel Arndt, Matthias Maier, 2013, 2014
    */
-  template <typename Container>
+  template <typename CONTAINER>
   void
   collect_periodic_faces
-  (const Container         &container,
+  (const CONTAINER          &container,
    const types::boundary_id b_id,
    const int                direction,
-   std::vector<PeriodicFacePair<typename Container::cell_iterator> > &matched_pairs,
-   const dealii::Tensor<1,Container::space_dimension> &offset = dealii::Tensor<1,Container::space_dimension>());
+   std::vector<PeriodicFacePair<typename CONTAINER::cell_iterator> > &matched_pairs,
+   const dealii::Tensor<1,CONTAINER::space_dimension> &offset = dealii::Tensor<1,CONTAINER::space_dimension>(),
+   const FullMatrix<double> &matrix = FullMatrix<double>(),
+   const std::vector<unsigned int> &first_vector_components = std::vector<unsigned int>());
+
 
   /*@}*/
   /**
index 56aaa3616ebce58c027189a9ae9153db92124c97..cf52dcac5e25b8b7d362322d956b73ed7bed9834 100644 (file)
@@ -13,6 +13,7 @@
 //
 // ---------------------------------------------------------------------
 
+#include <deal.II/base/std_cxx1x/array.h>
 #include <deal.II/base/thread_management.h>
 #include <deal.II/base/table.h>
 #include <deal.II/base/template_constraints.h>
@@ -1711,7 +1712,7 @@ namespace DoFTools
         {
           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),
+          Assert(face_1->get_fe(face_1_index) == face_2->get_fe(face_2_index),
                  ExcMessage ("Matching periodic cells need to use the same finite element"));
 
           const FiniteElement<dim, spacedim> &fe = face_1->get_fe(face_1_index);
@@ -1834,8 +1835,8 @@ namespace DoFTools
                       const unsigned int j =
                         cell_to_rotated_face_index[fe.face_to_cell_index(identity_constraint_target,
                                                                          0, /* It doesn't really matter, just assume
-                           * we're on the first face...
-                           */
+                                                                             * we're on the first face...
+                                                                             */
                                                                          face_orientation, face_flip, face_rotation)];
 
                       // if the two aren't already identity constrained (whichever way
@@ -1856,10 +1857,8 @@ namespace DoFTools
                           // 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(jj, 0, /* It doesn't really matter, just assume
-                               * we're on the first face...
-                               */
-                                                                             face_orientation, face_flip, face_rotation)];
+                            cell_to_rotated_face_index[fe.face_to_cell_index
+                                                       (jj, 0, face_orientation, face_flip, face_rotation)];
 
                           // And finally constrain the two DoFs respecting component_mask:
                           if (transformation(i,jj) != 0)
@@ -1873,6 +1872,9 @@ namespace DoFTools
   }
 
 
+  // Implementation of the low level interface:
+
+
   template <typename FaceIterator>
   void
   make_periodicity_constraints (const FaceIterator                          &face_1,
@@ -1881,7 +1883,9 @@ namespace DoFTools
                                 const ComponentMask                         &component_mask,
                                 const bool                                   face_orientation,
                                 const bool                                   face_flip,
-                                const bool                                   face_rotation)
+                                const bool                                   face_rotation,
+                                const FullMatrix<double>                    &matrix,
+                                const std::vector<unsigned int>             &first_vector_components)
   {
     static const int dim = FaceIterator::AccessorType::dimension;
 
@@ -1907,6 +1911,24 @@ namespace DoFTools
     Assert(face_1->at_boundary() && face_2->at_boundary(),
            ExcMessage ("Faces for periodicity constraints must be on the boundary"));
 
+    Assert(first_vector_components.empty()
+           || (matrix.size(0) == (int)dim && matrix.size(1) == (int)dim),
+           ExcMessage ("first_vector_components is nonempty, so matrix must be a rotation "
+                       "matrix exactly of size dim"));
+
+#ifdef DEBUG
+    if (!face_1->has_children())
+      {
+        Assert(face_1->n_active_fe_indices()==1, ExcInternalError());
+        const unsigned int n_dofs = face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face;
+        Assert(first_vector_components.size() != 0 ||
+               ( (matrix.size(0) == 0 && matrix.size(1) == 0) ||
+                 (matrix.size(0) == n_dofs && matrix.size(1) == n_dofs) ||
+                 (matrix.size(0) == (int)dim && matrix.size(1) == (int)dim) ),
+               ExcMessage ("first_vector_components is empty, so matrix must have either "
+                           "size 0 or dim or the size must be equal to the number of DoFs on the face"));
+      }
+#endif
 
     // A lookup table on how to go through the child faces depending on the
     // orientation:
@@ -1967,30 +1989,156 @@ namespace DoFTools
                                           component_mask,
                                           face_orientation,
                                           face_flip,
-                                          face_rotation);
+                                          face_rotation,
+                                          matrix,
+                                          first_vector_components);
           }
       }
     else
       // otherwise at least one of the two faces is active and
       // we need to enter the constraints
       {
+        // Build up the transformation matrix:
+        FullMatrix<double> transformation;
+        const unsigned int n_dofs = face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face;
+
+        if (matrix.size(0) == n_dofs)
+          transformation = matrix;
+        else if (!first_vector_components.empty() && matrix.size(0) == (int)dim)
+          {
+            // The matrix describes a rotation and we have to build a
+            // transformation matrix, we assume that for a 0° rotation
+            // we would have to build the identity matrix
+
+            const FiniteElement<dim> &fe1
+              = face_1->get_fe(face_1->nth_active_fe_index(0));
+
+            Quadrature<dim-1> quadrature (fe1.get_unit_face_support_points());
+
+            FEFaceValues<dim> fe_face_values
+            (fe1, quadrature, update_q_points);
+
+            // have an array that stores the location of each vector-dof tuple
+            // we want to rotate.
+            typedef std_cxx1x::array<unsigned int, dim> DoFTuple;
+
+            FullMatrix<double> rot_transformation
+              = IdentityMatrix(face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face);
+
+            for (unsigned int i=0; i<fe1.dofs_per_face; ++i)
+              {
+                std::vector<unsigned int>::const_iterator comp_it
+                  = std::find (first_vector_components.begin(),
+                               first_vector_components.end(),
+                               fe1.face_system_to_component_index(i).first);
+                if (comp_it != first_vector_components.end())
+                  {
+                    const unsigned int first_vector_component = *comp_it;
+
+                    // find corresponding other components of vector
+                    DoFTuple vector_dofs;
+                    vector_dofs[0] = i;
+
+                    Assert(*comp_it+dim<=fe1.n_components(),
+                           ExcMessage("Error: the finite element does not have enough components "
+                                      "to define rotated periodic boundaries."));
+
+                    for (unsigned int k=0; k<fe1.dofs_per_face; ++k)
+                      if ((k != i)
+                          &&
+                          (quadrature.point(k) == quadrature.point(i))
+                          &&
+                          (fe1.face_system_to_component_index(k).first >=
+                           first_vector_component)
+                          &&
+                          (fe1.face_system_to_component_index(k).first <
+                           first_vector_component + dim))
+                        vector_dofs[fe1.face_system_to_component_index(k).first -
+                                    first_vector_component]
+                          = k;
+
+                    for (int i=0; i<dim; ++i)
+                      {
+                        rot_transformation[vector_dofs[i]][vector_dofs[i]]=0.;
+                        for (int j=0; j<dim; ++j)
+                          rot_transformation[vector_dofs[i]][vector_dofs[j]]=matrix[i][j];
+                      }
+                  }
+              }
+
+            transformation = rot_transformation;
+          }
+        else
+          transformation = IdentityMatrix(n_dofs);
+
         if (face_2->has_children() == false)
-          set_periodicity_constraints(face_2, face_1,
-                                      FullMatrix<double>(IdentityMatrix(face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face)),
-                                      constraint_matrix,
-                                      component_mask,
-                                      face_orientation, face_flip, face_rotation);
+          {
+            FullMatrix<double> inverse(transformation.size(0));
+            inverse.invert(transformation);
+            set_periodicity_constraints(face_2, face_1,
+                                        inverse,
+                                        constraint_matrix,
+                                        component_mask,
+                                        face_orientation, face_flip, face_rotation);
+          }
         else
-          set_periodicity_constraints(face_1, face_2,
-                                      FullMatrix<double>(IdentityMatrix(face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face)),
-                                      constraint_matrix,
-                                      component_mask,
-                                      face_orientation, face_flip, face_rotation);
+          {
+            set_periodicity_constraints(face_1, face_2,
+                                        transformation,
+                                        constraint_matrix,
+                                        component_mask,
+                                        face_orientation, face_flip, face_rotation);
+          }
       }
   }
 
 
 
+  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();
+
+    // Loop over all periodic faces...
+    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]);
+
+        Assert(face_1->at_boundary() && face_2->at_boundary(),
+               ExcInternalError());
+
+        Assert (face_1 != face_2,
+                ExcInternalError());
+
+        // ... and apply the low level make_periodicity_constraints function to
+        // every matching pair:
+        make_periodicity_constraints(face_1,
+                                     face_2,
+                                     constraint_matrix,
+                                     component_mask,
+                                     it->orientation[0],
+                                     it->orientation[1],
+                                     it->orientation[2],
+                                     it->matrix,
+                                     it->first_vector_components);
+      }
+  }
+
+
+  // High level interface variants:
+
+
   template<typename DH>
   void
   make_periodicity_constraints (const DH                       &dof_handler,
@@ -2093,47 +2241,6 @@ 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();
-
-
-    // And apply the low level make_periodicity_constraints function to
-    // every matching pair:
-    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]);
-
-        Assert(face_1->at_boundary() && face_2->at_boundary(),
-               ExcInternalError());
-
-        Assert (face_1 != face_2,
-                ExcInternalError());
-
-        make_periodicity_constraints(face_1,
-                                     face_2,
-                                     constraint_matrix,
-                                     component_mask,
-                                     it->orientation[0],
-                                     it->orientation[1],
-                                     it->orientation[2]);
-      }
-  }
-
-
-
   namespace internal
   {
     namespace Assembler
index 4b2e4e8e21bd19b7c1429374e5fae8991863bebe..3017a6d0b0fd93f1c7a60274a63e81aeb3e90d42 100644 (file)
@@ -30,7 +30,17 @@ for (DH : DOFHANDLERS; deal_II_dimension : DIMENSIONS)
                                           const DH::face_iterator &,
                                           dealii::ConstraintMatrix &,
                                           const ComponentMask &,
-                                          bool, bool, bool);
+                                          bool, bool, bool,
+                                          const FullMatrix<double> &,
+                                          const std::vector<unsigned int> &);
+
+  template
+  void
+  DoFTools::make_periodicity_constraints<DH>
+  (const std::vector<GridTools::PeriodicFacePair<DH::cell_iterator> > &,
+   dealii::ConstraintMatrix &,
+   const ComponentMask &);
+
 
   template
   void
@@ -67,14 +77,6 @@ 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<DH::cell_iterator> > &,
-   dealii::ConstraintMatrix &,
-   const ComponentMask &);
-  
 #endif
 }
 
index 47570805b758ee5925fd3774f13afdce9dae769a..72a48e641d1261751afac7d9d34b43143ce3e5b3 100644 (file)
@@ -2367,14 +2367,15 @@ next_cell:
    * An orthogonal equality test for points:
    *
    * point1 and point2 are considered equal, if
-   *    (point1 + offset) - point2
+   *   matrix.(point1 + offset) - point2
    * is parallel to the unit vector in <direction>
    */
   template<int spacedim>
   inline bool orthogonal_equality (const dealii::Point<spacedim> &point1,
                                    const dealii::Point<spacedim> &point2,
                                    const int                     direction,
-                                   const dealii::Tensor<1,spacedim> &offset)
+                                   const dealii::Tensor<1,spacedim> &offset,
+                                   const FullMatrix<double>      &matrix)
   {
     Assert (0<=direction && direction<spacedim,
             ExcIndexRange (direction, 0, spacedim));
@@ -2384,9 +2385,15 @@ next_cell:
         if (i == direction)
           continue;
 
-        if (fabs(point1(i) + offset[i] - point2(i)) > 1.e-10)
+        double transformed_p1_comp=0.;
+
+        for (int j = 0; j < spacedim; ++j)
+          transformed_p1_comp += matrix(i,j)*point1(j) + offset[j];
+
+        if (fabs(transformed_p1_comp - point2(i)) > 1.e-10)
           return false;
       }
+
     return true;
   }
 
@@ -2475,7 +2482,8 @@ next_cell:
                        const FaceIterator &face1,
                        const FaceIterator &face2,
                        const int          direction,
-                       const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset)
+                       const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset,
+                       const FullMatrix<double> &matrix)
   {
     static const int dim = FaceIterator::AccessorType::dimension;
 
@@ -2495,7 +2503,7 @@ next_cell:
              it++)
           {
             if (orthogonal_equality(face1->vertex(i),face2->vertex(*it),
-                                    direction, offset))
+                                    direction, offset, matrix))
               {
                 matching[i] = *it;
                 face2_vertices.erase(it);
@@ -2518,11 +2526,12 @@ next_cell:
   orthogonal_equality (const FaceIterator &face1,
                        const FaceIterator &face2,
                        const int          direction,
-                       const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset)
+                       const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset,
+                       const FullMatrix<double> &matrix)
   {
     // Call the function above with a dummy orientation array
     std::bitset<3> dummy;
-    return orthogonal_equality (dummy, face1, face2, direction, offset);
+    return orthogonal_equality (dummy, face1, face2, direction, offset, matrix);
   }
 
 
@@ -2535,9 +2544,11 @@ next_cell:
   match_periodic_face_pairs
   (std::set<std::pair<CellIterator, unsigned int> > &pairs1,
    std::set<std::pair<typename identity<CellIterator>::type, unsigned int> > &pairs2,
-   const int direction,
-   std::vector<PeriodicFacePair<CellIterator> > &matched_pairs,
-   const dealii::Tensor<1,CellIterator::AccessorType::space_dimension> &offset)
+   const int                                        direction,
+   std::vector<PeriodicFacePair<CellIterator> >     &matched_pairs,
+   const dealii::Tensor<1,CellIterator::AccessorType::space_dimension> &offset,
+   const FullMatrix<double>                         &matrix,
+   const std::vector<unsigned int>                  &first_vector_components)
   {
     static const int space_dim = CellIterator::AccessorType::space_dimension;
     Assert (0<=direction && direction<space_dim,
@@ -2563,13 +2574,14 @@ next_cell:
             if (GridTools::orthogonal_equality(orientation,
                                                cell1->face(face_idx1),
                                                cell2->face(face_idx2),
-                                               direction, offset))
+                                               direction, offset,
+                                               matrix))
               {
                 // 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};
+                = {{cell1, cell2},{face_idx1, face_idx2}, orientation, matrix, first_vector_components};
                 matched_pairs.push_back(matched_face);
                 pairs2.erase(it2);
                 ++n_matches;
@@ -2588,12 +2600,14 @@ next_cell:
   template<typename CONTAINER>
   void
   collect_periodic_faces
-  (const CONTAINER          &container,
-   const types::boundary_id b_id1,
-   const types::boundary_id b_id2,
-   const int                direction,
+  (const CONTAINER                            &container,
+   const types::boundary_id                   b_id1,
+   const types::boundary_id                   b_id2,
+   const int                                  direction,
    std::vector<PeriodicFacePair<typename CONTAINER::cell_iterator> > &matched_pairs,
-   const dealii::Tensor<1,CONTAINER::space_dimension> &offset)
+   const Tensor<1,CONTAINER::space_dimension> &offset,
+   const FullMatrix<double>                   &matrix,
+   const std::vector<unsigned int>            &first_vector_components)
   {
     static const int dim = CONTAINER::dimension;
     static const int space_dim = CONTAINER::space_dimension;
@@ -2632,7 +2646,8 @@ next_cell:
             ExcMessage ("Unmatched faces on periodic boundaries"));
 
     // and call match_periodic_face_pairs that does the actual matching:
-    match_periodic_face_pairs(pairs1, pairs2, direction, matched_pairs, offset);
+    match_periodic_face_pairs(pairs1, pairs2, direction, matched_pairs,
+                              offset, matrix, first_vector_components);
   }
 
 
@@ -2640,11 +2655,13 @@ next_cell:
   template<typename CONTAINER>
   void
   collect_periodic_faces
-  (const CONTAINER          &container,
-   const types::boundary_id b_id,
-   const int                direction,
+  (const CONTAINER                            &container,
+   const types::boundary_id                   b_id,
+   const int                                  direction,
    std::vector<PeriodicFacePair<typename CONTAINER::cell_iterator> > &matched_pairs,
-   const dealii::Tensor<1,CONTAINER::space_dimension> &offset)
+   const Tensor<1,CONTAINER::space_dimension> &offset,
+   const FullMatrix<double>                   &matrix,
+   const std::vector<unsigned int>            &first_vector_components)
   {
     static const int dim = CONTAINER::dimension;
     static const int space_dim = CONTAINER::space_dimension;
@@ -2690,7 +2707,8 @@ next_cell:
 #endif
 
     // and call match_periodic_face_pairs that does the actual matching:
-    match_periodic_face_pairs(pairs1, pairs2, direction, matched_pairs, offset);
+    match_periodic_face_pairs(pairs1, pairs2, direction, matched_pairs,
+                              offset, matrix, first_vector_components);
 
 #ifdef DEBUG
     //check for standard orientation
index b949d29f2d19b3a86d882bbfc30a4dc4005c35ee..1c2ee8d9ec836b5890e7cef96d363db986aeb6f5 100644 (file)
@@ -302,26 +302,30 @@ for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II
                                                        const X::active_face_iterator&,
                                                        const X::active_face_iterator&,
                                                        const int,
-                                                       const Tensor<1,deal_II_space_dimension> &);
+                                                       const Tensor<1,deal_II_space_dimension> &,
+                                                       const FullMatrix<double> &);
 
     template
     bool orthogonal_equality<X::face_iterator> (std::bitset<3> &,
                                                 const X::face_iterator&,
                                                 const X::face_iterator&,
                                                 const int,
-                                                const Tensor<1,deal_II_space_dimension> &);
+                                                const Tensor<1,deal_II_space_dimension> &,
+                                                const FullMatrix<double> &);
 
     template
     bool orthogonal_equality<X::active_face_iterator> (const X::active_face_iterator&,
                                                        const X::active_face_iterator&,
                                                        const int,
-                                                       const Tensor<1,deal_II_space_dimension> &);
+                                                       const Tensor<1,deal_II_space_dimension> &,
+                                                       const FullMatrix<double> &);
 
     template
     bool orthogonal_equality<X::face_iterator> (const X::face_iterator&,
                                                 const X::face_iterator&,
                                                 const int,
-                                                const Tensor<1,deal_II_space_dimension> &);
+                                                const Tensor<1,deal_II_space_dimension> &,
+                                                const FullMatrix<double> &);
 
     #if deal_II_dimension >= 2
 
@@ -331,14 +335,18 @@ for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II
                                       const types::boundary_id,
                                       const int,
                                       std::vector<PeriodicFacePair<X::cell_iterator> > &,
-                                      const Tensor<1,X::space_dimension> &);
+                                      const Tensor<1,X::space_dimension> &,
+                                      const FullMatrix<double> &,
+                                      const std::vector<unsigned int> &);
 
       template
       void collect_periodic_faces<X> (const X &,
                                       const types::boundary_id,
                                       const int,
                                       std::vector<PeriodicFacePair<X::cell_iterator> > &,
-                                      const Tensor<1,X::space_dimension> &);
+                                      const Tensor<1,X::space_dimension> &,
+                                      const FullMatrix<double> &,
+                                      const std::vector<unsigned int> &);
 
     #endif
 
@@ -352,7 +360,6 @@ for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS
    #if deal_II_dimension >= 2
 
      namespace GridTools \{
-
       template
       void
       collect_periodic_faces<parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension> >
@@ -361,7 +368,9 @@ for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS
                                   const types::boundary_id,
                                   const int,
                                   std::vector<PeriodicFacePair<parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension>::cell_iterator> > &,
-                                  const Tensor<1,parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension>::space_dimension> &);
+                                  const Tensor<1,parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension>::space_dimension> &,
+                                  const FullMatrix<double> &,
+                                  const std::vector<unsigned int> &);
 
       template
       void
@@ -370,8 +379,9 @@ for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS
                                   const types::boundary_id,
                                   const int,
                                   std::vector<PeriodicFacePair<parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension>::cell_iterator> > &,
-                                  const Tensor<1,parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension>::space_dimension> &);
-
+                                  const Tensor<1,parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension>::space_dimension> &,
+                                  const FullMatrix<double> &,
+                                  const std::vector<unsigned int> &);
      \}
    #endif
 #endif

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.