]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Restructure tensor-product kernels via 1d interpolants 15891/head
authorMartin Kronbichler <martin.kronbichler@uni-a.de>
Thu, 17 Aug 2023 21:11:16 +0000 (23:11 +0200)
committerMartin Kronbichler <martin.kronbichler@uni-a.de>
Fri, 18 Aug 2023 11:14:18 +0000 (13:14 +0200)
include/deal.II/matrix_free/tensor_product_kernels.h

index 98680df2e15f2f22070176071bfd2caad778b20c..5abcf6507e7428d05b062db8a0015f0055962f47 100644 (file)
@@ -95,1699 +95,961 @@ namespace internal
 
 
   /**
-   * Generic evaluator framework that valuates the given shape data in general
-   * dimensions using the tensor product form. Depending on the particular
-   * layout in the matrix entries, this corresponds to a usual matrix-matrix
-   * product or a matrix-matrix product including some symmetries.
-   *
-   * @tparam variant Variant of evaluation used for creating template
-   *                 specializations
-   * @tparam dim Dimension of the function
-   * @tparam n_rows Number of rows in the transformation matrix, which corresponds
-   *                to the number of 1d shape functions in the usual tensor
-   *                contraction setting
-   * @tparam n_columns Number of columns in the transformation matrix, which
-   *                   corresponds to the number of 1d shape functions in the
-   *                   usual tensor contraction setting
-   * @tparam Number Abstract number type for input and output arrays
-   * @tparam Number2 Abstract number type for coefficient arrays (defaults to
-   *                 same type as the input/output arrays); must implement
-   *                 operator* with Number to be valid
-   */
-  template <EvaluatorVariant variant,
-            int              dim,
-            int              n_rows,
-            int              n_columns,
-            typename Number,
-            typename Number2 = Number>
-  struct EvaluatorTensorProduct
-  {};
-
-  /**
-   * Evaluator framework for anisotropic polynomial spaces that valuates the
-   * given shape data in general dimensions using the tensor product form.
-   *
-   * @tparam variant Variant of evaluation used for creating template
-   *                 specializations
-   * @tparam dim Dimension of the function
-   * @tparam n_rows Number of rows in the transformation matrix, which corresponds
-   *                to the number of 1d shape functions in the usual tensor
-   *                contraction setting
-   * @tparam n_columns Number of columns in the transformation matrix, which
-   *                   corresponds to the number of 1d shape functions in the
-   *                   usual tensor contraction setting
-   * @tparam Number Abstract number type for input and output arrays
-   * @tparam Number2 Abstract number type for coefficient arrays (defaults to
-   *                 same type as the input/output arrays); must implement
-   *                 operator* with Number to be valid
-   * @tparam normal_dir Indicates the direction of the continuous component for the
-   *                    Raviart-Thomas space in terms of the normal onto the
-   * face, e.g 0 if the  is in x-direction, 1 if in y-direction, and 2 if in
-   * z-direction.
-   */
-  template <EvaluatorVariant variant,
-            int              dim,
-            int              n_rows,
-            int              n_columns,
-            int              normal_dir,
-            typename Number,
-            typename Number2 = Number>
-  struct EvaluatorTensorProductAnisotropic
-  {};
-
-
-
-  /**
-   * Internal evaluator for shape function in arbitrary dimension using the
-   * tensor product form of the basis functions.
-   *
-   * @tparam dim Space dimension in which this class is applied
-   * @tparam n_rows Number of rows in the transformation matrix, which corresponds
-   *                to the number of 1d shape functions in the usual tensor
-   *                contraction setting
-   * @tparam n_columns Number of columns in the transformation matrix, which
-   *                   corresponds to the number of 1d shape functions in the
-   *                   usual tensor contraction setting
-   * @tparam Number Abstract number type for input and output arrays
-   * @tparam Number2 Abstract number type for coefficient arrays (defaults to
-   *                 same type as the input/output arrays); must implement
-   *                 operator* with Number and produce Number as an output to
-   *                 be a valid type
-   */
-  template <int dim,
-            int n_rows,
-            int n_columns,
-            typename Number,
-            typename Number2>
-  struct EvaluatorTensorProduct<evaluate_general,
-                                dim,
-                                n_rows,
-                                n_columns,
-                                Number,
-                                Number2>
-  {
-    static constexpr unsigned int n_rows_of_product =
-      Utilities::pow(n_rows, dim);
-    static constexpr unsigned int n_columns_of_product =
-      Utilities::pow(n_columns, dim);
-
-    /**
-     * Empty constructor. Does nothing. Be careful when using 'values' and
-     * related methods because they need to be filled with the other pointer
-     */
-    EvaluatorTensorProduct()
-      : shape_values(nullptr)
-      , shape_gradients(nullptr)
-      , shape_hessians(nullptr)
-    {}
-
-    /**
-     * Constructor, taking the data from ShapeInfo
-     */
-    EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
-                           const AlignedVector<Number2> &shape_gradients,
-                           const AlignedVector<Number2> &shape_hessians,
-                           const unsigned int            dummy1 = 0,
-                           const unsigned int            dummy2 = 0)
-      : shape_values(shape_values.begin())
-      , shape_gradients(shape_gradients.begin())
-      , shape_hessians(shape_hessians.begin())
-    {
-      // We can enter this function either for the apply() path that has
-      // n_rows * n_columns entries or for the apply_face() path that only has
-      // n_rows * 3 entries in the array. Since we cannot decide about the use
-      // we must allow for both here.
-      Assert(shape_values.empty() ||
-               shape_values.size() == n_rows * n_columns ||
-               shape_values.size() == 3 * n_rows,
-             ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
-      Assert(shape_gradients.empty() ||
-               shape_gradients.size() == n_rows * n_columns,
-             ExcDimensionMismatch(shape_gradients.size(), n_rows * n_columns));
-      Assert(shape_hessians.empty() ||
-               shape_hessians.size() == n_rows * n_columns,
-             ExcDimensionMismatch(shape_hessians.size(), n_rows * n_columns));
-      (void)dummy1;
-      (void)dummy2;
-    }
-
-    /**
-     * Constructor, taking the data from ShapeInfo via raw pointers
-     */
-    EvaluatorTensorProduct(const Number2 *    shape_values,
-                           const Number2 *    shape_gradients,
-                           const Number2 *    shape_hessians,
-                           const unsigned int dummy1 = 0,
-                           const unsigned int dummy2 = 0)
-      : shape_values(shape_values)
-      , shape_gradients(shape_gradients)
-      , shape_hessians(shape_hessians)
-    {
-      (void)dummy1;
-      (void)dummy2;
-    }
-
-    template <int direction, bool contract_over_rows, bool add>
-    void
-    values(const Number in[], Number out[]) const
-    {
-      apply<direction, contract_over_rows, add>(shape_values, in, out);
-    }
-
-    template <int direction, bool contract_over_rows, bool add>
-    void
-    gradients(const Number in[], Number out[]) const
-    {
-      apply<direction, contract_over_rows, add>(shape_gradients, in, out);
-    }
-
-    template <int direction, bool contract_over_rows, bool add>
-    void
-    hessians(const Number in[], Number out[]) const
-    {
-      apply<direction, contract_over_rows, add>(shape_hessians, in, out);
-    }
-
-    template <int direction, bool contract_over_rows, bool add>
-    void
-    values_one_line(const Number in[], Number out[]) const
-    {
-      Assert(shape_values != nullptr, ExcNotInitialized());
-      apply<direction, contract_over_rows, add, true>(shape_values, in, out);
-    }
-
-    template <int direction, bool contract_over_rows, bool add>
-    void
-    gradients_one_line(const Number in[], Number out[]) const
-    {
-      Assert(shape_gradients != nullptr, ExcNotInitialized());
-      apply<direction, contract_over_rows, add, true>(shape_gradients, in, out);
-    }
-
-    template <int direction, bool contract_over_rows, bool add>
-    void
-    hessians_one_line(const Number in[], Number out[]) const
-    {
-      Assert(shape_hessians != nullptr, ExcNotInitialized());
-      apply<direction, contract_over_rows, add, true>(shape_hessians, in, out);
-    }
-
-    /**
-     * This function applies the tensor product kernel, corresponding to a
-     * multiplication of 1d stripes, along the given @p direction of the tensor
-     * data in the input array. This function allows the @p in and @p out
-     * arrays to alias for the case n_rows == n_columns, i.e., it is safe to
-     * perform the contraction in place where @p in and @p out point to the
-     * same address. For the case n_rows != n_columns, the output is in general
-     * not correct.
-     *
-     * @tparam direction Direction that is evaluated
-     * @tparam contract_over_rows If true, the tensor contraction sums
-     *                            over the rows in the given @p shape_data
-     *                            array, otherwise it sums over the columns
-     * @tparam add If true, the result is added to the output vector, else
-     *             the computed values overwrite the content in the output
-     * @tparam one_line If true, the kernel is only applied along a single 1d
-     *                  stripe within a dim-dimensional tensor, not the full
-     *                  n_rows^dim points as in the @p false case.
-     *
-     * @param shape_data Transformation matrix with @p n_rows rows and
-     *                   @p n_columns columns, stored in row-major format
-     * @param in Pointer to the start of the input data vector
-     * @param out Pointer to the start of the output data vector
-     */
-    template <int  direction,
-              bool contract_over_rows,
-              bool add,
-              bool one_line = false>
-    static void
-    apply(const Number2 *DEAL_II_RESTRICT shape_data,
-          const Number *                  in,
-          Number *                        out);
-
-    /**
-     * This function applies the tensor product operation to produce face values
-     * from cell values. As opposed to the apply method, this method assumes
-     * that the directions orthogonal to the face have n_rows degrees of
-     * freedom per direction and not n_columns for those directions lower than
-     * the one currently applied. In other words, apply_face() must be called
-     * before calling any interpolation within the face.
-     *
-     * @tparam face_direction Direction of the normal vector (0=x, 1=y, etc)
-     * @tparam contract_onto_face If true, the input vector is of size n_rows^dim
-     *                            and interpolation into n_rows^(dim-1) points
-     *                            is performed. This is a typical scenario in
-     *                            FEFaceEvaluation::evaluate() calls. If false,
-     *                            data from n_rows^(dim-1) points is expanded
-     *                            into the n_rows^dim points of the higher-
-     *                            dimensional data array. Derivatives in the
-     *                            case contract_onto_face==false are summed
-     *                            together
-     * @tparam add If true, the result is added to the output vector, else
-     *             the computed values overwrite the content in the output
-     * @tparam max_derivative Sets the number of derivatives that should be
-     *             computed. 0 means only values, 1 means values and first
-     *             derivatives, 2 second derivates. Note that all the
-     *             derivatives access the data in @p shape_values passed to
-     *             the constructor of the class
-     *
-     * @param in address of the input data vector
-     * @param out address of the output data vector
-     */
-    template <int  face_direction,
-              bool contract_onto_face,
-              bool add,
-              int  max_derivative>
-    void
-    apply_face(const Number *DEAL_II_RESTRICT in,
-               Number *DEAL_II_RESTRICT       out) const;
-
-  private:
-    const Number2 *shape_values;
-    const Number2 *shape_gradients;
-    const Number2 *shape_hessians;
-  };
-
-
-
-  template <int dim,
-            int n_rows,
-            int n_columns,
-            typename Number,
-            typename Number2>
-  template <int direction, bool contract_over_rows, bool add, bool one_line>
-  inline void
-  EvaluatorTensorProduct<evaluate_general,
-                         dim,
-                         n_rows,
-                         n_columns,
-                         Number,
-                         Number2>::apply(const Number2 *DEAL_II_RESTRICT
-                                                       shape_data,
-                                         const Number *in,
-                                         Number *      out)
-  {
-    static_assert(one_line == false || direction == dim - 1,
-                  "Single-line evaluation only works for direction=dim-1.");
-    Assert(shape_data != nullptr,
-           ExcMessage(
-             "The given array shape_data must not be the null pointer!"));
-    Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
-             in != out,
-           ExcMessage("In-place operation only supported for "
-                      "n_rows==n_columns or single-line interpolation"));
-    AssertIndexRange(direction, dim);
-    constexpr int mm = contract_over_rows ? n_rows : n_columns,
-                  nn = contract_over_rows ? n_columns : n_rows;
-
-    constexpr int stride    = Utilities::pow(n_columns, direction);
-    constexpr int n_blocks1 = one_line ? 1 : stride;
-    constexpr int n_blocks2 =
-      Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
-
-    for (int i2 = 0; i2 < n_blocks2; ++i2)
-      {
-        for (int i1 = 0; i1 < n_blocks1; ++i1)
-          {
-            Number x[mm];
-            for (int i = 0; i < mm; ++i)
-              x[i] = in[stride * i];
-            for (int col = 0; col < nn; ++col)
-              {
-                Number2 val0;
-                if (contract_over_rows == true)
-                  val0 = shape_data[col];
-                else
-                  val0 = shape_data[col * n_columns];
-                Number res0 = val0 * x[0];
-                for (int i = 1; i < mm; ++i)
-                  {
-                    if (contract_over_rows == true)
-                      val0 = shape_data[i * n_columns + col];
-                    else
-                      val0 = shape_data[col * n_columns + i];
-                    res0 += val0 * x[i];
-                  }
-                if (add)
-                  out[stride * col] += res0;
-                else
-                  out[stride * col] = res0;
-              }
-
-            if (one_line == false)
-              {
-                ++in;
-                ++out;
-              }
-          }
-        if (one_line == false)
-          {
-            in += stride * (mm - 1);
-            out += stride * (nn - 1);
-          }
-      }
-  }
-
-
-
-  template <int dim,
-            int n_rows,
-            int n_columns,
-            typename Number,
-            typename Number2>
-  template <int  face_direction,
-            bool contract_onto_face,
-            bool add,
-            int  max_derivative>
-  inline void
-  EvaluatorTensorProduct<evaluate_general,
-                         dim,
-                         n_rows,
-                         n_columns,
-                         Number,
-                         Number2>::apply_face(const Number *DEAL_II_RESTRICT in,
-                                              Number *DEAL_II_RESTRICT
-                                                out) const
-  {
-    Assert(dim > 0, ExcMessage("Only dim=1,2,3 supported"));
-    static_assert(max_derivative >= 0 && max_derivative < 3,
-                  "Only derivative orders 0-2 implemented");
-    Assert(shape_values != nullptr,
-           ExcMessage(
-             "The given array shape_values must not be the null pointer."));
-
-    constexpr int n_blocks1 = (dim > 1 ? n_rows : 1);
-    constexpr int n_blocks2 = (dim > 2 ? n_rows : 1);
-
-    AssertIndexRange(face_direction, dim);
-    constexpr int in_stride  = Utilities::pow(n_rows, face_direction);
-    constexpr int out_stride = Utilities::pow(n_rows, dim - 1);
-    const Number2 *DEAL_II_RESTRICT shape_values = this->shape_values;
-
-    for (int i2 = 0; i2 < n_blocks2; ++i2)
-      {
-        for (int i1 = 0; i1 < n_blocks1; ++i1)
-          {
-            if (contract_onto_face == true)
-              {
-                Number res0 = shape_values[0] * in[0];
-                Number res1, res2;
-                if (max_derivative > 0)
-                  res1 = shape_values[n_rows] * in[0];
-                if (max_derivative > 1)
-                  res2 = shape_values[2 * n_rows] * in[0];
-                for (int ind = 1; ind < n_rows; ++ind)
-                  {
-                    res0 += shape_values[ind] * in[in_stride * ind];
-                    if (max_derivative > 0)
-                      res1 += shape_values[ind + n_rows] * in[in_stride * ind];
-                    if (max_derivative > 1)
-                      res2 +=
-                        shape_values[ind + 2 * n_rows] * in[in_stride * ind];
-                  }
-                if (add)
-                  {
-                    out[0] += res0;
-                    if (max_derivative > 0)
-                      out[out_stride] += res1;
-                    if (max_derivative > 1)
-                      out[2 * out_stride] += res2;
-                  }
-                else
-                  {
-                    out[0] = res0;
-                    if (max_derivative > 0)
-                      out[out_stride] = res1;
-                    if (max_derivative > 1)
-                      out[2 * out_stride] = res2;
-                  }
-              }
-            else
-              {
-                for (int col = 0; col < n_rows; ++col)
-                  {
-                    if (add)
-                      out[col * in_stride] += shape_values[col] * in[0];
-                    else
-                      out[col * in_stride] = shape_values[col] * in[0];
-                    if (max_derivative > 0)
-                      out[col * in_stride] +=
-                        shape_values[col + n_rows] * in[out_stride];
-                    if (max_derivative > 1)
-                      out[col * in_stride] +=
-                        shape_values[col + 2 * n_rows] * in[2 * out_stride];
-                  }
-              }
-
-            // increment: in regular case, just go to the next point in
-            // x-direction. If we are at the end of one chunk in x-dir, need
-            // to jump over to the next layer in z-direction
-            switch (face_direction)
-              {
-                case 0:
-                  in += contract_onto_face ? n_rows : 1;
-                  out += contract_onto_face ? 1 : n_rows;
-                  break;
-                case 1:
-                  ++in;
-                  ++out;
-                  // faces 2 and 3 in 3d use local coordinate system zx, which
-                  // is the other way around compared to the tensor
-                  // product. Need to take that into account.
-                  if (dim == 3)
-                    {
-                      if (contract_onto_face)
-                        out += n_rows - 1;
-                      else
-                        in += n_rows - 1;
-                    }
-                  break;
-                case 2:
-                  ++in;
-                  ++out;
-                  break;
-                default:
-                  Assert(false, ExcNotImplemented());
-              }
-          }
-
-        // adjust for local coordinate system zx
-        if (face_direction == 1 && dim == 3)
-          {
-            if (contract_onto_face)
-              {
-                in += n_rows * (n_rows - 1);
-                out -= n_rows * n_rows - 1;
-              }
-            else
-              {
-                out += n_rows * (n_rows - 1);
-                in -= n_rows * n_rows - 1;
-              }
-          }
-      }
-  }
-
-
-
-  /**
-   * Internal evaluator for shape function using the tensor product form
-   * of the basis functions. The same as the other templated class but
-   * without making use of template arguments and variable loop bounds
-   * instead.
-   *
-   * @tparam dim Space dimension in which this class is applied
-   * @tparam Number Abstract number type for input and output arrays
-   * @tparam Number2 Abstract number type for coefficient arrays (defaults to
-   *                 same type as the input/output arrays); must implement
-   *                 operator* with Number and produce Number as an output to
-   *                 be a valid type
-   */
-  template <int dim, typename Number, typename Number2>
-  struct EvaluatorTensorProduct<evaluate_general, dim, 0, 0, Number, Number2>
-  {
-    static constexpr unsigned int n_rows_of_product =
-      numbers::invalid_unsigned_int;
-    static constexpr unsigned int n_columns_of_product =
-      numbers::invalid_unsigned_int;
-
-    /**
-     * Empty constructor. Does nothing. Be careful when using 'values' and
-     * related methods because they need to be filled with the other constructor
-     */
-    EvaluatorTensorProduct()
-      : shape_values(nullptr)
-      , shape_gradients(nullptr)
-      , shape_hessians(nullptr)
-      , n_rows(numbers::invalid_unsigned_int)
-      , n_columns(numbers::invalid_unsigned_int)
-    {}
-
-    /**
-     * Constructor, taking the data from ShapeInfo
-     */
-    EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
-                           const AlignedVector<Number2> &shape_gradients,
-                           const AlignedVector<Number2> &shape_hessians,
-                           const unsigned int            n_rows,
-                           const unsigned int            n_columns)
-      : shape_values(shape_values.begin())
-      , shape_gradients(shape_gradients.begin())
-      , shape_hessians(shape_hessians.begin())
-      , n_rows(n_rows)
-      , n_columns(n_columns)
-    {
-      // We can enter this function either for the apply() path that has
-      // n_rows * n_columns entries or for the apply_face() path that only has
-      // n_rows * 3 entries in the array. Since we cannot decide about the use
-      // we must allow for both here.
-      Assert(shape_values.empty() ||
-               shape_values.size() == n_rows * n_columns ||
-               shape_values.size() == n_rows * 3,
-             ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
-      Assert(shape_gradients.empty() ||
-               shape_gradients.size() == n_rows * n_columns,
-             ExcDimensionMismatch(shape_gradients.size(), n_rows * n_columns));
-      Assert(shape_hessians.empty() ||
-               shape_hessians.size() == n_rows * n_columns,
-             ExcDimensionMismatch(shape_hessians.size(), n_rows * n_columns));
-    }
-
-    /**
-     * Constructor, taking the data from ShapeInfo
-     */
-    EvaluatorTensorProduct(const Number2 *    shape_values,
-                           const Number2 *    shape_gradients,
-                           const Number2 *    shape_hessians,
-                           const unsigned int n_rows,
-                           const unsigned int n_columns)
-      : shape_values(shape_values)
-      , shape_gradients(shape_gradients)
-      , shape_hessians(shape_hessians)
-      , n_rows(n_rows)
-      , n_columns(n_columns)
-    {}
-
-    template <int direction, bool contract_over_rows, bool add>
-    void
-    values(const Number *in, Number *out) const
-    {
-      apply<direction, contract_over_rows, add>(shape_values, in, out);
-    }
-
-    template <int direction, bool contract_over_rows, bool add>
-    void
-    gradients(const Number *in, Number *out) const
-    {
-      apply<direction, contract_over_rows, add>(shape_gradients, in, out);
-    }
-
-    template <int direction, bool contract_over_rows, bool add>
-    void
-    hessians(const Number *in, Number *out) const
-    {
-      apply<direction, contract_over_rows, add>(shape_hessians, in, out);
-    }
-
-    template <int direction, bool contract_over_rows, bool add>
-    void
-    values_one_line(const Number in[], Number out[]) const
-    {
-      Assert(shape_values != nullptr, ExcNotInitialized());
-      apply<direction, contract_over_rows, add, true>(shape_values, in, out);
-    }
-
-    template <int direction, bool contract_over_rows, bool add>
-    void
-    gradients_one_line(const Number in[], Number out[]) const
-    {
-      Assert(shape_gradients != nullptr, ExcNotInitialized());
-      apply<direction, contract_over_rows, add, true>(shape_gradients, in, out);
-    }
-
-    template <int direction, bool contract_over_rows, bool add>
-    void
-    hessians_one_line(const Number in[], Number out[]) const
-    {
-      Assert(shape_hessians != nullptr, ExcNotInitialized());
-      apply<direction, contract_over_rows, add, true>(shape_hessians, in, out);
-    }
-
-    template <int  direction,
-              bool contract_over_rows,
-              bool add,
-              bool one_line = false>
-    void
-    apply(const Number2 *DEAL_II_RESTRICT shape_data,
-          const Number *                  in,
-          Number *                        out) const;
-
-    template <int  face_direction,
-              bool contract_onto_face,
-              bool add,
-              int  max_derivative>
-    void
-    apply_face(const Number *DEAL_II_RESTRICT in,
-               Number *DEAL_II_RESTRICT       out) const;
-
-    const Number2 *    shape_values;
-    const Number2 *    shape_gradients;
-    const Number2 *    shape_hessians;
-    const unsigned int n_rows;
-    const unsigned int n_columns;
-  };
-
-
-
-  template <int dim, typename Number, typename Number2>
-  template <int direction, bool contract_over_rows, bool add, bool one_line>
-  inline void
-  EvaluatorTensorProduct<evaluate_general, dim, 0, 0, Number, Number2>::apply(
-    const Number2 *DEAL_II_RESTRICT shape_data,
-    const Number *                  in,
-    Number *                        out) const
-  {
-    static_assert(one_line == false || direction == dim - 1,
-                  "Single-line evaluation only works for direction=dim-1.");
-    Assert(shape_data != nullptr,
-           ExcMessage(
-             "The given array shape_data must not be the null pointer!"));
-    Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
-             in != out,
-           ExcMessage("In-place operation only supported for "
-                      "n_rows==n_columns or single-line interpolation"));
-    AssertIndexRange(direction, dim);
-    const int mm = contract_over_rows ? n_rows : n_columns,
-              nn = contract_over_rows ? n_columns : n_rows;
-
-    const int stride =
-      direction == 0 ? 1 : Utilities::fixed_power<direction>(n_columns);
-    const int n_blocks1 = one_line ? 1 : stride;
-    const int n_blocks2 = direction >= dim - 1 ?
-                            1 :
-                            Utilities::fixed_power<dim - direction - 1>(n_rows);
-    Assert(n_rows <= 128, ExcNotImplemented());
-
-    // specialization for n_rows = 2 that manually unrolls the innermost loop
-    // to make the operation perform better (not completely as good as the
-    // templated one, but much better than the generic version down below,
-    // because the loop over col can be more effectively unrolled by the
-    // compiler)
-    if (contract_over_rows && n_rows == 2)
-      {
-        const Number2 *shape_data_1 = shape_data + n_columns;
-        for (int i2 = 0; i2 < n_blocks2; ++i2)
-          {
-            for (int i1 = 0; i1 < n_blocks1; ++i1)
-              {
-                const Number x0 = in[0], x1 = in[stride];
-                for (int col = 0; col < nn; ++col)
-                  {
-                    const Number result =
-                      shape_data[col] * x0 + shape_data_1[col] * x1;
-                    if (add)
-                      out[stride * col] += result;
-                    else
-                      out[stride * col] = result;
-                  }
-
-                if (one_line == false)
-                  {
-                    ++in;
-                    ++out;
-                  }
-              }
-            if (one_line == false)
-              {
-                in += stride * (mm - 1);
-                out += stride * (nn - 1);
-              }
-          }
-      }
-    // specialization for n = 3
-    else if (contract_over_rows && n_rows == 3)
-      {
-        const Number2 *shape_data_1 = shape_data + n_columns;
-        const Number2 *shape_data_2 = shape_data + 2 * n_columns;
-        for (int i2 = 0; i2 < n_blocks2; ++i2)
-          {
-            for (int i1 = 0; i1 < n_blocks1; ++i1)
-              {
-                const Number x0 = in[0], x1 = in[stride], x2 = in[2 * stride];
-                for (int col = 0; col < nn; ++col)
-                  {
-                    const Number result = shape_data[col] * x0 +
-                                          shape_data_1[col] * x1 +
-                                          shape_data_2[col] * x2;
-                    if (add)
-                      out[stride * col] += result;
-                    else
-                      out[stride * col] = result;
-                  }
-
-                if (one_line == false)
-                  {
-                    ++in;
-                    ++out;
-                  }
-              }
-            if (one_line == false)
-              {
-                in += stride * (mm - 1);
-                out += stride * (nn - 1);
-              }
-          }
-      }
-    // general loop for all other cases
-    else
-      for (int i2 = 0; i2 < n_blocks2; ++i2)
-        {
-          for (int i1 = 0; i1 < n_blocks1; ++i1)
-            {
-              Number x[129];
-              for (int i = 0; i < mm; ++i)
-                x[i] = in[stride * i];
-              for (int col = 0; col < nn; ++col)
-                {
-                  Number2 val0;
-                  if (contract_over_rows == true)
-                    val0 = shape_data[col];
-                  else
-                    val0 = shape_data[col * n_columns];
-                  Number res0 = val0 * x[0];
-                  for (int i = 1; i < mm; ++i)
-                    {
-                      if (contract_over_rows == true)
-                        val0 = shape_data[i * n_columns + col];
-                      else
-                        val0 = shape_data[col * n_columns + i];
-                      res0 += val0 * x[i];
-                    }
-                  if (add)
-                    out[stride * col] += res0;
-                  else
-                    out[stride * col] = res0;
-                }
-
-              if (one_line == false)
-                {
-                  ++in;
-                  ++out;
-                }
-            }
-          if (one_line == false)
-            {
-              in += stride * (mm - 1);
-              out += stride * (nn - 1);
-            }
-        }
-  }
-
-
-
-  template <int dim, typename Number, typename Number2>
-  template <int  face_direction,
-            bool contract_onto_face,
-            bool add,
-            int  max_derivative>
-  inline void
-  EvaluatorTensorProduct<evaluate_general, dim, 0, 0, Number, Number2>::
-    apply_face(const Number *DEAL_II_RESTRICT in,
-               Number *DEAL_II_RESTRICT       out) const
-  {
-    Assert(shape_values != nullptr,
-           ExcMessage(
-             "The given array shape_data must not be the null pointer!"));
-    static_assert(dim > 0 && dim < 4, "Only dim=1,2,3 supported");
-    const int n_blocks1 = dim > 1 ? n_rows : 1;
-    const int n_blocks2 = dim > 2 ? n_rows : 1;
-
-    AssertIndexRange(face_direction, dim);
-    const int in_stride =
-      face_direction > 0 ? Utilities::fixed_power<face_direction>(n_rows) : 1;
-    const int out_stride =
-      dim > 1 ? Utilities::fixed_power<dim - 1>(n_rows) : 1;
-
-    for (int i2 = 0; i2 < n_blocks2; ++i2)
-      {
-        for (int i1 = 0; i1 < n_blocks1; ++i1)
-          {
-            if (contract_onto_face == true)
-              {
-                Number res0 = shape_values[0] * in[0];
-                Number res1, res2;
-                if (max_derivative > 0)
-                  res1 = shape_values[n_rows] * in[0];
-                if (max_derivative > 1)
-                  res2 = shape_values[2 * n_rows] * in[0];
-                for (unsigned int ind = 1; ind < n_rows; ++ind)
-                  {
-                    res0 += shape_values[ind] * in[in_stride * ind];
-                    if (max_derivative > 0)
-                      res1 += shape_values[ind + n_rows] * in[in_stride * ind];
-                    if (max_derivative > 1)
-                      res2 +=
-                        shape_values[ind + 2 * n_rows] * in[in_stride * ind];
-                  }
-                if (add)
-                  {
-                    out[0] += res0;
-                    if (max_derivative > 0)
-                      out[out_stride] += res1;
-                    if (max_derivative > 1)
-                      out[2 * out_stride] += res2;
-                  }
-                else
-                  {
-                    out[0] = res0;
-                    if (max_derivative > 0)
-                      out[out_stride] = res1;
-                    if (max_derivative > 1)
-                      out[2 * out_stride] = res2;
-                  }
-              }
-            else
-              {
-                for (unsigned int col = 0; col < n_rows; ++col)
-                  {
-                    if (add)
-                      out[col * in_stride] += shape_values[col] * in[0];
-                    else
-                      out[col * in_stride] = shape_values[col] * in[0];
-                    if (max_derivative > 0)
-                      out[col * in_stride] +=
-                        shape_values[col + n_rows] * in[out_stride];
-                    if (max_derivative > 1)
-                      out[col * in_stride] +=
-                        shape_values[col + 2 * n_rows] * in[2 * out_stride];
-                  }
-              }
-
-            // increment: in regular case, just go to the next point in
-            // x-direction. If we are at the end of one chunk in x-dir, need
-            // to jump over to the next layer in z-direction
-            switch (face_direction)
-              {
-                case 0:
-                  in += contract_onto_face ? n_rows : 1;
-                  out += contract_onto_face ? 1 : n_rows;
-                  break;
-                case 1:
-                  ++in;
-                  ++out;
-                  // faces 2 and 3 in 3d use local coordinate system zx, which
-                  // is the other way around compared to the tensor
-                  // product. Need to take that into account.
-                  if (dim == 3)
-                    {
-                      if (contract_onto_face)
-                        out += n_rows - 1;
-                      else
-                        in += n_rows - 1;
-                    }
-                  break;
-                case 2:
-                  ++in;
-                  ++out;
-                  break;
-                default:
-                  Assert(false, ExcNotImplemented());
-              }
+   * One-dimensional kernel for use by the generic tensor product
+   * interpolation as provided by the class EvaluatorTensorProduct,
+   * implementing a matrix-vector product along this dimension, controlled by
+   * the number of rows and columns and the stride in the input and output
+   * arrays, which are embedded into some lexicographic ordering of unknowns
+   * in a tensor-product arrangement.
+   *
+   * Besides this generic function for templated loop lengths, there are
+   * several specializations of this class to account for run-time matrix
+   * sizes as well as some symmetries that reduce the data access or
+   * arithmetic operations. The specializations are technically realized by
+   * conditional function overloading with std::enable_if_t based on the first
+   * template parameter.
+   */
+  template <EvaluatorVariant  variant,
+            EvaluatorQuantity quantity,
+            int               n_rows,
+            int               n_columns,
+            int               stride_in,
+            int               stride_out,
+            bool              transpose_matrix,
+            bool              add,
+            typename Number,
+            typename Number2>
+  std::enable_if_t<(variant == evaluate_general), void>
+  apply_matrix_vector_product(const Number2 *matrix,
+                              const Number * in,
+                              Number *       out)
+  {
+    // We can only statically assert that one argument is non-zero because
+    // face evaluation might instantiate some functions, so we need to use the
+    // run-time assert to verify that we do not end up involuntarily.
+    static_assert(n_rows > 0 || n_columns > 0,
+                  "Specialization only for n_rows, n_columns > 0");
+    Assert(n_rows > 0 && n_columns > 0,
+           ExcInternalError("The evaluation needs n_rows, n_columns > 0, but " +
+                            std::to_string(n_rows) + ", " +
+                            std::to_string(n_columns) + " was passed!"));
+    static_assert(quantity == EvaluatorQuantity::value,
+                  "This function should only use EvaluatorQuantity::value");
+
+    constexpr int mm = transpose_matrix ? n_rows : n_columns,
+                  nn = transpose_matrix ? n_columns : n_rows;
+
+    std::array<Number, mm> x;
+    for (int i = 0; i < mm; ++i)
+      x[i] = in[stride_in * i];
+    for (int col = 0; col < nn; ++col)
+      {
+        Number res0;
+        if (transpose_matrix == true)
+          {
+            res0 = matrix[col] * x[0];
+            for (int i = 1; i < mm; ++i)
+              res0 += matrix[i * n_columns + col] * x[i];
           }
-        if (face_direction == 1 && dim == 3)
+        else
           {
-            // adjust for local coordinate system zx
-            if (contract_onto_face)
-              {
-                in += n_rows * (n_rows - 1);
-                out -= n_rows * n_rows - 1;
-              }
-            else
-              {
-                out += n_rows * (n_rows - 1);
-                in -= n_rows * n_rows - 1;
-              }
+            res0 = matrix[col * n_columns] * x[0];
+            for (int i = 1; i < mm; ++i)
+              res0 += matrix[col * n_columns + i] * x[i];
           }
+        if (add)
+          out[stride_out * col] += res0;
+        else
+          out[stride_out * col] = res0;
       }
   }
 
 
 
   /**
-   * Internal evaluator for 1d-3d shape function using the tensor product form
-   * of the basis functions. This class specializes the general application of
-   * tensor-product based elements for "symmetric" finite elements, i.e., when
-   * the shape functions are symmetric about 0.5 and the quadrature points
-   * are, too.
-   *
-   * @tparam dim Space dimension in which this class is applied
-   * @tparam n_rows Number of rows in the transformation matrix, which corresponds
-   *                to the number of 1d shape functions in the usual tensor
-   *                contraction setting
-   * @tparam n_columns Number of columns in the transformation matrix, which
-   *                   corresponds to the number of 1d shape functions in the
-   *                   usual tensor contraction setting
-   * @tparam Number Abstract number type for input and output arrays
-   * @tparam Number2 Abstract number type for coefficient arrays (defaults to
-   *                 same type as the input/output arrays); must implement
-   *                 operator* with Number and produce Number as an output to
-   *                 be a valid type
+   * Specialization of the matrix-vector kernel for run-time loop bounds in
+   * the generic evaluator.
    */
-  template <int dim,
-            int n_rows,
-            int n_columns,
+  template <EvaluatorVariant  variant,
+            EvaluatorQuantity quantity,
+            bool              transpose_matrix,
+            bool              add,
             typename Number,
             typename Number2>
-  struct EvaluatorTensorProduct<evaluate_symmetric,
-                                dim,
-                                n_rows,
-                                n_columns,
-                                Number,
-                                Number2>
+  std::enable_if_t<(variant == evaluate_general), void>
+  apply_matrix_vector_product(const Number2 *matrix,
+                              const Number * in,
+                              Number *       out,
+                              const int      n_rows,
+                              const int      n_columns,
+                              const int      stride_in,
+                              const int      stride_out)
   {
-    static constexpr unsigned int n_rows_of_product =
-      Utilities::pow(n_rows, dim);
-    static constexpr unsigned int n_columns_of_product =
-      Utilities::pow(n_columns, dim);
-
-    /**
-     * Constructor, taking the data from ShapeInfo
-     */
-    EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
-                           const AlignedVector<Number2> &shape_gradients,
-                           const AlignedVector<Number2> &shape_hessians,
-                           const unsigned int            dummy1 = 0,
-                           const unsigned int            dummy2 = 0)
-      : shape_values(shape_values.begin())
-      , shape_gradients(shape_gradients.begin())
-      , shape_hessians(shape_hessians.begin())
-    {
-      Assert(shape_values.empty() || shape_values.size() == n_rows * n_columns,
-             ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
-      Assert(shape_gradients.empty() ||
-               shape_gradients.size() == n_rows * n_columns,
-             ExcDimensionMismatch(shape_gradients.size(), n_rows * n_columns));
-      Assert(shape_hessians.empty() ||
-               shape_hessians.size() == n_rows * n_columns,
-             ExcDimensionMismatch(shape_hessians.size(), n_rows * n_columns));
-      (void)dummy1;
-      (void)dummy2;
-    }
-
-    template <int direction, bool contract_over_rows, bool add>
-    void
-    values(const Number in[], Number out[]) const;
+    const int mm = transpose_matrix ? n_rows : n_columns,
+              nn = transpose_matrix ? n_columns : n_rows;
+    Assert(n_rows <= 128, ExcNotImplemented());
+    Assert(n_rows > 0 && n_columns > 0,
+           ExcInternalError("Empty evaluation task!"));
+    Assert(n_rows > 0 && n_columns > 0,
+           ExcInternalError("The evaluation needs n_rows, n_columns > 0, but " +
+                            std::to_string(n_rows) + ", " +
+                            std::to_string(n_columns) + " was passed!"));
 
-    template <int direction, bool contract_over_rows, bool add>
-    void
-    gradients(const Number in[], Number out[]) const;
+    static_assert(quantity == EvaluatorQuantity::value,
+                  "This function should only use EvaluatorQuantity::value");
 
-    template <int direction, bool contract_over_rows, bool add>
-    void
-    hessians(const Number in[], Number out[]) const;
+    // specialization for n_rows = 2 that manually unrolls the innermost loop
+    // to make the operation perform better (not completely as good as the
+    // templated one, but much better than the generic version down below,
+    // because the loop over col can be more effectively unrolled by the
+    // compiler)
+    if (transpose_matrix && n_rows == 2)
+      {
+        const Number2 *matrix_1 = matrix + n_columns;
+        const Number   x0 = in[0], x1 = in[stride_in];
+        for (int col = 0; col < nn; ++col)
+          {
+            const Number result = matrix[col] * x0 + matrix_1[col] * x1;
+            if (add)
+              out[stride_out * col] += result;
+            else
+              out[stride_out * col] = result;
+          }
+      }
+    else if (transpose_matrix && n_rows == 3)
+      {
+        const Number2 *matrix_1 = matrix + n_columns;
+        const Number2 *matrix_2 = matrix_1 + n_columns;
+        const Number   x0 = in[0], x1 = in[stride_in], x2 = in[2 * stride_in];
+        for (int col = 0; col < nn; ++col)
+          {
+            const Number result =
+              matrix[col] * x0 + matrix_1[col] * x1 + matrix_2[col] * x2;
+            if (add)
+              out[stride_out * col] += result;
+            else
+              out[stride_out * col] = result;
+          }
+      }
+    else
+      {
+        std::array<Number, 129> x;
+        for (int i = 0; i < mm; ++i)
+          x[i] = in[stride_in * i];
 
-  private:
-    const Number2 *shape_values;
-    const Number2 *shape_gradients;
-    const Number2 *shape_hessians;
-  };
+        Number res0;
+        for (int col = 0; col < nn; ++col)
+          {
+            if (transpose_matrix == true)
+              {
+                res0 = matrix[col] * x[0];
+                for (int i = 1; i < mm; ++i)
+                  res0 += matrix[i * n_columns + col] * x[i];
+              }
+            else
+              {
+                res0 = matrix[col * n_columns] * x[0];
+                for (int i = 1; i < mm; ++i)
+                  res0 += matrix[col * n_columns + i] * x[i];
+              }
+            if (add)
+              out[stride_out * col] += res0;
+            else
+              out[stride_out * col] = res0;
+          }
+      }
+  }
 
 
 
-  // In this case, the 1d shape values read (sorted lexicographically, rows
-  // run over 1d dofs, columns over quadrature points):
-  // Q2 --> [ 0.687  0 -0.087 ]
-  //        [ 0.4    1  0.4   ]
-  //        [-0.087  0  0.687 ]
-  // Q3 --> [ 0.66   0.003  0.002  0.049 ]
-  //        [ 0.521  1.005 -0.01  -0.230 ]
-  //        [-0.230 -0.01   1.005  0.521 ]
-  //        [ 0.049  0.002  0.003  0.66  ]
-  // Q4 --> [ 0.658  0.022  0 -0.007 -0.032 ]
-  //        [ 0.608  1.059  0  0.039  0.176 ]
-  //        [-0.409 -0.113  1 -0.113 -0.409 ]
-  //        [ 0.176  0.039  0  1.059  0.608 ]
-  //        [-0.032 -0.007  0  0.022  0.658 ]
-  //
-  // In these matrices, we want to use avoid computations involving zeros and
-  // ones and in addition use the symmetry in entries to reduce the number of
-  // read operations.
-  template <int dim,
-            int n_rows,
-            int n_columns,
+  /**
+   * Internal evaluator specialized for "symmetric" finite elements, i.e.,
+   * when the shape functions and quadrature points are symmetric about the
+   * middle point, making the matrix entries the same when starting to read in
+   * the (1,1) entry forward compared to the (N,N) entry backward.
+   */
+  template <EvaluatorVariant  variant,
+            EvaluatorQuantity quantity,
+            int               n_rows,
+            int               n_columns,
+            int               stride_in,
+            int               stride_out,
+            bool              transpose_matrix,
+            bool              add,
             typename Number,
             typename Number2>
-  template <int direction, bool contract_over_rows, bool add>
-  inline void
-  EvaluatorTensorProduct<evaluate_symmetric,
-                         dim,
-                         n_rows,
-                         n_columns,
-                         Number,
-                         Number2>::values(const Number in[], Number out[]) const
+  std::enable_if_t<(variant == evaluate_symmetric), void>
+  apply_matrix_vector_product(const Number2 *matrix,
+                              const Number * in,
+                              Number *       out)
   {
-    Assert(shape_values != nullptr, ExcNotInitialized());
-    AssertIndexRange(direction, dim);
-    constexpr int mm     = contract_over_rows ? n_rows : n_columns,
-                  nn     = contract_over_rows ? n_columns : n_rows;
+    // We can only statically assert that one argument is non-zero because
+    // face evaluation might instantiate some functions, so we need to use the
+    // run-time assert to verify that we do not end up involuntarily.
+    static_assert(n_rows > 0 || n_columns > 0,
+                  "Specialization only for n_rows, n_columns > 0");
+    Assert(n_rows > 0 && n_columns > 0,
+           ExcInternalError("The evaluation needs n_rows, n_columns > 0, but " +
+                            std::to_string(n_rows) + ", " +
+                            std::to_string(n_columns) + " was passed!"));
+
+    constexpr int mm     = transpose_matrix ? n_rows : n_columns,
+                  nn     = transpose_matrix ? n_columns : n_rows;
     constexpr int n_cols = nn / 2;
     constexpr int mid    = mm / 2;
 
-    constexpr int stride    = Utilities::pow(n_columns, direction);
-    constexpr int n_blocks1 = stride;
-    constexpr int n_blocks2 =
-      Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
+    std::array<Number, mm> x;
+    for (int i = 0; i < mm; ++i)
+      x[i] = in[stride_in * i];
 
-    for (int i2 = 0; i2 < n_blocks2; ++i2)
+    if (quantity == EvaluatorQuantity::value)
       {
-        for (int i1 = 0; i1 < n_blocks1; ++i1)
+        // In this case, the 1d shape values read (sorted lexicographically,
+        // rows run over 1d dofs, columns over quadrature points):
+        // Q2 --> [ 0.687  0 -0.087 ]
+        //        [ 0.4    1  0.4   ]
+        //        [-0.087  0  0.687 ]
+        // Q3 --> [ 0.66   0.003  0.002  0.049 ]
+        //        [ 0.521  1.005 -0.01  -0.230 ]
+        //        [-0.230 -0.01   1.005  0.521 ]
+        //        [ 0.049  0.002  0.003  0.66  ]
+        // Q4 --> [ 0.658  0.022  0 -0.007 -0.032 ]
+        //        [ 0.608  1.059  0  0.039  0.176 ]
+        //        [-0.409 -0.113  1 -0.113 -0.409 ]
+        //        [ 0.176  0.039  0  1.059  0.608 ]
+        //        [-0.032 -0.007  0  0.022  0.658 ]
+        //
+        // In these matrices, we want to use avoid computations involving
+        // zeros and ones and use the symmetry in entries starting from (1,1)
+        // forward and (N,N) backward, respectively to reduce the number of
+        // read operations.
+        for (int col = 0; col < n_cols; ++col)
           {
-            for (int col = 0; col < n_cols; ++col)
+            Number2 val0, val1;
+            Number  res0, res1;
+            if (transpose_matrix == true)
               {
-                Number2 val0, val1;
-                Number  in0, in1, res0, res1;
-                if (contract_over_rows == true)
-                  {
-                    val0 = shape_values[col];
-                    val1 = shape_values[nn - 1 - col];
-                  }
-                else
-                  {
-                    val0 = shape_values[col * n_columns];
-                    val1 = shape_values[(col + 1) * n_columns - 1];
-                  }
-                if (mid > 0)
-                  {
-                    in0  = in[0];
-                    in1  = in[stride * (mm - 1)];
-                    res0 = val0 * in0;
-                    res1 = val1 * in0;
-                    res0 += val1 * in1;
-                    res1 += val0 * in1;
-                    for (int ind = 1; ind < mid; ++ind)
-                      {
-                        if (contract_over_rows == true)
-                          {
-                            val0 = shape_values[ind * n_columns + col];
-                            val1 = shape_values[ind * n_columns + nn - 1 - col];
-                          }
-                        else
-                          {
-                            val0 = shape_values[col * n_columns + ind];
-                            val1 =
-                              shape_values[(col + 1) * n_columns - 1 - ind];
-                          }
-                        in0 = in[stride * ind];
-                        in1 = in[stride * (mm - 1 - ind)];
-                        res0 += val0 * in0;
-                        res1 += val1 * in0;
-                        res0 += val1 * in1;
-                        res1 += val0 * in1;
-                      }
-                  }
-                else
-                  res0 = res1 = Number();
-                if (contract_over_rows == true)
+                val0 = matrix[col];
+                val1 = matrix[nn - 1 - col];
+              }
+            else
+              {
+                val0 = matrix[col * n_columns];
+                val1 = matrix[(col + 1) * n_columns - 1];
+              }
+            if (mid > 0)
+              {
+                res0 = val0 * x[0];
+                res1 = val1 * x[0];
+                res0 += val1 * x[mm - 1];
+                res1 += val0 * x[mm - 1];
+                for (int ind = 1; ind < mid; ++ind)
                   {
-                    if (mm % 2 == 1)
+                    if (transpose_matrix == true)
                       {
-                        val0 = shape_values[mid * n_columns + col];
-                        in1  = val0 * in[stride * mid];
-                        res0 += in1;
-                        res1 += in1;
+                        val0 = matrix[ind * n_columns + col];
+                        val1 = matrix[ind * n_columns + nn - 1 - col];
                       }
-                  }
-                else
-                  {
-                    if (mm % 2 == 1 && nn % 2 == 0)
+                    else
                       {
-                        val0 = shape_values[col * n_columns + mid];
-                        in1  = val0 * in[stride * mid];
-                        res0 += in1;
-                        res1 += in1;
+                        val0 = matrix[col * n_columns + ind];
+                        val1 = matrix[(col + 1) * n_columns - 1 - ind];
                       }
+                    res0 += val0 * x[ind];
+                    res1 += val1 * x[ind];
+                    res0 += val1 * x[mm - 1 - ind];
+                    res1 += val0 * x[mm - 1 - ind];
                   }
-                if (add)
+              }
+            else
+              res0 = res1 = Number();
+            if (transpose_matrix == true)
+              {
+                if (mm % 2 == 1)
                   {
-                    out[stride * col] += res0;
-                    out[stride * (nn - 1 - col)] += res1;
+                    const Number tmp = matrix[mid * n_columns + col] * x[mid];
+                    res0 += tmp;
+                    res1 += tmp;
                   }
-                else
+              }
+            else
+              {
+                if (mm % 2 == 1 && nn % 2 == 0)
                   {
-                    out[stride * col]            = res0;
-                    out[stride * (nn - 1 - col)] = res1;
+                    const Number tmp = matrix[col * n_columns + mid] * x[mid];
+                    res0 += tmp;
+                    res1 += tmp;
                   }
               }
-            if (contract_over_rows == true && nn % 2 == 1 && mm % 2 == 1)
+            if (add)
               {
-                if (add)
-                  out[stride * n_cols] += in[stride * mid];
-                else
-                  out[stride * n_cols] = in[stride * mid];
+                out[stride_out * col] += res0;
+                out[stride_out * (nn - 1 - col)] += res1;
+              }
+            else
+              {
+                out[stride_out * col]            = res0;
+                out[stride_out * (nn - 1 - col)] = res1;
               }
-            else if (contract_over_rows == true && nn % 2 == 1)
+          }
+        if (transpose_matrix == true && nn % 2 == 1 && mm % 2 == 1)
+          {
+            if (add)
+              out[stride_out * n_cols] += x[mid];
+            else
+              out[stride_out * n_cols] = x[mid];
+          }
+        else if (transpose_matrix == true && nn % 2 == 1)
+          {
+            Number res0;
+            if (mid > 0)
               {
-                Number  res0;
-                Number2 val0 = shape_values[n_cols];
-                if (mid > 0)
+                res0 = matrix[n_cols] * (x[0] + x[mm - 1]);
+                for (int ind = 1; ind < mid; ++ind)
                   {
-                    res0 = val0 * (in[0] + in[stride * (mm - 1)]);
-                    for (int ind = 1; ind < mid; ++ind)
-                      {
-                        val0 = shape_values[ind * n_columns + n_cols];
-                        res0 += val0 * (in[stride * ind] +
-                                        in[stride * (mm - 1 - ind)]);
-                      }
+                    const Number2 val0 = matrix[ind * n_columns + n_cols];
+                    res0 += val0 * (x[ind] + in[mm - 1 - ind]);
                   }
-                else
-                  res0 = Number();
-                if (add)
-                  out[stride * n_cols] += res0;
-                else
-                  out[stride * n_cols] = res0;
               }
-            else if (contract_over_rows == false && nn % 2 == 1)
+            else
+              res0 = Number();
+            if (add)
+              out[stride_out * n_cols] += res0;
+            else
+              out[stride_out * n_cols] = res0;
+          }
+        else if (transpose_matrix == false && nn % 2 == 1)
+          {
+            Number res0;
+            if (mid > 0)
+              {
+                res0 = matrix[n_cols * n_columns] * (x[0] + x[mm - 1]);
+                for (int ind = 1; ind < mid; ++ind)
+                  {
+                    const Number2 val0 = matrix[n_cols * n_columns + ind];
+                    res0 += val0 * (x[ind] + x[mm - 1 - ind]);
+                    ;
+                  }
+                if (mm % 2)
+                  res0 += x[mid];
+              }
+            else
+              res0 = in[0];
+            if (add)
+              out[stride_out * n_cols] += res0;
+            else
+              out[stride_out * n_cols] = res0;
+          }
+      }
+    else if (quantity == EvaluatorQuantity::gradient)
+      {
+        // For the specialized loop used for gradient computations we again
+        // exploit symmetries according to the following entries (sorted
+        // lexicographically, rows run over 1d dofs, columns over quadrature
+        // points):
+        // Q2 --> [-2.549 -1  0.549 ]
+        //        [ 3.098  0 -3.098 ]
+        //        [-0.549  1  2.549 ]
+        // Q3 --> [-4.315 -1.03  0.5  -0.44  ]
+        //        [ 6.07  -1.44 -2.97  2.196 ]
+        //        [-2.196  2.97  1.44 -6.07  ]
+        //        [ 0.44  -0.5   1.03  4.315 ]
+        // Q4 --> [-6.316 -1.3    0.333 -0.353  0.413 ]
+        //        [10.111 -2.76  -2.667  2.066 -2.306 ]
+        //        [-5.688  5.773  0     -5.773  5.688 ]
+        //        [ 2.306 -2.066  2.667  2.76 -10.111 ]
+        //        [-0.413  0.353 -0.333 -0.353  0.413 ]
+        for (int col = 0; col < n_cols; ++col)
+          {
+            Number2 val0, val1;
+            Number  res0, res1;
+            if (transpose_matrix == true)
+              {
+                val0 = matrix[col];
+                val1 = matrix[nn - 1 - col];
+              }
+            else
+              {
+                val0 = matrix[col * n_columns];
+                val1 = matrix[(nn - col - 1) * n_columns];
+              }
+            if (mid > 0)
               {
-                Number res0;
-                if (mid > 0)
+                res0 = val0 * x[0];
+                res1 = val1 * x[0];
+                res0 -= val1 * x[mm - 1];
+                res1 -= val0 * x[mm - 1];
+                for (int ind = 1; ind < mid; ++ind)
                   {
-                    Number2 val0 = shape_values[n_cols * n_columns];
-                    res0         = val0 * (in[0] + in[stride * (mm - 1)]);
-                    for (int ind = 1; ind < mid; ++ind)
+                    if (transpose_matrix == true)
+                      {
+                        val0 = matrix[ind * n_columns + col];
+                        val1 = matrix[ind * n_columns + nn - 1 - col];
+                      }
+                    else
                       {
-                        val0       = shape_values[n_cols * n_columns + ind];
-                        Number in1 = val0 * (in[stride * ind] +
-                                             in[stride * (mm - 1 - ind)]);
-                        res0 += in1;
+                        val0 = matrix[col * n_columns + ind];
+                        val1 = matrix[(nn - col - 1) * n_columns + ind];
                       }
-                    if (mm % 2)
-                      res0 += in[stride * mid];
+                    res0 += val0 * x[ind];
+                    res1 += val1 * x[ind];
+                    res0 -= val1 * x[mm - 1 - ind];
+                    res1 -= val0 * x[mm - 1 - ind];
                   }
+              }
+            else
+              res0 = res1 = Number();
+            if (mm % 2 == 1)
+              {
+                if (transpose_matrix == true)
+                  val0 = matrix[mid * n_columns + col];
                 else
-                  res0 = in[0];
-                if (add)
-                  out[stride * n_cols] += res0;
+                  val0 = matrix[col * n_columns + mid];
+                const Number tmp = val0 * x[mid];
+                res0 += tmp;
+                res1 -= tmp;
+              }
+            if (add)
+              {
+                out[stride_out * col] += res0;
+                out[stride_out * (nn - 1 - col)] += res1;
+              }
+            else
+              {
+                out[stride_out * col]            = res0;
+                out[stride_out * (nn - 1 - col)] = res1;
+              }
+          }
+        if (nn % 2 == 1)
+          {
+            Number2 val0;
+            Number  res0;
+            if (transpose_matrix == true)
+              val0 = matrix[n_cols];
+            else
+              val0 = matrix[n_cols * n_columns];
+            res0 = val0 * (x[0] - x[mm - 1]);
+            for (int ind = 1; ind < mid; ++ind)
+              {
+                if (transpose_matrix == true)
+                  val0 = matrix[ind * n_columns + n_cols];
                 else
-                  out[stride * n_cols] = res0;
+                  val0 = matrix[n_cols * n_columns + ind];
+                Number in1 = val0 * (x[ind] - x[mm - 1 - ind]);
+                res0 += in1;
               }
-
-            ++in;
-            ++out;
+            if (add)
+              out[stride_out * n_cols] += res0;
+            else
+              out[stride_out * n_cols] = res0;
           }
-        in += stride * (mm - 1);
-        out += stride * (nn - 1);
       }
-  }
-
-
-
-  // For the specialized loop used for the gradient computation in
-  // here, the 1d shape values read (sorted lexicographically, rows
-  // run over 1d dofs, columns over quadrature points):
-  // Q2 --> [-2.549 -1  0.549 ]
-  //        [ 3.098  0 -3.098 ]
-  //        [-0.549  1  2.549 ]
-  // Q3 --> [-4.315 -1.03  0.5  -0.44  ]
-  //        [ 6.07  -1.44 -2.97  2.196 ]
-  //        [-2.196  2.97  1.44 -6.07  ]
-  //        [ 0.44  -0.5   1.03  4.315 ]
-  // Q4 --> [-6.316 -1.3    0.333 -0.353  0.413 ]
-  //        [10.111 -2.76  -2.667  2.066 -2.306 ]
-  //        [-5.688  5.773  0     -5.773  5.688 ]
-  //        [ 2.306 -2.066  2.667  2.76 -10.111 ]
-  //        [-0.413  0.353 -0.333 -0.353  0.413 ]
-  //
-  // In these matrices, we want to use avoid computations involving
-  // zeros and ones and in addition use the symmetry in entries to
-  // reduce the number of read operations.
-  template <int dim,
-            int n_rows,
-            int n_columns,
-            typename Number,
-            typename Number2>
-  template <int direction, bool contract_over_rows, bool add>
-  inline void
-  EvaluatorTensorProduct<evaluate_symmetric,
-                         dim,
-                         n_rows,
-                         n_columns,
-                         Number,
-                         Number2>::gradients(const Number in[],
-                                             Number       out[]) const
-  {
-    Assert(shape_gradients != nullptr, ExcNotInitialized());
-    AssertIndexRange(direction, dim);
-    constexpr int mm     = contract_over_rows ? n_rows : n_columns,
-                  nn     = contract_over_rows ? n_columns : n_rows;
-    constexpr int n_cols = nn / 2;
-    constexpr int mid    = mm / 2;
-
-    constexpr int stride    = Utilities::pow(n_columns, direction);
-    constexpr int n_blocks1 = stride;
-    constexpr int n_blocks2 =
-      Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
-
-    for (int i2 = 0; i2 < n_blocks2; ++i2)
+    else
       {
-        for (int i1 = 0; i1 < n_blocks1; ++i1)
+        // Hessians are almost the same as values, apart from some missing '1'
+        // entries
+        for (int col = 0; col < n_cols; ++col)
           {
-            for (int col = 0; col < n_cols; ++col)
+            Number2 val0, val1;
+            Number  res0, res1;
+            if (transpose_matrix == true)
               {
-                Number2 val0, val1;
-                Number  in0, in1, res0, res1;
-                if (contract_over_rows == true)
-                  {
-                    val0 = shape_gradients[col];
-                    val1 = shape_gradients[nn - 1 - col];
-                  }
-                else
-                  {
-                    val0 = shape_gradients[col * n_columns];
-                    val1 = shape_gradients[(nn - col - 1) * n_columns];
-                  }
-                if (mid > 0)
+                val0 = matrix[col];
+                val1 = matrix[nn - 1 - col];
+              }
+            else
+              {
+                val0 = matrix[col * n_columns];
+                val1 = matrix[(col + 1) * n_columns - 1];
+              }
+            if (mid > 0)
+              {
+                res0 = val0 * x[0];
+                res1 = val1 * x[0];
+                res0 += val1 * x[mm - 1];
+                res1 += val0 * x[mm - 1];
+                for (int ind = 1; ind < mid; ++ind)
                   {
-                    in0  = in[0];
-                    in1  = in[stride * (mm - 1)];
-                    res0 = val0 * in0;
-                    res1 = val1 * in0;
-                    res0 -= val1 * in1;
-                    res1 -= val0 * in1;
-                    for (int ind = 1; ind < mid; ++ind)
+                    if (transpose_matrix == true)
                       {
-                        if (contract_over_rows == true)
-                          {
-                            val0 = shape_gradients[ind * n_columns + col];
-                            val1 =
-                              shape_gradients[ind * n_columns + nn - 1 - col];
-                          }
-                        else
-                          {
-                            val0 = shape_gradients[col * n_columns + ind];
-                            val1 =
-                              shape_gradients[(nn - col - 1) * n_columns + ind];
-                          }
-                        in0 = in[stride * ind];
-                        in1 = in[stride * (mm - 1 - ind)];
-                        res0 += val0 * in0;
-                        res1 += val1 * in0;
-                        res0 -= val1 * in1;
-                        res1 -= val0 * in1;
+                        val0 = matrix[ind * n_columns + col];
+                        val1 = matrix[ind * n_columns + nn - 1 - col];
                       }
-                  }
-                else
-                  res0 = res1 = Number();
-                if (mm % 2 == 1)
-                  {
-                    if (contract_over_rows == true)
-                      val0 = shape_gradients[mid * n_columns + col];
                     else
-                      val0 = shape_gradients[col * n_columns + mid];
-                    in1 = val0 * in[stride * mid];
-                    res0 += in1;
-                    res1 -= in1;
-                  }
-                if (add)
-                  {
-                    out[stride * col] += res0;
-                    out[stride * (nn - 1 - col)] += res1;
-                  }
-                else
-                  {
-                    out[stride * col]            = res0;
-                    out[stride * (nn - 1 - col)] = res1;
+                      {
+                        val0 = matrix[col * n_columns + ind];
+                        val1 = matrix[(col + 1) * n_columns - 1 - ind];
+                      }
+                    res0 += val0 * x[ind];
+                    res1 += val1 * x[ind];
+                    res0 += val1 * x[mm - 1 - ind];
+                    res1 += val0 * x[mm - 1 - ind];
                   }
               }
-            if (nn % 2 == 1)
+            else
+              res0 = res1 = Number();
+            if (mm % 2 == 1)
               {
-                Number2 val0;
-                Number  res0;
-                if (contract_over_rows == true)
-                  val0 = shape_gradients[n_cols];
+                if (transpose_matrix == true)
+                  val0 = matrix[mid * n_columns + col];
                 else
-                  val0 = shape_gradients[n_cols * n_columns];
-                res0 = val0 * (in[0] - in[stride * (mm - 1)]);
+                  val0 = matrix[col * n_columns + mid];
+                const Number tmp = val0 * x[mid];
+                res0 += tmp;
+                res1 += tmp;
+              }
+            if (add)
+              {
+                out[stride_out * col] += res0;
+                out[stride_out * (nn - 1 - col)] += res1;
+              }
+            else
+              {
+                out[stride_out * col]            = res0;
+                out[stride_out * (nn - 1 - col)] = res1;
+              }
+          }
+        if (nn % 2 == 1)
+          {
+            Number2 val0;
+            Number  res0;
+            if (transpose_matrix == true)
+              val0 = matrix[n_cols];
+            else
+              val0 = matrix[n_cols * n_columns];
+            if (mid > 0)
+              {
+                res0 = val0 * (x[0] + x[mm - 1]);
                 for (int ind = 1; ind < mid; ++ind)
                   {
-                    if (contract_over_rows == true)
-                      val0 = shape_gradients[ind * n_columns + n_cols];
+                    if (transpose_matrix == true)
+                      val0 = matrix[ind * n_columns + n_cols];
                     else
-                      val0 = shape_gradients[n_cols * n_columns + ind];
-                    Number in1 =
-                      val0 * (in[stride * ind] - in[stride * (mm - 1 - ind)]);
+                      val0 = matrix[n_cols * n_columns + ind];
+                    Number in1 = val0 * (x[ind] + x[mm - 1 - ind]);
                     res0 += in1;
                   }
-                if (add)
-                  out[stride * n_cols] += res0;
+              }
+            else
+              res0 = Number();
+            if (mm % 2 == 1)
+              {
+                if (transpose_matrix == true)
+                  val0 = matrix[mid * n_columns + n_cols];
                 else
-                  out[stride * n_cols] = res0;
+                  val0 = matrix[n_cols * n_columns + mid];
+                res0 += val0 * x[mid];
               }
-
-            ++in;
-            ++out;
+            if (add)
+              out[stride_out * n_cols] += res0;
+            else
+              out[stride_out * n_cols] = res0;
           }
-        in += stride * (mm - 1);
-        out += stride * (nn - 1);
       }
   }
 
 
 
-  // evaluates the given shape data in 1d-3d using the tensor product
-  // form assuming the symmetries of unit cell shape hessians for
-  // finite elements in FEEvaluation
-  template <int dim,
-            int n_rows,
-            int n_columns,
+  /**
+   * Internal evaluator specialized for "symmetric" finite elements in the
+   * evenodd matrix format.
+   *
+   * This function implements a different approach to the symmetric case for
+   * values, gradients, and Hessians as in the above matrices: It is possible
+   * to reduce the cost per dimension from N^2 to N^2/2, where N is the number
+   * of 1d dofs (there are only N^2/2 different entries in the shape matrix,
+   * so this is plausible). The approach is based on the idea of applying the
+   * operator on the even and odd part of the input vectors separately, given
+   * that the basis of shape functions evaluated at quadrature points is
+   * symmetric. This method is presented e.g. in the book "Implementing
+   * Spectral Methods for Partial Differential Equations" by David A. Kopriva,
+   * Springer, 2009, section 3.5.3 (Even-Odd-Decomposition). Even though the
+   * experiments in the book say that the method is not efficient for N<20, it
+   * is more efficient in the context where the loop bounds are compile-time
+   * constants (templates).
+   */
+  template <EvaluatorVariant  variant,
+            EvaluatorQuantity quantity,
+            int               n_rows_static,
+            int               n_columns_static,
+            int               stride_in_static,
+            int               stride_out_static,
+            bool              transpose_matrix,
+            bool              add,
             typename Number,
             typename Number2>
-  template <int direction, bool contract_over_rows, bool add>
-  inline void
-  EvaluatorTensorProduct<evaluate_symmetric,
-                         dim,
-                         n_rows,
-                         n_columns,
-                         Number,
-                         Number2>::hessians(const Number in[],
-                                            Number       out[]) const
+  std::enable_if_t<(variant == evaluate_evenodd), void>
+  apply_matrix_vector_product(const Number2 *matrix,
+                              const Number * in,
+                              Number *       out,
+                              int            n_rows_runtime     = 0,
+                              int            n_columns_runtime  = 0,
+                              int            stride_in_runtime  = 0,
+                              int            stride_out_runtime = 0)
   {
-    Assert(shape_hessians != nullptr, ExcNotInitialized());
-    AssertIndexRange(direction, dim);
-    constexpr int mm     = contract_over_rows ? n_rows : n_columns;
-    constexpr int nn     = contract_over_rows ? n_columns : n_rows;
-    constexpr int n_cols = nn / 2;
-    constexpr int mid    = mm / 2;
+    const int n_rows = n_rows_static == 0 ? n_rows_runtime : n_rows_static;
+    const int n_columns =
+      n_rows_static == 0 ? n_columns_runtime : n_columns_static;
+    const int stride_in =
+      n_rows_static == 0 ? stride_in_runtime : stride_in_static;
+    const int stride_out =
+      n_rows_static == 0 ? stride_out_runtime : stride_out_static;
+
+    Assert(n_rows > 0 && n_columns > 0,
+           ExcInternalError("The evaluation needs n_rows, n_columns > 0, but " +
+                            std::to_string(n_rows) + ", " +
+                            std::to_string(n_columns) + " was passed!"));
+
+    const int mm     = transpose_matrix ? n_rows : n_columns,
+              nn     = transpose_matrix ? n_columns : n_rows;
+    const int n_cols = nn / 2;
+    const int mid    = mm / 2;
+
+    constexpr int max_mid = 16; // for non-templated execution
+    constexpr int static_mid =
+      n_rows_static == 0 ? 1 : (transpose_matrix ? n_rows : n_columns) / 2;
+    const int offset = (n_columns + 1) / 2;
 
-    constexpr int stride    = Utilities::pow(n_columns, direction);
-    constexpr int n_blocks1 = stride;
-    constexpr int n_blocks2 =
-      Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
+    Assert((n_rows_static != 0 && n_columns_static != 0) || mid <= max_mid,
+           ExcNotImplemented());
 
-    for (int i2 = 0; i2 < n_blocks2; ++i2)
+    std::array<Number, (n_rows_static != 0 ? static_mid : max_mid)> xp, xm;
+    for (int i = 0; i < mid; ++i)
       {
-        for (int i1 = 0; i1 < n_blocks1; ++i1)
+        if (transpose_matrix == true && quantity == EvaluatorQuantity::gradient)
+          {
+            xp[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
+            xm[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
+          }
+        else
+          {
+            xp[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
+            xm[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
+          }
+      }
+    Number xmid = in[stride_in * mid];
+    for (int col = 0; col < n_cols; ++col)
+      {
+        Number r0, r1;
+        if (mid > 0)
           {
-            for (int col = 0; col < n_cols; ++col)
+            if (transpose_matrix == true)
               {
-                Number2 val0, val1;
-                Number  in0, in1, res0, res1;
-                if (contract_over_rows == true)
-                  {
-                    val0 = shape_hessians[col];
-                    val1 = shape_hessians[nn - 1 - col];
-                  }
-                else
-                  {
-                    val0 = shape_hessians[col * n_columns];
-                    val1 = shape_hessians[(col + 1) * n_columns - 1];
-                  }
-                if (mid > 0)
-                  {
-                    in0  = in[0];
-                    in1  = in[stride * (mm - 1)];
-                    res0 = val0 * in0;
-                    res1 = val1 * in0;
-                    res0 += val1 * in1;
-                    res1 += val0 * in1;
-                    for (int ind = 1; ind < mid; ++ind)
-                      {
-                        if (contract_over_rows == true)
-                          {
-                            val0 = shape_hessians[ind * n_columns + col];
-                            val1 =
-                              shape_hessians[ind * n_columns + nn - 1 - col];
-                          }
-                        else
-                          {
-                            val0 = shape_hessians[col * n_columns + ind];
-                            val1 =
-                              shape_hessians[(col + 1) * n_columns - 1 - ind];
-                          }
-                        in0 = in[stride * ind];
-                        in1 = in[stride * (mm - 1 - ind)];
-                        res0 += val0 * in0;
-                        res1 += val1 * in0;
-                        res0 += val1 * in1;
-                        res1 += val0 * in1;
-                      }
-                  }
-                else
-                  res0 = res1 = Number();
-                if (mm % 2 == 1)
-                  {
-                    if (contract_over_rows == true)
-                      val0 = shape_hessians[mid * n_columns + col];
-                    else
-                      val0 = shape_hessians[col * n_columns + mid];
-                    in1 = val0 * in[stride * mid];
-                    res0 += in1;
-                    res1 += in1;
-                  }
-                if (add)
-                  {
-                    out[stride * col] += res0;
-                    out[stride * (nn - 1 - col)] += res1;
-                  }
-                else
-                  {
-                    out[stride * col]            = res0;
-                    out[stride * (nn - 1 - col)] = res1;
-                  }
+                r0 = matrix[col] * xp[0];
+                r1 = matrix[(n_rows - 1) * offset + col] * xm[0];
               }
-            if (nn % 2 == 1)
+            else
               {
-                Number2 val0;
-                Number  res0;
-                if (contract_over_rows == true)
-                  val0 = shape_hessians[n_cols];
-                else
-                  val0 = shape_hessians[n_cols * n_columns];
-                if (mid > 0)
+                r0 = matrix[col * offset] * xp[0];
+                r1 = matrix[(n_rows - 1 - col) * offset] * xm[0];
+              }
+            for (int ind = 1; ind < mid; ++ind)
+              {
+                if (transpose_matrix == true)
                   {
-                    res0 = val0 * (in[0] + in[stride * (mm - 1)]);
-                    for (int ind = 1; ind < mid; ++ind)
-                      {
-                        if (contract_over_rows == true)
-                          val0 = shape_hessians[ind * n_columns + n_cols];
-                        else
-                          val0 = shape_hessians[n_cols * n_columns + ind];
-                        Number in1 = val0 * (in[stride * ind] +
-                                             in[stride * (mm - 1 - ind)]);
-                        res0 += in1;
-                      }
+                    r0 += matrix[ind * offset + col] * xp[ind];
+                    r1 += matrix[(n_rows - 1 - ind) * offset + col] * xm[ind];
                   }
                 else
-                  res0 = Number();
-                if (mm % 2 == 1)
                   {
-                    if (contract_over_rows == true)
-                      val0 = shape_hessians[mid * n_columns + n_cols];
-                    else
-                      val0 = shape_hessians[n_cols * n_columns + mid];
-                    res0 += val0 * in[stride * mid];
+                    r0 += matrix[col * offset + ind] * xp[ind];
+                    r1 += matrix[(n_rows - 1 - col) * offset + ind] * xm[ind];
                   }
-                if (add)
-                  out[stride * n_cols] += res0;
-                else
-                  out[stride * n_cols] = res0;
               }
+          }
+        else
+          r0 = r1 = Number();
+        if (mm % 2 == 1 && transpose_matrix == true)
+          {
+            if (quantity == EvaluatorQuantity::gradient)
+              r1 += matrix[mid * offset + col] * xmid;
+            else
+              r0 += matrix[mid * offset + col] * xmid;
+          }
+        else if (mm % 2 == 1 &&
+                 (nn % 2 == 0 || quantity != EvaluatorQuantity::value ||
+                  mm == 3))
+          r0 += matrix[col * offset + mid] * xmid;
+
+        if (add)
+          {
+            out[stride_out * col] += r0 + r1;
+            if (quantity == EvaluatorQuantity::gradient &&
+                transpose_matrix == false)
+              out[stride_out * (nn - 1 - col)] += r1 - r0;
+            else
+              out[stride_out * (nn - 1 - col)] += r0 - r1;
+          }
+        else
+          {
+            out[stride_out * col] = r0 + r1;
+            if (quantity == EvaluatorQuantity::gradient &&
+                transpose_matrix == false)
+              out[stride_out * (nn - 1 - col)] = r1 - r0;
+            else
+              out[stride_out * (nn - 1 - col)] = r0 - r1;
+          }
+      }
+    if (quantity == EvaluatorQuantity::value && transpose_matrix == true &&
+        nn % 2 == 1 && mm % 2 == 1 && mm > 3)
+      {
+        if (add)
+          out[stride_out * n_cols] += matrix[mid * offset + n_cols] * xmid;
+        else
+          out[stride_out * n_cols] = matrix[mid * offset + n_cols] * xmid;
+      }
+    else if (transpose_matrix == true && nn % 2 == 1)
+      {
+        Number r0;
+        if (mid > 0)
+          {
+            r0 = matrix[n_cols] * xp[0];
+            for (int ind = 1; ind < mid; ++ind)
+              r0 += matrix[ind * offset + n_cols] * xp[ind];
+          }
+        else
+          r0 = Number();
+        if (quantity != EvaluatorQuantity::gradient && mm % 2 == 1)
+          r0 += matrix[mid * offset + n_cols] * xmid;
 
-            ++in;
-            ++out;
+        if (add)
+          out[stride_out * n_cols] += r0;
+        else
+          out[stride_out * n_cols] = r0;
+      }
+    else if (transpose_matrix == false && nn % 2 == 1)
+      {
+        Number r0;
+        if (mid > 0)
+          {
+            if (quantity == EvaluatorQuantity::gradient)
+              {
+                r0 = matrix[n_cols * offset] * xm[0];
+                for (int ind = 1; ind < mid; ++ind)
+                  r0 += matrix[n_cols * offset + ind] * xm[ind];
+              }
+            else
+              {
+                r0 = matrix[n_cols * offset] * xp[0];
+                for (int ind = 1; ind < mid; ++ind)
+                  r0 += matrix[n_cols * offset + ind] * xp[ind];
+              }
           }
-        in += stride * (mm - 1);
-        out += stride * (nn - 1);
+        else
+          r0 = Number();
+
+        if (quantity != EvaluatorQuantity::gradient && mm % 2 == 1)
+          r0 += matrix[n_cols * offset + mid] * xmid;
+
+        if (add)
+          out[stride_out * n_cols] += r0;
+        else
+          out[stride_out * n_cols] = r0;
       }
   }
 
 
 
-  template <int dim,
-            int n_rows_static,
-            int n_columns_static,
+  /**
+   * Internal evaluator specialized for "symmetric" finite elements in the
+   * evenodd matrix format with run-time bounds.
+   */
+  template <EvaluatorVariant  variant,
+            EvaluatorQuantity quantity,
+            bool              transpose_matrix,
+            bool              add,
             typename Number,
-            typename Number2,
-            int  direction,
-            bool contract_over_rows,
-            bool add,
-            int  type,
-            bool one_line>
-  inline void
-  even_odd_apply(const int                       n_rows_in,
-                 const int                       n_columns_in,
-                 const Number2 *DEAL_II_RESTRICT shapes,
-                 const Number *                  in,
-                 Number *                        out)
+            typename Number2>
+  std::enable_if_t<(variant == evaluate_evenodd), void>
+  apply_matrix_vector_product(const Number2 *matrix,
+                              const Number * in,
+                              Number *       out,
+                              int            n_rows,
+                              int            n_columns,
+                              int            stride_in,
+                              int            stride_out)
   {
-    static_assert(type < 3, "Only three variants type=0,1,2 implemented");
-    static_assert(one_line == false || direction == dim - 1,
-                  "Single-line evaluation only works for direction=dim-1.");
-
-    const int n_rows = n_rows_static == -1 ? n_rows_in : n_rows_static;
-    const int n_columns =
-      n_columns_static == -1 ? n_columns_in : n_columns_static;
-
-    Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
-             in != out,
-           ExcMessage("In-place operation only supported for "
-                      "n_rows==n_columns or single-line interpolation"));
+    apply_matrix_vector_product<evaluate_evenodd,
+                                quantity,
+                                0,
+                                0,
+                                0,
+                                0,
+                                transpose_matrix,
+                                add>(
+      matrix, in, out, n_rows, n_columns, stride_in, stride_out);
+  }
 
-    // We cannot statically assert that direction is less than dim, so must do
-    // an additional dynamic check
-    AssertIndexRange(direction, dim);
 
-    const int     nn = contract_over_rows ? n_columns : n_rows;
-    const int     mm = contract_over_rows ? n_rows : n_columns;
-    constexpr int mm_static =
-      contract_over_rows ? n_rows_static : n_columns_static;
-    const int     n_cols     = nn / 2;
-    const int     mid        = mm / 2;
-    constexpr int mid_static = mm_static / 2;
-    constexpr int max_mid    = 15; // for non-templated execution
 
-    Assert((n_rows_static != -1 && n_columns_static != -1) || mid <= max_mid,
-           ExcNotImplemented());
+  /**
+   * Internal evaluator specialized for "symmetric" finite elements in the
+   * symmetric_hierarchical matrix format.
+   *
+   * This class implements an approach similar to the even-odd decomposition
+   * but with a different type of symmetry. In this case, we assume that a
+   * single shape function already shows the symmetry over the quadrature
+   * points, rather than the complete basis that is considered in the even-odd
+   * case. In particular, we assume that the shape functions are ordered as in
+   * the Legendre basis, with symmetric shape functions in the even slots
+   * (rows of the values array) and point-symmetric in the odd slots. Like the
+   * even-odd decomposition, the number of operations are N^2/2 rather than
+   * N^2 FMAs (fused multiply-add), where N is the number of 1d dofs. The
+   * difference is in the way the input and output quantities are symmetrized.
+   */
+  template <EvaluatorVariant  variant,
+            EvaluatorQuantity quantity,
+            int               n_rows,
+            int               n_columns,
+            int               stride_in,
+            int               stride_out,
+            bool              transpose_matrix,
+            bool              add,
+            typename Number,
+            typename Number2>
+  std::enable_if_t<(variant == evaluate_symmetric_hierarchical), void>
+  apply_matrix_vector_product(const Number2 *matrix,
+                              const Number * in,
+                              Number *       out)
+  {
+    static_assert(n_rows > 0 && n_columns > 0,
+                  "Specialization requires n_rows, n_columns > 0");
 
-    const int stride    = Utilities::pow(n_columns, direction);
-    const int n_blocks1 = one_line ? 1 : stride;
-    const int n_blocks2 =
-      Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
+    constexpr bool evaluate_antisymmetric =
+      (quantity == EvaluatorQuantity::gradient);
 
-    const int offset = (n_columns + 1) / 2;
+    constexpr int mm     = transpose_matrix ? n_rows : n_columns,
+                  nn     = transpose_matrix ? n_columns : n_rows;
+    constexpr int n_cols = nn / 2;
+    constexpr int mid    = mm / 2;
 
-    // this code may look very inefficient at first sight due to the many
-    // different cases with if's at the innermost loop part, but all of the
-    // conditionals can be evaluated at compile time because they are
-    // templates, so the compiler should optimize everything away
-    for (int i2 = 0; i2 < n_blocks2; ++i2)
+    if (transpose_matrix)
       {
-        for (int i1 = 0; i1 < n_blocks1; ++i1)
+        std::array<Number, mm> x;
+        for (unsigned int i = 0; i < mm; ++i)
+          x[i] = in[stride_in * i];
+        for (unsigned int col = 0; col < n_cols; ++col)
           {
-            constexpr unsigned int mid_size =
-              (n_rows_static == -1 || n_columns_static == -1) ?
-                max_mid :
-                (mid_static > 0 ? mid_static : 1);
-            Number xp[mid_size], xm[mid_size];
-            for (int i = 0; i < mid; ++i)
+            Number r0, r1;
+            if (mid > 0)
               {
-                if (contract_over_rows == true && type == 1)
-                  {
-                    xp[i] = in[stride * i] - in[stride * (mm - 1 - i)];
-                    xm[i] = in[stride * i] + in[stride * (mm - 1 - i)];
-                  }
-                else
+                r0 = matrix[col] * x[0];
+                r1 = matrix[col + n_columns] * x[1];
+                for (unsigned int ind = 1; ind < mid; ++ind)
                   {
-                    xp[i] = in[stride * i] + in[stride * (mm - 1 - i)];
-                    xm[i] = in[stride * i] - in[stride * (mm - 1 - i)];
+                    r0 += matrix[col + 2 * ind * n_columns] * x[2 * ind];
+                    r1 +=
+                      matrix[col + (2 * ind + 1) * n_columns] * x[2 * ind + 1];
                   }
               }
-            Number xmid = in[stride * mid];
-            for (int col = 0; col < n_cols; ++col)
+            else
+              r0 = r1 = Number();
+            if (mm % 2 == 1)
+              r0 += matrix[col + (mm - 1) * n_columns] * x[mm - 1];
+            if (add)
               {
-                Number r0, r1;
-                if (mid > 0)
-                  {
-                    if (contract_over_rows == true)
-                      {
-                        r0 = shapes[col] * xp[0];
-                        r1 = shapes[(n_rows - 1) * offset + col] * xm[0];
-                      }
-                    else
-                      {
-                        r0 = shapes[col * offset] * xp[0];
-                        r1 = shapes[(n_rows - 1 - col) * offset] * xm[0];
-                      }
-                    for (int ind = 1; ind < mid; ++ind)
-                      {
-                        if (contract_over_rows == true)
-                          {
-                            r0 += shapes[ind * offset + col] * xp[ind];
-                            r1 += shapes[(n_rows - 1 - ind) * offset + col] *
-                                  xm[ind];
-                          }
-                        else
-                          {
-                            r0 += shapes[col * offset + ind] * xp[ind];
-                            r1 += shapes[(n_rows - 1 - col) * offset + ind] *
-                                  xm[ind];
-                          }
-                      }
-                  }
+                out[stride_out * col] += r0 + r1;
+                if (evaluate_antisymmetric)
+                  out[stride_out * (nn - 1 - col)] += r1 - r0;
                 else
-                  r0 = r1 = Number();
-                if (mm % 2 == 1 && contract_over_rows == true)
-                  {
-                    if (type == 1)
-                      r1 += shapes[mid * offset + col] * xmid;
-                    else
-                      r0 += shapes[mid * offset + col] * xmid;
-                  }
-                else if (mm % 2 == 1 && (nn % 2 == 0 || type > 0 || mm == 3))
-                  r0 += shapes[col * offset + mid] * xmid;
-
-                if (add)
-                  {
-                    out[stride * col] += r0 + r1;
-                    if (type == 1 && contract_over_rows == false)
-                      out[stride * (nn - 1 - col)] += r1 - r0;
-                    else
-                      out[stride * (nn - 1 - col)] += r0 - r1;
-                  }
+                  out[stride_out * (nn - 1 - col)] += r0 - r1;
+              }
+            else
+              {
+                out[stride_out * col] = r0 + r1;
+                if (evaluate_antisymmetric)
+                  out[stride_out * (nn - 1 - col)] = r1 - r0;
                 else
+                  out[stride_out * (nn - 1 - col)] = r0 - r1;
+              }
+          }
+        if (nn % 2 == 1)
+          {
+            Number             r0;
+            const unsigned int shift = evaluate_antisymmetric ? 1 : 0;
+            if (mid > 0)
+              {
+                r0 = matrix[n_cols + shift * n_columns] * x[shift];
+                for (unsigned int ind = 1; ind < mid; ++ind)
+                  r0 += matrix[n_cols + (2 * ind + shift) * n_columns] *
+                        x[2 * ind + shift];
+              }
+            else
+              r0 = 0;
+            if (!evaluate_antisymmetric && mm % 2 == 1)
+              r0 += matrix[n_cols + (mm - 1) * n_columns] * x[mm - 1];
+            if (add)
+              out[stride_out * n_cols] += r0;
+            else
+              out[stride_out * n_cols] = r0;
+          }
+      }
+    else
+      {
+        std::array<Number, mid + 1> xp, xm;
+        for (int i = 0; i < mid; ++i)
+          if (!evaluate_antisymmetric)
+            {
+              xp[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
+              xm[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
+            }
+          else
+            {
+              xp[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
+              xm[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
+            }
+        if (mm % 2 == 1)
+          xp[mid] = in[stride_in * mid];
+        for (unsigned int col = 0; col < n_cols; ++col)
+          {
+            Number r0, r1;
+            if (mid > 0)
+              {
+                r0 = matrix[2 * col * n_columns] * xp[0];
+                r1 = matrix[(2 * col + 1) * n_columns] * xm[0];
+                for (unsigned int ind = 1; ind < mid; ++ind)
                   {
-                    out[stride * col] = r0 + r1;
-                    if (type == 1 && contract_over_rows == false)
-                      out[stride * (nn - 1 - col)] = r1 - r0;
-                    else
-                      out[stride * (nn - 1 - col)] = r0 - r1;
+                    r0 += matrix[2 * col * n_columns + ind] * xp[ind];
+                    r1 += matrix[(2 * col + 1) * n_columns + ind] * xm[ind];
                   }
               }
-            if (type == 0 && contract_over_rows == true && nn % 2 == 1 &&
-                mm % 2 == 1 && mm > 3)
-              {
-                if (add)
-                  out[stride * n_cols] += shapes[mid * offset + n_cols] * xmid;
-                else
-                  out[stride * n_cols] = shapes[mid * offset + n_cols] * xmid;
-              }
-            else if (contract_over_rows == true && nn % 2 == 1)
+            else
+              r0 = r1 = Number();
+            if (mm % 2 == 1)
               {
-                Number r0;
-                if (mid > 0)
-                  {
-                    r0 = shapes[n_cols] * xp[0];
-                    for (int ind = 1; ind < mid; ++ind)
-                      r0 += shapes[ind * offset + n_cols] * xp[ind];
-                  }
-                else
-                  r0 = Number();
-                if (type != 1 && mm % 2 == 1)
-                  r0 += shapes[mid * offset + n_cols] * xmid;
-
-                if (add)
-                  out[stride * n_cols] += r0;
+                if (evaluate_antisymmetric)
+                  r1 += matrix[(2 * col + 1) * n_columns + mid] * xp[mid];
                 else
-                  out[stride * n_cols] = r0;
+                  r0 += matrix[2 * col * n_columns + mid] * xp[mid];
               }
-            else if (contract_over_rows == false && nn % 2 == 1)
+            if (add)
               {
-                Number r0;
-                if (mid > 0)
-                  {
-                    if (type == 1)
-                      {
-                        r0 = shapes[n_cols * offset] * xm[0];
-                        for (int ind = 1; ind < mid; ++ind)
-                          r0 += shapes[n_cols * offset + ind] * xm[ind];
-                      }
-                    else
-                      {
-                        r0 = shapes[n_cols * offset] * xp[0];
-                        for (int ind = 1; ind < mid; ++ind)
-                          r0 += shapes[n_cols * offset + ind] * xp[ind];
-                      }
-                  }
-                else
-                  r0 = Number();
-
-                if ((type == 0 || type == 2) && mm % 2 == 1)
-                  r0 += shapes[n_cols * offset + mid] * xmid;
-
-                if (add)
-                  out[stride * n_cols] += r0;
-                else
-                  out[stride * n_cols] = r0;
+                out[stride_out * (2 * col)] += r0;
+                out[stride_out * (2 * col + 1)] += r1;
               }
-            if (one_line == false)
+            else
               {
-                in += 1;
-                out += 1;
+                out[stride_out * (2 * col)]     = r0;
+                out[stride_out * (2 * col + 1)] = r1;
               }
           }
-        if (one_line == false)
+        if (nn % 2 == 1)
           {
-            in += stride * (mm - 1);
-            out += stride * (nn - 1);
+            Number r0;
+            if (mid > 0)
+              {
+                r0 = matrix[(nn - 1) * n_columns] * xp[0];
+                for (unsigned int ind = 1; ind < mid; ++ind)
+                  r0 += matrix[(nn - 1) * n_columns + ind] * xp[ind];
+              }
+            else
+              r0 = Number();
+            if (mm % 2 == 1 && !evaluate_antisymmetric)
+              r0 += matrix[(nn - 1) * n_columns + mid] * xp[mid];
+            if (add)
+              out[stride_out * (nn - 1)] += r0;
+            else
+              out[stride_out * (nn - 1)] = r0;
           }
       }
   }
@@ -1795,22 +1057,13 @@ namespace internal
 
 
   /**
-   * Internal evaluator for 1d-3d shape function using the tensor product form
-   * of the basis functions.
-   *
-   * This class implements a different approach to the symmetric case for
-   * values, gradients, and Hessians also treated with the above functions: It
-   * is possible to reduce the cost per dimension from N^2 to N^2/2, where N
-   * is the number of 1d dofs (there are only N^2/2 different entries in the
-   * shape matrix, so this is plausible). The approach is based on the idea of
-   * applying the operator on the even and odd part of the input vectors
-   * separately, given that the shape functions evaluated on quadrature points
-   * are symmetric. This method is presented e.g. in the book "Implementing
-   * Spectral Methods for Partial Differential Equations" by David A. Kopriva,
-   * Springer, 2009, section 3.5.3 (Even-Odd-Decomposition). Even though the
-   * experiments in the book say that the method is not efficient for N<20, it
-   * is more efficient in the context where the loop bounds are compile-time
-   * constants (templates).
+   * Generic evaluator framework that valuates the given shape data in general
+   * dimensions using the tensor product form. Depending on the particular
+   * layout in the matrix entries, this corresponds to a usual matrix-matrix
+   * product or a matrix-matrix product including some symmetries. The actual
+   * work is implemented by functions of type apply_matrix_vector_product
+   * working on a single dimension, controlled by suitable strides, using the
+   * kernel specified via variant.
    *
    * @tparam dim Space dimension in which this class is applied
    * @tparam n_rows Number of rows in the transformation matrix, which corresponds
@@ -1825,17 +1078,13 @@ namespace internal
    *                 operator* with Number and produce Number as an output to
    *                 be a valid type
    */
-  template <int dim,
-            int n_rows,
-            int n_columns,
+  template <EvaluatorVariant variant,
+            int              dim,
+            int              n_rows,
+            int              n_columns,
             typename Number,
-            typename Number2>
-  struct EvaluatorTensorProduct<evaluate_evenodd,
-                                dim,
-                                n_rows,
-                                n_columns,
-                                Number,
-                                Number2>
+            typename Number2 = Number>
+  struct EvaluatorTensorProduct
   {
     static constexpr unsigned int n_rows_of_product =
       Utilities::pow(n_rows, dim);
@@ -1844,8 +1093,7 @@ namespace internal
 
     /**
      * Empty constructor. Does nothing. Be careful when using 'values' and
-     * related methods because they need to be filled with the other
-     * constructor passing in at least an array for the values.
+     * related methods because they need to be filled with the other pointer
      */
     EvaluatorTensorProduct()
       : shape_values(nullptr)
@@ -1854,38 +1102,62 @@ namespace internal
     {}
 
     /**
-     * Constructor, taking the data from ShapeInfo (using the even-odd
-     * variants stored there)
-     */
-    EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values)
-      : shape_values(shape_values.begin())
-      , shape_gradients(nullptr)
-      , shape_hessians(nullptr)
-    {
-      AssertDimension(shape_values.size(), n_rows * ((n_columns + 1) / 2));
-    }
-
-    /**
-     * Constructor, taking the data from ShapeInfo (using the even-odd
-     * variants stored there)
+     * Constructor, taking the data from ShapeInfo
      */
     EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
                            const AlignedVector<Number2> &shape_gradients,
                            const AlignedVector<Number2> &shape_hessians,
-                           const unsigned int            dummy1 = 0,
-                           const unsigned int            dummy2 = 0)
+                           const unsigned int = 0,
+                           const unsigned int = 0)
       : shape_values(shape_values.begin())
       , shape_gradients(shape_gradients.begin())
       , shape_hessians(shape_hessians.begin())
     {
-      // In this function, we allow for dummy pointers if some of values,
-      // gradients or hessians should not be computed
-      if (!shape_values.empty())
-        AssertDimension(shape_values.size(), n_rows * ((n_columns + 1) / 2));
-      if (!shape_gradients.empty())
-        AssertDimension(shape_gradients.size(), n_rows * ((n_columns + 1) / 2));
-      if (!shape_hessians.empty())
-        AssertDimension(shape_hessians.size(), n_rows * ((n_columns + 1) / 2));
+      if (variant == evaluate_evenodd)
+        {
+          if (!shape_values.empty())
+            AssertDimension(shape_values.size(),
+                            n_rows * ((n_columns + 1) / 2));
+          if (!shape_gradients.empty())
+            AssertDimension(shape_gradients.size(),
+                            n_rows * ((n_columns + 1) / 2));
+          if (!shape_hessians.empty())
+            AssertDimension(shape_hessians.size(),
+                            n_rows * ((n_columns + 1) / 2));
+        }
+      else
+        {
+          // We can enter this function either for the apply() path that has
+          // n_rows * n_columns entries or for the apply_face() path that only
+          // has n_rows * 3 entries in the array. Since we cannot decide about
+          // the use we must allow for both here.
+          Assert(shape_values.empty() ||
+                   shape_values.size() == n_rows * n_columns ||
+                   shape_values.size() == 3 * n_rows,
+                 ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
+          Assert(shape_gradients.empty() ||
+                   shape_gradients.size() == n_rows * n_columns,
+                 ExcDimensionMismatch(shape_gradients.size(),
+                                      n_rows * n_columns));
+          Assert(shape_hessians.empty() ||
+                   shape_hessians.size() == n_rows * n_columns,
+                 ExcDimensionMismatch(shape_hessians.size(),
+                                      n_rows * n_columns));
+        }
+    }
+
+    /**
+     * Constructor, taking the data from ShapeInfo via raw pointers
+     */
+    EvaluatorTensorProduct(const Number2 *    shape_values,
+                           const Number2 *    shape_gradients,
+                           const Number2 *    shape_hessians,
+                           const unsigned int dummy1 = 0,
+                           const unsigned int dummy2 = 0)
+      : shape_values(shape_values)
+      , shape_gradients(shape_gradients)
+      , shape_hessians(shape_hessians)
+    {
       (void)dummy1;
       (void)dummy2;
     }
@@ -1894,24 +1166,35 @@ namespace internal
     void
     values(const Number in[], Number out[]) const
     {
-      Assert(shape_values != nullptr, ExcNotInitialized());
-      apply<direction, contract_over_rows, add, 0>(shape_values, in, out);
+      apply<direction,
+            contract_over_rows,
+            add,
+            false,
+            EvaluatorQuantity::value>(shape_values, in, out);
     }
 
     template <int direction, bool contract_over_rows, bool add>
     void
     gradients(const Number in[], Number out[]) const
     {
-      Assert(shape_gradients != nullptr, ExcNotInitialized());
-      apply<direction, contract_over_rows, add, 1>(shape_gradients, in, out);
+      constexpr EvaluatorQuantity gradient_type =
+        (variant == evaluate_general ? EvaluatorQuantity::value :
+                                       EvaluatorQuantity::gradient);
+      apply<direction, contract_over_rows, add, false, gradient_type>(
+        shape_gradients, in, out);
     }
 
     template <int direction, bool contract_over_rows, bool add>
     void
     hessians(const Number in[], Number out[]) const
     {
-      Assert(shape_hessians != nullptr, ExcNotInitialized());
-      apply<direction, contract_over_rows, add, 2>(shape_hessians, in, out);
+      constexpr EvaluatorQuantity hessian_type =
+        (((variant == evaluate_general) |
+          (variant == evaluate_symmetric_hierarchical)) ?
+           EvaluatorQuantity::value :
+           EvaluatorQuantity::hessian);
+      apply<direction, contract_over_rows, add, false, hessian_type>(
+        shape_hessians, in, out);
     }
 
     template <int direction, bool contract_over_rows, bool add>
@@ -1919,7 +1202,8 @@ namespace internal
     values_one_line(const Number in[], Number out[]) const
     {
       Assert(shape_values != nullptr, ExcNotInitialized());
-      apply<direction, contract_over_rows, add, 0, true>(shape_values, in, out);
+      apply<direction, contract_over_rows, add, true, EvaluatorQuantity::value>(
+        shape_values, in, out);
     }
 
     template <int direction, bool contract_over_rows, bool add>
@@ -1927,9 +1211,11 @@ namespace internal
     gradients_one_line(const Number in[], Number out[]) const
     {
       Assert(shape_gradients != nullptr, ExcNotInitialized());
-      apply<direction, contract_over_rows, add, 1, true>(shape_gradients,
-                                                         in,
-                                                         out);
+      constexpr EvaluatorQuantity gradient_type =
+        (variant == evaluate_general ? EvaluatorQuantity::value :
+                                       EvaluatorQuantity::gradient);
+      apply<direction, contract_over_rows, add, true, gradient_type>(
+        shape_gradients, in, out);
     }
 
     template <int direction, bool contract_over_rows, bool add>
@@ -1937,9 +1223,13 @@ namespace internal
     hessians_one_line(const Number in[], Number out[]) const
     {
       Assert(shape_hessians != nullptr, ExcNotInitialized());
-      apply<direction, contract_over_rows, add, 2, true>(shape_hessians,
-                                                         in,
-                                                         out);
+      constexpr EvaluatorQuantity hessian_type =
+        (((variant == evaluate_general) |
+          (variant == evaluate_symmetric_hierarchical)) ?
+           EvaluatorQuantity::value :
+           EvaluatorQuantity::hessian);
+      apply<direction, contract_over_rows, add, true, hessian_type>(
+        shape_hessians, in, out);
     }
 
     /**
@@ -1948,8 +1238,8 @@ namespace internal
      * data in the input array. This function allows the @p in and @p out
      * arrays to alias for the case n_rows == n_columns, i.e., it is safe to
      * perform the contraction in place where @p in and @p out point to the
-     * same address. For the case n_rows != n_columns, the output is only
-     * correct if @p one_line is set to true.
+     * same address. For the case n_rows != n_columns, the output is in general
+     * not correct.
      *
      * @tparam direction Direction that is evaluated
      * @tparam contract_over_rows If true, the tensor contraction sums
@@ -1957,10 +1247,6 @@ namespace internal
      *                            array, otherwise it sums over the columns
      * @tparam add If true, the result is added to the output vector, else
      *             the computed values overwrite the content in the output
-     * @tparam type Determines whether to use the symmetries appearing in
-     *              shape values (type=0), shape gradients (type=1) or
-     *              second derivatives (type=2, similar to type 0 but
-     *              without two additional zero entries)
      * @tparam one_line If true, the kernel is only applied along a single 1d
      *                  stripe within a dim-dimensional tensor, not the full
      *                  n_rows^dim points as in the @p false case.
@@ -1973,24 +1259,49 @@ namespace internal
     template <int  direction,
               bool contract_over_rows,
               bool add,
-              int  type,
-              bool one_line = false>
+              bool one_line     = false,
+              EvaluatorQuantity = EvaluatorQuantity::value>
     static void
     apply(const Number2 *DEAL_II_RESTRICT shape_data,
           const Number *                  in,
-          Number *                        out)
-    {
-      even_odd_apply<dim,
-                     n_rows,
-                     n_columns,
-                     Number,
-                     Number2,
-                     direction,
-                     contract_over_rows,
-                     add,
-                     type,
-                     one_line>(n_rows, n_columns, shape_data, in, out);
-    }
+          Number *                        out);
+
+    /**
+     * This function applies the tensor product operation to produce face values
+     * from cell values. As opposed to the apply method, this method assumes
+     * that the directions orthogonal to the face have n_rows degrees of
+     * freedom per direction and not n_columns for those directions lower than
+     * the one currently applied. In other words, apply_face() must be called
+     * before calling any interpolation within the face.
+     *
+     * @tparam face_direction Direction of the normal vector (0=x, 1=y, etc)
+     * @tparam contract_onto_face If true, the input vector is of size n_rows^dim
+     *                            and interpolation into n_rows^(dim-1) points
+     *                            is performed. This is a typical scenario in
+     *                            FEFaceEvaluation::evaluate() calls. If false,
+     *                            data from n_rows^(dim-1) points is expanded
+     *                            into the n_rows^dim points of the higher-
+     *                            dimensional data array. Derivatives in the
+     *                            case contract_onto_face==false are summed
+     *                            together
+     * @tparam add If true, the result is added to the output vector, else
+     *             the computed values overwrite the content in the output
+     * @tparam max_derivative Sets the number of derivatives that should be
+     *             computed. 0 means only values, 1 means values and first
+     *             derivatives, 2 second derivates. Note that all the
+     *             derivatives access the data in @p shape_values passed to
+     *             the constructor of the class
+     *
+     * @param in address of the input data vector
+     * @param out address of the output data vector
+     */
+    template <int  face_direction,
+              bool contract_onto_face,
+              bool add,
+              int  max_derivative>
+    void
+    apply_face(const Number *DEAL_II_RESTRICT in,
+               Number *DEAL_II_RESTRICT       out) const;
 
   private:
     const Number2 *shape_values;
@@ -1999,243 +1310,336 @@ namespace internal
   };
 
 
-  /**
-   * Internal evaluator for shape function using the tensor product form
-   * of the basis functions. The same as the other templated class but
-   * without making use of template arguments and variable loop bounds
-   * instead.
-   */
-  template <int dim, typename Number, typename Number2>
-  struct EvaluatorTensorProduct<evaluate_evenodd, dim, 0, 0, Number, Number2>
+
+  template <EvaluatorVariant variant,
+            int              dim,
+            int              n_rows,
+            int              n_columns,
+            typename Number,
+            typename Number2>
+  template <int               direction,
+            bool              contract_over_rows,
+            bool              add,
+            bool              one_line,
+            EvaluatorQuantity quantity>
+  inline void
+  EvaluatorTensorProduct<variant, dim, n_rows, n_columns, Number, Number2>::
+    apply(const Number2 *DEAL_II_RESTRICT shape_data,
+          const Number *                  in,
+          Number *                        out)
   {
-    EvaluatorTensorProduct()
-      : shape_values(nullptr)
-      , shape_gradients(nullptr)
-      , shape_hessians(nullptr)
-      , n_rows(numbers::invalid_unsigned_int)
-      , n_columns(numbers::invalid_unsigned_int)
-    {}
+    static_assert(one_line == false || direction == dim - 1,
+                  "Single-line evaluation only works for direction=dim-1.");
+    Assert(shape_data != nullptr,
+           ExcMessage(
+             "The given array shape_data must not be the null pointer!"));
+    Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
+             in != out,
+           ExcMessage("In-place operation only supported for "
+                      "n_rows==n_columns or single-line interpolation"));
+    AssertIndexRange(direction, dim);
+    constexpr int mm = contract_over_rows ? n_rows : n_columns,
+                  nn = contract_over_rows ? n_columns : n_rows;
 
-    EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
-                           const unsigned int            n_rows    = 0,
-                           const unsigned int            n_columns = 0)
-      : shape_values(shape_values.begin())
-      , shape_gradients(nullptr)
-      , shape_hessians(nullptr)
-      , n_rows(n_rows)
-      , n_columns(n_columns)
-    {
-      AssertDimension(shape_values.size(), n_rows * ((n_columns + 1) / 2));
-    }
+    constexpr int stride    = Utilities::pow(n_columns, direction);
+    constexpr int n_blocks1 = one_line ? 1 : stride;
+    constexpr int n_blocks2 =
+      Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
 
-    EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
-                           const AlignedVector<Number2> &shape_gradients,
-                           const AlignedVector<Number2> &shape_hessians,
-                           const unsigned int            n_rows    = 0,
-                           const unsigned int            n_columns = 0)
-      : shape_values(shape_values.begin())
-      , shape_gradients(shape_gradients.begin())
-      , shape_hessians(shape_hessians.begin())
-      , n_rows(n_rows)
-      , n_columns(n_columns)
-    {
-      if (!shape_values.empty())
-        AssertDimension(shape_values.size(), n_rows * ((n_columns + 1) / 2));
-      if (!shape_gradients.empty())
-        AssertDimension(shape_gradients.size(), n_rows * ((n_columns + 1) / 2));
-      if (!shape_hessians.empty())
-        AssertDimension(shape_hessians.size(), n_rows * ((n_columns + 1) / 2));
-    }
+    for (int i2 = 0; i2 < n_blocks2; ++i2)
+      {
+        for (int i1 = 0; i1 < n_blocks1; ++i1)
+          {
+            apply_matrix_vector_product<variant,
+                                        quantity,
+                                        n_rows,
+                                        n_columns,
+                                        stride,
+                                        stride,
+                                        contract_over_rows,
+                                        add>(shape_data, in, out);
 
-    template <int direction, bool contract_over_rows, bool add>
-    void
-    values(const Number in[], Number out[]) const
-    {
-      Assert(shape_values != nullptr, ExcNotInitialized());
-      apply<direction, contract_over_rows, add, 0>(shape_values, in, out);
-    }
+            if (one_line == false)
+              {
+                ++in;
+                ++out;
+              }
+          }
+        if (one_line == false)
+          {
+            in += stride * (mm - 1);
+            out += stride * (nn - 1);
+          }
+      }
+  }
 
-    template <int direction, bool contract_over_rows, bool add>
-    void
-    gradients(const Number in[], Number out[]) const
-    {
-      Assert(shape_gradients != nullptr, ExcNotInitialized());
-      apply<direction, contract_over_rows, add, 1>(shape_gradients, in, out);
-    }
 
-    template <int direction, bool contract_over_rows, bool add>
-    void
-    hessians(const Number in[], Number out[]) const
-    {
-      Assert(shape_hessians != nullptr, ExcNotInitialized());
-      apply<direction, contract_over_rows, add, 2>(shape_hessians, in, out);
-    }
 
-    template <int direction, bool contract_over_rows, bool add>
-    void
-    values_one_line(const Number in[], Number out[]) const
-    {
-      Assert(shape_values != nullptr, ExcNotInitialized());
-      apply<direction, contract_over_rows, add, 0, true>(shape_values, in, out);
-    }
+  template <EvaluatorVariant variant,
+            int              dim,
+            int              n_rows,
+            int              n_columns,
+            typename Number,
+            typename Number2>
+  template <int  face_direction,
+            bool contract_onto_face,
+            bool add,
+            int  max_derivative>
+  inline void
+  EvaluatorTensorProduct<variant, dim, n_rows, n_columns, Number, Number2>::
+    apply_face(const Number *DEAL_II_RESTRICT in,
+               Number *DEAL_II_RESTRICT       out) const
+  {
+    Assert(dim > 0, ExcMessage("Only dim=1,2,3 supported"));
+    static_assert(max_derivative >= 0 && max_derivative < 3,
+                  "Only derivative orders 0-2 implemented");
+    Assert(shape_values != nullptr,
+           ExcMessage(
+             "The given array shape_values must not be the null pointer."));
 
-    template <int direction, bool contract_over_rows, bool add>
-    void
-    gradients_one_line(const Number in[], Number out[]) const
-    {
-      Assert(shape_gradients != nullptr, ExcNotInitialized());
-      apply<direction, contract_over_rows, add, 1, true>(shape_gradients,
-                                                         in,
-                                                         out);
-    }
+    constexpr int n_blocks1 = (dim > 1 ? n_rows : 1);
+    constexpr int n_blocks2 = (dim > 2 ? n_rows : 1);
 
-    template <int direction, bool contract_over_rows, bool add>
-    void
-    hessians_one_line(const Number in[], Number out[]) const
-    {
-      Assert(shape_hessians != nullptr, ExcNotInitialized());
-      apply<direction, contract_over_rows, add, 2, true>(shape_hessians,
-                                                         in,
-                                                         out);
-    }
+    AssertIndexRange(face_direction, dim);
+    constexpr int in_stride  = Utilities::pow(n_rows, face_direction);
+    constexpr int out_stride = Utilities::pow(n_rows, dim - 1);
+    const Number2 *DEAL_II_RESTRICT shape_values = this->shape_values;
 
-    template <int  direction,
-              bool contract_over_rows,
-              bool add,
-              int  type,
-              bool one_line = false>
-    void
-    apply(const Number2 *DEAL_II_RESTRICT shape_data,
-          const Number *                  in,
-          Number *                        out) const
-    {
-      even_odd_apply<dim,
-                     -1,
-                     -1,
-                     Number,
-                     Number2,
-                     direction,
-                     contract_over_rows,
-                     add,
-                     type,
-                     one_line>(n_rows, n_columns, shape_data, in, out);
-    }
+    for (int i2 = 0; i2 < n_blocks2; ++i2)
+      {
+        for (int i1 = 0; i1 < n_blocks1; ++i1)
+          {
+            if (contract_onto_face == true)
+              {
+                Number res0 = shape_values[0] * in[0];
+                Number res1, res2;
+                if (max_derivative > 0)
+                  res1 = shape_values[n_rows] * in[0];
+                if (max_derivative > 1)
+                  res2 = shape_values[2 * n_rows] * in[0];
+                for (int ind = 1; ind < n_rows; ++ind)
+                  {
+                    res0 += shape_values[ind] * in[in_stride * ind];
+                    if (max_derivative > 0)
+                      res1 += shape_values[ind + n_rows] * in[in_stride * ind];
+                    if (max_derivative > 1)
+                      res2 +=
+                        shape_values[ind + 2 * n_rows] * in[in_stride * ind];
+                  }
+                if (add)
+                  {
+                    out[0] += res0;
+                    if (max_derivative > 0)
+                      out[out_stride] += res1;
+                    if (max_derivative > 1)
+                      out[2 * out_stride] += res2;
+                  }
+                else
+                  {
+                    out[0] = res0;
+                    if (max_derivative > 0)
+                      out[out_stride] = res1;
+                    if (max_derivative > 1)
+                      out[2 * out_stride] = res2;
+                  }
+              }
+            else
+              {
+                for (int col = 0; col < n_rows; ++col)
+                  {
+                    if (add)
+                      out[col * in_stride] += shape_values[col] * in[0];
+                    else
+                      out[col * in_stride] = shape_values[col] * in[0];
+                    if (max_derivative > 0)
+                      out[col * in_stride] +=
+                        shape_values[col + n_rows] * in[out_stride];
+                    if (max_derivative > 1)
+                      out[col * in_stride] +=
+                        shape_values[col + 2 * n_rows] * in[2 * out_stride];
+                  }
+              }
 
-  private:
-    const Number2 *    shape_values;
-    const Number2 *    shape_gradients;
-    const Number2 *    shape_hessians;
-    const unsigned int n_rows;
-    const unsigned int n_columns;
-  };
+            // increment: in regular case, just go to the next point in
+            // x-direction. If we are at the end of one chunk in x-dir, need
+            // to jump over to the next layer in z-direction
+            switch (face_direction)
+              {
+                case 0:
+                  in += contract_onto_face ? n_rows : 1;
+                  out += contract_onto_face ? 1 : n_rows;
+                  break;
+                case 1:
+                  ++in;
+                  ++out;
+                  // faces 2 and 3 in 3d use local coordinate system zx, which
+                  // is the other way around compared to the tensor
+                  // product. Need to take that into account.
+                  if (dim == 3)
+                    {
+                      if (contract_onto_face)
+                        out += n_rows - 1;
+                      else
+                        in += n_rows - 1;
+                    }
+                  break;
+                case 2:
+                  ++in;
+                  ++out;
+                  break;
+                default:
+                  Assert(false, ExcNotImplemented());
+              }
+          }
+
+        // adjust for local coordinate system zx
+        if (face_direction == 1 && dim == 3)
+          {
+            if (contract_onto_face)
+              {
+                in += n_rows * (n_rows - 1);
+                out -= n_rows * n_rows - 1;
+              }
+            else
+              {
+                out += n_rows * (n_rows - 1);
+                in -= n_rows * n_rows - 1;
+              }
+          }
+      }
+  }
 
 
 
   /**
-   * Internal evaluator for 1d-3d shape function using the tensor product form
-   * of the basis functions.
-   *
-   * This class implements an approach similar to the even-odd decomposition
-   * but with a different type of symmetry. In this case, we assume that a
-   * single shape function already shows the symmetry over the quadrature
-   * points, rather than the complete basis that is considered in the even-odd
-   * case. In particular, we assume that the shape functions are ordered as in
-   * the Legendre basis, with symmetric shape functions in the even slots
-   * (rows of the values array) and point-symmetric in the odd slots. Like the
-   * even-odd decomposition, the number of operations are N^2/2 rather than
-   * N^2 FMAs (fused multiply-add), where N is the number of 1d dofs. The
-   * difference is in the way the input and output quantities are symmetrized.
+   * Internal evaluator for shape function using the tensor product form
+   * of the basis functions. The same as the other templated class but
+   * without making use of template arguments and variable loop bounds
+   * instead.
    *
    * @tparam dim Space dimension in which this class is applied
-   * @tparam n_rows Number of rows in the transformation matrix, which corresponds
-   *                to the number of 1d shape functions in the usual tensor
-   *                contraction setting
-   * @tparam n_columns Number of columns in the transformation matrix, which
-   *                   corresponds to the number of 1d shape functions in the
-   *                   usual tensor contraction setting
    * @tparam Number Abstract number type for input and output arrays
    * @tparam Number2 Abstract number type for coefficient arrays (defaults to
    *                 same type as the input/output arrays); must implement
    *                 operator* with Number and produce Number as an output to
    *                 be a valid type
    */
-  template <int dim,
-            int n_rows,
-            int n_columns,
+  template <EvaluatorVariant variant,
+            int              dim,
             typename Number,
             typename Number2>
-  struct EvaluatorTensorProduct<evaluate_symmetric_hierarchical,
-                                dim,
-                                n_rows,
-                                n_columns,
-                                Number,
-                                Number2>
+  struct EvaluatorTensorProduct<variant, dim, 0, 0, Number, Number2>
   {
     static constexpr unsigned int n_rows_of_product =
-      Utilities::pow(n_rows, dim);
+      numbers::invalid_unsigned_int;
     static constexpr unsigned int n_columns_of_product =
-      Utilities::pow(n_columns, dim);
+      numbers::invalid_unsigned_int;
 
     /**
      * Empty constructor. Does nothing. Be careful when using 'values' and
-     * related methods because they need to be filled with the other
-     * constructor passing in at least an array for the values.
+     * related methods because they need to be filled with the other constructor
      */
     EvaluatorTensorProduct()
       : shape_values(nullptr)
       , shape_gradients(nullptr)
       , shape_hessians(nullptr)
+      , n_rows(numbers::invalid_unsigned_int)
+      , n_columns(numbers::invalid_unsigned_int)
     {}
 
     /**
-     * Constructor, taking the data from ShapeInfo (using the even-odd
-     * variants stored there)
-     */
-    EvaluatorTensorProduct(const AlignedVector<Number> &shape_values)
-      : shape_values(shape_values.begin())
-      , shape_gradients(nullptr)
-      , shape_hessians(nullptr)
-    {}
-
-    /**
-     * Constructor, taking the data from ShapeInfo (using the even-odd
-     * variants stored there)
+     * Constructor, taking the data from ShapeInfo
      */
     EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
                            const AlignedVector<Number2> &shape_gradients,
                            const AlignedVector<Number2> &shape_hessians,
-                           const unsigned int            dummy1 = 0,
-                           const unsigned int            dummy2 = 0)
+                           const unsigned int            n_rows    = 0,
+                           const unsigned int            n_columns = 0)
       : shape_values(shape_values.begin())
       , shape_gradients(shape_gradients.begin())
       , shape_hessians(shape_hessians.begin())
+      , n_rows(n_rows)
+      , n_columns(n_columns)
     {
-      (void)dummy1;
-      (void)dummy2;
+      if (variant == evaluate_evenodd)
+        {
+          if (!shape_values.empty())
+            AssertDimension(shape_values.size(),
+                            n_rows * ((n_columns + 1) / 2));
+          if (!shape_gradients.empty())
+            AssertDimension(shape_gradients.size(),
+                            n_rows * ((n_columns + 1) / 2));
+          if (!shape_hessians.empty())
+            AssertDimension(shape_hessians.size(),
+                            n_rows * ((n_columns + 1) / 2));
+        }
+      else
+        {
+          // We can enter this function either for the apply() path that has
+          // n_rows * n_columns entries or for the apply_face() path that only
+          // has n_rows * 3 entries in the array. Since we cannot decide about
+          // the use we must allow for both here.
+          Assert(shape_values.empty() ||
+                   shape_values.size() == n_rows * n_columns ||
+                   shape_values.size() == n_rows * 3,
+                 ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
+          Assert(shape_gradients.empty() ||
+                   shape_gradients.size() == n_rows * n_columns,
+                 ExcDimensionMismatch(shape_gradients.size(),
+                                      n_rows * n_columns));
+          Assert(shape_hessians.empty() ||
+                   shape_hessians.size() == n_rows * n_columns,
+                 ExcDimensionMismatch(shape_hessians.size(),
+                                      n_rows * n_columns));
+        }
     }
 
+    /**
+     * Constructor, taking the data from ShapeInfo
+     */
+    EvaluatorTensorProduct(const Number2 *    shape_values,
+                           const Number2 *    shape_gradients,
+                           const Number2 *    shape_hessians,
+                           const unsigned int n_rows    = 0,
+                           const unsigned int n_columns = 0)
+      : shape_values(shape_values)
+      , shape_gradients(shape_gradients)
+      , shape_hessians(shape_hessians)
+      , n_rows(n_rows)
+      , n_columns(n_columns)
+    {}
+
     template <int direction, bool contract_over_rows, bool add>
     void
-    values(const Number in[], Number out[]) const
+    values(const Number *in, Number *out) const
     {
-      Assert(shape_values != nullptr, ExcNotInitialized());
-      apply<direction, contract_over_rows, add, 0>(shape_values, in, out);
+      apply<direction,
+            contract_over_rows,
+            add,
+            false,
+            EvaluatorQuantity::value>(shape_values, in, out);
     }
 
     template <int direction, bool contract_over_rows, bool add>
     void
-    gradients(const Number in[], Number out[]) const
+    gradients(const Number *in, Number *out) const
     {
-      Assert(shape_gradients != nullptr, ExcNotInitialized());
-      apply<direction, contract_over_rows, add, 1>(shape_gradients, in, out);
+      constexpr EvaluatorQuantity gradient_type =
+        (variant != evaluate_evenodd ? EvaluatorQuantity::value :
+                                       EvaluatorQuantity::gradient);
+      apply<direction, contract_over_rows, add, false, gradient_type>(
+        shape_gradients, in, out);
     }
 
     template <int direction, bool contract_over_rows, bool add>
     void
-    hessians(const Number in[], Number out[]) const
+    hessians(const Number *in, Number *out) const
     {
-      Assert(shape_hessians != nullptr, ExcNotInitialized());
-      apply<direction, contract_over_rows, add, 0>(shape_hessians, in, out);
+      constexpr EvaluatorQuantity hessian_type =
+        (variant != evaluate_evenodd ? EvaluatorQuantity::value :
+                                       EvaluatorQuantity::hessian);
+      apply<direction, contract_over_rows, add, false, hessian_type>(
+        shape_hessians, in, out);
     }
 
     template <int direction, bool contract_over_rows, bool add>
@@ -2243,274 +1647,297 @@ namespace internal
     values_one_line(const Number in[], Number out[]) const
     {
       Assert(shape_values != nullptr, ExcNotInitialized());
-      apply<direction, contract_over_rows, add, 0, true>(shape_values, in, out);
-    }
-
-    template <int direction, bool contract_over_rows, bool add>
-    void
-    gradients_one_line(const Number in[], Number out[]) const
-    {
-      Assert(shape_gradients != nullptr, ExcNotInitialized());
-      apply<direction, contract_over_rows, add, 1, true>(shape_gradients,
-                                                         in,
-                                                         out);
+      apply<direction, contract_over_rows, add, true, EvaluatorQuantity::value>(
+        shape_values, in, out);
     }
 
     template <int direction, bool contract_over_rows, bool add>
-    void
-    hessians_one_line(const Number in[], Number out[]) const
-    {
-      Assert(shape_hessians != nullptr, ExcNotInitialized());
-      apply<direction, contract_over_rows, add, 0, true>(shape_hessians,
-                                                         in,
-                                                         out);
-    }
-
-    /**
-     * This function applies the tensor product kernel, corresponding to a
-     * multiplication of 1d stripes, along the given @p direction of the tensor
-     * data in the input array. This function allows the @p in and @p out
-     * arrays to alias for the case n_rows == n_columns, i.e., it is safe to
-     * perform the contraction in place where @p in and @p out point to the
-     * same address. For the case n_rows != n_columns, the output is only
-     * correct if @p one_line is set to true.
-     *
-     * @tparam direction Direction that is evaluated
-     * @tparam contract_over_rows If true, the tensor contraction sums
-     *                            over the rows in the given @p shape_data
-     *                            array, otherwise it sums over the columns
-     * @tparam add If true, the result is added to the output vector, else
-     *             the computed values overwrite the content in the output
-     * @tparam type Determines whether the evaluation is symmetric in even
-     *              rows (type=0) or odd rows (type=1) of @p shape_data and
-     *              skew-symmetric in odd rows (type=0) or even rows (type=1)
-     * @tparam one_line If true, the kernel is only applied along a single 1d
-     *                  stripe within a dim-dimensional tensor, not the full
-     *                  n_rows^dim points as in the @p false case.
-     *
-     * @param shape_data Transformation matrix with @p n_rows rows and
-     *                   @p n_columns columns, stored in row-major format
-     * @param in Pointer to the start of the input data vector
-     * @param out Pointer to the start of the output data vector
-     */
-    template <int  direction,
-              bool contract_over_rows,
-              bool add,
-              int  type,
-              bool one_line = false>
-    static void
+    void
+    gradients_one_line(const Number in[], Number out[]) const
+    {
+      Assert(shape_gradients != nullptr, ExcNotInitialized());
+      constexpr EvaluatorQuantity gradient_type =
+        (variant != evaluate_evenodd ? EvaluatorQuantity::value :
+                                       EvaluatorQuantity::gradient);
+      apply<direction, contract_over_rows, add, true, gradient_type>(
+        shape_gradients, in, out);
+    }
+
+    template <int direction, bool contract_over_rows, bool add>
+    void
+    hessians_one_line(const Number in[], Number out[]) const
+    {
+      Assert(shape_hessians != nullptr, ExcNotInitialized());
+      constexpr EvaluatorQuantity hessian_type =
+        (variant != evaluate_evenodd ? EvaluatorQuantity::value :
+                                       EvaluatorQuantity::hessian);
+      apply<direction, contract_over_rows, add, true, hessian_type>(
+        shape_hessians, in, out);
+    }
+
+    template <int               direction,
+              bool              contract_over_rows,
+              bool              add,
+              bool              one_line = false,
+              EvaluatorQuantity quantity = EvaluatorQuantity::value>
+    void
     apply(const Number2 *DEAL_II_RESTRICT shape_data,
           const Number *                  in,
-          Number *                        out);
+          Number *                        out) const;
 
-  private:
-    const Number2 *shape_values;
-    const Number2 *shape_gradients;
-    const Number2 *shape_hessians;
+    template <int  face_direction,
+              bool contract_onto_face,
+              bool add,
+              int  max_derivative>
+    void
+    apply_face(const Number *DEAL_II_RESTRICT in,
+               Number *DEAL_II_RESTRICT       out) const;
+
+    const Number2 *    shape_values;
+    const Number2 *    shape_gradients;
+    const Number2 *    shape_hessians;
+    const unsigned int n_rows;
+    const unsigned int n_columns;
   };
 
 
 
-  template <int dim,
-            int n_rows,
-            int n_columns,
+  template <EvaluatorVariant variant,
+            int              dim,
             typename Number,
             typename Number2>
-  template <int  direction,
-            bool contract_over_rows,
-            bool add,
-            int  type,
-            bool one_line>
+  template <int               direction,
+            bool              contract_over_rows,
+            bool              add,
+            bool              one_line,
+            EvaluatorQuantity quantity>
   inline void
-  EvaluatorTensorProduct<evaluate_symmetric_hierarchical,
-                         dim,
-                         n_rows,
-                         n_columns,
-                         Number,
-                         Number2>::apply(const Number2 *DEAL_II_RESTRICT shapes,
-                                         const Number *                  in,
-                                         Number *                        out)
+  EvaluatorTensorProduct<variant, dim, 0, 0, Number, Number2>::apply(
+    const Number2 *DEAL_II_RESTRICT shape_data,
+    const Number *                  in,
+    Number *                        out) const
   {
     static_assert(one_line == false || direction == dim - 1,
                   "Single-line evaluation only works for direction=dim-1.");
-    static_assert(
-      type == 0 || type == 1,
-      "Only types 0 and 1 implemented for evaluate_symmetric_hierarchical.");
+    Assert(shape_data != nullptr,
+           ExcMessage(
+             "The given array shape_data must not be the null pointer!"));
     Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
              in != out,
            ExcMessage("In-place operation only supported for "
                       "n_rows==n_columns or single-line interpolation"));
-
-    // We cannot statically assert that direction is less than dim, so must do
-    // an additional dynamic check
     AssertIndexRange(direction, dim);
+    const int mm = contract_over_rows ? n_rows : n_columns,
+              nn = contract_over_rows ? n_columns : n_rows;
 
-    constexpr int nn     = contract_over_rows ? n_columns : n_rows;
-    constexpr int mm     = contract_over_rows ? n_rows : n_columns;
-    constexpr int n_cols = nn / 2;
-    constexpr int mid    = mm / 2;
+    const int stride =
+      direction == 0 ? 1 : Utilities::fixed_power<direction>(n_columns);
+    const int n_blocks1 = one_line ? 1 : stride;
+    const int n_blocks2 = direction >= dim - 1 ?
+                            1 :
+                            Utilities::fixed_power<dim - direction - 1>(n_rows);
+    Assert(n_rows <= 128, ExcNotImplemented());
 
-    constexpr int stride    = Utilities::pow(n_columns, direction);
-    constexpr int n_blocks1 = one_line ? 1 : stride;
-    constexpr int n_blocks2 =
-      Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
+    for (int i2 = 0; i2 < n_blocks2; ++i2)
+      {
+        for (int i1 = 0; i1 < n_blocks1; ++i1)
+          {
+            // the empty template case can only run the general evaluator or
+            // evenodd
+            constexpr EvaluatorVariant restricted_variant =
+              variant == evaluate_evenodd ? evaluate_evenodd : evaluate_general;
+            apply_matrix_vector_product<restricted_variant,
+                                        quantity,
+                                        contract_over_rows,
+                                        add>(
+              shape_data, in, out, n_rows, n_columns, stride, stride);
+
+            if (one_line == false)
+              {
+                ++in;
+                ++out;
+              }
+          }
+        if (one_line == false)
+          {
+            in += stride * (mm - 1);
+            out += stride * (nn - 1);
+          }
+      }
+  }
+
+
+
+  template <EvaluatorVariant variant,
+            int              dim,
+            typename Number,
+            typename Number2>
+  template <int  face_direction,
+            bool contract_onto_face,
+            bool add,
+            int  max_derivative>
+  inline void
+  EvaluatorTensorProduct<variant, dim, 0, 0, Number, Number2>::apply_face(
+    const Number *DEAL_II_RESTRICT in,
+    Number *DEAL_II_RESTRICT       out) const
+  {
+    Assert(shape_values != nullptr,
+           ExcMessage(
+             "The given array shape_data must not be the null pointer!"));
+    static_assert(dim > 0 && dim < 4, "Only dim=1,2,3 supported");
+    const int n_blocks1 = dim > 1 ? n_rows : 1;
+    const int n_blocks2 = dim > 2 ? n_rows : 1;
+
+    AssertIndexRange(face_direction, dim);
+    const int in_stride =
+      face_direction > 0 ? Utilities::fixed_power<face_direction>(n_rows) : 1;
+    const int out_stride =
+      dim > 1 ? Utilities::fixed_power<dim - 1>(n_rows) : 1;
 
-    // this code may look very inefficient at first sight due to the many
-    // different cases with if's at the innermost loop part, but all of the
-    // conditionals can be evaluated at compile time because they are
-    // templates, so the compiler should optimize everything away
     for (int i2 = 0; i2 < n_blocks2; ++i2)
       {
         for (int i1 = 0; i1 < n_blocks1; ++i1)
           {
-            if (contract_over_rows)
+            if (contract_onto_face == true)
               {
-                Number x[mm];
-                for (unsigned int i = 0; i < mm; ++i)
-                  x[i] = in[stride * i];
-                for (unsigned int col = 0; col < n_cols; ++col)
+                Number res0 = shape_values[0] * in[0];
+                Number res1, res2;
+                if (max_derivative > 0)
+                  res1 = shape_values[n_rows] * in[0];
+                if (max_derivative > 1)
+                  res2 = shape_values[2 * n_rows] * in[0];
+                for (unsigned int ind = 1; ind < n_rows; ++ind)
                   {
-                    Number r0, r1;
-                    if (mid > 0)
-                      {
-                        r0 = shapes[col] * x[0];
-                        r1 = shapes[col + n_columns] * x[1];
-                        for (unsigned int ind = 1; ind < mid; ++ind)
-                          {
-                            r0 +=
-                              shapes[col + 2 * ind * n_columns] * x[2 * ind];
-                            r1 += shapes[col + (2 * ind + 1) * n_columns] *
-                                  x[2 * ind + 1];
-                          }
-                      }
-                    else
-                      r0 = r1 = Number();
-                    if (mm % 2 == 1)
-                      r0 += shapes[col + (mm - 1) * n_columns] * x[mm - 1];
-                    if (add)
-                      {
-                        out[stride * col] += r0 + r1;
-                        if (type == 1)
-                          out[stride * (nn - 1 - col)] += r1 - r0;
-                        else
-                          out[stride * (nn - 1 - col)] += r0 - r1;
-                      }
-                    else
-                      {
-                        out[stride * col] = r0 + r1;
-                        if (type == 1)
-                          out[stride * (nn - 1 - col)] = r1 - r0;
-                        else
-                          out[stride * (nn - 1 - col)] = r0 - r1;
-                      }
+                    res0 += shape_values[ind] * in[in_stride * ind];
+                    if (max_derivative > 0)
+                      res1 += shape_values[ind + n_rows] * in[in_stride * ind];
+                    if (max_derivative > 1)
+                      res2 +=
+                        shape_values[ind + 2 * n_rows] * in[in_stride * ind];
+                  }
+                if (add)
+                  {
+                    out[0] += res0;
+                    if (max_derivative > 0)
+                      out[out_stride] += res1;
+                    if (max_derivative > 1)
+                      out[2 * out_stride] += res2;
                   }
-                if (nn % 2 == 1)
+                else
                   {
-                    Number             r0;
-                    const unsigned int shift = type == 1 ? 1 : 0;
-                    if (mid > 0)
-                      {
-                        r0 = shapes[n_cols + shift * n_columns] * x[shift];
-                        for (unsigned int ind = 1; ind < mid; ++ind)
-                          r0 += shapes[n_cols + (2 * ind + shift) * n_columns] *
-                                x[2 * ind + shift];
-                      }
-                    else
-                      r0 = 0;
-                    if (type != 1 && mm % 2 == 1)
-                      r0 += shapes[n_cols + (mm - 1) * n_columns] * x[mm - 1];
-                    if (add)
-                      out[stride * n_cols] += r0;
-                    else
-                      out[stride * n_cols] = r0;
+                    out[0] = res0;
+                    if (max_derivative > 0)
+                      out[out_stride] = res1;
+                    if (max_derivative > 1)
+                      out[2 * out_stride] = res2;
                   }
               }
             else
               {
-                Number xp[mid + 1], xm[mid > 0 ? mid : 1];
-                for (int i = 0; i < mid; ++i)
-                  if (type == 0)
-                    {
-                      xp[i] = in[stride * i] + in[stride * (mm - 1 - i)];
-                      xm[i] = in[stride * i] - in[stride * (mm - 1 - i)];
-                    }
-                  else
-                    {
-                      xp[i] = in[stride * i] - in[stride * (mm - 1 - i)];
-                      xm[i] = in[stride * i] + in[stride * (mm - 1 - i)];
-                    }
-                if (mm % 2 == 1)
-                  xp[mid] = in[stride * mid];
-                for (unsigned int col = 0; col < n_cols; ++col)
-                  {
-                    Number r0, r1;
-                    if (mid > 0)
-                      {
-                        r0 = shapes[2 * col * n_columns] * xp[0];
-                        r1 = shapes[(2 * col + 1) * n_columns] * xm[0];
-                        for (unsigned int ind = 1; ind < mid; ++ind)
-                          {
-                            r0 += shapes[2 * col * n_columns + ind] * xp[ind];
-                            r1 +=
-                              shapes[(2 * col + 1) * n_columns + ind] * xm[ind];
-                          }
-                      }
-                    else
-                      r0 = r1 = Number();
-                    if (mm % 2 == 1)
-                      {
-                        if (type == 1)
-                          r1 +=
-                            shapes[(2 * col + 1) * n_columns + mid] * xp[mid];
-                        else
-                          r0 += shapes[2 * col * n_columns + mid] * xp[mid];
-                      }
-                    if (add)
-                      {
-                        out[stride * (2 * col)] += r0;
-                        out[stride * (2 * col + 1)] += r1;
-                      }
-                    else
-                      {
-                        out[stride * (2 * col)]     = r0;
-                        out[stride * (2 * col + 1)] = r1;
-                      }
-                  }
-                if (nn % 2 == 1)
+                for (unsigned int col = 0; col < n_rows; ++col)
                   {
-                    Number r0;
-                    if (mid > 0)
-                      {
-                        r0 = shapes[(nn - 1) * n_columns] * xp[0];
-                        for (unsigned int ind = 1; ind < mid; ++ind)
-                          r0 += shapes[(nn - 1) * n_columns + ind] * xp[ind];
-                      }
-                    else
-                      r0 = Number();
-                    if (mm % 2 == 1 && type == 0)
-                      r0 += shapes[(nn - 1) * n_columns + mid] * xp[mid];
                     if (add)
-                      out[stride * (nn - 1)] += r0;
+                      out[col * in_stride] += shape_values[col] * in[0];
                     else
-                      out[stride * (nn - 1)] = r0;
+                      out[col * in_stride] = shape_values[col] * in[0];
+                    if (max_derivative > 0)
+                      out[col * in_stride] +=
+                        shape_values[col + n_rows] * in[out_stride];
+                    if (max_derivative > 1)
+                      out[col * in_stride] +=
+                        shape_values[col + 2 * n_rows] * in[2 * out_stride];
                   }
               }
-            if (one_line == false)
+
+            // increment: in regular case, just go to the next point in
+            // x-direction. If we are at the end of one chunk in x-dir, need
+            // to jump over to the next layer in z-direction
+            switch (face_direction)
               {
-                in += 1;
-                out += 1;
+                case 0:
+                  in += contract_onto_face ? n_rows : 1;
+                  out += contract_onto_face ? 1 : n_rows;
+                  break;
+                case 1:
+                  ++in;
+                  ++out;
+                  // faces 2 and 3 in 3d use local coordinate system zx, which
+                  // is the other way around compared to the tensor
+                  // product. Need to take that into account.
+                  if (dim == 3)
+                    {
+                      if (contract_onto_face)
+                        out += n_rows - 1;
+                      else
+                        in += n_rows - 1;
+                    }
+                  break;
+                case 2:
+                  ++in;
+                  ++out;
+                  break;
+                default:
+                  Assert(false, ExcNotImplemented());
               }
           }
-        if (one_line == false)
+        if (face_direction == 1 && dim == 3)
           {
-            in += stride * (mm - 1);
-            out += stride * (nn - 1);
+            // adjust for local coordinate system zx
+            if (contract_onto_face)
+              {
+                in += n_rows * (n_rows - 1);
+                out -= n_rows * n_rows - 1;
+              }
+            else
+              {
+                out += n_rows * (n_rows - 1);
+                in -= n_rows * n_rows - 1;
+              }
           }
       }
   }
 
 
 
+  /**
+   * Generic evaluator framework that valuates the given shape data in general
+   * dimensions using the tensor product form. Depending on the particular
+   * layout in the matrix entries, this corresponds to a usual matrix-matrix
+   * product or a matrix-matrix product including some symmetries. The actual
+   * work is implemented by functions of type apply_matrix_vector_product
+   * working on a single dimension, controlled by suitable strides, using the
+   * kernel specified via variant.
+   *
+   * @tparam variant Variant of evaluation used for creating template
+   *                 specializations
+   * @tparam dim Dimension of the function
+   * @tparam n_rows Number of rows in the transformation matrix, which corresponds
+   *                to the number of 1d shape functions in the usual tensor
+   *                contraction setting
+   * @tparam n_columns Number of columns in the transformation matrix, which
+   *                   corresponds to the number of 1d shape functions in the
+   *                   usual tensor contraction setting
+   * @tparam Number Abstract number type for input and output arrays
+   * @tparam Number2 Abstract number type for coefficient arrays (defaults to
+   *                 same type as the input/output arrays); must implement
+   *                 operator* with Number to be valid
+   * @tparam normal_dir Indicates the direction of the continuous component for the
+   *                    Raviart-Thomas space in terms of the normal onto the
+   * face, e.g 0 if the  is in x-direction, 1 if in y-direction, and 2 if in
+   * z-direction.
+   */
+  template <EvaluatorVariant variant,
+            int              dim,
+            int              n_rows,
+            int              n_columns,
+            int              normal_dir,
+            typename Number,
+            typename Number2 = Number>
+  struct EvaluatorTensorProductAnisotropic
+  {};
+
+
+
   /**
    * Internal evaluator for shape function in 2d and 3d using the
    * tensor product form of the anisotropic basis functions of the

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.