From 1b387c723f9d58186cbac38686f59ceab07f1a5d Mon Sep 17 00:00:00 2001
From: Christoph Lehmann <christoph.lehmann@ufz.de>
Date: Tue, 19 Apr 2022 10:39:22 +0200
Subject: [PATCH] [T] Extended quadrature tests to all schemes and orders
---
.../MathLib/TestGaussLegendreIntegration.cpp | 637 +++++++-----------
1 file changed, 241 insertions(+), 396 deletions(-)
diff --git a/Tests/MathLib/TestGaussLegendreIntegration.cpp b/Tests/MathLib/TestGaussLegendreIntegration.cpp
index 8b8898264d4..653ac0bd376 100644
--- a/Tests/MathLib/TestGaussLegendreIntegration.cpp
+++ b/Tests/MathLib/TestGaussLegendreIntegration.cpp
@@ -9,6 +9,7 @@
#include <gtest/gtest.h>
+#include <algorithm>
#include <cmath>
#include <limits>
@@ -17,10 +18,17 @@
#include "MeshLib/Elements/Element.h"
#include "MeshLib/IO/VtkIO/VtuInterface.h"
#include "MeshLib/MeshSubset.h"
+
+#ifdef USE_PETSC
+#include "MeshLib/NodePartitionedMesh.h"
+#endif
+
#include "NumLib/DOF/LocalToGlobalIndexMap.h"
#include "NumLib/Fem/InitShapeMatrices.h"
#include "NumLib/Fem/ShapeMatrixPolicy.h"
+#include "Polynomials.h"
#include "ProcessLib/Utils/CreateLocalAssemblers.h"
+#include "Tests/MeshLib/UnitCube.h"
#include "Tests/VectorUtils.h"
namespace GaussLegendreTest
@@ -47,7 +55,7 @@ static std::array<double, 3> interpolateNodeCoordinates(
return res;
}
-class LocalAssemblerDataInterface
+class LocalAssemblerInterface
{
public:
using Function = std::function<double(std::array<double, 3> const&)>;
@@ -55,20 +63,20 @@ public:
virtual double integrate(Function const& f,
unsigned const integration_order) const = 0;
- virtual ~LocalAssemblerDataInterface() = default;
+ virtual ~LocalAssemblerInterface() = default;
};
template <typename ShapeFunction, typename IntegrationMethod, int GlobalDim>
-class LocalAssemblerData : public LocalAssemblerDataInterface
+class LocalAssembler : public LocalAssemblerInterface
{
using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
public:
- LocalAssemblerData(MeshLib::Element const& e,
- std::size_t const /*local_matrix_size*/,
- bool is_axially_symmetric,
- unsigned const /*integration_order*/)
+ LocalAssembler(MeshLib::Element const& e,
+ std::size_t const /*local_matrix_size*/,
+ bool is_axially_symmetric,
+ unsigned const /*integration_order*/)
: _e(e)
{
if (is_axially_symmetric)
@@ -91,16 +99,12 @@ public:
for (unsigned ip = 0; ip < sms.size(); ++ip)
{
- // We need to assert that detJ is constant in each element in order
- // to be sure that in every element that we are integrating over the
- // effective polynomial degree is the one we expect.
- EXPECT_NEAR(sms[0].detJ, sms[ip].detJ,
- std::numeric_limits<double>::epsilon());
auto const& N = sms[ip].N;
+ auto const weight =
+ integration_method.getWeightedPoint(ip).getWeight();
auto const coords = interpolateNodeCoordinates(_e, N);
auto const function_value = f(coords);
- integral += function_value * sms[ip].detJ *
- integration_method.getWeightedPoint(ip).getWeight();
+ integral += function_value * sms[ip].detJ * weight;
}
return integral;
@@ -113,8 +117,6 @@ private:
class IntegrationTestProcess
{
public:
- using LocalAssembler = LocalAssemblerDataInterface;
-
IntegrationTestProcess(MeshLib::Mesh const& mesh,
unsigned const integration_order)
: _integration_order(integration_order),
@@ -124,15 +126,14 @@ public:
_mesh_subset_all_nodes};
_dof_table = std::make_unique<NumLib::LocalToGlobalIndexMap>(
- std::move(all_mesh_subsets), NumLib::ComponentOrder::BY_COMPONENT);
+ std::move(all_mesh_subsets), NumLib::ComponentOrder::BY_LOCATION);
- // createAssemblers(mesh);
- ProcessLib::createLocalAssemblers<LocalAssemblerData>(
+ ProcessLib::createLocalAssemblers<LocalAssembler>(
mesh.getDimension(), mesh.getElements(), *_dof_table,
_local_assemblers, mesh.isAxiallySymmetric(), _integration_order);
}
- double integrate(LocalAssembler::Function const& f) const
+ double integrate(LocalAssemblerInterface::Function const& f) const
{
double integral = 0;
for (auto const& la : _local_assemblers)
@@ -149,473 +150,317 @@ private:
MeshLib::MeshSubset _mesh_subset_all_nodes;
std::unique_ptr<NumLib::LocalToGlobalIndexMap> _dof_table;
- std::vector<std::unique_ptr<LocalAssembler>> _local_assemblers;
+ std::vector<std::unique_ptr<LocalAssemblerInterface>> _local_assemblers;
};
-struct FBase
-{
-private:
- static std::vector<double> initCoeffs(std::size_t const num_coeffs)
- {
- std::vector<double> cs(num_coeffs);
- fillVectorRandomly(cs);
- return cs;
- }
-
-public:
- explicit FBase(std::size_t const num_coeffs)
- : coeffs(initCoeffs(num_coeffs))
- {
- }
+} // namespace GaussLegendreTest
- virtual double operator()(
- std::array<double, 3> const& /*coords*/) const = 0;
- virtual double getAnalyticalIntegralOverUnitCube() const = 0;
+/* *****************************************************************************
+ *
+ * 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
+ * certain degree.
+ *
+ * The coefficients of the tested polynomials are chosen randomly.
+ *
+ **************************************************************************** */
- LocalAssemblerDataInterface::Function getClosure() const
+// Test fixture.
+class MathLibIntegrationGaussLegendreTestBase : public ::testing::Test
+{
+protected:
+ // Up to the polynomial order returned by this function the numerical
+ // integration should be exact in the base case test.
+ static unsigned maxExactPolynomialOrderBase(
+ unsigned const integration_order)
{
- return [this](std::array<double, 3> const& coords)
- { return this->operator()(coords); };
+ // Base case is defined only for polynomials up to order 2.
+ return std::min(2u, 2 * integration_order - 1);
}
- virtual ~FBase() = default;
-
- std::vector<double> const coeffs;
-};
-
-struct FConst final : FBase
-{
- FConst() : FBase(1) {}
-
- double operator()(std::array<double, 3> const& /*unused*/) const override
+ // Up to the polynomial order returned by this function the numerical
+ // integration should be exact in the separable polynomial test.
+ static unsigned maxExactPolynomialOrderSeparable(
+ unsigned const integration_order)
{
- return coeffs[0];
+ return 2 * integration_order - 1;
}
- double getAnalyticalIntegralOverUnitCube() const override
+ // Up to the polynomial order returned by this function the numerical
+ // integration should be exact in the nonseparable polynomial test.
+ static unsigned maxExactPolynomialOrderNonseparable(
+ unsigned const integration_order)
{
- return coeffs[0];
+ return 2 * integration_order - 1;
}
};
-struct FLin final : FBase
+template <typename MeshElementType>
+class MathLibIntegrationGaussLegendreTest
+ : public MathLibIntegrationGaussLegendreTestBase
{
- FLin() : FBase(4) {}
-
- 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];
- return coeffs[0] + coeffs[1] * x + coeffs[2] * y + coeffs[3] * z;
- }
-
- double getAnalyticalIntegralOverUnitCube() const override
- {
- return coeffs[0];
- }
};
-struct FQuad final : FBase
+// Specialization for triangles.
+template <>
+class MathLibIntegrationGaussLegendreTest<MeshLib::Tri>
+ : public MathLibIntegrationGaussLegendreTestBase
{
- FQuad() : FBase(10) {}
-
- double operator()(std::array<double, 3> const& coords) const override
+protected:
+ static unsigned maxExactPolynomialOrderSeparable(
+ unsigned const integration_order)
{
- auto const x = coords[0];
- auto const y = coords[1];
- auto const z = coords[2];
- return coeffs[0] + coeffs[1] * x + coeffs[2] * y + coeffs[3] * z +
- coeffs[4] * x * y + coeffs[5] * y * z + coeffs[6] * z * x +
- coeffs[7] * x * x + coeffs[8] * y * y + coeffs[9] * z * z;
+ switch (integration_order)
+ {
+ case 1:
+ return 0; // constants
+ case 2:
+ return 1; // most complicated monomial is x*y*z
+ case 3:
+ return 1; // most complicated monomial is x*y*z
+ case 4:
+ return 2; // most complicated monomial is x^2 * y^2 * z^2
+ }
+
+ OGS_FATAL("Unsupported integration order {}", integration_order);
}
- double getAnalyticalIntegralOverUnitCube() const override
+ static unsigned maxExactPolynomialOrderNonseparable(
+ unsigned const integration_order)
{
- double const a = -.5;
- double const b = .5;
-
- double const a3 = a * a * a;
- double const b3 = b * b * b;
+ switch (integration_order)
+ {
+ case 1:
+ return 1; // at most linear affine functions
+ case 2:
+ return 3; // most complicated monomial is x^3 or x*y*z
+ case 3:
+ return 3; // most complicated monomial is x^3 or x*y*z
+ case 4:
+ return 5; // most complicated monomial is x^5 or x^2 * y^2 * z
+ }
- return coeffs[0] + (coeffs[7] + coeffs[8] + coeffs[9]) * (b3 - a3) / 3.;
+ OGS_FATAL("Unsupported integration order {}", integration_order);
}
};
-struct F3DSeparablePolynomial final : FBase
+// Specialization for tetrahedra.
+template <>
+class MathLibIntegrationGaussLegendreTest<MeshLib::Tet>
+ : public MathLibIntegrationGaussLegendreTestBase
{
- explicit F3DSeparablePolynomial(unsigned polynomial_degree)
- : FBase(3 * polynomial_degree + 3), _degree(polynomial_degree)
- {
- }
-
- // f(x, y, z) = g(x) * h(y) * i(z)
- double operator()(std::array<double, 3> const& coords) const override
+protected:
+ static unsigned maxExactPolynomialOrderSeparable(
+ unsigned const integration_order)
{
- double res = 1.0;
- for (unsigned d : {0, 1, 2})
+ switch (integration_order)
{
- auto const x = coords[d];
-
- double poly = 0.0;
- for (unsigned n = 0; n <= _degree; ++n)
- {
- poly += coeffs[n + d * (_degree + 1)] * std::pow(x, n);
- }
-
- res *= poly;
+ case 1:
+ return 0; // constants
+ case 2:
+ return 1; // most complicated monomial is x*y*z
+ case 3:
+ return 1; // most complicated monomial is x*y*z
+ case 4:
+ return 1; // most complicated monomial is x*y*z
}
- return res;
+ OGS_FATAL("Unsupported integration order {}", integration_order);
}
- // [ F(x, y, z) ]_a^b = [ G(x) ]_a^b * [ H(y) ]_a^b * [ I(z) ]_a^b
- double getAnalyticalIntegralOverUnitCube() const override
+ static unsigned maxExactPolynomialOrderNonseparable(
+ unsigned const integration_order)
{
- double const a = -.5;
- double const b = .5;
-
- double res = 1.0;
- for (unsigned d : {0, 1, 2})
+ switch (integration_order)
{
- double poly = 0.0;
- for (unsigned n = 0; n <= _degree; ++n)
- {
- poly += coeffs[n + d * (_degree + 1)] *
- (std::pow(b, n + 1) - std::pow(a, n + 1)) / (n + 1);
- }
-
- res *= poly;
+ case 1:
+ return 1; // at most linear affine functions
+ case 2:
+ return 3; // most complicated monomial is x^3 or x*y*z
+ case 3:
+ return 5; // most complicated monomial is x^5 or x^2 * y^2 * z
+ case 4:
+ return 5; // most complicated monomial is x^5 or x^2 * y^2 * z
}
- return res;
+ OGS_FATAL("Unsupported integration order {}", integration_order);
}
-
-private:
- unsigned const _degree;
};
-unsigned binomial_coefficient(unsigned n, unsigned k)
+// Specialization for prisms.
+template <>
+class MathLibIntegrationGaussLegendreTest<MeshLib::Prism>
+ : public MathLibIntegrationGaussLegendreTestBase
{
- EXPECT_GE(n, k);
- unsigned res = 1;
-
- for (unsigned i = n; i > k; --i)
+protected:
+ static unsigned maxExactPolynomialOrderSeparable(
+ unsigned const integration_order)
{
- res *= i;
- }
+ switch (integration_order)
+ {
+ case 1:
+ return 0; // constants
+ case 2:
+ return 1; // most complicated monomial is x*y*z
+ case 3:
+ return 2; // most complicated monomial is x^2 * y^2 * z^2
+ case 4:
+ return 2; // most complicated monomial is x^2 * y^2 * z^2
+ }
- for (unsigned i = n - k; i > 0; --i)
- {
- res /= i;
+ OGS_FATAL("Unsupported integration order {}", integration_order);
}
- 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 }
- explicit F3DNonseparablePolynomial(unsigned polynomial_degree)
- : FBase(binomial_coefficient(4 + polynomial_degree - 1, 4 - 1)),
- _degree(polynomial_degree)
+ static unsigned maxExactPolynomialOrderNonseparable(
+ unsigned const integration_order)
{
+ switch (integration_order)
+ {
+ case 1:
+ return 1; // at most linear affine functions
+ case 2:
+ return 3; // most complicated monomial is x^3 or x*y*z
+ case 3:
+ return 5; // most complicated monomial is x^5 or x^2 * y^2 * z
+ case 4:
+ return 5; // most complicated monomial is x^5 or x^2 * y^2 * z
+ }
+
+ OGS_FATAL("Unsupported integration order {}", integration_order);
}
+};
- double operator()(std::array<double, 3> const& coords) const override
+// Specialization for prisms.
+template <>
+class MathLibIntegrationGaussLegendreTest<MeshLib::Pyramid>
+ : public MathLibIntegrationGaussLegendreTestBase
+{
+protected:
+ static unsigned maxExactPolynomialOrderSeparable(
+ unsigned const integration_order)
{
- 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)
+ // Note: In pyramids det J changes over the element. So effectively we
+ // are integrating a polynomial of higher degree than indicated here.
+ switch (integration_order)
{
- 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;
- }
- }
+ case 1:
+ return 1; // most complicated monomial is x*y*z
+ case 2:
+ return 1; // most complicated monomial is x*y*z
+ case 3:
+ return 3; // most complicated monomial is x^3 * y^3 * z^3
+ case 4:
+ return 3; // most complicated monomial is x^3 * y^3 * z^3
}
- EXPECT_EQ(coeffs.size(), index);
-
- return res;
+ OGS_FATAL("Unsupported integration order {}", integration_order);
}
- double getAnalyticalIntegralOverUnitCube() const override
+ static unsigned maxExactPolynomialOrderNonseparable(
+ unsigned const integration_order)
{
- 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)
+ // Note: In pyramids det J changes over the element. So effectively we
+ // are integrating a polynomial of higher degree than indicated here.
+ switch (integration_order)
{
- 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;
- }
- }
+ case 1:
+ return 1; // at most linear affine functions
+ case 2:
+ return 3; // most complicated monomial is x^3 or x*y*z
+ case 3:
+ return 3; // most complicated monomial is x^3 or x*y*z
+ case 4:
+ return 3; // most complicated monomial is x^3 or x*y*z
}
- EXPECT_EQ(coeffs.size(), index);
-
- return res;
+ OGS_FATAL("Unsupported integration order {}", integration_order);
}
-
-private:
- unsigned const _degree;
};
-std::unique_ptr<FBase> getF(unsigned polynomial_order)
-{
- std::vector<double> coeffs;
-
- switch (polynomial_order)
- {
- case 0:
- return std::make_unique<FConst>();
- case 1:
- return std::make_unique<FLin>();
- case 2:
- return std::make_unique<FQuad>();
- }
-
- OGS_FATAL("unsupported polynomial order: {:d}.", polynomial_order);
-}
-
-} // namespace GaussLegendreTest
+using MeshElementTypes =
+ ::testing::Types<MeshLib::Line, MeshLib::Quad, MeshLib::Hex, MeshLib::Tri,
+ MeshLib::Tet, MeshLib::Prism, MeshLib::Pyramid>;
-/* *****************************************************************************
- *
- * 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
- * certain degree.
- *
- * The coefficients of the tested polynomials are chosen randomly.
- *
- **************************************************************************** */
+TYPED_TEST_SUITE(MathLibIntegrationGaussLegendreTest, MeshElementTypes);
-static double const eps = 10 * std::numeric_limits<double>::epsilon();
+template <typename MeshElementType, typename MaxExactPolynomialOrderFct,
+ typename PolynomialFactory>
+static void mathLibIntegrationGaussLegendreTestImpl(
+ MaxExactPolynomialOrderFct const& max_exact_polynomial_order_fct,
+ PolynomialFactory const& polynomial_factory)
+{
+ auto const eps = 10 * std::numeric_limits<double>::epsilon();
+ auto const mesh_ptr = createUnitCube<MeshElementType>();
-/* The tests in this file fundamentally rely on a mesh being read from a vtu
- * file, and a DOF table being computed for that mesh. Since our PETSc build
- * only works with node partitioned meshes, it cannot digest the read meshes.
- */
#ifndef USE_PETSC
-#define OGS_DONT_TEST_THIS_IF_PETSC(group, test_case) TEST(group, test_case)
+ auto const& mesh = *mesh_ptr;
#else
-#define OGS_DONT_TEST_THIS_IF_PETSC(group, test_case) \
- TEST(group, DISABLED_##test_case)
+ MeshLib::NodePartitionedMesh const mesh{*mesh_ptr};
#endif
-OGS_DONT_TEST_THIS_IF_PETSC(MathLib, IntegrationGaussLegendreTet)
-{
- std::unique_ptr<MeshLib::Mesh> mesh_tet(
- MeshLib::IO::VtuInterface::readVTUFile(
- TestInfoLib::TestInfo::data_path + "/MathLib/unit_cube_tet.vtu"));
-
- for (unsigned integration_order : {1, 2, 3, 4})
- {
- DBUG("\n==== integration order: {:d}.\n", integration_order);
- GaussLegendreTest::IntegrationTestProcess pcs_tet(*mesh_tet,
- integration_order);
-
- for (unsigned polynomial_order : {0, 1, 2})
- {
- if (polynomial_order > 2 * integration_order - 1)
- {
- break;
- }
-
- DBUG(" == polynomial order: {:d}.", polynomial_order);
- auto f = GaussLegendreTest::getF(polynomial_order);
-
- auto const integral_tet = pcs_tet.integrate(f->getClosure());
- EXPECT_NEAR(f->getAnalyticalIntegralOverUnitCube(), integral_tet,
- eps);
- }
- }
-}
-
-OGS_DONT_TEST_THIS_IF_PETSC(MathLib, IntegrationGaussLegendreHex)
-{
- std::unique_ptr<MeshLib::Mesh> mesh_hex(
- MeshLib::IO::VtuInterface::readVTUFile(
- TestInfoLib::TestInfo::data_path + "/MathLib/unit_cube_hex.vtu"));
-
- for (unsigned integration_order : {1, 2, 3})
- {
- DBUG("\n==== integration order: {:d}.\n", integration_order);
- GaussLegendreTest::IntegrationTestProcess pcs_hex(*mesh_hex,
- integration_order);
-
- for (unsigned polynomial_order : {0, 1, 2})
- {
- if (polynomial_order > 2 * integration_order - 1)
- {
- break;
- }
-
- DBUG(" == polynomial order: {:d}.", polynomial_order);
- auto f = GaussLegendreTest::getF(polynomial_order);
-
- auto const integral_hex = pcs_hex.integrate(f->getClosure());
- EXPECT_NEAR(f->getAnalyticalIntegralOverUnitCube(), integral_hex,
- eps);
- }
- }
-}
-
-// This test is disabled, because the polynomials involved are too complicated
-// to be exactly integrated over tetrahedra using Gauss-Legendre quadrature
-OGS_DONT_TEST_THIS_IF_PETSC(
- MathLib, DISABLED_IntegrationGaussLegendreTetSeparablePolynomial)
-{
- std::unique_ptr<MeshLib::Mesh> mesh_tet(
- MeshLib::IO::VtuInterface::readVTUFile(
- TestInfoLib::TestInfo::data_path + "/MathLib/unit_cube_tet.vtu"));
-
for (unsigned integration_order : {1, 2, 3, 4})
{
- DBUG("\n==== integration order: {:d}.\n", integration_order);
- GaussLegendreTest::IntegrationTestProcess pcs_tet(*mesh_tet,
- integration_order);
+ GaussLegendreTest::IntegrationTestProcess pcs(mesh, integration_order);
for (unsigned polynomial_order = 0;
- // Gauss-Legendre integration is exact up to this order!
- polynomial_order < 2 * integration_order;
+ polynomial_order <=
+ max_exact_polynomial_order_fct(integration_order);
++polynomial_order)
{
- DBUG(" == polynomial order: {:d}.", polynomial_order);
- GaussLegendreTest::F3DSeparablePolynomial f(polynomial_order);
-
- auto const integral_tet = pcs_tet.integrate(f.getClosure());
- EXPECT_NEAR(f.getAnalyticalIntegralOverUnitCube(), integral_tet,
- eps);
+ auto const f = polynomial_factory(polynomial_order);
+
+ auto const analytical_integral =
+ f->getAnalyticalIntegralOverUnitCube();
+ auto const numerical_integral = pcs.integrate(*f);
+ EXPECT_NEAR(analytical_integral, numerical_integral, eps)
+ << "integration order: " << integration_order
+ << "\npolynomial order: " << polynomial_order //
+ << "\nf: " << *f;
}
}
}
-OGS_DONT_TEST_THIS_IF_PETSC(MathLib,
- IntegrationGaussLegendreHexSeparablePolynomial)
+TYPED_TEST(MathLibIntegrationGaussLegendreTest, Base)
{
- std::unique_ptr<MeshLib::Mesh> mesh_hex(
- MeshLib::IO::VtuInterface::readVTUFile(
- TestInfoLib::TestInfo::data_path + "/MathLib/unit_cube_hex.vtu"));
-
- for (unsigned integration_order : {1, 2, 3, 4})
+ using MeshElementType = TypeParam;
+ auto const polynomial_factory = [](unsigned const polynomial_order)
{
- DBUG("\n==== integration order: {:d}.\n", integration_order);
- GaussLegendreTest::IntegrationTestProcess 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: {:d}.", polynomial_order);
- GaussLegendreTest::F3DSeparablePolynomial f(polynomial_order);
+ auto constexpr dim = MeshElementType::dimension;
+ return TestPolynomials::getF<dim>(polynomial_order);
+ };
- auto const integral_hex = pcs_hex.integrate(f.getClosure());
- EXPECT_NEAR(f.getAnalyticalIntegralOverUnitCube(), integral_hex,
- eps);
- }
- }
+ mathLibIntegrationGaussLegendreTestImpl<MeshElementType>(
+ this->maxExactPolynomialOrderBase, polynomial_factory);
}
-OGS_DONT_TEST_THIS_IF_PETSC(MathLib,
- IntegrationGaussLegendreTetNonSeparablePolynomial)
+TYPED_TEST(MathLibIntegrationGaussLegendreTest, SeparablePolynomial)
{
- std::unique_ptr<MeshLib::Mesh> mesh_tet(
- MeshLib::IO::VtuInterface::readVTUFile(
- TestInfoLib::TestInfo::data_path + "/MathLib/unit_cube_tet.vtu"));
-
- // quadrature must be exact up to the given polynomial degrees
- // for integration order n = 1, 2, 3 it is 2*n-1 = 1, 3, 5
- // for integration order 4 it is 5, i.e. the same as for 3rd integration
- // order
- const std::vector<unsigned> maximum_polynomial_degrees{0, 1, 3, 5, 5};
+ using MeshElementType = TypeParam;
- for (unsigned integration_order : {1, 2, 3, 4})
+ auto const polynomial_factory = [](unsigned const polynomial_order)
{
- DBUG("\n==== integration order: {:d}.\n", integration_order);
- GaussLegendreTest::IntegrationTestProcess pcs_tet(*mesh_tet,
- integration_order);
+ auto constexpr dim = MeshElementType::dimension;
+ return std::make_unique<TestPolynomials::FSeparablePolynomial<dim>>(
+ polynomial_order);
+ };
- for (unsigned polynomial_order = 0;
- polynomial_order <= maximum_polynomial_degrees[integration_order];
- ++polynomial_order)
- {
- DBUG(" == polynomial order: {:d}.", polynomial_order);
- GaussLegendreTest::F3DNonseparablePolynomial f(polynomial_order);
-
- auto const integral_tet = pcs_tet.integrate(f.getClosure());
- EXPECT_NEAR(f.getAnalyticalIntegralOverUnitCube(), integral_tet,
- eps);
- }
- }
+ mathLibIntegrationGaussLegendreTestImpl<MeshElementType>(
+ this->maxExactPolynomialOrderSeparable, polynomial_factory);
}
-OGS_DONT_TEST_THIS_IF_PETSC(MathLib,
- IntegrationGaussLegendreHexNonSeparablePolynomial)
+TYPED_TEST(MathLibIntegrationGaussLegendreTest, NonseparablePolynomial)
{
- std::unique_ptr<MeshLib::Mesh> mesh_hex(
- MeshLib::IO::VtuInterface::readVTUFile(
- TestInfoLib::TestInfo::data_path + "/MathLib/unit_cube_hex.vtu"));
+ using MeshElementType = TypeParam;
- for (unsigned integration_order : {1, 2, 3, 4})
+ auto const polynomial_factory = [](unsigned const polynomial_order)
{
- DBUG("\n==== integration order: {:d}.\n", integration_order);
- GaussLegendreTest::IntegrationTestProcess pcs_hex(*mesh_hex,
- integration_order);
+ auto constexpr dim = MeshElementType::dimension;
+ return std::make_unique<TestPolynomials::FNonseparablePolynomial<dim>>(
+ polynomial_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: {:d}.", polynomial_order);
- GaussLegendreTest::F3DNonseparablePolynomial f(polynomial_order);
-
- auto const integral_hex = pcs_hex.integrate(f.getClosure());
- EXPECT_NEAR(f.getAnalyticalIntegralOverUnitCube(), integral_hex,
- eps);
- }
- }
+ mathLibIntegrationGaussLegendreTestImpl<MeshElementType>(
+ this->maxExactPolynomialOrderNonseparable, polynomial_factory);
}
-
-#undef OGS_DONT_TEST_THIS_IF_PETSC
--
GitLab