From a5f3b79f089478d25aaa1b1bd2150028fe467e6a Mon Sep 17 00:00:00 2001
From: hartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Fri, 15 Oct 2004 07:00:24 +0000
Subject: [PATCH] Extension of CylinderBoundary class to represent circular
 tubes also along the y- and z-axis.

git-svn-id: https://svn.dealii.org/trunk@9710 0785d39b-7218-0410-832d-ea1e28bc413d
---
 .../deal.II/include/grid/tria_boundary_lib.h  | 32 ++++++++++---
 .../deal.II/source/grid/tria_boundary_lib.cc  | 47 ++++++++++---------
 deal.II/doc/news/c-5.0.html                   |  9 ++++
 3 files changed, 61 insertions(+), 27 deletions(-)

diff --git a/deal.II/deal.II/include/grid/tria_boundary_lib.h b/deal.II/deal.II/include/grid/tria_boundary_lib.h
index 8b2c7c9216..4d0676c3ad 100644
--- a/deal.II/deal.II/include/grid/tria_boundary_lib.h
+++ b/deal.II/deal.II/include/grid/tria_boundary_lib.h
@@ -20,10 +20,11 @@
 
 /**
  * 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 HyperBallBoundary,
- * new points are projected by dividing the straight line between the
- * old two points and adjusting the radius in the @p yz-plane.
+ * points are projected on a circular tube along the <tt>x-</tt>,
+ * <tt>y-</tt> or <tt>z</tt>-axis. The radius of the tube can be
+ * set. Similar to HyperBallBoundary, new points are projected by
+ * dividing the straight line between the old two points and adjusting
+ * the radius in the @p yz-plane (xz-plane or xy-plane, respectively).
  *
  * This class was developed to be used in conjunction with the
  * @p cylinder function of GridGenerator. It should be used for
@@ -40,9 +41,16 @@ class CylinderBoundary : public StraightBoundary<dim>
 {
   public:
 				     /**
-				      * Constructor
+				      * Constructor. Per default
+				      * circular tube along the x-axis
+				      * (<tt>axis=0<\tt>). Choose
+				      * <tt>axis=1<\tt> or
+				      * <tt>axis=2<\tt> for a tube
+				      * along the y- or z-axis,
+				      * respectively.
 				      */
-    CylinderBoundary (const double     radius = 1.0);
+    CylinderBoundary (const double       radius = 1.0,
+		      const unsigned int axis   = 0);
 
 				     /**
 				      * Refer to the general documentation of
@@ -120,6 +128,18 @@ class CylinderBoundary : public StraightBoundary<dim>
 				      */
     const double radius;
 
+				     /**
+				      * Denotes the axis of the
+				      * circular
+				      * tube. <tt>axis=0<\tt>,
+				      * <tt>axis=1<\tt> and
+				      * <tt>axis=2<\tt> denote the
+				      * <tt>x-axis</tt>,
+				      * <tt>y-axis</tt> and
+				      * <tt>z-axis</tt>, respectively.
+				      */
+    const unsigned int axis;
+
   private:
 
 				     /**
diff --git a/deal.II/deal.II/source/grid/tria_boundary_lib.cc b/deal.II/deal.II/source/grid/tria_boundary_lib.cc
index b61bb505d6..97b6a9bdf2 100644
--- a/deal.II/deal.II/source/grid/tria_boundary_lib.cc
+++ b/deal.II/deal.II/source/grid/tria_boundary_lib.cc
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003 by the deal.II authors
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -21,8 +21,10 @@
 
 
 template <int dim>
-CylinderBoundary<dim>::CylinderBoundary (const double radius) :
-		radius(radius)
+CylinderBoundary<dim>::CylinderBoundary (const double radius,
+					 const unsigned int axis) :
+		radius(radius),
+		axis(axis)
 {}
 
 
@@ -34,12 +36,13 @@ CylinderBoundary<dim>::get_new_point_on_line (const typename Triangulation<dim>:
   Point<dim> middle = StraightBoundary<dim>::get_new_point_on_line (line);
 				   // project to boundary
   if (dim>=3
-      && line->vertex(0).square()-line->vertex(0)(0)*line->vertex(0)(0) >= radius*radius-1.e-12
-      && line->vertex(1).square()-line->vertex(1)(0)*line->vertex(1)(0) >= radius*radius-1.e-12)
+      && line->vertex(0).square()-line->vertex(0)(axis)*line->vertex(0)(axis) >= radius*radius-1.e-10
+      && line->vertex(1).square()-line->vertex(1)(axis)*line->vertex(1)(axis) >= radius*radius-1.e-10)
     {
-      const double f = radius / std::sqrt(middle.square()-middle(0)*middle(0));
-      for (unsigned int i=1;i<dim;++i)
-	middle(i) *= f;
+      const double f = radius / std::sqrt(middle.square()-middle(axis)*middle(axis));
+      for (unsigned int i=0; i<dim; ++i)
+	if (i!=axis)
+	  middle(i) *= f;
     }
   return middle;
 }
@@ -55,15 +58,16 @@ get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const
   Point<3> middle = StraightBoundary<3>::get_new_point_on_quad (quad);
   
 				   // project to boundary
-  if (quad->vertex(0).square()-quad->vertex(0)(0)*quad->vertex(0)(0) >= radius*radius-1.e-12
-      && quad->vertex(1).square()-quad->vertex(1)(0)*quad->vertex(1)(0) >= radius*radius-1.e-12
-      && quad->vertex(2).square()-quad->vertex(2)(0)*quad->vertex(2)(0) >= radius*radius-1.e-12
-      && quad->vertex(3).square()-quad->vertex(3)(0)*quad->vertex(3)(0) >= radius*radius-1.e-12)
+  if (quad->vertex(0).square()-quad->vertex(0)(axis)*quad->vertex(0)(axis) >= radius*radius-1.e-10
+      && quad->vertex(1).square()-quad->vertex(1)(axis)*quad->vertex(1)(axis) >= radius*radius-1.e-10
+      && quad->vertex(2).square()-quad->vertex(2)(axis)*quad->vertex(2)(axis) >= radius*radius-1.e-10
+      && quad->vertex(3).square()-quad->vertex(3)(axis)*quad->vertex(3)(axis) >= radius*radius-1.e-10)
     
     {
-      const double f = radius / std::sqrt(middle.square()-middle(0)*middle(0));
-      for (unsigned int i=1;i<3;++i)
-	middle(i) *= f;
+      const double f = radius / std::sqrt(middle.square()-middle(axis)*middle(axis));
+      for (unsigned int i=0; i<3; ++i)
+	if (i!=axis)
+	  middle(i) *= f;
     }
   return middle;
 }
@@ -112,8 +116,8 @@ CylinderBoundary<dim>::get_intermediate_points_between_points (
   ds /= n+1;
 
   bool scale = (dim>=3
-		&& v0.square()-v0(0)*v0(0) >= radius*radius-1.e-12
-		&& v1.square()-v1(0)*v1(0) >= radius*radius-1.e-12);
+		&& v0.square()-v0(axis)*v0(axis) >= radius*radius-1.e-10
+		&& v1.square()-v1(axis)*v1(axis) >= radius*radius-1.e-10);
   
   for (unsigned int i=0; i<n; ++i)
     {
@@ -124,9 +128,10 @@ CylinderBoundary<dim>::get_intermediate_points_between_points (
 
       if (scale)
 	{
-	  const double f = radius / std::sqrt(points[i].square()-points[i](0)*points[i](0));
-	  for (unsigned int d=1;d<dim;++d)
-	    points[i](d) *= f;
+	  const double f = radius / std::sqrt(points[i].square()-points[i](axis)*points[i](axis));
+	  for (unsigned int d=0; d<dim; ++d)
+	    if (i!=axis)
+	      points[i](d) *= f;
 	}
     }
 }
@@ -196,7 +201,7 @@ get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &face,
   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.;
+      face_vertex_normals[vertex][axis] = 0.;
     }
 }
 
diff --git a/deal.II/doc/news/c-5.0.html b/deal.II/doc/news/c-5.0.html
index f78d446a8c..38f023982f 100644
--- a/deal.II/doc/news/c-5.0.html
+++ b/deal.II/doc/news/c-5.0.html
@@ -359,6 +359,15 @@ inconvenience this causes.
 <h3>deal.II</h3>
 
 <ol>
+  <li> <p>
+       Improved: The <code class="class">CylinderBoundary</code>
+       represented the hull of a circular tubes along the x-axis. It
+       is now extended to allow for circular tubes also along the y-
+       and z-axis.
+       <br>
+       (RH 2004/10/15)
+       </p>
+
   <li> <p>
        Fixed: The <code class="class">ConstraintMatrix</code> class had some
        algorithms that were linear in the number of constraints. Since these
-- 
2.39.5