diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 4e39374a30ca179ae879bf522b604c5ddc96598b..52a41bcb69fc8181fce73b92059b49cf1cccf314 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -164,6 +164,7 @@ add_executable (mpi_unit_tests
test_DiscreteFunctionP0.cpp
test_DiscreteFunctionP0Vector.cpp
test_DiscreteFunctionVectorIntegrator.cpp
+ test_DiscreteFunctionVectorIntegratorByZone.cpp
test_DiscreteFunctionVectorInterpoler.cpp
test_DiscreteFunctionVectorInterpolerByZone.cpp
test_EmbeddedIDiscreteFunctionMathFunctions.cpp
diff --git a/tests/test_DiscreteFunctionVectorIntegratorByZone.cpp b/tests/test_DiscreteFunctionVectorIntegratorByZone.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..98111d4c6b7f1dd4788745db480e72c06c083945
--- /dev/null
+++ b/tests/test_DiscreteFunctionVectorIntegratorByZone.cpp
@@ -0,0 +1,440 @@
+#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/IntegrateCellValue.hpp>
+#include <language/utils/PugsFunctionAdapter.hpp>
+#include <language/utils/SymbolTable.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+#include <analysis/GaussQuadratureDescriptor.hpp>
+#include <analysis/IQuadratureDescriptor.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshCellZone.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <mesh/NamedZoneDescriptor.hpp>
+
+#include <scheme/DiscreteFunctionDescriptorP0Vector.hpp>
+#include <scheme/DiscreteFunctionP0Vector.hpp>
+#include <scheme/DiscreteFunctionVectorIntegrator.hpp>
+
+#include <pegtl/string_input.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("DiscreteFunctionVectorIntegratorByZone", "[scheme]")
+{
+ auto same_cell_value = [](const CellValue<const double>& fi, const size_t i, const auto& f) -> bool {
+ for (CellId cell_id = 0; cell_id < fi.numberOfItems(); ++cell_id) {
+ if (fi[cell_id] != f[cell_id][i]) {
+ return false;
+ }
+ }
+
+ return true;
+ };
+
+ auto register_function = [](const TAO_PEGTL_NAMESPACE::position& position,
+ const std::shared_ptr<SymbolTable>& symbol_table, const std::string& name,
+ std::vector<FunctionSymbolId>& function_id_list) {
+ auto [i_symbol, found] = symbol_table->find(name, 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);
+ function_id_list.push_back(function_symbol_id);
+ };
+
+ SECTION("1D")
+ {
+ constexpr size_t Dimension = 1;
+
+ std::shared_ptr<const IQuadratureDescriptor> quadrature_descriptor = std::make_shared<GaussQuadratureDescriptor>(3);
+
+ auto mesh_1d = MeshDataBaseForTests::get().unordered1DMesh();
+
+ std::vector<std::shared_ptr<const IZoneDescriptor>> zone_list;
+ zone_list.push_back(std::make_shared<NamedZoneDescriptor>("LEFT"));
+
+ auto mesh_cell_zone = getMeshCellZone(*mesh_1d, *zone_list[0]);
+ auto zone_cell_list = mesh_cell_zone.cellList();
+
+ std::string_view data = R"(
+import math;
+let B_scalar_non_linear_1d: R^1 -> B, x -> (exp(2 * x[0]) + 3 < 3.3);
+let N_scalar_non_linear_1d: R^1 -> N, x -> floor(3 * x[0] * x[0] + 2);
+let Z_scalar_non_linear_1d: R^1 -> Z, x -> floor(exp(2 * x[0]) - 1);
+let R_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);
+
+ ASTModulesImporter{*ast};
+ ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+ ASTSymbolTableBuilder{*ast};
+ ASTNodeDataTypeBuilder{*ast};
+
+ ASTNodeTypeCleaner<language::var_declaration>{*ast};
+ ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+ ASTNodeExpressionBuilder{*ast};
+
+ 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::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+ std::vector<FunctionSymbolId> function_id_list;
+ register_function(position, symbol_table, "B_scalar_non_linear_1d", function_id_list);
+ register_function(position, symbol_table, "N_scalar_non_linear_1d", function_id_list);
+ register_function(position, symbol_table, "Z_scalar_non_linear_1d", function_id_list);
+ register_function(position, symbol_table, "R_scalar_non_linear_1d", function_id_list);
+
+ DiscreteFunctionVectorIntegrator integrator(mesh_1d, zone_list, quadrature_descriptor,
+ std::make_shared<DiscreteFunctionDescriptorP0Vector>(),
+ function_id_list);
+ std::shared_ptr discrete_function = integrator.integrate();
+
+ size_t i = 0;
+
+ {
+ Array<double> array =
+ IntegrateCellValue<double(TinyVector<1>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_1d,
+ zone_cell_list);
+
+ CellValue<double> cell_value{mesh_1d->connectivity()};
+ cell_value.fill(0);
+
+ parallel_for(
+ zone_cell_list.size(), PUGS_LAMBDA(const size_t j) {
+ const CellId cell_id = zone_cell_list[j];
+ cell_value[cell_id] = array[j];
+ });
+
+ REQUIRE(same_cell_value(cell_value, i++,
+ dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+ }
+
+ {
+ Array<double> array =
+ IntegrateCellValue<double(TinyVector<1>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_1d,
+ zone_cell_list);
+
+ CellValue<double> cell_value{mesh_1d->connectivity()};
+ cell_value.fill(0);
+
+ parallel_for(
+ zone_cell_list.size(), PUGS_LAMBDA(const size_t j) {
+ const CellId cell_id = zone_cell_list[j];
+ cell_value[cell_id] = array[j];
+ });
+
+ REQUIRE(same_cell_value(cell_value, i++,
+ dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+ }
+
+ {
+ Array<double> array =
+ IntegrateCellValue<double(TinyVector<1>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_1d,
+ zone_cell_list);
+
+ CellValue<double> cell_value{mesh_1d->connectivity()};
+ cell_value.fill(0);
+
+ parallel_for(
+ zone_cell_list.size(), PUGS_LAMBDA(const size_t j) {
+ const CellId cell_id = zone_cell_list[j];
+ cell_value[cell_id] = array[j];
+ });
+
+ REQUIRE(same_cell_value(cell_value, i++,
+ dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+ }
+
+ {
+ Array<double> array =
+ IntegrateCellValue<double(TinyVector<1>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_1d,
+ zone_cell_list);
+
+ CellValue<double> cell_value{mesh_1d->connectivity()};
+ cell_value.fill(0);
+
+ parallel_for(
+ zone_cell_list.size(), PUGS_LAMBDA(const size_t j) {
+ const CellId cell_id = zone_cell_list[j];
+ cell_value[cell_id] = array[j];
+ });
+
+ REQUIRE(same_cell_value(cell_value, i++,
+ dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+ }
+
+ REQUIRE(i == function_id_list.size());
+ }
+
+ SECTION("2D")
+ {
+ constexpr size_t Dimension = 2;
+
+ std::shared_ptr<const IQuadratureDescriptor> quadrature_descriptor = std::make_shared<GaussQuadratureDescriptor>(3);
+
+ auto mesh_2d = MeshDataBaseForTests::get().hybrid2DMesh();
+
+ std::vector<std::shared_ptr<const IZoneDescriptor>> zone_list;
+ zone_list.push_back(std::make_shared<NamedZoneDescriptor>("LEFT"));
+
+ auto mesh_cell_zone = getMeshCellZone(*mesh_2d, *zone_list[0]);
+ auto zone_cell_list = mesh_cell_zone.cellList();
+
+ std::string_view data = R"(
+import math;
+let B_scalar_non_linear_2d: R^2 -> B, x -> (exp(2 * x[0]) + 3 > 3.3);
+let N_scalar_non_linear_2d: R^2 -> N, x -> floor(3 * (x[0] * x[1]) * (x[0] * x[1]) + 2);
+let Z_scalar_non_linear_2d: R^2 -> Z, x -> floor(exp(2 * x[1]) - 1);
+let R_scalar_non_linear_2d: R^2 -> R, x -> 2 * exp(x[0] + 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};
+
+ 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::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+ std::vector<FunctionSymbolId> function_id_list;
+ register_function(position, symbol_table, "B_scalar_non_linear_2d", function_id_list);
+ register_function(position, symbol_table, "N_scalar_non_linear_2d", function_id_list);
+ register_function(position, symbol_table, "Z_scalar_non_linear_2d", function_id_list);
+ register_function(position, symbol_table, "R_scalar_non_linear_2d", function_id_list);
+
+ DiscreteFunctionVectorIntegrator integrator(mesh_2d, zone_list, quadrature_descriptor,
+ std::make_shared<DiscreteFunctionDescriptorP0Vector>(),
+ function_id_list);
+ std::shared_ptr discrete_function = integrator.integrate();
+
+ size_t i = 0;
+
+ {
+ Array<double> array =
+ IntegrateCellValue<double(TinyVector<2>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_2d,
+ zone_cell_list);
+
+ CellValue<double> cell_value{mesh_2d->connectivity()};
+ cell_value.fill(0);
+
+ parallel_for(
+ zone_cell_list.size(), PUGS_LAMBDA(const size_t j) {
+ const CellId cell_id = zone_cell_list[j];
+ cell_value[cell_id] = array[j];
+ });
+
+ REQUIRE(same_cell_value(cell_value, i++,
+ dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+ }
+
+ {
+ Array<double> array =
+ IntegrateCellValue<double(TinyVector<2>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_2d,
+ zone_cell_list);
+
+ CellValue<double> cell_value{mesh_2d->connectivity()};
+ cell_value.fill(0);
+
+ parallel_for(
+ zone_cell_list.size(), PUGS_LAMBDA(const size_t j) {
+ const CellId cell_id = zone_cell_list[j];
+ cell_value[cell_id] = array[j];
+ });
+
+ REQUIRE(same_cell_value(cell_value, i++,
+ dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+ }
+
+ {
+ Array<double> array =
+ IntegrateCellValue<double(TinyVector<2>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_2d,
+ zone_cell_list);
+
+ CellValue<double> cell_value{mesh_2d->connectivity()};
+ cell_value.fill(0);
+
+ parallel_for(
+ zone_cell_list.size(), PUGS_LAMBDA(const size_t j) {
+ const CellId cell_id = zone_cell_list[j];
+ cell_value[cell_id] = array[j];
+ });
+
+ REQUIRE(same_cell_value(cell_value, i++,
+ dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+ }
+
+ {
+ Array<double> array =
+ IntegrateCellValue<double(TinyVector<2>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_2d,
+ zone_cell_list);
+
+ CellValue<double> cell_value{mesh_2d->connectivity()};
+ cell_value.fill(0);
+
+ parallel_for(
+ zone_cell_list.size(), PUGS_LAMBDA(const size_t j) {
+ const CellId cell_id = zone_cell_list[j];
+ cell_value[cell_id] = array[j];
+ });
+
+ REQUIRE(same_cell_value(cell_value, i++,
+ dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+ }
+
+ REQUIRE(i == function_id_list.size());
+ }
+
+ SECTION("3D")
+ {
+ constexpr size_t Dimension = 3;
+
+ std::shared_ptr<const IQuadratureDescriptor> quadrature_descriptor = std::make_shared<GaussQuadratureDescriptor>(3);
+
+ auto mesh_3d = MeshDataBaseForTests::get().hybrid3DMesh();
+
+ std::vector<std::shared_ptr<const IZoneDescriptor>> zone_list;
+ zone_list.push_back(std::make_shared<NamedZoneDescriptor>("LEFT"));
+
+ auto mesh_cell_zone = getMeshCellZone(*mesh_3d, *zone_list[0]);
+ auto zone_cell_list = mesh_cell_zone.cellList();
+
+ std::string_view data = R"(
+import math;
+let B_scalar_non_linear_3d: R^3 -> B, x -> (exp(2 * x[0] + x[2]) + 3 > 3.4);
+let N_scalar_non_linear_3d: R^3 -> N, x -> floor(3 * (x[0] * x[1]) * (x[0] * x[1]) + 2);
+let Z_scalar_non_linear_3d: R^3 -> Z, x -> floor(exp(2 * x[1]) - x[2]);
+let R_scalar_non_linear_3d: R^3 -> R, x -> 2 * exp(x[0] + x[1]) + 3 * 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};
+
+ 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::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+ std::vector<FunctionSymbolId> function_id_list;
+ register_function(position, symbol_table, "B_scalar_non_linear_3d", function_id_list);
+ register_function(position, symbol_table, "N_scalar_non_linear_3d", function_id_list);
+ register_function(position, symbol_table, "Z_scalar_non_linear_3d", function_id_list);
+ register_function(position, symbol_table, "R_scalar_non_linear_3d", function_id_list);
+
+ DiscreteFunctionVectorIntegrator integrator(mesh_3d, zone_list, quadrature_descriptor,
+ std::make_shared<DiscreteFunctionDescriptorP0Vector>(),
+ function_id_list);
+ std::shared_ptr discrete_function = integrator.integrate();
+
+ size_t i = 0;
+
+ {
+ Array<double> array =
+ IntegrateCellValue<double(TinyVector<3>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_3d,
+ zone_cell_list);
+
+ CellValue<double> cell_value{mesh_3d->connectivity()};
+ cell_value.fill(0);
+
+ parallel_for(
+ zone_cell_list.size(), PUGS_LAMBDA(const size_t j) {
+ const CellId cell_id = zone_cell_list[j];
+ cell_value[cell_id] = array[j];
+ });
+
+ REQUIRE(same_cell_value(cell_value, i++,
+ dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+ }
+
+ {
+ Array<double> array =
+ IntegrateCellValue<double(TinyVector<3>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_3d,
+ zone_cell_list);
+
+ CellValue<double> cell_value{mesh_3d->connectivity()};
+ cell_value.fill(0);
+
+ parallel_for(
+ zone_cell_list.size(), PUGS_LAMBDA(const size_t j) {
+ const CellId cell_id = zone_cell_list[j];
+ cell_value[cell_id] = array[j];
+ });
+
+ REQUIRE(same_cell_value(cell_value, i++,
+ dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+ }
+
+ {
+ Array<double> array =
+ IntegrateCellValue<double(TinyVector<3>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_3d,
+ zone_cell_list);
+
+ CellValue<double> cell_value{mesh_3d->connectivity()};
+ cell_value.fill(0);
+
+ parallel_for(
+ zone_cell_list.size(), PUGS_LAMBDA(const size_t j) {
+ const CellId cell_id = zone_cell_list[j];
+ cell_value[cell_id] = array[j];
+ });
+
+ REQUIRE(same_cell_value(cell_value, i++,
+ dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+ }
+
+ {
+ Array<double> array =
+ IntegrateCellValue<double(TinyVector<3>)>::integrate(function_id_list[i], *quadrature_descriptor, *mesh_3d,
+ zone_cell_list);
+
+ CellValue<double> cell_value{mesh_3d->connectivity()};
+ cell_value.fill(0);
+
+ parallel_for(
+ zone_cell_list.size(), PUGS_LAMBDA(const size_t j) {
+ const CellId cell_id = zone_cell_list[j];
+ cell_value[cell_id] = array[j];
+ });
+
+ REQUIRE(same_cell_value(cell_value, i++,
+ dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+ }
+
+ REQUIRE(i == function_id_list.size());
+ }
+}