#include <base/config.h>
#include <base/exceptions.h>
#include <base/subscriptor.h>
+#include <base/std_cxx1x/shared_ptr.h>
#include <vector>
* 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
* Add a second polynomial.
*/
Polynomial<number>& operator += (const Polynomial<number>& p);
-
+
/**
* 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
static
std::vector<Polynomial<number> >
generate_complete_basis (const unsigned int degree);
-
+
private:
/**
* Needed by constructor.
static std::vector<number> make_vector(unsigned int n,
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
+ static
void
compute_coefficients (const unsigned int n,
const unsigned int support_point,
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
/**
* 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
static
std::vector<Polynomial<double> >
generate_complete_basis (const unsigned int degree);
-
+
private:
/**
* Compute coefficients recursively.
*/
static const std::vector<double> &
get_coefficients (const unsigned int p);
-
+
static std::vector<const std::vector<double> *> recursive_coefficients;
- };
+ };
}
/** @} */
/* -------------------------- inline functions --------------------- */
-namespace Polynomials
+namespace Polynomials
{
template <typename number>
inline
- Polynomial<number>::Polynomial ()
+ Polynomial<number>::Polynomial ()
{}
-
+
template <typename number>
inline
unsigned int
// ------------------ class Legendre --------------- //
-//TODO:[?] This class leaks memory, but only at the very end of a program.
-// Since it expands the Legendre<number>::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<const std::vector<double> *>
- Legendre::recursive_coefficients(20,
- static_cast<const std::vector<double>*>(0));
- std::vector<const std::vector<double> *>
- Legendre::shifted_coefficients(20,
- static_cast<const std::vector<double>*>(0));
+ std::vector<std_cxx1x::shared_ptr<const std::vector<double> > >
+ Legendre::recursive_coefficients(20);
+ std::vector<std_cxx1x::shared_ptr<const std::vector<double> > >
+ Legendre::shifted_coefficients(20);
Legendre::Legendre (const unsigned int k)
// no, then generate the
// respective coefficients
{
- recursive_coefficients.resize (k+1, 0);
+ recursive_coefficients.resize (k+1);
if (k<=1)
{
(*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<const std::vector<double> >(c0);
+ recursive_coefficients[1] =
+ std_cxx1x::shared_ptr<const std::vector<double> >(c1);
// Compute polynomials
// orthogonal on [0,1]
Polynomial<double>::shift<SHIFT_TYPE> (*c1, -1.);
Polynomial<double>::scale(*c1, 2.);
Polynomial<double>::multiply(*c1, std::sqrt(3.));
- shifted_coefficients[0]=c0;
- shifted_coefficients[1]=c1;
+ shifted_coefficients[0]=std_cxx1x::shared_ptr<const std::vector<double> >(c0);
+ shifted_coefficients[1]=std_cxx1x::shared_ptr<const std::vector<double> >(c1);
}
else
{
// created vector to the
// const pointer in the
// coefficients array
- recursive_coefficients[k] = ck;
+ recursive_coefficients[k] =
+ std_cxx1x::shared_ptr<const std::vector<double> >(ck);
// and compute the
// coefficients for [0,1]
ck = new std::vector<double>(*ck);
Polynomial<double>::shift<SHIFT_TYPE> (*ck, -1.);
Polynomial<double>::scale(*ck, 2.);
Polynomial<double>::multiply(*ck, std::sqrt(2.*k+1.));
- shifted_coefficients[k] = ck;
+ shifted_coefficients[k] =
+ std_cxx1x::shared_ptr<const std::vector<double> >(ck);
};
};
}