]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Move the internal FEValuesViews functions into a new file.
authorDavid Wells <drwells@email.unc.edu>
Fri, 18 Aug 2023 15:59:34 +0000 (11:59 -0400)
committerDavid Wells <drwells@email.unc.edu>
Sat, 19 Aug 2023 15:04:03 +0000 (11:04 -0400)
Most of the compilation work for FEValues is actually the inner functions used
by FEValuesViews objects. Explicitly instantiating these functions in their own
translation unit significantly lowers compilation time. Here are the numbers for
release mode:

master:
- fe_values_views: 7:58 wall time, 4.5 GB max RSS

feature:
- fe_values_views: 0:32 wall time, 1.6 GB max RSS
- fe_values_views_internal: 2:58 wall time, 2.8 GB max RSS

I had to adjust the type signatures (on the first ArrayView argument) to avoid
instantiating things for 'const float' rather than float et al.

include/deal.II/fe/fe_values_views.h
include/deal.II/fe/fe_values_views_internal.h [new file with mode: 0644]
source/fe/CMakeLists.txt
source/fe/fe_values_views.cc
source/fe/fe_values_views_internal.cc [new file with mode: 0644]
source/fe/fe_values_views_internal.inst.in [new file with mode: 0644]

index 5c03c9d15370ddc27ee90c04048f0b3ec212900a..b8459e04cf356a7680e902b58e0ad94f25399201 100644 (file)
 
 #include <deal.II/base/config.h>
 
-#include <deal.II/base/derivative_form.h>
 #include <deal.II/base/exceptions.h>
-#include <deal.II/base/point.h>
-#include <deal.II/base/quadrature.h>
-#include <deal.II/base/std_cxx20/iota_view.h>
-#include <deal.II/base/subscriptor.h>
+#include <deal.II/base/smartpointer.h>
 #include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/base/tensor.h>
 
-#include <deal.II/dofs/dof_accessor.h>
-#include <deal.II/dofs/dof_handler.h>
-
-#include <deal.II/fe/fe.h>
 #include <deal.II/fe/fe_update_flags.h>
 #include <deal.II/fe/fe_values_extractors.h>
-#include <deal.II/fe/mapping.h>
-#include <deal.II/fe/mapping_related_data.h>
-
-#include <deal.II/grid/tria.h>
-#include <deal.II/grid/tria_iterator.h>
-
-#include <deal.II/hp/q_collection.h>
 
 #include <deal.II/lac/read_vector.h>
 
-#include <algorithm>
-#include <memory>
 #include <type_traits>
+#include <vector>
 
 DEAL_II_NAMESPACE_OPEN
 
diff --git a/include/deal.II/fe/fe_values_views_internal.h b/include/deal.II/fe/fe_values_views_internal.h
new file mode 100644 (file)
index 0000000..1c53cad
--- /dev/null
@@ -0,0 +1,250 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 2023 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_fe_values_views_internal_h
+#define dealii_fe_values_views_internal_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/array_view.h>
+#include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/base/table.h>
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/fe/fe_values_views.h>
+
+#include <vector>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace FEValuesViews
+{
+  /**
+   * Internal namespace for the utility functions used to actually compute
+   * values by FEValuesViews
+   */
+  namespace internal
+  {
+    // Given values of degrees of freedom, evaluate the
+    // values/gradients/... at quadrature points
+
+    // ------------------------- scalar functions --------------------------
+
+    /**
+     * Compute function values for Scalars.
+     */
+    template <int dim, int spacedim, typename Number>
+    void
+    do_function_values(
+      const ArrayView<const Number> &dof_values,
+      const Table<2, double> &       shape_values,
+      const std::vector<typename Scalar<dim, spacedim>::ShapeFunctionData>
+        &shape_function_data,
+      std::vector<typename ProductType<Number, double>::type> &values);
+
+    /**
+     * same code for gradient and Hessian, template argument 'order' to give
+     * the order of the derivative (= rank of gradient/Hessian tensor)
+     */
+    template <int order, int dim, int spacedim, typename Number>
+    void
+    do_function_derivatives(
+      const ArrayView<const Number> &                  dof_values,
+      const Table<2, dealii::Tensor<order, spacedim>> &shape_derivatives,
+      const std::vector<typename Scalar<dim, spacedim>::ShapeFunctionData>
+        &shape_function_data,
+      std::vector<
+        typename ProductType<Number, dealii::Tensor<order, spacedim>>::type>
+        &derivatives);
+
+    /**
+     * Compute Laplacian values for Scalars.
+     */
+    template <int dim, int spacedim, typename Number>
+    void
+    do_function_laplacians(
+      const ArrayView<const Number> &              dof_values,
+      const Table<2, dealii::Tensor<2, spacedim>> &shape_hessians,
+      const std::vector<typename Scalar<dim, spacedim>::ShapeFunctionData>
+        &shape_function_data,
+      std::vector<typename Scalar<dim, spacedim>::
+                    template solution_laplacian_type<Number>> &laplacians);
+
+    // ----------------------------- vector part ---------------------------
+
+    /**
+     * Compute function values for Vectors.
+     */
+    template <int dim, int spacedim, typename Number>
+    void
+    do_function_values(
+      const ArrayView<const Number> &dof_values,
+      const Table<2, double> &       shape_values,
+      const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
+        &shape_function_data,
+      std::vector<
+        typename ProductType<Number, dealii::Tensor<1, spacedim>>::type>
+        &values);
+
+    /**
+     * same code for gradient and Hessian, template argument 'order' to give
+     * the order of the derivative (= rank of gradient/Hessian tensor)
+     */
+    template <int order, int dim, int spacedim, typename Number>
+    void
+    do_function_derivatives(
+      const ArrayView<const Number> &                  dof_values,
+      const Table<2, dealii::Tensor<order, spacedim>> &shape_derivatives,
+      const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
+        &shape_function_data,
+      std::vector<
+        typename ProductType<Number, dealii::Tensor<order + 1, spacedim>>::type>
+        &derivatives);
+
+    /**
+     * Compute gradient values for Vectors.
+     */
+    template <int dim, int spacedim, typename Number>
+    void
+    do_function_symmetric_gradients(
+      const ArrayView<const Number> &              dof_values,
+      const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients,
+      const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
+        &shape_function_data,
+      std::vector<
+        typename ProductType<Number,
+                             dealii::SymmetricTensor<2, spacedim>>::type>
+        &symmetric_gradients);
+
+    /**
+     * Compute divergence values for Vectors.
+     */
+    template <int dim, int spacedim, typename Number>
+    void
+    do_function_divergences(
+      const ArrayView<const Number> &              dof_values,
+      const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients,
+      const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
+        &shape_function_data,
+      std::vector<typename Vector<dim, spacedim>::
+                    template solution_divergence_type<Number>> &divergences);
+
+    /**
+     * Compute curl values for Vectors.
+     */
+    template <int dim, int spacedim, typename Number>
+    void
+    do_function_curls(
+      const ArrayView<const Number> &              dof_values,
+      const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients,
+      const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
+        &shape_function_data,
+      std::vector<typename ProductType<
+        Number,
+        typename dealii::internal::CurlType<spacedim>::type>::type> &curls);
+
+    /**
+     * Compute Laplacian values for Vectors.
+     */
+    template <int dim, int spacedim, typename Number>
+    void
+    do_function_laplacians(
+      const ArrayView<const Number> &              dof_values,
+      const Table<2, dealii::Tensor<2, spacedim>> &shape_hessians,
+      const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
+        &shape_function_data,
+      std::vector<typename Vector<dim, spacedim>::
+                    template solution_laplacian_type<Number>> &laplacians);
+
+    // ---------------------- symmetric tensor part ------------------------
+
+    /**
+     * Compute values for symmetric tensors.
+     */
+    template <int dim, int spacedim, typename Number>
+    void
+    do_function_values(
+      const ArrayView<const Number> & dof_values,
+      const dealii::Table<2, double> &shape_values,
+      const std::vector<
+        typename SymmetricTensor<2, dim, spacedim>::ShapeFunctionData>
+        &shape_function_data,
+      std::vector<
+        typename ProductType<Number,
+                             dealii::SymmetricTensor<2, spacedim>>::type>
+        &values);
+
+    /**
+     * Compute divergence values for symmetric tensors.
+     */
+    template <int dim, int spacedim, typename Number>
+    void
+    do_function_divergences(
+      const ArrayView<const Number> &              dof_values,
+      const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients,
+      const std::vector<
+        typename SymmetricTensor<2, dim, spacedim>::ShapeFunctionData>
+        &shape_function_data,
+      std::vector<typename SymmetricTensor<2, dim, spacedim>::
+                    template solution_divergence_type<Number>> &divergences);
+
+    // ---------------------- non-symmetric tensor part ------------------------
+
+    /**
+     * Compute values for nonsymmetric tensors.
+     */
+    template <int dim, int spacedim, typename Number>
+    void
+    do_function_values(
+      const ArrayView<const Number> & dof_values,
+      const dealii::Table<2, double> &shape_values,
+      const std::vector<typename Tensor<2, dim, spacedim>::ShapeFunctionData>
+        &shape_function_data,
+      std::vector<
+        typename ProductType<Number, dealii::Tensor<2, spacedim>>::type>
+        &values);
+
+    /**
+     * Compute divergence values for nonsymmetric tensors.
+     */
+    template <int dim, int spacedim, typename Number>
+    void
+    do_function_divergences(
+      const ArrayView<const Number> &              dof_values,
+      const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients,
+      const std::vector<typename Tensor<2, dim, spacedim>::ShapeFunctionData>
+        &shape_function_data,
+      std::vector<typename Tensor<2, dim, spacedim>::
+                    template solution_divergence_type<Number>> &divergences);
+
+    /**
+     * Compute gradient values for nonsymmetric tensors.
+     */
+    template <int dim, int spacedim, typename Number>
+    void
+    do_function_gradients(
+      const ArrayView<const Number> &              dof_values,
+      const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients,
+      const std::vector<typename Tensor<2, dim, spacedim>::ShapeFunctionData>
+        &shape_function_data,
+      std::vector<typename Tensor<2, dim, spacedim>::
+                    template solution_gradient_type<Number>> &gradients);
+  } // end of namespace internal
+} // namespace FEValuesViews
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index 57b6cc663e0c5990291a018ff3c4c0b970607e19..361a293e3523227e439210298b89eada48a1d589 100644 (file)
@@ -70,6 +70,7 @@ set(_separate_src
   fe_values.cc
   fe_values_base.cc
   fe_values_views.cc
+  fe_values_views_internal.cc
   mapping_fe_field.cc
   mapping_fe_field_inst2.cc
   fe_tools.cc
@@ -127,6 +128,7 @@ set(_inst
   fe_trace.inst.in
   fe_values_base.inst.in
   fe_values_views.inst.in
+  fe_values_views_internal.inst.in
   fe_values.inst.in
   fe_wedge_p.inst.in
   mapping_c1.inst.in
index a620da229bf853790f3e354dacbea1784f50b8ae..8326181698304c15a588706e1fa3a6046d50710d 100644 (file)
 // ---------------------------------------------------------------------
 
 #include <deal.II/base/array_view.h>
-#include <deal.II/base/memory_consumption.h>
-#include <deal.II/base/multithread_info.h>
 #include <deal.II/base/numbers.h>
-#include <deal.II/base/quadrature.h>
-#include <deal.II/base/signaling_nan.h>
-#include <deal.II/base/thread_management.h>
 
 #include <deal.II/differentiation/ad.h>
 
-#include <deal.II/dofs/dof_accessor.h>
-
 #include <deal.II/fe/fe.h>
 #include <deal.II/fe/fe_values_base.h>
 #include <deal.II/fe/fe_values_views.h>
-#include <deal.II/fe/mapping.h>
-
-#include <deal.II/grid/tria_accessor.h>
-#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/fe/fe_values_views_internal.h>
 
 #include <deal.II/lac/vector.h>
 
-#include <boost/container/small_vector.hpp>
-
-#include <iomanip>
-#include <memory>
-#include <type_traits>
-
 DEAL_II_NAMESPACE_OPEN
 
 
 namespace internal
 {
-  template <int dim, int spacedim>
-  inline std::vector<unsigned int>
-  make_shape_function_to_row_table(const FiniteElement<dim, spacedim> &fe)
-  {
-    std::vector<unsigned int> shape_function_to_row_table(
-      fe.n_dofs_per_cell() * fe.n_components(), numbers::invalid_unsigned_int);
-    unsigned int row = 0;
-    for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
-      {
-        // loop over all components that are nonzero for this particular
-        // shape function. if a component is zero then we leave the
-        // value in the table unchanged (at the invalid value)
-        // otherwise it is mapped to the next free entry
-        unsigned int nth_nonzero_component = 0;
-        for (unsigned int c = 0; c < fe.n_components(); ++c)
-          if (fe.get_nonzero_components(i)[c] == true)
-            {
-              shape_function_to_row_table[i * fe.n_components() + c] =
-                row + nth_nonzero_component;
-              ++nth_nonzero_component;
-            }
-        row += fe.n_nonzero_components(i);
-      }
-
-    return shape_function_to_row_table;
-  }
-
   namespace
   {
-    // Check to see if a DoF value is zero, implying that subsequent operations
-    // with the value have no effect.
-    template <typename Number, typename T = void>
-    struct CheckForZero
-    {
-      static bool
-      value(const Number &value)
-      {
-        return value == dealii::internal::NumberType<Number>::value(0.0);
-      }
-    };
-
-    // For auto-differentiable numbers, the fact that a DoF value is zero
-    // does not imply that its derivatives are zero as well. So we
-    // can't filter by value for these number types.
-    // Note that we also want to avoid actually checking the value itself,
-    // since some AD numbers are not contextually convertible to booleans.
-    template <typename Number>
-    struct CheckForZero<
-      Number,
-      std::enable_if_t<Differentiation::AD::is_ad_number<Number>::value>>
+    template <int dim, int spacedim>
+    inline std::vector<unsigned int>
+    make_shape_function_to_row_table(const FiniteElement<dim, spacedim> &fe)
     {
-      static bool
-      value(const Number & /*value*/)
-      {
-        return false;
-      }
-    };
+      std::vector<unsigned int> shape_function_to_row_table(
+        fe.n_dofs_per_cell() * fe.n_components(),
+        numbers::invalid_unsigned_int);
+      unsigned int row = 0;
+      for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
+        {
+          // loop over all components that are nonzero for this particular
+          // shape function. if a component is zero then we leave the
+          // value in the table unchanged (at the invalid value)
+          // otherwise it is mapped to the next free entry
+          unsigned int nth_nonzero_component = 0;
+          for (unsigned int c = 0; c < fe.n_components(); ++c)
+            if (fe.get_nonzero_components(i)[c] == true)
+              {
+                shape_function_to_row_table[i * fe.n_components() + c] =
+                  row + nth_nonzero_component;
+                ++nth_nonzero_component;
+              }
+          row += fe.n_nonzero_components(i);
+        }
+
+      return shape_function_to_row_table;
+    }
   } // namespace
 } // namespace internal
 
@@ -400,1116 +356,6 @@ namespace FEValuesViews
 
 
 
-  namespace internal
-  {
-    // Given values of degrees of freedom, evaluate the
-    // values/gradients/... at quadrature points
-
-    // ------------------------- scalar functions --------------------------
-    template <int dim, int spacedim, typename Number>
-    void
-    do_function_values(
-      const ArrayView<Number> &dof_values,
-      const Table<2, double> & shape_values,
-      const std::vector<typename Scalar<dim, spacedim>::ShapeFunctionData>
-        &shape_function_data,
-      std::vector<typename ProductType<Number, double>::type> &values)
-    {
-      const unsigned int dofs_per_cell       = dof_values.size();
-      const unsigned int n_quadrature_points = values.size();
-
-      std::fill(values.begin(),
-                values.end(),
-                dealii::internal::NumberType<Number>::value(0.0));
-
-      for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
-           ++shape_function)
-        if (shape_function_data[shape_function]
-              .is_nonzero_shape_function_component)
-          {
-            const Number &value = dof_values[shape_function];
-            // For auto-differentiable numbers, the fact that a DoF value is
-            // zero does not imply that its derivatives are zero as well. So we
-            // can't filter by value for these number types.
-            if (dealii::internal::CheckForZero<Number>::value(value) == true)
-              continue;
-
-            const double *shape_value_ptr =
-              &shape_values(shape_function_data[shape_function].row_index, 0);
-            for (unsigned int q_point = 0; q_point < n_quadrature_points;
-                 ++q_point)
-              values[q_point] += value * (*shape_value_ptr++);
-          }
-    }
-
-
-
-    // same code for gradient and Hessian, template argument 'order' to give
-    // the order of the derivative (= rank of gradient/Hessian tensor)
-    template <int order, int dim, int spacedim, typename Number>
-    void
-    do_function_derivatives(
-      const ArrayView<Number> &                        dof_values,
-      const Table<2, dealii::Tensor<order, spacedim>> &shape_derivatives,
-      const std::vector<typename Scalar<dim, spacedim>::ShapeFunctionData>
-        &shape_function_data,
-      std::vector<
-        typename ProductType<Number, dealii::Tensor<order, spacedim>>::type>
-        &derivatives)
-    {
-      const unsigned int dofs_per_cell       = dof_values.size();
-      const unsigned int n_quadrature_points = derivatives.size();
-
-      std::fill(
-        derivatives.begin(),
-        derivatives.end(),
-        typename ProductType<Number, dealii::Tensor<order, spacedim>>::type());
-
-      for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
-           ++shape_function)
-        if (shape_function_data[shape_function]
-              .is_nonzero_shape_function_component)
-          {
-            const Number &value = dof_values[shape_function];
-            // For auto-differentiable numbers, the fact that a DoF value is
-            // zero does not imply that its derivatives are zero as well. So we
-            // can't filter by value for these number types.
-            if (dealii::internal::CheckForZero<Number>::value(value) == true)
-              continue;
-
-            const dealii::Tensor<order, spacedim> *shape_derivative_ptr =
-              &shape_derivatives[shape_function_data[shape_function].row_index]
-                                [0];
-            for (unsigned int q_point = 0; q_point < n_quadrature_points;
-                 ++q_point)
-              derivatives[q_point] += value * (*shape_derivative_ptr++);
-          }
-    }
-
-
-
-    template <int dim, int spacedim, typename Number>
-    void
-    do_function_laplacians(
-      const ArrayView<Number> &                    dof_values,
-      const Table<2, dealii::Tensor<2, spacedim>> &shape_hessians,
-      const std::vector<typename Scalar<dim, spacedim>::ShapeFunctionData>
-        &shape_function_data,
-      std::vector<typename Scalar<dim, spacedim>::
-                    template solution_laplacian_type<Number>> &laplacians)
-    {
-      const unsigned int dofs_per_cell       = dof_values.size();
-      const unsigned int n_quadrature_points = laplacians.size();
-
-      std::fill(
-        laplacians.begin(),
-        laplacians.end(),
-        typename Scalar<dim,
-                        spacedim>::template solution_laplacian_type<Number>());
-
-      for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
-           ++shape_function)
-        if (shape_function_data[shape_function]
-              .is_nonzero_shape_function_component)
-          {
-            const Number &value = dof_values[shape_function];
-            // For auto-differentiable numbers, the fact that a DoF value is
-            // zero does not imply that its derivatives are zero as well. So we
-            // can't filter by value for these number types.
-            if (dealii::internal::CheckForZero<Number>::value(value) == true)
-              continue;
-
-            const dealii::Tensor<2, spacedim> *shape_hessian_ptr =
-              &shape_hessians[shape_function_data[shape_function].row_index][0];
-            for (unsigned int q_point = 0; q_point < n_quadrature_points;
-                 ++q_point)
-              laplacians[q_point] += value * trace(*shape_hessian_ptr++);
-          }
-    }
-
-
-
-    // ----------------------------- vector part ---------------------------
-
-    template <int dim, int spacedim, typename Number>
-    void
-    do_function_values(
-      const ArrayView<Number> &dof_values,
-      const Table<2, double> & shape_values,
-      const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
-        &shape_function_data,
-      std::vector<
-        typename ProductType<Number, dealii::Tensor<1, spacedim>>::type>
-        &values)
-    {
-      const unsigned int dofs_per_cell       = dof_values.size();
-      const unsigned int n_quadrature_points = values.size();
-
-      std::fill(
-        values.begin(),
-        values.end(),
-        typename ProductType<Number, dealii::Tensor<1, spacedim>>::type());
-
-      for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
-           ++shape_function)
-        {
-          const int snc =
-            shape_function_data[shape_function].single_nonzero_component;
-
-          if (snc == -2)
-            // shape function is zero for the selected components
-            continue;
-
-          const Number &value = dof_values[shape_function];
-          // For auto-differentiable numbers, the fact that a DoF value is zero
-          // does not imply that its derivatives are zero as well. So we
-          // can't filter by value for these number types.
-          if (dealii::internal::CheckForZero<Number>::value(value) == true)
-            continue;
-
-          if (snc != -1)
-            {
-              const unsigned int comp = shape_function_data[shape_function]
-                                          .single_nonzero_component_index;
-              const double *shape_value_ptr = &shape_values(snc, 0);
-              for (unsigned int q_point = 0; q_point < n_quadrature_points;
-                   ++q_point)
-                values[q_point][comp] += value * (*shape_value_ptr++);
-            }
-          else
-            for (unsigned int d = 0; d < spacedim; ++d)
-              if (shape_function_data[shape_function]
-                    .is_nonzero_shape_function_component[d])
-                {
-                  const double *shape_value_ptr = &shape_values(
-                    shape_function_data[shape_function].row_index[d], 0);
-                  for (unsigned int q_point = 0; q_point < n_quadrature_points;
-                       ++q_point)
-                    values[q_point][d] += value * (*shape_value_ptr++);
-                }
-        }
-    }
-
-
-
-    template <int order, int dim, int spacedim, typename Number>
-    void
-    do_function_derivatives(
-      const ArrayView<Number> &                        dof_values,
-      const Table<2, dealii::Tensor<order, spacedim>> &shape_derivatives,
-      const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
-        &shape_function_data,
-      std::vector<
-        typename ProductType<Number, dealii::Tensor<order + 1, spacedim>>::type>
-        &derivatives)
-    {
-      const unsigned int dofs_per_cell       = dof_values.size();
-      const unsigned int n_quadrature_points = derivatives.size();
-
-      std::fill(
-        derivatives.begin(),
-        derivatives.end(),
-        typename ProductType<Number,
-                             dealii::Tensor<order + 1, spacedim>>::type());
-
-      for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
-           ++shape_function)
-        {
-          const int snc =
-            shape_function_data[shape_function].single_nonzero_component;
-
-          if (snc == -2)
-            // shape function is zero for the selected components
-            continue;
-
-          const Number &value = dof_values[shape_function];
-          // For auto-differentiable numbers, the fact that a DoF value is zero
-          // does not imply that its derivatives are zero as well. So we
-          // can't filter by value for these number types.
-          if (dealii::internal::CheckForZero<Number>::value(value) == true)
-            continue;
-
-          if (snc != -1)
-            {
-              const unsigned int comp = shape_function_data[shape_function]
-                                          .single_nonzero_component_index;
-              const dealii::Tensor<order, spacedim> *shape_derivative_ptr =
-                &shape_derivatives[snc][0];
-              for (unsigned int q_point = 0; q_point < n_quadrature_points;
-                   ++q_point)
-                derivatives[q_point][comp] += value * (*shape_derivative_ptr++);
-            }
-          else
-            for (unsigned int d = 0; d < spacedim; ++d)
-              if (shape_function_data[shape_function]
-                    .is_nonzero_shape_function_component[d])
-                {
-                  const dealii::Tensor<order, spacedim> *shape_derivative_ptr =
-                    &shape_derivatives[shape_function_data[shape_function]
-                                         .row_index[d]][0];
-                  for (unsigned int q_point = 0; q_point < n_quadrature_points;
-                       ++q_point)
-                    derivatives[q_point][d] +=
-                      value * (*shape_derivative_ptr++);
-                }
-        }
-    }
-
-
-
-    template <int dim, int spacedim, typename Number>
-    void
-    do_function_symmetric_gradients(
-      const ArrayView<Number> &                    dof_values,
-      const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients,
-      const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
-        &shape_function_data,
-      std::vector<
-        typename ProductType<Number,
-                             dealii::SymmetricTensor<2, spacedim>>::type>
-        &symmetric_gradients)
-    {
-      const unsigned int dofs_per_cell       = dof_values.size();
-      const unsigned int n_quadrature_points = symmetric_gradients.size();
-
-      std::fill(
-        symmetric_gradients.begin(),
-        symmetric_gradients.end(),
-        typename ProductType<Number,
-                             dealii::SymmetricTensor<2, spacedim>>::type());
-
-      for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
-           ++shape_function)
-        {
-          const int snc =
-            shape_function_data[shape_function].single_nonzero_component;
-
-          if (snc == -2)
-            // shape function is zero for the selected components
-            continue;
-
-          const Number &value = dof_values[shape_function];
-          // For auto-differentiable numbers, the fact that a DoF value is zero
-          // does not imply that its derivatives are zero as well. So we
-          // can't filter by value for these number types.
-          if (dealii::internal::CheckForZero<Number>::value(value) == true)
-            continue;
-
-          if (snc != -1)
-            {
-              const unsigned int comp = shape_function_data[shape_function]
-                                          .single_nonzero_component_index;
-              const dealii::Tensor<1, spacedim> *shape_gradient_ptr =
-                &shape_gradients[snc][0];
-              for (unsigned int q_point = 0; q_point < n_quadrature_points;
-                   ++q_point)
-                symmetric_gradients[q_point] +=
-                  value * dealii::SymmetricTensor<2, spacedim>(
-                            symmetrize_single_row(comp, *shape_gradient_ptr++));
-            }
-          else
-            for (unsigned int q_point = 0; q_point < n_quadrature_points;
-                 ++q_point)
-              {
-                typename ProductType<Number, dealii::Tensor<2, spacedim>>::type
-                  grad;
-                for (unsigned int d = 0; d < spacedim; ++d)
-                  if (shape_function_data[shape_function]
-                        .is_nonzero_shape_function_component[d])
-                    grad[d] =
-                      value *
-                      shape_gradients[shape_function_data[shape_function]
-                                        .row_index[d]][q_point];
-                symmetric_gradients[q_point] += symmetrize(grad);
-              }
-        }
-    }
-
-
-
-    template <int dim, int spacedim, typename Number>
-    void
-    do_function_divergences(
-      const ArrayView<Number> &                    dof_values,
-      const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients,
-      const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
-        &shape_function_data,
-      std::vector<typename Vector<dim, spacedim>::
-                    template solution_divergence_type<Number>> &divergences)
-    {
-      const unsigned int dofs_per_cell       = dof_values.size();
-      const unsigned int n_quadrature_points = divergences.size();
-
-      std::fill(
-        divergences.begin(),
-        divergences.end(),
-        typename Vector<dim,
-                        spacedim>::template solution_divergence_type<Number>());
-
-      for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
-           ++shape_function)
-        {
-          const int snc =
-            shape_function_data[shape_function].single_nonzero_component;
-
-          if (snc == -2)
-            // shape function is zero for the selected components
-            continue;
-
-          const Number &value = dof_values[shape_function];
-          // For auto-differentiable numbers, the fact that a DoF value is zero
-          // does not imply that its derivatives are zero as well. So we
-          // can't filter by value for these number types.
-          if (dealii::internal::CheckForZero<Number>::value(value) == true)
-            continue;
-
-          if (snc != -1)
-            {
-              const unsigned int comp = shape_function_data[shape_function]
-                                          .single_nonzero_component_index;
-              const dealii::Tensor<1, spacedim> *shape_gradient_ptr =
-                &shape_gradients[snc][0];
-              for (unsigned int q_point = 0; q_point < n_quadrature_points;
-                   ++q_point)
-                divergences[q_point] += value * (*shape_gradient_ptr++)[comp];
-            }
-          else
-            for (unsigned int d = 0; d < spacedim; ++d)
-              if (shape_function_data[shape_function]
-                    .is_nonzero_shape_function_component[d])
-                {
-                  const dealii::Tensor<1, spacedim> *shape_gradient_ptr =
-                    &shape_gradients[shape_function_data[shape_function]
-                                       .row_index[d]][0];
-                  for (unsigned int q_point = 0; q_point < n_quadrature_points;
-                       ++q_point)
-                    divergences[q_point] += value * (*shape_gradient_ptr++)[d];
-                }
-        }
-    }
-
-
-
-    template <int dim, int spacedim, typename Number>
-    void
-    do_function_curls(
-      const ArrayView<Number> &                    dof_values,
-      const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients,
-      const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
-        &shape_function_data,
-      std::vector<typename ProductType<
-        Number,
-        typename dealii::internal::CurlType<spacedim>::type>::type> &curls)
-    {
-      const unsigned int dofs_per_cell       = dof_values.size();
-      const unsigned int n_quadrature_points = curls.size();
-
-      std::fill(curls.begin(),
-                curls.end(),
-                typename ProductType<
-                  Number,
-                  typename dealii::internal::CurlType<spacedim>::type>::type());
-
-      switch (spacedim)
-        {
-          case 1:
-            {
-              Assert(false,
-                     ExcMessage(
-                       "Computing the curl in 1d is not a useful operation"));
-              break;
-            }
-
-          case 2:
-            {
-              for (unsigned int shape_function = 0;
-                   shape_function < dofs_per_cell;
-                   ++shape_function)
-                {
-                  const int snc = shape_function_data[shape_function]
-                                    .single_nonzero_component;
-
-                  if (snc == -2)
-                    // shape function is zero for the selected components
-                    continue;
-
-                  const Number &value = dof_values[shape_function];
-                  // For auto-differentiable numbers, the fact that a DoF value
-                  // is zero does not imply that its derivatives are zero as
-                  // well. So we can't filter by value for these number types.
-                  if (dealii::internal::CheckForZero<Number>::value(value) ==
-                      true)
-                    continue;
-
-                  if (snc != -1)
-                    {
-                      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];
-                    }
-                  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])
-                        {
-                          const dealii::Tensor<1,
-                                               spacedim> *shape_gradient_ptr =
-                            &shape_gradients[shape_function_data[shape_function]
-                                               .row_index[0]][0];
-
-                          for (unsigned int q_point = 0;
-                               q_point < n_quadrature_points;
-                               ++q_point)
-                            curls[q_point][0] -=
-                              value * (*shape_gradient_ptr++)[1];
-                        }
-
-                      if (shape_function_data[shape_function]
-                            .is_nonzero_shape_function_component[1])
-                        {
-                          const dealii::Tensor<1,
-                                               spacedim> *shape_gradient_ptr =
-                            &shape_gradients[shape_function_data[shape_function]
-                                               .row_index[1]][0];
-
-                          for (unsigned int q_point = 0;
-                               q_point < n_quadrature_points;
-                               ++q_point)
-                            curls[q_point][0] +=
-                              value * (*shape_gradient_ptr++)[0];
-                        }
-                    }
-                }
-              break;
-            }
-
-          case 3:
-            {
-              for (unsigned int shape_function = 0;
-                   shape_function < dofs_per_cell;
-                   ++shape_function)
-                {
-                  const int snc = shape_function_data[shape_function]
-                                    .single_nonzero_component;
-
-                  if (snc == -2)
-                    // shape function is zero for the selected components
-                    continue;
-
-                  const Number &value = dof_values[shape_function];
-                  // For auto-differentiable numbers, the fact that a DoF value
-                  // is zero does not imply that its derivatives are zero as
-                  // well. So we can't filter by value for these number types.
-                  if (dealii::internal::CheckForZero<Number>::value(value) ==
-                      true)
-                    continue;
-
-                  if (snc != -1)
-                    {
-                      const dealii::Tensor<1, spacedim> *shape_gradient_ptr =
-                        &shape_gradients[snc][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][1] +=
-                                    value * (*shape_gradient_ptr)[2];
-                                  curls[q_point][2] -=
-                                    value * (*shape_gradient_ptr++)[1];
-                                }
-
-                              break;
-                            }
-
-                          case 1:
-                            {
-                              for (unsigned int q_point = 0;
-                                   q_point < n_quadrature_points;
-                                   ++q_point)
-                                {
-                                  curls[q_point][0] -=
-                                    value * (*shape_gradient_ptr)[2];
-                                  curls[q_point][2] +=
-                                    value * (*shape_gradient_ptr++)[0];
-                                }
-
-                              break;
-                            }
-
-                          case 2:
-                            {
-                              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])
-                        {
-                          const dealii::Tensor<1,
-                                               spacedim> *shape_gradient_ptr =
-                            &shape_gradients[shape_function_data[shape_function]
-                                               .row_index[0]][0];
-
-                          for (unsigned int q_point = 0;
-                               q_point < n_quadrature_points;
-                               ++q_point)
-                            {
-                              curls[q_point][1] +=
-                                value * (*shape_gradient_ptr)[2];
-                              curls[q_point][2] -=
-                                value * (*shape_gradient_ptr++)[1];
-                            }
-                        }
-
-                      if (shape_function_data[shape_function]
-                            .is_nonzero_shape_function_component[1])
-                        {
-                          const dealii::Tensor<1,
-                                               spacedim> *shape_gradient_ptr =
-                            &shape_gradients[shape_function_data[shape_function]
-                                               .row_index[1]][0];
-
-                          for (unsigned int q_point = 0;
-                               q_point < n_quadrature_points;
-                               ++q_point)
-                            {
-                              curls[q_point][0] -=
-                                value * (*shape_gradient_ptr)[2];
-                              curls[q_point][2] +=
-                                value * (*shape_gradient_ptr++)[0];
-                            }
-                        }
-
-                      if (shape_function_data[shape_function]
-                            .is_nonzero_shape_function_component[2])
-                        {
-                          const dealii::Tensor<1,
-                                               spacedim> *shape_gradient_ptr =
-                            &shape_gradients[shape_function_data[shape_function]
-                                               .row_index[2]][0];
-
-                          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];
-                            }
-                        }
-                    }
-                }
-            }
-        }
-    }
-
-
-
-    template <int dim, int spacedim, typename Number>
-    void
-    do_function_laplacians(
-      const ArrayView<Number> &                    dof_values,
-      const Table<2, dealii::Tensor<2, spacedim>> &shape_hessians,
-      const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
-        &shape_function_data,
-      std::vector<typename Vector<dim, spacedim>::
-                    template solution_laplacian_type<Number>> &laplacians)
-    {
-      const unsigned int dofs_per_cell       = dof_values.size();
-      const unsigned int n_quadrature_points = laplacians.size();
-
-      std::fill(
-        laplacians.begin(),
-        laplacians.end(),
-        typename Vector<dim,
-                        spacedim>::template solution_laplacian_type<Number>());
-
-      for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
-           ++shape_function)
-        {
-          const int snc =
-            shape_function_data[shape_function].single_nonzero_component;
-
-          if (snc == -2)
-            // shape function is zero for the selected components
-            continue;
-
-          const Number &value = dof_values[shape_function];
-          // For auto-differentiable numbers, the fact that a DoF value is zero
-          // does not imply that its derivatives are zero as well. So we
-          // can't filter by value for these number types.
-          if (dealii::internal::CheckForZero<Number>::value(value) == true)
-            continue;
-
-          if (snc != -1)
-            {
-              const unsigned int comp = shape_function_data[shape_function]
-                                          .single_nonzero_component_index;
-              const dealii::Tensor<2, spacedim> *shape_hessian_ptr =
-                &shape_hessians[snc][0];
-              for (unsigned int q_point = 0; q_point < n_quadrature_points;
-                   ++q_point)
-                laplacians[q_point][comp] +=
-                  value * trace(*shape_hessian_ptr++);
-            }
-          else
-            for (unsigned int d = 0; d < spacedim; ++d)
-              if (shape_function_data[shape_function]
-                    .is_nonzero_shape_function_component[d])
-                {
-                  const dealii::Tensor<2, spacedim> *shape_hessian_ptr =
-                    &shape_hessians[shape_function_data[shape_function]
-                                      .row_index[d]][0];
-                  for (unsigned int q_point = 0; q_point < n_quadrature_points;
-                       ++q_point)
-                    laplacians[q_point][d] +=
-                      value * trace(*shape_hessian_ptr++);
-                }
-        }
-    }
-
-
-
-    // ---------------------- symmetric tensor part ------------------------
-
-    template <int dim, int spacedim, typename Number>
-    void
-    do_function_values(
-      const ArrayView<Number> &       dof_values,
-      const dealii::Table<2, double> &shape_values,
-      const std::vector<
-        typename SymmetricTensor<2, dim, spacedim>::ShapeFunctionData>
-        &shape_function_data,
-      std::vector<
-        typename ProductType<Number,
-                             dealii::SymmetricTensor<2, spacedim>>::type>
-        &values)
-    {
-      const unsigned int dofs_per_cell       = dof_values.size();
-      const unsigned int n_quadrature_points = values.size();
-
-      std::fill(
-        values.begin(),
-        values.end(),
-        typename ProductType<Number,
-                             dealii::SymmetricTensor<2, spacedim>>::type());
-
-      for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
-           ++shape_function)
-        {
-          const int snc =
-            shape_function_data[shape_function].single_nonzero_component;
-
-          if (snc == -2)
-            // shape function is zero for the selected components
-            continue;
-
-          const Number &value = dof_values[shape_function];
-          // For auto-differentiable numbers, the fact that a DoF value is zero
-          // does not imply that its derivatives are zero as well. So we
-          // can't filter by value for these number types.
-          if (dealii::internal::CheckForZero<Number>::value(value) == true)
-            continue;
-
-          if (snc != -1)
-            {
-              const TableIndices<2> comp = dealii::
-                SymmetricTensor<2, spacedim>::unrolled_to_component_indices(
-                  shape_function_data[shape_function]
-                    .single_nonzero_component_index);
-              const double *shape_value_ptr = &shape_values(snc, 0);
-              for (unsigned int q_point = 0; q_point < n_quadrature_points;
-                   ++q_point)
-                values[q_point][comp] += value * (*shape_value_ptr++);
-            }
-          else
-            for (unsigned int d = 0;
-                 d <
-                 dealii::SymmetricTensor<2, spacedim>::n_independent_components;
-                 ++d)
-              if (shape_function_data[shape_function]
-                    .is_nonzero_shape_function_component[d])
-                {
-                  const TableIndices<2> comp =
-                    dealii::SymmetricTensor<2, spacedim>::
-                      unrolled_to_component_indices(d);
-                  const double *shape_value_ptr = &shape_values(
-                    shape_function_data[shape_function].row_index[d], 0);
-                  for (unsigned int q_point = 0; q_point < n_quadrature_points;
-                       ++q_point)
-                    values[q_point][comp] += value * (*shape_value_ptr++);
-                }
-        }
-    }
-
-
-
-    template <int dim, int spacedim, typename Number>
-    void
-    do_function_divergences(
-      const ArrayView<Number> &                    dof_values,
-      const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients,
-      const std::vector<
-        typename SymmetricTensor<2, dim, spacedim>::ShapeFunctionData>
-        &shape_function_data,
-      std::vector<typename SymmetricTensor<2, dim, spacedim>::
-                    template solution_divergence_type<Number>> &divergences)
-    {
-      const unsigned int dofs_per_cell       = dof_values.size();
-      const unsigned int n_quadrature_points = divergences.size();
-
-      std::fill(divergences.begin(),
-                divergences.end(),
-                typename SymmetricTensor<2, dim, spacedim>::
-                  template solution_divergence_type<Number>());
-
-      for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
-           ++shape_function)
-        {
-          const int snc =
-            shape_function_data[shape_function].single_nonzero_component;
-
-          if (snc == -2)
-            // shape function is zero for the selected components
-            continue;
-
-          const Number &value = dof_values[shape_function];
-          // For auto-differentiable numbers, the fact that a DoF value is zero
-          // does not imply that its derivatives are zero as well. So we
-          // can't filter by value for these number types.
-          if (dealii::internal::CheckForZero<Number>::value(value) == true)
-            continue;
-
-          if (snc != -1)
-            {
-              const unsigned int comp = shape_function_data[shape_function]
-                                          .single_nonzero_component_index;
-
-              const dealii::Tensor<1, spacedim> *shape_gradient_ptr =
-                &shape_gradients[snc][0];
-
-              const unsigned int ii = dealii::SymmetricTensor<2, spacedim>::
-                unrolled_to_component_indices(comp)[0];
-              const unsigned int jj = dealii::SymmetricTensor<2, spacedim>::
-                unrolled_to_component_indices(comp)[1];
-
-              for (unsigned int q_point = 0; q_point < n_quadrature_points;
-                   ++q_point, ++shape_gradient_ptr)
-                {
-                  divergences[q_point][ii] += value * (*shape_gradient_ptr)[jj];
-
-                  if (ii != jj)
-                    divergences[q_point][jj] +=
-                      value * (*shape_gradient_ptr)[ii];
-                }
-            }
-          else
-            {
-              for (unsigned int d = 0;
-                   d <
-                   dealii::SymmetricTensor<2,
-                                           spacedim>::n_independent_components;
-                   ++d)
-                if (shape_function_data[shape_function]
-                      .is_nonzero_shape_function_component[d])
-                  {
-                    Assert(false, ExcNotImplemented());
-
-                    // the following implementation needs to be looked over -- I
-                    // think it can't be right, because we are in a case where
-                    // there is no single nonzero component
-                    //
-                    // the following is not implemented! we need to consider the
-                    // interplay between multiple non-zero entries in shape
-                    // function and the representation as a symmetric
-                    // second-order tensor
-                    const unsigned int comp =
-                      shape_function_data[shape_function]
-                        .single_nonzero_component_index;
-
-                    const dealii::Tensor<1, spacedim> *shape_gradient_ptr =
-                      &shape_gradients[shape_function_data[shape_function]
-                                         .row_index[d]][0];
-                    for (unsigned int q_point = 0;
-                         q_point < n_quadrature_points;
-                         ++q_point, ++shape_gradient_ptr)
-                      {
-                        for (unsigned int j = 0; j < spacedim; ++j)
-                          {
-                            const unsigned int vector_component =
-                              dealii::SymmetricTensor<2, spacedim>::
-                                component_to_unrolled_index(
-                                  TableIndices<2>(comp, j));
-                            divergences[q_point][vector_component] +=
-                              value * (*shape_gradient_ptr++)[j];
-                          }
-                      }
-                  }
-            }
-        }
-    }
-
-    // ---------------------- non-symmetric tensor part ------------------------
-
-    template <int dim, int spacedim, typename Number>
-    void
-    do_function_values(
-      const ArrayView<Number> &       dof_values,
-      const dealii::Table<2, double> &shape_values,
-      const std::vector<typename Tensor<2, dim, spacedim>::ShapeFunctionData>
-        &shape_function_data,
-      std::vector<
-        typename ProductType<Number, dealii::Tensor<2, spacedim>>::type>
-        &values)
-    {
-      const unsigned int dofs_per_cell       = dof_values.size();
-      const unsigned int n_quadrature_points = values.size();
-
-      std::fill(
-        values.begin(),
-        values.end(),
-        typename ProductType<Number, dealii::Tensor<2, spacedim>>::type());
-
-      for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
-           ++shape_function)
-        {
-          const int snc =
-            shape_function_data[shape_function].single_nonzero_component;
-
-          if (snc == -2)
-            // shape function is zero for the selected components
-            continue;
-
-          const Number &value = dof_values[shape_function];
-          // For auto-differentiable numbers, the fact that a DoF value is zero
-          // does not imply that its derivatives are zero as well. So we
-          // can't filter by value for these number types.
-          if (dealii::internal::CheckForZero<Number>::value(value) == true)
-            continue;
-
-          if (snc != -1)
-            {
-              const unsigned int comp = shape_function_data[shape_function]
-                                          .single_nonzero_component_index;
-
-              const TableIndices<2> indices =
-                dealii::Tensor<2, spacedim>::unrolled_to_component_indices(
-                  comp);
-
-              const double *shape_value_ptr = &shape_values(snc, 0);
-              for (unsigned int q_point = 0; q_point < n_quadrature_points;
-                   ++q_point)
-                values[q_point][indices] += value * (*shape_value_ptr++);
-            }
-          else
-            for (unsigned int d = 0; d < dim * dim; ++d)
-              if (shape_function_data[shape_function]
-                    .is_nonzero_shape_function_component[d])
-                {
-                  const TableIndices<2> indices =
-                    dealii::Tensor<2, spacedim>::unrolled_to_component_indices(
-                      d);
-
-                  const double *shape_value_ptr = &shape_values(
-                    shape_function_data[shape_function].row_index[d], 0);
-                  for (unsigned int q_point = 0; q_point < n_quadrature_points;
-                       ++q_point)
-                    values[q_point][indices] += value * (*shape_value_ptr++);
-                }
-        }
-    }
-
-
-
-    template <int dim, int spacedim, typename Number>
-    void
-    do_function_divergences(
-      const ArrayView<Number> &                    dof_values,
-      const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients,
-      const std::vector<typename Tensor<2, dim, spacedim>::ShapeFunctionData>
-        &shape_function_data,
-      std::vector<typename Tensor<2, dim, spacedim>::
-                    template solution_divergence_type<Number>> &divergences)
-    {
-      const unsigned int dofs_per_cell       = dof_values.size();
-      const unsigned int n_quadrature_points = divergences.size();
-
-      std::fill(
-        divergences.begin(),
-        divergences.end(),
-        typename Tensor<2, dim, spacedim>::template solution_divergence_type<
-          Number>());
-
-      for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
-           ++shape_function)
-        {
-          const int snc =
-            shape_function_data[shape_function].single_nonzero_component;
-
-          if (snc == -2)
-            // shape function is zero for the selected components
-            continue;
-
-          const Number &value = dof_values[shape_function];
-          // For auto-differentiable numbers, the fact that a DoF value is zero
-          // does not imply that its derivatives are zero as well. So we
-          // can't filter by value for these number types.
-          if (dealii::internal::CheckForZero<Number>::value(value) == true)
-            continue;
-
-          if (snc != -1)
-            {
-              const unsigned int comp = shape_function_data[shape_function]
-                                          .single_nonzero_component_index;
-
-              const dealii::Tensor<1, spacedim> *shape_gradient_ptr =
-                &shape_gradients[snc][0];
-
-              const TableIndices<2> indices =
-                dealii::Tensor<2, spacedim>::unrolled_to_component_indices(
-                  comp);
-              const unsigned int ii = indices[0];
-              const unsigned int jj = indices[1];
-
-              for (unsigned int q_point = 0; q_point < n_quadrature_points;
-                   ++q_point, ++shape_gradient_ptr)
-                {
-                  divergences[q_point][ii] += value * (*shape_gradient_ptr)[jj];
-                }
-            }
-          else
-            {
-              for (unsigned int d = 0; d < dim * dim; ++d)
-                if (shape_function_data[shape_function]
-                      .is_nonzero_shape_function_component[d])
-                  {
-                    Assert(false, ExcNotImplemented());
-                  }
-            }
-        }
-    }
-
-
-
-    template <int dim, int spacedim, typename Number>
-    void
-    do_function_gradients(
-      const ArrayView<Number> &                    dof_values,
-      const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients,
-      const std::vector<typename Tensor<2, dim, spacedim>::ShapeFunctionData>
-        &shape_function_data,
-      std::vector<typename Tensor<2, dim, spacedim>::
-                    template solution_gradient_type<Number>> &gradients)
-    {
-      const unsigned int dofs_per_cell       = dof_values.size();
-      const unsigned int n_quadrature_points = gradients.size();
-
-      std::fill(
-        gradients.begin(),
-        gradients.end(),
-        typename Tensor<2, dim, spacedim>::template solution_gradient_type<
-          Number>());
-
-      for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
-           ++shape_function)
-        {
-          const int snc =
-            shape_function_data[shape_function].single_nonzero_component;
-
-          if (snc == -2)
-            // shape function is zero for the selected components
-            continue;
-
-          const Number &value = dof_values[shape_function];
-          // For auto-differentiable numbers, the fact that a DoF value is zero
-          // does not imply that its derivatives are zero as well. So we
-          // can't filter by value for these number types.
-          if (dealii::internal::CheckForZero<Number>::value(value) == true)
-            continue;
-
-          if (snc != -1)
-            {
-              const unsigned int comp = shape_function_data[shape_function]
-                                          .single_nonzero_component_index;
-
-              const dealii::Tensor<1, spacedim> *shape_gradient_ptr =
-                &shape_gradients[snc][0];
-
-              const TableIndices<2> indices =
-                dealii::Tensor<2, spacedim>::unrolled_to_component_indices(
-                  comp);
-              const unsigned int ii = indices[0];
-              const unsigned int jj = indices[1];
-
-              for (unsigned int q_point = 0; q_point < n_quadrature_points;
-                   ++q_point, ++shape_gradient_ptr)
-                {
-                  gradients[q_point][ii][jj] += value * (*shape_gradient_ptr);
-                }
-            }
-          else
-            {
-              for (unsigned int d = 0; d < dim * dim; ++d)
-                if (shape_function_data[shape_function]
-                      .is_nonzero_shape_function_component[d])
-                  {
-                    Assert(false, ExcNotImplemented());
-                  }
-            }
-        }
-    }
-
-  } // end of namespace internal
-
-
-
   template <int dim, int spacedim>
   template <typename Number>
   void
@@ -1529,7 +375,7 @@ namespace FEValuesViews
     fe_values->present_cell.get_interpolated_dof_values(fe_function,
                                                         dof_values);
     internal::do_function_values<dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_values,
       shape_function_data,
       values);
@@ -1553,7 +399,7 @@ namespace FEValuesViews
     AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
 
     internal::do_function_values<dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_values,
       shape_function_data,
       values);
@@ -1581,7 +427,7 @@ namespace FEValuesViews
     fe_values->present_cell.get_interpolated_dof_values(fe_function,
                                                         dof_values);
     internal::do_function_derivatives<1, dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_gradients,
       shape_function_data,
       gradients);
@@ -1605,7 +451,7 @@ namespace FEValuesViews
     AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
 
     internal::do_function_derivatives<1, dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_gradients,
       shape_function_data,
       gradients);
@@ -1633,7 +479,7 @@ namespace FEValuesViews
     fe_values->present_cell.get_interpolated_dof_values(fe_function,
                                                         dof_values);
     internal::do_function_derivatives<2, dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_hessians,
       shape_function_data,
       hessians);
@@ -1657,7 +503,7 @@ namespace FEValuesViews
     AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
 
     internal::do_function_derivatives<2, dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_hessians,
       shape_function_data,
       hessians);
@@ -1685,7 +531,7 @@ namespace FEValuesViews
     fe_values->present_cell.get_interpolated_dof_values(fe_function,
                                                         dof_values);
     internal::do_function_laplacians<dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_hessians,
       shape_function_data,
       laplacians);
@@ -1709,7 +555,7 @@ namespace FEValuesViews
     AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
 
     internal::do_function_laplacians<dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_hessians,
       shape_function_data,
       laplacians);
@@ -1738,7 +584,7 @@ namespace FEValuesViews
     fe_values->present_cell.get_interpolated_dof_values(fe_function,
                                                         dof_values);
     internal::do_function_derivatives<3, dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_3rd_derivatives,
       shape_function_data,
       third_derivatives);
@@ -1763,7 +609,7 @@ namespace FEValuesViews
     AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
 
     internal::do_function_derivatives<3, dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_3rd_derivatives,
       shape_function_data,
       third_derivatives);
@@ -1789,7 +635,7 @@ namespace FEValuesViews
     fe_values->present_cell.get_interpolated_dof_values(fe_function,
                                                         dof_values);
     internal::do_function_values<dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_values,
       shape_function_data,
       values);
@@ -1813,7 +659,7 @@ namespace FEValuesViews
     AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
 
     internal::do_function_values<dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_values,
       shape_function_data,
       values);
@@ -1841,7 +687,7 @@ namespace FEValuesViews
     fe_values->present_cell.get_interpolated_dof_values(fe_function,
                                                         dof_values);
     internal::do_function_derivatives<1, dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_gradients,
       shape_function_data,
       gradients);
@@ -1865,7 +711,7 @@ namespace FEValuesViews
     AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
 
     internal::do_function_derivatives<1, dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_gradients,
       shape_function_data,
       gradients);
@@ -1894,7 +740,7 @@ namespace FEValuesViews
     fe_values->present_cell.get_interpolated_dof_values(fe_function,
                                                         dof_values);
     internal::do_function_symmetric_gradients<dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_gradients,
       shape_function_data,
       symmetric_gradients);
@@ -1919,7 +765,7 @@ namespace FEValuesViews
     AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
 
     internal::do_function_symmetric_gradients<dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_gradients,
       shape_function_data,
       symmetric_gradients);
@@ -1948,7 +794,7 @@ namespace FEValuesViews
     fe_values->present_cell.get_interpolated_dof_values(fe_function,
                                                         dof_values);
     internal::do_function_divergences<dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_gradients,
       shape_function_data,
       divergences);
@@ -1972,7 +818,7 @@ namespace FEValuesViews
     AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
 
     internal::do_function_divergences<dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_gradients,
       shape_function_data,
       divergences);
@@ -2000,7 +846,7 @@ namespace FEValuesViews
     fe_values->present_cell.get_interpolated_dof_values(fe_function,
                                                         dof_values);
     internal::do_function_curls<dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_gradients,
       shape_function_data,
       curls);
@@ -2024,7 +870,7 @@ namespace FEValuesViews
     AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
 
     internal::do_function_curls<dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_gradients,
       shape_function_data,
       curls);
@@ -2052,7 +898,7 @@ namespace FEValuesViews
     fe_values->present_cell.get_interpolated_dof_values(fe_function,
                                                         dof_values);
     internal::do_function_derivatives<2, dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_hessians,
       shape_function_data,
       hessians);
@@ -2076,7 +922,7 @@ namespace FEValuesViews
     AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
 
     internal::do_function_derivatives<2, dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_hessians,
       shape_function_data,
       hessians);
@@ -2109,7 +955,7 @@ namespace FEValuesViews
     fe_values->present_cell.get_interpolated_dof_values(fe_function,
                                                         dof_values);
     internal::do_function_laplacians<dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_hessians,
       shape_function_data,
       laplacians);
@@ -2136,7 +982,7 @@ namespace FEValuesViews
     AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
 
     internal::do_function_laplacians<dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_hessians,
       shape_function_data,
       laplacians);
@@ -2165,7 +1011,7 @@ namespace FEValuesViews
     fe_values->present_cell.get_interpolated_dof_values(fe_function,
                                                         dof_values);
     internal::do_function_derivatives<3, dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_3rd_derivatives,
       shape_function_data,
       third_derivatives);
@@ -2190,7 +1036,7 @@ namespace FEValuesViews
     AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
 
     internal::do_function_derivatives<3, dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_3rd_derivatives,
       shape_function_data,
       third_derivatives);
@@ -2216,7 +1062,7 @@ namespace FEValuesViews
     fe_values->present_cell.get_interpolated_dof_values(fe_function,
                                                         dof_values);
     internal::do_function_values<dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_values,
       shape_function_data,
       values);
@@ -2240,7 +1086,7 @@ namespace FEValuesViews
     AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
 
     internal::do_function_values<dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_values,
       shape_function_data,
       values);
@@ -2269,7 +1115,7 @@ namespace FEValuesViews
     fe_values->present_cell.get_interpolated_dof_values(fe_function,
                                                         dof_values);
     internal::do_function_divergences<dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_gradients,
       shape_function_data,
       divergences);
@@ -2294,7 +1140,7 @@ namespace FEValuesViews
     AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
 
     internal::do_function_divergences<dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_gradients,
       shape_function_data,
       divergences);
@@ -2320,7 +1166,7 @@ namespace FEValuesViews
     fe_values->present_cell.get_interpolated_dof_values(fe_function,
                                                         dof_values);
     internal::do_function_values<dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_values,
       shape_function_data,
       values);
@@ -2344,7 +1190,7 @@ namespace FEValuesViews
     AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
 
     internal::do_function_values<dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_values,
       shape_function_data,
       values);
@@ -2373,7 +1219,7 @@ namespace FEValuesViews
     fe_values->present_cell.get_interpolated_dof_values(fe_function,
                                                         dof_values);
     internal::do_function_divergences<dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_gradients,
       shape_function_data,
       divergences);
@@ -2397,7 +1243,7 @@ namespace FEValuesViews
     AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
 
     internal::do_function_divergences<dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_gradients,
       shape_function_data,
       divergences);
@@ -2426,7 +1272,7 @@ namespace FEValuesViews
     fe_values->present_cell.get_interpolated_dof_values(fe_function,
                                                         dof_values);
     internal::do_function_gradients<dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_gradients,
       shape_function_data,
       gradients);
@@ -2450,7 +1296,7 @@ namespace FEValuesViews
     AssertDimension(dof_values.size(), fe_values->dofs_per_cell);
 
     internal::do_function_gradients<dim, spacedim>(
-      make_array_view(dof_values.begin(), dof_values.end()),
+      make_const_array_view(dof_values),
       fe_values->finite_element_output.shape_gradients,
       shape_function_data,
       gradients);
diff --git a/source/fe/fe_values_views_internal.cc b/source/fe/fe_values_views_internal.cc
new file mode 100644 (file)
index 0000000..4a038e2
--- /dev/null
@@ -0,0 +1,1167 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 2023 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/array_view.h>
+#include <deal.II/base/numbers.h>
+
+#include <deal.II/differentiation/ad.h>
+
+#include <deal.II/fe/fe_values_views_internal.h>
+
+#include <type_traits>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace FEValuesViews
+{
+  namespace internal
+  {
+    namespace
+    {
+      // Check to see if a DoF value is zero, implying that subsequent
+      // operations with the value have no effect.
+      template <typename Number, typename T = void>
+      struct CheckForZero
+      {
+        static bool
+        value(const Number &value)
+        {
+          return value == dealii::internal::NumberType<Number>::value(0.0);
+        }
+      };
+
+      // For auto-differentiable numbers, the fact that a DoF value is zero
+      // does not imply that its derivatives are zero as well. So we
+      // can't filter by value for these number types.
+      // Note that we also want to avoid actually checking the value itself,
+      // since some AD numbers are not contextually convertible to booleans.
+      template <typename Number>
+      struct CheckForZero<
+        Number,
+        std::enable_if_t<Differentiation::AD::is_ad_number<Number>::value>>
+      {
+        static bool
+        value(const Number & /*value*/)
+        {
+          return false;
+        }
+      };
+    } // namespace
+
+    template <int dim, int spacedim, typename Number>
+    void
+    do_function_values(
+      const ArrayView<const Number> &dof_values,
+      const Table<2, double> &       shape_values,
+      const std::vector<typename Scalar<dim, spacedim>::ShapeFunctionData>
+        &shape_function_data,
+      std::vector<typename ProductType<Number, double>::type> &values)
+    {
+      const unsigned int dofs_per_cell       = dof_values.size();
+      const unsigned int n_quadrature_points = values.size();
+
+      std::fill(values.begin(),
+                values.end(),
+                dealii::internal::NumberType<Number>::value(0.0));
+
+      for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
+           ++shape_function)
+        if (shape_function_data[shape_function]
+              .is_nonzero_shape_function_component)
+          {
+            const Number &value = dof_values[shape_function];
+            // For auto-differentiable numbers, the fact that a DoF value is
+            // zero does not imply that its derivatives are zero as well. So we
+            // can't filter by value for these number types.
+            if (CheckForZero<Number>::value(value) == true)
+              continue;
+
+            const double *shape_value_ptr =
+              &shape_values(shape_function_data[shape_function].row_index, 0);
+            for (unsigned int q_point = 0; q_point < n_quadrature_points;
+                 ++q_point)
+              values[q_point] += value * (*shape_value_ptr++);
+          }
+    }
+
+
+
+    template <int order, int dim, int spacedim, typename Number>
+    void
+    do_function_derivatives(
+      const ArrayView<const Number> &                  dof_values,
+      const Table<2, dealii::Tensor<order, spacedim>> &shape_derivatives,
+      const std::vector<typename Scalar<dim, spacedim>::ShapeFunctionData>
+        &shape_function_data,
+      std::vector<
+        typename ProductType<Number, dealii::Tensor<order, spacedim>>::type>
+        &derivatives)
+    {
+      const unsigned int dofs_per_cell       = dof_values.size();
+      const unsigned int n_quadrature_points = derivatives.size();
+
+      std::fill(
+        derivatives.begin(),
+        derivatives.end(),
+        typename ProductType<Number, dealii::Tensor<order, spacedim>>::type());
+
+      for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
+           ++shape_function)
+        if (shape_function_data[shape_function]
+              .is_nonzero_shape_function_component)
+          {
+            const Number &value = dof_values[shape_function];
+            // For auto-differentiable numbers, the fact that a DoF value is
+            // zero does not imply that its derivatives are zero as well. So we
+            // can't filter by value for these number types.
+            if (CheckForZero<Number>::value(value) == true)
+              continue;
+
+            const dealii::Tensor<order, spacedim> *shape_derivative_ptr =
+              &shape_derivatives[shape_function_data[shape_function].row_index]
+                                [0];
+            for (unsigned int q_point = 0; q_point < n_quadrature_points;
+                 ++q_point)
+              derivatives[q_point] += value * (*shape_derivative_ptr++);
+          }
+    }
+
+
+
+    template <int dim, int spacedim, typename Number>
+    void
+    do_function_laplacians(
+      const ArrayView<const Number> &              dof_values,
+      const Table<2, dealii::Tensor<2, spacedim>> &shape_hessians,
+      const std::vector<typename Scalar<dim, spacedim>::ShapeFunctionData>
+        &shape_function_data,
+      std::vector<typename Scalar<dim, spacedim>::
+                    template solution_laplacian_type<Number>> &laplacians)
+    {
+      const unsigned int dofs_per_cell       = dof_values.size();
+      const unsigned int n_quadrature_points = laplacians.size();
+
+      std::fill(
+        laplacians.begin(),
+        laplacians.end(),
+        typename Scalar<dim,
+                        spacedim>::template solution_laplacian_type<Number>());
+
+      for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
+           ++shape_function)
+        if (shape_function_data[shape_function]
+              .is_nonzero_shape_function_component)
+          {
+            const Number &value = dof_values[shape_function];
+            // For auto-differentiable numbers, the fact that a DoF value is
+            // zero does not imply that its derivatives are zero as well. So we
+            // can't filter by value for these number types.
+            if (CheckForZero<Number>::value(value) == true)
+              continue;
+
+            const dealii::Tensor<2, spacedim> *shape_hessian_ptr =
+              &shape_hessians[shape_function_data[shape_function].row_index][0];
+            for (unsigned int q_point = 0; q_point < n_quadrature_points;
+                 ++q_point)
+              laplacians[q_point] += value * trace(*shape_hessian_ptr++);
+          }
+    }
+
+
+
+    // ----------------------------- vector part ---------------------------
+
+    template <int dim, int spacedim, typename Number>
+    void
+    do_function_values(
+      const ArrayView<const Number> &dof_values,
+      const Table<2, double> &       shape_values,
+      const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
+        &shape_function_data,
+      std::vector<
+        typename ProductType<Number, dealii::Tensor<1, spacedim>>::type>
+        &values)
+    {
+      const unsigned int dofs_per_cell       = dof_values.size();
+      const unsigned int n_quadrature_points = values.size();
+
+      std::fill(
+        values.begin(),
+        values.end(),
+        typename ProductType<Number, dealii::Tensor<1, spacedim>>::type());
+
+      for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
+           ++shape_function)
+        {
+          const int snc =
+            shape_function_data[shape_function].single_nonzero_component;
+
+          if (snc == -2)
+            // shape function is zero for the selected components
+            continue;
+
+          const Number &value = dof_values[shape_function];
+          // For auto-differentiable numbers, the fact that a DoF value is zero
+          // does not imply that its derivatives are zero as well. So we
+          // can't filter by value for these number types.
+          if (CheckForZero<Number>::value(value) == true)
+            continue;
+
+          if (snc != -1)
+            {
+              const unsigned int comp = shape_function_data[shape_function]
+                                          .single_nonzero_component_index;
+              const double *shape_value_ptr = &shape_values(snc, 0);
+              for (unsigned int q_point = 0; q_point < n_quadrature_points;
+                   ++q_point)
+                values[q_point][comp] += value * (*shape_value_ptr++);
+            }
+          else
+            for (unsigned int d = 0; d < spacedim; ++d)
+              if (shape_function_data[shape_function]
+                    .is_nonzero_shape_function_component[d])
+                {
+                  const double *shape_value_ptr = &shape_values(
+                    shape_function_data[shape_function].row_index[d], 0);
+                  for (unsigned int q_point = 0; q_point < n_quadrature_points;
+                       ++q_point)
+                    values[q_point][d] += value * (*shape_value_ptr++);
+                }
+        }
+    }
+
+
+
+    template <int order, int dim, int spacedim, typename Number>
+    void
+    do_function_derivatives(
+      const ArrayView<const Number> &                  dof_values,
+      const Table<2, dealii::Tensor<order, spacedim>> &shape_derivatives,
+      const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
+        &shape_function_data,
+      std::vector<
+        typename ProductType<Number, dealii::Tensor<order + 1, spacedim>>::type>
+        &derivatives)
+    {
+      const unsigned int dofs_per_cell       = dof_values.size();
+      const unsigned int n_quadrature_points = derivatives.size();
+
+      std::fill(
+        derivatives.begin(),
+        derivatives.end(),
+        typename ProductType<Number,
+                             dealii::Tensor<order + 1, spacedim>>::type());
+
+      for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
+           ++shape_function)
+        {
+          const int snc =
+            shape_function_data[shape_function].single_nonzero_component;
+
+          if (snc == -2)
+            // shape function is zero for the selected components
+            continue;
+
+          const Number &value = dof_values[shape_function];
+          // For auto-differentiable numbers, the fact that a DoF value is zero
+          // does not imply that its derivatives are zero as well. So we
+          // can't filter by value for these number types.
+          if (CheckForZero<Number>::value(value) == true)
+            continue;
+
+          if (snc != -1)
+            {
+              const unsigned int comp = shape_function_data[shape_function]
+                                          .single_nonzero_component_index;
+              const dealii::Tensor<order, spacedim> *shape_derivative_ptr =
+                &shape_derivatives[snc][0];
+              for (unsigned int q_point = 0; q_point < n_quadrature_points;
+                   ++q_point)
+                derivatives[q_point][comp] += value * (*shape_derivative_ptr++);
+            }
+          else
+            for (unsigned int d = 0; d < spacedim; ++d)
+              if (shape_function_data[shape_function]
+                    .is_nonzero_shape_function_component[d])
+                {
+                  const dealii::Tensor<order, spacedim> *shape_derivative_ptr =
+                    &shape_derivatives[shape_function_data[shape_function]
+                                         .row_index[d]][0];
+                  for (unsigned int q_point = 0; q_point < n_quadrature_points;
+                       ++q_point)
+                    derivatives[q_point][d] +=
+                      value * (*shape_derivative_ptr++);
+                }
+        }
+    }
+
+
+
+    template <int dim, int spacedim, typename Number>
+    void
+    do_function_symmetric_gradients(
+      const ArrayView<const Number> &              dof_values,
+      const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients,
+      const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
+        &shape_function_data,
+      std::vector<
+        typename ProductType<Number,
+                             dealii::SymmetricTensor<2, spacedim>>::type>
+        &symmetric_gradients)
+    {
+      const unsigned int dofs_per_cell       = dof_values.size();
+      const unsigned int n_quadrature_points = symmetric_gradients.size();
+
+      std::fill(
+        symmetric_gradients.begin(),
+        symmetric_gradients.end(),
+        typename ProductType<Number,
+                             dealii::SymmetricTensor<2, spacedim>>::type());
+
+      for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
+           ++shape_function)
+        {
+          const int snc =
+            shape_function_data[shape_function].single_nonzero_component;
+
+          if (snc == -2)
+            // shape function is zero for the selected components
+            continue;
+
+          const Number &value = dof_values[shape_function];
+          // For auto-differentiable numbers, the fact that a DoF value is zero
+          // does not imply that its derivatives are zero as well. So we
+          // can't filter by value for these number types.
+          if (CheckForZero<Number>::value(value) == true)
+            continue;
+
+          if (snc != -1)
+            {
+              const unsigned int comp = shape_function_data[shape_function]
+                                          .single_nonzero_component_index;
+              const dealii::Tensor<1, spacedim> *shape_gradient_ptr =
+                &shape_gradients[snc][0];
+              for (unsigned int q_point = 0; q_point < n_quadrature_points;
+                   ++q_point)
+                symmetric_gradients[q_point] +=
+                  value * dealii::SymmetricTensor<2, spacedim>(
+                            symmetrize_single_row(comp, *shape_gradient_ptr++));
+            }
+          else
+            for (unsigned int q_point = 0; q_point < n_quadrature_points;
+                 ++q_point)
+              {
+                typename ProductType<Number, dealii::Tensor<2, spacedim>>::type
+                  grad;
+                for (unsigned int d = 0; d < spacedim; ++d)
+                  if (shape_function_data[shape_function]
+                        .is_nonzero_shape_function_component[d])
+                    grad[d] =
+                      value *
+                      shape_gradients[shape_function_data[shape_function]
+                                        .row_index[d]][q_point];
+                symmetric_gradients[q_point] += symmetrize(grad);
+              }
+        }
+    }
+
+
+
+    template <int dim, int spacedim, typename Number>
+    void
+    do_function_divergences(
+      const ArrayView<const Number> &              dof_values,
+      const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients,
+      const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
+        &shape_function_data,
+      std::vector<typename Vector<dim, spacedim>::
+                    template solution_divergence_type<Number>> &divergences)
+    {
+      const unsigned int dofs_per_cell       = dof_values.size();
+      const unsigned int n_quadrature_points = divergences.size();
+
+      std::fill(
+        divergences.begin(),
+        divergences.end(),
+        typename Vector<dim,
+                        spacedim>::template solution_divergence_type<Number>());
+
+      for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
+           ++shape_function)
+        {
+          const int snc =
+            shape_function_data[shape_function].single_nonzero_component;
+
+          if (snc == -2)
+            // shape function is zero for the selected components
+            continue;
+
+          const Number &value = dof_values[shape_function];
+          // For auto-differentiable numbers, the fact that a DoF value is zero
+          // does not imply that its derivatives are zero as well. So we
+          // can't filter by value for these number types.
+          if (CheckForZero<Number>::value(value) == true)
+            continue;
+
+          if (snc != -1)
+            {
+              const unsigned int comp = shape_function_data[shape_function]
+                                          .single_nonzero_component_index;
+              const dealii::Tensor<1, spacedim> *shape_gradient_ptr =
+                &shape_gradients[snc][0];
+              for (unsigned int q_point = 0; q_point < n_quadrature_points;
+                   ++q_point)
+                divergences[q_point] += value * (*shape_gradient_ptr++)[comp];
+            }
+          else
+            for (unsigned int d = 0; d < spacedim; ++d)
+              if (shape_function_data[shape_function]
+                    .is_nonzero_shape_function_component[d])
+                {
+                  const dealii::Tensor<1, spacedim> *shape_gradient_ptr =
+                    &shape_gradients[shape_function_data[shape_function]
+                                       .row_index[d]][0];
+                  for (unsigned int q_point = 0; q_point < n_quadrature_points;
+                       ++q_point)
+                    divergences[q_point] += value * (*shape_gradient_ptr++)[d];
+                }
+        }
+    }
+
+
+
+    template <int dim, int spacedim, typename Number>
+    void
+    do_function_curls(
+      const ArrayView<const Number> &              dof_values,
+      const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients,
+      const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
+        &shape_function_data,
+      std::vector<typename ProductType<
+        Number,
+        typename dealii::internal::CurlType<spacedim>::type>::type> &curls)
+    {
+      const unsigned int dofs_per_cell       = dof_values.size();
+      const unsigned int n_quadrature_points = curls.size();
+
+      std::fill(curls.begin(),
+                curls.end(),
+                typename ProductType<
+                  Number,
+                  typename dealii::internal::CurlType<spacedim>::type>::type());
+
+      switch (spacedim)
+        {
+          case 1:
+            {
+              Assert(false,
+                     ExcMessage(
+                       "Computing the curl in 1d is not a useful operation"));
+              break;
+            }
+
+          case 2:
+            {
+              for (unsigned int shape_function = 0;
+                   shape_function < dofs_per_cell;
+                   ++shape_function)
+                {
+                  const int snc = shape_function_data[shape_function]
+                                    .single_nonzero_component;
+
+                  if (snc == -2)
+                    // shape function is zero for the selected components
+                    continue;
+
+                  const Number &value = dof_values[shape_function];
+                  // For auto-differentiable numbers, the fact that a DoF value
+                  // is zero does not imply that its derivatives are zero as
+                  // well. So we can't filter by value for these number types.
+                  if (CheckForZero<Number>::value(value) == true)
+                    continue;
+
+                  if (snc != -1)
+                    {
+                      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];
+                    }
+                  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])
+                        {
+                          const dealii::Tensor<1,
+                                               spacedim> *shape_gradient_ptr =
+                            &shape_gradients[shape_function_data[shape_function]
+                                               .row_index[0]][0];
+
+                          for (unsigned int q_point = 0;
+                               q_point < n_quadrature_points;
+                               ++q_point)
+                            curls[q_point][0] -=
+                              value * (*shape_gradient_ptr++)[1];
+                        }
+
+                      if (shape_function_data[shape_function]
+                            .is_nonzero_shape_function_component[1])
+                        {
+                          const dealii::Tensor<1,
+                                               spacedim> *shape_gradient_ptr =
+                            &shape_gradients[shape_function_data[shape_function]
+                                               .row_index[1]][0];
+
+                          for (unsigned int q_point = 0;
+                               q_point < n_quadrature_points;
+                               ++q_point)
+                            curls[q_point][0] +=
+                              value * (*shape_gradient_ptr++)[0];
+                        }
+                    }
+                }
+              break;
+            }
+
+          case 3:
+            {
+              for (unsigned int shape_function = 0;
+                   shape_function < dofs_per_cell;
+                   ++shape_function)
+                {
+                  const int snc = shape_function_data[shape_function]
+                                    .single_nonzero_component;
+
+                  if (snc == -2)
+                    // shape function is zero for the selected components
+                    continue;
+
+                  const Number &value = dof_values[shape_function];
+                  // For auto-differentiable numbers, the fact that a DoF value
+                  // is zero does not imply that its derivatives are zero as
+                  // well. So we can't filter by value for these number types.
+                  if (CheckForZero<Number>::value(value) == true)
+                    continue;
+
+                  if (snc != -1)
+                    {
+                      const dealii::Tensor<1, spacedim> *shape_gradient_ptr =
+                        &shape_gradients[snc][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][1] +=
+                                    value * (*shape_gradient_ptr)[2];
+                                  curls[q_point][2] -=
+                                    value * (*shape_gradient_ptr++)[1];
+                                }
+
+                              break;
+                            }
+
+                          case 1:
+                            {
+                              for (unsigned int q_point = 0;
+                                   q_point < n_quadrature_points;
+                                   ++q_point)
+                                {
+                                  curls[q_point][0] -=
+                                    value * (*shape_gradient_ptr)[2];
+                                  curls[q_point][2] +=
+                                    value * (*shape_gradient_ptr++)[0];
+                                }
+
+                              break;
+                            }
+
+                          case 2:
+                            {
+                              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])
+                        {
+                          const dealii::Tensor<1,
+                                               spacedim> *shape_gradient_ptr =
+                            &shape_gradients[shape_function_data[shape_function]
+                                               .row_index[0]][0];
+
+                          for (unsigned int q_point = 0;
+                               q_point < n_quadrature_points;
+                               ++q_point)
+                            {
+                              curls[q_point][1] +=
+                                value * (*shape_gradient_ptr)[2];
+                              curls[q_point][2] -=
+                                value * (*shape_gradient_ptr++)[1];
+                            }
+                        }
+
+                      if (shape_function_data[shape_function]
+                            .is_nonzero_shape_function_component[1])
+                        {
+                          const dealii::Tensor<1,
+                                               spacedim> *shape_gradient_ptr =
+                            &shape_gradients[shape_function_data[shape_function]
+                                               .row_index[1]][0];
+
+                          for (unsigned int q_point = 0;
+                               q_point < n_quadrature_points;
+                               ++q_point)
+                            {
+                              curls[q_point][0] -=
+                                value * (*shape_gradient_ptr)[2];
+                              curls[q_point][2] +=
+                                value * (*shape_gradient_ptr++)[0];
+                            }
+                        }
+
+                      if (shape_function_data[shape_function]
+                            .is_nonzero_shape_function_component[2])
+                        {
+                          const dealii::Tensor<1,
+                                               spacedim> *shape_gradient_ptr =
+                            &shape_gradients[shape_function_data[shape_function]
+                                               .row_index[2]][0];
+
+                          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];
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+
+
+    template <int dim, int spacedim, typename Number>
+    void
+    do_function_laplacians(
+      const ArrayView<const Number> &              dof_values,
+      const Table<2, dealii::Tensor<2, spacedim>> &shape_hessians,
+      const std::vector<typename Vector<dim, spacedim>::ShapeFunctionData>
+        &shape_function_data,
+      std::vector<typename Vector<dim, spacedim>::
+                    template solution_laplacian_type<Number>> &laplacians)
+    {
+      const unsigned int dofs_per_cell       = dof_values.size();
+      const unsigned int n_quadrature_points = laplacians.size();
+
+      std::fill(
+        laplacians.begin(),
+        laplacians.end(),
+        typename Vector<dim,
+                        spacedim>::template solution_laplacian_type<Number>());
+
+      for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
+           ++shape_function)
+        {
+          const int snc =
+            shape_function_data[shape_function].single_nonzero_component;
+
+          if (snc == -2)
+            // shape function is zero for the selected components
+            continue;
+
+          const Number &value = dof_values[shape_function];
+          // For auto-differentiable numbers, the fact that a DoF value is zero
+          // does not imply that its derivatives are zero as well. So we
+          // can't filter by value for these number types.
+          if (CheckForZero<Number>::value(value) == true)
+            continue;
+
+          if (snc != -1)
+            {
+              const unsigned int comp = shape_function_data[shape_function]
+                                          .single_nonzero_component_index;
+              const dealii::Tensor<2, spacedim> *shape_hessian_ptr =
+                &shape_hessians[snc][0];
+              for (unsigned int q_point = 0; q_point < n_quadrature_points;
+                   ++q_point)
+                laplacians[q_point][comp] +=
+                  value * trace(*shape_hessian_ptr++);
+            }
+          else
+            for (unsigned int d = 0; d < spacedim; ++d)
+              if (shape_function_data[shape_function]
+                    .is_nonzero_shape_function_component[d])
+                {
+                  const dealii::Tensor<2, spacedim> *shape_hessian_ptr =
+                    &shape_hessians[shape_function_data[shape_function]
+                                      .row_index[d]][0];
+                  for (unsigned int q_point = 0; q_point < n_quadrature_points;
+                       ++q_point)
+                    laplacians[q_point][d] +=
+                      value * trace(*shape_hessian_ptr++);
+                }
+        }
+    }
+
+
+
+    // ---------------------- symmetric tensor part ------------------------
+
+    template <int dim, int spacedim, typename Number>
+    void
+    do_function_values(
+      const ArrayView<const Number> & dof_values,
+      const dealii::Table<2, double> &shape_values,
+      const std::vector<
+        typename SymmetricTensor<2, dim, spacedim>::ShapeFunctionData>
+        &shape_function_data,
+      std::vector<
+        typename ProductType<Number,
+                             dealii::SymmetricTensor<2, spacedim>>::type>
+        &values)
+    {
+      const unsigned int dofs_per_cell       = dof_values.size();
+      const unsigned int n_quadrature_points = values.size();
+
+      std::fill(
+        values.begin(),
+        values.end(),
+        typename ProductType<Number,
+                             dealii::SymmetricTensor<2, spacedim>>::type());
+
+      for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
+           ++shape_function)
+        {
+          const int snc =
+            shape_function_data[shape_function].single_nonzero_component;
+
+          if (snc == -2)
+            // shape function is zero for the selected components
+            continue;
+
+          const Number &value = dof_values[shape_function];
+          // For auto-differentiable numbers, the fact that a DoF value is zero
+          // does not imply that its derivatives are zero as well. So we
+          // can't filter by value for these number types.
+          if (CheckForZero<Number>::value(value) == true)
+            continue;
+
+          if (snc != -1)
+            {
+              const TableIndices<2> comp = dealii::
+                SymmetricTensor<2, spacedim>::unrolled_to_component_indices(
+                  shape_function_data[shape_function]
+                    .single_nonzero_component_index);
+              const double *shape_value_ptr = &shape_values(snc, 0);
+              for (unsigned int q_point = 0; q_point < n_quadrature_points;
+                   ++q_point)
+                values[q_point][comp] += value * (*shape_value_ptr++);
+            }
+          else
+            for (unsigned int d = 0;
+                 d <
+                 dealii::SymmetricTensor<2, spacedim>::n_independent_components;
+                 ++d)
+              if (shape_function_data[shape_function]
+                    .is_nonzero_shape_function_component[d])
+                {
+                  const TableIndices<2> comp =
+                    dealii::SymmetricTensor<2, spacedim>::
+                      unrolled_to_component_indices(d);
+                  const double *shape_value_ptr = &shape_values(
+                    shape_function_data[shape_function].row_index[d], 0);
+                  for (unsigned int q_point = 0; q_point < n_quadrature_points;
+                       ++q_point)
+                    values[q_point][comp] += value * (*shape_value_ptr++);
+                }
+        }
+    }
+
+
+
+    template <int dim, int spacedim, typename Number>
+    void
+    do_function_divergences(
+      const ArrayView<const Number> &              dof_values,
+      const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients,
+      const std::vector<
+        typename SymmetricTensor<2, dim, spacedim>::ShapeFunctionData>
+        &shape_function_data,
+      std::vector<typename SymmetricTensor<2, dim, spacedim>::
+                    template solution_divergence_type<Number>> &divergences)
+    {
+      const unsigned int dofs_per_cell       = dof_values.size();
+      const unsigned int n_quadrature_points = divergences.size();
+
+      std::fill(divergences.begin(),
+                divergences.end(),
+                typename SymmetricTensor<2, dim, spacedim>::
+                  template solution_divergence_type<Number>());
+
+      for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
+           ++shape_function)
+        {
+          const int snc =
+            shape_function_data[shape_function].single_nonzero_component;
+
+          if (snc == -2)
+            // shape function is zero for the selected components
+            continue;
+
+          const Number &value = dof_values[shape_function];
+          // For auto-differentiable numbers, the fact that a DoF value is zero
+          // does not imply that its derivatives are zero as well. So we
+          // can't filter by value for these number types.
+          if (CheckForZero<Number>::value(value) == true)
+            continue;
+
+          if (snc != -1)
+            {
+              const unsigned int comp = shape_function_data[shape_function]
+                                          .single_nonzero_component_index;
+
+              const dealii::Tensor<1, spacedim> *shape_gradient_ptr =
+                &shape_gradients[snc][0];
+
+              const unsigned int ii = dealii::SymmetricTensor<2, spacedim>::
+                unrolled_to_component_indices(comp)[0];
+              const unsigned int jj = dealii::SymmetricTensor<2, spacedim>::
+                unrolled_to_component_indices(comp)[1];
+
+              for (unsigned int q_point = 0; q_point < n_quadrature_points;
+                   ++q_point, ++shape_gradient_ptr)
+                {
+                  divergences[q_point][ii] += value * (*shape_gradient_ptr)[jj];
+
+                  if (ii != jj)
+                    divergences[q_point][jj] +=
+                      value * (*shape_gradient_ptr)[ii];
+                }
+            }
+          else
+            {
+              for (unsigned int d = 0;
+                   d <
+                   dealii::SymmetricTensor<2,
+                                           spacedim>::n_independent_components;
+                   ++d)
+                if (shape_function_data[shape_function]
+                      .is_nonzero_shape_function_component[d])
+                  {
+                    Assert(false, ExcNotImplemented());
+
+                    // the following implementation needs to be looked over -- I
+                    // think it can't be right, because we are in a case where
+                    // there is no single nonzero component
+                    //
+                    // the following is not implemented! we need to consider the
+                    // interplay between multiple non-zero entries in shape
+                    // function and the representation as a symmetric
+                    // second-order tensor
+                    const unsigned int comp =
+                      shape_function_data[shape_function]
+                        .single_nonzero_component_index;
+
+                    const dealii::Tensor<1, spacedim> *shape_gradient_ptr =
+                      &shape_gradients[shape_function_data[shape_function]
+                                         .row_index[d]][0];
+                    for (unsigned int q_point = 0;
+                         q_point < n_quadrature_points;
+                         ++q_point, ++shape_gradient_ptr)
+                      {
+                        for (unsigned int j = 0; j < spacedim; ++j)
+                          {
+                            const unsigned int vector_component =
+                              dealii::SymmetricTensor<2, spacedim>::
+                                component_to_unrolled_index(
+                                  TableIndices<2>(comp, j));
+                            divergences[q_point][vector_component] +=
+                              value * (*shape_gradient_ptr++)[j];
+                          }
+                      }
+                  }
+            }
+        }
+    }
+
+    // ---------------------- non-symmetric tensor part ------------------------
+
+    template <int dim, int spacedim, typename Number>
+    void
+    do_function_values(
+      const ArrayView<const Number> & dof_values,
+      const dealii::Table<2, double> &shape_values,
+      const std::vector<typename Tensor<2, dim, spacedim>::ShapeFunctionData>
+        &shape_function_data,
+      std::vector<
+        typename ProductType<Number, dealii::Tensor<2, spacedim>>::type>
+        &values)
+    {
+      const unsigned int dofs_per_cell       = dof_values.size();
+      const unsigned int n_quadrature_points = values.size();
+
+      std::fill(
+        values.begin(),
+        values.end(),
+        typename ProductType<Number, dealii::Tensor<2, spacedim>>::type());
+
+      for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
+           ++shape_function)
+        {
+          const int snc =
+            shape_function_data[shape_function].single_nonzero_component;
+
+          if (snc == -2)
+            // shape function is zero for the selected components
+            continue;
+
+          const Number &value = dof_values[shape_function];
+          // For auto-differentiable numbers, the fact that a DoF value is zero
+          // does not imply that its derivatives are zero as well. So we
+          // can't filter by value for these number types.
+          if (CheckForZero<Number>::value(value) == true)
+            continue;
+
+          if (snc != -1)
+            {
+              const unsigned int comp = shape_function_data[shape_function]
+                                          .single_nonzero_component_index;
+
+              const TableIndices<2> indices =
+                dealii::Tensor<2, spacedim>::unrolled_to_component_indices(
+                  comp);
+
+              const double *shape_value_ptr = &shape_values(snc, 0);
+              for (unsigned int q_point = 0; q_point < n_quadrature_points;
+                   ++q_point)
+                values[q_point][indices] += value * (*shape_value_ptr++);
+            }
+          else
+            for (unsigned int d = 0; d < dim * dim; ++d)
+              if (shape_function_data[shape_function]
+                    .is_nonzero_shape_function_component[d])
+                {
+                  const TableIndices<2> indices =
+                    dealii::Tensor<2, spacedim>::unrolled_to_component_indices(
+                      d);
+
+                  const double *shape_value_ptr = &shape_values(
+                    shape_function_data[shape_function].row_index[d], 0);
+                  for (unsigned int q_point = 0; q_point < n_quadrature_points;
+                       ++q_point)
+                    values[q_point][indices] += value * (*shape_value_ptr++);
+                }
+        }
+    }
+
+
+
+    template <int dim, int spacedim, typename Number>
+    void
+    do_function_divergences(
+      const ArrayView<const Number> &              dof_values,
+      const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients,
+      const std::vector<typename Tensor<2, dim, spacedim>::ShapeFunctionData>
+        &shape_function_data,
+      std::vector<typename Tensor<2, dim, spacedim>::
+                    template solution_divergence_type<Number>> &divergences)
+    {
+      const unsigned int dofs_per_cell       = dof_values.size();
+      const unsigned int n_quadrature_points = divergences.size();
+
+      std::fill(
+        divergences.begin(),
+        divergences.end(),
+        typename Tensor<2, dim, spacedim>::template solution_divergence_type<
+          Number>());
+
+      for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
+           ++shape_function)
+        {
+          const int snc =
+            shape_function_data[shape_function].single_nonzero_component;
+
+          if (snc == -2)
+            // shape function is zero for the selected components
+            continue;
+
+          const Number &value = dof_values[shape_function];
+          // For auto-differentiable numbers, the fact that a DoF value is zero
+          // does not imply that its derivatives are zero as well. So we
+          // can't filter by value for these number types.
+          if (CheckForZero<Number>::value(value) == true)
+            continue;
+
+          if (snc != -1)
+            {
+              const unsigned int comp = shape_function_data[shape_function]
+                                          .single_nonzero_component_index;
+
+              const dealii::Tensor<1, spacedim> *shape_gradient_ptr =
+                &shape_gradients[snc][0];
+
+              const TableIndices<2> indices =
+                dealii::Tensor<2, spacedim>::unrolled_to_component_indices(
+                  comp);
+              const unsigned int ii = indices[0];
+              const unsigned int jj = indices[1];
+
+              for (unsigned int q_point = 0; q_point < n_quadrature_points;
+                   ++q_point, ++shape_gradient_ptr)
+                {
+                  divergences[q_point][ii] += value * (*shape_gradient_ptr)[jj];
+                }
+            }
+          else
+            {
+              for (unsigned int d = 0; d < dim * dim; ++d)
+                if (shape_function_data[shape_function]
+                      .is_nonzero_shape_function_component[d])
+                  {
+                    Assert(false, ExcNotImplemented());
+                  }
+            }
+        }
+    }
+
+
+
+    template <int dim, int spacedim, typename Number>
+    void
+    do_function_gradients(
+      const ArrayView<const Number> &              dof_values,
+      const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients,
+      const std::vector<typename Tensor<2, dim, spacedim>::ShapeFunctionData>
+        &shape_function_data,
+      std::vector<typename Tensor<2, dim, spacedim>::
+                    template solution_gradient_type<Number>> &gradients)
+    {
+      const unsigned int dofs_per_cell       = dof_values.size();
+      const unsigned int n_quadrature_points = gradients.size();
+
+      std::fill(
+        gradients.begin(),
+        gradients.end(),
+        typename Tensor<2, dim, spacedim>::template solution_gradient_type<
+          Number>());
+
+      for (unsigned int shape_function = 0; shape_function < dofs_per_cell;
+           ++shape_function)
+        {
+          const int snc =
+            shape_function_data[shape_function].single_nonzero_component;
+
+          if (snc == -2)
+            // shape function is zero for the selected components
+            continue;
+
+          const Number &value = dof_values[shape_function];
+          // For auto-differentiable numbers, the fact that a DoF value is zero
+          // does not imply that its derivatives are zero as well. So we
+          // can't filter by value for these number types.
+          if (CheckForZero<Number>::value(value) == true)
+            continue;
+
+          if (snc != -1)
+            {
+              const unsigned int comp = shape_function_data[shape_function]
+                                          .single_nonzero_component_index;
+
+              const dealii::Tensor<1, spacedim> *shape_gradient_ptr =
+                &shape_gradients[snc][0];
+
+              const TableIndices<2> indices =
+                dealii::Tensor<2, spacedim>::unrolled_to_component_indices(
+                  comp);
+              const unsigned int ii = indices[0];
+              const unsigned int jj = indices[1];
+
+              for (unsigned int q_point = 0; q_point < n_quadrature_points;
+                   ++q_point, ++shape_gradient_ptr)
+                {
+                  gradients[q_point][ii][jj] += value * (*shape_gradient_ptr);
+                }
+            }
+          else
+            {
+              for (unsigned int d = 0; d < dim * dim; ++d)
+                if (shape_function_data[shape_function]
+                      .is_nonzero_shape_function_component[d])
+                  {
+                    Assert(false, ExcNotImplemented());
+                  }
+            }
+        }
+    }
+  } // end of namespace internal
+} // namespace FEValuesViews
+
+
+
+/*------------------------------- Explicit Instantiations -------------*/
+
+#include "fe_values_views_internal.inst"
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/source/fe/fe_values_views_internal.inst.in b/source/fe/fe_values_views_internal.inst.in
new file mode 100644 (file)
index 0000000..2f4cb3a
--- /dev/null
@@ -0,0 +1,204 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 2022 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#ifndef DOXYGEN
+
+for (S : ALL_SCALAR_TYPES; deal_II_dimension : DIMENSIONS;
+     deal_II_space_dimension : SPACE_DIMENSIONS)
+  {
+#  if deal_II_dimension <= deal_II_space_dimension
+    namespace FEValuesViews
+
+    {
+      namespace internal
+
+      {
+        template void
+        do_function_values<deal_II_dimension, deal_II_space_dimension, S>(
+          const ArrayView<const S> &,
+          const Table<2, double> &,
+          const std::vector<Scalar<deal_II_dimension,
+                                   deal_II_space_dimension>::ShapeFunctionData>
+            &,
+          std::vector<ProductType<S, double>::type> &);
+
+
+        template void
+        do_function_laplacians<deal_II_dimension, deal_II_space_dimension, S>(
+          const ArrayView<const S> &,
+          const Table<2, dealii::Tensor<2, deal_II_space_dimension>> &,
+          const std::vector<Scalar<deal_II_dimension,
+                                   deal_II_space_dimension>::ShapeFunctionData>
+            &,
+          std::vector<Scalar<deal_II_dimension, deal_II_space_dimension>::
+                        template solution_laplacian_type<S>> &);
+
+        template void
+        do_function_values<deal_II_dimension, deal_II_space_dimension, S>(
+          const ArrayView<const S> &,
+          const Table<2, double> &,
+          const std::vector<Vector<deal_II_dimension,
+                                   deal_II_space_dimension>::ShapeFunctionData>
+            &,
+          std::vector<
+            ProductType<S, dealii::Tensor<1, deal_II_space_dimension>>::type>
+            &);
+
+        template void
+        do_function_symmetric_gradients<deal_II_dimension,
+                                        deal_II_space_dimension,
+                                        S>(
+          const ArrayView<const S> &,
+          const Table<2, dealii::Tensor<1, deal_II_space_dimension>> &,
+          const std::vector<Vector<deal_II_dimension,
+                                   deal_II_space_dimension>::ShapeFunctionData>
+            &,
+          std::vector<ProductType<
+            S,
+            dealii::SymmetricTensor<2, deal_II_space_dimension>>::type> &);
+
+        template void
+        do_function_divergences<deal_II_dimension, deal_II_space_dimension, S>(
+          const ArrayView<const S> &,
+          const Table<2, dealii::Tensor<1, deal_II_space_dimension>> &,
+          const std::vector<Vector<deal_II_dimension,
+                                   deal_II_space_dimension>::ShapeFunctionData>
+            &,
+          std::vector<Vector<deal_II_dimension, deal_II_space_dimension>::
+                        template solution_divergence_type<S>> &);
+
+        template void
+        do_function_curls<deal_II_dimension, deal_II_space_dimension, S>(
+          const ArrayView<const S> &,
+          const Table<2, dealii::Tensor<1, deal_II_space_dimension>> &,
+          const std::vector<Vector<deal_II_dimension,
+                                   deal_II_space_dimension>::ShapeFunctionData>
+            &,
+          std::vector<ProductType<
+            S,
+            dealii::internal::CurlType<deal_II_space_dimension>::type>::type>
+            &);
+
+        template void
+        do_function_laplacians<deal_II_dimension, deal_II_space_dimension, S>(
+          const ArrayView<const S> &,
+          const Table<2, dealii::Tensor<2, deal_II_space_dimension>> &,
+          const std::vector<Vector<deal_II_dimension,
+                                   deal_II_space_dimension>::ShapeFunctionData>
+            &,
+          std::vector<Vector<deal_II_dimension, deal_II_space_dimension>::
+                        template solution_laplacian_type<S>> &);
+
+        template void
+        do_function_values<deal_II_dimension, deal_II_space_dimension, S>(
+          const ArrayView<const S> &,
+          const dealii::Table<2, double> &,
+          const std::vector<
+            SymmetricTensor<2, deal_II_dimension, deal_II_space_dimension>::
+              ShapeFunctionData> &,
+          std::vector<ProductType<
+            S,
+            dealii::SymmetricTensor<2, deal_II_space_dimension>>::type> &);
+
+        template void
+        do_function_divergences<deal_II_dimension, deal_II_space_dimension, S>(
+          const ArrayView<const S> &,
+          const Table<2, dealii::Tensor<1, deal_II_space_dimension>> &,
+          const std::vector<
+            SymmetricTensor<2, deal_II_dimension, deal_II_space_dimension>::
+              ShapeFunctionData> &,
+          std::vector<
+            SymmetricTensor<2, deal_II_dimension, deal_II_space_dimension>::
+              template solution_divergence_type<S>> &);
+
+        template void
+        do_function_values<deal_II_dimension, deal_II_space_dimension, S>(
+          const ArrayView<const S> &,
+          const dealii::Table<2, double> &,
+          const std::vector<
+            Tensor<2, deal_II_dimension, deal_II_space_dimension>::
+              ShapeFunctionData> &,
+          std::vector<
+            ProductType<S, dealii::Tensor<2, deal_II_space_dimension>>::type>
+            &);
+
+        template void
+        do_function_divergences<deal_II_dimension, deal_II_space_dimension, S>(
+          const ArrayView<const S> &,
+          const Table<2, dealii::Tensor<1, deal_II_space_dimension>> &,
+          const std::vector<
+            Tensor<2, deal_II_dimension, deal_II_space_dimension>::
+              ShapeFunctionData> &,
+          std::vector<Tensor<2, deal_II_dimension, deal_II_space_dimension>::
+                        template solution_divergence_type<S>> &);
+
+        template void
+        do_function_gradients<deal_II_dimension, deal_II_space_dimension, S>(
+          const ArrayView<const S> &,
+          const Table<2, dealii::Tensor<1, deal_II_space_dimension>> &,
+          const std::vector<
+            Tensor<2, deal_II_dimension, deal_II_space_dimension>::
+              ShapeFunctionData> &,
+          std::vector<Tensor<2, deal_II_dimension, deal_II_space_dimension>::
+                        template solution_gradient_type<S>> &);
+
+      \}
+    \}
+#  endif
+  }
+
+for (S : ALL_SCALAR_TYPES; deal_II_dimension : DIMENSIONS; ORDER : DIMENSIONS;
+     deal_II_space_dimension : SPACE_DIMENSIONS)
+  {
+#  if deal_II_dimension <= deal_II_space_dimension
+    namespace FEValuesViews
+    \{
+      namespace internal
+      \{
+        template void
+        do_function_derivatives<ORDER,
+                                deal_II_dimension,
+                                deal_II_space_dimension,
+                                S>(
+          const ArrayView<const S> &,
+          const Table<2, dealii::Tensor<ORDER, deal_II_space_dimension>> &,
+          const std::vector<Scalar<deal_II_dimension,
+                                   deal_II_space_dimension>::ShapeFunctionData>
+            &,
+          std::vector<
+            ProductType<S,
+                        dealii::Tensor<ORDER, deal_II_space_dimension>>::type>
+            &);
+
+        template void
+        do_function_derivatives<ORDER,
+                                deal_II_dimension,
+                                deal_II_space_dimension,
+                                S>(
+          const ArrayView<const S> &,
+          const Table<2, dealii::Tensor<ORDER, deal_II_space_dimension>> &,
+          const std::vector<Vector<deal_II_dimension,
+                                   deal_II_space_dimension>::ShapeFunctionData>
+            &,
+          std::vector<ProductType<
+            S,
+            dealii::Tensor<ORDER + 1, deal_II_space_dimension>>::type> &);
+      \}
+    \}
+#  endif
+  }
+
+
+#endif

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.