]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fix ghost layer on MG levels for periodic boundary conditions
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Sun, 18 Oct 2015 07:43:47 +0000 (09:43 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 19 Oct 2015 13:41:33 +0000 (15:41 +0200)
include/deal.II/distributed/tria.h
source/distributed/tria.cc
source/dofs/dof_handler_policy.cc

index 291101d33349e8a9cc3a03402abfdc60ae3e343f..229c5e864d3e5f8c308dd62ad6e0388f868cf2fb 100644 (file)
@@ -915,12 +915,21 @@ namespace parallel
        * grid cells are badly ordered this may mean that individual parts of
        * the forest stored on a local machine may be split across coarse grid
        * cells that are not geometrically close. Consequently, we apply a
-       * Cuthill-McKee preordering to ensure that the part of the forest
-       * stored by p4est is located on geometrically close coarse grid cells.
+       * hierarchical preordering according to
+       * SparsityTools::reorder_hierarchical() to ensure that the part of the
+       * forest stored by p4est is located on geometrically close coarse grid
+       * cells.
        */
       std::vector<types::global_dof_index> coarse_cell_to_p4est_tree_permutation;
       std::vector<types::global_dof_index> p4est_tree_to_coarse_cell_permutation;
 
+      /**
+       * If add_periodicity() is called, this variable stores the given
+       * periodic face pairs on level 0 for later access during the
+       * identification of ghost cells for the multigrid hierarchy.
+       */
+      std::vector<GridTools::PeriodicFacePair<cell_iterator> > periodic_face_pairs_level_0;
+
       /**
        * Return a pointer to the p4est tree that belongs to the given
        * dealii_coarse_cell_index()
@@ -957,8 +966,8 @@ namespace parallel
       void attach_mesh_data();
 
       /**
-       * fills a map that, for each vertex, lists all the processors whose
-       * subdomains are adjacent to that vertex.  Used by
+       * Fills a map that, for each vertex, lists all the processors whose
+       * subdomains are adjacent to that vertex. Used by
        * DoFHandler::Policy::ParallelDistributed.
        */
       void
@@ -966,6 +975,29 @@ namespace parallel
       (std::map<unsigned int, std::set<dealii::types::subdomain_id> >
        &vertices_with_ghost_neighbors);
 
+      /**
+       * Fills a map that, for each vertex, lists all the processors whose
+       * subdomains are adjacent to that vertex on the given level for the
+       * multigrid hierarchy. Used by DoFHandler::Policy::ParallelDistributed.
+       */
+      void
+      fill_level_vertices_with_ghost_neighbors
+      (const unsigned int level,
+       std::map<unsigned int, std::set<dealii::types::subdomain_id> >
+       &vertices_with_ghost_neighbors);
+
+      /**
+       * This method returns a bit vector of length tria.n_vertices()
+       * indicating the locally active vertices on a level, i.e., the vertices
+       * touched by the locally owned level cells for use in geometric
+       * multigrid (possibly including the vertices due to periodic boundary
+       * conditions) are marked by true.
+       *
+       * Used by DoFHandler::Policy::ParallelDistributed.
+       */
+      std::vector<bool>
+      mark_locally_active_vertices_on_level(const unsigned int level) const;
+
       template <int, int> friend class dealii::internal::DoFHandler::Policy::ParallelDistributed;
     };
 
@@ -1068,6 +1100,23 @@ namespace parallel
       (std::map<unsigned int, std::set<dealii::types::subdomain_id> >
        &vertices_with_ghost_neighbors);
 
+      /**
+       * Like above, this method, which is only implemented for dim = 2 or 3,
+       * needs a stub because it is used in dof_handler_policy.cc
+       */
+      void
+      fill_level_vertices_with_ghost_neighbors
+      (const unsigned int level,
+       std::map<unsigned int, std::set<dealii::types::subdomain_id> >
+       &vertices_with_ghost_neighbors);
+
+      /**
+       * Like above, this method, which is only implemented for dim = 2 or 3,
+       * needs a stub because it is used in dof_handler_policy.cc
+       */
+      std::vector<bool>
+      mark_locally_active_vertices_on_level(const unsigned int level) const;
+
     };
   }
 }
index 09a7961a93d4a5d8eed92eba90f2710207e2649d..5c92934683e83681f5cd732abfa7e80490a4f91c 100644 (file)
@@ -2585,6 +2585,8 @@ namespace parallel
       coarse_cell_to_p4est_tree_permutation.resize (0);
       p4est_tree_to_coarse_cell_permutation.resize (0);
 
+      periodic_face_pairs_level_0.clear();
+
       dealii::Triangulation<dim,spacedim>::clear ();
 
       this->update_number_cache ();
@@ -3337,7 +3339,7 @@ namespace parallel
           // -1. Note that we do not fill other cells we could figure out the
           // same way, because we might accidentally set an id for a cell that
           // is not a ghost cell.
-          for (unsigned int lvl=this->n_levels(); lvl>0;)
+          for (unsigned int lvl=this->n_levels(); lvl>0; )
             {
               --lvl;
               for (typename Triangulation<dim,spacedim>::cell_iterator cell = this->begin(lvl); cell!=this->end(lvl); ++cell)
@@ -3356,20 +3358,8 @@ namespace parallel
           //step 2: make sure all the neighbors to our level_cells exist. Need
           //to look up in p4est...
           std::vector<std::vector<bool> > marked_vertices(this->n_levels());
-
-          for (unsigned int lvl=0; lvl<this->n_levels(); ++lvl)
-            {
-              marked_vertices[lvl].resize(this->n_vertices(), false);
-
-              for (typename dealii::Triangulation<dim,spacedim>::cell_iterator
-                   cell = this->begin(lvl);
-                   cell != this->end(lvl); ++cell)
-                if (cell->level_subdomain_id() == this->locally_owned_subdomain())
-                  for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
-                    marked_vertices[lvl][cell->vertex_index(v)] = true;
-            }
-
-
+          for (unsigned int lvl=0; lvl < this->n_levels(); ++lvl)
+            marked_vertices[lvl] = mark_locally_active_vertices_on_level(lvl);
 
           for (typename Triangulation<dim,spacedim>::cell_iterator cell = this->begin(0); cell!=this->end(0); ++cell)
             {
@@ -4195,6 +4185,171 @@ namespace parallel
     }
 
 
+
+    namespace
+    {
+      /**
+       * ensures that if one of the two vertices on a periodic face is marked
+       * as active (i.e., belonging to an owned level cell), also the other
+       * one is active
+       */
+      template <typename ITERATOR>
+      void
+      mark_periodic_vertices_recursively(const GridTools::PeriodicFacePair<ITERATOR> &periodic,
+                                         const int target_level,
+                                         std::vector<bool> &active_vertices_on_level)
+      {
+        if (periodic.cell[0]->level() > target_level)
+          return;
+
+        // for hanging nodes there is nothing to do since we are interested in
+        // the connections on the same level...
+        if (periodic.cell[0]->level() < target_level &&
+            periodic.cell[0]->has_children() &&
+            periodic.cell[1]->has_children())
+          {
+            // copy orientations etc. from parent to child
+            GridTools::PeriodicFacePair<ITERATOR> periodic_child = periodic;
+
+            // find appropriate pairs of child elements
+            for (unsigned int cf=0; cf<periodic.cell[0]->face(periodic.face_idx[0])->n_children(); ++cf)
+              {
+                unsigned int c=0;
+                for (; c<periodic.cell[0]->n_children(); ++c)
+                  {
+                    if (periodic.cell[0]->child(c)->face(periodic.face_idx[0]) ==
+                        periodic.cell[0]->face(periodic.face_idx[0])->child(cf))
+                      break;
+                  }
+                Assert(c < periodic.cell[0]->n_children(),
+                       ExcMessage("Face child not found"));
+                periodic_child.cell[0] = periodic.cell[0]->child(c);
+                periodic_child.face_idx[0] = periodic.face_idx[0];
+
+                c=0;
+                for (; c<periodic.cell[1]->n_children(); ++c)
+                  {
+                    if (periodic.cell[1]->child(c)->face(periodic.face_idx[1]) ==
+                        periodic.cell[1]->face(periodic.face_idx[1])->child(cf))
+                      break;
+                  }
+                Assert(c < periodic.cell[1]->n_children(),
+                       ExcMessage("Face child not found"));
+                periodic_child.cell[1] = periodic.cell[1]->child(c);
+                periodic_child.face_idx[1] = periodic.face_idx[1];
+
+                // recursive call into children
+                mark_periodic_vertices_recursively (periodic_child, target_level,
+                                                    active_vertices_on_level);
+              }
+            return;
+          }
+
+        if (periodic.cell[0]->level() != target_level)
+          return;
+
+        // TODO: fix non-standard orientation
+        Assert(periodic.orientation[0] == true &&
+               periodic.orientation[1] == false &&
+               periodic.orientation[2] == false,
+               ExcNotImplemented());
+
+        for (unsigned int v=0; v<GeometryInfo<ITERATOR::AccessorType::dimension-1>::vertices_per_cell; ++v)
+          if (active_vertices_on_level[periodic.cell[0]->face(periodic.face_idx[0])->vertex_index(v)] ||
+              active_vertices_on_level[periodic.cell[1]->face(periodic.face_idx[1])->vertex_index(v)])
+            active_vertices_on_level[periodic.cell[0]->face(periodic.face_idx[0])->vertex_index(v)]
+              = active_vertices_on_level[periodic.cell[1]->face(periodic.face_idx[1])->vertex_index(v)]
+                = true;
+      }
+
+
+
+      /**
+       * ensures that always both vertices over a periodic face are identified
+       * together
+       */
+      template <typename ITERATOR>
+      void
+      set_periodic_ghost_neighbors_recursively(const GridTools::PeriodicFacePair<ITERATOR> &periodic,
+                                               const int target_level,
+                                               std::map<unsigned int, std::set<dealii::types::subdomain_id> > &vertices_with_ghost_neighbors)
+      {
+        if (periodic.cell[0]->level() > target_level)
+          return;
+
+        // for hanging nodes there is nothing to do since we are interested in
+        // the connections on the same level...
+        if (periodic.cell[0]->level() < target_level &&
+            periodic.cell[0]->has_children() &&
+            periodic.cell[1]->has_children())
+          {
+            // copy orientations etc. from parent to child
+            GridTools::PeriodicFacePair<ITERATOR> periodic_child = periodic;
+
+            // find appropriate pairs of child elements
+            for (unsigned int cf=0; cf<periodic.cell[0]->face(periodic.face_idx[0])->n_children(); ++cf)
+              {
+                unsigned int c=0;
+                for (; c<periodic.cell[0]->n_children(); ++c)
+                  {
+                    if (periodic.cell[0]->child(c)->face(periodic.face_idx[0]) ==
+                        periodic.cell[0]->face(periodic.face_idx[0])->child(cf))
+                      break;
+                  }
+                Assert(c < periodic.cell[0]->n_children(),
+                       ExcMessage("Face child not found"));
+                periodic_child.cell[0] = periodic.cell[0]->child(c);
+                periodic_child.face_idx[0] = periodic.face_idx[0];
+
+                c=0;
+                for (; c<periodic.cell[1]->n_children(); ++c)
+                  {
+                    if (periodic.cell[1]->child(c)->face(periodic.face_idx[1]) ==
+                        periodic.cell[1]->face(periodic.face_idx[1])->child(cf))
+                      break;
+                  }
+                Assert(c < periodic.cell[1]->n_children(),
+                       ExcMessage("Face child not found"));
+                periodic_child.cell[1] = periodic.cell[1]->child(c);
+                periodic_child.face_idx[1] = periodic.face_idx[1];
+
+                // recursive call into children
+                set_periodic_ghost_neighbors_recursively (periodic_child, target_level,
+                                                          vertices_with_ghost_neighbors);
+              }
+            return;
+          }
+
+        if (periodic.cell[0]->level() != target_level)
+          return;
+
+        // TODO: fix non-standard orientation
+        Assert(periodic.orientation[0] == true &&
+               periodic.orientation[1] == false &&
+               periodic.orientation[2] == false,
+               ExcNotImplemented());
+
+        for (unsigned int v=0; v<GeometryInfo<ITERATOR::AccessorType::dimension-1>::vertices_per_cell; ++v)
+          {
+            const unsigned int
+            idx0 = periodic.cell[0]->face(periodic.face_idx[0])->vertex_index(v),
+            idx1 = periodic.cell[1]->face(periodic.face_idx[1])->vertex_index(v);
+            if (vertices_with_ghost_neighbors.find(idx0) !=
+                vertices_with_ghost_neighbors.end())
+              vertices_with_ghost_neighbors[idx1].
+              insert(vertices_with_ghost_neighbors[idx0].begin(),
+                     vertices_with_ghost_neighbors[idx0].end());
+            if (vertices_with_ghost_neighbors.find(idx1) !=
+                vertices_with_ghost_neighbors.end())
+              vertices_with_ghost_neighbors[idx0].
+              insert(vertices_with_ghost_neighbors[idx1].begin(),
+                     vertices_with_ghost_neighbors[idx1].end());
+          }
+      }
+    }
+
+
+
     /**
      * Determine the neighboring subdomains that are adjacent to each vertex.
      * This is achieved via the p4est_iterate/p8est_iterate tool
@@ -4240,6 +4395,58 @@ namespace parallel
     }
 
 
+
+    /**
+     * Determine the neighboring subdomains that are adjacent to each vertex
+     * on the given multigrid level
+     */
+    template <int dim, int spacedim>
+    void
+    Triangulation<dim,spacedim>::
+    fill_level_vertices_with_ghost_neighbors
+    (const unsigned int level,
+     std::map<unsigned int, std::set<dealii::types::subdomain_id> >
+     &vertices_with_ghost_neighbors)
+    {
+      const std::vector<bool> locally_active_vertices =
+        mark_locally_active_vertices_on_level(level);
+      for (cell_iterator cell = this->begin(level); cell != this->end(level); ++cell)
+        if (cell->level_subdomain_id() != dealii::numbers::artificial_subdomain_id
+            && cell->level_subdomain_id() != this->locally_owned_subdomain())
+          for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+            if (locally_active_vertices[cell->vertex_index(v)])
+              vertices_with_ghost_neighbors[cell->vertex_index(v)]
+              .insert (cell->level_subdomain_id());
+
+      for (unsigned int i=0; i<periodic_face_pairs_level_0.size(); ++i)
+        set_periodic_ghost_neighbors_recursively(periodic_face_pairs_level_0[i],
+                                                 level, vertices_with_ghost_neighbors);
+    }
+
+
+
+    template<int dim, int spacedim>
+    std::vector<bool>
+    Triangulation<dim,spacedim>
+    ::mark_locally_active_vertices_on_level (const unsigned int level) const
+    {
+      Assert (dim>1, ExcNotImplemented());
+
+      std::vector<bool> marked_vertices(this->n_vertices(), false);
+      for (cell_iterator cell = this->begin(level); cell != this->end(level); ++cell)
+        if (cell->level_subdomain_id() == this->locally_owned_subdomain())
+          for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+            marked_vertices[cell->vertex_index(v)] = true;
+
+      for (unsigned int i=0; i<periodic_face_pairs_level_0.size(); ++i)
+        mark_periodic_vertices_recursively(periodic_face_pairs_level_0[i],
+                                           level, marked_vertices);
+
+      return marked_vertices;
+    }
+
+
+
     template<int dim, int spacedim>
     void
     Triangulation<dim,spacedim>::add_periodicity
@@ -4377,6 +4584,10 @@ namespace parallel
           AssertThrow (false, ExcInternalError());
         }
 
+      periodic_face_pairs_level_0.insert(periodic_face_pairs_level_0.end(),
+                                         periodicity_vector.begin(),
+                                         periodicity_vector.end());
+
 #else
       Assert(false, ExcMessage ("Need p4est version >= 0.3.4.1!"));
 #endif
@@ -4603,6 +4814,26 @@ namespace parallel
     {
       Assert (false, ExcNotImplemented());
     }
+
+    template <int spacedim>
+    void
+    Triangulation<1,spacedim>::
+    fill_level_vertices_with_ghost_neighbors
+    (const unsigned int /*level*/,
+     std::map<unsigned int, std::set<dealii::types::subdomain_id> >
+     &/*vertices_with_ghost_neighbors*/)
+    {
+      Assert (false, ExcNotImplemented());
+    }
+
+    template <int spacedim>
+    std::vector<bool>
+    Triangulation<1,spacedim>::
+    mark_locally_active_vertices_on_level (const unsigned int) const
+    {
+      Assert (false, ExcNotImplemented());
+      return std::vector<bool>();
+    }
   }
 }
 
index 7dd0650180ac0cdc10a0cdceb847e0fce33fcda3..dda7e89b4fa62c2ccef4299ffd10f7898d31e612 100644 (file)
@@ -1658,52 +1658,6 @@ namespace internal
 
         }
 
-        /**
-         * Return a vector in which,
-         * for every vertex index, we
-         * mark whether a locally owned
-         * cell is adjacent.
-         */
-        template <int dim, int spacedim>
-        std::vector<bool>
-        mark_locally_active_vertices (const parallel::distributed::Triangulation<dim,spacedim> &triangulation)
-        {
-          std::vector<bool> locally_active_vertices (triangulation.n_vertices(),
-                                                     false);
-
-          for (typename dealii::Triangulation<dim,spacedim>::active_cell_iterator
-               cell = triangulation.begin_active();
-               cell != triangulation.end(); ++cell)
-            if (cell->subdomain_id() == triangulation.locally_owned_subdomain())
-              for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
-                locally_active_vertices[cell->vertex_index(v)] = true;
-
-          return locally_active_vertices;
-        }
-
-        /**
-         * Return a vector in which,
-         * for every vertex index, we
-         * mark whether a locally owned
-         * cell is adjacent.
-         */
-        template <int dim, int spacedim>
-        std::vector<bool>
-        mark_locally_active_vertices_on_level (const parallel::distributed::Triangulation<dim,spacedim> &triangulation, unsigned int level)
-        {
-          std::vector<bool> locally_active_vertices (triangulation.n_vertices(),
-                                                     false);
-
-          if (level<triangulation.n_levels())
-            for (typename dealii::Triangulation<dim,spacedim>::cell_iterator
-                 cell = triangulation.begin(level);
-                 cell != triangulation.end(level); ++cell)
-              if (cell->level_subdomain_id() == triangulation.locally_owned_subdomain())
-                for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
-                  locally_active_vertices[cell->vertex_index(v)] = true;
-
-          return locally_active_vertices;
-        }
 
 
         template <int spacedim>
@@ -1763,7 +1717,6 @@ namespace internal
             }
 
 
-
           //sending
           std::vector<std::vector<char> > sendbuffers (needs_to_get_cells.size());
           std::vector<std::vector<char> >::iterator buffer = sendbuffers.begin();
@@ -1997,7 +1950,6 @@ namespace internal
               }
 
 
-
           //sending
           std::vector<std::vector<char> > sendbuffers (needs_to_get_cells.size());
           std::vector<std::vector<char> >::iterator buffer = sendbuffers.begin();
@@ -2397,6 +2349,7 @@ namespace internal
       }
 
 
+
       template <int dim, int spacedim>
       void
       ParallelDistributed<dim, spacedim>::
@@ -2571,34 +2524,21 @@ namespace internal
                     cell->set_user_flag();
               }
 
-            //mark the vertices we are interested
-            //in, i.e. belonging to own and marked cells
+            //mark the vertices we are interested in, i.e. belonging to own
+            //and marked cells
             const std::vector<bool> locally_active_vertices
-              = mark_locally_active_vertices_on_level (*tr, level);
+              = tr->mark_locally_active_vertices_on_level (level);
 
-            // add each ghostcells'
-            // subdomain to the vertex and
-            // keep track of interesting
-            // neighbors
+            // add each ghostcells' subdomain to the vertex and keep track of
+            // interesting neighbors
             std::map<unsigned int, std::set<dealii::types::subdomain_id> >
             vertices_with_ghost_neighbors;
-            if (level < tr->n_levels())
-              for (typename DoFHandler<dim,spacedim>::level_cell_iterator
-                   cell = dof_handler.begin(level);
-                   cell != dof_handler.end(level); ++cell)
-                if (cell->level_subdomain_id() != dealii::numbers::artificial_subdomain_id
-                    && cell->level_subdomain_id() != tr->locally_owned_subdomain())
-                  for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
-                    if (locally_active_vertices[cell->vertex_index(v)])
-                      vertices_with_ghost_neighbors[cell->vertex_index(v)]
-                      .insert (cell->level_subdomain_id());
-
-
-            /* Send and receive cells. After this,
-            only the local cells are marked,
-            that received new data. This has to
-            be communicated in a second
-            communication step. */
+            tr->fill_level_vertices_with_ghost_neighbors(level,
+                                                         vertices_with_ghost_neighbors);
+
+            // Send and receive cells. After this, only the local cells are
+            // marked, that received new data. This has to be communicated in
+            // a second communication step.
 
             communicate_mg_dof_indices_on_marked_cells( dof_handler,
                                                         vertices_with_ghost_neighbors,

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.