diff --git a/src/analysis/SimplicialGaussLegendreQuadrature.cpp b/src/analysis/SimplicialGaussLegendreQuadrature.cpp
index c6143bdec0f2013660378d36625788e21408c7a7..48843282c27922109de9716778a8eee478883f9f 100644
--- a/src/analysis/SimplicialGaussLegendreQuadrature.cpp
+++ b/src/analysis/SimplicialGaussLegendreQuadrature.cpp
@@ -196,5 +196,273 @@ template <>
void
SimplicialGaussLegendreQuadrature<3>::_buildPointAndWeightLists()
{
- throw NotImplementedError("gauss legendre quadrature on 3d simplex");
+ switch (m_order) {
+ case 0:
+ case 1: {
+ constexpr size_t nb_points = 1;
+ SmallArray<TinyVector<3>> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ point_list[0] = {-0.5, -0.5, -0.5};
+ weight_list[0] = 1.33333333333333;
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 2: {
+ constexpr size_t nb_points = 4;
+ SmallArray<TinyVector<3>> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ point_list[0] = {-0.7236067977, -0.7236067977, -0.7236067977};
+ point_list[1] = {+0.1708203932, -0.7236067977, -0.7236067977};
+ point_list[2] = {-0.7236067977, +0.1708203932, -0.7236067977};
+ point_list[3] = {-0.7236067977, -0.7236067977, +0.1708203932};
+
+ weight_list[0] = 0.333333333333333;
+ weight_list[1] = 0.333333333333333;
+ weight_list[2] = 0.333333333333333;
+ weight_list[3] = 0.333333333333333;
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 3: {
+ constexpr size_t nb_points = 5;
+ SmallArray<TinyVector<3>> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ point_list[0] = {-0.5000000000, -0.5000000000, -0.5000000000};
+ point_list[1] = {-0.66666666666666666666, -0.66666666666666666666, -0.66666666666666666666};
+ point_list[2] = {-0.66666666666666666666, -0.66666666666666666666, +0.0000000000};
+ point_list[3] = {-0.66666666666666666666, +0.0000000000, -0.66666666666666666666};
+ point_list[4] = {+0.0000000000, -0.66666666666666666666, -0.66666666666666666666};
+
+ weight_list[0] = -1.066666666666666;
+ weight_list[1] = +0.600000000000000;
+ weight_list[2] = +0.600000000000000;
+ weight_list[3] = +0.600000000000000;
+ weight_list[4] = +0.600000000000000;
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 4: {
+ constexpr size_t nb_points = 11;
+ SmallArray<TinyVector<3>> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ point_list[0] = {-0.5000000000, -0.5000000000, -0.5000000000};
+ point_list[1] = {-0.8571428571, -0.8571428571, -0.8571428571};
+ point_list[2] = {-0.8571428571, -0.8571428571, +0.5714285714};
+ point_list[3] = {-0.8571428571, +0.5714285714, -0.8571428571};
+ point_list[4] = {+0.5714285714, -0.8571428571, -0.8571428571};
+ point_list[5] = {-0.2011928476, -0.2011928476, -0.7988071523};
+ point_list[6] = {-0.2011928476, -0.7988071523, -0.2011928476};
+ point_list[7] = {-0.7988071523, -0.2011928476, -0.2011928476};
+ point_list[8] = {-0.2011928476, -0.7988071523, -0.7988071523};
+ point_list[9] = {-0.7988071523, -0.2011928476, -0.7988071523};
+ point_list[10] = {-0.7988071523, -0.7988071523, -0.2011928476};
+
+ weight_list[0] = -0.105244444444444;
+ weight_list[1] = +0.060977777777777;
+ weight_list[2] = +0.060977777777777;
+ weight_list[3] = +0.060977777777777;
+ weight_list[4] = +0.060977777777777;
+ weight_list[5] = +0.199111111111111;
+ weight_list[6] = +0.199111111111111;
+ weight_list[7] = +0.199111111111111;
+ weight_list[8] = +0.199111111111111;
+ weight_list[9] = +0.199111111111111;
+ weight_list[10] = +0.199111111111111;
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 5: {
+ constexpr size_t nb_points = 14;
+ SmallArray<TinyVector<3>> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ point_list[0] = {-0.8145294993, -0.8145294993, -0.8145294993};
+ point_list[1] = {+0.4435884981, -0.8145294993, -0.8145294993};
+ point_list[2] = {-0.8145294993, +0.4435884981, -0.8145294993};
+ point_list[3] = {-0.8145294993, -0.8145294993, +0.4435884981};
+ point_list[4] = {-0.3782281614, -0.3782281614, -0.3782281614};
+ point_list[5] = {-0.8653155155, -0.3782281614, -0.3782281614};
+ point_list[6] = {-0.3782281614, -0.8653155155, -0.3782281614};
+ point_list[7] = {-0.3782281614, -0.3782281614, -0.8653155155};
+ point_list[8] = {-0.0910074082, -0.0910074082, -0.9089925917};
+ point_list[9] = {-0.0910074082, -0.9089925917, -0.0910074082};
+ point_list[10] = {-0.9089925917, -0.0910074082, -0.0910074082};
+ point_list[11] = {-0.0910074082, -0.9089925917, -0.9089925917};
+ point_list[12] = {-0.9089925917, -0.0910074082, -0.9089925917};
+ point_list[13] = {-0.9089925917, -0.9089925917, -0.0910074082};
+
+ weight_list[0] = 0.0979907241;
+ weight_list[1] = 0.0979907241;
+ weight_list[2] = 0.0979907241;
+ weight_list[3] = 0.0979907241;
+ weight_list[4] = 0.1502505676;
+ weight_list[5] = 0.1502505676;
+ weight_list[6] = 0.1502505676;
+ weight_list[7] = 0.1502505676;
+ weight_list[8] = 0.0567280277;
+ weight_list[9] = 0.0567280277;
+ weight_list[10] = 0.0567280277;
+ weight_list[11] = 0.0567280277;
+ weight_list[12] = 0.0567280277;
+ weight_list[13] = 0.0567280277;
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 6: {
+ constexpr size_t nb_points = 24;
+ SmallArray<TinyVector<3>> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ point_list[0] = {-0.5707942574, -0.5707942574, -0.5707942574};
+ point_list[1] = {-0.2876172275, -0.5707942574, -0.5707942574};
+ point_list[2] = {-0.5707942574, -0.2876172275, -0.5707942574};
+ point_list[3] = {-0.5707942574, -0.5707942574, -0.2876172275};
+ point_list[4] = {-0.9186520829, -0.9186520829, -0.9186520829};
+ point_list[5] = {+0.7559562487, -0.9186520829, -0.9186520829};
+ point_list[6] = {-0.9186520829, +0.7559562487, -0.9186520829};
+ point_list[7] = {-0.9186520829, -0.9186520829, +0.7559562487};
+ point_list[8] = {-0.3553242197, -0.3553242197, -0.3553242197};
+ point_list[9] = {-0.9340273408, -0.3553242197, -0.3553242197};
+ point_list[10] = {-0.3553242197, -0.9340273408, -0.3553242197};
+ point_list[11] = {-0.3553242197, -0.3553242197, -0.9340273408};
+ point_list[12] = {-0.8726779962, -0.8726779962, -0.4606553370};
+ point_list[13] = {-0.8726779962, -0.4606553370, -0.8726779962};
+ point_list[14] = {-0.8726779962, -0.8726779962, +0.2060113295};
+ point_list[15] = {-0.8726779962, +0.2060113295, -0.8726779962};
+ point_list[16] = {-0.8726779962, -0.4606553370, +0.2060113295};
+ point_list[17] = {-0.8726779962, +0.2060113295, -0.4606553370};
+ point_list[18] = {-0.4606553370, -0.8726779962, -0.8726779962};
+ point_list[19] = {-0.4606553370, -0.8726779962, +0.2060113295};
+ point_list[20] = {-0.4606553370, +0.2060113295, -0.8726779962};
+ point_list[21] = {+0.2060113295, -0.8726779962, -0.4606553370};
+ point_list[22] = {+0.2060113295, -0.8726779962, -0.8726779962};
+ point_list[23] = {+0.2060113295, -0.4606553370, -0.8726779962};
+
+ weight_list[0] = 0.0532303336;
+ weight_list[1] = 0.0532303336;
+ weight_list[2] = 0.0532303336;
+ weight_list[3] = 0.0532303336;
+ weight_list[4] = 0.0134362814;
+ weight_list[5] = 0.0134362814;
+ weight_list[6] = 0.0134362814;
+ weight_list[7] = 0.0134362814;
+ weight_list[8] = 0.0738095753;
+ weight_list[9] = 0.0738095753;
+ weight_list[10] = 0.0738095753;
+ weight_list[11] = 0.0738095753;
+ weight_list[12] = 0.0642857142;
+ weight_list[13] = 0.0642857142;
+ weight_list[14] = 0.0642857142;
+ weight_list[15] = 0.0642857142;
+ weight_list[16] = 0.0642857142;
+ weight_list[17] = 0.0642857142;
+ weight_list[18] = 0.0642857142;
+ weight_list[19] = 0.0642857142;
+ weight_list[20] = 0.0642857142;
+ weight_list[21] = 0.0642857142;
+ weight_list[22] = 0.0642857142;
+ weight_list[23] = 0.0642857142;
+ weight_list[23] = 0.0642857142;
+ weight_list[23] = 0.0642857142;
+ weight_list[23] = 0.0642857142;
+ weight_list[23] = 0.0642857142;
+ weight_list[23] = 0.0642857142;
+ weight_list[23] = 0.0642857142;
+ weight_list[23] = 0.0642857142;
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 7: {
+ constexpr size_t nb_points = 31;
+ SmallArray<TinyVector<3>> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ point_list[0] = {+0.0000000000, +0.0000000000, -1.0000000000};
+ point_list[1] = {+0.0000000000, -1.0000000000, +0.0000000000};
+ point_list[2] = {-1.0000000000, +0.0000000000, +0.0000000000};
+ point_list[3] = {-1.0000000000, -1.0000000000, +0.0000000000};
+ point_list[4] = {-1.0000000000, +0.0000000000, -1.0000000000};
+ point_list[5] = {+0.0000000000, -1.0000000000, -1.0000000000};
+ point_list[6] = {-0.5000000000, -0.5000000000, -0.5000000000};
+ point_list[7] = {-0.8435736153, -0.8435736153, -0.8435736153};
+ point_list[8] = {-0.8435736153, -0.8435736153, +0.5307208460};
+ point_list[9] = {-0.8435736153, +0.5307208460, -0.8435736153};
+ point_list[10] = {+0.5307208460, -0.8435736153, -0.8435736153};
+ point_list[11] = {-0.7563135666, -0.7563135666, -0.7563135666};
+ point_list[12] = {-0.7563135666, -0.7563135666, +0.2689407000};
+ point_list[13] = {-0.7563135666, +0.2689407000, -0.7563135666};
+ point_list[14] = {+0.2689407000, -0.7563135666, -0.7563135666};
+ point_list[15] = {-0.3349216711, -0.3349216711, -0.3349216711};
+ point_list[16] = {-0.3349216711, -0.3349216711, -0.9952349866};
+ point_list[17] = {-0.3349216711, -0.9952349866, -0.3349216711};
+ point_list[18] = {-0.9952349866, -0.3349216711, -0.3349216711};
+ point_list[19] = {-0.8000000000, -0.8000000000, -0.6000000000};
+ point_list[20] = {-0.8000000000, -0.6000000000, -0.8000000000};
+ point_list[21] = {-0.8000000000, -0.8000000000, +0.2000000000};
+ point_list[22] = {-0.8000000000, +0.2000000000, -0.8000000000};
+ point_list[23] = {-0.8000000000, -0.6000000000, +0.2000000000};
+ point_list[24] = {-0.8000000000, +0.2000000000, -0.6000000000};
+ point_list[25] = {-0.6000000000, -0.8000000000, -0.8000000000};
+ point_list[26] = {-0.6000000000, -0.8000000000, +0.2000000000};
+ point_list[27] = {-0.6000000000, +0.2000000000, -0.8000000000};
+ point_list[28] = {+0.2000000000, -0.8000000000, -0.6000000000};
+ point_list[29] = {+0.2000000000, -0.8000000000, -0.8000000000};
+ point_list[30] = {+0.2000000000, -0.6000000000, -0.8000000000};
+
+ weight_list[0] = +0.0077601410;
+ weight_list[1] = +0.0077601410;
+ weight_list[2] = +0.0077601410;
+ weight_list[3] = +0.0077601410;
+ weight_list[4] = +0.0077601410;
+ weight_list[5] = +0.0077601410;
+ weight_list[6] = +0.1461137877;
+ weight_list[7] = +0.0847995321;
+ weight_list[8] = +0.0847995321;
+ weight_list[9] = +0.0847995321;
+ weight_list[10] = +0.0847995321;
+ weight_list[11] = -0.5001419209;
+ weight_list[12] = -0.5001419209;
+ weight_list[13] = -0.5001419209;
+ weight_list[14] = -0.5001419209;
+ weight_list[15] = +0.0391314021;
+ weight_list[16] = +0.0391314021;
+ weight_list[17] = +0.0391314021;
+ weight_list[18] = +0.0391314021;
+ weight_list[19] = +0.2204585537;
+ weight_list[20] = +0.2204585537;
+ weight_list[21] = +0.2204585537;
+ weight_list[22] = +0.2204585537;
+ weight_list[23] = +0.2204585537;
+ weight_list[24] = +0.2204585537;
+ weight_list[25] = +0.2204585537;
+ weight_list[26] = +0.2204585537;
+ weight_list[27] = +0.2204585537;
+ weight_list[28] = +0.2204585537;
+ weight_list[29] = +0.2204585537;
+ weight_list[30] = +0.2204585537;
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ default: {
+ throw NormalError("Gauss-Legendre quadrature formulae handle orders up to 7 on tetrahedra");
+ }
+ }
}
diff --git a/tests/test_SimplicialGaussLegendreQuadrature.cpp b/tests/test_SimplicialGaussLegendreQuadrature.cpp
index 9517a9f1fd2326424fac92ffbbdf91b1f0ea4efa..a4334e6ad958835441972a70c99bb571875a64bb 100644
--- a/tests/test_SimplicialGaussLegendreQuadrature.cpp
+++ b/tests/test_SimplicialGaussLegendreQuadrature.cpp
@@ -725,4 +725,225 @@ TEST_CASE("SimplicialGaussLegendreQuadrature", "[analysis]")
"error: Gauss-Legendre quadrature formulae handle orders up to 7 on triangles");
}
}
+
+ SECTION("3D")
+ {
+ auto integrate = [](auto f, auto quadrature_formula) {
+ auto point_list = quadrature_formula.pointList();
+ auto weight_list = quadrature_formula.weightList();
+
+ auto value = weight_list[0] * f(point_list[0]);
+ for (size_t i = 1; i < weight_list.size(); ++i) {
+ value += weight_list[i] * f(point_list[i]);
+ }
+
+ return value;
+ };
+
+ auto integrate_on_tetra = [](auto f, auto quadrature_formula, const std::array<TinyVector<3>, 4>& tetrahedron) {
+ AffineTransformation<3> t{tetrahedron};
+
+ auto point_list = quadrature_formula.pointList();
+ auto weight_list = quadrature_formula.weightList();
+
+ auto value = weight_list[0] * f(t(point_list[0]));
+ for (size_t i = 1; i < weight_list.size(); ++i) {
+ value += weight_list[i] * f(t(point_list[i]));
+ }
+
+ return t.jacobianDeterminant() * value;
+ };
+
+ auto get_order = [&integrate, &integrate_on_tetra](auto f, auto quadrature_formula, const double exact_value) {
+ using R3 = TinyVector<3>;
+
+ const R3 A{-1, -1, -1};
+ const R3 B{+1, -1, -1};
+ const R3 C{-1, +1, -1};
+ const R3 D{-1, -1, +1};
+
+ const R3 M_AB = 0.5 * (A + B);
+ const R3 M_AC = 0.5 * (A + C);
+ const R3 M_AD = 0.5 * (A + D);
+ const R3 M_BC = 0.5 * (B + C);
+ const R3 M_BD = 0.5 * (B + D);
+ const R3 M_CD = 0.5 * (C + D);
+
+ const double int_T_hat = integrate(f, quadrature_formula);
+ const double int_refined //
+ = integrate_on_tetra(f, quadrature_formula, {A, M_AB, M_AC, M_AD}) +
+ integrate_on_tetra(f, quadrature_formula, {B, M_AB, M_BD, M_BC}) +
+ integrate_on_tetra(f, quadrature_formula, {C, M_AC, M_BC, M_CD}) +
+ integrate_on_tetra(f, quadrature_formula, {D, M_AD, M_CD, M_BD}) +
+ integrate_on_tetra(f, quadrature_formula, {M_BD, M_AC, M_AD, M_CD}) +
+ integrate_on_tetra(f, quadrature_formula, {M_AB, M_AC, M_AD, M_BD}) +
+ integrate_on_tetra(f, quadrature_formula, {M_BC, M_AB, M_BD, M_AC}) +
+ integrate_on_tetra(f, quadrature_formula, {M_BC, M_AC, M_BD, M_CD});
+
+ return -std::log((int_refined - exact_value) / (int_T_hat - exact_value)) / std::log(2);
+ };
+
+ auto p0 = [](const TinyVector<3>&) { return 4; };
+ auto p1 = [](const TinyVector<3>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return 2 * x + 3 * y + z - 1;
+ };
+ auto p2 = [&p1](const TinyVector<3>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return p1(X) * (2.5 * x - 3 * y + z + 3);
+ };
+ auto p3 = [&p2](const TinyVector<3>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return p2(X) * (3 * x + 2 * y - 3 * z - 1);
+ };
+ auto p4 = [&p3](const TinyVector<3>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return p3(X) * (2 * x - 0.5 * y - 1.3 * z + 1);
+ };
+ auto p5 = [&p4](const TinyVector<3>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return p4(X) * (-0.1 * x + 1.3 * y - 3 * z + 1);
+ };
+ auto p6 = [&p5](const TinyVector<3>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return p5(X) * 7875. / 143443 * (2 * x - y + 4 * z + 1);
+ };
+ auto p7 = [&p6](const TinyVector<3>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return p6(X) * (0.7 * x - 2.7 * y + 1.3 * z - 2);
+ };
+ auto p8 = [&p7](const TinyVector<3>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return p7(X) * (0.3 * x + 1.2 * y - 0.7 * z + 0.2);
+ };
+
+ SECTION("order 2")
+ {
+ SimplicialGaussLegendreQuadrature<3> l1(1);
+
+ REQUIRE(l1.numberOfPoints() == 1);
+
+ REQUIRE(integrate(p0, l1) == Catch::Approx(16. / 3));
+ REQUIRE(integrate(p1, l1) == Catch::Approx(-16. / 3));
+ REQUIRE(integrate(p2, l1) != Catch::Approx(-47. / 3));
+
+ REQUIRE(get_order(p2, l1, -47. / 3) >= Catch::Approx(2));
+ }
+
+ SECTION("order 3")
+ {
+ SimplicialGaussLegendreQuadrature<3> l2(2);
+
+ REQUIRE(l2.numberOfPoints() == 4);
+
+ REQUIRE(integrate(p0, l2) == Catch::Approx(16. / 3));
+ REQUIRE(integrate(p1, l2) == Catch::Approx(-16. / 3));
+ REQUIRE(integrate(p2, l2) == Catch::Approx(-47. / 3));
+ REQUIRE(integrate(p3, l2) != Catch::Approx(557. / 15));
+
+ REQUIRE(get_order(p3, l2, 557. / 15) >= Catch::Approx(3));
+ }
+
+ SECTION("order 4")
+ {
+ SimplicialGaussLegendreQuadrature<3> l3(3);
+
+ REQUIRE(l3.numberOfPoints() == 5);
+
+ REQUIRE(integrate(p0, l3) == Catch::Approx(16. / 3));
+ REQUIRE(integrate(p1, l3) == Catch::Approx(-16. / 3));
+ REQUIRE(integrate(p2, l3) == Catch::Approx(-47. / 3));
+ REQUIRE(integrate(p3, l3) == Catch::Approx(557. / 15));
+ REQUIRE(integrate(p4, l3) != Catch::Approx(4499. / 1575));
+ }
+
+ SECTION("order 5")
+ {
+ SimplicialGaussLegendreQuadrature<3> l4(4);
+
+ REQUIRE(l4.numberOfPoints() == 11);
+
+ REQUIRE(integrate(p0, l4) == Catch::Approx(16. / 3));
+ REQUIRE(integrate(p1, l4) == Catch::Approx(-16. / 3));
+ REQUIRE(integrate(p2, l4) == Catch::Approx(-47. / 3));
+ REQUIRE(integrate(p3, l4) == Catch::Approx(557. / 15));
+ REQUIRE(integrate(p4, l4) == Catch::Approx(4499. / 1575));
+ REQUIRE(integrate(p5, l4) != Catch::Approx(143443. / 7875));
+
+ REQUIRE(get_order(p5, l4, 143443. / 7875) >= Catch::Approx(5));
+ }
+
+ SECTION("order 6")
+ {
+ SimplicialGaussLegendreQuadrature<3> l5(5);
+
+ REQUIRE(l5.numberOfPoints() == 14);
+
+ REQUIRE(integrate(p0, l5) == Catch::Approx(16. / 3));
+ REQUIRE(integrate(p1, l5) == Catch::Approx(-16. / 3));
+ REQUIRE(integrate(p2, l5) == Catch::Approx(-47. / 3));
+ REQUIRE(integrate(p3, l5) == Catch::Approx(557. / 15));
+ REQUIRE(integrate(p4, l5) == Catch::Approx(4499. / 1575));
+ REQUIRE(integrate(p5, l5) == Catch::Approx(143443. / 7875));
+ REQUIRE(integrate(p6, l5) != Catch::Approx(-75773. / 2581974));
+ }
+
+ SECTION("order 7")
+ {
+ SimplicialGaussLegendreQuadrature<3> l6(6);
+
+ REQUIRE(l6.numberOfPoints() == 24);
+
+ REQUIRE(integrate(p0, l6) == Catch::Approx(16. / 3));
+ REQUIRE(integrate(p1, l6) == Catch::Approx(-16. / 3));
+ REQUIRE(integrate(p2, l6) == Catch::Approx(-47. / 3));
+ REQUIRE(integrate(p3, l6) == Catch::Approx(557. / 15));
+ REQUIRE(integrate(p4, l6) == Catch::Approx(4499. / 1575));
+ REQUIRE(integrate(p5, l6) == Catch::Approx(143443. / 7875));
+ REQUIRE(integrate(p6, l6) == Catch::Approx(-75773. / 2581974));
+ REQUIRE(integrate(p7, l6) != Catch::Approx(86951548. / 32274675));
+
+ REQUIRE(get_order(p7, l6, 86951548. / 32274675) >= Catch::Approx(7));
+ }
+
+ SECTION("order 8")
+ {
+ SimplicialGaussLegendreQuadrature<3> l7(7);
+
+ REQUIRE(l7.numberOfPoints() == 31);
+
+ REQUIRE(integrate(p0, l7) == Catch::Approx(16. / 3));
+ REQUIRE(integrate(p1, l7) == Catch::Approx(-16. / 3));
+ REQUIRE(integrate(p2, l7) == Catch::Approx(-47. / 3));
+ REQUIRE(integrate(p3, l7) == Catch::Approx(557. / 15));
+ REQUIRE(integrate(p4, l7) == Catch::Approx(4499. / 1575));
+ REQUIRE(integrate(p5, l7) == Catch::Approx(143443. / 7875));
+ REQUIRE(integrate(p6, l7) == Catch::Approx(-75773. / 2581974));
+ REQUIRE(integrate(p7, l7) == Catch::Approx(86951548. / 32274675));
+
+ REQUIRE(get_order(p8, l7, -863556317. / 322746750) == Catch::Approx(8).margin(0.01));
+ }
+
+ SECTION("not implemented formulae")
+ {
+ REQUIRE_THROWS_WITH(SimplicialGaussLegendreQuadrature<3>(8),
+ "error: Gauss-Legendre quadrature formulae handle orders up to 7 on tetrahedra");
+ }
+ }
}