From b3efbf84539383fd232c9a9b84cf8d315e207386 Mon Sep 17 00:00:00 2001 From: maier Date: Sat, 1 Dec 2012 17:17:25 +0000 Subject: [PATCH] Bugfix: Make DoFTools::make_periodicity_constraints orientation aware 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 | 382 +++++++++++++---------- deal.II/source/dofs/dof_tools.cc | 189 ++++++++--- deal.II/source/dofs/dof_tools.inst.in | 5 +- 3 files changed, 354 insertions(+), 222 deletions(-) diff --git a/deal.II/include/deal.II/dofs/dof_tools.h b/deal.II/include/deal.II/dofs/dof_tools.h index c39b590c47..5a9b032536 100644 --- a/deal.II/include/deal.II/dofs/dof_tools.h +++ b/deal.II/include/deal.II/dofs/dof_tools.h @@ -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 i on @p face_1 gets + * constraint to the DoF with (local face) index i 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: face_orientation must always be true, + * face_rotation is always false, and face_flip has the + * meaning of line_flip; this implies e.g. for Q1: + * + * @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 void make_periodicity_constraints (const FaceIterator &face_1, const typename identity::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 - * 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 (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 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 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()); //@} diff --git a/deal.II/source/dofs/dof_tools.cc b/deal.II/source/dofs/dof_tools.cc index b9da248100..4e28cc84a7 100644 --- a/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/source/dofs/dof_tools.cc @@ -3290,23 +3290,56 @@ namespace DoFTools } + template void make_periodicity_constraints (const FaceIterator &face_1, const typename identity::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::max_children_per_face && @@ -3315,29 +3348,44 @@ namespace DoFTools for (unsigned int i = 0; i < GeometryInfo::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 &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 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 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 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 - 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::iterator it = matched_cells.begin(); + typedef typename DH::face_iterator FaceIterator; + typedef std::map > > 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]); } } diff --git a/deal.II/source/dofs/dof_tools.inst.in b/deal.II/source/dofs/dof_tools.inst.in index 1b9a510043..7aa2cb456a 100644 --- a/deal.II/source/dofs/dof_tools.inst.in +++ b/deal.II/source/dofs/dof_tools.inst.in @@ -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> &, -- 2.39.5