]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Added IntegratedLegendreSZ which implements the integrated Legendre polynomials descr...
authorRoss Kynch <rkynch@gmail.com>
Tue, 30 May 2017 21:54:53 +0000 (22:54 +0100)
committerRoss Kynch <rkynch@gmail.com>
Wed, 7 Jun 2017 16:46:24 +0000 (17:46 +0100)
doc/news/changes/minor/20170606Kynch [new file with mode: 0644]
include/deal.II/base/polynomials_integrated_legendre_sz.h [new file with mode: 0644]
source/base/CMakeLists.txt
source/base/polynomials_integrated_legendre_sz.cc [new file with mode: 0644]
tests/base/polynomials_integrated_legendre_sz.cc [new file with mode: 0644]
tests/base/polynomials_integrated_legendre_sz.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20170606Kynch b/doc/news/changes/minor/20170606Kynch
new file mode 100644 (file)
index 0000000..a21d70e
--- /dev/null
@@ -0,0 +1,5 @@
+New: Added a new polynomial class IntegratedLegendreSZ. This implements
+the integrated Legendre polynomials described in the 2006 PhD thesis
+of Sabine Zaglmayr.
+<br>
+(Ross Kynch, 2017/06/06)
diff --git a/include/deal.II/base/polynomials_integrated_legendre_sz.h b/include/deal.II/base/polynomials_integrated_legendre_sz.h
new file mode 100644 (file)
index 0000000..e55e789
--- /dev/null
@@ -0,0 +1,92 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 - 2017 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii__polynomials_integrated_legendre_sz_h
+#define dealii__polynomials_integrated_legendre_sz_h
+
+#include <deal.II/base/config.h>
+#include <deal.II/base/polynomial.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/thread_management.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+/**
+ * Class implementing the integrated Legendre polynomials described in the PhD thesis of Sabine Zaglmayer.
+ *
+ * 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:
+ *
+ * - $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:
+ *
+ * $L_{n} = \frac{l_{n} - l_{n-2}}{2n-1)}$
+ *
+ */
+class IntegratedLegendreSZ : public Polynomials::Polynomial<double>
+{
+public:
+  /**
+   * Constructor generating the coefficient of the polynomials up to degree p.
+   */
+  IntegratedLegendreSZ (const unsigned int p);
+
+
+  /**
+   * 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.
+   */
+  static const std::vector<double> &get_coefficients (const unsigned int k);
+};
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index 4b8d56ade97c9404e74e94291373b6a335a86505..4625cb74183eead0cac8544f0d0b6f1129010475 100644 (file)
@@ -51,6 +51,7 @@ SET(_src
   polynomials_bernstein.cc
   polynomials_bdm.cc
   polynomials_nedelec.cc
+  polynomials_integrated_legendre_sz.cc
   polynomial_space.cc
   polynomials_p.cc
   polynomials_piecewise.cc
diff --git a/source/base/polynomials_integrated_legendre_sz.cc b/source/base/polynomials_integrated_legendre_sz.cc
new file mode 100644 (file)
index 0000000..21df96a
--- /dev/null
@@ -0,0 +1,146 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 - 2017 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+#include <deal.II/base/polynomials_integrated_legendre_sz.h>
+
+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))
+{}
+
+
+
+void IntegratedLegendreSZ::compute_coefficients (const unsigned int k_)
+{
+  unsigned int k = k_;
+
+  // 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);
+
+  // The first 2 coefficients are hard-coded
+  if (k==0)   k=1;
+
+
+  // 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
+    {
+      // 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);
+        }
+    }
+}
+
+
+
+const std::vector<double> &IntegratedLegendreSZ::get_coefficients (const unsigned int k)
+{
+  // first make sure the coefficients get computed if so necessary
+  compute_coefficients (k);
+
+  // 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];
+}
+
+
+
+std::vector<Polynomials::Polynomial<double> >
+IntegratedLegendreSZ::generate_complete_basis (const unsigned int degree)
+{
+  std::vector<Polynomials::Polynomial<double> > v;
+  v.reserve(degree + 1);
+  for (unsigned int i=0; i<=degree; ++i)
+    {
+      v.push_back (IntegratedLegendreSZ(i));
+    }
+  return v;
+}
+
+
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/tests/base/polynomials_integrated_legendre_sz.cc b/tests/base/polynomials_integrated_legendre_sz.cc
new file mode 100644 (file)
index 0000000..b027bfe
--- /dev/null
@@ -0,0 +1,70 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// This tests the stability of the polynomial evaluation of IntegratedLegendreSZ.
+
+#include "../tests.h"
+#include <iomanip>
+#include <fstream>
+#include <cmath>
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/polynomial.h>
+#include <deal.II/base/polynomials_integrated_legendre_sz.h>
+#include <deal.II/base/quadrature_lib.h>
+
+
+using namespace Polynomials;
+
+
+void check_at_one (const std::vector<Polynomial<double> > &p)
+{
+  // Ignore first two polynomials as the integrated Legendre polynomials are only defined
+  // for degree > 1, it is only added to maintain the recursive relation.
+
+  deallog << "Function value of polynomial at right end point: ";
+  for (unsigned int i=2; i<p.size(); ++i)
+    {
+      deallog << '.';
+      const double y = p[i].value(1.);
+      if (std::fabs(y) > 1e-13)
+        deallog << "Error1  lg y=" << std::log10(std::fabs(y)) << std::endl;
+    }
+  deallog << std::endl;
+}
+
+
+
+void
+check_poly (const unsigned int n)
+{
+  deallog << "Degree: " << n+1 << std::endl;
+  std::vector<Polynomial<double> > p = IntegratedLegendreSZ::generate_complete_basis(n);
+  check_at_one (p);
+  deallog << std::endl;
+}
+
+
+
+int main()
+{
+  std::ofstream logfile("output");
+  deallog << std::setprecision(3);
+  deallog.attach(logfile);
+  deallog.threshold_double(1.e-10);
+
+  check_poly(25);
+}
diff --git a/tests/base/polynomials_integrated_legendre_sz.output b/tests/base/polynomials_integrated_legendre_sz.output
new file mode 100644 (file)
index 0000000..8bc294c
--- /dev/null
@@ -0,0 +1,4 @@
+
+DEAL::Degree: 26
+DEAL::Function value of polynomial at right end point: ........................
+DEAL::

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.