From d4cb9852e2adaad7435072a3f7f0e58e48a75415 Mon Sep 17 00:00:00 2001
From: heister <heister@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Tue, 4 Feb 2014 18:43:20 +0000
Subject: [PATCH] make GridTools::transform work with hanging nodes

git-svn-id: https://svn.dealii.org/trunk@32402 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/doc/news/changes.h                |  5 ++
 deal.II/include/deal.II/grid/grid_tools.h | 67 +++++++++++++++++------
 2 files changed, 56 insertions(+), 16 deletions(-)

diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index be608b1720..6aba67d558 100644
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -124,6 +124,11 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li> Changed: GridTools::transform() can now deal with meshes with hanging nodes.
+  <br>
+  (Timo Heister, 2014/02/04)
+  </li>
+
   <li>Fixed: Calling FEValuesViews::Vector::get_function_curls() computed
   wrong results in some cases (see https://code.google.com/p/dealii/issues/detail?id=182).
   This is now fixed.
diff --git a/deal.II/include/deal.II/grid/grid_tools.h b/deal.II/include/deal.II/grid/grid_tools.h
index 5f43b3af9e..e8de5d95ab 100644
--- a/deal.II/include/deal.II/grid/grid_tools.h
+++ b/deal.II/include/deal.II/grid/grid_tools.h
@@ -171,15 +171,7 @@ namespace GridTools
   /**
    * Transform the vertices of the given
    * triangulation by applying the
-   * function object provided as first argument to all its vertices. Since
-   * the internal consistency of a
-   * triangulation can only be guaranteed
-   * if the transformation is applied to
-   * the vertices of only one level of
-   * hierarchically refined cells, this
-   * function may only be used if all cells
-   * of the triangulation are on the same
-   * refinement level.
+   * function object provided as first argument to all its vertices.
    *
    * The transformation given as
    * argument is used to transform
@@ -194,6 +186,13 @@ namespace GridTools
    * return value have to be of
    * type <tt>Point<spacedim></tt>.
    *
+   * Note: if you are using a parallel::distributed::Triangulation you will have
+   * hanging nodes in your local Triangulation even if your "global" mesh has
+   * no hanging nodes. This will cause issues with wrong positioning of hanging
+   * nodes in ghost cells. The active cells will be correct, but keep in mind that
+   * computations like KellyErrorEstimator will give wrong answers. A safe bet is
+   * to use this function prior to any refinement in parallel.
+   *
    * This function is used in the
    * "Possibilities for extensions" section
    * of step-38. It is also used in step-49.
@@ -1204,13 +1203,6 @@ namespace GridTools
   void transform (const Predicate    &predicate,
                   Triangulation<dim, spacedim> &triangulation)
   {
-    // ensure that all the cells of the
-    // triangulation are on the same level
-    Assert (triangulation.n_levels() ==
-            static_cast<unsigned int>(triangulation.begin_active()->level()+1),
-            ExcMessage ("Not all cells of this triangulation are at the same "
-                        "refinement level, as is required for this function."));
-
     std::vector<bool> treated_vertices (triangulation.n_vertices(),
                                         false);
 
@@ -1232,6 +1224,49 @@ namespace GridTools
             // and mark it as treated
             treated_vertices[cell->vertex_index(v)] = true;
           };
+
+
+    // now fix any vertices on hanging nodes so that we don't create any holes
+    if (dim==2)
+      {
+        typename Triangulation<dim,spacedim>::active_cell_iterator
+        cell = triangulation.begin_active(),
+        endc = triangulation.end();
+        for (; cell!=endc; ++cell)
+          for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+            if (cell->face(face)->has_children() &&
+                !cell->face(face)->at_boundary())
+              { // this line has children
+                cell->face(face)->child(0)->vertex(1)
+                                       = (cell->face(face)->vertex(0) +
+                                           cell->face(face)->vertex(1)) / 2;
+              }
+      }
+    else if (dim==3)
+      {
+        typename Triangulation<dim,spacedim>::active_cell_iterator
+        cell = triangulation.begin_active(),
+        endc = triangulation.end();
+        for (; cell!=endc; ++cell)
+          for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+            if (cell->face(face)->has_children() &&
+                !cell->face(face)->at_boundary())
+              { // this face has hanging nodes
+                cell->face(face)->child(0)->vertex(1)
+                    = (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) / 2.0;
+                cell->face(face)->child(0)->vertex(2)
+                    = (cell->face(face)->vertex(0) + cell->face(face)->vertex(2)) / 2.0;
+                cell->face(face)->child(1)->vertex(3)
+                    = (cell->face(face)->vertex(1) + cell->face(face)->vertex(3)) / 2.0;
+                cell->face(face)->child(2)->vertex(3)
+                    = (cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) / 2.0;
+
+                // center of the face
+                cell->face(face)->child(0)->vertex(3)
+                    = (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)
+                        + cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) / 4.0;
+              }
+      }
   }
 
 
-- 
2.39.5