From a43b25333fad390e83c875c33096af4348d8b82d Mon Sep 17 00:00:00 2001
From: Christoph Lehmann <christoph.lehmann@ufz.de>
Date: Tue, 15 Aug 2017 13:59:03 +0200
Subject: [PATCH] [T] added a further class of polynomials
---
Tests/MathLib/TestIntegration2.cpp | 174 ++++++++++++++++++++++++++++-
1 file changed, 173 insertions(+), 1 deletion(-)
diff --git a/Tests/MathLib/TestIntegration2.cpp b/Tests/MathLib/TestIntegration2.cpp
index 1bca6fc49a4..709cf181d72 100644
--- a/Tests/MathLib/TestIntegration2.cpp
+++ b/Tests/MathLib/TestIntegration2.cpp
@@ -296,6 +296,108 @@ private:
unsigned const _degree;
};
+unsigned binomial_coefficient(unsigned n, unsigned k)
+{
+ EXPECT_GE(n, k);
+ unsigned res = 1;
+
+ for (unsigned i = n; i > k; --i)
+ {
+ res *= i;
+ }
+
+ for (unsigned i = n - k; i > 0; --i)
+ {
+ res /= i;
+ }
+
+ return res;
+}
+
+/* This function is a polynomial where for each monomial a_ijk x^i y^j z^k
+ * holds: i + j + k <= n, where n is the overall polynomial degree
+ */
+struct F3DNonseparablePolynomial final : FBase
+{
+ // The number of coefficients/monomials are obtained as follows: Compute the
+ // number of combinations with repititions when drawing
+ // polynomial_degree times from the set { x, y, z, 1 }
+ F3DNonseparablePolynomial(unsigned polynomial_degree)
+ : FBase(binomial_coefficient(4 + polynomial_degree - 1, 4 - 1)),
+ _degree(polynomial_degree)
+ {
+ }
+
+ double operator()(std::array<double, 3> const& coords) const override
+ {
+ auto const x = coords[0];
+ auto const y = coords[1];
+ auto const z = coords[2];
+
+ double res = 0.0;
+ unsigned index = 0;
+
+ for (unsigned x_deg = 0; x_deg <= _degree; ++x_deg)
+ {
+ for (unsigned y_deg = 0; x_deg + y_deg <= _degree; ++y_deg)
+ {
+ for (unsigned z_deg = 0; x_deg + y_deg + z_deg <= _degree;
+ ++z_deg)
+ {
+ EXPECT_GT(coeffs.size(), index);
+
+ res += coeffs[index] * std::pow(x, x_deg) *
+ std::pow(y, y_deg) * std::pow(z, z_deg);
+
+ ++index;
+ }
+ }
+ }
+
+ EXPECT_EQ(coeffs.size(), index);
+
+ return res;
+ }
+
+ double getAnalyticalIntegralOverUnitCube() const override
+ {
+ double const a = -.5;
+ double const b = .5;
+
+ double res = 0.0;
+ unsigned index = 0;
+
+ for (unsigned x_deg = 0; x_deg <= _degree; ++x_deg)
+ {
+ for (unsigned y_deg = 0; x_deg + y_deg <= _degree; ++y_deg)
+ {
+ for (unsigned z_deg = 0; x_deg + y_deg + z_deg <= _degree;
+ ++z_deg)
+ {
+ EXPECT_GT(coeffs.size(), index);
+
+ res += coeffs[index] *
+ (std::pow(b, x_deg + 1) - std::pow(a, x_deg + 1)) /
+ (x_deg + 1) *
+ (std::pow(b, y_deg + 1) - std::pow(a, y_deg + 1)) /
+ (y_deg + 1) *
+ (std::pow(b, z_deg + 1) - std::pow(a, z_deg + 1)) /
+ (z_deg + 1);
+
+ ++index;
+ }
+ }
+ }
+
+ EXPECT_EQ(coeffs.size(), index);
+
+ return res;
+ }
+
+private:
+ unsigned const _degree;
+};
+
std::unique_ptr<FBase> getF(unsigned polynomial_order)
{
std::vector<double> coeffs;
@@ -313,6 +415,18 @@ std::unique_ptr<FBase> getF(unsigned polynomial_order)
OGS_FATAL("unsupported polynomial order: %d.", polynomial_order);
}
+/* *****************************************************************************
+ *
+ * The idea behind the tests in this file is to integrate polynomials of
+ * different degree over the unit cube.
+ *
+ * Gauss-Legendre integration should be able to exactly integrate those up to a
+ * certian degree.
+ *
+ * The coefficients of the tested polynomials are chosen randomly.
+ *
+ **************************************************************************** */
+
TEST(MathLib, IntegrationGaussLegendreTet)
{
auto const eps = 2 * std::numeric_limits<double>::epsilon();
@@ -367,7 +481,9 @@ TEST(MathLib, IntegrationGaussLegendreHex)
}
}
-TEST(MathLib, IntegrationGaussLegendreTetSeparablePolynomial)
+// This test is disabled, because the polynomials involved are too complicated
+// to be exactly integrated over tetrahedra using Gauss-Legendre quadrature
+TEST(MathLib, DISABLED_IntegrationGaussLegendreTetSeparablePolynomial)
{
auto const eps = 2 * std::numeric_limits<double>::epsilon();
@@ -422,3 +538,59 @@ TEST(MathLib, IntegrationGaussLegendreHexSeparablePolynomial)
}
}
}
+
+TEST(MathLib, IntegrationGaussLegendreTetNonSeparablePolynomial)
+{
+ auto const eps = 2 * std::numeric_limits<double>::epsilon();
+
+ std::unique_ptr<MeshLib::Mesh> mesh_tet(
+ MeshLib::IO::VtuInterface::readVTUFile(BaseLib::BuildInfo::data_path +
+ "/MathLib/unit_cube_tet.vtu"));
+
+ for (unsigned integration_order : {1, 2, 3})
+ {
+ DBUG("\n==== integration order: %u.\n", integration_order);
+ TestProcess pcs_tet(*mesh_tet, integration_order);
+
+ for (unsigned polynomial_order = 0;
+ // Gauss-Legendre integration is exact up to this order!
+ polynomial_order < 2 * integration_order;
+ ++polynomial_order)
+ {
+ DBUG(" == polynomial order: %u.", polynomial_order);
+ F3DNonseparablePolynomial f(polynomial_order);
+
+ auto const integral_tet = pcs_tet.integrate(f.getClosure());
+ EXPECT_NEAR(f.getAnalyticalIntegralOverUnitCube(), integral_tet,
+ eps);
+ }
+ }
+}
+
+TEST(MathLib, IntegrationGaussLegendreHexNonSeparablePolynomial)
+{
+ auto const eps = 10 * std::numeric_limits<double>::epsilon();
+
+ std::unique_ptr<MeshLib::Mesh> mesh_hex(
+ MeshLib::IO::VtuInterface::readVTUFile(BaseLib::BuildInfo::data_path +
+ "/MathLib/unit_cube_hex.vtu"));
+
+ for (unsigned integration_order : {1, 2, 3, 4})
+ {
+ DBUG("\n==== integration order: %u.\n", integration_order);
+ TestProcess pcs_hex(*mesh_hex, integration_order);
+
+ for (unsigned polynomial_order = 0;
+ // Gauss-Legendre integration is exact up to this order!
+ polynomial_order < 2 * integration_order;
+ ++polynomial_order)
+ {
+ DBUG(" == polynomial order: %u.", polynomial_order);
+ F3DNonseparablePolynomial f(polynomial_order);
+
+ auto const integral_hex = pcs_hex.integrate(f.getClosure());
+ EXPECT_NEAR(f.getAnalyticalIntegralOverUnitCube(), integral_hex,
+ eps);
+ }
+ }
+}
--
GitLab