[MaL] fixed Gauss quadrature for Tet 2rd order

MathLib/Integration/GaussLegendreTet.cpp

@@ -25,37 +25,41 @@ GaussLegendreTet<2>::X = {{ {{1./4., 1./4., 1./4.}},
                             {{1./6., 1./6., 1./2.}} }};
 double const GaussLegendreTet<2>::W[5] = {-2./15., 0.075, 0.075, 0.075, 0.075};
 
+std::array<std::array<double, 3>, GaussLegendreTet<3>::NPoints> initGLTet3X()
+{
+    // Cf. Gellert, M., Harbord, R., 1991. Moderate degree cubature formulas for
+    // 3-D tetrahedral finite-element approximations. Communications in Applied
+    // Numerical Methods 7, 487–495. doi:10.1002/cnm.1630070609
+    const double a = 0.0673422422100983;
+    const double b = 0.3108859192633005;
+    const double c = 0.7217942490673264;
+    const double d = 0.0927352503108912;
+    const double e = 0.4544962958743506;
+    const double f = 0.045503704125649;
+
+    return {{{a, b, b},
+             {b, a, b},
+             {b, b, a},
+             {b, b, b},
+             {c, d, d},
+             {d, c, d},
+             {d, d, c},
+             {d, d, d},
+             {e, e, f},
+             {e, f, e},
+             {e, f, f},
+             {f, e, e},
+             {f, e, f},
+             {f, f, e}}};
+}
+
 const std::array<std::array<double, 3>, GaussLegendreTet<3>::NPoints>
-GaussLegendreTet<3>::X = {{ {{1./2., 1./2., 1./2.}},
-                            {{0.09197107805272303, 0.09197107805272303, 0.09197107805272303}},
-                            {{0.72408676584183096, 0.09197107805272303, 0.09197107805272303}},
-                            {{0.09197107805272303, 0.72408676584183096, 0.09197107805272303}},
-                            {{0.09197107805272303, 0.09197107805272303, 0.72408676584183096}},
-                            {{0.44364916731037080, 0.05635083268962915, 0.05635083268962915}},
-                            {{0.05635083268962915, 0.44364916731037080, 0.05635083268962915}},
-                            {{0.05635083268962915, 0.05635083268962915, 0.44364916731037080}},
-                            {{0.05635083268962915, 0.44364916731037080, 0.44364916731037080}},
-                            {{0.44364916731037080, 0.05635083268962915, 0.44364916731037080}},
-                            {{0.44364916731037080, 0.44364916731037080, 0.05635083268962915}},
-                            {{0.31979362782962989, 0.31979362782962989, 0.31979362782962989}},
-                            {{0.04061911651111023, 0.31979362782962989, 0.31979362782962989}},
-                            {{0.31979362782962989, 0.04061911651111023, 0.31979362782962989}},
-                            {{0.31979362782962989, 0.31979362782962989, 0.04061911651111023}} }};
-double const GaussLegendreTet<3>::W[15] = {0.019753086419753086,
-                                           0.011989513963169772,
-                                           0.011989513963169772,
-                                           0.011989513963169772,
-                                           0.011989513963169772,
-                                           0.008818342151675485,
-                                           0.008818342151675485,
-                                           0.008818342151675485,
-                                           0.008818342151675485,
-                                           0.008818342151675485,
-                                           0.008818342151675485,
-                                           0.011511367871045397,
-                                           0.011511367871045397,
-                                           0.011511367871045397,
-                                           0.011511367871045397
-                                           };
+    GaussLegendreTet<3>::X = initGLTet3X();
+
+static const double p = 0.1126879257180162 / 6.;
+static const double q = 0.0734930431163619 / 6.;
+static const double r = 0.0425460207770812 / 6.;
 
+double const GaussLegendreTet<3>::W[GaussLegendreTet<3>::NPoints] = {
+    p, p, p, p, q, q, q, q, r, r, r, r, r, r};
 }

MathLib/Integration/GaussLegendreTet.h

@@ -38,7 +38,7 @@ struct GaussLegendreTet<2> {
 template <>
 struct GaussLegendreTet<3> {
     static MATHLIB_EXPORT const unsigned Order = 3;
-    static MATHLIB_EXPORT const unsigned NPoints = 15;
+    static MATHLIB_EXPORT const unsigned NPoints = 14;
     static MATHLIB_EXPORT const std::array<std::array<double, 3>, NPoints> X;
     static MATHLIB_EXPORT const double W[NPoints];
 };