From a0ef66438537149d9794f5a370d8e74e770af10f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Tue, 9 Nov 2021 19:10:37 +0100
Subject: [PATCH] Improve testing of Q1Transformation
These additional tests exhibited an issue wrt the Jacobian calculation
in the 3D case.
---
src/geometry/Q1Transformation.hpp | 14 ++--
tests/test_Q1Transformation.cpp | 124 +++++++++++++++++++++++-------
2 files changed, 105 insertions(+), 33 deletions(-)
diff --git a/src/geometry/Q1Transformation.hpp b/src/geometry/Q1Transformation.hpp
index bd23e5095..55e87bcbb 100644
--- a/src/geometry/Q1Transformation.hpp
+++ b/src/geometry/Q1Transformation.hpp
@@ -43,8 +43,11 @@ class Q1Transformation
const double& x = X[0];
const double& y = X[1];
- const TinyMatrix<Dimension, Dimension> J = {T(0, 0) + T(0, 2) * y, T(0, 1) + T(0, 2) * x, //
- T(1, 0) + T(1, 2) * y, T(1, 1) + T(1, 2) * x};
+ const TinyMatrix<Dimension, Dimension> J = {T(0, 0) + T(0, 2) * y, //
+ T(0, 1) + T(0, 2) * x,
+ //
+ T(1, 0) + T(1, 2) * y, //
+ T(1, 1) + T(1, 2) * x};
return det(J);
} else {
static_assert(Dimension == 3, "invalid dimension");
@@ -54,16 +57,17 @@ class Q1Transformation
const double& z = X[2];
const TinyMatrix<Dimension, Dimension> J = {T(0, 0) + T(0, 3) * y + T(0, 5) * z + T(0, 6) * y * z,
- T(0, 1) + T(0, 3) * x + T(0, 4) * z + T(0, 6) * x * y,
+ T(0, 1) + T(0, 3) * x + T(0, 4) * z + T(0, 6) * x * z,
T(0, 2) + T(0, 4) * y + T(0, 5) * x + T(0, 6) * x * y,
//
T(1, 0) + T(1, 3) * y + T(1, 5) * z + T(1, 6) * y * z,
- T(1, 1) + T(1, 3) * x + T(1, 4) * z + T(1, 6) * x * y,
+ T(1, 1) + T(1, 3) * x + T(1, 4) * z + T(1, 6) * x * z,
T(1, 2) + T(1, 4) * y + T(1, 5) * x + T(1, 6) * x * y,
//
T(2, 0) + T(2, 3) * y + T(2, 5) * z + T(2, 6) * y * z,
- T(2, 1) + T(2, 3) * x + T(2, 4) * z + T(2, 6) * x * y,
+ T(2, 1) + T(2, 3) * x + T(2, 4) * z + T(2, 6) * x * z,
T(2, 2) + T(2, 4) * y + T(2, 5) * x + T(2, 6) * x * y};
+
return det(J);
}
}
diff --git a/tests/test_Q1Transformation.cpp b/tests/test_Q1Transformation.cpp
index f42a60b4b..680b53e22 100644
--- a/tests/test_Q1Transformation.cpp
+++ b/tests/test_Q1Transformation.cpp
@@ -32,34 +32,60 @@ TEST_CASE("Q1Transformation", "[geometry]")
{
using R2 = TinyVector<2>;
+ const R2 a_hat = {-1, -1};
+ const R2 b_hat = {+1, -1};
+ const R2 c_hat = {+1, +1};
+ const R2 d_hat = {-1, +1};
+
+ const R2 m_hat = zero;
+
const R2 a = {0, 0};
const R2 b = {8, -2};
const R2 c = {12, 7};
const R2 d = {3, 7};
+ const R2 m = 0.25 * (a + b + c + d);
+
const std::array points = {a, b, c, d};
const Q1Transformation<2> t(points);
SECTION("values")
{
- REQUIRE(t({-1, -1})[0] == Catch::Approx(a[0]));
- REQUIRE(t({-1, -1})[1] == Catch::Approx(a[1]));
+ REQUIRE(t(a_hat)[0] == Catch::Approx(a[0]));
+ REQUIRE(t(a_hat)[1] == Catch::Approx(a[1]));
- REQUIRE(t({+1, -1})[0] == Catch::Approx(b[0]));
- REQUIRE(t({+1, -1})[1] == Catch::Approx(b[1]));
+ REQUIRE(t(b_hat)[0] == Catch::Approx(b[0]));
+ REQUIRE(t(b_hat)[1] == Catch::Approx(b[1]));
- REQUIRE(t({+1, +1})[0] == Catch::Approx(c[0]));
- REQUIRE(t({+1, +1})[1] == Catch::Approx(c[1]));
+ REQUIRE(t(c_hat)[0] == Catch::Approx(c[0]));
+ REQUIRE(t(c_hat)[1] == Catch::Approx(c[1]));
- REQUIRE(t({-1, +1})[0] == Catch::Approx(d[0]));
- REQUIRE(t({-1, +1})[1] == Catch::Approx(d[1]));
+ REQUIRE(t(d_hat)[0] == Catch::Approx(d[0]));
+ REQUIRE(t(d_hat)[1] == Catch::Approx(d[1]));
- REQUIRE(t({0, 0})[0] == Catch::Approx(0.25 * (a + b + c + d)[0]));
- REQUIRE(t({0, 0})[1] == Catch::Approx(0.25 * (a + b + c + d)[1]));
+ REQUIRE(t(m_hat)[0] == Catch::Approx(m[0]));
+ REQUIRE(t(m_hat)[1] == Catch::Approx(m[1]));
}
SECTION("Jacobian determinant")
{
+ SECTION("at points")
+ {
+ auto detJ = [](const R2& X) {
+ const double x = X[0];
+ const double y = X[1];
+
+ return 0.25 * ((0.5 * x + 4) * (y + 17) - 0.5 * (x + 7) * (y - 1));
+ };
+
+ 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(m_hat) == Catch::Approx(detJ(m_hat)));
+ }
+
SECTION("Gauss order 1")
{
// One point is enough in 2d
@@ -69,20 +95,28 @@ TEST_CASE("Q1Transformation", "[geometry]")
surface += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i));
}
- // 71.5 is actually the volume of the hexahedron
+ // 71.5 is actually the surface of the quadrangle
REQUIRE(surface == Catch::Approx(71.5));
}
- SECTION("Gauss order 7")
+ SECTION("Gauss order 3")
{
- TensorialGaussLegendreQuadrature<2> gauss(7);
- double surface = 0;
+ auto p = [](const R2& X) {
+ const double& x = X[0];
+ const double& y = X[1];
+
+ return (2 * x + 3 * y + 2) * (x - 2 * y - 1);
+ };
+
+ // Jacbian determinant is a degree 1 polynomial, so the
+ // following formula is required to reach exactness
+ TensorialGaussLegendreQuadrature<2> gauss(3);
+ double integral = 0;
for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
- surface += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i));
+ integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(gauss.point(i));
}
- // 71.5 is actually the volume of the hexahedron
- REQUIRE(surface == Catch::Approx(71.5));
+ REQUIRE(integral == Catch::Approx(-479. / 2));
}
}
}
@@ -102,7 +136,7 @@ TEST_CASE("Q1Transformation", "[geometry]")
const R3 m_hat = zero;
- const R3 a = {1, 2, 0};
+ const R3 a = {0, 0, 0};
const R3 b = {3, 1, 3};
const R3 c = {2, 5, 2};
const R3 d = {0, 3, 1};
@@ -132,29 +166,63 @@ TEST_CASE("Q1Transformation", "[geometry]")
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")
{
- // At least 2^3 points is are required in 2d
- TensorialGaussLegendreQuadrature<3> gauss(3);
+ // The jacobian determinant is of maximal degree 2 in each variable
+ TensorialGaussLegendreQuadrature<3> gauss(2);
double volume = 0;
for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
volume += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i));
}
- // 22.5 is actually the volume of the hexahedron
- REQUIRE(volume == Catch::Approx(22.5));
+ // 353. / 12 is actually the volume of the hexahedron
+ REQUIRE(volume == Catch::Approx(353. / 12));
}
- SECTION("Gauss order 7")
+ SECTION("Gauss order 4")
{
- TensorialGaussLegendreQuadrature<3> gauss(7);
- double volume = 0;
+ 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;
+ };
+
+ // Jacbian determinant is a degree 2 polynomial, so the
+ // following formula is required to reach exactness
+ TensorialGaussLegendreQuadrature<3> gauss(4);
+ double integral = 0;
for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
- volume += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i));
+ integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(gauss.point(i));
}
- // 22.5 is actually the volume of the hexahedron
- REQUIRE(volume == Catch::Approx(22.5));
+ REQUIRE(integral == Catch::Approx(8795. / 72));
}
}
}
--
GitLab