From d566de88673fb5151c5f78d8bf5f0b2073d70787 Mon Sep 17 00:00:00 2001
From: Rene Gassmoeller <rene.gassmoeller@mailbox.org>
Date: Wed, 12 Oct 2016 14:09:20 -0600
Subject: [PATCH] Provide paths for static and dynamic memory.

---
 include/deal.II/base/polynomial.h            |  8 ++---
 include/deal.II/base/polynomials_piecewise.h |  8 ++---
 source/base/polynomial.cc                    | 30 +++++++++---------
 source/base/polynomials_piecewise.cc         | 16 +++++-----
 source/base/tensor_product_polynomials.cc    | 32 ++++++++++++++------
 5 files changed, 53 insertions(+), 41 deletions(-)

diff --git a/include/deal.II/base/polynomial.h b/include/deal.II/base/polynomial.h
index 22c60b963b..8e4f22a162 100644
--- a/include/deal.II/base/polynomial.h
+++ b/include/deal.II/base/polynomial.h
@@ -105,16 +105,16 @@ 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>x</tt>.  <tt>values[i], i=0,...,n_derivatives</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.
+     * determined by @p n_derivatives and @p values has to provide sufficient
+     * space for @p n_derivatives + 1 values.
      *
      * This function uses the Horner scheme for numerical stability of the
      * evaluation.
      */
     void value (const number         x,
-                const unsigned int values_size,
+                const unsigned int n_derivatives,
                 number *values) const;
 
     /**
diff --git a/include/deal.II/base/polynomials_piecewise.h b/include/deal.II/base/polynomials_piecewise.h
index e070407690..a79933c0ed 100644
--- a/include/deal.II/base/polynomials_piecewise.h
+++ b/include/deal.II/base/polynomials_piecewise.h
@@ -98,10 +98,10 @@ 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>x</tt>.  <tt>values[i], i=0,...,n_derivatives</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.
+     * determined by @p n_derivatives and @p values has to provide sufficient
+     * space for @p n_derivatives + 1 values.
      *
      * Note that all the derivatives evaluate to zero at the border between
      * intervals (assuming exact arithmetics) in the interior of the unit
@@ -112,7 +112,7 @@ namespace Polynomials
      * make sense.
      */
     void value (const number         x,
-                const unsigned int values_size,
+                const unsigned int n_derivatives,
                 number *values) const;
 
     /**
diff --git a/source/base/polynomial.cc b/source/base/polynomial.cc
index 32b03950eb..cedc1f4936 100644
--- a/source/base/polynomial.cc
+++ b/source/base/polynomial.cc
@@ -103,7 +103,7 @@ namespace Polynomials
   {
     Assert (values.size() > 0, ExcZero());
 
-    value(x,values.size(),&values[0]);
+    value(x,values.size()-1,&values[0]);
   }
 
 
@@ -111,10 +111,10 @@ namespace Polynomials
   template <typename number>
   void
   Polynomial<number>::value (const number         x,
-                             const unsigned int values_size,
+                             const unsigned int n_derivatives,
                              number *values) const
   {
-    Assert(values_size > 0, ExcZero());
+    Assert(n_derivatives >= 0, ExcZero());
 
     // evaluate Lagrange polynomial and derivatives
     if (in_lagrange_product_form == true)
@@ -123,11 +123,11 @@ namespace Polynomials
         // form (x-x_1)*(x-x_2)*...*(x-x_n), expand the derivatives like
         // automatic differentiation does.
         const unsigned int n_supp = lagrange_support_points.size();
-        switch (values_size)
+        switch (n_derivatives)
           {
           default:
             values[0] = 1;
-            for (unsigned int d=1; d<values_size; ++d)
+            for (unsigned int d=1; d<=n_derivatives; ++d)
               values[d] = 0;
             for (unsigned int i=0; i<n_supp; ++i)
               {
@@ -139,7 +139,7 @@ namespace Polynomials
                 // i.e., expand value v and derivative one). since we reuse a
                 // value from the next lower derivative from the steps before,
                 // need to start from the highest derivative
-                for (unsigned int k=values_size-1; k>0; --k)
+                for (unsigned int k=n_derivatives; k>0; --k)
                   values[k] = (values[k] * v + values[k-1]);
                 values[0] *= v;
               }
@@ -153,7 +153,7 @@ namespace Polynomials
             // p^(n)(x)/k! into the actual form of the derivative
             {
               number k_faculty = 1;
-              for (unsigned int k=0; k<values_size; ++k)
+              for (unsigned int k=0; k<=n_derivatives; ++k)
                 {
                   values[k] *= k_faculty * lagrange_weight;
                   k_faculty *= static_cast<number>(k+1);
@@ -161,10 +161,10 @@ namespace Polynomials
             }
             break;
 
-          // manually implement size 1 (values only), size 2 (value + first
-          // derivative), and size 3 (up to second derivative) since they
+          // manually implement case 0 (values only), case 1 (value + first
+          // derivative), and case 2 (up to second derivative) since they
           // might be called often. then, we can unroll the loop.
-          case 1:
+          case 0:
             values[0] = 1;
             for (unsigned int i=0; i<n_supp; ++i)
               {
@@ -174,7 +174,7 @@ namespace Polynomials
             values[0] *= lagrange_weight;
             break;
 
-          case 2:
+          case 1:
             values[0] = 1;
             values[1] = 0;
             for (unsigned int i=0; i<n_supp; ++i)
@@ -187,7 +187,7 @@ namespace Polynomials
             values[1] *= lagrange_weight;
             break;
 
-          case 3:
+          case 2:
             values[0] = 1;
             values[1] = 0;
             values[2] = 0;
@@ -211,7 +211,7 @@ namespace Polynomials
     // if we only need the value, then call the other function since that is
     // significantly faster (there is no need to allocate and free memory,
     // which is really expensive compared to all the other operations!)
-    if (values_size == 1)
+    if (n_derivatives == 0)
       {
         values[0] = value(x);
         return;
@@ -226,7 +226,7 @@ namespace Polynomials
     // loop over all requested derivatives. note that derivatives @p{j>m} are
     // necessarily zero, as they differentiate the polynomial more often than
     // the highest power is
-    const unsigned int min_valuessize_m=std::min(values_size, m);
+    const unsigned int min_valuessize_m=std::min(n_derivatives+1, m);
     for (unsigned int j=0; j<min_valuessize_m; ++j)
       {
         for (int k=m-2; k>=static_cast<int>(j); --k)
@@ -237,7 +237,7 @@ namespace Polynomials
       }
 
     // fill higher derivatives by zero
-    for (unsigned int j=min_valuessize_m; j<values_size; ++j)
+    for (unsigned int j=min_valuessize_m; j<=n_derivatives; ++j)
       values[j] = 0;
   }
 
diff --git a/source/base/polynomials_piecewise.cc b/source/base/polynomials_piecewise.cc
index 614b4b3d85..cd082ea688 100644
--- a/source/base/polynomials_piecewise.cc
+++ b/source/base/polynomials_piecewise.cc
@@ -47,7 +47,7 @@ namespace Polynomials
   {
     Assert (values.size() > 0, ExcZero());
 
-    value(x,values.size(),&values[0]);
+    value(x,values.size()-1,&values[0]);
   }
 
 
@@ -55,10 +55,10 @@ namespace Polynomials
   template <typename number>
   void
   PiecewisePolynomial<number>::value (const number         x,
-                                      const unsigned int values_size,
+                                      const unsigned int n_derivatives,
                                       number *values) const
   {
-    Assert (values_size > 0, ExcZero());
+    Assert (n_derivatives >= 0, ExcZero());
 
     // shift polynomial if necessary
     number y = x;
@@ -72,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<=n_derivatives; ++k)
                   values[k] = 0;
                 return;
               }
@@ -89,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<=n_derivatives; ++k)
                   values[k] = 0;
                 return;
               }
@@ -106,16 +106,16 @@ namespace Polynomials
              (interval < n_intervals-1 || derivative_change_sign == -1.)))
           {
             values[0] = value(x);
-            for (unsigned int d=1; d<values_size; ++d)
+            for (unsigned int d=1; d<=n_derivatives; ++d)
               values[d] = 0;
             return;
           }
       }
 
-    polynomial.value(y, values_size, values);
+    polynomial.value(y, n_derivatives, values);
 
     // change sign if necessary
-    for (unsigned int j=1; j<values_size; j+=2)
+    for (unsigned int j=1; j<=n_derivatives; j+=2)
       values[j] *= derivative_change_sign;
   }
 
diff --git a/source/base/tensor_product_polynomials.cc b/source/base/tensor_product_polynomials.cc
index fa382021db..5729912339 100644
--- a/source/base/tensor_product_polynomials.cc
+++ b/source/base/tensor_product_polynomials.cc
@@ -17,7 +17,6 @@
 #include <deal.II/base/polynomials_piecewise.h>
 #include <deal.II/base/exceptions.h>
 #include <deal.II/base/table.h>
-
 #include <deal.II/base/std_cxx11/array.h>
 
 DEAL_II_NAMESPACE_OPEN
@@ -313,6 +312,9 @@ compute (const Point<dim>            &p,
                     case 2:
                       index = polynomial*factor + f;
                       break;
+                    default:
+                      ExcNotImplemented();
+                      break;
                     }
 
                   const unsigned int i = index_map_inverse[index];
@@ -326,17 +328,27 @@ compute (const Point<dim>            &p,
 
   // Compute the values (and derivatives, if
   // necessary) of all polynomials at this
-  // 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
+  // evaluation point. To avoid expensive memory allocation,
+  // we use a small amount of memory on the stack, and store the
+  // result in an array (that has enough
   // fields for any evaluation of values and
-  // 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());
+  // derivatives, up to the 4th derivative, for up to 20 polynomials).
+  // If someone uses a larger number of
+  // polynomials, we need to allocate more memory on the heap.
+  std_cxx11::array<std_cxx11::array<double,5>, dim> *v;
+  std_cxx11::array<std_cxx11::array<std_cxx11::array<double,5>, dim>, 20> small_array;
+  std::vector<std_cxx11::array<std_cxx11::array<double,5>, dim> > large_array;
+
+  const unsigned int n_polynomials = polynomials.size();
+  if (n_polynomials > 20)
+    {
+      large_array.resize(n_polynomials);
+      v = &large_array[0];
+    }
+  else
+    v = &small_array[0];
 
-  for (unsigned int i=0; i<polynomials.size(); ++i)
+  for (unsigned int i=0; i<n_polynomials; ++i)
     for (unsigned int d=0; d<dim; ++d)
       polynomials[i].value(p(d), n_values_and_derivatives, &v[i][d][0]);
 
-- 
2.39.5