NGSolve  4.9
fem/recursive_pol_tet.hpp
00001 namespace ngfem
00002 {
00003 
00004   class TetShapesInnerLegendre
00005   {
00006   public:
00007     template <typename Sx, typename Sy, typename Sz, typename T>
00008     static int Calc (int n, Sx x, Sy y, Sz z, T & values)
00009     {
00010       int ii = 0, i, j, k;
00011       ArrayMem<Sx, 20> polx(n+1), poly(n+1), polz(n+1);
00012     
00013       ScaledLegendrePolynomial (n-4, x, 1-y-z, polx);
00014       ScaledLegendrePolynomial (n-4, 2*y - (1-z) , (1-z),poly);  //SZ (1-z)-1, poly);
00015       LegendrePolynomial (n-4, 2*z-1, polz);
00016       Sx bub = (1-x-y-z) * (1+x-y-z) * y * z ; 
00017   
00018       for (i = 0; i <= n-4; i++)
00019         for (j = 0; j <= n-4-i; j++)
00020           for (k = 0; k <= n-4-i-j; k++)
00021             values[ii++] = bub * polx[i] * poly[j] * polz[k];
00022     
00023       return ii;
00024     }
00025 
00026     template <typename Sx, typename Sy, typename Sz, typename T1, typename T2, typename T3>
00027     static void CalcSplitted (int n, Sx x, Sy y, Sz z, T1 & val1, T2 & val2, T3 & val3)
00028     {
00029       Sx bub1 = sqr (1-y-z) - sqr(x);
00030       ScaledLegendrePolynomialMult (n-4, x, 1-y-z, bub1, val1);
00031       ScaledLegendrePolynomialMult (n-4, 2*y - (1-z) , (1-z), y, val2); 
00032 
00033       // LegendrePolynomialMult (n-4, 2*z-1, z, val3);
00034       LegendrePolynomial leg;
00035       leg.EvalMult (n-4, 2*z-1, z, val3);
00036 
00037       /*
00038         ScaledLegendrePolynomial (n-4, x, 1-y-z, val1);
00039         // Sx bub1 = (1-x-y-z) * (1+x-y-z);
00040         Sx bub1 = sqr (1-y-z) - sqr(x);
00041         for (int i = 0; i <= n-4; i++)
00042         val1[i] *= bub1;
00043 
00044         ScaledLegendrePolynomial (n-4, 2*y - (1-z) , (1-z), val2); 
00045         for (int i = 0; i <= n-4; i++)
00046         val2[i] *= y;
00047 
00048         LegendrePolynomial (n-4, 2*z-1, val3);
00049         for (int i = 0; i <= n-4; i++)
00050         val3[i] *= z;
00051       */
00052     }
00053   };
00054 
00055 
00056 
00057   class TetShapesInnerJacobi
00058   {
00059   public:
00060     template <typename Sx, typename Sy, typename Sz, typename T>
00061     static int Calc (int n, Sx x, Sy y, Sz z, T & values)
00062     {
00063       int ii = 0, j, k;
00064       ArrayMem<Sx, 20> polx(n+1), poly(n+1), polz(n+1);
00065 
00066       Sx bub = y * z * (1-x-y-z) * (1+x-y-z);
00067       ScaledJacobiPolynomial (n-4, x, (1-y-z), 2, 2, polx);
00068     
00069       for (int ix = 0; ix <= n-4; ix++)
00070         {
00071           ScaledJacobiPolynomial (n-4, (2*y-1+z),(1-z), 2*ix+5, 2, poly);
00072           JacobiPolynomial (n-4, 2*z-1, 2*ix+5, 2, polz);
00073         
00074           for (j = 0; j <= n-4-ix; j++)
00075             for (k = 0; k <= n-4-ix-j; k++)
00076               values[ii++] = bub * polx[ix] * poly[j] * polz[k];
00077         }
00078       return ii;
00079     }
00080   };
00081 
00082   class TetShapesFaceLegendre
00083   {
00084   public:
00085     template <typename Sx, typename Sy, typename Sz, typename T>
00086     static int Calc (int n, Sx x, Sy y, Sz z, T & values)
00087     {
00088       ArrayMem<Sx, 20> polx(n+1), poly(n+1);
00089     
00090       ScaledLegendrePolynomial (n-3, x, 1-y-z, polx);
00091       ScaledLegendrePolynomial (n-3, 2*y-(1-z),(1-z), poly);
00092       Sx bub = (1-x-y-z) * (1+x-y-z)*y; 
00093 
00094       int ii = 0;
00095       for (int i = 0; i <= n-3; i++)
00096         for (int j = 0; j <= n-3-i; j++)
00097           values[ii++] = bub * polx[i] * poly[j];
00098  
00099       return ii;
00100     }
00101     template <typename Sx, typename Sy, typename Sz, typename T>
00102     static void CalcSplitted (int n, Sx x, Sy y, Sz z, T & val1, T & val2)
00103     {
00104       Sx bub1 = sqr (1-y-z) - sqr (x); 
00105       ScaledLegendrePolynomialMult (n-3, x, 1-y-z, bub1, val1);
00106       ScaledLegendrePolynomialMult (n-3, 2*y-(1-z),(1-z), y, val2);
00107       /*
00108         ScaledLegendrePolynomial (n-3, x, 1-y-z, val1);
00109 
00110         // Sx bub1 = (1-x-y-z) * (1+x-y-z);
00111         Sx bub1 = sqr (1-y-z) - sqr (x); 
00112         for (int i = 0; i <= n-3; i++)
00113         val1[i] *= bub1;
00114 
00115         ScaledLegendrePolynomial (n-3, 2*y-(1-z),(1-z), val2);
00116         for (int i = 0; i <= n-3; i++)
00117         val2[i] *= y;
00118       */
00119     }
00120   };
00121 
00122 
00123 
00124 
00125   class TetShapesFaceJacobi
00126   {
00127   public:
00128     template <typename Sx, typename Sy, typename Sz, typename T>
00129     static int Calc (int n, Sx x, Sy y, Sz z, T & values)
00130     {
00131       int ii = 0;
00132       ArrayMem<Sx, 20> polx(n+1), poly(n+1);
00133 
00134       Sx bub = y * (1-x-y-z) * (1+x-y-z);
00135       ScaledJacobiPolynomial (n-3, x, 1-y-z, 2, 2, polx);
00136     
00137       for (int ix = 0; ix <= n-3; ix++)
00138         {
00139           ScaledJacobiPolynomial (n-3, 2*y-1+z, 1-z, 2*ix+5, 2, poly);
00140           for (int j = 0; j <= n-3-ix; j++)
00141             values[ii++] = bub * polx[ix] * poly[j];
00142         }
00143       return ii;
00144     }
00145   };
00146 
00147 
00148   class TetShapesFaceOpt1
00149   {
00150   public:
00151     template <typename Sx, typename Sy, typename Sz, typename T>
00152     static int Calc (int n, Sx x, Sy y, Sz z, T & values)
00153     {
00154       int nd = Calc1 (n, x, y, z, values);
00155       ArrayMem<Sx, 100> hvalues(nd);
00156 
00157       Sx lam1 = 0.5 * (1 + x - y - z);
00158       Sx lam2 = 0.5 * (1 - x - y - z);
00159       Sx lam3 = y;
00160       Sx lam4 = z;
00161 
00162       Sx hlam1 = 0;
00163       Sx hlam2 = lam2;
00164       Sx hlam3 = lam3;
00165       Sx hlam4 = lam4 + lam1;
00166 
00167       Sx hx = hlam1 - hlam2;
00168       Sx hy = hlam3;
00169       Sx hz = hlam4;
00170 
00171       Sx frac;
00172       if (hlam4 < 1e-12)
00173         frac = 0.0;
00174       else
00175         frac = lam4 / hlam4;
00176 
00177       Calc1 (n, hx, hy, hz, hvalues);
00178       for (int i = 0; i < nd; i++)
00179         values[i] -= frac * hvalues[i];
00180 
00181       return nd;
00182     }
00183   
00184     template <typename Sx, typename Sy, typename Sz, typename T>
00185     static int Calc1 (int n, Sx x, Sy y, Sz z, T & values)
00186     {
00187       int nd = Calc2 (n, x, y, z, values);
00188       ArrayMem<Sx, 100> hvalues(nd);
00189 
00190       Sx lam1 = 0.5 * (1 + x - y - z);
00191       Sx lam2 = 0.5 * (1 - x - y - z);
00192       Sx lam3 = y;
00193       Sx lam4 = z;
00194 
00195       Sx hlam1 = lam1;
00196       Sx hlam2 = 0;
00197       Sx hlam3 = lam3;
00198       Sx hlam4 = lam4 + lam2;
00199 
00200       Sx hx = hlam1 - hlam2;
00201       Sx hy = hlam3;
00202       Sx hz = hlam4;
00203 
00204       Sx frac;
00205       if (hlam4 < 1e-12)
00206         frac = 0.0;
00207       else
00208         frac = lam4 / hlam4;
00209 
00210       Calc2 (n, hx, hy, hz, hvalues);
00211       for (int i = 0; i < nd; i++)
00212         values[i] -= frac * hvalues[i];
00213     
00214       return nd;
00215     }
00216 
00217     template <typename Sx, typename Sy, typename Sz, typename T>
00218     static int Calc2 (int n, Sx x, Sy y, Sz z, T & values)
00219     {
00220       Sx * hp = &values[0];
00221       int nd = Calc3 (n, x, y, z, hp);
00222       ArrayMem<Sx, 100> hvalues(nd);
00223 
00224       Sx lam1 = 0.5 * (1 + x - y - z);
00225       Sx lam2 = 0.5 * (1 - x - y - z);
00226       Sx lam3 = y;
00227       Sx lam4 = z;
00228 
00229       Sx hlam1 = lam1;
00230       Sx hlam2 = lam2;
00231       Sx hlam3 = 0;
00232       Sx hlam4 = lam4 + lam3;
00233 
00234       Sx hx = hlam1 - hlam2;
00235       Sx hy = hlam3;
00236       Sx hz = hlam4;
00237 
00238       Sx frac;
00239       if (hlam4 < 1e-12)
00240         frac = 0.0;
00241       else
00242         frac = lam4 / hlam4;
00243 
00244       hp = &hvalues[0];
00245       Calc3 (n, hx, hy, hz, hp);
00246       for (int i = 0; i < nd; i++)
00247         values[i] -= frac * hvalues[i];
00248     
00249       return nd;
00250     }
00251 
00252 
00253 
00254 
00255 
00256 
00257     template <typename Sx, typename Sy, typename Sz, typename T>
00258     static int Calc3 (int n, Sx x, Sy y, Sz z, T & values)
00259     {
00260       int ii = 0;
00261       ArrayMem<Sx, 20> polx(n+1), poly(n+1);
00262 
00263       const IntegrationRule & rule = SelectIntegrationRule (ET_TRIG, n+2);
00264 
00265       for (int ix = 0; ix <= n-3; ix++)
00266         for (int j = 0; j <= n-3-ix; j++)
00267           values[ii++] = 0;
00268       for (int i = 0; i < rule.GetNIP(); i++)
00269         {
00270           ii = 0;
00271           const IntegrationPoint & ip = rule[i];
00272 
00273           Sx hx = x + z * (-1+2*ip(0)+ip(1));
00274           Sy hy = y + z * ip(1);
00275 
00276           Sx bub = hy * (1-hx-hy) * (1+hx-hy);
00277           ScaledJacobiPolynomial (n-3, hx, 1-hy, 2, 2, polx);
00278 
00279           Sx fac = 2 * bub * ip.Weight(); //  / (z*z);
00280 
00281           for (int ix = 0; ix <= n-3; ix++)
00282             {
00283               ScaledJacobiPolynomial (n-3, 2*hy-1, 1, 2*ix+5, 2, poly);
00284               for (int j = 0; j <= n-3-ix; j++)
00285                 values[ii++] += fac * polx[ix] * poly[j];
00286             }
00287         }
00288       return ii;
00289     }
00290   };
00291 
00292 
00293 
00294 
00295   class TetShapesFaceOpt2
00296   {
00297   public:
00298     template <typename Sx, typename Sy, typename Sz, typename T>
00299     static int Calc (int n, Sx x, Sy y, Sz z, T & values)
00300     {
00301       int ii = 0, i, j;
00302       ArrayMem<Sx, 20> polx(n+1), poly(n+1);
00303 
00304       const IntegrationRule & rule = SelectIntegrationRule (ET_TRIG, n+2);
00305 
00306       for (int ix = 0; ix <= n-3; ix++)
00307         for (j = 0; j <= n-3-ix; j++)
00308           values[ii++] = 0;
00309       for (i = 0; i < rule.GetNIP(); i++)
00310         {
00311           ii = 0;
00312           const IntegrationPoint & ip = rule[i];
00313 
00314           Sx hx = x + z * (-1+2*ip(0)+ip(1));
00315           Sy hy = y + z * ip(1);
00316 
00317           //Sx bub = hy * (1-hx-hy) * (1+hx-hy);
00318           ScaledJacobiPolynomial (n-3, hx, 1-hy, 2, 2, polx);
00319 
00320           Sx fac = 2  * ip.Weight(); //  / (z*z);
00321 
00322           for (int ix = 0; ix <= n-3; ix++)
00323             {
00324               ScaledJacobiPolynomial (n-3, 2*hy-1, 1, 2*ix+5, 2, poly);
00325               for (j = 0; j <= n-3-ix; j++)
00326                 values[ii++] += fac * polx[ix] * poly[j];
00327             }
00328         }
00329       int jj = ii;
00330       ArrayMem<Sx, 100> hvalues(jj);
00331       for (int i = 0; i < ii; i++)
00332         hvalues[i] = values[i];
00333       ii = 0; jj = 0;
00334       for (int k = 0; k <= n-3; k++)
00335         {
00336           for (int m = 0; m <= n-3-k; m++)
00337             {
00338               values[ii++] = y*(1-x-y-z)*(1+x-y-z)*hvalues[jj++];
00339             }
00340         }
00341 
00342       return ii;
00343     }
00344   };
00345 
00346 
00347 
00348 
00349 
00350 
00351 
00352 
00353 
00354 
00355   /*
00357  virtual void CalcShape (const IntegrationPoint & ip,
00358  FlatVector<> shape) const
00359  {
00360  double lam1 = ip(0);
00361  double lam2 = 1-ip(0)-ip(1);
00362  double lam3 = ip(1);
00363  double bub = lam1*lam2*lam3;
00364 
00365  ArrayMem<double, 50> polx(order-2);
00366  ArrayMem<double, 50> poly(order-2);
00367 
00368  double fac = 1;
00369  for (int i = 0; i <= order-3; i++)
00370  {
00371  polx[i] *= fac;
00372  fac *= (1-lam3);
00373  }
00374 
00375  int ii = 0;
00376  for (int ix = 0; ix <= order-3; ix++)
00377  {
00378  JacobiPolynomial (order-3, lam3-(lam1+lam2), 2*ix+5, 2, poly);
00379  for (int iy = 0; iy <= order-3-ix; iy++)
00380  shape(ii++) = bub*polx[ix] * poly[iy];
00381  }
00382  }
00383 
00384   */
00385 
00386 
00387 }