]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implement HalfHyperShellBoundary and the corresponding grid generation function.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 8 Aug 2000 11:42:59 +0000 (11:42 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 8 Aug 2000 11:42:59 +0000 (11:42 +0000)
git-svn-id: https://svn.dealii.org/trunk@3234 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

index f37d336a49f158b0fffe53d82d99203120b23da4..e46d51dbaae38d0c134475e401f1a5cad8ad1f9d 100644 (file)
@@ -41,7 +41,7 @@ template <int dim> class Triangulation;
  *      $[-1,1]^2 - [0,1]^2$. To create a hyper-L in one dimension results in
  *      an error. The function is also implemented for three space dimensions.
  *
- *    @item Hyper balls:
+ *    @item Hyper ball:
  *      You get the circle or ball (or generalized: hyperball) around origin
  *      @p{p} and with radius @p{r} by calling
  *      @ref{GridGenerator}@p{::hyper_ball (tria, p, r)}. The circle is triangulated
@@ -75,7 +75,13 @@ template <int dim> class Triangulation;
  *      suitable boundary class is provided as @ref{HyperSphereBoundary}
  *      in the library.
  *
- * @item Slit domain: The slit domain is a variant of the hyper cube
+ *    @item Half shells: This is a variant of the preceeding domain,
+ *      where only one half of the shell is taken, defined by the first
+ *      coordinate having non-negative values. This class is useful in
+ *      computations with rotational symmetry, as then a 3d shell can be
+ *      modeled by a 2d half shell.
+ *
+ *    @item Slit domain: The slit domain is a variant of the hyper cube
  *      domain. In two spatial dimensions, it is a square into which a slit
  *      is sawed; if the initial square is though to be composed of four
  *      smaller squares, then two of them are not connected even though
@@ -221,6 +227,49 @@ class GridGenerator
                             const double        outer_radius,
                             const unsigned int  n_cells = 0);
 
+
+                                    /**
+                                     * Produce a half hyper-shell,
+                                     * i.e. the space between two
+                                     * circles in two space
+                                     * dimensions and the region
+                                     * between two spheres in 3d,
+                                     * with given inner and outer
+                                     * radius and a given number of
+                                     * elements for this initial
+                                     * triangulation.  However,
+                                     * opposed to the previous
+                                     * function, it does not produce
+                                     * a whole shell, but only one
+                                     * half of it, namely that part
+                                     * for which the first component
+                                     * is restricted to non-negative
+                                     * values. The purpose of this
+                                     * class is to enable
+                                     * computations for solutions
+                                     * which have rotational
+                                     * symmetry, in which case the
+                                     * half shell in 2d represents a
+                                     * shell in 3d.
+                                     *
+                                     * If the number of
+                                     * initial cells is zero (as is
+                                     * the default), then it is
+                                     * computed adaptively such that
+                                     * the resulting elements have
+                                     * the least aspect ratio.
+                                     *
+                                     * The triangulation needs to be
+                                     * void upon calling this
+                                     * function.
+                                     */
+    template <int dim>
+    static void half_hyper_shell (Triangulation<dim> &tria,
+                                 const Point<dim>   &center,
+                                 const double        inner_radius,
+                                 const double        outer_radius,
+                                 const unsigned int  n_cells = 0);
+    
                                     /**
                                      * Exception
                                      */
index 70a8695cd61458f6a10764fc03df111853180e36..501272620b6c994e5b5e1e4e18e7d32fc12997b9 100644 (file)
@@ -33,7 +33,8 @@
  * @author Wolfgang Bangerth, 1998
  */
 template <int dim>
-class HyperBallBoundary : public StraightBoundary<dim> {
+class HyperBallBoundary : public StraightBoundary<dim>
+{
   public:
                                     /**
                                      * Constructor
@@ -80,16 +81,20 @@ class HyperBallBoundary : public StraightBoundary<dim> {
 };
 
 
+
 /**
  * Variant of @p{HyperBallBoundary} which denotes a half hyper ball
  * where the first coordinate is restricted to the range $x>=0$ (or
  * $x>=center(0)$). In two dimensions, this equals the right half
- * circle, in three space dimensions it is a half ball.
+ * circle, in three space dimensions it is a half ball. This class
+ * might be useful for computations with rotational symmetry, where
+ * one dimension is the radius from the axis of rotation.
  *
  * @author Wolfgang Bangerth, 1999
  */
 template <int dim>
-class HalfHyperBallBoundary : public HyperBallBoundary<dim> {
+class HalfHyperBallBoundary : public HyperBallBoundary<dim>
+{
   public:
                                     /**
                                      * Constructor
@@ -115,6 +120,7 @@ class HalfHyperBallBoundary : public HyperBallBoundary<dim> {
 };
 
 
+
 /**
  * Class describing the boundaries of a hyper shell. Only the center
  * of the two spheres needs to be given, the radii of inner and outer
@@ -154,4 +160,47 @@ class HyperShellBoundary : public StraightBoundary<dim>
 };
 
 
+
+/**
+ * Variant of @p{HyperShellBoundary} which denotes a half hyper shell
+ * where the first coordinate is restricted to the range $x>=0$ (or
+ * $x>=center(0)$). In two dimensions, this equals the right half arc,
+ * in three space dimensions it is a half shell. This class might be
+ * useful for computations with rotational symmetry, where one
+ * dimension is the radius from the axis of rotation.
+ *
+ * @author Wolfgang Bangerth, 2000
+ */
+template <int dim>
+class HalfHyperShellBoundary : public HyperShellBoundary<dim> 
+{
+  public:
+                                    /**
+                                     * Constructor. The center of the
+                                     * spheres defaults to the
+                                     * origin.
+                                     */
+    HalfHyperShellBoundary (const Point<dim> &center = Point<dim>());
+    
+                                    /**
+                                     * Construct a new point on a line.
+                                     */
+    virtual Point<dim>
+    get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const;  
+    
+                                    /**
+                                     * Construct a new point on a quad.
+                                     */
+    virtual Point<dim>
+    get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const;
+    
+  private:
+                                    /**
+                                     * Store the center of the spheres.
+                                     */
+    const Point<dim> center;
+};
+
+
+
 #endif
index b44f8eadb0bc4ca3273f19525fbc5532c441155c..8356e2f4a280e83ce31255c4e967624591738c7f 100644 (file)
@@ -35,6 +35,7 @@ void GridGenerator::hyper_cube<> (Triangulation<1> &tria,
 };
 
 
+
 template <>
 void GridGenerator::hyper_cube_slit<> (Triangulation<1> &,
                                       const double,
@@ -43,6 +44,7 @@ void GridGenerator::hyper_cube_slit<> (Triangulation<1> &,
 };
 
 
+
 template <>
 void GridGenerator::hyper_L<> (Triangulation<1> &,
                               const double,
@@ -51,6 +53,7 @@ void GridGenerator::hyper_L<> (Triangulation<1> &,
 };
 
 
+
 template <>
 void GridGenerator::hyper_ball<> (Triangulation<1> &,
                                  const Point<1> &,
@@ -59,6 +62,7 @@ void GridGenerator::hyper_ball<> (Triangulation<1> &,
 };
 
 
+
 template <>
 void GridGenerator::hyper_shell<> (Triangulation<1> &,
                                   const Point<1> &,
@@ -71,6 +75,7 @@ void GridGenerator::hyper_shell<> (Triangulation<1> &,
 #endif
 
 
+
 #if deal_II_dimension == 2
 
 template <>
@@ -136,7 +141,9 @@ void GridGenerator::enclosed_hyper_cube (Triangulation<2> &tria,
                             cells,
                             SubCellData());       // no boundary information
 }
-  
+
+
+
 template <>
 void GridGenerator::hyper_cube_slit<> (Triangulation<2> &tria,
                                       const double left,
@@ -169,6 +176,7 @@ void GridGenerator::hyper_cube_slit<> (Triangulation<2> &tria,
 };
 
 
+
 template <>
 void GridGenerator::hyper_L<> (Triangulation<2> &tria,
                               const double a,
@@ -201,6 +209,7 @@ void GridGenerator::hyper_L<> (Triangulation<2> &tria,
 };
 
 
+
 template <>
 void GridGenerator::hyper_ball<> (Triangulation<2> &tria,
                                  const Point<2>    &p,
@@ -238,6 +247,7 @@ void GridGenerator::hyper_ball<> (Triangulation<2> &tria,
 };
 
 
+
 template <>
 void GridGenerator::hyper_shell<> (Triangulation<2>   &tria,
                                   const Point<2>     &center,
@@ -297,6 +307,76 @@ void GridGenerator::hyper_shell<> (Triangulation<2>   &tria,
 };
 
 
+
+template <>
+void GridGenerator::half_hyper_shell<> (Triangulation<2>   &tria,
+                                       const Point<2>     &center,
+                                       const double        inner_radius,
+                                       const double        outer_radius,
+                                       const unsigned int  n_cells)
+{
+  Assert ((inner_radius > 0) && (inner_radius < outer_radius),
+         ExcInvalidRadii ());
+  
+  const double pi     = 3.14159265359;
+                                  // determine the number of cells
+                                  // for the grid. if not provided by
+                                  // the user determine it such that
+                                  // the length of each cell on the
+                                  // median (in the middle between
+                                  // the two circles) is equal to its
+                                  // radial extent (which is the
+                                  // difference between the two
+                                  // radii)
+  const unsigned int N = (n_cells == 0 ?
+                         static_cast<unsigned int>
+                         (ceil((pi* (outer_radius + inner_radius)/2) /
+                               (outer_radius - inner_radius))) :
+                         n_cells);
+
+                                  // set up N+1 vertices on the
+                                  // outer and N+1 vertices on
+                                  // the inner circle. the
+                                  // first N+1 ones are on the
+                                  // outer one, and all are
+                                  // numbered counter-clockwise
+  vector<Point<2> > vertices(2*(N+1));
+  for (unsigned int i=0; i<=N; ++i)
+    {
+                                  // enforce that the x-coordinates
+                                  // of the first and last point of
+                                  // each half-circle are exactly
+                                  // zero (contrary to what we may
+                                  // compute using the imprecise
+                                  // value of pi)
+      vertices[i] =  Point<2>( ( (i==0) || (i==N) ?
+                                0 :
+                                cos(pi * i/N - pi/2) ),
+                              sin(pi * i/N - pi/2)) * outer_radius;
+      vertices[i+N+1] = vertices[i] * (inner_radius/outer_radius);
+
+      vertices[i]     += center;
+      vertices[i+N+1] += center;
+    };
+
+
+  vector<CellData<2> > cells (N, CellData<2>());
+       
+  for (unsigned int i=0; i<N; ++i) 
+    {
+      cells[i].vertices[0] = i;
+      cells[i].vertices[1] = (i+1)%(N+1);
+      cells[i].vertices[2] = N+1+((i+1)%(N+1));
+      cells[i].vertices[3] = N+1+i;
+           
+      cells[i].material_id = 0;
+    };
+  
+  tria.create_triangulation (vertices, cells, SubCellData());
+};
+
+
+
 #endif
 
 
@@ -327,6 +407,7 @@ void GridGenerator::hyper_cube<> (Triangulation<3> &tria,
 };
 
 
+
 template <>
 void GridGenerator::hyper_cube_slit<> (Triangulation<3> &,
                                       const double,
@@ -335,6 +416,7 @@ void GridGenerator::hyper_cube_slit<> (Triangulation<3> &,
 };
 
 
+
 template<>
 void GridGenerator::enclosed_hyper_cube (Triangulation<3> &tria,
                                         const double l,
index 945d4c6fab2e7264fbe68a0a16334b5b2d1b84d0..67eff897c54f62a4a872ee87485c4c7506227a17 100644 (file)
@@ -26,6 +26,7 @@ HyperBallBoundary<dim>::HyperBallBoundary (const Point<dim> p,
 {};
 
 
+
 template <int dim>
 Point<dim>
 HyperBallBoundary<dim>::get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const
@@ -41,6 +42,7 @@ HyperBallBoundary<dim>::get_new_point_on_line (const typename Triangulation<dim>
 };
 
 
+
 #if deal_II_dimension == 1
 
 template <>
@@ -55,6 +57,7 @@ get_new_point_on_quad (const Triangulation<1>::quad_iterator &) const
 #endif
 
 
+
 template <int dim>
 Point<dim>
 HyperBallBoundary<dim>::
@@ -71,6 +74,7 @@ get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) c
 };
 
 
+
 template <int dim>
 Point<dim>
 HyperBallBoundary<dim>::get_center () const 
@@ -79,6 +83,7 @@ HyperBallBoundary<dim>::get_center () const
 };
 
 
+
 template <int dim>
 double
 HyperBallBoundary<dim>::get_radius () const 
@@ -87,6 +92,9 @@ HyperBallBoundary<dim>::get_radius () const
 };
 
 
+/* ---------------------------------------------------------------------- */
+
+
 template <int dim>
 HalfHyperBallBoundary<dim>::HalfHyperBallBoundary (const Point<dim> center,
                                                   const double     radius) :
@@ -94,6 +102,7 @@ HalfHyperBallBoundary<dim>::HalfHyperBallBoundary (const Point<dim> center,
 {};
 
 
+
 template <int dim>
 Point<dim>
 HalfHyperBallBoundary<dim>::
@@ -107,6 +116,7 @@ get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) c
 };
 
 
+
 #if deal_II_dimension == 1
 
 template <>
@@ -121,6 +131,7 @@ get_new_point_on_quad (const Triangulation<1>::quad_iterator &) const
 #endif
 
 
+
 template <int dim>
 Point<dim>
 HalfHyperBallBoundary<dim>::
@@ -134,12 +145,18 @@ get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) c
 };
 
 
+
+/* ---------------------------------------------------------------------- */
+
+
+
 template <int dim>
 HyperShellBoundary<dim>::HyperShellBoundary (const Point<dim> &center) :
                center (center) 
 {};
 
 
+
 template <int dim>
 Point<dim>
 HyperShellBoundary<dim>::
@@ -161,6 +178,7 @@ get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) c
 };
 
 
+
 #if deal_II_dimension == 1
 
 template <>
@@ -175,6 +193,7 @@ get_new_point_on_quad (const Triangulation<1>::quad_iterator &) const
 #endif
 
 
+
 template <int dim>
 Point<dim>
 HyperShellBoundary<dim>::
@@ -196,7 +215,80 @@ get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) c
 };
 
 
+/* ---------------------------------------------------------------------- */
+
+
+
+
+template <int dim>
+HalfHyperShellBoundary<dim>::HalfHyperShellBoundary (const Point<dim> &center) :
+               HyperShellBoundary<dim> (center) 
+{};
+
+
+
+template <int dim>
+Point<dim>
+HalfHyperShellBoundary<dim>::
+get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const 
+{
+                                  // first check whether the two end
+                                  // points of the line are on the
+                                  // axis of symmetry. if so, then
+                                  // return the mid point
+  if ((line->vertex(0)(0) == center(0)) &&
+      (line->vertex(1)(0) == center(0)))
+    return (line->vertex(0) + line->vertex(1))/2;
+  
+
+                                  // otherwise we are on the outer or
+                                  // inner part of the shell. proceed
+                                  // as in the base class
+  return HyperShellBoundary<dim>::get_new_point_on_line (line);
+};
+
+
+
+#if deal_II_dimension == 1
+
+template <>
+Point<1>
+HalfHyperShellBoundary<1>::
+get_new_point_on_quad (const Triangulation<1>::quad_iterator &) const
+{
+  Assert (false, ExcInternalError());
+  return Point<1>();
+};
+
+#endif
+
+
+
+template <int dim>
+Point<dim>
+HalfHyperShellBoundary<dim>::
+get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const
+{
+                                  // same thing as for the new point
+                                  // on the line
+  if ((quad->vertex(0)(0) == center(0)) &&
+      (quad->vertex(1)(0) == center(0)) &&
+      (quad->vertex(2)(0) == center(0)) &&
+      (quad->vertex(3)(0) == center(0)))
+    return (quad->vertex(0) + quad->vertex(1) +
+           quad->vertex(2) + quad->vertex(3)   )/4;
+  
+
+                                  // otherwise we are on the outer or
+                                  // inner part of the shell. proceed
+                                  // as in the base class
+  return HyperShellBoundary<dim>::get_new_point_on_quad (quad);
+};
+
+
+
 // explicit instantiations
 template class HyperBallBoundary<deal_II_dimension>;
 template class HalfHyperBallBoundary<deal_II_dimension>;
 template class HyperShellBoundary<deal_II_dimension>;
+template class HalfHyperShellBoundary<deal_II_dimension>;

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.