From 6f564f5b39ade2b2c2b0dab1008eaa5640d707fb Mon Sep 17 00:00:00 2001 From: maier Date: Sat, 1 Dec 2012 17:16:43 +0000 Subject: [PATCH] Bugfixes: - Refactor GridTools::collect_periodic_cell_pairs to GridTools::collect_periodic_face_pairs - Make GridTools::collect_periodix_face_pairs orientation (orientation, flip, rotation) aware git-svn-id: https://svn.dealii.org/trunk@27715 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/grid/grid_tools.h | 178 +++++++---- deal.II/source/grid/grid_tools.cc | 355 ++++++++++++++++------ deal.II/source/grid/grid_tools.inst.in | 67 +++- 3 files changed, 431 insertions(+), 169 deletions(-) diff --git a/deal.II/include/deal.II/grid/grid_tools.h b/deal.II/include/deal.II/grid/grid_tools.h index 3387e7f8fe..8be9ef13e9 100644 --- a/deal.II/include/deal.II/grid/grid_tools.h +++ b/deal.II/include/deal.II/grid/grid_tools.h @@ -20,6 +20,7 @@ #include #include +#include #include DEAL_II_NAMESPACE_OPEN @@ -278,7 +279,7 @@ namespace GridTools * function which computes all * active vertex patches by looping * over cells. - * + * * Find and return a vector of * iterators to active cells that * surround a given vertex @p vertex. @@ -914,71 +915,136 @@ namespace GridTools = std::set()); + /** - * This function loops over all cells - * from @p begin to @p end and collects a - * set of periodic cell pairs for - * @p direction: + * An orthogonal equality test for faces. + * + * @p face1 and @p face2 are considered equal, if a one to one matching + * between its vertices can be achieved via an orthogonal equality + * relation: Two vertices v_1 and v_2 are considered + * equal, if (v_1 + offset) - v_2 is parallel to the unit + * vector in @p direction. + * + * If the matching was successful, the _relative_ orientation of @p face1 + * with respect to @p face2 is returned in the bitset @p orientation, + * where + * @code + * orientation[0] -> face_orientation + * orientation[1] -> face_flip + * orientation[2] -> face_rotation + * @endcode + * + * In 2D face_orientation is always true, + * face_rotation is always false, and face_flip has the + * meaning of line_flip. More precisely in 3d: + * + * face_orientation: + * true if @p face1 and @p face2 have the same orientation. + * Otherwise, the vertex indices of @p face1 match the vertex indices of + * @p face2 in the following manner: + * + * @code + * face1: face2: + * + * 1 - 3 2 - 3 + * | | <--> | | + * 0 - 2 0 - 1 + * @endcode * - * 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. + * face_flip: + * true if the matched vertices are rotated by 180 degrees: * - * This function tries to match - * all cells with faces belonging to the - * left' boundary with the cells with faces - * belonging to the 'right' boundary - * with the help of an - * orthogonal equality relation: - * Two cells do match if the vertices of - * the respective boundary faces - * can be transformed into each other by - * parallel translation in @p direction. + * @code + * face1: face2: * - * 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)$. + * 1 - 0 2 - 3 + * | | <--> | | + * 3 - 2 0 - 1 + * @endcode + * + * face_rotation: true if the matched vertices are + * rotated by 90 degrees counterclockwise: + * + * @code + * face1: face2: + * + * 0 - 2 2 - 3 + * | | <--> | | + * 1 - 3 0 - 1 + * @endcode + * + * and any combination of that... + * More information on the topic can be found in the + * @ref GlossFaceOrientation "glossary" article. + * + * @author Matthias Maier, 2012 */ - template - std::map - collect_periodic_cell_pairs (const CellIterator &begin, - const typename identity::type &end, - const types::boundary_id boundary_component, + template + bool + orthogonal_equality (std::bitset<3> &orientation, + const FaceIterator &face1, + const FaceIterator &face2, + const int direction, + const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset); + + + /** + * Same function as above, but doesn't return the actual orientation + */ + template + bool + orthogonal_equality (const FaceIterator &face1, + const FaceIterator &face2, + const int direction, + const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset); + + + /** + * This function loops over all faces from @p begin to @p end and + * collects a set of periodic face pairs for @p direction: + * + * 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 orthogonal_equality. + * + * 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, int direction, - const dealii::Tensor<1,CellIterator::AccessorType::space_dimension> - &offset = dealii::Tensor<1,CellIterator::AccessorType::space_dimension>()); + const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset); /** - * Same as above but matches all active - * cells, i.e. this function calls the - * above function with - * @p dof_handler.begin_active() and - * @p dof_handler.end() as the first two - * arguments. + * Same function as above, but accepts a Triangulation or DoFHandler + * object @p dof_handler instead of an explicit face iterator range. + * + * This function will collect periodic face pairs on the highest (i.e. + * coarsest) mesh level. + * + * @author Matthias Maier, 2012 */ template - std::map - collect_periodic_cell_pairs (const DH &dof_handler, - const types::boundary_id boundary_component, - int direction, - const dealii::Tensor<1,DH::space_dimension> - &offset = Tensor<1,DH::space_dimension>()); + std::map > > + collect_periodic_face_pairs (const DH &dof_handler, /*TODO: Name*/ + const types::boundary_id b_id1, + const types::boundary_id b_id2, + int direction, + const dealii::Tensor<1,DH::space_dimension> &offset); + /** @@ -1030,7 +1096,7 @@ namespace GridTools << " is not used in the given triangulation"); -} +} /*namespace GridTools*/ diff --git a/deal.II/source/grid/grid_tools.cc b/deal.II/source/grid/grid_tools.cc index 40b99fc361..a01fff9f02 100644 --- a/deal.II/source/grid/grid_tools.cc +++ b/deal.II/source/grid/grid_tools.cc @@ -12,6 +12,7 @@ //--------------------------------------------------------------------------- +#include #include #include #include @@ -2159,142 +2160,294 @@ next_cell: return surface_to_volume_mapping; } - // Internally used in - // collect_periodic_cell_pairs. - template - inline bool orthogonal_equality (const CellIterator &cell1, - const CellIterator &cell2, - const int direction, - const dealii::Tensor<1,CellIterator::AccessorType::space_dimension> - &offset) + + + /* + * Internally used in orthogonal_equality + * + * An orthogonal equality test for points: + * + * point1 and point2 are considered equal, if + * (point1 + offset) - point2 + * is parallel to the unit vector in + */ + template + inline bool orthogonal_equality (const dealii::Point &point1, + const dealii::Point &point2, + const int direction, + const dealii::Tensor<1,spacedim> &offset) { - // An orthogonal equality test for - // cells: - // Two cells are equal if the - // corresponding vertices of the - // boundary faces can be transformed - // into each other parallel to the - // given direction. - // It is assumed that - // cell1->face(2*direction) and - // cell2->face(2*direction+1) are the - // corresponding boundary faces. - // The optional argument offset will - // be added to each vertex of cell1 - - static const int dim = CellIterator::AccessorType::dimension; - static const int space_dim = CellIterator::AccessorType::space_dimension; - Assert(dim == space_dim, - ExcNotImplemented()); - // ... otherwise we would need two - // directions: One for selecting the - // faces, thus living in dim, - // the other direction for selecting the - // direction for the orthogonal - // projection, thus living in - // space_dim. - - for (int i = 0; i < space_dim; ++i) + Assert (0<=direction && direction::vertices_per_face; ++v) - if (fabs ( cell1->face(2*direction)->vertex(v)(i) - + offset[i] - - cell2->face(2*direction+1)->vertex(v)(i)) - > 1.e-10) - return false; + if (fabs(point1(i) + offset[i] - point2(i)) > 1.e-10) + return false; } return true; } - template - std::map - collect_periodic_cell_pairs (const CellIterator &begin, - const typename identity::type &end, - const types::boundary_id boundary_component, - int direction, - const dealii::Tensor<1,CellIterator::AccessorType::space_dimension> - &offset) + + /* + * Internally used in orthogonal_equality + * + * A lookup table to transform vertex matchings to orientation flags of + * the form (face_orientation, face_flip, face_rotation) + * + * See the comment on the next function as well as the detailed + * documentation of make_periodicity_constraints and + * collect_periodic_face_pairs for details + */ + template struct OrientationLookupTable {}; + + template<> struct OrientationLookupTable<1> { - static const int space_dim = CellIterator::AccessorType::space_dimension; - Assert (0<=direction && direction::vertices_per_face> MATCH_T; + static inline std::bitset<3> lookup (const MATCH_T &) + { + // The 1D case is trivial + return 4; // [true ,false,false] + } + }; + + template<> struct OrientationLookupTable<2> + { + typedef std_cxx1x::array::vertices_per_face> MATCH_T; + static inline std::bitset<3> lookup (const MATCH_T &matching) + { + // In 2D matching faces (=lines) results in two cases: Either + // they are aligned or flipped. We store this "line_flip" + // property somewhat sloppy as "face_flip" + // (always: face_orientation = true, face_rotation = false) + + static const MATCH_T m_tff = {{ 0 , 1 }}; + if (matching == m_tff) return 1; // [true ,false,false] + static const MATCH_T m_ttf = {{ 1 , 0 }}; + if (matching == m_ttf) return 3; // [true ,true ,false] + AssertThrow(false, ExcInternalError()); + } + }; + + template<> struct OrientationLookupTable<3> + { + typedef std_cxx1x::array::vertices_per_face> MATCH_T; + static inline std::bitset<3> lookup (const MATCH_T &matching) + { + // The full fledged 3D case. *Yay* + // See the documentation in include/deal.II/base/geometry_info.h + // as well as the actual implementation in source/grid/tria.cc + // for more details... + + static const MATCH_T m_tff = {{ 0 , 1 , 2 , 3 }}; + if (matching == m_tff) return 1; // [true ,false,false] + static const MATCH_T m_tft = {{ 1 , 3 , 0 , 2 }}; + if (matching == m_tft) return 5; // [true ,false,true ] + static const MATCH_T m_ttf = {{ 3 , 2 , 1 , 0 }}; + if (matching == m_ttf) return 3; // [true ,true ,false] + static const MATCH_T m_ttt = {{ 2 , 0 , 3 , 1 }}; + if (matching == m_ttt) return 7; // [true ,true ,true ] + static const MATCH_T m_fff = {{ 0 , 2 , 1 , 3 }}; + if (matching == m_fff) return 0; // [false,false,false] + static const MATCH_T m_fft = {{ 2 , 3 , 0 , 1 }}; + if (matching == m_fft) return 4; // [false,false,true ] + static const MATCH_T m_ftf = {{ 3 , 1 , 2 , 0 }}; + if (matching == m_ftf) return 2; // [false,true ,false] + static const MATCH_T m_ftt = {{ 1 , 0 , 3 , 2 }}; + if (matching == m_ftt) return 6; // [false,true ,true ] + AssertThrow(false, ExcInternalError()); + } + }; + - std::set cells1; - std::set cells2; - // collect boundary cells, i.e. cells - // with face(2*direction), resp. - // face(2*direction+1) being on the - // boundary and having the correct - // color: - CellIterator cell = begin; - for (; cell!= end; ++cell) + template + inline bool + orthogonal_equality (std::bitset<3> &orientation, + const FaceIterator &face1, + const FaceIterator &face2, + const int direction, + const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset) + { + static const int dim = FaceIterator::AccessorType::dimension; + + // Do a full matching of the face vertices: + + std_cxx1x:: + array::vertices_per_face> matching; + + std::set face2_vertices; + for (unsigned int i = 0; i < GeometryInfo::vertices_per_face; ++i) + face2_vertices.insert(i); + + for (unsigned int i = 0; i < GeometryInfo::vertices_per_face; ++i) { - if (cell->face(2*direction)->at_boundary() && - cell->face(2*direction)->boundary_indicator() == boundary_component) + for (std::set::iterator it = face2_vertices.begin(); + it != face2_vertices.end(); + it++) { - cells1.insert(cell); - } - if (cell->face(2*direction+1)->at_boundary() && - cell->face(2*direction+1)->boundary_indicator() == boundary_component) - { - cells2.insert(cell); + if (orthogonal_equality(face1->vertex(i),face2->vertex(*it), + direction, offset)) + { + matching[i] = *it; + face2_vertices.erase(it); + break; // jump out of the innermost loop + } } } - Assert (cells1.size() == cells2.size(), - ExcMessage ("Unmatched cells on periodic boundaries")); + // And finally, a lookup to determine the ordering bitmask: + if (face2_vertices.empty()) + orientation = OrientationLookupTable::lookup(matching); + + return face2_vertices.empty(); + } + + + + template + inline bool + orthogonal_equality (const FaceIterator &face1, + const FaceIterator &face2, + const int direction, + const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset) + { + // Call the function above with a dummy orientation array + std::bitset<3> dummy; + return orthogonal_equality (dummy, face1, face2, direction, offset); + } + + - std::map matched_cells; + /* + * Internally used in 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) + { + static const int space_dim = FaceIterator::AccessorType::space_dimension; + Assert (0<=direction && direction > ResultPair; + std::map matched_faces; - // Match with a complexity of O(n^2). - // This could be improved... - typedef typename std::set::const_iterator SetIterator; - for (SetIterator it1 = cells1.begin(); it1 != cells1.end(); ++it1) + // Match with a complexity of O(n^2). This could be improved... + std::bitset<3> orientation; + typedef typename std::set::const_iterator SetIterator; + for (SetIterator it1 = faces1.begin(); it1 != faces1.end(); ++it1) { - for (SetIterator it2 = cells2.begin(); it2 != cells2.end(); ++it2) + for (SetIterator it2 = faces2.begin(); it2 != faces2.end(); ++it2) { - if (orthogonal_equality(*it1, *it2, direction, offset)) + if (GridTools::orthogonal_equality(orientation, *it1, *it2, + direction, offset)) { - // We have a match, so insert the - // matching pairs and remove the matched - // cell in cells2 to speed up the + // We have a match, so insert the matching pairs and + // remove the matched cell in faces2 to speed up the // matching: - matched_cells[*it1] = *it2; - cells2.erase(it2); + matched_faces[*it1] = std::make_pair(*it2, orientation); + faces2.erase(it2); break; } } } - Assert (matched_cells.size() == cells1.size() && - cells2.size() == 0, - ExcMessage ("Unmatched cells on periodic boundaries")); + AssertThrow (matched_faces.size() == faces1.size() && faces2.size() == 0, + ExcMessage ("Unmatched faces on periodic boundaries")); - return matched_cells; + return matched_faces; } + + 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, + int direction, + const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset) + { + static const int space_dim = FaceIterator::AccessorType::space_dimension; + Assert (0<=direction && direction faces1; + std::set faces2; + + for (FaceIterator face = begin; face!= end; ++face) + { + 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")); + + // and call match_periodic_face_pairs that does the actual matching: + + return match_periodic_face_pairs(faces1, faces2, direction, offset); + } + + + template - std::map - collect_periodic_cell_pairs (const DH &dof_handler, - const types::boundary_id boundary_component, - int direction, - const dealii::Tensor<1,DH::space_dimension> - &offset) + std::map > > + collect_periodic_face_pairs (const DH &dof_handler, + const types::boundary_id b_id1, + const types::boundary_id b_id2, + int direction, + const dealii::Tensor<1,DH::space_dimension> &offset) { - return collect_periodic_cell_pairs - (dof_handler.begin_active(), - dof_handler.end(), - boundary_component, - direction, - 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) + { + 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")); + + // and call match_periodic_face_pairs that does the actual matching: + + return match_periodic_face_pairs(faces1, faces2, direction, offset); } diff --git a/deal.II/source/grid/grid_tools.inst.in b/deal.II/source/grid/grid_tools.inst.in index d8c6f0e51a..fe5ce669bb 100644 --- a/deal.II/source/grid/grid_tools.inst.in +++ b/deal.II/source/grid/grid_tools.inst.in @@ -184,22 +184,65 @@ for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II namespace GridTools \{ template - std::map - collect_periodic_cell_pairs(const X::cell_iterator &, - const X::cell_iterator &, - const types::boundary_id, - int, - const Tensor<1,deal_II_space_dimension> &); + bool orthogonal_equality (std::bitset<3> &, + const X::active_face_iterator&, + const X::active_face_iterator&, + const int, + const Tensor<1,deal_II_space_dimension> &); template - std::map - collect_periodic_cell_pairs (const X&, - const types::boundary_id, - int, - const Tensor<1,deal_II_space_dimension> &); + bool orthogonal_equality (std::bitset<3> &, + const X::face_iterator&, + const X::face_iterator&, + const int, + const Tensor<1,deal_II_space_dimension> &); + template + bool orthogonal_equality (const X::active_face_iterator&, + const X::active_face_iterator&, + const int, + const Tensor<1,deal_II_space_dimension> &); + + template + bool orthogonal_equality (const X::face_iterator&, + const X::face_iterator&, + const int, + const Tensor<1,deal_II_space_dimension> &); + + #if deal_II_dimension >= 2 + + template + std::map > > + collect_periodic_face_pairs + (const X::face_iterator &, + const X::face_iterator &, + const types::boundary_id, + const types::boundary_id, + int, + const Tensor<1,deal_II_space_dimension> &); + + template + std::map > > + collect_periodic_face_pairs + (const X::active_face_iterator &, + const X::active_face_iterator &, + const types::boundary_id, + const types::boundary_id, + int, + const Tensor<1,deal_II_space_dimension> &); + + template + std::map > > + collect_periodic_face_pairs + (const X &, + const types::boundary_id, + const types::boundary_id, + int, + const Tensor<1,deal_II_space_dimension> &); + + #endif \} - #endif +#endif } -- 2.39.5