From 39ecc278652d70b92848a4df7ee542b6139d5b78 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Wed, 30 Mar 2022 21:35:21 +0200
Subject: [PATCH] Add an optional zone(s) argument to interpolate functions
The values are set to 0 in cells that are not in the zone(s)
---
src/language/modules/SchemeModule.cpp | 50 +
src/scheme/DiscreteFunctionInterpoler.cpp | 87 +-
src/scheme/DiscreteFunctionInterpoler.hpp | 18 +
.../DiscreteFunctionVectorInterpoler.cpp | 84 +-
.../DiscreteFunctionVectorInterpoler.hpp | 19 +
tests/CMakeLists.txt | 2 +
.../test_DiscreteFunctionInterpolerByZone.cpp | 1074 +++++++++++++++++
...DiscreteFunctionVectorInterpolerByZone.cpp | 511 ++++++++
8 files changed, 1832 insertions(+), 13 deletions(-)
create mode 100644 tests/test_DiscreteFunctionInterpolerByZone.cpp
create mode 100644 tests/test_DiscreteFunctionVectorInterpolerByZone.cpp
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index c724e1674..e4b52dc22 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -11,6 +11,7 @@
#include <language/utils/TypeDescriptor.hpp>
#include <mesh/Connectivity.hpp>
#include <mesh/IBoundaryDescriptor.hpp>
+#include <mesh/IZoneDescriptor.hpp>
#include <mesh/Mesh.hpp>
#include <mesh/MeshData.hpp>
#include <mesh/MeshDataManager.hpp>
@@ -158,6 +159,55 @@ SchemeModule::SchemeModule()
));
+ this->_addBuiltinFunction("interpolate",
+ std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
+ const IDiscreteFunction>(std::shared_ptr<const IMesh>,
+ const std::vector<std::shared_ptr<const IZoneDescriptor>>&,
+ std::shared_ptr<const IDiscreteFunctionDescriptor>,
+ const std::vector<FunctionSymbolId>&)>>(
+ [](std::shared_ptr<const IMesh> mesh,
+ const std::vector<std::shared_ptr<const IZoneDescriptor>>& interpolation_zone_list,
+ std::shared_ptr<const IDiscreteFunctionDescriptor> discrete_function_descriptor,
+ const std::vector<FunctionSymbolId>& function_id_list)
+ -> std::shared_ptr<const IDiscreteFunction> {
+ return DiscreteFunctionVectorInterpoler{mesh, interpolation_zone_list,
+ discrete_function_descriptor, function_id_list}
+ .interpolate();
+ }
+
+ ));
+
+ this->_addBuiltinFunction("interpolate",
+ std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
+ const IDiscreteFunction>(std::shared_ptr<const IMesh>,
+ const std::vector<std::shared_ptr<const IZoneDescriptor>>&,
+ std::shared_ptr<const IDiscreteFunctionDescriptor>,
+ const FunctionSymbolId&)>>(
+ [](std::shared_ptr<const IMesh> mesh,
+ const std::vector<std::shared_ptr<const IZoneDescriptor>>& interpolation_zone_list,
+ std::shared_ptr<const IDiscreteFunctionDescriptor> discrete_function_descriptor,
+ const FunctionSymbolId& function_id) -> std::shared_ptr<const IDiscreteFunction> {
+ switch (discrete_function_descriptor->type()) {
+ case DiscreteFunctionType::P0: {
+ return DiscreteFunctionInterpoler{mesh, interpolation_zone_list,
+ discrete_function_descriptor, function_id}
+ .interpolate();
+ }
+ case DiscreteFunctionType::P0Vector: {
+ return DiscreteFunctionVectorInterpoler{mesh,
+ interpolation_zone_list,
+ discrete_function_descriptor,
+ {function_id}}
+ .interpolate();
+ }
+ default: {
+ throw NormalError("invalid function descriptor type");
+ }
+ }
+ }
+
+ ));
+
this->_addBuiltinFunction("randomizeMesh",
std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
const IMesh>(std::shared_ptr<const IMesh>,
diff --git a/src/scheme/DiscreteFunctionInterpoler.cpp b/src/scheme/DiscreteFunctionInterpoler.cpp
index c39b6ee9d..daf7646ee 100644
--- a/src/scheme/DiscreteFunctionInterpoler.cpp
+++ b/src/scheme/DiscreteFunctionInterpoler.cpp
@@ -1,20 +1,82 @@
#include <scheme/DiscreteFunctionInterpoler.hpp>
#include <language/utils/InterpolateItemValue.hpp>
+#include <mesh/MeshCellZone.hpp>
#include <scheme/DiscreteFunctionP0.hpp>
#include <utils/Exceptions.hpp>
template <size_t Dimension, typename DataType, typename ValueType>
std::shared_ptr<IDiscreteFunction>
-DiscreteFunctionInterpoler::_interpolate() const
+DiscreteFunctionInterpoler::_interpolateOnZoneList() const
+{
+ static_assert(std::is_convertible_v<DataType, ValueType>);
+ Assert(m_zone_list.size() > 0, "no zone list provided");
+
+ std::shared_ptr p_mesh = std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(m_mesh);
+ using MeshDataType = MeshData<Dimension>;
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*p_mesh);
+
+ CellValue<bool> is_in_zone{p_mesh->connectivity()};
+ is_in_zone.fill(false);
+
+ size_t number_of_cells = 0;
+ for (const auto& zone : m_zone_list) {
+ auto mesh_cell_zone = getMeshCellZone(*p_mesh, *zone);
+ const auto& cell_list = mesh_cell_zone.cellList();
+ for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
+ const CellId cell_id = cell_list[i_cell];
+ if (is_in_zone[cell_id]) {
+ std::ostringstream os;
+ os << "cell " << cell_id << " (number " << p_mesh->connectivity().cellNumber()[cell_id]
+ << ") is present multiple times in zone list";
+ throw NormalError(os.str());
+ }
+ ++number_of_cells;
+ is_in_zone[cell_id] = true;
+ }
+ }
+
+ Array<CellId> zone_cell_list{number_of_cells};
+ {
+ size_t i_cell = 0;
+ for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
+ if (is_in_zone[cell_id]) {
+ zone_cell_list[i_cell++] = cell_id;
+ }
+ }
+ }
+
+ Array<const DataType> array =
+ InterpolateItemValue<DataType(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(m_function_id,
+ mesh_data.xj(),
+ zone_cell_list);
+
+ CellValue<ValueType> cell_value{p_mesh->connectivity()};
+ if constexpr (is_tiny_vector_v<ValueType> or is_tiny_matrix_v<ValueType>) {
+ cell_value.fill(zero);
+ } else {
+ cell_value.fill(0);
+ }
+
+ parallel_for(
+ zone_cell_list.size(), PUGS_LAMBDA(const size_t i_cell) { cell_value[zone_cell_list[i_cell]] = array[i_cell]; });
+
+ return std::make_shared<DiscreteFunctionP0<Dimension, ValueType>>(p_mesh, cell_value);
+}
+
+template <size_t Dimension, typename DataType, typename ValueType>
+std::shared_ptr<IDiscreteFunction>
+DiscreteFunctionInterpoler::_interpolateGlobally() const
{
- std::shared_ptr mesh = std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(m_mesh);
+ Assert(m_zone_list.size() == 0, "invalid call when zones are defined");
+
+ std::shared_ptr p_mesh = std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(m_mesh);
using MeshDataType = MeshData<Dimension>;
- MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*p_mesh);
if constexpr (std::is_same_v<DataType, ValueType>) {
return std::make_shared<
- DiscreteFunctionP0<Dimension, ValueType>>(mesh,
+ DiscreteFunctionP0<Dimension, ValueType>>(p_mesh,
InterpolateItemValue<DataType(TinyVector<Dimension>)>::
template interpolate<ItemType::cell>(m_function_id, mesh_data.xj()));
} else {
@@ -24,12 +86,23 @@ DiscreteFunctionInterpoler::_interpolate() const
InterpolateItemValue<DataType(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(m_function_id,
mesh_data.xj());
- CellValue<ValueType> cell_value{mesh->connectivity()};
+ CellValue<ValueType> cell_value{p_mesh->connectivity()};
parallel_for(
- mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { cell_value[cell_id] = cell_data[cell_id]; });
+ p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { cell_value[cell_id] = cell_data[cell_id]; });
- return std::make_shared<DiscreteFunctionP0<Dimension, ValueType>>(mesh, cell_value);
+ return std::make_shared<DiscreteFunctionP0<Dimension, ValueType>>(p_mesh, cell_value);
+ }
+}
+
+template <size_t Dimension, typename DataType, typename ValueType>
+std::shared_ptr<IDiscreteFunction>
+DiscreteFunctionInterpoler::_interpolate() const
+{
+ if (m_zone_list.size() == 0) {
+ return this->_interpolateGlobally<Dimension, DataType, ValueType>();
+ } else {
+ return this->_interpolateOnZoneList<Dimension, DataType, ValueType>();
}
}
diff --git a/src/scheme/DiscreteFunctionInterpoler.hpp b/src/scheme/DiscreteFunctionInterpoler.hpp
index e5724e3fe..295c2fd9a 100644
--- a/src/scheme/DiscreteFunctionInterpoler.hpp
+++ b/src/scheme/DiscreteFunctionInterpoler.hpp
@@ -3,6 +3,7 @@
#include <language/utils/FunctionSymbolId.hpp>
#include <mesh/IMesh.hpp>
+#include <mesh/IZoneDescriptor.hpp>
#include <scheme/IDiscreteFunction.hpp>
#include <scheme/IDiscreteFunctionDescriptor.hpp>
@@ -12,9 +13,16 @@ class DiscreteFunctionInterpoler
{
private:
std::shared_ptr<const IMesh> m_mesh;
+ const std::vector<std::shared_ptr<const IZoneDescriptor>> m_zone_list;
std::shared_ptr<const IDiscreteFunctionDescriptor> m_discrete_function_descriptor;
const FunctionSymbolId m_function_id;
+ template <size_t Dimension, typename DataType, typename ValueType = DataType>
+ std::shared_ptr<IDiscreteFunction> _interpolateOnZoneList() const;
+
+ template <size_t Dimension, typename DataType, typename ValueType = DataType>
+ std::shared_ptr<IDiscreteFunction> _interpolateGlobally() const;
+
template <size_t Dimension, typename DataType, typename ValueType = DataType>
std::shared_ptr<IDiscreteFunction> _interpolate() const;
@@ -30,6 +38,16 @@ class DiscreteFunctionInterpoler
: m_mesh{mesh}, m_discrete_function_descriptor{discrete_function_descriptor}, m_function_id{function_id}
{}
+ DiscreteFunctionInterpoler(const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const IZoneDescriptor>>& zone_list,
+ const std::shared_ptr<const IDiscreteFunctionDescriptor>& discrete_function_descriptor,
+ const FunctionSymbolId& function_id)
+ : m_mesh{mesh},
+ m_zone_list{zone_list},
+ m_discrete_function_descriptor{discrete_function_descriptor},
+ m_function_id{function_id}
+ {}
+
DiscreteFunctionInterpoler(const DiscreteFunctionInterpoler&) = delete;
DiscreteFunctionInterpoler(DiscreteFunctionInterpoler&&) = delete;
diff --git a/src/scheme/DiscreteFunctionVectorInterpoler.cpp b/src/scheme/DiscreteFunctionVectorInterpoler.cpp
index 98d6513c9..128c735e4 100644
--- a/src/scheme/DiscreteFunctionVectorInterpoler.cpp
+++ b/src/scheme/DiscreteFunctionVectorInterpoler.cpp
@@ -1,22 +1,94 @@
#include <scheme/DiscreteFunctionVectorInterpoler.hpp>
#include <language/utils/InterpolateItemArray.hpp>
+#include <mesh/MeshCellZone.hpp>
#include <scheme/DiscreteFunctionP0Vector.hpp>
#include <utils/Exceptions.hpp>
template <size_t Dimension, typename DataType>
std::shared_ptr<IDiscreteFunction>
-DiscreteFunctionVectorInterpoler::_interpolate() const
+DiscreteFunctionVectorInterpoler::_interpolateOnZoneList() const
{
- std::shared_ptr mesh = std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(m_mesh);
+ Assert(m_zone_list.size() > 0, "no zone list provided");
+ std::shared_ptr p_mesh = std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(m_mesh);
using MeshDataType = MeshData<Dimension>;
- MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*p_mesh);
+
+ CellValue<bool> is_in_zone{p_mesh->connectivity()};
+ is_in_zone.fill(false);
+
+ size_t number_of_cells = 0;
+ for (const auto& zone : m_zone_list) {
+ auto mesh_cell_zone = getMeshCellZone(*p_mesh, *zone);
+ const auto& cell_list = mesh_cell_zone.cellList();
+ for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
+ const CellId cell_id = cell_list[i_cell];
+ if (is_in_zone[cell_id]) {
+ std::ostringstream os;
+ os << "cell " << cell_id << " (number " << p_mesh->connectivity().cellNumber()[cell_id]
+ << ") is present multiple times in zone list";
+ throw NormalError(os.str());
+ }
+ ++number_of_cells;
+ is_in_zone[cell_id] = true;
+ }
+ }
+
+ Array<CellId> zone_cell_list{number_of_cells};
+ {
+ size_t i_cell = 0;
+ for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
+ if (is_in_zone[cell_id]) {
+ zone_cell_list[i_cell++] = cell_id;
+ }
+ }
+ }
+
+ Table<const DataType> table =
+ InterpolateItemArray<DataType(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(m_function_id_list,
+ mesh_data.xj(),
+ zone_cell_list);
+
+ CellArray<DataType> cell_array{p_mesh->connectivity(), m_function_id_list.size()};
+ cell_array.fill(0);
+
+ parallel_for(
+ zone_cell_list.size(), PUGS_LAMBDA(const size_t i_cell) {
+ for (size_t i = 0; i < table.numberOfColumns(); ++i) {
+ cell_array[zone_cell_list[i_cell]][i] = table(i_cell, i);
+ }
+ });
+
+ return std::make_shared<DiscreteFunctionP0Vector<Dimension, DataType>>(p_mesh, cell_array);
+}
+
+template <size_t Dimension, typename DataType>
+std::shared_ptr<IDiscreteFunction>
+DiscreteFunctionVectorInterpoler::_interpolateGlobally() const
+{
+ Assert(m_zone_list.size() == 0, "invalid call when zones are defined");
+
+ std::shared_ptr p_mesh = std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(m_mesh);
+
+ using MeshDataType = MeshData<Dimension>;
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*p_mesh);
return std::make_shared<
- DiscreteFunctionP0Vector<Dimension, DataType>>(mesh, InterpolateItemArray<DataType(TinyVector<Dimension>)>::
- template interpolate<ItemType::cell>(m_function_id_list,
- mesh_data.xj()));
+ DiscreteFunctionP0Vector<Dimension, DataType>>(p_mesh, InterpolateItemArray<DataType(TinyVector<Dimension>)>::
+ template interpolate<ItemType::cell>(m_function_id_list,
+ mesh_data.xj()));
+}
+
+template <size_t Dimension, typename DataType>
+std::shared_ptr<IDiscreteFunction>
+DiscreteFunctionVectorInterpoler::_interpolate() const
+{
+ if (m_zone_list.size() == 0) {
+ return this->_interpolateGlobally<Dimension, DataType>();
+ } else {
+ return this->_interpolateOnZoneList<Dimension, DataType>();
+ }
}
template <size_t Dimension>
diff --git a/src/scheme/DiscreteFunctionVectorInterpoler.hpp b/src/scheme/DiscreteFunctionVectorInterpoler.hpp
index 093f290c4..4bd917264 100644
--- a/src/scheme/DiscreteFunctionVectorInterpoler.hpp
+++ b/src/scheme/DiscreteFunctionVectorInterpoler.hpp
@@ -3,6 +3,7 @@
#include <language/utils/FunctionSymbolId.hpp>
#include <mesh/IMesh.hpp>
+#include <mesh/IZoneDescriptor.hpp>
#include <scheme/IDiscreteFunction.hpp>
#include <scheme/IDiscreteFunctionDescriptor.hpp>
@@ -13,9 +14,16 @@ class DiscreteFunctionVectorInterpoler
{
private:
std::shared_ptr<const IMesh> m_mesh;
+ const std::vector<std::shared_ptr<const IZoneDescriptor>> m_zone_list;
std::shared_ptr<const IDiscreteFunctionDescriptor> m_discrete_function_descriptor;
const std::vector<FunctionSymbolId> m_function_id_list;
+ template <size_t Dimension, typename DataType>
+ std::shared_ptr<IDiscreteFunction> _interpolateOnZoneList() const;
+
+ template <size_t Dimension, typename DataType>
+ std::shared_ptr<IDiscreteFunction> _interpolateGlobally() const;
+
template <size_t Dimension, typename DataType>
std::shared_ptr<IDiscreteFunction> _interpolate() const;
@@ -32,6 +40,17 @@ class DiscreteFunctionVectorInterpoler
: m_mesh{mesh}, m_discrete_function_descriptor{discrete_function_descriptor}, m_function_id_list{function_id_list}
{}
+ DiscreteFunctionVectorInterpoler(
+ const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const IZoneDescriptor>>& zone_list,
+ const std::shared_ptr<const IDiscreteFunctionDescriptor>& discrete_function_descriptor,
+ const std::vector<FunctionSymbolId>& function_id_list)
+ : m_mesh{mesh},
+ m_zone_list{zone_list},
+ m_discrete_function_descriptor{discrete_function_descriptor},
+ m_function_id_list{function_id_list}
+ {}
+
DiscreteFunctionVectorInterpoler(const DiscreteFunctionVectorInterpoler&) = delete;
DiscreteFunctionVectorInterpoler(DiscreteFunctionVectorInterpoler&&) = delete;
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 37a87d841..78997684c 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -158,9 +158,11 @@ add_executable (mpi_unit_tests
mpi_test_main.cpp
test_Connectivity.cpp
test_DiscreteFunctionInterpoler.cpp
+ test_DiscreteFunctionInterpolerByZone.cpp
test_DiscreteFunctionP0.cpp
test_DiscreteFunctionP0Vector.cpp
test_DiscreteFunctionVectorInterpoler.cpp
+ test_DiscreteFunctionVectorInterpolerByZone.cpp
test_EmbeddedIDiscreteFunctionMathFunctions.cpp
test_EmbeddedIDiscreteFunctionOperators.cpp
test_InterpolateItemArray.cpp
diff --git a/tests/test_DiscreteFunctionInterpolerByZone.cpp b/tests/test_DiscreteFunctionInterpolerByZone.cpp
new file mode 100644
index 000000000..4ad5f464f
--- /dev/null
+++ b/tests/test_DiscreteFunctionInterpolerByZone.cpp
@@ -0,0 +1,1074 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <language/ast/ASTBuilder.hpp>
+#include <language/ast/ASTModulesImporter.hpp>
+#include <language/ast/ASTNodeDataTypeBuilder.hpp>
+#include <language/ast/ASTNodeExpressionBuilder.hpp>
+#include <language/ast/ASTNodeFunctionEvaluationExpressionBuilder.hpp>
+#include <language/ast/ASTNodeFunctionExpressionBuilder.hpp>
+#include <language/ast/ASTNodeTypeCleaner.hpp>
+#include <language/ast/ASTSymbolTableBuilder.hpp>
+#include <language/utils/PugsFunctionAdapter.hpp>
+#include <language/utils/SymbolTable.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshCellZone.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <mesh/NamedZoneDescriptor.hpp>
+
+#include <scheme/DiscreteFunctionDescriptorP0.hpp>
+#include <scheme/DiscreteFunctionInterpoler.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
+
+#include <pegtl/string_input.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("DiscreteFunctionInterpolerByZone", "[scheme]")
+{
+ auto same_cell_value = [](auto f, auto g) -> bool {
+ using ItemIdType = typename decltype(f)::index_type;
+ for (ItemIdType item_id = 0; item_id < f.numberOfItems(); ++item_id) {
+ if (f[item_id] != g[item_id]) {
+ return false;
+ }
+ }
+
+ return true;
+ };
+
+ SECTION("1D")
+ {
+ constexpr size_t Dimension = 1;
+
+ auto mesh_1d = MeshDataBaseForTests::get().unordered1DMesh();
+
+ std::vector<std::shared_ptr<const IZoneDescriptor>> zone_list;
+ zone_list.push_back(std::make_shared<NamedZoneDescriptor>("LEFT"));
+
+ auto mesh_cell_zone = getMeshCellZone(*mesh_1d, *zone_list[0]);
+ CellValue<bool> is_cell_in_zone{mesh_1d->connectivity()};
+ is_cell_in_zone.fill(false);
+ auto zone_cell_list = mesh_cell_zone.cellList();
+ for (size_t i_cell = 0; i_cell < zone_cell_list.size(); ++i_cell) {
+ is_cell_in_zone[zone_cell_list[i_cell]] = true;
+ }
+
+ auto xj = MeshDataManager::instance().getMeshData(*mesh_1d).xj();
+
+ std::string_view data = R"(
+import math;
+let B_scalar_non_linear_1d: R^1 -> B, x -> (exp(2 * x[0]) + 3 < 3.3);
+let N_scalar_non_linear_1d: R^1 -> N, x -> floor(3 * x[0] * x[0] + 2);
+let Z_scalar_non_linear_1d: R^1 -> Z, x -> floor(exp(2 * x[0]) - 1);
+let R_scalar_non_linear_1d: R^1 -> R, x -> 2 * exp(x[0]) + 3;
+let R1_non_linear_1d: R^1 -> R^1, x -> 2 * exp(x[0]);
+let R2_non_linear_1d: R^1 -> R^2, x -> (2 * exp(x[0]), -3*x[0]);
+let R3_non_linear_1d: R^1 -> R^3, x -> (2 * exp(x[0]) + 3, x[0] - 2, 3);
+let R1x1_non_linear_1d: R^1 -> R^1x1, x -> (2 * exp(x[0]) * sin(x[0]) + 3);
+let R2x2_non_linear_1d: R^1 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[0]) + 3, sin(x[0] - 2 * x[0]), 3, x[0] * x[0]);
+let R3x3_non_linear_1d: R^1 -> R^3x3, x -> (2 * exp(x[0]) * sin(x[0]) + 3, sin(x[0] - 2 * x[0]), 3, x[0] * x[0], -4*x[0], 2*x[0]+1, 3, -6*x[0], exp(x[0]));
+)";
+ TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+ auto ast = ASTBuilder::build(input);
+
+ ASTModulesImporter{*ast};
+ ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+ ASTSymbolTableBuilder{*ast};
+ ASTNodeDataTypeBuilder{*ast};
+
+ ASTNodeTypeCleaner<language::var_declaration>{*ast};
+ ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+ ASTNodeExpressionBuilder{*ast};
+
+ std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+ TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+ position.byte = data.size(); // ensure that variables are declared at this point
+
+ SECTION("B_scalar_non_linear_1d")
+ {
+ auto [i_symbol, found] = symbol_table->find("B_scalar_non_linear_1d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<double> cell_value{mesh_1d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = std::exp(2 * x[0]) + 3 < 3.3;
+ } else {
+ cell_value[cell_id] = false;
+ }
+ });
+
+ DiscreteFunctionInterpoler interpoler(mesh_1d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+ function_symbol_id);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ REQUIRE(
+ same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+ }
+
+ SECTION("N_scalar_non_linear_1d")
+ {
+ auto [i_symbol, found] = symbol_table->find("N_scalar_non_linear_1d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<double> cell_value{mesh_1d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = std::floor(3 * x[0] * x[0] + 2);
+ } else {
+ cell_value[cell_id] = 0;
+ }
+ });
+
+ DiscreteFunctionInterpoler interpoler(mesh_1d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+ function_symbol_id);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ REQUIRE(
+ same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+ }
+
+ SECTION("Z_scalar_non_linear_1d")
+ {
+ auto [i_symbol, found] = symbol_table->find("Z_scalar_non_linear_1d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<double> cell_value{mesh_1d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = std::floor(std::exp(2 * x[0]) - 1);
+ } else {
+ cell_value[cell_id] = 0;
+ }
+ });
+
+ DiscreteFunctionInterpoler interpoler(mesh_1d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+ function_symbol_id);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ REQUIRE(
+ same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+ }
+
+ SECTION("R_scalar_non_linear_1d")
+ {
+ auto [i_symbol, found] = symbol_table->find("R_scalar_non_linear_1d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<double> cell_value{mesh_1d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = 2 * std::exp(x[0]) + 3;
+ } else {
+ cell_value[cell_id] = 0;
+ }
+ });
+
+ DiscreteFunctionInterpoler interpoler(mesh_1d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+ function_symbol_id);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ REQUIRE(
+ same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+ }
+
+ SECTION("R1_non_linear_1d")
+ {
+ using DataType = TinyVector<1>;
+
+ auto [i_symbol, found] = symbol_table->find("R1_non_linear_1d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<DataType> cell_value{mesh_1d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = DataType{2 * std::exp(x[0])};
+ } else {
+ cell_value[cell_id] = zero;
+ }
+ });
+
+ DiscreteFunctionInterpoler interpoler(mesh_1d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+ function_symbol_id);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ REQUIRE(
+ same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+ }
+
+ SECTION("R2_non_linear_1d")
+ {
+ using DataType = TinyVector<2>;
+
+ auto [i_symbol, found] = symbol_table->find("R2_non_linear_1d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<DataType> cell_value{mesh_1d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = DataType{2 * std::exp(x[0]), -3 * x[0]};
+ } else {
+ cell_value[cell_id] = zero;
+ }
+ });
+
+ DiscreteFunctionInterpoler interpoler(mesh_1d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+ function_symbol_id);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ REQUIRE(
+ same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+ }
+
+ SECTION("R3_non_linear_1d")
+ {
+ using DataType = TinyVector<3>;
+
+ auto [i_symbol, found] = symbol_table->find("R3_non_linear_1d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<DataType> cell_value{mesh_1d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = DataType{2 * std::exp(x[0]) + 3, x[0] - 2, 3};
+ } else {
+ cell_value[cell_id] = zero;
+ }
+ });
+
+ DiscreteFunctionInterpoler interpoler(mesh_1d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+ function_symbol_id);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ REQUIRE(
+ same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+ }
+
+ SECTION("R1x1_non_linear_1d")
+ {
+ using DataType = TinyMatrix<1>;
+
+ auto [i_symbol, found] = symbol_table->find("R1x1_non_linear_1d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<DataType> cell_value{mesh_1d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = DataType{2 * std::exp(x[0]) * std::sin(x[0]) + 3};
+ } else {
+ cell_value[cell_id] = zero;
+ }
+ });
+
+ DiscreteFunctionInterpoler interpoler(mesh_1d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+ function_symbol_id);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ REQUIRE(
+ same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+ }
+
+ SECTION("R2x2_non_linear_1d")
+ {
+ using DataType = TinyMatrix<2>;
+
+ auto [i_symbol, found] = symbol_table->find("R2x2_non_linear_1d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<DataType> cell_value{mesh_1d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] =
+ DataType{2 * std::exp(x[0]) * std::sin(x[0]) + 3, std::sin(x[0] - 2 * x[0]), 3, x[0] * x[0]};
+ } else {
+ cell_value[cell_id] = zero;
+ }
+ });
+
+ DiscreteFunctionInterpoler interpoler(mesh_1d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+ function_symbol_id);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ REQUIRE(
+ same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+ }
+
+ SECTION("R3x3_non_linear_1d")
+ {
+ using DataType = TinyMatrix<3>;
+
+ auto [i_symbol, found] = symbol_table->find("R3x3_non_linear_1d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<DataType> cell_value{mesh_1d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = DataType{2 * exp(x[0]) * std::sin(x[0]) + 3,
+ std::sin(x[0] - 2 * x[0]),
+ 3,
+ x[0] * x[0],
+ -4 * x[0],
+ 2 * x[0] + 1,
+ 3,
+ -6 * x[0],
+ std::exp(x[0])};
+ } else {
+ cell_value[cell_id] = zero;
+ }
+ });
+
+ DiscreteFunctionInterpoler interpoler(mesh_1d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+ function_symbol_id);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ REQUIRE(
+ same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+ }
+ }
+
+ SECTION("2D")
+ {
+ constexpr size_t Dimension = 2;
+
+ auto mesh_2d = MeshDataBaseForTests::get().hybrid2DMesh();
+
+ std::vector<std::shared_ptr<const IZoneDescriptor>> zone_list;
+ zone_list.push_back(std::make_shared<NamedZoneDescriptor>("LEFT"));
+
+ auto mesh_cell_zone = getMeshCellZone(*mesh_2d, *zone_list[0]);
+ CellValue<bool> is_cell_in_zone{mesh_2d->connectivity()};
+ is_cell_in_zone.fill(false);
+ auto zone_cell_list = mesh_cell_zone.cellList();
+ for (size_t i_cell = 0; i_cell < zone_cell_list.size(); ++i_cell) {
+ is_cell_in_zone[zone_cell_list[i_cell]] = true;
+ }
+
+ auto xj = MeshDataManager::instance().getMeshData(*mesh_2d).xj();
+
+ std::string_view data = R"(
+import math;
+let B_scalar_non_linear_2d: R^2 -> B, x -> (exp(2 * x[0])< 2*x[1]);
+let N_scalar_non_linear_2d: R^2 -> N, x -> floor(3 * (x[0] + x[1]) * (x[0] + x[1]) + 2);
+let Z_scalar_non_linear_2d: R^2 -> Z, x -> floor(exp(2 * x[0]) - 3 * x[1]);
+let R_scalar_non_linear_2d: R^2 -> R, x -> 2 * exp(x[0]) + 3 * x[1];
+let R1_non_linear_2d: R^2 -> R^1, x -> 2 * exp(x[0]);
+let R2_non_linear_2d: R^2 -> R^2, x -> (2 * exp(x[0]), -3*x[1]);
+let R3_non_linear_2d: R^2 -> R^3, x -> (2 * exp(x[0]) + 3, x[1] - 2, 3);
+let R1x1_non_linear_2d: R^2 -> R^1x1, x -> (2 * exp(x[0]) * sin(x[1]) + 3);
+let R2x2_non_linear_2d: R^2 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3, sin(x[1] - 2 * x[0]), 3, x[1] * x[0]);
+let R3x3_non_linear_2d: R^2 -> R^3x3, x -> (2 * exp(x[0]) * sin(x[1]) + 3, sin(x[1] - 2 * x[0]), 3, x[1] * x[0], -4*x[1], 2*x[0]+1, 3, -6*x[0], exp(x[1]));
+)";
+ TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+ auto ast = ASTBuilder::build(input);
+
+ ASTModulesImporter{*ast};
+ ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+ ASTSymbolTableBuilder{*ast};
+ ASTNodeDataTypeBuilder{*ast};
+
+ ASTNodeTypeCleaner<language::var_declaration>{*ast};
+ ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+ ASTNodeExpressionBuilder{*ast};
+
+ std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+ TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+ position.byte = data.size(); // ensure that variables are declared at this point
+
+ SECTION("B_scalar_non_linear_2d")
+ {
+ auto [i_symbol, found] = symbol_table->find("B_scalar_non_linear_2d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<double> cell_value{mesh_2d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = std::exp(2 * x[0]) < 2 * x[1];
+ } else {
+ cell_value[cell_id] = false;
+ }
+ });
+
+ DiscreteFunctionInterpoler interpoler(mesh_2d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+ function_symbol_id);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ REQUIRE(
+ same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+ }
+
+ SECTION("N_scalar_non_linear_2d")
+ {
+ auto [i_symbol, found] = symbol_table->find("N_scalar_non_linear_2d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<double> cell_value{mesh_2d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = std::floor(3 * (x[0] + x[1]) * (x[0] + x[1]) + 2);
+ } else {
+ cell_value[cell_id] = 0;
+ }
+ });
+
+ DiscreteFunctionInterpoler interpoler(mesh_2d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+ function_symbol_id);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ REQUIRE(
+ same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+ }
+
+ SECTION("Z_scalar_non_linear_2d")
+ {
+ auto [i_symbol, found] = symbol_table->find("Z_scalar_non_linear_2d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<double> cell_value{mesh_2d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = std::floor(std::exp(2 * x[0]) - 3 * x[1]);
+ } else {
+ cell_value[cell_id] = 0;
+ }
+ });
+
+ DiscreteFunctionInterpoler interpoler(mesh_2d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+ function_symbol_id);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ REQUIRE(
+ same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+ }
+
+ SECTION("R_scalar_non_linear_2d")
+ {
+ auto [i_symbol, found] = symbol_table->find("R_scalar_non_linear_2d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<double> cell_value{mesh_2d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = 2 * std::exp(x[0]) + 3 * x[1];
+ } else {
+ cell_value[cell_id] = 0;
+ }
+ });
+
+ DiscreteFunctionInterpoler interpoler(mesh_2d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+ function_symbol_id);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ REQUIRE(
+ same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+ }
+
+ SECTION("R1_non_linear_2d")
+ {
+ using DataType = TinyVector<1>;
+
+ auto [i_symbol, found] = symbol_table->find("R1_non_linear_2d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<DataType> cell_value{mesh_2d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = DataType{2 * std::exp(x[0])};
+ } else {
+ cell_value[cell_id] = zero;
+ }
+ });
+
+ DiscreteFunctionInterpoler interpoler(mesh_2d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+ function_symbol_id);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ REQUIRE(
+ same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+ }
+
+ SECTION("R2_non_linear_2d")
+ {
+ using DataType = TinyVector<2>;
+
+ auto [i_symbol, found] = symbol_table->find("R2_non_linear_2d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<DataType> cell_value{mesh_2d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = DataType{2 * std::exp(x[0]), -3 * x[1]};
+ } else {
+ cell_value[cell_id] = zero;
+ }
+ });
+
+ DiscreteFunctionInterpoler interpoler(mesh_2d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+ function_symbol_id);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ REQUIRE(
+ same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+ }
+
+ SECTION("R3_non_linear_2d")
+ {
+ using DataType = TinyVector<3>;
+
+ auto [i_symbol, found] = symbol_table->find("R3_non_linear_2d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<DataType> cell_value{mesh_2d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = DataType{2 * std::exp(x[0]) + 3, x[1] - 2, 3};
+ } else {
+ cell_value[cell_id] = zero;
+ }
+ });
+
+ DiscreteFunctionInterpoler interpoler(mesh_2d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+ function_symbol_id);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ REQUIRE(
+ same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+ }
+
+ SECTION("R1x1_non_linear_2d")
+ {
+ using DataType = TinyMatrix<1>;
+
+ auto [i_symbol, found] = symbol_table->find("R1x1_non_linear_2d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<DataType> cell_value{mesh_2d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3};
+ } else {
+ cell_value[cell_id] = zero;
+ }
+ });
+
+ DiscreteFunctionInterpoler interpoler(mesh_2d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+ function_symbol_id);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ REQUIRE(
+ same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+ }
+
+ SECTION("R2x2_non_linear_2d")
+ {
+ using DataType = TinyMatrix<2>;
+
+ auto [i_symbol, found] = symbol_table->find("R2x2_non_linear_2d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<DataType> cell_value{mesh_2d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] =
+ DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3, std::sin(x[1] - 2 * x[0]), 3, x[1] * x[0]};
+ } else {
+ cell_value[cell_id] = zero;
+ }
+ });
+
+ DiscreteFunctionInterpoler interpoler(mesh_2d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+ function_symbol_id);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ REQUIRE(
+ same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+ }
+
+ SECTION("R3x3_non_linear_2d")
+ {
+ using DataType = TinyMatrix<3>;
+
+ auto [i_symbol, found] = symbol_table->find("R3x3_non_linear_2d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<DataType> cell_value{mesh_2d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+
+ cell_value[cell_id] = DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3,
+ std::sin(x[1] - 2 * x[0]),
+ 3,
+ x[1] * x[0],
+ -4 * x[1],
+ 2 * x[0] + 1,
+ 3,
+ -6 * x[0],
+ std::exp(x[1])};
+ } else {
+ cell_value[cell_id] = zero;
+ }
+ });
+
+ DiscreteFunctionInterpoler interpoler(mesh_2d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+ function_symbol_id);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ REQUIRE(
+ same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+ }
+ }
+
+ SECTION("3D")
+ {
+ constexpr size_t Dimension = 3;
+
+ auto mesh_3d = MeshDataBaseForTests::get().hybrid3DMesh();
+
+ std::vector<std::shared_ptr<const IZoneDescriptor>> zone_list;
+ zone_list.push_back(std::make_shared<NamedZoneDescriptor>("LEFT"));
+
+ auto mesh_cell_zone = getMeshCellZone(*mesh_3d, *zone_list[0]);
+ CellValue<bool> is_cell_in_zone{mesh_3d->connectivity()};
+ is_cell_in_zone.fill(false);
+ auto zone_cell_list = mesh_cell_zone.cellList();
+ for (size_t i_cell = 0; i_cell < zone_cell_list.size(); ++i_cell) {
+ is_cell_in_zone[zone_cell_list[i_cell]] = true;
+ }
+
+ auto xj = MeshDataManager::instance().getMeshData(*mesh_3d).xj();
+
+ std::string_view data = R"(
+import math;
+let B_scalar_non_linear_3d: R^3 -> B, x -> (exp(2 * x[0])< 2*x[1]+x[2]);
+let N_scalar_non_linear_3d: R^3 -> N, x -> floor(3 * (x[0] + x[1]) * (x[0] + x[1]) + x[2] * x[2]);
+let Z_scalar_non_linear_3d: R^3 -> Z, x -> floor(exp(2 * x[0]) - 3 * x[1] + x[2]);
+let R_scalar_non_linear_3d: R^3 -> R, x -> 2 * exp(x[0]+x[2]) + 3 * x[1];
+let R1_non_linear_3d: R^3 -> R^1, x -> 2 * exp(x[0])+sin(x[1] + x[2]);
+let R2_non_linear_3d: R^3 -> R^2, x -> (2 * exp(x[0]), -3*x[1] * x[2]);
+let R3_non_linear_3d: R^3 -> R^3, x -> (2 * exp(x[0]) + 3, x[1] - 2, 3 * x[2]);
+let R1x1_non_linear_3d: R^3 -> R^1x1, x -> (2 * exp(x[0]) * sin(x[1]) + 3 * x[2]);
+let R2x2_non_linear_3d: R^3 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3, sin(x[2] - 2 * x[0]), 3, x[1] * x[0] - x[2]);
+let R3x3_non_linear_3d: R^3 -> R^3x3, x -> (2 * exp(x[0]) * sin(x[1]) + 3, sin(x[1] - 2 * x[2]), 3, x[1] * x[2], -4*x[1], 2*x[2]+1, 3, -6*x[2], exp(x[1] + x[2]));
+)";
+ TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+ auto ast = ASTBuilder::build(input);
+
+ ASTModulesImporter{*ast};
+ ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+ ASTSymbolTableBuilder{*ast};
+ ASTNodeDataTypeBuilder{*ast};
+
+ ASTNodeTypeCleaner<language::var_declaration>{*ast};
+ ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+ ASTNodeExpressionBuilder{*ast};
+
+ std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+ TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+ position.byte = data.size(); // ensure that variables are declared at this point
+
+ SECTION("B_scalar_non_linear_3d")
+ {
+ auto [i_symbol, found] = symbol_table->find("B_scalar_non_linear_3d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<double> cell_value{mesh_3d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = std::exp(2 * x[0]) < 2 * x[1] + x[2];
+ } else {
+ cell_value[cell_id] = false;
+ }
+ });
+
+ DiscreteFunctionInterpoler interpoler(mesh_3d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+ function_symbol_id);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ REQUIRE(
+ same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+ }
+
+ SECTION("N_scalar_non_linear_3d")
+ {
+ auto [i_symbol, found] = symbol_table->find("N_scalar_non_linear_3d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<double> cell_value{mesh_3d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = std::floor(3 * (x[0] + x[1]) * (x[0] + x[1]) + x[2] * x[2]);
+ } else {
+ cell_value[cell_id] = 0;
+ }
+ });
+
+ DiscreteFunctionInterpoler interpoler(mesh_3d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+ function_symbol_id);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ REQUIRE(
+ same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+ }
+
+ SECTION("Z_scalar_non_linear_3d")
+ {
+ auto [i_symbol, found] = symbol_table->find("Z_scalar_non_linear_3d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<double> cell_value{mesh_3d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = std::floor(std::exp(2 * x[0]) - 3 * x[1] + x[2]);
+ } else {
+ cell_value[cell_id] = 0;
+ }
+ });
+
+ DiscreteFunctionInterpoler interpoler(mesh_3d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+ function_symbol_id);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ REQUIRE(
+ same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+ }
+
+ SECTION("R_scalar_non_linear_3d")
+ {
+ auto [i_symbol, found] = symbol_table->find("R_scalar_non_linear_3d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<double> cell_value{mesh_3d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = 2 * std::exp(x[0] + x[2]) + 3 * x[1];
+ } else {
+ cell_value[cell_id] = 0;
+ }
+ });
+
+ DiscreteFunctionInterpoler interpoler(mesh_3d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+ function_symbol_id);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ REQUIRE(
+ same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+ }
+
+ SECTION("R1_non_linear_3d")
+ {
+ using DataType = TinyVector<1>;
+
+ auto [i_symbol, found] = symbol_table->find("R1_non_linear_3d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<DataType> cell_value{mesh_3d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = DataType{2 * std::exp(x[0]) + std::sin(x[1] + x[2])};
+ } else {
+ cell_value[cell_id] = zero;
+ }
+ });
+
+ DiscreteFunctionInterpoler interpoler(mesh_3d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+ function_symbol_id);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ REQUIRE(
+ same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+ }
+
+ SECTION("R2_non_linear_3d")
+ {
+ using DataType = TinyVector<2>;
+
+ auto [i_symbol, found] = symbol_table->find("R2_non_linear_3d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<DataType> cell_value{mesh_3d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = DataType{2 * std::exp(x[0]), -3 * x[1] * x[2]};
+ } else {
+ cell_value[cell_id] = zero;
+ }
+ });
+
+ DiscreteFunctionInterpoler interpoler(mesh_3d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+ function_symbol_id);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ REQUIRE(
+ same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+ }
+
+ SECTION("R3_non_linear_3d")
+ {
+ using DataType = TinyVector<3>;
+
+ auto [i_symbol, found] = symbol_table->find("R3_non_linear_3d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<DataType> cell_value{mesh_3d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = DataType{2 * std::exp(x[0]) + 3, x[1] - 2, 3 * x[2]};
+ } else {
+ cell_value[cell_id] = zero;
+ }
+ });
+
+ DiscreteFunctionInterpoler interpoler(mesh_3d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+ function_symbol_id);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ REQUIRE(
+ same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+ }
+
+ SECTION("R1x1_non_linear_3d")
+ {
+ using DataType = TinyMatrix<1>;
+
+ auto [i_symbol, found] = symbol_table->find("R1x1_non_linear_3d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<DataType> cell_value{mesh_3d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3 * x[2]};
+ } else {
+ cell_value[cell_id] = zero;
+ }
+ });
+
+ DiscreteFunctionInterpoler interpoler(mesh_3d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+ function_symbol_id);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ REQUIRE(
+ same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+ }
+
+ SECTION("R2x2_non_linear_3d")
+ {
+ using DataType = TinyMatrix<2>;
+
+ auto [i_symbol, found] = symbol_table->find("R2x2_non_linear_3d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<DataType> cell_value{mesh_3d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] =
+ DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3, std::sin(x[2] - 2 * x[0]), 3, x[1] * x[0] - x[2]};
+ } else {
+ cell_value[cell_id] = zero;
+ }
+ });
+
+ DiscreteFunctionInterpoler interpoler(mesh_3d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+ function_symbol_id);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ REQUIRE(
+ same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+ }
+
+ SECTION("R3x3_non_linear_3d")
+ {
+ using DataType = TinyMatrix<3>;
+
+ auto [i_symbol, found] = symbol_table->find("R3x3_non_linear_3d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<DataType> cell_value{mesh_3d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+
+ cell_value[cell_id] = DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3,
+ std::sin(x[1] - 2 * x[2]),
+ 3,
+ x[1] * x[2],
+ -4 * x[1],
+ 2 * x[2] + 1,
+ 3,
+ -6 * x[2],
+ std::exp(x[1] + x[2])};
+ } else {
+ cell_value[cell_id] = zero;
+ }
+ });
+
+ DiscreteFunctionInterpoler interpoler(mesh_3d, zone_list, std::make_shared<DiscreteFunctionDescriptorP0>(),
+ function_symbol_id);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ REQUIRE(
+ same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+ }
+ }
+}
diff --git a/tests/test_DiscreteFunctionVectorInterpolerByZone.cpp b/tests/test_DiscreteFunctionVectorInterpolerByZone.cpp
new file mode 100644
index 000000000..172c68d4f
--- /dev/null
+++ b/tests/test_DiscreteFunctionVectorInterpolerByZone.cpp
@@ -0,0 +1,511 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <language/ast/ASTBuilder.hpp>
+#include <language/ast/ASTModulesImporter.hpp>
+#include <language/ast/ASTNodeDataTypeBuilder.hpp>
+#include <language/ast/ASTNodeExpressionBuilder.hpp>
+#include <language/ast/ASTNodeFunctionEvaluationExpressionBuilder.hpp>
+#include <language/ast/ASTNodeFunctionExpressionBuilder.hpp>
+#include <language/ast/ASTNodeTypeCleaner.hpp>
+#include <language/ast/ASTSymbolTableBuilder.hpp>
+#include <language/utils/PugsFunctionAdapter.hpp>
+#include <language/utils/SymbolTable.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshCellZone.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <mesh/NamedZoneDescriptor.hpp>
+
+#include <scheme/DiscreteFunctionDescriptorP0Vector.hpp>
+#include <scheme/DiscreteFunctionP0Vector.hpp>
+#include <scheme/DiscreteFunctionVectorInterpoler.hpp>
+
+#include <pegtl/string_input.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("DiscreteFunctionVectorInterpolerByZone", "[scheme]")
+{
+ auto same_cell_value = [](const CellValue<const double>& fi, const size_t i, const auto& f) -> bool {
+ for (CellId cell_id = 0; cell_id < fi.numberOfItems(); ++cell_id) {
+ if (fi[cell_id] != f[cell_id][i]) {
+ return false;
+ }
+ }
+
+ return true;
+ };
+
+ auto register_function = [](const TAO_PEGTL_NAMESPACE::position& position,
+ const std::shared_ptr<SymbolTable>& symbol_table, const std::string& name,
+ std::vector<FunctionSymbolId>& function_id_list) {
+ auto [i_symbol, found] = symbol_table->find(name, 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);
+ function_id_list.push_back(function_symbol_id);
+ };
+
+ SECTION("1D")
+ {
+ constexpr size_t Dimension = 1;
+
+ auto mesh_1d = MeshDataBaseForTests::get().unordered1DMesh();
+
+ std::vector<std::shared_ptr<const IZoneDescriptor>> zone_list;
+ zone_list.push_back(std::make_shared<NamedZoneDescriptor>("LEFT"));
+
+ auto mesh_cell_zone = getMeshCellZone(*mesh_1d, *zone_list[0]);
+ CellValue<bool> is_cell_in_zone{mesh_1d->connectivity()};
+ is_cell_in_zone.fill(false);
+ auto zone_cell_list = mesh_cell_zone.cellList();
+ for (size_t i_cell = 0; i_cell < zone_cell_list.size(); ++i_cell) {
+ is_cell_in_zone[zone_cell_list[i_cell]] = true;
+ }
+
+ auto xj = MeshDataManager::instance().getMeshData(*mesh_1d).xj();
+
+ std::string_view data = R"(
+import math;
+let B_scalar_non_linear_1d: R^1 -> B, x -> (exp(2 * x[0]) + 3 > 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;
+)";
+ 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};
+
+ TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+ position.byte = data.size(); // ensure that variables are declared at this point
+
+ std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+ std::vector<FunctionSymbolId> function_id_list;
+ register_function(position, symbol_table, "B_scalar_non_linear_1d", function_id_list);
+ register_function(position, symbol_table, "N_scalar_non_linear_1d", function_id_list);
+ register_function(position, symbol_table, "Z_scalar_non_linear_1d", function_id_list);
+ register_function(position, symbol_table, "R_scalar_non_linear_1d", function_id_list);
+
+ DiscreteFunctionVectorInterpoler interpoler(mesh_1d, zone_list,
+ std::make_shared<DiscreteFunctionDescriptorP0Vector>(),
+ function_id_list);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ size_t i = 0;
+
+ {
+ CellValue<double> cell_value{mesh_1d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = std::exp(2 * x[0]) + 3 > 4;
+ } else {
+ cell_value[cell_id] = 0;
+ }
+ });
+
+ REQUIRE(same_cell_value(cell_value, i++,
+ dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+ }
+
+ {
+ CellValue<double> cell_value{mesh_1d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = std::floor(3 * x[0] * x[0] + 2);
+ } else {
+ cell_value[cell_id] = 0;
+ }
+ });
+
+ REQUIRE(same_cell_value(cell_value, i++,
+ dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+ }
+
+ {
+ CellValue<double> cell_value{mesh_1d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = std::floor(std::exp(2 * x[0]) - 1);
+ } else {
+ cell_value[cell_id] = 0;
+ }
+ });
+
+ REQUIRE(same_cell_value(cell_value, i++,
+ dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+ }
+
+ {
+ CellValue<double> cell_value{mesh_1d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = 2 * std::exp(x[0]) + 3;
+ } else {
+ cell_value[cell_id] = 0;
+ }
+ });
+
+ REQUIRE(same_cell_value(cell_value, i++,
+ dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+ }
+
+ REQUIRE(i == function_id_list.size());
+ }
+
+ SECTION("2D")
+ {
+ constexpr size_t Dimension = 2;
+
+ auto mesh_2d = MeshDataBaseForTests::get().hybrid2DMesh();
+
+ std::vector<std::shared_ptr<const IZoneDescriptor>> zone_list;
+ zone_list.push_back(std::make_shared<NamedZoneDescriptor>("LEFT"));
+
+ auto mesh_cell_zone = getMeshCellZone(*mesh_2d, *zone_list[0]);
+ CellValue<bool> is_cell_in_zone{mesh_2d->connectivity()};
+ is_cell_in_zone.fill(false);
+ auto zone_cell_list = mesh_cell_zone.cellList();
+ for (size_t i_cell = 0; i_cell < zone_cell_list.size(); ++i_cell) {
+ is_cell_in_zone[zone_cell_list[i_cell]] = true;
+ }
+
+ auto xj = MeshDataManager::instance().getMeshData(*mesh_2d).xj();
+
+ std::string_view data = R"(
+import math;
+let B_scalar_non_linear_2d: R^2 -> B, x -> (exp(2 * x[0]) + 3 > 4);
+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[1]) - 1);
+let R_scalar_non_linear_2d: R^2 -> R, x -> 2 * exp(x[0] + 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};
+
+ TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+ position.byte = data.size(); // ensure that variables are declared at this point
+
+ std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+ std::vector<FunctionSymbolId> function_id_list;
+ register_function(position, symbol_table, "B_scalar_non_linear_2d", function_id_list);
+ register_function(position, symbol_table, "N_scalar_non_linear_2d", function_id_list);
+ register_function(position, symbol_table, "Z_scalar_non_linear_2d", function_id_list);
+ register_function(position, symbol_table, "R_scalar_non_linear_2d", function_id_list);
+
+ DiscreteFunctionVectorInterpoler interpoler(mesh_2d, zone_list,
+ std::make_shared<DiscreteFunctionDescriptorP0Vector>(),
+ function_id_list);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ size_t i = 0;
+
+ {
+ CellValue<double> cell_value{mesh_2d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = std::exp(2 * x[0]) + 3 > 4;
+ } else {
+ cell_value[cell_id] = 0;
+ }
+ });
+
+ REQUIRE(same_cell_value(cell_value, i++,
+ dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+ }
+
+ {
+ CellValue<double> cell_value{mesh_2d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = std::floor(3 * (x[0] * x[1]) * (x[0] * x[1]) + 2);
+ } else {
+ cell_value[cell_id] = 0;
+ }
+ });
+
+ REQUIRE(same_cell_value(cell_value, i++,
+ dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+ }
+
+ {
+ CellValue<double> cell_value{mesh_2d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = std::floor(std::exp(2 * x[1]) - 1);
+ } else {
+ cell_value[cell_id] = 0;
+ }
+ });
+
+ REQUIRE(same_cell_value(cell_value, i++,
+ dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+ }
+
+ {
+ CellValue<double> cell_value{mesh_2d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = 2 * std::exp(x[0] + x[1]) + 3;
+ } else {
+ cell_value[cell_id] = 0;
+ }
+ });
+
+ REQUIRE(same_cell_value(cell_value, i++,
+ dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+ }
+
+ REQUIRE(i == function_id_list.size());
+ }
+
+ SECTION("3D")
+ {
+ constexpr size_t Dimension = 3;
+
+ auto mesh_3d = MeshDataBaseForTests::get().hybrid3DMesh();
+
+ std::vector<std::shared_ptr<const IZoneDescriptor>> zone_list;
+ zone_list.push_back(std::make_shared<NamedZoneDescriptor>("LEFT"));
+
+ auto mesh_cell_zone = getMeshCellZone(*mesh_3d, *zone_list[0]);
+ CellValue<bool> is_cell_in_zone{mesh_3d->connectivity()};
+ is_cell_in_zone.fill(false);
+ auto zone_cell_list = mesh_cell_zone.cellList();
+ for (size_t i_cell = 0; i_cell < zone_cell_list.size(); ++i_cell) {
+ is_cell_in_zone[zone_cell_list[i_cell]] = true;
+ }
+
+ auto xj = MeshDataManager::instance().getMeshData(*mesh_3d).xj();
+
+ std::string_view data = R"(
+import math;
+let B_scalar_non_linear_3d: R^3 -> B, x -> (exp(2 * x[0] + x[2]) + 3 > 4);
+let N_scalar_non_linear_3d: R^3 -> N, x -> floor(3 * (x[0] * x[1]) * (x[0] * x[1]) + 2);
+let Z_scalar_non_linear_3d: R^3 -> Z, x -> floor(exp(2 * x[1]) - x[2]);
+let R_scalar_non_linear_3d: R^3 -> R, x -> 2 * exp(x[0] + x[1]) + 3 * 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};
+
+ TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+ position.byte = data.size(); // ensure that variables are declared at this point
+
+ std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+ std::vector<FunctionSymbolId> function_id_list;
+ register_function(position, symbol_table, "B_scalar_non_linear_3d", function_id_list);
+ register_function(position, symbol_table, "N_scalar_non_linear_3d", function_id_list);
+ register_function(position, symbol_table, "Z_scalar_non_linear_3d", function_id_list);
+ register_function(position, symbol_table, "R_scalar_non_linear_3d", function_id_list);
+
+ DiscreteFunctionVectorInterpoler interpoler(mesh_3d, zone_list,
+ std::make_shared<DiscreteFunctionDescriptorP0Vector>(),
+ function_id_list);
+ std::shared_ptr discrete_function = interpoler.interpolate();
+
+ size_t i = 0;
+
+ {
+ CellValue<double> cell_value{mesh_3d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = std::exp(2 * x[0] + x[2]) + 3 > 4;
+ } else {
+ cell_value[cell_id] = 0;
+ }
+ });
+
+ REQUIRE(same_cell_value(cell_value, i++,
+ dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+ }
+
+ {
+ CellValue<double> cell_value{mesh_3d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = std::floor(3 * (x[0] * x[1]) * (x[0] * x[1]) + 2);
+ } else {
+ cell_value[cell_id] = 0;
+ }
+ });
+
+ REQUIRE(same_cell_value(cell_value, i++,
+ dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+ }
+
+ {
+ CellValue<double> cell_value{mesh_3d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = std::floor(std::exp(2 * x[1]) - x[2]);
+ } else {
+ cell_value[cell_id] = 0;
+ }
+ });
+
+ REQUIRE(same_cell_value(cell_value, i++,
+ dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+ }
+
+ {
+ CellValue<double> cell_value{mesh_3d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ if (is_cell_in_zone[cell_id]) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = 2 * std::exp(x[0] + x[1]) + 3 * x[2];
+ } else {
+ cell_value[cell_id] = 0;
+ }
+ });
+
+ REQUIRE(same_cell_value(cell_value, i++,
+ dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*discrete_function)));
+ }
+
+ REQUIRE(i == function_id_list.size());
+ }
+
+ SECTION("errors")
+ {
+ std::array mesh_list = MeshDataBaseForTests::get().all3DMeshes();
+
+ for (auto named_mesh : mesh_list) {
+ SECTION(named_mesh.name())
+ {
+ auto mesh_3d = named_mesh.mesh();
+
+ auto xj = MeshDataManager::instance().getMeshData(*mesh_3d).xj();
+
+ std::string_view data = R"(
+import math;
+let B_scalar_non_linear_2d: R^2 -> B, x -> (exp(2 * x[0] + x[1]) + 3 > 4);
+let N_scalar_non_linear_1d: R^1 -> N, x -> floor(3 * x[0] * x[0] + 2);
+let Z_scalar_non_linear_3d: R^3 -> Z, x -> floor(exp(2 * x[1]) - x[2]);
+let R_scalar_non_linear_3d: R^3 -> R, x -> 2 * exp(x[0] + x[1]) + 3 * x[2];
+let R2_scalar_non_linear_3d: R^3 -> R^2, x -> (2 * exp(x[0] + x[1]) + 3 * x[2], x[0] - 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};
+
+ TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+ position.byte = data.size(); // ensure that variables are declared at this point
+
+ std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+ SECTION("invalid function type")
+ {
+ std::vector<FunctionSymbolId> function_id_list;
+ register_function(position, symbol_table, "B_scalar_non_linear_2d", function_id_list);
+ register_function(position, symbol_table, "N_scalar_non_linear_1d", function_id_list);
+ register_function(position, symbol_table, "Z_scalar_non_linear_3d", function_id_list);
+ register_function(position, symbol_table, "R_scalar_non_linear_3d", function_id_list);
+
+ DiscreteFunctionVectorInterpoler interpoler(mesh_3d, std::make_shared<DiscreteFunctionDescriptorP0Vector>(),
+ function_id_list);
+
+ const std::string error_msg = R"(error: invalid function type
+note: expecting R^3 -> R
+note: provided function B_scalar_non_linear_2d: R^2 -> B)";
+
+ REQUIRE_THROWS_WITH(interpoler.interpolate(), error_msg);
+ }
+
+ SECTION("invalid value type")
+ {
+ std::vector<FunctionSymbolId> function_id_list;
+ register_function(position, symbol_table, "Z_scalar_non_linear_3d", function_id_list);
+ register_function(position, symbol_table, "R_scalar_non_linear_3d", function_id_list);
+ register_function(position, symbol_table, "R2_scalar_non_linear_3d", function_id_list);
+
+ DiscreteFunctionVectorInterpoler interpoler(mesh_3d, std::make_shared<DiscreteFunctionDescriptorP0Vector>(),
+ function_id_list);
+
+ const std::string error_msg = R"(error: vector functions require scalar value type.
+Invalid interpolation value type: R^2)";
+
+ REQUIRE_THROWS_WITH(interpoler.interpolate(), error_msg);
+ }
+
+ SECTION("invalid discrete function type")
+ {
+ const std::string error_msg = "error: invalid discrete function type for vector interpolation";
+
+ DiscreteFunctionVectorInterpoler interpoler{mesh_3d, std::make_shared<DiscreteFunctionDescriptorP0>(), {}};
+ REQUIRE_THROWS_WITH(interpoler.interpolate(), error_msg);
+ }
+ }
+ }
+ }
+}
--
GitLab