From: David Wells <drwells@email.unc.edu>
Date: Sat, 13 May 2023 13:45:02 +0000 (-0400)
Subject: Make FEValues::get_function_gradients() use ReadVector.
X-Git-Tag: relicensing~820^2~6
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=7fbe3db5c56d0ddb358123234483260d49b3186f;p=dealii.git

Make FEValues::get_function_gradients() use ReadVector.
---

diff --git a/include/deal.II/fe/fe_values.h b/include/deal.II/fe/fe_values.h
index 9dd445848c..218b391153 100644
--- a/include/deal.II/fe/fe_values.h
+++ b/include/deal.II/fe/fe_values.h
@@ -484,12 +484,11 @@ namespace FEValuesViews
      *
      * @dealiiRequiresUpdateFlags{update_gradients}
      */
-    template <class InputVector>
+    template <typename Number>
     void
     get_function_gradients(
-      const InputVector &fe_function,
-      std::vector<solution_gradient_type<typename InputVector::value_type>>
-        &gradients) const;
+      const ReadVector<Number> &                   fe_function,
+      std::vector<solution_gradient_type<Number>> &gradients) const;
 
     /**
      * This function relates to get_function_gradients() in the same way
@@ -1171,12 +1170,11 @@ namespace FEValuesViews
      *
      * @dealiiRequiresUpdateFlags{update_gradients}
      */
-    template <class InputVector>
+    template <typename Number>
     void
     get_function_gradients(
-      const InputVector &fe_function,
-      std::vector<solution_gradient_type<typename InputVector::value_type>>
-        &gradients) const;
+      const ReadVector<Number> &                   fe_function,
+      std::vector<solution_gradient_type<Number>> &gradients) const;
 
     /**
      * This function relates to get_function_gradients() in the same way
@@ -2157,12 +2155,11 @@ namespace FEValuesViews
      *
      * @dealiiRequiresUpdateFlags{update_gradients}
      */
-    template <class InputVector>
+    template <typename Number>
     void
     get_function_gradients(
-      const InputVector &fe_function,
-      std::vector<solution_gradient_type<typename InputVector::value_type>>
-        &gradients) const;
+      const ReadVector<Number> &                   fe_function,
+      std::vector<solution_gradient_type<Number>> &gradients) const;
 
     /**
      * This function relates to get_function_gradients() in the same way
@@ -2887,12 +2884,11 @@ public:
    *
    * @dealiiRequiresUpdateFlags{update_gradients}
    */
-  template <class InputVector>
+  template <typename Number>
   void
   get_function_gradients(
-    const InputVector &fe_function,
-    std::vector<Tensor<1, spacedim, typename InputVector::value_type>>
-      &gradients) const;
+    const ReadVector<Number> &                fe_function,
+    std::vector<Tensor<1, spacedim, Number>> &gradients) const;
 
   /**
    * This function does the same as the other get_function_gradients(), but
@@ -2910,13 +2906,11 @@ public:
    *
    * @dealiiRequiresUpdateFlags{update_gradients}
    */
-  template <class InputVector>
+  template <typename Number>
   void
   get_function_gradients(
-    const InputVector &fe_function,
-    std::vector<
-      std::vector<Tensor<1, spacedim, typename InputVector::value_type>>>
-      &gradients) const;
+    const ReadVector<Number> &                             fe_function,
+    std::vector<std::vector<Tensor<1, spacedim, Number>>> &gradients) const;
 
   /**
    * This function relates to the first of the get_function_gradients() function
@@ -2926,13 +2920,12 @@ public:
    *
    * @dealiiRequiresUpdateFlags{update_gradients}
    */
-  template <class InputVector>
+  template <typename Number>
   void
   get_function_gradients(
-    const InputVector &                             fe_function,
+    const ReadVector<Number> &                      fe_function,
     const ArrayView<const types::global_dof_index> &indices,
-    std::vector<Tensor<1, spacedim, typename InputVector::value_type>>
-      &gradients) const;
+    std::vector<Tensor<1, spacedim, Number>> &      gradients) const;
 
   /**
    * This function relates to the first of the get_function_gradients() function
@@ -2942,14 +2935,12 @@ public:
    *
    * @dealiiRequiresUpdateFlags{update_gradients}
    */
-  template <class InputVector>
+  template <typename Number>
   void
   get_function_gradients(
-    const InputVector &                             fe_function,
-    const ArrayView<const types::global_dof_index> &indices,
-    ArrayView<
-      std::vector<Tensor<1, spacedim, typename InputVector::value_type>>>
-               gradients,
+    const ReadVector<Number> &                          fe_function,
+    const ArrayView<const types::global_dof_index> &    indices,
+    ArrayView<std::vector<Tensor<1, spacedim, Number>>> gradients,
     const bool quadrature_points_fastest = false) const;
 
   /** @} */
diff --git a/source/fe/fe_values.cc b/source/fe/fe_values.cc
index 2f2117f344..7df91e0af0 100644
--- a/source/fe/fe_values.cc
+++ b/source/fe/fe_values.cc
@@ -1591,12 +1591,11 @@ namespace FEValuesViews
 
 
   template <int dim, int spacedim>
-  template <class InputVector>
+  template <typename Number>
   void
   Scalar<dim, spacedim>::get_function_gradients(
-    const InputVector &fe_function,
-    std::vector<solution_gradient_type<typename InputVector::value_type>>
-      &gradients) const
+    const ReadVector<Number> &                   fe_function,
+    std::vector<solution_gradient_type<Number>> &gradients) const
   {
     Assert(fe_values->update_flags & update_gradients,
            (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField(
@@ -1607,8 +1606,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<1, dim, spacedim>(
@@ -1621,7 +1619,7 @@ namespace FEValuesViews
 
 
   template <int dim, int spacedim>
-  template <class InputVector>
+  template <typename InputVector>
   void
   Scalar<dim, spacedim>::get_function_gradients_from_local_dof_values(
     const InputVector &dof_values,
@@ -1859,12 +1857,11 @@ namespace FEValuesViews
 
 
   template <int dim, int spacedim>
-  template <class InputVector>
+  template <typename Number>
   void
   Vector<dim, spacedim>::get_function_gradients(
-    const InputVector &fe_function,
-    std::vector<solution_gradient_type<typename InputVector::value_type>>
-      &gradients) const
+    const ReadVector<Number> &                   fe_function,
+    std::vector<solution_gradient_type<Number>> &gradients) const
   {
     Assert(fe_values->update_flags & update_gradients,
            (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField(
@@ -1875,8 +1872,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<1, dim, spacedim>(
@@ -1889,7 +1885,7 @@ namespace FEValuesViews
 
 
   template <int dim, int spacedim>
-  template <class InputVector>
+  template <typename InputVector>
   void
   Vector<dim, spacedim>::get_function_gradients_from_local_dof_values(
     const InputVector &dof_values,
@@ -2461,12 +2457,11 @@ namespace FEValuesViews
 
 
   template <int dim, int spacedim>
-  template <class InputVector>
+  template <typename Number>
   void
   Tensor<2, dim, spacedim>::get_function_gradients(
-    const InputVector &fe_function,
-    std::vector<solution_gradient_type<typename InputVector::value_type>>
-      &gradients) const
+    const ReadVector<Number> &                   fe_function,
+    std::vector<solution_gradient_type<Number>> &gradients) const
   {
     Assert(fe_values->update_flags & update_gradients,
            (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField(
@@ -2478,8 +2473,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_gradients<dim, spacedim>(
@@ -3384,14 +3378,12 @@ FEValuesBase<dim, spacedim>::get_function_values(
 
 
 template <int dim, int spacedim>
-template <class InputVector>
+template <typename Number>
 void
 FEValuesBase<dim, spacedim>::get_function_gradients(
-  const InputVector &fe_function,
-  std::vector<Tensor<1, spacedim, typename InputVector::value_type>> &gradients)
-  const
+  const ReadVector<Number> &                fe_function,
+  std::vector<Tensor<1, spacedim, Number>> &gradients) const
 {
-  using Number = typename InputVector::value_type;
   Assert(this->update_flags & update_gradients,
          ExcAccessToUninitializedField("update_gradients"));
   AssertDimension(fe->n_components(), 1);
@@ -3410,25 +3402,22 @@ FEValuesBase<dim, spacedim>::get_function_gradients(
 
 
 template <int dim, int spacedim>
-template <class InputVector>
+template <typename Number>
 void
 FEValuesBase<dim, spacedim>::get_function_gradients(
-  const InputVector &                             fe_function,
+  const ReadVector<Number> &                      fe_function,
   const ArrayView<const types::global_dof_index> &indices,
-  std::vector<Tensor<1, spacedim, typename InputVector::value_type>> &gradients)
-  const
+  std::vector<Tensor<1, spacedim, Number>> &      gradients) const
 {
-  using Number = typename InputVector::value_type;
   Assert(this->update_flags & update_gradients,
          ExcAccessToUninitializedField("update_gradients"));
   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_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_gradients,
                                     gradients);
 }
@@ -3436,15 +3425,12 @@ FEValuesBase<dim, spacedim>::get_function_gradients(
 
 
 template <int dim, int spacedim>
-template <class InputVector>
+template <typename Number>
 void
 FEValuesBase<dim, spacedim>::get_function_gradients(
-  const InputVector &fe_function,
-  std::vector<
-    std::vector<Tensor<1, spacedim, typename InputVector::value_type>>>
-    &gradients) const
+  const ReadVector<Number> &                             fe_function,
+  std::vector<std::vector<Tensor<1, spacedim, Number>>> &gradients) const
 {
-  using Number = typename InputVector::value_type;
   Assert(this->update_flags & update_gradients,
          ExcAccessToUninitializedField("update_gradients"));
   Assert(present_cell.is_initialized(), ExcNotReinited());
@@ -3464,16 +3450,14 @@ FEValuesBase<dim, spacedim>::get_function_gradients(
 
 
 template <int dim, int spacedim>
-template <class InputVector>
+template <typename Number>
 void
 FEValuesBase<dim, spacedim>::get_function_gradients(
-  const InputVector &                             fe_function,
-  const ArrayView<const types::global_dof_index> &indices,
-  ArrayView<std::vector<Tensor<1, spacedim, typename InputVector::value_type>>>
-             gradients,
+  const ReadVector<Number> &                          fe_function,
+  const ArrayView<const types::global_dof_index> &    indices,
+  ArrayView<std::vector<Tensor<1, spacedim, Number>>> gradients,
   const bool quadrature_points_fastest) 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,
@@ -3481,9 +3465,9 @@ FEValuesBase<dim, spacedim>::get_function_gradients(
   Assert(this->update_flags & update_gradients,
          ExcAccessToUninitializedField("update_gradients"));
 
-  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]);
+  boost::container::small_vector<Number, 200> dof_values(dofs_per_cell);
+  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_gradients,
diff --git a/source/fe/fe_values.inst.in b/source/fe/fe_values.inst.in
index 54ba212b8f..6fb81be92f 100644
--- a/source/fe/fe_values.inst.in
+++ b/source/fe/fe_values.inst.in
@@ -145,6 +145,54 @@ for (S : REAL_AND_COMPLEX_SCALARS; deal_II_dimension : DIMENSIONS;
                              const ArrayView<const types::global_dof_index> &,
                              ArrayView<std::vector<S>>,
                              bool) const;
+
+    template void
+    FEValuesViews::Scalar<deal_II_dimension, deal_II_space_dimension>::
+      get_function_gradients<S>(
+        const ReadVector<S> &,
+        std::vector<
+          ProductType<S, dealii::Tensor<1, deal_II_space_dimension>>::type> &)
+        const;
+
+    template void
+    FEValuesViews::Vector<deal_II_dimension, deal_II_space_dimension>::
+      get_function_gradients<S>(
+        const ReadVector<S> &,
+        std::vector<
+          ProductType<S, dealii::Tensor<2, deal_II_space_dimension>>::type> &)
+        const;
+
+    template void
+    FEValuesViews::Tensor<2, deal_II_dimension, deal_II_space_dimension>::
+      get_function_gradients<S>(
+        const ReadVector<S> &,
+        std::vector<
+          ProductType<S, dealii::Tensor<3, deal_II_space_dimension>>::type> &)
+        const;
+
+    template void FEValuesBase<deal_II_dimension, deal_II_space_dimension>::
+      get_function_gradients<S>(
+        const ReadVector<S> &,
+        std::vector<dealii::Tensor<1, deal_II_space_dimension, S>> &) const;
+
+    template void FEValuesBase<deal_II_dimension, deal_II_space_dimension>::
+      get_function_gradients<S>(
+        const ReadVector<S> &,
+        const ArrayView<const types::global_dof_index> &,
+        std::vector<dealii::Tensor<1, deal_II_space_dimension, S>> &) const;
+
+    template void FEValuesBase<deal_II_dimension, deal_II_space_dimension>::
+      get_function_gradients<S>(
+        const ReadVector<S> &,
+        std::vector<std::vector<dealii::Tensor<1, deal_II_space_dimension, S>>>
+          &) const;
+
+    template void FEValuesBase<deal_II_dimension, deal_II_space_dimension>::
+      get_function_gradients<S>(
+        const ReadVector<S> &,
+        const ArrayView<const types::global_dof_index> &,
+        ArrayView<std::vector<dealii::Tensor<1, deal_II_space_dimension, S>>>,
+        bool) const;
 #  endif
   }
 
@@ -155,14 +203,6 @@ for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS;
   {
 #  if deal_II_dimension <= deal_II_space_dimension
     template void
-    FEValuesViews::Scalar<deal_II_dimension, deal_II_space_dimension>::
-      get_function_gradients<dealii::VEC>(
-        const dealii::VEC &,
-        std::vector<
-          ProductType<dealii::VEC::value_type,
-                      dealii::Tensor<1, deal_II_space_dimension>>::type> &)
-        const;
-    template void
     FEValuesViews::Scalar<deal_II_dimension, deal_II_space_dimension>::
       get_function_hessians<dealii::VEC>(
         const dealii::VEC &,
@@ -186,14 +226,6 @@ for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS;
         const;
 
     template void
-    FEValuesViews::Vector<deal_II_dimension, deal_II_space_dimension>::
-      get_function_gradients<dealii::VEC>(
-        const dealii::VEC &,
-        std::vector<
-          ProductType<dealii::VEC::value_type,
-                      dealii::Tensor<2, 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 &,
@@ -254,14 +286,6 @@ for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS;
           ProductType<dealii::VEC::value_type,
                       dealii::Tensor<1, deal_II_space_dimension>>::type> &)
         const;
-    template void
-    FEValuesViews::Tensor<2, deal_II_dimension, deal_II_space_dimension>::
-      get_function_gradients<dealii::VEC>(
-        const dealii::VEC &,
-        std::vector<
-          ProductType<dealii::VEC::value_type,
-                      dealii::Tensor<3, deal_II_space_dimension>>::type> &)
-        const;
 #  endif
   }
 
@@ -387,32 +411,6 @@ for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS;
   {
 #  if deal_II_dimension <= deal_II_space_dimension
 
-    template void FEValuesBase<deal_II_dimension, deal_II_space_dimension>::
-      get_function_gradients<VEC>(
-        const VEC &,
-        std::vector<dealii::Tensor<1, deal_II_space_dimension, VEC::value_type>>
-          &) const;
-    template void FEValuesBase<deal_II_dimension, deal_II_space_dimension>::
-      get_function_gradients<VEC>(
-        const VEC &,
-        const ArrayView<const types::global_dof_index> &,
-        std::vector<dealii::Tensor<1, deal_II_space_dimension, VEC::value_type>>
-          &) const;
-
-    template void FEValuesBase<deal_II_dimension, deal_II_space_dimension>::
-      get_function_gradients<VEC>(
-        const VEC &,
-        std::vector<std::vector<
-          dealii::Tensor<1, deal_II_space_dimension, VEC::value_type>>> &)
-        const;
-    template void FEValuesBase<deal_II_dimension, deal_II_space_dimension>::
-      get_function_gradients<VEC>(
-        const VEC &,
-        const ArrayView<const types::global_dof_index> &,
-        ArrayView<std::vector<
-          dealii::Tensor<1, deal_II_space_dimension, VEC::value_type>>>,
-        bool) const;
-
     template void FEValuesBase<deal_II_dimension, deal_II_space_dimension>::
       get_function_hessians<VEC>(
         const VEC &,