include/deal.II/base/polynomial.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: polynomial.h 25345 2012-03-31 08:37:04Z bangerth @f$
00003 //
00004 //    Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2009, 2010, 2011, 2012 by the deal.II authors
00005 //
00006 //    This file is subject to QPL and may not be  distributed
00007 //    without copyright and license information. Please refer
00008 //    to the file deal.II/doc/license.html for the  text  and
00009 //    further information on this license.
00010 //
00011 //---------------------------------------------------------------------------
00012 #ifndef __deal2__polynomial_h
00013 #define __deal2__polynomial_h
00014 
00015 
00016 
00017 #include <deal.II/base/config.h>
00018 #include <deal.II/base/exceptions.h>
00019 #include <deal.II/base/subscriptor.h>
00020 #include <deal.II/base/point.h>
00021 #include <deal.II/base/std_cxx1x/shared_ptr.h>
00022 
00023 #include <vector>
00024 
00025 DEAL_II_NAMESPACE_OPEN
00026 
00036 namespace Polynomials
00037 {
00038 
00048   template <typename number>
00049   class Polynomial : public Subscriptor
00050   {
00051     public:
00069       Polynomial (const std::vector<number> &coefficients);
00070 
00075       Polynomial (const unsigned int n);
00076 
00087     Polynomial (const std::vector<Point<1> > &lagrange_support_points,
00088                 const unsigned int            evaluation_point);
00089 
00094       Polynomial ();
00095 
00104       number value (const number x) const;
00105 
00122       void value (const number         x,
00123                   std::vector<number> &values) const;
00124 
00134       unsigned int degree () const;
00135 
00148       void scale (const number factor);
00149 
00178       template <typename number2>
00179       void shift (const number2 offset);
00180 
00185       Polynomial<number> derivative () const;
00186 
00193       Polynomial<number> primitive () const;
00194 
00198       Polynomial<number>& operator *= (const double s);
00199 
00203       Polynomial<number>& operator *= (const Polynomial<number>& p);
00204 
00208       Polynomial<number>& operator += (const Polynomial<number>& p);
00209 
00213       Polynomial<number>& operator -= (const Polynomial<number>& p);
00214 
00218       bool operator == (const Polynomial<number> & p)  const;
00219 
00223       void print(std::ostream& out) const;
00224 
00229       template <class Archive>
00230       void serialize (Archive & ar, const unsigned int version);
00231 
00232     protected:
00233 
00238       static void scale (std::vector<number> &coefficients,
00239                          const number         factor);
00240 
00245       template <typename number2>
00246       static void shift (std::vector<number> &coefficients,
00247                          const number2        shift);
00248 
00252       static void multiply (std::vector<number>& coefficients,
00253                             const number factor);
00254 
00262       void transform_into_standard_form ();
00263 
00276       std::vector<number> coefficients;
00277 
00284       bool in_lagrange_product_form;
00285 
00292       std::vector<number> lagrange_support_points;
00293 
00300       number lagrange_weight;
00301   };
00302 
00303 
00310   template <typename number>
00311   class Monomial : public Polynomial<number>
00312   {
00313     public:
00320       Monomial(const unsigned int n,
00321                const double coefficient = 1.);
00322 
00335       static
00336       std::vector<Polynomial<number> >
00337       generate_complete_basis (const unsigned int degree);
00338 
00339     private:
00343       static std::vector<number> make_vector(unsigned int n,
00344                                              const double coefficient);
00345   };
00346 
00347 
00366   class LagrangeEquidistant: public Polynomial<double>
00367   {
00368     public:
00378       LagrangeEquidistant (const unsigned int n,
00379                            const unsigned int support_point);
00380 
00397       static
00398       std::vector<Polynomial<double> >
00399       generate_complete_basis (const unsigned int degree);
00400 
00401     private:
00402 
00411       static
00412       void
00413       compute_coefficients (const unsigned int n,
00414                             const unsigned int support_point,
00415                             std::vector<double>& a);
00416   };
00417 
00418 
00429   std::vector<Polynomial<double> >
00430   generate_complete_Lagrange_basis (const std::vector<Point<1> >& points);
00431 
00432 
00433 
00448   class Legendre : public Polynomial<double>
00449   {
00450     public:
00455       Legendre (const unsigned int p);
00456 
00469       static
00470       std::vector<Polynomial<double> >
00471       generate_complete_basis (const unsigned int degree);
00472 
00473     private:
00477       static std::vector<std_cxx1x::shared_ptr<const std::vector<double> > > shifted_coefficients;
00478 
00492       static std::vector<std_cxx1x::shared_ptr<const std::vector<double> > > recursive_coefficients;
00493 
00506       static void compute_coefficients (const unsigned int p);
00507 
00515       static const std::vector<double> &
00516       get_coefficients (const unsigned int k);
00517   };
00518 
00540   class Lobatto : public Polynomial<double>
00541   {
00542     public:
00549       Lobatto (const unsigned int p = 0);
00550 
00558       static std::vector<Polynomial<double> >
00559       generate_complete_basis (const unsigned int p);
00560 
00561     private:
00565       std::vector<double> compute_coefficients (const unsigned int p);
00566   };
00567 
00568 
00569 
00610   class Hierarchical : public Polynomial<double>
00611   {
00612     public:
00619       Hierarchical (const unsigned int p);
00620 
00640       static
00641       std::vector<Polynomial<double> >
00642       generate_complete_basis (const unsigned int degree);
00643 
00644     private:
00648       static void compute_coefficients (const unsigned int p);
00649 
00657      static const std::vector<double> &
00658      get_coefficients (const unsigned int p);
00659 
00673     static std::vector<std_cxx1x::shared_ptr<const std::vector<double> > > recursive_coefficients;
00674    };
00675 }
00676 
00677 
00678 
00681 /* -------------------------- inline functions --------------------- */
00682 
00683 namespace Polynomials
00684 {
00685   template <typename number>
00686   inline
00687   Polynomial<number>::Polynomial ()
00688     :
00689     in_lagrange_product_form (false),
00690     lagrange_weight          (1.)
00691   {}
00692 
00693 
00694 
00695   template <typename number>
00696   inline
00697   unsigned int
00698   Polynomial<number>::degree () const
00699   {
00700     if (in_lagrange_product_form == true)
00701       {
00702         return lagrange_support_points.size();
00703       }
00704     else
00705       {
00706         Assert (coefficients.size()>0, ExcEmptyObject());
00707         return coefficients.size() - 1;
00708       }
00709   }
00710 
00711 
00712 
00713   template <typename number>
00714   inline
00715   number
00716   Polynomial<number>::value (const number x) const
00717   {
00718     if (in_lagrange_product_form == false)
00719       {
00720         Assert (coefficients.size() > 0, ExcEmptyObject());
00721 
00722                                      // Horner scheme
00723         const unsigned int m=coefficients.size();
00724         number value = coefficients.back();
00725         for (int k=m-2; k>=0; --k)
00726           value = value*x + coefficients[k];
00727         return value;
00728       }
00729     else
00730       {
00731                                 // direct evaluation of Lagrange polynomial
00732         const unsigned int m = lagrange_support_points.size();
00733         number value = 1.;
00734         for (unsigned int j=0; j<m; ++j)
00735           value *= x-lagrange_support_points[j];
00736         value *= lagrange_weight;
00737         return value;
00738       }
00739   }
00740 
00741 
00742 
00743   template <typename number>
00744   template <class Archive>
00745   inline
00746   void
00747   Polynomial<number>::serialize (Archive & ar, const unsigned int)
00748   {
00749                                      // forward to serialization
00750                                      // function in the base class.
00751     ar & static_cast<Subscriptor &>(*this);
00752     ar & coefficients;
00753     ar & in_lagrange_product_form;
00754     ar & lagrange_support_points;
00755     ar & lagrange_weight;
00756   }
00757 
00758 }
00759 DEAL_II_NAMESPACE_CLOSE
00760 
00761 #endif
00762 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

deal.II documentation generated on Wed May 23 2012 06:07:31 by doxygen 1.7.3