From fc605df00bf62afb5a741d67e002f2173da98d23 Mon Sep 17 00:00:00 2001 From: Ross Kynch Date: Tue, 30 May 2017 22:54:53 +0100 Subject: [PATCH] Added IntegratedLegendreSZ which implements the integrated Legendre polynomials described in the PhD thesis of Sabine Zaglmayr. --- doc/news/changes/minor/20170606Kynch | 5 + .../base/polynomials_integrated_legendre_sz.h | 92 +++++++++++ source/base/CMakeLists.txt | 1 + .../polynomials_integrated_legendre_sz.cc | 146 ++++++++++++++++++ .../polynomials_integrated_legendre_sz.cc | 70 +++++++++ .../polynomials_integrated_legendre_sz.output | 4 + 6 files changed, 318 insertions(+) create mode 100644 doc/news/changes/minor/20170606Kynch create mode 100644 include/deal.II/base/polynomials_integrated_legendre_sz.h create mode 100644 source/base/polynomials_integrated_legendre_sz.cc create mode 100644 tests/base/polynomials_integrated_legendre_sz.cc create mode 100644 tests/base/polynomials_integrated_legendre_sz.output diff --git a/doc/news/changes/minor/20170606Kynch b/doc/news/changes/minor/20170606Kynch new file mode 100644 index 0000000000..a21d70e467 --- /dev/null +++ b/doc/news/changes/minor/20170606Kynch @@ -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. +
+(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 index 0000000000..e55e789486 --- /dev/null +++ b/include/deal.II/base/polynomials_integrated_legendre_sz.h @@ -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 +#include +#include +#include +#include + + +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 +{ +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> 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>> 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 &get_coefficients (const unsigned int k); +}; + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/base/CMakeLists.txt b/source/base/CMakeLists.txt index 4b8d56ade9..4625cb7418 100644 --- a/source/base/CMakeLists.txt +++ b/source/base/CMakeLists.txt @@ -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 index 0000000000..21df96a94c --- /dev/null +++ b/source/base/polynomials_integrated_legendre_sz.cc @@ -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_NAMESPACE_OPEN + +// Reserve space for polynomials up to degree 19. +std::vector>> IntegratedLegendreSZ::recursive_coefficients(20); + +// Define the static mutex member. +Threads::Mutex IntegratedLegendreSZ::coefficients_lock; + + +IntegratedLegendreSZ::IntegratedLegendreSZ (const unsigned int k) + : + Polynomials::Polynomial (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 >()))) + // 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 *c0 = new std::vector(1); + (*c0)[0] = -1.; + + std::vector *c1 = new std::vector(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 >(c0); + recursive_coefficients[1] = std::shared_ptr >(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 *ck = new std::vector(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; isize(); i++) + { + (*ck)[i] *=a; + } + + // finally assign the newly created vector to the const pointer in the/ coefficients array + recursive_coefficients[k] = std::shared_ptr >(ck); + } + } +} + + + +const std::vector &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 > +IntegratedLegendreSZ::generate_complete_basis (const unsigned int degree) +{ + std::vector > 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 index 0000000000..b027bfea72 --- /dev/null +++ b/tests/base/polynomials_integrated_legendre_sz.cc @@ -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 +#include +#include + +#include +#include +#include +#include + + +using namespace Polynomials; + + +void check_at_one (const std::vector > &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 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 > 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 index 0000000000..8bc294c7bd --- /dev/null +++ b/tests/base/polynomials_integrated_legendre_sz.output @@ -0,0 +1,4 @@ + +DEAL::Degree: 26 +DEAL::Function value of polynomial at right end point: ........................ +DEAL:: -- 2.39.5