From: Wolfgang Bangerth <bangerth@math.tamu.edu>
Date: Wed, 3 May 2006 15:04:28 +0000 (+0000)
Subject: Add more tests that all fail for various reasons at the moment.
X-Git-Tag: v8.0.0~11787
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=437d7fae5d64df7ee7a5132fe0017470a612e49c;p=dealii.git

Add more tests that all fail for various reasons at the moment.


git-svn-id: https://svn.dealii.org/trunk@13018 0785d39b-7218-0410-832d-ea1e28bc413d
---

diff --git a/tests/hp/step-10.cc b/tests/hp/step-10.cc
new file mode 100644
index 0000000000..733b981986
--- /dev/null
+++ b/tests/hp/step-10.cc
@@ -0,0 +1,236 @@
+//----------------------------  step-10.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2005, 2006 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  step-10.cc  ---------------------------
+
+
+// a hp-ified version of step-10
+
+
+#include <base/logstream.h>
+#include <fstream>
+std::ofstream logfile("step-10/output");
+
+
+#include <base/quadrature_lib.h>
+#include <base/convergence_table.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_boundary_lib.h>
+#include <grid/tria.h>
+#include <grid/grid_out.h>
+#include <dofs/hp_dof_handler.h>
+#include <dofs/dof_accessor.h>
+#include <fe/fe_q.h>
+#include <fe/hp_fe_values.h>
+
+#include <fe/mapping_q.h>
+
+#include <iostream>
+#include <fstream>
+#include <cmath>
+
+const long double pi = 3.141592653589793238462643;
+
+
+
+template <int dim>
+void gnuplot_output()
+{
+  deallog << "Output of grids into gnuplot files:" << std::endl
+	    << "===================================" << std::endl;
+
+  Triangulation<dim> triangulation;
+  GridGenerator::hyper_ball (triangulation);
+  static const HyperBallBoundary<dim> boundary;
+  triangulation.set_boundary (0, boundary);
+
+  for (unsigned int refinement=0; refinement<2;
+       ++refinement, triangulation.refine_global(1))
+    {
+      deallog << "Refinement level: " << refinement << std::endl;
+
+      std::string filename_base = "ball";
+      filename_base += '0'+refinement;
+
+      for (unsigned int degree=1; degree<4; ++degree)
+	{
+	  deallog << "Degree = " << degree << std::endl;
+
+	  const MappingQ<dim> mapping (degree);
+
+
+	  GridOut grid_out;
+	  GridOutFlags::Gnuplot gnuplot_flags(false, 30);
+	  grid_out.set_flags(gnuplot_flags);
+  
+	  std::string filename = filename_base+"_mapping_q";
+	  filename += ('0'+degree);
+	  filename += ".dat";
+	  std::ofstream gnuplot_file (filename.c_str());
+
+	  grid_out.write_gnuplot (triangulation, gnuplot_file, &mapping);
+	}
+      deallog << std::endl;
+    }
+}
+
+template <int dim>
+void compute_pi_by_area ()
+{
+  deallog << "Computation of Pi by the area:" << std::endl
+	    << "==============================" << std::endl;
+
+  const hp::QCollection<dim> quadrature(QGauss<dim>(4));
+
+  for (unsigned int degree=1; degree<5; ++degree)
+    {
+      deallog << "Degree = " << degree << std::endl;
+
+      Triangulation<dim> triangulation;
+      GridGenerator::hyper_ball (triangulation);
+  
+      static const HyperBallBoundary<dim> boundary;
+      triangulation.set_boundary (0, boundary);
+
+      MappingQ<dim> m(degree);
+      const hp::MappingCollection<dim> mapping (m);
+
+      const hp::FECollection<dim>     dummy_fe (FE_Q<dim>(1));
+
+      hp::DoFHandler<dim> dof_handler (triangulation);
+
+      hp::FEValues<dim> x_fe_values (mapping, dummy_fe, quadrature,
+                               update_JxW_values);
+
+      ConvergenceTable table;
+
+      for (unsigned int refinement=0; refinement<6;
+	   ++refinement, triangulation.refine_global (1))
+	{
+	  table.add_value("cells", triangulation.n_active_cells());
+
+	  dof_handler.distribute_dofs (dummy_fe);
+
+	  long double area = 0;
+
+	  typename hp::DoFHandler<dim>::active_cell_iterator
+	    cell = dof_handler.begin_active(),
+	    endc = dof_handler.end();
+	  for (; cell!=endc; ++cell)
+	    {
+	      x_fe_values.reinit (cell);
+	      const FEValues<dim> &fe_values = x_fe_values.get_present_fe_values();
+	      for (unsigned int i=0; i<fe_values.n_quadrature_points; ++i)
+		area += fe_values.JxW (i);
+	    };
+
+	  table.add_value("eval.pi", static_cast<double> (area));
+	  table.add_value("error",   static_cast<double> (std::fabs(area-pi)));
+	};
+
+      table.omit_column_from_convergence_rate_evaluation("cells");
+      table.omit_column_from_convergence_rate_evaluation("eval.pi");
+      table.evaluate_all_convergence_rates(ConvergenceTable::reduction_rate_log2);
+
+      table.set_precision("eval.pi", 16);
+      table.set_scientific("error", true);
+
+      table.write_text(deallog.get_file_stream());
+
+      deallog << std::endl;
+    };
+}
+
+
+template <int dim>
+void compute_pi_by_perimeter ()
+{
+  deallog << "Computation of Pi by the perimeter:" << std::endl
+	    << "===================================" << std::endl;
+
+  const hp::QCollection<dim-1> quadrature(4);
+
+  for (unsigned int degree=1; degree<5; ++degree)
+    {
+      deallog << "Degree = " << degree << std::endl;
+      Triangulation<dim> triangulation;
+      GridGenerator::hyper_ball (triangulation);
+  
+      static const HyperBallBoundary<dim> boundary;
+      triangulation.set_boundary (0, boundary);
+
+      const MappingQ<dim> m (degree);
+      const hp::MappingCollection<dim> mapping(m);
+      const hp::FECollection<dim>     fe (FE_Q<dim>(1));
+
+      hp::DoFHandler<dim> dof_handler (triangulation);
+
+      hp::FEFaceValues<dim> x_fe_face_values (mapping, fe, quadrature,
+                                        update_JxW_values);
+      ConvergenceTable table;
+
+      for (unsigned int refinement=0; refinement<6;
+	   ++refinement, triangulation.refine_global (1))
+	{
+	  table.add_value("cells", triangulation.n_active_cells());
+
+	  dof_handler.distribute_dofs (fe);
+
+	  typename hp::DoFHandler<dim>::active_cell_iterator
+	    cell = dof_handler.begin_active(),
+	    endc = dof_handler.end();
+	  long double perimeter = 0;
+	  for (; cell!=endc; ++cell)
+	    for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+	      if (cell->face(face_no)->at_boundary())
+		{
+		  x_fe_face_values.reinit (cell, face_no);
+		  const FEFaceValues<dim> &fe_face_values
+		    = x_fe_face_values.get_present_fe_values ();
+		  
+		  for (unsigned int i=0; i<fe_face_values.n_quadrature_points; ++i)
+		    perimeter += fe_face_values.JxW (i);
+		};
+	  table.add_value("eval.pi", static_cast<double> (perimeter/2.));
+	  table.add_value("error",   static_cast<double> (std::fabs(perimeter/2.-pi)));
+	};
+
+      table.omit_column_from_convergence_rate_evaluation("cells");
+      table.omit_column_from_convergence_rate_evaluation("eval.pi");
+      table.evaluate_all_convergence_rates(ConvergenceTable::reduction_rate_log2);
+
+      table.set_precision("eval.pi", 16);
+      table.set_scientific("error", true);
+
+      table.write_text(deallog.get_file_stream());
+
+      deallog << std::endl;
+    };
+}
+
+
+int main () 
+{
+  logfile.precision(2);
+  
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);  
+
+  gnuplot_output<2>();
+
+  compute_pi_by_area<2> ();
+  compute_pi_by_perimeter<2> ();
+  
+  return 0;
+}
diff --git a/tests/hp/step-6.cc b/tests/hp/step-6.cc
new file mode 100644
index 0000000000..2731346c4d
--- /dev/null
+++ b/tests/hp/step-6.cc
@@ -0,0 +1,415 @@
+//----------------------------  step-6.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2005, 2006 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  step-6.cc  ---------------------------
+
+
+// a hp-ified version of step-6
+
+
+#include <base/logstream.h>
+#include <fstream>
+std::ofstream logfile("step-6/output");
+
+
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+#include <base/logstream.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <lac/sparse_matrix.h>
+#include <lac/solver_cg.h>
+#include <lac/precondition.h>
+#include <grid/tria.h>
+#include <dofs/hp_dof_handler.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_boundary_lib.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <fe/hp_fe_values.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+
+#include <fstream>
+#include <iostream>
+
+#include <fe/fe_q.h>
+#include <grid/grid_out.h>
+
+
+#include <dofs/dof_constraints.h>
+
+#include <grid/grid_refinement.h>
+
+#include <numerics/error_estimator.h>
+
+
+
+
+template <int dim>
+class LaplaceProblem 
+{
+  public:
+    LaplaceProblem ();
+    ~LaplaceProblem ();
+
+    void run ();
+    
+  private:
+    void setup_system ();
+    void assemble_system ();
+    void solve ();
+    void refine_grid ();
+    void output_results (const unsigned int cycle) const;
+
+    Triangulation<dim>   triangulation;
+
+    hp::DoFHandler<dim>      dof_handler;
+    hp::FECollection<dim>            fe;
+
+    ConstraintMatrix     hanging_node_constraints;
+
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+
+    Vector<double>       solution;
+    Vector<double>       system_rhs;
+};
+
+
+
+
+template <int dim>
+class Coefficient : public Function<dim> 
+{
+  public:
+    Coefficient () : Function<dim>() {};
+    
+    virtual double value (const Point<dim>   &p,
+			  const unsigned int  component = 0) const;
+    
+    virtual void value_list (const std::vector<Point<dim> > &points,
+			     std::vector<double>            &values,
+			     const unsigned int              component = 0) const;
+};
+
+
+
+template <int dim>
+double Coefficient<dim>::value (const Point<dim> &p,
+				const unsigned int) const 
+{
+  if (p.square() < 0.5*0.5)
+    return 20;
+  else
+    return 1;
+}
+
+
+
+template <int dim>
+void Coefficient<dim>::value_list (const std::vector<Point<dim> > &points,
+				   std::vector<double>            &values,
+				   const unsigned int              component) const 
+{
+  const unsigned int n_points = points.size();
+
+  Assert (values.size() == n_points, 
+	  ExcDimensionMismatch (values.size(), n_points));
+  
+  Assert (component == 0, 
+	  ExcIndexRange (component, 0, 1));
+  
+  for (unsigned int i=0; i<n_points; ++i)
+    {
+      if (points[i].square() < 0.5*0.5)
+	values[i] = 20;
+      else
+	values[i] = 1;
+    }
+}
+
+
+
+
+template <int dim>
+LaplaceProblem<dim>::LaplaceProblem () :
+		dof_handler (triangulation),
+                fe (FE_Q<dim>(2))
+{}
+
+
+
+template <int dim>
+LaplaceProblem<dim>::~LaplaceProblem () 
+{
+  dof_handler.clear ();
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::setup_system ()
+{
+  dof_handler.distribute_dofs (fe);
+
+  sparsity_pattern.reinit (dof_handler.n_dofs(),
+			   dof_handler.n_dofs(),
+			   dof_handler.max_couplings_between_dofs());
+  DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
+
+  solution.reinit (dof_handler.n_dofs());
+  system_rhs.reinit (dof_handler.n_dofs());
+
+  
+  hanging_node_constraints.clear ();
+  DoFTools::make_hanging_node_constraints (dof_handler,
+					   hanging_node_constraints);
+
+  hanging_node_constraints.close ();
+
+  hanging_node_constraints.condense (sparsity_pattern);
+
+  sparsity_pattern.compress();
+
+  system_matrix.reinit (sparsity_pattern);
+}
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_system () 
+{  
+  const hp::QCollection<dim>  quadrature_formula(3);
+
+  hp::FEValues<dim> x_fe_values (fe, quadrature_formula, 
+			   update_values    |  update_gradients |
+			   update_q_points  |  update_JxW_values);
+
+  const unsigned int   dofs_per_cell = fe[0].dofs_per_cell;
+  const unsigned int   n_q_points    = quadrature_formula[0].n_quadrature_points;
+
+  FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+  Vector<double>       cell_rhs (dofs_per_cell);
+
+  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+
+  const Coefficient<dim> coefficient;
+  std::vector<double>    coefficient_values (n_q_points);
+
+  typename hp::DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
+    {
+      cell_matrix = 0;
+      cell_rhs = 0;
+
+      x_fe_values.reinit (cell);
+      const FEValues<dim> &fe_values = x_fe_values.get_present_fe_values();
+      
+      coefficient.value_list (fe_values.get_quadrature_points(),
+			      coefficient_values);
+      
+      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+	for (unsigned int i=0; i<dofs_per_cell; ++i)
+	  {
+	    for (unsigned int j=0; j<dofs_per_cell; ++j)
+	      cell_matrix(i,j) += (coefficient_values[q_point] *
+				   fe_values.shape_grad(i,q_point) *
+				   fe_values.shape_grad(j,q_point) *
+				   fe_values.JxW(q_point));
+
+	    cell_rhs(i) += (fe_values.shape_value(i,q_point) *
+			    1.0 *
+			    fe_values.JxW(q_point));
+	  }
+
+      cell->get_dof_indices (local_dof_indices);
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+	{
+	  for (unsigned int j=0; j<dofs_per_cell; ++j)
+	    system_matrix.add (local_dof_indices[i],
+			       local_dof_indices[j],
+			       cell_matrix(i,j));
+	  
+	  system_rhs(local_dof_indices[i]) += cell_rhs(i);
+	}
+    }
+
+  hanging_node_constraints.condense (system_matrix);
+  hanging_node_constraints.condense (system_rhs);
+
+  std::map<unsigned int,double> boundary_values;
+  VectorTools::interpolate_boundary_values (dof_handler,
+					    0,
+					    ZeroFunction<dim>(),
+					    boundary_values);
+  MatrixTools::apply_boundary_values (boundary_values,
+				      system_matrix,
+				      solution,
+				      system_rhs);
+}
+
+
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::solve () 
+{
+  SolverControl           solver_control (1000, 1e-12);
+  SolverCG<>              cg (solver_control);
+
+  PreconditionSSOR<> preconditioner;
+  preconditioner.initialize(system_matrix, 1.2);
+
+  cg.solve (system_matrix, solution, system_rhs,
+	    preconditioner);
+
+  hanging_node_constraints.distribute (solution);
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::refine_grid ()
+{
+  Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+
+  KellyErrorEstimator<dim>::estimate (dof_handler,
+				      QGauss<dim-1>(3),
+				      typename FunctionMap<dim>::type(),
+				      solution,
+				      estimated_error_per_cell);
+
+  GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+						   estimated_error_per_cell,
+						   0.3, 0.03);
+
+  triangulation.execute_coarsening_and_refinement ();
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
+{
+  Assert (cycle < 10, ExcNotImplemented());
+
+  std::string filename = "grid-";
+  filename += ('0' + cycle);
+  filename += ".eps";
+  
+  std::ofstream output (filename.c_str());
+
+  GridOut grid_out;
+  grid_out.write_eps (triangulation, output);
+}
+
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::run () 
+{
+  for (unsigned int cycle=0; cycle<8; ++cycle)
+    {
+      deallog << "Cycle " << cycle << ':' << std::endl;
+
+      if (cycle == 0)
+	{
+	  GridGenerator::hyper_ball (triangulation);
+
+	  static const HyperBallBoundary<dim> boundary;
+	  triangulation.set_boundary (0, boundary);
+
+	  triangulation.refine_global (1);
+	}
+      else
+	refine_grid ();
+      
+
+      deallog << "   Number of active cells:       "
+		<< triangulation.n_active_cells()
+		<< std::endl;
+
+      setup_system ();
+
+      deallog << "   Number of degrees of freedom: "
+		<< dof_handler.n_dofs()
+		<< std::endl;
+      
+      assemble_system ();
+      solve ();
+      output_results (cycle);
+    }
+
+  DataOutBase::EpsFlags eps_flags;
+  eps_flags.z_scaling = 4;
+  
+  DataOut<dim,hp::DoFHandler<dim> > data_out;
+  data_out.set_flags (eps_flags);
+
+  data_out.attach_dof_handler (dof_handler);
+  data_out.add_data_vector (solution, "solution");
+  data_out.build_patches ();
+  
+  std::ofstream output ("final-solution.eps");
+  data_out.write_eps (output);
+}
+
+
+
+int main () 
+{
+  logfile.precision(2);
+  
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);  
+
+  try
+    {
+      deallog.depth_console (0);
+
+      LaplaceProblem<2> laplace_problem_2d;
+      laplace_problem_2d.run ();
+    }
+  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;
+}
diff --git a/tests/hp/step-7.cc b/tests/hp/step-7.cc
new file mode 100644
index 0000000000..e053f7c411
--- /dev/null
+++ b/tests/hp/step-7.cc
@@ -0,0 +1,742 @@
+//----------------------------  step-7.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2005, 2006 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  step-7.cc  ---------------------------
+
+
+// a hp-ified version of step-7
+
+
+#include <base/logstream.h>
+#include <fstream>
+std::ofstream logfile("step-7/output");
+
+
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+#include <base/logstream.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <lac/sparse_matrix.h>
+#include <lac/solver_cg.h>
+#include <lac/precondition.h>
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/grid_refinement.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_boundary_lib.h>
+#include <dofs/hp_dof_handler.h>
+#include <dofs/dof_constraints.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <fe/fe_q.h>
+#include <numerics/matrices.h>
+#include <numerics/error_estimator.h>
+#include <numerics/data_out.h>
+
+#include <dofs/dof_renumbering.h>
+#include <base/smartpointer.h>
+#include <numerics/vectors.h>
+#include <base/convergence_table.h>
+#include <fe/hp_fe_values.h>
+
+#include <typeinfo>
+#include <fstream>
+#include <iostream>
+
+
+
+template <int dim>
+class SolutionBase 
+{
+  protected:
+    static const unsigned int n_source_centers = 3;    
+    static const Point<dim>   source_centers[n_source_centers];
+    static const double       width;
+};
+
+
+template <>
+const Point<1>
+SolutionBase<1>::source_centers[SolutionBase<1>::n_source_centers]
+= { Point<1>(-1.0 / 3.0), 
+    Point<1>(0.0), 
+    Point<1>(+1.0 / 3.0)   };
+
+template <>
+const Point<2>
+SolutionBase<2>::source_centers[SolutionBase<2>::n_source_centers]
+= { Point<2>(-0.5, +0.5), 
+    Point<2>(-0.5, -0.5), 
+    Point<2>(+0.5, -0.5)   };
+
+template <int dim>
+const double SolutionBase<dim>::width = 1./3.;
+
+
+
+template <int dim>
+class Solution : public Function<dim>,
+		 protected SolutionBase<dim>
+{
+  public:
+    Solution () : Function<dim>() {};
+    
+    virtual double value (const Point<dim>   &p,
+			  const unsigned int  component = 0) const;
+    
+    virtual Tensor<1,dim> gradient (const Point<dim>   &p,
+				    const unsigned int  component = 0) const;
+};
+
+
+template <int dim>
+double Solution<dim>::value (const Point<dim>   &p,
+			     const unsigned int) const
+{
+  double return_value = 0;
+  for (unsigned int i=0; i<this->n_source_centers; ++i)
+    {
+      const Point<dim> x_minus_xi = p - this->source_centers[i];
+      return_value += std::exp(-x_minus_xi.square() /
+			       (this->width * this->width));
+    }
+  
+  return return_value;
+}
+
+
+template <int dim>
+Tensor<1,dim> Solution<dim>::gradient (const Point<dim>   &p,
+				       const unsigned int) const
+{
+  Tensor<1,dim> return_value;
+
+  for (unsigned int i=0; i<this->n_source_centers; ++i)
+    {
+      const Point<dim> x_minus_xi = p - this->source_centers[i];
+      
+      return_value += (-2 / (this->width * this->width) *
+		       std::exp(-x_minus_xi.square() /
+				(this->width * this->width)) *
+		       x_minus_xi);
+    }
+  
+  return return_value;
+}
+
+
+
+template <int dim>
+class RightHandSide : public Function<dim>,
+		      protected SolutionBase<dim>
+{
+  public:
+    RightHandSide () : Function<dim>() {};
+    
+    virtual double value (const Point<dim>   &p,
+			  const unsigned int  component = 0) const;
+};
+
+
+template <int dim>
+double RightHandSide<dim>::value (const Point<dim>   &p,
+				  const unsigned int) const
+{
+  double return_value = 0;
+  for (unsigned int i=0; i<this->n_source_centers; ++i)
+    {
+      const Point<dim> x_minus_xi = p - this->source_centers[i];
+      
+      return_value += ((2*dim - 4*x_minus_xi.square()/
+			(this->width * this->width)) / 
+		       (this->width * this->width) *
+		       std::exp(-x_minus_xi.square() /
+				(this->width * this->width)));
+      return_value += std::exp(-x_minus_xi.square() /
+			       (this->width * this->width));
+    }
+  
+  return return_value;
+}
+
+
+
+template <int dim>
+class HelmholtzProblem 
+{
+  public:
+    enum RefinementMode {
+	  global_refinement, adaptive_refinement
+    };
+    
+    HelmholtzProblem (const hp::FECollection<dim> &fe,
+                      const RefinementMode      refinement_mode);
+
+    ~HelmholtzProblem ();
+
+    void run ();
+    
+  private:
+    void setup_system ();
+    void assemble_system ();
+    void solve ();
+    void refine_grid ();
+    void process_solution (const unsigned int cycle);
+
+    Triangulation<dim>                      triangulation;
+    hp::DoFHandler<dim>                         dof_handler;
+
+    SmartPointer<const hp::FECollection<dim> > fe;
+
+    ConstraintMatrix                        hanging_node_constraints;
+
+    SparsityPattern                         sparsity_pattern;
+    SparseMatrix<double>                    system_matrix;
+
+    Vector<double>                          solution;
+    Vector<double>                          system_rhs;
+
+    const RefinementMode                    refinement_mode;
+
+    ConvergenceTable                        convergence_table;
+};
+
+
+
+
+template <int dim>
+HelmholtzProblem<dim>::HelmholtzProblem (const hp::FECollection<dim> &fe,
+                                         const RefinementMode refinement_mode) :
+		dof_handler (triangulation),
+		fe (&fe),
+		refinement_mode (refinement_mode)
+{}
+
+
+
+template <int dim>
+HelmholtzProblem<dim>::~HelmholtzProblem () 
+{
+  dof_handler.clear ();
+}
+
+
+
+template <int dim>
+void HelmholtzProblem<dim>::setup_system ()
+{
+  dof_handler.distribute_dofs (*fe);
+  DoFRenumbering::Cuthill_McKee (dof_handler);
+
+  hanging_node_constraints.clear ();
+  DoFTools::make_hanging_node_constraints (dof_handler,
+					   hanging_node_constraints);
+  hanging_node_constraints.close ();
+
+  sparsity_pattern.reinit (dof_handler.n_dofs(),
+			   dof_handler.n_dofs(),
+			   dof_handler.max_couplings_between_dofs());
+  DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
+  hanging_node_constraints.condense (sparsity_pattern);
+  sparsity_pattern.compress();
+
+  system_matrix.reinit (sparsity_pattern);
+
+  solution.reinit (dof_handler.n_dofs());
+  system_rhs.reinit (dof_handler.n_dofs());
+}
+
+
+
+template <int dim>
+void HelmholtzProblem<dim>::assemble_system () 
+{  
+  hp::QCollection<dim>   quadrature_formula(3);
+  hp::QCollection<dim-1> face_quadrature_formula(3);
+
+  const unsigned int n_q_points    = quadrature_formula[0].n_quadrature_points;
+  const unsigned int n_face_q_points = face_quadrature_formula[0].n_quadrature_points;
+
+  const unsigned int dofs_per_cell = (*fe)[0].dofs_per_cell;
+
+  FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+  Vector<double>       cell_rhs (dofs_per_cell);
+
+  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+  
+  hp::FEValues<dim>  x_fe_values (*fe, quadrature_formula, 
+			    update_values   | update_gradients |
+                            update_q_points | update_JxW_values);
+
+  hp::FEFaceValues<dim> x_fe_face_values (*fe, face_quadrature_formula, 
+				    update_values         | update_q_points  |
+                                    update_normal_vectors | update_JxW_values);
+
+  const RightHandSide<dim> right_hand_side;
+  std::vector<double>  rhs_values (n_q_points);
+
+  const Solution<dim> exact_solution;
+  
+  typename hp::DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
+    {
+      cell_matrix = 0;
+      cell_rhs = 0;
+
+      x_fe_values.reinit (cell);
+      const FEValues<dim> &fe_values = x_fe_values.get_present_fe_values();
+
+      right_hand_side.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)
+	  {
+	    for (unsigned int j=0; j<dofs_per_cell; ++j)
+	      cell_matrix(i,j) += ((fe_values.shape_grad(i,q_point) *
+				    fe_values.shape_grad(j,q_point)
+                                    +
+                                    fe_values.shape_value(i,q_point) *
+				    fe_values.shape_value(j,q_point)) *
+                                   fe_values.JxW(q_point));
+
+	    cell_rhs(i) += (fe_values.shape_value(i,q_point) *
+			    rhs_values [q_point] *
+			    fe_values.JxW(q_point));
+	  }
+
+      for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+	if (cell->face(face)->at_boundary()
+            &&
+            (cell->face(face)->boundary_indicator() == 1))
+	  {
+	    x_fe_face_values.reinit (cell, face);
+	    const FEFaceValues<dim> &fe_face_values
+	      = x_fe_face_values.get_present_fe_values ();
+
+	    for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
+	      {
+		const double neumann_value
+		  = (exact_solution.gradient (fe_face_values.quadrature_point(q_point)) *
+		     fe_face_values.normal_vector(q_point));
+
+		for (unsigned int i=0; i<dofs_per_cell; ++i)
+		  cell_rhs(i) += (neumann_value *
+				  fe_face_values.shape_value(i,q_point) *
+				  fe_face_values.JxW(q_point));
+	      }
+	  }
+
+      cell->get_dof_indices (local_dof_indices);
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+	{
+	  for (unsigned int j=0; j<dofs_per_cell; ++j)
+	    system_matrix.add (local_dof_indices[i],
+			       local_dof_indices[j],
+			       cell_matrix(i,j));
+	  
+	  system_rhs(local_dof_indices[i]) += cell_rhs(i);
+	}
+    }
+
+  hanging_node_constraints.condense (system_matrix);
+  hanging_node_constraints.condense (system_rhs);
+
+  std::map<unsigned int,double> boundary_values;
+  VectorTools::interpolate_boundary_values (dof_handler,
+					    0,
+					    Solution<dim>(),
+					    boundary_values);
+  MatrixTools::apply_boundary_values (boundary_values,
+				      system_matrix,
+				      solution,
+				      system_rhs);
+}
+
+
+
+template <int dim>
+void HelmholtzProblem<dim>::solve () 
+{
+  SolverControl           solver_control (1000, 1e-12);
+  SolverCG<>              cg (solver_control);
+
+  PreconditionSSOR<> preconditioner;
+  preconditioner.initialize(system_matrix, 1.2);
+
+  cg.solve (system_matrix, solution, system_rhs,
+	    preconditioner);
+
+  hanging_node_constraints.distribute (solution);
+}
+
+
+
+template <int dim>
+void HelmholtzProblem<dim>::refine_grid ()
+{
+  switch (refinement_mode) 
+    {
+      case global_refinement:
+      {
+	triangulation.refine_global (1);
+	break;
+      }
+
+      case adaptive_refinement:
+      {
+	Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+
+	typename FunctionMap<dim>::type neumann_boundary;
+	KellyErrorEstimator<dim>::estimate (dof_handler,
+					    hp::QCollection<dim-1>(3),
+					    neumann_boundary,
+					    solution,
+					    estimated_error_per_cell);
+
+	GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+							 estimated_error_per_cell,
+							 0.3, 0.03);
+	
+	triangulation.execute_coarsening_and_refinement ();
+
+	break;
+      }
+
+      default:
+      {
+        Assert (false, ExcNotImplemented());
+      }
+    }
+}
+
+
+
+template <int dim>
+void HelmholtzProblem<dim>::process_solution (const unsigned int cycle)
+{
+  Vector<float> difference_per_cell (triangulation.n_active_cells());
+  VectorTools::integrate_difference (dof_handler,
+				     solution,
+				     Solution<dim>(),
+				     difference_per_cell,
+				     hp::QCollection<dim>(3),
+				     VectorTools::L2_norm);
+  const double L2_error = difference_per_cell.l2_norm();
+
+  VectorTools::integrate_difference (dof_handler,
+				     solution,
+				     Solution<dim>(),
+				     difference_per_cell,
+				     hp::QCollection<dim>(3),
+				     VectorTools::H1_seminorm);
+  const double H1_error = difference_per_cell.l2_norm();
+
+  const QTrapez<1>     q_trapez;
+  const QIterated<dim> q_iterated (q_trapez, 5);
+  VectorTools::integrate_difference (dof_handler,
+				     solution,
+				     Solution<dim>(),
+				     difference_per_cell,
+				     q_iterated,
+				     VectorTools::Linfty_norm);
+  const double Linfty_error = difference_per_cell.linfty_norm();
+
+  const unsigned int n_active_cells=triangulation.n_active_cells();
+  const unsigned int n_dofs=dof_handler.n_dofs();
+  
+  deallog << "Cycle " << cycle << ':' 
+	    << std::endl
+	    << "   Number of active cells:       "
+	    << n_active_cells
+	    << std::endl
+	    << "   Number of degrees of freedom: "
+	    << n_dofs
+	    << std::endl;
+
+  convergence_table.add_value("cycle", cycle);
+  convergence_table.add_value("cells", n_active_cells);
+  convergence_table.add_value("dofs", n_dofs);
+  convergence_table.add_value("L2", L2_error);
+  convergence_table.add_value("H1", H1_error);
+  convergence_table.add_value("Linfty", Linfty_error);
+}
+
+
+
+template <int dim>
+void HelmholtzProblem<dim>::run () 
+{
+  for (unsigned int cycle=0; cycle<7; ++cycle)
+    {
+      if (cycle == 0)
+	{
+	  GridGenerator::hyper_cube (triangulation, -1, 1);
+	  triangulation.refine_global (1);
+
+          typename Triangulation<dim>::cell_iterator
+            cell = triangulation.begin (),
+            endc = triangulation.end();
+	  for (; cell!=endc; ++cell)
+	    for (unsigned int face=0;
+                 face<GeometryInfo<dim>::faces_per_cell;
+                 ++face)
+	      if ((cell->face(face)->center()(0) == -1)
+		  ||
+		  (cell->face(face)->center()(1) == -1))
+		cell->face(face)->set_boundary_indicator (1);
+	}
+      else
+        refine_grid ();
+      
+
+      setup_system ();
+      
+      assemble_system ();
+      solve ();
+
+      process_solution (cycle);
+    }
+  
+  
+  std::string gmv_filename;
+  switch (refinement_mode)
+    {
+      case global_refinement:
+	    gmv_filename = "solution-global";
+	    break;
+      case adaptive_refinement:
+	    gmv_filename = "solution-adaptive";
+	    break;
+      default:
+	    Assert (false, ExcNotImplemented());
+    }
+
+  switch ((*fe)[0].degree)
+    {
+      case 1:
+	    gmv_filename += "-q1";
+	    break;
+      case 2:
+	    gmv_filename += "-q2";
+	    break;
+
+      default:
+	    Assert (false, ExcNotImplemented());
+    }
+
+  gmv_filename += ".gmv";    
+  std::ofstream output (gmv_filename.c_str());
+
+  DataOut<dim,hp::DoFHandler<dim> > data_out;
+  data_out.attach_dof_handler (dof_handler);
+  data_out.add_data_vector (solution, "solution");
+
+  data_out.build_patches ((*fe)[0].degree);
+  data_out.write_gmv (output);
+
+  
+  
+  convergence_table.set_precision("L2", 3);
+  convergence_table.set_precision("H1", 3);
+  convergence_table.set_precision("Linfty", 3);
+
+  convergence_table.set_scientific("L2", true);
+  convergence_table.set_scientific("H1", true);
+  convergence_table.set_scientific("Linfty", true);
+
+  convergence_table.set_tex_caption("cells", "\\# cells");
+  convergence_table.set_tex_caption("dofs", "\\# dofs");
+  convergence_table.set_tex_caption("L2", "$L^2$-error");
+  convergence_table.set_tex_caption("H1", "$H^1$-error");
+  convergence_table.set_tex_caption("Linfty", "$L^\\infty$-error");
+
+  convergence_table.set_tex_format("cells", "r");
+  convergence_table.set_tex_format("dofs", "r");
+
+  deallog << std::endl;
+  convergence_table.write_text(deallog);
+  
+  std::string error_filename = "error";
+  switch (refinement_mode)
+    {
+      case global_refinement:
+            error_filename += "-global";
+            break;
+      case adaptive_refinement:
+            error_filename += "-adaptive";
+            break;
+      default:
+            Assert (false, ExcNotImplemented());
+    }
+
+  switch ((*fe)[0].degree)
+    {
+      case 1:
+            error_filename += "-q1";
+            break;
+      case 2:
+            error_filename += "-q2";
+            break;
+      default:
+            Assert (false, ExcNotImplemented());
+    }
+  
+  error_filename += ".tex";
+  std::ofstream error_table_file(error_filename.c_str());
+  
+  convergence_table.write_tex(error_table_file);
+
+  
+
+  if (refinement_mode==global_refinement)
+    {
+      convergence_table.add_column_to_supercolumn("cycle", "n cells");
+      convergence_table.add_column_to_supercolumn("cells", "n cells");
+      
+      std::vector<std::string> new_order;
+      new_order.push_back("n cells");
+      new_order.push_back("H1");
+      new_order.push_back("L2");
+      convergence_table.set_column_order (new_order);
+
+      convergence_table
+	.evaluate_convergence_rates("L2", ConvergenceTable::reduction_rate);
+      convergence_table
+	.evaluate_convergence_rates("L2", ConvergenceTable::reduction_rate_log2);
+      convergence_table
+	.evaluate_convergence_rates("H1", ConvergenceTable::reduction_rate_log2);
+
+      deallog << std::endl;
+      convergence_table.write_text(deallog);
+
+      std::string conv_filename = "convergence";
+      switch (refinement_mode)
+	{
+	  case global_refinement:
+		conv_filename += "-global";
+		break;
+	  case adaptive_refinement:
+		conv_filename += "-adaptive";
+		break;
+	  default:
+		Assert (false, ExcNotImplemented());
+	}
+      switch ((*fe)[0].degree)
+	{
+	  case 1:
+		conv_filename += "-q1";
+		break;
+	  case 2:
+		conv_filename += "-q2";
+		break;
+	  default:
+		Assert (false, ExcNotImplemented());
+	}
+      conv_filename += ".tex";
+
+      std::ofstream table_file(conv_filename.c_str());
+      convergence_table.write_tex(table_file);
+    }
+}
+
+
+int main () 
+{
+  logfile.precision(2);
+  
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);  
+
+  const unsigned int dim = 2;
+
+  try
+    {
+      deallog.depth_console (0);
+      
+      
+      {
+	deallog << "Solving with Q1 elements, adaptive refinement" << std::endl
+		  << "=============================================" << std::endl
+		  << std::endl;
+	
+	hp::FECollection<dim> fe(FE_Q<dim>(1));
+	HelmholtzProblem<dim>
+	  helmholtz_problem_2d (fe, HelmholtzProblem<dim>::adaptive_refinement);
+
+	helmholtz_problem_2d.run ();
+
+	deallog << std::endl;
+      }
+	
+      {
+	deallog << "Solving with Q1 elements, global refinement" << std::endl
+		  << "===========================================" << std::endl
+		  << std::endl;
+	
+	hp::FECollection<dim> fe(FE_Q<dim>(1));
+	HelmholtzProblem<dim>
+	  helmholtz_problem_2d (fe, HelmholtzProblem<dim>::global_refinement);
+
+	helmholtz_problem_2d.run ();
+
+	deallog << std::endl;
+      }
+       
+      {
+	deallog << "Solving with Q2 elements, global refinement" << std::endl
+		  << "===========================================" << std::endl
+		  << std::endl;
+	
+	hp::FECollection<dim> fe(FE_Q<dim>(2));
+	HelmholtzProblem<dim>
+	  helmholtz_problem_2d (fe, HelmholtzProblem<dim>::global_refinement);
+
+	helmholtz_problem_2d.run ();
+
+	deallog << std::endl;
+      }
+       
+    }
+  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;
+}
+
+
+template const double SolutionBase<2>::width;
diff --git a/tests/hp/step-8.cc b/tests/hp/step-8.cc
new file mode 100644
index 0000000000..4744612bbb
--- /dev/null
+++ b/tests/hp/step-8.cc
@@ -0,0 +1,468 @@
+//----------------------------  step-8.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2005, 2006 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  step-8.cc  ---------------------------
+
+
+// a hp-ified version of step-8
+
+
+#include <base/logstream.h>
+#include <fstream>
+std::ofstream logfile("step-8/output");
+
+
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+#include <base/logstream.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <lac/sparse_matrix.h>
+#include <lac/solver_cg.h>
+#include <lac/precondition.h>
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/grid_refinement.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_boundary_lib.h>
+#include <dofs/hp_dof_handler.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <fe/hp_fe_values.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+#include <dofs/dof_constraints.h>
+#include <numerics/error_estimator.h>
+
+#include <fe/fe_system.h>
+#include <fe/fe_q.h>
+
+#include <fstream>
+#include <iostream>
+
+
+
+template <int dim>
+class ElasticProblem 
+{
+  public:
+    ElasticProblem ();
+    ~ElasticProblem ();
+    void run ();
+    
+  private:
+    void setup_system ();
+    void assemble_system ();
+    void solve ();
+    void refine_grid ();
+    void output_results (const unsigned int cycle) const;
+
+    Triangulation<dim>   triangulation;
+    hp::DoFHandler<dim>      dof_handler;
+
+    hp::FECollection<dim>        fe;
+
+    ConstraintMatrix     hanging_node_constraints;
+
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+
+    Vector<double>       solution;
+    Vector<double>       system_rhs;
+};
+
+
+
+template <int dim>
+class RightHandSide :  public Function<dim> 
+{
+  public:
+    RightHandSide ();
+    
+    virtual void vector_value (const Point<dim> &p,
+			       Vector<double>   &values) const;
+
+    virtual void vector_value_list (const std::vector<Point<dim> > &points,
+				    std::vector<Vector<double> >   &value_list) const;
+};
+
+
+template <int dim>
+RightHandSide<dim>::RightHandSide ()
+		:
+		Function<dim> (dim)
+{}
+
+
+template <int dim>
+inline
+void RightHandSide<dim>::vector_value (const Point<dim> &p,
+				       Vector<double>   &values) const 
+{
+  Assert (values.size() == dim, 
+	  ExcDimensionMismatch (values.size(), dim));
+  Assert (dim >= 2, ExcNotImplemented());
+  
+  Point<dim> point_1, point_2;
+  point_1(0) = 0.5;
+  point_2(0) = -0.5;
+  
+  if (((p-point_1).square() < 0.2*0.2) ||
+      ((p-point_2).square() < 0.2*0.2))
+    values(0) = 1;
+  else
+    values(0) = 0;
+  
+  if (p.square() < 0.2*0.2)
+    values(1) = 1;
+  else
+    values(1) = 0;    
+}
+
+
+
+template <int dim>
+void RightHandSide<dim>::vector_value_list (const std::vector<Point<dim> > &points,
+					    std::vector<Vector<double> >   &value_list) const 
+{
+  Assert (value_list.size() == points.size(), 
+	  ExcDimensionMismatch (value_list.size(), points.size()));
+
+  const unsigned int n_points = points.size();
+
+  for (unsigned int p=0; p<n_points; ++p)
+    RightHandSide<dim>::vector_value (points[p],
+				      value_list[p]);
+}
+
+
+
+
+
+
+template <int dim>
+ElasticProblem<dim>::ElasticProblem ()
+		:
+		dof_handler (triangulation),
+		fe (FESystem<dim>(FE_Q<dim>(1), dim))
+{}
+
+
+
+
+template <int dim>
+ElasticProblem<dim>::~ElasticProblem () 
+{
+  dof_handler.clear ();
+}
+
+
+
+template <int dim>
+void ElasticProblem<dim>::setup_system ()
+{
+  dof_handler.distribute_dofs (fe);
+  hanging_node_constraints.clear ();
+  DoFTools::make_hanging_node_constraints (dof_handler,
+					   hanging_node_constraints);
+  hanging_node_constraints.close ();
+  sparsity_pattern.reinit (dof_handler.n_dofs(),
+			   dof_handler.n_dofs(),
+			   dof_handler.max_couplings_between_dofs());
+  DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
+
+  hanging_node_constraints.condense (sparsity_pattern);
+
+  sparsity_pattern.compress();
+
+  system_matrix.reinit (sparsity_pattern);
+
+  solution.reinit (dof_handler.n_dofs());
+  system_rhs.reinit (dof_handler.n_dofs());
+}
+
+
+
+template <int dim>
+void ElasticProblem<dim>::assemble_system () 
+{  
+  hp::QCollection<dim>  quadrature_formula(2);
+
+  hp::FEValues<dim> x_fe_values (fe, quadrature_formula, 
+			   update_values   | update_gradients |
+                           update_q_points | update_JxW_values);
+
+  const unsigned int   dofs_per_cell = fe[0].dofs_per_cell;
+  const unsigned int   n_q_points    = quadrature_formula[0].n_quadrature_points;
+
+  FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+  Vector<double>       cell_rhs (dofs_per_cell);
+
+  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+
+  std::vector<double>     lambda_values (n_q_points);
+  std::vector<double>     mu_values (n_q_points);
+
+  ConstantFunction<dim> lambda(1.), mu(1.);
+
+  RightHandSide<dim>      right_hand_side;
+  std::vector<Vector<double> > rhs_values (n_q_points,
+					   Vector<double>(dim));
+
+
+  typename hp::DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
+						 endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
+    {
+      cell_matrix = 0;
+      cell_rhs = 0;
+
+      x_fe_values.reinit (cell);
+      const FEValues<dim> &fe_values = x_fe_values.get_present_fe_values();
+      
+      lambda.value_list (fe_values.get_quadrature_points(), lambda_values);
+      mu.value_list     (fe_values.get_quadrature_points(), mu_values);
+
+      right_hand_side.vector_value_list (fe_values.get_quadrature_points(),
+					 rhs_values);
+      
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+	{
+	  const unsigned int 
+	    component_i = fe[0].system_to_component_index(i).first;
+	  
+	  for (unsigned int j=0; j<dofs_per_cell; ++j) 
+	    {
+	      const unsigned int 
+		component_j = fe[0].system_to_component_index(j).first;
+	      
+	      for (unsigned int q_point=0; q_point<n_q_points;
+		   ++q_point)
+		{
+		  cell_matrix(i,j) 
+		    += 
+		    (
+		      (fe_values.shape_grad(i,q_point)[component_i] *
+		       fe_values.shape_grad(j,q_point)[component_j] *
+		       lambda_values[q_point])
+		      +
+		      (fe_values.shape_grad(i,q_point)[component_j] *
+		       fe_values.shape_grad(j,q_point)[component_i] *
+		       mu_values[q_point])
+		      +
+		      ((component_i == component_j) ?
+		       (fe_values.shape_grad(i,q_point) *
+			fe_values.shape_grad(j,q_point) *
+			mu_values[q_point])  :
+		       0)
+		    )
+		    *
+		    fe_values.JxW(q_point);
+		}
+	    }
+	}
+
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+	{
+	  const unsigned int 
+	    component_i = fe[0].system_to_component_index(i).first;
+	  
+	  for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+	    cell_rhs(i) += fe_values.shape_value(i,q_point) *
+			   rhs_values[q_point](component_i) *
+			   fe_values.JxW(q_point);
+	}
+
+      cell->get_dof_indices (local_dof_indices);
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+	{
+	  for (unsigned int j=0; j<dofs_per_cell; ++j)
+	    system_matrix.add (local_dof_indices[i],
+			       local_dof_indices[j],
+			       cell_matrix(i,j));
+	  
+	  system_rhs(local_dof_indices[i]) += cell_rhs(i);
+	}
+    }
+
+  hanging_node_constraints.condense (system_matrix);
+  hanging_node_constraints.condense (system_rhs);
+
+  std::map<unsigned int,double> boundary_values;
+  VectorTools::interpolate_boundary_values (dof_handler,
+					    0,
+					    ZeroFunction<dim>(dim),
+					    boundary_values);
+  MatrixTools::apply_boundary_values (boundary_values,
+				      system_matrix,
+				      solution,
+				      system_rhs);
+}
+
+
+
+
+template <int dim>
+void ElasticProblem<dim>::solve () 
+{
+  SolverControl           solver_control (1000, 1e-12);
+  SolverCG<>              cg (solver_control);
+
+  PreconditionSSOR<> preconditioner;
+  preconditioner.initialize(system_matrix, 1.2);
+
+  cg.solve (system_matrix, solution, system_rhs,
+	    preconditioner);
+
+  hanging_node_constraints.distribute (solution);
+}
+
+
+
+template <int dim>
+void ElasticProblem<dim>::refine_grid ()
+{
+  Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+
+  typename FunctionMap<dim>::type neumann_boundary;
+  KellyErrorEstimator<dim>::estimate (dof_handler,
+				      hp::QCollection<dim-1>(2),
+				      neumann_boundary,
+				      solution,
+				      estimated_error_per_cell);
+
+  GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+						   estimated_error_per_cell,
+						   0.3, 0.03);
+
+  triangulation.execute_coarsening_and_refinement ();
+}
+
+
+
+template <int dim>
+void ElasticProblem<dim>::output_results (const unsigned int cycle) const
+{
+  std::string filename = "solution-";
+  filename += ('0' + cycle);
+  Assert (cycle < 10, ExcInternalError());
+  
+  filename += ".gmv";
+  std::ofstream output (filename.c_str());
+
+  DataOut<dim,hp::DoFHandler<dim> > data_out;
+  data_out.attach_dof_handler (dof_handler);
+
+ 
+
+  std::vector<std::string> solution_names;
+  switch (dim)
+    {
+      case 1:
+	    solution_names.push_back ("displacement");
+	    break;
+      case 2:
+	    solution_names.push_back ("x_displacement");	    
+	    solution_names.push_back ("y_displacement");
+	    break;
+      case 3:
+	    solution_names.push_back ("x_displacement");	    
+	    solution_names.push_back ("y_displacement");
+	    solution_names.push_back ("z_displacement");
+	    break;
+      default:
+	    Assert (false, ExcNotImplemented());
+    }
+	     
+  data_out.add_data_vector (solution, solution_names);
+  data_out.build_patches ();
+  data_out.write_gmv (output);
+}
+
+
+
+
+template <int dim>
+void ElasticProblem<dim>::run () 
+{
+  for (unsigned int cycle=0; cycle<8; ++cycle)
+    {
+      deallog << "Cycle " << cycle << ':' << std::endl;
+
+      if (cycle == 0)
+	{
+	  GridGenerator::hyper_cube (triangulation, -1, 1);
+	  triangulation.refine_global (2);
+	}
+      else
+	refine_grid ();
+
+      deallog << "   Number of active cells:       "
+		<< triangulation.n_active_cells()
+		<< std::endl;
+
+      setup_system ();
+
+      deallog << "   Number of degrees of freedom: "
+		<< dof_handler.n_dofs()
+		<< std::endl;
+      
+      assemble_system ();
+      solve ();
+      output_results (cycle);
+    }
+}
+
+
+int main () 
+{
+  logfile.precision(2);
+  
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);  
+
+  try
+    {
+      deallog.depth_console (0);
+
+      ElasticProblem<2> elastic_problem_2d;
+      elastic_problem_2d.run ();
+    }
+  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;
+}
diff --git a/tests/hp/step-9.cc b/tests/hp/step-9.cc
new file mode 100644
index 0000000000..c5db4b5e32
--- /dev/null
+++ b/tests/hp/step-9.cc
@@ -0,0 +1,761 @@
+//----------------------------  step-9.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2005, 2006 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  step-9.cc  ---------------------------
+
+
+// a hp-ified version of step-9
+
+
+#include <base/logstream.h>
+#include <fstream>
+std::ofstream logfile("step-9/output");
+
+
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+#include <base/logstream.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <lac/sparse_matrix.h>
+#include <lac/solver_bicgstab.h>
+#include <lac/precondition.h>
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/grid_refinement.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_boundary_lib.h>
+#include <dofs/hp_dof_handler.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <dofs/dof_constraints.h>
+#include <fe/hp_fe_values.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+#include <fe/fe_q.h>
+#include <grid/grid_out.h>
+
+#include <base/thread_management.h>
+#include <base/multithread_info.h>
+
+#include <base/tensor_function.h>
+
+#include <numerics/error_estimator.h>
+
+#include <fstream>
+#include <iostream>
+
+
+
+#ifndef M_PI
+#  define	M_PI		3.14159265358979323846
+#endif
+
+
+template <int dim>
+class AdvectionProblem 
+{
+  public:
+    AdvectionProblem ();
+    ~AdvectionProblem ();
+    void run ();
+    
+  private:
+    void setup_system ();
+    void assemble_system ();
+    void assemble_system_interval (const typename hp::DoFHandler<dim>::active_cell_iterator &begin,
+			 	   const typename hp::DoFHandler<dim>::active_cell_iterator &end);
+    
+    void solve ();
+    void refine_grid ();
+    void output_results (const unsigned int cycle) const;
+
+    Triangulation<dim>   triangulation;
+    hp::DoFHandler<dim>      dof_handler;
+
+    hp::FECollection<dim>            fe;
+
+    ConstraintMatrix     hanging_node_constraints;
+
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+
+    Vector<double>       solution;
+    Vector<double>       system_rhs;
+
+    Threads::ThreadMutex     assembler_lock;
+};
+
+
+
+
+template <int dim>
+class AdvectionField : public TensorFunction<1,dim>
+{
+  public:
+    AdvectionField () : TensorFunction<1,dim> () {};
+    
+    virtual Tensor<1,dim> value (const Point<dim> &p) const;
+    
+    virtual void value_list (const std::vector<Point<dim> > &points,
+			     std::vector<Tensor<1,dim> >    &values) const;
+
+    DeclException2 (ExcDimensionMismatch,
+		    unsigned int, unsigned int,
+		    << "The vector has size " << arg1 << " but should have "
+		    << arg2 << " elements.");
+};
+
+
+
+template <int dim>
+Tensor<1,dim> 
+AdvectionField<dim>::value (const Point<dim> &p) const 
+{
+  Point<dim> value;
+  value[0] = 2;
+  for (unsigned int i=1; i<dim; ++i)
+    value[i] = 1+0.8*std::sin(8*M_PI*p[0]);
+
+  return value;
+}
+
+
+
+template <int dim>
+void
+AdvectionField<dim>::value_list (const std::vector<Point<dim> > &points,
+				 std::vector<Tensor<1,dim> >    &values) const 
+{
+  Assert (values.size() == points.size(),
+	  ExcDimensionMismatch (values.size(), points.size()));
+  
+  for (unsigned int i=0; i<points.size(); ++i)
+    values[i] = AdvectionField<dim>::value (points[i]);
+}
+
+
+
+
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+  public:
+    RightHandSide () : Function<dim>() {};
+    
+    virtual double value (const Point<dim>   &p,
+			  const unsigned int  component = 0) const;
+    
+    virtual void value_list (const std::vector<Point<dim> > &points,
+			     std::vector<double>            &values,
+			     const unsigned int              component = 0) const;
+    
+  private:
+    static const Point<dim> center_point;
+};
+
+
+template <>
+const Point<1> RightHandSide<1>::center_point = Point<1> (-0.75);
+
+template <>
+const Point<2> RightHandSide<2>::center_point = Point<2> (-0.75, -0.75);
+
+template <>
+const Point<3> RightHandSide<3>::center_point = Point<3> (-0.75, -0.75, -0.75);
+
+
+
+template <int dim>
+double
+RightHandSide<dim>::value (const Point<dim>   &p,
+			   const unsigned int  component) const 
+{
+  Assert (component == 0, ExcIndexRange (component, 0, 1));
+  const double diameter = 0.1;
+  return ( (p-center_point).square() < diameter*diameter ?
+	   .1/std::pow(diameter,dim) :
+	   0);
+}
+
+
+
+template <int dim>
+void
+RightHandSide<dim>::value_list (const std::vector<Point<dim> > &points,
+				std::vector<double>            &values,
+				const unsigned int              component) const 
+{
+  Assert (values.size() == points.size(),
+	  ExcDimensionMismatch (values.size(), points.size()));
+  
+  for (unsigned int i=0; i<points.size(); ++i)
+    values[i] = RightHandSide<dim>::value (points[i], component);
+}
+
+
+
+template <int dim>
+class BoundaryValues : public Function<dim>
+{
+  public:
+    BoundaryValues () : Function<dim>() {};
+
+    virtual double value (const Point<dim>   &p,
+			  const unsigned int  component = 0) const;
+    
+    virtual void value_list (const std::vector<Point<dim> > &points,
+			     std::vector<double>            &values,
+			     const unsigned int              component = 0) const;
+};
+
+
+
+template <int dim>
+double
+BoundaryValues<dim>::value (const Point<dim>   &p,
+			    const unsigned int  component) const 
+{
+  Assert (component == 0, ExcIndexRange (component, 0, 1));
+
+  const double sine_term = std::sin(16*M_PI*std::sqrt(p.square()));
+  const double weight    = std::exp(-5*p.square()) / std::exp(-5.);
+  return sine_term * weight;
+}
+
+
+
+template <int dim>
+void
+BoundaryValues<dim>::value_list (const std::vector<Point<dim> > &points,
+				 std::vector<double>            &values,
+				 const unsigned int              component) const 
+{
+  Assert (values.size() == points.size(),
+	  ExcDimensionMismatch (values.size(), points.size()));
+  
+  for (unsigned int i=0; i<points.size(); ++i)
+    values[i] = BoundaryValues<dim>::value (points[i], component);
+}
+
+
+
+
+class GradientEstimation
+{
+  public:
+    template <int dim>
+    static void estimate (const hp::DoFHandler<dim> &dof,
+			  const Vector<double>  &solution,
+			  Vector<float>         &error_per_cell);
+
+    DeclException2 (ExcInvalidVectorLength,
+		    int, int,
+		    << "Vector has length " << arg1 << ", but should have "
+		    << arg2);
+    DeclException0 (ExcInsufficientDirections);
+
+  private:
+    typedef std::pair<unsigned int,unsigned int> IndexInterval;
+
+    template <int dim>
+    static void estimate_interval (const hp::DoFHandler<dim> &dof,
+				   const Vector<double>  &solution,
+				   const IndexInterval   &index_interval,
+				   Vector<float>         &error_per_cell);    
+};
+
+
+
+
+
+template <int dim>
+AdvectionProblem<dim>::AdvectionProblem () :
+		dof_handler (triangulation),
+		fe(FE_Q<dim>(1))
+{}
+
+
+
+template <int dim>
+AdvectionProblem<dim>::~AdvectionProblem () 
+{
+  dof_handler.clear ();
+}
+
+
+
+template <int dim>
+void AdvectionProblem<dim>::setup_system ()
+{
+  dof_handler.distribute_dofs (fe);
+
+  hanging_node_constraints.clear ();
+  DoFTools::make_hanging_node_constraints (dof_handler,
+					   hanging_node_constraints);
+  hanging_node_constraints.close ();
+
+  sparsity_pattern.reinit (dof_handler.n_dofs(),
+			   dof_handler.n_dofs(),
+			   dof_handler.max_couplings_between_dofs());
+  DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
+
+  hanging_node_constraints.condense (sparsity_pattern);
+
+  sparsity_pattern.compress();
+
+  system_matrix.reinit (sparsity_pattern);
+
+  solution.reinit (dof_handler.n_dofs());
+  system_rhs.reinit (dof_handler.n_dofs());
+}
+
+
+
+template <int dim>
+void AdvectionProblem<dim>::assemble_system () 
+{
+  const unsigned int n_threads = multithread_info.n_default_threads;
+  Threads::ThreadGroup<> threads;
+
+  typedef typename hp::DoFHandler<dim>::active_cell_iterator active_cell_iterator;
+  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);
+
+  for (unsigned int thread=0; thread<n_threads; ++thread)
+    threads += Threads::spawn (*this, &AdvectionProblem<dim>::assemble_system_interval)
+               (thread_ranges[thread].first,
+                thread_ranges[thread].second);
+
+  threads.join_all ();  
+
+
+  hanging_node_constraints.condense (system_matrix);
+  hanging_node_constraints.condense (system_rhs);
+}
+
+
+ 
+template <int dim>
+void
+AdvectionProblem<dim>::
+assemble_system_interval (const typename hp::DoFHandler<dim>::active_cell_iterator &begin,
+			  const typename hp::DoFHandler<dim>::active_cell_iterator &end) 
+{
+  const AdvectionField<dim> advection_field;
+  const RightHandSide<dim>  right_hand_side;
+  const BoundaryValues<dim> boundary_values;
+  
+  hp::QCollection<dim>   quadrature_formula(QGauss<dim>(2));
+  hp::QCollection<dim-1> face_quadrature_formula(QGauss<dim-1>(2));
+  
+  hp::FEValues<dim> x_fe_values (fe, quadrature_formula, 
+			   update_values   | update_gradients |
+                           update_q_points | update_JxW_values);
+  hp::FEFaceValues<dim> x_fe_face_values (fe, face_quadrature_formula,
+				    update_values     | update_q_points   |
+                                    update_JxW_values | update_normal_vectors);
+
+  const unsigned int   dofs_per_cell   = fe[0].dofs_per_cell;
+  const unsigned int   n_q_points      = quadrature_formula[0].n_quadrature_points;
+  const unsigned int   n_face_q_points = face_quadrature_formula[0].n_quadrature_points;
+
+  FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+  Vector<double>       cell_rhs (dofs_per_cell);
+
+  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+
+  std::vector<double>         rhs_values (n_q_points);
+  std::vector<Tensor<1,dim> > advection_directions (n_q_points);
+  std::vector<double>         face_boundary_values (n_face_q_points);
+  std::vector<Tensor<1,dim> > face_advection_directions (n_face_q_points);
+
+  typename hp::DoFHandler<dim>::active_cell_iterator cell;
+  for (cell=begin; cell!=end; ++cell)
+    {
+      cell_matrix = 0;
+      cell_rhs = 0;
+
+      x_fe_values.reinit (cell);
+      const FEValues<dim> &fe_values = x_fe_values.get_present_fe_values();
+
+      advection_field.value_list (fe_values.get_quadrature_points(),
+				  advection_directions);
+      right_hand_side.value_list (fe_values.get_quadrature_points(),
+				  rhs_values);
+
+      const double delta = 0.1 * cell->diameter ();
+
+      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) += ((advection_directions[q_point] *
+				    fe_values.shape_grad(j,q_point)   *
+				    (fe_values.shape_value(i,q_point) +
+				     delta *
+				     (advection_directions[q_point] *
+				      fe_values.shape_grad(i,q_point)))) *
+				   fe_values.JxW(q_point));
+
+	    cell_rhs(i) += ((fe_values.shape_value(i,q_point) +
+			     delta *
+			     (advection_directions[q_point] *
+			      fe_values.shape_grad(i,q_point))        ) *
+			    rhs_values[i] *
+			    fe_values.JxW (q_point));
+	  };
+
+      for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+	if (cell->face(face)->at_boundary())
+	  {
+	    x_fe_face_values.reinit (cell, face);
+	    const FEFaceValues<dim> &fe_face_values
+	      = x_fe_face_values.get_present_fe_values ();
+	    
+	    boundary_values.value_list (fe_face_values.get_quadrature_points(),
+					face_boundary_values);
+	    advection_field.value_list (fe_face_values.get_quadrature_points(),
+					face_advection_directions);
+	    
+	    for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
+	      if (fe_face_values.normal_vector(q_point) *
+		  face_advection_directions[q_point]
+		  < 0)
+		for (unsigned int i=0; i<dofs_per_cell; ++i)
+		  {
+		    for (unsigned int j=0; j<dofs_per_cell; ++j)
+		      cell_matrix(i,j) -= (face_advection_directions[q_point] *
+					   fe_face_values.normal_vector(q_point) *
+					   fe_face_values.shape_value(i,q_point) *
+					   fe_face_values.shape_value(j,q_point) *
+					   fe_face_values.JxW(q_point));
+		    
+		    cell_rhs(i) -= (face_advection_directions[q_point] *
+				    fe_face_values.normal_vector(q_point) *
+				    face_boundary_values[q_point]         *
+				    fe_face_values.shape_value(i,q_point) *
+				    fe_face_values.JxW(q_point));
+		  };
+	  };
+      
+
+      cell->get_dof_indices (local_dof_indices);
+
+      assembler_lock.acquire ();
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+	{
+	  for (unsigned int j=0; j<dofs_per_cell; ++j)
+	    system_matrix.add (local_dof_indices[i],
+			       local_dof_indices[j],
+			       cell_matrix(i,j));
+	  
+	  system_rhs(local_dof_indices[i]) += cell_rhs(i);
+	};
+      assembler_lock.release ();
+    };
+}
+
+
+
+template <int dim>
+void AdvectionProblem<dim>::solve () 
+{
+  SolverControl           solver_control (1000, 1e-12);
+  SolverBicgstab<>        bicgstab (solver_control);
+
+  PreconditionJacobi<> preconditioner;
+  preconditioner.initialize(system_matrix, 1.0);
+
+  bicgstab.solve (system_matrix, solution, system_rhs,
+		  preconditioner);
+
+  hanging_node_constraints.distribute (solution);
+}
+
+
+template <int dim>
+void AdvectionProblem<dim>::refine_grid ()
+{
+  Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+
+  GradientEstimation::estimate (dof_handler,
+				solution,
+				estimated_error_per_cell);
+
+  GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+						   estimated_error_per_cell,
+						   0.5, 0.03);
+
+  triangulation.execute_coarsening_and_refinement ();
+}
+
+
+
+template <int dim>
+void AdvectionProblem<dim>::output_results (const unsigned int cycle) const
+{
+  std::string filename = "grid-";
+  filename += ('0' + cycle);
+  Assert (cycle < 10, ExcInternalError());
+  
+  filename += ".eps";
+  std::ofstream output (filename.c_str());
+
+  GridOut grid_out;
+  grid_out.write_eps (triangulation, output);
+}
+
+
+template <int dim>
+void AdvectionProblem<dim>::run () 
+{
+  for (unsigned int cycle=0; cycle<6; ++cycle)
+    {
+      deallog << "Cycle " << cycle << ':' << std::endl;
+
+      if (cycle == 0)
+	{
+	  GridGenerator::hyper_cube (triangulation, -1, 1);
+	  triangulation.refine_global (4);
+	}
+      else
+	{
+	  refine_grid ();
+	};
+      
+
+      deallog << "   Number of active cells:       "
+		<< triangulation.n_active_cells()
+		<< std::endl;
+
+      setup_system ();
+
+      deallog << "   Number of degrees of freedom: "
+		<< dof_handler.n_dofs()
+		<< std::endl;
+      
+      assemble_system ();
+      solve ();
+      output_results (cycle);
+    };
+
+  DataOut<dim,hp::DoFHandler<dim> > data_out;
+  data_out.attach_dof_handler (dof_handler);
+  data_out.add_data_vector (solution, "solution");
+  data_out.build_patches ();
+  
+  std::ofstream output ("final-solution.gmv");
+  data_out.write_gmv (output);
+}
+
+
+
+
+template <int dim>
+void 
+GradientEstimation::estimate (const hp::DoFHandler<dim> &dof_handler,
+			      const Vector<double>  &solution,
+			      Vector<float>         &error_per_cell)
+{
+  Assert (error_per_cell.size() == dof_handler.get_tria().n_active_cells(),
+	  ExcInvalidVectorLength (error_per_cell.size(),
+				  dof_handler.get_tria().n_active_cells()));
+
+  const unsigned int n_threads = multithread_info.n_default_threads;
+  std::vector<IndexInterval> index_intervals
+    = Threads::split_interval (0, dof_handler.get_tria().n_active_cells(),
+			       n_threads);
+
+  Threads::ThreadGroup<> threads;
+  void (*estimate_interval_ptr) (const hp::DoFHandler<dim> &,
+				 const Vector<double> &,
+				 const IndexInterval &,
+				 Vector<float> &)
+    = &GradientEstimation::template estimate_interval<dim>;
+  for (unsigned int i=0; i<n_threads; ++i)
+    threads += Threads::spawn (estimate_interval_ptr)(dof_handler, solution,
+                                                      index_intervals[i],
+                                                      error_per_cell);
+  threads.join_all ();
+}
+
+
+template <int dim>
+void 
+GradientEstimation::estimate_interval (const hp::DoFHandler<dim> &dof_handler,
+				       const Vector<double>  &solution,
+				       const IndexInterval   &index_interval,
+				       Vector<float>         &error_per_cell)
+{
+  QMidpoint<dim> xmidpoint;
+  hp::QCollection<dim> midpoint_rule (xmidpoint);
+  hp::FEValues<dim>  x_fe_midpoint_value (dof_handler.get_fe(),
+				    midpoint_rule,
+				    update_values | update_q_points);
+  
+  Tensor<2,dim> Y;
+
+  typename hp::DoFHandler<dim>::active_cell_iterator cell, endc;
+
+  cell = dof_handler.begin_active();
+  advance (cell, static_cast<signed int>(index_interval.first));
+
+  endc = dof_handler.begin_active();
+  advance (endc, static_cast<signed int>(index_interval.second));
+
+  Vector<float>::iterator
+    error_on_this_cell = error_per_cell.begin() + index_interval.first;
+  
+
+  std::vector<typename hp::DoFHandler<dim>::active_cell_iterator> active_neighbors;
+  active_neighbors.reserve (GeometryInfo<dim>::faces_per_cell *
+			    GeometryInfo<dim>::subfaces_per_face);
+
+  for (; cell!=endc; ++cell, ++error_on_this_cell)
+    {
+      x_fe_midpoint_value.reinit (cell);
+      const FEValues<dim> &fe_midpoint_value = x_fe_midpoint_value.get_present_fe_values ();
+      
+      Y.clear ();
+
+      Tensor<1,dim> projected_gradient;
+
+
+      active_neighbors.clear ();
+      for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+	if (! cell->at_boundary(face_no))
+	  {
+	    const typename hp::DoFHandler<dim>::face_iterator 
+	      face = cell->face(face_no);
+	    const typename hp::DoFHandler<dim>::cell_iterator 
+	      neighbor = cell->neighbor(face_no);
+
+	    if (neighbor->active())
+	      active_neighbors.push_back (neighbor);
+	    else
+	      {
+		if (dim == 1)
+		  {
+		    typename hp::DoFHandler<dim>::cell_iterator
+		      neighbor_child = neighbor;
+		    while (neighbor_child->has_children())
+		      neighbor_child = neighbor_child->child (face_no==0 ? 1 : 0);
+		    
+		    Assert (neighbor_child->neighbor(face_no==0 ? 1 : 0)==cell,
+			    ExcInternalError());
+		    
+		    active_neighbors.push_back (neighbor_child);
+		  }
+		else
+		  for (unsigned int subface_no=0; subface_no<face->n_children(); ++subface_no)
+		    active_neighbors.push_back (
+		      cell->neighbor_child_on_subface(face_no, subface_no));
+	      };
+	  };
+
+      const Point<dim> this_center = fe_midpoint_value.quadrature_point(0);
+
+      std::vector<double> this_midpoint_value(1);
+      fe_midpoint_value.get_function_values (solution, this_midpoint_value);
+		
+
+      std::vector<double> neighbor_midpoint_value(1);
+      typename std::vector<typename hp::DoFHandler<dim>::active_cell_iterator>::const_iterator
+	neighbor_ptr = active_neighbors.begin();
+      for (; neighbor_ptr!=active_neighbors.end(); ++neighbor_ptr)
+	{
+	  const typename hp::DoFHandler<dim>::active_cell_iterator
+	    neighbor = *neighbor_ptr;
+	    
+	  x_fe_midpoint_value.reinit (neighbor);
+	  const FEValues<dim> &fe_midpoint_value = x_fe_midpoint_value.get_present_fe_values ();
+	  const Point<dim> neighbor_center = fe_midpoint_value.quadrature_point(0);
+
+	  fe_midpoint_value.get_function_values (solution,
+                                                 neighbor_midpoint_value);
+
+	  Point<dim>   y        = neighbor_center - this_center;
+	  const double distance = std::sqrt(y.square());
+	  y /= distance;
+	  
+	  for (unsigned int i=0; i<dim; ++i)
+	    for (unsigned int j=0; j<dim; ++j)
+	      Y[i][j] += y[i] * y[j];
+	  
+	  projected_gradient += (neighbor_midpoint_value[0] -
+				 this_midpoint_value[0]) /
+				distance *
+				y;
+	};
+
+      AssertThrow (determinant(Y) != 0,
+		   ExcInsufficientDirections());
+
+      const Tensor<2,dim> Y_inverse = invert(Y);
+      
+      Point<dim> gradient;
+      contract (gradient, Y_inverse, projected_gradient);
+      
+      *error_on_this_cell = (std::pow(cell->diameter(),
+				      1+1.0*dim/2) *
+			     std::sqrt(gradient.square()));
+    };
+}
+
+
+
+int main () 
+{
+  logfile.precision(2);
+  
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);  
+
+  try
+    {
+      deallog.depth_console (0);
+
+      AdvectionProblem<2> advection_problem_2d;
+      advection_problem_2d.run ();
+    }
+  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;
+}