]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add the SPEC CPU2006 benchmark 447.dealII with 2 or 3 minor modifications to make...
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 25 Jan 2008 00:00:18 +0000 (00:00 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 25 Jan 2008 00:00:18 +0000 (00:00 +0000)
git-svn-id: https://svn.dealii.org/trunk@15681 0785d39b-7218-0410-832d-ea1e28bc413d

tests/benchmarks/spec2006-447.dealII.cc [new file with mode: 0644]

diff --git a/tests/benchmarks/spec2006-447.dealII.cc b/tests/benchmarks/spec2006-447.dealII.cc
new file mode 100644 (file)
index 0000000..617b60d
--- /dev/null
@@ -0,0 +1,4275 @@
+/* $Id: step-14.cc 2008 2006-01-23 23:53:16Z wolf $ */
+/* Author: Wolfgang Bangerth, ETH Zurich, 2002 */
+
+/*    $Id: step-14.cc 2008 2006-01-23 23:53:16Z wolf $       */
+/*    Version: $Name$                                          */
+/*                                                                */
+/*    Copyright (C) 2002, 2003, 2004, 2008 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 <base/quadrature_lib.h>
+#include <base/function.h>
+#include <base/logstream.h>
+#include <base/thread_management.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <lac/compressed_sparsity_pattern.h>
+#include <lac/sparse_matrix.h>
+#include <lac/solver_cg.h>
+#include <lac/vector_memory.h>
+#include <lac/precondition.h>
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/grid_out.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/grid_refinement.h>
+#include <grid/tria_boundary_lib.h>
+#include <grid/intergrid_map.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_constraints.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <fe/fe_q.h>
+#include <fe/fe_values.h>
+#include <fe/fe_tools.h>
+#include <fe/mapping_q.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+#include <numerics/error_estimator.h>
+#include <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<unsigned int> 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<unsigned int> 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,
+                                       QGauss3<dim-1>(),
+                                       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...
+    QGauss4<dim> quadrature;
+    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<unsigned int> 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,
+                                         QGauss3<dim-1>(),
+                                         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.