From 65bbc4d7ffe58c88285dfc1851a20772c978ea6d Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Sun, 13 Aug 2017 16:41:14 -0600 Subject: [PATCH] Express the ghost exchange of DoF indices using GridTools::exchange_cell_data_to_ghosts(). --- source/dofs/dof_handler_policy.cc | 466 +++++++++--------------------- 1 file changed, 141 insertions(+), 325 deletions(-) diff --git a/source/dofs/dof_handler_policy.cc b/source/dofs/dof_handler_policy.cc index d772c57e95..e4f88eb398 100644 --- a/source/dofs/dof_handler_policy.cc +++ b/source/dofs/dof_handler_policy.cc @@ -27,6 +27,7 @@ #include #include #include +#include DEAL_II_DISABLE_EXTRA_DIAGNOSTICS #include @@ -2765,94 +2766,6 @@ namespace internal - template - void - get_dof_indices_recursively (const typename parallel::distributed::Triangulation &tria, - const unsigned int tree_index, - const CellIteratorType &dealii_cell, - const typename dealii::internal::p4est::types::quadrant &p4est_cell, - const std::map > &vertices_with_ghost_neighbors, - std::map> &needs_to_get_cell) - { - // see if we have to recurse... - if (dealii_cell->has_children()) - { - typename dealii::internal::p4est::types::quadrant - p4est_child[GeometryInfo::max_children_per_cell]; - internal::p4est::init_quadrant_children(p4est_cell, p4est_child); - - - for (unsigned int c=0; c::max_children_per_cell; ++c) - get_dof_indices_recursively (tria, - tree_index, - dealii_cell->child(c), - p4est_child[c], - vertices_with_ghost_neighbors, - needs_to_get_cell); - } - else - { - // we're at a leaf cell. see if the cell is flagged as - // interesting. note that we have only flagged our own cells - // before - if (dealii_cell->user_flag_set() && !dealii_cell->is_ghost()) - { - Assert (dealii_cell->is_locally_owned(), ExcInternalError()); - - // check each vertex if it is interesting and push - // DoF indices if so - std::set send_to; - for (unsigned int v=0; v::vertices_per_cell; ++v) - { - const std::map >::const_iterator - neighbor_subdomains_of_vertex - = vertices_with_ghost_neighbors.find (dealii_cell->vertex_index(v)); - - if (neighbor_subdomains_of_vertex == - vertices_with_ghost_neighbors.end()) - continue; - - Assert(neighbor_subdomains_of_vertex->second.size()!=0, - ExcInternalError()); - - send_to.insert(neighbor_subdomains_of_vertex->second.begin(), - neighbor_subdomains_of_vertex->second.end()); - } - - if (send_to.size() > 0) - { - // this cell's dof_indices need to be sent to - // someone - std::vector - local_dof_indices (dealii_cell->get_fe().dofs_per_cell); - dealii_cell->get_dof_indices (local_dof_indices); - - for (const auto recipient_subdomain : send_to) - { - // get an iterator to what needs to be sent to - // that subdomain (if already exists), or create - // such an object - typename std::map>::iterator - p - = needs_to_get_cell.insert (std::make_pair(recipient_subdomain, - CellDataTransferBuffer())) - .first; - - p->second.tree_indices.push_back(tree_index); - p->second.quadrants.push_back(p4est_cell); - - p->second.dof_numbers_and_indices.push_back(dealii_cell->get_fe().dofs_per_cell); - p->second.dof_numbers_and_indices.insert(p->second.dof_numbers_and_indices.end(), - local_dof_indices.begin(), - local_dof_indices.end()); - - } - } - } - } - } - - template void @@ -3185,81 +3098,6 @@ namespace internal - template - void - set_dof_indices_recursively (const parallel::distributed::Triangulation &tria, - const typename dealii::internal::p4est::types::quadrant &p4est_cell, - const CellIteratorType &dealii_cell, - const typename dealii::internal::p4est::types::quadrant &quadrant, - dealii::types::global_dof_index *dofs) - { - if (internal::p4est::quadrant_is_equal(p4est_cell, quadrant)) - { - Assert(!dealii_cell->has_children(), ExcInternalError()); - Assert(dealii_cell->is_ghost(), ExcInternalError()); - - // update dof indices of cell - std::vector - dof_indices (dealii_cell->get_fe().dofs_per_cell); - dealii_cell->update_cell_dof_indices_cache(); - dealii_cell->get_dof_indices(dof_indices); - - bool complete = true; - for (unsigned int i=0; iuser_flag_set(), ExcInternalError()); - if (complete) - dealii_cell->clear_user_flag(); - - const_cast(dealii_cell)->set_dof_indices(dof_indices); - } - else if (dealii_cell->has_children() - && - internal::p4est::quadrant_is_ancestor (p4est_cell, quadrant)) - { - typename dealii::internal::p4est::types::quadrant - p4est_child[GeometryInfo::max_children_per_cell]; - internal::p4est::init_quadrant_children(p4est_cell, p4est_child); - - for (unsigned int c=0; c::max_children_per_cell; ++c) - set_dof_indices_recursively (tria, p4est_child[c], - dealii_cell->child(c), - quadrant, dofs); - } - } - - - /** * A function that communicates the DoF indices from that subset of * locally owned cells that have their user indices set to the @@ -3308,191 +3146,166 @@ namespace internal void communicate_dof_indices_on_marked_cells (const DoFHandlerType &dof_handler, - const std::map > &vertices_with_ghost_neighbors, - const std::vector &coarse_cell_to_p4est_tree_permutation, - const std::vector &p4est_tree_to_coarse_cell_permutation) + const std::map > &, + const std::vector &, + const std::vector &) { -#ifndef DEAL_II_WITH_P4EST +#ifndef DEAL_II_WITH_MPI (void)vertices_with_ghost_neighbors; Assert (false, ExcNotImplemented()); #else const unsigned int dim = DoFHandlerType::dimension; const unsigned int spacedim = DoFHandlerType::space_dimension; - const parallel::distributed::Triangulation *triangulation - = (dynamic_cast*> - (&dof_handler.get_triangulation())); - Assert (triangulation != nullptr, ExcInternalError()); - - // first, collect cells and their dof_indices to send to - // interested neighbors. the functions called in the - // following loop only work on cells that are locally owned - // and are adjacent to ghost cells and that are marked via - // the user flag. in the first phase, all locally owned - // cells fall in this category, but we remove some of - // these flags later in this function - typedef - std::map> - cellmap_t; - cellmap_t needs_to_get_cells; - - for (typename DoFHandlerType::level_cell_iterator - cell = dof_handler.begin(0); - cell != dof_handler.end(0); - ++cell) - { - typename dealii::internal::p4est::types::quadrant p4est_coarse_cell; - internal::p4est::init_coarse_quadrant(p4est_coarse_cell); - - get_dof_indices_recursively (*triangulation, - coarse_cell_to_p4est_tree_permutation[cell->index()], - cell, - p4est_coarse_cell, - vertices_with_ghost_neighbors, - needs_to_get_cells); - } - - - // step 1: send data to all owners of ghost cells - std::vector > sendbuffers (needs_to_get_cells.size()); - std::vector send_requests (needs_to_get_cells.size()); - - unsigned int idx = 0; - std::vector >::iterator buffer = sendbuffers.begin(); - for (typename cellmap_t::iterator it=needs_to_get_cells.begin(); - it!=needs_to_get_cells.end(); - ++it, ++buffer, ++idx) - { - Assert(it->second.tree_indices.size() > 0, - ExcInternalError()); - Assert(it->second.tree_indices.size() == it->second.quadrants.size(), - ExcInternalError()); - - // pack all the data into the buffer for this recipient - // and send it. keep data around till we can make sure - // that the packet has been received - *buffer = it->second.pack_data (); - const int ierr = MPI_Isend(&(*buffer)[0], buffer->size(), - MPI_BYTE, it->first, - 123, triangulation->get_communicator(), - &send_requests[idx]); - AssertThrowMPI(ierr); - } - - - // mark all of our own cells that miss some dof_data and collect - // the neighbors that are going to send stuff to us - std::set senders; + // define functions that pack data on cells that are ghost cells + // somewhere else, and unpack data on cells where we get information + // from elsewhere + auto pack + = [] (const typename DoFHandlerType::active_cell_iterator &cell) -> std::vector { - std::vector local_dof_indices; - typename DoFHandlerType::active_cell_iterator - cell, endc = dof_handler.end(); + Assert (cell->is_locally_owned(), ExcInternalError()); - for (cell = dof_handler.begin_active(); cell != endc; ++cell) - if (!cell->is_artificial()) - { - if (cell->is_ghost()) - { - if (cell->user_flag_set()) - senders.insert(cell->subdomain_id()); - } - else - { - // this is a locally owned cell. - // - // if we are in phase 1: if it is - // complete already, then we have sent its information - // to the other processors that need it above, and we don't - // need to do so again in phase 2. if it wasn't - // complete, then it will be completed below, and we - // will have to send the completed information a - // second time. - // - // if we are in phase 2, nothing we do matters, so - // we just run this loop anyway - local_dof_indices.resize (cell->get_fe().dofs_per_cell); - cell->get_dof_indices (local_dof_indices); - if (local_dof_indices.end() == - std::find (local_dof_indices.begin(), - local_dof_indices.end(), - numbers::invalid_dof_index)) - cell->clear_user_flag(); - } - } - } + // first see whether we need to do anything at all on this cell. + // this is determined by whether the user_flag is set on the + // cell that indicates whether the *complete* set of DoF indices + // had already been sent + if (cell->user_flag_set()) + { + // get dof indices for the current cell + std::vector local_dof_indices (cell->get_fe().dofs_per_cell); + cell->get_dof_indices (local_dof_indices); + // now see if there are dof indices that were previously + // unknown. this can only happen in phase 1, and in + // that case we know that the user flag must have been set + // + // in any case, if the cell *is* complete, we do not + // need to send the data any more in the next phase. indicate + // this by removing the user flag + if (std::find (local_dof_indices.begin(), + local_dof_indices.end(), + numbers::invalid_dof_index) + != + local_dof_indices.end()) + { + Assert (cell->user_flag_set(), ExcInternalError()); + } + else + cell->clear_user_flag(); - // step 2: receive ghost cell data - std::vector receive; - for (unsigned int i=0; iget_communicator(), &status); - AssertThrowMPI(ierr); - ierr = MPI_Get_count(&status, MPI_BYTE, &len); - AssertThrowMPI(ierr); - receive.resize(len); + return local_dof_indices; + } + else + { + // the fact that the user flag wasn't set means that there is + // nothing we need to send that hasn't been sent so far. + // so return an empty array, but also verify that indeed + // the cell is complete +#ifdef DEBUG + std::vector local_dof_indices (cell->get_fe().dofs_per_cell); + cell->get_dof_indices (local_dof_indices); - char *ptr = &receive[0]; - ierr = MPI_Recv(ptr, len, MPI_BYTE, status.MPI_SOURCE, status.MPI_TAG, - triangulation->get_communicator(), &status); - AssertThrowMPI(ierr); + Assert (std::find (local_dof_indices.begin(), + local_dof_indices.end(), + numbers::invalid_dof_index) + == + local_dof_indices.end(), + ExcInternalError()); +#endif + return std::vector(); + } + }; - CellDataTransferBuffer cell_data_transfer_buffer; - cell_data_transfer_buffer.unpack_data(receive); - unsigned int cells = cell_data_transfer_buffer.tree_indices.size(); - dealii::types::global_dof_index *dofs = cell_data_transfer_buffer.dof_numbers_and_indices.data(); + auto unpack + = [] (const typename DoFHandlerType::active_cell_iterator &cell, + const std::vector &dof_indices) -> void + { + // if we get an empty array, then the sending processor has + // determined that all data has already been sent. that can + // only happen if the complete data set had previously + // been sent, and in that case, (i) the user flag must have been + // deleted, and (ii) the current set of DoF indices must already + // be complete + // + // there is one curious corner case, namely if the cell has no DoFs + // (e.g., when using FE_Nothing), in which case there is nothing + // to do anyway (but we shouldn't fall into the Assert either) + if (dof_indices.empty()) + { + if (cell->get_fe().dofs_per_cell > 0) + Assert (cell->user_flag_set() == false, ExcInternalError()); +#ifdef DEBUG + std::vector local_dof_indices (cell->get_fe().dofs_per_cell); + cell->get_dof_indices (local_dof_indices); - // the dofs pointer contains for each cell the number of - // dofs on that cell (dofs[0]) followed by the dof - // indices themselves. this is repeated for all cells - // in the buffer - for (unsigned int c=0; cuser_flag_set(), ExcInternalError()); + } + else + cell->clear_user_flag(); + + // in any case, set the DoF indices on this cell. some + // of the ones we received may still be invalid because + // the sending processor did not know them yet, so we + // need to merge the ones we get with those that are + // already set here and may have already been known. for + // those that we already know *and* get, they must obviously + // agree + // + // before getting the local dof indices, we need to update the + // cell dof indices cache because we may have set dof indices + // on a neighboring ghost cell before this one, which may have + // affected the dof indices we know about the current cell + std::vector local_dof_indices (cell->get_fe().dofs_per_cell); + cell->update_cell_dof_indices_cache(); + cell->get_dof_indices (local_dof_indices); - typename dealii::internal::p4est::types::quadrant p4est_coarse_cell; - internal::p4est::init_coarse_quadrant(p4est_coarse_cell); + for (unsigned int i=0; iget_fe().dofs_per_cell==dofs[0], ExcInternalError()); + const_cast(cell)->set_dof_indices(local_dof_indices); + } + }; - set_dof_indices_recursively (*triangulation, - p4est_coarse_cell, - cell, - cell_data_transfer_buffer.quadrants[c], - (dofs+1)); - } - } + parallel::GridTools::exchange_cell_data_to_ghosts, DoFHandlerType> + (dof_handler, pack, unpack); // finally update the cell DoF indices caches to make sure // our internal data structures are consistent update_all_active_cell_dof_indices_caches (dof_handler); - // wait for all Isends to complete, so that we can safely destroy the - // buffers. - if (send_requests.size() > 0) - { - const int ierr = MPI_Waitall(send_requests.size(), &send_requests[0], - MPI_STATUSES_IGNORE); - AssertThrowMPI(ierr); - } - - - - // check that all messages got sent and received - Assert (Utilities::MPI::sum (needs_to_get_cells.size(), - triangulation->get_communicator()) - == - Utilities::MPI::sum (senders.size(), - triangulation->get_communicator()), - ExcInternalError()); - // have a barrier so that sends between two calls to this // function are not mixed up. // @@ -3506,6 +3319,9 @@ namespace internal // different tags for phase 1 and 2, but the cost of a // barrier is negligible compared to everything else we do // here + const parallel::distributed::Triangulation< dim, spacedim > *triangulation + = (dynamic_cast*> + (&dof_handler.get_triangulation())); const int ierr = MPI_Barrier(triangulation->get_communicator()); AssertThrowMPI(ierr); #endif -- 2.39.5