From dedcadc2b24099c68d1b2a2820febc0aa5b29280 Mon Sep 17 00:00:00 2001 From: MARMAJOU ISABELLE <id.counilh@wanadoo.fr> Date: Fri, 7 Mar 2025 16:22:36 +0100 Subject: [PATCH] Add more tests for IntegrateOnCells in the parabolic case --- tests/test_IntegrateOnCells_parabolic.cpp | 608 ++++++++++++++-------- 1 file changed, 379 insertions(+), 229 deletions(-) diff --git a/tests/test_IntegrateOnCells_parabolic.cpp b/tests/test_IntegrateOnCells_parabolic.cpp index 24ab9c75f..f003922ff 100644 --- a/tests/test_IntegrateOnCells_parabolic.cpp +++ b/tests/test_IntegrateOnCells_parabolic.cpp @@ -14,6 +14,7 @@ #include <language/utils/SymbolTable.hpp> #include <MeshDataBaseForTests.hpp> +#include <geometry/LineParabolicTransformation.hpp> #include <mesh/Connectivity.hpp> #include <mesh/DualMeshManager.hpp> #include <mesh/Mesh.hpp> @@ -59,373 +60,522 @@ TEST_CASE("IntegrateOnCells_parabolic", "[language]") return std::sqrt(error); }; - SECTION("Gauss quadrature") + SECTION("straight edges") { - auto quadrature_descriptor = GaussQuadratureDescriptor(12); - - SECTION("integrate on all cells") + SECTION("Gauss quadrature") { - SECTION("2D") + auto quadrature_descriptor = GaussQuadratureDescriptor(12); + + SECTION("integrate on all cells") { - constexpr size_t Dimension = 2; + SECTION("2D") + { + constexpr size_t Dimension = 2; - std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes(); + std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes(); - for (const auto& named_mesh : mesh_list) { - SECTION(named_mesh.name()) - { - auto mesh_2d_v = named_mesh.mesh(); - auto mesh_2d = mesh_2d_v->get<Mesh<2>>(); + for (const auto& named_mesh : mesh_list) { + SECTION(named_mesh.name()) + { + auto mesh_2d_v = named_mesh.mesh(); + auto mesh_2d = mesh_2d_v->get<Mesh<2>>(); - PolynomialMeshBuilder pb{mesh_2d_v, 2}; + PolynomialMeshBuilder pb{mesh_2d_v, 2}; - auto parabolic_mesh = pb.mesh()->get<PolynomialMesh<2>>(); + auto parabolic_mesh = pb.mesh()->get<PolynomialMesh<2>>(); - std::string_view data = R"( + std::string_view data = R"( import math; let scalar_2d: R^2 -> R, x -> 2*x[0]*x[1]+3; let R3_2d: R^2 -> R^3, x -> [2*x[0]*x[0] + x[1]*x[1] - 1, x[0]-2*x[1], 3]; let R2x2_2d: R^2 -> R^2x2, x -> [[2*x[0]*x[0]*x[1]+3, x[0]-2*x[1]], [3, x[0]*x[1]]]; )"; - TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"}; + TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"}; - auto ast = ASTBuilder::build(input); + auto ast = ASTBuilder::build(input); - ASTModulesImporter{*ast}; - ASTNodeTypeCleaner<language::import_instruction>{*ast}; + ASTModulesImporter{*ast}; + ASTNodeTypeCleaner<language::import_instruction>{*ast}; - ASTSymbolTableBuilder{*ast}; - ASTNodeDataTypeBuilder{*ast}; + ASTSymbolTableBuilder{*ast}; + ASTNodeDataTypeBuilder{*ast}; - ASTNodeTypeCleaner<language::var_declaration>{*ast}; - ASTNodeTypeCleaner<language::fct_declaration>{*ast}; - ASTNodeExpressionBuilder{*ast}; + ASTNodeTypeCleaner<language::var_declaration>{*ast}; + ASTNodeTypeCleaner<language::fct_declaration>{*ast}; + ASTNodeExpressionBuilder{*ast}; - std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table; + std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table; - // ensure that variables are declared at this point - TAO_PEGTL_NAMESPACE::position position{data.size(), 1, 1, "fixture"}; + // ensure that variables are declared at this point + TAO_PEGTL_NAMESPACE::position position{data.size(), 1, 1, "fixture"}; - SECTION("scalar 2d") - { - auto [i_symbol, found] = symbol_table->find("scalar_2d", position); - REQUIRE(found); - REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t); + SECTION("scalar 2d") + { + auto [i_symbol, found] = symbol_table->find("scalar_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); + FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table); - ASTNode::setStackDetails(false); + ASTNode::setStackDetails(false); - Array<double> integrate_value_polygonal(mesh_2d->numberOfCells()); - IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, - *mesh_2d, integrate_value_polygonal); + Array<double> integrate_value_polygonal(mesh_2d->numberOfCells()); + IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, + *mesh_2d, integrate_value_polygonal); - Array<double> integrate_value_parabolic(parabolic_mesh->numberOfCells()); - IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, - *parabolic_mesh, integrate_value_parabolic); + Array<double> integrate_value_parabolic(parabolic_mesh->numberOfCells()); + IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, + *parabolic_mesh, + integrate_value_parabolic); - ASTNode::setStackDetails(true); + ASTNode::setStackDetails(true); - REQUIRE(scalar_error(integrate_value_parabolic, integrate_value_polygonal) < 1E-12); - } + REQUIRE(scalar_error(integrate_value_parabolic, integrate_value_polygonal) < 1E-12); + } - SECTION("vector 2d") - { - using R3 = TinyVector<3>; - auto [i_symbol, found] = symbol_table->find("R3_2d", position); - REQUIRE(found); - REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t); + SECTION("vector 2d") + { + using R3 = TinyVector<3>; + auto [i_symbol, found] = symbol_table->find("R3_2d", position); + REQUIRE(found); + REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t); - FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table); + FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table); - ASTNode::setStackDetails(false); + ASTNode::setStackDetails(false); - Array<R3> integrate_value_polygonal(mesh_2d->numberOfCells()); - IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, - *mesh_2d, integrate_value_polygonal); + Array<R3> integrate_value_polygonal(mesh_2d->numberOfCells()); + IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, + *mesh_2d, integrate_value_polygonal); - Array<R3> integrate_value_parabolic(parabolic_mesh->numberOfCells()); - IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, - *parabolic_mesh, integrate_value_parabolic); + Array<R3> integrate_value_parabolic(parabolic_mesh->numberOfCells()); + IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, + *parabolic_mesh, integrate_value_parabolic); - ASTNode::setStackDetails(true); + ASTNode::setStackDetails(true); - REQUIRE(vector_error(integrate_value_parabolic, integrate_value_polygonal) < 1E-12); - } + REQUIRE(vector_error(integrate_value_parabolic, integrate_value_polygonal) < 1E-12); + } - SECTION("matrix 2d") - { - using R2x2 = TinyMatrix<2>; - auto [i_symbol, found] = symbol_table->find("R2x2_2d", position); - REQUIRE(found); - REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t); + SECTION("matrix 2d") + { + using R2x2 = TinyMatrix<2>; + auto [i_symbol, found] = symbol_table->find("R2x2_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); + FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table); - ASTNode::setStackDetails(false); + ASTNode::setStackDetails(false); - Array<R2x2> integrate_value_polygonal(mesh_2d->numberOfCells()); - IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, - *mesh_2d, integrate_value_polygonal); + Array<R2x2> integrate_value_polygonal(mesh_2d->numberOfCells()); + IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, + *mesh_2d, integrate_value_polygonal); - Array<R2x2> integrate_value_parabolic(parabolic_mesh->numberOfCells()); - IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, - *parabolic_mesh, integrate_value_parabolic); + Array<R2x2> integrate_value_parabolic(parabolic_mesh->numberOfCells()); + IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, + *parabolic_mesh, integrate_value_parabolic); - ASTNode::setStackDetails(true); + ASTNode::setStackDetails(true); - REQUIRE(matrix_error(integrate_value_parabolic, integrate_value_polygonal) < 1E-12); + REQUIRE(matrix_error(integrate_value_parabolic, integrate_value_polygonal) < 1E-12); + } } } } } } - } - SECTION("Gauss-Legendre quadrature") - { - auto quadrature_descriptor = GaussLegendreQuadratureDescriptor(12); - - SECTION("integrate on all cells") + SECTION("Gauss-Legendre quadrature") { - SECTION("2D") + auto quadrature_descriptor = GaussLegendreQuadratureDescriptor(12); + + SECTION("integrate on all cells") { - constexpr size_t Dimension = 2; + SECTION("2D") + { + constexpr size_t Dimension = 2; - std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes(); + std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes(); - for (const auto& named_mesh : mesh_list) { - SECTION(named_mesh.name()) - { - auto mesh_2d_v = named_mesh.mesh(); - auto mesh_2d = mesh_2d_v->get<Mesh<2>>(); + for (const auto& named_mesh : mesh_list) { + SECTION(named_mesh.name()) + { + auto mesh_2d_v = named_mesh.mesh(); + auto mesh_2d = mesh_2d_v->get<Mesh<2>>(); - PolynomialMeshBuilder pb{mesh_2d_v, 2}; + PolynomialMeshBuilder pb{mesh_2d_v, 2}; - auto parabolic_mesh = pb.mesh()->get<PolynomialMesh<2>>(); + auto parabolic_mesh = pb.mesh()->get<PolynomialMesh<2>>(); - std::string_view data = R"( + std::string_view data = R"( import math; let scalar_2d: R^2 -> R, x -> 2*x[0]*x[1]+3; let R3_2d: R^2 -> R^3, x -> [2*x[0]*x[0] + x[1]*x[1] - 1, x[0]-2*x[1], 3]; let R2x2_2d: R^2 -> R^2x2, x -> [[2*x[0]*x[0]*x[1]+3, x[0]-2*x[1]], [3, x[0]*x[1]]]; )"; - TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"}; + TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"}; - auto ast = ASTBuilder::build(input); + auto ast = ASTBuilder::build(input); - ASTModulesImporter{*ast}; - ASTNodeTypeCleaner<language::import_instruction>{*ast}; + ASTModulesImporter{*ast}; + ASTNodeTypeCleaner<language::import_instruction>{*ast}; - ASTSymbolTableBuilder{*ast}; - ASTNodeDataTypeBuilder{*ast}; + ASTSymbolTableBuilder{*ast}; + ASTNodeDataTypeBuilder{*ast}; - ASTNodeTypeCleaner<language::var_declaration>{*ast}; - ASTNodeTypeCleaner<language::fct_declaration>{*ast}; - ASTNodeExpressionBuilder{*ast}; + ASTNodeTypeCleaner<language::var_declaration>{*ast}; + ASTNodeTypeCleaner<language::fct_declaration>{*ast}; + ASTNodeExpressionBuilder{*ast}; - std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table; + std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table; - // ensure that variables are declared at this point - TAO_PEGTL_NAMESPACE::position position{data.size(), 1, 1, "fixture"}; + // ensure that variables are declared at this point + TAO_PEGTL_NAMESPACE::position position{data.size(), 1, 1, "fixture"}; - SECTION("scalar 2d") - { - auto [i_symbol, found] = symbol_table->find("scalar_2d", position); - REQUIRE(found); - REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t); + SECTION("scalar 2d") + { + auto [i_symbol, found] = symbol_table->find("scalar_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); + FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table); - ASTNode::setStackDetails(false); + ASTNode::setStackDetails(false); - Array<double> integrate_value_polygonal(mesh_2d->numberOfCells()); - IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, - *mesh_2d, integrate_value_polygonal); + Array<double> integrate_value_polygonal(mesh_2d->numberOfCells()); + IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, + *mesh_2d, integrate_value_polygonal); - Array<double> integrate_value_parabolic(parabolic_mesh->numberOfCells()); - IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, - *parabolic_mesh, integrate_value_parabolic); + Array<double> integrate_value_parabolic(parabolic_mesh->numberOfCells()); + IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, + *parabolic_mesh, + integrate_value_parabolic); - ASTNode::setStackDetails(true); + ASTNode::setStackDetails(true); - REQUIRE(scalar_error(integrate_value_parabolic, integrate_value_polygonal) < 1E-12); - } + REQUIRE(scalar_error(integrate_value_parabolic, integrate_value_polygonal) < 1E-12); + } - SECTION("vector 2d") - { - using R3 = TinyVector<3>; - auto [i_symbol, found] = symbol_table->find("R3_2d", position); - REQUIRE(found); - REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t); + SECTION("vector 2d") + { + using R3 = TinyVector<3>; + auto [i_symbol, found] = symbol_table->find("R3_2d", position); + REQUIRE(found); + REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t); - FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table); + FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table); - ASTNode::setStackDetails(false); + ASTNode::setStackDetails(false); - Array<R3> integrate_value_polygonal(mesh_2d->numberOfCells()); - IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, - *mesh_2d, integrate_value_polygonal); + Array<R3> integrate_value_polygonal(mesh_2d->numberOfCells()); + IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, + *mesh_2d, integrate_value_polygonal); - Array<R3> integrate_value_parabolic(parabolic_mesh->numberOfCells()); - IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, - *parabolic_mesh, integrate_value_parabolic); + Array<R3> integrate_value_parabolic(parabolic_mesh->numberOfCells()); + IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, + *parabolic_mesh, integrate_value_parabolic); - ASTNode::setStackDetails(true); + ASTNode::setStackDetails(true); - REQUIRE(vector_error(integrate_value_parabolic, integrate_value_polygonal) < 1E-12); - } + REQUIRE(vector_error(integrate_value_parabolic, integrate_value_polygonal) < 1E-12); + } - SECTION("matrix 2d") - { - using R2x2 = TinyMatrix<2>; - auto [i_symbol, found] = symbol_table->find("R2x2_2d", position); - REQUIRE(found); - REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t); + SECTION("matrix 2d") + { + using R2x2 = TinyMatrix<2>; + auto [i_symbol, found] = symbol_table->find("R2x2_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); + FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table); - ASTNode::setStackDetails(false); + ASTNode::setStackDetails(false); - Array<R2x2> integrate_value_polygonal(mesh_2d->numberOfCells()); - IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, - *mesh_2d, integrate_value_polygonal); + Array<R2x2> integrate_value_polygonal(mesh_2d->numberOfCells()); + IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, + *mesh_2d, integrate_value_polygonal); - Array<R2x2> integrate_value_parabolic(parabolic_mesh->numberOfCells()); - IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, - *parabolic_mesh, integrate_value_parabolic); + Array<R2x2> integrate_value_parabolic(parabolic_mesh->numberOfCells()); + IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, + *parabolic_mesh, integrate_value_parabolic); - ASTNode::setStackDetails(true); + ASTNode::setStackDetails(true); - REQUIRE(matrix_error(integrate_value_parabolic, integrate_value_polygonal) < 1E-12); + REQUIRE(matrix_error(integrate_value_parabolic, integrate_value_polygonal) < 1E-12); + } } } } } } - } - - SECTION("Gauss-Lobatto quadrature") - { - auto quadrature_descriptor = GaussLobattoQuadratureDescriptor(12); - SECTION("integrate on all cells") + SECTION("Gauss-Lobatto quadrature") { - SECTION("2D") + auto quadrature_descriptor = GaussLobattoQuadratureDescriptor(12); + + SECTION("integrate on all cells") { - constexpr size_t Dimension = 2; + SECTION("2D") + { + constexpr size_t Dimension = 2; - std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes(); + std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes(); - for (const auto& named_mesh : mesh_list) { - SECTION(named_mesh.name()) - { - auto mesh_2d_v = named_mesh.mesh(); - auto mesh_2d = mesh_2d_v->get<Mesh<2>>(); + for (const auto& named_mesh : mesh_list) { + SECTION(named_mesh.name()) + { + auto mesh_2d_v = named_mesh.mesh(); + auto mesh_2d = mesh_2d_v->get<Mesh<2>>(); - PolynomialMeshBuilder pb{mesh_2d_v, 2}; + PolynomialMeshBuilder pb{mesh_2d_v, 2}; - auto parabolic_mesh = pb.mesh()->get<PolynomialMesh<2>>(); + auto parabolic_mesh = pb.mesh()->get<PolynomialMesh<2>>(); - std::string_view data = R"( + std::string_view data = R"( import math; let scalar_2d: R^2 -> R, x -> 2*x[0]*x[1]+3; let R3_2d: R^2 -> R^3, x -> [2*x[0]*x[0] + x[1]*x[1] - 1, x[0]-2*x[1], 3]; let R2x2_2d: R^2 -> R^2x2, x -> [[2*x[0]*x[0]*x[1]+3, x[0]-2*x[1]], [3, x[0]*x[1]]]; )"; - TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"}; + TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"}; - auto ast = ASTBuilder::build(input); + auto ast = ASTBuilder::build(input); - ASTModulesImporter{*ast}; - ASTNodeTypeCleaner<language::import_instruction>{*ast}; + ASTModulesImporter{*ast}; + ASTNodeTypeCleaner<language::import_instruction>{*ast}; - ASTSymbolTableBuilder{*ast}; - ASTNodeDataTypeBuilder{*ast}; + ASTSymbolTableBuilder{*ast}; + ASTNodeDataTypeBuilder{*ast}; - ASTNodeTypeCleaner<language::var_declaration>{*ast}; - ASTNodeTypeCleaner<language::fct_declaration>{*ast}; - ASTNodeExpressionBuilder{*ast}; + ASTNodeTypeCleaner<language::var_declaration>{*ast}; + ASTNodeTypeCleaner<language::fct_declaration>{*ast}; + ASTNodeExpressionBuilder{*ast}; - std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table; + std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table; - // ensure that variables are declared at this point - TAO_PEGTL_NAMESPACE::position position{data.size(), 1, 1, "fixture"}; + // ensure that variables are declared at this point + TAO_PEGTL_NAMESPACE::position position{data.size(), 1, 1, "fixture"}; - SECTION("scalar 2d") - { - auto [i_symbol, found] = symbol_table->find("scalar_2d", position); - REQUIRE(found); - REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t); + SECTION("scalar 2d") + { + auto [i_symbol, found] = symbol_table->find("scalar_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); + + ASTNode::setStackDetails(false); + + Array<double> integrate_value_polygonal(mesh_2d->numberOfCells()); + IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, + *mesh_2d, integrate_value_polygonal); + + Array<double> integrate_value_parabolic(parabolic_mesh->numberOfCells()); + IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, + *parabolic_mesh, + integrate_value_parabolic); + + ASTNode::setStackDetails(true); - FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table); + REQUIRE(scalar_error(integrate_value_parabolic, integrate_value_polygonal) < 1E-12); + } - ASTNode::setStackDetails(false); + SECTION("vector 2d") + { + using R3 = TinyVector<3>; + auto [i_symbol, found] = symbol_table->find("R3_2d", position); + REQUIRE(found); + REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t); + + FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table); + + ASTNode::setStackDetails(false); + + Array<R3> integrate_value_polygonal(mesh_2d->numberOfCells()); + IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, + *mesh_2d, integrate_value_polygonal); + + Array<R3> integrate_value_parabolic(parabolic_mesh->numberOfCells()); + IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, + *parabolic_mesh, integrate_value_parabolic); - Array<double> integrate_value_polygonal(mesh_2d->numberOfCells()); - IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, + ASTNode::setStackDetails(true); + + REQUIRE(vector_error(integrate_value_parabolic, integrate_value_polygonal) < 1E-12); + } + + SECTION("matrix 2d") + { + using R2x2 = TinyMatrix<2>; + auto [i_symbol, found] = symbol_table->find("R2x2_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); + + ASTNode::setStackDetails(false); + + Array<R2x2> integrate_value_polygonal(mesh_2d->numberOfCells()); + IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, *mesh_2d, integrate_value_polygonal); - Array<double> integrate_value_parabolic(parabolic_mesh->numberOfCells()); - IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, + Array<R2x2> integrate_value_parabolic(parabolic_mesh->numberOfCells()); + IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, *parabolic_mesh, integrate_value_parabolic); - ASTNode::setStackDetails(true); + ASTNode::setStackDetails(true); - REQUIRE(scalar_error(integrate_value_parabolic, integrate_value_polygonal) < 1E-12); + REQUIRE(matrix_error(integrate_value_parabolic, integrate_value_polygonal) < 1E-12); + } } + } + } + } + } + } - SECTION("vector 2d") - { - using R3 = TinyVector<3>; - auto [i_symbol, found] = symbol_table->find("R3_2d", position); - REQUIRE(found); - REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t); + SECTION("parabolic edges") + { + constexpr size_t Dimension = 2; - FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table); + using R2 = TinyVector<Dimension>; - ASTNode::setStackDetails(false); + auto mesh_transform = [](const R2& x) -> R2 { + return R2{x[0] * (1 + 0.2 * std::sin(x[0] + x[1])), x[1] * (1 + 0.3 * std::sin(x[0] + x[1]))}; + }; - Array<R3> integrate_value_polygonal(mesh_2d->numberOfCells()); - IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, - *mesh_2d, integrate_value_polygonal); + std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes(); - Array<R3> integrate_value_parabolic(parabolic_mesh->numberOfCells()); - IntegrateOnCells<R3(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, - *parabolic_mesh, integrate_value_parabolic); + for (const auto& named_mesh : mesh_list) { + SECTION(named_mesh.name()) + { + auto mesh_2d_v = named_mesh.mesh(); + auto mesh_2d = mesh_2d_v->get<Mesh<2>>(); - ASTNode::setStackDetails(true); + PolynomialMeshBuilder pb{mesh_2d_v, 2}; - REQUIRE(vector_error(integrate_value_parabolic, integrate_value_polygonal) < 1E-12); + auto parabolic_mesh = pb.mesh()->get<PolynomialMesh<2>>(); + + auto xr = copy(parabolic_mesh->xr()); + auto xl = copy(parabolic_mesh->xl()); + + parallel_for( + parabolic_mesh->numberOfNodes(), + PUGS_LAMBDA(const NodeId& node_id) { xr[node_id] = mesh_transform(xr[node_id]); }); + + parallel_for( + parabolic_mesh->numberOfFaces(), PUGS_LAMBDA(const FaceId& face_id) { + for (size_t i = 0; i < xl.sizeOfArrays(); ++i) { + xl[face_id][i] = mesh_transform(xl[face_id][i]); } + }); - SECTION("matrix 2d") - { - using R2x2 = TinyMatrix<2>; - auto [i_symbol, found] = symbol_table->find("R2x2_2d", position); - REQUIRE(found); - REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t); + parabolic_mesh = std::make_shared<PolynomialMesh<2>>(parabolic_mesh->shared_connectivity(), xr, xl); - FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table); + auto Q = [](const R2& X) -> R2 { + const double& x = X[0]; + const double& y = X[1]; + return R2{3 + 2 * x - 3 * y + 5 * x * y, // + -4 - 3 * x + y * y + 2 * x * x * y}; + }; - ASTNode::setStackDetails(false); + FaceValue<double> face_integral(parabolic_mesh->connectivity()); + const auto face_to_node_matrix = parabolic_mesh->connectivity().faceToNodeMatrix(); - Array<R2x2> integrate_value_polygonal(mesh_2d->numberOfCells()); - IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, - *mesh_2d, integrate_value_polygonal); + QuadratureFormula<1> qf = QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(7)); - Array<R2x2> integrate_value_parabolic(parabolic_mesh->numberOfCells()); - IntegrateOnCells<R2x2(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, - *parabolic_mesh, integrate_value_parabolic); + parallel_for( + parabolic_mesh->numberOfFaces(), PUGS_LAMBDA(const FaceId& face_id) { + const NodeId node0 = face_to_node_matrix[face_id][0]; + const NodeId node1 = face_to_node_matrix[face_id][1]; - ASTNode::setStackDetails(true); + LineParabolicTransformation<Dimension> T{xr[node0], xl[face_id][0], xr[node1]}; - REQUIRE(matrix_error(integrate_value_parabolic, integrate_value_polygonal) < 1E-12); + double sum = 0; + for (size_t i = 0; i < qf.numberOfPoints(); ++i) { + const auto& x_hat = qf.point(i); + + const R2 Q_i = Q(T(x_hat)); + const R2 T_prime = T.velocity(x_hat); + + const R2 N_hat{T_prime[1], -T_prime[0]}; + + const double& omega_i = qf.weight(i); + + sum += omega_i * dot(Q_i, N_hat); } - } + + face_integral[face_id] = sum; + }); + + auto cell_face_is_reversed = parabolic_mesh->connectivity().cellFaceIsReversed(); + const auto cell_to_face_matrix = parabolic_mesh->connectivity().cellToFaceMatrix(); + + CellValue<double> integrate_ref{parabolic_mesh->connectivity()}; + + parallel_for( + parabolic_mesh->numberOfCells(), PUGS_LAMBDA(const CellId& cell_id) { + const auto cell_faces = cell_to_face_matrix[cell_id]; + const auto face_is_reversed = cell_face_is_reversed[cell_id]; + + double sum = 0; + for (size_t i_face = 0; i_face < cell_faces.size(); ++i_face) { + const FaceId face_id = cell_faces[i_face]; + const double sign = (face_is_reversed[i_face]) ? -1 : 1; + + sum += sign * face_integral[face_id]; + } + + integrate_ref[cell_id] = sum; + }); + + { + auto quadrature_descriptor = GaussQuadratureDescriptor(12); + + std::string_view data = R"( +import math; +let scalar_2d: R^2 -> R, x -> 2 + 7 * x[1] + 2 * x[0] * x[0]; +)"; + + TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"}; + + auto ast = ASTBuilder::build(input); + + ASTModulesImporter{*ast}; + ASTNodeTypeCleaner<language::import_instruction>{*ast}; + + ASTSymbolTableBuilder{*ast}; + ASTNodeDataTypeBuilder{*ast}; + + ASTNodeTypeCleaner<language::var_declaration>{*ast}; + ASTNodeTypeCleaner<language::fct_declaration>{*ast}; + ASTNodeExpressionBuilder{*ast}; + + std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table; + + // ensure that variables are declared at this point + TAO_PEGTL_NAMESPACE::position position{data.size(), 1, 1, "fixture"}; + + auto [i_symbol, found] = symbol_table->find("scalar_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> integrate_value(parabolic_mesh->numberOfCells()); + + ASTNode::setStackDetails(false); + + IntegrateOnCells<double(TinyVector<Dimension>)>::integrateTo(function_symbol_id, quadrature_descriptor, + *parabolic_mesh, integrate_value); + + ASTNode::setStackDetails(true); + + REQUIRE(scalar_error(integrate_ref.arrayView(), integrate_value) < 1E-12); } } } -- GitLab