]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add step-13 to have in CVS. Not yet finished, though.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 28 Mar 2002 08:06:46 +0000 (08:06 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 28 Mar 2002 08:06:46 +0000 (08:06 +0000)
git-svn-id: https://svn.dealii.org/trunk@5628 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-13/Makefile [new file with mode: 0644]
deal.II/examples/step-13/step-13.cc [new file with mode: 0644]

diff --git a/deal.II/examples/step-13/Makefile b/deal.II/examples/step-13/Makefile
new file mode 100644 (file)
index 0000000..c02ef0f
--- /dev/null
@@ -0,0 +1,159 @@
+# $Id$
+
+
+# For the small projects Makefile, you basically need to fill in only
+# four fields.
+#
+# The first is the name of the application. It is assumed that the
+# application name is the same as the base file name of the single C++
+# file from which the application is generated.
+target = $(basename $(shell echo step-*.cc))
+
+# The second field determines whether you want to run your program in
+# debug or optimized mode. The latter is significantly faster, but no
+# run-time checking of parameters and internal states is performed, so
+# you should set this value to `on' while you develop your program,
+# and to `off' when running production computations.
+debug-mode = on
+
+
+# As third field, we need to give the path to the top-level deal.II
+# directory. You need to adjust this to your needs. Since this path is
+# probably the most often needed one in the Makefile internals, it is
+# designated by a single-character variable, since that can be
+# reference using $D only, i.e. without the parentheses that are
+# required for most other parameters, as e.g. in $(target).
+D = ../../
+
+
+# The last field specifies the names of data and other files that
+# shall be deleted when calling `make clean'. Object and backup files,
+# executables and the like are removed anyway. Here, we give a list of
+# files in the various output formats that deal.II supports.
+clean-up-files = *gmv *gnuplot *gpl *eps *pov
+
+
+
+
+#
+#
+# Usually, you will not need to change something beyond this point.
+#
+#
+# The next statement tell the `make' program where to find the
+# deal.II top level directory and to include the file with the global
+# settings
+include $D/common/Make.global_options
+
+
+# Since the whole project consists of only one file, we need not
+# consider difficult dependencies. We only have to declare the
+# libraries which we want to link to the object file, and there need
+# to be two sets of libraries: one for the debug mode version of the
+# application and one for the optimized mode. Here we have selected
+# the versions for 2d. Note that the order in which the libraries are
+# given here is important and that your applications won't link
+# properly if they are given in another order.
+#
+# You may need to augment the lists of libraries when compiling your
+# program for other dimensions, or when using third party libraries
+libs.g   = $(lib-deal2-2d.g) \
+          $(lib-lac.g)      \
+           $(lib-base.g)
+libs.o   = $(lib-deal2-2d.o) \
+          $(lib-lac.o)      \
+           $(lib-base.o)
+
+
+# We now use the variable defined above which switch between debug and
+# optimized mode to select the correct compiler flags and the set of
+# libraries to link with. Included in the list of libraries is the
+# name of the object file which we will produce from the single C++
+# file. Note that by default we use the extension .go for object files
+# compiled in debug mode and .o for object files in optimized mode.
+ifeq ($(debug-mode),on)
+  libraries = $(target).go $(libs.g)
+  flags     = $(CXXFLAGS.g)
+else
+  libraries = $(target).o $(libs.o)
+  flags     = $(CXXFLAGS.o)
+endif
+
+
+# Now comes the first production rule: how to link the single object
+# file produced from the single C++ file into the executable. Since
+# this is the first rule in the Makefile, it is the one `make' selects
+# if you call it without arguments.
+$(target) : $(libraries)
+       @echo ============================ Linking $@
+       @$(CXX) -o $@ $^ $(LIBS) $(LDFLAGS)
+
+
+# To make running the application somewhat independent of the actual
+# program name, we usually declare a rule `run' which simply runs the
+# program. You can then run it by typing `make run'. This is also
+# useful if you want to call the executable with arguments which do
+# not change frequently. You may then want to add them to the
+# following rule:
+run: $(target)
+       @echo ============================ Running $<
+       @./$(target)
+
+
+# As a last rule to the `make' program, we define what to do when
+# cleaning up a directory. This usually involves deleting object files
+# and other automatically created files such as the executable itself,
+# backup files, and data files. Since the latter are not usually quite
+# diverse, you needed to declare them at the top of this file.
+clean:
+       -rm -f *.o *.go *~ Makefile.dep $(target) $(clean-up-files)
+
+
+# Since we have not yet stated how to make an object file from a C++
+# file, we should do so now. Since the many flags passed to the
+# compiler are usually not of much interest, we suppress the actual
+# command line using the `at' sign in the first column of the rules
+# and write the string indicating what we do instead.
+%.go : %.cc
+       @echo ==============debug========= $(<F)
+       @$(CXX) $(CXXFLAGS.g) -c $< -o $@
+%.o : %.cc
+       @echo ==============optimized===== $(<F)
+       @$(CXX) $(CXXFLAGS.o) -c $< -o $@
+
+
+# The following statement tells make that the rules `run' and `clean'
+# are not expected to produce files of the same name as Makefile rules
+# usually do.
+.PHONY: run clean
+
+
+# Finally there is a rule which you normally need not care much about:
+# since the executable depends on some include files from the library,
+# besides the C++ application file of course, it is necessary to
+# re-generate the executable when one of the files it depends on has
+# changed. The following rule to created a dependency file
+# `Makefile.dep', which `make' uses to determine when to regenerate
+# the executable. This file is automagically remade whenever needed,
+# i.e. whenever one of the cc-(include-path-base)/baseh-files changed. Make detects whether
+# to remake this file upon inclusion at the bottom of this file.
+#
+# The dependency file is created using a perl script.  Since the
+# script prefixes the output names by `lib(include-path-base)/baseo' or `lib(include-path-base)/basego' (it was
+# written for the sublibraries' Makefile), we have to strip that again
+# since object files are placed in the present directory for this
+# application. All these things are made in the next rule:
+Makefile.dep: $(target).cc Makefile \
+              $(shell echo $(include-path-base)/base/*.h    \
+                           $(include-path-lac)/lac/*.h      \
+                           $(include-path-deal2)/*/*.h)
+       @echo ============================ Remaking Makefile
+       @perl $D/common/scripts/make_dependencies.pl  $(INCLUDE) $(target).cc \
+               | perl -pi -e 's!lib/g?o/!!g;' \
+               > Makefile.dep
+
+# To make the dependencies known to `make', we finally have to include
+# them:
+include Makefile.dep
+
+
diff --git a/deal.II/examples/step-13/step-13.cc b/deal.II/examples/step-13/step-13.cc
new file mode 100644 (file)
index 0000000..70bc958
--- /dev/null
@@ -0,0 +1,1434 @@
+/* $Id$ */
+/* Author: Wolfgang Bangerth, University of Heidelberg, 2001, 2002 */
+
+/*    $Id$       */
+/*    Version: $Name$                                          */
+/*                                                                */
+/*    Copyright (C) 2001, 2002 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.                        */
+
+
+                                // 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
+                                // ``base'' - ``lac'' - ``grid'' -
+                                // ``dofs'' - ``fe'' - ``numerics''
+                                // (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/solver_cg.h>
+#include <lac/vector_memory.h>
+#include <lac/precondition.h>
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#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 <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>
+
+                                // Now for the C++ standard headers:
+#include <fstream>
+#include <list>
+
+                                // Just as in the step-5 example
+                                // program (see there for a lengthy
+                                // discussion of the subject), we
+                                // have to work around some
+                                // historical confusion with the
+                                // files declaring the stringstream
+                                // classes:
+#ifdef HAVE_STD_STRINGSTREAM
+#  include <sstream>
+#else
+#  include <strstream>
+#endif
+
+
+                                // @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
+                                // ``DoFHandler'' 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.
+                                //
+                                // From an abstract point of view, we
+                                // declare an abstract base class
+                                // that provides and evaluation
+                                // operator ``operator()'' which will
+                                // do the evaluation of the solution
+                                // (whatever derived classes might
+                                // consider an ``evaluation''). 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 ``operator()'' a
+                                // ``functor'' 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 solver
+                                // 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
+{
+
+                                  // Now for the abstract base class
+                                  // of evaluation classes: its main
+                                  // purpose is to declare a pure
+                                  // virtual function ``operator()''
+                                  // taking a ``DoFHandler'' 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 ();
+
+      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;
+  };
+
+
+                                  // @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
+                                  // ``operator()''. 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 ``DeclExceptionN''
+                                  // 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)
+  {};
+  
+
+
+                                  // 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
+                                            // happend to be located
+                                            // at the vertices, but
+                                            // are not associated
+                                            // with the vertices bur
+                                            // 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
+                                            // ``Assert
+                                            // (dof_handler.get_fe().dofs_per_vertex
+                                            // > 0,
+                                            // ExcNotImplemented())'',
+                                            // 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
+                                            // ``cell->vertex_dof_index(vertex,0)''
+                                            // would fail if we asked
+                                            // it to give us the DoF
+                                            // index of a vertex if
+                                            // there were none.
+                                            //
+                                            // We briefly note that
+                                            // this restriction on
+                                            // the allowed finite
+                                            // elements should be
+                                            // stated in the class
+                                            // documentation.
+           
+           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 computation were useless
+                                    // anyway. So make sure through
+                                    // the ``AssertThrow'' 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 ``throw''
+                                    // 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 ``catch'' clause in
+                                    // your main function (as this
+                                    // program has), you will catch
+                                    // all exceptions that are not
+                                    // caught somewhere 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));
+
+                                    // 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);
+  };
+
+
+
+
+                                  // @sect4{Generating output}
+
+                                  // A different, maybe slightly odd
+                                  // kind of ``evaluation'' of a
+                                  // solution is to output it to a
+                                  // file in a graphical
+                                  // format. Since in the evaluation
+                                  // functions we are given a
+                                  // ``DoFHandler'' 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 ``EvaluationBase'' base
+                                  // class, its main interface is the
+                                  // ``operator()''
+                                  // 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
+                                  // ``DataOutInterface'' class
+                                  // (which is a base class of
+                                  // ``DataOut'' through which we
+                                  // will access its fields) provides
+                                  // an enumeration field
+                                  // ``OutputFormat'', 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 ``ucd'',
+                                  // ``gnuplot'', ``povray'',
+                                  // ``eps'', ``gmv'', and ``vtk'',
+                                  // 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
+                                  // ``write'' function, which the
+                                  // branches to the
+                                  // ``write_gnuplot'',
+                                  // ``write_ucd'', 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
+                                  // ``DataOut::default_suffix''
+                                  // 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
+                                  // ``DataOut::write'' 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.
+                                  //
+                                  // The somewhat complicated use of
+                                  // the stringstream class,
+                                  // involving support from the
+                                  // preprocessor, as already
+                                  // explained in the step-5 example
+                                  // program.
+  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 ();
+  
+#ifdef HAVE_STD_STRINGSTREAM
+    std::ostringstream filename;
+#else
+    std::ostrstream filename;
+#endif
+    filename << output_name_base << "-"
+            << refinement_cycle
+            << data_out.default_suffix (output_format)
+            << std::ends;
+#ifdef HAVE_STD_STRINGSTREAM
+    std::ofstream out (filename.str().c_str());
+#else
+    std::ofstream out (filename.str());
+#endif
+    
+    data_out.write (out, output_format);
+  };
+
+
+                                  // In practical applications, one
+                                  // would add here a list of other
+                                  // possible evaluation classes,
+                                  // representing quantities of
+                                  // interest that one is interested
+                                  // in. For this examples, that much
+                                  // shall be sufficient, so we close
+                                  // the namespace.
+};
+
+  
+                                // @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.
+                                //
+                                // Since we have discussed Laplace
+                                // solvers already in considerable
+                                // detail in previous examples, the
+                                // 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 examples
+                                // 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
+{
+                                  // @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, solver a problem,
+                                  // postprocess the solution with a
+                                  // list of evaluation objects, and
+                                  // refine the grid,
+                                  // respectively. None of these
+                                  // functions actually does
+                                  // something itself.
+                                  //
+                                  // 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
+                                  // even derived from entirely
+                                  // abstract classes ``Object'',
+                                  // 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 then can 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 ``SmartPointer'' 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
+                                  // ``refine_grid'' function.
+  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;
+
+    protected:
+      const SmartPointer<Triangulation<dim> > triangulation;
+  };
+
+
+                                  // 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)
+  {};
+
+
+  template <int dim>
+  Base<dim>::~Base () 
+  {};
+  
+
+                                  // @sect3{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
+                                  // ``solve_problem'' and
+                                  // ``postprocess'' functions
+                                  // declared in the base class. It
+                                  // does not, however, implement the
+                                  // ``refine_grid'' method, as mesh
+                                  // refinement will be implemented
+                                  // in a number of derived classes.
+                                  //
+                                  // It also declares a new abstract
+                                  // virtual function,
+                                  // ``assemble_rhs'', 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.
+  template <int dim>
+  class Solver : public virtual Base<dim>
+  {
+    public:
+      Solver (Triangulation<dim>       &triangulation,
+             const FiniteElement<dim> &fe,
+             const Function<dim>      &boundary_values);
+      virtual ~Solver ();
+      virtual void solve_problem ();
+      virtual void postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const;
+
+    protected:
+      const SmartPointer<const FiniteElement<dim> >  fe;
+      DoFHandler<dim>                                dof_handler;
+      Vector<double>                                 solution;
+      const SmartPointer<const Function<dim> >       boundary_values;
+      
+      virtual void assemble_rhs (Vector<double> &rhs) const = 0;
+    
+    private:
+      struct LinearSystem
+      {
+         LinearSystem (const DoFHandler<dim> &dof_handler);
+
+         void solve (Vector<double> &solution) const;
+       
+         ConstraintMatrix     hanging_node_constraints;
+         SparsityPattern      sparsity_pattern;
+         SparseMatrix<double> matrix;
+         Vector<double>       rhs;
+      };
+
+      void assemble_linear_system (LinearSystem &linear_system);
+
+      void assemble_matrix (LinearSystem                                &linear_system,
+                           const DoFHandler<dim>::active_cell_iterator &begin_cell,
+                           const DoFHandler<dim>::active_cell_iterator &end_cell,
+                           Threads::ThreadMutex                        &mutex) const      ;
+  };
+
+
+
+
+  template <int dim>
+  Solver<dim>::Solver (Triangulation<dim>       &triangulation,
+                      const FiniteElement<dim> &fe,
+                      const Function<dim>      &boundary_values)
+                 :
+                 Base<dim> (triangulation),
+                 fe (&fe),
+                 dof_handler (triangulation),
+                 boundary_values (&boundary_values)
+  {};
+
+
+  template <int dim>
+  Solver<dim>::~Solver () 
+  {
+    dof_handler.clear ();
+  };
+
+
+
+  template <int dim>
+  void
+  Solver<dim>::solve_problem ()
+  {
+    dof_handler.distribute_dofs (*fe);
+    solution.reinit (dof_handler.n_dofs());
+
+    LinearSystem linear_system (dof_handler);
+    assemble_linear_system (linear_system);
+    linear_system.solve (solution);
+  };
+
+
+
+  template <int dim>
+  Solver<dim>::LinearSystem::
+  LinearSystem (const DoFHandler<dim> &dof_handler)
+  {
+    hanging_node_constraints.clear ();
+    DoFTools::make_hanging_node_constraints (dof_handler,
+                                            hanging_node_constraints);
+    hanging_node_constraints.close ();
+
+    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);
+
+    hanging_node_constraints.condense (sparsity_pattern);
+
+    sparsity_pattern.compress();
+
+    matrix.reinit (sparsity_pattern);
+    rhs.reinit (dof_handler.n_dofs());
+  };
+
+
+
+  template <int dim>
+  void
+  Solver<dim>::assemble_linear_system (LinearSystem &linear_system)
+  {
+    typedef typename DoFHandler<dim>::active_cell_iterator active_cell_iterator;
+
+    const unsigned int n_threads = multithread_info.n_default_threads;
+    std::vector<std::pair<active_cell_iterator,active_cell_iterator> >
+      thread_ranges 
+      = Threads::split_range<active_cell_iterator> (dof_handler.begin_active (),
+                                                   dof_handler.end (),
+                                                   n_threads);
+    Threads::ThreadMutex mutex;
+    Threads::ThreadManager thread_manager;
+    for (unsigned int thread=0; thread<n_threads; ++thread)
+      Threads::spawn (thread_manager,
+                     Threads::encapsulate(&Solver<dim>::assemble_matrix)
+                     .collect_args (this,
+                                    linear_system,
+                                    thread_ranges[thread].first,
+                                    thread_ranges[thread].second,
+                                    mutex));
+    assemble_rhs (linear_system.rhs);
+    linear_system.hanging_node_constraints.condense (linear_system.rhs);
+
+    thread_manager.wait ();
+    linear_system.hanging_node_constraints.condense (linear_system.matrix);
+
+    std::map<unsigned int,double> boundary_value_map;
+    VectorTools::interpolate_boundary_values (dof_handler,
+                                             0,
+                                             *boundary_values,
+                                             boundary_value_map);
+    MatrixTools::apply_boundary_values (boundary_value_map,
+                                       linear_system.matrix,
+                                       solution,
+                                       linear_system.rhs);
+
+  };
+
+  
+  template <int dim>
+  void
+  Solver<dim>::assemble_matrix (LinearSystem                                &linear_system,
+                               const DoFHandler<dim>::active_cell_iterator &begin_cell,
+                               const DoFHandler<dim>::active_cell_iterator &end_cell,
+                               Threads::ThreadMutex                        &mutex) const
+  {
+                                    //TODO: adaptive
+    QGauss4<dim>  quadrature_formula;
+
+    FEValues<dim> fe_values (*fe, quadrature_formula, 
+                            UpdateFlags(update_gradients |
+                                        update_JxW_values));
+
+    const unsigned int   dofs_per_cell = fe->dofs_per_cell;
+    const unsigned int   n_q_points    = quadrature_formula.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.clear ();
+
+       fe_values.reinit (cell);
+       const std::vector<std::vector<Tensor<1,dim> > >
+         & shape_grads  = fe_values.get_shape_grads();
+       const std::vector<double>
+         & JxW_values   = fe_values.get_JxW_values();
+
+       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) += (shape_grads[i][q_point] *
+                                  shape_grads[j][q_point] *
+                                  JxW_values[q_point]);
+
+
+       cell->get_dof_indices (local_dof_indices);
+       mutex.acquire ();
+       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));
+       mutex.release ();
+      };
+  };
+
+
+
+  template <int dim>
+  void
+  Solver<dim>::LinearSystem::solve (Vector<double> &solution) const
+  {
+    SolverControl           solver_control (1000, 1e-12);
+    PrimitiveVectorMemory<> vector_memory;
+    SolverCG<>              cg (solver_control, vector_memory);
+
+    PreconditionSSOR<> preconditioner;
+    preconditioner.initialize(matrix, 1.2);
+
+    cg.solve (matrix, solution, rhs, preconditioner);
+
+    hanging_node_constraints.distribute (solution);
+  };
+
+
+
+  template <int dim>
+  void
+  Solver<dim>::
+  postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const
+  {
+    postprocessor (dof_handler, solution);
+  };
+  
+
+//----------------------------------------------------------    
+
+  template <int dim>
+  class PrimalSolver : public Solver<dim>
+  {
+    public:
+      PrimalSolver (Triangulation<dim>       &triangulation,
+                   const FiniteElement<dim> &fe,
+                   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;
+  };
+
+
+
+  template <int dim>
+  PrimalSolver<dim>::
+  PrimalSolver (Triangulation<dim>       &triangulation,
+               const FiniteElement<dim> &fe,
+               const Function<dim>      &rhs_function,
+               const Function<dim>      &boundary_values)
+                 :
+                 Base<dim> (triangulation),
+                 Solver<dim> (triangulation, fe, boundary_values),
+                  rhs_function (&rhs_function)
+  {};
+
+
+
+  template <int dim>
+  void
+  PrimalSolver<dim>::
+  assemble_rhs (Vector<double> &rhs) const 
+  {
+                                    //TODO: adaptive
+    QGauss4<dim>  quadrature_formula;
+
+    FEValues<dim> fe_values (*fe, quadrature_formula, 
+                            UpdateFlags(update_values    |
+                                        update_q_points  |
+                                        update_JxW_values));
+
+    const unsigned int   dofs_per_cell = fe->dofs_per_cell;
+    const unsigned int   n_q_points    = quadrature_formula.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 = dof_handler.begin_active(),
+      endc = dof_handler.end();
+    for (; cell!=endc; ++cell)
+      {
+       cell_rhs.clear ();
+
+       fe_values.reinit (cell);
+       const FullMatrix<double> 
+         & shape_values = fe_values.get_shape_values();
+       const std::vector<double>
+         & JxW_values   = fe_values.get_JxW_values();
+       const std::vector<Point<dim> >
+         & q_points     = fe_values.get_quadrature_points();
+
+       rhs_function->value_list (q_points, rhs_values);
+      
+       for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+         for (unsigned int i=0; i<dofs_per_cell; ++i)
+           cell_rhs(i) += (shape_values (i,q_point) *
+                           rhs_values[q_point] *
+                           JxW_values[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);
+      };
+  };
+
+
+//----------------------------------------------------------    
+
+  template <int dim>
+  class RefinementKelly : public PrimalSolver<dim>
+  {
+    public:
+      RefinementKelly (Triangulation<dim>       &coarse_grid,
+                      const FiniteElement<dim> &fe,
+                      const Function<dim>      &rhs_function,
+                      const Function<dim>      &boundary_values);
+
+      virtual void refine_grid ();
+  };
+
+
+
+  template <int dim>
+  RefinementKelly<dim>::
+  RefinementKelly (Triangulation<dim>       &coarse_grid,
+                  const FiniteElement<dim> &fe,
+                  const Function<dim>      &rhs_function,
+                  const Function<dim>      &boundary_values)
+                 :
+                 Base<dim> (coarse_grid),
+    PrimalSolver<dim> (coarse_grid, fe, rhs_function, boundary_values)
+  {};
+
+
+
+  template <int dim>
+  void
+  RefinementKelly<dim>::refine_grid ()
+  {
+    Vector<float> estimated_error_per_cell (triangulation->n_active_cells());
+    KellyErrorEstimator<dim>::estimate (dof_handler,
+                                       QGauss3<dim-1>(),
+                                       typename FunctionMap<dim>::type(),
+                                       solution,
+                                       estimated_error_per_cell);
+    GridRefinement::refine_and_coarsen_fixed_number (*triangulation,
+                                                    estimated_error_per_cell,
+                                                    0.3, 0.03);
+    triangulation->execute_coarsening_and_refinement ();
+  };
+
+
+
+//----------------------------------------------------------    
+
+  template <int dim>
+  class RefinementGlobal : public PrimalSolver<dim>
+  {
+    public:
+      RefinementGlobal (Triangulation<dim>       &coarse_grid,
+                       const FiniteElement<dim> &fe,
+                       const Function<dim>      &rhs_function,
+                       const Function<dim>      &boundary_values);
+
+      virtual void refine_grid ();
+  };
+
+
+
+  template <int dim>
+  RefinementGlobal<dim>::
+  RefinementGlobal (Triangulation<dim>       &coarse_grid,
+                   const FiniteElement<dim> &fe,
+                   const Function<dim>      &rhs_function,
+                   const Function<dim>      &boundary_values)
+                 :
+                 Base<dim> (coarse_grid),
+    PrimalSolver<dim> (fe, rhs_function, boundary_values)
+  {};
+
+
+
+  template <int dim>
+  void
+  RefinementGlobal<dim>::refine_grid ()
+  {
+    triangulation->refine_global (1);
+  };
+};
+
+
+
+
+                                // @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 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 ``y''
+                                // replaced by ''z'' and so on. Given
+                                // this, the following two classes
+                                // are probably straightforward from
+                                // the previous examples.
+template <int dim>
+class Solution : public Function<dim>
+{
+  public:
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component) const;
+};
+
+
+template <int dim>
+double
+Solution<dim>::value (const Point<dim>   &p,
+                     const unsigned int  /*component*/) const
+{
+  double q = p(0);
+  for (unsigned int i=1; i<dim; ++i)
+    q += sin(10*p(i)+5*p(0)*p(0));
+  const double exponential = exp(q);
+  return exponential;
+};
+
+
+
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+  public:
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component) const;
+};
+
+
+template <int dim>
+double
+RightHandSide<dim>::value (const Point<dim>   &p,
+                          const unsigned int  /*component*/) const
+{
+  double q = p(0);
+  for (unsigned int i=1; i<dim; ++i)
+    q += sin(10*p(i)+5*p(0)*p(0));
+  const double u = exp(q);
+  double t1 = 1,
+        t2 = 0,
+        t3 = 0;
+  for (unsigned int i=1; i<dim; ++i)
+    {
+      t1 += cos(10*p(i)+5*p(0)*p(0)) * 10 * p(0);
+      t2 += 10*cos(10*p(i)+5*p(0)*p(0)) -
+           100*sin(10*p(i)+5*p(0)*p(0)) * p(0)*p(0);
+      t3 += 100*cos(10*p(i)+5*p(0)*p(0))*cos(10*p(i)+5*p(0)*p(0)) -
+           100*sin(10*p(i)+5*p(0)*p(0));
+    };
+  t1 = t1*t1;
+  
+  return -u*(t1+t2+t3);
+};
+
+
+
+                                // @sect3{The driver routines}
+
+
+template <int dim>
+void
+run_simulation (LaplaceSolver::Base<dim>                     &solver,
+               const std::list<Evaluation::EvaluationBase<dim> *> &postprocessor_list)
+{
+  const unsigned int max_steps = 10;
+  for (unsigned int step=0; step<max_steps; ++step)
+    {
+      std::cout << "Refinement cycle " << step << std::endl;
+      
+      solver.solve_problem ();
+
+      for (typename std::list<Evaluation::EvaluationBase<dim> *>::const_iterator
+            i = postprocessor_list.begin();
+          i != postprocessor_list.end(); ++i)
+       {
+         (*i)->set_refinement_cycle (step);
+         solver.postprocess (**i);
+       };
+
+      if (step!=max_steps-1)
+       solver.refine_grid ();
+    };
+};
+
+
+template <int dim>
+void solve_problem_kelly () 
+{      
+  Triangulation<dim> triangulation;
+  GridGenerator::hyper_cube (triangulation, -1, 1);
+  triangulation.refine_global (2);
+  FE_Q<dim> fe(1);
+  const RightHandSide<dim> rhs_function;
+  const Solution<dim>      boundary_values;
+      
+  LaplaceSolver::RefinementKelly<dim> kelly (triangulation, fe,
+                                            rhs_function,
+                                            boundary_values);
+  TableHandler results_table;
+  
+  Evaluation::PointValueEvaluation<dim>
+    postprocessor1 (Point<dim>(.5,.5), results_table);
+  Evaluation::SolutionOutput<dim>
+    postprocessor2 ("solution-kelly", DataOut<dim>::gnuplot);
+  std::list<Evaluation::EvaluationBase<dim> *> postprocessor_list;
+  postprocessor_list.push_back (&postprocessor1);
+  postprocessor_list.push_back (&postprocessor2);
+  
+  run_simulation (kelly, postprocessor_list);
+
+  results_table.write_text (std::cout);
+};
+
+
+    
+int main () 
+{
+  try
+    {
+      deallog.depth_console (0);
+
+      solve_problem_kelly<2> ();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+               << exc.what() << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      return 1;
+    }
+  catch (...) 
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      return 1;
+    };
+
+  return 0;
+};

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.