]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Do a patch that has happened on mainline but failed on this branch because of local...
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 20 Sep 2012 03:09:52 +0000 (03:09 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 20 Sep 2012 03:09:52 +0000 (03:09 +0000)
git-svn-id: https://svn.dealii.org/branches/branch_bigger_global_dof_indices_3@26545 0785d39b-7218-0410-832d-ea1e28bc413d

tests/benchmarks/old/spec2006-447.dealII.cc
tests/benchmarks/spec2006-447.dealII.cc [deleted file]

index 9fc45eac2ce92e9f327ef84b0a133cca447a9fc8..73ca71506a65abd1e1a647fdfeab650fe3a453d9 100644 (file)
@@ -107,7 +107,7 @@ class MassCoefficient : public Function<dim>
       }
 };
 
-    
+
 
                                 // @sect3{Evaluating the solution}
 
@@ -128,16 +128,16 @@ namespace Evaluation
 {
                                   // @sect4{The EvaluationBase class}
   template <int dim>
-  class EvaluationBase 
+  class EvaluationBase
   {
     public:
       virtual ~EvaluationBase ();
 
       void set_refinement_cycle (const unsigned int refinement_cycle);
-      
+
       virtual void operator () (const DoFHandler<dim> &dof_handler,
                                const Vector<double>  &solution) const = 0;
-    
+
       unsigned int refinement_cycle;
   };
 
@@ -145,9 +145,9 @@ namespace Evaluation
   template <int dim>
   EvaluationBase<dim>::~EvaluationBase ()
   {}
-  
 
-  
+
+
   template <int dim>
   void
   EvaluationBase<dim>::set_refinement_cycle (const unsigned int step)
@@ -162,15 +162,15 @@ namespace Evaluation
   {
     public:
       PointValueEvaluation (const Point<dim>   &evaluation_point);
-      
+
       virtual void operator () (const DoFHandler<dim> &dof_handler,
                                const Vector<double>  &solution) const;
-      
+
       DeclException1 (ExcEvaluationPointNotFound,
                      Point<dim>,
                      << "The evaluation point " << arg1
                      << " was not found among the vertices of the present grid.");
-    
+
       const Point<dim>  evaluation_point;
   };
 
@@ -181,14 +181,14 @@ namespace Evaluation
                  :
                  evaluation_point (evaluation_point)
   {}
-  
+
 
 
   template <int dim>
   void
   PointValueEvaluation<dim>::
   operator () (const DoFHandler<dim> &dof_handler,
-              const Vector<double>  &solution) const 
+              const Vector<double>  &solution) const
   {
     double point_value = 1e20;
 
@@ -246,15 +246,15 @@ namespace Evaluation
   {
     public:
       PointXDerivativeEvaluation (const Point<dim>   &evaluation_point);
-      
+
       virtual void operator () (const DoFHandler<dim> &dof_handler,
                                const Vector<double>  &solution) const;
-      
+
       DeclException1 (ExcEvaluationPointNotFound,
                      Point<dim>,
                      << "The evaluation point " << arg1
                      << " was not found among the vertices of the present grid.");
-    
+
       const Point<dim>  evaluation_point;
   };
 
@@ -265,7 +265,7 @@ namespace Evaluation
                  :
                  evaluation_point (evaluation_point)
   {}
-  
+
 
                                   // The more interesting things
                                   // happen inside the function doing
@@ -274,7 +274,7 @@ namespace Evaluation
   void
   PointXDerivativeEvaluation<dim>::
   operator () (const DoFHandler<dim> &dof_handler,
-              const Vector<double>  &solution) const 
+              const Vector<double>  &solution) const
   {
                                     // This time initialize the
                                     // return value with something
@@ -294,7 +294,7 @@ namespace Evaluation
                             update_gradients | update_q_points);
     std::vector<Tensor<1,dim> >
       solution_gradients (vertex_quadrature.n_quadrature_points);
-    
+
                                     // ...and next loop over all cells
                                     // and their vertices, and count
                                     // how often the vertex has been
@@ -429,7 +429,7 @@ namespace Evaluation
   }
 
 
-  
+
                                   // @sect4{The GridOutput class}
 
                                   // Since this program has a more
@@ -460,10 +460,10 @@ namespace Evaluation
   {
     public:
       GridOutput (const std::string &output_name_base);
-      
+
       virtual void operator () (const DoFHandler<dim> &dof_handler,
                                const Vector<double>  &solution) const;
-    
+
       const std::string output_name_base;
   };
 
@@ -474,7 +474,7 @@ namespace Evaluation
                  :
                  output_name_base (output_name_base)
   {}
-  
+
 
   template <int dim>
   void
@@ -496,12 +496,12 @@ namespace Evaluation
 #else
     std::ofstream out (filename.str());
 #endif
-    
+
     GridOut().write_eps (dof_handler.get_tria(), out);
   }
 }
 
-  
+
                                 // @sect3{The Laplace solver classes}
 
                                 // Next are the actual solver
@@ -518,8 +518,8 @@ namespace LaplaceSolver
                                   // which requires the class to be
                                   // known:
   template <int dim> class WeightedResidual;
-  
-  
+
+
                                   // @sect4{The Laplace solver base class}
 
                                   // This class is almost unchanged,
@@ -552,8 +552,8 @@ namespace LaplaceSolver
       virtual void set_refinement_cycle (const unsigned int cycle);
 
       virtual void output_solution () const = 0;
-      
-    
+
+
       const SmartPointer<Triangulation<dim> > triangulation;
 
       unsigned int refinement_cycle;
@@ -568,7 +568,7 @@ namespace LaplaceSolver
 
 
   template <int dim>
-  Base<dim>::~Base () 
+  Base<dim>::~Base ()
   {}
 
 
@@ -579,7 +579,7 @@ namespace LaplaceSolver
   {
     refinement_cycle = cycle;
   }
-  
+
 
                                   // @sect4{The Laplace Solver class}
 
@@ -593,7 +593,7 @@ namespace LaplaceSolver
       Solver (Triangulation<dim>       &triangulation,
              const FiniteElement<dim> &fe,
              const Quadrature<dim>    &quadrature,
-             const Quadrature<dim-1>  &face_quadrature,              
+             const Quadrature<dim-1>  &face_quadrature,
              const Function<dim>      &boundary_values);
       virtual
       ~Solver ();
@@ -609,24 +609,24 @@ namespace LaplaceSolver
       virtual
       unsigned int
       n_dofs () const;
-      
-    
+
+
       const SmartPointer<const FiniteElement<dim> >  fe;
       const SmartPointer<const Quadrature<dim> >     quadrature;
-      const SmartPointer<const Quadrature<dim-1> >   face_quadrature;      
+      const SmartPointer<const Quadrature<dim-1> >   face_quadrature;
       DoFHandler<dim>                                dof_handler;
       Vector<double>                                 solution;
       const SmartPointer<const Function<dim> >       boundary_values;
 
       virtual void assemble_rhs (Vector<double> &rhs) const = 0;
-    
-    
+
+
       struct LinearSystem
       {
          LinearSystem (const DoFHandler<dim> &dof_handler);
 
          void solve (Vector<double> &solution) const;
-       
+
          ConstraintMatrix     hanging_node_constraints;
          SparsityPattern      sparsity_pattern;
          SparseMatrix<double> matrix;
@@ -655,14 +655,14 @@ namespace LaplaceSolver
                  Base<dim> (triangulation),
                  fe (&fe),
                   quadrature (&quadrature),
-                  face_quadrature (&face_quadrature),    
+                  face_quadrature (&face_quadrature),
                  dof_handler (triangulation),
                  boundary_values (&boundary_values)
   {}
 
 
   template <int dim>
-  Solver<dim>::~Solver () 
+  Solver<dim>::~Solver ()
   {
     dof_handler.clear ();
   }
@@ -696,7 +696,7 @@ namespace LaplaceSolver
   {
     return dof_handler.n_dofs();
   }
-  
+
 
   template <int dim>
   void
@@ -708,7 +708,7 @@ namespace LaplaceSolver
 
     const unsigned int n_threads = multithread_info.n_default_threads;
     std::vector<std::pair<active_cell_iterator,active_cell_iterator> >
-      thread_ranges 
+      thread_ranges
       = Threads::split_range<active_cell_iterator> (dof_handler.begin_active (),
                                                    dof_handler.end (),
                                                    n_threads);
@@ -730,7 +730,7 @@ namespace LaplaceSolver
                                              0,
                                              *boundary_values,
                                              boundary_value_map);
-    
+
     threads.join_all ();
     linear_system.hanging_node_constraints.condense (linear_system.matrix);
 
@@ -749,7 +749,7 @@ namespace LaplaceSolver
                                Threads::ThreadMutex                                 &mutex) const
   {
     MappingQ<dim> mapping (4);
-    FEValues<dim> fe_values (mapping, *fe, *quadrature, 
+    FEValues<dim> fe_values (mapping, *fe, *quadrature,
                             UpdateFlags(update_gradients | update_values |
                                         update_q_points |
                                         update_JxW_values));
@@ -759,12 +759,12 @@ namespace LaplaceSolver
 
     FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
 
-    std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
 
     std::vector<double> laplace_coefficients (fe_values.n_quadrature_points);
     std::vector<double> mass_coefficients (fe_values.n_quadrature_points);
-    
-    
+
+
     for (typename DoFHandler<dim>::active_cell_iterator cell=begin_cell;
         cell!=end_cell; ++cell)
       {
@@ -776,8 +776,8 @@ namespace LaplaceSolver
                                              laplace_coefficients);
        MassCoefficient<dim>().value_list (fe_values.get_quadrature_points(),
                                           mass_coefficients);
-       
-       
+
+
        for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
          for (unsigned int i=0; i<dofs_per_cell; ++i)
            for (unsigned int j=0; j<dofs_per_cell; ++j)
@@ -812,7 +812,7 @@ namespace LaplaceSolver
     void (*mhnc_p) (const DoFHandler<dim> &,
                    ConstraintMatrix      &)
       = &DoFTools::make_hanging_node_constraints;
-    
+
     Threads::Thread<>
       mhnc_thread = Threads::spawn (mhnc_p)(dof_handler, hanging_node_constraints);
 
@@ -928,17 +928,17 @@ namespace LaplaceSolver
 
       virtual
       void solve_problem ();
-      
+
       virtual
       unsigned int n_dofs () const;
-      
+
       virtual
       void postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const;
 
       virtual
       void output_solution () const;
-      
-    
+
+
       const SmartPointer<const Function<dim> > rhs_function;
       virtual void assemble_rhs (Vector<double> &rhs) const;
 
@@ -1008,16 +1008,16 @@ namespace LaplaceSolver
   {
     abort ();
   }
-  
+
 
 
   template <int dim>
   void
   PrimalSolver<dim>::
-  assemble_rhs (Vector<double> &rhs) const 
+  assemble_rhs (Vector<double> &rhs) const
   {
     MappingQ<dim> mapping (4);
-    FEValues<dim> fe_values (mapping, *this->fe, *this->quadrature, 
+    FEValues<dim> fe_values (mapping, *this->fe, *this->quadrature,
                             UpdateFlags(update_values    |
                                         update_q_points  |
                                         update_JxW_values));
@@ -1027,7 +1027,7 @@ namespace LaplaceSolver
 
     Vector<double>       cell_rhs (dofs_per_cell);
     std::vector<double>  rhs_values (n_q_points);
-    std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
 
     typename DoFHandler<dim>::active_cell_iterator
       cell = this->dof_handler.begin_active(),
@@ -1040,7 +1040,7 @@ namespace LaplaceSolver
 
        rhs_function->value_list (fe_values.get_quadrature_points(),
                                  rhs_values);
-      
+
        for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
          for (unsigned int i=0; i<dofs_per_cell; ++i)
            cell_rhs(i) += (fe_values.shape_value(i,q_point) *
@@ -1191,7 +1191,7 @@ namespace LaplaceSolver
 
       virtual void refine_grid ();
 
-    
+
       const SmartPointer<const Function<dim> > weighting_function;
   };
 
@@ -1247,7 +1247,7 @@ namespace LaplaceSolver
     for (unsigned int cell_index=0; cell!=endc; ++cell, ++cell_index)
       estimated_error(cell_index)
        *= weighting_function->value (cell->center());
-    
+
     GridRefinement::refine_and_coarsen_fixed_number (*this->triangulation,
                                                     estimated_error,
                                                     0.2, 0.02);
@@ -1392,7 +1392,7 @@ namespace LaplaceSolver
 namespace Data
 {
                                   // @sect4{The SetUpBase and SetUp classes}
-  
+
                                   // Based on the above description,
                                   // the ``SetUpBase'' class then
                                   // looks as follows. To allow using
@@ -1436,12 +1436,12 @@ namespace Data
 
       virtual
       const Function<dim> &  get_right_hand_side () const;
-      
+
 
       virtual
       void create_coarse_grid (Triangulation<dim> &coarse_grid) const;
 
-    
+
       static const typename Traits::BoundaryValues boundary_values;
       static const typename Traits::RightHandSide  right_hand_side;
   };
@@ -1458,7 +1458,7 @@ namespace Data
                                   // functions:
   template <class Traits, int dim>
   const Function<dim> &
-  SetUp<Traits,dim>::get_boundary_values () const 
+  SetUp<Traits,dim>::get_boundary_values () const
   {
     return boundary_values;
   }
@@ -1466,7 +1466,7 @@ namespace Data
 
   template <class Traits, int dim>
   const Function<dim> &
-  SetUp<Traits,dim>::get_right_hand_side () const 
+  SetUp<Traits,dim>::get_right_hand_side () const
   {
     return right_hand_side;
   }
@@ -1475,11 +1475,11 @@ namespace Data
   template <class Traits, int dim>
   void
   SetUp<Traits,dim>::
-  create_coarse_grid (Triangulation<dim> &coarse_grid) const 
+  create_coarse_grid (Triangulation<dim> &coarse_grid) const
   {
     Traits::create_coarse_grid (coarse_grid);
   }
-  
+
 
                                   // @sect4{The CurvedRidges class}
 
@@ -1496,7 +1496,7 @@ namespace Data
       {
        public:
          BoundaryValues () : Function<dim> () {};
-         
+
          virtual double value (const Point<dim>   &p,
                                const unsigned int  component) const;
       };
@@ -1506,7 +1506,7 @@ namespace Data
       {
        public:
          RightHandSide () : Function<dim> () {};
-         
+
          virtual double value (const Point<dim>   &p,
                                const unsigned int  component) const;
       };
@@ -1515,8 +1515,8 @@ namespace Data
       void
       create_coarse_grid (Triangulation<dim> &coarse_grid);
   };
-  
-    
+
+
   template <int dim>
   double
   CurvedRidges<dim>::BoundaryValues::
@@ -1553,7 +1553,7 @@ namespace Data
              100*std::sin(10*p(i)+5*p(0)*p(0));
       };
     t1 = t1*t1;
-    
+
     return -u*(t1+t2+t3);
   }
 
@@ -1566,10 +1566,10 @@ namespace Data
     GridGenerator::hyper_cube (coarse_grid, -1, 1);
     coarse_grid.refine_global (2);
   }
-  
+
 
                                   // @sect4{The Exercise_2_3 class}
-  
+
                                   // This example program was written
                                   // while giving practical courses
                                   // for a lecture on adaptive finite
@@ -1611,7 +1611,7 @@ namespace Data
        public:
          RightHandSide () : ConstantFunction<dim> (1.) {};
       };
-      
+
                                       // Finally a function to
                                       // generate the coarse
                                       // grid. This is somewhat more
@@ -1651,7 +1651,7 @@ namespace Data
                                   // enough to do the creation by
                                   // hand, rather than using a mesh
                                   // generator.
-/*  
+/*
   template <>
   void
   Exercise_2_3<2>::
@@ -1686,27 +1686,27 @@ namespace Data
             Point<2> (0.,    -1.),
             Point<2> (+1./2, -1.),
             Point<2> (+1,    -1.),
-            
+
             Point<2> (-1.,   -1./2.),
             Point<2> (-1./2, -1./2.),
             Point<2> (0.,    -1./2.),
             Point<2> (+1./2, -1./2.),
             Point<2> (+1,    -1./2.),
-            
+
             Point<2> (-1.,   0.),
             Point<2> (-1./2, 0.),
             Point<2> (+1./2, 0.),
             Point<2> (+1,    0.),
-            
+
             Point<2> (-1.,   1./2.),
             Point<2> (-1./2, 1./2.),
             Point<2> (0.,    1./2.),
             Point<2> (+1./2, 1./2.),
             Point<2> (+1,    1./2.),
-            
+
             Point<2> (-1.,   1.),
             Point<2> (-1./2, 1.),
-            Point<2> (0.,    1.),                        
+            Point<2> (0.,    1.),
             Point<2> (+1./2, 1.),
             Point<2> (+1,    1.)    };
     const unsigned int
@@ -1750,7 +1750,7 @@ namespace Data
                                     // the material indicator to zero
                                     // for all the cells:
     std::vector<CellData<dim> > cells (n_cells, CellData<dim>());
-    for (unsigned int i=0; i<n_cells; ++i) 
+    for (unsigned int i=0; i<n_cells; ++i)
       {
        for (unsigned int j=0;
             j<GeometryInfo<dim>::vertices_per_cell;
@@ -1772,14 +1772,14 @@ namespace Data
     coarse_grid.create_triangulation (vertices,
                                      cells,
                                      SubCellData());
-    
+
                                     // And since we want that the
                                     // evaluation point (3/4,3/4) in
                                     // this example is a grid point,
                                     // we refine once globally:
     coarse_grid.refine_global (1);
   }
-*/  
+*/
 
 
   template <>
@@ -1792,7 +1792,7 @@ namespace Data
     coarse_grid.set_boundary (0, boundary);
     coarse_grid.refine_global (1);
   }
-  
+
 }
 
                                 // @sect4{Discussion}
@@ -1908,7 +1908,7 @@ namespace Data
 namespace DualFunctional
 {
                                   // @sect4{The DualFunctionalBase class}
-  
+
                                   // First start with the base class
                                   // for dual functionals. Since for
                                   // linear problems the
@@ -1930,7 +1930,7 @@ namespace DualFunctional
 
 
                                   // @sect4{The PointValueEvaluation class}
-  
+
                                   // As a first application, we
                                   // consider the functional
                                   // corresponding to the evaluation
@@ -1952,13 +1952,13 @@ namespace DualFunctional
       void
       assemble_rhs (const DoFHandler<dim> &dof_handler,
                    Vector<double>        &rhs) const;
-      
+
       DeclException1 (ExcEvaluationPointNotFound,
                      Point<dim>,
                      << "The evaluation point " << arg1
                      << " was not found among the vertices of the present grid.");
 
-    
+
       const Point<dim> evaluation_point;
   };
 
@@ -1969,7 +1969,7 @@ namespace DualFunctional
                  :
                  evaluation_point (evaluation_point)
   {}
-  
+
 
                                   // As for doing the main purpose of
                                   // the class, assembling the right
@@ -2042,7 +2042,7 @@ namespace DualFunctional
 
 
                                   // @sect4{The PointXDerivativeEvaluation class}
-  
+
                                   // As second application, we again
                                   // consider the evaluation of the
                                   // x-derivative of the solution at
@@ -2061,13 +2061,13 @@ namespace DualFunctional
       void
       assemble_rhs (const DoFHandler<dim> &dof_handler,
                    Vector<double>        &rhs) const;
-      
+
       DeclException1 (ExcEvaluationPointNotFound,
                      Point<dim>,
                      << "The evaluation point " << arg1
                      << " was not found among the vertices of the present grid.");
 
-    
+
       const Point<dim> evaluation_point;
   };
 
@@ -2078,7 +2078,7 @@ namespace DualFunctional
                  :
                  evaluation_point (evaluation_point)
   {}
-  
+
 
                                   // What is interesting is the
                                   // implementation of this
@@ -2130,7 +2130,7 @@ namespace DualFunctional
                                     // for the number of quadrature
                                     // points and shape functions...
     QGauss<dim> quadrature(4);
-    MappingQ<dim> mapping (4);    
+    MappingQ<dim> mapping (4);
     FEValues<dim>  fe_values (mapping, dof_handler.get_fe(), quadrature,
                              update_gradients |
                              update_q_points  |
@@ -2146,7 +2146,7 @@ namespace DualFunctional
                                     // shape functions at the
                                     // quadrature points:
     Vector<double> cell_rhs (dofs_per_cell);
-    std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
 
                                     // Finally have a variable in
                                     // which we will sum up the
@@ -2155,7 +2155,7 @@ namespace DualFunctional
                                     // integrating the unit functions
                                     // on these cells:
     double total_volume = 0;
-    
+
                                     // Then start the loop over all
                                     // cells, and select those cells
                                     // which are close enough to the
@@ -2178,7 +2178,7 @@ namespace DualFunctional
                                           // total area/volume.
          fe_values.reinit (cell);
          cell_rhs = 0;
-         
+
          for (unsigned int q=0; q<n_q_points; ++q)
            {
              for (unsigned int i=0; i<dofs_per_cell; ++i)
@@ -2217,7 +2217,7 @@ namespace DualFunctional
                                     // which we have integrated:
     rhs.scale (1./total_volume);
   }
-  
+
 
 }
 
@@ -2276,7 +2276,7 @@ namespace LaplaceSolver
       virtual
       void
       solve_problem ();
-      
+
       virtual
       unsigned int
       n_dofs () const;
@@ -2285,7 +2285,7 @@ namespace LaplaceSolver
       void
       postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const;
 
-    
+
       const SmartPointer<const DualFunctional::DualFunctionalBase<dim> > dual_functional;
       virtual void assemble_rhs (Vector<double> &rhs) const;
 
@@ -2340,13 +2340,13 @@ namespace LaplaceSolver
   {
     Solver<dim>::postprocess(postprocessor);
   }
-  
+
 
 
   template <int dim>
   void
   DualSolver<dim>::
-  assemble_rhs (Vector<double> &rhs) const 
+  assemble_rhs (Vector<double> &rhs) const
   {
     dual_functional->assemble_rhs (this->dof_handler, rhs);
   }
@@ -2392,7 +2392,7 @@ namespace LaplaceSolver
       virtual
       void
       postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const;
-      
+
       virtual
       unsigned int
       n_dofs () const;
@@ -2403,7 +2403,7 @@ namespace LaplaceSolver
       void
       output_solution () const;
 
-    
+
                                       // In the private section, we
                                       // have two functions that are
                                       // used to call the
@@ -2420,7 +2420,7 @@ namespace LaplaceSolver
                                       // to avoid that we have to
                                       // write this lengthy name
                                       // over and over again:
-                                       
+
       typedef
       typename DoFHandler<dim>::active_cell_iterator
       active_cell_iterator;
@@ -2549,8 +2549,8 @@ namespace LaplaceSolver
          const SmartPointer<const Function<dim> > right_hand_side;
 
          std::vector<double> cell_residual;
-         std::vector<double> rhs_values;         
-         std::vector<double> dual_weights;       
+         std::vector<double> rhs_values;
+         std::vector<double> dual_weights;
          typename std::vector<Tensor<2,dim> > cell_grad_grads;
          CellData (const FiniteElement<dim> &fe,
                    const Quadrature<dim>    &quadrature,
@@ -2565,14 +2565,14 @@ namespace LaplaceSolver
          FESubfaceValues<dim> fe_subface_values_cell;
 
          std::vector<double> jump_residual;
-         std::vector<double> dual_weights;       
+         std::vector<double> dual_weights;
          typename std::vector<Tensor<1,dim> > cell_grads;
          typename std::vector<Tensor<1,dim> > neighbor_grads;
          FaceData (const FiniteElement<dim> &fe,
                    const Quadrature<dim-1>  &face_quadrature);
       };
 
-      
+
 
                                       // Regarding the evaluation of
                                       // the error estimator, we have
@@ -2652,17 +2652,17 @@ namespace LaplaceSolver
                             update_q_points           |
                             update_JxW_values),
                  right_hand_side (&right_hand_side)
-  {  
+  {
     const unsigned int n_q_points
       = quadrature.n_quadrature_points;
-  
+
     cell_residual.resize(n_q_points);
-    rhs_values.resize(n_q_points);    
-    dual_weights.resize(n_q_points);    
+    rhs_values.resize(n_q_points);
+    dual_weights.resize(n_q_points);
     cell_grad_grads.resize(n_q_points);
   }
-  
-  
+
+
 
   template <int dim>
   WeightedResidual<dim>::FaceData::
@@ -2682,16 +2682,16 @@ namespace LaplaceSolver
                                           update_normal_vectors),
                  fe_subface_values_cell (mapping, fe, face_quadrature,
                                          update_gradients)
-  {  
+  {
     const unsigned int n_face_q_points
       = face_quadrature.n_quadrature_points;
-  
+
     jump_residual.resize(n_face_q_points);
-    dual_weights.resize(n_face_q_points);    
+    dual_weights.resize(n_face_q_points);
     cell_grads.resize(n_face_q_points);
     neighbor_grads.resize(n_face_q_points);
   }
-  
+
 
 
 
@@ -2735,7 +2735,7 @@ namespace LaplaceSolver
     threads.join_all ();
   }
 
-  
+
   template <int dim>
   void
   WeightedResidual<dim>::solve_primal_problem ()
@@ -2749,7 +2749,7 @@ namespace LaplaceSolver
   {
     DualSolver<dim>::solve_problem ();
   }
-  
+
 
   template <int dim>
   void
@@ -2758,8 +2758,8 @@ namespace LaplaceSolver
   {
     PrimalSolver<dim>::postprocess (postprocessor);
   }
-  
-  
+
+
   template <int dim>
   unsigned int
   WeightedResidual<dim>::n_dofs () const
@@ -2809,11 +2809,11 @@ namespace LaplaceSolver
                                     // of the error.
     GridRefinement::refine_and_coarsen_fixed_fraction (*this->triangulation,
                                                       error_indicators,
-                                                      0.5/std::sqrt(1.+step), 
+                                                      0.5/std::sqrt(1.+step),
                                                       0.2/std::sqrt(1.+step));
     this->triangulation->execute_coarsening_and_refinement ();
   }
-  
+
 
                                   // Since we want to output both the
                                   // primal and the dual solution, we
@@ -2905,7 +2905,7 @@ namespace LaplaceSolver
   {
     const PrimalSolver<dim> &primal_solver = *this;
     const DualSolver<dim>   &dual_solver   = *this;
-    
+
     ConstraintMatrix primal_hanging_node_constraints;
     DoFTools::make_hanging_node_constraints (primal_solver.dof_handler,
                                             primal_hanging_node_constraints);
@@ -2915,7 +2915,7 @@ namespace LaplaceSolver
                          dual_solver.solution,
                          primal_solver.dof_handler,
                          primal_hanging_node_constraints,
-                         dual_solution);    
+                         dual_solution);
 
                                     // approximate error, gradient,
                                     // and second derivative
@@ -2936,13 +2936,13 @@ namespace LaplaceSolver
                              primal_solver.dof_handler,
                              primal_solver.solution,
                              gradient_indicators);
-      
+
       DerivativeApproximation::
        approximate_second_derivative (mapping,
                                       primal_solver.dof_handler,
                                       primal_solver.solution,
                                       second_indicators);
-      
+
     }
                                     // distribute cell to dof vectors
     Vector<double> x_error_indicators (primal_solver.dof_handler.n_dofs());
@@ -2957,8 +2957,8 @@ namespace LaplaceSolver
     DoFTools::distribute_cell_to_dof_vector (primal_solver.dof_handler,
                                             second_indicators,
                                             x_second_indicators);
-    
-    
+
+
 
                                     // we generate too much output in
                                     // 3d. instead of doing it that
@@ -2996,20 +2996,20 @@ namespace LaplaceSolver
       {
        coarse_to_fine_map[cell]->get_interpolated_dof_values (primal_solver.solution,tmp);
        cell->set_dof_values (tmp, coarse_primal_solution);
-       
+
        coarse_to_fine_map[cell]->get_interpolated_dof_values (dual_solution,tmp);
        cell->set_dof_values (tmp, coarse_dual_solution);
-       
+
        coarse_to_fine_map[cell]->get_interpolated_dof_values (x_error_indicators,tmp);
        cell->set_dof_values (tmp, coarse_error_indicators);
-       
+
        coarse_to_fine_map[cell]->get_interpolated_dof_values (x_gradient_indicators,tmp);
        cell->set_dof_values (tmp, coarse_gradient_indicators);
-       
+
        coarse_to_fine_map[cell]->get_interpolated_dof_values (x_second_indicators,tmp);
        cell->set_dof_values (tmp, coarse_second_indicators);
       }
-       
+
     {
       DataOut<dim> data_out;
       data_out.attach_dof_handler (coarser_dof_handler);
@@ -3019,7 +3019,7 @@ namespace LaplaceSolver
       data_out.add_data_vector (coarse_gradient_indicators, "gradient");
       data_out.add_data_vector (coarse_second_indicators, "second_derivatives");
       data_out.build_patches ();
-  
+
 #ifdef HAVE_STD_STRINGSTREAM
       std::ostringstream filename;
 #else
@@ -3035,10 +3035,10 @@ namespace LaplaceSolver
 #else
       std::ofstream out (filename.str());
 #endif
-    
+
       data_out.write_gmv (out);
     }
-    
+
   }
 
 
@@ -3053,7 +3053,7 @@ namespace LaplaceSolver
                                   // functions that actually do the
                                   // work, and finally collects the
                                   // results.
-  
+
   template <int dim>
   void
   WeightedResidual<dim>::
@@ -3105,7 +3105,7 @@ namespace LaplaceSolver
                          dual_solver.dof_handler,
                          dual_hanging_node_constraints,
                          primal_solution);
-    
+
                                     // Then for computing the
                                     // interpolation of the
                                     // numerically approximated dual
@@ -3128,7 +3128,7 @@ namespace LaplaceSolver
                                       primal_solver.dof_handler,
                                       primal_hanging_node_constraints,
                                       dual_weights);
-    
+
                                     // Note that this could probably
                                     // have been more efficient since
                                     // those constraints have been
@@ -3139,7 +3139,7 @@ namespace LaplaceSolver
                                     // dual solution. We leave the
                                     // optimization of the program in
                                     // this respect as an exercise.
-    
+
                                     // Having computed the dual
                                     // weights we now proceed with
                                     // computing the cell and face
@@ -3206,7 +3206,7 @@ namespace LaplaceSolver
                   n_threads, i,
                   error_indicators,
                   face_integrals);
-    threads.join_all();    
+    threads.join_all();
 
                                     // Once the error contributions
                                     // are computed, sum them up. For
@@ -3222,7 +3222,7 @@ namespace LaplaceSolver
                                     // term, since the other half
                                     // will be taken by the
                                     // neighboring cell.
-    unsigned int present_cell=0;  
+    unsigned int present_cell=0;
     for (active_cell_iterator cell=dual_solver.dof_handler.begin_active();
         cell!=dual_solver.dof_handler.end();
         ++cell, ++present_cell)
@@ -3284,7 +3284,7 @@ namespace LaplaceSolver
                        *dual_solver.quadrature,
                        *primal_solver.rhs_function);
     FaceData face_data (*dual_solver.fe,
-                       *dual_solver.face_quadrature);    
+                       *dual_solver.face_quadrature);
 
                                     // Then calculate the start cell
                                     // for this thread. We let the
@@ -3320,7 +3320,7 @@ namespace LaplaceSolver
                                     // back right now
     if (cell == dual_solver.dof_handler.end())
       return;
-    
+
                                     // Next loop over all cells. The
                                     // check for loop end is done at
                                     // the end of the loop, along
@@ -3340,7 +3340,7 @@ namespace LaplaceSolver
                             dual_weights,
                             cell_data,
                             error_indicators);
-       
+
                                         // After computing the cell
                                         // terms, turn to the face
                                         // terms. For this, loop over
@@ -3364,12 +3364,12 @@ namespace LaplaceSolver
                                             // the list of faces with
                                             // a zero contribution to
                                             // the error.
-           if (cell->face(face_no)->at_boundary()) 
+           if (cell->face(face_no)->at_boundary())
              {
                face_integrals[cell->face(face_no)] = 0;
                continue;
              };
-           
+
                                             // Next, note that since
                                             // we want to compute the
                                             // jump terms on each
@@ -3413,7 +3413,7 @@ namespace LaplaceSolver
                                             // the coarse cell.
            if (cell->at_boundary(face_no) == false)
              if (cell->neighbor(face_no)->level() < cell->level())
-               continue;         
+               continue;
 
 
                                             // Now we know that we
@@ -3443,7 +3443,7 @@ namespace LaplaceSolver
                                           primal_solution,
                                           dual_weights,
                                           face_data,
-                                          face_integrals);       
+                                          face_integrals);
            else
              integrate_over_irregular_face (cell, face_no,
                                             primal_solution,
@@ -3518,7 +3518,7 @@ namespace LaplaceSolver
 
 
                                   // @sect4{Computing edge term error contributions - 1}
-  
+
                                   // On the other hand, computation
                                   // of the edge terms for the error
                                   // estimate is not so
@@ -3595,7 +3595,7 @@ namespace LaplaceSolver
                                     // that cell, and extract the
                                     // gradients on that cell:
     const active_cell_iterator neighbor = cell->neighbor(face_no);
-    face_data.fe_face_values_neighbor.reinit (neighbor, neighbor_neighbor);      
+    face_data.fe_face_values_neighbor.reinit (neighbor, neighbor_neighbor);
     face_data.fe_face_values_neighbor.get_function_grads (primal_solution,
                                                          face_data.neighbor_grads);
 
@@ -3614,7 +3614,7 @@ namespace LaplaceSolver
                                     // this face:
     face_data.fe_face_values_cell.get_function_values (dual_weights,
                                                       face_data.dual_weights);
-    
+
                                     // Finally, we have to compute
                                     // the sum over jump residuals,
                                     // dual weights, and quadrature
@@ -3657,7 +3657,7 @@ namespace LaplaceSolver
 
 
                                   // @sect4{Computing edge term error contributions - 2}
-  
+
                                   // We are still missing the case of
                                   // faces with hanging nodes. This
                                   // is what is covered in this
@@ -3680,7 +3680,7 @@ namespace LaplaceSolver
       n_q_points = face_data.fe_face_values_cell.n_quadrature_points;
 
     const typename DoFHandler<dim>::cell_iterator
-      neighbor = cell->neighbor(face_no);    
+      neighbor = cell->neighbor(face_no);
     Assert (neighbor.state() == IteratorState::valid,
            ExcInternalError());
     Assert (neighbor->has_children(),
@@ -3697,7 +3697,7 @@ namespace LaplaceSolver
                                     // the same.
     const unsigned int
       neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
-  
+
                                     // Then simply do everything we
                                     // did in the previous function
                                     // for one face for all the
@@ -3722,7 +3722,7 @@ namespace LaplaceSolver
                                             neighbor_neighbor);
        face_data.fe_face_values_neighbor.get_function_grads (primal_solution,
                                                              face_data.neighbor_grads);
-      
+
                                         // and finally building the
                                         // jump residuals. Since we
                                         // take the normal vector
@@ -3738,7 +3738,7 @@ namespace LaplaceSolver
                                         // Then get dual weights:
        face_data.fe_face_values_neighbor.get_function_values (dual_weights,
                                                               face_data.dual_weights);
-       
+
                                         // At last, sum up the
                                         // contribution of this
                                         // sub-face, and set it in
@@ -3767,21 +3767,21 @@ namespace LaplaceSolver
     typename DoFHandler<dim>::face_iterator face = cell->face(face_no);
     for (unsigned int subface_no=0;
         subface_no<GeometryInfo<dim>::subfaces_per_face;
-        ++subface_no) 
+        ++subface_no)
       {
        Assert (face_integrals.find(face->child(subface_no)) !=
                face_integrals.end(),
                ExcInternalError());
        Assert (face_integrals[face->child(subface_no)] != -1e20,
                ExcInternalError());
-      
+
        sum += face_integrals[face->child(subface_no)];
       };
                                     // Finally store the value with
                                     // the parent face.
     face_integrals[face] = sum;
   }
-  
+
 }
 
 
@@ -3837,7 +3837,7 @@ struct Framework
                                     // structure, these values are
                                     // all set to default values, for
                                     // simple use.
-    struct ProblemDescription 
+    struct ProblemDescription
     {
                                         // First allow for the
                                         // degrees of the piecewise
@@ -4006,7 +4006,7 @@ void Framework<dim>::run (const ProblemDescription &descriptor)
                                                      *descriptor.dual_functional);
        break;
       };
-       
+
       case ProblemDescription::global_refinement:
       {
        solver
@@ -4018,7 +4018,7 @@ void Framework<dim>::run (const ProblemDescription &descriptor)
                                                      descriptor.data->get_boundary_values());
        break;
       };
-       
+
       case ProblemDescription::kelly_indicator:
       {
        solver
@@ -4043,11 +4043,11 @@ void Framework<dim>::run (const ProblemDescription &descriptor)
                                                             *descriptor.kelly_weight);
        break;
       };
-           
+
       default:
            AssertThrow (false, ExcInternalError());
     };
-  
+
                                   // Now that all objects are in
                                   // place, run the main loop. The
                                   // stopping criterion is
@@ -4066,14 +4066,14 @@ void Framework<dim>::run (const ProblemDescription &descriptor)
     {
       std::cout << "Refinement cycle: "        << step
                << std::endl;
-           
+
       solver->set_refinement_cycle (step);
       solver->solve_problem ();
       solver->output_solution ();
 
       std::cout << "   Number of degrees of freedom="
                << solver->n_dofs() << std::endl;
-      
+
       for (typename EvaluatorList::const_iterator
             e = descriptor.evaluator_list.begin();
           e != descriptor.evaluator_list.end(); ++e)
@@ -4082,7 +4082,7 @@ void Framework<dim>::run (const ProblemDescription &descriptor)
          solver->postprocess (**e);
        };
 
-           
+
       if (solver->n_dofs() < descriptor.max_degrees_of_freedom)
        solver->refine_grid ();
       else
@@ -4111,7 +4111,7 @@ void Framework<dim>::run (const ProblemDescription &descriptor)
                                 // etc), and passes them packed into
                                 // a structure to the frame work
                                 // class above.
-int main (int argc, char **argv) 
+int main (int argc, char **argv)
 {
                                   // if no argument is given, then do 18
                                   // iterations
@@ -4134,10 +4134,10 @@ int main (int argc, char **argv)
                  << std::endl;
        exit (1);
       }
-  
-       
-  
-  
+
+
+
+
   deallog.depth_console (0);
   try
     {
@@ -4183,7 +4183,7 @@ int main (int argc, char **argv)
                                       // can also use
                                       // ``CurvedRidges<dim>'':
       descriptor.data = new Data::SetUp<Data::Exercise_2_3<dim>,dim> ();
-      
+
                                       // Next set first a dual
                                       // functional, then a list of
                                       // evaluation objects. We
@@ -4216,12 +4216,12 @@ int main (int argc, char **argv)
       const Point<dim> evaluation_point (0., 0., 0.);
       descriptor.dual_functional
        = new DualFunctional::PointValueEvaluation<dim> (evaluation_point);
-      
+
       Evaluation::PointValueEvaluation<dim>
        postprocessor1 (evaluation_point);
       Evaluation::GridOutput<dim>
        postprocessor2 ("grid");
-      
+
       descriptor.evaluator_list.push_back (&postprocessor1);
       descriptor.evaluator_list.push_back (&postprocessor2);
 
@@ -4236,7 +4236,7 @@ int main (int argc, char **argv)
 #else
       descriptor.max_degrees_of_freedom = 20000;
 #endif
-      
+
                                       // Finally pass the descriptor
                                       // object to a function that
                                       // runs the entire solution
@@ -4259,7 +4259,7 @@ int main (int argc, char **argv)
                << std::endl;
       return 1;
     }
-  catch (...) 
+  catch (...)
     {
       std::cerr << std::endl << std::endl
                << "----------------------------------------------------"
diff --git a/tests/benchmarks/spec2006-447.dealII.cc b/tests/benchmarks/spec2006-447.dealII.cc
deleted file mode 100644 (file)
index 7801034..0000000
+++ /dev/null
@@ -1,4275 +0,0 @@
-/* $Id$ */
-/* Author: Wolfgang Bangerth, ETH Zurich, 2002 */
-
-/*    $Id$       */
-/*    Version: $Name$                                          */
-/*                                                                */
-/*    Copyright (C) 2002, 2003, 2004, 2008, 2010 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.                        */
-
-
-                                // Start out with well known things...
-#include <deal.II/base/quadrature_lib.h>
-#include <deal.II/base/function.h>
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/thread_management.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/lac/compressed_sparsity_pattern.h>
-#include <deal.II/lac/sparse_matrix.h>
-#include <deal.II/lac/solver_cg.h>
-#include <deal.II/lac/vector_memory.h>
-#include <deal.II/lac/precondition.h>
-#include <deal.II/grid/tria.h>
-#include <deal.II/grid/grid_generator.h>
-#include <deal.II/grid/grid_out.h>
-#include <deal.II/grid/tria_accessor.h>
-#include <deal.II/grid/tria_iterator.h>
-#include <deal.II/grid/grid_refinement.h>
-#include <deal.II/grid/tria_boundary_lib.h>
-#include <deal.II/grid/intergrid_map.h>
-#include <deal.II/dofs/dof_handler.h>
-#include <deal.II/lac/constraint_matrix.h>
-#include <deal.II/dofs/dof_accessor.h>
-#include <deal.II/dofs/dof_tools.h>
-#include <deal.II/fe/fe_q.h>
-#include <deal.II/fe/fe_values.h>
-#include <deal.II/fe/fe_tools.h>
-#include <deal.II/fe/mapping_q.h>
-#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/numerics/matrix_tools.h>
-#include <deal.II/numerics/data_out.h>
-#include <deal.II/numerics/error_estimator.h>
-#include <deal.II/numerics/derivative_approximation.h>
-
-#include <iostream>
-#include <fstream>
-#include <list>
-#include <algorithm>
-#include <numeric>
-#include <stdlib.h>
-
-#ifdef HAVE_STD_STRINGSTREAM
-#  include <sstream>
-#else
-#  include <strstream>
-#endif
-
-using namespace dealii;
-
-
-/**
- * The number of iterations the program shall do. This is given as a
- * parameter to the executable.
- */
-unsigned int n_steps = 0;
-/**
- * The present step.
- */
-unsigned int step;
-
-                                /**
-                                 * Declare the coefficient in front
-                                 * of the Laplace operator of the
-                                 * equation.
-                                 */
-template <int dim>
-class LaplaceCoefficient : public Function<dim>
-{
-  public:
-    virtual double value (const Point<dim> &p,
-                         const unsigned int) const
-      {
-       return 1.+p*p;
-      }
-};
-
-
-                                /**
-                                 * Same for the coefficient in front
-                                 * of the mass term.
-                                 */
-template <int dim>
-class MassCoefficient : public Function<dim>
-{
-  public:
-    virtual double value (const Point<dim> &p,
-                         const unsigned int) const
-      {
-       double factor = 1;
-       for (unsigned int d=0; d<dim; ++d)
-         factor *= std::sin(4*p(d));
-       return 1.+factor/2;
-      }
-};
-
-    
-
-                                // @sect3{Evaluating the solution}
-
-                                // As mentioned in the introduction,
-                                // significant parts of the program
-                                // have simply been taken over from
-                                // the step-13 example program. We
-                                // therefore only comment on those
-                                // things that are new.
-                                //
-                                // First, the framework for
-                                // evaluation of solutions is
-                                // unchanged, i.e. the base class is
-                                // the same, and the class to
-                                // evaluate the solution at a grid
-                                // point is unchanged:
-namespace Evaluation
-{
-                                  // @sect4{The EvaluationBase class}
-  template <int dim>
-  class EvaluationBase 
-  {
-    public:
-      virtual ~EvaluationBase ();
-
-      void set_refinement_cycle (const unsigned int refinement_cycle);
-      
-      virtual void operator () (const DoFHandler<dim> &dof_handler,
-                               const Vector<double>  &solution) const = 0;
-    
-      unsigned int refinement_cycle;
-  };
-
-
-  template <int dim>
-  EvaluationBase<dim>::~EvaluationBase ()
-  {}
-  
-
-  
-  template <int dim>
-  void
-  EvaluationBase<dim>::set_refinement_cycle (const unsigned int step)
-  {
-    refinement_cycle = step;
-  }
-
-
-                                  // @sect4{The PointValueEvaluation class}
-  template <int dim>
-  class PointValueEvaluation : public EvaluationBase<dim>
-  {
-    public:
-      PointValueEvaluation (const Point<dim>   &evaluation_point);
-      
-      virtual void operator () (const DoFHandler<dim> &dof_handler,
-                               const Vector<double>  &solution) const;
-      
-      DeclException1 (ExcEvaluationPointNotFound,
-                     Point<dim>,
-                     << "The evaluation point " << arg1
-                     << " was not found among the vertices of the present grid.");
-    
-      const Point<dim>  evaluation_point;
-  };
-
-
-  template <int dim>
-  PointValueEvaluation<dim>::
-  PointValueEvaluation (const Point<dim>   &evaluation_point)
-                 :
-                 evaluation_point (evaluation_point)
-  {}
-  
-
-
-  template <int dim>
-  void
-  PointValueEvaluation<dim>::
-  operator () (const DoFHandler<dim> &dof_handler,
-              const Vector<double>  &solution) const 
-  {
-    double point_value = 1e20;
-
-    typename DoFHandler<dim>::active_cell_iterator
-      cell = dof_handler.begin_active(),
-      endc = dof_handler.end();
-    bool evaluation_point_found = false;
-    for (; (cell!=endc) && !evaluation_point_found; ++cell)
-      for (unsigned int vertex=0;
-          vertex<GeometryInfo<dim>::vertices_per_cell;
-          ++vertex)
-       if (cell->vertex(vertex).distance (evaluation_point)
-           <
-           cell->diameter() * 1e-8)
-         {
-           point_value = solution(cell->vertex_dof_index(vertex,0));
-
-           evaluation_point_found = true;
-           break;
-         };
-
-    AssertThrow (evaluation_point_found,
-                ExcEvaluationPointNotFound(evaluation_point));
-
-    std::cout << "   Point value=" << point_value
-             << std::endl;
-  }
-
-
-                                  // @sect4{The PointXDerivativeEvaluation class}
-
-                                  // Besides the class implementing
-                                  // the evaluation of the solution
-                                  // at one point, we here provide
-                                  // one which evaluates the gradient
-                                  // at a grid point. Since in
-                                  // general the gradient of a finite
-                                  // element function is not
-                                  // continuous at a vertex, we have
-                                  // to be a little bit more careful
-                                  // here. What we do is to loop over
-                                  // all cells, even if we have found
-                                  // the point already on one cell,
-                                  // and use the mean value of the
-                                  // gradient at the vertex taken
-                                  // from all adjacent cells.
-                                  //
-                                  // Given the interface of the
-                                  // ``PointValueEvaluation'' class,
-                                  // the declaration of this class
-                                  // provides little surprise, and
-                                  // neither does the constructor:
-  template <int dim>
-  class PointXDerivativeEvaluation : public EvaluationBase<dim>
-  {
-    public:
-      PointXDerivativeEvaluation (const Point<dim>   &evaluation_point);
-      
-      virtual void operator () (const DoFHandler<dim> &dof_handler,
-                               const Vector<double>  &solution) const;
-      
-      DeclException1 (ExcEvaluationPointNotFound,
-                     Point<dim>,
-                     << "The evaluation point " << arg1
-                     << " was not found among the vertices of the present grid.");
-    
-      const Point<dim>  evaluation_point;
-  };
-
-
-  template <int dim>
-  PointXDerivativeEvaluation<dim>::
-  PointXDerivativeEvaluation (const Point<dim>   &evaluation_point)
-                 :
-                 evaluation_point (evaluation_point)
-  {}
-  
-
-                                  // The more interesting things
-                                  // happen inside the function doing
-                                  // the actual evaluation:
-  template <int dim>
-  void
-  PointXDerivativeEvaluation<dim>::
-  operator () (const DoFHandler<dim> &dof_handler,
-              const Vector<double>  &solution) const 
-  {
-                                    // This time initialize the
-                                    // return value with something
-                                    // useful, since we will have to
-                                    // add up a number of
-                                    // contributions and take the
-                                    // mean value afterwards...
-    double point_derivative = 0;
-
-                                    // ...then have some objects of
-                                    // which the meaning wil become
-                                    // clear below...
-    QTrapez<dim>  vertex_quadrature;
-    MappingQ<dim> mapping (4);
-    FEValues<dim> fe_values (mapping, dof_handler.get_fe(),
-                            vertex_quadrature,
-                            update_gradients | update_q_points);
-    std::vector<Tensor<1,dim> >
-      solution_gradients (vertex_quadrature.n_quadrature_points);
-    
-                                    // ...and next loop over all cells
-                                    // and their vertices, and count
-                                    // how often the vertex has been
-                                    // found:
-    typename DoFHandler<dim>::active_cell_iterator
-      cell = dof_handler.begin_active(),
-      endc = dof_handler.end();
-    unsigned int evaluation_point_hits = 0;
-    for (; cell!=endc; ++cell)
-      for (unsigned int vertex=0;
-          vertex<GeometryInfo<dim>::vertices_per_cell;
-          ++vertex)
-       if (cell->vertex(vertex) == evaluation_point)
-         {
-                                            // Things are now no more
-                                            // as simple, since we
-                                            // can't get the gradient
-                                            // of the finite element
-                                            // field as before, where
-                                            // we simply had to pick
-                                            // one degree of freedom
-                                            // at a vertex.
-                                            //
-                                            // Rather, we have to
-                                            // evaluate the finite
-                                            // element field on this
-                                            // cell, and at a certain
-                                            // point. As you know,
-                                            // evaluating finite
-                                            // element fields at
-                                            // certain points is done
-                                            // through the
-                                            // ``FEValues'' class, so
-                                            // we use that. The
-                                            // question is: the
-                                            // ``FEValues'' object
-                                            // needs to be a given a
-                                            // quadrature formula and
-                                            // can then compute the
-                                            // values of finite
-                                            // element quantities at
-                                            // the quadrature
-                                            // points. Here, we don't
-                                            // want to do quadrature,
-                                            // we simply want to
-                                            // specify some points!
-                                            //
-                                            // Nevertheless, the same
-                                            // way is chosen: use a
-                                            // special quadrature
-                                            // rule with points at
-                                            // the vertices, since
-                                            // these are what we are
-                                            // interested in. The
-                                            // appropriate rule is
-                                            // the trapezoidal rule,
-                                            // so that is the reason
-                                            // why we used that one
-                                            // above.
-                                            //
-                                            // Thus: initialize the
-                                            // ``FEValues'' object on
-                                            // this cell,
-           fe_values.reinit (cell);
-                                            // and extract the
-                                            // gradients of the
-                                            // solution vector at the
-                                            // vertices:
-           fe_values.get_function_grads (solution,
-                                         solution_gradients);
-
-                                            // Now we have the
-                                            // gradients at all
-                                            // vertices, so pick out
-                                            // that one which belongs
-                                            // to the evaluation
-                                            // point (note that the
-                                            // order of vertices is
-                                            // not necessarily the
-                                            // same as that of the
-                                            // quadrature points):
-           unsigned int q_point = 0;
-           for (; q_point<solution_gradients.size(); ++q_point)
-             if (fe_values.quadrature_point(q_point) ==
-                 evaluation_point)
-               break;
-
-                                            // Check that the
-                                            // evaluation point was
-                                            // indeed found,
-           Assert (q_point < solution_gradients.size(),
-                   ExcInternalError());
-                                            // and if so take the
-                                            // x-derivative of the
-                                            // gradient there as the
-                                            // value which we are
-                                            // interested in, and
-                                            // increase the counter
-                                            // indicating how often
-                                            // we have added to that
-                                            // variable:
-           point_derivative += solution_gradients[q_point][0];
-           ++evaluation_point_hits;
-
-                                            // Finally break out of
-                                            // the innermost loop
-                                            // iterating over the
-                                            // vertices of the
-                                            // present cell, since if
-                                            // we have found the
-                                            // evaluation point at
-                                            // one vertex it cannot
-                                            // be at a following
-                                            // vertex as well:
-           break;
-         };
-
-                                    // Now we have looped over all
-                                    // cells and vertices, so check
-                                    // whether the point was found:
-    AssertThrow (evaluation_point_hits > 0,
-                ExcEvaluationPointNotFound(evaluation_point));
-
-                                    // We have simply summed up the
-                                    // contributions of all adjacent
-                                    // cells, so we still have to
-                                    // compute the mean value. Once
-                                    // this is done, report the status:
-    point_derivative /= evaluation_point_hits;
-    std::cout << "   Point x-derivative=" << point_derivative
-             << std::endl;
-  }
-
-
-  
-                                  // @sect4{The GridOutput class}
-
-                                  // Since this program has a more
-                                  // difficult structure (it computed
-                                  // a dual solution in addition to a
-                                  // primal one), writing out the
-                                  // solution is no more done by an
-                                  // evaluation object since we want
-                                  // to write both solutions at once
-                                  // into one file, and that requires
-                                  // some more information than
-                                  // available to the evaluation
-                                  // classes.
-                                  //
-                                  // However, we also want to look at
-                                  // the grids generated. This again
-                                  // can be done with one such
-                                  // class. Its structure is analog
-                                  // to the ``SolutionOutput'' class
-                                  // of the previous example program,
-                                  // so we do not discuss it here in
-                                  // more detail. Furthermore,
-                                  // everything that is used here has
-                                  // already been used in previous
-                                  // example programs.
-  template <int dim>
-  class GridOutput : public EvaluationBase<dim>
-  {
-    public:
-      GridOutput (const std::string &output_name_base);
-      
-      virtual void operator () (const DoFHandler<dim> &dof_handler,
-                               const Vector<double>  &solution) const;
-    
-      const std::string output_name_base;
-  };
-
-
-  template <int dim>
-  GridOutput<dim>::
-  GridOutput (const std::string &output_name_base)
-                 :
-                 output_name_base (output_name_base)
-  {}
-  
-
-  template <int dim>
-  void
-  GridOutput<dim>::operator () (const DoFHandler<dim> &dof_handler,
-                               const Vector<double>  &/*solution*/) const
-  {
-#ifdef HAVE_STD_STRINGSTREAM
-    std::ostringstream filename;
-#else
-    std::ostrstream filename;
-#endif
-    filename << "spec2006-447.dealII/"
-            << output_name_base << "-"
-            << this->refinement_cycle
-            << ".eps"
-            << std::ends;
-#ifdef HAVE_STD_STRINGSTREAM
-    std::ofstream out (filename.str().c_str());
-#else
-    std::ofstream out (filename.str());
-#endif
-    
-    GridOut().write_eps (dof_handler.get_tria(), out);
-  }
-}
-
-  
-                                // @sect3{The Laplace solver classes}
-
-                                // Next are the actual solver
-                                // classes. Again, we discuss only
-                                // the differences to the previous
-                                // program.
-namespace LaplaceSolver
-{
-                                  // Before everything else,
-                                  // forward-declare one class that
-                                  // we will have later, since we
-                                  // will want to make it a friend of
-                                  // some of the classes that follow,
-                                  // which requires the class to be
-                                  // known:
-  template <int dim> class WeightedResidual;
-  
-  
-                                  // @sect4{The Laplace solver base class}
-
-                                  // This class is almost unchanged,
-                                  // with the exception that it
-                                  // declares two more functions:
-                                  // ``output_solution'' will be used
-                                  // to generate output files from
-                                  // the actual solutions computed by
-                                  // derived classes, and the
-                                  // ``set_refinement_cycle''
-                                  // function by which the testing
-                                  // framework sets the number of the
-                                  // refinement cycle to a local
-                                  // variable in this class; this
-                                  // number is later used to generate
-                                  // filenames for the solution
-                                  // output.
-  template <int dim>
-  class Base
-  {
-    public:
-      Base (Triangulation<dim> &coarse_grid);
-      virtual ~Base ();
-
-      virtual void solve_problem () = 0;
-      virtual void postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const = 0;
-      virtual void refine_grid () = 0;
-      virtual unsigned int n_dofs () const = 0;
-
-      virtual void set_refinement_cycle (const unsigned int cycle);
-
-      virtual void output_solution () const = 0;
-      
-    
-      const SmartPointer<Triangulation<dim> > triangulation;
-
-      unsigned int refinement_cycle;
-  };
-
-
-  template <int dim>
-  Base<dim>::Base (Triangulation<dim> &coarse_grid)
-                 :
-                 triangulation (&coarse_grid)
-  {}
-
-
-  template <int dim>
-  Base<dim>::~Base () 
-  {}
-
-
-
-  template <int dim>
-  void
-  Base<dim>::set_refinement_cycle (const unsigned int cycle)
-  {
-    refinement_cycle = cycle;
-  }
-  
-
-                                  // @sect4{The Laplace Solver class}
-
-                                  // Likewise, the ``Solver'' class
-                                  // is entirely unchanged and will
-                                  // thus not be discussed.
-  template <int dim>
-  class Solver : public virtual Base<dim>
-  {
-    public:
-      Solver (Triangulation<dim>       &triangulation,
-             const FiniteElement<dim> &fe,
-             const Quadrature<dim>    &quadrature,
-             const Quadrature<dim-1>  &face_quadrature,              
-             const Function<dim>      &boundary_values);
-      virtual
-      ~Solver ();
-
-      virtual
-      void
-      solve_problem ();
-
-      virtual
-      void
-      postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const;
-
-      virtual
-      unsigned int
-      n_dofs () const;
-      
-    
-      const SmartPointer<const FiniteElement<dim> >  fe;
-      const SmartPointer<const Quadrature<dim> >     quadrature;
-      const SmartPointer<const Quadrature<dim-1> >   face_quadrature;      
-      DoFHandler<dim>                                dof_handler;
-      Vector<double>                                 solution;
-      const SmartPointer<const Function<dim> >       boundary_values;
-
-      virtual void assemble_rhs (Vector<double> &rhs) const = 0;
-    
-    
-      struct LinearSystem
-      {
-         LinearSystem (const DoFHandler<dim> &dof_handler);
-
-         void solve (Vector<double> &solution) const;
-       
-         ConstraintMatrix     hanging_node_constraints;
-         SparsityPattern      sparsity_pattern;
-         SparseMatrix<double> matrix;
-         Vector<double>       rhs;
-      };
-
-      void
-      assemble_linear_system (LinearSystem &linear_system);
-
-      void
-      assemble_matrix (LinearSystem                                         &linear_system,
-                      const typename DoFHandler<dim>::active_cell_iterator &begin_cell,
-                      const typename DoFHandler<dim>::active_cell_iterator &end_cell,
-                      Threads::ThreadMutex                                 &mutex) const;
-  };
-
-
-
-  template <int dim>
-  Solver<dim>::Solver (Triangulation<dim>       &triangulation,
-                      const FiniteElement<dim> &fe,
-                      const Quadrature<dim>    &quadrature,
-                      const Quadrature<dim-1>  &face_quadrature,
-                      const Function<dim>      &boundary_values)
-                 :
-                 Base<dim> (triangulation),
-                 fe (&fe),
-                  quadrature (&quadrature),
-                  face_quadrature (&face_quadrature),    
-                 dof_handler (triangulation),
-                 boundary_values (&boundary_values)
-  {}
-
-
-  template <int dim>
-  Solver<dim>::~Solver () 
-  {
-    dof_handler.clear ();
-  }
-
-
-  template <int dim>
-  void
-  Solver<dim>::solve_problem ()
-  {
-    dof_handler.distribute_dofs (*fe);
-    solution.reinit (dof_handler.n_dofs());
-
-    LinearSystem linear_system (dof_handler);
-    assemble_linear_system (linear_system);
-    linear_system.solve (solution);
-  }
-
-
-  template <int dim>
-  void
-  Solver<dim>::
-  postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const
-  {
-    postprocessor (dof_handler, solution);
-  }
-
-
-  template <int dim>
-  unsigned int
-  Solver<dim>::n_dofs () const
-  {
-    return dof_handler.n_dofs();
-  }
-  
-
-  template <int dim>
-  void
-  Solver<dim>::assemble_linear_system (LinearSystem &linear_system)
-  {
-    typedef
-      typename DoFHandler<dim>::active_cell_iterator
-      active_cell_iterator;
-
-    const unsigned int n_threads = multithread_info.n_default_threads;
-    std::vector<std::pair<active_cell_iterator,active_cell_iterator> >
-      thread_ranges 
-      = Threads::split_range<active_cell_iterator> (dof_handler.begin_active (),
-                                                   dof_handler.end (),
-                                                   n_threads);
-
-    Threads::ThreadMutex mutex;
-    Threads::ThreadGroup<> threads;
-    for (unsigned int thread=0; thread<n_threads; ++thread)
-      threads += Threads::spawn (*this, &Solver<dim>::assemble_matrix)
-                 (linear_system,
-                  thread_ranges[thread].first,
-                  thread_ranges[thread].second,
-                  mutex);
-
-    assemble_rhs (linear_system.rhs);
-    linear_system.hanging_node_constraints.condense (linear_system.rhs);
-
-    std::map<unsigned int,double> boundary_value_map;
-    VectorTools::interpolate_boundary_values (dof_handler,
-                                             0,
-                                             *boundary_values,
-                                             boundary_value_map);
-    
-    threads.join_all ();
-    linear_system.hanging_node_constraints.condense (linear_system.matrix);
-
-    MatrixTools::apply_boundary_values (boundary_value_map,
-                                       linear_system.matrix,
-                                       solution,
-                                       linear_system.rhs);
-  }
-
-
-  template <int dim>
-  void
-  Solver<dim>::assemble_matrix (LinearSystem                                         &linear_system,
-                               const typename DoFHandler<dim>::active_cell_iterator &begin_cell,
-                               const typename DoFHandler<dim>::active_cell_iterator &end_cell,
-                               Threads::ThreadMutex                                 &mutex) const
-  {
-    MappingQ<dim> mapping (4);
-    FEValues<dim> fe_values (mapping, *fe, *quadrature, 
-                            UpdateFlags(update_gradients | update_values |
-                                        update_q_points |
-                                        update_JxW_values));
-
-    const unsigned int   dofs_per_cell = fe->dofs_per_cell;
-    const unsigned int   n_q_points    = quadrature->n_quadrature_points;
-
-    FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
-
-    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
-
-    std::vector<double> laplace_coefficients (fe_values.n_quadrature_points);
-    std::vector<double> mass_coefficients (fe_values.n_quadrature_points);
-    
-    
-    for (typename DoFHandler<dim>::active_cell_iterator cell=begin_cell;
-        cell!=end_cell; ++cell)
-      {
-       cell_matrix = 0;
-
-       fe_values.reinit (cell);
-
-       LaplaceCoefficient<dim>().value_list (fe_values.get_quadrature_points(),
-                                             laplace_coefficients);
-       MassCoefficient<dim>().value_list (fe_values.get_quadrature_points(),
-                                          mass_coefficients);
-       
-       
-       for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
-         for (unsigned int i=0; i<dofs_per_cell; ++i)
-           for (unsigned int j=0; j<dofs_per_cell; ++j)
-             cell_matrix(i,j) += ((fe_values.shape_grad(i,q_point) *
-                                   fe_values.shape_grad(j,q_point) *
-                                   laplace_coefficients[q_point]
-                                   +
-                                   fe_values.shape_value(i,q_point) *
-                                   fe_values.shape_value(j,q_point) *
-                                   mass_coefficients[q_point]
-                                  ) *
-                                  fe_values.JxW(q_point));
-
-
-       cell->get_dof_indices (local_dof_indices);
-       Threads::ThreadMutex::ScopedLock lock (mutex);
-       for (unsigned int i=0; i<dofs_per_cell; ++i)
-         for (unsigned int j=0; j<dofs_per_cell; ++j)
-           linear_system.matrix.add (local_dof_indices[i],
-                                     local_dof_indices[j],
-                                     cell_matrix(i,j));
-      };
-  }
-
-
-  template <int dim>
-  Solver<dim>::LinearSystem::
-  LinearSystem (const DoFHandler<dim> &dof_handler)
-  {
-    hanging_node_constraints.clear ();
-
-    void (*mhnc_p) (const DoFHandler<dim> &,
-                   ConstraintMatrix      &)
-      = &DoFTools::make_hanging_node_constraints;
-    
-    Threads::Thread<>
-      mhnc_thread = Threads::spawn (mhnc_p)(dof_handler, hanging_node_constraints);
-
-                                    // make sparsity pattern. since
-                                    // in 3d the usual way just blows
-                                    // the roof w.r.t. memory
-                                    // consumption, use the detour
-                                    // via a compressed sparsity
-                                    // pattern that we later copy
-                                    // over
-    CompressedSparsityPattern csp(dof_handler.n_dofs(),
-                                 dof_handler.n_dofs());
-    DoFTools::make_sparsity_pattern (dof_handler, csp);
-
-    mhnc_thread.join ();
-    hanging_node_constraints.close ();
-    hanging_node_constraints.condense (csp);
-
-    sparsity_pattern.copy_from(csp);
-    matrix.reinit (sparsity_pattern);
-    rhs.reinit (dof_handler.n_dofs());
-  }
-
-
-
-  template <int dim>
-  void
-  Solver<dim>::LinearSystem::solve (Vector<double> &solution) const
-  {
-    SolverControl           solver_control (solution.size(), 1e-6);
-    PrimitiveVectorMemory<> vector_memory;
-    SolverCG<>              cg (solver_control, vector_memory);
-
-    PreconditionJacobi<> preconditioner;
-    preconditioner.initialize(matrix);
-
-    cg.solve (matrix, solution, rhs, preconditioner);
-
-    hanging_node_constraints.distribute (solution);
-  }
-
-
-
-
-                                  // @sect4{The PrimalSolver class}
-
-                                  // The ``PrimalSolver'' class is
-                                  // also mostly unchanged except for
-                                  // overloading the functions
-                                  // ``solve_problem'', ``n_dofs'',
-                                  // and ``postprocess'' of the base
-                                  // class, and implementing the
-                                  // ``output_solution''
-                                  // function. These overloaded
-                                  // functions do nothing particular
-                                  // besides calling the functions of
-                                  // the base class -- that seems
-                                  // superfluous, but works around a
-                                  // bug in a popular compiler which
-                                  // requires us to write such
-                                  // functions for the following
-                                  // scenario: Besides the
-                                  // ``PrimalSolver'' class, we will
-                                  // have a ``DualSolver'', both
-                                  // derived from ``Solver''. We will
-                                  // then have a final classes which
-                                  // derived from these two, which
-                                  // will then have two instances of
-                                  // the ``Solver'' class as its base
-                                  // classes. If we want, for
-                                  // example, the number of degrees
-                                  // of freedom of the primal solver,
-                                  // we would have to indicate this
-                                  // like so:
-                                  // ``PrimalSolver<dim>::n_dofs()''.
-                                  // However, the compiler does not
-                                  // accept this since the ``n_dofs''
-                                  // function is actually from a base
-                                  // class of the ``PrimalSolver''
-                                  // class, so we have to inject the
-                                  // name from the base to the
-                                  // derived class using these
-                                  // additional functions.
-                                  //
-                                  // Regarding the implementation of
-                                  // the ``output_solution''
-                                  // function, we keep the
-                                  // ``GlobalRefinement'' and
-                                  // ``RefinementKelly'' classes in
-                                  // this program, and they can then
-                                  // rely on the default
-                                  // implementation of this function
-                                  // which simply outputs the primal
-                                  // solution. The class implementing
-                                  // dual weighted error estimators
-                                  // will overload this function
-                                  // itself, to also output the dual
-                                  // solution.
-                                  //
-                                  // Except for this, the class is
-                                  // unchanged with respect to the
-                                  // previous example.
-  template <int dim>
-  class PrimalSolver : public Solver<dim>
-  {
-    public:
-      PrimalSolver (Triangulation<dim>       &triangulation,
-                   const FiniteElement<dim> &fe,
-                   const Quadrature<dim>    &quadrature,
-                   const Quadrature<dim-1>  &face_quadrature,
-                   const Function<dim>      &rhs_function,
-                   const Function<dim>      &boundary_values);
-
-      virtual
-      void solve_problem ();
-      
-      virtual
-      unsigned int n_dofs () const;
-      
-      virtual
-      void postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const;
-
-      virtual
-      void output_solution () const;
-      
-    
-      const SmartPointer<const Function<dim> > rhs_function;
-      virtual void assemble_rhs (Vector<double> &rhs) const;
-
-                                      // Now, in order to work around
-                                      // some problems in one of the
-                                      // compilers this library can
-                                      // be compiled with, we will
-                                      // have to use some
-                                      // workarounds. This will
-                                      // require that we declare a
-                                      // class that is actually
-                                      // derived from the present
-                                      // one, as a friend (strange as
-                                      // that seems). The full
-                                      // rationale will be explained
-                                      // below.
-      friend class WeightedResidual<dim>;
-  };
-
-
-  template <int dim>
-  PrimalSolver<dim>::
-  PrimalSolver (Triangulation<dim>       &triangulation,
-               const FiniteElement<dim> &fe,
-               const Quadrature<dim>    &quadrature,
-               const Quadrature<dim-1>  &face_quadrature,
-               const Function<dim>      &rhs_function,
-               const Function<dim>      &boundary_values)
-                 :
-                 Base<dim> (triangulation),
-                 Solver<dim> (triangulation, fe,
-                              quadrature, face_quadrature,
-                              boundary_values),
-                  rhs_function (&rhs_function)
-  {}
-
-
-  template <int dim>
-  void
-  PrimalSolver<dim>::solve_problem ()
-  {
-    Solver<dim>::solve_problem ();
-  }
-
-
-
-  template <int dim>
-  unsigned int
-  PrimalSolver<dim>::n_dofs() const
-  {
-    return Solver<dim>::n_dofs();
-  }
-
-
-  template <int dim>
-  void
-  PrimalSolver<dim>::
-  postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const
-  {
-    Solver<dim>::postprocess(postprocessor);
-  }
-
-
-  template <int dim>
-  void
-  PrimalSolver<dim>::output_solution () const
-  {
-    abort ();
-  }
-  
-
-
-  template <int dim>
-  void
-  PrimalSolver<dim>::
-  assemble_rhs (Vector<double> &rhs) const 
-  {
-    MappingQ<dim> mapping (4);
-    FEValues<dim> fe_values (mapping, *this->fe, *this->quadrature, 
-                            UpdateFlags(update_values    |
-                                        update_q_points  |
-                                        update_JxW_values));
-
-    const unsigned int   dofs_per_cell = this->fe->dofs_per_cell;
-    const unsigned int   n_q_points    = this->quadrature->n_quadrature_points;
-
-    Vector<double>       cell_rhs (dofs_per_cell);
-    std::vector<double>  rhs_values (n_q_points);
-    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
-
-    typename DoFHandler<dim>::active_cell_iterator
-      cell = this->dof_handler.begin_active(),
-      endc = this->dof_handler.end();
-    for (; cell!=endc; ++cell)
-      {
-       cell_rhs = 0;
-
-       fe_values.reinit (cell);
-
-       rhs_function->value_list (fe_values.get_quadrature_points(),
-                                 rhs_values);
-      
-       for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
-         for (unsigned int i=0; i<dofs_per_cell; ++i)
-           cell_rhs(i) += (fe_values.shape_value(i,q_point) *
-                           rhs_values[q_point] *
-                           fe_values.JxW(q_point));
-
-       cell->get_dof_indices (local_dof_indices);
-       for (unsigned int i=0; i<dofs_per_cell; ++i)
-         rhs(local_dof_indices[i]) += cell_rhs(i);
-      }
-  }
-
-
-                                  // @sect4{The RefinementGlobal and RefinementKelly classes}
-
-                                  // For the following two classes,
-                                  // the same applies as for most of
-                                  // the above: the class is taken
-                                  // from the previous example as-is:
-  template <int dim>
-  class RefinementGlobal : public PrimalSolver<dim>
-  {
-    public:
-      RefinementGlobal (Triangulation<dim>       &coarse_grid,
-                       const FiniteElement<dim> &fe,
-                       const Quadrature<dim>    &quadrature,
-                       const Quadrature<dim-1>  &face_quadrature,
-                       const Function<dim>      &rhs_function,
-                       const Function<dim>      &boundary_values);
-
-      virtual void refine_grid ();
-  };
-
-
-
-  template <int dim>
-  RefinementGlobal<dim>::
-  RefinementGlobal (Triangulation<dim>       &coarse_grid,
-                   const FiniteElement<dim> &fe,
-                   const Quadrature<dim>    &quadrature,
-                   const Quadrature<dim-1>  &face_quadrature,
-                   const Function<dim>      &rhs_function,
-                   const Function<dim>      &boundary_values)
-                 :
-                 Base<dim> (coarse_grid),
-                  PrimalSolver<dim> (coarse_grid, fe, quadrature,
-                                    face_quadrature, rhs_function,
-                                    boundary_values)
-  {}
-
-
-
-  template <int dim>
-  void
-  RefinementGlobal<dim>::refine_grid ()
-  {
-    this->triangulation->refine_global (1);
-  }
-
-
-
-  template <int dim>
-  class RefinementKelly : public PrimalSolver<dim>
-  {
-    public:
-      RefinementKelly (Triangulation<dim>       &coarse_grid,
-                      const FiniteElement<dim> &fe,
-                      const Quadrature<dim>    &quadrature,
-                      const Quadrature<dim-1>  &face_quadrature,
-                      const Function<dim>      &rhs_function,
-                      const Function<dim>      &boundary_values);
-
-      virtual void refine_grid ();
-  };
-
-
-
-  template <int dim>
-  RefinementKelly<dim>::
-  RefinementKelly (Triangulation<dim>       &coarse_grid,
-                  const FiniteElement<dim> &fe,
-                  const Quadrature<dim>    &quadrature,
-                  const Quadrature<dim-1>  &face_quadrature,
-                  const Function<dim>      &rhs_function,
-                  const Function<dim>      &boundary_values)
-                 :
-                 Base<dim> (coarse_grid),
-                  PrimalSolver<dim> (coarse_grid, fe, quadrature,
-                                    face_quadrature,
-                                    rhs_function, boundary_values)
-  {}
-
-
-
-  template <int dim>
-  void
-  RefinementKelly<dim>::refine_grid ()
-  {
-    Vector<float> estimated_error_per_cell (this->triangulation->n_active_cells());
-    KellyErrorEstimator<dim>::estimate (this->dof_handler,
-                                       QGauss<dim-1>(3),
-                                       typename FunctionMap<dim>::type(),
-                                       this->solution,
-                                       estimated_error_per_cell);
-    GridRefinement::refine_and_coarsen_fixed_number (*this->triangulation,
-                                                    estimated_error_per_cell,
-                                                    0.2, 0.02);
-    this->triangulation->execute_coarsening_and_refinement ();
-  }
-
-
-
-                                  // @sect4{The RefinementWeightedKelly class}
-
-                                  // This class is a variant of the
-                                  // previous one, in that it allows
-                                  // to weight the refinement
-                                  // indicators we get from the
-                                  // library's Kelly indicator by
-                                  // some function. We include this
-                                  // class since the goal of this
-                                  // example program is to
-                                  // demonstrate automatic refinement
-                                  // criteria even for complex output
-                                  // quantities such as point values
-                                  // or stresses. If we did not solve
-                                  // a dual problem and compute the
-                                  // weights thereof, we would
-                                  // probably be tempted to give a
-                                  // hand-crafted weighting to the
-                                  // indicators to account for the
-                                  // fact that we are going to
-                                  // evaluate these quantities. This
-                                  // class accepts such a weighting
-                                  // function as argument to its
-                                  // constructor:
-  template <int dim>
-  class RefinementWeightedKelly : public PrimalSolver<dim>
-  {
-    public:
-      RefinementWeightedKelly (Triangulation<dim>       &coarse_grid,
-                              const FiniteElement<dim> &fe,
-                              const Quadrature<dim>    &quadrature,
-                              const Quadrature<dim-1>  &face_quadrature,
-                              const Function<dim>      &rhs_function,
-                              const Function<dim>      &boundary_values,
-                              const Function<dim>      &weighting_function);
-
-      virtual void refine_grid ();
-
-    
-      const SmartPointer<const Function<dim> > weighting_function;
-  };
-
-
-
-  template <int dim>
-  RefinementWeightedKelly<dim>::
-  RefinementWeightedKelly (Triangulation<dim>       &coarse_grid,
-                          const FiniteElement<dim> &fe,
-                          const Quadrature<dim>    &quadrature,
-                          const Quadrature<dim-1>  &face_quadrature,
-                          const Function<dim>      &rhs_function,
-                          const Function<dim>      &boundary_values,
-                          const Function<dim>      &weighting_function)
-                 :
-                 Base<dim> (coarse_grid),
-                  PrimalSolver<dim> (coarse_grid, fe, quadrature,
-                                    face_quadrature,
-                                    rhs_function, boundary_values),
-                  weighting_function (&weighting_function)
-  {}
-
-
-
-                                  // Now, here comes the main
-                                  // function, including the
-                                  // weighting:
-  template <int dim>
-  void
-  RefinementWeightedKelly<dim>::refine_grid ()
-  {
-                                    // First compute some residual
-                                    // based error indicators for all
-                                    // cells by a method already
-                                    // implemented in the
-                                    // library. What exactly is
-                                    // computed can be read in the
-                                    // documentation of that class.
-    Vector<float> estimated_error (this->triangulation->n_active_cells());
-    KellyErrorEstimator<dim>::estimate (this->dof_handler,
-                                       *this->face_quadrature,
-                                       typename FunctionMap<dim>::type(),
-                                       this->solution,
-                                       estimated_error);
-
-                                    // Now we are going to weight
-                                    // these indicators by the value
-                                    // of the function given to the
-                                    // constructor:
-    typename DoFHandler<dim>::active_cell_iterator
-      cell = this->dof_handler.begin_active(),
-      endc = this->dof_handler.end();
-    for (unsigned int cell_index=0; cell!=endc; ++cell, ++cell_index)
-      estimated_error(cell_index)
-       *= weighting_function->value (cell->center());
-    
-    GridRefinement::refine_and_coarsen_fixed_number (*this->triangulation,
-                                                    estimated_error,
-                                                    0.2, 0.02);
-    this->triangulation->execute_coarsening_and_refinement ();
-  }
-
-}
-
-
-                                // @sect3{Equation data}
-                                //
-                                // In this example program, we work
-                                // with the same data sets as in the
-                                // previous one, but as it may so
-                                // happen that someone wants to run
-                                // the program with different
-                                // boundary values and right hand side
-                                // functions, or on a different grid,
-                                // we show a simple technique to do
-                                // exactly that. For more clarity, we
-                                // furthermore pack everything that
-                                // has to do with equation data into
-                                // a namespace of its own.
-                                //
-                                // The underlying assumption is that
-                                // this is a research program, and
-                                // that there we often have a number
-                                // of test cases that consist of a
-                                // domain, a right hand side,
-                                // boundary values, possibly a
-                                // specified coefficient, and a
-                                // number of other parameters. They
-                                // often vary all at the same time
-                                // when shifting from one example to
-                                // another. To make handling such
-                                // sets of problem description
-                                // parameters simple is the goal of
-                                // the following.
-                                //
-                                // Basically, the idea is this: let
-                                // us have a structure for each set
-                                // of data, in which we pack
-                                // everything that describes a test
-                                // case: here, these are two
-                                // subclasses, one called
-                                // ``BoundaryValues'' for the
-                                // boundary values of the exact
-                                // solution, and one called
-                                // ``RightHandSide'', and then a way
-                                // to generate the coarse grid. Since
-                                // the solution of the previous
-                                // example program looked like curved
-                                // ridges, we use this name here for
-                                // the enclosing class. Note that the
-                                // names of the two inner classes
-                                // have to be the same for all
-                                // enclosing test case classes, and
-                                // also that we have attached the
-                                // dimension template argument to the
-                                // enclosing class rather than to the
-                                // inner ones, to make further
-                                // processing simpler.  (From a
-                                // language viewpoint, a namespace
-                                // would be better to encapsulate
-                                // these inner classes, rather than a
-                                // structure. However, namespaces
-                                // cannot be given as template
-                                // arguments, so we use a structure
-                                // to allow a second object to select
-                                // from within its given
-                                // argument. The enclosing structure,
-                                // of course, has no member variables
-                                // apart from the classes it
-                                // declares, and a static function to
-                                // generate the coarse mesh; it will
-                                // in general never be instantiated.)
-                                //
-                                // The idea is then the following
-                                // (this is the right time to also
-                                // take a brief look at the code
-                                // below): we can generate objects
-                                // for boundary values and
-                                // right hand side by simply giving
-                                // the name of the outer class as a
-                                // template argument to a class which
-                                // we call here ``Data::SetUp'', and
-                                // it then creates objects for the
-                                // inner classes. In this case, to
-                                // get all that characterizes the
-                                // curved ridge solution, we would
-                                // simply generate an instance of
-                                // ``Data::SetUp<Data::CurvedRidge>'',
-                                // and everything we need to know
-                                // about the solution would be static
-                                // member variables and functions of
-                                // that object.
-                                //
-                                // This approach might seem like
-                                // overkill in this case, but will
-                                // become very handy once a certain
-                                // set up is not only characterized
-                                // by Dirichlet boundary values and a
-                                // right hand side function, but in
-                                // addition by material properties,
-                                // Neumann values, different boundary
-                                // descriptors, etc. In that case,
-                                // the ``SetUp'' class might consist
-                                // of a dozen or more objects, and
-                                // each descriptor class (like the
-                                // ``CurvedRidges'' class below)
-                                // would have to provide them. Then,
-                                // you will be happy to be able to
-                                // change from one set of data to
-                                // another by only changing the
-                                // template argument to the ``SetUp''
-                                // class at one place, rather than at
-                                // many.
-                                //
-                                // With this framework for different
-                                // test cases, we are almost
-                                // finished, but one thing remains:
-                                // by now we can select statically,
-                                // by changing one template argument,
-                                // which data set to choose. In order
-                                // to be able to do that dynamically,
-                                // i.e. at run time, we need a base
-                                // class. This we provide in the
-                                // obvious way, see below, with
-                                // virtual abstract functions. It
-                                // forces us to introduce a second
-                                // template parameter ``dim'' which
-                                // we need for the base class (which
-                                // could be avoided using some
-                                // template magic, but we omit that),
-                                // but that's all.
-                                //
-                                // Adding new testcases is now
-                                // simple, you don't have to touch
-                                // the framework classes, only a
-                                // structure like the
-                                // ``CurvedRidges'' one is needed.
-namespace Data
-{
-                                  // @sect4{The SetUpBase and SetUp classes}
-  
-                                  // Based on the above description,
-                                  // the ``SetUpBase'' class then
-                                  // looks as follows. To allow using
-                                  // the ``SmartPointer'' class with
-                                  // this class, we derived from the
-                                  // ``Subscriptor'' class.
-  template <int dim>
-  struct SetUpBase : public Subscriptor
-  {
-      virtual
-      const Function<dim> &  get_boundary_values () const = 0;
-
-      virtual
-      const Function<dim> &  get_right_hand_side () const = 0;
-
-      virtual
-      void create_coarse_grid (Triangulation<dim> &coarse_grid) const = 0;
-  };
-
-
-                                  // And now for the derived class
-                                  // that takes the template argument
-                                  // as explained above. For some
-                                  // reason, C++ requires us to
-                                  // define a constructor (which
-                                  // maybe empty), as otherwise a
-                                  // warning is generated that some
-                                  // data is not initialized.
-                                  //
-                                  // Here we pack the data elements
-                                  // into private variables, and
-                                  // allow access to them through the
-                                  // methods of the base class.
-  template <class Traits, int dim>
-  struct SetUp : public SetUpBase<dim>
-  {
-      SetUp () {};
-
-      virtual
-      const Function<dim> &  get_boundary_values () const;
-
-      virtual
-      const Function<dim> &  get_right_hand_side () const;
-      
-
-      virtual
-      void create_coarse_grid (Triangulation<dim> &coarse_grid) const;
-
-    
-      static const typename Traits::BoundaryValues boundary_values;
-      static const typename Traits::RightHandSide  right_hand_side;
-  };
-
-                                  // We have to provide definitions
-                                  // for the static member variables
-                                  // of the above class:
-  template <class Traits, int dim>
-  const typename Traits::BoundaryValues  SetUp<Traits,dim>::boundary_values;
-  template <class Traits, int dim>
-  const typename Traits::RightHandSide   SetUp<Traits,dim>::right_hand_side;
-
-                                  // And definitions of the member
-                                  // functions:
-  template <class Traits, int dim>
-  const Function<dim> &
-  SetUp<Traits,dim>::get_boundary_values () const 
-  {
-    return boundary_values;
-  }
-
-
-  template <class Traits, int dim>
-  const Function<dim> &
-  SetUp<Traits,dim>::get_right_hand_side () const 
-  {
-    return right_hand_side;
-  }
-
-
-  template <class Traits, int dim>
-  void
-  SetUp<Traits,dim>::
-  create_coarse_grid (Triangulation<dim> &coarse_grid) const 
-  {
-    Traits::create_coarse_grid (coarse_grid);
-  }
-  
-
-                                  // @sect4{The CurvedRidges class}
-
-                                  // The class that is used to
-                                  // describe the boundary values and
-                                  // right hand side of the ``curved
-                                  // ridge'' problem already used in
-                                  // the step-13 example program is
-                                  // then like so:
-  template <int dim>
-  struct CurvedRidges
-  {
-      class BoundaryValues : public Function<dim>
-      {
-       public:
-         BoundaryValues () : Function<dim> () {};
-         
-         virtual double value (const Point<dim>   &p,
-                               const unsigned int  component) const;
-      };
-
-
-      class RightHandSide : public Function<dim>
-      {
-       public:
-         RightHandSide () : Function<dim> () {};
-         
-         virtual double value (const Point<dim>   &p,
-                               const unsigned int  component) const;
-      };
-
-      static
-      void
-      create_coarse_grid (Triangulation<dim> &coarse_grid);
-  };
-  
-    
-  template <int dim>
-  double
-  CurvedRidges<dim>::BoundaryValues::
-  value (const Point<dim>   &p,
-        const unsigned int  /*component*/) const
-  {
-    double q = p(0);
-    for (unsigned int i=1; i<dim; ++i)
-      q += std::sin(10*p(i)+5*p(0)*p(0));
-    const double exponential = std::exp(q);
-    return exponential;
-  }
-
-
-
-  template <int dim>
-  double
-  CurvedRidges<dim>::RightHandSide::value (const Point<dim>   &p,
-                                          const unsigned int  /*component*/) const
-  {
-    double q = p(0);
-    for (unsigned int i=1; i<dim; ++i)
-      q += std::sin(10*p(i)+5*p(0)*p(0));
-    const double u = std::exp(q);
-    double t1 = 1,
-          t2 = 0,
-          t3 = 0;
-    for (unsigned int i=1; i<dim; ++i)
-      {
-       t1 += std::cos(10*p(i)+5*p(0)*p(0)) * 10 * p(0);
-       t2 += 10*std::cos(10*p(i)+5*p(0)*p(0)) -
-             100*std::sin(10*p(i)+5*p(0)*p(0)) * p(0)*p(0);
-       t3 += 100*std::cos(10*p(i)+5*p(0)*p(0))*std::cos(10*p(i)+5*p(0)*p(0)) -
-             100*std::sin(10*p(i)+5*p(0)*p(0));
-      };
-    t1 = t1*t1;
-    
-    return -u*(t1+t2+t3);
-  }
-
-
-  template <int dim>
-  void
-  CurvedRidges<dim>::
-  create_coarse_grid (Triangulation<dim> &coarse_grid)
-  {
-    GridGenerator::hyper_cube (coarse_grid, -1, 1);
-    coarse_grid.refine_global (2);
-  }
-  
-
-                                  // @sect4{The Exercise_2_3 class}
-  
-                                  // This example program was written
-                                  // while giving practical courses
-                                  // for a lecture on adaptive finite
-                                  // element methods and duality
-                                  // based error estimates. For these
-                                  // courses, we had one exercise,
-                                  // which required to solve the
-                                  // Laplace equation with constant
-                                  // right hand side on a square
-                                  // domain with a square hole in the
-                                  // center, and zero boundary
-                                  // values. Since the implementation
-                                  // of the properties of this
-                                  // problem is so particularly
-                                  // simple here, lets do it. As the
-                                  // number of the exercise was 2.3,
-                                  // we take the liberty to retain
-                                  // this name for the class as well.
-  template <int dim>
-  struct Exercise_2_3
-  {
-                                      // We need a class to denote
-                                      // the boundary values of the
-                                      // problem. In this case, this
-                                      // is simple: it's the zero
-                                      // function, so don't even
-                                      // declare a class, just a
-                                      // typedef:
-      typedef ZeroFunction<dim> BoundaryValues;
-
-                                      // Second, a class that denotes
-                                      // the right hand side. Since
-                                      // they are constant, just
-                                      // subclass the corresponding
-                                      // class of the library and be
-                                      // done:
-      class RightHandSide : public ConstantFunction<dim>
-      {
-       public:
-         RightHandSide () : ConstantFunction<dim> (1.) {};
-      };
-      
-                                      // Finally a function to
-                                      // generate the coarse
-                                      // grid. This is somewhat more
-                                      // complicated here, see
-                                      // immediately below.
-      static
-      void
-      create_coarse_grid (Triangulation<dim> &coarse_grid);
-  };
-
-
-                                  // As stated above, the grid for
-                                  // this example is the square
-                                  // [-1,1]^2 with the square
-                                  // [-1/2,1/2]^2 as hole in it. We
-                                  // create the coarse grid as 4
-                                  // times 4 cells with the middle
-                                  // four ones missing.
-                                  //
-                                  // Of course, the example has an
-                                  // extension to 3d, but since this
-                                  // function cannot be written in a
-                                  // dimension independent way we
-                                  // choose not to implement this
-                                  // here, but rather only specialize
-                                  // the template for dim=2. If you
-                                  // compile the program for 3d,
-                                  // you'll get a message from the
-                                  // linker that this function is not
-                                  // implemented for 3d, and needs to
-                                  // be provided.
-                                  //
-                                  // For the creation of this
-                                  // geometry, the library has no
-                                  // predefined method. In this case,
-                                  // the geometry is still simple
-                                  // enough to do the creation by
-                                  // hand, rather than using a mesh
-                                  // generator.
-/*  
-  template <>
-  void
-  Exercise_2_3<2>::
-  create_coarse_grid (Triangulation<2> &coarse_grid)
-  {
-                                    // First define the space
-                                    // dimension, to allow those
-                                    // parts of the function that are
-                                    // actually dimension independent
-                                    // to use this variable. That
-                                    // makes it simpler if you later
-                                    // takes this as a starting point
-                                    // to implement the 3d version.
-    const unsigned int dim = 2;
-
-                                    // Then have a list of
-                                    // vertices. Here, they are 24 (5
-                                    // times 5, with the middle one
-                                    // omitted). It is probably best
-                                    // to draw a sketch here. Note
-                                    // that we leave the number of
-                                    // vertices open at first, but
-                                    // then let the compiler compute
-                                    // this number afterwards. This
-                                    // reduces the possibility of
-                                    // having the dimension to large
-                                    // and leaving the last ones
-                                    // uninitialized.
-    static const Point<2> vertices_1[]
-      = {  Point<2> (-1.,   -1.),
-            Point<2> (-1./2, -1.),
-            Point<2> (0.,    -1.),
-            Point<2> (+1./2, -1.),
-            Point<2> (+1,    -1.),
-            
-            Point<2> (-1.,   -1./2.),
-            Point<2> (-1./2, -1./2.),
-            Point<2> (0.,    -1./2.),
-            Point<2> (+1./2, -1./2.),
-            Point<2> (+1,    -1./2.),
-            
-            Point<2> (-1.,   0.),
-            Point<2> (-1./2, 0.),
-            Point<2> (+1./2, 0.),
-            Point<2> (+1,    0.),
-            
-            Point<2> (-1.,   1./2.),
-            Point<2> (-1./2, 1./2.),
-            Point<2> (0.,    1./2.),
-            Point<2> (+1./2, 1./2.),
-            Point<2> (+1,    1./2.),
-            
-            Point<2> (-1.,   1.),
-            Point<2> (-1./2, 1.),
-            Point<2> (0.,    1.),                        
-            Point<2> (+1./2, 1.),
-            Point<2> (+1,    1.)    };
-    const unsigned int
-      n_vertices = sizeof(vertices_1) / sizeof(vertices_1[0]);
-
-                                    // From this static list of
-                                    // vertices, we generate an STL
-                                    // vector of the vertices, as
-                                    // this is the data type the
-                                    // library wants to see.
-    const std::vector<Point<dim> > vertices (&vertices_1[0],
-                                            &vertices_1[n_vertices]);
-
-                                    // Next, we have to define the
-                                    // cells and the vertices they
-                                    // contain. Here, we have 8
-                                    // vertices, but leave the number
-                                    // open and let it be computed
-                                    // afterwards:
-    static const int cell_vertices[][GeometryInfo<dim>::vertices_per_cell]
-      = {{0, 1, 6,5},
-        {1, 2, 7, 6},
-        {2, 3, 8, 7},
-        {3, 4, 9, 8},
-        {5, 6, 11, 10},
-        {8, 9, 13, 12},
-        {10, 11, 15, 14},
-        {12, 13, 18, 17},
-        {14, 15, 20, 19},
-        {15, 16, 21, 20},
-        {16, 17, 22, 21},
-        {17, 18, 23, 22}};
-    const unsigned int
-      n_cells = sizeof(cell_vertices) / sizeof(cell_vertices[0]);
-
-                                    // Again, we generate a C++
-                                    // vector type from this, but
-                                    // this time by looping over the
-                                    // cells (yes, this is
-                                    // boring). Additionally, we set
-                                    // the material indicator to zero
-                                    // for all the cells:
-    std::vector<CellData<dim> > cells (n_cells, CellData<dim>());
-    for (unsigned int i=0; i<n_cells; ++i) 
-      {
-       for (unsigned int j=0;
-            j<GeometryInfo<dim>::vertices_per_cell;
-            ++j)
-         cells[i].vertices[j] = cell_vertices[i][j];
-       cells[i].material_id = 0;
-      };
-
-                                    // Finally pass all this
-                                    // information to the library to
-                                    // generate a triangulation. The
-                                    // last parameter may be used to
-                                    // pass information about
-                                    // non-zero boundary indicators
-                                    // at certain faces of the
-                                    // triangulation to the library,
-                                    // but we don't want that here,
-                                    // so we give an empty object:
-    coarse_grid.create_triangulation (vertices,
-                                     cells,
-                                     SubCellData());
-    
-                                    // And since we want that the
-                                    // evaluation point (3/4,3/4) in
-                                    // this example is a grid point,
-                                    // we refine once globally:
-    coarse_grid.refine_global (1);
-  }
-*/  
-
-
-  template <>
-  void
-  Exercise_2_3<3>::
-  create_coarse_grid (Triangulation<3> &coarse_grid)
-  {
-    GridGenerator::hyper_ball (coarse_grid);
-    static HyperBallBoundary<3> boundary;
-    coarse_grid.set_boundary (0, boundary);
-    coarse_grid.refine_global (1);
-  }
-  
-}
-
-                                // @sect4{Discussion}
-                                //
-                                // As you have now read through this
-                                // framework, you may be wondering
-                                // why we have not chosen to
-                                // implement the classes implementing
-                                // a certain setup (like the
-                                // ``CurvedRidges'' class) directly
-                                // as classes derived from
-                                // ``Data::SetUpBase''. Indeed, we
-                                // could have done very well so. The
-                                // only reason is that then we would
-                                // have to have member variables for
-                                // the solution and right hand side
-                                // classes in the ``CurvedRidges''
-                                // class, as well as member functions
-                                // overloading the abstract functions
-                                // of the base class giving access to
-                                // these member variables. The
-                                // ``SetUp'' class has the sole
-                                // reason to relieve us from the need
-                                // to reiterate these member
-                                // variables and functions that would
-                                // be necessary in all such
-                                // classes. In some way, the template
-                                // mechanism here only provides a way
-                                // to have default implementations
-                                // for a number of functions that
-                                // depend on external quantities and
-                                // can thus not be provided using
-                                // normal virtual functions, at least
-                                // not without the help of templates.
-                                //
-                                // However, there might be good
-                                // reasons to actually implement
-                                // classes derived from
-                                // ``Data::SetUpBase'', for example
-                                // if the solution or right hand side
-                                // classes require constructors that
-                                // take arguments, which the
-                                // ``Data::SetUpBase'' class cannot
-                                // provide. In that case, subclassing
-                                // is a worthwhile strategy. Other
-                                // possibilities for special cases
-                                // are to derive from
-                                // ``Data::SetUp<SomeSetUp>'' where
-                                // ``SomeSetUp'' denotes a class, or
-                                // even to explicitly specialize
-                                // ``Data::SetUp<SomeSetUp>''. The
-                                // latter allows to transparently use
-                                // the way the ``SetUp'' class is
-                                // used for other set-ups, but with
-                                // special actions taken for special
-                                // arguments.
-                                //
-                                // A final observation favoring the
-                                // approach taken here is the
-                                // following: we have found numerous
-                                // times that when starting a
-                                // project, the number of parameters
-                                // (usually boundary values, right
-                                // hand side, coarse grid, just as
-                                // here) was small, and the number of
-                                // test cases was small as well. One
-                                // then starts out by handcoding them
-                                // into a number of ``switch''
-                                // statements. Over time, projects
-                                // grow, and so does the number of
-                                // test cases. The number of
-                                // ``switch'' statements grows with
-                                // that, and their length as well,
-                                // and one starts to find ways to
-                                // consider impossible examples where
-                                // domains, boundary values, and
-                                // right hand sides do not fit
-                                // together any more, and starts
-                                // loosing the overview over the
-                                // whole structure. Encapsulating
-                                // everything belonging to a certain
-                                // test case into a structure of its
-                                // own has proven worthwhile for
-                                // this, as it keeps everything that
-                                // belongs to one test case in one
-                                // place. Furthermore, it allows to
-                                // put these things all in one or
-                                // more files that are only devoted
-                                // to test cases and their data,
-                                // without having to bring their
-                                // actual implementation into contact
-                                // with the rest of the program.
-
-
-                                // @sect3{Dual functionals}
-
-                                // As with the other components of
-                                // the program, we put everything we
-                                // need to describe dual functionals
-                                // into a namespace of its own, and
-                                // define an abstract base class that
-                                // provides the interface the class
-                                // solving the dual problem needs for
-                                // its work.
-                                //
-                                // We will then implement two such
-                                // classes, for the evaluation of a
-                                // point value and of the derivative
-                                // of the solution at that point. For
-                                // these functionals we already have
-                                // the corresponding evaluation
-                                // objects, so they are comlementary.
-namespace DualFunctional
-{
-                                  // @sect4{The DualFunctionalBase class}
-  
-                                  // First start with the base class
-                                  // for dual functionals. Since for
-                                  // linear problems the
-                                  // characteristics of the dual
-                                  // problem play a role only in the
-                                  // right hand side, we only need to
-                                  // provide for a function that
-                                  // assembles the right hand side
-                                  // for a given discretization:
-  template <int dim>
-  class DualFunctionalBase : public Subscriptor
-  {
-    public:
-      virtual
-      void
-      assemble_rhs (const DoFHandler<dim> &dof_handler,
-                   Vector<double>        &rhs) const = 0;
-  };
-
-
-                                  // @sect4{The PointValueEvaluation class}
-  
-                                  // As a first application, we
-                                  // consider the functional
-                                  // corresponding to the evaluation
-                                  // of the solution's value at a
-                                  // given point which again we
-                                  // assume to be a vertex. Apart
-                                  // from the constructor that takes
-                                  // and stores the evaluation point,
-                                  // this class consists only of the
-                                  // function that implements
-                                  // assembling the right hand side.
-  template <int dim>
-  class PointValueEvaluation : public DualFunctionalBase<dim>
-  {
-    public:
-      PointValueEvaluation (const Point<dim> &evaluation_point);
-
-      virtual
-      void
-      assemble_rhs (const DoFHandler<dim> &dof_handler,
-                   Vector<double>        &rhs) const;
-      
-      DeclException1 (ExcEvaluationPointNotFound,
-                     Point<dim>,
-                     << "The evaluation point " << arg1
-                     << " was not found among the vertices of the present grid.");
-
-    
-      const Point<dim> evaluation_point;
-  };
-
-
-  template <int dim>
-  PointValueEvaluation<dim>::
-  PointValueEvaluation (const Point<dim> &evaluation_point)
-                 :
-                 evaluation_point (evaluation_point)
-  {}
-  
-
-                                  // As for doing the main purpose of
-                                  // the class, assembling the right
-                                  // hand side, let us first consider
-                                  // what is necessary: The right
-                                  // hand side of the dual problem is
-                                  // a vector of values J(phi_i),
-                                  // where J is the error functional,
-                                  // and phi_i is the i-th shape
-                                  // function. Here, J is the
-                                  // evaluation at the point x0,
-                                  // i.e. J(phi_i)=phi_i(x0).
-                                  //
-                                  // Now, we have assumed that the
-                                  // evaluation point is a
-                                  // vertex. Thus, for the usual
-                                  // finite elements we might be
-                                  // using in this program, we can
-                                  // take for granted that at such a
-                                  // point exactly one shape function
-                                  // is nonzero, and in particular
-                                  // has the value one. Thus, we set
-                                  // the right hand side vector to
-                                  // all-zeros, then seek for the
-                                  // shape function associated with
-                                  // that point and set the
-                                  // corresponding value of the right
-                                  // hand side vector to one:
-  template <int dim>
-  void
-  PointValueEvaluation<dim>::
-  assemble_rhs (const DoFHandler<dim> &dof_handler,
-               Vector<double>        &rhs) const
-  {
-                                    // So, first set everything to
-                                    // zeros...
-    rhs.reinit (dof_handler.n_dofs());
-
-                                    // ...then loop over cells and
-                                    // find the evaluation point
-                                    // among the vertices (or very
-                                    // close to a vertex, which may
-                                    // happen due to floating point
-                                    // round-off):
-    typename DoFHandler<dim>::active_cell_iterator
-      cell = dof_handler.begin_active(),
-      endc = dof_handler.end();
-    for (; cell!=endc; ++cell)
-      for (unsigned int vertex=0;
-          vertex<GeometryInfo<dim>::vertices_per_cell;
-          ++vertex)
-       if (cell->vertex(vertex).distance(evaluation_point)
-           < cell->diameter()*1e-8)
-         {
-                                            // Ok, found, so set
-                                            // corresponding entry,
-                                            // and leave function
-                                            // since we are finished:
-           rhs(cell->vertex_dof_index(vertex,0)) = 1;
-           return;
-         };
-
-                                    // Finally, a sanity check: if we
-                                    // somehow got here, then we must
-                                    // have missed the evaluation
-                                    // point, so raise an exception
-                                    // unconditionally:
-    AssertThrow (false, ExcEvaluationPointNotFound(evaluation_point));
-  }
-
-
-                                  // @sect4{The PointXDerivativeEvaluation class}
-  
-                                  // As second application, we again
-                                  // consider the evaluation of the
-                                  // x-derivative of the solution at
-                                  // one point. Again, the
-                                  // declaration of the class, and
-                                  // the implementation of its
-                                  // constructor is not too
-                                  // interesting:
-  template <int dim>
-  class PointXDerivativeEvaluation : public DualFunctionalBase<dim>
-  {
-    public:
-      PointXDerivativeEvaluation (const Point<dim> &evaluation_point);
-
-      virtual
-      void
-      assemble_rhs (const DoFHandler<dim> &dof_handler,
-                   Vector<double>        &rhs) const;
-      
-      DeclException1 (ExcEvaluationPointNotFound,
-                     Point<dim>,
-                     << "The evaluation point " << arg1
-                     << " was not found among the vertices of the present grid.");
-
-    
-      const Point<dim> evaluation_point;
-  };
-
-
-  template <int dim>
-  PointXDerivativeEvaluation<dim>::
-  PointXDerivativeEvaluation (const Point<dim> &evaluation_point)
-                 :
-                 evaluation_point (evaluation_point)
-  {}
-  
-
-                                  // What is interesting is the
-                                  // implementation of this
-                                  // functional: here,
-                                  // J(phi_i)=d/dx phi_i(x0).
-                                  //
-                                  // We could, as in the
-                                  // implementation of the respective
-                                  // evaluation object take the
-                                  // average of the gradients of each
-                                  // shape function phi_i at this
-                                  // evaluation point. However, we
-                                  // take a slightly different
-                                  // approach: we simply take the
-                                  // average over all cells that
-                                  // surround this point. The
-                                  // question which cells
-                                  // ``surrounds'' the evaluation
-                                  // point is made dependent on the
-                                  // mesh width by including those
-                                  // cells for which the distance of
-                                  // the cell's midpoint to the
-                                  // evaluation point is less than
-                                  // the cell's diameter.
-                                  //
-                                  // Taking the average of the
-                                  // gradient over the area/volume of
-                                  // these cells leads to a dual
-                                  // solution which is very close to
-                                  // the one which would result from
-                                  // the point evaluation of the
-                                  // gradient. It is simple to
-                                  // justify theoretically that this
-                                  // does not change the method
-                                  // significantly.
-  template <int dim>
-  void
-  PointXDerivativeEvaluation<dim>::
-  assemble_rhs (const DoFHandler<dim> &dof_handler,
-               Vector<double>        &rhs) const
-  {
-                                    // Again, first set all entries
-                                    // to zero:
-    rhs.reinit (dof_handler.n_dofs());
-
-                                    // Initialize a ``FEValues''
-                                    // object with a quadrature
-                                    // formula, have abbreviations
-                                    // for the number of quadrature
-                                    // points and shape functions...
-    QGauss<dim> quadrature(4);
-    MappingQ<dim> mapping (4);    
-    FEValues<dim>  fe_values (mapping, dof_handler.get_fe(), quadrature,
-                             update_gradients |
-                             update_q_points  |
-                             update_JxW_values);
-    const unsigned int n_q_points = fe_values.n_quadrature_points;
-    const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
-
-                                    // ...and have two objects that
-                                    // are used to store the global
-                                    // indices of the degrees of
-                                    // freedom on a cell, and the
-                                    // values of the gradients of the
-                                    // shape functions at the
-                                    // quadrature points:
-    Vector<double> cell_rhs (dofs_per_cell);
-    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
-
-                                    // Finally have a variable in
-                                    // which we will sum up the
-                                    // area/volume of the cells over
-                                    // which we integrate, by
-                                    // integrating the unit functions
-                                    // on these cells:
-    double total_volume = 0;
-    
-                                    // Then start the loop over all
-                                    // cells, and select those cells
-                                    // which are close enough to the
-                                    // evaluation point:
-    typename DoFHandler<dim>::active_cell_iterator
-      cell = dof_handler.begin_active(),
-      endc = dof_handler.end();
-    for (; cell!=endc; ++cell)
-      if (cell->center().distance(evaluation_point) <=
-         cell->diameter())
-       {
-                                          // If we have found such a
-                                          // cell, then initialize
-                                          // the ``FEValues'' object
-                                          // and integrate the
-                                          // x-component of the
-                                          // gradient of each shape
-                                          // function, as well as the
-                                          // unit function for the
-                                          // total area/volume.
-         fe_values.reinit (cell);
-         cell_rhs = 0;
-         
-         for (unsigned int q=0; q<n_q_points; ++q)
-           {
-             for (unsigned int i=0; i<dofs_per_cell; ++i)
-               cell_rhs(i) += fe_values.shape_grad(i,q)[0] *
-                              fe_values.JxW (q);
-             total_volume += fe_values.JxW (q);
-           };
-
-                                          // If we have the local
-                                          // contributions,
-                                          // distribute them to the
-                                          // global vector:
-         cell->get_dof_indices (local_dof_indices);
-         for (unsigned int i=0; i<dofs_per_cell; ++i)
-           rhs(local_dof_indices[i]) += cell_rhs(i);
-       };
-
-                                    // After we have looped over all
-                                    // cells, check whether we have
-                                    // found any at all, by making
-                                    // sure that their volume is
-                                    // non-zero. If not, then the
-                                    // results will be botched, as
-                                    // the right hand side should
-                                    // then still be zero, so throw
-                                    // an exception:
-    AssertThrow (total_volume > 0,
-                ExcEvaluationPointNotFound(evaluation_point));
-
-                                    // Finally, we have by now only
-                                    // integrated the gradients of
-                                    // the shape functions, not
-                                    // taking their mean value. We
-                                    // fix this by dividing by the
-                                    // measure of the volume over
-                                    // which we have integrated:
-    rhs.scale (1./total_volume);
-  }
-  
-
-}
-
-
-                                // @sect3{Extending the LaplaceSolver namespace}
-namespace LaplaceSolver
-{
-
-                                  // @sect4{The DualSolver class}
-
-                                  // In the same way as the
-                                  // ``PrimalSolver'' class above, we
-                                  // now implement a
-                                  // ``DualSolver''. It has all the
-                                  // same features, the only
-                                  // difference is that it does not
-                                  // take a function object denoting
-                                  // a right hand side object, but
-                                  // now takes a
-                                  // ``DualFunctionalBase'' object
-                                  // that will assemble the right
-                                  // hand side vector of the dual
-                                  // problem. The rest of the class
-                                  // is rather trivial.
-                                  //
-                                  // Since both primal and dual
-                                  // solver will use the same
-                                  // triangulation, but different
-                                  // discretizations, it now becomes
-                                  // clear why we have made the
-                                  // ``Base'' class a virtual one:
-                                  // since the final class will be
-                                  // derived from both
-                                  // ``PrimalSolver'' as well as
-                                  // ``DualSolver'', it would have
-                                  // two ``Base'' instances, would we
-                                  // not have marked the inheritance
-                                  // as virtual. Since in many
-                                  // applications the base class
-                                  // would store much more
-                                  // information than just the
-                                  // triangulation which needs to be
-                                  // shared between primal and dual
-                                  // solvers, we do not usually want
-                                  // to use two such base classes.
-  template <int dim>
-  class DualSolver : public Solver<dim>
-  {
-    public:
-      DualSolver (Triangulation<dim>       &triangulation,
-                 const FiniteElement<dim> &fe,
-                 const Quadrature<dim>    &quadrature,
-                 const Quadrature<dim-1>  &face_quadrature,
-                 const DualFunctional::DualFunctionalBase<dim> &dual_functional);
-
-      virtual
-      void
-      solve_problem ();
-      
-      virtual
-      unsigned int
-      n_dofs () const;
-
-      virtual
-      void
-      postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const;
-
-    
-      const SmartPointer<const DualFunctional::DualFunctionalBase<dim> > dual_functional;
-      virtual void assemble_rhs (Vector<double> &rhs) const;
-
-      static const ZeroFunction<dim> boundary_values;
-
-                                      // Same as above -- make a
-                                      // derived class a friend of
-                                      // this one:
-      friend class WeightedResidual<dim>;
-  };
-
-  template <int dim>
-  const ZeroFunction<dim> DualSolver<dim>::boundary_values;
-
-  template <int dim>
-  DualSolver<dim>::
-  DualSolver (Triangulation<dim>       &triangulation,
-             const FiniteElement<dim> &fe,
-             const Quadrature<dim>    &quadrature,
-             const Quadrature<dim-1>  &face_quadrature,
-             const DualFunctional::DualFunctionalBase<dim> &dual_functional)
-                 :
-                 Base<dim> (triangulation),
-                 Solver<dim> (triangulation, fe,
-                              quadrature, face_quadrature,
-                              boundary_values),
-                  dual_functional (&dual_functional)
-  {}
-
-
-  template <int dim>
-  void
-  DualSolver<dim>::solve_problem ()
-  {
-    Solver<dim>::solve_problem ();
-  }
-
-
-
-  template <int dim>
-  unsigned int
-  DualSolver<dim>::n_dofs() const
-  {
-    return Solver<dim>::n_dofs();
-  }
-
-
-  template <int dim>
-  void
-  DualSolver<dim>::
-  postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const
-  {
-    Solver<dim>::postprocess(postprocessor);
-  }
-  
-
-
-  template <int dim>
-  void
-  DualSolver<dim>::
-  assemble_rhs (Vector<double> &rhs) const 
-  {
-    dual_functional->assemble_rhs (this->dof_handler, rhs);
-  }
-
-
-                                  // @sect4{The WeightedResidual class}
-
-                                  // Here finally comes the main
-                                  // class of this program, the one
-                                  // that implements the dual
-                                  // weighted residual error
-                                  // estimator. It joins the primal
-                                  // and dual solver classes to use
-                                  // them for the computation of
-                                  // primal and dual solutions, and
-                                  // implements the error
-                                  // representation formula for use
-                                  // as error estimate and mesh
-                                  // refinement.
-                                  //
-                                  // The first few of the functions
-                                  // of this class are mostly
-                                  // overriders of the respective
-                                  // functions of the base class:
-  template <int dim>
-  class WeightedResidual : public PrimalSolver<dim>,
-                          public DualSolver<dim>
-  {
-    public:
-      WeightedResidual (Triangulation<dim>       &coarse_grid,
-                       const FiniteElement<dim> &primal_fe,
-                       const FiniteElement<dim> &dual_fe,
-                       const Quadrature<dim>    &quadrature,
-                       const Quadrature<dim-1>  &face_quadrature,
-                       const Function<dim>      &rhs_function,
-                       const Function<dim>      &boundary_values,
-                       const DualFunctional::DualFunctionalBase<dim> &dual_functional);
-
-      virtual
-      void
-      solve_problem ();
-
-      virtual
-      void
-      postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const;
-      
-      virtual
-      unsigned int
-      n_dofs () const;
-
-      virtual void refine_grid ();
-
-      virtual
-      void
-      output_solution () const;
-
-    
-                                      // In the private section, we
-                                      // have two functions that are
-                                      // used to call the
-                                      // ``solve_problem'' functions
-                                      // of the primal and dual base
-                                      // classes. These two functions
-                                      // will be called in parallel
-                                      // by the ``solve_problem''
-                                      // function of this class.
-      void solve_primal_problem ();
-      void solve_dual_problem ();
-                                      // Then declare abbreviations
-                                      // for active cell iterators,
-                                      // to avoid that we have to
-                                      // write this lengthy name
-                                      // over and over again:
-                                       
-      typedef
-      typename DoFHandler<dim>::active_cell_iterator
-      active_cell_iterator;
-
-                                      // Next, declare a data type
-                                      // that we will us to store the
-                                      // contribution of faces to the
-                                      // error estimator. The idea is
-                                      // that we can compute the face
-                                      // terms from each of the two
-                                      // cells to this face, as they
-                                      // are the same when viewed
-                                      // from both sides. What we
-                                      // will do is to compute them
-                                      // only once, based on some
-                                      // rules explained below which
-                                      // of the two adjacent cells
-                                      // will be in charge to do
-                                      // so. We then store the
-                                      // contribution of each face in
-                                      // a map mapping faces to their
-                                      // values, and only collect the
-                                      // contributions for each cell
-                                      // by looping over the cells a
-                                      // second time and grabbing the
-                                      // values from the map.
-                                      //
-                                      // The data type of this map is
-                                      // declared here:
-      typedef
-      typename std::map<typename DoFHandler<dim>::face_iterator,double>
-      FaceIntegrals;
-
-                                      // In the computation of the
-                                      // error estimates on cells and
-                                      // faces, we need a number of
-                                      // helper objects, such as
-                                      // ``FEValues'' and
-                                      // ``FEFaceValues'' functions,
-                                      // but also temporary objects
-                                      // storing the values and
-                                      // gradients of primal and dual
-                                      // solutions, for
-                                      // example. These fields are
-                                      // needed in the three
-                                      // functions that do the
-                                      // integration on cells, and
-                                      // regular and irregular faces,
-                                      // respectively.
-                                      //
-                                      // There are three reasonable
-                                      // ways to provide these
-                                      // fields: first, as local
-                                      // variables in the function
-                                      // that needs them; second, as
-                                      // member variables of this
-                                      // class; third, as arguments
-                                      // passed to that function.
-                                      //
-                                      // These three alternatives all
-                                      // have drawbacks: the third
-                                      // that their number is not
-                                      // neglectable and would make
-                                      // calling these functions a
-                                      // lengthy enterprise. The
-                                      // second has the drawback that
-                                      // it disallows
-                                      // parallelization, since the
-                                      // threads that will compute
-                                      // the error estimate have to
-                                      // have their own copies of
-                                      // these variables each, so
-                                      // member variables of the
-                                      // enclosing class will not
-                                      // work. The first approach,
-                                      // although straightforward,
-                                      // has a subtle but important
-                                      // drawback: we will call these
-                                      // functions over and over
-                                      // again, many thousand times
-                                      // maybe; it has now turned out
-                                      // that allocating vectors and
-                                      // other objects that need
-                                      // memory from the heap is an
-                                      // expensive business in terms
-                                      // of run-time, since memory
-                                      // allocation is expensive when
-                                      // several threads are
-                                      // involved. In our experience,
-                                      // more than 20 per cent of the
-                                      // total run time of error
-                                      // estimation functions are due
-                                      // to memory allocation, if
-                                      // done on a per-call level. It
-                                      // is thus significantly better
-                                      // to allocate the memory only
-                                      // once, and recycle the
-                                      // objects as often as
-                                      // possible.
-                                      //
-                                      // What to do? Our answer is to
-                                      // use a variant of the third
-                                      // strategy, namely generating
-                                      // these variables once in the
-                                      // main function of each
-                                      // thread, and passing them
-                                      // down to the functions that
-                                      // do the actual work. To avoid
-                                      // that we have to give these
-                                      // functions a dozen or so
-                                      // arguments, we pack all these
-                                      // variables into two
-                                      // structures, one which is
-                                      // used for the computations on
-                                      // cells, the other doing them
-                                      // on the faces. Instead of
-                                      // many individual objects, we
-                                      // will then only pass one such
-                                      // object to these functions,
-                                      // making their calling
-                                      // sequence simpler.
-      struct CellData
-      {
-         MappingQ<dim> mapping;
-         FEValues<dim>    fe_values;
-         const SmartPointer<const Function<dim> > right_hand_side;
-
-         std::vector<double> cell_residual;
-         std::vector<double> rhs_values;         
-         std::vector<double> dual_weights;       
-         typename std::vector<Tensor<2,dim> > cell_grad_grads;
-         CellData (const FiniteElement<dim> &fe,
-                   const Quadrature<dim>    &quadrature,
-                   const Function<dim>      &right_hand_side);
-      };
-
-      struct FaceData
-      {
-         MappingQ<dim> mapping;
-         FEFaceValues<dim>    fe_face_values_cell;
-         FEFaceValues<dim>    fe_face_values_neighbor;
-         FESubfaceValues<dim> fe_subface_values_cell;
-
-         std::vector<double> jump_residual;
-         std::vector<double> dual_weights;       
-         typename std::vector<Tensor<1,dim> > cell_grads;
-         typename std::vector<Tensor<1,dim> > neighbor_grads;
-         FaceData (const FiniteElement<dim> &fe,
-                   const Quadrature<dim-1>  &face_quadrature);
-      };
-
-      
-
-                                      // Regarding the evaluation of
-                                      // the error estimator, we have
-                                      // two driver functions that do
-                                      // this: the first is called to
-                                      // generate the cell-wise
-                                      // estimates, and splits up the
-                                      // task in a number of threads
-                                      // each of which work on a
-                                      // subset of the cells. The
-                                      // first function will run the
-                                      // second for each of these
-                                      // threads:
-      void estimate_error (Vector<float> &error_indicators) const;
-
-      void estimate_some (const Vector<double> &primal_solution,
-                         const Vector<double> &dual_weights,
-                         const unsigned int    n_threads,
-                         const unsigned int    this_thread,
-                         Vector<float>        &error_indicators,
-                         FaceIntegrals        &face_integrals) const;
-
-                                      // Then we have functions that
-                                      // do the actual integration of
-                                      // the error representation
-                                      // formula. They will treat the
-                                      // terms on the cell interiors,
-                                      // on those faces that have no
-                                      // hanging nodes, and on those
-                                      // faces with hanging nodes,
-                                      // respectively:
-      void
-      integrate_over_cell (const active_cell_iterator &cell,
-                          const unsigned int          cell_index,
-                          const Vector<double>       &primal_solution,
-                          const Vector<double>       &dual_weights,
-                          CellData                   &cell_data,
-                          Vector<float>              &error_indicators) const;
-
-      void
-      integrate_over_regular_face (const active_cell_iterator &cell,
-                                  const unsigned int          face_no,
-                                  const Vector<double>       &primal_solution,
-                                  const Vector<double>       &dual_weights,
-                                  FaceData                   &face_data,
-                                  FaceIntegrals              &face_integrals) const;
-      void
-      integrate_over_irregular_face (const active_cell_iterator &cell,
-                                    const unsigned int          face_no,
-                                    const Vector<double>       &primal_solution,
-                                    const Vector<double>       &dual_weights,
-                                    FaceData                   &face_data,
-                                    FaceIntegrals              &face_integrals) const;
-  };
-
-
-
-                                  // In the implementation of this
-                                  // class, we first have the
-                                  // constructors of the ``CellData''
-                                  // and ``FaceData'' member classes,
-                                  // and the ``WeightedResidual''
-                                  // constructor. They only
-                                  // initialize fields to their
-                                  // correct lengths, so we do not
-                                  // have to discuss them to length.
-  template <int dim>
-  WeightedResidual<dim>::CellData::
-  CellData (const FiniteElement<dim> &fe,
-           const Quadrature<dim>    &quadrature,
-           const Function<dim>      &right_hand_side)
-                 :
-                 mapping (4),
-                 fe_values (mapping, fe, quadrature,
-                            update_values             |
-                            update_second_derivatives |
-                            update_q_points           |
-                            update_JxW_values),
-                 right_hand_side (&right_hand_side)
-  {  
-    const unsigned int n_q_points
-      = quadrature.n_quadrature_points;
-  
-    cell_residual.resize(n_q_points);
-    rhs_values.resize(n_q_points);    
-    dual_weights.resize(n_q_points);    
-    cell_grad_grads.resize(n_q_points);
-  }
-  
-  
-
-  template <int dim>
-  WeightedResidual<dim>::FaceData::
-  FaceData (const FiniteElement<dim> &fe,
-           const Quadrature<dim-1>  &face_quadrature)
-                 :
-                 mapping (4),
-                 fe_face_values_cell (mapping, fe, face_quadrature,
-                                      update_values        |
-                                      update_gradients     |
-                                      update_JxW_values    |
-                                      update_normal_vectors),
-                 fe_face_values_neighbor (mapping, fe, face_quadrature,
-                                          update_values     |
-                                          update_gradients  |
-                                          update_JxW_values |
-                                          update_normal_vectors),
-                 fe_subface_values_cell (mapping, fe, face_quadrature,
-                                         update_gradients)
-  {  
-    const unsigned int n_face_q_points
-      = face_quadrature.n_quadrature_points;
-  
-    jump_residual.resize(n_face_q_points);
-    dual_weights.resize(n_face_q_points);    
-    cell_grads.resize(n_face_q_points);
-    neighbor_grads.resize(n_face_q_points);
-  }
-  
-
-
-
-  template <int dim>
-  WeightedResidual<dim>::
-  WeightedResidual (Triangulation<dim>       &coarse_grid,
-                   const FiniteElement<dim> &primal_fe,
-                   const FiniteElement<dim> &dual_fe,
-                   const Quadrature<dim>    &quadrature,
-                   const Quadrature<dim-1>  &face_quadrature,
-                   const Function<dim>      &rhs_function,
-                   const Function<dim>      &bv,
-                   const DualFunctional::DualFunctionalBase<dim> &dual_functional)
-                 :
-                 Base<dim> (coarse_grid),
-                  PrimalSolver<dim> (coarse_grid, primal_fe,
-                                    quadrature, face_quadrature,
-                                    rhs_function, bv),
-                  DualSolver<dim> (coarse_grid, dual_fe,
-                                  quadrature, face_quadrature,
-                                  dual_functional)
-  {}
-
-
-                                  // The next five functions are
-                                  // boring, as they simply relay
-                                  // their work to the base
-                                  // classes. The first calls the
-                                  // primal and dual solvers in
-                                  // parallel, while postprocessing
-                                  // the solution and retrieving the
-                                  // number of degrees of freedom is
-                                  // done by the primal class.
-  template <int dim>
-  void
-  WeightedResidual<dim>::solve_problem ()
-  {
-    Threads::ThreadGroup<> threads;
-    threads += Threads::spawn (*this, &WeightedResidual<dim>::solve_primal_problem)();
-    threads += Threads::spawn (*this, &WeightedResidual<dim>::solve_dual_problem)();
-    threads.join_all ();
-  }
-
-  
-  template <int dim>
-  void
-  WeightedResidual<dim>::solve_primal_problem ()
-  {
-    PrimalSolver<dim>::solve_problem ();
-  }
-
-  template <int dim>
-  void
-  WeightedResidual<dim>::solve_dual_problem ()
-  {
-    DualSolver<dim>::solve_problem ();
-  }
-  
-
-  template <int dim>
-  void
-  WeightedResidual<dim>::
-  postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const
-  {
-    PrimalSolver<dim>::postprocess (postprocessor);
-  }
-  
-  
-  template <int dim>
-  unsigned int
-  WeightedResidual<dim>::n_dofs () const
-  {
-    return PrimalSolver<dim>::n_dofs();
-  }
-
-
-
-                                  // Now, it is becoming more
-                                  // interesting: the ``refine_grid''
-                                  // function asks the error
-                                  // estimator to compute the
-                                  // cell-wise error indicators, then
-                                  // uses their absolute values for
-                                  // mesh refinement.
-  template <int dim>
-  void
-  WeightedResidual<dim>::refine_grid ()
-  {
-                                    // First call the function that
-                                    // computes the cell-wise and
-                                    // global error:
-    Vector<float> error_indicators (this->triangulation->n_active_cells());
-    estimate_error (error_indicators);
-
-                                    // Then note that marking cells
-                                    // for refinement or coarsening
-                                    // only works if all indicators
-                                    // are positive, to allow their
-                                    // comparison. Thus, drop the
-                                    // signs on all these indicators:
-    for (Vector<float>::iterator i=error_indicators.begin();
-        i != error_indicators.end(); ++i)
-      *i = std::fabs (*i);
-
-                                    // Finally, we can select between
-                                    // different strategies for
-                                    // refinement. The default here
-                                    // is to refine those cells with
-                                    // the largest error indicators
-                                    // that make up for a total of 80
-                                    // per cent of the error, while
-                                    // we coarsen those with the
-                                    // smallest indicators that make
-                                    // up for the bottom 2 per cent
-                                    // of the error.
-    GridRefinement::refine_and_coarsen_fixed_fraction (*this->triangulation,
-                                                      error_indicators,
-                                                      0.5/std::sqrt(1.+step), 
-                                                      0.2/std::sqrt(1.+step));
-    this->triangulation->execute_coarsening_and_refinement ();
-  }
-  
-
-                                  // Since we want to output both the
-                                  // primal and the dual solution, we
-                                  // overload the ``output_solution''
-                                  // function. The only interesting
-                                  // feature of this function is that
-                                  // the primal and dual solutions
-                                  // are defined on different finite
-                                  // element spaces, which is not the
-                                  // format the ``DataOut'' class
-                                  // expects. Thus, we have to
-                                  // transfer them to a common finite
-                                  // element space. Since we want the
-                                  // solutions only to see them
-                                  // qualitatively, we contend
-                                  // ourselves with interpolating the
-                                  // dual solution to the (smaller)
-                                  // primal space. For the
-                                  // interpolation, there is a
-                                  // library function, that takes a
-                                  // ``ConstraintMatrix'' object
-                                  // including the hanging node
-                                  // constraints. The rest is
-                                  // standard.
-                                  //
-                                  // There is, however, one
-                                  // work-around worth mentioning: in
-                                  // this function, as in a couple of
-                                  // following ones, we have to
-                                  // access the ``DoFHandler''
-                                  // objects and solutions of both
-                                  // the primal as well as of the
-                                  // dual solver. Since these are
-                                  // members of the ``Solver'' base
-                                  // class which exists twice in the
-                                  // class hierarchy leading to the
-                                  // present class (once as base
-                                  // class of the ``PrimalSolver''
-                                  // class, once as base class of the
-                                  // ``DualSolver'' class), we have
-                                  // to disambiguate accesses to them
-                                  // by telling the compiler a member
-                                  // of which of these two instances
-                                  // we want to access. The way to do
-                                  // this would be identify the
-                                  // member by pointing a path
-                                  // through the class hierarchy
-                                  // which disambiguates the base
-                                  // class, for example writing
-                                  // ``PrimalSolver::dof_handler'' to
-                                  // denote the member variable
-                                  // ``dof_handler'' from the
-                                  // ``Solver'' base class of the
-                                  // ``PrimalSolver''
-                                  // class. Unfortunately, this
-                                  // confuses gcc's version 2.96 (a
-                                  // version that was intended as a
-                                  // development snapshot, but
-                                  // delivered as system compiler by
-                                  // Red Hat in their 7.x releases)
-                                  // so much that it bails out and
-                                  // refuses to compile the code.
-                                  //
-                                  // Thus, we have to work around
-                                  // this problem. We do this by
-                                  // introducing references to the
-                                  // ``PrimalSolver'' and
-                                  // ``DualSolver'' components of the
-                                  // ``WeightedResidual'' object at
-                                  // the beginning of the
-                                  // function. Since each of these
-                                  // has an unambiguous base class
-                                  // ``Solver'', we can access the
-                                  // member variables we want through
-                                  // these references. However, we
-                                  // are now accessing protected
-                                  // member variables of these
-                                  // classes through a pointer other
-                                  // than the ``this'' pointer (in
-                                  // fact, this is of course the
-                                  // ``this'' pointer, but not
-                                  // explicitly). This finally is the
-                                  // reason why we had to declare the
-                                  // present class a friend of the
-                                  // classes we so access.
-  template <int dim>
-  void
-  WeightedResidual<dim>::output_solution () const
-  {
-    const PrimalSolver<dim> &primal_solver = *this;
-    const DualSolver<dim>   &dual_solver   = *this;
-    
-    ConstraintMatrix primal_hanging_node_constraints;
-    DoFTools::make_hanging_node_constraints (primal_solver.dof_handler,
-                                            primal_hanging_node_constraints);
-    primal_hanging_node_constraints.close();
-    Vector<double> dual_solution (primal_solver.dof_handler.n_dofs());
-    FETools::interpolate (dual_solver.dof_handler,
-                         dual_solver.solution,
-                         primal_solver.dof_handler,
-                         primal_hanging_node_constraints,
-                         dual_solution);    
-
-                                    // approximate error, gradient,
-                                    // and second derivative
-                                    // information as cell information
-    Vector<float> error_indicators (this->triangulation->n_active_cells());
-    Vector<float> gradient_indicators (this->triangulation->n_active_cells());
-    Vector<float> second_indicators (this->triangulation->n_active_cells());
-    {
-      MappingQ<dim> mapping(4);
-      KellyErrorEstimator<dim>::estimate (mapping, primal_solver.dof_handler,
-                                         QGauss<dim-1>(3),
-                                         typename FunctionMap<dim>::type(),
-                                         primal_solver.solution,
-                                         error_indicators);
-
-      DerivativeApproximation::
-       approximate_gradient (mapping,
-                             primal_solver.dof_handler,
-                             primal_solver.solution,
-                             gradient_indicators);
-      
-      DerivativeApproximation::
-       approximate_second_derivative (mapping,
-                                      primal_solver.dof_handler,
-                                      primal_solver.solution,
-                                      second_indicators);
-      
-    }
-                                    // distribute cell to dof vectors
-    Vector<double> x_error_indicators (primal_solver.dof_handler.n_dofs());
-    Vector<double> x_gradient_indicators (primal_solver.dof_handler.n_dofs());
-    Vector<double> x_second_indicators (primal_solver.dof_handler.n_dofs());
-    DoFTools::distribute_cell_to_dof_vector (primal_solver.dof_handler,
-                                            error_indicators,
-                                            x_error_indicators);
-    DoFTools::distribute_cell_to_dof_vector (primal_solver.dof_handler,
-                                            gradient_indicators,
-                                            x_gradient_indicators);
-    DoFTools::distribute_cell_to_dof_vector (primal_solver.dof_handler,
-                                            second_indicators,
-                                            x_second_indicators);
-    
-    
-
-                                    // we generate too much output in
-                                    // 3d. instead of doing it that
-                                    // way, simply generate a coarser
-                                    // mesh and output from there
-    Triangulation<dim> coarser_mesh;
-    coarser_mesh.copy_triangulation (*this->triangulation);
-    for (typename Triangulation<dim>::active_cell_iterator
-          cell = coarser_mesh.begin_active();
-        cell != coarser_mesh.end(); ++cell)
-      cell->set_coarsen_flag();
-    coarser_mesh.execute_coarsening_and_refinement ();
-
-                                    // next generate a DoF handler on
-                                    // that mesh and a map fron one
-                                    // to the other mesh
-    DoFHandler<dim> coarser_dof_handler (coarser_mesh);
-    coarser_dof_handler.distribute_dofs (primal_solver.dof_handler.get_fe());
-    InterGridMap<DoFHandler<dim> > coarse_to_fine_map;
-    coarse_to_fine_map.make_mapping (coarser_dof_handler,
-                                    primal_solver.dof_handler);
-
-                                    // finally we have to transfer
-                                    // the data vectors
-    Vector<double> coarse_primal_solution (coarser_dof_handler.n_dofs());
-    Vector<double> coarse_dual_solution (coarser_dof_handler.n_dofs());
-    Vector<double> coarse_error_indicators (coarser_dof_handler.n_dofs());
-    Vector<double> coarse_gradient_indicators (coarser_dof_handler.n_dofs());
-    Vector<double> coarse_second_indicators (coarser_dof_handler.n_dofs());
-
-    Vector<double> tmp (coarser_dof_handler.get_fe().dofs_per_cell);
-    for (typename DoFHandler<dim>::active_cell_iterator
-          cell = coarser_dof_handler.begin_active();
-        cell != coarser_dof_handler.end(); ++cell)
-      {
-       coarse_to_fine_map[cell]->get_interpolated_dof_values (primal_solver.solution,tmp);
-       cell->set_dof_values (tmp, coarse_primal_solution);
-       
-       coarse_to_fine_map[cell]->get_interpolated_dof_values (dual_solution,tmp);
-       cell->set_dof_values (tmp, coarse_dual_solution);
-       
-       coarse_to_fine_map[cell]->get_interpolated_dof_values (x_error_indicators,tmp);
-       cell->set_dof_values (tmp, coarse_error_indicators);
-       
-       coarse_to_fine_map[cell]->get_interpolated_dof_values (x_gradient_indicators,tmp);
-       cell->set_dof_values (tmp, coarse_gradient_indicators);
-       
-       coarse_to_fine_map[cell]->get_interpolated_dof_values (x_second_indicators,tmp);
-       cell->set_dof_values (tmp, coarse_second_indicators);
-      }
-       
-    {
-      DataOut<dim> data_out;
-      data_out.attach_dof_handler (coarser_dof_handler);
-      data_out.add_data_vector (coarse_primal_solution, "primal_solution");
-      data_out.add_data_vector (coarse_dual_solution, "dual_solution");
-      data_out.add_data_vector (coarse_error_indicators, "errors");
-      data_out.add_data_vector (coarse_gradient_indicators, "gradient");
-      data_out.add_data_vector (coarse_second_indicators, "second_derivatives");
-      data_out.build_patches ();
-  
-#ifdef HAVE_STD_STRINGSTREAM
-      std::ostringstream filename;
-#else
-      std::ostrstream filename;
-#endif
-      filename << "spec2006-447.dealII/"
-              << "solution-"
-              << this->refinement_cycle
-              << ".gmv"
-              << std::ends;
-#ifdef HAVE_STD_STRINGSTREAM
-      std::ofstream out (filename.str().c_str());
-#else
-      std::ofstream out (filename.str());
-#endif
-    
-      data_out.write_gmv (out);
-    }
-    
-  }
-
-
-                                  // @sect3{Estimating errors}
-
-                                  // @sect4{Error estimation driver functions}
-                                  //
-                                  // As for the actual computation of
-                                  // error estimates, let's start
-                                  // with the function that drives
-                                  // all this, i.e. calls those
-                                  // functions that actually do the
-                                  // work, and finally collects the
-                                  // results.
-  
-  template <int dim>
-  void
-  WeightedResidual<dim>::
-  estimate_error (Vector<float> &error_indicators) const
-  {
-    const PrimalSolver<dim> &primal_solver = *this;
-    const DualSolver<dim>   &dual_solver   = *this;
-
-                                    // The first task in computing
-                                    // the error is to set up vectors
-                                    // that denote the primal
-                                    // solution, and the weights
-                                    // (z-z_h)=(z-I_hz), both in the
-                                    // finite element space for which
-                                    // we have computed the dual
-                                    // solution. For this, we have to
-                                    // interpolate the primal
-                                    // solution to the dual finite
-                                    // element space, and to subtract
-                                    // the interpolation of the
-                                    // computed dual solution to the
-                                    // primal finite element
-                                    // space. Fortunately, the
-                                    // library provides functions for
-                                    // the interpolation into larger
-                                    // or smaller finite element
-                                    // spaces, so this is mostly
-                                    // obvious.
-                                    //
-                                    // First, let's do that for the
-                                    // primal solution: it is
-                                    // cell-wise interpolated into
-                                    // the finite element space in
-                                    // which we have solved the dual
-                                    // problem: But, again as in the
-                                    // ``WeightedResidual::output_solution''
-                                    // function we first need to
-                                    // create a ConstraintMatrix
-                                    // including the hanging node
-                                    // constraints, but this time of
-                                    // the dual finite element space.
-    ConstraintMatrix dual_hanging_node_constraints;
-    DoFTools::make_hanging_node_constraints (dual_solver.dof_handler,
-                                            dual_hanging_node_constraints);
-    dual_hanging_node_constraints.close();
-    Vector<double> primal_solution (dual_solver.dof_handler.n_dofs());
-    FETools::interpolate (primal_solver.dof_handler,
-                         primal_solver.solution,
-                         dual_solver.dof_handler,
-                         dual_hanging_node_constraints,
-                         primal_solution);
-    
-                                    // Then for computing the
-                                    // interpolation of the
-                                    // numerically approximated dual
-                                    // solution z into the finite
-                                    // element space of the primal
-                                    // solution and subtracting it
-                                    // from z: use the
-                                    // ``interpolate_difference''
-                                    // function, that gives (z-I_hz)
-                                    // in the element space of the
-                                    // dual solution.
-    ConstraintMatrix primal_hanging_node_constraints;
-    DoFTools::make_hanging_node_constraints (primal_solver.dof_handler,
-                                            primal_hanging_node_constraints);
-    primal_hanging_node_constraints.close();
-    Vector<double> dual_weights (dual_solver.dof_handler.n_dofs());
-    FETools::interpolation_difference (dual_solver.dof_handler,
-                                      dual_hanging_node_constraints,
-                                      dual_solver.solution,
-                                      primal_solver.dof_handler,
-                                      primal_hanging_node_constraints,
-                                      dual_weights);
-    
-                                    // Note that this could probably
-                                    // have been more efficient since
-                                    // those constraints have been
-                                    // used previously when
-                                    // assembling matrix and right
-                                    // hand side for the primal
-                                    // problem and writing out the
-                                    // dual solution. We leave the
-                                    // optimization of the program in
-                                    // this respect as an exercise.
-    
-                                    // Having computed the dual
-                                    // weights we now proceed with
-                                    // computing the cell and face
-                                    // residuals of the primal
-                                    // solution. First we set up a
-                                    // map between face iterators and
-                                    // their jump term contributions
-                                    // of faces to the error
-                                    // estimator. The reason is that
-                                    // we compute the jump terms only
-                                    // once, from one side of the
-                                    // face, and want to collect them
-                                    // only afterwards when looping
-                                    // over all cells a second time.
-                                    //
-                                    // We initialize this map already
-                                    // with a value of -1e20 for all
-                                    // faces, since this value will
-                                    // strike in the results if
-                                    // something should go wrong and
-                                    // we fail to compute the value
-                                    // for a face for some
-                                    // reason. Secondly, we
-                                    // initialize the map once before
-                                    // we branch to different threads
-                                    // since this way the map's
-                                    // structure is no more modified
-                                    // by the individual threads,
-                                    // only existing entries are set
-                                    // to new values. This relieves
-                                    // us from the necessity to
-                                    // synchronise the threads
-                                    // through a mutex each time they
-                                    // write to (and modify the
-                                    // structure of) this map.
-    FaceIntegrals face_integrals;
-    for (active_cell_iterator cell=dual_solver.dof_handler.begin_active();
-        cell!=dual_solver.dof_handler.end();
-        ++cell)
-      for (unsigned int face_no=0;
-          face_no<GeometryInfo<dim>::faces_per_cell;
-          ++face_no)
-       face_integrals[cell->face(face_no)] = -1e20;
-
-                                    // Then set up a vector with
-                                    // error indicators.  Reserve one
-                                    // slot for each cell and set it
-                                    // to zero.
-    error_indicators.reinit (dual_solver.dof_handler
-                            .get_tria().n_active_cells());
-
-                                    // Now start a number of threads
-                                    // which compute the error
-                                    // formula on parts of all the
-                                    // cells, and once they are all
-                                    // started wait until they have
-                                    // all finished:
-    const unsigned int n_threads = multithread_info.n_default_threads;
-    Threads::ThreadGroup<> threads;
-    for (unsigned int i=0; i<n_threads; ++i)
-      threads += Threads::spawn (*this, &WeightedResidual<dim>::estimate_some)
-                 (primal_solution,
-                  dual_weights,
-                  n_threads, i,
-                  error_indicators,
-                  face_integrals);
-    threads.join_all();    
-
-                                    // Once the error contributions
-                                    // are computed, sum them up. For
-                                    // this, note that the cell terms
-                                    // are already set, and that only
-                                    // the edge terms need to be
-                                    // collected. Thus, loop over all
-                                    // cells and their faces, make
-                                    // sure that the contributions of
-                                    // each of the faces are there,
-                                    // and add them up. Only take
-                                    // minus one half of the jump
-                                    // term, since the other half
-                                    // will be taken by the
-                                    // neighboring cell.
-    unsigned int present_cell=0;  
-    for (active_cell_iterator cell=dual_solver.dof_handler.begin_active();
-        cell!=dual_solver.dof_handler.end();
-        ++cell, ++present_cell)
-      for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell;
-          ++face_no)
-       {
-         Assert(face_integrals.find(cell->face(face_no)) !=
-                face_integrals.end(),
-                ExcInternalError());
-         error_indicators(present_cell)
-           -= 0.5*face_integrals[cell->face(face_no)];
-       };
-    std::cout << "   Estimated error="
-             << std::accumulate (error_indicators.begin(),
-                                 error_indicators.end(), 0.)
-             << std::endl;
-  }
-
-
-                                  // @sect4{Estimating on a subset of cells}
-
-                                  // Next we have the function that
-                                  // is called to estimate the error
-                                  // on a subset of cells. The
-                                  // function may be called multiply
-                                  // if the library was configured to
-                                  // use multi-threading. Here it
-                                  // goes:
-  template <int dim>
-  void
-  WeightedResidual<dim>::
-  estimate_some (const Vector<double> &primal_solution,
-                const Vector<double> &dual_weights,
-                const unsigned int    n_threads,
-                const unsigned int    this_thread,
-                Vector<float>        &error_indicators,
-                FaceIntegrals        &face_integrals) const
-  {
-    const PrimalSolver<dim> &primal_solver = *this;
-    const DualSolver<dim>   &dual_solver   = *this;
-
-                                    // At the beginning, we
-                                    // initialize two variables for
-                                    // each thread which may be
-                                    // running this function. The
-                                    // reason for these functions was
-                                    // discussed above, when the
-                                    // respective classes were
-                                    // discussed, so we here only
-                                    // point out that since they are
-                                    // local to the function that is
-                                    // spawned when running more than
-                                    // one thread, the data of these
-                                    // objects exists actually once
-                                    // per thread, so we don't have
-                                    // to take care about
-                                    // synchronising access to them.
-    CellData cell_data (*dual_solver.fe,
-                       *dual_solver.quadrature,
-                       *primal_solver.rhs_function);
-    FaceData face_data (*dual_solver.fe,
-                       *dual_solver.face_quadrature);    
-
-                                    // Then calculate the start cell
-                                    // for this thread. We let the
-                                    // different threads run on
-                                    // interleaved cells, i.e. for
-                                    // example if we have 4 threads,
-                                    // then the first thread treates
-                                    // cells 0, 4, 8, etc, while the
-                                    // second threads works on cells 1,
-                                    // 5, 9, and so on. The reason is
-                                    // that it takes vastly more time
-                                    // to work on cells with hanging
-                                    // nodes than on regular cells, but
-                                    // such cells are not evenly
-                                    // distributed across the range of
-                                    // cell iterators, so in order to
-                                    // have the different threads do
-                                    // approximately the same amount of
-                                    // work, we have to let them work
-                                    // interleaved to the effect of a
-                                    // pseudorandom distribution of the
-                                    // `hard' cells to the different
-                                    // threads.
-    active_cell_iterator cell=dual_solver.dof_handler.begin_active();
-    for (unsigned int t=0;
-        (t<this_thread) && (cell!=dual_solver.dof_handler.end());
-        ++t, ++cell);
-
-                                    // If there are no cells for this
-                                    // thread (for example if there
-                                    // are a total of less cells than
-                                    // there are threads), then go
-                                    // back right now
-    if (cell == dual_solver.dof_handler.end())
-      return;
-    
-                                    // Next loop over all cells. The
-                                    // check for loop end is done at
-                                    // the end of the loop, along
-                                    // with incrementing the loop
-                                    // index.
-    for (unsigned int cell_index=this_thread; true; )
-      {
-                                        // First task on each cell is
-                                        // to compute the cell
-                                        // residual contributions of
-                                        // this cell, and put them
-                                        // into the
-                                        // ``error_indicators''
-                                        // variable:
-       integrate_over_cell (cell, cell_index,
-                            primal_solution,
-                            dual_weights,
-                            cell_data,
-                            error_indicators);
-       
-                                        // After computing the cell
-                                        // terms, turn to the face
-                                        // terms. For this, loop over
-                                        // all faces of the present
-                                        // cell, and see whether
-                                        // something needs to be
-                                        // computed on it:
-       for (unsigned int face_no=0;
-            face_no<GeometryInfo<dim>::faces_per_cell;
-            ++face_no)
-         {
-                                            // First, if this face is
-                                            // part of the boundary,
-                                            // then there is nothing
-                                            // to do. However, to
-                                            // make things easier
-                                            // when summing up the
-                                            // contributions of the
-                                            // faces of cells, we
-                                            // enter this face into
-                                            // the list of faces with
-                                            // a zero contribution to
-                                            // the error.
-           if (cell->face(face_no)->at_boundary()) 
-             {
-               face_integrals[cell->face(face_no)] = 0;
-               continue;
-             };
-           
-                                            // Next, note that since
-                                            // we want to compute the
-                                            // jump terms on each
-                                            // face only once
-                                            // although we access it
-                                            // twice (if it is not at
-                                            // the boundary), we have
-                                            // to define some rules
-                                            // who is responsible for
-                                            // computing on a face:
-                                            //
-                                            // First, if the
-                                            // neighboring cell is on
-                                            // the same level as this
-                                            // one, i.e. neither
-                                            // further refined not
-                                            // coarser, then the one
-                                            // with the lower index
-                                            // within this level does
-                                            // the work. In other
-                                            // words: if the other
-                                            // one has a lower index,
-                                            // then skip work on this
-                                            // face:
-           if ((cell->neighbor(face_no)->has_children() == false) &&
-               (cell->neighbor(face_no)->level() == cell->level()) &&
-               (cell->neighbor(face_no)->index() < cell->index()))
-             continue;
-
-                                            // Likewise, we always
-                                            // work from the coarser
-                                            // cell if this and its
-                                            // neighbor differ in
-                                            // refinement. Thus, if
-                                            // the neighboring cell
-                                            // is less refined than
-                                            // the present one, then
-                                            // do nothing since we
-                                            // integrate over the
-                                            // subfaces when we visit
-                                            // the coarse cell.
-           if (cell->at_boundary(face_no) == false)
-             if (cell->neighbor(face_no)->level() < cell->level())
-               continue;         
-
-
-                                            // Now we know that we
-                                            // are in charge here, so
-                                            // actually compute the
-                                            // face jump terms. If
-                                            // the face is a regular
-                                            // one, i.e.  the other
-                                            // side's cell is neither
-                                            // coarser not finer than
-                                            // this cell, then call
-                                            // one function, and if
-                                            // the cell on the other
-                                            // side is further
-                                            // refined, then use
-                                            // another function. Note
-                                            // that the case that the
-                                            // cell on the other side
-                                            // is coarser cannot
-                                            // happen since we have
-                                            // decided above that we
-                                            // handle this case when
-                                            // we pass over that
-                                            // other cell.
-           if (cell->face(face_no)->has_children() == false)
-             integrate_over_regular_face (cell, face_no,
-                                          primal_solution,
-                                          dual_weights,
-                                          face_data,
-                                          face_integrals);       
-           else
-             integrate_over_irregular_face (cell, face_no,
-                                            primal_solution,
-                                            dual_weights,
-                                            face_data,
-                                            face_integrals);
-         };
-
-                                        // After computing the cell
-                                        // contributions and looping
-                                        // over the faces, go to the
-                                        // next cell for this
-                                        // thread. Note again that
-                                        // the cells for each of the
-                                        // threads are interleaved.
-                                        // If we are at the end of
-                                        // our workload, jump out
-                                        // of the loop.
-       for (unsigned int t=0;
-            ((t<n_threads) && (cell!=dual_solver.dof_handler.end()));
-            ++t, ++cell, ++cell_index);
-       if (cell == dual_solver.dof_handler.end())
-         break;
-      };
-  }
-
-
-                                  // @sect4{Computing cell term error contributions}
-
-                                  // As for the actual computation of
-                                  // the error contributions, first
-                                  // turn to the cell terms:
-  template <int dim>
-  void WeightedResidual<dim>::
-  integrate_over_cell (const active_cell_iterator &cell,
-                      const unsigned int          cell_index,
-                      const Vector<double>       &primal_solution,
-                      const Vector<double>       &dual_weights,
-                      CellData                   &cell_data,
-                      Vector<float>              &error_indicators) const
-  {
-                                    // The tasks to be done are what
-                                    // appears natural from looking
-                                    // at the error estimation
-                                    // formula: first compute the the
-                                    // right hand side and the
-                                    // Laplacian of the numerical
-                                    // solution at the quadrature
-                                    // points for the cell residual,
-    cell_data.fe_values.reinit (cell);
-    cell_data.right_hand_side
-      ->value_list (cell_data.fe_values.get_quadrature_points(),
-                   cell_data.rhs_values);
-    cell_data.fe_values.get_function_2nd_derivatives (primal_solution,
-                                                     cell_data.cell_grad_grads);
-
-                                    // ...then get the dual weights...
-    cell_data.fe_values.get_function_values (dual_weights,
-                                            cell_data.dual_weights);
-
-                                    // ...and finally build the sum
-                                    // over all quadrature points and
-                                    // store it with the present
-                                    // cell:
-    double sum = 0;
-    for (unsigned int p=0; p<cell_data.fe_values.n_quadrature_points; ++p)
-      sum += ((cell_data.rhs_values[p]+trace(cell_data.cell_grad_grads[p])) *
-             cell_data.dual_weights[p] *
-             cell_data.fe_values.JxW (p));
-    error_indicators(cell_index) += sum;
-  }
-
-
-                                  // @sect4{Computing edge term error contributions - 1}
-  
-                                  // On the other hand, computation
-                                  // of the edge terms for the error
-                                  // estimate is not so
-                                  // simple. First, we have to
-                                  // distinguish between faces with
-                                  // and without hanging
-                                  // nodes. Because it is the simple
-                                  // case, we first consider the case
-                                  // without hanging nodes on a face
-                                  // (let's call this the `regular'
-                                  // case):
-  template <int dim>
-  void WeightedResidual<dim>::
-  integrate_over_regular_face (const active_cell_iterator &cell,
-                              const unsigned int          face_no,
-                              const Vector<double>       &primal_solution,
-                              const Vector<double>       &dual_weights,
-                              FaceData                   &face_data,
-                              FaceIntegrals              &face_integrals) const
-  {
-    const unsigned int
-      n_q_points = face_data.fe_face_values_cell.n_quadrature_points;
-
-                                    // The first step is to get the
-                                    // values of the gradients at the
-                                    // quadrature points of the
-                                    // finite element field on the
-                                    // present cell. For this,
-                                    // initialize the
-                                    // ``FEFaceValues'' object
-                                    // corresponding to this side of
-                                    // the face, and extract the
-                                    // gradients using that
-                                    // object.
-    face_data.fe_face_values_cell.reinit (cell, face_no);
-    face_data.fe_face_values_cell.get_function_grads (primal_solution,
-                                                     face_data.cell_grads);
-
-                                    // The second step is then to
-                                    // extract the gradients of the
-                                    // finite element solution at the
-                                    // quadrature points on the other
-                                    // side of the face, i.e. from
-                                    // the neighboring cell.
-                                    //
-                                    // For this, do a sanity check
-                                    // before: make sure that the
-                                    // neigbor actually exists (yes,
-                                    // we should not have come here
-                                    // if the neighbor did not exist,
-                                    // but in complicated software
-                                    // there are bugs, so better
-                                    // check this), and if this is
-                                    // not the case throw an error.
-    Assert (cell->neighbor(face_no).state() == IteratorState::valid,
-           ExcInternalError());
-                                    // If we have that, then we need
-                                    // to find out with which face of
-                                    // the neighboring cell we have
-                                    // to work, i.e. the
-                                    // ``home-many''the neighbor the
-                                    // present cell is of the cell
-                                    // behind the present face. For
-                                    // this, there is a function, and
-                                    // we put the result into a
-                                    // variable with the name
-                                    // ``neighbor_neighbor'':
-    const unsigned int
-      neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
-                                    // Then define an abbreviation
-                                    // for the neigbor cell,
-                                    // initialize the
-                                    // ``FEFaceValues'' object on
-                                    // that cell, and extract the
-                                    // gradients on that cell:
-    const active_cell_iterator neighbor = cell->neighbor(face_no);
-    face_data.fe_face_values_neighbor.reinit (neighbor, neighbor_neighbor);      
-    face_data.fe_face_values_neighbor.get_function_grads (primal_solution,
-                                                         face_data.neighbor_grads);
-
-                                    // Now that we have the gradients
-                                    // on this and the neighboring
-                                    // cell, compute the jump
-                                    // residual by multiplying the
-                                    // jump in the gradient with the
-                                    // normal vector:
-    for (unsigned int p=0; p<n_q_points; ++p)
-      face_data.jump_residual[p]
-       = ((face_data.cell_grads[p] - face_data.neighbor_grads[p]) *
-          face_data.fe_face_values_cell.normal_vector(p));
-
-                                    // Next get the dual weights for
-                                    // this face:
-    face_data.fe_face_values_cell.get_function_values (dual_weights,
-                                                      face_data.dual_weights);
-    
-                                    // Finally, we have to compute
-                                    // the sum over jump residuals,
-                                    // dual weights, and quadrature
-                                    // weights, to get the result for
-                                    // this face:
-    double face_integral = 0;
-    for (unsigned int p=0; p<n_q_points; ++p)
-      face_integral += (face_data.jump_residual[p] *
-                       face_data.dual_weights[p]  *
-                       face_data.fe_face_values_cell.JxW(p));
-
-                                    // Double check that the element
-                                    // already exists and that it was
-                                    // not already written to...
-    Assert (face_integrals.find (cell->face(face_no)) != face_integrals.end(),
-           ExcInternalError());
-    Assert (face_integrals[cell->face(face_no)] == -1e20,
-           ExcInternalError());
-
-                                    // ...then store computed value
-                                    // at assigned location. Note
-                                    // that the stored value does not
-                                    // contain the factor 1/2 that
-                                    // appears in the error
-                                    // representation. The reason is
-                                    // that the term actually does
-                                    // not have this factor if we
-                                    // loop over all faces in the
-                                    // triangulation, but only
-                                    // appears if we write it as a
-                                    // sum over all cells and all
-                                    // faces of each cell; we thus
-                                    // visit the same face twice. We
-                                    // take account of this by using
-                                    // this factor -1/2 later, when we
-                                    // sum up the contributions for
-                                    // each cell individually.
-    face_integrals[cell->face(face_no)] = face_integral;
-  }
-
-
-                                  // @sect4{Computing edge term error contributions - 2}
-  
-                                  // We are still missing the case of
-                                  // faces with hanging nodes. This
-                                  // is what is covered in this
-                                  // function:
-  template <int dim>
-  void WeightedResidual<dim>::
-  integrate_over_irregular_face (const active_cell_iterator &cell,
-                                const unsigned int          face_no,
-                                const Vector<double>       &primal_solution,
-                                const Vector<double>       &dual_weights,
-                                FaceData                   &face_data,
-                                FaceIntegrals              &face_integrals) const
-  {
-                                    // First again two abbreviations,
-                                    // and some consistency checks
-                                    // whether the function is called
-                                    // only on faces for which it is
-                                    // supposed to be called:
-    const unsigned int
-      n_q_points = face_data.fe_face_values_cell.n_quadrature_points;
-
-    const typename DoFHandler<dim>::cell_iterator
-      neighbor = cell->neighbor(face_no);    
-    Assert (neighbor.state() == IteratorState::valid,
-           ExcInternalError());
-    Assert (neighbor->has_children(),
-           ExcInternalError());
-
-                                    // Then find out which neighbor
-                                    // the present cell is of the
-                                    // adjacent cell. Note that we
-                                    // will operator on the children
-                                    // of this adjacent cell, but
-                                    // that their orientation is the
-                                    // same as that of their mother,
-                                    // i.e. the neigbor direction is
-                                    // the same.
-    const unsigned int
-      neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
-  
-                                    // Then simply do everything we
-                                    // did in the previous function
-                                    // for one face for all the
-                                    // sub-faces now:
-    for (unsigned int subface_no=0;
-        subface_no<GeometryInfo<dim>::subfaces_per_face;
-        ++subface_no)
-      {
-       const active_cell_iterator neighbor_child
-          = cell->neighbor_child_on_subface (face_no, subface_no);
-
-                                        // Now start the work by
-                                        // again getting the gradient
-                                        // of the solution first at
-                                        // this side of the
-                                        // interface,
-       face_data.fe_subface_values_cell.reinit (cell, face_no, subface_no);
-       face_data.fe_subface_values_cell.get_function_grads (primal_solution,
-                                                            face_data.cell_grads);
-                                        // then at the other side,
-       face_data.fe_face_values_neighbor.reinit (neighbor_child,
-                                            neighbor_neighbor);
-       face_data.fe_face_values_neighbor.get_function_grads (primal_solution,
-                                                             face_data.neighbor_grads);
-      
-                                        // and finally building the
-                                        // jump residuals. Since we
-                                        // take the normal vector
-                                        // from the other cell this
-                                        // time, revert the sign of
-                                        // the first term compared to
-                                        // the other function:
-       for (unsigned int p=0; p<n_q_points; ++p)
-         face_data.jump_residual[p]
-            = ((face_data.neighbor_grads[p] - face_data.cell_grads[p]) *
-               face_data.fe_face_values_neighbor.normal_vector(p));
-
-                                        // Then get dual weights:
-       face_data.fe_face_values_neighbor.get_function_values (dual_weights,
-                                                              face_data.dual_weights);
-       
-                                        // At last, sum up the
-                                        // contribution of this
-                                        // sub-face, and set it in
-                                        // the global map:
-       double face_integral = 0;
-       for (unsigned int p=0; p<n_q_points; ++p)
-         face_integral += (face_data.jump_residual[p] *
-                           face_data.dual_weights[p] *
-                           face_data.fe_face_values_neighbor.JxW(p));
-       face_integrals[neighbor_child->face(neighbor_neighbor)]
-         = face_integral;
-      };
-
-                                    // Once the contributions of all
-                                    // sub-faces are computed, loop
-                                    // over all sub-faces to collect
-                                    // and store them with the mother
-                                    // face for simple use when later
-                                    // collecting the error terms of
-                                    // cells. Again make safety
-                                    // checks that the entries for
-                                    // the sub-faces have been
-                                    // computed and do not carry an
-                                    // invalid value.
-    double sum = 0;
-    typename DoFHandler<dim>::face_iterator face = cell->face(face_no);
-    for (unsigned int subface_no=0;
-        subface_no<GeometryInfo<dim>::subfaces_per_face;
-        ++subface_no) 
-      {
-       Assert (face_integrals.find(face->child(subface_no)) !=
-               face_integrals.end(),
-               ExcInternalError());
-       Assert (face_integrals[face->child(subface_no)] != -1e20,
-               ExcInternalError());
-      
-       sum += face_integrals[face->child(subface_no)];
-      };
-                                    // Finally store the value with
-                                    // the parent face.
-    face_integrals[face] = sum;
-  }
-  
-}
-
-
-                                // @sect3{A simulation framework}
-
-                                // In the previous example program,
-                                // we have had two functions that
-                                // were used to drive the process of
-                                // solving on subsequently finer
-                                // grids. We extend this here to
-                                // allow for a number of parameters
-                                // to be passed to these functions,
-                                // and put all of that into framework
-                                // class.
-                                //
-                                // You will have noted that this
-                                // program is built up of a number of
-                                // small parts (evaluation functions,
-                                // solver classes implementing
-                                // various refinement methods,
-                                // different dual functionals,
-                                // different problem and data
-                                // descriptions), which makes the
-                                // program relatively simple to
-                                // extend, but also allows to solve a
-                                // large number of different problems
-                                // by replacing one part by
-                                // another. We reflect this
-                                // flexibility by declaring a
-                                // structure in the following
-                                // framework class that holds a
-                                // number of parameters that may be
-                                // set to test various combinations
-                                // of the parts of this program, and
-                                // which can be used to test it at
-                                // various problems and
-                                // discretizations in a simple way.
-template <int dim>
-struct Framework
-{
-  public:
-                                    // First, we declare two
-                                    // abbreviations for simple use
-                                    // of the respective data types:
-    typedef Evaluation::EvaluationBase<dim> Evaluator;
-    typedef std::list<Evaluator*>           EvaluatorList;
-
-
-                                    // Then we have the structure
-                                    // which declares all the
-                                    // parameters that may be set. In
-                                    // the default constructor of the
-                                    // structure, these values are
-                                    // all set to default values, for
-                                    // simple use.
-    struct ProblemDescription 
-    {
-                                        // First allow for the
-                                        // degrees of the piecewise
-                                        // polynomials by which the
-                                        // primal and dual problems
-                                        // will be discretized. They
-                                        // default to (bi-,
-                                        // tri-)linear ansatz
-                                        // functions for the primal,
-                                        // and (bi-, tri-)quadratic
-                                        // ones for the dual
-                                        // problem. If a refinement
-                                        // criterion is chosen that
-                                        // does not need the solution
-                                        // of a dual problem, the
-                                        // value of the dual finite
-                                        // element degree is of
-                                        // course ignored.
-       unsigned int primal_fe_degree;
-       unsigned int dual_fe_degree;
-
-                                        // Then have an object that
-                                        // describes the problem
-                                        // type, i.e. right hand
-                                        // side, domain, boundary
-                                        // values, etc. The pointer
-                                        // needed here defaults to
-                                        // the Null pointer, i.e. you
-                                        // will have to set it in
-                                        // actual instances of this
-                                        // object to make it useful.
-       SmartPointer<const Data::SetUpBase<dim> > data;
-
-                                        // Since we allow to use
-                                        // different refinement
-                                        // criteria (global
-                                        // refinement, refinement by
-                                        // the Kelly error indicator,
-                                        // possibly with a weight,
-                                        // and using the dual
-                                        // estimator), define a
-                                        // number of enumeration
-                                        // values, and subsequently a
-                                        // variable of that type. It
-                                        // will default to
-                                        // ``dual_weighted_error_estimator''.
-       enum RefinementCriterion {
-             dual_weighted_error_estimator,
-             global_refinement,
-             kelly_indicator,
-             weighted_kelly_indicator
-       };
-
-       RefinementCriterion refinement_criterion;
-
-                                        // Next, an object that
-                                        // describes the dual
-                                        // functional. It is only
-                                        // needed if the dual
-                                        // weighted residual
-                                        // refinement is chosen, and
-                                        // also defaults to a Null
-                                        // pointer.
-       SmartPointer<const DualFunctional::DualFunctionalBase<dim> > dual_functional;
-
-                                        // Then a list of evaluation
-                                        // objects. Its default value
-                                        // is empty, i.e. no
-                                        // evaluation objects.
-       EvaluatorList evaluator_list;
-
-                                        // Next to last, a function
-                                        // that is used as a weight
-                                        // to the
-                                        // ``RefinementWeightedKelly''
-                                        // class. The default value
-                                        // of this pointer is zero,
-                                        // but you have to set it to
-                                        // some other value if you
-                                        // want to use the
-                                        // ``weighted_kelly_indicator''
-                                        // refinement criterion.
-       SmartPointer<const Function<dim> > kelly_weight;
-
-                                        // Finally, we have a
-                                        // variable that denotes the
-                                        // maximum number of degrees
-                                        // of freedom we allow for
-                                        // the (primal)
-                                        // discretization. If it is
-                                        // exceeded, we stop the
-                                        // process of solving and
-                                        // intermittend mesh
-                                        // refinement. Its default
-                                        // value is 20,000.
-       unsigned int max_degrees_of_freedom;
-
-                                        // Finally the default
-                                        // constructor of this class:
-       ProblemDescription ();
-    };
-
-                                    // The driver framework class
-                                    // only has one method which
-                                    // calls solver and mesh
-                                    // refinement intermittently, and
-                                    // does some other small tasks in
-                                    // between. Since it does not
-                                    // need data besides the
-                                    // parameters given to it, we
-                                    // make it static:
-    static void run (const ProblemDescription &descriptor);
-};
-
-
-                                // As for the implementation, first
-                                // the constructor of the parameter
-                                // object, setting all values to
-                                // their defaults:
-template <int dim>
-Framework<dim>::ProblemDescription::ProblemDescription ()
-               :
-               primal_fe_degree (1),
-               dual_fe_degree (2),
-               refinement_criterion (dual_weighted_error_estimator),
-               max_degrees_of_freedom (1000)
-{}
-
-
-
-                                // Then the function which drives the
-                                // whole process:
-template <int dim>
-void Framework<dim>::run (const ProblemDescription &descriptor)
-{
-                                  // First create a triangulation
-                                  // from the given data object,
-  Triangulation<dim>
-    triangulation (Triangulation<dim>::smoothing_on_refinement);
-  descriptor.data->create_coarse_grid (triangulation);
-
-                                  // then a set of finite elements
-                                  // and appropriate quadrature
-                                  // formula:
-  const FE_Q<dim>     primal_fe(descriptor.primal_fe_degree);
-  const FE_Q<dim>     dual_fe(descriptor.dual_fe_degree);
-  const QGauss<dim>   quadrature(descriptor.dual_fe_degree+1);
-  const QGauss<dim-1> face_quadrature(descriptor.dual_fe_degree+1);
-
-                                  // Next, select one of the classes
-                                  // implementing different
-                                  // refinement criteria.
-  LaplaceSolver::Base<dim> * solver = 0;
-  switch (descriptor.refinement_criterion)
-    {
-      case ProblemDescription::dual_weighted_error_estimator:
-      {
-       solver
-         = new LaplaceSolver::WeightedResidual<dim> (triangulation,
-                                                     primal_fe,
-                                                     dual_fe,
-                                                     quadrature,
-                                                     face_quadrature,
-                                                     descriptor.data->get_right_hand_side(),
-                                                     descriptor.data->get_boundary_values(),
-                                                     *descriptor.dual_functional);
-       break;
-      };
-       
-      case ProblemDescription::global_refinement:
-      {
-       solver
-         = new LaplaceSolver::RefinementGlobal<dim> (triangulation,
-                                                     primal_fe,
-                                                     quadrature,
-                                                     face_quadrature,
-                                                     descriptor.data->get_right_hand_side(),
-                                                     descriptor.data->get_boundary_values());
-       break;
-      };
-       
-      case ProblemDescription::kelly_indicator:
-      {
-       solver
-         = new LaplaceSolver::RefinementKelly<dim> (triangulation,
-                                                    primal_fe,
-                                                    quadrature,
-                                                    face_quadrature,
-                                                    descriptor.data->get_right_hand_side(),
-                                                    descriptor.data->get_boundary_values());
-       break;
-      };
-
-      case ProblemDescription::weighted_kelly_indicator:
-      {
-       solver
-         = new LaplaceSolver::RefinementWeightedKelly<dim> (triangulation,
-                                                            primal_fe,
-                                                            quadrature,
-                                                            face_quadrature,
-                                                            descriptor.data->get_right_hand_side(),
-                                                            descriptor.data->get_boundary_values(),
-                                                            *descriptor.kelly_weight);
-       break;
-      };
-           
-      default:
-           AssertThrow (false, ExcInternalError());
-    };
-  
-                                  // Now that all objects are in
-                                  // place, run the main loop. The
-                                  // stopping criterion is
-                                  // implemented at the bottom of the
-                                  // loop.
-                                  //
-                                  // In the loop, first set the new
-                                  // cycle number, then solve the
-                                  // problem, output its solution(s),
-                                  // apply the evaluation objects to
-                                  // it, then decide whether we want
-                                  // to refine the mesh further and
-                                  // solve again on this mesh, or
-                                  // jump out of the loop.
-  for (step=0; step<=n_steps; ++step)
-    {
-      std::cout << "Refinement cycle: "        << step
-               << std::endl;
-           
-      solver->set_refinement_cycle (step);
-      solver->solve_problem ();
-      solver->output_solution ();
-
-      std::cout << "   Number of degrees of freedom="
-               << solver->n_dofs() << std::endl;
-      
-      for (typename EvaluatorList::const_iterator
-            e = descriptor.evaluator_list.begin();
-          e != descriptor.evaluator_list.end(); ++e)
-       {
-         (*e)->set_refinement_cycle (step);
-         solver->postprocess (**e);
-       };
-
-           
-      if (solver->n_dofs() < descriptor.max_degrees_of_freedom)
-       solver->refine_grid ();
-      else
-       break;
-    };
-
-                                  // After the loop has run, clean up
-                                  // the screen, and delete objects
-                                  // no more needed:
-  std::cout << std::endl;
-  delete solver;
-  solver = 0;
-}
-
-
-
-
-                                // @sect3{The main function}
-
-                                // Here finally comes the main
-                                // function. It drives the whole
-                                // process by specifying a set of
-                                // parameters to be used for the
-                                // simulation (polynomial degrees,
-                                // evaluation and dual functionals,
-                                // etc), and passes them packed into
-                                // a structure to the frame work
-                                // class above.
-int main (int argc, char **argv) 
-{
-                                  // if no argument is given, then do 18
-                                  // iterations
-  if (argc == 1)
-    n_steps = 18;
-  else
-    if (argc == 2)
-      {
-       n_steps = atoi(argv[1]);
-       if ((n_steps==0) || (n_steps>100))
-         {
-           std::cout << "Please call this program with an argument in the range 1..100"
-                     << std::endl;
-           exit (1);
-         }
-      }
-    else
-      {
-       std::cout << "Please call this program with a single argument in the range 1..100"
-                 << std::endl;
-       exit (1);
-      }
-  
-       
-  
-  
-  deallog.depth_console (0);
-  try
-    {
-                                      // Describe the problem we want
-                                      // to solve here by passing a
-                                      // descriptor object to the
-                                      // function doing the rest of
-                                      // the work:
-      const unsigned int dim = 3;
-      Framework<dim>::ProblemDescription descriptor;
-
-                                      // First set the refinement
-                                      // criterion we wish to use:
-      descriptor.refinement_criterion
-       = Framework<dim>::ProblemDescription::dual_weighted_error_estimator;
-                                      // Here, we could as well have
-                                      // used ``global_refinement''
-                                      // or
-                                      // ``weighted_kelly_indicator''. Note
-                                      // that the information given
-                                      // about dual finite elements,
-                                      // dual functional, etc is only
-                                      // important for the given
-                                      // choice of refinement
-                                      // criterion, and is ignored
-                                      // otherwise.
-
-                                      // Then set the polynomial
-                                      // degrees of primal and dual
-                                      // problem. We choose here
-                                      // bi-linear and bi-quadratic
-                                      // ones:
-      descriptor.primal_fe_degree = 1;
-      descriptor.dual_fe_degree   = 2;
-
-                                      // Then set the description of
-                                      // the test case, i.e. domain,
-                                      // boundary values, and right
-                                      // hand side. These are
-                                      // prepackaged in classes. We
-                                      // take here the description of
-                                      // ``Exercise_2_3'', but you
-                                      // can also use
-                                      // ``CurvedRidges<dim>'':
-      descriptor.data = new Data::SetUp<Data::Exercise_2_3<dim>,dim> ();
-      
-                                      // Next set first a dual
-                                      // functional, then a list of
-                                      // evaluation objects. We
-                                      // choose as default the
-                                      // evaluation of the
-                                      // value at an
-                                      // evaluation point,
-                                      // represented by the classes
-                                      // ``PointValueEvaluation''
-                                      // in the namespaces of
-                                      // evaluation and dual
-                                      // functional classes. You can
-                                      // also set the
-                                      // ``PointXDerivativeEvaluation''
-                                      // classes for the x-derivative
-                                      // instead of the value
-                                      // at the evaluation point.
-                                      //
-                                      // Note that dual functional
-                                      // and evaluation objects
-                                      // should match. However, you
-                                      // can give as many evaluation
-                                      // functionals as you want, so
-                                      // you can have both point
-                                      // value and derivative
-                                      // evaluated after each step.
-                                      // One such additional
-                                      // evaluation is to output the
-                                      // grid in each step.
-      const Point<dim> evaluation_point (0., 0., 0.);
-      descriptor.dual_functional
-       = new DualFunctional::PointValueEvaluation<dim> (evaluation_point);
-      
-      Evaluation::PointValueEvaluation<dim>
-       postprocessor1 (evaluation_point);
-      Evaluation::GridOutput<dim>
-       postprocessor2 ("grid");
-      
-      descriptor.evaluator_list.push_back (&postprocessor1);
-      descriptor.evaluator_list.push_back (&postprocessor2);
-
-                                      // Set the maximal number of
-                                      // degrees of freedom after
-                                      // which we want the program to
-                                      // stop refining the mesh
-                                      // further:
-#if defined(SPEC_CPU)
-                                      // raise from 20000 to 30000. (jfk p6f)
-      descriptor.max_degrees_of_freedom = 30000;
-#else
-      descriptor.max_degrees_of_freedom = 20000;
-#endif
-      
-                                      // Finally pass the descriptor
-                                      // object to a function that
-                                      // runs the entire solution
-                                      // with it:
-      Framework<dim>::run (descriptor);
-    }
-
-                                  // Catch exceptions to give
-                                  // information about things that
-                                  // failed:
-  catch (std::exception &exc)
-    {
-      std::cerr << std::endl << std::endl
-               << "----------------------------------------------------"
-               << std::endl;
-      std::cerr << "Exception on processing: " << std::endl
-               << exc.what() << std::endl
-               << "Aborting!" << std::endl
-               << "----------------------------------------------------"
-               << std::endl;
-      return 1;
-    }
-  catch (...) 
-    {
-      std::cerr << std::endl << std::endl
-               << "----------------------------------------------------"
-               << std::endl;
-      std::cerr << "Unknown exception!" << std::endl
-               << "Aborting!" << std::endl
-               << "----------------------------------------------------"
-               << std::endl;
-      return 1;
-    };
-
-  return 0;
-}

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.