From 4987b9227b23b7abbf89bb7c8b8021e81e546fc8 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Tue, 22 Mar 2016 15:46:11 +0100 Subject: [PATCH] Store orientation as well --- include/deal.II/grid/tria.h | 3 +- source/grid/tria.cc | 95 +++++++++++++++++-------------------- 2 files changed, 46 insertions(+), 52 deletions(-) diff --git a/include/deal.II/grid/tria.h b/include/deal.II/grid/tria.h index 83665a8c2c..a20531642d 100644 --- a/include/deal.II/grid/tria.h +++ b/include/deal.II/grid/tria.h @@ -38,6 +38,7 @@ DEAL_II_ENABLE_EXTRA_DIAGNOSTICS #include #include #include +#include DEAL_II_NAMESPACE_OPEN @@ -1481,7 +1482,7 @@ public: /** * If add_periodicity() is called, this variable stores the active periodic face pairs. */ - std::map, std::pair > periodic_face_map; + std::map, std::pair, std::bitset<3> > > periodic_face_map; /** * Make the dimension available in function templates. diff --git a/source/grid/tria.cc b/source/grid/tria.cc index c7760c1c8f..4e06ab94fe 100644 --- a/source/grid/tria.cc +++ b/source/grid/tria.cc @@ -914,14 +914,18 @@ namespace (const typename Triangulation::cell_iterator &cell_1, const typename Triangulation::cell_iterator &cell_2, unsigned int n_face_1, unsigned int n_face_2, - const bool face_orientation, const bool face_flip, const bool face_rotation, + const std::bitset<3> &orientation, typename std::map::cell_iterator, unsigned int>, - std::pair::cell_iterator, unsigned int> > &periodic_face_map) + std::pair::cell_iterator, unsigned int>, std::bitset<3> > > &periodic_face_map) { typedef typename Triangulation::face_iterator FaceIterator; const FaceIterator face_1 = cell_1->face(n_face_1); const FaceIterator face_2 = cell_2->face(n_face_2); + const bool face_orientation = orientation[0]; + const bool face_flip = orientation[1]; + const bool face_rotation = orientation[2]; + 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")); @@ -938,7 +942,9 @@ namespace typedef std::pair::cell_iterator, unsigned int> CellFace; const CellFace cell_face_1 (cell_1, n_face_1); const CellFace cell_face_2 (cell_2, n_face_2); - const std::pair periodic_faces (cell_face_1, cell_face_2); + const std::pair > cell_face_orientation_2 (cell_face_2, orientation); + + const std::pair > > periodic_faces (cell_face_1, cell_face_orientation_2); // Only one periodic neighbor is allowed Assert(periodic_face_map.count(cell_face_1) == 0, ExcInternalError()); @@ -1003,18 +1009,22 @@ namespace // find subcell ids that belong to the subface indices unsigned int child_cell_1 - = GeometryInfo::child_cell_on_face(cell_1->refinement_case(), n_face_1, i, face_orientation, - face_flip, face_rotation, face_1->refinement_case()); + = GeometryInfo::child_cell_on_face + (cell_1->refinement_case(), n_face_1, i, cell_1->face_orientation(n_face_1), + cell_1->face_flip(n_face_1), cell_1->face_rotation(n_face_1), face_1->refinement_case()); unsigned int child_cell_2 - = GeometryInfo::child_cell_on_face(cell_2->refinement_case(), n_face_2, j, face_orientation, - face_flip, face_rotation, face_2->refinement_case()); + = GeometryInfo::child_cell_on_face + (cell_2->refinement_case(), n_face_2, j, cell_2->face_orientation(n_face_2), + cell_2->face_flip(n_face_2), cell_2->face_rotation(n_face_2), face_2->refinement_case()); + + Assert(cell_1->child(child_cell_1)->face(n_face_1) == face_1->child(i), ExcInternalError()); + Assert(cell_2->child(child_cell_2)->face(n_face_2) == face_2->child(j), ExcInternalError()); // precondition: subcell has the same orientation as cell (so that the face numbers coincide) // recursive call update_periodic_face_map_recursively (cell_1->child(child_cell_1), cell_2->child(child_cell_2), - n_face_1, n_face_2, face_orientation, face_flip, face_rotation, - periodic_face_map); + n_face_1, n_face_2, orientation, periodic_face_map); } } else //only face_1 has children @@ -1023,14 +1033,13 @@ namespace { // find subcell ids that belong to the subface indices unsigned int child_cell_1 - = GeometryInfo::child_cell_on_face(cell_1->refinement_case(), n_face_1, i, face_orientation, - face_flip, face_rotation, face_1->refinement_case()); + = GeometryInfo::child_cell_on_face(cell_1->refinement_case(), n_face_1, i, cell_1->face_orientation(n_face_1), + cell_1->face_flip(n_face_1), cell_1->face_rotation(n_face_1), face_1->refinement_case()); // recursive call update_periodic_face_map_recursively (cell_1->child(child_cell_1), cell_2, - n_face_1, n_face_2, face_orientation, face_flip, face_rotation, - periodic_face_map); + n_face_1, n_face_2, orientation, periodic_face_map); } } @@ -11732,22 +11741,7 @@ Triangulation::add_periodicity periodicity_vector.end()); //Now initialize periodic_face_map - typename std::vector >::const_iterator it; - for (it=periodicity_vector.begin(); it!=periodicity_vector.end(); ++it) - { - const CellFace cell_face_1 (it->cell[0], it->face_idx[0]); - const CellFace cell_face_2 (it->cell[1], it->face_idx[1]); - - //Only one periodic neighbor is allowed - Assert(periodic_face_map.count(cell_face_1)==0, ExcInternalError()); - Assert(periodic_face_map.count(cell_face_2)==0, ExcInternalError()); - - const std::pair periodic_faces_1 (cell_face_1, cell_face_2); - const std::pair periodic_faces_2 (cell_face_2, cell_face_1); - - periodic_face_map.insert(periodic_faces_1); - periodic_face_map.insert(periodic_faces_2); - } + update_periodic_face_map(); } @@ -11830,38 +11824,37 @@ Triangulation::update_periodic_face_map () { update_periodic_face_map_recursively (it->cell[0], it->cell[1], it->face_idx[0], it->face_idx[1], - it->orientation[0], it->orientation[1], it->orientation[2], - periodic_face_map); + it->orientation, periodic_face_map); + //for the other way, we need to invert the orientation - bool orientation, flip, rotation; - orientation = it->orientation[0]; - rotation = it->orientation[2]; - if (it->orientation[0] && it->orientation [2]) - flip = !it->orientation[1]; - else - flip = it->orientation[1]; + std::bitset<3> inverted_orientation; + { + bool orientation, flip, rotation; + orientation = it->orientation[0]; + rotation = it->orientation[2]; + flip = orientation ? rotation ^ it->orientation[1] : it->orientation[1]; + inverted_orientation[0] = orientation; + inverted_orientation[1] = flip; + inverted_orientation[2] = rotation; + } update_periodic_face_map_recursively (it->cell[1], it->cell[0], it->face_idx[1], it->face_idx[0], - orientation, flip, rotation, - periodic_face_map); + inverted_orientation, periodic_face_map); } //check consistency - std::cout << "check consistency" << std::endl; - //cell->neighbor(neighbor)->neighbor_child_on_subface(face_no, subface_no)==cell. typename std::map, - std::pair >::const_iterator it_test; + std::pair, std::bitset<3> > >::const_iterator it_test; for (it_test=periodic_face_map.begin(); it_test!=periodic_face_map.end(); ++it_test) { - if (!it_test->second.first->has_children()) + const Triangulation::cell_iterator cell_1 = it_test->first.first; + const Triangulation::cell_iterator cell_2 = it_test->second.first.first; + if (cell_1->level() == cell_2->level()) { - if (!periodic_face_map[it_test->second].first->has_children()) - { - // if the neighbor is active as well, the same pair - // order swapped has to be in the map - Assert(periodic_face_map[it_test->second] == it_test->first, - ExcInternalError()); - } + // if both cells have the same neighbor, then the same pair + // order swapped has to be in the map + Assert (periodic_face_map[it_test->second.first].first == it_test->first, + ExcInternalError()); } } } -- 2.39.5