From: guido Date: Mon, 28 May 2001 11:27:22 +0000 (+0000) Subject: Legendre polynomials fixed X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=813096040b2a2ea71d9f7764135b479df8140f21;p=dealii-svn.git Legendre polynomials fixed git-svn-id: https://svn.dealii.org/trunk@4734 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/include/base/polynomial.h b/deal.II/base/include/base/polynomial.h index 5153321895..b520904bbc 100644 --- a/deal.II/base/include/base/polynomial.h +++ b/deal.II/base/include/base/polynomial.h @@ -43,7 +43,7 @@ class Polynomial : public Subscriptor * of the array denotes the * constant term, the second the * linear one, and so on. The - * order of the polynomial + * degree of the polynomial * represented by this object is * thus the number of elements in * the @p{coefficient} array @@ -85,6 +85,17 @@ class Polynomial : public Subscriptor 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; + /** * Exception */ @@ -214,5 +225,14 @@ class Legendre : public Polynomial }; +template +inline +unsigned int +Polynomial::degree () const +{ + Assert (coefficients.size()>0, ExcVoidPolynomial()); + return coefficients.size() - 1; +} + #endif diff --git a/deal.II/base/source/legendre.cc b/deal.II/base/source/legendre.cc index cc3c975d44..81c74a4a41 100644 --- a/deal.II/base/source/legendre.cc +++ b/deal.II/base/source/legendre.cc @@ -91,15 +91,17 @@ Legendre::compute_coefficients (const unsigned int k) std::vector *ck = new std::vector(k+1); - const number a = 1./k+1; - const number b = a*(2*k+1); + const number a = 1./(k); + const number b = a*(2*k-1); + const number c = a*(k-1); (*ck)[k] = b*(*coefficients[k-1])[k-1]; (*ck)[k-1] = b*(*coefficients[k-1])[k-2]; for (unsigned int i=1 ; i<= k-2 ; ++i) (*ck)[i] = b*(*coefficients[k-1])[i-1] - - k*a*(*coefficients[k-2])[i]; - (*ck)[0] = -k*a*(*coefficients[k-2])[0]; + -c*(*coefficients[k-2])[i]; + + (*ck)[0] = -c*(*coefficients[k-2])[0]; // finally assign the newly // created vector to the