]> https://gitweb.dealii.org/ - dealii.git/commitdiff
ConeBoundary and GridGenerator::truncated_cone. Implementation by Markus Buerg.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 17 Aug 2009 17:26:47 +0000 (17:26 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 17 Aug 2009 17:26:47 +0000 (17:26 +0000)
git-svn-id: https://svn.dealii.org/trunk@19288 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/grid/grid_generator.h
deal.II/deal.II/include/grid/tria_boundary_lib.h
deal.II/deal.II/source/grid/grid_generator.cc
deal.II/deal.II/source/grid/tria_boundary_lib.cc
deal.II/doc/authors.html
deal.II/doc/news/changes.h

index 33b297fe7569809ad7684509892b1b3ca4786df3..39ef2464203baa06f7c485dfdc6b52956e515b9e 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors
+//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -37,7 +37,7 @@ template <typename number> class SparseMatrix;
  * 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
@@ -147,7 +147,7 @@ class GridGenerator
                                      * adds 2<sup>i</sup>. For
                                      * instance, the center point
                                      * (1,-1,1) yields a material id 5.
-                                     * 
+                                     *
                                      * @note The triangulation needs to be
                                      * void upon calling this
                                      * function.
@@ -166,7 +166,7 @@ class GridGenerator
                                      * dimension @p i,
                                      * <tt>repetitions[i]</tt> cells are
                                      * generated.
-                                     * 
+                                     *
                                      * To get cells with an aspect
                                      * ratio different from that of
                                      * the domain, use different
@@ -179,7 +179,7 @@ class GridGenerator
                                      * of integers denoting the
                                      * number of subdivisions in each
                                      * coordinate direction.
-                                     * 
+                                     *
                                      * If the @p colorize flag is
                                      * set, the
                                      * @p boundary_indicators of the
@@ -279,7 +279,7 @@ class GridGenerator
                                const Point<dim>                         &p,
                                const Table<dim,unsigned char>           &material_id,
                                const bool                               colorize=false);
-    
+
                                     /**
                                      * A parallelogram. The first
                                      * corner point is the
@@ -304,8 +304,8 @@ class GridGenerator
     parallelogram(Triangulation<dim>&  tria,
                  const Tensor<2,dim>& corners,
                  const bool           colorize=false);
-                 
-                             
+
+
                                     /**
                                      * Hypercube with a layer of
                                      * hypercubes around it. The
@@ -339,7 +339,7 @@ class GridGenerator
                                     const double      right= 1.,
                                     const double      thickness = 1.,
                                     const bool        colorize = false);
-    
+
                                     /**
                                      * Initialize the given
                                      * triangulation with a
@@ -366,7 +366,7 @@ class GridGenerator
                                      * @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>(),
@@ -381,14 +381,14 @@ class GridGenerator
                                      * perpendicular to the
                                      * <i>x</i>-axis.
                                      *
-                                     * The boundary indicators for the final 
+                                     * The boundary indicators for the final
                                      * triangulation are 0 for the curved boundary and
-                                     * 1 for the cut plane. 
+                                     * 1 for the cut plane.
                                      *
                                      * The appropriate
-                                     * boundary class is 
+                                     * boundary class is
                                      * HalfHyperBallBoundary, or HyperBallBoundary.
-                                     * 
+                                     *
                                      * @note The triangulation needs to be
                                      * void upon calling this
                                      * function.
@@ -429,7 +429,55 @@ class GridGenerator
     static 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 <tt>x=-half_length</tt>
+                                      * to <tt>x=half_length</tt> and
+                                      * its projection into the @p
+                                      * yz-plane is a circle of radius
+                                      * @p radius_0 at
+                                      * <tt>x=-half_length</tt> and a
+                                      * circle of radius @p radius_1
+                                      * at <tt>x=+half_length</tt>.
+                                      * In between the radius is
+                                      * linearly decreasing.
+                                      *
+                                      * In two dimensions, the cone is
+                                      * a trapezoid from
+                                      * <tt>x=-half_length</tt> to
+                                      * <tt>x=+half_length</tt> and
+                                      * from <tt>y=-radius</tt> to
+                                      * <tt>y=radius_0</tt> at
+                                      * <tt>x=-half_length</tt> and
+                                      * from <tt>y=-radius_1</tt> to
+                                      * <tt>y=radius_1</tt> at
+                                      * <tt>x=+half_length</tt>.  In
+                                      * between the range of
+                                      * <tt>y</tt> is linearly
+                                      * decreasing.
+                                     *
+                                      * The boundaries are colored
+                                      * according to the following
+                                      * scheme: 0 for the hull of the
+                                      * cone, 1 for the left hand
+                                      * face and 2 for the right hand
+                                      * face.
+                                      *
+                                      * @note The triangulation needs to be
+                                      * void upon calling this
+                                      * function.
+                                     *
+                                     * @author Markus B&uuml;rg, 2009
+                                      */
+    template <int dim>
+    static void
+    truncated_cone (Triangulation<dim> &tria,
+                   const double radius_0 = 1.0,
+                   const double radius_1 = 0.5,
+                   const double half_length = 1.0);
+
                                     /**
                                      * Initialize the given
                                      * triangulation with a hyper-L
@@ -460,7 +508,7 @@ class GridGenerator
     static void hyper_L (Triangulation<dim> &tria,
                         const double      left = -1.,
                         const double      right= 1.);
-    
+
                                      /**
                                      * Initialize the given
                                      * Triangulation with a hypercube
@@ -481,7 +529,7 @@ class GridGenerator
                                      * <i>z</i>-direction, such that
                                      * a plane cuts the lower half of
                                      * a rectangle in two.
-                                     
+
                                      * This function is declared to
                                      * exist for triangulations of
                                      * all space dimensions, but
@@ -497,7 +545,7 @@ class GridGenerator
                                 const double      left = 0.,
                                 const double      right= 1.,
                                 const bool colorize = false);
-    
+
                                     /**
                                      * Produce a hyper-shell,
                                      * the region between two
@@ -550,7 +598,7 @@ class GridGenerator
                             const double        outer_radius,
                             const unsigned int  n_cells = 0,
                             bool colorize = false);
-    
+
                                     /**
                                      * Produce a half hyper-shell,
                                      * i.e. the space between two
@@ -595,7 +643,7 @@ class GridGenerator
                                  const double        inner_radius,
                                  const double        outer_radius,
                                  const unsigned int  n_cells = 0);
-    
+
                                     /**
                                      * Produce a domain that is the space
                                      * between two cylinders in 3d, with
@@ -626,7 +674,7 @@ class GridGenerator
                                 const unsigned int  n_radial_cells = 0,
                                 const unsigned int  n_axial_cells = 0);
 
-                                    /** 
+                                    /**
                                      * This class produces a square
                                      * on the <i>xy</i>-plane with a
                                      * circular hole in the middle,
@@ -634,29 +682,29 @@ class GridGenerator
                                      * (only in 3d).
                                      *
                                      *  @image html cubes_hole.png
-                                     *  
+                                     *
                                      * It is implemented in 2d and
                                      * 3d, and takes the following
                                      * arguments:
-                                     * 
+                                     *
                                      * @arg @p inner_radius: size of the
-                                      *    internal hole 
+                                      *    internal hole
                                      * @arg @p  outer_radius: size of the
                                       *    biggest enclosed cylinder
                                      * @arg @p L: extension on the @p z-direction
                                      * @arg @p repetitions: number of subdivisions
                                      *      along the @p z-direction
-                                     * @arg @p colorize: wether to assign different 
+                                     * @arg @p colorize: wether to assign different
                                      *     boundary indicators to different faces.
-                                     *    The colors are given in lexicographic 
-                                     *    ordering for the flat faces (0 to 3 in 2d, 
-                                     *    0 to 5 in 3d) plus the curved hole 
+                                     *    The colors are given in lexicographic
+                                     *    ordering for the flat faces (0 to 3 in 2d,
+                                     *    0 to 5 in 3d) plus the curved hole
                                      *    (4 in 2d, and 6 in 3d).
-                                     *    If @p colorize is set to false, then flat faces 
+                                     *    If @p colorize is set to false, then flat faces
                                      *    get the number 0 and the hole gets number 1.
                                      */
     template<int dim>
-    static void hyper_cube_with_cylindrical_hole (Triangulation<dim> &triangulation, 
+    static void hyper_cube_with_cylindrical_hole (Triangulation<dim> &triangulation,
                                                const double inner_radius = .25,
                                                const double outer_radius = .5,
                                                const double L = .5,
@@ -680,7 +728,7 @@ class GridGenerator
                         const unsigned int   n_rotations,
                         const double         R,
                         const double         r);
-    
+
                                     /**
                                      * This function transformes the
                                      * @p Triangulation @p tria
@@ -723,7 +771,7 @@ class GridGenerator
                                      */
     DeclException1 (ExcInvalidRepetitionsDimension,
                     int,
-                    << "The vector of repetitions  must have " 
+                    << "The vector of repetitions  must have "
                    << arg1 <<" elements.");
 
   private:
@@ -766,7 +814,7 @@ class GridGenerator
                          const Point<dim>& center,
                          const double inner_radius,
                          const double outer_radius);
-    
+
                                     /**
                                      * Solve the Laplace equation for
                                      * @p laplace_transformation
index bf4ad320c8ed974970d1b3d188033c377f808fb3..62c9c8bae2d826a6aa734f3a0a0721dafe42f03b 100644 (file)
@@ -142,8 +142,8 @@ class CylinderBoundary : public StraightBoundary<dim>
                                      * see below, flag is set true.
                                      */
     DeclException0 (ExcRadiusNotSet);
-    
-    
+
+
   protected:
                                     /**
                                      * Radius of the cylinder.
@@ -188,10 +188,163 @@ class CylinderBoundary : public StraightBoundary<dim>
 
 
 /**
- *   Specialisation of Boundary<dim>, which places the new point on
+ * Boundary object for the hull of a (truncated) cone with two
+ * different radii at the two ends. If one radius is chosen to be 0
+ * the object describes the boundary of a cone. In three dimensions,
+ * points are projected on an arbitrarily oriented (truncated) cone
+ * described by the two endpoints and the corresponding radii. Similar
+ * to HyperBallBoundary, new points are projected by dividing the
+ * straight line between the old two points and adjusting the radius
+ * from the axis.
+ *
+ * This class is derived from StraightBoundary rather than from
+ * Boundary, which would seem natural, since this way we can use the
+ * StraightBoundary::in_between() function.
+ *
+ * @author Markus B&uuml;rg, 2009
+ */
+template <int dim>
+class ConeBoundary : public StraightBoundary<dim>
+{
+  public:
+                                    /**
+                                     * Constructor. Here the boundary
+                                     * object is constructed. The
+                                     * points <tt>x_0</tt> and
+                                     * <tt>x_1</tt> describe the
+                                     * starting and ending points of
+                                     * the axis of the (cutted)
+                                     * cone. <tt>radius_0</tt>
+                                     * denotes the radius
+                                     * corresponding to <tt>x_0</tt>
+                                     * and <tt>radius_1</tt> the one
+                                     * corresponding to <tt>x_1</tt>.
+                                     */
+    ConeBoundary (const double radius_0,
+                 const double radius_1,
+                 const Point<dim> x_0,
+                 const Point<dim> x_1);
+
+                                    /**
+                                     * Return the radius of the
+                                     * (cutted) cone at given point
+                                     * <tt>x</tt> on the axis.
+                                     */
+    double get_radius (const Point<dim> x) const;
+
+                                    /**
+                                     * Refer to the general
+                                     * documentation of this class
+                                     * and the documentation of the
+                                     * base class.
+                                     */
+    virtual
+    Point<dim>
+    get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const;
+
+                                    /**
+                                     * Refer to the general
+                                     * documentation of this class
+                                     * and the documentation of the
+                                     * base class.
+                                     */
+    virtual
+    Point<dim>
+    get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const;
+
+                                    /**
+                                     * Refer to the general
+                                     * documentation of this class
+                                     * and the documentation of the
+                                     * base class.
+                                     *
+                                     * Calls @p
+                                     * get_intermediate_points_between_points.
+                                     */
+    virtual
+    void
+    get_intermediate_points_on_line (const typename Triangulation<dim>::line_iterator &line,
+                                    std::vector<Point<dim> > &points) const;
+
+                                    /**
+                                     * Refer to the general
+                                     * documentation of this class
+                                     * and the documentation of the
+                                     * base class.
+                                     *
+                                     * Only implemented for
+                                     * <tt>dim=3</tt> and for
+                                     * <tt>points.size()==1</tt>.
+                                     */
+    virtual
+    void
+    get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterator &quad,
+                                    std::vector<Point<dim> > &points) const;
+
+                                    /**
+                                     * Compute the normals to the
+                                     * boundary at the vertices of
+                                     * the given face.
+                                     *
+                                     * Refer to the general
+                                     * documentation of this class
+                                     * and the documentation of the
+                                     * base class.
+                                     */
+    virtual
+    void
+    get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &face,
+                            typename Boundary<dim>::FaceVertexNormals &face_vertex_normals) const;
+
+  protected:
+                                    /**
+                                     * First radius of the (cutted)
+                                     * cone.
+                                     */
+    const double radius_0;
+
+                                    /**
+                                     * Second radius of the (cutted)
+                                     * cone.
+                                     */
+    const double radius_1;
+
+                                    /**
+                                     * Starting point of the axis.
+                                     */
+    const Point<dim> x_0;
+
+                                    /**
+                                     * Ending point of the axis.
+                                     */
+    const Point<dim> x_1;
+
+  private:
+                                    /**
+                                     * Called by @p
+                                     * get_intermediate_points_on_line
+                                     * and by @p
+                                     * get_intermediate_points_on_quad.
+                                     *
+                                     * Refer to the general
+                                     * documentation of @p
+                                     * get_intermediate_points_on_line
+                                     * in the documentation of the
+                                     * base class.
+                                     */
+    void
+    get_intermediate_points_between_points (const Point<dim> &p0,
+                                           const Point<dim> &p1,
+                                           std::vector<Point<dim> > &points) const;
+};
+
+
+
+/**
+ *   Specialization of Boundary<dim>, which places the new point on
  *   the boundary of a ball in arbitrary dimension. It works by projecting
  *   the point in the middle of the old points onto the ball. The middle is
- *   defined as the arithmetic mean of the points. 
+ *   defined as the arithmetic mean of the points.
  *
  *   The center of the ball and its radius may be given upon construction of
  *   an object of this type. They default to the origin and a radius of 1.0.
@@ -287,10 +440,10 @@ class HyperBallBoundary : public StraightBoundary<dim,spacedim>
                                      * see below, flag is set true.
                                      */
     DeclException0 (ExcRadiusNotSet);
-    
-    
+
+
   protected:
-    
+
                                     /**
                                      * Center point of the hyperball.
                                      */
@@ -335,7 +488,7 @@ class HyperBallBoundary : public StraightBoundary<dim,spacedim>
                                      * base class.
                                      */
     void get_intermediate_points_between_points (const Point<spacedim> &p0, const Point<spacedim> &p1,
-                                                std::vector<Point<spacedim> > &points) const;    
+                                                std::vector<Point<spacedim> > &points) const;
 };
 
 
@@ -464,7 +617,7 @@ class HyperShellBoundary : public HyperBallBoundary<dim>
  * @author Wolfgang Bangerth, 2000, 2009
  */
 template <int dim>
-class HalfHyperShellBoundary : public HyperShellBoundary<dim> 
+class HalfHyperShellBoundary : public HyperShellBoundary<dim>
 {
   public:
                                     /**
@@ -482,13 +635,13 @@ class HalfHyperShellBoundary : public HyperShellBoundary<dim>
     HalfHyperShellBoundary (const Point<dim> &center = Point<dim>(),
                            const double inner_radius = -1,
                            const double outer_radius = -1);
-    
+
                                     /**
                                      * Construct a new point on a line.
                                      */
     virtual Point<dim>
-    get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const;  
-    
+    get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const;
+
                                     /**
                                      * Construct a new point on a quad.
                                      */
index 3d1135478a778d30c05fadf5e9ca83a67f4223fd..586f3dcd26b36c248d44b9bb98408c1febcec3a9 100644 (file)
@@ -70,7 +70,7 @@ namespace
        Point<3>( 0, 0,-1),
        Point<3>( 0, 0, 1)
   };
-  
+
 #endif
 }
 
@@ -86,11 +86,11 @@ GridGenerator::hyper_rectangle (Triangulation<dim,spacedim> &tria,
                                   // 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)
     {
@@ -101,24 +101,24 @@ GridGenerator::hyper_rectangle (Triangulation<dim,spacedim> &tria,
       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 ());
@@ -174,7 +174,7 @@ void GridGenerator::hyper_cube (Triangulation<dim,spacedim> &tria,
 {
   Assert (left < right,
          ExcMessage ("Invalid left and right bounds of hypercube"));
-  
+
   Point<spacedim> p1;
   Point<spacedim> p2;
 
@@ -204,12 +204,12 @@ GridGenerator::moebius (
   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."));
-    
-  
+
+
   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;
-  
+
   for (unsigned int i=0; i<n_cells; ++i)
     for (unsigned int j=0; j<4; ++j)
       {
@@ -219,7 +219,7 @@ GridGenerator::moebius (
       }
 
   unsigned int offset=0;
-  
+
   std::vector<CellData<dim> > cells (n_cells);
   for (unsigned int i=0; i<n_cells; ++i)
     {
@@ -233,7 +233,7 @@ GridGenerator::moebius (
       offset+=4;
       cells[i].material_id=0;
     }
-  
+
                                   // 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;
@@ -257,7 +257,7 @@ GridGenerator::parallelogram (
   const bool      colorize)
 {
   Assert (dim==2, ExcNotImplemented());
-  
+
   std::vector<Point<dim> > vertices (GeometryInfo<dim>::vertices_per_cell);
 
   vertices[1] = corners[0];
@@ -288,7 +288,7 @@ GridGenerator::subdivided_hyper_cube (Triangulation<dim> &tria,
   Assert (repetitions >= 1, ExcInvalidRepetitions(repetitions));
   Assert (left < right,
          ExcMessage ("Invalid left and right bounds of hypercube"));
-  
+
                                    // first generate the necessary
                                    // points
   const double delta = (right-left)/repetitions;
@@ -340,7 +340,7 @@ GridGenerator::subdivided_hyper_cube (Triangulation<dim> &tria,
            cells[x].material_id = 0;
          }
        break;
-       
+
       case 2:
        cells.resize (repetitions*repetitions);
        for (unsigned int y=0; y<repetitions; ++y)
@@ -386,7 +386,7 @@ GridGenerator::subdivided_hyper_cube (Triangulation<dim> &tria,
             Assert (false, ExcNotImplemented());
     }
 
-  tria.create_triangulation (points, cells, SubCellData());  
+  tria.create_triangulation (points, cells, SubCellData());
 }
 
 
@@ -402,14 +402,14 @@ GridGenerator::subdivided_hyper_rectangle (
 {
                                   // contributed by Joerg R. Weimar
                                   // (j.weimar@jweimar.de) 2003
-  Assert(repetitions.size() == dim, 
+  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));
@@ -419,14 +419,14 @@ GridGenerator::subdivided_hyper_rectangle (
                                   // 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]));
-               
+
       delta[i] = (p2[i]-p1[i])/repetitions[i];
     }
+
                                    // then generate the necessary
                                    // points
   std::vector<Point<dim> > points;
@@ -473,7 +473,7 @@ GridGenerator::subdivided_hyper_rectangle (
           }
         break;
       }
-      
+
       case 2:
       {
         cells.resize (repetitions[1]*repetitions[0]);
@@ -489,12 +489,12 @@ GridGenerator::subdivided_hyper_rectangle (
             }
         break;
       }
-      
+
       case 3:
       {
         const unsigned int n_x  = (repetitions[0]+1);
         const unsigned int n_xy = (repetitions[0]+1)*(repetitions[1]+1);
-        
+
         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)
@@ -513,14 +513,14 @@ GridGenerator::subdivided_hyper_rectangle (
                 cells[c].material_id = 0;
               }
         break;
-        
+
       }
 
       default:
             Assert (false, ExcNotImplemented());
     }
 
-  tria.create_triangulation (points, cells, SubCellData());  
+  tria.create_triangulation (points, cells, SubCellData());
 
   if (colorize)
     {
@@ -528,7 +528,7 @@ GridGenerator::subdivided_hyper_rectangle (
                                       // 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.
@@ -536,7 +536,7 @@ GridGenerator::subdivided_hyper_rectangle (
         = 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);
@@ -557,7 +557,7 @@ subdivided_hyper_rectangle(Triangulation<dim>              &tria,
                                   // contributed by Joerg R. Weimar
                                   // (j.weimar@jweimar.de) 2003
                                   // modified by Yaqi Wang 2006
-  Assert(step_sz.size() == dim, 
+  Assert(step_sz.size() == dim,
         ExcInvalidRepetitionsDimension(dim));
 
 
@@ -582,13 +582,13 @@ subdivided_hyper_rectangle(Triangulation<dim>              &tria,
        }
 
       double x = 0;
-      for (unsigned int j=0; j<step_sizes.at(i).size(); j++) 
+      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) );
     }
 
+
                                    // then generate the necessary
                                    // points
   std::vector<Point<dim> > points;
@@ -598,7 +598,7 @@ subdivided_hyper_rectangle(Triangulation<dim>              &tria,
       {
        double x=0;
        for (unsigned int i=0; ; ++i)
-         { 
+         {
            points.push_back (Point<dim> (p1[0]+x));
 
                                             // form partial sums. in
@@ -608,12 +608,12 @@ subdivided_hyper_rectangle(Triangulation<dim>              &tria,
                                             // and exit early instead
            if (i == step_sizes[0].size())
              break;
-           
+
            x += step_sizes[0][i];
          }
        break;
       }
-      
+
       case 2:
       {
        double y=0;
@@ -626,13 +626,13 @@ subdivided_hyper_rectangle(Triangulation<dim>              &tria,
                                              p1[1]+y));
                if (i == step_sizes[0].size())
                  break;
-               
+
                x += step_sizes[0][i];
              }
-           
+
            if (j == step_sizes[1].size())
              break;
-           
+
            y += step_sizes[1][j];
          }
        break;
@@ -654,24 +654,24 @@ subdivided_hyper_rectangle(Triangulation<dim>              &tria,
                                                  p1[2]+z));
                    if (i == step_sizes[0].size())
                      break;
-                   
+
                    x += step_sizes[0][i];
                  }
 
                if (j == step_sizes[1].size())
                  break;
-               
+
                y += step_sizes[1][j];
              }
-           
+
            if (k == step_sizes[2].size())
              break;
-           
+
            z += step_sizes[2][k];
          }
        break;
       }
-      
+
       default:
             Assert (false, ExcNotImplemented());
     }
@@ -692,7 +692,7 @@ subdivided_hyper_rectangle(Triangulation<dim>              &tria,
           }
         break;
       }
-      
+
       case 2:
       {
         cells.resize (step_sizes[1].size()*step_sizes[0].size());
@@ -708,12 +708,12 @@ subdivided_hyper_rectangle(Triangulation<dim>              &tria,
             }
         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);
-        
+
         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)
@@ -732,14 +732,14 @@ subdivided_hyper_rectangle(Triangulation<dim>              &tria,
                 cells[c].material_id = 0;
               }
         break;
-        
+
       }
 
       default:
             Assert (false, ExcNotImplemented());
     }
 
-  tria.create_triangulation (points, cells, SubCellData());  
+  tria.create_triangulation (points, cells, SubCellData());
 
   if (colorize)
     {
@@ -747,7 +747,7 @@ subdivided_hyper_rectangle(Triangulation<dim>              &tria,
                                       // 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.
@@ -758,7 +758,7 @@ subdivided_hyper_rectangle(Triangulation<dim>              &tria,
                             *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);
@@ -779,12 +779,12 @@ subdivided_hyper_rectangle (Triangulation<1>&                             tria,
                             const bool                                    colorize)
 {
                                   // contributed by Yaqi Wang 2006
-  Assert(spacing.size() == 1, 
+  Assert(spacing.size() == 1,
         ExcInvalidRepetitionsDimension(1));
 
   const unsigned int n_cells = material_id.size(0);
 
-  Assert(spacing[0].size() == n_cells, 
+  Assert(spacing[0].size() == n_cells,
         ExcInvalidRepetitionsDimension(1));
 
   double delta = std::numeric_limits<double>::max();
@@ -792,7 +792,7 @@ subdivided_hyper_rectangle (Triangulation<1>&                             tria,
     Assert (spacing[0][i] >= 0, ExcInvalidRepetitions(-1));
     delta = std::min (delta, spacing[0][i]);
   }
-  
+
                                    // generate the necessary points
   std::vector<Point<1> > points;
   double ax = p[0];
@@ -820,7 +820,7 @@ subdivided_hyper_rectangle (Triangulation<1>&                             tria,
   SubCellData t;
   GridTools::delete_unused_vertices (points, cells, t);
 
-  tria.create_triangulation (points, cells, t);  
+  tria.create_triangulation (points, cells, t);
 
                                    // set boundary indicator
   if (colorize)
@@ -841,7 +841,7 @@ subdivided_hyper_rectangle (Triangulation<2>&                         tria,
                             const bool                                    colorize)
 {
                                   // contributed by Yaqi Wang 2006
-  Assert(spacing.size() == 2, 
+  Assert(spacing.size() == 2,
         ExcInvalidRepetitionsDimension(2));
 
   std::vector<unsigned int> repetitions(2);
@@ -856,10 +856,10 @@ subdivided_hyper_rectangle (Triangulation<2>&                         tria,
          Assert (spacing[i][j] >= 0, ExcInvalidRepetitions(-1));
          delta = std::min (delta, spacing[i][j]);
        }
-      Assert(material_id.size(i) == repetitions[i], 
+      Assert(material_id.size(i) == repetitions[i],
             ExcInvalidRepetitionsDimension(i));
     }
+
                                    // generate the necessary points
   std::vector<Point<2> > points;
   double ay = p[1];
@@ -901,7 +901,7 @@ subdivided_hyper_rectangle (Triangulation<2>&                         tria,
   SubCellData t;
   GridTools::delete_unused_vertices (points, cells, t);
 
-  tria.create_triangulation (points, cells, t);  
+  tria.create_triangulation (points, cells, t);
 
                                     // set boundary indicator
   if (colorize)
@@ -943,8 +943,8 @@ subdivided_hyper_rectangle (Triangulation<3>&                           tria,
 {
                                   // contributed by Yaqi Wang 2006
   const unsigned int dim = 3;
-  
-  Assert(spacing.size() == dim, 
+
+  Assert(spacing.size() == dim,
         ExcInvalidRepetitionsDimension(dim));
 
   std::vector<unsigned int> repetitions(dim);
@@ -959,7 +959,7 @@ subdivided_hyper_rectangle (Triangulation<3>&                           tria,
          Assert (spacing[i][j] >= 0, ExcInvalidRepetitions(-1));
          delta = std::min (delta, spacing[i][j]);
        }
-      Assert(material_id.size(i) == repetitions[i], 
+      Assert(material_id.size(i) == repetitions[i],
             ExcInvalidRepetitionsDimension(i));
     }
 
@@ -1018,7 +1018,7 @@ subdivided_hyper_rectangle (Triangulation<3>&                           tria,
   SubCellData t;
   GridTools::delete_unused_vertices (points, cells, t);
 
-  tria.create_triangulation (points, cells, t);  
+  tria.create_triangulation (points, cells, t);
 
                                   // set boundary indicator
   if (colorize)
@@ -1077,14 +1077,14 @@ GridGenerator::colorize_subdivided_hyper_rectangle (Triangulation<dim> &tria,
                                                    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)
@@ -1110,7 +1110,7 @@ GridGenerator::colorize_subdivided_hyper_rectangle (Triangulation<dim> &tria,
                                             // but we could not find
                                             // on which boundary.
            Assert (false, ExcInternalError());
-         
+
        }
     }
   for (typename Triangulation<dim>::cell_iterator cell = tria.begin();
@@ -1188,6 +1188,18 @@ void GridGenerator::cylinder (Triangulation<dim> &,
 
 
 
+// Implementation for 1D only
+template <int dim>
+void GridGenerator::truncated_cone (Triangulation<dim> &,
+                                   const double,
+                                   const double,
+                                   const double)
+{
+  Assert (false, ExcNotImplemented());
+}
+
+
+
 // Implementation for 1D only
 template <int dim>
 void GridGenerator::hyper_shell (Triangulation<dim> &,
@@ -1276,7 +1288,7 @@ void GridGenerator::enclosed_hyper_cube (Triangulation<dim> &tria,
                                       1, 0, 2,
                                       9, 8,10
   };
-  
+
   std::vector<CellData<dim> > cells(9);
   k = 0;
   for (unsigned int i0=0;i0<3;++i0)
@@ -1342,6 +1354,45 @@ GridGenerator::hyper_cube_slit (Triangulation<dim> &tria,
 }
 
 
+
+template <int dim>
+void GridGenerator::truncated_cone (Triangulation<dim> &triangulation,
+                                   const double radius_0,
+                                   const double radius_1,
+                                   const double half_length)
+{
+  Point<dim> vertices_tmp[4];
+
+  vertices_tmp[0] = Point<dim> (-half_length, -radius_0);
+  vertices_tmp[1] = Point<dim> (half_length, -radius_1);
+  vertices_tmp[2] = Point<dim> (-half_length, radius_0);
+  vertices_tmp[3] = Point<dim> (half_length, radius_1);
+
+  const std::vector<Point<dim> > vertices (&vertices_tmp[0], &vertices_tmp[4]);
+  unsigned int cell_vertices[1][GeometryInfo<dim>::vertices_per_cell];
+
+  for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell; ++i)
+    cell_vertices[0][i] = i;
+
+  std::vector<CellData<dim> > cells (1, CellData<dim> ());
+
+  for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell; ++i)
+    cells[0].vertices[i] = cell_vertices[0][i];
+
+  cells[0].material_id = 0;
+  triangulation.create_triangulation (vertices, cells, SubCellData ());
+
+  typename Triangulation<dim>::cell_iterator cell = triangulation.begin ();
+
+  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 <int dim>
@@ -1363,14 +1414,14 @@ GridGenerator::hyper_L (Triangulation<dim> &tria,
                                   {3, 4, 6, 7}};
 
   std::vector<CellData<dim> > cells (3, CellData<dim>());
-  
-  for (unsigned int i=0; i<3; ++i) 
+
+  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<dim> >(&vertices[0], &vertices[8]),
     cells,
@@ -1398,7 +1449,7 @@ GridGenerator::hyper_ball (Triangulation<dim> &tria,
                                 p+Point<dim>(+1,+1)*(radius/std::sqrt(2.0)*a),
                                 p+Point<dim>(-1,+1)*(radius/std::sqrt(2.0)),
                                 p+Point<dim>(+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},
@@ -1406,14 +1457,14 @@ GridGenerator::hyper_ball (Triangulation<dim> &tria,
                                   {6, 4, 7, 5}};
 
   std::vector<CellData<dim> > cells (5, CellData<dim>());
-  
-  for (unsigned int i=0; i<5; ++i) 
+
+  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;
     };
-  
+
   tria.create_triangulation (
     std::vector<Point<dim> >(&vertices[0], &vertices[8]),
     cells,
@@ -1433,9 +1484,9 @@ void GridGenerator::hyper_shell (Triangulation<dim>   &tria,
 {
   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
@@ -1469,17 +1520,17 @@ void GridGenerator::hyper_shell (Triangulation<dim>   &tria,
     };
 
   std::vector<CellData<dim> > cells (N, CellData<dim>());
-       
-  for (unsigned int i=0; i<N; ++i) 
+
+  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);
-           
+
       cells[i].material_id = 0;
     };
-  
+
   tria.create_triangulation (
     vertices, cells, SubCellData());
 
@@ -1535,7 +1586,7 @@ GridGenerator::cylinder (Triangulation<dim> &tria,
                break;
          default:
                f->set_boundary_indicator(0);
-               break;      
+               break;
        }
       ++f;
     }
@@ -1574,21 +1625,21 @@ GridGenerator::half_hyper_ball (Triangulation<dim> &tria,
                                 p+Point<dim>(+1,+1)*(radius/std::sqrt(2.0)*a),
                                 p+Point<dim>(0,+1)*radius,
                                 p+Point<dim>(+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<dim> > cells (4, CellData<dim>());
-  
-  for (unsigned int i=0; i<4; ++i) 
+
+  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<dim> >(&vertices[0], &vertices[8]),
     cells,
@@ -1626,7 +1677,7 @@ GridGenerator::half_hyper_shell (Triangulation<dim>   &tria,
 {
   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
@@ -1670,17 +1721,17 @@ GridGenerator::half_hyper_shell (Triangulation<dim>   &tria,
 
 
   std::vector<CellData<dim> > cells (N, CellData<dim>());
-       
-  for (unsigned int i=0; i<N; ++i) 
+
+  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));
-           
+
       cells[i].material_id = 0;
     };
-  
+
   tria.create_triangulation (vertices, cells, SubCellData());
 }
 
@@ -1701,7 +1752,7 @@ void GridGenerator::hyper_cube_slit (Triangulation<dim>& tria,
 {
   const double rl2=(right+left)/2;
   const double len = (right-left)/2.;
-  
+
   const Point<dim> vertices[20] = {
        Point<dim>(left, left , -len/2.),
        Point<dim>(rl2,  left , -len/2.),
@@ -1784,7 +1835,7 @@ void GridGenerator::enclosed_hyper_cube (Triangulation<dim> &tria,
        33,32,34,
        41,40,42
   };
-  
+
   std::vector<CellData<dim> > cells(27);
   k = 0;
   for (unsigned int z=0;z<3;++z)
@@ -1811,6 +1862,80 @@ void GridGenerator::enclosed_hyper_cube (Triangulation<dim> &tria,
 
 
 
+template <int dim>
+void GridGenerator::truncated_cone (Triangulation<dim> &triangulation,
+                                   const double radius_0,
+                                   const double radius_1,
+                                   const double half_length)
+{
+  Assert (dim <= 3, ExcNotImplemented());
+                                  // Determine number of cells and vertices
+  const unsigned int n_cells = std::floor (half_length / std::max (radius_0, radius_1) + 0.5);
+  const unsigned int n_vertices = 4 * (n_cells + 1);
+  std::vector<Point<dim> > vertices_tmp(n_vertices);
+
+  vertices_tmp[0] = Point<dim> (-half_length, 0, -radius_0);
+  vertices_tmp[1] = Point<dim> (-half_length, radius_0, 0);
+  vertices_tmp[2] = Point<dim> (-half_length, -radius_0, 0);
+  vertices_tmp[3] = Point<dim> (-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<dim> (dx, 0, 0.5 * (radius_0 - radius_1) * dx / half_length);
+    vertices_tmp[4 * i + 5]
+      = vertices_tmp[4 * i + 1] +
+      Point<dim> (dx, 0.5 * (radius_1 - radius_0) * dx / half_length, 0);
+    vertices_tmp[4 * i + 6]
+      = vertices_tmp[4 * i + 2] +
+      Point<dim> (dx, 0.5 * (radius_0 - radius_1) * dx / half_length, 0);
+    vertices_tmp[4 * i + 7]
+      = vertices_tmp[4 * i + 3] +
+      Point<dim> (dx, 0, 0.5 * (radius_1 - radius_0) * dx / half_length);
+  }
+
+  const std::vector<Point<dim> > vertices (&vertices_tmp[0], &vertices_tmp[n_vertices]);
+  Table<2,unsigned int> cell_vertices(n_cells,GeometryInfo<dim>::vertices_per_cell);
+
+  for (unsigned int i = 0; i < n_cells; ++i)
+    for (unsigned int j = 0; j < GeometryInfo<dim>::vertices_per_cell; ++j)
+      cell_vertices[i][j] = 4 * i + j;
+
+  std::vector<CellData<dim> > cells (n_cells, CellData<dim> ());
+
+  for (unsigned int i = 0; i < n_cells; ++i) {
+    for (unsigned int j = 0; j < GeometryInfo<dim>::vertices_per_cell; ++j)
+      cells[i].vertices[j] = cell_vertices[i][j];
+
+    cells[i].material_id = 0;
+  }
+
+  triangulation.create_triangulation (vertices, cells, SubCellData ());
+
+  for (typename Triangulation<dim>::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);
+    }
+
+    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 = 0; i < 4; ++i)
+      cell->face (i)->set_boundary_indicator (0);
+  }
+}
+
+
 // Implementation for 3D only
 template <int dim>
 void
@@ -1862,8 +1987,8 @@ GridGenerator::hyper_L (Triangulation<dim> &tria,
                                   {12, 13, 21, 22, 15, 16, 24, 25}};
 
   std::vector<CellData<dim> > cells (7, CellData<dim>());
-  
-  for (unsigned int i=0; i<7; ++i) 
+
+  for (unsigned int i=0; i<7; ++i)
     {
       for (unsigned int j=0; j<8; ++j)
        cells[i].vertices[j] = cell_vertices[i][j];
@@ -1923,10 +2048,10 @@ GridGenerator::hyper_ball (Triangulation<dim> &tria,
                                         {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<dim> > cells (n_cells, CellData<dim>());
-  
-  for (unsigned int i=0; i<n_cells; ++i) 
+
+  for (unsigned int i=0; i<n_cells; ++i)
     {
       for (unsigned int j=0; j<GeometryInfo<dim>::vertices_per_cell; ++j)
        cells[i].vertices[j] = cell_vertices[i][j];
@@ -1949,7 +2074,7 @@ GridGenerator::cylinder (Triangulation<dim> &tria,
                         const double half_length)
 {
   Assert (dim <= 3, ExcNotImplemented());
-  
+
                                   // Copy the base from hyper_ball<dim>
                                   // and transform it to yz
   const double d = radius/std::sqrt(2.0);
@@ -1987,7 +2112,7 @@ GridGenerator::cylinder (Triangulation<dim> &tria,
       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},
@@ -1998,10 +2123,10 @@ GridGenerator::cylinder (Triangulation<dim> &tria,
   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<dim> > cells (10, CellData<dim>());
-  
-  for (unsigned int i=0; i<10; ++i) 
+
+  for (unsigned int i=0; i<10; ++i)
     {
       for (unsigned int j=0; j<8; ++j)
        cells[i].vertices[j] = cell_vertices[i][j];
@@ -2025,7 +2150,7 @@ GridGenerator::cylinder (Triangulation<dim> &tria,
                                   // above
   typename Triangulation<dim>::cell_iterator cell = tria.begin();
   typename Triangulation<dim>::cell_iterator end = tria.end();
-  
+
   for (; cell != end; ++cell)
     for (unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i)
       if (cell->at_boundary(i))
@@ -2105,7 +2230,7 @@ GridGenerator::half_hyper_ball (Triangulation<dim>& tria,
 
     std::vector<CellData<dim> > cells (6, CellData<dim>());
 
-    for (unsigned int i=0; i<6; ++i) 
+    for (unsigned int i=0; i<6; ++i)
     {
        for (unsigned int j=0; j<8; ++j)
            cells[i].vertices[j] = cell_vertices[i][j];
@@ -2131,7 +2256,7 @@ GridGenerator::half_hyper_ball (Triangulation<dim>& tria,
            // element
            if (cell->face(i)->center()(0) < center(0)+1.e-5) {
                cell->face(i)->set_boundary_indicator(1);
-               for (unsigned int j=0;j<GeometryInfo<dim>::lines_per_face;++j) 
+               for (unsigned int j=0;j<GeometryInfo<dim>::lines_per_face;++j)
                    cell->face(i)->line(j)->set_boundary_indicator(1);
            }
        }
@@ -2144,7 +2269,7 @@ GridGenerator::half_hyper_ball (Triangulation<dim>& tria,
            // If the center is not on the plane x=0, this is a curvilinear
            // element
            if (cell->face(i)->center()(0) > center(0)+1.e-5) {
-               for (unsigned int j=0;j<GeometryInfo<dim>::lines_per_face;++j) 
+               for (unsigned int j=0;j<GeometryInfo<dim>::lines_per_face;++j)
                    cell->face(i)->line(j)->set_boundary_indicator(0);
            }
        }
@@ -2168,7 +2293,7 @@ void GridGenerator::hyper_shell (Triangulation<dim>& tria,
   const double orad = outer_radius/std::sqrt(3.0);
   std::vector<Point<dim> > vertices;
   std::vector<CellData<dim> > cells;
-  
+
                                   // Start with the shell bounded by
                                   // two nested cubes
   if (n <= 6)
@@ -2177,7 +2302,7 @@ void GridGenerator::hyper_shell (Triangulation<dim>& tria,
        vertices.push_back(p+hexahedron[i]*irad);
       for (unsigned int i=0;i<8;++i)
        vertices.push_back(p+hexahedron[i]*orad);
-      
+
       const unsigned int n_cells = 6;
       const int cell_vertices[n_cells][8] =
        {{8, 9, 10, 11, 0, 1, 2, 3}, // bottom
@@ -2186,10 +2311,10 @@ void GridGenerator::hyper_shell (Triangulation<dim>& tria,
         {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
-      
+
       cells.resize(n_cells, CellData<dim>());
-      
-      for (unsigned int i=0; i<n_cells; ++i) 
+
+      for (unsigned int i=0; i<n_cells; ++i)
        {
          for (unsigned int j=0; j<GeometryInfo<dim>::vertices_per_cell; ++j)
            cells[i].vertices[j] = cell_vertices[i][j];
@@ -2224,10 +2349,10 @@ void GridGenerator::hyper_shell (Triangulation<dim>& tria,
         {  2,  8, 11,  6},
         { 12,  0,  2,  8},
         {  1, 10, 12,  0}};
-      
+
       cells.resize(n_cells, CellData<dim>());
-      
-      for (unsigned int i=0; i<n_cells; ++i) 
+
+      for (unsigned int i=0; i<n_cells; ++i)
        {
          for (unsigned int j=0; j<4; ++j)
            {
@@ -2241,11 +2366,11 @@ void GridGenerator::hyper_shell (Triangulation<dim>& tria,
     {
       Assert(false, ExcIndexRange(n, 1, 7));
     }
-  
+
   tria.create_triangulation (vertices, cells,
                             SubCellData());       // no boundary
                                                   // information
-  
+
   if (colorize)
     colorize_hyper_shell(tria, p, inner_radius, outer_radius);
 }
@@ -2271,7 +2396,7 @@ GridGenerator::colorize_hyper_shell (
       (++cell)->face(2)->set_boundary_indicator(1);
       (++cell)->face(0)->set_boundary_indicator(1);
       (++cell)->face(2)->set_boundary_indicator(1);
-      (++cell)->face(0)->set_boundary_indicator(1);      
+      (++cell)->face(0)->set_boundary_indicator(1);
     }
   else
                                     // For higher polyhedra, this is regular.
@@ -2337,7 +2462,7 @@ GridGenerator::half_hyper_shell (Triangulation<dim>& tria,
 
       std::vector<CellData<dim> > cells (5, CellData<dim>());
 
-      for (unsigned int i=0; i<5; ++i) 
+      for (unsigned int i=0; i<5; ++i)
        {
          for (unsigned int j=0; j<8; ++j)
            cells[i].vertices[j] = cell_vertices[i][j];
@@ -2369,9 +2494,9 @@ void GridGenerator::cylinder_shell (Triangulation<dim>   &tria,
 {
   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
@@ -2416,11 +2541,11 @@ void GridGenerator::cylinder_shell (Triangulation<dim>   &tria,
                           j*length/N_z);
         vertices_3d.push_back (v);
       }
-                            
+
   std::vector<CellData<dim> > cells (N_r*N_z, CellData<dim>());
-  
+
   for (unsigned int j=0; j<N_z; ++j)
-    for (unsigned int i=0; i<N_r; ++i) 
+    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;
@@ -2431,10 +2556,10 @@ void GridGenerator::cylinder_shell (Triangulation<dim>   &tria,
         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;
       }
-  
+
   tria.create_triangulation (
     vertices_3d, cells, SubCellData());
 }
@@ -2459,13 +2584,13 @@ void GridGenerator::laplace_solve (const SparseMatrix<double> &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);
-  
+
   Vector<double> f(n_dofs);
-  
+
   SF.add_constraints(m);
   SF.apply_constraints (f, true);
   solver.solve(SF, u, f, PF);
@@ -2491,21 +2616,21 @@ void GridGenerator::laplace_transformation (Triangulation<dim> &tria,
 {
                                   // first provide everything that is
                                   // needed for solving a Laplace
-                                  // equation.  
-  MappingQ1<dim> mapping_q1;  
+                                  // equation.
+  MappingQ1<dim> mapping_q1;
   FE_Q<dim> q1(1);
-  
+
   DoFHandler<dim> dof_handler(tria);
   dof_handler.distribute_dofs(q1);
   SparsityPattern sparsity_pattern (dof_handler.n_dofs (), dof_handler.n_dofs (),
                                    dof_handler.max_couplings_between_dofs());
   DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
   sparsity_pattern.compress ();
-  
+
   SparseMatrix<double> S(sparsity_pattern);
-  
+
   QGauss4<dim> quadrature;
-  
+
   MatrixCreator::create_laplace_matrix(mapping_q1, dof_handler, quadrature, S);
 
                                   // set up the boundary values for
@@ -2539,20 +2664,20 @@ void GridGenerator::laplace_transformation (Triangulation<dim> &tria,
                }
          }
     }
-  
+
                                   // 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
@@ -2589,7 +2714,7 @@ void GridGenerator::hyper_cube_with_cylindrical_hole (Triangulation<dim> &,
 #if deal_II_dimension == 2
 
 template<int dim>
-void GridGenerator::hyper_cube_with_cylindrical_hole (Triangulation<dim> &triangulation, 
+void GridGenerator::hyper_cube_with_cylindrical_hole (Triangulation<dim> &triangulation,
                                   const double inner_radius,
                                   const double outer_radius,
                                   const double, // width,
@@ -2632,7 +2757,7 @@ void GridGenerator::hyper_cube_with_cylindrical_hole (Triangulation<dim> &triang
            }
          }
        }
-      } 
+      }
   }
   double eps = 1e-3 * outer_radius;
   cell = triangulation.begin_active();
@@ -2650,7 +2775,7 @@ void GridGenerator::hyper_cube_with_cylindrical_hole (Triangulation<dim> &triang
            cell->face(f)->set_boundary_indicator(2);
          else if(std::abs(dy - outer_radius) < eps)
            cell->face(f)->set_boundary_indicator(3);
-         else 
+         else
            cell->face(f)->set_boundary_indicator(4);
        } else {
          double d = (cell->face(f)->center() - center).norm();
@@ -2660,7 +2785,7 @@ void GridGenerator::hyper_cube_with_cylindrical_hole (Triangulation<dim> &triang
            cell->face(f)->set_boundary_indicator(0);
        }
       }
-   }   
+   }
 }
 
 #endif
@@ -2668,7 +2793,7 @@ void GridGenerator::hyper_cube_with_cylindrical_hole (Triangulation<dim> &triang
 #if deal_II_dimension == 3
 
 template<int dim>
-void GridGenerator::hyper_cube_with_cylindrical_hole(Triangulation<dim> &triangulation, 
+void GridGenerator::hyper_cube_with_cylindrical_hole(Triangulation<dim> &triangulation,
                                                 const double inner_radius,
                                                 const double outer_radius,
                                                 const double L,
@@ -2677,15 +2802,15 @@ void GridGenerator::hyper_cube_with_cylindrical_hole(Triangulation<dim> &triangu
 {
   Assert(inner_radius < outer_radius,
         ExcMessage("outer_radius has to be bigger than inner_radius."));
-  Assert(L > 0, 
+  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);
-  
+
   typename Triangulation<dim>::active_cell_iterator
     cell = triangulation.begin_active(),
     endc = triangulation.end();
@@ -2718,7 +2843,7 @@ void GridGenerator::hyper_cube_with_cylindrical_hole(Triangulation<dim> &triangu
            }
          }
        }
-      } 
+      }
   }
   double eps = 1e-3 * outer_radius;
   cell = triangulation.begin_active();
@@ -2728,26 +2853,26 @@ void GridGenerator::hyper_cube_with_cylindrical_hole(Triangulation<dim> &triangu
        double dx = cell->face(f)->center()(0);
        double dy = cell->face(f)->center()(1);
        double dz = cell->face(f)->center()(2);
-       
+
        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 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)
@@ -2766,7 +2891,7 @@ void GridGenerator::hyper_cube_with_cylindrical_hole(Triangulation<dim> &triangu
            cell->face(f)->set_boundary_indicator(0);
        }
       }
-   }   
+   }
 }
 
 #endif
@@ -2789,14 +2914,14 @@ GridGenerator::hyper_rectangle<deal_II_dimension> (
   const bool);
 
 template void
-GridGenerator::subdivided_hyper_rectangle<deal_II_dimension> 
+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>&, bool);
 
 template void
-GridGenerator::subdivided_hyper_rectangle<deal_II_dimension> 
+GridGenerator::subdivided_hyper_rectangle<deal_II_dimension>
 (Triangulation<deal_II_dimension> &,
  const std::vector<std::vector<double> >&,
  const Point<deal_II_dimension>&,
@@ -2822,6 +2947,10 @@ GridGenerator::cylinder (
   Triangulation<deal_II_dimension>&, double, double);
 
 
+template void
+GridGenerator::truncated_cone (
+  Triangulation<deal_II_dimension>&, double, double, double);
+
 template void
 GridGenerator::hyper_L (
   Triangulation<deal_II_dimension>&, double, double);
@@ -2851,7 +2980,7 @@ GridGenerator::half_hyper_shell (
   const Point<deal_II_dimension>&, double, double, unsigned int);
 
 
-template void 
+template void
 GridGenerator::hyper_cube_with_cylindrical_hole (
   Triangulation<deal_II_dimension> &,
   const double, const double, const double, const unsigned int, bool);
index 803508b3b79f40a63617e819fbb2096f46c5f1d6..871bf7d7b19451cf8c6a76864c0d374e2550c0de 100644 (file)
@@ -96,7 +96,7 @@ CylinderBoundary<3>::
 get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const
 {
   const Point<3> middle = StraightBoundary<3>::get_new_point_on_quad (quad);
-  
+
                                   // same algorithm as above
   const unsigned int dim = 3;
   const Point<dim> vector_from_axis = (middle-point_on_axis) -
@@ -135,7 +135,7 @@ CylinderBoundary<dim>::get_intermediate_points_on_line (
     get_intermediate_points_between_points(line->vertex(0), line->vertex(1), points);
 }
 
-  
+
 template <int dim>
 void
 CylinderBoundary<dim>::get_intermediate_points_between_points (
@@ -181,18 +181,18 @@ CylinderBoundary<3>::get_intermediate_points_on_quad (
     {
       unsigned int m=static_cast<unsigned int> (std::sqrt(static_cast<double>(points.size())));
       Assert(points.size()==m*m, ExcInternalError());
-      
+
       std::vector<Point<3> > lp0(m);
       std::vector<Point<3> > lp1(m);
 
       get_intermediate_points_on_line(quad->line(0), lp0);
       get_intermediate_points_on_line(quad->line(1), lp1);
-      
+
       std::vector<Point<3> > lps(m);
       for (unsigned int i=0; i<m; ++i)
        {
          get_intermediate_points_between_points(lp0[i], lp1[i], lps);
-         
+
          for (unsigned int j=0; j<m; ++j)
            points[i*m+j]=lps[j];
        }
@@ -239,7 +239,7 @@ get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &face,
 
       const Point<dim> vector_from_axis = (vertex-point_on_axis) -
                                          ((vertex-point_on_axis) * direction) * direction;
-      
+
       face_vertex_normals[v] = (vector_from_axis / vector_from_axis.norm());
     }
 }
@@ -248,12 +248,257 @@ get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &face,
 
 template <int dim>
 double
-CylinderBoundary<dim>::get_radius () const 
+CylinderBoundary<dim>::get_radius () const
 {
   return radius;
 }
 
 
+//======================================================================//
+
+template<int dim>
+ConeBoundary<dim>::ConeBoundary (const double radius_0,
+                                const double radius_1,
+                                const Point<dim> x_0,
+                                const Point<dim> x_1)
+               :
+               radius_0 (radius_0),
+               radius_1 (radius_1),
+               x_0 (x_0),
+               x_1 (x_1)
+{}
+
+
+
+template<int dim>
+double ConeBoundary<dim>::get_radius (Point<dim> x) const
+{
+  for (unsigned int i = 0; i < dim; ++i)
+    if ((x_1 (i) - x_0 (i)) != 0)
+      return (radius_1 - radius_0) * (x (i) - x_0 (i)) / (x_1 (i) - x_0 (i)) + radius_0;
+
+  return 0;
+}
+
+
+
+template<int dim>
+void
+ConeBoundary<dim>::
+get_intermediate_points_between_points (const Point<dim> &p0,
+                                       const Point<dim> &p1,
+                                       std::vector<Point<dim> > &points) const
+{
+  const unsigned int n = points.size ();
+  const Point<dim> axis = x_1 - x_0;
+  const Point<dim> dx = (p1 - p0) / (n + 1);
+
+  Assert (n > 0, ExcInternalError ());
+
+  for (unsigned int i = 0; i < n; ++i)
+    {
+                                      // Compute the current point.
+      const Point<dim> x_i = p0 + (i + 1) * dx;
+                                      // To project this point on the
+                                      // boundary of the cone we first
+                                      // compute the orthogonal
+                                      // projection of this point onto
+                                      // the axis of the cone.
+      const double c = (x_i - x_0) * axis / axis.square ();
+      const Point<dim> x_ip = x_0 + c * axis;
+                                      // If the middle point is located
+                                      // on the axis hand over the
+                                      // current point itself, else
+                                      // hand over the projection of it
+                                      // on the boundary.
+      if (x_ip.norm () <= 1e-10 * x_i.norm ())
+       points[i] = x_i;
+      else {
+                                        // Compute the projection of
+                                        // the middle point on the
+                                        // boundary of the cone.
+       points[i] = x_ip + get_radius (x_ip) *  (x_i - x_ip) / (x_i - x_ip).norm ();
+      }
+    }
+}
+
+template<int dim>
+Point<dim>
+ConeBoundary<dim>::
+get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const
+{
+  const Point<dim> axis = x_1 - x_0;
+                                  // Compute the middle point of the line.
+  const Point<dim> middle = StraightBoundary<dim>::get_new_point_on_line (line);
+                                  // To project it on the boundary of
+                                  // the cone we first compute the
+                                  // orthogonal projection of the
+                                  // middle point onto the axis of
+                                  // the cone.
+  const double c = (middle - x_0) * axis / axis.square ();
+  const Point<dim> middle_p = x_0 + c * axis;
+                                  // If the middle point is located
+                                  // on the axis return the middle
+                                  // point itself, else return the
+                                  // projection of it on the
+                                  // boundary.
+  if (middle_p.norm () <= 1e-10 * middle.norm ())
+    return middle;
+  else {
+                                    // Compute the projection of the
+                                    // middle point on the boundary
+                                    // of the cone.
+    return middle_p + get_radius (middle_p) * (middle - middle_p) / (middle - middle_p).norm ();
+  }
+}
+
+
+
+template<>
+Point<3>
+ConeBoundary<3>::get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const
+{
+  const Point<3> axis = x_1 - x_0;
+                                  // Compute the middle point of the
+                                  // quad.
+  const Point<3> middle = StraightBoundary<3>::get_new_point_on_quad (quad);
+                                  // Same algorithm as above: To
+                                  // project it on the boundary of
+                                  // the cone we first compute the
+                                  // orthogonal projection of the
+                                  // middle point onto the axis of
+                                  // the cone.
+  const double c = (middle - x_0) * axis / axis.square ();
+  const Point<3> middle_p = x_0 + c * axis;
+                                  // If the middle point is located
+                                  // on the axis return the middle
+                                  // point itself, else return the
+                                  // projection of it on the
+                                  // boundary.
+  if (middle_p.norm () <= 1e-10 * middle.norm ())
+    return middle;
+  else {
+                                    // Compute the projection of the
+                                    // middle point on the boundary
+                                    // of the cone.
+    return middle_p + get_radius (middle_p) * (middle - middle_p) / (middle - middle_p).norm ();
+  }
+}
+
+
+
+template<int dim>
+Point<dim>
+ConeBoundary<dim>::
+get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &) const
+{
+  Assert (false, ExcImpossibleInDim (dim));
+
+  return Point<dim>();
+}
+
+
+
+template<int dim>
+void
+ConeBoundary<dim>::
+get_intermediate_points_on_line (const typename Triangulation<dim>::line_iterator &line,
+                                std::vector<Point<dim> > &points) const
+{
+  if (points.size () == 1)
+    points[0] = get_new_point_on_line (line);
+  else
+    get_intermediate_points_between_points (line->vertex (0), line->vertex (1), points);
+}
+
+
+#if deal_II_dimension == 3
+
+template<>
+void
+ConeBoundary<3>::
+get_intermediate_points_on_quad (const Triangulation<3>::quad_iterator &quad,
+                                std::vector<Point<3> > &points) const
+{
+  if (points.size () == 1)
+    points[0] = get_new_point_on_quad (quad);
+  else
+    {
+      unsigned int n = static_cast<unsigned int> (std::sqrt (static_cast<double> (points.size ())));
+
+      Assert (points.size () == n * n, ExcInternalError ());
+
+      std::vector<Point<3> > points_line_0 (n);
+      std::vector<Point<3> > points_line_1 (n);
+
+      get_intermediate_points_on_line (quad->line (0), points_line_0);
+      get_intermediate_points_on_line (quad->line (1), points_line_1);
+
+      std::vector<Point<3> > points_line_segment (n);
+
+      for (unsigned int i = 0; i < n; ++i)
+       {
+         get_intermediate_points_between_points (points_line_0[i],
+                                                 points_line_1[i],
+                                                 points_line_segment);
+
+         for (unsigned int j = 0; j < n; ++j)
+           points[i * n + j] = points_line_segment[j];
+       }
+    }
+}
+
+#endif
+
+template <int dim>
+void
+ConeBoundary<dim>::
+get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterator &,
+                                std::vector<Point<dim> > &) const
+{
+  Assert (false, ExcImpossibleInDim (dim));
+}
+
+
+#if deal_II_dimension == 1
+
+template<>
+void
+ConeBoundary<1>::
+get_normals_at_vertices (const Triangulation<1>::face_iterator &,
+                        Boundary<1>::FaceVertexNormals &) const
+{
+  Assert (false, ExcImpossibleInDim (1));
+}
+
+#endif
+
+template<int dim>
+void
+ConeBoundary<dim>::
+get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &face,
+                        typename Boundary<dim>::FaceVertexNormals &face_vertex_normals) const
+{
+  const Point<dim> axis = x_1 - x_0;
+
+  for (unsigned int vertex = 0; vertex < GeometryInfo<dim>::vertices_per_cell; ++vertex)
+    {
+                                      // Compute the orthogonal
+                                      // projection of the vertex onto
+                                      // the axis of the cone.
+      const double c = (face->vertex (vertex) - x_0) * axis / axis.square ();
+      const Point<dim> vertex_p = x_0 + c * axis;
+                                      // Then compute the vector
+                                      // pointing from the point
+                                      // <tt>vertex_p</tt> on the axis
+                                      // to the vertex.
+      const Point<dim> axis_to_vertex = face->vertex (vertex) - vertex_p;
+
+      face_vertex_normals[vertex] = axis_to_vertex / axis_to_vertex.norm ();
+    }
+}
+
+
 //======================================================================//
 
 template <int dim, int spacedim>
@@ -269,7 +514,7 @@ Point<spacedim>
 HyperBallBoundary<dim,spacedim>::get_new_point_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line) const
 {
   Point<spacedim> middle = StraightBoundary<dim,spacedim>::get_new_point_on_line (line);
-  
+
   middle -= center;
 
   double r=0;
@@ -282,7 +527,7 @@ HyperBallBoundary<dim,spacedim>::get_new_point_on_line (const typename Triangula
     r=radius;
                                   // project to boundary
   middle *= r / std::sqrt(middle.square());
-  
+
   middle += center;
   return middle;
 }
@@ -319,7 +564,7 @@ HyperBallBoundary<dim,spacedim>::
 get_new_point_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad) const
 {
   Point<spacedim> middle = StraightBoundary<dim,spacedim>::get_new_point_on_quad (quad);
-  
+
   middle -= center;
 
   double r=0;
@@ -332,7 +577,7 @@ get_new_point_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator
     r=radius;
                                   // project to boundary
   middle *= r / std::sqrt(middle.square());
-  
+
   middle += center;
   return middle;
 }
@@ -364,7 +609,7 @@ HyperBallBoundary<dim,spacedim>::get_intermediate_points_on_line (
     get_intermediate_points_between_points(line->vertex(0), line->vertex(1), points);
 }
 
-  
+
 template <int dim, int spacedim>
 void
 HyperBallBoundary<dim,spacedim>::get_intermediate_points_between_points (
@@ -377,7 +622,7 @@ HyperBallBoundary<dim,spacedim>::get_intermediate_points_between_points (
   const Point<spacedim> v0=p0-center,
                   v1=p1-center;
   const double length=std::sqrt((v1-v0).square());
-  
+
   double eps=1e-12;
   double r=0;
   if (compute_radius_automatically)
@@ -388,17 +633,17 @@ HyperBallBoundary<dim,spacedim>::get_intermediate_points_between_points (
   else
     r=radius;
 
-  
+
   const double r2=r*r;
   Assert(std::fabs(v0.square()-r2)<eps*r2, ExcInternalError());
   Assert(std::fabs(v1.square()-r2)<eps*r2, ExcInternalError());
-  
+
   const double alpha=std::acos((v0*v1)/std::sqrt(v0.square()*v1.square()));
   const double d_alpha=alpha/(n+1);
   const Point<spacedim> pm=0.5*(v0+v1);
-  
+
   const double h=std::sqrt(pm.square());
-  
+
   double beta=0;
   unsigned int left_index=0, right_index=0;
   if ((n+1)%2==0)
@@ -419,7 +664,7 @@ HyperBallBoundary<dim,spacedim>::get_intermediate_points_between_points (
       right_index=n/2;
       left_index=n/2-1;
     }
-      
+
                                   // n even:  m=n/2,
                                   // n odd:   m=(n-1)/2
   const unsigned int m=n/2;
@@ -427,12 +672,12 @@ HyperBallBoundary<dim,spacedim>::get_intermediate_points_between_points (
     {
       Assert(right_index<n, ExcInternalError());
       Assert(left_index+1>=1, ExcInternalError());
-         
+
       double d=h*std::tan(beta);
       points[right_index]=pm+d/length*(v1-v0);
       points[left_index]=pm-d/length*(v1-v0);
     }
-      
+
 
                                   // project the points from the
                                   // straight line to the
@@ -458,18 +703,18 @@ HyperBallBoundary<3>::get_intermediate_points_on_quad (
     {
       unsigned int m=static_cast<unsigned int> (std::sqrt(static_cast<double>(points.size())));
       Assert(points.size()==m*m, ExcInternalError());
-      
+
       std::vector<Point<3> > lp0(m);
       std::vector<Point<3> > lp1(m);
-      
+
       get_intermediate_points_on_line(quad->line(0), lp0);
       get_intermediate_points_on_line(quad->line(1), lp1);
-      
+
       std::vector<Point<3> > lps(m);
       for (unsigned int i=0; i<m; ++i)
        {
          get_intermediate_points_between_points(lp0[i], lp1[i], lps);
-         
+
          for (unsigned int j=0; j<m; ++j)
            points[i*m+j]=lps[j];
        }
@@ -527,7 +772,7 @@ get_normals_at_vertices (const typename Triangulation<dim,spacedim>::face_iterat
 
 template <int dim, int spacedim>
 Point<spacedim>
-HyperBallBoundary<dim,spacedim>::get_center () const 
+HyperBallBoundary<dim,spacedim>::get_center () const
 {
   return center;
 }
@@ -536,7 +781,7 @@ HyperBallBoundary<dim,spacedim>::get_center () const
 
 template <int dim, int spacedim>
 double
-HyperBallBoundary<dim,spacedim>::get_radius () const 
+HyperBallBoundary<dim,spacedim>::get_radius () const
 {
   Assert(!compute_radius_automatically, ExcRadiusNotSet());
   return radius;
@@ -557,7 +802,7 @@ HalfHyperBallBoundary<dim>::HalfHyperBallBoundary (const Point<dim> center,
 template <int dim>
 Point<dim>
 HalfHyperBallBoundary<dim>::
-get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const 
+get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const
 {
                                   // check whether center of object is
                                   // at x==0, since then it belongs
@@ -577,7 +822,7 @@ get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) c
 template <>
 Point<1>
 HalfHyperBallBoundary<1>::
-get_new_point_on_quad (const Triangulation<1>::quad_iterator &) const 
+get_new_point_on_quad (const Triangulation<1>::quad_iterator &) const
 {
   Assert (false, ExcInternalError());
   return Point<1>();
@@ -590,7 +835,7 @@ get_new_point_on_quad (const Triangulation<1>::quad_iterator &) const
 template <int dim>
 Point<dim>
 HalfHyperBallBoundary<dim>::
-get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const 
+get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const
 {
   const Point<dim> quad_center = quad->center();
   if (quad_center(0) == this->center(0))
@@ -728,7 +973,7 @@ HalfHyperShellBoundary<dim>::HalfHyperShellBoundary (const Point<dim> &center,
 template <int dim>
 Point<dim>
 HalfHyperShellBoundary<dim>::
-get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const 
+get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const
 {
   switch (dim)
     {
@@ -748,7 +993,7 @@ get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) c
                                           // as in the base class
          return HyperShellBoundary<dim>::get_new_point_on_line (line);
       }
-      
+
                                             // in 3d, a line is a straight
                                             // line if it is on the symmetry
                                             // plane and if not both of its
@@ -756,7 +1001,7 @@ get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) c
                                             // inner or outer sphere
       case 3:
       {
-       
+
            if (((line->vertex(0)(0) == this->center(0))
                 &&
                 (line->vertex(1)(0) == this->center(0)))
@@ -783,7 +1028,7 @@ get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) c
       default:
            Assert (false, ExcNotImplemented());
     }
-  
+
   return Point<dim>();
 }
 
@@ -822,7 +1067,7 @@ get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) c
                                      quad->vertex(2) + quad->vertex(3)   )/4;
       const Point<dim> quad_center_offset = quad_center - this->center;
 
-      
+
       if (std::fabs (quad->line(0)->center().distance(this->center) -
                     quad->line(1)->center().distance(this->center))
          < 1e-12 * outer_radius)
@@ -830,7 +1075,7 @@ get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) c
                                           // lines 0 and 1 are radial
          const double needed_radius
            = quad->line(0)->center().distance(this->center);
-         
+
          return (this->center +
                  quad_center_offset/quad_center_offset.norm() * needed_radius);
        }
@@ -841,14 +1086,14 @@ get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) c
                                           // lines 2 and 3 are radial
          const double needed_radius
            = quad->line(2)->center().distance(this->center);
-         
+
          return (this->center +
                  quad_center_offset/quad_center_offset.norm() * needed_radius);
        }
       else
        Assert (false, ExcInternalError());
     }
-  
+
                                   // otherwise we are on the outer or
                                   // inner part of the shell. proceed
                                   // as in the base class
@@ -881,7 +1126,7 @@ get_intermediate_points_on_line (const typename Triangulation<dim>::line_iterato
                                           // as in the base class
          HyperShellBoundary<dim>::get_intermediate_points_on_line (line, points);
       }
-      
+
                                             // in 3d, a line is a straight
                                             // line if it is on the symmetry
                                             // plane and if not both of its
@@ -889,7 +1134,7 @@ get_intermediate_points_on_line (const typename Triangulation<dim>::line_iterato
                                             // inner or outer sphere
       case 3:
       {
-       
+
            if (((line->vertex(0)(0) == this->center(0))
                 &&
                 (line->vertex(1)(0) == this->center(0)))
@@ -927,7 +1172,7 @@ get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterato
                                 std::vector<Point<dim> > &points) const
 {
   Assert (dim < 3, ExcNotImplemented());
-  
+
                                   // check whether center of object is
                                   // at x==0, since then it belongs
                                   // to the plane part of the
@@ -986,8 +1231,9 @@ get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &face,
 
 
 // explicit instantiations
-template class HyperBallBoundary<deal_II_dimension>;
 template class CylinderBoundary<deal_II_dimension>;
+template class ConeBoundary<deal_II_dimension>;
+template class HyperBallBoundary<deal_II_dimension>;
 template class HalfHyperBallBoundary<deal_II_dimension>;
 template class HyperShellBoundary<deal_II_dimension>;
 template class HalfHyperShellBoundary<deal_II_dimension>;
index 8a817b40c25fc1e84b537279b7a382ace90ad805..ffc62c805f30080273a5f125c6bbfff584517795 100644 (file)
@@ -60,6 +60,9 @@ contact with them, send email to <tt>authors</tt> at <tt>dealii.org</tt>.
 <li><em>Michael Anderson:</em>
     Linear complexity grid reordering algorithm in 2d and 3d.
 
+<li><em>Markus B&uuml;rg:</em>
+    Conical meshes and geometries.
+
 <li><em>John Burnell:</em>
     Configuration on Microsoft Windows systems.
 
@@ -81,7 +84,7 @@ contact with them, send email to <tt>authors</tt> at <tt>dealii.org</tt>.
     cubit journal file to export to ucd mesh format; &nbsp;
     FEFieldFunction and ParsedFunction classes; &nbsp;
     work on the codimension-one meshes, DoFHandler, and finite
-    elements; &nbsp; 
+    elements; &nbsp;
     singular integration; &nbsp;
     Step-34 tutorial program; &nbsp;
     random bug fixes and enhancements.
@@ -170,7 +173,7 @@ contact with them, send email to <tt>authors</tt> at <tt>dealii.org</tt>.
 
 <li><em>Yaqi Wang:</em>
     The step-28 tutorial program;
-    &nbsp; 
+    &nbsp;
     some grid generation functions.
 
 <li><em>Joshua White:</em>
index b5944929ba955639811227c8908c9457549dbd0d..e9e97e08a0fd59cdd550a22ac21ea56002bdb545 100644 (file)
@@ -359,6 +359,15 @@ inconvenience this causes.
 <h3>deal.II</h3>
 
 <ol>
+  <li>
+  <p>
+  New: The function GridGenerator::truncated_cone() and the class ConeBoundary
+  can now be used to describe conical objects.
+  <br>
+  (Markus B&uuml;rg 2009/08/17)
+  </p>
+  </li>
+
   <li>
   <p>
   New: Instead of asking face or edge iterators for their boundary indicator

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.