Clarify the docs of the hierarchical polynomials, and change what generate_complete_b...

git-svn-id: https://svn.dealii.org/trunk@7414 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/base/include/base/polynomial.h b/deal.II/base/include/base/polynomial.h
index a5bbe5f7f710e3c24648cec3318e6e07f1b14022..ee6f29dd397fedce158c7104909d888dbba65b35 100644 (file)
--- a/deal.II/base/include/base/polynomial.h
+++ b/deal.II/base/include/base/polynomial.h
@@ -239,7 +239,7 @@ namespace Polynomials
                                         * space of polynomials up to the
                                         * given degree. The polynomials
                                         * are generated by calling the
-                                        * destructor of this class with
+                                        * constructor of this class with
                                         * the same degree but support
                                         * point running from zero to
                                         * @p{degree}. This function may
@@ -366,6 +366,22 @@ namespace Polynomials
  * 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{p} will generate the
+ * following: if @p{p==0}, then the resulting polynomial is the linear
+ * function associated with the left vertex, if @p{p==1} the one
+ * associated with the right vertex. For higher values of @p{p}, you
+ * get the polynomial of degree @p{p} that is orthogonal to all
+ * previous ones. Note that for @p{p==0} you therefore do @em{not}
+ * get a polynomial of degree zero, but one of degree one. This is to
+ * allow generating a complete basis for polynomial spaces, by just
+ * iterating over the indices given to the constructor.
+ *
+ * On the other hand, the function @p{generate_complete_basis} creates
+ * a complete basis of given degree. In order to be consistent with
+ * the concept of a polynomial degree, if the given argument is zero,
+ * it does @em{not} return the linear polynomial described above, but
+ * rather a constant polynomial.
+ *
  * @author Brian Carnes, 2002
  */
   template <typename number>
@@ -374,19 +390,28 @@ namespace Polynomials
     public:
                                      /**
                                      * Constructor for polynomial of
-                                     * order @p{p}.
+                                     * order @p{p}. There is an
+                                     * exception for @p{p==0}, see
+                                     * the general documentation.
                                      */
       Hierarchical (const unsigned int p);
 
                                     /**
-                                     * Return a vector of Hierarchical
-                                     * polynomial objects of orders
-                                     * zero through @p{degree}, which
-                                     * then spans the full space of
-                                     * polynomials up to the given
-                                     * degree. This function may be
+                                     * Return a vector of
+                                     * Hierarchical polynomial
+                                     * objects of orders zero through
+                                     * @p{degree}, which then spans
+                                     * the full space of polynomials
+                                     * up to the given degree. Note
+                                     * that there is an exception if
+                                     * the given @p{degree} equals
+                                     * zero, see the general
+                                     * documentation of this class.
+                                     *
+                                     * This function may be
                                      * used to initialize the
-                                     * @ref{TensorProductPolynomials}
+                                     * @ref{TensorProductPolynomials},
+                                     * @ref{AnisotropicPolynomials},
                                      * and @ref{PolynomialSpace}
                                      * classes.
                                      */

diff --git a/deal.II/base/source/polynomial.cc b/deal.II/base/source/polynomial.cc
index 38ac842782dd212426e080fdad46d286b5768038..373485c3c9697bd036f995fe40bce336b87ab2b2 100644 (file)
--- a/deal.II/base/source/polynomial.cc
+++ b/deal.II/base/source/polynomial.cc
@@ -34,6 +34,7 @@ namespace
 }
 
 
+
 namespace Polynomials
 {
 
@@ -41,7 +42,8 @@ namespace Polynomials
 
 
   template <typename number>
-  Polynomial<number>::Polynomial (const std::vector<number> &a):
+  Polynomial<number>::Polynomial (const std::vector<number> &a)
+                  :
                   coefficients(a)
   {}
 
@@ -119,8 +121,8 @@ namespace Polynomials
 
   template <typename number>
   void
-  Polynomial<number>::scale(std::vector<number>& coefficients,
-                            const number factor)
+  Polynomial<number>::scale(std::vector<number> &coefficients,
+                            const number         factor)
   {
     double f = 1.;
     for (typename std::vector<number>::iterator c = coefficients.begin();
@@ -841,11 +843,26 @@ namespace Polynomials
   std::vector<Polynomial<number> >
   Hierarchical<number>::generate_complete_basis (const unsigned int degree)
   {
-    std::vector<Polynomial<double> > v;
-    v.reserve(degree+1);
-    for (unsigned int i=0; i<=degree; ++i)
-      v.push_back (Hierarchical<double>(i));
-    return v;
+    if (degree==0)
+                                       // create constant
+                                       // polynomial. note that we
+                                       // can't use the other branch
+                                       // of the if-statement, since
+                                       // calling the constructor of
+                                       // this class with argument
+                                       // zero does _not_ create the
+                                       // constant polynomial, but
+                                       // rather 1-x
+      return std::vector<Polynomial<double> >
+        (1, Polynomial<double> (std::vector<double> (1,1.)));
+    else
+      {
+        std::vector<Polynomial<double> > v;
+        v.reserve(degree+1);
+        for (unsigned int i=0; i<=degree; ++i)
+          v.push_back (Hierarchical<double>(i));
+        return v;
+      }
   }
 }