From: Ross Kynch Date: Thu, 8 Jun 2017 22:41:00 +0000 (+0100) Subject: Fixed comment spacing in header. Fixed some typos in documentation. Replaced use... X-Git-Tag: v9.0.0-rc1~1517^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9cc781dd4168ad9a0471db574070fe42e0678e2d;p=dealii.git Fixed comment spacing in header. Fixed some typos in documentation. Replaced use of pre-computed coefficients and replaced with a simple on-the-fly calculation. --- diff --git a/include/deal.II/base/polynomials_integrated_legendre_sz.h b/include/deal.II/base/polynomials_integrated_legendre_sz.h index e55e789486..55c8eb9181 100644 --- a/include/deal.II/base/polynomials_integrated_legendre_sz.h +++ b/include/deal.II/base/polynomials_integrated_legendre_sz.h @@ -26,18 +26,21 @@ 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 { 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> 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. + * Main function to compute the co-efficients of the polynomial at degree p. */ - static const std::vector &get_coefficients (const unsigned int k); + static const std::vector get_coefficients (const unsigned int k); }; DEAL_II_NAMESPACE_CLOSE diff --git a/source/base/polynomials_integrated_legendre_sz.cc b/source/base/polynomials_integrated_legendre_sz.cc index 21df96a94c..2873ee64e9 100644 --- a/source/base/polynomials_integrated_legendre_sz.cc +++ b/source/base/polynomials_integrated_legendre_sz.cc @@ -18,13 +18,6 @@ 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)) @@ -32,99 +25,48 @@ IntegratedLegendreSZ::IntegratedLegendreSZ (const unsigned int k) -void IntegratedLegendreSZ::compute_coefficients (const unsigned int k_) +const std::vector IntegratedLegendreSZ::get_coefficients (const unsigned int k) { - unsigned int k = k_; + std::vector 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 coefficients_km2 = get_coefficients(k-2); + std::vector 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 >()))) - // 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 *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); - } + 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 &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