From f683b264d52318678baa821e3e1ed2b58aec29bf Mon Sep 17 00:00:00 2001 From: Markus Buerg Date: Wed, 15 Sep 2010 15:00:03 +0000 Subject: [PATCH] Lost patches recovered. git-svn-id: https://svn.dealii.org/trunk@21990 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/source/polynomial.cc | 306 +++++++++++++++++------------- 1 file changed, 171 insertions(+), 135 deletions(-) diff --git a/deal.II/base/source/polynomial.cc b/deal.II/base/source/polynomial.cc index a84e6379b3..3ee279fae3 100644 --- a/deal.II/base/source/polynomial.cc +++ b/deal.II/base/source/polynomial.cc @@ -52,6 +52,14 @@ namespace Polynomials + template + Polynomial::Polynomial (const unsigned int n) + : + coefficients(n+1, 0.) + {} + + + template void Polynomial::value (const number x, @@ -157,7 +165,7 @@ namespace Polynomials Polynomial& Polynomial::operator *= (const Polynomial& p) { - // Degree of the product + // Degree of the product unsigned int new_degree = this->degree() + p.degree(); std::vector new_coefficients(new_degree+1, 0.); @@ -175,8 +183,8 @@ namespace Polynomials Polynomial& Polynomial::operator += (const Polynomial& p) { - // if necessary expand the number - // of coefficients we store + // if necessary expand the number + // of coefficients we store if (p.coefficients.size() > coefficients.size()) coefficients.resize (p.coefficients.size(), 0.); @@ -191,8 +199,8 @@ namespace Polynomials Polynomial& Polynomial::operator -= (const Polynomial& p) { - // if necessary expand the number - // of coefficients we store + // if necessary expand the number + // of coefficients we store if (p.coefficients.size() > coefficients.size()) coefficients.resize (p.coefficients.size(), 0.); @@ -329,7 +337,7 @@ namespace Polynomials template std::vector Monomial::make_vector(unsigned int n, - double coefficient) + double coefficient) { std::vector result(n+1, 0.); result[n] = coefficient; @@ -339,8 +347,8 @@ namespace Polynomials template Monomial::Monomial (unsigned int n, - double coefficient) - : Polynomial(make_vector(n, coefficient)) + double coefficient) + : Polynomial(make_vector(n, coefficient)) {} @@ -358,17 +366,47 @@ namespace Polynomials LagrangeEquidistant::LagrangeEquidistant (const unsigned int n, const unsigned int support_point) - : - Polynomial(compute_coefficients(n,support_point)) - {} + { + if (n <= 10) + { + this->coefficients.resize(n+1); + compute_coefficients(n, support_point, this->coefficients); + } + else + { + // We have precomputed tables + // up to degree 10. For + // higher order, we have to + // compute by hand. + + // Start with the constant one + this->coefficients.resize(1); + this->coefficients[0] = 1.; + // Then compute the Lagrange + // polynomial as the product + // of linear factors + std::vector two (2, 1.); + for (unsigned int k=0;k<=n;++k) + { + if (k != support_point) + { + two[0] = -1.*k/n; + Polynomial factor(two); + factor.scale(1.*n/(support_point - k)); + (*this) *= factor; + } + } + } + } - std::vector + + void LagrangeEquidistant::compute_coefficients (const unsigned int n, - const unsigned int support_point) + const unsigned int support_point, + std::vector& a) { - std::vector a (n+1); Assert(support_point >& points) { std::vector > p(points.size()); - // polynomials are built as - // products of linear - // factors. The coefficient in - // front of the linear term is - // always 1. + // polynomials are built as + // products of linear + // factors. The coefficient in + // front of the linear term is + // always 1. std::vector linear(2, 1.); - // We start with a constant polynomial + // We start with a constant polynomial std::vector one(1, 1.); for (unsigned int i=0;i(one); - for (unsigned int k=0;k factor(linear); - factor *= 1./(points[i](0)-points[k](0)); - p[i] *= factor; - } + // Construct interpolation formula + p[i] = Polynomial(one); + for (unsigned int k=0;k factor(linear); + factor *= 1./(points[i](0)-points[k](0)); + p[i] *= factor; + } } - // Since the previous operation - // is subject to round-off error - // amplification, we correct the - // polynomials here. + // Since the previous operation + // is subject to round-off error + // amplification, we correct the + // polynomials here. for (unsigned int i=0;i q = p[k]; - if (i==k) - { - q *= 1.-value; - p[i] += q; - } - else - { - q *= -value; - p[i] += q; - } - } + for (unsigned int k=0;k q = p[k]; + if (i==k) + { + q *= 1.-value; + p[i] += q; + } + else + { + q *= -value; + p[i] += q; + } + } } return p; @@ -666,22 +702,12 @@ namespace Polynomials // ------------------ class Legendre --------------- // -//TODO:[?] This class leaks memory, but only at the very end of a program. -// Since it expands the Legendre::coefficients array, the elements -// of this static variable are not destroyed at the end of the program -// run. While this is not a problem (since the returned memory could -// not be used anyway then), it is a little confusing when looking at -// a memory checker such as "purify". Maybe, this should be handled somehow -// to avoid this confusion in future. - // Reserve space for polynomials up to degree 19. Should be sufficient // for the start. - std::vector *> - Legendre::recursive_coefficients(20, - static_cast*>(0)); - std::vector *> - Legendre::shifted_coefficients(20, - static_cast*>(0)); + std::vector > > + Legendre::recursive_coefficients(20); + std::vector > > + Legendre::shifted_coefficients(20); Legendre::Legendre (const unsigned int k) @@ -721,11 +747,12 @@ namespace Polynomials // already exist? if ((recursive_coefficients.size() < k+1) || ((recursive_coefficients.size() >= k+1) && - (recursive_coefficients[k] == 0))) + (recursive_coefficients[k] == + std_cxx1x::shared_ptr >()))) // no, then generate the // respective coefficients { - recursive_coefficients.resize (k+1, 0); + recursive_coefficients.resize (k+1); if (k<=1) { @@ -745,9 +772,16 @@ namespace Polynomials (*c1)[1] = 1.; // now make these arrays - // const - recursive_coefficients[0] = c0; - recursive_coefficients[1] = c1; + // 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_cxx1x::shared_ptr >(c0); + recursive_coefficients[1] = + std_cxx1x::shared_ptr >(c1); + // Compute polynomials // orthogonal on [0,1] c0 = new std::vector(*c0); @@ -758,8 +792,8 @@ namespace Polynomials Polynomial::shift (*c1, -1.); Polynomial::scale(*c1, 2.); Polynomial::multiply(*c1, std::sqrt(3.)); - shifted_coefficients[0]=c0; - shifted_coefficients[1]=c1; + shifted_coefficients[0]=std_cxx1x::shared_ptr >(c0); + shifted_coefficients[1]=std_cxx1x::shared_ptr >(c1); } else { @@ -793,14 +827,16 @@ namespace Polynomials // created vector to the // const pointer in the // coefficients array - recursive_coefficients[k] = ck; + recursive_coefficients[k] = + std_cxx1x::shared_ptr >(ck); // and compute the // coefficients for [0,1] ck = new std::vector(*ck); Polynomial::shift (*ck, -1.); Polynomial::scale(*ck, 2.); Polynomial::multiply(*ck, std::sqrt(2.*k+1.)); - shifted_coefficients[k] = ck; + shifted_coefficients[k] = + std_cxx1x::shared_ptr >(ck); }; }; } @@ -948,15 +984,15 @@ std::vector > Lobatto::generate_complete_basis (const unsigne // check: does the information // already exist? if ( (recursive_coefficients.size() < k+1) || - ((recursive_coefficients.size() >= k+1) && + ((recursive_coefficients.size() >= k+1) && (recursive_coefficients[k] == 0)) ) - // no, then generate the - // respective coefficients + // no, then generate the + // respective coefficients { - recursive_coefficients.resize (k+1, 0); + recursive_coefficients.resize (k+1, 0); - if (k<=1) - { + if (k<=1) + { // create coefficients // vectors for k=0 and k=1 // @@ -965,37 +1001,37 @@ std::vector > Lobatto::generate_complete_basis (const unsigne // later assign it to the // coefficients array to // make it const - std::vector *c0 = new std::vector(2); - (*c0)[0] = 1.; - (*c0)[1] = -1.; + std::vector *c0 = new std::vector(2); + (*c0)[0] = 1.; + (*c0)[1] = -1.; - std::vector *c1 = new std::vector(2); - (*c1)[0] = 0.; - (*c1)[1] = 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; - } - else if (k==2) - { + recursive_coefficients[0] = c0; + recursive_coefficients[1] = c1; + } + else if (k==2) + { coefficients_lock.release (); - compute_coefficients(1); - coefficients_lock.acquire (); + compute_coefficients(1); + coefficients_lock.acquire (); - std::vector *c2 = new std::vector(3); + std::vector *c2 = new std::vector(3); - const double a = 1.; //1./8.; + const double a = 1.; //1./8.; - (*c2)[0] = 0.*a; - (*c2)[1] = -4.*a; - (*c2)[2] = 4.*a; + (*c2)[0] = 0.*a; + (*c2)[1] = -4.*a; + (*c2)[2] = 4.*a; - recursive_coefficients[2] = c2; - } - else - { + recursive_coefficients[2] = c2; + } + else + { // for larger numbers, // compute the coefficients // recursively. to do so, @@ -1004,39 +1040,39 @@ std::vector > Lobatto::generate_complete_basis (const unsigne // allow the called // function to acquire it // itself - coefficients_lock.release (); - compute_coefficients(k-1); - coefficients_lock.acquire (); + coefficients_lock.release (); + compute_coefficients(k-1); + coefficients_lock.acquire (); - std::vector *ck = new std::vector(k+1); + std::vector *ck = new std::vector(k+1); - const double a = 1.; //1./(2.*k); + const double a = 1.; //1./(2.*k); - (*ck)[0] = - a*(*recursive_coefficients[k-1])[0]; + (*ck)[0] = - a*(*recursive_coefficients[k-1])[0]; - for (unsigned int i=1; i<=k-1; ++i) - (*ck)[i] = a*( 2.*(*recursive_coefficients[k-1])[i-1] - - (*recursive_coefficients[k-1])[i] ); + for (unsigned int i=1; i<=k-1; ++i) + (*ck)[i] = a*( 2.*(*recursive_coefficients[k-1])[i-1] + - (*recursive_coefficients[k-1])[i] ); - (*ck)[k] = a*2.*(*recursive_coefficients[k-1])[k-1]; - // for even degrees, we need - // to add a multiple of - // basis fcn phi_2 - if ( (k%2) == 0 ) - { - double b = 1.; //8.; - //for (unsigned int i=1; i<=k; i++) - // b /= 2.*i; + (*ck)[k] = a*2.*(*recursive_coefficients[k-1])[k-1]; + // for even degrees, we need + // to add a multiple of + // basis fcn phi_2 + if ( (k%2) == 0 ) + { + double b = 1.; //8.; + //for (unsigned int i=1; i<=k; i++) + // b /= 2.*i; - (*ck)[1] += b*(*recursive_coefficients[2])[1]; - (*ck)[2] += b*(*recursive_coefficients[2])[2]; - } + (*ck)[1] += b*(*recursive_coefficients[2])[1]; + (*ck)[2] += b*(*recursive_coefficients[2])[2]; + } // finally assign the newly // created vector to the // const pointer in the // coefficients array - recursive_coefficients[k] = ck; - }; + recursive_coefficients[k] = ck; + }; }; } @@ -1045,21 +1081,21 @@ std::vector > Lobatto::generate_complete_basis (const unsigne const std::vector & Hierarchical::get_coefficients (const unsigned int k) { - // first make sure the coefficients - // get computed if so necessary + // 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 + // then get a pointer to the array + // of coefficients. do that in a MT // safe way coefficients_lock.acquire (); const std::vector *p = recursive_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 the object pointed + // to. since this object does not + // change any more once computed, + // this is MT safe return *p; } @@ -1113,4 +1149,4 @@ namespace Polynomials template class Monomial; } -DEAL_II_NAMESPACE_CLOSE +DEAL_II_NAMESPACE_CLOSE \ No newline at end of file -- 2.39.5