]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add DoFTools::make_periodicity_constraints
authorMatthias Maier <tamiko@kyomu.43-1.org>
Mon, 21 May 2012 22:22:59 +0000 (22:22 +0000)
committerMatthias Maier <tamiko@kyomu.43-1.org>
Mon, 21 May 2012 22:22:59 +0000 (22:22 +0000)
git-svn-id: https://svn.dealii.org/trunk@25533 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 3db11e7135fd8be027c2e63f81eeff27943189f5..00728695da2b4e269fa56da9226584815505f115 100644 (file)
@@ -926,6 +926,189 @@ namespace DoFTools
                                  ConstraintMatrix &constraints);
                                    //@}
 
+
+                                   /**
+                                    * @name Periodic Boundary Conditions
+                                    * @{
+                                    */
+
+                                   /**
+                                    * 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 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 (i.e. an
+                                    * empty array), all components are
+                                    * constrainted. 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.
+                                    */
+  template<typename FaceIterator>
+  void
+  make_periodicity_constraints (const FaceIterator                          &face_1,
+                                const typename identity<FaceIterator>::type &face_2,
+                                dealii::ConstraintMatrix                    &constraint_matrix,
+                                const std::vector<bool>                     &component_mask = std::vector<bool>());
+
+
+                                   /**
+                                    * 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 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 (i.e. an
+                                    * empty array), all components are
+                                    * constrainted. 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.
+                                    */
+  template<typename DH>
+  void
+  make_periodicity_constraints (const DH                   &dof_handler,
+                                const types::boundary_id_t boundary_component,
+                                const int                  direction,
+                                dealii::ConstraintMatrix   &constraint_matrix,
+                                const std::vector<bool>    &component_mask = std::vector<bool>());
+
+
+                                   /**
+                                    * 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)$.
+                                    */
+  template<typename DH>
+  void
+  make_periodicity_constraints (const DH                       &dof_handler,
+                                const types::boundary_id_t     boundary_component,
+                                const int                      direction,
+                                dealii::Tensor<1,DH::space_dimension>
+                                                               &offset,
+                                dealii::ConstraintMatrix       &constraint_matrix,
+                                const std::vector<bool>        &component_mask = std::vector<bool>());
+                                   //@}
+
+
                                    /**
                                     * Take a vector of values which live on
                                     * cells (e.g. an error per cell) and
index b6efd20146b87cea6012d84b8e93a22851d2f5f0..ced4db6d8ea661f09161a7f00533076fc3bed537 100644 (file)
@@ -3298,6 +3298,150 @@ 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 std::vector<bool>                     &component_mask = std::vector<bool>())
+  {
+    static const int dim = FaceIterator::AccessorType::dimension;
+
+    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 applu make_periodicty_constrains
+                                   // recursively:
+    if (face_1->has_children() && face_2->has_children()) {
+      Assert(face_1->n_children() == GeometryInfo<dim>::max_children_per_face &&
+             face_2->n_children() == GeometryInfo<dim>::max_children_per_face,
+             ExcNotImplemented());
+
+      for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i) {
+        make_periodicity_constraints (face_1->child(i),
+                                      face_2->child(i),
+                                      constraint_matrix,
+                                      component_mask);
+      }
+      return;
+    }
+                                   // .. 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 ..
+    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"));
+
+    const dealii::FiniteElement<dim> &fe = face_1->get_fe(face_1_index);
+
+    Assert (component_mask.size() == 0 ||
+            component_mask.size() == fe.n_components(),
+            ExcMessage ("The number of components in the mask has to be either "
+                        "zero or equal to the number of components in the finite "
+                        "element."));
+
+    const unsigned int dofs_per_face = fe.dofs_per_face;
+
+    std::vector<unsigned int> dofs_1(dofs_per_face);
+    std::vector<unsigned int> dofs_2(dofs_per_face);
+
+    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):
+    for (unsigned int i = 0; i < dofs_per_face; ++i) {
+      if (component_mask.size() == 0 ||
+          component_mask[fe.face_system_to_component_index(i).first]) {
+        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);
+        }
+      }
+    }
+  }
+
+
+  template<typename DH>
+  void
+  make_periodicity_constraints (const DH                       &dof_handler,
+                                const types::boundary_id_t     boundary_component,
+                                const int                      direction,
+                                dealii::ConstraintMatrix       &constraint_matrix,
+                                const std::vector<bool>        &component_mask = std::vector<bool>())
+  {
+    Tensor<1,DH::space_dimension> dummy;
+    make_periodicity_constraints (dof_handler,
+                                  boundary_component,
+                                  direction,
+                                  dummy,
+                                  constraint_matrix,
+                                  component_mask);
+  }
+
+
+  template<typename DH>
+  void
+  make_periodicity_constraints (const DH                       &dof_handler,
+                                const types::boundary_id_t     boundary_component,
+                                const int                      direction,
+                                dealii::Tensor<1,DH::space_dimension>
+                                                               &offset,
+                                dealii::ConstraintMatrix       &constraint_matrix,
+                                const std::vector<bool>        &component_mask = std::vector<bool>())
+  {
+    static const int space_dim = DH::space_dimension;
+    Assert (0<=direction && direction<space_dim,
+            ExcIndexRange (direction, 0, space_dim));
+
+    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 (auto 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);
+
+      Assert(face_1->at_boundary() && face_2->at_boundary(),
+             ExcInternalError());
+
+      Assert (face_1->boundary_indicator() == boundary_component &&
+              face_2->boundary_indicator() == boundary_component,
+              ExcInternalError());
+
+      make_periodicity_constraints(face_1,
+                                   face_2,
+                                   constraint_matrix,
+                                   component_mask);
+    }
+  }
+
+
 
   namespace internal
   {
index f0b3fe8f53ec044ea5ba7015e194b90e17e97e29..ad03511dc23e44526f91f559c05f73657966d4da 100644 (file)
@@ -316,6 +316,36 @@ for (DH : DOFHANDLERS; deal_II_dimension : DIMENSIONS)
 }
 
 
+
+for (DH : DOFHANDLERS; deal_II_dimension : DIMENSIONS)
+{
+  template
+  void
+  DoFTools::make_periodicity_constraints (const DH::face_iterator &,
+                                          const DH::face_iterator &,
+                                          dealii::ConstraintMatrix &,
+                                          const std::vector<bool> &);
+
+  template
+  void
+  DoFTools::make_periodicity_constraints(const DH &,
+                                         const types::boundary_id_t,
+                                         const int,
+                                         dealii::ConstraintMatrix &,
+                                         const std::vector<bool> &);
+
+  template
+  void
+  DoFTools::make_periodicity_constraints(const DH &,
+                                         const types::boundary_id_t,
+                                         const int,
+                                         dealii::Tensor<1,DH::space_dimension> &,
+                                         dealii::ConstraintMatrix &,
+                                         const std::vector<bool> &);
+}
+
+
+
 for (deal_II_dimension : DIMENSIONS)
 {
   template

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.