From bcbd9e9f96830fc9ef8cd3a0ae3415f5f79f9cf5 Mon Sep 17 00:00:00 2001
From: David Wells <drwells@email.unc.edu>
Date: Fri, 2 Jun 2023 10:48:10 -0400
Subject: [PATCH] Make all higher derivative functions of FEValues use
 ReadVector.

---
 include/deal.II/fe/fe_values.h | 182 +++++++++-----------
 source/fe/fe_values.cc         | 303 ++++++++++++---------------------
 source/fe/fe_values.inst.in    | 165 ++++++++----------
 3 files changed, 263 insertions(+), 387 deletions(-)

diff --git a/include/deal.II/fe/fe_values.h b/include/deal.II/fe/fe_values.h
index 218b391153..b34b0c5e70 100644
--- a/include/deal.II/fe/fe_values.h
+++ b/include/deal.II/fe/fe_values.h
@@ -520,12 +520,11 @@ namespace FEValuesViews
      *
      * @dealiiRequiresUpdateFlags{update_hessians}
      */
-    template <class InputVector>
+    template <typename Number>
     void
     get_function_hessians(
-      const InputVector &fe_function,
-      std::vector<solution_hessian_type<typename InputVector::value_type>>
-        &hessians) const;
+      const ReadVector<Number> &                  fe_function,
+      std::vector<solution_hessian_type<Number>> &hessians) const;
 
     /**
      * This function relates to get_function_hessians() in the same way
@@ -559,12 +558,11 @@ namespace FEValuesViews
      *
      * @dealiiRequiresUpdateFlags{update_hessians}
      */
-    template <class InputVector>
+    template <typename Number>
     void
     get_function_laplacians(
-      const InputVector &fe_function,
-      std::vector<solution_laplacian_type<typename InputVector::value_type>>
-        &laplacians) const;
+      const ReadVector<Number> &                    fe_function,
+      std::vector<solution_laplacian_type<Number>> &laplacians) const;
 
     /**
      * This function relates to get_function_laplacians() in the same way
@@ -598,13 +596,12 @@ namespace FEValuesViews
      *
      * @dealiiRequiresUpdateFlags{update_third_derivatives}
      */
-    template <class InputVector>
+    template <typename Number>
     void
     get_function_third_derivatives(
-      const InputVector &fe_function,
-      std::vector<
-        solution_third_derivative_type<typename InputVector::value_type>>
-        &third_derivatives) const;
+      const ReadVector<Number> &                           fe_function,
+      std::vector<solution_third_derivative_type<Number>> &third_derivatives)
+      const;
 
     /**
      * This function relates to get_function_third_derivatives() in the same way
@@ -1212,12 +1209,11 @@ namespace FEValuesViews
      *
      * @dealiiRequiresUpdateFlags{update_gradients}
      */
-    template <class InputVector>
+    template <typename Number>
     void
     get_function_symmetric_gradients(
-      const InputVector &fe_function,
-      std::vector<
-        solution_symmetric_gradient_type<typename InputVector::value_type>>
+      const ReadVector<Number> &fe_function,
+      std::vector<solution_symmetric_gradient_type<Number>>
         &symmetric_gradients) const;
 
     /**
@@ -1252,12 +1248,11 @@ namespace FEValuesViews
      *
      * @dealiiRequiresUpdateFlags{update_gradients}
      */
-    template <class InputVector>
+    template <typename Number>
     void
     get_function_divergences(
-      const InputVector &fe_function,
-      std::vector<solution_divergence_type<typename InputVector::value_type>>
-        &divergences) const;
+      const ReadVector<Number> &                     fe_function,
+      std::vector<solution_divergence_type<Number>> &divergences) const;
 
     /**
      * This function relates to get_function_divergences() in the same way
@@ -1290,12 +1285,10 @@ namespace FEValuesViews
      *
      * @dealiiRequiresUpdateFlags{update_gradients}
      */
-    template <class InputVector>
+    template <typename Number>
     void
-    get_function_curls(
-      const InputVector &fe_function,
-      std::vector<solution_curl_type<typename InputVector::value_type>> &curls)
-      const;
+    get_function_curls(const ReadVector<Number> &               fe_function,
+                       std::vector<solution_curl_type<Number>> &curls) const;
 
     /**
      * This function relates to get_function_curls() in the same way
@@ -1327,12 +1320,11 @@ namespace FEValuesViews
      *
      * @dealiiRequiresUpdateFlags{update_hessians}
      */
-    template <class InputVector>
+    template <typename Number>
     void
     get_function_hessians(
-      const InputVector &fe_function,
-      std::vector<solution_hessian_type<typename InputVector::value_type>>
-        &hessians) const;
+      const ReadVector<Number> &                  fe_function,
+      std::vector<solution_hessian_type<Number>> &hessians) const;
 
     /**
      * This function relates to get_function_hessians() in the same way
@@ -1365,12 +1357,11 @@ namespace FEValuesViews
      *
      * @dealiiRequiresUpdateFlags{update_hessians}
      */
-    template <class InputVector>
+    template <typename Number>
     void
     get_function_laplacians(
-      const InputVector &fe_function,
-      std::vector<solution_laplacian_type<typename InputVector::value_type>>
-        &laplacians) const;
+      const ReadVector<Number> &                    fe_function,
+      std::vector<solution_laplacian_type<Number>> &laplacians) const;
 
     /**
      * This function relates to get_function_laplacians() in the same way
@@ -1403,13 +1394,12 @@ namespace FEValuesViews
      *
      * @dealiiRequiresUpdateFlags{update_third_derivatives}
      */
-    template <class InputVector>
+    template <typename Number>
     void
     get_function_third_derivatives(
-      const InputVector &fe_function,
-      std::vector<
-        solution_third_derivative_type<typename InputVector::value_type>>
-        &third_derivatives) const;
+      const ReadVector<Number> &                           fe_function,
+      std::vector<solution_third_derivative_type<Number>> &third_derivatives)
+      const;
 
     /**
      * This function relates to get_function_third_derivatives() in the same way
@@ -1745,12 +1735,11 @@ namespace FEValuesViews
      *
      * @dealiiRequiresUpdateFlags{update_gradients}
      */
-    template <class InputVector>
+    template <typename Number>
     void
     get_function_divergences(
-      const InputVector &fe_function,
-      std::vector<solution_divergence_type<typename InputVector::value_type>>
-        &divergences) const;
+      const ReadVector<Number> &                     fe_function,
+      std::vector<solution_divergence_type<Number>> &divergences) const;
 
     /**
      * This function relates to get_function_divergences() in the same way
@@ -2119,12 +2108,11 @@ namespace FEValuesViews
      *
      * @dealiiRequiresUpdateFlags{update_gradients}
      */
-    template <class InputVector>
+    template <typename Number>
     void
     get_function_divergences(
-      const InputVector &fe_function,
-      std::vector<solution_divergence_type<typename InputVector::value_type>>
-        &divergences) const;
+      const ReadVector<Number> &                     fe_function,
+      std::vector<solution_divergence_type<Number>> &divergences) const;
 
     /**
      * This function relates to get_function_divergences() in the same way
@@ -2986,12 +2974,11 @@ public:
    *
    * @dealiiRequiresUpdateFlags{update_hessians}
    */
-  template <class InputVector>
+  template <typename Number>
   void
   get_function_hessians(
-    const InputVector &fe_function,
-    std::vector<Tensor<2, spacedim, typename InputVector::value_type>>
-      &hessians) const;
+    const ReadVector<Number> &                fe_function,
+    std::vector<Tensor<2, spacedim, Number>> &hessians) const;
 
   /**
    * This function does the same as the other get_function_hessians(), but
@@ -3010,13 +2997,11 @@ public:
    *
    * @dealiiRequiresUpdateFlags{update_hessians}
    */
-  template <class InputVector>
+  template <typename Number>
   void
   get_function_hessians(
-    const InputVector &fe_function,
-    std::vector<
-      std::vector<Tensor<2, spacedim, typename InputVector::value_type>>>
-      &        hessians,
+    const ReadVector<Number> &                             fe_function,
+    std::vector<std::vector<Tensor<2, spacedim, Number>>> &hessians,
     const bool quadrature_points_fastest = false) const;
 
   /**
@@ -3027,13 +3012,12 @@ public:
    *
    * @dealiiRequiresUpdateFlags{update_hessians}
    */
-  template <class InputVector>
+  template <typename Number>
   void
   get_function_hessians(
-    const InputVector &                             fe_function,
+    const ReadVector<Number> &                      fe_function,
     const ArrayView<const types::global_dof_index> &indices,
-    std::vector<Tensor<2, spacedim, typename InputVector::value_type>>
-      &hessians) const;
+    std::vector<Tensor<2, spacedim, Number>> &      hessians) const;
 
   /**
    * This function relates to the first of the get_function_hessians() function
@@ -3043,14 +3027,12 @@ public:
    *
    * @dealiiRequiresUpdateFlags{update_hessians}
    */
-  template <class InputVector>
+  template <typename Number>
   void
   get_function_hessians(
-    const InputVector &                             fe_function,
-    const ArrayView<const types::global_dof_index> &indices,
-    ArrayView<
-      std::vector<Tensor<2, spacedim, typename InputVector::value_type>>>
-               hessians,
+    const ReadVector<Number> &                          fe_function,
+    const ArrayView<const types::global_dof_index> &    indices,
+    ArrayView<std::vector<Tensor<2, spacedim, Number>>> hessians,
     const bool quadrature_points_fastest = false) const;
 
   /**
@@ -3093,11 +3075,10 @@ public:
    *
    * @dealiiRequiresUpdateFlags{update_hessians}
    */
-  template <class InputVector>
+  template <typename Number>
   void
-  get_function_laplacians(
-    const InputVector &                            fe_function,
-    std::vector<typename InputVector::value_type> &laplacians) const;
+  get_function_laplacians(const ReadVector<Number> &fe_function,
+                          std::vector<Number> &     laplacians) const;
 
   /**
    * This function does the same as the other get_function_laplacians(), but
@@ -3118,11 +3099,10 @@ public:
    *
    * @dealiiRequiresUpdateFlags{update_hessians}
    */
-  template <class InputVector>
+  template <typename Number>
   void
-  get_function_laplacians(
-    const InputVector &                                    fe_function,
-    std::vector<Vector<typename InputVector::value_type>> &laplacians) const;
+  get_function_laplacians(const ReadVector<Number> &   fe_function,
+                          std::vector<Vector<Number>> &laplacians) const;
 
   /**
    * This function relates to the first of the get_function_laplacians()
@@ -3132,12 +3112,12 @@ public:
    *
    * @dealiiRequiresUpdateFlags{update_hessians}
    */
-  template <class InputVector>
+  template <typename Number>
   void
   get_function_laplacians(
-    const InputVector &                             fe_function,
+    const ReadVector<Number> &                      fe_function,
     const ArrayView<const types::global_dof_index> &indices,
-    std::vector<typename InputVector::value_type> & laplacians) const;
+    std::vector<Number> &                           laplacians) const;
 
   /**
    * This function relates to the first of the get_function_laplacians()
@@ -3147,12 +3127,12 @@ public:
    *
    * @dealiiRequiresUpdateFlags{update_hessians}
    */
-  template <class InputVector>
+  template <typename Number>
   void
   get_function_laplacians(
-    const InputVector &                                    fe_function,
-    const ArrayView<const types::global_dof_index> &       indices,
-    std::vector<Vector<typename InputVector::value_type>> &laplacians) const;
+    const ReadVector<Number> &                      fe_function,
+    const ArrayView<const types::global_dof_index> &indices,
+    std::vector<Vector<Number>> &                   laplacians) const;
 
   /**
    * This function relates to the first of the get_function_laplacians()
@@ -3162,12 +3142,12 @@ public:
    *
    * @dealiiRequiresUpdateFlags{update_hessians}
    */
-  template <class InputVector>
+  template <typename Number>
   void
   get_function_laplacians(
-    const InputVector &                                         fe_function,
-    const ArrayView<const types::global_dof_index> &            indices,
-    std::vector<std::vector<typename InputVector::value_type>> &laplacians,
+    const ReadVector<Number> &                      fe_function,
+    const ArrayView<const types::global_dof_index> &indices,
+    std::vector<std::vector<Number>> &              laplacians,
     const bool quadrature_points_fastest = false) const;
 
   /** @} */
@@ -3212,12 +3192,11 @@ public:
    *
    * @dealiiRequiresUpdateFlags{update_3rd_derivatives}
    */
-  template <class InputVector>
+  template <typename Number>
   void
   get_function_third_derivatives(
-    const InputVector &fe_function,
-    std::vector<Tensor<3, spacedim, typename InputVector::value_type>>
-      &third_derivatives) const;
+    const ReadVector<Number> &                fe_function,
+    std::vector<Tensor<3, spacedim, Number>> &third_derivatives) const;
 
   /**
    * This function does the same as the other
@@ -3237,13 +3216,11 @@ public:
    *
    * @dealiiRequiresUpdateFlags{update_3rd_derivatives}
    */
-  template <class InputVector>
+  template <typename Number>
   void
   get_function_third_derivatives(
-    const InputVector &fe_function,
-    std::vector<
-      std::vector<Tensor<3, spacedim, typename InputVector::value_type>>>
-      &        third_derivatives,
+    const ReadVector<Number> &                             fe_function,
+    std::vector<std::vector<Tensor<3, spacedim, Number>>> &third_derivatives,
     const bool quadrature_points_fastest = false) const;
 
   /**
@@ -3254,13 +3231,12 @@ public:
    *
    * @dealiiRequiresUpdateFlags{update_3rd_derivatives}
    */
-  template <class InputVector>
+  template <typename Number>
   void
   get_function_third_derivatives(
-    const InputVector &                             fe_function,
+    const ReadVector<Number> &                      fe_function,
     const ArrayView<const types::global_dof_index> &indices,
-    std::vector<Tensor<3, spacedim, typename InputVector::value_type>>
-      &third_derivatives) const;
+    std::vector<Tensor<3, spacedim, Number>> &      third_derivatives) const;
 
   /**
    * This function relates to the first of the get_function_third_derivatives()
@@ -3270,14 +3246,12 @@ public:
    *
    * @dealiiRequiresUpdateFlags{update_3rd_derivatives}
    */
-  template <class InputVector>
+  template <typename Number>
   void
   get_function_third_derivatives(
-    const InputVector &                             fe_function,
-    const ArrayView<const types::global_dof_index> &indices,
-    ArrayView<
-      std::vector<Tensor<3, spacedim, typename InputVector::value_type>>>
-               third_derivatives,
+    const ReadVector<Number> &                          fe_function,
+    const ArrayView<const types::global_dof_index> &    indices,
+    ArrayView<std::vector<Tensor<3, spacedim, Number>>> third_derivatives,
     const bool quadrature_points_fastest = false) const;
   /** @} */
 
diff --git a/source/fe/fe_values.cc b/source/fe/fe_values.cc
index 7df91e0af0..ffe9371e6b 100644
--- a/source/fe/fe_values.cc
+++ b/source/fe/fe_values.cc
@@ -32,18 +32,7 @@
 #include <deal.II/grid/tria_accessor.h>
 #include <deal.II/grid/tria_iterator.h>
 
-#include <deal.II/lac/block_vector.h>
-#include <deal.II/lac/la_parallel_block_vector.h>
-#include <deal.II/lac/la_parallel_vector.h>
-#include <deal.II/lac/la_vector.h>
-#include <deal.II/lac/petsc_block_vector.h>
-#include <deal.II/lac/petsc_vector.h>
-#include <deal.II/lac/trilinos_epetra_vector.h>
-#include <deal.II/lac/trilinos_parallel_block_vector.h>
-#include <deal.II/lac/trilinos_tpetra_vector.h>
-#include <deal.II/lac/trilinos_vector.h>
 #include <deal.II/lac/vector.h>
-#include <deal.II/lac/vector_element_access.h>
 
 #include <boost/container/small_vector.hpp>
 
@@ -56,25 +45,6 @@ DEAL_II_NAMESPACE_OPEN
 
 namespace internal
 {
-  template <class VectorType>
-  typename VectorType::value_type inline get_vector_element(
-    const VectorType &            vector,
-    const types::global_dof_index cell_number)
-  {
-    return internal::ElementAccess<VectorType>::get(vector, cell_number);
-  }
-
-
-
-  IndexSet::value_type inline get_vector_element(
-    const IndexSet &              is,
-    const types::global_dof_index cell_number)
-  {
-    return (is.is_element(cell_number) ? 1 : 0);
-  }
-
-
-
   template <int dim, int spacedim>
   inline std::vector<unsigned int>
   make_shape_function_to_row_table(const FiniteElement<dim, spacedim> &fe)
@@ -1643,12 +1613,11 @@ namespace FEValuesViews
 
 
   template <int dim, int spacedim>
-  template <class InputVector>
+  template <typename Number>
   void
   Scalar<dim, spacedim>::get_function_hessians(
-    const InputVector &fe_function,
-    std::vector<solution_hessian_type<typename InputVector::value_type>>
-      &hessians) const
+    const ReadVector<Number> &                  fe_function,
+    std::vector<solution_hessian_type<Number>> &hessians) const
   {
     Assert(fe_values->update_flags & update_hessians,
            (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField(
@@ -1659,8 +1628,7 @@ namespace FEValuesViews
                     fe_values->present_cell.n_dofs_for_dof_handler());
 
     // get function values of dofs on this cell
-    dealii::Vector<typename InputVector::value_type> dof_values(
-      fe_values->dofs_per_cell);
+    dealii::Vector<Number> dof_values(fe_values->dofs_per_cell);
     fe_values->present_cell.get_interpolated_dof_values(fe_function,
                                                         dof_values);
     internal::do_function_derivatives<2, dim, spacedim>(
@@ -1697,12 +1665,11 @@ namespace FEValuesViews
 
 
   template <int dim, int spacedim>
-  template <class InputVector>
+  template <typename Number>
   void
   Scalar<dim, spacedim>::get_function_laplacians(
-    const InputVector &fe_function,
-    std::vector<solution_laplacian_type<typename InputVector::value_type>>
-      &laplacians) const
+    const ReadVector<Number> &                    fe_function,
+    std::vector<solution_laplacian_type<Number>> &laplacians) const
   {
     Assert(fe_values->update_flags & update_hessians,
            (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField(
@@ -1713,8 +1680,7 @@ namespace FEValuesViews
                     fe_values->present_cell.n_dofs_for_dof_handler());
 
     // get function values of dofs on this cell
-    dealii::Vector<typename InputVector::value_type> dof_values(
-      fe_values->dofs_per_cell);
+    dealii::Vector<Number> dof_values(fe_values->dofs_per_cell);
     fe_values->present_cell.get_interpolated_dof_values(fe_function,
                                                         dof_values);
     internal::do_function_laplacians<dim, spacedim>(
@@ -1751,13 +1717,12 @@ namespace FEValuesViews
 
 
   template <int dim, int spacedim>
-  template <class InputVector>
+  template <typename Number>
   void
   Scalar<dim, spacedim>::get_function_third_derivatives(
-    const InputVector &fe_function,
-    std::vector<
-      solution_third_derivative_type<typename InputVector::value_type>>
-      &third_derivatives) const
+    const ReadVector<Number> &                           fe_function,
+    std::vector<solution_third_derivative_type<Number>> &third_derivatives)
+    const
   {
     Assert(fe_values->update_flags & update_3rd_derivatives,
            (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField(
@@ -1768,8 +1733,7 @@ namespace FEValuesViews
                     fe_values->present_cell.n_dofs_for_dof_handler());
 
     // get function values of dofs on this cell
-    dealii::Vector<typename InputVector::value_type> dof_values(
-      fe_values->dofs_per_cell);
+    dealii::Vector<Number> dof_values(fe_values->dofs_per_cell);
     fe_values->present_cell.get_interpolated_dof_values(fe_function,
                                                         dof_values);
     internal::do_function_derivatives<3, dim, spacedim>(
@@ -1909,13 +1873,12 @@ namespace FEValuesViews
 
 
   template <int dim, int spacedim>
-  template <class InputVector>
+  template <typename Number>
   void
   Vector<dim, spacedim>::get_function_symmetric_gradients(
-    const InputVector &fe_function,
-    std::vector<
-      solution_symmetric_gradient_type<typename InputVector::value_type>>
-      &symmetric_gradients) const
+    const ReadVector<Number> &                             fe_function,
+    std::vector<solution_symmetric_gradient_type<Number>> &symmetric_gradients)
+    const
   {
     Assert(fe_values->update_flags & update_gradients,
            (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField(
@@ -1926,8 +1889,7 @@ namespace FEValuesViews
                     fe_values->present_cell.n_dofs_for_dof_handler());
 
     // get function values of dofs on this cell
-    dealii::Vector<typename InputVector::value_type> dof_values(
-      fe_values->dofs_per_cell);
+    dealii::Vector<Number> dof_values(fe_values->dofs_per_cell);
     fe_values->present_cell.get_interpolated_dof_values(fe_function,
                                                         dof_values);
     internal::do_function_symmetric_gradients<dim, spacedim>(
@@ -1965,12 +1927,11 @@ namespace FEValuesViews
 
 
   template <int dim, int spacedim>
-  template <class InputVector>
+  template <typename Number>
   void
   Vector<dim, spacedim>::get_function_divergences(
-    const InputVector &fe_function,
-    std::vector<solution_divergence_type<typename InputVector::value_type>>
-      &divergences) const
+    const ReadVector<Number> &                     fe_function,
+    std::vector<solution_divergence_type<Number>> &divergences) const
   {
     Assert(fe_values->update_flags & update_gradients,
            (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField(
@@ -1982,8 +1943,7 @@ namespace FEValuesViews
 
     // get function values of dofs
     // on this cell
-    dealii::Vector<typename InputVector::value_type> dof_values(
-      fe_values->dofs_per_cell);
+    dealii::Vector<Number> dof_values(fe_values->dofs_per_cell);
     fe_values->present_cell.get_interpolated_dof_values(fe_function,
                                                         dof_values);
     internal::do_function_divergences<dim, spacedim>(
@@ -2020,12 +1980,11 @@ namespace FEValuesViews
 
 
   template <int dim, int spacedim>
-  template <class InputVector>
+  template <typename Number>
   void
   Vector<dim, spacedim>::get_function_curls(
-    const InputVector &fe_function,
-    std::vector<solution_curl_type<typename InputVector::value_type>> &curls)
-    const
+    const ReadVector<Number> &               fe_function,
+    std::vector<solution_curl_type<Number>> &curls) const
   {
     Assert(fe_values->update_flags & update_gradients,
            (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField(
@@ -2036,8 +1995,7 @@ namespace FEValuesViews
                     fe_values->present_cell.n_dofs_for_dof_handler());
 
     // get function values of dofs on this cell
-    dealii::Vector<typename InputVector::value_type> dof_values(
-      fe_values->dofs_per_cell);
+    dealii::Vector<Number> dof_values(fe_values->dofs_per_cell);
     fe_values->present_cell.get_interpolated_dof_values(fe_function,
                                                         dof_values);
     internal::do_function_curls<dim, spacedim>(
@@ -2074,12 +2032,11 @@ namespace FEValuesViews
 
 
   template <int dim, int spacedim>
-  template <class InputVector>
+  template <typename Number>
   void
   Vector<dim, spacedim>::get_function_hessians(
-    const InputVector &fe_function,
-    std::vector<solution_hessian_type<typename InputVector::value_type>>
-      &hessians) const
+    const ReadVector<Number> &                  fe_function,
+    std::vector<solution_hessian_type<Number>> &hessians) const
   {
     Assert(fe_values->update_flags & update_hessians,
            (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField(
@@ -2090,8 +2047,7 @@ namespace FEValuesViews
                     fe_values->present_cell.n_dofs_for_dof_handler());
 
     // get function values of dofs on this cell
-    dealii::Vector<typename InputVector::value_type> dof_values(
-      fe_values->dofs_per_cell);
+    dealii::Vector<Number> dof_values(fe_values->dofs_per_cell);
     fe_values->present_cell.get_interpolated_dof_values(fe_function,
                                                         dof_values);
     internal::do_function_derivatives<2, dim, spacedim>(
@@ -2128,12 +2084,11 @@ namespace FEValuesViews
 
 
   template <int dim, int spacedim>
-  template <class InputVector>
+  template <typename Number>
   void
   Vector<dim, spacedim>::get_function_laplacians(
-    const InputVector &fe_function,
-    std::vector<solution_value_type<typename InputVector::value_type>>
-      &laplacians) const
+    const ReadVector<Number> &                fe_function,
+    std::vector<solution_value_type<Number>> &laplacians) const
   {
     Assert(fe_values->update_flags & update_hessians,
            (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField(
@@ -2149,8 +2104,7 @@ namespace FEValuesViews
                            fe_values->present_cell.n_dofs_for_dof_handler()));
 
     // get function values of dofs on this cell
-    dealii::Vector<typename InputVector::value_type> dof_values(
-      fe_values->dofs_per_cell);
+    dealii::Vector<Number> dof_values(fe_values->dofs_per_cell);
     fe_values->present_cell.get_interpolated_dof_values(fe_function,
                                                         dof_values);
     internal::do_function_laplacians<dim, spacedim>(
@@ -2190,13 +2144,12 @@ namespace FEValuesViews
 
 
   template <int dim, int spacedim>
-  template <class InputVector>
+  template <typename Number>
   void
   Vector<dim, spacedim>::get_function_third_derivatives(
-    const InputVector &fe_function,
-    std::vector<
-      solution_third_derivative_type<typename InputVector::value_type>>
-      &third_derivatives) const
+    const ReadVector<Number> &                           fe_function,
+    std::vector<solution_third_derivative_type<Number>> &third_derivatives)
+    const
   {
     Assert(fe_values->update_flags & update_3rd_derivatives,
            (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField(
@@ -2207,8 +2160,7 @@ namespace FEValuesViews
                     fe_values->present_cell.n_dofs_for_dof_handler());
 
     // get function values of dofs on this cell
-    dealii::Vector<typename InputVector::value_type> dof_values(
-      fe_values->dofs_per_cell);
+    dealii::Vector<Number> dof_values(fe_values->dofs_per_cell);
     fe_values->present_cell.get_interpolated_dof_values(fe_function,
                                                         dof_values);
     internal::do_function_derivatives<3, dim, spacedim>(
@@ -2296,12 +2248,11 @@ namespace FEValuesViews
 
 
   template <int dim, int spacedim>
-  template <class InputVector>
+  template <typename Number>
   void
   SymmetricTensor<2, dim, spacedim>::get_function_divergences(
-    const InputVector &fe_function,
-    std::vector<solution_divergence_type<typename InputVector::value_type>>
-      &divergences) const
+    const ReadVector<Number> &                     fe_function,
+    std::vector<solution_divergence_type<Number>> &divergences) const
   {
     Assert(fe_values->update_flags & update_gradients,
            (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField(
@@ -2313,8 +2264,7 @@ namespace FEValuesViews
 
     // get function values of dofs
     // on this cell
-    dealii::Vector<typename InputVector::value_type> dof_values(
-      fe_values->dofs_per_cell);
+    dealii::Vector<Number> dof_values(fe_values->dofs_per_cell);
     fe_values->present_cell.get_interpolated_dof_values(fe_function,
                                                         dof_values);
     internal::do_function_divergences<dim, spacedim>(
@@ -2402,12 +2352,11 @@ namespace FEValuesViews
 
 
   template <int dim, int spacedim>
-  template <class InputVector>
+  template <typename Number>
   void
   Tensor<2, dim, spacedim>::get_function_divergences(
-    const InputVector &fe_function,
-    std::vector<solution_divergence_type<typename InputVector::value_type>>
-      &divergences) const
+    const ReadVector<Number> &                     fe_function,
+    std::vector<solution_divergence_type<Number>> &divergences) const
   {
     Assert(fe_values->update_flags & update_gradients,
            (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField(
@@ -2419,8 +2368,7 @@ namespace FEValuesViews
 
     // get function values of dofs
     // on this cell
-    dealii::Vector<typename InputVector::value_type> dof_values(
-      fe_values->dofs_per_cell);
+    dealii::Vector<Number> dof_values(fe_values->dofs_per_cell);
     fe_values->present_cell.get_interpolated_dof_values(fe_function,
                                                         dof_values);
     internal::do_function_divergences<dim, spacedim>(
@@ -3469,7 +3417,7 @@ FEValuesBase<dim, spacedim>::get_function_gradients(
   auto view = make_array_view(dof_values.begin(), dof_values.end());
   fe_function.extract_subvector_to(indices, view);
   internal::do_function_derivatives(
-    make_array_view(dof_values.begin(), dof_values.end()),
+    view,
     this->finite_element_output.shape_gradients,
     *fe,
     this->finite_element_output.shape_function_to_row_table,
@@ -3481,14 +3429,12 @@ FEValuesBase<dim, spacedim>::get_function_gradients(
 
 
 template <int dim, int spacedim>
-template <class InputVector>
+template <typename Number>
 void
 FEValuesBase<dim, spacedim>::get_function_hessians(
-  const InputVector &fe_function,
-  std::vector<Tensor<2, spacedim, typename InputVector::value_type>> &hessians)
-  const
+  const ReadVector<Number> &                fe_function,
+  std::vector<Tensor<2, spacedim, Number>> &hessians) const
 {
-  using Number = typename InputVector::value_type;
   AssertDimension(fe->n_components(), 1);
   Assert(this->update_flags & update_hessians,
          ExcAccessToUninitializedField("update_hessians"));
@@ -3507,25 +3453,22 @@ FEValuesBase<dim, spacedim>::get_function_hessians(
 
 
 template <int dim, int spacedim>
-template <class InputVector>
+template <typename Number>
 void
 FEValuesBase<dim, spacedim>::get_function_hessians(
-  const InputVector &                             fe_function,
+  const ReadVector<Number> &                      fe_function,
   const ArrayView<const types::global_dof_index> &indices,
-  std::vector<Tensor<2, spacedim, typename InputVector::value_type>> &hessians)
-  const
+  std::vector<Tensor<2, spacedim, Number>> &      hessians) const
 {
-  using Number = typename InputVector::value_type;
   Assert(this->update_flags & update_hessians,
          ExcAccessToUninitializedField("update_hessians"));
   AssertDimension(fe_function.size(), present_cell.n_dofs_for_dof_handler());
   AssertDimension(indices.size(), dofs_per_cell);
 
   boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
-  for (unsigned int i = 0; i < dofs_per_cell; ++i)
-    dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
-  internal::do_function_derivatives(make_array_view(dof_values.begin(),
-                                                    dof_values.end()),
+  auto view = make_array_view(dof_values.begin(), dof_values.end());
+  fe_function.extract_subvector_to(indices, view);
+  internal::do_function_derivatives(view,
                                     this->finite_element_output.shape_hessians,
                                     hessians);
 }
@@ -3533,16 +3476,13 @@ FEValuesBase<dim, spacedim>::get_function_hessians(
 
 
 template <int dim, int spacedim>
-template <class InputVector>
+template <typename Number>
 void
 FEValuesBase<dim, spacedim>::get_function_hessians(
-  const InputVector &fe_function,
-  std::vector<
-    std::vector<Tensor<2, spacedim, typename InputVector::value_type>>>
-    &        hessians,
+  const ReadVector<Number> &                             fe_function,
+  std::vector<std::vector<Tensor<2, spacedim, Number>>> &hessians,
   const bool quadrature_points_fastest) const
 {
-  using Number = typename InputVector::value_type;
   Assert(this->update_flags & update_hessians,
          ExcAccessToUninitializedField("update_hessians"));
   Assert(present_cell.is_initialized(), ExcNotReinited());
@@ -3563,26 +3503,24 @@ FEValuesBase<dim, spacedim>::get_function_hessians(
 
 
 template <int dim, int spacedim>
-template <class InputVector>
+template <typename Number>
 void
 FEValuesBase<dim, spacedim>::get_function_hessians(
-  const InputVector &                             fe_function,
-  const ArrayView<const types::global_dof_index> &indices,
-  ArrayView<std::vector<Tensor<2, spacedim, typename InputVector::value_type>>>
-             hessians,
+  const ReadVector<Number> &                          fe_function,
+  const ArrayView<const types::global_dof_index> &    indices,
+  ArrayView<std::vector<Tensor<2, spacedim, Number>>> hessians,
   const bool quadrature_points_fastest) const
 {
-  using Number = typename InputVector::value_type;
   Assert(this->update_flags & update_hessians,
          ExcAccessToUninitializedField("update_hessians"));
   Assert(indices.size() % dofs_per_cell == 0,
          ExcNotMultiple(indices.size(), dofs_per_cell));
 
   boost::container::small_vector<Number, 200> dof_values(indices.size());
-  for (unsigned int i = 0; i < indices.size(); ++i)
-    dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
+  auto view = make_array_view(dof_values.begin(), dof_values.end());
+  fe_function.extract_subvector_to(indices, view);
   internal::do_function_derivatives(
-    make_array_view(dof_values.begin(), dof_values.end()),
+    view,
     this->finite_element_output.shape_hessians,
     *fe,
     this->finite_element_output.shape_function_to_row_table,
@@ -3594,13 +3532,12 @@ FEValuesBase<dim, spacedim>::get_function_hessians(
 
 
 template <int dim, int spacedim>
-template <class InputVector>
+template <typename Number>
 void
 FEValuesBase<dim, spacedim>::get_function_laplacians(
-  const InputVector &                            fe_function,
-  std::vector<typename InputVector::value_type> &laplacians) const
+  const ReadVector<Number> &fe_function,
+  std::vector<Number> &     laplacians) const
 {
-  using Number = typename InputVector::value_type;
   Assert(this->update_flags & update_hessians,
          ExcAccessToUninitializedField("update_hessians"));
   AssertDimension(fe->n_components(), 1);
@@ -3619,24 +3556,22 @@ FEValuesBase<dim, spacedim>::get_function_laplacians(
 
 
 template <int dim, int spacedim>
-template <class InputVector>
+template <typename Number>
 void
 FEValuesBase<dim, spacedim>::get_function_laplacians(
-  const InputVector &                             fe_function,
+  const ReadVector<Number> &                      fe_function,
   const ArrayView<const types::global_dof_index> &indices,
-  std::vector<typename InputVector::value_type> & laplacians) const
+  std::vector<Number> &                           laplacians) const
 {
-  using Number = typename InputVector::value_type;
   Assert(this->update_flags & update_hessians,
          ExcAccessToUninitializedField("update_hessians"));
   AssertDimension(fe->n_components(), 1);
   AssertDimension(indices.size(), dofs_per_cell);
 
   boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
-  for (unsigned int i = 0; i < dofs_per_cell; ++i)
-    dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
-  internal::do_function_laplacians(make_array_view(dof_values.begin(),
-                                                   dof_values.end()),
+  auto view = make_array_view(dof_values.begin(), dof_values.end());
+  fe_function.extract_subvector_to(indices, view);
+  internal::do_function_laplacians(view,
                                    this->finite_element_output.shape_hessians,
                                    laplacians);
 }
@@ -3644,13 +3579,12 @@ FEValuesBase<dim, spacedim>::get_function_laplacians(
 
 
 template <int dim, int spacedim>
-template <class InputVector>
+template <typename Number>
 void
 FEValuesBase<dim, spacedim>::get_function_laplacians(
-  const InputVector &                                    fe_function,
-  std::vector<Vector<typename InputVector::value_type>> &laplacians) const
+  const ReadVector<Number> &   fe_function,
+  std::vector<Vector<Number>> &laplacians) const
 {
-  using Number = typename InputVector::value_type;
   Assert(present_cell.is_initialized(), ExcNotReinited());
   Assert(this->update_flags & update_hessians,
          ExcAccessToUninitializedField("update_hessians"));
@@ -3670,14 +3604,13 @@ FEValuesBase<dim, spacedim>::get_function_laplacians(
 
 
 template <int dim, int spacedim>
-template <class InputVector>
+template <typename Number>
 void
 FEValuesBase<dim, spacedim>::get_function_laplacians(
-  const InputVector &                                    fe_function,
-  const ArrayView<const types::global_dof_index> &       indices,
-  std::vector<Vector<typename InputVector::value_type>> &laplacians) const
+  const ReadVector<Number> &                      fe_function,
+  const ArrayView<const types::global_dof_index> &indices,
+  std::vector<Vector<Number>> &                   laplacians) const
 {
-  using Number = typename InputVector::value_type;
   // 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,
@@ -3686,10 +3619,10 @@ FEValuesBase<dim, spacedim>::get_function_laplacians(
          ExcAccessToUninitializedField("update_hessians"));
 
   boost::container::small_vector<Number, 200> dof_values(indices.size());
-  for (unsigned int i = 0; i < indices.size(); ++i)
-    dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
+  auto view = make_array_view(dof_values.begin(), dof_values.end());
+  fe_function.extract_subvector_to(indices, view);
   internal::do_function_laplacians(
-    make_array_view(dof_values.begin(), dof_values.end()),
+    view,
     this->finite_element_output.shape_hessians,
     *fe,
     this->finite_element_output.shape_function_to_row_table,
@@ -3701,25 +3634,24 @@ FEValuesBase<dim, spacedim>::get_function_laplacians(
 
 
 template <int dim, int spacedim>
-template <class InputVector>
+template <typename Number>
 void
 FEValuesBase<dim, spacedim>::get_function_laplacians(
-  const InputVector &                                         fe_function,
-  const ArrayView<const types::global_dof_index> &            indices,
-  std::vector<std::vector<typename InputVector::value_type>> &laplacians,
+  const ReadVector<Number> &                      fe_function,
+  const ArrayView<const types::global_dof_index> &indices,
+  std::vector<std::vector<Number>> &              laplacians,
   const bool quadrature_points_fastest) const
 {
-  using Number = typename InputVector::value_type;
   Assert(indices.size() % dofs_per_cell == 0,
          ExcNotMultiple(indices.size(), dofs_per_cell));
   Assert(this->update_flags & update_hessians,
          ExcAccessToUninitializedField("update_hessians"));
 
   boost::container::small_vector<Number, 200> dof_values(indices.size());
-  for (unsigned int i = 0; i < indices.size(); ++i)
-    dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
+  auto view = make_array_view(dof_values.begin(), dof_values.end());
+  fe_function.extract_subvector_to(indices, view);
   internal::do_function_laplacians(
-    make_array_view(dof_values.begin(), dof_values.end()),
+    view,
     this->finite_element_output.shape_hessians,
     *fe,
     this->finite_element_output.shape_function_to_row_table,
@@ -3731,14 +3663,12 @@ FEValuesBase<dim, spacedim>::get_function_laplacians(
 
 
 template <int dim, int spacedim>
-template <class InputVector>
+template <typename Number>
 void
 FEValuesBase<dim, spacedim>::get_function_third_derivatives(
-  const InputVector &fe_function,
-  std::vector<Tensor<3, spacedim, typename InputVector::value_type>>
-    &third_derivatives) const
+  const ReadVector<Number> &                fe_function,
+  std::vector<Tensor<3, spacedim, Number>> &third_derivatives) const
 {
-  using Number = typename InputVector::value_type;
   AssertDimension(fe->n_components(), 1);
   Assert(this->update_flags & update_3rd_derivatives,
          ExcAccessToUninitializedField("update_3rd_derivatives"));
@@ -3757,42 +3687,35 @@ FEValuesBase<dim, spacedim>::get_function_third_derivatives(
 
 
 template <int dim, int spacedim>
-template <class InputVector>
+template <typename Number>
 void
 FEValuesBase<dim, spacedim>::get_function_third_derivatives(
-  const InputVector &                             fe_function,
+  const ReadVector<Number> &                      fe_function,
   const ArrayView<const types::global_dof_index> &indices,
-  std::vector<Tensor<3, spacedim, typename InputVector::value_type>>
-    &third_derivatives) const
+  std::vector<Tensor<3, spacedim, Number>> &      third_derivatives) const
 {
-  using Number = typename InputVector::value_type;
   Assert(this->update_flags & update_3rd_derivatives,
          ExcAccessToUninitializedField("update_3rd_derivatives"));
   AssertDimension(fe_function.size(), present_cell.n_dofs_for_dof_handler());
   AssertDimension(indices.size(), dofs_per_cell);
 
   boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
-  for (unsigned int i = 0; i < dofs_per_cell; ++i)
-    dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
+  auto view = make_array_view(dof_values.begin(), dof_values.end());
+  fe_function.extract_subvector_to(indices, view);
   internal::do_function_derivatives(
-    make_array_view(dof_values.begin(), dof_values.end()),
-    this->finite_element_output.shape_3rd_derivatives,
-    third_derivatives);
+    view, this->finite_element_output.shape_3rd_derivatives, third_derivatives);
 }
 
 
 
 template <int dim, int spacedim>
-template <class InputVector>
+template <typename Number>
 void
 FEValuesBase<dim, spacedim>::get_function_third_derivatives(
-  const InputVector &fe_function,
-  std::vector<
-    std::vector<Tensor<3, spacedim, typename InputVector::value_type>>>
-    &        third_derivatives,
+  const ReadVector<Number> &                             fe_function,
+  std::vector<std::vector<Tensor<3, spacedim, Number>>> &third_derivatives,
   const bool quadrature_points_fastest) const
 {
-  using Number = typename InputVector::value_type;
   Assert(this->update_flags & update_3rd_derivatives,
          ExcAccessToUninitializedField("update_3rd_derivatives"));
   Assert(present_cell.is_initialized(), ExcNotReinited());
@@ -3813,26 +3736,24 @@ FEValuesBase<dim, spacedim>::get_function_third_derivatives(
 
 
 template <int dim, int spacedim>
-template <class InputVector>
+template <typename Number>
 void
 FEValuesBase<dim, spacedim>::get_function_third_derivatives(
-  const InputVector &                             fe_function,
-  const ArrayView<const types::global_dof_index> &indices,
-  ArrayView<std::vector<Tensor<3, spacedim, typename InputVector::value_type>>>
-             third_derivatives,
+  const ReadVector<Number> &                          fe_function,
+  const ArrayView<const types::global_dof_index> &    indices,
+  ArrayView<std::vector<Tensor<3, spacedim, Number>>> third_derivatives,
   const bool quadrature_points_fastest) const
 {
-  using Number = typename InputVector::value_type;
   Assert(this->update_flags & update_3rd_derivatives,
          ExcAccessToUninitializedField("update_3rd_derivatives"));
   Assert(indices.size() % dofs_per_cell == 0,
          ExcNotMultiple(indices.size(), dofs_per_cell));
 
   boost::container::small_vector<Number, 200> dof_values(indices.size());
-  for (unsigned int i = 0; i < indices.size(); ++i)
-    dof_values[i] = internal::get_vector_element(fe_function, indices[i]);
+  auto view = make_array_view(dof_values.begin(), dof_values.end());
+  fe_function.extract_subvector_to(indices, view);
   internal::do_function_derivatives(
-    make_array_view(dof_values.begin(), dof_values.end()),
+    view,
     this->finite_element_output.shape_3rd_derivatives,
     *fe,
     this->finite_element_output.shape_function_to_row_table,
diff --git a/source/fe/fe_values.inst.in b/source/fe/fe_values.inst.in
index 6fb81be92f..d75f616bd7 100644
--- a/source/fe/fe_values.inst.in
+++ b/source/fe/fe_values.inst.in
@@ -198,93 +198,82 @@ for (S : REAL_AND_COMPLEX_SCALARS; deal_II_dimension : DIMENSIONS;
 
 
 
-for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS;
+for (S : REAL_AND_COMPLEX_SCALARS; deal_II_dimension : DIMENSIONS;
      deal_II_space_dimension : SPACE_DIMENSIONS)
   {
 #  if deal_II_dimension <= deal_II_space_dimension
     template void
     FEValuesViews::Scalar<deal_II_dimension, deal_II_space_dimension>::
-      get_function_hessians<dealii::VEC>(
-        const dealii::VEC &,
+      get_function_hessians<S>(
+        const dealii::ReadVector<S> &,
         std::vector<
-          ProductType<dealii::VEC::value_type,
-                      dealii::Tensor<2, deal_II_space_dimension>>::type> &)
+          ProductType<S, dealii::Tensor<2, deal_II_space_dimension>>::type> &)
         const;
     template void
     FEValuesViews::Scalar<deal_II_dimension, deal_II_space_dimension>::
-      get_function_laplacians<dealii::VEC>(
-        const dealii::VEC &,
-        std::vector<ProductType<dealii::VEC::value_type, double>::type> &)
+      get_function_laplacians<S>(const dealii::ReadVector<S> &,
+                                 std::vector<ProductType<S, double>::type> &)
         const;
     template void
     FEValuesViews::Scalar<deal_II_dimension, deal_II_space_dimension>::
-      get_function_third_derivatives<dealii::VEC>(
-        const dealii::VEC &,
+      get_function_third_derivatives<S>(
+        const dealii::ReadVector<S> &,
         std::vector<
-          ProductType<dealii::VEC::value_type,
-                      dealii::Tensor<3, deal_II_space_dimension>>::type> &)
+          ProductType<S, dealii::Tensor<3, deal_II_space_dimension>>::type> &)
         const;
 
     template void
     FEValuesViews::Vector<deal_II_dimension, deal_II_space_dimension>::
-      get_function_symmetric_gradients<dealii::VEC>(
-        const dealii::VEC &,
+      get_function_symmetric_gradients<S>(
+        const dealii::ReadVector<S> &,
         std::vector<ProductType<
-          dealii::VEC::value_type,
+          S,
           dealii::SymmetricTensor<2, deal_II_space_dimension>>::type> &) const;
+    template void FEValuesViews::
+      Vector<deal_II_dimension, deal_II_space_dimension>::get_function_curls<S>(
+        const dealii::ReadVector<S> &,
+        std::vector<ProductType<S, curl_type>::type> &) const;
     template void
     FEValuesViews::Vector<deal_II_dimension, deal_II_space_dimension>::
-      get_function_curls<dealii::VEC>(
-        const dealii::VEC &,
-        std::vector<ProductType<dealii::VEC::value_type, curl_type>::type> &)
-        const;
-    template void
-    FEValuesViews::Vector<deal_II_dimension, deal_II_space_dimension>::
-      get_function_divergences<dealii::VEC>(
-        const dealii::VEC &,
-        std::vector<ProductType<dealii::VEC::value_type, divergence_type>::type>
-          &) const;
+      get_function_divergences<S>(
+        const dealii::ReadVector<S> &,
+        std::vector<ProductType<S, divergence_type>::type> &) const;
     template void
     FEValuesViews::Vector<deal_II_dimension, deal_II_space_dimension>::
-      get_function_hessians<dealii::VEC>(
-        const dealii::VEC &,
+      get_function_hessians<S>(
+        const dealii::ReadVector<S> &,
         std::vector<
-          ProductType<dealii::VEC::value_type,
-                      dealii::Tensor<3, deal_II_space_dimension>>::type> &)
+          ProductType<S, dealii::Tensor<3, deal_II_space_dimension>>::type> &)
         const;
     template void
     FEValuesViews::Vector<deal_II_dimension, deal_II_space_dimension>::
-      get_function_laplacians<dealii::VEC>(
-        const dealii::VEC &,
+      get_function_laplacians<S>(
+        const dealii::ReadVector<S> &,
         std::vector<
-          ProductType<dealii::VEC::value_type,
-                      dealii::Tensor<1, deal_II_space_dimension>>::type> &)
+          ProductType<S, dealii::Tensor<1, deal_II_space_dimension>>::type> &)
         const;
     template void
     FEValuesViews::Vector<deal_II_dimension, deal_II_space_dimension>::
-      get_function_third_derivatives<dealii::VEC>(
-        const dealii::VEC &,
+      get_function_third_derivatives<S>(
+        const dealii::ReadVector<S> &,
         std::vector<
-          ProductType<dealii::VEC::value_type,
-                      dealii::Tensor<4, deal_II_space_dimension>>::type> &)
+          ProductType<S, dealii::Tensor<4, deal_II_space_dimension>>::type> &)
         const;
 
     template void FEValuesViews::
       SymmetricTensor<2, deal_II_dimension, deal_II_space_dimension>::
-        get_function_divergences<dealii::VEC>(
-          const dealii::VEC &,
+        get_function_divergences<S>(
+          const dealii::ReadVector<S> &,
           std::vector<
-            ProductType<dealii::VEC::value_type,
-                        dealii::Tensor<1, deal_II_space_dimension>>::type> &)
+            ProductType<S, dealii::Tensor<1, deal_II_space_dimension>>::type> &)
           const;
 
     template void
     FEValuesViews::Tensor<2, deal_II_dimension, deal_II_space_dimension>::
-      get_function_divergences<dealii::VEC>(
-        const dealii::VEC &,
+      get_function_divergences<S>(
+        const dealii::ReadVector<S> &,
         std::vector<
-          ProductType<dealii::VEC::value_type,
-                      dealii::Tensor<1, deal_II_space_dimension>>::type> &)
+          ProductType<S, dealii::Tensor<1, deal_II_space_dimension>>::type> &)
         const;
 #  endif
   }
@@ -406,88 +395,80 @@ for (VEC : GENERAL_CONTAINER_TYPES; Number : ALL_SCALAR_TYPES;
   }
 
 
-for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS;
+for (S : REAL_AND_COMPLEX_SCALARS; deal_II_dimension : DIMENSIONS;
      deal_II_space_dimension : SPACE_DIMENSIONS)
   {
 #  if deal_II_dimension <= deal_II_space_dimension
 
     template void FEValuesBase<deal_II_dimension, deal_II_space_dimension>::
-      get_function_hessians<VEC>(
-        const VEC &,
-        std::vector<dealii::Tensor<2, deal_II_space_dimension, VEC::value_type>>
-          &) const;
+      get_function_hessians<S>(
+        const ReadVector<S> &,
+        std::vector<dealii::Tensor<2, deal_II_space_dimension, S>> &) const;
     template void FEValuesBase<deal_II_dimension, deal_II_space_dimension>::
-      get_function_hessians<VEC>(
-        const VEC &,
+      get_function_hessians<S>(
+        const ReadVector<S> &,
         const ArrayView<const types::global_dof_index> &,
-        std::vector<dealii::Tensor<2, deal_II_space_dimension, VEC::value_type>>
-          &) const;
+        std::vector<dealii::Tensor<2, deal_II_space_dimension, S>> &) const;
 
     template void FEValuesBase<deal_II_dimension, deal_II_space_dimension>::
-      get_function_hessians<VEC>(
-        const VEC &,
-        std::vector<std::vector<
-          dealii::Tensor<2, deal_II_space_dimension, VEC::value_type>>> &,
+      get_function_hessians<S>(
+        const ReadVector<S> &,
+        std::vector<std::vector<dealii::Tensor<2, deal_II_space_dimension, S>>>
+          &,
         bool) const;
     template void FEValuesBase<deal_II_dimension, deal_II_space_dimension>::
-      get_function_hessians<VEC>(
-        const VEC &,
+      get_function_hessians<S>(
+        const ReadVector<S> &,
         const ArrayView<const types::global_dof_index> &,
-        ArrayView<std::vector<
-          dealii::Tensor<2, deal_II_space_dimension, VEC::value_type>>>,
+        ArrayView<std::vector<dealii::Tensor<2, deal_II_space_dimension, S>>>,
         bool) const;
 
     template void FEValuesBase<deal_II_dimension, deal_II_space_dimension>::
-      get_function_laplacians<VEC>(const VEC &, std::vector<VEC::value_type> &)
-        const;
+      get_function_laplacians<S>(const ReadVector<S> &, std::vector<S> &) const;
     template void FEValuesBase<deal_II_dimension, deal_II_space_dimension>::
-      get_function_laplacians<VEC>(
-        const VEC &,
+      get_function_laplacians<S>(
+        const ReadVector<S> &,
         const ArrayView<const types::global_dof_index> &,
-        std::vector<VEC::value_type> &) const;
+        std::vector<S> &) const;
 
     template void FEValuesBase<deal_II_dimension, deal_II_space_dimension>::
-      get_function_laplacians<VEC>(const VEC &,
-                                   std::vector<Vector<VEC::value_type>> &)
-        const;
+      get_function_laplacians<S>(const ReadVector<S> &,
+                                 std::vector<Vector<S>> &) const;
 
     template void FEValuesBase<deal_II_dimension, deal_II_space_dimension>::
-      get_function_laplacians<VEC>(
-        const VEC &,
+      get_function_laplacians<S>(
+        const ReadVector<S> &,
         const ArrayView<const types::global_dof_index> &,
-        std::vector<Vector<VEC::value_type>> &) const;
+        std::vector<Vector<S>> &) const;
 
     template void FEValuesBase<deal_II_dimension, deal_II_space_dimension>::
-      get_function_laplacians<VEC>(
-        const VEC &,
+      get_function_laplacians<S>(
+        const ReadVector<S> &,
         const ArrayView<const types::global_dof_index> &,
-        std::vector<std::vector<VEC::value_type>> &,
+        std::vector<std::vector<S>> &,
         bool) const;
 
     template void FEValuesBase<deal_II_dimension, deal_II_space_dimension>::
-      get_function_third_derivatives<VEC>(
-        const VEC &,
-        std::vector<dealii::Tensor<3, deal_II_space_dimension, VEC::value_type>>
-          &) const;
+      get_function_third_derivatives<S>(
+        const ReadVector<S> &,
+        std::vector<dealii::Tensor<3, deal_II_space_dimension, S>> &) const;
     template void FEValuesBase<deal_II_dimension, deal_II_space_dimension>::
-      get_function_third_derivatives<VEC>(
-        const VEC &,
+      get_function_third_derivatives<S>(
+        const ReadVector<S> &,
         const ArrayView<const types::global_dof_index> &,
-        std::vector<dealii::Tensor<3, deal_II_space_dimension, VEC::value_type>>
-          &) const;
+        std::vector<dealii::Tensor<3, deal_II_space_dimension, S>> &) const;
 
     template void FEValuesBase<deal_II_dimension, deal_II_space_dimension>::
-      get_function_third_derivatives<VEC>(
-        const VEC &,
-        std::vector<std::vector<
-          dealii::Tensor<3, deal_II_space_dimension, VEC::value_type>>> &,
+      get_function_third_derivatives<S>(
+        const ReadVector<S> &,
+        std::vector<std::vector<dealii::Tensor<3, deal_II_space_dimension, S>>>
+          &,
         bool) const;
     template void FEValuesBase<deal_II_dimension, deal_II_space_dimension>::
-      get_function_third_derivatives<VEC>(
-        const VEC &,
+      get_function_third_derivatives<S>(
+        const ReadVector<S> &,
         const ArrayView<const types::global_dof_index> &,
-        ArrayView<std::vector<
-          dealii::Tensor<3, deal_II_space_dimension, VEC::value_type>>>,
+        ArrayView<std::vector<dealii::Tensor<3, deal_II_space_dimension, S>>>,
         bool) const;
 
 #  endif
-- 
2.39.5