From: hartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Tue, 9 Mar 2004 14:45:44 +0000 (+0000)
Subject: Implementation of a complete polynomial space of order k based on monomials. Makes... 
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b22e7ea16beacc51054bc6477a960b616c469b0d;p=dealii-svn.git

Implementation of a complete polynomial space of order k based on monomials. Makes use of the PolynomialSpace class and imposes an ordering of the monomials which ensures identical ordering of Subspaces.


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

diff --git a/deal.II/base/include/base/polynomials_p.h b/deal.II/base/include/base/polynomials_p.h
new file mode 100644
index 0000000000..ee110547fc
--- /dev/null
+++ b/deal.II/base/include/base/polynomials_p.h
@@ -0,0 +1,113 @@
+//-------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2004 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//-------------------------------------------------------------------
+#ifndef __deal2__polynomials_P_h
+#define __deal2__polynomials_P_h
+
+
+#include <base/config.h>
+#include <base/exceptions.h>
+#include <base/tensor.h>
+#include <base/point.h>
+#include <base/polynomial.h>
+#include <base/polynomial_space.h>
+#include <base/table.h>
+
+#include <vector>
+
+
+/**
+ * @brief The complete polynomial space of order <tt>k</tt> based on
+ * the monomials.
+ *
+ * This class implements the polynomial space of order <tt>k</tt>
+ * based on the monomials ${1,x,x^2,...}$. I.e. in <tt>d</tt>
+ * dimensions it constructs all polynomials of the form $\prod_{i=1}^d
+ * x_i^{n_i}$, where $\sum_i n_i\leq k$. The base polynomials are
+ * given a specific ordering, e.g. in 2 dimensions:
+ * ${1,x,y,xy,x^2,y^2,x^2y,xy^2,x^3,y^3,...}$. The ordering of the
+ * monomials in $P_k1$ matches the ordering of the monomials in $P_k2$
+ * for $k2>k1$.
+ *
+ * @author Ralf Hartmann, 2004
+ */
+template <int dim>
+class PolynomialsP: public PolynomialSpace<dim>
+{
+  public:
+				     /**
+				      * Constructor. Creates all basis
+				      * functions of $P_k$.
+				      * @arg k: the degree of the
+				      * polynomial space
+				      */
+    PolynomialsP (const unsigned int k);
+
+				     /**
+				      * Returns the degree <tt>k</tt>
+				      * of the polynomial space
+				      * <tt>P_k</tt>.
+				      *
+				      * Note, that this number is
+				      * <tt>PolynomialSpace::degree()-1</tt>,
+				      * compare definition in
+				      * @ref{PolynomialSpace}.
+				      */
+    unsigned int degree() const;
+
+				     /**
+				      * For the <tt>n</tt>th
+				      * polynomial $p_n(x,y,z)=x^i y^j
+				      * z^k$ this function gives the
+				      * degrees i,j,k in the x,y,z
+				      * directions.
+				      */
+    void directional_degrees(unsigned int n,
+			     unsigned int (&degrees)[dim]) const;
+    
+  private:
+
+				     /**
+				      * Fills the <tt>index_map</tt>.
+				      */
+    void create_polynomial_ordering(vector<unsigned int> &index_map) const;
+
+				     /**
+				      * Degree <tt>k<tt> of the
+				      * polynomial space $P_k$,
+				      * i.e. the number <tt>k<tt>
+				      * which was given to the
+				      * constructor.
+				      */
+    const unsigned int k;
+};
+
+
+
+template <int dim>
+inline unsigned int
+PolynomialsP<dim>::degree() const
+{
+  return k;
+}
+
+
+template <int dim>
+inline void
+PolynomialsP<dim>::directional_degrees(unsigned int n,
+				       unsigned int (&degrees)[dim]) const
+{
+  compute_index(n,degrees);
+}
+
+
+#endif
diff --git a/deal.II/base/source/polynomials_p.cc b/deal.II/base/source/polynomials_p.cc
new file mode 100644
index 0000000000..e2ce67889f
--- /dev/null
+++ b/deal.II/base/source/polynomials_p.cc
@@ -0,0 +1,98 @@
+//-------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2004 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//-------------------------------------------------------------------
+
+#include <base/polynomials_p.h>
+
+
+template <int dim>
+PolynomialsP<dim>::PolynomialsP (const unsigned int k)
+		:
+		PolynomialSpace<dim>(Polynomials::Monomial<double>::generate_complete_basis(k)),
+		k(k)
+{
+  vector<unsigned int> index_map(n());
+  create_polynomial_ordering(index_map);
+  set_polynomial_ordering(index_map);
+}
+
+
+template <>
+void PolynomialsP<1>::create_polynomial_ordering(
+  vector<unsigned int> &index_map) const
+{
+  Assert(index_map.size()==n(), ExcDimensionMismatch(index_map.size(), n()));
+
+				   // identity
+  for (unsigned int i=0; i<n(); ++i)
+    index_map[i]=i;
+}
+
+
+
+const unsigned int imap2[6][21]=
+{{0},
+ {0,1,2},
+ {0,1,3,4,2,5},
+ {0,1,4,5,2,7,6,8,3,9},
+ {0,1,5,6,2,9,7,10,3,12,11,8,13,4,14},
+ {0,1,6,7,2,11,8,12,3,15,13,9,16,4,18,14,17,10,19,5,20}
+};
+
+
+template <>
+void PolynomialsP<2>::create_polynomial_ordering(
+  vector<unsigned int> &index_map) const
+{
+  Assert(index_map.size()==n(), ExcDimensionMismatch(index_map.size(), n()));
+  Assert(k<=5, ExcNotImplemented());
+  
+				   // Given the number i of the
+				   // polynomial in
+				   // $1,x,y,xy,x2,y2,...$,
+				   // index_map[i] gives the number of
+				   // the polynomial in
+				   // PolynomialSpace.
+  for (unsigned int i=0; i<n(); ++i)
+    index_map[i]=imap2[k][i];
+}
+
+
+const unsigned int imap3[4][20]=
+{{0},
+ {0,1,2,3},
+ {0,1,3,6,4,7,8,2,5,9},
+ {0,1,4,10,5,11,13,2,7,16,14,6,12,8,15,17,18,3,9,19}
+};
+
+template <>
+void PolynomialsP<3>::create_polynomial_ordering(
+  vector<unsigned int> &index_map) const
+{
+  Assert(index_map.size()==n(), ExcDimensionMismatch(index_map.size(), n()));
+  Assert(k<=3, ExcNotImplemented());
+  
+				   // Given the number i of the
+				   // polynomial in
+				   // $1,x,y,xy,x2,y2,...$,
+				   // index_map[i] gives the number of
+				   // the polynomial in
+				   // PolynomialSpace.
+  for (unsigned int i=0; i<n(); ++i)
+    index_map[i]=imap3[k][i];
+}
+
+
+
+template class PolynomialsP<1>;
+template class PolynomialsP<2>;
+template class PolynomialsP<3>;