From eaaa319c58057a3b49fb3fe994ed65aa0d176cf1 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Thu, 25 Feb 2016 19:51:49 +0100 Subject: [PATCH] Implement updating periodic_face_map --- include/deal.II/distributed/tria.h | 7 - include/deal.II/grid/tria.h | 20 ++- source/distributed/tria.cc | 22 +-- source/grid/tria.cc | 238 ++++++++++++++++++++++++++++- 4 files changed, 264 insertions(+), 23 deletions(-) diff --git a/include/deal.II/distributed/tria.h b/include/deal.II/distributed/tria.h index 95e5859e65..4ed9440966 100644 --- a/include/deal.II/distributed/tria.h +++ b/include/deal.II/distributed/tria.h @@ -878,13 +878,6 @@ namespace parallel std::vector coarse_cell_to_p4est_tree_permutation; std::vector p4est_tree_to_coarse_cell_permutation; - /** - * If add_periodicity() is called, this variable stores the given - * periodic face pairs on level 0 for later access during the - * identification of ghost cells for the multigrid hierarchy. - */ - std::vector > periodic_face_pairs_level_0; - /** * Return a pointer to the p4est tree that belongs to the given * dealii_coarse_cell_index() diff --git a/include/deal.II/grid/tria.h b/include/deal.II/grid/tria.h index 0678a6abfb..bd6ece4cf1 100644 --- a/include/deal.II/grid/tria.h +++ b/include/deal.II/grid/tria.h @@ -46,6 +46,11 @@ template class Boundary; template class StraightBoundary; template class Manifold; +namespace GridTools +{ + template struct PeriodicFacePair; +} + template class TriaAccessor; template class TriaAccessor<0,1,spacedim>; template class TriaAccessorBase; @@ -1468,15 +1473,15 @@ public: /** * If add_periodicity() is called, this variable stores the given * periodic face pairs on level 0 for later access during the - * identification of ghost cells for the multigrid hierarchy. + * identification of ghost cells for the multigrid hierarchy and for + * setting up the periodic_face_map. */ std::vector > periodic_face_pairs_level_0; /** * If add_periodicity() is called, this variable stores the active periodic face pairs. - * The variable is updated on refinement. */ - std::map, GridTools::PeriodicFacePair > periodic_face_map; + std::map, std::pair > periodic_face_map; /** * Make the dimension available in function templates. @@ -2978,7 +2983,7 @@ public: * Declare the (coarse) face pairs given in the argument of this function * as periodic. This way it it possible to obtain neighbors across periodic * boundaries. - * This function initializes ActivePeriodicFacePairs and on refinement this + * This function initializes periodic_face_map that stores the active periodic faces. * * The vector can be filled by the function * GridTools::collect_periodic_faces. @@ -3100,6 +3105,13 @@ protected: const unsigned int magic_number2, std::istream &in); + /** + * Recreate information about periodic neighbors from periodic_face_pairs_level_0. + * Only faces that belong to locally owned or ghosted cells are stored. + */ + void update_periodic_face_map (); + + private: /** * @name Cell iterator functions for internal use diff --git a/source/distributed/tria.cc b/source/distributed/tria.cc index 3423854a29..0bdab3a7b8 100644 --- a/source/distributed/tria.cc +++ b/source/distributed/tria.cc @@ -2244,6 +2244,7 @@ namespace parallel } this->update_number_cache (); + Triangulation::update_periodic_face_map(); } @@ -2582,8 +2583,6 @@ namespace parallel coarse_cell_to_p4est_tree_permutation.resize (0); p4est_tree_to_coarse_cell_permutation.resize (0); - periodic_face_pairs_level_0.clear(); - dealii::Triangulation::clear (); this->update_number_cache (); @@ -2774,6 +2773,7 @@ namespace parallel } this->update_number_cache (); + Triangulation::update_periodic_face_map(); } @@ -3207,7 +3207,7 @@ namespace parallel template bool enforce_mesh_balance_over_periodic_boundaries (Triangulation &tria, - const std::vector::cell_iterator> > periodic_face_pairs_level_0) + const std::vector::cell_iterator> > &periodic_face_pairs_level_0) { if (periodic_face_pairs_level_0.empty()) return false; @@ -3361,7 +3361,7 @@ namespace parallel if (this->smooth_grid & dealii::Triangulation::limit_level_difference_at_vertices) mesh_changed = enforce_mesh_balance_over_periodic_boundaries(*this, - periodic_face_pairs_level_0); + this->periodic_face_pairs_level_0); } while (mesh_changed); @@ -3894,8 +3894,8 @@ namespace parallel refinement_in_progress = false; - this->update_number_cache (); + Triangulation::update_periodic_face_map(); } template @@ -3965,6 +3965,7 @@ namespace parallel // update how many cells, edges, etc, we store locally this->update_number_cache (); + Triangulation::update_periodic_face_map(); } @@ -4693,8 +4694,8 @@ namespace parallel vertices_with_ghost_neighbors[cell->vertex_index(v)] .insert (cell->level_subdomain_id()); - for (unsigned int i=0; iperiodic_face_pairs_level_0.size(); ++i) + set_periodic_ghost_neighbors_recursively(this->periodic_face_pairs_level_0[i], level, vertices_with_ghost_neighbors); } @@ -4715,8 +4716,8 @@ namespace parallel for (unsigned int v=0; v::vertices_per_cell; ++v) marked_vertices[cell->vertex_index(v)] = true; - for (unsigned int i=0; iperiodic_face_pairs_level_0.size(); ++i) + mark_periodic_vertices_recursively(this->periodic_face_pairs_level_0[i], level, marked_vertices); return marked_vertices; @@ -4862,7 +4863,7 @@ namespace parallel } //finally call the base class for storing the periodicity information - Triangulation::add_periodicity(periodicity_vector); + dealii::Triangulation::add_periodicity(periodicity_vector); #else Assert(false, ExcMessage ("Need p4est version >= 0.3.4.1!")); #endif @@ -4964,6 +4965,7 @@ namespace parallel } this->update_number_cache (); + Triangulation::update_periodic_face_map(); } diff --git a/source/grid/tria.cc b/source/grid/tria.cc index 78c065079c..69852420dc 100644 --- a/source/grid/tria.cc +++ b/source/grid/tria.cc @@ -1,4 +1,5 @@ // --------------------------------------------------------------------- + // // Copyright (C) 1998 - 2016 by the deal.II authors // @@ -907,6 +908,158 @@ namespace } } + + template + void update_periodic_face_map_recursively + (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, + typename std::map::cell_iterator, unsigned int>, + std::pair::cell_iterator, unsigned int> > &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); + + 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 != face_2, + ExcMessage ("face_1 and face_2 are equal!")); + + Assert(face_1->at_boundary() &&face_2->at_boundary(), + ExcMessage ("Periodic faces must be on the boundary")); + + // A lookup table on how to go through the child cells depending on the + // orientation: + // see Documentation of GeometryInfo for details + + 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 + } + } + }; + + if (cell_1->has_children()) + { + if (cell_2->has_children()) + { + // In the case that both faces have children, we loop over all + // children and apply update_periodic_face_map_recursively recursively: + + Assert(face_1->n_children() == GeometryInfo::max_children_per_face + && face_2->n_children() == GeometryInfo::max_children_per_face, + ExcNotImplemented()); + + for (unsigned int i = 0; i < GeometryInfo::max_children_per_face; ++i) + { + // Lookup the index for the second face + unsigned int j = 0; + 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()); + } + + // 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()); + 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()); + + // 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); + } + } + else //only face_1 has children + { + for (unsigned int i = 0; i < GeometryInfo::max_children_per_face; ++i) + { + // 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()); + + // 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); + } + } + + } + else + { + // Otherwise at least one of the two faces is active and + // we need to do some work and enter the periodic face pairs! + + //periodic boundary can only be on the boundary of the domain + Assert(cell_1->at_boundary(n_face_1), ExcInternalError()); + Assert(cell_2->at_boundary(n_face_2), ExcInternalError()); + + //Ignore artificial cells, cell_2 may not be active + if (!cell_1->is_artificial()) + { + // if cell_2 is active it must not be artificial + if (!cell_2->has_children() && cell_2->is_artificial()) + { + std::cout<<"The artificial neighbor of the cell with center " <center() + <<" is the cell with center "<center()<is_artificial(), ExcInternalError()); + } + if (cell_2->has_children() || !cell_2->is_artificial()) + { + // insert periodic face pair for both cells + 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); + + // Only one periodic neighbor is allowed + Assert(periodic_face_map.count(cell_face_1) == 0, ExcInternalError()); + periodic_face_map.insert(periodic_faces); + } + } + } + } + }// end of anonymous namespace @@ -8823,6 +8976,8 @@ void Triangulation::clear () { clear_despite_subscriptions(); signals.clear(); + periodic_face_pairs_level_0.clear(); + periodic_face_map.clear(); } @@ -11588,14 +11743,36 @@ Triangulation::get_triangulation () const template -virtual void -Triangulation::add_periodicity() +void +Triangulation::add_periodicity +(const std::vector > & + periodicity_vector) { + typedef std::pair CellFace; periodic_face_pairs_level_0.insert(periodic_face_pairs_level_0.end(), periodicity_vector.begin(), 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); + } } + template void Triangulation::execute_coarsening_and_refinement () @@ -11634,6 +11811,14 @@ Triangulation::execute_coarsening_and_refinement () AssertThrow (cells_with_distorted_children.distorted_cells.size() == 0, cells_with_distorted_children); + + // For parallel::distributed::Triangulations we update + // periodic_face_map later. At this point cells might be + // artificial which prevents us from storing the correct faces + const parallel::Triangulation< dim, spacedim > *distributed_triangulation + = dynamic_cast *> (this); + if (!distributed_triangulation) + update_periodic_face_map(); } @@ -11656,6 +11841,55 @@ Triangulation::reset_active_cell_indices () } +template +void +Triangulation::update_periodic_face_map () +{ + //first empty the currently stored objects + periodic_face_map.clear(); + + typename std::vector::cell_iterator> >::const_iterator it; + for (it=periodic_face_pairs_level_0.begin(); it!=periodic_face_pairs_level_0.end(); ++it) + { + 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); + //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]; + 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); + } + + //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; + for (it_test=periodic_face_map.begin(); it_test!=periodic_face_map.end(); ++it_test) + { + if (!it_test->second.first->has_children()) + { + 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()); + } + } + } +} + + template void -- 2.39.5