From: kronbichler Date: Wed, 3 Feb 2010 15:59:43 +0000 (+0000) Subject: Fix a memory leak. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=eb4481d215f3d36fce71ab3fab6e564733bc0cc6;p=dealii-svn.git Fix a memory leak. git-svn-id: https://svn.dealii.org/trunk@20489 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/include/base/polynomial.h b/deal.II/base/include/base/polynomial.h index 8eda982788..9a840f0cad 100644 --- a/deal.II/base/include/base/polynomial.h +++ b/deal.II/base/include/base/polynomial.h @@ -404,6 +404,11 @@ namespace Polynomials */ Legendre (const unsigned int p); + /** + * Destructor. + */ + ~Legendre(); + /** * Return a vector of Legendre * polynomial objects of degrees diff --git a/deal.II/base/source/polynomial.cc b/deal.II/base/source/polynomial.cc index e7611a19b2..d4259788ca 100644 --- a/deal.II/base/source/polynomial.cc +++ b/deal.II/base/source/polynomial.cc @@ -1,5 +1,5 @@ //--------------------------------------------------------------------------- -// $Id$ +// $Id$ // Version: $Name$ // // Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2009 by the deal.II authors @@ -31,7 +31,7 @@ DEAL_II_NAMESPACE_OPEN // to be a problem since we only need it on very rare occasions. if // someone finds this is a bottleneck, feel free to replace it by a // more fine-grained solution -namespace +namespace { Threads::ThreadMutex coefficients_lock; } @@ -68,8 +68,8 @@ namespace Polynomials Assert (coefficients.size() > 0, ExcEmptyObject()); Assert (values.size() > 0, ExcZero()); const unsigned int values_size=values.size(); - - + + // if we only need the value, then // call the other function since // that is significantly faster @@ -124,7 +124,7 @@ namespace Polynomials { *c *= f; f *= factor; - } + } } @@ -160,7 +160,7 @@ namespace Polynomials return *this; } - + template Polynomial& Polynomial::operator *= (const Polynomial& p) @@ -169,16 +169,16 @@ namespace Polynomials unsigned int new_degree = this->degree() + p.degree(); std::vector new_coefficients(new_degree+1, 0.); - + for (unsigned int i=0; icoefficients.size(); ++j) new_coefficients[i+j] += this->coefficients[j]*p.coefficients[i]; this->coefficients = new_coefficients; - + return *this; } - + template Polynomial& Polynomial::operator += (const Polynomial& p) @@ -194,7 +194,7 @@ namespace Polynomials return *this; } - + template Polynomial& Polynomial::operator -= (const Polynomial& p) @@ -210,13 +210,13 @@ namespace Polynomials return *this; } - + 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 " @@ -228,12 +228,12 @@ namespace Polynomials // args. note that this code is // actually unreachable coefficients[0] = offset; -#else +#else // 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. @@ -250,7 +250,7 @@ namespace Polynomials // needed and computed // successively. number2 offset_power = offset; - + // Compute (x+offset)^d // and modify all values c_k // with k Polynomial Polynomial::derivative () const @@ -305,7 +305,7 @@ namespace Polynomials return Polynomial (newcoefficients); } - + template Polynomial @@ -318,7 +318,7 @@ namespace Polynomials return Polynomial (newcoefficients); } - + template void @@ -343,15 +343,15 @@ namespace Polynomials result[n] = coefficient; return result; } - - + + template Monomial::Monomial (unsigned int n, double coefficient) : Polynomial(make_vector(n, coefficient)) {} - - + + template std::vector > Monomial::generate_complete_basis (const unsigned int degree) @@ -378,7 +378,7 @@ namespace Polynomials // 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.; @@ -387,7 +387,7 @@ namespace Polynomials // polynomial as the product // of linear factors std::vector two (2, 1.); - + for (unsigned int k=0;k<=n;++k) { if (k != support_point) @@ -397,10 +397,10 @@ namespace Polynomials factor.scale(1.*n/(support_point - k)); (*this) *= factor; } - } + } } } - + void LagrangeEquidistant::compute_coefficients (const unsigned int n, @@ -413,7 +413,7 @@ namespace Polynomials Assert(support_point > Lagrange::generate_complete_basis (const std::vector >& points) { @@ -657,7 +657,7 @@ namespace Polynomials std::vector linear(2, 1.); // We start with a constant polynomial std::vector one(1, 1.); - + for (unsigned int i=0;i *c0 = 0; + std::swap (recursive_coefficients[i], c0); + delete c0; + } + for (unsigned int i=0; i *c0 = 0; + std::swap (shifted_coefficients[i], c0); + delete c0; + } + } + + + void Legendre::compute_coefficients (const unsigned int k_) { @@ -740,7 +763,7 @@ namespace Polynomials #else typedef long double SHIFT_TYPE; #endif - + unsigned int k = k_; // first make sure that no other @@ -762,7 +785,7 @@ namespace Polynomials // respective coefficients { recursive_coefficients.resize (k+1, 0); - + if (k<=1) { // create coefficients @@ -782,20 +805,21 @@ namespace Polynomials // now make these arrays // const - recursive_coefficients[0] = c0; - recursive_coefficients[1] = c1; + 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, -1.); - Polynomial::scale(*c0, 2.); - Polynomial::shift (*c1, -1.); - Polynomial::scale(*c1, 2.); - Polynomial::multiply(*c1, std::sqrt(3.)); - shifted_coefficients[0]=c0; - shifted_coefficients[1]=c1; + std::vector *d0 = new std::vector(*c0); + std::vector *d1 = new std::vector(*c1); + + Polynomial::shift (*d0, -1.); + Polynomial::scale(*d0, 2.); + Polynomial::shift (*d1, -1.); + Polynomial::scale(*d1, 2.); + Polynomial::multiply(*d1, std::sqrt(3.)); + shifted_coefficients[0]=d0; + shifted_coefficients[1]=d1; } else { @@ -812,11 +836,11 @@ namespace Polynomials coefficients_lock.acquire (); std::vector *ck = new std::vector(k+1); - + const double a = 1./(k); const double b = a*(2*k-1); const double 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) @@ -902,20 +926,20 @@ namespace Polynomials // until we quit this function Threads::ThreadMutex::ScopedLock lock(coefficients_lock); - // The first 2 coefficients + // 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.size() >= k+1) && (recursive_coefficients[k] == 0)) ) // no, then generate the // respective coefficients { recursive_coefficients.resize (k+1, 0); - + if (k<=1) { // create coefficients @@ -952,7 +976,7 @@ namespace Polynomials (*c2)[0] = 0.*a; (*c2)[1] = -4.*a; (*c2)[2] = 4.*a; - + recursive_coefficients[2] = c2; } else @@ -970,15 +994,15 @@ namespace Polynomials coefficients_lock.acquire (); std::vector *ck = new std::vector(k+1); - + const double a = 1.; //1./(2.*k); (*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] ); - + (*ck)[k] = a*2.*(*recursive_coefficients[k-1])[k-1]; // for even degrees, we need // to add a multiple of @@ -991,7 +1015,7 @@ namespace Polynomials (*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