From 7b920016cb84462d52756682d89337175a3bfb2f Mon Sep 17 00:00:00 2001 From: maier Date: Fri, 21 Feb 2014 10:59:28 +0000 Subject: [PATCH] Reorganize make_periodicity_constraints, introduce a rotation matrix and augment documentation git-svn-id: https://svn.dealii.org/branches/branch_periodic_bc@32524 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 7 + deal.II/include/deal.II/dofs/dof_tools.h | 164 +++++++++--------- deal.II/include/deal.II/grid/grid_tools.h | 49 +++++- deal.II/source/dofs/dof_tools_constraints.cc | 98 ++++++----- .../source/dofs/dof_tools_constraints.inst.in | 19 +- deal.II/source/grid/grid_tools.cc | 15 +- deal.II/source/grid/grid_tools.inst.in | 12 +- 7 files changed, 213 insertions(+), 151 deletions(-) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 986abc2d37..4cc794bbc9 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -135,6 +135,13 @@ inconvenience this causes.

Specific improvements

    +
  1. New: DoFTools::make_periodicity_constraints and + GridTools::collect_periodic_faces now take a rotation matrix as optional + argument that describes how to rotate vector valued components of an + FESystem. +
    + (Daniel Arndt, Matthias Maier, 2014/02/21) +
  2. Improved: DoFRenumbering::Cuthill_McKee can now run with distributed triangulations with the renumbering only done within each processor's subdomain. diff --git a/deal.II/include/deal.II/dofs/dof_tools.h b/deal.II/include/deal.II/dofs/dof_tools.h index 22a6a74efa..f695233853 100644 --- a/deal.II/include/deal.II/dofs/dof_tools.h +++ b/deal.II/include/deal.II/dofs/dof_tools.h @@ -746,10 +746,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 @@ -843,15 +843,17 @@ namespace DoFTools * and any combination of that... * @endcode * - * More information on the topic can be found in the - * @ref GlossFaceOrientation "glossary" article. + * The optional parameter @p rotation is a rotation matrix that is + * applied to all vector valued components of an FESystem on @p face_1 + * prior to periodically constraining them to the DoFs of @p face_2. * - * @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 void @@ -862,7 +864,38 @@ 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 &rotation = IdentityMatrix(FaceIterator::space_dimension)); + + + + /** + * 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 + void + make_periodicity_constraints + (const std::vector > + &periodic_faces, + dealii::ConstraintMatrix &constraint_matrix, + const ComponentMask &component_mask = ComponentMask()); @@ -871,8 +904,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 @@ -902,12 +934,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 */ @@ -922,35 +956,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 - 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 @@ -967,12 +972,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 void @@ -984,53 +991,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 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 void make_periodicity_constraints - (const std::vector > - &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; //@} diff --git a/deal.II/include/deal.II/grid/grid_tools.h b/deal.II/include/deal.II/grid/grid_tools.h index a439b3da7b..f7246afe1e 100644 --- a/deal.II/include/deal.II/grid/grid_tools.h +++ b/deal.II/include/deal.II/grid/grid_tools.h @@ -988,19 +988,42 @@ namespace GridTools const std::set &boundary_ids = std::set()); + /** - * 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 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 rotation matrix that describes how vector valued DoFs of the first face + * should be rotated prior to constraining to the DoFs of the second + * face, see make_periodicity_constraints. + */ + FullMatrix rotation; }; + /** * An orthogonal equality test for faces. * @@ -1110,6 +1133,11 @@ 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 rotation can be specified that + * describes how vector valued DoFs of the first boundary @p b_id1 should + * be rotated prior to constraining to the DoFs of the second boundary @p + * b_id2, see make_periodicity_constraints. + * * @note The created std::vector can be used in * DoFTools::make_periodicity_constraints and in * parallel::distributed::Triangulation::add_periodicity to enforce @@ -1120,17 +1148,18 @@ 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 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 > &matched_pairs, - const dealii::Tensor<1,CONTAINER::space_dimension> &offset = dealii::Tensor<1,CONTAINER::space_dimension>()); + const Tensor<1,CONTAINER::space_dimension> &offset = dealii::Tensor<1,CONTAINER::space_dimension>(), + const FullMatrix &rotation = IdentityMatrix(CONTAINER::space_dimension)); /** @@ -1151,16 +1180,18 @@ 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 void collect_periodic_faces - (const CONTAINER &container, + (const CONTAINER &container, const types::boundary_id b_id, const int direction, std::vector > &matched_pairs, - const dealii::Tensor<1,CONTAINER::space_dimension> &offset = dealii::Tensor<1,CONTAINER::space_dimension>()); + const dealii::Tensor<1,CONTAINER::space_dimension> &offset = dealii::Tensor<1,CONTAINER::space_dimension>(), + const FullMatrix &rotation = IdentityMatrix(CONTAINER::space_dimension)); + /** diff --git a/deal.II/source/dofs/dof_tools_constraints.cc b/deal.II/source/dofs/dof_tools_constraints.cc index 890f2ac036..132396d24b 100644 --- a/deal.II/source/dofs/dof_tools_constraints.cc +++ b/deal.II/source/dofs/dof_tools_constraints.cc @@ -1875,6 +1875,9 @@ namespace DoFTools } + // Implementation of the low level interface: + + template void make_periodicity_constraints (const FaceIterator &face_1, @@ -1883,7 +1886,8 @@ namespace DoFTools const ComponentMask &component_mask, const bool face_orientation, const bool face_flip, - const bool face_rotation) + const bool face_rotation, + const FullMatrix &rotation) { static const int dim = FaceIterator::AccessorType::dimension; @@ -1969,7 +1973,8 @@ namespace DoFTools component_mask, face_orientation, face_flip, - face_rotation); + face_rotation, + rotation); } } else @@ -1978,13 +1983,13 @@ namespace DoFTools { if (face_2->has_children() == false) set_periodicity_constraints(face_2, face_1, - FullMatrix(IdentityMatrix(face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face)), + FullMatrix(IdentityMatrix(face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face)), // TODO constraint_matrix, component_mask, face_orientation, face_flip, face_rotation); else set_periodicity_constraints(face_1, face_2, - FullMatrix(IdentityMatrix(face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face)), + FullMatrix(IdentityMatrix(face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face)), // TODO constraint_matrix, component_mask, face_orientation, face_flip, face_rotation); @@ -1993,6 +1998,50 @@ namespace DoFTools + template + void + make_periodicity_constraints + (const std::vector > + &periodic_faces, + dealii::ConstraintMatrix &constraint_matrix, + const ComponentMask &component_mask) + { + typedef std::vector > + 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->rotation); + } + } + + + // High level interface variants: + + template void make_periodicity_constraints (const DH &dof_handler, @@ -2095,47 +2144,6 @@ namespace DoFTools - template - void - make_periodicity_constraints - (const std::vector > - &periodic_faces, - dealii::ConstraintMatrix &constraint_matrix, - const ComponentMask &component_mask) - { - typedef std::vector > - 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 diff --git a/deal.II/source/dofs/dof_tools_constraints.inst.in b/deal.II/source/dofs/dof_tools_constraints.inst.in index d61208de26..450332f5d0 100644 --- a/deal.II/source/dofs/dof_tools_constraints.inst.in +++ b/deal.II/source/dofs/dof_tools_constraints.inst.in @@ -31,7 +31,16 @@ for (DH : DOFHANDLERS; deal_II_dimension : DIMENSIONS) const DH::face_iterator &, dealii::ConstraintMatrix &, const ComponentMask &, - bool, bool, bool); + bool, bool, bool, + const FullMatrix &); + + template + void + DoFTools::make_periodicity_constraints + (const std::vector > &, + dealii::ConstraintMatrix &, + const ComponentMask &); + template void @@ -68,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 - (const std::vector > &, - dealii::ConstraintMatrix &, - const ComponentMask &); - #endif } diff --git a/deal.II/source/grid/grid_tools.cc b/deal.II/source/grid/grid_tools.cc index e96dab3f2f..c5ada5eac1 100644 --- a/deal.II/source/grid/grid_tools.cc +++ b/deal.II/source/grid/grid_tools.cc @@ -2381,7 +2381,8 @@ next_cell: std::set::type, unsigned int> > &pairs2, const int direction, std::vector > &matched_pairs, - const dealii::Tensor<1,CellIterator::AccessorType::space_dimension> &offset) + const dealii::Tensor<1,CellIterator::AccessorType::space_dimension> &offset, + const FullMatrix &rotation) { static const int space_dim = CellIterator::AccessorType::space_dimension; Assert (0<=direction && direction matched_face - = {{cell1, cell2},{face_idx1, face_idx2}, orientation}; + = {{cell1, cell2},{face_idx1, face_idx2}, orientation, rotation}; matched_pairs.push_back(matched_face); pairs2.erase(it2); ++n_matches; @@ -2437,7 +2438,8 @@ next_cell: const types::boundary_id b_id2, const int direction, std::vector > &matched_pairs, - const dealii::Tensor<1,CONTAINER::space_dimension> &offset) + const Tensor<1,CONTAINER::space_dimension> &offset, + const FullMatrix &rotation) { static const int dim = CONTAINER::dimension; static const int space_dim = CONTAINER::space_dimension; @@ -2476,7 +2478,7 @@ 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, rotation); } @@ -2488,7 +2490,8 @@ next_cell: const types::boundary_id b_id, const int direction, std::vector > &matched_pairs, - const dealii::Tensor<1,CONTAINER::space_dimension> &offset) + const Tensor<1,CONTAINER::space_dimension> &offset, + const FullMatrix &rotation) { static const int dim = CONTAINER::dimension; static const int space_dim = CONTAINER::space_dimension; @@ -2534,7 +2537,7 @@ 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, rotation); #ifdef DEBUG //check for standard orientation diff --git a/deal.II/source/grid/grid_tools.inst.in b/deal.II/source/grid/grid_tools.inst.in index a85fc66fb5..8586f204b9 100644 --- a/deal.II/source/grid/grid_tools.inst.in +++ b/deal.II/source/grid/grid_tools.inst.in @@ -320,14 +320,16 @@ for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II const types::boundary_id, const int, std::vector > &, - const Tensor<1,X::space_dimension> &); + const Tensor<1,X::space_dimension> &, + const FullMatrix &); template void collect_periodic_faces (const X &, const types::boundary_id, const int, std::vector > &, - const Tensor<1,X::space_dimension> &); + const Tensor<1,X::space_dimension> &, + const FullMatrix &); #endif @@ -350,7 +352,8 @@ for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS const types::boundary_id, const int, std::vector::cell_iterator> > &, - const Tensor<1,parallel::distributed::Triangulation::space_dimension> &); + const Tensor<1,parallel::distributed::Triangulation::space_dimension> &, + const FullMatrix &); template void @@ -359,7 +362,8 @@ for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS const types::boundary_id, const int, std::vector::cell_iterator> > &, - const Tensor<1,parallel::distributed::Triangulation::space_dimension> &); + const Tensor<1,parallel::distributed::Triangulation::space_dimension> &, + const FullMatrix &); \} #endif -- 2.39.5