From: Rene Gassmoeller <rene.gassmoeller@mailbox.org>
Date: Thu, 15 Sep 2016 01:19:18 +0000 (-0600)
Subject: Provide shortcut and allocate on stack
X-Git-Tag: v8.5.0-rc1~577^2~1
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=42956ac212c15d80de460db322e409065bd92401;p=dealii.git

Provide shortcut and allocate on stack
---

diff --git a/include/deal.II/base/polynomial.h b/include/deal.II/base/polynomial.h
index 01888a2292..22c60b963b 100644
--- a/include/deal.II/base/polynomial.h
+++ b/include/deal.II/base/polynomial.h
@@ -103,6 +103,20 @@ namespace Polynomials
     void value (const number         x,
                 std::vector<number> &values) const;
 
+    /**
+     * Return the values and the derivatives of the Polynomial at point
+     * <tt>x</tt>.  <tt>values[i], i=0,...,values_size-1</tt> includes the
+     * <tt>i</tt>th derivative. The number of derivatives to be computed is
+     * determined by @p values_size and @p values has to provide sufficient
+     * space for @p values_size values.
+     *
+     * This function uses the Horner scheme for numerical stability of the
+     * evaluation.
+     */
+    void value (const number         x,
+                const unsigned int values_size,
+                number *values) const;
+
     /**
      * Degree of the polynomial. This is the degree reflected by the number of
      * coefficients provided by the constructor. Leading non-zero coefficients
diff --git a/include/deal.II/base/polynomials_piecewise.h b/include/deal.II/base/polynomials_piecewise.h
index 3746c1e48d..e070407690 100644
--- a/include/deal.II/base/polynomials_piecewise.h
+++ b/include/deal.II/base/polynomials_piecewise.h
@@ -83,7 +83,7 @@ namespace Polynomials
      * Return the values and the derivatives of the Polynomial at point
      * <tt>x</tt>.  <tt>values[i], i=0,...,values.size()-1</tt> includes the
      * <tt>i</tt>th derivative. The number of derivatives to be computed is
-     * thus determined by the size of the array passed.
+     * thus determined by the size of the vector passed.
      *
      * Note that all the derivatives evaluate to zero at the border between
      * intervals (assuming exact arithmetics) in the interior of the unit
@@ -96,6 +96,25 @@ namespace Polynomials
     void value (const number         x,
                 std::vector<number> &values) const;
 
+    /**
+     * Return the values and the derivatives of the Polynomial at point
+     * <tt>x</tt>.  <tt>values[i], i=0,...,values_size-1</tt> includes the
+     * <tt>i</tt>th derivative.The number of derivatives to be computed is
+     * determined by @p values_size and @p values has to provide sufficient
+     * space for @p values_size values.
+     *
+     * Note that all the derivatives evaluate to zero at the border between
+     * intervals (assuming exact arithmetics) in the interior of the unit
+     * interval, as there is no unique gradient value in that case for a
+     * piecewise polynomial. This is not always desired (e.g., when evaluating
+     * jumps of gradients on the element boundary), but it is the user's
+     * responsibility to avoid evaluation at these points when it does not
+     * make sense.
+     */
+    void value (const number         x,
+                const unsigned int values_size,
+                number *values) const;
+
     /**
      * Degree of the polynomial. This is the degree of the underlying base
      * polynomial.
diff --git a/source/base/polynomial.cc b/source/base/polynomial.cc
index cdea89d303..32b03950eb 100644
--- a/source/base/polynomial.cc
+++ b/source/base/polynomial.cc
@@ -102,7 +102,19 @@ namespace Polynomials
                              std::vector<number> &values) const
   {
     Assert (values.size() > 0, ExcZero());
-    const unsigned int values_size=values.size();
+
+    value(x,values.size(),&values[0]);
+  }
+
+
+
+  template <typename number>
+  void
+  Polynomial<number>::value (const number         x,
+                             const unsigned int values_size,
+                             number *values) const
+  {
+    Assert(values_size > 0, ExcZero());
 
     // evaluate Lagrange polynomial and derivatives
     if (in_lagrange_product_form == true)
diff --git a/source/base/polynomials_piecewise.cc b/source/base/polynomials_piecewise.cc
index 25656864f9..614b4b3d85 100644
--- a/source/base/polynomials_piecewise.cc
+++ b/source/base/polynomials_piecewise.cc
@@ -46,7 +46,19 @@ namespace Polynomials
                                       std::vector<number> &values) const
   {
     Assert (values.size() > 0, ExcZero());
-    const unsigned int values_size=values.size();
+
+    value(x,values.size(),&values[0]);
+  }
+
+
+
+  template <typename number>
+  void
+  PiecewisePolynomial<number>::value (const number         x,
+                                      const unsigned int values_size,
+                                      number *values) const
+  {
+    Assert (values_size > 0, ExcZero());
 
     // shift polynomial if necessary
     number y = x;
@@ -60,7 +72,7 @@ namespace Polynomials
             const double offset = step * interval;
             if (x<offset || x>offset+step+step)
               {
-                for (unsigned int k=0; k<values.size(); ++k)
+                for (unsigned int k=0; k<values_size; ++k)
                   values[k] = 0;
                 return;
               }
@@ -77,7 +89,7 @@ namespace Polynomials
             const double offset = step * interval;
             if (x<offset || x>offset+step)
               {
-                for (unsigned int k=0; k<values.size(); ++k)
+                for (unsigned int k=0; k<values_size; ++k)
                   values[k] = 0;
                 return;
               }
@@ -100,7 +112,7 @@ namespace Polynomials
           }
       }
 
-    polynomial.value(y, values);
+    polynomial.value(y, values_size, values);
 
     // change sign if necessary
     for (unsigned int j=1; j<values_size; j+=2)
diff --git a/source/base/tensor_product_polynomials.cc b/source/base/tensor_product_polynomials.cc
index 9ee229c6b2..fa382021db 100644
--- a/source/base/tensor_product_polynomials.cc
+++ b/source/base/tensor_product_polynomials.cc
@@ -18,6 +18,8 @@
 #include <deal.II/base/exceptions.h>
 #include <deal.II/base/table.h>
 
+#include <deal.II/base/std_cxx11/array.h>
+
 DEAL_II_NAMESPACE_OPEN
 
 
@@ -254,7 +256,7 @@ compute (const Point<dim>            &p,
   Assert (fourth_derivatives.size()==n_tensor_pols|| fourth_derivatives.size()==0,
           ExcDimensionMismatch2(fourth_derivatives.size(), n_tensor_pols, 0));
 
-  const bool update_values     = (values.size() == n_tensor_pols),
+  const bool update_values     = (values.size()==n_tensor_pols),
              update_grads      = (grads.size()==n_tensor_pols),
              update_grad_grads = (grad_grads.size()==n_tensor_pols),
              update_3rd_derivatives      = (third_derivatives.size()==n_tensor_pols),
@@ -275,26 +277,68 @@ compute (const Point<dim>            &p,
   if (update_4th_derivatives)
     n_values_and_derivatives = 5;
 
+  // Provide a shortcut if only values are requested. For this case usually the
+  // temporary memory allocation below does not pay off. Also we loop over all
+  // tensor polynomials in a peculiar way to avoid all dynamic memory allocation.
+  // We need to evaluate the polynomial value n_polynomials*dim times, and
+  // need to compute n_polynomials^dim tensor polynomials. Therefore we can split
+  // the loop over the tensor polynomials into one loop over n_polynomials,
+  // evaluate this polynomial for each dimension and multiply it with the
+  // associated n_polynomials^(dim-1) tensor polynomials before moving to the
+  // next polynomial.
+  if (n_values_and_derivatives == 1)
+    {
+      const unsigned int n_polynomials = polynomials.size();
+      const unsigned int factor = n_tensor_pols/n_polynomials;
+
+      std::fill(values.begin(),values.end(),1.0);
+
+      for (unsigned int polynomial=0; polynomial<n_polynomials; ++polynomial)
+        {
+          for (unsigned int d=0; d<dim; ++d)
+            {
+              const double value = polynomials[polynomial].value(p(d));
+
+              for (unsigned int f=0; f<factor; ++f)
+                {
+                  unsigned int index = 0;
+                  switch (d)
+                    {
+                    case 0:
+                      index = polynomial+f*n_polynomials;
+                      break;
+                    case 1:
+                      index = polynomial*n_polynomials + (f/n_polynomials)*n_polynomials*n_polynomials + f%n_polynomials;
+                      break;
+                    case 2:
+                      index = polynomial*factor + f;
+                      break;
+                    }
+
+                  const unsigned int i = index_map_inverse[index];
+                  values[i] *= value;
+                }
+            }
+        }
+      return;
+    }
+
 
-  // compute the values (and derivatives, if
+  // Compute the values (and derivatives, if
   // necessary) of all polynomials at this
-  // evaluation point. to avoid many
-  // reallocation, use one std::vector for
-  // polynomial evaluation and store the
-  // result as Tensor<1,5> (that has enough
+  // evaluation point. To avoid expensive memory allocations,
+  // use alloca to allocate a (small) amount of memory
+  // on the stack and store the
+  // result in a vector of arrays (that has enough
   // fields for any evaluation of values and
-  // derivatives, up to the 4th derivative)
-  Table<2,Tensor<1,5> > v(dim, polynomials.size());
-  {
-    std::vector<double> tmp (n_values_and_derivatives);
+  // derivatives, up to the 4th derivative).
+  std_cxx11::array<std_cxx11::array<double, 5>, dim> *v =
+    (std_cxx11::array<std_cxx11::array<double, 5>, dim> *)
+    alloca(sizeof(std_cxx11::array<std_cxx11::array<double, 5>, dim>) * polynomials.size());
+
+  for (unsigned int i=0; i<polynomials.size(); ++i)
     for (unsigned int d=0; d<dim; ++d)
-      for (unsigned int i=0; i<polynomials.size(); ++i)
-        {
-          polynomials[i].value(p(d), tmp);
-          for (unsigned int e=0; e<n_values_and_derivatives; ++e)
-            v(d,i)[e] = tmp[e];
-        };
-  }
+      polynomials[i].value(p(d), n_values_and_derivatives, &v[i][d][0]);
 
   for (unsigned int i=0; i<n_tensor_pols; ++i)
     {
@@ -309,7 +353,7 @@ compute (const Point<dim>            &p,
         {
           values[i] = 1;
           for (unsigned int x=0; x<dim; ++x)
-            values[i] *= v(x,indices[x])[0];
+            values[i] *= v[indices[x]][x][0];
         }
 
       if (update_grads)
@@ -317,7 +361,7 @@ compute (const Point<dim>            &p,
           {
             grads[i][d] = 1.;
             for (unsigned int x=0; x<dim; ++x)
-              grads[i][d] *= v(x,indices[x])[d==x];
+              grads[i][d] *= v[indices[x]][x][d==x];
           }
 
       if (update_grad_grads)
@@ -332,7 +376,7 @@ compute (const Point<dim>            &p,
                   if (d2==x) ++derivative;
 
                   grad_grads[i][d1][d2]
-                  *= v(x,indices[x])[derivative];
+                  *= v[indices[x]][x][derivative];
                 }
             }
 
@@ -350,7 +394,7 @@ compute (const Point<dim>            &p,
                     if (d3==x) ++derivative;
 
                     third_derivatives[i][d1][d2][d3]
-                    *= v(x,indices[x])[derivative];
+                    *= v[indices[x]][x][derivative];
                   }
               }
 
@@ -370,7 +414,7 @@ compute (const Point<dim>            &p,
                       if (d4==x) ++derivative;
 
                       fourth_derivatives[i][d1][d2][d3][d4]
-                      *= v(x,indices[x])[derivative];
+                      *= v[indices[x]][x][derivative];
                     }
                 }
     }