NGSolve  4.9
fem/recursive_pol_trig.hpp
00001 namespace ngfem
00002 {
00003 
00008   class TrigShapesInnerLegendre
00009   {
00010   public:
00012     template <typename Sx, typename Sy, typename T>
00013     static void CalcSplitted (int n, Sx x, Sy y, T & val1, T & val2)
00014     {
00015       LegendrePolynomial leg;
00016 
00017       Sx bub1 = (1-x-y)*(1+x-y);
00018       // ScaledLegendrePolynomialMult (n-3, x, 1-y, bub1, val1); 
00019       leg.EvalScaledMult (n-3, x, 1-y, bub1, val1); 
00020 
00021       //LegendrePolynomialMult (n-3, 2*y-1, y, val2);
00022       leg.EvalMult (n-3, 2*y-1, y, val2);
00023 
00024       /*
00025         ScaledLegendrePolynomial (n-3, x, 1-y, val1); 
00026         for (int i = 0; i <= n-3; i++)
00027         val1[i] *= bub1;
00028 
00029         LegendrePolynomial (n-3, 2*y-1, val2);
00030         for (int i = 0; i <= n-3; i++)
00031         val2[i] *= y;
00032       */
00033     }
00034 
00035 
00036     template <int n, typename Sx, typename Sy, typename T>
00037     static void CalcSplitted (Sx x, Sy y, T & val1, T & val2)
00038     {
00039       Sx bub1 = (1-x-y)*(1+x-y);
00040       // ScaledLegendrePolynomialMult (n-3, x, 1-y, bub1, val1); 
00041       // LegendrePolynomialFO<n-3>::EvalScaledMult (x, 1-y, bub1, val1); 
00042       LegendrePolynomial leg;
00043       leg.EvalScaledMultFO<n-3> (x, 1-y, bub1, val1);
00044 
00045       // LegendrePolynomialMult (n-3, 2*y-1, y, val2);
00046       // LegendrePolynomialFO<n-3>::EvalMult (2*y-1, y, val2);
00047       leg.EvalMultFO<n-3> (2*y-1, y, val2);
00048     }
00049 
00050 
00052     template < typename Sx, typename Sy, typename T>
00053     static int Calc (int n, Sx x, Sy y, T & values)
00054     {
00055       ArrayMem<Sx, 20> polx(n-2), poly(n-2);
00056 
00057       ScaledLegendrePolynomial (n-3, x, 1-y, polx);
00058       LegendrePolynomial (n-3, 2*y-1, poly);
00059       Sx bub = y * (1-x-y) * (1+x-y);
00060 
00061       int ii = 0;
00062       for (int i = 0; i <= n-3; i++)
00063         for (int j = 0; j <= n-3-i; j++)
00064           values[ii++] = bub * polx[i] * poly[j];
00065 
00066       return ii;
00067     }
00068 
00069     template <int n, typename Sx, typename Sy, typename T>
00070     static int Calc (Sx x, Sy y, T & values)
00071     {
00072       // ArrayMem<Sx, 20> polx(n-2), poly(n-2);
00073       Sx polx[n], poly[n];
00074 
00075       /*
00076         ScaledLegendrePolynomial (n-3, x, 1-y, polx);
00077         LegendrePolynomial (n-3, 2*y-1, poly);
00078         Sx bub = y * (1-x-y) * (1+x-y);
00079       */
00080       CalcSplitted<n> (x, y, polx, poly);
00081       int ii = 0;
00082       for (int i = 0; i <= n-3; i++)
00083         for (int j = 0; j <= n-3-i; j++)
00084           values[ii++] = /* bub * */ polx[i] * poly[j];
00085 
00086       return ii;
00087     }
00088 
00089 
00090 
00091 
00092 
00094     template <typename Sx, typename Sy, typename Sf, typename T>
00095     static int CalcMult (int n, Sx x, Sy y, Sf & fac, T & values)
00096     {
00097       ArrayMem<Sx, 20> polx(n-2), poly(n-2);
00098 
00099       ScaledLegendrePolynomial (n-3, x, 1-y, polx);
00100       LegendrePolynomial (n-3, 2*y-1, poly);
00101 
00102       Sx bub = fac * y * (1-x-y) * (1+x-y);
00103 
00104       int ii = 0;
00105       for (int i = 0; i <= n-3; i++)
00106         for (int j = 0; j <= n-3-i; j++)
00107           values[ii++] = bub * polx[i] * poly[j];
00108 
00109       return ii;
00110     }
00111 
00112 
00113   };
00114 
00115 
00120   class TrigShapesInnerJacobi
00121   {
00122   public:
00124     template <typename Sx, typename Sy, typename T>
00125     static int Calc (int n, Sx x, Sy y, T & values)
00126     {
00127       int ii = 0;
00128       ArrayMem<Sx, 20> polx(n+1), poly(n+1);
00129 
00130       Sx bub = y * (1-x-y) * (1+x-y);
00131       ScaledJacobiPolynomial (n-3, x, 1-y, 2, 2, polx);
00132 
00133       for (int ix = 0; ix <= n-3; ix++)
00134         {
00135           JacobiPolynomial (n-3, 2*y-1, 2*ix+5, 2, poly);
00136           for (int j = 0; j <= n-3-ix; j++)
00137             values[ii++] = bub * polx[ix] * poly[j];
00138         }
00139       return ii;
00140     }
00141 
00142 
00143   };
00144 
00145 
00146   //*************************************MONOMIAL EXTENSION***********************************************************
00147 
00161   class TrigExtensionMonomial
00162   {
00163   public:
00165     template <class Sx, class Sy, class T>
00166     inline static int CalcTrigExt (int n, Sx x, Sy y, T & values)
00167     {
00168       Sy fy = (1-y)*(1-y);
00169       Sx p3 = 0;
00170       Sx p2 = -1;
00171       Sx p1 = x;
00172 
00173       for (int j=2; j<=n; j++)
00174         {
00175           p3=p2; p2=p1;
00176           p1=( (2*j-3) * x * p2 - (j-3) * fy * p3) / j;
00177           values[j-2] = p1;
00178         }
00179       return n-1;
00180     }
00181 
00182 
00184     template <class T>
00185     inline static int CalcTrigExtDeriv (int n, double x, double y, T & values)
00186     {
00187       double fy = (1-y)*(1-y);
00188       double p3 = 0, p3x = 0, p3y = 0;
00189       double p2 = -1, p2x = 0, p2y = 0;
00190       double p1 = x, p1x = 1, p1y = 0;
00191 
00192       for (int j=2; j<=n; j++)
00193         {
00194           p3=p2; p3x = p2x; p3y = p2y;
00195           p2=p1; p2x = p1x; p2y = p1y;
00196           double c1 = (2.0*j-3) / j;
00197           double c2 = (j-3.0) / j;
00198 
00199           p1  = c1 * x * p2 - c2 * fy * p3;
00200           p1x = c1 * p2 + c1 * x * p2x - c2 * fy * p3x;
00201           p1y = c1 * x * p2y - (c2 * 2 * (y-1) * p3 + c2 * fy * p3y);
00202           values (j-2, 0) = p1x;
00203           values (j-2, 1) = p1y;
00204         }
00205       return n-1;
00206     }
00207 
00209     template <class Sx, class T>
00210     inline static int Calc (int n, Sx x, T & values)
00211     {
00212       Sx p3 = 0;
00213       Sx p2 = -1;
00214       Sx p1 = x;
00215 
00216       for (int j=2; j<=n; j++)
00217         {
00218           p3=p2; p2=p1;
00219           p1=( (2*j-3) * x * p2 - (j-3) * p3) / j;
00220           values[j-2] = p1;
00221         }
00222       return n-1;
00223 
00224     }
00225 
00227     template <class T>
00228     inline static int CalcDeriv (int n, double x, T & values)
00229     {
00230       double p1 = 1.0, p2 = 0.0, p3;
00231 
00232       for (int j=1; j<=n-1; j++)
00233         {
00234           p3 = p2; p2 = p1;
00235           p1 = ((2.0*j-1.0)*x*p2 - (j-1.0)*p3) / j;
00236           values[j-1] = p1;
00237         }
00238       return n-1;
00239   
00240     }
00241   };
00242 
00243 
00244   // ***********************  OPTIMAL EXTENSION  *****************************************
00258   class TrigExtensionOptimal
00259   {
00260     enum { SIZE = 1000 };
00261     static double coefs[SIZE][5];
00262     static bool initialized;
00263   public:
00264 
00265     TrigExtensionOptimal ();
00266 
00267     template <class Sx, class Sy, class T>
00268     inline static int CalcTrigExt (int p, Sx x, Sy y, T & values)
00269     {
00270       Sx p1(0.0), p2(0.0), p3(0.0), p4(0.0), p5(0.0);
00271       Sx bub = (1.+x-y)*(x+y-1);
00272       switch(p)
00273         {
00274         default:
00275         case 5:
00276 
00277           {
00278             p5 = -1./24. * bub*x * (-9.+21.*x*x-14.*y+35.*y*y);
00279             values[3] = p5;
00280           }
00281 
00282         case 4:
00283 
00284           {
00285             p4 = -0.125 * bub * (-1.+5.*x*x-2.*y+3.*y*y);
00286             values[2] = p4;
00287           }
00288 
00289         case 3:
00290 
00291           {
00292             p3 = -0.5 * bub * x;
00293             values[1] = p3;
00294           }
00295 
00296         case 2:
00297 
00298           {
00299             p2 = -0.5 * bub;
00300             values[0] = p2;
00301           }
00302         case 1:
00303         case 0:
00304           ;
00305         }
00306 
00307       for (int j = 6; j <= p; j++)
00308         {
00309           p1 = 
00310             -coefs[j-6][0] * p2 
00311             + coefs[j-6][1] * x * p3 
00312             - (coefs[j-6][2]+coefs[j-6][3]*(x*x-y*y)) * p4 
00313             + coefs[j-6][4] * x * p5;
00314           p2 = p3; p3 = p4; p4 = p5; p5 = p1;
00315           values[j-2] = p1;
00316         }
00317 
00318       return p-1;
00319     }
00320 
00321     template <class T>
00322     inline static int CalcTrigExtDeriv (int n, double x, double y, T & values)
00323     {
00324       ArrayMem<AutoDiff<2>,10> ad_values(n-1);
00325       AutoDiff<2> ad_x(x, 0);
00326       AutoDiff<2> ad_y(y, 1);
00327 
00328       CalcTrigExt (n, ad_x, ad_y, ad_values);
00329 
00330       for (int i = 0; i < n-1; i++)
00331         for (int j = 0; j < 2; j++)
00332           values(i,j) = ad_values[i].DValue(j);
00333       return n-1;
00334     }
00335 
00336 
00337     template <class Sx, class T>
00338     inline static int Calc (int n, Sx x, T & values)
00339     {
00340       Sx y = 0.0;
00341       return CalcTrigExt (n, x, y, values);
00342     }
00343 
00344     template <class T>
00345     inline static int CalcDeriv (int n, double x, T & values)
00346     {
00347       return CalcTrigExtDeriv (n, x, 0.0, values);
00348     }
00349   };
00350 
00351 
00352 
00353   // ************************ MINIMAL EXTENSION *************************
00354 
00355 
00359   class TrigExtensionMin
00360   {
00361   public:
00362     template <class Sx, class Sy, class T>
00363     inline static int CalcTrigExt (int n, Sx x, Sy y, T & values)
00364     {
00365 #ifdef ABC
00366       //    TrigExtensionOptimal::CalcTrigExt (n, x, y, values);
00367 
00368       ArrayMem<Sx, 20> polx(n+1);
00369 
00370       const IntegrationRule & rule =
00371         GetIntegrationRules().SelectIntegrationRule (ET_SEGM, n+2);
00372 
00373       Sx bub = (1-x-y) * (1+x-y);
00374       for (int i = 0; i <= n-2; i++)
00375         values[i] = 0;
00376 
00377       for (int i = 0; i < rule.GetNIP(); i++)
00378         {
00379           const IntegrationPoint & ip = rule.GetIP(i);
00380 
00381           Sx hx = x + y * (2.0*ip(0)-1.0);
00382                          
00383           JacobiPolynomial (n-2, hx, 1, 1, polx);
00384           // LegendrePolynomial (n-2, hx, polx);
00385 
00386           Sx fac = bub * ip.Weight(); //  / (z*z);
00387 
00388           for (int j = 0; j <= n-2; j++)
00389             values[j] += fac * polx[j];
00390         }
00391 
00392 
00393 
00394       // IntegratedLegendreMonomialExt::CalcTrigExt (n, x, y, values);
00395 
00396       /*
00397         Sy fy = (1-y)*(1-y);
00398         Sx p3 = 0;
00399         Sx p2 = -1;
00400         Sx p1 = x;
00401 
00402         for (int j=2; j<=n; j++)
00403         {
00404         p3=p2; p2=p1;
00405         p1=( (2*j-3) * x * p2 - (j-3) * fy * p3) / j;
00406         values[j-2] = p1;
00407         }
00408       */
00409 
00410       /*
00411         Sx fac = 1;
00412         for (int j = n-2; j >= 0; j--)
00413         {
00414         values[j] = fac * values[j];
00415         fac = fac * (1-y);
00416         }
00417       */
00418 
00419       for (int j = 0; j < n-3; j++)
00420         values[j] *= VertexExtensionOptimal<2>::Calc (n-2-j, 1-y);
00421 
00422       return n-1;
00423 #endif
00424 
00425       if (y == 0.0) y += 1e-8;
00426       ArrayMem<Sx, 20> pol1(n+2), pol2(n+2);
00427       JacobiPolynomial (n+1, x-y, 1, 1, pol1);
00428       JacobiPolynomial (n+1, x+y, 1, 1, pol2);
00429       for (int i = 0; i <= n-2; i++)
00430         values[i] = (pol2[i+1]-pol1[i+1]) / (2*y);
00431 
00432       /*
00433         ScaledJacobiPolynomial (n-2, x, 1-y, 2, 2, values);
00434         for (int i = 0; i < n-1; i++)
00435         values[i] *= 0.5 * (i+4) * (1-x*x);
00436       */
00437 
00438 
00439 
00440       LowEnergyVertexPolynomials2D (n, 1-2*y, pol1);
00441       for (int i = 0; i <= n-2; i++)
00442         values[i] *= pol1[n-2-i] * (1-x-y) * (1+x-y);
00443 
00444       return n-1;
00445     }
00446 
00447     template <class T>
00448     inline static int CalcTrigExtDeriv (int n, double x, double y, T & values)
00449     {
00450       ArrayMem<AutoDiff<2>,10> ad_values(n-1);
00451       AutoDiff<2> ad_x(x, 0);
00452       AutoDiff<2> ad_y(y, 1);
00453 
00454       CalcTrigExt (n, ad_x, ad_y, ad_values);
00455 
00456       for (int i = 0; i < n-1; i++)
00457         for (int j = 0; j < 2; j++)
00458           values(i,j) = ad_values[i].DValue(j);
00459       return n-1;
00460     }
00461 
00462 
00463 
00464     template <class Sx, class T>
00465     inline static int Calc (int n, Sx x, T & values)
00466     {
00467       /*
00468         Sx y = 0.0;
00469         CalcTrigExt (n, x, y, values);
00470       */
00471 
00472       JacobiPolynomial (n-2, x, 2, 2, values);
00473       for (int i = 0; i < n-1; i++)
00474         values[i] *= 0.5 * (i+4) * (1-x*x);
00475 
00476       return n-1;
00477 
00478     }
00479 
00480     template <class T>
00481     inline static int CalcDeriv (int n, double x, T & values)
00482     {
00483       ArrayMem<AutoDiff<1>,10> ad_values(n-1);
00484       AutoDiff<1> ad_x(x, 0);
00485 
00486       Calc (n, ad_x, ad_values);
00487 
00488       for (int i = 0; i < n-1; i++)
00489         values[i] = ad_values[i].DValue(0);
00490       return n-1;
00491 
00492       // return CalcTrigExtDeriv (n, x, 0.0, values);
00493     }
00494   };
00495 
00496 
00497 }