diff --git a/src/analysis/EconomicalGaussQuadrature.cpp b/src/analysis/EconomicalGaussQuadrature.cpp
index efc83e10070ea381662c8a69fc46eacebb946cba..50385d1787fcd0adaed988e9efd2ef0a69732b40 100644
--- a/src/analysis/EconomicalGaussQuadrature.cpp
+++ b/src/analysis/EconomicalGaussQuadrature.cpp
@@ -313,212 +313,235 @@ template <>
void
EconomicalGaussQuadrature<3>::_buildPointAndWeightLists()
{
- using R3 = TinyVector<3>;
- using QuadratureGroup = std::array<double, 4>;
+ using R3 = TinyVector<3>;
- auto fill_quadrature_points = [](auto quadrature_group_list, auto& point_list, auto& weight_list) {
+ struct Descriptor
+ {
+ int id;
+ double weight;
+ std::vector<double> lambda_list;
+ };
+
+ auto fill_quadrature_points = [](auto descriptor_list, auto& point_list, auto& weight_list) {
Assert(point_list.size() == weight_list.size());
size_t k = 0;
- for (size_t i = 0; i < quadrature_group_list.size(); ++i) {
- const auto& [w, xi, eta, zeta] = quadrature_group_list[i];
- if (xi == eta) {
- if (xi == zeta) {
- if (xi == 0) {
- // xi = eta =zeta = 0
- point_list[k] = {0, 0, 0};
- weight_list[k] = w;
-
- k += 1;
- } else {
- // xi = eta = zeta != 0
- point_list[k + 0] = {+xi, +eta, +zeta};
- point_list[k + 1] = {-xi, +eta, +zeta};
- point_list[k + 2] = {+xi, -eta, +zeta};
- point_list[k + 3] = {-xi, -eta, +zeta};
- point_list[k + 4] = {+xi, +eta, -zeta};
- point_list[k + 5] = {-xi, +eta, -zeta};
- point_list[k + 6] = {+xi, -eta, -zeta};
- point_list[k + 7] = {-xi, -eta, -zeta};
-
- for (size_t i = 0; i < 8; ++i) {
- weight_list[k + i] = w;
- }
-
- k += 8;
- }
- } else {
- if (zeta == 0) {
- // xi = eta and zeta = 0
- point_list[k + 0] = {+xi, +eta, 0};
- point_list[k + 1] = {-xi, +eta, 0};
- point_list[k + 2] = {+xi, -eta, 0};
- point_list[k + 3] = {-xi, -eta, 0};
- point_list[k + 4] = {+xi, 0, +eta};
- point_list[k + 5] = {-xi, 0, +eta};
- point_list[k + 6] = {+xi, 0, -eta};
- point_list[k + 7] = {-xi, 0, -eta};
- point_list[k + 8] = {0, +xi, +eta};
- point_list[k + 9] = {0, -xi, +eta};
- point_list[k + 10] = {0, +xi, -eta};
- point_list[k + 11] = {0, -xi, -eta};
-
- for (size_t i = 0; i < 12; ++i) {
- weight_list[k + i] = w;
- }
-
- k += 12;
- } else {
- // xi = eta != zeta (!=0)
- point_list[k + 0] = {+xi, +xi, +zeta};
- point_list[k + 1] = {-xi, +xi, +zeta};
- point_list[k + 2] = {+xi, -xi, +zeta};
- point_list[k + 3] = {-xi, -xi, +zeta};
- point_list[k + 4] = {+xi, +xi, -zeta};
- point_list[k + 5] = {-xi, +xi, -zeta};
- point_list[k + 6] = {+xi, -xi, -zeta};
- point_list[k + 7] = {-xi, -xi, -zeta};
- point_list[k + 8] = {+xi, +zeta, +xi};
- point_list[k + 9] = {-xi, +zeta, +xi};
- point_list[k + 10] = {+xi, +zeta, -xi};
- point_list[k + 11] = {-xi, +zeta, -xi};
- point_list[k + 12] = {+xi, -zeta, +xi};
- point_list[k + 13] = {-xi, -zeta, +xi};
- point_list[k + 14] = {+xi, -zeta, -xi};
- point_list[k + 15] = {-xi, -zeta, -xi};
- point_list[k + 16] = {+zeta, +xi, +xi};
- point_list[k + 17] = {+zeta, -xi, +xi};
- point_list[k + 18] = {+zeta, +xi, -xi};
- point_list[k + 19] = {+zeta, -xi, -xi};
- point_list[k + 20] = {-zeta, +xi, +xi};
- point_list[k + 21] = {-zeta, -xi, +xi};
- point_list[k + 22] = {-zeta, +xi, -xi};
- point_list[k + 23] = {-zeta, -xi, -xi};
-
- for (size_t i = 0; i < 24; ++i) {
- weight_list[k + i] = w;
- }
-
- k += 24;
- }
+ for (size_t i = 0; i < descriptor_list.size(); ++i) {
+ const auto [id, unit_weight, position_list] = descriptor_list[i];
+
+ const double w = 8. * unit_weight;
+
+ switch (id) {
+ case 1: {
+ Assert(position_list.size() == 0);
+ point_list[k] = {0, 0, 0};
+ weight_list[k] = w;
+ ++k;
+ break;
+ }
+ case 2: {
+ Assert(position_list.size() == 1);
+ const double a = position_list[0];
+
+ point_list[k + 0] = {+a, 0, 0};
+ point_list[k + 1] = {-a, 0, 0};
+ point_list[k + 2] = {0, +a, 0};
+ point_list[k + 3] = {0, -a, 0};
+ point_list[k + 4] = {0, 0, +a};
+ point_list[k + 5] = {0, 0, -a};
+
+ for (size_t l = 0; l < 6; ++l) {
+ weight_list[k + l] = w;
}
- } else {
- if (zeta == 0) {
- if (eta == 0) {
- // (xi != 0) and eta = zeta = 0
- point_list[k + 0] = {+xi, 0, 0};
- point_list[k + 1] = {-xi, 0, 0};
- point_list[k + 2] = {0, +xi, 0};
- point_list[k + 3] = {0, -xi, 0};
- point_list[k + 4] = {0, 0, +xi};
- point_list[k + 5] = {0, 0, -xi};
-
- for (size_t i = 0; i < 6; ++i) {
- weight_list[k + i] = w;
- }
-
- k += 6;
-
- } else {
- // (xi != eta) and zeta = 0
- point_list[k + 0] = {+xi, +eta, 0};
- point_list[k + 1] = {-xi, +eta, 0};
- point_list[k + 2] = {+xi, -eta, 0};
- point_list[k + 3] = {-xi, -eta, 0};
- point_list[k + 4] = {+xi, 0, +eta};
- point_list[k + 5] = {-xi, 0, +eta};
- point_list[k + 6] = {+xi, 0, -eta};
- point_list[k + 7] = {-xi, 0, -eta};
- point_list[k + 8] = {0, +xi, +eta};
- point_list[k + 9] = {0, -xi, +eta};
- point_list[k + 10] = {0, +xi, -eta};
- point_list[k + 11] = {0, -xi, -eta};
-
- point_list[k + 12] = {+eta, +xi, 0};
- point_list[k + 13] = {-eta, +xi, 0};
- point_list[k + 14] = {+eta, -xi, 0};
- point_list[k + 15] = {-eta, -xi, 0};
- point_list[k + 16] = {+eta, 0, +xi};
- point_list[k + 17] = {-eta, 0, +xi};
- point_list[k + 18] = {+eta, 0, -xi};
- point_list[k + 19] = {-eta, 0, -xi};
- point_list[k + 20] = {0, +eta, +xi};
- point_list[k + 21] = {0, -eta, +xi};
- point_list[k + 22] = {0, +eta, -xi};
- point_list[k + 23] = {0, -eta, -xi};
-
- for (size_t i = 0; i < 24; ++i) {
- weight_list[k + i] = w;
- }
-
- k += 24;
- }
- } else {
- // xi != eta != zeta
- point_list[k + 0] = {+xi, +eta, +zeta};
- point_list[k + 1] = {-xi, +eta, +zeta};
- point_list[k + 2] = {+xi, -eta, +zeta};
- point_list[k + 3] = {-xi, -eta, +zeta};
- point_list[k + 4] = {+xi, +eta, -zeta};
- point_list[k + 5] = {-xi, +eta, -zeta};
- point_list[k + 6] = {+xi, -eta, -zeta};
- point_list[k + 7] = {-xi, -eta, -zeta};
-
- point_list[k + 8] = {+xi, +zeta, +eta};
- point_list[k + 9] = {-xi, +zeta, +eta};
- point_list[k + 10] = {+xi, +zeta, -eta};
- point_list[k + 11] = {-xi, +zeta, -eta};
- point_list[k + 12] = {+xi, -zeta, +eta};
- point_list[k + 13] = {-xi, -zeta, +eta};
- point_list[k + 14] = {+xi, -zeta, -eta};
- point_list[k + 15] = {-xi, -zeta, -eta};
-
- point_list[k + 16] = {+eta, +xi, +zeta};
- point_list[k + 17] = {+eta, -xi, +zeta};
- point_list[k + 18] = {-eta, +xi, +zeta};
- point_list[k + 19] = {-eta, -xi, +zeta};
- point_list[k + 20] = {+eta, +xi, -zeta};
- point_list[k + 21] = {+eta, -xi, -zeta};
- point_list[k + 22] = {-eta, +xi, -zeta};
- point_list[k + 23] = {-eta, -xi, -zeta};
-
- point_list[k + 24] = {+eta, +zeta, +xi};
- point_list[k + 25] = {+eta, +zeta, -xi};
- point_list[k + 26] = {-eta, +zeta, +xi};
- point_list[k + 27] = {-eta, +zeta, -xi};
- point_list[k + 28] = {+eta, -zeta, +xi};
- point_list[k + 29] = {+eta, -zeta, -xi};
- point_list[k + 30] = {-eta, -zeta, +xi};
- point_list[k + 31] = {-eta, -zeta, -xi};
-
- point_list[k + 32] = {+zeta, +xi, +eta};
- point_list[k + 33] = {+zeta, -xi, +eta};
- point_list[k + 34] = {+zeta, +xi, -eta};
- point_list[k + 35] = {+zeta, -xi, -eta};
- point_list[k + 36] = {-zeta, +xi, +eta};
- point_list[k + 37] = {-zeta, -xi, +eta};
- point_list[k + 38] = {-zeta, +xi, -eta};
- point_list[k + 39] = {-zeta, -xi, -eta};
-
- point_list[k + 40] = {+zeta, +eta, +xi};
- point_list[k + 41] = {+zeta, +eta, -xi};
- point_list[k + 42] = {+zeta, -eta, +xi};
- point_list[k + 43] = {+zeta, -eta, -xi};
- point_list[k + 44] = {-zeta, +eta, +xi};
- point_list[k + 45] = {-zeta, +eta, -xi};
- point_list[k + 46] = {-zeta, -eta, +xi};
- point_list[k + 47] = {-zeta, -eta, -xi};
-
- for (size_t i = 0; i < 48; ++i) {
- weight_list[k + i] = w;
- }
- k += 48;
+ k += 6;
+ break;
+ }
+ case 3: {
+ Assert(position_list.size() == 1);
+ const double a = position_list[0];
+
+ point_list[k + 0] = {+a, +a, +a};
+ point_list[k + 1] = {+a, +a, -a};
+ point_list[k + 2] = {+a, -a, +a};
+ point_list[k + 3] = {+a, -a, -a};
+ point_list[k + 4] = {-a, +a, +a};
+ point_list[k + 5] = {-a, +a, -a};
+ point_list[k + 6] = {-a, -a, +a};
+ point_list[k + 7] = {-a, -a, -a};
+
+ for (size_t l = 0; l < 8; ++l) {
+ weight_list[k + l] = w;
}
+
+ k += 8;
+ break;
}
- }
+ case 4: {
+ Assert(position_list.size() == 1);
+ const double a = position_list[0];
+
+ point_list[k + 0] = {+a, +a, 0};
+ point_list[k + 1] = {+a, -a, 0};
+ point_list[k + 2] = {-a, +a, 0};
+ point_list[k + 3] = {-a, -a, 0};
+ point_list[k + 4] = {+a, 0, +a};
+ point_list[k + 5] = {+a, 0, -a};
+ point_list[k + 6] = {-a, 0, +a};
+ point_list[k + 7] = {-a, 0, -a};
+ point_list[k + 8] = {0, +a, +a};
+ point_list[k + 9] = {0, +a, -a};
+ point_list[k + 10] = {0, -a, +a};
+ point_list[k + 11] = {0, -a, -a};
+
+ for (size_t l = 0; l < 12; ++l) {
+ weight_list[k + l] = w;
+ }
- Assert(k == point_list.size(), "invalid number of quadrature points");
+ k += 12;
+ break;
+ }
+ case 5: {
+ Assert(position_list.size() == 2);
+ const double a = position_list[0];
+ const double b = position_list[1];
+
+ point_list[k + 0] = {+a, +b, 0};
+ point_list[k + 1] = {+a, -b, 0};
+ point_list[k + 2] = {-a, +b, 0};
+ point_list[k + 3] = {-a, -b, 0};
+ point_list[k + 4] = {+b, +a, 0};
+ point_list[k + 5] = {-b, +a, 0};
+ point_list[k + 6] = {+b, -a, 0};
+ point_list[k + 7] = {-b, -a, 0};
+ point_list[k + 8] = {+a, 0, +b};
+ point_list[k + 9] = {+a, 0, -b};
+ point_list[k + 10] = {-a, 0, +b};
+ point_list[k + 11] = {-a, 0, -b};
+ point_list[k + 12] = {+b, 0, +a};
+ point_list[k + 13] = {-b, 0, +a};
+ point_list[k + 14] = {+b, 0, -a};
+ point_list[k + 15] = {-b, 0, -a};
+ point_list[k + 16] = {0, +a, +b};
+ point_list[k + 17] = {0, +a, -b};
+ point_list[k + 18] = {0, -a, +b};
+ point_list[k + 19] = {0, -a, -b};
+ point_list[k + 20] = {0, +b, +a};
+ point_list[k + 21] = {0, -b, +a};
+ point_list[k + 22] = {0, +b, -a};
+ point_list[k + 23] = {0, -b, -a};
+
+ for (size_t l = 0; l < 24; ++l) {
+ weight_list[k + l] = w;
+ }
+
+ k += 24;
+ break;
+ }
+ case 6: {
+ Assert(position_list.size() == 2);
+ const double a = position_list[0];
+ const double b = position_list[1];
+
+ point_list[k + 0] = {+a, +a, +b};
+ point_list[k + 1] = {+a, +a, -b};
+ point_list[k + 2] = {+a, -a, +b};
+ point_list[k + 3] = {+a, -a, -b};
+ point_list[k + 4] = {-a, +a, +b};
+ point_list[k + 5] = {-a, +a, -b};
+ point_list[k + 6] = {-a, -a, +b};
+ point_list[k + 7] = {-a, -a, -b};
+ point_list[k + 8] = {+a, +b, +a};
+ point_list[k + 9] = {+a, -b, +a};
+ point_list[k + 10] = {+a, +b, -a};
+ point_list[k + 11] = {+a, -b, -a};
+ point_list[k + 12] = {-a, +b, +a};
+ point_list[k + 13] = {-a, -b, +a};
+ point_list[k + 14] = {-a, +b, -a};
+ point_list[k + 15] = {-a, -b, -a};
+ point_list[k + 16] = {+b, +a, +a};
+ point_list[k + 17] = {-b, +a, +a};
+ point_list[k + 18] = {+b, +a, -a};
+ point_list[k + 19] = {-b, +a, -a};
+ point_list[k + 20] = {+b, -a, +a};
+ point_list[k + 21] = {-b, -a, +a};
+ point_list[k + 22] = {+b, -a, -a};
+ point_list[k + 23] = {-b, -a, -a};
+
+ for (size_t l = 0; l < 24; ++l) {
+ weight_list[k + l] = w;
+ }
+
+ k += 24;
+ break;
+ }
+ case 7: {
+ Assert(position_list.size() == 3);
+ const double a = position_list[0];
+ const double b = position_list[1];
+ const double c = position_list[2];
+
+ point_list[k + 0] = {+a, +b, +c};
+ point_list[k + 1] = {+a, +b, -c};
+ point_list[k + 2] = {+a, -b, +c};
+ point_list[k + 3] = {+a, -b, -c};
+ point_list[k + 4] = {-a, +b, +c};
+ point_list[k + 5] = {-a, +b, -c};
+ point_list[k + 6] = {-a, -b, +c};
+ point_list[k + 7] = {-a, -b, -c};
+ point_list[k + 8] = {+a, +c, +b};
+ point_list[k + 9] = {+a, -c, +b};
+ point_list[k + 10] = {+a, +c, -b};
+ point_list[k + 11] = {+a, -c, -b};
+ point_list[k + 12] = {-a, +c, +b};
+ point_list[k + 13] = {-a, -c, +b};
+ point_list[k + 14] = {-a, +c, -b};
+ point_list[k + 15] = {-a, -c, -b};
+ point_list[k + 16] = {+b, +a, +c};
+ point_list[k + 17] = {+b, +a, -c};
+ point_list[k + 18] = {-b, +a, +c};
+ point_list[k + 19] = {-b, +a, -c};
+ point_list[k + 20] = {+b, -a, +c};
+ point_list[k + 21] = {+b, -a, -c};
+ point_list[k + 22] = {-b, -a, +c};
+ point_list[k + 23] = {-b, -a, -c};
+ point_list[k + 24] = {+b, +c, +a};
+ point_list[k + 25] = {+b, -c, +a};
+ point_list[k + 26] = {-b, +c, +a};
+ point_list[k + 27] = {-b, -c, +a};
+ point_list[k + 28] = {+b, +c, -a};
+ point_list[k + 29] = {+b, -c, -a};
+ point_list[k + 30] = {-b, +c, -a};
+ point_list[k + 31] = {-b, -c, -a};
+ point_list[k + 32] = {+c, +a, +b};
+ point_list[k + 33] = {-c, +a, +b};
+ point_list[k + 34] = {+c, +a, -b};
+ point_list[k + 35] = {-c, +a, -b};
+ point_list[k + 36] = {+c, -a, +b};
+ point_list[k + 37] = {-c, -a, +b};
+ point_list[k + 38] = {+c, -a, -b};
+ point_list[k + 39] = {-c, -a, -b};
+ point_list[k + 40] = {+c, +b, +a};
+ point_list[k + 41] = {-c, +b, +a};
+ point_list[k + 42] = {+c, -b, +a};
+ point_list[k + 43] = {-c, -b, +a};
+ point_list[k + 44] = {+c, +b, -a};
+ point_list[k + 45] = {-c, +b, -a};
+ point_list[k + 46] = {+c, -b, -a};
+ point_list[k + 47] = {-c, -b, -a};
+
+ for (size_t l = 0; l < 48; ++l) {
+ weight_list[k + l] = w;
+ }
+
+ k += 48;
+ break;
+ }
+ default: {
+ throw UnexpectedError("invalid quadrature id");
+ }
+ }
+ }
};
switch (m_degree) {
@@ -528,10 +551,10 @@ EconomicalGaussQuadrature<3>::_buildPointAndWeightLists()
SmallArray<R3> point_list(nb_points);
SmallArray<double> weight_list(nb_points);
- std::array<QuadratureGroup, 1> quadrature_group_list = //
- {QuadratureGroup{8.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}};
+ std::array descriptor_list = //
+ {Descriptor{1, 1.000000000000000e+00, {}}};
- fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+ fill_quadrature_points(descriptor_list, point_list, weight_list);
m_point_list = point_list;
m_weight_list = weight_list;
@@ -543,10 +566,10 @@ EconomicalGaussQuadrature<3>::_buildPointAndWeightLists()
SmallArray<R3> point_list(nb_points);
SmallArray<double> weight_list(nb_points);
- std::array<QuadratureGroup, 1> quadrature_group_list = //
- {QuadratureGroup{1.33333333333333, 1.000000000000000, 0.00000000000000, 0.00000000000000}};
+ std::array descriptor_list = //
+ {Descriptor{2, 1.666666666666667e-01, {+1.000000000000000e+00}}};
- fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+ fill_quadrature_points(descriptor_list, point_list, weight_list);
m_point_list = point_list;
m_weight_list = weight_list;
@@ -558,11 +581,11 @@ EconomicalGaussQuadrature<3>::_buildPointAndWeightLists()
SmallArray<R3> point_list(nb_points);
SmallArray<double> weight_list(nb_points);
- std::array quadrature_group_list = //
- {QuadratureGroup{8.86426592797784E-01, 0.795822425754221, 0.000000000000000, 0.000000000000000},
- QuadratureGroup{3.35180055401662E-01, 0.758786910639328, 0.758786910639328, 0.758786910639328}};
+ std::array descriptor_list = //
+ {Descriptor{2, 1.108033240997230e-01, {-7.958224257542215e-01}},
+ Descriptor{3, 4.189750692520776e-02, {+7.587869106393281e-01}}};
- fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+ fill_quadrature_points(descriptor_list, point_list, weight_list);
m_point_list = point_list;
m_weight_list = weight_list;
@@ -570,17 +593,17 @@ EconomicalGaussQuadrature<3>::_buildPointAndWeightLists()
}
case 6:
case 7: {
- constexpr size_t nb_points = 27;
+ constexpr size_t nb_points = 34;
SmallArray<R3> point_list(nb_points);
SmallArray<double> weight_list(nb_points);
- std::array quadrature_group_list = //
- {QuadratureGroup{7.88073482744211E-01, 0.000000000000000, 0.000000000000000, 0.000000000000000},
- QuadratureGroup{4.99369002307720E-01, 0.848418011472252, 0.000000000000000, 0.000000000000000},
- QuadratureGroup{4.78508449425127E-01, 0.652816472101691, 0.652816472101691, 0.652816472101691},
- QuadratureGroup{3.23037423340374E-02, 1.106412898626718, 1.106412898626718, 0.000000000000000}};
+ std::array descriptor_list = //
+ {Descriptor{2, 3.400458474980031e-02, {-9.388053721060176e-01}},
+ Descriptor{3, 2.529672440135186e-02, {+7.463336100128160e-01}},
+ Descriptor{3, 5.417063533485642e-02, {+4.101308983320144e-01}},
+ Descriptor{4, 1.335280113429432e-02, {+9.064606901228371e-01}}};
- fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+ fill_quadrature_points(descriptor_list, point_list, weight_list);
m_point_list = point_list;
m_weight_list = weight_list;
@@ -588,19 +611,18 @@ EconomicalGaussQuadrature<3>::_buildPointAndWeightLists()
}
case 8:
case 9: {
- constexpr size_t nb_points = 53;
+ constexpr size_t nb_points = 58;
SmallArray<R3> point_list(nb_points);
SmallArray<double> weight_list(nb_points);
- std::array quadrature_group_list = //
- {QuadratureGroup{+5.88405321380412e-01, 0.000000000000000, 0.000000000000000, 0.000000000000000},
- QuadratureGroup{-1.52097068487023e-01, 1.064082230328777, 0.000000000000000, 0.000000000000000},
- QuadratureGroup{+3.69012523996709e-01, 0.905830033000216, 0.000000000000000, 0.000000000000000},
- QuadratureGroup{+1.04007450974435e-01, 0.817286490798906, 0.817286490798906, 0.817286490798906},
- QuadratureGroup{+3.80660357224238e-01, 0.501292956337400, 0.501292956337400, 0.501292956337400},
- QuadratureGroup{+9.30316449988371e-02, 1.017168937265364, 0.650007853956632, 0.000000000000000}};
+ std::array descriptor_list = //
+ {Descriptor{2, 5.415937446870682e-02, {-6.136814695917090e-01}},
+ Descriptor{3, 6.268599412418628e-03, {+8.700997846619759e-01}},
+ Descriptor{3, 2.485747976800294e-02, {+5.641108070200300e-01}},
+ Descriptor{4, 1.147372576702221e-02, {+8.776871232576783e-01}},
+ Descriptor{6, 1.201460043917167e-02, {+4.322679026308622e-01, +9.385304218646717e-01}}};
- fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+ fill_quadrature_points(descriptor_list, point_list, weight_list);
m_point_list = point_list;
m_weight_list = weight_list;
@@ -608,22 +630,20 @@ EconomicalGaussQuadrature<3>::_buildPointAndWeightLists()
}
case 10:
case 11: {
- constexpr size_t nb_points = 89;
+ constexpr size_t nb_points = 90;
SmallArray<R3> point_list(nb_points);
SmallArray<double> weight_list(nb_points);
- std::array quadrature_group_list = //
- {QuadratureGroup{-7.29841346362385e-01, 0.000000000000000, 0.000000000000000, 0.000000000000000},
- QuadratureGroup{-2.95542895787525e-01, 1.044892251909493, 0.000000000000000, 0.000000000000000},
- QuadratureGroup{+7.52817194936725e-02, 1.128826539781045, 0.000000000000000, 0.000000000000000},
- QuadratureGroup{+2.88865143683389e-01, 0.965386862750349, 0.000000000000000, 0.000000000000000},
- QuadratureGroup{+3.93648801431496e-01, 0.398089853865140, 0.000000000000000, 0.000000000000000},
- QuadratureGroup{+1.41560420881811e-02, 0.888733872960214, 0.888733872960214, 0.888733872960214},
- QuadratureGroup{+2.32318705969944e-01, 0.554161513623834, 0.554161513623834, 0.554161513623834},
- QuadratureGroup{+1.18935816908526e-01, 0.920049645350334, 0.544394886959702, 0.000000000000000},
- QuadratureGroup{+4.70861309652733e-02, 0.904534033733291, 0.904534033733291, 0.496742636335202}};
+ std::array descriptor_list = //
+ {Descriptor{2, 2.961692452402056e-02, {+7.221330388744185e-01}},
+ Descriptor{3, 7.751183787523073e-03, {+8.094882019630989e-01}},
+ Descriptor{3, 2.222100810905048e-02, {-5.336540088804971e-01}},
+ Descriptor{3, 1.698870214455200e-02, {+2.807725866512744e-01}},
+ Descriptor{4, 1.602728177693560e-02, {+8.039334672152845e-01}},
+ Descriptor{6, 1.618821190032323e-03, {+9.800994910090715e-01, -5.307838311938264e-01}},
+ Descriptor{6, 8.976342110119554e-03, {+4.056859801950964e-01, +9.545832189295661e-01}}};
- fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+ fill_quadrature_points(descriptor_list, point_list, weight_list);
m_point_list = point_list;
m_weight_list = weight_list;
@@ -631,24 +651,22 @@ EconomicalGaussQuadrature<3>::_buildPointAndWeightLists()
}
case 12:
case 13: {
- constexpr size_t nb_points = 151;
+ constexpr size_t nb_points = 154;
SmallArray<R3> point_list(nb_points);
SmallArray<double> weight_list(nb_points);
- std::array quadrature_group_list = //
- {QuadratureGroup{+1.20671303848159e-01, 0.000000000000000, 0.000000000000000, 0.000000000000000},
- QuadratureGroup{+8.27155895731376e-02, 0.957697664927095, 0.000000000000000, 0.000000000000000},
- QuadratureGroup{-2.00338143812084e-01, 0.799168794176284, 0.000000000000000, 0.000000000000000},
- QuadratureGroup{+1.19014821875897e-01, 0.478190772481903, 0.000000000000000, 0.000000000000000},
- QuadratureGroup{+3.36821300476000e-02, 0.863636428036401, 0.863636428036401, 0.863636428036401},
- QuadratureGroup{+9.91539279234894e-02, 0.625548513457035, 0.625548513457035, 0.625548513457035},
- QuadratureGroup{+8.83375530928565e-02, 0.367859232104514, 0.367859232104514, 0.367859232104514},
- QuadratureGroup{+3.89329599510302e-02, 0.869451316547931, 0.869451316547931, 0.000000000000000},
- QuadratureGroup{+1.54603530012075e-01, 0.767220340241960, 0.354778664757558, 0.000000000000000},
- QuadratureGroup{+3.80152555202625e-03, 1.012782207207705, 1.012782207207705, 0.537923170210472},
- QuadratureGroup{+3.81806114347457e-02, 0.944392525641525, 0.688378084301308, 0.372427809150424}};
+ std::array descriptor_list = //
+ {Descriptor{2, 1.739128248747042e-02, {-6.344099893707986e-01}},
+ Descriptor{3, 1.838657780061070e-03, {+9.078540479628353e-01}},
+ Descriptor{3, 1.129102483142921e-02, {-2.362052584516462e-01}},
+ Descriptor{4, 5.975759635289482e-03, {+7.375746343132366e-01}},
+ Descriptor{5, 8.691231524417118e-03, {+4.607161476964481e-01, +9.333452030108366e-01}},
+ Descriptor{6, 1.061455871034256e-02, {+3.733841515799997e-01, +6.642163895931538e-01}},
+ Descriptor{6, 1.690358916666084e-03, {+9.495110814642814e-01, +2.576763878569212e-01}},
+ Descriptor{6, 2.284240552839863e-03, {+6.656694827601780e-01, +9.948700428018972e-01}},
+ Descriptor{6, 6.674015652391932e-03, {-8.051009105032320e-01, -4.521405059548300e-01}}};
- fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+ fill_quadrature_points(descriptor_list, point_list, weight_list);
m_point_list = point_list;
m_weight_list = weight_list;
@@ -656,27 +674,26 @@ EconomicalGaussQuadrature<3>::_buildPointAndWeightLists()
}
case 14:
case 15: {
- constexpr size_t nb_points = 235;
+ constexpr size_t nb_points = 256;
SmallArray<R3> point_list(nb_points);
SmallArray<double> weight_list(nb_points);
- std::array quadrature_group_list = //
- {QuadratureGroup{+3.47325185603722e-02, 0.000000000000000, 0.000000000000000, 0.000000000000000},
- QuadratureGroup{+6.06208354441015e-02, 0.955412599148645, 0.000000000000000, 0.000000000000000},
- QuadratureGroup{-9.26612667312444e-02, 0.851429576022840, 0.000000000000000, 0.000000000000000},
- QuadratureGroup{+7.67301289360290e-02, 0.339914995753808, 0.000000000000000, 0.000000000000000},
- QuadratureGroup{+3.54673063527669e-05, 1.122550051350066, 1.122550051350066, 1.122550051350066},
- QuadratureGroup{+6.13391727954448e-02, 0.683321227375114, 0.683321227375114, 0.683321227375114},
- QuadratureGroup{+1.23406122314294e-01, 0.398466735636029, 0.398466735636029, 0.398466735636029},
- QuadratureGroup{+5.04204099161268e-03, 1.031795356557355, 0.870614907692603, 0.000000000000000},
- QuadratureGroup{+7.27045726671404e-02, 0.859754056650996, 0.417236245841784, 0.000000000000000},
- QuadratureGroup{+5.55784362060740e-02, 0.677373166268794, 0.234354064403370, 0.000000000000000},
- QuadratureGroup{+1.51521479476415e-02, 0.931914101559387, 0.931914101559387, 0.710695749151767},
- QuadratureGroup{+2.70620502572634e-02, 0.871331764342344, 0.871331764342344, 0.277186077584005},
- QuadratureGroup{+4.18922675266216e-02, 0.693611775581488, 0.693611775581488, 0.258751904925954},
- QuadratureGroup{+2.08443087896894e-02, 0.959113457650743, 0.622683549553047, 0.406412429958087}};
-
- fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+ std::array descriptor_list = //
+ {Descriptor{2, 3.340305010694517e-03, {+9.239155149504112e-01}},
+ Descriptor{2, 5.018552337877279e-03, {-2.173948613666498e-01}},
+ Descriptor{3, 9.082758817224784e-03, {+3.384731869494254e-01}},
+ Descriptor{3, 2.936408738387597e-03, {+7.930891310393379e-01}},
+ Descriptor{4, 9.318238647048784e-03, {+6.965519027266717e-01}},
+ Descriptor{5, 6.990786235751096e-03, {+2.642667592628472e-01, +5.849899486873243e-01}},
+ Descriptor{5, 2.841005747490632e-03, {+5.529505913370675e-01, +9.828048727801773e-01}},
+ Descriptor{6, 5.309574072330840e-03, {+2.720325890860906e-01, +8.780894076733160e-01}},
+ Descriptor{6, 3.570894536567292e-03, {+6.186478600855522e-01, +9.399813353075436e-01}},
+ Descriptor{6, 3.505872179585672e-03, {+8.631544625530380e-01, +3.042225951752605e-01}},
+ Descriptor{6, 6.665199051138311e-03, {+4.616057726468051e-01, +7.371361183496182e-01}},
+ Descriptor{6, 7.781518386120085e-04, {+9.603141805416257e-01, +7.670146080341008e-01}},
+ Descriptor{7, 6.249800796596737e-04, {+2.802518397917547e-01, +8.885938329579797e-01, +9.942601945848643e-01}}};
+
+ fill_quadrature_points(descriptor_list, point_list, weight_list);
m_point_list = point_list;
m_weight_list = weight_list;
@@ -684,30 +701,29 @@ EconomicalGaussQuadrature<3>::_buildPointAndWeightLists()
}
case 16:
case 17: {
- constexpr size_t nb_points = 307;
+ constexpr size_t nb_points = 346;
SmallArray<R3> point_list(nb_points);
SmallArray<double> weight_list(nb_points);
- std::array quadrature_group_list = //
- {QuadratureGroup{+1.27375023441727e-01, 0.000000000000000, 0.000000000000000, 0.000000000000000},
- QuadratureGroup{-6.41337796713059e-03, 1.086788874756602, 0.000000000000000, 0.000000000000000},
- QuadratureGroup{+5.09259994505536e-02, 0.922008663747819, 0.000000000000000, 0.000000000000000},
- QuadratureGroup{-1.13461087562205e-01, 0.641870173734099, 0.000000000000000, 0.000000000000000},
- QuadratureGroup{+1.11151619355886e-02, 0.891607017924525, 0.891607017924525, 0.891607017924525},
- QuadratureGroup{+3.71864272269592e-02, 0.734488336558120, 0.734488336558120, 0.734488336558120},
- QuadratureGroup{+7.72186740324679e-02, 0.518312376398533, 0.518312376398533, 0.518312376398533},
- QuadratureGroup{+2.13441739735667e-03, 0.989312892186012, 0.989312892186012, 0.000000000000000},
- QuadratureGroup{+8.08047120120055e-02, 0.364474560584746, 0.364474560584746, 0.000000000000000},
- QuadratureGroup{+1.59035001385071e-02, 0.948247130849454, 0.780127958605308, 0.000000000000000},
- QuadratureGroup{+7.65449188677983e-03, 1.045702906353664, 0.348818446769418, 0.000000000000000},
- QuadratureGroup{+8.55284062610809e-02, 0.660531051473545, 0.241437081938889, 0.000000000000000},
- QuadratureGroup{+1.97991253834944e-03, 0.990502583320736, 0.990502583320736, 0.757226655088057},
- QuadratureGroup{+8.14314291841830e-03, 0.926857250551066, 0.926857250551066, 0.315061197287156},
- QuadratureGroup{+3.57632627556923e-02, 0.710700783405400, 0.710700783405400, 0.208469339137832},
- QuadratureGroup{+1.55500497623716e-02, 0.957873006911044, 0.780040197522315, 0.527111921367462},
- QuadratureGroup{+3.79403443748498e-02, 0.880972079148717, 0.513930736899189, 0.269433491346948}};
-
- fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+ std::array descriptor_list = //
+ {Descriptor{2, 9.614592128929721e-03, {+4.853150302302675e-01}},
+ Descriptor{3, 4.393827602159037e-03, {+7.406788158740985e-01}},
+ Descriptor{3, 4.964725194523337e-05, {+9.999995835822517e-01}},
+ Descriptor{4, 2.506284464832531e-03, {+7.546268052131762e-01}},
+ Descriptor{4, 6.133232482070283e-03, {+5.496818141686911e-01}},
+ Descriptor{4, 3.052513689971173e-04, {+9.996685037422811e-01}},
+ Descriptor{5, 2.825437733033075e-03, {+6.157290991684734e-01, +9.530419407049950e-01}},
+ Descriptor{5, 4.817704598799924e-03, {+2.756812731993550e-01, +8.171240681070916e-01}},
+ Descriptor{6, 2.283142338953617e-03, {+2.553429540372821e-01, +9.697740446619386e-01}},
+ Descriptor{6, 1.177437487278679e-03, {+9.352363848980443e-01, +7.467233074183380e-01}},
+ Descriptor{6, 6.184060750627916e-03, {+3.430591461388653e-01, +6.186639578085585e-01}},
+ Descriptor{6, 2.674655837593652e-03, {+2.555841601150721e-01, -1.540702494043778e-01}},
+ Descriptor{6, 3.419177677578799e-03, {+5.416211416468879e-01, +8.956569517854811e-01}},
+ Descriptor{6, 2.334617556615003e-03, {+8.874561091440486e-01, +2.414649778934788e-01}},
+ Descriptor{7, 1.052790164887146e-03, {+8.014733742030056e-01, +9.887847181449660e-01, +4.683593080344387e-01}},
+ Descriptor{7, 2.743830940763945e-03, {+6.099140660050306e-01, +7.868230999431114e-01, +3.205108203811766e-01}}};
+
+ fill_quadrature_points(descriptor_list, point_list, weight_list);
m_point_list = point_list;
m_weight_list = weight_list;
@@ -715,44 +731,77 @@ EconomicalGaussQuadrature<3>::_buildPointAndWeightLists()
}
case 18:
case 19: {
- constexpr size_t nb_points = 435;
+ constexpr size_t nb_points = 454;
SmallArray<R3> point_list(nb_points);
SmallArray<double> weight_list(nb_points);
- std::array quadrature_group_list = //
- {QuadratureGroup{-2.96256080038908e-01, 0.000000000000000, 0.000000000000000, 0.000000000000000},
- QuadratureGroup{-3.18213148997343e-02, 0.950006122416513, 0.000000000000000, 0.000000000000000},
- QuadratureGroup{+8.47491973828124e-02, 0.914965429121309, 0.000000000000000, 0.000000000000000},
- QuadratureGroup{-1.86007315013108e-01, 0.822592418605602, 0.000000000000000, 0.000000000000000},
- QuadratureGroup{+1.19234266342962e-01, 0.652114589201603, 0.000000000000000, 0.000000000000000},
- QuadratureGroup{+9.27104797995355e-02, 0.157242622070114, 0.000000000000000, 0.000000000000000},
- QuadratureGroup{+6.00225327364294e-02, 1.001394961101233, 1.001394961101233, 1.001394961101233},
- QuadratureGroup{+6.01997623450650e-03, 0.949635648752703, 0.949635648752703, 0.949635648752703},
- QuadratureGroup{+3.56895101325082e-02, 0.594117176975313, 0.594117176975313, 0.594117176975313},
- QuadratureGroup{+8.17900786633523e-02, 0.326702572971493, 0.326702572971493, 0.326702572971493},
- QuadratureGroup{-3.59448438162817e-02, 1.003460985202348, 1.003460985202348, 0.000000000000000},
- QuadratureGroup{-6.27840019030977e-02, 0.893211074627396, 0.893211074627396, 0.000000000000000},
- QuadratureGroup{+5.12421306349611e-02, 0.488224948121536, 0.488224948121536, 0.000000000000000},
- QuadratureGroup{+2.09987036073207e-02, 1.010406331416018, 0.988440552758072, 0.000000000000000},
- QuadratureGroup{+1.29343518268887e-02, 0.989413303740966, 0.384711778165505, 0.000000000000000},
- QuadratureGroup{+5.08521957924782e-02, 0.907323294169450, 0.833136115100468, 0.000000000000000},
- QuadratureGroup{+6.25369819391845e-02, 0.831834985558718, 0.220277012103345, 0.000000000000000},
- QuadratureGroup{-2.31758216427620e-02, 1.012434400912337, 1.012434400912337, 0.978322332670673},
- QuadratureGroup{+9.71367824650409e-03, 0.940302135044088, 0.940302135044088, 0.343148808079936},
- QuadratureGroup{+2.66634824440037e-03, 1.036512025425065, 1.036512025425065, 0.933501389902834},
- QuadratureGroup{+2.38739030828427e-02, 0.827104008787951, 0.827104008787951, 0.595273916358421},
- QuadratureGroup{+4.49207551990765e-03, 0.853891500333261, 0.990340675813395, 0.703401752648627},
- QuadratureGroup{+2.05647395630476e-02, 0.946965100974417, 0.623429573213959, 0.361644847592954},
- QuadratureGroup{+3.90081809778401e-02, 0.754378926224437, 0.538990774213829, 0.258107746839302}};
+ std::array descriptor_list = //
+ {Descriptor{2, 2.057371005731085e-03, {+2.357950388975374e-01}},
+ Descriptor{3, 9.660877915728369e-04, {+8.884950400440639e-01}},
+ Descriptor{3, 1.635901446799172e-03, {+6.691818469930052e-01}},
+ Descriptor{5, 8.499912813105116e-04, {+7.123100397018546e-01, +9.982977792850621e-01}},
+ Descriptor{5, 3.595003721736422e-03, {+3.997141162742279e-01, +1.720241294475643e-01}},
+ Descriptor{5, 3.836063429622649e-03, {+8.765805264721680e-01, -4.319090088935915e-01}},
+ Descriptor{6, 5.044128485798918e-03, {-3.101731974045544e-01, +5.411043641351523e-01}},
+ Descriptor{6, 1.951830532879285e-04, {+9.841852549757261e-01, +8.206649948892372e-01}},
+ Descriptor{6, 1.758454840661237e-03, {+8.443872335315854e-01, -1.473766326322086e-01}},
+ Descriptor{6, 3.887320262932200e-03, {+4.875993998852213e-01, +7.674914745033637e-01}},
+ Descriptor{6, 1.887401083675721e-03, {+7.641318497967444e-01, +3.863910365012072e-01}},
+ Descriptor{6, 1.699165183521457e-03, {+7.112537388224306e-01, +9.032012949929131e-01}},
+ Descriptor{6, 7.588378288377641e-04, {+9.551489495939087e-01, +2.226160644570115e-01}},
+ Descriptor{6, 1.455762954103218e-03, {+1.696919803869971e-01, +9.651628675592381e-01}},
+ Descriptor{6, 4.007905075757146e-03, {+6.085335258834057e-01, +9.412726803541549e-02}},
+ Descriptor{6, 1.399186075584291e-03, {+4.265517787687217e-01, +9.784781131096041e-01}},
+ Descriptor{6, 4.099280788926685e-03, {+1.836920702608942e-01, +7.457754226757105e-01}},
+ Descriptor{7, 1.951861257384991e-03, {+6.535330886596633e-01, +9.133523126050157e-01, +2.888629020786040e-01}},
+ Descriptor{7, 9.537937942918817e-04, {+8.475635090984786e-01, +9.748482139426561e-01, +5.534277540813458e-01}}};
+
+ fill_quadrature_points(descriptor_list, point_list, weight_list);
- fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 20:
+ case 21: {
+ constexpr size_t nb_points = 580;
+ SmallArray<R3> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ std::array descriptor_list = //
+ {Descriptor{2, 3.007917561622623e-03, {-6.174957125232848e-01}},
+ Descriptor{2, 6.157551855986320e-03, {-2.913675401437903e-01}},
+ Descriptor{3, 7.451864938264940e-04, {+8.622583182294242e-01}},
+ Descriptor{3, 3.473649875824585e-03, {+5.666798116095041e-01}},
+ Descriptor{5, 5.231065457280251e-04, {+7.962288576876380e-01, +9.999999871114571e-01}},
+ Descriptor{5, 4.721261422431887e-03, {+6.935661772320008e-01, +3.463806774517028e-01}},
+ Descriptor{5, 2.337222373763527e-03, {+8.718362020420203e-01, -6.252215000307171e-01}},
+ Descriptor{6, 3.841754898565957e-03, {-4.164545409594995e-01, +8.053481425029153e-01}},
+ Descriptor{6, 1.507767719634259e-04, {+9.738519593091541e-01, +8.745003988374457e-01}},
+ Descriptor{6, 3.214079409518292e-03, {+7.290234505608810e-01, +2.451675433295293e-01}},
+ Descriptor{6, 2.245446622613030e-03, {+6.496117102813809e-01, +8.441031264282526e-01}},
+ Descriptor{6, 1.041367673151780e-03, {+8.919707516952021e-01, +1.301935885566839e-01}},
+ Descriptor{6, 5.920000965024342e-04, {+7.833271437290619e-01, +9.532187072534166e-01}},
+ Descriptor{6, 4.100172337698149e-04, {+9.699411873288849e-01, +2.710159328732420e-01}},
+ Descriptor{6, 7.465850926621290e-04, {+1.720339772207737e-01, +9.847690037797490e-01}},
+ Descriptor{6, 3.743289683970358e-03, {+5.380985926191851e-01, -2.207631346215846e-01}},
+ Descriptor{6, 4.534357504253366e-04, {+5.536849035893876e-01, +9.989728978231427e-01}},
+ Descriptor{6, 4.182867578178122e-03, {+2.235580139510551e-01, +4.418577456047065e-01}},
+ Descriptor{6, 7.389797381822837e-04, {+8.700625815926416e-01, +4.988147208937697e-01}},
+ Descriptor{7, 1.371282553387760e-03, {+7.501439986013229e-01, +9.439853549264192e-01, +3.972114022576276e-01}},
+ Descriptor{7, 3.675472168789246e-04, {+9.025858122958752e-01, +9.846237882037482e-01, +6.307126444914427e-01}},
+ Descriptor{7, 1.258019483515844e-03, {+2.217794596478607e-01, +9.822849317288322e-02, +8.721988578583585e-01}},
+ Descriptor{7, 1.516565561694638e-03, {+9.568422563542535e-01, +2.081392891496346e-01, +4.958409266028638e-01}}};
+
+ fill_quadrature_points(descriptor_list, point_list, weight_list);
m_point_list = point_list;
m_weight_list = weight_list;
break;
}
+
default: {
- throw NormalError("Gauss quadrature formulae handle degrees up to 19 on cubes");
+ throw NormalError("Gauss quadrature formulae handle degrees up to 21 on cubes");
}
}
}
diff --git a/src/analysis/EconomicalGaussQuadrature.hpp b/src/analysis/EconomicalGaussQuadrature.hpp
index 601c76f2206b6a5191d81f0c6b888019677d0f1c..1353c1cb051670a01f0702388f878363f794852f 100644
--- a/src/analysis/EconomicalGaussQuadrature.hpp
+++ b/src/analysis/EconomicalGaussQuadrature.hpp
@@ -7,13 +7,13 @@
* Defines Economical Gauss quadrature on the reference element
* \f$]-1,1[^d\f$, where \f$d\f$ denotes the \a Dimension.
*
- * \note formulae are given by D. A. DUNAVANT
+ * \note formulae are provided by
*
* In 2d: ‘Economical symmetrical quadrature rules for complete
- * polynomials over a square domain‘ (1985).
+ * polynomials over a square domain‘ D. A. Dunavant (1985).
*
- * In 3d: 'Efficient symmetrical cubature rules for complete
- * polynomials of high degree over the unit cube' (1986).
+ * In 3d: 'High-order symmetric cubature rules for tetrahedra and
+ * pyramids' Jan Jaśkowiec & N. Sukumar (2020).
*/
template <size_t Dimension>
class EconomicalGaussQuadrature final : public QuadratureForumla<Dimension>
diff --git a/tests/test_EconomicalGaussQuadrature.cpp b/tests/test_EconomicalGaussQuadrature.cpp
index 2f20ba78941909cb8fc944f069350e6f9d755643..5c51921d4650a2a12b103e02c12b8b0fd1cbb611 100644
--- a/tests/test_EconomicalGaussQuadrature.cpp
+++ b/tests/test_EconomicalGaussQuadrature.cpp
@@ -1153,6 +1153,18 @@ TEST_CASE("EconomicalGaussQuadrature", "[analysis]")
const double z = X[2];
return p19(X) * (0.3 * x - 0.7 * y - 0.8 * z + 0.7);
};
+ auto p21 = [&p20](const TinyVector<3>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return p20(X) * (-0.9 * x + 0.2 * y + 0.5 * z - 0.6);
+ };
+ auto p22 = [&p21](const TinyVector<3>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return p21(X) * (0.3 * x - 0.6 * y - 0.7 * z + 0.2);
+ };
SECTION("degree 0 and 1")
{
@@ -1212,7 +1224,7 @@ TEST_CASE("EconomicalGaussQuadrature", "[analysis]")
REQUIRE(same_formulas(l6, l7));
- REQUIRE(l7.numberOfPoints() == 27);
+ REQUIRE(l7.numberOfPoints() == 34);
REQUIRE(integrate(p0, l7) == Catch::Approx(32));
REQUIRE(integrate(p1, l7) == Catch::Approx(-8));
@@ -1234,7 +1246,7 @@ TEST_CASE("EconomicalGaussQuadrature", "[analysis]")
REQUIRE(same_formulas(l8, l9));
- REQUIRE(l9.numberOfPoints() == 53);
+ REQUIRE(l9.numberOfPoints() == 58);
REQUIRE(integrate(p0, l9) == Catch::Approx(32));
REQUIRE(integrate(p1, l9) == Catch::Approx(-8));
@@ -1258,7 +1270,7 @@ TEST_CASE("EconomicalGaussQuadrature", "[analysis]")
REQUIRE(same_formulas(l10, l11));
- REQUIRE(l11.numberOfPoints() == 89);
+ REQUIRE(l11.numberOfPoints() == 90);
REQUIRE(integrate(p0, l11) == Catch::Approx(32));
REQUIRE(integrate(p1, l11) == Catch::Approx(-8));
@@ -1284,7 +1296,7 @@ TEST_CASE("EconomicalGaussQuadrature", "[analysis]")
REQUIRE(same_formulas(l12, l13));
- REQUIRE(l13.numberOfPoints() == 151);
+ REQUIRE(l13.numberOfPoints() == 154);
REQUIRE(integrate(p0, l13) == Catch::Approx(32));
REQUIRE(integrate(p1, l13) == Catch::Approx(-8));
@@ -1312,7 +1324,7 @@ TEST_CASE("EconomicalGaussQuadrature", "[analysis]")
REQUIRE(same_formulas(l14, l15));
- REQUIRE(l15.numberOfPoints() == 235);
+ REQUIRE(l15.numberOfPoints() == 256);
REQUIRE(integrate(p0, l15) == Catch::Approx(32));
REQUIRE(integrate(p1, l15) == Catch::Approx(-8));
@@ -1342,7 +1354,7 @@ TEST_CASE("EconomicalGaussQuadrature", "[analysis]")
REQUIRE(same_formulas(l16, l17));
- REQUIRE(l17.numberOfPoints() == 307);
+ REQUIRE(l17.numberOfPoints() == 346);
REQUIRE(integrate(p0, l17) == Catch::Approx(32));
REQUIRE(integrate(p1, l17) == Catch::Approx(-8));
@@ -1374,7 +1386,7 @@ TEST_CASE("EconomicalGaussQuadrature", "[analysis]")
REQUIRE(same_formulas(l18, l19));
- REQUIRE(l19.numberOfPoints() == 435);
+ REQUIRE(l19.numberOfPoints() == 454);
REQUIRE(integrate(p0, l19) == Catch::Approx(32));
REQUIRE(integrate(p1, l19) == Catch::Approx(-8));
@@ -1401,10 +1413,48 @@ TEST_CASE("EconomicalGaussQuadrature", "[analysis]")
REQUIRE(get_order(p20, l19, -11898262429946164522483495921. / 836985568090546875000000000.) == Catch::Approx(20));
}
+ SECTION("degree 20 and 21")
+ {
+ EconomicalGaussQuadrature<3> l20(20);
+ EconomicalGaussQuadrature<3> l21(21);
+
+ REQUIRE(same_formulas(l20, l21));
+
+ REQUIRE(l21.numberOfPoints() == 580);
+
+ REQUIRE(integrate(p0, l21) == Catch::Approx(32));
+ REQUIRE(integrate(p1, l21) == Catch::Approx(-8));
+ REQUIRE(integrate(p2, l21) == Catch::Approx(-32));
+ REQUIRE(integrate(p3, l21) == Catch::Approx(108));
+ REQUIRE(integrate(p4, l21) == Catch::Approx(868. / 75));
+ REQUIRE(integrate(p5, l21) == Catch::Approx(11176. / 225));
+ REQUIRE(integrate(p6, l21) == Catch::Approx(-34781. / 430329));
+ REQUIRE(integrate(p7, l21) == Catch::Approx(-37338109. / 4303290));
+ REQUIRE(integrate(p8, l21) == Catch::Approx(422437099. / 21516450));
+ REQUIRE(integrate(p9, l21) == Catch::Approx(-7745999747. / 358607500));
+ REQUIRE(integrate(p10, l21) == Catch::Approx(-564286973089. / 14792559375));
+ REQUIRE(integrate(p11, l21) == Catch::Approx(5047102242313. / 59170237500));
+ REQUIRE(integrate(p12, l21) == Catch::Approx(41226980237154884. / 504796088671875));
+ REQUIRE(integrate(p13, l21) == Catch::Approx(-2061220959094693133. / 26922458062500000));
+ REQUIRE(integrate(p14, l21) == Catch::Approx(9658346476244058223. / 134612290312500000));
+ REQUIRE(integrate(p15, l21) == Catch::Approx(-74949066496612419191. / 673061451562500000.));
+ REQUIRE(integrate(p16, l21) == Catch::Approx(-46230170725718969828017. / 228840893531250000000.));
+ REQUIRE(integrate(p17, l21) == Catch::Approx(2946902633453348474221. / 190700744609375000000.));
+ REQUIRE(integrate(p18, l21) == Catch::Approx(904723313909284441962799. / 47555998186962890625000.));
+ REQUIRE(integrate(p19, l21) == Catch::Approx(-91477977618647751170958517. / 22826879129742187500000000.));
+ REQUIRE(integrate(p20, l21) == Catch::Approx(-11898262429946164522483495921. / 836985568090546875000000000.));
+ REQUIRE(integrate(p21, l21) == Catch::Approx(7694483338683700814691225463. / 224192562881396484375000000.));
+ REQUIRE(integrate(p22, l21) !=
+ Catch::Approx(10761048146311587678467825954981. / 481266701652064453125000000000.));
+
+ REQUIRE(get_order(p22, l21, 10761048146311587678467825954981. / 481266701652064453125000000000.) ==
+ Catch::Approx(22));
+ }
+
SECTION("not implemented formulae")
{
- REQUIRE_THROWS_WITH(EconomicalGaussQuadrature<3>(20),
- "error: Gauss quadrature formulae handle degrees up to 19 on cubes");
+ REQUIRE_THROWS_WITH(EconomicalGaussQuadrature<3>(22),
+ "error: Gauss quadrature formulae handle degrees up to 21 on cubes");
}
}