From bb66aff13a6e8217f71241fce89b32bafa81cab4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Thu, 21 Oct 2021 16:03:16 +0200
Subject: [PATCH] Add economical Gauss quadrature for square and cube (up to
degree 7)
---
src/analysis/CMakeLists.txt | 1 +
src/analysis/EconomicalGaussQuadrature.cpp | 276 +++++++
src/analysis/EconomicalGaussQuadrature.hpp | 31 +
tests/CMakeLists.txt | 1 +
tests/test_EconomicalGaussQuadrature.cpp | 877 +++++++++++++++++++++
5 files changed, 1186 insertions(+)
create mode 100644 src/analysis/EconomicalGaussQuadrature.cpp
create mode 100644 src/analysis/EconomicalGaussQuadrature.hpp
create mode 100644 tests/test_EconomicalGaussQuadrature.cpp
diff --git a/src/analysis/CMakeLists.txt b/src/analysis/CMakeLists.txt
index bb934db37..5e30087d9 100644
--- a/src/analysis/CMakeLists.txt
+++ b/src/analysis/CMakeLists.txt
@@ -2,6 +2,7 @@
add_library(
PugsAnalysis
+ EconomicalGaussQuadrature.cpp
SimplicialGaussQuadrature.cpp
TensorialGaussLegendreQuadrature.cpp
TensorialGaussLobattoQuadrature.cpp)
diff --git a/src/analysis/EconomicalGaussQuadrature.cpp b/src/analysis/EconomicalGaussQuadrature.cpp
new file mode 100644
index 000000000..27ada266d
--- /dev/null
+++ b/src/analysis/EconomicalGaussQuadrature.cpp
@@ -0,0 +1,276 @@
+#include <analysis/EconomicalGaussQuadrature.hpp>
+#include <analysis/TensorialGaussLegendreQuadrature.hpp>
+#include <utils/Exceptions.hpp>
+
+template <>
+void
+EconomicalGaussQuadrature<1>::_buildPointAndWeightLists()
+{
+ TensorialGaussLegendreQuadrature<1> gauss(m_degree);
+
+ m_point_list = gauss.pointList();
+ m_weight_list = gauss.weightList();
+}
+
+template <>
+void
+EconomicalGaussQuadrature<2>::_buildPointAndWeightLists()
+{
+ switch (m_degree) {
+ case 0:
+ case 1: {
+ constexpr size_t nb_points = 1;
+ SmallArray<TinyVector<2>> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ point_list[0] = {0, 0};
+ weight_list[0] = 4;
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 2:
+ case 3: {
+ constexpr size_t nb_points = 4;
+ SmallArray<TinyVector<2>> 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};
+
+ weight_list[0] = 1;
+ weight_list[1] = 1;
+ weight_list[2] = 1;
+ weight_list[3] = 1;
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 4:
+ case 5: {
+ constexpr size_t nb_points = 8;
+ SmallArray<TinyVector<2>> 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;
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 6:
+ case 7: {
+ constexpr size_t nb_points = 12;
+ SmallArray<TinyVector<2>> 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;
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ default: {
+ throw NormalError("Gauss quadrature formulae handle degrees up to 7 on quadrangles");
+ }
+ }
+}
+
+template <>
+void
+EconomicalGaussQuadrature<3>::_buildPointAndWeightLists()
+{
+ switch (m_degree) {
+ 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, 0, 0};
+
+ weight_list[0] = 8;
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 2:
+ case 3: {
+ constexpr size_t nb_points = 6;
+ SmallArray<TinyVector<3>> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ point_list[0] = {+1.0000000000, +0.0000000000, +0.0000000000};
+ point_list[1] = {-1.0000000000, +0.0000000000, +0.0000000000};
+ point_list[2] = {+0.0000000000, +1.0000000000, +0.0000000000};
+ point_list[3] = {+0.0000000000, -1.0000000000, +0.0000000000};
+ point_list[4] = {+0.0000000000, +0.0000000000, +1.0000000000};
+ point_list[5] = {+0.0000000000, +0.0000000000, -1.0000000000};
+
+ weight_list[0] = 1.3333333333;
+ weight_list[1] = 1.3333333333;
+ weight_list[2] = 1.3333333333;
+ weight_list[3] = 1.3333333333;
+ weight_list[4] = 1.3333333333;
+ weight_list[5] = 1.3333333333;
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 4:
+ 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.7958224257, +0.0000000000, +0.0000000000};
+ point_list[1] = {-0.7958224257, +0.0000000000, +0.0000000000};
+ point_list[2] = {+0.0000000000, +0.7958224257, +0.0000000000};
+ point_list[3] = {+0.0000000000, -0.7958224257, +0.0000000000};
+ point_list[4] = {+0.0000000000, +0.0000000000, +0.7958224257};
+ point_list[5] = {+0.0000000000, +0.0000000000, -0.7958224257};
+ point_list[6] = {+0.7587869106, +0.7587869106, +0.7587869106};
+ point_list[7] = {+0.7587869106, -0.7587869106, +0.7587869106};
+ point_list[8] = {+0.7587869106, +0.7587869106, -0.7587869106};
+ point_list[9] = {+0.7587869106, -0.7587869106, -0.7587869106};
+ point_list[10] = {-0.7587869106, +0.7587869106, +0.7587869106};
+ point_list[11] = {-0.7587869106, -0.7587869106, +0.7587869106};
+ point_list[12] = {-0.7587869106, +0.7587869106, -0.7587869106};
+ point_list[13] = {-0.7587869106, -0.7587869106, -0.7587869106};
+
+ weight_list[0] = 0.8864265927;
+ weight_list[1] = 0.8864265927;
+ weight_list[2] = 0.8864265927;
+ weight_list[3] = 0.8864265927;
+ weight_list[4] = 0.8864265927;
+ weight_list[5] = 0.8864265927;
+ weight_list[6] = 0.3351800554;
+ weight_list[7] = 0.3351800554;
+ weight_list[8] = 0.3351800554;
+ weight_list[9] = 0.3351800554;
+ weight_list[10] = 0.3351800554;
+ weight_list[11] = 0.3351800554;
+ weight_list[12] = 0.3351800554;
+ weight_list[13] = 0.3351800554;
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 6:
+ case 7: {
+ constexpr size_t nb_points = 27;
+ SmallArray<TinyVector<3>> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ point_list[0] = {+0.0000000000, +0.0000000000, +0.0000000000};
+ point_list[1] = {+0.8484180014, +0.0000000000, +0.0000000000};
+ point_list[2] = {-0.8484180014, +0.0000000000, +0.0000000000};
+ point_list[3] = {+0.0000000000, +0.8484180014, +0.0000000000};
+ point_list[4] = {+0.0000000000, -0.8484180014, +0.0000000000};
+ point_list[5] = {+0.0000000000, +0.0000000000, +0.8484180014};
+ point_list[6] = {+0.0000000000, +0.0000000000, -0.8484180014};
+ point_list[7] = {+0.6528164721, +0.6528164721, +0.6528164721};
+ point_list[8] = {+0.6528164721, -0.6528164721, +0.6528164721};
+ point_list[9] = {+0.6528164721, +0.6528164721, -0.6528164721};
+ point_list[10] = {+0.6528164721, -0.6528164721, -0.6528164721};
+ point_list[11] = {-0.6528164721, +0.6528164721, +0.6528164721};
+ point_list[12] = {-0.6528164721, -0.6528164721, +0.6528164721};
+ point_list[13] = {-0.6528164721, +0.6528164721, -0.6528164721};
+ point_list[14] = {-0.6528164721, -0.6528164721, -0.6528164721};
+ point_list[15] = {+0.0000000000, +1.1064128986, +1.1064128986};
+ point_list[16] = {+0.0000000000, -1.1064128986, +1.1064128986};
+ point_list[17] = {+0.0000000000, +1.1064128986, -1.1064128986};
+ point_list[18] = {+0.0000000000, -1.1064128986, -1.1064128986};
+ point_list[19] = {+1.1064128986, +0.0000000000, +1.1064128986};
+ point_list[20] = {-1.1064128986, +0.0000000000, +1.1064128986};
+ point_list[21] = {+1.1064128986, +0.0000000000, -1.1064128986};
+ point_list[22] = {-1.1064128986, +0.0000000000, -1.1064128986};
+ point_list[23] = {+1.1064128986, +1.1064128986, +0.0000000000};
+ point_list[24] = {-1.1064128986, +1.1064128986, +0.0000000000};
+ point_list[25] = {+1.1064128986, -1.1064128986, +0.0000000000};
+ point_list[26] = {-1.1064128986, -1.1064128986, +0.0000000000};
+
+ weight_list[0] = 0.7880734827;
+ weight_list[1] = 0.4993690023;
+ weight_list[2] = 0.4993690023;
+ weight_list[3] = 0.4993690023;
+ weight_list[4] = 0.4993690023;
+ weight_list[5] = 0.4993690023;
+ weight_list[6] = 0.4993690023;
+ weight_list[7] = 0.4785084494;
+ weight_list[8] = 0.4785084494;
+ weight_list[9] = 0.4785084494;
+ weight_list[10] = 0.4785084494;
+ weight_list[11] = 0.4785084494;
+ weight_list[12] = 0.4785084494;
+ weight_list[13] = 0.4785084494;
+ weight_list[14] = 0.4785084494;
+ weight_list[15] = 0.0323037423;
+ weight_list[16] = 0.0323037423;
+ weight_list[17] = 0.0323037423;
+ weight_list[18] = 0.0323037423;
+ weight_list[19] = 0.0323037423;
+ weight_list[20] = 0.0323037423;
+ weight_list[21] = 0.0323037423;
+ weight_list[22] = 0.0323037423;
+ weight_list[23] = 0.0323037423;
+ weight_list[24] = 0.0323037423;
+ weight_list[25] = 0.0323037423;
+ weight_list[26] = 0.0323037423;
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ default: {
+ throw NormalError("Gauss quadrature formulae handle degrees up to 7 on hexahedra");
+ }
+ }
+}
diff --git a/src/analysis/EconomicalGaussQuadrature.hpp b/src/analysis/EconomicalGaussQuadrature.hpp
new file mode 100644
index 000000000..0862d0a66
--- /dev/null
+++ b/src/analysis/EconomicalGaussQuadrature.hpp
@@ -0,0 +1,31 @@
+#ifndef ECONOMICAL_GAUSS_QUADRATURE_HPP
+#define ECONOMICAL_GAUSS_QUADRATURE_HPP
+
+#include <analysis/QuadratureFormula.hpp>
+
+/**
+ * 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]
+ */
+template <size_t Dimension>
+class EconomicalGaussQuadrature final : public QuadratureForumla<Dimension>
+{
+ private:
+ void _buildPointAndWeightLists();
+
+ public:
+ EconomicalGaussQuadrature(EconomicalGaussQuadrature&&) = default;
+ EconomicalGaussQuadrature(const EconomicalGaussQuadrature&) = default;
+
+ explicit EconomicalGaussQuadrature(const size_t degree) : QuadratureForumla<Dimension>(degree)
+ {
+ this->_buildPointAndWeightLists();
+ }
+
+ EconomicalGaussQuadrature() = delete;
+ ~EconomicalGaussQuadrature() = default;
+};
+
+#endif // ECONOMICAL_GAUSS_QUADRATURE_HPP
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index b188acce3..a55587ae5 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -68,6 +68,7 @@ add_executable (unit_tests
test_DiscreteFunctionType.cpp
test_DiscreteFunctionUtils.cpp
test_DoWhileProcessor.cpp
+ test_EconomicalGaussQuadrature.cpp
test_EigenvalueSolver.cpp
test_EmbeddedData.cpp
test_EmbeddedIDiscreteFunctionUtils.cpp
diff --git a/tests/test_EconomicalGaussQuadrature.cpp b/tests/test_EconomicalGaussQuadrature.cpp
new file mode 100644
index 000000000..1da84db8d
--- /dev/null
+++ b/tests/test_EconomicalGaussQuadrature.cpp
@@ -0,0 +1,877 @@
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <analysis/EconomicalGaussQuadrature.hpp>
+#include <geometry/AffineTransformation.hpp>
+#include <utils/Exceptions.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("EconomicalGaussQuadrature", "[analysis]")
+{
+ auto same_formulas = [](auto quadrature_formula0, auto quadrature_formula1) {
+ auto point_list0 = quadrature_formula0.pointList();
+ auto weight_list0 = quadrature_formula0.weightList();
+
+ auto point_list1 = quadrature_formula1.pointList();
+ auto weight_list1 = quadrature_formula1.weightList();
+
+ bool is_same = true;
+
+ for (size_t i = 0; i < point_list0.size(); ++i) {
+ if (point_list0[i] != point_list1[i]) {
+ return false;
+ }
+ }
+ for (size_t i = 0; i < weight_list0.size(); ++i) {
+ if (weight_list0[i] != weight_list1[i]) {
+ return false;
+ }
+ }
+
+ return is_same;
+ };
+
+ SECTION("1D")
+ {
+ auto is_symmetric_formula = [](auto quadrature_formula) {
+ auto point_list = quadrature_formula.pointList();
+ auto weight_list = quadrature_formula.weightList();
+
+ bool is_symmetric = true;
+
+ const size_t middle_index = point_list.size() / 2;
+ for (size_t i = 0; i <= middle_index; ++i) {
+ if (point_list[i] != -point_list[point_list.size() - 1 - i]) {
+ return false;
+ }
+ }
+ for (size_t i = 0; i <= middle_index; ++i) {
+ if (weight_list[i] != weight_list[point_list.size() - 1 - i]) {
+ return false;
+ }
+ }
+
+ return is_symmetric;
+ };
+
+ auto integrate = [](auto f, auto quadrature_formula, const double a, const double b) {
+ auto point_list = quadrature_formula.pointList();
+ auto weight_list = quadrature_formula.weightList();
+
+ double alpha = 0.5 * (b - a);
+ double beta = 0.5 * (a + b);
+
+ auto x = [&alpha, &beta](auto x_hat) { return alpha * x_hat + beta; };
+
+ auto value = weight_list[0] * f(x(point_list[0]));
+ for (size_t i = 1; i < weight_list.size(); ++i) {
+ value += weight_list[i] * f(x(point_list[i]));
+ }
+
+ return alpha * value;
+ };
+
+ auto get_order = [&integrate](auto f, auto quadrature_formula, const double a, const double b,
+ const double exact_value) {
+ double int_ab = integrate(f, quadrature_formula, a, b);
+ double int_first_half = integrate(f, quadrature_formula, a, 0.5 * (a + b));
+ double int_second_half = integrate(f, quadrature_formula, 0.5 * (a + b), b);
+
+ return -std::log((int_first_half + int_second_half - exact_value) / (int_ab - exact_value)) / std::log(2);
+ };
+
+ auto p0 = [](const TinyVector<1>&) { return 1; };
+ auto p1 = [&p0](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 2 * x + p0(X);
+ };
+ auto p2 = [&p1](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 3 * x * x + p1(X);
+ };
+ auto p3 = [&p2](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 4 * std::pow(x, 3) + p2(X);
+ };
+ auto p4 = [&p3](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 5 * std::pow(x, 4) + p3(X);
+ };
+ auto p5 = [&p4](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 6 * std::pow(x, 5) + p4(X);
+ };
+ auto p6 = [&p5](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 7 * std::pow(x, 6) + p5(X);
+ };
+ auto p7 = [&p6](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 8 * std::pow(x, 7) + p6(X);
+ };
+ auto p8 = [&p7](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 9 * std::pow(x, 8) + p7(X);
+ };
+ auto p9 = [&p8](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 10 * std::pow(x, 9) + p8(X);
+ };
+ auto p10 = [&p9](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 11 * std::pow(x, 10) + p9(X);
+ };
+ auto p11 = [&p10](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 12 * std::pow(x, 11) + p10(X);
+ };
+ auto p12 = [&p11](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 13 * std::pow(x, 12) + p11(X);
+ };
+ auto p13 = [&p12](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 14 * std::pow(x, 13) + p12(X);
+ };
+ auto p14 = [&p13](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 15 * std::pow(x, 14) + p13(X);
+ };
+ auto p15 = [&p14](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 16 * std::pow(x, 15) + p14(X);
+ };
+ auto p16 = [&p15](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 17 * std::pow(x, 16) + p15(X);
+ };
+ auto p17 = [&p16](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 18 * std::pow(x, 17) + p16(X);
+ };
+ auto p18 = [&p17](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 19 * std::pow(x, 18) + p17(X);
+ };
+ auto p19 = [&p18](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 20 * std::pow(x, 19) + p18(X);
+ };
+ auto p20 = [&p19](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 21 * std::pow(x, 20) + p19(X);
+ };
+ auto p21 = [&p20](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 22 * std::pow(x, 21) + p20(X);
+ };
+ auto p22 = [&p21](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 23 * std::pow(x, 22) + p21(X);
+ };
+ auto p23 = [&p22](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 24 * std::pow(x, 23) + p22(X);
+ };
+ auto p24 = [&p23](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 25 * std::pow(x, 24) + p23(X);
+ };
+
+ SECTION("degree 0 and 1")
+ {
+ EconomicalGaussQuadrature<1> l0(0);
+ EconomicalGaussQuadrature<1> l1(1);
+
+ REQUIRE(same_formulas(l0, l1));
+ REQUIRE(is_symmetric_formula(l1));
+
+ REQUIRE(l1.numberOfPoints() == 1);
+
+ REQUIRE(integrate(p0, l1, 0.5, 1) == Catch::Approx(0.5));
+ REQUIRE(integrate(p1, l1, 0, 1) == Catch::Approx(2));
+ REQUIRE(integrate(p2, l1, 0, 1) != Catch::Approx(3));
+
+ REQUIRE(get_order(p2, l1, -1, 1, 4) == Catch::Approx(2));
+ }
+
+ SECTION("degree 2 and 3")
+ {
+ EconomicalGaussQuadrature<1> l2(2);
+ EconomicalGaussQuadrature<1> l3(3);
+
+ REQUIRE(same_formulas(l2, l3));
+ REQUIRE(is_symmetric_formula(l3));
+
+ REQUIRE(l3.numberOfPoints() == 2);
+
+ REQUIRE(integrate(p0, l3, 0.5, 1) == Catch::Approx(0.5));
+ REQUIRE(integrate(p1, l3, 0, 1) == Catch::Approx(2));
+ REQUIRE(integrate(p2, l3, 0, 1) == Catch::Approx(3));
+ REQUIRE(integrate(p3, l3, 0, 1) == Catch::Approx(4));
+ REQUIRE(integrate(p4, l3, 0, 1) != Catch::Approx(5));
+
+ REQUIRE(get_order(p4, l3, -1, 1, 6) == Catch::Approx(4));
+ }
+
+ SECTION("degree 4 and 5")
+ {
+ EconomicalGaussQuadrature<1> l4(4);
+ EconomicalGaussQuadrature<1> l5(5);
+
+ REQUIRE(same_formulas(l4, l5));
+ REQUIRE(is_symmetric_formula(l5));
+
+ REQUIRE(l5.numberOfPoints() == 3);
+
+ REQUIRE(integrate(p0, l5, 0.5, 1) == Catch::Approx(0.5));
+ REQUIRE(integrate(p1, l5, 0, 1) == Catch::Approx(2));
+ REQUIRE(integrate(p2, l5, 0, 1) == Catch::Approx(3));
+ REQUIRE(integrate(p3, l5, 0, 1) == Catch::Approx(4));
+ REQUIRE(integrate(p4, l5, 0, 1) == Catch::Approx(5));
+ REQUIRE(integrate(p5, l5, 0, 1) == Catch::Approx(6));
+ REQUIRE(integrate(p6, l5, 0, 1) != Catch::Approx(7));
+
+ REQUIRE(get_order(p6, l5, -1, 1, 8) == Catch::Approx(6));
+ }
+
+ SECTION("degree 6 and 7")
+ {
+ EconomicalGaussQuadrature<1> l6(6);
+ EconomicalGaussQuadrature<1> l7(7);
+
+ REQUIRE(same_formulas(l6, l7));
+ REQUIRE(is_symmetric_formula(l7));
+
+ REQUIRE(l7.numberOfPoints() == 4);
+
+ REQUIRE(integrate(p0, l7, 0.5, 1) == Catch::Approx(0.5));
+ REQUIRE(integrate(p1, l7, 0, 1) == Catch::Approx(2));
+ REQUIRE(integrate(p2, l7, 0, 1) == Catch::Approx(3));
+ REQUIRE(integrate(p3, l7, 0, 1) == Catch::Approx(4));
+ REQUIRE(integrate(p4, l7, 0, 1) == Catch::Approx(5));
+ REQUIRE(integrate(p5, l7, 0, 1) == Catch::Approx(6));
+ REQUIRE(integrate(p6, l7, 0, 1) == Catch::Approx(7));
+ REQUIRE(integrate(p7, l7, 0, 1) == Catch::Approx(8));
+ REQUIRE(integrate(p8, l7, 0, 1) != Catch::Approx(9));
+
+ REQUIRE(get_order(p8, l7, -1, 1, 10) == Catch::Approx(8));
+ }
+
+ SECTION("degree 8 and 9")
+ {
+ EconomicalGaussQuadrature<1> l8(8);
+ EconomicalGaussQuadrature<1> l9(9);
+
+ REQUIRE(same_formulas(l8, l9));
+ REQUIRE(is_symmetric_formula(l9));
+
+ REQUIRE(l9.numberOfPoints() == 5);
+
+ REQUIRE(integrate(p0, l9, 0.5, 1) == Catch::Approx(0.5));
+ REQUIRE(integrate(p1, l9, 0, 1) == Catch::Approx(2));
+ REQUIRE(integrate(p2, l9, 0, 1) == Catch::Approx(3));
+ REQUIRE(integrate(p3, l9, 0, 1) == Catch::Approx(4));
+ REQUIRE(integrate(p4, l9, 0, 1) == Catch::Approx(5));
+ REQUIRE(integrate(p5, l9, 0, 1) == Catch::Approx(6));
+ REQUIRE(integrate(p6, l9, 0, 1) == Catch::Approx(7));
+ REQUIRE(integrate(p7, l9, 0, 1) == Catch::Approx(8));
+ REQUIRE(integrate(p8, l9, 0, 1) == Catch::Approx(9));
+ REQUIRE(integrate(p9, l9, 0, 1) == Catch::Approx(10));
+ REQUIRE(integrate(p10, l9, 0, 1) != Catch::Approx(11).epsilon(1E-13));
+
+ REQUIRE(get_order(p10, l9, -1, 1, 12) == Catch::Approx(10));
+ }
+
+ SECTION("degree 10 and 11")
+ {
+ EconomicalGaussQuadrature<1> l10(10);
+ EconomicalGaussQuadrature<1> l11(11);
+
+ REQUIRE(same_formulas(l10, l11));
+ REQUIRE(is_symmetric_formula(l11));
+
+ REQUIRE(l11.numberOfPoints() == 6);
+
+ REQUIRE(integrate(p0, l11, 0.5, 1) == Catch::Approx(0.5));
+ REQUIRE(integrate(p1, l11, 0, 1) == Catch::Approx(2));
+ REQUIRE(integrate(p2, l11, 0, 1) == Catch::Approx(3));
+ REQUIRE(integrate(p3, l11, 0, 1) == Catch::Approx(4));
+ REQUIRE(integrate(p4, l11, 0, 1) == Catch::Approx(5));
+ REQUIRE(integrate(p5, l11, 0, 1) == Catch::Approx(6));
+ REQUIRE(integrate(p6, l11, 0, 1) == Catch::Approx(7));
+ REQUIRE(integrate(p7, l11, 0, 1) == Catch::Approx(8));
+ REQUIRE(integrate(p8, l11, 0, 1) == Catch::Approx(9));
+ REQUIRE(integrate(p9, l11, 0, 1) == Catch::Approx(10));
+ REQUIRE(integrate(p10, l11, 0, 1) == Catch::Approx(11));
+ REQUIRE(integrate(p11, l11, 0, 1) == Catch::Approx(12));
+ REQUIRE(integrate(p12, l11, 0, 1) != Catch::Approx(13).epsilon(1E-13));
+
+ REQUIRE(get_order(p12, l11, -1, 1, 14) == Catch::Approx(12));
+ }
+
+ SECTION("degree 12 and 13")
+ {
+ EconomicalGaussQuadrature<1> l12(12);
+ EconomicalGaussQuadrature<1> l13(13);
+
+ REQUIRE(same_formulas(l12, l13));
+ REQUIRE(is_symmetric_formula(l13));
+
+ REQUIRE(l13.numberOfPoints() == 7);
+
+ REQUIRE(integrate(p0, l13, 0.5, 1) == Catch::Approx(0.5));
+ REQUIRE(integrate(p1, l13, 0, 1) == Catch::Approx(2));
+ REQUIRE(integrate(p2, l13, 0, 1) == Catch::Approx(3));
+ REQUIRE(integrate(p3, l13, 0, 1) == Catch::Approx(4));
+ REQUIRE(integrate(p4, l13, 0, 1) == Catch::Approx(5));
+ REQUIRE(integrate(p5, l13, 0, 1) == Catch::Approx(6));
+ REQUIRE(integrate(p6, l13, 0, 1) == Catch::Approx(7));
+ REQUIRE(integrate(p7, l13, 0, 1) == Catch::Approx(8));
+ REQUIRE(integrate(p8, l13, 0, 1) == Catch::Approx(9));
+ REQUIRE(integrate(p9, l13, 0, 1) == Catch::Approx(10));
+ REQUIRE(integrate(p10, l13, 0, 1) == Catch::Approx(11));
+ REQUIRE(integrate(p11, l13, 0, 1) == Catch::Approx(12));
+ REQUIRE(integrate(p12, l13, 0, 1) == Catch::Approx(13));
+ REQUIRE(integrate(p13, l13, 0, 1) == Catch::Approx(14));
+ REQUIRE(integrate(p14, l13, 0, 1) != Catch::Approx(15).epsilon(1E-13));
+
+ REQUIRE(get_order(p14, l13, -1, 1, 16) == Catch::Approx(14));
+ }
+
+ SECTION("degree 14 and 15")
+ {
+ EconomicalGaussQuadrature<1> l14(14);
+ EconomicalGaussQuadrature<1> l15(15);
+
+ REQUIRE(same_formulas(l14, l15));
+ REQUIRE(is_symmetric_formula(l15));
+
+ REQUIRE(l15.numberOfPoints() == 8);
+
+ REQUIRE(integrate(p0, l15, 0.5, 1) == Catch::Approx(0.5));
+ REQUIRE(integrate(p1, l15, 0, 1) == Catch::Approx(2));
+ REQUIRE(integrate(p2, l15, 0, 1) == Catch::Approx(3));
+ REQUIRE(integrate(p3, l15, 0, 1) == Catch::Approx(4));
+ REQUIRE(integrate(p4, l15, 0, 1) == Catch::Approx(5));
+ REQUIRE(integrate(p5, l15, 0, 1) == Catch::Approx(6));
+ REQUIRE(integrate(p6, l15, 0, 1) == Catch::Approx(7));
+ REQUIRE(integrate(p7, l15, 0, 1) == Catch::Approx(8));
+ REQUIRE(integrate(p8, l15, 0, 1) == Catch::Approx(9));
+ REQUIRE(integrate(p9, l15, 0, 1) == Catch::Approx(10));
+ REQUIRE(integrate(p10, l15, 0, 1) == Catch::Approx(11));
+ REQUIRE(integrate(p11, l15, 0, 1) == Catch::Approx(12));
+ REQUIRE(integrate(p12, l15, 0, 1) == Catch::Approx(13));
+ REQUIRE(integrate(p13, l15, 0, 1) == Catch::Approx(14));
+ REQUIRE(integrate(p14, l15, 0, 1) == Catch::Approx(15));
+ REQUIRE(integrate(p15, l15, 0, 1) == Catch::Approx(16));
+ REQUIRE(integrate(p16, l15, 0, 1) != Catch::Approx(17).epsilon(1E-13));
+
+ REQUIRE(get_order(p16, l15, -1, 1, 18) == Catch::Approx(16));
+ }
+
+ SECTION("degree 16 and 17")
+ {
+ EconomicalGaussQuadrature<1> l16(16);
+ EconomicalGaussQuadrature<1> l17(17);
+
+ REQUIRE(same_formulas(l16, l17));
+ REQUIRE(is_symmetric_formula(l17));
+
+ REQUIRE(l17.numberOfPoints() == 9);
+
+ REQUIRE(integrate(p0, l17, 0.5, 1) == Catch::Approx(0.5));
+ REQUIRE(integrate(p1, l17, 0, 1) == Catch::Approx(2));
+ REQUIRE(integrate(p2, l17, 0, 1) == Catch::Approx(3));
+ REQUIRE(integrate(p3, l17, 0, 1) == Catch::Approx(4));
+ REQUIRE(integrate(p4, l17, 0, 1) == Catch::Approx(5));
+ REQUIRE(integrate(p5, l17, 0, 1) == Catch::Approx(6));
+ REQUIRE(integrate(p6, l17, 0, 1) == Catch::Approx(7));
+ REQUIRE(integrate(p7, l17, 0, 1) == Catch::Approx(8));
+ REQUIRE(integrate(p8, l17, 0, 1) == Catch::Approx(9));
+ REQUIRE(integrate(p9, l17, 0, 1) == Catch::Approx(10));
+ REQUIRE(integrate(p10, l17, 0, 1) == Catch::Approx(11));
+ REQUIRE(integrate(p11, l17, 0, 1) == Catch::Approx(12));
+ REQUIRE(integrate(p12, l17, 0, 1) == Catch::Approx(13));
+ REQUIRE(integrate(p13, l17, 0, 1) == Catch::Approx(14));
+ REQUIRE(integrate(p14, l17, 0, 1) == Catch::Approx(15));
+ REQUIRE(integrate(p15, l17, 0, 1) == Catch::Approx(16));
+ REQUIRE(integrate(p16, l17, 0, 1) == Catch::Approx(17));
+ REQUIRE(integrate(p17, l17, 0, 1) == Catch::Approx(18));
+ REQUIRE(integrate(p18, l17, 0, 1) != Catch::Approx(19).epsilon(1E-13));
+
+ REQUIRE(get_order(p18, l17, -1, 1, 20) == Catch::Approx(18));
+ }
+
+ SECTION("degree 18 and 19")
+ {
+ EconomicalGaussQuadrature<1> l18(18);
+ EconomicalGaussQuadrature<1> l19(19);
+
+ REQUIRE(same_formulas(l18, l19));
+ REQUIRE(is_symmetric_formula(l19));
+
+ REQUIRE(l19.numberOfPoints() == 10);
+
+ REQUIRE(integrate(p0, l19, 0.5, 1) == Catch::Approx(0.5));
+ REQUIRE(integrate(p1, l19, 0, 1) == Catch::Approx(2));
+ REQUIRE(integrate(p2, l19, 0, 1) == Catch::Approx(3));
+ REQUIRE(integrate(p3, l19, 0, 1) == Catch::Approx(4));
+ REQUIRE(integrate(p4, l19, 0, 1) == Catch::Approx(5));
+ REQUIRE(integrate(p5, l19, 0, 1) == Catch::Approx(6));
+ REQUIRE(integrate(p6, l19, 0, 1) == Catch::Approx(7));
+ REQUIRE(integrate(p7, l19, 0, 1) == Catch::Approx(8));
+ REQUIRE(integrate(p8, l19, 0, 1) == Catch::Approx(9));
+ REQUIRE(integrate(p9, l19, 0, 1) == Catch::Approx(10));
+ REQUIRE(integrate(p10, l19, 0, 1) == Catch::Approx(11));
+ REQUIRE(integrate(p11, l19, 0, 1) == Catch::Approx(12));
+ REQUIRE(integrate(p12, l19, 0, 1) == Catch::Approx(13));
+ REQUIRE(integrate(p13, l19, 0, 1) == Catch::Approx(14));
+ REQUIRE(integrate(p14, l19, 0, 1) == Catch::Approx(15));
+ REQUIRE(integrate(p15, l19, 0, 1) == Catch::Approx(16));
+ REQUIRE(integrate(p16, l19, 0, 1) == Catch::Approx(17));
+ REQUIRE(integrate(p17, l19, 0, 1) == Catch::Approx(18));
+ REQUIRE(integrate(p18, l19, 0, 1) == Catch::Approx(19));
+ REQUIRE(integrate(p19, l19, 0, 1) == Catch::Approx(20));
+ REQUIRE(integrate(p20, l19, 0, 1) != Catch::Approx(21).epsilon(1E-13));
+
+ REQUIRE(get_order(p20, l19, -1, 1, 22) == Catch::Approx(20));
+ }
+
+ SECTION("degree 20 and 21")
+ {
+ EconomicalGaussQuadrature<1> l20(20);
+ EconomicalGaussQuadrature<1> l21(21);
+
+ REQUIRE(same_formulas(l20, l21));
+ REQUIRE(is_symmetric_formula(l21));
+
+ REQUIRE(l21.numberOfPoints() == 11);
+
+ REQUIRE(integrate(p0, l21, 0.5, 1) == Catch::Approx(0.5));
+ REQUIRE(integrate(p1, l21, 0, 1) == Catch::Approx(2));
+ REQUIRE(integrate(p2, l21, 0, 1) == Catch::Approx(3));
+ REQUIRE(integrate(p3, l21, 0, 1) == Catch::Approx(4));
+ REQUIRE(integrate(p4, l21, 0, 1) == Catch::Approx(5));
+ REQUIRE(integrate(p5, l21, 0, 1) == Catch::Approx(6));
+ REQUIRE(integrate(p6, l21, 0, 1) == Catch::Approx(7));
+ REQUIRE(integrate(p7, l21, 0, 1) == Catch::Approx(8));
+ REQUIRE(integrate(p8, l21, 0, 1) == Catch::Approx(9));
+ REQUIRE(integrate(p9, l21, 0, 1) == Catch::Approx(10));
+ REQUIRE(integrate(p10, l21, 0, 1) == Catch::Approx(11));
+ REQUIRE(integrate(p11, l21, 0, 1) == Catch::Approx(12));
+ REQUIRE(integrate(p12, l21, 0, 1) == Catch::Approx(13));
+ REQUIRE(integrate(p13, l21, 0, 1) == Catch::Approx(14));
+ REQUIRE(integrate(p14, l21, 0, 1) == Catch::Approx(15));
+ REQUIRE(integrate(p15, l21, 0, 1) == Catch::Approx(16));
+ REQUIRE(integrate(p16, l21, 0, 1) == Catch::Approx(17));
+ REQUIRE(integrate(p17, l21, 0, 1) == Catch::Approx(18));
+ REQUIRE(integrate(p18, l21, 0, 1) == Catch::Approx(19));
+ REQUIRE(integrate(p19, l21, 0, 1) == Catch::Approx(20));
+ REQUIRE(integrate(p20, l21, 0, 1) == Catch::Approx(21));
+ REQUIRE(integrate(p21, l21, 0, 1) == Catch::Approx(22));
+ REQUIRE(integrate(p22, l21, 0, 1) != Catch::Approx(23).epsilon(1E-14));
+
+ REQUIRE(get_order(p22, l21, -1, 1, 24) == Catch::Approx(22).margin(0.02));
+ }
+
+ SECTION("degree 22 and 23")
+ {
+ EconomicalGaussQuadrature<1> l22(22);
+ EconomicalGaussQuadrature<1> l23(23);
+
+ REQUIRE(same_formulas(l22, l23));
+ REQUIRE(is_symmetric_formula(l23));
+
+ REQUIRE(l23.numberOfPoints() == 12);
+
+ REQUIRE(integrate(p0, l23, 0.5, 1) == Catch::Approx(0.5));
+ REQUIRE(integrate(p1, l23, 0, 1) == Catch::Approx(2));
+ REQUIRE(integrate(p2, l23, 0, 1) == Catch::Approx(3));
+ REQUIRE(integrate(p3, l23, 0, 1) == Catch::Approx(4));
+ REQUIRE(integrate(p4, l23, 0, 1) == Catch::Approx(5));
+ REQUIRE(integrate(p5, l23, 0, 1) == Catch::Approx(6));
+ REQUIRE(integrate(p6, l23, 0, 1) == Catch::Approx(7));
+ REQUIRE(integrate(p7, l23, 0, 1) == Catch::Approx(8));
+ REQUIRE(integrate(p8, l23, 0, 1) == Catch::Approx(9));
+ REQUIRE(integrate(p9, l23, 0, 1) == Catch::Approx(10));
+ REQUIRE(integrate(p10, l23, 0, 1) == Catch::Approx(11));
+ REQUIRE(integrate(p11, l23, 0, 1) == Catch::Approx(12));
+ REQUIRE(integrate(p12, l23, 0, 1) == Catch::Approx(13));
+ REQUIRE(integrate(p13, l23, 0, 1) == Catch::Approx(14));
+ REQUIRE(integrate(p14, l23, 0, 1) == Catch::Approx(15));
+ REQUIRE(integrate(p15, l23, 0, 1) == Catch::Approx(16));
+ REQUIRE(integrate(p16, l23, 0, 1) == Catch::Approx(17));
+ REQUIRE(integrate(p17, l23, 0, 1) == Catch::Approx(18));
+ REQUIRE(integrate(p18, l23, 0, 1) == Catch::Approx(19));
+ REQUIRE(integrate(p19, l23, 0, 1) == Catch::Approx(20));
+ REQUIRE(integrate(p20, l23, 0, 1) == Catch::Approx(21));
+ REQUIRE(integrate(p21, l23, 0, 1) == Catch::Approx(22));
+ REQUIRE(integrate(p22, l23, 0, 1) == Catch::Approx(23));
+ REQUIRE(integrate(p23, l23, 0, 1) == Catch::Approx(24));
+ REQUIRE(integrate(p24, l23, 0, 1) != Catch::Approx(25).epsilon(1E-15));
+
+ REQUIRE(get_order(p24, l23, -1, 1, 26) == Catch::Approx(24).margin(0.05));
+ }
+
+ SECTION("not implemented formulae")
+ {
+ REQUIRE_THROWS_WITH(EconomicalGaussQuadrature<1>(24),
+ "error: Gauss-Legendre quadrature formulae handle degrees up to 23");
+ }
+ }
+
+ SECTION("2D")
+ {
+ 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_rectangle = [](auto f, auto quadrature_formula, const std::array<TinyVector<2>, 3>& triangle) {
+ AffineTransformation<2> t{triangle};
+
+ 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_rectangle](auto f, auto quadrature_formula, const double exact_value) {
+ using R2 = TinyVector<2>;
+ const double int_K_hat = integrate(f, quadrature_formula);
+ const double int_refined //
+ = integrate_on_rectangle(f, quadrature_formula, {R2{-1, -1}, R2{0, -1}, R2{-1, 0}}) +
+ integrate_on_rectangle(f, quadrature_formula, {R2{0, -1}, R2{1, -1}, R2{0, 0}}) +
+ integrate_on_rectangle(f, quadrature_formula, {R2{-1, 0}, R2{0, 0}, R2{-1, 1}}) +
+ integrate_on_rectangle(f, quadrature_formula, {R2{0, 0}, R2{1, 0}, R2{0, 1}});
+
+ return -std::log((int_refined - exact_value) / (int_K_hat - exact_value)) / std::log(2);
+ };
+
+ auto p0 = [](const TinyVector<2>&) { return 2; };
+ auto p1 = [](const TinyVector<2>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ return 2 * x + 3 * y + 1;
+ };
+ auto p2 = [&p1](const TinyVector<2>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ return p1(X) * (2.5 * x - 3 * y + 3);
+ };
+ auto p3 = [&p2](const TinyVector<2>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ return p2(X) * (-1.5 * x + 3 * y - 3);
+ };
+ auto p4 = [&p3](const TinyVector<2>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ return p3(X) * (x + y + 1);
+ };
+ auto p5 = [&p4](const TinyVector<2>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ return p4(X) * (-0.2 * x - 1.3 * y - 0.7);
+ };
+ auto p6 = [&p5](const TinyVector<2>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ return p5(X) * (3 * x - 2 * y + 3);
+ };
+ auto p7 = [&p6](const TinyVector<2>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ return p6(X) * (-2 * x + 4 * y - 7);
+ };
+ auto p8 = [&p7](const TinyVector<2>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ return p7(X) * (2 * x - 3 * y - 3);
+ };
+
+ SECTION("degree 0 and 1")
+ {
+ EconomicalGaussQuadrature<2> l1(1);
+
+ REQUIRE(l1.numberOfPoints() == 1);
+
+ REQUIRE(integrate(p0, l1) == Catch::Approx(8));
+ REQUIRE(integrate(p1, l1) == Catch::Approx(4));
+ REQUIRE(integrate(p2, l1) != Catch::Approx(20. / 3));
+
+ REQUIRE(get_order(p2, l1, 20. / 3) == Catch::Approx(2));
+ }
+
+ SECTION("degree 2 and 3")
+ {
+ EconomicalGaussQuadrature<2> l2(2);
+ EconomicalGaussQuadrature<2> l3(3);
+
+ REQUIRE(same_formulas(l2, l3));
+
+ REQUIRE(l3.numberOfPoints() == 4);
+
+ REQUIRE(integrate(p0, l3) == Catch::Approx(8));
+ REQUIRE(integrate(p1, l3) == Catch::Approx(4));
+ REQUIRE(integrate(p2, l3) == Catch::Approx(20. / 3));
+ REQUIRE(integrate(p3, l3) == Catch::Approx(-13));
+ REQUIRE(integrate(p4, l3) != Catch::Approx(-1184. / 15));
+
+ REQUIRE(get_order(p4, l3, -1184. / 15) == Catch::Approx(4));
+ }
+
+ SECTION("degree 4 and 5")
+ {
+ EconomicalGaussQuadrature<2> l4(4);
+ EconomicalGaussQuadrature<2> l5(5);
+
+ REQUIRE(same_formulas(l4, l5));
+
+ REQUIRE(l5.numberOfPoints() == 8);
+
+ REQUIRE(integrate(p0, l5) == Catch::Approx(8));
+ REQUIRE(integrate(p1, l5) == Catch::Approx(4));
+ REQUIRE(integrate(p2, l5) == Catch::Approx(20. / 3));
+ REQUIRE(integrate(p3, l5) == Catch::Approx(-13));
+ REQUIRE(integrate(p4, l5) == Catch::Approx(-1184. / 15));
+ REQUIRE(integrate(p5, l5) == Catch::Approx(1971. / 25));
+ REQUIRE(integrate(p6, l5) != Catch::Approx(60441. / 175));
+
+ REQUIRE(get_order(p6, l5, 60441. / 175) == Catch::Approx(6));
+ }
+
+ SECTION("degree 6 and 7")
+ {
+ EconomicalGaussQuadrature<2> l6(6);
+ EconomicalGaussQuadrature<2> l7(7);
+
+ REQUIRE(same_formulas(l6, l7));
+
+ REQUIRE(l7.numberOfPoints() == 12);
+
+ REQUIRE(integrate(p0, l7) == Catch::Approx(8));
+ REQUIRE(integrate(p1, l7) == Catch::Approx(4));
+ REQUIRE(integrate(p2, l7) == Catch::Approx(20. / 3));
+ REQUIRE(integrate(p3, l7) == Catch::Approx(-13));
+ REQUIRE(integrate(p4, l7) == Catch::Approx(-1184. / 15));
+ REQUIRE(integrate(p5, l7) == Catch::Approx(1971. / 25));
+ REQUIRE(integrate(p6, l7) == Catch::Approx(60441. / 175));
+ REQUIRE(integrate(p7, l7) == Catch::Approx(-1307119. / 525));
+ REQUIRE(integrate(p8, l7) != Catch::Approx(957697. / 175));
+
+ REQUIRE(get_order(p8, l7, 957697. / 175) == Catch::Approx(8));
+ }
+
+ SECTION("not implemented formulae")
+ {
+ REQUIRE_THROWS_WITH(EconomicalGaussQuadrature<2>(8),
+ "error: Gauss quadrature formulae handle degrees up to 7 on quadrangles");
+ }
+ }
+
+ 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_brick = [](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_brick](auto f, auto quadrature_formula, const double exact_value) {
+ using R3 = TinyVector<3>;
+
+ const double int_K_hat = integrate(f, quadrature_formula);
+ const double int_refined //
+ = integrate_on_brick(f, quadrature_formula, {R3{-1, -1, -1}, R3{+0, -1, -1}, R3{-1, +0, -1}, R3{-1, -1, +0}}) +
+ integrate_on_brick(f, quadrature_formula, {R3{+0, -1, -1}, R3{+1, -1, -1}, R3{+0, +0, -1}, R3{+0, -1, +0}}) +
+ integrate_on_brick(f, quadrature_formula, {R3{-1, +0, -1}, R3{+0, +0, -1}, R3{-1, +1, -1}, R3{-1, +0, +0}}) +
+ integrate_on_brick(f, quadrature_formula, {R3{+0, +0, -1}, R3{+1, +0, -1}, R3{+0, +1, -1}, R3{+0, +0, +0}}) +
+ integrate_on_brick(f, quadrature_formula, {R3{-1, -1, +0}, R3{+0, -1, +0}, R3{-1, +0, +0}, R3{-1, -1, +1}}) +
+ integrate_on_brick(f, quadrature_formula, {R3{+0, -1, +0}, R3{+1, -1, +0}, R3{+0, +0, +0}, R3{+0, -1, +1}}) +
+ integrate_on_brick(f, quadrature_formula, {R3{-1, +0, +0}, R3{+0, +0, +0}, R3{-1, +1, +0}, R3{-1, +0, +1}}) +
+ integrate_on_brick(f, quadrature_formula, {R3{+0, +0, +0}, R3{+1, +0, +0}, R3{+0, +1, +0}, R3{+0, +0, +1}});
+
+ return -std::log((int_refined - exact_value) / (int_K_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("degree 0 and 1")
+ {
+ EconomicalGaussQuadrature<3> l1(1);
+
+ REQUIRE(l1.numberOfPoints() == 1);
+
+ REQUIRE(integrate(p0, l1) == Catch::Approx(32));
+ REQUIRE(integrate(p1, l1) == Catch::Approx(-8));
+ REQUIRE(integrate(p2, l1) != Catch::Approx(-32));
+
+ REQUIRE(get_order(p2, l1, -32) == Catch::Approx(2));
+ }
+
+ SECTION("degree 2 and 3")
+ {
+ EconomicalGaussQuadrature<3> l2(2);
+ EconomicalGaussQuadrature<3> l3(2);
+
+ REQUIRE(same_formulas(l2, l3));
+
+ REQUIRE(l3.numberOfPoints() == 6);
+
+ REQUIRE(integrate(p0, l3) == Catch::Approx(32));
+ REQUIRE(integrate(p1, l3) == Catch::Approx(-8));
+ REQUIRE(integrate(p2, l3) == Catch::Approx(-32));
+ REQUIRE(integrate(p3, l3) == Catch::Approx(108));
+ REQUIRE(integrate(p4, l3) != Catch::Approx(868. / 75));
+
+ REQUIRE(get_order(p4, l3, 868. / 75) == Catch::Approx(4));
+ }
+
+ SECTION("degree 4 and 5")
+ {
+ EconomicalGaussQuadrature<3> l4(4);
+ EconomicalGaussQuadrature<3> l5(5);
+
+ REQUIRE(same_formulas(l4, l5));
+
+ REQUIRE(l5.numberOfPoints() == 14);
+
+ REQUIRE(integrate(p0, l5) == Catch::Approx(32));
+ REQUIRE(integrate(p1, l5) == Catch::Approx(-8));
+ REQUIRE(integrate(p2, l5) == Catch::Approx(-32));
+ REQUIRE(integrate(p3, l5) == Catch::Approx(108));
+ REQUIRE(integrate(p4, l5) == Catch::Approx(868. / 75));
+ REQUIRE(integrate(p5, l5) == Catch::Approx(11176. / 225));
+ REQUIRE(integrate(p6, l5) != Catch::Approx(-34781. / 430329));
+
+ REQUIRE(get_order(p6, l5, -34781. / 430329) == Catch::Approx(6));
+ }
+
+ SECTION("degree 6 and 7")
+ {
+ EconomicalGaussQuadrature<3> l6(6);
+ EconomicalGaussQuadrature<3> l7(7);
+
+ REQUIRE(same_formulas(l6, l7));
+
+ REQUIRE(l7.numberOfPoints() == 27);
+
+ REQUIRE(integrate(p0, l7) == Catch::Approx(32));
+ REQUIRE(integrate(p1, l7) == Catch::Approx(-8));
+ REQUIRE(integrate(p2, l7) == Catch::Approx(-32));
+ REQUIRE(integrate(p3, l7) == Catch::Approx(108));
+ REQUIRE(integrate(p4, l7) == Catch::Approx(868. / 75));
+ REQUIRE(integrate(p5, l7) == Catch::Approx(11176. / 225));
+ REQUIRE(integrate(p6, l7) == Catch::Approx(-34781. / 430329));
+ REQUIRE(integrate(p7, l7) == Catch::Approx(-37338109. / 4303290));
+ REQUIRE(integrate(p8, l7) != Catch::Approx(422437099. / 21516450));
+
+ REQUIRE(get_order(p8, l7, 422437099. / 21516450) == Catch::Approx(8));
+ }
+
+ SECTION("not implemented formulae")
+ {
+ REQUIRE_THROWS_WITH(EconomicalGaussQuadrature<3>(8),
+ "error: Gauss quadrature formulae handle degrees up to 7 on hexahedra");
+ }
+ }
+
+ SECTION("Access functions")
+ {
+ EconomicalGaussQuadrature<3> quadrature_formula(7);
+ auto point_list = quadrature_formula.pointList();
+ auto weight_list = quadrature_formula.weightList();
+
+ REQUIRE(point_list.size() == quadrature_formula.numberOfPoints());
+ REQUIRE(weight_list.size() == quadrature_formula.numberOfPoints());
+
+ for (size_t i = 0; i < quadrature_formula.numberOfPoints(); ++i) {
+ REQUIRE(&point_list[i] == &quadrature_formula.point(i));
+ REQUIRE(&weight_list[i] == &quadrature_formula.weight(i));
+ }
+ }
+}
--
GitLab