]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add more tests that all fail for various reasons at the moment.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 3 May 2006 15:04:28 +0000 (15:04 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 3 May 2006 15:04:28 +0000 (15:04 +0000)
git-svn-id: https://svn.dealii.org/trunk@13018 0785d39b-7218-0410-832d-ea1e28bc413d

tests/hp/step-10.cc [new file with mode: 0644]
tests/hp/step-6.cc [new file with mode: 0644]
tests/hp/step-7.cc [new file with mode: 0644]
tests/hp/step-8.cc [new file with mode: 0644]
tests/hp/step-9.cc [new file with mode: 0644]

diff --git a/tests/hp/step-10.cc b/tests/hp/step-10.cc
new file mode 100644 (file)
index 0000000..733b981
--- /dev/null
@@ -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 (file)
index 0000000..2731346
--- /dev/null
@@ -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 (file)
index 0000000..e053f7c
--- /dev/null
@@ -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 (file)
index 0000000..4744612
--- /dev/null
@@ -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 (file)
index 0000000..c5db4b5
--- /dev/null
@@ -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;
+}

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.