]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Part 2
authorbuerg <buerg@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 15 Sep 2010 19:58:41 +0000 (19:58 +0000)
committerbuerg <buerg@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 15 Sep 2010 19:58:41 +0000 (19:58 +0000)
git-svn-id: https://svn.dealii.org/trunk@21991 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/polynomial.h

index 7c5889558b106da046590b33da5f820f40f01346..274c311cb8b926bb1c19572f24577551dddf81db 100644 (file)
@@ -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 <base/config.h>
 #include <base/exceptions.h>
 #include <base/subscriptor.h>
+#include <base/std_cxx1x/shared_ptr.h>
 
 #include <vector>
 
@@ -68,12 +69,18 @@ namespace Polynomials
                                         */
       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.
@@ -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 <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.
                                         */
@@ -220,7 +227,7 @@ namespace Polynomials
                                         */
       static void multiply (std::vector<number>& coefficients,
                             const number factor);
-    
+   
                                        /**
                                         * Coefficients of the polynomial
                                         * $\sum_i a_i x^i$. This vector
@@ -248,14 +255,14 @@ namespace Polynomials
        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
@@ -272,15 +279,15 @@ namespace Polynomials
       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
@@ -334,7 +341,7 @@ namespace Polynomials
       static
       std::vector<Polynomial<double> >
       generate_complete_basis (const unsigned int degree);
-    
+   
     private:
 
                                        /**
@@ -345,10 +352,11 @@ namespace Polynomials
                                         * 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);
   };
 
 /**
@@ -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<Polynomial<double> >
       generate_complete_basis (const std::vector<Point<1> >& points);
   };
-  
-  
-  
 /**
  * Legendre polynomials of arbitrary degree on <tt>[0,1]</tt>.
  *
@@ -412,30 +420,33 @@ namespace Polynomials
       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
@@ -471,52 +482,52 @@ namespace Polynomials
   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
@@ -541,49 +552,49 @@ namespace Polynomials
   {
     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);
  
@@ -595,13 +606,13 @@ namespace Polynomials
 
 /* -------------------------- inline functions --------------------- */
 
-namespace Polynomials 
+namespace Polynomials
 {
   template <typename number>
   inline
-  Polynomial<number>::Polynomial () 
+  Polynomial<number>::Polynomial ()
   {}
-  
   template <typename number>
   inline
   unsigned int

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.