]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Replace the hp adaptation of step-12 by a mixture between step-13 and step-14. It...
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 13 Dec 2006 05:44:44 +0000 (05:44 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 13 Dec 2006 05:44:44 +0000 (05:44 +0000)
git-svn-id: https://svn.dealii.org/trunk@14229 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 37cc153d3255d33a42bdc55e2ca9401ea56d53e4..8cca3c9577b35a3199e0372279e863565ffd1068 100644 (file)
 /* $Id$ */
-/* Authors: Ralf Hartmann, University of Heidelberg, 2001 */
-/*          hp-Version, Oliver Kayser-Herold, TU-Braunschweig 2005 */
+/* Author: Wolfgang Bangerth, University of Heidelberg, 2001, 2002 */
 
 /*    $Id$       */
 /*    Version: $Name$                                          */
 /*                                                                */
-/*    Copyright (C) 2001-2006 by the deal.II authors              */
+/*    Copyright (C) 2001, 2002, 2003, 2004, 2006 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.                        */
 
-                                // The first few files have already
-                                // been covered in previous examples
-                                // and will thus not be further
-                                // commented on.
+
+                                // As in all programs, we start with
+                                // a list of include files from the
+                                // library, and as usual they are in
+                                // the standard order which is
+                                // <code>base</code> -- <code>lac</code> -- <code>grid</code> --
+                                // <code>dofs</code> -- <code>fe</code> -- <code>numerics</code>
+                                // (as each of these categories
+                                // roughly builds upon previous
+                                // ones), then C++ standard headers:
 #include <base/quadrature_lib.h>
 #include <base/function.h>
+#include <base/logstream.h>
+#include <base/table_handler.h>
+#include <base/thread_management.h>
 #include <lac/vector.h>
+#include <lac/full_matrix.h>
 #include <lac/sparse_matrix.h>
-#include <lac/vector_memory.h>
-#include <lac/solver_gmres.h>
+#include <lac/solver_cg.h>
 #include <lac/precondition.h>
-#include <lac/sparse_ilu.h>
 #include <grid/tria.h>
 #include <grid/grid_generator.h>
-#include <grid/grid_out.h>
-#include <grid/grid_refinement.h>
 #include <grid/tria_accessor.h>
 #include <grid/tria_iterator.h>
-#include <fe/fe_values.h>
-                                 // Instead of the usual DoFHandler
-                                 // the hp::DoFHandler class has to be
-                                 // used to gain access to the
-                                 // hp-Functionality. The hpDoFHandler
-                                 // provides essentially the same
-                                 // interface as the standard DoFHandler.
-#include <dofs/hp_dof_handler.h>
-                                 // Due to the implementation of the
-                                 // hp-Method, the DoFAccessor classes
-                                 // stay the same as for the DoFHandler,
-                                 // and hence also require the same
-                                 // include files.
+#include <grid/grid_refinement.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_constraints.h>
 #include <dofs/dof_accessor.h>
 #include <dofs/dof_tools.h>
-#include <dofs/dof_renumbering.h>
+#include <fe/fe_q.h>
+#include <fe/fe_values.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
 #include <numerics/data_out.h>
+#include <numerics/error_estimator.h>
 
-#include <fe/mapping_q1.h>
-#include <fe/fe_dgq.h>
-
-                                // We are going to use gradients as
-                                // refinement indicator.
-#include <numerics/derivative_approximation.h>
-                                // Finally we do some time comparison
-                                // using the <code>Timer</code> class.
-#include <base/timer.h>
-
-                                // And this again is C++:
+                                // Now for the C++ standard headers:
 #include <iostream>
 #include <fstream>
-
-                                // Now some additional utility classes will
-                                 // be necessary. As the name "collection"
-                                 // suggests, these are essentially
-                                 // container classes which store the
-                                 // quadrature and finite element objects,
-                                 // for the different polynomial degrees.
-#include <fe/fe_collection.h>
-#include <fe/q_collection.h>
-
-                                // The first of the following
-                                 // two files provides the hp::FEValues
-                                 // class, which implements the same
-                                 // functionality as the FEValues class
-                                 // with the difference that it takes
-                                 // "collection" objects instead of
-                                 // single finite element or quadrature
-                                 // objects. I.e. instead of the
-                                 // usual Quadrature object a
-                                 // QCollection object is needed.
-#include <fe/hp_fe_values.h>
-
-                                 // A compressed sparsity pattern is
-                                 // not an explicit prerequisite for the
-                                 // use of the hp-functionality. But as
-                                 // a standard sparsity pattern has to
-                                 // based on a bandwidth estimate for the
-                                 // highest polynomial degree, it would
-                                 // be prohibitively bad. Therefore,
-                                 // the recommended way is to explicitly
-                                 // build a compressed sparsity pattern before
-                                 // creating the matrices.
-#include <lac/compressed_sparsity_pattern.h>
-
+#include <list>
+#include <sstream>
 
                                 // The last step is as in all
                                 // previous programs:
 using namespace dealii;
 
-
-                                // @sect3{Equation data}
+                                // @sect3{Evaluation of the solution}
+
+                                // As for the program itself, we
+                                // first define classes that evaluate
+                                // the solutions of a Laplace
+                                // equation. In fact, they can
+                                // evaluate every kind of solution,
+                                // as long as it is described by a
+                                // <code>DoFHandler</code> object, and a
+                                // solution vector. We define them
+                                // here first, even before the
+                                // classes that actually generate the
+                                // solution to be evaluated, since we
+                                // need to declare an abstract base
+                                // class that the solver classes can
+                                // refer to.
                                 //
-                                // First we define the classes
-                                // representing the equation-specific
-                                // functions. Both classes, <code>RHS</code>
-                                // and <code>BoundaryValues</code>, are
-                                // derived from the <code>Function</code>
-                                // class. Only the <code>value_list</code>
-                                // function are implemented because
-                                // only lists of function values are
-                                // computed rather than single
-                                // values.
-template <int dim>
-class RHS:  public Function<dim>
+                                // From an abstract point of view, we
+                                // declare a pure base class
+                                // that provides an evaluation
+                                // operator <code>operator()</code> which will
+                                // do the evaluation of the solution
+                                // (whatever derived classes might
+                                // consider an <code>evaluation</code>). Since
+                                // this is the only real function of
+                                // this base class (except for some
+                                // bookkeeping machinery), one
+                                // usually terms such a class that
+                                // only has an <code>operator()</code> a
+                                // <code>functor</code> in C++ terminology,
+                                // since it is used just like a
+                                // function object.
+                                //
+                                // Objects of this functor type will
+                                // then later be passed to the solver
+                                // object, which applies it to the
+                                // solution just computed. The
+                                // evaluation objects may then
+                                // extract any quantity they like
+                                // from the solution. The advantage
+                                // of putting these evaluation
+                                // functions into a separate
+                                // hierarchy of classes is that by
+                                // design they cannot use the
+                                // internals of the solver object and
+                                // are therefore independent of
+                                // changes to the way the solver
+                                // works. Furthermore, it is trivial
+                                // to write another evaluation class
+                                // without modifying the solver
+                                // class, which speeds up programming
+                                // (not being able to use internals
+                                // of another class also means that
+                                // you do not have to worry about
+                                // them -- programming evaluators is
+                                // usually a rather quickly done
+                                // task), as well as compilation (if
+                                // solver and evaluation classes are
+                                // put into different files: the
+                                // solver only needs to see the
+                                // declaration of the abstract base
+                                // class, and therefore does not need
+                                // to be recompiled upon addition of
+                                // a new evaluation class, or
+                                // modification of an old one).
+                                // On a related note, you can reuse
+                                // the evaluation classes for other
+                                // projects, solving different
+                                // equations.
+                                //
+                                // In order to improve separation of
+                                // code into different modules, we
+                                // put the evaluation classes into a
+                                // namespace of their own. This makes
+                                // it easier to actually solve
+                                // different equations in the same
+                                // program, by assembling it from
+                                // existing building blocks. The
+                                // reason for this is that classes
+                                // for similar purposes tend to have
+                                // the same name, although they were
+                                // developed in different
+                                // contexts. In order to be able to
+                                // use them together in one program,
+                                // it is necessary that they are
+                                // placed in different
+                                // namespaces. This we do here:
+namespace Evaluation
 {
-  public:
-    virtual void value_list (const std::vector<Point<dim> > &points,
-                            std::vector<double> &values,
-                            const unsigned int component=0) const;
-};
 
+                                  // Now for the abstract base class
+                                  // of evaluation classes: its main
+                                  // purpose is to declare a pure
+                                  // virtual function <code>operator()</code>
+                                  // taking a <code>DoFHandler</code> object,
+                                  // and the solution vector. In
+                                  // order to be able to use pointers
+                                  // to this base class only, it also
+                                  // has to declare a virtual
+                                  // destructor, which however does
+                                  // nothing. Besides this, it only
+                                  // provides for a little bit of
+                                  // bookkeeping: since we usually
+                                  // want to evaluate solutions on
+                                  // subsequent refinement levels, we
+                                  // store the number of the present
+                                  // refinement cycle, and provide a
+                                  // function to change this number.
+  template <int dim>
+  class EvaluationBase 
+  {
+    public:
+      virtual ~EvaluationBase ();
 
-template <int dim>
-class BoundaryValues:  public Function<dim>
-{
-  public:
-    virtual void value_list (const std::vector<Point<dim> > &points,
-                            std::vector<double> &values,
-                            const unsigned int component=0) const;
-};
+      void set_refinement_cycle (const unsigned int refinement_cycle);
+      
+      virtual void operator () (const DoFHandler<dim> &dof_handler,
+                               const Vector<double>  &solution) const = 0;
+    protected:
+      unsigned int refinement_cycle;
+  };
+
+
+                                  // After the declaration has been
+                                  // discussed above, the
+                                  // implementation is rather
+                                  // straightforward:
+  template <int dim>
+  EvaluationBase<dim>::~EvaluationBase ()
+  {}
+  
 
+  
+  template <int dim>
+  void
+  EvaluationBase<dim>::set_refinement_cycle (const unsigned int step)
+  {
+    refinement_cycle = step;
+  }
 
-                                // The class <code>Beta</code> represents the
-                                // vector valued flow field of the
-                                // linear transport equation and is
-                                // not derived from the <code>Function</code>
-                                // class as we prefer to get function
-                                // values of type <code>Point</code> rather
-                                // than of type
-                                // <code>Vector@<double@></code>. This, because
-                                // there exist scalar products
-                                // between <code>Point</code> and <code>Point</code> as
-                                // well as between <code>Point</code> and
-                                // <code>Tensor</code>, simplifying terms like
-                                // $\beta\cdot n$ and
-                                // $\beta\cdot\nabla v$.
-                                 //
-                                 // An unnecessary empty constructor
-                                 // is added to the class to work
-                                 // around a bug in Compaq's cxx
-                                 // compiler which otherwise reports
-                                 // an error about an omitted
-                                 // initializer for an object of
-                                 // this class further down.
-template <int dim>
-class Beta
-{
-  public:
-    Beta () {};
-    void value_list (const std::vector<Point<dim> > &points,
-                    std::vector<Point<dim> > &values) const;
-};
 
+                                  // @sect4{%Point evaluation}
+
+                                  // The next thing is to implement
+                                  // actual evaluation classes. As
+                                  // noted in the introduction, we'd
+                                  // like to extract a point value
+                                  // from the solution, so the first
+                                  // class does this in its
+                                  // <code>operator()</code>. The actual point
+                                  // is given to this class through
+                                  // the constructor, as well as a
+                                  // table object into which it will
+                                  // put its findings.
+                                  //
+                                  // Finding out the value of a
+                                  // finite element field at an
+                                  // arbitrary point is rather
+                                  // difficult, if we cannot rely on
+                                  // knowing the actual finite
+                                  // element used, since then we
+                                  // cannot, for example, interpolate
+                                  // between nodes. For simplicity,
+                                  // we therefore assume here that
+                                  // the point at which we want to
+                                  // evaluate the field is actually a
+                                  // node. If, in the process of
+                                  // evaluating the solution, we find
+                                  // that we did not encounter this
+                                  // point upon looping over all
+                                  // vertices, we then have to throw
+                                  // an exception in order to signal
+                                  // to the calling functions that
+                                  // something has gone wrong, rather
+                                  // than silently ignore this error.
+                                  //
+                                  // In the step-9 example program,
+                                  // we have already seen how such an
+                                  // exception class can be declared,
+                                  // using the <code>DeclExceptionN</code>
+                                  // macros. We use this mechanism
+                                  // here again.
+                                  //
+                                  // From this, the actual
+                                  // declaration of this class should
+                                  // be evident. Note that of course
+                                  // even if we do not list a
+                                  // destructor explicitely, an
+                                  // implicit destructor is generated
+                                  // from the compiler, and it is
+                                  // virtual just as the one of the
+                                  // base class.
+  template <int dim>
+  class PointValueEvaluation : public EvaluationBase<dim>
+  {
+    public:
+      PointValueEvaluation (const Point<dim>   &evaluation_point,
+                           TableHandler       &results_table);
+      
+      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.");
+    private:
+      const Point<dim>  evaluation_point;
+      TableHandler     &results_table;
+  };
+
+
+                                  // As for the definition, the
+                                  // constructor is trivial, just
+                                  // taking data and storing it in
+                                  // object-local ones:
+  template <int dim>
+  PointValueEvaluation<dim>::
+  PointValueEvaluation (const Point<dim>   &evaluation_point,
+                       TableHandler       &results_table)
+                 :
+                 evaluation_point (evaluation_point),
+                 results_table (results_table)
+  {}
+  
 
-                                // The implementation of the
-                                // <code>value_list</code> functions of these
-                                // classes are rather simple.  For
-                                // simplicity the right hand side is
-                                // set to be zero but will be
-                                // assembled anyway.
-template <int dim>
-void RHS<dim>::value_list(const std::vector<Point<dim> > &points,
-                         std::vector<double> &values,
-                         const unsigned int) const
-{
-                                  // Usually we check whether input
-                                  // parameter have the right sizes.
-  Assert(values.size()==points.size(),
-        ExcDimensionMismatch(values.size(),points.size()));
 
-  for (unsigned int i=0; i<values.size(); ++i)
-    values[i]=0;
-}
+                                  // Now for the function that is
+                                  // mainly of interest in this
+                                  // class, the computation of the
+                                  // point value:
+  template <int dim>
+  void
+  PointValueEvaluation<dim>::
+  operator () (const DoFHandler<dim> &dof_handler,
+              const Vector<double>  &solution) const 
+  {
+                                    // First allocate a variable that
+                                    // will hold the point
+                                    // value. Initialize it with a
+                                    // value that is clearly bogus,
+                                    // so that if we fail to set it
+                                    // to a reasonable value, we will
+                                    // note at once. This may not be
+                                    // necessary in a function as
+                                    // small as this one, since we
+                                    // can easily see all possible
+                                    // paths of execution here, but
+                                    // it proved to be helpful for
+                                    // more complex cases, and so we
+                                    // employ this strategy here as
+                                    // well.
+    double point_value = 1e20;
+
+                                    // Then loop over all cells and
+                                    // all their vertices, and check
+                                    // whether a vertex matches the
+                                    // evaluation point. If this is
+                                    // the case, then extract the
+                                    // point value, set a flag that
+                                    // we have found the point of
+                                    // interest, and exit the loop.
+    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) == evaluation_point)
+         {
+                                            // In order to extract
+                                            // the point value from
+                                            // the global solution
+                                            // vector, pick that
+                                            // component that belongs
+                                            // to the vertex of
+                                            // interest, and, in case
+                                            // the solution is
+                                            // vector-valued, take
+                                            // the first component of
+                                            // it:
+           point_value = solution(cell->vertex_dof_index(vertex,0));
+                                            // Note that by this we
+                                            // have made an
+                                            // assumption that is not
+                                            // valid always and
+                                            // should be documented
+                                            // in the class
+                                            // declaration if this
+                                            // were code for a real
+                                            // application rather
+                                            // than a tutorial
+                                            // program: we assume
+                                            // that the finite
+                                            // element used for the
+                                            // solution we try to
+                                            // evaluate actually has
+                                            // degrees of freedom
+                                            // associated with
+                                            // vertices. This, for
+                                            // example, does not hold
+                                            // for discontinuous
+                                            // elements, were the
+                                            // support points for the
+                                            // shape functions
+                                            // happen to be located
+                                            // at the vertices, but
+                                            // are not associated
+                                            // with the vertices but
+                                            // rather with the cell
+                                            // interior, since
+                                            // association with
+                                            // vertices would imply
+                                            // continuity there. It
+                                            // would also not hold
+                                            // for edge oriented
+                                            // elements, and the
+                                            // like.
+                                            //
+                                            // Ideally, we would
+                                            // check this at the
+                                            // beginning of the
+                                            // function, for example
+                                            // by a statement like
+                                            // <code>Assert
+                                            // (dof_handler.get_fe().dofs_per_vertex
+                                            // @> 0,
+                                            // ExcNotImplemented())</code>,
+                                            // which should make it
+                                            // quite clear what is
+                                            // going wrong when the
+                                            // exception is
+                                            // triggered. In this
+                                            // case, we omit it
+                                            // (which is indeed bad
+                                            // style), but knowing
+                                            // that that does not
+                                            // hurt here, since the
+                                            // statement
+                                            // <code>cell-@>vertex_dof_index(vertex,0)</code>
+                                            // would fail if we asked
+                                            // it to give us the DoF
+                                            // index of a vertex if
+                                            // there were none.
+                                            //
+                                            // We stress again that
+                                            // this restriction on
+                                            // the allowed finite
+                                            // elements should be
+                                            // stated in the class
+                                            // documentation.
+
+                                            // Since we found the
+                                            // right point, we now
+                                            // set the respective
+                                            // flag and exit the
+                                            // innermost loop. The
+                                            // outer loop will the
+                                            // also be terminated due
+                                            // to the set flag.
+           evaluation_point_found = true;
+           break;
+         };
+
+                                    // Finally, we'd like to make
+                                    // sure that we have indeed found
+                                    // the evaluation point, since if
+                                    // that were not so we could not
+                                    // give a reasonable value of the
+                                    // solution there and the rest of
+                                    // the computations were useless
+                                    // anyway. So make sure through
+                                    // the <code>AssertThrow</code> macro
+                                    // already used in the step-9
+                                    // program that we have indeed
+                                    // found this point. If this is
+                                    // not so, the macro throws an
+                                    // exception of the type that is
+                                    // given to it as second
+                                    // argument, but compared to a
+                                    // straightforward <code>throw</code>
+                                    // statement, it fills the
+                                    // exception object with a set of
+                                    // additional information, for
+                                    // example the source file and
+                                    // line number where the
+                                    // exception was generated, and
+                                    // the condition that failed. If
+                                    // you have a <code>catch</code> clause in
+                                    // your main function (as this
+                                    // program has), you will catch
+                                    // all exceptions that are not
+                                    // caught somewhere in between
+                                    // and thus already handled, and
+                                    // this additional information
+                                    // will help you find out what
+                                    // happened and where it went
+                                    // wrong.
+    AssertThrow (evaluation_point_found,
+                ExcEvaluationPointNotFound(evaluation_point));
+                                    // Note that we have used the
+                                    // <code>Assert</code> macro in other
+                                    // example programs as well. It
+                                    // differed from the
+                                    // <code>AssertThrow</code> macro used
+                                    // here in that it simply aborts
+                                    // the program, rather than
+                                    // throwing an exception, and
+                                    // that it did so only in debug
+                                    // mode. It was the right macro
+                                    // to use to check about the size
+                                    // of vectors passed as arguments
+                                    // to functions, and the like.
+                                    //
+                                    // However, here the situation is
+                                    // different: whether we find the
+                                    // evaluation point or not may
+                                    // change from refinement to
+                                    // refinement (for example, if
+                                    // the four cells around point
+                                    // are coarsened away, then the
+                                    // point may vanish after
+                                    // refinement and
+                                    // coarsening). This is something
+                                    // that cannot be predicted from
+                                    // a few number of runs of the
+                                    // program in debug mode, but
+                                    // should be checked always, also
+                                    // in production runs. Thus the
+                                    // use of the <code>AssertThrow</code>
+                                    // macro here.
+    
+                                    // Now, if we are sure that we
+                                    // have found the evaluation
+                                    // point, we can add the results
+                                    // into the table of results:
+    results_table.add_value ("DoFs", dof_handler.n_dofs());
+    results_table.add_value ("u(x_0)", point_value);
+  }
 
 
-                                // The flow field is chosen to be
-                                // circular, counterclockwise, and with
-                                // the origin as midpoint.
-template <int dim>
-void Beta<dim>::value_list(const std::vector<Point<dim> > &points,
-                          std::vector<Point<dim> > &values) const
-{
-  Assert(values.size()==points.size(),
-        ExcDimensionMismatch(values.size(),points.size()));
 
-  for (unsigned int i=0; i<points.size(); ++i)
-    {
-      const Point<dim> &p=points[i];
-      Point<dim> &beta=values[i];
 
-      beta(0) = -p(1);
-      beta(1) = p(0);
-      beta /= std::sqrt(beta.square());
-    }
-}
+                                  // @sect4{Generating output}
+
+                                  // A different, maybe slightly odd
+                                  // kind of <code>evaluation</code> of a
+                                  // solution is to output it to a
+                                  // file in a graphical
+                                  // format. Since in the evaluation
+                                  // functions we are given a
+                                  // <code>DoFHandler</code> object and the
+                                  // solution vector, we have all we
+                                  // need to do this, so we can do it
+                                  // in an evaluation class. The
+                                  // reason for actually doing so
+                                  // instead of putting it into the
+                                  // class that computed the solution
+                                  // is that this way we have more
+                                  // flexibility: if we choose to
+                                  // only output certain aspects of
+                                  // it, or not output it at all. In
+                                  // any case, we do not need to
+                                  // modify the solver class, we just
+                                  // have to modify one of the
+                                  // modules out of which we build
+                                  // this program. This form of
+                                  // encapsulation, as above, helps
+                                  // us to keep each part of the
+                                  // program rather simple as the
+                                  // interfaces are kept simple, and
+                                  // no access to hidden data is
+                                  // possible.
+                                  //
+                                  // Since this class which generates
+                                  // the output is derived from the
+                                  // common <code>EvaluationBase</code> base
+                                  // class, its main interface is the
+                                  // <code>operator()</code>
+                                  // function. Furthermore, it has a
+                                  // constructor taking a string that
+                                  // will be used as the base part of
+                                  // the file name to which output
+                                  // will be sent (we will augment it
+                                  // by a number indicating the
+                                  // number of the refinement cycle
+                                  // -- the base class has this
+                                  // information at hand --, and a
+                                  // suffix), and the constructor
+                                  // also takes a value that
+                                  // indicates which format is
+                                  // requested, i.e. for which
+                                  // graphics program we shall
+                                  // generate output (from this we
+                                  // will then also generate the
+                                  // suffix of the filename to which
+                                  // we write).
+                                  //
+                                  // Regarding the output format, the
+                                  // <code>DataOutInterface</code> class
+                                  // (which is a base class of
+                                  // <code>DataOut</code> through which we
+                                  // will access its fields) provides
+                                  // an enumeration field
+                                  // <code>OutputFormat</code>, which lists
+                                  // names for all supported output
+                                  // formats. At the time of writing
+                                  // of this program, the supported
+                                  // graphics formats are represented
+                                  // by the enum values <code>ucd</code>,
+                                  // <code>gnuplot</code>, <code>povray</code>,
+                                  // <code>eps</code>, <code>gmv</code>, <code>tecplot</code>,
+                                  // <code>tecplot_binary</code>, <code>dx</code>, and
+                                  // <code>vtk</code>, but this list will
+                                  // certainly grow over time. Now,
+                                  // within various functions of that
+                                  // base class, you can use values
+                                  // of this type to get information
+                                  // about these graphics formats
+                                  // (for example the default suffix
+                                  // used for files of each format),
+                                  // and you can call a generic
+                                  // <code>write</code> function, which then
+                                  // branches to the
+                                  // <code>write_gnuplot</code>,
+                                  // <code>write_ucd</code>, etc functions
+                                  // which we have used in previous
+                                  // examples already, based on the
+                                  // value of a second argument given
+                                  // to it denoting the required
+                                  // output format. This mechanism
+                                  // makes it simple to write an
+                                  // extensible program that can
+                                  // decide which output format to
+                                  // use at runtime, and it also
+                                  // makes it rather simple to write
+                                  // the program in a way such that
+                                  // it takes advantage of newly
+                                  // implemented output formats,
+                                  // without the need to change the
+                                  // application program.
+                                  //
+                                  // Of these two fields, the base
+                                  // name and the output format
+                                  // descriptor, the constructor
+                                  // takes values and stores them for
+                                  // later use by the actual
+                                  // evaluation function.
+  template <int dim>
+  class SolutionOutput : public EvaluationBase<dim>
+  {
+    public:
+      SolutionOutput (const std::string                         &output_name_base,
+                     const typename DataOut<dim>::OutputFormat  output_format);
+      
+      virtual void operator () (const DoFHandler<dim> &dof_handler,
+                               const Vector<double>  &solution) const;
+    private:
+      const std::string                         output_name_base;
+      const typename DataOut<dim>::OutputFormat output_format;
+  };
+
+
+  template <int dim>
+  SolutionOutput<dim>::
+  SolutionOutput (const std::string                         &output_name_base,
+                 const typename DataOut<dim>::OutputFormat  output_format)
+                 :
+                 output_name_base (output_name_base),
+                 output_format (output_format)
+  {}
+  
+
+                                  // After the description above, the
+                                  // function generating the actual
+                                  // output is now relatively
+                                  // straightforward. The only
+                                  // particularly interesting feature
+                                  // over previous example programs
+                                  // is the use of the
+                                  // <code>DataOut::default_suffix</code>
+                                  // function, returning the usual
+                                  // suffix for files of a given
+                                  // format (e.g. ".eps" for
+                                  // encapsulated postscript files,
+                                  // ".gnuplot" for Gnuplot files),
+                                  // and of the generic
+                                  // <code>DataOut::write</code> function with
+                                  // a second argument, which
+                                  // branches to the actual output
+                                  // functions for the different
+                                  // graphics formats, based on the
+                                  // value of the format descriptor
+                                  // passed as second argument.
+                                  //
+                                  // Also note that we have to prefix
+                                  // <code>this-@></code> to access a member
+                                  // variable of the template
+                                  // dependent base class. The reason
+                                  // here, and further down in the
+                                  // program is the same as the one
+                                  // described in the step-7 example
+                                  // program (look for <code>two-stage
+                                  // name lookup</code> there).
+  template <int dim>
+  void
+  SolutionOutput<dim>::operator () (const DoFHandler<dim> &dof_handler,
+                                   const Vector<double>  &solution) const
+  {
+    DataOut<dim> data_out;
+    data_out.attach_dof_handler (dof_handler);
+    data_out.add_data_vector (solution, "solution");
+    data_out.build_patches ();
+  
+    std::ostringstream filename;
+    filename << output_name_base << "-"
+            << this->refinement_cycle
+            << data_out.default_suffix (output_format)
+            << std::ends;
+    std::ofstream out (filename.str().c_str());
+    
+    data_out.write (out, output_format);
+  }
 
 
-                                // Hence the inflow boundary of the
-                                // unit square [0,1]^2 are the right
-                                // and the lower boundaries. We
-                                // prescribe discontinuous boundary
-                                // values 1 and 0 on the x-axis and
-                                // value 0 on the right boundary. The
-                                // values of this function on the
-                                // outflow boundaries will not be
-                                // used within the DG scheme.
-template <int dim>
-void BoundaryValues<dim>::value_list(const std::vector<Point<dim> > &points,
-                                      std::vector<double> &values,
-                                      const unsigned int) const
-{
-  Assert(values.size()==points.size(),
-        ExcDimensionMismatch(values.size(),points.size()));
 
-  for (unsigned int i=0; i<values.size(); ++i)
-    {
-      if (points[i](0)<0.5)
-       values[i]=1.;
-      else
-       values[i]=0.;
-    }
+                                  // @sect4{Other evaluations}
+  
+                                  // In practical applications, one
+                                  // would add here a list of other
+                                  // possible evaluation classes,
+                                  // representing quantities that one
+                                  // may be interested in. For this
+                                  // example, that much shall be
+                                  // sufficient, so we close the
+                                  // namespace.
 }
 
-
-                                // @sect3{Class: DGTransportEquation}
+  
+                                // @sect3{The Laplace solver classes}
+
+                                // After defining what we want to
+                                // know of the solution, we should
+                                // now care how to get at it. We will
+                                // pack everything we need into a
+                                // namespace of its own, for much the
+                                // same reasons as for the
+                                // evaluations above.
                                 //
-                                // Next we define the
-                                // equation-dependent and
-                                // DG-method-dependent class
-                                // <code>DGTransportEquation</code>. Its
-                                // member functions were already
-                                // mentioned in the Introduction and
-                                // will be explained
-                                // below. Furthermore it includes
-                                // objects of the previously defined
-                                // <code>Beta</code>, <code>RHS</code> and
-                                // <code>BoundaryValues</code> function
-                                // classes.
-template <int dim>
-class DGTransportEquation
+                                // Since we have discussed Laplace
+                                // solvers already in considerable
+                                // detail in previous examples, there
+                                // is not much new stuff
+                                // following. Rather, we have to a
+                                // great extent cannibalized previous
+                                // examples and put them, in slightly
+                                // different form, into this example
+                                // program. We will therefore mostly
+                                // be concerned with discussing the
+                                // differences to previous examples.
+                                //
+                                // Basically, as already said in the
+                                // introduction, the lack of new
+                                // stuff in this example is
+                                // deliberate, as it is more to
+                                // demonstrate software design
+                                // practices, rather than
+                                // mathematics. The emphasis in
+                                // explanations below will therefore
+                                // be more on the actual
+                                // implementation.
+namespace LaplaceSolver
 {
-  public:
-    DGTransportEquation();
-
-    void assemble_cell_term(const FEValues<dim>& fe_v,
-                           FullMatrix<double> &u_v_matrix,
-                           Vector<double> &cell_vector) const;
-    
-    void assemble_boundary_term(const FEFaceValues<dim>& fe_v,
-                               FullMatrix<double> &u_v_matrix,
-                               Vector<double> &cell_vector) const;
-    
-    void assemble_face_term1(const FEFaceValuesBase<dim>& fe_v,
-                            const FEFaceValuesBase<dim>& fe_v_neighbor,
-                            FullMatrix<double> &u_v_matrix,
-                            FullMatrix<double> &un_v_matrix) const;
-
-    void assemble_face_term2(const FEFaceValuesBase<dim>& fe_v,
-                            const FEFaceValuesBase<dim>& fe_v_neighbor,
-                            FullMatrix<double> &u_v_matrix,
-                            FullMatrix<double> &un_v_matrix,
-                            FullMatrix<double> &u_vn_matrix,
-                            FullMatrix<double> &un_vn_matrix) const;
-  private:
-    const Beta<dim> beta_function;
-    const RHS<dim> rhs_function;
-    const BoundaryValues<dim> boundary_function;
-};
+                                  // @sect4{An abstract base class}
+
+                                  // In defining a Laplace solver, we
+                                  // start out by declaring an
+                                  // abstract base class, that has no
+                                  // functionality itself except for
+                                  // taking and storing a pointer to
+                                  // the triangulation to be used
+                                  // later.
+                                  //
+                                  // This base class is very general,
+                                  // and could as well be used for
+                                  // any other stationary problem. It
+                                  // provides declarations of
+                                  // functions that shall, in derived
+                                  // classes, solve a problem,
+                                  // postprocess the solution with a
+                                  // list of evaluation objects, and
+                                  // refine the grid,
+                                  // respectively. None of these
+                                  // functions actually does
+                                  // something itself in the base
+                                  // class.
+                                  //
+                                  // Due to the lack of actual
+                                  // functionality, the programming
+                                  // style of declaring very abstract
+                                  // base classes reminds of the
+                                  // style used in Smalltalk or Java
+                                  // programs, where all classes are
+                                  // derived from entirely abstract
+                                  // classes <code>Object</code>, even number
+                                  // representations. The author
+                                  // admits that he does not
+                                  // particularly like the use of
+                                  // such a style in C++, as it puts
+                                  // style over reason. Furthermore,
+                                  // it promotes the use of virtual
+                                  // functions for everything (for
+                                  // example, in Java, all functions
+                                  // are virtual per se), which,
+                                  // however, has proven to be rather
+                                  // inefficient in many applications
+                                  // where functions are often only
+                                  // accessing data, not doing
+                                  // computations, and therefore
+                                  // quickly return; the overhead of
+                                  // virtual functions can then be
+                                  // significant. The opinion of the
+                                  // author is to have abstract base
+                                  // classes wherever at least some
+                                  // part of the code of actual
+                                  // implementations can be shared
+                                  // and thus separated into the base
+                                  // class.
+                                  //
+                                  // Besides all these theoretical
+                                  // questions, we here have a good
+                                  // reason, which will become
+                                  // clearer to the reader
+                                  // below. Basically, we want to be
+                                  // able to have a family of
+                                  // different Laplace solvers that
+                                  // differ so much that no larger
+                                  // common subset of functionality
+                                  // could be found. We therefore
+                                  // just declare such an abstract
+                                  // base class, taking a pointer to
+                                  // a triangulation in the
+                                  // constructor and storing it
+                                  // henceforth. Since this
+                                  // triangulation will be used
+                                  // throughout all computations, we
+                                  // have to make sure that the
+                                  // triangulation exists until the
+                                  // destructor exits. We do this by
+                                  // keeping a <code>SmartPointer</code> to
+                                  // this triangulation, which uses a
+                                  // counter in the triangulation
+                                  // class to denote the fact that
+                                  // there is still an object out
+                                  // there using this triangulation,
+                                  // thus leading to an abort in case
+                                  // the triangulation is attempted
+                                  // to be destructed while this
+                                  // object still uses it.
+                                  //
+                                  // Note that while the pointer
+                                  // itself is declared constant
+                                  // (i.e. throughout the lifetime of
+                                  // this object, the pointer points
+                                  // to the same object), it is not
+                                  // declared as a pointer to a
+                                  // constant triangulation. In fact,
+                                  // by this we allow that derived
+                                  // classes refine or coarsen the
+                                  // triangulation within the
+                                  // <code>refine_grid</code> function.
+                                  //
+                                  // Finally, we have a function
+                                  // <code>n_dofs</code> is only a tool for
+                                  // the driver functions to decide
+                                  // whether we want to go on with
+                                  // mesh refinement or not. It
+                                  // returns the number of degrees of
+                                  // freedom the present simulation
+                                  // has.
+  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;
+      
+    protected:
+      const SmartPointer<Triangulation<dim> > triangulation;
+  };
 
 
-template <int dim>
-DGTransportEquation<dim>::DGTransportEquation ()
-               :
-               beta_function (),
-               rhs_function (),
-               boundary_function ()
-{}
+                                  // The implementation of the only
+                                  // two non-abstract functions is
+                                  // then rather boring:
+  template <int dim>
+  Base<dim>::Base (Triangulation<dim> &coarse_grid)
+                 :
+                 triangulation (&coarse_grid)
+  {}
 
 
-                                // @sect4{Function: assemble_cell_term}
-                                //
-                                // The <code>assemble_cell_term</code>
-                                // function assembles the cell terms
-                                // of the discretization.
-                                // <code>u_v_matrix</code> is a cell matrix,
-                                // i.e. for a DG method of degree 1,
-                                // it is of size 4 times 4, and
-                                // <code>cell_vector</code> is of size 4.
-                                // When this function is invoked,
-                                // <code>fe_v</code> is already reinit'ed with the
-                                // current cell before and includes
-                                // all shape values needed.
-template <int dim>
-void DGTransportEquation<dim>::assemble_cell_term(
-  const FEValues<dim> &fe_v,
-  FullMatrix<double> &u_v_matrix,
-  Vector<double> &cell_vector) const
-{
-                                  // First we ask <code>fe_v</code> for the
-                                  // shape gradients, shape values and
-                                  // quadrature weights,
-  const std::vector<double> &JxW = fe_v.get_JxW_values ();
-
-                                  // Then the flow field beta and the
-                                  // <code>rhs_function</code> are evaluated at
-                                  // the quadrature points,
-  std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
-  std::vector<double> rhs (fe_v.n_quadrature_points);
-  
-  beta_function.value_list (fe_v.get_quadrature_points(), beta);
-  rhs_function.value_list (fe_v.get_quadrature_points(), rhs);
+  template <int dim>
+  Base<dim>::~Base () 
+  {}
   
-                                  // and the cell matrix and cell
-                                  // vector are assembled due to the
-                                  // terms $-(u,\beta\cdot\nabla
-                                  // v)_K$ and $(f,v)_K$.
-  for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
-    for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i) 
+
+                                  // @sect4{A general solver class}
+
+                                  // Following now the main class
+                                  // that implements assembling the
+                                  // matrix of the linear system,
+                                  // solving it, and calling the
+                                  // postprocessor objects on the
+                                  // solution. It implements the
+                                  // <code>solve_problem</code> and
+                                  // <code>postprocess</code> functions
+                                  // declared in the base class. It
+                                  // does not, however, implement the
+                                  // <code>refine_grid</code> method, as mesh
+                                  // refinement will be implemented
+                                  // in a number of derived classes.
+                                  //
+                                  // It also declares a new abstract
+                                  // virtual function,
+                                  // <code>assemble_rhs</code>, that needs to
+                                  // be overloaded in subclasses. The
+                                  // reason is that we will implement
+                                  // two different classes that will
+                                  // implement different methods to
+                                  // assemble the right hand side
+                                  // vector. This function might also
+                                  // be interesting in cases where
+                                  // the right hand side depends not
+                                  // simply on a continuous function,
+                                  // but on something else as well,
+                                  // for example the solution of
+                                  // another discretized problem,
+                                  // etc. The latter happens
+                                  // frequently in non-linear
+                                  // problems.
+                                  //
+                                  // As we mentioned previously, the
+                                  // actual content of this class is
+                                  // not new, but a mixture of
+                                  // various techniques already used
+                                  // in previous examples. We will
+                                  // therefore not discuss them in
+                                  // detail, but refer the reader to
+                                  // these programs.
+                                  //
+                                  // Basically, in a few words, the
+                                  // constructor of this class takes
+                                  // pointers to a triangulation, a
+                                  // finite element, and a function
+                                  // object representing the boundary
+                                  // values. These are either passed
+                                  // down to the base class's
+                                  // constructor, or are stored and
+                                  // used to generate a
+                                  // <code>DoFHandler</code> object
+                                  // later. Since finite elements and
+                                  // quadrature formula should match,
+                                  // it is also passed a quadrature
+                                  // object.
+                                  //
+                                  // The <code>solve_problem</code> sets up
+                                  // the data structures for the
+                                  // actual solution, calls the
+                                  // functions to assemble the linear
+                                  // system, and solves it.
+                                  //
+                                  // The <code>postprocess</code> function
+                                  // finally takes an evaluation
+                                  // object and applies it to the
+                                  // computed solution.
+                                  //
+                                  // The <code>n_dofs</code> function finally
+                                  // implements the pure virtual
+                                  // function of the base class.
+  template <int dim>
+  class Solver : public virtual Base<dim>
+  {
+    public:
+      Solver (Triangulation<dim>       &triangulation,
+             const FiniteElement<dim> &fe,
+             const Quadrature<dim>    &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;
+      
+                                      // In the protected section of
+                                      // this class, we first have a
+                                      // number of member variables,
+                                      // of which the use should be
+                                      // clear from the previous
+                                      // examples:
+    protected:
+      const SmartPointer<const FiniteElement<dim> >  fe;
+      const SmartPointer<const Quadrature<dim> >     quadrature;
+      DoFHandler<dim>                                dof_handler;
+      Vector<double>                                 solution;
+      const SmartPointer<const Function<dim> >       boundary_values;
+
+                                      // Then we declare an abstract
+                                      // function that will be used
+                                      // to assemble the right hand
+                                      // side. As explained above,
+                                      // there are various cases for
+                                      // which this action differs
+                                      // strongly in what is
+                                      // necessary, so we defer this
+                                      // to derived classes:
+      virtual void assemble_rhs (Vector<double> &rhs) const = 0;
+    
+                                      // Next, in the private
+                                      // section, we have a small
+                                      // class which represents an
+                                      // entire linear system, i.e. a
+                                      // matrix, a right hand side,
+                                      // and a solution vector, as
+                                      // well as the constraints that
+                                      // are applied to it, such as
+                                      // those due to hanging
+                                      // nodes. Its constructor
+                                      // initializes the various
+                                      // subobjects, and there is a
+                                      // function that implements a
+                                      // conjugate gradient method as
+                                      // solver.
+    private:
+      struct LinearSystem
       {
-       for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
-         u_v_matrix(i,j) -= beta[point]*fe_v.shape_grad(i,point)*
-                             fe_v.shape_value(j,point) *
-                             JxW[point];
+         LinearSystem (const DoFHandler<dim> &dof_handler);
+
+         void solve (Vector<double> &solution) const;
        
-       cell_vector(i) += rhs[point] * fe_v.shape_value(i,point) * JxW[point];
-      }
-}
+         ConstraintMatrix     hanging_node_constraints;
+         SparsityPattern      sparsity_pattern;
+         SparseMatrix<double> matrix;
+         Vector<double>       rhs;
+      };
+
+                                      // Finally, there is a pair of
+                                      // functions which will be used
+                                      // to assemble the actual
+                                      // system matrix. It calls the
+                                      // virtual function assembling
+                                      // the right hand side, and
+                                      // installs a number threads
+                                      // each running the second
+                                      // function which assembles
+                                      // part of the system
+                                      // matrix. The mechanism for
+                                      // doing so is the same as in
+                                      // the step-9 example program.
+      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;
+  };
+
+
+
+                                  // Now here comes the constructor
+                                  // of the class. It does not do
+                                  // much except store pointers to
+                                  // the objects given, and generate
+                                  // <code>DoFHandler</code> object
+                                  // initialized with the given
+                                  // pointer to a triangulation. This
+                                  // causes the DoF handler to store
+                                  // that pointer, but does not
+                                  // already generate a finite
+                                  // element numbering (we only ask
+                                  // for that in the
+                                  // <code>solve_problem</code> function).
+  template <int dim>
+  Solver<dim>::Solver (Triangulation<dim>       &triangulation,
+                      const FiniteElement<dim> &fe,
+                      const Quadrature<dim>    &quadrature,
+                      const Function<dim>      &boundary_values)
+                 :
+                 Base<dim> (triangulation),
+                 fe (&fe),
+                  quadrature (&quadrature),
+                 dof_handler (triangulation),
+                 boundary_values (&boundary_values)
+  {}
+
+
+                                  // The destructor is simple, it
+                                  // only clears the information
+                                  // stored in the DoF handler object
+                                  // to release the memory.
+  template <int dim>
+  Solver<dim>::~Solver () 
+  {
+    dof_handler.clear ();
+  }
 
 
-                                // @sect4{Function: assemble_boundary_term}
-                                //
-                                // The <code>assemble_boundary_term</code>
-                                // function assembles the face terms
-                                // at boundary faces.  When this
-                                // function is invoked, <code>fe_v</code> is
-                                // already reinit'ed with the current
-                                // cell and current face. Hence it
-                                // provides the shape values on that
-                                // boundary face.
-template <int dim>
-void DGTransportEquation<dim>::assemble_boundary_term(
-  const FEFaceValues<dim>& fe_v,    
-  FullMatrix<double> &u_v_matrix,
-  Vector<double> &cell_vector) const
-{
-                                  // Again, as in the previous
-                                  // function, we ask the <code>FEValues</code>
-                                  // object for the shape values and
-                                  // the quadrature weights
-  const std::vector<double> &JxW = fe_v.get_JxW_values ();
-                                  // but here also for the normals.
-  const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
-
-                                  // We evaluate the flow field
-                                  // and the boundary values at the
-                                  // quadrature points.
-  std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
-  std::vector<double> g(fe_v.n_quadrature_points);
-  
-  beta_function.value_list (fe_v.get_quadrature_points(), beta);
-  boundary_function.value_list (fe_v.get_quadrature_points(), g);
-
-                                  // Then we assemble cell vector and
-                                  // cell matrix according to the DG
-                                  // method given in the
-                                  // introduction.
-  for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
-    {
-      const double beta_n=beta[point] * normals[point];      
-                                        // We assemble the term
-                                        // $(\beta\cdot n
-                                        // u,v)_{\partial K_+}$,
-      if (beta_n>0)
-       for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
-         for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
-           u_v_matrix(i,j) += beta_n *
-                              fe_v.shape_value(j,point) *
-                              fe_v.shape_value(i,point) *
-                              JxW[point];
-      else
-                                        // and the term $(\beta\cdot
-                                        // n g,v)_{\partial
-                                        // K_-\cap\partial\Omega}$,
-       for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
-         cell_vector(i) -= beta_n *
-                           g[point] *
-                           fe_v.shape_value(i,point) *
-                           JxW[point];
-    }
-}
+                                  // The next function is the one
+                                  // which delegates the main work in
+                                  // solving the problem: it sets up
+                                  // the DoF handler object with the
+                                  // finite element given to the
+                                  // constructor of this object, the
+                                  // creates an object that denotes
+                                  // the linear system (i.e. the
+                                  // matrix, the right hand side
+                                  // vector, and the solution
+                                  // vector), calls the function to
+                                  // assemble it, and finally solves
+                                  // it:
+  template <int dim>
+  void
+  Solver<dim>::solve_problem ()
+  {
+    dof_handler.distribute_dofs (*fe);
 
+    std::cout << "Number of degrees of freedom: "
+              << dof_handler.n_dofs()
+              << std::endl;
+    
+    solution.reinit (dof_handler.n_dofs());
 
-                                // @sect4{Function: assemble_face_term1}
-                                //
-                                // The <code>assemble_face_term1</code>
-                                // function assembles the face terms
-                                // corresponding to the first version
-                                // of the DG method, cf. above. For
-                                // that case, the face terms are
-                                // given as a sum of integrals over
-                                // all cell boundaries.
-                                //
-                                // When this function is invoked,
-                                // <code>fe_v</code> and <code>fe_v_neighbor</code> are
-                                // already reinit'ed with the current
-                                // cell and the neighoring cell,
-                                // respectively, as well as with the
-                                // current face. Hence they provide
-                                // the inner and outer shape values
-                                // on the face.
-                                //
-                                // In addition to the cell matrix
-                                // <code>u_v_matrix</code> this function has
-                                // got a new argument
-                                // <code>un_v_matrix</code>, that stores
-                                // contributions to the system matrix
-                                // that are based on outer values of
-                                // u, see $\hat u_h$ in the
-                                // introduction, and inner values of
-                                // v, see $v_h$. Here we note that
-                                // <code>un</code> is the short notation for
-                                // <code>u_neighbor</code> and represents
-                                // $\hat u_h$.
-template <int dim>
-void DGTransportEquation<dim>::assemble_face_term1(
-  const FEFaceValuesBase<dim>& fe_v,
-  const FEFaceValuesBase<dim>& fe_v_neighbor,      
-  FullMatrix<double> &u_v_matrix,
-  FullMatrix<double> &un_v_matrix) const
-{
-                                  // Again, as in the previous
-                                  // function, we ask the FEValues
-                                  // objects for the shape values,
-                                  // the quadrature weights and the
-                                  // normals
-  const std::vector<double> &JxW = fe_v.get_JxW_values ();
-  const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
-
-                                  // and we evaluate the flow field
-                                  // at the quadrature points.
-  std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
-  beta_function.value_list (fe_v.get_quadrature_points(), beta);
-
-                                  // Then we assemble the cell
-                                  // matrices according to the DG
-                                  // method given in the
-                                  // introduction.
-  for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
-    {
-      const double beta_n=beta[point] * normals[point];
-                                        // We assemble the term
-                                        // $(\beta\cdot n
-                                        // u,v)_{\partial K_+}$,
-      if (beta_n>0)
-       for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
-         for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
-           u_v_matrix(i,j) += beta_n *
-                              fe_v.shape_value(j,point) *
-                              fe_v.shape_value(i,point) *
-                              JxW[point];
-      else
-                                        // and the
-                                        // term $(\beta\cdot n
-                                        // \hat u,v)_{\partial
-                                        // K_-}$.
-       for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
-         for (unsigned int k=0; k<fe_v_neighbor.dofs_per_cell; ++k)
-           un_v_matrix(i,k) += beta_n *
-                               fe_v_neighbor.shape_value(k,point) *
-                               fe_v.shape_value(i,point) *
-                               JxW[point];
-    }
-}
+    LinearSystem linear_system (dof_handler);
 
+    std::cout << "Number of constraints       : "
+              << linear_system.hanging_node_constraints.n_constraints()
+              << std::endl;
 
-                                // @sect4{Function: assemble_face_term2}
-                                //
-                                // Now we look at the
-                                // <code>assemble_face_term2</code> function
-                                // that assembles the face terms
-                                // corresponding to the second
-                                // version of the DG method,
-                                // cf. above. For that case the face
-                                // terms are given as a sum of
-                                // integrals over all faces.  Here we
-                                // need two additional cell matrices
-                                // <code>u_vn_matrix</code> and
-                                // <code>un_vn_matrix</code> that will store
-                                // contributions due to terms
-                                // involving u and vn as well as un
-                                // and vn.
-template <int dim>
-void DGTransportEquation<dim>::assemble_face_term2(
-  const FEFaceValuesBase<dim>& fe_v,
-  const FEFaceValuesBase<dim>& fe_v_neighbor,
-  FullMatrix<double> &u_v_matrix,
-  FullMatrix<double> &un_v_matrix,
-  FullMatrix<double> &u_vn_matrix,
-  FullMatrix<double> &un_vn_matrix) const
-{
-                                  // the first few lines are the same
-  const std::vector<double> &JxW = fe_v.get_JxW_values ();
-  const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
+    assemble_linear_system (linear_system);
+    linear_system.solve (solution);
+  }
 
-  std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
-  
-  beta_function.value_list (fe_v.get_quadrature_points(), beta);
 
-  for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
-    {
-      const double beta_n=beta[point] * normals[point];
-      if (beta_n>0)
-       {
-                                          // This terms we've already seen.
-         for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
-           for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
-             u_v_matrix(i,j) += beta_n *
-                                fe_v.shape_value(j,point) *
-                                fe_v.shape_value(i,point) *
-                                JxW[point];
-
-                                          // We additionally assemble
-                                          // the term $(\beta\cdot n
-                                          // u,\hat v)_{\partial
-                                          // K_+}$,
-         for (unsigned int k=0; k<fe_v_neighbor.dofs_per_cell; ++k)
-           for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
-             u_vn_matrix(k,j) -= beta_n *
-                                 fe_v.shape_value(j,point) *
-                                 fe_v_neighbor.shape_value(k,point) *
-                                 JxW[point];
-       }
-      else
-       {
-                                          // This one we've already
-                                          // seen, too.
-         for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
-           for (unsigned int l=0; l<fe_v_neighbor.dofs_per_cell; ++l)
-             un_v_matrix(i,l) += beta_n *
-                                 fe_v_neighbor.shape_value(l,point) *
-                                 fe_v.shape_value(i,point) *
-                                 JxW[point];
-
-                                          // And this is another new
-                                          // one: $(\beta\cdot n \hat
-                                          // u,\hat v)_{\partial
-                                          // K_-}$.
-         for (unsigned int k=0; k<fe_v_neighbor.dofs_per_cell; ++k)
-           for (unsigned int l=0; l<fe_v_neighbor.dofs_per_cell; ++l)
-             un_vn_matrix(k,l) -= beta_n *
-                                  fe_v_neighbor.shape_value(l,point) *
-                                  fe_v_neighbor.shape_value(k,point) *
-                                  JxW[point];
-       }
-    }
-}
+                                  // As stated above, the
+                                  // <code>postprocess</code> function takes
+                                  // an evaluation object, and
+                                  // applies it to the computed
+                                  // solution. This function may be
+                                  // called multiply, once for each
+                                  // evaluation of the solution which
+                                  // the user required.
+  template <int dim>
+  void
+  Solver<dim>::
+  postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const
+  {
+    postprocessor (dof_handler, solution);
+  }
 
 
-                                // @sect3{Class: DGMethod}
-                                //
-                                // After these preparations, we
-                                // proceed with the main part of this
-                                // program. The main class, here
-                                // called <code>DGMethod</code> is basically
-                                // the main class of step 6. One of
-                                // the differences is that there's no
-                                // ConstraintMatrix object. This is,
-                                // because there are no hanging node
-                                // constraints in DG discretizations.
-template <int dim>
-class DGMethod
-{
-  public:
-    DGMethod ();
-    ~DGMethod ();
+                                  // The <code>n_dofs</code> function should
+                                  // be self-explanatory:
+  template <int dim>
+  unsigned int
+  Solver<dim>::n_dofs () const
+  {
+    return dof_handler.n_dofs();
+  }
+  
 
-    void run ();
-    
-  private:
-    void setup_system ();
-    void assemble_system1 ();
-    void assemble_system2 ();
-    void solve (Vector<double> &solution);
-    void refine_grid ();
-    void output_results (const unsigned int cycle) const;
-
-    Triangulation<dim>   triangulation;
+                                  // The following function assembles
+                                  // matrix and right hand side of
+                                  // the linear system to be solved
+                                  // in each step. It goes along the
+                                  // same lines as used in previous
+                                  // examples, so we explain it only
+                                  // briefly:
+  template <int dim>
+  void
+  Solver<dim>::assemble_linear_system (LinearSystem &linear_system)
+  {
+                                    // First define a convenience
+                                    // abbreviation for these lengthy
+                                    // iterator names...
+    typedef
+      typename DoFHandler<dim>::active_cell_iterator
+      active_cell_iterator;
+
+                                    // ... and use it to split up the
+                                    // set of cells into a number of
+                                    // pieces of equal size. The
+                                    // number of blocks is set to the
+                                    // default number of threads to
+                                    // be used, which by default is
+                                    // set to the number of
+                                    // processors found in your
+                                    // computer at startup of the
+                                    // program:
+    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);
+
+                                    // These ranges are then assigned
+                                    // to a number of threads which
+                                    // we create next. Each will
+                                    // assemble the local cell
+                                    // matrices on the assigned
+                                    // cells, and fill the matrix
+                                    // object with it. Since there is
+                                    // need for synchronization when
+                                    // filling the same matrix from
+                                    // different threads, we need a
+                                    // mutex here:
+    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);
+
+                                    // While the spawned threads
+                                    // assemble the system matrix, we
+                                    // can already compute the right
+                                    // hand side vector in the main
+                                    // thread, and condense away the
+                                    // constraints due to hanging
+                                    // nodes:
+    assemble_rhs (linear_system.rhs);
+    linear_system.hanging_node_constraints.condense (linear_system.rhs);
+
+                                    // And while we're already at it
+                                    // to compute things in parallel,
+                                    // interpolating boundary values
+                                    // is one more thing that can be
+                                    // done independently, so we do
+                                    // it here:
+    std::map<unsigned int,double> boundary_value_map;
+    VectorTools::interpolate_boundary_values (dof_handler,
+                                             0,
+                                             *boundary_values,
+                                             boundary_value_map);
     
-                                    // In contrast to the example code
-                                     // of step-12, this time DG elements
-                                    // of different degree will be used.
-                                    // The different FiniteElement
-                                     // objects for the different polynomial
-                                     // degrees will be stored in the
-                                     // fe_collection object.
-    hp::FECollection<dim>    fe_collection;
-
-                                    // As already mentioned, the
-                                     // standard DoFHandler has to be
-                                     // replaced by a hp::DoFHandler.
-    hp::DoFHandler<dim>    dof_handler;
-
-    SparsityPattern sparsity;
-    SparseMatrix<double> system_matrix;
-
-                                    // We define the quadrature
-                                    // formulae for the cell and the
-                                    // face terms of the
-                                    // discretization.
-                                     // Clearly the hp-Method requires
-                                     // a complete set of quadrature
-                                     // rules for each polynomial
-                                     // degree which will be used in the
-                                     // computations.
-    hp::QCollection<dim>   quadratures;
-    hp::QCollection<dim-1> face_quadratures;
     
-                                    // And there are two solution
-                                    // vectors, that store the
-                                    // solutions to the problems
-                                    // corresponding to the two
-                                    // different assembling routines
-                                    // <code>assemble_system1</code> and
-                                    // <code>assemble_system2</code>;
-    Vector<double>       solution1;
-    Vector<double>       solution2;
-    Vector<double>       right_hand_side;
+                                    // If this is done, wait for the
+                                    // matrix assembling threads, and
+                                    // condense the constraints in
+                                    // the matrix as well:
+    threads.join_all ();
+    linear_system.hanging_node_constraints.condense (linear_system.matrix);
+
+                                    // Now that we have the linear
+                                    // system, we can also treat
+                                    // boundary values, which need to
+                                    // be eliminated from both the
+                                    // matrix and the right hand
+                                    // side:
+    MatrixTools::apply_boundary_values (boundary_value_map,
+                                       linear_system.matrix,
+                                       solution,
+                                       linear_system.rhs);
+
+  }
+
+
+                                  // The second of this pair of
+                                  // functions takes a range of cell
+                                  // iterators, and assembles the
+                                  // system matrix on this part of
+                                  // the domain. Since it's actions
+                                  // have all been explained in
+                                  // previous programs, we do not
+                                  // comment on it any more, except
+                                  // for one pointe below.
+  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
+  {
+    FEValues<dim> fe_values (*fe, *quadrature, 
+                            update_gradients | 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);
+
+    for (typename DoFHandler<dim>::active_cell_iterator cell=begin_cell;
+        cell!=end_cell; ++cell)
+      {
+       cell_matrix = 0;
+
+       fe_values.reinit (cell);
+
+       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) *
+                                  fe_values.JxW(q_point));
+
+
+       cell->get_dof_indices (local_dof_indices);
+
+                                         // In the step-9 program, we
+                                         // have shown that you have
+                                         // to use the mutex to lock
+                                         // the matrix when copying
+                                         // the elements from the
+                                         // local to the global
+                                         // matrix. This was necessary
+                                         // to avoid that two threads
+                                         // access it at the same
+                                         // time, eventually
+                                         // overwriting their
+                                         // respective
+                                         // work. Previously, we have
+                                         // used the <code>acquire</code> and
+                                         // <code>release</code> functions of
+                                         // the mutex to lock and
+                                         // unlock the mutex,
+                                         // respectively. While this
+                                         // is valid, there is one
+                                         // possible catch: if between
+                                         // the locking operation and
+                                         // the unlocking operation an
+                                         // exception is thrown, the
+                                         // mutex remains in the
+                                         // locked state, and in some
+                                         // cases this might lead to
+                                         // deadlocks. A similar
+                                         // situation arises, when one
+                                         // changes the code to have a
+                                         // return statement somewhere
+                                         // in the middle of the
+                                         // locked block, and forgets
+                                         // that before we call
+                                         // <code>return</code>, we also have
+                                         // to unlock the mutex. This
+                                         // all is not be a problem
+                                         // here, but we want to show
+                                         // the general technique to
+                                         // cope with these problems
+                                         // nevertheless: have an
+                                         // object that upon
+                                         // initialization (i.e. in
+                                         // its constructor) locks the
+                                         // mutex, and on running the
+                                         // destructor unlocks it
+                                         // again. This is called the
+                                         // <code>scoped lock</code> pattern
+                                         // (apparently invented by
+                                         // Doug Schmidt originally),
+                                         // and it works because
+                                         // destructors of local
+                                         // objects are also run when
+                                         // we exit the function
+                                         // either through a
+                                         // <code>return</code> statement, or
+                                         // when an exception is
+                                         // raised. Thus, it is
+                                         // guaranteed that the mutex
+                                         // will always be unlocked
+                                         // when we exit this part of
+                                         // the program, whether the
+                                         // operation completed
+                                         // successfully or not,
+                                         // whether the exit path was
+                                         // something we implemented
+                                         // willfully or whether the
+                                         // function was exited by an
+                                         // exception that we did not
+                                         // forsee.
+                                         //
+                                         // deal.II implements the
+                                         // scoped locking pattern in
+                                         // the
+                                         // ThreadMutex::ScopedLock
+                                         // class: it takes the mutex
+                                         // in the constructor and
+                                         // locks it; in its
+                                         // destructor, it unlocks it
+                                         // again. So here is how it
+                                         // is used:
+        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));
+                                         // Here, at the brace, the
+                                         // current scope ends, so the
+                                         // <code>lock</code> variable goes out
+                                         // of existence and its
+                                         // destructor the mutex is
+                                         // unlocked.
+      };
+  }
+
+
+                                  // Now for the functions that
+                                  // implement actions in the linear
+                                  // system class. First, the
+                                  // constructor initializes all data
+                                  // elements to their correct sizes,
+                                  // and sets up a number of
+                                  // additional data structures, such
+                                  // as constraints due to hanging
+                                  // nodes. Since setting up the
+                                  // hanging nodes and finding out
+                                  // about the nonzero elements of
+                                  // the matrix is independent, we do
+                                  // that in parallel (if the library
+                                  // was configured to use
+                                  // concurrency, at least;
+                                  // otherwise, the actions are
+                                  // performed sequentially). Note
+                                  // that we spawn only one thread,
+                                  // and do the second action in the
+                                  // main thread. Since only one
+                                  // thread is generated, we don't
+                                  // use the <code>Threads::ThreadGroup</code>
+                                  // class here, but rather use the
+                                  // one created thread object
+                                  // directly to wait for this
+                                  // particular thread's exit.
+                                  //
+                                  // Note that taking up the address
+                                  // of the
+                                  // <code>DoFTools::make_hanging_node_constraints</code>
+                                  // function is a little tricky,
+                                  // since there are actually three
+                                  // of them, one for each supported
+                                  // space dimension. Taking
+                                  // addresses of overloaded
+                                  // functions is somewhat
+                                  // complicated in C++, since the
+                                  // address-of operator <code>&</code> in
+                                  // that case returns more like a
+                                  // set of values (the addresses of
+                                  // all functions with that name),
+                                  // and selecting the right one is
+                                  // then the next step. If the
+                                  // context dictates which one to
+                                  // take (for example by assigning
+                                  // to a function pointer of known
+                                  // type), then the compiler can do
+                                  // that by itself, but if this set
+                                  // of pointers shall be given as
+                                  // the argument to a function that
+                                  // takes a template, the compiler
+                                  // could choose all without having
+                                  // a preference for one. We
+                                  // therefore have to make it clear
+                                  // to the compiler which one we
+                                  // would like to have; for this, we
+                                  // could use a cast, but for more
+                                  // clarity, we assign it to a
+                                  // temporary <code>mhnc_p</code> (short for
+                                  // <code>pointer to
+                                  // make_hanging_node_constraints</code>)
+                                  // with the right type, and using
+                                  // this pointer instead.
+  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;
     
-                                    // Finally this class includes an
-                                    // object of the
-                                    // DGTransportEquations class
-                                    // described above.
-    const DGTransportEquation<dim> dg;
-};
+    Threads::Thread<>
+      mhnc_thread = Threads::spawn (mhnc_p)(dof_handler,
+                                            hanging_node_constraints);
+
+    sparsity_pattern.reinit (dof_handler.n_dofs(),
+                            dof_handler.n_dofs(),
+                            dof_handler.max_couplings_between_dofs());
+    DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
+
+                                    // Wait until the
+                                    // <code>hanging_node_constraints</code>
+                                    // object is fully set up, then
+                                    // close it and use it to
+                                    // condense the sparsity pattern:
+    mhnc_thread.join ();
+    hanging_node_constraints.close ();
+    hanging_node_constraints.condense (sparsity_pattern);
+
+                                    // Finally, close the sparsity
+                                    // pattern, initialize the
+                                    // matrix, and set the right hand
+                                    // side vector to the right size.
+    sparsity_pattern.compress();
+    matrix.reinit (sparsity_pattern);
+    rhs.reinit (dof_handler.n_dofs());
+  }
 
 
-template <int dim>
-DGMethod<dim>::DGMethod ()
-               :
-               dof_handler (triangulation),
-               dg ()
-{
-                                  // Change here for hp
-                                  // methods of
-                                  // different maximum degrees.
-  const unsigned int hp_degree = 5;
-  for (unsigned int i = 0; i < hp_degree; ++i)
-    {
-      fe_collection.push_back (FE_DGQ<dim> (i));
-      quadratures.push_back (QGauss<dim> (i+2));
-      face_quadratures.push_back (QGauss<dim-1> (i+2));
-    }
-}
 
+                                  // The second function of this
+                                  // class simply solves the linear
+                                  // system by a preconditioned
+                                  // conjugate gradient method. This
+                                  // has been extensively discussed
+                                  // before, so we don't dwell into
+                                  // it any more.
+  template <int dim>
+  void
+  Solver<dim>::LinearSystem::solve (Vector<double> &solution) const
+  {
+    SolverControl           solver_control (1000, 1e-12);
+    SolverCG<>              cg (solver_control);
 
-template <int dim>
-DGMethod<dim>::~DGMethod () 
-{
-  dof_handler.clear ();
-}
+    PreconditionSSOR<> preconditioner;
+    preconditioner.initialize(matrix, 1.2);
 
+    cg.solve (matrix, solution, rhs, preconditioner);
 
-template <int dim>
-void DGMethod<dim>::setup_system ()
-{
-                                  // First we need to distribute the
-                                  // DoFs.
-  dof_handler.distribute_dofs (fe_collection);
-                                  // In order to get a good
-                                  // preconditioner, the degrees of
-                                  // freedom should be ordered in
-                                  // downstream direction.  First, we
-                                  // initalize a vector fairly close
-                                  // to the real vector field; since
-                                  // this is for preconditioning
-                                  // only, a rough approximation is
-                                  // sufficient.
-  Point<dim> sorting_direction;
-  for (unsigned int d=0;d<dim;++d)
-    sorting_direction(d) = 1.;
-                                  // Now do the sorting of the
-                                  // degrees of freedom.
-  DoFRenumbering::downstream_dg(dof_handler, sorting_direction);
-  
-                                  // The DoFs of a cell are coupled
-                                  // with all DoFs of all neighboring
-                                  // cells.  Therefore the maximum
-                                  // number of matrix entries per row
-                                  // is needed when all neighbors of
-                                  // a cell are once more refined
-                                  // than the cell under
-                                  // consideration.
-  CompressedSparsityPattern compressed_pattern (dof_handler.n_dofs ());
-  DoFTools::make_sparsity_pattern (dof_handler, compressed_pattern);
-  DoFTools::make_flux_sparsity_pattern (dof_handler, compressed_pattern);
-
-  sparsity.copy_from(compressed_pattern);
-  system_matrix.reinit (sparsity);
-
-  solution1.reinit (dof_handler.n_dofs());
-  solution2.reinit (dof_handler.n_dofs());
-  right_hand_side.reinit (dof_handler.n_dofs());
-}
+    hanging_node_constraints.distribute (solution);
+  }
 
 
-                                // @sect4{Function: assemble_system1}
-                                //
-                                // We proceed with the
-                                // <code>assemble_system1</code> function that
-                                // implements the DG discretization
-                                // in its first version. This
-                                // function repeatedly calls the
-                                // <code>assemble_cell_term</code>,
-                                // <code>assemble_boundary_term</code> and
-                                // <code>assemble_face_term1</code> functions
-                                // of the <code>DGTransportEquation</code>
-                                // object.  The
-                                // <code>assemble_boundary_term</code> covers
-                                // the first case mentioned in the
-                                // introduction.
-                                //
-                                // 1. face is at boundary
-                                //
-                                // This function takes a
-                                // <code>FEFaceValues</code> object as
-                                // argument.  In contrast to that
-                                // <code>assemble_face_term1</code>
-                                // takes two <code>FEFaceValuesBase</code>
-                                // objects; one for the shape
-                                // functions on the current cell and
-                                // the other for shape functions on
-                                // the neighboring cell under
-                                // consideration. Both objects are
-                                // either of class <code>FEFaceValues</code>
-                                // or of class <code>FESubfaceValues</code>
-                                // (both derived from
-                                // <code>FEFaceValuesBase</code>) according to
-                                // the remaining cases mentioned
-                                // in the introduction:
-                                //
-                                // 2. neighboring cell is finer
-                                // (current cell: <code>FESubfaceValues</code>,
-                                // neighboring cell: <code>FEFaceValues</code>);
-                                //
-                                // 3. neighboring cell is of the same
-                                // refinement level (both, current
-                                // and neighboring cell:
-                                // <code>FEFaceValues</code>);
-                                //
-                                // 4. neighboring cell is coarser
-                                // (current cell: <code>FEFaceValues</code>,
-                                // neighboring cell:
-                                // <code>FESubfaceValues</code>).
-                                //
-                                // If we considered globally refined
-                                // meshes then only case 3 would
-                                // occur. But as we consider also
-                                // locally refined meshes we need to
-                                // distinguish all four cases making
-                                // the following assembling function
-                                // a bit longish.
-template <int dim>
-void DGMethod<dim>::assemble_system1 () 
-{
-                                  // First we create the
-                                  // <code>UpdateFlags</code> for the
-                                  // <code>FEValues</code> and the
-                                  // <code>FEFaceValues</code> objects.
-  const UpdateFlags update_flags = update_values
-                                   | update_gradients
-                                   | update_q_points
-                                   | update_JxW_values;
-
-                                  // Note, that on faces we do not
-                                  // need gradients but we need
-                                  // normal vectors.
-  const UpdateFlags face_update_flags = update_values
-                                        | update_q_points
-                                        | update_JxW_values
-                                        | update_normal_vectors;
-  
-                                  // On the neighboring cell we only
-                                  // need the shape values. Given a
-                                  // specific face, the quadrature
-                                  // points and `JxW values' are the
-                                  // same as for the current cells,
-                                  // the normal vectors are known to
-                                  // be the negative of the normal
-                                  // vectors of the current cell.
-  const UpdateFlags neighbor_face_update_flags = update_values;
-   
-                                  // Then we create the <code>FEValues</code>
-                                  // object. Here, we use the default
-                                  // MappingQ1. different mapping
-                                  // create a MappingCollection first
-                                  // and call the respective
-                                  // hp::FEValues constructor.
-  hp::FEValues<dim> fe_v_x (fe_collection, quadratures, update_flags);
-  
-                                  // Similarly we create the
-                                  // <code>FEFaceValues</code> and
-                                  // <code>FESubfaceValues</code> objects for
-                                  // both, the current and the
-                                  // neighboring cell. Within the
-                                  // following nested loop over all
-                                  // cells and all faces of the cell
-                                  // they will be reinited to the
-                                  // current cell and the face (and
-                                  // subface) number.
-  hp::FEFaceValues<dim> fe_v_face_x (
-    fe_collection, face_quadratures, face_update_flags);
-  hp::FESubfaceValues<dim> fe_v_subface_x (
-    fe_collection, face_quadratures, face_update_flags);
-  hp::FEFaceValues<dim> fe_v_face_neighbor_x (
-    fe_collection, face_quadratures, neighbor_face_update_flags);
-  hp::FESubfaceValues<dim> fe_v_subface_neighbor_x (
-    fe_collection, face_quadratures, neighbor_face_update_flags);
-
-                                  // Now we create the cell matrices
-                                  // and vectors. Here we need two
-                                  // cell matrices, both for face
-                                  // terms that include test
-                                  // functions <code>v</code> (shape functions
-                                  // of the current cell). To be more
-                                  // precise, the first matrix will
-                                  // include the `u and v terms' and
-                                  // the second that will include the
-                                  // `un and v terms'. Here we recall
-                                  // the convention that `un' is
-                                  // the shortcut for `u_neighbor'
-                                  // and represents the $u_hat$, see
-                                  // introduction.
-  const unsigned int max_dofs_per_cell = fe_collection.max_dofs_per_cell ();
-
-  FullMatrix<double> u_v_matrix (max_dofs_per_cell, max_dofs_per_cell);
-  FullMatrix<double> un_v_matrix (max_dofs_per_cell, max_dofs_per_cell);
-  Vector<double>  cell_vector (max_dofs_per_cell);
-
-                                  // Furthermore we need some cell
-                                  // iterators.
-  typename hp::DoFHandler<dim>::active_cell_iterator
-    cell = dof_handler.begin_active(),
-    endc = dof_handler.end();
-
-                                  // Now we start the loop over all
-                                  // active cells.
-  for (;cell!=endc; ++cell) 
-    {
-      const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
-      std::vector<unsigned int> dofs (dofs_per_cell);
-      std::vector<unsigned int> dofs_neighbor;
-
-                                      // In the
-                                      // <code>assemble_face_term1</code>
-                                      // function contributions to
-                                      // the cell matrices and the
-                                      // cell vector are only
-                                      // ADDED. Therefore on each
-                                      // cell we need to reset the
-                                      // <code>u_v_matrix</code> and
-                                      // <code>cell_vector</code> to zero,
-                                      // before assembling the cell terms.
-      u_v_matrix = 0;
-      cell_vector = 0;
-      
-                                      // Now we reinit the <code>FEValues</code>
-                                      // object for the current cell
-      fe_v_x.reinit (cell);
-
-                                      // and call the function
-                                      // that assembles the cell
-                                      // terms. The first argument is
-                                      // the <code>FEValues</code> that was
-                                      // previously reinit'ed on the
-                                      // current cell.
-      dg.assemble_cell_term(fe_v_x.get_present_fe_values (),
-                           u_v_matrix,
-                           cell_vector);
-
-                                      // As in previous examples the
-                                      // vector `dofs' includes the
-                                      // dof_indices.
-      dofs.resize (dofs_per_cell);
-      cell->get_dof_indices (dofs);
-
-                                      // This is the start of the
-                                      // nested loop over all faces.
-      for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
-       {
-                                          // First we set the face
-                                          // iterator
-         typename hp::DoFHandler<dim>::face_iterator face=cell->face(face_no);
-         
-                                          // and clear the
-                                          // <code>un_v_matrix</code> on each
-                                          // face.
-         un_v_matrix = 0;
-
-                                          // Now we distinguish the
-                                          // four different cases in
-                                          // the ordering mentioned
-                                          // above. We start with
-                                          // faces belonging to the
-                                          // boundary of the domain.
-         if (face->at_boundary())
-           {
-                                              // We reinit the
-                                              // <code>FEFaceValues</code>
-                                              // object to the
-                                              // current face
-             fe_v_face_x.reinit (cell, face_no);
-
-                                              // and assemble the
-                                              // corresponding face
-                                              // terms.
-             dg.assemble_boundary_term(fe_v_face_x.get_present_fe_values (),
-                                       u_v_matrix,
-                                       cell_vector);
-           }
-         else
-           {
-                                              // Now we are not on
-                                              // the boundary of the
-                                              // domain, therefore
-                                              // there must exist a
-                                              // neighboring cell.
-             typename hp::DoFHandler<dim>::cell_iterator neighbor=
-               cell->neighbor(face_no);
-             
-                                              // We proceed with the
-                                              // second and most
-                                              // complicated case:
-                                              // the neighboring cell
-                                              // is more refined than
-                                              // the current cell. As
-                                              // in deal.II
-                                              // neighboring cells
-                                              // are restricted to
-                                              // have a level
-                                              // difference of not
-                                              // more than one, the
-                                              // neighboring cell is
-                                              // known to be at most
-                                              // ONCE more refined
-                                              // than the current
-                                              // cell. Furthermore
-                                              // also the face is
-                                              // more refined,
-                                              // i.e. it has
-                                              // children. Here we
-                                              // note that the
-                                              // following part of
-                                              // code will not work
-                                              // for <code>dim==1</code>.
-             if (face->has_children())
-               {
-                                                  // First we store
-                                                  // which number the
-                                                  // current cell has
-                                                  // in the list of
-                                                  // neighbors of the
-                                                  // neighboring
-                                                  // cell. Hence,
-                                                  // neighbor-@>neighbor(neighbor2)
-                                                  // equals the
-                                                  // current cell
-                                                  // <code>cell</code>.
-                 const unsigned int neighbor2=
-                   cell->neighbor_of_neighbor(face_no);
-                 
-                 
-                                                  // We loop over
-                                                  // subfaces
-                 for (unsigned int subface_no=0;
-                      subface_no<face->n_children(); ++subface_no)
-                   {
-                                                      // and set the
-                                                      // cell
-                                                      // iterator
-                                                      // <code>neighbor_child</code>
-                                                      // to the cell
-                                                      // placed
-                                                      // `behind' the
-                                                      // current
-                                                      // subface.
-                     typename hp::DoFHandler<dim>::active_cell_iterator neighbor_child=
-                       neighbor->child(GeometryInfo<dim>::
-                                       child_cell_on_face(neighbor2,subface_no));
-
-                                                      // an additional speciality
-                                                      // for the hp method appears
-                                                      // on the faces. To get an
-                                                      // efficient assembly, the
-                                                      // lowest order but 
-                                                      // sufficient quadrature
-                                                      // rule should be used. Hence
-                                                      // the face quadrature rule of the
-                                                      // higher order element
-                                                      // will be used.
-                     const unsigned int quadrature_index = 
-                       std::max (neighbor_child->active_fe_index (),
-                                  cell->active_fe_index ());
-                       
-
-                                                      // As these are
-                                                      // quite
-                                                      // complicated
-                                                      // indirections
-                                                      // which one
-                                                      // does not
-                                                      // usually get
-                                                      // right at
-                                                      // first
-                                                      // attempt we
-                                                      // check for
-                                                      // the internal
-                                                      // consistency.
-                     Assert (neighbor_child->face(neighbor2) == face->child(subface_no),
-                             ExcInternalError());
-                     Assert (!neighbor_child->has_children(), ExcInternalError());
-
-                                                      // We need to
-                                                      // reset the
-                                                      // <code>un_v_matrix</code>
-                                                      // on each
-                                                      // subface
-                                                      // because on
-                                                      // each subface
-                                                      // the <code>un</code>
-                                                      // belong to
-                                                      // different
-                                                      // neighboring
-                                                      // cells.
-                     un_v_matrix = 0;
-                     
-                                                      // As already
-                                                      // mentioned
-                                                      // above for
-                                                      // the current
-                                                      // case (case
-                                                      // 2) we employ
-                                                      // the
-                                                      // <code>FESubfaceValues</code>
-                                                      // of the
-                                                      // current
-                                                      // cell (here
-                                                      // reinited for
-                                                      // the current
-                                                      // cell, face
-                                                      // and subface)
-                                                      // and we
-                                                      // employ the
-                                                      // FEFaceValues
-                                                      // of the
-                                                      // neighboring
-                                                      // child cell.
-                     fe_v_subface_x.reinit (cell, face_no, subface_no, quadrature_index);
-                     fe_v_face_neighbor_x.reinit (neighbor_child, neighbor2, quadrature_index);
-
-                     dg.assemble_face_term1(fe_v_subface_x.get_present_fe_values (),
-                                            fe_v_face_neighbor_x.get_present_fe_values (),
-                                            u_v_matrix,
-                                            un_v_matrix);
-                     
-                                                      // Then we get
-                                                      // the dof
-                                                      // indices of
-                                                      // the
-                                                      // neighbor_child
-                                                      // cell
-                     dofs_neighbor.resize (neighbor_child->get_fe().dofs_per_cell);
-                     neighbor_child->get_dof_indices (dofs_neighbor);
-                                                               
-                                                      // and
-                                                      // distribute
-                                                      // <code>un_v_matrix</code>
-                                                      // to the
-                                                      // system_matrix
-                     for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
-                       for (unsigned int k=0; k<neighbor_child->get_fe().dofs_per_cell; ++k)
-                         system_matrix.add(dofs[i], dofs_neighbor[k],
-                                           un_v_matrix(i,k));
-                   }
-                                                  // End of <code>if
-                                                  // (face-@>has_children())</code>
-               }
-             else
-               {
-                                                  // We proceed with
-                                                  // case 3,
-                                                  // i.e. neighboring
-                                                  // cell is of the
-                                                  // same refinement
-                                                  // level as the
-                                                  // current cell.
-                 if (neighbor->level() == cell->level()) 
-                   {
-                                                      // Like before
-                                                      // we store
-                                                      // which number
-                                                      // the current
-                                                      // cell has in
-                                                      // the list of
-                                                      // neighbors of
-                                                      // the
-                                                      // neighboring
-                                                      // cell.
-                     const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no);
-
-                                                      // Like before. Use
-                                                      // quadrature rule
-                                                       // of higher order
-                                                       // cell.
-                     const unsigned int quadrature_index = 
-                       std::max (neighbor->active_fe_index (),
-                                  cell->active_fe_index ());
-
-                                                      // We reinit
-                                                      // the
-                                                      // <code>FEFaceValues</code>
-                                                      // of the
-                                                      // current and
-                                                      // neighboring
-                                                      // cell to the
-                                                      // current face
-                                                      // and assemble
-                                                      // the
-                                                      // corresponding
-                                                      // face terms.
-                     fe_v_face_x.reinit (cell, face_no, quadrature_index);
-                     fe_v_face_neighbor_x.reinit (neighbor, neighbor2, quadrature_index);
-                     
-                     dg.assemble_face_term1(fe_v_face_x.get_present_fe_values (),
-                                            fe_v_face_neighbor_x.get_present_fe_values (),
-                                            u_v_matrix,
-                                            un_v_matrix);
-                                                      // End of <code>if
-                                                      // (neighbor-@>level()
-                                                      // ==
-                                                      // cell-@>level())</code>
-                   }
-                 else
-                   {
-                                                      // Finally we
-                                                      // consider
-                                                      // case 4. When
-                                                      // the
-                                                      // neighboring
-                                                      // cell is not
-                                                      // finer and
-                                                      // not of the
-                                                      // same
-                                                      // refinement
-                                                      // level as the
-                                                      // current cell
-                                                      // it must be
-                                                      // coarser.
-                     Assert(neighbor->level() < cell->level(), ExcInternalError());
-
-                                                      // Find out the
-                                                      // how many'th
-                                                      // face_no and
-                                                      // subface_no
-                                                      // the current
-                                                      // face is
-                                                      // w.r.t. the
-                                                      // neighboring
-                                                      // cell.
-                     const std::pair<unsigned int, unsigned int> faceno_subfaceno=
-                       cell->neighbor_of_coarser_neighbor(face_no);
-                     const unsigned int neighbor_face_no=faceno_subfaceno.first,
-                                     neighbor_subface_no=faceno_subfaceno.second;
-
-                     Assert (neighbor->neighbor(neighbor_face_no)
-                             ->child(GeometryInfo<dim>::child_cell_on_face(
-                               face_no,neighbor_subface_no)) == cell, ExcInternalError());
-
-
-                                                      // Like before. Use
-                                                      // quadrature rule
-                                                       // of higher order
-                                                       // cell.
-                     const unsigned int quadrature_index = 
-                       std::max (neighbor->active_fe_index (),
-                                  cell->active_fe_index ());
-
-                                                      // Reinit the
-                                                      // appropriate
-                                                      // <code>FEFaceValues</code>
-                                                      // and assemble
-                                                      // the face
-                                                      // terms.
-                     fe_v_face_x.reinit (cell, face_no, quadrature_index);
-                     fe_v_subface_neighbor_x.reinit (neighbor, neighbor_face_no,
-                                                     neighbor_subface_no, quadrature_index);
-                     
-                     dg.assemble_face_term1(fe_v_face_x.get_present_fe_values (),
-                                            fe_v_subface_neighbor_x.get_present_fe_values (),
-                                            u_v_matrix,
-                                            un_v_matrix);
-                   }
-
-                                                  // Now we get the
-                                                  // dof indices of
-                                                  // the
-                                                  // <code>neighbor_child</code>
-                                                  // cell,
-                 dofs_neighbor.resize (neighbor->get_fe().dofs_per_cell);
-                 neighbor->get_dof_indices (dofs_neighbor);
-
-                                                  // and distribute the
-                                                  // <code>un_v_matrix</code>.
-                 for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
-                   for (unsigned int k=0; k<neighbor->get_fe().dofs_per_cell; ++k)
-                     system_matrix.add(dofs[i], dofs_neighbor[k],
-                                       un_v_matrix(i,k));
-               }
-                                              // End of <code>face not at boundary</code>:
-           }
-                                          // End of loop over all faces:
-       }
-      
-                                      // Finally we distribute the
-                                      // <code>u_v_matrix</code>
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-       for (unsigned int j=0; j<dofs_per_cell; ++j)
-         system_matrix.add(dofs[i], dofs[j], u_v_matrix(i,j));
+
+
+                                  // @sect4{A primal solver}
+
+                                  // In the previous section, a base
+                                  // class for Laplace solvers was
+                                  // implemented, that lacked the
+                                  // functionality to assemble the
+                                  // right hand side vector, however,
+                                  // for reasons that were explained
+                                  // there. Now we implement a
+                                  // corresponding class that can do
+                                  // this for the case that the right
+                                  // hand side of a problem is given
+                                  // as a function object.
+                                  //
+                                  // The actions of the class are
+                                  // rather what you have seen
+                                  // already in previous examples
+                                  // already, so a brief explanation
+                                  // should suffice: the constructor
+                                  // takes the same data as does that
+                                  // of the underlying class (to
+                                  // which it passes all information)
+                                  // except for one function object
+                                  // that denotes the right hand side
+                                  // of the problem. A pointer to
+                                  // this object is stored (again as
+                                  // a <code>SmartPointer</code>, in order to
+                                  // make sure that the function
+                                  // object is not deleted as long as
+                                  // it is still used by this class).
+                                  //
+                                  // The only functional part of this
+                                  // class is the <code>assemble_rhs</code>
+                                  // method that does what its name
+                                  // suggests.
+  template <int dim>
+  class PrimalSolver : public Solver<dim>
+  {
+    public:
+      PrimalSolver (Triangulation<dim>       &triangulation,
+                   const FiniteElement<dim> &fe,
+                   const Quadrature<dim>    &quadrature,
+                   const Function<dim>      &rhs_function,
+                   const Function<dim>      &boundary_values);
+    protected:
+      const SmartPointer<const Function<dim> > rhs_function;
+      virtual void assemble_rhs (Vector<double> &rhs) const;
+  };
+
+
+                                  // The constructor of this class
+                                  // basically does what it is
+                                  // announced to do above...
+  template <int dim>
+  PrimalSolver<dim>::
+  PrimalSolver (Triangulation<dim>       &triangulation,
+               const FiniteElement<dim> &fe,
+               const Quadrature<dim>    &quadrature,
+               const Function<dim>      &rhs_function,
+               const Function<dim>      &boundary_values)
+                 :
+                 Base<dim> (triangulation),
+                 Solver<dim> (triangulation, fe,
+                              quadrature, boundary_values),
+                  rhs_function (&rhs_function)
+  {}
+
+
+
+                                  // ... as does the <code>assemble_rhs</code>
+                                  // function. Since this is
+                                  // explained in several of the
+                                  // previous example programs, we
+                                  // leave it at that.
+  template <int dim>
+  void
+  PrimalSolver<dim>::
+  assemble_rhs (Vector<double> &rhs) const 
+  {
+    FEValues<dim> fe_values (*this->fe, *this->quadrature, 
+                            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);
       
-                                      // and the cell vector.
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-       right_hand_side(dofs[i]) += cell_vector(i);
-    }
-}
+       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{Function: assemble_system2}
-                                //
-                                // We proceed with the
-                                // <code>assemble_system2</code> function that
-                                // implements the DG discretization
-                                // in its second version. This
-                                // function is very similar to the
-                                // <code>assemble_system1</code>
-                                // function. Therefore, here we only
-                                // discuss the differences between
-                                // the two functions. This function
-                                // repeatedly calls the
-                                // <code>assemble_face_term2</code> function
-                                // of the DGTransportEquation object,
-                                // that assembles the face terms
-                                // written as a sum of integrals over
-                                // all faces. Therefore, we need to
-                                // make sure that each face is
-                                // treated only once. This is achieved
-                                // by introducing the rule:
-                                // 
-                                // a) If the current and the
-                                // neighboring cells are of the same
-                                // refinement level we access and
-                                // treat the face from the cell with
-                                // lower index.
-                                //
-                                // b) If the two cells are of
-                                // different refinement levels we
-                                // access and treat the face from the
-                                // coarser cell.
-                                //
-                                // Due to rule b) we do not need to
-                                // consider case 4 (neighboring cell
-                                // is coarser) any more.
+                                  // @sect4{Local refinement by the Kelly error indicator}
+
+                                  // The second class implementing
+                                  // refinement strategies uses the
+                                  // Kelly refinemet indicator used
+                                  // in various example programs
+                                  // before. Since this indicator is
+                                  // already implemented in a class
+                                  // of its own inside the deal.II
+                                  // library, there is not much t do
+                                  // here except cal the function
+                                  // computing the indicator, then
+                                  // using it to select a number of
+                                  // cells for refinement and
+                                  // coarsening, and refinement the
+                                  // mesh accordingly.
+                                  //
+                                  // Again, this should now be
+                                  // sufficiently standard to allow
+                                  // the omission of further
+                                  // comments.
+  template <int dim>
+  class RefinementKelly : public PrimalSolver<dim>
+  {
+    public:
+      RefinementKelly (Triangulation<dim>       &coarse_grid,
+                      const FiniteElement<dim> &fe,
+                      const Quadrature<dim>    &quadrature,
+                      const Function<dim>      &rhs_function,
+                      const Function<dim>      &boundary_values);
 
-template <int dim>
-void DGMethod<dim>::assemble_system2 () 
-{
-  const UpdateFlags update_flags = update_values
-                                   | update_gradients
-                                   | update_q_points
-                                   | update_JxW_values;
-  
-  const UpdateFlags face_update_flags = update_values
-                                        | update_q_points
-                                        | update_JxW_values
-                                        | update_normal_vectors;
-  
-  const UpdateFlags neighbor_face_update_flags = update_values;
-
-                                  // Here we do not need
-                                  // <code>fe_v_face_neighbor</code> as case 4
-                                  // does not occur.
-  hp::FEValues<dim> fe_v_x (
-    fe_collection, quadratures, update_flags);
-  hp::FEFaceValues<dim> fe_v_face_x (
-    fe_collection, face_quadratures, face_update_flags);
-  hp::FESubfaceValues<dim> fe_v_subface_x (
-    fe_collection, face_quadratures, face_update_flags);
-  hp::FEFaceValues<dim> fe_v_face_neighbor_x (
-    fe_collection, face_quadratures, neighbor_face_update_flags);
-
-  const unsigned int max_dofs_per_cell = fe_collection.max_dofs_per_cell ();
-
-  FullMatrix<double> u_v_matrix (max_dofs_per_cell, max_dofs_per_cell);
-  FullMatrix<double> un_v_matrix (max_dofs_per_cell, max_dofs_per_cell);
-  
-                                  // Additionally we need the
-                                  // following two cell matrices,
-                                  // both for face term that include
-                                  // test function <code>vn</code> (shape
-                                  // functions of the neighboring
-                                  // cell). To be more precise, the
-                                  // first matrix will include the `u
-                                  // and vn terms' and the second
-                                  // that will include the `un and vn
-                                  // terms'.
-  FullMatrix<double> u_vn_matrix (max_dofs_per_cell, max_dofs_per_cell);
-  FullMatrix<double> un_vn_matrix (max_dofs_per_cell, max_dofs_per_cell);
-  
-  Vector<double>  cell_vector (max_dofs_per_cell);
-
-                                  // The following lines are roughly
-                                  // the same as in the previous
-                                  // function.
-  typename hp::DoFHandler<dim>::active_cell_iterator
-    cell = dof_handler.begin_active(),
-    endc = dof_handler.end();
-  for (;cell!=endc; ++cell) 
-    {
-      const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
-      std::vector<unsigned int> dofs (dofs_per_cell);
-      std::vector<unsigned int> dofs_neighbor (dofs_per_cell);
+      virtual void refine_grid ();
+  };
 
-      u_v_matrix = 0;
-      cell_vector = 0;
 
-      fe_v_x.reinit (cell);
 
-      dg.assemble_cell_term(fe_v_x.get_present_fe_values (),
-                           u_v_matrix,
-                           cell_vector);
-      
-      cell->get_dof_indices (dofs);
+  template <int dim>
+  RefinementKelly<dim>::
+  RefinementKelly (Triangulation<dim>       &coarse_grid,
+                  const FiniteElement<dim> &fe,
+                  const Quadrature<dim>    &quadrature,
+                  const Function<dim>      &rhs_function,
+                  const Function<dim>      &boundary_values)
+                 :
+                 Base<dim> (coarse_grid),
+                  PrimalSolver<dim> (coarse_grid, fe, quadrature,
+                                    rhs_function, boundary_values)
+  {}
+
+
+
+  template <int dim>
+  void
+  RefinementKelly<dim>::refine_grid ()
+  {
+    Vector<float> estimated_error_per_cell (this->triangulation->n_active_cells());
+    KellyErrorEstimator<dim>::estimate (this->dof_handler,
+                                       QGauss<dim-1>(3),
+                                       typename FunctionMap<dim>::type(),
+                                       this->solution,
+                                       estimated_error_per_cell);
+    GridRefinement::refine_and_coarsen_fixed_number (*this->triangulation,
+                                                    estimated_error_per_cell,
+                                                    0.3, 0.03);
+    this->triangulation->execute_coarsening_and_refinement ();
+  }
 
-      for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
-       {
-         typename hp::DoFHandler<dim>::face_iterator face=
-           cell->face(face_no);
-
-                                          // Case 1:
-         if (face->at_boundary())
-           {
-             fe_v_face_x.reinit (cell, face_no);
-
-             dg.assemble_boundary_term(fe_v_face_x.get_present_fe_values (),
-                                       u_v_matrix,
-                                       cell_vector);
-           }
-         else
-           {
-             Assert (cell->neighbor(face_no).state() == IteratorState::valid,
-                     ExcInternalError());
-             typename hp::DoFHandler<dim>::cell_iterator neighbor=
-               cell->neighbor(face_no);
-
-             const unsigned int dofs_on_neighbor = neighbor->get_fe().dofs_per_cell;
-
-                                              // Case 2:
-             if (face->has_children())
-               {
-                 const unsigned int neighbor2=
-                   cell->neighbor_of_neighbor(face_no);
-                 
-                 for (unsigned int subface_no=0;
-                      subface_no<face->n_children(); ++subface_no)
-                   {
-                     typename hp::DoFHandler<dim>::cell_iterator neighbor_child=
-                       neighbor->child(GeometryInfo<dim>::child_cell_on_face(
-                         neighbor2,subface_no));
-                     const unsigned int dofs_on_neighbor_child = neighbor_child->get_fe().dofs_per_cell;
-
-                     const unsigned int quadrature_index = 
-                       std::max (neighbor_child->active_fe_index (),
-                                  cell->active_fe_index ());
-
-                     Assert (neighbor_child->face(neighbor2) == face->child(subface_no),
-                             ExcInternalError());
-                     Assert (!neighbor_child->has_children(), ExcInternalError());
-                     
-                     un_v_matrix = 0;
-                     u_vn_matrix = 0;
-                     un_vn_matrix = 0;
-                     
-                     fe_v_subface_x.reinit (cell, face_no, subface_no, quadrature_index);
-                     fe_v_face_neighbor_x.reinit (neighbor_child, neighbor2, quadrature_index);
-
-                     dg.assemble_face_term2(fe_v_subface_x.get_present_fe_values (),
-                                            fe_v_face_neighbor_x.get_present_fe_values (),
-                                            u_v_matrix,
-                                            un_v_matrix,
-                                            u_vn_matrix,
-                                            un_vn_matrix);
-
-                     dofs_neighbor.resize (dofs_on_neighbor_child);
-                     neighbor_child->get_dof_indices (dofs_neighbor);
-
-                     for (unsigned int i=0; i<dofs_per_cell; ++i)
-                       for (unsigned int j=0; j<dofs_on_neighbor_child; ++j)
-                         system_matrix.add(dofs[i], dofs_neighbor[j],
-                                           un_v_matrix(i,j));
-
-                     for (unsigned int i=0; i<dofs_on_neighbor_child; ++i)
-                       for (unsigned int j=0; j<dofs_per_cell; ++j)
-                         system_matrix.add(dofs_neighbor[i], dofs[j],
-                                           u_vn_matrix(i,j));
-
-                     for (unsigned int i=0; i<dofs_on_neighbor_child; ++i)
-                       for (unsigned int j=0; j<dofs_on_neighbor_child; ++j)
-                         system_matrix.add(dofs_neighbor[i], dofs_neighbor[j],
-                                           un_vn_matrix(i,j));
-                   }
-               }
-             else
-               {
-                                                  // Case 3, with the
-                                                  // additional rule
-                                                  // a)
-                 if (neighbor->level() == cell->level() &&
-                     neighbor->index() > cell->index()) 
-                   {
-                     const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no);
-                     
-                     const unsigned int quadrature_index = 
-                       std::max (neighbor->active_fe_index (),
-                                  cell->active_fe_index ());
-
-                     un_v_matrix = 0;
-                     u_vn_matrix = 0;
-                     un_vn_matrix = 0;
-                     
-                     fe_v_face_x.reinit (cell, face_no, quadrature_index);
-                     fe_v_face_neighbor_x.reinit (neighbor, neighbor2, quadrature_index);
-
-                     dg.assemble_face_term2(fe_v_face_x.get_present_fe_values (),
-                                            fe_v_face_neighbor_x.get_present_fe_values (),
-                                            u_v_matrix,
-                                            un_v_matrix,
-                                            u_vn_matrix,
-                                            un_vn_matrix);
-
-                     dofs_neighbor.resize (dofs_on_neighbor);
-                     neighbor->get_dof_indices (dofs_neighbor);
-
-                     for (unsigned int i=0; i<dofs_per_cell; ++i)
-                       for (unsigned int j=0; j<dofs_on_neighbor; ++j)
-                         system_matrix.add(dofs[i], dofs_neighbor[j],
-                                           un_v_matrix(i,j));
-
-                     for (unsigned int i=0; i<dofs_on_neighbor; ++i)
-                       for (unsigned int j=0; j<dofs_per_cell; ++j)
-                         system_matrix.add(dofs_neighbor[i], dofs[j],
-                                           u_vn_matrix(i,j));
-
-                     for (unsigned int i=0; i<dofs_on_neighbor; ++i)
-                       for (unsigned int j=0; j<dofs_on_neighbor; ++j)
-                         system_matrix.add(dofs_neighbor[i], dofs_neighbor[j],
-                                           un_vn_matrix(i,j));
-
-                   }
-
-                                                  // Due to rule b)
-                                                  // we do not need
-                                                  // to consider case
-                                                  // 4.
-               }
-           }
-       }
-      
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-       for (unsigned int j=0; j<dofs_per_cell; ++j)
-         system_matrix.add(dofs[i], dofs[j], u_v_matrix(i,j));
-      
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-       right_hand_side(dofs[i]) += cell_vector(i);
-    }
 }
 
 
-                                // @sect3{All the rest}
-                                //
-                                // First, we have to solve the
-                                // discrete system. Since we solve a
-                                // transport equation, the matrix is
-                                // nonsymmetric. Hence, we use a
-                                // GMRES solver.
+
+
+                                // @sect3{Equation data}
+
+                                // As this is one more academic
+                                // example, we'd like to compare
+                                // exact and computed solution
+                                // against each other. For this, we
+                                // need to declare function classes
+                                // representing the exact solution
+                                // (for comparison and for the
+                                // Dirichlet boundary values), as
+                                // well as a class that denotes the
+                                // right hand side of the equation
+                                // (this is simply the Laplace
+                                // operator applied to the exact
+                                // solution we'd like to recover).
                                 //
-                                // For a preconditioner, we use the
-                                // ILU method. Since we already
-                                // sorted the degrees of freedom in
-                                // downwind direction, this should be
-                                // quite efficient. Actually, a block
-                                // Gauss-Seidel method would be an
-                                // exact solver, but it has not been
-                                // implemented for variable block
-                                // sizes yet.
+                                // For this example, let us choose as
+                                // exact solution the function
+                                // $u(x,y)=exp(x+sin(10y+5x^2))$. In more
+                                // than two dimensions, simply repeat
+                                // the sine-factor with <code>y</code>
+                                // replaced by <code>z</code> and so on. Given
+                                // this, the following two classes
+                                // are probably straightforward from
+                                // the previous examples.
                                 //
+                                // As in previous examples, the C++
+                                // language forces us to declare and
+                                // define a constructor to the
+                                // following classes even though they
+                                // are empty. This is due to the fact
+                                // that the base class has no default
+                                // constructor (i.e. one without
+                                // arguments), even though it has a
+                                // constructor which has default
+                                // values for all arguments.
 template <int dim>
-void DGMethod<dim>::solve (Vector<double> &solution) 
+class Solution : public Function<dim>
 {
-  SolverControl solver_control (10000, 1e-12, false, true);
-  SolverGMRES<Vector<double> > solver (solver_control);
-                                  // Initialize the ILU
-                                  // preconditioner. We decide for
-                                  // two additional off diagonals in
-                                  // order to enhance its
-                                  // performance.
-  SparseILU<double>::AdditionalData data(0., 2);
-  SparseILU<double> preconditioner;
-  preconditioner.initialize (system_matrix, data);
-                                  // Then solve the system:
-  solver.solve (system_matrix, solution, right_hand_side,
-             preconditioner);
-}
+  public:
+    Solution () : Function<dim> () {};
+    
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component) const;
+};
 
 
-                                // We refine the grid according to a
-                                // very simple refinement criterion,
-                                // namely an approximation to the
-                                // gradient of the solution. As here
-                                // we consider the DG(1) method
-                                // (i.e. we use piecewise bilinear
-                                // shape functions) we could simply
-                                // compute the gradients on each
-                                // cell. But we do not want to base
-                                // our refinement indicator on the
-                                // gradients on each cell only, but
-                                // want to base them also on jumps of
-                                // the discontinuous solution
-                                // function over faces between
-                                // neighboring cells. The simpliest
-                                // way of doing that is to compute
-                                // approximative gradients by
-                                // difference quotients including the
-                                // cell under consideration and its
-                                // neighbors. This is done by the
-                                // <code>DerivativeApproximation</code> class
-                                // that computes the approximate
-                                // gradients in a way similar to the
-                                // <code>GradientEstimation</code> described
-                                // in Step 9 of this tutorial. In
-                                // fact, the
-                                // <code>DerivativeApproximation</code> class
-                                // was developed following the
-                                // <code>GradientEstimation</code> class of
-                                // Step 9. Relating to the
-                                // discussion in Step 9, here we
-                                // consider $h^{1+d/2}|\nabla_h
-                                // u_h|$. Futhermore we note that we
-                                // do not consider approximate second
-                                // derivatives because solutions to
-                                // the linear advection equation are
-                                // in general not in $H^2$ but in $H^1$
-                                // (to be more precise, in $H^1_\beta$)
-                                // only.
 template <int dim>
-void DGMethod<dim>::refine_grid ()
+double
+Solution<dim>::value (const Point<dim>   &p,
+                     const unsigned int  /*component*/) const
 {
-                                  // The <code>DerivativeApproximation</code>
-                                  // class computes the gradients to
-                                  // float precision. This is
-                                  // sufficient as they are
-                                  // approximate and serve as
-                                  // refinement indicators only.
-  Vector<float> gradient_indicator (triangulation.n_active_cells());
-
-                                  // Now the approximate gradients
-                                  // are computed  
-  DerivativeApproximation::approximate_gradient (dof_handler,
-                                                solution2,
-                                                gradient_indicator);
-  
-                                  // and they are cell-wise scaled by
-                                  // the factor $h^{1+d/2}$
-  typename hp::DoFHandler<dim>::active_cell_iterator
-     cell = dof_handler.begin_active(),
-    endc = dof_handler.end();
-
-  for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no)
-    gradient_indicator(cell_no)*=std::pow(cell->diameter(), 1+1.0*dim/2);
-
-                                  // Finally they serve as refinement
-                                  // indicator.
-  GridRefinement::refine_and_coarsen_fixed_number (triangulation,
-                                                  gradient_indicator,
-                                                  0.3, 0.1);
-
-                                   // Simple heuristic hp-Refinement. 
-                                   // As the indicator marks nonsmooth
-                                   // regions, p-refine all non marked 
-                                   // regions, while the marked
-                                   // regions clearly deserve an 
-                                   // h-refinement.
-  cell = dof_handler.begin_active ();
-  for (; cell!=endc; ++cell)
-    if (!cell->refine_flag_set ()
-       &&
-       (cell->active_fe_index() < fe_collection.size()-1))
-      cell->set_active_fe_index (cell->active_fe_index () + 1);
-
-  triangulation.execute_coarsening_and_refinement ();
+  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;
 }
 
 
-                                // The output of this program
-                                // consists of eps-files of the
-                                // adaptively refined grids and the
-                                // numerical solutions given in
-                                // gnuplot format. This was covered
-                                // in previous examples and will not
-                                // be further commented on.
+
 template <int dim>
-void DGMethod<dim>::output_results (const unsigned int cycle) const
+class RightHandSide : public Function<dim>
 {
-                                  // Write the grid in eps format.
-  std::string filename = "grid-";
-  filename += ('0' + cycle);
-  Assert (cycle < 10, ExcInternalError());
-  
-  filename += ".eps";
-  deallog << "Writing grid to <" << filename << ">..." << std::endl;
-  std::ofstream eps_output (filename.c_str());
-
-  GridOut grid_out;
-  grid_out.write_eps (triangulation, eps_output);
+  public:
+    RightHandSide () : Function<dim> () {};
+    
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component) const;
+};
 
-  Vector<double> active_fe_indices (triangulation.n_active_cells());
-  {
-    unsigned int index = 0;
-    for (typename hp::DoFHandler<dim>::active_cell_iterator
-          cell = dof_handler.begin_active();
-        cell != dof_handler.end(); ++cell, ++index)
-      active_fe_indices(index) = cell->active_fe_index ();
-  }
-  
-  
-                                  // Output of the solution in
-                                  // gnuplot format.
-  filename = "sol-";
-  filename += ('0' + cycle);
-  Assert (cycle < 10, ExcInternalError());
-  
-  //  filename += ".gnuplot";
-  filename += ".gmv";
-  deallog << "Writing solution to <" << filename << ">..."
-           << std::endl << std::endl;
-  std::ofstream gnuplot_output (filename.c_str());
-  
-  DataOut<dim, hp::DoFHandler<dim> > data_out;
-  data_out.attach_dof_handler (dof_handler);
-  data_out.add_data_vector (solution2, "u");
-  data_out.add_data_vector (active_fe_indices, "fe_index");
 
-  data_out.build_patches (4);
-  
-//  data_out.write_gnuplot(gnuplot_output);
-  data_out.write_gmv(gnuplot_output);
+template <int dim>
+double
+RightHandSide<dim>::value (const Point<dim>   &/*p*/,
+                          const unsigned int  /*component*/) const
+{
+  return 1.;
 }
 
 
-                                // The following <code>run</code> function is
-                                // similar to previous examples. The
-                                // only difference is that the
-                                // problem is assembled and solved
-                                // twice on each refinement step;
-                                // first by <code>assemble_system1</code> that
-                                // implements the first version and
-                                // then by <code>assemble_system2</code> that
-                                // implements the second version of
-                                // writing the DG
-                                // discretization. Furthermore the
-                                // time needed by each of the two
-                                // assembling routines is measured.
+
+                                // @sect3{The driver routines}
+
+                                // What is now missing are only the
+                                // functions that actually select the
+                                // various options, and run the
+                                // simulation on successively finer
+                                // grids to monitor the progress as
+                                // the mesh is refined.
+                                //
+                                // This we do in the following
+                                // function: it takes a solver
+                                // object, and a list of
+                                // postprocessing (evaluation)
+                                // objects, and runs them with
+                                // intermittent mesh refinement:
 template <int dim>
-void DGMethod<dim>::run () 
+void
+run_simulation (LaplaceSolver::Base<dim>                     &solver,
+               const std::list<Evaluation::EvaluationBase<dim> *> &postprocessor_list)
 {
-  for (unsigned int cycle=0; cycle<7; ++cycle)
+                                  // We will give an indicator of the
+                                  // step we are presently computing,
+                                  // in order to keep the user
+                                  // informed that something is still
+                                  // happening, and that the program
+                                  // is not in an endless loop. This
+                                  // is the head of this status line:
+  std::cout << "Refinement cycle: ";
+
+                                  // Then start a loop which only
+                                  // terminates once the number of
+                                  // degrees of freedom is larger
+                                  // than 20,000 (you may of course
+                                  // change this limit, if you need
+                                  // more -- or less -- accuracy from
+                                  // your program).
+  for (unsigned int step=0; true; ++step)
     {
-      deallog << "Cycle " << cycle << ':' << std::endl;
-
-      if (cycle == 0)
+                                      // Then give the <code>alive</code>
+                                      // indication for this
+                                      // iteration. Note that the
+                                      // <code>std::flush</code> is needed to
+                                      // have the text actually
+                                      // appear on the screen, rather
+                                      // than only in some buffer
+                                      // that is only flushed the
+                                      // next time we issue an
+                                      // end-line.
+      std::cout << step << " " << std::flush;
+
+                                      // Now solve the problem on the
+                                      // present grid, and run the
+                                      // evaluators on it. The long
+                                      // type name of iterators into
+                                      // the list is a little
+                                      // annoying, but could be
+                                      // shortened by a typedef, if
+                                      // so desired.
+      solver.solve_problem ();
+
+      for (typename std::list<Evaluation::EvaluationBase<dim> *>::const_iterator
+            i = postprocessor_list.begin();
+          i != postprocessor_list.end(); ++i)
        {
-         GridGenerator::hyper_cube (triangulation);
+         (*i)->set_refinement_cycle (step);
+         solver.postprocess (**i);
+       };
 
-         triangulation.refine_global (1);
-       }
-      else
-       refine_grid ();
-      
 
-      deallog << "   Number of active cells:       "
-               << triangulation.n_active_cells()
-               << std::endl;
+                                      // Now check whether more
+                                      // iterations are required, or
+                                      // whether the loop shall be
+                                      // ended:
+      if (solver.n_dofs() < 20000)
+       solver.refine_grid ();
+      else
+       break;
+    };
 
-      setup_system ();
+                                  // Finally end the line in which we
+                                  // displayed status reports:
+  std::cout << std::endl;
+}
 
-      deallog << "   Number of degrees of freedom: "
-               << dof_handler.n_dofs()
-               << std::endl;
 
-                                      // The constructor of the Timer
-                                      // class automatically starts
-                                      // the time measurement.
-      Timer assemble_timer;
-                                      // First assembling routine.
-      assemble_system1 ();
-                                      // The operator () accesses the
-                                      // current time without
-                                      // disturbing the time
-                                      // measurement.
-      deallog << "Time of assemble_system1: "
-               << assemble_timer()
-               << std::endl;
-      solve (solution1);
-
-
-                                      // As preparation for the
-                                      // second assembling routine we
-                                      // reinit the system matrix, the
-                                      // right hand side vector and
-                                      // the Timer object.
-      system_matrix = 0;
-      right_hand_side = 0;
-      assemble_timer.reset();
-
-                                      // We start the Timer,
-      assemble_timer.start();
-                                      // call the second assembling routine
-      assemble_system2 ();
-                                      // and access the current time.
-      deallog << "Time of assemble_system2: "
-               << assemble_timer()
-               << std::endl;
-      solve (solution2);
-
-                                      // To make sure that both
-                                      // versions of the DG method
-                                      // yield the same
-                                      // discretization and hence the
-                                      // same solution we check the
-                                      // two solutions for equality.
-      solution1-=solution2;
-
-      const double difference=solution1.linfty_norm();
-      if (difference>1e-12)
-       deallog << "solution1 and solution2 differ!!" << std::endl;
-      else
-       deallog << "solution1 and solution2 coincide." << std::endl;
-       
-                                      // Finally we perform the
-                                      // output.
-      output_results (cycle);
+void
+create_coarse_grid (Triangulation<2> &coarse_grid)
+{
+  const unsigned int dim = 2;
+  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]);
+  const std::vector<Point<dim> > vertices (&vertices_1[0],
+                                           &vertices_1[n_vertices]);
+  static const int cell_vertices[][GeometryInfo<dim>::vertices_per_cell]
+    = {{0, 1, 5, 6},
+       {1, 2, 6, 7},
+       {2, 3, 7, 8},
+       {3, 4, 8, 9},
+       {5, 6, 10, 11},
+       {8, 9, 12, 13},
+       {10, 11, 14, 15},
+       {12, 13, 17, 18},
+       {14, 15, 19, 20},
+       {15, 16, 20, 21},
+       {16, 17, 21, 22},
+       {17, 18, 22, 23}};
+  const unsigned int
+    n_cells = sizeof(cell_vertices) / sizeof(cell_vertices[0]);
+
+  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;
     }
+
+  coarse_grid.create_triangulation (vertices,
+                                    cells,
+                                    SubCellData());
+  coarse_grid.refine_global (1);
+}
+
+
+                                // The final function is one which
+                                // takes the name of a solver
+                                // (presently "kelly" and "global"
+                                // are allowed), creates a solver
+                                // object out of it using a coarse
+                                // grid (in this case the ubiquitous
+                                // unit square) and a finite element
+                                // object (here the likewise
+                                // ubiquitous bilinear one), and uses
+                                // that solver to ask for the
+                                // solution of the problem on a
+                                // sequence of successively refined
+                                // grids.
+                                //
+                                // The function also sets up two of
+                                // evaluation functions, one
+                                // evaluating the solution at the
+                                // point (0.5,0.5), the other writing
+                                // out the solution to a file.
+template <int dim>
+void solve_problem () 
+{
+  Triangulation<dim> triangulation;
+  create_coarse_grid (triangulation);
+
+  const FE_Q<dim>          fe(1);
+  const QGauss<dim>       quadrature(4);
+  const RightHandSide<dim> rhs_function;
+  const ZeroFunction<dim>      boundary_values;
+
+                                  // Create a solver object of the
+                                  // kind indicated by the argument
+                                  // to this function. If the name is
+                                  // not recognized, throw an
+                                  // exception!
+  LaplaceSolver::RefinementKelly<dim> solver (triangulation, fe,
+                                              quadrature,
+                                              rhs_function,
+                                              boundary_values);
+
+                                  // Next create a table object in
+                                  // which the values of the
+                                  // numerical solution at the point
+                                  // (0.5,0.5) will be stored, and
+                                  // create a respective evaluation
+                                  // object:
+  TableHandler results_table;
+  Evaluation::PointValueEvaluation<dim>
+    postprocessor1 (Point<dim>(0.5,0.5), results_table);
+
+                                  // Also generate an evaluator which
+                                  // writes out the solution:
+  Evaluation::SolutionOutput<dim>
+    postprocessor2 (std::string("solution"),
+                   DataOut<dim>::vtk);
+
+                                  // Take these two evaluation
+                                  // objects and put them in a
+                                  // list...
+  std::list<Evaluation::EvaluationBase<dim> *> postprocessor_list;
+  postprocessor_list.push_back (&postprocessor1);
+  postprocessor_list.push_back (&postprocessor2);
+
+                                  // ... which we can then pass on to
+                                  // the function that actually runs
+                                  // the simulation on successively
+                                  // refined grids:
+  run_simulation (solver, postprocessor_list);
+
+                                  // When this all is done, write out
+                                  // the results of the point
+                                  // evaluations, and finally delete
+                                  // the solver object:
+  results_table.write_text (std::cout);
+
+                                  // And one blank line after all
+                                  // results:
+  std::cout << std::endl;
 }
 
-                                // The following <code>main</code> function is
-                                // similar to previous examples and
-                                // need not to be commented on.
+
+
+                                // There is not much to say about the
+                                // main function. It follows the same
+                                // pattern as in all previous
+                                // examples, with attempts to catch
+                                // thrown exceptions, and displaying
+                                // as much information as possible if
+                                // we should get some. The rest is
+                                // self-explanatory.
 int main () 
 {
   try
     {
-         DGMethod<2> dgmethod;
-         dgmethod.run ();
+      deallog.depth_console (0);
+
+      solve_problem<2> ();
     }
   catch (std::exception &exc)
     {
@@ -1805,8 +2053,6 @@ int main ()
                << 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.