Add Prism conform transformation and the associated tests

src/geometry/PrismTransformation.hpp

+#ifndef PRISM_TRANSFORMATION_HPP
+#define PRISM_TRANSFORMATION_HPP
+
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+
+#include <array>
+
+class PrismTransformation
+{
+ private:
+  constexpr static size_t Dimension      = 3;
+  constexpr static size_t NumberOfPoints = 6;
+
+  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[2], x[1] * x[2]};
+    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 double& z = X[2];
+
+    const TinyMatrix<Dimension, Dimension> J = {T(0, 0) + T(0, 3) * z,   //
+                                                T(0, 1) + T(0, 4) * z,   //
+                                                T(0, 2) + T(0, 3) * x + T(0, 4) * y,
+                                                //
+                                                T(1, 0) + T(1, 3) * z,   //
+                                                T(1, 1) + T(1, 4) * z,   //
+                                                T(1, 2) + T(1, 3) * x + T(1, 4) * y,
+                                                //
+                                                T(2, 0) + T(2, 3) * z,   //
+                                                T(2, 1) + T(2, 4) * z,   //
+                                                T(2, 2) + T(2, 3) * x + T(2, 4) * y};
+    return det(J);
+  }
+
+  PUGS_INLINE
+  PrismTransformation(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];
+    const TinyVector<Dimension>& f = image_nodes[5];
+
+    for (size_t i = 0; i < Dimension; ++i) {
+      m_coefficients(i, 0) = 0.5 * (b[i] - a[i] + e[i] - d[i]);
+      m_coefficients(i, 1) = 0.5 * (c[i] + f[i] - a[i] - d[i]);
+      m_coefficients(i, 2) = 0.5 * (d[i] - a[i]);
+      m_coefficients(i, 3) = 0.5 * (a[i] - b[i] + e[i] - d[i]);
+      m_coefficients(i, 4) = 0.5 * (f[i] - c[i] + a[i] - d[i]);
+    }
+
+    m_shift = 0.5 * (a + d);
+  }
+
+  ~PrismTransformation() = default;
+};
+
+#endif   // PRISM_TRANSFORMATION_HPP

tests/CMakeLists.txt

@@ -96,6 +96,7 @@ add_executable (unit_tests
   test_ParseError.cpp
   test_PETScUtils.cpp
   test_PrismGaussQuadrature.cpp
+  test_PrismTransformation.cpp
   test_PugsAssert.cpp
   test_PugsFunctionAdapter.cpp
   test_PugsUtils.cpp

tests/test_PrismTransformation.cpp

+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <analysis/PrismGaussQuadrature.hpp>
+#include <geometry/PrismTransformation.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("PrismTransformation", "[geometry]")
+{
+  using R3 = TinyVector<3>;
+
+  const R3 a_hat = {0, 0, -1};
+  const R3 b_hat = {1, 0, -1};
+  const R3 c_hat = {0, 1, -1};
+  const R3 d_hat = {0, 0, +1};
+  const R3 e_hat = {1, 0, +1};
+  const R3 f_hat = {0, 1, +1};
+
+  const R3 m_hat = {1. / 3, 1. / 3, 0};
+
+  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 std::array points = {a, b, c, d, e, f};
+  const PrismTransformation 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));
+    REQUIRE(l2Norm(t(f_hat) - f) == Catch::Approx(0));
+
+    R3 m = (1. / 6) * (a + b + c + d + e + f);
+
+    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];
+        const double& z = X[2];
+
+        return ((2 * y + 0.5 * x + 0.5) * (z + 2) - (y - 0.5 * x - 0.5) * (2 * z + 4)) +
+               (1.5 - 0.5 * z) * ((2 * y + 0.5 * x + 0.5) * (0.5 - 2.5 * z) - (0.5 - 2.5 * y) * (2 * z + 4)) +
+               (0.5 * z + 3.5) * ((0.5 - 2.5 * y) * (z + 2) - (y - 0.5 * x - 0.5) * (0.5 - 2.5 * z));
+      };
+
+      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(f_hat) == Catch::Approx(detJ(f_hat)));
+
+      REQUIRE(t.jacobianDeterminant(m_hat) == Catch::Approx(detJ(m_hat)));
+    }
+
+    SECTION("volume calculation")
+    {
+      // due to the z component of the jacobian determinant, degree 3
+      // polynomials must be exactly integrated
+      PrismGaussQuadrature gauss(3);
+      double volume = 0;
+      for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
+        volume += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i));
+      }
+
+      // 11/2 is actually the volume of the prism
+      REQUIRE(volume == Catch::Approx(11. / 2));
+    }
+
+    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;
+      };
+
+      // 5 is the minimum quadrature rule to integrate the polynomial on the prism
+      PrismGaussQuadrature gauss(5);
+      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(17627. / 720));
+    }
+  }
+}