]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Cylinder domains
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 7 Dec 2001 14:44:22 +0000 (14:44 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 7 Dec 2001 14:44:22 +0000 (14:44 +0000)
git-svn-id: https://svn.dealii.org/trunk@5326 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/grid_out.cc
deal.II/deal.II/source/grid/tria_boundary_lib.cc

index 1b010905c5a6379ef2df3bdb22aad06c7525168f..0e45017fc1115b8781cc188757e9c84ac0f033fb 100644 (file)
@@ -251,6 +251,47 @@ class GridGenerator
                            const Point<3>   &center = Point<3>(),
                            const double      radius = 1.);
 
+                                    /**
+                                     * Create a cylinder around the
+                                     * x-axis.  The cylinder extends
+                                     * from @p{x=-half_length} to
+                                     * @p{x=+half_length} and its
+                                     * projection into the
+                                     * @p{yz}-plane is a circle of
+                                     * radius @p{radius}.
+                                     *
+                                     * The boundaries are colored
+                                     * according to the following
+                                     * scheme: 0 for the hull of the
+                                     * cylinder, 1 for the left hand
+                                     * face and 2 for the right hand
+                                     * face.
+                                     */
+    static void cylinder (Triangulation<3> &tria,
+                         const double      radius = 1.,
+                         const double      half_length = 1.);
+
+                                    /**
+                                     * Projection of the
+                                     * three-dimensional cylinder
+                                     * into the @p{xy}-plane.
+                                     * Therefore, this is simply a square.
+                                     *
+                                     * Coloring is like in 3D.
+                                     */
+    static void cylinder (Triangulation<2> &tria,
+                         const double      radius = 1.,
+                         const double      half_length = 1.);
+
+                                    /**
+                                     * Non-implemented dummy for compilation purposes.
+                                     */
+    static void cylinder (Triangulation<1> &tria,
+                         const double      radius,
+                         const double      half_length);
+
+
+
                                     /**
                                      * Initialize the given
                                      * triangulation with a hyper-L
index 68ea8c8100ea261d9482fab5999b2ec23444a0f7..11144ef60d0dd284f3fe48801cb0df60faa9e9e7 100644 (file)
 #include <grid/tria_boundary.h>
 
 
+/**
+ * Boundary object for the hull of a cylinder.  In three dimensions,
+ * points are projected on a circular tube along the @p{x}-axis. The
+ * radius of the tube can be set. Similar to @ref{HyperBallBoundary},
+ * new points are projected by dividing the straight line between the
+ * old two points and adjusting the radius in the @p{yz}-plane.
+ *
+ * This class was developed to be used in conjunction with the
+ * @p{cylinder} function of @ref{GridGenerator}. It should be used for
+ * the hull of the cylinder only (boundary indicator 0).
+ *
+ *   This class is derived from @ref{StraightBoundary} rather than from
+ *   @ref{Boundary}, which would seem natural, since this way we can use the
+ *   @ref{StraightBoundary}@p{<dim>::in_between(neighbors)} function.
+ *
+ *   @author Guido Kanschat, 2001
+ */
+template <int dim>
+class CylinderBoundary : public StraightBoundary<dim>
+{
+  public:
+                                    /**
+                                     * Constructor
+                                     */
+    CylinderBoundary (const double     radius = 1.0);
+
+                                    /**
+                                     * 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,
+                                    typename std::vector<Point<dim> > &points) const;
+
+                                    /**
+                                     * Refer to the general
+                                     * documentation of this class
+                                     * and the documentation of the
+                                     * base class.
+                                     *
+                                     * Only implemented for @p{dim=3}
+                                     * and for @p{points.size()==1}.
+                                     */
+    virtual void
+    get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterator &quad,
+                                    typename 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;
+
+                                    /**
+                                     * Return the radius of the cylinder.
+                                     */
+    double get_radius () const;
+
+                                    /**
+                                     * Exception. Thrown by the
+                                     * @p{get_radius} if the
+                                     * @p{compute_radius_automatically},
+                                     * see below, flag is set true.
+                                     */
+    DeclException0 (ExcRadiusNotSet);
+    
+    
+  protected:
+                                    /**
+                                     * Radius of the cylinder.
+                                     */
+    const double radius;
+
+  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,
+                                                typename std::vector<Point<dim> > &points) const;    
+};
+
+
+
 /**
  *   Specialisation of @ref{Boundary}<dim>, which places the new point on
  *   the boundary of a ball in arbitrary dimension. It works by projecting
index a0f3e382a08ca3037293cabd86cc40b796915a03..b9130cf1212a08dc9950922e54c27067ef868432 100644 (file)
@@ -133,6 +133,7 @@ GridGenerator::colorize_hyper_rectangle (Triangulation<dim> &tria)
            cell->face(3)->set_boundary_indicator (1);
            cell->face(4)->set_boundary_indicator (5);
            cell->face(5)->set_boundary_indicator (0);
+           break;
       default:
            Assert(false, ExcNotImplemented());
     };
@@ -186,6 +187,25 @@ void GridGenerator::hyper_ball (Triangulation<1> &,
 
 
 
+void
+GridGenerator::cylinder (Triangulation<2> &,
+                        const double,
+                        const double)
+{
+  Assert (false, ExcInternalError());  
+}
+
+
+
+void GridGenerator::cylinder (Triangulation<1> &,
+                             const double,
+                             const double)
+{
+  Assert (false, ExcInternalError());
+};
+
+
+
 void GridGenerator::hyper_shell (Triangulation<1> &,
                                 const Point<1> &,
                                 const double,
@@ -413,6 +433,38 @@ void GridGenerator::hyper_shell (Triangulation<2>   &tria,
 
 
 
+void
+GridGenerator::cylinder (Triangulation<2> &tria,
+                        const double radius,
+                        const double half_length)
+{
+  Point<2> p1 (-half_length, -radius);
+  Point<2> p2 (half_length, radius);
+
+  hyper_rectangle(tria, p1, p2, true);
+
+  Triangulation<2>::face_iterator f = tria.begin_face();
+  Triangulation<2>::face_iterator end = tria.end_face();
+  while (f != end)
+    {
+      switch (f->boundary_indicator())
+       {
+         case 0:
+           f->set_boundary_indicator(1);
+           break;
+         case 1:
+           f->set_boundary_indicator(2);
+           break;
+         default:
+           f->set_boundary_indicator(0);
+           break;          
+       }
+      ++f;
+    }
+}
+
+
+
 void
 GridGenerator::half_hyper_shell (Triangulation<2>   &tria,
                                 const Point<2>     &center,
@@ -697,6 +749,81 @@ GridGenerator::hyper_ball (Triangulation<3> &tria,
 
 
 
+void
+GridGenerator::cylinder (Triangulation<3> &tria,
+                        const double radius,
+                        const double half_length)
+{
+                                  // Copy the base from hyper_ball<2>
+                                  // and transform it to yz
+  const double d = radius/std::sqrt(2.0);
+  const double a = d/(1+std::sqrt(2.0));
+  const Point<3> vertices[24] = {
+    Point<3>(-half_length, -d,-d),
+      Point<3>(-half_length,  d,-d),
+      Point<3>(-half_length, -a,-a),
+      Point<3>(-half_length,  a,-a),
+      Point<3>(-half_length, -a, a),
+      Point<3>(-half_length,  a, a),
+      Point<3>(-half_length, -d, d),
+      Point<3>(-half_length,  d, d),
+    Point<3>(0, -d,-d),
+      Point<3>(0,  d,-d),
+      Point<3>(0, -a,-a),
+      Point<3>(0,  a,-a),
+      Point<3>(0, -a, a),
+      Point<3>(0,  a, a),
+      Point<3>(0, -d, d),
+      Point<3>(0,  d, d),
+    Point<3>(half_length, -d,-d),
+      Point<3>(half_length,  d,-d),
+      Point<3>(half_length, -a,-a),
+      Point<3>(half_length,  a,-a),
+      Point<3>(half_length, -a, a),
+      Point<3>(half_length,  a, a),
+      Point<3>(half_length, -d, d),
+      Point<3>(half_length,  d, d),
+      };
+  
+  int cell_vertices[10][8] = {
+       {0,1,3,2,8,9,11,10},
+       {0,2,4,6,8,10,12,14},
+       {2,3,5,4,10,11,13,12},
+       {1,7,5,3,9,15,13,11},
+       {6,4,5,7,14,12,13,15}
+  };
+  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<3> > cells (10, CellData<3>());
+  
+  for (unsigned int i=0; i<10; ++i) 
+    {
+      for (unsigned int j=0; j<8; ++j)
+       cells[i].vertices[j] = cell_vertices[i][j];
+      cells[i].material_id = 0;
+    };
+  
+  tria.create_triangulation (std::vector<Point<3> >(&vertices[0], &vertices[24]),
+                            cells,
+                            SubCellData());       // no boundary information
+
+  Triangulation<3>::cell_iterator cell = tria.begin();
+  Triangulation<3>::cell_iterator end = tria.end();
+  
+  while (cell != end)
+    {
+      if (cell->face(0)->boundary_indicator() != 255)
+       cell->face(0)->set_boundary_indicator(1);
+      if (cell->face(1)->boundary_indicator() != 255)
+       cell->face(1)->set_boundary_indicator(2);
+      ++cell;
+    }
+};
+
+
+
 void GridGenerator::hyper_shell (Triangulation<3>   &,
                                 const Point<3>     &,
                                 const double        ,
index 968a2ffdd46bd73dc08e45b0b28221c82d44c7e4..09706c59087990259ddc59f8214b3d164b619867 100644 (file)
@@ -31,85 +31,7 @@ template <int dim>
 void GridOut::write_dx (const Triangulation<dim> &tria,
                        std::ostream             &out) 
 {
-  AssertThrow (out, ExcIO());
-
-                                  // get the positions of the
-                                  // vertices and whether they are
-                                  // used.
-  const std::vector<Point<dim> > &vertices    = tria.get_vertices();
-  const std::vector<bool>        &vertex_used = tria.get_used_vertices();
-  
-  const unsigned int n_vertices = tria.n_used_vertices();
-
-  typename Triangulation<dim>::active_cell_iterator       cell=tria.begin_active();
-  const typename Triangulation<dim>::active_cell_iterator endc=tria.end();
-
-                                  // start with ucd data
-  out << n_vertices << ' '
-      << tria.n_active_cells() + (ucd_flags.write_faces ?
-                                 n_boundary_faces(tria) :
-                                 0)
-      << " 0 0 0"                  // no data
-      << std::endl;
-
-                                  // actually write the vertices.
-                                  // note that we shall number them
-                                  // with first index 1 instead of 0
-  for (unsigned int i=0; i<vertices.size(); ++i)
-    if (vertex_used[i])
-      {
-       out << i+1                 // vertex index
-           << "  "
-           << vertices[i];
-       for (unsigned int d=dim+1; d<=3; ++d)
-         out << " 0";             // fill with zeroes
-       out << std::endl;
-      };
-       
-                                  // write cells. Enumerate cells
-                                  // consecutively, starting with 1
-  unsigned int cell_index=1;
-  for (cell=tria.begin_active();
-       cell!=endc; ++cell, ++cell_index)
-    {
-      out << cell_index << ' '
-         << static_cast<unsigned int>(cell->material_id())
-         << " ";
-      switch (dim) 
-       {
-         case 1:  out << "line    "; break;
-         case 2:  out << "quad    "; break;
-         case 3:  out << "hex     "; break;
-         default:
-               Assert (false, ExcNotImplemented());
-       };
-
-                                      // it follows a list of the
-                                      // vertices of each cell. in 1d
-                                      // this is simply a list of the
-                                      // two vertices, in 2d its counter
-                                      // clockwise, as usual in this
-                                      // library. in 3d, the same applies
-                                      // (special thanks to AVS for
-                                      // numbering their vertices in a
-                                      // way compatible to deal.II!)
-                                      //
-                                      // technical reference:
-                                      // AVS Developer's Guide, Release 4,
-                                      // May, 1992, p. E6
-                                      //
-                                      // note: vertex numbers are 1-base
-      for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell;
-          ++vertex)
-       out << cell->vertex_index(vertex)+1 << ' ';
-      out << std::endl;
-    };
-
-                                  // write faces with non-zero boundary
-                                  // indicator
-  if (ucd_flags.write_faces)
-    write_ucd_faces (tria, cell_index, out);
-    
+  Assert(false, ExcNotImplemented());
   AssertThrow (out, ExcIO());
 };
 
index 27b39a1f43fe2193e7df0ddf5c6f7a6f748ae1d6..678821e3f48dc4bca8859d6eab11d3c9221f5cbc 100644 (file)
 #include <cmath>
 
 
+template <int dim>
+CylinderBoundary<dim>::CylinderBoundary (const double radius) :
+               radius(radius)
+{};
+
+
+
+template <int dim>
+Point<dim>
+CylinderBoundary<dim>::get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const
+{
+  Point<dim> middle = StraightBoundary<dim>::get_new_point_on_line (line);
+    
+                                  // project to boundary
+  if (dim>=3)
+    middle *= radius / std::sqrt(middle.square()-middle(0)*middle(0));
+  return middle;
+};
+
+
+template <int dim>
+Point<dim>
+CylinderBoundary<dim>::
+get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const
+{
+  Point<dim> middle = StraightBoundary<dim>::get_new_point_on_quad (quad);
+  
+                                  // project to boundary
+  if (dim>=3)
+    middle *= radius / std::sqrt(middle.square()-middle(0)*middle(0));
+  return middle;
+};
+
+
+template <int dim>
+void
+CylinderBoundary<dim>::get_intermediate_points_on_line (
+  const typename Triangulation<dim>::line_iterator &line,
+  typename 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);
+}
+
+  
+template <int dim>
+void
+CylinderBoundary<dim>::get_intermediate_points_between_points (
+  const Point<dim> &v0, const Point<dim> &v1,
+  typename std::vector<Point<dim> > &points) const
+{
+  const unsigned int n=points.size();
+  Assert(n>0, ExcInternalError());
+                                  // Do a simple linear interpolation
+                                  // followed by projection. If this
+                                  // is not sufficient, try to
+                                  // understand the sophisticated
+                                  // code in HyperBall later.
+  Point<dim> ds = v1-v0;
+  ds /= n+1;
+  
+  for (unsigned int i=0; i<n; ++i)
+    {
+      if (i==0)
+       points[i] = v0+ds;
+      else
+       points[i] = points[i-1]+ds;
+      
+      points[i] *= radius / std::sqrt(points[i].square()-points[i](0)*points[i](0));
+    }
+}
+
+
+#if deal_II_dimension == 3
+
+template <>
+void
+CylinderBoundary<3>::get_intermediate_points_on_quad (
+  const Triangulation<3>::quad_iterator &quad,
+  std::vector<Point<3> > &points) const
+{
+  unsigned int m=static_cast<unsigned int> (std::sqrt(static_cast<double>(points.size())));
+  Assert(points.size()==m*m, ExcInternalError());
+
+  std::vector<Point<3> > lp3(m);
+  std::vector<Point<3> > lp1(m);
+  
+  get_intermediate_points_on_line(quad->line(3), lp3);
+  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(lp3[i], lp1[i], lps);
+      
+      for (unsigned int j=0; j<m; ++j)
+       points[i*m+j]=lps[j];
+    }
+}
+
+#endif
+
+
+template <int dim>
+void
+CylinderBoundary<dim>::get_intermediate_points_on_quad (
+  const typename Triangulation<dim>::quad_iterator &,
+  typename std::vector<Point<dim> > &) const
+{
+  Assert(false, Boundary<dim>::ExcFunctionNotUseful(dim));
+}
+
+
+
+#if deal_II_dimension == 1
+
+template <>
+void
+CylinderBoundary<1>::
+get_normals_at_vertices (const Triangulation<1>::face_iterator &,
+                        Boundary<1>::FaceVertexNormals &) const
+{
+  Assert (false, Boundary<1>::ExcFunctionNotUseful(1));
+};
+
+#endif
+
+
+template <int dim>
+void
+CylinderBoundary<dim>::
+get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &face,
+                        typename Boundary<dim>::FaceVertexNormals &face_vertex_normals) const
+{
+  for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_face; ++vertex)
+    {
+      face_vertex_normals[vertex] = face->vertex(vertex);
+      face_vertex_normals[vertex][0] = 0.;
+    }
+};
+
+
+
+template <int dim>
+double
+CylinderBoundary<dim>::get_radius () const 
+{
+  return radius;
+};
+
+
+//======================================================================//
+
 template <int dim>
 HyperBallBoundary<dim>::HyperBallBoundary (const Point<dim> p,
                                           const double     radius) :
@@ -597,6 +752,7 @@ 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 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.