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

src/geometry/AffineTransformation.hpp

+#ifndef AFFINE_TRANSFORMATION_HPP
+#define AFFINE_TRANSFORMATION_HPP
+
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+
+#include <array>
+
+template <size_t Dimension>
+class AffineTransformation
+{
+ private:
+  constexpr static size_t NumberOfPoints = Dimension + 1;
+
+  TinyMatrix<Dimension> m_jacobian;
+  TinyVector<Dimension> m_shift;
+  double m_jacobian_determinant;
+
+ public:
+  PUGS_INLINE
+  TinyVector<Dimension>
+  operator()(const TinyVector<Dimension>& x) const
+  {
+    return m_jacobian * x + m_shift;
+  }
+
+  double
+  jacobianDeterminant() const
+  {
+    return m_jacobian_determinant;
+  }
+
+  PUGS_INLINE
+  AffineTransformation(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_jacobian(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];
+
+      for (size_t i = 0; i < Dimension; ++i) {
+        m_jacobian(i, 0) = 0.5 * (b[i] - a[i]);
+        m_jacobian(i, 1) = 0.5 * (c[i] - a[i]);
+
+        m_shift[i] = 0.5 * (b[i] + c[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];
+
+      for (size_t i = 0; i < Dimension; ++i) {
+        m_jacobian(i, 0) = 0.5 * (b[i] - a[i]);
+        m_jacobian(i, 1) = 0.5 * (c[i] - a[i]);
+        m_jacobian(i, 2) = 0.5 * (d[i] - a[i]);
+
+        m_shift[i] = 0.5 * (b[i] + c[i] + d[i] - a[i]);
+      }
+    }
+
+    m_jacobian_determinant = det(m_jacobian);
+  }
+
+  ~AffineTransformation() = default;
+};
+
+#endif   // AFFINE_TRANSFORMATION_HPP

tests/CMakeLists.txt

@@ -9,6 +9,7 @@ add_executable (unit_tests
   test_AffectationProcessor.cpp
   test_AffectationToStringProcessor.cpp
   test_AffectationToTupleProcessor.cpp
+  test_AffineTransformation.cpp
   test_Array.cpp
   test_ArraySubscriptProcessor.cpp
   test_ASTBuilder.cpp

tests/test_AffineTransformation.cpp

+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <geometry/AffineTransformation.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("AffineTransformation", "[geometry]")
+{
+  SECTION("1D")
+  {
+    using R1 = TinyVector<1>;
+
+    const R1 a = 0.3;
+    const R1 b = 2.7;
+
+    const std::array points = {a, b};
+
+    const AffineTransformation<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() == Catch::Approx(1.2));
+  }
+
+  SECTION("2D")
+  {
+    using R2 = TinyVector<2>;
+
+    const R2 a = {1, 2};
+    const R2 b = {3, 1};
+    const R2 c = {2, 5};
+
+    const std::array points = {a, b, c};
+    const AffineTransformation<2> t(points);
+
+    REQUIRE(t({-1, -1})[0] == Catch::Approx(1));
+    REQUIRE(t({-1, -1})[1] == Catch::Approx(2));
+
+    REQUIRE(t({+1, -1})[0] == Catch::Approx(3));
+    REQUIRE(t({+1, -1})[1] == Catch::Approx(1));
+
+    REQUIRE(t({-1, +1})[0] == Catch::Approx(2));
+    REQUIRE(t({-1, +1})[1] == Catch::Approx(5));
+
+    REQUIRE(t({-1. / 3, -1. / 3})[0] == Catch::Approx(2));
+    REQUIRE(t({-1. / 3, -1. / 3})[1] == Catch::Approx(8. / 3));
+
+    REQUIRE(t.jacobianDeterminant() == Catch::Approx(7. / 4));
+  }
+
+  SECTION("3D")
+  {
+    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 std::array points = {a, b, c, d};
+    const AffineTransformation<3> t(points);
+
+    REQUIRE(t({-1, -1, -1})[0] == Catch::Approx(1));
+    REQUIRE(t({-1, -1, -1})[1] == Catch::Approx(2));
+    REQUIRE(t({-1, -1, -1})[2] == Catch::Approx(1));
+
+    REQUIRE(t({+1, -1, -1})[0] == Catch::Approx(3));
+    REQUIRE(t({+1, -1, -1})[1] == Catch::Approx(1));
+    REQUIRE(t({+1, -1, -1})[2] == Catch::Approx(3));
+
+    REQUIRE(t({-1, +1, -1})[0] == Catch::Approx(2));
+    REQUIRE(t({-1, +1, -1})[1] == Catch::Approx(5));
+    REQUIRE(t({-1, +1, -1})[2] == Catch::Approx(2));
+
+    REQUIRE(t({-1, -1, +1})[0] == Catch::Approx(2));
+    REQUIRE(t({-1, -1, +1})[1] == Catch::Approx(3));
+    REQUIRE(t({-1, -1, +1})[2] == Catch::Approx(4));
+
+    REQUIRE(t({-0.5, -0.5, -0.5})[0] == Catch::Approx(2));
+    REQUIRE(t({-0.5, -0.5, -0.5})[1] == Catch::Approx(11. / 4));
+    REQUIRE(t({-0.5, -0.5, -0.5})[2] == Catch::Approx(2.5));
+
+    REQUIRE(t.jacobianDeterminant() == Catch::Approx(7. / 4));
+  }
+}