From: guido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Wed, 19 Dec 2001 17:29:01 +0000 (+0000)
Subject: cylinder fixed
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=78b6a0062881fe2c9b8c4e378c782373a7625024;p=dealii-svn.git

cylinder fixed


git-svn-id: https://svn.dealii.org/trunk@5341 0785d39b-7218-0410-832d-ea1e28bc413d
---

diff --git a/deal.II/deal.II/source/grid/grid_generator.cc b/deal.II/deal.II/source/grid/grid_generator.cc
index b9130cf121..79d6c6f64e 100644
--- a/deal.II/deal.II/source/grid/grid_generator.cc
+++ b/deal.II/deal.II/source/grid/grid_generator.cc
@@ -758,32 +758,39 @@ GridGenerator::cylinder (Triangulation<3> &tria,
 				   // 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),
+  Point<3> vertices[24] = {
+    Point<3>(-d, -half_length,-d),
+      Point<3>( d, -half_length,-d),
+      Point<3>(-a, -half_length,-a),
+      Point<3>( a, -half_length,-a),
+      Point<3>(-a, -half_length, a),
+      Point<3>( a, -half_length, a),
+      Point<3>(-d, -half_length, d),
+      Point<3>( d, -half_length, d),
+      Point<3>(-d, 0,-d),
+      Point<3>( d, 0,-d),
+      Point<3>(-a, 0,-a),
+      Point<3>( a, 0,-a),
+      Point<3>(-a, 0, a),
+      Point<3>( a, 0, a),
+      Point<3>(-d, 0, d),
+      Point<3>( d, 0, d),
+      Point<3>(-d, half_length,-d),
+      Point<3>( d, half_length,-d),
+      Point<3>(-a, half_length,-a),
+      Point<3>( a, half_length,-a),
+      Point<3>(-a, half_length, a),
+      Point<3>( a, half_length, a),
+      Point<3>(-d, half_length, d),
+      Point<3>( d, half_length, d),
       };
+				   // Turn cylinder such that y->x
+  for (unsigned int i=0;i<24;++i)
+    {
+      const double h = vertices[i](1);
+      vertices[i](1) = -vertices[i](0);
+      vertices[i](0) = h;
+    }
   
   int cell_vertices[10][8] = {
 	{0,1,3,2,8,9,11,10},
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 678821e3f4..fa443d675b 100644
--- a/deal.II/deal.II/source/grid/tria_boundary_lib.cc
+++ b/deal.II/deal.II/source/grid/tria_boundary_lib.cc
@@ -32,10 +32,15 @@ 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));
+  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)
+    {
+      const double f = radius / std::sqrt(middle.square()-middle(0)*middle(0));
+      for (unsigned int i=1;i<dim;++i)
+	middle(i) *= f;
+    }
   return middle;
 };
 
@@ -48,8 +53,17 @@ get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) c
   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));
+  if (dim>=3
+      && 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)
+      
+    {
+      const double f = radius / std::sqrt(middle.square()-middle(0)*middle(0));
+      for (unsigned int i=1;i<dim;++i)
+	middle(i) *= f;
+    }
   return middle;
 };
 
@@ -82,6 +96,10 @@ CylinderBoundary<dim>::get_intermediate_points_between_points (
 				   // code in HyperBall later.
   Point<dim> ds = v1-v0;
   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);
   
   for (unsigned int i=0; i<n; ++i)
     {
@@ -89,8 +107,13 @@ CylinderBoundary<dim>::get_intermediate_points_between_points (
 	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 (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;
+	}
     }
 }