From cd8646f2016b78a9145b9d15a28937adb76cf035 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Fri, 8 Oct 2021 12:10:29 +0200
Subject: [PATCH] Start writing of simplicial Gauss-Lengendre quadrature
formulae
Only sets 1d formulae which is trivially set as the tensorial ones
---
src/analysis/CMakeLists.txt | 1 +
.../SimplicialGaussLegendreQuadrature.cpp | 27 +
.../SimplicialGaussLegendreQuadrature.hpp | 33 ++
tests/CMakeLists.txt | 1 +
...test_SimplicialGaussLegendreQuadrature.cpp | 524 ++++++++++++++++++
5 files changed, 586 insertions(+)
create mode 100644 src/analysis/SimplicialGaussLegendreQuadrature.cpp
create mode 100644 src/analysis/SimplicialGaussLegendreQuadrature.hpp
create mode 100644 tests/test_SimplicialGaussLegendreQuadrature.cpp
diff --git a/src/analysis/CMakeLists.txt b/src/analysis/CMakeLists.txt
index a99e6a339..833370313 100644
--- a/src/analysis/CMakeLists.txt
+++ b/src/analysis/CMakeLists.txt
@@ -2,5 +2,6 @@
add_library(
PugsAnalysis
+ SimplicialGaussLegendreQuadrature.cpp
TensorialGaussLegendreQuadrature.cpp
TensorialGaussLobattoQuadrature.cpp)
diff --git a/src/analysis/SimplicialGaussLegendreQuadrature.cpp b/src/analysis/SimplicialGaussLegendreQuadrature.cpp
new file mode 100644
index 000000000..d09e9a020
--- /dev/null
+++ b/src/analysis/SimplicialGaussLegendreQuadrature.cpp
@@ -0,0 +1,27 @@
+#include <analysis/SimplicialGaussLegendreQuadrature.hpp>
+#include <analysis/TensorialGaussLegendreQuadrature.hpp>
+#include <utils/Exceptions.hpp>
+
+template <>
+void
+SimplicialGaussLegendreQuadrature<1>::_buildPointAndWeightLists()
+{
+ TensorialGaussLegendreQuadrature<1> gauss(m_order);
+
+ m_point_list = gauss.pointList();
+ m_weight_list = gauss.weightList();
+}
+
+template <>
+void
+SimplicialGaussLegendreQuadrature<2>::_buildPointAndWeightLists()
+{
+ throw NotImplementedError("gauss legendre quadrature on 2d simplex");
+}
+
+template <>
+void
+SimplicialGaussLegendreQuadrature<3>::_buildPointAndWeightLists()
+{
+ throw NotImplementedError("gauss legendre quadrature on 3d simplex");
+}
diff --git a/src/analysis/SimplicialGaussLegendreQuadrature.hpp b/src/analysis/SimplicialGaussLegendreQuadrature.hpp
new file mode 100644
index 000000000..824eb0995
--- /dev/null
+++ b/src/analysis/SimplicialGaussLegendreQuadrature.hpp
@@ -0,0 +1,33 @@
+#ifndef SIMPLICIAL_GAUSS_LEGENDRE_QUADRATURE_HPP
+#define SIMPLICIAL_GAUSS_LEGENDRE_QUADRATURE_HPP
+
+#include <analysis/QuadratureFormula.hpp>
+
+/**
+ * Defines Gauss Legendre quadrature on the reference simplex element
+ * - in 1d, the segment \f$]-1,1[\f$
+ * - in 2d, the triangle defined by \f$(-1,-1)\f$, \f$(1,-1)\f$ and \f$(-1,1)\f$
+ * - in 3d, the tetrahedron joining \f$(-1,-1-1)\f$, \f$(1,-1,-1)\f$, \f$(-1,1,-1)\f$ and \f$(-1,1-,1)\f$
+ *
+ * \note formulae are extracted from High-order Finite Element Method [2004 - Chapman & Hall]
+ */
+template <size_t Dimension>
+class SimplicialGaussLegendreQuadrature final : public QuadratureForumla<Dimension>
+{
+ private:
+ void _buildPointAndWeightLists();
+
+ public:
+ SimplicialGaussLegendreQuadrature(SimplicialGaussLegendreQuadrature&&) = default;
+ SimplicialGaussLegendreQuadrature(const SimplicialGaussLegendreQuadrature&) = default;
+
+ explicit SimplicialGaussLegendreQuadrature(const size_t order) : QuadratureForumla<Dimension>(order)
+ {
+ this->_buildPointAndWeightLists();
+ }
+
+ SimplicialGaussLegendreQuadrature() = delete;
+ ~SimplicialGaussLegendreQuadrature() = default;
+};
+
+#endif // SIMPLICIAL_GAUSS_LEGENDRE_QUADRATURE_HPP
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 7ecbb1f77..69b85561e 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -99,6 +99,7 @@ add_executable (unit_tests
test_PugsUtils.cpp
test_Q1Transformation.cpp
test_RevisionInfo.cpp
+ test_SimplicialGaussLegendreQuadrature.cpp
test_SmallArray.cpp
test_SmallVector.cpp
test_SymbolTable.cpp
diff --git a/tests/test_SimplicialGaussLegendreQuadrature.cpp b/tests/test_SimplicialGaussLegendreQuadrature.cpp
new file mode 100644
index 000000000..80fd1c7e0
--- /dev/null
+++ b/tests/test_SimplicialGaussLegendreQuadrature.cpp
@@ -0,0 +1,524 @@
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <analysis/SimplicialGaussLegendreQuadrature.hpp>
+#include <utils/Exceptions.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("SimplicialGaussLegendreQuadrature", "[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("order 0 and 1")
+ {
+ SimplicialGaussLegendreQuadrature<1> l0(0);
+ SimplicialGaussLegendreQuadrature<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("order 2 and 3")
+ {
+ SimplicialGaussLegendreQuadrature<1> l2(2);
+ SimplicialGaussLegendreQuadrature<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("order 4 and 5")
+ {
+ SimplicialGaussLegendreQuadrature<1> l4(4);
+ SimplicialGaussLegendreQuadrature<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("order 6 and 7")
+ {
+ SimplicialGaussLegendreQuadrature<1> l6(6);
+ SimplicialGaussLegendreQuadrature<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("order 8 and 9")
+ {
+ SimplicialGaussLegendreQuadrature<1> l8(8);
+ SimplicialGaussLegendreQuadrature<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("order 10 and 11")
+ {
+ SimplicialGaussLegendreQuadrature<1> l10(10);
+ SimplicialGaussLegendreQuadrature<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("order 12 and 13")
+ {
+ SimplicialGaussLegendreQuadrature<1> l12(12);
+ SimplicialGaussLegendreQuadrature<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("order 14 and 15")
+ {
+ SimplicialGaussLegendreQuadrature<1> l14(14);
+ SimplicialGaussLegendreQuadrature<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("order 16 and 17")
+ {
+ SimplicialGaussLegendreQuadrature<1> l16(16);
+ SimplicialGaussLegendreQuadrature<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("order 18 and 19")
+ {
+ SimplicialGaussLegendreQuadrature<1> l18(18);
+ SimplicialGaussLegendreQuadrature<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("order 20 and 21")
+ {
+ SimplicialGaussLegendreQuadrature<1> l20(20);
+ SimplicialGaussLegendreQuadrature<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("order 22 and 23")
+ {
+ SimplicialGaussLegendreQuadrature<1> l22(22);
+ SimplicialGaussLegendreQuadrature<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(SimplicialGaussLegendreQuadrature<1>(24),
+ "error: Gauss-Legendre quadrature formulae handle orders up to 23");
+ }
+ }
+}
--
GitLab