From 85df947d3c2b08fc562b0c7d8720cb67097e2249 Mon Sep 17 00:00:00 2001 From: wolf Date: Tue, 15 Oct 2002 16:44:06 +0000 Subject: [PATCH] Move all 1d polynomials into a namespace Polynomials. git-svn-id: https://svn.dealii.org/trunk@6654 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/polynomial.h | 516 ++++---- deal.II/base/include/base/polynomial_space.h | 2 +- .../include/base/tensor_product_polynomials.h | 2 +- deal.II/base/source/polynomial.cc | 1158 +++++++++-------- deal.II/deal.II/include/fe/mapping_q.h | 1 - deal.II/deal.II/source/fe/fe_dgp.cc | 2 +- .../deal.II/source/fe/fe_dgp_nonparametric.cc | 2 +- deal.II/deal.II/source/fe/fe_dgq.cc | 4 +- deal.II/deal.II/source/fe/fe_q.cc | 10 +- deal.II/deal.II/source/fe/mapping_q.cc | 5 +- deal.II/doc/news/2002/c-3-4.html | 7 + tests/base/polynomial1d.cc | 3 + tests/base/polynomial_test.cc | 3 + 13 files changed, 874 insertions(+), 841 deletions(-) diff --git a/deal.II/base/include/base/polynomial.h b/deal.II/base/include/base/polynomial.h index e5eadbab0b..5eb0359512 100644 --- a/deal.II/base/include/base/polynomial.h +++ b/deal.II/base/include/base/polynomial.h @@ -22,6 +22,13 @@ #include +/** + * A namespace in which classes relating to the description of + * 1d polynomial spaces are declared. + */ +namespace Polynomials +{ + /** * Base class for all 1D polynomials. A polynomial is represented in * this class by its coefficients, which are set through the @@ -31,164 +38,164 @@ * * @author Ralf Hartmann, Guido Kanschat, 2000 */ -template -class Polynomial : public Subscriptor -{ - public: - /** - * Constructor. The coefficients - * of the polynomial are passed - * as arguments, and denote the - * polynomial @p{\sum_i a[i] - * x^i}, i.e. the first element - * of the array denotes the - * constant term, the second the - * linear one, and so on. The - * degree of the polynomial - * represented by this object is - * thus the number of elements in - * the @p{coefficient} array - * minus one. - */ - Polynomial (const std::vector &coefficients); + template + class Polynomial : public Subscriptor + { + public: + /** + * Constructor. The coefficients + * of the polynomial are passed + * as arguments, and denote the + * polynomial @p{\sum_i a[i] + * x^i}, i.e. the first element + * of the array denotes the + * constant term, the second the + * linear one, and so on. The + * degree of the polynomial + * represented by this object is + * thus the number of elements in + * the @p{coefficient} array + * minus one. + */ + Polynomial (const std::vector &coefficients); - /** - * Return the value of this - * polynomial at the given point. - * - * This function uses the Horner - * scheme for numerical stability - * of the evaluation. - */ - number value (const number x) const; + /** + * Return the value of this + * polynomial at the given point. + * + * This function uses the Horner + * scheme for numerical stability + * of the evaluation. + */ + number value (const number x) const; - /** - * Return the values and the - * derivatives of the - * @p{Polynomial} at point @p{x}. - * @p{values[i], - * i=0,...,values.size()-1} - * includes the @p{i}th - * derivative. The number of - * derivatives to be computed is - * thus determined by the size of - * the array passed. - * - * This function uses the Horner - * scheme for numerical stability - * of the evaluation. - */ - void value (const number x, - std::vector &values) const; + /** + * Return the values and the + * derivatives of the + * @p{Polynomial} at point @p{x}. + * @p{values[i], + * i=0,...,values.size()-1} + * includes the @p{i}th + * derivative. The number of + * derivatives to be computed is + * thus determined by the size of + * the array passed. + * + * This function uses the Horner + * scheme for numerical stability + * of the evaluation. + */ + void value (const number x, + std::vector &values) const; - /** - * Degree of the polynomial. This - * is the degree reflected by the - * number of coefficients - * provided by the - * constructor. Leading non-zero - * coefficients are not treated - * separately. - */ - unsigned int degree () const; + /** + * Degree of the polynomial. This + * is the degree reflected by the + * number of coefficients + * provided by the + * constructor. Leading non-zero + * coefficients are not treated + * separately. + */ + unsigned int degree () const; - /** - * Scale the abscissa of the - * polynomial. Given the - * polynomial $p(t)$ and the - * scaling $t = ax$, then the - * result of this operation is - * the polynomial $q$, such that - * $q(x) = p(t)$. - * - * The operation is performed in - * place. - */ - void scale (const number factor); + /** + * Scale the abscissa of the + * polynomial. Given the + * polynomial $p(t)$ and the + * scaling $t = ax$, then the + * result of this operation is + * the polynomial $q$, such that + * $q(x) = p(t)$. + * + * The operation is performed in + * place. + */ + void scale (const number factor); - /** - * Shift the abscissa oft the - * polynomial. Given the - * polynomial $p(t)$ and the - * shift $t = x + a$, then the - * result of this operation is - * the polynomial $q$, such that - * $q(x) = p(t)$. - * - * The template parameter allows - * to compute the new - * coefficients with higher - * accuracy, since all - * computations are performed - * with type @p{number2}. This - * may be necessary, since this - * operation involves a big - * number of additions. On a Sun - * Sparc Ultra with Solaris 2.8, - * the difference between - * @p{double} and @p{long double} - * was not significant, though. - * - * The operation is performed in - * place, i.e. the coefficients - * of the present object are - * changed. - */ - template - void shift (const number2 offset); + /** + * Shift the abscissa oft the + * polynomial. Given the + * polynomial $p(t)$ and the + * shift $t = x + a$, then the + * result of this operation is + * the polynomial $q$, such that + * $q(x) = p(t)$. + * + * The template parameter allows + * to compute the new + * coefficients with higher + * accuracy, since all + * computations are performed + * with type @p{number2}. This + * may be necessary, since this + * operation involves a big + * number of additions. On a Sun + * Sparc Ultra with Solaris 2.8, + * the difference between + * @p{double} and @p{long double} + * was not significant, though. + * + * The operation is performed in + * place, i.e. the coefficients + * of the present object are + * changed. + */ + template + void shift (const number2 offset); - /** - * Print coefficients. - */ - void print(std::ostream& out) const; + /** + * Print coefficients. + */ + void print(std::ostream& out) const; - /** - * Exception - */ - DeclException0 (ExcEmptyArray); + /** + * Exception + */ + DeclException0 (ExcEmptyArray); - /** - * Exception - */ - DeclException0 (ExcVoidPolynomial); + /** + * Exception + */ + DeclException0 (ExcVoidPolynomial); - protected: + protected: - /** - * This function performs the - * actual scaling. - */ - static void scale (std::vector &coefficients, - const number factor); + /** + * This function performs the + * actual scaling. + */ + static void scale (std::vector &coefficients, + const number factor); - /** - * This function performs the - * actual shift - */ - template - static void shift (std::vector &coefficients, - const number2 shift); + /** + * This function performs the + * actual shift + */ + template + static void shift (std::vector &coefficients, + const number2 shift); - /** - * Multiply polynomial by a factor. - */ - static void multiply (std::vector& coefficients, - const number factor); + /** + * Multiply polynomial by a factor. + */ + static void multiply (std::vector& coefficients, + const number factor); - /** - * Coefficients of the polynomial - * $\sum_i a_i x^i$. This vector - * is filled by the constructor - * of this class and may be - * passed down by derived - * classes. - * - * This vector cannot be constant - * since we want to allow copying - * of polynomials. - */ - std::vector coefficients; -}; + /** + * Coefficients of the polynomial + * $\sum_i a_i x^i$. This vector + * is filled by the constructor + * of this class and may be + * passed down by derived + * classes. + * + * This vector cannot be constant + * since we want to allow copying + * of polynomials. + */ + std::vector coefficients; + }; @@ -210,56 +217,56 @@ class Polynomial : public Subscriptor * * @author Ralf Hartmann, 2000 */ -class LagrangeEquidistant: public Polynomial -{ - public: - /** - * Constructor. Takes the order - * @p{n} of the Lagrangian - * polynom and the index - * @p{support_point} of the - * support point. Fills the - * @p{coefficients} of the base - * class @p{Polynomial}. - */ - LagrangeEquidistant (const unsigned int n, - const unsigned int support_point); + class LagrangeEquidistant: public Polynomial + { + public: + /** + * Constructor. Takes the order + * @p{n} of the Lagrangian + * polynom and the index + * @p{support_point} of the + * support point. Fills the + * @p{coefficients} of the base + * class @p{Polynomial}. + */ + LagrangeEquidistant (const unsigned int n, + const unsigned int support_point); - /** - * Return a vector of polynomial - * objects of order @p{degree}, - * which then spans the full - * space of polynomials up to the - * given degree. The polynomials - * are generated by calling the - * destructor of this class with - * the same degree but support - * point running from zero to - * @p{degree}. This function may - * be used to initialize the - * @ref{TensorProductPolynomials} - * and @ref{PolynomialSpace} - * classes. - */ - static - std::vector > - generate_complete_basis (const unsigned int degree); + /** + * Return a vector of polynomial + * objects of order @p{degree}, + * which then spans the full + * space of polynomials up to the + * given degree. The polynomials + * are generated by calling the + * destructor of this class with + * the same degree but support + * point running from zero to + * @p{degree}. This function may + * be used to initialize the + * @ref{TensorProductPolynomials} + * and @ref{PolynomialSpace} + * classes. + */ + static + std::vector > + generate_complete_basis (const unsigned int degree); - private: + private: - /** - * Computes the @p{coefficients} - * of the base class - * @p{Polynomial}. This function - * is @p{static} to allow to be - * called in the - * constructor. - */ - static - std::vector - compute_coefficients (const unsigned int n, - const unsigned int support_point); -}; + /** + * Computes the @p{coefficients} + * of the base class + * @p{Polynomial}. This function + * is @p{static} to allow to be + * called in the + * constructor. + */ + static + std::vector + compute_coefficients (const unsigned int n, + const unsigned int support_point); + }; /** @@ -274,76 +281,81 @@ class LagrangeEquidistant: public Polynomial * * @author Guido Kanschat, 2000 */ -template -class Legendre : public Polynomial -{ - public: - /** - * Constructor for polynomial of - * order @p{k}. - */ - Legendre (const unsigned int k); + template + class Legendre : public Polynomial + { + public: + /** + * Constructor for polynomial of + * order @p{k}. + */ + Legendre (const unsigned int k); - /** - * Return a vector of Legendre - * polynomial objects of orders - * zero through @p{degree}, which - * then spans the full space of - * polynomials up to the given - * degree. This function may be - * used to initialize the - * @ref{TensorProductPolynomials} - * and @ref{PolynomialSpace} - * classes. - */ - static - std::vector > - generate_complete_basis (const unsigned int degree); + /** + * Return a vector of Legendre + * polynomial objects of orders + * zero through @p{degree}, which + * then spans the full space of + * polynomials up to the given + * degree. This function may be + * used to initialize the + * @ref{TensorProductPolynomials} + * and @ref{PolynomialSpace} + * classes. + */ + static + std::vector > + generate_complete_basis (const unsigned int degree); - private: - /** - * Coefficients for the interval $[0,1]$. - */ - static std::vector *> shifted_coefficients; + private: + /** + * Coefficients for the interval $[0,1]$. + */ + static std::vector *> shifted_coefficients; - /** - * 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. - */ - static std::vector *> recursive_coefficients; + /** + * 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. + */ + static std::vector *> recursive_coefficients; - /** - * Compute coefficients recursively. - */ - static void compute_coefficients (const unsigned int k); + /** + * Compute coefficients recursively. + */ + static void compute_coefficients (const unsigned int k); - /** - * Get coefficients for - * constructor. This way, it can - * use the non-standard - * constructor of - * @ref{Polynomial}. - */ - static const std::vector & - get_coefficients (const unsigned int k); -}; + /** + * Get coefficients for + * constructor. This way, it can + * use the non-standard + * constructor of + * @ref{Polynomial}. + */ + static const std::vector & + get_coefficients (const unsigned int k); + }; + +} /* -------------------------- inline functions --------------------- */ -template -inline -unsigned int -Polynomial::degree () const +namespace Polynomials { - Assert (coefficients.size()>0, ExcVoidPolynomial()); - return coefficients.size() - 1; + template + inline + unsigned int + Polynomial::degree () const + { + Assert (coefficients.size()>0, ExcVoidPolynomial()); + return coefficients.size() - 1; + } } #endif diff --git a/deal.II/base/include/base/polynomial_space.h b/deal.II/base/include/base/polynomial_space.h index a2b2cc98dd..c26f2cfcfb 100644 --- a/deal.II/base/include/base/polynomial_space.h +++ b/deal.II/base/include/base/polynomial_space.h @@ -143,7 +143,7 @@ class PolynomialSpace * polynomials given to the * constructor. */ - const std::vector > polynomials; + const std::vector > polynomials; /** * Store the precomputed value diff --git a/deal.II/base/include/base/tensor_product_polynomials.h b/deal.II/base/include/base/tensor_product_polynomials.h index 06801ee2aa..9a737445d8 100644 --- a/deal.II/base/include/base/tensor_product_polynomials.h +++ b/deal.II/base/include/base/tensor_product_polynomials.h @@ -182,7 +182,7 @@ class TensorProductPolynomials * polynomials given to the * constructor. */ - const std::vector > polynomials; + const std::vector > polynomials; /** * Number of tensor product diff --git a/deal.II/base/source/polynomial.cc b/deal.II/base/source/polynomial.cc index bb12b4452c..61b866cf3a 100644 --- a/deal.II/base/source/polynomial.cc +++ b/deal.II/base/source/polynomial.cc @@ -33,464 +33,467 @@ namespace }; +namespace Polynomials +{ + // -------------------- class Polynomial ---------------- // -template -Polynomial::Polynomial (const std::vector &a): - coefficients(a) -{} + template + Polynomial::Polynomial (const std::vector &a): + coefficients(a) + {} -template -number -Polynomial::value (const number x) const -{ - Assert (coefficients.size() > 0, ExcVoidPolynomial()); - const unsigned int m=coefficients.size(); + template + number + Polynomial::value (const number x) const + { + Assert (coefficients.size() > 0, ExcVoidPolynomial()); + const unsigned int m=coefficients.size(); - // Horner scheme - number value = coefficients.back(); - for (int k=m-2; k>=0; --k) - value = value*x + coefficients[k]; + // Horner scheme + number value = coefficients.back(); + for (int k=m-2; k>=0; --k) + value = value*x + coefficients[k]; - return value; -} + return value; + } -template -void -Polynomial::value (const number x, - std::vector &values) const -{ - Assert (coefficients.size() > 0, ExcVoidPolynomial()); - Assert (values.size() > 0, ExcEmptyArray()); - const unsigned int values_size=values.size(); + template + void + Polynomial::value (const number x, + std::vector &values) const + { + Assert (coefficients.size() > 0, ExcVoidPolynomial()); + Assert (values.size() > 0, ExcEmptyArray()); + const unsigned int values_size=values.size(); - // if we only need the value, then - // call the other function since - // that is significantly faster - // (there is no need to allocate - // and free memory, which is really - // expensive compared to all the - // other operations!) - if (values_size == 1) - { - values[0] = value(x); - return; - }; - - // if there are derivatives needed, - // then do it properly by the - // full Horner scheme - const unsigned int m=coefficients.size(); - std::vector a(coefficients); - unsigned int j_faculty=1; - - // loop over all requested - // derivatives. note that - // derivatives @p{j>m} are - // necessarily zero, as they - // differentiate the polynomial - // more often than the highest - // power is - const unsigned int min_valuessize_m=std::min(values_size, m); - for (unsigned int j=0; j=static_cast(j); --k) - a[k]+=x*a[k+1]; - values[j]=j_faculty*a[j]; - - j_faculty*=j+1; - } - - // fill higher derivatives by zero - for (unsigned int j=min_valuessize_m; j a(coefficients); + unsigned int j_faculty=1; + + // loop over all requested + // derivatives. note that + // derivatives @p{j>m} are + // necessarily zero, as they + // differentiate the polynomial + // more often than the highest + // power is + const unsigned int min_valuessize_m=std::min(values_size, m); + for (unsigned int j=0; j=static_cast(j); --k) + a[k]+=x*a[k+1]; + values[j]=j_faculty*a[j]; + + j_faculty*=j+1; + } + // fill higher derivatives by zero + for (unsigned int j=min_valuessize_m; j -void -Polynomial::scale(std::vector& coefficients, - const number factor) -{ - double f = 1.; - for (typename std::vector::iterator c = coefficients.begin(); - c != coefficients.end(); ++c) - { - *c *= f; - f *= factor; - } -} + template + void + Polynomial::scale(std::vector& coefficients, + const number factor) + { + double f = 1.; + for (typename std::vector::iterator c = coefficients.begin(); + c != coefficients.end(); ++c) + { + *c *= f; + f *= factor; + } + } -template -void -Polynomial::scale(const number factor) -{ - scale (coefficients, factor); -} + template + void + Polynomial::scale(const number factor) + { + scale (coefficients, factor); + } -template -void -Polynomial::multiply(std::vector& coefficients, - const number factor) -{ - for (typename std::vector::iterator c = coefficients.begin(); - c != coefficients.end(); ++c) - *c *= factor; -} + + template + void + Polynomial::multiply(std::vector &coefficients, + const number factor) + { + for (typename std::vector::iterator c = coefficients.begin(); + c != coefficients.end(); ++c) + *c *= factor; + } -template -template -void -Polynomial::shift(std::vector& coefficients, - const number2 offset) -{ + template + template + void + Polynomial::shift(std::vector& coefficients, + const number2 offset) + { #ifdef DEAL_II_LONG_DOUBLE_LOOP_BUG - AssertThrow (false, - ExcMessage("Sorry, but the compiler you are using has a bug that disallows " - "compilation of this function, so you cannot use it. Read more " - "about the bug and when it occurs in the aclocal.m4 file in the " - "top level directory (watch for the string " - "DEAL_II_LONG_DOUBLE_LOOP_BUG)")); - // calm down warning for unused - // args. note that this code is - // actually unreachable - coefficients[0] = offset; + AssertThrow (false, + ExcMessage("Sorry, but the compiler you are using has a bug that disallows " + "compilation of this function, so you cannot use it. Read more " + "about the bug and when it occurs in the aclocal.m4 file in the " + "top level directory (watch for the string " + "DEAL_II_LONG_DOUBLE_LOOP_BUG)")); + // calm down warning for unused + // args. note that this code is + // actually unreachable + coefficients[0] = offset; #else - // Copy coefficients to a vector of - // accuracy given by the argument - std::vector new_coefficients(coefficients.begin(), - coefficients.end()); + // Copy coefficients to a vector of + // accuracy given by the argument + std::vector new_coefficients(coefficients.begin(), + coefficients.end()); - // Traverse all coefficients from - // c_1. c_0 will be modified by - // higher degrees, only. - for (unsigned int d=1; d -template -void -Polynomial::shift(const number2 offset) -{ - shift(coefficients, offset); -} + template + template + void + Polynomial::shift(const number2 offset) + { + shift(coefficients, offset); + } -template -void -Polynomial::print(std::ostream& out) const -{ - for (int i=degree();i>=0;--i) - { - out << coefficients[i] << " x^" << i << std::endl; - } -} + template + void + Polynomial::print(std::ostream& out) const + { + for (int i=degree();i>=0;--i) + { + out << coefficients[i] << " x^" << i << std::endl; + } + } // ------------------ class LagrangeEquidistant --------------- // -LagrangeEquidistant::LagrangeEquidistant (const unsigned int n, - const unsigned int support_point): - Polynomial(compute_coefficients(n,support_point)) -{} + LagrangeEquidistant::LagrangeEquidistant (const unsigned int n, + const unsigned int support_point): + Polynomial(compute_coefficients(n,support_point)) + {} -std::vector -LagrangeEquidistant::compute_coefficients (const unsigned int n, - const unsigned int support_point) -{ - std::vector a (n+1); - Assert(support_point + LagrangeEquidistant::compute_coefficients (const unsigned int n, + const unsigned int support_point) + { + std::vector a (n+1); + Assert(support_point > -LagrangeEquidistant:: -generate_complete_basis (const unsigned int degree) -{ - if (degree==0) - // create constant polynomial - return std::vector > - (1, Polynomial (std::vector (1,1.))); - else - { - // create array of Lagrange - // polynomials - std::vector > v; - for (unsigned int i=0; i<=degree; ++i) - v.push_back(LagrangeEquidistant(degree,i)); - return v; - }; -}; + return a; + } + + + std::vector > + LagrangeEquidistant:: + generate_complete_basis (const unsigned int degree) + { + if (degree==0) + // create constant polynomial + return std::vector > + (1, Polynomial (std::vector (1,1.))); + else + { + // create array of Lagrange + // polynomials + std::vector > v; + for (unsigned int i=0; i<=degree; ++i) + v.push_back(LagrangeEquidistant(degree,i)); + return v; + }; + }; @@ -507,14 +510,14 @@ generate_complete_basis (const unsigned int degree) // Reserve space for polynomials up to degree 19. Should be sufficient // for the start. -template -std::vector *> -Legendre::recursive_coefficients( - 20, static_cast*>(0)); -template -std::vector *> -Legendre::shifted_coefficients( - 20, static_cast*>(0)); + template + std::vector *> + Legendre::recursive_coefficients( + 20, static_cast*>(0)); + template + std::vector *> + Legendre::shifted_coefficients( + 20, static_cast*>(0)); @@ -525,169 +528,174 @@ Legendre::shifted_coefficients( #define SHIFT_TYPE long double #endif -template -void -Legendre::compute_coefficients (const unsigned int k_) -{ - unsigned int k = k_; - - // first make sure that no other - // thread intercepts the operation - // of this function - coefficients_lock.acquire (); - - // 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] == 0))) - // no, then generate the - // respective coefficients - { - recursive_coefficients.resize (k+1, 0); + template + void + Legendre::compute_coefficients (const unsigned int k_) + { + unsigned int k = k_; + + // first make sure that no other + // thread intercepts the operation + // of this function + coefficients_lock.acquire (); + + // 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] == 0))) + // no, then generate the + // respective coefficients + { + recursive_coefficients.resize (k+1, 0); - if (k<=1) - { - // create coefficients - // vectors for k=0 and k=1 - // - // allocate the respective - // amount of memory and - // 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 - recursive_coefficients[0] = c0; - recursive_coefficients[1] = c1; - // Compute polynomials - // orthogonal on [0,1] - c0 = new std::vector(*c0); - c1 = new std::vector(*c1); + if (k<=1) + { + // create coefficients + // vectors for k=0 and k=1 + // + // allocate the respective + // amount of memory and + // 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 + recursive_coefficients[0] = c0; + recursive_coefficients[1] = c1; + // Compute polynomials + // orthogonal on [0,1] + c0 = new std::vector(*c0); + c1 = new std::vector(*c1); - Polynomial::shift(*c0, (SHIFT_TYPE) -1.); - Polynomial::scale(*c0, 2.); - Polynomial::shift(*c1, (SHIFT_TYPE) -1.); - Polynomial::scale(*c1, 2.); - Polynomial::multiply(*c1, std::sqrt(3.)); - shifted_coefficients[0]=c0; - shifted_coefficients[1]=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); + Polynomial::shift(*c0, (SHIFT_TYPE) -1.); + Polynomial::scale(*c0, 2.); + Polynomial::shift(*c1, (SHIFT_TYPE) -1.); + Polynomial::scale(*c1, 2.); + Polynomial::multiply(*c1, std::sqrt(3.)); + shifted_coefficients[0]=c0; + shifted_coefficients[1]=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 number a = 1./(k); - const number b = a*(2*k-1); - const number c = a*(k-1); + const number a = 1./(k); + const number b = a*(2*k-1); + const number c = a*(k-1); - (*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]; - - // finally assign the newly - // created vector to the - // const pointer in the - // coefficients array - recursive_coefficients[k] = ck; - // and compute the - // coefficients for [0,1] - ck = new std::vector(*ck); - shift(*ck,(SHIFT_TYPE) -1.); - Polynomial::scale(*ck, 2.); - Polynomial::multiply(*ck, std::sqrt(2.*k+1.)); - shifted_coefficients[k] = ck; - }; - }; - - // now, everything is done, so - // release the lock again - coefficients_lock.release (); + (*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]; + + // finally assign the newly + // created vector to the + // const pointer in the + // coefficients array + recursive_coefficients[k] = ck; + // and compute the + // coefficients for [0,1] + ck = new std::vector(*ck); + shift(*ck,(SHIFT_TYPE) -1.); + Polynomial::scale(*ck, 2.); + Polynomial::multiply(*ck, std::sqrt(2.*k+1.)); + shifted_coefficients[k] = ck; + }; + }; + + // now, everything is done, so + // release the lock again + coefficients_lock.release (); + } + + + + template + const std::vector & + Legendre::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 + coefficients_lock.acquire (); + const std::vector *p = shifted_coefficients[k]; + coefficients_lock.release (); + + // return the object pointed + // to. since this object does not + // change any more once computed, + // this is MT safe + return *p; + } + + + + template + Legendre::Legendre (const unsigned int k) + : + Polynomial (get_coefficients(k)) + {} + + + + template + std::vector > + Legendre::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 (Legendre(i)); + return v; + }; + } +// ------------------ explicit instantiations --------------- // -template -const std::vector & -Legendre::get_coefficients (const unsigned int k) +namespace Polynomials { - // 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 - coefficients_lock.acquire (); - const std::vector *p = shifted_coefficients[k]; - coefficients_lock.release (); - - // return the object pointed - // to. since this object does not - // change any more once computed, - // this is MT safe - return *p; + template class Polynomial; + template class Polynomial; + template class Polynomial; + + template void Polynomial::shift(const float offset); + template void Polynomial::shift(const double offset); + template void Polynomial::shift(const double offset); + template void Polynomial::shift(const long double offset); + template void Polynomial::shift(const long double offset); + template void Polynomial::shift(const long double offset); + + template class Legendre; } - - - -template -Legendre::Legendre (const unsigned int k) - : - Polynomial (get_coefficients(k)) -{} - - - -template -std::vector > -Legendre::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 (Legendre(i)); - return v; -}; - - -// ------------------ explicit instantiations --------------- // - -template class Polynomial; -template class Polynomial; -template class Polynomial; - -template void Polynomial::shift(const float offset); -template void Polynomial::shift(const double offset); -template void Polynomial::shift(const double offset); -template void Polynomial::shift(const long double offset); -template void Polynomial::shift(const long double offset); -template void Polynomial::shift(const long double offset); - -template class Legendre; diff --git a/deal.II/deal.II/include/fe/mapping_q.h b/deal.II/deal.II/include/fe/mapping_q.h index dbc26e4a1e..8757e60985 100644 --- a/deal.II/deal.II/include/fe/mapping_q.h +++ b/deal.II/deal.II/include/fe/mapping_q.h @@ -19,7 +19,6 @@ #include template class TensorProductPolynomials; -class LagrangeEquidistant; diff --git a/deal.II/deal.II/source/fe/fe_dgp.cc b/deal.II/deal.II/source/fe/fe_dgp.cc index 7018ccf84c..8d6e8f4540 100644 --- a/deal.II/deal.II/source/fe/fe_dgp.cc +++ b/deal.II/deal.II/source/fe/fe_dgp.cc @@ -30,7 +30,7 @@ FE_DGP::FE_DGP (const unsigned int degree) std::vector >(FiniteElementData(get_dpo_vector(degree),1).dofs_per_cell, std::vector(1,true))), degree(degree), - polynomial_space (Legendre::generate_complete_basis(degree)) + polynomial_space (Polynomials::Legendre::generate_complete_basis(degree)) { // if defined, copy over matrices // from precomputed arrays diff --git a/deal.II/deal.II/source/fe/fe_dgp_nonparametric.cc b/deal.II/deal.II/source/fe/fe_dgp_nonparametric.cc index 2dd287b739..b6367ac44b 100644 --- a/deal.II/deal.II/source/fe/fe_dgp_nonparametric.cc +++ b/deal.II/deal.II/source/fe/fe_dgp_nonparametric.cc @@ -30,7 +30,7 @@ FE_DGPNonparametric::FE_DGPNonparametric (const unsigned int degree) std::vector >(FiniteElementData(get_dpo_vector(degree),1).dofs_per_cell, std::vector(1,true))), degree(degree), - polynomial_space (Legendre::generate_complete_basis(degree)) + polynomial_space (Polynomials::Legendre::generate_complete_basis(degree)) { // if defined, copy over matrices // from precomputed arrays diff --git a/deal.II/deal.II/source/fe/fe_dgq.cc b/deal.II/deal.II/source/fe/fe_dgq.cc index 9f7f38eb8a..1702fbb5c6 100644 --- a/deal.II/deal.II/source/fe/fe_dgq.cc +++ b/deal.II/deal.II/source/fe/fe_dgq.cc @@ -32,8 +32,8 @@ FE_DGQ::FE_DGQ (const unsigned int degree) true), std::vector >(FiniteElementData(get_dpo_vector(degree),1).dofs_per_cell, std::vector(1,true))), - degree(degree), - polynomial_space (LagrangeEquidistant::generate_complete_basis(degree)) + degree(degree), + polynomial_space (Polynomials::LagrangeEquidistant::generate_complete_basis(degree)) { // generate permutation/rotation // index sets to generate some diff --git a/deal.II/deal.II/source/fe/fe_q.cc b/deal.II/deal.II/source/fe/fe_q.cc index 450063a9e5..6bcdfbbf19 100644 --- a/deal.II/deal.II/source/fe/fe_q.cc +++ b/deal.II/deal.II/source/fe/fe_q.cc @@ -31,11 +31,11 @@ FE_Q::FE_Q (const unsigned int degree) false), std::vector >(FiniteElementData(get_dpo_vector(degree),1).dofs_per_cell, std::vector(1,true))), - degree(degree), - renumber(this->dofs_per_cell, 0), - renumber_inverse(this->dofs_per_cell, 0), - face_renumber(this->dofs_per_face, 0), - polynomial_space(LagrangeEquidistant::generate_complete_basis(degree)) + degree(degree), + renumber(this->dofs_per_cell, 0), + renumber_inverse(this->dofs_per_cell, 0), + face_renumber(this->dofs_per_face, 0), + polynomial_space(Polynomials::LagrangeEquidistant::generate_complete_basis(degree)) { // do some internal book-keeping on // cells and faces. if in 1d, the diff --git a/deal.II/deal.II/source/fe/mapping_q.cc b/deal.II/deal.II/source/fe/mapping_q.cc index b0d5ac2ab7..f90805632e 100644 --- a/deal.II/deal.II/source/fe/mapping_q.cc +++ b/deal.II/deal.II/source/fe/mapping_q.cc @@ -13,6 +13,7 @@ #include #include +#include #include #include #include @@ -112,9 +113,9 @@ MappingQ::MappingQ (const unsigned int p) // polynomials used as shape // functions for the Qp mapping of // cells at the boundary. - std::vector v; + std::vector v; for (unsigned int i=0; i<=degree; ++i) - v.push_back(LagrangeEquidistant(degree,i)); + v.push_back(Polynomials::LagrangeEquidistant(degree,i)); tensor_pols = new TensorProductPolynomials (v); Assert (n_shape_functions==tensor_pols->n(), diff --git a/deal.II/doc/news/2002/c-3-4.html b/deal.II/doc/news/2002/c-3-4.html index a43abbcc52..774be5c151 100644 --- a/deal.II/doc/news/2002/c-3-4.html +++ b/deal.II/doc/news/2002/c-3-4.html @@ -167,6 +167,13 @@ contributor's names are abbreviated by WB (Wolfgang Bangerth), GK

base

    +
  1. Changed: Because they became too many, the classes describing 1d + polynomials are now in a namespace + Polynomials. +
    + (WB 2002/10/14) +

    +
  2. Changed: When an exception is thrown but not caught in a sub-thread, this exception is not passed to the main thread by the operating system. Rather, if the exception is not caught from the function that diff --git a/tests/base/polynomial1d.cc b/tests/base/polynomial1d.cc index 02d8e2e2b3..3a56a68eb0 100644 --- a/tests/base/polynomial1d.cc +++ b/tests/base/polynomial1d.cc @@ -20,6 +20,9 @@ #include +using namespace Polynomials; + + double scalar_product (const Polynomial& p1, const Polynomial& p2) { diff --git a/tests/base/polynomial_test.cc b/tests/base/polynomial_test.cc index 82e00eb057..585e659ef0 100644 --- a/tests/base/polynomial_test.cc +++ b/tests/base/polynomial_test.cc @@ -20,6 +20,9 @@ #include +using namespace Polynomials; + + template void check_poly(const Point& x, const POLY& p) -- 2.39.5