Add P3 transformation for triangles and related unit tests

src/geometry/TriangleP3Transformation.hpp

+#ifndef TRIANGLE_P3_TRANSFORMATION_HPP
+#define TRIANGLE_P3_TRANSFORMATION_HPP
+
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+
+template <size_t GivenDimension>
+class TriangleP3Transformation;
+
+template <>
+class TriangleP3Transformation<2>
+{
+ public:
+  constexpr static size_t Dimension      = 2;
+  constexpr static size_t NumberOfPoints = 10;
+
+ private:
+  // Points are defined as follows (image of the unit triangle with a right angle in a0)
+  //
+  //   a2
+  //   |  `
+  //  a220  a122
+  //   |     `
+  //  a200  a012    a112
+  //   |         `
+  //   a0- a001- a001 -a1
+
+  TinyMatrix<Dimension, NumberOfPoints> m_points_images;
+
+  PUGS_INLINE
+  double
+  _lambda0(const TinyVector<Dimension>& X_hat) const
+  {
+    const double& x_hat = X_hat[0];
+    const double& y_hat = X_hat[1];
+
+    return 1 - x_hat - y_hat;
+  }
+
+  PUGS_INLINE
+  double
+  _lambda1(const TinyVector<Dimension>& X_hat) const
+  {
+    const double& x_hat = X_hat[0];
+
+    return x_hat;
+  }
+
+  PUGS_INLINE
+  double
+  _lambda2(const TinyVector<Dimension>& X_hat) const
+  {
+    const double& y_hat = X_hat[1];
+
+    return y_hat;
+  }
+
+  PUGS_INLINE
+  double
+  _dx_lambda0(const TinyVector<Dimension>&) const
+  {
+    return -1;
+  }
+
+  PUGS_INLINE
+  double
+  _dy_lambda0(const TinyVector<Dimension>&) const
+  {
+    return -1;
+  }
+
+  PUGS_INLINE
+  double
+  _dx_lambda1(const TinyVector<Dimension>&) const
+  {
+    return 1;
+  }
+
+  PUGS_INLINE
+  double
+  _dy_lambda1(const TinyVector<Dimension>&) const
+  {
+    return 0;
+  }
+
+  PUGS_INLINE
+  double
+  _dx_lambda2(const TinyVector<Dimension>&) const
+  {
+    return 0;
+  }
+
+  PUGS_INLINE
+  double
+  _dy_lambda2(const TinyVector<Dimension>&) const
+  {
+    return 1;
+  }
+
+  // Nodal basis functions have the form lambda*(2*lambda -1)
+  PUGS_INLINE
+  double
+  _wa0(const TinyVector<Dimension>& X_hat) const
+  {
+    const double lambda0 = _lambda0(X_hat);
+
+    return 1 / 2. * lambda0 * (-1 + 3 * lambda0) * (-2 + 3 * lambda0);
+  }
+
+  PUGS_INLINE
+  double
+  _wa1(const TinyVector<Dimension>& X_hat) const
+  {
+    const double lambda1 = _lambda1(X_hat);
+    return 1 / 2. * lambda1 * (-1 + 3 * lambda1) * (-2 + 3 * lambda1);
+  }
+
+  PUGS_INLINE
+  double
+  _wa2(const TinyVector<Dimension>& X_hat) const
+  {
+    const double lambda2 = _lambda2(X_hat);
+
+    return 1 / 2. * lambda2 * (-1 + 3 * lambda2) * (-2 + 3 * lambda2);
+  }
+
+  // Edge basis functions have the form ???
+  PUGS_INLINE
+  double
+  _wa001(const TinyVector<Dimension>& X_hat) const
+  {
+    const double lambda0 = _lambda0(X_hat);
+    const double lambda1 = _lambda1(X_hat);
+
+    return 27. / 2 * lambda0 * lambda1 * (lambda0 - 1. / 3);
+  }
+  PUGS_INLINE
+  double
+  _wa011(const TinyVector<Dimension>& X_hat) const
+  {
+    const double lambda0 = _lambda0(X_hat);
+    const double lambda1 = _lambda1(X_hat);
+
+    return 27. / 2 * lambda0 * lambda1 * (lambda1 - 1. / 3);
+  }
+
+  PUGS_INLINE
+  double
+  _wa112(const TinyVector<Dimension>& X_hat) const
+  {
+    const double lambda1 = _lambda1(X_hat);
+    const double lambda2 = _lambda2(X_hat);
+    return 27. / 2 * lambda1 * lambda2 * (lambda1 - 1. / 3);
+  }
+
+  PUGS_INLINE
+  double
+  _wa122(const TinyVector<Dimension>& X_hat) const
+  {
+    const double lambda1 = _lambda1(X_hat);
+    const double lambda2 = _lambda2(X_hat);
+    return 27. / 2 * lambda1 * lambda2 * (lambda2 - 1. / 3);
+  }
+
+  PUGS_INLINE
+  double
+  _wa220(const TinyVector<Dimension>& X_hat) const
+  {
+    const double lambda0 = _lambda0(X_hat);
+    const double lambda2 = _lambda2(X_hat);
+    return 27. / 2 * lambda0 * lambda2 * (lambda2 - 1. / 3);
+  }
+
+  PUGS_INLINE
+  double
+  _wa200(const TinyVector<Dimension>& X_hat) const
+  {
+    const double lambda0 = _lambda0(X_hat);
+    const double lambda2 = _lambda2(X_hat);
+    return 27. / 2 * lambda0 * lambda2 * (lambda0 - 1. / 3);
+  }
+
+  PUGS_INLINE
+  double
+  _wa012(const TinyVector<Dimension>& X_hat) const
+  {
+    const double lambda0 = _lambda0(X_hat);
+    const double lambda1 = _lambda1(X_hat);
+    const double lambda2 = _lambda2(X_hat);
+
+    return 27. * lambda0 * lambda1 * lambda2;
+  }
+  /////////////////////
+  PUGS_INLINE
+  double
+  _dx_wa0(const TinyVector<Dimension>& X_hat) const
+  {
+    const double lambda0    = _lambda0(X_hat);
+    const double dx_lambda0 = _dx_lambda0(X_hat);
+    return 1. / 2 *
+           ((-1 + 3 * lambda0) * (-2 + 3 * lambda0) + 3 * lambda0 * (-2 + 3 * lambda0) +
+            3 * lambda0 * (-1 + 3 * lambda0)) *
+           dx_lambda0;
+  }
+
+  PUGS_INLINE
+  double
+  _dy_wa0(const TinyVector<Dimension>& X_hat) const
+  {
+    const double lambda0    = _lambda0(X_hat);
+    const double dy_lambda0 = _dy_lambda0(X_hat);
+    return 1. / 2 *
+           ((-1 + 3 * lambda0) * (-2 + 3 * lambda0) + 3 * lambda0 * (-2 + 3 * lambda0) +
+            3 * lambda0 * (-1 + 3 * lambda0)) *
+           dy_lambda0;
+  }
+
+  PUGS_INLINE
+  double
+  _dx_wa1(const TinyVector<Dimension>& X_hat) const
+  {
+    const double lambda1    = _lambda1(X_hat);
+    const double dx_lambda1 = _dx_lambda1(X_hat);
+    return 1. / 2 *
+           ((-1 + 3 * lambda1) * (-2 + 3 * lambda1) + 3 * lambda1 * (-2 + 3 * lambda1) +
+            3 * lambda1 * (-1 + 3 * lambda1)) *
+           dx_lambda1;
+  }
+  PUGS_INLINE
+  double
+  _dy_wa1(const TinyVector<Dimension>& X_hat) const
+  {
+    const double lambda1    = _lambda1(X_hat);
+    const double dy_lambda1 = _dy_lambda1(X_hat);
+    return 1. / 2 *
+           ((-1 + 3 * lambda1) * (-2 + 3 * lambda1) + 3 * lambda1 * (-2 + 3 * lambda1) +
+            3 * lambda1 * (-1 + 3 * lambda1)) *
+           dy_lambda1;
+  }
+  PUGS_INLINE
+  double
+  _dx_wa2(const TinyVector<Dimension>& X_hat) const
+  {
+    const double lambda2    = _lambda2(X_hat);
+    const double dx_lambda2 = _dx_lambda2(X_hat);
+    return 1. / 2 *
+           ((-1 + 3 * lambda2) * (-2 + 3 * lambda2) + 3 * lambda2 * (-2 + 3 * lambda2) +
+            3 * lambda2 * (-1 + 3 * lambda2)) *
+           dx_lambda2;
+  }
+  PUGS_INLINE
+  double
+  _dy_wa2(const TinyVector<Dimension>& X_hat) const
+  {
+    const double lambda2    = _lambda2(X_hat);
+    const double dy_lambda2 = _dy_lambda2(X_hat);
+    return 1. / 2 *
+           ((-1 + 3 * lambda2) * (-2 + 3 * lambda2) + 3 * lambda2 * (-2 + 3 * lambda2) +
+            3 * lambda2 * (-1 + 3 * lambda2)) *
+           dy_lambda2;
+  }
+  ///////////////////
+  PUGS_INLINE
+  double
+  _dx_wa001(const TinyVector<Dimension>& X_hat) const
+  {
+    const double lambda0    = _lambda0(X_hat);
+    const double dx_lambda0 = _dx_lambda0(X_hat);
+    const double lambda1    = _lambda1(X_hat);
+    const double dx_lambda1 = _dx_lambda1(X_hat);
+
+    return 27. / 2 * ((-1. / 3 + lambda0) * lambda0 * dx_lambda1 + lambda1 * (2 * lambda0 - 1. / 3) * dx_lambda0);
+  }
+  PUGS_INLINE
+  double
+  _dy_wa001(const TinyVector<Dimension>& X_hat) const
+  {
+    const double lambda0    = _lambda0(X_hat);
+    const double dy_lambda0 = _dy_lambda0(X_hat);
+    const double lambda1    = _lambda1(X_hat);
+    const double dy_lambda1 = _dy_lambda1(X_hat);
+
+    return 27. / 2 * ((-1. / 3 + lambda0) * lambda0 * dy_lambda1 + lambda1 * (2 * lambda0 - 1. / 3) * dy_lambda0);
+  }
+
+  PUGS_INLINE
+  double
+  _dx_wa011(const TinyVector<Dimension>& X_hat) const
+  {
+    const double lambda0    = _lambda0(X_hat);
+    const double dx_lambda0 = _dx_lambda0(X_hat);
+    const double lambda1    = _lambda1(X_hat);
+    const double dx_lambda1 = _dx_lambda1(X_hat);
+
+    return 27. / 2 * ((-1. / 3 + lambda1) * lambda1 * dx_lambda0 + lambda0 * (2 * lambda1 - 1. / 3) * dx_lambda1);
+  }
+
+  PUGS_INLINE
+  double
+  _dy_wa011(const TinyVector<Dimension>& X_hat) const
+  {
+    const double lambda0    = _lambda0(X_hat);
+    const double dy_lambda0 = _dy_lambda0(X_hat);
+    const double lambda1    = _lambda1(X_hat);
+    const double dy_lambda1 = _dy_lambda1(X_hat);
+
+    return 27. / 2 * ((-1. / 3 + lambda1) * lambda1 * dy_lambda0 + lambda0 * (2 * lambda1 - 1. / 3) * dy_lambda1);
+  }
+
+  PUGS_INLINE
+  double
+  _dx_wa122(const TinyVector<Dimension>& X_hat) const
+  {
+    const double lambda2    = _lambda2(X_hat);
+    const double dx_lambda2 = _dx_lambda2(X_hat);
+    const double lambda1    = _lambda1(X_hat);
+    const double dx_lambda1 = _dx_lambda1(X_hat);
+
+    return 27. / 2 * ((-1. / 3 + lambda2) * lambda2 * dx_lambda1 + lambda1 * (2 * lambda2 - 1. / 3) * dx_lambda2);
+  }
+
+  PUGS_INLINE
+  double
+  _dy_wa122(const TinyVector<Dimension>& X_hat) const
+  {
+    const double lambda2    = _lambda2(X_hat);
+    const double dy_lambda2 = _dy_lambda2(X_hat);
+    const double lambda1    = _lambda1(X_hat);
+    const double dy_lambda1 = _dy_lambda1(X_hat);
+
+    return 27. / 2 * ((-1. / 3 + lambda2) * lambda2 * dy_lambda1 + lambda1 * (2 * lambda2 - 1. / 3) * dy_lambda2);
+  }
+
+  PUGS_INLINE
+  double
+  _dx_wa112(const TinyVector<Dimension>& X_hat) const
+  {
+    const double lambda2    = _lambda2(X_hat);
+    const double dx_lambda2 = _dx_lambda2(X_hat);
+    const double lambda1    = _lambda1(X_hat);
+    const double dx_lambda1 = _dx_lambda1(X_hat);
+
+    return 27. / 2 * ((-1. / 3 + lambda1) * lambda1 * dx_lambda2 + lambda2 * (2 * lambda1 - 1. / 3) * dx_lambda1);
+  }
+
+  PUGS_INLINE
+  double
+  _dy_wa112(const TinyVector<Dimension>& X_hat) const
+  {
+    const double lambda2    = _lambda2(X_hat);
+    const double dy_lambda2 = _dy_lambda2(X_hat);
+    const double lambda1    = _lambda1(X_hat);
+    const double dy_lambda1 = _dy_lambda1(X_hat);
+
+    return 27. / 2 * ((-1. / 3 + lambda1) * lambda1 * dy_lambda2 + lambda2 * (2 * lambda1 - 1. / 3) * dy_lambda1);
+  }
+
+  PUGS_INLINE
+  double
+  _dx_wa220(const TinyVector<Dimension>& X_hat) const
+  {
+    const double lambda2    = _lambda2(X_hat);
+    const double dx_lambda2 = _dx_lambda2(X_hat);
+    const double lambda0    = _lambda0(X_hat);
+    const double dx_lambda0 = _dx_lambda0(X_hat);
+
+    return 27. / 2 * ((-1. / 3 + lambda2) * lambda2 * dx_lambda0 + lambda0 * (2 * lambda2 - 1. / 3) * dx_lambda2);
+  }
+
+  PUGS_INLINE
+  double
+  _dy_wa220(const TinyVector<Dimension>& X_hat) const
+  {
+    const double lambda2    = _lambda2(X_hat);
+    const double dy_lambda2 = _dy_lambda2(X_hat);
+    const double lambda0    = _lambda0(X_hat);
+    const double dy_lambda0 = _dy_lambda0(X_hat);
+
+    return 27. / 2 * ((-1. / 3 + lambda2) * lambda2 * dy_lambda0 + lambda0 * (2 * lambda2 - 1. / 3) * dy_lambda2);
+  }
+
+  PUGS_INLINE
+  double
+  _dx_wa200(const TinyVector<Dimension>& X_hat) const
+  {
+    const double lambda0    = _lambda0(X_hat);
+    const double dx_lambda0 = _dx_lambda0(X_hat);
+    const double lambda2    = _lambda2(X_hat);
+    const double dx_lambda2 = _dx_lambda2(X_hat);
+
+    return 27. / 2 * ((-1. / 3 + lambda0) * lambda0 * dx_lambda2 + lambda2 * (2 * lambda0 - 1. / 3) * dx_lambda0);
+  }
+  PUGS_INLINE
+  double
+  _dy_wa200(const TinyVector<Dimension>& X_hat) const
+  {
+    const double lambda0    = _lambda0(X_hat);
+    const double dy_lambda0 = _dy_lambda0(X_hat);
+    const double lambda2    = _lambda2(X_hat);
+    const double dy_lambda2 = _dy_lambda2(X_hat);
+
+    return 27. / 2 * ((-1. / 3 + lambda0) * lambda0 * dy_lambda2 + lambda2 * (2 * lambda0 - 1. / 3) * dy_lambda0);
+  }
+
+  PUGS_INLINE
+  double
+  _dx_wa012(const TinyVector<Dimension>& X_hat) const
+  {
+    const double lambda0    = _lambda0(X_hat);
+    const double dx_lambda0 = _dx_lambda0(X_hat);
+    const double lambda1    = _lambda1(X_hat);
+    const double dx_lambda1 = _dx_lambda1(X_hat);
+    const double lambda2    = _lambda2(X_hat);
+    const double dx_lambda2 = _dx_lambda2(X_hat);
+
+    return 27 * (dx_lambda0 * lambda1 * lambda2 + lambda0 * dx_lambda1 * lambda2 + lambda0 * lambda1 * dx_lambda2);
+  }
+
+  PUGS_INLINE
+  double
+  _dy_wa012(const TinyVector<Dimension>& X_hat) const
+  {
+    const double lambda0    = _lambda0(X_hat);
+    const double dy_lambda0 = _dy_lambda0(X_hat);
+    const double lambda1    = _lambda1(X_hat);
+    const double dy_lambda1 = _dy_lambda1(X_hat);
+    const double lambda2    = _lambda2(X_hat);
+    const double dy_lambda2 = _dy_lambda2(X_hat);
+
+    return 27 * (dy_lambda0 * lambda1 * lambda2 + lambda0 * dy_lambda1 * lambda2 + lambda0 * lambda1 * dy_lambda2);
+  }
+
+ public:
+  PUGS_INLINE
+  TinyVector<Dimension>
+  operator()(const TinyVector<2>& X_hat) const
+  {
+    TinyVector<NumberOfPoints> Wij{_wa0(X_hat),   _wa1(X_hat),   _wa2(X_hat),     //
+                                   _wa112(X_hat), _wa122(X_hat), _wa220(X_hat),   //
+                                   _wa200(X_hat), _wa001(X_hat), _wa011(X_hat),   //
+                                   _wa012(X_hat)};
+
+    return m_points_images * Wij;
+  }
+
+  double
+  jacobianDeterminant(const TinyVector<2>& X_hat) const
+  {
+    const TinyVector<NumberOfPoints> dx_Wij{_dx_wa0(X_hat),   _dx_wa1(X_hat),   _dx_wa2(X_hat),     //
+                                            _dx_wa112(X_hat), _dx_wa122(X_hat), _dx_wa220(X_hat),   //
+                                            _dx_wa200(X_hat), _dx_wa001(X_hat), _dx_wa011(X_hat),   //
+                                            _dx_wa012(X_hat)};
+
+    const TinyVector<NumberOfPoints> dy_Wij{_dy_wa0(X_hat),   _dy_wa1(X_hat),   _dy_wa2(X_hat),     //
+                                            _dy_wa112(X_hat), _dy_wa122(X_hat), _dy_wa220(X_hat),   //
+                                            _dy_wa200(X_hat), _dy_wa001(X_hat), _dy_wa011(X_hat),   //
+                                            _dy_wa012(X_hat)};
+
+    const TinyMatrix<Dimension, Dimension> J{m_points_images * dx_Wij,   //
+                                             m_points_images * dy_Wij};
+    return det(J);
+  }
+
+  PUGS_INLINE
+  TriangleP3Transformation(const TinyVector<Dimension>& a0,
+                           const TinyVector<Dimension>& a1,
+                           const TinyVector<Dimension>& a2,
+                           const TinyVector<Dimension>& a112,
+                           const TinyVector<Dimension>& a122,
+                           const TinyVector<Dimension>& a220,
+                           const TinyVector<Dimension>& a200,
+                           const TinyVector<Dimension>& a001,
+                           const TinyVector<Dimension>& a011,
+                           const TinyVector<Dimension>& a012)
+
+    : m_points_images{a0, a1, a2, a112, a122, a220, a200, a001, a011, a012}
+  {}
+
+  ~TriangleP3Transformation() = default;
+};
+
+#endif   // TRIANGLE_P3_TRANSFORMATION_HPP

tests/CMakeLists.txt

@@ -156,6 +156,7 @@ add_executable (unit_tests
   test_TinyVector.cpp
   test_TriangleGaussQuadrature.cpp
   test_TriangleP2Transformation.cpp
+  test_TriangleP3Transformation.cpp
   test_TriangleTransformation.cpp
   test_TupleToVectorProcessor.cpp
   test_UnaryExpressionProcessor.cpp

tests/test_TriangleP3Transformation.cpp

+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <analysis/GaussQuadratureDescriptor.hpp>
+#include <analysis/QuadratureManager.hpp>
+#include <geometry/TriangleP3Transformation.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("TriangleP3Transformation", "[geometry]")
+{
+  SECTION("2D")
+  {
+    using R2 = TinyVector<2>;
+
+    const R2 a0{1, 2};
+    const R2 a1{5, 1};
+    const R2 a2{2, 4};
+
+    const R2 a001{4. / 3, 1.6};
+    const R2 a011{8. / 3, 0.8};
+
+    const R2 a112{3.9, 1.9};
+    const R2 a122{3.2, 2.9};
+
+    const R2 a220{1.8, 3.2};
+    const R2 a200{3.2, 2.9};
+
+    const R2 a012{2.5, 2};
+
+    const TriangleP3Transformation<2> t(a0, a1, a2,   //
+                                        a112, a122,   //
+                                        a220, a200,   //
+                                        a001, a011,   //
+                                        a012);
+
+    const R2 a0_hat{0, 0};
+    const R2 a1_hat{1, 0};
+    const R2 a2_hat{0, 1};
+
+    const R2 a001_hat{1. / 3, 0};
+    const R2 a011_hat{2. / 3, 0};
+
+    const R2 a112_hat{2. / 3, 1. / 3};
+    const R2 a122_hat{1. / 3, 2. / 3};
+
+    const R2 a220_hat{0, 2. / 3};
+    const R2 a200_hat{0, 1. / 3};
+
+    const R2 a012_hat{1. / 3, 1. / 3};
+
+    CHECK(t(a0_hat)[0] == Catch::Approx(a0[0]));
+    CHECK(t(a0_hat)[1] == Catch::Approx(a0[1]));
+
+    CHECK(t(a1_hat)[0] == Catch::Approx(a1[0]));
+    CHECK(t(a1_hat)[1] == Catch::Approx(a1[1]));
+
+    CHECK(t(a2_hat)[0] == Catch::Approx(a2[0]));
+    CHECK(t(a2_hat)[1] == Catch::Approx(a2[1]));
+
+    CHECK(t(a001_hat)[0] == Catch::Approx(a001[0]));
+    CHECK(t(a001_hat)[1] == Catch::Approx(a001[1]));
+
+    CHECK(t(a011_hat)[0] == Catch::Approx(a011[0]));
+    CHECK(t(a011_hat)[1] == Catch::Approx(a011[1]));
+
+    CHECK(t(a122_hat)[0] == Catch::Approx(a122[0]));
+    CHECK(t(a122_hat)[1] == Catch::Approx(a122[1]));
+
+    CHECK(t(a112_hat)[0] == Catch::Approx(a112[0]));
+    CHECK(t(a112_hat)[1] == Catch::Approx(a112[1]));
+    //
+    CHECK(t(a220_hat)[0] == Catch::Approx(a220[0]));
+    CHECK(t(a220_hat)[1] == Catch::Approx(a220[1]));
+
+    CHECK(t(a200_hat)[0] == Catch::Approx(a200[0]));
+    CHECK(t(a200_hat)[1] == Catch::Approx(a200[1]));
+
+    CHECK(t(a012_hat)[0] == Catch::Approx(a012[0]));
+    CHECK(t(a012_hat)[1] == Catch::Approx(a012[1]));
+
+    // Values obtained using Maxima
+
+    const R2 X0_hat{0.2, 0.3};
+    REQUIRE(t(X0_hat)[0] == Catch::Approx(2.5464));
+    REQUIRE(t(X0_hat)[1] == Catch::Approx(2.31225));
+
+    const R2 X1_hat{0.6, 0.1};
+    REQUIRE(t(X1_hat)[0] == Catch::Approx(2.7548));
+    REQUIRE(t(X1_hat)[1] == Catch::Approx(1.10455));
+
+    const R2 X2_hat{0.8, 0.05};
+    REQUIRE(t(X2_hat)[0] == Catch::Approx(3.669725));
+    REQUIRE(t(X2_hat)[1] == Catch::Approx(0.8521937500000003));
+
+    auto detJ = [](const R2& X) {
+      const double x  = X[0];
+      const double y  = X[1];
+      const double x2 = x * x;
+      const double y2 = y * y;
+      const double x3 = x * x * x;
+      const double x4 = x * x * x * x;
+      const double y3 = y * y * y;
+      const double y4 = y * y * y * y;
+
+      return -414.315 * y4 - 1833.435 * x * y3 + 1150.4025 * y3 - 2720.9925 * x2 * y2 + 3589.7175 * x * y2 -
+             862.56 * y2 - 1598.94 * x3 * y + 3373.4475 * x2 * y - 1920.96 * x * y + 193.185 * y - 280.665 * x4 +
+             922.185 * x3 - 937.98 * x2 + 351.045 * x - 16.11;
+    };
+
+    REQUIRE(t.jacobianDeterminant(X0_hat) == Catch::Approx(detJ(X0_hat)));
+    REQUIRE(t.jacobianDeterminant(X1_hat) == Catch::Approx(detJ(X1_hat)));
+    REQUIRE(t.jacobianDeterminant(X2_hat) == Catch::Approx(detJ(X2_hat)));
+
+    QuadratureFormula<2> qf = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor(4));
+
+    double surface = 0;
+    for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+      surface += qf.weight(i) * t.jacobianDeterminant(qf.point(i));
+    }
+
+    // 113. / 30 is actually the surface of the P3-quadrangle
+    // (computed with Maxima)
+    REQUIRE(surface == Catch::Approx(42891. / 8000));
+  }
+}