From fe1a504a04dd625622ee61d570b6653196347fe7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Wed, 9 Feb 2022 14:16:18 +0100
Subject: [PATCH] Add tests for IntegrateCellValue

---
 src/language/utils/IntegrateCellValue.hpp |   5 +-
 tests/CMakeLists.txt                      |   1 +
 tests/test_IntegrateCellValue.cpp         | 451 ++++++++++++++++++++++
 tests/test_IntegrateCellValue.hpp         |   1 +
 tests/test_InterpolateItemArray.cpp       |   8 +-
 5 files changed, 459 insertions(+), 7 deletions(-)
 create mode 100644 tests/test_IntegrateCellValue.cpp
 create mode 100644 tests/test_IntegrateCellValue.hpp

diff --git a/src/language/utils/IntegrateCellValue.hpp b/src/language/utils/IntegrateCellValue.hpp
index c4eafc151..cd8bacc07 100644
--- a/src/language/utils/IntegrateCellValue.hpp
+++ b/src/language/utils/IntegrateCellValue.hpp
@@ -34,7 +34,7 @@ class IntegrateCellValue<OutputType(InputType)>
             const Array<const CellId>& list_of_cells)
   {
     return IntegrateOnCells<OutputType(const InputType)>::integrate(function_symbol_id, quadrature_descriptor, mesh,
-                                                                    list_of_cells);
+                                                                    Array<const CellId>{list_of_cells});
   }
 
   template <typename MeshType>
@@ -44,8 +44,7 @@ class IntegrateCellValue<OutputType(InputType)>
             const MeshType& mesh,
             const Array<CellId>& list_of_cells)
   {
-    return IntegrateOnCells<OutputType(const InputType)>::integrate(function_symbol_id, quadrature_descriptor, mesh,
-                                                                    Array<const CellId>{list_of_cells});
+    return integrate(function_symbol_id, quadrature_descriptor, mesh, Array<const CellId>{list_of_cells});
   }
 };
 
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index cf294f9ec..2e8f36589 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -98,6 +98,7 @@ add_executable (unit_tests
   test_GaussQuadratureDescriptor.cpp
   test_IfProcessor.cpp
   test_IncDecExpressionProcessor.cpp
+  test_IntegrateCellValue.cpp
   test_IntegrateOnCells.cpp
   test_INodeProcessor.cpp
   test_ItemId.cpp
diff --git a/tests/test_IntegrateCellValue.cpp b/tests/test_IntegrateCellValue.cpp
new file mode 100644
index 000000000..df35a215b
--- /dev/null
+++ b/tests/test_IntegrateCellValue.cpp
@@ -0,0 +1,451 @@
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <language/ast/ASTBuilder.hpp>
+#include <language/ast/ASTModulesImporter.hpp>
+#include <language/ast/ASTNodeDataTypeBuilder.hpp>
+#include <language/ast/ASTNodeExpressionBuilder.hpp>
+#include <language/ast/ASTNodeFunctionEvaluationExpressionBuilder.hpp>
+#include <language/ast/ASTNodeFunctionExpressionBuilder.hpp>
+#include <language/ast/ASTNodeTypeCleaner.hpp>
+#include <language/ast/ASTSymbolTableBuilder.hpp>
+#include <language/utils/PugsFunctionAdapter.hpp>
+#include <language/utils/SymbolTable.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/DualMeshManager.hpp>
+#include <mesh/Mesh.hpp>
+#include <scheme/CellIntegrator.hpp>
+
+#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
+#include <analysis/GaussLobattoQuadratureDescriptor.hpp>
+#include <analysis/GaussQuadratureDescriptor.hpp>
+
+#include <language/utils/IntegrateCellValue.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("IntegrateCellValue", "[language]")
+{
+  SECTION("integrate on all cells")
+  {
+    auto same_item_integral = [](auto f, auto g) -> bool {
+      using ItemIdType = typename decltype(f)::index_type;
+      for (ItemIdType item_id = 0; item_id < f.numberOfItems(); ++item_id) {
+        if (f[item_id] != g[item_id]) {
+          return false;
+        }
+      }
+
+      return true;
+    };
+
+    SECTION("1D")
+    {
+      constexpr size_t Dimension = 1;
+      auto quadrature_descriptor = GaussQuadratureDescriptor(3);
+
+      std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_1d = named_mesh.mesh();
+
+          std::string_view data = R"(
+import math;
+let R2x2_1d: R^1 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[0]) + 3, sin(x[0] - 2 * x[0]), 3, 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;
+
+          TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+          position.byte = data.size();   // ensure that variables are declared at this point
+
+          using R2x2             = TinyMatrix<2>;
+          auto [i_symbol, found] = symbol_table->find("R2x2_1d", 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);
+
+          CellValue<R2x2> cell_integral{mesh_1d->connectivity()};
+          auto f = [](const TinyVector<Dimension>& x) -> R2x2 {
+            return R2x2{2 * exp(x[0]) * sin(x[0]) + 3, sin(x[0] - 2 * x[0]), 3, x[0] * x[0]};
+          };
+          CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_1d, cell_integral);
+
+          CellValue<R2x2> integrate_value =
+            IntegrateCellValue<R2x2(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                       *mesh_1d);
+
+          REQUIRE(same_item_integral(cell_integral, integrate_value));
+        }
+      }
+    }
+
+    SECTION("2D")
+    {
+      constexpr size_t Dimension = 2;
+      auto quadrature_descriptor = GaussLobattoQuadratureDescriptor(3);
+
+      std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_2d = named_mesh.mesh();
+
+          std::string_view data = R"(
+import math;
+let R3_2d: R^2 -> R^3, x -> (2*exp(x[0])*sin(x[1])+3, x[0]-2*x[1], 3);
+)";
+
+          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;
+
+          TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+          position.byte = data.size();   // ensure that variables are declared at this point
+
+          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);
+
+          CellValue<R3> cell_integral{mesh_2d->connectivity()};
+          auto f = [](const TinyVector<Dimension>& x) -> R3 {
+            return R3{2 * exp(x[0]) * sin(x[1]) + 3, x[0] - 2 * x[1], 3};
+          };
+          CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_2d, cell_integral);
+
+          CellValue<R3> integrate_value =
+            IntegrateCellValue<R3(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                     *mesh_2d);
+
+          REQUIRE(same_item_integral(cell_integral, integrate_value));
+        }
+      }
+    }
+
+    SECTION("3D")
+    {
+      constexpr size_t Dimension = 3;
+      auto quadrature_descriptor = GaussLegendreQuadratureDescriptor(3);
+
+      using NamedMesh = MeshDataBaseForTests::NamedMesh<Dimension>;
+
+      std::vector<NamedMesh> mesh_list = [] {
+        std::vector<NamedMesh> extended_mesh_list;
+        std::array mesh_array = MeshDataBaseForTests::get().all3DMeshes();
+        for (size_t i = 0; i < mesh_array.size(); ++i) {
+          extended_mesh_list.push_back(MeshDataBaseForTests::get().all3DMeshes()[i]);
+        }
+        extended_mesh_list.push_back(NamedMesh("diamond dual", DualMeshManager::instance().getDiamondDualMesh(
+                                                                 *MeshDataBaseForTests::get().hybrid3DMesh())));
+        return extended_mesh_list;
+      }();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_3d = named_mesh.mesh();
+
+          std::string_view data = R"(
+import math;
+let scalar_3d: R^3 -> R, x -> 2 * exp(x[0]) * sin(x[1]) * x[2] + 3;
+)";
+
+          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;
+
+          TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+          position.byte = data.size();   // ensure that variables are declared at this point
+
+          auto [i_symbol, found] = symbol_table->find("scalar_3d", 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);
+
+          CellValue<double> cell_integral{mesh_3d->connectivity()};
+          auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * exp(x[0]) * sin(x[1]) * x[2] + 3; };
+          CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_3d, cell_integral);
+
+          CellValue<double> integrate_value =
+            IntegrateCellValue<double(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                         *mesh_3d);
+
+          REQUIRE(same_item_integral(cell_integral, integrate_value));
+        }
+      }
+    }
+  }
+
+  SECTION("integrate on cell list")
+  {
+    auto same_item_integral = [](auto f, auto g) -> bool {
+      using ItemIdType = typename decltype(g)::index_type;
+      for (ItemIdType item_id = 0; item_id < f.size(); ++item_id) {
+        if (f[item_id] != g[item_id]) {
+          return false;
+        }
+      }
+
+      return true;
+    };
+
+    SECTION("1D")
+    {
+      constexpr size_t Dimension = 1;
+      auto quadrature_descriptor = GaussLegendreQuadratureDescriptor(3);
+
+      std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_1d = named_mesh.mesh();
+          Array<CellId> cell_list{mesh_1d->numberOfCells() / 2 + mesh_1d->numberOfCells() % 2};
+
+          {
+            size_t k = 0;
+            for (CellId cell_id = 0; cell_id < mesh_1d->numberOfCells(); ++(++cell_id), ++k) {
+              cell_list[k] = cell_id;
+            }
+
+            REQUIRE(k == cell_list.size());
+          }
+
+          std::string_view data = R"(
+import math;
+let scalar_1d: R^1 -> R, x -> 2 * exp(x[0]) + 3;
+)";
+          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;
+
+          TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+          position.byte = data.size();   // ensure that variables are declared at this point
+
+          auto [i_symbol, found] = symbol_table->find("scalar_1d", 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);
+
+          auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * std::exp(x[0]) + 3; };
+
+          Array<const double> cell_integral = CellIntegrator::integrate(f, quadrature_descriptor, *mesh_1d, cell_list);
+          Array<const double> integrate_value =
+            IntegrateCellValue<double(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                         *mesh_1d, cell_list);
+
+          REQUIRE(same_item_integral(cell_integral, integrate_value));
+        }
+      }
+    }
+
+    SECTION("2D")
+    {
+      constexpr size_t Dimension = 2;
+      auto quadrature_descriptor = GaussLegendreQuadratureDescriptor(3);
+
+      std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_2d = named_mesh.mesh();
+
+          Array<CellId> cell_list{mesh_2d->numberOfCells() / 2 + mesh_2d->numberOfCells() % 2};
+
+          {
+            size_t k = 0;
+            for (CellId cell_id = 0; cell_id < mesh_2d->numberOfCells(); ++(++cell_id), ++k) {
+              cell_list[k] = cell_id;
+            }
+
+            REQUIRE(k == cell_list.size());
+          }
+
+          std::string_view data = R"(
+import math;
+let R3_2d: R^2 -> R^3, x -> (2*exp(x[0])*sin(x[1])+3, x[0]-2*x[1], 3);
+)";
+
+          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;
+
+          TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+          position.byte = data.size();   // ensure that variables are declared at this point
+
+          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);
+
+          auto f = [](const TinyVector<Dimension>& x) -> R3 {
+            return R3{2 * exp(x[0]) * sin(x[1]) + 3, x[0] - 2 * x[1], 3};
+          };
+
+          Array<const R3> cell_integral = CellIntegrator::integrate(f, quadrature_descriptor, *mesh_2d, cell_list);
+          Array<const R3> integrate_value =
+            IntegrateCellValue<R3(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                     *mesh_2d, cell_list);
+
+          REQUIRE(same_item_integral(cell_integral, integrate_value));
+        }
+      }
+    }
+
+    SECTION("3D")
+    {
+      constexpr size_t Dimension = 3;
+      auto quadrature_descriptor = GaussQuadratureDescriptor(3);
+
+      using NamedMesh = MeshDataBaseForTests::NamedMesh<Dimension>;
+
+      std::vector<NamedMesh> mesh_list = [] {
+        std::vector<NamedMesh> extended_mesh_list;
+        std::array mesh_array = MeshDataBaseForTests::get().all3DMeshes();
+        for (size_t i = 0; i < mesh_array.size(); ++i) {
+          extended_mesh_list.push_back(MeshDataBaseForTests::get().all3DMeshes()[i]);
+        }
+        extended_mesh_list.push_back(NamedMesh("diamond dual", DualMeshManager::instance().getDiamondDualMesh(
+                                                                 *MeshDataBaseForTests::get().hybrid3DMesh())));
+        return extended_mesh_list;
+      }();
+
+      for (auto named_mesh : mesh_list) {
+        SECTION(named_mesh.name())
+        {
+          auto mesh_3d = named_mesh.mesh();
+
+          Array<CellId> cell_list{mesh_3d->numberOfCells() / 2 + mesh_3d->numberOfCells() % 2};
+
+          {
+            size_t k = 0;
+            for (CellId cell_id = 0; cell_id < mesh_3d->numberOfCells(); ++(++cell_id), ++k) {
+              cell_list[k] = cell_id;
+            }
+
+            REQUIRE(k == cell_list.size());
+          }
+
+          std::string_view data = R"(
+import math;
+let R2x2_3d: R^3 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3 * cos(x[2]), sin(x[0] - 2 * x[1] * x[2]), 3, x[0] * x[1] * x[2]);
+)";
+
+          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;
+
+          TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+          position.byte = data.size();   // ensure that variables are declared at this point
+
+          using R2x2             = TinyMatrix<2>;
+          auto [i_symbol, found] = symbol_table->find("R2x2_3d", 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);
+
+          auto f = [](const TinyVector<Dimension>& x) -> R2x2 {
+            return R2x2{2 * exp(x[0]) * sin(x[1]) + 3 * cos(x[2]), sin(x[0] - 2 * x[1] * x[2]), 3, x[0] * x[1] * x[2]};
+          };
+
+          Array<const R2x2> cell_integral = CellIntegrator::integrate(f, quadrature_descriptor, *mesh_3d, cell_list);
+
+          Array<R2x2> integrate_value =
+            IntegrateCellValue<R2x2(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
+                                                                       *mesh_3d, cell_list);
+
+          REQUIRE(same_item_integral(cell_integral, integrate_value));
+        }
+      }
+    }
+  }
+}
diff --git a/tests/test_IntegrateCellValue.hpp b/tests/test_IntegrateCellValue.hpp
new file mode 100644
index 000000000..f9a9e32a9
--- /dev/null
+++ b/tests/test_IntegrateCellValue.hpp
@@ -0,0 +1 @@
+i_f_symboli_f_symboli_f_symbol
diff --git a/tests/test_InterpolateItemArray.cpp b/tests/test_InterpolateItemArray.cpp
index 2974e01d8..7fd968be9 100644
--- a/tests/test_InterpolateItemArray.cpp
+++ b/tests/test_InterpolateItemArray.cpp
@@ -297,10 +297,10 @@ let scalar_non_linear_3d: R^3 -> R, x -> 2 * exp(x[0]) * sin(x[1]) * x[2] + 3;
           }();
 
           std::string_view data = R"(
-  import math;
-  let scalar_affine_1d: R^1 -> R, x -> 2*x[0] + 2;
-  let scalar_non_linear_1d: R^1 -> R, x -> 2 * exp(x[0]) + 3;
-  )";
+import math;
+let scalar_affine_1d: R^1 -> R, x -> 2*x[0] + 2;
+let scalar_non_linear_1d: R^1 -> R, x -> 2 * exp(x[0]) + 3;
+)";
           TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
 
           auto ast = ASTBuilder::build(input);
-- 
GitLab