]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Change GridGenerator from a class to a namespace. Undo introduction of local typedef...
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 24 Jul 2013 12:58:44 +0000 (12:58 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 24 Jul 2013 12:58:44 +0000 (12:58 +0000)
git-svn-id: https://svn.dealii.org/trunk@30149 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/grid/grid_generator.h
deal.II/source/grid/grid_generator.cc
deal.II/source/grid/grid_generator.inst.in

index 59758c96011b47449da81db96d77db49b10b166d..3af5c44c42f6f9ee48c9ea01f423b97c2c0f3fcb 100644 (file)
@@ -44,7 +44,22 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
-
+  <li>Changed: GridGenerator used to be a class with only static members
+  but is now a namespace, like all other similar constructs in deal.II.
+  <br>
+  (Wolfgang Bangerth, 2013/07/24)
+  </li>
+
+  <li>Changed: In GridGenerator, several functions had erroneously been changed
+  to take an argument of type <code>size_type</code> rather than <code>unsigned
+  int</code>. <code>GridGenerator::size_type</code> was a typedef to
+  types::global_dof_index, which for most users was <code>unsigned int</code>
+  anyway, but could also be set to be a 64-bit integer type. In any case, the
+  change has been reverted and these functions take just a regular
+  <code>unsigned int</code> again.
+  <br>
+  (Wolfgang Bangerth, 2013/07/24)
+  </li>
 </ol>
 
 
index 39676698608ff92918e5df574b2d09817f0cf79f..58d4d3e2cf1087f04c0477fb844f58a5978dd9b0 100644 (file)
@@ -33,34 +33,25 @@ template <typename number> class SparseMatrix;
 
 
 /**
- * This class provides a collection of functions for generating basic
- * triangulations. Below, we try to provide some pictures in order to
- * illustrate at least the more complex ones.
+ * This namespace provides a collection of functions for generating
+ * triangulations for some basic geometries.
  *
  * Some of these functions receive a flag @p colorize. If this is
- * set, parts of the boundary receive different boundary numbers,
- * allowing them to be distinguished by application programs. See the
- * documentation of the functions for details.
- *
- * Additionally this class provides a function
- * (@p laplace_transformation) that smoothly transforms a grid
- * according to given new boundary points. This can be used to
- * transform (simple-shaped) grids to a more complicated ones, like a
- * shell onto a grid of an airfoil, for example.
- *
- * No meshes for the codimension one case are provided at the moment.
+ * set, parts of the boundary receive different boundary indicators
+ * (@ref GlossBoundaryIndicator),
+ * allowing them to be distinguished for the purpose of attaching geometry
+ * objects and evaluating different boundary conditions.
  *
+ * This namespace also provides a function
+ * GridGenerator::laplace_transformation that smoothly transforms a domain
+ * into another one. This can be used to
+ * transform basic geometries to more complicated ones, like a
+ * shell to a grid of an airfoil, for example.
  *
  * @ingroup grid
  */
-class GridGenerator
+namespace GridGenerator
 {
-public:
-  /**
-   * Declare type for number of cell.
-   */
-  typedef types::global_dof_index size_type;
-
   /**
    * Initialize the given triangulation with a hypercube (line in 1D, square
    * in 2D, etc) consisting of exactly one cell. The hypercube volume is the
@@ -79,9 +70,9 @@ public:
    * @note The triangulation needs to be void upon calling this function.
    */
   template <int dim, int spacedim>
-  static void hyper_cube (Triangulation<dim,spacedim>  &tria,
-                          const double        left = 0.,
-                          const double        right= 1.);
+  void hyper_cube (Triangulation<dim,spacedim>  &tria,
+                   const double        left = 0.,
+                   const double        right= 1.);
 
   /**
    * Same as hyper_cube(), but with the difference that not only one cell is
@@ -97,10 +88,10 @@ public:
    * @note The triangulation needs to be void upon calling this function.
    */
   template <int dim>
-  static void subdivided_hyper_cube (Triangulation<dim>  &tria,
-                                     const unsigned int  repetitions,
-                                     const double        left = 0.,
-                                     const double        right= 1.);
+  void subdivided_hyper_cube (Triangulation<dim>  &tria,
+                              const unsigned int  repetitions,
+                              const double        left = 0.,
+                              const double        right= 1.);
 
   /**
    * Create a coordinate-parallel brick from the two diagonally opposite
@@ -118,10 +109,10 @@ public:
    * @note The triangulation needs to be void upon calling this function.
    */
   template <int dim, int spacedim>
-  static void hyper_rectangle (Triangulation<dim,spacedim> &tria,
-                               const Point<spacedim>       &p1,
-                               const Point<spacedim>       &p2,
-                               const bool                  colorize = false);
+  void hyper_rectangle (Triangulation<dim,spacedim> &tria,
+                        const Point<spacedim>       &p1,
+                        const Point<spacedim>       &p2,
+                        const bool                  colorize = false);
 
   /**
    * Create a coordinate-parallel parallelepiped from the two diagonally
@@ -155,7 +146,6 @@ public:
    * program.
    */
   template <int dim>
-  static
   void
   subdivided_hyper_rectangle (Triangulation<dim>              &tria,
                               const std::vector<unsigned int> &repetitions,
@@ -179,7 +169,6 @@ public:
    * specified by the points @p p1 and @p p2.
    */
   template <int dim>
-  static
   void
   subdivided_hyper_rectangle(Triangulation<dim>                      &tria,
                              const std::vector<std::vector<double> > &step_sizes,
@@ -196,7 +185,6 @@ public:
    * i.e. the domain will have a void there.
    */
   template <int dim>
-  static
   void
   subdivided_hyper_rectangle (Triangulation<dim>                       &tria,
                               const std::vector< std::vector<double> > &spacing,
@@ -215,7 +203,6 @@ public:
    * @note The triangulation needs to be void upon calling this function.
    */
   template <int dim>
-  static
   void
   parallelogram(Triangulation<dim> &tria,
                 const Point<dim> (&corners)[dim],
@@ -225,7 +212,6 @@ public:
    * @deprecated Use the other function of same name.
    */
   template <int dim>
-  static
   void
   parallelogram(Triangulation<dim> &tria,
                 const Tensor<2,dim> &corners,
@@ -249,7 +235,6 @@ public:
    * function.
    */
   template <int dim>
-  static
   void
   parallelepiped (Triangulation<dim>  &tria,
                   const Point<dim>   (&corners) [dim],
@@ -268,10 +253,9 @@ public:
    * function.
    */
   template <int dim>
-  static
   void
   subdivided_parallelepiped (Triangulation<dim>  &tria,
-                             const size_type     n_subdivisions,
+                             const unsigned int   n_subdivisions,
                              const Point<dim>   (&corners) [dim],
                              const bool           colorize = false);
 
@@ -284,10 +268,9 @@ public:
    * function.
    */
   template <int dim>
-  static
   void
   subdivided_parallelepiped (Triangulation<dim>  &tria,
-                             const size_type    ( n_subdivisions) [dim],
+                             const unsigned int    ( n_subdivisions) [dim],
                              const Point<dim>   (&corners) [dim],
                              const bool           colorize = false);
 
@@ -307,11 +290,11 @@ public:
    * @note The triangulation needs to be void upon calling this function.
    */
   template <int dim>
-  static void enclosed_hyper_cube (Triangulation<dim> &tria,
-                                   const double      left = 0.,
-                                   const double      right= 1.,
-                                   const double      thickness = 1.,
-                                   const bool        colorize = false);
+  void enclosed_hyper_cube (Triangulation<dim> &tria,
+                            const double      left = 0.,
+                            const double      right= 1.,
+                            const double      thickness = 1.,
+                            const bool        colorize = false);
 
   /**
    * Initialize the given triangulation with a hyperball, i.e. a circle or a
@@ -328,9 +311,9 @@ public:
    * @note The triangulation needs to be void upon calling this function.
    */
   template <int dim>
-  static void hyper_ball (Triangulation<dim> &tria,
-                          const Point<dim>   &center = Point<dim>(),
-                          const double      radius = 1.);
+  void hyper_ball (Triangulation<dim> &tria,
+                   const Point<dim>   &center = Point<dim>(),
+                   const double      radius = 1.);
 
   /**
    * This class produces a half hyper-ball around <tt>center</tt>, which
@@ -346,9 +329,9 @@ public:
    * @note The triangulation needs to be void upon calling this function.
    */
   template <int dim>
-  static void half_hyper_ball (Triangulation<dim> &tria,
-                               const Point<dim>   &center = Point<dim>(),
-                               const double      radius = 1.);
+  void half_hyper_ball (Triangulation<dim> &tria,
+                        const Point<dim>   &center = Point<dim>(),
+                        const double      radius = 1.);
 
   /**
    * Create a cylinder around the x-axis.  The cylinder extends from
@@ -366,9 +349,9 @@ public:
    * @note The triangulation needs to be void upon calling this function.
    */
   template <int dim>
-  static void cylinder (Triangulation<dim> &tria,
-                        const double      radius = 1.,
-                        const double      half_length = 1.);
+  void cylinder (Triangulation<dim> &tria,
+                 const double      radius = 1.,
+                 const double      half_length = 1.);
 
   /**
    * Create a cutted cone around the x-axis.  The cone extends from
@@ -395,7 +378,7 @@ public:
    * @author Markus B&uuml;rg, 2009
    */
   template <int dim>
-  static void
+  void
   truncated_cone (Triangulation<dim> &tria,
                   const double radius_0 = 1.0,
                   const double radius_1 = 0.5,
@@ -417,9 +400,9 @@ public:
    * @note The triangulation needs to be void upon calling this function.
    */
   template <int dim>
-  static void hyper_L (Triangulation<dim> &tria,
-                       const double      left = -1.,
-                       const double      right= 1.);
+  void hyper_L (Triangulation<dim> &tria,
+                const double      left = -1.,
+                const double      right= 1.);
 
   /**
    * Initialize the given Triangulation with a hypercube with a slit. In each
@@ -437,10 +420,10 @@ public:
    * @note The triangulation needs to be void upon calling this function.
    */
   template <int dim>
-  static void hyper_cube_slit (Triangulation<dim> &tria,
-                               const double      left = 0.,
-                               const double      right= 1.,
-                               const bool colorize = false);
+  void hyper_cube_slit (Triangulation<dim> &tria,
+                        const double      left = 0.,
+                        const double      right= 1.,
+                        const bool colorize = false);
 
   /**
    * Produce a hyper-shell, the region between two spheres around
@@ -491,12 +474,12 @@ public:
    * @note The triangulation needs to be void upon calling this function.
    */
   template <int dim>
-  static void hyper_shell (Triangulation<dim>   &tria,
-                           const Point<dim>     &center,
-                           const double        inner_radius,
-                           const double        outer_radius,
-                           const size_type     n_cells = 0,
-                           bool colorize = false);
+  void hyper_shell (Triangulation<dim>   &tria,
+                    const Point<dim>     &center,
+                    const double        inner_radius,
+                    const double        outer_radius,
+                    const unsigned int     n_cells = 0,
+                    bool colorize = false);
 
   /**
    * Produce a half hyper-shell, i.e. the space between two circles in two
@@ -520,12 +503,12 @@ public:
    * @note The triangulation needs to be void upon calling this function.
    */
   template <int dim>
-  static void half_hyper_shell (Triangulation<dim>   &tria,
-                                const Point<dim>     &center,
-                                const double        inner_radius,
-                                const double        outer_radius,
-                                const size_type     n_cells = 0,
-                                const bool colorize = false);
+  void half_hyper_shell (Triangulation<dim>   &tria,
+                         const Point<dim>     &center,
+                         const double        inner_radius,
+                         const double        outer_radius,
+                         const unsigned int     n_cells = 0,
+                         const bool colorize = false);
 
 
   /**
@@ -549,12 +532,12 @@ public:
    * @note The triangulation needs to be void upon calling this function.
    */
   template <int dim>
-  static void quarter_hyper_shell (Triangulation<dim>   &tria,
-                                   const Point<dim>     &center,
-                                   const double        inner_radius,
-                                   const double        outer_radius,
-                                   const size_type     n_cells = 0,
-                                   const bool colorize = false);
+  void quarter_hyper_shell (Triangulation<dim>   &tria,
+                            const Point<dim>     &center,
+                            const double        inner_radius,
+                            const double        outer_radius,
+                            const unsigned int     n_cells = 0,
+                            const bool colorize = false);
 
   /**
    * Produce a domain that is the space between two cylinders in 3d, with
@@ -569,12 +552,12 @@ public:
    * @note The triangulation needs to be void upon calling this function.
    */
   template <int dim>
-  static void cylinder_shell (Triangulation<dim>   &tria,
-                              const double        length,
-                              const double        inner_radius,
-                              const double        outer_radius,
-                              const size_type     n_radial_cells = 0,
-                              const size_type     n_axial_cells = 0);
+  void cylinder_shell (Triangulation<dim>   &tria,
+                       const double        length,
+                       const double        inner_radius,
+                       const double        outer_radius,
+                       const unsigned int     n_radial_cells = 0,
+                       const unsigned int     n_axial_cells = 0);
 
 
 
@@ -592,9 +575,9 @@ public:
    * torus.
    */
 
-  static void torus (Triangulation<2,3>  &tria,
-                     const double         R,
-                     const double         r);
+  void torus (Triangulation<2,3>  &tria,
+              const double         R,
+              const double         r);
 
 
   /**
@@ -622,12 +605,12 @@ public:
    *    get the number 0 and the hole gets number 1.
    */
   template<int dim>
-  static void hyper_cube_with_cylindrical_hole (Triangulation<dim> &triangulation,
-                                                const double inner_radius = .25,
-                                                const double outer_radius = .5,
-                                                const double L = .5,
-                                                const size_type repetition = 1,
-                                                const bool colorize = false);
+  void hyper_cube_with_cylindrical_hole (Triangulation<dim> &triangulation,
+                                         const double inner_radius = .25,
+                                         const double outer_radius = .5,
+                                         const double L = .5,
+                                         const unsigned int repetition = 1,
+                                         const bool colorize = false);
 
   /**
    * Produce a ring of cells in 3D that is cut open, twisted and glued
@@ -639,11 +622,11 @@ public:
    * @param R           The radius of the circle, which forms the middle line of the torus containing the loop of cells. Must be greater than @p r.
    * @param r           The radius of the cylinder bend together as loop.
    */
-  static void moebius (Triangulation<3,3>  &tria,
-                       const size_type      n_cells,
-                       const unsigned int   n_rotations,
-                       const double         R,
-                       const double         r);
+  void moebius (Triangulation<3,3>  &tria,
+                const unsigned int      n_cells,
+                const unsigned int   n_rotations,
+                const double         R,
+                const double         r);
 
   /**
    * Given the two triangulations specified as the first two arguments, create
@@ -676,7 +659,6 @@ public:
    * GridTools::create_union_triangulation .
    */
   template <int dim, int spacedim>
-  static
   void
   merge_triangulations (const Triangulation<dim, spacedim> &triangulation_1,
                         const Triangulation<dim, spacedim> &triangulation_2,
@@ -690,10 +672,9 @@ public:
    * to the corresponding side walls in z direction. The bottom and top
    * get the next two free boundary indicators.
    */
-  static
   void
   extrude_triangulation(const Triangulation<2, 2> &input,
-                        const size_type n_slices,
+                        const unsigned int n_slices,
                         const double height,
                         Triangulation<3,3> &result);
 
@@ -712,9 +693,8 @@ public:
    * @deprecated This function has been moved to GridTools::laplace_transform
    */
   template <int dim>
-  static
   void laplace_transformation (Triangulation<dim> &tria,
-                               const std::map<size_type,Point<dim> > &new_points) DEAL_II_DEPRECATED;
+                               const std::map<unsigned int,Point<dim> > &new_points) DEAL_II_DEPRECATED;
 
   /**
    * Exception
@@ -734,67 +714,7 @@ public:
                   int,
                   << "The vector of repetitions  must have "
                   << arg1 <<" elements.");
-
-private:
-  /**
-   * Perform the action specified by the @p colorize flag of the
-   * hyper_rectangle() function of this class.
-   */
-  template <int dim, int spacedim>
-  static
-  void
-  colorize_hyper_rectangle (Triangulation<dim,spacedim> &tria);
-
-  /**
-   * Perform the action specified by the @p colorize flag of the
-   * subdivided_hyper_rectangle() function of this class. This function is
-   * singled out because it is dimension specific.
-   */
-  template <int dim>
-  static
-  void
-  colorize_subdivided_hyper_rectangle (Triangulation<dim> &tria,
-                                       const Point<dim>   &p1,
-                                       const Point<dim>   &p2,
-                                       const double        epsilon);
-
-  /**
-   * Assign boundary number zero to the inner shell boundary and 1 to the
-   * outer.
-   */
-  template<int dim>
-  static
-  void
-  colorize_hyper_shell (Triangulation<dim> &tria,
-                        const Point<dim> &center,
-                        const double inner_radius,
-                        const double outer_radius);
-
-
-  /**
-   * Assign boundary number zero the inner shell boundary, one to the outer
-   * shell boundary, two to the face with x=0, three to the face with y=0,
-   * four to the face with z=0.
-   */
-  template<int dim>
-  static
-  void
-  colorize_quarter_hyper_shell(Triangulation<dim> &tria,
-                               const Point<dim> &center,
-                               const double inner_radius,
-                               const double outer_radius);
-
-  /**
-   * Solve the Laplace equation for @p laplace_transformation function for one
-   * of the @p dim space dimensions. Externalized into a function of its own
-   * in order to allow parallel execution.
-   */
-  static
-  void
-  laplace_solve (const SparseMatrix<double>          &S,
-                 const std::map<size_type,double> &m,
-                 Vector<double>                      &u);
-};
+}
 
 
 
index e834ca633467cd8a19b156e2a6a46535d03e7bf9..443cc9d4795a26d1f28629c524bb703a749a9290 100644 (file)
@@ -46,3784 +46,3793 @@ DEAL_II_NAMESPACE_OPEN
 
 
 
-namespace
+namespace GridGenerator
 {
-  // Corner points of the cube [-1,1]^3
-  const Point<3> hexahedron[8] =
+  namespace
   {
-    Point<3>(-1,-1,-1),
-    Point<3>(+1,-1,-1),
-    Point<3>(-1,+1,-1),
-    Point<3>(+1,+1,-1),
-    Point<3>(-1,-1,+1),
-    Point<3>(+1,-1,+1),
-    Point<3>(-1,+1,+1),
-    Point<3>(+1,+1,+1)
-  };
-
-  // Octahedron inscribed in the cube
-  // [-1,1]^3
-  const Point<3> octahedron[6] =
-  {
-    Point<3>(-1, 0, 0),
-    Point<3>( 1, 0, 0),
-    Point<3>( 0,-1, 0),
-    Point<3>( 0, 1, 0),
-    Point<3>( 0, 0,-1),
-    Point<3>( 0, 0, 1)
-  };
-}
-
+    // Corner points of the cube [-1,1]^3
+    const Point<3> hexahedron[8] =
+    {
+      Point<3>(-1,-1,-1),
+      Point<3>(+1,-1,-1),
+      Point<3>(-1,+1,-1),
+      Point<3>(+1,+1,-1),
+      Point<3>(-1,-1,+1),
+      Point<3>(+1,-1,+1),
+      Point<3>(-1,+1,+1),
+      Point<3>(+1,+1,+1)
+    };
 
-template <int dim, int spacedim>
-void
-GridGenerator::hyper_rectangle (Triangulation<dim,spacedim> &tria,
-                                const Point<spacedim>   &p_1,
-                                const Point<spacedim>   &p_2,
-                                const bool          colorize)
-{
-  // First, normalize input such that
-  // p1 is lower in all coordinate directions.
-  Point<spacedim> p1(p_1);
-  Point<spacedim> p2(p_2);
+    // Octahedron inscribed in the cube
+    // [-1,1]^3
+    const Point<3> octahedron[6] =
+    {
+      Point<3>(-1, 0, 0),
+      Point<3>( 1, 0, 0),
+      Point<3>( 0,-1, 0),
+      Point<3>( 0, 1, 0),
+      Point<3>( 0, 0,-1),
+      Point<3>( 0, 0, 1)
+    };
 
-  for (unsigned int i=0; i<spacedim; ++i)
-    if (p1(i) > p2(i))
-      std::swap (p1(i), p2(i));
 
-  std::vector<Point<spacedim> > vertices (GeometryInfo<dim>::vertices_per_cell);
-  switch (dim)
+    /**
+     * Perform the action specified by the @p colorize flag of the
+     * hyper_rectangle() function of this class.
+     */
+    template <int dim, int spacedim>
+    void
+    colorize_hyper_rectangle (Triangulation<dim,spacedim> &tria)
     {
-    case 1:
-      vertices[0] = p1;
-      vertices[1] = p2;
-      break;
-    case 2:
-      vertices[0] = vertices[1] = p1;
-      vertices[2] = vertices[3] = p2;
-
-      vertices[1](0) = p2(0);
-      vertices[2](0) = p1(0);
-      break;
-    case 3:
-      vertices[0] = vertices[1] = vertices[2] = vertices[3] = p1;
-      vertices[4] = vertices[5] = vertices[6] = vertices[7] = p2;
+      // there is nothing to do in 1d
+      if (dim > 1)
+        {
+          // there is only one cell, so
+          // simple task
+          const typename Triangulation<dim,spacedim>::cell_iterator
+          cell = tria.begin();
+          for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+            cell->face(f)->set_boundary_indicator (f);
+        }
+    }
 
-      vertices[1](0) = p2(0);
-      vertices[2](1) = p2(1);
-      vertices[3](0) = p2(0);
-      vertices[3](1) = p2(1);
 
-      vertices[4](0) = p1(0);
-      vertices[4](1) = p1(1);
-      vertices[5](1) = p1(1);
-      vertices[6](0) = p1(0);
 
-      break;
-    default:
-      Assert (false, ExcNotImplemented ());
+    void
+    colorize_subdivided_hyper_rectangle (Triangulation<1> &tria,
+                                         const Point<1> &,
+                                         const Point<1> &,
+                                         const double)
+    {
+      for (Triangulation<1>::cell_iterator cell = tria.begin();
+           cell != tria.end(); ++cell)
+        if (cell->center()(0) > 0)
+          cell->set_material_id(1);
+      // boundary indicators are set to
+      // 0 (left) and 1 (right) by default.
     }
 
-  // Prepare cell data
-  std::vector<CellData<dim> > cells (1);
-  for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
-    cells[0].vertices[i] = i;
-  cells[0].material_id = 0;
 
-  tria.create_triangulation (vertices, cells, SubCellData());
 
-  // Assign boundary indicators
-  if (colorize)
-    colorize_hyper_rectangle (tria);
-}
+    template <int dim>
+    void
+    colorize_subdivided_hyper_rectangle (Triangulation<dim> &tria,
+                                         const Point<dim>   &p1,
+                                         const Point<dim>   &p2,
+                                         const double        epsilon)
+    {
 
+      // run through all faces and check
+      // if one of their center coordinates matches
+      // one of the corner points. Comparisons
+      // are made using an epsilon which
+      // should be smaller than the smallest cell
+      // diameter.
 
+      typename Triangulation<dim>::face_iterator face = tria.begin_face(),
+                                                 endface = tria.end_face();
+      for (; face!=endface; ++face)
+        {
+          if (face->boundary_indicator() == 0)
+            {
+              const Point<dim> center (face->center());
+              if (std::abs(center(0)-p1[0]) < epsilon)
+                face->set_boundary_indicator(0);
+              else if (std::abs(center(0) - p2[0]) < epsilon)
+                face->set_boundary_indicator(1);
+              else if (dim > 1 && std::abs(center(1) - p1[1]) < epsilon)
+                face->set_boundary_indicator(2);
+              else if (dim > 1 && std::abs(center(1) - p2[1]) < epsilon)
+                face->set_boundary_indicator(3);
+              else if (dim > 2 && std::abs(center(2) - p1[2]) < epsilon)
+                face->set_boundary_indicator(4);
+              else if (dim > 2 && std::abs(center(2) - p2[2]) < epsilon)
+                face->set_boundary_indicator(5);
+              else
+                // triangulation says it
+                // is on the boundary,
+                // but we could not find
+                // on which boundary.
+                Assert (false, ExcInternalError());
 
-template <int dim, int spacedim>
-void
-GridGenerator::colorize_hyper_rectangle (Triangulation<dim,spacedim> &tria)
-{
-  // there is nothing to do in 1d
-  if (dim > 1)
-    {
-      // there is only one cell, so
-      // simple task
-      const typename Triangulation<dim,spacedim>::cell_iterator
-      cell = tria.begin();
-      for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
-        cell->face(f)->set_boundary_indicator (f);
+            }
+        }
+      for (typename Triangulation<dim>::cell_iterator cell = tria.begin();
+           cell != tria.end(); ++cell)
+        {
+          char id = 0;
+          for (unsigned int d=0; d<dim; ++d)
+            if (cell->center()(d) > 0) id += 1 << d;
+          cell->set_material_id(id);
+        }
     }
-}
-
 
 
-template <int dim, int spacedim>
-void GridGenerator::hyper_cube (Triangulation<dim,spacedim> &tria,
-                                const double                 left,
-                                const double                 right)
-{
-  Assert (left < right,
-          ExcMessage ("Invalid left-to-right bounds of hypercube"));
-
-  Point<spacedim> p1;
-  Point<spacedim> p2;
+    /**
+     * Assign boundary number zero to the inner shell boundary and 1 to the
+     * outer.
+     */
+    void colorize_hyper_shell (Triangulation<1> &,
+                               const Point<1> &,
+                               const double,
+                               const double)
+    {
+      Assert (false, ExcNotImplemented());
+    }
 
-  p1(spacedim-1) = 0;
-  p2(spacedim-1) = 0;
 
-  for (unsigned int i=0; i<dim; ++i)
+    /**
+     * Assign boundary number zero to the inner shell boundary and 1 to the
+     * outer.
+     */
+    void
+    colorize_hyper_shell (
+      Triangulation<2> &tria,
+      const Point<2> &, const double, const double)
     {
-      p1(i) = left;
-      p2(i) = right;
+      // In spite of receiving geometrical
+      // data, we do this only based on
+      // topology.
+
+      // For the mesh based on  cube,
+      // this is highly irregular
+      for (Triangulation<2>::cell_iterator cell = tria.begin();
+           cell != tria.end(); ++cell)
+        {
+          Assert (cell->face(2)->at_boundary(), ExcInternalError());
+          cell->face(2)->set_boundary_indicator(1);
+        }
     }
-  hyper_rectangle (tria, p1, p2);
-}
 
 
+    /**
+     * Assign boundary number zero to the inner shell boundary and 1 to the
+     * outer.
+     */
+    void
+    colorize_hyper_shell (Triangulation<3> &tria,
+                          const Point<3> &,
+                          const double,
+                          const double)
+    {
+      // the following uses a good amount
+      // of knowledge about the
+      // orientation of cells. this is
+      // probably not good style...
+      if (tria.n_cells() == 6)
+        {
+          Triangulation<3>::cell_iterator cell = tria.begin();
 
-void
-GridGenerator::moebius (
-  Triangulation<3>  &tria,
-  const size_type      n_cells,
-  const unsigned int   n_rotations,
-  const double         R,
-  const double         r)
-{
-  const unsigned int dim=3;
-  Assert (n_cells>4, ExcMessage("More than 4 cells are needed to create a moebius grid."));
-  Assert (r>0 && R>0, ExcMessage("Outer and inner radius must be positive."));
-  Assert (R>r, ExcMessage("Outer radius must be greater than inner radius."));
+          cell->face(4)->set_boundary_indicator(1);
+          Assert (cell->face(4)->at_boundary(), ExcInternalError());
 
+          (++cell)->face(2)->set_boundary_indicator(1);
+          Assert (cell->face(2)->at_boundary(), ExcInternalError());
 
-  std::vector<Point<dim> > vertices (4*n_cells);
-  double beta_step=n_rotations*numbers::PI/2.0/n_cells;
-  double alpha_step=2.0*numbers::PI/n_cells;
+          (++cell)->face(2)->set_boundary_indicator(1);
+          Assert (cell->face(2)->at_boundary(), ExcInternalError());
 
-  for (size_type i=0; i<n_cells; ++i)
-    for (unsigned int j=0; j<4; ++j)
-      {
-        vertices[4*i+j][0]=R*std::cos(i*alpha_step)+r*std::cos(i*beta_step+j*numbers::PI/2.0)*std::cos(i*alpha_step);
-        vertices[4*i+j][1]=R*std::sin(i*alpha_step)+r*std::cos(i*beta_step+j*numbers::PI/2.0)*std::sin(i*alpha_step);
-        vertices[4*i+j][2]=r*std::sin(i*beta_step+j*numbers::PI/2.0);
-      }
+          (++cell)->face(0)->set_boundary_indicator(1);
+          Assert (cell->face(0)->at_boundary(), ExcInternalError());
 
-  size_type offset=0;
+          (++cell)->face(2)->set_boundary_indicator(1);
+          Assert (cell->face(2)->at_boundary(), ExcInternalError());
 
-  std::vector<CellData<dim> > cells (n_cells);
-  for (size_type i=0; i<n_cells; ++i)
-    {
-      for (unsigned int j=0; j<2; ++j)
+          (++cell)->face(0)->set_boundary_indicator(1);
+          Assert (cell->face(0)->at_boundary(), ExcInternalError());
+        }
+      else if (tria.n_cells() == 12)
+        {
+          // again use some internal
+          // knowledge
+          for (Triangulation<3>::cell_iterator cell = tria.begin();
+               cell != tria.end(); ++cell)
+            {
+              Assert (cell->face(5)->at_boundary(), ExcInternalError());
+              cell->face(5)->set_boundary_indicator(1);
+            }
+        }
+      else if (tria.n_cells() == 96)
         {
-          cells[i].vertices[0+4*j]=offset+0+4*j;
-          cells[i].vertices[1+4*j]=offset+3+4*j;
-          cells[i].vertices[2+4*j]=offset+2+4*j;
-          cells[i].vertices[3+4*j]=offset+1+4*j;
+          // the 96-cell hypershell is
+          // based on a once refined
+          // 12-cell mesh. consequently,
+          // since the outer faces all
+          // are face_no==5 above, so
+          // they are here (unless they
+          // are in the interior). Use
+          // this to assign boundary
+          // indicators, but also make
+          // sure that we encounter
+          // exactly 48 such faces
+          unsigned int count = 0;
+          for (Triangulation<3>::cell_iterator cell = tria.begin();
+               cell != tria.end(); ++cell)
+            if (cell->face(5)->at_boundary())
+              {
+                cell->face(5)->set_boundary_indicator(1);
+                ++count;
+              }
+          Assert (count == 48, ExcInternalError());
         }
-      offset+=4;
-      cells[i].material_id=0;
+      else
+        Assert (false, ExcNotImplemented());
     }
 
-  // now correct the last four vertices
-  cells[n_cells-1].vertices[4]=(0+n_rotations)%4;
-  cells[n_cells-1].vertices[5]=(3+n_rotations)%4;
-  cells[n_cells-1].vertices[6]=(2+n_rotations)%4;
-  cells[n_cells-1].vertices[7]=(1+n_rotations)%4;
 
-  GridReordering<dim>::invert_all_cells_of_negative_grid(vertices,cells);
-  tria.create_triangulation_compatibility (vertices, cells, SubCellData());
-}
 
+    /**
+     * Assign boundary number zero the inner shell boundary, one to the outer
+     * shell boundary, two to the face with x=0, three to the face with y=0,
+     * four to the face with z=0.
+     */
+    void
+    colorize_quarter_hyper_shell(Triangulation<3> &tria,
+                                 const Point<3> &center,
+                                 const double inner_radius,
+                                 const double outer_radius)
+    {
+      if (tria.n_cells() != 3)
+        AssertThrow (false, ExcNotImplemented());
 
+      double middle = (outer_radius-inner_radius)/2e0 + inner_radius;
+      double eps = 1e-3*middle;
+      Triangulation<3>::cell_iterator cell = tria.begin();
 
-void
-GridGenerator::torus (Triangulation<2,3>  &tria,
-                      const double         R,
-                      const double         r)
-{
-  Assert (R>r, ExcMessage("Outer radius must be greater than inner radius."));
-
-  const unsigned int dim=2;
-  const unsigned int spacedim=3;
-  std::vector<Point<spacedim> > vertices (16);
-
-  vertices[0]=Point<spacedim>(R-r,0,0);
-  vertices[1]=Point<spacedim>(R,-r,0);
-  vertices[2]=Point<spacedim>(R+r,0,0);
-  vertices[3]=Point<spacedim>(R, r,0);
-  vertices[4]=Point<spacedim>(0,0,R-r);
-  vertices[5]=Point<spacedim>(0,-r,R);
-  vertices[6]=Point<spacedim>(0,0,R+r);
-  vertices[7]=Point<spacedim>(0,r,R);
-  vertices[8]=Point<spacedim>(-(R-r),0,0);
-  vertices[9]=Point<spacedim>(-R,-r,0);
-  vertices[10]=Point<spacedim>(-(R+r),0,0);
-  vertices[11]=Point<spacedim>(-R, r,0);
-  vertices[12]=Point<spacedim>(0,0,-(R-r));
-  vertices[13]=Point<spacedim>(0,-r,-R);
-  vertices[14]=Point<spacedim>(0,0,-(R+r));
-  vertices[15]=Point<spacedim>(0,r,-R);
-
-  std::vector<CellData<dim> > cells (16);
-  //Right Hand Orientation
-  cells[0].vertices[0] =  0;
-  cells[0].vertices[1] =  4;
-  cells[0].vertices[2] =  7;
-  cells[0].vertices[3] =  3;
-  cells[0].material_id = 0;
-
-  cells[1].vertices[0] =  1;
-  cells[1].vertices[1] =  5;
-  cells[1].vertices[2] =  4;
-  cells[1].vertices[3] =  0;
-  cells[1].material_id = 0;
-
-  cells[2].vertices[0] =  2;
-  cells[2].vertices[1] =  6;
-  cells[2].vertices[2] =  5;
-  cells[2].vertices[3] =  1;
-  cells[2].material_id = 0;
-
-  cells[3].vertices[0] =  3;
-  cells[3].vertices[1] =  7;
-  cells[3].vertices[2] =  6;
-  cells[3].vertices[3] =  2;
-  cells[3].material_id = 0;
-
-  cells[4].vertices[0] =  4;
-  cells[4].vertices[1] =  8;
-  cells[4].vertices[2] =  11;
-  cells[4].vertices[3] =  7;
-  cells[4].material_id = 0;
-
-  cells[5].vertices[0] =  5;
-  cells[5].vertices[1] =  9;
-  cells[5].vertices[2] =  8;
-  cells[5].vertices[3] =  4;
-  cells[5].material_id = 0;
-
-  cells[6].vertices[0] =  6;
-  cells[6].vertices[1] =  10;
-  cells[6].vertices[2] =  9;
-  cells[6].vertices[3] =  5;
-  cells[6].material_id = 0;
-
-  cells[7].vertices[0] =  7;
-  cells[7].vertices[1] =  11;
-  cells[7].vertices[2] =  10;
-  cells[7].vertices[3] =  6;
-  cells[7].material_id = 0;
-
-  cells[8].vertices[0] =  8;
-  cells[8].vertices[1] =  12;
-  cells[8].vertices[2] =  15;
-  cells[8].vertices[3] =  11;
-  cells[8].material_id = 0;
-
-  cells[9].vertices[0] =  9;
-  cells[9].vertices[1] =  13;
-  cells[9].vertices[2] =  12;
-  cells[9].vertices[3] =  8;
-  cells[9].material_id = 0;
-
-  cells[10].vertices[0] =  10;
-  cells[10].vertices[1] =  14;
-  cells[10].vertices[2] =  13;
-  cells[10].vertices[3] =  9;
-  cells[10].material_id = 0;
-
-  cells[11].vertices[0] =  11;
-  cells[11].vertices[1] =  15;
-  cells[11].vertices[2] =  14;
-  cells[11].vertices[3] =  10;
-  cells[11].material_id = 0;
-
-  cells[12].vertices[0] =  12;
-  cells[12].vertices[1] =  0;
-  cells[12].vertices[2] =  3;
-  cells[12].vertices[3] =  15;
-  cells[12].material_id = 0;
-
-  cells[13].vertices[0] =  13;
-  cells[13].vertices[1] =  1;
-  cells[13].vertices[2] =  0;
-  cells[13].vertices[3] =  12;
-  cells[13].material_id = 0;
-
-  cells[14].vertices[0] =  14;
-  cells[14].vertices[1] =  2;
-  cells[14].vertices[2] =  1;
-  cells[14].vertices[3] =  13;
-  cells[14].material_id = 0;
-
-  cells[15].vertices[0] =  15;
-  cells[15].vertices[1] =  3;
-  cells[15].vertices[2] =  2;
-  cells[15].vertices[3] =  14;
-  cells[15].material_id = 0;
-
-  // Must call this to be able to create a
-  // correct triangulation in dealii, read
-  // GridReordering<> doc
-  GridReordering<dim,spacedim>::reorder_cells (cells);
-  tria.create_triangulation_compatibility (vertices, cells, SubCellData());
-}
+      for (; cell!=tria.end(); ++cell)
+        for (unsigned int f=0; f<GeometryInfo<3>::faces_per_cell; ++f)
+          {
+            if (!cell->face(f)->at_boundary())
+              continue;
 
+            double radius = cell->face(f)->center().norm() - center.norm();
+            if (std::fabs(cell->face(f)->center()(0)) < eps ) // x = 0 set boundary 2
+              {
+                cell->face(f)->set_boundary_indicator(2);
+                for (unsigned int j=0; j<GeometryInfo<3>::lines_per_face; ++j)
+                  if (cell->face(f)->line(j)->at_boundary())
+                    if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) > eps)
+                      cell->face(f)->line(j)->set_boundary_indicator(2);
+              }
+            else if (std::fabs(cell->face(f)->center()(1)) < eps) // y = 0 set boundary 3
+              {
+                cell->face(f)->set_boundary_indicator(3);
+                for (unsigned int j=0; j<GeometryInfo<3>::lines_per_face; ++j)
+                  if (cell->face(f)->line(j)->at_boundary())
+                    if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) > eps)
+                      cell->face(f)->line(j)->set_boundary_indicator(3);
+              }
+            else if (std::fabs(cell->face(f)->center()(2)) < eps ) // z = 0 set boundary 4
+              {
+                cell->face(f)->set_boundary_indicator(4);
+                for (unsigned int j=0; j<GeometryInfo<3>::lines_per_face; ++j)
+                  if (cell->face(f)->line(j)->at_boundary())
+                    if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) > eps)
+                      cell->face(f)->line(j)->set_boundary_indicator(4);
+              }
+            else if (radius < middle) // inner radius set boundary 0
+              {
+                cell->face(f)->set_boundary_indicator(0);
+                for (unsigned int j=0; j<GeometryInfo<3>::lines_per_face; ++j)
+                  if (cell->face(f)->line(j)->at_boundary())
+                    if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) < eps)
+                      cell->face(f)->line(j)->set_boundary_indicator(0);
+              }
+            else if (radius > middle) // outer radius set boundary 1
+              {
+                cell->face(f)->set_boundary_indicator(1);
+                for (unsigned int j=0; j<GeometryInfo<3>::lines_per_face; ++j)
+                  if (cell->face(f)->line(j)->at_boundary())
+                    if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) < eps)
+                      cell->face(f)->line(j)->set_boundary_indicator(1);
+              }
+            else
+              AssertThrow (false, ExcInternalError());
+          }
+    }
 
-// Implementation for 2D only
-template<>
-void
-GridGenerator::parallelogram (Triangulation<2>  &tria,
-                              const Point<2> (&corners)[2],
-                              const bool         colorize)
-{
-  std::vector<Point<2> > vertices (GeometryInfo<2>::vertices_per_cell);
-
-  vertices[1] = corners[0];
-  vertices[2] = corners[1];
-  vertices[3] = vertices[1] + vertices[2];
-  // Prepare cell data
-  std::vector<CellData<2> > cells (1);
-  for (unsigned int i=0; i<GeometryInfo<2>::vertices_per_cell; ++i)
-    cells[0].vertices[i] = i;
-  cells[0].material_id = 0;
-
-  tria.create_triangulation (vertices, cells, SubCellData());
-
-  // Assign boundary indicators
-  if (colorize)
-    colorize_hyper_rectangle (tria);
-}
+  }
 
 
-template<>
-void
-GridGenerator::parallelogram (Triangulation<2>  &tria,
-                              const Tensor<2,2> &corners,
-                              const bool         colorize)
-{
-  // simply pass everything to the other function of same name
-  const Point<2> x[2] = { corners[0], corners[1] };
-  parallelogram (tria, x, colorize);
-}
+  template <int dim, int spacedim>
+  void
+  hyper_rectangle (Triangulation<dim,spacedim> &tria,
+                   const Point<spacedim>   &p_1,
+                   const Point<spacedim>   &p_2,
+                   const bool          colorize)
+  {
+    // First, normalize input such that
+    // p1 is lower in all coordinate directions.
+    Point<spacedim> p1(p_1);
+    Point<spacedim> p2(p_2);
 
+    for (unsigned int i=0; i<spacedim; ++i)
+      if (p1(i) > p2(i))
+        std::swap (p1(i), p2(i));
 
+    std::vector<Point<spacedim> > vertices (GeometryInfo<dim>::vertices_per_cell);
+    switch (dim)
+      {
+      case 1:
+        vertices[0] = p1;
+        vertices[1] = p2;
+        break;
+      case 2:
+        vertices[0] = vertices[1] = p1;
+        vertices[2] = vertices[3] = p2;
+
+        vertices[1](0) = p2(0);
+        vertices[2](0) = p1(0);
+        break;
+      case 3:
+        vertices[0] = vertices[1] = vertices[2] = vertices[3] = p1;
+        vertices[4] = vertices[5] = vertices[6] = vertices[7] = p2;
+
+        vertices[1](0) = p2(0);
+        vertices[2](1) = p2(1);
+        vertices[3](0) = p2(0);
+        vertices[3](1) = p2(1);
+
+        vertices[4](0) = p1(0);
+        vertices[4](1) = p1(1);
+        vertices[5](1) = p1(1);
+        vertices[6](0) = p1(0);
+
+        break;
+      default:
+        Assert (false, ExcNotImplemented ());
+      }
 
-// Parallelepiped implementation in 1d, 2d, and 3d. @note The
-// implementation in 1d is similar to hyper_rectangle(), and in 2d is
-// similar to parallelogram().
-//
-// The GridReordering::reorder_grid is made use of towards the end of
-// this function. Thus the triangulation is explicitly constructed for
-// all dim here since it is slightly different in that respect
-// (cf. hyper_rectangle(), parallelogram()).
-template<int dim>
-void
-GridGenerator::parallelepiped (Triangulation<dim>  &tria,
-                               const Point<dim>   (&corners) [dim],
-                               const bool           colorize)
-{
-  // Check that none of the user defined vertices overlap
-  for (unsigned int i=0; i<dim; ++i)
-    for (unsigned int j=i+1; j<dim; ++j)
-      Assert ((corners[i]!=corners[j]),
-              ExcMessage ("Invalid distance between corner points of parallelepiped."));
+    // Prepare cell data
+    std::vector<CellData<dim> > cells (1);
+    for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+      cells[0].vertices[i] = i;
+    cells[0].material_id = 0;
 
-  // Note: vertex[0] is the origin and is initialised as so here:
-  std::vector<Point<dim> > vertices (GeometryInfo<dim>::vertices_per_cell);
+    tria.create_triangulation (vertices, cells, SubCellData());
 
-  switch (dim)
-    {
-      // A line (1d parallelepiped)
-    case 1:
-      vertices[1] = corners[0];
-      break;
+    // Assign boundary indicators
+    if (colorize)
+      colorize_hyper_rectangle (tria);
+  }
 
-      // A parallelogram (2d parallelepiped)
-    case 2:
-      // assign corners to vertices:
-      vertices[1] = corners[0];
-      vertices[2] = corners[1];
 
-      // compose the remaining vertex:
-      vertices[3] = vertices[1] + vertices[2];
-      break;
+  template <int dim, int spacedim>
+  void hyper_cube (Triangulation<dim,spacedim> &tria,
+                   const double                 left,
+                   const double                 right)
+  {
+    Assert (left < right,
+            ExcMessage ("Invalid left-to-right bounds of hypercube"));
 
-      // A parallelepiped (3d parallelepiped)
-    case 3:
-      // assign corners to vertices:
-      vertices[1] = corners[0];
-      vertices[2] = corners[1];
-      vertices[4] = corners[2];
-
-      // compose the remaining vertices:
-      vertices[3] = vertices[1] + vertices[2];
-      vertices[5] = vertices[1] + vertices[4];
-      vertices[6] = vertices[2] + vertices[4];
-      vertices[7] = vertices[1] + vertices[2] + vertices[4];
-      break;
+    Point<spacedim> p1;
+    Point<spacedim> p2;
 
-    default:
-      Assert (false, ExcNotImplemented());
-    }
+    p1(spacedim-1) = 0;
+    p2(spacedim-1) = 0;
 
-  // Prepare cell data and wipe material identity
-  std::vector<CellData<dim> > cells (1);
-  for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
-    cells[0].vertices[i] = i;
-  cells[0].material_id = 0;
+    for (unsigned int i=0; i<dim; ++i)
+      {
+        p1(i) = left;
+        p2(i) = right;
+      }
+    hyper_rectangle (tria, p1, p2);
+  }
 
-  // Check ordering of vertices and create triangulation
-  GridReordering<dim>::reorder_cells (cells);
-  tria.create_triangulation (vertices, cells, SubCellData());
 
-  // Finally assign boundary indicators according to hyper_rectangle
-  if (colorize)
-    colorize_hyper_rectangle (tria);
-}
 
-template<int dim>
-void
-GridGenerator::subdivided_parallelepiped (Triangulation<dim>  &tria,
-                                          const size_type      n_subdivisions,
-                                          const Point<dim>   (&corners) [dim],
-                                          const bool           colorize)
-{
-  // Equalise number of subdivisions in each dim-direction, heir
-  // validity will be checked later
-  size_type (n_subdivisions_) [dim];
-  for (unsigned int i=0; i<dim; ++i)
-    n_subdivisions_[i] = n_subdivisions;
-
-  // and call the function below
-  GridGenerator::subdivided_parallelepiped (tria, n_subdivisions_,
-                                            corners,
-                                            colorize);
-}
+  void
+  moebius (Triangulation<3>  &tria,
+           const unsigned int      n_cells,
+           const unsigned int   n_rotations,
+           const double         R,
+           const double         r)
+  {
+    const unsigned int dim=3;
+    Assert (n_cells>4, ExcMessage("More than 4 cells are needed to create a moebius grid."));
+    Assert (r>0 && R>0, ExcMessage("Outer and inner radius must be positive."));
+    Assert (R>r, ExcMessage("Outer radius must be greater than inner radius."));
 
-template<int dim>
-void
-GridGenerator::subdivided_parallelepiped (Triangulation<dim>  &tria,
-                                          const size_type    ( n_subdivisions) [dim],
-                                          const Point<dim>   (&corners) [dim],
-                                          const bool           colorize)
-{
-  // Zero n_subdivisions is the origin only, which makes no sense
-  for (unsigned int i=0; i<dim; ++i)
-    Assert (n_subdivisions[i]>0, ExcInvalidRepetitions(n_subdivisions[i]));
 
-  // Check corners do not overlap (unique)
-  for (unsigned int i=0; i<dim; ++i)
-    for (unsigned int j=i+1; j<dim; ++j)
-      Assert ((corners[i]!=corners[j]),
-              ExcMessage ("Invalid distance between corner points of parallelepiped."));
+    std::vector<Point<dim> > vertices (4*n_cells);
+    double beta_step=n_rotations*numbers::PI/2.0/n_cells;
+    double alpha_step=2.0*numbers::PI/n_cells;
 
-  // Create a list of points
-  Point<dim> (delta) [dim];
+    for (unsigned int i=0; i<n_cells; ++i)
+      for (unsigned int j=0; j<4; ++j)
+        {
+          vertices[4*i+j][0]=R*std::cos(i*alpha_step)+r*std::cos(i*beta_step+j*numbers::PI/2.0)*std::cos(i*alpha_step);
+          vertices[4*i+j][1]=R*std::sin(i*alpha_step)+r*std::cos(i*beta_step+j*numbers::PI/2.0)*std::sin(i*alpha_step);
+          vertices[4*i+j][2]=r*std::sin(i*beta_step+j*numbers::PI/2.0);
+        }
 
-  for (unsigned int i=0; i<dim; ++i)
-    delta[i] = corners[i]/n_subdivisions[i];
-  std::vector<Point<dim> > points;
+    unsigned int offset=0;
 
-  switch (dim)
-    {
-    case 1:
-      for (size_type x=0; x<=n_subdivisions[0]; ++x)
-        points.push_back (Point<dim> (x*delta[0]));
-      break;
+    std::vector<CellData<dim> > cells (n_cells);
+    for (unsigned int i=0; i<n_cells; ++i)
+      {
+        for (unsigned int j=0; j<2; ++j)
+          {
+            cells[i].vertices[0+4*j]=offset+0+4*j;
+            cells[i].vertices[1+4*j]=offset+3+4*j;
+            cells[i].vertices[2+4*j]=offset+2+4*j;
+            cells[i].vertices[3+4*j]=offset+1+4*j;
+          }
+        offset+=4;
+        cells[i].material_id=0;
+      }
 
-    case 2:
-      for (size_type y=0; y<=n_subdivisions[1]; ++y)
-        for (size_type x=0; x<=n_subdivisions[0]; ++x)
-          points.push_back (Point<dim> (x*delta[0] + y*delta[1]));
-      break;
+    // now correct the last four vertices
+    cells[n_cells-1].vertices[4]=(0+n_rotations)%4;
+    cells[n_cells-1].vertices[5]=(3+n_rotations)%4;
+    cells[n_cells-1].vertices[6]=(2+n_rotations)%4;
+    cells[n_cells-1].vertices[7]=(1+n_rotations)%4;
 
-    case 3:
-      for (size_type z=0; z<=n_subdivisions[2]; ++z)
-        for (size_type y=0; y<=n_subdivisions[1]; ++y)
-          for (size_type x=0; x<=n_subdivisions[0]; ++x)
-            points.push_back (Point<dim> (x*delta[0] + y*delta[1] + z*delta[2]));
-      break;
+    GridReordering<dim>::invert_all_cells_of_negative_grid(vertices,cells);
+    tria.create_triangulation_compatibility (vertices, cells, SubCellData());
+  }
 
-    default:
-      Assert (false, ExcNotImplemented());
-    }
 
-  // Prepare cell data
-  size_type n_cells = 1;
-  for (unsigned int i=0; i<dim; ++i)
-    n_cells *= n_subdivisions[i];
-  std::vector<CellData<dim> > cells (n_cells);
 
-  // Create fixed ordering of
-  switch (dim)
-    {
-    case 1:
-      for (size_type x=0; x<n_subdivisions[0]; ++x)
-        {
-          cells[x].vertices[0] = x;
-          cells[x].vertices[1] = x+1;
+  void
+  torus (Triangulation<2,3>  &tria,
+         const double         R,
+         const double         r)
+  {
+    Assert (R>r, ExcMessage("Outer radius must be greater than inner radius."));
+
+    const unsigned int dim=2;
+    const unsigned int spacedim=3;
+    std::vector<Point<spacedim> > vertices (16);
+
+    vertices[0]=Point<spacedim>(R-r,0,0);
+    vertices[1]=Point<spacedim>(R,-r,0);
+    vertices[2]=Point<spacedim>(R+r,0,0);
+    vertices[3]=Point<spacedim>(R, r,0);
+    vertices[4]=Point<spacedim>(0,0,R-r);
+    vertices[5]=Point<spacedim>(0,-r,R);
+    vertices[6]=Point<spacedim>(0,0,R+r);
+    vertices[7]=Point<spacedim>(0,r,R);
+    vertices[8]=Point<spacedim>(-(R-r),0,0);
+    vertices[9]=Point<spacedim>(-R,-r,0);
+    vertices[10]=Point<spacedim>(-(R+r),0,0);
+    vertices[11]=Point<spacedim>(-R, r,0);
+    vertices[12]=Point<spacedim>(0,0,-(R-r));
+    vertices[13]=Point<spacedim>(0,-r,-R);
+    vertices[14]=Point<spacedim>(0,0,-(R+r));
+    vertices[15]=Point<spacedim>(0,r,-R);
+
+    std::vector<CellData<dim> > cells (16);
+    //Right Hand Orientation
+    cells[0].vertices[0] =  0;
+    cells[0].vertices[1] =  4;
+    cells[0].vertices[2] =  7;
+    cells[0].vertices[3] =  3;
+    cells[0].material_id = 0;
+
+    cells[1].vertices[0] =  1;
+    cells[1].vertices[1] =  5;
+    cells[1].vertices[2] =  4;
+    cells[1].vertices[3] =  0;
+    cells[1].material_id = 0;
+
+    cells[2].vertices[0] =  2;
+    cells[2].vertices[1] =  6;
+    cells[2].vertices[2] =  5;
+    cells[2].vertices[3] =  1;
+    cells[2].material_id = 0;
+
+    cells[3].vertices[0] =  3;
+    cells[3].vertices[1] =  7;
+    cells[3].vertices[2] =  6;
+    cells[3].vertices[3] =  2;
+    cells[3].material_id = 0;
+
+    cells[4].vertices[0] =  4;
+    cells[4].vertices[1] =  8;
+    cells[4].vertices[2] =  11;
+    cells[4].vertices[3] =  7;
+    cells[4].material_id = 0;
+
+    cells[5].vertices[0] =  5;
+    cells[5].vertices[1] =  9;
+    cells[5].vertices[2] =  8;
+    cells[5].vertices[3] =  4;
+    cells[5].material_id = 0;
+
+    cells[6].vertices[0] =  6;
+    cells[6].vertices[1] =  10;
+    cells[6].vertices[2] =  9;
+    cells[6].vertices[3] =  5;
+    cells[6].material_id = 0;
+
+    cells[7].vertices[0] =  7;
+    cells[7].vertices[1] =  11;
+    cells[7].vertices[2] =  10;
+    cells[7].vertices[3] =  6;
+    cells[7].material_id = 0;
+
+    cells[8].vertices[0] =  8;
+    cells[8].vertices[1] =  12;
+    cells[8].vertices[2] =  15;
+    cells[8].vertices[3] =  11;
+    cells[8].material_id = 0;
+
+    cells[9].vertices[0] =  9;
+    cells[9].vertices[1] =  13;
+    cells[9].vertices[2] =  12;
+    cells[9].vertices[3] =  8;
+    cells[9].material_id = 0;
+
+    cells[10].vertices[0] =  10;
+    cells[10].vertices[1] =  14;
+    cells[10].vertices[2] =  13;
+    cells[10].vertices[3] =  9;
+    cells[10].material_id = 0;
+
+    cells[11].vertices[0] =  11;
+    cells[11].vertices[1] =  15;
+    cells[11].vertices[2] =  14;
+    cells[11].vertices[3] =  10;
+    cells[11].material_id = 0;
+
+    cells[12].vertices[0] =  12;
+    cells[12].vertices[1] =  0;
+    cells[12].vertices[2] =  3;
+    cells[12].vertices[3] =  15;
+    cells[12].material_id = 0;
+
+    cells[13].vertices[0] =  13;
+    cells[13].vertices[1] =  1;
+    cells[13].vertices[2] =  0;
+    cells[13].vertices[3] =  12;
+    cells[13].material_id = 0;
+
+    cells[14].vertices[0] =  14;
+    cells[14].vertices[1] =  2;
+    cells[14].vertices[2] =  1;
+    cells[14].vertices[3] =  13;
+    cells[14].material_id = 0;
+
+    cells[15].vertices[0] =  15;
+    cells[15].vertices[1] =  3;
+    cells[15].vertices[2] =  2;
+    cells[15].vertices[3] =  14;
+    cells[15].material_id = 0;
+
+    // Must call this to be able to create a
+    // correct triangulation in dealii, read
+    // GridReordering<> doc
+    GridReordering<dim,spacedim>::reorder_cells (cells);
+    tria.create_triangulation_compatibility (vertices, cells, SubCellData());
+  }
 
-          // wipe material id
-          cells[x].material_id = 0;
-        }
-      break;
 
-    case 2:
-    {
-      // Shorthand
-      const size_type n_dy = n_subdivisions[1];
-      const size_type n_dx = n_subdivisions[0];
+// Implementation for 2D only
+  template<>
+  void
+  parallelogram (Triangulation<2>  &tria,
+                 const Point<2> (&corners)[2],
+                 const bool         colorize)
+  {
+    std::vector<Point<2> > vertices (GeometryInfo<2>::vertices_per_cell);
+
+    vertices[1] = corners[0];
+    vertices[2] = corners[1];
+    vertices[3] = vertices[1] + vertices[2];
+    // Prepare cell data
+    std::vector<CellData<2> > cells (1);
+    for (unsigned int i=0; i<GeometryInfo<2>::vertices_per_cell; ++i)
+      cells[0].vertices[i] = i;
+    cells[0].material_id = 0;
+
+    tria.create_triangulation (vertices, cells, SubCellData());
+
+    // Assign boundary indicators
+    if (colorize)
+      colorize_hyper_rectangle (tria);
+  }
+
+
+  template<>
+  void
+  parallelogram (Triangulation<2>  &tria,
+                 const Tensor<2,2> &corners,
+                 const bool         colorize)
+  {
+    // simply pass everything to the other function of same name
+    const Point<2> x[2] = { corners[0], corners[1] };
+    parallelogram (tria, x, colorize);
+  }
+
+
+
+// Parallelepiped implementation in 1d, 2d, and 3d. @note The
+// implementation in 1d is similar to hyper_rectangle(), and in 2d is
+// similar to parallelogram().
+//
+// The GridReordering::reorder_grid is made use of towards the end of
+// this function. Thus the triangulation is explicitly constructed for
+// all dim here since it is slightly different in that respect
+// (cf. hyper_rectangle(), parallelogram()).
+  template<int dim>
+  void
+  parallelepiped (Triangulation<dim>  &tria,
+                  const Point<dim>   (&corners) [dim],
+                  const bool           colorize)
+  {
+    // Check that none of the user defined vertices overlap
+    for (unsigned int i=0; i<dim; ++i)
+      for (unsigned int j=i+1; j<dim; ++j)
+        Assert ((corners[i]!=corners[j]),
+                ExcMessage ("Invalid distance between corner points of parallelepiped."));
+
+    // Note: vertex[0] is the origin and is initialised as so here:
+    std::vector<Point<dim> > vertices (GeometryInfo<dim>::vertices_per_cell);
+
+    switch (dim)
+      {
+        // A line (1d parallelepiped)
+      case 1:
+        vertices[1] = corners[0];
+        break;
+
+        // A parallelogram (2d parallelepiped)
+      case 2:
+        // assign corners to vertices:
+        vertices[1] = corners[0];
+        vertices[2] = corners[1];
+
+        // compose the remaining vertex:
+        vertices[3] = vertices[1] + vertices[2];
+        break;
+
+        // A parallelepiped (3d parallelepiped)
+      case 3:
+        // assign corners to vertices:
+        vertices[1] = corners[0];
+        vertices[2] = corners[1];
+        vertices[4] = corners[2];
+
+        // compose the remaining vertices:
+        vertices[3] = vertices[1] + vertices[2];
+        vertices[5] = vertices[1] + vertices[4];
+        vertices[6] = vertices[2] + vertices[4];
+        vertices[7] = vertices[1] + vertices[2] + vertices[4];
+        break;
+
+      default:
+        Assert (false, ExcNotImplemented());
+      }
+
+    // Prepare cell data and wipe material identity
+    std::vector<CellData<dim> > cells (1);
+    for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+      cells[0].vertices[i] = i;
+    cells[0].material_id = 0;
+
+    // Check ordering of vertices and create triangulation
+    GridReordering<dim>::reorder_cells (cells);
+    tria.create_triangulation (vertices, cells, SubCellData());
+
+    // Finally assign boundary indicators according to hyper_rectangle
+    if (colorize)
+      colorize_hyper_rectangle (tria);
+  }
+
+  template<int dim>
+  void
+  subdivided_parallelepiped (Triangulation<dim>  &tria,
+                             const unsigned int      n_subdivisions,
+                             const Point<dim>   (&corners) [dim],
+                             const bool           colorize)
+  {
+    // Equalise number of subdivisions in each dim-direction, heir
+    // validity will be checked later
+    unsigned int (n_subdivisions_) [dim];
+    for (unsigned int i=0; i<dim; ++i)
+      n_subdivisions_[i] = n_subdivisions;
+
+    // and call the function below
+    subdivided_parallelepiped (tria, n_subdivisions_,
+                               corners,
+                               colorize);
+  }
+
+  template<int dim>
+  void
+  subdivided_parallelepiped (Triangulation<dim>  &tria,
+                             const unsigned int    ( n_subdivisions) [dim],
+                             const Point<dim>   (&corners) [dim],
+                             const bool           colorize)
+  {
+    // Zero n_subdivisions is the origin only, which makes no sense
+    for (unsigned int i=0; i<dim; ++i)
+      Assert (n_subdivisions[i]>0, ExcInvalidRepetitions(n_subdivisions[i]));
 
-      for (size_type y=0; y<n_dy; ++y)
-        for (size_type x=0; x<n_dx; ++x)
+    // Check corners do not overlap (unique)
+    for (unsigned int i=0; i<dim; ++i)
+      for (unsigned int j=i+1; j<dim; ++j)
+        Assert ((corners[i]!=corners[j]),
+                ExcMessage ("Invalid distance between corner points of parallelepiped."));
+
+    // Create a list of points
+    Point<dim> (delta) [dim];
+
+    for (unsigned int i=0; i<dim; ++i)
+      delta[i] = corners[i]/n_subdivisions[i];
+    std::vector<Point<dim> > points;
+
+    switch (dim)
+      {
+      case 1:
+        for (unsigned int x=0; x<=n_subdivisions[0]; ++x)
+          points.push_back (Point<dim> (x*delta[0]));
+        break;
+
+      case 2:
+        for (unsigned int y=0; y<=n_subdivisions[1]; ++y)
+          for (unsigned int x=0; x<=n_subdivisions[0]; ++x)
+            points.push_back (Point<dim> (x*delta[0] + y*delta[1]));
+        break;
+
+      case 3:
+        for (unsigned int z=0; z<=n_subdivisions[2]; ++z)
+          for (unsigned int y=0; y<=n_subdivisions[1]; ++y)
+            for (unsigned int x=0; x<=n_subdivisions[0]; ++x)
+              points.push_back (Point<dim> (x*delta[0] + y*delta[1] + z*delta[2]));
+        break;
+
+      default:
+        Assert (false, ExcNotImplemented());
+      }
+
+    // Prepare cell data
+    unsigned int n_cells = 1;
+    for (unsigned int i=0; i<dim; ++i)
+      n_cells *= n_subdivisions[i];
+    std::vector<CellData<dim> > cells (n_cells);
+
+    // Create fixed ordering of
+    switch (dim)
+      {
+      case 1:
+        for (unsigned int x=0; x<n_subdivisions[0]; ++x)
           {
-            const size_type    c = y*n_dx         + x;
-            cells[c].vertices[0] = y*(n_dx+1)     + x;
-            cells[c].vertices[1] = y*(n_dx+1)     + x+1;
-            cells[c].vertices[2] = (y+1)*(n_dx+1) + x;
-            cells[c].vertices[3] = (y+1)*(n_dx+1) + x+1;
+            cells[x].vertices[0] = x;
+            cells[x].vertices[1] = x+1;
 
             // wipe material id
-            cells[c].material_id = 0;
+            cells[x].material_id = 0;
           }
-    }
-    break;
+        break;
 
-    case 3:
-    {
-      // Shorthand
-      const size_type n_dz = n_subdivisions[2];
-      const size_type n_dy = n_subdivisions[1];
-      const size_type n_dx = n_subdivisions[0];
-
-      for (size_type z=0; z<n_dz; ++z)
-        for (size_type y=0; y<n_dy; ++y)
-          for (size_type x=0; x<n_dx; ++x)
-            {
-              const size_type    c = z*n_dy*n_dx             + y*n_dx         + x;
+      case 2:
+      {
+        // Shorthand
+        const unsigned int n_dy = n_subdivisions[1];
+        const unsigned int n_dx = n_subdivisions[0];
 
-              cells[c].vertices[0] = z*(n_dy+1)*(n_dx+1)     + y*(n_dx+1)     + x;
-              cells[c].vertices[1] = z*(n_dy+1)*(n_dx+1)     + y*(n_dx+1)     + x+1;
-              cells[c].vertices[2] = z*(n_dy+1)*(n_dx+1)     + (y+1)*(n_dx+1) + x;
-              cells[c].vertices[3] = z*(n_dy+1)*(n_dx+1)     + (y+1)*(n_dx+1) + x+1;
-              cells[c].vertices[4] = (z+1)*(n_dy+1)*(n_dx+1) + y*(n_dx+1)     + x;
-              cells[c].vertices[5] = (z+1)*(n_dy+1)*(n_dx+1) + y*(n_dx+1)     + x+1;
-              cells[c].vertices[6] = (z+1)*(n_dy+1)*(n_dx+1) + (y+1)*(n_dx+1) + x;
-              cells[c].vertices[7] = (z+1)*(n_dy+1)*(n_dx+1) + (y+1)*(n_dx+1) + x+1;
+        for (unsigned int y=0; y<n_dy; ++y)
+          for (unsigned int x=0; x<n_dx; ++x)
+            {
+              const unsigned int    c = y*n_dx         + x;
+              cells[c].vertices[0] = y*(n_dx+1)     + x;
+              cells[c].vertices[1] = y*(n_dx+1)     + x+1;
+              cells[c].vertices[2] = (y+1)*(n_dx+1) + x;
+              cells[c].vertices[3] = (y+1)*(n_dx+1) + x+1;
 
               // wipe material id
               cells[c].material_id = 0;
             }
+      }
       break;
-    }
 
-    default:
-      Assert (false, ExcNotImplemented());
-    }
+      case 3:
+      {
+        // Shorthand
+        const unsigned int n_dz = n_subdivisions[2];
+        const unsigned int n_dy = n_subdivisions[1];
+        const unsigned int n_dx = n_subdivisions[0];
+
+        for (unsigned int z=0; z<n_dz; ++z)
+          for (unsigned int y=0; y<n_dy; ++y)
+            for (unsigned int x=0; x<n_dx; ++x)
+              {
+                const unsigned int    c = z*n_dy*n_dx             + y*n_dx         + x;
+
+                cells[c].vertices[0] = z*(n_dy+1)*(n_dx+1)     + y*(n_dx+1)     + x;
+                cells[c].vertices[1] = z*(n_dy+1)*(n_dx+1)     + y*(n_dx+1)     + x+1;
+                cells[c].vertices[2] = z*(n_dy+1)*(n_dx+1)     + (y+1)*(n_dx+1) + x;
+                cells[c].vertices[3] = z*(n_dy+1)*(n_dx+1)     + (y+1)*(n_dx+1) + x+1;
+                cells[c].vertices[4] = (z+1)*(n_dy+1)*(n_dx+1) + y*(n_dx+1)     + x;
+                cells[c].vertices[5] = (z+1)*(n_dy+1)*(n_dx+1) + y*(n_dx+1)     + x+1;
+                cells[c].vertices[6] = (z+1)*(n_dy+1)*(n_dx+1) + (y+1)*(n_dx+1) + x;
+                cells[c].vertices[7] = (z+1)*(n_dy+1)*(n_dx+1) + (y+1)*(n_dx+1) + x+1;
+
+                // wipe material id
+                cells[c].material_id = 0;
+              }
+        break;
+      }
 
-  // Create triangulation
-  tria.create_triangulation (points, cells, SubCellData());
+      default:
+        Assert (false, ExcNotImplemented());
+      }
 
-  // Finally assign boundary indicators according to hyper_rectangle
-  if (colorize)
-    colorize_hyper_rectangle (tria);
-}
+    // Create triangulation
+    tria.create_triangulation (points, cells, SubCellData());
 
+    // Finally assign boundary indicators according to hyper_rectangle
+    if (colorize)
+      colorize_hyper_rectangle (tria);
+  }
 
-template <int dim>
-void
-GridGenerator::subdivided_hyper_cube (Triangulation<dim> &tria,
-                                      const unsigned int  repetitions,
-                                      const double        left,
-                                      const double        right)
-{
-  Assert (repetitions >= 1, ExcInvalidRepetitions(repetitions));
-  Assert (left < right,
-          ExcMessage ("Invalid left-to-right bounds of hypercube"));
-
-  // first generate the necessary
-  // points
-  const double delta = (right-left)/repetitions;
-  std::vector<Point<dim> > points;
-  switch (dim)
-    {
-    case 1:
-      for (unsigned int x=0; x<=repetitions; ++x)
-        points.push_back (Point<dim> (left+x*delta));
-      break;
 
-    case 2:
-      for (unsigned int y=0; y<=repetitions; ++y)
+  template <int dim>
+  void
+  subdivided_hyper_cube (Triangulation<dim> &tria,
+                         const unsigned int  repetitions,
+                         const double        left,
+                         const double        right)
+  {
+    Assert (repetitions >= 1, ExcInvalidRepetitions(repetitions));
+    Assert (left < right,
+            ExcMessage ("Invalid left-to-right bounds of hypercube"));
+
+    // first generate the necessary
+    // points
+    const double delta = (right-left)/repetitions;
+    std::vector<Point<dim> > points;
+    switch (dim)
+      {
+      case 1:
         for (unsigned int x=0; x<=repetitions; ++x)
-          points.push_back (Point<dim> (left+x*delta,
-                                        left+y*delta));
-      break;
+          points.push_back (Point<dim> (left+x*delta));
+        break;
 
-    case 3:
-      for (unsigned int z=0; z<=repetitions; ++z)
+      case 2:
         for (unsigned int y=0; y<=repetitions; ++y)
           for (unsigned int x=0; x<=repetitions; ++x)
             points.push_back (Point<dim> (left+x*delta,
-                                          left+y*delta,
-                                          left+z*delta));
-      break;
-
-    default:
-      Assert (false, ExcNotImplemented());
-    }
-
-  // next create the cells
-  // Prepare cell data
-  std::vector<CellData<dim> > cells;
-  // Define these as abbreviations
-  // for the step sizes below. The
-  // number of points in a single
-  // direction is repetitions+1
-  const unsigned int dy = repetitions+1;
-  const unsigned int dz = dy*dy;
-  switch (dim)
-    {
-    case 1:
-      cells.resize (repetitions);
-      for (unsigned int x=0; x<repetitions; ++x)
-        {
-          cells[x].vertices[0] = x;
-          cells[x].vertices[1] = x+1;
-          cells[x].material_id = 0;
-        }
-      break;
+                                          left+y*delta));
+        break;
+
+      case 3:
+        for (unsigned int z=0; z<=repetitions; ++z)
+          for (unsigned int y=0; y<=repetitions; ++y)
+            for (unsigned int x=0; x<=repetitions; ++x)
+              points.push_back (Point<dim> (left+x*delta,
+                                            left+y*delta,
+                                            left+z*delta));
+        break;
+
+      default:
+        Assert (false, ExcNotImplemented());
+      }
 
-    case 2:
-      cells.resize (repetitions*repetitions);
-      for (unsigned int y=0; y<repetitions; ++y)
+    // next create the cells
+    // Prepare cell data
+    std::vector<CellData<dim> > cells;
+    // Define these as abbreviations
+    // for the step sizes below. The
+    // number of points in a single
+    // direction is repetitions+1
+    const unsigned int dy = repetitions+1;
+    const unsigned int dz = dy*dy;
+    switch (dim)
+      {
+      case 1:
+        cells.resize (repetitions);
         for (unsigned int x=0; x<repetitions; ++x)
           {
-            const unsigned int c = x  +y*repetitions;
-            cells[c].vertices[0] = x  +y*dy;
-            cells[c].vertices[1] = x+1+y*dy;
-            cells[c].vertices[2] = x  +(y+1)*dy;
-            cells[c].vertices[3] = x+1+(y+1)*dy;
-            cells[c].material_id = 0;
+            cells[x].vertices[0] = x;
+            cells[x].vertices[1] = x+1;
+            cells[x].material_id = 0;
           }
-      break;
+        break;
 
-    case 3:
-      cells.resize (repetitions*repetitions*repetitions);
-      for (unsigned int z=0; z<repetitions; ++z)
+      case 2:
+        cells.resize (repetitions*repetitions);
         for (unsigned int y=0; y<repetitions; ++y)
           for (unsigned int x=0; x<repetitions; ++x)
             {
-              const unsigned int c = x+y*repetitions
-                                     +z*repetitions*repetitions;
-              cells[c].vertices[0] = x  +y*dy    +z*dz;
-              cells[c].vertices[1] = x+1+y*dy    +z*dz;
-              cells[c].vertices[2] = x  +(y+1)*dy+z*dz;
-              cells[c].vertices[3] = x+1+(y+1)*dy+z*dz;
-              cells[c].vertices[4] = x  +y*dy    +(z+1)*dz;
-              cells[c].vertices[5] = x+1+y*dy    +(z+1)*dz;
-              cells[c].vertices[6] = x  +(y+1)*dy+(z+1)*dz;
-              cells[c].vertices[7] = x+1+(y+1)*dy+(z+1)*dz;
+              const unsigned int c = x  +y*repetitions;
+              cells[c].vertices[0] = x  +y*dy;
+              cells[c].vertices[1] = x+1+y*dy;
+              cells[c].vertices[2] = x  +(y+1)*dy;
+              cells[c].vertices[3] = x+1+(y+1)*dy;
               cells[c].material_id = 0;
             }
-      break;
+        break;
 
-    default:
-      // should be trivial to
-      // do for 3d as well, but
-      // am too tired at this
-      // point of the night to
-      // do that...
-      //
-      // contributions are welcome!
-      Assert (false, ExcNotImplemented());
-    }
+      case 3:
+        cells.resize (repetitions*repetitions*repetitions);
+        for (unsigned int z=0; z<repetitions; ++z)
+          for (unsigned int y=0; y<repetitions; ++y)
+            for (unsigned int x=0; x<repetitions; ++x)
+              {
+                const unsigned int c = x+y*repetitions
+                                       +z*repetitions*repetitions;
+                cells[c].vertices[0] = x  +y*dy    +z*dz;
+                cells[c].vertices[1] = x+1+y*dy    +z*dz;
+                cells[c].vertices[2] = x  +(y+1)*dy+z*dz;
+                cells[c].vertices[3] = x+1+(y+1)*dy+z*dz;
+                cells[c].vertices[4] = x  +y*dy    +(z+1)*dz;
+                cells[c].vertices[5] = x+1+y*dy    +(z+1)*dz;
+                cells[c].vertices[6] = x  +(y+1)*dy+(z+1)*dz;
+                cells[c].vertices[7] = x+1+(y+1)*dy+(z+1)*dz;
+                cells[c].material_id = 0;
+              }
+        break;
+
+      default:
+        // should be trivial to
+        // do for 3d as well, but
+        // am too tired at this
+        // point of the night to
+        // do that...
+        //
+        // contributions are welcome!
+        Assert (false, ExcNotImplemented());
+      }
 
-  tria.create_triangulation (points, cells, SubCellData());
-}
+    tria.create_triangulation (points, cells, SubCellData());
+  }
 
 
 
-template <int dim>
-void
-GridGenerator::subdivided_hyper_rectangle (
-  Triangulation<dim>              &tria,
-  const std::vector<unsigned int> &repetitions,
-  const Point<dim>                &p_1,
-  const Point<dim>                &p_2,
-  const bool                       colorize)
-{
-  // contributed by Joerg R. Weimar
-  // (j.weimar@jweimar.de) 2003
-  Assert(repetitions.size() == dim,
-         ExcInvalidRepetitionsDimension(dim));
-  // First, normalize input such that
-  // p1 is lower in all coordinate
-  // directions.
-  Point<dim> p1(p_1);
-  Point<dim> p2(p_2);
-
-  for (unsigned int i=0; i<dim; ++i)
-    if (p1(i) > p2(i))
-      std::swap (p1(i), p2(i));
-
-  // then check that all repetitions
-  // are >= 1, and calculate deltas
-  // convert repetitions from double
-  // to int by taking the ceiling.
-  Point<dim> delta;
-
-  for (unsigned int i=0; i<dim; ++i)
-    {
-      Assert (repetitions[i] >= 1, ExcInvalidRepetitions(repetitions[i]));
+  template <int dim>
+  void
+  subdivided_hyper_rectangle (
+    Triangulation<dim>              &tria,
+    const std::vector<unsigned int> &repetitions,
+    const Point<dim>                &p_1,
+    const Point<dim>                &p_2,
+    const bool                       colorize)
+  {
+    // contributed by Joerg R. Weimar
+    // (j.weimar@jweimar.de) 2003
+    Assert(repetitions.size() == dim,
+           ExcInvalidRepetitionsDimension(dim));
+    // First, normalize input such that
+    // p1 is lower in all coordinate
+    // directions.
+    Point<dim> p1(p_1);
+    Point<dim> p2(p_2);
+
+    for (unsigned int i=0; i<dim; ++i)
+      if (p1(i) > p2(i))
+        std::swap (p1(i), p2(i));
 
-      delta[i] = (p2[i]-p1[i])/repetitions[i];
-    }
+    // then check that all repetitions
+    // are >= 1, and calculate deltas
+    // convert repetitions from double
+    // to int by taking the ceiling.
+    Point<dim> delta;
 
-  // then generate the necessary
-  // points
-  std::vector<Point<dim> > points;
-  switch (dim)
-    {
-    case 1:
-      for (unsigned int x=0; x<=repetitions[0]; ++x)
-        points.push_back (Point<dim> (p1[0]+x*delta[0]));
-      break;
+    for (unsigned int i=0; i<dim; ++i)
+      {
+        Assert (repetitions[i] >= 1, ExcInvalidRepetitions(repetitions[i]));
 
-    case 2:
-      for (unsigned int y=0; y<=repetitions[1]; ++y)
+        delta[i] = (p2[i]-p1[i])/repetitions[i];
+      }
+
+    // then generate the necessary
+    // points
+    std::vector<Point<dim> > points;
+    switch (dim)
+      {
+      case 1:
         for (unsigned int x=0; x<=repetitions[0]; ++x)
-          points.push_back (Point<dim> (p1[0]+x*delta[0],
-                                        p1[1]+y*delta[1]));
-      break;
+          points.push_back (Point<dim> (p1[0]+x*delta[0]));
+        break;
 
-    case 3:
-      for (unsigned int z=0; z<=repetitions[2]; ++z)
+      case 2:
         for (unsigned int y=0; y<=repetitions[1]; ++y)
           for (unsigned int x=0; x<=repetitions[0]; ++x)
             points.push_back (Point<dim> (p1[0]+x*delta[0],
-                                          p1[1]+y*delta[1],
-                                          p1[2]+z*delta[2]));
-      break;
-
-    default:
-      Assert (false, ExcNotImplemented());
-    }
-
-  // next create the cells
-  // Prepare cell data
-  std::vector<CellData<dim> > cells;
-  switch (dim)
-    {
-    case 1:
-    {
-      cells.resize (repetitions[0]);
-      for (unsigned int x=0; x<repetitions[0]; ++x)
-        {
-          cells[x].vertices[0] = x;
-          cells[x].vertices[1] = x+1;
-          cells[x].material_id = 0;
-        }
-      break;
-    }
+                                          p1[1]+y*delta[1]));
+        break;
+
+      case 3:
+        for (unsigned int z=0; z<=repetitions[2]; ++z)
+          for (unsigned int y=0; y<=repetitions[1]; ++y)
+            for (unsigned int x=0; x<=repetitions[0]; ++x)
+              points.push_back (Point<dim> (p1[0]+x*delta[0],
+                                            p1[1]+y*delta[1],
+                                            p1[2]+z*delta[2]));
+        break;
+
+      default:
+        Assert (false, ExcNotImplemented());
+      }
 
-    case 2:
-    {
-      cells.resize (repetitions[1]*repetitions[0]);
-      for (unsigned int y=0; y<repetitions[1]; ++y)
+    // next create the cells
+    // Prepare cell data
+    std::vector<CellData<dim> > cells;
+    switch (dim)
+      {
+      case 1:
+      {
+        cells.resize (repetitions[0]);
         for (unsigned int x=0; x<repetitions[0]; ++x)
           {
-            const unsigned int c = x+y*repetitions[0];
-            cells[c].vertices[0] = y*(repetitions[0]+1)+x;
-            cells[c].vertices[1] = y*(repetitions[0]+1)+x+1;
-            cells[c].vertices[2] = (y+1)*(repetitions[0]+1)+x;
-            cells[c].vertices[3] = (y+1)*(repetitions[0]+1)+x+1;
-            cells[c].material_id = 0;
+            cells[x].vertices[0] = x;
+            cells[x].vertices[1] = x+1;
+            cells[x].material_id = 0;
           }
-      break;
-    }
-
-    case 3:
-    {
-      const unsigned int n_x  = (repetitions[0]+1);
-      const unsigned int n_xy = (repetitions[0]+1)*(repetitions[1]+1);
+        break;
+      }
 
-      cells.resize (repetitions[2]*repetitions[1]*repetitions[0]);
-      for (unsigned int z=0; z<repetitions[2]; ++z)
+      case 2:
+      {
+        cells.resize (repetitions[1]*repetitions[0]);
         for (unsigned int y=0; y<repetitions[1]; ++y)
           for (unsigned int x=0; x<repetitions[0]; ++x)
             {
-              const unsigned int c = x+y*repetitions[0] +
-                                     z*repetitions[0]*repetitions[1];
-              cells[c].vertices[0] = z*n_xy + y*n_x + x;
-              cells[c].vertices[1] = z*n_xy + y*n_x + x+1;
-              cells[c].vertices[2] = z*n_xy + (y+1)*n_x + x;
-              cells[c].vertices[3] = z*n_xy + (y+1)*n_x + x+1;
-              cells[c].vertices[4] = (z+1)*n_xy + y*n_x + x;
-              cells[c].vertices[5] = (z+1)*n_xy + y*n_x + x+1;
-              cells[c].vertices[6] = (z+1)*n_xy + (y+1)*n_x + x;
-              cells[c].vertices[7] = (z+1)*n_xy + (y+1)*n_x + x+1;
+              const unsigned int c = x+y*repetitions[0];
+              cells[c].vertices[0] = y*(repetitions[0]+1)+x;
+              cells[c].vertices[1] = y*(repetitions[0]+1)+x+1;
+              cells[c].vertices[2] = (y+1)*(repetitions[0]+1)+x;
+              cells[c].vertices[3] = (y+1)*(repetitions[0]+1)+x+1;
               cells[c].material_id = 0;
             }
-      break;
+        break;
+      }
 
-    }
+      case 3:
+      {
+        const unsigned int n_x  = (repetitions[0]+1);
+        const unsigned int n_xy = (repetitions[0]+1)*(repetitions[1]+1);
 
-    default:
-      Assert (false, ExcNotImplemented());
-    }
+        cells.resize (repetitions[2]*repetitions[1]*repetitions[0]);
+        for (unsigned int z=0; z<repetitions[2]; ++z)
+          for (unsigned int y=0; y<repetitions[1]; ++y)
+            for (unsigned int x=0; x<repetitions[0]; ++x)
+              {
+                const unsigned int c = x+y*repetitions[0] +
+                                       z*repetitions[0]*repetitions[1];
+                cells[c].vertices[0] = z*n_xy + y*n_x + x;
+                cells[c].vertices[1] = z*n_xy + y*n_x + x+1;
+                cells[c].vertices[2] = z*n_xy + (y+1)*n_x + x;
+                cells[c].vertices[3] = z*n_xy + (y+1)*n_x + x+1;
+                cells[c].vertices[4] = (z+1)*n_xy + y*n_x + x;
+                cells[c].vertices[5] = (z+1)*n_xy + y*n_x + x+1;
+                cells[c].vertices[6] = (z+1)*n_xy + (y+1)*n_x + x;
+                cells[c].vertices[7] = (z+1)*n_xy + (y+1)*n_x + x+1;
+                cells[c].material_id = 0;
+              }
+        break;
 
-  tria.create_triangulation (points, cells, SubCellData());
+      }
 
-  if (colorize)
-    {
-      // to colorize, run through all
-      // faces of all cells and set
-      // boundary indicator to the
-      // correct value if it was 0.
-
-      // use a large epsilon to
-      // compare numbers to avoid
-      // roundoff problems.
-      const double epsilon
-        = 0.01 * *std::min_element (&delta[0], &delta[0]+dim);
-      Assert (epsilon > 0,
-              ExcMessage ("The distance between corner points must be positive."))
-
-      // actual code is external since
-      // 1-D is different from 2/3D.
-      colorize_subdivided_hyper_rectangle (tria, p1, p2, epsilon);
-    }
-}
+      default:
+        Assert (false, ExcNotImplemented());
+      }
 
+    tria.create_triangulation (points, cells, SubCellData());
 
+    if (colorize)
+      {
+        // to colorize, run through all
+        // faces of all cells and set
+        // boundary indicator to the
+        // correct value if it was 0.
+
+        // use a large epsilon to
+        // compare numbers to avoid
+        // roundoff problems.
+        const double epsilon
+          = 0.01 * *std::min_element (&delta[0], &delta[0]+dim);
+        Assert (epsilon > 0,
+                ExcMessage ("The distance between corner points must be positive."))
+
+        // actual code is external since
+        // 1-D is different from 2/3D.
+        colorize_subdivided_hyper_rectangle (tria, p1, p2, epsilon);
+      }
+  }
 
-template <int dim>
-void
-GridGenerator::subdivided_hyper_rectangle(
-  Triangulation<dim>              &tria,
-  const std::vector<std::vector<double> > &step_sz,
-  const Point<dim>                &p_1,
-  const Point<dim>                &p_2,
-  const bool                       colorize)
-{
-  // contributed by Joerg R. Weimar
-  // (j.weimar@jweimar.de) 2003
-  // modified by Yaqi Wang 2006
-  Assert(step_sz.size() == dim,
-         ExcInvalidRepetitionsDimension(dim));
-
-
-  // First, normalize input such that
-  // p1 is lower in all coordinate
-  // directions.
-
-  // and check the consistency of
-  // step sizes, i.e. that they all
-  // add up to the sizes specified by
-  // p_1 and p_2
-  Point<dim> p1(p_1);
-  Point<dim> p2(p_2);
-  std::vector< std::vector<double> > step_sizes(step_sz);
-
-  for (unsigned int i=0; i<dim; ++i)
-    {
-      if (p1(i) > p2(i))
-        {
-          std::swap (p1(i), p2(i));
-          std::reverse (step_sizes[i].begin(), step_sizes[i].end());
-        }
 
-      double x = 0;
-      for (size_type j=0; j<step_sizes.at(i).size(); j++)
-        x += step_sizes[i][j];
-      Assert(std::fabs(x - (p2(i)-p1(i))) <= 1e-12*std::fabs(x),
-             ExcInvalidRepetitions (i) );
-    }
 
+  template <int dim>
+  void
+  subdivided_hyper_rectangle(
+    Triangulation<dim>              &tria,
+    const std::vector<std::vector<double> > &step_sz,
+    const Point<dim>                &p_1,
+    const Point<dim>                &p_2,
+    const bool                       colorize)
+  {
+    // contributed by Joerg R. Weimar
+    // (j.weimar@jweimar.de) 2003
+    // modified by Yaqi Wang 2006
+    Assert(step_sz.size() == dim,
+           ExcInvalidRepetitionsDimension(dim));
+
+
+    // First, normalize input such that
+    // p1 is lower in all coordinate
+    // directions.
+
+    // and check the consistency of
+    // step sizes, i.e. that they all
+    // add up to the sizes specified by
+    // p_1 and p_2
+    Point<dim> p1(p_1);
+    Point<dim> p2(p_2);
+    std::vector< std::vector<double> > step_sizes(step_sz);
+
+    for (unsigned int i=0; i<dim; ++i)
+      {
+        if (p1(i) > p2(i))
+          {
+            std::swap (p1(i), p2(i));
+            std::reverse (step_sizes[i].begin(), step_sizes[i].end());
+          }
 
-  // then generate the necessary
-  // points
-  std::vector<Point<dim> > points;
-  switch (dim)
-    {
-    case 1:
-    {
-      double x=0;
-      for (size_type i=0; ; ++i)
-        {
-          points.push_back (Point<dim> (p1[0]+x));
-
-          // form partial sums. in
-          // the last run through
-          // avoid accessing
-          // non-existent values
-          // and exit early instead
-          if (i == step_sizes[0].size())
-            break;
+        double x = 0;
+        for (unsigned int j=0; j<step_sizes.at(i).size(); j++)
+          x += step_sizes[i][j];
+        Assert(std::fabs(x - (p2(i)-p1(i))) <= 1e-12*std::fabs(x),
+               ExcInvalidRepetitions (i) );
+      }
 
-          x += step_sizes[0][i];
-        }
-      break;
-    }
 
-    case 2:
-    {
-      double y=0;
-      for (size_type j=0; ; ++j)
-        {
-          double x=0;
-          for (size_type i=0; ; ++i)
-            {
-              points.push_back (Point<dim> (p1[0]+x,
-                                            p1[1]+y));
-              if (i == step_sizes[0].size())
-                break;
+    // then generate the necessary
+    // points
+    std::vector<Point<dim> > points;
+    switch (dim)
+      {
+      case 1:
+      {
+        double x=0;
+        for (unsigned int i=0; ; ++i)
+          {
+            points.push_back (Point<dim> (p1[0]+x));
 
-              x += step_sizes[0][i];
-            }
+            // form partial sums. in
+            // the last run through
+            // avoid accessing
+            // non-existent values
+            // and exit early instead
+            if (i == step_sizes[0].size())
+              break;
 
-          if (j == step_sizes[1].size())
-            break;
+            x += step_sizes[0][i];
+          }
+        break;
+      }
 
-          y += step_sizes[1][j];
-        }
-      break;
+      case 2:
+      {
+        double y=0;
+        for (unsigned int j=0; ; ++j)
+          {
+            double x=0;
+            for (unsigned int i=0; ; ++i)
+              {
+                points.push_back (Point<dim> (p1[0]+x,
+                                              p1[1]+y));
+                if (i == step_sizes[0].size())
+                  break;
 
-    }
-    case 3:
-    {
-      double z=0;
-      for (size_type k=0; ; ++k)
-        {
-          double y=0;
-          for (size_type j=0; ; ++j)
-            {
-              double x=0;
-              for (size_type i=0; ; ++i)
-                {
-                  points.push_back (Point<dim> (p1[0]+x,
-                                                p1[1]+y,
-                                                p1[2]+z));
-                  if (i == step_sizes[0].size())
-                    break;
+                x += step_sizes[0][i];
+              }
 
-                  x += step_sizes[0][i];
-                }
+            if (j == step_sizes[1].size())
+              break;
 
-              if (j == step_sizes[1].size())
-                break;
+            y += step_sizes[1][j];
+          }
+        break;
 
-              y += step_sizes[1][j];
-            }
+      }
+      case 3:
+      {
+        double z=0;
+        for (unsigned int k=0; ; ++k)
+          {
+            double y=0;
+            for (unsigned int j=0; ; ++j)
+              {
+                double x=0;
+                for (unsigned int i=0; ; ++i)
+                  {
+                    points.push_back (Point<dim> (p1[0]+x,
+                                                  p1[1]+y,
+                                                  p1[2]+z));
+                    if (i == step_sizes[0].size())
+                      break;
 
-          if (k == step_sizes[2].size())
-            break;
+                    x += step_sizes[0][i];
+                  }
 
-          z += step_sizes[2][k];
-        }
-      break;
-    }
+                if (j == step_sizes[1].size())
+                  break;
 
-    default:
-      Assert (false, ExcNotImplemented());
-    }
+                y += step_sizes[1][j];
+              }
 
-  // next create the cells
-  // Prepare cell data
-  std::vector<CellData<dim> > cells;
-  switch (dim)
-    {
-    case 1:
-    {
-      cells.resize (step_sizes[0].size());
-      for (size_type x=0; x<step_sizes[0].size(); ++x)
-        {
-          cells[x].vertices[0] = x;
-          cells[x].vertices[1] = x+1;
-          cells[x].material_id = 0;
-        }
-      break;
-    }
+            if (k == step_sizes[2].size())
+              break;
 
-    case 2:
-    {
-      cells.resize (step_sizes[1].size()*step_sizes[0].size());
-      for (size_type y=0; y<step_sizes[1].size(); ++y)
-        for (size_type x=0; x<step_sizes[0].size(); ++x)
-          {
-            const size_type c = x+y*step_sizes[0].size();
-            cells[c].vertices[0] = y*(step_sizes[0].size()+1)+x;
-            cells[c].vertices[1] = y*(step_sizes[0].size()+1)+x+1;
-            cells[c].vertices[2] = (y+1)*(step_sizes[0].size()+1)+x;
-            cells[c].vertices[3] = (y+1)*(step_sizes[0].size()+1)+x+1;
-            cells[c].material_id = 0;
+            z += step_sizes[2][k];
           }
-      break;
-    }
+        break;
+      }
 
-    case 3:
-    {
-      const size_type n_x  = (step_sizes[0].size()+1);
-      const size_type n_xy = (step_sizes[0].size()+1)*(step_sizes[1].size()+1);
+      default:
+        Assert (false, ExcNotImplemented());
+      }
+
+    // next create the cells
+    // Prepare cell data
+    std::vector<CellData<dim> > cells;
+    switch (dim)
+      {
+      case 1:
+      {
+        cells.resize (step_sizes[0].size());
+        for (unsigned int x=0; x<step_sizes[0].size(); ++x)
+          {
+            cells[x].vertices[0] = x;
+            cells[x].vertices[1] = x+1;
+            cells[x].material_id = 0;
+          }
+        break;
+      }
 
-      cells.resize (step_sizes[2].size()*step_sizes[1].size()*step_sizes[0].size());
-      for (size_type z=0; z<step_sizes[2].size(); ++z)
-        for (size_type y=0; y<step_sizes[1].size(); ++y)
-          for (size_type x=0; x<step_sizes[0].size(); ++x)
+      case 2:
+      {
+        cells.resize (step_sizes[1].size()*step_sizes[0].size());
+        for (unsigned int y=0; y<step_sizes[1].size(); ++y)
+          for (unsigned int x=0; x<step_sizes[0].size(); ++x)
             {
-              const size_type    c = x+y*step_sizes[0].size() +
-                                     z*step_sizes[0].size()*step_sizes[1].size();
-              cells[c].vertices[0] = z*n_xy + y*n_x + x;
-              cells[c].vertices[1] = z*n_xy + y*n_x + x+1;
-              cells[c].vertices[2] = z*n_xy + (y+1)*n_x + x;
-              cells[c].vertices[3] = z*n_xy + (y+1)*n_x + x+1;
-              cells[c].vertices[4] = (z+1)*n_xy + y*n_x + x;
-              cells[c].vertices[5] = (z+1)*n_xy + y*n_x + x+1;
-              cells[c].vertices[6] = (z+1)*n_xy + (y+1)*n_x + x;
-              cells[c].vertices[7] = (z+1)*n_xy + (y+1)*n_x + x+1;
+              const unsigned int c = x+y*step_sizes[0].size();
+              cells[c].vertices[0] = y*(step_sizes[0].size()+1)+x;
+              cells[c].vertices[1] = y*(step_sizes[0].size()+1)+x+1;
+              cells[c].vertices[2] = (y+1)*(step_sizes[0].size()+1)+x;
+              cells[c].vertices[3] = (y+1)*(step_sizes[0].size()+1)+x+1;
               cells[c].material_id = 0;
             }
-      break;
+        break;
+      }
 
-    }
+      case 3:
+      {
+        const unsigned int n_x  = (step_sizes[0].size()+1);
+        const unsigned int n_xy = (step_sizes[0].size()+1)*(step_sizes[1].size()+1);
 
-    default:
-      Assert (false, ExcNotImplemented());
-    }
+        cells.resize (step_sizes[2].size()*step_sizes[1].size()*step_sizes[0].size());
+        for (unsigned int z=0; z<step_sizes[2].size(); ++z)
+          for (unsigned int y=0; y<step_sizes[1].size(); ++y)
+            for (unsigned int x=0; x<step_sizes[0].size(); ++x)
+              {
+                const unsigned int    c = x+y*step_sizes[0].size() +
+                                          z*step_sizes[0].size()*step_sizes[1].size();
+                cells[c].vertices[0] = z*n_xy + y*n_x + x;
+                cells[c].vertices[1] = z*n_xy + y*n_x + x+1;
+                cells[c].vertices[2] = z*n_xy + (y+1)*n_x + x;
+                cells[c].vertices[3] = z*n_xy + (y+1)*n_x + x+1;
+                cells[c].vertices[4] = (z+1)*n_xy + y*n_x + x;
+                cells[c].vertices[5] = (z+1)*n_xy + y*n_x + x+1;
+                cells[c].vertices[6] = (z+1)*n_xy + (y+1)*n_x + x;
+                cells[c].vertices[7] = (z+1)*n_xy + (y+1)*n_x + x+1;
+                cells[c].material_id = 0;
+              }
+        break;
 
-  tria.create_triangulation (points, cells, SubCellData());
+      }
 
-  if (colorize)
-    {
-      // to colorize, run through all
-      // faces of all cells and set
-      // boundary indicator to the
-      // correct value if it was 0.
-
-      // use a large epsilon to
-      // compare numbers to avoid
-      // roundoff problems.
-      double min_size = *std::min_element (step_sizes[0].begin(),
-                                           step_sizes[0].end());
-      for (unsigned int i=1; i<dim; ++i)
-        min_size = std::min (min_size,
-                             *std::min_element (step_sizes[i].begin(),
-                                                step_sizes[i].end()));
-      const double epsilon = 0.01 * min_size;
-
-      // actual code is external since
-      // 1-D is different from 2/3D.
-      colorize_subdivided_hyper_rectangle (tria, p1, p2, epsilon);
-    }
-}
+      default:
+        Assert (false, ExcNotImplemented());
+      }
 
+    tria.create_triangulation (points, cells, SubCellData());
 
+    if (colorize)
+      {
+        // to colorize, run through all
+        // faces of all cells and set
+        // boundary indicator to the
+        // correct value if it was 0.
+
+        // use a large epsilon to
+        // compare numbers to avoid
+        // roundoff problems.
+        double min_size = *std::min_element (step_sizes[0].begin(),
+                                             step_sizes[0].end());
+        for (unsigned int i=1; i<dim; ++i)
+          min_size = std::min (min_size,
+                               *std::min_element (step_sizes[i].begin(),
+                                                  step_sizes[i].end()));
+        const double epsilon = 0.01 * min_size;
+
+        // actual code is external since
+        // 1-D is different from 2/3D.
+        colorize_subdivided_hyper_rectangle (tria, p1, p2, epsilon);
+      }
+  }
 
-template <>
-void
-GridGenerator::subdivided_hyper_rectangle (
-  Triangulation<1>                             &tria,
-  const std::vector< std::vector<double> >     &spacing,
-  const Point<1>                               &p,
-  const Table<1,types::material_id>                 &material_id,
-  const bool                                    colorize)
-{
-  // contributed by Yaqi Wang 2006
-  Assert(spacing.size() == 1,
-         ExcInvalidRepetitionsDimension(1));
 
-  const size_type n_cells = material_id.size(0);
 
-  Assert(spacing[0].size() == n_cells,
-         ExcInvalidRepetitionsDimension(1));
+  template <>
+  void
+  subdivided_hyper_rectangle (
+    Triangulation<1>                             &tria,
+    const std::vector< std::vector<double> >     &spacing,
+    const Point<1>                               &p,
+    const Table<1,types::material_id>                 &material_id,
+    const bool                                    colorize)
+  {
+    // contributed by Yaqi Wang 2006
+    Assert(spacing.size() == 1,
+           ExcInvalidRepetitionsDimension(1));
 
-  double delta = std::numeric_limits<double>::max();
-  for (size_type i=0; i<n_cells; i++)
-    {
-      Assert (spacing[0][i] >= 0, ExcInvalidRepetitions(-1));
-      delta = std::min (delta, spacing[0][i]);
-    }
+    const unsigned int n_cells = material_id.size(0);
 
-  // generate the necessary points
-  std::vector<Point<1> > points;
-  double ax = p[0];
-  for (size_type x=0; x<=n_cells; ++x)
-    {
-      points.push_back (Point<1> (ax));
-      if (x<n_cells)
-        ax += spacing[0][x];
-    }
-  // create the cells
-  size_type n_val_cells = 0;
-  for (size_type i=0; i<n_cells; i++)
-    if (material_id[i]!=numbers::invalid_material_id) n_val_cells++;
+    Assert(spacing[0].size() == n_cells,
+           ExcInvalidRepetitionsDimension(1));
 
-  std::vector<CellData<1> > cells(n_val_cells);
-  size_type id = 0;
-  for (size_type x=0; x<n_cells; ++x)
-    if (material_id[x] != numbers::invalid_material_id)
+    double delta = std::numeric_limits<double>::max();
+    for (unsigned int i=0; i<n_cells; i++)
       {
-        cells[id].vertices[0] = x;
-        cells[id].vertices[1] = x+1;
-        cells[id].material_id = material_id[x];
-        id++;
+        Assert (spacing[0][i] >= 0, ExcInvalidRepetitions(-1));
+        delta = std::min (delta, spacing[0][i]);
       }
-  // create triangulation
-  SubCellData t;
-  GridTools::delete_unused_vertices (points, cells, t);
 
-  tria.create_triangulation (points, cells, t);
+    // generate the necessary points
+    std::vector<Point<1> > points;
+    double ax = p[0];
+    for (unsigned int x=0; x<=n_cells; ++x)
+      {
+        points.push_back (Point<1> (ax));
+        if (x<n_cells)
+          ax += spacing[0][x];
+      }
+    // create the cells
+    unsigned int n_val_cells = 0;
+    for (unsigned int i=0; i<n_cells; i++)
+      if (material_id[i]!=numbers::invalid_material_id) n_val_cells++;
+
+    std::vector<CellData<1> > cells(n_val_cells);
+    unsigned int id = 0;
+    for (unsigned int x=0; x<n_cells; ++x)
+      if (material_id[x] != numbers::invalid_material_id)
+        {
+          cells[id].vertices[0] = x;
+          cells[id].vertices[1] = x+1;
+          cells[id].material_id = material_id[x];
+          id++;
+        }
+    // create triangulation
+    SubCellData t;
+    GridTools::delete_unused_vertices (points, cells, t);
 
-  // set boundary indicator
-  if (colorize)
-    Assert (false, ExcNotImplemented());
-}
+    tria.create_triangulation (points, cells, t);
 
+    // set boundary indicator
+    if (colorize)
+      Assert (false, ExcNotImplemented());
+  }
 
-template <>
-void
-GridGenerator::subdivided_hyper_rectangle (
-  Triangulation<2>                         &tria,
-  const std::vector< std::vector<double> >     &spacing,
-  const Point<2>                               &p,
-  const Table<2,types::material_id>          &material_id,
-  const bool                                    colorize)
-{
-  // contributed by Yaqi Wang 2006
-  Assert(spacing.size() == 2,
-         ExcInvalidRepetitionsDimension(2));
-
-  std::vector<size_type> repetitions(2);
-  size_type n_cells = 1;
-  double delta = std::numeric_limits<double>::max();
-  for (unsigned int i=0; i<2; i++)
-    {
-      repetitions[i] = spacing[i].size();
-      n_cells *= repetitions[i];
-      for (size_type j=0; j<repetitions[i]; j++)
-        {
-          Assert (spacing[i][j] >= 0, ExcInvalidRepetitions(-1));
-          delta = std::min (delta, spacing[i][j]);
-        }
-      Assert(material_id.size(i) == repetitions[i],
-             ExcInvalidRepetitionsDimension(i));
-    }
 
-  // generate the necessary points
-  std::vector<Point<2> > points;
-  double ay = p[1];
-  for (size_type y=0; y<=repetitions[1]; ++y)
-    {
-      double ax = p[0];
-      for (size_type x=0; x<=repetitions[0]; ++x)
-        {
-          points.push_back (Point<2> (ax,ay));
-          if (x<repetitions[0])
-            ax += spacing[0][x];
-        }
-      if (y<repetitions[1])
-        ay += spacing[1][y];
-    }
+  template <>
+  void
+  subdivided_hyper_rectangle (
+    Triangulation<2>                         &tria,
+    const std::vector< std::vector<double> >     &spacing,
+    const Point<2>                               &p,
+    const Table<2,types::material_id>          &material_id,
+    const bool                                    colorize)
+  {
+    // contributed by Yaqi Wang 2006
+    Assert(spacing.size() == 2,
+           ExcInvalidRepetitionsDimension(2));
+
+    std::vector<unsigned int> repetitions(2);
+    unsigned int n_cells = 1;
+    double delta = std::numeric_limits<double>::max();
+    for (unsigned int i=0; i<2; i++)
+      {
+        repetitions[i] = spacing[i].size();
+        n_cells *= repetitions[i];
+        for (unsigned int j=0; j<repetitions[i]; j++)
+          {
+            Assert (spacing[i][j] >= 0, ExcInvalidRepetitions(-1));
+            delta = std::min (delta, spacing[i][j]);
+          }
+        Assert(material_id.size(i) == repetitions[i],
+               ExcInvalidRepetitionsDimension(i));
+      }
 
-  // create the cells
-  size_type n_val_cells = 0;
-  for (size_type i=0; i<material_id.size(0); i++)
-    for (size_type j=0; j<material_id.size(1); j++)
-      if (material_id[i][j] != numbers::invalid_material_id)
-        n_val_cells++;
-
-  std::vector<CellData<2> > cells(n_val_cells);
-  size_type id = 0;
-  for (size_type y=0; y<repetitions[1]; ++y)
-    for (size_type x=0; x<repetitions[0]; ++x)
-      if (material_id[x][y]!=numbers::invalid_material_id)
-        {
-          cells[id].vertices[0] = y*(repetitions[0]+1)+x;
-          cells[id].vertices[1] = y*(repetitions[0]+1)+x+1;
-          cells[id].vertices[2] = (y+1)*(repetitions[0]+1)+x;
-          cells[id].vertices[3] = (y+1)*(repetitions[0]+1)+x+1;
-          cells[id].material_id = material_id[x][y];
-          id++;
-        }
+    // generate the necessary points
+    std::vector<Point<2> > points;
+    double ay = p[1];
+    for (unsigned int y=0; y<=repetitions[1]; ++y)
+      {
+        double ax = p[0];
+        for (unsigned int x=0; x<=repetitions[0]; ++x)
+          {
+            points.push_back (Point<2> (ax,ay));
+            if (x<repetitions[0])
+              ax += spacing[0][x];
+          }
+        if (y<repetitions[1])
+          ay += spacing[1][y];
+      }
 
-  // create triangulation
-  SubCellData t;
-  GridTools::delete_unused_vertices (points, cells, t);
+    // create the cells
+    unsigned int n_val_cells = 0;
+    for (unsigned int i=0; i<material_id.size(0); i++)
+      for (unsigned int j=0; j<material_id.size(1); j++)
+        if (material_id[i][j] != numbers::invalid_material_id)
+          n_val_cells++;
 
-  tria.create_triangulation (points, cells, t);
+    std::vector<CellData<2> > cells(n_val_cells);
+    unsigned int id = 0;
+    for (unsigned int y=0; y<repetitions[1]; ++y)
+      for (unsigned int x=0; x<repetitions[0]; ++x)
+        if (material_id[x][y]!=numbers::invalid_material_id)
+          {
+            cells[id].vertices[0] = y*(repetitions[0]+1)+x;
+            cells[id].vertices[1] = y*(repetitions[0]+1)+x+1;
+            cells[id].vertices[2] = (y+1)*(repetitions[0]+1)+x;
+            cells[id].vertices[3] = (y+1)*(repetitions[0]+1)+x+1;
+            cells[id].material_id = material_id[x][y];
+            id++;
+          }
 
-  // set boundary indicator
-  if (colorize)
-    {
-      double eps = 0.01 * delta;
-      Triangulation<2>::cell_iterator cell = tria.begin(),
-                                      endc = tria.end();
-      for (; cell !=endc; ++cell)
-        {
-          Point<2> cell_center = cell->center();
-          for (unsigned int f=0; f<GeometryInfo<2>::faces_per_cell; ++f)
-            if (cell->face(f)->boundary_indicator() == 0)
-              {
-                Point<2> face_center = cell->face(f)->center();
-                for (unsigned int i=0; i<2; ++i)
-                  {
-                    if (face_center[i]<cell_center[i]-eps)
-                      cell->face(f)->set_boundary_indicator(i*2);
-                    if (face_center[i]>cell_center[i]+eps)
-                      cell->face(f)->set_boundary_indicator(i*2+1);
-                  }
-              }
-        }
-    }
-}
+    // create triangulation
+    SubCellData t;
+    GridTools::delete_unused_vertices (points, cells, t);
 
+    tria.create_triangulation (points, cells, t);
 
-template <>
-void
-GridGenerator::subdivided_hyper_rectangle (
-  Triangulation<3>                           &tria,
-  const std::vector< std::vector<double> >     &spacing,
-  const Point<3>                             &p,
-  const Table<3,types::material_id>               &material_id,
-  const bool                                    colorize)
-{
-  // contributed by Yaqi Wang 2006
-  const unsigned int dim = 3;
+    // set boundary indicator
+    if (colorize)
+      {
+        double eps = 0.01 * delta;
+        Triangulation<2>::cell_iterator cell = tria.begin(),
+                                        endc = tria.end();
+        for (; cell !=endc; ++cell)
+          {
+            Point<2> cell_center = cell->center();
+            for (unsigned int f=0; f<GeometryInfo<2>::faces_per_cell; ++f)
+              if (cell->face(f)->boundary_indicator() == 0)
+                {
+                  Point<2> face_center = cell->face(f)->center();
+                  for (unsigned int i=0; i<2; ++i)
+                    {
+                      if (face_center[i]<cell_center[i]-eps)
+                        cell->face(f)->set_boundary_indicator(i*2);
+                      if (face_center[i]>cell_center[i]+eps)
+                        cell->face(f)->set_boundary_indicator(i*2+1);
+                    }
+                }
+          }
+      }
+  }
 
-  Assert(spacing.size() == dim,
-         ExcInvalidRepetitionsDimension(dim));
 
-  std::vector<size_type > repetitions(dim);
-  size_type n_cells = 1;
-  double delta = std::numeric_limits<double>::max();
-  for (unsigned int i=0; i<dim; i++)
-    {
-      repetitions[i] = spacing[i].size();
-      n_cells *= repetitions[i];
-      for (size_type j=0; j<repetitions[i]; j++)
-        {
-          Assert (spacing[i][j] >= 0, ExcInvalidRepetitions(-1));
-          delta = std::min (delta, spacing[i][j]);
-        }
-      Assert(material_id.size(i) == repetitions[i],
-             ExcInvalidRepetitionsDimension(i));
-    }
+  template <>
+  void
+  subdivided_hyper_rectangle (
+    Triangulation<3>                           &tria,
+    const std::vector< std::vector<double> >     &spacing,
+    const Point<3>                             &p,
+    const Table<3,types::material_id>               &material_id,
+    const bool                                    colorize)
+  {
+    // contributed by Yaqi Wang 2006
+    const unsigned int dim = 3;
 
-  // generate the necessary points
-  std::vector<Point<dim> > points;
-  double az = p[2];
-  for (size_type z=0; z<=repetitions[2]; ++z)
-    {
-      double ay = p[1];
-      for (size_type y=0; y<=repetitions[1]; ++y)
-        {
-          double ax = p[0];
-          for (size_type x=0; x<=repetitions[0]; ++x)
+    Assert(spacing.size() == dim,
+           ExcInvalidRepetitionsDimension(dim));
+
+    std::vector<unsigned int > repetitions(dim);
+    unsigned int n_cells = 1;
+    double delta = std::numeric_limits<double>::max();
+    for (unsigned int i=0; i<dim; i++)
+      {
+        repetitions[i] = spacing[i].size();
+        n_cells *= repetitions[i];
+        for (unsigned int j=0; j<repetitions[i]; j++)
+          {
+            Assert (spacing[i][j] >= 0, ExcInvalidRepetitions(-1));
+            delta = std::min (delta, spacing[i][j]);
+          }
+        Assert(material_id.size(i) == repetitions[i],
+               ExcInvalidRepetitionsDimension(i));
+      }
+
+    // generate the necessary points
+    std::vector<Point<dim> > points;
+    double az = p[2];
+    for (unsigned int z=0; z<=repetitions[2]; ++z)
+      {
+        double ay = p[1];
+        for (unsigned int y=0; y<=repetitions[1]; ++y)
+          {
+            double ax = p[0];
+            for (unsigned int x=0; x<=repetitions[0]; ++x)
+              {
+                points.push_back (Point<dim> (ax,ay,az));
+                if (x<repetitions[0])
+                  ax += spacing[0][x];
+              }
+            if (y<repetitions[1])
+              ay += spacing[1][y];
+          }
+        if (z<repetitions[2])
+          az += spacing[2][z];
+      }
+
+    // create the cells
+    unsigned int n_val_cells = 0;
+    for (unsigned int i=0; i<material_id.size(0); i++)
+      for (unsigned int j=0; j<material_id.size(1); j++)
+        for (unsigned int k=0; k<material_id.size(2); k++)
+          if (material_id[i][j][k]!=numbers::invalid_material_id)
+            n_val_cells++;
+
+    std::vector<CellData<dim> > cells(n_val_cells);
+    unsigned int id = 0;
+    const unsigned int n_x  = (repetitions[0]+1);
+    const unsigned int n_xy = (repetitions[0]+1)*(repetitions[1]+1);
+    for (unsigned int z=0; z<repetitions[2]; ++z)
+      for (unsigned int y=0; y<repetitions[1]; ++y)
+        for (unsigned int x=0; x<repetitions[0]; ++x)
+          if (material_id[x][y][z]!=numbers::invalid_material_id)
             {
-              points.push_back (Point<dim> (ax,ay,az));
-              if (x<repetitions[0])
-                ax += spacing[0][x];
+              cells[id].vertices[0] = z*n_xy + y*n_x + x;
+              cells[id].vertices[1] = z*n_xy + y*n_x + x+1;
+              cells[id].vertices[2] = z*n_xy + (y+1)*n_x + x;
+              cells[id].vertices[3] = z*n_xy + (y+1)*n_x + x+1;
+              cells[id].vertices[4] = (z+1)*n_xy + y*n_x + x;
+              cells[id].vertices[5] = (z+1)*n_xy + y*n_x + x+1;
+              cells[id].vertices[6] = (z+1)*n_xy + (y+1)*n_x + x;
+              cells[id].vertices[7] = (z+1)*n_xy + (y+1)*n_x + x+1;
+              cells[id].material_id = material_id[x][y][z];
+              id++;
             }
-          if (y<repetitions[1])
-            ay += spacing[1][y];
-        }
-      if (z<repetitions[2])
-        az += spacing[2][z];
-    }
 
-  // create the cells
-  size_type n_val_cells = 0;
-  for (size_type i=0; i<material_id.size(0); i++)
-    for (size_type j=0; j<material_id.size(1); j++)
-      for (size_type k=0; k<material_id.size(2); k++)
-        if (material_id[i][j][k]!=numbers::invalid_material_id)
-          n_val_cells++;
+    // create triangulation
+    SubCellData t;
+    GridTools::delete_unused_vertices (points, cells, t);
 
-  std::vector<CellData<dim> > cells(n_val_cells);
-  size_type id = 0;
-  const size_type n_x  = (repetitions[0]+1);
-  const size_type n_xy = (repetitions[0]+1)*(repetitions[1]+1);
-  for (size_type z=0; z<repetitions[2]; ++z)
-    for (size_type y=0; y<repetitions[1]; ++y)
-      for (size_type x=0; x<repetitions[0]; ++x)
-        if (material_id[x][y][z]!=numbers::invalid_material_id)
+    tria.create_triangulation (points, cells, t);
+
+    // set boundary indicator
+    if (colorize && dim>1)
+      {
+        double eps = 0.01 * delta;
+        Triangulation<dim>::cell_iterator cell = tria.begin(),
+                                          endc = tria.end();
+        for (; cell !=endc; ++cell)
           {
-            cells[id].vertices[0] = z*n_xy + y*n_x + x;
-            cells[id].vertices[1] = z*n_xy + y*n_x + x+1;
-            cells[id].vertices[2] = z*n_xy + (y+1)*n_x + x;
-            cells[id].vertices[3] = z*n_xy + (y+1)*n_x + x+1;
-            cells[id].vertices[4] = (z+1)*n_xy + y*n_x + x;
-            cells[id].vertices[5] = (z+1)*n_xy + y*n_x + x+1;
-            cells[id].vertices[6] = (z+1)*n_xy + (y+1)*n_x + x;
-            cells[id].vertices[7] = (z+1)*n_xy + (y+1)*n_x + x+1;
-            cells[id].material_id = material_id[x][y][z];
-            id++;
+            Point<dim> cell_center = cell->center();
+            for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+              if (cell->face(f)->boundary_indicator() == 0)
+                {
+                  Point<dim> face_center = cell->face(f)->center();
+                  for (unsigned int i=0; i<dim; ++i)
+                    {
+                      if (face_center[i]<cell_center[i]-eps)
+                        cell->face(f)->set_boundary_indicator(i*2);
+                      if (face_center[i]>cell_center[i]+eps)
+                        cell->face(f)->set_boundary_indicator(i*2+1);
+                    }
+                }
           }
+      }
+  }
 
-  // create triangulation
-  SubCellData t;
-  GridTools::delete_unused_vertices (points, cells, t);
 
-  tria.create_triangulation (points, cells, t);
 
-  // set boundary indicator
-  if (colorize && dim>1)
-    {
-      double eps = 0.01 * delta;
-      Triangulation<dim>::cell_iterator cell = tria.begin(),
-                                        endc = tria.end();
-      for (; cell !=endc; ++cell)
-        {
-          Point<dim> cell_center = cell->center();
-          for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
-            if (cell->face(f)->boundary_indicator() == 0)
-              {
-                Point<dim> face_center = cell->face(f)->center();
-                for (unsigned int i=0; i<dim; ++i)
-                  {
-                    if (face_center[i]<cell_center[i]-eps)
-                      cell->face(f)->set_boundary_indicator(i*2);
-                    if (face_center[i]>cell_center[i]+eps)
-                      cell->face(f)->set_boundary_indicator(i*2+1);
-                  }
-              }
-        }
-    }
-}
+  template <>
+  void hyper_cube_slit (Triangulation<1> &,
+                        const double,
+                        const double,
+                        const bool)
+  {
+    Assert (false, ExcNotImplemented());
+  }
 
 
 
-template <>
-void
-GridGenerator::colorize_subdivided_hyper_rectangle (
-  Triangulation<1> &tria,
-  const Point<1> &,
-  const Point<1> &,
-  const double)
-{
-  for (Triangulation<1>::cell_iterator cell = tria.begin();
-       cell != tria.end(); ++cell)
-    if (cell->center()(0) > 0)
-      cell->set_material_id(1);
-  // boundary indicators are set to
-  // 0 (left) and 1 (right) by default.
-}
+  template <>
+  void enclosed_hyper_cube (Triangulation<1> &,
+                            const double,
+                            const double,
+                            const double,
+                            const bool)
+  {
+    Assert (false, ExcNotImplemented());
+  }
 
 
 
-template <int dim>
-void
-GridGenerator::colorize_subdivided_hyper_rectangle (Triangulation<dim> &tria,
-                                                    const Point<dim>   &p1,
-                                                    const Point<dim>   &p2,
-                                                    const double        epsilon)
-{
+  template <>
+  void hyper_L (Triangulation<1> &,
+                const double,
+                const double)
+  {
+    Assert (false, ExcNotImplemented());
+  }
 
-  // run through all faces and check
-  // if one of their center coordinates matches
-  // one of the corner points. Comparisons
-  // are made using an epsilon which
-  // should be smaller than the smallest cell
-  // diameter.
 
-  typename Triangulation<dim>::face_iterator face = tria.begin_face(),
-                                             endface = tria.end_face();
-  for (; face!=endface; ++face)
-    {
-      if (face->boundary_indicator() == 0)
-        {
-          const Point<dim> center (face->center());
-          if (std::abs(center(0)-p1[0]) < epsilon)
-            face->set_boundary_indicator(0);
-          else if (std::abs(center(0) - p2[0]) < epsilon)
-            face->set_boundary_indicator(1);
-          else if (dim > 1 && std::abs(center(1) - p1[1]) < epsilon)
-            face->set_boundary_indicator(2);
-          else if (dim > 1 && std::abs(center(1) - p2[1]) < epsilon)
-            face->set_boundary_indicator(3);
-          else if (dim > 2 && std::abs(center(2) - p1[2]) < epsilon)
-            face->set_boundary_indicator(4);
-          else if (dim > 2 && std::abs(center(2) - p2[2]) < epsilon)
-            face->set_boundary_indicator(5);
-          else
-            // triangulation says it
-            // is on the boundary,
-            // but we could not find
-            // on which boundary.
-            Assert (false, ExcInternalError());
 
-        }
-    }
-  for (typename Triangulation<dim>::cell_iterator cell = tria.begin();
-       cell != tria.end(); ++cell)
-    {
-      char id = 0;
-      for (unsigned int d=0; d<dim; ++d)
-        if (cell->center()(d) > 0) id += 1 << d;
-      cell->set_material_id(id);
-    }
-}
+  template <>
+  void hyper_ball (Triangulation<1> &,
+                   const Point<1> &,
+                   const double)
+  {
+    Assert (false, ExcNotImplemented());
+  }
 
 
-template <>
-void GridGenerator::hyper_cube_slit (Triangulation<1> &,
-                                     const double,
-                                     const double,
-                                     const bool)
-{
-  Assert (false, ExcNotImplemented());
-}
 
+  template <>
+  void cylinder (Triangulation<1> &,
+                 const double,
+                 const double)
+  {
+    Assert (false, ExcNotImplemented());
+  }
 
 
-template <>
-void GridGenerator::enclosed_hyper_cube (Triangulation<1> &,
-                                         const double,
-                                         const double,
-                                         const double,
-                                         const bool)
-{
-  Assert (false, ExcNotImplemented());
-}
 
+  template <>
+  void truncated_cone (Triangulation<1> &,
+                       const double,
+                       const double,
+                       const double)
+  {
+    Assert (false, ExcNotImplemented());
+  }
 
 
-template <>
-void GridGenerator::hyper_L (Triangulation<1> &,
-                             const double,
-                             const double)
-{
-  Assert (false, ExcNotImplemented());
-}
 
+  template <>
+  void hyper_shell (Triangulation<1> &,
+                    const Point<1> &,
+                    const double,
+                    const double,
+                    const unsigned int ,
+                    const bool)
+  {
+    Assert (false, ExcNotImplemented());
+  }
 
 
-template <>
-void GridGenerator::hyper_ball (Triangulation<1> &,
-                                const Point<1> &,
-                                const double)
-{
-  Assert (false, ExcNotImplemented());
-}
+  template <>
+  void cylinder_shell (Triangulation<1> &,
+                       const double,
+                       const double,
+                       const double,
+                       const unsigned int ,
+                       const unsigned int )
+  {
+    Assert (false, ExcNotImplemented());
+  }
 
 
+  template <>
+  void
+  half_hyper_ball (Triangulation<1> &,
+                   const Point<1> &,
+                   const double)
+  {
+    Assert (false, ExcNotImplemented());
+  }
 
-template <>
-void GridGenerator::cylinder (Triangulation<1> &,
-                              const double,
-                              const double)
-{
-  Assert (false, ExcNotImplemented());
-}
 
+  template <>
+  void
+  half_hyper_shell (Triangulation<1> &,
+                    const Point<1> &,
+                    const double,
+                    const double,
+                    const unsigned int ,
+                    const bool)
+  {
+    Assert (false, ExcNotImplemented());
+  }
+
+  template <>
+  void quarter_hyper_shell (Triangulation<1> &,
+                            const Point<1> &,
+                            const double,
+                            const double,
+                            const unsigned int ,
+                            const bool)
+  {
+    Assert (false, ExcNotImplemented());
+  }
+
+  template <>
+  void enclosed_hyper_cube (Triangulation<2> &tria,
+                            const double        left,
+                            const double        right,
+                            const double        thickness,
+                            const bool          colorize)
+  {
+    Assert(left<right,
+           ExcMessage ("Invalid left-to-right bounds of enclosed hypercube"));
+
+    std::vector<Point<2> > vertices(16);
+    double coords[4];
+    coords[0] = left-thickness;
+    coords[1] = left;
+    coords[2] = right;
+    coords[3] = right+thickness;
+
+    unsigned int k=0;
+    for (unsigned int i0=0; i0<4; ++i0)
+      for (unsigned int i1=0; i1<4; ++i1)
+        vertices[k++] = Point<2>(coords[i1], coords[i0]);
+
+    const types::material_id materials[9] = { 5, 4, 6,
+                                              1, 0, 2,
+                                              9, 8,10
+                                            };
+
+    std::vector<CellData<2> > cells(9);
+    k = 0;
+    for (unsigned int i0=0; i0<3; ++i0)
+      for (unsigned int i1=0; i1<3; ++i1)
+        {
+          cells[k].vertices[0] = i1+4*i0;
+          cells[k].vertices[1] = i1+4*i0+1;
+          cells[k].vertices[2] = i1+4*i0+4;
+          cells[k].vertices[3] = i1+4*i0+5;
+          if (colorize)
+            cells[k].material_id = materials[k];
+          ++k;
+        }
+    tria.create_triangulation (vertices,
+                               cells,
+                               SubCellData());       // no boundary information
+  }
 
 
-template <>
-void GridGenerator::truncated_cone (Triangulation<1> &,
-                                    const double,
-                                    const double,
-                                    const double)
-{
-  Assert (false, ExcNotImplemented());
-}
 
+// Implementation for 2D only
+  template <>
+  void
+  hyper_cube_slit (Triangulation<2> &tria,
+                   const double left,
+                   const double right,
+                   const bool colorize)
+  {
+    const double rl2=(right+left)/2;
+    const Point<2> vertices[10] = { Point<2>(left, left ),
+                                    Point<2>(rl2,  left ),
+                                    Point<2>(rl2,  rl2  ),
+                                    Point<2>(left, rl2  ),
+                                    Point<2>(right,left ),
+                                    Point<2>(right,rl2  ),
+                                    Point<2>(rl2,  right),
+                                    Point<2>(left, right),
+                                    Point<2>(right,right),
+                                    Point<2>(rl2,  left )
+                                  };
+    const int cell_vertices[4][4] = { { 0,1,3,2 },
+      { 9,4,2,5 },
+      { 3,2,7,6 },
+      { 2,5,6,8 }
+    };
+    std::vector<CellData<2> > cells (4, CellData<2>());
+    for (unsigned int i=0; i<4; ++i)
+      {
+        for (unsigned int j=0; j<4; ++j)
+          cells[i].vertices[j] = cell_vertices[i][j];
+        cells[i].material_id = 0;
+      };
+    tria.create_triangulation (
+      std::vector<Point<2> >(&vertices[0], &vertices[10]),
+      cells,
+      SubCellData());       // no boundary information
 
+    if (colorize)
+      {
+        Triangulation<2>::cell_iterator cell = tria.begin();
+        cell->face(1)->set_boundary_indicator(1);
+        ++cell;
+        cell->face(3)->set_boundary_indicator(2);
+      }
+  }
 
-template <>
-void GridGenerator::hyper_shell (Triangulation<1> &,
-                                 const Point<1> &,
-                                 const double,
-                                 const double,
-                                 const size_type ,
-                                 const bool)
-{
-  Assert (false, ExcNotImplemented());
-}
 
 
-template <>
-void GridGenerator::colorize_hyper_shell (Triangulation<1> &,
-                                          const Point<1> &,
-                                          const double,
-                                          const double)
-{
-  Assert (false, ExcNotImplemented());
-}
+  template <>
+  void truncated_cone (Triangulation<2> &triangulation,
+                       const double radius_0,
+                       const double radius_1,
+                       const double half_length)
+  {
+    Point<2> vertices_tmp[4];
 
+    vertices_tmp[0] = Point<2> (-half_length, -radius_0);
+    vertices_tmp[1] = Point<2> (half_length, -radius_1);
+    vertices_tmp[2] = Point<2> (-half_length, radius_0);
+    vertices_tmp[3] = Point<2> (half_length, radius_1);
 
-template <>
-void GridGenerator::cylinder_shell (Triangulation<1> &,
-                                    const double,
-                                    const double,
-                                    const double,
-                                    const size_type ,
-                                    const size_type )
-{
-  Assert (false, ExcNotImplemented());
-}
+    const std::vector<Point<2> > vertices (&vertices_tmp[0], &vertices_tmp[4]);
+    unsigned int cell_vertices[1][GeometryInfo<2>::vertices_per_cell];
 
+    for (unsigned int i = 0; i < GeometryInfo<2>::vertices_per_cell; ++i)
+      cell_vertices[0][i] = i;
 
-template <>
-void
-GridGenerator::half_hyper_ball (Triangulation<1> &,
-                                const Point<1> &,
-                                const double)
-{
-  Assert (false, ExcNotImplemented());
-}
+    std::vector<CellData<2> > cells (1, CellData<2> ());
 
+    for (unsigned int i = 0; i < GeometryInfo<2>::vertices_per_cell; ++i)
+      cells[0].vertices[i] = cell_vertices[0][i];
 
-template <>
-void
-GridGenerator::half_hyper_shell (Triangulation<1> &,
-                                 const Point<1> &,
-                                 const double,
-                                 const double,
-                                 const size_type ,
-                                 const bool)
-{
-  Assert (false, ExcNotImplemented());
-}
+    cells[0].material_id = 0;
+    triangulation.create_triangulation (vertices, cells, SubCellData ());
 
-template <>
-void GridGenerator::quarter_hyper_shell (Triangulation<1> &,
-                                         const Point<1> &,
-                                         const double,
-                                         const double,
-                                         const size_type ,
-                                         const bool)
-{
-  Assert (false, ExcNotImplemented());
-}
+    Triangulation<2>::cell_iterator cell = triangulation.begin ();
 
-template <>
-void GridGenerator::enclosed_hyper_cube (Triangulation<2> &tria,
-                                         const double        left,
-                                         const double        right,
-                                         const double        thickness,
-                                         const bool          colorize)
-{
-  Assert(left<right,
-         ExcMessage ("Invalid left-to-right bounds of enclosed hypercube"));
-
-  std::vector<Point<2> > vertices(16);
-  double coords[4];
-  coords[0] = left-thickness;
-  coords[1] = left;
-  coords[2] = right;
-  coords[3] = right+thickness;
-
-  unsigned int k=0;
-  for (unsigned int i0=0; i0<4; ++i0)
-    for (unsigned int i1=0; i1<4; ++i1)
-      vertices[k++] = Point<2>(coords[i1], coords[i0]);
-
-  const types::material_id materials[9] = { 5, 4, 6,
-                                            1, 0, 2,
-                                            9, 8,10
-                                          };
-
-  std::vector<CellData<2> > cells(9);
-  k = 0;
-  for (unsigned int i0=0; i0<3; ++i0)
-    for (unsigned int i1=0; i1<3; ++i1)
-      {
-        cells[k].vertices[0] = i1+4*i0;
-        cells[k].vertices[1] = i1+4*i0+1;
-        cells[k].vertices[2] = i1+4*i0+4;
-        cells[k].vertices[3] = i1+4*i0+5;
-        if (colorize)
-          cells[k].material_id = materials[k];
-        ++k;
-      }
-  tria.create_triangulation (vertices,
-                             cells,
-                             SubCellData());       // no boundary information
-}
+    cell->face (0)->set_boundary_indicator (1);
+    cell->face (1)->set_boundary_indicator (2);
+
+    for (unsigned int i = 2; i < 4; ++i)
+      cell->face (i)->set_boundary_indicator (0);
+  }
 
 
 
+//TODO: Colorize edges as circumference, left and right radius
 // Implementation for 2D only
-template <>
-void
-GridGenerator::hyper_cube_slit (Triangulation<2> &tria,
-                                const double left,
-                                const double right,
-                                const bool colorize)
-{
-  const double rl2=(right+left)/2;
-  const Point<2> vertices[10] = { Point<2>(left, left ),
-                                  Point<2>(rl2,  left ),
-                                  Point<2>(rl2,  rl2  ),
-                                  Point<2>(left, rl2  ),
-                                  Point<2>(right,left ),
-                                  Point<2>(right,rl2  ),
-                                  Point<2>(rl2,  right),
-                                  Point<2>(left, right),
-                                  Point<2>(right,right),
-                                  Point<2>(rl2,  left )
-                                };
-  const int cell_vertices[4][4] = { { 0,1,3,2 },
-    { 9,4,2,5 },
-    { 3,2,7,6 },
-    { 2,5,6,8 }
-  };
-  std::vector<CellData<2> > cells (4, CellData<2>());
-  for (unsigned int i=0; i<4; ++i)
-    {
-      for (unsigned int j=0; j<4; ++j)
-        cells[i].vertices[j] = cell_vertices[i][j];
-      cells[i].material_id = 0;
+  template <>
+  void
+  hyper_L (Triangulation<2> &tria,
+           const double a,
+           const double b)
+  {
+    const Point<2> vertices[8] = { Point<2> (a,a),
+                                   Point<2> ((a+b)/2,a),
+                                   Point<2> (b,a),
+                                   Point<2> (a,(a+b)/2),
+                                   Point<2> ((a+b)/2,(a+b)/2),
+                                   Point<2> (b,(a+b)/2),
+                                   Point<2> (a,b),
+                                   Point<2> ((a+b)/2,b)
+                                 };
+    const int cell_vertices[3][4] = {{0, 1, 3, 4},
+      {1, 2, 4, 5},
+      {3, 4, 6, 7}
     };
-  tria.create_triangulation (
-    std::vector<Point<2> >(&vertices[0], &vertices[10]),
-    cells,
-    SubCellData());       // no boundary information
 
-  if (colorize)
-    {
-      Triangulation<2>::cell_iterator cell = tria.begin();
-      cell->face(1)->set_boundary_indicator(1);
-      ++cell;
-      cell->face(3)->set_boundary_indicator(2);
-    }
-}
+    std::vector<CellData<2> > cells (3, CellData<2>());
 
+    for (unsigned int i=0; i<3; ++i)
+      {
+        for (unsigned int j=0; j<4; ++j)
+          cells[i].vertices[j] = cell_vertices[i][j];
+        cells[i].material_id = 0;
+      };
 
+    tria.create_triangulation (
+      std::vector<Point<2> >(&vertices[0], &vertices[8]),
+      cells,
+      SubCellData());       // no boundary information
+  }
 
-template <>
-void GridGenerator::truncated_cone (Triangulation<2> &triangulation,
-                                    const double radius_0,
-                                    const double radius_1,
-                                    const double half_length)
-{
-  Point<2> vertices_tmp[4];
 
-  vertices_tmp[0] = Point<2> (-half_length, -radius_0);
-  vertices_tmp[1] = Point<2> (half_length, -radius_1);
-  vertices_tmp[2] = Point<2> (-half_length, radius_0);
-  vertices_tmp[3] = Point<2> (half_length, radius_1);
 
-  const std::vector<Point<2> > vertices (&vertices_tmp[0], &vertices_tmp[4]);
-  unsigned int cell_vertices[1][GeometryInfo<2>::vertices_per_cell];
+// Implementation for 2D only
+  template <>
+  void
+  hyper_ball (Triangulation<2> &tria,
+              const Point<2>   &p,
+              const double      radius)
+  {
+    // equilibrate cell sizes at
+    // transition from the inner part
+    // to the radial cells
+    const double a = 1./(1+std::sqrt(2.0));
+    const Point<2> vertices[8] = { p+Point<2>(-1,-1) *(radius/std::sqrt(2.0)),
+                                   p+Point<2>(+1,-1) *(radius/std::sqrt(2.0)),
+                                   p+Point<2>(-1,-1) *(radius/std::sqrt(2.0)*a),
+                                   p+Point<2>(+1,-1) *(radius/std::sqrt(2.0)*a),
+                                   p+Point<2>(-1,+1) *(radius/std::sqrt(2.0)*a),
+                                   p+Point<2>(+1,+1) *(radius/std::sqrt(2.0)*a),
+                                   p+Point<2>(-1,+1) *(radius/std::sqrt(2.0)),
+                                   p+Point<2>(+1,+1) *(radius/std::sqrt(2.0))
+                                 };
+
+    const int cell_vertices[5][4] = {{0, 1, 2, 3},
+      {0, 2, 6, 4},
+      {2, 3, 4, 5},
+      {1, 7, 3, 5},
+      {6, 4, 7, 5}
+    };
 
-  for (unsigned int i = 0; i < GeometryInfo<2>::vertices_per_cell; ++i)
-    cell_vertices[0][i] = i;
+    std::vector<CellData<2> > cells (5, CellData<2>());
 
-  std::vector<CellData<2> > cells (1, CellData<2> ());
+    for (unsigned int i=0; i<5; ++i)
+      {
+        for (unsigned int j=0; j<4; ++j)
+          cells[i].vertices[j] = cell_vertices[i][j];
+        cells[i].material_id = 0;
+      };
 
-  for (unsigned int i = 0; i < GeometryInfo<2>::vertices_per_cell; ++i)
-    cells[0].vertices[i] = cell_vertices[0][i];
+    tria.create_triangulation (
+      std::vector<Point<2> >(&vertices[0], &vertices[8]),
+      cells,
+      SubCellData());       // no boundary information
+  }
 
-  cells[0].material_id = 0;
-  triangulation.create_triangulation (vertices, cells, SubCellData ());
 
-  Triangulation<2>::cell_iterator cell = triangulation.begin ();
 
-  cell->face (0)->set_boundary_indicator (1);
-  cell->face (1)->set_boundary_indicator (2);
+  template <>
+  void hyper_shell (Triangulation<2> &tria,
+                    const Point<2>   &center,
+                    const double      inner_radius,
+                    const double      outer_radius,
+                    const unsigned int   n_cells,
+                    const bool colorize)
+  {
+    Assert ((inner_radius > 0) && (inner_radius < outer_radius),
+            ExcInvalidRadii ());
+
+    const double pi = numbers::PI;
+
+    // determine the number of cells
+    // for the grid. if not provided by
+    // the user determine it such that
+    // the length of each cell on the
+    // median (in the middle between
+    // the two circles) is equal to its
+    // radial extent (which is the
+    // difference between the two
+    // radii)
+    const unsigned int N = (n_cells == 0 ?
+                            static_cast<unsigned int>
+                            (std::ceil((2*pi* (outer_radius + inner_radius)/2) /
+                                       (outer_radius - inner_radius))) :
+                            n_cells);
+
+    // set up N vertices on the
+    // outer and N vertices on
+    // the inner circle. the
+    // first N ones are on the
+    // outer one, and all are
+    // numbered counter-clockwise
+    std::vector<Point<2> > vertices(2*N);
+    for (unsigned int i=0; i<N; ++i)
+      {
+        vertices[i] = Point<2>( std::cos(2*pi * i/N),
+                                std::sin(2*pi * i/N)) * outer_radius;
+        vertices[i+N] = vertices[i] * (inner_radius/outer_radius);
 
-  for (unsigned int i = 2; i < 4; ++i)
-    cell->face (i)->set_boundary_indicator (0);
-}
+        vertices[i]   += center;
+        vertices[i+N] += center;
+      };
 
+    std::vector<CellData<2> > cells (N, CellData<2>());
 
+    for (unsigned int i=0; i<N; ++i)
+      {
+        cells[i].vertices[0] = i;
+        cells[i].vertices[1] = (i+1)%N;
+        cells[i].vertices[2] = N+i;
+        cells[i].vertices[3] = N+((i+1)%N);
 
-//TODO: Colorize edges as circumference, left and right radius
-// Implementation for 2D only
-template <>
-void
-GridGenerator::hyper_L (Triangulation<2> &tria,
-                        const double a,
-                        const double b)
-{
-  const Point<2> vertices[8] = { Point<2> (a,a),
-                                 Point<2> ((a+b)/2,a),
-                                 Point<2> (b,a),
-                                 Point<2> (a,(a+b)/2),
-                                 Point<2> ((a+b)/2,(a+b)/2),
-                                 Point<2> (b,(a+b)/2),
-                                 Point<2> (a,b),
-                                 Point<2> ((a+b)/2,b)
-                               };
-  const int cell_vertices[3][4] = {{0, 1, 3, 4},
-    {1, 2, 4, 5},
-    {3, 4, 6, 7}
-  };
-
-  std::vector<CellData<2> > cells (3, CellData<2>());
-
-  for (unsigned int i=0; i<3; ++i)
-    {
-      for (unsigned int j=0; j<4; ++j)
-        cells[i].vertices[j] = cell_vertices[i][j];
-      cells[i].material_id = 0;
-    };
+        cells[i].material_id = 0;
+      };
 
-  tria.create_triangulation (
-    std::vector<Point<2> >(&vertices[0], &vertices[8]),
-    cells,
-    SubCellData());       // no boundary information
-}
+    tria.create_triangulation (
+      vertices, cells, SubCellData());
 
+    if (colorize)
+      colorize_hyper_shell(tria, center, inner_radius, outer_radius);
+  }
 
 
 // Implementation for 2D only
-template <>
-void
-GridGenerator::hyper_ball (Triangulation<2> &tria,
-                           const Point<2>   &p,
-                           const double      radius)
-{
-  // equilibrate cell sizes at
-  // transition from the inner part
-  // to the radial cells
-  const double a = 1./(1+std::sqrt(2.0));
-  const Point<2> vertices[8] = { p+Point<2>(-1,-1) *(radius/std::sqrt(2.0)),
-                                 p+Point<2>(+1,-1) *(radius/std::sqrt(2.0)),
-                                 p+Point<2>(-1,-1) *(radius/std::sqrt(2.0)*a),
-                                 p+Point<2>(+1,-1) *(radius/std::sqrt(2.0)*a),
-                                 p+Point<2>(-1,+1) *(radius/std::sqrt(2.0)*a),
-                                 p+Point<2>(+1,+1) *(radius/std::sqrt(2.0)*a),
-                                 p+Point<2>(-1,+1) *(radius/std::sqrt(2.0)),
-                                 p+Point<2>(+1,+1) *(radius/std::sqrt(2.0))
-                               };
-
-  const int cell_vertices[5][4] = {{0, 1, 2, 3},
-    {0, 2, 6, 4},
-    {2, 3, 4, 5},
-    {1, 7, 3, 5},
-    {6, 4, 7, 5}
-  };
-
-  std::vector<CellData<2> > cells (5, CellData<2>());
-
-  for (unsigned int i=0; i<5; ++i)
-    {
-      for (unsigned int j=0; j<4; ++j)
-        cells[i].vertices[j] = cell_vertices[i][j];
-      cells[i].material_id = 0;
-    };
+  template <>
+  void
+  cylinder (Triangulation<2> &tria,
+            const double radius,
+            const double half_length)
+  {
+    Point<2> p1 (-half_length, -radius);
+    Point<2> p2 (half_length, radius);
 
-  tria.create_triangulation (
-    std::vector<Point<2> >(&vertices[0], &vertices[8]),
-    cells,
-    SubCellData());       // no boundary information
-}
+    hyper_rectangle(tria, p1, p2, true);
+
+    Triangulation<2>::face_iterator f = tria.begin_face();
+    Triangulation<2>::face_iterator end = tria.end_face();
+    while (f != end)
+      {
+        switch (f->boundary_indicator())
+          {
+          case 0:
+            f->set_boundary_indicator(1);
+            break;
+          case 1:
+            f->set_boundary_indicator(2);
+            break;
+          default:
+            f->set_boundary_indicator(0);
+            break;
+          }
+        ++f;
+      }
+  }
 
 
 
 // Implementation for 2D only
-template<>
-void
-GridGenerator::colorize_hyper_shell (
-  Triangulation<2> &tria,
-  const Point<2> &, const double, const double)
-{
-  // In spite of receiving geometrical
-  // data, we do this only based on
-  // topology.
-
-  // For the mesh based on  cube,
-  // this is highly irregular
-  for (Triangulation<2>::cell_iterator cell = tria.begin();
-       cell != tria.end(); ++cell)
-    {
-      Assert (cell->face(2)->at_boundary(), ExcInternalError());
-      cell->face(2)->set_boundary_indicator(1);
-    }
-}
-
+  template <>
+  void cylinder_shell (Triangulation<2> &,
+                       const double,
+                       const double,
+                       const double,
+                       const unsigned int,
+                       const unsigned int)
+  {
+    Assert (false, ExcNotImplemented());
+  }
 
-template <>
-void GridGenerator::hyper_shell (Triangulation<2> &tria,
-                                 const Point<2>   &center,
-                                 const double      inner_radius,
-                                 const double      outer_radius,
-                                 const size_type   n_cells,
-                                 const bool colorize)
-{
-  Assert ((inner_radius > 0) && (inner_radius < outer_radius),
-          ExcInvalidRadii ());
-
-  const double pi = numbers::PI;
-
-  // determine the number of cells
-  // for the grid. if not provided by
-  // the user determine it such that
-  // the length of each cell on the
-  // median (in the middle between
-  // the two circles) is equal to its
-  // radial extent (which is the
-  // difference between the two
-  // radii)
-  const size_type N = (n_cells == 0 ?
-                       static_cast<size_type>
-                       (std::ceil((2*pi* (outer_radius + inner_radius)/2) /
-                                  (outer_radius - inner_radius))) :
-                       n_cells);
-
-  // set up N vertices on the
-  // outer and N vertices on
-  // the inner circle. the
-  // first N ones are on the
-  // outer one, and all are
-  // numbered counter-clockwise
-  std::vector<Point<2> > vertices(2*N);
-  for (size_type i=0; i<N; ++i)
-    {
-      vertices[i] = Point<2>( std::cos(2*pi * i/N),
-                              std::sin(2*pi * i/N)) * outer_radius;
-      vertices[i+N] = vertices[i] * (inner_radius/outer_radius);
 
-      vertices[i]   += center;
-      vertices[i+N] += center;
+  template <>
+  void
+  half_hyper_ball (Triangulation<2> &tria,
+                   const Point<2>   &p,
+                   const double      radius)
+  {
+    // equilibrate cell sizes at
+    // transition from the inner part
+    // to the radial cells
+    const double a = 1./(1+std::sqrt(2.0));
+    const Point<2> vertices[8] = { p+Point<2>(0,-1) *radius,
+                                   p+Point<2>(+1,-1) *(radius/std::sqrt(2.0)),
+                                   p+Point<2>(0,-1) *(radius/std::sqrt(2.0)*a),
+                                   p+Point<2>(+1,-1) *(radius/std::sqrt(2.0)*a),
+                                   p+Point<2>(0,+1) *(radius/std::sqrt(2.0)*a),
+                                   p+Point<2>(+1,+1) *(radius/std::sqrt(2.0)*a),
+                                   p+Point<2>(0,+1) *radius,
+                                   p+Point<2>(+1,+1) *(radius/std::sqrt(2.0))
+                                 };
+
+    const int cell_vertices[5][4] = {{0, 1, 2, 3},
+      {2, 3, 4, 5},
+      {1, 7, 3, 5},
+      {6, 4, 7, 5}
     };
 
-  std::vector<CellData<2> > cells (N, CellData<2>());
-
-  for (size_type i=0; i<N; ++i)
-    {
-      cells[i].vertices[0] = i;
-      cells[i].vertices[1] = (i+1)%N;
-      cells[i].vertices[2] = N+i;
-      cells[i].vertices[3] = N+((i+1)%N);
-
-      cells[i].material_id = 0;
-    };
+    std::vector<CellData<2> > cells (4, CellData<2>());
 
-  tria.create_triangulation (
-    vertices, cells, SubCellData());
+    for (unsigned int i=0; i<4; ++i)
+      {
+        for (unsigned int j=0; j<4; ++j)
+          cells[i].vertices[j] = cell_vertices[i][j];
+        cells[i].material_id = 0;
+      };
 
-  if (colorize)
-    colorize_hyper_shell(tria, center, inner_radius, outer_radius);
-}
+    tria.create_triangulation (
+      std::vector<Point<2> >(&vertices[0], &vertices[8]),
+      cells,
+      SubCellData());       // no boundary information
 
+    Triangulation<2>::cell_iterator cell = tria.begin();
+    Triangulation<2>::cell_iterator end = tria.end();
 
-// Implementation for 2D only
-template <>
-void
-GridGenerator::cylinder (Triangulation<2> &tria,
-                         const double radius,
-                         const double half_length)
-{
-  Point<2> p1 (-half_length, -radius);
-  Point<2> p2 (half_length, radius);
 
-  hyper_rectangle(tria, p1, p2, true);
+    while (cell != end)
+      {
+        for (unsigned int i=0; i<GeometryInfo<2>::faces_per_cell; ++i)
+          {
+            if (cell->face(i)->boundary_indicator() == numbers::internal_face_boundary_id)
+              continue;
 
-  Triangulation<2>::face_iterator f = tria.begin_face();
-  Triangulation<2>::face_iterator end = tria.end_face();
-  while (f != end)
-    {
-      switch (f->boundary_indicator())
-        {
-        case 0:
-          f->set_boundary_indicator(1);
-          break;
-        case 1:
-          f->set_boundary_indicator(2);
-          break;
-        default:
-          f->set_boundary_indicator(0);
-          break;
-        }
-      ++f;
-    }
-}
+            // If x is zero, then this is part of the plane
+            if (cell->face(i)->center()(0) < p(0)+1.e-5)
+              cell->face(i)->set_boundary_indicator(1);
+          }
+        ++cell;
+      }
+  }
 
 
 
 // Implementation for 2D only
-template <>
-void GridGenerator::cylinder_shell (Triangulation<2> &,
-                                    const double,
-                                    const double,
-                                    const double,
-                                    const size_type,
-                                    const size_type)
-{
-  Assert (false, ExcNotImplemented());
-}
+  template <>
+  void
+  half_hyper_shell (Triangulation<2> &tria,
+                    const Point<2>   &center,
+                    const double      inner_radius,
+                    const double      outer_radius,
+                    const unsigned int   n_cells,
+                    const bool colorize)
+  {
+    Assert ((inner_radius > 0) && (inner_radius < outer_radius),
+            ExcInvalidRadii ());
+
+    const double pi     = numbers::PI;
+    // determine the number of cells
+    // for the grid. if not provided by
+    // the user determine it such that
+    // the length of each cell on the
+    // median (in the middle between
+    // the two circles) is equal to its
+    // radial extent (which is the
+    // difference between the two
+    // radii)
+    const unsigned int N = (n_cells == 0 ?
+                            static_cast<unsigned int>
+                            (std::ceil((pi* (outer_radius + inner_radius)/2) /
+                                       (outer_radius - inner_radius))) :
+                            n_cells);
+
+    // set up N+1 vertices on the
+    // outer and N+1 vertices on
+    // the inner circle. the
+    // first N+1 ones are on the
+    // outer one, and all are
+    // numbered counter-clockwise
+    std::vector<Point<2> > vertices(2*(N+1));
+    for (unsigned int i=0; i<=N; ++i)
+      {
+        // enforce that the x-coordinates
+        // of the first and last point of
+        // each half-circle are exactly
+        // zero (contrary to what we may
+        // compute using the imprecise
+        // value of pi)
+        vertices[i] =  Point<2>( ( (i==0) || (i==N) ?
+                                   0 :
+                                   std::cos(pi * i/N - pi/2) ),
+                                 std::sin(pi * i/N - pi/2)) * outer_radius;
+        vertices[i+N+1] = vertices[i] * (inner_radius/outer_radius);
+
+        vertices[i]     += center;
+        vertices[i+N+1] += center;
+      };
 
 
-template <>
-void
-GridGenerator::half_hyper_ball (Triangulation<2> &tria,
-                                const Point<2>   &p,
-                                const double      radius)
-{
-  // equilibrate cell sizes at
-  // transition from the inner part
-  // to the radial cells
-  const double a = 1./(1+std::sqrt(2.0));
-  const Point<2> vertices[8] = { p+Point<2>(0,-1) *radius,
-                                 p+Point<2>(+1,-1) *(radius/std::sqrt(2.0)),
-                                 p+Point<2>(0,-1) *(radius/std::sqrt(2.0)*a),
-                                 p+Point<2>(+1,-1) *(radius/std::sqrt(2.0)*a),
-                                 p+Point<2>(0,+1) *(radius/std::sqrt(2.0)*a),
-                                 p+Point<2>(+1,+1) *(radius/std::sqrt(2.0)*a),
-                                 p+Point<2>(0,+1) *radius,
-                                 p+Point<2>(+1,+1) *(radius/std::sqrt(2.0))
-                               };
-
-  const int cell_vertices[5][4] = {{0, 1, 2, 3},
-    {2, 3, 4, 5},
-    {1, 7, 3, 5},
-    {6, 4, 7, 5}
-  };
-
-  std::vector<CellData<2> > cells (4, CellData<2>());
-
-  for (unsigned int i=0; i<4; ++i)
-    {
-      for (unsigned int j=0; j<4; ++j)
-        cells[i].vertices[j] = cell_vertices[i][j];
-      cells[i].material_id = 0;
-    };
+    std::vector<CellData<2> > cells (N, CellData<2>());
 
-  tria.create_triangulation (
-    std::vector<Point<2> >(&vertices[0], &vertices[8]),
-    cells,
-    SubCellData());       // no boundary information
+    for (unsigned int i=0; i<N; ++i)
+      {
+        cells[i].vertices[0] = i;
+        cells[i].vertices[1] = (i+1)%(N+1);
+        cells[i].vertices[2] = N+1+i;
+        cells[i].vertices[3] = N+1+((i+1)%(N+1));
 
-  Triangulation<2>::cell_iterator cell = tria.begin();
-  Triangulation<2>::cell_iterator end = tria.end();
+        cells[i].material_id = 0;
+      };
 
+    tria.create_triangulation (vertices, cells, SubCellData());
 
-  while (cell != end)
-    {
-      for (unsigned int i=0; i<GeometryInfo<2>::faces_per_cell; ++i)
-        {
-          if (cell->face(i)->boundary_indicator() == numbers::internal_face_boundary_id)
-            continue;
+    if (colorize)
+      {
+        Triangulation<2>::cell_iterator cell = tria.begin();
+        for (; cell!=tria.end(); ++cell)
+          {
+            cell->face(2)->set_boundary_indicator(1);
+          }
+        tria.begin()->face(0)->set_boundary_indicator(3);
 
-          // If x is zero, then this is part of the plane
-          if (cell->face(i)->center()(0) < p(0)+1.e-5)
-            cell->face(i)->set_boundary_indicator(1);
-        }
-      ++cell;
-    }
-}
+        tria.last()->face(1)->set_boundary_indicator(2);
+      }
+  }
 
 
+  template <>
+  void quarter_hyper_shell (Triangulation<2> &tria,
+                            const Point<2>   &center,
+                            const double      inner_radius,
+                            const double      outer_radius,
+                            const unsigned int   n_cells,
+                            const bool colorize)
+  {
+    Assert ((inner_radius > 0) && (inner_radius < outer_radius),
+            ExcInvalidRadii ());
+
+    const double pi     = numbers::PI;
+    // determine the number of cells
+    // for the grid. if not provided by
+    // the user determine it such that
+    // the length of each cell on the
+    // median (in the middle between
+    // the two circles) is equal to its
+    // radial extent (which is the
+    // difference between the two
+    // radii)
+    const unsigned int N = (n_cells == 0 ?
+                            static_cast<unsigned int>
+                            (std::ceil((pi* (outer_radius + inner_radius)/4) /
+                                       (outer_radius - inner_radius))) :
+                            n_cells);
+
+    // set up N+1 vertices on the
+    // outer and N+1 vertices on
+    // the inner circle. the
+    // first N+1 ones are on the
+    // outer one, and all are
+    // numbered counter-clockwise
+    std::vector<Point<2> > vertices(2*(N+1));
+    for (unsigned int i=0; i<=N; ++i)
+      {
+        // enforce that the x-coordinates
+        // of the last point is exactly
+        // zero (contrary to what we may
+        // compute using the imprecise
+        // value of pi)
+        vertices[i] =  Point<2>( ( (i==N) ?
+                                   0 :
+                                   std::cos(pi * i/N/2) ),
+                                 std::sin(pi * i/N/2)) * outer_radius;
+        vertices[i+N+1] = vertices[i] * (inner_radius/outer_radius);
+
+        vertices[i]     += center;
+        vertices[i+N+1] += center;
+      };
 
-// Implementation for 2D only
-template <>
-void
-GridGenerator::half_hyper_shell (Triangulation<2> &tria,
-                                 const Point<2>   &center,
-                                 const double      inner_radius,
-                                 const double      outer_radius,
-                                 const size_type   n_cells,
-                                 const bool colorize)
-{
-  Assert ((inner_radius > 0) && (inner_radius < outer_radius),
-          ExcInvalidRadii ());
-
-  const double pi     = numbers::PI;
-  // determine the number of cells
-  // for the grid. if not provided by
-  // the user determine it such that
-  // the length of each cell on the
-  // median (in the middle between
-  // the two circles) is equal to its
-  // radial extent (which is the
-  // difference between the two
-  // radii)
-  const size_type N = (n_cells == 0 ?
-                       static_cast<size_type>
-                       (std::ceil((pi* (outer_radius + inner_radius)/2) /
-                                  (outer_radius - inner_radius))) :
-                       n_cells);
-
-  // set up N+1 vertices on the
-  // outer and N+1 vertices on
-  // the inner circle. the
-  // first N+1 ones are on the
-  // outer one, and all are
-  // numbered counter-clockwise
-  std::vector<Point<2> > vertices(2*(N+1));
-  for (size_type i=0; i<=N; ++i)
-    {
-      // enforce that the x-coordinates
-      // of the first and last point of
-      // each half-circle are exactly
-      // zero (contrary to what we may
-      // compute using the imprecise
-      // value of pi)
-      vertices[i] =  Point<2>( ( (i==0) || (i==N) ?
-                                 0 :
-                                 std::cos(pi * i/N - pi/2) ),
-                               std::sin(pi * i/N - pi/2)) * outer_radius;
-      vertices[i+N+1] = vertices[i] * (inner_radius/outer_radius);
-
-      vertices[i]     += center;
-      vertices[i+N+1] += center;
-    };
 
+    std::vector<CellData<2> > cells (N, CellData<2>());
 
-  std::vector<CellData<2> > cells (N, CellData<2>());
+    for (unsigned int i=0; i<N; ++i)
+      {
+        cells[i].vertices[0] = i;
+        cells[i].vertices[1] = (i+1)%(N+1);
+        cells[i].vertices[2] = N+1+i;
+        cells[i].vertices[3] = N+1+((i+1)%(N+1));
 
-  for (size_type i=0; i<N; ++i)
-    {
-      cells[i].vertices[0] = i;
-      cells[i].vertices[1] = (i+1)%(N+1);
-      cells[i].vertices[2] = N+1+i;
-      cells[i].vertices[3] = N+1+((i+1)%(N+1));
+        cells[i].material_id = 0;
+      };
 
-      cells[i].material_id = 0;
-    };
+    tria.create_triangulation (vertices, cells, SubCellData());
 
-  tria.create_triangulation (vertices, cells, SubCellData());
+    if (colorize)
+      {
+        Triangulation<2>::cell_iterator cell = tria.begin();
+        for (; cell!=tria.end(); ++cell)
+          {
+            cell->face(2)->set_boundary_indicator(1);
+          }
+        tria.begin()->face(0)->set_boundary_indicator(3);
 
-  if (colorize)
-    {
-      Triangulation<2>::cell_iterator cell = tria.begin();
-      for (; cell!=tria.end(); ++cell)
-        {
-          cell->face(2)->set_boundary_indicator(1);
-        }
-      tria.begin()->face(0)->set_boundary_indicator(3);
+        tria.last()->face(1)->set_boundary_indicator(2);
+      }
+  }
 
-      tria.last()->face(1)->set_boundary_indicator(2);
-    }
-}
 
 
-template <>
-void GridGenerator::quarter_hyper_shell (Triangulation<2> &tria,
-                                         const Point<2>   &center,
-                                         const double      inner_radius,
-                                         const double      outer_radius,
-                                         const size_type   n_cells,
-                                         const bool colorize)
-{
-  Assert ((inner_radius > 0) && (inner_radius < outer_radius),
-          ExcInvalidRadii ());
-
-  const double pi     = numbers::PI;
-  // determine the number of cells
-  // for the grid. if not provided by
-  // the user determine it such that
-  // the length of each cell on the
-  // median (in the middle between
-  // the two circles) is equal to its
-  // radial extent (which is the
-  // difference between the two
-  // radii)
-  const size_type N = (n_cells == 0 ?
-                       static_cast<size_type>
-                       (std::ceil((pi* (outer_radius + inner_radius)/4) /
-                                  (outer_radius - inner_radius))) :
-                       n_cells);
-
-  // set up N+1 vertices on the
-  // outer and N+1 vertices on
-  // the inner circle. the
-  // first N+1 ones are on the
-  // outer one, and all are
-  // numbered counter-clockwise
-  std::vector<Point<2> > vertices(2*(N+1));
-  for (size_type i=0; i<=N; ++i)
-    {
-      // enforce that the x-coordinates
-      // of the last point is exactly
-      // zero (contrary to what we may
-      // compute using the imprecise
-      // value of pi)
-      vertices[i] =  Point<2>( ( (i==N) ?
-                                 0 :
-                                 std::cos(pi * i/N/2) ),
-                               std::sin(pi * i/N/2)) * outer_radius;
-      vertices[i+N+1] = vertices[i] * (inner_radius/outer_radius);
-
-      vertices[i]     += center;
-      vertices[i+N+1] += center;
+// Implementation for 3D only
+  template <>
+  void hyper_cube_slit (Triangulation<3> &tria,
+                        const double left,
+                        const double right,
+                        const bool colorize)
+  {
+    const double rl2=(right+left)/2;
+    const double len = (right-left)/2.;
+
+    const Point<3> vertices[20] =
+    {
+      Point<3>(left, left , -len/2.),
+      Point<3>(rl2,  left , -len/2.),
+      Point<3>(rl2,  rl2  , -len/2.),
+      Point<3>(left, rl2  , -len/2.),
+      Point<3>(right,left , -len/2.),
+      Point<3>(right,rl2  , -len/2.),
+      Point<3>(rl2,  right, -len/2.),
+      Point<3>(left, right, -len/2.),
+      Point<3>(right,right, -len/2.),
+      Point<3>(rl2,  left , -len/2.),
+      Point<3>(left, left , len/2.),
+      Point<3>(rl2,  left , len/2.),
+      Point<3>(rl2,  rl2  , len/2.),
+      Point<3>(left, rl2  , len/2.),
+      Point<3>(right,left , len/2.),
+      Point<3>(right,rl2  , len/2.),
+      Point<3>(rl2,  right, len/2.),
+      Point<3>(left, right, len/2.),
+      Point<3>(right,right, len/2.),
+      Point<3>(rl2,  left , len/2.)
     };
-
-
-  std::vector<CellData<2> > cells (N, CellData<2>());
-
-  for (size_type i=0; i<N; ++i)
-    {
-      cells[i].vertices[0] = i;
-      cells[i].vertices[1] = (i+1)%(N+1);
-      cells[i].vertices[2] = N+1+i;
-      cells[i].vertices[3] = N+1+((i+1)%(N+1));
-
-      cells[i].material_id = 0;
+    const int cell_vertices[4][8] = { { 0,1,3,2, 10, 11, 13, 12 },
+      { 9,4,2,5, 19,14, 12, 15 },
+      { 3,2,7,6,13,12,17,16 },
+      { 2,5,6,8,12,15,16,18 }
     };
+    std::vector<CellData<3> > cells (4, CellData<3>());
+    for (unsigned int i=0; i<4; ++i)
+      {
+        for (unsigned int j=0; j<8; ++j)
+          cells[i].vertices[j] = cell_vertices[i][j];
+        cells[i].material_id = 0;
+      };
+    tria.create_triangulation (
+      std::vector<Point<3> >(&vertices[0], &vertices[20]),
+      cells,
+      SubCellData());       // no boundary information
 
-  tria.create_triangulation (vertices, cells, SubCellData());
-
-  if (colorize)
-    {
-      Triangulation<2>::cell_iterator cell = tria.begin();
-      for (; cell!=tria.end(); ++cell)
-        {
-          cell->face(2)->set_boundary_indicator(1);
-        }
-      tria.begin()->face(0)->set_boundary_indicator(3);
-
-      tria.last()->face(1)->set_boundary_indicator(2);
-    }
-}
+    if (colorize)
+      {
+        Assert(false, ExcNotImplemented());
+        Triangulation<3>::cell_iterator cell = tria.begin();
+        cell->face(1)->set_boundary_indicator(1);
+        ++cell;
+        cell->face(3)->set_boundary_indicator(2);
+      }
+  }
 
 
 
 // Implementation for 3D only
-template <>
-void GridGenerator::hyper_cube_slit (Triangulation<3> &tria,
-                                     const double left,
-                                     const double right,
-                                     const bool colorize)
-{
-  const double rl2=(right+left)/2;
-  const double len = (right-left)/2.;
-
-  const Point<3> vertices[20] =
+  template <>
+  void enclosed_hyper_cube (Triangulation<3> &tria,
+                            const double        left,
+                            const double        right,
+                            const double        thickness,
+                            const bool          colorize)
   {
-    Point<3>(left, left , -len/2.),
-    Point<3>(rl2,  left , -len/2.),
-    Point<3>(rl2,  rl2  , -len/2.),
-    Point<3>(left, rl2  , -len/2.),
-    Point<3>(right,left , -len/2.),
-    Point<3>(right,rl2  , -len/2.),
-    Point<3>(rl2,  right, -len/2.),
-    Point<3>(left, right, -len/2.),
-    Point<3>(right,right, -len/2.),
-    Point<3>(rl2,  left , -len/2.),
-    Point<3>(left, left , len/2.),
-    Point<3>(rl2,  left , len/2.),
-    Point<3>(rl2,  rl2  , len/2.),
-    Point<3>(left, rl2  , len/2.),
-    Point<3>(right,left , len/2.),
-    Point<3>(right,rl2  , len/2.),
-    Point<3>(rl2,  right, len/2.),
-    Point<3>(left, right, len/2.),
-    Point<3>(right,right, len/2.),
-    Point<3>(rl2,  left , len/2.)
-  };
-  const int cell_vertices[4][8] = { { 0,1,3,2, 10, 11, 13, 12 },
-    { 9,4,2,5, 19,14, 12, 15 },
-    { 3,2,7,6,13,12,17,16 },
-    { 2,5,6,8,12,15,16,18 }
-  };
-  std::vector<CellData<3> > cells (4, CellData<3>());
-  for (unsigned int i=0; i<4; ++i)
-    {
-      for (unsigned int j=0; j<8; ++j)
-        cells[i].vertices[j] = cell_vertices[i][j];
-      cells[i].material_id = 0;
+    Assert(left<right,
+           ExcMessage ("Invalid left-to-right bounds of enclosed hypercube"));
+
+    std::vector<Point<3> > vertices(64);
+    double coords[4];
+    coords[0] = left-thickness;
+    coords[1] = left;
+    coords[2] = right;
+    coords[3] = right+thickness;
+
+    unsigned int k=0;
+    for (unsigned int z=0; z<4; ++z)
+      for (unsigned int y=0; y<4; ++y)
+        for (unsigned int x=0; x<4; ++x)
+          vertices[k++] = Point<3>(coords[x], coords[y], coords[z]);
+
+    const types::material_id materials[27] =
+    {
+      21,20,22,
+      17,16,18,
+      25,24,26,
+      5 , 4, 6,
+      1 , 0, 2,
+      9 , 8,10,
+      37,36,38,
+      33,32,34,
+      41,40,42
     };
-  tria.create_triangulation (
-    std::vector<Point<3> >(&vertices[0], &vertices[20]),
-    cells,
-    SubCellData());       // no boundary information
 
-  if (colorize)
-    {
-      Assert(false, ExcNotImplemented());
-      Triangulation<3>::cell_iterator cell = tria.begin();
-      cell->face(1)->set_boundary_indicator(1);
-      ++cell;
-      cell->face(3)->set_boundary_indicator(2);
-    }
-}
+    std::vector<CellData<3> > cells(27);
+    k = 0;
+    for (unsigned int z=0; z<3; ++z)
+      for (unsigned int y=0; y<3; ++y)
+        for (unsigned int x=0; x<3; ++x)
+          {
+            cells[k].vertices[0] = x+4*y+16*z;
+            cells[k].vertices[1] = x+4*y+16*z+1;
+            cells[k].vertices[2] = x+4*y+16*z+4;
+            cells[k].vertices[3] = x+4*y+16*z+5;
+            cells[k].vertices[4] = x+4*y+16*z+16;
+            cells[k].vertices[5] = x+4*y+16*z+17;
+            cells[k].vertices[6] = x+4*y+16*z+20;
+            cells[k].vertices[7] = x+4*y+16*z+21;
+            if (colorize)
+              cells[k].material_id = materials[k];
+            ++k;
+          }
+    tria.create_triangulation (
+      vertices,
+      cells,
+      SubCellData());       // no boundary information
+  }
 
 
 
-// Implementation for 3D only
-template <>
-void GridGenerator::enclosed_hyper_cube (Triangulation<3> &tria,
-                                         const double        left,
-                                         const double        right,
-                                         const double        thickness,
-                                         const bool          colorize)
-{
-  Assert(left<right,
-         ExcMessage ("Invalid left-to-right bounds of enclosed hypercube"));
-
-  std::vector<Point<3> > vertices(64);
-  double coords[4];
-  coords[0] = left-thickness;
-  coords[1] = left;
-  coords[2] = right;
-  coords[3] = right+thickness;
-
-  unsigned int k=0;
-  for (unsigned int z=0; z<4; ++z)
-    for (unsigned int y=0; y<4; ++y)
-      for (unsigned int x=0; x<4; ++x)
-        vertices[k++] = Point<3>(coords[x], coords[y], coords[z]);
-
-  const types::material_id materials[27] =
+  template <>
+  void truncated_cone (Triangulation<3> &triangulation,
+                       const double radius_0,
+                       const double radius_1,
+                       const double half_length)
   {
-    21,20,22,
-    17,16,18,
-    25,24,26,
-    5 , 4, 6,
-    1 , 0, 2,
-    9 , 8,10,
-    37,36,38,
-    33,32,34,
-    41,40,42
-  };
-
-  std::vector<CellData<3> > cells(27);
-  k = 0;
-  for (unsigned int z=0; z<3; ++z)
-    for (unsigned int y=0; y<3; ++y)
-      for (unsigned int x=0; x<3; ++x)
-        {
-          cells[k].vertices[0] = x+4*y+16*z;
-          cells[k].vertices[1] = x+4*y+16*z+1;
-          cells[k].vertices[2] = x+4*y+16*z+4;
-          cells[k].vertices[3] = x+4*y+16*z+5;
-          cells[k].vertices[4] = x+4*y+16*z+16;
-          cells[k].vertices[5] = x+4*y+16*z+17;
-          cells[k].vertices[6] = x+4*y+16*z+20;
-          cells[k].vertices[7] = x+4*y+16*z+21;
-          if (colorize)
-            cells[k].material_id = materials[k];
-          ++k;
-        }
-  tria.create_triangulation (
-    vertices,
-    cells,
-    SubCellData());       // no boundary information
-}
-
-
-
-template <>
-void GridGenerator::truncated_cone (Triangulation<3> &triangulation,
-                                    const double radius_0,
-                                    const double radius_1,
-                                    const double half_length)
-{
-  // Determine number of cells and vertices
-  const size_type
-  n_cells = static_cast<size_type>(std::floor (half_length /
-                                               std::max (radius_0,
-                                                         radius_1) +
-                                               0.5));
-  const size_type n_vertices = 4 * (n_cells + 1);
-  std::vector<Point<3> > vertices_tmp(n_vertices);
-
-  vertices_tmp[0] = Point<3> (-half_length, 0, -radius_0);
-  vertices_tmp[1] = Point<3> (-half_length, radius_0, 0);
-  vertices_tmp[2] = Point<3> (-half_length, -radius_0, 0);
-  vertices_tmp[3] = Point<3> (-half_length, 0, radius_0);
-
-  const double dx = 2 * half_length / n_cells;
-
-  for (size_type i = 0; i < n_cells; ++i)
-    {
-      vertices_tmp[4 * (i + 1)]
-        = vertices_tmp[4 * i] +
-          Point<3> (dx, 0, 0.5 * (radius_0 - radius_1) * dx / half_length);
-      vertices_tmp[4 * i + 5]
-        = vertices_tmp[4 * i + 1] +
-          Point<3> (dx, 0.5 * (radius_1 - radius_0) * dx / half_length, 0);
-      vertices_tmp[4 * i + 6]
-        = vertices_tmp[4 * i + 2] +
-          Point<3> (dx, 0.5 * (radius_0 - radius_1) * dx / half_length, 0);
-      vertices_tmp[4 * i + 7]
-        = vertices_tmp[4 * i + 3] +
-          Point<3> (dx, 0, 0.5 * (radius_1 - radius_0) * dx / half_length);
-    }
+    // Determine number of cells and vertices
+    const unsigned int
+    n_cells = static_cast<unsigned int>(std::floor (half_length /
+                                                    std::max (radius_0,
+                                                              radius_1) +
+                                                    0.5));
+    const unsigned int n_vertices = 4 * (n_cells + 1);
+    std::vector<Point<3> > vertices_tmp(n_vertices);
+
+    vertices_tmp[0] = Point<3> (-half_length, 0, -radius_0);
+    vertices_tmp[1] = Point<3> (-half_length, radius_0, 0);
+    vertices_tmp[2] = Point<3> (-half_length, -radius_0, 0);
+    vertices_tmp[3] = Point<3> (-half_length, 0, radius_0);
+
+    const double dx = 2 * half_length / n_cells;
+
+    for (unsigned int i = 0; i < n_cells; ++i)
+      {
+        vertices_tmp[4 * (i + 1)]
+          = vertices_tmp[4 * i] +
+            Point<3> (dx, 0, 0.5 * (radius_0 - radius_1) * dx / half_length);
+        vertices_tmp[4 * i + 5]
+          = vertices_tmp[4 * i + 1] +
+            Point<3> (dx, 0.5 * (radius_1 - radius_0) * dx / half_length, 0);
+        vertices_tmp[4 * i + 6]
+          = vertices_tmp[4 * i + 2] +
+            Point<3> (dx, 0.5 * (radius_0 - radius_1) * dx / half_length, 0);
+        vertices_tmp[4 * i + 7]
+          = vertices_tmp[4 * i + 3] +
+            Point<3> (dx, 0, 0.5 * (radius_1 - radius_0) * dx / half_length);
+      }
 
-  const std::vector<Point<3> > vertices (vertices_tmp.begin(),
-                                         vertices_tmp.end());
-  Table<2,unsigned int> cell_vertices(n_cells,GeometryInfo<3>::vertices_per_cell);
+    const std::vector<Point<3> > vertices (vertices_tmp.begin(),
+                                           vertices_tmp.end());
+    Table<2,unsigned int> cell_vertices(n_cells,GeometryInfo<3>::vertices_per_cell);
 
-  for (size_type i = 0; i < n_cells; ++i)
-    for (unsigned int j = 0; j < GeometryInfo<3>::vertices_per_cell; ++j)
-      cell_vertices[i][j] = 4 * i + j;
+    for (unsigned int i = 0; i < n_cells; ++i)
+      for (unsigned int j = 0; j < GeometryInfo<3>::vertices_per_cell; ++j)
+        cell_vertices[i][j] = 4 * i + j;
 
-  std::vector<CellData<3> > cells (n_cells, CellData<3> ());
+    std::vector<CellData<3> > cells (n_cells, CellData<3> ());
 
-  for (size_type i = 0; i < n_cells; ++i)
-    {
-      for (unsigned int j = 0; j < GeometryInfo<3>::vertices_per_cell; ++j)
-        cells[i].vertices[j] = cell_vertices[i][j];
+    for (unsigned int i = 0; i < n_cells; ++i)
+      {
+        for (unsigned int j = 0; j < GeometryInfo<3>::vertices_per_cell; ++j)
+          cells[i].vertices[j] = cell_vertices[i][j];
 
-      cells[i].material_id = 0;
-    }
+        cells[i].material_id = 0;
+      }
 
-  triangulation.create_triangulation (vertices, cells, SubCellData ());
+    triangulation.create_triangulation (vertices, cells, SubCellData ());
 
-  for (Triangulation<3>::cell_iterator cell = triangulation.begin ();
-       cell != triangulation.end (); ++cell)
-    {
-      if (cell->vertex (0) (0) == -half_length)
-        {
-          cell->face (4)->set_boundary_indicator (1);
+    for (Triangulation<3>::cell_iterator cell = triangulation.begin ();
+         cell != triangulation.end (); ++cell)
+      {
+        if (cell->vertex (0) (0) == -half_length)
+          {
+            cell->face (4)->set_boundary_indicator (1);
 
-          for (unsigned int i = 0; i < 4; ++i)
-            cell->line (i)->set_boundary_indicator (0);
-        }
+            for (unsigned int i = 0; i < 4; ++i)
+              cell->line (i)->set_boundary_indicator (0);
+          }
 
-      if (cell->vertex (4) (0) == half_length)
-        {
-          cell->face (5)->set_boundary_indicator (2);
+        if (cell->vertex (4) (0) == half_length)
+          {
+            cell->face (5)->set_boundary_indicator (2);
 
-          for (unsigned int i = 4; i < 8; ++i)
-            cell->line (i)->set_boundary_indicator (0);
-        }
+            for (unsigned int i = 4; i < 8; ++i)
+              cell->line (i)->set_boundary_indicator (0);
+          }
 
-      for (unsigned int i = 0; i < 4; ++i)
-        cell->face (i)->set_boundary_indicator (0);
-    }
-}
+        for (unsigned int i = 0; i < 4; ++i)
+          cell->face (i)->set_boundary_indicator (0);
+      }
+  }
 
 
 // Implementation for 3D only
-template <>
-void
-GridGenerator::hyper_L (Triangulation<3> &tria,
-                        const double      a,
-                        const double      b)
-{
-  // we slice out the top back right
-  // part of the cube
-  const Point<3> vertices[26]
-  =
+  template <>
+  void
+  hyper_L (Triangulation<3> &tria,
+           const double      a,
+           const double      b)
   {
-    // front face of the big cube
-    Point<3> (a,      a,a),
-    Point<3> ((a+b)/2,a,a),
-    Point<3> (b,      a,a),
-    Point<3> (a,      a,(a+b)/2),
-    Point<3> ((a+b)/2,a,(a+b)/2),
-    Point<3> (b,      a,(a+b)/2),
-    Point<3> (a,      a,b),
-    Point<3> ((a+b)/2,a,b),
-    Point<3> (b,      a,b),
-    // middle face of the big cube
-    Point<3> (a,      (a+b)/2,a),
-    Point<3> ((a+b)/2,(a+b)/2,a),
-    Point<3> (b,      (a+b)/2,a),
-    Point<3> (a,      (a+b)/2,(a+b)/2),
-    Point<3> ((a+b)/2,(a+b)/2,(a+b)/2),
-    Point<3> (b,      (a+b)/2,(a+b)/2),
-    Point<3> (a,      (a+b)/2,b),
-    Point<3> ((a+b)/2,(a+b)/2,b),
-    Point<3> (b,      (a+b)/2,b),
-    // back face of the big cube
-    // last (top right) point is missing
-    Point<3> (a,      b,a),
-    Point<3> ((a+b)/2,b,a),
-    Point<3> (b,      b,a),
-    Point<3> (a,      b,(a+b)/2),
-    Point<3> ((a+b)/2,b,(a+b)/2),
-    Point<3> (b,      b,(a+b)/2),
-    Point<3> (a,      b,b),
-    Point<3> ((a+b)/2,b,b)
-  };
-  const int cell_vertices[7][8] = {{0, 1, 9, 10, 3, 4, 12, 13},
-    {1, 2, 10, 11, 4, 5, 13, 14},
-    {3, 4, 12, 13, 6, 7, 15, 16},
-    {4, 5, 13, 14, 7, 8, 16, 17},
-    {9, 10, 18, 19, 12, 13, 21, 22},
-    {10, 11, 19, 20, 13, 14, 22, 23},
-    {12, 13, 21, 22, 15, 16, 24, 25}
-  };
-
-  std::vector<CellData<3> > cells (7, CellData<3>());
-
-  for (unsigned int i=0; i<7; ++i)
-    {
-      for (unsigned int j=0; j<8; ++j)
-        cells[i].vertices[j] = cell_vertices[i][j];
-      cells[i].material_id = 0;
+    // we slice out the top back right
+    // part of the cube
+    const Point<3> vertices[26]
+    =
+    {
+      // front face of the big cube
+      Point<3> (a,      a,a),
+      Point<3> ((a+b)/2,a,a),
+      Point<3> (b,      a,a),
+      Point<3> (a,      a,(a+b)/2),
+      Point<3> ((a+b)/2,a,(a+b)/2),
+      Point<3> (b,      a,(a+b)/2),
+      Point<3> (a,      a,b),
+      Point<3> ((a+b)/2,a,b),
+      Point<3> (b,      a,b),
+      // middle face of the big cube
+      Point<3> (a,      (a+b)/2,a),
+      Point<3> ((a+b)/2,(a+b)/2,a),
+      Point<3> (b,      (a+b)/2,a),
+      Point<3> (a,      (a+b)/2,(a+b)/2),
+      Point<3> ((a+b)/2,(a+b)/2,(a+b)/2),
+      Point<3> (b,      (a+b)/2,(a+b)/2),
+      Point<3> (a,      (a+b)/2,b),
+      Point<3> ((a+b)/2,(a+b)/2,b),
+      Point<3> (b,      (a+b)/2,b),
+      // back face of the big cube
+      // last (top right) point is missing
+      Point<3> (a,      b,a),
+      Point<3> ((a+b)/2,b,a),
+      Point<3> (b,      b,a),
+      Point<3> (a,      b,(a+b)/2),
+      Point<3> ((a+b)/2,b,(a+b)/2),
+      Point<3> (b,      b,(a+b)/2),
+      Point<3> (a,      b,b),
+      Point<3> ((a+b)/2,b,b)
+    };
+    const int cell_vertices[7][8] = {{0, 1, 9, 10, 3, 4, 12, 13},
+      {1, 2, 10, 11, 4, 5, 13, 14},
+      {3, 4, 12, 13, 6, 7, 15, 16},
+      {4, 5, 13, 14, 7, 8, 16, 17},
+      {9, 10, 18, 19, 12, 13, 21, 22},
+      {10, 11, 19, 20, 13, 14, 22, 23},
+      {12, 13, 21, 22, 15, 16, 24, 25}
     };
 
-  tria.create_triangulation (
-    std::vector<Point<3> >(&vertices[0], &vertices[26]),
-    cells,
-    SubCellData());       // no boundary information
-}
+    std::vector<CellData<3> > cells (7, CellData<3>());
+
+    for (unsigned int i=0; i<7; ++i)
+      {
+        for (unsigned int j=0; j<8; ++j)
+          cells[i].vertices[j] = cell_vertices[i][j];
+        cells[i].material_id = 0;
+      };
+
+    tria.create_triangulation (
+      std::vector<Point<3> >(&vertices[0], &vertices[26]),
+      cells,
+      SubCellData());       // no boundary information
+  }
 
 
 
 // Implementation for 3D only
-template <>
-void
-GridGenerator::hyper_ball (Triangulation<3> &tria,
-                           const Point<3>   &p,
-                           const double radius)
-{
-  const double a = 1./(1+std::sqrt(3.0)); // equilibrate cell sizes at transition
-  // from the inner part to the radial
-  // cells
-  const unsigned int n_vertices = 16;
-  const Point<3> vertices[n_vertices]
-  =
+  template <>
+  void
+  hyper_ball (Triangulation<3> &tria,
+              const Point<3>   &p,
+              const double radius)
   {
-    // first the vertices of the inner
-    // cell
-    p+Point<3>(-1,-1,-1) *(radius/std::sqrt(3.0)*a),
-    p+Point<3>(+1,-1,-1) *(radius/std::sqrt(3.0)*a),
-    p+Point<3>(+1,-1,+1) *(radius/std::sqrt(3.0)*a),
-    p+Point<3>(-1,-1,+1) *(radius/std::sqrt(3.0)*a),
-    p+Point<3>(-1,+1,-1) *(radius/std::sqrt(3.0)*a),
-    p+Point<3>(+1,+1,-1) *(radius/std::sqrt(3.0)*a),
-    p+Point<3>(+1,+1,+1) *(radius/std::sqrt(3.0)*a),
-    p+Point<3>(-1,+1,+1) *(radius/std::sqrt(3.0)*a),
-    // now the eight vertices at
-    // the outer sphere
-    p+Point<3>(-1,-1,-1) *(radius/std::sqrt(3.0)),
-    p+Point<3>(+1,-1,-1) *(radius/std::sqrt(3.0)),
-    p+Point<3>(+1,-1,+1) *(radius/std::sqrt(3.0)),
-    p+Point<3>(-1,-1,+1) *(radius/std::sqrt(3.0)),
-    p+Point<3>(-1,+1,-1) *(radius/std::sqrt(3.0)),
-    p+Point<3>(+1,+1,-1) *(radius/std::sqrt(3.0)),
-    p+Point<3>(+1,+1,+1) *(radius/std::sqrt(3.0)),
-    p+Point<3>(-1,+1,+1) *(radius/std::sqrt(3.0)),
-  };
-
-  // one needs to draw the seven cubes to
-  // understand what's going on here
-  const unsigned int n_cells = 7;
-  const int cell_vertices[n_cells][8] = {{0, 1, 4, 5, 3, 2, 7, 6}, // center
-    {8, 9, 12, 13, 0, 1, 4, 5}, // bottom
-    {9, 13, 1, 5, 10, 14, 2, 6}, // right
-    {11, 10, 3, 2, 15, 14, 7, 6}, // top
-    {8, 0, 12, 4, 11, 3, 15, 7}, // left
-    {8, 9, 0, 1, 11, 10, 3, 2}, // front
-    {12, 4, 13, 5, 15, 7, 14, 6}
-  }; // back
-
-  std::vector<CellData<3> > cells (n_cells, CellData<3>());
-
-  for (unsigned int i=0; i<n_cells; ++i)
-    {
-      for (unsigned int j=0; j<GeometryInfo<3>::vertices_per_cell; ++j)
-        cells[i].vertices[j] = cell_vertices[i][j];
-      cells[i].material_id = 0;
+    const double a = 1./(1+std::sqrt(3.0)); // equilibrate cell sizes at transition
+    // from the inner part to the radial
+    // cells
+    const unsigned int n_vertices = 16;
+    const Point<3> vertices[n_vertices]
+    =
+    {
+      // first the vertices of the inner
+      // cell
+      p+Point<3>(-1,-1,-1) *(radius/std::sqrt(3.0)*a),
+      p+Point<3>(+1,-1,-1) *(radius/std::sqrt(3.0)*a),
+      p+Point<3>(+1,-1,+1) *(radius/std::sqrt(3.0)*a),
+      p+Point<3>(-1,-1,+1) *(radius/std::sqrt(3.0)*a),
+      p+Point<3>(-1,+1,-1) *(radius/std::sqrt(3.0)*a),
+      p+Point<3>(+1,+1,-1) *(radius/std::sqrt(3.0)*a),
+      p+Point<3>(+1,+1,+1) *(radius/std::sqrt(3.0)*a),
+      p+Point<3>(-1,+1,+1) *(radius/std::sqrt(3.0)*a),
+      // now the eight vertices at
+      // the outer sphere
+      p+Point<3>(-1,-1,-1) *(radius/std::sqrt(3.0)),
+      p+Point<3>(+1,-1,-1) *(radius/std::sqrt(3.0)),
+      p+Point<3>(+1,-1,+1) *(radius/std::sqrt(3.0)),
+      p+Point<3>(-1,-1,+1) *(radius/std::sqrt(3.0)),
+      p+Point<3>(-1,+1,-1) *(radius/std::sqrt(3.0)),
+      p+Point<3>(+1,+1,-1) *(radius/std::sqrt(3.0)),
+      p+Point<3>(+1,+1,+1) *(radius/std::sqrt(3.0)),
+      p+Point<3>(-1,+1,+1) *(radius/std::sqrt(3.0)),
     };
 
-  tria.create_triangulation (
-    std::vector<Point<3> >(&vertices[0], &vertices[n_vertices]),
-    cells,
-    SubCellData());       // no boundary information
-}
+    // one needs to draw the seven cubes to
+    // understand what's going on here
+    const unsigned int n_cells = 7;
+    const int cell_vertices[n_cells][8] = {{0, 1, 4, 5, 3, 2, 7, 6}, // center
+      {8, 9, 12, 13, 0, 1, 4, 5}, // bottom
+      {9, 13, 1, 5, 10, 14, 2, 6}, // right
+      {11, 10, 3, 2, 15, 14, 7, 6}, // top
+      {8, 0, 12, 4, 11, 3, 15, 7}, // left
+      {8, 9, 0, 1, 11, 10, 3, 2}, // front
+      {12, 4, 13, 5, 15, 7, 14, 6}
+    }; // back
+
+    std::vector<CellData<3> > cells (n_cells, CellData<3>());
+
+    for (unsigned int i=0; i<n_cells; ++i)
+      {
+        for (unsigned int j=0; j<GeometryInfo<3>::vertices_per_cell; ++j)
+          cells[i].vertices[j] = cell_vertices[i][j];
+        cells[i].material_id = 0;
+      };
+
+    tria.create_triangulation (
+      std::vector<Point<3> >(&vertices[0], &vertices[n_vertices]),
+      cells,
+      SubCellData());       // no boundary information
+  }
 
 
 
 // Implementation for 3D only
-template <>
-void
-GridGenerator::cylinder (Triangulation<3> &tria,
-                         const double radius,
-                         const double half_length)
-{
-  // Copy the base from hyper_ball<3>
-  // and transform it to yz
-  const double d = radius/std::sqrt(2.0);
-  const double a = d/(1+std::sqrt(2.0));
-  Point<3> vertices[24] =
+  template <>
+  void
+  cylinder (Triangulation<3> &tria,
+            const double radius,
+            const double half_length)
   {
-    Point<3>(-d, -half_length,-d),
-    Point<3>( d, -half_length,-d),
-    Point<3>(-a, -half_length,-a),
-    Point<3>( a, -half_length,-a),
-    Point<3>(-a, -half_length, a),
-    Point<3>( a, -half_length, a),
-    Point<3>(-d, -half_length, d),
-    Point<3>( d, -half_length, d),
-    Point<3>(-d, 0,-d),
-    Point<3>( d, 0,-d),
-    Point<3>(-a, 0,-a),
-    Point<3>( a, 0,-a),
-    Point<3>(-a, 0, a),
-    Point<3>( a, 0, a),
-    Point<3>(-d, 0, d),
-    Point<3>( d, 0, d),
-    Point<3>(-d, half_length,-d),
-    Point<3>( d, half_length,-d),
-    Point<3>(-a, half_length,-a),
-    Point<3>( a, half_length,-a),
-    Point<3>(-a, half_length, a),
-    Point<3>( a, half_length, a),
-    Point<3>(-d, half_length, d),
-    Point<3>( d, half_length, d),
-  };
-  // Turn cylinder such that y->x
-  for (unsigned int i=0; i<24; ++i)
-    {
-      const double h = vertices[i](1);
-      vertices[i](1) = -vertices[i](0);
-      vertices[i](0) = h;
-    }
+    // Copy the base from hyper_ball<3>
+    // and transform it to yz
+    const double d = radius/std::sqrt(2.0);
+    const double a = d/(1+std::sqrt(2.0));
+    Point<3> vertices[24] =
+    {
+      Point<3>(-d, -half_length,-d),
+      Point<3>( d, -half_length,-d),
+      Point<3>(-a, -half_length,-a),
+      Point<3>( a, -half_length,-a),
+      Point<3>(-a, -half_length, a),
+      Point<3>( a, -half_length, a),
+      Point<3>(-d, -half_length, d),
+      Point<3>( d, -half_length, d),
+      Point<3>(-d, 0,-d),
+      Point<3>( d, 0,-d),
+      Point<3>(-a, 0,-a),
+      Point<3>( a, 0,-a),
+      Point<3>(-a, 0, a),
+      Point<3>( a, 0, a),
+      Point<3>(-d, 0, d),
+      Point<3>( d, 0, d),
+      Point<3>(-d, half_length,-d),
+      Point<3>( d, half_length,-d),
+      Point<3>(-a, half_length,-a),
+      Point<3>( a, half_length,-a),
+      Point<3>(-a, half_length, a),
+      Point<3>( a, half_length, a),
+      Point<3>(-d, half_length, d),
+      Point<3>( d, half_length, d),
+    };
+    // Turn cylinder such that y->x
+    for (unsigned int i=0; i<24; ++i)
+      {
+        const double h = vertices[i](1);
+        vertices[i](1) = -vertices[i](0);
+        vertices[i](0) = h;
+      }
 
-  int cell_vertices[10][8] =
-  {
-    {0, 1, 8, 9, 2, 3, 10, 11},
-    {0, 2, 8, 10, 6, 4, 14, 12},
-    {2, 3, 10, 11, 4, 5, 12, 13},
-    {1, 7, 9, 15, 3, 5, 11, 13},
-    {6, 4, 14, 12, 7, 5, 15, 13}
-  };
-  for (unsigned int i=0; i<5; ++i)
-    for (unsigned int j=0; j<8; ++j)
-      cell_vertices[i+5][j] = cell_vertices[i][j]+8;
-
-  std::vector<CellData<3> > cells (10, CellData<3>());
-
-  for (unsigned int i=0; i<10; ++i)
+    int cell_vertices[10][8] =
     {
-      for (unsigned int j=0; j<8; ++j)
-        cells[i].vertices[j] = cell_vertices[i][j];
-      cells[i].material_id = 0;
+      {0, 1, 8, 9, 2, 3, 10, 11},
+      {0, 2, 8, 10, 6, 4, 14, 12},
+      {2, 3, 10, 11, 4, 5, 12, 13},
+      {1, 7, 9, 15, 3, 5, 11, 13},
+      {6, 4, 14, 12, 7, 5, 15, 13}
     };
+    for (unsigned int i=0; i<5; ++i)
+      for (unsigned int j=0; j<8; ++j)
+        cell_vertices[i+5][j] = cell_vertices[i][j]+8;
 
-  tria.create_triangulation (
-    std::vector<Point<3> >(&vertices[0], &vertices[24]),
-    cells,
-    SubCellData());       // no boundary information
-
-  // set boundary indicators for the
-  // faces at the ends to 1 and 2,
-  // respectively. note that we also
-  // have to deal with those lines
-  // that are purely in the interior
-  // of the ends. we determine whether
-  // an edge is purely in the
-  // interior if one of its vertices
-  // is at coordinates '+-a' as set
-  // above
-  Triangulation<3>::cell_iterator cell = tria.begin();
-  Triangulation<3>::cell_iterator end = tria.end();
-
-  for (; cell != end; ++cell)
-    for (unsigned int i=0; i<GeometryInfo<3>::faces_per_cell; ++i)
-      if (cell->at_boundary(i))
-        {
-          if (cell->face(i)->center()(0) > half_length-1.e-5)
-            {
-              cell->face(i)->set_boundary_indicator(2);
-
-              for (unsigned int e=0; e<GeometryInfo<3>::lines_per_face; ++e)
-                if ((std::fabs(cell->face(i)->line(e)->vertex(0)[1]) == a) ||
-                    (std::fabs(cell->face(i)->line(e)->vertex(0)[2]) == a) ||
-                    (std::fabs(cell->face(i)->line(e)->vertex(1)[1]) == a) ||
-                    (std::fabs(cell->face(i)->line(e)->vertex(1)[2]) == a))
-                  cell->face(i)->line(e)->set_boundary_indicator(2);
-            }
-          else if (cell->face(i)->center()(0) < -half_length+1.e-5)
-            {
-              cell->face(i)->set_boundary_indicator(1);
+    std::vector<CellData<3> > cells (10, CellData<3>());
 
-              for (unsigned int e=0; e<GeometryInfo<3>::lines_per_face; ++e)
-                if ((std::fabs(cell->face(i)->line(e)->vertex(0)[1]) == a) ||
-                    (std::fabs(cell->face(i)->line(e)->vertex(0)[2]) == a) ||
-                    (std::fabs(cell->face(i)->line(e)->vertex(1)[1]) == a) ||
-                    (std::fabs(cell->face(i)->line(e)->vertex(1)[2]) == a))
-                  cell->face(i)->line(e)->set_boundary_indicator(1);
-            }
-        }
-}
+    for (unsigned int i=0; i<10; ++i)
+      {
+        for (unsigned int j=0; j<8; ++j)
+          cells[i].vertices[j] = cell_vertices[i][j];
+        cells[i].material_id = 0;
+      };
+
+    tria.create_triangulation (
+      std::vector<Point<3> >(&vertices[0], &vertices[24]),
+      cells,
+      SubCellData());       // no boundary information
+
+    // set boundary indicators for the
+    // faces at the ends to 1 and 2,
+    // respectively. note that we also
+    // have to deal with those lines
+    // that are purely in the interior
+    // of the ends. we determine whether
+    // an edge is purely in the
+    // interior if one of its vertices
+    // is at coordinates '+-a' as set
+    // above
+    Triangulation<3>::cell_iterator cell = tria.begin();
+    Triangulation<3>::cell_iterator end = tria.end();
+
+    for (; cell != end; ++cell)
+      for (unsigned int i=0; i<GeometryInfo<3>::faces_per_cell; ++i)
+        if (cell->at_boundary(i))
+          {
+            if (cell->face(i)->center()(0) > half_length-1.e-5)
+              {
+                cell->face(i)->set_boundary_indicator(2);
+
+                for (unsigned int e=0; e<GeometryInfo<3>::lines_per_face; ++e)
+                  if ((std::fabs(cell->face(i)->line(e)->vertex(0)[1]) == a) ||
+                      (std::fabs(cell->face(i)->line(e)->vertex(0)[2]) == a) ||
+                      (std::fabs(cell->face(i)->line(e)->vertex(1)[1]) == a) ||
+                      (std::fabs(cell->face(i)->line(e)->vertex(1)[2]) == a))
+                    cell->face(i)->line(e)->set_boundary_indicator(2);
+              }
+            else if (cell->face(i)->center()(0) < -half_length+1.e-5)
+              {
+                cell->face(i)->set_boundary_indicator(1);
+
+                for (unsigned int e=0; e<GeometryInfo<3>::lines_per_face; ++e)
+                  if ((std::fabs(cell->face(i)->line(e)->vertex(0)[1]) == a) ||
+                      (std::fabs(cell->face(i)->line(e)->vertex(0)[2]) == a) ||
+                      (std::fabs(cell->face(i)->line(e)->vertex(1)[1]) == a) ||
+                      (std::fabs(cell->face(i)->line(e)->vertex(1)[2]) == a))
+                    cell->face(i)->line(e)->set_boundary_indicator(1);
+              }
+          }
+  }
 
 
 
 // Implementation for 3D only
-template <>
-void
-GridGenerator::half_hyper_ball (Triangulation<3> &tria,
-                                const Point<3> &center,
-                                const double radius)
-{
-  // These are for the two lower squares
-  const double d = radius/std::sqrt(2.0);
-  const double a = d/(1+std::sqrt(2.0));
-  // These are for the two upper square
-  const double b = a/2.0;
-  const double c = d/2.0;
-  // And so are these
-  const double hb = radius*std::sqrt(3.0)/4.0;
-  const double hc = radius*std::sqrt(3.0)/2.0;
-
-  Point<3> vertices[16] =
+  template <>
+  void
+  half_hyper_ball (Triangulation<3> &tria,
+                   const Point<3> &center,
+                   const double radius)
   {
-    center+Point<3>( 0,  d, -d),
-    center+Point<3>( 0, -d, -d),
-    center+Point<3>( 0,  a, -a),
-    center+Point<3>( 0, -a, -a),
-    center+Point<3>( 0,  a,  a),
-    center+Point<3>( 0, -a,  a),
-    center+Point<3>( 0,  d,  d),
-    center+Point<3>( 0, -d,  d),
-
-    center+Point<3>(hc,  c, -c),
-    center+Point<3>(hc, -c, -c),
-    center+Point<3>(hb,  b, -b),
-    center+Point<3>(hb, -b, -b),
-    center+Point<3>(hb,  b,  b),
-    center+Point<3>(hb, -b,  b),
-    center+Point<3>(hc,  c,  c),
-    center+Point<3>(hc, -c,  c),
-  };
-
-  int cell_vertices[6][8] =
-  {
-    {0, 1, 8, 9, 2, 3, 10, 11},
-    {0, 2, 8, 10, 6, 4, 14, 12},
-    {2, 3, 10, 11, 4, 5, 12, 13},
-    {1, 7, 9, 15, 3, 5, 11, 13},
-    {6, 4, 14, 12, 7, 5, 15, 13},
-    {8, 10, 9, 11, 14, 12, 15, 13}
-  };
-
-  std::vector<CellData<3> > cells (6, CellData<3>());
+    // These are for the two lower squares
+    const double d = radius/std::sqrt(2.0);
+    const double a = d/(1+std::sqrt(2.0));
+    // These are for the two upper square
+    const double b = a/2.0;
+    const double c = d/2.0;
+    // And so are these
+    const double hb = radius*std::sqrt(3.0)/4.0;
+    const double hc = radius*std::sqrt(3.0)/2.0;
+
+    Point<3> vertices[16] =
+    {
+      center+Point<3>( 0,  d, -d),
+      center+Point<3>( 0, -d, -d),
+      center+Point<3>( 0,  a, -a),
+      center+Point<3>( 0, -a, -a),
+      center+Point<3>( 0,  a,  a),
+      center+Point<3>( 0, -a,  a),
+      center+Point<3>( 0,  d,  d),
+      center+Point<3>( 0, -d,  d),
+
+      center+Point<3>(hc,  c, -c),
+      center+Point<3>(hc, -c, -c),
+      center+Point<3>(hb,  b, -b),
+      center+Point<3>(hb, -b, -b),
+      center+Point<3>(hb,  b,  b),
+      center+Point<3>(hb, -b,  b),
+      center+Point<3>(hc,  c,  c),
+      center+Point<3>(hc, -c,  c),
+    };
 
-  for (unsigned int i=0; i<6; ++i)
+    int cell_vertices[6][8] =
     {
-      for (unsigned int j=0; j<8; ++j)
-        cells[i].vertices[j] = cell_vertices[i][j];
-      cells[i].material_id = 0;
+      {0, 1, 8, 9, 2, 3, 10, 11},
+      {0, 2, 8, 10, 6, 4, 14, 12},
+      {2, 3, 10, 11, 4, 5, 12, 13},
+      {1, 7, 9, 15, 3, 5, 11, 13},
+      {6, 4, 14, 12, 7, 5, 15, 13},
+      {8, 10, 9, 11, 14, 12, 15, 13}
     };
 
-  tria.create_triangulation (
-    std::vector<Point<3> >(&vertices[0], &vertices[16]),
-    cells,
-    SubCellData());       // no boundary information
+    std::vector<CellData<3> > cells (6, CellData<3>());
 
-  Triangulation<3>::cell_iterator cell = tria.begin();
-  Triangulation<3>::cell_iterator end = tria.end();
-
-  // go over all faces. for the ones on the flat face, set boundary
-  // indicator for face and edges to one; the rest will remain at
-  // zero but we have to pay attention to those edges that are
-  // at the perimeter of the flat face since they should not be
-  // set to one
-  while (cell != end)
-    {
-      for (unsigned int i=0; i<GeometryInfo<3>::faces_per_cell; ++i)
-        {
-          if (!cell->at_boundary(i))
-            continue;
+    for (unsigned int i=0; i<6; ++i)
+      {
+        for (unsigned int j=0; j<8; ++j)
+          cells[i].vertices[j] = cell_vertices[i][j];
+        cells[i].material_id = 0;
+      };
 
-          // If the center is on the plane x=0, this is a planar element. set
-          // its boundary indicator. also set the boundary indicators of the
-          // bounding faces unless both vertices are on the perimeter
-          if (cell->face(i)->center()(0) < center(0)+1.e-5*radius)
-            {
-              cell->face(i)->set_boundary_indicator(1);
-              for (unsigned int j=0; j<GeometryInfo<3>::lines_per_face; ++j)
-                {
-                  const Point<3> vertices[2]
-                    = { cell->face(i)->line(j)->vertex(0),
-                        cell->face(i)->line(j)->vertex(1)
-                      };
-                  if ((std::fabs(vertices[0].distance(center)-radius) >
-                       1e-5*radius)
-                      ||
-                      (std::fabs(vertices[1].distance(center)-radius) >
-                       1e-5*radius))
-                    cell->face(i)->line(j)->set_boundary_indicator(1);
-                }
-            }
-        }
-      ++cell;
-    }
-}
+    tria.create_triangulation (
+      std::vector<Point<3> >(&vertices[0], &vertices[16]),
+      cells,
+      SubCellData());       // no boundary information
 
-// Implementation for 3D only
-template<>
-void
-GridGenerator::
-colorize_hyper_shell (Triangulation<3> &tria,
-                      const Point<3> &,
-                      const double,
-                      const double)
-{
-  // the following uses a good amount
-  // of knowledge about the
-  // orientation of cells. this is
-  // probably not good style...
-  if (tria.n_cells() == 6)
-    {
-      Triangulation<3>::cell_iterator cell = tria.begin();
+    Triangulation<3>::cell_iterator cell = tria.begin();
+    Triangulation<3>::cell_iterator end = tria.end();
 
-      cell->face(4)->set_boundary_indicator(1);
-      Assert (cell->face(4)->at_boundary(), ExcInternalError());
-
-      (++cell)->face(2)->set_boundary_indicator(1);
-      Assert (cell->face(2)->at_boundary(), ExcInternalError());
-
-      (++cell)->face(2)->set_boundary_indicator(1);
-      Assert (cell->face(2)->at_boundary(), ExcInternalError());
-
-      (++cell)->face(0)->set_boundary_indicator(1);
-      Assert (cell->face(0)->at_boundary(), ExcInternalError());
-
-      (++cell)->face(2)->set_boundary_indicator(1);
-      Assert (cell->face(2)->at_boundary(), ExcInternalError());
-
-      (++cell)->face(0)->set_boundary_indicator(1);
-      Assert (cell->face(0)->at_boundary(), ExcInternalError());
-    }
-  else if (tria.n_cells() == 12)
-    {
-      // again use some internal
-      // knowledge
-      for (Triangulation<3>::cell_iterator cell = tria.begin();
-           cell != tria.end(); ++cell)
-        {
-          Assert (cell->face(5)->at_boundary(), ExcInternalError());
-          cell->face(5)->set_boundary_indicator(1);
-        }
-    }
-  else if (tria.n_cells() == 96)
-    {
-      // the 96-cell hypershell is
-      // based on a once refined
-      // 12-cell mesh. consequently,
-      // since the outer faces all
-      // are face_no==5 above, so
-      // they are here (unless they
-      // are in the interior). Use
-      // this to assign boundary
-      // indicators, but also make
-      // sure that we encounter
-      // exactly 48 such faces
-      unsigned int count = 0;
-      for (Triangulation<3>::cell_iterator cell = tria.begin();
-           cell != tria.end(); ++cell)
-        if (cell->face(5)->at_boundary())
+    // go over all faces. for the ones on the flat face, set boundary
+    // indicator for face and edges to one; the rest will remain at
+    // zero but we have to pay attention to those edges that are
+    // at the perimeter of the flat face since they should not be
+    // set to one
+    while (cell != end)
+      {
+        for (unsigned int i=0; i<GeometryInfo<3>::faces_per_cell; ++i)
           {
-            cell->face(5)->set_boundary_indicator(1);
-            ++count;
-          }
-      Assert (count == 48, ExcInternalError());
-    }
-  else
-    Assert (false, ExcNotImplemented());
-}
-
+            if (!cell->at_boundary(i))
+              continue;
 
+            // If the center is on the plane x=0, this is a planar element. set
+            // its boundary indicator. also set the boundary indicators of the
+            // bounding faces unless both vertices are on the perimeter
+            if (cell->face(i)->center()(0) < center(0)+1.e-5*radius)
+              {
+                cell->face(i)->set_boundary_indicator(1);
+                for (unsigned int j=0; j<GeometryInfo<3>::lines_per_face; ++j)
+                  {
+                    const Point<3> vertices[2]
+                      = { cell->face(i)->line(j)->vertex(0),
+                          cell->face(i)->line(j)->vertex(1)
+                        };
+                    if ((std::fabs(vertices[0].distance(center)-radius) >
+                         1e-5*radius)
+                        ||
+                        (std::fabs(vertices[1].distance(center)-radius) >
+                         1e-5*radius))
+                      cell->face(i)->line(j)->set_boundary_indicator(1);
+                  }
+              }
+          }
+        ++cell;
+      }
+  }
 
-template <>
-void
-GridGenerator::hyper_shell (Triangulation<3> &tria,
-                            const Point<3> &p,
-                            const double inner_radius,
-                            const double outer_radius,
-                            const size_type n,
-                            const bool colorize)
-{
-  Assert ((inner_radius > 0) && (inner_radius < outer_radius),
-          ExcInvalidRadii ());
 
-  const double irad = inner_radius/std::sqrt(3.0);
-  const double orad = outer_radius/std::sqrt(3.0);
-  std::vector<Point<3> > vertices;
-  std::vector<CellData<3> > cells;
+  template <>
+  void
+  hyper_shell (Triangulation<3> &tria,
+               const Point<3> &p,
+               const double inner_radius,
+               const double outer_radius,
+               const unsigned int n,
+               const bool colorize)
+  {
+    Assert ((inner_radius > 0) && (inner_radius < outer_radius),
+            ExcInvalidRadii ());
 
-  // Start with the shell bounded by
-  // two nested cubes
-  if (n == 6)
-    {
-      for (unsigned int i=0; i<8; ++i)
-        vertices.push_back(p+hexahedron[i]*irad);
-      for (unsigned int i=0; i<8; ++i)
-        vertices.push_back(p+hexahedron[i]*orad);
+    const double irad = inner_radius/std::sqrt(3.0);
+    const double orad = outer_radius/std::sqrt(3.0);
+    std::vector<Point<3> > vertices;
+    std::vector<CellData<3> > cells;
 
-      const unsigned int n_cells = 6;
-      const int cell_vertices[n_cells][8] =
+    // Start with the shell bounded by
+    // two nested cubes
+    if (n == 6)
       {
-        {8, 9, 10, 11, 0, 1, 2, 3}, // bottom
-        {9, 11, 1, 3, 13, 15, 5, 7}, // right
-        {12, 13, 4, 5, 14, 15, 6, 7}, // top
-        {8, 0, 10, 2, 12, 4, 14, 6}, // left
-        {8, 9, 0, 1, 12, 13, 4, 5}, // front
-        {10, 2, 11, 3, 14, 6, 15, 7}
-      }; // back
+        for (unsigned int i=0; i<8; ++i)
+          vertices.push_back(p+hexahedron[i]*irad);
+        for (unsigned int i=0; i<8; ++i)
+          vertices.push_back(p+hexahedron[i]*orad);
 
-      cells.resize(n_cells, CellData<3>());
-
-      for (unsigned int i=0; i<n_cells; ++i)
+        const unsigned int n_cells = 6;
+        const int cell_vertices[n_cells][8] =
         {
-          for (unsigned int j=0; j<GeometryInfo<3>::vertices_per_cell; ++j)
-            cells[i].vertices[j] = cell_vertices[i][j];
-          cells[i].material_id = 0;
-        }
+          {8, 9, 10, 11, 0, 1, 2, 3}, // bottom
+          {9, 11, 1, 3, 13, 15, 5, 7}, // right
+          {12, 13, 4, 5, 14, 15, 6, 7}, // top
+          {8, 0, 10, 2, 12, 4, 14, 6}, // left
+          {8, 9, 0, 1, 12, 13, 4, 5}, // front
+          {10, 2, 11, 3, 14, 6, 15, 7}
+        }; // back
 
-      tria.create_triangulation (vertices, cells, SubCellData());
-    }
-  // A more regular subdivision can
-  // be obtained by two nested
-  // rhombic dodecahedra
-  else if (n == 12)
-    {
-      for (unsigned int i=0; i<8; ++i)
-        vertices.push_back(p+hexahedron[i]*irad);
-      for (unsigned int i=0; i<6; ++i)
-        vertices.push_back(p+octahedron[i]*inner_radius);
-      for (unsigned int i=0; i<8; ++i)
-        vertices.push_back(p+hexahedron[i]*orad);
-      for (unsigned int i=0; i<6; ++i)
-        vertices.push_back(p+octahedron[i]*outer_radius);
-
-      const unsigned int n_cells = 12;
-      const unsigned int rhombi[n_cells][4] =
-      {
-        { 10,  4,  0,  8},
-        {  4, 13,  8,  6},
-        { 10,  5,  4, 13},
-        {  1,  9, 10,  5},
-        {  9,  7,  5, 13},
-        {  7, 11, 13,  6},
-        {  9,  3,  7, 11},
-        {  1, 12,  9,  3},
-        { 12,  2,  3, 11},
-        {  2,  8, 11,  6},
-        { 12,  0,  2,  8},
-        {  1, 10, 12,  0}
-      };
+        cells.resize(n_cells, CellData<3>());
 
-      cells.resize(n_cells, CellData<3>());
+        for (unsigned int i=0; i<n_cells; ++i)
+          {
+            for (unsigned int j=0; j<GeometryInfo<3>::vertices_per_cell; ++j)
+              cells[i].vertices[j] = cell_vertices[i][j];
+            cells[i].material_id = 0;
+          }
 
-      for (unsigned int i=0; i<n_cells; ++i)
+        tria.create_triangulation (vertices, cells, SubCellData());
+      }
+    // A more regular subdivision can
+    // be obtained by two nested
+    // rhombic dodecahedra
+    else if (n == 12)
+      {
+        for (unsigned int i=0; i<8; ++i)
+          vertices.push_back(p+hexahedron[i]*irad);
+        for (unsigned int i=0; i<6; ++i)
+          vertices.push_back(p+octahedron[i]*inner_radius);
+        for (unsigned int i=0; i<8; ++i)
+          vertices.push_back(p+hexahedron[i]*orad);
+        for (unsigned int i=0; i<6; ++i)
+          vertices.push_back(p+octahedron[i]*outer_radius);
+
+        const unsigned int n_cells = 12;
+        const unsigned int rhombi[n_cells][4] =
         {
-          for (unsigned int j=0; j<4; ++j)
-            {
-              cells[i].vertices[j  ] = rhombi[i][j];
-              cells[i].vertices[j+4] = rhombi[i][j] + 14;
-            }
-          cells[i].material_id = 0;
-        }
+          { 10,  4,  0,  8},
+          {  4, 13,  8,  6},
+          { 10,  5,  4, 13},
+          {  1,  9, 10,  5},
+          {  9,  7,  5, 13},
+          {  7, 11, 13,  6},
+          {  9,  3,  7, 11},
+          {  1, 12,  9,  3},
+          { 12,  2,  3, 11},
+          {  2,  8, 11,  6},
+          { 12,  0,  2,  8},
+          {  1, 10, 12,  0}
+        };
 
-      tria.create_triangulation (vertices, cells, SubCellData());
-    }
-  else if (n == 96)
-    {
-      // create a triangulation based on the
-      // 12-cell one where we refine the mesh
-      // once and then re-arrange all
-      // interior nodes so that the mesh is
-      // the least distorted
-      HyperShellBoundary<3> boundary (p);
-      Triangulation<3> tmp;
-      GridGenerator::hyper_shell (tmp, p, inner_radius, outer_radius, 12);
-      tmp.set_boundary(0, boundary);
-      tmp.set_boundary(1, boundary);
-      tmp.refine_global (1);
-
-      // let's determine the distance at
-      // which the interior nodes should be
-      // from the center. let's say we
-      // measure distances in multiples of
-      // outer_radius and call
-      // r=inner_radius.
-      //
-      // then note
-      // that we now have 48 faces on the
-      // inner and 48 on the outer sphere,
-      // each with an area of approximately
-      // 4*pi/48*r^2 and 4*pi/48, for
-      // a face edge length of approximately
-      // sqrt(pi/12)*r and sqrt(pi/12)
-      //
-      // let's say we put the interior nodes
-      // at a distance rho, then a measure of
-      // deformation for the inner cells
-      // would be
-      //   di=max(sqrt(pi/12)*r/(rho-r),
-      //          (rho-r)/sqrt(pi/12)/r)
-      // and for the outer cells
-      //   do=max(sqrt(pi/12)/(1-rho),
-      //          (1-rho)/sqrt(pi/12))
-      //
-      // we now seek a rho so that the
-      // deformation of cells on the inside
-      // and outside is equal. there are in
-      // principle four possibilities for one
-      // of the branches of do== one of the
-      // branches of di, though not all of
-      // them satisfy do==di, of
-      // course. however, we are not
-      // interested in cases where the inner
-      // cell is long and skinny and the
-      // outer one tall -- yes, they have the
-      // same aspect ratio, but in different
-      // space directions.
-      //
-      // so it only boils down to the
-      // following two possibilities: the
-      // first branch of each max(.,.)
-      // functions are equal, or the second
-      // one are. on the other hand, since
-      // they two branches are reciprocals of
-      // each other, if one pair of branches
-      // is equal, so is the other
-      //
-      // this yields the following equation
-      // for rho:
-      //   sqrt(pi/12)*r/(rho-r)
-      //   == sqrt(pi/12)/(1-rho)
-      // with solution rho=2r/(1+r)
-      const double r = inner_radius / outer_radius;
-      const double rho = 2*r/(1+r);
-
-      // then this is the distance of the
-      // interior nodes from the center:
-      const double middle_radius = rho * outer_radius;
-
-      // mark vertices we've already moved or
-      // that we want to ignore: we don't
-      // want to move vertices at the inner
-      // or outer boundaries
-      std::vector<bool> vertex_already_treated (tmp.n_vertices(), false);
-      for (Triangulation<3>::active_cell_iterator cell = tmp.begin_active();
-           cell != tmp.end(); ++cell)
-        for (unsigned int f=0; f<GeometryInfo<3>::faces_per_cell; ++f)
-          if (cell->at_boundary(f))
-            for (unsigned int v=0; v<GeometryInfo<3>::vertices_per_face; ++v)
-              vertex_already_treated[cell->face(f)->vertex_index(v)] = true;
-
-      // now move the remaining vertices
-      for (Triangulation<3>::active_cell_iterator cell = tmp.begin_active();
-           cell != tmp.end(); ++cell)
-        for (unsigned int v=0; v<GeometryInfo<3>::vertices_per_cell; ++v)
-          if (vertex_already_treated[cell->vertex_index(v)] == false)
-            {
-              // this is a new interior
-              // vertex. mesh refinement may
-              // have placed it at a number
-              // of places in radial
-              // direction and oftentimes not
-              // in a particularly good
-              // one. move it to halfway
-              // between inner and outer
-              // sphere
-              const Point<3> old_distance = cell->vertex(v) - p;
-              const double old_radius = cell->vertex(v).distance(p);
-              cell->vertex(v) = p + old_distance * (middle_radius / old_radius);
-
-              vertex_already_treated[cell->vertex_index(v)] = true;
-            }
+        cells.resize(n_cells, CellData<3>());
 
-      // now copy the resulting level 1 cells
-      // into the new triangulation,
-      cells.resize(tmp.n_active_cells(), CellData<3>());
+        for (unsigned int i=0; i<n_cells; ++i)
+          {
+            for (unsigned int j=0; j<4; ++j)
+              {
+                cells[i].vertices[j  ] = rhombi[i][j];
+                cells[i].vertices[j+4] = rhombi[i][j] + 14;
+              }
+            cells[i].material_id = 0;
+          }
 
-      unsigned int index = 0;
-      for (Triangulation<3>::active_cell_iterator cell = tmp.begin_active();
-           cell != tmp.end(); ++cell, ++index)
-        {
+        tria.create_triangulation (vertices, cells, SubCellData());
+      }
+    else if (n == 96)
+      {
+        // create a triangulation based on the
+        // 12-cell one where we refine the mesh
+        // once and then re-arrange all
+        // interior nodes so that the mesh is
+        // the least distorted
+        HyperShellBoundary<3> boundary (p);
+        Triangulation<3> tmp;
+        hyper_shell (tmp, p, inner_radius, outer_radius, 12);
+        tmp.set_boundary(0, boundary);
+        tmp.set_boundary(1, boundary);
+        tmp.refine_global (1);
+
+        // let's determine the distance at
+        // which the interior nodes should be
+        // from the center. let's say we
+        // measure distances in multiples of
+        // outer_radius and call
+        // r=inner_radius.
+        //
+        // then note
+        // that we now have 48 faces on the
+        // inner and 48 on the outer sphere,
+        // each with an area of approximately
+        // 4*pi/48*r^2 and 4*pi/48, for
+        // a face edge length of approximately
+        // sqrt(pi/12)*r and sqrt(pi/12)
+        //
+        // let's say we put the interior nodes
+        // at a distance rho, then a measure of
+        // deformation for the inner cells
+        // would be
+        //   di=max(sqrt(pi/12)*r/(rho-r),
+        //          (rho-r)/sqrt(pi/12)/r)
+        // and for the outer cells
+        //   do=max(sqrt(pi/12)/(1-rho),
+        //          (1-rho)/sqrt(pi/12))
+        //
+        // we now seek a rho so that the
+        // deformation of cells on the inside
+        // and outside is equal. there are in
+        // principle four possibilities for one
+        // of the branches of do== one of the
+        // branches of di, though not all of
+        // them satisfy do==di, of
+        // course. however, we are not
+        // interested in cases where the inner
+        // cell is long and skinny and the
+        // outer one tall -- yes, they have the
+        // same aspect ratio, but in different
+        // space directions.
+        //
+        // so it only boils down to the
+        // following two possibilities: the
+        // first branch of each max(.,.)
+        // functions are equal, or the second
+        // one are. on the other hand, since
+        // they two branches are reciprocals of
+        // each other, if one pair of branches
+        // is equal, so is the other
+        //
+        // this yields the following equation
+        // for rho:
+        //   sqrt(pi/12)*r/(rho-r)
+        //   == sqrt(pi/12)/(1-rho)
+        // with solution rho=2r/(1+r)
+        const double r = inner_radius / outer_radius;
+        const double rho = 2*r/(1+r);
+
+        // then this is the distance of the
+        // interior nodes from the center:
+        const double middle_radius = rho * outer_radius;
+
+        // mark vertices we've already moved or
+        // that we want to ignore: we don't
+        // want to move vertices at the inner
+        // or outer boundaries
+        std::vector<bool> vertex_already_treated (tmp.n_vertices(), false);
+        for (Triangulation<3>::active_cell_iterator cell = tmp.begin_active();
+             cell != tmp.end(); ++cell)
+          for (unsigned int f=0; f<GeometryInfo<3>::faces_per_cell; ++f)
+            if (cell->at_boundary(f))
+              for (unsigned int v=0; v<GeometryInfo<3>::vertices_per_face; ++v)
+                vertex_already_treated[cell->face(f)->vertex_index(v)] = true;
+
+        // now move the remaining vertices
+        for (Triangulation<3>::active_cell_iterator cell = tmp.begin_active();
+             cell != tmp.end(); ++cell)
           for (unsigned int v=0; v<GeometryInfo<3>::vertices_per_cell; ++v)
-            cells[index].vertices[v] = cell->vertex_index(v);
-          cells[index].material_id = 0;
-        }
-
-      tria.create_triangulation (tmp.get_vertices(), cells, SubCellData());
-    }
-  else
-    {
-      Assert(false, ExcMessage ("Invalid number of coarse mesh cells."));
-    }
+            if (vertex_already_treated[cell->vertex_index(v)] == false)
+              {
+                // this is a new interior
+                // vertex. mesh refinement may
+                // have placed it at a number
+                // of places in radial
+                // direction and oftentimes not
+                // in a particularly good
+                // one. move it to halfway
+                // between inner and outer
+                // sphere
+                const Point<3> old_distance = cell->vertex(v) - p;
+                const double old_radius = cell->vertex(v).distance(p);
+                cell->vertex(v) = p + old_distance * (middle_radius / old_radius);
+
+                vertex_already_treated[cell->vertex_index(v)] = true;
+              }
 
-  if (colorize)
-    colorize_hyper_shell(tria, p, inner_radius, outer_radius);
-}
+        // now copy the resulting level 1 cells
+        // into the new triangulation,
+        cells.resize(tmp.n_active_cells(), CellData<3>());
 
+        unsigned int index = 0;
+        for (Triangulation<3>::active_cell_iterator cell = tmp.begin_active();
+             cell != tmp.end(); ++cell, ++index)
+          {
+            for (unsigned int v=0; v<GeometryInfo<3>::vertices_per_cell; ++v)
+              cells[index].vertices[v] = cell->vertex_index(v);
+            cells[index].material_id = 0;
+          }
 
+        tria.create_triangulation (tmp.get_vertices(), cells, SubCellData());
+      }
+    else
+      {
+        Assert(false, ExcMessage ("Invalid number of coarse mesh cells."));
+      }
 
+    if (colorize)
+      colorize_hyper_shell(tria, p, inner_radius, outer_radius);
+  }
 
-// Implementation for 3D only
-template <>
-void
-GridGenerator::half_hyper_shell (Triangulation<3> &tria,
-                                 const Point<3> &center,
-                                 const double inner_radius,
-                                 const double outer_radius,
-                                 const size_type n,
-                                 const bool colorize)
-{
-  Assert ((inner_radius > 0) && (inner_radius < outer_radius),
-          ExcInvalidRadii ());
-  Assert(colorize == false, ExcNotImplemented());
 
-  if (n <= 5)
-    {
-      // These are for the two lower squares
-      const double d = outer_radius/std::sqrt(2.0);
-      const double a = inner_radius/std::sqrt(2.0);
-      // These are for the two upper square
-      const double b = a/2.0;
-      const double c = d/2.0;
-      // And so are these
-      const double hb = inner_radius*std::sqrt(3.0)/2.0;
-      const double hc = outer_radius*std::sqrt(3.0)/2.0;
-
-      Point<3> vertices[16] =
-      {
-        center+Point<3>( 0,  d, -d),
-        center+Point<3>( 0, -d, -d),
-        center+Point<3>( 0,  a, -a),
-        center+Point<3>( 0, -a, -a),
-        center+Point<3>( 0,  a,  a),
-        center+Point<3>( 0, -a,  a),
-        center+Point<3>( 0,  d,  d),
-        center+Point<3>( 0, -d,  d),
-
-        center+Point<3>(hc,  c, -c),
-        center+Point<3>(hc, -c, -c),
-        center+Point<3>(hb,  b, -b),
-        center+Point<3>(hb, -b, -b),
-        center+Point<3>(hb,  b,  b),
-        center+Point<3>(hb, -b,  b),
-        center+Point<3>(hc,  c,  c),
-        center+Point<3>(hc, -c,  c),
-      };
 
-      int cell_vertices[5][8] =
-      {
-        {0, 1, 8, 9, 2, 3, 10, 11},
-        {0, 2, 8, 10, 6, 4, 14, 12},
-        {1, 7, 9, 15, 3, 5, 11, 13},
-        {6, 4, 14, 12, 7, 5, 15, 13},
-        {8, 10, 9, 11, 14, 12, 15, 13}
-      };
 
-      std::vector<CellData<3> > cells (5, CellData<3>());
+// Implementation for 3D only
+  template <>
+  void
+  half_hyper_shell (Triangulation<3> &tria,
+                    const Point<3> &center,
+                    const double inner_radius,
+                    const double outer_radius,
+                    const unsigned int n,
+                    const bool colorize)
+  {
+    Assert ((inner_radius > 0) && (inner_radius < outer_radius),
+            ExcInvalidRadii ());
+    Assert(colorize == false, ExcNotImplemented());
 
-      for (unsigned int i=0; i<5; ++i)
+    if (n <= 5)
+      {
+        // These are for the two lower squares
+        const double d = outer_radius/std::sqrt(2.0);
+        const double a = inner_radius/std::sqrt(2.0);
+        // These are for the two upper square
+        const double b = a/2.0;
+        const double c = d/2.0;
+        // And so are these
+        const double hb = inner_radius*std::sqrt(3.0)/2.0;
+        const double hc = outer_radius*std::sqrt(3.0)/2.0;
+
+        Point<3> vertices[16] =
         {
-          for (unsigned int j=0; j<8; ++j)
-            cells[i].vertices[j] = cell_vertices[i][j];
-          cells[i].material_id = 0;
+          center+Point<3>( 0,  d, -d),
+          center+Point<3>( 0, -d, -d),
+          center+Point<3>( 0,  a, -a),
+          center+Point<3>( 0, -a, -a),
+          center+Point<3>( 0,  a,  a),
+          center+Point<3>( 0, -a,  a),
+          center+Point<3>( 0,  d,  d),
+          center+Point<3>( 0, -d,  d),
+
+          center+Point<3>(hc,  c, -c),
+          center+Point<3>(hc, -c, -c),
+          center+Point<3>(hb,  b, -b),
+          center+Point<3>(hb, -b, -b),
+          center+Point<3>(hb,  b,  b),
+          center+Point<3>(hb, -b,  b),
+          center+Point<3>(hc,  c,  c),
+          center+Point<3>(hc, -c,  c),
         };
 
-      tria.create_triangulation (
-        std::vector<Point<3> >(&vertices[0], &vertices[16]),
-        cells,
-        SubCellData());       // no boundary information
-    }
-  else
-    {
-      Assert(false, ExcIndexRange(n, 0, 5));
-    }
+        int cell_vertices[5][8] =
+        {
+          {0, 1, 8, 9, 2, 3, 10, 11},
+          {0, 2, 8, 10, 6, 4, 14, 12},
+          {1, 7, 9, 15, 3, 5, 11, 13},
+          {6, 4, 14, 12, 7, 5, 15, 13},
+          {8, 10, 9, 11, 14, 12, 15, 13}
+        };
 
-}
+        std::vector<CellData<3> > cells (5, CellData<3>());
 
-// Implementation for 3D only
-template <>
-void
-GridGenerator::colorize_quarter_hyper_shell(Triangulation<3> &tria,
-                                            const Point<3> &center,
-                                            const double inner_radius,
-                                            const double outer_radius)
-{
+        for (unsigned int i=0; i<5; ++i)
+          {
+            for (unsigned int j=0; j<8; ++j)
+              cells[i].vertices[j] = cell_vertices[i][j];
+            cells[i].material_id = 0;
+          };
+
+        tria.create_triangulation (
+          std::vector<Point<3> >(&vertices[0], &vertices[16]),
+          cells,
+          SubCellData());       // no boundary information
+      }
+    else
+      {
+        Assert(false, ExcIndexRange(n, 0, 5));
+      }
 
-  if (tria.n_cells() != 3)
-    AssertThrow (false, ExcNotImplemented());
+  }
 
-  double middle = (outer_radius-inner_radius)/2e0 + inner_radius;
-  double eps = 1e-3*middle;
-  Triangulation<3>::cell_iterator cell = tria.begin();
 
-  for (; cell!=tria.end(); ++cell)
-    for (unsigned int f=0; f<GeometryInfo<3>::faces_per_cell; ++f)
+// Implementation for 3D only
+  template <>
+  void quarter_hyper_shell (Triangulation<3> &tria,
+                            const Point<3> &center,
+                            const double inner_radius,
+                            const double outer_radius,
+                            const unsigned int n,
+                            const bool colorize)
+  {
+    Assert ((inner_radius > 0) && (inner_radius < outer_radius),
+            ExcInvalidRadii ());
+    if (n == 0 || n == 3)
       {
-        if (!cell->face(f)->at_boundary())
-          continue;
+        const double a = inner_radius*std::sqrt(2.0)/2e0;
+        const double b = outer_radius*std::sqrt(2.0)/2e0;
+        const double c = a*std::sqrt(3.0)/2e0;
+        const double d = b*std::sqrt(3.0)/2e0;
+        const double e = outer_radius/2e0;
+        const double h = inner_radius/2e0;
+
+        std::vector<Point<3> > vertices;
+
+        vertices.push_back (center+Point<3>( 0,  inner_radius, 0)); //0
+        vertices.push_back (center+Point<3>( a,  a, 0));                  //1
+        vertices.push_back (center+Point<3>( b,  b, 0));                  //2
+        vertices.push_back (center+Point<3>( 0, outer_radius, 0));        //3
+        vertices.push_back (center+Point<3>( 0, a , a));                  //4
+        vertices.push_back (center+Point<3>( c, c , h));                  //5
+        vertices.push_back (center+Point<3>( d, d , e));                  //6
+        vertices.push_back (center+Point<3>( 0, b , b));                  //7
+        vertices.push_back (center+Point<3>( inner_radius, 0 , 0));       //8
+        vertices.push_back (center+Point<3>( outer_radius, 0 , 0));       //9
+        vertices.push_back (center+Point<3>( a, 0 , a));  //10
+        vertices.push_back (center+Point<3>( b, 0 , b));  //11
+        vertices.push_back (center+Point<3>( 0, 0 , inner_radius));       //12
+        vertices.push_back (center+Point<3>( 0, 0 , outer_radius));       //13
+
+        const int cell_vertices[3][8] =
+        {
+          {0, 1, 3, 2, 4, 5, 7, 6},
+          {1, 8, 2, 9, 5, 10, 6, 11},
+          {4, 5, 7, 6, 12, 10, 13, 11},
+        };
+        std::vector<CellData<3> > cells(3);
 
-        double radius = cell->face(f)->center().norm() - center.norm();
-        if (std::fabs(cell->face(f)->center()(0)) < eps ) // x = 0 set boundary 2
-          {
-            cell->face(f)->set_boundary_indicator(2);
-            for (unsigned int j=0; j<GeometryInfo<3>::lines_per_face; ++j)
-              if (cell->face(f)->line(j)->at_boundary())
-                if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) > eps)
-                  cell->face(f)->line(j)->set_boundary_indicator(2);
-          }
-        else if (std::fabs(cell->face(f)->center()(1)) < eps) // y = 0 set boundary 3
-          {
-            cell->face(f)->set_boundary_indicator(3);
-            for (unsigned int j=0; j<GeometryInfo<3>::lines_per_face; ++j)
-              if (cell->face(f)->line(j)->at_boundary())
-                if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) > eps)
-                  cell->face(f)->line(j)->set_boundary_indicator(3);
-          }
-        else if (std::fabs(cell->face(f)->center()(2)) < eps ) // z = 0 set boundary 4
-          {
-            cell->face(f)->set_boundary_indicator(4);
-            for (unsigned int j=0; j<GeometryInfo<3>::lines_per_face; ++j)
-              if (cell->face(f)->line(j)->at_boundary())
-                if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) > eps)
-                  cell->face(f)->line(j)->set_boundary_indicator(4);
-          }
-        else if (radius < middle) // inner radius set boundary 0
+        for (unsigned int i=0; i<3; ++i)
           {
-            cell->face(f)->set_boundary_indicator(0);
-            for (unsigned int j=0; j<GeometryInfo<3>::lines_per_face; ++j)
-              if (cell->face(f)->line(j)->at_boundary())
-                if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) < eps)
-                  cell->face(f)->line(j)->set_boundary_indicator(0);
+            for (unsigned int j=0; j<8; ++j)
+              cells[i].vertices[j] = cell_vertices[i][j];
+            cells[i].material_id = 0;
           }
-        else if (radius > middle) // outer radius set boundary 1
-          {
-            cell->face(f)->set_boundary_indicator(1);
-            for (unsigned int j=0; j<GeometryInfo<3>::lines_per_face; ++j)
-              if (cell->face(f)->line(j)->at_boundary())
-                if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) < eps)
-                  cell->face(f)->line(j)->set_boundary_indicator(1);
-          }
-        else
-          AssertThrow (false, ExcInternalError());
+
+        tria.create_triangulation ( vertices, cells, SubCellData());       // no boundary information
+      }
+    else
+      {
+        AssertThrow(false, ExcNotImplemented());
       }
-}
 
+    if (colorize)
+      colorize_quarter_hyper_shell(tria, center, inner_radius, outer_radius);
+  }
 
 
 // Implementation for 3D only
-template <>
-void GridGenerator::quarter_hyper_shell (Triangulation<3> &tria,
-                                         const Point<3> &center,
-                                         const double inner_radius,
-                                         const double outer_radius,
-                                         const size_type n,
-                                         const bool colorize)
-{
-  Assert ((inner_radius > 0) && (inner_radius < outer_radius),
-          ExcInvalidRadii ());
-  if (n == 0 || n == 3)
-    {
-      const double a = inner_radius*std::sqrt(2.0)/2e0;
-      const double b = outer_radius*std::sqrt(2.0)/2e0;
-      const double c = a*std::sqrt(3.0)/2e0;
-      const double d = b*std::sqrt(3.0)/2e0;
-      const double e = outer_radius/2e0;
-      const double h = inner_radius/2e0;
-
-      std::vector<Point<3> > vertices;
-
-      vertices.push_back (center+Point<3>( 0,  inner_radius, 0)); //0
-      vertices.push_back (center+Point<3>( a,  a, 0));                  //1
-      vertices.push_back (center+Point<3>( b,  b, 0));                  //2
-      vertices.push_back (center+Point<3>( 0, outer_radius, 0));        //3
-      vertices.push_back (center+Point<3>( 0, a , a));                  //4
-      vertices.push_back (center+Point<3>( c, c , h));                  //5
-      vertices.push_back (center+Point<3>( d, d , e));                  //6
-      vertices.push_back (center+Point<3>( 0, b , b));                  //7
-      vertices.push_back (center+Point<3>( inner_radius, 0 , 0));       //8
-      vertices.push_back (center+Point<3>( outer_radius, 0 , 0));       //9
-      vertices.push_back (center+Point<3>( a, 0 , a));  //10
-      vertices.push_back (center+Point<3>( b, 0 , b));  //11
-      vertices.push_back (center+Point<3>( 0, 0 , inner_radius));       //12
-      vertices.push_back (center+Point<3>( 0, 0 , outer_radius));       //13
-
-      const int cell_vertices[3][8] =
-      {
-        {0, 1, 3, 2, 4, 5, 7, 6},
-        {1, 8, 2, 9, 5, 10, 6, 11},
-        {4, 5, 7, 6, 12, 10, 13, 11},
+  template <>
+  void cylinder_shell (Triangulation<3> &tria,
+                       const double      length,
+                       const double      inner_radius,
+                       const double      outer_radius,
+                       const unsigned int   n_radial_cells,
+                       const unsigned int   n_axial_cells)
+  {
+    Assert ((inner_radius > 0) && (inner_radius < outer_radius),
+            ExcInvalidRadii ());
+
+    const double pi = numbers::PI;
+
+    // determine the number of cells
+    // for the grid. if not provided by
+    // the user determine it such that
+    // the length of each cell on the
+    // median (in the middle between
+    // the two circles) is equal to its
+    // radial extent (which is the
+    // difference between the two
+    // radii)
+    const unsigned int N_r = (n_radial_cells == 0 ?
+                              static_cast<unsigned int>
+                              (std::ceil((2*pi* (outer_radius + inner_radius)/2) /
+                                         (outer_radius - inner_radius))) :
+                              n_radial_cells);
+    const unsigned int N_z = (n_axial_cells == 0 ?
+                              static_cast<unsigned int>
+                              (std::ceil (length /
+                                          (2*pi*(outer_radius + inner_radius)/2/N_r))) :
+                              n_axial_cells);
+
+    // set up N vertices on the
+    // outer and N vertices on
+    // the inner circle. the
+    // first N ones are on the
+    // outer one, and all are
+    // numbered counter-clockwise
+    std::vector<Point<2> > vertices_2d(2*N_r);
+    for (unsigned int i=0; i<N_r; ++i)
+      {
+        vertices_2d[i] = Point<2>( std::cos(2*pi * i/N_r),
+                                   std::sin(2*pi * i/N_r)) * outer_radius;
+        vertices_2d[i+N_r] = vertices_2d[i] * (inner_radius/outer_radius);
       };
-      std::vector<CellData<3> > cells(3);
 
-      for (unsigned int i=0; i<3; ++i)
+    std::vector<Point<3> > vertices_3d;
+    vertices_3d.reserve (2*N_r*(N_z+1));
+    for (unsigned int j=0; j<=N_z; ++j)
+      for (unsigned int i=0; i<2*N_r; ++i)
         {
-          for (unsigned int j=0; j<8; ++j)
-            cells[i].vertices[j] = cell_vertices[i][j];
-          cells[i].material_id = 0;
+          const Point<3> v (vertices_2d[i][0],
+                            vertices_2d[i][1],
+                            j*length/N_z);
+          vertices_3d.push_back (v);
         }
 
-      tria.create_triangulation ( vertices, cells, SubCellData());       // no boundary information
-    }
-  else
-    {
-      AssertThrow(false, ExcNotImplemented());
-    }
-
-  if (colorize)
-    colorize_quarter_hyper_shell(tria, center, inner_radius, outer_radius);
-}
+    std::vector<CellData<3> > cells (N_r*N_z, CellData<3>());
 
+    for (unsigned int j=0; j<N_z; ++j)
+      for (unsigned int i=0; i<N_r; ++i)
+        {
+          cells[i+j*N_r].vertices[0] = i + (j+1)*2*N_r;
+          cells[i+j*N_r].vertices[1] = (i+1)%N_r + (j+1)*2*N_r;
+          cells[i+j*N_r].vertices[2] = i + j*2*N_r;
+          cells[i+j*N_r].vertices[3] = (i+1)%N_r + j*2*N_r;
 
-// Implementation for 3D only
-template <>
-void GridGenerator::cylinder_shell (Triangulation<3> &tria,
-                                    const double      length,
-                                    const double      inner_radius,
-                                    const double      outer_radius,
-                                    const size_type   n_radial_cells,
-                                    const size_type   n_axial_cells)
-{
-  Assert ((inner_radius > 0) && (inner_radius < outer_radius),
-          ExcInvalidRadii ());
-
-  const double pi = numbers::PI;
-
-  // determine the number of cells
-  // for the grid. if not provided by
-  // the user determine it such that
-  // the length of each cell on the
-  // median (in the middle between
-  // the two circles) is equal to its
-  // radial extent (which is the
-  // difference between the two
-  // radii)
-  const size_type N_r = (n_radial_cells == 0 ?
-                         static_cast<size_type>
-                         (std::ceil((2*pi* (outer_radius + inner_radius)/2) /
-                                    (outer_radius - inner_radius))) :
-                         n_radial_cells);
-  const size_type N_z = (n_axial_cells == 0 ?
-                         static_cast<size_type>
-                         (std::ceil (length /
-                                     (2*pi*(outer_radius + inner_radius)/2/N_r))) :
-                         n_axial_cells);
-
-  // set up N vertices on the
-  // outer and N vertices on
-  // the inner circle. the
-  // first N ones are on the
-  // outer one, and all are
-  // numbered counter-clockwise
-  std::vector<Point<2> > vertices_2d(2*N_r);
-  for (size_type i=0; i<N_r; ++i)
-    {
-      vertices_2d[i] = Point<2>( std::cos(2*pi * i/N_r),
-                                 std::sin(2*pi * i/N_r)) * outer_radius;
-      vertices_2d[i+N_r] = vertices_2d[i] * (inner_radius/outer_radius);
-    };
+          cells[i+j*N_r].vertices[4] = N_r+i + (j+1)*2*N_r;
+          cells[i+j*N_r].vertices[5] = N_r+((i+1)%N_r) + (j+1)*2*N_r;
+          cells[i+j*N_r].vertices[6] = N_r+i + j*2*N_r;
+          cells[i+j*N_r].vertices[7] = N_r+((i+1)%N_r) + j*2*N_r;
 
-  std::vector<Point<3> > vertices_3d;
-  vertices_3d.reserve (2*N_r*(N_z+1));
-  for (size_type j=0; j<=N_z; ++j)
-    for (size_type i=0; i<2*N_r; ++i)
-      {
-        const Point<3> v (vertices_2d[i][0],
-                          vertices_2d[i][1],
-                          j*length/N_z);
-        vertices_3d.push_back (v);
-      }
+          cells[i+j*N_r].material_id = 0;
+        }
 
-  std::vector<CellData<3> > cells (N_r*N_z, CellData<3>());
+    tria.create_triangulation (
+      vertices_3d, cells, SubCellData());
+  }
 
-  for (size_type j=0; j<N_z; ++j)
-    for (size_type i=0; i<N_r; ++i)
-      {
-        cells[i+j*N_r].vertices[0] = i + (j+1)*2*N_r;
-        cells[i+j*N_r].vertices[1] = (i+1)%N_r + (j+1)*2*N_r;
-        cells[i+j*N_r].vertices[2] = i + j*2*N_r;
-        cells[i+j*N_r].vertices[3] = (i+1)%N_r + j*2*N_r;
 
-        cells[i+j*N_r].vertices[4] = N_r+i + (j+1)*2*N_r;
-        cells[i+j*N_r].vertices[5] = N_r+((i+1)%N_r) + (j+1)*2*N_r;
-        cells[i+j*N_r].vertices[6] = N_r+i + j*2*N_r;
-        cells[i+j*N_r].vertices[7] = N_r+((i+1)%N_r) + j*2*N_r;
 
-        cells[i+j*N_r].material_id = 0;
+  template <int dim, int spacedim>
+  void
+  merge_triangulations (const Triangulation<dim, spacedim> &triangulation_1,
+                        const Triangulation<dim, spacedim> &triangulation_2,
+                        Triangulation<dim, spacedim>       &result)
+  {
+    Assert (triangulation_1.n_levels() == 1,
+            ExcMessage ("The input triangulations must be coarse meshes."));
+    Assert (triangulation_2.n_levels() == 1,
+            ExcMessage ("The input triangulations must be coarse meshes."));
+
+    // get the union of the set of vertices
+    std::vector<Point<spacedim> > vertices = triangulation_1.get_vertices();
+    vertices.insert (vertices.end(),
+                     triangulation_2.get_vertices().begin(),
+                     triangulation_2.get_vertices().end());
+
+    // now form the union of the set of cells
+    std::vector<CellData<dim> > cells;
+    cells.reserve (triangulation_1.n_cells() + triangulation_2.n_cells());
+    for (typename Triangulation<dim,spacedim>::cell_iterator
+         cell = triangulation_1.begin(); cell != triangulation_1.end(); ++cell)
+      {
+        CellData<dim> this_cell;
+        for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+          this_cell.vertices[v] = cell->vertex_index(v);
+        this_cell.material_id = cell->material_id();
+        cells.push_back (this_cell);
       }
 
-  tria.create_triangulation (
-    vertices_3d, cells, SubCellData());
-}
-
-
-
-template <int dim, int spacedim>
-void
-GridGenerator::
-merge_triangulations (const Triangulation<dim, spacedim> &triangulation_1,
-                      const Triangulation<dim, spacedim> &triangulation_2,
-                      Triangulation<dim, spacedim>       &result)
-{
-  Assert (triangulation_1.n_levels() == 1,
-          ExcMessage ("The input triangulations must be coarse meshes."));
-  Assert (triangulation_2.n_levels() == 1,
-          ExcMessage ("The input triangulations must be coarse meshes."));
-
-  // get the union of the set of vertices
-  std::vector<Point<spacedim> > vertices = triangulation_1.get_vertices();
-  vertices.insert (vertices.end(),
-                   triangulation_2.get_vertices().begin(),
-                   triangulation_2.get_vertices().end());
-
-  // now form the union of the set of cells
-  std::vector<CellData<dim> > cells;
-  cells.reserve (triangulation_1.n_cells() + triangulation_2.n_cells());
-  for (typename Triangulation<dim,spacedim>::cell_iterator
-       cell = triangulation_1.begin(); cell != triangulation_1.end(); ++cell)
-    {
-      CellData<dim> this_cell;
-      for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
-        this_cell.vertices[v] = cell->vertex_index(v);
-      this_cell.material_id = cell->material_id();
-      cells.push_back (this_cell);
-    }
-
-  // now do the same for the other other mesh. note that we have to
-  // translate the vertex indices
-  for (typename Triangulation<dim,spacedim>::cell_iterator
-       cell = triangulation_2.begin(); cell != triangulation_2.end(); ++cell)
-    {
-      CellData<dim> this_cell;
-      for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
-        this_cell.vertices[v] = cell->vertex_index(v) + triangulation_1.n_vertices();
-      this_cell.material_id = cell->material_id();
-      cells.push_back (this_cell);
-    }
-
-  // throw out duplicated vertices from the two meshes
-  // and create the triangulation
-  SubCellData subcell_data;
-  std::vector<unsigned int> considered_vertices;
-  GridTools::delete_duplicated_vertices (vertices, cells, subcell_data, considered_vertices);
-  result.clear ();
-  result.create_triangulation (vertices, cells, subcell_data);
-}
-
-void
-GridGenerator::
-extrude_triangulation(const Triangulation<2, 2> &input,
-                      const size_type n_slices,
-                      const double height,
-                      Triangulation<3,3> &result)
-{
-  Assert (input.n_levels() == 1,
-          ExcMessage ("The input triangulations must be coarse meshes."));
-  Assert(result.n_cells()==0, ExcMessage("resultin Triangulation need to be empty upon calling extrude_triangulation."));
-  Assert(height>0, ExcMessage("The height in extrude_triangulation needs to be positive."));
-  Assert(n_slices>=2, ExcMessage("The number of slices in extrude_triangulation needs to be at least 2."));
-
-  std::vector<Point<3> > points(n_slices*input.n_vertices());
-  std::vector<CellData<3> > cells;
-  cells.reserve((n_slices-1)*input.n_active_cells());
+    // now do the same for the other other mesh. note that we have to
+    // translate the vertex indices
+    for (typename Triangulation<dim,spacedim>::cell_iterator
+         cell = triangulation_2.begin(); cell != triangulation_2.end(); ++cell)
+      {
+        CellData<dim> this_cell;
+        for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+          this_cell.vertices[v] = cell->vertex_index(v) + triangulation_1.n_vertices();
+        this_cell.material_id = cell->material_id();
+        cells.push_back (this_cell);
+      }
 
-  for (size_type slice=0; slice<n_slices; ++slice)
-    {
-      for (size_type i=0; i<input.n_vertices(); ++i)
+    // throw out duplicated vertices from the two meshes
+    // and create the triangulation
+    SubCellData subcell_data;
+    std::vector<unsigned int> considered_vertices;
+    GridTools::delete_duplicated_vertices (vertices, cells, subcell_data, considered_vertices);
+    result.clear ();
+    result.create_triangulation (vertices, cells, subcell_data);
+  }
+
+
+  void
+  extrude_triangulation(const Triangulation<2, 2> &input,
+                        const unsigned int n_slices,
+                        const double height,
+                        Triangulation<3,3> &result)
+  {
+    Assert (input.n_levels() == 1,
+            ExcMessage ("The input triangulations must be coarse meshes."));
+    Assert(result.n_cells()==0, ExcMessage("resultin Triangulation need to be empty upon calling extrude_triangulation."));
+    Assert(height>0, ExcMessage("The height in extrude_triangulation needs to be positive."));
+    Assert(n_slices>=2, ExcMessage("The number of slices in extrude_triangulation needs to be at least 2."));
 
-        {
-          const Point<2> &v = input.get_vertices()[i];
-          points[i+slice*input.n_vertices()](0) = v(0);
-          points[i+slice*input.n_vertices()](1) = v(1);
-          points[i+slice*input.n_vertices()](2) = height * slice / (n_slices-1);
-        }
-    }
+    std::vector<Point<3> > points(n_slices*input.n_vertices());
+    std::vector<CellData<3> > cells;
+    cells.reserve((n_slices-1)*input.n_active_cells());
 
-  for (Triangulation<2,2>::cell_iterator
-       cell = input.begin(); cell != input.end(); ++cell)
-    {
-      for (size_type slice=0; slice<n_slices-1; ++slice)
-        {
-          CellData<3> this_cell;
-          for (unsigned int v=0; v<GeometryInfo<2>::vertices_per_cell; ++v)
-            {
-              this_cell.vertices[v]
-                = cell->vertex_index(v)+slice*input.n_vertices();
-              this_cell.vertices[v+GeometryInfo<2>::vertices_per_cell]
-                = cell->vertex_index(v)+(slice+1)*input.n_vertices();
-            }
+    for (unsigned int slice=0; slice<n_slices; ++slice)
+      {
+        for (unsigned int i=0; i<input.n_vertices(); ++i)
 
-          this_cell.material_id = cell->material_id();
-          cells.push_back(this_cell);
-        }
-    }
+          {
+            const Point<2> &v = input.get_vertices()[i];
+            points[i+slice*input.n_vertices()](0) = v(0);
+            points[i+slice*input.n_vertices()](1) = v(1);
+            points[i+slice*input.n_vertices()](2) = height * slice / (n_slices-1);
+          }
+      }
 
-  SubCellData s;
-  types::boundary_id bid=0;
-  s.boundary_quads.reserve(input.n_active_lines()*(n_slices-1) + input.n_active_cells()*2);
-  for (Triangulation<2,2>::cell_iterator
-       cell = input.begin(); cell != input.end(); ++cell)
-    {
-      CellData<2> quad;
-      for (unsigned int f=0; f<4; ++f)
-        if (cell->at_boundary(f))
+    for (Triangulation<2,2>::cell_iterator
+         cell = input.begin(); cell != input.end(); ++cell)
+      {
+        for (unsigned int slice=0; slice<n_slices-1; ++slice)
           {
-            quad.boundary_id = cell->face(f)->boundary_indicator();
-            bid = std::max(bid, quad.boundary_id);
-            for (size_type slice=0; slice<n_slices-1; ++slice)
+            CellData<3> this_cell;
+            for (unsigned int v=0; v<GeometryInfo<2>::vertices_per_cell; ++v)
               {
-                quad.vertices[0] = cell->face(f)->vertex_index(0)+slice*input.n_vertices();
-                quad.vertices[1] = cell->face(f)->vertex_index(1)+slice*input.n_vertices();
-                quad.vertices[2] = cell->face(f)->vertex_index(0)+(slice+1)*input.n_vertices();
-                quad.vertices[3] = cell->face(f)->vertex_index(1)+(slice+1)*input.n_vertices();
-                s.boundary_quads.push_back(quad);
+                this_cell.vertices[v]
+                  = cell->vertex_index(v)+slice*input.n_vertices();
+                this_cell.vertices[v+GeometryInfo<2>::vertices_per_cell]
+                  = cell->vertex_index(v)+(slice+1)*input.n_vertices();
               }
-          }
-    }
 
-  for (Triangulation<2,2>::cell_iterator
-       cell = input.begin(); cell != input.end(); ++cell)
-    {
-      CellData<2> quad;
-      quad.boundary_id = bid + 1;
-      quad.vertices[0] = cell->vertex_index(0);
-      quad.vertices[1] = cell->vertex_index(1);
-      quad.vertices[2] = cell->vertex_index(2);
-      quad.vertices[3] = cell->vertex_index(3);
-      s.boundary_quads.push_back(quad);
-
-      quad.boundary_id = bid + 2;
-      for (int i=0; i<4; ++i)
-        quad.vertices[i] += (n_slices-1)*input.n_vertices();
-      s.boundary_quads.push_back(quad);
-    }
+            this_cell.material_id = cell->material_id();
+            cells.push_back(this_cell);
+          }
+      }
 
+    SubCellData s;
+    types::boundary_id bid=0;
+    s.boundary_quads.reserve(input.n_active_lines()*(n_slices-1) + input.n_active_cells()*2);
+    for (Triangulation<2,2>::cell_iterator
+         cell = input.begin(); cell != input.end(); ++cell)
+      {
+        CellData<2> quad;
+        for (unsigned int f=0; f<4; ++f)
+          if (cell->at_boundary(f))
+            {
+              quad.boundary_id = cell->face(f)->boundary_indicator();
+              bid = std::max(bid, quad.boundary_id);
+              for (unsigned int slice=0; slice<n_slices-1; ++slice)
+                {
+                  quad.vertices[0] = cell->face(f)->vertex_index(0)+slice*input.n_vertices();
+                  quad.vertices[1] = cell->face(f)->vertex_index(1)+slice*input.n_vertices();
+                  quad.vertices[2] = cell->face(f)->vertex_index(0)+(slice+1)*input.n_vertices();
+                  quad.vertices[3] = cell->face(f)->vertex_index(1)+(slice+1)*input.n_vertices();
+                  s.boundary_quads.push_back(quad);
+                }
+            }
+      }
 
+    for (Triangulation<2,2>::cell_iterator
+         cell = input.begin(); cell != input.end(); ++cell)
+      {
+        CellData<2> quad;
+        quad.boundary_id = bid + 1;
+        quad.vertices[0] = cell->vertex_index(0);
+        quad.vertices[1] = cell->vertex_index(1);
+        quad.vertices[2] = cell->vertex_index(2);
+        quad.vertices[3] = cell->vertex_index(3);
+        s.boundary_quads.push_back(quad);
+
+        quad.boundary_id = bid + 2;
+        for (int i=0; i<4; ++i)
+          quad.vertices[i] += (n_slices-1)*input.n_vertices();
+        s.boundary_quads.push_back(quad);
+      }
 
-  result.create_triangulation (
-    points,
-    cells,
-    s);
-}
+    result.create_triangulation (points,
+                                 cells,
+                                 s);
+  }
 
 
-// make the following function inline. this is so that it becomes marked
-// internal/weak for the linker and we don't get multiply defined errors
-// when linking with more than one dimension at a time. Usually we used
-// the trick of putting these functions in a .all_dimensions.cc file, but
-// this is not necessary here as this is an internal only function.
-inline
-void GridGenerator::laplace_solve (const SparseMatrix<double> &S,
-                                   const std::map<size_type,double> &m,
-                                   Vector<double> &u)
-{
-  const size_type n_dofs=S.n();
-  FilteredMatrix<Vector<double> > SF (S);
-  PreconditionJacobi<SparseMatrix<double> > prec;
-  prec.initialize(S, 1.2);
-  FilteredMatrix<Vector<double> > PF (prec);
+  /**
+   * Solve the Laplace equation for @p laplace_transformation function for one
+   * of the @p dim space dimensions. Factorized into a function of its own
+   * in order to allow parallel execution.
+   */
+  void laplace_solve (const SparseMatrix<double> &S,
+                      const std::map<unsigned int,double> &m,
+                      Vector<double> &u)
+  {
+    const unsigned int n_dofs=S.n();
+    FilteredMatrix<Vector<double> > SF (S);
+    PreconditionJacobi<SparseMatrix<double> > prec;
+    prec.initialize(S, 1.2);
+    FilteredMatrix<Vector<double> > PF (prec);
 
-  SolverControl control (n_dofs, 1.e-10, false, false);
-  GrowingVectorMemory<Vector<double> > mem;
-  SolverCG<Vector<double> > solver (control, mem);
+    SolverControl control (n_dofs, 1.e-10, false, false);
+    GrowingVectorMemory<Vector<double> > mem;
+    SolverCG<Vector<double> > solver (control, mem);
 
-  Vector<double> f(n_dofs);
+    Vector<double> f(n_dofs);
 
-  SF.add_constraints(m);
-  SF.apply_constraints (f, true);
-  solver.solve(SF, u, f, PF);
-}
+    SF.add_constraints(m);
+    SF.apply_constraints (f, true);
+    solver.solve(SF, u, f, PF);
+  }
 
 
 // Implementation for 1D only
-template <>
-void GridGenerator::laplace_transformation (Triangulation<1> &,
-                                            const std::map<size_type,Point<1> > &)
-{
-  Assert(false, ExcNotImplemented());
-}
+  template <>
+  void laplace_transformation (Triangulation<1> &,
+                               const std::map<unsigned int,Point<1> > &)
+  {
+    Assert(false, ExcNotImplemented());
+  }
 
 
 // Implementation for dimensions except 1
-template <int dim>
-void GridGenerator::laplace_transformation (Triangulation<dim> &tria,
-                                            const std::map<size_type,Point<dim> > &new_points)
-{
-  // first provide everything that is
-  // needed for solving a Laplace
-  // equation.
-  MappingQ1<dim> mapping_q1;
-  FE_Q<dim> q1(1);
-
-  DoFHandler<dim> dof_handler(tria);
-  dof_handler.distribute_dofs(q1);
-
-  CompressedSparsityPattern c_sparsity_pattern (dof_handler.n_dofs (),
-                                                dof_handler.n_dofs ());
-  DoFTools::make_sparsity_pattern (dof_handler, c_sparsity_pattern);
-  c_sparsity_pattern.compress ();
-
-  SparsityPattern sparsity_pattern;
-  sparsity_pattern.copy_from (c_sparsity_pattern);
-  sparsity_pattern.compress ();
-
-  SparseMatrix<double> S(sparsity_pattern);
-
-  QGauss<dim> quadrature(4);
-
-  MatrixCreator::create_laplace_matrix(mapping_q1, dof_handler, quadrature, S);
-
-  // set up the boundary values for
-  // the laplace problem
-  std::vector<std::map<size_type,double> > m(dim);
-  typename std::map<size_type,Point<dim> >::const_iterator map_iter;
-  typename std::map<size_type,Point<dim> >::const_iterator map_end=new_points.end();
-
-  // fill these maps using the data
-  // given by new_points
-  typename DoFHandler<dim>::cell_iterator cell=dof_handler.begin_active(),
-                                          endc=dof_handler.end();
-  typename DoFHandler<dim>::face_iterator face;
-  for (; cell!=endc; ++cell)
-    {
-      if (cell->at_boundary())
-        for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
-          {
-            face=cell->face(face_no);
-            if (face->at_boundary())
-              for (unsigned int vertex_no=0;
-                   vertex_no<GeometryInfo<dim>::vertices_per_face; ++vertex_no)
-                {
-                  const size_type vertex_index=face->vertex_index(vertex_no);
-                  map_iter=new_points.find(vertex_index);
-
-                  if (map_iter!=map_end)
-                    for (unsigned int i=0; i<dim; ++i)
-                      m[i].insert(std::pair<size_type,double> (
-                                    face->vertex_dof_index(vertex_no, 0), map_iter->second(i)));
-                }
-          }
-    }
+  template <int dim>
+  void laplace_transformation (Triangulation<dim> &tria,
+                               const std::map<unsigned int,Point<dim> > &new_points)
+  {
+    // first provide everything that is
+    // needed for solving a Laplace
+    // equation.
+    MappingQ1<dim> mapping_q1;
+    FE_Q<dim> q1(1);
+
+    DoFHandler<dim> dof_handler(tria);
+    dof_handler.distribute_dofs(q1);
+
+    CompressedSparsityPattern c_sparsity_pattern (dof_handler.n_dofs (),
+                                                  dof_handler.n_dofs ());
+    DoFTools::make_sparsity_pattern (dof_handler, c_sparsity_pattern);
+    c_sparsity_pattern.compress ();
+
+    SparsityPattern sparsity_pattern;
+    sparsity_pattern.copy_from (c_sparsity_pattern);
+    sparsity_pattern.compress ();
+
+    SparseMatrix<double> S(sparsity_pattern);
+
+    QGauss<dim> quadrature(4);
+
+    MatrixCreator::create_laplace_matrix(mapping_q1, dof_handler, quadrature, S);
+
+    // set up the boundary values for
+    // the laplace problem
+    std::vector<std::map<unsigned int,double> > m(dim);
+    typename std::map<unsigned int,Point<dim> >::const_iterator map_iter;
+    typename std::map<unsigned int,Point<dim> >::const_iterator map_end=new_points.end();
+
+    // fill these maps using the data
+    // given by new_points
+    typename DoFHandler<dim>::cell_iterator cell=dof_handler.begin_active(),
+                                            endc=dof_handler.end();
+    typename DoFHandler<dim>::face_iterator face;
+    for (; cell!=endc; ++cell)
+      {
+        if (cell->at_boundary())
+          for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+            {
+              face=cell->face(face_no);
+              if (face->at_boundary())
+                for (unsigned int vertex_no=0;
+                     vertex_no<GeometryInfo<dim>::vertices_per_face; ++vertex_no)
+                  {
+                    const unsigned int vertex_index=face->vertex_index(vertex_no);
+                    map_iter=new_points.find(vertex_index);
 
-  // solve the dim problems with
-  // different right hand sides.
-  Vector<double> us[dim];
-  for (unsigned int i=0; i<dim; ++i)
-    us[i].reinit (dof_handler.n_dofs());
-
-  // solve linear systems in parallel
-  Threads::TaskGroup<> tasks;
-  for (unsigned int i=0; i<dim; ++i)
-    tasks += Threads::new_task (&GridGenerator::laplace_solve,
-                                S, m[i], us[i]);
-  tasks.join_all ();
-
-  // change the coordinates of the
-  // points of the triangulation
-  // according to the computed values
-  for (cell=dof_handler.begin_active(); cell!=endc; ++cell)
-    for (size_type vertex_no=0;
-         vertex_no<GeometryInfo<dim>::vertices_per_cell; ++vertex_no)
-      {
-        Point<dim> &v=cell->vertex(vertex_no);
-        const size_type dof_index=cell->vertex_dof_index(vertex_no, 0);
-        for (unsigned int i=0; i<dim; ++i)
-          v(i)=us[i](dof_index);
+                    if (map_iter!=map_end)
+                      for (unsigned int i=0; i<dim; ++i)
+                        m[i].insert(std::pair<unsigned int,double> (
+                                      face->vertex_dof_index(vertex_no, 0), map_iter->second(i)));
+                  }
+            }
       }
-}
 
+    // solve the dim problems with
+    // different right hand sides.
+    Vector<double> us[dim];
+    for (unsigned int i=0; i<dim; ++i)
+      us[i].reinit (dof_handler.n_dofs());
+
+    // solve linear systems in parallel
+    Threads::TaskGroup<> tasks;
+    for (unsigned int i=0; i<dim; ++i)
+      tasks += Threads::new_task (&laplace_solve,
+                                  S, m[i], us[i]);
+    tasks.join_all ();
+
+    // change the coordinates of the
+    // points of the triangulation
+    // according to the computed values
+    for (cell=dof_handler.begin_active(); cell!=endc; ++cell)
+      for (unsigned int vertex_no=0;
+           vertex_no<GeometryInfo<dim>::vertices_per_cell; ++vertex_no)
+        {
+          Point<dim> &v=cell->vertex(vertex_no);
+          const unsigned int dof_index=cell->vertex_dof_index(vertex_no, 0);
+          for (unsigned int i=0; i<dim; ++i)
+            v(i)=us[i](dof_index);
+        }
+  }
 
 
-template <>
-void GridGenerator::hyper_cube_with_cylindrical_hole (Triangulation<1> &,
-                                                      const double,
-                                                      const double,
-                                                      const double,
-                                                      const size_type,
-                                                      bool)
-{
-  Assert(false, ExcNotImplemented());
-}
 
+  template <>
+  void hyper_cube_with_cylindrical_hole (Triangulation<1> &,
+                                         const double,
+                                         const double,
+                                         const double,
+                                         const unsigned int,
+                                         bool)
+  {
+    Assert(false, ExcNotImplemented());
+  }
 
 
-template <>
-void
-GridGenerator::hyper_cube_with_cylindrical_hole (Triangulation<2> &triangulation,
-                                                 const double inner_radius,
-                                                 const double outer_radius,
-                                                 const double, // width,
-                                                 const size_type, // width_repetition,
-                                                 bool colorize)
-{
-  const int dim = 2;
-
-  Assert(inner_radius < outer_radius,
-         ExcMessage("outer_radius has to be bigger than inner_radius."));
-
-  Point<dim> center;
-  // We create an hyper_shell in two dimensions, and then we modify it.
-  GridGenerator::hyper_shell (triangulation,
-                              center, inner_radius, outer_radius,
-                              8);
-  Triangulation<dim>::active_cell_iterator
-  cell = triangulation.begin_active(),
-  endc = triangulation.end();
-  std::vector<bool> treated_vertices(triangulation.n_vertices(), false);
-  for (; cell != endc; ++cell)
-    {
-      for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
-        if (cell->face(f)->at_boundary())
-          {
-            for (unsigned int v=0; v < GeometryInfo<dim>::vertices_per_face; ++v)
-              {
-                size_type vv = cell->face(f)->vertex_index(v);
-                if (treated_vertices[vv] == false)
-                  {
-                    treated_vertices[vv] = true;
-                    switch (vv)
-                      {
-                      case 1:
-                        cell->face(f)->vertex(v) = center+Point<dim>(outer_radius,outer_radius);
-                        break;
-                      case 3:
-                        cell->face(f)->vertex(v) = center+Point<dim>(-outer_radius,outer_radius);
-                        break;
-                      case 5:
-                        cell->face(f)->vertex(v) = center+Point<dim>(-outer_radius,-outer_radius);
-                        break;
-                      case 7:
-                        cell->face(f)->vertex(v) = center+Point<dim>(outer_radius,-outer_radius);
-                      default:
-                        break;
-                      }
-                  }
-              }
-          }
-    }
-  double eps = 1e-3 * outer_radius;
-  cell = triangulation.begin_active();
-  for (; cell != endc; ++cell)
-    {
-      for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
-        if (cell->face(f)->at_boundary())
-          {
-            double dx = cell->face(f)->center()(0) - center(0);
-            double dy = cell->face(f)->center()(1) - center(1);
-            if (colorize)
-              {
-                if (std::abs(dx + outer_radius) < eps)
-                  cell->face(f)->set_boundary_indicator(0);
-                else if (std::abs(dx - outer_radius) < eps)
-                  cell->face(f)->set_boundary_indicator(1);
-                else if (std::abs(dy + outer_radius) < eps)
-                  cell->face(f)->set_boundary_indicator(2);
-                else if (std::abs(dy - outer_radius) < eps)
-                  cell->face(f)->set_boundary_indicator(3);
-                else
-                  cell->face(f)->set_boundary_indicator(4);
-              }
-            else
-              {
-                double d = (cell->face(f)->center() - center).norm();
-                if (d-inner_radius < 0)
-                  cell->face(f)->set_boundary_indicator(1);
-                else
-                  cell->face(f)->set_boundary_indicator(0);
-              }
-          }
-    }
-}
 
-template <>
-void GridGenerator::hyper_cube_with_cylindrical_hole(Triangulation<3> &triangulation,
-                                                     const double inner_radius,
-                                                     const double outer_radius,
-                                                     const double L,
-                                                     const size_type Nz,
-                                                     bool colorize)
-{
-  const int dim = 3;
-
-  Assert(inner_radius < outer_radius,
-         ExcMessage("outer_radius has to be bigger than inner_radius."));
-  Assert(L > 0,
-         ExcMessage("Must give positive extension L"));
-  Assert(Nz >= 1, ExcLowerRange(1, Nz));
-
-  GridGenerator::cylinder_shell (triangulation,
-                                 L, inner_radius, outer_radius,
-                                 8,
-                                 Nz);
-
-  Triangulation<dim>::active_cell_iterator
-  cell = triangulation.begin_active(),
-  endc = triangulation.end();
-  std::vector<bool> treated_vertices(triangulation.n_vertices(), false);
-  for (; cell != endc; ++cell)
-    {
-      for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
-        if (cell->face(f)->at_boundary())
-          {
-            for (unsigned int v=0; v < GeometryInfo<dim>::vertices_per_face; ++v)
-              {
-                size_type vv = cell->face(f)->vertex_index(v);
-                if (treated_vertices[vv] == false)
-                  {
-                    treated_vertices[vv] = true;
-                    for (size_type i=0; i<=Nz; ++i)
-                      {
-                        double d = ((double) i)*L/((double) Nz);
-                        switch (vv-i*16)
-                          {
-                          case 1:
-                            cell->face(f)->vertex(v) = Point<dim>(outer_radius,outer_radius,d);
-                            break;
-                          case 3:
-                            cell->face(f)->vertex(v) = Point<dim>(-outer_radius,outer_radius,d);
-                            break;
-                          case 5:
-                            cell->face(f)->vertex(v) = Point<dim>(-outer_radius,-outer_radius,d);
-                            break;
-                          case 7:
-                            cell->face(f)->vertex(v) = Point<dim>(outer_radius,-outer_radius,d);
-                            break;
-                          default:
-                            break;
-                          }
-                      }
-                  }
-              }
-          }
-    }
-  double eps = 1e-3 * outer_radius;
-  cell = triangulation.begin_active();
-  for (; cell != endc; ++cell)
-    {
-      for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
-        if (cell->face(f)->at_boundary())
-          {
-            double dx = cell->face(f)->center()(0);
-            double dy = cell->face(f)->center()(1);
-            double dz = cell->face(f)->center()(2);
+  template <>
+  void
+  hyper_cube_with_cylindrical_hole (Triangulation<2> &triangulation,
+                                    const double inner_radius,
+                                    const double outer_radius,
+                                    const double, // width,
+                                    const unsigned int, // width_repetition,
+                                    bool colorize)
+  {
+    const int dim = 2;
+
+    Assert(inner_radius < outer_radius,
+           ExcMessage("outer_radius has to be bigger than inner_radius."));
+
+    Point<dim> center;
+    // We create an hyper_shell in two dimensions, and then we modify it.
+    hyper_shell (triangulation,
+                 center, inner_radius, outer_radius,
+                 8);
+    Triangulation<dim>::active_cell_iterator
+    cell = triangulation.begin_active(),
+    endc = triangulation.end();
+    std::vector<bool> treated_vertices(triangulation.n_vertices(), false);
+    for (; cell != endc; ++cell)
+      {
+        for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+          if (cell->face(f)->at_boundary())
+            {
+              for (unsigned int v=0; v < GeometryInfo<dim>::vertices_per_face; ++v)
+                {
+                  unsigned int vv = cell->face(f)->vertex_index(v);
+                  if (treated_vertices[vv] == false)
+                    {
+                      treated_vertices[vv] = true;
+                      switch (vv)
+                        {
+                        case 1:
+                          cell->face(f)->vertex(v) = center+Point<dim>(outer_radius,outer_radius);
+                          break;
+                        case 3:
+                          cell->face(f)->vertex(v) = center+Point<dim>(-outer_radius,outer_radius);
+                          break;
+                        case 5:
+                          cell->face(f)->vertex(v) = center+Point<dim>(-outer_radius,-outer_radius);
+                          break;
+                        case 7:
+                          cell->face(f)->vertex(v) = center+Point<dim>(outer_radius,-outer_radius);
+                        default:
+                          break;
+                        }
+                    }
+                }
+            }
+      }
+    double eps = 1e-3 * outer_radius;
+    cell = triangulation.begin_active();
+    for (; cell != endc; ++cell)
+      {
+        for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+          if (cell->face(f)->at_boundary())
+            {
+              double dx = cell->face(f)->center()(0) - center(0);
+              double dy = cell->face(f)->center()(1) - center(1);
+              if (colorize)
+                {
+                  if (std::abs(dx + outer_radius) < eps)
+                    cell->face(f)->set_boundary_indicator(0);
+                  else if (std::abs(dx - outer_radius) < eps)
+                    cell->face(f)->set_boundary_indicator(1);
+                  else if (std::abs(dy + outer_radius) < eps)
+                    cell->face(f)->set_boundary_indicator(2);
+                  else if (std::abs(dy - outer_radius) < eps)
+                    cell->face(f)->set_boundary_indicator(3);
+                  else
+                    cell->face(f)->set_boundary_indicator(4);
+                }
+              else
+                {
+                  double d = (cell->face(f)->center() - center).norm();
+                  if (d-inner_radius < 0)
+                    cell->face(f)->set_boundary_indicator(1);
+                  else
+                    cell->face(f)->set_boundary_indicator(0);
+                }
+            }
+      }
+  }
 
-            if (colorize)
-              {
-                if (std::abs(dx + outer_radius) < eps)
-                  cell->face(f)->set_boundary_indicator(0);
 
-                else if (std::abs(dx - outer_radius) < eps)
-                  cell->face(f)->set_boundary_indicator(1);
 
-                else if (std::abs(dy + outer_radius) < eps)
-                  cell->face(f)->set_boundary_indicator(2);
+  template <>
+  void hyper_cube_with_cylindrical_hole(Triangulation<3> &triangulation,
+                                        const double inner_radius,
+                                        const double outer_radius,
+                                        const double L,
+                                        const unsigned int Nz,
+                                        bool colorize)
+  {
+    const int dim = 3;
+
+    Assert(inner_radius < outer_radius,
+           ExcMessage("outer_radius has to be bigger than inner_radius."));
+    Assert(L > 0,
+           ExcMessage("Must give positive extension L"));
+    Assert(Nz >= 1, ExcLowerRange(1, Nz));
+
+    cylinder_shell (triangulation,
+                    L, inner_radius, outer_radius,
+                    8,
+                    Nz);
+
+    Triangulation<dim>::active_cell_iterator
+    cell = triangulation.begin_active(),
+    endc = triangulation.end();
+    std::vector<bool> treated_vertices(triangulation.n_vertices(), false);
+    for (; cell != endc; ++cell)
+      {
+        for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+          if (cell->face(f)->at_boundary())
+            {
+              for (unsigned int v=0; v < GeometryInfo<dim>::vertices_per_face; ++v)
+                {
+                  unsigned int vv = cell->face(f)->vertex_index(v);
+                  if (treated_vertices[vv] == false)
+                    {
+                      treated_vertices[vv] = true;
+                      for (unsigned int i=0; i<=Nz; ++i)
+                        {
+                          double d = ((double) i)*L/((double) Nz);
+                          switch (vv-i*16)
+                            {
+                            case 1:
+                              cell->face(f)->vertex(v) = Point<dim>(outer_radius,outer_radius,d);
+                              break;
+                            case 3:
+                              cell->face(f)->vertex(v) = Point<dim>(-outer_radius,outer_radius,d);
+                              break;
+                            case 5:
+                              cell->face(f)->vertex(v) = Point<dim>(-outer_radius,-outer_radius,d);
+                              break;
+                            case 7:
+                              cell->face(f)->vertex(v) = Point<dim>(outer_radius,-outer_radius,d);
+                              break;
+                            default:
+                              break;
+                            }
+                        }
+                    }
+                }
+            }
+      }
+    double eps = 1e-3 * outer_radius;
+    cell = triangulation.begin_active();
+    for (; cell != endc; ++cell)
+      {
+        for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+          if (cell->face(f)->at_boundary())
+            {
+              double dx = cell->face(f)->center()(0);
+              double dy = cell->face(f)->center()(1);
+              double dz = cell->face(f)->center()(2);
 
-                else if (std::abs(dy - outer_radius) < eps)
-                  cell->face(f)->set_boundary_indicator(3);
+              if (colorize)
+                {
+                  if (std::abs(dx + outer_radius) < eps)
+                    cell->face(f)->set_boundary_indicator(0);
 
-                else if (std::abs(dz) < eps)
-                  cell->face(f)->set_boundary_indicator(4);
+                  else if (std::abs(dx - outer_radius) < eps)
+                    cell->face(f)->set_boundary_indicator(1);
 
-                else if (std::abs(dz - L) < eps)
-                  cell->face(f)->set_boundary_indicator(5);
+                  else if (std::abs(dy + outer_radius) < eps)
+                    cell->face(f)->set_boundary_indicator(2);
 
-                else
-                  {
-                    cell->face(f)->set_boundary_indicator(6);
-                    for (unsigned int l=0; l<GeometryInfo<dim>::lines_per_face; ++l)
-                      cell->face(f)->line(l)->set_boundary_indicator(6);
-                  }
+                  else if (std::abs(dy - outer_radius) < eps)
+                    cell->face(f)->set_boundary_indicator(3);
 
-              }
-            else
-              {
-                Point<dim> c = cell->face(f)->center();
-                c(2) = 0;
-                double d = c.norm();
-                if (d-inner_radius < 0)
-                  {
-                    cell->face(f)->set_boundary_indicator(1);
-                    for (unsigned int l=0; l<GeometryInfo<dim>::lines_per_face; ++l)
-                      cell->face(f)->line(l)->set_boundary_indicator(1);
-                  }
-                else
-                  cell->face(f)->set_boundary_indicator(0);
-              }
-          }
-    }
-}
+                  else if (std::abs(dz) < eps)
+                    cell->face(f)->set_boundary_indicator(4);
+
+                  else if (std::abs(dz - L) < eps)
+                    cell->face(f)->set_boundary_indicator(5);
+
+                  else
+                    {
+                      cell->face(f)->set_boundary_indicator(6);
+                      for (unsigned int l=0; l<GeometryInfo<dim>::lines_per_face; ++l)
+                        cell->face(f)->line(l)->set_boundary_indicator(6);
+                    }
 
+                }
+              else
+                {
+                  Point<dim> c = cell->face(f)->center();
+                  c(2) = 0;
+                  double d = c.norm();
+                  if (d-inner_radius < 0)
+                    {
+                      cell->face(f)->set_boundary_indicator(1);
+                      for (unsigned int l=0; l<GeometryInfo<dim>::lines_per_face; ++l)
+                        cell->face(f)->line(l)->set_boundary_indicator(1);
+                    }
+                  else
+                    cell->face(f)->set_boundary_indicator(0);
+                }
+            }
+      }
+  }
+}
 
 // explicit instantiations
 #include "grid_generator.inst"
index fadcea0bc79024596630f0fb46e782b97f07c0ff..39a7641c44f2ae312d006f5e6d2c6d7829e61cdd 100644 (file)
 
 for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
   {
+namespace GridGenerator 
+\{
 #if deal_II_dimension <= deal_II_space_dimension
     template void
-      GridGenerator::hyper_rectangle<deal_II_dimension, deal_II_space_dimension> (
+      hyper_rectangle<deal_II_dimension, deal_II_space_dimension> (
        Triangulation<deal_II_dimension, deal_II_space_dimension> &,
        const Point<deal_II_space_dimension>&, const Point<deal_II_space_dimension>&,
        const bool);
     
     template void
-      GridGenerator::colorize_hyper_rectangle<deal_II_dimension, deal_II_space_dimension> (
+      colorize_hyper_rectangle<deal_II_dimension, deal_II_space_dimension> (
        Triangulation<deal_II_dimension, deal_II_space_dimension> &);
     
     template void
-      GridGenerator::hyper_cube<deal_II_dimension, deal_II_space_dimension> (
+      hyper_cube<deal_II_dimension, deal_II_space_dimension> (
        Triangulation<deal_II_dimension, deal_II_space_dimension> &, const double, const double);
     
     template
       void
-      GridGenerator::merge_triangulations
+      merge_triangulations
       (const Triangulation<deal_II_dimension,deal_II_space_dimension> &triangulation_1,
        const Triangulation<deal_II_dimension,deal_II_space_dimension> &triangulation_2,
        Triangulation<deal_II_dimension,deal_II_space_dimension>       &result);
       
 #endif
+\}  
   }
 
 
 for (deal_II_dimension : DIMENSIONS)
 {
-  
+namespace GridGenerator \{  
   template void
-    GridGenerator::subdivided_hyper_cube<deal_II_dimension> (
+    subdivided_hyper_cube<deal_II_dimension> (
       Triangulation<deal_II_dimension> &,
       const unsigned int, const double, const double);
   
   template void
-    GridGenerator::subdivided_hyper_rectangle<deal_II_dimension>
+    subdivided_hyper_rectangle<deal_II_dimension>
     (Triangulation<deal_II_dimension> &,
      const std::vector<unsigned int>&,
      const Point<deal_II_dimension>&,
@@ -61,7 +64,7 @@ for (deal_II_dimension : DIMENSIONS)
   
   template
     void
-    GridGenerator::subdivided_hyper_rectangle(
+    subdivided_hyper_rectangle(
       Triangulation<deal_II_dimension>              &,
       const std::vector<std::vector<double> > &,
       const Point<deal_II_dimension>                &,
@@ -69,32 +72,31 @@ for (deal_II_dimension : DIMENSIONS)
       const bool                       );
   
     template void
-      GridGenerator::parallelepiped<deal_II_dimension> (
+      parallelepiped<deal_II_dimension> (
        Triangulation<deal_II_dimension>&,
        const Point<deal_II_dimension> (&) [deal_II_dimension], 
        const bool);  
   
     template void
-      GridGenerator::subdivided_parallelepiped<deal_II_dimension> (
+      subdivided_parallelepiped<deal_II_dimension> (
        Triangulation<deal_II_dimension>&,
-       const size_type,
+       const unsigned int,
        const Point<deal_II_dimension> (&) [deal_II_dimension], 
        const bool);  
 
     template void
-      GridGenerator::subdivided_parallelepiped<deal_II_dimension> (
+      subdivided_parallelepiped<deal_II_dimension> (
        Triangulation<deal_II_dimension>&,
-       const size_type [deal_II_dimension], 
+       const unsigned int [deal_II_dimension], 
        const Point<deal_II_dimension> (&) [deal_II_dimension], 
        const bool);  
   
 #if deal_II_dimension > 1
   template void
-    GridGenerator::
     laplace_transformation<deal_II_dimension> (Triangulation<deal_II_dimension> &,
                                               const std::map<types::global_dof_index,Point<deal_II_dimension> > &);
 #endif
   
-  
+\}  
  }
 

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.