]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fixed comment spacing in header. Fixed some typos in documentation. Replaced use... 4496/head
authorRoss Kynch <rkynch@gmail.com>
Thu, 8 Jun 2017 22:41:00 +0000 (23:41 +0100)
committerRoss Kynch <rkynch@gmail.com>
Fri, 9 Jun 2017 00:14:51 +0000 (01:14 +0100)
include/deal.II/base/polynomials_integrated_legendre_sz.h
source/base/polynomials_integrated_legendre_sz.cc

index e55e789486a1572b580dd9c6ee3c01101d28327e..55c8eb9181111639a0da04f8ce7225fcfeb0d6a6 100644 (file)
 DEAL_II_NAMESPACE_OPEN
 
 /**
- * Class implementing the integrated Legendre polynomials described in the PhD thesis of Sabine Zaglmayer.
+ * Class implementing the integrated Legendre polynomials described in the PhD
+ * thesis of Sabine Zaglmayr.
  *
- * This class was written based upon the existing deal.II Legendre class as a base, but with the coefficents adjusted
- * so that the recursive formula is for the integrated Legendre polynomials described in the PhD thesis of
- * Sabine Zaglmayer. The polynomials can be generated recursively from:
+ * This class was written based upon the existing deal.II Legendre class as a
+ * base, but with the coefficents adjusted so that the recursive formula is for
+ * the integrated Legendre polynomials described in the PhD thesis of Sabine
+ * Zaglmayr. The polynomials can be generated recursively from:
  *
  * - $L_{0}(x) = -1$ (added so that it can be generated recursively from 0)
  * - $L_{1}(x) = x$
  * - $L_{2}(x) = \frac{(x^2 - 1)}{2}$
  * - $(n+1)L_{n+1} = (2n-1)L_{n} - (n-2)L_{n-1}$.
  *
- * However, it is also possible to generate them directly from the Legendre polynomials:
+ * However, it is also possible to generate them directly from the Legendre
+ * polynomials:
  *
  * $L_{n} = \frac{l_{n} - l_{n-2}}{2n-1)}$
  *
@@ -46,45 +49,21 @@ class IntegratedLegendreSZ : public Polynomials::Polynomial<double>
 {
 public:
   /**
-   * Constructor generating the coefficient of the polynomials up to degree p.
+   * Constructor generating the coefficients of the polynomials at degree p.
    */
   IntegratedLegendreSZ (const unsigned int p);
 
-
   /**
-   * Returns the complete set of Integrated Legendre polynomials up to the given degree.
+   * Returns the complete set of Integrated Legendre polynomials up to the
+   * given degree.
    */
   static std::vector<Polynomials::Polynomial<double>> generate_complete_basis (const unsigned int degree);
 
-
 private:
   /**
-   * Lock that guarantees that at most one thread is changing and accessing the recursive_coefficients array.
-   */
-  static Threads::Mutex coefficients_lock;
-
-
-  /**
-   * Vector with already computed coefficients. For each degree of the
-   * polynomial, we keep one pointer to the list of coefficients; we do so
-   * rather than keeping a vector of vectors in order to simplify
-   * programming multithread-safe. In order to avoid memory leak, we use a
-   * shared_ptr in order to correctly free the memory of the vectors when
-   * the global destructor is called.
-   */
-  static std::vector<std::shared_ptr<const std::vector<double>>> recursive_coefficients;
-
-
-  /**
-   * Main function to compute the co-efficients of the polyonial at degree p.
-   */
-  static void compute_coefficients (const unsigned int p);
-
-
-  /**
-   * Get coefficients for constructor.
+   * Main function to compute the co-efficients of the polynomial at degree p.
    */
-  static const std::vector<double> &get_coefficients (const unsigned int k);
+  static const std::vector<double> get_coefficients (const unsigned int k);
 };
 
 DEAL_II_NAMESPACE_CLOSE
index 21df96a94c27822e9f60e02438b28cf554e5cff5..2873ee64e92cc7e79f129819b2a6afe6125a54e1 100644 (file)
 
 DEAL_II_NAMESPACE_OPEN
 
-// Reserve space for polynomials up to degree 19.
-std::vector<std::shared_ptr<const std::vector<double>>> IntegratedLegendreSZ::recursive_coefficients(20);
-
-// Define the static mutex member.
-Threads::Mutex IntegratedLegendreSZ::coefficients_lock;
-
-
 IntegratedLegendreSZ::IntegratedLegendreSZ (const unsigned int k)
   :
   Polynomials::Polynomial<double> (get_coefficients(k))
@@ -32,99 +25,48 @@ IntegratedLegendreSZ::IntegratedLegendreSZ (const unsigned int k)
 
 
 
-void IntegratedLegendreSZ::compute_coefficients (const unsigned int k_)
+const std::vector<double> IntegratedLegendreSZ::get_coefficients (const unsigned int k)
 {
-  unsigned int k = k_;
+  std::vector<double> coefficients(k+1);
 
-  // first make sure that no other thread intercepts the operation of this function;
-  // for this, acquire the lock until we quit this function
-  Threads::Mutex::ScopedLock lock(coefficients_lock);
+  // first two polynomials are hard-coded:
+  if (k==0)
+    {
+      coefficients[0] = -1.;
+      return coefficients;
+    }
+  else if (k==1)
+    {
+      coefficients[0] = 0.;
+      coefficients[1] = 1.;
+      return coefficients;
+    }
 
-  // The first 2 coefficients are hard-coded
-  if (k==0)   k=1;
+  // General formula is:
+  // k*L_{k}(x) = (2*k-3)*x*L_{k-1} - (k-3)*L_{k-2}.
+  std::vector<double> coefficients_km2 = get_coefficients(k-2);
+  std::vector<double> coefficients_km1 = get_coefficients(k-1);
 
+  const double a = 1.0 / k;
+  const double b = 2.0*k - 3.0;
+  const double c = k - 3.0;
 
-  // check: does the information already exist?
-  if ((recursive_coefficients.size() < k+1) ||
-      ((recursive_coefficients.size() >= k+1) &&
-       (recursive_coefficients[k] == std::shared_ptr<const std::vector<double> >())))
-    // no, then generate the respective coefficients
+  // To maintain stability, delay the division (multiplication by a) until the end.
+  for (unsigned int i=1; i<=k-2; i++)
     {
-      // make sure that there is enough space in the array for the coefficients,
-      // so we have to resize it to size k+1
-
-      // but it's more complicated than that: we call this function recursively, so if we simply
-      // resize it to k+1 here, then compute the coefficients for degree k-1 by calling this
-      // function recursively, then it will reset the size to k -- not enough for what we want to do below. the
-      // solution therefore is to only resize the size if we are going to *increase* it
-      if (recursive_coefficients.size() < k+1)
-        {
-          recursive_coefficients.resize (k+1);
-        }
-      if (k<=1)
-        {
-          // create coefficients vectors for k=0 and k=1
-          //
-          // allocate the respective later assign it to the coefficients array to make it const
-          std::vector<double> *c0 = new std::vector<double>(1);
-          (*c0)[0] = -1.;
-
-          std::vector<double> *c1 = new std::vector<double>(2);
-          (*c1)[0] = 0.;
-          (*c1)[1] = 1.;
-
-          // now make these arrays const. use shared_ptr for recursive_coefficients because
-          // that avoids a memory leak that would appear if we used plain pointers.
-          recursive_coefficients[0] = std::shared_ptr<const std::vector<double> >(c0);
-          recursive_coefficients[1] = std::shared_ptr<const std::vector<double> >(c1);
-
-        }
-      else
-        {
-          // for larger numbers, compute the coefficients recursively. to do so, we have to release the
-          // lock temporarily to allow the called function to acquire it itself
-          coefficients_lock.release ();
-          compute_coefficients(k-1);
-          coefficients_lock.acquire ();
-
-          std::vector<double> *ck = new std::vector<double>(k+1);
-
-          const double a = 1.0 / k;
-          const double b = 2.0*k - 3.0;
-          const double c = k - 3.0;
-
-          // To maintain stability, delay the division (multiplication by a) until the end.
-
-          (*ck)[k]   = b*(*recursive_coefficients[k-1])[k-1];
-          (*ck)[k-1] = b*(*recursive_coefficients[k-1])[k-2];
-          for (unsigned int i=1; i<= k-2 ; ++i)
-            {
-              (*ck)[i] = b*(*recursive_coefficients[k-1])[i-1] - c*(*recursive_coefficients[k-2])[i];
-            }
-
-          (*ck)[0] = -c*(*recursive_coefficients[k-2])[0];
-
-          for (unsigned int i=0; i<ck->size(); i++)
-            {
-              (*ck)[i] *=a;
-            }
-
-          // finally assign the newly created vector to the const pointer in the/ coefficients array
-          recursive_coefficients[k] = std::shared_ptr<const std::vector<double> >(ck);
-        }
+      coefficients[i] = b*coefficients_km1[i-1] - c*coefficients_km2[i];
     }
-}
-
 
+  coefficients[0] = -c*coefficients_km2[0];
+  coefficients[k]   = b*coefficients_km1[k-1];
+  coefficients[k-1] = b*coefficients_km1[k-2];
 
-const std::vector<double> &IntegratedLegendreSZ::get_coefficients (const unsigned int k)
-{
-  // first make sure the coefficients get computed if so necessary
-  compute_coefficients (k);
+  for (unsigned int i=0; i<coefficients.size(); i++)
+    {
+      coefficients[i] *= a;
+    }
 
-  // then get a pointer to the array of coefficients. do that in a MT safe way
-  Threads::Mutex::ScopedLock lock (coefficients_lock);
-  return *recursive_coefficients[k];
+  return coefficients;
 }
 
 

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.