diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 2e8f36589f0f340dc9ca682d64ae1a52df760e33..ea44f119d6096cbd1b304c39c00400eefec16ca2 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_IntegrateCellArray.cpp
test_IntegrateCellValue.cpp
test_IntegrateOnCells.cpp
test_INodeProcessor.cpp
diff --git a/tests/test_IntegrateCellArray.cpp b/tests/test_IntegrateCellArray.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b6b8db8c8ca7bc0aa895f8be6b7d5f8dbccd89ff
--- /dev/null
+++ b/tests/test_IntegrateCellArray.cpp
@@ -0,0 +1,626 @@
+#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/IntegrateCellArray.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("IntegrateCellArray", "[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) {
+ for (size_t i = 0; i < f.sizeOfArrays(); ++i) {
+ if (f[item_id][i] != g[item_id][i]) {
+ 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;
+import math;
+let f: R^1 -> R, x -> 2*x[0] + 2;
+let g: 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
+
+ std::vector<FunctionSymbolId> function_symbol_id_list;
+
+ {
+ auto [i_symbol, found] = symbol_table->find("f", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ function_symbol_id_list.push_back(
+ FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+ }
+
+ {
+ auto [i_symbol, found] = symbol_table->find("g", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ function_symbol_id_list.push_back(
+ FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+ }
+
+ CellArray<double> cell_integral_array{mesh_1d->connectivity(), 2};
+
+ {
+ CellValue<double> cell_f_integral{mesh_1d->connectivity()};
+ auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * x[0] + 2; };
+ CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_1d, cell_f_integral);
+
+ parallel_for(
+ mesh_1d->numberOfCells(),
+ PUGS_LAMBDA(const CellId cell_id) { cell_integral_array[cell_id][0] = cell_f_integral[cell_id]; });
+
+ CellValue<double> cell_g_integral{mesh_1d->connectivity()};
+ auto g = [](const TinyVector<Dimension>& x) -> double { return 2 * exp(x[0]) + 3; };
+ CellIntegrator::integrateTo(g, quadrature_descriptor, *mesh_1d, cell_g_integral);
+
+ parallel_for(
+ mesh_1d->numberOfCells(),
+ PUGS_LAMBDA(const CellId cell_id) { cell_integral_array[cell_id][1] = cell_g_integral[cell_id]; });
+ }
+
+ CellArray<double> integrate_array =
+ IntegrateCellArray<double(TinyVector<Dimension>)>::integrate(function_symbol_id_list, quadrature_descriptor,
+ *mesh_1d);
+
+ REQUIRE(same_item_integral(cell_integral_array, integrate_array));
+ }
+ }
+ }
+
+ 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 f: R^2 -> R, x -> 2*x[0] + 3*x[1] + 2;
+let g: R^2 -> R, x -> 2*exp(x[0])*sin(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
+
+ std::vector<FunctionSymbolId> function_symbol_id_list;
+
+ {
+ auto [i_symbol, found] = symbol_table->find("f", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ function_symbol_id_list.push_back(
+ FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+ }
+
+ {
+ auto [i_symbol, found] = symbol_table->find("g", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ function_symbol_id_list.push_back(
+ FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+ }
+
+ CellArray<double> cell_integral_array{mesh_2d->connectivity(), 2};
+
+ {
+ CellValue<double> cell_f_integral{mesh_2d->connectivity()};
+ auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * x[0] + 3 * x[1] + 2; };
+ CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_2d, cell_f_integral);
+
+ parallel_for(
+ mesh_2d->numberOfCells(),
+ PUGS_LAMBDA(const CellId cell_id) { cell_integral_array[cell_id][0] = cell_f_integral[cell_id]; });
+
+ CellValue<double> cell_g_integral{mesh_2d->connectivity()};
+ auto g = [](const TinyVector<Dimension>& x) -> double { return 2 * exp(x[0]) * sin(x[1]) + 3; };
+ CellIntegrator::integrateTo(g, quadrature_descriptor, *mesh_2d, cell_g_integral);
+
+ parallel_for(
+ mesh_2d->numberOfCells(),
+ PUGS_LAMBDA(const CellId cell_id) { cell_integral_array[cell_id][1] = cell_g_integral[cell_id]; });
+ }
+
+ CellArray<double> integrate_array =
+ IntegrateCellArray<double(TinyVector<Dimension>)>::integrate(function_symbol_id_list, quadrature_descriptor,
+ *mesh_2d);
+
+ REQUIRE(same_item_integral(cell_integral_array, integrate_array));
+ }
+ }
+ }
+
+ 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 f: R^3 -> R, x -> 2 * x[0] + 3 * x[1] + 2 * x[2] - 1;
+let g: 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
+
+ std::vector<FunctionSymbolId> function_symbol_id_list;
+
+ {
+ auto [i_symbol, found] = symbol_table->find("f", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ function_symbol_id_list.push_back(
+ FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+ }
+
+ {
+ auto [i_symbol, found] = symbol_table->find("g", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ function_symbol_id_list.push_back(
+ FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+ }
+
+ CellArray<double> cell_integral_array{mesh_3d->connectivity(), 2};
+
+ {
+ CellValue<double> cell_f_integral{mesh_3d->connectivity()};
+ auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * x[0] + 3 * x[1] + 2 * x[2] - 1; };
+ CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_3d, cell_f_integral);
+
+ parallel_for(
+ mesh_3d->numberOfCells(),
+ PUGS_LAMBDA(const CellId cell_id) { cell_integral_array[cell_id][0] = cell_f_integral[cell_id]; });
+
+ CellValue<double> cell_g_integral{mesh_3d->connectivity()};
+ auto g = [](const TinyVector<Dimension>& x) -> double { return 2 * exp(x[0]) * sin(x[1]) * x[2] + 3; };
+ CellIntegrator::integrateTo(g, quadrature_descriptor, *mesh_3d, cell_g_integral);
+
+ parallel_for(
+ mesh_3d->numberOfCells(),
+ PUGS_LAMBDA(const CellId cell_id) { cell_integral_array[cell_id][1] = cell_g_integral[cell_id]; });
+ }
+
+ CellArray<double> integrate_array =
+ IntegrateCellArray<double(TinyVector<Dimension>)>::integrate(function_symbol_id_list, quadrature_descriptor,
+ *mesh_3d);
+
+ REQUIRE(same_item_integral(cell_integral_array, integrate_array));
+ }
+ }
+ }
+ }
+
+ SECTION("integrate on cell list")
+ {
+ auto same_item_integral = [](auto f, auto g) -> bool {
+ using ItemIdType = typename decltype(f)::index_type;
+ for (ItemIdType item_id = 0; item_id < f.numberOfRows(); ++item_id) {
+ for (size_t i = 0; i < f.numberOfColumns(); ++i) {
+ if (f[item_id][i] != g[item_id][i]) {
+ 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 f: R^1 -> R, x -> 2*x[0] + 2;
+let g: 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
+
+ std::vector<FunctionSymbolId> function_symbol_id_list;
+
+ {
+ auto [i_symbol, found] = symbol_table->find("f", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ function_symbol_id_list.push_back(
+ FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+ }
+
+ {
+ auto [i_symbol, found] = symbol_table->find("g", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ function_symbol_id_list.push_back(
+ FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+ }
+
+ Table<double> cell_integral_array{cell_list.size(), 2};
+
+ {
+ auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * x[0] + 2; };
+ Array<double> cell_f_integral = CellIntegrator::integrate(f, quadrature_descriptor, *mesh_1d, cell_list);
+
+ parallel_for(
+ cell_integral_array.numberOfRows(),
+ PUGS_LAMBDA(const CellId cell_id) { cell_integral_array[cell_id][0] = cell_f_integral[cell_id]; });
+
+ auto g = [](const TinyVector<Dimension>& x) -> double { return 2 * exp(x[0]) + 3; };
+ Array<double> cell_g_integral = CellIntegrator::integrate(g, quadrature_descriptor, *mesh_1d, cell_list);
+
+ parallel_for(
+ cell_integral_array.numberOfRows(),
+ PUGS_LAMBDA(const CellId cell_id) { cell_integral_array[cell_id][1] = cell_g_integral[cell_id]; });
+ }
+
+ Table<const double> integrate_value =
+ IntegrateCellArray<double(TinyVector<Dimension>)>::integrate(function_symbol_id_list, quadrature_descriptor,
+ *mesh_1d, cell_list);
+
+ REQUIRE(same_item_integral(cell_integral_array, 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 f: R^2 -> R, x -> 2*x[0] + 3*x[1] + 2;
+let g: R^2 -> R, x -> 2*exp(x[0])*sin(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
+
+ std::vector<FunctionSymbolId> function_symbol_id_list;
+
+ {
+ auto [i_symbol, found] = symbol_table->find("f", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ function_symbol_id_list.push_back(
+ FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+ }
+
+ {
+ auto [i_symbol, found] = symbol_table->find("g", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ function_symbol_id_list.push_back(
+ FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+ }
+
+ Table<double> cell_integral_array{cell_list.size(), 2};
+
+ {
+ auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * x[0] + 3 * x[1] + 2; };
+ Array<double> cell_f_integral = CellIntegrator::integrate(f, quadrature_descriptor, *mesh_2d, cell_list);
+
+ parallel_for(
+ cell_integral_array.numberOfRows(),
+ PUGS_LAMBDA(const CellId cell_id) { cell_integral_array[cell_id][0] = cell_f_integral[cell_id]; });
+
+ auto g = [](const TinyVector<Dimension>& x) -> double { return 2 * exp(x[0]) * sin(x[1]) + 3; };
+ Array<double> cell_g_integral = CellIntegrator::integrate(g, quadrature_descriptor, *mesh_2d, cell_list);
+
+ parallel_for(
+ cell_integral_array.numberOfRows(),
+ PUGS_LAMBDA(const CellId cell_id) { cell_integral_array[cell_id][1] = cell_g_integral[cell_id]; });
+ }
+
+ Table<const double> integrate_value =
+ IntegrateCellArray<double(TinyVector<Dimension>)>::integrate(function_symbol_id_list, quadrature_descriptor,
+ *mesh_2d, cell_list);
+
+ REQUIRE(same_item_integral(cell_integral_array, 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 f: R^3 -> R, x -> 2 * x[0] + 3 * x[1] + 2 * x[2] - 1;
+let g: 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
+
+ std::vector<FunctionSymbolId> function_symbol_id_list;
+
+ {
+ auto [i_symbol, found] = symbol_table->find("f", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ function_symbol_id_list.push_back(
+ FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+ }
+
+ {
+ auto [i_symbol, found] = symbol_table->find("g", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ function_symbol_id_list.push_back(
+ FunctionSymbolId(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table));
+ }
+
+ Table<double> cell_integral_array{cell_list.size(), 2};
+
+ {
+ auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * x[0] + 3 * x[1] + 2 * x[2] - 1; };
+ Array<double> cell_f_integral = CellIntegrator::integrate(f, quadrature_descriptor, *mesh_3d, cell_list);
+
+ parallel_for(
+ cell_integral_array.numberOfRows(),
+ PUGS_LAMBDA(const CellId cell_id) { cell_integral_array[cell_id][0] = cell_f_integral[cell_id]; });
+
+ auto g = [](const TinyVector<Dimension>& x) -> double { return 2 * exp(x[0]) * sin(x[1]) * x[2] + 3; };
+ Array<double> cell_g_integral = CellIntegrator::integrate(g, quadrature_descriptor, *mesh_3d, cell_list);
+
+ parallel_for(
+ cell_integral_array.numberOfRows(),
+ PUGS_LAMBDA(const CellId cell_id) { cell_integral_array[cell_id][1] = cell_g_integral[cell_id]; });
+ }
+
+ Table<const double> integrate_value =
+ IntegrateCellArray<double(TinyVector<Dimension>)>::integrate(function_symbol_id_list, quadrature_descriptor,
+ *mesh_3d, cell_list);
+
+ REQUIRE(same_item_integral(cell_integral_array, integrate_value));
+ }
+ }
+ }
+ }
+}