From: bangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Sun, 1 Jul 2007 16:48:31 +0000 (+0000)
Subject: Remove some of the stuff we won't need any more for the tutorial program
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6a9daecff988525f5a67c22c253107f8b7915f17;p=dealii-svn.git

Remove some of the stuff we won't need any more for the tutorial program

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

diff --git a/deal.II/examples/step-27/step-27.cc b/deal.II/examples/step-27/step-27.cc
index 03c210f5c6..4c4d2d2319 100644
--- a/deal.II/examples/step-27/step-27.cc
+++ b/deal.II/examples/step-27/step-27.cc
@@ -1,10 +1,10 @@
 /* $Id$ */
-/* Author: Wolfgang Bangerth, University of Heidelberg, 2000 */
+/* Author: Wolfgang Bangerth, Texas A&M University, 2006, 2007 */
 
 /*    $Id$       */
 /*    Version: $Name$                                          */
 /*                                                                */
-/*    Copyright (C) 2000, 2001, 2002, 2003, 2004, 2006, 2007 by the deal.II authors */
+/*    Copyright (C) 2006, 2007 by the deal.II authors */
 /*                                                                */
 /*    This file is subject to QPL and may not be  distributed     */
 /*    without copyright and license information. Please refer     */
@@ -25,7 +25,7 @@
 #include <lac/vector.h>
 #include <lac/full_matrix.h>
 #include <lac/sparse_matrix.h>
-#include <lac/compressed_sparsity_pattern.h>
+#include <lac/compressed_set_sparsity_pattern.h>
 #include <lac/solver_cg.h>
 #include <lac/precondition.h>
 #include <grid/tria.h>
@@ -42,15 +42,15 @@
 #include <numerics/matrices.h>
 #include <numerics/data_out.h>
 
-#include <fstream>
-#include <iostream>
-
 #include <fe/fe_q.h>
 #include <grid/grid_out.h>
 #include <dofs/dof_constraints.h>
 #include <grid/grid_refinement.h>
 #include <numerics/error_estimator.h>
 
+#include <fstream>
+#include <iostream>
+
 #include <complex>
 
 
@@ -62,7 +62,7 @@ template <int dim>
 class LaplaceProblem 
 {
   public:
-    LaplaceProblem (const bool condense_glob = true);
+    LaplaceProblem ();
     ~LaplaceProblem ();
 
     void run ();
@@ -91,10 +91,7 @@ class LaplaceProblem
     Vector<double>       solution;
     Vector<double>       system_rhs;
 
-    Timer distr, condense, hang, assemble, solver;
-
     const unsigned int max_degree;
-    const bool condense_global;
 };
 
 
@@ -138,10 +135,10 @@ RightHandSide<dim>::value (const Point<dim>   &p,
 
 
 template <int dim>
-LaplaceProblem<dim>::LaplaceProblem (bool condense_glob) :
-  dof_handler (triangulation),
-  max_degree (dim == 2 ? 7 : 5),
-  condense_global (condense_glob)
+LaplaceProblem<dim>::LaplaceProblem ()
+		:
+		dof_handler (triangulation),
+		max_degree (dim == 2 ? 7 : 5)
 {
   for (unsigned int degree=2; degree<=max_degree; ++degree)
     {
@@ -161,61 +158,22 @@ LaplaceProblem<dim>::~LaplaceProblem ()
 template <int dim>
 void LaplaceProblem<dim>::setup_system ()
 {
-  distr.reset();
-  distr.start();
   dof_handler.distribute_dofs (fe_collection);
-  distr.stop();
 
   solution.reinit (dof_handler.n_dofs());
   system_rhs.reinit (dof_handler.n_dofs());
 
   hanging_node_constraints.clear ();
-  hang.reset();
-  hang.start();
   DoFTools::make_hanging_node_constraints (dof_handler,
 					   hanging_node_constraints);
 
   hanging_node_constraints.close ();
-  hang.stop();
 
-  if (dim < 3)
-    {
-      sparsity_pattern.reinit (dof_handler.n_dofs(),
-			       dof_handler.n_dofs(),
-			       dof_handler.max_couplings_between_dofs());
-      condense.reset();
-      if (condense_global)
-	{
-	  DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
-	  condense.start();
-	  hanging_node_constraints.condense (sparsity_pattern);
-	  condense.stop();
-	}
-      else
-	DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern,
-					 hanging_node_constraints);
-
-      sparsity_pattern.compress();      
-    }
-  else
-    {
-      CompressedSparsityPattern csp (dof_handler.n_dofs(),
-				     dof_handler.n_dofs());
-      condense.reset();
-      if (condense_global)
-	{
-	  DoFTools::make_sparsity_pattern (dof_handler, csp);
-
-	  condense.start();
-	  hanging_node_constraints.condense (csp);
-	  condense.stop();
-	}
-      else
-	DoFTools::make_sparsity_pattern (dof_handler, csp,
-					 hanging_node_constraints);
-	
-      sparsity_pattern.copy_from (csp);
-    }
+  CompressedSetSparsityPattern csp (dof_handler.n_dofs(),
+				    dof_handler.n_dofs());
+  DoFTools::make_sparsity_pattern (dof_handler, csp,
+				   hanging_node_constraints);
+  sparsity_pattern.copy_from (csp);
 
   system_matrix.reinit (sparsity_pattern);
 }
@@ -223,9 +181,6 @@ void LaplaceProblem<dim>::setup_system ()
 template <int dim>
 void LaplaceProblem<dim>::assemble_system () 
 {
-  assemble.reset ();
-  assemble.start ();
-  
   hp::FEValues<dim> hp_fe_values (fe_collection,
 				  quadrature_collection, 
 				  update_values    |  update_gradients |
@@ -270,40 +225,15 @@ void LaplaceProblem<dim>::assemble_system ()
 
       cell->get_dof_indices (local_dof_indices);
 
-      if (condense_global)
-	{
-	  for (unsigned int i=0; i<dofs_per_cell; ++i)
-	    {
-	      for (unsigned int j=0; j<dofs_per_cell; ++j)
-		system_matrix.add (local_dof_indices[i],
-				   local_dof_indices[j],
-				   cell_matrix(i,j));
-	      
-	      system_rhs(local_dof_indices[i]) += cell_rhs(i);
-	    }
-	}
-      else
-	{
-	  hanging_node_constraints
-	    .distribute_local_to_global (cell_matrix,
-					 local_dof_indices,
-					 system_matrix);
+      hanging_node_constraints
+	.distribute_local_to_global (cell_matrix,
+				     local_dof_indices,
+				     system_matrix);
 	  
-	  hanging_node_constraints
-	    .distribute_local_to_global (cell_rhs,
-					 local_dof_indices,
-					 system_rhs);
-	}
-    }
-
-  assemble.stop();
-  
-  if (condense_global)
-    {
-      condense.start();
-      hanging_node_constraints.condense (system_matrix);
-      hanging_node_constraints.condense (system_rhs);
-      condense.stop();
+      hanging_node_constraints
+	.distribute_local_to_global (cell_rhs,
+				     local_dof_indices,
+				     system_rhs);
     }
 
   std::map<unsigned int,double> boundary_values;
@@ -590,22 +520,13 @@ void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
     Vector<float> smoothness_indicators (triangulation.n_active_cells());
     estimate_smoothness (smoothness_indicators);
 
-    Vector<double> smoothness_field (dof_handler.n_dofs());
-    DoFTools::distribute_cell_to_dof_vector (dof_handler,
-					     smoothness_indicators,
-					     smoothness_field);
-
     Vector<float> fe_indices (triangulation.n_active_cells());
     {
       typename hp::DoFHandler<dim>::active_cell_iterator
 	cell = dof_handler.begin_active(),
 	endc = dof_handler.end();
       for (unsigned int index=0; cell!=endc; ++cell, ++index)
-	{
-	  
-	  fe_indices(index) = cell->active_fe_index();
-//	  smoothness_indicators(index) *= std::sqrt(cell->diameter());
-	}
+	fe_indices(index) = cell->active_fe_index();
     }
     
     const std::string filename = "solution-" +
@@ -615,8 +536,7 @@ void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
 
     data_out.attach_dof_handler (dof_handler);
     data_out.add_data_vector (solution, "solution");
-    data_out.add_data_vector (smoothness_indicators, "smoothness1");
-    data_out.add_data_vector (smoothness_field, "smoothness2");
+    data_out.add_data_vector (smoothness_indicators, "smoothness");
     data_out.add_data_vector (fe_indices, "fe_index");
     data_out.build_patches ();
   
@@ -698,586 +618,6 @@ void LaplaceProblem<2>::create_coarse_grid ()
 
 
 
-namespace BreastPhantom
-{
-  
-
-				   // 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,
-      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)
-	{
-	  CellData<3> cell;
-	  for (unsigned int j=0; j<4; ++j)
-	    {
-	      cell.vertices[j]   = connectivity[i][j]+first_upper_vertex;
-	      cell.vertices[j+4] = connectivity[i][j];
-	    }
-	  cell.material_id = 0;
-	  cells.push_back (cell);
-	}
-
-				       // mark the 8 vertical edges with
-				       // the correct boundary indicator
-      for (unsigned int i=0; i<8; ++i)
-	{
-	  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);
-	}
-
-				       // 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
-    {
-				       // 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]);
-	  }
-      }
-    }
-  
-
-				     // 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);
-}
-
-
 
 template <int dim>
 void LaplaceProblem<dim>::run () 
@@ -1296,9 +636,6 @@ void LaplaceProblem<dim>::run ()
 		<< triangulation.n_active_cells()
 		<< std::endl;
 
-      Timer all;
-      all.reset();
-      all.start();
       setup_system ();
 
       std::cout << "   Number of degrees of freedom: "
@@ -1311,20 +648,7 @@ void LaplaceProblem<dim>::run ()
       assemble_system ();
       
 
-      solver.reset();
-      solver.start();
       solve ();
-      solver.stop();
-      
-      all.stop();
-
-      std::cout << "   All: " << all()
-		<< ", distr: " << distr()
-		<< ", hang: " << hang()
-		<< ", condense: " << condense()
-		<< ", assemble: " << assemble()
-		<< ", solver: " << solver()
-		<< std::endl;
       
       output_results (cycle);
     }
@@ -1336,7 +660,7 @@ int main ()
     {
       deallog.depth_console (0);
 
-      LaplaceProblem<3> laplace_problem (false);
+      LaplaceProblem<2> laplace_problem;
       laplace_problem.run ();
     }
   catch (std::exception &exc)