]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Provide function to evaluate polynomial on array of positions 13798/head
authorMartin Kronbichler <martin.kronbichler@uni-a.de>
Mon, 23 May 2022 18:45:20 +0000 (20:45 +0200)
committerMartin Kronbichler <martin.kronbichler@uni-a.de>
Wed, 25 May 2022 13:16:10 +0000 (15:16 +0200)
include/deal.II/base/polynomial.h
include/deal.II/matrix_free/tensor_product_kernels.h

index 24e9a634183b84e61533e916493ec50c84d464b6..984e120c15889ffbee0cc389199a632b7c5177ff 100644 (file)
@@ -24,6 +24,7 @@
 #include <deal.II/base/point.h>
 #include <deal.II/base/subscriptor.h>
 
+#include <array>
 #include <memory>
 #include <vector>
 
@@ -142,6 +143,24 @@ namespace Polynomials
           const unsigned int n_derivatives,
           Number2 *          values) const;
 
+    /**
+     * Similar to the function above, but evaluate the polynomials on several
+     * positions at once, as described by the array argument @p points. This
+     * function is can be faster than the other function when the same
+     * polynomial should be evaluated on several positions at once, e.g., the
+     * x,y,z coordinates of a point for tensor-product polynomials.
+     *
+     * The template type `Number2` must implement arithmetic
+     * operations such as additions or multiplication with the type
+     * `number` of the polynomial, and must be convertible from
+     * `number` by `operator=`.
+     */
+    template <std::size_t n_entries, typename Number2>
+    void
+    values_of_array(const std::array<Number2, n_entries> &points,
+                    const unsigned int                    n_derivatives,
+                    std::array<Number2, n_entries> *      values) const;
+
     /**
      * Degree of the polynomial. This is the degree reflected by the number of
      * coefficients provided by the constructor. Leading non-zero coefficients
@@ -827,6 +846,21 @@ namespace Polynomials
   Polynomial<number>::value(const Number2      x,
                             const unsigned int n_derivatives,
                             Number2 *          values) const
+  {
+    values_of_array(std::array<Number2, 1ul>{{x}},
+                    n_derivatives,
+                    reinterpret_cast<std::array<Number2, 1ul> *>(values));
+  }
+
+
+
+  template <typename number>
+  template <std::size_t n_entries, typename Number2>
+  inline void
+  Polynomial<number>::values_of_array(
+    const std::array<Number2, n_entries> &x,
+    const unsigned int                    n_derivatives,
+    std::array<Number2, n_entries> *      values) const
   {
     // evaluate Lagrange polynomial and derivatives
     if (in_lagrange_product_form == true)
@@ -839,12 +873,16 @@ namespace Polynomials
         switch (n_derivatives)
           {
             default:
-              values[0] = 1.;
-              for (unsigned int d = 1; d <= n_derivatives; ++d)
-                values[d] = 0.;
+              for (unsigned int e = 0; e < n_entries; ++e)
+                values[0][e] = 1.;
+              for (unsigned int k = 1; k <= n_derivatives; ++k)
+                for (unsigned int e = 0; e < n_entries; ++e)
+                  values[k][e] = 0.;
               for (unsigned int i = 0; i < n_supp; ++i)
                 {
-                  const Number2 v = x - lagrange_support_points[i];
+                  std::array<Number2, n_entries> v = x;
+                  for (unsigned int e = 0; e < n_entries; ++e)
+                    v[e] -= lagrange_support_points[i];
 
                   // multiply by (x-x_i) and compute action on all derivatives,
                   // too (inspired from automatic differentiation: implement the
@@ -853,8 +891,10 @@ namespace Polynomials
                   // value from the next lower derivative from the steps before,
                   // need to start from the highest derivative
                   for (unsigned int k = n_derivatives; k > 0; --k)
-                    values[k] = (values[k] * v + values[k - 1]);
-                  values[0] *= v;
+                    for (unsigned int e = 0; e < n_entries; ++e)
+                      values[k][e] = (values[k][e] * v[e] + values[k - 1][e]);
+                  for (unsigned int e = 0; e < n_entries; ++e)
+                    values[0][e] *= v[e];
                 }
               // finally, multiply by the weight in the Lagrange
               // denominator. Could be done instead of setting values[0] = 1
@@ -868,7 +908,8 @@ namespace Polynomials
                 number k_factorial = 1;
                 for (unsigned int k = 0; k <= n_derivatives; ++k)
                   {
-                    values[k] *= k_factorial * weight;
+                    for (unsigned int e = 0; e < n_entries; ++e)
+                      values[k][e] *= k_factorial * weight;
                     k_factorial *= static_cast<number>(k + 1);
                   }
               }
@@ -881,46 +922,60 @@ namespace Polynomials
             // compiler with the pointer aliasing analysis.
             case 0:
               {
-                Number2 value = 1.;
+                std::array<Number2, n_entries> value;
+                for (unsigned int e = 0; e < n_entries; ++e)
+                  value[e] = 1.;
                 for (unsigned int i = 0; i < n_supp; ++i)
-                  {
-                    const Number2 v = x - lagrange_support_points[i];
-                    value *= v;
-                  }
-                values[0] = weight * value;
+                  for (unsigned int e = 0; e < n_entries; ++e)
+                    value[e] *= (x[e] - lagrange_support_points[i]);
+
+                for (unsigned int e = 0; e < n_entries; ++e)
+                  values[0][e] = weight * value[e];
                 break;
               }
 
             case 1:
               {
-                Number2 value      = 1.;
-                Number2 derivative = 0.;
+                std::array<Number2, n_entries> value, derivative = {};
+                for (unsigned int e = 0; e < n_entries; ++e)
+                  value[e] = 1.;
                 for (unsigned int i = 0; i < n_supp; ++i)
+                  for (unsigned int e = 0; e < n_entries; ++e)
+                    {
+                      const Number2 v = x[e] - lagrange_support_points[i];
+                      derivative[e]   = derivative[e] * v + value[e];
+                      value[e] *= v;
+                    }
+
+                for (unsigned int e = 0; e < n_entries; ++e)
                   {
-                    const Number2 v = x - lagrange_support_points[i];
-                    derivative      = derivative * v + value;
-                    value *= v;
+                    values[0][e] = weight * value[e];
+                    values[1][e] = weight * derivative[e];
                   }
-                values[0] = weight * value;
-                values[1] = weight * derivative;
                 break;
               }
 
             case 2:
               {
-                Number2 value      = 1.;
-                Number2 derivative = 0.;
-                Number2 second     = 0.;
+                std::array<Number2, n_entries> value, derivative = {},
+                                                      second = {};
+                for (unsigned int e = 0; e < n_entries; ++e)
+                  value[e] = 1.;
                 for (unsigned int i = 0; i < n_supp; ++i)
+                  for (unsigned int e = 0; e < n_entries; ++e)
+                    {
+                      const Number2 v = x[e] - lagrange_support_points[i];
+                      second[e]       = second[e] * v + derivative[e];
+                      derivative[e]   = derivative[e] * v + value[e];
+                      value[e] *= v;
+                    }
+
+                for (unsigned int e = 0; e < n_entries; ++e)
                   {
-                    const Number2 v = x - lagrange_support_points[i];
-                    second          = second * v + derivative;
-                    derivative      = derivative * v + value;
-                    value *= v;
+                    values[0][e] = weight * value[e];
+                    values[1][e] = weight * derivative[e];
+                    values[2][e] = static_cast<number>(2) * weight * second[e];
                   }
-                values[0] = weight * value;
-                values[1] = weight * derivative;
-                values[2] = static_cast<number>(2) * weight * second;
                 break;
               }
           }
@@ -931,9 +986,12 @@ namespace Polynomials
 
     // if derivatives are needed, then do it properly by the full
     // Horner scheme
-    const unsigned int   m = coefficients.size();
-    std::vector<Number2> a(coefficients.size());
-    std::copy(coefficients.begin(), coefficients.end(), a.begin());
+    const unsigned int                          m = coefficients.size();
+    std::vector<std::array<Number2, n_entries>> a(coefficients.size());
+    for (unsigned int i = 0; i < coefficients.size(); ++i)
+      for (unsigned int e = 0; e < n_entries; ++e)
+        a[i][e] = coefficients[i];
+
     unsigned int j_factorial = 1;
 
     // loop over all requested derivatives. note that derivatives @p{j>m} are
@@ -943,15 +1001,18 @@ namespace Polynomials
     for (unsigned int j = 0; j < min_valuessize_m; ++j)
       {
         for (int k = m - 2; k >= static_cast<int>(j); --k)
-          a[k] += x * a[k + 1];
-        values[j] = static_cast<number>(j_factorial) * a[j];
+          for (unsigned int e = 0; e < n_entries; ++e)
+            a[k][e] += x[e] * a[k + 1][e];
+        for (unsigned int e = 0; e < n_entries; ++e)
+          values[j][e] = static_cast<number>(j_factorial) * a[j][e];
 
         j_factorial *= j + 1;
       }
 
     // fill higher derivatives by zero
     for (unsigned int j = min_valuessize_m; j <= n_derivatives; ++j)
-      values[j] = 0.;
+      for (unsigned int e = 0; e < n_entries; ++e)
+        values[j][e] = 0.;
   }
 
 
index 3f5233b42ab644c014207a7d06bd97ab71737518..50f31668005e027de460f673158f2e3845ac1132 100644 (file)
@@ -20,6 +20,7 @@
 #include <deal.II/base/config.h>
 
 #include <deal.II/base/aligned_vector.h>
+#include <deal.II/base/ndarray.h>
 #include <deal.II/base/polynomial.h>
 #include <deal.II/base/utilities.h>
 
@@ -2961,12 +2962,14 @@ namespace internal
       }
 
     AssertIndexRange(n_shapes, 200);
-    std::array<Number2, 2 * dim * 200> shapes;
+    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)
-      for (int i = 0; i < n_shapes; ++i)
-        poly[i].value(p[d], 1, shapes.data() + 2 * (d * n_shapes + i));
+      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
@@ -2984,22 +2987,22 @@ namespace internal
             if (renumber.empty())
               for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
                 {
-                  value += shapes[2 * i0] * values[i];
-                  deriv += shapes[2 * i0 + 1] * values[i];
+                  value += shapes[i0][0][0] * values[i];
+                  deriv += shapes[i0][1][0] * values[i];
                 }
             else
               for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
                 {
-                  value += shapes[2 * i0] * values[renumber[i]];
-                  deriv += shapes[2 * i0 + 1] * values[renumber[i]];
+                  value += shapes[i0][0][0] * values[renumber[i]];
+                  deriv += shapes[i0][1][0] * values[renumber[i]];
                 }
 
             // Interpolation + derivative in y direction
             if (dim > 1)
               {
-                value_y += value * shapes[2 * n_shapes + 2 * i1];
-                deriv_x += deriv * shapes[2 * n_shapes + 2 * i1];
-                deriv_y += value * shapes[2 * n_shapes + 2 * i1 + 1];
+                value_y += value * shapes[i1][0][1];
+                deriv_x += deriv * shapes[i1][0][1];
+                deriv_y += value * shapes[i1][1][1];
               }
             else
               {
@@ -3010,10 +3013,10 @@ namespace internal
         if (dim == 3)
           {
             // Interpolation + derivative in z direction
-            result.first += value_y * shapes[4 * n_shapes + 2 * i2];
-            result.second[0] += deriv_x * shapes[4 * n_shapes + 2 * i2];
-            result.second[1] += deriv_y * shapes[4 * n_shapes + 2 * i2];
-            result.second[2] += value_y * shapes[4 * n_shapes + 2 * i2 + 1];
+            result.first += value_y * shapes[i2][0][2];
+            result.second[0] += deriv_x * shapes[i2][0][2];
+            result.second[1] += deriv_y * shapes[i2][0][2];
+            result.second[2] += value_y * shapes[i2][1][2];
           }
         else if (dim == 2)
           {
@@ -3050,12 +3053,14 @@ namespace internal
            ExcDimensionMismatch(renumber.size(), values.size()));
 
     AssertIndexRange(n_shapes, 200);
-    std::array<Number2, 3 * dim * 200> shapes;
+    dealii::ndarray<Number2, 200, 3, dim> shapes;
 
     // Evaluate 1D polynomials and their derivatives
+    std::array<Number2, dim> point;
     for (unsigned int d = 0; d < dim; ++d)
-      for (int i = 0; i < n_shapes; ++i)
-        poly[i].value(p[d], 2, shapes.data() + 3 * (d * n_shapes + i));
+      point[d] = p[d];
+    for (int i = 0; i < n_shapes; ++i)
+      poly[i].values_of_array(point, 2, &shapes[i][0]);
 
     // Go through the tensor product of shape functions and interpolate
     // with optimal algorithm
@@ -3074,16 +3079,16 @@ namespace internal
             if (renumber.empty())
               for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
                 {
-                  value += shapes[3 * i0] * values[i];
-                  deriv_1 += shapes[3 * i0 + 1] * values[i];
-                  deriv_2 += shapes[3 * i0 + 2] * values[i];
+                  value += shapes[i0][0][0] * values[i];
+                  deriv_1 += shapes[i0][1][0] * values[i];
+                  deriv_2 += shapes[i0][2][0] * values[i];
                 }
             else
               for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
                 {
-                  value += shapes[3 * i0] * values[renumber[i]];
-                  deriv_1 += shapes[3 * i0 + 1] * values[renumber[i]];
-                  deriv_2 += shapes[3 * i0 + 2] * values[renumber[i]];
+                  value += shapes[i0][0][0] * values[renumber[i]];
+                  deriv_1 += shapes[i0][1][0] * values[renumber[i]];
+                  deriv_2 += shapes[i0][2][0] * values[renumber[i]];
                 }
 
             // Interpolation + derivative in y direction
@@ -3091,13 +3096,13 @@ namespace internal
               {
                 if (dim > 2)
                   {
-                    value_y += value * shapes[3 * n_shapes + 3 * i1];
-                    deriv_x += deriv_1 * shapes[3 * n_shapes + 3 * i1];
-                    deriv_y += value * shapes[3 * n_shapes + 3 * i1 + 1];
+                    value_y += value * shapes[i1][0][1];
+                    deriv_x += deriv_1 * shapes[i1][0][1];
+                    deriv_y += value * shapes[i1][1][1];
                   }
-                deriv_xx += deriv_2 * shapes[3 * n_shapes + 3 * i1];
-                deriv_xy += deriv_1 * shapes[3 * n_shapes + 3 * i1 + 1];
-                deriv_yy += value * shapes[3 * n_shapes + 3 * i1 + 2];
+                deriv_xx += deriv_2 * shapes[i1][0][1];
+                deriv_xy += deriv_1 * shapes[i1][1][1];
+                deriv_yy += value * shapes[i1][2][1];
               }
             else
               {
@@ -3107,12 +3112,12 @@ namespace internal
         if (dim == 3)
           {
             // Interpolation + derivative in z direction
-            result[0][0] += deriv_xx * shapes[6 * n_shapes + 3 * i2];
-            result[0][1] += deriv_xy * shapes[6 * n_shapes + 3 * i2];
-            result[0][2] += deriv_x * shapes[6 * n_shapes + 3 * i2 + 1];
-            result[1][1] += deriv_yy * shapes[6 * n_shapes + 3 * i2];
-            result[1][2] += deriv_y * shapes[6 * n_shapes + 3 * i2 + 1];
-            result[2][2] += value_y * shapes[6 * n_shapes + 3 * i2 + 2];
+            result[0][0] += deriv_xx * shapes[i2][0][2];
+            result[0][1] += deriv_xy * shapes[i2][0][2];
+            result[0][2] += deriv_x * shapes[i2][1][2];
+            result[1][1] += deriv_yy * shapes[i2][0][2];
+            result[1][2] += deriv_y * shapes[i2][1][2];
+            result[2][2] += value_y * shapes[i2][2][2];
           }
         else if (dim == 2)
           {
@@ -3149,42 +3154,43 @@ namespace internal
            ExcDimensionMismatch(renumber.size(), values.size()));
 
     AssertIndexRange(n_shapes, 200);
-    std::array<Number, 2 * dim * 200> shapes;
+    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)
-      for (int i = 0; i < n_shapes; ++i)
-        poly[i].value(p[d], 1, shapes.data() + 2 * (d * n_shapes + i));
+      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
     for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
       {
         const Number2 test_value_z =
-          dim > 2 ? (value * shapes[4 * n_shapes + 2 * i2] +
-                     gradient[2] * shapes[4 * n_shapes + 2 * i2 + 1]) :
-                    value;
+          dim > 2 ?
+            (value * shapes[i2][0][2] + gradient[2] * shapes[i2][1][2]) :
+            value;
         const Number2 test_grad_x =
-          dim > 2 ? gradient[0] * shapes[4 * n_shapes + 2 * i2] : gradient[0];
-        const Number2 test_grad_y =
-          dim > 2 ? gradient[1] * shapes[4 * n_shapes + 2 * i2] :
-                    (dim > 1 ? gradient[1] : Number2());
+          dim > 2 ? gradient[0] * shapes[i2][0][2] : gradient[0];
+        const Number2 test_grad_y = dim > 2 ?
+                                      gradient[1] * shapes[i2][0][2] :
+                                      (dim > 1 ? gradient[1] : Number2());
         for (int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1)
           {
-            const Number2 test_value_y =
-              dim > 1 ? (test_value_z * shapes[2 * n_shapes + 2 * i1] +
-                         test_grad_y * shapes[2 * n_shapes + 2 * i1 + 1]) :
-                        test_value_z;
+            const Number2 test_value_y = dim > 1 ?
+                                           (test_value_z * shapes[i1][0][1] +
+                                            test_grad_y * shapes[i1][1][1]) :
+                                           test_value_z;
             const Number2 test_grad_xy =
-              dim > 1 ? test_grad_x * shapes[2 * n_shapes + 2 * i1] :
-                        test_grad_x;
+              dim > 1 ? test_grad_x * shapes[i1][0][1] : test_grad_x;
             if (renumber.empty())
               for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
-                values[i] += shapes[2 * i0] * test_value_y +
-                             shapes[2 * i0 + 1] * test_grad_xy;
+                values[i] += shapes[i0][0][0] * test_value_y +
+                             shapes[i0][1][0] * test_grad_xy;
             else
               for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
-                values[renumber[i]] += shapes[2 * i0] * test_value_y +
-                                       shapes[2 * i0 + 1] * test_grad_xy;
+                values[renumber[i]] += shapes[i0][0][0] * test_value_y +
+                                       shapes[i0][1][0] * test_grad_xy;
           }
       }
   }

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.