]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Check in a more interesting 3d geometry
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 25 May 2007 20:03:21 +0000 (20:03 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 25 May 2007 20:03:21 +0000 (20:03 +0000)
git-svn-id: https://svn.dealii.org/trunk@14708 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-27/step-27.cc

index 05db6bac683e1c73874a8faefd385391c9e01e28..7e28967828fef330e218c14ecad5a1a3592d52b4 100644 (file)
@@ -34,6 +34,7 @@
 #include <grid/tria_accessor.h>
 #include <grid/tria_iterator.h>
 #include <grid/tria_boundary_lib.h>
+#include <grid/grid_reordering.h>
 #include <dofs/dof_accessor.h>
 #include <dofs/dof_tools.h>
 #include <hp/fe_values.h>
@@ -112,10 +113,23 @@ double
 RightHandSide<dim>::value (const Point<dim>   &p,
                           const unsigned int  /*component*/) const
 {
-  double product = 1;
-  for (unsigned int d=0; d<dim; ++d)
-    product *= (p[d]+1);
-  return product;
+  switch (dim)
+    {
+      case 2:
+      {
+       double product = 1;
+       for (unsigned int d=0; d<dim; ++d)
+         product *= (p[d]+1);
+       return product;
+      }
+      
+      case 3:
+           return (p[0]>std::fabs(p[1]) ? 1 : 0);
+           
+      default:
+           Assert (false, ExcNotImplemented());
+    }
+  return 0.;
 }
 
 
@@ -124,7 +138,7 @@ RightHandSide<dim>::value (const Point<dim>   &p,
 template <int dim>
 LaplaceProblem<dim>::LaplaceProblem () :
   dof_handler (triangulation),
-  max_degree (dim == 2 ? 7 : 4)
+  max_degree (dim == 2 ? 7 : 5)
 {
   for (unsigned int degree=2; degree<=max_degree; ++degree)
     {
@@ -555,7 +569,7 @@ void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
     
     const std::string filename = "solution-" +
                                 Utilities::int_to_string (cycle, 2) +
-                                ".vtk";
+                                ".gmv";
     DataOut<dim,hp::DoFHandler<dim> > data_out;
 
     data_out.attach_dof_handler (dof_handler);
@@ -566,7 +580,7 @@ void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
     data_out.build_patches ();
   
     std::ofstream output (filename.c_str());
-    data_out.write_vtk (output);
+    data_out.write_gmv (output);
   }
 }
 
@@ -643,99 +657,584 @@ void LaplaceProblem<2>::create_coarse_grid ()
 
 
 
-template <>
-void LaplaceProblem<3>::create_coarse_grid ()
+namespace BreastPhantom
 {
-  const unsigned int dim = 3;
-
-  //  GridGenerator::hyper_cube (triangulation);
-  //  triangulation.refine_global (1);
-
-  // Create a hollow cube, in analogy to the 2D example.
-  // The grid generation is done in two steps. First the
-  // cell data is created on a uniform grid. In the
-  // second step, the unused vertices are removed.
-  const unsigned char hollow [4][4] = {{1,1,1,1},
-                                      {1,0,0,1},
-                                      {1,0,0,1},
-                                      {1,1,1,1}};
-  const unsigned char solid [4][4] = {{1,1,1,1},
-                                     {1,1,1,1},
-                                     {1,1,1,1},
-                                     {1,1,1,1}};
-  const unsigned char (*layers[4])[4][4] = {&solid, &hollow, &hollow, &solid};
-
-  std::vector<CellData<dim> > cells;
-  std::vector<bool> vertex_used (5*5*5, false);
-
-  for (unsigned int zc = 0; zc < 4; ++zc)
-    for (unsigned int yc = 0; yc < 4; ++yc)
-      for (unsigned int xc = 0; xc < 4; ++xc)
+  
+
+                                  // Radius of the sphere of the breast
+                                  // phantom geometry
+  static const double hemisphere_radius = 5;
+
+                                  // Radius of the disk underneath
+  static const double bottom_disk_radius = 10;
+
+                                  // Bottom z-coordinate of the disk
+                                  // underneath
+  const double bottom_disk_floor = -3;
+                                  // Top z-coordinate of the disk
+                                  // underneath
+  const double bottom_disk_ceil = -.5;
+
+                                  // radius of the inner set of cells of
+                                  // the sphere geometry
+  const double interior_hemisphere_radius
+  = hemisphere_radius/(1.+std::sqrt(2.0));
+      
+  template <int dim>
+  class SphereBoundary : public HyperBallBoundary<dim> 
+  {
+    public:
+      SphereBoundary () 
+                     :
+                     HyperBallBoundary<dim> (Point<dim>(), hemisphere_radius) 
+       {}
+  };
+  
+  
+  template <int dim>
+  class CylinderBoundary : public StraightBoundary<dim>
+  {
+    public:
+      typedef
+      typename Triangulation<dim>::line_iterator
+      line_iterator;
+      
+      typedef
+      typename Triangulation<dim>::quad_iterator
+      quad_iterator;
+      
+      typedef
+      typename Triangulation<dim>::face_iterator
+      face_iterator;
+
+                                      /**
+                                       * Constructor.
+                                       */
+      CylinderBoundary (const double radius)
+                     :
+                     radius (radius)
+       {}
+          
+          
+      virtual Point<dim>
+      get_new_point_on_line (const line_iterator &line) const;
+
+      virtual Point<dim>
+      get_new_point_on_quad (const quad_iterator &quad) const;
+
+      virtual void
+      get_intermediate_points_on_line (const line_iterator &line,
+                                      std::vector<Point<dim> > &points) const;
+
+      virtual void
+      get_intermediate_points_on_quad (const quad_iterator &quad,
+                                      std::vector<Point<dim> > &points) const;
+
+      virtual void
+      get_normals_at_vertices (const face_iterator &face,
+                              typename Boundary<dim>::FaceVertexNormals &face_vertex_normals) const;
+
+    private:
+      const double radius;
+          
+      void
+      get_intermediate_points_between_points (const Point<dim> &p0, const Point<dim> &p1,
+                                             std::vector<Point<dim> > &points) const;    
+  };
+
+  
+  template <>
+  Point<3>
+  CylinderBoundary<3>::
+  get_new_point_on_line (const line_iterator &line) const
+  {
+    const Point<3> middle = StraightBoundary<3>::get_new_point_on_line (line);
+                                    // project to boundary
+    Point<3> p(middle[0], middle[1], 0);
+    p *= radius/std::sqrt(p.square());
+
+    return Point<3> (p[0], p[1], middle[2]);
+  }
+
+
+  template<>
+  Point<3>
+  CylinderBoundary<3>::
+  get_new_point_on_quad (const quad_iterator &quad) const
+  {
+    Point<3> middle = StraightBoundary<3>::get_new_point_on_quad (quad);
+  
+                                    // project to boundary
+    Point<3> p(middle[0], middle[1], 0);
+    p *= radius/std::sqrt(p.square());
+
+    return Point<3> (p[0], p[1], middle[2]);
+  }
+
+
+  template <int dim>
+  void
+  CylinderBoundary<dim>::
+  get_intermediate_points_on_line (const line_iterator &line,
+                                  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> &,
+                                         const Point<dim> &,
+                                         std::vector<Point<dim> > &) const
+  {
+    Assert (false, ExcNotImplemented());
+  }
+
+
+  template <>
+  void
+  CylinderBoundary<3>::
+  get_intermediate_points_on_quad (const Triangulation<3>::quad_iterator &,
+                                  std::vector<Point<3> > &) const
+  {
+    Assert (false, ExcNotImplemented());
+  }
+
+
+  template <int dim>
+  void
+  CylinderBoundary<dim>::
+  get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &,
+                          typename Boundary<dim>::FaceVertexNormals &) const
+  {
+    Assert (false, ExcNotImplemented());
+  }
+
+
+  void
+  create_coarse_grid (Triangulation<3> &coarse_grid)
+  {
+    const unsigned int dim = 3;
+  
+    std::vector<Point<dim> >    vertices;
+    std::vector<CellData<dim> > cells;
+    SubCellData                 sub_cell_data;
+
+    const unsigned char
+      bottom_cylinder_boundary_id  = 10,
+      middle_cylinder_boundary_id  = 11,
+      spherical_boundary_id = 12,
+      all_other_boundary_id = 13,
+      straight_nondirichlet_boundary = 14;
+        
+        
+                                    // first build up the cells of the
+                                    // cylinder
+    {
+                                      // the vertices in each plane of
+                                      // the cylinder are located on
+                                      // three concentric rings of radii
+                                      // interior_hemisphere_radius,
+                                      // hemisphere_radius, and
+                                      // bottom_disk_radius,
+                                      // respectively. first generate
+                                      // these three rings
+      const Point<3> ring_points[8] = { Point<3>(-1,0,0),
+                                       Point<3>(-1,-1,0) / std::sqrt(2.),
+                                       Point<3>(0,-1,0),
+                                       Point<3>(+1,-1,0) / std::sqrt(2.),
+                                       Point<3>(+1,0,0),
+                                       Point<3>(+1,+1,0) / std::sqrt(2.),
+                                       Point<3>(0,+1,0),
+                                       Point<3>(-1,+1,0) / std::sqrt(2.) };
+
+                                      // first the point in the middle
+                                      // and the rest of those on the
+                                      // upper surface
+      vertices.push_back (Point<3>(0,0,bottom_disk_ceil));
+      for (unsigned int ring=0; ring<3; ++ring)
+       for (unsigned int i=0; i<8; ++i)
+         vertices.push_back (ring_points[i] * (ring == 0 ?
+                                               interior_hemisphere_radius :
+                                               (ring == 1 ? hemisphere_radius :
+                                                bottom_disk_radius))
+                             +
+                             Point<3>(0,0,bottom_disk_ceil));
+
+                                      // then points on lower surface
+      vertices.push_back (Point<3>(0,0,bottom_disk_floor));
+      for (unsigned int ring=0; ring<3; ++ring)
+       for (unsigned int i=0; i<8; ++i)
+         vertices.push_back (ring_points[i] * (ring == 0 ?
+                                               interior_hemisphere_radius :
+                                               (ring == 1 ?
+                                                hemisphere_radius :
+                                                bottom_disk_radius))
+                             +
+                             Point<3>(0,0,bottom_disk_floor));
+
+      const unsigned int n_vertices_per_surface = 25;
+      Assert (vertices.size() == n_vertices_per_surface*2,
+             ExcInternalError());
+    
+                                      // next create cells from these
+                                      // vertices. only store the
+                                      // vertices of the upper surface,
+                                      // the lower ones are the same
+                                      // +12
+      {
+       const unsigned int connectivity[20][4]
+         = { { 1, 2, 3, 0 },  // four cells in the center
+             { 3, 4, 5, 0 },
+             { 0, 5, 6, 7 },
+             { 1, 0, 7, 8 },            
+          
+             { 9, 10, 2, 1 },  // eight cells of inner ring
+             { 10, 11, 3, 2 },
+             { 11, 12, 4, 3 },
+             { 4, 12, 13, 5 },
+             { 5, 13, 14, 6 },
+             { 6, 14, 15, 7 },
+             { 8, 7, 15, 16 },
+             { 9, 1, 8, 16 },
+
+             { 17, 18, 10, 9 },  // eight cells of outer ring
+             { 18, 19, 11, 10 },
+             { 19, 20, 12, 11 },
+             { 12, 20, 21, 13 },
+             { 13, 21, 22, 14 },
+             { 14, 22, 23, 15 },
+             { 16, 15, 23, 24 },
+             { 17, 9, 16, 24 }   };
+
+                                        // now create cells out of this
+       for (unsigned int i=0; i<20; ++i)
+         {
+           CellData<3> cell;
+           for (unsigned int j=0; j<4; ++j)
+             {
+               cell.vertices[j]   = connectivity[i][j];
+               cell.vertices[j+4] = connectivity[i][j]+n_vertices_per_surface;
+             }
+           cell.material_id = 0;
+           cells.push_back (cell);
+         }
+      }
+    
+                                      // associate edges and faces on the
+                                      // outer boundary with boundary
+                                      // indicator of the cylinder
+                                      // boundary indicator. do this the
+                                      // same way as above, just this
+                                      // time with faces (edges follow
+                                      // from this immediately. some
+                                      // edges are duplicated since they
+                                      // belong to more than one cell,
+                                      // but that doesn't harm us here)
+      {
+       const unsigned int connectivity[8][2]
+         = { { 17,18 }, { 18, 19 }, { 19, 20 }, { 20, 21 },
+             { 21,22 }, { 22, 23 }, { 23, 24 }, { 24, 17 }};
+
+       for (unsigned int i=0; i<8; ++i)
+         {
+           const CellData<2> face = 
+             { { connectivity[i][0]+n_vertices_per_surface,
+                 connectivity[i][1]+n_vertices_per_surface,
+                 connectivity[i][1],
+                 connectivity[i][0] },
+               bottom_cylinder_boundary_id };
+           sub_cell_data.boundary_quads.push_back (face);
+
+           const CellData<1> edges[4] = 
+             { { { connectivity[i][0],    connectivity[i][1]    },
+                 bottom_cylinder_boundary_id },
+               { { connectivity[i][0]+n_vertices_per_surface,
+                   connectivity[i][1]+n_vertices_per_surface },
+                 bottom_cylinder_boundary_id },
+               { { connectivity[i][0]+n_vertices_per_surface,
+                   connectivity[i][0]    },
+                 bottom_cylinder_boundary_id },
+               { { connectivity[i][1]+n_vertices_per_surface,
+                   connectivity[i][1]    },
+                 bottom_cylinder_boundary_id } };
+           for (unsigned int i=0; i<4; ++i)
+             sub_cell_data.boundary_lines.push_back (edges[i]);
+         }
+      }
+    }
+
+                                    // next build up the middle ring. for
+                                    // this, copy the first 17 vertices
+                                    // up to z=0
+    {
+      const unsigned int first_upper_vertex = vertices.size();
+          
+      for (unsigned int i=0; i<17; ++i)
+       vertices.push_back (Point<3>(vertices[i][0], vertices[i][1], 0));
+
+                                      // next create cells from these
+                                      // vertices. only store the
+                                      // vertices of the lower surface,
+                                      // the lower ones are the same
+                                      // +first_upper_vertex
+      const unsigned int connectivity[12][4]
+       = { { 1, 2, 3, 0 },  // four cells in the center
+           { 3, 4, 5, 0 },
+           { 0, 5, 6, 7 },
+           { 1, 0, 7, 8 },            
+          
+           { 9, 10, 2, 1 },  // eight cells of ring
+           { 10, 11, 3, 2 },
+           { 11, 12, 4, 3 },
+           { 4, 12, 13, 5 },
+           { 5, 13, 14, 6 },
+           { 6, 14, 15, 7 },
+           { 8, 7, 15, 16 },
+           { 9, 1, 8, 16 }};
+                                      // now create cells out of this
+      for (unsigned int i=0; i<12; ++i)
        {
-         // Check if we have to create a cell
-         if ((*layers[zc])[xc][yc] == 1)
+         CellData<3> cell;
+         for (unsigned int j=0; j<4; ++j)
            {
-             const unsigned int z_vert = 25;
-             const unsigned int y_vert = 5;
-             unsigned int zoffs = zc * z_vert;
-             unsigned int yoffs = yc * y_vert;
-             unsigned int base_vert = zoffs + yoffs + xc;
-
-             CellData<dim> cell;
-             cell.vertices[0] = base_vert;
-             cell.vertices[1] = cell.vertices[0] + 1;
-             cell.vertices[2] = cell.vertices[0] + y_vert;
-             cell.vertices[3] = cell.vertices[1] + y_vert;
-             cell.vertices[4] = cell.vertices[0] + z_vert;
-             cell.vertices[5] = cell.vertices[1] + z_vert;
-             cell.vertices[6] = cell.vertices[2] + z_vert;
-             cell.vertices[7] = cell.vertices[3] + z_vert;
-             cell.material_id = 0;
-             cells.push_back (cell);
-
-             // Now add entries to the list of used
-             // vertices.
-             for (unsigned int i = 0; i < 8; ++i)
-               vertex_used[cell.vertices[i]] = true;
+             cell.vertices[j]   = connectivity[i][j]+first_upper_vertex;
+             cell.vertices[j+4] = connectivity[i][j];
            }
+         cell.material_id = 0;
+         cells.push_back (cell);
        }
 
-  // Now create vertices and renumber stuff;
-  std::vector<Point<dim> > vertices;
-  std::vector<unsigned int> vert_renumber (5*5*5, 0);
-  const double scale = 0.5;
-  unsigned int v_indx = 0;
-
-  for (int zv = 0; zv < 5; ++zv)
-    for (int yv = 0; yv < 5; ++yv)
-      for (int xv = 0; xv < 5; ++xv)
+                                      // mark the 8 vertical edges with
+                                      // the correct boundary indicator
+      for (unsigned int i=0; i<8; ++i)
        {
-         Point<dim> p_new ((double)(xv-2) * scale,
-                           (double)(yv-2) * scale,
-                           (double)(zv-2) * scale);
-
-         if (vertex_used[v_indx])
-           {
-             vert_renumber[v_indx] = vertices.size ();
-             vertices.push_back (p_new);
-           }
-         v_indx++;
+         const CellData<1> edge = { { 9, 9+first_upper_vertex },
+                                    middle_cylinder_boundary_id };
+         sub_cell_data.boundary_lines.push_back (edge);
+       }
+                                      // likewise with the 8 tangential
+                                      // edges on the lower disk. the
+                                      // edges at the interface between
+                                      // the middle disk and the
+                                      // hemisphere are handled by the
+                                      // hemisphere boundary
+      for (unsigned int i=0; i<8; ++i)
+       {
+         const CellData<1> edge = { { 9+i, 9+(i+1)%8},
+                                    middle_cylinder_boundary_id };
+         sub_cell_data.boundary_lines.push_back (edge);
        }
 
-  // Finally renumber the vertex indices in the cells
-  std::vector<CellData<dim> >::iterator cell_iterator;
-  for (cell_iterator = cells.begin (); cell_iterator != cells.end ();
-       ++cell_iterator)
+                                      // then assign face indicators
+      for (unsigned int i=0; i<8; ++i)
+       {
+         const CellData<2> face = { { 9+i,
+                                      9+(i+1)%8,
+                                      9+(i+1)%8+first_upper_vertex,
+                                      9+i+first_upper_vertex},
+                                    middle_cylinder_boundary_id };
+         sub_cell_data.boundary_quads.push_back (face);
+       }
+    }
+        
+                                    // the final part is setting the
+                                    // half-sphere on top of this
     {
-      for (unsigned int i = 0; i < 8; ++i)
-       cell_iterator->vertices[i] = 
-         vert_renumber[cell_iterator->vertices[i]];
+                                      // add four cubes to the top of
+                                      // the inner four cells, as well
+                                      // as 8 to their outside
+      {
+                                        // mirror the first nine vertices
+                                        // above the surface, and scale
+                                        // them to a certain distance
+                                        // outward
+       const double rx = hemisphere_radius / (1+std::sqrt(3.0));
+       for (unsigned int i=0; i<9; ++i)
+         {
+           const Point<3> p (vertices[i][0],
+                             vertices[i][1],
+                             i == 0 ?
+                             1
+                             :
+                             std::max(std::fabs(vertices[i][0]),
+                                      std::fabs(vertices[i][1])));
+           vertices.push_back (p / std::sqrt(p.square()) * rx);
+         }
+       Assert (vertices.size() == 76, ExcInternalError());
+
+                                        // same with the next ring of
+                                        // vertices, except that they
+                                        // go to hemisphere_radius
+       for (unsigned int i=9; i<17; ++i)
+         {
+           Point<3> p (vertices[i][0],
+                       vertices[i][1],
+                       std::max(std::fabs(vertices[i][0]),
+                                std::fabs(vertices[i][1])));
+           vertices.push_back (p / std::sqrt(p.square()) *
+                               hemisphere_radius);
+         }
+       Assert (vertices.size() == 84, ExcInternalError());
+      
+                                        // make 12 cells out of this
+       const unsigned int connectivity[12][4]
+         = { { 1, 2, 3, 0 },  // four cells in the center
+             { 3, 4, 5, 0 },
+             { 0, 5, 6, 7 },
+             { 1, 0, 7, 8 },
+
+             { 9, 10, 2, 1 },  // eight cells of inner ring
+             { 10, 11, 3, 2 },
+             { 11, 12, 4, 3 },
+             { 4, 12, 13, 5 },
+             { 5, 13, 14, 6 },
+             { 6, 14, 15, 7 },
+             { 8, 7, 15, 16 },
+             { 9, 1, 8, 16 },
+         };
+
+       for (unsigned int i=0; i<12; ++i)
+         {
+           CellData<3> cell;
+           for (unsigned int j=0; j<4; ++j)
+             {
+               cell.vertices[j]   = connectivity[i][j]+67;
+               cell.vertices[j+4] = connectivity[i][j]+50;
+             }
+           cell.material_id = 0;
+           cells.push_back (cell);
+         }
+      }
+
+                                      // assign boundary indicators to
+                                      // the faces and edges of these
+                                      // cells
+      {
+                                        // these are the numbers of the
+                                        // vertices on the top surface
+                                        // of the cylinder, with one
+                                        // "wrap-around":
+       const unsigned int vertices[9] = 
+         { 9, 10, 11, 12, 13, 14, 15, 16, 9 };
+                                        // their counter-parts are the
+                                        // same +67
+       for (unsigned int i=0; i<8; ++i)
+         {
+                                            // generate a face
+           const CellData<2> face = 
+             { { vertices[i]+50,   vertices[i+1]+50,
+                 vertices[i+1]+67, vertices[i]+67 },
+               spherical_boundary_id };
+           sub_cell_data.boundary_quads.push_back (face);
+
+                                            // same for the faces
+           const CellData<1> edges[4] =
+             { { { vertices[i]+50,   vertices[i+1]+50 },
+                 spherical_boundary_id },
+               { { vertices[i]+67,   vertices[i+1]+67 },
+                 spherical_boundary_id },
+               { { vertices[i]+50,   vertices[i]+67   },
+                 spherical_boundary_id },
+               { { vertices[i+1]+50, vertices[i+1]+67 },
+                 spherical_boundary_id } };
+           for (unsigned int j=0; j<4; ++j)
+             sub_cell_data.boundary_lines.push_back (edges[j]);
+         }
+      }  
+
+
+                                      // finally top the building
+                                      // with four closing cells and
+                                      // the vertex at the top
+      {
+       vertices.push_back (Point<3> (0,0,hemisphere_radius));
+
+       const unsigned int connectivity[4][8]
+         = { { 59, 60, 61, 67,   51, 52, 53, 50 },
+             { 61, 62, 63, 67,   53, 54, 55, 50 },
+             { 67, 63, 64, 65,   50, 55, 56, 57 },
+             { 59, 67, 65, 66,   51, 50, 57, 58 }};
+      
+       for (unsigned int i=0; i<4; ++i)
+         {
+           CellData<3> cell;
+           for (unsigned int j=0; j<8; ++j)
+             cell.vertices[j]   = connectivity[i][j]+17;
+           cell.material_id   = 0;
+           cells.push_back (cell);
+         }
+
+                                        // generate boundary
+                                        // information for these cells,
+                                        // too
+       for (unsigned int i=0; i<4; ++i)
+         {
+           const CellData<2> face = 
+             { { connectivity[i][0]+17, connectivity[i][1]+17,
+                 connectivity[i][2]+17, connectivity[i][3]+17 },
+               spherical_boundary_id };
+           sub_cell_data.boundary_quads.push_back (face);
+
+           const CellData<1> edges[4] =
+             { { { connectivity[i][0]+17, connectivity[i][1]+17 },
+                 spherical_boundary_id },
+               { { connectivity[i][1]+17, connectivity[i][2]+17 },
+                 spherical_boundary_id },
+               { { connectivity[i][2]+17, connectivity[i][3]+17 },
+                 spherical_boundary_id },
+               { { connectivity[i][3]+17, connectivity[i][0]+17 },
+                 spherical_boundary_id } };
+           for (unsigned int j=0; j<4; ++j)
+             sub_cell_data.boundary_lines.push_back (edges[j]);
+         }
+      }
     }
+  
 
-  // Now create triangulation
-  triangulation.create_triangulation (vertices,
-                                     cells,
-                                     SubCellData());
+                                    // finally generate a triangulation
+                                    // out of this
+    GridReordering<3>::reorder_cells (cells);
+    coarse_grid.create_triangulation_compatibility (vertices, cells,
+                                                   sub_cell_data);
+
+                                    // then associate boundary objects
+                                    // with the different boundary
+                                    // indicators
+    static const CylinderBoundary<3>
+      bottom_cylinder_boundary (bottom_disk_radius);
+    static const CylinderBoundary<3>
+      middle_cylinder_boundary (hemisphere_radius);
+    static const SphereBoundary<3> sphere_boundary;
+
+    coarse_grid.set_boundary (bottom_cylinder_boundary_id,
+                             bottom_cylinder_boundary);
+    coarse_grid.set_boundary (middle_cylinder_boundary_id,
+                             middle_cylinder_boundary);
+    coarse_grid.set_boundary (spherical_boundary_id,
+                             sphere_boundary);
+
+    for (Triangulation<dim>::active_cell_iterator cell=coarse_grid.begin_active();
+        cell != coarse_grid.end(); ++cell)
+      for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+       if ((cell->face(f)->boundary_indicator() == 0)
+           &&
+           (cell->face(f)->center()[2] >= (bottom_disk_floor+bottom_disk_ceil)/2))
+         cell->face(f)->set_boundary_indicator(straight_nondirichlet_boundary);
+  }
+}
+
+
+
+template <>
+void LaplaceProblem<3>::create_coarse_grid ()
+{
+  BreastPhantom::create_coarse_grid (triangulation);
 }
 
 
@@ -799,7 +1298,7 @@ int main ()
     {
       deallog.depth_console (0);
 
-      LaplaceProblem<2> laplace_problem;
+      LaplaceProblem<3> laplace_problem;
       laplace_problem.run ();
     }
   catch (std::exception &exc)

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.