diff --git a/src/analysis/EconomicalGaussQuadrature.cpp b/src/analysis/EconomicalGaussQuadrature.cpp
index 27ada266d21ce6cdb2628b00faad7c17e7cae3f6..c3033f89b1bc50b752eaf800f622b9bf6dc0c91b 100644
--- a/src/analysis/EconomicalGaussQuadrature.cpp
+++ b/src/analysis/EconomicalGaussQuadrature.cpp
@@ -16,15 +16,87 @@ template <>
void
EconomicalGaussQuadrature<2>::_buildPointAndWeightLists()
{
+ using R2 = TinyVector<2>;
+ using QuadratureGroup = std::array<double, 3>;
+
+ auto fill_quadrature_points = [](auto quadrature_group_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) {
+ auto [w, xi, eta] = quadrature_group_list[i];
+ if (xi != eta) {
+ Assert(xi != 0, "quadrature group was not defined correctly");
+ if (eta == 0) {
+ point_list[k + 0] = {xi, eta};
+ point_list[k + 1] = {-xi, eta};
+ point_list[k + 2] = {eta, xi};
+ point_list[k + 3] = {eta, -xi};
+
+ weight_list[k + 0] = w;
+ weight_list[k + 1] = w;
+ weight_list[k + 2] = w;
+ weight_list[k + 3] = w;
+
+ k += 4;
+ } else {
+ point_list[k + 0] = {xi, eta};
+ point_list[k + 1] = {-xi, eta};
+ point_list[k + 2] = {xi, -eta};
+ point_list[k + 3] = {-xi, -eta};
+ point_list[k + 4] = {eta, xi};
+ point_list[k + 5] = {-eta, xi};
+ point_list[k + 6] = {eta, -xi};
+ point_list[k + 7] = {-eta, -xi};
+
+ weight_list[k + 0] = w;
+ weight_list[k + 1] = w;
+ weight_list[k + 2] = w;
+ weight_list[k + 3] = w;
+ weight_list[k + 4] = w;
+ weight_list[k + 5] = w;
+ weight_list[k + 6] = w;
+ weight_list[k + 7] = w;
+
+ k += 8;
+ }
+ } else {
+ if (xi == 0) {
+ point_list[k] = {xi, eta};
+
+ weight_list[k] = w;
+
+ k += 1;
+ } else {
+ point_list[k + 0] = {xi, eta};
+ point_list[k + 1] = {-xi, eta};
+ point_list[k + 2] = {xi, -eta};
+ point_list[k + 3] = {-xi, -eta};
+
+ weight_list[k + 0] = w;
+ weight_list[k + 1] = w;
+ weight_list[k + 2] = w;
+ weight_list[k + 3] = w;
+
+ k += 4;
+ }
+ }
+ }
+
+ Assert(k == point_list.size(), "invalid number of quadrature points");
+ };
+
switch (m_degree) {
case 0:
case 1: {
constexpr size_t nb_points = 1;
- SmallArray<TinyVector<2>> point_list(nb_points);
+ SmallArray<R2> point_list(nb_points);
SmallArray<double> weight_list(nb_points);
- point_list[0] = {0, 0};
- weight_list[0] = 4;
+ std::array<QuadratureGroup, 1> quadrature_group_list = //
+ {QuadratureGroup{4.000000000000000, 0.000000000000000, 0.000000000000000}};
+
+ fill_quadrature_points(quadrature_group_list, point_list, weight_list);
m_point_list = point_list;
m_weight_list = weight_list;
@@ -33,18 +105,13 @@ EconomicalGaussQuadrature<2>::_buildPointAndWeightLists()
case 2:
case 3: {
constexpr size_t nb_points = 4;
- SmallArray<TinyVector<2>> point_list(nb_points);
+ SmallArray<R2> point_list(nb_points);
SmallArray<double> weight_list(nb_points);
- point_list[0] = {+0.577350269189626, +0.577350269189626};
- point_list[1] = {+0.577350269189626, -0.577350269189626};
- point_list[2] = {-0.577350269189626, +0.577350269189626};
- point_list[3] = {-0.577350269189626, -0.577350269189626};
+ std::array<QuadratureGroup, 1> quadrature_group_list = //
+ {QuadratureGroup{1.000000000000000, 0.577350269189626, 0.577350269189626}};
- weight_list[0] = 1;
- weight_list[1] = 1;
- weight_list[2] = 1;
- weight_list[3] = 1;
+ fill_quadrature_points(quadrature_group_list, point_list, weight_list);
m_point_list = point_list;
m_weight_list = weight_list;
@@ -53,26 +120,14 @@ EconomicalGaussQuadrature<2>::_buildPointAndWeightLists()
case 4:
case 5: {
constexpr size_t nb_points = 8;
- SmallArray<TinyVector<2>> point_list(nb_points);
+ SmallArray<R2> point_list(nb_points);
SmallArray<double> weight_list(nb_points);
- point_list[0] = {+0.683130051063973, +0.000000000000000};
- point_list[1] = {-0.683130051063973, +0.000000000000000};
- point_list[2] = {+0.000000000000000, +0.683130051063973};
- point_list[3] = {+0.000000000000000, -0.683130051063973};
- point_list[4] = {+0.881917103688197, +0.881917103688197};
- point_list[5] = {+0.881917103688197, -0.881917103688197};
- point_list[6] = {-0.881917103688197, +0.881917103688197};
- point_list[7] = {-0.881917103688197, -0.881917103688197};
-
- weight_list[0] = 0.816326530612245;
- weight_list[1] = 0.816326530612245;
- weight_list[2] = 0.816326530612245;
- weight_list[3] = 0.816326530612245;
- weight_list[4] = 0.183673469387755;
- weight_list[5] = 0.183673469387755;
- weight_list[6] = 0.183673469387755;
- weight_list[7] = 0.183673469387755;
+ std::array quadrature_group_list = //
+ {QuadratureGroup{0.816326530612245, 0.683130051063973, 0.000000000000000},
+ QuadratureGroup{0.183673469387755, 0.881917103688197, 0.881917103688197}};
+
+ fill_quadrature_points(quadrature_group_list, point_list, weight_list);
m_point_list = point_list;
m_weight_list = weight_list;
@@ -81,41 +136,182 @@ EconomicalGaussQuadrature<2>::_buildPointAndWeightLists()
case 6:
case 7: {
constexpr size_t nb_points = 12;
- SmallArray<TinyVector<2>> point_list(nb_points);
+ SmallArray<R2> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ std::array quadrature_group_list = //
+ {QuadratureGroup{0.241975308641975, 0.925820099772551, 0.000000000000000},
+ QuadratureGroup{0.237431774690630, 0.805979782918599, 0.805979782918599},
+ QuadratureGroup{0.520592916667394, 0.380554433208316, 0.380554433208316}};
+
+ fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 8:
+ case 9: {
+ constexpr size_t nb_points = 20;
+ SmallArray<R2> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ std::array quadrature_group_list = //
+ {QuadratureGroup{0.018475842507491, 1.121225763866564, 0.000000000000000},
+ QuadratureGroup{0.390052939160735, 0.451773049920657, 0.000000000000000},
+ QuadratureGroup{0.083095178026482, 0.891849420851512, 0.891849420851512},
+ QuadratureGroup{0.254188020152646, 0.824396370149276, 0.411623426336542}};
+
+ fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 10:
+ case 11: {
+ constexpr size_t nb_points = 25;
+ SmallArray<R2> point_list(nb_points);
SmallArray<double> weight_list(nb_points);
- point_list[0] = {+0.925820099772551, +0.000000000000000};
- point_list[1] = {-0.925820099772551, +0.000000000000000};
- point_list[2] = {+0.000000000000000, +0.925820099772551};
- point_list[3] = {+0.000000000000000, -0.925820099772551};
- point_list[4] = {+0.805979782918599, +0.805979782918599};
- point_list[5] = {+0.805979782918599, -0.805979782918599};
- point_list[6] = {-0.805979782918599, +0.805979782918599};
- point_list[7] = {-0.805979782918599, -0.805979782918599};
- point_list[8] = {+0.380554433208316, +0.380554433208316};
- point_list[9] = {+0.380554433208316, -0.380554433208316};
- point_list[10] = {-0.380554433208316, +0.380554433208316};
- point_list[11] = {-0.380554433208316, -0.380554433208316};
-
- weight_list[0] = 0.241975308641975;
- weight_list[1] = 0.241975308641975;
- weight_list[2] = 0.241975308641975;
- weight_list[3] = 0.241975308641975;
- weight_list[4] = 0.237431774690630;
- weight_list[5] = 0.237431774690630;
- weight_list[6] = 0.237431774690630;
- weight_list[7] = 0.237431774690630;
- weight_list[8] = 0.520592916667394;
- weight_list[9] = 0.520592916667394;
- weight_list[10] = 0.520592916667394;
- weight_list[11] = 0.520592916667394;
+ std::array quadrature_group_list = //
+ {QuadratureGroup{0.365379525585903, 0.000000000000000, 0.000000000000000},
+ QuadratureGroup{0.027756165564204, 1.044402915409813, 0.000000000000000},
+ QuadratureGroup{0.244272057751754, 0.769799068396649, 0.000000000000000},
+ QuadratureGroup{0.034265103851229, 0.935787012440540, 0.935787012440540},
+ QuadratureGroup{0.308993036133713, 0.413491953449114, 0.413491953449114},
+ QuadratureGroup{0.146684377651312, 0.883025508525690, 0.575653595840465}};
+
+ fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 12:
+ case 13: {
+ constexpr size_t nb_points = 36;
+ SmallArray<R2> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ std::array quadrature_group_list = //
+ {QuadratureGroup{0.005656169693764, 1.086056158573971, 0.000000000000000},
+ QuadratureGroup{0.192443867470396, 0.658208197042585, 0.000000000000000},
+ QuadratureGroup{0.005166832979773, 1.001300602991729, 1.001300602991729},
+ QuadratureGroup{0.200302559622138, 0.584636168775946, 0.584636168775946},
+ QuadratureGroup{0.228125175912536, 0.246795612720261, 0.246795612720261},
+ QuadratureGroup{0.117496926974491, 0.900258815287201, 0.304720678579870},
+ QuadratureGroup{0.066655770186205, 0.929866705560780, 0.745052720131169}};
+
+ fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 14:
+ case 15: {
+ constexpr size_t nb_points = 45;
+ SmallArray<R2> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ std::array quadrature_group_list = //
+ {QuadratureGroup{-0.001768979827207, 0.000000000000000, 0.000000000000000},
+ QuadratureGroup{+0.012816726617512, 1.027314357719367, 0.000000000000000},
+ QuadratureGroup{+0.119897873101347, 0.856766776147643, 0.000000000000000},
+ QuadratureGroup{+0.210885452208801, 0.327332998189723, 0.000000000000000},
+ QuadratureGroup{+0.006392720128215, 0.967223740028505, 0.967223740028505},
+ QuadratureGroup{+0.104415680788580, 0.732168901749711, 0.732168901749711},
+ QuadratureGroup{+0.168053047203816, 0.621974427996805, 0.321696694921009},
+ QuadratureGroup{+0.076169694452294, 0.928618480068352, 0.455124178121179},
+ QuadratureGroup{+0.028794154400064, 0.960457474887516, 0.809863684081217}};
+
+ fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 16:
+ case 17: {
+ constexpr size_t nb_points = 60;
+ SmallArray<R2> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ std::array quadrature_group_list = //
+ {QuadratureGroup{0.020614915919991, 0.989353074512600, 0.000000000000000},
+ QuadratureGroup{0.128025716179910, 0.376285207157973, 0.000000000000000},
+ QuadratureGroup{0.005511739534032, 0.978848279262233, 0.978848279262233},
+ QuadratureGroup{0.039207712457142, 0.885794729164116, 0.885794729164116},
+ QuadratureGroup{0.076396945079863, 0.171756123838348, 0.171756123838348},
+ QuadratureGroup{0.141513729949972, 0.590499273806002, 0.319505036634574},
+ QuadratureGroup{0.083903279363798, 0.799079131916863, 0.597972451929457},
+ QuadratureGroup{0.060394163649685, 0.803743962958745, 0.058344481776551},
+ QuadratureGroup{0.057387752969213, 0.936506276127495, 0.347386316166203},
+ QuadratureGroup{0.021922559481864, 0.981321179805452, 0.706000287798646}};
+
+ fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 18:
+ case 19: {
+ constexpr size_t nb_points = 72;
+ SmallArray<R2> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ std::array quadrature_group_list = //
+ {QuadratureGroup{0.038205406871462, 0.943962831808239, 0.000000000000000},
+ QuadratureGroup{0.135368502976521, 0.536918434376013, 0.000000000000000},
+ QuadratureGroup{0.005773708558664, 0.973981076394170, 0.973981076394170},
+ QuadratureGroup{0.067460759759473, 0.742995535327609, 0.742995535327609},
+ QuadratureGroup{0.140899115227892, 0.285010052188916, 0.285010052188916},
+ QuadratureGroup{0.047466627685662, 0.068354569272491, 0.068354569272491},
+ QuadratureGroup{0.078619467342982, 0.802952004398543, 0.203345534163332},
+ QuadratureGroup{0.094979169511394, 0.634244672807882, 0.426572172992877},
+ QuadratureGroup{0.022331162356015, 0.978350706908227, 0.295830776620995},
+ QuadratureGroup{0.055594877793785, 0.901672714410389, 0.541983037327871},
+ QuadratureGroup{0.006049054506376, 1.007018449383116, 0.669414798783936},
+ QuadratureGroup{0.024839207949609, 0.945161453573471, 0.829501421477824}};
+
+ 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 = 88;
+ SmallArray<R2> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ std::array quadrature_group_list = //
+ {QuadratureGroup{0.019503841092684, 0.980883148832881, 0.000000000000000},
+ QuadratureGroup{0.089012127744268, 0.678152700336576, 0.000000000000000},
+ QuadratureGroup{0.114568584702749, 0.240599282275864, 0.000000000000000},
+ QuadratureGroup{0.007463627359106, 0.965176994929162, 0.965176994929162},
+ QuadratureGroup{0.050585943594705, 0.749698539312765, 0.749698539312765},
+ QuadratureGroup{0.074613865184212, 0.568983925500818, 0.568983925500818},
+ QuadratureGroup{0.023501091310143, 0.971086142843168, 0.355832132274584},
+ QuadratureGroup{0.011588562644144, 0.983453947854968, 0.645588139196562},
+ QuadratureGroup{0.023073245798171, 0.933927707027213, 0.821920249234369},
+ QuadratureGroup{0.001570221774472, 1.014086498915039, 0.862185099566557},
+ QuadratureGroup{0.049102258016277, 0.877914842155496, 0.168914072450263},
+ QuadratureGroup{0.042512352239126, 0.882246882640128, 0.568113580166780},
+ QuadratureGroup{0.067270936863160, 0.741324453314596, 0.371360260002223},
+ QuadratureGroup{0.103507336515645, 0.469570217710647, 0.237333359193547}};
+
+ fill_quadrature_points(quadrature_group_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 7 on quadrangles");
+ throw NormalError("Gauss quadrature formulae handle degrees up to 21 on squares");
}
}
}
@@ -270,7 +466,7 @@ EconomicalGaussQuadrature<3>::_buildPointAndWeightLists()
break;
}
default: {
- throw NormalError("Gauss quadrature formulae handle degrees up to 7 on hexahedra");
+ throw NormalError("Gauss quadrature formulae handle degrees up to 7 on cubes");
}
}
}
diff --git a/src/analysis/EconomicalGaussQuadrature.hpp b/src/analysis/EconomicalGaussQuadrature.hpp
index 0862d0a66486d23c99ae9a0525302673bb2e99b4..601c76f2206b6a5191d81f0c6b888019677d0f1c 100644
--- a/src/analysis/EconomicalGaussQuadrature.hpp
+++ b/src/analysis/EconomicalGaussQuadrature.hpp
@@ -7,7 +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 extracted from High-order Finite Element Method [2004 - Chapman & Hall]
+ * \note formulae are given by D. A. DUNAVANT
+ *
+ * In 2d: ‘Economical symmetrical quadrature rules for complete
+ * polynomials over a square domain‘ (1985).
+ *
+ * In 3d: 'Efficient symmetrical cubature rules for complete
+ * polynomials of high degree over the unit cube' (1986).
*/
template <size_t Dimension>
class EconomicalGaussQuadrature final : public QuadratureForumla<Dimension>
diff --git a/tests/test_EconomicalGaussQuadrature.cpp b/tests/test_EconomicalGaussQuadrature.cpp
index 1da84db8d901c65e4ac092bcb5595b504fdb284d..18b7143cc42b685274a7ddeb4f1045ac824fdbec 100644
--- a/tests/test_EconomicalGaussQuadrature.cpp
+++ b/tests/test_EconomicalGaussQuadrature.cpp
@@ -604,6 +604,76 @@ TEST_CASE("EconomicalGaussQuadrature", "[analysis]")
const double y = X[1];
return p7(X) * (2 * x - 3 * y - 3);
};
+ auto p9 = [&p8](const TinyVector<2>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ return p8(X) * (x + 2 * y - 1.7);
+ };
+ auto p10 = [&p9](const TinyVector<2>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ return p9(X) * (-1.3 * x - 1.7 * y + 1.3);
+ };
+ auto p11 = [&p10](const TinyVector<2>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ return p10(X) * (0.8 * x - 3.1 * y + 0.6);
+ };
+ auto p12 = [&p11](const TinyVector<2>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ return p11(X) * (1.8 * x + 1.3 * y - 0.3);
+ };
+ auto p13 = [&p12](const TinyVector<2>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ return p12(X) * (-0.9 * x + 1.1 * y - 0.6);
+ };
+ auto p14 = [&p13](const TinyVector<2>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ return p13(X) * (0.6 * x + 0.3 * y + 1.1);
+ };
+ auto p15 = [&p14](const TinyVector<2>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ return p14(X) * (0.5 * x - 0.4 * y - 1.1);
+ };
+ auto p16 = [&p15](const TinyVector<2>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ return p15(X) * (-0.3 * x - 0.3 * y - 0.2);
+ };
+ auto p17 = [&p16](const TinyVector<2>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ return p16(X) * (-0.1 * x - 0.4 * y - 0.3);
+ };
+ auto p18 = [&p17](const TinyVector<2>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ return p17(X) * (0.2 * x + 0.3 * y + 0.3);
+ };
+ auto p19 = [&p18](const TinyVector<2>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ return p18(X) * (2.1 * x + 3.3 * y - 0.3);
+ };
+ auto p20 = [&p19](const TinyVector<2>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ return p19(X) * (1.2 * x - 2.1 * y + 0.6);
+ };
+ auto p21 = [&p20](const TinyVector<2>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ return p20(X) * (-1.3 * x - 1.2 * y + 0.7);
+ };
+ auto p22 = [&p21](const TinyVector<2>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ return p21(X) * (1.1 * x - 2.2 * y - 0.3);
+ };
SECTION("degree 0 and 1")
{
@@ -678,10 +748,221 @@ TEST_CASE("EconomicalGaussQuadrature", "[analysis]")
REQUIRE(get_order(p8, l7, 957697. / 175) == Catch::Approx(8));
}
+ SECTION("degree 8 and 9")
+ {
+ EconomicalGaussQuadrature<2> l8(8);
+ EconomicalGaussQuadrature<2> l9(9);
+
+ REQUIRE(same_formulas(l8, l9));
+
+ REQUIRE(l9.numberOfPoints() == 20);
+
+ REQUIRE(integrate(p0, l9) == Catch::Approx(8));
+ REQUIRE(integrate(p1, l9) == Catch::Approx(4));
+ REQUIRE(integrate(p2, l9) == Catch::Approx(20. / 3));
+ REQUIRE(integrate(p3, l9) == Catch::Approx(-13));
+ REQUIRE(integrate(p4, l9) == Catch::Approx(-1184. / 15));
+ REQUIRE(integrate(p5, l9) == Catch::Approx(1971. / 25));
+ REQUIRE(integrate(p6, l9) == Catch::Approx(60441. / 175));
+ REQUIRE(integrate(p7, l9) == Catch::Approx(-1307119. / 525));
+ REQUIRE(integrate(p8, l9) == Catch::Approx(957697. / 175));
+ REQUIRE(integrate(p9, l9) == Catch::Approx(-196981. / 5250));
+ REQUIRE(integrate(p10, l9) != Catch::Approx(78447601. / 577500));
+
+ REQUIRE(get_order(p10, l9, 78447601. / 577500) == Catch::Approx(10));
+ }
+
+ SECTION("degree 10 and 11")
+ {
+ EconomicalGaussQuadrature<2> l10(10);
+ EconomicalGaussQuadrature<2> l11(11);
+
+ REQUIRE(same_formulas(l10, l11));
+
+ REQUIRE(l11.numberOfPoints() == 25);
+
+ REQUIRE(integrate(p0, l11) == Catch::Approx(8));
+ REQUIRE(integrate(p1, l11) == Catch::Approx(4));
+ REQUIRE(integrate(p2, l11) == Catch::Approx(20. / 3));
+ REQUIRE(integrate(p3, l11) == Catch::Approx(-13));
+ REQUIRE(integrate(p4, l11) == Catch::Approx(-1184. / 15));
+ REQUIRE(integrate(p5, l11) == Catch::Approx(1971. / 25));
+ REQUIRE(integrate(p6, l11) == Catch::Approx(60441. / 175));
+ REQUIRE(integrate(p7, l11) == Catch::Approx(-1307119. / 525));
+ REQUIRE(integrate(p8, l11) == Catch::Approx(957697. / 175));
+ REQUIRE(integrate(p9, l11) == Catch::Approx(-196981. / 5250));
+ REQUIRE(integrate(p10, l11) == Catch::Approx(78447601. / 577500));
+ REQUIRE(integrate(p11, l11) == Catch::Approx(673235482069. / 86625000));
+ REQUIRE(integrate(p12, l11) != Catch::Approx(-4092600398251. / 303187500));
+
+ REQUIRE(get_order(p12, l11, -4092600398251. / 303187500) == Catch::Approx(12));
+ }
+
+ SECTION("degree 12 and 13")
+ {
+ EconomicalGaussQuadrature<2> l12(12);
+ EconomicalGaussQuadrature<2> l13(13);
+
+ REQUIRE(same_formulas(l12, l13));
+
+ REQUIRE(l13.numberOfPoints() == 36);
+
+ REQUIRE(integrate(p0, l13) == Catch::Approx(8));
+ REQUIRE(integrate(p1, l13) == Catch::Approx(4));
+ REQUIRE(integrate(p2, l13) == Catch::Approx(20. / 3));
+ REQUIRE(integrate(p3, l13) == Catch::Approx(-13));
+ REQUIRE(integrate(p4, l13) == Catch::Approx(-1184. / 15));
+ REQUIRE(integrate(p5, l13) == Catch::Approx(1971. / 25));
+ REQUIRE(integrate(p6, l13) == Catch::Approx(60441. / 175));
+ REQUIRE(integrate(p7, l13) == Catch::Approx(-1307119. / 525));
+ REQUIRE(integrate(p8, l13) == Catch::Approx(957697. / 175));
+ REQUIRE(integrate(p9, l13) == Catch::Approx(-196981. / 5250));
+ REQUIRE(integrate(p10, l13) == Catch::Approx(78447601. / 577500));
+ REQUIRE(integrate(p11, l13) == Catch::Approx(673235482069. / 86625000));
+ REQUIRE(integrate(p12, l13) == Catch::Approx(-4092600398251. / 303187500));
+ REQUIRE(integrate(p13, l13) == Catch::Approx(77231697272647. / 5053125000));
+ REQUIRE(integrate(p14, l13) != Catch::Approx(46574962939049. / 9384375000));
+
+ REQUIRE(get_order(p14, l13, 46574962939049. / 9384375000) == Catch::Approx(14));
+ }
+
+ SECTION("degree 14 and 15")
+ {
+ EconomicalGaussQuadrature<2> l14(14);
+ EconomicalGaussQuadrature<2> l15(15);
+
+ REQUIRE(same_formulas(l14, l15));
+
+ REQUIRE(l15.numberOfPoints() == 45);
+
+ REQUIRE(integrate(p0, l15) == Catch::Approx(8));
+ REQUIRE(integrate(p1, l15) == Catch::Approx(4));
+ REQUIRE(integrate(p2, l15) == Catch::Approx(20. / 3));
+ REQUIRE(integrate(p3, l15) == Catch::Approx(-13));
+ REQUIRE(integrate(p4, l15) == Catch::Approx(-1184. / 15));
+ REQUIRE(integrate(p5, l15) == Catch::Approx(1971. / 25));
+ REQUIRE(integrate(p6, l15) == Catch::Approx(60441. / 175));
+ REQUIRE(integrate(p7, l15) == Catch::Approx(-1307119. / 525));
+ REQUIRE(integrate(p8, l15) == Catch::Approx(957697. / 175));
+ REQUIRE(integrate(p9, l15) == Catch::Approx(-196981. / 5250));
+ REQUIRE(integrate(p10, l15) == Catch::Approx(78447601. / 577500));
+ REQUIRE(integrate(p11, l15) == Catch::Approx(673235482069. / 86625000));
+ REQUIRE(integrate(p12, l15) == Catch::Approx(-4092600398251. / 303187500));
+ REQUIRE(integrate(p13, l15) == Catch::Approx(77231697272647. / 5053125000));
+ REQUIRE(integrate(p14, l15) == Catch::Approx(46574962939049. / 9384375000));
+ REQUIRE(integrate(p15, l15) == Catch::Approx(-4818864487842259. / 1094843750000));
+ REQUIRE(integrate(p16, l15) != Catch::Approx(-89885697514686141. / 23265429687500));
+
+ REQUIRE(get_order(p16, l15, -89885697514686141. / 23265429687500) == Catch::Approx(16));
+ }
+
+ SECTION("degree 16 and 17")
+ {
+ EconomicalGaussQuadrature<2> l16(16);
+ EconomicalGaussQuadrature<2> l17(17);
+
+ REQUIRE(same_formulas(l16, l17));
+
+ REQUIRE(l17.numberOfPoints() == 60);
+
+ REQUIRE(integrate(p0, l17) == Catch::Approx(8));
+ REQUIRE(integrate(p1, l17) == Catch::Approx(4));
+ REQUIRE(integrate(p2, l17) == Catch::Approx(20. / 3));
+ REQUIRE(integrate(p3, l17) == Catch::Approx(-13));
+ REQUIRE(integrate(p4, l17) == Catch::Approx(-1184. / 15));
+ REQUIRE(integrate(p5, l17) == Catch::Approx(1971. / 25));
+ REQUIRE(integrate(p6, l17) == Catch::Approx(60441. / 175));
+ REQUIRE(integrate(p7, l17) == Catch::Approx(-1307119. / 525));
+ REQUIRE(integrate(p8, l17) == Catch::Approx(957697. / 175));
+ REQUIRE(integrate(p9, l17) == Catch::Approx(-196981. / 5250));
+ REQUIRE(integrate(p10, l17) == Catch::Approx(78447601. / 577500));
+ REQUIRE(integrate(p11, l17) == Catch::Approx(673235482069. / 86625000));
+ REQUIRE(integrate(p12, l17) == Catch::Approx(-4092600398251. / 303187500));
+ REQUIRE(integrate(p13, l17) == Catch::Approx(77231697272647. / 5053125000));
+ REQUIRE(integrate(p14, l17) == Catch::Approx(46574962939049. / 9384375000));
+ REQUIRE(integrate(p15, l17) == Catch::Approx(-4818864487842259. / 1094843750000));
+ REQUIRE(integrate(p16, l17) == Catch::Approx(-89885697514686141. / 23265429687500));
+ REQUIRE(integrate(p17, l17) == Catch::Approx(16706355156097391. / 12085937500000));
+ REQUIRE(integrate(p18, l17) != Catch::Approx(495866230514635109957. / 397838847656250000.).epsilon(1E-8));
+
+ REQUIRE(get_order(p18, l17, 495866230514635109957. / 397838847656250000.) == Catch::Approx(18).margin(0.001));
+ }
+
+ SECTION("degree 18 and 19")
+ {
+ EconomicalGaussQuadrature<2> l18(18);
+ EconomicalGaussQuadrature<2> l19(19);
+
+ REQUIRE(same_formulas(l18, l19));
+
+ REQUIRE(l19.numberOfPoints() == 72);
+
+ REQUIRE(integrate(p0, l19) == Catch::Approx(8));
+ REQUIRE(integrate(p1, l19) == Catch::Approx(4));
+ REQUIRE(integrate(p2, l19) == Catch::Approx(20. / 3));
+ REQUIRE(integrate(p3, l19) == Catch::Approx(-13));
+ REQUIRE(integrate(p4, l19) == Catch::Approx(-1184. / 15));
+ REQUIRE(integrate(p5, l19) == Catch::Approx(1971. / 25));
+ REQUIRE(integrate(p6, l19) == Catch::Approx(60441. / 175));
+ REQUIRE(integrate(p7, l19) == Catch::Approx(-1307119. / 525));
+ REQUIRE(integrate(p8, l19) == Catch::Approx(957697. / 175));
+ REQUIRE(integrate(p9, l19) == Catch::Approx(-196981. / 5250));
+ REQUIRE(integrate(p10, l19) == Catch::Approx(78447601. / 577500));
+ REQUIRE(integrate(p11, l19) == Catch::Approx(673235482069. / 86625000));
+ REQUIRE(integrate(p12, l19) == Catch::Approx(-4092600398251. / 303187500));
+ REQUIRE(integrate(p13, l19) == Catch::Approx(77231697272647. / 5053125000));
+ REQUIRE(integrate(p14, l19) == Catch::Approx(46574962939049. / 9384375000));
+ REQUIRE(integrate(p15, l19) == Catch::Approx(-4818864487842259. / 1094843750000));
+ REQUIRE(integrate(p16, l19) == Catch::Approx(-89885697514686141. / 23265429687500));
+ REQUIRE(integrate(p17, l19) == Catch::Approx(16706355156097391. / 12085937500000));
+ REQUIRE(integrate(p18, l19) == Catch::Approx(495866230514635109957. / 397838847656250000.));
+ REQUIRE(integrate(p19, l19) == Catch::Approx(703712939580204375319. / 132612949218750000.));
+ REQUIRE(integrate(p20, l19) != Catch::Approx(-891851528496270127477. / 884086328125000000.));
+
+ REQUIRE(get_order(p20, l19, -891851528496270127477. / 884086328125000000.) == Catch::Approx(20).margin(0.001));
+ }
+
+ SECTION("degree 20 and 21")
+ {
+ EconomicalGaussQuadrature<2> l20(20);
+ EconomicalGaussQuadrature<2> l21(21);
+
+ REQUIRE(same_formulas(l20, l21));
+
+ REQUIRE(l21.numberOfPoints() == 88);
+
+ REQUIRE(integrate(p0, l21) == Catch::Approx(8));
+ REQUIRE(integrate(p1, l21) == Catch::Approx(4));
+ REQUIRE(integrate(p2, l21) == Catch::Approx(20. / 3));
+ REQUIRE(integrate(p3, l21) == Catch::Approx(-13));
+ REQUIRE(integrate(p4, l21) == Catch::Approx(-1184. / 15));
+ REQUIRE(integrate(p5, l21) == Catch::Approx(1971. / 25));
+ REQUIRE(integrate(p6, l21) == Catch::Approx(60441. / 175));
+ REQUIRE(integrate(p7, l21) == Catch::Approx(-1307119. / 525));
+ REQUIRE(integrate(p8, l21) == Catch::Approx(957697. / 175));
+ REQUIRE(integrate(p9, l21) == Catch::Approx(-196981. / 5250));
+ REQUIRE(integrate(p10, l21) == Catch::Approx(78447601. / 577500));
+ REQUIRE(integrate(p11, l21) == Catch::Approx(673235482069. / 86625000));
+ REQUIRE(integrate(p12, l21) == Catch::Approx(-4092600398251. / 303187500));
+ REQUIRE(integrate(p13, l21) == Catch::Approx(77231697272647. / 5053125000));
+ REQUIRE(integrate(p14, l21) == Catch::Approx(46574962939049. / 9384375000));
+ REQUIRE(integrate(p15, l21) == Catch::Approx(-4818864487842259. / 1094843750000));
+ REQUIRE(integrate(p16, l21) == Catch::Approx(-89885697514686141. / 23265429687500));
+ REQUIRE(integrate(p17, l21) == Catch::Approx(16706355156097391. / 12085937500000));
+ REQUIRE(integrate(p18, l21) == Catch::Approx(495866230514635109957. / 397838847656250000.));
+ REQUIRE(integrate(p19, l21) == Catch::Approx(703712939580204375319. / 132612949218750000.));
+ REQUIRE(integrate(p20, l21) == Catch::Approx(-891851528496270127477. / 884086328125000000.));
+ REQUIRE(integrate(p21, l21) == Catch::Approx(31710268999580650802107. / 66306474609375000000.));
+ REQUIRE(integrate(p22, l21) != Catch::Approx(-1827205780869627586647799. / 726213769531250000000.).epsilon(1E-8));
+
+ REQUIRE(get_order(p22, l21, -1827205780869627586647799. / 726213769531250000000.) ==
+ Catch::Approx(22).margin(0.01));
+ }
+
SECTION("not implemented formulae")
{
- REQUIRE_THROWS_WITH(EconomicalGaussQuadrature<2>(8),
- "error: Gauss quadrature formulae handle degrees up to 7 on quadrangles");
+ REQUIRE_THROWS_WITH(EconomicalGaussQuadrature<2>(22),
+ "error: Gauss quadrature formulae handle degrees up to 21 on squares");
}
}
@@ -856,7 +1137,7 @@ TEST_CASE("EconomicalGaussQuadrature", "[analysis]")
SECTION("not implemented formulae")
{
REQUIRE_THROWS_WITH(EconomicalGaussQuadrature<3>(8),
- "error: Gauss quadrature formulae handle degrees up to 7 on hexahedra");
+ "error: Gauss quadrature formulae handle degrees up to 7 on cubes");
}
}