]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Split tensor product points and single point kernels 16736/head
authorMaximilian Bergbauer <maximilian.bergbauer@tum.de>
Fri, 8 Mar 2024 08:22:01 +0000 (09:22 +0100)
committerMaximilian Bergbauer <maximilian.bergbauer@tum.de>
Fri, 8 Mar 2024 12:04:49 +0000 (13:04 +0100)
include/deal.II/fe/mapping_q_internal.h
include/deal.II/matrix_free/fe_point_evaluation.h
include/deal.II/matrix_free/tensor_product_kernels.h
include/deal.II/matrix_free/tensor_product_point_kernels.h [new file with mode: 0644]

index 77bd61c7e936060b9268d5d616a890e449e8dbe7..5bc48229fa717efb277b1f918af6e6ccd873b170 100644 (file)
@@ -37,7 +37,7 @@
 #include <deal.II/matrix_free/fe_evaluation_data.h>
 #include <deal.II/matrix_free/mapping_info_storage.h>
 #include <deal.II/matrix_free/shape_info.h>
-#include <deal.II/matrix_free/tensor_product_kernels.h>
+#include <deal.II/matrix_free/tensor_product_point_kernels.h>
 
 #include <array>
 #include <limits>
index 84303226aa4b3188416ba51d3c609958d7bba319..e7f142616318e8db1790774582d6154c1ea7a8da 100644 (file)
@@ -31,7 +31,7 @@
 #include <deal.II/matrix_free/evaluation_kernels_face.h>
 #include <deal.II/matrix_free/mapping_info_storage.h>
 #include <deal.II/matrix_free/shape_info.h>
-#include <deal.II/matrix_free/tensor_product_kernels.h>
+#include <deal.II/matrix_free/tensor_product_point_kernels.h>
 
 #include <deal.II/non_matching/mapping_info.h>
 
index 90cb883ba2b397589f79252be13f81b192ef7743..4f2319420bde35b3ae99cffb1d763b999d834aa9 100644 (file)
@@ -2014,1572 +2014,6 @@ namespace internal
       }
   }
 
-
-
-  /**
-   * Struct to avoid using Tensor<1, dim, Point<dim2>> in
-   * evaluate_tensor_product_value_and_gradient because a Point cannot be used
-   * within Tensor. Instead, a specialization of this struct upcasts the point
-   * to a Tensor<1,dim>.
-   */
-  template <typename Number, typename Number2>
-  struct ProductTypeNoPoint
-  {
-    using type = typename ProductType<Number, Number2>::type;
-  };
-
-  template <int dim, typename Number, typename Number2>
-  struct ProductTypeNoPoint<Point<dim, Number>, Number2>
-  {
-    using type = typename ProductType<Tensor<1, dim, Number>, Number2>::type;
-  };
-
-
-
-  /**
-   * Computes the values and derivatives of the 1d polynomials @p poly at the
-   * specified point @p p and stores it in @p shapes.
-   */
-  template <int dim, typename Number>
-  inline void
-  compute_values_of_array(
-    dealii::ndarray<Number, 2, dim>                    *shapes,
-    const std::vector<Polynomials::Polynomial<double>> &poly,
-    const Point<dim, Number>                           &p,
-    const unsigned int                                  derivative = 1)
-  {
-    const int n_shapes = poly.size();
-
-    // Evaluate 1d polynomials and their derivatives
-    std::array<Number, dim> point;
-    for (unsigned int d = 0; d < dim; ++d)
-      point[d] = p[d];
-    for (int i = 0; i < n_shapes; ++i)
-      poly[i].values_of_array(point, derivative, shapes[i].data());
-  }
-
-
-
-  /**
-   * Specialization of above function for dim = 0. Should not be called.
-   */
-  template <typename Number>
-  inline void
-  compute_values_of_array(dealii::ndarray<Number, 2, 0> *,
-                          const std::vector<Polynomials::Polynomial<double>> &,
-                          const Point<0, Number> &,
-                          const unsigned int)
-  {
-    DEAL_II_ASSERT_UNREACHABLE();
-  }
-
-
-
-  /**
-   * Interpolate inner dimensions of tensor product shape functions.
-   */
-  template <int dim,
-            int length,
-            typename Number2,
-            typename Number,
-            int  n_values    = 1,
-            bool do_renumber = true,
-            int  stride      = 1>
-  inline
-#ifndef DEBUG
-    DEAL_II_ALWAYS_INLINE
-#endif
-      std::array<typename ProductTypeNoPoint<Number, Number2>::type,
-                 2 + n_values>
-      do_interpolate_xy(const Number                           *values,
-                        const std::vector<unsigned int>        &renumber,
-                        const dealii::ndarray<Number2, 2, dim> *shapes,
-                        const int n_shapes_runtime,
-                        int      &i)
-  {
-    static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented");
-    static_assert(1 <= n_values && n_values <= 2,
-                  "Only n_values=1,2 implemented");
-
-    const int n_shapes = length > 0 ? length : n_shapes_runtime;
-
-    // If n_values > 1, we want to interpolate from a second array,
-    // placed in the same array immediately after the main data. This
-    // is used to interpolate normal derivatives onto faces.
-    const Number *values_2 =
-      n_values > 1 ?
-        values + stride * (length > 0 ?
-                             Utilities::pow(length, dim) :
-                             Utilities::fixed_power<dim>(n_shapes_runtime)) :
-        nullptr;
-    using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
-    std::array<Number3, 2 + n_values> result = {};
-    for (int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1)
-      {
-        // Interpolation + derivative x direction
-        std::array<Number3, 1 + n_values> inner_result = {};
-
-        // Distinguish the inner loop based on whether we have a
-        // renumbering or not
-        if (do_renumber && !renumber.empty())
-          for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
-            {
-              // gradient
-              inner_result[0] +=
-                shapes[i0][1][0] * values[renumber[i] * stride];
-              // values
-              inner_result[1] +=
-                shapes[i0][0][0] * values[renumber[i] * stride];
-              if (n_values > 1)
-                inner_result[2] +=
-                  shapes[i0][0][0] * values_2[renumber[i] * stride];
-            }
-        else
-          for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
-            {
-              // gradient
-              inner_result[0] += shapes[i0][1][0] * values[i * stride];
-              // values
-              inner_result[1] += shapes[i0][0][0] * values[i * stride];
-              if (n_values > 1)
-                inner_result[2] += shapes[i0][0][0] * values_2[i * stride];
-            }
-
-        if (dim > 1)
-          {
-            // Interpolation + derivative in y direction
-            // gradient
-            result[0] += inner_result[0] * shapes[i1][0][1];
-            result[1] += inner_result[1] * shapes[i1][1][1];
-            // values
-            result[2] += inner_result[1] * shapes[i1][0][1];
-            if (n_values > 1)
-              result[3] += inner_result[2] * shapes[i1][0][1];
-          }
-        else
-          {
-            // gradient
-            result[0] = inner_result[0];
-            // values
-            result[1] = inner_result[1];
-            if (n_values > 1)
-              result[2] = inner_result[2];
-          }
-      }
-    return result;
-  }
-
-
-
-  /**
-   * Interpolates the values and gradients into the points specified in
-   * @p compute_values_of_array() with help of the precomputed @p shapes.
-   */
-  template <int dim,
-            typename Number,
-            typename Number2,
-            int  n_values    = 1,
-            bool do_renumber = true,
-            int  stride      = 1>
-  inline std::array<typename ProductTypeNoPoint<Number, Number2>::type,
-                    dim + n_values>
-  evaluate_tensor_product_value_and_gradient_shapes(
-    const dealii::ndarray<Number2, 2, dim> *shapes,
-    const int                               n_shapes,
-    const Number                           *values,
-    const std::vector<unsigned int>        &renumber = {})
-  {
-    static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented");
-    static_assert(1 <= n_values && n_values <= 2,
-                  "Only n_values=1,2 implemented");
-
-    using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
-
-    std::array<Number3, dim + n_values> result = {};
-    if (dim == 0)
-      {
-        // We only need the interpolation of the value and normal derivatives on
-        // faces of a 1d element. As the interpolation is the value at the
-        // point, simply set the result vector accordingly.
-        result[0] = values[0];
-        if (n_values > 1)
-          result[1] = values[1 * stride];
-        return result;
-      }
-
-    // Go through the tensor product of shape functions and interpolate
-    // with optimal algorithm
-    for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
-      {
-        std::array<Number3, 2 + n_values> inner_result;
-        // Generate separate code with known loop bounds for the most common
-        // cases
-        if (n_shapes == 2)
-          inner_result =
-            do_interpolate_xy<dim,
-                              2,
-                              Number2,
-                              Number,
-                              n_values,
-                              do_renumber,
-                              stride>(values, renumber, shapes, n_shapes, i);
-        else if (n_shapes == 3)
-          inner_result =
-            do_interpolate_xy<dim,
-                              3,
-                              Number2,
-                              Number,
-                              n_values,
-                              do_renumber,
-                              stride>(values, renumber, shapes, n_shapes, i);
-        else if (n_shapes == 4)
-          inner_result =
-            do_interpolate_xy<dim,
-                              4,
-                              Number2,
-                              Number,
-                              n_values,
-                              do_renumber,
-                              stride>(values, renumber, shapes, n_shapes, i);
-        else if (n_shapes == 5)
-          inner_result =
-            do_interpolate_xy<dim,
-                              5,
-                              Number2,
-                              Number,
-                              n_values,
-                              do_renumber,
-                              stride>(values, renumber, shapes, n_shapes, i);
-        else if (n_shapes == 6)
-          inner_result =
-            do_interpolate_xy<dim,
-                              6,
-                              Number2,
-                              Number,
-                              n_values,
-                              do_renumber,
-                              stride>(values, renumber, shapes, n_shapes, i);
-        else
-          inner_result =
-            do_interpolate_xy<dim,
-                              -1,
-                              Number2,
-                              Number,
-                              n_values,
-                              do_renumber,
-                              stride>(values, renumber, shapes, n_shapes, i);
-        if (dim == 3)
-          {
-            // derivative + interpolation in z direction
-            // gradient
-            result[0] += inner_result[0] * shapes[i2][0][2];
-            result[1] += inner_result[1] * shapes[i2][0][2];
-            result[2] += inner_result[2] * shapes[i2][1][2];
-            // values
-            result[3] += inner_result[2] * shapes[i2][0][2];
-            if (n_values > 1)
-              result[4] += inner_result[3] * shapes[i2][0][2];
-          }
-        else if (dim == 2)
-          {
-            // gradient
-            result[0] = inner_result[0];
-            result[1] = inner_result[1];
-            // values
-            result[2] = inner_result[2];
-            if (n_values > 1)
-              result[3] = inner_result[3];
-          }
-        else
-          {
-            // gradient
-            result[0] = inner_result[0];
-            // values
-            result[1] = inner_result[1];
-            if (n_values > 1)
-              result[2] = inner_result[2];
-          }
-      }
-
-    return result;
-  }
-
-
-
-  /**
-   * Specializes @p evaluate_tensor_product_value_and_gradient() for linear
-   * polynomials which massively reduces the necessary instructions.
-   */
-  template <int dim,
-            typename Number,
-            typename Number2,
-            int n_values = 1,
-            int stride   = 1>
-  inline std::array<typename ProductTypeNoPoint<Number, Number2>::type,
-                    dim + n_values>
-  evaluate_tensor_product_value_and_gradient_linear(
-    const Number              *values,
-    const Point<dim, Number2> &p)
-  {
-    static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented");
-    static_assert(1 <= n_values && n_values <= 2,
-                  "Only n_values=1,2 implemented");
-
-    using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
-
-    // If n_values > 1, we want to interpolate from a second array,
-    // placed in the same array immediately after the main data. This
-    // is used to interpolate normal derivatives onto faces.
-
-    std::array<Number3, dim + n_values> result;
-    if (dim == 0)
-      {
-        // we only need the value on faces of a 1d element
-        result[0] = values[0];
-        if (n_values > 1)
-          result[1] = values[1 * stride];
-      }
-    else if (dim == 1)
-      {
-        // gradient
-        result[0] = Number3(values[stride] - values[0]);
-        // values
-        result[1] = Number3(values[0]) + p[0] * result[0];
-        if (n_values > 1)
-          result[2] = Number3(values[2 * stride]) +
-                      p[0] * (values[3 * stride] - values[2 * stride]);
-      }
-    else if (dim == 2)
-      {
-        const Number3 val10 = Number3(values[stride] - values[0]);
-        const Number3 val32 = Number3(values[3 * stride] - values[2 * stride]);
-        const Number3 tmp0  = Number3(values[0]) + p[0] * val10;
-        const Number3 tmp1  = Number3(values[2 * stride]) + p[0] * val32;
-
-        // gradient
-        result[0] = val10 + p[1] * (val32 - val10);
-        result[1] = tmp1 - tmp0;
-
-        // values
-        result[2] = tmp0 + p[1] * result[1];
-
-        if (n_values > 1)
-          {
-            const Number3 tmp0_2 =
-              Number3(values[4 * stride]) +
-              p[0] * (values[5 * stride] - values[4 * stride]);
-            const Number3 tmp1_2 =
-              Number3(values[6 * stride]) +
-              p[0] * (values[7 * stride] - values[6 * stride]);
-            result[3] = tmp0_2 + p[1] * (tmp1_2 - tmp0_2);
-          }
-      }
-    else if (dim == 3)
-      {
-        const Number3 val10 = Number3(values[stride] - values[0]);
-        const Number3 val32 = Number3(values[3 * stride] - values[2 * stride]);
-        const Number3 tmp0  = Number3(values[0]) + p[0] * val10;
-        const Number3 tmp1  = Number3(values[2 * stride]) + p[0] * val32;
-        const Number3 tmp10 = tmp1 - tmp0;
-        const Number3 tmpy0 = tmp0 + p[1] * tmp10;
-
-        const Number3 val54 = Number3(values[5 * stride] - values[4 * stride]);
-        const Number3 val76 = Number3(values[7 * stride] - values[6 * stride]);
-        const Number3 tmp2  = Number3(values[4 * stride]) + p[0] * val54;
-        const Number3 tmp3  = Number3(values[6 * stride]) + p[0] * val76;
-        const Number3 tmp32 = tmp3 - tmp2;
-        const Number3 tmpy1 = tmp2 + p[1] * tmp32;
-
-        // gradient
-        result[2]           = tmpy1 - tmpy0;
-        result[1]           = tmp10 + p[2] * (tmp32 - tmp10);
-        const Number3 tmpz0 = val10 + p[1] * (val32 - val10);
-        result[0] = tmpz0 + p[2] * (val54 + p[1] * (val76 - val54) - tmpz0);
-
-        // value
-        result[3] = tmpy0 + p[2] * result[2];
-        Assert(n_values == 1, ExcNotImplemented());
-      }
-
-    return result;
-  }
-
-
-
-  /**
-   * Compute the polynomial interpolation of a tensor product shape function
-   * $\varphi_i$ given a vector of coefficients $u_i$ in the form
-   * $u_h(\mathbf{x}) = \sum_{i=1}^{k^d} \varphi_i(\mathbf{x}) u_i$. The shape
-   * functions $\varphi_i(\mathbf{x}) =
-   * \prod_{d=1}^{\text{dim}}\varphi_{i_d}^\text{1d}(x_d)$ represent a tensor
-   * product. The function returns a pair with the value of the interpolation
-   * as the first component and the gradient in reference coordinates as the
-   * second component. Note that for compound types (e.g. the `values` field
-   * begin a Point<spacedim> argument), the components of the gradient are
-   * sorted as Tensor<1, dim, Tensor<1, spacedim>> with the derivatives
-   * as the first index; this is a consequence of the generic arguments in the
-   * function.
-   *
-   * @param poly The underlying one-dimensional polynomial basis
-   * $\{\varphi^{1d}_{i_1}\}$ given as a vector of polynomials.
-   *
-   * @param values The expansion coefficients $u_i$ of type `Number` in
-   * the polynomial interpolation. The coefficients can be simply `double`
-   * variables but e.g. also Point<spacedim> in case they define arithmetic
-   * operations with the type `Number2`.
-   *
-   * @param p The position in reference coordinates where the interpolation
-   * should be evaluated.
-   *
-   * @param d_linear Flag to specify whether a d-linear (linear in 1d,
-   * bi-linear in 2d, tri-linear in 3d) interpolation should be made, which
-   * allows to unroll loops and considerably speed up evaluation.
-   *
-   * @param renumber Optional parameter to specify a renumbering in the
-   * coefficient vector, assuming that `values[renumber[i]]` returns
-   * the lexicographic (tensor product) entry of the coefficients. If the
-   * vector is entry, the values are assumed to be sorted lexicographically.
-   */
-  template <int dim, typename Number, typename Number2>
-  inline std::pair<
-    typename ProductTypeNoPoint<Number, Number2>::type,
-    Tensor<1, dim, typename ProductTypeNoPoint<Number, Number2>::type>>
-  evaluate_tensor_product_value_and_gradient(
-    const std::vector<Polynomials::Polynomial<double>> &poly,
-    const std::vector<Number>                          &values,
-    const Point<dim, Number2>                          &p,
-    const bool                                          d_linear = false,
-    const std::vector<unsigned int>                    &renumber = {})
-  {
-    using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
-
-    std::array<Number3, dim + 1> result;
-    if (d_linear)
-      {
-        result =
-          evaluate_tensor_product_value_and_gradient_linear(values.data(), p);
-      }
-    else
-      {
-        AssertIndexRange(poly.size(), 200);
-        std::array<dealii::ndarray<Number2, 2, dim>, 200> shapes;
-        compute_values_of_array(shapes.data(), poly, p);
-        result = evaluate_tensor_product_value_and_gradient_shapes<dim,
-                                                                   Number,
-                                                                   Number2>(
-          shapes.data(), poly.size(), values.data(), renumber);
-      }
-    return std::make_pair(result[dim],
-                          Tensor<1, dim, Number3>(
-                            ArrayView<Number3>(result.data(), dim)));
-  }
-
-
-
-  template <int dim,
-            int length,
-            typename Number2,
-            typename Number,
-            bool do_renumber = true,
-            int  stride      = 1>
-  inline
-#ifndef DEBUG
-    DEAL_II_ALWAYS_INLINE
-#endif
-    typename ProductTypeNoPoint<Number, Number2>::type
-    do_interpolate_xy_value(const Number                           *values,
-                            const std::vector<unsigned int>        &renumber,
-                            const dealii::ndarray<Number2, 2, dim> *shapes,
-                            const int n_shapes_runtime,
-                            int      &i)
-  {
-    const int n_shapes = length > 0 ? length : n_shapes_runtime;
-    using Number3      = typename ProductTypeNoPoint<Number, Number2>::type;
-    Number3 result     = {};
-    for (int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1)
-      {
-        // Interpolation x direction
-        Number3 value = {};
-
-        // Distinguish the inner loop based on whether we have a
-        // renumbering or not
-        if (do_renumber && !renumber.empty())
-          for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
-            value += shapes[i0][0][0] * values[renumber[i] * stride];
-        else
-          for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
-            value += shapes[i0][0][0] * values[i * stride];
-
-        if (dim > 1)
-          result += value * shapes[i1][0][1];
-        else
-          result = value;
-      }
-    return result;
-  }
-
-
-
-  template <int dim,
-            typename Number,
-            typename Number2,
-            bool do_renumber = true,
-            int  stride      = 1>
-  inline typename ProductTypeNoPoint<Number, Number2>::type
-  evaluate_tensor_product_value_shapes(
-    const dealii::ndarray<Number2, 2, dim> *shapes,
-    const int                               n_shapes,
-    const Number                           *values,
-    const std::vector<unsigned int>        &renumber = {})
-  {
-    static_assert(dim >= 0 && dim <= 3, "Only dim=0,1,2,3 implemented");
-
-    // we only need the value on faces of a 1d element
-    if (dim == 0)
-      {
-        return values[0];
-      }
-
-    using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
-
-    // Go through the tensor product of shape functions and interpolate
-    // with optimal algorithm
-    Number3 result = {};
-    for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
-      {
-        Number3 inner_result;
-        // Generate separate code with known loop bounds for the most common
-        // cases
-        if (n_shapes == 2)
-          inner_result = do_interpolate_xy_value<dim,
-                                                 2,
-                                                 Number2,
-                                                 Number,
-                                                 do_renumber,
-                                                 stride>(
-            values, renumber, shapes, n_shapes, i);
-        else if (n_shapes == 3)
-          inner_result = do_interpolate_xy_value<dim,
-                                                 3,
-                                                 Number2,
-                                                 Number,
-                                                 do_renumber,
-                                                 stride>(
-            values, renumber, shapes, n_shapes, i);
-        else if (n_shapes == 4)
-          inner_result = do_interpolate_xy_value<dim,
-                                                 4,
-                                                 Number2,
-                                                 Number,
-                                                 do_renumber,
-                                                 stride>(
-            values, renumber, shapes, n_shapes, i);
-        else if (n_shapes == 5)
-          inner_result = do_interpolate_xy_value<dim,
-                                                 5,
-                                                 Number2,
-                                                 Number,
-                                                 do_renumber,
-                                                 stride>(
-            values, renumber, shapes, n_shapes, i);
-        else if (n_shapes == 6)
-          inner_result = do_interpolate_xy_value<dim,
-                                                 6,
-                                                 Number2,
-                                                 Number,
-                                                 do_renumber,
-                                                 stride>(
-            values, renumber, shapes, n_shapes, i);
-        else
-          inner_result = do_interpolate_xy_value<dim,
-                                                 -1,
-                                                 Number2,
-                                                 Number,
-                                                 do_renumber,
-                                                 stride>(
-            values, renumber, shapes, n_shapes, i);
-        if (dim == 3)
-          {
-            // Interpolation + derivative in z direction
-            result += inner_result * shapes[i2][0][2];
-          }
-        else
-          {
-            result = inner_result;
-          }
-      }
-
-    return result;
-  }
-
-
-
-  template <int dim, typename Number, typename Number2, int stride = 1>
-  inline typename ProductTypeNoPoint<Number, Number2>::type
-  evaluate_tensor_product_value_linear(const Number              *values,
-                                       const Point<dim, Number2> &p)
-  {
-    static_assert(dim >= 0 && dim <= 3, "Only dim=0,1,2,3 implemented");
-
-    using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
-
-    if (dim == 0)
-      {
-        // we only need the value on faces of a 1d element
-        return values[0];
-      }
-    else if (dim == 1)
-      {
-        return Number3(values[0]) + p[0] * Number3(values[stride] - values[0]);
-      }
-    else if (dim == 2)
-      {
-        const Number3 val10 = Number3(values[stride] - values[0]);
-        const Number3 val32 = Number3(values[3 * stride] - values[2 * stride]);
-        const Number3 tmp0  = Number3(values[0]) + p[0] * val10;
-        const Number3 tmp1  = Number3(values[2 * stride]) + p[0] * val32;
-        return tmp0 + p[1] * (tmp1 - tmp0);
-      }
-    else if (dim == 3)
-      {
-        const Number3 val10 = Number3(values[stride] - values[0]);
-        const Number3 val32 = Number3(values[3 * stride] - values[2 * stride]);
-        const Number3 tmp0  = Number3(values[0]) + p[0] * val10;
-        const Number3 tmp1  = Number3(values[2 * stride]) + p[0] * val32;
-        const Number3 tmpy0 = tmp0 + p[1] * (tmp1 - tmp0);
-
-        const Number3 val54 = Number3(values[5 * stride] - values[4 * stride]);
-        const Number3 val76 = Number3(values[7 * stride] - values[6 * stride]);
-        const Number3 tmp2  = Number3(values[4 * stride]) + p[0] * val54;
-        const Number3 tmp3  = Number3(values[6 * stride]) + p[0] * val76;
-        const Number3 tmpy1 = tmp2 + p[1] * (tmp3 - tmp2);
-
-        return tmpy0 + p[2] * (tmpy1 - tmpy0);
-      }
-
-    // work around a compile error: missing return statement
-    return Number3();
-  }
-
-
-
-  template <int dim, typename Number, typename Number2>
-  inline typename ProductTypeNoPoint<Number, Number2>::type
-  evaluate_tensor_product_value(
-    const std::vector<Polynomials::Polynomial<double>> &poly,
-    const std::vector<Number>                          &values,
-    const Point<dim, Number2>                          &p,
-    const bool                                          d_linear = false,
-    const std::vector<unsigned int>                    &renumber = {})
-  {
-    typename ProductTypeNoPoint<Number, Number2>::type result;
-    if (d_linear)
-      {
-        result = evaluate_tensor_product_value_linear(values.data(), p);
-      }
-    else
-      {
-        AssertIndexRange(poly.size(), 200);
-        std::array<dealii::ndarray<Number2, 2, dim>, 200> shapes;
-        const int                n_shapes = poly.size();
-        std::array<Number2, dim> point;
-        for (unsigned int d = 0; d < dim; ++d)
-          point[d] = p[d];
-        for (int i = 0; i < n_shapes; ++i)
-          poly[i].values_of_array(point, 0, shapes[i].data());
-        result = evaluate_tensor_product_value_shapes<dim, Number, Number2>(
-          shapes.data(), n_shapes, values.data(), renumber);
-      }
-    return result;
-  }
-
-
-
-  /**
-   * This function computes derivatives of arbitrary orders in 1d, returning a
-   * Tensor with the respective derivative
-   */
-  template <int derivative_order, typename Number, typename Number2>
-  inline Tensor<1, 1, typename ProductTypeNoPoint<Number, Number2>::type>
-  evaluate_tensor_product_higher_derivatives(
-    const std::vector<Polynomials::Polynomial<double>> &poly,
-    const std::vector<Number>                          &values,
-    const Point<1, Number2>                            &p,
-    const std::vector<unsigned int>                    &renumber = {})
-  {
-    using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
-
-    const int n_shapes = poly.size();
-    AssertDimension(n_shapes, values.size());
-    Assert(renumber.empty() || renumber.size() == values.size(),
-           ExcDimensionMismatch(renumber.size(), values.size()));
-
-    std::array<Number2, derivative_order + 1> shapes;
-    Tensor<1, 1, Number3>                     result;
-    if (renumber.empty())
-      for (int i = 0; i < n_shapes; ++i)
-        {
-          poly[i].value(p[0], derivative_order, shapes.data());
-          result[0] += shapes[derivative_order] * values[i];
-        }
-    else
-      for (int i = 0; i < n_shapes; ++i)
-        {
-          poly[i].value(p[0], derivative_order, shapes.data());
-          result[0] += shapes[derivative_order] * values[renumber[i]];
-        }
-    return result;
-  }
-
-
-
-  /**
-   * This function computes derivatives of arbitrary orders in 2d, returning a
-   * Tensor with the respective derivatives
-   */
-  template <int derivative_order, typename Number, typename Number2>
-  inline Tensor<1,
-                derivative_order + 1,
-                typename ProductTypeNoPoint<Number, Number2>::type>
-  evaluate_tensor_product_higher_derivatives(
-    const std::vector<Polynomials::Polynomial<double>> &poly,
-    const std::vector<Number>                          &values,
-    const Point<2, Number2>                            &p,
-    const std::vector<unsigned int>                    &renumber = {})
-  {
-    using Number3     = typename ProductTypeNoPoint<Number, Number2>::type;
-    constexpr int dim = 2;
-
-    const int n_shapes = poly.size();
-    AssertDimension(Utilities::pow(n_shapes, 2), values.size());
-    Assert(renumber.empty() || renumber.size() == values.size(),
-           ExcDimensionMismatch(renumber.size(), values.size()));
-
-    AssertIndexRange(n_shapes, 100);
-    dealii::ndarray<Number2, 100, derivative_order + 1, dim> shapes;
-    // Evaluate 1d polynomials and their derivatives
-    std::array<Number2, dim> point;
-    for (unsigned int d = 0; d < dim; ++d)
-      point[d] = p[d];
-    for (int i = 0; i < n_shapes; ++i)
-      poly[i].values_of_array(point, derivative_order, &shapes[i][0]);
-
-    Tensor<1, derivative_order + 1, Number3> result;
-    for (int i1 = 0, i = 0; i1 < n_shapes; ++i1)
-      {
-        Tensor<1, derivative_order + 1, Number3> result_x;
-        if (renumber.empty())
-          for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
-            for (unsigned int d = 0; d <= derivative_order; ++d)
-              result_x[d] += shapes[i0][d][0] * values[i];
-        else
-          for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
-            for (unsigned int d = 0; d <= derivative_order; ++d)
-              result_x[d] += shapes[i0][d][0] * values[renumber[i]];
-
-        for (unsigned int d = 0; d <= derivative_order; ++d)
-          result[d] += shapes[i1][d][1] * result_x[derivative_order - d];
-      }
-    return result;
-  }
-
-
-
-  /**
-   * This function computes derivatives of arbitrary orders in 3d, returning a
-   * Tensor with the respective derivatives
-   */
-  template <int derivative_order, typename Number, typename Number2>
-  inline Tensor<1,
-                ((derivative_order + 1) * (derivative_order + 2)) / 2,
-                typename ProductTypeNoPoint<Number, Number2>::type>
-  evaluate_tensor_product_higher_derivatives(
-    const std::vector<Polynomials::Polynomial<double>> &poly,
-    const std::vector<Number>                          &values,
-    const Point<3, Number2>                            &p,
-    const std::vector<unsigned int>                    &renumber = {})
-  {
-    using Number3     = typename ProductTypeNoPoint<Number, Number2>::type;
-    constexpr int dim = 3;
-    constexpr int n_derivatives =
-      ((derivative_order + 1) * (derivative_order + 2)) / 2;
-
-    const int n_shapes = poly.size();
-    AssertDimension(Utilities::pow(n_shapes, 3), values.size());
-    Assert(renumber.empty() || renumber.size() == values.size(),
-           ExcDimensionMismatch(renumber.size(), values.size()));
-
-    AssertIndexRange(n_shapes, 100);
-    dealii::ndarray<Number2, 100, derivative_order + 1, dim> shapes;
-    // Evaluate 1d polynomials and their derivatives
-    std::array<Number2, dim> point;
-    for (unsigned int d = 0; d < dim; ++d)
-      point[d] = p[d];
-    for (int i = 0; i < n_shapes; ++i)
-      poly[i].values_of_array(point, derivative_order, &shapes[i][0]);
-
-    Tensor<1, n_derivatives, Number3> result;
-    for (int i2 = 0, i = 0; i2 < n_shapes; ++i2)
-      {
-        Tensor<1, n_derivatives, Number3> result_xy;
-        for (int i1 = 0; i1 < n_shapes; ++i1)
-          {
-            // apply x derivatives
-            Tensor<1, derivative_order + 1, Number3> result_x;
-            if (renumber.empty())
-              for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
-                for (unsigned int d = 0; d <= derivative_order; ++d)
-                  result_x[d] += shapes[i0][d][0] * values[i];
-            else
-              for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
-                for (unsigned int d = 0; d <= derivative_order; ++d)
-                  result_x[d] += shapes[i0][d][0] * values[renumber[i]];
-
-            // multiply by y derivatives, sorting them in upper triangular
-            // matrix, starting with highest global derivative order,
-            // decreasing the combined order of xy derivatives by one in each
-            // row, to be combined with z derivatives in the next step
-            for (unsigned int d = 0, c = 0; d <= derivative_order; ++d)
-              for (unsigned int e = d; e <= derivative_order; ++e, ++c)
-                result_xy[c] +=
-                  shapes[i1][e - d][1] * result_x[derivative_order - e];
-          }
-
-        // multiply by z derivatives, starting with highest x derivative
-        for (unsigned int d = 0, c = 0; d <= derivative_order; ++d)
-          for (unsigned int e = d; e <= derivative_order; ++e, ++c)
-            result[c] += shapes[i2][d][2] * result_xy[c];
-      }
-    return result;
-  }
-
-
-
-  template <int dim, typename Number, typename Number2>
-  SymmetricTensor<2, dim, typename ProductTypeNoPoint<Number, Number2>::type>
-  evaluate_tensor_product_hessian(
-    const std::vector<Polynomials::Polynomial<double>> &poly,
-    const std::vector<Number>                          &values,
-    const Point<dim, Number2>                          &p,
-    const std::vector<unsigned int>                    &renumber = {})
-  {
-    static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented");
-
-    const auto hessian =
-      evaluate_tensor_product_higher_derivatives<2>(poly, values, p, renumber);
-
-    using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
-    SymmetricTensor<2, dim, Number3> result;
-    if (dim == 1)
-      result[0][0] = hessian[0];
-    else if (dim >= 2)
-      {
-        // derivatives in Hessian are xx, xy, yy, xz, yz, zz, so must re-order
-        // them for 3D
-        for (unsigned int d = 0, c = 0; d < 2; ++d)
-          for (unsigned int e = d; e < 2; ++e, ++c)
-            result[d][e] = hessian[c];
-        if (dim == 3)
-          {
-            for (unsigned int d = 0; d < 2; ++d)
-              result[d][2] = hessian[3 + d];
-            result[2][2] = hessian[5];
-          }
-      }
-
-    return result;
-  }
-
-
-
-  /**
-   * Test inner dimensions of tensor product shape functions and accumulate.
-   */
-  template <int dim,
-            int length,
-            typename Number2,
-            typename Number,
-            bool add,
-            int  n_values = 1>
-  inline
-#ifndef DEBUG
-    DEAL_II_ALWAYS_INLINE
-#endif
-    void
-    do_apply_test_functions_xy(
-      Number2                                 *values,
-      const dealii::ndarray<Number, 2, dim>   *shapes,
-      const std::array<Number2, 2 + n_values> &test_grads_value,
-      const int                                n_shapes_runtime,
-      int                                     &i)
-  {
-    static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented");
-    static_assert(1 <= n_values && n_values <= 2,
-                  "Only n_values=1,2 implemented");
-
-    // Note that 'add' is a template argument, so the compiler will remove
-    // these checks
-    if (length > 0)
-      {
-        constexpr unsigned int         array_size = length > 0 ? length : 1;
-        std::array<Number, array_size> shape_values_x;
-        std::array<Number, array_size> shape_derivs_x;
-        for (unsigned int j = 0; j < array_size; ++j)
-          {
-            shape_values_x[j] = shapes[j][0][0];
-            shape_derivs_x[j] = shapes[j][1][0];
-          }
-        for (int i1 = 0; i1 < (dim > 1 ? length : 1); ++i1)
-          {
-            const Number2 test_value_y =
-              dim > 1 ? (test_grads_value[2] * shapes[i1][0][1] +
-                         test_grads_value[1] * shapes[i1][1][1]) :
-                        test_grads_value[2];
-            const Number2 test_grad_xy =
-              dim > 1 ? test_grads_value[0] * shapes[i1][0][1] :
-                        test_grads_value[0];
-            Number2 test_value_y_2;
-            if (n_values > 1)
-              test_value_y_2 = dim > 1 ?
-                                 test_grads_value[3] * shapes[i1][0][1] :
-                                 test_grads_value[3];
-
-            Number2 *values_ptr = values + i + i1 * length;
-            Number2 *values_ptr_2 =
-              n_values > 1 ? values_ptr + Utilities::pow(length, dim) : nullptr;
-            for (int i0 = 0; i0 < length; ++i0)
-              {
-                if (add)
-                  values_ptr[i0] += shape_values_x[i0] * test_value_y;
-                else
-                  values_ptr[i0] = shape_values_x[i0] * test_value_y;
-                values_ptr[i0] += shape_derivs_x[i0] * test_grad_xy;
-                if (n_values > 1)
-                  {
-                    if (add)
-                      values_ptr_2[i0] += shape_values_x[i0] * test_value_y_2;
-                    else
-                      values_ptr_2[i0] = shape_values_x[i0] * test_value_y_2;
-                  }
-              }
-          }
-        i += (dim > 1 ? length * length : length);
-      }
-    else
-      {
-        for (int i1 = 0; i1 < (dim > 1 ? n_shapes_runtime : 1); ++i1)
-          {
-            const Number2 test_value_y =
-              dim > 1 ? (test_grads_value[2] * shapes[i1][0][1] +
-                         test_grads_value[1] * shapes[i1][1][1]) :
-                        test_grads_value[2];
-            const Number2 test_grad_xy =
-              dim > 1 ? test_grads_value[0] * shapes[i1][0][1] :
-                        test_grads_value[0];
-            Number2 test_value_y_2;
-            if (n_values > 1)
-              test_value_y_2 = dim > 1 ?
-                                 test_grads_value[3] * shapes[i1][0][1] :
-                                 test_grads_value[3];
-
-            Number2 *values_ptr = values + i + i1 * n_shapes_runtime;
-            Number2 *values_ptr_2 =
-              n_values > 1 ?
-                values_ptr + Utilities::fixed_power<dim>(n_shapes_runtime) :
-                nullptr;
-            for (int i0 = 0; i0 < n_shapes_runtime; ++i0)
-              {
-                if (add)
-                  values_ptr[i0] += shapes[i0][0][0] * test_value_y;
-                else
-                  values_ptr[i0] = shapes[i0][0][0] * test_value_y;
-                values_ptr[i0] += shapes[i0][1][0] * test_grad_xy;
-                if (n_values > 1)
-                  {
-                    if (add)
-                      values_ptr_2[i0] += shapes[i0][0][0] * test_value_y_2;
-                    else
-                      values_ptr_2[i0] = shapes[i0][0][0] * test_value_y_2;
-                  }
-              }
-          }
-        i += (dim > 1 ? n_shapes_runtime * n_shapes_runtime : n_shapes_runtime);
-      }
-  }
-
-
-
-  /**
-   * Same as evaluate_tensor_product_value_and_gradient_shapes() but for
-   * integration.
-   */
-  template <int dim,
-            typename Number,
-            typename Number2,
-            bool add,
-            int  n_values = 1>
-  inline void
-  integrate_add_tensor_product_value_and_gradient_shapes(
-    const dealii::ndarray<Number, 2, dim> *shapes,
-    const int                              n_shapes,
-    const Number2                         *value,
-    const Tensor<1, dim, Number2>         &gradient,
-    Number2                               *values)
-  {
-    static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented");
-    static_assert(1 <= n_values && n_values <= 2,
-                  "Only n_values=1,2 implemented");
-
-    // Note that 'add' is a template argument, so the compiler will remove
-    // these checks
-    if (dim == 0)
-      {
-        if (add)
-          values[0] += value[0];
-        else
-          values[0] = value[0];
-        if (n_values > 1)
-          {
-            if (add)
-              values[1] += value[1];
-            else
-              values[1] = value[1];
-          }
-        return;
-      }
-
-    // Implement the transpose of the function above
-    // as in evaluate, use `int` type to produce better code in this context
-    std::array<Number2, 2 + n_values> test_grads_value;
-    for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
-      {
-        // test grad x
-        test_grads_value[0] =
-          dim > 2 ? gradient[0] * shapes[i2][0][2] : gradient[0];
-        // test grad y
-        test_grads_value[1] = dim > 2 ? gradient[1] * shapes[i2][0][2] :
-                                        (dim > 1 ? gradient[1] : Number2());
-        // test value z
-        test_grads_value[2] =
-          dim > 2 ?
-            (value[0] * shapes[i2][0][2] + gradient[2] * shapes[i2][1][2]) :
-            value[0];
-
-        if (n_values > 1)
-          test_grads_value[3] =
-            dim > 2 ? value[1] * shapes[i2][0][2] : value[1];
-        // Generate separate code with known loop bounds for the most common
-        // cases
-        if (n_shapes == 2)
-          do_apply_test_functions_xy<dim, 2, Number2, Number, add, n_values>(
-            values, shapes, test_grads_value, n_shapes, i);
-        else if (n_shapes == 3)
-          do_apply_test_functions_xy<dim, 3, Number2, Number, add, n_values>(
-            values, shapes, test_grads_value, n_shapes, i);
-        else if (n_shapes == 4)
-          do_apply_test_functions_xy<dim, 4, Number2, Number, add, n_values>(
-            values, shapes, test_grads_value, n_shapes, i);
-        else if (n_shapes == 5)
-          do_apply_test_functions_xy<dim, 5, Number2, Number, add, n_values>(
-            values, shapes, test_grads_value, n_shapes, i);
-        else if (n_shapes == 6)
-          do_apply_test_functions_xy<dim, 6, Number2, Number, add, n_values>(
-            values, shapes, test_grads_value, n_shapes, i);
-        else
-          do_apply_test_functions_xy<dim, -1, Number2, Number, add, n_values>(
-            values, shapes, test_grads_value, n_shapes, i);
-      }
-  }
-
-
-
-  /**
-   * Specializes @p integrate_add_tensor_product_value_and_gradient_shapes() for linear
-   * polynomials which massively reduces the necessary instructions.
-   */
-  template <int dim,
-            typename Number,
-            typename Number2,
-            bool add,
-            int  n_values = 1>
-  inline void
-  integrate_add_tensor_product_value_and_gradient_linear(
-    const Number2                 *value,
-    const Tensor<1, dim, Number2> &gradient,
-    Number2                       *values,
-    const Point<dim, Number>      &p)
-  {
-    static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented");
-    static_assert(1 <= n_values && n_values <= 2,
-                  "Only n_values=1,2 implemented");
-
-    // Note that 'add' is a template argument, so the compiler will remove
-    // these checks
-    if (dim == 0)
-      {
-        if (add)
-          values[0] += value[0];
-        else
-          values[0] = value[0];
-        if (n_values > 1)
-          {
-            if (add)
-              values[1] += value[1];
-            else
-              values[1] = value[1];
-          }
-      }
-    else if (dim == 1)
-      {
-        const Number2 difference = value[0] * p[0] + gradient[0];
-        if (add)
-          {
-            values[0] += value[0] - difference;
-            values[1] += difference;
-          }
-        else
-          {
-            values[0] = value[0] - difference;
-            values[1] = difference;
-          }
-        if (n_values > 1)
-          {
-            const Number2 product = value[1] * p[0];
-            if (add)
-              {
-                values[2] += value[1] - product;
-                values[3] += product;
-              }
-            else
-              {
-                values[2] = value[1] - product;
-                values[3] = product;
-              }
-          }
-      }
-    else if (dim == 2)
-      {
-        const Number2 test_value_y1 = value[0] * p[1] + gradient[1];
-        const Number2 test_value_y0 = value[0] - test_value_y1;
-        const Number2 test_grad_xy1 = gradient[0] * p[1];
-        const Number2 test_grad_xy0 = gradient[0] - test_grad_xy1;
-        const Number2 value0        = p[0] * test_value_y0 + test_grad_xy0;
-        const Number2 value1        = p[0] * test_value_y1 + test_grad_xy1;
-
-        if (add)
-          {
-            values[0] += test_value_y0 - value0;
-            values[1] += value0;
-            values[2] += test_value_y1 - value1;
-            values[3] += value1;
-          }
-        else
-          {
-            values[0] = test_value_y0 - value0;
-            values[1] = value0;
-            values[2] = test_value_y1 - value1;
-            values[3] = value1;
-          }
-
-        if (n_values > 1)
-          {
-            const Number2 test_value_y1_2 = value[1] * p[1];
-            const Number2 test_value_y0_2 = value[1] - test_value_y1_2;
-            const Number2 value0_2        = p[0] * test_value_y0_2;
-            const Number2 value1_2        = p[0] * test_value_y1_2;
-
-            if (add)
-              {
-                values[4] += test_value_y0_2 - value0_2;
-                values[5] += value0_2;
-                values[6] += test_value_y1_2 - value1_2;
-                values[7] += value1_2;
-              }
-            else
-              {
-                values[4] = test_value_y0_2 - value0_2;
-                values[5] = value0_2;
-                values[6] = test_value_y1_2 - value1_2;
-                values[7] = value1_2;
-              }
-          }
-      }
-    else if (dim == 3)
-      {
-        Assert(n_values == 1, ExcNotImplemented());
-
-        const Number2 test_value_z1 = value[0] * p[2] + gradient[2];
-        const Number2 test_value_z0 = value[0] - test_value_z1;
-        const Number2 test_grad_x1  = gradient[0] * p[2];
-        const Number2 test_grad_x0  = gradient[0] - test_grad_x1;
-        const Number2 test_grad_y1  = gradient[1] * p[2];
-        const Number2 test_grad_y0  = gradient[1] - test_grad_y1;
-
-        const Number2 test_value_y01 = test_value_z0 * p[1] + test_grad_y0;
-        const Number2 test_value_y00 = test_value_z0 - test_value_y01;
-        const Number2 test_grad_xy01 = test_grad_x0 * p[1];
-        const Number2 test_grad_xy00 = test_grad_x0 - test_grad_xy01;
-        const Number2 test_value_y11 = test_value_z1 * p[1] + test_grad_y1;
-        const Number2 test_value_y10 = test_value_z1 - test_value_y11;
-        const Number2 test_grad_xy11 = test_grad_x1 * p[1];
-        const Number2 test_grad_xy10 = test_grad_x1 - test_grad_xy11;
-
-        const Number2 value00 = p[0] * test_value_y00 + test_grad_xy00;
-        const Number2 value01 = p[0] * test_value_y01 + test_grad_xy01;
-        const Number2 value10 = p[0] * test_value_y10 + test_grad_xy10;
-        const Number2 value11 = p[0] * test_value_y11 + test_grad_xy11;
-
-        if (add)
-          {
-            values[0] += test_value_y00 - value00;
-            values[1] += value00;
-            values[2] += test_value_y01 - value01;
-            values[3] += value01;
-            values[4] += test_value_y10 - value10;
-            values[5] += value10;
-            values[6] += test_value_y11 - value11;
-            values[7] += value11;
-          }
-        else
-          {
-            values[0] = test_value_y00 - value00;
-            values[1] = value00;
-            values[2] = test_value_y01 - value01;
-            values[3] = value01;
-            values[4] = test_value_y10 - value10;
-            values[5] = value10;
-            values[6] = test_value_y11 - value11;
-            values[7] = value11;
-          }
-      }
-  }
-
-
-
-  /**
-   * Calls the correct @p integrate_add_tensor_product_value_and_gradient_...()
-   * function depending on if values should be added to or set and if
-   * polynomials are linear.
-   */
-  template <bool is_linear,
-            int  dim,
-            typename Number,
-            typename Number2,
-            int n_values = 1>
-  inline void
-  integrate_tensor_product_value_and_gradient(
-    const dealii::ndarray<Number, 2, dim> *shapes,
-    const unsigned int                     n_shapes,
-    const Number2                         *value,
-    const Tensor<1, dim, Number2>         &gradient,
-    Number2                               *values,
-    const Point<dim, Number>              &p,
-    const bool                             do_add)
-  {
-    if (do_add)
-      {
-        if (is_linear)
-          internal::integrate_add_tensor_product_value_and_gradient_linear<
-            dim,
-            Number,
-            Number2,
-            true,
-            n_values>(value, gradient, values, p);
-        else
-          internal::integrate_add_tensor_product_value_and_gradient_shapes<
-            dim,
-            Number,
-            Number2,
-            true,
-            n_values>(shapes, n_shapes, value, gradient, values);
-      }
-    else
-      {
-        if (is_linear)
-          internal::integrate_add_tensor_product_value_and_gradient_linear<
-            dim,
-            Number,
-            Number2,
-            false,
-            n_values>(value, gradient, values, p);
-        else
-          internal::integrate_add_tensor_product_value_and_gradient_shapes<
-            dim,
-            Number,
-            Number2,
-            false,
-            n_values>(shapes, n_shapes, value, gradient, values);
-      }
-  }
-
-
-
-  /**
-   * Test inner dimensions of tensor product shape functions and accumulate.
-   */
-  template <int dim, int length, typename Number2, typename Number, bool add>
-  inline
-#ifndef DEBUG
-    DEAL_II_ALWAYS_INLINE
-#endif
-    void
-    do_apply_test_functions_xy_value(
-      Number2                               *values,
-      const dealii::ndarray<Number, 2, dim> *shapes,
-      const Number2                         &test_value,
-      const int                              n_shapes_runtime,
-      int                                   &i)
-  {
-    if (length > 0)
-      {
-        constexpr unsigned int         array_size = length > 0 ? length : 1;
-        std::array<Number, array_size> shape_values_x;
-        for (unsigned int i1 = 0; i1 < array_size; ++i1)
-          shape_values_x[i1] = shapes[i1][0][0];
-        for (unsigned int i1 = 0; i1 < (dim > 1 ? length : 1); ++i1)
-          {
-            const Number2 test_value_y =
-              dim > 1 ? test_value * shapes[i1][0][1] : test_value;
-
-            Number2 *values_ptr = values + i + i1 * length;
-            for (unsigned int i0 = 0; i0 < length; ++i0)
-              {
-                if (add)
-                  values_ptr[i0] += shape_values_x[i0] * test_value_y;
-                else
-                  values_ptr[i0] = shape_values_x[i0] * test_value_y;
-              }
-          }
-        i += (dim > 1 ? length * length : length);
-      }
-    else
-      {
-        for (int i1 = 0; i1 < (dim > 1 ? n_shapes_runtime : 1); ++i1)
-          {
-            const Number2 test_value_y =
-              dim > 1 ? test_value * shapes[i1][0][1] : test_value;
-
-            Number2 *values_ptr = values + i + i1 * n_shapes_runtime;
-            for (int i0 = 0; i0 < n_shapes_runtime; ++i0)
-              {
-                if (add)
-                  values_ptr[i0] += shapes[i0][0][0] * test_value_y;
-                else
-                  values_ptr[i0] = shapes[i0][0][0] * test_value_y;
-              }
-          }
-        i += (dim > 1 ? n_shapes_runtime * n_shapes_runtime : n_shapes_runtime);
-      }
-  }
-
-
-
-  /**
-   * Same as evaluate_tensor_product_value_shapes() but for integration.
-   */
-  template <int dim, typename Number, typename Number2, bool add>
-  inline void
-  integrate_add_tensor_product_value_shapes(
-    const dealii::ndarray<Number, 2, dim> *shapes,
-    const int                              n_shapes,
-    const Number2                         &value,
-    Number2                               *values)
-  {
-    static_assert(dim >= 0 && dim <= 3, "Only dim=0,1,2,3 implemented");
-
-    // as in evaluate, use `int` type to produce better code in this context
-
-    if (dim == 0)
-      {
-        if (add)
-          values[0] += value;
-        else
-          values[0] = value;
-        return;
-      }
-
-    // Implement the transpose of the function above
-    Number2 test_value;
-    for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
-      {
-        // test value z
-        test_value = dim > 2 ? value * shapes[i2][0][2] : value;
-
-        // Generate separate code with known loop bounds for the most common
-        // cases
-        if (n_shapes == 2)
-          do_apply_test_functions_xy_value<dim, 2, Number2, Number, add>(
-            values, shapes, test_value, n_shapes, i);
-        else if (n_shapes == 3)
-          do_apply_test_functions_xy_value<dim, 3, Number2, Number, add>(
-            values, shapes, test_value, n_shapes, i);
-        else if (n_shapes == 4)
-          do_apply_test_functions_xy_value<dim, 4, Number2, Number, add>(
-            values, shapes, test_value, n_shapes, i);
-        else if (n_shapes == 5)
-          do_apply_test_functions_xy_value<dim, 5, Number2, Number, add>(
-            values, shapes, test_value, n_shapes, i);
-        else if (n_shapes == 6)
-          do_apply_test_functions_xy_value<dim, 6, Number2, Number, add>(
-            values, shapes, test_value, n_shapes, i);
-        else
-          do_apply_test_functions_xy_value<dim, -1, Number2, Number, add>(
-            values, shapes, test_value, n_shapes, i);
-      }
-  }
-
-
-
-  /**
-   * Specializes @p integrate_tensor_product_value_shapes() for linear
-   * polynomials which massively reduces the necessary instructions.
-   */
-  template <int dim, typename Number, typename Number2, bool add>
-  inline void
-  integrate_add_tensor_product_value_linear(const Number2            &value,
-                                            Number2                  *values,
-                                            const Point<dim, Number> &p)
-  {
-    static_assert(dim >= 0 && dim <= 3, "Only dim=0,1,2,3 implemented");
-
-    if (dim == 0)
-      {
-        if (add)
-          values[0] += value;
-        else
-          values[0] = value;
-      }
-    else if (dim == 1)
-      {
-        const auto x0 = 1. - p[0], x1 = p[0];
-
-        if (add)
-          {
-            values[0] += value * x0;
-            values[1] += value * x1;
-          }
-        else
-          {
-            values[0] = value * x0;
-            values[1] = value * x1;
-          }
-      }
-    else if (dim == 2)
-      {
-        const auto x0 = 1. - p[0], x1 = p[0], y0 = 1. - p[1], y1 = p[1];
-
-        const auto test_value_y0 = value * y0;
-        const auto test_value_y1 = value * y1;
-
-        if (add)
-          {
-            values[0] += x0 * test_value_y0;
-            values[1] += x1 * test_value_y0;
-            values[2] += x0 * test_value_y1;
-            values[3] += x1 * test_value_y1;
-          }
-        else
-          {
-            values[0] = x0 * test_value_y0;
-            values[1] = x1 * test_value_y0;
-            values[2] = x0 * test_value_y1;
-            values[3] = x1 * test_value_y1;
-          }
-      }
-    else if (dim == 3)
-      {
-        const auto x0 = 1. - p[0], x1 = p[0], y0 = 1. - p[1], y1 = p[1],
-                   z0 = 1. - p[2], z1 = p[2];
-
-        const auto test_value_z0 = value * z0;
-        const auto test_value_z1 = value * z1;
-
-        const auto test_value_y00 = test_value_z0 * y0;
-        const auto test_value_y01 = test_value_z0 * y1;
-        const auto test_value_y10 = test_value_z1 * y0;
-        const auto test_value_y11 = test_value_z1 * y1;
-
-        if (add)
-          {
-            values[0] += x0 * test_value_y00;
-            values[1] += x1 * test_value_y00;
-            values[2] += x0 * test_value_y01;
-            values[3] += x1 * test_value_y01;
-            values[4] += x0 * test_value_y10;
-            values[5] += x1 * test_value_y10;
-            values[6] += x0 * test_value_y11;
-            values[7] += x1 * test_value_y11;
-          }
-        else
-          {
-            values[0] = x0 * test_value_y00;
-            values[1] = x1 * test_value_y00;
-            values[2] = x0 * test_value_y01;
-            values[3] = x1 * test_value_y01;
-            values[4] = x0 * test_value_y10;
-            values[5] = x1 * test_value_y10;
-            values[6] = x0 * test_value_y11;
-            values[7] = x1 * test_value_y11;
-          }
-      }
-  }
-
-
-
-  /**
-   * Calls the correct @p integrate_add_tensor_product_value_...()
-   * function depending on if values should be added to or set and if
-   * polynomials are linear.
-   */
-  template <bool is_linear, int dim, typename Number, typename Number2>
-  inline void
-  integrate_tensor_product_value(const dealii::ndarray<Number, 2, dim> *shapes,
-                                 const unsigned int        n_shapes,
-                                 const Number2            &value,
-                                 Number2                  *values,
-                                 const Point<dim, Number> &p,
-                                 const bool                do_add)
-  {
-    if (do_add)
-      {
-        if (is_linear)
-          internal::integrate_add_tensor_product_value_linear<dim,
-                                                              Number,
-                                                              Number2,
-                                                              true>(value,
-                                                                    values,
-                                                                    p);
-        else
-          internal::integrate_add_tensor_product_value_shapes<dim,
-                                                              Number,
-                                                              Number2,
-                                                              true>(shapes,
-                                                                    n_shapes,
-                                                                    value,
-                                                                    values);
-      }
-    else
-      {
-        if (is_linear)
-          internal::integrate_add_tensor_product_value_linear<dim,
-                                                              Number,
-                                                              Number2,
-                                                              false>(value,
-                                                                     values,
-                                                                     p);
-        else
-          internal::integrate_add_tensor_product_value_shapes<dim,
-                                                              Number,
-                                                              Number2,
-                                                              false>(shapes,
-                                                                     n_shapes,
-                                                                     value,
-                                                                     values);
-      }
-  }
-
-
-
   template <int dim, int n_points_1d_template, typename Number>
   inline void
   weight_fe_q_dofs_by_entity(const Number      *weights,
diff --git a/include/deal.II/matrix_free/tensor_product_point_kernels.h b/include/deal.II/matrix_free/tensor_product_point_kernels.h
new file mode 100644 (file)
index 0000000..4e4e76c
--- /dev/null
@@ -0,0 +1,1601 @@
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2020 - 2024 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+
+#ifndef dealii_matrix_free_tensor_product_point_kernels_h
+#define dealii_matrix_free_tensor_product_point_kernels_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/aligned_vector.h>
+#include <deal.II/base/ndarray.h>
+#include <deal.II/base/polynomial.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/matrix_free/shape_info.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+
+namespace internal
+{
+  /**
+   * Struct to avoid using Tensor<1, dim, Point<dim2>> in
+   * evaluate_tensor_product_value_and_gradient because a Point cannot be used
+   * within Tensor. Instead, a specialization of this struct upcasts the point
+   * to a Tensor<1,dim>.
+   */
+  template <typename Number, typename Number2>
+  struct ProductTypeNoPoint
+  {
+    using type = typename ProductType<Number, Number2>::type;
+  };
+
+  template <int dim, typename Number, typename Number2>
+  struct ProductTypeNoPoint<Point<dim, Number>, Number2>
+  {
+    using type = typename ProductType<Tensor<1, dim, Number>, Number2>::type;
+  };
+
+
+
+  /**
+   * Computes the values and derivatives of the 1d polynomials @p poly at the
+   * specified point @p p and stores it in @p shapes.
+   */
+  template <int dim, typename Number>
+  inline void
+  compute_values_of_array(
+    dealii::ndarray<Number, 2, dim>                    *shapes,
+    const std::vector<Polynomials::Polynomial<double>> &poly,
+    const Point<dim, Number>                           &p,
+    const unsigned int                                  derivative = 1)
+  {
+    const int n_shapes = poly.size();
+
+    // Evaluate 1d polynomials and their derivatives
+    std::array<Number, dim> point;
+    for (unsigned int d = 0; d < dim; ++d)
+      point[d] = p[d];
+    for (int i = 0; i < n_shapes; ++i)
+      poly[i].values_of_array(point, derivative, shapes[i].data());
+  }
+
+
+
+  /**
+   * Specialization of above function for dim = 0. Should not be called.
+   */
+  template <typename Number>
+  inline void
+  compute_values_of_array(dealii::ndarray<Number, 2, 0> *,
+                          const std::vector<Polynomials::Polynomial<double>> &,
+                          const Point<0, Number> &,
+                          const unsigned int)
+  {
+    DEAL_II_ASSERT_UNREACHABLE();
+  }
+
+
+
+  /**
+   * Interpolate inner dimensions of tensor product shape functions.
+   */
+  template <int dim,
+            int length,
+            typename Number2,
+            typename Number,
+            int  n_values    = 1,
+            bool do_renumber = true,
+            int  stride      = 1>
+  inline
+#ifndef DEBUG
+    DEAL_II_ALWAYS_INLINE
+#endif
+      std::array<typename ProductTypeNoPoint<Number, Number2>::type,
+                 2 + n_values>
+      do_interpolate_xy(const Number                           *values,
+                        const std::vector<unsigned int>        &renumber,
+                        const dealii::ndarray<Number2, 2, dim> *shapes,
+                        const int n_shapes_runtime,
+                        int      &i)
+  {
+    static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented");
+    static_assert(1 <= n_values && n_values <= 2,
+                  "Only n_values=1,2 implemented");
+
+    const int n_shapes = length > 0 ? length : n_shapes_runtime;
+
+    // If n_values > 1, we want to interpolate from a second array,
+    // placed in the same array immediately after the main data. This
+    // is used to interpolate normal derivatives onto faces.
+    const Number *values_2 =
+      n_values > 1 ?
+        values + stride * (length > 0 ?
+                             Utilities::pow(length, dim) :
+                             Utilities::fixed_power<dim>(n_shapes_runtime)) :
+        nullptr;
+    using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
+    std::array<Number3, 2 + n_values> result = {};
+    for (int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1)
+      {
+        // Interpolation + derivative x direction
+        std::array<Number3, 1 + n_values> inner_result = {};
+
+        // Distinguish the inner loop based on whether we have a
+        // renumbering or not
+        if (do_renumber && !renumber.empty())
+          for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
+            {
+              // gradient
+              inner_result[0] +=
+                shapes[i0][1][0] * values[renumber[i] * stride];
+              // values
+              inner_result[1] +=
+                shapes[i0][0][0] * values[renumber[i] * stride];
+              if (n_values > 1)
+                inner_result[2] +=
+                  shapes[i0][0][0] * values_2[renumber[i] * stride];
+            }
+        else
+          for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
+            {
+              // gradient
+              inner_result[0] += shapes[i0][1][0] * values[i * stride];
+              // values
+              inner_result[1] += shapes[i0][0][0] * values[i * stride];
+              if (n_values > 1)
+                inner_result[2] += shapes[i0][0][0] * values_2[i * stride];
+            }
+
+        if (dim > 1)
+          {
+            // Interpolation + derivative in y direction
+            // gradient
+            result[0] += inner_result[0] * shapes[i1][0][1];
+            result[1] += inner_result[1] * shapes[i1][1][1];
+            // values
+            result[2] += inner_result[1] * shapes[i1][0][1];
+            if (n_values > 1)
+              result[3] += inner_result[2] * shapes[i1][0][1];
+          }
+        else
+          {
+            // gradient
+            result[0] = inner_result[0];
+            // values
+            result[1] = inner_result[1];
+            if (n_values > 1)
+              result[2] = inner_result[2];
+          }
+      }
+    return result;
+  }
+
+
+
+  /**
+   * Interpolates the values and gradients into the points specified in
+   * @p compute_values_of_array() with help of the precomputed @p shapes.
+   */
+  template <int dim,
+            typename Number,
+            typename Number2,
+            int  n_values    = 1,
+            bool do_renumber = true,
+            int  stride      = 1>
+  inline std::array<typename ProductTypeNoPoint<Number, Number2>::type,
+                    dim + n_values>
+  evaluate_tensor_product_value_and_gradient_shapes(
+    const dealii::ndarray<Number2, 2, dim> *shapes,
+    const int                               n_shapes,
+    const Number                           *values,
+    const std::vector<unsigned int>        &renumber = {})
+  {
+    static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented");
+    static_assert(1 <= n_values && n_values <= 2,
+                  "Only n_values=1,2 implemented");
+
+    using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
+
+    std::array<Number3, dim + n_values> result = {};
+    if (dim == 0)
+      {
+        // We only need the interpolation of the value and normal derivatives on
+        // faces of a 1d element. As the interpolation is the value at the
+        // point, simply set the result vector accordingly.
+        result[0] = values[0];
+        if (n_values > 1)
+          result[1] = values[1 * stride];
+        return result;
+      }
+
+    // Go through the tensor product of shape functions and interpolate
+    // with optimal algorithm
+    for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
+      {
+        std::array<Number3, 2 + n_values> inner_result;
+        // Generate separate code with known loop bounds for the most common
+        // cases
+        if (n_shapes == 2)
+          inner_result =
+            do_interpolate_xy<dim,
+                              2,
+                              Number2,
+                              Number,
+                              n_values,
+                              do_renumber,
+                              stride>(values, renumber, shapes, n_shapes, i);
+        else if (n_shapes == 3)
+          inner_result =
+            do_interpolate_xy<dim,
+                              3,
+                              Number2,
+                              Number,
+                              n_values,
+                              do_renumber,
+                              stride>(values, renumber, shapes, n_shapes, i);
+        else if (n_shapes == 4)
+          inner_result =
+            do_interpolate_xy<dim,
+                              4,
+                              Number2,
+                              Number,
+                              n_values,
+                              do_renumber,
+                              stride>(values, renumber, shapes, n_shapes, i);
+        else if (n_shapes == 5)
+          inner_result =
+            do_interpolate_xy<dim,
+                              5,
+                              Number2,
+                              Number,
+                              n_values,
+                              do_renumber,
+                              stride>(values, renumber, shapes, n_shapes, i);
+        else if (n_shapes == 6)
+          inner_result =
+            do_interpolate_xy<dim,
+                              6,
+                              Number2,
+                              Number,
+                              n_values,
+                              do_renumber,
+                              stride>(values, renumber, shapes, n_shapes, i);
+        else
+          inner_result =
+            do_interpolate_xy<dim,
+                              -1,
+                              Number2,
+                              Number,
+                              n_values,
+                              do_renumber,
+                              stride>(values, renumber, shapes, n_shapes, i);
+        if (dim == 3)
+          {
+            // derivative + interpolation in z direction
+            // gradient
+            result[0] += inner_result[0] * shapes[i2][0][2];
+            result[1] += inner_result[1] * shapes[i2][0][2];
+            result[2] += inner_result[2] * shapes[i2][1][2];
+            // values
+            result[3] += inner_result[2] * shapes[i2][0][2];
+            if (n_values > 1)
+              result[4] += inner_result[3] * shapes[i2][0][2];
+          }
+        else if (dim == 2)
+          {
+            // gradient
+            result[0] = inner_result[0];
+            result[1] = inner_result[1];
+            // values
+            result[2] = inner_result[2];
+            if (n_values > 1)
+              result[3] = inner_result[3];
+          }
+        else
+          {
+            // gradient
+            result[0] = inner_result[0];
+            // values
+            result[1] = inner_result[1];
+            if (n_values > 1)
+              result[2] = inner_result[2];
+          }
+      }
+
+    return result;
+  }
+
+
+
+  /**
+   * Specializes @p evaluate_tensor_product_value_and_gradient() for linear
+   * polynomials which massively reduces the necessary instructions.
+   */
+  template <int dim,
+            typename Number,
+            typename Number2,
+            int n_values = 1,
+            int stride   = 1>
+  inline std::array<typename ProductTypeNoPoint<Number, Number2>::type,
+                    dim + n_values>
+  evaluate_tensor_product_value_and_gradient_linear(
+    const Number              *values,
+    const Point<dim, Number2> &p)
+  {
+    static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented");
+    static_assert(1 <= n_values && n_values <= 2,
+                  "Only n_values=1,2 implemented");
+
+    using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
+
+    // If n_values > 1, we want to interpolate from a second array,
+    // placed in the same array immediately after the main data. This
+    // is used to interpolate normal derivatives onto faces.
+
+    std::array<Number3, dim + n_values> result;
+    if (dim == 0)
+      {
+        // we only need the value on faces of a 1d element
+        result[0] = values[0];
+        if (n_values > 1)
+          result[1] = values[1 * stride];
+      }
+    else if (dim == 1)
+      {
+        // gradient
+        result[0] = Number3(values[stride] - values[0]);
+        // values
+        result[1] = Number3(values[0]) + p[0] * result[0];
+        if (n_values > 1)
+          result[2] = Number3(values[2 * stride]) +
+                      p[0] * (values[3 * stride] - values[2 * stride]);
+      }
+    else if (dim == 2)
+      {
+        const Number3 val10 = Number3(values[stride] - values[0]);
+        const Number3 val32 = Number3(values[3 * stride] - values[2 * stride]);
+        const Number3 tmp0  = Number3(values[0]) + p[0] * val10;
+        const Number3 tmp1  = Number3(values[2 * stride]) + p[0] * val32;
+
+        // gradient
+        result[0] = val10 + p[1] * (val32 - val10);
+        result[1] = tmp1 - tmp0;
+
+        // values
+        result[2] = tmp0 + p[1] * result[1];
+
+        if (n_values > 1)
+          {
+            const Number3 tmp0_2 =
+              Number3(values[4 * stride]) +
+              p[0] * (values[5 * stride] - values[4 * stride]);
+            const Number3 tmp1_2 =
+              Number3(values[6 * stride]) +
+              p[0] * (values[7 * stride] - values[6 * stride]);
+            result[3] = tmp0_2 + p[1] * (tmp1_2 - tmp0_2);
+          }
+      }
+    else if (dim == 3)
+      {
+        const Number3 val10 = Number3(values[stride] - values[0]);
+        const Number3 val32 = Number3(values[3 * stride] - values[2 * stride]);
+        const Number3 tmp0  = Number3(values[0]) + p[0] * val10;
+        const Number3 tmp1  = Number3(values[2 * stride]) + p[0] * val32;
+        const Number3 tmp10 = tmp1 - tmp0;
+        const Number3 tmpy0 = tmp0 + p[1] * tmp10;
+
+        const Number3 val54 = Number3(values[5 * stride] - values[4 * stride]);
+        const Number3 val76 = Number3(values[7 * stride] - values[6 * stride]);
+        const Number3 tmp2  = Number3(values[4 * stride]) + p[0] * val54;
+        const Number3 tmp3  = Number3(values[6 * stride]) + p[0] * val76;
+        const Number3 tmp32 = tmp3 - tmp2;
+        const Number3 tmpy1 = tmp2 + p[1] * tmp32;
+
+        // gradient
+        result[2]           = tmpy1 - tmpy0;
+        result[1]           = tmp10 + p[2] * (tmp32 - tmp10);
+        const Number3 tmpz0 = val10 + p[1] * (val32 - val10);
+        result[0] = tmpz0 + p[2] * (val54 + p[1] * (val76 - val54) - tmpz0);
+
+        // value
+        result[3] = tmpy0 + p[2] * result[2];
+        Assert(n_values == 1, ExcNotImplemented());
+      }
+
+    return result;
+  }
+
+
+
+  /**
+   * Compute the polynomial interpolation of a tensor product shape function
+   * $\varphi_i$ given a vector of coefficients $u_i$ in the form
+   * $u_h(\mathbf{x}) = \sum_{i=1}^{k^d} \varphi_i(\mathbf{x}) u_i$. The shape
+   * functions $\varphi_i(\mathbf{x}) =
+   * \prod_{d=1}^{\text{dim}}\varphi_{i_d}^\text{1d}(x_d)$ represent a tensor
+   * product. The function returns a pair with the value of the interpolation
+   * as the first component and the gradient in reference coordinates as the
+   * second component. Note that for compound types (e.g. the `values` field
+   * begin a Point<spacedim> argument), the components of the gradient are
+   * sorted as Tensor<1, dim, Tensor<1, spacedim>> with the derivatives
+   * as the first index; this is a consequence of the generic arguments in the
+   * function.
+   *
+   * @param poly The underlying one-dimensional polynomial basis
+   * $\{\varphi^{1d}_{i_1}\}$ given as a vector of polynomials.
+   *
+   * @param values The expansion coefficients $u_i$ of type `Number` in
+   * the polynomial interpolation. The coefficients can be simply `double`
+   * variables but e.g. also Point<spacedim> in case they define arithmetic
+   * operations with the type `Number2`.
+   *
+   * @param p The position in reference coordinates where the interpolation
+   * should be evaluated.
+   *
+   * @param d_linear Flag to specify whether a d-linear (linear in 1d,
+   * bi-linear in 2d, tri-linear in 3d) interpolation should be made, which
+   * allows to unroll loops and considerably speed up evaluation.
+   *
+   * @param renumber Optional parameter to specify a renumbering in the
+   * coefficient vector, assuming that `values[renumber[i]]` returns
+   * the lexicographic (tensor product) entry of the coefficients. If the
+   * vector is entry, the values are assumed to be sorted lexicographically.
+   */
+  template <int dim, typename Number, typename Number2>
+  inline std::pair<
+    typename ProductTypeNoPoint<Number, Number2>::type,
+    Tensor<1, dim, typename ProductTypeNoPoint<Number, Number2>::type>>
+  evaluate_tensor_product_value_and_gradient(
+    const std::vector<Polynomials::Polynomial<double>> &poly,
+    const std::vector<Number>                          &values,
+    const Point<dim, Number2>                          &p,
+    const bool                                          d_linear = false,
+    const std::vector<unsigned int>                    &renumber = {})
+  {
+    using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
+
+    std::array<Number3, dim + 1> result;
+    if (d_linear)
+      {
+        result =
+          evaluate_tensor_product_value_and_gradient_linear(values.data(), p);
+      }
+    else
+      {
+        AssertIndexRange(poly.size(), 200);
+        std::array<dealii::ndarray<Number2, 2, dim>, 200> shapes;
+        compute_values_of_array(shapes.data(), poly, p);
+        result = evaluate_tensor_product_value_and_gradient_shapes<dim,
+                                                                   Number,
+                                                                   Number2>(
+          shapes.data(), poly.size(), values.data(), renumber);
+      }
+    return std::make_pair(result[dim],
+                          Tensor<1, dim, Number3>(
+                            ArrayView<Number3>(result.data(), dim)));
+  }
+
+
+
+  template <int dim,
+            int length,
+            typename Number2,
+            typename Number,
+            bool do_renumber = true,
+            int  stride      = 1>
+  inline
+#ifndef DEBUG
+    DEAL_II_ALWAYS_INLINE
+#endif
+    typename ProductTypeNoPoint<Number, Number2>::type
+    do_interpolate_xy_value(const Number                           *values,
+                            const std::vector<unsigned int>        &renumber,
+                            const dealii::ndarray<Number2, 2, dim> *shapes,
+                            const int n_shapes_runtime,
+                            int      &i)
+  {
+    const int n_shapes = length > 0 ? length : n_shapes_runtime;
+    using Number3      = typename ProductTypeNoPoint<Number, Number2>::type;
+    Number3 result     = {};
+    for (int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1)
+      {
+        // Interpolation x direction
+        Number3 value = {};
+
+        // Distinguish the inner loop based on whether we have a
+        // renumbering or not
+        if (do_renumber && !renumber.empty())
+          for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
+            value += shapes[i0][0][0] * values[renumber[i] * stride];
+        else
+          for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
+            value += shapes[i0][0][0] * values[i * stride];
+
+        if (dim > 1)
+          result += value * shapes[i1][0][1];
+        else
+          result = value;
+      }
+    return result;
+  }
+
+
+
+  template <int dim,
+            typename Number,
+            typename Number2,
+            bool do_renumber = true,
+            int  stride      = 1>
+  inline typename ProductTypeNoPoint<Number, Number2>::type
+  evaluate_tensor_product_value_shapes(
+    const dealii::ndarray<Number2, 2, dim> *shapes,
+    const int                               n_shapes,
+    const Number                           *values,
+    const std::vector<unsigned int>        &renumber = {})
+  {
+    static_assert(dim >= 0 && dim <= 3, "Only dim=0,1,2,3 implemented");
+
+    // we only need the value on faces of a 1d element
+    if (dim == 0)
+      {
+        return values[0];
+      }
+
+    using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
+
+    // Go through the tensor product of shape functions and interpolate
+    // with optimal algorithm
+    Number3 result = {};
+    for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
+      {
+        Number3 inner_result;
+        // Generate separate code with known loop bounds for the most common
+        // cases
+        if (n_shapes == 2)
+          inner_result = do_interpolate_xy_value<dim,
+                                                 2,
+                                                 Number2,
+                                                 Number,
+                                                 do_renumber,
+                                                 stride>(
+            values, renumber, shapes, n_shapes, i);
+        else if (n_shapes == 3)
+          inner_result = do_interpolate_xy_value<dim,
+                                                 3,
+                                                 Number2,
+                                                 Number,
+                                                 do_renumber,
+                                                 stride>(
+            values, renumber, shapes, n_shapes, i);
+        else if (n_shapes == 4)
+          inner_result = do_interpolate_xy_value<dim,
+                                                 4,
+                                                 Number2,
+                                                 Number,
+                                                 do_renumber,
+                                                 stride>(
+            values, renumber, shapes, n_shapes, i);
+        else if (n_shapes == 5)
+          inner_result = do_interpolate_xy_value<dim,
+                                                 5,
+                                                 Number2,
+                                                 Number,
+                                                 do_renumber,
+                                                 stride>(
+            values, renumber, shapes, n_shapes, i);
+        else if (n_shapes == 6)
+          inner_result = do_interpolate_xy_value<dim,
+                                                 6,
+                                                 Number2,
+                                                 Number,
+                                                 do_renumber,
+                                                 stride>(
+            values, renumber, shapes, n_shapes, i);
+        else
+          inner_result = do_interpolate_xy_value<dim,
+                                                 -1,
+                                                 Number2,
+                                                 Number,
+                                                 do_renumber,
+                                                 stride>(
+            values, renumber, shapes, n_shapes, i);
+        if (dim == 3)
+          {
+            // Interpolation + derivative in z direction
+            result += inner_result * shapes[i2][0][2];
+          }
+        else
+          {
+            result = inner_result;
+          }
+      }
+
+    return result;
+  }
+
+
+
+  template <int dim, typename Number, typename Number2, int stride = 1>
+  inline typename ProductTypeNoPoint<Number, Number2>::type
+  evaluate_tensor_product_value_linear(const Number              *values,
+                                       const Point<dim, Number2> &p)
+  {
+    static_assert(dim >= 0 && dim <= 3, "Only dim=0,1,2,3 implemented");
+
+    using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
+
+    if (dim == 0)
+      {
+        // we only need the value on faces of a 1d element
+        return values[0];
+      }
+    else if (dim == 1)
+      {
+        return Number3(values[0]) + p[0] * Number3(values[stride] - values[0]);
+      }
+    else if (dim == 2)
+      {
+        const Number3 val10 = Number3(values[stride] - values[0]);
+        const Number3 val32 = Number3(values[3 * stride] - values[2 * stride]);
+        const Number3 tmp0  = Number3(values[0]) + p[0] * val10;
+        const Number3 tmp1  = Number3(values[2 * stride]) + p[0] * val32;
+        return tmp0 + p[1] * (tmp1 - tmp0);
+      }
+    else if (dim == 3)
+      {
+        const Number3 val10 = Number3(values[stride] - values[0]);
+        const Number3 val32 = Number3(values[3 * stride] - values[2 * stride]);
+        const Number3 tmp0  = Number3(values[0]) + p[0] * val10;
+        const Number3 tmp1  = Number3(values[2 * stride]) + p[0] * val32;
+        const Number3 tmpy0 = tmp0 + p[1] * (tmp1 - tmp0);
+
+        const Number3 val54 = Number3(values[5 * stride] - values[4 * stride]);
+        const Number3 val76 = Number3(values[7 * stride] - values[6 * stride]);
+        const Number3 tmp2  = Number3(values[4 * stride]) + p[0] * val54;
+        const Number3 tmp3  = Number3(values[6 * stride]) + p[0] * val76;
+        const Number3 tmpy1 = tmp2 + p[1] * (tmp3 - tmp2);
+
+        return tmpy0 + p[2] * (tmpy1 - tmpy0);
+      }
+
+    // work around a compile error: missing return statement
+    return Number3();
+  }
+
+
+
+  template <int dim, typename Number, typename Number2>
+  inline typename ProductTypeNoPoint<Number, Number2>::type
+  evaluate_tensor_product_value(
+    const std::vector<Polynomials::Polynomial<double>> &poly,
+    const std::vector<Number>                          &values,
+    const Point<dim, Number2>                          &p,
+    const bool                                          d_linear = false,
+    const std::vector<unsigned int>                    &renumber = {})
+  {
+    typename ProductTypeNoPoint<Number, Number2>::type result;
+    if (d_linear)
+      {
+        result = evaluate_tensor_product_value_linear(values.data(), p);
+      }
+    else
+      {
+        AssertIndexRange(poly.size(), 200);
+        std::array<dealii::ndarray<Number2, 2, dim>, 200> shapes;
+        const int                n_shapes = poly.size();
+        std::array<Number2, dim> point;
+        for (unsigned int d = 0; d < dim; ++d)
+          point[d] = p[d];
+        for (int i = 0; i < n_shapes; ++i)
+          poly[i].values_of_array(point, 0, shapes[i].data());
+        result = evaluate_tensor_product_value_shapes<dim, Number, Number2>(
+          shapes.data(), n_shapes, values.data(), renumber);
+      }
+    return result;
+  }
+
+
+
+  /**
+   * This function computes derivatives of arbitrary orders in 1d, returning a
+   * Tensor with the respective derivative
+   */
+  template <int derivative_order, typename Number, typename Number2>
+  inline Tensor<1, 1, typename ProductTypeNoPoint<Number, Number2>::type>
+  evaluate_tensor_product_higher_derivatives(
+    const std::vector<Polynomials::Polynomial<double>> &poly,
+    const std::vector<Number>                          &values,
+    const Point<1, Number2>                            &p,
+    const std::vector<unsigned int>                    &renumber = {})
+  {
+    using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
+
+    const int n_shapes = poly.size();
+    AssertDimension(n_shapes, values.size());
+    Assert(renumber.empty() || renumber.size() == values.size(),
+           ExcDimensionMismatch(renumber.size(), values.size()));
+
+    std::array<Number2, derivative_order + 1> shapes;
+    Tensor<1, 1, Number3>                     result;
+    if (renumber.empty())
+      for (int i = 0; i < n_shapes; ++i)
+        {
+          poly[i].value(p[0], derivative_order, shapes.data());
+          result[0] += shapes[derivative_order] * values[i];
+        }
+    else
+      for (int i = 0; i < n_shapes; ++i)
+        {
+          poly[i].value(p[0], derivative_order, shapes.data());
+          result[0] += shapes[derivative_order] * values[renumber[i]];
+        }
+    return result;
+  }
+
+
+
+  /**
+   * This function computes derivatives of arbitrary orders in 2d, returning a
+   * Tensor with the respective derivatives
+   */
+  template <int derivative_order, typename Number, typename Number2>
+  inline Tensor<1,
+                derivative_order + 1,
+                typename ProductTypeNoPoint<Number, Number2>::type>
+  evaluate_tensor_product_higher_derivatives(
+    const std::vector<Polynomials::Polynomial<double>> &poly,
+    const std::vector<Number>                          &values,
+    const Point<2, Number2>                            &p,
+    const std::vector<unsigned int>                    &renumber = {})
+  {
+    using Number3     = typename ProductTypeNoPoint<Number, Number2>::type;
+    constexpr int dim = 2;
+
+    const int n_shapes = poly.size();
+    AssertDimension(Utilities::pow(n_shapes, 2), values.size());
+    Assert(renumber.empty() || renumber.size() == values.size(),
+           ExcDimensionMismatch(renumber.size(), values.size()));
+
+    AssertIndexRange(n_shapes, 100);
+    dealii::ndarray<Number2, 100, derivative_order + 1, dim> shapes;
+    // Evaluate 1d polynomials and their derivatives
+    std::array<Number2, dim> point;
+    for (unsigned int d = 0; d < dim; ++d)
+      point[d] = p[d];
+    for (int i = 0; i < n_shapes; ++i)
+      poly[i].values_of_array(point, derivative_order, &shapes[i][0]);
+
+    Tensor<1, derivative_order + 1, Number3> result;
+    for (int i1 = 0, i = 0; i1 < n_shapes; ++i1)
+      {
+        Tensor<1, derivative_order + 1, Number3> result_x;
+        if (renumber.empty())
+          for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
+            for (unsigned int d = 0; d <= derivative_order; ++d)
+              result_x[d] += shapes[i0][d][0] * values[i];
+        else
+          for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
+            for (unsigned int d = 0; d <= derivative_order; ++d)
+              result_x[d] += shapes[i0][d][0] * values[renumber[i]];
+
+        for (unsigned int d = 0; d <= derivative_order; ++d)
+          result[d] += shapes[i1][d][1] * result_x[derivative_order - d];
+      }
+    return result;
+  }
+
+
+
+  /**
+   * This function computes derivatives of arbitrary orders in 3d, returning a
+   * Tensor with the respective derivatives
+   */
+  template <int derivative_order, typename Number, typename Number2>
+  inline Tensor<1,
+                ((derivative_order + 1) * (derivative_order + 2)) / 2,
+                typename ProductTypeNoPoint<Number, Number2>::type>
+  evaluate_tensor_product_higher_derivatives(
+    const std::vector<Polynomials::Polynomial<double>> &poly,
+    const std::vector<Number>                          &values,
+    const Point<3, Number2>                            &p,
+    const std::vector<unsigned int>                    &renumber = {})
+  {
+    using Number3     = typename ProductTypeNoPoint<Number, Number2>::type;
+    constexpr int dim = 3;
+    constexpr int n_derivatives =
+      ((derivative_order + 1) * (derivative_order + 2)) / 2;
+
+    const int n_shapes = poly.size();
+    AssertDimension(Utilities::pow(n_shapes, 3), values.size());
+    Assert(renumber.empty() || renumber.size() == values.size(),
+           ExcDimensionMismatch(renumber.size(), values.size()));
+
+    AssertIndexRange(n_shapes, 100);
+    dealii::ndarray<Number2, 100, derivative_order + 1, dim> shapes;
+    // Evaluate 1d polynomials and their derivatives
+    std::array<Number2, dim> point;
+    for (unsigned int d = 0; d < dim; ++d)
+      point[d] = p[d];
+    for (int i = 0; i < n_shapes; ++i)
+      poly[i].values_of_array(point, derivative_order, &shapes[i][0]);
+
+    Tensor<1, n_derivatives, Number3> result;
+    for (int i2 = 0, i = 0; i2 < n_shapes; ++i2)
+      {
+        Tensor<1, n_derivatives, Number3> result_xy;
+        for (int i1 = 0; i1 < n_shapes; ++i1)
+          {
+            // apply x derivatives
+            Tensor<1, derivative_order + 1, Number3> result_x;
+            if (renumber.empty())
+              for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
+                for (unsigned int d = 0; d <= derivative_order; ++d)
+                  result_x[d] += shapes[i0][d][0] * values[i];
+            else
+              for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
+                for (unsigned int d = 0; d <= derivative_order; ++d)
+                  result_x[d] += shapes[i0][d][0] * values[renumber[i]];
+
+            // multiply by y derivatives, sorting them in upper triangular
+            // matrix, starting with highest global derivative order,
+            // decreasing the combined order of xy derivatives by one in each
+            // row, to be combined with z derivatives in the next step
+            for (unsigned int d = 0, c = 0; d <= derivative_order; ++d)
+              for (unsigned int e = d; e <= derivative_order; ++e, ++c)
+                result_xy[c] +=
+                  shapes[i1][e - d][1] * result_x[derivative_order - e];
+          }
+
+        // multiply by z derivatives, starting with highest x derivative
+        for (unsigned int d = 0, c = 0; d <= derivative_order; ++d)
+          for (unsigned int e = d; e <= derivative_order; ++e, ++c)
+            result[c] += shapes[i2][d][2] * result_xy[c];
+      }
+    return result;
+  }
+
+
+
+  template <int dim, typename Number, typename Number2>
+  SymmetricTensor<2, dim, typename ProductTypeNoPoint<Number, Number2>::type>
+  evaluate_tensor_product_hessian(
+    const std::vector<Polynomials::Polynomial<double>> &poly,
+    const std::vector<Number>                          &values,
+    const Point<dim, Number2>                          &p,
+    const std::vector<unsigned int>                    &renumber = {})
+  {
+    static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented");
+
+    const auto hessian =
+      evaluate_tensor_product_higher_derivatives<2>(poly, values, p, renumber);
+
+    using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
+    SymmetricTensor<2, dim, Number3> result;
+    if (dim == 1)
+      result[0][0] = hessian[0];
+    else if (dim >= 2)
+      {
+        // derivatives in Hessian are xx, xy, yy, xz, yz, zz, so must re-order
+        // them for 3D
+        for (unsigned int d = 0, c = 0; d < 2; ++d)
+          for (unsigned int e = d; e < 2; ++e, ++c)
+            result[d][e] = hessian[c];
+        if (dim == 3)
+          {
+            for (unsigned int d = 0; d < 2; ++d)
+              result[d][2] = hessian[3 + d];
+            result[2][2] = hessian[5];
+          }
+      }
+
+    return result;
+  }
+
+
+
+  /**
+   * Test inner dimensions of tensor product shape functions and accumulate.
+   */
+  template <int dim,
+            int length,
+            typename Number2,
+            typename Number,
+            bool add,
+            int  n_values = 1>
+  inline
+#ifndef DEBUG
+    DEAL_II_ALWAYS_INLINE
+#endif
+    void
+    do_apply_test_functions_xy(
+      Number2                                 *values,
+      const dealii::ndarray<Number, 2, dim>   *shapes,
+      const std::array<Number2, 2 + n_values> &test_grads_value,
+      const int                                n_shapes_runtime,
+      int                                     &i)
+  {
+    static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented");
+    static_assert(1 <= n_values && n_values <= 2,
+                  "Only n_values=1,2 implemented");
+
+    // Note that 'add' is a template argument, so the compiler will remove
+    // these checks
+    if (length > 0)
+      {
+        constexpr unsigned int         array_size = length > 0 ? length : 1;
+        std::array<Number, array_size> shape_values_x;
+        std::array<Number, array_size> shape_derivs_x;
+        for (unsigned int j = 0; j < array_size; ++j)
+          {
+            shape_values_x[j] = shapes[j][0][0];
+            shape_derivs_x[j] = shapes[j][1][0];
+          }
+        for (int i1 = 0; i1 < (dim > 1 ? length : 1); ++i1)
+          {
+            const Number2 test_value_y =
+              dim > 1 ? (test_grads_value[2] * shapes[i1][0][1] +
+                         test_grads_value[1] * shapes[i1][1][1]) :
+                        test_grads_value[2];
+            const Number2 test_grad_xy =
+              dim > 1 ? test_grads_value[0] * shapes[i1][0][1] :
+                        test_grads_value[0];
+            Number2 test_value_y_2;
+            if (n_values > 1)
+              test_value_y_2 = dim > 1 ?
+                                 test_grads_value[3] * shapes[i1][0][1] :
+                                 test_grads_value[3];
+
+            Number2 *values_ptr = values + i + i1 * length;
+            Number2 *values_ptr_2 =
+              n_values > 1 ? values_ptr + Utilities::pow(length, dim) : nullptr;
+            for (int i0 = 0; i0 < length; ++i0)
+              {
+                if (add)
+                  values_ptr[i0] += shape_values_x[i0] * test_value_y;
+                else
+                  values_ptr[i0] = shape_values_x[i0] * test_value_y;
+                values_ptr[i0] += shape_derivs_x[i0] * test_grad_xy;
+                if (n_values > 1)
+                  {
+                    if (add)
+                      values_ptr_2[i0] += shape_values_x[i0] * test_value_y_2;
+                    else
+                      values_ptr_2[i0] = shape_values_x[i0] * test_value_y_2;
+                  }
+              }
+          }
+        i += (dim > 1 ? length * length : length);
+      }
+    else
+      {
+        for (int i1 = 0; i1 < (dim > 1 ? n_shapes_runtime : 1); ++i1)
+          {
+            const Number2 test_value_y =
+              dim > 1 ? (test_grads_value[2] * shapes[i1][0][1] +
+                         test_grads_value[1] * shapes[i1][1][1]) :
+                        test_grads_value[2];
+            const Number2 test_grad_xy =
+              dim > 1 ? test_grads_value[0] * shapes[i1][0][1] :
+                        test_grads_value[0];
+            Number2 test_value_y_2;
+            if (n_values > 1)
+              test_value_y_2 = dim > 1 ?
+                                 test_grads_value[3] * shapes[i1][0][1] :
+                                 test_grads_value[3];
+
+            Number2 *values_ptr = values + i + i1 * n_shapes_runtime;
+            Number2 *values_ptr_2 =
+              n_values > 1 ?
+                values_ptr + Utilities::fixed_power<dim>(n_shapes_runtime) :
+                nullptr;
+            for (int i0 = 0; i0 < n_shapes_runtime; ++i0)
+              {
+                if (add)
+                  values_ptr[i0] += shapes[i0][0][0] * test_value_y;
+                else
+                  values_ptr[i0] = shapes[i0][0][0] * test_value_y;
+                values_ptr[i0] += shapes[i0][1][0] * test_grad_xy;
+                if (n_values > 1)
+                  {
+                    if (add)
+                      values_ptr_2[i0] += shapes[i0][0][0] * test_value_y_2;
+                    else
+                      values_ptr_2[i0] = shapes[i0][0][0] * test_value_y_2;
+                  }
+              }
+          }
+        i += (dim > 1 ? n_shapes_runtime * n_shapes_runtime : n_shapes_runtime);
+      }
+  }
+
+
+
+  /**
+   * Same as evaluate_tensor_product_value_and_gradient_shapes() but for
+   * integration.
+   */
+  template <int dim,
+            typename Number,
+            typename Number2,
+            bool add,
+            int  n_values = 1>
+  inline void
+  integrate_add_tensor_product_value_and_gradient_shapes(
+    const dealii::ndarray<Number, 2, dim> *shapes,
+    const int                              n_shapes,
+    const Number2                         *value,
+    const Tensor<1, dim, Number2>         &gradient,
+    Number2                               *values)
+  {
+    static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented");
+    static_assert(1 <= n_values && n_values <= 2,
+                  "Only n_values=1,2 implemented");
+
+    // Note that 'add' is a template argument, so the compiler will remove
+    // these checks
+    if (dim == 0)
+      {
+        if (add)
+          values[0] += value[0];
+        else
+          values[0] = value[0];
+        if (n_values > 1)
+          {
+            if (add)
+              values[1] += value[1];
+            else
+              values[1] = value[1];
+          }
+        return;
+      }
+
+    // Implement the transpose of the function above
+    // as in evaluate, use `int` type to produce better code in this context
+    std::array<Number2, 2 + n_values> test_grads_value;
+    for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
+      {
+        // test grad x
+        test_grads_value[0] =
+          dim > 2 ? gradient[0] * shapes[i2][0][2] : gradient[0];
+        // test grad y
+        test_grads_value[1] = dim > 2 ? gradient[1] * shapes[i2][0][2] :
+                                        (dim > 1 ? gradient[1] : Number2());
+        // test value z
+        test_grads_value[2] =
+          dim > 2 ?
+            (value[0] * shapes[i2][0][2] + gradient[2] * shapes[i2][1][2]) :
+            value[0];
+
+        if (n_values > 1)
+          test_grads_value[3] =
+            dim > 2 ? value[1] * shapes[i2][0][2] : value[1];
+        // Generate separate code with known loop bounds for the most common
+        // cases
+        if (n_shapes == 2)
+          do_apply_test_functions_xy<dim, 2, Number2, Number, add, n_values>(
+            values, shapes, test_grads_value, n_shapes, i);
+        else if (n_shapes == 3)
+          do_apply_test_functions_xy<dim, 3, Number2, Number, add, n_values>(
+            values, shapes, test_grads_value, n_shapes, i);
+        else if (n_shapes == 4)
+          do_apply_test_functions_xy<dim, 4, Number2, Number, add, n_values>(
+            values, shapes, test_grads_value, n_shapes, i);
+        else if (n_shapes == 5)
+          do_apply_test_functions_xy<dim, 5, Number2, Number, add, n_values>(
+            values, shapes, test_grads_value, n_shapes, i);
+        else if (n_shapes == 6)
+          do_apply_test_functions_xy<dim, 6, Number2, Number, add, n_values>(
+            values, shapes, test_grads_value, n_shapes, i);
+        else
+          do_apply_test_functions_xy<dim, -1, Number2, Number, add, n_values>(
+            values, shapes, test_grads_value, n_shapes, i);
+      }
+  }
+
+
+
+  /**
+   * Specializes @p integrate_add_tensor_product_value_and_gradient_shapes() for linear
+   * polynomials which massively reduces the necessary instructions.
+   */
+  template <int dim,
+            typename Number,
+            typename Number2,
+            bool add,
+            int  n_values = 1>
+  inline void
+  integrate_add_tensor_product_value_and_gradient_linear(
+    const Number2                 *value,
+    const Tensor<1, dim, Number2> &gradient,
+    Number2                       *values,
+    const Point<dim, Number>      &p)
+  {
+    static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented");
+    static_assert(1 <= n_values && n_values <= 2,
+                  "Only n_values=1,2 implemented");
+
+    // Note that 'add' is a template argument, so the compiler will remove
+    // these checks
+    if (dim == 0)
+      {
+        if (add)
+          values[0] += value[0];
+        else
+          values[0] = value[0];
+        if (n_values > 1)
+          {
+            if (add)
+              values[1] += value[1];
+            else
+              values[1] = value[1];
+          }
+      }
+    else if (dim == 1)
+      {
+        const Number2 difference = value[0] * p[0] + gradient[0];
+        if (add)
+          {
+            values[0] += value[0] - difference;
+            values[1] += difference;
+          }
+        else
+          {
+            values[0] = value[0] - difference;
+            values[1] = difference;
+          }
+        if (n_values > 1)
+          {
+            const Number2 product = value[1] * p[0];
+            if (add)
+              {
+                values[2] += value[1] - product;
+                values[3] += product;
+              }
+            else
+              {
+                values[2] = value[1] - product;
+                values[3] = product;
+              }
+          }
+      }
+    else if (dim == 2)
+      {
+        const Number2 test_value_y1 = value[0] * p[1] + gradient[1];
+        const Number2 test_value_y0 = value[0] - test_value_y1;
+        const Number2 test_grad_xy1 = gradient[0] * p[1];
+        const Number2 test_grad_xy0 = gradient[0] - test_grad_xy1;
+        const Number2 value0        = p[0] * test_value_y0 + test_grad_xy0;
+        const Number2 value1        = p[0] * test_value_y1 + test_grad_xy1;
+
+        if (add)
+          {
+            values[0] += test_value_y0 - value0;
+            values[1] += value0;
+            values[2] += test_value_y1 - value1;
+            values[3] += value1;
+          }
+        else
+          {
+            values[0] = test_value_y0 - value0;
+            values[1] = value0;
+            values[2] = test_value_y1 - value1;
+            values[3] = value1;
+          }
+
+        if (n_values > 1)
+          {
+            const Number2 test_value_y1_2 = value[1] * p[1];
+            const Number2 test_value_y0_2 = value[1] - test_value_y1_2;
+            const Number2 value0_2        = p[0] * test_value_y0_2;
+            const Number2 value1_2        = p[0] * test_value_y1_2;
+
+            if (add)
+              {
+                values[4] += test_value_y0_2 - value0_2;
+                values[5] += value0_2;
+                values[6] += test_value_y1_2 - value1_2;
+                values[7] += value1_2;
+              }
+            else
+              {
+                values[4] = test_value_y0_2 - value0_2;
+                values[5] = value0_2;
+                values[6] = test_value_y1_2 - value1_2;
+                values[7] = value1_2;
+              }
+          }
+      }
+    else if (dim == 3)
+      {
+        Assert(n_values == 1, ExcNotImplemented());
+
+        const Number2 test_value_z1 = value[0] * p[2] + gradient[2];
+        const Number2 test_value_z0 = value[0] - test_value_z1;
+        const Number2 test_grad_x1  = gradient[0] * p[2];
+        const Number2 test_grad_x0  = gradient[0] - test_grad_x1;
+        const Number2 test_grad_y1  = gradient[1] * p[2];
+        const Number2 test_grad_y0  = gradient[1] - test_grad_y1;
+
+        const Number2 test_value_y01 = test_value_z0 * p[1] + test_grad_y0;
+        const Number2 test_value_y00 = test_value_z0 - test_value_y01;
+        const Number2 test_grad_xy01 = test_grad_x0 * p[1];
+        const Number2 test_grad_xy00 = test_grad_x0 - test_grad_xy01;
+        const Number2 test_value_y11 = test_value_z1 * p[1] + test_grad_y1;
+        const Number2 test_value_y10 = test_value_z1 - test_value_y11;
+        const Number2 test_grad_xy11 = test_grad_x1 * p[1];
+        const Number2 test_grad_xy10 = test_grad_x1 - test_grad_xy11;
+
+        const Number2 value00 = p[0] * test_value_y00 + test_grad_xy00;
+        const Number2 value01 = p[0] * test_value_y01 + test_grad_xy01;
+        const Number2 value10 = p[0] * test_value_y10 + test_grad_xy10;
+        const Number2 value11 = p[0] * test_value_y11 + test_grad_xy11;
+
+        if (add)
+          {
+            values[0] += test_value_y00 - value00;
+            values[1] += value00;
+            values[2] += test_value_y01 - value01;
+            values[3] += value01;
+            values[4] += test_value_y10 - value10;
+            values[5] += value10;
+            values[6] += test_value_y11 - value11;
+            values[7] += value11;
+          }
+        else
+          {
+            values[0] = test_value_y00 - value00;
+            values[1] = value00;
+            values[2] = test_value_y01 - value01;
+            values[3] = value01;
+            values[4] = test_value_y10 - value10;
+            values[5] = value10;
+            values[6] = test_value_y11 - value11;
+            values[7] = value11;
+          }
+      }
+  }
+
+
+
+  /**
+   * Calls the correct @p integrate_add_tensor_product_value_and_gradient_...()
+   * function depending on if values should be added to or set and if
+   * polynomials are linear.
+   */
+  template <bool is_linear,
+            int  dim,
+            typename Number,
+            typename Number2,
+            int n_values = 1>
+  inline void
+  integrate_tensor_product_value_and_gradient(
+    const dealii::ndarray<Number, 2, dim> *shapes,
+    const unsigned int                     n_shapes,
+    const Number2                         *value,
+    const Tensor<1, dim, Number2>         &gradient,
+    Number2                               *values,
+    const Point<dim, Number>              &p,
+    const bool                             do_add)
+  {
+    if (do_add)
+      {
+        if (is_linear)
+          internal::integrate_add_tensor_product_value_and_gradient_linear<
+            dim,
+            Number,
+            Number2,
+            true,
+            n_values>(value, gradient, values, p);
+        else
+          internal::integrate_add_tensor_product_value_and_gradient_shapes<
+            dim,
+            Number,
+            Number2,
+            true,
+            n_values>(shapes, n_shapes, value, gradient, values);
+      }
+    else
+      {
+        if (is_linear)
+          internal::integrate_add_tensor_product_value_and_gradient_linear<
+            dim,
+            Number,
+            Number2,
+            false,
+            n_values>(value, gradient, values, p);
+        else
+          internal::integrate_add_tensor_product_value_and_gradient_shapes<
+            dim,
+            Number,
+            Number2,
+            false,
+            n_values>(shapes, n_shapes, value, gradient, values);
+      }
+  }
+
+
+
+  /**
+   * Test inner dimensions of tensor product shape functions and accumulate.
+   */
+  template <int dim, int length, typename Number2, typename Number, bool add>
+  inline
+#ifndef DEBUG
+    DEAL_II_ALWAYS_INLINE
+#endif
+    void
+    do_apply_test_functions_xy_value(
+      Number2                               *values,
+      const dealii::ndarray<Number, 2, dim> *shapes,
+      const Number2                         &test_value,
+      const int                              n_shapes_runtime,
+      int                                   &i)
+  {
+    if (length > 0)
+      {
+        constexpr unsigned int         array_size = length > 0 ? length : 1;
+        std::array<Number, array_size> shape_values_x;
+        for (unsigned int i1 = 0; i1 < array_size; ++i1)
+          shape_values_x[i1] = shapes[i1][0][0];
+        for (unsigned int i1 = 0; i1 < (dim > 1 ? length : 1); ++i1)
+          {
+            const Number2 test_value_y =
+              dim > 1 ? test_value * shapes[i1][0][1] : test_value;
+
+            Number2 *values_ptr = values + i + i1 * length;
+            for (unsigned int i0 = 0; i0 < length; ++i0)
+              {
+                if (add)
+                  values_ptr[i0] += shape_values_x[i0] * test_value_y;
+                else
+                  values_ptr[i0] = shape_values_x[i0] * test_value_y;
+              }
+          }
+        i += (dim > 1 ? length * length : length);
+      }
+    else
+      {
+        for (int i1 = 0; i1 < (dim > 1 ? n_shapes_runtime : 1); ++i1)
+          {
+            const Number2 test_value_y =
+              dim > 1 ? test_value * shapes[i1][0][1] : test_value;
+
+            Number2 *values_ptr = values + i + i1 * n_shapes_runtime;
+            for (int i0 = 0; i0 < n_shapes_runtime; ++i0)
+              {
+                if (add)
+                  values_ptr[i0] += shapes[i0][0][0] * test_value_y;
+                else
+                  values_ptr[i0] = shapes[i0][0][0] * test_value_y;
+              }
+          }
+        i += (dim > 1 ? n_shapes_runtime * n_shapes_runtime : n_shapes_runtime);
+      }
+  }
+
+
+
+  /**
+   * Same as evaluate_tensor_product_value_shapes() but for integration.
+   */
+  template <int dim, typename Number, typename Number2, bool add>
+  inline void
+  integrate_add_tensor_product_value_shapes(
+    const dealii::ndarray<Number, 2, dim> *shapes,
+    const int                              n_shapes,
+    const Number2                         &value,
+    Number2                               *values)
+  {
+    static_assert(dim >= 0 && dim <= 3, "Only dim=0,1,2,3 implemented");
+
+    // as in evaluate, use `int` type to produce better code in this context
+
+    if (dim == 0)
+      {
+        if (add)
+          values[0] += value;
+        else
+          values[0] = value;
+        return;
+      }
+
+    // Implement the transpose of the function above
+    Number2 test_value;
+    for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
+      {
+        // test value z
+        test_value = dim > 2 ? value * shapes[i2][0][2] : value;
+
+        // Generate separate code with known loop bounds for the most common
+        // cases
+        if (n_shapes == 2)
+          do_apply_test_functions_xy_value<dim, 2, Number2, Number, add>(
+            values, shapes, test_value, n_shapes, i);
+        else if (n_shapes == 3)
+          do_apply_test_functions_xy_value<dim, 3, Number2, Number, add>(
+            values, shapes, test_value, n_shapes, i);
+        else if (n_shapes == 4)
+          do_apply_test_functions_xy_value<dim, 4, Number2, Number, add>(
+            values, shapes, test_value, n_shapes, i);
+        else if (n_shapes == 5)
+          do_apply_test_functions_xy_value<dim, 5, Number2, Number, add>(
+            values, shapes, test_value, n_shapes, i);
+        else if (n_shapes == 6)
+          do_apply_test_functions_xy_value<dim, 6, Number2, Number, add>(
+            values, shapes, test_value, n_shapes, i);
+        else
+          do_apply_test_functions_xy_value<dim, -1, Number2, Number, add>(
+            values, shapes, test_value, n_shapes, i);
+      }
+  }
+
+
+
+  /**
+   * Specializes @p integrate_tensor_product_value_shapes() for linear
+   * polynomials which massively reduces the necessary instructions.
+   */
+  template <int dim, typename Number, typename Number2, bool add>
+  inline void
+  integrate_add_tensor_product_value_linear(const Number2            &value,
+                                            Number2                  *values,
+                                            const Point<dim, Number> &p)
+  {
+    static_assert(dim >= 0 && dim <= 3, "Only dim=0,1,2,3 implemented");
+
+    if (dim == 0)
+      {
+        if (add)
+          values[0] += value;
+        else
+          values[0] = value;
+      }
+    else if (dim == 1)
+      {
+        const auto x0 = 1. - p[0], x1 = p[0];
+
+        if (add)
+          {
+            values[0] += value * x0;
+            values[1] += value * x1;
+          }
+        else
+          {
+            values[0] = value * x0;
+            values[1] = value * x1;
+          }
+      }
+    else if (dim == 2)
+      {
+        const auto x0 = 1. - p[0], x1 = p[0], y0 = 1. - p[1], y1 = p[1];
+
+        const auto test_value_y0 = value * y0;
+        const auto test_value_y1 = value * y1;
+
+        if (add)
+          {
+            values[0] += x0 * test_value_y0;
+            values[1] += x1 * test_value_y0;
+            values[2] += x0 * test_value_y1;
+            values[3] += x1 * test_value_y1;
+          }
+        else
+          {
+            values[0] = x0 * test_value_y0;
+            values[1] = x1 * test_value_y0;
+            values[2] = x0 * test_value_y1;
+            values[3] = x1 * test_value_y1;
+          }
+      }
+    else if (dim == 3)
+      {
+        const auto x0 = 1. - p[0], x1 = p[0], y0 = 1. - p[1], y1 = p[1],
+                   z0 = 1. - p[2], z1 = p[2];
+
+        const auto test_value_z0 = value * z0;
+        const auto test_value_z1 = value * z1;
+
+        const auto test_value_y00 = test_value_z0 * y0;
+        const auto test_value_y01 = test_value_z0 * y1;
+        const auto test_value_y10 = test_value_z1 * y0;
+        const auto test_value_y11 = test_value_z1 * y1;
+
+        if (add)
+          {
+            values[0] += x0 * test_value_y00;
+            values[1] += x1 * test_value_y00;
+            values[2] += x0 * test_value_y01;
+            values[3] += x1 * test_value_y01;
+            values[4] += x0 * test_value_y10;
+            values[5] += x1 * test_value_y10;
+            values[6] += x0 * test_value_y11;
+            values[7] += x1 * test_value_y11;
+          }
+        else
+          {
+            values[0] = x0 * test_value_y00;
+            values[1] = x1 * test_value_y00;
+            values[2] = x0 * test_value_y01;
+            values[3] = x1 * test_value_y01;
+            values[4] = x0 * test_value_y10;
+            values[5] = x1 * test_value_y10;
+            values[6] = x0 * test_value_y11;
+            values[7] = x1 * test_value_y11;
+          }
+      }
+  }
+
+
+
+  /**
+   * Calls the correct @p integrate_add_tensor_product_value_...()
+   * function depending on if values should be added to or set and if
+   * polynomials are linear.
+   */
+  template <bool is_linear, int dim, typename Number, typename Number2>
+  inline void
+  integrate_tensor_product_value(const dealii::ndarray<Number, 2, dim> *shapes,
+                                 const unsigned int        n_shapes,
+                                 const Number2            &value,
+                                 Number2                  *values,
+                                 const Point<dim, Number> &p,
+                                 const bool                do_add)
+  {
+    if (do_add)
+      {
+        if (is_linear)
+          internal::integrate_add_tensor_product_value_linear<dim,
+                                                              Number,
+                                                              Number2,
+                                                              true>(value,
+                                                                    values,
+                                                                    p);
+        else
+          internal::integrate_add_tensor_product_value_shapes<dim,
+                                                              Number,
+                                                              Number2,
+                                                              true>(shapes,
+                                                                    n_shapes,
+                                                                    value,
+                                                                    values);
+      }
+    else
+      {
+        if (is_linear)
+          internal::integrate_add_tensor_product_value_linear<dim,
+                                                              Number,
+                                                              Number2,
+                                                              false>(value,
+                                                                     values,
+                                                                     p);
+        else
+          internal::integrate_add_tensor_product_value_shapes<dim,
+                                                              Number,
+                                                              Number2,
+                                                              false>(shapes,
+                                                                     n_shapes,
+                                                                     value,
+                                                                     values);
+      }
+  }
+} // end of namespace internal
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#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.