// $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
#include <base/config.h>
#include <base/exceptions.h>
#include <base/subscriptor.h>
+#include <base/std_cxx1x/shared_ptr.h>
#include <vector>
*/
Polynomial (const std::vector<number> &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.
* of the evaluation.
*/
number value (const number x) const;
-
+
/**
* Return the values and the
* derivatives of the
template <typename number2>
void shift (const number2 offset);
- /**
- * Compute the derivative of a
- * polynomial.
- */
+ /**
+ * Compute the derivative of a
+ * polynomial.
+ */
Polynomial<number> 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<number> primitive () const;
- /**
- * Multiply with a scalar.
- */
+ /**
+ * Multiply with a scalar.
+ */
Polynomial<number>& operator *= (const double s);
- /**
- * Multiply with another polynomial.
- */
+ /**
+ * Multiply with another polynomial.
+ */
Polynomial<number>& operator *= (const Polynomial<number>& p);
- /**
- * Add a second polynomial.
- */
+ /**
+ * Add a second polynomial.
+ */
Polynomial<number>& operator += (const Polynomial<number>& p);
-
- /**
- * Subtract a second polynomial.
- */
+
+ /**
+ * Subtract a second polynomial.
+ */
Polynomial<number>& operator -= (const Polynomial<number>& p);
-
+
/**
* Print coefficients.
*/
*/
static void multiply (std::vector<number>& coefficients,
const number factor);
-
+
/**
* Coefficients of the polynomial
* $\sum_i a_i x^i$. This vector
public Polynomial<number>
{
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
static
std::vector<Polynomial<number> >
generate_complete_basis (const unsigned int degree);
-
+
private:
- /**
- * Needed by constructor.
- */
+ /**
+ * Needed by constructor.
+ */
static std::vector<number> make_vector(unsigned int n,
- const double coefficient);
+ const double coefficient);
};
-
+
/**
* Lagrange polynomials with equidistant interpolation points in
static
std::vector<Polynomial<double> >
generate_complete_basis (const unsigned int degree);
-
+
private:
/**
* called in the
* constructor.
*/
- static
- std::vector<double>
+ static
+ void
compute_coefficients (const unsigned int n,
- const unsigned int support_point);
+ const unsigned int support_point,
+ std::vector<double>& a);
};
/**
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<Polynomial<double> >
generate_complete_basis (const std::vector<Point<1> >& points);
};
-
-
-
+
+
+
/**
* Legendre polynomials of arbitrary degree on <tt>[0,1]</tt>.
*
static
std::vector<Polynomial<double> >
generate_complete_basis (const unsigned int degree);
-
+
private:
/**
* Coefficients for the interval $[0,1]$.
*/
- static std::vector<const std::vector<double> *> shifted_coefficients;
-
+ static std::vector<std_cxx1x::shared_ptr<const std::vector<double> > > 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<const std::vector<double> *> recursive_coefficients;
-
+ static std::vector<std_cxx1x::shared_ptr<const std::vector<double> > > recursive_coefficients;
+
/**
* Compute coefficients recursively.
*/
static void compute_coefficients (const unsigned int p);
-
+
/**
* Get coefficients for
* constructor. This way, it can
class Lobatto : public Polynomial<double>
{
public:
- /**
- * Constructor for polynomial of degree
- * <tt>p</tt>. There is an exception
- * for <tt>p==0</tt>, see the general
- * documentation.
- */
+ /**
+ * Constructor for polynomial of degree
+ * <tt>p</tt>. There is an exception
+ * for <tt>p==0</tt>, see the general
+ * documentation.
+ */
Lobatto (const unsigned int p = 0);
- /**
- * Return the polynomials with index
- * <tt>0</tt> up to
- * <tt>degree</tt>. There is an
- * exception for <tt>p==0</tt>, see the
- * general documentation.
- */
+ /**
+ * Return the polynomials with index
+ * <tt>0</tt> up to
+ * <tt>degree</tt>. There is an
+ * exception for <tt>p==0</tt>, see the
+ * general documentation.
+ */
static std::vector<Polynomial<double> >
generate_complete_basis (const unsigned int p);
private:
- /**
- * Compute coefficients recursively.
- */
+ /**
+ * Compute coefficients recursively.
+ */
std::vector<double> compute_coefficients (const unsigned int p);
};
-
+
/**
* Hierarchical polynomials of arbitrary degree on <tt>[0,1]</tt>.
*
- * When Constructing a Hierarchical polynomial of degree <tt>p</tt>,
+ * When Constructing a Hierarchical polynomial of degree <tt>p</tt>,
* 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 <tt>p</tt> will generate the
{
public:
/**
- * Constructor for polynomial of
- * degree <tt>p</tt>. There is an
- * exception for <tt>p==0</tt>, see
- * the general documentation.
- */
+ * Constructor for polynomial of
+ * degree <tt>p</tt>. There is an
+ * exception for <tt>p==0</tt>, see
+ * the general documentation.
+ */
Hierarchical (const unsigned int p);
- /**
- * Return a vector of
- * Hierarchical polynomial
- * objects of degrees zero through
- * <tt>degree</tt>, which then spans
- * the full space of polynomials
- * up to the given degree. Note
- * that there is an exception if
- * the given <tt>degree</tt> 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
+ * <tt>degree</tt>, which then spans
+ * the full space of polynomials
+ * up to the given degree. Note
+ * that there is an exception if
+ * the given <tt>degree</tt> 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<Polynomial<double> >
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<double> &
get_coefficients (const unsigned int p);
/* -------------------------- inline functions --------------------- */
-namespace Polynomials
+namespace Polynomials
{
template <typename number>
inline
- Polynomial<number>::Polynomial ()
+ Polynomial<number>::Polynomial ()
{}
-
+
template <typename number>
inline
unsigned int