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 b9130cf1212a08dc9950922e54c27067ef868432..79d6c6f64e6d48cc1aae2b76da45934d2afe8e96 100644 (file)
--- 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 678821e3f48dc4bca8859d6eab11d3c9221f5cbc..fa443d675b239abc2654233c7ead7bfef2eb4863 100644 (file)
--- 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;
+       }
     }
 }