From 51876e698544d8cabc69a214a6a9207e316eff82 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Sun, 18 Oct 2015 09:43:47 +0200 Subject: [PATCH] Fix ghost layer on MG levels for periodic boundary conditions --- include/deal.II/distributed/tria.h | 57 ++++++- source/distributed/tria.cc | 261 +++++++++++++++++++++++++++-- source/dofs/dof_handler_policy.cc | 84 ++-------- 3 files changed, 311 insertions(+), 91 deletions(-) diff --git a/include/deal.II/distributed/tria.h b/include/deal.II/distributed/tria.h index 291101d333..229c5e864d 100644 --- a/include/deal.II/distributed/tria.h +++ b/include/deal.II/distributed/tria.h @@ -915,12 +915,21 @@ namespace parallel * grid cells are badly ordered this may mean that individual parts of * the forest stored on a local machine may be split across coarse grid * cells that are not geometrically close. Consequently, we apply a - * Cuthill-McKee preordering to ensure that the part of the forest - * stored by p4est is located on geometrically close coarse grid cells. + * hierarchical preordering according to + * SparsityTools::reorder_hierarchical() to ensure that the part of the + * forest stored by p4est is located on geometrically close coarse grid + * cells. */ 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() @@ -957,8 +966,8 @@ namespace parallel void attach_mesh_data(); /** - * fills a map that, for each vertex, lists all the processors whose - * subdomains are adjacent to that vertex. Used by + * Fills a map that, for each vertex, lists all the processors whose + * subdomains are adjacent to that vertex. Used by * DoFHandler::Policy::ParallelDistributed. */ void @@ -966,6 +975,29 @@ namespace parallel (std::map > &vertices_with_ghost_neighbors); + /** + * Fills a map that, for each vertex, lists all the processors whose + * subdomains are adjacent to that vertex on the given level for the + * multigrid hierarchy. Used by DoFHandler::Policy::ParallelDistributed. + */ + void + fill_level_vertices_with_ghost_neighbors + (const unsigned int level, + std::map > + &vertices_with_ghost_neighbors); + + /** + * This method returns a bit vector of length tria.n_vertices() + * indicating the locally active vertices on a level, i.e., the vertices + * touched by the locally owned level cells for use in geometric + * multigrid (possibly including the vertices due to periodic boundary + * conditions) are marked by true. + * + * Used by DoFHandler::Policy::ParallelDistributed. + */ + std::vector + mark_locally_active_vertices_on_level(const unsigned int level) const; + template friend class dealii::internal::DoFHandler::Policy::ParallelDistributed; }; @@ -1068,6 +1100,23 @@ namespace parallel (std::map > &vertices_with_ghost_neighbors); + /** + * Like above, this method, which is only implemented for dim = 2 or 3, + * needs a stub because it is used in dof_handler_policy.cc + */ + void + fill_level_vertices_with_ghost_neighbors + (const unsigned int level, + std::map > + &vertices_with_ghost_neighbors); + + /** + * Like above, this method, which is only implemented for dim = 2 or 3, + * needs a stub because it is used in dof_handler_policy.cc + */ + std::vector + mark_locally_active_vertices_on_level(const unsigned int level) const; + }; } } diff --git a/source/distributed/tria.cc b/source/distributed/tria.cc index 09a7961a93..5c92934683 100644 --- a/source/distributed/tria.cc +++ b/source/distributed/tria.cc @@ -2585,6 +2585,8 @@ 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 (); @@ -3337,7 +3339,7 @@ namespace parallel // -1. Note that we do not fill other cells we could figure out the // same way, because we might accidentally set an id for a cell that // is not a ghost cell. - for (unsigned int lvl=this->n_levels(); lvl>0;) + for (unsigned int lvl=this->n_levels(); lvl>0; ) { --lvl; for (typename Triangulation::cell_iterator cell = this->begin(lvl); cell!=this->end(lvl); ++cell) @@ -3356,20 +3358,8 @@ namespace parallel //step 2: make sure all the neighbors to our level_cells exist. Need //to look up in p4est... std::vector > marked_vertices(this->n_levels()); - - for (unsigned int lvl=0; lvln_levels(); ++lvl) - { - marked_vertices[lvl].resize(this->n_vertices(), false); - - for (typename dealii::Triangulation::cell_iterator - cell = this->begin(lvl); - cell != this->end(lvl); ++cell) - if (cell->level_subdomain_id() == this->locally_owned_subdomain()) - for (unsigned int v=0; v::vertices_per_cell; ++v) - marked_vertices[lvl][cell->vertex_index(v)] = true; - } - - + for (unsigned int lvl=0; lvl < this->n_levels(); ++lvl) + marked_vertices[lvl] = mark_locally_active_vertices_on_level(lvl); for (typename Triangulation::cell_iterator cell = this->begin(0); cell!=this->end(0); ++cell) { @@ -4195,6 +4185,171 @@ namespace parallel } + + namespace + { + /** + * ensures 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 + */ + template + void + mark_periodic_vertices_recursively(const GridTools::PeriodicFacePair &periodic, + const int target_level, + std::vector &active_vertices_on_level) + { + if (periodic.cell[0]->level() > target_level) + return; + + // for hanging nodes there is nothing to do since we are interested in + // the connections on the same level... + if (periodic.cell[0]->level() < target_level && + periodic.cell[0]->has_children() && + periodic.cell[1]->has_children()) + { + // copy orientations etc. from parent to child + GridTools::PeriodicFacePair periodic_child = periodic; + + // find appropriate pairs of child elements + for (unsigned int cf=0; cfface(periodic.face_idx[0])->n_children(); ++cf) + { + unsigned int c=0; + for (; cn_children(); ++c) + { + if (periodic.cell[0]->child(c)->face(periodic.face_idx[0]) == + periodic.cell[0]->face(periodic.face_idx[0])->child(cf)) + break; + } + Assert(c < periodic.cell[0]->n_children(), + ExcMessage("Face child not found")); + periodic_child.cell[0] = periodic.cell[0]->child(c); + periodic_child.face_idx[0] = periodic.face_idx[0]; + + c=0; + for (; cn_children(); ++c) + { + if (periodic.cell[1]->child(c)->face(periodic.face_idx[1]) == + periodic.cell[1]->face(periodic.face_idx[1])->child(cf)) + break; + } + Assert(c < periodic.cell[1]->n_children(), + ExcMessage("Face child not found")); + periodic_child.cell[1] = periodic.cell[1]->child(c); + periodic_child.face_idx[1] = periodic.face_idx[1]; + + // recursive call into children + mark_periodic_vertices_recursively (periodic_child, target_level, + active_vertices_on_level); + } + return; + } + + if (periodic.cell[0]->level() != target_level) + return; + + // TODO: fix non-standard orientation + Assert(periodic.orientation[0] == true && + periodic.orientation[1] == false && + periodic.orientation[2] == false, + ExcNotImplemented()); + + for (unsigned int v=0; v::vertices_per_cell; ++v) + if (active_vertices_on_level[periodic.cell[0]->face(periodic.face_idx[0])->vertex_index(v)] || + active_vertices_on_level[periodic.cell[1]->face(periodic.face_idx[1])->vertex_index(v)]) + active_vertices_on_level[periodic.cell[0]->face(periodic.face_idx[0])->vertex_index(v)] + = active_vertices_on_level[periodic.cell[1]->face(periodic.face_idx[1])->vertex_index(v)] + = true; + } + + + + /** + * ensures that always both vertices over a periodic face are identified + * together + */ + template + void + set_periodic_ghost_neighbors_recursively(const GridTools::PeriodicFacePair &periodic, + const int target_level, + std::map > &vertices_with_ghost_neighbors) + { + if (periodic.cell[0]->level() > target_level) + return; + + // for hanging nodes there is nothing to do since we are interested in + // the connections on the same level... + if (periodic.cell[0]->level() < target_level && + periodic.cell[0]->has_children() && + periodic.cell[1]->has_children()) + { + // copy orientations etc. from parent to child + GridTools::PeriodicFacePair periodic_child = periodic; + + // find appropriate pairs of child elements + for (unsigned int cf=0; cfface(periodic.face_idx[0])->n_children(); ++cf) + { + unsigned int c=0; + for (; cn_children(); ++c) + { + if (periodic.cell[0]->child(c)->face(periodic.face_idx[0]) == + periodic.cell[0]->face(periodic.face_idx[0])->child(cf)) + break; + } + Assert(c < periodic.cell[0]->n_children(), + ExcMessage("Face child not found")); + periodic_child.cell[0] = periodic.cell[0]->child(c); + periodic_child.face_idx[0] = periodic.face_idx[0]; + + c=0; + for (; cn_children(); ++c) + { + if (periodic.cell[1]->child(c)->face(periodic.face_idx[1]) == + periodic.cell[1]->face(periodic.face_idx[1])->child(cf)) + break; + } + Assert(c < periodic.cell[1]->n_children(), + ExcMessage("Face child not found")); + periodic_child.cell[1] = periodic.cell[1]->child(c); + periodic_child.face_idx[1] = periodic.face_idx[1]; + + // recursive call into children + set_periodic_ghost_neighbors_recursively (periodic_child, target_level, + vertices_with_ghost_neighbors); + } + return; + } + + if (periodic.cell[0]->level() != target_level) + return; + + // TODO: fix non-standard orientation + Assert(periodic.orientation[0] == true && + periodic.orientation[1] == false && + periodic.orientation[2] == false, + ExcNotImplemented()); + + for (unsigned int v=0; v::vertices_per_cell; ++v) + { + const unsigned int + idx0 = periodic.cell[0]->face(periodic.face_idx[0])->vertex_index(v), + idx1 = periodic.cell[1]->face(periodic.face_idx[1])->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()); + } + } + } + + + /** * Determine the neighboring subdomains that are adjacent to each vertex. * This is achieved via the p4est_iterate/p8est_iterate tool @@ -4240,6 +4395,58 @@ namespace parallel } + + /** + * Determine the neighboring subdomains that are adjacent to each vertex + * on the given multigrid level + */ + template + void + Triangulation:: + fill_level_vertices_with_ghost_neighbors + (const unsigned int level, + std::map > + &vertices_with_ghost_neighbors) + { + const std::vector locally_active_vertices = + mark_locally_active_vertices_on_level(level); + for (cell_iterator cell = this->begin(level); cell != this->end(level); ++cell) + if (cell->level_subdomain_id() != dealii::numbers::artificial_subdomain_id + && cell->level_subdomain_id() != this->locally_owned_subdomain()) + for (unsigned int v=0; v::vertices_per_cell; ++v) + if (locally_active_vertices[cell->vertex_index(v)]) + vertices_with_ghost_neighbors[cell->vertex_index(v)] + .insert (cell->level_subdomain_id()); + + for (unsigned int i=0; i + std::vector + Triangulation + ::mark_locally_active_vertices_on_level (const unsigned int level) const + { + Assert (dim>1, ExcNotImplemented()); + + std::vector marked_vertices(this->n_vertices(), false); + for (cell_iterator cell = this->begin(level); cell != this->end(level); ++cell) + if (cell->level_subdomain_id() == this->locally_owned_subdomain()) + for (unsigned int v=0; v::vertices_per_cell; ++v) + marked_vertices[cell->vertex_index(v)] = true; + + for (unsigned int i=0; i void Triangulation::add_periodicity @@ -4377,6 +4584,10 @@ namespace parallel AssertThrow (false, ExcInternalError()); } + periodic_face_pairs_level_0.insert(periodic_face_pairs_level_0.end(), + periodicity_vector.begin(), + periodicity_vector.end()); + #else Assert(false, ExcMessage ("Need p4est version >= 0.3.4.1!")); #endif @@ -4603,6 +4814,26 @@ namespace parallel { Assert (false, ExcNotImplemented()); } + + template + void + Triangulation<1,spacedim>:: + fill_level_vertices_with_ghost_neighbors + (const unsigned int /*level*/, + std::map > + &/*vertices_with_ghost_neighbors*/) + { + Assert (false, ExcNotImplemented()); + } + + template + std::vector + Triangulation<1,spacedim>:: + mark_locally_active_vertices_on_level (const unsigned int) const + { + Assert (false, ExcNotImplemented()); + return std::vector(); + } } } diff --git a/source/dofs/dof_handler_policy.cc b/source/dofs/dof_handler_policy.cc index 7dd0650180..dda7e89b4f 100644 --- a/source/dofs/dof_handler_policy.cc +++ b/source/dofs/dof_handler_policy.cc @@ -1658,52 +1658,6 @@ namespace internal } - /** - * Return a vector in which, - * for every vertex index, we - * mark whether a locally owned - * cell is adjacent. - */ - template - std::vector - mark_locally_active_vertices (const parallel::distributed::Triangulation &triangulation) - { - std::vector locally_active_vertices (triangulation.n_vertices(), - false); - - for (typename dealii::Triangulation::active_cell_iterator - cell = triangulation.begin_active(); - cell != triangulation.end(); ++cell) - if (cell->subdomain_id() == triangulation.locally_owned_subdomain()) - for (unsigned int v=0; v::vertices_per_cell; ++v) - locally_active_vertices[cell->vertex_index(v)] = true; - - return locally_active_vertices; - } - - /** - * Return a vector in which, - * for every vertex index, we - * mark whether a locally owned - * cell is adjacent. - */ - template - std::vector - mark_locally_active_vertices_on_level (const parallel::distributed::Triangulation &triangulation, unsigned int level) - { - std::vector locally_active_vertices (triangulation.n_vertices(), - false); - - if (level::cell_iterator - cell = triangulation.begin(level); - cell != triangulation.end(level); ++cell) - if (cell->level_subdomain_id() == triangulation.locally_owned_subdomain()) - for (unsigned int v=0; v::vertices_per_cell; ++v) - locally_active_vertices[cell->vertex_index(v)] = true; - - return locally_active_vertices; - } template @@ -1763,7 +1717,6 @@ namespace internal } - //sending std::vector > sendbuffers (needs_to_get_cells.size()); std::vector >::iterator buffer = sendbuffers.begin(); @@ -1997,7 +1950,6 @@ namespace internal } - //sending std::vector > sendbuffers (needs_to_get_cells.size()); std::vector >::iterator buffer = sendbuffers.begin(); @@ -2397,6 +2349,7 @@ namespace internal } + template void ParallelDistributed:: @@ -2571,34 +2524,21 @@ namespace internal cell->set_user_flag(); } - //mark the vertices we are interested - //in, i.e. belonging to own and marked cells + //mark the vertices we are interested in, i.e. belonging to own + //and marked cells const std::vector locally_active_vertices - = mark_locally_active_vertices_on_level (*tr, level); + = tr->mark_locally_active_vertices_on_level (level); - // add each ghostcells' - // subdomain to the vertex and - // keep track of interesting - // neighbors + // add each ghostcells' subdomain to the vertex and keep track of + // interesting neighbors std::map > vertices_with_ghost_neighbors; - if (level < tr->n_levels()) - for (typename DoFHandler::level_cell_iterator - cell = dof_handler.begin(level); - cell != dof_handler.end(level); ++cell) - if (cell->level_subdomain_id() != dealii::numbers::artificial_subdomain_id - && cell->level_subdomain_id() != tr->locally_owned_subdomain()) - for (unsigned int v=0; v::vertices_per_cell; ++v) - if (locally_active_vertices[cell->vertex_index(v)]) - vertices_with_ghost_neighbors[cell->vertex_index(v)] - .insert (cell->level_subdomain_id()); - - - /* Send and receive cells. After this, - only the local cells are marked, - that received new data. This has to - be communicated in a second - communication step. */ + tr->fill_level_vertices_with_ghost_neighbors(level, + vertices_with_ghost_neighbors); + + // Send and receive cells. After this, only the local cells are + // marked, that received new data. This has to be communicated in + // a second communication step. communicate_mg_dof_indices_on_marked_cells( dof_handler, vertices_with_ghost_neighbors, -- 2.39.5