Add pyramid conform transformation and associated tests

src/geometry/PyramidTransformation.hpp

+#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

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

tests/test_PyramidTransformation.cpp

+#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));
+    }
+  }
+}