00001
00002
00003
00004
00005
00006
00007
00008
00009
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
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
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
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
00750
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