]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Reimplementation of Q1-Q4 LagrangePolynomials. New implementation of Q5-Q10 cvs diff...
authorRalf Hartmann <Ralf.Hartmann@dlr.de>
Tue, 23 Jan 2001 17:15:51 +0000 (17:15 +0000)
committerRalf Hartmann <Ralf.Hartmann@dlr.de>
Tue, 23 Jan 2001 17:15:51 +0000 (17:15 +0000)
git-svn-id: https://svn.dealii.org/trunk@3770 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/source/polynomial.cc

index 7ea2e58a8f4fdb4d99b67a6e497fc7848febc2b9..9dd7a5bd358c632ef28ac10709f464813d65046c 100644 (file)
@@ -104,130 +104,190 @@ LagrangeEquidistant::compute_coefficients (const unsigned int n,
   std::vector<double> a (n+1);
   Assert(support_point<n+1, ExcIndexRange(support_point, 0, n+1));
 
+  unsigned int n_functions=n+1;
+  Assert(support_point<n_functions,
+        ExcIndexRange(support_point, 0, n_functions));
+  double const *x=0;
+  
   switch (n)
     {
-      case 0:
-           switch (support_point)
-             {
-               case 0:
-                     a[0]=1.;
-                     break;
-               default:
-                     Assert(false, ExcInvalidSupportPoint(support_point));
-             }
-           break;
       case 1:
-           switch (support_point)
-             {
-               case 0:
-                     a[0]=1.;
-                     a[1]=-1.;
-                     break;
-               case 1:
-                     a[0]=0.;
-                     a[1]=1.;
-                     break;
-               default:
-                     Assert(false, ExcInvalidSupportPoint(support_point));
-             }
-           break;
+      {
+       static const double x1[4]=
+       { 1.0, -1.0, 0.0, 1.0 };
+       x=&x1[0];
+       break;  
+      }
       case 2:
-           switch (support_point)
-             {
-               case 0:
-                     a[0]=1.;
-                     a[1]=-3.;
-                     a[2]=2.;
-                     break;
-               case 1:
-                     a[0]=0.;
-                     a[1]=4.;                
-                     a[2]=-4.;
-                     break;
-               case 2:
-                     a[0]=0.;
-                     a[1]=-1.;               
-                     a[2]=2.;
-                     break;
-               default:
-                     Assert(false, ExcInvalidSupportPoint(support_point));
-             }
-           break;
+      {
+       static const double x2[9]=
+       { 1.0, -3.0, 2.0, 0.0, 4.0, -4.0, 0.0, -1.0, 2.0 };
+       x=&x2[0];
+       break;
+      }
       case 3:
-           switch (support_point)
-             {
-               case 0:
-                     a[0]=1.0;
-                     a[1]=-11.0/2.0;
-                     a[2]=9.0;
-                     a[3]=-9.0/2.0;
-                     break;
-               case 1:
-                     a[0]=0.;
-                     a[1]=9.0;               
-                     a[2]=-45.0/2.0;
-                     a[3]=27.0/2.0;
-                     break;
-               case 2:
-                     a[0]=0.;
-                     a[1]=-9.0/2.0;                  
-                     a[2]=18.0;
-                     a[3]=-27.0/2.0;
-                     break;
-               case 3:
-                     a[0]=0.;
-                     a[1]=1.;                
-                     a[2]=-9.0/2.0;
-                     a[3]=9.0/2.0;
-                     break;
-               default:
-                     Assert(false, ExcInvalidSupportPoint(support_point));
-             }
-           break;
+      {
+       static const double x3[16]=
+       { 1.0, -11.0/2.0, 9.0, -9.0/2.0, 0.0, 9.0, -45.0/2.0, 27.0/2.0, 0.0,
+         -9.0/2.0, 18.0, -27.0/2.0, 0.0, 1.0, -9.0/2.0, 9.0/2.0};
+       x=&x3[0];
+       break;
+      }
       case 4:
-           switch (support_point)
-             {
-               case 0:
-                     a[0]=1.;
-                     a[1]=-25.0/3.0;
-                     a[2]=70.0/3.0;
-                     a[3]=-80.0/3.0;
-                     a[4]=32.0/3.0;
-                     break;
-               case 1:
-                     a[0]=0.;
-                     a[1]=16.0;
-                     a[2]=-208.0/3.0;
-                     a[3]=96.0;
-                     a[4]=-128.0/3.0;
-                     break;
-               case 2:
-                     a[0]=0.;
-                     a[1]=-12.0;
-                     a[2]=76.0;
-                     a[3]=-128.0;
-                     a[4]=64.0;
-                     break;
-               case 3:
-                     a[0]=0.;
-                     a[1]=16.0/3.0;
-                     a[2]=-112.0/3.0;
-                     a[3]=224.0/3.0;
-                     a[4]=-128.0/3.0;
-                     break;
-               case 4:
-                     a[0]=0.;
-                     a[1]=-1.;
-                     a[2]=22.0/3.0;
-                     a[3]=-16.0;
-                     a[4]=32.0/3.0;
-                     break;
-               default:
-                     Assert(false, ExcInvalidSupportPoint(support_point));
-             }
-           break;
+      {
+       static const double x4[25]=
+       { 1.0, -25.0/3.0, 70.0/3.0, -80.0/3.0, 32.0/3.0, 0.0, 16.0,
+         -208.0/3.0, 96.0, -128.0/3.0, 0.0, -12.0, 76.0, -128.0, 64.0, 0.0,
+         16.0/3.0, -112.0/3.0, 224.0/3.0, -128.0/3.0, 0.0, -1.0, 22.0/3.0,
+         -16.0, 32.0/3.0};
+       x=&x4[0];
+       break;  
+      }
+      case 5:
+      {
+       static const double x5[36]=
+       { 1.0, -137.0/12.0, 375.0/8.0, -2125.0/24.0, 625.0/8.0, -625.0/24.0,
+         0.0, 25.0, -1925.0/12.0, 8875.0/24.0, -4375.0/12.0, 3125.0/24.0,
+         0.0, -25.0, 2675.0/12.0, -7375.0/12.0, 8125.0/12.0, -3125.0/12.0,
+         0.0, 50.0/3.0, -325.0/2.0, 6125.0/12.0, -625.0, 3125.0/12.0, 0.0,
+         -25.0/4.0, 1525.0/24.0, -5125.0/24.0, 6875.0/24.0, -3125.0/24.0,
+         0.0, 1.0, -125.0/12.0, 875.0/24.0, -625.0/12.0, 625.0/24.0};
+       x=&x5[0];
+       break;
+      }
+      case 6:
+      {
+       static const double x6[49]=
+       { 1.0, -147.0/10.0, 406.0/5.0, -441.0/2.0, 315.0, -1134.0/5.0,
+         324.0/5.0, 0.0, 36.0, -1566.0/5.0, 1044.0, -1674.0, 1296.0,
+         -1944.0/5.0, 0.0, -45.0, 1053.0/2.0, -4149.0/2.0, 3699.0, -3078.0,
+         972.0, 0.0, 40.0, -508.0, 2232.0, -4356.0, 3888.0, -1296.0, 0.0,
+         -45.0/2.0, 297.0, -2763.0/2.0, 2889.0, -2754.0, 972.0, 0.0,
+         36.0/5.0, -486.0/5.0, 468.0, -1026.0, 5184.0/5.0, -1944.0/5.0, 0.0,
+         -1.0, 137.0/10.0, -135.0/2.0, 153.0, -162.0, 324.0/5.0};
+       x=&x6[0];
+       break;
+      }
+      case 7:
+      {
+       static const double x7[64]=
+       { 1.0, -363.0/20.0, 22981.0/180.0, -331681.0/720.0, 16807.0/18.0,
+         -386561.0/360.0, 117649.0/180.0, -117649.0/720.0, 0.0, 49.0,
+         -10927.0/20.0, 109417.0/45.0, -88837.0/16.0, 991613.0/144.0,
+         -352947.0/80.0, 823543.0/720.0, 0.0, -147.0/2.0, 43071.0/40.0,
+         -1347647.0/240.0, 170471.0/12.0, -151263.0/8.0, 1529437.0/120.0,
+         -823543.0/240.0, 0.0, 245.0/3.0, -46501.0/36.0, 133427.0/18.0,
+         -2926819.0/144.0, 4151329.0/144.0, -2941225.0/144.0,
+         823543.0/144.0, 0.0, -245.0/4.0, 2009.0/2.0, -872935.0/144.0,
+         52822.0/3.0, -1899191.0/72.0, 117649.0/6.0, -823543.0/144.0, 0.0,
+         147.0/5.0, -9849.0/20.0, 45962.0/15.0, -444185.0/48.0,
+         1159683.0/80.0, -2705927.0/240.0, 823543.0/240.0, 0.0, -49.0/6.0,
+         49931.0/360.0, -634207.0/720.0, 98441.0/36.0, -319333.0/72.0,
+         1294139.0/360.0, -823543.0/720.0, 0.0, 1.0, -343.0/20.0,
+         9947.0/90.0, -16807.0/48.0, 84035.0/144.0, -117649.0/240.0,
+         117649.0/720.0};
+       x=&x7[0];
+       break;
+      }
+      case 8:
+      {
+       static const double x8[81]=
+       { 1.0, -761.0/35.0, 59062.0/315.0, -4272.0/5.0, 34208.0/15.0,
+         -18432.0/5.0, 53248.0/15.0, -65536.0/35.0, 131072.0/315.0, 0.0,
+         64.0, -30784.0/35.0, 44672.0/9.0, -673792.0/45.0, 235520.0/9.0,
+         -1196032.0/45.0, 131072.0/9.0, -1048576.0/315.0, 0.0, -112.0,
+         9936.0/5.0, -587296.0/45.0, 1956992.0/45.0, -733184.0/9.0,
+         3915776.0/45.0, -2228224.0/45.0, 524288.0/45.0, 0.0, 448.0/3.0,
+         -128192.0/45.0, 102016.0/5.0, -1097728.0/15.0, 145408.0,
+         -2441216.0/15.0, 1441792.0/15.0, -1048576.0/45.0, 0.0, -140.0,
+         2764.0, -186496.0/9.0, 703552.0/9.0, -1466368.0/9.0, 1712128.0/9.0,
+         -1048576.0/9.0, 262144.0/9.0, 0.0, 448.0/5.0, -9024.0/5.0,
+         626048.0/45.0, -2443264.0/45.0, 5285888.0/45.0, -6406144.0/45.0,
+         4063232.0/45.0, -1048576.0/45.0, 0.0, -112.0/3.0, 34288.0/45.0,
+         -5984.0, 358784.0/15.0, -53248.0, 999424.0/15.0, -131072.0/3.0,
+         524288.0/45.0, 0.0, 64.0/7.0, -6592.0/35.0, 67456.0/45.0,
+         -274432.0/45.0, 124928.0/9.0, -802816.0/45.0, 3801088.0/315.0,
+         -1048576.0/315.0, 0.0, -1.0, 726.0/35.0, -7504.0/45.0, 30944.0/45.0,
+         -14336.0/9.0, 94208.0/45.0, -65536.0/45.0, 131072.0/315.0};
+       x=&x8[0];
+       break;
+      }
+      case 9:
+      {
+       static const double x9[100]=
+       { 1.0, -7129.0/280.0, 58635.0/224.0, -40707.0/28.0, 623295.0/128.0,
+         -6589431.0/640.0, 885735.0/64.0, -5137263.0/448.0, 4782969.0/896.0,
+         -4782969.0/4480.0, 0.0, 81.0, -373329.0/280.0, 10307331.0/1120.0,
+         -5589243.0/160.0, 51221727.0/640.0, -4546773.0/40.0,
+         31355019.0/320.0, -52612659.0/1120.0, 43046721.0/4480.0, 0.0,
+         -162.0, 475389.0/140.0, -15190173.0/560.0, 18152829.0/160.0,
+         -44529507.0/160.0, 33244587.0/80.0, -3720087.0/10.0,
+         205667667.0/1120.0, -43046721.0/1120.0, 0.0, 252.0, -56601.0/10.0,
+         1959363.0/40.0, -8776431.0/40.0, 91020753.0/160.0,
+         -71035947.0/80.0, 16474671.0/20.0, -33480783.0/80.0,
+         14348907.0/160.0, 0.0, -567.0/2.0, 526419.0/80.0, -4752351.0/80.0,
+         89119521.0/320.0, -241241409.0/320.0, 195629337.0/160.0,
+         -187598673.0/160.0, 196101729.0/320.0, -43046721.0/320.0, 0.0,
+         1134.0/5.0, -21465.0/4.0, 795339.0/16.0, -3844017.0/16.0,
+         215023653.0/320.0, -18009945.0/16.0, 35606547.0/32.0,
+         -4782969.0/8.0, 43046721.0/320.0, 0.0, -126.0, 60381.0/20.0,
+         -2276289.0/80.0, 22480173.0/160.0, -64448703.0/160.0,
+         55447011.0/80.0, -28166373.0/40.0, 62178597.0/160.0,
+         -14348907.0/160.0, 0.0, 324.0/7.0, -78327.0/70.0, 2989629.0/280.0,
+         -2142531.0/40.0, 25043337.0/160.0, -22025277.0/80.0,
+         80247591.0/280.0, -90876411.0/560.0, 43046721.0/1120.0, 0.0,
+         -81.0/8.0, 275967.0/1120.0, -1328967.0/560.0, 7712091.0/640.0,
+         -22878207.0/640.0, 20490003.0/320.0, -21789081.0/320.0,
+         176969853.0/4480.0, -43046721.0/4480.0, 0.0, 1.0, -6849.0/280.0,
+         265779.0/1120.0, -194643.0/160.0, 2337903.0/640.0, -531441.0/80.0,
+         2302911.0/320.0, -4782969.0/1120.0, 4782969.0/4480.0};
+       x=&x9[0];
+       break;
+      }
+      case 10:
+      {
+       static const double x10[121]=
+       { 1.0, -7381.0/252.0, 177133.0/504.0, -10511875.0/4536.0,
+         42711625.0/4536.0, -5369375.0/216.0, 4695625.0/108.0,
+         -9453125.0/189.0, 6875000.0/189.0, -8593750.0/567.0,
+         1562500.0/567.0, 0.0, 100.0, -121525.0/63.0, 1997825.0/126.0,
+         -82992625.0/1134.0, 3775625.0/18.0, -20965625.0/54.0,
+         4187500.0/9.0, -65937500.0/189.0, 3125000.0/21.0,
+         -15625000.0/567.0, 0.0, -225.0, 153025.0/28.0, -2898075.0/56.0,
+         33095875.0/126.0, -57981875.0/72.0, 56396875.0/36.0,
+         -17546875.0/9.0, 94843750.0/63.0, -41406250.0/63.0, 7812500.0/63.0,
+         0.0, 400.0, -654100.0/63.0, 20028950.0/189.0, -108434750.0/189.0,
+         16686250.0/9.0, -33868750.0/9.0, 43625000.0/9.0, -242500000.0/63.0,
+         325000000.0/189.0, -62500000.0/189.0, 0.0, -525.0, 168775.0/12.0,
+         -1792225.0/12.0, 91073375.0/108.0, -102070625.0/36.0,
+         107321875.0/18.0, -71281250.0/9.0, 19375000.0/3.0, -26562500.0/9.0,
+         15625000.0/27.0, 0.0, 504.0, -13754.0, 149625.0, -7818625.0/9.0,
+         27074375.0/9.0, -58608125.0/9.0, 80000000.0/9.0, -66875000.0/9.0,
+         31250000.0/9.0, -6250000.0/9.0, 0.0, -350.0, 174025.0/18.0,
+         -11544725.0/108.0, 34178875.0/54.0, -80666875.0/36.0,
+         89384375.0/18.0, -62468750.0/9.0, 5937500.0, -76562500.0/27.0,
+         15625000.0/27.0, 0.0, 1200.0/7.0, -100300.0/21.0, 1121950.0/21.0,
+         -60659750.0/189.0, 10401250.0/9.0, -7831250.0/3.0,
+         234625000.0/63.0, -205000000.0/63.0, 100000000.0/63.0,
+         -62500000.0/189.0, 0.0, -225.0/4.0, 88325.0/56.0, -996675.0/56.0,
+         54486625.0/504.0, -28405625.0/72.0, 32584375.0/36.0,
+         -11828125.0/9.0, 73750000.0/63.0, -36718750.0/63.0, 7812500.0/63.0,
+         0.0, 100.0/9.0, -6575.0/21.0, 4033825.0/1134.0, -24717625.0/1134.0,
+         4341875.0/54.0, -10090625.0/54.0, 7437500.0/27.0,
+         -47187500.0/189.0, 71875000.0/567.0, -15625000.0/567.0, 0.0, -1.0,
+         7129.0/252.0, -162875.0/504.0, 1130750.0/567.0, -59375.0/8.0,
+         1883125.0/108.0, -78125.0/3.0, 4531250.0/189.0, -781250.0/63.0,
+         1562500.0/567.0};
+       x=&x10[0];
+       break;
+      }
       default:
            Assert(false, ExcNotImplemented());
     }
+
+  Assert(x!=0, ExcInternalError());
+  for (unsigned int i=0; i<n_functions; ++i)
+    a[i]=x[support_point*n_functions+i];
+  
   return a;
 }

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.