NGSolve
4.9
|
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