From 7878df99e97c19cd1e6a0e9c785b489e70981cb3 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Tue, 22 Mar 2016 17:44:45 +0100 Subject: [PATCH] Use periodic_face_map also for the multigrid structures in distributed/tria.h, compare --- include/deal.II/distributed/tria.h | 4 +- include/deal.II/grid/tria.h | 33 +++-- source/distributed/tria.cc | 226 ++++++++++++++++++++++++++--- source/grid/tria.cc | 17 ++- tests/mpi/periodicity_04.cc | 7 +- 5 files changed, 244 insertions(+), 43 deletions(-) diff --git a/include/deal.II/distributed/tria.h b/include/deal.II/distributed/tria.h index 4ed9440966..3a12ebf34f 100644 --- a/include/deal.II/distributed/tria.h +++ b/include/deal.II/distributed/tria.h @@ -946,7 +946,7 @@ namespace parallel */ void fill_level_vertices_with_ghost_neighbors - (const unsigned int level, + (const int level, std::map > &vertices_with_ghost_neighbors); @@ -960,7 +960,7 @@ namespace parallel * Used by DoFHandler::Policy::ParallelDistributed. */ std::vector - mark_locally_active_vertices_on_level(const unsigned int level) const; + mark_locally_active_vertices_on_level(const int level) const; template friend class dealii::internal::DoFHandler::Policy::ParallelDistributed; }; diff --git a/include/deal.II/grid/tria.h b/include/deal.II/grid/tria.h index a20531642d..7ecbe2acf9 100644 --- a/include/deal.II/grid/tria.h +++ b/include/deal.II/grid/tria.h @@ -1471,19 +1471,6 @@ public: distorted_cells; }; - /** - * 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 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. - */ - std::map, std::pair, std::bitset<3> > > periodic_face_map; - /** * Make the dimension available in function templates. */ @@ -2984,7 +2971,6 @@ 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 periodic_face_map that stores the active periodic faces. * * The vector can be filled by the function * GridTools::collect_periodic_faces. @@ -3000,6 +2986,12 @@ public: add_periodicity (const std::vector > &); + /** + * Return the periodic_face_map. + */ + const std::map,std::pair, std::bitset<3> > > & + get_periodic_face_map() const; + BOOST_SERIALIZATION_SPLIT_MEMBER() @@ -3113,6 +3105,19 @@ protected: private: + /** + * 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 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. + */ + std::map, std::pair, std::bitset<3> > > periodic_face_map; + /** * @name Cell iterator functions for internal use * @{ diff --git a/source/distributed/tria.cc b/source/distributed/tria.cc index 0bdab3a7b8..2e052d49e2 100644 --- a/source/distributed/tria.cc +++ b/source/distributed/tria.cc @@ -2244,7 +2244,7 @@ namespace parallel } this->update_number_cache (); - Triangulation::update_periodic_face_map(); + this->update_periodic_face_map(); } @@ -2773,7 +2773,7 @@ namespace parallel } this->update_number_cache (); - Triangulation::update_periodic_face_map(); + this->update_periodic_face_map(); } @@ -3197,6 +3197,13 @@ namespace parallel topological_vertex_numbering[periodic.cell[0]->face(periodic.face_idx[0])->vertex_index(vface0)] = topological_vertex_numbering[periodic.cell[1]->face(periodic.face_idx[1])->vertex_index(v)] = min_index; + + /* std::cout << "Old identified vertices " << periodic.cell[0]->face(periodic.face_idx[0])->vertex_index(vface0) + << " at " << periodic.cell[0]->face(periodic.face_idx[0])->vertex(vface0) + << " and " << periodic.cell[1]->face(periodic.face_idx[1])->vertex_index(v) + << " at " << periodic.cell[1]->face(periodic.face_idx[1])->vertex(v) + << std::endl;*/ + } } @@ -3206,10 +3213,9 @@ namespace parallel // artificial cell layer (the active cells are taken care of by p4est) template bool enforce_mesh_balance_over_periodic_boundaries - (Triangulation &tria, - const std::vector::cell_iterator> > &periodic_face_pairs_level_0) + (Triangulation &tria) { - if (periodic_face_pairs_level_0.empty()) + if (tria.get_periodic_face_map().size()==0) return false; std::vector flags_before[2]; @@ -3217,14 +3223,100 @@ namespace parallel tria.save_refine_flags (flags_before[1]); std::vector topological_vertex_numbering(tria.n_vertices()); + std::vector topological_vertex_numbering_cmp(tria.n_vertices()); for (unsigned int i=0; i::cell_iterator cell_iterator; + typename std::map, + std::pair, std::bitset<3> > >::const_iterator it; + for (it = tria.get_periodic_face_map().begin(); it!= tria.get_periodic_face_map().end(); ++it) + { + const cell_iterator &cell_1 = it->first.first; + const unsigned int face_no_1 = it->first.second; + const cell_iterator &cell_2 = it->second.first.first; + const unsigned int face_no_2 = it->second.first.second; + const std::bitset<3> face_orientation = it->second.second; + + if (cell_1->level() == cell_2->level()) + { + for (unsigned int v=0; v::vertices_per_cell; ++v) + { + // take possible non-standard orientation of face on cell[0] into + // account + const unsigned int vface0 = + GeometryInfo::standard_to_real_face_vertex(v,face_orientation[0], + face_orientation[1], + face_orientation[2]); + const unsigned int vi0 = topological_vertex_numbering[cell_1->face(face_no_1)->vertex_index(vface0)]; + const unsigned int vi1 = topological_vertex_numbering[cell_2->face(face_no_2)->vertex_index(v)]; + const unsigned int min_index = std::min(vi0, vi1); + topological_vertex_numbering[cell_1->face(face_no_1)->vertex_index(vface0)] + = topological_vertex_numbering[cell_2->face(face_no_2)->vertex_index(v)] + = min_index; + /* std::cout << "Newly identified vertices " << cell_1->face(face_no_1)->vertex_index(vface0) + << " at " << cell_1->face(face_no_1)->vertex(vface0) + << " and " << cell_2->face(face_no_2)->vertex_index(v) + << " at " << cell_2->face(face_no_2)->vertex(v) + << std::endl;*/ + } + } + } + + // There must not be any chains! + for (unsigned int i=0; idealii::Triangulation::prepare_coarsening_and_refinement(); - + this->update_periodic_face_map(); // enforce 2:1 mesh balance over periodic boundaries if (this->smooth_grid & dealii::Triangulation::limit_level_difference_at_vertices) - mesh_changed = enforce_mesh_balance_over_periodic_boundaries(*this, - this->periodic_face_pairs_level_0); + mesh_changed = enforce_mesh_balance_over_periodic_boundaries(*this); } while (mesh_changed); @@ -3895,7 +3986,7 @@ namespace parallel refinement_in_progress = false; this->update_number_cache (); - Triangulation::update_periodic_face_map(); + this->update_periodic_face_map(); } template @@ -3965,7 +4056,7 @@ namespace parallel // update how many cells, edges, etc, we store locally this->update_number_cache (); - Triangulation::update_periodic_face_map(); + this->update_periodic_face_map(); } @@ -4678,7 +4769,7 @@ namespace parallel void Triangulation:: fill_level_vertices_with_ghost_neighbors - (const unsigned int level, + (const int level, std::map > &vertices_with_ghost_neighbors) { @@ -4694,9 +4785,72 @@ namespace parallel vertices_with_ghost_neighbors[cell->vertex_index(v)] .insert (cell->level_subdomain_id()); + std::map > + vertices_with_ghost_neighbors_cmp = vertices_with_ghost_neighbors; + + 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); + level, vertices_with_ghost_neighbors_cmp); + + //now for the vertices on periodic faces + typename std::map, + std::pair, std::bitset<3> > >::const_iterator it; + + for (it = this->get_periodic_face_map().begin(); it!= this->get_periodic_face_map().end(); ++it) + { + const cell_iterator &cell_1 = it->first.first; + const unsigned int face_no_1 = it->first.second; + const cell_iterator &cell_2 = it->second.first.first; + const unsigned int face_no_2 = it->second.first.second; + const std::bitset<3> face_orientation = it->second.second; + + if (cell_1->level() == level && + cell_2->level() == level) + { + for (unsigned int v=0; v::vertices_per_cell; ++v) + { + // take possible non-standard orientation of faces into account + const unsigned int vface0 = + GeometryInfo::standard_to_real_face_vertex(v,face_orientation[0], + face_orientation[1], + face_orientation[2]); + const unsigned int idx0 = cell_1->face(face_no_1)->vertex_index(vface0); + const unsigned int idx1 = cell_2->face(face_no_2)->vertex_index(v); + if (vertices_with_ghost_neighbors.find(idx0) != vertices_with_ghost_neighbors.end()) + vertices_with_ghost_neighbors[idx1].insert(vertices_with_ghost_neighbors[idx0].begin(), + vertices_with_ghost_neighbors[idx0].end()); + if (vertices_with_ghost_neighbors.find(idx1) != vertices_with_ghost_neighbors.end()) + vertices_with_ghost_neighbors[idx0].insert(vertices_with_ghost_neighbors[idx1].begin(), + vertices_with_ghost_neighbors[idx1].end()); + } + } + } + + if (vertices_with_ghost_neighbors != vertices_with_ghost_neighbors_cmp) + { + std::cout << "Print vertices_with_ghost_neighbors: " << std::endl; + std::map >::const_iterator it; + for (it = vertices_with_ghost_neighbors.begin(); it!=vertices_with_ghost_neighbors.end(); ++it) + { + std::cout << "Vertex " << it->first << " is on subdomains "; + std::set::const_iterator it2; + for (it2 = it->second.begin(); it2!=it->second.end(); ++it2) + std::cout << *it2 << " "; + std::cout << std::endl; + } + + std::cout << "Print vertices_with_ghost_neighbors_cmp: " << std::endl; + for (it = vertices_with_ghost_neighbors_cmp.begin(); it!=vertices_with_ghost_neighbors_cmp.end(); ++it) + { + std::cout << "Vertex " << it->first << " is on subdomains "; + std::set::const_iterator it2; + for (it2 = it->second.begin(); it2!=it->second.end(); ++it2) + std::cout << *it2 << " "; + std::cout << std::endl; + } + Assert(false, ExcInternalError()); + } } @@ -4704,7 +4858,7 @@ namespace parallel template std::vector Triangulation - ::mark_locally_active_vertices_on_level (const unsigned int level) const + ::mark_locally_active_vertices_on_level (const int level) const { Assert (dim>1, ExcNotImplemented()); @@ -4716,10 +4870,48 @@ namespace parallel for (unsigned int v=0; v::vertices_per_cell; ++v) marked_vertices[cell->vertex_index(v)] = true; + std::vector marked_vertices_cmp = marked_vertices; + 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); + level, marked_vertices_cmp); + + /** + * ensure that if one of the two vertices on a periodic face is marked + * as active (i.e., belonging to an owned level cell), also the other + * one is active + */ + typename std::map, + std::pair, std::bitset<3> > >::const_iterator it; + + for (it = this->get_periodic_face_map().begin(); it!= this->get_periodic_face_map().end(); ++it) + { + const cell_iterator &cell_1 = it->first.first; + const unsigned int face_no_1 = it->first.second; + const cell_iterator &cell_2 = it->second.first.first; + const unsigned int face_no_2 = it->second.first.second; + const std::bitset<3> &face_orientation = it->second.second; + + if (cell_1->level() == level && + cell_2->level() == level) + { + for (unsigned int v=0; v::vertices_per_cell; ++v) + { + // take possible non-standard orientation of faces into account + const unsigned int vface0 = + GeometryInfo::standard_to_real_face_vertex(v,face_orientation[0], + face_orientation[1], + face_orientation[2]); + if (marked_vertices[cell_1->face(face_no_1)->vertex_index(vface0)] || + marked_vertices[cell_2->face(face_no_2)->vertex_index(v)]) + marked_vertices[cell_1->face(face_no_1)->vertex_index(vface0)] + = marked_vertices[cell_2->face(face_no_2)->vertex_index(v)] + = true; + } + } + } + Assert(marked_vertices == marked_vertices_cmp, ExcInternalError()); return marked_vertices; } @@ -4965,7 +5157,7 @@ namespace parallel } this->update_number_cache (); - Triangulation::update_periodic_face_map(); + this->update_periodic_face_map(); } diff --git a/source/grid/tria.cc b/source/grid/tria.cc index 4e06ab94fe..99a3e47e25 100644 --- a/source/grid/tria.cc +++ b/source/grid/tria.cc @@ -11744,6 +11744,14 @@ Triangulation::add_periodicity update_periodic_face_map(); } +template +const typename std::map::cell_iterator, unsigned int>, + std::pair::cell_iterator, unsigned int>, std::bitset<3> > > & + Triangulation::get_periodic_face_map() const +{ + return periodic_face_map; +} + template void @@ -11784,12 +11792,7 @@ 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. - const parallel::Triangulation< dim, spacedim > *distributed_triangulation - = dynamic_cast *> (this); - if (!distributed_triangulation) - update_periodic_face_map(); + update_periodic_face_map(); } @@ -11819,7 +11822,7 @@ Triangulation::update_periodic_face_map () //first empty the currently stored objects periodic_face_map.clear(); - typename std::vector::cell_iterator> >::const_iterator it; + typename std::vector >::const_iterator it; for (it=periodic_face_pairs_level_0.begin(); it!=periodic_face_pairs_level_0.end(); ++it) { update_periodic_face_map_recursively diff --git a/tests/mpi/periodicity_04.cc b/tests/mpi/periodicity_04.cc index 7211af7135..844355036c 100644 --- a/tests/mpi/periodicity_04.cc +++ b/tests/mpi/periodicity_04.cc @@ -241,7 +241,7 @@ void check triangulation.execute_coarsening_and_refinement(); typedef std::pair::cell_iterator, unsigned int> CellFace; - const typename std::map > > &face_map = triangulation.periodic_face_map; + const typename std::map > > &face_map = triangulation.get_periodic_face_map(); typename std::map > >::const_iterator it; int sum_of_pairs_local = face_map.size(); int sum_of_pairs_global; @@ -263,8 +263,9 @@ void check { std::cout << "face_center_1: " << face_center_1 << std::endl; std::cout << "face_center_2: " << face_center_2 << std::endl; - for (auto it = triangulation.periodic_face_map.begin(); - it != triangulation.periodic_face_map.end(); ++it) + typename std::map > >::const_iterator it; + for (it = triangulation.get_periodic_face_map().begin(); + it != triangulation.get_periodic_face_map().end(); ++it) { std::cout << "The cell with center " << it->first.first->center() << " has on face " << it->first.second -- 2.39.5