From f5b25ea3ea1117daab733a8015dc88b1e5e8b9f4 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Sat, 8 Nov 2014 15:09:48 -0600 Subject: [PATCH] Move communicate_locally_moved_vertices to parallel::distributed::Triangulation. --- include/deal.II/distributed/tria.h | 85 ++++++ include/deal.II/grid/grid_tools.h | 20 -- source/distributed/tria.cc | 452 ++++++++++++++++++++++++++++ source/grid/grid_tools.cc | 460 +---------------------------- source/grid/grid_tools.inst.in | 4 - 5 files changed, 538 insertions(+), 483 deletions(-) diff --git a/include/deal.II/distributed/tria.h b/include/deal.II/distributed/tria.h index f0b515b2ef..4fb06a477e 100644 --- a/include/deal.II/distributed/tria.h +++ b/include/deal.II/distributed/tria.h @@ -468,6 +468,55 @@ namespace parallel */ virtual void execute_coarsening_and_refinement (); + /** + * When vertices have been moved locally, for example using code + * like + * @code + * cell->vertex(0) = new_location; + * @endcode + * then this function can be used to update the location of + * vertices between MPI processes. + * + * All the vertices that have been moved and might be in the ghost layer + * of a process have to be reported in the @p vertex_locally_moved argument. + * This ensures that that part of the information that has to be send + * between processes is actually sent. Additionally, it is quite important that vertices + * on the boundary between processes are reported on exactly one process + * (e.g. the one with the highest id). Otherwise we could expect + * undesirable results if multiple processes move a vertex differently. + * A typical strategy is to let processor $i$ move those vertices that + * are adjacent to cells whose owners include processor $i$ but no + * other processor $j$ with $j &vertex_locally_moved); + /** * Return the subdomain id of those cells that are owned by the * current processor. All cells in the triangulation that do not @@ -958,6 +1007,42 @@ namespace parallel const std::vector & get_p4est_tree_to_coarse_cell_permutation() const; + /** + * When vertices have been moved locally, for example using code + * like + * @code + * cell->vertex(0) = new_location; + * @endcode + * then this function can be used to update the location of + * vertices between MPI processes. + * + * All the vertices that have been moved and might be in the ghost layer + * of a process have to be reported in the @p vertex_locally_moved argument. + * This ensures that that part of the information that has to be send + * between processes is actually sent. Additionally, it is quite important that vertices + * on the boundary between processes are reported on exactly one process + * (e.g. the one with the highest id). Otherwise we could expect + * undesirable results if multiple processes move a vertex differently. + * A typical strategy is to let processor $i$ move those vertices that + * are adjacent to cells whose owners include processor $i$ but no + * other processor $j$ with $j &vertex_locally_moved); + /** * Return the subdomain id of those cells that are owned by the * current processor. All cells in the triangulation that do not diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index 15a06de672..eb00e598fc 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -373,26 +373,6 @@ namespace GridTools Triangulation &triangulation, const bool keep_boundary=true); - /** - * When vertices have been moved locally, this function updates the - * location of the vertices in the ghost layer on all MPI processes. - * - * All the vertices that have been moved and might be in the ghost layer - * of a process have to be reported to @p vertex_locally_moved. - * This ensures that we reduce the information that has to be send - * between processes. Additionally, it is quite important that vertices - * on the boundary between processes are reported on exactly one process - * (e.g. the one with the highest id). Otherwise we could expect - * undesirable results if multiple processes move a vertex differently. - * - * This functions needs not to be called if all MPI processes know - * how the position of all the vertices changed. - */ - template - void - communicate_locally_moved_vertices (parallel::distributed::Triangulation &triangulation, - const std::vector &vertex_locally_moved); - /*@}*/ /** * @name Finding cells and vertices of a triangulation diff --git a/source/distributed/tria.cc b/source/distributed/tria.cc index 07f4a6986b..643e3e3a64 100644 --- a/source/distributed/tria.cc +++ b/source/distributed/tria.cc @@ -2099,6 +2099,309 @@ namespace parallel } + // This anonymous namespace contains utility for + // the function Triangulation::communicate_locally_moved_vertices + namespace CommunicateLocallyMovedVertices + { + namespace + { + /** + * A list of tree+quadrant and their vertex indices. + * The bool vector describes which vertices are of interest + * and should be set on the receiving processes. + */ + template + struct CellInfo + { + // store all the tree_indices we send/receive consecutively (n_cells entries) + std::vector tree_index; + // store all the quadrants we send/receive consecutively (n_cells entries) + std::vector::quadrant> quadrants; + // store for each cell the number of vertices we send/receive + // and then the vertex indices (for each cell: n_vertices+1 entries) + std::vector vertex_indices; + // store for each cell the vertices we send/receive + // (for each cell n_vertices entries) + std::vector > vertices; + // for receiving and unpacking data we need to store pointers to the + // first vertex and vertex_index on each cell additionally + // both vectors have as many entries as there are cells + std::vector first_vertex_indices; + std::vector* > first_vertices; + + unsigned int bytes_for_buffer () const + { + return (sizeof(unsigned int) + + tree_index.size() * sizeof(unsigned int) + + quadrants.size() * sizeof(typename dealii::internal::p4est + ::types::quadrant) + + vertices.size() * sizeof(dealii::Point)) + + vertex_indices.size() * sizeof(unsigned int); + } + + void pack_data (std::vector &buffer) const + { + buffer.resize(bytes_for_buffer()); + + char *ptr = &buffer[0]; + + const unsigned int num_cells = tree_index.size(); + std::memcpy(ptr, &num_cells, sizeof(unsigned int)); + ptr += sizeof(unsigned int); + + std::memcpy(ptr, + &tree_index[0], + num_cells*sizeof(unsigned int)); + ptr += num_cells*sizeof(unsigned int); + + std::memcpy(ptr, + &quadrants[0], + num_cells * sizeof(typename dealii::internal::p4est:: + types::quadrant)); + ptr += num_cells*sizeof(typename dealii::internal::p4est::types:: + quadrant); + + std::memcpy(ptr, + &vertex_indices[0], + vertex_indices.size() * sizeof(unsigned int)); + ptr += vertex_indices.size() * sizeof(unsigned int); + + std::memcpy(ptr, + &vertices[0], + vertices.size() * sizeof(dealii::Point)); + ptr += vertices.size() * sizeof(dealii::Point); + + Assert (ptr == &buffer[0]+buffer.size(), + ExcInternalError()); + + } + + void unpack_data (const std::vector &buffer) + { + const char *ptr = &buffer[0]; + unsigned int cells; + memcpy(&cells, ptr, sizeof(unsigned int)); + ptr += sizeof(unsigned int); + + tree_index.resize(cells); + memcpy(&tree_index[0],ptr,sizeof(unsigned int)*cells); + ptr += sizeof(unsigned int)*cells; + + quadrants.resize(cells); + memcpy(&quadrants[0],ptr, + sizeof(typename dealii::internal::p4est::types::quadrant)*cells); + ptr += sizeof(typename dealii::internal::p4est::types::quadrant)*cells; + + vertex_indices.clear(); + first_vertex_indices.resize(cells); + std::vector n_vertices_on_cell(cells); + std::vector first_indices (cells); + for (unsigned int c=0; c(ptr); + first_indices[c] = vertex_indices.size(); + vertex_indices.push_back(*vertex_index); + n_vertices_on_cell[c] = *vertex_index; + ptr += sizeof(unsigned int); + // Now copy all the 'real' vertex_indices + vertex_indices.resize(vertex_indices.size() + n_vertices_on_cell[c]); + memcpy(&vertex_indices[vertex_indices.size() - n_vertices_on_cell[c]], + ptr, n_vertices_on_cell[c]*sizeof(unsigned int)); + ptr += n_vertices_on_cell[c]*sizeof(unsigned int); + } + for (unsigned int c=0; c *const vertex + = reinterpret_cast * const>(ptr); + first_indices[c] = vertices.size(); + vertices.push_back(*vertex); + ptr += sizeof(dealii::Point); + vertices.resize(vertices.size() + n_vertices_on_cell[c]-1); + memcpy(&vertices[vertices.size() - (n_vertices_on_cell[c]-1)], + ptr, (n_vertices_on_cell[c]-1)*sizeof(dealii::Point)); + ptr += (n_vertices_on_cell[c]-1)*sizeof(dealii::Point); + } + for (unsigned int c=0; c + void + fill_vertices_recursively (const typename parallel::distributed::Triangulation &tria, + const unsigned int tree_index, + const typename Triangulation::cell_iterator &dealii_cell, + const typename dealii::internal::p4est::types::quadrant &p4est_cell, + const std::map > &vertices_with_ghost_neighbors, + const std::vector &vertex_locally_moved, + 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]; + dealii::internal::p4est::init_quadrant_children(p4est_cell, p4est_child); + + + for (unsigned int c=0; c::max_children_per_cell; ++c) + fill_vertices_recursively(tria, + tree_index, + dealii_cell->child(c), + p4est_child[c], + vertices_with_ghost_neighbors, + vertex_locally_moved, + needs_to_get_cell); + return; + } + + // We're at a leaf cell. + // If one of the vertices of the cell is interesting, + // sent all moved vertices of the cell. + if (dealii_cell->is_locally_owned()) + { + 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()) + { + 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) + { + std::vector vertex_indices; + std::vector > local_vertices; + for (unsigned int v=0; v::vertices_per_cell; ++v) + if (vertex_locally_moved[dealii_cell->vertex_index(v)]) + { + vertex_indices.push_back(v); + local_vertices.push_back(dealii_cell->vertex(v)); + } + + for (std::set::iterator it=send_to.begin(); + it!=send_to.end(); ++it) + { + const dealii::types::subdomain_id subdomain = *it; + + // 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(subdomain, + CellInfo())) + .first; + + p->second.tree_index.push_back(tree_index); + p->second.quadrants.push_back(p4est_cell); + + p->second.vertex_indices.push_back(vertex_indices.size()); + p->second.vertex_indices.insert(p->second.vertex_indices.end(), + vertex_indices.begin(), + vertex_indices.end()); + + p->second.vertices.insert(p->second.vertices.end(), + local_vertices.begin(), + local_vertices.end()); + } + } + } + } + + + + // After the cell data has been received this function is responsible + // for moving the vertices in the corresponding ghost layer locally. + // As in fill_vertices_recursively for each dealii cell on the + // coarsest level the corresponding p4est_cell has to be provided + // when calling this function. By recursing through through all + // children we consider each active cell. + // Additionally, we need to give a pointer to the first vertex indices + // and vertices. Since the first information saved in vertex_indices + // is the number of vertices this all the information we need. + template + void + set_vertices_recursively ( + const parallel::distributed::Triangulation &tria, + const typename dealii::internal::p4est::types::quadrant &p4est_cell, + const typename Triangulation::cell_iterator &dealii_cell, + const typename dealii::internal::p4est::types::quadrant &quadrant, + const dealii::Point *const vertices, + const unsigned int *const vertex_indices) + { + if (dealii::internal::p4est::quadrant_is_equal(p4est_cell, quadrant)) + { + Assert(!dealii_cell->is_artificial(), ExcInternalError()); + Assert(!dealii_cell->has_children(), ExcInternalError()); + Assert(!dealii_cell->is_locally_owned(), ExcInternalError()); + + const unsigned int n_vertices = vertex_indices[0]; + + // update dof indices of cell + for (unsigned int i=0; ivertex(vertex_indices[i+1]) = vertices[i]; + + return; + } + + if (! dealii_cell->has_children()) + return; + + if (! dealii::internal::p4est::quadrant_is_ancestor (p4est_cell, quadrant)) + return; + + typename dealii::internal::p4est::types::quadrant + p4est_child[GeometryInfo::max_children_per_cell]; + dealii::internal::p4est::init_quadrant_children(p4est_cell, p4est_child); + + for (unsigned int c=0; c::max_children_per_cell; ++c) + set_vertices_recursively (tria, p4est_child[c], + dealii_cell->child(c), + quadrant, vertices, + vertex_indices); + } + } + } + + template void @@ -3182,6 +3485,146 @@ namespace parallel + + + + template + void + Triangulation:: + communicate_locally_moved_vertices (const std::vector &vertex_locally_moved) + { + Assert (vertex_locally_moved.size() == this->n_vertices(), + ExcDimensionMismatch(vertex_locally_moved.size(), + this->n_vertices())); + + std::map > + vertices_with_ghost_neighbors; + + // First find out which process should receive which vertices... + for (typename Triangulation::active_cell_iterator + cell=this->begin_active(); cell!=this->end(); + ++cell) + if (cell->is_ghost()) + for (unsigned int vertex_no=0; + vertex_no::vertices_per_cell; ++vertex_no) + { + const unsigned int global_vertex_no = cell->vertex_index(vertex_no); + vertices_with_ghost_neighbors[global_vertex_no].insert + (cell->subdomain_id()); + } + + // now collect cells and their vertices + // for the interested neighbors + typedef + std::map > cellmap_t; + cellmap_t needs_to_get_cells; + + for (typename Triangulation::cell_iterator + cell = this->begin(0); + cell != this->end(0); + ++cell) + { + typename dealii::internal::p4est::types::quadrant p4est_coarse_cell; + dealii::internal::p4est::init_coarse_quadrant(p4est_coarse_cell); + + CommunicateLocallyMovedVertices::fill_vertices_recursively + (*this, + this->get_coarse_cell_to_p4est_tree_permutation()[cell->index()], + cell, + p4est_coarse_cell, + vertices_with_ghost_neighbors, + vertex_locally_moved, + needs_to_get_cells); + } + + // sending + std::vector > sendbuffers (needs_to_get_cells.size()); + std::vector >::iterator buffer = sendbuffers.begin(); + std::vector requests (needs_to_get_cells.size()); + std::vector destinations; + + unsigned int idx=0; + + for (typename cellmap_t::iterator it=needs_to_get_cells.begin(); + it!=needs_to_get_cells.end(); + ++it, ++buffer, ++idx) + { + const unsigned int num_cells = it->second.tree_index.size(); + destinations.push_back(it->first); + + Assert(num_cells==it->second.quadrants.size(), ExcInternalError()); + Assert(num_cells>0, 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 + it->second.pack_data (*buffer); + MPI_Isend(&(*buffer)[0], buffer->size(), + MPI_BYTE, it->first, + 123, this->get_communicator(), &requests[idx]); + } + + Assert(destinations.size()==needs_to_get_cells.size(), ExcInternalError()); + + // collect the neighbors + // that are going to send stuff to us + const std::vector senders + = Utilities::MPI::compute_point_to_point_communication_pattern + (this->get_communicator(), destinations); + + // receive ghostcelldata + std::vector receive; + CommunicateLocallyMovedVertices::CellInfo cellinfo; + for (unsigned int i=0; iget_communicator(), &status); + MPI_Get_count(&status, MPI_BYTE, &len); + receive.resize(len); + + char *ptr = &receive[0]; + MPI_Recv(ptr, len, MPI_BYTE, status.MPI_SOURCE, status.MPI_TAG, + this->get_communicator(), &status); + + cellinfo.unpack_data(receive); + const unsigned int cells = cellinfo.tree_index.size(); + for (unsigned int c=0; c::cell_iterator + cell (this, + 0, + this->get_p4est_tree_to_coarse_cell_permutation()[cellinfo.tree_index[c]]); + + typename dealii::internal::p4est::types::quadrant p4est_coarse_cell; + dealii::internal::p4est::init_coarse_quadrant(p4est_coarse_cell); + + CommunicateLocallyMovedVertices::set_vertices_recursively (*this, + p4est_coarse_cell, + cell, + cellinfo.quadrants[c], + cellinfo.first_vertices[c], + cellinfo.first_vertex_indices[c]); + } + } + + // complete all sends, so that we can + // safely destroy the buffers. + if (requests.size() > 0) + MPI_Waitall(requests.size(), &requests[0], MPI_STATUSES_IGNORE); + + //check all msgs got sent and received + Assert(Utilities::MPI::sum(needs_to_get_cells.size(), this->get_communicator()) + == Utilities::MPI::sum(senders.size(), this->get_communicator()), + ExcInternalError()); + } + + + template types::subdomain_id Triangulation::locally_owned_subdomain () const @@ -3911,6 +4354,15 @@ namespace parallel + template + void + Triangulation<1,spacedim>::communicate_locally_moved_vertices + (const std::vector &vertex_locally_moved) + { + Assert (false, ExcNotImplemented()); + } + + template types::subdomain_id Triangulation<1,spacedim>::locally_owned_subdomain () const diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index d66d12b7f0..b736aadf55 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -95,308 +95,6 @@ namespace GridTools } -#ifdef DEAL_II_WITH_P4EST - // This anonymous namespace contains utility for - // the function GridTools::communicate_locally_moved_vertices - namespace - { - /** - * A list of tree+quadrant and their vertex indices. - * The bool vector describes which vertices are of interest - * and should be set on the receiving processes. - */ - template - struct CellInfo - { - // store all the tree_indices we send/receive consecutively (n_cells entries) - std::vector tree_index; - // store all the quadrants we send/receive consecutively (n_cells entries) - std::vector::quadrant> quadrants; - // store for each cell the number of vertices we send/receive - // and then the vertex indices (for each cell: n_vertices+1 entries) - std::vector vertex_indices; - // store for each cell the vertices we send/receive - // (for each cell n_vertices entries) - std::vector > vertices; - // for receiving and unpacking data we need to store pointers to the - // first vertex and vertex_index on each cell additionally - // both vectors have as many entries as there are cells - std::vector first_vertex_indices; - std::vector* > first_vertices; - - unsigned int bytes_for_buffer () const - { - return (sizeof(unsigned int) + - tree_index.size() * sizeof(unsigned int) + - quadrants.size() * sizeof(typename dealii::internal::p4est - ::types::quadrant) + - vertices.size() * sizeof(dealii::Point)) + - vertex_indices.size() * sizeof(unsigned int); - } - - void pack_data (std::vector &buffer) const - { - buffer.resize(bytes_for_buffer()); - - char *ptr = &buffer[0]; - - const unsigned int num_cells = tree_index.size(); - std::memcpy(ptr, &num_cells, sizeof(unsigned int)); - ptr += sizeof(unsigned int); - - std::memcpy(ptr, - &tree_index[0], - num_cells*sizeof(unsigned int)); - ptr += num_cells*sizeof(unsigned int); - - std::memcpy(ptr, - &quadrants[0], - num_cells * sizeof(typename dealii::internal::p4est:: - types::quadrant)); - ptr += num_cells*sizeof(typename dealii::internal::p4est::types:: - quadrant); - - std::memcpy(ptr, - &vertex_indices[0], - vertex_indices.size() * sizeof(unsigned int)); - ptr += vertex_indices.size() * sizeof(unsigned int); - - std::memcpy(ptr, - &vertices[0], - vertices.size() * sizeof(dealii::Point)); - ptr += vertices.size() * sizeof(dealii::Point); - - Assert (ptr == &buffer[0]+buffer.size(), - ExcInternalError()); - - } - - void unpack_data (const std::vector &buffer) - { - const char *ptr = &buffer[0]; - unsigned int cells; - memcpy(&cells, ptr, sizeof(unsigned int)); - ptr += sizeof(unsigned int); - - tree_index.resize(cells); - memcpy(&tree_index[0],ptr,sizeof(unsigned int)*cells); - ptr += sizeof(unsigned int)*cells; - - quadrants.resize(cells); - memcpy(&quadrants[0],ptr, - sizeof(typename dealii::internal::p4est::types::quadrant)*cells); - ptr += sizeof(typename dealii::internal::p4est::types::quadrant)*cells; - - vertex_indices.clear(); - first_vertex_indices.resize(cells); - std::vector n_vertices_on_cell(cells); - std::vector first_indices (cells); - for (unsigned int c=0; c(ptr); - first_indices[c] = vertex_indices.size(); - vertex_indices.push_back(*vertex_index); - n_vertices_on_cell[c] = *vertex_index; - ptr += sizeof(unsigned int); - // Now copy all the 'real' vertex_indices - vertex_indices.resize(vertex_indices.size() + n_vertices_on_cell[c]); - memcpy(&vertex_indices[vertex_indices.size() - n_vertices_on_cell[c]], - ptr, n_vertices_on_cell[c]*sizeof(unsigned int)); - ptr += n_vertices_on_cell[c]*sizeof(unsigned int); - } - for (unsigned int c=0; c *const vertex - = reinterpret_cast * const>(ptr); - first_indices[c] = vertices.size(); - vertices.push_back(*vertex); - ptr += sizeof(dealii::Point); - vertices.resize(vertices.size() + n_vertices_on_cell[c]-1); - memcpy(&vertices[vertices.size() - (n_vertices_on_cell[c]-1)], - ptr, (n_vertices_on_cell[c]-1)*sizeof(dealii::Point)); - ptr += (n_vertices_on_cell[c]-1)*sizeof(dealii::Point); - } - for (unsigned int c=0; c - void - fill_vertices_recursively (const typename parallel::distributed::Triangulation &tria, - const unsigned int tree_index, - const typename Triangulation::cell_iterator &dealii_cell, - const typename dealii::internal::p4est::types::quadrant &p4est_cell, - const std::map > &vertices_with_ghost_neighbors, - const std::vector &vertex_locally_moved, - 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) - fill_vertices_recursively(tria, - tree_index, - dealii_cell->child(c), - p4est_child[c], - vertices_with_ghost_neighbors, - vertex_locally_moved, - needs_to_get_cell); - return; - } - - // We're at a leaf cell. - // If one of the vertices of the cell is interesting, - // sent all moved vertices of the cell. - if (dealii_cell->is_locally_owned()) - { - 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()) - { - 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) - { - std::vector vertex_indices; - std::vector > local_vertices; - for (unsigned int v=0; v::vertices_per_cell; ++v) - if (vertex_locally_moved[dealii_cell->vertex_index(v)]) - { - vertex_indices.push_back(v); - local_vertices.push_back(dealii_cell->vertex(v)); - } - - for (std::set::iterator it=send_to.begin(); - it!=send_to.end(); ++it) - { - const dealii::types::subdomain_id subdomain = *it; - - // 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(subdomain, - CellInfo())) - .first; - - p->second.tree_index.push_back(tree_index); - p->second.quadrants.push_back(p4est_cell); - - p->second.vertex_indices.push_back(vertex_indices.size()); - p->second.vertex_indices.insert(p->second.vertex_indices.end(), - vertex_indices.begin(), - vertex_indices.end()); - - p->second.vertices.insert(p->second.vertices.end(), - local_vertices.begin(), - local_vertices.end()); - } - } - } - } - - - - // After the cell data has been received this function is responsible - // for moving the vertices in the corresponding ghost layer locally. - // As in fill_vertices_recursively for each dealii cell on the - // coarsest level the corresponding p4est_cell has to be provided - // when calling this function. By recursing through through all - // children we consider each active cell. - // Additionally, we need to give a pointer to the first vertex indices - // and vertices. Since the first information saved in vertex_indices - // is the number of vertices this all the information we need. - template - void - set_vertices_recursively ( - const parallel::distributed::Triangulation &tria, - const typename dealii::internal::p4est::types::quadrant &p4est_cell, - const typename Triangulation::cell_iterator &dealii_cell, - const typename dealii::internal::p4est::types::quadrant &quadrant, - const dealii::Point *const vertices, - const unsigned int *const vertex_indices) - { - if (internal::p4est::quadrant_is_equal(p4est_cell, quadrant)) - { - Assert(!dealii_cell->is_artificial(), ExcInternalError()); - Assert(!dealii_cell->has_children(), ExcInternalError()); - Assert(!dealii_cell->is_locally_owned(), ExcInternalError()); - - const unsigned int n_vertices = vertex_indices[0]; - - // update dof indices of cell - for (unsigned int i=0; ivertex(vertex_indices[i+1]) = vertices[i]; - - return; - } - - if (! dealii_cell->has_children()) - return; - - if (! internal::p4est::quadrant_is_ancestor (p4est_cell, quadrant)) - return; - - 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_vertices_recursively (tria, p4est_child[c], - dealii_cell->child(c), - quadrant, vertices, - vertex_indices); - } - } -#endif //DEAL_II_WITH_P4EST - - template double @@ -937,162 +635,6 @@ namespace GridTools - template - void - communicate_locally_moved_vertices - (parallel::distributed::Triangulation<1,spacedim> &triangulation, - const std::vector &vertex_locally_moved) - { - Assert (false, ExcNotImplemented()); - } - - - - template - void - communicate_locally_moved_vertices - (parallel::distributed::Triangulation &triangulation, - const std::vector &vertex_locally_moved) - { -#ifndef DEAL_II_WITH_P4EST - (void)vertex_locally_moved; - Assert (false, ExcNotImplemented()); -#else - - Assert (vertex_locally_moved.size() == triangulation.n_vertices(), - ExcDimensionMismatch(vertex_locally_moved.size(), - triangulation.n_vertices())); - - std::map > - vertices_with_ghost_neighbors; - - // First find out which process should receive which vertices... - for (typename Triangulation::active_cell_iterator - cell=triangulation.begin_active(); cell!=triangulation.end(); - ++cell) - if (cell->is_ghost()) - for (unsigned int vertex_no=0; - vertex_no::vertices_per_cell; ++vertex_no) - { - const unsigned int global_vertex_no = cell->vertex_index(vertex_no); - vertices_with_ghost_neighbors[global_vertex_no].insert - (cell->subdomain_id()); - } - - // now collect cells and their vertices - // for the interested neighbors - typedef - std::map > cellmap_t; - cellmap_t needs_to_get_cells; - - for (typename Triangulation::cell_iterator - cell = triangulation.begin(0); - cell != triangulation.end(0); - ++cell) - { - typename dealii::internal::p4est::types::quadrant p4est_coarse_cell; - internal::p4est::init_coarse_quadrant(p4est_coarse_cell); - - fill_vertices_recursively - (triangulation, - triangulation.get_coarse_cell_to_p4est_tree_permutation()[cell->index()], - cell, - p4est_coarse_cell, - vertices_with_ghost_neighbors, - vertex_locally_moved, - needs_to_get_cells); - } - - //sending - std::vector > sendbuffers (needs_to_get_cells.size()); - std::vector >::iterator buffer = sendbuffers.begin(); - std::vector requests (needs_to_get_cells.size()); - std::vector destinations; - - unsigned int idx=0; - - for (typename cellmap_t::iterator it=needs_to_get_cells.begin(); - it!=needs_to_get_cells.end(); - ++it, ++buffer, ++idx) - { - const unsigned int num_cells = it->second.tree_index.size(); - destinations.push_back(it->first); - - Assert(num_cells==it->second.quadrants.size(), ExcInternalError()); - Assert(num_cells>0, 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 - it->second.pack_data (*buffer); - MPI_Isend(&(*buffer)[0], buffer->size(), - MPI_BYTE, it->first, - 123, triangulation.get_communicator(), &requests[idx]); - } - - Assert(destinations.size()==needs_to_get_cells.size(), ExcInternalError()); - - // collect the neighbors - // that are going to send stuff to us - const std::vector senders - = Utilities::MPI::compute_point_to_point_communication_pattern - (triangulation.get_communicator(), destinations); - - // receive ghostcelldata - std::vector receive; - CellInfo cellinfo; - for (unsigned int i=0; i::cell_iterator - cell (&triangulation, - 0, - triangulation.get_p4est_tree_to_coarse_cell_permutation()[cellinfo.tree_index[c]]); - - typename dealii::internal::p4est::types::quadrant p4est_coarse_cell; - internal::p4est::init_coarse_quadrant(p4est_coarse_cell); - - set_vertices_recursively (triangulation, - p4est_coarse_cell, - cell, - cellinfo.quadrants[c], - cellinfo.first_vertices[c], - cellinfo.first_vertex_indices[c]); - } - } - - // complete all sends, so that we can - // safely destroy the buffers. - if (requests.size() > 0) - MPI_Waitall(requests.size(), &requests[0], MPI_STATUSES_IGNORE); - - //check all msgs 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()); - -#endif //DEAL_II_WITH_P4EST - } - - - /** * Distort a triangulation in * some random way. @@ -1224,7 +766,7 @@ namespace GridTools } if (distributed) - communicate_locally_moved_vertices(*tr, vertex_moved); + tr->communicate_locally_moved_vertices(vertex_moved); // Correct hanging nodes if necessary diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in index 91e6ef835c..a11c8cbf72 100644 --- a/source/grid/grid_tools.inst.in +++ b/source/grid/grid_tools.inst.in @@ -185,10 +185,6 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS Triangulation &, const bool); - template - void communicate_locally_moved_vertices (parallel::distributed::Triangulation &, - const std::vector &); - template void get_face_connectivity_of_cells (const Triangulation &triangulation, -- 2.39.5