From 462b8877383eb6e10f0b2d7dbe5cca47b95b4770 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Mon, 4 Oct 2021 12:14:36 +0200
Subject: [PATCH] Fix quadrature tests and add and order evaluation
functionality
---
tests/test_FunctionEvaluation.cpp | 38 +++++++++++++++++++++++++++----
1 file changed, 33 insertions(+), 5 deletions(-)
diff --git a/tests/test_FunctionEvaluation.cpp b/tests/test_FunctionEvaluation.cpp
index 88c9e552a..7893ba07e 100644
--- a/tests/test_FunctionEvaluation.cpp
+++ b/tests/test_FunctionEvaluation.cpp
@@ -38,16 +38,29 @@ TEST_CASE("QuadratureFormula", "[analysis]")
auto point_list = quadrature_formula.pointList();
auto weight_list = quadrature_formula.weightList();
- auto x = [&a, &b](auto x_hat) { return 0.5 * (b - a) * (x_hat + 1); };
+ 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 0.5 * (b - a) * value;
+
+ 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 3; };
+ auto p0 = [](const TinyVector<1>&) { return 1; };
auto p1 = [](const TinyVector<1>& X) {
const double x = X[0];
return 2 * x + 1;
@@ -56,6 +69,14 @@ TEST_CASE("QuadratureFormula", "[analysis]")
const double x = X[0];
return 3 * x * x + 2 * x + 1;
};
+ auto p3 = [](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 4 * std::pow(x, 3) + 3 * x * x + 2 * x + 1;
+ };
+ auto p4 = [](const TinyVector<1>& X) {
+ const double x = X[0];
+ return 5 * std::pow(x, 4) + 4 * std::pow(x, 3) + 3 * x * x + 2 * x + 1;
+ };
SECTION("exact integration")
{
@@ -67,8 +88,11 @@ TEST_CASE("QuadratureFormula", "[analysis]")
REQUIRE(l0.numberOfPoints() == 1);
REQUIRE(l1.numberOfPoints() == 1);
- REQUIRE(integrate(p0, l1, 0.5, 1) == Catch::Approx(1.5));
+ 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, 0, 1, 3) == Catch::Approx(2));
}
SECTION("order 2 and 3")
@@ -79,9 +103,13 @@ TEST_CASE("QuadratureFormula", "[analysis]")
REQUIRE(l2.numberOfPoints() == 2);
REQUIRE(l3.numberOfPoints() == 2);
- REQUIRE(integrate(p0, l3, 0.5, 1) == Catch::Approx(1.5));
+ 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, 0, 1, 5) == Catch::Approx(4));
}
}
}
--
GitLab