]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Include Joergs subdivided_hyper_rectangle.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 9 Sep 2003 14:11:48 +0000 (14:11 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 9 Sep 2003 14:11:48 +0000 (14:11 +0000)
git-svn-id: https://svn.dealii.org/trunk@7966 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/grid/grid_generator.h
deal.II/deal.II/source/grid/grid_generator.cc

index 526716da81c22854a56285e43f3f61834039cdeb..2e2e7233ebaceae1543fcef6d2d346ae8f06b863 100644 (file)
@@ -49,6 +49,12 @@ template <typename number> class SparseMatrix;
  *      @ref{GridGenerator}@p{::hyper_rectangle (tria, p1, p2)}, with two
  *      opposite corner points @p{p1} and @p{p2}.
  *
+ *    @item Rectangular coordinate-parallel domains with different numbers of
+ *     cells in each coordinate direction are generated by 
+ *      @ref{GridGenerator}@p{::subdivided_hyper_rectangle 
+ *     (tria, int[dim] repetitions, p1, p2, colorize)}, with two
+ *      opposite corner points @p{p1} and @p{p2}.
+ *
  *    @item Generalized L-shape domain:
  *      using the @ref{GridGenerator}@p{::hyper_L (tria, a,b)} function produces
  *      the hypercube with the interval $[a,b]$ without the hypercube
@@ -84,7 +90,7 @@ template <typename number> class SparseMatrix;
  *      the triangulation to be refined at the outer boundaries. The class
  *      @ref{HalfHyperBallBoundary} will provide a boundary object.
  *
-*    @item Hyper shell: A hyper shell is the region between two hyper
+ *    @item Hyper shell: A hyper shell is the region between two hyper
  *      sphere with the same origin. Therefore, it is a ring in two
  *      spatial dimensions. To triangulation it, call the function
  *      @ref{GridGenerator}@p{::hyper_shell (tria, origin, inner_radius, outer_radius, N)},
@@ -131,7 +137,7 @@ template <typename number> class SparseMatrix;
  * transform (simple-shaped) grids to a more complicated ones, like a
  * shell onto a grid of an airfoil, for example.
  *
- * @author Wolfgang Bangerth, Ralf Hartmann, Guido Kanschat, Stefan Nauber, 1998, 1999, 2000, 2001, 2002, 2003.
+ * @author Wolfgang Bangerth, Ralf Hartmann, Guido Kanschat, Stefan Nauber, Joerg Weimar, 1998, 1999, 2000, 2001, 2002, 2003.
  */
 class GridGenerator
 {
@@ -205,6 +211,49 @@ class GridGenerator
                                 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},
+                                     * @p{repetitions[i]} cells are
+                                     * generated.
+                                     * 
+                                     * To get cells with an aspect
+                                     * ratio different from that of
+                                     * the recatangle, use different
+                                     * numbers of subdivisions in
+                                     * different coordinate
+                                     * directions. The minimum number
+                                     * of subdivisions in each
+                                     * direction is
+                                     * 1. @p{repetitions} is a list
+                                     * of integers denoting the
+                                     * number of subdivisions in each
+                                     * coordinate direction.
+                                     * 
+                                     * If the @p{colorize} flag is
+                                     * set, the
+                                     * @p{boundary_indicator}s of the
+                                     * surfaces are assigned, such
+                                     * that the lower one in
+                                     * @p{x}-direction is 0, the
+                                     * upper one is 1. The indicators
+                                     * for the surfaces in
+                                     * @p{y}-direction are 2 and 3,
+                                     * the ones for @p{z} are 4 and
+                                     * 5.
+                                     */
+    template <int dim>
+    static
+    void
+    subdivided_hyper_rectangle (Triangulation<dim>              &tria,
+                               const std::vector<unsigned int> &repetitions,
+                               const Point<dim>                &p1,
+                               const Point<dim>                &p2,
+                               const bool                       colorize=false);
+
                                     /**
                                      * Hypercube with a layer of
                                      * hypercubes around it. The
@@ -591,6 +640,13 @@ class GridGenerator
                     int,
                     << "The number of repetitions " << arg1
                     << " must be >=1.");
+                                    /**
+                                     * Exception
+                                     */
+    DeclException1 (ExcInvalidRepetitionsDimension,
+                    int,
+                    << "The vector of repetitions  must have " 
+                   << arg1 <<" elements.");
 
   private:
                                     /**
@@ -608,6 +664,35 @@ class GridGenerator
                                      */
     static void colorize_hyper_rectangle (Triangulation<1> &tria);
     
+                                    /**
+                                     * Perform the action specified
+                                     * by the @p{colorize} flag of
+                                     * the
+                                     * @ref{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);
+
+                                    /**
+                                     * Declaration of the above
+                                     * function for 1d.
+                                     */
+    static
+    void
+    colorize_subdivided_hyper_rectangle (Triangulation<1> &tria,
+                                        const Point<1>   &p1,
+                                        const Point<1>   &p2,
+                                        const double      epsilon);
+
                                     /**
                                      * Solve the Laplace equation for
                                      * @p{laplace_transformation}
@@ -617,9 +702,11 @@ class GridGenerator
                                      * a function of its own in order
                                      * to allow parallel execution.
                                      */
-    static void laplace_solve (const SparseMatrix<double>          &S,
-                              const std::map<unsigned int,double> &m,
-                              Vector<double>                      &u);
+    static
+    void
+    laplace_solve (const SparseMatrix<double>          &S,
+                  const std::map<unsigned int,double> &m,
+                  Vector<double>                      &u);
 };
 
 
index 8645257c57ff18e786c1b4f42bdd2911901901bd..f55c3224d5d5223e160df9dcf9139ea9fe5f08f3 100644 (file)
@@ -243,6 +243,203 @@ GridGenerator::subdivided_hyper_cube (Triangulation<dim> &tria,
 }
 
 
+
+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.
+  double delta[dim];
+  
+  for (unsigned int i=0; i<dim; ++i)
+    {
+      Assert (repetitions[i] >= 1, ExcInvalidRepetitions(repetitions[i]));
+               
+      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]));
+            break;
+
+      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]));
+            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());
+    }
+
+                                   // 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;
+
+      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];
+                  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+1;
+                  cells[c].vertices[3] = (y+1)*(repetitions[0]+1)+x;
+                  cells[c].material_id = 0;
+                }
+            break;
+
+      default:                                 
+                                            // copied from hyper_rectangle:
+                                             // 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());  
+
+  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 epsilon = 0.01*delta[0];
+    if (dim > 1)
+      epsilon = std::min(epsilon,0.01*delta[1]);
+    
+    if (dim > 2)
+      epsilon = std::min(epsilon,0.01*delta[2]);
+    
+                                    // actual code is external since
+                                    // 1-D is different from 2/3D.
+    colorize_subdivided_hyper_rectangle (tria,p1,p2,epsilon);
+  }
+}
+
+
+#if deal_II_dimension == 1
+
+void
+GridGenerator::colorize_subdivided_hyper_rectangle (Triangulation<1> &,
+                                                   const Point<1>   &,
+                                                   const Point<1>   &,
+                                                   const double      )
+{
+                                  // nothing to do in 1d
+                                  // boundary indicators are set to
+                                  // 0 (left) and 1 (right) by default.
+}
+
+#else
+
+template <int dim>
+void
+GridGenerator::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>::raw_face_iterator face = tria.begin_raw_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());
+         
+       }
+    }
+}
+
+#endif
+
+
+
 #if deal_II_dimension == 1
 
 void GridGenerator::hyper_cube_slit (Triangulation<1> &,
@@ -348,15 +545,15 @@ GridGenerator::hyper_cube_slit (Triangulation<2> &tria,
 {
   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 ) };
+                                 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,2,3 },
                                    { 9,4,5,2 },
                                    { 3,2,6,7 },
@@ -390,13 +587,13 @@ GridGenerator::hyper_L (Triangulation<2> &tria,
 {
   const unsigned int dim=2;
   const Point<dim> vertices[8] = { Point<dim> (a,a),
-                                    Point<dim> ((a+b)/2,a),
-                                    Point<dim> (b,a),
-                                    Point<dim> (a,(a+b)/2),
-                                    Point<dim> ((a+b)/2,(a+b)/2),
-                                    Point<dim> (b,(a+b)/2),
-                                    Point<dim> (a,b),
-                                    Point<dim> ((a+b)/2,b)  };
+                                  Point<dim> ((a+b)/2,a),
+                                  Point<dim> (b,a),
+                                  Point<dim> (a,(a+b)/2),
+                                  Point<dim> ((a+b)/2,(a+b)/2),
+                                  Point<dim> (b,(a+b)/2),
+                                  Point<dim> (a,b),
+                                  Point<dim> ((a+b)/2,b)  };
   const int cell_vertices[3][4] = {{0, 1, 4, 3},
                                   {1, 2, 5, 4},
                                   {3, 4, 7, 6}};
@@ -427,13 +624,13 @@ GridGenerator::hyper_ball (Triangulation<2> &tria,
                                   // 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)) };
+                                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, 3, 2},
                                   {0, 2, 4, 6},
@@ -534,14 +731,14 @@ GridGenerator::cylinder (Triangulation<2> &tria,
       switch (f->boundary_indicator())
        {
          case 0:
-           f->set_boundary_indicator(1);
-           break;
+               f->set_boundary_indicator(1);
+               break;
          case 1:
-           f->set_boundary_indicator(2);
-           break;
+               f->set_boundary_indicator(2);
+               break;
          default:
-           f->set_boundary_indicator(0);
-           break;          
+               f->set_boundary_indicator(0);
+               break;      
        }
       ++f;
     }
@@ -559,13 +756,13 @@ GridGenerator::half_hyper_ball (Triangulation<2> &tria,
                                   // 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)) };
+                                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, 3, 2},
                                   {2, 3, 5, 4},
@@ -739,37 +936,37 @@ GridGenerator::hyper_L (Triangulation<3> &tria,
                                   // part of the cube
   const Point<dim> vertices[26]
     = {
-                                    // front face of the big cube
-      Point<dim> (a,      a,a),
-      Point<dim> ((a+b)/2,a,a),
-      Point<dim> (b,      a,a),
-      Point<dim> (a,      a,(a+b)/2),
-      Point<dim> ((a+b)/2,a,(a+b)/2),
-      Point<dim> (b,      a,(a+b)/2),
-      Point<dim> (a,      a,b),
-      Point<dim> ((a+b)/2,a,b),
-      Point<dim> (b,      a,b),
-                                      // middle face of the big cube
-      Point<dim> (a,      (a+b)/2,a),
-      Point<dim> ((a+b)/2,(a+b)/2,a),
-      Point<dim> (b,      (a+b)/2,a),
-      Point<dim> (a,      (a+b)/2,(a+b)/2),
-      Point<dim> ((a+b)/2,(a+b)/2,(a+b)/2),
-      Point<dim> (b,      (a+b)/2,(a+b)/2),
-      Point<dim> (a,      (a+b)/2,b),
-      Point<dim> ((a+b)/2,(a+b)/2,b),
-      Point<dim> (b,      (a+b)/2,b),
-                                      // back face of the big cube
-                                      // last (top right) point is missing
-      Point<dim> (a,      b,a),
-      Point<dim> ((a+b)/2,b,a),
-      Point<dim> (b,      b,a),
-      Point<dim> (a,      b,(a+b)/2),
-      Point<dim> ((a+b)/2,b,(a+b)/2),
-      Point<dim> (b,      b,(a+b)/2),
-      Point<dim> (a,      b,b),
-      Point<dim> ((a+b)/2,b,b)
-      };
+                                          // front face of the big cube
+         Point<dim> (a,      a,a),
+         Point<dim> ((a+b)/2,a,a),
+         Point<dim> (b,      a,a),
+         Point<dim> (a,      a,(a+b)/2),
+         Point<dim> ((a+b)/2,a,(a+b)/2),
+         Point<dim> (b,      a,(a+b)/2),
+         Point<dim> (a,      a,b),
+         Point<dim> ((a+b)/2,a,b),
+         Point<dim> (b,      a,b),
+                                          // middle face of the big cube
+         Point<dim> (a,      (a+b)/2,a),
+         Point<dim> ((a+b)/2,(a+b)/2,a),
+         Point<dim> (b,      (a+b)/2,a),
+         Point<dim> (a,      (a+b)/2,(a+b)/2),
+         Point<dim> ((a+b)/2,(a+b)/2,(a+b)/2),
+         Point<dim> (b,      (a+b)/2,(a+b)/2),
+         Point<dim> (a,      (a+b)/2,b),
+         Point<dim> ((a+b)/2,(a+b)/2,b),
+         Point<dim> (b,      (a+b)/2,b),
+                                          // back face of the big cube
+                                          // last (top right) point is missing
+         Point<dim> (a,      b,a),
+         Point<dim> ((a+b)/2,b,a),
+         Point<dim> (b,      b,a),
+         Point<dim> (a,      b,(a+b)/2),
+         Point<dim> ((a+b)/2,b,(a+b)/2),
+         Point<dim> (b,      b,(a+b)/2),
+         Point<dim> (a,      b,b),
+         Point<dim> ((a+b)/2,b,b)
+    };
   const int cell_vertices[7][8] = {{0, 1, 4, 3, 9, 10, 13, 12},
                                   {1, 2, 5, 4, 10, 11, 14, 13},
                                   {3, 4, 7, 6, 12, 13, 16, 15},
@@ -805,27 +1002,27 @@ GridGenerator::hyper_ball (Triangulation<3> &tria,
   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)),
-      };
+                                          // 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
@@ -866,31 +1063,31 @@ GridGenerator::cylinder (Triangulation<3> &tria,
   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),
-      };
+       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)
     {
@@ -995,7 +1192,7 @@ void GridGenerator::laplace_transformation (Triangulation<dim> &tria,
                                   // fill these maps using the data
                                   // given by new_points
   typename DoFHandler<dim>::cell_iterator cell=dof_handler.begin_active(),
-                                endc=dof_handler.end();
+                                         endc=dof_handler.end();
   typename DoFHandler<dim>::face_iterator face;
   for (; cell!=endc; ++cell)
     {
@@ -1013,7 +1210,7 @@ void GridGenerator::laplace_transformation (Triangulation<dim> &tria,
                  
                  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)));
+                                 face->vertex_dof_index(vertex_no, 0), map_iter->second(i)));
                }
          }
     }
@@ -1088,6 +1285,14 @@ GridGenerator::subdivided_hyper_cube<deal_II_dimension> (Triangulation<deal_II_d
                                                          const double,
                                                          const double);
 
+template void
+GridGenerator::subdivided_hyper_rectangle<deal_II_dimension> 
+(Triangulation<deal_II_dimension> &,
+ const std::vector<unsigned int>&,
+ const Point<deal_II_dimension>&,
+ const Point<deal_II_dimension>&,
+ const bool);
+
 
 #if deal_II_dimension != 1
 template void

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.