Add Q1 transformation (R^d -> R^d)

src/geometry/Q1Transformation.hpp

+#ifndef Q1_TRANSFORMATION_HPP
+#define Q1_TRANSFORMATION_HPP
+
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+
+#include <array>
+
+template <size_t Dimension>
+class Q1Transformation
+{
+ private:
+  constexpr static size_t NumberOfPoints = std::pow(2, Dimension);
+
+  TinyMatrix<Dimension, NumberOfPoints - 1> m_coefficients;
+  TinyVector<Dimension> m_shift;
+
+ public:
+  PUGS_INLINE
+  TinyVector<Dimension>
+  operator()(const TinyVector<Dimension>& x) const
+  {
+    if constexpr (Dimension == 1) {
+      return m_coefficients * x + m_shift;
+    } else if constexpr (Dimension == 2) {
+      const TinyVector<NumberOfPoints - 1> X = {x[0], x[1], x[0] * x[1]};
+      return m_coefficients * X + m_shift;
+    } else {
+      static_assert(Dimension == 3, "invalid dimension");
+      const TinyVector<NumberOfPoints - 1> X =
+        {x[0], x[1], x[2], x[0] * x[1], x[1] * x[2], x[0] * x[2], x[0] * x[1] * x[2]};
+      return m_coefficients * X + m_shift;
+    }
+  }
+
+  PUGS_INLINE double
+  jacobianDeterminant([[maybe_unused]] const TinyVector<Dimension>& X) const
+  {
+    if constexpr (Dimension == 1) {
+      return m_coefficients(0, 0);
+    } else if constexpr (Dimension == 2) {
+      const auto& T   = m_coefficients;
+      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};
+      return det(J);
+    } else {
+      static_assert(Dimension == 3, "invalid dimension");
+      const auto& T   = m_coefficients;
+      const double& x = X[0];
+      const double& y = X[1];
+      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, 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, 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, 2) + T(2, 4) * y + T(2, 5) * x + T(2, 6) * x * y};
+      return det(J);
+    }
+  }
+
+  PUGS_INLINE
+  Q1Transformation(const std::array<TinyVector<Dimension>, NumberOfPoints>& image_nodes)
+  {
+    static_assert(Dimension >= 1 and Dimension <= 3, "invalid dimension");
+    if constexpr (Dimension == 1) {
+      const TinyVector<Dimension>& a = image_nodes[0];
+      const TinyVector<Dimension>& b = image_nodes[1];
+
+      m_coefficients(0, 0) = 0.5 * (b[0] - a[0]);
+
+      m_shift[0] = 0.5 * (a[0] + b[0]);
+    } else if constexpr (Dimension == 2) {
+      const TinyVector<Dimension>& a = image_nodes[0];
+      const TinyVector<Dimension>& b = image_nodes[1];
+      const TinyVector<Dimension>& c = image_nodes[2];
+      const TinyVector<Dimension>& d = image_nodes[3];
+
+      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]);
+      }
+    } else {
+      static_assert(Dimension == 3);
+
+      const TinyVector<Dimension>& a = image_nodes[0];
+      const TinyVector<Dimension>& b = image_nodes[1];
+      const TinyVector<Dimension>& c = image_nodes[2];
+      const TinyVector<Dimension>& d = image_nodes[3];
+      const TinyVector<Dimension>& e = image_nodes[4];
+      const TinyVector<Dimension>& f = image_nodes[5];
+      const TinyVector<Dimension>& g = image_nodes[6];
+      const TinyVector<Dimension>& h = image_nodes[7];
+
+      for (size_t i = 0; i < Dimension; ++i) {
+        m_coefficients(i, 0) = 0.125 * (-a[i] + b[i] + c[i] - d[i] - e[i] + f[i] + g[i] - h[i]);
+        m_coefficients(i, 1) = 0.125 * (-a[i] - b[i] + c[i] + d[i] - e[i] - f[i] + g[i] + h[i]);
+        m_coefficients(i, 2) = 0.125 * (-a[i] - b[i] - c[i] - d[i] + e[i] + f[i] + g[i] + h[i]);
+        m_coefficients(i, 3) = 0.125 * (+a[i] - b[i] + c[i] - d[i] + e[i] - f[i] + g[i] - h[i]);
+        m_coefficients(i, 4) = 0.125 * (+a[i] + b[i] - c[i] - d[i] - e[i] - f[i] + g[i] + h[i]);
+        m_coefficients(i, 5) = 0.125 * (+a[i] - b[i] - c[i] + d[i] - e[i] + f[i] + g[i] - h[i]);
+        m_coefficients(i, 6) = 0.125 * (-a[i] + b[i] - c[i] + d[i] + e[i] - f[i] + g[i] - h[i]);
+
+        m_shift[i] = 0.125 * (a[i] + b[i] + c[i] + d[i] + e[i] + f[i] + g[i] + h[i]);
+      }
+    }
+  }
+
+  ~Q1Transformation() = default;
+};
+
+#endif   // Q1_TRANSFORMATION_HPP

tests/CMakeLists.txt

@@ -99,6 +99,7 @@ add_executable (unit_tests
   test_PugsAssert.cpp
   test_PugsFunctionAdapter.cpp
   test_PugsUtils.cpp
+  test_Q1Transformation.cpp
   test_RevisionInfo.cpp
   test_SmallArray.cpp
   test_SmallVector.cpp

tests/test_Q1Transformation.cpp

+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <analysis/GaussLegendreQuadrature.hpp>
+#include <geometry/Q1Transformation.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("Q1Transformation", "[geometry]")
+{
+  SECTION("1D")
+  {
+    using R1 = TinyVector<1>;
+
+    const R1 a = 0.3;
+    const R1 b = 2.7;
+
+    const std::array points = {a, b};
+
+    const Q1Transformation<1> t(points);
+
+    REQUIRE(t(-1)[0] == Catch::Approx(0.3));
+    REQUIRE(t(0)[0] == Catch::Approx(1.5));
+    REQUIRE(t(1)[0] == Catch::Approx(2.7));
+
+    REQUIRE(t.jacobianDeterminant(R1{-1}) == Catch::Approx(1.2));
+    REQUIRE(t.jacobianDeterminant(R1{0}) == Catch::Approx(1.2));
+    REQUIRE(t.jacobianDeterminant(R1{3}) == Catch::Approx(1.2));
+  }
+
+  SECTION("2D")
+  {
+    using R2 = TinyVector<2>;
+
+    const R2 a = {0, 0};
+    const R2 b = {8, -2};
+    const R2 c = {12, 7};
+    const R2 d = {3, 7};
+
+    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({+1, -1})[0] == Catch::Approx(b[0]));
+      REQUIRE(t({+1, -1})[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({-1, +1})[0] == Catch::Approx(d[0]));
+      REQUIRE(t({-1, +1})[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]));
+    }
+
+    SECTION("Jacobian determinant")
+    {
+      SECTION("Gauss order 1")
+      {
+        // One point is enough in 2d
+        GaussLegendreQuadrature<2> gauss(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 volume of the hexahedron
+        REQUIRE(surface == Catch::Approx(71.5));
+      }
+
+      SECTION("Gauss order 7")
+      {
+        GaussLegendreQuadrature<2> gauss(7);
+        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 volume of the hexahedron
+        REQUIRE(surface == Catch::Approx(71.5));
+      }
+    }
+  }
+
+  SECTION("3D")
+  {
+    using R3 = TinyVector<3>;
+
+    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 = {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 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 std::array points = {a, b, c, d, e, f, g, h};
+    const Q1Transformation<3> t(points);
+
+    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("Gauss order 3")
+      {
+        // At least 2^3 points is are required in 2d
+        GaussLegendreQuadrature<3> gauss(3);
+        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));
+      }
+
+      SECTION("Gauss order 7")
+      {
+        GaussLegendreQuadrature<3> gauss(7);
+        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));
+      }
+    }
+  }
+}