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)}$
*
{
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
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_)
+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;
}
projects / dealii.git / commitdiff