From: heister <heister@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Wed, 19 Dec 2012 17:52:41 +0000 (+0000)
Subject: merge from upstream
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d251ee32a0bb0ddd40ca0e513bda8d1142c340f1;p=dealii-svn.git

merge from upstream

git-svn-id: https://svn.dealii.org/branches/branch_bigger_global_dof_indices_3@27830 0785d39b-7218-0410-832d-ea1e28bc413d
---

d251ee32a0bb0ddd40ca0e513bda8d1142c340f1
diff --cc deal.II/include/deal.II/lac/constraint_matrix.h
index e665929d2b,103433b068..03143da3cc
--- a/deal.II/include/deal.II/lac/constraint_matrix.h
+++ b/deal.II/include/deal.II/lac/constraint_matrix.h
@@@ -1903,9 -1903,18 +1903,18 @@@ ConstraintMatrix::set_inhomogeneity (co
  
  
  
+ inline
+ unsigned int
+ ConstraintMatrix::n_constraints () const
+ {
+   return lines.size();
+ }
+ 
+ 
+ 
  inline
  bool
 -ConstraintMatrix::is_constrained (const unsigned int index) const
 +ConstraintMatrix::is_constrained (const types::global_dof_index index) const
  {
    const unsigned int line_index = calculate_line_index(index);
    return ((line_index < lines_cache.size())
diff --cc deal.II/include/deal.II/lac/sparse_matrix.h
index 2d77baf4b0,9092530630..add3c01108
--- a/deal.II/include/deal.II/lac/sparse_matrix.h
+++ b/deal.II/include/deal.II/lac/sparse_matrix.h
@@@ -610,24 -484,17 +484,17 @@@ public
    explicit SparseMatrix (const SparsityPattern &sparsity);
  
    /**
-    * Copy constructor: initialize
-    * the matrix with the identity
-    * matrix. This constructor will
-    * throw an exception if the
-    * sizes of the sparsity pattern
-    * and the identity matrix do not
-    * coincide, or if the sparsity
-    * pattern does not provide for
-    * nonzero entries on the entire
-    * diagonal.
+    * Copy constructor: initialize the matrix with the identity matrix. This
+    * constructor will throw an exception if the sizes of the sparsity pattern
+    * and the identity matrix do not coincide, or if the sparsity pattern does
+    * not provide for nonzero entries on the entire diagonal.
     */
    SparseMatrix (const SparsityPattern &sparsity,
 -                const IdentityMatrix &id);
 +                const IdentityMatrix  &id);
  
    /**
-    * Destructor. Free all memory, but do not
-    * release the memory of the sparsity
-    * structure.
+    * Destructor. Free all memory, but do not release the memory of the
+    * sparsity structure.
     */
    virtual ~SparseMatrix ();
  
@@@ -653,34 -510,21 +510,21 @@@
    SparseMatrix<number> &operator = (const SparseMatrix<number> &);
  
    /**
-    * Copy operator: initialize
-    * the matrix with the identity
-    * matrix. This operator will
-    * throw an exception if the
-    * sizes of the sparsity pattern
-    * and the identity matrix do not
-    * coincide, or if the sparsity
-    * pattern does not provide for
-    * nonzero entries on the entire
-    * diagonal.
+    * Copy operator: initialize the matrix with the identity matrix. This
+    * operator will throw an exception if the sizes of the sparsity pattern and
+    * the identity matrix do not coincide, or if the sparsity pattern does not
+    * provide for nonzero entries on the entire diagonal.
     */
    SparseMatrix<number> &
 -  operator= (const IdentityMatrix &id);
 +  operator= (const IdentityMatrix  &id);
  
    /**
-    * This operator assigns a scalar to
-    * a matrix. Since this does usually
-    * not make much sense (should we set
-    * all matrix entries to this value?
-    * Only the nonzero entries of the
-    * sparsity pattern?), this operation
-    * is only allowed if the actual
-    * value to be assigned is zero. This
-    * operator only exists to allow for
-    * the obvious notation
-    * <tt>matrix=0</tt>, which sets all
-    * elements of the matrix to zero,
-    * but keep the sparsity pattern
+    * This operator assigns a scalar to a matrix. Since this does usually not
+    * make much sense (should we set all matrix entries to this value?  Only
+    * the nonzero entries of the sparsity pattern?), this operation is only
+    * allowed if the actual value to be assigned is zero. This operator only
+    * exists to allow for the obvious notation <tt>matrix=0</tt>, which sets
+    * all elements of the matrix to zero, but keep the sparsity pattern
     * previously used.
     */
    SparseMatrix &operator = (const double d);
diff --cc deal.II/include/deal.II/lac/sparsity_pattern.h
index 0eb8828807,fe87da3f19..5c13c278e0
--- a/deal.II/include/deal.II/lac/sparsity_pattern.h
+++ b/deal.II/include/deal.II/lac/sparsity_pattern.h
@@@ -496,46 -385,26 +385,26 @@@ public
    /**
     * Make a copy with extra off-diagonals.
     *
-    * This constructs objects intended for
-    * the application of the ILU(n)-method
-    * or other incomplete decompositions.
-    * Therefore, additional to the original
-    * entry structure, space for
-    * <tt>extra_off_diagonals</tt>
-    * side-diagonals is provided on both
-    * sides of the main diagonal.
-    *
-    * <tt>max_per_row</tt> is the
-    * maximum number of nonzero
-    * elements per row which this
-    * structure is to hold. It is
-    * assumed that this number is
-    * sufficiently large to
-    * accommodate both the elements
-    * in <tt>original</tt> as well
-    * as the new off-diagonal
-    * elements created by this
-    * constructor. You will usually
-    * want to give the same number
-    * as you gave for
-    * <tt>original</tt> plus the
-    * number of side diagonals times
-    * two. You may however give a
-    * larger value if you wish to
-    * add further nonzero entries
-    * for the decomposition based on
-    * other criteria than their
-    * being on side-diagonals.
-    *
-    * This function requires that
-    * <tt>original</tt> refers to a
-    * quadratic matrix structure.
-    * It must be compressed. The
-    * matrix structure is not
-    * compressed after this function
-    * finishes.
+    * This constructs objects intended for the application of the ILU(n)-method
+    * or other incomplete decompositions.  Therefore, additional to the
+    * original entry structure, space for <tt>extra_off_diagonals</tt>
+    * side-diagonals is provided on both sides of the main diagonal.
+    *
+    * <tt>max_per_row</tt> is the maximum number of nonzero elements per row
+    * which this structure is to hold. It is assumed that this number is
+    * sufficiently large to accommodate both the elements in <tt>original</tt>
+    * as well as the new off-diagonal elements created by this constructor. You
+    * will usually want to give the same number as you gave for
+    * <tt>original</tt> plus the number of side diagonals times two. You may
+    * however give a larger value if you wish to add further nonzero entries
+    * for the decomposition based on other criteria than their being on
+    * side-diagonals.
+    *
+    * This function requires that <tt>original</tt> refers to a quadratic
+    * matrix structure.  It must be compressed. The matrix structure is not
+    * compressed after this function finishes.
     */
 -  SparsityPattern (const SparsityPattern &original,
 +  SparsityPattern (const SparsityPattern  &original,
                     const unsigned int      max_per_row,
                     const unsigned int      extra_off_diagonals);
  
diff --cc deal.II/source/dofs/dof_tools.cc
index 30697f7ee0,619d478814..fcfea2260e
--- a/deal.II/source/dofs/dof_tools.cc
+++ b/deal.II/source/dofs/dof_tools.cc
@@@ -382,19 -367,14 +367,14 @@@ namespace DoFTool
        };
  #endif
  
 -    std::vector<unsigned int> dofs_on_this_face;
 +    std::vector<types::global_dof_index> dofs_on_this_face;
      dofs_on_this_face.reserve (max_dofs_per_face(dof));
  
-     // loop over all faces to check
-     // whether they are at a
-     // boundary. note that we need not
-     // take special care of single
-     // lines (using
-     // @p{cell->has_boundary_lines}),
-     // since we do not support
-     // boundaries of dimension dim-2,
-     // and so every boundary line is
-     // also part of a boundary face.
+     // loop over all faces to check whether they are at a boundary. note
+     // that we need not take special care of single lines (using
+     // @p{cell->has_boundary_lines}), since we do not support boundaries of
+     // dimension dim-2, and so every boundary line is also part of a
+     // boundary face.
      typename DH::active_cell_iterator cell = dof.begin_active(),
                                        endc = dof.end();
      for (; cell!=endc; ++cell)
@@@ -445,10 -422,9 +422,9 @@@
                cell = cell->child(direction);
  
              const unsigned int dofs_per_vertex = cell->get_fe().dofs_per_vertex;
 -            std::vector<unsigned int> boundary_dof_boundary_indices (dofs_per_vertex);
 +            std::vector<types::global_dof_index> boundary_dof_boundary_indices (dofs_per_vertex);
  
-             // next get boundary mapped dof
-             // indices of boundary dofs
+             // next get boundary mapped dof indices of boundary dofs
              for (unsigned int i=0; i<dofs_per_vertex; ++i)
                boundary_dof_boundary_indices[i]
                  = dof_to_boundary_mapping[cell->vertex_dof_index(direction,i)];
@@@ -1217,10 -1172,9 +1172,9 @@@
    {
      namespace
      {
-       inline
-       bool
+       inline bool
        check_master_dof_list (const FullMatrix<double> &face_interpolation_matrix,
 -                             const std::vector<unsigned int> &master_dof_list)
 +                             const std::vector<types::global_dof_index> &master_dof_list)
        {
          const unsigned int N = master_dof_list.size();
  
@@@ -1387,43 -1302,28 +1302,28 @@@
                  (fe2.dofs_per_quad <= fe1.dofs_per_quad),
                  ExcInternalError());
  
-         // the idea here is to designate as
-         // many DoFs in fe1 per object
-         // (vertex, line, quad) as master as
-         // there are such dofs in fe2
-         // (indices are int, because we want
-         // to avoid the 'unsigned int < 0 is
-         // always false warning for the cases
-         // at the bottom in 1d and 2d)
+         // the idea here is to designate as many DoFs in fe1 per object
+         // (vertex, line, quad) as master as there are such dofs in fe2
+         // (indices are int, because we want to avoid the 'unsigned int < 0
+         // is always false warning for the cases at the bottom in 1d and
+         // 2d)
          //
-         // as mentioned in the paper, it is
-         // not always easy to find a set of
-         // master dofs that produces an
-         // invertible matrix. to this end, we
-         // check in each step whether the
-         // matrix is still invertible and
-         // simply discard this dof if the
-         // matrix is not invertible anymore.
+         // as mentioned in the paper, it is not always easy to find a set
+         // of master dofs that produces an invertible matrix. to this end,
+         // we check in each step whether the matrix is still invertible and
+         // simply discard this dof if the matrix is not invertible anymore.
          //
-         // the cases where we did have
-         // trouble in the past were with
-         // adding more quad dofs when Q3 and
-         // Q4 elements meet at a refined face
-         // in 3d (see the hp/crash_12 test
-         // that tests that we can do exactly
-         // this, and failed before we had
-         // code to compensate for this
-         // case). the other case are system
-         // elements: if we have say a Q1Q2 vs
-         // a Q2Q3 element, then we can't just
-         // take all master dofs on a line
-         // from a single base element, since
-         // the shape functions of that base
-         // element are independent of that of
-         // the other one. this latter case
-         // shows up when running
-         // hp/hp_constraints_q_system_06
+         // the cases where we did have trouble in the past were with adding
+         // more quad dofs when Q3 and Q4 elements meet at a refined face in
+         // 3d (see the hp/crash_12 test that tests that we can do exactly
+         // this, and failed before we had code to compensate for this
+         // case). the other case are system elements: if we have say a Q1Q2
+         // vs a Q2Q3 element, then we can't just take all master dofs on a
+         // line from a single base element, since the shape functions of
+         // that base element are independent of that of the other one. this
+         // latter case shows up when running hp/hp_constraints_q_system_06
+ 
 -        std::vector<unsigned int> master_dof_list;
 +        std::vector<types::global_dof_index> master_dof_list;
          unsigned int index = 0;
          for (int v=0;
               v<static_cast<signed int>(GeometryInfo<dim>::vertices_per_face);
@@@ -1523,10 -1411,9 +1411,9 @@@
          AssertDimension (index, fe1.dofs_per_face);
          AssertDimension (master_dof_list.size(), fe2.dofs_per_face);
  
-         // finally copy the list into the
-         // mask
+         // finally copy the list into the mask
          std::fill (master_dof_mask.begin(), master_dof_mask.end(), false);
 -        for (std::vector<unsigned int>::const_iterator i=master_dof_list.begin();
 +        for (std::vector<types::global_dof_index>::const_iterator i=master_dof_list.begin();
               i!=master_dof_list.end(); ++i)
            master_dof_mask[*i] = true;
        }
@@@ -1802,24 -1661,18 +1661,18 @@@
  
  
        /**
-        * Copy constraints into a constraint
-        * matrix object.
+        * Copy constraints into a constraint matrix object.
         *
-        * This function removes zero
-        * constraints and those, which
-        * constrain a DoF which was
-        * already eliminated in one of
-        * the previous steps of the hp
-        * hanging node procedure.
+        * This function removes zero constraints and those, which constrain
+        * a DoF which was already eliminated in one of the previous steps of
+        * the hp hanging node procedure.
         *
-        * It also suppresses very small
-        * entries in the constraint matrix to
-        * avoid making the sparsity pattern
-        * fuller than necessary.
+        * It also suppresses very small entries in the constraint matrix to
+        * avoid making the sparsity pattern fuller than necessary.
         */
        void
 -      filter_constraints (const std::vector<unsigned int> &master_dofs,
 -                          const std::vector<unsigned int> &slave_dofs,
 +      filter_constraints (const std::vector<types::global_dof_index> &master_dofs,
 +                          const std::vector<types::global_dof_index> &slave_dofs,
                            const FullMatrix<double> &face_constraints,
                            ConstraintMatrix &constraints)
        {
@@@ -1947,28 -1781,26 +1781,26 @@@
  
  
      void
 -    make_hp_hanging_node_constraints (const dealii::hp::DoFHandler<1> & /*dof_handler*/,
 -                                      ConstraintMatrix &        /*constraints*/)
 +    make_hp_hanging_node_constraints (const dealii::hp::DoFHandler<1> &/*dof_handler*/,
 +                                      ConstraintMatrix        &/*constraints*/)
      {
-       // we may have to compute
-       // constraints for
-       // vertices. gotta think about
+       // we may have to compute constraints for vertices. gotta think about
        // that a bit more
- //TODO[WB]: think about what to do here...
+ 
+       //TODO[WB]: think about what to do here...
      }
  
  
  
      void
 -    make_oldstyle_hanging_node_constraints (const dealii::hp::DoFHandler<1> & /*dof_handler*/,
 -                                            ConstraintMatrix &        /*constraints*/,
 +    make_oldstyle_hanging_node_constraints (const dealii::hp::DoFHandler<1> &/*dof_handler*/,
 +                                            ConstraintMatrix        &/*constraints*/,
                                              dealii::internal::int2type<1>)
      {
-       // we may have to compute
-       // constraints for
-       // vertices. gotta think about
+       // we may have to compute constraints for vertices. gotta think about
        // that a bit more
- //TODO[WB]: think about what to do here...
+ 
+       //TODO[WB]: think about what to do here...
      }
  
  
@@@ -2057,25 -1886,18 +1886,18 @@@
  
        const unsigned int spacedim = DH::space_dimension;
  
 -      std::vector<unsigned int> dofs_on_mother;
 -      std::vector<unsigned int> dofs_on_children;
 +      std::vector<types::global_dof_index> dofs_on_mother;
 +      std::vector<types::global_dof_index> dofs_on_children;
  
-       // loop over all lines; only on
-       // lines there can be constraints.
-       // We do so by looping over all
-       // active cells and checking
-       // whether any of the faces are
-       // refined which can only be from
-       // the neighboring cell because
-       // this one is active. In that
-       // case, the face is subject to
+       // loop over all lines; only on lines there can be constraints. We do
+       // so by looping over all active cells and checking whether any of
+       // the faces are refined which can only be from the neighboring cell
+       // because this one is active. In that case, the face is subject to
        // constraints
        //
-       // note that even though we may
-       // visit a face twice if the
-       // neighboring cells are equally
-       // refined, we can only visit each
-       // face with hanging nodes once
+       // note that even though we may visit a face twice if the neighboring
+       // cells are equally refined, we can only visit each face with
+       // hanging nodes once
        typename DH::active_cell_iterator cell = dof_handler.begin_active(),
                                          endc = dof_handler.end();
        for (; cell!=endc; ++cell)
@@@ -2211,25 -2010,18 +2010,18 @@@
      {
        const unsigned int dim = 3;
  
 -      std::vector<unsigned int> dofs_on_mother;
 -      std::vector<unsigned int> dofs_on_children;
 +      std::vector<types::global_dof_index> dofs_on_mother;
 +      std::vector<types::global_dof_index> dofs_on_children;
  
-       // loop over all quads; only on
-       // quads there can be constraints.
-       // We do so by looping over all
-       // active cells and checking
-       // whether any of the faces are
-       // refined which can only be from
-       // the neighboring cell because
-       // this one is active. In that
-       // case, the face is subject to
+       // loop over all quads; only on quads there can be constraints. We do
+       // so by looping over all active cells and checking whether any of
+       // the faces are refined which can only be from the neighboring cell
+       // because this one is active. In that case, the face is subject to
        // constraints
        //
-       // note that even though we may
-       // visit a face twice if the
-       // neighboring cells are equally
-       // refined, we can only visit each
-       // face with hanging nodes once
+       // note that even though we may visit a face twice if the neighboring
+       // cells are equally refined, we can only visit each face with
+       // hanging nodes once
        typename DH::active_cell_iterator cell = dof_handler.begin_active(),
                                          endc = dof_handler.end();
        for (; cell!=endc; ++cell)
@@@ -2459,30 -2221,21 +2221,21 @@@
        const unsigned int spacedim = DH::space_dimension;
  
  
-       // a matrix to be used for
-       // constraints below. declared
-       // here and simply resized down
-       // below to avoid permanent
-       // re-allocation of memory
+       // a matrix to be used for constraints below. declared here and
+       // simply resized down below to avoid permanent re-allocation of
+       // memory
        FullMatrix<double> constraint_matrix;
  
-       // similarly have arrays that
-       // will hold master and slave
-       // dof numbers, as well as a
-       // scratch array needed for the
-       // complicated case below
+       // similarly have arrays that will hold master and slave dof numbers,
+       // as well as a scratch array needed for the complicated case below
 -      std::vector<unsigned int> master_dofs;
 -      std::vector<unsigned int> slave_dofs;
 -      std::vector<unsigned int> scratch_dofs;
 +      std::vector<types::global_dof_index> master_dofs;
 +      std::vector<types::global_dof_index> slave_dofs;
 +      std::vector<types::global_dof_index> scratch_dofs;
  
-       // caches for the face and
-       // subface interpolation
-       // matrices between different
-       // (or the same) finite
-       // elements. we compute them
-       // only once, namely the first
-       // time they are needed, and
-       // then just reuse them
+       // caches for the face and subface interpolation matrices between
+       // different (or the same) finite elements. we compute them only
+       // once, namely the first time they are needed, and then just reuse
+       // them
        Table<2,std_cxx1x::shared_ptr<FullMatrix<double> > >
        face_interpolation_matrices (n_finite_elements (dof_handler),
                                     n_finite_elements (dof_handler));
@@@ -3836,9 -3318,8 +3318,8 @@@
                    ExcInternalError());
          }
  
-       // then loop over all cells and do
-       // the work
+       // then loop over all cells and do the work
 -      std::vector<unsigned int> indices;
 +      std::vector<types::global_dof_index> indices;
        for (typename DH::active_cell_iterator c=dof.begin_active();
             c!=dof.end(); ++ c)
          if (c->is_locally_owned())
@@@ -4112,10 -3573,9 +3573,9 @@@
      for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
        local_selected_dofs[i] = component_mask[local_component_asssociation[i]];
  
-     // then loop over all cells and do
-     // work
+     // then loop over all cells and do work
 -    std::vector<unsigned int> indices(fe.dofs_per_cell);
 -    typename MGDoFHandler<dim,spacedim>::cell_iterator c;
 +    std::vector<types::global_dof_index> indices(fe.dofs_per_cell);
 +    typename DH::cell_iterator c;
      for (c = dof.begin(level) ; c != dof.end(level) ; ++ c)
        {
          c->get_mg_dof_indices(indices);
@@@ -4199,26 -3656,20 +3656,20 @@@
           (component_mask.n_selected_components(n_components(dof_handler)) !=
            n_components(dof_handler)));
  
 -    std::vector<unsigned int> face_dof_indices;
 +    std::vector<types::global_dof_index> face_dof_indices;
      face_dof_indices.reserve (max_dofs_per_face(dof_handler));
  
-     // now loop over all cells and
-     // check whether their faces are at
-     // the boundary. note that we need
-     // not take special care of single
-     // lines being at the boundary
-     // (using
-     // @p{cell->has_boundary_lines}),
-     // since we do not support
-     // boundaries of dimension dim-2,
-     // and so every isolated boundary
-     // line is also part of a boundary
-     // face which we will be visiting
-     // sooner or later
+     // now loop over all cells and check whether their faces are at the
+     // boundary. note that we need not take special care of single lines
+     // being at the boundary (using @p{cell->has_boundary_lines}), since we
+     // do not support boundaries of dimension dim-2, and so every isolated
+     // boundary line is also part of a boundary face which we will be
+     // visiting sooner or later
      for (typename DH::active_cell_iterator cell=dof_handler.begin_active();
           cell!=dof_handler.end(); ++cell)
-       // only work on cells that are either locally owned or at
-       // least ghost cells
+ 
+       // only work on cells that are either locally owned or at least ghost
+       // cells
        if (cell->is_artificial() == false)
          for (unsigned int face=0;
               face<GeometryInfo<DH::dimension>::faces_per_cell; ++face)
@@@ -4318,26 -3755,18 +3755,18 @@@
      const bool check_vector_component
        = (component_mask.represents_the_all_selected_mask() == false);
  
-     // clear and reset array by default
-     // values
+     // clear and reset array by default values
      selected_dofs.clear ();
      selected_dofs.resize (dof_handler.n_dofs(), false);
 -    std::vector<unsigned int> cell_dof_indices;
 +    std::vector<types::global_dof_index> cell_dof_indices;
      cell_dof_indices.reserve (max_dofs_per_cell(dof_handler));
  
-     // now loop over all cells and
-     // check whether their faces are at
-     // the boundary. note that we need
-     // not take special care of single
-     // lines being at the boundary
-     // (using
-     // @p{cell->has_boundary_lines}),
-     // since we do not support
-     // boundaries of dimension dim-2,
-     // and so every isolated boundary
-     // line is also part of a boundary
-     // face which we will be visiting
-     // sooner or later
+     // now loop over all cells and check whether their faces are at the
+     // boundary. note that we need not take special care of single lines
+     // being at the boundary (using @p{cell->has_boundary_lines}), since we
+     // do not support boundaries of dimension dim-2, and so every isolated
+     // boundary line is also part of a boundary face which we will be
+     // visiting sooner or later
      for (typename DH::active_cell_iterator cell=dof_handler.begin_active();
           cell!=dof_handler.end(); ++cell)
        for (unsigned int face=0;
@@@ -4537,12 -3955,11 +3955,11 @@@
      // preset all values by false
      std::fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false);
  
 -    std::vector<unsigned int> local_dof_indices;
 +    std::vector<types::global_dof_index> local_dof_indices;
      local_dof_indices.reserve (max_dofs_per_cell(dof_handler));
  
-     // this function is similar to the
-     // make_sparsity_pattern function,
-     // see there for more information
+     // this function is similar to the make_sparsity_pattern function, see
+     // there for more information
      typename DH::active_cell_iterator
      cell = dof_handler.begin_active(),
      endc = dof_handler.end();
@@@ -4579,14 -3996,11 +3996,11 @@@
      // collect all the locally owned dofs
      dof_set = dof_handler.locally_owned_dofs();
  
-     // add the DoF on the adjacent ghost cells
-     // to the IndexSet, cache them in a
-     // set. need to check each dof manually
-     // because we can't be sure that the dof
-     // range of locally_owned_dofs is really
-     // contiguous.
+     // add the DoF on the adjacent ghost cells to the IndexSet, cache them
+     // in a set. need to check each dof manually because we can't be sure
+     // that the dof range of locally_owned_dofs is really contiguous.
 -    std::vector<unsigned int> dof_indices;
 -    std::set<unsigned int> global_dof_indices;
 +    std::vector<types::global_dof_index> dof_indices;
 +    std::set<types::global_dof_index> global_dof_indices;
  
      typename DH::active_cell_iterator cell = dof_handler.begin_active(),
                                        endc = dof_handler.end();
@@@ -4618,14 -4032,11 +4032,11 @@@
      // collect all the locally owned dofs
      dof_set = dof_handler.locally_owned_dofs();
  
-     // add the DoF on the adjacent ghost cells
-     // to the IndexSet, cache them in a
-     // set. need to check each dof manually
-     // because we can't be sure that the dof
-     // range of locally_owned_dofs is really
-     // contiguous.
+     // add the DoF on the adjacent ghost cells to the IndexSet, cache them
+     // in a set. need to check each dof manually because we can't be sure
+     // that the dof range of locally_owned_dofs is really contiguous.
 -    std::vector<unsigned int> dof_indices;
 -    std::set<unsigned int> global_dof_indices;
 +    std::vector<types::global_dof_index> dof_indices;
 +    std::set<types::global_dof_index> global_dof_indices;
  
      typename DH::active_cell_iterator cell = dof_handler.begin_active(),
                                        endc = dof_handler.end();
@@@ -4733,19 -4142,15 +4142,15 @@@
      std::fill_n (subdomain_association.begin(), dof_handler.n_dofs(),
                   types::invalid_subdomain_id);
  
 -    std::vector<unsigned int> local_dof_indices;
 +    std::vector<types::global_dof_index> local_dof_indices;
      local_dof_indices.reserve (max_dofs_per_cell(dof_handler));
  
-     // pseudo-randomly assign variables
-     // which lie on the interface
-     // between subdomains to each of
-     // the two or more
+     // pseudo-randomly assign variables which lie on the interface between
+     // subdomains to each of the two or more
      bool coin_flip = true;
  
-     // loop over all cells and record
-     // which subdomain a DoF belongs
-     // to. toss a coin in case it is on
-     // an interface
+     // loop over all cells and record which subdomain a DoF belongs to.
+     // toss a coin in case it is on an interface
      typename DH::active_cell_iterator
      cell = dof_handler.begin_active(),
      endc = dof_handler.end();
@@@ -4824,18 -4223,15 +4223,15 @@@
  
      IndexSet index_set (dof_handler.n_dofs());
  
 -    std::vector<unsigned int> local_dof_indices;
 +    std::vector<types::global_dof_index> local_dof_indices;
      local_dof_indices.reserve (max_dofs_per_cell(dof_handler));
  
-     // first generate an unsorted list of all
-     // indices which we fill from the back. could
-     // also insert them directly into the
-     // IndexSet, but that inserts indices in the
-     // middle, which is an O(n^2) algorithm and
-     // hence too expensive. Could also use
-     // std::set, but that is in general more
-     // expensive than a vector
+     // first generate an unsorted list of all indices which we fill from
+     // the back. could also insert them directly into the IndexSet, but
+     // that inserts indices in the middle, which is an O(n^2) algorithm and
+     // hence too expensive. Could also use std::set, but that is in general
+     // more expensive than a vector
 -    std::vector<unsigned int> subdomain_indices;
 +    std::vector<types::global_dof_index> subdomain_indices;
  
      typename DH::active_cell_iterator
      cell = dof_handler.begin_active(),
@@@ -4965,14 -4361,14 +4361,14 @@@
      void
      resolve_components (const hp::FECollection<dim,spacedim> &fe_collection,
                          const std::vector<unsigned char> &dofs_by_component,
 -                        const std::vector<unsigned int> &target_component,
 +                        const std::vector<unsigned int>  &target_component,
                          const bool                        only_once,
 -                        std::vector<unsigned int>        &dofs_per_component,
 +                        std::vector<types::global_dof_index>        &dofs_per_component,
                          unsigned int                     &component)
      {
-       // assert that all elements in the collection have the same
-       // structure (base elements and multiplicity, components per base
-       // element) and then simply call the function above
+       // assert that all elements in the collection have the same structure
+       // (base elements and multiplicity, components per base element) and
+       // then simply call the function above
        for (unsigned int fe=1; fe<fe_collection.size(); ++fe)
          {
            Assert (fe_collection[fe].n_components() == fe_collection[0].n_components(),
@@@ -5587,13 -4903,11 +4903,11 @@@
          unsigned int n_parameters_on_fine_grid=0;
          if (true)
            {
-             // have a flag for each dof on
-             // the fine grid and set it
-             // to true if this is an
-             // interesting dof. finally count
-             // how many true's there
+             // have a flag for each dof on the fine grid and set it to true
+             // if this is an interesting dof. finally count how many true's
+             // there
              std::vector<bool> dof_is_interesting (fine_grid.n_dofs(), false);
 -            std::vector<unsigned int>  local_dof_indices (fine_fe.dofs_per_cell);
 +            std::vector<types::global_dof_index>  local_dof_indices (fine_fe.dofs_per_cell);
  
              for (typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator
                   cell=fine_grid.begin_active();
diff --cc deal.II/source/fe/fe_values.cc
index a1f22a30b2,5a63aa75f2..028fccb660
--- a/deal.II/source/fe/fe_values.cc
+++ b/deal.II/source/fe/fe_values.cc
@@@ -2036,85 -2371,34 +2371,34 @@@ template <int dim, int spacedim
  template <class InputVector, typename number>
  void FEValuesBase<dim,spacedim>::get_function_values (
    const InputVector &fe_function,
 -  const VectorSlice<const std::vector<unsigned int> > &indices,
 +  const VectorSlice<const std::vector<types::global_dof_index> > &indices,
    std::vector<Vector<number> > &values) const
  {
-   // One value per quadrature point
-   Assert (n_quadrature_points == values.size(),
-           ExcDimensionMismatch(values.size(), n_quadrature_points));
- 
-   const unsigned int n_components = fe->n_components();
- 
-   // Size of indices must be a
-   // multiple of dofs_per_cell such
-   // that an integer number of
-   // function values is generated in
-   // each point.
+   // Size of indices must be a multiple of dofs_per_cell such that an integer
+   // number of function values is generated in each point.
    Assert (indices.size() % dofs_per_cell == 0,
            ExcNotMultiple(indices.size(), dofs_per_cell));
- 
-   // The number of components of the
-   // result may be a multiple of the
-   // number of components of the
-   // finite element
-   const unsigned int result_components = indices.size() * n_components / dofs_per_cell;
- 
-   for (unsigned int i=0; i<values.size(); ++i)
-     Assert (values[i].size() == result_components,
-             ExcDimensionMismatch(values[i].size(), result_components));
- 
-   // If the result has more
-   // components than the finite
-   // element, we need this number for
-   // loops filling all components
-   const unsigned int component_multiple = result_components / n_components;
- 
    Assert (this->update_flags & update_values, ExcAccessToUninitializedField());
  
-   // initialize with zero
-   for (unsigned int i=0; i<values.size(); ++i)
-     std::fill_n (values[i].begin(), values[i].size(), 0);
- 
-   // add up contributions of trial
-   // functions. now check whether the
-   // shape function is primitive or
-   // not. if it is, then set its only
-   // non-zero component, otherwise
-   // loop over components
-   for (unsigned int mc = 0; mc < component_multiple; ++mc)
-     for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
-       {
-         const double value = get_vector_element (fe_function, indices[shape_func+mc*dofs_per_cell]);
-         if (value == 0.)
-           continue;
- 
-         if (fe->is_primitive(shape_func))
-           {
-             const unsigned int comp = fe->system_to_component_index(shape_func).first
-                                       + mc * n_components;
-             const unsigned int
-             row = this->shape_function_to_row_table[shape_func * fe->n_components() + comp];
- 
-             const double *shape_value_ptr = &this->shape_values(row, 0);
-             for (unsigned int point=0; point<n_quadrature_points; ++point)
-               values[point](comp) += value * *shape_value_ptr++;
-           }
-         else
-           for (unsigned int c=0; c<n_components; ++c)
-             {
-               if (fe->get_nonzero_components(shape_func)[c] == false)
-                 continue;
- 
-               const unsigned int
-               row = this->shape_function_to_row_table[shape_func * fe->n_components() + c];
- 
-               const double *shape_value_ptr = &this->shape_values(row, 0);
-               const unsigned int comp = c + mc * n_components;
- 
-               for (unsigned int point=0; point<n_quadrature_points; ++point)
-                 values[point](comp) += value * *shape_value_ptr++;
-             }
-       }
+   VectorSlice<std::vector<Vector<number> > > val(values);
+   if (indices.size() <= 100)
+     {
+       double dof_values[100];
+       for (unsigned int i=0; i<dofs_per_cell; ++i)
+         dof_values[i] = get_vector_element (fe_function, indices[i]);
+       internal::do_function_values(&dof_values[0], this->shape_values, *fe,
+                                    this->shape_function_to_row_table, val,
+                                    false, indices.size()/dofs_per_cell);
+     }
+   else
+     {
+       Vector<double> dof_values(100);
+       for (unsigned int i=0; i<dofs_per_cell; ++i)
+         dof_values[i] = get_vector_element (fe_function, indices[i]);
+       internal::do_function_values(dof_values.begin(), this->shape_values, *fe,
+                                    this->shape_function_to_row_table, val,
+                                    false, indices.size()/dofs_per_cell);
+     }
  }
  
  
@@@ -2587,41 -2581,28 +2581,28 @@@ template <int dim, int spacedim
  template <class InputVector>
  void FEValuesBase<dim,spacedim>::get_function_hessians (
    const InputVector &fe_function,
 -  const VectorSlice<const std::vector<unsigned int> > &indices,
 +  const VectorSlice<const std::vector<types::global_dof_index> > &indices,
    std::vector<Tensor<2,spacedim> > &hessians) const
  {
-   Assert (this->update_flags & update_second_derivatives, ExcAccessToUninitializedField());
-   // This function fills a single
-   // component only
-   Assert (fe->n_components() == 1,
-           ExcDimensionMismatch(fe->n_components(), 1));
-   // One index for each dof
-   Assert (indices.size() == dofs_per_cell,
-           ExcDimensionMismatch(indices.size(), dofs_per_cell));
-   // This vector has one entry for
-   // each quadrature point
-   Assert (hessians.size() == n_quadrature_points,
-           ExcDimensionMismatch(hessians.size(), n_quadrature_points));
- 
-   // initialize with zero
-   std::fill_n (hessians.begin(), n_quadrature_points, Tensor<2,spacedim>());
- 
-   // add up contributions of trial
-   // functions. note that here we
-   // deal with scalar finite
-   // elements, so no need to check
-   // for non-primitivity of shape
-   // functions
-   for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
+   Assert (this->update_flags & update_second_derivatives,
+           ExcAccessToUninitializedField());
+   AssertDimension (fe_function.size(), present_cell->n_dofs_for_dof_handler());
+   AssertDimension (indices.size(), dofs_per_cell);
+   if (dofs_per_cell <= 100)
+     {
+       double dof_values[100];
+       for (unsigned int i=0; i<dofs_per_cell; ++i)
+         dof_values[i] = get_vector_element (fe_function, indices[i]);
+       internal::do_function_derivatives(&dof_values[0], this->shape_hessians,
+                                         hessians);
+     }
+   else
      {
-       const double value = get_vector_element (fe_function, indices[shape_func]);
-       if (value == 0.)
-         continue;
- 
-       const Tensor<2,spacedim> *shape_hessians_ptr
-         = &this->shape_hessians[shape_func][0];
-       for (unsigned int point=0; point<n_quadrature_points; ++point)
-         hessians[point] += value * *shape_hessians_ptr++;
+       Vector<double> dof_values(dofs_per_cell);
+       for (unsigned int i=0; i<dofs_per_cell; ++i)
+         dof_values[i] = get_vector_element (fe_function, indices[i]);
+       internal::do_function_derivatives(dof_values.begin(), this->shape_hessians,
+                                         hessians);
      }
  }
  
@@@ -3000,87 -2745,34 +2745,34 @@@ template <int dim, int spacedim
  template <class InputVector, typename number>
  void FEValuesBase<dim,spacedim>::get_function_laplacians (
    const InputVector &fe_function,
 -  const VectorSlice<const std::vector<unsigned int> > &indices,
 +  const VectorSlice<const std::vector<types::global_dof_index> > &indices,
    std::vector<Vector<number> > &laplacians) const
  {
-   // One value per quadrature point
-   Assert (n_quadrature_points == laplacians.size(),
-           ExcDimensionMismatch(laplacians.size(), n_quadrature_points));
- 
-   const unsigned int n_components = fe->n_components();
- 
-   // Size of indices must be a
-   // multiple of dofs_per_cell such
-   // that an integer number of
-   // function values is generated in
-   // each point.
+   // Size of indices must be a multiple of dofs_per_cell such that an integer
+   // number of function values is generated in each point.
    Assert (indices.size() % dofs_per_cell == 0,
            ExcNotMultiple(indices.size(), dofs_per_cell));
- 
-   // The number of components of the
-   // result may be a multiple of the
-   // number of components of the
-   // finite element
-   const unsigned int result_components = indices.size() * n_components / dofs_per_cell;
- 
-   for (unsigned int i=0; i<laplacians.size(); ++i)
-     Assert (laplacians[i].size() == result_components,
-             ExcDimensionMismatch(laplacians[i].size(), result_components));
- 
-   // If the result has more
-   // components than the finite
-   // element, we need this number for
-   // loops filling all components
-   const unsigned int component_multiple = result_components / n_components;
- 
    Assert (this->update_flags & update_hessians, ExcAccessToUninitializedField());
- 
-   // initialize with zero
-   for (unsigned int i=0; i<laplacians.size(); ++i)
-     std::fill_n (laplacians[i].begin(), laplacians[i].size(), 0);
- 
-   // add up contributions of trial
-   // functions. now check whether the
-   // shape function is primitive or
-   // not. if it is, then set its only
-   // non-zero component, otherwise
-   // loop over components
-   for (unsigned int mc = 0; mc < component_multiple; ++mc)
-     for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
-       {
-         const double value = get_vector_element (fe_function, indices[shape_func+mc*dofs_per_cell]);
-         if (value == 0.)
-           continue;
- 
-         if (fe->is_primitive(shape_func))
-           {
-             const unsigned int comp = fe->system_to_component_index(shape_func).first
-                                       + mc * n_components;
-             const unsigned int
-             row = this->shape_function_to_row_table[shape_func * fe->n_components() + comp];
- 
-             const Tensor<2,spacedim> *shape_hessian_ptr
-               = &this->shape_hessians[row][0];
-             for (unsigned int point=0; point<n_quadrature_points; ++point)
-               laplacians[point](comp) += value * trace(*shape_hessian_ptr++);
-           }
-         else
-           for (unsigned int c=0; c<n_components; ++c)
-             {
-               if (fe->get_nonzero_components(shape_func)[c] == false)
-                 continue;
- 
-               const unsigned int
-               row = this->shape_function_to_row_table[shape_func * fe->n_components() + c];
- 
-               const Tensor<2,spacedim> *shape_hessian_ptr
-                 = &this->shape_hessians[row][0];
-               const unsigned int comp = c + mc * n_components;
- 
-               for (unsigned int point=0; point<n_quadrature_points; ++point)
-                 laplacians[point](comp) += value * trace(*shape_hessian_ptr++);
-             }
-       }
+   if (indices.size() <= 100)
+     {
+       double dof_values[100];
+       for (unsigned int i=0; i<indices.size(); ++i)
+         dof_values[i] = get_vector_element (fe_function, indices[i]);
+       internal::do_function_laplacians(&dof_values[0], this->shape_hessians,
+                                        *fe, this->shape_function_to_row_table,
+                                        laplacians, false,
+                                        indices.size()/dofs_per_cell);
+     }
+   else
+     {
+       Vector<double> dof_values(indices.size());
+       for (unsigned int i=0; i<indices.size(); ++i)
+         dof_values[i] = get_vector_element (fe_function, indices[i]);
+       internal::do_function_laplacians(dof_values.begin(),this->shape_hessians,
+                                        *fe, this->shape_function_to_row_table,
+                                        laplacians, false,
+                                        indices.size()/dofs_per_cell);
+     }
  }