diff --git a/src/geometry/SquareTransformation.hpp b/src/geometry/SquareTransformation.hpp
index 4d2906bb514fe90137f7270404c2f5b5f54d24a2..d79ae12c96c696f6be0bcf6245ede0b16c5f2a83 100644
--- a/src/geometry/SquareTransformation.hpp
+++ b/src/geometry/SquareTransformation.hpp
@@ -6,7 +6,11 @@
#include <array>
-class SquareTransformation
+template <size_t GivenDimension>
+class SquareTransformation;
+
+template <>
+class SquareTransformation<2>
{
public:
constexpr static size_t Dimension = 2;
@@ -58,4 +62,62 @@ class SquareTransformation
~SquareTransformation() = default;
};
+template <size_t GivenDimension>
+class SquareTransformation
+{
+ public:
+ constexpr static size_t Dimension = GivenDimension;
+ static_assert(Dimension == 3, "Square transformation is defined in dimension 2 or 3");
+
+ constexpr static size_t NumberOfPoints = 4;
+
+ private:
+ TinyMatrix<Dimension, NumberOfPoints - 1> m_coefficients;
+ TinyVector<Dimension> m_shift;
+
+ public:
+ PUGS_INLINE
+ TinyVector<Dimension>
+ operator()(const TinyVector<2>& x) const
+ {
+ const TinyVector<NumberOfPoints - 1> X = {x[0], x[1], x[0] * x[1]};
+ return m_coefficients * X + m_shift;
+ }
+
+ PUGS_INLINE double
+ areaVariationNorm(const TinyVector<2>& X) const
+ {
+ const auto& T = m_coefficients;
+ const double& x = X[0];
+ const double& y = X[1];
+
+ const TinyVector<Dimension> dxT = {T(0, 0) + T(0, 2) * y, //
+ T(1, 0) + T(1, 2) * y, //
+ T(2, 0) + T(2, 2) * y};
+
+ const TinyVector<Dimension> dyT = {T(0, 1) + T(0, 2) * x, //
+ T(1, 1) + T(1, 2) * x, //
+ T(2, 1) + T(2, 2) * x};
+
+ return l2Norm(crossProduct(dxT, dyT));
+ }
+
+ PUGS_INLINE
+ SquareTransformation(const TinyVector<Dimension>& a,
+ const TinyVector<Dimension>& b,
+ const TinyVector<Dimension>& c,
+ const TinyVector<Dimension>& d)
+ {
+ for (size_t i = 0; i < Dimension; ++i) {
+ m_coefficients(i, 0) = 0.25 * (-a[i] + b[i] + c[i] - d[i]);
+ m_coefficients(i, 1) = 0.25 * (-a[i] - b[i] + c[i] + d[i]);
+ m_coefficients(i, 2) = 0.25 * (+a[i] - b[i] + c[i] - d[i]);
+
+ m_shift[i] = 0.25 * (a[i] + b[i] + c[i] + d[i]);
+ }
+ }
+
+ ~SquareTransformation() = default;
+};
+
#endif // SQUARE_TRANSFORMATION_HPP
diff --git a/src/language/utils/IntegrateOnCells.hpp b/src/language/utils/IntegrateOnCells.hpp
index 216483dc42530998aefbe3418aaad3b7d634be63..669190ca6b31895a66165671e08a192fcbda4b3f 100644
--- a/src/language/utils/IntegrateOnCells.hpp
+++ b/src/language/utils/IntegrateOnCells.hpp
@@ -158,8 +158,8 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
} else if constexpr (Dimension == 2) {
switch (cell_type[cell_id]) {
case CellType::Triangle: {
- const SquareTransformation T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]],
- xr[cell_to_node[2]]);
+ const SquareTransformation<2> T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]],
+ xr[cell_to_node[2]]);
for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
const auto xi = qf.point(i_point);
@@ -170,8 +170,8 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
break;
}
case CellType::Quadrangle: {
- const SquareTransformation T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]],
- xr[cell_to_node[3]]);
+ const SquareTransformation<2> T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]],
+ xr[cell_to_node[3]]);
for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
const auto xi = qf.point(i_point);
@@ -411,8 +411,8 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
break;
}
case CellType::Quadrangle: {
- const SquareTransformation T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]],
- xr[cell_to_node[3]]);
+ const SquareTransformation<2> T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]],
+ xr[cell_to_node[3]]);
const auto qf = QuadratureManager::instance().getSquareFormula(quadrature_descriptor);
for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
diff --git a/tests/test_SquareTransformation.cpp b/tests/test_SquareTransformation.cpp
index 19523959948a9ae3980e5df9f6efb4dcb0e12734..fa12d756f27945d0189f19785c7c5bc1caccb0d5 100644
--- a/tests/test_SquareTransformation.cpp
+++ b/tests/test_SquareTransformation.cpp
@@ -9,96 +9,204 @@
TEST_CASE("SquareTransformation", "[geometry]")
{
- using R2 = TinyVector<2>;
+ SECTION("2D")
+ {
+ 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 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 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 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 R2 m = 0.25 * (a + b + c + d);
- const SquareTransformation t(a, b, c, d);
+ const SquareTransformation<2> t(a, b, c, d);
- SECTION("values")
- {
- REQUIRE(t(a_hat)[0] == Catch::Approx(a[0]));
- REQUIRE(t(a_hat)[1] == Catch::Approx(a[1]));
+ SECTION("values")
+ {
+ REQUIRE(t(a_hat)[0] == Catch::Approx(a[0]));
+ REQUIRE(t(a_hat)[1] == Catch::Approx(a[1]));
- REQUIRE(t(b_hat)[0] == Catch::Approx(b[0]));
- REQUIRE(t(b_hat)[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(c_hat)[0] == Catch::Approx(c[0]));
- REQUIRE(t(c_hat)[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(d_hat)[0] == Catch::Approx(d[0]));
- REQUIRE(t(d_hat)[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(m_hat)[0] == Catch::Approx(m[0]));
- REQUIRE(t(m_hat)[1] == Catch::Approx(m[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")
+ SECTION("Jacobian determinant")
{
- auto detJ = [](const R2& X) {
- const double x = X[0];
- const double y = X[1];
+ 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));
- };
+ 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(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
+ const QuadratureFormula<2>& gauss =
+ QuadratureManager::instance().getSquareFormula(GaussLegendreQuadratureDescriptor(1));
+
+ double surface = 0;
+ for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
+ surface += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i));
+ }
+
+ // 71.5 is actually the surface of the quadrangle
+ REQUIRE(surface == Catch::Approx(71.5));
+ }
- REQUIRE(t.jacobianDeterminant(m_hat) == Catch::Approx(detJ(m_hat)));
+ SECTION("Gauss order 3")
+ {
+ 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
+ const QuadratureFormula<2>& gauss =
+ QuadratureManager::instance().getSquareFormula(GaussLegendreQuadratureDescriptor(3));
+
+ double integral = 0;
+ for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
+ integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(gauss.point(i));
+ }
+
+ REQUIRE(integral == Catch::Approx(-479. / 2));
+ }
}
+ }
+
+ SECTION("3D")
+ {
+ 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;
- SECTION("Gauss order 1")
+ using R3 = TinyVector<3>;
+
+ const R3 a = {0, 0, -1};
+ const R3 b = {8, -2, 3};
+ const R3 c = {12, 7, 2};
+ const R3 d = {3, 7, 1};
+
+ const R3 m = 0.25 * (a + b + c + d);
+
+ const SquareTransformation<3> t(a, b, c, d);
+
+ SECTION("values")
{
- // One point is enough in 2d
- const QuadratureFormula<2>& gauss =
- QuadratureManager::instance().getSquareFormula(GaussLegendreQuadratureDescriptor(1));
+ REQUIRE(t(a_hat)[0] == Catch::Approx(a[0]));
+ REQUIRE(t(a_hat)[1] == Catch::Approx(a[1]));
+ REQUIRE(t(a_hat)[2] == Catch::Approx(a[2]));
- double surface = 0;
- for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
- surface += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i));
- }
+ REQUIRE(t(b_hat)[0] == Catch::Approx(b[0]));
+ REQUIRE(t(b_hat)[1] == Catch::Approx(b[1]));
+ REQUIRE(t(b_hat)[2] == Catch::Approx(b[2]));
- // 71.5 is actually the surface of the quadrangle
- REQUIRE(surface == Catch::Approx(71.5));
+ REQUIRE(t(c_hat)[0] == Catch::Approx(c[0]));
+ REQUIRE(t(c_hat)[1] == Catch::Approx(c[1]));
+ REQUIRE(t(c_hat)[2] == Catch::Approx(c[2]));
+
+ REQUIRE(t(d_hat)[0] == Catch::Approx(d[0]));
+ REQUIRE(t(d_hat)[1] == Catch::Approx(d[1]));
+ REQUIRE(t(d_hat)[2] == Catch::Approx(d[2]));
+
+ REQUIRE(t(m_hat)[0] == Catch::Approx(m[0]));
+ REQUIRE(t(m_hat)[1] == Catch::Approx(m[1]));
+ REQUIRE(t(m_hat)[2] == Catch::Approx(m[2]));
}
- SECTION("Gauss order 3")
+ SECTION("Area variation norm")
{
- auto p = [](const R2& X) {
- const double& x = X[0];
- const double& y = X[1];
+ auto area_variation_norm = [&](const R2& X) {
+ const double x = X[0];
+ const double y = X[1];
+
+ const R3 J1 = 0.25 * (-a + b + c - d);
+ const R3 J2 = 0.25 * (-a - b + c + d);
+ const R3 J3 = 0.25 * (a - b + c - d);
- return (2 * x + 3 * y + 2) * (x - 2 * y - 1);
+ return l2Norm(crossProduct(J1 + y * J3, J2 + x * J3));
};
- // Jacbian determinant is a degree 1 polynomial, so the
- // following formula is required to reach exactness
- const QuadratureFormula<2>& gauss =
- QuadratureManager::instance().getSquareFormula(GaussLegendreQuadratureDescriptor(3));
+ SECTION("at points")
+ {
+ REQUIRE(t.areaVariationNorm(a_hat) == Catch::Approx(area_variation_norm(a_hat)));
+ REQUIRE(t.areaVariationNorm(b_hat) == Catch::Approx(area_variation_norm(b_hat)));
+ REQUIRE(t.areaVariationNorm(c_hat) == Catch::Approx(area_variation_norm(c_hat)));
+ REQUIRE(t.areaVariationNorm(d_hat) == Catch::Approx(area_variation_norm(d_hat)));
- double integral = 0;
- for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
- integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(gauss.point(i));
+ REQUIRE(t.areaVariationNorm(m_hat) == Catch::Approx(area_variation_norm(m_hat)));
}
- REQUIRE(integral == Catch::Approx(-479. / 2));
+ 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 - 3 * y * y + y + 2 * z * z - 0.5 * z + 2;
+ };
+
+ SECTION("Gauss order 3")
+ {
+ // Due to the area variation term (square root of a
+ // polynomial), Gauss-like quadrature cannot be exact for
+ // general 3d quadrangle.
+ //
+ // Moreover, even the surface of general 3d quadrangle cannot
+ // be computed exactly.
+ //
+ // So for this test we integrate the following function:
+ // f(x,y,z) := p(x,y,z) * |dA(t^-1(x,y,z))|
+ //
+ // dA denotes the area change on the reference element. This
+ // function can be exactly integrated since computing the
+ // integral on the reference quadrangle, one gets a dA^2,
+ // which is a polynomial.
+ const QuadratureFormula<2>& gauss =
+ QuadratureManager::instance().getSquareFormula(GaussLegendreQuadratureDescriptor(4));
+
+ double integral = 0;
+ for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
+ const double fX = (t.areaVariationNorm(gauss.point(i)) * p(t(gauss.point(i))));
+ integral += gauss.weight(i) * t.areaVariationNorm(gauss.point(i)) * fX;
+ }
+
+ REQUIRE(integral == Catch::Approx(60519715. / 576));
+ }
}
}
}