]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Remove old benchmarks that haven't been looked at in years. 436/head
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 13 Jan 2015 20:16:41 +0000 (14:16 -0600)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 13 Jan 2015 20:16:41 +0000 (14:16 -0600)
They almost certainly also don't build any more.

tests/benchmarks/old/Makefile [deleted file]
tests/benchmarks/old/benchmark.cc [deleted file]
tests/benchmarks/old/dof_handler_size.cc [deleted file]
tests/benchmarks/old/dof_handler_timing.cc [deleted file]
tests/benchmarks/old/dof_handler_timing_01.cc [deleted file]
tests/benchmarks/old/patches/library/2007-05-30 [deleted file]
tests/benchmarks/old/performance.cc [deleted file]
tests/benchmarks/old/performance.pl [deleted file]
tests/benchmarks/old/spec2006-447.dealII.cc [deleted file]

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

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

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.