Add more tests for IntegrateOnCells for cubic cells

tests/test_IntegrateOnCells_cubic.cpp

@@ -14,6 +14,7 @@
 #include <language/utils/SymbolTable.hpp>
 
 #include <MeshDataBaseForTests.hpp>
+#include <geometry/LineCubicTransformation.hpp>
 #include <mesh/Connectivity.hpp>
 #include <mesh/DualMeshManager.hpp>
 #include <mesh/Mesh.hpp>
@@ -59,6 +60,8 @@ TEST_CASE("IntegrateOnCells_cubic", "[language]")
     return std::sqrt(error);
   };
 
+  SECTION("straight edges")
+  {
     SECTION("Gauss quadrature")
     {
       auto quadrature_descriptor = GaussQuadratureDescriptor(20);
@@ -431,3 +434,147 @@ let R2x2_2d: R^2 -> R^2x2, x -> [[2*x[0]*x[0]*x[1]+3, x[0]-2*x[1]], [3, x[0]*x[1
       }
     }
   }
+
+  SECTION("cubic edges")
+  {
+    constexpr size_t Dimension = 2;
+
+    using R2 = TinyVector<Dimension>;
+
+    auto mesh_transform = [](const R2& x) -> R2 {
+      return R2{x[0] * (1 + 0.2 * std::sin(x[0] + x[1])), x[1] * (1 + 0.3 * std::sin(x[0] + x[1]))};
+    };
+
+    std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+    for (const auto& named_mesh : mesh_list) {
+      SECTION(named_mesh.name())
+      {
+        auto mesh_2d_v = named_mesh.mesh();
+        auto mesh_2d   = mesh_2d_v->get<Mesh<2>>();
+
+        PolynomialMeshBuilder pb{mesh_2d_v, 3};
+
+        auto cubic_mesh = pb.mesh()->get<PolynomialMesh<2>>();
+
+        auto xr = copy(cubic_mesh->xr());
+        auto xl = copy(cubic_mesh->xl());
+
+        parallel_for(
+          cubic_mesh->numberOfNodes(),
+          PUGS_LAMBDA(const NodeId& node_id) { xr[node_id] = mesh_transform(xr[node_id]); });
+
+        parallel_for(
+          cubic_mesh->numberOfFaces(), PUGS_LAMBDA(const FaceId& face_id) {
+            for (size_t i = 0; i < xl.sizeOfArrays(); ++i) {
+              xl[face_id][i] = mesh_transform(xl[face_id][i]);
+            }
+          });
+
+        cubic_mesh = std::make_shared<PolynomialMesh<2>>(cubic_mesh->shared_connectivity(), xr, xl);
+
+        auto Q = [](const R2& X) -> R2 {
+          const double& x = X[0];
+          const double& y = X[1];
+          return R2{3 + 2 * x - 3 * y + 5 * x * y,   //
+                    -4 - 3 * x + y * y + 2 * x * x * y};
+        };
+
+        FaceValue<double> face_integral(cubic_mesh->connectivity());
+        const auto face_to_node_matrix = cubic_mesh->connectivity().faceToNodeMatrix();
+
+        QuadratureFormula<1> qf = QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(11));
+
+        parallel_for(
+          cubic_mesh->numberOfFaces(), PUGS_LAMBDA(const FaceId& face_id) {
+            const NodeId node0 = face_to_node_matrix[face_id][0];
+            const NodeId node1 = face_to_node_matrix[face_id][1];
+
+            LineCubicTransformation<Dimension> T{xr[node0], xl[face_id][0], xl[face_id][1], xr[node1]};
+
+            double sum = 0;
+            for (size_t i = 0; i < qf.numberOfPoints(); ++i) {
+              const auto& x_hat = qf.point(i);
+
+              const R2 Q_i     = Q(T(x_hat));
+              const R2 T_prime = T.velocity(x_hat);
+
+              const R2 N_hat{T_prime[1], -T_prime[0]};
+
+              const double& omega_i = qf.weight(i);
+
+              sum += omega_i * dot(Q_i, N_hat);
+            }
+
+            face_integral[face_id] = sum;
+          });
+
+        auto cell_face_is_reversed     = cubic_mesh->connectivity().cellFaceIsReversed();
+        const auto cell_to_face_matrix = cubic_mesh->connectivity().cellToFaceMatrix();
+
+        CellValue<double> integrate_ref{cubic_mesh->connectivity()};
+
+        parallel_for(
+          cubic_mesh->numberOfCells(), PUGS_LAMBDA(const CellId& cell_id) {
+            const auto cell_faces       = cell_to_face_matrix[cell_id];
+            const auto face_is_reversed = cell_face_is_reversed[cell_id];
+
+            double sum = 0;
+            for (size_t i_face = 0; i_face < cell_faces.size(); ++i_face) {
+              const FaceId face_id = cell_faces[i_face];
+              const double sign    = (face_is_reversed[i_face]) ? -1 : 1;
+
+              sum += sign * face_integral[face_id];
+            }
+
+            integrate_ref[cell_id] = sum;
+          });
+
+        {
+          auto quadrature_descriptor = GaussQuadratureDescriptor(18);
+
+          std::string_view data = R"(
+import math;
+let scalar_2d: R^2 -> R, x -> 2 + 7 * x[1] + 2 * x[0] * x[0];
+)";
+
+          TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+          auto ast = ASTBuilder::build(input);
+
+          ASTModulesImporter{*ast};
+          ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+          ASTSymbolTableBuilder{*ast};
+          ASTNodeDataTypeBuilder{*ast};
+
+          ASTNodeTypeCleaner<language::var_declaration>{*ast};
+          ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+          ASTNodeExpressionBuilder{*ast};
+
+          std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+          // ensure that variables are declared at this point
+          TAO_PEGTL_NAMESPACE::position position{data.size(), 1, 1, "fixture"};
+
+          auto [i_symbol, found] = symbol_table->find("scalar_2d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          Array<double> integrate_value(cubic_mesh->numberOfCells());
+
+          ASTNode::setStackDetails(false);
+
+          IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                       *cubic_mesh, integrate_value);
+
+          ASTNode::setStackDetails(true);
+
+          REQUIRE(scalar_error(integrate_ref.arrayView(), integrate_value) < 1E-12);
+        }
+      }
+    }
+  }
+}