Select Git revision
test_CubeTransformation.cpp
Stéphane Del Pino authored
It avoids some spurious automatic conversion (especially from arithmetic types to TinyVector<1> or TinyMatrix<1>) and could reduce stress on the compiler.
test_CubeTransformation.cpp 9.60 KiB
#include <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp>
#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
#include <analysis/GaussLobattoQuadratureDescriptor.hpp>
#include <analysis/GaussQuadratureDescriptor.hpp>
#include <analysis/QuadratureManager.hpp>
#include <geometry/CubeTransformation.hpp>
// clazy:excludeall=non-pod-global-static
TEST_CASE("CubeTransformation", "[geometry]")
{
using R3 = TinyVector<3>;
SECTION("invertible transformation")
{
const R3 a_hat{-1, -1, -1};
const R3 b_hat{+1, -1, -1};
const R3 c_hat{+1, +1, -1};
const R3 d_hat{-1, +1, -1};
const R3 e_hat{-1, -1, +1};
const R3 f_hat{+1, -1, +1};
const R3 g_hat{+1, +1, +1};
const R3 h_hat{-1, +1, +1};
const R3 m_hat = zero;
const R3 a{0, 0, 0};
const R3 b{3, 1, 3};
const R3 c{2, 5, 2};
const R3 d{0, 3, 1};
const R3 e{1, 2, 5};
const R3 f{3, 1, 7};
const R3 g{2, 5, 5};
const R3 h{0, 3, 6};
const R3 m = 0.125 * (a + b + c + d + e + f + g + h);
const CubeTransformation t(a, b, c, d, e, f, g, h);
SECTION("values")
{
REQUIRE(l2Norm(t(a_hat) - a) == Catch::Approx(0));
REQUIRE(l2Norm(t(b_hat) - b) == Catch::Approx(0));
REQUIRE(l2Norm(t(c_hat) - c) == Catch::Approx(0));
REQUIRE(l2Norm(t(d_hat) - d) == Catch::Approx(0));
REQUIRE(l2Norm(t(e_hat) - e) == Catch::Approx(0));
REQUIRE(l2Norm(t(f_hat) - f) == Catch::Approx(0));
REQUIRE(l2Norm(t(g_hat) - g) == Catch::Approx(0));
REQUIRE(l2Norm(t(h_hat) - h) == Catch::Approx(0));
REQUIRE(l2Norm(t(m_hat) - m) == Catch::Approx(0));
}
SECTION("Jacobian determinant")
{
SECTION("at points")
{
auto detJ = [](const R3& X) {
const double x = X[0];
const double y = X[1];
const double z = X[2];
return -(3 * x * y * y * z + 9 * y * y * z - 6 * x * x * y * z + 6 * x * y * z - 96 * y * z + 14 * x * x * z -
49 * x * z + 119 * z - 21 * x * y * y + 21 * y * y + 90 * x * y - 10 * y + 40 * x * x - 109 * x -
491) /
128;
};
REQUIRE(t.jacobianDeterminant(a_hat) == Catch::Approx(detJ(a_hat)));
REQUIRE(t.jacobianDeterminant(b_hat) == Catch::Approx(detJ(b_hat)));
REQUIRE(t.jacobianDeterminant(c_hat) == Catch::Approx(detJ(c_hat)));
REQUIRE(t.jacobianDeterminant(d_hat) == Catch::Approx(detJ(d_hat)));
REQUIRE(t.jacobianDeterminant(e_hat) == Catch::Approx(detJ(e_hat)));
REQUIRE(t.jacobianDeterminant(f_hat) == Catch::Approx(detJ(f_hat)));
REQUIRE(t.jacobianDeterminant(g_hat) == Catch::Approx(detJ(g_hat)));
REQUIRE(t.jacobianDeterminant(h_hat) == Catch::Approx(detJ(h_hat)));
REQUIRE(t.jacobianDeterminant(m_hat) == Catch::Approx(detJ(m_hat)));
}
SECTION("Gauss order 3")
{
// The jacobian determinant is of maximal degree 2 in each variable
const QuadratureFormula<3>& gauss =
QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor(2));
double volume = 0;
for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
volume += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i));
}
// 353. / 12 is actually the volume of the hexahedron
REQUIRE(volume == Catch::Approx(353. / 12));
}
auto p = [](const R3& X) {
const double& x = X[0];
const double& y = X[1];
const double& z = X[2];
return 2 * x * x + 3 * y * z + z * z + 2 * x - 3 * y + 0.5 * z + 3;
};
SECTION("Gauss order 6")
{
// Jacbian determinant is a degree 2 polynomial, so the
// following formula is required to reach exactness
const QuadratureFormula<3>& gauss = QuadratureManager::instance().getCubeFormula(GaussQuadratureDescriptor(6));
double integral = 0;
for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(t(gauss.point(i)));
}
REQUIRE(integral == Catch::Approx(488879. / 360));
}
SECTION("Gauss-Legendre order 4")
{
// Jacbian determinant is a degree 2 polynomial, so the
// following formula is required to reach exactness
const QuadratureFormula<3>& gauss =
QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor(4));
double integral = 0;
for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(t(gauss.point(i)));
}
REQUIRE(integral == Catch::Approx(488879. / 360));
}
SECTION("Gauss-Lobatto order 4")
{
// Jacbian determinant is a degree 2 polynomial, so the
// following formula is required to reach exactness
const QuadratureFormula<3>& gauss =
QuadratureManager::instance().getCubeFormula(GaussLobattoQuadratureDescriptor(4));
double integral = 0;
for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(t(gauss.point(i)));
}
REQUIRE(integral == Catch::Approx(488879. / 360));
}
}
}
SECTION("degenerate to tetrahedron")
{
using R3 = TinyVector<3>;
const R3 a{1, 2, 1};
const R3 b{3, 1, 3};
const R3 c{2, 5, 2};
const R3 d{2, 3, 4};
const CubeTransformation t(a, b, c, c, d, d, d, d);
auto p = [](const R3& X) {
const double x = X[0];
const double y = X[1];
const double z = X[2];
return 2 * x * x + 3 * x * y + y * y + 3 * y - z * z + 2 * x * z + 2 * z;
};
SECTION("Polynomial Gauss integral")
{
QuadratureFormula<3> qf = QuadratureManager::instance().getCubeFormula(GaussQuadratureDescriptor(6));
double sum = 0;
for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
const R3 xi = qf.point(i);
sum += qf.weight(i) * t.jacobianDeterminant(xi) * p(t(xi));
}
REQUIRE(sum == Catch::Approx(231. / 2));
}
SECTION("Polynomial Gauss-Lobatto integral")
{
QuadratureFormula<3> qf = QuadratureManager::instance().getCubeFormula(GaussLobattoQuadratureDescriptor(4));
double sum = 0;
for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
const R3 xi = qf.point(i);
sum += qf.weight(i) * t.jacobianDeterminant(xi) * p(t(xi));
}
REQUIRE(sum == Catch::Approx(231. / 2));
}
SECTION("Polynomial Gauss-Legendre integral")
{
QuadratureFormula<3> qf = QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor(4));
double sum = 0;
for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
const R3 xi = qf.point(i);
sum += qf.weight(i) * t.jacobianDeterminant(xi) * p(t(xi));
}
REQUIRE(sum == Catch::Approx(231. / 2));
}
}
SECTION("degenerate to prism") {
const R3 a{1, 2, 0};
const R3 b{3, 1, 3};
const R3 c{2, 5, 2};
const R3 d{0, 3, 1};
const R3 e{1, 2, 5};
const R3 f{3, 1, 7};
const CubeTransformation t(a, b, c, c, d, e, f, f);
auto p = [](const R3& X) {
const double x = X[0];
const double y = X[1];
const double z = X[2];
return 3 * x * x + 2 * y * y + 3 * z * z + 4 * x + 3 * y + 2 * z + 1;
};
SECTION("Polynomial Gauss integral")
{
// 8 is the minimum quadrature rule to integrate the polynomial
// on the prism due to the jacobian determinant expression
const QuadratureFormula<3>& gauss = QuadratureManager::instance().getCubeFormula(GaussQuadratureDescriptor(8));
double integral = 0;
for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(t(gauss.point(i)));
}
REQUIRE(integral == Catch::Approx(30377. / 90));
}
SECTION("Polynomial Gauss-Legendre integral")
{
const QuadratureFormula<3>& gauss =
QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor(4));
double integral = 0;
for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(t(gauss.point(i)));
}
REQUIRE(integral == Catch::Approx(30377. / 90));
}
SECTION("Polynomial Gauss-Lobatto integral")
{
const QuadratureFormula<3>& gauss =
QuadratureManager::instance().getCubeFormula(GaussLobattoQuadratureDescriptor(4));
double integral = 0;
for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(t(gauss.point(i)));
}
REQUIRE(integral == Catch::Approx(30377. / 90));
}
}
SECTION("degenerate to pyramid")
{
const R3 a{1, 2, 0};
const R3 b{3, 1, 3};
const R3 c{2, 5, 2};
const R3 d{0, 3, 1};
const R3 e{1, 2, 5};
const CubeTransformation t(a, b, c, d, e, e, e, e);
auto p = [](const R3& X) { const double x = X[0];
const double y = X[1];
const double z = X[2];
return 3 * x * x + 2 * y * y + 3 * z * z + 4 * x + 3 * y + 2 * z + 1;
};
// 4 is the minimum quadrature rule to integrate the polynomial on the pyramid
const QuadratureFormula<3>& gauss =
QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor(20));
double integral = 0;
for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(t(gauss.point(i)));
}
REQUIRE(integral == Catch::Approx(7238. / 15));
}
}