From 59056e60b93fab45ce65ae00389b543da39d3b88 Mon Sep 17 00:00:00 2001
From: young <young@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Fri, 14 Feb 2014 13:35:20 +0000
Subject: [PATCH] Blindly patch in fe_values from Denis.

git-svn-id: https://svn.dealii.org/branches/branch_petscscalar_complex@32486 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/include/deal.II/fe/fe_values.h     |  36 +--
 deal.II/source/fe/fe_values.cc             | 327 ++++++++++-----------
 deal.II/source/fe/fe_values.decl.1.inst.in |   4 +-
 deal.II/source/fe/fe_values.decl.2.inst.in |   4 +-
 deal.II/source/fe/fe_values.impl.1.inst.in |  32 +-
 deal.II/source/fe/fe_values.impl.2.inst.in |  28 +-
 deal.II/source/fe/fe_values.inst.in        |  15 +-
 7 files changed, 233 insertions(+), 213 deletions(-)

diff --git a/deal.II/include/deal.II/fe/fe_values.h b/deal.II/include/deal.II/fe/fe_values.h
index 0802018b05..eeb297f369 100644
--- a/deal.II/include/deal.II/fe/fe_values.h
+++ b/deal.II/include/deal.II/fe/fe_values.h
@@ -1709,10 +1709,10 @@ public:
    * Since this function allows for fairly general combinations of argument
    * sizes, be aware that the checks on the arguments may not detect errors.
    */
-  template <class InputVector>
+  template <class InputVector, typename number>
   void get_function_values (const InputVector &fe_function,
                             const VectorSlice<const std::vector<types::global_dof_index> > &indices,
-                            VectorSlice<std::vector<std::vector<double> > > values,
+                            VectorSlice<std::vector<std::vector<number> > > values,
                             const bool quadrature_points_fastest) const;
 
   //@}
@@ -1751,9 +1751,9 @@ public:
    * type IndexSet, then the function is represented as one that is either
    * zero or one, depending on whether a DoF index is in the set or not.
    */
-  template <class InputVector>
+  template <class InputVector, typename number>
   void get_function_gradients (const InputVector      &fe_function,
-                               std::vector<Tensor<1,spacedim> > &gradients) const;
+                               std::vector<Tensor<1,spacedim,number> > &gradients) const;
 
   /**
    * This function does the same as the other get_function_gradients(), but
@@ -1769,27 +1769,27 @@ public:
    * derivative in coordinate direction $d$ of the $c$th vector component of
    * the vector field at quadrature point $q$ of the current cell.
    */
-  template <class InputVector>
+  template <class InputVector, typename number>
   void get_function_gradients (const InputVector               &fe_function,
-                               std::vector<std::vector<Tensor<1,spacedim> > > &gradients) const;
+                               std::vector<std::vector<Tensor<1,spacedim,number> > > &gradients) const;
 
   /**
    * Function gradient access with more flexibility. see get_function_values()
    * with corresponding arguments.
    */
-  template <class InputVector>
+  template <class InputVector, typename number>
   void get_function_gradients (const InputVector &fe_function,
                                const VectorSlice<const std::vector<types::global_dof_index> > &indices,
-                               std::vector<Tensor<1,spacedim> > &gradients) const;
+                               std::vector<Tensor<1,spacedim,number> > &gradients) const;
 
   /**
    * Function gradient access with more flexibility. see get_function_values()
    * with corresponding arguments.
    */
-  template <class InputVector>
+  template <class InputVector, typename number>
   void get_function_gradients (const InputVector &fe_function,
                                const VectorSlice<const std::vector<types::global_dof_index> > &indices,
-                               VectorSlice<std::vector<std::vector<Tensor<1,spacedim> > > > gradients,
+                               VectorSlice<std::vector<std::vector<Tensor<1,spacedim,number> > > > gradients,
                                bool quadrature_points_fastest = false) const;
 
   /**
@@ -1860,10 +1860,10 @@ public:
    * type IndexSet, then the function is represented as one that is either
    * zero or one, depending on whether a DoF index is in the set or not.
    */
-  template <class InputVector>
+  template <class InputVector, typename number>
   void
   get_function_hessians (const InputVector &fe_function,
-                         std::vector<Tensor<2,spacedim> > &hessians) const;
+                         std::vector<Tensor<2,spacedim,number> > &hessians) const;
 
   /**
    * This function does the same as the other get_function_hessians(), but
@@ -1880,31 +1880,31 @@ public:
    * component of the vector field at quadrature point $q$ of the current
    * cell.
    */
-  template <class InputVector>
+  template <class InputVector, typename number>
   void
   get_function_hessians (const InputVector      &fe_function,
-                         std::vector<std::vector<Tensor<2,spacedim> > > &hessians,
+                         std::vector<std::vector<Tensor<2,spacedim,number> > > &hessians,
                          bool quadrature_points_fastest = false) const;
 
   /**
    * Access to the second derivatives of a function with more flexibility. see
    * get_function_values() with corresponding arguments.
    */
-  template <class InputVector>
+  template <class InputVector, typename number>
   void get_function_hessians (
     const InputVector &fe_function,
     const VectorSlice<const std::vector<types::global_dof_index> > &indices,
-    std::vector<Tensor<2,spacedim> > &hessians) const;
+    std::vector<Tensor<2,spacedim,number> > &hessians) const;
 
   /**
    * Access to the second derivatives of a function with more flexibility. see
    * get_function_values() with corresponding arguments.
    */
-  template <class InputVector>
+  template <class InputVector, typename number>
   void get_function_hessians (
     const InputVector &fe_function,
     const VectorSlice<const std::vector<types::global_dof_index> > &indices,
-    VectorSlice<std::vector<std::vector<Tensor<2,spacedim> > > > hessians,
+    VectorSlice<std::vector<std::vector<Tensor<2,spacedim,number> > > > hessians,
     bool quadrature_points_fastest = false) const;
 
   /**
diff --git a/deal.II/source/fe/fe_values.cc b/deal.II/source/fe/fe_values.cc
index 5a99c08bce..f565011ec8 100644
--- a/deal.II/source/fe/fe_values.cc
+++ b/deal.II/source/fe/fe_values.cc
@@ -1,7 +1,7 @@
 // ---------------------------------------------------------------------
 // $Id$
 //
-// Copyright (C) 1998 - 2014 by the deal.II authors
+// Copyright (C) 1998 - 2013 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
@@ -800,22 +800,24 @@ namespace FEValuesViews
                   const dealii::Tensor<1, spacedim> *shape_gradient_ptr =
                     &shape_gradients[snc][0];
 
-                  Assert (shape_function_data[shape_function].single_nonzero_component >= 0,
-                          ExcInternalError());
-                  // we're in 2d, so the formula for the curl is simple:
-                  if (shape_function_data[shape_function].single_nonzero_component_index == 0)
-                    for (unsigned int q_point = 0;
-                        q_point < n_quadrature_points; ++q_point)
-                      curls[q_point][0] -= value * (*shape_gradient_ptr++)[1];
-                  else
-                    for (unsigned int q_point = 0;
-                        q_point < n_quadrature_points; ++q_point)
-                      curls[q_point][0] += value * (*shape_gradient_ptr++)[0];
+                  switch (shape_function_data[shape_function].single_nonzero_component_index)
+                    {
+                    case 0:
+                    {
+                      for (unsigned int q_point = 0;
+                           q_point < n_quadrature_points; ++q_point)
+                        curls[q_point][0] -= value * (*shape_gradient_ptr++)[1];
+
+                      break;
+                    }
+
+                    default:
+                      for (unsigned int q_point = 0;
+                           q_point < n_quadrature_points; ++q_point)
+                        curls[q_point][0] += value * (*shape_gradient_ptr)[0];
+                    }
                 }
               else
-                // we have multiple non-zero components in the shape functions. not
-                // all of them must necessarily be within the 2-component window
-                // this FEValuesViews::Vector object considers, however.
                 {
                   if (shape_function_data[shape_function].is_nonzero_shape_function_component[0])
                     {
@@ -885,26 +887,17 @@ namespace FEValuesViews
                       break;
                     }
 
-                    case 2:
-                    {
+                    default:
                       for (unsigned int q_point = 0;
                            q_point < n_quadrature_points; ++q_point)
                         {
                           curls[q_point][0] += value * (*shape_gradient_ptr)[1];
                           curls[q_point][1] -= value * (*shape_gradient_ptr++)[0];
                         }
-                      break;
-                    }
-
-                    default:
-                      Assert (false, ExcInternalError());
                     }
                 }
 
               else
-                // we have multiple non-zero components in the shape functions. not
-                // all of them must necessarily be within the 3-component window
-                // this FEValuesViews::Vector object considers, however.
                 {
                   if (shape_function_data[shape_function].is_nonzero_shape_function_component[0])
                     {
@@ -1007,7 +1000,7 @@ namespace FEValuesViews
     template <int dim, int spacedim>
     void
     do_function_values (const ::dealii::Vector<double> &dof_values,
-                        const dealii::Table<2,double>          &shape_values,
+                        const Table<2,double>          &shape_values,
                         const std::vector<typename SymmetricTensor<2,dim,spacedim>::ShapeFunctionData> &shape_function_data,
                         std::vector<dealii::SymmetricTensor<2,spacedim> > &values)
     {
@@ -1149,7 +1142,7 @@ namespace FEValuesViews
     template <int dim, int spacedim>
     void
     do_function_values (const ::dealii::Vector<double> &dof_values,
-                        const dealii::Table<2,double>          &shape_values,
+                        const Table<2,double>          &shape_values,
                         const std::vector<typename Tensor<2,dim,spacedim>::ShapeFunctionData> &shape_function_data,
                         std::vector<dealii::Tensor<2,spacedim> > &values)
     {
@@ -1768,10 +1761,10 @@ public:
   /// @p get_interpolated_dof_values
   /// of the iterator with the
   /// given arguments.
-  virtual
-  void
-  get_interpolated_dof_values (const IndexSet &in,
-                               Vector<double> &out) const = 0;
+//  virtual
+//  void
+//  get_interpolated_dof_values (const IndexSet &in,
+//                               Vector<double> &out) const = 0;
 };
 
 
@@ -1837,10 +1830,10 @@ public:
   /// @p get_interpolated_dof_values
   /// of the iterator with the
   /// given arguments.
-  virtual
-  void
-  get_interpolated_dof_values (const IndexSet &in,
-                               Vector<double> &out) const;
+//  virtual
+//  void
+//  get_interpolated_dof_values (const IndexSet &in,
+//                               Vector<double> &out) const;
 
 private:
   /**
@@ -1940,10 +1933,10 @@ public:
   /// @p get_interpolated_dof_values
   /// of the iterator with the
   /// given arguments.
-  virtual
-  void
-  get_interpolated_dof_values (const IndexSet &in,
-                               Vector<double> &out) const;
+//  virtual
+//  void
+//  get_interpolated_dof_values (const IndexSet &in,
+//                               Vector<double> &out) const;
 
 private:
   /**
@@ -2002,21 +1995,21 @@ FEValuesBase<dim,spacedim>::CellIterator<CI>::n_dofs_for_dof_handler () const
 #include "fe_values.impl.1.inst"
 
 
-template <int dim, int spacedim>
-template <typename CI>
-void
-FEValuesBase<dim,spacedim>::CellIterator<CI>::
-get_interpolated_dof_values (const IndexSet &in,
-                             Vector<double> &out) const
-{
-  Assert (cell->has_children() == false, ExcNotImplemented());
-
-  std::vector<types::global_dof_index> dof_indices (cell->get_fe().dofs_per_cell);
-  cell->get_dof_indices (dof_indices);
-
-  for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
-    out[i] = (in.is_element (dof_indices[i]) ? 1 : 0);
-}
+//template <int dim, int spacedim>
+//template <typename CI>
+//void
+//FEValuesBase<dim,spacedim>::CellIterator<CI>::
+//get_interpolated_dof_values (const IndexSet &in,
+//                             Vector<double> &out) const
+//{
+//  Assert (cell->has_children() == false, ExcNotImplemented());
+//
+//  std::vector<types::global_dof_index> dof_indices (cell->get_fe().dofs_per_cell);
+//  cell->get_dof_indices (dof_indices);
+//
+//  for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
+//    out[i] = (in.is_element (dof_indices[i]) ? 1 : 0);
+//}
 
 
 /* ---------------- FEValuesBase<dim,spacedim>::TriaCellIterator --------- */
@@ -2063,14 +2056,14 @@ FEValuesBase<dim,spacedim>::TriaCellIterator::n_dofs_for_dof_handler () const
 #include "fe_values.impl.2.inst"
 
 
-template <int dim, int spacedim>
-void
-FEValuesBase<dim,spacedim>::TriaCellIterator::
-get_interpolated_dof_values (const IndexSet &,
-                             Vector<double> &) const
-{
-  Assert (false, ExcMessage (message_string));
-}
+//template <int dim, int spacedim>
+//void
+//FEValuesBase<dim,spacedim>::TriaCellIterator::
+//get_interpolated_dof_values (const IndexSet &,
+//                             Vector<double> &) const
+//{
+//  Assert (false, ExcMessage (message_string));
+//}
 
 
 
@@ -2205,10 +2198,10 @@ namespace internal
   // compilation and reduces the size of the final file since all the
   // different global vectors get channeled through the same code.
 
-  template <typename Number>
+  template <typename Number, typename Number2>
   void
-  do_function_values (const double          *dof_values_ptr,
-                      const dealii::Table<2,double> &shape_values,
+  do_function_values (const Number2         *dof_values_ptr,
+                      const Table<2,double> &shape_values,
                       std::vector<Number>   &values)
   {
     // scalar finite elements, so shape_values.size() == dofs_per_cell
@@ -2229,8 +2222,8 @@ namespace internal
     // the shape_values data stored contiguously
     for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
       {
-        const double value = dof_values_ptr[shape_func];
-        if (value == 0.)
+        const Number2 value = dof_values_ptr[shape_func];
+        if (value == Number2())
           continue;
 
         const double *shape_value_ptr = &shape_values(shape_func, 0);
@@ -2239,28 +2232,19 @@ namespace internal
       }
   }
 
-  template <int dim, int spacedim, typename VectorType>
+  template <int dim, int spacedim, typename VectorType, typename Number>
   void
-  do_function_values (const double                      *dof_values_ptr,
-                      const dealii::Table<2,double>             &shape_values,
+  do_function_values (const Number                      *dof_values_ptr,
+                      const Table<2,double>             &shape_values,
                       const FiniteElement<dim,spacedim> &fe,
                       const std::vector<unsigned int> &shape_function_to_row_table,
                       VectorSlice<std::vector<VectorType> > &values,
                       const bool quadrature_points_fastest  = false,
                       const unsigned int component_multiple = 1)
   {
-    // initialize with zero
-    for (unsigned int i=0; i<values.size(); ++i)
-      std::fill_n (values[i].begin(), values[i].size(),
-                   typename VectorType::value_type());
-
-    // see if there the current cell has DoFs at all, and if not
-    // then there is nothing else to do.
     const unsigned int dofs_per_cell = fe.dofs_per_cell;
-    if (dofs_per_cell == 0)
-      return;
-    
-    const unsigned int n_quadrature_points = shape_values.n_cols();
+    const unsigned int n_quadrature_points = dofs_per_cell > 0 ?
+                                             shape_values.n_cols() : 0;
     const unsigned int n_components = fe.n_components();
 
     // Assert that we can write all components into the result vectors
@@ -2278,14 +2262,19 @@ namespace internal
           AssertDimension (values[i].size(), result_components);
       }
 
+    // initialize with zero
+    for (unsigned int i=0; i<values.size(); ++i)
+      std::fill_n (values[i].begin(), values[i].size(),
+                   typename VectorType::value_type());
+
     // 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 = dof_values_ptr[shape_func+mc*dofs_per_cell];
-          if (value == 0.)
+          const Number value = dof_values_ptr[shape_func+mc*dofs_per_cell];
+          if (value == Number())
             continue;
 
           if (fe.is_primitive(shape_func))
@@ -2336,11 +2325,11 @@ namespace internal
 
   // use the same implementation for gradients and Hessians, distinguish them
   // by the rank of the tensors
-  template <int order, int spacedim>
+  template <int order, int spacedim, typename Number>
   void
-  do_function_derivatives (const double                     *dof_values_ptr,
+  do_function_derivatives (const Number                     *dof_values_ptr,
                            const std::vector<std::vector<Tensor<order,spacedim> > > &shape_derivatives,
-                           std::vector<Tensor<order,spacedim> > &derivatives)
+                           std::vector<Tensor<order,spacedim,Number> > &derivatives)
   {
     const unsigned int dofs_per_cell = shape_derivatives.size();
     const unsigned int n_quadrature_points = dofs_per_cell > 0 ?
@@ -2348,7 +2337,7 @@ namespace internal
     AssertDimension(derivatives.size(), n_quadrature_points);
 
     // initialize with zero
-    std::fill_n (derivatives.begin(), n_quadrature_points, Tensor<order,spacedim>());
+    std::fill_n (derivatives.begin(), n_quadrature_points, Tensor<order,spacedim,Number>());
 
     // add up contributions of trial functions. note that here we deal with
     // scalar finite elements, so no need to check for non-primitivity of
@@ -2359,8 +2348,8 @@ namespace internal
     // access the shape_gradients/hessians data stored contiguously
     for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
       {
-        const double value = dof_values_ptr[shape_func];
-        if (value == 0.)
+        const Number value = dof_values_ptr[shape_func];
+        if (value == Number())
           continue;
 
         const Tensor<order,spacedim> *shape_derivative_ptr
@@ -2370,29 +2359,19 @@ namespace internal
       }
   }
 
-  template <int order, int dim, int spacedim>
+  template <int order, int dim, int spacedim, typename Number>
   void
-  do_function_derivatives (const double                      *dof_values_ptr,
+  do_function_derivatives (const Number                      *dof_values_ptr,
                            const std::vector<std::vector<Tensor<order,spacedim> > > &shape_derivatives,
                            const FiniteElement<dim,spacedim> &fe,
                            const std::vector<unsigned int> &shape_function_to_row_table,
-                           VectorSlice<std::vector<std::vector<Tensor<order,spacedim> > > > &derivatives,
+                           VectorSlice<std::vector<std::vector<Tensor<order,spacedim,Number> > > > &derivatives,
                            const bool quadrature_points_fastest  = false,
                            const unsigned int component_multiple = 1)
   {
-    // initialize with zero
-    for (unsigned int i=0; i<derivatives.size(); ++i)
-      std::fill_n (derivatives[i].begin(), derivatives[i].size(),
-                   Tensor<order,spacedim>());
-
-    // see if there the current cell has DoFs at all, and if not
-    // then there is nothing else to do.
     const unsigned int dofs_per_cell = fe.dofs_per_cell;
-    if (dofs_per_cell == 0)
-      return;
-
-
-    const unsigned int n_quadrature_points = shape_derivatives[0].size();
+    const unsigned int n_quadrature_points = dofs_per_cell > 0 ?
+                                             shape_derivatives[0].size() : 0;
     const unsigned int n_components = fe.n_components();
 
     // Assert that we can write all components into the result vectors
@@ -2410,14 +2389,19 @@ namespace internal
           AssertDimension (derivatives[i].size(), result_components);
       }
 
+    // initialize with zero
+    for (unsigned int i=0; i<derivatives.size(); ++i)
+      std::fill_n (derivatives[i].begin(), derivatives[i].size(),
+                   Tensor<order,spacedim,Number>());
+
     // 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 = dof_values_ptr[shape_func+mc*dofs_per_cell];
-          if (value == 0.)
+          const Number value = dof_values_ptr[shape_func+mc*dofs_per_cell];
+          if (value == Number())
             continue;
 
           if (fe.is_primitive(shape_func))
@@ -2461,9 +2445,9 @@ namespace internal
         }
   }
 
-  template <int spacedim, typename Number>
+  template <int spacedim, typename Number, typename Number2>
   void
-  do_function_laplacians (const double        *dof_values_ptr,
+  do_function_laplacians (const Number2        *dof_values_ptr,
                           const std::vector<std::vector<Tensor<2,spacedim> > > &shape_hessians,
                           std::vector<Number> &laplacians)
   {
@@ -2480,8 +2464,8 @@ namespace internal
     // the trace of the Hessian.
     for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
       {
-        const double value = dof_values_ptr[shape_func];
-        if (value == 0.)
+        const Number2 value = dof_values_ptr[shape_func];
+        if (value == Number2())
           continue;
 
         const Tensor<2,spacedim> *shape_hessian_ptr
@@ -2491,9 +2475,9 @@ namespace internal
       }
   }
 
-  template <int dim, int spacedim, typename VectorType>
+  template <int dim, int spacedim, typename VectorType, typename Number>
   void
-  do_function_laplacians (const double                    *dof_values_ptr,
+  do_function_laplacians (const Number                    *dof_values_ptr,
                           const std::vector<std::vector<Tensor<2,spacedim> > > &shape_hessians,
                           const FiniteElement<dim,spacedim> &fe,
                           const std::vector<unsigned int> &shape_function_to_row_table,
@@ -2501,19 +2485,9 @@ namespace internal
                           const bool quadrature_points_fastest  = false,
                           const unsigned int component_multiple = 1)
   {
-    // initialize with zero
-    for (unsigned int i=0; i<laplacians.size(); ++i)
-      std::fill_n (laplacians[i].begin(), laplacians[i].size(),
-                   typename VectorType::value_type());
-
-    // see if there the current cell has DoFs at all, and if not
-    // then there is nothing else to do.
     const unsigned int dofs_per_cell = fe.dofs_per_cell;
-    if (dofs_per_cell == 0)
-      return;
-
-
-    const unsigned int n_quadrature_points = shape_hessians[0].size();
+    const unsigned int n_quadrature_points = dofs_per_cell > 0 ?
+                                             shape_hessians[0].size() : 0;
     const unsigned int n_components = fe.n_components();
 
     // Assert that we can write all components into the result vectors
@@ -2531,14 +2505,19 @@ namespace internal
           AssertDimension (laplacians[i].size(), result_components);
       }
 
+    // initialize with zero
+    for (unsigned int i=0; i<laplacians.size(); ++i)
+      std::fill_n (laplacians[i].begin(), laplacians[i].size(),
+                   typename VectorType::value_type());
+
     // 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 = dof_values_ptr[shape_func+mc*dofs_per_cell];
-          if (value == 0.)
+          const Number value = dof_values_ptr[shape_func+mc*dofs_per_cell];
+          if (value == Number())
             continue;
 
           if (fe.is_primitive(shape_func))
@@ -2606,7 +2585,7 @@ void FEValuesBase<dim,spacedim>::get_function_values (
                    present_cell->n_dofs_for_dof_handler());
 
   // get function values of dofs on this cell
-  Vector<double> dof_values (dofs_per_cell);
+  Vector<number> dof_values (dofs_per_cell);
   present_cell->get_interpolated_dof_values(fe_function, dof_values);
   internal::do_function_values (dof_values.begin(), this->shape_values,
                                 values);
@@ -2629,14 +2608,14 @@ void FEValuesBase<dim,spacedim>::get_function_values (
   // avoid allocation when the local size is small enough
   if (dofs_per_cell <= 100)
     {
-      double dof_values[100];
+	  number 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, values);
     }
   else
     {
-      Vector<double> dof_values(dofs_per_cell);
+      Vector<number> 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_values(dof_values.begin(), this->shape_values,
@@ -2660,7 +2639,7 @@ void FEValuesBase<dim,spacedim>::get_function_values (
   AssertDimension (fe_function.size(), present_cell->n_dofs_for_dof_handler());
 
   // get function values of dofs on this cell
-  Vector<double> dof_values (dofs_per_cell);
+  Vector<number> dof_values (dofs_per_cell);
   present_cell->get_interpolated_dof_values(fe_function, dof_values);
   VectorSlice<std::vector<Vector<number> > > val(values);
   internal::do_function_values(dof_values.begin(), this->shape_values, *fe,
@@ -2686,7 +2665,7 @@ void FEValuesBase<dim,spacedim>::get_function_values (
   VectorSlice<std::vector<Vector<number> > > val(values);
   if (indices.size() <= 100)
     {
-      double dof_values[100];
+	  number 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,
@@ -2695,7 +2674,7 @@ void FEValuesBase<dim,spacedim>::get_function_values (
     }
   else
     {
-      Vector<double> dof_values(100);
+      Vector<number> 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,
@@ -2707,11 +2686,11 @@ void FEValuesBase<dim,spacedim>::get_function_values (
 
 
 template <int dim, int spacedim>
-template <class InputVector>
+template <class InputVector, typename number>
 void FEValuesBase<dim,spacedim>::get_function_values (
   const InputVector &fe_function,
   const VectorSlice<const std::vector<types::global_dof_index> > &indices,
-  VectorSlice<std::vector<std::vector<double> > > values,
+  VectorSlice<std::vector<std::vector<number> > > values,
   bool quadrature_points_fastest) const
 {
   Assert (this->update_flags & update_values,
@@ -2724,7 +2703,7 @@ void FEValuesBase<dim,spacedim>::get_function_values (
 
   if (indices.size() <= 100)
     {
-      double dof_values[100];
+	  number 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_values(&dof_values[0], this->shape_values, *fe,
@@ -2734,7 +2713,7 @@ void FEValuesBase<dim,spacedim>::get_function_values (
     }
   else
     {
-      Vector<double> dof_values(indices.size());
+      Vector<number> 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_values(dof_values.begin(), this->shape_values, *fe,
@@ -2747,11 +2726,11 @@ void FEValuesBase<dim,spacedim>::get_function_values (
 
 
 template <int dim, int spacedim>
-template <class InputVector>
+template <class InputVector, typename number>
 void
 FEValuesBase<dim,spacedim>::get_function_gradients (
   const InputVector           &fe_function,
-  std::vector<Tensor<1,spacedim> > &gradients) const
+  std::vector<Tensor<1,spacedim,number> > &gradients) const
 {
   Assert (this->update_flags & update_gradients,
           ExcAccessToUninitializedField("update_gradients"));
@@ -2761,7 +2740,7 @@ FEValuesBase<dim,spacedim>::get_function_gradients (
   AssertDimension (fe_function.size(), present_cell->n_dofs_for_dof_handler());
 
   // get function values of dofs on this cell
-  Vector<double> dof_values (dofs_per_cell);
+  Vector<number> dof_values (dofs_per_cell);
   present_cell->get_interpolated_dof_values(fe_function, dof_values);
   internal::do_function_derivatives(dof_values.begin(), this->shape_gradients,
                                     gradients);
@@ -2770,11 +2749,11 @@ FEValuesBase<dim,spacedim>::get_function_gradients (
 
 
 template <int dim, int spacedim>
-template <class InputVector>
+template <class InputVector, typename number>
 void FEValuesBase<dim,spacedim>::get_function_gradients (
   const InputVector &fe_function,
   const VectorSlice<const std::vector<types::global_dof_index> > &indices,
-  std::vector<Tensor<1,spacedim> > &gradients) const
+  std::vector<Tensor<1,spacedim,number> > &gradients) const
 {
   Assert (this->update_flags & update_gradients,
           ExcAccessToUninitializedField("update_gradients"));
@@ -2782,7 +2761,7 @@ void FEValuesBase<dim,spacedim>::get_function_gradients (
   AssertDimension (indices.size(), dofs_per_cell);
   if (dofs_per_cell <= 100)
     {
-      double dof_values[100];
+      number 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_gradients,
@@ -2790,7 +2769,7 @@ void FEValuesBase<dim,spacedim>::get_function_gradients (
     }
   else
     {
-      Vector<double> dof_values(dofs_per_cell);
+      Vector<number> 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_gradients,
@@ -2802,11 +2781,11 @@ void FEValuesBase<dim,spacedim>::get_function_gradients (
 
 
 template <int dim, int spacedim>
-template <class InputVector>
+template <class InputVector, typename number>
 void
 FEValuesBase<dim,spacedim>::get_function_gradients (
   const InputVector                              &fe_function,
-  std::vector<std::vector<Tensor<1,spacedim> > > &gradients) const
+  std::vector<std::vector<Tensor<1,spacedim,number> > > &gradients) const
 {
   Assert (this->update_flags & update_gradients,
           ExcAccessToUninitializedField("update_gradients"));
@@ -2815,9 +2794,9 @@ FEValuesBase<dim,spacedim>::get_function_gradients (
   AssertDimension (fe_function.size(), present_cell->n_dofs_for_dof_handler());
 
   // get function values of dofs on this cell
-  Vector<double> dof_values (dofs_per_cell);
+  Vector<number> dof_values (dofs_per_cell);
   present_cell->get_interpolated_dof_values(fe_function, dof_values);
-  VectorSlice<std::vector<std::vector<Tensor<1,spacedim> > > > grads(gradients);
+  VectorSlice<std::vector<std::vector<Tensor<1,spacedim,number> > > > grads(gradients);
   internal::do_function_derivatives(dof_values.begin(), this->shape_gradients,
                                     *fe, this->shape_function_to_row_table,
                                     grads);
@@ -2826,11 +2805,11 @@ FEValuesBase<dim,spacedim>::get_function_gradients (
 
 
 template <int dim, int spacedim>
-template <class InputVector>
+template <class InputVector, typename number>
 void FEValuesBase<dim,spacedim>::get_function_gradients (
   const InputVector &fe_function,
   const VectorSlice<const std::vector<types::global_dof_index> > &indices,
-  VectorSlice<std::vector<std::vector<Tensor<1,spacedim> > > > gradients,
+  VectorSlice<std::vector<std::vector<Tensor<1,spacedim,number> > > > gradients,
   bool quadrature_points_fastest) const
 {
   // Size of indices must be a multiple of dofs_per_cell such that an integer
@@ -2842,7 +2821,7 @@ void FEValuesBase<dim,spacedim>::get_function_gradients (
 
   if (indices.size() <= 100)
     {
-      double dof_values[100];
+      number 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_derivatives(&dof_values[0], this->shape_gradients,
@@ -2852,7 +2831,7 @@ void FEValuesBase<dim,spacedim>::get_function_gradients (
     }
   else
     {
-      Vector<double> dof_values(indices.size());
+      Vector<number> 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_derivatives(dof_values.begin(),this->shape_gradients,
@@ -2865,11 +2844,11 @@ void FEValuesBase<dim,spacedim>::get_function_gradients (
 
 
 template <int dim, int spacedim>
-template <class InputVector>
+template <class InputVector, typename number>
 void
 FEValuesBase<dim,spacedim>::
 get_function_hessians (const InputVector                &fe_function,
-                       std::vector<Tensor<2,spacedim> > &hessians) const
+                       std::vector<Tensor<2,spacedim,number> > &hessians) const
 {
   AssertDimension (fe->n_components(), 1);
   Assert (this->update_flags & update_hessians,
@@ -2879,7 +2858,7 @@ get_function_hessians (const InputVector                &fe_function,
   AssertDimension (fe_function.size(), present_cell->n_dofs_for_dof_handler());
 
   // get function values of dofs on this cell
-  Vector<double> dof_values (dofs_per_cell);
+  Vector<number> dof_values (dofs_per_cell);
   present_cell->get_interpolated_dof_values(fe_function, dof_values);
   internal::do_function_derivatives(dof_values.begin(), this->shape_hessians,
                                     hessians);
@@ -2888,11 +2867,11 @@ get_function_hessians (const InputVector                &fe_function,
 
 
 template <int dim, int spacedim>
-template <class InputVector>
+template <class InputVector, typename number>
 void FEValuesBase<dim,spacedim>::get_function_hessians (
   const InputVector &fe_function,
   const VectorSlice<const std::vector<types::global_dof_index> > &indices,
-  std::vector<Tensor<2,spacedim> > &hessians) const
+  std::vector<Tensor<2,spacedim,number> > &hessians) const
 {
   Assert (this->update_flags & update_hessians,
           ExcAccessToUninitializedField("update_hessians"));
@@ -2900,7 +2879,7 @@ void FEValuesBase<dim,spacedim>::get_function_hessians (
   AssertDimension (indices.size(), dofs_per_cell);
   if (dofs_per_cell <= 100)
     {
-      double dof_values[100];
+      number 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,
@@ -2908,7 +2887,7 @@ void FEValuesBase<dim,spacedim>::get_function_hessians (
     }
   else
     {
-      Vector<double> dof_values(dofs_per_cell);
+      Vector<number> 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,
@@ -2920,11 +2899,11 @@ void FEValuesBase<dim,spacedim>::get_function_hessians (
 
 
 template <int dim, int spacedim>
-template <class InputVector>
+template <class InputVector, typename number>
 void
 FEValuesBase<dim,spacedim>::
 get_function_hessians (const InputVector                         &fe_function,
-                       std::vector<std::vector<Tensor<2,spacedim> > > &hessians,
+                       std::vector<std::vector<Tensor<2,spacedim,number> > > &hessians,
                        bool quadrature_points_fastest) const
 {
   Assert (this->update_flags & update_hessians,
@@ -2934,9 +2913,9 @@ get_function_hessians (const InputVector                         &fe_function,
   AssertDimension (fe_function.size(), present_cell->n_dofs_for_dof_handler());
 
   // get function values of dofs on this cell
-  Vector<double> dof_values (dofs_per_cell);
+  Vector<number> dof_values (dofs_per_cell);
   present_cell->get_interpolated_dof_values(fe_function, dof_values);
-  VectorSlice<std::vector<std::vector<Tensor<2,spacedim> > > > hes(hessians);
+  VectorSlice<std::vector<std::vector<Tensor<2,spacedim,number> > > > hes(hessians);
   internal::do_function_derivatives(dof_values.begin(), this->shape_hessians,
                                     *fe, this->shape_function_to_row_table,
                                     hes, quadrature_points_fastest);
@@ -2945,11 +2924,11 @@ get_function_hessians (const InputVector                         &fe_function,
 
 
 template <int dim, int spacedim>
-template <class InputVector>
+template <class InputVector, typename number>
 void FEValuesBase<dim, spacedim>::get_function_hessians (
   const InputVector &fe_function,
   const VectorSlice<const std::vector<types::global_dof_index> > &indices,
-  VectorSlice<std::vector<std::vector<Tensor<2,spacedim> > > > hessians,
+  VectorSlice<std::vector<std::vector<Tensor<2,spacedim,number> > > > hessians,
   bool quadrature_points_fastest) const
 {
   Assert (this->update_flags & update_hessians,
@@ -2958,7 +2937,7 @@ void FEValuesBase<dim, spacedim>::get_function_hessians (
           ExcNotMultiple(indices.size(), dofs_per_cell));
   if (indices.size() <= 100)
     {
-      double dof_values[100];
+      number 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_derivatives(&dof_values[0], this->shape_hessians,
@@ -2968,7 +2947,7 @@ void FEValuesBase<dim, spacedim>::get_function_hessians (
     }
   else
     {
-      Vector<double> dof_values(indices.size());
+      Vector<number> 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_derivatives(dof_values.begin(),this->shape_hessians,
@@ -2994,7 +2973,7 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
   AssertDimension (fe_function.size(), present_cell->n_dofs_for_dof_handler());
 
   // get function values of dofs on this cell
-  Vector<double> dof_values (dofs_per_cell);
+  Vector<number> dof_values (dofs_per_cell);
   present_cell->get_interpolated_dof_values(fe_function, dof_values);
   internal::do_function_laplacians(dof_values.begin(), this->shape_hessians,
                                    laplacians);
@@ -3015,7 +2994,7 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
   AssertDimension (indices.size(), dofs_per_cell);
   if (dofs_per_cell <= 100)
     {
-      double dof_values[100];
+      number 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_laplacians(&dof_values[0], this->shape_hessians,
@@ -3023,7 +3002,7 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
     }
   else
     {
-      Vector<double> dof_values(dofs_per_cell);
+      Vector<number> 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_laplacians(dof_values.begin(), this->shape_hessians,
@@ -3046,7 +3025,7 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
   AssertDimension (fe_function.size(), present_cell->n_dofs_for_dof_handler());
 
   // get function values of dofs on this cell
-  Vector<double> dof_values (dofs_per_cell);
+  Vector<number> dof_values (dofs_per_cell);
   present_cell->get_interpolated_dof_values(fe_function, dof_values);
   internal::do_function_laplacians(dof_values.begin(), this->shape_hessians,
                                    *fe, this->shape_function_to_row_table,
@@ -3070,7 +3049,7 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
           ExcAccessToUninitializedField("update_hessians"));
   if (indices.size() <= 100)
     {
-      double dof_values[100];
+      number 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,
@@ -3080,7 +3059,7 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
     }
   else
     {
-      Vector<double> dof_values(indices.size());
+      Vector<number> 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,
@@ -3106,7 +3085,7 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
           ExcAccessToUninitializedField("update_hessians"));
   if (indices.size() <= 100)
     {
-      double dof_values[100];
+      number 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,
@@ -3116,7 +3095,7 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
     }
   else
     {
-      Vector<double> dof_values(indices.size());
+      Vector<number> 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,
diff --git a/deal.II/source/fe/fe_values.decl.1.inst.in b/deal.II/source/fe/fe_values.decl.1.inst.in
index 6cbab13c8f..99a4856d08 100644
--- a/deal.II/source/fe/fe_values.decl.1.inst.in
+++ b/deal.II/source/fe/fe_values.decl.1.inst.in
@@ -19,7 +19,7 @@
 // Declarations of member functions of FEValuesBase::CellIteratorBase
 // and derived classes
 
-for (VEC : REAL_SERIAL_VECTORS)
+for (VEC : SERIAL_VECTORS_REAL; S: REAL_SCALARS)
   {
 				      /// Call
 				      /// @p get_interpolated_dof_values
@@ -28,7 +28,7 @@ for (VEC : REAL_SERIAL_VECTORS)
     virtual
     void
     get_interpolated_dof_values (const VEC &in,
-				 Vector<double> &out) const = 0;
+				 Vector<S> &out) const = 0;
   }
   
 for (VEC : SERIAL_VECTORS_COMPLEX)
diff --git a/deal.II/source/fe/fe_values.decl.2.inst.in b/deal.II/source/fe/fe_values.decl.2.inst.in
index f676fad32a..6a26e8e72f 100644
--- a/deal.II/source/fe/fe_values.decl.2.inst.in
+++ b/deal.II/source/fe/fe_values.decl.2.inst.in
@@ -19,7 +19,7 @@
 // Declarations of member functions of FEValuesBase::CellIteratorBase
 // and derived classes
 
-for (VEC : REAL_SERIAL_VECTORS)
+for (VEC : SERIAL_VECTORS_REAL; S: REAL_SCALARS)
   {
 				      /// Call
 				      /// @p get_interpolated_dof_values
@@ -28,7 +28,7 @@ for (VEC : REAL_SERIAL_VECTORS)
     virtual
     void
     get_interpolated_dof_values (const VEC      &in,
-				 Vector<double> &out) const;
+				 Vector<S> &out) const;
   }  
   
 for (VEC : SERIAL_VECTORS_COMPLEX)
diff --git a/deal.II/source/fe/fe_values.impl.1.inst.in b/deal.II/source/fe/fe_values.impl.1.inst.in
index 731d6ccb5a..96d4ee678b 100644
--- a/deal.II/source/fe/fe_values.impl.1.inst.in
+++ b/deal.II/source/fe/fe_values.impl.1.inst.in
@@ -15,29 +15,49 @@
 // ---------------------------------------------------------------------
 
 
-for (VEC : REAL_SERIAL_VECTORS)
+for (VEC : SERIAL_VECTORS_REAL; S: REAL_SCALARS)
   {
     template <int dim, int spacedim>
     template <typename CI>
     void
     FEValuesBase<dim,spacedim>::CellIterator<CI>::
     get_interpolated_dof_values (const VEC      &in,
-			         Vector<double> &out) const
+			         Vector<S> &out) const
     \{
       cell->get_interpolated_dof_values (in, out);
     \}
   }
-
-for (VEC : COMPLEX_SERIAL_VECTORS)
+  
+for (VEC : SERIAL_VECTORS_COMPLEX)
   {
     template <int dim, int spacedim>
     template <typename CI>
     void
     FEValuesBase<dim,spacedim>::CellIterator<CI>::
-    get_interpolated_dof_values (const VEC                     &in,
+    get_interpolated_dof_values (const VEC      &in,
 			         Vector<std::complex<double> > &out) const
     \{
       cell->get_interpolated_dof_values (in, out);
     \}
-  }
+  }  
+
+  
+ for (S: SCALARS)
+  {
+  
+  template <int dim, int spacedim>
+  template <typename CI>
+  void
+  FEValuesBase<dim,spacedim>::CellIterator<CI>::
+  get_interpolated_dof_values (const IndexSet &in,
+                               Vector<S> &out) const
+  \{
+    Assert (cell->has_children() == false, ExcNotImplemented());
+
+    std::vector<types::global_dof_index> dof_indices (cell->get_fe().dofs_per_cell);
+    cell->get_dof_indices (dof_indices);
 
+    for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
+      out[i] = (in.is_element (dof_indices[i]) ? 1 : 0);
+  \}
+  } 
\ No newline at end of file
diff --git a/deal.II/source/fe/fe_values.impl.2.inst.in b/deal.II/source/fe/fe_values.impl.2.inst.in
index 4209718bf6..368c147d48 100644
--- a/deal.II/source/fe/fe_values.impl.2.inst.in
+++ b/deal.II/source/fe/fe_values.impl.2.inst.in
@@ -15,28 +15,38 @@
 // ---------------------------------------------------------------------
 
 
-for (VEC : REAL_SERIAL_VECTORS)
+for (VEC : SERIAL_VECTORS_REAL; S: REAL_SCALARS)
   {
     template <int dim, int spacedim>
     void
     FEValuesBase<dim,spacedim>::TriaCellIterator::
-    get_interpolated_dof_values (const VEC      &,
-    			         Vector<double> &) const
+    get_interpolated_dof_values (const VEC &,
+    			         Vector<S> &) const
     \{
       Assert (false, ExcMessage (message_string));
     \}
   }
-
-
-for (VEC : COMPLEX_SERIAL_VECTORS)
+  
+for (VEC : SERIAL_VECTORS_COMPLEX)
   {
     template <int dim, int spacedim>
     void
     FEValuesBase<dim,spacedim>::TriaCellIterator::
-    get_interpolated_dof_values (const VEC                     &,
+    get_interpolated_dof_values (const VEC &,
     			         Vector<std::complex<double> > &) const
     \{
       Assert (false, ExcMessage (message_string));
     \}
-  }
-
+  }  
+  
+for (S: SCALARS)
+  {
+   template <int dim, int spacedim>
+   void
+   FEValuesBase<dim,spacedim>::TriaCellIterator::
+   get_interpolated_dof_values (const IndexSet &,
+                                Vector<S> &) const
+   \{
+      Assert (false, ExcMessage (message_string));
+   \}
+  }  
diff --git a/deal.II/source/fe/fe_values.inst.in b/deal.II/source/fe/fe_values.inst.in
index d0201ca59d..be28a92084 100644
--- a/deal.II/source/fe/fe_values.inst.in
+++ b/deal.II/source/fe/fe_values.inst.in
@@ -64,7 +64,7 @@ for (dof_handler : DOFHANDLER_TEMPLATES; deal_II_dimension : DIMENSIONS; deal_II
 
 
 
-for (VEC : REAL_SERIAL_VECTORS; deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
+for (VEC : SERIAL_VECTORS_REAL; deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
   {
 #if deal_II_dimension <= deal_II_space_dimension
 #if (deal_II_space_dimension == DIM_A) || (deal_II_space_dimension == DIM_B)
@@ -137,7 +137,7 @@ for (VEC : REAL_SERIAL_VECTORS; deal_II_dimension : DIMENSIONS; deal_II_space_di
 
 
 
-for (VEC : REAL_SERIAL_VECTORS; deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
+for (VEC : SERIAL_VECTORS_REAL; deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
   {
 #if deal_II_dimension <= deal_II_space_dimension
 #if (deal_II_space_dimension == DIM_A) || (deal_II_space_dimension == DIM_B)
@@ -254,6 +254,17 @@ for (VEC : REAL_SERIAL_VECTORS; deal_II_dimension : DIMENSIONS; deal_II_space_di
 #endif
   }
 
+for (VEC : SERIAL_VECTORS_COMPLEX; deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
+  {
+#if deal_II_dimension <= deal_II_space_dimension
+#if (deal_II_space_dimension == DIM_A) || (deal_II_space_dimension == DIM_B)
+
+//    template
+//      void FEValuesBase<deal_II_dimension,deal_II_space_dimension>::get_function_values<VEC>
+//      (const VEC&, std::vector<std::complex<double> >&) const;
+#endif
+#endif      
+  }
 
 // instantiations for VEC=IndexSet
 
-- 
2.39.5