diff --git a/MathLib/Integration/GaussLegendreTet.cpp b/MathLib/Integration/GaussLegendreTet.cpp
index b2a7fdeaccd43928c840f52cccfbf5d0b3507b42..1faa686b15a83eed368b71e2ee6ebfb0563ece6c 100644
--- a/MathLib/Integration/GaussLegendreTet.cpp
+++ b/MathLib/Integration/GaussLegendreTet.cpp
@@ -11,10 +11,11 @@
namespace MathLib
{
-
-template <> const std::array<std::array<double, 3>, GaussLegendreTet<1>::NPoints>
-GaussLegendreTet<1>::X = {{{{1./3., 1./3., 1./3.}}}};
-template <> double const GaussLegendreTet<1>::W[1] = {1.0};
+template <>
+const std::array<std::array<double, 3>, GaussLegendreTet<1>::NPoints>
+ GaussLegendreTet<1>::X = {{{{1. / 4., 1. / 4., 1. / 4.}}}};
+template <>
+double const GaussLegendreTet<1>::W[1] = {1. / 6.};
const std::array<std::array<double, 3>, GaussLegendreTet<2>::NPoints>
GaussLegendreTet<2>::X = {{ {{1./4., 1./4., 1./4.}},
@@ -24,37 +25,41 @@ GaussLegendreTet<2>::X = {{ {{1./4., 1./4., 1./4.}},
{{1./6., 1./6., 1./2.}} }};
double const GaussLegendreTet<2>::W[5] = {-2./15., 0.075, 0.075, 0.075, 0.075};
+std::array<std::array<double, 3>, GaussLegendreTet<3>::NPoints> initGLTet3X()
+{
+ // Cf. Gellert, M., Harbord, R., 1991. Moderate degree cubature formulas for
+ // 3-D tetrahedral finite-element approximations. Communications in Applied
+ // Numerical Methods 7, 487–495. doi:10.1002/cnm.1630070609
+ const double a = 0.0673422422100983;
+ const double b = 0.3108859192633005;
+ const double c = 0.7217942490673264;
+ const double d = 0.0927352503108912;
+ const double e = 0.4544962958743506;
+ const double f = 0.045503704125649;
+
+ return {{{a, b, b},
+ {b, a, b},
+ {b, b, a},
+ {b, b, b},
+ {c, d, d},
+ {d, c, d},
+ {d, d, c},
+ {d, d, d},
+ {e, e, f},
+ {e, f, e},
+ {e, f, f},
+ {f, e, e},
+ {f, e, f},
+ {f, f, e}}};
+}
+
const std::array<std::array<double, 3>, GaussLegendreTet<3>::NPoints>
-GaussLegendreTet<3>::X = {{ {{1./2., 1./2., 1./2.}},
- {{0.09197107805272303, 0.09197107805272303, 0.09197107805272303}},
- {{0.72408676584183096, 0.09197107805272303, 0.09197107805272303}},
- {{0.09197107805272303, 0.72408676584183096, 0.09197107805272303}},
- {{0.09197107805272303, 0.09197107805272303, 0.72408676584183096}},
- {{0.44364916731037080, 0.05635083268962915, 0.05635083268962915}},
- {{0.05635083268962915, 0.44364916731037080, 0.05635083268962915}},
- {{0.05635083268962915, 0.05635083268962915, 0.44364916731037080}},
- {{0.05635083268962915, 0.44364916731037080, 0.44364916731037080}},
- {{0.44364916731037080, 0.05635083268962915, 0.44364916731037080}},
- {{0.44364916731037080, 0.44364916731037080, 0.05635083268962915}},
- {{0.31979362782962989, 0.31979362782962989, 0.31979362782962989}},
- {{0.04061911651111023, 0.31979362782962989, 0.31979362782962989}},
- {{0.31979362782962989, 0.04061911651111023, 0.31979362782962989}},
- {{0.31979362782962989, 0.31979362782962989, 0.04061911651111023}} }};
-double const GaussLegendreTet<3>::W[15] = {0.019753086419753086,
- 0.011989513963169772,
- 0.011989513963169772,
- 0.011989513963169772,
- 0.011989513963169772,
- 0.008818342151675485,
- 0.008818342151675485,
- 0.008818342151675485,
- 0.008818342151675485,
- 0.008818342151675485,
- 0.008818342151675485,
- 0.011511367871045397,
- 0.011511367871045397,
- 0.011511367871045397,
- 0.011511367871045397
- };
+ GaussLegendreTet<3>::X = initGLTet3X();
+
+static const double p = 0.1126879257180162 / 6.;
+static const double q = 0.0734930431163619 / 6.;
+static const double r = 0.0425460207770812 / 6.;
+double const GaussLegendreTet<3>::W[GaussLegendreTet<3>::NPoints] = {
+ p, p, p, p, q, q, q, q, r, r, r, r, r, r};
}
diff --git a/MathLib/Integration/GaussLegendreTet.h b/MathLib/Integration/GaussLegendreTet.h
index 4a6ac7e6c039bc04e4496ace6c46b0e2b142aac2..110635849f5bfdc34aedc4d5e56c83c374310272 100644
--- a/MathLib/Integration/GaussLegendreTet.h
+++ b/MathLib/Integration/GaussLegendreTet.h
@@ -38,7 +38,7 @@ struct GaussLegendreTet<2> {
template <>
struct GaussLegendreTet<3> {
static MATHLIB_EXPORT const unsigned Order = 3;
- static MATHLIB_EXPORT const unsigned NPoints = 15;
+ static MATHLIB_EXPORT const unsigned NPoints = 14;
static MATHLIB_EXPORT const std::array<std::array<double, 3>, NPoints> X;
static MATHLIB_EXPORT const double W[NPoints];
};
diff --git a/Tests/Data b/Tests/Data
index 611c3253b6254fe1ad1e53cb72b8b124ff11c16e..bb88a0209350e88e857f8e95b8845ab82ddd78c0 160000
--- a/Tests/Data
+++ b/Tests/Data
@@ -1 +1 @@
-Subproject commit 611c3253b6254fe1ad1e53cb72b8b124ff11c16e
+Subproject commit bb88a0209350e88e857f8e95b8845ab82ddd78c0
diff --git a/Tests/MathLib/TestGaussLegendreIntegration.cpp b/Tests/MathLib/TestGaussLegendreIntegration.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..abab841932bb84609bacb7937b91312856ea0506
--- /dev/null
+++ b/Tests/MathLib/TestGaussLegendreIntegration.cpp
@@ -0,0 +1,613 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/LICENSE.txt
+ */
+
+#include <gtest/gtest.h>
+
+#include <cmath>
+#include <limits>
+
+#include "BaseLib/BuildInfo.h"
+
+#include "MathLib/Integration/GaussLegendreTet.h"
+
+#include "MeshLib/Elements/Element.h"
+#include "MeshLib/IO/VtkIO/VtuInterface.h"
+#include "MeshLib/MeshSubsets.h"
+
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+
+#include "ProcessLib/Utils/CreateLocalAssemblers.h"
+#include "ProcessLib/Utils/InitShapeMatrices.h"
+
+#include "Tests/VectorUtils.h"
+
+namespace GaussLegendreTest
+{
+template <typename Shp>
+static std::array<double, 3> interpolateNodeCoordinates(
+ MeshLib::Element const& e, Shp const& N)
+{
+ std::array<double, 3> res;
+
+ auto const* const nodes = e.getNodes();
+ auto node_coords = N;
+
+ for (std::size_t d = 0; d < 3; ++d)
+ {
+ for (unsigned ip = 0; ip < N.size(); ++ip)
+ {
+ node_coords[ip] = (*nodes[ip])[d];
+ }
+
+ res[d] = N.dot(node_coords);
+ }
+
+ return res;
+}
+
+class LocalAssemblerDataInterface
+{
+public:
+ using Function = std::function<double(std::array<double, 3> const&)>;
+
+ virtual double integrate(Function const& f,
+ unsigned const integration_order) const = 0;
+
+ virtual ~LocalAssemblerDataInterface() = default;
+};
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+class LocalAssemblerData : public LocalAssemblerDataInterface
+{
+ 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*/)
+ : _e(e)
+ {
+ if (is_axially_symmetric)
+ OGS_FATAL("Only testing Cartesian meshes!");
+ }
+
+ double integrate(const Function& f,
+ unsigned const integration_order) const override
+ {
+ double integral = 0;
+
+ auto const sms =
+ ProcessLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
+ IntegrationMethod, GlobalDim>(
+ _e, false /*is_axially_symmetric*/,
+ IntegrationMethod{integration_order});
+ IntegrationMethod integration_method{integration_order};
+
+ 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 coords = interpolateNodeCoordinates(_e, N);
+ auto const function_value = f(coords);
+ integral += function_value * sms[ip].detJ *
+ integration_method.getWeightedPoint(ip).getWeight();
+ }
+
+ return integral;
+ }
+
+private:
+ MeshLib::Element const& _e;
+};
+
+class IntegrationTestProcess
+{
+public:
+ using LocalAssembler = LocalAssemblerDataInterface;
+
+ IntegrationTestProcess(MeshLib::Mesh const& mesh,
+ unsigned const integration_order)
+ : _integration_order(integration_order),
+ _mesh_subset_all_nodes(mesh, &mesh.getNodes())
+ {
+ std::vector<MeshLib::MeshSubsets> all_mesh_subsets;
+ all_mesh_subsets.emplace_back(&_mesh_subset_all_nodes);
+
+ _dof_table = std::make_unique<NumLib::LocalToGlobalIndexMap>(
+ std::move(all_mesh_subsets), NumLib::ComponentOrder::BY_COMPONENT);
+
+ // createAssemblers(mesh);
+ ProcessLib::createLocalAssemblers<LocalAssemblerData>(
+ mesh.getDimension(), mesh.getElements(), *_dof_table, 1,
+ _local_assemblers, mesh.isAxiallySymmetric(), _integration_order);
+ }
+
+ double integrate(LocalAssembler::Function const& f) const
+ {
+ double integral = 0;
+ for (auto const& la : _local_assemblers)
+ {
+ integral += la->integrate(f, _integration_order);
+ }
+
+ return integral;
+ }
+
+private:
+ unsigned const _integration_order;
+
+ MeshLib::MeshSubset _mesh_subset_all_nodes;
+ std::unique_ptr<NumLib::LocalToGlobalIndexMap> _dof_table;
+
+ std::vector<std::unique_ptr<LocalAssembler>> _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:
+ 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;
+
+ LocalAssemblerDataInterface::Function getClosure() const
+ {
+ return [this](std::array<double, 3> const& coords) {
+ return this->operator()(coords);
+ };
+ }
+
+ virtual ~FBase() = default;
+
+ std::vector<double> const coeffs;
+};
+
+struct FConst final : FBase
+{
+ FConst() : FBase(1) {}
+
+ double operator()(std::array<double, 3> const&) const override
+ {
+ return coeffs[0];
+ }
+
+ double getAnalyticalIntegralOverUnitCube() const override
+ {
+ return coeffs[0];
+ }
+};
+
+struct FLin final : FBase
+{
+ 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
+{
+ FQuad() : FBase(10) {}
+
+ 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 +
+ coeffs[4] * x * y + coeffs[5] * y * z + coeffs[6] * z * x +
+ coeffs[7] * x * x + coeffs[8] * y * y + coeffs[9] * z * z;
+ }
+
+ double getAnalyticalIntegralOverUnitCube() const override
+ {
+ double const a = -.5;
+ double const b = .5;
+
+ double const a3 = a * a * a;
+ double const b3 = b * b * b;
+
+ return coeffs[0] + (coeffs[7] + coeffs[8] + coeffs[9]) * (b3 - a3) / 3.;
+ }
+};
+
+struct F3DSeparablePolynomial final : FBase
+{
+ 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
+ {
+ double res = 1.0;
+ for (unsigned d : {0, 1, 2})
+ {
+ 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;
+ }
+
+ return res;
+ }
+
+ // [ F(x, y, z) ]_a^b = [ G(x) ]_a^b * [ H(y) ]_a^b * [ I(z) ]_a^b
+ double getAnalyticalIntegralOverUnitCube() const override
+ {
+ double const a = -.5;
+ double const b = .5;
+
+ double res = 1.0;
+ for (unsigned d : {0, 1, 2})
+ {
+ 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;
+ }
+
+ return res;
+ }
+
+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;
+
+ 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
+
+/* *****************************************************************************
+ *
+ * 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.
+ *
+ **************************************************************************** */
+
+static double const eps = 10 * std::numeric_limits<double>::epsilon();
+
+/* 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)
+#else
+#define OGS_DONT_TEST_THIS_IF_PETSC(group, test_case) \
+ TEST(group, DISABLED_##test_case)
+#endif
+
+OGS_DONT_TEST_THIS_IF_PETSC(MathLib, IntegrationGaussLegendreTet)
+{
+ 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);
+ 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: %u.", 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(BaseLib::BuildInfo::data_path +
+ "/MathLib/unit_cube_hex.vtu"));
+
+ for (unsigned integration_order : {1, 2, 3})
+ {
+ DBUG("\n==== integration order: %u.\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: %u.", 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(BaseLib::BuildInfo::data_path +
+ "/MathLib/unit_cube_tet.vtu"));
+
+ for (unsigned integration_order : {1, 2, 3})
+ {
+ DBUG("\n==== integration order: %u.\n", integration_order);
+ GaussLegendreTest::IntegrationTestProcess 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);
+ GaussLegendreTest::F3DSeparablePolynomial 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,
+ IntegrationGaussLegendreHexSeparablePolynomial)
+{
+ 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);
+ 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: %u.", polynomial_order);
+ GaussLegendreTest::F3DSeparablePolynomial f(polynomial_order);
+
+ auto const integral_hex = pcs_hex.integrate(f.getClosure());
+ EXPECT_NEAR(f.getAnalyticalIntegralOverUnitCube(), integral_hex,
+ eps);
+ }
+ }
+}
+
+OGS_DONT_TEST_THIS_IF_PETSC(MathLib,
+ IntegrationGaussLegendreTetNonSeparablePolynomial)
+{
+ 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);
+ GaussLegendreTest::IntegrationTestProcess 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);
+ 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,
+ IntegrationGaussLegendreHexNonSeparablePolynomial)
+{
+ 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);
+ 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: %u.", 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
diff --git a/Tests/NumLib/TestExtrapolation.cpp b/Tests/NumLib/TestExtrapolation.cpp
index 075556d43ce5aec93cdad88f370e124ccfd2264e..10a82220df115efb9d7bf886d8357dca789626be 100644
--- a/Tests/NumLib/TestExtrapolation.cpp
+++ b/Tests/NumLib/TestExtrapolation.cpp
@@ -11,8 +11,8 @@
#include "MathLib/LinAlg/LinAlg.h"
-#include "MeshLib/MeshGenerators/MeshGenerator.h"
#include "MeshLib/IO/writeMeshToFile.h"
+#include "MeshLib/MeshGenerators/MeshGenerator.h"
#include "NumLib/DOF/DOFTableUtil.h"
#include "NumLib/DOF/MatrixProviderUser.h"
@@ -24,13 +24,13 @@
#include "NumLib/Function/Interpolation.h"
#include "NumLib/NumericsConfig.h"
-#include "ProcessLib/Utils/LocalDataInitializer.h"
#include "ProcessLib/Utils/CreateLocalAssemblers.h"
#include "ProcessLib/Utils/InitShapeMatrices.h"
+#include "ProcessLib/Utils/LocalDataInitializer.h"
#include "Tests/VectorUtils.h"
-namespace
+namespace ExtrapolationTest
{
template <typename ShapeMatrices>
void interpolateNodalValuesToIntegrationPoints(
@@ -39,15 +39,13 @@ void interpolateNodalValuesToIntegrationPoints(
shape_matrices,
std::vector<double>& interpolated_values)
{
- for (unsigned ip=0; ip<shape_matrices.size(); ++ip)
+ for (unsigned ip = 0; ip < shape_matrices.size(); ++ip)
{
NumLib::shapeFunctionInterpolate(
local_nodal_values, shape_matrices[ip].N, interpolated_values[ip]);
}
}
-} // anonymous namespace
-
class LocalAssemblerDataInterface : public NumLib::ExtrapolatableElement
{
public:
@@ -96,8 +94,8 @@ public:
{
}
- Eigen::Map<const Eigen::RowVectorXd>
- getShapeMatrix(const unsigned integration_point) const override
+ Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
+ const unsigned integration_point) const override
{
auto const& N = _shape_matrices[integration_point].N;
@@ -129,7 +127,7 @@ public:
void interpolateNodalValuesToIntegrationPoints(
std::vector<double> const& local_nodal_values) override
{
- ::interpolateNodalValuesToIntegrationPoints(
+ ExtrapolationTest::interpolateNodalValuesToIntegrationPoints(
local_nodal_values, _shape_matrices, _int_pt_values);
}
@@ -139,7 +137,7 @@ private:
std::vector<double> _int_pt_values;
};
-class TestProcess
+class ExtrapolationTestProcess
{
public:
using LocalAssembler = LocalAssemblerDataInterface;
@@ -148,9 +146,10 @@ public:
using ExtrapolatorImplementation =
NumLib::LocalLinearLeastSquaresExtrapolator;
- TestProcess(MeshLib::Mesh const& mesh, unsigned const integration_order)
- : _integration_order(integration_order)
- , _mesh_subset_all_nodes(mesh, &mesh.getNodes())
+ ExtrapolationTestProcess(MeshLib::Mesh const& mesh,
+ unsigned const integration_order)
+ : _integration_order(integration_order),
+ _mesh_subset_all_nodes(mesh, &mesh.getNodes())
{
std::vector<MeshLib::MeshSubsets> all_mesh_subsets;
all_mesh_subsets.emplace_back(&_mesh_subset_all_nodes);
@@ -211,13 +210,14 @@ private:
std::unique_ptr<ExtrapolatorInterface> _extrapolator;
};
-void extrapolate(TestProcess const& pcs, IntegrationPointValuesMethod method,
+void extrapolate(ExtrapolationTestProcess const& pcs,
+ IntegrationPointValuesMethod method,
GlobalVector const& expected_extrapolated_global_nodal_values,
std::size_t const nnodes, std::size_t const nelements)
{
namespace LinAlg = MathLib::LinAlg;
- auto const tolerance_dx = 30.0 * std::numeric_limits<double>::epsilon();
+ auto const tolerance_dx = 30.0 * std::numeric_limits<double>::epsilon();
auto const tolerance_res = 15.0 * std::numeric_limits<double>::epsilon();
const double t = 0.0;
@@ -227,7 +227,7 @@ void extrapolate(TestProcess const& pcs, IntegrationPointValuesMethod method,
auto const& x_extra = *result.first;
auto const& residual = *result.second;
- ASSERT_EQ(nnodes, x_extra.size());
+ ASSERT_EQ(nnodes, x_extra.size());
ASSERT_EQ(nelements, residual.size());
auto const res_norm = LinAlg::normMax(residual);
@@ -235,14 +235,16 @@ void extrapolate(TestProcess const& pcs, IntegrationPointValuesMethod method,
EXPECT_GT(tolerance_res, res_norm);
auto delta_x = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
- expected_extrapolated_global_nodal_values);
- LinAlg::axpy(*delta_x, -1.0, x_extra); // delta_x = x_expected - x_extra
+ expected_extrapolated_global_nodal_values);
+ LinAlg::axpy(*delta_x, -1.0, x_extra); // delta_x = x_expected - x_extra
auto const dx_norm = LinAlg::normMax(*delta_x);
DBUG("maximum norm of delta x: %g", dx_norm);
EXPECT_GT(tolerance_dx, dx_norm);
}
+} // anonymous namespace
+
#ifndef USE_PETSC
TEST(NumLib, Extrapolation)
#else
@@ -264,18 +266,20 @@ TEST(NumLib, DISABLED_Extrapolation)
// generate mesh
std::unique_ptr<MeshLib::Mesh> mesh(
- MeshLib::MeshGenerator::generateRegularHexMesh(
- mesh_length, mesh_elements_in_each_direction));
+ MeshLib::MeshGenerator::generateRegularHexMesh(
+ mesh_length, mesh_elements_in_each_direction));
for (unsigned integration_order : {2, 3, 4})
{
namespace LinAlg = MathLib::LinAlg;
- auto const nnodes = mesh->getNumberOfNodes();
+ auto const nnodes = mesh->getNumberOfNodes();
auto const nelements = mesh->getNumberOfElements();
- DBUG("number of nodes: %lu, number of elements: %lu", nnodes, nelements);
+ DBUG("number of nodes: %lu, number of elements: %lu", nnodes,
+ nelements);
- TestProcess pcs(*mesh, integration_order);
+ ExtrapolationTest::ExtrapolationTestProcess pcs(*mesh,
+ integration_order);
// generate random nodal values
MathLib::MatrixSpecifications spec{nnodes, nnodes, nullptr, nullptr};
@@ -287,16 +291,20 @@ TEST(NumLib, DISABLED_Extrapolation)
// test extrapolation of a quantity that is stored in the local
// assembler
- extrapolate(pcs, &LocalAssemblerDataInterface::getStoredQuantity, *x,
- nnodes, nelements);
+ ExtrapolationTest::extrapolate(
+ pcs,
+ &ExtrapolationTest::LocalAssemblerDataInterface::getStoredQuantity,
+ *x, nnodes, nelements);
// expect 2*x as extraplation result for derived quantity
auto two_x = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(*x);
- LinAlg::axpy(*two_x, 1.0, *x); // two_x = x + x
+ LinAlg::axpy(*two_x, 1.0, *x); // two_x = x + x
// test extrapolation of a quantity that is derived from some
// integration point values
- extrapolate(pcs, &LocalAssemblerDataInterface::getDerivedQuantity,
- *two_x, nnodes, nelements);
+ ExtrapolationTest::extrapolate(
+ pcs,
+ &ExtrapolationTest::LocalAssemblerDataInterface::getDerivedQuantity,
+ *two_x, nnodes, nelements);
}
}
diff --git a/Tests/VectorUtils.h b/Tests/VectorUtils.h
index d8c9adc0527080dcb4f9b403e28ad85be350bbb1..673f2083565ea4de47e6dfb6494ffcf6fde40744 100644
--- a/Tests/VectorUtils.h
+++ b/Tests/VectorUtils.h
@@ -28,3 +28,15 @@ void fillVectorRandomly(Vector& x)
finalizeVectorAssembly(x);
#endif
}
+
+inline void fillVectorRandomly(std::vector<double>& x)
+{
+ std::random_device rd;
+ std::mt19937 random_number_generator(rd());
+ std::uniform_real_distribution<double> rnd;
+
+ for (auto& value : x)
+ {
+ value = rnd(random_number_generator);
+ }
+}