]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Move slightly forward.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 22 Apr 2002 08:33:43 +0000 (08:33 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 22 Apr 2002 08:33:43 +0000 (08:33 +0000)
git-svn-id: https://svn.dealii.org/trunk@5701 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-14/Makefile
deal.II/examples/step-14/step-14.cc

index c02ef0f148efae8f641362e505aadb2afadd4f50..2780e004977524ac0e3a81b80576d25b417f120a 100644 (file)
@@ -148,8 +148,8 @@ Makefile.dep: $(target).cc Makefile \
                            $(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;' \
+       @$(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
index 0283a294e00f3f450def86a037e3c37518f4a681..67ee9296b03002ba1b8c9454c0c4e4fed7d80481 100644 (file)
@@ -34,6 +34,7 @@
 #include <dofs/dof_tools.h>
 #include <fe/fe_q.h>
 #include <fe/fe_values.h>
+#include <fe/fe_tools.h>
 #include <numerics/vectors.h>
 #include <numerics/matrices.h>
 #include <numerics/data_out.h>
@@ -42,6 +43,8 @@
 #include <iostream>
 #include <fstream>
 #include <list>
+#include <algorithm>
+#include <numeric>
 
 #ifdef HAVE_STD_STRINGSTREAM
 #  include <sstream>
@@ -146,66 +149,6 @@ namespace Evaluation
     results_table.add_value ("u(x_0)", point_value);
   };
 
-
-
-
-
-  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)
-  {};
-  
-
-  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);
-  };
-
-
-
-  
 };
 
   
@@ -224,9 +167,15 @@ namespace LaplaceSolver
       virtual void postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const = 0;
       virtual void refine_grid () = 0;
       virtual unsigned int n_dofs () const = 0;
+
+      virtual void set_refinement_cycle (const unsigned int cycle);
+
+      virtual void output_solution () const = 0;
       
     protected:
       const SmartPointer<Triangulation<dim> > triangulation;
+
+      unsigned int refinement_cycle;
   };
 
 
@@ -240,6 +189,15 @@ namespace LaplaceSolver
   template <int dim>
   Base<dim>::~Base () 
   {};
+
+
+
+  template <int dim>
+  void
+  Base<dim>::set_refinement_cycle (const unsigned int cycle)
+  {
+    refinement_cycle = cycle;
+  };
   
 
 
@@ -250,6 +208,7 @@ namespace LaplaceSolver
       Solver (Triangulation<dim>       &triangulation,
              const FiniteElement<dim> &fe,
              const Quadrature<dim>    &quadrature,
+             const Quadrature<dim-1>  &face_quadrature,              
              const Function<dim>      &boundary_values);
       virtual
       ~Solver ();
@@ -269,6 +228,7 @@ namespace LaplaceSolver
     protected:
       const SmartPointer<const FiniteElement<dim> >  fe;
       const SmartPointer<const Quadrature<dim> >     quadrature;
+      const SmartPointer<const Quadrature<dim-1> >   face_quadrature;      
       DoFHandler<dim>                                dof_handler;
       Vector<double>                                 solution;
       const SmartPointer<const Function<dim> >       boundary_values;
@@ -304,11 +264,13 @@ namespace LaplaceSolver
   Solver<dim>::Solver (Triangulation<dim>       &triangulation,
                       const FiniteElement<dim> &fe,
                       const Quadrature<dim>    &quadrature,
+                      const Quadrature<dim-1>  &face_quadrature,
                       const Function<dim>      &boundary_values)
                  :
                  Base<dim> (triangulation),
                  fe (&fe),
                   quadrature (&quadrature),
+                  face_quadrature (&face_quadrature),    
                  dof_handler (triangulation),
                  boundary_values (&boundary_values)
   {};
@@ -506,6 +468,7 @@ namespace LaplaceSolver
       PrimalSolver (Triangulation<dim>       &triangulation,
                    const FiniteElement<dim> &fe,
                    const Quadrature<dim>    &quadrature,
+                   const Quadrature<dim-1>  &face_quadrature,
                    const Function<dim>      &rhs_function,
                    const Function<dim>      &boundary_values);
 
@@ -533,12 +496,14 @@ namespace LaplaceSolver
   PrimalSolver (Triangulation<dim>       &triangulation,
                const FiniteElement<dim> &fe,
                const Quadrature<dim>    &quadrature,
+               const Quadrature<dim-1>  &face_quadrature,
                const Function<dim>      &rhs_function,
                const Function<dim>      &boundary_values)
                  :
                  Base<dim> (triangulation),
                  Solver<dim> (triangulation, fe,
-                              quadrature, boundary_values),
+                              quadrature, face_quadrature,
+                              boundary_values),
                   rhs_function (&rhs_function)
   {};
 
@@ -664,6 +629,7 @@ namespace LaplaceSolver
       RefinementKelly (Triangulation<dim>       &coarse_grid,
                       const FiniteElement<dim> &fe,
                       const Quadrature<dim>    &quadrature,
+                      const Quadrature<dim-1>  &face_quadrature,
                       const Function<dim>      &rhs_function,
                       const Function<dim>      &boundary_values);
 
@@ -677,11 +643,13 @@ namespace LaplaceSolver
   RefinementKelly (Triangulation<dim>       &coarse_grid,
                   const FiniteElement<dim> &fe,
                   const Quadrature<dim>    &quadrature,
+                  const Quadrature<dim-1>  &face_quadrature,
                   const Function<dim>      &rhs_function,
                   const Function<dim>      &boundary_values)
                  :
                  Base<dim> (coarse_grid),
                   PrimalSolver<dim> (coarse_grid, fe, quadrature,
+                                    face_quadrature,
                                     rhs_function, boundary_values)
   {};
 
@@ -722,11 +690,15 @@ 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;
+//    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;
+  double s = 1;
+  for (unsigned int i=0; i<dim; ++i)
+    s *= sin(3.1415926536*p(i));
+  return s;
 };
 
 
@@ -747,24 +719,28 @@ 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;
+//    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);
+//    return -u*(t1+t2+t3);
+  double s = 1;
+  for (unsigned int i=0; i<dim; ++i)
+    s *= sin(3.1415926536*p(i));
+  return dim*3.1415926536*3.1415926536*s;  
 };
 
 
@@ -793,6 +769,11 @@ namespace DualFunctional
       void
       assemble_rhs (const DoFHandler<dim> &dof_handler,
                    Vector<double>        &rhs) const;
+      DeclException1 (ExcEvaluationPointNotFound,
+                     Point<dim>,
+                     << "The evaluation point " << arg1
+                     << " was not found among the vertices of the present grid.");
+
     protected:
       const Point<dim> evaluation_point;
   };
@@ -846,6 +827,7 @@ namespace LaplaceSolver
       DualSolver (Triangulation<dim>       &triangulation,
                  const FiniteElement<dim> &fe,
                  const Quadrature<dim>    &quadrature,
+                 const Quadrature<dim-1>  &face_quadrature,
                  const DualFunctional::DualFunctionalBase<dim> &dual_functional);
 
                                       // XXX
@@ -868,17 +850,21 @@ namespace LaplaceSolver
       static const ZeroFunction<dim> boundary_values;
   };
 
+  template <int dim>
+  const ZeroFunction<dim> DualSolver<dim>::boundary_values;
 
   template <int dim>
   DualSolver<dim>::
   DualSolver (Triangulation<dim>       &triangulation,
              const FiniteElement<dim> &fe,
              const Quadrature<dim>    &quadrature,
+             const Quadrature<dim-1>  &face_quadrature,
              const DualFunctional::DualFunctionalBase<dim> &dual_functional)
                  :
                  Base<dim> (triangulation),
                  Solver<dim> (triangulation, fe,
-                              quadrature, boundary_values),
+                              quadrature, face_quadrature,
+                              boundary_values),
                   dual_functional (&dual_functional)
   {};
 
@@ -925,8 +911,10 @@ namespace LaplaceSolver
   {
     public:
       WeightedResidual (Triangulation<dim>       &coarse_grid,
-                       const FiniteElement<dim> &fe,
+                       const FiniteElement<dim> &primal_fe,
+                       const FiniteElement<dim> &dual_fe,
                        const Quadrature<dim>    &quadrature,
+                       const Quadrature<dim-1>  &face_quadrature,
                        const Function<dim>      &rhs_function,
                        const Function<dim>      &boundary_values,
                        const DualFunctional::DualFunctionalBase<dim> &dual_functional);
@@ -943,25 +931,248 @@ namespace LaplaceSolver
       unsigned int
       n_dofs () const;
 
-      virtual void refine_grid () {};
+      virtual void refine_grid ();
+
+      virtual
+      void
+      output_solution () const;
+
+    private:
+
+                                      /**
+                                       * Declare a data type to
+                                       * represent the mapping between
+                                       * faces and integrated jumps of
+                                       * gradients of each of the
+                                       * solution vectors. See the
+                                       * general documentation of this
+                                       * class for more information.
+                                       */
+      typedef typename std::pair<double,typename DoFHandler<dim>::active_cell_iterator> FaceEntry;
+      typedef typename std::map<typename DoFHandler<dim>::face_iterator,FaceEntry> FaceIntegrals;
+
+
+                                      /**
+                                       * Redeclare an active cell iterator.
+                                       * This is simply for convenience.
+                                       */
+      typedef typename DoFHandler<dim>::active_cell_iterator active_cell_iterator;
+
+                                      /**
+                                     * All data needed by the several
+                                     * functions of the error
+                                     * estimator is gathered in this
+                                     * struct. It is passed as a
+                                     * reference to the separate
+                                     * functions in this class.
+                                     *
+                                     * The reason for invention of
+                                     * this object is two-fold:
+                                     * first, class member data is
+                                     * not possible because no real
+                                     * object is created (all
+                                     * functions are @p{static}),
+                                     * which is a historical
+                                     * reason. Second, if we don't
+                                     * collect the data the various
+                                     * functions need somewhere at a
+                                     * central place, that would mean
+                                     * that the functions would have
+                                     * to allocate them upon
+                                     * need. However, then some
+                                     * variables would be allocated
+                                     * over and over again, which can
+                                     * take a significant amount of
+                                     * time (10-20 per cent) and most
+                                     * importantly, memory allocation
+                                     * requires synchronisation in
+                                     * multithreaded mode. While that
+                                     * is done by the C++ library and
+                                     * has not to be handcoded, it
+                                     * nevertheless seriously damages
+                                     * the ability to efficiently run
+                                     * the functions of this class in
+                                     * parallel, since they are quite
+                                     * often blocked by these
+                                     * synchronisation points.
+                                     *
+                                     * Thus, every thread gets an
+                                     * instance of this class to work
+                                     * with and needs not allocate
+                                     * memory itself, or synchronise
+                                     * with other threads.
+                                     */
+      struct FaceData
+      {
+         FEFaceValues<dim>    fe_face_values_cell;
+         FEFaceValues<dim>    fe_face_values_neighbor;
+         FESubfaceValues<dim> fe_subface_values_cell;
+
+         std::vector<double> jump_residual;
+         std::vector<double> dual_weights;       
+         typename std::vector<Tensor<1,dim> > cell_grads;
+         typename std::vector<Tensor<1,dim> > neighbor_grads;
+         FaceData (const FiniteElement<dim> &dof_handler,
+                   const Quadrature<dim-1>  &face_quadrature);
+      };
+
+      struct CellData
+      {
+         FEValues<dim>    fe_values;
+         const SmartPointer<const Function<dim> > right_hand_side;
+
+         std::vector<double> cell_residual;
+         std::vector<double> rhs_values;         
+         std::vector<double> dual_weights;       
+         typename std::vector<Tensor<2,dim> > cell_grad_grads;
+         CellData (const FiniteElement<dim> &dof_handler,
+                   const Quadrature<dim>    &quadrature,
+                   const Function<dim>      &right_hand_side);
+      };
+      
+
+
+      void estimate_error (Vector<float> &error_indicators) const;
+
+      void estimate_some (const Vector<double> &primal_solution,
+                         const Vector<double> &dual_weights,
+                         const unsigned int    n_threads,
+                         const unsigned int    this_thread,
+                         Vector<float>        &error_indicators,
+                         FaceIntegrals        &face_integrals) const;
+
+      void
+      integrate_over_cell (const active_cell_iterator &cell,
+                          const unsigned int          cell_index,
+                          const Vector<double>       &primal_solution,
+                          const Vector<double>       &dual_weights,
+                          CellData                   &cell_data,
+                          Vector<float>              &error_indicators) const;
+      
+                                      /**
+                                       * Actually do the computation on
+                                       * a face which has no hanging
+                                       * nodes (it is regular), i.e.
+                                       * either on the other side there
+                                       * is nirvana (face is at
+                                       * boundary), or the other side's
+                                       * refinement level is the same
+                                       * as that of this side, then
+                                       * handle the integration of
+                                       * these both cases together.
+                                       *
+                                       * The meaning of the parameters
+                                       * becomes clear when looking at
+                                       * the source code. This function
+                                       * is only externalized from
+                                       * @p{estimate_error} to avoid
+                                       * ending up with a function of
+                                       * 500 lines of code.
+                                       */
+      void
+      integrate_over_regular_face (const active_cell_iterator &cell,
+                                  const unsigned int          face_no,
+                                  const Vector<double>       &primal_solution,
+                                  const Vector<double>       &dual_weights,
+                                  FaceData                   &face_data,
+                                  FaceIntegrals              &face_integrals) const;
+      
+
+                                      /**
+                                       * The same applies as for the
+                                       * function above, except that
+                                       * integration is over face
+                                       * @p{face_no} of @p{cell}, where
+                                       * the respective neighbor is
+                                       * refined, so that the
+                                       * integration is a bit more
+                                       * complex.
+                                       */
+      void
+      integrate_over_irregular_face (const active_cell_iterator &cell,
+                                    const unsigned int          face_no,
+                                    const Vector<double>       &primal_solution,
+                                    const Vector<double>       &dual_weights,
+                                    FaceData                   &face_data,
+                                    FaceIntegrals              &face_integrals) const;
   };
 
 
 
+
+
+  template <int dim>
+  WeightedResidual<dim>::FaceData::
+  FaceData (const FiniteElement<dim> &fe,
+           const Quadrature<dim-1>  &face_quadrature)
+                 :
+                 fe_face_values_cell (fe, face_quadrature,
+                                      update_values        |
+                                      update_gradients     |
+                                      update_JxW_values    |
+                                      update_normal_vectors),
+                 fe_face_values_neighbor (fe, face_quadrature,
+                                          update_values     |
+                                          update_gradients  |
+                                          update_JxW_values |
+                                          update_normal_vectors),
+                 fe_subface_values_cell (fe, face_quadrature,
+                                         update_gradients)
+  {  
+    const unsigned int n_face_q_points
+      = face_quadrature.n_quadrature_points;
+  
+    jump_residual.resize(n_face_q_points);
+    dual_weights.resize(n_face_q_points);    
+    cell_grads.resize(n_face_q_points);
+    neighbor_grads.resize(n_face_q_points);
+  };
+  
+
+
+  template <int dim>
+  WeightedResidual<dim>::CellData::
+  CellData (const FiniteElement<dim> &fe,
+           const Quadrature<dim>    &quadrature,
+           const Function<dim>      &right_hand_side)
+                 :
+                 fe_values (fe, quadrature,
+                            update_values             |
+                            update_second_derivatives |
+                            update_q_points           |
+                            update_JxW_values),
+                 right_hand_side (&right_hand_side)
+  {  
+    const unsigned int n_q_points
+      = quadrature.n_quadrature_points;
+  
+    cell_residual.resize(n_q_points);
+    rhs_values.resize(n_q_points);    
+    dual_weights.resize(n_q_points);    
+    cell_grad_grads.resize(n_q_points);
+  };
+  
+  
+
+
   template <int dim>
   WeightedResidual<dim>::
   WeightedResidual (Triangulation<dim>       &coarse_grid,
-                   const FiniteElement<dim> &fe,
+                   const FiniteElement<dim> &primal_fe,
+                   const FiniteElement<dim> &dual_fe,
                    const Quadrature<dim>    &quadrature,
+                   const Quadrature<dim-1>  &face_quadrature,
                    const Function<dim>      &rhs_function,
                    const Function<dim>      &boundary_values,
                    const DualFunctional::DualFunctionalBase<dim> &dual_functional)
                  :
                  Base<dim> (coarse_grid),
-                  PrimalSolver<dim> (coarse_grid, fe, quadrature,
+                  PrimalSolver<dim> (coarse_grid, primal_fe,
+                                    quadrature, face_quadrature,
                                     rhs_function, boundary_values),
-    DualSolver<dim> (coarse_grid, fe, quadrature,
-                    dual_functional)
+                  DualSolver<dim> (coarse_grid, dual_fe,
+                                  quadrature, face_quadrature,
+                                  dual_functional)
   {};
 
 
@@ -988,10 +1199,688 @@ namespace LaplaceSolver
   {
     return PrimalSolver<dim>::n_dofs();
   };
+
+
+
+  template <int dim>
+  void
+  WeightedResidual<dim>::refine_grid ()
+  {
+    Vector<float> error_indicators (triangulation->n_active_cells());
+    estimate_error (error_indicators);
+    std::cout << "Estimated error="
+             << std::accumulate (error_indicators.begin(),
+                                 error_indicators.end(), 0.)
+             << std::endl;
+    DataOut<dim> data_out;
+    ofstream x("x");
+    Vector<double> xe (error_indicators.begin(),
+                      error_indicators.end());
+    data_out.attach_dof_handler (DualSolver<dim>::dof_handler);
+    data_out.add_data_vector (xe, "e");
+    data_out.build_patches ();
+    data_out.write_gnuplot (x);
+    
+    std::transform (error_indicators.begin(),
+                   error_indicators.end(),
+                   error_indicators.begin(),
+                   &fabs);
+    GridRefinement::refine_and_coarsen_fixed_number (*triangulation,
+                                                    error_indicators,
+                                                    0.3, 0.03);
+    triangulation->execute_coarsening_and_refinement ();
+  };
+  
+
+  template <int dim>
+  void
+  WeightedResidual<dim>::output_solution () const
+  {
+    for (unsigned int solution_type=0; solution_type<2; ++solution_type)
+      {
+       DataOut<dim> data_out;
+
+       switch (solution_type)
+         {
+           case 0:
+                 data_out.attach_dof_handler (PrimalSolver<dim>::dof_handler);
+                 data_out.add_data_vector (PrimalSolver<dim>::solution,
+                                           "primal_solution");
+                 break;
+           case 1:
+                 data_out.attach_dof_handler (DualSolver<dim>::dof_handler);
+                 data_out.add_data_vector (DualSolver<dim>::solution,
+                                           "dual_solution");
+                 break;
+           default:
+                 Assert (false, ExcInternalError());
+         };
+       data_out.build_patches ();
+  
+#ifdef HAVE_STD_STRINGSTREAM
+       std::ostringstream filename;
+#else
+       std::ostrstream filename;
+#endif
+       filename << "solution-"
+                << (solution_type == 0 ?
+                    "primal-" : "dual-")
+                << refinement_cycle
+                << ".gnuplot"
+                << std::ends;
+#ifdef HAVE_STD_STRINGSTREAM
+       std::ofstream out (filename.str().c_str());
+#else
+       std::ofstream out (filename.str());
+#endif
+    
+       data_out.write (out, DataOut<dim>::gnuplot);
+      };
+  };
+
+
+  template <int dim>
+  void
+  WeightedResidual<dim>::
+  estimate_error (Vector<float> &error_indicators) const
+  {
+                                    // The first task in computing
+                                    // the error is to set up vectors
+                                    // that denote the primal
+                                    // solution, and the weights
+                                    // (z-z_h)=(z-I_hz), both in the
+                                    // finite element space for which
+                                    // we have computed the dual
+                                    // solution. For this, we have to
+                                    // interpolate the primal
+                                    // solution to the dual finite
+                                    // element space, and to subtract
+                                    // the interpolation of the
+                                    // computed dual solution to the
+                                    // primal finite element
+                                    // space. Fortunately, the
+                                    // library provides functions for
+                                    // these two actions.
+    Vector<double> primal_solution (DualSolver<dim>::dof_handler.n_dofs());
+    FETools::interpolate (PrimalSolver<dim>::dof_handler,
+                         PrimalSolver<dim>::solution,
+                         DualSolver<dim>::dof_handler,
+                         primal_solution);    
+    Vector<double> dual_weights (DualSolver<dim>::dof_handler.n_dofs());
+    FETools::interpolation_difference (DualSolver<dim>::dof_handler,
+                                      DualSolver<dim>::solution,
+                                      *PrimalSolver<dim>::fe,
+                                      dual_weights);
+    
+    
+                                    // Map of integrals indexed by
+                                    // the corresponding face. In this map
+                                    // we store the integrated jump of the
+                                    // gradient for each face.
+                                    // At the end of the function, we again
+                                    // loop over the cells and collect the
+                                    // contributions of the different faces
+                                    // of the cell.
+                                    // 
+                                    // The initial values for all faces
+                                    // are set to -1e20. It would cost
+                                    // a lot of time to synchronise the
+                                    // initialisation (i.e. the
+                                    // creation of new keys) of the map
+                                    // in multithreaded mode. Negative
+                                    // value indicates that the face
+                                    // has not yet been processed.
+    FaceIntegrals face_integrals;
+    for (active_cell_iterator cell=DualSolver<dim>::dof_handler.begin_active();
+        cell!=DualSolver<dim>::dof_handler.end();
+        ++cell)
+      for (unsigned int face_no=0;
+          face_no<GeometryInfo<dim>::faces_per_cell;
+          ++face_no)
+       face_integrals[cell->face(face_no)].first = -1e20;
+
+                                    // reserve one slot for each cell
+                                    // and set it to zero
+    error_indicators.reinit (DualSolver<dim>::dof_handler
+                            .get_tria().n_active_cells());
+
+
+                                  // all the data needed in the error
+                                  // estimator by each of the threads
+                                  // is gathered in the following
+                                  // stuctures
+                                  //
+                                  // note that if no component mask
+                                  // was given, then treat all
+                                  // components
+    const unsigned int n_threads = multithread_info.n_default_threads;
+  
+                                    // split all cells into threads if
+                                    // multithreading is used and run
+                                    // the whole thing
+    Threads::ThreadManager thread_manager;
+    for (unsigned int i=0; i<n_threads; ++i)
+      Threads::spawn (thread_manager,
+                     Threads::encapsulate (&WeightedResidual<dim>::
+                                           estimate_some)
+                     .collect_args (this,
+                                    primal_solution,
+                                    dual_weights,
+                                    n_threads, i,
+                                    error_indicators,
+                                    face_integrals));
+    thread_manager.wait();
+  
+                                    // finally add up the
+                                    // contributions of the faces for
+                                    // each cell
+  
+    unsigned int present_cell=0;  
+    for (active_cell_iterator cell=DualSolver<dim>::dof_handler.begin_active();
+        cell!=DualSolver<dim>::dof_handler.end();
+        ++cell, ++present_cell)
+      {
+                                        // loop over all faces of this cell
+       for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell;
+            ++face_no)
+         {
+           Assert(face_integrals.find(cell->face(face_no)) != face_integrals.end(),
+                  ExcInternalError());
+           if (true || (face_integrals[cell->face(face_no)].second
+               ==
+                        cell))
+             error_indicators(present_cell)
+               += -0.5*face_integrals[cell->face(face_no)].first;
+           else
+             error_indicators(present_cell)
+               -= -0.5*face_integrals[cell->face(face_no)].first;
+         };
+      };
+  };
+
+
+  template <int dim>
+  void
+  WeightedResidual<dim>::
+  estimate_some (const Vector<double> &primal_solution,
+                const Vector<double> &dual_weights,
+                const unsigned int    n_threads,
+                const unsigned int    this_thread,
+                Vector<float>        &error_indicators,
+                FaceIntegrals        &face_integrals) const
+  {
+    FaceData face_data (*DualSolver<dim>::fe,
+                       *DualSolver<dim>::face_quadrature);    
+    CellData cell_data (*DualSolver<dim>::fe,
+                       *DualSolver<dim>::quadrature,
+                       *PrimalSolver<dim>::rhs_function);
+
+                                  // First calculate the start cell
+                                  // for this thread. We let the
+                                  // different threads run on
+                                  // interleaved cells, i.e. for
+                                  // example if we have 4 threads,
+                                  // then the first thread treates
+                                  // cells 0, 4, 8, etc, while the
+                                  // second threads works on cells 1,
+                                  // 5, 9, and so on. The reason is
+                                  // that it takes vastly more time
+                                  // to work on cells with hanging
+                                  // nodes than on regular cells, but
+                                  // such cells are not evenly
+                                  // distributed across the range of
+                                  // cell iterators, so in order to
+                                  // have the different threads do
+                                  // approximately the same amount of
+                                  // work, we have to let them work
+                                  // interleaved to the effect of a
+                                  // pseudorandom distribution of the
+                                  // `hard' cells to the different
+                                  // threads.
+    active_cell_iterator cell=DualSolver<dim>::dof_handler.begin_active();
+    for (unsigned int t=0;
+        (t<this_thread) && (cell!=DualSolver<dim>::dof_handler.end());
+        ++t, ++cell);
+
+  
+                                    // Then loop over all cells. The
+                                    // check for loop end is done at
+                                    // the end of the loop, along
+                                    // with incrementing the loop
+                                    // index.
+    for (unsigned int cell_index=this_thread; true; )
+      {
+       
+       integrate_over_cell (cell, cell_index,
+                            primal_solution,
+                            dual_weights,
+                            cell_data,
+                            error_indicators);
+       
+                                        // After computing the cell
+                                        // terms, turn to the face
+                                        // terms. For this, loop over
+                                        // all faces of the present
+                                        // cell, and see whether
+                                        // something needs to be
+                                        // computed on it:
+       for (unsigned int face_no=0;
+            face_no<GeometryInfo<dim>::faces_per_cell;
+            ++face_no)
+         {
+                                            // First, if this face is
+                                            // part of the boundary,
+                                            // then there is nothing
+                                            // to do. However, to
+                                            // make things easier
+                                            // when summing up the
+                                            // contributions of the
+                                            // faces of cells, we
+                                            // enter this face into
+                                            // the list of faces with
+                                            // a zero contribution to
+                                            // the error, and also
+                                            // mark the cell on which
+                                            // we computed this
+                                            // value.
+           if (cell->face(face_no)->at_boundary()) 
+             {
+               face_integrals[cell->face(face_no)].first = 0;
+               face_integrals[cell->face(face_no)].second = cell;
+               continue;
+             };
+           
+                                            // Next, note that since
+                                            // we want to compute the
+                                            // jump terms on each
+                                            // face only once,
+                                            // although we access it
+                                            // twice if it is not at
+                                            // the boundary, we have
+                                            // to define some rules
+                                            // who is responsible for
+                                            // computing on a face:
+                                            //
+                                            // First, if the
+                                            // neighboring cell is on
+                                            // the same level as this
+                                            // one, i.e. neither
+                                            // further refined not
+                                            // coarser, then the one
+                                            // with the lower index
+                                            // within this level does
+                                            // the work. In other
+                                            // words: if the other
+                                            // one has a lower index,
+                                            // then skip work on this
+                                            // face:
+           if ((cell->neighbor(face_no)->has_children() == false) &&
+               (cell->neighbor(face_no)->level() == cell->level()) &&
+               (cell->neighbor(face_no)->index() < cell->index()))
+             continue;
+
+                                            // Likewise, we always
+                                            // work from the coarser
+                                            // cell if this and its
+                                            // neighbor differ in
+                                            // refinement. Thus, if
+                                            // the neighboring cell
+                                            // is less refined than
+                                            // the present one, then
+                                            // do nothing since we
+                                            // integrate over the
+                                            // subfaces when we visit
+                                            // the coarse cell.
+           if (cell->at_boundary(face_no) == false)
+             if (cell->neighbor(face_no)->level() < cell->level())
+               continue;         
+
+
+                                            // Now we know that we
+                                            // are in charge here, so
+                                            // actually compute the
+                                            // face jump terms. If
+                                            // the face is a regular
+                                            // one, i.e.  the other
+                                            // side's cell is neither
+                                            // coarser not finer than
+                                            // this cell, then call
+                                            // one function, and if
+                                            // the cell on the other
+                                            // side is further
+                                            // refined, then use
+                                            // another function. Note
+                                            // that the case that the
+                                            // cell on the other side
+                                            // is coarser cannot
+                                            // happen since we have
+                                            // decided above that we
+                                            // handle this case when
+                                            // we pass over that
+                                            // other cell.
+           if (cell->face(face_no)->has_children() == false)
+             integrate_over_regular_face (cell, face_no,
+                                          primal_solution,
+                                          dual_weights,
+                                          face_data,
+                                          face_integrals);       
+           else
+             integrate_over_irregular_face (cell, face_no,
+                                            primal_solution,
+                                            dual_weights,
+                                            face_data,
+                                            face_integrals);
+         };
+
+                                        // After computing the cell
+                                        // contributions and looping
+                                        // over the faces, go to the
+                                        // next cell for this
+                                        // thread. Note again that
+                                        // the cells for each of the
+                                        // threads are
+                                        // interleaved. If we are at
+                                        // the end of our workload,
+                                        // jump out of the loop.
+       for (unsigned int t=0;
+            ((t<n_threads) && (cell!=DualSolver<dim>::dof_handler.end()));
+            ++t, ++cell, ++cell_index);
+       if (cell == DualSolver<dim>::dof_handler.end())
+         break;
+      };    
+  };
+
+
+  template <int dim>
+  void WeightedResidual<dim>::
+  integrate_over_cell (const active_cell_iterator &cell,
+                      const unsigned int          cell_index,
+                      const Vector<double>       &primal_solution,
+                      const Vector<double>       &dual_weights,
+                      CellData                   &cell_data,
+                      Vector<float>              &error_indicators) const
+  {
+    cell_data.fe_values.reinit (cell);
+    cell_data.right_hand_side
+      ->value_list (cell_data.fe_values.get_quadrature_points(),
+                   cell_data.rhs_values);
+    cell_data.fe_values.get_function_2nd_derivatives (primal_solution,
+                                                     cell_data.cell_grad_grads);
+    cell_data.fe_values.get_function_values (dual_weights,
+                                            cell_data.dual_weights);
+    double sum = 0;
+    for (unsigned int p=0; p<cell_data.fe_values.n_quadrature_points; ++p)
+      sum += ((cell_data.rhs_values[p]+trace(cell_data.cell_grad_grads[p])) *
+             cell_data.dual_weights[p] *
+             cell_data.fe_values.JxW (p));
+    error_indicators(cell_index) += sum;
+  };
   
+  template <int dim>
+  void WeightedResidual<dim>::
+  integrate_over_regular_face (const active_cell_iterator &cell,
+                              const unsigned int          face_no,
+                              const Vector<double>       &primal_solution,
+                              const Vector<double>       &dual_weights,
+                              FaceData                   &face_data,
+                              FaceIntegrals              &face_integrals) const
+  {
+    const unsigned int
+      n_q_points = face_data.fe_face_values_cell.n_quadrature_points;
+
+                                    // The first step is to get the
+                                    // values of the gradients at the
+                                    // quadrature points of the
+                                    // finite element field on the
+                                    // present cell. For this,
+                                    // initialize the
+                                    // ``FEFaceValues'' object
+                                    // corresponding to this side of
+                                    // the face, and extract the
+                                    // gradients using that
+                                    // object.
+    face_data.fe_face_values_cell.reinit (cell, face_no);
+    face_data.fe_face_values_cell.get_function_grads (primal_solution,
+                                                     face_data.cell_grads);
+
+                                    // The second step is then to
+                                    // extract the gradients of the
+                                    // finite element solution at the
+                                    // quadrature points on the other
+                                    // side of the face, i.e. from
+                                    // the neighboring cell.
+                                    //
+                                    // For this, do a sanity check
+                                    // before: make sure that the
+                                    // neigbor actually exists (yes,
+                                    // we should not have come here
+                                    // if the neighbor did not exist,
+                                    // but in complicated software
+                                    // there are bugs, so better
+                                    // check this), and if this is
+                                    // not the case throw an error.
+    Assert (cell->neighbor(face_no).state() == IteratorState::valid,
+           ExcInternalError());
+                                    // If we have that, then we need
+                                    // to find out with which face of
+                                    // the neighboring cell we have
+                                    // to work, i.e. the
+                                    // ``home-many''the neighbor the
+                                    // present cell is of the cell
+                                    // behind the present face. For
+                                    // this, there is a function, and
+                                    // we put the result into a
+                                    // variable with the name
+                                    // ``neighbor_neighbor'':
+    const unsigned int
+      neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
+                                    // Then define an abbreviation
+                                    // for the neigbor cell,
+                                    // initialize the
+                                    // ``FEFaceValues'' object on
+                                    // that cell, and extract the
+                                    // gradients on that cell:
+    const active_cell_iterator neighbor = cell->neighbor(face_no);
+    face_data.fe_face_values_neighbor.reinit (neighbor, neighbor_neighbor);      
+    face_data.fe_face_values_neighbor.get_function_grads (primal_solution,
+                                                         face_data.neighbor_grads);
+
+                                    // Now that we have the gradients
+                                    // on this and the neighboring
+                                    // cell, compute the jump
+                                    // residual by multiplying the
+                                    // jump in the gradient with the
+                                    // normal vector:
+    for (unsigned int p=0; p<n_q_points; ++p)
+      face_data.jump_residual[p]
+       = ((face_data.neighbor_grads[p] - face_data.cell_grads[p]) *
+          face_data.fe_face_values_cell.normal_vector(p));
+
+                                    // Next get the dual weights for
+                                    // this face:
+    face_data.fe_face_values_cell.get_function_values (dual_weights,
+                                                      face_data.dual_weights);
+    
+                                    // Finally, we have to compute
+                                    // the sum over jump residuals,
+                                    // dual weights, and quadrature
+                                    // weights, to get the result for
+                                    // this face:
+    double face_integral = 0;
+    for (unsigned int p=0; p<n_q_points; ++p)
+      face_integral += (face_data.jump_residual[p] *
+                       face_data.dual_weights[p]  *
+                       face_data.fe_face_values_cell.JxW(p));
+
+                                    // Double check that the element
+                                    // already exists and that it was
+                                    // not already written to...
+    Assert (face_integrals.find (cell->face(face_no)) != face_integrals.end(),
+           ExcInternalError());
+    Assert (face_integrals[cell->face(face_no)].first == -1e20,
+           ExcInternalError());
+
+                                    // ...then store computed value
+                                    // at assigned location:
+    face_integrals[cell->face(face_no)].first = face_integral;
+    face_integrals[cell->face(face_no)].second = cell;
+  };
+
+
+
+  template <int dim>
+  void WeightedResidual<dim>::
+  integrate_over_irregular_face (const active_cell_iterator &cell,
+                                const unsigned int          face_no,
+                                const Vector<double>       &primal_solution,
+                                const Vector<double>       &dual_weights,
+                                FaceData                   &face_data,
+                                FaceIntegrals              &face_integrals) const
+  {
+                                    // First again two abbreviations,
+                                    // and some consistency checks
+                                    // whether the function is called
+                                    // only on faces for which it is
+                                    // supposed to be called:
+    const unsigned int
+      n_q_points = face_data.fe_face_values_cell.n_quadrature_points;
+
+    const typename DoFHandler<dim>::cell_iterator
+      neighbor = cell->neighbor(face_no);    
+    Assert (neighbor.state() == IteratorState::valid,
+           ExcInternalError());
+    Assert (neighbor->has_children(),
+           ExcInternalError());
+
+                                    // Then find out which neighbor
+                                    // the present cell is of the
+                                    // adjacent cell. Note that we
+                                    // will operator on the children
+                                    // of this adjacent cell, but
+                                    // that their orientation is the
+                                    // same as that of their mother,
+                                    // i.e. the neigbor direction is
+                                    // the same.
+    const unsigned int
+      neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
   
-  template class WeightedResidual<2>;
+                                    // Then simply do everything we
+                                    // did in the previous function
+                                    // for one face for all the
+                                    // sub-faces now:
+    for (unsigned int subface_no=0;
+        subface_no<GeometryInfo<dim>::subfaces_per_face;
+        ++subface_no)
+      {
+                                        // Start with some checks
+                                        // again: get an iterator
+                                        // pointing to the cell
+                                        // behind the present subface
+                                        // and check whether its face
+                                        // is a subface of the one we
+                                        // are considering. If that
+                                        // were not the case, then
+                                        // there would be either a
+                                        // bug in the
+                                        // ``neighbor_neighbor''
+                                        // function called above, or
+                                        // -- worse -- some function
+                                        // in the library did not
+                                        // keep to some underlying
+                                        // assumptions about cells,
+                                        // their children, and their
+                                        // faces. In any case, even
+                                        // though this assertion
+                                        // should not be triggered,
+                                        // it does not harm to be
+                                        // cautious, and in optimized
+                                        // mode computations the
+                                        // assertion will be removed
+                                        // anyway.
+       const active_cell_iterator neighbor_child
+         = neighbor->child(GeometryInfo<dim>::
+                           child_cell_on_face(neighbor_neighbor,
+                                              subface_no));
+       Assert (neighbor_child->face(neighbor_neighbor) ==
+               cell->face(face_no)->child(subface_no),
+               ExcInternalError());
+
+                                        // Now start the work by
+                                        // again getting the gradient
+                                        // of the solution first at
+                                        // this side of the
+                                        // interface,
+       face_data.fe_subface_values_cell.reinit (cell, face_no, subface_no);
+       face_data.fe_subface_values_cell.get_function_grads (primal_solution,
+                                                            face_data.cell_grads);
+                                        // then at the other side,
+       face_data.fe_face_values_neighbor.reinit (neighbor_child,
+                                            neighbor_neighbor);
+       face_data.fe_face_values_neighbor.get_function_grads (primal_solution,
+                                                             face_data.neighbor_grads);
+      
+                                        // and finally building the
+                                        // jump residuals. Since we
+                                        // take the normal vector
+                                        // from the other cell this
+                                        // time, revert the sign of
+                                        // the first term compared to
+                                        // the other function:
+       for (unsigned int p=0; p<n_q_points; ++p)
+         face_data.jump_residual[p]
+            = ((face_data.cell_grads[p] - face_data.neighbor_grads[p]) *
+               face_data.fe_face_values_neighbor.normal_vector(p));
+
+                                        // Then get dual weights:
+       face_data.fe_face_values_neighbor.get_function_values (dual_weights,
+                                                              face_data.dual_weights);
+       
+                                        // At last, sum up the
+                                        // contribution of this
+                                        // sub-face, and set it in
+                                        // the global map:
+       double face_integral = 0;
+       for (unsigned int p=0; p<n_q_points; ++p)
+         face_integral += (face_data.jump_residual[p] *
+                           face_data.dual_weights[p] *
+                           face_data.fe_face_values_neighbor.JxW(p));
+       face_integrals[neighbor_child->face(neighbor_neighbor)].first
+         = face_integral;
+       face_integrals[neighbor_child->face(neighbor_neighbor)].second
+         = cell;
+      };
 
+                                    // Once the contributions of all
+                                    // sub-faces are computed, loop
+                                    // over all sub-faces to collect
+                                    // and store them with the mother
+                                    // face for simple use when later
+                                    // collecting the error terms of
+                                    // cells. Again make safety
+                                    // checks that the entries for
+                                    // the sub-faces have been
+                                    // computed and do not carry an
+                                    // invalid value.
+    double sum = 0;
+    typename DoFHandler<dim>::face_iterator face = cell->face(face_no);
+    for (unsigned int subface_no=0;
+        subface_no<GeometryInfo<dim>::subfaces_per_face;
+        ++subface_no) 
+      {
+       Assert (face_integrals.find(face->child(subface_no)) !=
+               face_integrals.end(),
+               ExcInternalError());
+       Assert (face_integrals[face->child(subface_no)].first != -1e20,
+               ExcInternalError());
+      
+       sum += face_integrals[face->child(subface_no)].first;
+      };
+                                    // Finally store the value with
+                                    // the parent face.
+    face_integrals[face].first = sum;
+    face_integrals[face].second = cell;    
+  };
+  
 };
 
 
@@ -1006,9 +1895,13 @@ run_simulation (LaplaceSolver::Base<dim>                     &solver,
 
   for (unsigned int step=0; true; ++step)
     {
-      std::cout << step << " " << std::flush;
+      std::cout << step << " Solving "
+               << solver.n_dofs()
+               << std::endl;
 
+      solver.set_refinement_cycle (step);
       solver.solve_problem ();
+      solver.output_solution ();
 
       for (typename std::list<Evaluation::EvaluationBase<dim> *>::const_iterator
             i = postprocessor_list.begin();
@@ -1019,7 +1912,7 @@ run_simulation (LaplaceSolver::Base<dim>                     &solver,
        };
 
 
-      if (solver.n_dofs() < 20000)
+      if (solver.n_dofs() < 5000)
        solver.refine_grid ();
       else
        break;
@@ -1030,6 +1923,52 @@ run_simulation (LaplaceSolver::Base<dim>                     &solver,
 
 
 
+void
+create_triangulation (Triangulation<2> &tria)
+{
+  const Point<2>
+    vertices[16] = {  Point<2> (-1.,   -1.),
+                     Point<2> (-1./3, -1.),
+                     Point<2> (+1./3, -1.),
+                     Point<2> (+1,    -1.),
+                     Point<2> (-1.,   -1./3.),
+                     Point<2> (-1./3, -1./3.),
+                     Point<2> (+1./3, -1./3.),
+                     Point<2> (+1,    -1./3.),
+                     Point<2> (-1.,   1./3.),
+                     Point<2> (-1./3, 1./3.),
+                     Point<2> (+1./3, 1./3.),
+                     Point<2> (+1,    1./3.),
+                     Point<2> (-1.,   1.),
+                     Point<2> (-1./3, 1.),
+                     Point<2> (+1./3, 1.),
+                     Point<2> (+1,    1.)    };
+  
+  const int cell_vertices[8][4] = {{0, 1, 5, 4},
+                                  {1, 2, 6, 5},
+                                  {2, 3, 7, 6},
+                                  {4, 5, 9, 8},
+                                  {6, 7, 11, 10},
+                                  {8,9,13,12},
+                                  {9,10,14,13},
+                                  {10,11,15,14}};
+
+  std::vector<CellData<2> > cells (8, CellData<2>());
+  
+  for (unsigned int i=0; i<8; ++i) 
+    {
+      for (unsigned int j=0; j<4; ++j)
+       cells[i].vertices[j] = cell_vertices[i][j];
+      cells[i].material_id = 0;
+    };
+  
+  tria.create_triangulation (std::vector<Point<2> >(&vertices[0], &vertices[16]),
+                            cells,
+                            SubCellData());       // no boundary information
+};
+
+
+
 template <int dim>
 void solve_problem (const std::string &solver_name) 
 {
@@ -1039,38 +1978,36 @@ void solve_problem (const std::string &solver_name)
            << std::string (header.size(), '-') << std::endl;
 
   Triangulation<dim> triangulation;
+//  create_triangulation (triangulation);
   GridGenerator::hyper_cube (triangulation, -1, 1);
   triangulation.refine_global (2);
-  const FE_Q<dim>          fe(1);
+  const FE_Q<dim>          primal_fe(1);
+  const FE_Q<dim>          dual_fe(2);
   const QGauss4<dim>       quadrature;
+  const QGauss4<dim-1>     face_quadrature;  
   const RightHandSide<dim> rhs_function;
   const Solution<dim>      boundary_values;
 
+  const Point<dim> evaluation_point(0.5,0.5);
+  const DualFunctional::PointValueEvaluation<dim>
+    dual_functional (evaluation_point);
+  
   LaplaceSolver::Base<dim> * solver = 0;
-  if (solver_name == "global")
-    solver = new LaplaceSolver::RefinementGlobal<dim> (triangulation, fe,
-                                                      quadrature,
-                                                      rhs_function,
-                                                      boundary_values);
-  else if (solver_name == "kelly")
-    solver = new LaplaceSolver::RefinementKelly<dim> (triangulation, fe,
-                                                     quadrature,
-                                                     rhs_function,
-                                                     boundary_values);
-  else
-    AssertThrow (false, ExcNotImplemented());
+  solver = new LaplaceSolver::WeightedResidual<dim> (triangulation,
+                                                    primal_fe,
+                                                    dual_fe,
+                                                    quadrature,
+                                                    face_quadrature,
+                                                    rhs_function,
+                                                    boundary_values,
+                                                    dual_functional);
 
   TableHandler results_table;
   Evaluation::PointValueEvaluation<dim>
     postprocessor1 (Point<dim>(0.5,0.5), results_table);
 
-  Evaluation::SolutionOutput<dim>
-    postprocessor2 (std::string("solution-")+solver_name,
-                   DataOut<dim>::gnuplot);
-
   std::list<Evaluation::EvaluationBase<dim> *> postprocessor_list;
   postprocessor_list.push_back (&postprocessor1);
-  postprocessor_list.push_back (&postprocessor2);
 
   run_simulation (*solver, postprocessor_list);
 
@@ -1089,7 +2026,7 @@ int main ()
       deallog.depth_console (0);
 
       solve_problem<2> ("global");
-      solve_problem<2> ("kelly");      
+//      solve_problem<2> ("kelly");      
     }
   catch (std::exception &exc)
     {

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.