]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Compute values_of_array only once
authorMaximilian Bergbauer <bergbauer@lnm.mw.tum.de>
Fri, 24 Mar 2023 22:06:41 +0000 (23:06 +0100)
committerMaximilian Bergbauer <bergbauer@lnm.mw.tum.de>
Thu, 30 Mar 2023 18:45:34 +0000 (20:45 +0200)
include/deal.II/matrix_free/fe_point_evaluation.h
include/deal.II/matrix_free/tensor_product_kernels.h
include/deal.II/non_matching/mapping_info.h

index 03e7592746a14f69bd2943d135ffa8250584de2f..40b65738b1486a50f575ed014f32d06fdc74bfcb 100644 (file)
@@ -649,6 +649,14 @@ private:
   void
   setup(const unsigned int first_selected_component);
 
+  /**
+   * Shared functionality of all @p reinit() functions. Resizes data fields and
+   * precomputes the @p shapes vector, holding the evaluation of 1D basis
+   * functions of tensor product polynomials, if necessary.
+   */
+  void
+  do_reinit();
+
   /**
    * Number of quadrature points of the current cell/face.
    */
@@ -775,6 +783,12 @@ private:
    * Bool indicating if class is reinitialized and data vectors a resized.
    */
   bool is_reinitialized;
+
+  /**
+   * Vector containing tensor product shape functions evaluated (during
+   * reinit()) at the vectorized unit points.
+   */
+  AlignedVector<dealii::ndarray<VectorizedArray<Number>, 2, dim>> shapes;
 };
 
 // ----------------------- template and inline function ----------------------
@@ -830,6 +844,8 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::setup(
   AssertIndexRange(first_selected_component + n_components,
                    fe->n_components() + 1);
 
+  shapes.reserve(100);
+
   bool         same_base_element   = true;
   unsigned int base_element_number = 0;
   component_in_base_element        = 0;
@@ -914,14 +930,7 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::reinit(
       fe_values->reinit(cell);
     }
 
-  const_cast<unsigned int &>(n_q_points) = unit_points.size();
-
-  if (update_flags & update_values)
-    values.resize(n_q_points, numbers::signaling_nan<value_type>());
-  if (update_flags & update_gradients)
-    gradients.resize(n_q_points, numbers::signaling_nan<gradient_type>());
-
-  is_reinitialized = true;
+  do_reinit();
 }
 
 
@@ -934,16 +943,7 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::reinit(
   current_cell_index  = cell_index;
   current_face_number = numbers::invalid_unsigned_int;
 
-  const_cast<unsigned int &>(n_q_points) =
-    mapping_info->get_unit_points(current_cell_index, current_face_number)
-      .size();
-
-  if (update_flags & update_values)
-    values.resize(n_q_points, numbers::signaling_nan<value_type>());
-  if (update_flags & update_gradients)
-    gradients.resize(n_q_points, numbers::signaling_nan<gradient_type>());
-
-  is_reinitialized = true;
+  do_reinit();
 }
 
 
@@ -957,15 +957,49 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::reinit(
   current_cell_index  = cell_index;
   current_face_number = face_number;
 
-  const_cast<unsigned int &>(n_q_points) =
-    mapping_info->get_unit_points(current_cell_index, current_face_number)
-      .size();
+  do_reinit();
+}
+
+
+
+template <int n_components, int dim, int spacedim, typename Number>
+void
+FEPointEvaluation<n_components, dim, spacedim, Number>::do_reinit()
+{
+  const auto unit_points =
+    mapping_info->get_unit_points(current_cell_index, current_face_number);
+
+  const_cast<unsigned int &>(n_q_points) = unit_points.size();
 
   if (update_flags & update_values)
     values.resize(n_q_points, numbers::signaling_nan<value_type>());
   if (update_flags & update_gradients)
     gradients.resize(n_q_points, numbers::signaling_nan<gradient_type>());
 
+  if (!polynomials_are_hat_functions)
+    {
+      const std::size_t n_points = unit_points.size();
+      const std::size_t n_lanes  = VectorizedArray<Number>::size();
+      const std::size_t n_batches =
+        n_points / n_lanes + (n_points % n_lanes > 0 ? 1 : 0);
+      const std::size_t n_shapes = poly.size();
+      shapes.resize_fast(n_batches * n_shapes);
+      for (unsigned int i = 0, qb = 0; i < n_points; i += n_lanes, ++qb)
+        {
+          // convert to vectorized format
+          Point<dim, VectorizedArray<Number>> vectorized_points;
+          for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
+            for (unsigned int d = 0; d < dim; ++d)
+              vectorized_points[d][j] = unit_points[i + j][d];
+
+          auto view =
+            make_array_view(shapes.begin() + qb * n_shapes,
+                            shapes.begin() + (qb * n_shapes + n_shapes));
+
+          internal::compute_values_of_array(view, poly, vectorized_points);
+        }
+    }
+
   is_reinitialized = true;
 }
 
@@ -977,6 +1011,9 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::evaluate(
   const ArrayView<const Number> &         solution_values,
   const EvaluationFlags::EvaluationFlags &evaluation_flag)
 {
+  if (!is_reinitialized)
+    reinit(numbers::invalid_unsigned_int);
+
   if (n_q_points == 0)
     return;
 
@@ -986,14 +1023,6 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::evaluate(
       fast_path)
     {
       // fast path with tensor product evaluation
-      const auto unit_points =
-        mapping_info->get_unit_points(current_cell_index, current_face_number);
-
-      // we need to call reinit() here if we reuse the same MappingInfo object
-      // for several FEPointEvaluation objects to resize the data fields
-      if (!is_reinitialized || n_q_points != unit_points.size())
-        reinit(numbers::invalid_unsigned_int);
-
       if (solution_renumbered.size() != dofs_per_component)
         solution_renumbered.resize(dofs_per_component);
       for (unsigned int comp = 0; comp < n_components; ++comp)
@@ -1011,23 +1040,41 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::evaluate(
       unit_gradients.resize(n_q_points,
                             numbers::signaling_nan<gradient_type>());
 
+      const auto unit_points =
+        mapping_info->get_unit_points(current_cell_index, current_face_number);
+      const auto &mapping_data =
+        mapping_info->get_mapping_data(current_cell_index, current_face_number);
+
       const std::size_t n_points = unit_points.size();
       const std::size_t n_lanes  = VectorizedArray<Number>::size();
-      for (unsigned int i = 0; i < n_points; i += n_lanes)
+      for (unsigned int i = 0, qb = 0; i < n_points; i += n_lanes, ++qb)
         {
-          // convert to vectorized format
-          Point<dim, VectorizedArray<Number>> vectorized_points;
-          for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
-            for (unsigned int d = 0; d < dim; ++d)
-              vectorized_points[d][j] = unit_points[i + j][d];
-
           // compute
-          const auto val_and_grad =
-            internal::evaluate_tensor_product_value_and_gradient(
-              poly,
-              solution_renumbered,
-              vectorized_points,
-              polynomials_are_hat_functions);
+          const unsigned int n_shapes     = poly.size();
+          const auto         val_and_grad = [&]() {
+            if (polynomials_are_hat_functions)
+              {
+                // convert to vectorized format
+                Point<dim, VectorizedArray<Number>> vectorized_points;
+                for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
+                  for (unsigned int d = 0; d < dim; ++d)
+                    vectorized_points[d][j] = unit_points[i + j][d];
+
+                return internal::
+                  evaluate_tensor_product_value_and_gradient_linear(
+                    poly, solution_renumbered, vectorized_points);
+              }
+            else
+              return internal::
+                evaluate_tensor_product_value_and_gradient_shapes<
+                  dim,
+                  value_type,
+                  VectorizedArray<Number>>(
+                  make_array_view(shapes.begin() + qb * n_shapes,
+                                  shapes.begin() + (qb * n_shapes + n_shapes)),
+                  poly.size(),
+                  solution_renumbered);
+          }();
 
           // convert back to standard format
           if (evaluation_flag & EvaluationFlags::values)
@@ -1048,9 +1095,6 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::evaluate(
                     Number>::set_gradient(val_and_grad.second,
                                           j,
                                           unit_gradients[i + j]);
-                  const auto &mapping_data =
-                    mapping_info->get_mapping_data(current_cell_index,
-                                                   current_face_number);
                   gradients[i + j] = apply_transformation(
                     mapping_data.inverse_jacobians[i + j].transpose(),
                     unit_gradients[i + j]);
@@ -1134,13 +1178,6 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::integrate(
       fast_path)
     {
       // fast path with tensor product integration
-      const auto unit_points =
-        mapping_info->get_unit_points(current_cell_index, current_face_number);
-
-      // we need to call reinit() here if we reuse the same MappingInfo object
-      // for several FEPointEvaluation objects to resize the data fields
-      if (!is_reinitialized || n_q_points != unit_points.size())
-        reinit(numbers::invalid_unsigned_int);
 
       if (integration_flags & EvaluationFlags::values)
         AssertIndexRange(n_q_points, values.size() + 1);
@@ -1156,16 +1193,15 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::integrate(
           n_components,
           VectorizedArray<Number>>::value_type());
 
+      const auto unit_points =
+        mapping_info->get_unit_points(current_cell_index, current_face_number);
+      const auto &mapping_data =
+        mapping_info->get_mapping_data(current_cell_index, current_face_number);
+
       const std::size_t n_points = unit_points.size();
       const std::size_t n_lanes  = VectorizedArray<Number>::size();
-      for (unsigned int i = 0; i < n_points; i += n_lanes)
+      for (unsigned int i = 0, qb = 0; i < n_points; i += n_lanes, ++qb)
         {
-          // convert to vectorized format
-          Point<dim, VectorizedArray<Number>> vectorized_points;
-          for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
-            for (unsigned int d = 0; d < dim; ++d)
-              vectorized_points[d][j] = unit_points[i + j][d];
-
           typename internal::ProductTypeNoPoint<value_type,
                                                 VectorizedArray<Number>>::type
             value = {};
@@ -1184,9 +1220,6 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::integrate(
           if (integration_flags & EvaluationFlags::gradients)
             for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
               {
-                const auto &mapping_data =
-                  mapping_info->get_mapping_data(current_cell_index,
-                                                 current_face_number);
                 gradients[i + j] =
                   apply_transformation(mapping_data.inverse_jacobians[i + j],
                                        gradients[i + j]);
@@ -1196,12 +1229,34 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::integrate(
               }
 
           // compute
-          internal::integrate_add_tensor_product_value_and_gradient(
-            poly,
-            value,
-            gradient,
-            vectorized_points,
-            solution_renumbered_vectorized);
+          const unsigned int n_shapes = poly.size();
+          if (polynomials_are_hat_functions)
+            {
+              // convert to vectorized format
+              Point<dim, VectorizedArray<Number>> vectorized_points;
+              for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
+                for (unsigned int d = 0; d < dim; ++d)
+                  vectorized_points[d][j] = unit_points[i + j][d];
+
+              internal::integrate_add_tensor_product_value_and_gradient_linear(
+                poly,
+                value,
+                gradient,
+                solution_renumbered_vectorized,
+                vectorized_points);
+            }
+          else
+            internal::integrate_add_tensor_product_value_and_gradient_shapes<
+              dim,
+              VectorizedArray<Number>,
+              typename internal::
+                ProductTypeNoPoint<value_type, VectorizedArray<Number>>::type>(
+              make_array_view(shapes.begin() + qb * n_shapes,
+                              shapes.begin() + (qb * n_shapes + n_shapes)),
+              n_shapes,
+              value,
+              gradient,
+              solution_renumbered_vectorized);
         }
 
       // add between the lanes and write into the result
index 18f0faddee36026876977a6169add760e5b96982..74b7279fe29f60456628d6a7bf0b9154b0443179 100644 (file)
@@ -2982,14 +2982,37 @@ namespace internal
 
 
 
+  /**
+   * Computes the values and derivatives of the 1d polynomials @p poly at the
+   * specified point @p p and stores it in @p shapes.
+   */
+  template <int dim, typename Number>
+  inline void
+  compute_values_of_array(
+    ArrayView<dealii::ndarray<Number, 2, dim>> &        shapes,
+    const std::vector<Polynomials::Polynomial<double>> &poly,
+    const Point<dim, Number> &                          p)
+  {
+    const int n_shapes = poly.size();
+
+    // Evaluate 1d polynomials and their derivatives
+    std::array<Number, dim> point;
+    for (unsigned int d = 0; d < dim; ++d)
+      point[d] = p[d];
+    for (int i = 0; i < n_shapes; ++i)
+      poly[i].values_of_array(point, 1, &shapes[i][0]);
+  }
+
+
+
   /**
    * Interpolate inner dimensions of tensor product shape functions.
    */
   template <int dim, int length, typename Number2, typename Number>
   inline std::array<typename ProductTypeNoPoint<Number, Number2>::type, 3>
-  do_interpolate_xy(const std::vector<Number> &                  values,
-                    const std::vector<unsigned int> &            renumber,
-                    const dealii::ndarray<Number2, 200, 2, dim> &shapes,
+  do_interpolate_xy(const std::vector<Number> &                        values,
+                    const std::vector<unsigned int> &                  renumber,
+                    const ArrayView<dealii::ndarray<Number2, 2, dim>> &shapes,
                     const int n_shapes_runtime,
                     int &     i)
   {
@@ -3035,122 +3058,27 @@ namespace internal
 
 
   /**
-   * Compute the polynomial interpolation of a tensor product shape function
-   * $\varphi_i$ given a vector of coefficients $u_i$ in the form
-   * $u_h(\mathbf{x}) = \sum_{i=1}^{k^d} \varphi_i(\mathbf{x}) u_i$. The shape
-   * functions $\varphi_i(\mathbf{x}) =
-   * \prod_{d=1}^{\text{dim}}\varphi_{i_d}^\text{1d}(x_d)$ represent a tensor
-   * product. The function returns a pair with the value of the interpolation
-   * as the first component and the gradient in reference coordinates as the
-   * second component. Note that for compound types (e.g. the `values` field
-   * begin a Point<spacedim> argument), the components of the gradient are
-   * sorted as Tensor<1, dim, Tensor<1, spacedim>> with the derivatives
-   * as the first index; this is a consequence of the generic arguments in the
-   * function.
-   *
-   * @param poly The underlying one-dimensional polynomial basis
-   * $\{\varphi^{1d}_{i_1}\}$ given as a vector of polynomials.
-   *
-   * @param values The expansion coefficients $u_i$ of type `Number` in
-   * the polynomial interpolation. The coefficients can be simply `double`
-   * variables but e.g. also Point<spacedim> in case they define arithmetic
-   * operations with the type `Number2`.
-   *
-   * @param p The position in reference coordinates where the interpolation
-   * should be evaluated.
-   *
-   * @param d_linear Flag to specify whether a d-linear (linear in 1d,
-   * bi-linear in 2d, tri-linear in 3d) interpolation should be made, which
-   * allows to unroll loops and considerably speed up evaluation.
-   *
-   * @param renumber Optional parameter to specify a renumbering in the
-   * coefficient vector, assuming that `values[renumber[i]]` returns
-   * the lexicographic (tensor product) entry of the coefficients. If the
-   * vector is entry, the values are assumed to be sorted lexicographically.
+   * Interpolates the values and gradients into the points specified in
+   * @p compute_values_of_array() with help of the precomputed @p shapes.
    */
   template <int dim, typename Number, typename Number2>
   inline std::pair<
     typename ProductTypeNoPoint<Number, Number2>::type,
     Tensor<1, dim, typename ProductTypeNoPoint<Number, Number2>::type>>
-  evaluate_tensor_product_value_and_gradient(
-    const std::vector<Polynomials::Polynomial<double>> &poly,
-    const std::vector<Number> &                         values,
-    const Point<dim, Number2> &                         p,
-    const bool                                          d_linear = false,
-    const std::vector<unsigned int> &                   renumber = {})
+  evaluate_tensor_product_value_and_gradient_shapes(
+    const ArrayView<dealii::ndarray<Number2, 2, dim>> &shapes,
+    const int                                          n_shapes,
+    const std::vector<Number> &                        values,
+    const std::vector<unsigned int> &                  renumber = {})
   {
     static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented");
 
     using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
 
-    // use `int` type for this variable and the loops below to inform the
-    // compiler that the loops below will never overflow, which allows it to
-    // generate more optimized code for the variable loop bounds in the
-    // present context
-    const int n_shapes = poly.size();
     AssertDimension(Utilities::pow(n_shapes, dim), values.size());
     Assert(renumber.empty() || renumber.size() == values.size(),
            ExcDimensionMismatch(renumber.size(), values.size()));
 
-    // shortcut for linear interpolation to speed up evaluation
-    if (d_linear)
-      {
-        AssertDimension(poly.size(), 2);
-        for (unsigned int i = 0; i < renumber.size(); ++i)
-          AssertDimension(renumber[i], i);
-
-        if (dim == 1)
-          {
-            Tensor<1, dim, Number3> derivative;
-            derivative[0] = values[1] - values[0];
-            return std::make_pair((1. - p[0]) * values[0] + p[0] * values[1],
-                                  derivative);
-          }
-        else if (dim == 2)
-          {
-            const Number2           x0 = 1. - p[0], x1 = p[0];
-            const Number3           tmp0   = x0 * values[0] + x1 * values[1];
-            const Number3           tmp1   = x0 * values[2] + x1 * values[3];
-            const Number3           mapped = (1. - p[1]) * tmp0 + p[1] * tmp1;
-            Tensor<1, dim, Number3> derivative;
-            derivative[0] = (1. - p[1]) * (values[1] - values[0]) +
-                            p[1] * (values[3] - values[2]);
-            derivative[1] = tmp1 - tmp0;
-            return std::make_pair(mapped, derivative);
-          }
-        else if (dim == 3)
-          {
-            const Number2 x0 = 1. - p[0], x1 = p[0], y0 = 1. - p[1], y1 = p[1],
-                          z0 = 1. - p[2], z1 = p[2];
-            const Number3           tmp0   = x0 * values[0] + x1 * values[1];
-            const Number3           tmp1   = x0 * values[2] + x1 * values[3];
-            const Number3           tmpy0  = y0 * tmp0 + y1 * tmp1;
-            const Number3           tmp2   = x0 * values[4] + x1 * values[5];
-            const Number3           tmp3   = x0 * values[6] + x1 * values[7];
-            const Number3           tmpy1  = y0 * tmp2 + y1 * tmp3;
-            const Number3           mapped = z0 * tmpy0 + z1 * tmpy1;
-            Tensor<1, dim, Number3> derivative;
-            derivative[2] = tmpy1 - tmpy0;
-            derivative[1] = z0 * (tmp1 - tmp0) + z1 * (tmp3 - tmp2);
-            derivative[0] =
-              z0 *
-                (y0 * (values[1] - values[0]) + y1 * (values[3] - values[2])) +
-              z1 *
-                (y0 * (values[5] - values[4]) + y1 * (values[7] - values[6]));
-            return std::make_pair(mapped, derivative);
-          }
-      }
-
-    AssertIndexRange(n_shapes, 200);
-    dealii::ndarray<Number2, 200, 2, dim> shapes;
-
-    // Evaluate 1d polynomials and their derivatives
-    std::array<Number2, dim> point;
-    for (unsigned int d = 0; d < dim; ++d)
-      point[d] = p[d];
-    for (int i = 0; i < n_shapes; ++i)
-      poly[i].values_of_array(point, 1, &shapes[i][0]);
-
     // Go through the tensor product of shape functions and interpolate
     // with optimal algorithm
     std::pair<Number3, Tensor<1, dim, Number3>> result = {};
@@ -3199,6 +3127,147 @@ namespace internal
 
 
 
+  /**
+   * Specializes @p evaluate_tensor_product_value_and_gradient() for linear
+   * polynomials which massively reduces the necessary instructions.
+   */
+  template <int dim, typename Number, typename Number2>
+  inline std::pair<
+    typename ProductTypeNoPoint<Number, Number2>::type,
+    Tensor<1, dim, typename ProductTypeNoPoint<Number, Number2>::type>>
+  evaluate_tensor_product_value_and_gradient_linear(
+    const std::vector<Polynomials::Polynomial<double>> &poly,
+    const std::vector<Number> &                         values,
+    const Point<dim, Number2> &                         p,
+    const std::vector<unsigned int> &                   renumber = {})
+  {
+    (void)poly;
+    static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented");
+
+    using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
+
+    AssertDimension(Utilities::pow(poly.size(), dim), values.size());
+    Assert(renumber.empty() || renumber.size() == values.size(),
+           ExcDimensionMismatch(renumber.size(), values.size()));
+
+    AssertDimension(poly.size(), 2);
+    for (unsigned int i = 0; i < renumber.size(); ++i)
+      AssertDimension(renumber[i], i);
+
+    if (dim == 1)
+      {
+        Tensor<1, dim, Number3> derivative;
+        derivative[0] = values[1] - values[0];
+        return std::make_pair((1. - p[0]) * values[0] + p[0] * values[1],
+                              derivative);
+      }
+    else if (dim == 2)
+      {
+        const Number2           x0 = 1. - p[0], x1 = p[0];
+        const Number3           tmp0   = x0 * values[0] + x1 * values[1];
+        const Number3           tmp1   = x0 * values[2] + x1 * values[3];
+        const Number3           mapped = (1. - p[1]) * tmp0 + p[1] * tmp1;
+        Tensor<1, dim, Number3> derivative;
+        derivative[0] = (1. - p[1]) * (values[1] - values[0]) +
+                        p[1] * (values[3] - values[2]);
+        derivative[1] = tmp1 - tmp0;
+        return std::make_pair(mapped, derivative);
+      }
+    else if (dim == 3)
+      {
+        const Number2 x0 = 1. - p[0], x1 = p[0], y0 = 1. - p[1], y1 = p[1],
+                      z0 = 1. - p[2], z1 = p[2];
+        const Number3           tmp0   = x0 * values[0] + x1 * values[1];
+        const Number3           tmp1   = x0 * values[2] + x1 * values[3];
+        const Number3           tmpy0  = y0 * tmp0 + y1 * tmp1;
+        const Number3           tmp2   = x0 * values[4] + x1 * values[5];
+        const Number3           tmp3   = x0 * values[6] + x1 * values[7];
+        const Number3           tmpy1  = y0 * tmp2 + y1 * tmp3;
+        const Number3           mapped = z0 * tmpy0 + z1 * tmpy1;
+        Tensor<1, dim, Number3> derivative;
+        derivative[2] = tmpy1 - tmpy0;
+        derivative[1] = z0 * (tmp1 - tmp0) + z1 * (tmp3 - tmp2);
+        derivative[0] =
+          z0 * (y0 * (values[1] - values[0]) + y1 * (values[3] - values[2])) +
+          z1 * (y0 * (values[5] - values[4]) + y1 * (values[7] - values[6]));
+        return std::make_pair(mapped, derivative);
+      }
+
+    // work around a compile error: missing return statement
+    return std::make_pair(Number3(), Tensor<1, dim, Number3>());
+  }
+
+
+
+  /**
+   * Compute the polynomial interpolation of a tensor product shape function
+   * $\varphi_i$ given a vector of coefficients $u_i$ in the form
+   * $u_h(\mathbf{x}) = \sum_{i=1}^{k^d} \varphi_i(\mathbf{x}) u_i$. The shape
+   * functions $\varphi_i(\mathbf{x}) =
+   * \prod_{d=1}^{\text{dim}}\varphi_{i_d}^\text{1d}(x_d)$ represent a tensor
+   * product. The function returns a pair with the value of the interpolation
+   * as the first component and the gradient in reference coordinates as the
+   * second component. Note that for compound types (e.g. the `values` field
+   * begin a Point<spacedim> argument), the components of the gradient are
+   * sorted as Tensor<1, dim, Tensor<1, spacedim>> with the derivatives
+   * as the first index; this is a consequence of the generic arguments in the
+   * function.
+   *
+   * @param poly The underlying one-dimensional polynomial basis
+   * $\{\varphi^{1d}_{i_1}\}$ given as a vector of polynomials.
+   *
+   * @param values The expansion coefficients $u_i$ of type `Number` in
+   * the polynomial interpolation. The coefficients can be simply `double`
+   * variables but e.g. also Point<spacedim> in case they define arithmetic
+   * operations with the type `Number2`.
+   *
+   * @param p The position in reference coordinates where the interpolation
+   * should be evaluated.
+   *
+   * @param d_linear Flag to specify whether a d-linear (linear in 1d,
+   * bi-linear in 2d, tri-linear in 3d) interpolation should be made, which
+   * allows to unroll loops and considerably speed up evaluation.
+   *
+   * @param renumber Optional parameter to specify a renumbering in the
+   * coefficient vector, assuming that `values[renumber[i]]` returns
+   * the lexicographic (tensor product) entry of the coefficients. If the
+   * vector is entry, the values are assumed to be sorted lexicographically.
+   */
+  template <int dim, typename Number, typename Number2>
+  inline std::pair<
+    typename ProductTypeNoPoint<Number, Number2>::type,
+    Tensor<1, dim, typename ProductTypeNoPoint<Number, Number2>::type>>
+  evaluate_tensor_product_value_and_gradient(
+    const std::vector<Polynomials::Polynomial<double>> &poly,
+    const std::vector<Number> &                         values,
+    const Point<dim, Number2> &                         p,
+    const bool                                          d_linear = false,
+    const std::vector<unsigned int> &                   renumber = {})
+  {
+    if (d_linear)
+      {
+        return evaluate_tensor_product_value_and_gradient_linear(poly,
+                                                                 values,
+                                                                 p,
+                                                                 renumber);
+      }
+    else
+      {
+        std::array<dealii::ndarray<Number2, 2, dim>, 200> shapes;
+
+        auto view = make_array_view(shapes);
+
+        compute_values_of_array(view, poly, p);
+
+        return evaluate_tensor_product_value_and_gradient_shapes<dim,
+                                                                 Number,
+                                                                 Number2>(
+          view, poly.size(), values, renumber);
+      }
+  }
+
+
+
   template <int dim, typename Number, typename Number2>
   SymmetricTensor<2, dim, typename ProductTypeNoPoint<Number, Number2>::type>
   evaluate_tensor_product_hessian(
@@ -3305,12 +3374,13 @@ namespace internal
    */
   template <int dim, int length, typename Number2, typename Number>
   inline void
-  do_apply_test_functions_xy(AlignedVector<Number2> &         values,
-                             const std::vector<unsigned int> &renumber,
-                             const dealii::ndarray<Number, 200, 2, dim> &shapes,
-                             const std::array<Number2, 3> &test_grads_value,
-                             const int                     n_shapes_runtime,
-                             int &                         i)
+  do_apply_test_functions_xy(
+    AlignedVector<Number2> &                          values,
+    const std::vector<unsigned int> &                 renumber,
+    const ArrayView<dealii::ndarray<Number, 2, dim>> &shapes,
+    const std::array<Number2, 3> &                    test_grads_value,
+    const int                                         n_shapes_runtime,
+    int &                                             i)
   {
     const int n_shapes = length > 0 ? length : n_shapes_runtime;
     for (int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1)
@@ -3344,32 +3414,21 @@ namespace internal
    */
   template <int dim, typename Number, typename Number2>
   inline void
-  integrate_add_tensor_product_value_and_gradient(
-    const std::vector<Polynomials::Polynomial<double>> &poly,
-    const Number2 &                                     value,
-    const Tensor<1, dim, Number2> &                     gradient,
-    const Point<dim, Number> &                          p,
-    AlignedVector<Number2> &                            values,
-    const std::vector<unsigned int> &                   renumber = {})
+  integrate_add_tensor_product_value_and_gradient_shapes(
+    const ArrayView<dealii::ndarray<Number, 2, dim>> &shapes,
+    const int                                         n_shapes,
+    const Number2 &                                   value,
+    const Tensor<1, dim, Number2> &                   gradient,
+    AlignedVector<Number2> &                          values,
+    const std::vector<unsigned int> &                 renumber = {})
   {
     static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented");
 
     // as in evaluate, use `int` type to produce better code in this context
-    const int n_shapes = poly.size();
     AssertDimension(Utilities::pow(n_shapes, dim), values.size());
     Assert(renumber.empty() || renumber.size() == values.size(),
            ExcDimensionMismatch(renumber.size(), values.size()));
 
-    AssertIndexRange(n_shapes, 200);
-    dealii::ndarray<Number, 200, 2, dim> shapes;
-
-    // Evaluate 1d polynomials and their derivatives
-    std::array<Number, dim> point;
-    for (unsigned int d = 0; d < dim; ++d)
-      point[d] = p[d];
-    for (int i = 0; i < n_shapes; ++i)
-      poly[i].values_of_array(point, 1, &shapes[i][0]);
-
     // Implement the transpose of the function above
     std::array<Number2, 3> test_grads_value;
     for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
@@ -3411,6 +3470,86 @@ namespace internal
 
 
 
+  /**
+   * Specializes @p evaluate_tensor_product_value_and_gradient() for linear
+   * polynomials which massively reduces the necessary instructions.
+   */
+  template <int dim, typename Number, typename Number2>
+  inline void
+  integrate_add_tensor_product_value_and_gradient_linear(
+    const std::vector<Polynomials::Polynomial<double>> &poly,
+    const Number2 &                                     value,
+    const Tensor<1, dim, Number2> &                     gradient,
+    AlignedVector<Number2> &                            values,
+    const Point<dim, Number> &                          p,
+    const std::vector<unsigned int> &                   renumber = {})
+  {
+    (void)poly;
+    static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented");
+
+    AssertDimension(Utilities::pow(poly.size(), dim), values.size());
+    Assert(renumber.empty() || renumber.size() == values.size(),
+           ExcDimensionMismatch(renumber.size(), values.size()));
+
+    AssertDimension(poly.size(), 2);
+    for (unsigned int i = 0; i < renumber.size(); ++i)
+      AssertDimension(renumber[i], i);
+
+    if (dim == 1)
+      {
+        const auto x0 = 1. - p[0], x1 = p[0];
+
+        values[0] = value * x0 - gradient[0];
+        values[1] = value * x1 + gradient[0];
+      }
+    else if (dim == 2)
+      {
+        const auto x0 = 1. - p[0], x1 = p[0], y0 = 1. - p[1], y1 = p[1];
+
+        const auto test_value_y0 = value * y0 - gradient[1];
+        const auto test_grad_xy0 = gradient[0] * y0;
+        const auto test_value_y1 = value * y1 + gradient[1];
+        const auto test_grad_xy1 = gradient[0] * y1;
+
+        values[0] += x0 * test_value_y0 - test_grad_xy0;
+        values[1] += x1 * test_value_y0 + test_grad_xy0;
+        values[2] += x0 * test_value_y1 - test_grad_xy1;
+        values[3] += x1 * test_value_y1 + test_grad_xy1;
+      }
+    else if (dim == 3)
+      {
+        const auto x0 = 1. - p[0], x1 = p[0], y0 = 1. - p[1], y1 = p[1],
+                   z0 = 1. - p[2], z1 = p[2];
+
+        const auto test_value_z0 = value * z0 - gradient[2];
+        const auto test_grad_x0  = gradient[0] * z0;
+        const auto test_grad_y0  = gradient[1] * z0;
+        const auto test_value_z1 = value * z1 + gradient[2];
+        const auto test_grad_x1  = gradient[0] * z1;
+        const auto test_grad_y1  = gradient[1] * z1;
+
+        const auto test_value_y00 = test_value_z0 * y0 - test_grad_y0;
+        const auto test_grad_xy00 = test_grad_x0 * y0;
+        const auto test_value_y01 = test_value_z0 * y1 + test_grad_y0;
+        const auto test_grad_xy01 = test_grad_x0 * y1;
+        const auto test_value_y10 = test_value_z1 * y0 - test_grad_y1;
+        const auto test_grad_xy10 = test_grad_x1 * y0;
+        const auto test_value_y11 = test_value_z1 * y1 + test_grad_y1;
+        const auto test_grad_xy11 = test_grad_x1 * y1;
+
+        values[0] += x0 * test_value_y00 - test_grad_xy00;
+        values[1] += x1 * test_value_y00 + test_grad_xy00;
+        values[2] += x0 * test_value_y01 - test_grad_xy01;
+        values[3] += x1 * test_value_y01 + test_grad_xy01;
+        values[4] += x0 * test_value_y10 - test_grad_xy10;
+        values[5] += x1 * test_value_y10 + test_grad_xy10;
+        values[6] += x0 * test_value_y11 - test_grad_xy11;
+        values[7] += x1 * test_value_y11 + test_grad_xy11;
+      }
+  }
+
+
+
   template <int dim, int n_points_1d_template, typename Number>
   inline void
   weight_fe_q_dofs_by_entity(const Number *     weights,
index 21347f05303a723e89f7bfb8e181139f78c44860..63ad40f87d07b7c4a68d2c367d6e4d184d22e8a8 100644 (file)
@@ -616,7 +616,7 @@ namespace NonMatching
 
 
   template <int dim, int spacedim>
-  const ArrayView<const Point<dim>>
+  inline const ArrayView<const Point<dim>>
   MappingInfo<dim, spacedim>::get_unit_points(
     const unsigned int cell_index,
     const unsigned int face_number) const
@@ -703,7 +703,7 @@ namespace NonMatching
 
 
   template <int dim, int spacedim>
-  const typename MappingInfo<dim, spacedim>::MappingData &
+  inline const typename MappingInfo<dim, spacedim>::MappingData &
   MappingInfo<dim, spacedim>::get_mapping_data(
     const unsigned int cell_index,
     const unsigned int face_number) const

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.