From b9f78d92708d992bd1b225c19bb2853fa5a38d76 Mon Sep 17 00:00:00 2001
From: Daniel Arndt <d.arndt@math.uni-goettingen.de>
Date: Tue, 5 May 2015 10:30:07 +0200
Subject: [PATCH] Fixed colorization for GridGenerator::half_hyper_shell

---
 doc/news/changes.h                    |  5 ++++
 include/deal.II/grid/grid_generator.h |  6 ++---
 source/grid/grid_generator.cc         | 33 +++++++++++++++++++++++++--
 3 files changed, 39 insertions(+), 5 deletions(-)

diff --git a/doc/news/changes.h b/doc/news/changes.h
index 822adabc74..e0a089473e 100644
--- a/doc/news/changes.h
+++ b/doc/news/changes.h
@@ -472,6 +472,11 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li> Fixed: GridGenerator::half_hyper_shell can now be colorized.
+  <br>
+  (Daniel Arndt, 2015/05/05)
+  </li>
+
   <li> New: The function VectorTools::point_gradient has been added to compute
   the gradient of a given FE function.
   <br>
diff --git a/include/deal.II/grid/grid_generator.h b/include/deal.II/grid/grid_generator.h
index 3f51f26a31..0fbc8f959f 100644
--- a/include/deal.II/grid/grid_generator.h
+++ b/include/deal.II/grid/grid_generator.h
@@ -506,9 +506,9 @@ namespace GridGenerator
    * computed adaptively such that the resulting elements have the least
    * aspect ratio.
    *
-   * If colorize is set to true, the inner, outer, left, and right boundary
-   * get indicator 0, 1, 2, and 3, respectively. Otherwise all indicators are
-   * set to 0.
+   * If colorize is set to true, the inner, outer, and the part of the boundary
+   * where $x=0$, get indicator 0, 1, and 2, respectively. Otherwise all
+   * indicators are set to 0.
    *
    * @note The triangulation needs to be void upon calling this function.
    */
diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc
index 4c2439ecde..7e7fa1641c 100644
--- a/source/grid/grid_generator.cc
+++ b/source/grid/grid_generator.cc
@@ -3026,10 +3026,8 @@ namespace GridGenerator
                     const unsigned int n,
                     const bool colorize)
   {
-    (void)colorize;
     Assert ((inner_radius > 0) && (inner_radius < outer_radius),
             ExcInvalidRadii ());
-    Assert(colorize == false, ExcNotImplemented());
 
     if (n <= 5)
       {
@@ -3091,7 +3089,38 @@ namespace GridGenerator
       {
         Assert(false, ExcIndexRange(n, 0, 5));
       }
+    if (colorize)
+      {
+        // We want to use a standard boundary description where
+        // the boundary is not curved. Hence set boundary id 2 to
+        // to all faces in a first step.
+        typename Triangulation<3>::cell_iterator cell = tria.begin();
+        for (; cell!=tria.end(); ++cell)
+          for (unsigned int i=0; i<GeometryInfo<3>::faces_per_cell; ++i)
+            if (cell->at_boundary(i))
+              cell->face(i)->set_all_boundary_ids(2);
+
+        // Next look for the curved boundaries. If the first component is
+        // not near to zero, the boundary is curved. Then decide whether the center
+        // is nearer to the inner or outer boundary to set the correct boundary id.
+        for (cell=tria.begin(); cell!=tria.end(); ++cell)
+          for (unsigned int i=0; i<GeometryInfo<3>::faces_per_cell; ++i)
+            if (cell->at_boundary(i))
+              {
+                const typename Triangulation<3>::face_iterator face
+                  = cell->face(i);
 
+                const Point<3> face_center (face->center());
+                if (std::abs((face_center(0)-center(0)) > 1.e-3))
+                  {
+                    if (std::abs((face_center-center).norm()-inner_radius) <
+                        std::abs((face_center-center).norm()-outer_radius))
+                      face->set_all_boundary_ids(0);
+                    else
+                      face->set_all_boundary_ids(1);
+                  }
+              }
+      }
   }
 
 
-- 
2.39.5