From 8fba3d253db76f1c37db70aa2d9b0167b96f623f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Thu, 7 Oct 2021 18:53:23 +0200
Subject: [PATCH] Add affine transformation (R^d -> R^d)
---
src/geometry/AffineTransformation.hpp | 78 ++++++++++++++++++++++++
tests/CMakeLists.txt | 1 +
tests/test_AffineTransformation.cpp | 88 +++++++++++++++++++++++++++
3 files changed, 167 insertions(+)
create mode 100644 src/geometry/AffineTransformation.hpp
create mode 100644 tests/test_AffineTransformation.cpp
diff --git a/src/geometry/AffineTransformation.hpp b/src/geometry/AffineTransformation.hpp
new file mode 100644
index 000000000..e7179ef5d
--- /dev/null
+++ b/src/geometry/AffineTransformation.hpp
@@ -0,0 +1,78 @@
+#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
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index b5d57ecad..5d0fa65b1 100644
--- a/tests/CMakeLists.txt
+++ b/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
diff --git a/tests/test_AffineTransformation.cpp b/tests/test_AffineTransformation.cpp
new file mode 100644
index 000000000..cce135741
--- /dev/null
+++ b/tests/test_AffineTransformation.cpp
@@ -0,0 +1,88 @@
+#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));
+ }
+}
--
GitLab