From af2477fff3b6a99491a82ec7820ab0aada3c76be Mon Sep 17 00:00:00 2001
From: wolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Mon, 19 Apr 2004 21:09:40 +0000
Subject: [PATCH] Implement a cylinder shell in 3d.

git-svn-id: https://svn.dealii.org/trunk@9066 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/deal.II/include/grid/grid_generator.h |  22 +++
 deal.II/deal.II/source/grid/grid_generator.cc |  80 ++++++++
 tests/bits/Makefile                           |   3 +-
 tests/bits/cylinder_shell_01.cc               |  44 +++++
 .../bits/cylinder_shell_01.output             | 175 ++++++++++++++++++
 5 files changed, 323 insertions(+), 1 deletion(-)
 create mode 100644 tests/bits/cylinder_shell_01.cc
 create mode 100644 tests/results/i686-pc-linux-gnu+gcc3.2/bits/cylinder_shell_01.output

diff --git a/deal.II/deal.II/include/grid/grid_generator.h b/deal.II/deal.II/include/grid/grid_generator.h
index 206d93254a..bd07bcf03e 100644
--- a/deal.II/deal.II/include/grid/grid_generator.h
+++ b/deal.II/deal.II/include/grid/grid_generator.h
@@ -536,6 +536,28 @@ class GridGenerator
 			     const double        outer_radius,
 			     const unsigned int  n_cells = 0);
 
+				     /**
+				      * Produce a domain that is the space
+				      * between two cylinders in 3d, with
+				      * given length, inner and outer radius
+				      * and a given number of elements for
+				      * this initial triangulation. If @p
+				      * n_radial_cells is zero (as is the
+				      * default), then it is computed
+				      * adaptively such that the resulting
+				      * elements have the least aspect
+				      * ratio. The same holds for @p
+				      * n_axial_cells.
+				      *
+				      * The triangulation needs to be void
+				      * upon calling this function.
+				      */
+    static void cylinder_shell (Triangulation<3>   &tria,
+                                const double        length,
+                                const double        inner_radius,
+                                const double        outer_radius,
+                                const unsigned int  n_radial_cells = 0,
+                                const unsigned int  n_axial_cells = 0);
 
 				     /**
 				      * This class produces a half hyper-ball,
diff --git a/deal.II/deal.II/source/grid/grid_generator.cc b/deal.II/deal.II/source/grid/grid_generator.cc
index 55e392be21..8b1b7a1006 100644
--- a/deal.II/deal.II/source/grid/grid_generator.cc
+++ b/deal.II/deal.II/source/grid/grid_generator.cc
@@ -1172,6 +1172,86 @@ void GridGenerator::hyper_shell (Triangulation<3>   &,
 }
 
 
+
+void GridGenerator::cylinder_shell (Triangulation<3>   &tria,
+                                    const double        length,
+                                    const double        inner_radius,
+                                    const double        outer_radius,
+                                    const unsigned int  n_radial_cells,
+                                    const unsigned int  n_axial_cells)
+{
+  Assert ((inner_radius > 0) && (inner_radius < outer_radius),
+	  ExcInvalidRadii ());
+  
+  const double pi = deal_II_numbers::PI;
+  
+				   // 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_r = (n_radial_cells == 0 ?
+                            static_cast<unsigned int>
+                            (std::ceil((2*pi* (outer_radius + inner_radius)/2) /
+                                       (outer_radius - inner_radius))) :
+                            n_radial_cells);
+  const unsigned int N_z = (n_axial_cells == 0 ?
+                            static_cast<unsigned int>
+                            (std::ceil (length /
+                                        (2*pi*(outer_radius + inner_radius)/2/N_r))) :
+                            n_axial_cells);
+
+				   // set up N vertices on the
+				   // outer and N vertices on
+				   // the inner circle. the
+				   // first N ones are on the
+				   // outer one, and all are
+				   // numbered counter-clockwise
+  std::vector<Point<2> > vertices_2d(2*N_r);
+  for (unsigned int i=0; i<N_r; ++i)
+    {
+      vertices_2d[i] = Point<2>( std::cos(2*pi * i/N_r),
+                                 std::sin(2*pi * i/N_r)) * outer_radius;
+      vertices_2d[i+N_r] = vertices_2d[i] * (inner_radius/outer_radius);
+    };
+
+  std::vector<Point<3> > vertices_3d;
+  vertices_3d.reserve (2*N_r*(N_z+1));
+  for (unsigned int j=0; j<=N_z; ++j)
+    for (unsigned int i=0; i<2*N_r; ++i)
+      {
+        const Point<3> v (vertices_2d[i][0],
+                          vertices_2d[i][1],
+                          j*length/N_z);
+        vertices_3d.push_back (v);
+      }
+                            
+  std::vector<CellData<3> > cells (N_r*N_z, CellData<3>());
+  
+  for (unsigned int j=0; j<N_z; ++j)
+    for (unsigned int i=0; i<N_r; ++i) 
+      {
+        cells[i+j*N_r].vertices[0] = i + j*2*N_r;
+        cells[i+j*N_r].vertices[1] = (i+1)%N_r + j*2*N_r;
+        cells[i+j*N_r].vertices[2] = N_r+((i+1)%N_r) + j*2*N_r;
+        cells[i+j*N_r].vertices[3] = N_r+i + j*2*N_r;
+        
+        cells[i+j*N_r].vertices[4] = i + (j+1)*2*N_r;
+        cells[i+j*N_r].vertices[5] = (i+1)%N_r + (j+1)*2*N_r;
+        cells[i+j*N_r].vertices[6] = N_r+((i+1)%N_r) + (j+1)*2*N_r;
+        cells[i+j*N_r].vertices[7] = N_r+i + (j+1)*2*N_r;
+        
+        cells[i+j*N_r].material_id = 0;
+      }
+  
+  tria.create_triangulation (vertices_3d, cells, SubCellData());
+}
+
+
 #endif
 
 
diff --git a/tests/bits/Makefile b/tests/bits/Makefile
index 37e954b547..1726bb96c3 100644
--- a/tests/bits/Makefile
+++ b/tests/bits/Makefile
@@ -57,7 +57,8 @@ tests_x = anna_? \
 	  subdomain_ids_* \
 	  dof_constraints_* \
 	  apply_boundary_values_* \
-	  error_estimator_*
+	  error_estimator_* \
+	  cylinder_shell_*
 
 ifeq ($(USE_CONTRIB_PETSC),yes)
   tests_x += petsc_*
diff --git a/tests/bits/cylinder_shell_01.cc b/tests/bits/cylinder_shell_01.cc
new file mode 100644
index 0000000000..9517325cb1
--- /dev/null
+++ b/tests/bits/cylinder_shell_01.cc
@@ -0,0 +1,44 @@
+//----------------------------  cylinder_shell_01.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 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
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  cylinder_shell_01.cc  ---------------------------
+
+// test that the grid generated by GridGenerator::cylinder_shell<3> works as
+// expected
+
+
+#include "../tests.h"
+#include <base/logstream.h>
+#include <grid/tria.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_accessor.h>
+#include <grid/grid_generator.h>
+#include <fstream>
+
+    
+
+int main () 
+{
+  std::ofstream logfile("cylinder_shell_01.output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  logfile.precision (2);
+
+                                   // generate a hyperball in 3d
+  Triangulation<3> tria;
+  GridGenerator::cylinder_shell (tria, 1, .8, 1);
+
+                                   // make sure that all cells have positive
+                                   // volume
+  for (Triangulation<3>::active_cell_iterator cell=tria.begin_active();
+       cell!=tria.end(); ++cell)
+    deallog << cell << ' ' << cell->measure () << std::endl;
+}
diff --git a/tests/results/i686-pc-linux-gnu+gcc3.2/bits/cylinder_shell_01.output b/tests/results/i686-pc-linux-gnu+gcc3.2/bits/cylinder_shell_01.output
new file mode 100644
index 0000000000..2bb791c4b2
--- /dev/null
+++ b/tests/results/i686-pc-linux-gnu+gcc3.2/bits/cylinder_shell_01.output
@@ -0,0 +1,175 @@
+
+DEAL::0.0 0.0064
+DEAL::0.1 0.0064
+DEAL::0.2 0.0064
+DEAL::0.3 0.0064
+DEAL::0.4 0.0064
+DEAL::0.5 0.0064
+DEAL::0.6 0.0064
+DEAL::0.7 0.0064
+DEAL::0.8 0.0064
+DEAL::0.9 0.0064
+DEAL::0.10 0.0064
+DEAL::0.11 0.0064
+DEAL::0.12 0.0064
+DEAL::0.13 0.0064
+DEAL::0.14 0.0064
+DEAL::0.15 0.0064
+DEAL::0.16 0.0064
+DEAL::0.17 0.0064
+DEAL::0.18 0.0064
+DEAL::0.19 0.0064
+DEAL::0.20 0.0064
+DEAL::0.21 0.0064
+DEAL::0.22 0.0064
+DEAL::0.23 0.0064
+DEAL::0.24 0.0064
+DEAL::0.25 0.0064
+DEAL::0.26 0.0064
+DEAL::0.27 0.0064
+DEAL::0.28 0.0064
+DEAL::0.29 0.0064
+DEAL::0.30 0.0064
+DEAL::0.31 0.0064
+DEAL::0.32 0.0064
+DEAL::0.33 0.0064
+DEAL::0.34 0.0064
+DEAL::0.35 0.0064
+DEAL::0.36 0.0064
+DEAL::0.37 0.0064
+DEAL::0.38 0.0064
+DEAL::0.39 0.0064
+DEAL::0.40 0.0064
+DEAL::0.41 0.0064
+DEAL::0.42 0.0064
+DEAL::0.43 0.0064
+DEAL::0.44 0.0064
+DEAL::0.45 0.0064
+DEAL::0.46 0.0064
+DEAL::0.47 0.0064
+DEAL::0.48 0.0064
+DEAL::0.49 0.0064
+DEAL::0.50 0.0064
+DEAL::0.51 0.0064
+DEAL::0.52 0.0064
+DEAL::0.53 0.0064
+DEAL::0.54 0.0064
+DEAL::0.55 0.0064
+DEAL::0.56 0.0064
+DEAL::0.57 0.0064
+DEAL::0.58 0.0064
+DEAL::0.59 0.0064
+DEAL::0.60 0.0064
+DEAL::0.61 0.0064
+DEAL::0.62 0.0064
+DEAL::0.63 0.0064
+DEAL::0.64 0.0064
+DEAL::0.65 0.0064
+DEAL::0.66 0.0064
+DEAL::0.67 0.0064
+DEAL::0.68 0.0064
+DEAL::0.69 0.0064
+DEAL::0.70 0.0064
+DEAL::0.71 0.0064
+DEAL::0.72 0.0064
+DEAL::0.73 0.0064
+DEAL::0.74 0.0064
+DEAL::0.75 0.0064
+DEAL::0.76 0.0064
+DEAL::0.77 0.0064
+DEAL::0.78 0.0064
+DEAL::0.79 0.0064
+DEAL::0.80 0.0064
+DEAL::0.81 0.0064
+DEAL::0.82 0.0064
+DEAL::0.83 0.0064
+DEAL::0.84 0.0064
+DEAL::0.85 0.0064
+DEAL::0.86 0.0064
+DEAL::0.87 0.0064
+DEAL::0.88 0.0064
+DEAL::0.89 0.0064
+DEAL::0.90 0.0064
+DEAL::0.91 0.0064
+DEAL::0.92 0.0064
+DEAL::0.93 0.0064
+DEAL::0.94 0.0064
+DEAL::0.95 0.0064
+DEAL::0.96 0.0064
+DEAL::0.97 0.0064
+DEAL::0.98 0.0064
+DEAL::0.99 0.0064
+DEAL::0.100 0.0064
+DEAL::0.101 0.0064
+DEAL::0.102 0.0064
+DEAL::0.103 0.0064
+DEAL::0.104 0.0064
+DEAL::0.105 0.0064
+DEAL::0.106 0.0064
+DEAL::0.107 0.0064
+DEAL::0.108 0.0064
+DEAL::0.109 0.0064
+DEAL::0.110 0.0064
+DEAL::0.111 0.0064
+DEAL::0.112 0.0064
+DEAL::0.113 0.0064
+DEAL::0.114 0.0064
+DEAL::0.115 0.0064
+DEAL::0.116 0.0064
+DEAL::0.117 0.0064
+DEAL::0.118 0.0064
+DEAL::0.119 0.0064
+DEAL::0.120 0.0064
+DEAL::0.121 0.0064
+DEAL::0.122 0.0064
+DEAL::0.123 0.0064
+DEAL::0.124 0.0064
+DEAL::0.125 0.0064
+DEAL::0.126 0.0064
+DEAL::0.127 0.0064
+DEAL::0.128 0.0064
+DEAL::0.129 0.0064
+DEAL::0.130 0.0064
+DEAL::0.131 0.0064
+DEAL::0.132 0.0064
+DEAL::0.133 0.0064
+DEAL::0.134 0.0064
+DEAL::0.135 0.0064
+DEAL::0.136 0.0064
+DEAL::0.137 0.0064
+DEAL::0.138 0.0064
+DEAL::0.139 0.0064
+DEAL::0.140 0.0064
+DEAL::0.141 0.0064
+DEAL::0.142 0.0064
+DEAL::0.143 0.0064
+DEAL::0.144 0.0064
+DEAL::0.145 0.0064
+DEAL::0.146 0.0064
+DEAL::0.147 0.0064
+DEAL::0.148 0.0064
+DEAL::0.149 0.0064
+DEAL::0.150 0.0064
+DEAL::0.151 0.0064
+DEAL::0.152 0.0064
+DEAL::0.153 0.0064
+DEAL::0.154 0.0064
+DEAL::0.155 0.0064
+DEAL::0.156 0.0064
+DEAL::0.157 0.0064
+DEAL::0.158 0.0064
+DEAL::0.159 0.0064
+DEAL::0.160 0.0064
+DEAL::0.161 0.0064
+DEAL::0.162 0.0064
+DEAL::0.163 0.0064
+DEAL::0.164 0.0064
+DEAL::0.165 0.0064
+DEAL::0.166 0.0064
+DEAL::0.167 0.0064
+DEAL::0.168 0.0064
+DEAL::0.169 0.0064
+DEAL::0.170 0.0064
+DEAL::0.171 0.0064
+DEAL::0.172 0.0064
+DEAL::0.173 0.0064
-- 
2.39.5