From: Daniel Arndt Date: Mon, 16 Sep 2013 17:02:05 +0000 (+0000) Subject: Changed the data type for periodic face pairs X-Git-Tag: v8.1.0~805 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d5fbf81519f3626d9730ecce88f9e086e4a7f69a;p=dealii.git Changed the data type for periodic face pairs git-svn-id: https://svn.dealii.org/trunk@30732 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/distributed/tria.h b/deal.II/include/deal.II/distributed/tria.h index 95c34afc1d..5d647258f1 100644 --- a/deal.II/include/deal.II/distributed/tria.h +++ b/deal.II/include/deal.II/distributed/tria.h @@ -158,7 +158,8 @@ namespace internal } } - +//forward declaration of the data type for periodic face pairs +namespace GridTools {template struct PeriodicFacePair;} namespace parallel { @@ -700,15 +701,36 @@ namespace parallel get_p4est_tree_to_coarse_cell_permutation() const; + /** * Join faces in the p4est forest due to periodic boundary conditions. * + * The vector can be filled by the function + * GridTools::collect_periodic_faces. + * + * @todo At the moment just default orientation is implemented. + * + * @note Before this function can be used the triangulation has to be + * initialized and must not be refined. + * Calling this function more than once is possible, but not recommended: + * The function destroys and rebuilds the p4est forest each time it is called. + */ + void + add_periodicity + (const std::vector >&); + + /** + * Same as the function above, but takes a different argument. + * * The entries in the std::vector should have the form * std_cxx1x::tuple. * * The vector can be filled by the function - * DoFTools::identify_periodic_face_pairs. - * + * GridTools::identify_periodic_face_pairs. + * + * @note This function can only be used if the faces are in + * default orientation. + * * @note Before this function can be used the triangulation has to be * initialized and must not be refined. * Calling this function more than once is possible, but not recommended: @@ -720,6 +742,7 @@ namespace parallel cell_iterator, unsigned int> >&); + private: /** * MPI communicator to be diff --git a/deal.II/include/deal.II/dofs/dof_tools.h b/deal.II/include/deal.II/dofs/dof_tools.h index e4d237faf3..9b39197104 100644 --- a/deal.II/include/deal.II/dofs/dof_tools.h +++ b/deal.II/include/deal.II/dofs/dof_tools.h @@ -53,6 +53,7 @@ class ConstraintMatrix; template class InterGridMap; template class Mapping; +namespace GridTools {template struct PeriodicFacePair;} //TODO: map_support_points_to_dofs should generate a multimap, rather than just a map, since several dofs may be located at the same support point @@ -1087,13 +1088,15 @@ 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 = ComponentMask(), - const bool face_orientation = true, - const bool face_flip = false, - const bool face_rotation = false); + make_periodicity_constraints + (const FaceIterator &face_1, + const typename identity::type &face_2, + dealii::ConstraintMatrix &constraint_matrix, + const ComponentMask &component_mask = ComponentMask(), + const bool face_orientation = true, + const bool face_flip = false, + const bool face_rotation = false); + /** @@ -1143,12 +1146,13 @@ namespace DoFTools */ 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::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()); /** @@ -1167,17 +1171,18 @@ namespace DoFTools * * @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators" * - * @author Matthias Maier, 2012 + * @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()); + 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()); /** @@ -1204,11 +1209,12 @@ namespace DoFTools */ template void - make_periodicity_constraints (const DH &dof_handler, - const types::boundary_id b_id, - const int direction, - dealii::ConstraintMatrix &constraint_matrix, - const ComponentMask &component_mask = ComponentMask()); + make_periodicity_constraints + (const DH &dof_handler, + const types::boundary_id b_id, + const int direction, + dealii::ConstraintMatrix &constraint_matrix, + const ComponentMask &component_mask = ComponentMask()); /** @@ -1233,12 +1239,31 @@ namespace DoFTools */ template void - make_periodicity_constraints (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()); + make_periodicity_constraints + (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()); + + /** + * Same as above but the periodicity information is given by + * @p periodic_faces. This std::vector 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. + */ + template + void + make_periodicity_constraints + (const std::vector > + &periodic_faces, + dealii::ConstraintMatrix &constraint_matrix, + const ComponentMask &component_mask = ComponentMask()); //@} diff --git a/deal.II/include/deal.II/grid/grid_tools.h b/deal.II/include/deal.II/grid/grid_tools.h index 83e40cf654..fc43a3cd46 100644 --- a/deal.II/include/deal.II/grid/grid_tools.h +++ b/deal.II/include/deal.II/grid/grid_tools.h @@ -990,7 +990,18 @@ 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 + */ + template + struct PeriodicFacePair + { + CellIterator cell[2]; + unsigned int face_idx[2]; + std::bitset<3> orientation; + }; /** * An orthogonal equality test for faces. @@ -1076,8 +1087,8 @@ namespace GridTools /** - * This function loops over all faces from @p begin to @p end and - * collects a set of periodic face pairs for @p direction: + * This function will collect periodic face pairs on the highest (i.e. + * coarsest) mesh level. * * Define a 'first' boundary as all boundary faces having boundary_id * @p b_id1 and a 'second' boundary consisting of all faces belonging @@ -1087,45 +1098,33 @@ namespace GridTools * boundary with faces belonging to the second boundary with the help * of orthogonal_equality(). * - * The bitset that is returned together with the second face encodes the + * The bitset that is returned inside of PeriodicFacePair encodes the * _relative_ orientation of the first face with respect to the second * face, see the documentation of orthogonal_equality for further details. * + * The @p direction refers to the space direction in which periodicity + * is enforced. + * * 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. This can * be used to implement conditions such as $u(0,y)=u(1,y+1)$. * - * @author Matthias Maier, 2012 - */ - template - std::map > > - collect_periodic_face_pairs (const FaceIterator &begin, - const typename identity::type &end, - const types::boundary_id b_id1, - const types::boundary_id b_id2, - const int direction, - const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset); - - - /** - * Same function as above, but accepts a Triangulation or DoFHandler - * object @p container (a container is a collection of objects, here a - * collection of cells) instead of an explicit face iterator range. - * - * This function will collect periodic face pairs on the highest (i.e. - * coarsest) mesh level. + * @note The created std::vector can be used in + * DoFTools::make_periodicity_constraints and in + * parallel::distributed::Triangulation::add_periodicity to enforce + * periodicity algebraically. * - * @author Matthias Maier, 2012 + * @author Daniel Arndt, 2013 */ template - std::map > > - collect_periodic_face_pairs (const DH &container, - const types::boundary_id b_id1, - const types::boundary_id b_id2, - const int direction, - const dealii::Tensor<1,DH::space_dimension> &offset); - + std::vector > + collect_periodic_faces + (const DH &dof_handler, + const types::boundary_id b_id1, + const types::boundary_id b_id2, + const unsigned int direction, + const dealii::Tensor<1,DH::space_dimension> &offset); /** * This compatibility version of collect_periodic_face_pairs only works @@ -1145,40 +1144,101 @@ namespace GridTools * meshes with cells not in @ref GlossFaceOrientation * "standard orientation". * + * @author Daniel Arndt, 2013 + */ + template + std::vector > + collect_periodic_faces + (const DH &dof_handler, + const types::boundary_id b_id, + const unsigned int direction, + const dealii::Tensor<1,DH::space_dimension> &offset); + + /** + * This function does the same as collect_periodic_faces but returns a + * different data type. + * + * @author Matthias Maier, 2012 + * + * @note The returned data type is not compatible with + * DoFTools::make_periodicity_constraints + * + * @deprecated + */ + template + std::map > > + collect_periodic_face_pairs + (const DH &container, + const types::boundary_id b_id1, + const types::boundary_id b_id2, + int direction, + const dealii::Tensor<1,DH::space_dimension> &offset) DEAL_II_DEPRECATED; + + /** + * This compatibility version of collect_periodic_face_pairs only works + * on grids with cells in @ref GlossFaceOrientation "standard orientation". + * * @author Matthias Maier, 2012 + * + * @note The returned data type is not compatible with + * DoFTools::make_periodicity_constraints + * + * @deprecated */ template std::map - collect_periodic_face_pairs (const DH &dof_handler, /*TODO: Name*/ - const types::boundary_id b_id, - const int direction, - const dealii::Tensor<1,DH::space_dimension> &offset); + collect_periodic_face_pairs + (const DH &dof_handler, /*TODO: Name*/ + const types::boundary_id b_id, + int direction, + const dealii::Tensor<1,DH::space_dimension> &offset) DEAL_II_DEPRECATED; /** - * Add periodicity information to the @p periodicity_vector for the - * boundaries with boundary id @p b_id1 and @p b_id2 in cartesian - * direction @p direction. + * This version loops over all faces from @p begin to @p end + * instead of accepting a DoFHandler or a Triangulation. * - * This function tries to match all faces belonging to the first - * boundary with faces belonging to the second boundary - * by comparing the center of the cell faces. To find the correct - * corresponding faces, the direction argument indicates in which - * cartesian direction periodicity should be set. - * ((0,1,2) -> (x,y,z)-direction) + * @author Matthias Maier, 2012 + * + * @note This function cannot produce the return as the other + * collect_periodic_faces functions. + * + * @deprecated + */ + template + std::map > > + collect_periodic_face_pairs + (const FaceIterator &begin, + const typename identity::type &end, + const types::boundary_id b_id1, + const types::boundary_id b_id2, + const int direction, + const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset) + DEAL_II_DEPRECATED; + + /** + * This function does the same as collect_periodic_faces but returns a + * different data type. + * + * @author Daniel Arndt, 2013 * - * The output of this function can be used in + * @note The returned data type is not compatible with + * DoFTools::make_periodicity_constraints, but with a version of * parallel::distributed::Triangulation::add_periodicity + * + * @note Use collect_periodic_faces instead. + * + * @deprecated */ template void identify_periodic_face_pairs - (const DH &dof_handler, - const types::boundary_id b_id1, - const types::boundary_id b_id2, - const unsigned int direction, - std::vector > - &periodicity_vector); + (const DH &dof_handler, + const types::boundary_id b_id1, + const types::boundary_id b_id2, + const unsigned int direction, + std::vector > + &periodicity_vector) DEAL_II_DEPRECATED; /** diff --git a/deal.II/source/distributed/tria.cc b/deal.II/source/distributed/tria.cc index ce4cbc6e97..b036739d6f 100644 --- a/deal.II/source/distributed/tria.cc +++ b/deal.II/source/distributed/tria.cc @@ -23,8 +23,8 @@ #include #include #include -#include #include +#include #ifdef DEAL_II_WITH_P4EST # include @@ -3436,39 +3436,30 @@ namespace parallel return mpi_communicator; } - - template + template void Triangulation::add_periodicity - (const std::vector >& - periodicity_vector) + (const std::vector >& + periodicity_vector) { -#if DEAL_II_P4EST_VERSION_GTE(0,3,4,1) + #if DEAL_II_P4EST_VERSION_GTE(0,3,4,1) Assert (triangulation_has_content == true, ExcMessage ("The triangulation is empty!")); Assert (this->n_levels() == 1, ExcMessage ("The triangulation is refined!")); - - typename std::vector - >::const_iterator periodic_tuple; - periodic_tuple = periodicity_vector.begin(); - - typename std::vector - >::const_iterator periodic_end; - periodic_end = periodicity_vector.end(); - - for (; periodic_tuple > + FaceVector; + typename FaceVector::const_iterator it, periodic_end; + it = periodicity_vector.begin(); + + for (; it(*periodic_tuple); - const cell_iterator second_cell=std_cxx1x::get<2>(*periodic_tuple); - const unsigned int face_right=std_cxx1x::get<3>(*periodic_tuple); - const unsigned int face_left=std_cxx1x::get<1>(*periodic_tuple); - + const cell_iterator first_cell = it->cell[0]; + const cell_iterator second_cell = it->cell[1]; + const unsigned int face_left = it->face_idx[0]; + const unsigned int face_right = it->face_idx[1]; + //respective cells of the matching faces in p4est const unsigned int tree_left = coarse_cell_to_p4est_tree_permutation[std::distance(this->begin(), @@ -3477,6 +3468,12 @@ namespace parallel = coarse_cell_to_p4est_tree_permutation[std::distance(this->begin(), second_cell)]; + //TODO Add support for non default orientation. + Assert(it->orientation == 1, + ExcMessage("Found a face match with non standard orientation. " + "This function is only suitable for meshes with " + "cells in default orientation")); + dealii::internal::p4est::functions:: connectivity_join_faces (connectivity, tree_left, @@ -3484,6 +3481,7 @@ namespace parallel face_left, face_right, /* orientation */ 0); + /* The orientation parameter above describes the difference between * the cell on the left and the cell on the right would number of the * corners of the face. In the periodic domains most users will want, @@ -3495,25 +3493,60 @@ namespace parallel * date. */ } - - + + Assert(dealii::internal::p4est::functions::connectivity_is_valid - (connectivity) == 1, - ExcInternalError()); - - // now create a forest out of the connectivity data structure + (connectivity) == 1, ExcInternalError()); + + // now create a forest out of the connectivity data structure dealii::internal::p4est::functions::destroy (parallel_forest); parallel_forest = dealii::internal::p4est::functions:: new_forest (mpi_communicator, connectivity, - /* minimum initial number of quadrants per tree */ 0, - /* minimum level of upfront refinement */ 0, - /* use uniform upfront refinement */ 1, - /* user_data_size = */ 0, - /* user_data_constructor = */ NULL, - /* user_pointer */ this); + /* minimum initial number of quadrants per tree */ 0, + /* minimum level of upfront refinement */ 0, + /* use uniform upfront refinement */ 1, + /* user_data_size = */ 0, + /* user_data_constructor = */ NULL, + /* user_pointer */ this); + + #else + Assert(false, ExcMessage ("Need p4est version >= 0.3.4.1!")); + #endif + } + + + template + void + Triangulation::add_periodicity + (const std::vector >& + periodicity_vector) + { +#if DEAL_II_P4EST_VERSION_GTE(0,3,4,1) + typedef std::vector > + FaceVector; + + typename FaceVector::const_iterator it, end_periodic; + it = periodicity_vector.begin(); + end_periodic = periodicity_vector.end(); + + std::vector > periodic_faces; + + for(; it!=end_periodic; ++it) + { + const cell_iterator cell1 = std::get<0> (*it); + const cell_iterator cell2 = std::get<2> (*it); + const unsigned int face_idx1 = std::get<1> (*it); + const unsigned int face_idx2 = std::get<3> (*it); + const GridTools::PeriodicFacePair matched_face + = {{cell1, cell2},{face_idx1, face_idx2}, 1}; + periodic_faces.push_back(matched_face); + } + add_periodicity(periodic_faces); #else Assert(false, ExcMessage ("Need p4est version >= 0.3.4.1!")); #endif diff --git a/deal.II/source/dofs/dof_tools_constraints.cc b/deal.II/source/dofs/dof_tools_constraints.cc index 402557bee3..adaf93b7f0 100644 --- a/deal.II/source/dofs/dof_tools_constraints.cc +++ b/deal.II/source/dofs/dof_tools_constraints.cc @@ -1970,7 +1970,7 @@ namespace DoFTools } - + template void make_periodicity_constraints (const DH &dof_handler, @@ -2010,24 +2010,23 @@ namespace DoFTools ExcMessage ("The boundary indicators b_id1 and b_id2 must be" "different to denote different boundaries.")); - typedef typename DH::face_iterator FaceIterator; - typedef std::map > > FaceMap; + typedef std::vector > FaceVector; // Collect matching periodic cells on the coarsest level: - FaceMap matched_cells = - GridTools::collect_periodic_face_pairs(dof_handler, - b_id1, b_id2, - direction, offset); + FaceVector matched_cells = + GridTools::collect_periodic_faces(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(); + for (typename FaceVector::iterator it = matched_cells.begin(); it != matched_cells.end(); ++it) { typedef typename DH::face_iterator FaceIterator; - const FaceIterator &face_1 = it->first; - const FaceIterator &face_2 = it->second.first; - const std::bitset<3> &orientation = it->second.second; + const FaceIterator &face_1 = it->cell[0]->face(it->face_idx[0]); + const FaceIterator &face_2 = it->cell[1]->face(it->face_idx[1]); + const std::bitset<3> &orientation = it->orientation; Assert(face_1->at_boundary() && face_2->at_boundary(), ExcInternalError()); @@ -2088,23 +2087,22 @@ namespace DoFTools Assert(dim == space_dim, ExcNotImplemented()); - typedef typename DH::face_iterator FaceIterator; - typedef std::map FaceMap; + typedef std::vector > FaceVector; // Collect matching periodic cells on the coarsest level: - FaceMap matched_cells = - GridTools::collect_periodic_face_pairs(dof_handler, - b_id, - direction, offset); + const FaceVector matched_cells = + GridTools::collect_periodic_faces(dof_handler, b_id, + direction, offset); // And apply the low level make_periodicity_constraints function to // every matching pair: - for (typename FaceMap::iterator it = matched_cells.begin(); + for (typename FaceVector::const_iterator it = matched_cells.begin(); it != matched_cells.end(); ++it) { typedef typename DH::face_iterator FaceIterator; - const FaceIterator &face_1 = it->first; - const FaceIterator &face_2 = it->second; + 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()); @@ -2126,6 +2124,38 @@ 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(); + + 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]); + + make_periodicity_constraints(face_1, + face_2, + constraint_matrix, + component_mask, + it->orientation[0], + it->orientation[1], + it->orientation[2]); + } + } + + + namespace internal { namespace diff --git a/deal.II/source/dofs/dof_tools_constraints.inst.in b/deal.II/source/dofs/dof_tools_constraints.inst.in index 246f53aa02..5c9bdd7c88 100644 --- a/deal.II/source/dofs/dof_tools_constraints.inst.in +++ b/deal.II/source/dofs/dof_tools_constraints.inst.in @@ -68,6 +68,14 @@ 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 eeeb0dc6ad..292b113d62 100644 --- a/deal.II/source/grid/grid_tools.cc +++ b/deal.II/source/grid/grid_tools.cc @@ -2217,7 +2217,6 @@ next_cell: } - /* * Internally used in orthogonal_equality * @@ -2226,7 +2225,7 @@ next_cell: * * See the comment on the next function as well as the detailed * documentation of make_periodicity_constraints and - * collect_periodic_face_pairs for details + * collect_periodic_faces for details */ template struct OrientationLookupTable {}; @@ -2355,14 +2354,82 @@ next_cell: /* - * Internally used in collect_periodic_face_pairs + * Internally used in collect_periodic_faces + */ + template + std::vector > + match_periodic_face_pairs + (std::set > + &pairs1, + std::set::type, unsigned int> > + &pairs2, + int direction, + const dealii::Tensor<1,CellIterator::AccessorType::space_dimension> + &offset) + { + static const int space_dim = CellIterator::AccessorType::space_dimension; + Assert (0<=direction && direction > matched_faces; + + // Match with a complexity of O(n^2). This could be improved... + std::bitset<3> orientation; + typedef typename std::set + >::const_iterator PairIterator; + for (PairIterator it1 = pairs1.begin(); it1 != pairs1.end(); ++it1) + { + for (PairIterator it2 = pairs2.begin(); it2 != pairs2.end(); ++it2) + { + const CellIterator cell1 = it1->first; + const CellIterator cell2 = it2->first; + const unsigned int face_idx1 = it1->second; + const unsigned int face_idx2 = it2->second; + if (GridTools::orthogonal_equality(orientation, + cell1->face(face_idx1), + cell2->face(face_idx2), + direction, offset)) + { + // We have a match, so insert the matching pairs and + // remove the matched cell in pairs2 to speed up the + // matching: + const PeriodicFacePair matched_face + = {{cell1, cell2},{face_idx1, face_idx2}, orientation}; + matched_faces.push_back(matched_face); + pairs2.erase(it2); + break; + } + } + } + + AssertThrow (matched_faces.size() == pairs1.size() && pairs2.size() == 0, + ExcMessage ("Unmatched faces on periodic boundaries")); + + return matched_faces; + } + + /* Deprecated version of the function above with different return value. + * It is used the deprecated collect_periodic_face_pairs. */ template std::map > > - match_periodic_face_pairs (std::set &faces1, /* not const! */ - std::set::type> &faces2, /* not const! */ - int direction, - const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset) + match_periodic_face_pairs + (std::set &faces1, /* not const! */ + std::set::type> &faces2, /* not const! */ + int direction, + const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset) + DEAL_II_DEPRECATED; + + template + std::map > > + match_periodic_face_pairs + (std::set &faces1, /* not const! */ + std::set::type> &faces2, /* not const! */ + int direction, + const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset) { static const int space_dim = FaceIterator::AccessorType::space_dimension; Assert (0<=direction && direction + std::vector > + collect_periodic_faces + (const DH &dof_handler, + const types::boundary_id b_id1, + const types::boundary_id b_id2, + const unsigned int direction, + const dealii::Tensor<1,DH::space_dimension> &offset) + { + static const unsigned int dim = DH::dimension; + static const unsigned int space_dim = DH::space_dimension; + Assert (0<=direction && direction > pairs1; + std::set > pairs2; + + for (typename DH::cell_iterator cell = dof_handler.begin(0); + cell != dof_handler.end(0); ++cell) + { + for (unsigned int i = 0; i < GeometryInfo::faces_per_cell; ++i) + { + const typename DH::face_iterator face = cell->face(i); + if (face->at_boundary() && face->boundary_indicator() == b_id1) + { + const std::pair pair1 + = std::make_pair(cell, i); + pairs1.insert(pair1); + } + + if (face->at_boundary() && face->boundary_indicator() == b_id2) + { + const std::pair pair2 + = std::make_pair(cell, i); + pairs2.insert(pair2); + } + } + } + + Assert (pairs1.size() == pairs2.size(), + ExcMessage ("Unmatched faces on periodic boundaries")); + + // and call match_periodic_face_pairs that does the actual matching: + return match_periodic_face_pairs(pairs1, pairs2, direction, offset); + } + + template + std::vector > + collect_periodic_faces (const DH &dof_handler, + const types::boundary_id b_id, + const unsigned int direction, + const dealii::Tensor<1,DH::space_dimension> &offset) + { + static const unsigned int dim = DH::dimension; + static const unsigned int space_dim = DH::space_dimension; + Assert (0<=direction && direction > pairs1; + std::set > pairs2; + + for (typename DH::cell_iterator cell = dof_handler.begin(0); + cell != dof_handler.end(0); ++cell) + { + const typename DH::face_iterator face_1 = cell->face(2*direction); + const typename DH::face_iterator face_2 = cell->face(2*direction+1); + + if (face_1->at_boundary() && face_1->boundary_indicator() == b_id) + { + const std::pair pair1 + = std::make_pair(cell, 2*direction); + pairs1.insert(pair1); + } + + if (face_2->at_boundary() && face_2->boundary_indicator() == b_id) + { + const std::pair pair1 + = std::make_pair(cell, 2*direction+1); + pairs1.insert(pair1); + } + } + + Assert (pairs1.size() == pairs2.size(), + ExcMessage ("Unmatched faces on periodic boundaries")); + + // and call match_periodic_face_pairs that does the actual matching: + + typedef std::vector > + FaceVector; + + FaceVector matching = match_periodic_face_pairs(pairs1, pairs2, + direction, offset); + + for (typename FaceVector::iterator pairing = matching.begin(); + pairing != matching.end(); ++pairing) + { + Assert(pairing->orientation == 1, + ExcMessage("Found a face match with non standard orientation. " + "This function is only suitable for meshes with cells " + "in default orientation")); + } + + return matching; + } template std::map > > @@ -2436,8 +2616,6 @@ next_cell: return match_periodic_face_pairs(faces1, faces2, direction, offset); } - - template std::map > > collect_periodic_face_pairs (const DH &dof_handler, @@ -2446,41 +2624,26 @@ next_cell: int direction, const dealii::Tensor<1,DH::space_dimension> &offset) { - static const int dim = DH::dimension; - static const int space_dim = DH::space_dimension; - Assert (0<=direction && direction > FaceVector; - std::set faces1; - std::set faces2; + const FaceVector face_vector + = collect_periodic_faces (dof_handler, b_id1, b_id2, direction, offset); - for (typename DH::cell_iterator cell = dof_handler.begin(0); - cell != dof_handler.end(0); ++cell) - { - for (unsigned int i = 0; i < GeometryInfo::faces_per_cell; ++i) - { - const typename DH::face_iterator face = cell->face(i); - - if (face->at_boundary() && face->boundary_indicator() == b_id1) - faces1.insert(face); - - if (face->at_boundary() && face->boundary_indicator() == b_id2) - faces2.insert(face); - } - } - - Assert (faces1.size() == faces2.size(), - ExcMessage ("Unmatched faces on periodic boundaries")); + std::map > > + return_value; + for(typename FaceVector::const_iterator it = face_vector.begin(); + it != face_vector.end(); ++it) + { + const typename DH::face_iterator face1 = it->cell[0]->face(it->face_idx[0]); + const typename DH::face_iterator face2 = it->cell[1]->face(it->face_idx[1]); + return_value[face1] = std::make_pair(face2, it->orientation); + } - // and call match_periodic_face_pairs that does the actual matching: - return match_periodic_face_pairs(faces1, faces2, direction, offset); + return return_value; } - - + template std::map collect_periodic_face_pairs (const DH &dof_handler, @@ -2488,59 +2651,24 @@ next_cell: int direction, const dealii::Tensor<1,DH::space_dimension> &offset) { - static const int dim = DH::dimension; - static const int space_dim = DH::space_dimension; - Assert (0<=direction && direction faces1; - std::set faces2; - - for (typename DH::cell_iterator cell = dof_handler.begin(0); - cell != dof_handler.end(0); ++cell) - { - const typename DH::face_iterator face_1 = cell->face(2*direction); - const typename DH::face_iterator face_2 = cell->face(2*direction+1); - - if (face_1->at_boundary() && face_1->boundary_indicator() == b_id) - faces1.insert(face_1); - - if (face_2->at_boundary() && face_2->boundary_indicator() == b_id) - faces2.insert(face_2); - } - - Assert (faces1.size() == faces2.size(), - ExcMessage ("Unmatched faces on periodic boundaries")); - - // and call match_periodic_face_pairs that does the actual matching: - - typedef std::map > > FaceMap; - FaceMap matching = match_periodic_face_pairs(faces1, faces2, direction, offset); - - std::map - return_value; - - for (typename FaceMap::iterator pairing = matching.begin(); - pairing != matching.end(); ++pairing) - { - Assert(pairing->second.second == 1, - ExcMessage("Found a face match with non standard orientation. " - "This function is only suitable for meshes with cells " - "in default orientation")); - - return_value[pairing->first] = pairing->second.first; - } - + typedef std::vector > FaceVector; + + const FaceVector face_vector + = collect_periodic_faces (dof_handler, b_id, direction, offset); + + std::map return_value; + for(typename FaceVector::const_iterator it = face_vector.begin(); + it != face_vector.end(); ++it) + { + const typename DH::face_iterator face1 = it->cell[0]->face(it->face_idx[0]); + const typename DH::face_iterator face2 = it->cell[1]->face(it->face_idx[1]); + return_value[face1] = face2; + } + return return_value; } - + template void identify_periodic_face_pairs @@ -2552,72 +2680,26 @@ next_cell: typename DH::cell_iterator, unsigned int> > &periodicity_vector) { - static const unsigned int dim = DH::dimension; - static const unsigned int spacedim = DH::space_dimension; - - Assert (0<=direction && direction, - Point > face_locations; - - // Collect faces with boundary_indicator b_id1 - typename DH::cell_iterator cell = dof_handler.begin(); - typename DH::cell_iterator endc = dof_handler.end(); - for (; cell != endc; ++cell) - for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) - if(cell->face(f)->boundary_indicator() == b_id1) - { - Point face_center (cell->face(f)->center()); - face_center(direction)=0.; - const std::pair - cell_face_pair = std::make_pair(cell, f); - face_locations[cell_face_pair]=face_center; - } - - // Match faces with boundary_indicator b_id2 to the ones in - //face_locations - cell = dof_handler.begin(); - for (; cell != endc; ++cell) - for(unsigned int f=0;f::faces_per_cell;++f) - if(cell->face(f)->boundary_indicator() == b_id2) - { - typename std::map, - Point >::iterator p - = face_locations.begin(); - - Point center2 (cell->face(f)->center()); - center2(direction)=0.; - - for (; p != face_locations.end(); ++p) - { - if (center2.distance(p->second) < 1e-4*cell->face(f)->diameter()) - { - const std_cxx1x::tuple - periodic_tuple (p->first.first, - p->first.second, - cell, - f); - - periodicity_vector.push_back(periodic_tuple); - - face_locations.erase(p); - break; - } - Assert (p != face_locations.end(), - ExcMessage ("No corresponding face was found!")); - } - } + typedef std::vector > + FaceVector; + const FaceVector periodic_faces + = collect_periodic_faces(dof_handler, b_id1, b_id2, direction, + dealii::Tensor<1,DH::space_dimension> ()); + + typename FaceVector::const_iterator it, end_faces; + it = periodic_faces.begin(); + end_faces = periodic_faces.end(); + + for(; it!=end_faces; ++it) + { - Assert (face_locations.size() == 0, - ExcMessage ("There are unmatched faces!")); + const std_cxx1x::tuple + periodic_tuple (it->cell[0], it->face_idx[0], + it->cell[1], it->face_idx[1]); + + periodicity_vector.push_back(periodic_tuple); + } } diff --git a/deal.II/source/grid/grid_tools.inst.in b/deal.II/source/grid/grid_tools.inst.in index fdaecb93c5..465931c953 100644 --- a/deal.II/source/grid/grid_tools.inst.in +++ b/deal.II/source/grid/grid_tools.inst.in @@ -234,6 +234,23 @@ for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II #if deal_II_dimension >= 2 + template + std::vector > + collect_periodic_faces + (const X &, + const types::boundary_id, + const types::boundary_id, + unsigned int, + const Tensor<1,deal_II_space_dimension> &); + + template + std::vector > + collect_periodic_faces + (const X &, + const types::boundary_id, + unsigned int, + const Tensor<1,deal_II_space_dimension> &); + template std::map > > collect_periodic_face_pairs @@ -243,7 +260,6 @@ for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II const types::boundary_id, int, const Tensor<1,deal_II_space_dimension> &); - template std::map > > collect_periodic_face_pairs @@ -270,7 +286,7 @@ for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II const types::boundary_id, int, const Tensor<1,deal_II_space_dimension> &); - + template void identify_periodic_face_pairs