From 2039cd5a0172c8b1be1de9240887ff0631b61982 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 18:31:04 +0200
Subject: [PATCH] Add Gauss Legendre formulae up to order 7 on reference
triangle
---
.../SimplicialGaussLegendreQuadrature.cpp | 175 ++++++++++++-
...test_SimplicialGaussLegendreQuadrature.cpp | 240 ++++++++++++++++--
2 files changed, 396 insertions(+), 19 deletions(-)
diff --git a/src/analysis/SimplicialGaussLegendreQuadrature.cpp b/src/analysis/SimplicialGaussLegendreQuadrature.cpp
index d09e9a020..c6143bdec 100644
--- a/src/analysis/SimplicialGaussLegendreQuadrature.cpp
+++ b/src/analysis/SimplicialGaussLegendreQuadrature.cpp
@@ -16,7 +16,180 @@ template <>
void
SimplicialGaussLegendreQuadrature<2>::_buildPointAndWeightLists()
{
- throw NotImplementedError("gauss legendre quadrature on 2d simplex");
+ switch (m_order) {
+ 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.333333333333333, -0.333333333333333};
+ weight_list[0] = 2;
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 2: {
+ constexpr size_t nb_points = 3;
+ SmallArray<TinyVector<2>> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ point_list[0] = {-0.666666666666667, -0.666666666666667};
+ point_list[1] = {-0.666666666666667, +0.333333333333333};
+ point_list[2] = {+0.333333333333333, -0.666666666666667};
+
+ weight_list[0] = 0.666666666666667;
+ weight_list[1] = 0.666666666666667;
+ weight_list[2] = 0.666666666666667;
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ 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.333333333333333, -0.333333333333333};
+ point_list[1] = {-0.600000000000000, -0.600000000000000};
+ point_list[2] = {-0.600000000000000, +0.200000000000000};
+ point_list[3] = {+0.200000000000000, -0.600000000000000};
+
+ weight_list[0] = -1.125;
+ weight_list[1] = 1.041666666666667;
+ weight_list[2] = 1.041666666666667;
+ weight_list[3] = 1.041666666666667;
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 4: {
+ constexpr size_t nb_points = 6;
+ SmallArray<TinyVector<2>> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ point_list[0] = {-0.108103018168070, -0.108103018168070};
+ point_list[1] = {-0.108103018168070, -0.783793963663860};
+ point_list[2] = {-0.783793963663860, -0.108103018168070};
+ point_list[3] = {-0.816847572980458, -0.816847572980458};
+ point_list[4] = {-0.816847572980458, +0.633695145960918};
+ point_list[5] = {+0.633695145960918, -0.816847572980458};
+
+ weight_list[0] = 0.446763179356022;
+ weight_list[1] = 0.446763179356022;
+ weight_list[2] = 0.446763179356022;
+ weight_list[3] = 0.219903487310644;
+ weight_list[4] = 0.219903487310644;
+ weight_list[5] = 0.219903487310644;
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 5: {
+ constexpr size_t nb_points = 7;
+ SmallArray<TinyVector<2>> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ point_list[0] = {-0.333333333333333, -0.333333333333333};
+ point_list[1] = {-0.059715871789770, -0.059715871789770};
+ point_list[2] = {-0.059715871789770, -0.880568256420460};
+ point_list[3] = {-0.880568256420460, -0.059715871789770};
+ point_list[4] = {-0.797426985353088, -0.797426985353088};
+ point_list[5] = {-0.797426985353088, +0.594853970706174};
+ point_list[6] = {+0.594853970706174, -0.797426985353088};
+
+ weight_list[0] = 0.45;
+ weight_list[1] = 0.264788305577012;
+ weight_list[2] = 0.264788305577012;
+ weight_list[3] = 0.264788305577012;
+ weight_list[4] = 0.251878361089654;
+ weight_list[5] = 0.251878361089654;
+ weight_list[6] = 0.251878361089654;
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 6: {
+ constexpr size_t nb_points = 12;
+ SmallArray<TinyVector<2>> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ point_list[0] = {-0.501426509658180, -0.501426509658180};
+ point_list[1] = {-0.501426509658180, +0.002853019316358};
+ point_list[2] = {+0.002853019316358, -0.501426509658180};
+ point_list[3] = {-0.873821971016996, -0.873821971016996};
+ point_list[4] = {-0.873821971016996, +0.747643942033992};
+ point_list[5] = {+0.747643942033992, -0.873821971016996};
+ point_list[6] = {-0.379295097932432, +0.273004998242798};
+ point_list[7] = {+0.273004998242798, -0.893709900310366};
+ point_list[8] = {-0.893709900310366, -0.379295097932432};
+ point_list[9] = {-0.379295097932432, -0.893709900310366};
+ point_list[10] = {+0.273004998242798, -0.379295097932432};
+ point_list[11] = {-0.893709900310366, +0.273004998242798};
+
+ weight_list[0] = 0.233572551452758;
+ weight_list[1] = 0.233572551452758;
+ weight_list[2] = 0.233572551452758;
+ weight_list[3] = 0.101689812740414;
+ weight_list[4] = 0.101689812740414;
+ weight_list[5] = 0.101689812740414;
+ weight_list[6] = 0.165702151236748;
+ weight_list[7] = 0.165702151236748;
+ weight_list[8] = 0.165702151236748;
+ weight_list[9] = 0.165702151236748;
+ weight_list[10] = 0.165702151236748;
+ weight_list[11] = 0.165702151236748;
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 7: {
+ constexpr size_t nb_points = 13;
+ SmallArray<TinyVector<2>> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ point_list[0] = {-0.333333333333333, -0.333333333333333};
+ point_list[1] = {-0.479308067841920, -0.479308067841920};
+ point_list[2] = {-0.479308067841920, -0.041383864316160};
+ point_list[3] = {-0.041383864316160, -0.479308067841920};
+ point_list[4] = {-0.869739794195568, -0.869739794195568};
+ point_list[5] = {-0.869739794195568, +0.739479588391136};
+ point_list[6] = {+0.739479588391136, -0.869739794195568};
+ point_list[7] = {-0.374269007990252, +0.276888377139620};
+ point_list[8] = {+0.276888377139620, -0.902619369149368};
+ point_list[9] = {-0.902619369149368, -0.374269007990252};
+ point_list[10] = {-0.374269007990252, -0.902619369149368};
+ point_list[11] = {+0.276888377139620, -0.374269007990252};
+ point_list[12] = {-0.902619369149368, +0.276888377139620};
+
+ weight_list[0] = -0.299140088935364;
+ weight_list[1] = +0.351230514866416;
+ weight_list[2] = +0.351230514866416;
+ weight_list[3] = +0.351230514866416;
+ weight_list[4] = +0.106694471217676;
+ weight_list[5] = +0.106694471217676;
+ weight_list[6] = +0.106694471217676;
+ weight_list[7] = +0.154227521780514;
+ weight_list[8] = +0.154227521780514;
+ weight_list[9] = +0.154227521780514;
+ weight_list[10] = +0.154227521780514;
+ weight_list[11] = +0.154227521780514;
+ weight_list[12] = +0.154227521780514;
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ default: {
+ throw NormalError("Gauss-Legendre quadrature formulae handle orders up to 7 on triangles");
+ }
+ }
}
template <>
diff --git a/tests/test_SimplicialGaussLegendreQuadrature.cpp b/tests/test_SimplicialGaussLegendreQuadrature.cpp
index 80fd1c7e0..9517a9f1f 100644
--- a/tests/test_SimplicialGaussLegendreQuadrature.cpp
+++ b/tests/test_SimplicialGaussLegendreQuadrature.cpp
@@ -3,37 +3,38 @@
#include <catch2/matchers/catch_matchers_all.hpp>
#include <analysis/SimplicialGaussLegendreQuadrature.hpp>
+#include <geometry/AffineTransformation.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();
+ SECTION("1D")
+ {
+ 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();
+ auto point_list1 = quadrature_formula1.pointList();
+ auto weight_list1 = quadrature_formula1.weightList();
- bool is_same = true;
+ 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 < 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;
+ for (size_t i = 0; i < weight_list0.size(); ++i) {
+ if (weight_list0[i] != weight_list1[i]) {
+ return false;
+ }
}
- }
- return is_same;
- };
+ return is_same;
+ };
- SECTION("1D")
- {
auto is_symmetric_formula = [](auto quadrature_formula) {
auto point_list = quadrature_formula.pointList();
auto weight_list = quadrature_formula.weightList();
@@ -521,4 +522,207 @@ TEST_CASE("SimplicialGaussLegendreQuadrature", "[analysis]")
"error: Gauss-Legendre quadrature formulae handle orders 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_triangle = [](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_triangle](auto f, auto quadrature_formula, const double exact_value) {
+ using R2 = TinyVector<2>;
+ const double int_T_hat = integrate(f, quadrature_formula);
+ const double int_refined //
+ = integrate_on_triangle(f, quadrature_formula, {R2{-1, -1}, R2{0, -1}, R2{-1, 0}}) +
+ integrate_on_triangle(f, quadrature_formula, {R2{0, -1}, R2{0, 0}, R2{-1, 0}}) +
+ integrate_on_triangle(f, quadrature_formula, {R2{0, -1}, R2{1, -1}, R2{0, 0}}) +
+ integrate_on_triangle(f, quadrature_formula, {R2{-1, 0}, R2{0, 0}, R2{-1, 1}});
+
+ return -std::log((int_refined - exact_value) / (int_T_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("order 1")
+ {
+ SimplicialGaussLegendreQuadrature<2> l1(1);
+
+ REQUIRE(l1.numberOfPoints() == 1);
+
+ REQUIRE(integrate(p0, l1) == Catch::Approx(4));
+ REQUIRE(integrate(p1, l1) == Catch::Approx(-4. / 3));
+ REQUIRE(integrate(p2, l1) != Catch::Approx(-19. / 3));
+
+ REQUIRE(get_order(p2, l1, -19. / 3) == Catch::Approx(2));
+ }
+
+ SECTION("order 2")
+ {
+ SimplicialGaussLegendreQuadrature<2> l2(2);
+
+ REQUIRE(l2.numberOfPoints() == 3);
+
+ REQUIRE(integrate(p0, l2) == Catch::Approx(4));
+ REQUIRE(integrate(p1, l2) == Catch::Approx(-4. / 3));
+ REQUIRE(integrate(p2, l2) == Catch::Approx(-19. / 3));
+ REQUIRE(integrate(p3, l2) != Catch::Approx(146. / 5));
+
+ // In a weird way, the formula achieves 4th order on this test
+ REQUIRE(get_order(p3, l2, 146. / 5) == Catch::Approx(4));
+ }
+
+ SECTION("order 3")
+ {
+ SimplicialGaussLegendreQuadrature<2> l3(3);
+
+ REQUIRE(l3.numberOfPoints() == 4);
+
+ REQUIRE(integrate(p0, l3) == Catch::Approx(4));
+ REQUIRE(integrate(p1, l3) == Catch::Approx(-4. / 3));
+ REQUIRE(integrate(p2, l3) == Catch::Approx(-19. / 3));
+ REQUIRE(integrate(p3, l3) == Catch::Approx(146. / 5));
+ REQUIRE(integrate(p4, l3) != Catch::Approx(-25. / 6));
+
+ REQUIRE(get_order(p4, l3, -25. / 6) == Catch::Approx(4));
+ }
+
+ SECTION("order 4")
+ {
+ SimplicialGaussLegendreQuadrature<2> l4(4);
+
+ REQUIRE(l4.numberOfPoints() == 6);
+
+ REQUIRE(integrate(p0, l4) == Catch::Approx(4));
+ REQUIRE(integrate(p1, l4) == Catch::Approx(-4. / 3));
+ REQUIRE(integrate(p2, l4) == Catch::Approx(-19. / 3));
+ REQUIRE(integrate(p3, l4) == Catch::Approx(146. / 5));
+ REQUIRE(integrate(p4, l4) == Catch::Approx(-25. / 6));
+ REQUIRE(integrate(p5, l4) != Catch::Approx(-17. / 10));
+
+ // In a weird way, the formula achieves 6th order on this test
+ REQUIRE(get_order(p5, l4, -17. / 10) == Catch::Approx(6));
+ }
+
+ SECTION("order 5")
+ {
+ SimplicialGaussLegendreQuadrature<2> l5(5);
+
+ REQUIRE(l5.numberOfPoints() == 7);
+
+ REQUIRE(integrate(p0, l5) == Catch::Approx(4));
+ REQUIRE(integrate(p1, l5) == Catch::Approx(-4. / 3));
+ REQUIRE(integrate(p2, l5) == Catch::Approx(-19. / 3));
+ REQUIRE(integrate(p3, l5) == Catch::Approx(146. / 5));
+ REQUIRE(integrate(p4, l5) == Catch::Approx(-25. / 6));
+ REQUIRE(integrate(p5, l5) == Catch::Approx(-17. / 10));
+ REQUIRE(integrate(p6, l5) != Catch::Approx(-197. / 175));
+
+ REQUIRE(get_order(p6, l5, -197. / 175) == Catch::Approx(6));
+ }
+
+ SECTION("order 6")
+ {
+ SimplicialGaussLegendreQuadrature<2> l6(6);
+
+ REQUIRE(l6.numberOfPoints() == 12);
+
+ REQUIRE(integrate(p0, l6) == Catch::Approx(4));
+ REQUIRE(integrate(p1, l6) == Catch::Approx(-4. / 3));
+ REQUIRE(integrate(p2, l6) == Catch::Approx(-19. / 3));
+ REQUIRE(integrate(p3, l6) == Catch::Approx(146. / 5));
+ REQUIRE(integrate(p4, l6) == Catch::Approx(-25. / 6));
+ REQUIRE(integrate(p5, l6) == Catch::Approx(-17. / 10));
+ REQUIRE(integrate(p6, l6) == Catch::Approx(-197. / 175));
+ REQUIRE(integrate(p7, l6) != Catch::Approx(-4507. / 1575));
+
+ REQUIRE(get_order(p7, l6, -4507. / 1575) == Catch::Approx(8));
+ }
+
+ SECTION("order 7")
+ {
+ SimplicialGaussLegendreQuadrature<2> l7(7);
+
+ REQUIRE(l7.numberOfPoints() == 13);
+
+ REQUIRE(integrate(p0, l7) == Catch::Approx(4));
+ REQUIRE(integrate(p1, l7) == Catch::Approx(-4. / 3));
+ REQUIRE(integrate(p2, l7) == Catch::Approx(-19. / 3));
+ REQUIRE(integrate(p3, l7) == Catch::Approx(146. / 5));
+ REQUIRE(integrate(p4, l7) == Catch::Approx(-25. / 6));
+ REQUIRE(integrate(p5, l7) == Catch::Approx(-17. / 10));
+ REQUIRE(integrate(p6, l7) == Catch::Approx(-197. / 175));
+ REQUIRE(integrate(p7, l7) == Catch::Approx(-4507. / 1575));
+ REQUIRE(integrate(p8, l7) != Catch::Approx(-47734. / 315));
+
+ REQUIRE(get_order(p8, l7, -47734. / 315) == Catch::Approx(8));
+ }
+
+ SECTION("not implemented formulae")
+ {
+ REQUIRE_THROWS_WITH(SimplicialGaussLegendreQuadrature<2>(8),
+ "error: Gauss-Legendre quadrature formulae handle orders up to 7 on triangles");
+ }
+ }
}
--
GitLab