From: heister <heister@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Wed, 4 Sep 2013 15:28:48 +0000 (+0000)
Subject: fix bug in parallel MG: need to import additional level ghost neighbors which are... 
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=037b4fafbe726d66ab7c7a62538d686c2d956a53;p=dealii-svn.git

fix bug in parallel MG: need to import additional level ghost neighbors which are coarser

git-svn-id: https://svn.dealii.org/trunk@30595 0785d39b-7218-0410-832d-ea1e28bc413d
---

diff --git a/deal.II/source/distributed/tria.cc b/deal.II/source/distributed/tria.cc
index d4b8767c96..4269403891 100644
--- a/deal.II/source/distributed/tria.cc
+++ b/deal.II/source/distributed/tria.cc
@@ -2795,7 +2795,32 @@ namespace parallel
                 }
             }
 
+          //step 4: Special case: on each level we need all the face neighbors of our own level cells
+          // these are normally on the same level, unless the neighbor is active and coarser. It
+          // can end up on a different processor. Luckily, the level_subdomain_id can be figured out
+          // without communication, because the cell is active (and so level_subdomain_id=subdomain_id):
+          for (typename Triangulation<dim,spacedim>::cell_iterator cell = this->begin(); cell!=this->end(); ++cell)
+            {
+              if (cell->level_subdomain_id()!=this->locally_owned_subdomain())
+                continue;
 
+              for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+                {
+                  if (cell->face(f)->at_boundary())
+                    continue;
+                  if (cell->neighbor(f)->level() < cell->level()
+                      &&
+                      cell->neighbor(f)->level_subdomain_id()!=this->locally_owned_subdomain())
+                    {
+                      Assert(cell->neighbor(f)->active(), ExcInternalError());
+                      Assert(cell->neighbor(f)->subdomain_id() != numbers::artificial_subdomain_id, ExcInternalError());
+                      Assert(cell->neighbor(f)->level_subdomain_id() == numbers::artificial_subdomain_id
+                           || cell->neighbor(f)->level_subdomain_id() == cell->neighbor(f)->subdomain_id(), ExcInternalError());
+                      cell->neighbor(f)->set_level_subdomain_id(cell->neighbor(f)->subdomain_id());
+                    }
+                }
+
+            }
 
         }
 
diff --git a/deal.II/source/dofs/dof_handler_policy.cc b/deal.II/source/dofs/dof_handler_policy.cc
index 30a9c96ecb..496a9e16f4 100644
--- a/deal.II/source/dofs/dof_handler_policy.cc
+++ b/deal.II/source/dofs/dof_handler_policy.cc
@@ -1212,6 +1212,40 @@ namespace internal
                     }
                 }
 
+              //additionally (multigrid only), we can have the case that children of our neighbor
+              //have us as a neighbor. In this case we and the children are active.
+              if (dealii_cell->active())
+                {
+                  for (unsigned int f=0;f<GeometryInfo<dim>::faces_per_cell;++f)
+                    {
+                      if (dealii_cell->at_boundary(f))
+                        continue;
+                      typename DoFHandler<dim,spacedim>::level_cell_iterator neighbor = dealii_cell->neighbor(f);
+                      if (!neighbor->has_children())
+                        continue;
+
+                      for (unsigned int subface=0; subface<GeometryInfo<dim>::max_children_per_face; ++subface)
+                        {
+                          typename DoFHandler<dim,spacedim>::level_cell_iterator child = dealii_cell->neighbor_child_on_subface(f,subface);
+                          dealii::types::subdomain_id dest = child->subdomain_id();
+                          Assert(dest != dealii::numbers::artificial_subdomain_id, ExcInternalError());
+                          if (dest != tria.locally_owned_subdomain())
+                            {
+                              if (send_to.find(dest) == send_to.end())
+                                {
+                                  std::cout << "********" << std::endl;
+                                }
+                              send_to.insert(dest);
+
+                            }
+                        }
+
+                    }
+
+                }
+
+
+              // send if we have something to send
               if (send_to.size() > 0)
                 {
                   // this cell's dof_indices
@@ -2320,7 +2354,8 @@ namespace internal
               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())
+                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)]
@@ -2364,7 +2399,7 @@ namespace internal
                                      local_dof_indices.end(),
                                      DoFHandler<dim,spacedim>::invalid_dof_index))
                         {
-                          Assert(false, ExcMessage ("not all dofs got distributed!"));
+                          Assert(false, ExcMessage ("not all DoFs got distributed!"));
                         }
                     }
               }