]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Merge up from the upstream branch.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 23 Aug 2012 03:58:24 +0000 (03:58 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 23 Aug 2012 03:58:24 +0000 (03:58 +0000)
git-svn-id: https://svn.dealii.org/branches/branch_bigger_global_dof_indices_3@26096 0785d39b-7218-0410-832d-ea1e28bc413d

25 files changed:
1  2 
deal.II/contrib/boost-1.49.0/include/boost/fusion/tuple/detail/preprocessed/tuple.hpp
deal.II/examples/step-42/step-42.cc
deal.II/include/deal.II/base/types.h
deal.II/include/deal.II/dofs/dof_handler.h
deal.II/include/deal.II/dofs/dof_renumbering.h
deal.II/include/deal.II/dofs/dof_tools.h
deal.II/include/deal.II/hp/dof_handler.h
deal.II/include/deal.II/lac/block_vector.h
deal.II/include/deal.II/lac/block_vector_base.h
deal.II/include/deal.II/multigrid/mg_transfer.h
deal.II/include/deal.II/multigrid/mg_transfer.templates.h
deal.II/include/deal.II/numerics/vector_tools.h
deal.II/include/deal.II/numerics/vector_tools.templates.h
deal.II/source/distributed/tria.cc
deal.II/source/dofs/dof_handler.cc
deal.II/source/dofs/dof_handler_policy.cc
deal.II/source/dofs/dof_renumbering.cc
deal.II/source/dofs/dof_renumbering.inst.in
deal.II/source/dofs/dof_tools.cc
deal.II/source/dofs/dof_tools.inst.in
deal.II/source/fe/fe_tools.cc
deal.II/source/hp/dof_handler.cc
deal.II/source/multigrid/mg_tools.cc
deal.II/source/multigrid/mg_tools.inst.in
tests/fe/abf_01.cc

index 8f9b94a5cd83c8c0fa1f97d5ecd846f04fb53088,8f9b94a5cd83c8c0fa1f97d5ecd846f04fb53088..e1cd851940ba3919ec0bdab15a5bda293392e333
@@@ -18,4 -18,4 +18,4 @@@
  #include <boost/fusion/tuple/detail/preprocessed/tuple50.hpp>
  #else
  #error "FUSION_MAX_VECTOR_SIZE out of bounds for preprocessed headers"
--#endif
++#endif
index 8fbb364edc9fed1f12585bf17590975885f325b7,3a9a898ef80aac1250d5466b0b3e2649ff1c51a4..54bdc047368196cf3307a5da0ecc72088b2c9176
@@@ -433,1053 -476,921 +476,920 @@@ namespace Step4
    }
  
    template <int dim>
-   void Obstacle<dim>::vector_value (const Point<dim> &p,
-                                   Vector<double> &values) const
+   void PlasticityContactProblem<dim>::make_grid ()
    {
-     for (unsigned int c=0; c<this->n_components; ++c)
-       values(c) = Obstacle<dim>::value (p, c);
+     std::vector<unsigned int> repet(3);
+     repet[0] = 1;//20;
+     repet[1] = 1;
+     repet[2] = 1;
+     Point<dim> p1 (0,0,0);
+     Point<dim> p2 (1.0, 1.0, 1.0);
+     GridGenerator::subdivided_hyper_rectangle (triangulation, repet, p1, p2);
+     Triangulation<3>::active_cell_iterator
+       cell = triangulation.begin_active(),
+       endc = triangulation.end();
+     /* boundary_indicators:
+               _______
+              /  9    /|
+             /______ / |
+           8|       | 8|
+            |   8   | /
+            |_______|/
+                6
+      */
+     for (; cell!=endc; ++cell)
+       for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+         {
+           if (cell->face (face)->center ()[2] == p2(2))
+             cell->face (face)->set_boundary_indicator (9);
+           if (cell->face (face)->center ()[0] == p1(0) ||
+               cell->face (face)->center ()[0] == p2(0) ||
+               cell->face (face)->center ()[1] == p1(1) ||
+               cell->face (face)->center ()[1] == p2(1))
+             cell->face (face)->set_boundary_indicator (8);
+           if (cell->face (face)->center ()[2] == p1(2))
+             cell->face (face)->set_boundary_indicator (6);
+         }
+     triangulation.refine_global (n_refinements_global);
+     // Lokale Verfeinerung des Gitters
+     for (int step=0; step<n_refinements_local; ++step)
+       {
+         cell = triangulation.begin_active();  // Iterator ueber alle Zellen
+         for (; cell!=endc; ++cell)
+            for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+            {
+   //     if (cell->face (face)->at_boundary()
+   //         && cell->face (face)->boundary_indicator () == 9)
+   //       {
+   //         cell->set_refine_flag ();
+   //         break;
+   //       }
+   //     else if (cell->level () == n_refinements + n_refinements_local - 1)
+   //       {
+   //         cell->set_refine_flag ();
+   //         break;
+   //       }
+   //     if (cell->face (face)->at_boundary()
+   //         && cell->face (face)->boundary_indicator () == 9)
+   //       {
+   //         if (cell->face (face)->vertex (0)(0) <= 0.7 &&
+   //             cell->face (face)->vertex (1)(0) >= 0.3 &&
+   //             cell->face (face)->vertex (0)(1) <= 0.875 &&
+   //             cell->face (face)->vertex (2)(1) >= 0.125)
+   //                   {
+   //             cell->set_refine_flag ();
+   //             break;
+   //           }
+   //       }
+               if (step == 0 &&
+                   cell->center ()(2) < n_refinements_local*9.0/64.0)
+                {
+                  cell->set_refine_flag ();
+                  break;
+                }
+           };
+         triangulation.execute_coarsening_and_refinement ();
+       };
    }
- }
  
+   template <int dim>
+   void PlasticityContactProblem<dim>::setup_system ()
+   {
+     // setup dofs
+     {
+       dof_handler.distribute_dofs (fe);
  
-                                  // @sect3{Implementation of the <code>Step4</code> class}
-                                  // Next for the implementation of the class
-                                  // template that makes use of the functions
-                                  // above. As before, we will write everything
- template <int dim>
- Step4<dim>::Step4 (int _n_refinements_global, int _n_refinements_local)
-   :
-   n_refinements_global (_n_refinements_global),
-   n_refinements_local (_n_refinements_local),
-   mpi_communicator (MPI_COMM_WORLD),
-   triangulation (mpi_communicator),
-   fe (FE_Q<dim>(1), dim),
-   dof_handler (triangulation),
-   pcout (std::cout,
-        (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)),
-   sigma_0 (400),
-   gamma (1.e-2),
-   e_modul (2.0e5),
-   nu (0.3)
- {
-   // double _E, double _nu, double _sigma_0, double _gamma
-   plast_lin_hard.reset (new ConstitutiveLaw<dim> (e_modul, nu, sigma_0, gamma, mpi_communicator, pcout));
- }
+       locally_owned_dofs = dof_handler.locally_owned_dofs ();
+       locally_relevant_dofs.clear();
+       DoFTools::extract_locally_relevant_dofs (dof_handler,
+                                                locally_relevant_dofs);
+     }
  
- template <int dim>
- void Step4<dim>::make_grid ()
- {
-   std::vector<unsigned int> repet(3);
-   repet[0] = 1;//20;
-   repet[1] = 1;
-   repet[2] = 1;
-   Point<dim> p1 (0,0,0);
-   Point<dim> p2 (1.0, 1.0, 1.0);
-   GridGenerator::subdivided_hyper_rectangle (triangulation, repet, p1, p2);
-   Triangulation<3>::active_cell_iterator
-     cell = triangulation.begin_active(),
-     endc = triangulation.end();
-   /* boundary_indicators:
-             _______
-            /  9    /|
-           /______ / |
-         8|       | 8|
-          |   8   | /
-          |_______|/
-              6
-    */
-   for (; cell!=endc; ++cell)
-     for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-       {
-       if (cell->face (face)->center ()[2] == p2(2))
-         cell->face (face)->set_boundary_indicator (9);
-       if (cell->face (face)->center ()[0] == p1(0) ||
-           cell->face (face)->center ()[0] == p2(0) ||
-           cell->face (face)->center ()[1] == p1(1) ||
-           cell->face (face)->center ()[1] == p2(1))
-         cell->face (face)->set_boundary_indicator (8);
-       if (cell->face (face)->center ()[2] == p1(2))
-         cell->face (face)->set_boundary_indicator (6);
-       }
+     // setup hanging nodes and dirichlet constraints
+     {
+       constraints_hanging_nodes.reinit (locally_relevant_dofs);
+       DoFTools::make_hanging_node_constraints (dof_handler,
+                                                constraints_hanging_nodes);
+       constraints_hanging_nodes.close ();
+       pcout << "Number of active cells: "
+                  << triangulation.n_active_cells()
+                  << std::endl
+                  << "Total number of cells: "
+                  << triangulation.n_cells()
+                  << std::endl
+                  << "Number of degrees of freedom: "
+                  << dof_handler.n_dofs ()
+                  << std::endl;
+       dirichlet_constraints ();
+     }
  
-   triangulation.refine_global (n_refinements_global);
+     // Initialzation for matrices and vectors
+     {
+       solution.reinit (locally_relevant_dofs, mpi_communicator);
+       system_rhs_newton.reinit (locally_owned_dofs, mpi_communicator);
+       old_solution.reinit (system_rhs_newton);
+       resid_vector.reinit (system_rhs_newton);
+       diag_mass_matrix_vector.reinit (system_rhs_newton);
+       active_set.set_size (locally_relevant_dofs.size ());
+     }
  
-   // Lokale Verfeinerung des Gitters
-   for (int step=0; step<n_refinements_local; ++step)
+     // setup sparsity pattern
      {
-       cell = triangulation.begin_active();  // Iterator ueber alle Zellen
-       for (; cell!=endc; ++cell)
-          for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-          {
- //       if (cell->face (face)->at_boundary()
- //           && cell->face (face)->boundary_indicator () == 9)
- //         {
- //           cell->set_refine_flag ();
- //           break;
- //         }
- //       else if (cell->level () == n_refinements + n_refinements_local - 1)
- //         {
- //           cell->set_refine_flag ();
- //           break;
- //         }
- //       if (cell->face (face)->at_boundary()
- //           && cell->face (face)->boundary_indicator () == 9)
- //         {
- //           if (cell->face (face)->vertex (0)(0) <= 0.7 &&
- //               cell->face (face)->vertex (1)(0) >= 0.3 &&
- //               cell->face (face)->vertex (0)(1) <= 0.875 &&
- //               cell->face (face)->vertex (2)(1) >= 0.125)
- //             {
- //               cell->set_refine_flag ();
- //               break;
- //             }
- //         }
-           
-           if (step == 0 &&
-               cell->center ()(2) < n_refinements_local*9.0/64.0)
-            {
-              cell->set_refine_flag ();
-              break;
-            }     
-         };
-       triangulation.execute_coarsening_and_refinement ();
-     };
- }
+       TrilinosWrappers::SparsityPattern sp (locally_owned_dofs,
+                                             mpi_communicator);
  
- template <int dim>
- void Step4<dim>::setup_system ()
- {
-   // setup dofs
-   {
-     dof_handler.distribute_dofs (fe);
+       DoFTools::make_sparsity_pattern (dof_handler, sp, constraints_dirichlet_hanging_nodes, false,
+                                        Utilities::MPI::this_mpi_process(mpi_communicator));
+       sp.compress();
+       system_matrix_newton.reinit (sp);
  
-     locally_owned_dofs = dof_handler.locally_owned_dofs ();
-     locally_relevant_dofs.clear();
-     DoFTools::extract_locally_relevant_dofs (dof_handler,
-                                            locally_relevant_dofs);
+       TrilinosWrappers::SparseMatrix mass_matrix;
+       mass_matrix.reinit (sp);
+       assemble_mass_matrix_diagonal (mass_matrix);
+       const unsigned int
+         start = (system_rhs_newton.local_range().first),
+         end   = (system_rhs_newton.local_range().second);
+       for (unsigned int j=start; j<end; j++)
+         diag_mass_matrix_vector (j) = mass_matrix.diag_element (j);
+       number_iterations = 0;
+       diag_mass_matrix_vector.compress ();
+     }
    }
  
-   // setup hanging nodes and dirichlet constraints
+   template <int dim>
+   void PlasticityContactProblem<dim>::assemble_nl_system (TrilinosWrappers::MPI::Vector &u)
    {
-     // constraints_hanging_nodes.clear ();
-     constraints_hanging_nodes.reinit (locally_relevant_dofs);
-     DoFTools::make_hanging_node_constraints (dof_handler,
-                                            constraints_hanging_nodes);
-     constraints_hanging_nodes.close ();
-     pcout << "Number of active cells: "
-              << triangulation.n_active_cells()
-              << std::endl
-              << "Total number of cells: "
-              << triangulation.n_cells()
-              << std::endl
-              << "Number of degrees of freedom: "
-              << dof_handler.n_dofs ()
-              << std::endl;
-     dirichlet_constraints ();
+     QGauss<dim>  quadrature_formula(2);
+     QGauss<dim-1>  face_quadrature_formula(2);
+     FEValues<dim> fe_values (fe, quadrature_formula,
+                              UpdateFlags(update_values    |
+                                          update_gradients |
+                                          update_q_points  |
+                                          update_JxW_values));
+     FEFaceValues<dim> fe_values_face (fe, face_quadrature_formula,
+                                       update_values   | update_quadrature_points |
+                                       update_JxW_values);
+     const unsigned int   dofs_per_cell = fe.dofs_per_cell;
+     const unsigned int   n_q_points    = quadrature_formula.size ();
+     const unsigned int   n_face_q_points = face_quadrature_formula.size();
+     const EquationData::RightHandSide<dim> right_hand_side;
+     std::vector<Vector<double> > right_hand_side_values (n_q_points,
+                                                          Vector<double>(dim));
+     std::vector<Vector<double> > right_hand_side_values_face (n_face_q_points,
+                                                               Vector<double>(dim));
+     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<types::global_dof_index> local_dof_indices (dofs_per_cell);
+     typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
+                                                    endc = dof_handler.end();
+     const FEValuesExtractors::Vector displacement (0);
+     TrilinosWrappers::MPI::Vector   test_rhs(solution);
+     const double kappa = 1.0;
+     for (; cell!=endc; ++cell)
+       if (cell->is_locally_owned())
+         {
+           fe_values.reinit (cell);
+           cell_matrix = 0;
+           cell_rhs = 0;
+           right_hand_side.vector_value_list (fe_values.get_quadrature_points(),
+                                              right_hand_side_values);
+           std::vector<SymmetricTensor<2,dim> > strain_tensor (n_q_points);
+           fe_values[displacement].get_function_symmetric_gradients (u, strain_tensor);
+           for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+             {
+               SymmetricTensor<4,dim> stress_strain_tensor_linearized;
+               SymmetricTensor<4,dim> stress_strain_tensor;
+               SymmetricTensor<2,dim> stress_tensor;
+               plast_lin_hard->linearized_plast_linear_hardening (stress_strain_tensor_linearized,
+                                                                  stress_strain_tensor,
+                                                                  strain_tensor[q_point]);
+               //      if (q_point == 0)
+               //      std::cout<< stress_strain_tensor_linearized <<std::endl;
+               //      std::cout<< stress_strain_tensor <<std::endl;
+               for (unsigned int i=0; i<dofs_per_cell; ++i)
+                 {
+                   stress_tensor = stress_strain_tensor_linearized * plast_lin_hard->get_strain(fe_values, i, q_point);
+                   for (unsigned int j=0; j<dofs_per_cell; ++j)
+                     {
+                       cell_matrix(i,j) += (stress_tensor *
+                                            plast_lin_hard->get_strain(fe_values, j, q_point) *
+                                            fe_values.JxW (q_point));
+                     }
+                   // the linearized part a(v^i;v^i,v) of the rhs
+                   cell_rhs(i) += (stress_tensor *
+                                   strain_tensor[q_point] *
+                                   fe_values.JxW (q_point));
+                   // the residual part a(v^i;v) of the rhs
+                   cell_rhs(i) -= (strain_tensor[q_point] * stress_strain_tensor *
+                                   plast_lin_hard->get_strain(fe_values, i, q_point) *
+                                   fe_values.JxW (q_point));
+                   // the residual part F(v) of the rhs
+                   Tensor<1,dim> rhs_values;
+                   rhs_values = 0;
+                   cell_rhs(i) += (fe_values[displacement].value (i, q_point) *
+                                   rhs_values *
+                                   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 () == 9)
+                 {
+                   fe_values_face.reinit (cell, face);
+                   right_hand_side.vector_value_list (fe_values_face.get_quadrature_points(),
+                                                      right_hand_side_values_face);
+                   for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
+                     {
+                       Tensor<1,dim> rhs_values;
+                       rhs_values = 0;
+                       for (unsigned int i=0; i<dofs_per_cell; ++i)
+                         cell_rhs(i) += (fe_values_face[displacement].value (i, q_point) *
+                                         rhs_values *
+                                         fe_values_face.JxW (q_point));
+                     }
+                 }
+             }
+           cell->get_dof_indices (local_dof_indices);
+           constraints.distribute_local_to_global (cell_matrix, cell_rhs,
+                                                   local_dof_indices,
+                                                   system_matrix_newton, system_rhs_newton, true);
+         };
+     system_matrix_newton.compress ();
+     system_rhs_newton.compress (Add);
    }
  
-   // Initialzation for matrices and vectors
+   template <int dim>
+   void PlasticityContactProblem<dim>::residual_nl_system (TrilinosWrappers::MPI::Vector &u)
    {
-     solution.reinit (locally_relevant_dofs, mpi_communicator);
-     system_rhs_newton.reinit (locally_owned_dofs, mpi_communicator);
-     old_solution.reinit (system_rhs_newton);
-     resid_vector.reinit (system_rhs_newton);
-     diag_mass_matrix_vector.reinit (system_rhs_newton);
-     active_set.set_size (locally_relevant_dofs.size ());
+     QGauss<dim>  quadrature_formula(2);
+     QGauss<dim-1> face_quadrature_formula(2);
+     FEValues<dim> fe_values (fe, quadrature_formula,
+                              UpdateFlags(update_values    |
+                                          update_gradients |
+                                          update_q_points  |
+                                          update_JxW_values));
+     FEFaceValues<dim> fe_values_face (fe, face_quadrature_formula,
+                                       update_values   | update_quadrature_points |
+                                       update_JxW_values);
+     const unsigned int   dofs_per_cell = fe.dofs_per_cell;
+     const unsigned int   n_q_points    = quadrature_formula.size ();
+     const unsigned int   n_face_q_points = face_quadrature_formula.size();
+     const EquationData::RightHandSide<dim> right_hand_side;
+     std::vector<Vector<double> > right_hand_side_values (n_q_points,
+                                                          Vector<double>(dim));
+     std::vector<Vector<double> > right_hand_side_values_face (n_face_q_points,
+                                                               Vector<double>(dim));
+     Vector<double>       cell_rhs (dofs_per_cell);
+     std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+     const FEValuesExtractors::Vector displacement (0);
+     typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
+                                                    endc = dof_handler.end();
+     unsigned int elast_points = 0;
+     unsigned int plast_points = 0;
+     double       yield = 0;
+     unsigned int cell_number = 0;
+     for (; cell!=endc; ++cell)
+       if (cell->is_locally_owned())
+         {
+           fe_values.reinit (cell);
+           cell_rhs = 0;
+           right_hand_side.vector_value_list (fe_values.get_quadrature_points(),
+                                              right_hand_side_values);
+           std::vector<SymmetricTensor<2,dim> > strain_tensor (n_q_points);
+           fe_values[displacement].get_function_symmetric_gradients (u, strain_tensor);
+           for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+             {
+               SymmetricTensor<4,dim> stress_strain_tensor;
+               SymmetricTensor<2,dim> stress_tensor;
+               plast_lin_hard->plast_linear_hardening (stress_strain_tensor, strain_tensor[q_point],
+                                                       elast_points, plast_points, yield);
+               for (unsigned int i=0; i<dofs_per_cell; ++i)
+                 {
+                   cell_rhs(i) -= (strain_tensor[q_point] * stress_strain_tensor * //(stress_tensor) *
+                                   plast_lin_hard->get_strain(fe_values, i, q_point) *
+                                   fe_values.JxW (q_point));
+                   Tensor<1,dim> rhs_values;
+                   rhs_values = 0;
+                   cell_rhs(i) += ((fe_values[displacement].value (i, q_point) *
+                                    rhs_values) *
+                                   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 () == 9)
+                 {
+                   fe_values_face.reinit (cell, face);
+                   right_hand_side.vector_value_list (fe_values_face.get_quadrature_points(),
+                                                      right_hand_side_values_face);
+                   for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
+                     {
+                       Tensor<1,dim> rhs_values;
+                       rhs_values = 0;
+                       for (unsigned int i=0; i<dofs_per_cell; ++i)
+                         cell_rhs(i) += (fe_values_face[displacement].value (i, q_point) *
+                                         rhs_values *
+                                         fe_values_face.JxW (q_point));
+                     }
+                 }
+             }
+           cell->get_dof_indices (local_dof_indices);
+           constraints_dirichlet_hanging_nodes.distribute_local_to_global (cell_rhs,
+                                                                           local_dof_indices,
+                                                                           system_rhs_newton);
+           cell_number += 1;
+         };
+     system_rhs_newton.compress ();
+     unsigned int sum_elast_points = Utilities::MPI::sum(elast_points, mpi_communicator);
+     unsigned int sum_plast_points = Utilities::MPI::sum(plast_points, mpi_communicator);
+     pcout<< "Elast-Points = " << sum_elast_points <<std::endl;
+     pcout<< "Plast-Points = " << sum_plast_points <<std::endl;
    }
  
-   // setup sparsity pattern
+   template <int dim>
+   void PlasticityContactProblem<dim>::assemble_mass_matrix_diagonal (TrilinosWrappers::SparseMatrix &mass_matrix)
    {
-     TrilinosWrappers::SparsityPattern sp (locally_owned_dofs,
-                                         mpi_communicator);
+     QTrapez<dim-1>  face_quadrature_formula;
  
-     DoFTools::make_sparsity_pattern (dof_handler, sp, constraints_dirichlet_hanging_nodes, false,
-                                    Utilities::MPI::this_mpi_process(mpi_communicator));
+     FEFaceValues<dim> fe_values_face (fe, face_quadrature_formula,
+                                       update_values   |
+                                       update_quadrature_points |
+                                       update_JxW_values);
  
-     sp.compress();
+     const unsigned int   dofs_per_cell      = fe.dofs_per_cell;
+     const unsigned int   n_face_q_points    = face_quadrature_formula.size();
  
-     system_matrix_newton.reinit (sp);
+     FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
  
-     mass_matrix.reinit (sp);
-   }    
+     std::vector<unsigned int> local_dof_indices (dofs_per_cell);
  
-   assemble_mass_matrix ();
-   const unsigned int
-     start = (system_rhs_newton.local_range().first),
-     end   = (system_rhs_newton.local_range().second);
-   for (unsigned int j=start; j<end; j++)
-     diag_mass_matrix_vector (j) = mass_matrix.diag_element (j);
-   number_iterations = 0;
-   
-   diag_mass_matrix_vector.compress ();
- }
+     const FEValuesExtractors::Vector displacement (0);
  
- template <int dim>
- void Step4<dim>::assemble_mass_matrix ()
- {
-   QTrapez<dim-1>  face_quadrature_formula;
+     typename DoFHandler<dim>::active_cell_iterator
+       cell = dof_handler.begin_active(),
+       endc = dof_handler.end();
  
-   FEFaceValues<dim> fe_values_face (fe, face_quadrature_formula,
-                                     update_values   | update_quadrature_points | update_JxW_values);
+     for (; cell!=endc; ++cell)
+       if (cell->is_locally_owned())
+         for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+           if (cell->face (face)->at_boundary()
+               && cell->face (face)->boundary_indicator () == 9)
+             {
+               fe_values_face.reinit (cell, face);
+               cell_matrix = 0;
  
-   const unsigned int   dofs_per_cell      = fe.dofs_per_cell;
-   const unsigned int   dofs_per_face      = fe.dofs_per_face;
-   const unsigned int   n_face_q_points    = face_quadrature_formula.size();
+               for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
+                 for (unsigned int i=0; i<dofs_per_cell; ++i)
+                   cell_matrix(i,i) += (fe_values_face[displacement].value (i, q_point) *
+                                        fe_values_face[displacement].value (i, q_point) *
+                                        fe_values_face.JxW (q_point));
  
-   FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+               cell->get_dof_indices (local_dof_indices);
  
-   std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+               constraints_dirichlet_hanging_nodes.distribute_local_to_global (cell_matrix,
+                                                                               local_dof_indices,
+                                                                               mass_matrix);
+             }
  
-   const FEValuesExtractors::Vector displacement (0);
+     mass_matrix.compress ();
+   }
  
-   typename DoFHandler<dim>::active_cell_iterator
-     cell = dof_handler.begin_active(),
-     endc = dof_handler.end();
+                                    // @sect4{PlasticityContactProblem::update_solution_and_constraints}
  
-   for (; cell!=endc; ++cell)
-     if (cell->is_locally_owned())
-       for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-       if (cell->face (face)->at_boundary()
-           && cell->face (face)->boundary_indicator () == 9)
-         {
-           fe_values_face.reinit (cell, face);
-           cell_matrix = 0;
-           for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
-             for (unsigned int i=0; i<dofs_per_cell; ++i)
-               cell_matrix(i,i) += (fe_values_face[displacement].value (i, q_point) *
-                                    fe_values_face[displacement].value (i, q_point) *
-                                    fe_values_face.JxW (q_point));
-           cell->get_dof_indices (local_dof_indices);
-           constraints_dirichlet_hanging_nodes.distribute_local_to_global (cell_matrix,
-                                                                           local_dof_indices,
-                                                                           mass_matrix);
-         }
-   mass_matrix.compress ();
- }
+                                    // Projection and updating of the active set
+                                    // for the dofs which penetrates the obstacle.
+   template <int dim>
+   void PlasticityContactProblem<dim>::update_solution_and_constraints ()
+   {
+     clock_t                        start_proj, end_proj;
  
- template <int dim>
- void Step4<dim>::assemble_nl_system (TrilinosWrappers::MPI::Vector &u)
- {
-   QGauss<dim>  quadrature_formula(2);
-   QGauss<dim-1>  face_quadrature_formula(2);
+     const EquationData::Obstacle<dim>     obstacle;
+     std::vector<bool>                     vertex_touched (dof_handler.n_dofs (), false);
  
-   FEValues<dim> fe_values (fe, quadrature_formula,
-                          UpdateFlags(update_values    |
-                                      update_gradients |
-                                      update_q_points  |
-                                      update_JxW_values));
+     typename DoFHandler<dim>::active_cell_iterator
+     cell = dof_handler.begin_active(),
+     endc = dof_handler.end();
  
-   FEFaceValues<dim> fe_values_face (fe, face_quadrature_formula,
-                                   update_values   | update_quadrature_points |
-                                   update_JxW_values);
+     TrilinosWrappers::MPI::Vector         distributed_solution (system_rhs_newton);
+     distributed_solution = solution;
+     TrilinosWrappers::MPI::Vector         lambda (solution);
+     lambda = resid_vector;
+     TrilinosWrappers::MPI::Vector         diag_mass_matrix_vector_relevant (solution);
+     diag_mass_matrix_vector_relevant = diag_mass_matrix_vector;
+     constraints.reinit(locally_relevant_dofs);
+     active_set.clear ();
+     IndexSet     active_set_locally_owned;
+     active_set_locally_owned.set_size (locally_owned_dofs.size ());
+     const double c = 100.0*e_modul;
+     for (; cell!=endc; ++cell)
+       if (cell->is_locally_owned())
+         for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+           if (cell->face (face)->at_boundary()
+               && cell->face (face)->boundary_indicator () == 9)
+             for (unsigned int v=0; v<GeometryInfo<dim-1>::vertices_per_cell; ++v)
+               {
+                 unsigned int index_z = cell->face (face)->vertex_dof_index (v,2);
+                 if (vertex_touched[cell->face (face)->vertex_index(v)] == false)
+                   vertex_touched[cell->face (face)->vertex_index(v)] = true;
+                 else
+                   continue;
+                 // the local row where
+                 Point<dim> point (cell->face (face)->vertex (v)[0],/* + solution (index_x),*/
+                                   cell->face (face)->vertex (v)[1],
+                                   cell->face (face)->vertex (v)[2]);
+                 double obstacle_value = obstacle.value (point, 2);
+                 double solution_index_z = solution (index_z);
+                 double gap = obstacle_value - point (2);
+                 if (lambda (index_z) +
+                     c *
+                     diag_mass_matrix_vector_relevant (index_z) *
+                     (solution_index_z - gap)
+                     > 0)
+                   {
+                     constraints.add_line (index_z);
+                     constraints.set_inhomogeneity (index_z, gap);
+                     distributed_solution (index_z) = gap;
+                     if (locally_relevant_dofs.is_element (index_z))
+                       active_set.add_index (index_z);
+                     if (locally_owned_dofs.is_element (index_z))
+                       active_set_locally_owned.add_index (index_z);
+                     // std::cout<< index_z << ", "
+                     //           << "Error: " << lambda (index_z) +
+                     //   diag_mass_matrix_vector_relevant (index_z)*c*(solution_index_z - gap)
+                     //           << ", " << lambda (index_z)
+                     //           << ", " << diag_mass_matrix_vector_relevant (index_z)
+                     //           << ", " << obstacle_value
+                     //           << ", " << solution_index_z
+                     //           <<std::endl;
+                   }
+               }
+     distributed_solution.compress(Insert);
+     unsigned int sum_contact_constraints = Utilities::MPI::sum(active_set_locally_owned.n_elements (),
+                                                                mpi_communicator);
+     pcout << "Number of Contact-Constaints: " << sum_contact_constraints <<std::endl;
+     solution = distributed_solution;
+     constraints.close ();
+     const ConstraintMatrix::MergeConflictBehavior
+       merge_conflict_behavior = ConstraintMatrix::left_object_wins;
+     constraints.merge (constraints_dirichlet_hanging_nodes, merge_conflict_behavior);
+   }
  
-   const unsigned int   dofs_per_cell = fe.dofs_per_cell;
-   const unsigned int   n_q_points    = quadrature_formula.size ();
-   const unsigned int   n_face_q_points = face_quadrature_formula.size();
+   template <int dim>
+   void PlasticityContactProblem<dim>::dirichlet_constraints ()
+   {
+     /* boundary_indicators:
+               _______
+              /  9    /|
+             /______ / |
+           8|       | 8|
+            |   8   | /
+            |_______|/
+                6
+      */
+     constraints_dirichlet_hanging_nodes.reinit (locally_relevant_dofs);
+     constraints_dirichlet_hanging_nodes.merge (constraints_hanging_nodes);
+     std::vector<bool> component_mask (dim, true);
+     component_mask[0] = true;
+     component_mask[1] = true;
+     component_mask[2] = true;
+     VectorTools::interpolate_boundary_values (dof_handler,
+                                               6,
+                                               EquationData::BoundaryValues<dim>(),
+                                               constraints_dirichlet_hanging_nodes,
+                                               component_mask);
+     component_mask[0] = true;
+     component_mask[1] = true;
+     component_mask[2] = false;
+     VectorTools::interpolate_boundary_values (dof_handler,
+                                               8,
+                                               EquationData::BoundaryValues<dim>(),
+                                               constraints_dirichlet_hanging_nodes,
+                                               component_mask);
+     constraints_dirichlet_hanging_nodes.close ();
+   }
  
-   const EquationData::RightHandSide<dim> right_hand_side;
-   std::vector<Vector<double> > right_hand_side_values (n_q_points,
-                                                        Vector<double>(dim));
-   std::vector<Vector<double> > right_hand_side_values_face (n_face_q_points,
-                                                           Vector<double>(dim));
+   template <int dim>
+   void PlasticityContactProblem<dim>::solve ()
+   {
+     pcout << "Solving ..." << std::endl;
+     Timer t;
  
-   FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
-   Vector<double>       cell_rhs (dofs_per_cell);
+     TrilinosWrappers::MPI::Vector    distributed_solution (system_rhs_newton);
+     distributed_solution = solution;
  
-   std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+     constraints_hanging_nodes.set_zero (distributed_solution);
  
-   typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
-                                                endc = dof_handler.end();
+     // Solving iterative
  
-   const FEValuesExtractors::Vector displacement (0);
+     MPI_Barrier (mpi_communicator);
+     t.restart();
  
-   TrilinosWrappers::MPI::Vector   test_rhs(solution);
-   const double kappa = 1.0;
-   for (; cell!=endc; ++cell)
-     if (cell->is_locally_owned())
-       {
-       fe_values.reinit (cell);
-       cell_matrix = 0;
-       cell_rhs = 0;
-       right_hand_side.vector_value_list (fe_values.get_quadrature_points(),
-                                          right_hand_side_values);
-       std::vector<SymmetricTensor<2,dim> > strain_tensor (n_q_points);
-       fe_values[displacement].get_function_symmetric_gradients (u, strain_tensor);
-       for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
-         {
-           SymmetricTensor<4,dim> stress_strain_tensor_linearized;
-           SymmetricTensor<4,dim> stress_strain_tensor;
-           SymmetricTensor<2,dim> stress_tensor;
-           plast_lin_hard->linearized_plast_linear_hardening (stress_strain_tensor_linearized,
-                                                              stress_strain_tensor,
-                                                              strain_tensor[q_point]);
-           //          if (q_point == 0)
-           //          std::cout<< stress_strain_tensor_linearized <<std::endl;
-           //          std::cout<< stress_strain_tensor <<std::endl;
-           for (unsigned int i=0; i<dofs_per_cell; ++i)
-             {
-               stress_tensor = stress_strain_tensor_linearized * plast_lin_hard->get_strain(fe_values, i, q_point);
-               for (unsigned int j=0; j<dofs_per_cell; ++j)
-                 {
-                   cell_matrix(i,j) += (stress_tensor *
-                                        plast_lin_hard->get_strain(fe_values, j, q_point) *
-                                        fe_values.JxW (q_point));
-                 }
-               // the linearized part a(v^i;v^i,v) of the rhs
-               cell_rhs(i) += (stress_tensor *
-                               strain_tensor[q_point] *
-                               fe_values.JxW (q_point));
-               // the residual part a(v^i;v) of the rhs
-               cell_rhs(i) -= (strain_tensor[q_point] * stress_strain_tensor *
-                               plast_lin_hard->get_strain(fe_values, i, q_point) *
-                               fe_values.JxW (q_point));
-               // the residual part F(v) of the rhs
-               Tensor<1,dim> rhs_values;
-               rhs_values = 0;
-               cell_rhs(i) += (fe_values[displacement].value (i, q_point) *
-                               rhs_values *
-                               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 () == 9)
-             {
-               fe_values_face.reinit (cell, face);
-               right_hand_side.vector_value_list (fe_values_face.get_quadrature_points(),
-                                                  right_hand_side_values_face);
-               for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
-                 {
-                   Tensor<1,dim> rhs_values;
-                   rhs_values = 0;
-                   for (unsigned int i=0; i<dofs_per_cell; ++i)
-                     cell_rhs(i) += (fe_values_face[displacement].value (i, q_point) *
-                                     rhs_values *
-                                     fe_values_face.JxW (q_point));
-                 }
-             }
-         }
-       cell->get_dof_indices (local_dof_indices);
-       constraints.distribute_local_to_global (cell_matrix, cell_rhs,
-                                               local_dof_indices,
-                                               system_matrix_newton, system_rhs_newton, true);
-       };
+     preconditioner_u.initialize (system_matrix_newton, additional_data);
  
-   system_matrix_newton.compress ();
-   system_rhs_newton.compress (Add);
- }
+     MPI_Barrier (mpi_communicator);
+     t.stop();
+     if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
+     run_time[6] += t.wall_time();
  
- template <int dim>
- void Step4<dim>::residual_nl_system (TrilinosWrappers::MPI::Vector &u,
-                                    Vector<double>                &sigma_eff_vector)
- {
-   QGauss<dim>  quadrature_formula(2);
-   QGauss<dim-1> face_quadrature_formula(2);
-   FEValues<dim> fe_values (fe, quadrature_formula,
-                          UpdateFlags(update_values    |
-                                      update_gradients |
-                                      update_q_points  |
-                                      update_JxW_values));
-   FEFaceValues<dim> fe_values_face (fe, face_quadrature_formula,
-                                   update_values   | update_quadrature_points |
-                                   update_JxW_values);
-   const unsigned int   dofs_per_cell = fe.dofs_per_cell;
-   const unsigned int   n_q_points    = quadrature_formula.size ();
-   const unsigned int   n_face_q_points = face_quadrature_formula.size();
-   const EquationData::RightHandSide<dim> right_hand_side;
-   std::vector<Vector<double> > right_hand_side_values (n_q_points,
-                                                        Vector<double>(dim));
-   std::vector<Vector<double> > right_hand_side_values_face (n_face_q_points,
-                                                           Vector<double>(dim));
-   Vector<double>       cell_rhs (dofs_per_cell);
-   Vector<double>       cell_sigma_eff (dofs_per_cell);
-   std::vector<unsigned int> local_dof_indices (dofs_per_cell);
-   const FEValuesExtractors::Vector displacement (0);
-   typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
-                                                endc = dof_handler.end();
-   unsigned int elast_points = 0;
-   unsigned int plast_points = 0;
-   double       sigma_eff = 0;
-   double       yield = 0;
-   unsigned int cell_number = 0;
-   for (; cell!=endc; ++cell)
-     if (cell->is_locally_owned())
-       {
-       fe_values.reinit (cell);
-       cell_rhs = 0;
-       right_hand_side.vector_value_list (fe_values.get_quadrature_points(),
-                                          right_hand_side_values);
-       std::vector<SymmetricTensor<2,dim> > strain_tensor (n_q_points);
-       fe_values[displacement].get_function_symmetric_gradients (u, strain_tensor);
-       for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
-         {
-           SymmetricTensor<4,dim> stress_strain_tensor;
-           SymmetricTensor<2,dim> stress_tensor;
-           plast_lin_hard->plast_linear_hardening (stress_strain_tensor, strain_tensor[q_point],
-                                                   elast_points, plast_points, sigma_eff, yield);
-           // sigma_eff_vector (cell_number) += sigma_eff;
-           sigma_eff_vector (cell_number) += yield;
-           /*  if (q_point == 0)
-               std::cout<< stress_strain_tensor <<std::endl;*/
-           for (unsigned int i=0; i<dofs_per_cell; ++i)
-             {
-               cell_rhs(i) -= (strain_tensor[q_point] * stress_strain_tensor * //(stress_tensor) *
-                               plast_lin_hard->get_strain(fe_values, i, q_point) *
-                               fe_values.JxW (q_point));
-               Tensor<1,dim> rhs_values;
-               rhs_values = 0;
-               cell_rhs(i) += ((fe_values[displacement].value (i, q_point) *
-                                rhs_values) *
-                               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 () == 9)
-             {
-               fe_values_face.reinit (cell, face);
-               right_hand_side.vector_value_list (fe_values_face.get_quadrature_points(),
-                                                  right_hand_side_values_face);
-               for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
-                 {
-                   Tensor<1,dim> rhs_values;
-                   rhs_values = 0;
-                   for (unsigned int i=0; i<dofs_per_cell; ++i)
-                     cell_rhs(i) += (fe_values_face[displacement].value (i, q_point) *
-                                     rhs_values *
-                                     fe_values_face.JxW (q_point));
-                 }
-             }
-         }
-       cell->get_dof_indices (local_dof_indices);
-       constraints_dirichlet_hanging_nodes.distribute_local_to_global (cell_rhs,
-                                                                       local_dof_indices,
-                                                                       system_rhs_newton);
-       sigma_eff_vector(cell_number) /= n_q_points;
-       cell_number += 1;
-       };
+     MPI_Barrier (mpi_communicator);
+     t.restart();
  
-   system_rhs_newton.compress ();
+   //   ReductionControl                 reduction_control (10000, 1e-15, 1e-4);
+   //   SolverCG<TrilinosWrappers::MPI::Vector>
+   //     solver (reduction_control, mpi_communicator);
+   //   solver.solve (system_matrix_newton, distributed_solution, system_rhs_newton, preconditioner_u);
  
-   unsigned int sum_elast_points = Utilities::MPI::sum(elast_points, mpi_communicator);
-   unsigned int sum_plast_points = Utilities::MPI::sum(plast_points, mpi_communicator);
-   pcout<< "Elast-Points = " << sum_elast_points <<std::endl;
-   pcout<< "Plast-Points = " << sum_plast_points <<std::endl;
- }
+     PrimitiveVectorMemory<TrilinosWrappers::MPI::Vector> mem;
+     TrilinosWrappers::MPI::Vector    tmp (system_rhs_newton);
+     const double solver_tolerance = 1e-4 *
+           system_matrix_newton.residual (tmp, distributed_solution, system_rhs_newton);
+     SolverControl solver_control (system_matrix_newton.m(), solver_tolerance);
+     SolverFGMRES<TrilinosWrappers::MPI::Vector>
+        solver(solver_control, mem,
+        SolverFGMRES<TrilinosWrappers::MPI::Vector>::
+        AdditionalData(30, true));
+     solver.solve(system_matrix_newton, distributed_solution, system_rhs_newton, preconditioner_u);
  
-                                  // @sect4{Step4::projection_active_set}
+     pcout << "Initial error: " << solver_control.initial_value() <<std::endl;
+     pcout << "   " << solver_control.last_step()
+               << " CG iterations needed to obtain convergence with an error: "
+               <<  solver_control.last_value()
+               << std::endl;
  
-                                // Projection and updating of the active set
-                                  // for the dofs which penetrates the obstacle.
- template <int dim>
- void Step4<dim>::projection_active_set ()
- {
-   clock_t                        start_proj, end_proj;
-   
-   const EquationData::Obstacle<dim>     obstacle;
-   std::vector<bool>                     vertex_touched (dof_handler.n_dofs (), false);
-   
-   typename DoFHandler<dim>::active_cell_iterator
-   cell = dof_handler.begin_active(),
-   endc = dof_handler.end();
-   TrilinosWrappers::MPI::Vector         distributed_solution (system_rhs_newton);
-   distributed_solution = solution;
-   TrilinosWrappers::MPI::Vector         lambda (solution);
-   lambda = resid_vector;
-   TrilinosWrappers::MPI::Vector         diag_mass_matrix_vector_relevant (solution);
-   diag_mass_matrix_vector_relevant = diag_mass_matrix_vector;
-     
-   constraints.reinit(locally_relevant_dofs);
-   active_set.clear ();
-   IndexSet     active_set_locally_owned;
-   active_set_locally_owned.set_size (locally_owned_dofs.size ());
-   const double c = 100.0*e_modul;
-   for (; cell!=endc; ++cell)
-     if (cell->is_locally_owned())
-       for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-       if (cell->face (face)->at_boundary()
-           && cell->face (face)->boundary_indicator () == 9)
-         for (unsigned int v=0; v<GeometryInfo<dim-1>::vertices_per_cell; ++v)
-           {
-             unsigned int index_z = cell->face (face)->vertex_dof_index (v,2);
-             if (vertex_touched[cell->face (face)->vertex_index(v)] == false)
-               vertex_touched[cell->face (face)->vertex_index(v)] = true;
-             else
-               continue;
-             // the local row where
-             Point<dim> point (cell->face (face)->vertex (v)[0],/* + solution (index_x),*/
-                               cell->face (face)->vertex (v)[1],
-                               cell->face (face)->vertex (v)[2]);
-             double obstacle_value = obstacle.value (point, 2);
-             double solution_index_z = solution (index_z);
-             double gap = obstacle_value - point (2);
-             if (lambda (index_z) +
-                 c*diag_mass_matrix_vector_relevant (index_z)*(solution_index_z - gap) > 0)
-               {
-                 constraints.add_line (index_z);
-                 constraints.set_inhomogeneity (index_z, gap);
-                 distributed_solution (index_z) = gap;
-                 if (locally_relevant_dofs.is_element (index_z))
-                   active_set.add_index (index_z);
-                 if (locally_owned_dofs.is_element (index_z))
-                   active_set_locally_owned.add_index (index_z);
-                 // std::cout<< index_z << ", "
-                 //       << "Error: " << lambda (index_z) +
-                 //   diag_mass_matrix_vector_relevant (index_z)*c*(solution_index_z - gap)
-                 //       << ", " << lambda (index_z)
-                 //       << ", " << diag_mass_matrix_vector_relevant (index_z)
-                 //       << ", " << obstacle_value
-                 //       << ", " << solution_index_z
-                 //       <<std::endl;
-               }
-           }
-   distributed_solution.compress(Insert);
-   unsigned int sum_contact_constraints = Utilities::MPI::sum(active_set_locally_owned.n_elements (),
-                                                            mpi_communicator);
-   pcout << "Number of Contact-Constaints: " << sum_contact_constraints <<std::endl;
-   solution = distributed_solution;
-   constraints.close ();
-   const ConstraintMatrix::MergeConflictBehavior
-     merge_conflict_behavior = ConstraintMatrix::left_object_wins;
-   constraints.merge (constraints_dirichlet_hanging_nodes, merge_conflict_behavior);
- }
+     MPI_Barrier (mpi_communicator);
+     t.stop();
+     if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
+     run_time[7] += t.wall_time();
  
- template <int dim>
- void Step4<dim>::dirichlet_constraints ()
- {
-   /* boundary_indicators:
-             _______
-            /  9    /|
-           /______ / |
-         8|       | 8|
-          |   8   | /
-          |_______|/
-              6
-    */
-   constraints_dirichlet_hanging_nodes.reinit (locally_relevant_dofs);
-   constraints_dirichlet_hanging_nodes.merge (constraints_hanging_nodes);
-   std::vector<bool> component_mask (dim, true);
-   component_mask[0] = true;
-   component_mask[1] = true;
-   component_mask[2] = true;
-   VectorTools::interpolate_boundary_values (dof_handler,
-                                           6,
-                                           EquationData::BoundaryValues<dim>(),
-                                           constraints_dirichlet_hanging_nodes,
-                                           component_mask);
-   component_mask[0] = true;
-   component_mask[1] = true;
-   component_mask[2] = false;
-   VectorTools::interpolate_boundary_values (dof_handler,
-                                           8,
-                                           EquationData::BoundaryValues<dim>(),
-                                           constraints_dirichlet_hanging_nodes,
-                                           component_mask);
-   constraints_dirichlet_hanging_nodes.close ();
- }
+     number_iterations += solver_control.last_step();
  
- template <int dim>
- void Step4<dim>::solve ()
- {
-   pcout << "Solving ..." << std::endl;
-   Timer t;
-   TrilinosWrappers::MPI::Vector    distributed_solution (system_rhs_newton);
-   distributed_solution = solution;
-   
-   constraints_hanging_nodes.set_zero (distributed_solution);
-   // Solving iterative
-   MPI_Barrier (mpi_communicator);
-   t.restart();
-   preconditioner_u.initialize (system_matrix_newton, additional_data);
-   
-   MPI_Barrier (mpi_communicator);
-   t.stop();
-   if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
-   run_time[6] += t.wall_time();
-   
-   MPI_Barrier (mpi_communicator);
-   t.restart();
- //   ReductionControl                 reduction_control (10000, 1e-15, 1e-4);
- //   SolverCG<TrilinosWrappers::MPI::Vector>
- //     solver (reduction_control, mpi_communicator);
- //   solver.solve (system_matrix_newton, distributed_solution, system_rhs_newton, preconditioner_u);
-   PrimitiveVectorMemory<TrilinosWrappers::MPI::Vector> mem;
-   TrilinosWrappers::MPI::Vector    tmp (system_rhs_newton);
-   const double solver_tolerance = 1e-4 *
-       system_matrix_newton.residual (tmp, distributed_solution, system_rhs_newton);
-   SolverControl solver_control (system_matrix_newton.m(), solver_tolerance);
-   SolverFGMRES<TrilinosWrappers::MPI::Vector>
-      solver(solver_control, mem,
-      SolverFGMRES<TrilinosWrappers::MPI::Vector>::
-      AdditionalData(30, true));
-   solver.solve(system_matrix_newton, distributed_solution, system_rhs_newton, preconditioner_u);          
-   pcout << "Initial error: " << solver_control.initial_value() <<std::endl;
-   pcout << "   " << solver_control.last_step()
-           << " CG iterations needed to obtain convergence with an error: "
-           <<  solver_control.last_value()
-           << std::endl;
-   MPI_Barrier (mpi_communicator);
-   t.stop();
-   if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
-   run_time[7] += t.wall_time();
-   number_iterations += solver_control.last_step();
-   constraints.distribute (distributed_solution);
-   solution = distributed_solution;
- }
- template <int dim>
- void Step4<dim>::solve_newton ()
- {
-   double                         resid=0;
-   double                         resid_old=100000;
-   TrilinosWrappers::MPI::Vector  res (system_rhs_newton);
-   TrilinosWrappers::MPI::Vector  tmp_vector (system_rhs_newton);
-   Timer                          t;
-   std::vector<std::vector<bool> > constant_modes;
-   std::vector<bool>  components (dim,true);
-   components[dim] = false;
-   DoFTools::extract_constant_modes (dof_handler, components,
-                                   constant_modes);
-   additional_data.elliptic = true;
-   additional_data.n_cycles = 1;
-   additional_data.w_cycle = false;
-   additional_data.output_details = false;
-   additional_data.smoother_sweeps = 2;
-   additional_data.aggregation_threshold = 1e-2;
-   IndexSet                            active_set_old (active_set);
-   Vector<double>                      sigma_eff_vector;
-   sigma_eff_vector.reinit (triangulation.n_active_cells());  
-   unsigned int j = 0;
-   unsigned int number_assemble_system = 0;
-   for (; j<=100;j++)
-     {
-       pcout<< " " <<std::endl;
-       pcout<< j << ". Iteration of the inexact Newton-method." <<std::endl;
-       pcout<< "Update of active set" <<std::endl;
-       
-       MPI_Barrier (mpi_communicator);
-       t.restart(); 
-       
-       projection_active_set ();
-       
-       MPI_Barrier (mpi_communicator);
-       t.stop();
-       if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
-         run_time[5] += t.wall_time();
-       
-       pcout<< "Assembling ... " <<std::endl;
-       MPI_Barrier (mpi_communicator);
-       t.restart(); 
-       system_matrix_newton = 0;
-       system_rhs_newton = 0;
-       assemble_nl_system (solution);  //compute Newton-Matrix
-       MPI_Barrier (mpi_communicator);
-       t.stop();
-       if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
-         run_time[1] += t.wall_time();
-       number_assemble_system += 1;
-       MPI_Barrier (mpi_communicator);
-       t.restart();
-       solve ();
-       MPI_Barrier (mpi_communicator);
-       t.stop();
-       if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
-         run_time[2] += t.wall_time();
-       
-       TrilinosWrappers::MPI::Vector    distributed_solution (system_rhs_newton);
-       distributed_solution = solution;
-       int damped = 0;
-       tmp_vector = old_solution;
-       double a = 0;
-       for (unsigned int i=0; (i<10)&&(!damped); i++)
-       {
-         a=pow(0.5, static_cast<double>(i));
-         old_solution = tmp_vector;
-         old_solution.sadd(1-a,a, distributed_solution);
-         MPI_Barrier (mpi_communicator);
-           t.restart();
-         system_rhs_newton = 0;
-         sigma_eff_vector = 0;
-         solution = old_solution;
-         residual_nl_system (solution, sigma_eff_vector);
-         res = system_rhs_newton;
-         const unsigned int
-           start_res     = (res.local_range().first),
-           end_res       = (res.local_range().second);
-         for (unsigned int n=start_res; n<end_res; ++n)
-           if (constraints.is_inhomogeneously_constrained (n))
-             {
-               // pcout<< i << ". " << constraints.get_inhomogeneity (n)
-               //       << ". " << res (n)
-               //       << ", start = " << start_res
-               //       << ", end = " << end_res
-               //       <<std::endl;
-               res(n) = 0;
-             }
-         resid = res.l2_norm ();
-         pcout<< "Residual: " << resid <<std::endl;
-         if (resid<resid_old)
-           {
-             pcout<< "Newton-damping parameter alpha = " << a <<std::endl;
-             damped=1;
-           }
-         MPI_Barrier (mpi_communicator);
-           t.stop();
-           if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
-             run_time[3] += t.wall_time();
-       }
-       if (resid<1e-8)
-       {
-         pcout<< "Inexact Newton-method stopped with residual = " << resid <<std::endl;
-         pcout<< "Number of Assembling systems = " << number_assemble_system <<std::endl;
-         break;
-       }
-       resid_old=resid;
-       resid_vector = system_rhs_newton;
-       if (active_set == active_set_old && resid < 1e-10)
-       break;
-       active_set_old = active_set;
-     } // End of active-set-loop
-   
-   pcout<< "Creating output." <<std::endl;
-   MPI_Barrier (mpi_communicator);
-   t.restart();
-   std::ostringstream filename_solution;
-   filename_solution << "solution";
-   // filename_solution << "solution_";
-   // filename_solution << k;
-   output_results (filename_solution.str ());
-   // output_results (sigma_eff_vector, "sigma_eff");
-   MPI_Barrier (mpi_communicator);
-   t.stop();
-   if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
-     run_time[4] += t.wall_time();
-   pcout<< "Number of Solver-Iterations = " << number_iterations <<std::endl;
-   pcout<< "%%%%%% Rechenzeit make grid and setup = " << run_time[0] <<std::endl;
-   pcout<< "%%%%%% Rechenzeit projection active set = " << run_time[5] <<std::endl;
-   pcout<< "%%%%%% Rechenzeit assemble system = " << run_time[1] <<std::endl;
-   pcout<< "%%%%%% Rechenzeit solve system = " << run_time[2] <<std::endl;
-   pcout<< "%%%%%% Rechenzeit preconditioner = " << run_time[6] <<std::endl;
-   pcout<< "%%%%%% Rechenzeit solve with CG = " << run_time[7] <<std::endl;
-   pcout<< "%%%%%% Rechenzeit error and lambda = " << run_time[3] <<std::endl;
-   pcout<< "%%%%%% Rechenzeit output = " << run_time[4] <<std::endl;
- }
+     constraints.distribute (distributed_solution);
  
- template <int dim>
- void Step4<dim>::output_results (const std::string& title) const
- {
-   move_mesh (solution);
+     solution = distributed_solution;
+   }
  
-   TrilinosWrappers::MPI::Vector         lambda (solution);
-   lambda = resid_vector;
+   template <int dim>
+   void PlasticityContactProblem<dim>::solve_newton ()
+   {
+     double                         resid=0;
+     double                         resid_old=100000;
+     TrilinosWrappers::MPI::Vector  res (system_rhs_newton);
+     TrilinosWrappers::MPI::Vector  tmp_vector (system_rhs_newton);
+     Timer                          t;
+     std::vector<std::vector<bool> > constant_modes;
+     std::vector<bool>  components (dim,true);
+     components[dim] = false;
+     DoFTools::extract_constant_modes (dof_handler, components,
+                                       constant_modes);
+     additional_data.elliptic = true;
+     additional_data.n_cycles = 1;
+     additional_data.w_cycle = false;
+     additional_data.output_details = false;
+     additional_data.smoother_sweeps = 2;
+     additional_data.aggregation_threshold = 1e-2;
+     IndexSet                            active_set_old (active_set);
+     unsigned int j = 0;
+     unsigned int number_assemble_system = 0;
+     for (; j<=100;j++)
+       {
+         pcout<< " " <<std::endl;
+         pcout<< j << ". Iteration of the inexact Newton-method." <<std::endl;
+         pcout<< "Update of active set" <<std::endl;
+         MPI_Barrier (mpi_communicator);
+         t.restart();
+         update_solution_and_constraints ();
+         MPI_Barrier (mpi_communicator);
+         t.stop();
+         if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
+           run_time[5] += t.wall_time();
+         pcout<< "Assembling ... " <<std::endl;
+         MPI_Barrier (mpi_communicator);
+         t.restart();
+         system_matrix_newton = 0;
+         system_rhs_newton = 0;
+         assemble_nl_system (solution);  //compute Newton-Matrix
+         MPI_Barrier (mpi_communicator);
+         t.stop();
+         if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
+           run_time[1] += t.wall_time();
+         number_assemble_system += 1;
+         MPI_Barrier (mpi_communicator);
+         t.restart();
+         solve ();
+         MPI_Barrier (mpi_communicator);
+         t.stop();
+         if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
+           run_time[2] += t.wall_time();
+         TrilinosWrappers::MPI::Vector    distributed_solution (system_rhs_newton);
+         distributed_solution = solution;
+         int damped = 0;
+         tmp_vector = old_solution;
+         double a = 0;
+         for (unsigned int i=0; (i<10)&&(!damped); i++)
+           {
+             a=pow(0.5, static_cast<double>(i));
+             old_solution = tmp_vector;
+             old_solution.sadd(1-a,a, distributed_solution);
+             MPI_Barrier (mpi_communicator);
+             t.restart();
+             system_rhs_newton = 0;
+             solution = old_solution;
+             residual_nl_system (solution);
+             res = system_rhs_newton;
+             const unsigned int
+               start_res     = (res.local_range().first),
+               end_res       = (res.local_range().second);
+             for (unsigned int n=start_res; n<end_res; ++n)
+               if (constraints.is_inhomogeneously_constrained (n))
+                 {
+                   // pcout<< i << ". " << constraints.get_inhomogeneity (n)
+                   //   << ". " << res (n)
+                   //   << ", start = " << start_res
+                   //   << ", end = " << end_res
+                   //   <<std::endl;
+                   res(n) = 0;
+                 }
+             resid = res.l2_norm ();
+             pcout<< "Residual: " << resid <<std::endl;
+             if (resid<resid_old)
+               {
+                 pcout<< "Newton-damping parameter alpha = " << a <<std::endl;
+                 damped=1;
+               }
+             MPI_Barrier (mpi_communicator);
+             t.stop();
+             if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
+               run_time[3] += t.wall_time();
+           }
+         if (resid<1e-8)
+           {
+             pcout<< "Inexact Newton-method stopped with residual = " << resid <<std::endl;
+             pcout<< "Number of Assembling systems = " << number_assemble_system <<std::endl;
+             break;
+           }
+         resid_old=resid;
+         resid_vector = system_rhs_newton;
+         if (active_set == active_set_old && resid < 1e-10)
+           break;
+         active_set_old = active_set;
+       } // End of active-set-loop
+     pcout<< "Creating output." <<std::endl;
+     MPI_Barrier (mpi_communicator);
+     t.restart();
+     std::ostringstream filename_solution;
+     filename_solution << "solution";
+     // filename_solution << "solution_";
+     // filename_solution << k;
+     output_results (filename_solution.str ());
+     MPI_Barrier (mpi_communicator);
+     t.stop();
+     if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
+       run_time[4] += t.wall_time();
+     pcout<< "Number of Solver-Iterations = " << number_iterations <<std::endl;
+     pcout<< "%%%%%% Rechenzeit make grid and setup = " << run_time[0] <<std::endl;
+     pcout<< "%%%%%% Rechenzeit projection active set = " << run_time[5] <<std::endl;
+     pcout<< "%%%%%% Rechenzeit assemble system = " << run_time[1] <<std::endl;
+     pcout<< "%%%%%% Rechenzeit solve system = " << run_time[2] <<std::endl;
+     pcout<< "%%%%%% Rechenzeit preconditioner = " << run_time[6] <<std::endl;
+     pcout<< "%%%%%% Rechenzeit solve with CG = " << run_time[7] <<std::endl;
+     pcout<< "%%%%%% Rechenzeit error and lambda = " << run_time[3] <<std::endl;
+     pcout<< "%%%%%% Rechenzeit output = " << run_time[4] <<std::endl;
+   }
  
-   DataOut<dim> data_out;
+   template <int dim>
+   void PlasticityContactProblem<dim>::output_results (const std::string& title) const
+   {
+     move_mesh (solution);
  
-   data_out.attach_dof_handler (dof_handler);
+     TrilinosWrappers::MPI::Vector         lambda (solution);
+     lambda = resid_vector;
  
-   const std::vector<DataComponentInterpretation::DataComponentInterpretation>
-     data_component_interpretation
-     (dim, DataComponentInterpretation::component_is_part_of_vector);
-   data_out.add_data_vector (solution, std::vector<std::string>(dim, "Displacement"),
-                           DataOut<dim>::type_dof_data,
-                           data_component_interpretation);
-   data_out.add_data_vector (lambda, std::vector<std::string>(dim, "Residual"),
-                           DataOut<dim>::type_dof_data,
-                           data_component_interpretation);
-   data_out.add_data_vector (active_set, std::vector<std::string>(dim, "ActiveSet"),
-                           DataOut<dim>::type_dof_data,
-                           data_component_interpretation);
+     DataOut<dim> data_out;
  
-   Vector<float> subdomain (triangulation.n_active_cells());
-   for (unsigned int i=0; i<subdomain.size(); ++i)
-     subdomain(i) = triangulation.locally_owned_subdomain();
-   data_out.add_data_vector (subdomain, "subdomain");
+     data_out.attach_dof_handler (dof_handler);
  
-   data_out.build_patches ();
+     const std::vector<DataComponentInterpretation::DataComponentInterpretation>
+       data_component_interpretation
+       (dim, DataComponentInterpretation::component_is_part_of_vector);
+     data_out.add_data_vector (solution, std::vector<std::string>(dim, "Displacement"),
+                               DataOut<dim>::type_dof_data,
+                               data_component_interpretation);
+     data_out.add_data_vector (lambda, std::vector<std::string>(dim, "Residual"),
+                               DataOut<dim>::type_dof_data,
+                               data_component_interpretation);
+     data_out.add_data_vector (active_set, std::vector<std::string>(dim, "ActiveSet"),
+                               DataOut<dim>::type_dof_data,
+                               data_component_interpretation);
  
-   const std::string filename = (title + "-" +
-                               Utilities::int_to_string
-                               (triangulation.locally_owned_subdomain(), 4));
+     Vector<float> subdomain (triangulation.n_active_cells());
+     for (unsigned int i=0; i<subdomain.size(); ++i)
+       subdomain(i) = triangulation.locally_owned_subdomain();
+     data_out.add_data_vector (subdomain, "subdomain");
  
-   std::ofstream output_vtu ((filename + ".vtu").c_str ());
-   data_out.write_vtu (output_vtu);
+     data_out.build_patches ();
  
-   if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
-     {
-       std::vector<std::string> filenames;
-       for (unsigned int i=0;
-          i<Utilities::MPI::n_mpi_processes(mpi_communicator);
-          ++i)
-       filenames.push_back ("solution-" +
-                            Utilities::int_to_string (i, 4) +
-                            ".vtu");
-       std::ofstream master_output ((filename + ".pvtu").c_str());
-       data_out.write_pvtu_record (master_output, filenames);
-     }
+     const std::string filename = (title + "-" +
+                                   Utilities::int_to_string
+                                   (triangulation.locally_owned_subdomain(), 4));
  
-   TrilinosWrappers::MPI::Vector  tmp (solution);
-   tmp *= -1;
-   move_mesh (tmp);
- }
+     std::ofstream output_vtu ((filename + ".vtu").c_str ());
+     data_out.write_vtu (output_vtu);
  
- template <int dim>
- void Step4<dim>::move_mesh (const TrilinosWrappers::MPI::Vector &_complete_displacement) const
- {
-   pcout<< "Moving mesh." <<std::endl;
-   std::vector<bool> vertex_touched (triangulation.n_vertices(),
-                                   false);
-   for (typename DoFHandler<dim>::active_cell_iterator
-        cell = dof_handler.begin_active ();
-        cell != dof_handler.end(); ++cell)
-     if (cell->is_locally_owned())
-       for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
-       {
-         if (vertex_touched[cell->vertex_index(v)] == false)
-           {
-             vertex_touched[cell->vertex_index(v)] = true;
-             Point<dim> vertex_displacement;
-             for (unsigned int d=0; d<dim; ++d)
-               {
-                 if (_complete_displacement(cell->vertex_dof_index(v,d)) != 0)
-                   vertex_displacement[d]
-                     = _complete_displacement(cell->vertex_dof_index(v,d));
-               }
-             cell->vertex(v) += vertex_displacement;
-           }
-       }
- }
+     if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
+       {
+         std::vector<std::string> filenames;
+         for (unsigned int i=0;
+              i<Utilities::MPI::n_mpi_processes(mpi_communicator);
+              ++i)
+           filenames.push_back ("solution-" +
+                                Utilities::int_to_string (i, 4) +
+                                ".vtu");
+         std::ofstream master_output ((filename + ".pvtu").c_str());
+         data_out.write_pvtu_record (master_output, filenames);
+       }
  
- template <int dim>
- void Step4<dim>::output_results (TrilinosWrappers::MPI::Vector vector,
-                                const std::string& title) const
- {
-   DataOut<dim> data_out;
+     TrilinosWrappers::MPI::Vector  tmp (solution);
+     tmp *= -1;
+     move_mesh (tmp);
+   }
  
-   data_out.attach_dof_handler (dof_handler);
-   data_out.add_data_vector (vector, "vector_to_plot");
+   template <int dim>
+   void PlasticityContactProblem<dim>::move_mesh (const TrilinosWrappers::MPI::Vector &_complete_displacement) const
+   {
+     pcout<< "Moving mesh." <<std::endl;
+     std::vector<bool> vertex_touched (triangulation.n_vertices(),
+                                       false);
+     for (typename DoFHandler<dim>::active_cell_iterator
+            cell = dof_handler.begin_active ();
+          cell != dof_handler.end(); ++cell)
+       if (cell->is_locally_owned())
+         for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+           {
+             if (vertex_touched[cell->vertex_index(v)] == false)
+               {
+                 vertex_touched[cell->vertex_index(v)] = true;
+                 Point<dim> vertex_displacement;
+                 for (unsigned int d=0; d<dim; ++d)
+                   {
+                     if (_complete_displacement(cell->vertex_dof_index(v,d)) != 0)
+                       vertex_displacement[d]
+                         = _complete_displacement(cell->vertex_dof_index(v,d));
+                   }
+                 cell->vertex(v) += vertex_displacement;
+               }
+           }
+   }
  
-   data_out.build_patches ();
+   template <int dim>
+   void PlasticityContactProblem<dim>::output_results (TrilinosWrappers::MPI::Vector vector,
+                                    const std::string& title) const
+   {
+     DataOut<dim> data_out;
  
-   std::ofstream output_vtk (dim == 2 ?
-                           (title + ".vtk").c_str () :
-                           (title + ".vtk").c_str ());
-   data_out.write_vtk (output_vtk);
- }
+     data_out.attach_dof_handler (dof_handler);
+     data_out.add_data_vector (vector, "vector_to_plot");
  
- template <int dim>
- void Step4<dim>::output_results (Vector<double> vector, const std::string& title) const
- {
-   DataOut<dim> data_out;
+     data_out.build_patches ();
  
-   data_out.attach_dof_handler (dof_handler);
-   data_out.add_data_vector (vector, "vector_to_plot");
+     std::ofstream output_vtk (dim == 2 ?
+                               (title + ".vtk").c_str () :
+                               (title + ".vtk").c_str ());
+     data_out.write_vtk (output_vtk);
+   }
  
-   data_out.build_patches ();
+   template <int dim>
+   void PlasticityContactProblem<dim>::run ()
+   {
+     pcout << "Solving problem in " << dim << " space dimensions." << std::endl;
 -  
 +
-   std::ofstream output_vtk (dim == 2 ?
-                           (title + ".vtk").c_str () :
-                           (title + ".vtk").c_str ());
-   data_out.write_vtk (output_vtk);
- }
+     run_time.resize (8);
 -  
 +
- template <int dim>
- void Step4<dim>::run ()
- {
-   pcout << "Solving problem in " << dim << " space dimensions." << std::endl;
+     clock_t     start, end;
 -  
 +
-   run_time.resize (8);
+     start = clock();
+     make_grid();
+     //  mesh_surface ();
  
-   clock_t     start, end;
+     setup_system ();
 -    
 +
-   start = clock();
-   make_grid();
-   //  mesh_surface ();
-   
-   setup_system (); 
-   
-   end = clock();
-   run_time[0] = (double)(end-start)/CLOCKS_PER_SEC;
+     end = clock();
+     run_time[0] = (double)(end-start)/CLOCKS_PER_SEC;
  
-   solve_newton ();
+     solve_newton ();
+   }
  }
  
                                   // @sect3{The <code>main</code> function}
  
-                                // And this is the main function. It also
-                                // looks mostly like in step-3, but if you
-                                // look at the code below, note how we first
-                                // create a variable of type
-                                // <code>Step4@<2@></code> (forcing
-                                // the compiler to compile the class template
-                                // with <code>dim</code> replaced by
-                                // <code>2</code>) and run a 2d simulation,
-                                // and then we do the whole thing over in 3d.
-                                //
-                                // In practice, this is probably not what you
-                                // would do very frequently (you probably
-                                // either want to solve a 2d problem, or one
-                                // in 3d, but not both at the same
-                                // time). However, it demonstrates the
-                                // mechanism by which we can simply change
-                                // which dimension we want in a single place,
-                                // and thereby force the compiler to
-                                // recompile the dimension independent class
-                                // templates for the dimension we
-                                // request. The emphasis here lies on the
-                                // fact that we only need to change a single
-                                // place. This makes it rather trivial to
-                                // debug the program in 2d where computations
-                                // are fast, and then switch a single place
-                                // to a 3 to run the much more computing
-                                // intensive program in 3d for `real'
-                                // computations.
-                                //
-                                // Each of the two blocks is enclosed in
-                                // braces to make sure that the
-                                // <code>laplace_problem_2d</code> variable
-                                // goes out of scope (and releases the memory
-                                // it holds) before we move on to allocate
-                                // memory for the 3d case. Without the
-                                // additional braces, the
-                                // <code>laplace_problem_2d</code> variable
-                                // would only be destroyed at the end of the
-                                // function, i.e. after running the 3d
-                                // problem, and would needlessly hog memory
-                                // while the 3d run could actually use it.
-                                  //
-                                  // Finally, the first line of the function is
-                                  // used to suppress some output.  Remember
-                                  // that in the previous example, we had the
-                                  // output from the linear solvers about the
-                                  // starting residual and the number of the
-                                  // iteration where convergence was
-                                  // detected. This can be suppressed through
-                                  // the <code>deallog.depth_console(0)</code>
-                                  // call.
-                                  //
-                                  // The rationale here is the following: the
-                                  // deallog (i.e. deal-log, not de-allog)
-                                  // variable represents a stream to which some
-                                  // parts of the library write output. It
-                                  // redirects this output to the console and
-                                  // if required to a file. The output is
-                                  // nested in a way so that each function can
-                                  // use a prefix string (separated by colons)
-                                  // for each line of output; if it calls
-                                  // another function, that may also use its
-                                  // prefix which is then printed after the one
-                                  // of the calling function. Since output from
-                                  // functions which are nested deep below is
-                                  // usually not as important as top-level
-                                  // output, you can give the deallog variable
-                                  // a maximal depth of nested output for
-                                  // output to console and file. The depth zero
-                                  // which we gave here means that no output is
-                                  // written. By changing it you can get more
-                                  // information about the innards of the
-                                  // library.
  int main (int argc, char *argv[])
  {
+   using namespace dealii;
+   using namespace Step42;
    deallog.depth_console (0);
  
    clock_t     start, end;
index 80997dd2c1f21854328946867e996bb82c936af7,be5ecda5383f4e5b5041baa8b4da414f91a77a97..9d6bdabbbbcc27ef5aecace97abbde968b2a011c
@@@ -63,60 -69,20 +69,60 @@@ namespace type
                                      * the @ref distributed module for
                                      * more information.
                                      */
-   const unsigned int artificial_subdomain_id = static_cast<subdomain_id_t>(-2);
+   const unsigned int artificial_subdomain_id = static_cast<subdomain_id>(-2);
  
 -                                   /**
 -                                    * The type used to denote global dof
 -                                    * indices.
 -                                    */
 +
 +#ifdef DEAL_II_USE_LARGE_INDEX_TYPE
 +                                 /**
 +                                  * The type used for global indices of
 +                                  * degrees of freedom. While in sequential
 +                                  * computations the 4 billion indices of
 +                                  * 32-bit unsigned integers is plenty,
 +                                  * parallel computations using the
 +                                  * parallel::distributed::Triangulation
 +                                  * class can overflow this number and we
 +                                  * need a bigger index space.
 +                                  *
 +                                  * The data type always indicates an
 +                                  * unsigned integer type.
 +                                  */
 +  typedef unsigned long long int global_dof_index;
 +
 +                                 /**
 +                                  * An identifier that denotes the MPI type
 +                                  * associated with types::global_dof_index.
 +                                  */
 +#  define DEAL_II_DOF_INDEX_MPI_TYPE MPI_UNSIGNED_LONG_LONG
 +#else
 +                                 /**
 +                                  * The type used for global indices of
 +                                  * degrees of freedom. While in sequential
 +                                  * computations the 4 billion indices of
 +                                  * 32-bit unsigned integers is plenty,
 +                                  * parallel computations using the
 +                                  * parallel::distributed::Triangulation
 +                                  * class can overflow this number and we
 +                                  * need a bigger index space.
 +                                  *
 +                                  * The data type always indicates an
 +                                  * unsigned integer type.
 +                                  */
    typedef unsigned int global_dof_index;
  
 -                                   /**
 -                                    * An invalid value for indices of degrees
 -                                    * of freedom.
 -                                    */
 +                                 /**
 +                                  * An identifier that denotes the MPI type
 +                                  * associated with types::global_dof_index.
 +                                  */
 +#  define DEAL_II_DOF_INDEX_MPI_TYPE MPI_UNSIGNED
 +#endif
 +
 +                                 /**
 +                                  * An invalid value for indices of degrees
 +                                  * of freedom.
 +                                  */
    const global_dof_index invalid_dof_index = static_cast<global_dof_index>(-1);
  
 +
            /**
             * The type used to denote boundary indicators associated with every
             * piece of the boundary and, in the case of meshes that describe
             * Invalid material_id which we need in several places as a default value.
             * We assume that all material_ids lie in the range [0, invalid_material_id).
             */
-   const material_id_t invalid_material_id = static_cast<material_id_t>(-1);
+   const material_id invalid_material_id = static_cast<material_id>(-1);
  
 +
 +
  }
  
  
index 6ae454508ed069fd0c00be916395de6cbe49d882,b236cdea757acb5654cb529a106d4b77170a4dc8..d352b8d2f5a474b5b218f90eddb8bdd2a0c8f1c8
@@@ -318,9 -318,7 +318,9 @@@ class DoFHandler  :  public Subscripto
                                        * object to the finite element.
                                        */
      virtual void distribute_dofs (const FiniteElement<dim,spacedim> &fe,
 -                                  const unsigned int        offset = 0);
 +                                                        const types::global_dof_index     offset = 0);
-     
++
 +    virtual void distribute_mg_dofs (const FiniteElement<dim, spacedim>& fe, const types::global_dof_index offset = 0);
  
                                       /**
                                        * After distribute_dofs() with
                                        * DoFs which are constrained by
                                        * hanging nodes, see @ref constraints.
                                        */
 -    unsigned int n_dofs () const;
 +    types::global_dof_index n_dofs () const;
-     
++
 +    types::global_dof_index n_dofs (const unsigned int level) const;
  
                                       /**
                                        * Return the number of degrees of freedom
                                        * indicators under
                                        * consideration.
                                        */
 -    unsigned int
 +    types::global_dof_index
-     n_boundary_dofs (const std::set<types::boundary_id_t> &boundary_indicators) const;
+     n_boundary_dofs (const std::set<types::boundary_id> &boundary_indicators) const;
  
                                       /**
                                        * Access to an object informing
                                        */
      DoFHandler & operator = (const DoFHandler &);
  
-       
 +    class MGVertexDoFs {
 +      private:
 +        unsigned int coarsest_level;
 +        unsigned int finest_level;
 +        types::global_dof_index* indices;
 +        types::global_dof_index* indices_offset;
-     
++
 +      public:
 +        DeclException0 (ExcNoMemory);
 +        MGVertexDoFs ();
 +        ~MGVertexDoFs ();
 +        unsigned int get_coarsest_level () const;
 +        unsigned int get_finest_level () const;
 +        types::global_dof_index get_index (const unsigned int level, const unsigned int dof_number) const;
 +        void init (const unsigned int coarsest_level, const unsigned int finest_level, const unsigned int dofs_per_vertex);
 +        void set_index (const unsigned int level, const unsigned int dof_number, const types::global_dof_index index);
 +    };
++
 +    void clear_mg_space ();
 +
                                       /**
                                        * Free all used memory.
                                        */
      void clear_space ();
-     
 +    void reserve_space ();
-     
++
 +    template <int structdim>
 +    types::global_dof_index get_dof_index (const unsigned int obj_level, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index) const;
 +
 +    template<int structdim>
 +    void set_dof_index (const unsigned int obj_level, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index) const;
 +
                                       /**
                                        * Space to store the DoF numbers
                                        * for the different
                                        * Triangulation objects.
                                        */
      std::vector<dealii::internal::DoFHandler::DoFLevel<dim>*> levels;
-     
 +    std::vector<dealii::internal::DoFHandler::DoFLevel<dim>*> mg_levels;
 +
                                       /**
                                        * Space to store DoF numbers of
                                        * faces. They are not stored in
                                        * array.
                                        */
      dealii::internal::DoFHandler::DoFFaces<dim> *faces;
-     
 +    dealii::internal::DoFHandler::DoFFaces<dim>* mg_faces;
-     
++
 +    std::vector<MGVertexDoFs> mg_vertex_dofs;
-     
++
 +    std::vector<types::global_dof_index> mg_used_dofs;
 +
                                       /**
                                        * Make accessor objects friends.
                                        */
  
  #ifndef DOXYGEN
  
 -template <> unsigned int DoFHandler<1>::n_boundary_dofs () const;
 -template <> unsigned int DoFHandler<1>::n_boundary_dofs (const FunctionMap &) const;
 -template <> unsigned int DoFHandler<1>::n_boundary_dofs (const std::set<types::boundary_id> &) const;
 +template <> types::global_dof_index DoFHandler<1>::n_boundary_dofs () const;
 +template <> types::global_dof_index DoFHandler<1>::n_boundary_dofs (const FunctionMap &) const;
- template <> types::global_dof_index DoFHandler<1>::n_boundary_dofs (const std::set<types::boundary_id_t> &) const;
++template <> types::global_dof_index DoFHandler<1>::n_boundary_dofs (const std::set<types::boundary_id> &) const;
  
  
  /* ----------------------- Inline functions ---------------------------------- */
index fd7eda09b9a9a34e7597775b835ae01456e4c2fa,22a9d90963ddc6e80939107d7b8dfdf1847dc8d7..b6f7b2a5311cb65ced485353adcf87fed79367b5
@@@ -493,10 -493,10 +493,10 @@@ namespace DoFTool
    template <class DH, class SparsityPattern>
    void
    make_sparsity_pattern (const DH               &dof,
 -                         SparsityPattern        &sparsity_pattern,
 -                         const ConstraintMatrix &constraints = ConstraintMatrix(),
 -                         const bool              keep_constrained_dofs = true,
 -                         const types::subdomain_id subdomain_id = types::invalid_subdomain_id);
 +                                   SparsityPattern        &sparsity_pattern,
 +                                   const ConstraintMatrix &constraints = ConstraintMatrix(),
 +                                   const bool              keep_constrained_dofs = true,
-                                    const types::subdomain_id_t subdomain_id = types::invalid_subdomain_id);
++                                   const types::subdomain_id subdomain_id = types::invalid_subdomain_id);
  
                                     /**
                                      * Locate non-zero entries for
    template <class DH, class SparsityPattern>
    void
    make_sparsity_pattern (const DH                 &dof,
 -                         const Table<2, Coupling> &coupling,
 -                         SparsityPattern          &sparsity_pattern,
 -                         const ConstraintMatrix   &constraints = ConstraintMatrix(),
 -                         const bool                keep_constrained_dofs = true,
 -                         const types::subdomain_id subdomain_id = types::invalid_subdomain_id);
 +                                   const Table<2, Coupling> &coupling,
 +                                   SparsityPattern          &sparsity_pattern,
 +                                   const ConstraintMatrix   &constraints = ConstraintMatrix(),
 +                                   const bool                keep_constrained_dofs = true,
-                                    const types::subdomain_id_t subdomain_id = types::invalid_subdomain_id);
++                                   const types::subdomain_id subdomain_id = types::invalid_subdomain_id);
  
                                     /**
                                      * @deprecated This is the old
    template<class DH, class SparsityPattern>
    void
    make_flux_sparsity_pattern (const DH        &dof_handler,
 -                              SparsityPattern &sparsity_pattern,
 -                              const ConstraintMatrix   &constraints,
 -                              const bool                keep_constrained_dofs = true,
 -                              const types::subdomain_id  subdomain_id = numbers::invalid_unsigned_int);
 +                                        SparsityPattern &sparsity_pattern,
 +                                        const ConstraintMatrix   &constraints,
 +                                        const bool                keep_constrained_dofs = true,
-                                         const types::subdomain_id_t  subdomain_id = numbers::invalid_unsigned_int);
++                                        const types::subdomain_id  subdomain_id = numbers::invalid_unsigned_int);
  
                                     /**
                                      * This function does the same as
    extract_boundary_dofs (const DH                   &dof_handler,
                           const std::vector<bool>    &component_mask,
                           std::vector<bool>          &selected_dofs,
-                          const std::set<types::boundary_id_t> &boundary_indicators = std::set<types::boundary_id_t>());
+                          const std::set<types::boundary_id> &boundary_indicators = std::set<types::boundary_id>());
  
 +  /**
 +   * This function does the same as the previous one but it
 +   * returns its result as an IndexSet rather than a std::vector@<bool@>.
 +   * Thus, it can also be called for DoFHandler objects that are
 +   * defined on parallel::distributed::Triangulation objects.
 +   *
 +   * @note If the DoFHandler object is indeed defined on a
 +   * parallel::distributed::Triangulation, then the @p selected_dofs
 +   * index set will contain only those degrees of freedom on the
 +   * boundary that belong to the locally relevant set (see
 +   * @ref GlossLocallyRelevantDof "locally relevant DoFs").
 +   *
 +   * @param dof_handler The object that describes which degrees of freedom
 +   *          live on which cell
 +   * @param component_mask A mask denoting the vector components of the
 +   *          finite element that should be considered (see also
 +   *          @ref GlossComponentMask).
 +   * @param selected_dofs The IndexSet object that is returned and that
 +   *          will contain the indices of degrees of freedom that are
 +   *          located on the boundary (and correspond to the selected
 +   *          vector components and boundary indicators, depending on
 +   *          the values of the @p component_mask and @p boundary_indicators
 +   *          arguments).
 +   * @param boundary_indicators If empty, this function extracts the
 +   *          indices of the degrees of freedom for all parts of the boundary.
 +   *          If it is a non-empty list, then the function only considers
 +   *          boundary faces with the boundary indicators listed in this
 +   *          argument.
 +   */
 +  template <class DH>
 +  void
 +  extract_boundary_dofs (const DH                   &dof_handler,
 +                         const std::vector<bool>    &component_mask,
 +                         IndexSet                   &selected_dofs,
 +                         const std::set<types::boundary_id_t> &boundary_indicators = std::set<types::boundary_id_t>());
 +
    /**
     * This function does the same as the previous one but it
     * returns its result as an IndexSet rather than a std::vector@<bool@>.
    template <class DH>
    void
    extract_dofs_with_support_on_boundary (const DH                   &dof_handler,
 -                                         const std::vector<bool>    &component_mask,
 -                                         std::vector<bool>          &selected_dofs,
 -                                         const std::set<types::boundary_id> &boundary_indicators = std::set<types::boundary_id>());
 +                                                           const std::vector<bool>    &component_mask,
 +                                                           std::vector<bool>          &selected_dofs,
-                                                            const std::set<types::boundary_id_t> &boundary_indicators = std::set<types::boundary_id_t>());
++                                                           const std::set<types::boundary_id> &boundary_indicators = std::set<types::boundary_id>());
  
                                     /**
                                      * @name Hanging Nodes
    template <class DH>
    void
    extract_subdomain_dofs (const DH           &dof_handler,
-                                     const types::subdomain_id_t subdomain_id,
 -                          const types::subdomain_id subdomain_id,
 -                          std::vector<bool>  &selected_dofs);
++                                    const types::subdomain_id subdomain_id,
 +                                    std::vector<bool>  &selected_dofs);
  
  
                                     /**
    template <class DH>
    void
    get_subdomain_association (const DH                  &dof_handler,
-                                        std::vector<types::subdomain_id_t> &subdomain);
 -                             std::vector<types::subdomain_id> &subdomain);
++                                       std::vector<types::subdomain_id> &subdomain);
  
                                     /**
                                      * Count how many degrees of freedom are
    template <class DH>
    unsigned int
    count_dofs_with_subdomain_association (const DH           &dof_handler,
-                                                            const types::subdomain_id_t subdomain);
 -                                         const types::subdomain_id subdomain);
++                                                           const types::subdomain_id subdomain);
  
                                     /**
                                      * Count how many degrees of freedom are
    template <class DH>
    void
    count_dofs_with_subdomain_association (const DH           &dof_handler,
-                                                            const types::subdomain_id_t subdomain,
 -                                         const types::subdomain_id subdomain,
 -                                         std::vector<unsigned int> &n_dofs_on_subdomain);
++                                                           const types::subdomain_id subdomain,
 +                                                           std::vector<unsigned int> &n_dofs_on_subdomain);
  
                                     /**
                                      * Return a set of indices that denotes
    template <class DH>
    IndexSet
    dof_indices_with_subdomain_association (const DH           &dof_handler,
-                                                             const types::subdomain_id_t subdomain);
 -                                          const types::subdomain_id subdomain);
 -
++                                                            const types::subdomain_id subdomain);
                                     // @}
                                     /**
                                      * @name Dof indices for patches
    template <class DH>
    void
    map_dof_to_boundary_indices (const DH                      &dof_handler,
-                                          const std::set<types::boundary_id_t> &boundary_indicators,
 -                               const std::set<types::boundary_id> &boundary_indicators,
 -                               std::vector<unsigned int>     &mapping);
++                                         const std::set<types::boundary_id> &boundary_indicators,
 +                                         std::vector<types::global_dof_index>     &mapping);
  
                                     /**
                                      * Return a list of support
index 3876c81c8d732430513bc1d4718c923df3807bc8,c15ffc2a426757e213463c92ec5f89262453b97b..81e9ee98649bec5bfdf247faa7e6fcea861805f3
@@@ -863,8 -862,8 +863,8 @@@ namespace h
                                          * indicators under
                                          * consideration.
                                          */
 -      unsigned int
 +      types::global_dof_index
-       n_boundary_dofs (const std::set<types::boundary_id_t> &boundary_indicators) const;
+       n_boundary_dofs (const std::set<types::boundary_id> &boundary_indicators) const;
  
                                         /**
                                          * Return the number of
                                          * copy such an object.
                                          */
        DoFHandler & operator = (const DoFHandler &);
-       
 +      class MGVertexDoFs {
 +        public:
 +          MGVertexDoFs ();
 +          ~MGVertexDoFs ();
 +          types::global_dof_index get_index (const unsigned int level, const unsigned int dof_number) const;
 +          void set_index (const unsigned int level, const unsigned int dof_number, const types::global_dof_index index);
 +      };
 +
                                         /**
                                          * Free all used memory.
                                          */
        void clear_space ();
-       
 +      template<int structdim>
 +      types::global_dof_index get_dof_index (const unsigned int obj_level, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index) const;
-       
++
 +      template<int structdim>
 +      void set_dof_index (const unsigned int obj_level, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index) const;
 +
                                         /**
                                          *  Create default tables for
                                          *  the active_fe_indices in
                                          * the present class.
                                          */
        std::vector<unsigned int>      vertex_dofs_offsets;
-       
 +      std::vector<MGVertexDoFs> mg_vertex_dofs;
 +
                                         /**
                                          * Array to store the
                                          * information, if a cell on
      return *tria;
    }
  
-   types::global_dof_index DoFHandler<dim, spacedim>::MGVertexDoFs::get_index (const unsigned int level, const unsigned int dof_number) const {
 +  template<int dim, int spacedim>
 +  inline
-   
++  DoFHandler<dim, spacedim>::MGVertexDoFs::MGVertexDoFs()
++  {
++    Assert (false, ExcNotImplemented ());
++  }
++
++  template<int dim, int spacedim>
++  inline
++  DoFHandler<dim, spacedim>::MGVertexDoFs::~MGVertexDoFs()
++  {
++    Assert (false, ExcNotImplemented ());
++  }
++
++  template<int dim, int spacedim>
++  inline
++  types::global_dof_index DoFHandler<dim, spacedim>::MGVertexDoFs::get_index (const unsigned int,
++                                                                            const unsigned int) const
++  {
 +    Assert (false, ExcNotImplemented ());
 +    return invalid_dof_index;
 +  }
-   void DoFHandler<dim, spacedim>::MGVertexDoFs::set_index (const unsigned int level, const unsigned int dof_number, const types::global_dof_index index) {
 +  template<int dim, int spacedim>
 +  inline
++  void DoFHandler<dim, spacedim>::MGVertexDoFs::set_index (const unsigned int,
++                                                         const unsigned int,
++                                                         types::global_dof_index)
++  {
 +    Assert (false, ExcNotImplemented ());
 +  }
  
  
  #endif
index bcd44fdf6e3b621e080560e3177da0b703e61dd6,94b82212dc953478024c25355a4acaaeb5f39ff1..7e51e549cd52635fa46eff839e2f7a7dd00df471
@@@ -356,16 -215,11 +356,11 @@@ MGTransferPrebuilt<VECTOR>::memory_cons
  {
    std::size_t result = sizeof(*this);
    result += sizeof(unsigned int) * sizes.size();
- #ifdef DEAL_PREFER_MATRIX_EZ
-   std::vector<std_cxx1x::shared_ptr<SparseMatrixEZ<double> > >::const_iterator m;
-   const std::vector<std_cxx1x::shared_ptr<SparseMatrixEZ<double> > >::const_iterator end = prolongation_matrices.end();
-   for (m = prolongation_matrices.begin(); m != end ; ++m)
-     result += *m->memory_consumption();
- #else
    for (unsigned int i=0;i<prolongation_matrices.size();++i)
      result += prolongation_matrices[i]->memory_consumption()
 -              + prolongation_sparsities[i]->memory_consumption();
 +                + prolongation_sparsities[i]->memory_consumption();
- #endif
    return result;
  }
  
index 0b962ce8039a3a04ddca85177ebf4db9d7a3814f,f6a1387560ce0487deb6e9a6a955923f61beed3e..0c85f89126c7e146e02053b8f182b1dc78520fb0
@@@ -767,9 -767,9 +767,9 @@@ namespace VectorTool
    void
    interpolate_boundary_values (const Mapping<DH::dimension,DH::space_dimension>            &mapping,
                                 const DH                 &dof,
-                                const types::boundary_id_t            boundary_component,
+                                const types::boundary_id            boundary_component,
                                 const Function<DH::space_dimension>           &boundary_function,
 -                               std::map<unsigned int,double> &boundary_values,
 +                               std::map<types::global_dof_index,double> &boundary_values,
                                 const std::vector<bool>       &component_mask = std::vector<bool>());
  
                                     /**
    template <class DH>
    void
    interpolate_boundary_values (const DH        &dof,
-                                const types::boundary_id_t            boundary_component,
+                                const types::boundary_id            boundary_component,
                                 const Function<DH::space_dimension>           &boundary_function,
 -                               std::map<unsigned int,double> &boundary_values,
 +                               std::map<types::global_dof_index,double> &boundary_values,
                                 const std::vector<bool>       &component_mask = std::vector<bool>());
  
  
index 09c182635e25bc6b3cf7b7b17133b3b3b3400bec,2df4c30f05ff8302998d7219d82f3be636275f1f..61a33e0d60c5b3900095725433c55d6fdf785de7
@@@ -1941,9 -1942,9 +1942,9 @@@ namespace VectorTool
    void
    interpolate_boundary_values (const Mapping<DH::dimension, DH::space_dimension>            &mapping,
                                 const DH                 &dof,
-                                const types::boundary_id_t            boundary_component,
+                                const types::boundary_id            boundary_component,
                                 const Function<DH::space_dimension>           &boundary_function,
 -                               std::map<unsigned int,double> &boundary_values,
 +                               std::map<types::global_dof_index,double> &boundary_values,
                                 const std::vector<bool>       &component_mask)
    {
      typename FunctionMap<DH::space_dimension>::type function_map;
    template <class DH>
    void
    interpolate_boundary_values (const DH                 &dof,
-                                const types::boundary_id_t            boundary_component,
+                                const types::boundary_id            boundary_component,
                                 const Function<DH::space_dimension>           &boundary_function,
 -                               std::map<unsigned int,double> &boundary_values,
 +                               std::map<types::global_dof_index,double> &boundary_values,
                                 const std::vector<bool>       &component_mask)
    {
      Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
      else
        AssertDimension (dof.get_fe().n_components(), component_mapping.size());
  
 -    std::vector<unsigned int> dof_to_boundary_mapping;
 +    std::vector<types::global_dof_index> dof_to_boundary_mapping;
-     std::set<types::boundary_id_t> selected_boundary_components;
+     std::set<types::boundary_id> selected_boundary_components;
      for (typename FunctionMap<spacedim>::type::const_iterator i=boundary_functions.begin();
           i!=boundary_functions.end(); ++i)
        selected_boundary_components.insert (i->first);
Simple merge
index dc84a2fbd25f4039457a2d30278c4d1a803af098,9e466cd2d0bb311e284a669b5f3f255daf9a91d7..4be24de818feb783d705743f0668a9bda7a65b8b
@@@ -257,427 -257,99 +257,427 @@@ namespace interna
                                            * was selected when calling
                                            * @p distribute_dofs the last time.
                                            */
 -        template <int spacedim>
 -        static
 -        void reserve_space (DoFHandler<1,spacedim> &dof_handler)
 -          {
 -            dof_handler.vertex_dofs
 -              .resize(dof_handler.tria->n_vertices() *
 -                      dof_handler.selected_fe->dofs_per_vertex,
 -                      DoFHandler<1,spacedim>::invalid_dof_index);
 -
 -            for (unsigned int i=0; i<dof_handler.tria->n_levels(); ++i)
 -              {
 -                dof_handler.levels
 -                  .push_back (new internal::DoFHandler::DoFLevel<1>);
 -
 -                dof_handler.levels.back()->lines.dofs
 -                  .resize (dof_handler.tria->n_raw_lines(i) *
 -                           dof_handler.selected_fe->dofs_per_line,
 -                           DoFHandler<1,spacedim>::invalid_dof_index);
 -
 -                dof_handler.levels.back()->cell_dof_indices_cache
 -                  .resize (dof_handler.tria->n_raw_lines(i) *
 -                           dof_handler.selected_fe->dofs_per_cell,
 -                           DoFHandler<1,spacedim>::invalid_dof_index);
 -              }
 -          }
 -
 -
 -        template <int spacedim>
 -        static
 -        void reserve_space (DoFHandler<2,spacedim> &dof_handler)
 -          {
 -            dof_handler.vertex_dofs
 -              .resize(dof_handler.tria->n_vertices() *
 -                      dof_handler.selected_fe->dofs_per_vertex,
 -                      DoFHandler<2,spacedim>::invalid_dof_index);
 -
 -            for (unsigned int i=0; i<dof_handler.tria->n_levels(); ++i)
 -              {
 -                dof_handler.levels.push_back (new internal::DoFHandler::DoFLevel<2>);
 -
 -                dof_handler.levels.back()->quads.dofs
 -                  .resize (dof_handler.tria->n_raw_quads(i) *
 -                           dof_handler.selected_fe->dofs_per_quad,
 -                           DoFHandler<2,spacedim>::invalid_dof_index);
 -
 -                dof_handler.levels.back()->cell_dof_indices_cache
 -                  .resize (dof_handler.tria->n_raw_quads(i) *
 -                           dof_handler.selected_fe->dofs_per_cell,
 -                           DoFHandler<2,spacedim>::invalid_dof_index);
 -              }
 -
 -            dof_handler.faces = new internal::DoFHandler::DoFFaces<2>;
 -            dof_handler.faces->lines.dofs
 -              .resize (dof_handler.tria->n_raw_lines() *
 -                       dof_handler.selected_fe->dofs_per_line,
 -                       DoFHandler<2,spacedim>::invalid_dof_index);
 -          }
 -
 -
 -        template <int spacedim>
 -        static
 -        void reserve_space (DoFHandler<3,spacedim> &dof_handler)
 -          {
 -            dof_handler.vertex_dofs
 -              .resize(dof_handler.tria->n_vertices() *
 -                      dof_handler.selected_fe->dofs_per_vertex,
 -                      DoFHandler<3,spacedim>::invalid_dof_index);
 -
 -            for (unsigned int i=0; i<dof_handler.tria->n_levels(); ++i)
 -              {
 -                dof_handler.levels.push_back (new internal::DoFHandler::DoFLevel<3>);
 -
 -                dof_handler.levels.back()->hexes.dofs
 -                  .resize (dof_handler.tria->n_raw_hexs(i) *
 -                           dof_handler.selected_fe->dofs_per_hex,
 -                           DoFHandler<3,spacedim>::invalid_dof_index);
 -
 -                dof_handler.levels.back()->cell_dof_indices_cache
 -                  .resize (dof_handler.tria->n_raw_hexs(i) *
 -                           dof_handler.selected_fe->dofs_per_cell,
 -                           DoFHandler<3,spacedim>::invalid_dof_index);
 -              }
 -            dof_handler.faces = new internal::DoFHandler::DoFFaces<3>;
 -
 -            dof_handler.faces->lines.dofs
 -              .resize (dof_handler.tria->n_raw_lines() *
 -                       dof_handler.selected_fe->dofs_per_line,
 -                       DoFHandler<3,spacedim>::invalid_dof_index);
 -            dof_handler.faces->quads.dofs
 -              .resize (dof_handler.tria->n_raw_quads() *
 -                       dof_handler.selected_fe->dofs_per_quad,
 -                       DoFHandler<3,spacedim>::invalid_dof_index);
 -          }
 +          template <int spacedim>
 +          static
 +          void reserve_space (DoFHandler<1,spacedim> &dof_handler)
 +            {
 +              dof_handler.vertex_dofs
 +                .resize(dof_handler.tria->n_vertices() *
 +                            dof_handler.selected_fe->dofs_per_vertex,
 +                            DoFHandler<1,spacedim>::invalid_dof_index);
 +
 +              for (unsigned int i=0; i<dof_handler.tria->n_levels(); ++i)
 +                {
 +                      dof_handler.levels
 +                        .push_back (new internal::DoFHandler::DoFLevel<1>);
 +
 +                      dof_handler.levels.back()->lines.dofs
 +                        .resize (dof_handler.tria->n_raw_lines(i) *
 +                                     dof_handler.selected_fe->dofs_per_line,
 +                                     DoFHandler<1,spacedim>::invalid_dof_index);
 +
 +                      dof_handler.levels.back()->cell_dof_indices_cache
 +                        .resize (dof_handler.tria->n_raw_lines(i) *
 +                                     dof_handler.selected_fe->dofs_per_cell,
 +                                     DoFHandler<1,spacedim>::invalid_dof_index);
 +                }
 +            }
 +
 +
 +          template <int spacedim>
 +          static
 +          void reserve_space (DoFHandler<2,spacedim> &dof_handler)
 +            {
 +              dof_handler.vertex_dofs
 +                .resize(dof_handler.tria->n_vertices() *
 +                            dof_handler.selected_fe->dofs_per_vertex,
 +                            DoFHandler<2,spacedim>::invalid_dof_index);
 +
 +              for (unsigned int i=0; i<dof_handler.tria->n_levels(); ++i)
 +                {
 +                      dof_handler.levels.push_back (new internal::DoFHandler::DoFLevel<2>);
 +
 +                      dof_handler.levels.back()->quads.dofs
 +                        .resize (dof_handler.tria->n_raw_quads(i) *
 +                                     dof_handler.selected_fe->dofs_per_quad,
 +                                     DoFHandler<2,spacedim>::invalid_dof_index);
 +
 +                      dof_handler.levels.back()->cell_dof_indices_cache
 +                        .resize (dof_handler.tria->n_raw_quads(i) *
 +                                     dof_handler.selected_fe->dofs_per_cell,
 +                                     DoFHandler<2,spacedim>::invalid_dof_index);
 +                }
 +
 +              dof_handler.faces = new internal::DoFHandler::DoFFaces<2>;
 +              dof_handler.faces->lines.dofs
 +                .resize (dof_handler.tria->n_raw_lines() *
 +                             dof_handler.selected_fe->dofs_per_line,
 +                             DoFHandler<2,spacedim>::invalid_dof_index);
 +            }
 +
 +
 +          template <int spacedim>
 +          static
 +          void reserve_space (DoFHandler<3,spacedim> &dof_handler)
 +            {
 +              dof_handler.vertex_dofs
 +                .resize(dof_handler.tria->n_vertices() *
 +                            dof_handler.selected_fe->dofs_per_vertex,
 +                            DoFHandler<3,spacedim>::invalid_dof_index);
 +
 +              for (unsigned int i=0; i<dof_handler.tria->n_levels(); ++i)
 +                {
 +                      dof_handler.levels.push_back (new internal::DoFHandler::DoFLevel<3>);
 +
 +                      dof_handler.levels.back()->hexes.dofs
 +                        .resize (dof_handler.tria->n_raw_hexs(i) *
 +                                     dof_handler.selected_fe->dofs_per_hex,
 +                                     DoFHandler<3,spacedim>::invalid_dof_index);
 +
 +                      dof_handler.levels.back()->cell_dof_indices_cache
 +                        .resize (dof_handler.tria->n_raw_hexs(i) *
 +                                     dof_handler.selected_fe->dofs_per_cell,
 +                                     DoFHandler<3,spacedim>::invalid_dof_index);
 +                }
 +              dof_handler.faces = new internal::DoFHandler::DoFFaces<3>;
 +
 +              dof_handler.faces->lines.dofs
 +                .resize (dof_handler.tria->n_raw_lines() *
 +                             dof_handler.selected_fe->dofs_per_line,
 +                             DoFHandler<3,spacedim>::invalid_dof_index);
 +              dof_handler.faces->quads.dofs
 +                .resize (dof_handler.tria->n_raw_quads() *
 +                             dof_handler.selected_fe->dofs_per_quad,
 +                             DoFHandler<3,spacedim>::invalid_dof_index);
 +            }
-         
++
 +          template<int spacedim>
 +          static
 +          void reserve_space_mg (DoFHandler<1, spacedim>& dof_handler) {
 +            Assert (dof_handler.get_tria ().n_levels () > 0, ExcMessage ("Invalid triangulation"));
 +            dof_handler.clear_mg_space ();
-           
++
 +            const dealii::Triangulation<1, spacedim>& tria = dof_handler.get_tria ();
 +            const unsigned int& dofs_per_line = dof_handler.get_fe ().dofs_per_line;
 +          const unsigned int& n_levels = tria.n_levels ();
-           
++
 +            for (unsigned int i = 0; i < n_levels; ++i) {
 +              dof_handler.mg_levels.push_back (new internal::DoFHandler::DoFLevel<1>);
 +              dof_handler.mg_levels.back ()->lines.dofs = std::vector<unsigned int> (tria.n_raw_lines (i) * dofs_per_line, DoFHandler<1>::invalid_dof_index);
 +            }
-           
++
 +            const unsigned int& n_vertices = tria.n_vertices ();
-           
++
 +            dof_handler.mg_vertex_dofs.resize (n_vertices);
-           
++
 +            std::vector<unsigned int> max_level (n_vertices, 0);
 +            std::vector<unsigned int> min_level (n_vertices, n_levels);
-           
++
 +            for (typename dealii::Triangulation<1, spacedim>::cell_iterator cell = tria.begin (); cell != tria.end (); ++cell) {
 +              const unsigned int level = cell->level ();
-             
++
 +              for (unsigned int vertex = 0; vertex < GeometryInfo<1>::vertices_per_cell; ++vertex) {
 +                const unsigned int vertex_index = cell->vertex_index (vertex);
-               
++
 +                if (min_level[vertex_index] > level)
 +                  min_level[vertex_index] = level;
-               
++
 +                if (max_level[vertex_index] < level)
 +                  max_level[vertex_index] = level;
 +              }
 +            }
-           
++
 +            for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
 +              if (tria.vertex_used (vertex)) {
 +                Assert (min_level[vertex] < n_levels, ExcInternalError ());
 +                Assert (max_level[vertex] >= min_level[vertex], ExcInternalError ());
 +                dof_handler.mg_vertex_dofs[vertex].init (min_level[vertex], max_level[vertex], dof_handler.get_fe ().dofs_per_vertex);
 +              }
-             
++
 +              else {
 +                Assert (min_level[vertex] == n_levels, ExcInternalError ());
 +                Assert (max_level[vertex] == 0, ExcInternalError ());
 +                dof_handler.mg_vertex_dofs[vertex].init (1, 0, 0);
 +              }
 +          }
-         
++
 +          template<int spacedim>
 +          static
 +          void reserve_space_mg (DoFHandler<2, spacedim>& dof_handler) {
 +            Assert (dof_handler.get_tria ().n_levels () > 0, ExcMessage ("Invalid triangulation"));
 +            dof_handler.clear_mg_space ();
-           
++
 +            const dealii::FiniteElement<2, spacedim>& fe = dof_handler.get_fe ();
 +            const dealii::Triangulation<2, spacedim>& tria = dof_handler.get_tria ();
 +            const unsigned int& n_levels = tria.n_levels ();
-           
++
 +            for (unsigned int i = 0; i < n_levels; ++i) {
 +              dof_handler.mg_levels.push_back (new internal::DoFHandler::DoFLevel<2>);
 +              dof_handler.mg_levels.back ()->quads.dofs = std::vector<unsigned int> (tria.n_raw_quads (i) * fe.dofs_per_quad, DoFHandler<2>::invalid_dof_index);
 +            }
-           
++
 +            dof_handler.mg_faces = new internal::DoFHandler::DoFFaces<2>;
 +            dof_handler.mg_faces->lines.dofs = std::vector<unsigned int> (tria.n_raw_lines () * fe.dofs_per_line, DoFHandler<2>::invalid_dof_index);
-           
++
 +            const unsigned int& n_vertices = tria.n_vertices ();
-           
++
 +            dof_handler.mg_vertex_dofs.resize (n_vertices);
-           
++
 +            std::vector<unsigned int> max_level (n_vertices, 0);
 +            std::vector<unsigned int> min_level (n_vertices, n_levels);
-           
++
 +            for (typename dealii::Triangulation<2, spacedim>::cell_iterator cell = tria.begin (); cell != tria.end (); ++cell) {
 +              const unsigned int level = cell->level ();
-             
++
 +              for (unsigned int vertex = 0; vertex < GeometryInfo<2>::vertices_per_cell; ++vertex) {
 +                const unsigned int vertex_index = cell->vertex_index (vertex);
-               
++
 +                if (min_level[vertex_index] > level)
 +                  min_level[vertex_index] = level;
-               
++
 +                if (max_level[vertex_index] < level)
 +                  max_level[vertex_index] = level;
 +              }
 +            }
-           
++
 +            for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
 +              if (tria.vertex_used (vertex)) {
 +                Assert (min_level[vertex] < n_levels, ExcInternalError ());
 +                Assert (max_level[vertex] >= min_level[vertex], ExcInternalError ());
 +                dof_handler.mg_vertex_dofs[vertex].init (min_level[vertex], max_level[vertex], fe.dofs_per_vertex);
 +              }
-             
++
 +              else {
 +                Assert (min_level[vertex] == n_levels, ExcInternalError ());
 +                Assert (max_level[vertex] == 0, ExcInternalError ());
 +                dof_handler.mg_vertex_dofs[vertex].init (1, 0, 0);
 +              }
 +          }
-         
++
 +          template<int spacedim>
 +          static
 +          void reserve_space_mg (DoFHandler<3, spacedim>& dof_handler) {
 +            Assert (dof_handler.get_tria ().n_levels () > 0, ExcMessage ("Invalid triangulation"));
 +            dof_handler.clear_mg_space ();
-           
++
 +            const dealii::FiniteElement<3, spacedim>& fe = dof_handler.get_fe ();
 +            const dealii::Triangulation<3, spacedim>& tria = dof_handler.get_tria ();
 +            const unsigned int& n_levels = tria.n_levels ();
-           
++
 +            for (unsigned int i = 0; i < n_levels; ++i) {
 +              dof_handler.mg_levels.push_back (new internal::DoFHandler::DoFLevel<3>);
 +              dof_handler.mg_levels.back ()->hexes.dofs = std::vector<unsigned int> (tria.n_raw_hexs (i) * fe.dofs_per_hex, DoFHandler<3>::invalid_dof_index);
 +            }
-           
++
 +            dof_handler.mg_faces = new internal::DoFHandler::DoFFaces<3>;
 +            dof_handler.mg_faces->lines.dofs = std::vector<unsigned int> (tria.n_raw_lines () * fe.dofs_per_line, DoFHandler<3>::invalid_dof_index);
 +            dof_handler.mg_faces->quads.dofs = std::vector<unsigned int> (tria.n_raw_quads () * fe.dofs_per_quad, DoFHandler<3>::invalid_dof_index);
-           
++
 +            const unsigned int& n_vertices = tria.n_vertices ();
-           
++
 +            dof_handler.mg_vertex_dofs.resize (n_vertices);
-           
++
 +            std::vector<unsigned int> max_level (n_vertices, 0);
 +            std::vector<unsigned int> min_level (n_vertices, n_levels);
-           
++
 +            for (typename dealii::Triangulation<3, spacedim>::cell_iterator cell = tria.begin (); cell != tria.end (); ++cell) {
 +              const unsigned int level = cell->level ();
-             
++
 +              for (unsigned int vertex = 0; vertex < GeometryInfo<3>::vertices_per_cell; ++vertex) {
 +                const unsigned int vertex_index = cell->vertex_index (vertex);
-               
++
 +                if (min_level[vertex_index] > level)
 +                  min_level[vertex_index] = level;
-               
++
 +                if (max_level[vertex_index] < level)
 +                  max_level[vertex_index] = level;
 +              }
 +            }
-           
++
 +            for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
 +              if (tria.vertex_used (vertex)) {
 +                Assert (min_level[vertex] < n_levels, ExcInternalError ());
 +                Assert (max_level[vertex] >= min_level[vertex], ExcInternalError ());
 +                dof_handler.mg_vertex_dofs[vertex].init (min_level[vertex], max_level[vertex], fe.dofs_per_vertex);
 +              }
-             
++
 +              else {
 +                Assert (min_level[vertex] == n_levels, ExcInternalError ());
 +                Assert (max_level[vertex] == 0, ExcInternalError ());
 +                dof_handler.mg_vertex_dofs[vertex].init (1, 0, 0);
 +              }
 +          }
-         
++
 +          template<int spacedim>
 +          static
 +          types::global_dof_index distribute_dofs_on_cell (typename DoFHandler<1, spacedim>::cell_iterator& cell, types::global_dof_index next_free_dof) {
 +                const FiniteElement<1, spacedim>& fe = cell->get_fe ();
-               
++
 +            if (fe.dofs_per_vertex > 0)
 +              for (unsigned int vertex = 0; vertex < GeometryInfo<1>::vertices_per_cell; ++vertex) {
 +                typename DoFHandler<1, spacedim>::cell_iterator neighbor = cell->neighbor (vertex);
-               
++
 +                if (neighbor.state () == IteratorState::valid)
 +                  if (neighbor->user_flag_set () && (neighbor->level () == cell->level ())) {
 +                    if (vertex == 0)
 +                      for (unsigned int dof = 0; dof < fe.dofs_per_vertex; ++dof)
 +                        cell->set_mg_vertex_dof_index (cell->level (), 0, dof, neighbor->mg_vertex_dof_index (cell->level (), 1, dof));
-                   
++
 +                    else
 +                      for (unsigned int dof = 0; dof < fe.dofs_per_vertex; ++dof)
 +                        cell->set_mg_vertex_dof_index (cell->level (), 1, dof, neighbor->mg_vertex_dof_index (cell->level (), 0, dof));
-                   
++
 +                    continue;
 +                  }
-               
++
 +                for (unsigned int dof = 0; dof < fe.dofs_per_vertex; ++dof)
 +                  cell->set_mg_vertex_dof_index (cell->level (), vertex, dof, next_free_dof++);
 +              }
-           
++
 +            if (fe.dofs_per_line > 0)
 +              for (unsigned int dof = 0; dof < fe.dofs_per_line; ++dof)
 +                cell->set_mg_dof_index (cell->level (), dof, next_free_dof++);
-           
++
 +            cell->set_user_flag ();
 +            return next_free_dof;
 +          }
-         
++
 +          template<int spacedim>
 +          static
 +          types::global_dof_index distribute_dofs_on_cell (typename DoFHandler<2, spacedim>::cell_iterator& cell, types::global_dof_index next_free_dof) {
 +            const FiniteElement<2, spacedim>& fe = cell->get_fe ();
-           
++
 +            if (fe.dofs_per_vertex > 0)
 +              for (unsigned int vertex = 0; vertex < GeometryInfo<2>::vertices_per_cell; ++vertex)
 +                if (cell->mg_vertex_dof_index (cell->level (), vertex, 0) == DoFHandler<2>::invalid_dof_index)
 +                  for (unsigned int dof = 0; dof < fe.dofs_per_vertex; ++dof)
 +                    cell->set_mg_vertex_dof_index (cell->level (), vertex, dof, next_free_dof++);
-           
++
 +            if (fe.dofs_per_line > 0)
 +              for (unsigned int face = 0; face < GeometryInfo<2>::faces_per_cell; ++face) {
 +                typename DoFHandler<2, spacedim>::line_iterator line = cell->line (face);
-               
++
 +                if (line->mg_dof_index (cell->level (), 0) == DoFHandler<2>::invalid_dof_index)
 +                  for (unsigned int dof = 0; dof < fe.dofs_per_line; ++dof)
 +                    line->set_mg_dof_index (cell->level (), dof, next_free_dof++);
 +              }
-           
++
 +            if (fe.dofs_per_quad > 0)
 +              for (unsigned int dof = 0; dof < fe.dofs_per_quad; ++dof)
 +                cell->set_mg_dof_index (cell->level (), dof, next_free_dof++);
-           
++
 +            cell->set_user_flag ();
 +            return next_free_dof;
 +          }
-         
++
 +          template<int spacedim>
 +          static
 +          types::global_dof_index distribute_dofs_on_cell (typename DoFHandler<3, spacedim>::cell_iterator& cell, types::global_dof_index next_free_dof) {
 +            const FiniteElement<3, spacedim>& fe = cell->get_fe ();
-           
++
 +            if (fe.dofs_per_vertex > 0)
 +              for (unsigned int vertex = 0; vertex < GeometryInfo<3>::vertices_per_cell; ++vertex)
 +                if (cell->mg_vertex_dof_index (cell->level (), vertex, 0) == DoFHandler<3>::invalid_dof_index)
 +                  for (unsigned int dof = 0; dof < fe.dofs_per_vertex; ++dof)
 +                    cell->set_mg_vertex_dof_index (cell->level (), vertex, dof, next_free_dof++);
-           
++
 +            if (fe.dofs_per_line > 0)
 +              for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_cell; ++line) {
 +                typename DoFHandler<3, spacedim>::line_iterator line_it = cell->line (line);
-               
++
 +                if (line_it->mg_dof_index (cell->level (), 0) == DoFHandler<3>::invalid_dof_index)
 +                  for (unsigned int dof = 0; dof < fe.dofs_per_line; ++dof)
 +                    line_it->set_mg_dof_index (cell->level (), dof, next_free_dof++);
 +              }
-           
++
 +            if (fe.dofs_per_quad > 0)
 +              for (unsigned int face = 0; face < GeometryInfo<3>::quads_per_cell; ++face) {
 +                typename DoFHandler<3, spacedim>::quad_iterator quad = cell->quad (face);
-               
++
 +                if (quad->mg_dof_index (cell->level (), 0) == DoFHandler<3>::invalid_dof_index)
 +                  for (unsigned int dof = 0; dof < fe.dofs_per_quad; ++dof)
 +                    quad->set_mg_dof_index (cell->level (), dof, next_free_dof++);
 +              }
-           
++
 +            if (fe.dofs_per_hex > 0)
 +              for (unsigned int dof = 0; dof < fe.dofs_per_hex; ++dof)
 +                cell->set_mg_dof_index (cell->level (), dof, next_free_dof++);
-           
++
 +            cell->set_user_flag ();
 +            return next_free_dof;
 +          }
-         
++
 +          template<int spacedim>
 +          static
 +          types::global_dof_index get_dof_index (const DoFHandler<1, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<1>& mg_level, internal::DoFHandler::DoFFaces<1>&, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const int2type<1>) {
 +            return mg_level.lines.get_dof_index (dof_handler, obj_index, fe_index, local_index);
 +          }
-         
++
 +          template<int spacedim>
 +          static
 +          types::global_dof_index get_dof_index (const DoFHandler<2, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<2>&, internal::DoFHandler::DoFFaces<2>& mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const int2type<1>) {
 +            return mg_faces.lines.get_dof_index (dof_handler, obj_index, fe_index, local_index);
 +          }
-         
++
 +          template<int spacedim>
 +          static
 +          types::global_dof_index get_dof_index (const DoFHandler<2, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<2>& mg_level, internal::DoFHandler::DoFFaces<2>&, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const int2type<2>) {
 +            return mg_level.quads.get_dof_index (dof_handler, obj_index, fe_index, local_index);
 +          }
-         
++
 +          template<int spacedim>
 +          static
 +          types::global_dof_index get_dof_index (const DoFHandler<3, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<3>&, internal::DoFHandler::DoFFaces<3>& mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const int2type<1>) {
 +            return mg_faces.lines.get_dof_index (dof_handler, obj_index, fe_index, local_index);
 +          }
-         
++
 +          template<int spacedim>
 +          static
 +          types::global_dof_index get_dof_index (const DoFHandler<3, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<3>&, internal::DoFHandler::DoFFaces<3>& mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const int2type<2>) {
 +            return mg_faces.quads.get_dof_index (dof_handler, obj_index, fe_index, local_index);
 +          }
-         
++
 +          template<int spacedim>
 +          static
 +          types::global_dof_index get_dof_index (const DoFHandler<3, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<3>& mg_level, internal::DoFHandler::DoFFaces<3>&, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const int2type<3>) {
 +            return mg_level.hexes.get_dof_index (dof_handler, obj_index, fe_index, local_index);
 +          }
-         
++
 +          template<int spacedim>
 +          static
 +          void set_dof_index (const DoFHandler<1, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<1>& mg_level, internal::DoFHandler::DoFFaces<1>&, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const int2type<1>) {
 +            mg_level.lines.set_dof_index (dof_handler, obj_index, fe_index, local_index, global_index);
 +          }
-         
++
 +          template<int spacedim>
 +          static
 +          void set_dof_index (const DoFHandler<2, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<2>&, internal::DoFHandler::DoFFaces<2>& mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const int2type<1>) {
 +            mg_faces.lines.set_dof_index (dof_handler, obj_index, fe_index, local_index, global_index);
 +          }
-         
++
 +          template<int spacedim>
 +          static
 +          void set_dof_index (const DoFHandler<2, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<2>& mg_level, internal::DoFHandler::DoFFaces<2>&, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const int2type<2>) {
 +            mg_level.quads.set_dof_index (dof_handler, obj_index, fe_index, local_index, global_index);
 +          }
-         
++
 +          template<int spacedim>
 +          static
 +          void set_dof_index (const DoFHandler<3, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<3>&, internal::DoFHandler::DoFFaces<3>& mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const int2type<1>) {
 +            mg_faces.lines.set_dof_index (dof_handler, obj_index, fe_index, local_index, global_index);
 +          }
-         
++
 +          template<int spacedim>
 +          static
 +          void set_dof_index (const DoFHandler<3, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<3>&, internal::DoFHandler::DoFFaces<3>& mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const int2type<2>) {
 +            mg_faces.quads.set_dof_index (dof_handler, obj_index, fe_index, local_index, global_index);
 +          }
-         
++
 +          template<int spacedim>
 +          static
 +          void set_dof_index (const DoFHandler<3, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<3>& mg_level, internal::DoFHandler::DoFFaces<3>&, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const int2type<3>) {
 +            mg_level.hexes.set_dof_index (dof_handler, obj_index, fe_index, local_index, global_index);
 +          }
      };
    }
  }
@@@ -1819,15 -1489,15 +1819,15 @@@ types::global_dof_index DoFHandler<1>::
  
  
  template <>
- types::global_dof_index DoFHandler<1>::n_boundary_dofs (const std::set<types::boundary_id_t> &boundary_indicators) const
 -unsigned int DoFHandler<1>::n_boundary_dofs (const std::set<types::boundary_id> &boundary_indicators) const
++types::global_dof_index DoFHandler<1>::n_boundary_dofs (const std::set<types::boundary_id> &boundary_indicators) const
  {
                                     // check that only boundary
                                     // indicators 0 and 1 are allowed
                                     // in 1d
-   for (std::set<types::boundary_id_t>::const_iterator i=boundary_indicators.begin();
+   for (std::set<types::boundary_id>::const_iterator i=boundary_indicators.begin();
         i!=boundary_indicators.end(); ++i)
      Assert ((*i == 0) || (*i == 1),
 -            ExcInvalidBoundaryIndicator());
 +              ExcInvalidBoundaryIndicator());
  
    return boundary_indicators.size()*get_fe().dofs_per_vertex;
  }
@@@ -1858,15 -1528,15 +1858,15 @@@ types::global_dof_index DoFHandler<1,2>
  
  
  template <>
- types::global_dof_index DoFHandler<1,2>::n_boundary_dofs (const std::set<types::boundary_id_t> &boundary_indicators) const
 -unsigned int DoFHandler<1,2>::n_boundary_dofs (const std::set<types::boundary_id> &boundary_indicators) const
++types::global_dof_index DoFHandler<1,2>::n_boundary_dofs (const std::set<types::boundary_id> &boundary_indicators) const
  {
                                     // check that only boundary
                                     // indicators 0 and 1 are allowed
                                     // in 1d
-   for (std::set<types::boundary_id_t>::const_iterator i=boundary_indicators.begin();
+   for (std::set<types::boundary_id>::const_iterator i=boundary_indicators.begin();
         i!=boundary_indicators.end(); ++i)
      Assert ((*i == 0) || (*i == 1),
 -            ExcInvalidBoundaryIndicator());
 +              ExcInvalidBoundaryIndicator());
  
    return boundary_indicators.size()*get_fe().dofs_per_vertex;
  }
@@@ -1932,11 -1602,11 +1932,11 @@@ DoFHandler<dim,spacedim>::n_boundary_do
  
  
  template<int dim, int spacedim>
 -unsigned int
 +types::global_dof_index
- DoFHandler<dim,spacedim>::n_boundary_dofs (const std::set<types::boundary_id_t> &boundary_indicators) const
+ DoFHandler<dim,spacedim>::n_boundary_dofs (const std::set<types::boundary_id> &boundary_indicators) const
  {
    Assert (boundary_indicators.find (types::internal_face_boundary_id) == boundary_indicators.end(),
 -          ExcInvalidBoundaryIndicator());
 +            ExcInvalidBoundaryIndicator());
  
    std::set<int> boundary_dofs;
  
@@@ -1964,25 -1634,16 +1964,25 @@@ std::size_
  DoFHandler<dim,spacedim>::memory_consumption () const
  {
    std::size_t mem = (MemoryConsumption::memory_consumption (tria) +
 -                     MemoryConsumption::memory_consumption (selected_fe) +
 -                     MemoryConsumption::memory_consumption (block_info_object) +
 -                     MemoryConsumption::memory_consumption (levels) +
 -                     MemoryConsumption::memory_consumption (*faces) +
 -                     MemoryConsumption::memory_consumption (faces) +
 -                     sizeof (number_cache) +
 -                     MemoryConsumption::memory_consumption (vertex_dofs));
 +                           MemoryConsumption::memory_consumption (selected_fe) +
 +                           MemoryConsumption::memory_consumption (block_info_object) +
 +                           MemoryConsumption::memory_consumption (levels) +
 +                           MemoryConsumption::memory_consumption (*faces) +
 +                           MemoryConsumption::memory_consumption (faces) +
 +                           sizeof (number_cache) +
 +                           MemoryConsumption::memory_consumption (vertex_dofs));
    for (unsigned int i=0; i<levels.size(); ++i)
      mem += MemoryConsumption::memory_consumption (*levels[i]);
-   
 +  for (unsigned int level = 0; level < mg_levels.size (); ++level)
 +    mem += mg_levels[level]->memory_consumption ();
-   
++
 +  mem += MemoryConsumption::memory_consumption (*mg_faces);
-   
++
 +  for (unsigned int i = 0; i < mg_vertex_dofs.size (); ++i)
 +    mem += sizeof (MGVertexDoFs) + (1 + mg_vertex_dofs[i].get_finest_level () - mg_vertex_dofs[i].get_coarsest_level ()) * sizeof (types::global_dof_index);
-   
++
 +  mem += MemoryConsumption::memory_consumption (mg_used_dofs);
    return mem;
  }
  
@@@ -2027,53 -1688,6 +2027,53 @@@ distribute_dofs (const FiniteElement<di
  }
  
  
-   
 +template<int dim, int spacedim>
 +void DoFHandler<dim, spacedim>::distribute_mg_dofs (const FiniteElement<dim, spacedim>& fe, const types::global_dof_index offset) {
 +  Assert ((dynamic_cast<const parallel::distributed::Triangulation<dim, spacedim>*> (&*tria) == 0), ExcMessage ("Invalid triangulation"));
 +  distribute_dofs (fe);
 +  reserve_space ();
-   
++
 +  const unsigned int& n_levels = (*tria).n_levels ();
-   
++
 +  mg_used_dofs.resize (n_levels, 0);
-   
++
 +  std::vector<bool> user_flags;
-   
++
 +  (*tria).save_user_flags (user_flags);
 +  const_cast<Triangulation<dim, spacedim>&>(*tria).clear_user_flags ();
-     
++
 +  for (unsigned int level = 0; level < n_levels; ++level) {
 +    types::global_dof_index next_free_dof = offset;
-     
++
 +    for (cell_iterator cell = begin (level); cell != end (level); ++cell)
 +      next_free_dof = internal::DoFHandler::Implementation::distribute_dofs_on_cell<spacedim> (cell, next_free_dof);
-   
++
 +    mg_used_dofs[level] = next_free_dof;
 +  }
-   
++
 +  const_cast<Triangulation<dim, spacedim>&>(*tria).load_user_flags (user_flags);
 +  block_info_object.initialize (*this, true, true);
 +}
 +
 +template<int dim, int spacedim>
 +void DoFHandler<dim, spacedim>::reserve_space () {
 +  internal::DoFHandler::Implementation::reserve_space_mg (*this);
 +}
 +
 +template<int dim, int spacedim>
 +void DoFHandler<dim, spacedim>::clear_mg_space () {
 +  for (unsigned int i = 0; i < mg_levels.size (); ++i)
 +    delete mg_levels[i];
-   
++
 +  mg_levels.clear ();
 +  delete mg_faces;
 +  mg_faces = NULL;
-   
++
 +  std::vector<MGVertexDoFs> tmp;
++
 +  std::swap (mg_vertex_dofs, tmp);
 +}
 +
  
  template<int dim, int spacedim>
  void DoFHandler<dim,spacedim>::initialize_local_block_info ()
@@@ -2196,72 -1809,6 +2196,72 @@@ void DoFHandler<dim,spacedim>::clear_sp
    std::swap (vertex_dofs, tmp);
  }
  
-   
 +template<int dim, int spacedim>
 +template<int structdim>
 +types::global_dof_index DoFHandler<dim, spacedim>::get_dof_index (const unsigned int obj_level, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index) const {
 +  return internal::DoFHandler::Implementation::get_dof_index (*this, *this->mg_levels[obj_level], *this->mg_faces, obj_index, fe_index, local_index, internal::int2type<structdim> ());
 +}
 +
 +template<int dim, int spacedim>
 +template<int structdim>
 +void DoFHandler<dim, spacedim>::set_dof_index (const unsigned int obj_level, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index) const {
 +  internal::DoFHandler::Implementation::set_dof_index (*this, *this->mg_levels[obj_level], *this->mg_faces, obj_index, fe_index, local_index, global_index, internal::int2type<structdim> ());
 +}
 +
 +template<int dim, int spacedim>
 +DoFHandler<dim, spacedim>::MGVertexDoFs::MGVertexDoFs (): coarsest_level (numbers::invalid_unsigned_int), finest_level (0), indices (0), indices_offset (0) {
 +}
 +
 +template<int dim, int spacedim>
 +DoFHandler<dim, spacedim>::MGVertexDoFs::~MGVertexDoFs () {
 +  delete[] indices;
 +  delete[] indices_offset;
 +}
 +
 +template<int dim, int spacedim>
 +void DoFHandler<dim, spacedim>::MGVertexDoFs::init (const unsigned int cl, const unsigned int fl, const unsigned int dofs_per_vertex) {
 +  if (indices != 0) {
 +    delete[] indices;
 +    indices = 0;
 +  }
-   
++
 +  if (indices_offset != 0) {
 +    delete[] indices_offset;
 +    indices_offset = 0;
 +  }
-   
++
 +  coarsest_level = cl;
 +  finest_level = fl;
-   
++
 +  if (cl > fl)
 +    return;
-   
++
 +  const unsigned int n_levels = finest_level - coarsest_level + 1;
 +  const unsigned int n_indices = n_levels * dofs_per_vertex;
-   
++
 +  indices = new unsigned int[n_indices];
 +  Assert (indices != 0, ExcNoMemory ());
-   
++
 +  for (unsigned int i = 0; i < n_indices; ++i)
 +    indices[i] = DoFHandler<dim, spacedim>::invalid_dof_index;
-   
++
 +  indices_offset = new unsigned int[n_levels];
 +  Assert (indices != 0, ExcNoMemory ());
++
 +  for (unsigned int i = 0; i < n_levels; ++i)
 +    indices_offset[i] = i * dofs_per_vertex;
 +}
 +
 +template<int dim, int spacedim>
 +unsigned int DoFHandler<dim, spacedim>::MGVertexDoFs::get_coarsest_level () const {
 +  return coarsest_level;
 +}
 +
 +template<int dim, int spacedim>
 +unsigned int DoFHandler<dim, spacedim>::MGVertexDoFs::get_finest_level () const {
 +  return finest_level;
 +}
 +
  
  /*-------------- Explicit Instantiations -------------------------------*/
  #include "dof_handler.inst"
index 4e12ba94c2e3138d3f26c4a8df8fedd465973e9b,5c8185b46ff8fc9f87f6a5c5a200e85e1573fb76..ea53851a2d1f6071992b4450ad71c2382f73709e
@@@ -223,9 -224,9 +223,9 @@@ namespace interna
                                              */
            template <int dim, int spacedim>
            static
 -          unsigned int
 -          distribute_dofs (const unsigned int        offset,
 +        types::global_dof_index 
 +          distribute_dofs (const types::global_dof_index offset,
-                            const types::subdomain_id_t subdomain_id,
+                            const types::subdomain_id subdomain_id,
                             DoFHandler<dim,spacedim> &dof_handler)
              {
                const dealii::Triangulation<dim,spacedim> & tria
                                             // mark all own cells, that miss some
                                             // dof_data and collect the neighbors
                                             // that are going to send stuff to us
-           std::set<dealii::types::subdomain_id_t> senders;
+           std::set<dealii::types::subdomain_id> senders;
            {
 -            std::vector<unsigned int> local_dof_indices;
 +            std::vector<types::global_dof_index> local_dof_indices;
              typename DoFHandler<dim,spacedim>::active_cell_iterator
                cell, endc = dof_handler.end();
  
index fad2b36a3af39ba64b24f94a1df1f8c2e4049fae,99bec3b82325ce10f5ca1b0845ff9161c0cc9b0e..f74aba9bd0b797973aef2923da70e405e813cf37
@@@ -1038,9 -1037,24 +1038,24 @@@ namespace DoFRenumberin
  
  
  
+   template <class DH>
+   void
+   sort_selected_dofs_back (DH&                      dof_handler,
+                            const std::vector<bool>& selected_dofs,
+                            const unsigned int       level)
+   {
+     std::vector<unsigned int> renumbering(dof_handler.n_dofs(level),
+                                           DH::invalid_dof_index);
+     compute_sort_selected_dofs_back(renumbering, dof_handler, selected_dofs, level);
+     dof_handler.renumber_dofs(level, renumbering);
+   }
    template <class DH>
    void
 -  compute_sort_selected_dofs_back (std::vector<unsigned int>& new_indices,
 +  compute_sort_selected_dofs_back (std::vector<types::global_dof_index>& new_indices,
                                     const DH&                  dof_handler,
                                     const std::vector<bool>&   selected_dofs)
    {
index 85c08f4f0e42a6bd4869fa55e641eb9664be0bda,e23dd691a8399107fb8ca80851567498f4c01fec..d42c9c7f9676840e9d896cf88550061cbfb30c9c
@@@ -57,10 -57,10 +57,10 @@@ namespace DoFTool
    template <class DH, class SparsityPattern>
    void
    make_sparsity_pattern (const DH               &dof,
 -                         SparsityPattern        &sparsity,
 -                         const ConstraintMatrix &constraints,
 -                         const bool              keep_constrained_dofs,
 -                         const types::subdomain_id subdomain_id)
 +                                   SparsityPattern        &sparsity,
 +                                   const ConstraintMatrix &constraints,
 +                                   const bool              keep_constrained_dofs,
-                                    const types::subdomain_id_t subdomain_id)
++                                   const types::subdomain_id subdomain_id)
    {
      const unsigned int n_dofs = dof.n_dofs();
  
    template <class DH, class SparsityPattern>
    void
    make_sparsity_pattern (const DH                &dof,
 -                         const Table<2,Coupling> &couplings,
 -                         SparsityPattern         &sparsity,
 -                         const ConstraintMatrix  &constraints,
 -                         const bool               keep_constrained_dofs,
 -                         const types::subdomain_id subdomain_id)
 +                                   const Table<2,Coupling> &couplings,
 +                                   SparsityPattern         &sparsity,
 +                                   const ConstraintMatrix  &constraints,
 +                                   const bool               keep_constrained_dofs,
-                                    const types::subdomain_id_t subdomain_id)
++                                   const types::subdomain_id subdomain_id)
    {
      const unsigned int n_dofs = dof.n_dofs();
  
    template <class DH, class SparsityPattern>
    void
    make_flux_sparsity_pattern (const DH        &dof,
 -                              SparsityPattern &sparsity,
 -                              const ConstraintMatrix &constraints,
 -                              const bool              keep_constrained_dofs,
 +                                        SparsityPattern &sparsity,
 +                                        const ConstraintMatrix &constraints,
 +                                        const bool              keep_constrained_dofs,
-                               const types::subdomain_id_t subdomain_id)
+                               const types::subdomain_id subdomain_id)
    {
      const unsigned int n_dofs = dof.n_dofs();
  
    template <class DH>
    void
    extract_dofs_with_support_on_boundary (const DH                      &dof_handler,
 -                                         const std::vector<bool>       &component_select,
 -                                         std::vector<bool>             &selected_dofs,
 -                                         const std::set<types::boundary_id> &boundary_indicators)
 +                                                           const std::vector<bool>       &component_select,
 +                                                           std::vector<bool>             &selected_dofs,
-                                                            const std::set<types::boundary_id_t> &boundary_indicators)
++                                                           const std::set<types::boundary_id> &boundary_indicators)
    {
      AssertDimension (component_select.size(), n_components(dof_handler));
      Assert (boundary_indicators.find (types::internal_face_boundary_id) == boundary_indicators.end(),
    template <class DH>
    void
    extract_subdomain_dofs (const DH                   &dof_handler,
-                                     const types::subdomain_id_t subdomain_id,
 -                          const types::subdomain_id subdomain_id,
 -                          std::vector<bool>          &selected_dofs)
++                                    const types::subdomain_id subdomain_id,
 +                                    std::vector<bool>          &selected_dofs)
    {
      Assert(selected_dofs.size() == dof_handler.n_dofs(),
 -           ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs()));
 +             ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs()));
  
                                       // preset all values by false
      std::fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false);
    template <class DH>
    void
    get_subdomain_association (const DH                  &dof_handler,
-                                        std::vector<types::subdomain_id_t> &subdomain_association)
 -                             std::vector<types::subdomain_id> &subdomain_association)
++                                       std::vector<types::subdomain_id> &subdomain_association)
    {
                                       // if the Triangulation is distributed, the
                                       // only thing we can usefully ask is for
        endc = dof_handler.end();
      for (; cell!=endc; ++cell)
        {
 -        Assert (cell->is_artificial() == false,
 -                ExcMessage ("You can't call this function for meshes that "
 -                            "have artificial cells."));
 +          Assert (cell->is_artificial() == false,
 +                      ExcMessage ("You can't call this function for meshes that "
 +                                      "have artificial cells."));
  
-           const types::subdomain_id_t subdomain_id = cell->subdomain_id();
 -        const types::subdomain_id subdomain_id = cell->subdomain_id();
 -        const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
 -        local_dof_indices.resize (dofs_per_cell);
 -        cell->get_dof_indices (local_dof_indices);
++          const types::subdomain_id subdomain_id = cell->subdomain_id();
 +          const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
 +          local_dof_indices.resize (dofs_per_cell);
 +          cell->get_dof_indices (local_dof_indices);
  
                                           // set subdomain ids. if dofs
                                           // already have their values
    template <class DH>
    unsigned int
    count_dofs_with_subdomain_association (const DH           &dof_handler,
-                                                            const types::subdomain_id_t subdomain)
 -                                         const types::subdomain_id subdomain)
++                                                           const types::subdomain_id subdomain)
    {
-     std::vector<types::subdomain_id_t> subdomain_association (dof_handler.n_dofs());
+     std::vector<types::subdomain_id> subdomain_association (dof_handler.n_dofs());
      get_subdomain_association (dof_handler, subdomain_association);
  
      return std::count (subdomain_association.begin(),
    template <class DH>
    IndexSet
    dof_indices_with_subdomain_association (const DH           &dof_handler,
-                                                             const types::subdomain_id_t subdomain)
 -                                          const types::subdomain_id subdomain)
++                                                            const types::subdomain_id subdomain)
    {
  
                                       // If we have a distributed::Triangulation only
    template <class DH>
    void
    count_dofs_with_subdomain_association (const DH           &dof_handler,
-                                                            const types::subdomain_id_t subdomain,
 -                                         const types::subdomain_id subdomain,
 -                                         std::vector<unsigned int> &n_dofs_on_subdomain)
++                                                           const types::subdomain_id subdomain,
 +                                                           std::vector<unsigned int> &n_dofs_on_subdomain)
    {
      Assert (n_dofs_on_subdomain.size() == dof_handler.get_fe().n_components(),
 -            ExcDimensionMismatch (n_dofs_on_subdomain.size(),
 -                                  dof_handler.get_fe().n_components()));
 +              ExcDimensionMismatch (n_dofs_on_subdomain.size(),
 +                                                dof_handler.get_fe().n_components()));
      std::fill (n_dofs_on_subdomain.begin(), n_dofs_on_subdomain.end(), 0);
  
                                       // in debug mode, make sure that there are
  #if defined(DEAL_II_USE_P4EST) && defined(DEAL_II_COMPILER_SUPPORTS_MPI)
      const unsigned int dim = DH::dimension;
      const unsigned int spacedim = DH::space_dimension;
-     
      if (const parallel::distributed::Triangulation<dim,spacedim> * tria
 -        = (dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>
 -           (&dof_handler.get_tria())))
 +          = (dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>
 +             (&dof_handler.get_tria())))
        {
 -        std::vector<unsigned int> local_dof_count = dofs_per_component;
 +          std::vector<unsigned int> local_dof_count = dofs_per_component;
  
 -        MPI_Allreduce ( &local_dof_count[0], &dofs_per_component[0], n_target_components,
 -                        MPI_UNSIGNED, MPI_SUM, tria->get_communicator());
 +          MPI_Allreduce ( &local_dof_count[0], &dofs_per_component[0], n_target_components,
 +                                  MPI_UNSIGNED, MPI_SUM, tria->get_communicator());
        }
  #endif
    }
  #endif
  
  
 -        return n_parameters_on_fine_grid;
 +          return n_parameters_on_fine_grid;
        }
-       
-       
      }
    }
  
    template <class DH>
    void map_dof_to_boundary_indices (
      const DH                      &dof_handler,
-     const std::set<types::boundary_id_t> &boundary_indicators,
+     const std::set<types::boundary_id> &boundary_indicators,
 -    std::vector<unsigned int>     &mapping)
 +    std::vector<types::global_dof_index>     &mapping)
    {
      Assert (&dof_handler.get_fe() != 0, ExcNoFESelected());
      Assert (boundary_indicators.find (types::internal_face_boundary_id) == boundary_indicators.end(),
  
                  local_dof_indices.resize(cell->get_fe().dofs_per_cell);
                  cell->get_dof_indices(local_dof_indices);
-  
                  const std::vector<Point<spacedim> > & points =
 -                    fe_values.get_quadrature_points();
 +                  fe_values.get_quadrature_points();
                  for (unsigned int i = 0; i < cell->get_fe().dofs_per_cell; ++i)
                    // insert the values into the map
                    support_points[local_dof_indices[i]] = points[i];
index a629dcc61ec96fb1d650fa90900b460354f4866d,f1ecde11e7a751498893f0eb046a3f16e15f3b0f..4f3d9e3eea1ba6173d6cdfeea4ede36d212e30af
@@@ -604,14 -600,14 +604,14 @@@ templat
  unsigned int
  DoFTools::count_dofs_with_subdomain_association<hp::DoFHandler<deal_II_dimension,deal_II_dimension+1> >
  (const hp::DoFHandler<deal_II_dimension,deal_II_dimension+1> &,
-  const types::subdomain_id_t);
+  const types::subdomain_id);
 -
 + 
  template
  IndexSet
  DoFTools::dof_indices_with_subdomain_association<hp::DoFHandler<deal_II_dimension,deal_II_dimension+1> >
  (const hp::DoFHandler<deal_II_dimension,deal_II_dimension+1> &,
-  const types::subdomain_id_t);
+  const types::subdomain_id);
 -
 + 
  template
  void
  DoFTools::get_subdomain_association<DoFHandler<deal_II_dimension,deal_II_dimension+1> >
Simple merge
index 6362c78416f9f461d8040821143f4da5a2c05643,75f26683925ed98b433c3c4fd83e6c72d86e04a5..6db6f20fa79595dca694523e85d6b97cc1b29f84
@@@ -2745,7 -2745,7 +2745,7 @@@ namespace h
  
  
    template <>
-   types::global_dof_index DoFHandler<1>::n_boundary_dofs (const std::set<types::boundary_id_t> &boundary_indicators) const
 -  unsigned int DoFHandler<1>::n_boundary_dofs (const std::set<types::boundary_id> &boundary_indicators) const
++  types::global_dof_index DoFHandler<1>::n_boundary_dofs (const std::set<types::boundary_id> &boundary_indicators) const
    {
      Assert (finite_elements != 0, ExcNoFESelected());
  
    }
  
    template <>
-   types::global_dof_index DoFHandler<1,2>::n_boundary_dofs (const std::set<types::boundary_id_t> &) const
 -  unsigned int DoFHandler<1,2>::n_boundary_dofs (const std::set<types::boundary_id> &) const
++  types::global_dof_index DoFHandler<1,2>::n_boundary_dofs (const std::set<types::boundary_id> &) const
    {
      Assert(false,ExcNotImplemented());
      return 0;
@@@ -2820,7 -2819,7 +2820,7 @@@ types::global_dof_index DoFHandler<1,3>
    }
  
    template <>
-   types::global_dof_index DoFHandler<1,3>::n_boundary_dofs (const std::set<types::boundary_id_t> &) const
 -  unsigned int DoFHandler<1,3>::n_boundary_dofs (const std::set<types::boundary_id> &) const
++  types::global_dof_index DoFHandler<1,3>::n_boundary_dofs (const std::set<types::boundary_id> &) const
    {
      Assert(false,ExcNotImplemented());
      return 0;
  
  
    template<int dim, int spacedim>
 -  unsigned int
 +  types::global_dof_index
-   DoFHandler<dim,spacedim>::n_boundary_dofs (const std::set<types::boundary_id_t> &boundary_indicators) const
+   DoFHandler<dim,spacedim>::n_boundary_dofs (const std::set<types::boundary_id> &boundary_indicators) const
    {
      Assert (finite_elements != 0, ExcNoFESelected());
      Assert (boundary_indicators.find (types::internal_face_boundary_id) == boundary_indicators.end(),
  
  
    template <>
-   types::global_dof_index DoFHandler<2,3>::n_boundary_dofs (const std::set<types::boundary_id_t> &) const
 -  unsigned int DoFHandler<2,3>::n_boundary_dofs (const std::set<types::boundary_id> &) const
++  types::global_dof_index DoFHandler<2,3>::n_boundary_dofs (const std::set<types::boundary_id> &) const
    {
      Assert(false,ExcNotImplemented());
      return 0;
@@@ -4143,31 -4140,6 +4143,32 @@@ template <
    }
  
  
-   DoFHandler<dim,spacedim>::get_dof_index (const unsigned int       obj_level,
-                                            const unsigned int       obj_index,
-                                            const unsigned int       fe_index,
-                                            const unsigned int       local_index) const
 +  template <int dim, int spacedim>
 +  template <int structdim>
 +  types::global_dof_index
-     
++  DoFHandler<dim,spacedim>::get_dof_index (const unsigned int,
++                                           const unsigned int,
++                                           const unsigned int,
++                                           const unsigned int) const
 +  {
-   DoFHandler<dim,spacedim>::set_dof_index (const unsigned int obj_level,
-                                            const unsigned int obj_index,
-                                            const unsigned int fe_index,
-                                            const unsigned int local_index,
-                                            const types::global_dof_index global_index) const
++    Assert (false, ExcNotImplemented());
++    return numbers::invalid_unsigned_int;
 +  }
 +
 +
 +  template <int dim, int spacedim>
 +  template <int structdim>
 +  void
-     
++  DoFHandler<dim,spacedim>::set_dof_index (const unsigned int,
++                                           const unsigned int,
++                                           const unsigned int,
++                                           const unsigned int,
++                                           const types::global_dof_index) const
 +  {
++    Assert (false, ExcNotImplemented());
 +  }
 +
 +
    template<int dim, int spacedim>
    void DoFHandler<dim,spacedim>::clear_space ()
    {
index d60d62e20a2d842314b3e6c619743b01cb128757,411b60a7e024012f628ae8e031e995659389affc..c81d6320326b5bb475e593663a542d6b2a469d99
@@@ -867,9 -866,9 +867,9 @@@ namespace MGTool
    {
      const FiniteElement<dim>& fe = dof.get_fe();
      const unsigned int n_comp = fe.n_components();
-  
      Assert ((level>=1) && (level<dof.get_tria().n_levels()),
 -            ExcIndexRange(level, 1, dof.get_tria().n_levels()));
 +              ExcIndexRange(level, 1, dof.get_tria().n_levels()));
  
      const unsigned int fine_dofs = dof.n_dofs(level);
      const unsigned int coarse_dofs = dof.n_dofs(level-1);
      if (component_mask.size() == 0)
        {
          typename MGDoFHandler<dim,spacedim>::cell_iterator
 -          cell = dof.begin(),
 -          endc = dof.end();
 +            cell = dof.begin(),
 +            endc = dof.end();
          for (; cell!=endc; ++cell)
 -          {
 -            const FiniteElement<dim> &fe = cell->get_fe();
 -            const unsigned int level = cell->level();
 -            local_dofs.resize(fe.dofs_per_face);
 -
 -            for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
 -                 ++face_no)
 -              if (cell->at_boundary(face_no) == true)
 -                {
 -                  const typename MGDoFHandler<dim,spacedim>::face_iterator
 -                    face = cell->face(face_no);
 -                  const types::boundary_id bi = face->boundary_indicator();
 +            {
 +              const FiniteElement<dim> &fe = cell->get_fe();
 +              const unsigned int level = cell->level();
 +              local_dofs.resize(fe.dofs_per_face);
 +
 +              for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
 +                   ++face_no)
 +                if (cell->at_boundary(face_no) == true)
 +                  {
 +                        const typename MGDoFHandler<dim,spacedim>::face_iterator
 +                          face = cell->face(face_no);
-                         const types::boundary_id_t bi = face->boundary_indicator();
++                        const types::boundary_id bi = face->boundary_indicator();
                                                     // Face is listed in
                                                     // boundary map
 -                  if (function_map.find(bi) != function_map.end())
 -                    {
 -                      face->get_mg_dof_indices(level, local_dofs);
 -                      for (unsigned int i=0;i<fe.dofs_per_face;++i)
 -                        boundary_indices[level].insert(local_dofs[i]);
 -                    }
 -                }
 -          }
 +                        if (function_map.find(bi) != function_map.end())
 +                          {
 +                            face->get_mg_dof_indices(level, local_dofs);
 +                            for (unsigned int i=0;i<fe.dofs_per_face;++i)
 +                              boundary_indices[level].insert(local_dofs[i]);
 +                          }
 +                  }
 +            }
        }
      else
        {
                                                 // that are non-zero for the
                                                 // components we are interested in,
                                                 // are in fact primitive
 -              for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
 -                {
 -                  const std::vector<bool> &nonzero_component_array
 -                    = cell->get_fe().get_nonzero_components (i);
 -                  for (unsigned int c=0; c<n_components; ++c)
 -                    if ((nonzero_component_array[c] == true)
 -                        &&
 -                        (component_mask[c] == true))
 -                      Assert (cell->get_fe().is_primitive (i),
 -                              ExcMessage ("This function can only deal with requested boundary "
 -                                          "values that correspond to primitive (scalar) base "
 -                                          "elements"));
 -                }
 -
 -              typename MGDoFHandler<dim,spacedim>::face_iterator face = cell->face(face_no);
 -              const types::boundary_id boundary_component = face->boundary_indicator();
 -              if (function_map.find(boundary_component) != function_map.end())
 +                for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
 +                  {
 +                        const std::vector<bool> &nonzero_component_array
 +                          = cell->get_fe().get_nonzero_components (i);
 +                        for (unsigned int c=0; c<n_components; ++c)
 +                          if ((nonzero_component_array[c] == true)
 +                              &&
 +                              (component_mask[c] == true))
 +                            Assert (cell->get_fe().is_primitive (i),
 +                                        ExcMessage ("This function can only deal with requested boundary "
 +                                                            "values that correspond to primitive (scalar) base "
 +                                                            "elements"));
 +                  }
 +
 +                typename MGDoFHandler<dim,spacedim>::face_iterator face = cell->face(face_no);
-                 const types::boundary_id_t boundary_component = face->boundary_indicator();
++                const types::boundary_id boundary_component = face->boundary_indicator();
 +                if (function_map.find(boundary_component) != function_map.end())
                                                   // face is of the right component
 -                {
 +                  {
                                                     // get indices, physical location and
                                                     // boundary values of dofs on this
                                                     // face
          std::fill (cell_dofs.begin(), cell_dofs.end(), false);
  
          for (unsigned int face_nr=0; face_nr<GeometryInfo<dim>::faces_per_cell; ++face_nr)
 -          {
 -            const typename DoFHandler<dim,spacedim>::face_iterator face = cell->face(face_nr);
 -            if (!face->at_boundary())
 -              {
 +            {
 +              const typename DoFHandler<dim,spacedim>::face_iterator face = cell->face(face_nr);
 +              if (!face->at_boundary())
 +                {
                                                   //interior face
 -                const typename MGDoFHandler<dim>::cell_iterator
 -                  neighbor = cell->neighbor(face_nr);
 +                  const typename MGDoFHandler<dim>::cell_iterator
 +                        neighbor = cell->neighbor(face_nr);
  
 -                if (neighbor->level() < cell->level())
 -                  {
 -                    for (unsigned int j=0; j<dofs_per_face; ++j)
 -                      cell_dofs[fe.face_to_cell_index(j,face_nr)] = true;
 +                  if (neighbor->level() < cell->level())
 +                        {
 +                          for (unsigned int j=0; j<dofs_per_face; ++j)
 +                            cell_dofs[fe.face_to_cell_index(j,face_nr)] = true;
  
 -                  }
 -              }
 -          }
 +                        }
 +                }
 +            }
-  
          const unsigned int level = cell->level();
          cell->get_mg_dof_indices (local_dof_indices);
  
    }
  
  
 -template <int dim, int spacedim>
 -void
 -extract_non_interface_dofs (const MGDoFHandler<dim,spacedim> &mg_dof_handler,
 -                            std::vector<std::set<unsigned int> >  &non_interface_dofs)
 -{
 -  Assert (non_interface_dofs.size() == mg_dof_handler.get_tria().n_levels(),
 -          ExcDimensionMismatch (non_interface_dofs.size(),
 -                                mg_dof_handler.get_tria().n_levels()));
 +  template <int dim, int spacedim>
 +  void
 +  extract_non_interface_dofs (const MGDoFHandler<dim,spacedim> &mg_dof_handler,
 +                            std::vector<std::set<unsigned int> >  &non_interface_dofs)
 +  {
 +    Assert (non_interface_dofs.size() == mg_dof_handler.get_tria().n_levels(),
 +              ExcDimensionMismatch (non_interface_dofs.size(),
 +                                                mg_dof_handler.get_tria().n_levels()));
  
 -  const FiniteElement<dim,spacedim> &fe = mg_dof_handler.get_fe();
 +    const FiniteElement<dim,spacedim> &fe = mg_dof_handler.get_fe();
  
 -  const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
 -  const unsigned int   dofs_per_face   = fe.dofs_per_face;
 +    const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
 +    const unsigned int   dofs_per_face   = fe.dofs_per_face;
  
 -  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 -  std::vector<bool> cell_dofs(dofs_per_cell, false);
 -  std::vector<bool> cell_dofs_interface(dofs_per_cell, false);
 +    std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 +    std::vector<bool> cell_dofs(dofs_per_cell, false);
 +    std::vector<bool> cell_dofs_interface(dofs_per_cell, false);
  
 -  typename MGDoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(),
 -                                            endc = mg_dof_handler.end();
 +    typename MGDoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(),
 +                                                                endc = mg_dof_handler.end();
  
  
 -  for (; cell!=endc; ++cell)
 -    {
 -      std::fill (cell_dofs.begin(), cell_dofs.end(), false);
 -      std::fill (cell_dofs_interface.begin(), cell_dofs_interface.end(), false);
 +    for (; cell!=endc; ++cell)
 +      {
 +        std::fill (cell_dofs.begin(), cell_dofs.end(), false);
 +        std::fill (cell_dofs_interface.begin(), cell_dofs_interface.end(), false);
-  
 -      for (unsigned int face_nr=0; face_nr<GeometryInfo<dim>::faces_per_cell; ++face_nr)
 -        {
 -          const typename DoFHandler<dim,spacedim>::face_iterator face = cell->face(face_nr);
 -          if (!face->at_boundary())
 -            {
 +        for (unsigned int face_nr=0; face_nr<GeometryInfo<dim>::faces_per_cell; ++face_nr)
 +            {
 +              const typename DoFHandler<dim,spacedim>::face_iterator face = cell->face(face_nr);
 +              if (!face->at_boundary())
 +                {
                                                 //interior face
 -              const typename MGDoFHandler<dim>::cell_iterator
 -                neighbor = cell->neighbor(face_nr);
 -
 -              if ((neighbor->level() < cell->level()))
 +                  const typename MGDoFHandler<dim>::cell_iterator
 +                        neighbor = cell->neighbor(face_nr);
 +
 +                  if ((neighbor->level() < cell->level()))
 +                        {
 +                          for (unsigned int j=0; j<dofs_per_face; ++j)
 +                            cell_dofs_interface[fe.face_to_cell_index(j,face_nr)] = true;
 +                        }
 +                else
                  {
                    for (unsigned int j=0; j<dofs_per_face; ++j)
 -                    cell_dofs_interface[fe.face_to_cell_index(j,face_nr)] = true;
 +                          cell_dofs[fe.face_to_cell_index(j,face_nr)] = true;
                  }
 -              else
 +                }
 +            else
                {
 -                for (unsigned int j=0; j<dofs_per_face; ++j)
 -                    cell_dofs[fe.face_to_cell_index(j,face_nr)] = true;
 -              }
 -            }
 -          else
 -          {
                                                 //boundary face
                  for (unsigned int j=0; j<dofs_per_face; ++j)
 -                    cell_dofs[fe.face_to_cell_index(j,face_nr)] = true;
 -          }
 -        }
 +                        cell_dofs[fe.face_to_cell_index(j,face_nr)] = true;
 +              }
 +            }
  
 -      const unsigned int level = cell->level();
 -      cell->get_mg_dof_indices (local_dof_indices);
 +        const unsigned int level = cell->level();
 +        cell->get_mg_dof_indices (local_dof_indices);
  
 -      for(unsigned int i=0; i<dofs_per_cell; ++i)
 -          if (cell_dofs[i] && !cell_dofs_interface[i])
 -            non_interface_dofs[level].insert(local_dof_indices[i]);
 -    }
 -}
 +        for(unsigned int i=0; i<dofs_per_cell; ++i)
 +            if (cell_dofs[i] && !cell_dofs_interface[i])
 +              non_interface_dofs[level].insert(local_dof_indices[i]);
 +      }
 +  }
  
  
    template <int dim, int spacedim>
index 0f78643c9768b7b93348afc683fc47fe8b416709,b4d74691ffdbdb655ff97181be452cc7819c9659..ad42c3ce05076364a0f3dd1199cf082be8e77dc5
@@@ -77,15 -77,15 +77,18 @@@ for (deal_II_dimension : DIMENSIONS
          const Table<2,DoFTools::Coupling>&, const Table<2,DoFTools::Coupling>&);
  #endif
  
 -      template void count_dofs_per_component<deal_II_dimension> (
 +      template void count_dofs_per_component (
        const MGDoFHandler<deal_II_dimension>&, std::vector<std::vector<unsigned int> >&,
        bool, std::vector<unsigned int>);
 -      template void count_dofs_per_component<deal_II_dimension> (
 +      template void count_dofs_per_component (
        const MGDoFHandler<deal_II_dimension>&, std::vector<std::vector<unsigned int> >&,
        std::vector<unsigned int>);
 -      template void count_dofs_per_block<deal_II_dimension> (
 +      template void count_dofs_per_block (
        const MGDoFHandler<deal_II_dimension>&, std::vector<std::vector<unsigned int> >&,
        std::vector<unsigned int>);
++      template void count_dofs_per_block (
++      const DoFHandler<deal_II_dimension>&, std::vector<std::vector<unsigned int> >&,
++      std::vector<unsigned int>);
  
  #if deal_II_dimension > 1
        template void make_boundary_list(
  #endif
  
  #if deal_II_dimension == 3
 -      template void count_dofs_per_block<1,3> (
 +      template void count_dofs_per_block (
        const MGDoFHandler<1,3>&,
        std::vector<std::vector<unsigned int> >&, std::vector<unsigned int>);
++
++      template void count_dofs_per_block (
++      const DoFHandler<1,3>&,
++      std::vector<std::vector<unsigned int> >&, std::vector<unsigned int>);
++      template void count_dofs_per_block (
++      const DoFHandler<2,3>&,
++      std::vector<std::vector<unsigned int> >&, std::vector<unsigned int>);
++      template void count_dofs_per_block (
++      const DoFHandler<1,2>&,
++      std::vector<std::vector<unsigned int> >&, std::vector<unsigned int>);
  #endif
        \}
    }
Simple merge

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.