From: Daniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Date: Tue, 22 Mar 2016 15:35:35 +0000 (+0100)
Subject: Use only periodioc face_map for the multigrid structures in distributed/tria.h
X-Git-Tag: v8.5.0-rc1~1151^2~3
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b4ba1602dfdfe72579542651b86703c00b755fd1;p=dealii.git

Use only periodioc face_map for the multigrid structures in distributed/tria.h
---

diff --git a/source/distributed/tria.cc b/source/distributed/tria.cc
index 2e052d49e2..605000a0bc 100644
--- a/source/distributed/tria.cc
+++ b/source/distributed/tria.cc
@@ -3124,91 +3124,6 @@ namespace parallel
 
     namespace
     {
-      // this function combines vertices that have different locations (and
-      // thus, different vertex_index) but represent the same topological
-      // entity over periodic boundaries. The vector
-      // topological_vertex_numbering contains a linear map from 0 to
-      // n_vertices at input and at output relates periodic vertices with only
-      // one vertex index. The output is used to always identify the same
-      // vertex according to the periodicity, e.g. when finding the maximum
-      // cell level around a vertex.
-      //
-      // Example: On a 3D cell with vertices numbered from 0 to 7 and periodic
-      // boundary conditions in x direction, the vector
-      // topological_vertex_numbering will contain the numbers
-      // {0,0,2,2,4,4,6,6} (because the vertex pairs {0,1}, {2,3}, {4,5},
-      // {6,7} belong together, respectively). If periodicity is set in x and
-      // z direction, the output is {0,0,2,2,0,0,2,2}, and if periodicity is
-      // in all directions, the output is simply {0,0,0,0,0,0,0,0}.
-      template <typename ITERATOR>
-      void
-      identify_periodic_vertices_recursively(const GridTools::PeriodicFacePair<ITERATOR> &periodic,
-                                             std::vector<unsigned int> &topological_vertex_numbering)
-      {
-        const unsigned int dim = ITERATOR::AccessorType::dimension;
-
-        // for hanging nodes we will consider all necessary coupling already
-        // on the parent level, so we just need to consider neighbors of the
-        // same level
-        if (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)
-              {
-                const unsigned int child_index_0 =
-                  GeometryInfo<dim>::child_cell_on_face(periodic.cell[0]->refinement_case(),
-                                                        periodic.face_idx[0], cf,
-                                                        periodic.orientation[0],
-                                                        periodic.orientation[1],
-                                                        periodic.orientation[2]);
-                periodic_child.cell[0] = periodic.cell[0]->child(child_index_0);
-                periodic_child.face_idx[0] = periodic.face_idx[0];
-
-                // the second face is in standard orientation in terms of the
-                // periodic face pair
-                const unsigned int child_index_1 =
-                  GeometryInfo<dim>::child_cell_on_face(periodic.cell[1]->refinement_case(),
-                                                        periodic.face_idx[1], cf);
-
-                periodic_child.cell[1] = periodic.cell[1]->child(child_index_1);
-                periodic_child.face_idx[1] = periodic.face_idx[1];
-
-                // recursive call into children
-                identify_periodic_vertices_recursively (periodic_child,
-                                                        topological_vertex_numbering);
-              }
-          }
-
-        for (unsigned int v=0; v<GeometryInfo<dim-1>::vertices_per_cell; ++v)
-          {
-            // take possible non-standard orientation of face on cell[0] into
-            // account
-            const unsigned int vface0 =
-              GeometryInfo<dim>::standard_to_real_face_vertex(v,periodic.orientation[0],
-                                                              periodic.orientation[1],
-                                                              periodic.orientation[2]);
-            const unsigned int vi0 = topological_vertex_numbering[periodic.cell[0]->face(periodic.face_idx[0])->vertex_index(vface0)];
-            const unsigned int vi1 = topological_vertex_numbering[periodic.cell[1]->face(periodic.face_idx[1])->vertex_index(v)];
-            const unsigned int min_index = std::min(vi0, vi1);
-            topological_vertex_numbering[periodic.cell[0]->face(periodic.face_idx[0])->vertex_index(vface0)]
-              = topological_vertex_numbering[periodic.cell[1]->face(periodic.face_idx[1])->vertex_index(v)]
-                = min_index;
-
-            /*  std::cout << "Old identified vertices " << periodic.cell[0]->face(periodic.face_idx[0])->vertex_index(vface0)
-                        << " at "                 << periodic.cell[0]->face(periodic.face_idx[0])->vertex(vface0)
-                        << " and "                << periodic.cell[1]->face(periodic.face_idx[1])->vertex_index(v)
-                        << " at "                 << periodic.cell[1]->face(periodic.face_idx[1])->vertex(v)
-                        << std::endl;*/
-
-          }
-      }
-
-
-
       // ensures the 2:1 mesh balance for periodic boundary conditions in the
       // artificial cell layer (the active cells are taken care of by p4est)
       template <int dim, int spacedim>
@@ -3223,18 +3138,8 @@ namespace parallel
         tria.save_refine_flags (flags_before[1]);
 
         std::vector<unsigned int> topological_vertex_numbering(tria.n_vertices());
-        std::vector<unsigned int> topological_vertex_numbering_cmp(tria.n_vertices());
         for (unsigned int i=0; i<topological_vertex_numbering.size(); ++i)
-          {
-            topological_vertex_numbering[i] = i;
-            topological_vertex_numbering_cmp[i] = i;
-          }
-        //std::cout << std::endl;
-        for (unsigned int i=0; i<periodic_face_pairs_level_0.size(); ++i)
-          {
-            identify_periodic_vertices_recursively(periodic_face_pairs_level_0[i],
-                                                   topological_vertex_numbering_cmp);
-          }
+          topological_vertex_numbering[i] = i;
         // combine vertices that have different locations (and thus, different
         // vertex_index) but represent the same topological entity over periodic
         // boundaries. The vector topological_vertex_numbering contains a linear
@@ -3277,11 +3182,6 @@ namespace parallel
                     topological_vertex_numbering[cell_1->face(face_no_1)->vertex_index(vface0)]
                       = topological_vertex_numbering[cell_2->face(face_no_2)->vertex_index(v)]
                         = min_index;
-                    /*  std::cout << "Newly identified vertices " << cell_1->face(face_no_1)->vertex_index(vface0)
-                              << " at "                 << cell_1->face(face_no_1)->vertex(vface0)
-                              << " and "                << cell_2->face(face_no_2)->vertex_index(v)
-                              << " at "                 << cell_2->face(face_no_2)->vertex(v)
-                              << std::endl;*/
                   }
               }
           }
@@ -3295,27 +3195,6 @@ namespace parallel
                      ExcInternalError());
           }
 
-        for (unsigned int i=0; i<topological_vertex_numbering_cmp.size(); ++i)
-          {
-            const unsigned int j = topological_vertex_numbering_cmp[i];
-            if (j != i)
-              Assert(topological_vertex_numbering_cmp[j] == j,
-                     ExcInternalError());
-          }
-
-        if (topological_vertex_numbering != topological_vertex_numbering_cmp)
-          {
-            for (unsigned int i=0; i<topological_vertex_numbering.size(); ++i)
-              {
-                if (topological_vertex_numbering[i] != topological_vertex_numbering_cmp[i])
-                  std::cout << i
-                            << " old: " << topological_vertex_numbering_cmp[i]
-                            << " new: " << topological_vertex_numbering[i]
-                            << std::endl;
-              }
-            Assert( false, ExcInternalError());
-          }
-
 
         // this code is replicated from grid/tria.cc but using an indirection
         // for periodic boundary conditions
@@ -4554,167 +4433,6 @@ 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;
-
-        const unsigned int dim = ITERATOR::AccessorType::dimension;
-        // 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)
-              {
-                const unsigned int child_index_0 =
-                  GeometryInfo<dim>::child_cell_on_face(periodic.cell[0]->refinement_case(),
-                                                        periodic.face_idx[0], cf,
-                                                        periodic.orientation[0],
-                                                        periodic.orientation[1],
-                                                        periodic.orientation[2]);
-                periodic_child.cell[0] = periodic.cell[0]->child(child_index_0);
-                periodic_child.face_idx[0] = periodic.face_idx[0];
-
-                // the second face is in standard orientation in terms of the
-                // periodic face pair
-                const unsigned int child_index_1 =
-                  GeometryInfo<dim>::child_cell_on_face(periodic.cell[1]->refinement_case(),
-                                                        periodic.face_idx[1], cf);
-
-                periodic_child.cell[1] = periodic.cell[1]->child(child_index_1);
-                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;
-
-        for (unsigned int v=0; v<GeometryInfo<dim-1>::vertices_per_cell; ++v)
-          {
-            // take possible non-standard orientation of face on cell[0] into
-            // account
-            const unsigned int vface0 =
-              GeometryInfo<dim>::standard_to_real_face_vertex(v,periodic.orientation[0],
-                                                              periodic.orientation[1],
-                                                              periodic.orientation[2]);
-            if (active_vertices_on_level[periodic.cell[0]->face(periodic.face_idx[0])->vertex_index(vface0)] ||
-                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(vface0)]
-                = 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
@@ -4785,14 +4503,6 @@ namespace parallel
               vertices_with_ghost_neighbors[cell->vertex_index(v)]
               .insert (cell->level_subdomain_id());
 
-      std::map<unsigned int, std::set<dealii::types::subdomain_id> >
-      vertices_with_ghost_neighbors_cmp = vertices_with_ghost_neighbors;
-
-
-      for (unsigned int i=0; i<this->periodic_face_pairs_level_0.size(); ++i)
-        set_periodic_ghost_neighbors_recursively(this->periodic_face_pairs_level_0[i],
-                                                 level, vertices_with_ghost_neighbors_cmp);
-
       //now for the vertices on periodic faces
       typename std::map<std::pair<cell_iterator, unsigned int>,
                std::pair<std::pair<cell_iterator,unsigned int>, std::bitset<3> > >::const_iterator it;
@@ -4826,31 +4536,6 @@ namespace parallel
                 }
             }
         }
-
-      if (vertices_with_ghost_neighbors != vertices_with_ghost_neighbors_cmp)
-        {
-          std::cout << "Print vertices_with_ghost_neighbors: " << std::endl;
-          std::map<unsigned int, std::set<dealii::types::subdomain_id> >::const_iterator it;
-          for (it = vertices_with_ghost_neighbors.begin(); it!=vertices_with_ghost_neighbors.end(); ++it)
-            {
-              std::cout << "Vertex " << it->first << " is on subdomains ";
-              std::set<dealii::types::subdomain_id>::const_iterator it2;
-              for (it2 = it->second.begin(); it2!=it->second.end(); ++it2)
-                std::cout << *it2 << " ";
-              std::cout << std::endl;
-            }
-
-          std::cout << "Print vertices_with_ghost_neighbors_cmp: " << std::endl;
-          for (it = vertices_with_ghost_neighbors_cmp.begin(); it!=vertices_with_ghost_neighbors_cmp.end(); ++it)
-            {
-              std::cout << "Vertex " << it->first << " is on subdomains ";
-              std::set<dealii::types::subdomain_id>::const_iterator it2;
-              for (it2 = it->second.begin(); it2!=it->second.end(); ++it2)
-                std::cout << *it2 << " ";
-              std::cout << std::endl;
-            }
-          Assert(false, ExcInternalError());
-        }
     }
 
 
@@ -4870,12 +4555,6 @@ namespace parallel
           for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
             marked_vertices[cell->vertex_index(v)] = true;
 
-      std::vector<bool> marked_vertices_cmp = marked_vertices;
-
-      for (unsigned int i=0; i<this->periodic_face_pairs_level_0.size(); ++i)
-        mark_periodic_vertices_recursively(this->periodic_face_pairs_level_0[i],
-                                           level, marked_vertices_cmp);
-
       /**
        * ensure 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
@@ -4911,7 +4590,6 @@ namespace parallel
             }
         }
 
-      Assert(marked_vertices == marked_vertices_cmp, ExcInternalError());
       return marked_vertices;
     }