Skip to content
Snippets Groups Projects
Commit a0ef6643 authored by Stéphane Del Pino's avatar Stéphane Del Pino
Browse files

Improve testing of Q1Transformation

These additional tests exhibited an issue wrt the Jacobian calculation
in the 3D case.
parent 48e599b9
No related branches found
No related tags found
1 merge request!124Add files for high order integration with quadratures
...@@ -43,8 +43,11 @@ class Q1Transformation ...@@ -43,8 +43,11 @@ class Q1Transformation
const double& x = X[0]; const double& x = X[0];
const double& y = X[1]; const double& y = X[1];
const TinyMatrix<Dimension, Dimension> J = {T(0, 0) + T(0, 2) * y, T(0, 1) + T(0, 2) * x, // const TinyMatrix<Dimension, Dimension> J = {T(0, 0) + T(0, 2) * y, //
T(1, 0) + T(1, 2) * y, T(1, 1) + T(1, 2) * x}; T(0, 1) + T(0, 2) * x,
//
T(1, 0) + T(1, 2) * y, //
T(1, 1) + T(1, 2) * x};
return det(J); return det(J);
} else { } else {
static_assert(Dimension == 3, "invalid dimension"); static_assert(Dimension == 3, "invalid dimension");
...@@ -54,16 +57,17 @@ class Q1Transformation ...@@ -54,16 +57,17 @@ class Q1Transformation
const double& z = X[2]; 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, 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(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, 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(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, 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}; T(2, 2) + T(2, 4) * y + T(2, 5) * x + T(2, 6) * x * y};
return det(J); return det(J);
} }
} }
......
...@@ -32,34 +32,60 @@ TEST_CASE("Q1Transformation", "[geometry]") ...@@ -32,34 +32,60 @@ TEST_CASE("Q1Transformation", "[geometry]")
{ {
using R2 = TinyVector<2>; 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 a = {0, 0};
const R2 b = {8, -2}; const R2 b = {8, -2};
const R2 c = {12, 7}; const R2 c = {12, 7};
const R2 d = {3, 7}; const R2 d = {3, 7};
const R2 m = 0.25 * (a + b + c + d);
const std::array points = {a, b, c, d}; const std::array points = {a, b, c, d};
const Q1Transformation<2> t(points); const Q1Transformation<2> t(points);
SECTION("values") SECTION("values")
{ {
REQUIRE(t({-1, -1})[0] == Catch::Approx(a[0])); REQUIRE(t(a_hat)[0] == Catch::Approx(a[0]));
REQUIRE(t({-1, -1})[1] == Catch::Approx(a[1])); REQUIRE(t(a_hat)[1] == Catch::Approx(a[1]));
REQUIRE(t({+1, -1})[0] == Catch::Approx(b[0])); REQUIRE(t(b_hat)[0] == Catch::Approx(b[0]));
REQUIRE(t({+1, -1})[1] == Catch::Approx(b[1])); REQUIRE(t(b_hat)[1] == Catch::Approx(b[1]));
REQUIRE(t({+1, +1})[0] == Catch::Approx(c[0])); REQUIRE(t(c_hat)[0] == Catch::Approx(c[0]));
REQUIRE(t({+1, +1})[1] == Catch::Approx(c[1])); REQUIRE(t(c_hat)[1] == Catch::Approx(c[1]));
REQUIRE(t({-1, +1})[0] == Catch::Approx(d[0])); REQUIRE(t(d_hat)[0] == Catch::Approx(d[0]));
REQUIRE(t({-1, +1})[1] == Catch::Approx(d[1])); 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(m_hat)[0] == Catch::Approx(m[0]));
REQUIRE(t({0, 0})[1] == Catch::Approx(0.25 * (a + b + c + d)[1])); REQUIRE(t(m_hat)[1] == Catch::Approx(m[1]));
} }
SECTION("Jacobian determinant") 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") SECTION("Gauss order 1")
{ {
// One point is enough in 2d // One point is enough in 2d
...@@ -69,20 +95,28 @@ TEST_CASE("Q1Transformation", "[geometry]") ...@@ -69,20 +95,28 @@ TEST_CASE("Q1Transformation", "[geometry]")
surface += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)); 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)); REQUIRE(surface == Catch::Approx(71.5));
} }
SECTION("Gauss order 7") SECTION("Gauss order 3")
{ {
TensorialGaussLegendreQuadrature<2> gauss(7); auto p = [](const R2& X) {
double surface = 0; 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) { 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(integral == Catch::Approx(-479. / 2));
REQUIRE(surface == Catch::Approx(71.5));
} }
} }
} }
...@@ -102,7 +136,7 @@ TEST_CASE("Q1Transformation", "[geometry]") ...@@ -102,7 +136,7 @@ TEST_CASE("Q1Transformation", "[geometry]")
const R3 m_hat = zero; 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 b = {3, 1, 3};
const R3 c = {2, 5, 2}; const R3 c = {2, 5, 2};
const R3 d = {0, 3, 1}; const R3 d = {0, 3, 1};
...@@ -132,29 +166,63 @@ TEST_CASE("Q1Transformation", "[geometry]") ...@@ -132,29 +166,63 @@ TEST_CASE("Q1Transformation", "[geometry]")
SECTION("Jacobian determinant") 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") SECTION("Gauss order 3")
{ {
// At least 2^3 points is are required in 2d // The jacobian determinant is of maximal degree 2 in each variable
TensorialGaussLegendreQuadrature<3> gauss(3); TensorialGaussLegendreQuadrature<3> gauss(2);
double volume = 0; double volume = 0;
for (size_t i = 0; i < gauss.numberOfPoints(); ++i) { for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
volume += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)); volume += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i));
} }
// 22.5 is actually the volume of the hexahedron // 353. / 12 is actually the volume of the hexahedron
REQUIRE(volume == Catch::Approx(22.5)); REQUIRE(volume == Catch::Approx(353. / 12));
} }
SECTION("Gauss order 7") SECTION("Gauss order 4")
{ {
TensorialGaussLegendreQuadrature<3> gauss(7); auto p = [](const R3& X) {
double volume = 0; 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) { 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(integral == Catch::Approx(8795. / 72));
REQUIRE(volume == Catch::Approx(22.5));
} }
} }
} }
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment