Skip to content
Snippets Groups Projects
Select Git revision
  • 4b59dbba0c955f9bca1b556a9d097d3a8f2172a7
  • develop default protected
  • feature/advection
  • feature/composite-scheme-other-fluxes
  • origin/stage/bouguettaia
  • save_clemence
  • feature/local-dt-fsi
  • feature/variational-hydro
  • feature/gmsh-reader
  • feature/reconstruction
  • feature/kinetic-schemes
  • feature/composite-scheme-sources
  • feature/serraille
  • feature/composite-scheme
  • hyperplastic
  • feature/polynomials
  • feature/gks
  • feature/implicit-solver-o2
  • feature/coupling_module
  • feature/implicit-solver
  • feature/merge-local-dt-fsi
  • v0.5.0 protected
  • v0.4.1 protected
  • v0.4.0 protected
  • v0.3.0 protected
  • v0.2.0 protected
  • v0.1.0 protected
  • Kidder
  • v0.0.4 protected
  • v0.0.3 protected
  • v0.0.2 protected
  • v0 protected
  • v0.0.1 protected
33 results

MeshDataBaseForTests.cpp

Blame
  • test_DiscreteFunctionInterpolerByZone.cpp 43.96 KiB
    #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/Mesh.hpp>
    #include <mesh/MeshCellZone.hpp>
    #include <mesh/MeshData.hpp>
    #include <mesh/MeshDataManager.hpp>
    #include <mesh/NamedZoneDescriptor.hpp>
    
    #include <scheme/DiscreteFunctionDescriptorP0.hpp>
    #include <scheme/DiscreteFunctionInterpoler.hpp>
    #include <scheme/DiscreteFunctionP0.hpp>
    
    #include <pegtl/string_input.hpp>
    
    // clazy:excludeall=non-pod-global-static
    
    TEST_CASE("DiscreteFunctionInterpolerByZone", "[scheme]")
    {
      auto same_cell_value = [](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 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]);
        CellValue<bool> is_cell_in_zone{mesh_1d->connectivity()};
        is_cell_in_zone.fill(false);
        auto zone_cell_list = mesh_cell_zone.cellList();
        for (size_t i_cell = 0; i_cell < zone_cell_list.size(); ++i_cell) {
          is_cell_in_zone[zone_cell_list[i_cell]] = true;
        }
    
        auto xj = MeshDataManager::instance().getMeshData(*mesh_1d).xj();
    
        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;
    let R1_non_linear_1d: R^1 -> R^1, x -> 2 * exp(x[0]);
    let R2_non_linear_1d: R^1 -> R^2, x -> (2 * exp(x[0]), -3*x[0]);
    let R3_non_linear_1d: R^1 -> R^3, x -> (2 * exp(x[0]) + 3, x[0] - 2, 3);
    let R1x1_non_linear_1d: R^1 -> R^1x1, x -> (2 * exp(x[0]) * sin(x[0]) + 3);
    let R2x2_non_linear_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]);
    let R3x3_non_linear_1d: R^1 -> R^3x3, x -> (2 * exp(x[0]) * sin(x[0]) + 3, sin(x[0] - 2 * x[0]), 3, x[0] * x[0], -4*x[0], 2*x[0]+1, 3, -6*x[0], exp(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
    
        SECTION("B_scalar_non_linear_1d")
        {
          auto [i_symbol, found] = symbol_table->find("B_scalar_non_linear_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<double> cell_value{mesh_1d->connectivity()};
          parallel_for(
            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
              if (is_cell_in_zone[cell_id]) {
                const TinyVector<Dimension>& x = xj[cell_id];
                cell_value[cell_id]            = std::exp(2 * x[0]) + 3 < 3.3;
              } else {
                cell_value[cell_id] = false;
              }
            });
    
          DiscreteFunctionInterpoler interpoler(mesh_1d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
                                                function_symbol_id);
          std::shared_ptr discrete_function = interpoler.interpolate();
    
          REQUIRE(
            same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
        }
    
        SECTION("N_scalar_non_linear_1d")
        {
          auto [i_symbol, found] = symbol_table->find("N_scalar_non_linear_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<double> cell_value{mesh_1d->connectivity()};
          parallel_for(
            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
              if (is_cell_in_zone[cell_id]) {
                const TinyVector<Dimension>& x = xj[cell_id];
                cell_value[cell_id]            = std::floor(3 * x[0] * x[0] + 2);
              } else {
                cell_value[cell_id] = 0;
              }
            });
    
          DiscreteFunctionInterpoler interpoler(mesh_1d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
                                                function_symbol_id);
          std::shared_ptr discrete_function = interpoler.interpolate();
    
          REQUIRE(
            same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
        }
    
        SECTION("Z_scalar_non_linear_1d")
        {
          auto [i_symbol, found] = symbol_table->find("Z_scalar_non_linear_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<double> cell_value{mesh_1d->connectivity()};
          parallel_for(
            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
              if (is_cell_in_zone[cell_id]) {
                const TinyVector<Dimension>& x = xj[cell_id];
                cell_value[cell_id]            = std::floor(std::exp(2 * x[0]) - 1);
              } else {
                cell_value[cell_id] = 0;
              }
            });
    
          DiscreteFunctionInterpoler interpoler(mesh_1d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
                                                function_symbol_id);
          std::shared_ptr discrete_function = interpoler.interpolate();
    
          REQUIRE(
            same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
        }
    
        SECTION("R_scalar_non_linear_1d")
        {
          auto [i_symbol, found] = symbol_table->find("R_scalar_non_linear_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<double> cell_value{mesh_1d->connectivity()};
          parallel_for(
            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
              if (is_cell_in_zone[cell_id]) {
                const TinyVector<Dimension>& x = xj[cell_id];
                cell_value[cell_id]            = 2 * std::exp(x[0]) + 3;
              } else {
                cell_value[cell_id] = 0;
              }
            });
    
          DiscreteFunctionInterpoler interpoler(mesh_1d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
                                                function_symbol_id);
          std::shared_ptr discrete_function = interpoler.interpolate();
    
          REQUIRE(
            same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
        }
    
        SECTION("R1_non_linear_1d")
        {
          using DataType = TinyVector<1>;
    
          auto [i_symbol, found] = symbol_table->find("R1_non_linear_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<DataType> cell_value{mesh_1d->connectivity()};
          parallel_for(
            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
              if (is_cell_in_zone[cell_id]) {
                const TinyVector<Dimension>& x = xj[cell_id];
                cell_value[cell_id]            = DataType{2 * std::exp(x[0])};
              } else {
                cell_value[cell_id] = zero;
              }
            });
    
          DiscreteFunctionInterpoler interpoler(mesh_1d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
                                                function_symbol_id);
          std::shared_ptr discrete_function = interpoler.interpolate();
    
          REQUIRE(
            same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
        }
    
        SECTION("R2_non_linear_1d")
        {
          using DataType = TinyVector<2>;
    
          auto [i_symbol, found] = symbol_table->find("R2_non_linear_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<DataType> cell_value{mesh_1d->connectivity()};
          parallel_for(
            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
              if (is_cell_in_zone[cell_id]) {
                const TinyVector<Dimension>& x = xj[cell_id];
                cell_value[cell_id]            = DataType{2 * std::exp(x[0]), -3 * x[0]};
              } else {
                cell_value[cell_id] = zero;
              }
            });
    
          DiscreteFunctionInterpoler interpoler(mesh_1d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
                                                function_symbol_id);
          std::shared_ptr discrete_function = interpoler.interpolate();
    
          REQUIRE(
            same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
        }
    
        SECTION("R3_non_linear_1d")
        {
          using DataType = TinyVector<3>;
    
          auto [i_symbol, found] = symbol_table->find("R3_non_linear_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<DataType> cell_value{mesh_1d->connectivity()};
          parallel_for(
            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
              if (is_cell_in_zone[cell_id]) {
                const TinyVector<Dimension>& x = xj[cell_id];
                cell_value[cell_id]            = DataType{2 * std::exp(x[0]) + 3, x[0] - 2, 3};
              } else {
                cell_value[cell_id] = zero;
              }
            });
    
          DiscreteFunctionInterpoler interpoler(mesh_1d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
                                                function_symbol_id);
          std::shared_ptr discrete_function = interpoler.interpolate();
    
          REQUIRE(
            same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
        }
    
        SECTION("R1x1_non_linear_1d")
        {
          using DataType = TinyMatrix<1>;
    
          auto [i_symbol, found] = symbol_table->find("R1x1_non_linear_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<DataType> cell_value{mesh_1d->connectivity()};
          parallel_for(
            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
              if (is_cell_in_zone[cell_id]) {
                const TinyVector<Dimension>& x = xj[cell_id];
                cell_value[cell_id]            = DataType{2 * std::exp(x[0]) * std::sin(x[0]) + 3};
              } else {
                cell_value[cell_id] = zero;
              }
            });
    
          DiscreteFunctionInterpoler interpoler(mesh_1d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
                                                function_symbol_id);
          std::shared_ptr discrete_function = interpoler.interpolate();
    
          REQUIRE(
            same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
        }
    
        SECTION("R2x2_non_linear_1d")
        {
          using DataType = TinyMatrix<2>;
    
          auto [i_symbol, found] = symbol_table->find("R2x2_non_linear_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<DataType> cell_value{mesh_1d->connectivity()};
          parallel_for(
            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
              if (is_cell_in_zone[cell_id]) {
                const TinyVector<Dimension>& x = xj[cell_id];
                cell_value[cell_id] =
                  DataType{2 * std::exp(x[0]) * std::sin(x[0]) + 3, std::sin(x[0] - 2 * x[0]), 3, x[0] * x[0]};
              } else {
                cell_value[cell_id] = zero;
              }
            });
    
          DiscreteFunctionInterpoler interpoler(mesh_1d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
                                                function_symbol_id);
          std::shared_ptr discrete_function = interpoler.interpolate();
    
          REQUIRE(
            same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
        }
    
        SECTION("R3x3_non_linear_1d")
        {
          using DataType = TinyMatrix<3>;
    
          auto [i_symbol, found] = symbol_table->find("R3x3_non_linear_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<DataType> cell_value{mesh_1d->connectivity()};
          parallel_for(
            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
              if (is_cell_in_zone[cell_id]) {
                const TinyVector<Dimension>& x = xj[cell_id];
                cell_value[cell_id]            = DataType{2 * exp(x[0]) * std::sin(x[0]) + 3,
                                               std::sin(x[0] - 2 * x[0]),
                                               3,
                                               x[0] * x[0],
                                               -4 * x[0],
                                               2 * x[0] + 1,
                                               3,
                                               -6 * x[0],
                                               std::exp(x[0])};
              } else {
                cell_value[cell_id] = zero;
              }
            });
    
          DiscreteFunctionInterpoler interpoler(mesh_1d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
                                                function_symbol_id);
          std::shared_ptr discrete_function = interpoler.interpolate();
    
          REQUIRE(
            same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
        }
      }
    
      SECTION("2D")
      {
        constexpr size_t Dimension = 2;
    
        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]);
        CellValue<bool> is_cell_in_zone{mesh_2d->connectivity()};
        is_cell_in_zone.fill(false);
        auto zone_cell_list = mesh_cell_zone.cellList();
        for (size_t i_cell = 0; i_cell < zone_cell_list.size(); ++i_cell) {
          is_cell_in_zone[zone_cell_list[i_cell]] = true;
        }
    
        auto xj = MeshDataManager::instance().getMeshData(*mesh_2d).xj();
    
        std::string_view data = R"(
    import math;
    let B_scalar_non_linear_2d: R^2 -> B, x -> (exp(2 * x[0])< 2*x[1]);
    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[0]) - 3 * x[1]);
    let R_scalar_non_linear_2d: R^2 -> R, x -> 2 * exp(x[0]) + 3 * x[1];
    let R1_non_linear_2d: R^2 -> R^1, x -> 2 * exp(x[0]);
    let R2_non_linear_2d: R^2 -> R^2, x -> (2 * exp(x[0]), -3*x[1]);
    let R3_non_linear_2d: R^2 -> R^3, x -> (2 * exp(x[0]) + 3, x[1] - 2, 3);
    let R1x1_non_linear_2d: R^2 -> R^1x1, x -> (2 * exp(x[0]) * sin(x[1]) + 3);
    let R2x2_non_linear_2d: R^2 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3, sin(x[1] - 2 * x[0]), 3, x[1] * x[0]);
    let R3x3_non_linear_2d: R^2 -> R^3x3, x -> (2 * exp(x[0]) * sin(x[1]) + 3, sin(x[1] - 2 * x[0]), 3, x[1] * x[0], -4*x[1], 2*x[0]+1, 3, -6*x[0], exp(x[1]));
    )";
        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
    
        SECTION("B_scalar_non_linear_2d")
        {
          auto [i_symbol, found] = symbol_table->find("B_scalar_non_linear_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<double> cell_value{mesh_2d->connectivity()};
          parallel_for(
            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
              if (is_cell_in_zone[cell_id]) {
                const TinyVector<Dimension>& x = xj[cell_id];
                cell_value[cell_id]            = std::exp(2 * x[0]) < 2 * x[1];
              } else {
                cell_value[cell_id] = false;
              }
            });
    
          DiscreteFunctionInterpoler interpoler(mesh_2d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
                                                function_symbol_id);
          std::shared_ptr discrete_function = interpoler.interpolate();
    
          REQUIRE(
            same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
        }
    
        SECTION("N_scalar_non_linear_2d")
        {
          auto [i_symbol, found] = symbol_table->find("N_scalar_non_linear_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<double> cell_value{mesh_2d->connectivity()};
          parallel_for(
            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
              if (is_cell_in_zone[cell_id]) {
                const TinyVector<Dimension>& x = xj[cell_id];
                cell_value[cell_id]            = std::floor(3 * (x[0] + x[1]) * (x[0] + x[1]) + 2);
              } else {
                cell_value[cell_id] = 0;
              }
            });
    
          DiscreteFunctionInterpoler interpoler(mesh_2d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
                                                function_symbol_id);
          std::shared_ptr discrete_function = interpoler.interpolate();
    
          REQUIRE(
            same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
        }
    
        SECTION("Z_scalar_non_linear_2d")
        {
          auto [i_symbol, found] = symbol_table->find("Z_scalar_non_linear_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<double> cell_value{mesh_2d->connectivity()};
          parallel_for(
            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
              if (is_cell_in_zone[cell_id]) {
                const TinyVector<Dimension>& x = xj[cell_id];
                cell_value[cell_id]            = std::floor(std::exp(2 * x[0]) - 3 * x[1]);
              } else {
                cell_value[cell_id] = 0;
              }
            });
    
          DiscreteFunctionInterpoler interpoler(mesh_2d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
                                                function_symbol_id);
          std::shared_ptr discrete_function = interpoler.interpolate();
    
          REQUIRE(
            same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
        }
    
        SECTION("R_scalar_non_linear_2d")
        {
          auto [i_symbol, found] = symbol_table->find("R_scalar_non_linear_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<double> cell_value{mesh_2d->connectivity()};
          parallel_for(
            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
              if (is_cell_in_zone[cell_id]) {
                const TinyVector<Dimension>& x = xj[cell_id];
                cell_value[cell_id]            = 2 * std::exp(x[0]) + 3 * x[1];
              } else {
                cell_value[cell_id] = 0;
              }
            });
    
          DiscreteFunctionInterpoler interpoler(mesh_2d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
                                                function_symbol_id);
          std::shared_ptr discrete_function = interpoler.interpolate();
    
          REQUIRE(
            same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
        }
    
        SECTION("R1_non_linear_2d")
        {
          using DataType = TinyVector<1>;
    
          auto [i_symbol, found] = symbol_table->find("R1_non_linear_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<DataType> cell_value{mesh_2d->connectivity()};
          parallel_for(
            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
              if (is_cell_in_zone[cell_id]) {
                const TinyVector<Dimension>& x = xj[cell_id];
                cell_value[cell_id]            = DataType{2 * std::exp(x[0])};
              } else {
                cell_value[cell_id] = zero;
              }
            });
    
          DiscreteFunctionInterpoler interpoler(mesh_2d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
                                                function_symbol_id);
          std::shared_ptr discrete_function = interpoler.interpolate();
    
          REQUIRE(
            same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
        }
    
        SECTION("R2_non_linear_2d")
        {
          using DataType = TinyVector<2>;
    
          auto [i_symbol, found] = symbol_table->find("R2_non_linear_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<DataType> cell_value{mesh_2d->connectivity()};
          parallel_for(
            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
              if (is_cell_in_zone[cell_id]) {
                const TinyVector<Dimension>& x = xj[cell_id];
                cell_value[cell_id]            = DataType{2 * std::exp(x[0]), -3 * x[1]};
              } else {
                cell_value[cell_id] = zero;
              }
            });
    
          DiscreteFunctionInterpoler interpoler(mesh_2d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
                                                function_symbol_id);
          std::shared_ptr discrete_function = interpoler.interpolate();
    
          REQUIRE(
            same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
        }
    
        SECTION("R3_non_linear_2d")
        {
          using DataType = TinyVector<3>;
    
          auto [i_symbol, found] = symbol_table->find("R3_non_linear_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<DataType> cell_value{mesh_2d->connectivity()};
          parallel_for(
            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
              if (is_cell_in_zone[cell_id]) {
                const TinyVector<Dimension>& x = xj[cell_id];
                cell_value[cell_id]            = DataType{2 * std::exp(x[0]) + 3, x[1] - 2, 3};
              } else {
                cell_value[cell_id] = zero;
              }
            });
    
          DiscreteFunctionInterpoler interpoler(mesh_2d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
                                                function_symbol_id);
          std::shared_ptr discrete_function = interpoler.interpolate();
    
          REQUIRE(
            same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
        }
    
        SECTION("R1x1_non_linear_2d")
        {
          using DataType = TinyMatrix<1>;
    
          auto [i_symbol, found] = symbol_table->find("R1x1_non_linear_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<DataType> cell_value{mesh_2d->connectivity()};
          parallel_for(
            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
              if (is_cell_in_zone[cell_id]) {
                const TinyVector<Dimension>& x = xj[cell_id];
                cell_value[cell_id]            = DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3};
              } else {
                cell_value[cell_id] = zero;
              }
            });
    
          DiscreteFunctionInterpoler interpoler(mesh_2d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
                                                function_symbol_id);
          std::shared_ptr discrete_function = interpoler.interpolate();
    
          REQUIRE(
            same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
        }
    
        SECTION("R2x2_non_linear_2d")
        {
          using DataType = TinyMatrix<2>;
    
          auto [i_symbol, found] = symbol_table->find("R2x2_non_linear_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<DataType> cell_value{mesh_2d->connectivity()};
          parallel_for(
            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
              if (is_cell_in_zone[cell_id]) {
                const TinyVector<Dimension>& x = xj[cell_id];
                cell_value[cell_id] =
                  DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3, std::sin(x[1] - 2 * x[0]), 3, x[1] * x[0]};
              } else {
                cell_value[cell_id] = zero;
              }
            });
    
          DiscreteFunctionInterpoler interpoler(mesh_2d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
                                                function_symbol_id);
          std::shared_ptr discrete_function = interpoler.interpolate();
    
          REQUIRE(
            same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
        }
    
        SECTION("R3x3_non_linear_2d")
        {
          using DataType = TinyMatrix<3>;
    
          auto [i_symbol, found] = symbol_table->find("R3x3_non_linear_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<DataType> cell_value{mesh_2d->connectivity()};
          parallel_for(
            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
              if (is_cell_in_zone[cell_id]) {
                const TinyVector<Dimension>& x = xj[cell_id];
    
                cell_value[cell_id] = DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3,
                                               std::sin(x[1] - 2 * x[0]),
                                               3,
                                               x[1] * x[0],
                                               -4 * x[1],
                                               2 * x[0] + 1,
                                               3,
                                               -6 * x[0],
                                               std::exp(x[1])};
              } else {
                cell_value[cell_id] = zero;
              }
            });
    
          DiscreteFunctionInterpoler interpoler(mesh_2d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
                                                function_symbol_id);
          std::shared_ptr discrete_function = interpoler.interpolate();
    
          REQUIRE(
            same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
        }
      }
    
      SECTION("3D")
      {
        constexpr size_t Dimension = 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]);
        CellValue<bool> is_cell_in_zone{mesh_3d->connectivity()};
        is_cell_in_zone.fill(false);
        auto zone_cell_list = mesh_cell_zone.cellList();
        for (size_t i_cell = 0; i_cell < zone_cell_list.size(); ++i_cell) {
          is_cell_in_zone[zone_cell_list[i_cell]] = true;
        }
    
        auto xj = MeshDataManager::instance().getMeshData(*mesh_3d).xj();
    
        std::string_view data = R"(
    import math;
    let B_scalar_non_linear_3d: R^3 -> B, x -> (exp(2 * x[0])< 2*x[1]+x[2]);
    let N_scalar_non_linear_3d: R^3 -> N, x -> floor(3 * (x[0] + x[1]) * (x[0] + x[1]) + x[2] * x[2]);
    let Z_scalar_non_linear_3d: R^3 -> Z, x -> floor(exp(2 * x[0]) - 3 * x[1] + x[2]);
    let R_scalar_non_linear_3d: R^3 -> R, x -> 2 * exp(x[0]+x[2]) + 3 * x[1];
    let R1_non_linear_3d: R^3 -> R^1, x -> 2 * exp(x[0])+sin(x[1] + x[2]);
    let R2_non_linear_3d: R^3 -> R^2, x -> (2 * exp(x[0]), -3*x[1] * x[2]);
    let R3_non_linear_3d: R^3 -> R^3, x -> (2 * exp(x[0]) + 3, x[1] - 2, 3 * x[2]);
    let R1x1_non_linear_3d: R^3 -> R^1x1, x -> (2 * exp(x[0]) * sin(x[1]) + 3 * x[2]);
    let R2x2_non_linear_3d: R^3 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3, sin(x[2] - 2 * x[0]), 3, x[1] * x[0] - x[2]);
    let R3x3_non_linear_3d: R^3 -> R^3x3, x -> (2 * exp(x[0]) * sin(x[1]) + 3, sin(x[1] - 2 * x[2]), 3, x[1] * x[2], -4*x[1], 2*x[2]+1, 3, -6*x[2], exp(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
    
        SECTION("B_scalar_non_linear_3d")
        {
          auto [i_symbol, found] = symbol_table->find("B_scalar_non_linear_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_value{mesh_3d->connectivity()};
          parallel_for(
            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
              if (is_cell_in_zone[cell_id]) {
                const TinyVector<Dimension>& x = xj[cell_id];
                cell_value[cell_id]            = std::exp(2 * x[0]) < 2 * x[1] + x[2];
              } else {
                cell_value[cell_id] = false;
              }
            });
    
          DiscreteFunctionInterpoler interpoler(mesh_3d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
                                                function_symbol_id);
          std::shared_ptr discrete_function = interpoler.interpolate();
    
          REQUIRE(
            same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
        }
    
        SECTION("N_scalar_non_linear_3d")
        {
          auto [i_symbol, found] = symbol_table->find("N_scalar_non_linear_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_value{mesh_3d->connectivity()};
          parallel_for(
            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
              if (is_cell_in_zone[cell_id]) {
                const TinyVector<Dimension>& x = xj[cell_id];
                cell_value[cell_id]            = std::floor(3 * (x[0] + x[1]) * (x[0] + x[1]) + x[2] * x[2]);
              } else {
                cell_value[cell_id] = 0;
              }
            });
    
          DiscreteFunctionInterpoler interpoler(mesh_3d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
                                                function_symbol_id);
          std::shared_ptr discrete_function = interpoler.interpolate();
    
          REQUIRE(
            same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
        }
    
        SECTION("Z_scalar_non_linear_3d")
        {
          auto [i_symbol, found] = symbol_table->find("Z_scalar_non_linear_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_value{mesh_3d->connectivity()};
          parallel_for(
            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
              if (is_cell_in_zone[cell_id]) {
                const TinyVector<Dimension>& x = xj[cell_id];
                cell_value[cell_id]            = std::floor(std::exp(2 * x[0]) - 3 * x[1] + x[2]);
              } else {
                cell_value[cell_id] = 0;
              }
            });
    
          DiscreteFunctionInterpoler interpoler(mesh_3d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
                                                function_symbol_id);
          std::shared_ptr discrete_function = interpoler.interpolate();
    
          REQUIRE(
            same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
        }
    
        SECTION("R_scalar_non_linear_3d")
        {
          auto [i_symbol, found] = symbol_table->find("R_scalar_non_linear_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_value{mesh_3d->connectivity()};
          parallel_for(
            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
              if (is_cell_in_zone[cell_id]) {
                const TinyVector<Dimension>& x = xj[cell_id];
                cell_value[cell_id]            = 2 * std::exp(x[0] + x[2]) + 3 * x[1];
              } else {
                cell_value[cell_id] = 0;
              }
            });
    
          DiscreteFunctionInterpoler interpoler(mesh_3d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
                                                function_symbol_id);
          std::shared_ptr discrete_function = interpoler.interpolate();
    
          REQUIRE(
            same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
        }
    
        SECTION("R1_non_linear_3d")
        {
          using DataType = TinyVector<1>;
    
          auto [i_symbol, found] = symbol_table->find("R1_non_linear_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<DataType> cell_value{mesh_3d->connectivity()};
          parallel_for(
            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
              if (is_cell_in_zone[cell_id]) {
                const TinyVector<Dimension>& x = xj[cell_id];
                cell_value[cell_id]            = DataType{2 * std::exp(x[0]) + std::sin(x[1] + x[2])};
              } else {
                cell_value[cell_id] = zero;
              }
            });
    
          DiscreteFunctionInterpoler interpoler(mesh_3d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
                                                function_symbol_id);
          std::shared_ptr discrete_function = interpoler.interpolate();
    
          REQUIRE(
            same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
        }
    
        SECTION("R2_non_linear_3d")
        {
          using DataType = TinyVector<2>;
    
          auto [i_symbol, found] = symbol_table->find("R2_non_linear_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<DataType> cell_value{mesh_3d->connectivity()};
          parallel_for(
            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
              if (is_cell_in_zone[cell_id]) {
                const TinyVector<Dimension>& x = xj[cell_id];
                cell_value[cell_id]            = DataType{2 * std::exp(x[0]), -3 * x[1] * x[2]};
              } else {
                cell_value[cell_id] = zero;
              }
            });
    
          DiscreteFunctionInterpoler interpoler(mesh_3d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
                                                function_symbol_id);
          std::shared_ptr discrete_function = interpoler.interpolate();
    
          REQUIRE(
            same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
        }
    
        SECTION("R3_non_linear_3d")
        {
          using DataType = TinyVector<3>;
    
          auto [i_symbol, found] = symbol_table->find("R3_non_linear_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<DataType> cell_value{mesh_3d->connectivity()};
          parallel_for(
            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
              if (is_cell_in_zone[cell_id]) {
                const TinyVector<Dimension>& x = xj[cell_id];
                cell_value[cell_id]            = DataType{2 * std::exp(x[0]) + 3, x[1] - 2, 3 * x[2]};
              } else {
                cell_value[cell_id] = zero;
              }
            });
    
          DiscreteFunctionInterpoler interpoler(mesh_3d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
                                                function_symbol_id);
          std::shared_ptr discrete_function = interpoler.interpolate();
    
          REQUIRE(
            same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
        }
    
        SECTION("R1x1_non_linear_3d")
        {
          using DataType = TinyMatrix<1>;
    
          auto [i_symbol, found] = symbol_table->find("R1x1_non_linear_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<DataType> cell_value{mesh_3d->connectivity()};
          parallel_for(
            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
              if (is_cell_in_zone[cell_id]) {
                const TinyVector<Dimension>& x = xj[cell_id];
                cell_value[cell_id]            = DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3 * x[2]};
              } else {
                cell_value[cell_id] = zero;
              }
            });
    
          DiscreteFunctionInterpoler interpoler(mesh_3d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
                                                function_symbol_id);
          std::shared_ptr discrete_function = interpoler.interpolate();
    
          REQUIRE(
            same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
        }
    
        SECTION("R2x2_non_linear_3d")
        {
          using DataType = TinyMatrix<2>;
    
          auto [i_symbol, found] = symbol_table->find("R2x2_non_linear_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<DataType> cell_value{mesh_3d->connectivity()};
          parallel_for(
            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
              if (is_cell_in_zone[cell_id]) {
                const TinyVector<Dimension>& x = xj[cell_id];
                cell_value[cell_id] =
                  DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3, std::sin(x[2] - 2 * x[0]), 3, x[1] * x[0] - x[2]};
              } else {
                cell_value[cell_id] = zero;
              }
            });
    
          DiscreteFunctionInterpoler interpoler(mesh_3d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
                                                function_symbol_id);
          std::shared_ptr discrete_function = interpoler.interpolate();
    
          REQUIRE(
            same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
        }
    
        SECTION("R3x3_non_linear_3d")
        {
          using DataType = TinyMatrix<3>;
    
          auto [i_symbol, found] = symbol_table->find("R3x3_non_linear_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<DataType> cell_value{mesh_3d->connectivity()};
          parallel_for(
            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
              if (is_cell_in_zone[cell_id]) {
                const TinyVector<Dimension>& x = xj[cell_id];
    
                cell_value[cell_id] = DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3,
                                               std::sin(x[1] - 2 * x[2]),
                                               3,
                                               x[1] * x[2],
                                               -4 * x[1],
                                               2 * x[2] + 1,
                                               3,
                                               -6 * x[2],
                                               std::exp(x[1] + x[2])};
              } else {
                cell_value[cell_id] = zero;
              }
            });
    
          DiscreteFunctionInterpoler interpoler(mesh_3d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
                                                function_symbol_id);
          std::shared_ptr discrete_function = interpoler.interpolate();
    
          REQUIRE(
            same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
        }
      }
    }