Select Git revision
ScalarDiamondScheme.hpp
test_IntegrateCellValue.cpp 16.43 KiB
#include <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_all.hpp>
#include <language/ast/ASTBuilder.hpp>
#include <language/ast/ASTModulesImporter.hpp>
#include <language/ast/ASTNodeDataTypeBuilder.hpp>
#include <language/ast/ASTNodeExpressionBuilder.hpp>
#include <language/ast/ASTNodeFunctionEvaluationExpressionBuilder.hpp>
#include <language/ast/ASTNodeFunctionExpressionBuilder.hpp>
#include <language/ast/ASTNodeTypeCleaner.hpp>
#include <language/ast/ASTSymbolTableBuilder.hpp>
#include <language/utils/PugsFunctionAdapter.hpp>
#include <language/utils/SymbolTable.hpp>
#include <MeshDataBaseForTests.hpp>
#include <mesh/Connectivity.hpp>
#include <mesh/DualMeshManager.hpp>
#include <mesh/Mesh.hpp>
#include <scheme/CellIntegrator.hpp>
#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
#include <analysis/GaussLobattoQuadratureDescriptor.hpp>
#include <analysis/GaussQuadratureDescriptor.hpp>
#include <language/utils/IntegrateCellValue.hpp>
// clazy:excludeall=non-pod-global-static
TEST_CASE("IntegrateCellValue", "[language]")
{
SECTION("integrate on all cells")
{
auto same_item_integral = [](auto f, auto g) -> bool {
using ItemIdType = typename decltype(f)::index_type;
for (ItemIdType item_id = 0; item_id < f.numberOfItems(); ++item_id) {
if (f[item_id] != g[item_id]) {
return false;
}
}
return true;
};
SECTION("1D")
{
constexpr size_t Dimension = 1;
auto quadrature_descriptor = GaussQuadratureDescriptor(3);
std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
for (const auto& named_mesh : mesh_list) {
SECTION(named_mesh.name())
{
auto mesh_1d = named_mesh.mesh();
std::string_view data = R"(
import math;
let R2x2_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]]];
)";
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"};
using R2x2 = TinyMatrix<2>;
auto [i_symbol, found] = symbol_table->find("R2x2_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<R2x2> cell_integral{mesh_1d->connectivity()};
auto f = [](const TinyVector<Dimension>& x) -> R2x2 {
return R2x2{2 * exp(x[0]) * sin(x[0]) + 3, sin(x[0] - 2 * x[0]), 3, x[0] * x[0]};
};
CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_1d, cell_integral);
CellValue<R2x2> integrate_value =
IntegrateCellValue<R2x2(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
*mesh_1d);
REQUIRE(same_item_integral(cell_integral, integrate_value));
}
}
}
SECTION("2D")
{
constexpr size_t Dimension = 2;
auto quadrature_descriptor = GaussLobattoQuadratureDescriptor(3);
std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
for (const auto& named_mesh : mesh_list) {
SECTION(named_mesh.name())
{
auto mesh_2d = named_mesh.mesh();
std::string_view data = R"(
import math;
let R3_2d: R^2 -> R^3, x -> [2*exp(x[0])*sin(x[1])+3, x[0]-2*x[1], 3];
)";
TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
auto ast = ASTBuilder::build(input);
ASTModulesImporter{*ast};
ASTNodeTypeCleaner<language::import_instruction>{*ast};
ASTSymbolTableBuilder{*ast};
ASTNodeDataTypeBuilder{*ast};
ASTNodeTypeCleaner<language::var_declaration>{*ast};
ASTNodeTypeCleaner<language::fct_declaration>{*ast};
ASTNodeExpressionBuilder{*ast};
std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
// ensure that variables are declared at this point
TAO_PEGTL_NAMESPACE::position position{data.size(), 1, 1, "fixture"};
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);
CellValue<R3> cell_integral{mesh_2d->connectivity()};
auto f = [](const TinyVector<Dimension>& x) -> R3 {
return R3{2 * exp(x[0]) * sin(x[1]) + 3, x[0] - 2 * x[1], 3};
};
CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_2d, cell_integral);
CellValue<R3> integrate_value =
IntegrateCellValue<R3(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
*mesh_2d);
REQUIRE(same_item_integral(cell_integral, integrate_value));
}
}
}
SECTION("3D")
{
constexpr size_t Dimension = 3;
auto quadrature_descriptor = GaussLegendreQuadratureDescriptor(3);
using NamedMesh = MeshDataBaseForTests::NamedMesh<Dimension>;
std::vector<NamedMesh> mesh_list = [] {
std::vector<NamedMesh> extended_mesh_list;
std::array mesh_array = MeshDataBaseForTests::get().all3DMeshes();
for (size_t i = 0; i < mesh_array.size(); ++i) {
extended_mesh_list.push_back(MeshDataBaseForTests::get().all3DMeshes()[i]);
}
extended_mesh_list.push_back(NamedMesh("diamond dual", DualMeshManager::instance().getDiamondDualMesh(
*MeshDataBaseForTests::get().hybrid3DMesh())));
return extended_mesh_list;
}();
for (const auto& named_mesh : mesh_list) {
SECTION(named_mesh.name())
{
auto mesh_3d = named_mesh.mesh();
std::string_view data = R"(
import math;
let scalar_3d: R^3 -> R, x -> 2 * exp(x[0]) * sin(x[1]) * x[2] + 3;
)";
TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
auto ast = ASTBuilder::build(input);
ASTModulesImporter{*ast};
ASTNodeTypeCleaner<language::import_instruction>{*ast};
ASTSymbolTableBuilder{*ast};
ASTNodeDataTypeBuilder{*ast};
ASTNodeTypeCleaner<language::var_declaration>{*ast};
ASTNodeTypeCleaner<language::fct_declaration>{*ast};
ASTNodeExpressionBuilder{*ast};
std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
// 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_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_integral{mesh_3d->connectivity()};
auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * exp(x[0]) * sin(x[1]) * x[2] + 3; };
CellIntegrator::integrateTo(f, quadrature_descriptor, *mesh_3d, cell_integral);
CellValue<double> integrate_value =
IntegrateCellValue<double(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
*mesh_3d);
REQUIRE(same_item_integral(cell_integral, integrate_value));
}
}
}
}
SECTION("integrate on cell list")
{
auto same_item_integral = [](auto f, auto g) -> bool {
using ItemIdType = typename decltype(g)::index_type;
for (ItemIdType item_id = 0; item_id < f.size(); ++item_id) {
if (f[item_id] != g[item_id]) {
return false;
}
}
return true;
};
SECTION("1D")
{
constexpr size_t Dimension = 1;
auto quadrature_descriptor = GaussLegendreQuadratureDescriptor(3);
std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
for (const auto& named_mesh : mesh_list) {
SECTION(named_mesh.name())
{
auto mesh_1d = named_mesh.mesh();
Array<CellId> cell_list{mesh_1d->numberOfCells() / 2 + mesh_1d->numberOfCells() % 2};
{
size_t k = 0;
for (CellId cell_id = 0; cell_id < mesh_1d->numberOfCells(); ++(++cell_id), ++k) {
cell_list[k] = cell_id;
}
REQUIRE(k == cell_list.size());
}
std::string_view data = R"(
import math;
let scalar_1d: R^1 -> R, x -> 2 * exp(x[0]) + 3;
)";
TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
auto ast = ASTBuilder::build(input);
ASTModulesImporter{*ast};
ASTNodeTypeCleaner<language::import_instruction>{*ast};
ASTSymbolTableBuilder{*ast};
ASTNodeDataTypeBuilder{*ast};
ASTNodeTypeCleaner<language::var_declaration>{*ast};
ASTNodeTypeCleaner<language::fct_declaration>{*ast};
ASTNodeExpressionBuilder{*ast};
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_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);
auto f = [](const TinyVector<Dimension>& x) -> double { return 2 * std::exp(x[0]) + 3; };
Array<const double> cell_integral = CellIntegrator::integrate(f, quadrature_descriptor, *mesh_1d, cell_list);
Array<const double> integrate_value =
IntegrateCellValue<double(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
*mesh_1d, cell_list);
REQUIRE(same_item_integral(cell_integral, integrate_value));
}
}
}
SECTION("2D")
{
constexpr size_t Dimension = 2;
auto quadrature_descriptor = GaussLegendreQuadratureDescriptor(3);
std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
for (const auto& named_mesh : mesh_list) {
SECTION(named_mesh.name())
{
auto mesh_2d = named_mesh.mesh();
Array<CellId> cell_list{mesh_2d->numberOfCells() / 2 + mesh_2d->numberOfCells() % 2};
{
size_t k = 0;
for (CellId cell_id = 0; cell_id < mesh_2d->numberOfCells(); ++(++cell_id), ++k) {
cell_list[k] = cell_id;
}
REQUIRE(k == cell_list.size());
}
std::string_view data = R"(
import math;
let R3_2d: R^2 -> R^3, x -> [2*exp(x[0])*sin(x[1])+3, x[0]-2*x[1], 3];
)";
TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
auto ast = ASTBuilder::build(input);
ASTModulesImporter{*ast};
ASTNodeTypeCleaner<language::import_instruction>{*ast};
ASTSymbolTableBuilder{*ast};
ASTNodeDataTypeBuilder{*ast};
ASTNodeTypeCleaner<language::var_declaration>{*ast};
ASTNodeTypeCleaner<language::fct_declaration>{*ast};
ASTNodeExpressionBuilder{*ast};
std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
// ensure that variables are declared at this point
TAO_PEGTL_NAMESPACE::position position{data.size(), 1, 1, "fixture"};
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);
auto f = [](const TinyVector<Dimension>& x) -> R3 {
return R3{2 * exp(x[0]) * sin(x[1]) + 3, x[0] - 2 * x[1], 3};
};
Array<const R3> cell_integral = CellIntegrator::integrate(f, quadrature_descriptor, *mesh_2d, cell_list);
Array<const R3> integrate_value =
IntegrateCellValue<R3(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
*mesh_2d, cell_list);
REQUIRE(same_item_integral(cell_integral, integrate_value));
}
}
}
SECTION("3D")
{
constexpr size_t Dimension = 3;
auto quadrature_descriptor = GaussQuadratureDescriptor(3);
using NamedMesh = MeshDataBaseForTests::NamedMesh<Dimension>;
std::vector<NamedMesh> mesh_list = [] {
std::vector<NamedMesh> extended_mesh_list;
std::array mesh_array = MeshDataBaseForTests::get().all3DMeshes();
for (size_t i = 0; i < mesh_array.size(); ++i) {
extended_mesh_list.push_back(MeshDataBaseForTests::get().all3DMeshes()[i]);
}
extended_mesh_list.push_back(NamedMesh("diamond dual", DualMeshManager::instance().getDiamondDualMesh(
*MeshDataBaseForTests::get().hybrid3DMesh())));
return extended_mesh_list;
}();
for (const auto& named_mesh : mesh_list) {
SECTION(named_mesh.name())
{
auto mesh_3d = named_mesh.mesh();
Array<CellId> cell_list{mesh_3d->numberOfCells() / 2 + mesh_3d->numberOfCells() % 2};
{
size_t k = 0;
for (CellId cell_id = 0; cell_id < mesh_3d->numberOfCells(); ++(++cell_id), ++k) {
cell_list[k] = cell_id;
}
REQUIRE(k == cell_list.size());
}
std::string_view data = R"(
import math;
let R2x2_3d: R^3 -> R^2x2, x -> [[2 * exp(x[0]) * sin(x[1]) + 3 * cos(x[2]), sin(x[0] - 2 * x[1] * x[2])], [3, x[0] * 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;
// ensure that variables are declared at this point
TAO_PEGTL_NAMESPACE::position position{data.size(), 1, 1, "fixture"};
using R2x2 = TinyMatrix<2>;
auto [i_symbol, found] = symbol_table->find("R2x2_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);
auto f = [](const TinyVector<Dimension>& x) -> R2x2 {
return R2x2{2 * exp(x[0]) * sin(x[1]) + 3 * cos(x[2]), sin(x[0] - 2 * x[1] * x[2]), 3, x[0] * x[1] * x[2]};
};
Array<const R2x2> cell_integral = CellIntegrator::integrate(f, quadrature_descriptor, *mesh_3d, cell_list);
Array<R2x2> integrate_value =
IntegrateCellValue<R2x2(TinyVector<Dimension>)>::integrate(function_symbol_id, quadrature_descriptor,
*mesh_3d, cell_list);
REQUIRE(same_item_integral(cell_integral, integrate_value));
}
}
}
}
}