diff --git a/src/geometry/Q1Transformation.hpp b/src/geometry/Q1Transformation.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..f31e67a2e1dfd74e70d49b3c6f47b020fd1dad38
--- /dev/null
+++ b/src/geometry/Q1Transformation.hpp
@@ -0,0 +1,124 @@
+#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
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 5d0fa65b16e918cbdb5e0d6e14351709c3c6bba9..69542a25317c0cbb4ea7f64964697cb948af042d 100644
--- a/tests/CMakeLists.txt
+++ b/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
diff --git a/tests/test_Q1Transformation.cpp b/tests/test_Q1Transformation.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b4e0314ef7c6c945e477e9bc426eeef1dafdc690
--- /dev/null
+++ b/tests/test_Q1Transformation.cpp
@@ -0,0 +1,161 @@
+#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));
+ }
+ }
+ }
+}