NGSolve  4.9
fem/generic_recpol.hpp
00001 #ifndef FILE_GENERIC_RECPOL
00002 #define FILE_GENERIC_RECPOL
00003 
00004 
00005 /*********************************************************************/
00006 /* File:   generic_recpol.hpp                                        */
00007 /* Author: Veronika Pillwein and Joachim Schöberl                    */
00008 /* Date:   5. Jan. 2005                                              */
00009 /*********************************************************************/
00010 
00011 namespace ngfem
00012 {
00013 
00014 
00025   class RecPol
00026   {
00027   protected:
00029     Array< double[3] > coefs;
00030   
00031   public:
00033     RecPol (int n) : coefs(n+1) 
00034     { 
00035       for (int i = 0; i <= n; i++)
00036         coefs[i][0] = coefs[i][1] = coefs[i][2] = 0.0;
00037     }
00038 
00040     int Order() const { return coefs.Size()-1; }
00041   
00043     template <class S, class T>
00044     void Evaluate (int n, S x, T & values) const
00045     {
00046       // T & values = const_cast<T&> (const_values);
00047       S v1, v2, v3;
00048 
00049       if (n < 0) return;
00050       v2 = values[0] = coefs[0][0];
00051 
00052       if (n == 0) return;
00053       v1 = values[1] = (coefs[1][0] + coefs[1][1] * x) * values[0];
00054 
00055       /*
00056       // hat nichts gebracht
00057       int i;
00058       for (i = 2; i < n; i+=2)
00059       {
00060       v2 *= coefs[i][2];
00061       v2 += (coefs[i][0] + coefs[i][1] * x) * v1;
00062       values[i] = v2;
00063 
00064       v1 *= coefs[i+1][2];
00065       v1 += (coefs[i+1][0] + coefs[i+1][1] * x) * v2;
00066       values[i+1] = v1;
00067       }
00068 
00069       if (i <= n)
00070       {
00071       v2 *= coefs[i][2];
00072       v2 += (coefs[i][0] + coefs[i][1] * x) * v1;
00073       values[i] = v2;
00074       }
00075       */
00076     
00077       for (int i = 2; i <= n; i++)
00078         {
00079           v3 = v2; v2 = v1;
00080           v1 = (coefs[i][0] + coefs[i][1] * x) * v2 + coefs[i][2] * v3;
00081           values[i] = v1;
00082         }
00083     }
00084 
00085 
00086     /*
00087       void Evaluate (int n, double x, double * values) const
00088       {
00089       double v1, v2, v3;
00090 
00091       if (n >= 0)
00092       v2 = values[0] = coefs[0][0];
00093       if (n >= 1)
00094       v1 = values[1] = (coefs[1][0] + coefs[1][1] * x) * values[0];
00095     
00096       for (int i = 2; i <= n; i++)
00097       {
00098       v3 = v2; v2 = v1;
00099       v1 = (coefs[i][0] + coefs[i][1] * x) * v2 + coefs[i][2] * v3;
00100       values[i] = v1;
00101       }
00102       }
00103     */
00104 
00111     void MultBubble ();
00112 
00119     void MultLinear (double a, double b);
00120 
00121   
00122     void Scale (double s)
00123     { coefs[0][0] *= s; }
00124 
00125 
00126     double A(int i) const { return coefs[i][0]; }
00127     double B(int i) const { return coefs[i][1]; }
00128     double C(int i) const { return coefs[i][2]; }
00129   };
00130 
00131 
00133   inline ostream & operator<< (ostream & ost, const RecPol & pol)
00134   {
00135     for (int i = 0; i <= pol.Order(); i++)
00136       ost << "a(" << i << ") = " << pol.A(i)
00137           << ", b(" << i << ") = " << pol.B(i)
00138           << ", c(" << i << ") = " << pol.C(i) << endl;
00139     return ost;
00140   }
00141 
00142 
00146   class RecPolMonomial : public RecPol
00147   {
00148   public:
00150     RecPolMonomial (int n)
00151       : RecPol (n)
00152     {
00153       coefs[0][0] = 1;
00154       for (int i = 1; i <= n; i++)
00155         coefs[i][1] = 1;
00156     }
00157   };
00158 
00159 
00163   class RecPolJacobi : public RecPol
00164   {
00165   public:
00166     RecPolJacobi (int n, double alpha, double beta)
00167       : RecPol (n)
00168     {
00169       coefs[0][0] = 1;
00170       coefs[1][0] = alpha+1 - 0.5 * (alpha + beta + 2);
00171       coefs[1][1] = 0.5 * (alpha+beta+2);
00172 
00173       for (int i = 2; i <= n; i++)
00174         {
00175           double fac = 1.0 / (2.0*i* (i+alpha+beta)*(2*i-2+alpha+beta));        
00176           coefs[i][0] = ( 2*i-1+alpha+beta ) * ( alpha*alpha-beta*beta ) * fac;
00177           coefs[i][1] = ( 2*i-2+alpha+beta ) * ( 2*i-1+alpha+beta ) * ( 2*i+alpha+beta ) * fac;
00178           coefs[i][2] = -2 * (i-1+alpha) * (i-1+beta) * (2*i+alpha+beta) * fac;
00179         }
00180     }
00181   };
00182 
00183 
00187   class RecPolLegendre : public RecPol
00188   {
00189   public:
00191     RecPolLegendre (int n)
00192       : RecPol (n)
00193     {
00194       coefs[0][0] = 1;
00195       coefs[1][1] = 1;
00196       for (int i = 2; i <= n; i++)
00197         {
00198           coefs[i][1] = (2.0*i-1)/i;
00199           coefs[i][2] = (1.0-i)/i;
00200         }
00201     }
00202   };
00203 
00204 
00205 
00206 
00207 
00208 
00209 
00210 
00211 
00212 
00214   extern void GenerateMatrix (const RecPol & pol1, const RecPol & pol2,
00215                               FlatVector<> pts, FlatVector<> coefs,
00216                               FlatMatrix<> mat);
00217 
00218 
00219 
00220 
00222   extern void GenerateMatrix (const RecPol & pol1, FlatMatrix<> pol1_vals,
00223                               const RecPol & pol2, FlatMatrix<> pol2_vals,
00224                               FlatVector<> coefs, FlatMatrix<> mat);
00225 
00226 
00227 
00228 }
00229 
00230 #endif
00231