]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Partly hp-ify these tests. Not all library calls work right now
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 31 Jul 2006 19:41:45 +0000 (19:41 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 31 Jul 2006 19:41:45 +0000 (19:41 +0000)
git-svn-id: https://svn.dealii.org/trunk@13548 0785d39b-7218-0410-832d-ea1e28bc413d

tests/hp/step-11.cc
tests/hp/step-12.cc
tests/hp/step-13.cc
tests/hp/step-14.cc
tests/hp/step-15.cc

index c2d88a1ea4cddcf76078f6ed411d304fe133834c..c1247d8948f1949375d823b18f3c7864633623e0 100644 (file)
@@ -37,7 +37,7 @@ std::ofstream logfile("step-11/output");
 #include <dofs/dof_accessor.h>
 #include <dofs/dof_tools.h>
 #include <fe/fe_q.h>
-#include <fe/fe_values.h>
+#include <fe/hp_fe_values.h>
 #include <fe/mapping_q.h>
 #include <numerics/vectors.h>
 #include <numerics/matrices.h>
@@ -63,9 +63,9 @@ class LaplaceProblem
     void solve ();
 
     Triangulation<dim>   triangulation;
-    FE_Q<dim>            fe;
-    DoFHandler<dim>      dof_handler;
-    MappingQ<dim>        mapping;
+    hp::FECollection<dim>            fe;
+    hp::DoFHandler<dim>      dof_handler;
+    hp::MappingCollection<dim>        mapping;
 
     SparsityPattern      sparsity_pattern;
     SparseMatrix<double> system_matrix;
@@ -81,9 +81,9 @@ class LaplaceProblem
 
 template <int dim>
 LaplaceProblem<dim>::LaplaceProblem (const unsigned int mapping_degree) :
-                fe (1),
+                fe (FE_Q<dim>(1)),
                dof_handler (triangulation),
-               mapping (mapping_degree)
+               mapping (MappingQ<dim>(mapping_degree))
 {
   deallog << "Using mapping with degree " << mapping_degree << ":"
            << std::endl
@@ -134,7 +134,7 @@ void LaplaceProblem<dim>::assemble_and_solve ()
 {
 
   const unsigned int gauss_degree
-    = std::max (static_cast<unsigned int>(std::ceil(1.*(mapping.get_degree()+1)/2)),
+    = std::max (static_cast<unsigned int>(std::ceil(1.*(static_cast<const MappingQ<dim>&>(mapping[0]).get_degree()+1)/2)),
                2U);
   MatrixTools::create_laplace_matrix (mapping, dof_handler,
                                      QGauss<dim>(gauss_degree),
index bb4b20849be1f8d8a867f05a640af4f0337a33ce..42b1e0fbecb19398e7b615a45ac99638f83e01c7 100644 (file)
@@ -29,7 +29,7 @@ std::ofstream logfile("step-12/output");
 #include <grid/grid_refinement.h>
 #include <grid/tria_accessor.h>
 #include <grid/tria_iterator.h>
-#include <fe/fe_values.h>
+#include <fe/hp_fe_values.h>
 #include <dofs/dof_handler.h>
 #include <dofs/dof_accessor.h>
 #include <dofs/dof_tools.h>
@@ -132,21 +132,23 @@ class DGTransportEquation
   public:
     DGTransportEquation();
 
-    void assemble_cell_term(const FEValues<dim>& fe_v,
+    void assemble_cell_term(const hp::FEValues<dim>& fe_v,
                            FullMatrix<double> &ui_vi_matrix,
                            Vector<double> &cell_vector) const;
     
-    void assemble_boundary_term(const FEFaceValues<dim>& fe_v,
+    void assemble_boundary_term(const hp::FEFaceValues<dim>& fe_v,
                                FullMatrix<double> &ui_vi_matrix,
                                Vector<double> &cell_vector) const;
-    
-    void assemble_face_term1(const FEFaceValuesBase<dim>& fe_v,
-                            const FEFaceValuesBase<dim>& fe_v_neighbor,
+
+    template <class X, class Y>
+    void assemble_face_term1(const X& fe_v,
+                            const Y& fe_v_neighbor,
                             FullMatrix<double> &ui_vi_matrix,
                             FullMatrix<double> &ue_vi_matrix) const;
 
-    void assemble_face_term2(const FEFaceValuesBase<dim>& fe_v,
-                            const FEFaceValuesBase<dim>& fe_v_neighbor,
+    template <class X, class Y>
+    void assemble_face_term2(const X& fe_v,
+                            const Y& fe_v_neighbor,
                             FullMatrix<double> &ui_vi_matrix,
                             FullMatrix<double> &ue_vi_matrix,
                             FullMatrix<double> &ui_ve_matrix,
@@ -169,149 +171,151 @@ DGTransportEquation<dim>::DGTransportEquation ()
 
 template <int dim>
 void DGTransportEquation<dim>::assemble_cell_term(
-  const FEValues<dim> &fe_v,
+  const hp::FEValues<dim> &fe_v,
   FullMatrix<double> &ui_vi_matrix,
   Vector<double> &cell_vector) const
 {
-  const std::vector<double> &JxW = fe_v.get_JxW_values ();
+  const std::vector<double> &JxW = fe_v.get_present_fe_values().get_JxW_values ();
 
-  std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
-  std::vector<double> rhs (fe_v.n_quadrature_points);
+  std::vector<Point<dim> > beta (fe_v.get_present_fe_values().n_quadrature_points);
+  std::vector<double> rhs (fe_v.get_present_fe_values().n_quadrature_points);
   
-  beta_function.value_list (fe_v.get_quadrature_points(), beta);
-  rhs_function.value_list (fe_v.get_quadrature_points(), rhs);
+  beta_function.value_list (fe_v.get_present_fe_values().get_quadrature_points(), beta);
+  rhs_function.value_list (fe_v.get_present_fe_values().get_quadrature_points(), rhs);
   
-  for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
-    for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i) 
+  for (unsigned int point=0; point<fe_v.get_present_fe_values().n_quadrature_points; ++point)
+    for (unsigned int i=0; i<fe_v.get_present_fe_values().dofs_per_cell; ++i) 
       {
-       for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
-         ui_vi_matrix(i,j) -= beta[point]*fe_v.shape_grad(i,point)*
-                             fe_v.shape_value(j,point) *
+       for (unsigned int j=0; j<fe_v.get_present_fe_values().dofs_per_cell; ++j)
+         ui_vi_matrix(i,j) -= beta[point]*fe_v.get_present_fe_values().shape_grad(i,point)*
+                             fe_v.get_present_fe_values().shape_value(j,point) *
                              JxW[point];
        
-       cell_vector(i) += rhs[point] * fe_v.shape_value(i,point) * JxW[point];
+       cell_vector(i) += rhs[point] * fe_v.get_present_fe_values().shape_value(i,point) * JxW[point];
       }
 }
 
 
 template <int dim>
 void DGTransportEquation<dim>::assemble_boundary_term(
-  const FEFaceValues<dim>& fe_v,    
+  const hp::FEFaceValues<dim>& fe_v,    
   FullMatrix<double> &ui_vi_matrix,
   Vector<double> &cell_vector) const
 {
-  const std::vector<double> &JxW = fe_v.get_JxW_values ();
-  const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
+  const std::vector<double> &JxW = fe_v.get_present_fe_values().get_JxW_values ();
+  const std::vector<Point<dim> > &normals = fe_v.get_present_fe_values().get_normal_vectors ();
 
-  std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
-  std::vector<double> g(fe_v.n_quadrature_points);
+  std::vector<Point<dim> > beta (fe_v.get_present_fe_values().n_quadrature_points);
+  std::vector<double> g(fe_v.get_present_fe_values().n_quadrature_points);
   
-  beta_function.value_list (fe_v.get_quadrature_points(), beta);
-  boundary_function.value_list (fe_v.get_quadrature_points(), g);
+  beta_function.value_list (fe_v.get_present_fe_values().get_quadrature_points(), beta);
+  boundary_function.value_list (fe_v.get_present_fe_values().get_quadrature_points(), g);
 
-  for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
+  for (unsigned int point=0; point<fe_v.get_present_fe_values().n_quadrature_points; ++point)
     {
       const double beta_n=beta[point] * normals[point];      
       if (beta_n>0)
-       for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
-         for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
+       for (unsigned int i=0; i<fe_v.get_present_fe_values().dofs_per_cell; ++i)
+         for (unsigned int j=0; j<fe_v.get_present_fe_values().dofs_per_cell; ++j)
            ui_vi_matrix(i,j) += beta_n *
-                              fe_v.shape_value(j,point) *
-                              fe_v.shape_value(i,point) *
+                              fe_v.get_present_fe_values().shape_value(j,point) *
+                              fe_v.get_present_fe_values().shape_value(i,point) *
                               JxW[point];
       else
-       for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
+       for (unsigned int i=0; i<fe_v.get_present_fe_values().dofs_per_cell; ++i)
          cell_vector(i) -= beta_n *
                            g[point] *
-                           fe_v.shape_value(i,point) *
+                           fe_v.get_present_fe_values().shape_value(i,point) *
                            JxW[point];
     }
 }
 
 
 template <int dim>
+template <class X, class Y>
 void DGTransportEquation<dim>::assemble_face_term1(
-  const FEFaceValuesBase<dim>& fe_v,
-  const FEFaceValuesBase<dim>& fe_v_neighbor,      
+  const X& fe_v,
+  const Y& fe_v_neighbor,      
   FullMatrix<double> &ui_vi_matrix,
   FullMatrix<double> &ue_vi_matrix) const
 {
-  const std::vector<double> &JxW = fe_v.get_JxW_values ();
-  const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
+  const std::vector<double> &JxW = fe_v.get_present_fe_values().get_JxW_values ();
+  const std::vector<Point<dim> > &normals = fe_v.get_present_fe_values().get_normal_vectors ();
 
-  std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
-  beta_function.value_list (fe_v.get_quadrature_points(), beta);
+  std::vector<Point<dim> > beta (fe_v.get_present_fe_values().n_quadrature_points);
+  beta_function.value_list (fe_v.get_present_fe_values().get_quadrature_points(), beta);
 
-  for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
+  for (unsigned int point=0; point<fe_v.get_present_fe_values().n_quadrature_points; ++point)
     {
       const double beta_n=beta[point] * normals[point];
       if (beta_n>0)
-       for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
-         for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
+       for (unsigned int i=0; i<fe_v.get_present_fe_values().dofs_per_cell; ++i)
+         for (unsigned int j=0; j<fe_v.get_present_fe_values().dofs_per_cell; ++j)
            ui_vi_matrix(i,j) += beta_n *
-                              fe_v.shape_value(j,point) *
-                              fe_v.shape_value(i,point) *
+                              fe_v.get_present_fe_values().shape_value(j,point) *
+                              fe_v.get_present_fe_values().shape_value(i,point) *
                               JxW[point];
       else
-       for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
-         for (unsigned int k=0; k<fe_v_neighbor.dofs_per_cell; ++k)
+       for (unsigned int i=0; i<fe_v.get_present_fe_values().dofs_per_cell; ++i)
+         for (unsigned int k=0; k<fe_v_neighbor.get_present_fe_values().dofs_per_cell; ++k)
            ue_vi_matrix(i,k) += beta_n *
-                               fe_v_neighbor.shape_value(k,point) *
-                               fe_v.shape_value(i,point) *
+                               fe_v_neighbor.get_present_fe_values().shape_value(k,point) *
+                               fe_v.get_present_fe_values().shape_value(i,point) *
                                JxW[point];
     }
 }
 
 
 template <int dim>
+template <class X, class Y>
 void DGTransportEquation<dim>::assemble_face_term2(
-  const FEFaceValuesBase<dim>& fe_v,
-  const FEFaceValuesBase<dim>& fe_v_neighbor,      
+  const X& fe_v,
+  const Y& fe_v_neighbor,      
   FullMatrix<double> &ui_vi_matrix,
   FullMatrix<double> &ue_vi_matrix,
   FullMatrix<double> &ui_ve_matrix,
   FullMatrix<double> &ue_ve_matrix) const
 {
-  const std::vector<double> &JxW = fe_v.get_JxW_values ();
-  const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
+  const std::vector<double> &JxW = fe_v.get_present_fe_values().get_JxW_values ();
+  const std::vector<Point<dim> > &normals = fe_v.get_present_fe_values().get_normal_vectors ();
 
-  std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
+  std::vector<Point<dim> > beta (fe_v.get_present_fe_values().n_quadrature_points);
   
-  beta_function.value_list (fe_v.get_quadrature_points(), beta);
+  beta_function.value_list (fe_v.get_present_fe_values().get_quadrature_points(), beta);
 
-  for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
+  for (unsigned int point=0; point<fe_v.get_present_fe_values().n_quadrature_points; ++point)
     {
       const double beta_n=beta[point] * normals[point];
       if (beta_n>0)
        {
-         for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
-           for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
+         for (unsigned int i=0; i<fe_v.get_present_fe_values().dofs_per_cell; ++i)
+           for (unsigned int j=0; j<fe_v.get_present_fe_values().dofs_per_cell; ++j)
              ui_vi_matrix(i,j) += beta_n *
-                                fe_v.shape_value(j,point) *
-                                fe_v.shape_value(i,point) *
+                                fe_v.get_present_fe_values().shape_value(j,point) *
+                                fe_v.get_present_fe_values().shape_value(i,point) *
                                 JxW[point];
 
-         for (unsigned int k=0; k<fe_v_neighbor.dofs_per_cell; ++k)
-           for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
+         for (unsigned int k=0; k<fe_v_neighbor.get_present_fe_values().dofs_per_cell; ++k)
+           for (unsigned int j=0; j<fe_v.get_present_fe_values().dofs_per_cell; ++j)
              ui_ve_matrix(k,j) -= beta_n *
-                                 fe_v.shape_value(j,point) *
-                                 fe_v_neighbor.shape_value(k,point) *
+                                 fe_v.get_present_fe_values().shape_value(j,point) *
+                                 fe_v_neighbor.get_present_fe_values().shape_value(k,point) *
                                  JxW[point];
        }
       else
        {
-         for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
-           for (unsigned int l=0; l<fe_v_neighbor.dofs_per_cell; ++l)
+         for (unsigned int i=0; i<fe_v.get_present_fe_values().dofs_per_cell; ++i)
+           for (unsigned int l=0; l<fe_v_neighbor.get_present_fe_values().dofs_per_cell; ++l)
              ue_vi_matrix(i,l) += beta_n *
-                                 fe_v_neighbor.shape_value(l,point) *
-                                 fe_v.shape_value(i,point) *
+                                 fe_v_neighbor.get_present_fe_values().shape_value(l,point) *
+                                 fe_v.get_present_fe_values().shape_value(i,point) *
                                  JxW[point];
 
-         for (unsigned int k=0; k<fe_v_neighbor.dofs_per_cell; ++k)
-           for (unsigned int l=0; l<fe_v_neighbor.dofs_per_cell; ++l)
+         for (unsigned int k=0; k<fe_v_neighbor.get_present_fe_values().dofs_per_cell; ++k)
+           for (unsigned int l=0; l<fe_v_neighbor.get_present_fe_values().dofs_per_cell; ++l)
              ue_ve_matrix(k,l) -= beta_n *
-                                  fe_v_neighbor.shape_value(l,point) *
-                                  fe_v_neighbor.shape_value(k,point) *
+                                  fe_v_neighbor.get_present_fe_values().shape_value(l,point) *
+                                  fe_v_neighbor.get_present_fe_values().shape_value(k,point) *
                                   JxW[point];
        }
     }
@@ -336,16 +340,16 @@ class DGMethod
     void output_results (const unsigned int cycle) const;
     
     Triangulation<dim>   triangulation;
-    const MappingQ1<dim> mapping;
+    const hp::MappingCollection<dim> mapping;
     
-    FE_DGQ<dim>          fe;
-    DoFHandler<dim>      dof_handler;
+    hp::FECollection<dim>          fe;
+    hp::DoFHandler<dim>      dof_handler;
 
     SparsityPattern      sparsity_pattern;
     SparseMatrix<double> system_matrix;
 
-    const QGauss<dim>   quadrature;
-    const QGauss<dim-1> face_quadrature;
+    const hp::QCollection<dim>   quadrature;
+    const hp::QCollection<dim-1> face_quadrature;
     
     Vector<double>       solution1;
     Vector<double>       solution2;
@@ -358,11 +362,11 @@ class DGMethod
 template <int dim>
 DGMethod<dim>::DGMethod ()
                :
-               mapping (),
-                fe (1),
+               mapping (MappingQ1<dim>()),
+                fe (FE_DGQ<dim>(1)),
                dof_handler (triangulation),
-               quadrature (4),
-               face_quadrature (4),
+               quadrature (QGauss<dim>(4)),
+               face_quadrature (QGauss<dim-1>(4)),
                dg ()
 {}
 
@@ -382,7 +386,7 @@ void DGMethod<dim>::setup_system ()
   sparsity_pattern.reinit (dof_handler.n_dofs(),
                           dof_handler.n_dofs(),
                           (GeometryInfo<dim>::faces_per_cell
-                           *GeometryInfo<dim>::subfaces_per_face+1)*fe.dofs_per_cell);
+                           *GeometryInfo<dim>::subfaces_per_face+1)*fe[0].dofs_per_cell);
   
   DoFTools::make_flux_sparsity_pattern (dof_handler, sparsity_pattern);
   
@@ -399,7 +403,7 @@ void DGMethod<dim>::setup_system ()
 template <int dim>
 void DGMethod<dim>::assemble_system1 () 
 {
-  const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
+  const unsigned int dofs_per_cell = dof_handler.get_fe()[0].dofs_per_cell;
   std::vector<unsigned int> dofs (dofs_per_cell);
   std::vector<unsigned int> dofs_neighbor (dofs_per_cell);
 
@@ -415,16 +419,16 @@ void DGMethod<dim>::assemble_system1 ()
   
   const UpdateFlags neighbor_face_update_flags = update_values;
    
-  FEValues<dim> fe_v (
+  hp::FEValues<dim> fe_v (
     mapping, fe, quadrature, update_flags);
   
-  FEFaceValues<dim> fe_v_face (
+  hp::FEFaceValues<dim> fe_v_face (
     mapping, fe, face_quadrature, face_update_flags);
-  FESubfaceValues<dim> fe_v_subface (
+  hp::FESubfaceValues<dim> fe_v_subface (
     mapping, fe, face_quadrature, face_update_flags);
-  FEFaceValues<dim> fe_v_face_neighbor (
+  hp::FEFaceValues<dim> fe_v_face_neighbor (
     mapping, fe, face_quadrature, neighbor_face_update_flags);
-  FESubfaceValues<dim> fe_v_subface_neighbor (
+  hp::FESubfaceValues<dim> fe_v_subface_neighbor (
     mapping, fe, face_quadrature, neighbor_face_update_flags);
 
   FullMatrix<double> ui_vi_matrix (dofs_per_cell, dofs_per_cell);
@@ -432,7 +436,7 @@ void DGMethod<dim>::assemble_system1 ()
 
   Vector<double>  cell_vector (dofs_per_cell);
 
-  typename DoFHandler<dim>::active_cell_iterator
+  typename hp::DoFHandler<dim>::active_cell_iterator
     cell = dof_handler.begin_active(),
     endc = dof_handler.end();
 
@@ -451,7 +455,7 @@ void DGMethod<dim>::assemble_system1 ()
 
       for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
        {
-         typename DoFHandler<dim>::face_iterator face=cell->face(face_no);
+         typename hp::DoFHandler<dim>::face_iterator face=cell->face(face_no);
          
          ue_vi_matrix = 0;
                  
@@ -465,7 +469,7 @@ void DGMethod<dim>::assemble_system1 ()
            }
          else
            {
-             typename DoFHandler<dim>::cell_iterator neighbor=
+             typename hp::DoFHandler<dim>::cell_iterator neighbor=
                cell->neighbor(face_no);;
 
              if (face->has_children())
@@ -477,7 +481,7 @@ void DGMethod<dim>::assemble_system1 ()
                  for (unsigned int subface_no=0;
                       subface_no<face->n_children(); ++subface_no)
                    {
-                     typename DoFHandler<dim>::active_cell_iterator
+                     typename hp::DoFHandler<dim>::active_cell_iterator
                         neighbor_child
                         = cell->neighbor_child_on_subface (face_no, subface_no);
                      
@@ -565,7 +569,7 @@ void DGMethod<dim>::assemble_system1 ()
 template <int dim>
 void DGMethod<dim>::assemble_system2 () 
 {
-  const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
+  const unsigned int dofs_per_cell = dof_handler.get_fe()[0].dofs_per_cell;
   std::vector<unsigned int> dofs (dofs_per_cell);
   std::vector<unsigned int> dofs_neighbor (dofs_per_cell);
 
@@ -581,13 +585,13 @@ void DGMethod<dim>::assemble_system2 ()
   
   const UpdateFlags neighbor_face_update_flags = update_values;
 
-  FEValues<dim> fe_v (
+  hp::FEValues<dim> fe_v (
     mapping, fe, quadrature, update_flags);
-  FEFaceValues<dim> fe_v_face (
+  hp::FEFaceValues<dim> fe_v_face (
     mapping, fe, face_quadrature, face_update_flags);
-  FESubfaceValues<dim> fe_v_subface (
+  hp::FESubfaceValues<dim> fe_v_subface (
     mapping, fe, face_quadrature, face_update_flags);
-  FEFaceValues<dim> fe_v_face_neighbor (
+  hp::FEFaceValues<dim> fe_v_face_neighbor (
     mapping, fe, face_quadrature, neighbor_face_update_flags);
 
 
@@ -599,7 +603,7 @@ void DGMethod<dim>::assemble_system2 ()
   
   Vector<double>  cell_vector (dofs_per_cell);
 
-  typename DoFHandler<dim>::active_cell_iterator
+  typename hp::DoFHandler<dim>::active_cell_iterator
     cell = dof_handler.begin_active(),
     endc = dof_handler.end();
   for (;cell!=endc; ++cell) 
@@ -617,7 +621,7 @@ void DGMethod<dim>::assemble_system2 ()
 
       for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
        {
-         typename DoFHandler<dim>::face_iterator face=
+         typename hp::DoFHandler<dim>::face_iterator face=
            cell->face(face_no);
 
          if (face->at_boundary())
@@ -632,7 +636,7 @@ void DGMethod<dim>::assemble_system2 ()
            {
              Assert (cell->neighbor(face_no).state() == IteratorState::valid,
                      ExcInternalError());
-             typename DoFHandler<dim>::cell_iterator neighbor=
+             typename hp::DoFHandler<dim>::cell_iterator neighbor=
                cell->neighbor(face_no);
              if (face->has_children())
                {
@@ -642,7 +646,7 @@ void DGMethod<dim>::assemble_system2 ()
                  for (unsigned int subface_no=0;
                       subface_no<face->n_children(); ++subface_no)
                    {
-                     typename DoFHandler<dim>::cell_iterator neighbor_child
+                     typename hp::DoFHandler<dim>::cell_iterator neighbor_child
                         = cell->neighbor_child_on_subface (face_no, subface_no);
                      Assert (neighbor_child->face(neighbor2) == face->child(subface_no),
                              ExcInternalError());
@@ -733,7 +737,7 @@ void DGMethod<dim>::solve (Vector<double> &solution)
 
   PreconditionBlockSSOR<SparseMatrix<double> > preconditioner;
 
-  preconditioner.initialize(system_matrix, fe.dofs_per_cell);
+  preconditioner.initialize(system_matrix, fe[0].dofs_per_cell);
 
   solver.solve (system_matrix, solution, right_hand_side,
                preconditioner);
@@ -750,7 +754,7 @@ void DGMethod<dim>::refine_grid ()
                                                 solution2,
                                                 gradient_indicator);
 
-  typename DoFHandler<dim>::active_cell_iterator
+  typename hp::DoFHandler<dim>::active_cell_iterator
     cell = dof_handler.begin_active(),
     endc = dof_handler.end();
   for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no)
@@ -785,7 +789,7 @@ void DGMethod<dim>::output_results (const unsigned int cycle) const
   deallog << "Writing solution to <" << filename << ">..."
            << std::endl << std::endl;
   
-  DataOut<dim> data_out;
+  DataOut<dim,hp::DoFHandler<dim> > data_out;
   data_out.attach_dof_handler (dof_handler);
   data_out.add_data_vector (solution2, "u");
 
index 2d80eb1603d8bfaa5b21839a7e13e01bb0b67d5b..8b8aea97c2040b17bd1d11ce904795b625d80ad1 100644 (file)
@@ -40,7 +40,7 @@ std::ofstream logfile("step-13/output");
 #include <dofs/dof_accessor.h>
 #include <dofs/dof_tools.h>
 #include <fe/fe_q.h>
-#include <fe/fe_values.h>
+#include <fe/hp_fe_values.h>
 #include <numerics/vectors.h>
 #include <numerics/matrices.h>
 #include <numerics/data_out.h>
@@ -64,7 +64,7 @@ namespace Evaluation
 
       void set_refinement_cycle (const unsigned int refinement_cycle);
       
-      virtual void operator () (const DoFHandler<dim> &dof_handler,
+      virtual void operator () (const hp::DoFHandler<dim> &dof_handler,
                                const Vector<double>  &solution) const = 0;
     protected:
       unsigned int refinement_cycle;
@@ -93,7 +93,7 @@ namespace Evaluation
       PointValueEvaluation (const Point<dim>   &evaluation_point,
                            TableHandler       &results_table);
       
-      virtual void operator () (const DoFHandler<dim> &dof_handler,
+      virtual void operator () (const hp::DoFHandler<dim> &dof_handler,
                                const Vector<double>  &solution) const;
       
       DeclException1 (ExcEvaluationPointNotFound,
@@ -120,12 +120,12 @@ namespace Evaluation
   template <int dim>
   void
   PointValueEvaluation<dim>::
-  operator () (const DoFHandler<dim> &dof_handler,
+  operator () (const hp::DoFHandler<dim> &dof_handler,
               const Vector<double>  &solution) const 
   {
     double point_value = 1e20;
 
-    typename DoFHandler<dim>::active_cell_iterator
+    typename hp::DoFHandler<dim>::active_cell_iterator
       cell = dof_handler.begin_active(),
       endc = dof_handler.end();
     bool evaluation_point_found = false;
@@ -135,7 +135,8 @@ namespace Evaluation
           ++vertex)
        if (cell->vertex(vertex) == evaluation_point)
          {
-           point_value = solution(cell->vertex_dof_index(vertex,0));
+           point_value = solution(cell->vertex_dof_index(vertex,0,
+                                                         cell->active_fe_index()));
 
            evaluation_point_found = true;
            break;
@@ -159,7 +160,7 @@ namespace Evaluation
       SolutionOutput (const std::string                         &output_name_base,
                      const typename DataOut<dim>::OutputFormat  output_format);
       
-      virtual void operator () (const DoFHandler<dim> &dof_handler,
+      virtual void operator () (const hp::DoFHandler<dim> &dof_handler,
                                const Vector<double>  &solution) const;
     private:
       const std::string                         output_name_base;
@@ -179,10 +180,10 @@ namespace Evaluation
 
   template <int dim>
   void
-  SolutionOutput<dim>::operator () (const DoFHandler<dim> &dof_handler,
+  SolutionOutput<dim>::operator () (const hp::DoFHandler<dim> &dof_handler,
                                    const Vector<double>  &solution) const
   {
-    DataOut<dim> data_out;
+    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 ();
@@ -241,8 +242,8 @@ namespace LaplaceSolver
   {
     public:
       Solver (Triangulation<dim>       &triangulation,
-             const FiniteElement<dim> &fe,
-             const Quadrature<dim>    &quadrature,
+             const hp::FECollection<dim> &fe,
+             const hp::QCollection<dim>    &quadrature,
              const Function<dim>      &boundary_values);
       virtual
       ~Solver ();
@@ -260,9 +261,9 @@ namespace LaplaceSolver
       n_dofs () const;
       
     protected:
-      const SmartPointer<const FiniteElement<dim> >  fe;
-      const SmartPointer<const Quadrature<dim> >     quadrature;
-      DoFHandler<dim>                                dof_handler;
+      const SmartPointer<const hp::FECollection<dim> >  fe;
+      const SmartPointer<const hp::QCollection<dim> >     quadrature;
+      hp::DoFHandler<dim>                                dof_handler;
       Vector<double>                                 solution;
       const SmartPointer<const Function<dim> >       boundary_values;
 
@@ -271,7 +272,7 @@ namespace LaplaceSolver
     private:
       struct LinearSystem
       {
-         LinearSystem (const DoFHandler<dim> &dof_handler);
+         LinearSystem (const hp::DoFHandler<dim> &dof_handler);
 
          void solve (Vector<double> &solution) const;
        
@@ -286,8 +287,8 @@ namespace LaplaceSolver
 
       void
       assemble_matrix (LinearSystem                                         &linear_system,
-                      const typename DoFHandler<dim>::active_cell_iterator &begin_cell,
-                      const typename DoFHandler<dim>::active_cell_iterator &end_cell,
+                      const typename hp::DoFHandler<dim>::active_cell_iterator &begin_cell,
+                      const typename hp::DoFHandler<dim>::active_cell_iterator &end_cell,
                       Threads::ThreadMutex                                 &mutex) const;
   };
 
@@ -295,8 +296,8 @@ namespace LaplaceSolver
 
   template <int dim>
   Solver<dim>::Solver (Triangulation<dim>       &triangulation,
-                      const FiniteElement<dim> &fe,
-                      const Quadrature<dim>    &quadrature,
+                      const hp::FECollection<dim> &fe,
+                      const hp::QCollection<dim>    &quadrature,
                       const Function<dim>      &boundary_values)
                  :
                  Base<dim> (triangulation),
@@ -349,7 +350,7 @@ namespace LaplaceSolver
   Solver<dim>::assemble_linear_system (LinearSystem &linear_system)
   {
     typedef
-      typename DoFHandler<dim>::active_cell_iterator
+      typename hp::DoFHandler<dim>::active_cell_iterator
       active_cell_iterator;
 
     const unsigned int n_threads = multithread_info.n_default_threads;
@@ -392,21 +393,21 @@ namespace LaplaceSolver
   template <int dim>
   void
   Solver<dim>::assemble_matrix (LinearSystem                                         &linear_system,
-                               const typename DoFHandler<dim>::active_cell_iterator &begin_cell,
-                               const typename DoFHandler<dim>::active_cell_iterator &end_cell,
+                               const typename hp::DoFHandler<dim>::active_cell_iterator &begin_cell,
+                               const typename hp::DoFHandler<dim>::active_cell_iterator &end_cell,
                                Threads::ThreadMutex                                 &mutex) const
   {
-    FEValues<dim> fe_values (*fe, *quadrature, 
+    hp::FEValues<dim> fe_values (*fe, *quadrature, 
                             update_gradients | update_JxW_values);
 
-    const unsigned int   dofs_per_cell = fe->dofs_per_cell;
-    const unsigned int   n_q_points    = quadrature->n_quadrature_points;
+    const unsigned int   dofs_per_cell = (*fe)[0].dofs_per_cell;
+    const unsigned int   n_q_points    = (*quadrature)[0].n_quadrature_points;
 
     FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
 
     std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 
-    for (typename DoFHandler<dim>::active_cell_iterator cell=begin_cell;
+    for (typename hp::DoFHandler<dim>::active_cell_iterator cell=begin_cell;
         cell!=end_cell; ++cell)
       {
        cell_matrix = 0;
@@ -416,9 +417,9 @@ namespace LaplaceSolver
        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.JxW(q_point));
+             cell_matrix(i,j) += (fe_values.get_present_fe_values().shape_grad(i,q_point) *
+                                  fe_values.get_present_fe_values().shape_grad(j,q_point) *
+                                  fe_values.get_present_fe_values().JxW(q_point));
 
 
        cell->get_dof_indices (local_dof_indices);
@@ -435,11 +436,11 @@ namespace LaplaceSolver
 
   template <int dim>
   Solver<dim>::LinearSystem::
-  LinearSystem (const DoFHandler<dim> &dof_handler)
+  LinearSystem (const hp::DoFHandler<dim> &dof_handler)
   {
     hanging_node_constraints.clear ();
 
-    void (*mhnc_p) (const DoFHandler<dim> &,
+    void (*mhnc_p) (const hp::DoFHandler<dim> &,
                    ConstraintMatrix      &)
       = &DoFTools::make_hanging_node_constraints;
     
@@ -487,8 +488,8 @@ namespace LaplaceSolver
   {
     public:
       PrimalSolver (Triangulation<dim>       &triangulation,
-                   const FiniteElement<dim> &fe,
-                   const Quadrature<dim>    &quadrature,
+                   const hp::FECollection<dim> &fe,
+                   const hp::QCollection<dim>    &quadrature,
                    const Function<dim>      &rhs_function,
                    const Function<dim>      &boundary_values);
     protected:
@@ -500,8 +501,8 @@ namespace LaplaceSolver
   template <int dim>
   PrimalSolver<dim>::
   PrimalSolver (Triangulation<dim>       &triangulation,
-               const FiniteElement<dim> &fe,
-               const Quadrature<dim>    &quadrature,
+               const hp::FECollection<dim> &fe,
+               const hp::QCollection<dim>    &quadrature,
                const Function<dim>      &rhs_function,
                const Function<dim>      &boundary_values)
                  :
@@ -518,32 +519,32 @@ namespace LaplaceSolver
   PrimalSolver<dim>::
   assemble_rhs (Vector<double> &rhs) const 
   {
-    FEValues<dim> fe_values (*this->fe, *this->quadrature, 
+    hp::FEValues<dim> fe_values (*this->fe, *this->quadrature, 
                             update_values | update_q_points  |
                              update_JxW_values);
 
-    const unsigned int   dofs_per_cell = this->fe->dofs_per_cell;
-    const unsigned int   n_q_points    = this->quadrature->n_quadrature_points;
+    const unsigned int   dofs_per_cell = (*this->fe)[0].dofs_per_cell;
+    const unsigned int   n_q_points    = (*this->quadrature)[0].n_quadrature_points;
 
     Vector<double>       cell_rhs (dofs_per_cell);
     std::vector<double>  rhs_values (n_q_points);
     std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 
-    typename DoFHandler<dim>::active_cell_iterator
+    typename hp::DoFHandler<dim>::active_cell_iterator
       cell = this->dof_handler.begin_active(),
       endc = this->dof_handler.end();
     for (; cell!=endc; ++cell)
       {
        cell_rhs = 0;
        fe_values.reinit (cell);
-       rhs_function->value_list (fe_values.get_quadrature_points(),
+       rhs_function->value_list (fe_values.get_present_fe_values().get_quadrature_points(),
                                  rhs_values);
       
        for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
          for (unsigned int i=0; i<dofs_per_cell; ++i)
-           cell_rhs(i) += (fe_values.shape_value(i,q_point) *
+           cell_rhs(i) += (fe_values.get_present_fe_values().shape_value(i,q_point) *
                            rhs_values[q_point] *
-                           fe_values.JxW(q_point));
+                           fe_values.get_present_fe_values().JxW(q_point));
 
        cell->get_dof_indices (local_dof_indices);
        for (unsigned int i=0; i<dofs_per_cell; ++i)
@@ -558,8 +559,8 @@ namespace LaplaceSolver
   {
     public:
       RefinementGlobal (Triangulation<dim>       &coarse_grid,
-                       const FiniteElement<dim> &fe,
-                       const Quadrature<dim>    &quadrature,
+                       const hp::FECollection<dim> &fe,
+                       const hp::QCollection<dim>    &quadrature,
                        const Function<dim>      &rhs_function,
                        const Function<dim>      &boundary_values);
 
@@ -571,8 +572,8 @@ namespace LaplaceSolver
   template <int dim>
   RefinementGlobal<dim>::
   RefinementGlobal (Triangulation<dim>       &coarse_grid,
-                   const FiniteElement<dim> &fe,
-                   const Quadrature<dim>    &quadrature,
+                   const hp::FECollection<dim> &fe,
+                   const hp::QCollection<dim>    &quadrature,
                    const Function<dim>      &rhs_function,
                    const Function<dim>      &boundary_values)
                  :
@@ -597,8 +598,8 @@ namespace LaplaceSolver
   {
     public:
       RefinementKelly (Triangulation<dim>       &coarse_grid,
-                      const FiniteElement<dim> &fe,
-                      const Quadrature<dim>    &quadrature,
+                      const hp::FECollection<dim> &fe,
+                      const hp::QCollection<dim>    &quadrature,
                       const Function<dim>      &rhs_function,
                       const Function<dim>      &boundary_values);
 
@@ -610,8 +611,8 @@ namespace LaplaceSolver
   template <int dim>
   RefinementKelly<dim>::
   RefinementKelly (Triangulation<dim>       &coarse_grid,
-                  const FiniteElement<dim> &fe,
-                  const Quadrature<dim>    &quadrature,
+                  const hp::FECollection<dim> &fe,
+                  const hp::QCollection<dim>    &quadrature,
                   const Function<dim>      &rhs_function,
                   const Function<dim>      &boundary_values)
                  :
@@ -752,8 +753,8 @@ void solve_problem (const std::string &solver_name)
   Triangulation<dim> triangulation;
   GridGenerator::hyper_cube (triangulation, -1, 1);
   triangulation.refine_global (2);
-  const FE_Q<dim>          fe(1);
-  const QGauss<dim>       quadrature(4);
+  const hp::FECollection<dim>          fe(FE_Q<dim>(1));
+  const hp::QCollection<dim>       quadrature(QGauss<dim>(4));
   const RightHandSide<dim> rhs_function;
   const Solution<dim>      boundary_values;
 
index 725cb2cbe2aaff80cc971cca5165628e2379e9d7..3c99840ca2124e9b357625dcdda8732f3e515ebc 100644 (file)
@@ -40,7 +40,7 @@ std::ofstream logfile("step-14/output");
 #include <dofs/dof_accessor.h>
 #include <dofs/dof_tools.h>
 #include <fe/fe_q.h>
-#include <fe/fe_values.h>
+#include <fe/hp_fe_values.h>
 #include <fe/fe_tools.h>
 #include <numerics/vectors.h>
 #include <numerics/matrices.h>
@@ -65,7 +65,7 @@ namespace Evaluation
 
       void set_refinement_cycle (const unsigned int refinement_cycle);
       
-      virtual void operator () (const DoFHandler<dim> &dof_handler,
+      virtual void operator () (const hp::DoFHandler<dim> &dof_handler,
                                const Vector<double>  &solution) const = 0;
     protected:
       unsigned int refinement_cycle;
@@ -92,7 +92,7 @@ namespace Evaluation
     public:
       PointValueEvaluation (const Point<dim>   &evaluation_point);
       
-      virtual void operator () (const DoFHandler<dim> &dof_handler,
+      virtual void operator () (const hp::DoFHandler<dim> &dof_handler,
                                const Vector<double>  &solution) const;
       
       DeclException1 (ExcEvaluationPointNotFound,
@@ -116,12 +116,12 @@ namespace Evaluation
   template <int dim>
   void
   PointValueEvaluation<dim>::
-  operator () (const DoFHandler<dim> &dof_handler,
+  operator () (const hp::DoFHandler<dim> &dof_handler,
               const Vector<double>  &solution) const 
   {
     double point_value = 1e20;
 
-    typename DoFHandler<dim>::active_cell_iterator
+    typename hp::DoFHandler<dim>::active_cell_iterator
       cell = dof_handler.begin_active(),
       endc = dof_handler.end();
     bool evaluation_point_found = false;
@@ -154,7 +154,7 @@ namespace Evaluation
     public:
       PointXDerivativeEvaluation (const Point<dim>   &evaluation_point);
       
-      virtual void operator () (const DoFHandler<dim> &dof_handler,
+      virtual void operator () (const hp::DoFHandler<dim> &dof_handler,
                                const Vector<double>  &solution) const;
       
       DeclException1 (ExcEvaluationPointNotFound,
@@ -177,19 +177,19 @@ namespace Evaluation
   template <int dim>
   void
   PointXDerivativeEvaluation<dim>::
-  operator () (const DoFHandler<dim> &dof_handler,
+  operator () (const hp::DoFHandler<dim> &dof_handler,
               const Vector<double>  &solution) const 
   {
     double point_derivative = 0;
 
     QTrapez<dim>  vertex_quadrature;
-    FEValues<dim> fe_values (dof_handler.get_fe(),
+    hp::FEValues<dim> fe_values (dof_handler.get_fe(),
                             vertex_quadrature,
                             update_gradients | update_q_points);
     std::vector<Tensor<1,dim> >
       solution_gradients (vertex_quadrature.n_quadrature_points);
     
-    typename DoFHandler<dim>::active_cell_iterator
+    typename hp::DoFHandler<dim>::active_cell_iterator
       cell = dof_handler.begin_active(),
       endc = dof_handler.end();
     unsigned int evaluation_point_hits = 0;
@@ -200,12 +200,12 @@ namespace Evaluation
        if (cell->vertex(vertex) == evaluation_point)
          {
            fe_values.reinit (cell);
-           fe_values.get_function_grads (solution,
+           fe_values.get_present_fe_values().get_function_grads (solution,
                                          solution_gradients);
 
            unsigned int q_point = 0;
            for (; q_point<solution_gradients.size(); ++q_point)
-             if (fe_values.quadrature_point(q_point) ==
+             if (fe_values.get_present_fe_values().quadrature_point(q_point) ==
                  evaluation_point)
                break;
 
@@ -234,7 +234,7 @@ namespace Evaluation
     public:
       GridOutput (const std::string &output_name_base);
       
-      virtual void operator () (const DoFHandler<dim> &dof_handler,
+      virtual void operator () (const hp::DoFHandler<dim> &dof_handler,
                                const Vector<double>  &solution) const;
     private:
       const std::string output_name_base;
@@ -251,7 +251,7 @@ namespace Evaluation
 
   template <int dim>
   void
-  GridOutput<dim>::operator () (const DoFHandler<dim> &dof_handler,
+  GridOutput<dim>::operator () (const hp::DoFHandler<dim> &dof_handler,
                                const Vector<double>  &/*solution*/) const
   {
     std::ostringstream filename;
@@ -323,9 +323,9 @@ namespace LaplaceSolver
   {
     public:
       Solver (Triangulation<dim>       &triangulation,
-             const FiniteElement<dim> &fe,
-             const Quadrature<dim>    &quadrature,
-             const Quadrature<dim-1>  &face_quadrature,              
+             const hp::FECollection<dim> &fe,
+             const hp::QCollection<dim>    &quadrature,
+             const hp::QCollection<dim-1>  &face_quadrature,         
              const Function<dim>      &boundary_values);
       virtual
       ~Solver ();
@@ -343,10 +343,10 @@ namespace LaplaceSolver
       n_dofs () const;
       
     protected:
-      const SmartPointer<const FiniteElement<dim> >  fe;
-      const SmartPointer<const Quadrature<dim> >     quadrature;
-      const SmartPointer<const Quadrature<dim-1> >   face_quadrature;      
-      DoFHandler<dim>                                dof_handler;
+      const SmartPointer<const hp::FECollection<dim> >  fe;
+      const SmartPointer<const hp::QCollection<dim> >     quadrature;
+      const SmartPointer<const hp::QCollection<dim-1> >   face_quadrature;      
+      hp::DoFHandler<dim>                                dof_handler;
       Vector<double>                                 solution;
       const SmartPointer<const Function<dim> >       boundary_values;
 
@@ -355,7 +355,7 @@ namespace LaplaceSolver
     private:
       struct LinearSystem
       {
-         LinearSystem (const DoFHandler<dim> &dof_handler);
+         LinearSystem (const hp::DoFHandler<dim> &dof_handler);
 
          void solve (Vector<double> &solution) const;
        
@@ -370,8 +370,8 @@ namespace LaplaceSolver
 
       void
       assemble_matrix (LinearSystem                                         &linear_system,
-                      const typename DoFHandler<dim>::active_cell_iterator &begin_cell,
-                      const typename DoFHandler<dim>::active_cell_iterator &end_cell,
+                      const typename hp::DoFHandler<dim>::active_cell_iterator &begin_cell,
+                      const typename hp::DoFHandler<dim>::active_cell_iterator &end_cell,
                       Threads::ThreadMutex                                 &mutex) const;
   };
 
@@ -379,9 +379,9 @@ namespace LaplaceSolver
 
   template <int dim>
   Solver<dim>::Solver (Triangulation<dim>       &triangulation,
-                      const FiniteElement<dim> &fe,
-                      const Quadrature<dim>    &quadrature,
-                      const Quadrature<dim-1>  &face_quadrature,
+                      const hp::FECollection<dim> &fe,
+                      const hp::QCollection<dim>    &quadrature,
+                      const hp::QCollection<dim-1>  &face_quadrature,
                       const Function<dim>      &boundary_values)
                  :
                  Base<dim> (triangulation),
@@ -435,7 +435,7 @@ namespace LaplaceSolver
   Solver<dim>::assemble_linear_system (LinearSystem &linear_system)
   {
     typedef
-      typename DoFHandler<dim>::active_cell_iterator
+      typename hp::DoFHandler<dim>::active_cell_iterator
       active_cell_iterator;
 
     const unsigned int n_threads = multithread_info.n_default_threads;
@@ -476,21 +476,21 @@ namespace LaplaceSolver
   template <int dim>
   void
   Solver<dim>::assemble_matrix (LinearSystem                                         &linear_system,
-                               const typename DoFHandler<dim>::active_cell_iterator &begin_cell,
-                               const typename DoFHandler<dim>::active_cell_iterator &end_cell,
+                               const typename hp::DoFHandler<dim>::active_cell_iterator &begin_cell,
+                               const typename hp::DoFHandler<dim>::active_cell_iterator &end_cell,
                                Threads::ThreadMutex                                 &mutex) const
   {
-    FEValues<dim> fe_values (*fe, *quadrature, 
+    hp::FEValues<dim> fe_values (*fe, *quadrature, 
                             update_gradients | update_JxW_values);
 
-    const unsigned int   dofs_per_cell = fe->dofs_per_cell;
-    const unsigned int   n_q_points    = quadrature->n_quadrature_points;
+    const unsigned int   dofs_per_cell = (*fe)[0].dofs_per_cell;
+    const unsigned int   n_q_points    = (*quadrature)[0].n_quadrature_points;
 
     FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
 
     std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 
-    for (typename DoFHandler<dim>::active_cell_iterator cell=begin_cell;
+    for (typename hp::DoFHandler<dim>::active_cell_iterator cell=begin_cell;
         cell!=end_cell; ++cell)
       {
        cell_matrix = 0;
@@ -500,9 +500,9 @@ namespace LaplaceSolver
        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.JxW(q_point));
+             cell_matrix(i,j) += (fe_values.get_present_fe_values().shape_grad(i,q_point) *
+                                  fe_values.get_present_fe_values().shape_grad(j,q_point) *
+                                  fe_values.get_present_fe_values().JxW(q_point));
 
 
        cell->get_dof_indices (local_dof_indices);
@@ -518,11 +518,11 @@ namespace LaplaceSolver
 
   template <int dim>
   Solver<dim>::LinearSystem::
-  LinearSystem (const DoFHandler<dim> &dof_handler)
+  LinearSystem (const hp::DoFHandler<dim> &dof_handler)
   {
     hanging_node_constraints.clear ();
 
-    void (*mhnc_p) (const DoFHandler<dim> &,
+    void (*mhnc_p) (const hp::DoFHandler<dim> &,
                    ConstraintMatrix      &)
       = &DoFTools::make_hanging_node_constraints;
     
@@ -569,9 +569,9 @@ namespace LaplaceSolver
   {
     public:
       PrimalSolver (Triangulation<dim>       &triangulation,
-                   const FiniteElement<dim> &fe,
-                   const Quadrature<dim>    &quadrature,
-                   const Quadrature<dim-1>  &face_quadrature,
+                   const hp::FECollection<dim> &fe,
+                   const hp::QCollection<dim>    &quadrature,
+                   const hp::QCollection<dim-1>  &face_quadrature,
                    const Function<dim>      &rhs_function,
                    const Function<dim>      &boundary_values);
 
@@ -598,9 +598,9 @@ namespace LaplaceSolver
   template <int dim>
   PrimalSolver<dim>::
   PrimalSolver (Triangulation<dim>       &triangulation,
-               const FiniteElement<dim> &fe,
-               const Quadrature<dim>    &quadrature,
-               const Quadrature<dim-1>  &face_quadrature,
+               const hp::FECollection<dim> &fe,
+               const hp::QCollection<dim>    &quadrature,
+               const hp::QCollection<dim-1>  &face_quadrature,
                const Function<dim>      &rhs_function,
                const Function<dim>      &boundary_values)
                  :
@@ -642,7 +642,7 @@ namespace LaplaceSolver
   void
   PrimalSolver<dim>::output_solution () const
   {
-    DataOut<dim> data_out;
+    DataOut<dim,hp::DoFHandler<dim> > data_out;
     data_out.attach_dof_handler (this->dof_handler);
     data_out.add_data_vector (this->solution, "solution");
     data_out.build_patches ();
@@ -662,18 +662,18 @@ namespace LaplaceSolver
   PrimalSolver<dim>::
   assemble_rhs (Vector<double> &rhs) const 
   {
-    FEValues<dim> fe_values (*this->fe, *this->quadrature, 
+    hp::FEValues<dim> fe_values (*this->fe, *this->quadrature, 
                             update_values  | update_q_points  |
                              update_JxW_values);
 
-    const unsigned int   dofs_per_cell = this->fe->dofs_per_cell;
-    const unsigned int   n_q_points    = this->quadrature->n_quadrature_points;
+    const unsigned int   dofs_per_cell = (*this->fe)[0].dofs_per_cell;
+    const unsigned int   n_q_points    = (*this->quadrature)[0].n_quadrature_points;
 
     Vector<double>       cell_rhs (dofs_per_cell);
     std::vector<double>  rhs_values (n_q_points);
     std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 
-    typename DoFHandler<dim>::active_cell_iterator
+    typename hp::DoFHandler<dim>::active_cell_iterator
       cell = this->dof_handler.begin_active(),
       endc = this->dof_handler.end();
     for (; cell!=endc; ++cell)
@@ -682,14 +682,14 @@ namespace LaplaceSolver
 
        fe_values.reinit (cell);
 
-       rhs_function->value_list (fe_values.get_quadrature_points(),
+       rhs_function->value_list (fe_values.get_present_fe_values().get_quadrature_points(),
                                  rhs_values);
       
        for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
          for (unsigned int i=0; i<dofs_per_cell; ++i)
-           cell_rhs(i) += (fe_values.shape_value(i,q_point) *
+           cell_rhs(i) += (fe_values.get_present_fe_values().shape_value(i,q_point) *
                            rhs_values[q_point] *
-                           fe_values.JxW(q_point));
+                           fe_values.get_present_fe_values().JxW(q_point));
 
        cell->get_dof_indices (local_dof_indices);
        for (unsigned int i=0; i<dofs_per_cell; ++i)
@@ -704,9 +704,9 @@ namespace LaplaceSolver
   {
     public:
       RefinementGlobal (Triangulation<dim>       &coarse_grid,
-                       const FiniteElement<dim> &fe,
-                       const Quadrature<dim>    &quadrature,
-                       const Quadrature<dim-1>  &face_quadrature,
+                       const hp::FECollection<dim> &fe,
+                       const hp::QCollection<dim>    &quadrature,
+                       const hp::QCollection<dim-1>  &face_quadrature,
                        const Function<dim>      &rhs_function,
                        const Function<dim>      &boundary_values);
 
@@ -718,9 +718,9 @@ namespace LaplaceSolver
   template <int dim>
   RefinementGlobal<dim>::
   RefinementGlobal (Triangulation<dim>       &coarse_grid,
-                   const FiniteElement<dim> &fe,
-                   const Quadrature<dim>    &quadrature,
-                   const Quadrature<dim-1>  &face_quadrature,
+                   const hp::FECollection<dim> &fe,
+                   const hp::QCollection<dim>    &quadrature,
+                   const hp::QCollection<dim-1>  &face_quadrature,
                    const Function<dim>      &rhs_function,
                    const Function<dim>      &boundary_values)
                  :
@@ -746,9 +746,9 @@ namespace LaplaceSolver
   {
     public:
       RefinementKelly (Triangulation<dim>       &coarse_grid,
-                      const FiniteElement<dim> &fe,
-                      const Quadrature<dim>    &quadrature,
-                      const Quadrature<dim-1>  &face_quadrature,
+                      const hp::FECollection<dim> &fe,
+                      const hp::QCollection<dim>    &quadrature,
+                      const hp::QCollection<dim-1>  &face_quadrature,
                       const Function<dim>      &rhs_function,
                       const Function<dim>      &boundary_values);
 
@@ -760,9 +760,9 @@ namespace LaplaceSolver
   template <int dim>
   RefinementKelly<dim>::
   RefinementKelly (Triangulation<dim>       &coarse_grid,
-                  const FiniteElement<dim> &fe,
-                  const Quadrature<dim>    &quadrature,
-                  const Quadrature<dim-1>  &face_quadrature,
+                  const hp::FECollection<dim> &fe,
+                  const hp::QCollection<dim>    &quadrature,
+                  const hp::QCollection<dim-1>  &face_quadrature,
                   const Function<dim>      &rhs_function,
                   const Function<dim>      &boundary_values)
                  :
@@ -798,9 +798,9 @@ namespace LaplaceSolver
   {
     public:
       RefinementWeightedKelly (Triangulation<dim>       &coarse_grid,
-                              const FiniteElement<dim> &fe,
-                              const Quadrature<dim>    &quadrature,
-                              const Quadrature<dim-1>  &face_quadrature,
+                              const hp::FECollection<dim> &fe,
+                              const hp::QCollection<dim>    &quadrature,
+                              const hp::QCollection<dim-1>  &face_quadrature,
                               const Function<dim>      &rhs_function,
                               const Function<dim>      &boundary_values,
                               const Function<dim>      &weighting_function);
@@ -816,9 +816,9 @@ namespace LaplaceSolver
   template <int dim>
   RefinementWeightedKelly<dim>::
   RefinementWeightedKelly (Triangulation<dim>       &coarse_grid,
-                          const FiniteElement<dim> &fe,
-                          const Quadrature<dim>    &quadrature,
-                          const Quadrature<dim-1>  &face_quadrature,
+                          const hp::FECollection<dim> &fe,
+                          const hp::QCollection<dim>    &quadrature,
+                          const hp::QCollection<dim-1>  &face_quadrature,
                           const Function<dim>      &rhs_function,
                           const Function<dim>      &boundary_values,
                           const Function<dim>      &weighting_function)
@@ -838,12 +838,12 @@ namespace LaplaceSolver
   {
     Vector<float> estimated_error (this->triangulation->n_active_cells());
     KellyErrorEstimator<dim>::estimate (this->dof_handler,
-                                       *this->face_quadrature,
+                                       (*this->face_quadrature)[0],
                                        typename FunctionMap<dim>::type(),
                                        this->solution,
                                        estimated_error);
 
-    typename DoFHandler<dim>::active_cell_iterator
+    typename hp::DoFHandler<dim>::active_cell_iterator
       cell = this->dof_handler.begin_active(),
       endc = this->dof_handler.end();
     for (unsigned int cell_index=0; cell!=endc; ++cell, ++cell_index)
@@ -1112,7 +1112,7 @@ namespace DualFunctional
     public:
       virtual
       void
-      assemble_rhs (const DoFHandler<dim> &dof_handler,
+      assemble_rhs (const hp::DoFHandler<dim> &dof_handler,
                    Vector<double>        &rhs) const = 0;
   };
 
@@ -1126,7 +1126,7 @@ namespace DualFunctional
 
       virtual
       void
-      assemble_rhs (const DoFHandler<dim> &dof_handler,
+      assemble_rhs (const hp::DoFHandler<dim> &dof_handler,
                    Vector<double>        &rhs) const;
       
       DeclException1 (ExcEvaluationPointNotFound,
@@ -1150,12 +1150,12 @@ namespace DualFunctional
   template <int dim>
   void
   PointValueEvaluation<dim>::
-  assemble_rhs (const DoFHandler<dim> &dof_handler,
+  assemble_rhs (const hp::DoFHandler<dim> &dof_handler,
                Vector<double>        &rhs) const
   {
     rhs.reinit (dof_handler.n_dofs());
 
-    typename DoFHandler<dim>::active_cell_iterator
+    typename hp::DoFHandler<dim>::active_cell_iterator
       cell = dof_handler.begin_active(),
       endc = dof_handler.end();
     for (; cell!=endc; ++cell)
@@ -1182,7 +1182,7 @@ namespace DualFunctional
 
       virtual
       void
-      assemble_rhs (const DoFHandler<dim> &dof_handler,
+      assemble_rhs (const hp::DoFHandler<dim> &dof_handler,
                    Vector<double>        &rhs) const;
       
       DeclException1 (ExcEvaluationPointNotFound,
@@ -1206,17 +1206,17 @@ namespace DualFunctional
   template <int dim>
   void
   PointXDerivativeEvaluation<dim>::
-  assemble_rhs (const DoFHandler<dim> &dof_handler,
+  assemble_rhs (const hp::DoFHandler<dim> &dof_handler,
                Vector<double>        &rhs) const
   {
     rhs.reinit (dof_handler.n_dofs());
 
     QGauss<dim> quadrature(4);
-    FEValues<dim>  fe_values (dof_handler.get_fe(), quadrature,
+    hp::FEValues<dim>  fe_values (dof_handler.get_fe(), quadrature,
                              update_gradients |
                              update_q_points  |
                              update_JxW_values);
-    const unsigned int n_q_points = fe_values.n_quadrature_points;
+    const unsigned int n_q_points = fe_values.get_present_fe_values().n_quadrature_points;
     const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
 
     Vector<double> cell_rhs (dofs_per_cell);
@@ -1224,7 +1224,7 @@ namespace DualFunctional
 
     double total_volume = 0;
     
-    typename DoFHandler<dim>::active_cell_iterator
+    typename hp::DoFHandler<dim>::active_cell_iterator
       cell = dof_handler.begin_active(),
       endc = dof_handler.end();
     for (; cell!=endc; ++cell)
@@ -1237,9 +1237,9 @@ namespace DualFunctional
          for (unsigned int q=0; q<n_q_points; ++q)
            {
              for (unsigned int i=0; i<dofs_per_cell; ++i)
-               cell_rhs(i) += fe_values.shape_grad(i,q)[0] *
-                              fe_values.JxW (q);
-             total_volume += fe_values.JxW (q);
+               cell_rhs(i) += fe_values.get_present_fe_values().shape_grad(i,q)[0] *
+                              fe_values.get_present_fe_values().JxW (q);
+             total_volume += fe_values.get_present_fe_values().JxW (q);
            };
 
          cell->get_dof_indices (local_dof_indices);
@@ -1266,9 +1266,9 @@ namespace LaplaceSolver
   {
     public:
       DualSolver (Triangulation<dim>       &triangulation,
-                 const FiniteElement<dim> &fe,
-                 const Quadrature<dim>    &quadrature,
-                 const Quadrature<dim-1>  &face_quadrature,
+                 const hp::FECollection<dim> &fe,
+                 const hp::QCollection<dim>    &quadrature,
+                 const hp::QCollection<dim-1>  &face_quadrature,
                  const DualFunctional::DualFunctionalBase<dim> &dual_functional);
 
       virtual
@@ -1298,9 +1298,9 @@ namespace LaplaceSolver
   template <int dim>
   DualSolver<dim>::
   DualSolver (Triangulation<dim>       &triangulation,
-             const FiniteElement<dim> &fe,
-             const Quadrature<dim>    &quadrature,
-             const Quadrature<dim-1>  &face_quadrature,
+             const hp::FECollection<dim> &fe,
+             const hp::QCollection<dim>    &quadrature,
+             const hp::QCollection<dim-1>  &face_quadrature,
              const DualFunctional::DualFunctionalBase<dim> &dual_functional)
                  :
                  Base<dim> (triangulation),
@@ -1354,10 +1354,10 @@ namespace LaplaceSolver
   {
     public:
       WeightedResidual (Triangulation<dim>       &coarse_grid,
-                       const FiniteElement<dim> &primal_fe,
-                       const FiniteElement<dim> &dual_fe,
-                       const Quadrature<dim>    &quadrature,
-                       const Quadrature<dim-1>  &face_quadrature,
+                       const hp::FECollection<dim> &primal_fe,
+                       const hp::FECollection<dim> &dual_fe,
+                       const hp::QCollection<dim>    &quadrature,
+                       const hp::QCollection<dim-1>  &face_quadrature,
                        const Function<dim>      &rhs_function,
                        const Function<dim>      &boundary_values,
                        const DualFunctional::DualFunctionalBase<dim> &dual_functional);
@@ -1385,39 +1385,39 @@ namespace LaplaceSolver
       void solve_dual_problem ();
                                        
       typedef
-      typename DoFHandler<dim>::active_cell_iterator
+      typename hp::DoFHandler<dim>::active_cell_iterator
       active_cell_iterator;
 
       typedef
-      typename std::map<typename DoFHandler<dim>::face_iterator,double>
+      typename std::map<typename hp::DoFHandler<dim>::face_iterator,double>
       FaceIntegrals;
 
       struct CellData
       {
-         FEValues<dim>    fe_values;
+         hp::FEValues<dim>    fe_values;
          const SmartPointer<const Function<dim> > right_hand_side;
 
          std::vector<double> cell_residual;
          std::vector<double> rhs_values;         
          std::vector<double> dual_weights;       
          typename std::vector<Tensor<2,dim> > cell_grad_grads;
-         CellData (const FiniteElement<dim> &fe,
-                   const Quadrature<dim>    &quadrature,
+         CellData (const hp::FECollection<dim> &fe,
+                   const hp::QCollection<dim>    &quadrature,
                    const Function<dim>      &right_hand_side);
       };
 
       struct FaceData
       {
-         FEFaceValues<dim>    fe_face_values_cell;
-         FEFaceValues<dim>    fe_face_values_neighbor;
-         FESubfaceValues<dim> fe_subface_values_cell;
+         hp::FEFaceValues<dim>    fe_face_values_cell;
+         hp::FEFaceValues<dim>    fe_face_values_neighbor;
+         hp::FESubfaceValues<dim> fe_subface_values_cell;
 
          std::vector<double> jump_residual;
          std::vector<double> dual_weights;       
          typename std::vector<Tensor<1,dim> > cell_grads;
          typename std::vector<Tensor<1,dim> > neighbor_grads;
-         FaceData (const FiniteElement<dim> &fe,
-                   const Quadrature<dim-1>  &face_quadrature);
+         FaceData (const hp::FECollection<dim> &fe,
+                   const hp::QCollection<dim-1>  &face_quadrature);
       };
 
       
@@ -1459,8 +1459,8 @@ namespace LaplaceSolver
 
   template <int dim>
   WeightedResidual<dim>::CellData::
-  CellData (const FiniteElement<dim> &fe,
-           const Quadrature<dim>    &quadrature,
+  CellData (const hp::FECollection<dim> &fe,
+           const hp::QCollection<dim>    &quadrature,
            const Function<dim>      &right_hand_side)
                  :
                  fe_values (fe, quadrature,
@@ -1471,7 +1471,7 @@ namespace LaplaceSolver
                  right_hand_side (&right_hand_side)
   {  
     const unsigned int n_q_points
-      = quadrature.n_quadrature_points;
+      = quadrature[0].n_quadrature_points;
   
     cell_residual.resize(n_q_points);
     rhs_values.resize(n_q_points);    
@@ -1483,8 +1483,8 @@ namespace LaplaceSolver
 
   template <int dim>
   WeightedResidual<dim>::FaceData::
-  FaceData (const FiniteElement<dim> &fe,
-           const Quadrature<dim-1>  &face_quadrature)
+  FaceData (const hp::FECollection<dim> &fe,
+           const hp::QCollection<dim-1>  &face_quadrature)
                  :
                  fe_face_values_cell (fe, face_quadrature,
                                       update_values        |
@@ -1500,7 +1500,7 @@ namespace LaplaceSolver
                                          update_gradients)
   {  
     const unsigned int n_face_q_points
-      = face_quadrature.n_quadrature_points;
+      = face_quadrature[0].n_quadrature_points;
   
     jump_residual.resize(n_face_q_points);
     dual_weights.resize(n_face_q_points);    
@@ -1514,10 +1514,10 @@ namespace LaplaceSolver
   template <int dim>
   WeightedResidual<dim>::
   WeightedResidual (Triangulation<dim>       &coarse_grid,
-                   const FiniteElement<dim> &primal_fe,
-                   const FiniteElement<dim> &dual_fe,
-                   const Quadrature<dim>    &quadrature,
-                   const Quadrature<dim-1>  &face_quadrature,
+                   const hp::FECollection<dim> &primal_fe,
+                   const hp::FECollection<dim> &dual_fe,
+                   const hp::QCollection<dim>    &quadrature,
+                   const hp::QCollection<dim-1>  &face_quadrature,
                    const Function<dim>      &rhs_function,
                    const Function<dim>      &bv,
                    const DualFunctional::DualFunctionalBase<dim> &dual_functional)
@@ -1612,7 +1612,7 @@ namespace LaplaceSolver
                          primal_hanging_node_constraints,
                          dual_solution);    
 
-    DataOut<dim> data_out;
+    DataOut<dim,hp::DoFHandler<dim> > data_out;
     data_out.attach_dof_handler (primal_solver.dof_handler);
 
     data_out.add_data_vector (primal_solver.solution,
@@ -1799,19 +1799,19 @@ namespace LaplaceSolver
   {
     cell_data.fe_values.reinit (cell);
     cell_data.right_hand_side
-      ->value_list (cell_data.fe_values.get_quadrature_points(),
+      ->value_list (cell_data.fe_values.get_present_fe_values().get_quadrature_points(),
                    cell_data.rhs_values);
-    cell_data.fe_values.get_function_2nd_derivatives (primal_solution,
+    cell_data.fe_values.get_present_fe_values().get_function_2nd_derivatives (primal_solution,
                                                      cell_data.cell_grad_grads);
 
-    cell_data.fe_values.get_function_values (dual_weights,
+    cell_data.fe_values.get_present_fe_values().get_function_values (dual_weights,
                                             cell_data.dual_weights);
 
     double sum = 0;
-    for (unsigned int p=0; p<cell_data.fe_values.n_quadrature_points; ++p)
+    for (unsigned int p=0; p<cell_data.fe_values.get_present_fe_values().n_quadrature_points; ++p)
       sum += ((cell_data.rhs_values[p]+trace(cell_data.cell_grad_grads[p])) *
              cell_data.dual_weights[p] *
-             cell_data.fe_values.JxW (p));
+             cell_data.fe_values.get_present_fe_values().JxW (p));
     error_indicators(cell_index) += sum;
   }
 
@@ -1827,10 +1827,10 @@ namespace LaplaceSolver
                               FaceIntegrals              &face_integrals) const
   {
     const unsigned int
-      n_q_points = face_data.fe_face_values_cell.n_quadrature_points;
+      n_q_points = face_data.fe_face_values_cell.get_present_fe_values().n_quadrature_points;
 
     face_data.fe_face_values_cell.reinit (cell, face_no);
-    face_data.fe_face_values_cell.get_function_grads (primal_solution,
+    face_data.fe_face_values_cell.get_present_fe_values().get_function_grads (primal_solution,
                                                      face_data.cell_grads);
 
     Assert (cell->neighbor(face_no).state() == IteratorState::valid,
@@ -1839,22 +1839,22 @@ namespace LaplaceSolver
       neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
     const active_cell_iterator neighbor = cell->neighbor(face_no);
     face_data.fe_face_values_neighbor.reinit (neighbor, neighbor_neighbor);      
-    face_data.fe_face_values_neighbor.get_function_grads (primal_solution,
+    face_data.fe_face_values_neighbor.get_present_fe_values().get_function_grads (primal_solution,
                                                          face_data.neighbor_grads);
 
     for (unsigned int p=0; p<n_q_points; ++p)
       face_data.jump_residual[p]
        = ((face_data.cell_grads[p] - face_data.neighbor_grads[p]) *
-          face_data.fe_face_values_cell.normal_vector(p));
+          face_data.fe_face_values_cell.get_present_fe_values().normal_vector(p));
 
-    face_data.fe_face_values_cell.get_function_values (dual_weights,
+    face_data.fe_face_values_cell.get_present_fe_values().get_function_values (dual_weights,
                                                       face_data.dual_weights);
     
     double face_integral = 0;
     for (unsigned int p=0; p<n_q_points; ++p)
       face_integral += (face_data.jump_residual[p] *
                        face_data.dual_weights[p]  *
-                       face_data.fe_face_values_cell.JxW(p));
+                       face_data.fe_face_values_cell.get_present_fe_values().JxW(p));
 
     Assert (face_integrals.find (cell->face(face_no)) != face_integrals.end(),
            ExcInternalError());
@@ -1876,11 +1876,11 @@ namespace LaplaceSolver
                                 FaceIntegrals              &face_integrals) const
   {
     const unsigned int
-      n_q_points = face_data.fe_face_values_cell.n_quadrature_points;
+      n_q_points = face_data.fe_face_values_cell.get_present_fe_values().n_quadrature_points;
 
-    const typename DoFHandler<dim>::face_iterator
+    const typename hp::DoFHandler<dim>::face_iterator
       face = cell->face(face_no);    
-    const typename DoFHandler<dim>::cell_iterator
+    const typename hp::DoFHandler<dim>::cell_iterator
       neighbor = cell->neighbor(face_no);    
     Assert (neighbor.state() == IteratorState::valid,
            ExcInternalError());
@@ -1900,26 +1900,26 @@ namespace LaplaceSolver
                ExcInternalError());
 
        face_data.fe_subface_values_cell.reinit (cell, face_no, subface_no);
-       face_data.fe_subface_values_cell.get_function_grads (primal_solution,
+       face_data.fe_subface_values_cell.get_present_fe_values().get_function_grads (primal_solution,
                                                             face_data.cell_grads);
        face_data.fe_face_values_neighbor.reinit (neighbor_child,
                                             neighbor_neighbor);
-       face_data.fe_face_values_neighbor.get_function_grads (primal_solution,
+       face_data.fe_face_values_neighbor.get_present_fe_values().get_function_grads (primal_solution,
                                                              face_data.neighbor_grads);
       
        for (unsigned int p=0; p<n_q_points; ++p)
          face_data.jump_residual[p]
             = ((face_data.neighbor_grads[p] - face_data.cell_grads[p]) *
-               face_data.fe_face_values_neighbor.normal_vector(p));
+               face_data.fe_face_values_neighbor.get_present_fe_values().normal_vector(p));
 
-       face_data.fe_face_values_neighbor.get_function_values (dual_weights,
+       face_data.fe_face_values_neighbor.get_present_fe_values().get_function_values (dual_weights,
                                                               face_data.dual_weights);
        
        double face_integral = 0;
        for (unsigned int p=0; p<n_q_points; ++p)
          face_integral += (face_data.jump_residual[p] *
                            face_data.dual_weights[p] *
-                           face_data.fe_face_values_neighbor.JxW(p));
+                           face_data.fe_face_values_neighbor.get_present_fe_values().JxW(p));
        face_integrals[neighbor_child->face(neighbor_neighbor)]
          = face_integral;
       };
@@ -2000,10 +2000,10 @@ void Framework<dim>::run (const ProblemDescription &descriptor)
     triangulation (Triangulation<dim>::smoothing_on_refinement);
   descriptor.data->create_coarse_grid (triangulation);
 
-  const FE_Q<dim>     primal_fe(descriptor.primal_fe_degree);
-  const FE_Q<dim>     dual_fe(descriptor.dual_fe_degree);
-  const QGauss<dim>   quadrature(descriptor.dual_fe_degree+1);
-  const QGauss<dim-1> face_quadrature(descriptor.dual_fe_degree+1);
+  const hp::FECollection<dim>     primal_fe(FE_Q<dim>(descriptor.primal_fe_degree));
+  const hp::FECollection<dim>     dual_fe(FE_Q<dim>(descriptor.dual_fe_degree));
+  const hp::QCollection<dim>   quadrature(QGauss<dim>(descriptor.dual_fe_degree+1));
+  const hp::QCollection<dim-1> face_quadrature(QGauss<dim-1>(descriptor.dual_fe_degree+1));
 
   LaplaceSolver::Base<dim> * solver = 0;
   switch (descriptor.refinement_criterion)
index ec7b7457f51f21fd792f08544338c4e6a2bdd7c9..49c949c68b27329ca4063569049c527cdabcc68e 100644 (file)
@@ -38,7 +38,8 @@ std::ofstream logfile("step-15/output");
 #include <dofs/dof_constraints.h>
 #include <dofs/dof_tools.h>
 #include <fe/fe_q.h>
-#include <fe/fe_values.h>
+#include <fe/hp_fe_values.h>
+#include <fe/q_collection.h>
 #include <numerics/vectors.h>
 #include <numerics/matrices.h>
 #include <numerics/data_out.h>
@@ -101,7 +102,7 @@ class MinimizationProblem
     void output_results () const;
     void refine_grid ();
 
-    static double energy (const DoFHandler<dim> &dof_handler,
+    static double energy (const hp::DoFHandler<dim> &dof_handler,
                           const Vector<double>  &function);
  
 
@@ -109,8 +110,8 @@ class MinimizationProblem
     
     Triangulation<dim>   triangulation;
 
-    FE_Q<dim>            fe;
-    DoFHandler<dim>      dof_handler;
+    hp::FECollection<dim>            fe;
+    hp::DoFHandler<dim>      dof_handler;
 
     ConstraintMatrix     hanging_node_constraints;
     
@@ -127,7 +128,7 @@ template <int dim>
 MinimizationProblem<dim>::MinimizationProblem (const unsigned int run_number)
                 :
                 run_number (run_number),
-                fe (1),
+                fe (FE_Q<dim>(1)),
                dof_handler (triangulation)
 {}
 
@@ -140,7 +141,7 @@ void MinimizationProblem<1>::initialize_solution ()
                             InitializationValues(),
                             present_solution);
 
-  DoFHandler<1>::cell_iterator cell;
+  hp::DoFHandler<1>::cell_iterator cell;
   cell = dof_handler.begin(0);
   while (cell->at_boundary(0) == false)
     cell = cell->neighbor(0);
@@ -183,13 +184,13 @@ void MinimizationProblem<dim>::assemble_step ()
   matrix.reinit (sparsity_pattern);
   residual.reinit (dof_handler.n_dofs());
 
-  QGauss<dim>  quadrature_formula(4);
-  FEValues<dim> fe_values (fe, quadrature_formula, 
+  hp::QCollection<dim>  quadrature_formula(QGauss<dim>(4));
+  hp::FEValues<dim> fe_values (fe, quadrature_formula, 
                           update_values   | update_gradients |
                            update_q_points | update_JxW_values);
 
-  const unsigned int dofs_per_cell = fe.dofs_per_cell;
-  const unsigned int n_q_points    = quadrature_formula.n_quadrature_points;
+  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);
@@ -199,7 +200,7 @@ void MinimizationProblem<dim>::assemble_step ()
   std::vector<double>         local_solution_values (n_q_points);
   std::vector<Tensor<1,dim> > local_solution_grads (n_q_points);
 
-  typename DoFHandler<dim>::active_cell_iterator
+  typename hp::DoFHandler<dim>::active_cell_iterator
     cell = dof_handler.begin_active(),
     endc = dof_handler.end();
   for (; cell!=endc; ++cell)
@@ -209,44 +210,44 @@ void MinimizationProblem<dim>::assemble_step ()
 
       fe_values.reinit (cell);
 
-      fe_values.get_function_values (present_solution,
+      fe_values.get_present_fe_values().get_function_values (present_solution,
                                      local_solution_values);
-      fe_values.get_function_grads (present_solution,
+      fe_values.get_present_fe_values().get_function_grads (present_solution,
                                     local_solution_grads);
 
       for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
         {
           const double u = local_solution_values[q_point],
-                       x = fe_values.quadrature_point(q_point)(0);
+                       x = fe_values.get_present_fe_values().quadrature_point(q_point)(0);
           const double x_minus_u3 = (x-std::pow(u,3));
           const Tensor<1,dim> u_prime = local_solution_grads[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.get_present_fe_values().shape_grad(i,q_point) *
+                    fe_values.get_present_fe_values().shape_grad(j,q_point) *
                     cell->diameter() *
                     cell->diameter()
                     +
-                    fe_values.shape_value(i,q_point) *
-                    fe_values.shape_value(j,q_point)) *
-                fe_values.JxW(q_point);
+                    fe_values.get_present_fe_values().shape_value(i,q_point) *
+                    fe_values.get_present_fe_values().shape_value(j,q_point)) *
+                fe_values.get_present_fe_values().JxW(q_point);
 
           for (unsigned int i=0; i<dofs_per_cell; ++i)
             cell_rhs(i) += -((6. * x_minus_u3 *
                               gradient_power (u_prime, 4) *
-                              fe_values.shape_value(i,q_point)
+                              fe_values.get_present_fe_values().shape_value(i,q_point)
                               *
                               (x_minus_u3 *
                                (u_prime * 
-                                fe_values.shape_grad(i,q_point))
+                                fe_values.get_present_fe_values().shape_grad(i,q_point))
                                -
                                (u_prime*u_prime) * u * u *
-                               fe_values.shape_value(i,q_point))
+                               fe_values.get_present_fe_values().shape_value(i,q_point))
                               )
                              *
-                             fe_values.JxW(q_point));
+                             fe_values.get_present_fe_values().JxW(q_point));
         }
       
       cell->get_dof_indices (local_dof_indices);
@@ -366,7 +367,7 @@ template <int dim>
 void
 MinimizationProblem<dim>::output_results () const
 {
-  DataOut<dim> data_out;
+  DataOut<dim,hp::DoFHandler<dim> > data_out;
   data_out.attach_dof_handler (dof_handler);
   data_out.add_data_vector (present_solution, "solution");
   data_out.build_patches ();
@@ -389,33 +390,34 @@ void MinimizationProblem<1>::refine_grid ()
   
   Vector<float> error_indicators (triangulation.n_active_cells());
 
-  QTrapez<dim> quadrature;
-  FEValues<dim> fe_values (fe, quadrature,
+  QTrapez<dim> q;
+  hp::QCollection<dim> quadrature(q);
+  hp::FEValues<dim> fe_values (fe, quadrature,
                            update_values   | update_gradients |
                            update_second_derivatives |
                            update_q_points | update_JxW_values);
 
-  FEValues<dim> neighbor_fe_values (fe, quadrature,
+  hp::FEValues<dim> neighbor_fe_values (fe, quadrature,
                                     update_gradients);
 
-  std::vector<double> local_values (quadrature.n_quadrature_points);
-  std::vector<Tensor<1,dim> > local_gradients (quadrature.n_quadrature_points);
-  std::vector<Tensor<2,dim> > local_2nd_derivs (quadrature.n_quadrature_points);
+  std::vector<double> local_values (quadrature[0].n_quadrature_points);
+  std::vector<Tensor<1,dim> > local_gradients (quadrature[0].n_quadrature_points);
+  std::vector<Tensor<2,dim> > local_2nd_derivs (quadrature[0].n_quadrature_points);
 
-  DoFHandler<dim>::active_cell_iterator
+  hp::DoFHandler<dim>::active_cell_iterator
     cell = dof_handler.begin_active (),
     endc = dof_handler.end ();
   for (unsigned int cell_index = 0; cell!=endc; ++cell, ++cell_index)
     {
       fe_values.reinit (cell);
-      fe_values.get_function_values (present_solution, local_values);
-      fe_values.get_function_grads (present_solution, local_gradients);
-      fe_values.get_function_2nd_derivatives (present_solution, local_2nd_derivs);
+      fe_values.get_present_fe_values().get_function_values (present_solution, local_values);
+      fe_values.get_present_fe_values().get_function_grads (present_solution, local_gradients);
+      fe_values.get_present_fe_values().get_function_2nd_derivatives (present_solution, local_2nd_derivs);
 
       double cell_residual_norm = 0;
-      for (unsigned int q=0; q<quadrature.n_quadrature_points; ++q)
+      for (unsigned int q=0; q<quadrature[0].n_quadrature_points; ++q)
         {
-          const double x             = fe_values.quadrature_point(q)[0];
+          const double x             = fe_values.get_present_fe_values().quadrature_point(q)[0];
           const double u             = local_values[q];
           const double u_prime       = local_gradients[q][0];
           const double u_doubleprime = local_2nd_derivs[q][0][0];
@@ -428,13 +430,13 @@ void MinimizationProblem<1>::refine_grid ()
                 2*u_prime*(1-3*u*u*u_prime)));
           
           cell_residual_norm += (local_residual_value * local_residual_value *
-                                 fe_values.JxW(q));
+                                 fe_values.get_present_fe_values().JxW(q));
         }
       error_indicators(cell_index) = cell_residual_norm *
                                      cell->diameter() * cell->diameter();
 
-      const double x_left  = fe_values.quadrature_point(0)[0];
-      const double x_right = fe_values.quadrature_point(1)[0];
+      const double x_left  = fe_values.get_present_fe_values().quadrature_point(0)[0];
+      const double x_right = fe_values.get_present_fe_values().quadrature_point(1)[0];
 
       Assert (x_left  == cell->vertex(0)[0], ExcInternalError());
       Assert (x_right == cell->vertex(1)[0], ExcInternalError());
@@ -447,12 +449,12 @@ void MinimizationProblem<1>::refine_grid ()
 
       if (cell->at_boundary(0) == false)
         {
-          DoFHandler<dim>::cell_iterator left_neighbor = cell->neighbor(0);
+          hp::DoFHandler<dim>::cell_iterator left_neighbor = cell->neighbor(0);
           while (left_neighbor->has_children())
             left_neighbor = left_neighbor->child(1);
 
           neighbor_fe_values.reinit (left_neighbor);
-          neighbor_fe_values.get_function_grads (present_solution, local_gradients);
+          neighbor_fe_values.get_present_fe_values().get_function_grads (present_solution, local_gradients);
 
           const double neighbor_u_prime_left = local_gradients[1][0];
 
@@ -465,12 +467,12 @@ void MinimizationProblem<1>::refine_grid ()
 
       if (cell->at_boundary(1) == false)
         {
-          DoFHandler<dim>::cell_iterator right_neighbor = cell->neighbor(1);
+          hp::DoFHandler<dim>::cell_iterator right_neighbor = cell->neighbor(1);
           while (right_neighbor->has_children())
             right_neighbor = right_neighbor->child(0);
           
           neighbor_fe_values.reinit (right_neighbor);
-          neighbor_fe_values.get_function_grads (present_solution, local_gradients);
+          neighbor_fe_values.get_present_fe_values().get_function_grads (present_solution, local_gradients);
 
           const double neighbor_u_prime_right = local_gradients[0][0];
 
@@ -508,39 +510,39 @@ void MinimizationProblem<1>::refine_grid ()
 
 template <int dim>
 double
-MinimizationProblem<dim>::energy (const DoFHandler<dim> &dof_handler,
+MinimizationProblem<dim>::energy (const hp::DoFHandler<dim> &dof_handler,
                                   const Vector<double>  &function)
 {
-  QGauss<dim>  quadrature_formula(4);
-  FEValues<dim> fe_values (dof_handler.get_fe(), quadrature_formula, 
+  hp::QCollection<dim>  quadrature_formula(QGauss<dim>(4));
+  hp::FEValues<dim> fe_values (dof_handler.get_fe(), quadrature_formula, 
                           update_values   | update_gradients |
                            update_q_points | update_JxW_values);
 
-  const unsigned int   n_q_points    = quadrature_formula.n_quadrature_points;
+  const unsigned int   n_q_points    = quadrature_formula[0].n_quadrature_points;
 
   std::vector<double>         local_solution_values (n_q_points);
   std::vector<Tensor<1,dim> > local_solution_grads (n_q_points);
 
   double energy = 0.;
 
-  typename DoFHandler<dim>::active_cell_iterator
+  typename hp::DoFHandler<dim>::active_cell_iterator
     cell = dof_handler.begin_active(),
     endc = dof_handler.end();
   for (; cell!=endc; ++cell)
     {
       fe_values.reinit (cell);
-      fe_values.get_function_values (function,
+      fe_values.get_present_fe_values().get_function_values (function,
                                      local_solution_values);
-      fe_values.get_function_grads (function,
+      fe_values.get_present_fe_values().get_function_grads (function,
                                     local_solution_grads);
 
       for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
-        energy += (std::pow (fe_values.quadrature_point(q_point)(0)
+        energy += (std::pow (fe_values.get_present_fe_values().quadrature_point(q_point)(0)
                              -
                              std::pow (local_solution_values[q_point], 3),
                              2) *
                    gradient_power (local_solution_grads[q_point], 6) *
-                   fe_values.JxW (q_point));
+                   fe_values.get_present_fe_values().JxW (q_point));
     }
 
   return energy;

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.