Commits

Lisandro Dalcin committed af17bf1

Update tabulation of quadrature rules

  • Participants
  • Parent commits 6077ab7

Comments (0)

Files changed (2)

src/petigaqdr.f90

      X(0) = 0.0_rk
      W(0) = 2.0_rk
   case (2) ! p <= 3
-     X(0) = -0.5773502691896257645091487805019576_rk ! 1/sqrt(3)
+     X(0) = -0.577350269189625764509148780501957455_rk ! 1/sqrt(3)
      X(1) = -X(0)
-     W(0) =  1.0_rk                                  ! 1
+     W(0) =  1.0_rk                                    ! 1
      W(1) =  W(0)
   case (3) ! p <= 5
-     X(0) = -0.7745966692414833770358530799564799_rk ! sqrt(3/5)
-     X(1) =  0.0_rk                                  ! 0
+     X(0) = -0.774596669241483377035853079956479922_rk ! sqrt(3/5)
+     X(1) =  0.0_rk                                    ! 0
      X(2) = -X(0)
-     W(0) =  0.5555555555555555555555555555555556_rk ! 5/9
-     W(1) =  0.8888888888888888888888888888888889_rk ! 8/9
+     W(0) =  0.555555555555555555555555555555555556_rk ! 5/9
+     W(1) =  0.888888888888888888888888888888888889_rk ! 8/9
      W(2) =  W(0)
   case (4) ! p <= 7
-     X(0) = -0.8611363115940525752239464888928094_rk ! sqrt((3+2*sqrt(6/5))/7)
-     X(1) = -0.3399810435848562648026657591032448_rk ! sqrt((3-2*sqrt(6/5))/7)
+     X(0) = -0.861136311594052575223946488892809506_rk ! sqrt((3+2*sqrt(6/5))/7)
+     X(1) = -0.339981043584856264802665759103244686_rk ! sqrt((3-2*sqrt(6/5))/7)
      X(2) = -X(1)
      X(3) = -X(0)
-     W(0) =  0.3478548451374538573730639492219994_rk ! (18-sqrt(30))/36
-     W(1) =  0.6521451548625461426269360507780006_rk ! (18+sqrt(30))/36
+     W(0) =  0.347854845137453857373063949221999408_rk ! (18-sqrt(30))/36
+     W(1) =  0.652145154862546142626936050778000592_rk ! (18+sqrt(30))/36
      W(2) =  W(1)
      W(3) =  W(0)
   case (5) ! p <= 9
-     X(0) = -0.9061798459386639927976268782993929_rk ! 1/3*sqrt(5+2*sqrt(10/7))
-     X(1) = -0.5384693101056830910363144207002086_rk ! 1/3*sqrt(5-2*sqrt(10/7))
-     X(2) =  0.0_rk                                  ! 0
+     X(0) = -0.906179845938663992797626878299392962_rk ! 1/3*sqrt(5+2*sqrt(10/7))
+     X(1) = -0.538469310105683091036314420700208806_rk ! 1/3*sqrt(5-2*sqrt(10/7))
+     X(2) =  0.0_rk                                    ! 0
      X(3) = -X(1)
      X(4) = -X(0)
-     W(0) =  0.2369268850561890875142640407199173_rk ! (322-13*sqrt(70))/900
-     W(1) =  0.4786286704993664680412915148356382_rk ! (322+13*sqrt(70))/900
-     W(2) =  0.5688888888888888888888888888888889_rk ! 128/225
+     W(0) =  0.236926885056189087514264040719917362_rk ! (322-13*sqrt(70))/900
+     W(1) =  0.478628670499366468041291514835638193_rk ! (322+13*sqrt(70))/900
+     W(2) =  0.568888888888888888888888888888888889_rk ! 128/225
      W(3) =  W(1)
      W(4) =  W(0)
   case (6) ! p <= 11
      W(0) =  1.0_rk
      W(1) =  W(0)
   case (3) ! p <= 3
-     X(0) = -1.0_rk                                   ! -1
-     X(1) =  0.0_rk                                   ! 0
+     X(0) = -1.0_rk                                    ! -1
+     X(1) =  0.0_rk                                    ! 0
      X(2) = -X(0)
-     W(0) =  0.33333333333333333333333333333333333_rk ! 1/3
-     W(1) =  1.33333333333333333333333333333333333_rk ! 4/3
+     W(0) =  0.333333333333333333333333333333333333_rk ! 1/3
+     W(1) =  1.333333333333333333333333333333333333_rk ! 4/3
      W(2) =  W(0)
   case (4) ! p <= 5
-     X(0) = -1.0_rk                                   ! -1
-     X(1) = -0.44721359549995793928183473374625525_rk ! 1/sqrt(5)
+     X(0) = -1.0_rk                                    ! -1
+     X(1) = -0.447213595499957939281834733746255246_rk ! 1/sqrt(5)
      X(2) = -X(1)
      X(3) = -X(0)
-     W(0) =  0.16666666666666666666666666666666667_rk ! 1/6
-     W(1) =  0.83333333333333333333333333333333333_rk ! 5/6
+     W(0) =  0.166666666666666666666666666666666667_rk ! 1/6
+     W(1) =  0.833333333333333333333333333333333343_rk ! 5/6
      W(2) =  W(1)
      W(3) =  W(0)
   case (5) ! p <= 7
-     X(0) = -1.0_rk                                   ! -1
-     X(1) = -0.65465367070797714379829245624685835_rk ! sqrt(3/7)
-     X(2) =  0.0_rk                                   ! 0
+     X(0) = -1.0_rk                                    ! -1
+     X(1) = -0.654653670707977143798292456246858356_rk ! sqrt(3/7)
+     X(2) =  0.0_rk                                    ! 0
      X(3) = -X(1)
      X(4) = -X(0)
-     W(0) =  0.1_rk                                   !  1/10
-     W(1) =  0.54444444444444444444444444444444444_rk ! 49/90
-     W(2) =  0.71111111111111111111111111111111111_rk ! 32/45
+     W(0) =  0.1_rk                                    !  1/10
+     W(1) =  0.544444444444444444444444444444444444_rk ! 49/90
+     W(2) =  0.711111111111111111111111111111111111_rk ! 32/45
      W(3) =  W(1)
      W(4) =  W(0)
   case default
   PetscFunctionReturn(0);
 }
 
+#if   defined(PETSC_USE_REAL_SINGLE)
+#define Q(constant) constant##f
+#elif defined(PETSC_USE_REAL_DOUBLE)
+#define Q(constant) constant
+#elif defined(PETSC_USE_REAL_LONG_DOUBLE)
+#define Q(constant) constant##L
+#elif   defined(PETSC_USE_REAL___FLOAT128)
+#define Q(constant) constant##Q
+#endif
+
 static PetscErrorCode GaussLegendreRule(PetscInt q, PetscReal X[], PetscReal W[])
 {
   switch (q)  {
   case (1): /* p = 1 */
-    X[0] = 0.0;
-    W[0] = 2.0;
+    X[0] =  Q(0.0);
+    W[0] =  Q(2.0);
     break;
   case (2): /* p = 3 */
-    X[0] = -0.5773502691896257645091487805019576; /* 1/sqrt(3) */
+    X[0] = -Q(0.577350269189625764509148780501957455); /* 1/sqrt(3) */
     X[1] = -X[0];
-    W[0] =  1.0;
+    W[0] =  Q(1.0);
     W[1] =  W[0];
     break;
   case (3): /* p = 5 */
-    X[0] = -0.7745966692414833770358530799564799; /* sqrt(3/5) */
-    X[1] =  0.0;
+    X[0] = -Q(0.774596669241483377035853079956479922); /* sqrt(3/5) */
+    X[1] =  Q(0.0);
     X[2] = -X[0];
-    W[0] =  0.5555555555555555555555555555555556; /* 5/9 */
-    W[1] =  0.8888888888888888888888888888888889; /* 8/9 */
+    W[0] =  Q(0.555555555555555555555555555555555556); /* 5/9 */
+    W[1] =  Q(0.888888888888888888888888888888888889); /* 8/9 */
     W[2] =  W[0];
     break;
   case (4): /* p = 7 */
-    X[0] = -0.8611363115940525752239464888928094; /* sqrt((3+2*sqrt(6/5))/7) */
-    X[1] = -0.3399810435848562648026657591032448; /* sqrt((3-2*sqrt(6/5))/7) */
+    X[0] = -Q(0.861136311594052575223946488892809506); /* sqrt((3+2*sqrt(6/5))/7) */
+    X[1] = -Q(0.339981043584856264802665759103244686); /* sqrt((3-2*sqrt(6/5))/7) */
     X[2] = -X[1];
     X[3] = -X[0];
-    W[0] =  0.3478548451374538573730639492219994; /* (18-sqrt(30))/36 */
-    W[1] =  0.6521451548625461426269360507780006; /* (18+sqrt(30))/36 */
+    W[0] =  Q(0.347854845137453857373063949221999408); /* (18-sqrt(30))/36 */
+    W[1] =  Q(0.652145154862546142626936050778000592); /* (18+sqrt(30))/36 */
     W[2] =  W[1];
     W[3] =  W[0];
     break;
   case (5): /* p = 9 */
-    X[0] = -0.9061798459386639927976268782993929; /* 1/3*sqrt(5+2*sqrt(10/7)) */
-    X[1] = -0.5384693101056830910363144207002086; /* 1/3*sqrt(5-2*sqrt(10/7)) */
-    X[2] =  0.0;
+    X[0] = -Q(0.906179845938663992797626878299392962); /* 1/3*sqrt(5+2*sqrt(10/7)) */
+    X[1] = -Q(0.538469310105683091036314420700208806); /* 1/3*sqrt(5-2*sqrt(10/7)) */
+    X[2] =  Q(0.0);
     X[3] = -X[1];
     X[4] = -X[0];
-    W[0] =  0.2369268850561890875142640407199173; /* (322-13*sqrt(70))/900 */
-    W[1] =  0.4786286704993664680412915148356382; /* (322+13*sqrt(70))/900 */
-    W[2] =  0.5688888888888888888888888888888889; /* 128/225 */
+    W[0] =  Q(0.236926885056189087514264040719917362); /* (322-13*sqrt(70))/900 */
+    W[1] =  Q(0.478628670499366468041291514835638193); /* (322+13*sqrt(70))/900 */
+    W[2] =  Q(0.568888888888888888888888888888888889); /* 128/225 */
     W[3] =  W[1];
     W[4] =  W[0];
     break;
   case (6): /* p = 11 */
-    X[0] = -0.9324695142031520278123015544939946; /* << NumericalDifferentialEquationAnalysis` */
-    X[1] = -0.6612093864662645136613995950199053; /* GaussianQuadratureWeights[6, -1, 1, 37]   */
-    X[2] = -0.2386191860831969086305017216807119;
+    X[0] = -Q(0.9324695142031520278123015544939946); /* << NumericalDifferentialEquationAnalysis` */
+    X[1] = -Q(0.6612093864662645136613995950199053); /* GaussianQuadratureWeights[6, -1, 1, 37]   */
+    X[2] = -Q(0.2386191860831969086305017216807119);
     X[3] = -X[2];
     X[4] = -X[1];
     X[5] = -X[0];
-    W[0] =  0.171324492379170345040296142172732894;
-    W[1] =  0.360761573048138607569833513837716112;
-    W[2] =  0.467913934572691047389870343989550995;
+    W[0] =  Q(0.171324492379170345040296142172732894);
+    W[1] =  Q(0.360761573048138607569833513837716112);
+    W[2] =  Q(0.467913934572691047389870343989550995);
     W[3] =  W[2];
     W[4] =  W[1];
     W[5] =  W[0];
     break;
   case (7): /* p = 13 */
-    X[0] = -0.9491079123427585245261896840478513; /* << NumericalDifferentialEquationAnalysis` */
-    X[1] = -0.7415311855993944398638647732807884; /* GaussianQuadratureWeights[7, -1, 1, 37]   */
-    X[2] = -0.4058451513773971669066064120769615;
-    X[3] =  0.0;
+    X[0] = -Q(0.9491079123427585245261896840478513); /* << NumericalDifferentialEquationAnalysis` */
+    X[1] = -Q(0.7415311855993944398638647732807884); /* GaussianQuadratureWeights[7, -1, 1, 37]   */
+    X[2] = -Q(0.4058451513773971669066064120769615);
+    X[3] =  Q(0.0);
     X[4] = -X[2];
     X[5] = -X[1];
     X[6] = -X[0];
-    W[0] =  0.129484966168869693270611432679082018;
-    W[1] =  0.279705391489276667901467771423779582;
-    W[2] =  0.381830050505118944950369775488975134;
-    W[3] =  0.417959183673469387755102040816326531;
+    W[0] =  Q(0.129484966168869693270611432679082018);
+    W[1] =  Q(0.279705391489276667901467771423779582);
+    W[2] =  Q(0.381830050505118944950369775488975134);
+    W[3] =  Q(0.417959183673469387755102040816326531);
     W[4] =  W[2];
     W[5] =  W[1];
     W[6] =  W[0];
     break;
   case (8): /* p = 15 */
-    X[0] = -0.9602898564975362316835608685694730; /* << NumericalDifferentialEquationAnalysis` */
-    X[1] = -0.7966664774136267395915539364758304; /* GaussianQuadratureWeights[8, -1, 1, 37]   */
-    X[2] = -0.5255324099163289858177390491892463;
-    X[3] = -0.1834346424956498049394761423601840;
+    X[0] = -Q(0.9602898564975362316835608685694730); /* << NumericalDifferentialEquationAnalysis` */
+    X[1] = -Q(0.7966664774136267395915539364758304); /* GaussianQuadratureWeights[8, -1, 1, 37]   */
+    X[2] = -Q(0.5255324099163289858177390491892463);
+    X[3] = -Q(0.1834346424956498049394761423601840);
     X[4] = -X[3];
     X[5] = -X[2];
     X[6] = -X[1];
     X[7] = -X[0];
-    W[0] =  0.101228536290376259152531354309962190;
-    W[1] =  0.222381034453374470544355994426240884;
-    W[2] =  0.313706645877887287337962201986601313;
-    W[3] =  0.362683783378361982965150449277195612;
+    W[0] =  Q(0.101228536290376259152531354309962190);
+    W[1] =  Q(0.222381034453374470544355994426240884);
+    W[2] =  Q(0.313706645877887287337962201986601313);
+    W[3] =  Q(0.362683783378361982965150449277195612);
     W[4] =  W[3];
     W[5] =  W[2];
     W[6] =  W[1];
     W[7] =  W[0];
     break;
   case (9): /* p = 17 */
-    X[0] = -0.9681602395076260898355762029036729; /* << NumericalDifferentialEquationAnalysis` */
-    X[1] = -0.8360311073266357942994297880697349; /* GaussianQuadratureWeights[9, -1, 1, 37]   */
-    X[2] = -0.6133714327005903973087020393414742;
-    X[3] = -0.3242534234038089290385380146433366;
-    X[4] =  0.0;
+    X[0] = -Q(0.9681602395076260898355762029036729); /* << NumericalDifferentialEquationAnalysis` */
+    X[1] = -Q(0.8360311073266357942994297880697349); /* GaussianQuadratureWeights[9, -1, 1, 37]   */
+    X[2] = -Q(0.6133714327005903973087020393414742);
+    X[3] = -Q(0.3242534234038089290385380146433366);
+    X[4] =  Q(0.0);
     X[5] = -X[3];
     X[6] = -X[2];
     X[7] = -X[1];
     X[8] = -X[0];
-    W[0] =  0.081274388361574411971892158110523651;
-    W[1] =  0.180648160694857404058472031242912810;
-    W[2] =  0.260610696402935462318742869418632850;
-    W[3] =  0.312347077040002840068630406584443666;
-    W[4] =  0.330239355001259763164525069286974049;
+    W[0] =  Q(0.081274388361574411971892158110523651);
+    W[1] =  Q(0.180648160694857404058472031242912810);
+    W[2] =  Q(0.260610696402935462318742869418632850);
+    W[3] =  Q(0.312347077040002840068630406584443666);
+    W[4] =  Q(0.330239355001259763164525069286974049);
     W[5] =  W[3];
     W[6] =  W[2];
     W[7] =  W[1];
     X[0] = -1.0;
     X[1] =  0.0;
     X[2] = -X[0];
-    W[0] =  0.3333333333333333333333333333333333; /* 1/3 */
-    W[1] =  1.3333333333333333333333333333333333; /* 4/3 */
+    W[0] =  0.33333333333333333333333333333333333; /* 1/3 */
+    W[1] =  1.33333333333333333333333333333333333; /* 4/3 */
     W[2] =  W[0];
     break;
   case (4): /* p = 5 */
     W[3] =  W[0];
     break;
   case (5): /* p = 7 */
-    X[0] = -1;
+    X[0] = -1.0;
     X[1] = -0.65465367070797714379829245624685835; /* sqrt(3/7) */
     X[2] =  0.0;
     X[3] = -X[1];