*/
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<i$; in other words, for vertices
+ * at the boundary of a subdomain, the processor with the lowest subdomain
+ * id "owns" a vertex.
+ *
+ * @note It only makes sense to move vertices that are either located
+ * on locally owned cells or on cells in the ghost layer. This is
+ * because you can be sure that these vertices indeed exist on the
+ * finest mesh aggregated over all processors, whereas vertices on
+ * artificial cells but not at least in the ghost layer may or may
+ * not exist on the globally finest mesh. Consequently, the
+ * @p vertex_locally_moved argument may not contain vertices that
+ * aren't at least on ghost cells.
+ *
+ * @note This function moves vertices in such a way that on every processor,
+ * the vertices of every locally owned and ghost cell is consistent with
+ * the corresponding location of these cells on other processors. On the
+ * other hand, the locations of artificial cells will in general be wrong
+ * since artificial cells may or may not exist on other processors
+ * and consequently it is not possible to determine their location
+ * in any way. This is not usually a problem since one never does anything
+ * on artificial cells. However, it may lead to problems if the mesh
+ * with moved vertices is refined in a later step. If that's what you
+ * want to do, the right way to do it is to save the offset applied to
+ * every vertex, call this function, and before refining or coarsening
+ * the mesh apply the opposite offset and call this function again.
+ *
+ * @see This function is used, for example, in GridTools::distort_random().
+ */
+ void
+ communicate_locally_moved_vertices (const std::vector<bool> &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
const std::vector<types::global_dof_index> &
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<i$; in other words, for vertices
+ * at the boundary of a subdomain, the processor with the lowest subdomain
+ * id "owns" a vertex.
+ *
+ * @note It only makes sense to move vertices that are either located
+ * on locally owned cells or on cells in the ghost layer. This is
+ * because you can be sure that these vertices indeed exist on the
+ * finest mesh aggregated over all processors, whereas vertices on
+ * artificial cells but not at least in the ghost layer may or may
+ * not exist on the globally finest mesh. Consequently, the
+ * @p vertex_locally_moved argument may not contain vertices that
+ * aren't at least on ghost cells.
+ *
+ * @see This function is used, for example, in GridTools::distort_random().
+ */
+ void
+ communicate_locally_moved_vertices (const std::vector<bool> &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
Triangulation<dim, spacedim> &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 <int dim, int spacedim>
- void
- communicate_locally_moved_vertices (parallel::distributed::Triangulation<dim, spacedim> &triangulation,
- const std::vector<bool> &vertex_locally_moved);
-
/*@}*/
/**
* @name Finding cells and vertices of a triangulation
}
+ // 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 <int dim, int spacedim>
+ struct CellInfo
+ {
+ // store all the tree_indices we send/receive consecutively (n_cells entries)
+ std::vector<unsigned int> tree_index;
+ // store all the quadrants we send/receive consecutively (n_cells entries)
+ std::vector<typename dealii::internal::p4est::types<dim>::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<unsigned int> vertex_indices;
+ // store for each cell the vertices we send/receive
+ // (for each cell n_vertices entries)
+ std::vector<dealii::Point<spacedim> > 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<unsigned int * > first_vertex_indices;
+ std::vector<dealii::Point<spacedim>* > 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<dim>::quadrant) +
+ vertices.size() * sizeof(dealii::Point<spacedim>)) +
+ vertex_indices.size() * sizeof(unsigned int);
+ }
+
+ void pack_data (std::vector<char> &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<dim>::quadrant));
+ ptr += num_cells*sizeof(typename dealii::internal::p4est::types<dim>::
+ 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<spacedim>));
+ ptr += vertices.size() * sizeof(dealii::Point<spacedim>);
+
+ Assert (ptr == &buffer[0]+buffer.size(),
+ ExcInternalError());
+
+ }
+
+ void unpack_data (const std::vector<char> &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<dim>::quadrant)*cells);
+ ptr += sizeof(typename dealii::internal::p4est::types<dim>::quadrant)*cells;
+
+ vertex_indices.clear();
+ first_vertex_indices.resize(cells);
+ std::vector<unsigned int> n_vertices_on_cell(cells);
+ std::vector<unsigned int> first_indices (cells);
+ for (unsigned int c=0; c<cells; ++c)
+ {
+ // The first 'vertex index' is the number of vertices.
+ // Additionally, we need to store the pointer to this
+ // vertex index with respect to the std::vector
+ const unsigned int *const vertex_index
+ = reinterpret_cast<const unsigned int *const>(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<cells; ++c)
+ first_vertex_indices[c] = &vertex_indices[first_indices[c]];
+
+ vertices.clear();
+ first_vertices.resize(cells);
+ for (unsigned int c=0; c<cells; ++c)
+ {
+ // We need to store a pointer to the first vertex.
+ const dealii::Point<spacedim> *const vertex
+ = reinterpret_cast<const dealii::Point<spacedim> * const>(ptr);
+ first_indices[c] = vertices.size();
+ vertices.push_back(*vertex);
+ ptr += sizeof(dealii::Point<spacedim>);
+ 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<spacedim>));
+ ptr += (n_vertices_on_cell[c]-1)*sizeof(dealii::Point<spacedim>);
+ }
+ for (unsigned int c=0; c<cells; ++c)
+ first_vertices[c] = &vertices[first_indices[c]];
+
+ Assert (ptr == &buffer[0]+buffer.size(),
+ ExcInternalError());
+ }
+ };
+
+
+
+ // This function is responsible for gathering the information
+ // we want to send to each process.
+ // For each dealii cell on the coarsest level the corresponding
+ // p4est_cell has to be provided when calling this function.
+ // By recursing through all children we consider each active cell.
+ // vertices_with_ghost_neighbors tells us which vertices
+ // are in the ghost layer and for which processes they might
+ // be interesting.
+ // Whether a vertex has actually been updated locally is
+ // stored in vertex_locally_moved. Only those are considered
+ // for sending.
+ // The gathered information is saved into needs_to_get_cell.
+ template <int dim, int spacedim>
+ void
+ fill_vertices_recursively (const typename parallel::distributed::Triangulation<dim,spacedim> &tria,
+ const unsigned int tree_index,
+ const typename Triangulation<dim,spacedim>::cell_iterator &dealii_cell,
+ const typename dealii::internal::p4est::types<dim>::quadrant &p4est_cell,
+ const std::map<unsigned int, std::set<dealii::types::subdomain_id> > &vertices_with_ghost_neighbors,
+ const std::vector<bool> &vertex_locally_moved,
+ std::map<dealii::types::subdomain_id, CellInfo<dim, spacedim> > &needs_to_get_cell)
+ {
+ // see if we have to
+ // recurse...
+ if (dealii_cell->has_children())
+ {
+ typename dealii::internal::p4est::types<dim>::quadrant
+ p4est_child[GeometryInfo<dim>::max_children_per_cell];
+ dealii::internal::p4est::init_quadrant_children<dim>(p4est_cell, p4est_child);
+
+
+ for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
+ fill_vertices_recursively<dim,spacedim>(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<dealii::types::subdomain_id> send_to;
+ for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+ {
+ const std::map<unsigned int, std::set<dealii::types::subdomain_id> >::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<unsigned int> vertex_indices;
+ std::vector<dealii::Point<spacedim> > local_vertices;
+ for (unsigned int v=0; v<GeometryInfo<dim>::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<dealii::types::subdomain_id>::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<dealii::types::subdomain_id, CellInfo<dim, spacedim> >::iterator
+ p
+ = needs_to_get_cell.insert (std::make_pair(subdomain,
+ CellInfo<dim,spacedim>()))
+ .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 <int dim, int spacedim>
+ void
+ set_vertices_recursively (
+ const parallel::distributed::Triangulation<dim,spacedim> &tria,
+ const typename dealii::internal::p4est::types<dim>::quadrant &p4est_cell,
+ const typename Triangulation<dim,spacedim>::cell_iterator &dealii_cell,
+ const typename dealii::internal::p4est::types<dim>::quadrant &quadrant,
+ const dealii::Point<spacedim> *const vertices,
+ const unsigned int *const vertex_indices)
+ {
+ if (dealii::internal::p4est::quadrant_is_equal<dim>(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; i<n_vertices; ++i)
+ dealii_cell->vertex(vertex_indices[i+1]) = vertices[i];
+
+ return;
+ }
+
+ if (! dealii_cell->has_children())
+ return;
+
+ if (! dealii::internal::p4est::quadrant_is_ancestor<dim> (p4est_cell, quadrant))
+ return;
+
+ typename dealii::internal::p4est::types<dim>::quadrant
+ p4est_child[GeometryInfo<dim>::max_children_per_cell];
+ dealii::internal::p4est::init_quadrant_children<dim>(p4est_cell, p4est_child);
+
+ for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
+ set_vertices_recursively<dim,spacedim> (tria, p4est_child[c],
+ dealii_cell->child(c),
+ quadrant, vertices,
+ vertex_indices);
+ }
+ }
+ }
+
+
template <int dim, int spacedim>
void
+
+
+
+ template <int dim, int spacedim>
+ void
+ Triangulation<dim,spacedim>::
+ communicate_locally_moved_vertices (const std::vector<bool> &vertex_locally_moved)
+ {
+ Assert (vertex_locally_moved.size() == this->n_vertices(),
+ ExcDimensionMismatch(vertex_locally_moved.size(),
+ this->n_vertices()));
+
+ std::map<unsigned int, std::set<dealii::types::subdomain_id> >
+ vertices_with_ghost_neighbors;
+
+ // First find out which process should receive which vertices...
+ for (typename Triangulation<dim,spacedim>::active_cell_iterator
+ cell=this->begin_active(); cell!=this->end();
+ ++cell)
+ if (cell->is_ghost())
+ for (unsigned int vertex_no=0;
+ vertex_no<GeometryInfo<dim>::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<dealii::types::subdomain_id, CommunicateLocallyMovedVertices::CellInfo<dim,spacedim> > cellmap_t;
+ cellmap_t needs_to_get_cells;
+
+ for (typename Triangulation<dim,spacedim>::cell_iterator
+ cell = this->begin(0);
+ cell != this->end(0);
+ ++cell)
+ {
+ typename dealii::internal::p4est::types<dim>::quadrant p4est_coarse_cell;
+ dealii::internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);
+
+ CommunicateLocallyMovedVertices::fill_vertices_recursively<dim,spacedim>
+ (*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<std::vector<char> > sendbuffers (needs_to_get_cells.size());
+ std::vector<std::vector<char> >::iterator buffer = sendbuffers.begin();
+ std::vector<MPI_Request> requests (needs_to_get_cells.size());
+ std::vector<unsigned int> 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<unsigned int> senders
+ = Utilities::MPI::compute_point_to_point_communication_pattern
+ (this->get_communicator(), destinations);
+
+ // receive ghostcelldata
+ std::vector<char> receive;
+ CommunicateLocallyMovedVertices::CellInfo<dim,spacedim> cellinfo;
+ for (unsigned int i=0; i<senders.size(); ++i)
+ {
+ MPI_Status status;
+ int len;
+ MPI_Probe(MPI_ANY_SOURCE, 123, this->get_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<cells; ++c)
+ {
+ typename dealii::parallel::distributed::Triangulation<dim,spacedim>::cell_iterator
+ cell (this,
+ 0,
+ this->get_p4est_tree_to_coarse_cell_permutation()[cellinfo.tree_index[c]]);
+
+ typename dealii::internal::p4est::types<dim>::quadrant p4est_coarse_cell;
+ dealii::internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);
+
+ CommunicateLocallyMovedVertices::set_vertices_recursively<dim,spacedim> (*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 <int dim, int spacedim>
types::subdomain_id
Triangulation<dim,spacedim>::locally_owned_subdomain () const
+ template <int spacedim>
+ void
+ Triangulation<1,spacedim>::communicate_locally_moved_vertices
+ (const std::vector<bool> &vertex_locally_moved)
+ {
+ Assert (false, ExcNotImplemented());
+ }
+
+
template <int spacedim>
types::subdomain_id
Triangulation<1,spacedim>::locally_owned_subdomain () const
}
-#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 <int dim, int spacedim>
- struct CellInfo
- {
- // store all the tree_indices we send/receive consecutively (n_cells entries)
- std::vector<unsigned int> tree_index;
- // store all the quadrants we send/receive consecutively (n_cells entries)
- std::vector<typename dealii::internal::p4est::types<dim>::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<unsigned int> vertex_indices;
- // store for each cell the vertices we send/receive
- // (for each cell n_vertices entries)
- std::vector<dealii::Point<spacedim> > 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<unsigned int * > first_vertex_indices;
- std::vector<dealii::Point<spacedim>* > 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<dim>::quadrant) +
- vertices.size() * sizeof(dealii::Point<spacedim>)) +
- vertex_indices.size() * sizeof(unsigned int);
- }
-
- void pack_data (std::vector<char> &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<dim>::quadrant));
- ptr += num_cells*sizeof(typename dealii::internal::p4est::types<dim>::
- 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<spacedim>));
- ptr += vertices.size() * sizeof(dealii::Point<spacedim>);
-
- Assert (ptr == &buffer[0]+buffer.size(),
- ExcInternalError());
-
- }
-
- void unpack_data (const std::vector<char> &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<dim>::quadrant)*cells);
- ptr += sizeof(typename dealii::internal::p4est::types<dim>::quadrant)*cells;
-
- vertex_indices.clear();
- first_vertex_indices.resize(cells);
- std::vector<unsigned int> n_vertices_on_cell(cells);
- std::vector<unsigned int> first_indices (cells);
- for (unsigned int c=0; c<cells; ++c)
- {
- // The first 'vertex index' is the number of vertices.
- // Additionally, we need to store the pointer to this
- // vertex index with respect to the std::vector
- const unsigned int *const vertex_index
- = reinterpret_cast<const unsigned int *const>(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<cells; ++c)
- first_vertex_indices[c] = &vertex_indices[first_indices[c]];
-
- vertices.clear();
- first_vertices.resize(cells);
- for (unsigned int c=0; c<cells; ++c)
- {
- // We need to store a pointer to the first vertex.
- const dealii::Point<spacedim> *const vertex
- = reinterpret_cast<const dealii::Point<spacedim> * const>(ptr);
- first_indices[c] = vertices.size();
- vertices.push_back(*vertex);
- ptr += sizeof(dealii::Point<spacedim>);
- 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<spacedim>));
- ptr += (n_vertices_on_cell[c]-1)*sizeof(dealii::Point<spacedim>);
- }
- for (unsigned int c=0; c<cells; ++c)
- first_vertices[c] = &vertices[first_indices[c]];
-
- Assert (ptr == &buffer[0]+buffer.size(),
- ExcInternalError());
- }
- };
-
-
-
- // This function is responsible for gathering the information
- // we want to send to each process.
- // For each dealii cell on the coarsest level the corresponding
- // p4est_cell has to be provided when calling this function.
- // By recursing through all children we consider each active cell.
- // vertices_with_ghost_neighbors tells us which vertices
- // are in the ghost layer and for which processes they might
- // be interesting.
- // Whether a vertex has actually been updated locally is
- // stored in vertex_locally_moved. Only those are considered
- // for sending.
- // The gathered information is saved into needs_to_get_cell.
- template <int dim, int spacedim>
- void
- fill_vertices_recursively (const typename parallel::distributed::Triangulation<dim,spacedim> &tria,
- const unsigned int tree_index,
- const typename Triangulation<dim,spacedim>::cell_iterator &dealii_cell,
- const typename dealii::internal::p4est::types<dim>::quadrant &p4est_cell,
- const std::map<unsigned int, std::set<dealii::types::subdomain_id> > &vertices_with_ghost_neighbors,
- const std::vector<bool> &vertex_locally_moved,
- std::map<dealii::types::subdomain_id, CellInfo<dim, spacedim> > &needs_to_get_cell)
- {
- // see if we have to
- // recurse...
- if (dealii_cell->has_children())
- {
- typename dealii::internal::p4est::types<dim>::quadrant
- p4est_child[GeometryInfo<dim>::max_children_per_cell];
- internal::p4est::init_quadrant_children<dim>(p4est_cell, p4est_child);
-
-
- for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
- fill_vertices_recursively<dim,spacedim>(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<dealii::types::subdomain_id> send_to;
- for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
- {
- const std::map<unsigned int, std::set<dealii::types::subdomain_id> >::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<unsigned int> vertex_indices;
- std::vector<dealii::Point<spacedim> > local_vertices;
- for (unsigned int v=0; v<GeometryInfo<dim>::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<dealii::types::subdomain_id>::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<dealii::types::subdomain_id, CellInfo<dim, spacedim> >::iterator
- p
- = needs_to_get_cell.insert (std::make_pair(subdomain,
- CellInfo<dim,spacedim>()))
- .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 <int dim, int spacedim>
- void
- set_vertices_recursively (
- const parallel::distributed::Triangulation<dim,spacedim> &tria,
- const typename dealii::internal::p4est::types<dim>::quadrant &p4est_cell,
- const typename Triangulation<dim,spacedim>::cell_iterator &dealii_cell,
- const typename dealii::internal::p4est::types<dim>::quadrant &quadrant,
- const dealii::Point<spacedim> *const vertices,
- const unsigned int *const vertex_indices)
- {
- if (internal::p4est::quadrant_is_equal<dim>(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; i<n_vertices; ++i)
- dealii_cell->vertex(vertex_indices[i+1]) = vertices[i];
-
- return;
- }
-
- if (! dealii_cell->has_children())
- return;
-
- if (! internal::p4est::quadrant_is_ancestor<dim> (p4est_cell, quadrant))
- return;
-
- typename dealii::internal::p4est::types<dim>::quadrant
- p4est_child[GeometryInfo<dim>::max_children_per_cell];
- internal::p4est::init_quadrant_children<dim>(p4est_cell, p4est_child);
-
- for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
- set_vertices_recursively<dim,spacedim> (tria, p4est_child[c],
- dealii_cell->child(c),
- quadrant, vertices,
- vertex_indices);
- }
- }
-#endif //DEAL_II_WITH_P4EST
-
-
template <int dim, int spacedim>
double
- template <int spacedim>
- void
- communicate_locally_moved_vertices
- (parallel::distributed::Triangulation<1,spacedim> &triangulation,
- const std::vector<bool> &vertex_locally_moved)
- {
- Assert (false, ExcNotImplemented());
- }
-
-
-
- template <int dim, int spacedim>
- void
- communicate_locally_moved_vertices
- (parallel::distributed::Triangulation<dim,spacedim> &triangulation,
- const std::vector<bool> &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<unsigned int, std::set<dealii::types::subdomain_id> >
- vertices_with_ghost_neighbors;
-
- // First find out which process should receive which vertices...
- for (typename Triangulation<dim,spacedim>::active_cell_iterator
- cell=triangulation.begin_active(); cell!=triangulation.end();
- ++cell)
- if (cell->is_ghost())
- for (unsigned int vertex_no=0;
- vertex_no<GeometryInfo<dim>::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<dealii::types::subdomain_id, CellInfo<dim,spacedim> > cellmap_t;
- cellmap_t needs_to_get_cells;
-
- for (typename Triangulation<dim,spacedim>::cell_iterator
- cell = triangulation.begin(0);
- cell != triangulation.end(0);
- ++cell)
- {
- typename dealii::internal::p4est::types<dim>::quadrant p4est_coarse_cell;
- internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);
-
- fill_vertices_recursively<dim,spacedim>
- (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<std::vector<char> > sendbuffers (needs_to_get_cells.size());
- std::vector<std::vector<char> >::iterator buffer = sendbuffers.begin();
- std::vector<MPI_Request> requests (needs_to_get_cells.size());
- std::vector<unsigned int> 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<unsigned int> senders
- = Utilities::MPI::compute_point_to_point_communication_pattern
- (triangulation.get_communicator(), destinations);
-
- // receive ghostcelldata
- std::vector<char> receive;
- CellInfo<dim,spacedim> cellinfo;
- for (unsigned int i=0; i<senders.size(); ++i)
- {
- MPI_Status status;
- int len;
- MPI_Probe(MPI_ANY_SOURCE, 123, triangulation.get_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,
- triangulation.get_communicator(), &status);
-
- cellinfo.unpack_data(receive);
- const unsigned int cells = cellinfo.tree_index.size();
- for (unsigned int c=0; c<cells; ++c)
- {
- typename parallel::distributed::Triangulation<dim,spacedim>::cell_iterator
- cell (&triangulation,
- 0,
- triangulation.get_p4est_tree_to_coarse_cell_permutation()[cellinfo.tree_index[c]]);
-
- typename dealii::internal::p4est::types<dim>::quadrant p4est_coarse_cell;
- internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);
-
- set_vertices_recursively<dim,spacedim> (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.
}
if (distributed)
- communicate_locally_moved_vertices(*tr, vertex_moved);
+ tr->communicate_locally_moved_vertices(vertex_moved);
// Correct hanging nodes if necessary
Triangulation<deal_II_dimension, deal_II_space_dimension> &,
const bool);
- template
- void communicate_locally_moved_vertices (parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension> &,
- const std::vector<bool> &);
-
template
void get_face_connectivity_of_cells
(const Triangulation<deal_II_dimension, deal_II_space_dimension> &triangulation,
projects / dealii.git / commitdiff