From 8298fa7e3a52fe940ea13b14975b4c3d7c924812 Mon Sep 17 00:00:00 2001 From: Markus Buerg Date: Wed, 15 Sep 2010 19:58:41 +0000 Subject: [PATCH] Part 2 git-svn-id: https://svn.dealii.org/trunk@21991 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/polynomial.h | 285 +++++++++++++------------ 1 file changed, 148 insertions(+), 137 deletions(-) diff --git a/deal.II/base/include/base/polynomial.h b/deal.II/base/include/base/polynomial.h index 7c5889558b..274c311cb8 100644 --- a/deal.II/base/include/base/polynomial.h +++ b/deal.II/base/include/base/polynomial.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2010 by the deal.II authors +// Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2009, 2010 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -18,6 +18,7 @@ #include #include #include +#include #include @@ -68,12 +69,18 @@ namespace Polynomials */ Polynomial (const std::vector &coefficients); - /** - * Default constructor creating - * an illegal object. - */ + /** + * Constructor creating a zero + * polynomial of degree @p n. + */ + Polynomial (const unsigned int n); + + /** + * Default constructor creating + * an illegal object. + */ Polynomial (); - + /** * Return the value of this * polynomial at the given point. @@ -83,7 +90,7 @@ namespace Polynomials * of the evaluation. */ number value (const number x) const; - + /** * Return the values and the * derivatives of the @@ -159,40 +166,40 @@ namespace Polynomials template void shift (const number2 offset); - /** - * Compute the derivative of a - * polynomial. - */ + /** + * Compute the derivative of a + * polynomial. + */ Polynomial derivative () const; - /** - * Compute the primitive of a - * polynomial. the coefficient - * of the zero order term of - * the polynomial is zero. - */ + /** + * Compute the primitive of a + * polynomial. the coefficient + * of the zero order term of + * the polynomial is zero. + */ Polynomial primitive () const; - /** - * Multiply with a scalar. - */ + /** + * Multiply with a scalar. + */ Polynomial& operator *= (const double s); - /** - * Multiply with another polynomial. - */ + /** + * Multiply with another polynomial. + */ Polynomial& operator *= (const Polynomial& p); - /** - * Add a second polynomial. - */ + /** + * Add a second polynomial. + */ Polynomial& operator += (const Polynomial& p); - - /** - * Subtract a second polynomial. - */ + + /** + * Subtract a second polynomial. + */ Polynomial& operator -= (const Polynomial& p); - + /** * Print coefficients. */ @@ -220,7 +227,7 @@ namespace Polynomials */ static void multiply (std::vector& coefficients, const number factor); - + /** * Coefficients of the polynomial * $\sum_i a_i x^i$. This vector @@ -248,14 +255,14 @@ namespace Polynomials public Polynomial { public: - /** - * Constructor, taking the - * degree of the monomial and - * an optional coefficient as - * arguments. - */ + /** + * Constructor, taking the + * degree of the monomial and + * an optional coefficient as + * arguments. + */ Monomial(const unsigned int n, - const double coefficient = 1.); + const double coefficient = 1.); /** * Return a vector of Monomial @@ -272,15 +279,15 @@ namespace Polynomials static std::vector > generate_complete_basis (const unsigned int degree); - + private: - /** - * Needed by constructor. - */ + /** + * Needed by constructor. + */ static std::vector make_vector(unsigned int n, - const double coefficient); + const double coefficient); }; - + /** * Lagrange polynomials with equidistant interpolation points in @@ -334,7 +341,7 @@ namespace Polynomials static std::vector > generate_complete_basis (const unsigned int degree); - + private: /** @@ -345,10 +352,11 @@ namespace Polynomials * called in the * constructor. */ - static - std::vector + static + void compute_coefficients (const unsigned int n, - const unsigned int support_point); + const unsigned int support_point, + std::vector& a); }; /** @@ -359,23 +367,23 @@ namespace Polynomials class Lagrange { public: - /** - * Given a set of points, this - * function returns all - * Lagrange polynomials for - * interpolation of these - * points. The number of - * polynomials is equal to the - * number of points and the - * maximum degree is one less. - */ + /** + * Given a set of points, this + * function returns all + * Lagrange polynomials for + * interpolation of these + * points. The number of + * polynomials is equal to the + * number of points and the + * maximum degree is one less. + */ static std::vector > generate_complete_basis (const std::vector >& points); }; - - - + + + /** * Legendre polynomials of arbitrary degree on [0,1]. * @@ -412,30 +420,33 @@ namespace Polynomials static std::vector > generate_complete_basis (const unsigned int degree); - + private: /** * Coefficients for the interval $[0,1]$. */ - static std::vector *> shifted_coefficients; - + 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 + * 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. + * 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 *> recursive_coefficients; - + static std::vector > > recursive_coefficients; + /** * Compute coefficients recursively. */ static void compute_coefficients (const unsigned int p); - + /** * Get coefficients for * constructor. This way, it can @@ -471,52 +482,52 @@ namespace Polynomials class Lobatto : public Polynomial { public: - /** - * Constructor for polynomial of degree - * p. There is an exception - * for p==0, see the general - * documentation. - */ + /** + * Constructor for polynomial of degree + * p. There is an exception + * for p==0, see the general + * documentation. + */ Lobatto (const unsigned int p = 0); - /** - * Return the polynomials with index - * 0 up to - * degree. There is an - * exception for p==0, see the - * general documentation. - */ + /** + * Return the polynomials with index + * 0 up to + * degree. There is an + * exception for p==0, see the + * general documentation. + */ static std::vector > generate_complete_basis (const unsigned int p); private: - /** - * Compute coefficients recursively. - */ + /** + * Compute coefficients recursively. + */ std::vector compute_coefficients (const unsigned int p); }; - + /** * Hierarchical polynomials of arbitrary degree on [0,1]. * - * When Constructing a Hierarchical polynomial of degree p, + * When Constructing a Hierarchical polynomial of degree p, * the coefficients will be computed by a recursion formula. The * coefficients are stored in a static data vector to be available * when needed next time. * - * These hierarchical polynomials are based on those of Demkowicz, Oden, + * These hierarchical polynomials are based on those of Demkowicz, Oden, * Rachowicz, and Hardy (CMAME 77 (1989) 79-112, Sec. 4). The first two - * polynomials are the standard linear shape functions given by + * polynomials are the standard linear shape functions given by * $\phi_{0}(x) = 1 - x$ and $\phi_{1}(x) = x$. For $l \geq 2$ * we use the definitions $\phi_{l}(x) = (2x-1)^l - 1, l = 2,4,6,...$ - * and $\phi_{l}(x) = (2x-1)^l - (2x-1), l = 3,5,7,...$. These satisfy the - * recursion relations $\phi_{l}(x) = (2x-1)\phi_{l-1}, l=3,5,7,...$ and - * $\phi_{l}(x) = (2x-1)\phi_{l-1} + \phi_{2}, l=4,6,8,...$. + * and $\phi_{l}(x) = (2x-1)^l - (2x-1), l = 3,5,7,...$. These satisfy the + * recursion relations $\phi_{l}(x) = (2x-1)\phi_{l-1}, l=3,5,7,...$ and + * $\phi_{l}(x) = (2x-1)\phi_{l-1} + \phi_{2}, l=4,6,8,...$. * - * The degrees of freedom are the values at the vertices and the + * The degrees of freedom are the values at the vertices and the * derivatives at the midpoint. Currently, we do not scale the - * polynomials in any way, although better conditioning of the + * polynomials in any way, although better conditioning of the * element stiffness matrix could possibly be achieved with scaling. * * Calling the constructor with a given index p will generate the @@ -541,49 +552,49 @@ namespace Polynomials { public: /** - * Constructor for polynomial of - * degree p. There is an - * exception for p==0, see - * the general documentation. - */ + * Constructor for polynomial of + * degree p. There is an + * exception for p==0, see + * the general documentation. + */ Hierarchical (const unsigned int p); - /** - * Return a vector of - * Hierarchical polynomial - * objects of degrees zero through - * degree, which then spans - * the full space of polynomials - * up to the given degree. Note - * that there is an exception if - * the given degree equals - * zero, see the general - * documentation of this class. - * - * This function may be - * used to initialize the - * TensorProductPolynomials, - * AnisotropicPolynomials, - * and PolynomialSpace - * classes. - */ + /** + * Return a vector of + * Hierarchical polynomial + * objects of degrees zero through + * degree, which then spans + * the full space of polynomials + * up to the given degree. Note + * that there is an exception if + * the given degree equals + * zero, see the general + * documentation of this class. + * + * This function may be + * used to initialize the + * TensorProductPolynomials, + * AnisotropicPolynomials, + * and PolynomialSpace + * classes. + */ static std::vector > generate_complete_basis (const unsigned int degree); - + private: - /** - * Compute coefficients recursively. - */ + /** + * Compute coefficients recursively. + */ static void compute_coefficients (const unsigned int p); - /** - * Get coefficients for - * constructor. This way, it can - * use the non-standard - * constructor of - * Polynomial. - */ + /** + * Get coefficients for + * constructor. This way, it can + * use the non-standard + * constructor of + * Polynomial. + */ static const std::vector & get_coefficients (const unsigned int p); @@ -595,13 +606,13 @@ namespace Polynomials /* -------------------------- inline functions --------------------- */ -namespace Polynomials +namespace Polynomials { template inline - Polynomial::Polynomial () + Polynomial::Polynomial () {} - + template inline unsigned int -- 2.39.5