Allow SquareTransformation in dimension 3

src/geometry/SquareTransformation.hpp

@@ -6,7 +6,11 @@
 
 #include <array>
 
-class SquareTransformation
+template <size_t GivenDimension>
+class SquareTransformation;
+
+template <>
+class SquareTransformation<2>
 {
  public:
   constexpr static size_t Dimension      = 2;
@@ -58,4 +62,62 @@ class SquareTransformation
   ~SquareTransformation() = default;
 };
 
+template <size_t GivenDimension>
+class SquareTransformation
+{
+ public:
+  constexpr static size_t Dimension = GivenDimension;
+  static_assert(Dimension == 3, "Square transformation is defined in dimension 2 or 3");
+
+  constexpr static size_t NumberOfPoints = 4;
+
+ private:
+  TinyMatrix<Dimension, NumberOfPoints - 1> m_coefficients;
+  TinyVector<Dimension> m_shift;
+
+ public:
+  PUGS_INLINE
+  TinyVector<Dimension>
+  operator()(const TinyVector<2>& x) const
+  {
+    const TinyVector<NumberOfPoints - 1> X = {x[0], x[1], x[0] * x[1]};
+    return m_coefficients * X + m_shift;
+  }
+
+  PUGS_INLINE double
+  areaVariationNorm(const TinyVector<2>& X) const
+  {
+    const auto& T   = m_coefficients;
+    const double& x = X[0];
+    const double& y = X[1];
+
+    const TinyVector<Dimension> dxT = {T(0, 0) + T(0, 2) * y,   //
+                                       T(1, 0) + T(1, 2) * y,   //
+                                       T(2, 0) + T(2, 2) * y};
+
+    const TinyVector<Dimension> dyT = {T(0, 1) + T(0, 2) * x,   //
+                                       T(1, 1) + T(1, 2) * x,   //
+                                       T(2, 1) + T(2, 2) * x};
+
+    return l2Norm(crossProduct(dxT, dyT));
+  }
+
+  PUGS_INLINE
+  SquareTransformation(const TinyVector<Dimension>& a,
+                       const TinyVector<Dimension>& b,
+                       const TinyVector<Dimension>& c,
+                       const TinyVector<Dimension>& d)
+  {
+    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]);
+    }
+  }
+
+  ~SquareTransformation() = default;
+};
+
 #endif   // SQUARE_TRANSFORMATION_HPP

src/language/utils/IntegrateOnCells.hpp

@@ -158,7 +158,7 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
       } else if constexpr (Dimension == 2) {
         switch (cell_type[cell_id]) {
         case CellType::Triangle: {
-          const SquareTransformation T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]],
+          const SquareTransformation<2> T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]],
                                           xr[cell_to_node[2]]);
 
           for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
@@ -170,7 +170,7 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
           break;
         }
         case CellType::Quadrangle: {
-          const SquareTransformation T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]],
+          const SquareTransformation<2> T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]],
                                           xr[cell_to_node[3]]);
 
           for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
@@ -411,7 +411,7 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
           break;
         }
         case CellType::Quadrangle: {
-          const SquareTransformation T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]],
+          const SquareTransformation<2> T(xr[cell_to_node[0]], xr[cell_to_node[1]], xr[cell_to_node[2]],
                                           xr[cell_to_node[3]]);
           const auto qf = QuadratureManager::instance().getSquareFormula(quadrature_descriptor);
 

tests/test_SquareTransformation.cpp

@@ -8,6 +8,8 @@
 // clazy:excludeall=non-pod-global-static
 
 TEST_CASE("SquareTransformation", "[geometry]")
+{
+  SECTION("2D")
   {
     using R2 = TinyVector<2>;
 
@@ -25,7 +27,7 @@ TEST_CASE("SquareTransformation", "[geometry]")
 
     const R2 m = 0.25 * (a + b + c + d);
 
-  const SquareTransformation t(a, b, c, d);
+    const SquareTransformation<2> t(a, b, c, d);
 
     SECTION("values")
     {
@@ -102,3 +104,109 @@ TEST_CASE("SquareTransformation", "[geometry]")
       }
     }
   }
+
+  SECTION("3D")
+  {
+    using R2 = TinyVector<2>;
+
+    const R2 a_hat = {-1, -1};
+    const R2 b_hat = {+1, -1};
+    const R2 c_hat = {+1, +1};
+    const R2 d_hat = {-1, +1};
+
+    const R2 m_hat = zero;
+
+    using R3 = TinyVector<3>;
+
+    const R3 a = {0, 0, -1};
+    const R3 b = {8, -2, 3};
+    const R3 c = {12, 7, 2};
+    const R3 d = {3, 7, 1};
+
+    const R3 m = 0.25 * (a + b + c + d);
+
+    const SquareTransformation<3> t(a, b, c, d);
+
+    SECTION("values")
+    {
+      REQUIRE(t(a_hat)[0] == Catch::Approx(a[0]));
+      REQUIRE(t(a_hat)[1] == Catch::Approx(a[1]));
+      REQUIRE(t(a_hat)[2] == Catch::Approx(a[2]));
+
+      REQUIRE(t(b_hat)[0] == Catch::Approx(b[0]));
+      REQUIRE(t(b_hat)[1] == Catch::Approx(b[1]));
+      REQUIRE(t(b_hat)[2] == Catch::Approx(b[2]));
+
+      REQUIRE(t(c_hat)[0] == Catch::Approx(c[0]));
+      REQUIRE(t(c_hat)[1] == Catch::Approx(c[1]));
+      REQUIRE(t(c_hat)[2] == Catch::Approx(c[2]));
+
+      REQUIRE(t(d_hat)[0] == Catch::Approx(d[0]));
+      REQUIRE(t(d_hat)[1] == Catch::Approx(d[1]));
+      REQUIRE(t(d_hat)[2] == Catch::Approx(d[2]));
+
+      REQUIRE(t(m_hat)[0] == Catch::Approx(m[0]));
+      REQUIRE(t(m_hat)[1] == Catch::Approx(m[1]));
+      REQUIRE(t(m_hat)[2] == Catch::Approx(m[2]));
+    }
+
+    SECTION("Area variation norm")
+    {
+      auto area_variation_norm = [&](const R2& X) {
+        const double x = X[0];
+        const double y = X[1];
+
+        const R3 J1 = 0.25 * (-a + b + c - d);
+        const R3 J2 = 0.25 * (-a - b + c + d);
+        const R3 J3 = 0.25 * (a - b + c - d);
+
+        return l2Norm(crossProduct(J1 + y * J3, J2 + x * J3));
+      };
+
+      SECTION("at points")
+      {
+        REQUIRE(t.areaVariationNorm(a_hat) == Catch::Approx(area_variation_norm(a_hat)));
+        REQUIRE(t.areaVariationNorm(b_hat) == Catch::Approx(area_variation_norm(b_hat)));
+        REQUIRE(t.areaVariationNorm(c_hat) == Catch::Approx(area_variation_norm(c_hat)));
+        REQUIRE(t.areaVariationNorm(d_hat) == Catch::Approx(area_variation_norm(d_hat)));
+
+        REQUIRE(t.areaVariationNorm(m_hat) == Catch::Approx(area_variation_norm(m_hat)));
+      }
+
+      auto p = [](const R3& X) {
+        const double x = X[0];
+        const double y = X[1];
+        const double z = X[2];
+        return 2 * x * x + 3 * x - 3 * y * y + y + 2 * z * z - 0.5 * z + 2;
+      };
+
+      SECTION("Gauss order 3")
+      {
+        // Due to the area variation term (square root of a
+        // polynomial), Gauss-like quadrature cannot be exact for
+        // general 3d quadrangle.
+        //
+        // Moreover, even the surface of general 3d quadrangle cannot
+        // be computed exactly.
+        //
+        // So for this test we integrate the following function:
+        // f(x,y,z) := p(x,y,z) * |dA(t^-1(x,y,z))|
+        //
+        // dA denotes the area change on the reference element. This
+        // function can be exactly integrated since computing the
+        // integral on the reference quadrangle, one gets a dA^2,
+        // which is a polynomial.
+        const QuadratureFormula<2>& gauss =
+          QuadratureManager::instance().getSquareFormula(GaussLegendreQuadratureDescriptor(4));
+
+        double integral = 0;
+        for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
+          const double fX = (t.areaVariationNorm(gauss.point(i)) * p(t(gauss.point(i))));
+          integral += gauss.weight(i) * t.areaVariationNorm(gauss.point(i)) * fX;
+        }
+
+        REQUIRE(integral == Catch::Approx(60519715. / 576));
+      }
+    }
+  }
+}