From 2c481fd3cb26c73012c00aaf815a9fa18070eb69 Mon Sep 17 00:00:00 2001 From: maier Date: Mon, 21 May 2012 22:22:59 +0000 Subject: [PATCH] Add DoFTools::make_periodicity_constraints git-svn-id: https://svn.dealii.org/trunk@25533 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/dofs/dof_tools.h | 183 +++++++++++++++++++++++ deal.II/source/dofs/dof_tools.cc | 144 ++++++++++++++++++ deal.II/source/dofs/dof_tools.inst.in | 30 ++++ 3 files changed, 357 insertions(+) diff --git a/deal.II/include/deal.II/dofs/dof_tools.h b/deal.II/include/deal.II/dofs/dof_tools.h index 3db11e7135..00728695da 100644 --- a/deal.II/include/deal.II/dofs/dof_tools.h +++ b/deal.II/include/deal.II/dofs/dof_tools.h @@ -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 + void + make_periodicity_constraints (const FaceIterator &face_1, + const typename identity::type &face_2, + dealii::ConstraintMatrix &constraint_matrix, + const std::vector &component_mask = std::vector()); + + + /** + * 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 + * face(2*direction)->at_boundary()==true. + * Analogously, a 'right' boundary + * consisting of all faces of @p + * boundary_component with unit normal + * in positive @p direction, + * face(2*direction+1)->at_boundary()==true. + * + * 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 + 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 &component_mask = std::vector()); + + + /** + * 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 + 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 &component_mask = std::vector()); + //@} + + /** * Take a vector of values which live on * cells (e.g. an error per cell) and diff --git a/deal.II/source/dofs/dof_tools.cc b/deal.II/source/dofs/dof_tools.cc index b6efd20146..ced4db6d8e 100644 --- a/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/source/dofs/dof_tools.cc @@ -3298,6 +3298,150 @@ namespace DoFTools } + template + void + make_periodicity_constraints (const FaceIterator &face_1, + const typename identity::type &face_2, + dealii::ConstraintMatrix &constraint_matrix, + const std::vector &component_mask = std::vector()) + { + 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::max_children_per_face && + face_2->n_children() == GeometryInfo::max_children_per_face, + ExcNotImplemented()); + + for (unsigned int i = 0; i < GeometryInfo::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 &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 dofs_1(dofs_per_face); + std::vector 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 + 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 &component_mask = std::vector()) + { + Tensor<1,DH::space_dimension> dummy; + make_periodicity_constraints (dof_handler, + boundary_component, + direction, + dummy, + constraint_matrix, + component_mask); + } + + + template + 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 &component_mask = std::vector()) + { + static const int space_dim = DH::space_dimension; + Assert (0<=direction && direction + 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 { diff --git a/deal.II/source/dofs/dof_tools.inst.in b/deal.II/source/dofs/dof_tools.inst.in index f0b3fe8f53..ad03511dc2 100644 --- a/deal.II/source/dofs/dof_tools.inst.in +++ b/deal.II/source/dofs/dof_tools.inst.in @@ -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 &); + + template + void + DoFTools::make_periodicity_constraints(const DH &, + const types::boundary_id_t, + const int, + dealii::ConstraintMatrix &, + const std::vector &); + + 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 &); +} + + + for (deal_II_dimension : DIMENSIONS) { template -- 2.39.5