From: young <young@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Sun, 13 Jan 2013 00:08:43 +0000 (+0000)
Subject: Generate a parallelepiped grid in 1d, 2d, and 3d.
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=1ea8bd41ff1e7d20b543d14aab9d00846e67568a;p=dealii-svn.git

Generate a parallelepiped grid in 1d, 2d, and 3d.

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

diff --git a/deal.II/include/deal.II/grid/grid_generator.h b/deal.II/include/deal.II/grid/grid_generator.h
index d9a15f9b5f..9f7b344ed2 100644
--- a/deal.II/include/deal.II/grid/grid_generator.h
+++ b/deal.II/include/deal.II/grid/grid_generator.h
@@ -303,11 +303,33 @@ public:
    * function.
    */
   template <int dim>
-  static void
-  parallelogram(Triangulation<dim> &tria,
-                const Tensor<2,dim> &corners,
-                const bool           colorize=false);
+  static 
+    void
+    parallelogram(Triangulation<dim> &tria,
+		  const Tensor<2,dim> &corners,
+		  const bool           colorize=false);
 
+  /**
+   * A parallelepiped. The first corner point is the origin. The
+   * <tt>dim</tt> adjacent points are the rank one subtensors of the
+   * tensor provided and additional points will be sums of those
+   * <tt>dim</tt> vectors. Thus in 1d this is a line, in 2d this is a
+   * parallelogram, and in 3d a parallelepiped. Colorizing is done
+   * according to hyper_rectangle().
+   *
+   * @note This function silently reorders the vertices on the cells
+   * to lexiographic ordering if they are not already ordered that way
+   * (see GridReordering::reorder_grid()).
+   *
+   * @note The triangulation needs to be void upon calling this
+   * function.
+   */
+  template <int dim>
+    static 
+    void
+    parallelepiped (Triangulation<dim>   &tria,
+		    const Tensor<2, dim> &corners,
+		    const bool            colorize = false);
 
   /**
    * Hypercube with a layer of
diff --git a/deal.II/source/grid/grid_generator.cc b/deal.II/source/grid/grid_generator.cc
index b976ca453a..2d3e21fab2 100644
--- a/deal.II/source/grid/grid_generator.cc
+++ b/deal.II/source/grid/grid_generator.cc
@@ -373,9 +373,9 @@ GridGenerator::torus (Triangulation<2,3> &tria,
 template<>
 void
 GridGenerator::parallelogram (
-  Triangulation<2> &tria,
+  Triangulation<2>  &tria,
   const Tensor<2,2> &corners,
-  const bool      colorize)
+  const bool         colorize)
 {
   std::vector<Point<2> > vertices (GeometryInfo<2>::vertices_per_cell);
 
@@ -397,6 +397,79 @@ GridGenerator::parallelogram (
 
 
 
+// Parallelepiped implementation in 1d, 2d, and 3d. @note The
+// implementation in 1d is equivalent to hyper_rectangle, and in 2d is
+// equivalent to parallelogram. Reordering of vertices to lexiographic
+// occurs if needed. For this reason the triangulation is explicitly
+// constructed for each dim here.
+template<int dim>
+void
+GridGenerator::parallelepiped (Triangulation<dim>   &tria,
+			       const Tensor<2, dim> &corners,
+			       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 indicatorsaccording to hyper_rectangle
+  if (colorize)
+    colorize_hyper_rectangle (tria);
+}
+
+
+
+
 template <int dim>
 void
 GridGenerator::subdivided_hyper_cube (Triangulation<dim> &tria,
diff --git a/deal.II/source/grid/grid_generator.inst.in b/deal.II/source/grid/grid_generator.inst.in
index f16f7005fd..0561f0a4a4 100644
--- a/deal.II/source/grid/grid_generator.inst.in
+++ b/deal.II/source/grid/grid_generator.inst.in
@@ -64,7 +64,12 @@ for (deal_II_dimension : DIMENSIONS)
       const Point<deal_II_dimension>                &,
       const bool                       );
   
-  
+
+    template void
+      GridGenerator::parallelepiped<deal_II_dimension> (
+	Triangulation<deal_II_dimension>               &,
+	const Tensor<2, deal_II_dimension>             &, 
+	const bool);  
   
   
 #if deal_II_dimension > 1