diff --git a/src/geometry/PyramidTransformation.hpp b/src/geometry/PyramidTransformation.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..19e7822ba3bc3e01da7be2840522e24707b68007
--- /dev/null
+++ b/src/geometry/PyramidTransformation.hpp
@@ -0,0 +1,72 @@
+#ifndef PYRAMID_TRANSFORMATION_HPP
+#define PYRAMID_TRANSFORMATION_HPP
+
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+
+#include <array>
+
+class PyramidTransformation
+{
+ private:
+ constexpr static size_t Dimension = 3;
+ constexpr static size_t NumberOfPoints = 5;
+
+ TinyMatrix<Dimension, NumberOfPoints - 1> m_coefficients;
+ TinyVector<Dimension> m_shift;
+
+ public:
+ PUGS_INLINE
+ TinyVector<Dimension>
+ operator()(const TinyVector<Dimension>& x) const
+ {
+ const TinyVector<NumberOfPoints - 1> X = {x[0], x[1], x[2], x[0] * x[1]};
+ return m_coefficients * X + m_shift;
+ }
+
+ double
+ jacobianDeterminant(const TinyVector<Dimension>& X) const
+ {
+ 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, 3) * y, //
+ T(0, 1) + T(0, 3) * x, //
+ T(0, 2),
+ //
+ T(1, 0) + T(1, 3) * y, //
+ T(1, 1) + T(1, 3) * x, //
+ T(1, 2),
+ //
+ T(2, 0) + T(2, 3) * y, //
+ T(2, 1) + T(2, 3) * x, //
+ T(2, 2)};
+ return det(J);
+ }
+
+ PUGS_INLINE
+ PyramidTransformation(const std::array<TinyVector<Dimension>, NumberOfPoints>& image_nodes)
+ {
+ static_assert(Dimension == 3, "invalid dimension");
+
+ 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];
+
+ m_shift = 0.25 * (a + b + c + d);
+
+ for (size_t i = 0; i < Dimension; ++i) {
+ m_coefficients(i, 0) = 0.5 * (b[i] + c[i]) - m_shift[i];
+ m_coefficients(i, 1) = 0.5 * (c[i] + d[i]) - m_shift[i];
+ m_coefficients(i, 2) = e[i] - m_shift[i];
+ m_coefficients(i, 3) = 0.5 * (a[i] + c[i]) - m_shift[i];
+ }
+ }
+
+ ~PyramidTransformation() = default;
+};
+
+#endif // PYRAMID_TRANSFORMATION_HPP
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index fadd4ec5d1138d1423c4e758e9cc59249684c8e7..ba8ceb6c279f71dddff992aa4aebac11a629f79c 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -101,6 +101,7 @@ add_executable (unit_tests
test_PugsFunctionAdapter.cpp
test_PugsUtils.cpp
test_PyramidGaussQuadrature.cpp
+ test_PyramidTransformation.cpp
test_Q1Transformation.cpp
test_RevisionInfo.cpp
test_SmallArray.cpp
diff --git a/tests/test_PyramidTransformation.cpp b/tests/test_PyramidTransformation.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..bfdcee769d97558a2ff8037752b8e060fb3a8c8b
--- /dev/null
+++ b/tests/test_PyramidTransformation.cpp
@@ -0,0 +1,95 @@
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <analysis/PyramidGaussQuadrature.hpp>
+#include <geometry/PyramidTransformation.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("PyramidTransformation", "[geometry]")
+{
+ using R3 = TinyVector<3>;
+
+ const R3 a_hat = {-1, -1, +0};
+ const R3 b_hat = {+1, -1, +0};
+ const R3 c_hat = {+1, +1, +0};
+ const R3 d_hat = {-1, +1, +0};
+ const R3 e_hat = {+0, +0, +1};
+
+ const R3 m_hat = {0, 0, 1. / 5};
+
+ 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 std::array points = {a, b, c, d, e};
+ const PyramidTransformation t(points);
+
+ SECTION("points")
+ {
+ 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));
+
+ R3 m = (1. / 5) * (a + b + c + d + e);
+ REQUIRE(l2Norm(t(m_hat) - m) == Catch::Approx(0).margin(1E-14));
+ }
+
+ SECTION("Jacobian determinant")
+ {
+ SECTION("at points")
+ {
+ auto detJ = [](const R3 X) {
+ const double& x = X[0];
+ const double& y = X[1];
+
+ return (43 * x + 13 * y + 93) / 16;
+ };
+
+ 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(e_hat) == Catch::Approx(detJ(e_hat)));
+
+ REQUIRE(t.jacobianDeterminant(m_hat) == Catch::Approx(detJ(m_hat)));
+ }
+
+ SECTION("volume calculation")
+ {
+ // The jacobian determinant is a degree 1 polynomial
+ PyramidGaussQuadrature gauss(1);
+ double volume = 0;
+ for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
+ volume += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i));
+ }
+
+ // 31 / 4 is actually the volume of the pyramid
+ REQUIRE(volume == Catch::Approx(31. / 4));
+ }
+
+ SECTION("exact polynomial integration")
+ {
+ auto p = [](const R3& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+
+ return 3 * x * x + 2 * y * y + 3 * z * z + 4 * x + 3 * y + 2 * z + 1;
+ };
+
+ // 3 is the minimum quadrature rule to integrate the polynomial on the pyramid
+ PyramidGaussQuadrature gauss(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(1513. / 60));
+ }
+ }
+}