]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Move communicate_locally_moved_vertices to parallel::distributed::Triangulation.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Sat, 8 Nov 2014 21:09:48 +0000 (15:09 -0600)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 12 Nov 2014 13:03:24 +0000 (07:03 -0600)
include/deal.II/distributed/tria.h
include/deal.II/grid/grid_tools.h
source/distributed/tria.cc
source/grid/grid_tools.cc
source/grid/grid_tools.inst.in

index f0b515b2ef0d4e18cc7450d6034bc33c95dab66c..4fb06a477e0dc1f75cbdd8d8236f8b1a62078ee7 100644 (file)
@@ -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<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
@@ -958,6 +1007,42 @@ namespace parallel
       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
index 15a06de672a909dd6c1a5d99d0d1255574ca8532..eb00e598fcacda7358e9cd3b5975076471fd5b6b 100644 (file)
@@ -373,26 +373,6 @@ namespace GridTools
                        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
index 07f4a6986b9e4f1aeefbcbe2da14b6f6a1abbd93..643e3e3a64a6aa0753b1d1d5e1788b95b4e61c72 100644 (file)
@@ -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 <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
@@ -3182,6 +3485,146 @@ namespace parallel
 
 
 
+
+
+
+    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
@@ -3911,6 +4354,15 @@ namespace parallel
 
 
 
+    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
index d66d12b7f0b75539ecfccb38cd486ac66048225b..b736aadf55cae0edaf487ef7900f934dc718fe6f 100644 (file)
@@ -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 <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
@@ -937,162 +635,6 @@ namespace GridTools
 
 
 
-  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.
@@ -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
index 91e6ef835c5410b160443f48ddf338233ddef3ef..a11c8cbf729e2b89d29e762068d215fb06cad897 100644 (file)
@@ -185,10 +185,6 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
                                     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,

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.