Skip to content
Snippets Groups Projects
Commit 1b387c72 authored by Christoph Lehmann's avatar Christoph Lehmann
Browse files

[T] Extended quadrature tests to all schemes and orders

parent 1b6cbed7
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Tests/MathLib/TestGaussLegendreIntegration.cpp
+ 241
− 396
View file @ 1b387c72
...@@ -9,6 +9,7 @@ ...@@ -9,6 +9,7 @@
#include <gtest/gtest.h> #include <gtest/gtest.h>
#include <algorithm>
#include <cmath> #include <cmath>
#include <limits> #include <limits>
...@@ -17,10 +18,17 @@ ...@@ -17,10 +18,17 @@
#include "MeshLib/Elements/Element.h" #include "MeshLib/Elements/Element.h"
#include "MeshLib/IO/VtkIO/VtuInterface.h" #include "MeshLib/IO/VtkIO/VtuInterface.h"
#include "MeshLib/MeshSubset.h" #include "MeshLib/MeshSubset.h"
#ifdef USE_PETSC
#include "MeshLib/NodePartitionedMesh.h"
#endif
#include "NumLib/DOF/LocalToGlobalIndexMap.h" #include "NumLib/DOF/LocalToGlobalIndexMap.h"
#include "NumLib/Fem/InitShapeMatrices.h" #include "NumLib/Fem/InitShapeMatrices.h"
#include "NumLib/Fem/ShapeMatrixPolicy.h" #include "NumLib/Fem/ShapeMatrixPolicy.h"
#include "Polynomials.h"
#include "ProcessLib/Utils/CreateLocalAssemblers.h" #include "ProcessLib/Utils/CreateLocalAssemblers.h"
#include "Tests/MeshLib/UnitCube.h"
#include "Tests/VectorUtils.h" #include "Tests/VectorUtils.h"
namespace GaussLegendreTest namespace GaussLegendreTest
...@@ -47,7 +55,7 @@ static std::array<double, 3> interpolateNodeCoordinates( ...@@ -47,7 +55,7 @@ static std::array<double, 3> interpolateNodeCoordinates(
return res; return res;
} }
class LocalAssemblerDataInterface class LocalAssemblerInterface
{ {
public: public:
using Function = std::function<double(std::array<double, 3> const&)>; using Function = std::function<double(std::array<double, 3> const&)>;
...@@ -55,20 +63,20 @@ public: ...@@ -55,20 +63,20 @@ public:
virtual double integrate(Function const& f, virtual double integrate(Function const& f,
unsigned const integration_order) const = 0; unsigned const integration_order) const = 0;
virtual ~LocalAssemblerDataInterface() = default; virtual ~LocalAssemblerInterface() = default;
}; };
template <typename ShapeFunction, typename IntegrationMethod, int GlobalDim> template <typename ShapeFunction, typename IntegrationMethod, int GlobalDim>
class LocalAssemblerData : public LocalAssemblerDataInterface class LocalAssembler : public LocalAssemblerInterface
{ {
using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>; using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices; using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
public: public:
LocalAssemblerData(MeshLib::Element const& e, LocalAssembler(MeshLib::Element const& e,
std::size_t const /*local_matrix_size*/, std::size_t const /*local_matrix_size*/,
bool is_axially_symmetric, bool is_axially_symmetric,
unsigned const /*integration_order*/) unsigned const /*integration_order*/)
: _e(e) : _e(e)
{ {
if (is_axially_symmetric) if (is_axially_symmetric)
...@@ -91,16 +99,12 @@ public: ...@@ -91,16 +99,12 @@ public:
for (unsigned ip = 0; ip < sms.size(); ++ip) 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& N = sms[ip].N;
auto const weight =
integration_method.getWeightedPoint(ip).getWeight();
auto const coords = interpolateNodeCoordinates(_e, N); auto const coords = interpolateNodeCoordinates(_e, N);
auto const function_value = f(coords); auto const function_value = f(coords);
integral += function_value * sms[ip].detJ * integral += function_value * sms[ip].detJ * weight;
integration_method.getWeightedPoint(ip).getWeight();
} }
return integral; return integral;
...@@ -113,8 +117,6 @@ private: ...@@ -113,8 +117,6 @@ private:
class IntegrationTestProcess class IntegrationTestProcess
{ {
public: public:
using LocalAssembler = LocalAssemblerDataInterface;
IntegrationTestProcess(MeshLib::Mesh const& mesh, IntegrationTestProcess(MeshLib::Mesh const& mesh,
unsigned const integration_order) unsigned const integration_order)
: _integration_order(integration_order), : _integration_order(integration_order),
...@@ -124,15 +126,14 @@ public: ...@@ -124,15 +126,14 @@ public:
_mesh_subset_all_nodes}; _mesh_subset_all_nodes};
_dof_table = std::make_unique<NumLib::LocalToGlobalIndexMap>( _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<LocalAssembler>(
ProcessLib::createLocalAssemblers<LocalAssemblerData>(
mesh.getDimension(), mesh.getElements(), *_dof_table, mesh.getDimension(), mesh.getElements(), *_dof_table,
_local_assemblers, mesh.isAxiallySymmetric(), _integration_order); _local_assemblers, mesh.isAxiallySymmetric(), _integration_order);
} }
double integrate(LocalAssembler::Function const& f) const double integrate(LocalAssemblerInterface::Function const& f) const
{ {
double integral = 0; double integral = 0;
for (auto const& la : _local_assemblers) for (auto const& la : _local_assemblers)
...@@ -149,473 +150,317 @@ private: ...@@ -149,473 +150,317 @@ private:
MeshLib::MeshSubset _mesh_subset_all_nodes; MeshLib::MeshSubset _mesh_subset_all_nodes;
std::unique_ptr<NumLib::LocalToGlobalIndexMap> _dof_table; 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 } // namespace GaussLegendreTest
{
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))
{
}
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) // Base case is defined only for polynomials up to order 2.
{ return this->operator()(coords); }; return std::min(2u, 2 * integration_order - 1);
} }
virtual ~FBase() = default; // Up to the polynomial order returned by this function the numerical
// integration should be exact in the separable polynomial test.
std::vector<double> const coeffs; static unsigned maxExactPolynomialOrderSeparable(
}; unsigned const integration_order)
struct FConst final : FBase
{
FConst() : FBase(1) {}
double operator()(std::array<double, 3> const& /*unused*/) const override
{ {
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) {} protected:
static unsigned maxExactPolynomialOrderSeparable(
double operator()(std::array<double, 3> const& coords) const override unsigned const integration_order)
{ {
auto const x = coords[0]; switch (integration_order)
auto const y = coords[1]; {
auto const z = coords[2]; case 1:
return coeffs[0] + coeffs[1] * x + coeffs[2] * y + coeffs[3] * z + return 0; // constants
coeffs[4] * x * y + coeffs[5] * y * z + coeffs[6] * z * x + case 2:
coeffs[7] * x * x + coeffs[8] * y * y + coeffs[9] * z * z; 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; switch (integration_order)
double const b = .5; {
case 1:
double const a3 = a * a * a; return 1; // at most linear affine functions
double const b3 = b * b * b; 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) protected:
: FBase(3 * polynomial_degree + 3), _degree(polynomial_degree) static unsigned maxExactPolynomialOrderSeparable(
{ unsigned const integration_order)
}
// f(x, y, z) = g(x) * h(y) * i(z)
double operator()(std::array<double, 3> const& coords) const override
{ {
double res = 1.0; switch (integration_order)
for (unsigned d : {0, 1, 2})
{ {
auto const x = coords[d]; case 1:
return 0; // constants
double poly = 0.0; case 2:
for (unsigned n = 0; n <= _degree; ++n) return 1; // most complicated monomial is x*y*z
{ case 3:
poly += coeffs[n + d * (_degree + 1)] * std::pow(x, n); return 1; // most complicated monomial is x*y*z
} case 4:
return 1; // most complicated monomial is x*y*z
res *= poly;
} }
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 static unsigned maxExactPolynomialOrderNonseparable(
double getAnalyticalIntegralOverUnitCube() const override unsigned const integration_order)
{ {
double const a = -.5; switch (integration_order)
double const b = .5;
double res = 1.0;
for (unsigned d : {0, 1, 2})
{ {
double poly = 0.0; case 1:
for (unsigned n = 0; n <= _degree; ++n) return 1; // at most linear affine functions
{ case 2:
poly += coeffs[n + d * (_degree + 1)] * return 3; // most complicated monomial is x^3 or x*y*z
(std::pow(b, n + 1) - std::pow(a, n + 1)) / (n + 1); case 3:
} return 5; // most complicated monomial is x^5 or x^2 * y^2 * z
case 4:
res *= poly; 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); protected:
unsigned res = 1; static unsigned maxExactPolynomialOrderSeparable(
unsigned const integration_order)
for (unsigned i = n; i > k; --i)
{ {
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) OGS_FATAL("Unsupported integration order {}", integration_order);
{
res /= i;
} }
return res; static unsigned maxExactPolynomialOrderNonseparable(
} unsigned const integration_order)
/* 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)
{ {
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]; // Note: In pyramids det J changes over the element. So effectively we
auto const y = coords[1]; // are integrating a polynomial of higher degree than indicated here.
auto const z = coords[2]; switch (integration_order)
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) case 1:
{ return 1; // most complicated monomial is x*y*z
for (unsigned z_deg = 0; x_deg + y_deg + z_deg <= _degree; case 2:
++z_deg) return 1; // most complicated monomial is x*y*z
{ case 3:
EXPECT_GT(coeffs.size(), index); return 3; // most complicated monomial is x^3 * y^3 * z^3
case 4:
res += coeffs[index] * std::pow(x, x_deg) * return 3; // most complicated monomial is x^3 * y^3 * z^3
std::pow(y, y_deg) * std::pow(z, z_deg);
++index;
}
}
} }
EXPECT_EQ(coeffs.size(), index); OGS_FATAL("Unsupported integration order {}", integration_order);
return res;
} }
double getAnalyticalIntegralOverUnitCube() const override static unsigned maxExactPolynomialOrderNonseparable(
unsigned const integration_order)
{ {
double const a = -.5; // Note: In pyramids det J changes over the element. So effectively we
double const b = .5; // are integrating a polynomial of higher degree than indicated here.
switch (integration_order)
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) case 1:
{ return 1; // at most linear affine functions
for (unsigned z_deg = 0; x_deg + y_deg + z_deg <= _degree; case 2:
++z_deg) return 3; // most complicated monomial is x^3 or x*y*z
{ case 3:
EXPECT_GT(coeffs.size(), index); return 3; // most complicated monomial is x^3 or x*y*z
case 4:
res += coeffs[index] * return 3; // most complicated monomial is x^3 or x*y*z
(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); OGS_FATAL("Unsupported integration order {}", integration_order);
return res;
} }
private:
unsigned const _degree;
}; };
std::unique_ptr<FBase> getF(unsigned polynomial_order) using MeshElementTypes =
{ ::testing::Types<MeshLib::Line, MeshLib::Quad, MeshLib::Hex, MeshLib::Tri,
std::vector<double> coeffs; MeshLib::Tet, MeshLib::Prism, MeshLib::Pyramid>;
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
/* ***************************************************************************** TYPED_TEST_SUITE(MathLibIntegrationGaussLegendreTest, MeshElementTypes);
*
* 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.
*
**************************************************************************** */
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 #ifndef USE_PETSC
#define OGS_DONT_TEST_THIS_IF_PETSC(group, test_case) TEST(group, test_case) auto const& mesh = *mesh_ptr;
#else #else
#define OGS_DONT_TEST_THIS_IF_PETSC(group, test_case) \ MeshLib::NodePartitionedMesh const mesh{*mesh_ptr};
TEST(group, DISABLED_##test_case)
#endif #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}) for (unsigned integration_order : {1, 2, 3, 4})
{ {
DBUG("\n==== integration order: {:d}.\n", integration_order); GaussLegendreTest::IntegrationTestProcess pcs(mesh, integration_order);
GaussLegendreTest::IntegrationTestProcess pcs_tet(*mesh_tet,
integration_order);
for (unsigned polynomial_order = 0; for (unsigned polynomial_order = 0;
// Gauss-Legendre integration is exact up to this order! polynomial_order <=
polynomial_order < 2 * integration_order; max_exact_polynomial_order_fct(integration_order);
++polynomial_order) ++polynomial_order)
{ {
DBUG(" == polynomial order: {:d}.", polynomial_order); auto const f = polynomial_factory(polynomial_order);
GaussLegendreTest::F3DSeparablePolynomial f(polynomial_order);
auto const analytical_integral =
auto const integral_tet = pcs_tet.integrate(f.getClosure()); f->getAnalyticalIntegralOverUnitCube();
EXPECT_NEAR(f.getAnalyticalIntegralOverUnitCube(), integral_tet, auto const numerical_integral = pcs.integrate(*f);
eps); EXPECT_NEAR(analytical_integral, numerical_integral, eps)
<< "integration order: " << integration_order
<< "\npolynomial order: " << polynomial_order //
<< "\nf: " << *f;
} }
} }
} }
OGS_DONT_TEST_THIS_IF_PETSC(MathLib, TYPED_TEST(MathLibIntegrationGaussLegendreTest, Base)
IntegrationGaussLegendreHexSeparablePolynomial)
{ {
std::unique_ptr<MeshLib::Mesh> mesh_hex( using MeshElementType = TypeParam;
MeshLib::IO::VtuInterface::readVTUFile( auto const polynomial_factory = [](unsigned const polynomial_order)
TestInfoLib::TestInfo::data_path + "/MathLib/unit_cube_hex.vtu"));
for (unsigned integration_order : {1, 2, 3, 4})
{ {
DBUG("\n==== integration order: {:d}.\n", integration_order); auto constexpr dim = MeshElementType::dimension;
GaussLegendreTest::IntegrationTestProcess pcs_hex(*mesh_hex, return TestPolynomials::getF<dim>(polynomial_order);
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 const integral_hex = pcs_hex.integrate(f.getClosure()); mathLibIntegrationGaussLegendreTestImpl<MeshElementType>(
EXPECT_NEAR(f.getAnalyticalIntegralOverUnitCube(), integral_hex, this->maxExactPolynomialOrderBase, polynomial_factory);
eps);
}
}
} }
OGS_DONT_TEST_THIS_IF_PETSC(MathLib, TYPED_TEST(MathLibIntegrationGaussLegendreTest, SeparablePolynomial)
IntegrationGaussLegendreTetNonSeparablePolynomial)
{ {
std::unique_ptr<MeshLib::Mesh> mesh_tet( using MeshElementType = TypeParam;
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};
for (unsigned integration_order : {1, 2, 3, 4}) auto const polynomial_factory = [](unsigned const polynomial_order)
{ {
DBUG("\n==== integration order: {:d}.\n", integration_order); auto constexpr dim = MeshElementType::dimension;
GaussLegendreTest::IntegrationTestProcess pcs_tet(*mesh_tet, return std::make_unique<TestPolynomials::FSeparablePolynomial<dim>>(
integration_order); polynomial_order);
};
for (unsigned polynomial_order = 0; mathLibIntegrationGaussLegendreTestImpl<MeshElementType>(
polynomial_order <= maximum_polynomial_degrees[integration_order]; this->maxExactPolynomialOrderSeparable, polynomial_factory);
++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);
}
}
} }
OGS_DONT_TEST_THIS_IF_PETSC(MathLib, TYPED_TEST(MathLibIntegrationGaussLegendreTest, NonseparablePolynomial)
IntegrationGaussLegendreHexNonSeparablePolynomial)
{ {
std::unique_ptr<MeshLib::Mesh> mesh_hex( using MeshElementType = TypeParam;
MeshLib::IO::VtuInterface::readVTUFile(
TestInfoLib::TestInfo::data_path + "/MathLib/unit_cube_hex.vtu"));
for (unsigned integration_order : {1, 2, 3, 4}) auto const polynomial_factory = [](unsigned const polynomial_order)
{ {
DBUG("\n==== integration order: {:d}.\n", integration_order); auto constexpr dim = MeshElementType::dimension;
GaussLegendreTest::IntegrationTestProcess pcs_hex(*mesh_hex, return std::make_unique<TestPolynomials::FNonseparablePolynomial<dim>>(
integration_order); polynomial_order);
};
for (unsigned polynomial_order = 0; mathLibIntegrationGaussLegendreTestImpl<MeshElementType>(
// Gauss-Legendre integration is exact up to this order! this->maxExactPolynomialOrderNonseparable, polynomial_factory);
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);
}
}
} }
#undef OGS_DONT_TEST_THIS_IF_PETSC
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment