diff --git a/src/scheme/DiscreteFunctionIntegrator.cpp b/src/scheme/DiscreteFunctionIntegrator.cpp
index fe0757cc0ef74e52b887e751fa678e7be6bdb186..c6a950a9ddc7f088996425da48a0affdcf61d233 100644
--- a/src/scheme/DiscreteFunctionIntegrator.cpp
+++ b/src/scheme/DiscreteFunctionIntegrator.cpp
@@ -73,26 +73,12 @@ DiscreteFunctionIntegrator::_integrateGlobally() const
using MeshType = Mesh<Connectivity<Dimension>>;
std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(m_mesh);
- if constexpr (std::is_same_v<DataType, ValueType>) {
- return std::make_shared<
- DiscreteFunctionP0<Dimension, ValueType>>(mesh,
- IntegrateCellValue<DataType(TinyVector<Dimension>)>::template integrate<
- MeshType>(m_function_id, *m_quadrature_descriptor, *mesh));
- } else {
- static_assert(std::is_convertible_v<DataType, ValueType>);
-
- CellValue<DataType> cell_data =
- IntegrateCellValue<DataType(TinyVector<Dimension>)>::template integrate<MeshType>(m_function_id,
- *m_quadrature_descriptor,
- *mesh);
-
- CellValue<ValueType> cell_value{mesh->connectivity()};
-
- parallel_for(
- mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { cell_value[cell_id] = cell_data[cell_id]; });
+ static_assert(std::is_convertible_v<DataType, ValueType>);
- return std::make_shared<DiscreteFunctionP0<Dimension, ValueType>>(mesh, cell_value);
- }
+ return std::make_shared<
+ DiscreteFunctionP0<Dimension, ValueType>>(mesh,
+ IntegrateCellValue<ValueType(TinyVector<Dimension>)>::template integrate<
+ MeshType>(m_function_id, *m_quadrature_descriptor, *mesh));
}
template <size_t Dimension, typename DataType, typename ValueType>
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 78997684c07a6ffb5457207c6ab4ffaf9d1992f5..2fd2c19e0bc610f8d82cc355ea1965eb486c7b3b 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -157,6 +157,7 @@ add_executable (unit_tests
add_executable (mpi_unit_tests
mpi_test_main.cpp
test_Connectivity.cpp
+ test_DiscreteFunctionIntegrator.cpp
test_DiscreteFunctionInterpoler.cpp
test_DiscreteFunctionInterpolerByZone.cpp
test_DiscreteFunctionP0.cpp
diff --git a/tests/test_DiscreteFunctionIntegrator.cpp b/tests/test_DiscreteFunctionIntegrator.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..465288f67d4efcd9a7a78eff53ab9649d381e572
--- /dev/null
+++ b/tests/test_DiscreteFunctionIntegrator.cpp
@@ -0,0 +1,781 @@
+#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/MeshData.hpp>
+#include <mesh/MeshDataManager.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("DiscreteFunctionIntegrator", "[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);
+
+ std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
+
+ for (auto named_mesh : mesh_list) {
+ SECTION(named_mesh.name())
+ {
+ auto mesh_1d = named_mesh.mesh();
+
+ std::string_view data = R"(
+import math;
+let B_scalar_non_linear_1d: R^1 -> B, x -> (exp(2 * x[0]) + 3 > 4);
+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 =
+ IntegrateCellValue<double(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_1d);
+
+ DiscreteFunctionIntegrator integrator(mesh_1d, 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);
+
+ CellValue<double> cell_value =
+ IntegrateCellValue<double(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_1d);
+
+ DiscreteFunctionIntegrator integrator(mesh_1d, 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);
+
+ CellValue<double> cell_value =
+ IntegrateCellValue<double(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_1d);
+
+ DiscreteFunctionIntegrator integrator(mesh_1d, 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);
+
+ CellValue<double> cell_value =
+ IntegrateCellValue<double(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_1d);
+
+ DiscreteFunctionIntegrator integrator(mesh_1d, 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);
+
+ CellValue<DataType> cell_value =
+ IntegrateCellValue<DataType(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor,
+ *mesh_1d);
+
+ DiscreteFunctionIntegrator integrator(mesh_1d, 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);
+
+ CellValue<DataType> cell_value =
+ IntegrateCellValue<DataType(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor,
+ *mesh_1d);
+
+ DiscreteFunctionIntegrator integrator(mesh_1d, 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);
+
+ CellValue<DataType> cell_value =
+ IntegrateCellValue<DataType(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor,
+ *mesh_1d);
+
+ DiscreteFunctionIntegrator integrator(mesh_1d, 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);
+
+ CellValue<DataType> cell_value =
+ IntegrateCellValue<DataType(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor,
+ *mesh_1d);
+
+ DiscreteFunctionIntegrator integrator(mesh_1d, 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);
+
+ CellValue<DataType> cell_value =
+ IntegrateCellValue<DataType(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor,
+ *mesh_1d);
+
+ DiscreteFunctionIntegrator integrator(mesh_1d, 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);
+
+ CellValue<DataType> cell_value =
+ IntegrateCellValue<DataType(TinyVector<1>)>::integrate(function_symbol_id, *quadrature_descriptor,
+ *mesh_1d);
+
+ DiscreteFunctionIntegrator integrator(mesh_1d, 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);
+
+ std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+ for (auto named_mesh : mesh_list) {
+ SECTION(named_mesh.name())
+ {
+ auto mesh_2d = named_mesh.mesh();
+
+ 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 =
+ IntegrateCellValue<double(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_2d);
+
+ DiscreteFunctionIntegrator integrator(mesh_2d, 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);
+
+ CellValue<double> cell_value =
+ IntegrateCellValue<double(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_2d);
+
+ DiscreteFunctionIntegrator integrator(mesh_2d, 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);
+
+ CellValue<double> cell_value =
+ IntegrateCellValue<double(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_2d);
+
+ DiscreteFunctionIntegrator integrator(mesh_2d, 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);
+
+ CellValue<double> cell_value =
+ IntegrateCellValue<double(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_2d);
+
+ DiscreteFunctionIntegrator integrator(mesh_2d, 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);
+
+ CellValue<DataType> cell_value =
+ IntegrateCellValue<DataType(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor,
+ *mesh_2d);
+
+ DiscreteFunctionIntegrator integrator(mesh_2d, 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);
+
+ CellValue<DataType> cell_value =
+ IntegrateCellValue<DataType(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor,
+ *mesh_2d);
+
+ DiscreteFunctionIntegrator integrator(mesh_2d, 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);
+
+ CellValue<DataType> cell_value =
+ IntegrateCellValue<DataType(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor,
+ *mesh_2d);
+
+ DiscreteFunctionIntegrator integrator(mesh_2d, 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);
+
+ CellValue<DataType> cell_value =
+ IntegrateCellValue<DataType(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor,
+ *mesh_2d);
+
+ DiscreteFunctionIntegrator integrator(mesh_2d, 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);
+
+ CellValue<DataType> cell_value =
+ IntegrateCellValue<DataType(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor,
+ *mesh_2d);
+
+ DiscreteFunctionIntegrator integrator(mesh_2d, 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);
+
+ CellValue<DataType> cell_value =
+ IntegrateCellValue<DataType(TinyVector<2>)>::integrate(function_symbol_id, *quadrature_descriptor,
+ *mesh_2d);
+
+ DiscreteFunctionIntegrator integrator(mesh_2d, 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);
+
+ std::array mesh_list = MeshDataBaseForTests::get().all3DMeshes();
+
+ for (auto named_mesh : mesh_list) {
+ SECTION(named_mesh.name())
+ {
+ auto mesh_3d = named_mesh.mesh();
+
+ 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 =
+ IntegrateCellValue<double(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_3d);
+
+ DiscreteFunctionIntegrator integrator(mesh_3d, 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);
+
+ CellValue<double> cell_value =
+ IntegrateCellValue<double(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_3d);
+
+ DiscreteFunctionIntegrator integrator(mesh_3d, 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);
+
+ CellValue<double> cell_value =
+ IntegrateCellValue<double(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_3d);
+
+ DiscreteFunctionIntegrator integrator(mesh_3d, 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);
+
+ CellValue<double> cell_value =
+ IntegrateCellValue<double(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor, *mesh_3d);
+
+ DiscreteFunctionIntegrator integrator(mesh_3d, 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);
+
+ CellValue<DataType> cell_value =
+ IntegrateCellValue<DataType(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor,
+ *mesh_3d);
+
+ DiscreteFunctionIntegrator integrator(mesh_3d, 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);
+
+ CellValue<DataType> cell_value =
+ IntegrateCellValue<DataType(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor,
+ *mesh_3d);
+
+ DiscreteFunctionIntegrator integrator(mesh_3d, 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);
+
+ CellValue<DataType> cell_value =
+ IntegrateCellValue<DataType(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor,
+ *mesh_3d);
+
+ DiscreteFunctionIntegrator integrator(mesh_3d, 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);
+
+ CellValue<DataType> cell_value =
+ IntegrateCellValue<DataType(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor,
+ *mesh_3d);
+
+ DiscreteFunctionIntegrator integrator(mesh_3d, 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);
+
+ CellValue<DataType> cell_value =
+ IntegrateCellValue<DataType(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor,
+ *mesh_3d);
+
+ DiscreteFunctionIntegrator integrator(mesh_3d, 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);
+
+ CellValue<DataType> cell_value =
+ IntegrateCellValue<DataType(TinyVector<3>)>::integrate(function_symbol_id, *quadrature_descriptor,
+ *mesh_3d);
+
+ DiscreteFunctionIntegrator integrator(mesh_3d, 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)));
+ }
+ }
+ }
+ }
+}