From a64dcc84b449deca42d815ba88e147ae9109d0a1 Mon Sep 17 00:00:00 2001
From: MARMAJOU ISABELLE <id.counilh@wanadoo.fr>
Date: Fri, 7 Mar 2025 16:35:17 +0100
Subject: [PATCH] Add more tests for IntegrateOnCells for cubic cells

---
 tests/test_IntegrateOnCells_cubic.cpp | 605 ++++++++++++++++----------
 1 file changed, 376 insertions(+), 229 deletions(-)

diff --git a/tests/test_IntegrateOnCells_cubic.cpp b/tests/test_IntegrateOnCells_cubic.cpp
index 28b60d607..aa378255a 100644
--- a/tests/test_IntegrateOnCells_cubic.cpp
+++ b/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,373 +60,519 @@ TEST_CASE("IntegrateOnCells_cubic", "[language]")
     return std::sqrt(error);
   };
 
-  SECTION("Gauss quadrature")
+  SECTION("straight edges")
   {
-    auto quadrature_descriptor = GaussQuadratureDescriptor(20);
-
-    SECTION("integrate on all cells")
+    SECTION("Gauss quadrature")
     {
-      SECTION("2D")
+      auto quadrature_descriptor = GaussQuadratureDescriptor(20);
+
+      SECTION("integrate on all cells")
       {
-        constexpr size_t Dimension = 2;
+        SECTION("2D")
+        {
+          constexpr size_t Dimension = 2;
 
-        std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+          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>>();
+          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};
+              PolynomialMeshBuilder pb{mesh_2d_v, 3};
 
-            auto cubic_mesh = pb.mesh()->get<PolynomialMesh<2>>();
+              auto cubic_mesh = pb.mesh()->get<PolynomialMesh<2>>();
 
-            std::string_view data = R"(
+              std::string_view data = R"(
 import math;
 let scalar_2d: R^2 -> R, x -> 2*x[0]*x[1]+3;
 let R3_2d: R^2 -> R^3, x -> [2*x[0]*x[0] + x[1]*x[1] - 1, x[0]-2*x[1], 3];
 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]]];
 )";
 
-            TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+              TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
 
-            auto ast = ASTBuilder::build(input);
+              auto ast = ASTBuilder::build(input);
 
-            ASTModulesImporter{*ast};
-            ASTNodeTypeCleaner<language::import_instruction>{*ast};
+              ASTModulesImporter{*ast};
+              ASTNodeTypeCleaner<language::import_instruction>{*ast};
 
-            ASTSymbolTableBuilder{*ast};
-            ASTNodeDataTypeBuilder{*ast};
+              ASTSymbolTableBuilder{*ast};
+              ASTNodeDataTypeBuilder{*ast};
 
-            ASTNodeTypeCleaner<language::var_declaration>{*ast};
-            ASTNodeTypeCleaner<language::fct_declaration>{*ast};
-            ASTNodeExpressionBuilder{*ast};
+              ASTNodeTypeCleaner<language::var_declaration>{*ast};
+              ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+              ASTNodeExpressionBuilder{*ast};
 
-            std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+              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"};
+              // ensure that variables are declared at this point
+              TAO_PEGTL_NAMESPACE::position position{data.size(), 1, 1, "fixture"};
 
-            SECTION("scalar 2d")
-            {
-              auto [i_symbol, found] = symbol_table->find("scalar_2d", position);
-              REQUIRE(found);
-              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+              SECTION("scalar 2d")
+              {
+                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);
+                FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
 
-              ASTNode::setStackDetails(false);
+                ASTNode::setStackDetails(false);
 
-              Array<double> integrate_value_polygonal(mesh_2d->numberOfCells());
-              IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
-                                                                           *mesh_2d, integrate_value_polygonal);
+                Array<double> integrate_value_polygonal(mesh_2d->numberOfCells());
+                IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                             *mesh_2d, integrate_value_polygonal);
 
-              Array<double> integrate_value_cubic(cubic_mesh->numberOfCells());
-              IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
-                                                                           *cubic_mesh, integrate_value_cubic);
+                Array<double> integrate_value_cubic(cubic_mesh->numberOfCells());
+                IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                             *cubic_mesh, integrate_value_cubic);
 
-              ASTNode::setStackDetails(true);
+                ASTNode::setStackDetails(true);
 
-              REQUIRE(scalar_error(integrate_value_cubic, integrate_value_polygonal) < 1E-12);
-            }
+                REQUIRE(scalar_error(integrate_value_cubic, integrate_value_polygonal) < 1E-12);
+              }
 
-            SECTION("vector 2d")
-            {
-              using R3               = TinyVector<3>;
-              auto [i_symbol, found] = symbol_table->find("R3_2d", position);
-              REQUIRE(found);
-              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+              SECTION("vector 2d")
+              {
+                using R3               = TinyVector<3>;
+                auto [i_symbol, found] = symbol_table->find("R3_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);
+                FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
 
-              ASTNode::setStackDetails(false);
+                ASTNode::setStackDetails(false);
 
-              Array<R3> integrate_value_polygonal(mesh_2d->numberOfCells());
-              IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
-                                                                       *mesh_2d, integrate_value_polygonal);
+                Array<R3> integrate_value_polygonal(mesh_2d->numberOfCells());
+                IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                         *mesh_2d, integrate_value_polygonal);
 
-              Array<R3> integrate_value_cubic(cubic_mesh->numberOfCells());
-              IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
-                                                                       *cubic_mesh, integrate_value_cubic);
+                Array<R3> integrate_value_cubic(cubic_mesh->numberOfCells());
+                IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                         *cubic_mesh, integrate_value_cubic);
 
-              ASTNode::setStackDetails(true);
+                ASTNode::setStackDetails(true);
 
-              REQUIRE(vector_error(integrate_value_cubic, integrate_value_polygonal) < 1E-12);
-            }
+                REQUIRE(vector_error(integrate_value_cubic, integrate_value_polygonal) < 1E-12);
+              }
 
-            SECTION("matrix 2d")
-            {
-              using R2x2             = TinyMatrix<2>;
-              auto [i_symbol, found] = symbol_table->find("R2x2_2d", position);
-              REQUIRE(found);
-              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+              SECTION("matrix 2d")
+              {
+                using R2x2             = TinyMatrix<2>;
+                auto [i_symbol, found] = symbol_table->find("R2x2_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);
+                FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
 
-              ASTNode::setStackDetails(false);
+                ASTNode::setStackDetails(false);
 
-              Array<R2x2> integrate_value_polygonal(mesh_2d->numberOfCells());
-              IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
-                                                                         *mesh_2d, integrate_value_polygonal);
+                Array<R2x2> integrate_value_polygonal(mesh_2d->numberOfCells());
+                IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                           *mesh_2d, integrate_value_polygonal);
 
-              Array<R2x2> integrate_value_cubic(cubic_mesh->numberOfCells());
-              IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
-                                                                         *cubic_mesh, integrate_value_cubic);
+                Array<R2x2> integrate_value_cubic(cubic_mesh->numberOfCells());
+                IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                           *cubic_mesh, integrate_value_cubic);
 
-              ASTNode::setStackDetails(true);
+                ASTNode::setStackDetails(true);
 
-              REQUIRE(matrix_error(integrate_value_cubic, integrate_value_polygonal) < 1E-12);
+                REQUIRE(matrix_error(integrate_value_cubic, integrate_value_polygonal) < 1E-12);
+              }
             }
           }
         }
       }
     }
-  }
 
-  SECTION("Gauss-Legendre quadrature")
-  {
-    auto quadrature_descriptor = GaussLegendreQuadratureDescriptor(20);
-
-    SECTION("integrate on all cells")
+    SECTION("Gauss-Legendre quadrature")
     {
-      SECTION("2D")
+      auto quadrature_descriptor = GaussLegendreQuadratureDescriptor(20);
+
+      SECTION("integrate on all cells")
       {
-        constexpr size_t Dimension = 2;
+        SECTION("2D")
+        {
+          constexpr size_t Dimension = 2;
 
-        std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+          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>>();
+          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};
+              PolynomialMeshBuilder pb{mesh_2d_v, 3};
 
-            auto cubic_mesh = pb.mesh()->get<PolynomialMesh<2>>();
+              auto cubic_mesh = pb.mesh()->get<PolynomialMesh<2>>();
 
-            std::string_view data = R"(
+              std::string_view data = R"(
 import math;
 let scalar_2d: R^2 -> R, x -> 2*x[0]*x[1]+3;
 let R3_2d: R^2 -> R^3, x -> [2*x[0]*x[0] + x[1]*x[1] - 1, x[0]-2*x[1], 3];
 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]]];
 )";
 
-            TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+              TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
 
-            auto ast = ASTBuilder::build(input);
+              auto ast = ASTBuilder::build(input);
 
-            ASTModulesImporter{*ast};
-            ASTNodeTypeCleaner<language::import_instruction>{*ast};
+              ASTModulesImporter{*ast};
+              ASTNodeTypeCleaner<language::import_instruction>{*ast};
 
-            ASTSymbolTableBuilder{*ast};
-            ASTNodeDataTypeBuilder{*ast};
+              ASTSymbolTableBuilder{*ast};
+              ASTNodeDataTypeBuilder{*ast};
 
-            ASTNodeTypeCleaner<language::var_declaration>{*ast};
-            ASTNodeTypeCleaner<language::fct_declaration>{*ast};
-            ASTNodeExpressionBuilder{*ast};
+              ASTNodeTypeCleaner<language::var_declaration>{*ast};
+              ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+              ASTNodeExpressionBuilder{*ast};
 
-            std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+              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"};
+              // ensure that variables are declared at this point
+              TAO_PEGTL_NAMESPACE::position position{data.size(), 1, 1, "fixture"};
 
-            SECTION("scalar 2d")
-            {
-              auto [i_symbol, found] = symbol_table->find("scalar_2d", position);
-              REQUIRE(found);
-              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+              SECTION("scalar 2d")
+              {
+                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);
+                FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
 
-              ASTNode::setStackDetails(false);
+                ASTNode::setStackDetails(false);
 
-              Array<double> integrate_value_polygonal(mesh_2d->numberOfCells());
-              IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
-                                                                           *mesh_2d, integrate_value_polygonal);
+                Array<double> integrate_value_polygonal(mesh_2d->numberOfCells());
+                IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                             *mesh_2d, integrate_value_polygonal);
 
-              Array<double> integrate_value_cubic(cubic_mesh->numberOfCells());
-              IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
-                                                                           *cubic_mesh, integrate_value_cubic);
+                Array<double> integrate_value_cubic(cubic_mesh->numberOfCells());
+                IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                             *cubic_mesh, integrate_value_cubic);
 
-              ASTNode::setStackDetails(true);
+                ASTNode::setStackDetails(true);
 
-              REQUIRE(scalar_error(integrate_value_cubic, integrate_value_polygonal) < 1E-12);
-            }
+                REQUIRE(scalar_error(integrate_value_cubic, integrate_value_polygonal) < 1E-12);
+              }
 
-            SECTION("vector 2d")
-            {
-              using R3               = TinyVector<3>;
-              auto [i_symbol, found] = symbol_table->find("R3_2d", position);
-              REQUIRE(found);
-              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+              SECTION("vector 2d")
+              {
+                using R3               = TinyVector<3>;
+                auto [i_symbol, found] = symbol_table->find("R3_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);
+                FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
 
-              ASTNode::setStackDetails(false);
+                ASTNode::setStackDetails(false);
 
-              Array<R3> integrate_value_polygonal(mesh_2d->numberOfCells());
-              IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
-                                                                       *mesh_2d, integrate_value_polygonal);
+                Array<R3> integrate_value_polygonal(mesh_2d->numberOfCells());
+                IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                         *mesh_2d, integrate_value_polygonal);
 
-              Array<R3> integrate_value_cubic(cubic_mesh->numberOfCells());
-              IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
-                                                                       *cubic_mesh, integrate_value_cubic);
+                Array<R3> integrate_value_cubic(cubic_mesh->numberOfCells());
+                IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                         *cubic_mesh, integrate_value_cubic);
 
-              ASTNode::setStackDetails(true);
+                ASTNode::setStackDetails(true);
 
-              REQUIRE(vector_error(integrate_value_cubic, integrate_value_polygonal) < 1E-12);
-            }
+                REQUIRE(vector_error(integrate_value_cubic, integrate_value_polygonal) < 1E-12);
+              }
 
-            SECTION("matrix 2d")
-            {
-              using R2x2             = TinyMatrix<2>;
-              auto [i_symbol, found] = symbol_table->find("R2x2_2d", position);
-              REQUIRE(found);
-              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+              SECTION("matrix 2d")
+              {
+                using R2x2             = TinyMatrix<2>;
+                auto [i_symbol, found] = symbol_table->find("R2x2_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);
+                FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
 
-              ASTNode::setStackDetails(false);
+                ASTNode::setStackDetails(false);
 
-              Array<R2x2> integrate_value_polygonal(mesh_2d->numberOfCells());
-              IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
-                                                                         *mesh_2d, integrate_value_polygonal);
+                Array<R2x2> integrate_value_polygonal(mesh_2d->numberOfCells());
+                IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                           *mesh_2d, integrate_value_polygonal);
 
-              Array<R2x2> integrate_value_cubic(cubic_mesh->numberOfCells());
-              IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
-                                                                         *cubic_mesh, integrate_value_cubic);
+                Array<R2x2> integrate_value_cubic(cubic_mesh->numberOfCells());
+                IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                           *cubic_mesh, integrate_value_cubic);
 
-              ASTNode::setStackDetails(true);
+                ASTNode::setStackDetails(true);
 
-              REQUIRE(matrix_error(integrate_value_cubic, integrate_value_polygonal) < 1E-12);
+                REQUIRE(matrix_error(integrate_value_cubic, integrate_value_polygonal) < 1E-12);
+              }
             }
           }
         }
       }
     }
-  }
-
-  SECTION("Gauss-Lobatto quadrature")
-  {
-    auto quadrature_descriptor = GaussLobattoQuadratureDescriptor(13);
 
-    SECTION("integrate on all cells")
+    SECTION("Gauss-Lobatto quadrature")
     {
-      SECTION("2D")
+      auto quadrature_descriptor = GaussLobattoQuadratureDescriptor(13);
+
+      SECTION("integrate on all cells")
       {
-        constexpr size_t Dimension = 2;
+        SECTION("2D")
+        {
+          constexpr size_t Dimension = 2;
 
-        std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+          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>>();
+          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};
+              PolynomialMeshBuilder pb{mesh_2d_v, 3};
 
-            auto cubic_mesh = pb.mesh()->get<PolynomialMesh<2>>();
+              auto cubic_mesh = pb.mesh()->get<PolynomialMesh<2>>();
 
-            std::string_view data = R"(
+              std::string_view data = R"(
 import math;
 let scalar_2d: R^2 -> R, x -> 2*x[0]*x[1]+3;
 let R3_2d: R^2 -> R^3, x -> [2*x[0]*x[0] + x[1]*x[1] - 1, x[0]-2*x[1], 3];
 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]]];
 )";
 
-            TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+              TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
 
-            auto ast = ASTBuilder::build(input);
+              auto ast = ASTBuilder::build(input);
 
-            ASTModulesImporter{*ast};
-            ASTNodeTypeCleaner<language::import_instruction>{*ast};
+              ASTModulesImporter{*ast};
+              ASTNodeTypeCleaner<language::import_instruction>{*ast};
 
-            ASTSymbolTableBuilder{*ast};
-            ASTNodeDataTypeBuilder{*ast};
+              ASTSymbolTableBuilder{*ast};
+              ASTNodeDataTypeBuilder{*ast};
 
-            ASTNodeTypeCleaner<language::var_declaration>{*ast};
-            ASTNodeTypeCleaner<language::fct_declaration>{*ast};
-            ASTNodeExpressionBuilder{*ast};
+              ASTNodeTypeCleaner<language::var_declaration>{*ast};
+              ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+              ASTNodeExpressionBuilder{*ast};
 
-            std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+              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"};
+              // ensure that variables are declared at this point
+              TAO_PEGTL_NAMESPACE::position position{data.size(), 1, 1, "fixture"};
 
-            SECTION("scalar 2d")
-            {
-              auto [i_symbol, found] = symbol_table->find("scalar_2d", position);
-              REQUIRE(found);
-              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+              SECTION("scalar 2d")
+              {
+                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);
+
+                ASTNode::setStackDetails(false);
+
+                Array<double> integrate_value_polygonal(mesh_2d->numberOfCells());
+                IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                             *mesh_2d, integrate_value_polygonal);
+
+                Array<double> integrate_value_cubic(cubic_mesh->numberOfCells());
+                IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                             *cubic_mesh, integrate_value_cubic);
+
+                ASTNode::setStackDetails(true);
 
-              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+                REQUIRE(scalar_error(integrate_value_cubic, integrate_value_polygonal) < 1E-12);
+              }
 
-              ASTNode::setStackDetails(false);
+              SECTION("vector 2d")
+              {
+                using R3               = TinyVector<3>;
+                auto [i_symbol, found] = symbol_table->find("R3_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);
+
+                ASTNode::setStackDetails(false);
+
+                Array<R3> integrate_value_polygonal(mesh_2d->numberOfCells());
+                IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                         *mesh_2d, integrate_value_polygonal);
+
+                Array<R3> integrate_value_cubic(cubic_mesh->numberOfCells());
+                IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                                                                         *cubic_mesh, integrate_value_cubic);
 
-              Array<double> integrate_value_polygonal(mesh_2d->numberOfCells());
-              IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                ASTNode::setStackDetails(true);
+
+                REQUIRE(vector_error(integrate_value_cubic, integrate_value_polygonal) < 1E-12);
+              }
+
+              SECTION("matrix 2d")
+              {
+                using R2x2             = TinyMatrix<2>;
+                auto [i_symbol, found] = symbol_table->find("R2x2_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);
+
+                ASTNode::setStackDetails(false);
+
+                Array<R2x2> integrate_value_polygonal(mesh_2d->numberOfCells());
+                IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
                                                                            *mesh_2d, integrate_value_polygonal);
 
-              Array<double> integrate_value_cubic(cubic_mesh->numberOfCells());
-              IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
+                Array<R2x2> integrate_value_cubic(cubic_mesh->numberOfCells());
+                IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
                                                                            *cubic_mesh, integrate_value_cubic);
 
-              ASTNode::setStackDetails(true);
+                ASTNode::setStackDetails(true);
 
-              REQUIRE(scalar_error(integrate_value_cubic, integrate_value_polygonal) < 1E-12);
+                REQUIRE(matrix_error(integrate_value_cubic, integrate_value_polygonal) < 1E-8);
+              }
             }
+          }
+        }
+      }
+    }
+  }
 
-            SECTION("vector 2d")
-            {
-              using R3               = TinyVector<3>;
-              auto [i_symbol, found] = symbol_table->find("R3_2d", position);
-              REQUIRE(found);
-              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+  SECTION("cubic edges")
+  {
+    constexpr size_t Dimension = 2;
 
-              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+    using R2 = TinyVector<Dimension>;
 
-              ASTNode::setStackDetails(false);
+    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]))};
+    };
 
-              Array<R3> integrate_value_polygonal(mesh_2d->numberOfCells());
-              IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
-                                                                       *mesh_2d, integrate_value_polygonal);
+    std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
 
-              Array<R3> integrate_value_cubic(cubic_mesh->numberOfCells());
-              IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
-                                                                       *cubic_mesh, integrate_value_cubic);
+    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>>();
 
-              ASTNode::setStackDetails(true);
+        PolynomialMeshBuilder pb{mesh_2d_v, 3};
 
-              REQUIRE(vector_error(integrate_value_cubic, integrate_value_polygonal) < 1E-12);
+        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]);
             }
+          });
 
-            SECTION("matrix 2d")
-            {
-              using R2x2             = TinyMatrix<2>;
-              auto [i_symbol, found] = symbol_table->find("R2x2_2d", position);
-              REQUIRE(found);
-              REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+        cubic_mesh = std::make_shared<PolynomialMesh<2>>(cubic_mesh->shared_connectivity(), xr, xl);
 
-              FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+        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};
+        };
 
-              ASTNode::setStackDetails(false);
+        FaceValue<double> face_integral(cubic_mesh->connectivity());
+        const auto face_to_node_matrix = cubic_mesh->connectivity().faceToNodeMatrix();
 
-              Array<R2x2> integrate_value_polygonal(mesh_2d->numberOfCells());
-              IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
-                                                                         *mesh_2d, integrate_value_polygonal);
+        QuadratureFormula<1> qf = QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(11));
 
-              Array<R2x2> integrate_value_cubic(cubic_mesh->numberOfCells());
-              IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor,
-                                                                         *cubic_mesh, integrate_value_cubic);
+        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];
 
-              ASTNode::setStackDetails(true);
+            LineCubicTransformation<Dimension> T{xr[node0], xl[face_id][0], xl[face_id][1], xr[node1]};
 
-              REQUIRE(matrix_error(integrate_value_cubic, integrate_value_polygonal) < 1E-8);
+            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);
         }
       }
     }
-- 
GitLab