From: Timo Heister <timo.heister@gmail.com>
Date: Wed, 2 Dec 2015 19:53:47 +0000 (-0500)
Subject: fix *_hyper_rectangle for codim and add new subdivided_parallelepiped
X-Git-Tag: v8.4.0-rc2~38^2~1
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=05bfbabc03edcdc9d2e4b73510359676998f3f04;p=dealii.git

fix *_hyper_rectangle for codim and add new subdivided_parallelepiped

- extend documentation what happens if dim!=spacedim
- add new subdivided_parallelepiped for codim1 and 2
- switch to Point<dim> in *_hyper_rectangle
---

diff --git a/include/deal.II/grid/grid_generator.h b/include/deal.II/grid/grid_generator.h
index 41f76f8877..f6d17c586b 100644
--- a/include/deal.II/grid/grid_generator.h
+++ b/include/deal.II/grid/grid_generator.h
@@ -18,6 +18,7 @@
 
 
 #include <deal.II/base/config.h>
+#include <deal.II/base/std_cxx11/array.h>
 #include <deal.II/base/exceptions.h>
 #include <deal.II/base/point.h>
 #include <deal.II/base/table.h>
@@ -48,8 +49,8 @@ namespace GridGenerator
 {
   /**
    * @name Creating meshes for basic geometries
-   * @{
    */
+  ///@{
 
   /**
    * Initialize the given triangulation with a hypercube (line in 1D, square
@@ -64,6 +65,11 @@ namespace GridGenerator
    *
    * @image html hyper_cubes.png
    *
+   * If @p dim < @p spacedim, this will create a @p dim dimensional object
+   * in the first @p dim coordinate directions embedded into the @p spacedim
+   * dimensional space with the remaining entries set to zero. For example,
+   * a <tt>Triangulation@<2,3@></tt> will be a square in the xy plane with z=0.
+   *
    * See also subdivided_hyper_cube() for a coarse mesh consisting of several
    * cells. See hyper_rectangle(), if different lengths in different ordinate
    * directions are required.
@@ -77,8 +83,7 @@ namespace GridGenerator
                    const bool                    colorize= false);
 
   /**
-   * \brief %Triangulation of a d-simplex with (d+1) vertices and mesh cells,
-   * resp.
+   * \brief %Triangulation of a d-simplex with (d+1) vertices and mesh cells.
    *
    * The @p vertices argument contains a vector with all d+1 vertices of the
    * simplex. They must be given in an order such that the vectors from the
@@ -90,6 +95,13 @@ namespace GridGenerator
    * @image html simplex_2d.png
    * @image html simplex_3d.png
    *
+   * @param tria The Triangulation to create. It needs to be empty upon calling this
+   * function.
+   *
+   * @param vertices The dim+1 corners of the simplex.
+   *
+   * @note Implemented for <tt>Triangulation@<2,2@></tt>, <tt>Triangulation@<3,3@></tt>.
+   *
    * @author Guido Kanschat
    * @date 2015
    */
@@ -103,10 +115,10 @@ namespace GridGenerator
    * cells. Thus, the number of cells filling the given volume is
    * <tt>repetitions<sup>dim</sup></tt>.
    *
-   * If spacedim=dim+1 the same mesh as in the case spacedim=dim is created,
-   * but the vertices have an additional coordinate =0. So, if dim=1 one
-   * obtains line along the x axis in the xy plane, and if dim=3 one obtains a
-   * square in lying in the xy plane in 3d space.
+   * If @p dim < @p spacedim, this will create a @p dim dimensional object
+   * in the first @p dim coordinate directions embedded into the @p spacedim
+   * dimensional space with the remaining entries set to zero. For example,
+   * a <tt>Triangulation@<2,3@></tt> will be a square in the xy plane with z=0.
    *
    * @note The triangulation needs to be void upon calling this function.
    */
@@ -129,23 +141,23 @@ namespace GridGenerator
    * 2<sup>i</sup>. For instance, the center point (1,-1,1) yields a material
    * id 5.
    *
-   * @note If spacedim>dim the same mesh as in the case spacedim=dim is
-   * created, but the vertices have all additional coordinates equal to the
-   * coordinates of p1, i.e., we generate a hyper_rectangle parallel to the
-   * x-axis (1d), or to the xy-axis passing through p1. The additional
-   * coordinates of the point p2 are ignored by this class in this case.
+   * If @p dim < @p spacedim, this will create a @p dim dimensional object
+   * in the first @p dim coordinate directions embedded into the @p spacedim
+   * dimensional space with the remaining entries set to zero. For example,
+   * a <tt>Triangulation@<2,3@></tt> will be a rectangle in the xy plane with z=0,
+   * defined by the two opposing corners @p p1 and @p p2.
    *
    * @note The triangulation needs to be void upon calling this function.
    */
   template <int dim, int spacedim>
   void hyper_rectangle (Triangulation<dim,spacedim> &tria,
-                        const Point<spacedim>       &p1,
-                        const Point<spacedim>       &p2,
+                        const Point<dim>            &p1,
+                        const Point<dim>            &p2,
                         const bool                  colorize = false);
 
   /**
    * Create a coordinate-parallel parallelepiped from the two diagonally
-   * opposite corner points @p p1 and @p p2. In dimension @p i,
+   * opposite corner points @p p1 and @p p2. In direction @p i,
    * <tt>repetitions[i]</tt> cells are generated.
    *
    * To get cells with an aspect ratio different from that of the domain, use
@@ -169,17 +181,34 @@ namespace GridGenerator
    * be true. That means the boundary indicator is 0 on the left and 1 on the
    * right.  See step-15 for details.
    *
-   * @note The triangulation needs to be void upon calling this function.
+   * If @p dim < @p spacedim, this will create a @p dim dimensional object
+   * in the first @p dim coordinate directions embedded into the @p spacedim
+   * dimensional space with the remaining entries set to zero. For example,
+   * a <tt>Triangulation@<2,3@></tt> will be a rectangle in the xy plane with z=0,
+   * defined by the two opposing corners @p p1 and @p p2.
    *
    * @note For an example of the use of this function see the step-28 tutorial
    * program.
+   *
+   * @param tria The Triangulation to create. It needs to be empty upon calling this
+   * function.
+   *
+   * @param repetitions A vector of dim positive values denoting the number of cells
+   * to generate in that direction.
+   *
+   * @param p1 First corner point.
+   *
+   * @param p2 Second corner opposite to @p p1.
+   *
+   * @param colorize Assign different boundary ids if set to true.
+   *
    */
   template <int dim, int spacedim>
   void
   subdivided_hyper_rectangle (Triangulation<dim,spacedim>     &tria,
                               const std::vector<unsigned int> &repetitions,
-                              const Point<spacedim>           &p1,
-                              const Point<spacedim>           &p2,
+                              const Point<dim>                &p1,
+                              const Point<dim>                &p2,
                               const bool                      colorize=false);
 
   /**
@@ -238,6 +267,15 @@ namespace GridGenerator
    *
    * @image html cheese_2d.png
    *
+   * If @p dim < @p spacedim, this will create a @p dim dimensional object
+   * in the first @p dim coordinate directions embedded into the @p spacedim
+   * dimensional space with the remaining entries set to zero.
+   *
+   * @param tria The Triangulation to create. It needs to be empty upon calling this
+   * function.
+   *
+   * @param holes Positive number of holes in each of the dim directions.
+
    * @author Guido Kanschat
    * @date 2015
    */
@@ -317,6 +355,32 @@ namespace GridGenerator
                              const Point<dim>   (&corners) [dim],
                              const bool           colorize = false);
 
+  /**
+   * A subdivided parallelepiped.
+   *
+   * @param tria The Triangulation to create. It needs to be empty upon calling this
+   * function.
+   *
+   * @param origin First corner of the parallelepiped.
+   *
+   * @param edges An array of @p dim tensors describing the length and direction of the edges
+   * from @p origin.
+   *
+   * @param subdivisions Number of subdivisions in each of the dim directions. Each
+   * entry must be positive. An empty vector is equivalent to one subdivision in
+   * each direction.
+   *
+   * @param colorize Assign different boundary ids if set to true.
+   *
+   * @note Implemented for all combinations of @p dim and @p spacedim.
+   */
+  template <int dim, int spacedim>
+  void
+  subdivided_parallelepiped (Triangulation<dim, spacedim>  &tria,
+                             const Point<spacedim> &origin,
+                             const std_cxx11::array<Tensor<1,spacedim>,dim> &edges,
+                             const std::vector<unsigned int> &subdivisions = std::vector<unsigned int>(),
+                             const bool colorize = false);
 
   /**
    * Hypercube with a layer of hypercubes around it. The first two parameters
@@ -351,6 +415,10 @@ namespace GridGenerator
    * This function is declared to exist for triangulations of all space
    * dimensions, but throws an error if called in 1d.
    *
+   * You should attach a SphericalManifold to the cells and faces for correct
+   * placement of vertices upon refinement and to be able to use higher order
+   * mappings.
+   *
    * @note The triangulation needs to be void upon calling this function.
    */
   template <int dim>
@@ -359,9 +427,14 @@ namespace GridGenerator
                    const double        radius = 1.);
 
   /**
-    * Creates an hyper sphere.
-    * This function is declared to exist for dim=1, spacedim =2, and
-    * dim=2, spacedim =3.
+    * Creates a hyper sphere, i.e., a surface of a ball in @p spacedim
+    * dimensions.
+    * This function only exists for dim+1=spacedim in 2 and 3 space
+    * dimensions.
+    *
+    * You should attach a SphericalManifold to the cells and faces for correct
+    * placement of vertices upon refinement and to be able to use higher order
+    * mappings.
     *
     * The following pictures are generated with:
     * @code
@@ -376,8 +449,8 @@ namespace GridGenerator
     * triangulation.refine_global(3);
     * @endcode
     *
-    * See the @ref manifold "documentation module on manifolds", for
-    * more detail.
+    * See the @ref manifold "documentation module on manifolds" for
+    * more details.
     *
     * @image html sphere.png
     * @image html sphere_section.png
@@ -753,14 +826,12 @@ namespace GridGenerator
                 const double        R,
                 const double        r);
 
-  /*
-   * @}
-   */
+  ///@}
 
   /**
    * @name Creating meshes from other meshes
-   * @{
    */
+  ///@{
 
   /**
    * Given the two triangulations specified as the first two arguments, create
@@ -918,16 +989,12 @@ namespace GridGenerator
   void flatten_triangulation(const Triangulation<dim,spacedim1> &in_tria,
                              Triangulation<dim,spacedim2> &out_tria);
 
-  /*
-   * @}
-   */
+  ///@}
 
   /**
-   * @name Creating lower-dimensional meshes from parts of higher-dimensional
-   * meshes
-   * @{
+   * @name Creating lower-dimensional meshes from parts of higher-dimensional meshes
    */
-
+  ///@{
 
 #ifdef _MSC_VER
   // Microsoft's VC++ has a bug where it doesn't want to recognize that
@@ -1028,14 +1095,14 @@ namespace GridGenerator
                              const std::set<types::boundary_id> &boundary_ids
                              = std::set<types::boundary_id>());
 
-  /*
-   * @}
-   */
+  ///@}
+
 
   /**
    * @name Exceptions
-   * @{
    */
+  ///@{
+
 
   /**
    * Exception
@@ -1056,9 +1123,7 @@ namespace GridGenerator
                   << "The vector of repetitions  must have "
                   << arg1 <<" elements.");
 
-  /*
-   * @}
-   */
+  ///@}
 }
 
 
diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc
index c8c1443d96..c73d78ae54 100644
--- a/source/grid/grid_generator.cc
+++ b/source/grid/grid_generator.cc
@@ -337,18 +337,19 @@ namespace GridGenerator
   template <int dim, int spacedim>
   void
   hyper_rectangle (Triangulation<dim,spacedim> &tria,
-                   const Point<spacedim>   &p_1,
-                   const Point<spacedim>   &p_2,
+                   const Point<dim>   &p_1,
+                   const Point<dim>   &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));
+    // First, extend dimensions from dim to spacedim and
+    // normalize such that p1 is lower in all coordinate
+    // directions. Additional entries will be 0.
+    Point<spacedim> p1, p2;
+    for (unsigned int i=0; i<dim; ++i)
+      {
+        p1(i) = std::min(p_1(i), p_2(i));
+        p2(i) = std::max(p_1(i), p_2(i));
+      }
 
     std::vector<Point<spacedim> > vertices (GeometryInfo<dim>::vertices_per_cell);
     switch (dim)
@@ -406,12 +407,7 @@ namespace GridGenerator
     Assert (left < right,
             ExcMessage ("Invalid left-to-right bounds of hypercube"));
 
-    Point<spacedim> p1;
-    Point<spacedim> p2;
-
-    p1(spacedim-1) = 0;
-    p2(spacedim-1) = 0;
-
+    Point<dim> p1, p2;
     for (unsigned int i=0; i<dim; ++i)
       {
         p1(i) = left;
@@ -746,22 +742,12 @@ namespace GridGenerator
                  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);
+    Point<2> origin;
+    std_cxx11::array<Tensor<1,2>,2> edges;
+    edges[0] = corners[0];
+    edges[1] = corners[1];
+    std::vector<unsigned int> subdivisions;
+    subdivided_parallelepiped(tria, origin, edges, subdivisions, colorize);
   }
 
 
@@ -801,14 +787,6 @@ namespace GridGenerator
                                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
   subdivided_parallelepiped (Triangulation<dim>  &tria,
@@ -820,41 +798,83 @@ namespace GridGenerator
                              const Point<dim>   (&corners) [dim],
                              const bool           colorize)
   {
-    // Zero n_subdivisions is the origin only, which makes no sense
+    Point<dim> origin;
+    std::vector<unsigned int> subdivisions;
+    std_cxx11::array<Tensor<1,dim>,dim> edges;
+    for (unsigned int i=0; i<dim; ++i)
+      {
+        subdivisions.push_back(n_subdivisions[i]);
+        edges[i] = corners[i];
+      }
+
+    subdivided_parallelepiped (tria, origin, edges, subdivisions, 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, int spacedim>
+  void
+  subdivided_parallelepiped (Triangulation<dim, spacedim>  &tria,
+                             const Point<spacedim> &origin,
+                             const std_cxx11::array<Tensor<1,spacedim>,dim> &edges,
+                             const std::vector<unsigned int> &subdivisions,
+                             const bool colorize)
+  {
+    if (subdivisions.size()==0)
+      {
+        std::vector<unsigned int> new_subdivisions(dim, 1);
+        subdivided_parallelepiped(tria, origin, edges, new_subdivisions, colorize);
+        return;
+      }
+
+    Assert(subdivisions.size()==dim, ExcMessage(""));
+
+    // check subdivisions
     for (unsigned int i=0; i<dim; ++i)
-      Assert (n_subdivisions[i]>0, ExcInvalidRepetitions(n_subdivisions[i]));
+      {
+        Assert (subdivisions[i]>0, ExcInvalidRepetitions(subdivisions[i]));
+        Assert (edges[i].norm()>0, ExcMessage("Edges in subdivided_parallelepiped() must not be degenerate."));
+      }
 
     // 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."));
+        Assert ((edges[i]!=edges[j]),
+                ExcMessage ("Degenerate edges of subdivided_parallelepiped encountered."));
 
     // 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;
+    std::vector<Point<spacedim> > points;
 
     switch (dim)
       {
       case 1:
-        for (unsigned int x=0; x<=n_subdivisions[0]; ++x)
-          points.push_back (double(x)*delta[0]);
+        for (unsigned int x=0; x<=subdivisions[0]; ++x)
+          points.push_back (origin + edges[0]/subdivisions[0]*x);
         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 (double(x)*delta[0] + double(y)*delta[1]);
+        for (unsigned int y=0; y<=subdivisions[1]; ++y)
+          for (unsigned int x=0; x<=subdivisions[0]; ++x)
+            points.push_back (origin
+                              + edges[0]/subdivisions[0]*x
+                              + edges[1]/subdivisions[1]*y);
         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 (double(x)*delta[0] + double(y)*delta[1] + double(z)*delta[2]);
+        for (unsigned int z=0; z<=subdivisions[2]; ++z)
+          for (unsigned int y=0; y<=subdivisions[1]; ++y)
+            for (unsigned int x=0; x<=subdivisions[0]; ++x)
+              points.push_back (
+                origin
+                + edges[0]/subdivisions[0]*x
+                + edges[1]/subdivisions[1]*y
+                + edges[2]/subdivisions[2]*z);
         break;
 
       default:
@@ -864,14 +884,14 @@ namespace GridGenerator
     // Prepare cell data
     unsigned int n_cells = 1;
     for (unsigned int i=0; i<dim; ++i)
-      n_cells *= n_subdivisions[i];
+      n_cells *= 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)
+        for (unsigned int x=0; x<subdivisions[0]; ++x)
           {
             cells[x].vertices[0] = x;
             cells[x].vertices[1] = x+1;
@@ -884,13 +904,13 @@ namespace GridGenerator
       case 2:
       {
         // Shorthand
-        const unsigned int n_dy = n_subdivisions[1];
-        const unsigned int n_dx = n_subdivisions[0];
+        const unsigned int n_dy = subdivisions[1];
+        const unsigned int n_dx = subdivisions[0];
 
         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;
+              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;
@@ -905,15 +925,15 @@ namespace GridGenerator
       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];
+        const unsigned int n_dz = subdivisions[2];
+        const unsigned int n_dy = subdivisions[1];
+        const unsigned int n_dx = 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;
+                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;
@@ -969,7 +989,7 @@ namespace GridGenerator
     Assert (left < right,
             ExcMessage ("Invalid left-to-right bounds of hypercube"));
 
-    Point<spacedim> p0,p1;
+    Point<dim> p0, p1;
     for (unsigned int i=0; i<dim; ++i)
       {
         p0[i] = left;
@@ -987,39 +1007,35 @@ namespace GridGenerator
   subdivided_hyper_rectangle (
     Triangulation<dim, spacedim>              &tria,
     const std::vector<unsigned int> &repetitions,
-    const Point<spacedim>                &p_1,
-    const Point<spacedim>                &p_2,
+    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<spacedim> p1(p_1);
-    Point<spacedim> p2(p_2);
 
+    // First, extend dimensions from dim to spacedim and
+    // normalize such that p1 is lower in all coordinate
+    // directions. Additional entries will be 0.
+    Point<spacedim> p1, p2;
     for (unsigned int i=0; i<dim; ++i)
-      if (p1(i) > p2(i))
-        std::swap (p1(i), p2(i));
+      {
+        p1(i) = std::min(p_1(i), p_2(i));
+        p2(i) = std::max(p_1(i), p_2(i));
+      }
 
-    // then check that all repetitions
-    // are >= 1, and calculate deltas
-    // convert repetitions from double
-    // to int by taking the ceiling.
+    // calculate deltas and validate input
     std::vector<Point<spacedim> > delta(dim);
-
     for (unsigned int i=0; i<dim; ++i)
       {
         Assert (repetitions[i] >= 1, ExcInvalidRepetitions(repetitions[i]));
 
         delta[i][i] = (p2[i]-p1[i])/repetitions[i];
+        Assert(delta[i][i]>0.0,
+               ExcMessage("The first dim entries of coordinates of p1 and p2 need to be different."));
       }
 
-    // then generate the necessary
-    // points
+    // then generate the points
     std::vector<Point<spacedim> > points;
     switch (dim)
       {
@@ -1048,7 +1064,6 @@ namespace GridGenerator
       }
 
     // next create the cells
-    // Prepare cell data
     std::vector<CellData<dim> > cells;
     switch (dim)
       {
@@ -1145,18 +1160,12 @@ namespace GridGenerator
     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
+    // 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
@@ -1367,7 +1376,6 @@ namespace GridGenerator
     const Table<1,types::material_id>                 &material_id,
     const bool                                    colorize)
   {
-    // contributed by Yaqi Wang 2006
     Assert(spacing.size() == 1,
            ExcInvalidRepetitionsDimension(1));
 
@@ -1428,7 +1436,6 @@ namespace GridGenerator
     const Table<2,types::material_id>          &material_id,
     const bool                                    colorize)
   {
-    // contributed by Yaqi Wang 2006
     Assert(spacing.size() == 2,
            ExcInvalidRepetitionsDimension(2));
 
@@ -1526,7 +1533,6 @@ namespace GridGenerator
     const Table<3,types::material_id>               &material_id,
     const bool                                    colorize)
   {
-    // contributed by Yaqi Wang 2006
     const unsigned int dim = 3;
 
     Assert(spacing.size() == dim,
diff --git a/source/grid/grid_generator.inst.in b/source/grid/grid_generator.inst.in
index c82e5a5f12..be673bd3c4 100644
--- a/source/grid/grid_generator.inst.in
+++ b/source/grid/grid_generator.inst.in
@@ -22,7 +22,7 @@ namespace GridGenerator
     template void
       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 Point<deal_II_dimension>&, const Point<deal_II_dimension>&,
 	const bool);
     
     template void
@@ -39,10 +39,18 @@ namespace GridGenerator
     subdivided_hyper_rectangle<deal_II_dimension, deal_II_space_dimension>
     (Triangulation<deal_II_dimension, deal_II_space_dimension> &,
      const std::vector<unsigned int>&,
-     const Point<deal_II_space_dimension>&,
-     const Point<deal_II_space_dimension>&,
+     const Point<deal_II_dimension>&,
+     const Point<deal_II_dimension>&,
      const bool);
 
+      template void
+      subdivided_parallelepiped<deal_II_dimension, deal_II_space_dimension>
+      (Triangulation<deal_II_dimension, deal_II_space_dimension> &,
+          const Point<deal_II_space_dimension> &,
+          const std_cxx11::array<Tensor<1,deal_II_space_dimension>,deal_II_dimension> &,
+                                 const std::vector<unsigned int> &,
+                                 const bool colorize);
+
   template void
     hyper_cross<deal_II_dimension, deal_II_space_dimension>
     (Triangulation<deal_II_dimension, deal_II_space_dimension> &,
@@ -147,6 +155,7 @@ namespace GridGenerator \{
 	const Point<deal_II_dimension> (&) [deal_II_dimension], 
 	const bool);
   #endif
+
   
 \}  
  }
diff --git a/tests/grid/grid_generator_08.cc b/tests/grid/grid_generator_08.cc
index 367f50dc29..4317a6ff8f 100644
--- a/tests/grid/grid_generator_08.cc
+++ b/tests/grid/grid_generator_08.cc
@@ -31,18 +31,15 @@
 template<int dim, int spacedim>
 void test(std::ostream &out)
 {
-  Point<spacedim> p1;
+  Point<dim> p1;
   p1[0] = 2.;
   if (dim>1) p1[1] = -1.;
-  if (dim>2) p1[2] = 0.;
-  Point<spacedim> p2;
+  Point<dim> p2;
   p2[0] = 3.;
   if (dim>1) p2[1] = 2.;
-  if (dim>2) p2[2] = 4.;
-  Point<spacedim> p3;
+  Point<dim> p3;
   p3[0] = 2.;
   if (dim>1) p3[1] = 1.;
-  if (dim>2) p3[2] = 4.;
 
   GridOut go;
   GridOut::OutputFormat format = GridOut::msh;
@@ -85,7 +82,6 @@ void test(std::ostream &out)
     }
 }
 
-
 int main()
 {
   std::ofstream logfile("output");
diff --git a/tests/grid/grid_parallelepiped_05.cc b/tests/grid/grid_parallelepiped_05.cc
new file mode 100644
index 0000000000..473236dee9
--- /dev/null
+++ b/tests/grid/grid_parallelepiped_05.cc
@@ -0,0 +1,116 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2013 - 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// test subdivided_parallelepiped
+
+#include "../tests.h"
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/logstream.h>
+#include <cmath>
+#include <cstdlib>
+
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+
+template <int dim>
+Point<dim> point(double x=0, double y=0, double z=0)
+{
+  Point<dim> p;
+  if (dim>0) p[0] = x;
+  if (dim>1) p[1] = y;
+  if (dim>2) p[2] = z;
+  return p;
+}
+
+
+// The simplest test case is to create a parallelepiped grid with a
+// number of subdivisions and output the result.
+template<int dim, int spacedim>
+void check (bool subdivide)
+{
+  deallog << "dim " << dim << " spacedim " << spacedim << std::endl;
+
+  Point<spacedim> origin;
+  for (unsigned int d=0; d<spacedim; ++d)
+    origin[d] = 0.1+d*1.0;
+
+  std_cxx11::array<Tensor<1,spacedim>,dim> edges;
+  switch (dim)
+    {
+    case 1:
+      edges[0] = point<spacedim>(0.5, 0.02, 0.03);
+      break;
+
+    case 2:
+      edges[0] = point<spacedim>(0.15, 0.50, 0.03);
+      edges[1] = point<spacedim>(0.70, 0.25, -0.01);
+      break;
+
+    case 3:
+      edges[0] = point<spacedim>(0.10, 0.50, 0.70);
+      edges[1] = point<spacedim>(1.50, 0.25, 0.70);
+      edges[2] = point<spacedim>(0.25, 0.50, 0.3);
+      break;
+
+    default:
+      Assert (false, ExcInternalError ());
+    }
+
+
+  std::vector<unsigned int> subdivisions;
+  if (subdivide)
+    {
+      subdivisions.resize(dim);
+      for (unsigned int d=0; d<dim; ++d)
+        subdivisions[d]=d+1;
+    }
+
+
+  bool colorize = false;
+
+  Triangulation<dim,spacedim> triangulation;
+  GridGenerator::subdivided_parallelepiped (triangulation, origin, edges,
+                                            subdivisions, colorize);
+
+  GridOut grid_out;
+
+  grid_out.write_gnuplot (triangulation, deallog.get_file_stream());
+}
+
+int main ()
+{
+  initlog();
+
+  check<1,1> (true);
+  check<1,2> (true);
+  check<1,3> (true);
+
+  check<2,2> (true);
+  check<2,3> (true);
+
+  check<3,3> (true);
+  check<3,3> (false);
+
+
+  deallog<< "OK" << std::endl;
+}
diff --git a/tests/grid/grid_parallelepiped_05.output b/tests/grid/grid_parallelepiped_05.output
new file mode 100644
index 0000000000..753252c3d3
--- /dev/null
+++ b/tests/grid/grid_parallelepiped_05.output
@@ -0,0 +1,217 @@
+
+DEAL::dim 1 spacedim 1
+0.100000 0 0
+0.600000 0 0
+
+
+DEAL::dim 1 spacedim 2
+0.100000 1.10000 0 0
+0.600000 1.12000 0 0
+
+
+DEAL::dim 1 spacedim 3
+0.100000 1.10000 2.10000 0 0
+0.600000 1.12000 2.13000 0 0
+
+
+DEAL::dim 2 spacedim 2
+0.100000 1.10000 0 0
+0.250000 1.60000 0 0
+0.600000 1.72500 0 0
+0.450000 1.22500 0 0
+0.100000 1.10000 0 0
+
+
+0.450000 1.22500 0 0
+0.600000 1.72500 0 0
+0.950000 1.85000 0 0
+0.800000 1.35000 0 0
+0.450000 1.22500 0 0
+
+
+DEAL::dim 2 spacedim 3
+0.100000 1.10000 2.10000 0 0
+0.250000 1.60000 2.13000 0 0
+0.600000 1.72500 2.12500 0 0
+0.450000 1.22500 2.09500 0 0
+0.100000 1.10000 2.10000 0 0
+
+
+0.450000 1.22500 2.09500 0 0
+0.600000 1.72500 2.12500 0 0
+0.950000 1.85000 2.12000 0 0
+0.800000 1.35000 2.09000 0 0
+0.450000 1.22500 2.09500 0 0
+
+
+DEAL::dim 3 spacedim 3
+0.100000 1.10000 2.10000 0 0
+0.200000 1.60000 2.80000 0 0
+0.283333 1.76667 2.90000 0 0
+0.183333 1.26667 2.20000 0 0
+0.100000 1.10000 2.10000 0 0
+
+0.850000 1.22500 2.45000 0 0
+0.950000 1.72500 3.15000 0 0
+1.03333 1.89167 3.25000 0 0
+0.933333 1.39167 2.55000 0 0
+0.850000 1.22500 2.45000 0 0
+
+0.100000 1.10000 2.10000 0 0
+0.850000 1.22500 2.45000 0 0
+
+0.200000 1.60000 2.80000 0 0
+0.950000 1.72500 3.15000 0 0
+
+0.283333 1.76667 2.90000 0 0
+1.03333 1.89167 3.25000 0 0
+
+0.183333 1.26667 2.20000 0 0
+0.933333 1.39167 2.55000 0 0
+
+0.850000 1.22500 2.45000 0 0
+0.950000 1.72500 3.15000 0 0
+1.03333 1.89167 3.25000 0 0
+0.933333 1.39167 2.55000 0 0
+0.850000 1.22500 2.45000 0 0
+
+1.60000 1.35000 2.80000 0 0
+1.70000 1.85000 3.50000 0 0
+1.78333 2.01667 3.60000 0 0
+1.68333 1.51667 2.90000 0 0
+1.60000 1.35000 2.80000 0 0
+
+0.850000 1.22500 2.45000 0 0
+1.60000 1.35000 2.80000 0 0
+
+0.950000 1.72500 3.15000 0 0
+1.70000 1.85000 3.50000 0 0
+
+1.03333 1.89167 3.25000 0 0
+1.78333 2.01667 3.60000 0 0
+
+0.933333 1.39167 2.55000 0 0
+1.68333 1.51667 2.90000 0 0
+
+0.183333 1.26667 2.20000 0 0
+0.283333 1.76667 2.90000 0 0
+0.366667 1.93333 3.00000 0 0
+0.266667 1.43333 2.30000 0 0
+0.183333 1.26667 2.20000 0 0
+
+0.933333 1.39167 2.55000 0 0
+1.03333 1.89167 3.25000 0 0
+1.11667 2.05833 3.35000 0 0
+1.01667 1.55833 2.65000 0 0
+0.933333 1.39167 2.55000 0 0
+
+0.183333 1.26667 2.20000 0 0
+0.933333 1.39167 2.55000 0 0
+
+0.283333 1.76667 2.90000 0 0
+1.03333 1.89167 3.25000 0 0
+
+0.366667 1.93333 3.00000 0 0
+1.11667 2.05833 3.35000 0 0
+
+0.266667 1.43333 2.30000 0 0
+1.01667 1.55833 2.65000 0 0
+
+0.933333 1.39167 2.55000 0 0
+1.03333 1.89167 3.25000 0 0
+1.11667 2.05833 3.35000 0 0
+1.01667 1.55833 2.65000 0 0
+0.933333 1.39167 2.55000 0 0
+
+1.68333 1.51667 2.90000 0 0
+1.78333 2.01667 3.60000 0 0
+1.86667 2.18333 3.70000 0 0
+1.76667 1.68333 3.00000 0 0
+1.68333 1.51667 2.90000 0 0
+
+0.933333 1.39167 2.55000 0 0
+1.68333 1.51667 2.90000 0 0
+
+1.03333 1.89167 3.25000 0 0
+1.78333 2.01667 3.60000 0 0
+
+1.11667 2.05833 3.35000 0 0
+1.86667 2.18333 3.70000 0 0
+
+1.01667 1.55833 2.65000 0 0
+1.76667 1.68333 3.00000 0 0
+
+0.266667 1.43333 2.30000 0 0
+0.366667 1.93333 3.00000 0 0
+0.450000 2.10000 3.10000 0 0
+0.350000 1.60000 2.40000 0 0
+0.266667 1.43333 2.30000 0 0
+
+1.01667 1.55833 2.65000 0 0
+1.11667 2.05833 3.35000 0 0
+1.20000 2.22500 3.45000 0 0
+1.10000 1.72500 2.75000 0 0
+1.01667 1.55833 2.65000 0 0
+
+0.266667 1.43333 2.30000 0 0
+1.01667 1.55833 2.65000 0 0
+
+0.366667 1.93333 3.00000 0 0
+1.11667 2.05833 3.35000 0 0
+
+0.450000 2.10000 3.10000 0 0
+1.20000 2.22500 3.45000 0 0
+
+0.350000 1.60000 2.40000 0 0
+1.10000 1.72500 2.75000 0 0
+
+1.01667 1.55833 2.65000 0 0
+1.11667 2.05833 3.35000 0 0
+1.20000 2.22500 3.45000 0 0
+1.10000 1.72500 2.75000 0 0
+1.01667 1.55833 2.65000 0 0
+
+1.76667 1.68333 3.00000 0 0
+1.86667 2.18333 3.70000 0 0
+1.95000 2.35000 3.80000 0 0
+1.85000 1.85000 3.10000 0 0
+1.76667 1.68333 3.00000 0 0
+
+1.01667 1.55833 2.65000 0 0
+1.76667 1.68333 3.00000 0 0
+
+1.11667 2.05833 3.35000 0 0
+1.86667 2.18333 3.70000 0 0
+
+1.20000 2.22500 3.45000 0 0
+1.95000 2.35000 3.80000 0 0
+
+1.10000 1.72500 2.75000 0 0
+1.85000 1.85000 3.10000 0 0
+
+DEAL::dim 3 spacedim 3
+0.100000 1.10000 2.10000 0 0
+0.200000 1.60000 2.80000 0 0
+0.450000 2.10000 3.10000 0 0
+0.350000 1.60000 2.40000 0 0
+0.100000 1.10000 2.10000 0 0
+
+1.60000 1.35000 2.80000 0 0
+1.70000 1.85000 3.50000 0 0
+1.95000 2.35000 3.80000 0 0
+1.85000 1.85000 3.10000 0 0
+1.60000 1.35000 2.80000 0 0
+
+0.100000 1.10000 2.10000 0 0
+1.60000 1.35000 2.80000 0 0
+
+0.200000 1.60000 2.80000 0 0
+1.70000 1.85000 3.50000 0 0
+
+0.450000 2.10000 3.10000 0 0
+1.95000 2.35000 3.80000 0 0
+
+0.350000 1.60000 2.40000 0 0
+1.85000 1.85000 3.10000 0 0
+
+DEAL::OK