diff --git a/src/scheme/DiscreteFunctionIntegrator.cpp b/src/scheme/DiscreteFunctionIntegrator.cpp index c6a950a9ddc7f088996425da48a0affdcf61d233..93c836461ce14ae05b51b879cf90525f73a793eb 100644 --- a/src/scheme/DiscreteFunctionIntegrator.cpp +++ b/src/scheme/DiscreteFunctionIntegrator.cpp @@ -46,10 +46,10 @@ DiscreteFunctionIntegrator::_integrateOnZoneList() const } } - Array<DataType> array = - IntegrateCellValue<DataType(TinyVector<Dimension>)>::template integrate<MeshType>(m_function_id, - *m_quadrature_descriptor, *p_mesh, - zone_cell_list); + Array<ValueType> array = + IntegrateCellValue<ValueType(TinyVector<Dimension>)>::template integrate<MeshType>(m_function_id, + *m_quadrature_descriptor, + *p_mesh, zone_cell_list); CellValue<ValueType> cell_value{p_mesh->connectivity()}; if constexpr (is_tiny_vector_v<ValueType> or is_tiny_matrix_v<ValueType>) { diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 2fd2c19e0bc610f8d82cc355ea1965eb486c7b3b..5035890f67986977626aaf13dbae968269aa66b6 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -158,6 +158,7 @@ add_executable (mpi_unit_tests mpi_test_main.cpp test_Connectivity.cpp test_DiscreteFunctionIntegrator.cpp + test_DiscreteFunctionIntegratorByZone.cpp test_DiscreteFunctionInterpoler.cpp test_DiscreteFunctionInterpolerByZone.cpp test_DiscreteFunctionP0.cpp diff --git a/tests/test_DiscreteFunctionIntegratorByZone.cpp b/tests/test_DiscreteFunctionIntegratorByZone.cpp new file mode 100644 index 0000000000000000000000000000000000000000..071f8cddb7e9f71e0b102bf361afb8715a7afe87 --- /dev/null +++ b/tests/test_DiscreteFunctionIntegratorByZone.cpp @@ -0,0 +1,1062 @@ +#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 <analysis/GaussQuadratureDescriptor.hpp> +#include <language/utils/IntegrateCellValue.hpp> +#include <scheme/DiscreteFunctionDescriptorP0.hpp> +#include <scheme/DiscreteFunctionIntegrator.hpp> +#include <scheme/DiscreteFunctionP0.hpp> + +#include <pegtl/string_input.hpp> + +// clazy:excludeall=non-pod-global-static + +TEST_CASE("DiscreteFunctionIntegratorByZone", "[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; + + 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; +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); + + Array<double> array = + IntegrateCellValue<double(TinyVector<1>)>::integrate(function_symbol_id, *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 i) { + const CellId cell_id = zone_cell_list[i]; + cell_value[cell_id] = array[i]; + }); + + DiscreteFunctionIntegrator integrator(mesh_1d, zone_list, quadrature_descriptor, function_symbol_id); + std::shared_ptr discrete_function = integrator.integrate(); + + 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); + + Array<double> array = + IntegrateCellValue<double(TinyVector<1>)>::integrate(function_symbol_id, *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 i) { + const CellId cell_id = zone_cell_list[i]; + cell_value[cell_id] = array[i]; + }); + + DiscreteFunctionIntegrator integrator(mesh_1d, zone_list, quadrature_descriptor, function_symbol_id); + std::shared_ptr discrete_function = integrator.integrate(); + + 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); + + Array<double> array = + IntegrateCellValue<double(TinyVector<1>)>::integrate(function_symbol_id, *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 i) { + const CellId cell_id = zone_cell_list[i]; + cell_value[cell_id] = array[i]; + }); + + DiscreteFunctionIntegrator integrator(mesh_1d, zone_list, quadrature_descriptor, function_symbol_id); + std::shared_ptr discrete_function = integrator.integrate(); + + 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); + + Array<double> array = + IntegrateCellValue<double(TinyVector<1>)>::integrate(function_symbol_id, *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 i) { + const CellId cell_id = zone_cell_list[i]; + cell_value[cell_id] = array[i]; + }); + + DiscreteFunctionIntegrator integrator(mesh_1d, zone_list, quadrature_descriptor, function_symbol_id); + std::shared_ptr discrete_function = integrator.integrate(); + + 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); + + Array<DataType> array = + IntegrateCellValue<DataType(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_1d, + zone_cell_list); + + CellValue<DataType> cell_value{mesh_1d->connectivity()}; + cell_value.fill(zero); + + parallel_for( + zone_cell_list.size(), PUGS_LAMBDA(const size_t i) { + const CellId cell_id = zone_cell_list[i]; + cell_value[cell_id] = array[i]; + }); + + DiscreteFunctionIntegrator integrator(mesh_1d, zone_list, quadrature_descriptor, function_symbol_id); + std::shared_ptr discrete_function = integrator.integrate(); + + 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); + + Array<DataType> array = + IntegrateCellValue<DataType(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_1d, + zone_cell_list); + + CellValue<DataType> cell_value{mesh_1d->connectivity()}; + cell_value.fill(zero); + + parallel_for( + zone_cell_list.size(), PUGS_LAMBDA(const size_t i) { + const CellId cell_id = zone_cell_list[i]; + cell_value[cell_id] = array[i]; + }); + + DiscreteFunctionIntegrator integrator(mesh_1d, zone_list, quadrature_descriptor, function_symbol_id); + std::shared_ptr discrete_function = integrator.integrate(); + + 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); + + Array<DataType> array = + IntegrateCellValue<DataType(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_1d, + zone_cell_list); + + CellValue<DataType> cell_value{mesh_1d->connectivity()}; + cell_value.fill(zero); + + parallel_for( + zone_cell_list.size(), PUGS_LAMBDA(const size_t i) { + const CellId cell_id = zone_cell_list[i]; + cell_value[cell_id] = array[i]; + }); + + DiscreteFunctionIntegrator integrator(mesh_1d, zone_list, quadrature_descriptor, function_symbol_id); + std::shared_ptr discrete_function = integrator.integrate(); + + 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); + + Array<DataType> array = + IntegrateCellValue<DataType(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_1d, + zone_cell_list); + + CellValue<DataType> cell_value{mesh_1d->connectivity()}; + cell_value.fill(zero); + + parallel_for( + zone_cell_list.size(), PUGS_LAMBDA(const size_t i) { + const CellId cell_id = zone_cell_list[i]; + cell_value[cell_id] = array[i]; + }); + + DiscreteFunctionIntegrator integrator(mesh_1d, zone_list, quadrature_descriptor, function_symbol_id); + std::shared_ptr discrete_function = integrator.integrate(); + + 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); + + Array<DataType> array = + IntegrateCellValue<DataType(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_1d, + zone_cell_list); + + CellValue<DataType> cell_value{mesh_1d->connectivity()}; + cell_value.fill(zero); + + parallel_for( + zone_cell_list.size(), PUGS_LAMBDA(const size_t i) { + const CellId cell_id = zone_cell_list[i]; + cell_value[cell_id] = array[i]; + }); + + DiscreteFunctionIntegrator integrator(mesh_1d, zone_list, quadrature_descriptor, function_symbol_id); + std::shared_ptr discrete_function = integrator.integrate(); + + 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); + + Array<DataType> array = + IntegrateCellValue<DataType(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_1d, + zone_cell_list); + + CellValue<DataType> cell_value{mesh_1d->connectivity()}; + cell_value.fill(zero); + + parallel_for( + zone_cell_list.size(), PUGS_LAMBDA(const size_t i) { + const CellId cell_id = zone_cell_list[i]; + cell_value[cell_id] = array[i]; + }); + + DiscreteFunctionIntegrator integrator(mesh_1d, zone_list, quadrature_descriptor, function_symbol_id); + std::shared_ptr discrete_function = integrator.integrate(); + + REQUIRE( + same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function))); + } + } + + 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])< 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); + + Array<double> array = + IntegrateCellValue<double(TinyVector<2>)>::integrate(function_symbol_id, *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 i) { + const CellId cell_id = zone_cell_list[i]; + cell_value[cell_id] = array[i]; + }); + + DiscreteFunctionIntegrator integrator(mesh_2d, zone_list, quadrature_descriptor, function_symbol_id); + std::shared_ptr discrete_function = integrator.integrate(); + + 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); + + Array<double> array = + IntegrateCellValue<double(TinyVector<2>)>::integrate(function_symbol_id, *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 i) { + const CellId cell_id = zone_cell_list[i]; + cell_value[cell_id] = array[i]; + }); + + DiscreteFunctionIntegrator integrator(mesh_2d, zone_list, quadrature_descriptor, function_symbol_id); + std::shared_ptr discrete_function = integrator.integrate(); + + 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); + + Array<double> array = + IntegrateCellValue<double(TinyVector<2>)>::integrate(function_symbol_id, *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 i) { + const CellId cell_id = zone_cell_list[i]; + cell_value[cell_id] = array[i]; + }); + + DiscreteFunctionIntegrator integrator(mesh_2d, zone_list, quadrature_descriptor, function_symbol_id); + std::shared_ptr discrete_function = integrator.integrate(); + + 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); + + Array<double> array = + IntegrateCellValue<double(TinyVector<2>)>::integrate(function_symbol_id, *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 i) { + const CellId cell_id = zone_cell_list[i]; + cell_value[cell_id] = array[i]; + }); + + DiscreteFunctionIntegrator integrator(mesh_2d, zone_list, quadrature_descriptor, function_symbol_id); + std::shared_ptr discrete_function = integrator.integrate(); + + 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); + + Array<DataType> array = + IntegrateCellValue<DataType(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_2d, + zone_cell_list); + + CellValue<DataType> cell_value{mesh_2d->connectivity()}; + cell_value.fill(zero); + + parallel_for( + zone_cell_list.size(), PUGS_LAMBDA(const size_t i) { + const CellId cell_id = zone_cell_list[i]; + cell_value[cell_id] = array[i]; + }); + + DiscreteFunctionIntegrator integrator(mesh_2d, zone_list, quadrature_descriptor, function_symbol_id); + std::shared_ptr discrete_function = integrator.integrate(); + + 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); + + Array<DataType> array = + IntegrateCellValue<DataType(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_2d, + zone_cell_list); + + CellValue<DataType> cell_value{mesh_2d->connectivity()}; + cell_value.fill(zero); + + parallel_for( + zone_cell_list.size(), PUGS_LAMBDA(const size_t i) { + const CellId cell_id = zone_cell_list[i]; + cell_value[cell_id] = array[i]; + }); + + DiscreteFunctionIntegrator integrator(mesh_2d, zone_list, quadrature_descriptor, function_symbol_id); + std::shared_ptr discrete_function = integrator.integrate(); + + 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); + + Array<DataType> array = + IntegrateCellValue<DataType(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_2d, + zone_cell_list); + + CellValue<DataType> cell_value{mesh_2d->connectivity()}; + cell_value.fill(zero); + + parallel_for( + zone_cell_list.size(), PUGS_LAMBDA(const size_t i) { + const CellId cell_id = zone_cell_list[i]; + cell_value[cell_id] = array[i]; + }); + + DiscreteFunctionIntegrator integrator(mesh_2d, zone_list, quadrature_descriptor, function_symbol_id); + std::shared_ptr discrete_function = integrator.integrate(); + + 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); + + Array<DataType> array = + IntegrateCellValue<DataType(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_2d, + zone_cell_list); + + CellValue<DataType> cell_value{mesh_2d->connectivity()}; + cell_value.fill(zero); + + parallel_for( + zone_cell_list.size(), PUGS_LAMBDA(const size_t i) { + const CellId cell_id = zone_cell_list[i]; + cell_value[cell_id] = array[i]; + }); + + DiscreteFunctionIntegrator integrator(mesh_2d, zone_list, quadrature_descriptor, function_symbol_id); + std::shared_ptr discrete_function = integrator.integrate(); + + 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); + + Array<DataType> array = + IntegrateCellValue<DataType(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_2d, + zone_cell_list); + + CellValue<DataType> cell_value{mesh_2d->connectivity()}; + cell_value.fill(zero); + + parallel_for( + zone_cell_list.size(), PUGS_LAMBDA(const size_t i) { + const CellId cell_id = zone_cell_list[i]; + cell_value[cell_id] = array[i]; + }); + + DiscreteFunctionIntegrator integrator(mesh_2d, zone_list, quadrature_descriptor, function_symbol_id); + std::shared_ptr discrete_function = integrator.integrate(); + + 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); + + Array<DataType> array = + IntegrateCellValue<DataType(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_2d, + zone_cell_list); + + CellValue<DataType> cell_value{mesh_2d->connectivity()}; + cell_value.fill(zero); + + parallel_for( + zone_cell_list.size(), PUGS_LAMBDA(const size_t i) { + const CellId cell_id = zone_cell_list[i]; + cell_value[cell_id] = array[i]; + }); + + DiscreteFunctionIntegrator integrator(mesh_2d, zone_list, quadrature_descriptor, function_symbol_id); + std::shared_ptr discrete_function = integrator.integrate(); + + REQUIRE( + same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function))); + } + } + + 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])< 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); + + Array<double> array = + IntegrateCellValue<double(TinyVector<3>)>::integrate(function_symbol_id, *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 i) { + const CellId cell_id = zone_cell_list[i]; + cell_value[cell_id] = array[i]; + }); + + DiscreteFunctionIntegrator integrator(mesh_3d, zone_list, quadrature_descriptor, function_symbol_id); + std::shared_ptr discrete_function = integrator.integrate(); + + 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); + + Array<double> array = + IntegrateCellValue<double(TinyVector<3>)>::integrate(function_symbol_id, *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 i) { + const CellId cell_id = zone_cell_list[i]; + cell_value[cell_id] = array[i]; + }); + + DiscreteFunctionIntegrator integrator(mesh_3d, zone_list, quadrature_descriptor, function_symbol_id); + std::shared_ptr discrete_function = integrator.integrate(); + + 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); + + Array<double> array = + IntegrateCellValue<double(TinyVector<3>)>::integrate(function_symbol_id, *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 i) { + const CellId cell_id = zone_cell_list[i]; + cell_value[cell_id] = array[i]; + }); + + DiscreteFunctionIntegrator integrator(mesh_3d, zone_list, quadrature_descriptor, function_symbol_id); + std::shared_ptr discrete_function = integrator.integrate(); + + 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); + + Array<double> array = + IntegrateCellValue<double(TinyVector<3>)>::integrate(function_symbol_id, *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 i) { + const CellId cell_id = zone_cell_list[i]; + cell_value[cell_id] = array[i]; + }); + + DiscreteFunctionIntegrator integrator(mesh_3d, zone_list, quadrature_descriptor, function_symbol_id); + std::shared_ptr discrete_function = integrator.integrate(); + + 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); + + Array<DataType> array = + IntegrateCellValue<DataType(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_3d, + zone_cell_list); + + CellValue<DataType> cell_value{mesh_3d->connectivity()}; + cell_value.fill(zero); + + parallel_for( + zone_cell_list.size(), PUGS_LAMBDA(const size_t i) { + const CellId cell_id = zone_cell_list[i]; + cell_value[cell_id] = array[i]; + }); + + DiscreteFunctionIntegrator integrator(mesh_3d, zone_list, quadrature_descriptor, function_symbol_id); + std::shared_ptr discrete_function = integrator.integrate(); + + 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); + + Array<DataType> array = + IntegrateCellValue<DataType(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_3d, + zone_cell_list); + + CellValue<DataType> cell_value{mesh_3d->connectivity()}; + cell_value.fill(zero); + + parallel_for( + zone_cell_list.size(), PUGS_LAMBDA(const size_t i) { + const CellId cell_id = zone_cell_list[i]; + cell_value[cell_id] = array[i]; + }); + + DiscreteFunctionIntegrator integrator(mesh_3d, zone_list, quadrature_descriptor, function_symbol_id); + std::shared_ptr discrete_function = integrator.integrate(); + + 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); + + Array<DataType> array = + IntegrateCellValue<DataType(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_3d, + zone_cell_list); + + CellValue<DataType> cell_value{mesh_3d->connectivity()}; + cell_value.fill(zero); + + parallel_for( + zone_cell_list.size(), PUGS_LAMBDA(const size_t i) { + const CellId cell_id = zone_cell_list[i]; + cell_value[cell_id] = array[i]; + }); + + DiscreteFunctionIntegrator integrator(mesh_3d, zone_list, quadrature_descriptor, function_symbol_id); + std::shared_ptr discrete_function = integrator.integrate(); + + 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); + + Array<DataType> array = + IntegrateCellValue<DataType(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_3d, + zone_cell_list); + + CellValue<DataType> cell_value{mesh_3d->connectivity()}; + cell_value.fill(zero); + + parallel_for( + zone_cell_list.size(), PUGS_LAMBDA(const size_t i) { + const CellId cell_id = zone_cell_list[i]; + cell_value[cell_id] = array[i]; + }); + + DiscreteFunctionIntegrator integrator(mesh_3d, zone_list, quadrature_descriptor, function_symbol_id); + std::shared_ptr discrete_function = integrator.integrate(); + + 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); + + Array<DataType> array = + IntegrateCellValue<DataType(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_3d, + zone_cell_list); + + CellValue<DataType> cell_value{mesh_3d->connectivity()}; + cell_value.fill(zero); + + parallel_for( + zone_cell_list.size(), PUGS_LAMBDA(const size_t i) { + const CellId cell_id = zone_cell_list[i]; + cell_value[cell_id] = array[i]; + }); + + DiscreteFunctionIntegrator integrator(mesh_3d, zone_list, quadrature_descriptor, function_symbol_id); + std::shared_ptr discrete_function = integrator.integrate(); + + 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); + + Array<DataType> array = + IntegrateCellValue<DataType(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_3d, + zone_cell_list); + + CellValue<DataType> cell_value{mesh_3d->connectivity()}; + cell_value.fill(zero); + + parallel_for( + zone_cell_list.size(), PUGS_LAMBDA(const size_t i) { + const CellId cell_id = zone_cell_list[i]; + cell_value[cell_id] = array[i]; + }); + + DiscreteFunctionIntegrator integrator(mesh_3d, zone_list, quadrature_descriptor, function_symbol_id); + std::shared_ptr discrete_function = integrator.integrate(); + + REQUIRE( + same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function))); + } + } +}