diff --git a/.gitignore b/.gitignore
index 3bb6d6cbcf2906a99d3302595285d50145395fe5..96214f7a18cd5cd618c43f91b02e940897be0c42 100644
--- a/.gitignore
+++ b/.gitignore
@@ -12,3 +12,4 @@ GPATH
GRTAGS
GTAGS
/.clangd/
+/.cache/
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 45c69b7a0d8744a1fde62b53170e7f3f8ef303c2..114ceecc5d49685a513708dced3dcffea83bff29 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -528,7 +528,9 @@ target_link_libraries(
PugsLanguageAlgorithms
PugsMesh
PugsAlgebra
+ PugsScheme
PugsUtils
+ PugsOutput
PugsLanguageUtils
kokkos
${PETSC_LIBRARIES}
@@ -551,7 +553,8 @@ if(${CMAKE_VERSION} VERSION_LESS "3.13.0")
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib)
else()
- install(TARGETS pugs PugsMesh PugsUtils PugsLanguage
+ install(TARGETS pugs PugsAlgebra PugsUtils PugsLanguage PugsLanguageAST PugsLanguageModules PugsLanguageAlgorithms PugsLanguageUtils PugsMesh PugsScheme PugsOutput
+
RUNTIME DESTINATION bin
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib)
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 1b21f79a7545c442be5448792132acbdf840947a..16904a7be0cddc9aa54206638ab4b5c759cd0de8 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -15,5 +15,8 @@ add_subdirectory(algebra)
# Pugs mesh
add_subdirectory(mesh)
+# Pugs mesh
+add_subdirectory(scheme)
+
# Pugs output
-#add_subdirectory(output)
+add_subdirectory(output)
diff --git a/src/language/algorithms/AcousticSolverAlgorithm.cpp b/src/language/algorithms/AcousticSolverAlgorithm.cpp
index 9811c90f0b0f2c52d5f72369efea4a4398f007e3..edce9e025f0feb4cd39acaadcb537451c9251065 100644
--- a/src/language/algorithms/AcousticSolverAlgorithm.cpp
+++ b/src/language/algorithms/AcousticSolverAlgorithm.cpp
@@ -1,12 +1,13 @@
#include <language/algorithms/AcousticSolverAlgorithm.hpp>
#include <language/utils/InterpolateItemValue.hpp>
-#include <output/VTKWriter.hpp>
#include <scheme/AcousticSolver.hpp>
#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
#include <scheme/IBoundaryDescriptor.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+#include <iomanip>
+
template <size_t Dimension>
AcousticSolverAlgorithm<Dimension>::AcousticSolverAlgorithm(
const AcousticSolverType& solver_type,
@@ -176,18 +177,9 @@ AcousticSolverAlgorithm<Dimension>::AcousticSolverAlgorithm(
mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { inv_mj[j] = 1. / mj[j]; });
}
- VTKWriter vtk_writer("mesh_" + std::to_string(Dimension), 0.01);
-
while ((t < tmax) and (iteration < itermax)) {
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
- vtk_writer.write(mesh,
- {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
- NamedItemValue{"coords", mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()},
- NamedItemValue{"cell_owner", mesh->connectivity().cellOwner()},
- NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}},
- t);
-
double dt = 0.4 * acoustic_solver.acoustic_dt(mesh_data.Vj(), cj);
if (t + dt > tmax) {
dt = tmax - t;
@@ -208,16 +200,6 @@ AcousticSolverAlgorithm<Dimension>::AcousticSolverAlgorithm(
std::cout << rang::style::bold << "Final time=" << rang::fgB::green << t << rang::style::reset << " reached after "
<< rang::fgB::cyan << iteration << rang::style::reset << rang::style::bold << " iterations"
<< rang::style::reset << '\n';
- {
- MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
-
- vtk_writer.write(mesh,
- {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
- NamedItemValue{"coords", mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()},
- NamedItemValue{"cell_owner", mesh->connectivity().cellOwner()},
- NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}},
- t, true); // forces last output
- }
}
template AcousticSolverAlgorithm<1>::AcousticSolverAlgorithm(
diff --git a/src/language/ast/ASTNodeAffectationExpressionBuilder.cpp b/src/language/ast/ASTNodeAffectationExpressionBuilder.cpp
index 884711c081ded349d32663f80180c06c4911bfba..a97c13443bd628a2226de54623c673d464916a0c 100644
--- a/src/language/ast/ASTNodeAffectationExpressionBuilder.cpp
+++ b/src/language/ast/ASTNodeAffectationExpressionBuilder.cpp
@@ -9,24 +9,24 @@
#include <language/utils/ParseError.hpp>
#include <utils/Exceptions.hpp>
-ASTNodeAffectationExpressionBuilder::ASTNodeAffectationExpressionBuilder(ASTNode& n)
+ASTNodeAffectationExpressionBuilder::ASTNodeAffectationExpressionBuilder(ASTNode& node)
{
- const ASTNodeDataType& target_data_type = n.children[0]->m_data_type;
- const ASTNodeDataType& source_data_type = n.children[1]->m_data_type;
+ const ASTNodeDataType& target_data_type = node.children[0]->m_data_type;
+ const ASTNodeDataType& source_data_type = node.children[1]->m_data_type;
const std::string affectation_name = [&] {
- if (n.is_type<language::eq_op>()) {
+ if (node.is_type<language::eq_op>()) {
return affectationMangler<language::eq_op>(target_data_type, source_data_type);
- } else if (n.is_type<language::multiplyeq_op>()) {
+ } else if (node.is_type<language::multiplyeq_op>()) {
return affectationMangler<language::multiplyeq_op>(target_data_type, source_data_type);
- } else if (n.is_type<language::divideeq_op>()) {
+ } else if (node.is_type<language::divideeq_op>()) {
return affectationMangler<language::divideeq_op>(target_data_type, source_data_type);
- } else if (n.is_type<language::pluseq_op>()) {
+ } else if (node.is_type<language::pluseq_op>()) {
return affectationMangler<language::pluseq_op>(target_data_type, source_data_type);
- } else if (n.is_type<language::minuseq_op>()) {
+ } else if (node.is_type<language::minuseq_op>()) {
return affectationMangler<language::minuseq_op>(target_data_type, source_data_type);
} else {
- throw ParseError("unexpected error: undefined affectation operator", std::vector{n.begin()});
+ throw ParseError("unexpected error: undefined affectation operator", std::vector{node.begin()});
}
}();
@@ -34,12 +34,15 @@ ASTNodeAffectationExpressionBuilder::ASTNodeAffectationExpressionBuilder(ASTNode
OperatorRepository::instance().getAffectationProcessorBuilder(affectation_name);
if (optional_processor_builder.has_value()) {
- n.m_node_processor = optional_processor_builder.value()->getNodeProcessor(n);
+ ASTNode& lhs_node = *node.children[0];
+ ASTNode& rhs_node = *node.children[1];
+
+ node.m_node_processor = optional_processor_builder.value()->getNodeProcessor(lhs_node, rhs_node);
} else {
std::ostringstream error_message;
error_message << "undefined affectation type: ";
error_message << rang::fgB::red << affectation_name << rang::fg::reset;
- throw ParseError(error_message.str(), std::vector{n.children[0]->begin()});
+ throw ParseError(error_message.str(), std::vector{node.begin()});
}
}
diff --git a/src/language/ast/ASTNodeDataTypeBuilder.cpp b/src/language/ast/ASTNodeDataTypeBuilder.cpp
index 018bc94f74160bfcd680df08ac9b2a08c54ba492..b5f65437ce1b45b4a403aaf079317de8f04c0de5 100644
--- a/src/language/ast/ASTNodeDataTypeBuilder.cpp
+++ b/src/language/ast/ASTNodeDataTypeBuilder.cpp
@@ -333,6 +333,22 @@ ASTNodeDataTypeBuilder::_buildNodeDataTypes(ASTNode& n) const
value_type = getMatrixDataType(image_node);
} else if (image_node.is_type<language::string_type>()) {
value_type = ASTNodeDataType::build<ASTNodeDataType::string_t>();
+ } else if (image_node.is_type<language::type_name_id>()) {
+ const std::string& type_name_id = image_node.string();
+
+ auto& symbol_table = *image_node.m_symbol_table;
+
+ const auto [i_type_symbol, found] = symbol_table.find(type_name_id, image_node.begin());
+ if (not found) {
+ throw ParseError("undefined type identifier", std::vector{image_node.begin()});
+ } else if (i_type_symbol->attributes().dataType() != ASTNodeDataType::type_name_id_t) {
+ std::ostringstream os;
+ os << "invalid type identifier, '" << type_name_id << "' was previously defined as a '"
+ << dataTypeName(i_type_symbol->attributes().dataType()) << '\'';
+ throw ParseError(os.str(), std::vector{image_node.begin()});
+ }
+
+ value_type = ASTNodeDataType::build<ASTNodeDataType::type_id_t>(type_name_id);
}
// LCOV_EXCL_START
@@ -348,7 +364,6 @@ ASTNodeDataTypeBuilder::_buildNodeDataTypes(ASTNode& n) const
}
n.m_data_type = ASTNodeDataType::build<ASTNodeDataType::typename_t>(
ASTNodeDataType::build<ASTNodeDataType::list_t>(sub_data_type_list));
-
} else if (n.is_type<language::for_post>() or n.is_type<language::for_init>() or
n.is_type<language::for_statement_block>()) {
n.m_data_type = ASTNodeDataType::build<ASTNodeDataType::void_t>();
@@ -361,13 +376,11 @@ ASTNodeDataTypeBuilder::_buildNodeDataTypes(ASTNode& n) const
const ASTNode& test_node = *n.children[0];
ASTNodeNaturalConversionChecker{test_node, ASTNodeDataType::build<ASTNodeDataType::bool_t>()};
-
} else if (n.is_type<language::do_while_statement>()) {
n.m_data_type = ASTNodeDataType::build<ASTNodeDataType::void_t>();
const ASTNode& test_node = *n.children[1];
ASTNodeNaturalConversionChecker{test_node, ASTNodeDataType::build<ASTNodeDataType::bool_t>()};
-
} else if (n.is_type<language::unary_not>() or n.is_type<language::unary_minus>()) {
auto& operator_repository = OperatorRepository::instance();
@@ -394,7 +407,6 @@ ASTNodeDataTypeBuilder::_buildNodeDataTypes(ASTNode& n) const
<< rang::style::reset;
throw ParseError(message.str(), n.begin());
}
-
} else if (n.is_type<language::unary_plusplus>() or n.is_type<language::unary_minusminus>() or
n.is_type<language::post_plusplus>() or n.is_type<language::post_minusminus>()) {
auto& operator_repository = OperatorRepository::instance();
@@ -428,7 +440,6 @@ ASTNodeDataTypeBuilder::_buildNodeDataTypes(ASTNode& n) const
<< rang::style::reset;
throw ParseError(message.str(), n.begin());
}
-
} else if (n.is_type<language::plus_op>() or n.is_type<language::minus_op>() or
n.is_type<language::multiply_op>() or n.is_type<language::divide_op>() or
n.is_type<language::lesser_op>() or n.is_type<language::lesser_or_eq_op>() or
@@ -496,7 +507,6 @@ ASTNodeDataTypeBuilder::_buildNodeDataTypes(ASTNode& n) const
<< "note: incompatible operand types " << dataTypeName(type_0) << " and " << dataTypeName(type_1);
throw ParseError(message.str(), n.begin());
}
-
} else if (n.is_type<language::function_evaluation>()) {
if (n.children[0]->m_data_type == ASTNodeDataType::function_t) {
const std::string& function_name = n.children[0]->string();
diff --git a/src/language/ast/ASTNodeDataTypeFlattener.cpp b/src/language/ast/ASTNodeDataTypeFlattener.cpp
index 51d055836393b1ee680b7d983756bee39020d992..6a06c72bb66926af546ff0a130315eacc8b54c9a 100644
--- a/src/language/ast/ASTNodeDataTypeFlattener.cpp
+++ b/src/language/ast/ASTNodeDataTypeFlattener.cpp
@@ -1,6 +1,7 @@
#include <language/ast/ASTNodeDataTypeFlattener.hpp>
#include <language/PEGGrammar.hpp>
+#include <language/utils/BuiltinFunctionEmbedder.hpp>
#include <language/utils/FunctionTable.hpp>
#include <language/utils/SymbolTable.hpp>
@@ -33,9 +34,20 @@ ASTNodeDataTypeFlattener::ASTNodeDataTypeFlattener(ASTNode& node, FlattenedDataT
}
break;
}
+ case ASTNodeDataType::builtin_function_t: {
+ uint64_t builtin_function_id = std::get<uint64_t>(i_function_symbol->attributes().value());
+ auto builtin_function_embedder = node.m_symbol_table->builtinFunctionEmbedderTable()[builtin_function_id];
+
+ const auto& compound_data_type = builtin_function_embedder->getReturnDataType();
+ for (auto data_type : compound_data_type.contentTypeList()) {
+ flattened_datatype_list.push_back({*data_type, node});
+ }
+
+ break;
+ }
// LCOV_EXCL_START
default: {
- throw ParseError("unexpected function type", node.begin());
+ throw ParseError{"unexpected function type", node.begin()};
}
// LCOV_EXCL_STOP
}
diff --git a/src/language/ast/ASTNodeListAffectationExpressionBuilder.cpp b/src/language/ast/ASTNodeListAffectationExpressionBuilder.cpp
index b135bbcd3f67127c23c044880d96ad28ee6527c5..6b53db72a1a257f74196a56f6e2c83694531fd5f 100644
--- a/src/language/ast/ASTNodeListAffectationExpressionBuilder.cpp
+++ b/src/language/ast/ASTNodeListAffectationExpressionBuilder.cpp
@@ -6,6 +6,9 @@
#include <language/utils/ASTNodeNaturalConversionChecker.hpp>
#include <language/utils/ParseError.hpp>
+#include <language/utils/AffectationMangler.hpp>
+#include <language/utils/OperatorRepository.hpp>
+
template <typename OperatorT>
void
ASTNodeListAffectationExpressionBuilder::_buildAffectationProcessor(
@@ -115,6 +118,25 @@ ASTNodeListAffectationExpressionBuilder::_buildAffectationProcessor(
}
};
+ auto add_affectation_processor_for_embedded_data = [&](const ASTNodeSubDataType& node_sub_data_type) {
+ if constexpr (std::is_same_v<OperatorT, language::eq_op>) {
+ switch (node_sub_data_type.m_data_type) {
+ case ASTNodeDataType::type_id_t: {
+ list_affectation_processor->template add<EmbeddedData, EmbeddedData>(value_node);
+ break;
+ }
+ // LCOV_EXCL_START
+ default: {
+ throw ParseError("unexpected error:invalid operand type for embedded data affectation",
+ std::vector{node_sub_data_type.m_parent_node.begin()});
+ }
+ // LCOV_EXCL_STOP
+ }
+ } else {
+ throw ParseError("unexpected error: undefined operator type for string affectation", std::vector{m_node.begin()});
+ }
+ };
+
auto add_affectation_processor_for_string_data = [&](const ASTNodeSubDataType& node_sub_data_type) {
if constexpr (std::is_same_v<OperatorT, language::eq_op>) {
switch (node_sub_data_type.m_data_type) {
@@ -192,6 +214,10 @@ ASTNodeListAffectationExpressionBuilder::_buildAffectationProcessor(
add_affectation_processor_for_data(double{}, node_sub_data_type);
break;
}
+ case ASTNodeDataType::type_id_t: {
+ add_affectation_processor_for_embedded_data(node_sub_data_type);
+ break;
+ }
case ASTNodeDataType::vector_t: {
switch (value_type.dimension()) {
case 1: {
@@ -251,7 +277,21 @@ ASTNodeListAffectationExpressionBuilder::_buildAffectationProcessor(
ASTNodeNaturalConversionChecker<AllowRToR1Conversion>(rhs_node_sub_data_type, value_node.m_data_type);
- add_affectation_processor_for_value(value_node.m_data_type, rhs_node_sub_data_type);
+ const std::string affectation_name =
+ affectationMangler<language::eq_op>(value_node.m_data_type, rhs_node_sub_data_type.m_data_type);
+
+ const auto& optional_processor_builder =
+ OperatorRepository::instance().getAffectationProcessorBuilder(affectation_name);
+
+ if (optional_processor_builder.has_value()) {
+ add_affectation_processor_for_value(value_node.m_data_type, rhs_node_sub_data_type);
+ } else {
+ std::ostringstream error_message;
+ error_message << "undefined affectation type: ";
+ error_message << rang::fgB::red << affectation_name << rang::fg::reset;
+
+ throw ParseError(error_message.str(), std::vector{m_node.children[0]->begin()});
+ }
}
template <typename OperatorT>
diff --git a/src/language/modules/CMakeLists.txt b/src/language/modules/CMakeLists.txt
index 020727478c5e64bac9ec0a780a0cfb80541fb28a..dfc777168d454337ffbd38ef084f37f50924b117 100644
--- a/src/language/modules/CMakeLists.txt
+++ b/src/language/modules/CMakeLists.txt
@@ -9,7 +9,7 @@ add_library(PugsLanguageModules
ModuleRepository.cpp
SchemeModule.cpp
UtilsModule.cpp
- VTKModule.cpp
+ WriterModule.cpp
)
diff --git a/src/language/modules/ModuleRepository.cpp b/src/language/modules/ModuleRepository.cpp
index 65bfcafd1f1c85f077a97f6e932fa97da0079537..44a437d04dc177cbfc6274420999c94fafaa22c9 100644
--- a/src/language/modules/ModuleRepository.cpp
+++ b/src/language/modules/ModuleRepository.cpp
@@ -7,7 +7,7 @@
#include <language/modules/MeshModule.hpp>
#include <language/modules/SchemeModule.hpp>
#include <language/modules/UtilsModule.hpp>
-#include <language/modules/VTKModule.hpp>
+#include <language/modules/WriterModule.hpp>
#include <language/utils/BasicAffectationRegistrerFor.hpp>
#include <language/utils/BuiltinFunctionEmbedder.hpp>
#include <language/utils/ParseError.hpp>
@@ -32,7 +32,7 @@ ModuleRepository::ModuleRepository()
this->_subscribe(std::make_unique<MeshModule>());
this->_subscribe(std::make_unique<SchemeModule>());
this->_subscribe(std::make_unique<UtilsModule>());
- this->_subscribe(std::make_unique<VTKModule>());
+ this->_subscribe(std::make_unique<WriterModule>());
}
template <typename NameEmbedderMapT, typename EmbedderTableT>
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 14f798eb6bac43b77fd98367a009fa19d5dcfcbd..8277ff9fc0ed6c2035520966bcdf304b794f4aeb 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -6,10 +6,14 @@
#include <language/utils/TypeDescriptor.hpp>
#include <mesh/Mesh.hpp>
#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
+#include <scheme/DiscreteFunctionDescriptorP0.hpp>
+#include <scheme/DiscreteFunctionInterpoler.hpp>
#include <scheme/FourierBoundaryConditionDescriptor.hpp>
#include <scheme/FreeBoundaryConditionDescriptor.hpp>
#include <scheme/IBoundaryConditionDescriptor.hpp>
#include <scheme/IBoundaryDescriptor.hpp>
+#include <scheme/IDiscreteFunction.hpp>
+#include <scheme/IDiscreteFunctionDescriptor.hpp>
#include <scheme/NamedBoundaryDescriptor.hpp>
#include <scheme/NeumannBoundaryConditionDescriptor.hpp>
#include <scheme/NumberedBoundaryDescriptor.hpp>
@@ -19,9 +23,33 @@
SchemeModule::SchemeModule()
{
+ this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const IDiscreteFunction>>);
+ this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const IDiscreteFunctionDescriptor>>);
+
this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const IBoundaryDescriptor>>);
this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const IBoundaryConditionDescriptor>>);
+ this->_addBuiltinFunction("P0", std::make_shared<
+ BuiltinFunctionEmbedder<std::shared_ptr<const IDiscreteFunctionDescriptor>()>>(
+ []() -> std::shared_ptr<const IDiscreteFunctionDescriptor> {
+ return std::make_shared<DiscreteFunctionDescriptorP0>();
+ }
+
+ ));
+
+ this->_addBuiltinFunction(
+ "interpolate",
+ std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
+ const IDiscreteFunction>(std::shared_ptr<const IMesh>, std::shared_ptr<const IDiscreteFunctionDescriptor>,
+ const FunctionSymbolId&)>>(
+ [](std::shared_ptr<const IMesh> mesh,
+ std::shared_ptr<const IDiscreteFunctionDescriptor> discrete_function_descriptor,
+ const FunctionSymbolId& function_id) -> std::shared_ptr<const IDiscreteFunction> {
+ return DiscreteFunctionInterpoler{mesh, discrete_function_descriptor, function_id}.interpolate();
+ }
+
+ ));
+
this->_addBuiltinFunction("boundaryName",
std::make_shared<
BuiltinFunctionEmbedder<std::shared_ptr<const IBoundaryDescriptor>(const std::string&)>>(
diff --git a/src/language/modules/SchemeModule.hpp b/src/language/modules/SchemeModule.hpp
index aa0f1322618956da9b3a4bc095c87f95ebf719fa..c49d1d4e8bc7fdf19430dac36d8bde4709eeda1f 100644
--- a/src/language/modules/SchemeModule.hpp
+++ b/src/language/modules/SchemeModule.hpp
@@ -20,6 +20,16 @@ template <>
inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const BasisType>> =
ASTNodeDataType::build<ASTNodeDataType::type_id_t>("basis_type");
+class IDiscreteFunction;
+template <>
+inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const IDiscreteFunction>> =
+ ASTNodeDataType::build<ASTNodeDataType::type_id_t>("Vh");
+
+class IDiscreteFunctionDescriptor;
+template <>
+inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const IDiscreteFunctionDescriptor>> =
+ ASTNodeDataType::build<ASTNodeDataType::type_id_t>("Vh_type");
+
class SchemeModule : public BuiltinModule
{
public:
diff --git a/src/language/modules/UtilsModule.cpp b/src/language/modules/UtilsModule.cpp
index 8b2271158c1548d05b39c379782abb3cc4754d66..b74ddc54eaab21ebf44dd53dbb585f4d0111c3fe 100644
--- a/src/language/modules/UtilsModule.cpp
+++ b/src/language/modules/UtilsModule.cpp
@@ -52,9 +52,7 @@ UtilsModule::UtilsModule()
std::make_shared<BuiltinFunctionEmbedder<std::string(const FunctionSymbolId&)>>(
[](const FunctionSymbolId& function_symbol_id) -> std::string {
- auto& function_table = function_symbol_id.symbolTable().functionTable();
-
- const auto& function_descriptor = function_table[function_symbol_id.id()];
+ const auto& function_descriptor = function_symbol_id.descriptor();
std::ostringstream os;
os << function_descriptor.name() << ": domain mapping\n";
diff --git a/src/language/modules/VTKModule.cpp b/src/language/modules/VTKModule.cpp
deleted file mode 100644
index ce84141ee866f7ac6fdbb8448db78d0d1b64de88..0000000000000000000000000000000000000000
--- a/src/language/modules/VTKModule.cpp
+++ /dev/null
@@ -1,61 +0,0 @@
-#include <language/modules/VTKModule.hpp>
-
-#include <language/utils/BuiltinFunctionEmbedder.hpp>
-#include <language/utils/TypeDescriptor.hpp>
-#include <mesh/Connectivity.hpp>
-#include <mesh/GmshReader.hpp>
-#include <mesh/Mesh.hpp>
-#include <output/VTKWriter.hpp>
-
-VTKModule::VTKModule()
-{
- this->_addBuiltinFunction("writeVTK",
- std::make_shared<
- BuiltinFunctionEmbedder<void(std::shared_ptr<const IMesh>, const std::string&)>>(
-
- [](std::shared_ptr<const IMesh> p_mesh, const std::string& filename) -> void {
- VTKWriter writer(filename, 0.1);
-
- static double time = 0;
-
- switch (p_mesh->dimension()) {
- case 1: {
- using MeshType = Mesh<Connectivity<1>>;
- const std::shared_ptr<const MeshType> mesh =
- std::dynamic_pointer_cast<const MeshType>(p_mesh);
-
- writer.write(mesh,
- OutputNamedItemValueSet{
- NamedItemValue{"cell_number", mesh->connectivity().cellNumber()}},
- time, true);
- break;
- }
- case 2: {
- using MeshType = Mesh<Connectivity<2>>;
- const std::shared_ptr<const MeshType> mesh =
- std::dynamic_pointer_cast<const MeshType>(p_mesh);
-
- writer.write(mesh,
- OutputNamedItemValueSet{
- NamedItemValue{"cell_number", mesh->connectivity().cellNumber()}},
- time, true);
- break;
- }
- case 3: {
- using MeshType = Mesh<Connectivity<3>>;
- const std::shared_ptr<const MeshType> mesh =
- std::dynamic_pointer_cast<const MeshType>(p_mesh);
-
- writer.write(mesh,
- OutputNamedItemValueSet{
- NamedItemValue{"cell_number", mesh->connectivity().cellNumber()}},
- time, true);
- break;
- }
- }
-
- time++;
- }
-
- ));
-}
diff --git a/src/language/modules/VTKModule.hpp b/src/language/modules/VTKModule.hpp
deleted file mode 100644
index 1578fad811b7dae1181053317e891dbb865ad035..0000000000000000000000000000000000000000
--- a/src/language/modules/VTKModule.hpp
+++ /dev/null
@@ -1,21 +0,0 @@
-#ifndef VTK_MODULE_HPP
-#define VTK_MODULE_HPP
-
-#include <language/modules/BuiltinModule.hpp>
-#include <utils/PugsMacros.hpp>
-
-class VTKModule : public BuiltinModule
-{
- public:
- std::string_view
- name() const final
- {
- return "vtk";
- }
-
- VTKModule();
-
- ~VTKModule() = default;
-};
-
-#endif // VTK_MODULE_HPP
diff --git a/src/language/modules/WriterModule.cpp b/src/language/modules/WriterModule.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..611952939906e0afce5ea4682428ec1eefd36eaa
--- /dev/null
+++ b/src/language/modules/WriterModule.cpp
@@ -0,0 +1,91 @@
+#include <language/modules/WriterModule.hpp>
+
+#include <language/utils/BuiltinFunctionEmbedder.hpp>
+#include <language/utils/TypeDescriptor.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/GmshReader.hpp>
+#include <mesh/Mesh.hpp>
+#include <output/GnuplotWriter.hpp>
+#include <output/IWriter.hpp>
+#include <output/NamedDiscreteFunction.hpp>
+#include <output/VTKWriter.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/IDiscreteFunction.hpp>
+#include <scheme/IDiscreteFunctionDescriptor.hpp>
+
+WriterModule::WriterModule()
+{
+ this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const NamedDiscreteFunction>>);
+
+ this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const IWriter>>);
+
+ this->_addBuiltinFunction("vtk_writer",
+ std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IWriter>(const std::string&,
+ const double&)>>(
+
+ [](const std::string& filename, const double& period) {
+ return std::make_shared<VTKWriter>(filename, period);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("gnuplot_writer",
+ std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IWriter>(const std::string&,
+ const double&)>>(
+
+ [](const std::string& filename, const double& period) {
+ return std::make_shared<GnuplotWriter>(filename, period);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("name_output",
+ std::make_shared<BuiltinFunctionEmbedder<
+ std::shared_ptr<const NamedDiscreteFunction>(std::shared_ptr<const IDiscreteFunction>,
+ const std::string&)>>(
+
+ [](std::shared_ptr<const IDiscreteFunction> discrete_function,
+ const std::string& name) -> std::shared_ptr<const NamedDiscreteFunction> {
+ return std::make_shared<const NamedDiscreteFunction>(discrete_function, name);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("write_mesh",
+ std::make_shared<BuiltinFunctionEmbedder<void(std::shared_ptr<const IWriter>,
+ std::shared_ptr<const IMesh>)>>(
+
+ [](std::shared_ptr<const IWriter> writer, std::shared_ptr<const IMesh> p_mesh) -> void {
+ writer->writeMesh(p_mesh);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("write",
+ std::make_shared<BuiltinFunctionEmbedder<
+ void(std::shared_ptr<const IWriter>,
+ const std::vector<std::shared_ptr<const NamedDiscreteFunction>>&, const double&)>>(
+
+ [](std::shared_ptr<const IWriter> writer,
+ const std::vector<std::shared_ptr<const NamedDiscreteFunction>>&
+ named_discrete_function_list,
+ const double& time) -> void {
+ writer->writeIfNeeded(named_discrete_function_list, time);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("force_write",
+ std::make_shared<BuiltinFunctionEmbedder<
+ void(std::shared_ptr<const IWriter>,
+ const std::vector<std::shared_ptr<const NamedDiscreteFunction>>&, const double&)>>(
+
+ [](std::shared_ptr<const IWriter> writer,
+ const std::vector<std::shared_ptr<const NamedDiscreteFunction>>&
+ named_discrete_function_list,
+ const double& time) -> void {
+ writer->writeForced(named_discrete_function_list, time);
+ }
+
+ ));
+}
diff --git a/src/language/modules/WriterModule.hpp b/src/language/modules/WriterModule.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..14a961d7a812b34aac1c5010ddda674749fdbfb8
--- /dev/null
+++ b/src/language/modules/WriterModule.hpp
@@ -0,0 +1,37 @@
+#ifndef WRITER_MODULE_HPP
+#define WRITER_MODULE_HPP
+
+#include <language/modules/BuiltinModule.hpp>
+#include <language/utils/ASTNodeDataTypeTraits.hpp>
+#include <utils/PugsMacros.hpp>
+
+class OutputNamedItemValueSet;
+class NamedDiscreteFunction;
+class IDiscreteFunction;
+
+#include <string>
+
+template <>
+inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const NamedDiscreteFunction>> =
+ ASTNodeDataType::build<ASTNodeDataType::type_id_t>("output");
+
+class IWriter;
+template <>
+inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const IWriter>> =
+ ASTNodeDataType::build<ASTNodeDataType::type_id_t>("writer");
+
+class WriterModule : public BuiltinModule
+{
+ public:
+ std::string_view
+ name() const final
+ {
+ return "writer";
+ }
+
+ WriterModule();
+
+ ~WriterModule() = default;
+};
+
+#endif // WRITER_MODULE_HPP
diff --git a/src/language/node_processor/AffectationProcessor.hpp b/src/language/node_processor/AffectationProcessor.hpp
index 959578619b023bc28236b1ae52141b018aad8c69..07a52de44e09cba57cf5b60aba65e844ff7faa94 100644
--- a/src/language/node_processor/AffectationProcessor.hpp
+++ b/src/language/node_processor/AffectationProcessor.hpp
@@ -121,7 +121,7 @@ class AffectationExecutor final : public IAffectationExecutor
}
} else {
if constexpr (std::is_same_v<OperatorT, language::eq_op>) {
- if constexpr (std::is_convertible_v<ValueT, DataT>) {
+ if constexpr (std::is_convertible_v<DataT, ValueT>) {
m_lhs = std::get<DataT>(rhs);
} else if constexpr (std::is_same_v<DataT, AggregateDataVariant>) {
const AggregateDataVariant& v = std::get<AggregateDataVariant>(rhs);
@@ -222,7 +222,7 @@ class MatrixComponentAffectationExecutor final : public IAffectationExecutor
}()};
public:
- MatrixComponentAffectationExecutor(ASTNode& node,
+ MatrixComponentAffectationExecutor(ASTNode& lhs_node,
ArrayT& lhs_array,
ASTNode& index0_expression,
ASTNode& index1_expression)
@@ -230,7 +230,7 @@ class MatrixComponentAffectationExecutor final : public IAffectationExecutor
{
// LCOV_EXCL_START
if constexpr (not m_is_defined) {
- throw ParseError("unexpected error: invalid operands to affectation expression", std::vector{node.begin()});
+ throw ParseError("unexpected error: invalid operands to affectation expression", std::vector{lhs_node.begin()});
}
// LCOV_EXCL_STOP
}
@@ -275,11 +275,8 @@ class MatrixComponentAffectationExecutor final : public IAffectationExecutor
}
} else {
if constexpr (std::is_same_v<OperatorT, language::eq_op>) {
- if constexpr (std::is_same_v<ValueT, DataT>) {
- m_lhs_array(index0_value, index1_value) = std::get<DataT>(rhs);
- } else {
- m_lhs_array(index0_value, index1_value) = static_cast<ValueT>(std::get<DataT>(rhs));
- }
+ static_assert(std::is_convertible_v<DataT, ValueT>, "unexpected types");
+ m_lhs_array(index0_value, index1_value) = static_cast<ValueT>(std::get<DataT>(rhs));
} else {
AffOp<OperatorT>().eval(m_lhs_array(index0_value, index1_value), std::get<DataT>(rhs));
}
@@ -354,11 +351,8 @@ class VectorComponentAffectationExecutor final : public IAffectationExecutor
}
} else {
if constexpr (std::is_same_v<OperatorT, language::eq_op>) {
- if constexpr (std::is_same_v<ValueT, DataT>) {
- m_lhs_array[index_value] = std::get<DataT>(rhs);
- } else {
- m_lhs_array[index_value] = static_cast<ValueT>(std::get<DataT>(rhs));
- }
+ static_assert(std::is_convertible_v<DataT, ValueT>, "incompatible data and value types");
+ m_lhs_array[index_value] = static_cast<ValueT>(std::get<DataT>(rhs));
} else {
AffOp<OperatorT>().eval(m_lhs_array[index_value], std::get<DataT>(rhs));
}
@@ -371,24 +365,16 @@ template <typename OperatorT, typename ValueT, typename DataT>
class AffectationProcessor final : public INodeProcessor
{
private:
- ASTNode& m_node;
+ ASTNode& m_rhs_node;
std::unique_ptr<IAffectationExecutor> m_affectation_executor;
- public:
- DataVariant
- execute(ExecutionPolicy& exec_policy)
+ std::unique_ptr<IAffectationExecutor>
+ _buildAffectationExecutor(ASTNode& lhs_node)
{
- m_affectation_executor->affect(exec_policy, m_node.children[1]->execute(exec_policy));
-
- return {};
- }
-
- AffectationProcessor(ASTNode& node) : m_node{node}
- {
- if (node.children[0]->is_type<language::name>()) {
- const std::string& symbol = m_node.children[0]->string();
- auto [i_symbol, found] = m_node.m_symbol_table->find(symbol, m_node.children[0]->begin());
+ if (lhs_node.is_type<language::name>()) {
+ const std::string& symbol = lhs_node.string();
+ auto [i_symbol, found] = lhs_node.m_symbol_table->find(symbol, lhs_node.begin());
Assert(found);
DataVariant& value = i_symbol->attributes().value();
@@ -397,23 +383,21 @@ class AffectationProcessor final : public INodeProcessor
}
using AffectationExecutorT = AffectationExecutor<OperatorT, ValueT, DataT>;
- m_affectation_executor = std::make_unique<AffectationExecutorT>(m_node, std::get<ValueT>(value));
- } else if (node.children[0]->is_type<language::subscript_expression>()) {
- auto& array_subscript_expression = *node.children[0];
-
- auto& array_expression = *array_subscript_expression.children[0];
+ return std::make_unique<AffectationExecutorT>(lhs_node, std::get<ValueT>(value));
+ } else if (lhs_node.is_type<language::subscript_expression>()) {
+ auto& array_expression = *lhs_node.children[0];
Assert(array_expression.is_type<language::name>());
const std::string& symbol = array_expression.string();
- auto [i_symbol, found] = m_node.m_symbol_table->find(symbol, array_subscript_expression.begin());
+ auto [i_symbol, found] = lhs_node.m_symbol_table->find(symbol, lhs_node.begin());
Assert(found);
DataVariant& value = i_symbol->attributes().value();
if (array_expression.m_data_type == ASTNodeDataType::vector_t) {
- Assert(array_subscript_expression.children.size() == 2);
+ Assert(lhs_node.children.size() == 2);
- auto& index_expression = *array_subscript_expression.children[1];
+ auto& index_expression = *lhs_node.children[1];
switch (array_expression.m_data_type.dimension()) {
case 1: {
@@ -422,9 +406,8 @@ class AffectationProcessor final : public INodeProcessor
value = ArrayTypeT{};
}
using AffectationExecutorT = VectorComponentAffectationExecutor<OperatorT, ArrayTypeT, ValueT, DataT>;
- m_affectation_executor =
- std::make_unique<AffectationExecutorT>(node, std::get<ArrayTypeT>(value), index_expression);
- break;
+
+ return std::make_unique<AffectationExecutorT>(lhs_node, std::get<ArrayTypeT>(value), index_expression);
}
case 2: {
using ArrayTypeT = TinyVector<2>;
@@ -432,9 +415,8 @@ class AffectationProcessor final : public INodeProcessor
value = ArrayTypeT{};
}
using AffectationExecutorT = VectorComponentAffectationExecutor<OperatorT, ArrayTypeT, ValueT, DataT>;
- m_affectation_executor =
- std::make_unique<AffectationExecutorT>(node, std::get<ArrayTypeT>(value), index_expression);
- break;
+
+ return std::make_unique<AffectationExecutorT>(lhs_node, std::get<ArrayTypeT>(value), index_expression);
}
case 3: {
using ArrayTypeT = TinyVector<3>;
@@ -442,23 +424,21 @@ class AffectationProcessor final : public INodeProcessor
value = ArrayTypeT{};
}
using AffectationExecutorT = VectorComponentAffectationExecutor<OperatorT, ArrayTypeT, ValueT, DataT>;
- m_affectation_executor =
- std::make_unique<AffectationExecutorT>(node, std::get<ArrayTypeT>(value), index_expression);
- break;
+
+ return std::make_unique<AffectationExecutorT>(lhs_node, std::get<ArrayTypeT>(value), index_expression);
}
// LCOV_EXCL_START
default: {
- throw ParseError("unexpected error: invalid vector dimension",
- std::vector{array_subscript_expression.begin()});
+ throw ParseError("unexpected error: invalid vector dimension", std::vector{lhs_node.begin()});
}
// LCOV_EXCL_STOP
}
} else if (array_expression.m_data_type == ASTNodeDataType::matrix_t) {
- Assert(array_subscript_expression.children.size() == 3);
+ Assert(lhs_node.children.size() == 3);
Assert(array_expression.m_data_type.nbRows() == array_expression.m_data_type.nbColumns());
- auto& index0_expression = *array_subscript_expression.children[1];
- auto& index1_expression = *array_subscript_expression.children[2];
+ auto& index0_expression = *lhs_node.children[1];
+ auto& index1_expression = *lhs_node.children[2];
switch (array_expression.m_data_type.nbRows()) {
case 1: {
@@ -467,9 +447,9 @@ class AffectationProcessor final : public INodeProcessor
value = ArrayTypeT{};
}
using AffectationExecutorT = MatrixComponentAffectationExecutor<OperatorT, ArrayTypeT, ValueT, DataT>;
- m_affectation_executor = std::make_unique<AffectationExecutorT>(node, std::get<ArrayTypeT>(value),
- index0_expression, index1_expression);
- break;
+
+ return std::make_unique<AffectationExecutorT>(lhs_node, std::get<ArrayTypeT>(value), index0_expression,
+ index1_expression);
}
case 2: {
using ArrayTypeT = TinyMatrix<2>;
@@ -477,9 +457,9 @@ class AffectationProcessor final : public INodeProcessor
value = ArrayTypeT{};
}
using AffectationExecutorT = MatrixComponentAffectationExecutor<OperatorT, ArrayTypeT, ValueT, DataT>;
- m_affectation_executor = std::make_unique<AffectationExecutorT>(node, std::get<ArrayTypeT>(value),
- index0_expression, index1_expression);
- break;
+
+ return std::make_unique<AffectationExecutorT>(lhs_node, std::get<ArrayTypeT>(value), index0_expression,
+ index1_expression);
}
case 3: {
using ArrayTypeT = TinyMatrix<3>;
@@ -487,14 +467,13 @@ class AffectationProcessor final : public INodeProcessor
value = ArrayTypeT{};
}
using AffectationExecutorT = MatrixComponentAffectationExecutor<OperatorT, ArrayTypeT, ValueT, DataT>;
- m_affectation_executor = std::make_unique<AffectationExecutorT>(node, std::get<ArrayTypeT>(value),
- index0_expression, index1_expression);
- break;
+
+ return std::make_unique<AffectationExecutorT>(lhs_node, std::get<ArrayTypeT>(value), index0_expression,
+ index1_expression);
}
// LCOV_EXCL_START
default: {
- throw ParseError("unexpected error: invalid vector dimension",
- std::vector{array_subscript_expression.begin()});
+ throw ParseError("unexpected error: invalid vector dimension", std::vector{lhs_node.begin()});
}
// LCOV_EXCL_STOP
}
@@ -505,25 +484,54 @@ class AffectationProcessor final : public INodeProcessor
}
} else {
// LCOV_EXCL_START
- throw ParseError("unexpected error: invalid lhs", std::vector{node.children[0]->begin()});
+ throw ParseError("unexpected error: invalid lhs", std::vector{lhs_node.begin()});
// LCOV_EXCL_STOP
}
}
+
+ public:
+ DataVariant
+ execute(ExecutionPolicy& exec_policy)
+ {
+ m_affectation_executor->affect(exec_policy, m_rhs_node.execute(exec_policy));
+
+ return {};
+ }
+
+ AffectationProcessor(ASTNode& lhs_node, ASTNode& rhs_node)
+ : m_rhs_node{rhs_node}, m_affectation_executor{this->_buildAffectationExecutor(lhs_node)}
+ {}
+};
+
+class AffectationToDataVariantProcessorBase : public INodeProcessor
+{
+ protected:
+ DataVariant* m_lhs;
+
+ public:
+ AffectationToDataVariantProcessorBase(ASTNode& lhs_node)
+ {
+ const std::string& symbol = lhs_node.string();
+ auto [i_symbol, found] = lhs_node.m_symbol_table->find(symbol, lhs_node.begin());
+ Assert(found);
+
+ m_lhs = &i_symbol->attributes().value();
+ }
+
+ virtual ~AffectationToDataVariantProcessorBase() = default;
};
template <typename OperatorT, typename ValueT>
-class AffectationToTinyVectorFromListProcessor final : public INodeProcessor
+class AffectationToTinyVectorFromListProcessor final : public AffectationToDataVariantProcessorBase
{
private:
- ASTNode& m_node;
-
- DataVariant* m_lhs;
+ ASTNode& m_rhs_node;
public:
DataVariant
execute(ExecutionPolicy& exec_policy)
{
- AggregateDataVariant children_values = std::get<AggregateDataVariant>(m_node.children[1]->execute(exec_policy));
+ AggregateDataVariant children_values = std::get<AggregateDataVariant>(m_rhs_node.execute(exec_policy));
static_assert(std::is_same_v<OperatorT, language::eq_op>, "forbidden affection operator for list to vectors");
@@ -537,7 +545,7 @@ class AffectationToTinyVectorFromListProcessor final : public INodeProcessor
v[i] = child_value;
} else {
// LCOV_EXCL_START
- throw ParseError("unexpected error: unexpected right hand side type in affectation", m_node.begin());
+ throw ParseError("unexpected error: unexpected right hand side type in affectation", m_rhs_node.begin());
// LCOV_EXCL_STOP
}
},
@@ -548,29 +556,22 @@ class AffectationToTinyVectorFromListProcessor final : public INodeProcessor
return {};
}
- AffectationToTinyVectorFromListProcessor(ASTNode& node) : m_node{node}
- {
- const std::string& symbol = m_node.children[0]->string();
- auto [i_symbol, found] = m_node.m_symbol_table->find(symbol, m_node.children[0]->begin());
- Assert(found);
-
- m_lhs = &i_symbol->attributes().value();
- }
+ AffectationToTinyVectorFromListProcessor(ASTNode& lhs_node, ASTNode& rhs_node)
+ : AffectationToDataVariantProcessorBase(lhs_node), m_rhs_node{rhs_node}
+ {}
};
template <typename OperatorT, typename ValueT>
-class AffectationToTinyMatrixFromListProcessor final : public INodeProcessor
+class AffectationToTinyMatrixFromListProcessor final : public AffectationToDataVariantProcessorBase
{
private:
- ASTNode& m_node;
-
- DataVariant* m_lhs;
+ ASTNode& m_rhs_node;
public:
DataVariant
execute(ExecutionPolicy& exec_policy)
{
- AggregateDataVariant children_values = std::get<AggregateDataVariant>(m_node.children[1]->execute(exec_policy));
+ AggregateDataVariant children_values = std::get<AggregateDataVariant>(m_rhs_node.execute(exec_policy));
static_assert(std::is_same_v<OperatorT, language::eq_op>, "forbidden affection operator for list to vectors");
@@ -585,7 +586,7 @@ class AffectationToTinyMatrixFromListProcessor final : public INodeProcessor
v(i, j) = child_value;
} else {
// LCOV_EXCL_START
- throw ParseError("unexpected error: unexpected right hand side type in affectation", m_node.begin());
+ throw ParseError("unexpected error: unexpected right hand side type in affectation", m_rhs_node.begin());
// LCOV_EXCL_STOP
}
},
@@ -597,29 +598,22 @@ class AffectationToTinyMatrixFromListProcessor final : public INodeProcessor
return {};
}
- AffectationToTinyMatrixFromListProcessor(ASTNode& node) : m_node{node}
- {
- const std::string& symbol = m_node.children[0]->string();
- auto [i_symbol, found] = m_node.m_symbol_table->find(symbol, m_node.children[0]->begin());
- Assert(found);
-
- m_lhs = &i_symbol->attributes().value();
- }
+ AffectationToTinyMatrixFromListProcessor(ASTNode& lhs_node, ASTNode& rhs_node)
+ : AffectationToDataVariantProcessorBase(lhs_node), m_rhs_node{rhs_node}
+ {}
};
template <typename ValueT>
-class AffectationToTupleProcessor final : public INodeProcessor
+class AffectationToTupleProcessor final : public AffectationToDataVariantProcessorBase
{
private:
- ASTNode& m_node;
-
- DataVariant* m_lhs;
+ ASTNode& m_rhs_node;
public:
DataVariant
execute(ExecutionPolicy& exec_policy)
{
- DataVariant value = m_node.children[1]->execute(exec_policy);
+ DataVariant value = m_rhs_node.execute(exec_policy);
std::visit(
[&](auto&& v) {
@@ -646,12 +640,12 @@ class AffectationToTupleProcessor final : public INodeProcessor
*m_lhs = std::vector<ValueT>{ValueT{zero}};
} else {
// LCOV_EXCL_START
- throw ParseError("unexpected error: unexpected right hand side type in affectation", m_node.begin());
+ throw ParseError("unexpected error: unexpected right hand side type in affectation", m_rhs_node.begin());
// LCOV_EXCL_STOP
}
} else {
// LCOV_EXCL_START
- throw ParseError("unexpected error: unexpected right hand side type in affectation", m_node.begin());
+ throw ParseError("unexpected error: unexpected right hand side type in affectation", m_rhs_node.begin());
// LCOV_EXCL_STOP
}
},
@@ -660,23 +654,16 @@ class AffectationToTupleProcessor final : public INodeProcessor
return {};
}
- AffectationToTupleProcessor(ASTNode& node) : m_node{node}
- {
- const std::string& symbol = m_node.children[0]->string();
- auto [i_symbol, found] = m_node.m_symbol_table->find(symbol, m_node.children[0]->begin());
- Assert(found);
-
- m_lhs = &i_symbol->attributes().value();
- }
+ AffectationToTupleProcessor(ASTNode& lhs_node, ASTNode& rhs_node)
+ : AffectationToDataVariantProcessorBase(lhs_node), m_rhs_node{rhs_node}
+ {}
};
template <typename ValueT>
-class AffectationToTupleFromListProcessor final : public INodeProcessor
+class AffectationToTupleFromListProcessor final : public AffectationToDataVariantProcessorBase
{
private:
- ASTNode& m_node;
-
- DataVariant* m_lhs;
+ ASTNode& m_rhs_node;
void
_copyAggregateDataVariant(const AggregateDataVariant& children_values)
@@ -711,7 +698,7 @@ class AffectationToTupleFromListProcessor final : public INodeProcessor
} else {
// LCOV_EXCL_START
throw ParseError("unexpected error: unexpected right hand side type in affectation",
- m_node.children[1]->children[i]->begin());
+ m_rhs_node.children[i]->begin());
// LCOV_EXCL_STOP
}
},
@@ -728,7 +715,7 @@ class AffectationToTupleFromListProcessor final : public INodeProcessor
} else {
// LCOV_EXCL_START
throw ParseError("unexpected error: unexpected right hand side type in affectation",
- m_node.children[1]->children[i]->begin());
+ m_rhs_node.children[i]->begin());
// LCOV_EXCL_STOP
}
} else if constexpr (is_tiny_matrix_v<ValueT>) {
@@ -745,7 +732,7 @@ class AffectationToTupleFromListProcessor final : public INodeProcessor
} else {
// LCOV_EXCL_START
throw ParseError("unexpected error: unexpected right hand side type in affectation",
- m_node.children[1]->children[i]->begin());
+ m_rhs_node.children[i]->begin());
// LCOV_EXCL_STOP
}
},
@@ -763,13 +750,13 @@ class AffectationToTupleFromListProcessor final : public INodeProcessor
} else {
// LCOV_EXCL_START
throw ParseError("unexpected error: unexpected right hand side type in affectation",
- m_node.children[1]->children[i]->begin());
+ m_rhs_node.children[i]->begin());
// LCOV_EXCL_STOP
}
} else {
// LCOV_EXCL_START
throw ParseError("unexpected error: unexpected right hand side type in affectation",
- m_node.children[1]->children[i]->begin());
+ m_rhs_node.children[i]->begin());
// LCOV_EXCL_STOP
}
},
@@ -805,8 +792,7 @@ class AffectationToTupleFromListProcessor final : public INodeProcessor
}
} else {
// LCOV_EXCL_START
- throw ParseError("unexpected error: unexpected right hand side type in tuple affectation",
- m_node.children[1]->begin());
+ throw ParseError("unexpected error: unexpected right hand side type in tuple affectation", m_rhs_node.begin());
// LCOV_EXCL_STOP
}
@@ -826,34 +812,23 @@ class AffectationToTupleFromListProcessor final : public INodeProcessor
this->_copyVector(value_list);
} else {
// LCOV_EXCL_START
- throw ParseError("unexpected error: invalid lhs (expecting list or tuple)",
- std::vector{m_node.children[1]->begin()});
+ throw ParseError("unexpected error: invalid lhs (expecting list or tuple)", std::vector{m_rhs_node.begin()});
// LCOV_EXCL_STOP
}
},
- m_node.children[1]->execute(exec_policy));
+ m_rhs_node.execute(exec_policy));
return {};
}
- AffectationToTupleFromListProcessor(ASTNode& node) : m_node{node}
- {
- const std::string& symbol = m_node.children[0]->string();
- auto [i_symbol, found] = m_node.m_symbol_table->find(symbol, m_node.children[0]->begin());
- Assert(found);
-
- m_lhs = &i_symbol->attributes().value();
- }
+ AffectationToTupleFromListProcessor(ASTNode& lhs_node, ASTNode& rhs_node)
+ : AffectationToDataVariantProcessorBase(lhs_node), m_rhs_node{rhs_node}
+ {}
};
template <typename ValueT>
-class AffectationFromZeroProcessor final : public INodeProcessor
+class AffectationFromZeroProcessor final : public AffectationToDataVariantProcessorBase
{
- private:
- ASTNode& m_node;
-
- DataVariant* m_lhs;
-
public:
DataVariant
execute(ExecutionPolicy&)
@@ -862,14 +837,7 @@ class AffectationFromZeroProcessor final : public INodeProcessor
return {};
}
- AffectationFromZeroProcessor(ASTNode& node) : m_node{node}
- {
- const std::string& symbol = m_node.children[0]->string();
- auto [i_symbol, found] = m_node.m_symbol_table->find(symbol, m_node.children[0]->begin());
- Assert(found);
-
- m_lhs = &i_symbol->attributes().value();
- }
+ AffectationFromZeroProcessor(ASTNode& lhs_node) : AffectationToDataVariantProcessorBase(lhs_node) {}
};
template <typename OperatorT>
diff --git a/src/language/utils/AffectationProcessorBuilder.hpp b/src/language/utils/AffectationProcessorBuilder.hpp
index 8d11599a01c67566e0f88fcb34cad581bc76de45..0e4b7ace65654b650989e817dbe2e4a545ae320f 100644
--- a/src/language/utils/AffectationProcessorBuilder.hpp
+++ b/src/language/utils/AffectationProcessorBuilder.hpp
@@ -15,18 +15,18 @@ class AffectationProcessorBuilder final : public IAffectationProcessorBuilder
public:
AffectationProcessorBuilder() = default;
std::unique_ptr<INodeProcessor>
- getNodeProcessor(ASTNode& node) const
+ getNodeProcessor(ASTNode& lhs_node, ASTNode& rhs_node) const final
{
if constexpr (std::is_same_v<ValueT, TinyVector<1>> and std::is_same_v<DataT, int64_t> and
std::is_same_v<OperatorT, language::eq_op>) {
// Special treatment for the case 0 -> R^1
- if ((node.children[1]->is_type<language::integer>()) and (std::stoi(node.children[1]->string()) == 0)) {
- return std::make_unique<AffectationFromZeroProcessor<ValueT>>(node);
+ if ((rhs_node.is_type<language::integer>()) and (std::stoi(rhs_node.string()) == 0)) {
+ return std::make_unique<AffectationFromZeroProcessor<ValueT>>(lhs_node);
} else {
- return std::make_unique<AffectationProcessor<OperatorT, ValueT, DataT>>(node);
+ return std::make_unique<AffectationProcessor<OperatorT, ValueT, DataT>>(lhs_node, rhs_node);
}
} else {
- return std::make_unique<AffectationProcessor<OperatorT, ValueT, DataT>>(node);
+ return std::make_unique<AffectationProcessor<OperatorT, ValueT, DataT>>(lhs_node, rhs_node);
}
}
};
@@ -37,9 +37,9 @@ class AffectationToTupleProcessorBuilder final : public IAffectationProcessorBui
public:
AffectationToTupleProcessorBuilder() = default;
std::unique_ptr<INodeProcessor>
- getNodeProcessor(ASTNode& node) const
+ getNodeProcessor(ASTNode& lhs_node, ASTNode& rhs_node) const final
{
- return std::make_unique<AffectationToTupleProcessor<ValueT>>(node);
+ return std::make_unique<AffectationToTupleProcessor<ValueT>>(lhs_node, rhs_node);
}
};
@@ -49,10 +49,10 @@ class AffectationToTupleFromListProcessorBuilder final : public IAffectationProc
public:
AffectationToTupleFromListProcessorBuilder() = default;
std::unique_ptr<INodeProcessor>
- getNodeProcessor(ASTNode& node) const
+ getNodeProcessor(ASTNode& lhs_node, ASTNode& rhs_node) const final
{
- ASTNodeNaturalConversionChecker(*node.children[1], node.children[0]->m_data_type);
- return std::make_unique<AffectationToTupleFromListProcessor<ValueT>>(node);
+ ASTNodeNaturalConversionChecker(rhs_node, lhs_node.m_data_type);
+ return std::make_unique<AffectationToTupleFromListProcessor<ValueT>>(lhs_node, rhs_node);
}
};
@@ -62,9 +62,9 @@ class AffectationToTinyVectorFromListProcessorBuilder final : public IAffectatio
public:
AffectationToTinyVectorFromListProcessorBuilder() = default;
std::unique_ptr<INodeProcessor>
- getNodeProcessor(ASTNode& node) const
+ getNodeProcessor(ASTNode& lhs_node, ASTNode& rhs_node) const final
{
- return std::make_unique<AffectationToTinyVectorFromListProcessor<OperatorT, ValueT>>(node);
+ return std::make_unique<AffectationToTinyVectorFromListProcessor<OperatorT, ValueT>>(lhs_node, rhs_node);
}
};
@@ -74,9 +74,9 @@ class AffectationToTinyMatrixFromListProcessorBuilder final : public IAffectatio
public:
AffectationToTinyMatrixFromListProcessorBuilder() = default;
std::unique_ptr<INodeProcessor>
- getNodeProcessor(ASTNode& node) const
+ getNodeProcessor(ASTNode& lhs_node, ASTNode& rhs_node) const final
{
- return std::make_unique<AffectationToTinyMatrixFromListProcessor<OperatorT, ValueT>>(node);
+ return std::make_unique<AffectationToTinyMatrixFromListProcessor<OperatorT, ValueT>>(lhs_node, rhs_node);
}
};
@@ -86,12 +86,12 @@ class AffectationFromZeroProcessorBuilder final : public IAffectationProcessorBu
public:
AffectationFromZeroProcessorBuilder() = default;
std::unique_ptr<INodeProcessor>
- getNodeProcessor(ASTNode& node) const
+ getNodeProcessor(ASTNode& lhs_node, ASTNode& rhs_node) const final
{
- if (std::stoi(node.children[1]->string()) == 0) {
- return std::make_unique<AffectationFromZeroProcessor<ValueT>>(node);
+ if (std::stoi(rhs_node.string()) == 0) {
+ return std::make_unique<AffectationFromZeroProcessor<ValueT>>(lhs_node);
} else {
- throw ParseError("invalid integral value (0 is the solely valid value)", std::vector{node.children[1]->begin()});
+ throw ParseError("invalid integral value (0 is the solely valid value)", std::vector{lhs_node.begin()});
}
}
};
diff --git a/src/language/utils/BuiltinFunctionEmbedder.hpp b/src/language/utils/BuiltinFunctionEmbedder.hpp
index afb991acadcb36236bcd55a34c7840ca399ab6f0..265379a368ef3b7f847cd8de2b763a70b525f8a4 100644
--- a/src/language/utils/BuiltinFunctionEmbedder.hpp
+++ b/src/language/utils/BuiltinFunctionEmbedder.hpp
@@ -130,23 +130,30 @@ class BuiltinFunctionEmbedder<FX(Args...)> : public IBuiltinFunctionEmbedder
return ast_node_data_type_from<T>;
}
- template <size_t I>
+ template <typename TupleT, size_t I>
PUGS_INLINE ASTNodeDataType
- _getOneParameterDataType(ArgsTuple& t) const
+ _getOneParameterDataType() const
{
- using ArgN_T = std::decay_t<decltype(std::get<I>(t))>;
+ using ArgN_T = std::decay_t<decltype(std::get<I>(TupleT{}))>;
return _getDataType<ArgN_T>();
}
template <size_t... I>
- PUGS_INLINE std::vector<ASTNodeDataType>
- _getParameterDataTypes(ArgsTuple t, std::index_sequence<I...>) const
+ PUGS_INLINE std::vector<ASTNodeDataType> _getParameterDataTypes(std::index_sequence<I...>) const
{
std::vector<ASTNodeDataType> parameter_type_list;
- (parameter_type_list.push_back(this->_getOneParameterDataType<I>(t)), ...);
+ (parameter_type_list.push_back(this->_getOneParameterDataType<ArgsTuple, I>()), ...);
return parameter_type_list;
}
+ template <size_t... I>
+ PUGS_INLINE std::vector<std::shared_ptr<const ASTNodeDataType>> _getCompoundDataTypes(std::index_sequence<I...>) const
+ {
+ std::vector<std::shared_ptr<const ASTNodeDataType>> compound_type_list;
+ (compound_type_list.push_back(std::make_shared<ASTNodeDataType>(this->_getOneParameterDataType<FX, I>())), ...);
+ return compound_type_list;
+ }
+
template <typename T>
PUGS_INLINE std::shared_ptr<IDataHandler>
_createHandler(std::shared_ptr<T> data) const
@@ -175,21 +182,52 @@ class BuiltinFunctionEmbedder<FX(Args...)> : public IBuiltinFunctionEmbedder
(_check_arg<I>(), ...);
}
+ template <typename ResultT>
+ PUGS_INLINE DataVariant
+ _resultToDataVariant(ResultT&& result) const
+ {
+ if constexpr (is_data_variant_v<std::decay_t<ResultT>>) {
+ return std::move(result);
+ } else {
+ return EmbeddedData(_createHandler(std::move(result)));
+ }
+ }
+
+ PUGS_INLINE
+ AggregateDataVariant
+ _applyToAggregate(const ArgsTuple& t) const
+ {
+ auto tuple_result = std::apply(m_f, t);
+ std::vector<DataVariant> vector_result;
+ vector_result.reserve(std::tuple_size_v<decltype(tuple_result)>);
+
+ std::apply([&](auto&&... result) { ((vector_result.emplace_back(_resultToDataVariant(result))), ...); },
+ tuple_result);
+
+ return vector_result;
+ }
+
public:
PUGS_INLINE ASTNodeDataType
getReturnDataType() const final
{
- return _getDataType<FX>();
+ if constexpr (is_std_tuple_v<FX>) {
+ constexpr size_t N = std::tuple_size_v<FX>;
+ using IndexSequence = std::make_index_sequence<N>;
+
+ return ASTNodeDataType::build<ASTNodeDataType::list_t>(this->_getCompoundDataTypes(IndexSequence{}));
+ } else {
+ return this->_getDataType<FX>();
+ }
}
PUGS_INLINE std::vector<ASTNodeDataType>
getParameterDataTypes() const final
{
- constexpr size_t N = std::tuple_size_v<ArgsTuple>;
- ArgsTuple t;
+ constexpr size_t N = std::tuple_size_v<ArgsTuple>;
using IndexSequence = std::make_index_sequence<N>;
- return this->_getParameterDataTypes(t, IndexSequence{});
+ return this->_getParameterDataTypes(IndexSequence{});
}
PUGS_INLINE size_t
@@ -198,7 +236,6 @@ class BuiltinFunctionEmbedder<FX(Args...)> : public IBuiltinFunctionEmbedder
return sizeof...(Args);
}
- public:
PUGS_INLINE
DataVariant
apply(const std::vector<DataVariant>& x) const final
@@ -210,6 +247,8 @@ class BuiltinFunctionEmbedder<FX(Args...)> : public IBuiltinFunctionEmbedder
this->_copyFromVector(t, x, IndexSequence{});
if constexpr (is_data_variant_v<FX>) {
return {std::apply(m_f, t)};
+ } else if constexpr (is_std_tuple_v<FX>) {
+ return this->_applyToAggregate(t);
} else if constexpr (std::is_same_v<FX, void>) {
std::apply(m_f, t);
return {};
@@ -245,6 +284,22 @@ class BuiltinFunctionEmbedder<FX(void)> : public IBuiltinFunctionEmbedder
private:
std::function<FX(void)> m_f;
+ template <typename TupleT, size_t I>
+ PUGS_INLINE ASTNodeDataType
+ _getOneParameterDataType() const
+ {
+ using ArgN_T = std::decay_t<decltype(std::get<I>(TupleT{}))>;
+ return _getDataType<ArgN_T>();
+ }
+
+ template <size_t... I>
+ PUGS_INLINE std::vector<std::shared_ptr<const ASTNodeDataType>> _getCompoundDataTypes(std::index_sequence<I...>) const
+ {
+ std::vector<std::shared_ptr<const ASTNodeDataType>> compound_type_list;
+ (compound_type_list.push_back(std::make_shared<ASTNodeDataType>(this->_getOneParameterDataType<FX, I>())), ...);
+ return compound_type_list;
+ }
+
template <typename T>
PUGS_INLINE std::shared_ptr<IDataHandler>
_createHandler(std::shared_ptr<T> data) const
@@ -260,11 +315,42 @@ class BuiltinFunctionEmbedder<FX(void)> : public IBuiltinFunctionEmbedder
return ast_node_data_type_from<T>;
}
+ template <typename ResultT>
+ PUGS_INLINE DataVariant
+ _resultToDataVariant(ResultT&& result) const
+ {
+ if constexpr (is_data_variant_v<std::decay_t<ResultT>>) {
+ return std::move(result);
+ } else {
+ return EmbeddedData(_createHandler(std::move(result)));
+ }
+ }
+
+ PUGS_INLINE
+ AggregateDataVariant
+ _applyToAggregate() const
+ {
+ auto tuple_result = m_f();
+ std::vector<DataVariant> vector_result;
+ vector_result.reserve(std::tuple_size_v<decltype(tuple_result)>);
+
+ std::apply([&](auto&&... result) { ((vector_result.emplace_back(_resultToDataVariant(result))), ...); },
+ tuple_result);
+
+ return vector_result;
+ }
+
public:
PUGS_INLINE ASTNodeDataType
getReturnDataType() const final
{
- return this->_getDataType<FX>();
+ if constexpr (is_std_tuple_v<FX>) {
+ constexpr size_t N = std::tuple_size_v<FX>;
+ using IndexSequence = std::make_index_sequence<N>;
+ return ASTNodeDataType::build<ASTNodeDataType::list_t>(this->_getCompoundDataTypes(IndexSequence{}));
+ } else {
+ return this->_getDataType<FX>();
+ }
}
PUGS_INLINE std::vector<ASTNodeDataType>
@@ -285,6 +371,8 @@ class BuiltinFunctionEmbedder<FX(void)> : public IBuiltinFunctionEmbedder
{
if constexpr (is_data_variant_v<FX>) {
return {m_f()};
+ } else if constexpr (is_std_tuple_v<FX>) {
+ return this->_applyToAggregate();
} else if constexpr (std::is_same_v<FX, void>) {
m_f();
return {};
diff --git a/src/language/utils/CMakeLists.txt b/src/language/utils/CMakeLists.txt
index 22c0ae52f184a7f18d852057cbe9b0616f1905e3..32a75c66f0dd73cc10c09ca2c6aecf640962c45c 100644
--- a/src/language/utils/CMakeLists.txt
+++ b/src/language/utils/CMakeLists.txt
@@ -22,6 +22,7 @@ add_library(PugsLanguageUtils
BinaryOperatorRegisterForZ.cpp
DataVariant.cpp
EmbeddedData.cpp
+ FunctionSymbolId.cpp
IncDecOperatorRegisterForN.cpp
IncDecOperatorRegisterForR.cpp
IncDecOperatorRegisterForZ.cpp
diff --git a/src/language/utils/FunctionSymbolId.cpp b/src/language/utils/FunctionSymbolId.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d45c35c2d3e71311128026c8779094e370eee22d
--- /dev/null
+++ b/src/language/utils/FunctionSymbolId.cpp
@@ -0,0 +1,8 @@
+#include <language/utils/FunctionSymbolId.hpp>
+#include <language/utils/SymbolTable.hpp>
+
+const FunctionDescriptor&
+FunctionSymbolId::descriptor() const
+{
+ return m_symbol_table->functionTable()[m_function_id];
+}
diff --git a/src/language/utils/FunctionSymbolId.hpp b/src/language/utils/FunctionSymbolId.hpp
index f96b783fa5bf3853507f4c650e02972faa5d6ebb..21496cb0e56e49f9c5a5d58d407a927f3d32cfd3 100644
--- a/src/language/utils/FunctionSymbolId.hpp
+++ b/src/language/utils/FunctionSymbolId.hpp
@@ -8,7 +8,9 @@
#include <iostream>
#include <memory>
+class FunctionDescriptor;
class SymbolTable;
+
class FunctionSymbolId
{
private:
@@ -29,6 +31,8 @@ class FunctionSymbolId
return *m_symbol_table;
}
+ [[nodiscard]] const FunctionDescriptor& descriptor() const;
+
friend std::ostream&
operator<<(std::ostream& os, const FunctionSymbolId& function_symbol_id)
{
diff --git a/src/language/utils/IAffectationProcessorBuilder.hpp b/src/language/utils/IAffectationProcessorBuilder.hpp
index 6d5f8540df9435bba4b86e6bee22f980b5727d7b..cc0ca85f67ece81b6c1313a13c429dfeb9718cd7 100644
--- a/src/language/utils/IAffectationProcessorBuilder.hpp
+++ b/src/language/utils/IAffectationProcessorBuilder.hpp
@@ -9,7 +9,7 @@ class INodeProcessor;
class IAffectationProcessorBuilder
{
public:
- virtual std::unique_ptr<INodeProcessor> getNodeProcessor(ASTNode& node) const = 0;
+ virtual std::unique_ptr<INodeProcessor> getNodeProcessor(ASTNode& lhs_node, ASTNode& rhs_node) const = 0;
virtual ~IAffectationProcessorBuilder() = default;
};
diff --git a/src/language/utils/PugsFunctionAdapter.hpp b/src/language/utils/PugsFunctionAdapter.hpp
index 9b57a908e6341291d14e5e8f9c80aa8a87275dcc..b6d3a259dbb676d5f1259c058e0ea1f36fb925fd 100644
--- a/src/language/utils/PugsFunctionAdapter.hpp
+++ b/src/language/utils/PugsFunctionAdapter.hpp
@@ -119,10 +119,10 @@ class PugsFunctionAdapter<OutputType(InputType...)>
[[nodiscard]] PUGS_INLINE static auto&
getFunctionExpression(const FunctionSymbolId& function_symbol_id)
{
- auto& function = function_symbol_id.symbolTable().functionTable()[function_symbol_id.id()];
- _checkFunction(function);
+ auto& function_descriptor = function_symbol_id.descriptor();
+ _checkFunction(function_descriptor);
- return *function.definitionNode().children[1];
+ return *function_descriptor.definitionNode().children[1];
}
[[nodiscard]] PUGS_INLINE static auto
@@ -209,7 +209,7 @@ class PugsFunctionAdapter<OutputType(InputType...)>
};
} else {
// If this point is reached must be a 0 vector
- return [](DataVariant &&) -> OutputType { return OutputType{ZeroType{}}; };
+ return [](DataVariant&&) -> OutputType { return OutputType{ZeroType{}}; };
}
}
case ASTNodeDataType::double_t: {
@@ -236,8 +236,8 @@ class PugsFunctionAdapter<OutputType(InputType...)>
AggregateDataVariant& v = std::get<AggregateDataVariant>(result);
OutputType x;
- for (size_t i = 0, l = 0; i < x.dimension(); ++i) {
- for (size_t j = 0; j < x.dimension(); ++j, ++l) {
+ for (size_t i = 0, l = 0; i < x.nbRows(); ++i) {
+ for (size_t j = 0; j < x.nbColumns(); ++j, ++l) {
std::visit(
[&](auto&& Aij) {
using Aij_T = std::decay_t<decltype(Aij)>;
@@ -288,7 +288,7 @@ class PugsFunctionAdapter<OutputType(InputType...)>
};
} else {
// If this point is reached must be a 0 matrix
- return [](DataVariant &&) -> OutputType { return OutputType{ZeroType{}}; };
+ return [](DataVariant&&) -> OutputType { return OutputType{ZeroType{}}; };
}
}
case ASTNodeDataType::double_t: {
diff --git a/src/language/utils/SymbolTable.hpp b/src/language/utils/SymbolTable.hpp
index 6aa089bd8de7aa5a3a6a6cf2a8115e77de25dbe8..772ac6564859dd402ba2086d4fed7dc0b4c64059 100644
--- a/src/language/utils/SymbolTable.hpp
+++ b/src/language/utils/SymbolTable.hpp
@@ -267,12 +267,15 @@ class SymbolTable
clearValues()
{
for (auto& symbol : m_symbol_list) {
- std::visit(
- [](auto&& value) {
- using T = std::decay_t<decltype(value)>;
- value = T{};
- },
- symbol.attributes().value());
+ // local functions must kept their values (id)
+ if (symbol.attributes().dataType() != ASTNodeDataType::function_t) {
+ std::visit(
+ [](auto&& value) {
+ using T = std::decay_t<decltype(value)>;
+ value = T{};
+ },
+ symbol.attributes().value());
+ }
}
}
diff --git a/src/output/CMakeLists.txt b/src/output/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..a3e6799657d2ed422fee5ab7018939ea080cf4c2
--- /dev/null
+++ b/src/output/CMakeLists.txt
@@ -0,0 +1,16 @@
+# ------------------- Source files --------------------
+
+add_library(
+ PugsOutput
+ GnuplotWriter.cpp
+ VTKWriter.cpp
+ WriterBase.cpp)
+
+# ------------------- Installation --------------------
+# temporary version workaround
+if(${CMAKE_VERSION} VERSION_LESS "3.13.0")
+ install(TARGETS PugsMesh
+ RUNTIME DESTINATION bin
+ LIBRARY DESTINATION lib
+ ARCHIVE DESTINATION lib)
+endif()
diff --git a/src/output/GnuplotWriter.cpp b/src/output/GnuplotWriter.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..436ae7ffa22d7f3aaa20cbd44e2331393958823c
--- /dev/null
+++ b/src/output/GnuplotWriter.cpp
@@ -0,0 +1,253 @@
+#include <output/GnuplotWriter.hpp>
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/ItemValue.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <utils/Messenger.hpp>
+#include <utils/PugsTraits.hpp>
+#include <utils/RevisionInfo.hpp>
+
+#include <utils/Demangle.hpp>
+
+#include <fstream>
+#include <iomanip>
+
+std::string
+GnuplotWriter::_getDateAndVersionComment() const
+{
+ std::ostringstream os;
+
+ std::time_t now = std::time(nullptr);
+ os << "# Generated by pugs: " << std::ctime(&now);
+ os << "# version: " << RevisionInfo::version() << '\n';
+ os << "# tag: " << RevisionInfo::gitTag() << '\n';
+ os << "# HEAD: " << RevisionInfo::gitHead() << '\n';
+ os << "# hash: " << RevisionInfo::gitHash() << " (" << ((RevisionInfo::gitIsClean()) ? "clean" : "dirty") << ")\n";
+ os << '\n';
+
+ return os.str();
+}
+
+std::string
+GnuplotWriter::_getFilename() const
+{
+ std::ostringstream sout;
+ sout << m_base_filename;
+ sout << '.' << std::setfill('0') << std::setw(4) << m_saved_times.size() << ".gnu";
+ return sout.str();
+}
+
+template <size_t Dimension>
+void
+GnuplotWriter::_writePreamble(const OutputNamedItemValueSet& output_named_item_value_set, std::ostream& fout) const
+{
+ fout << "# list of data\n";
+ fout << "# 1:x";
+ if constexpr (Dimension > 1) {
+ fout << " 2:y";
+ }
+ uint64_t i = Dimension + 1;
+ for (const auto& i_named_item_value : output_named_item_value_set) {
+ const std::string name = i_named_item_value.first;
+ const auto& item_value_variant = i_named_item_value.second;
+ std::visit(
+ [&](auto&& item_value) {
+ using ItemValueType = std::decay_t<decltype(item_value)>;
+ using DataType = std::decay_t<typename ItemValueType::data_type>;
+ if constexpr (std::is_arithmetic_v<DataType>) {
+ fout << ' ' << i++ << ':' << name;
+ } else if constexpr (is_tiny_vector_v<DataType>) {
+ for (size_t j = 0; j < DataType{}.dimension(); ++j) {
+ fout << ' ' << i++ << ':' << name << '[' << j << ']';
+ }
+ } else if constexpr (is_tiny_matrix_v<DataType>) {
+ for (size_t j = 0; j < DataType{}.nbRows(); ++j) {
+ for (size_t k = 0; k < DataType{}.nbColumns(); ++k) {
+ fout << ' ' << i++ << ':' << name << '(' << j << ',' << k << ')';
+ }
+ }
+ } else {
+ throw UnexpectedError("invalid data type");
+ }
+ },
+ item_value_variant);
+ }
+ fout << "\n\n";
+}
+
+template <typename DataType>
+void
+GnuplotWriter::_writeCellValue(const CellValue<DataType>& cell_value, CellId cell_id, std::ostream& fout) const
+{
+ const auto& value = cell_value[cell_id];
+ if constexpr (std::is_arithmetic_v<DataType>) {
+ fout << ' ' << value;
+ } else if constexpr (is_tiny_vector_v<std::decay_t<DataType>>) {
+ for (size_t i = 0; i < value.dimension(); ++i) {
+ fout << ' ' << value[i];
+ }
+ } else if constexpr (is_tiny_matrix_v<std::decay_t<DataType>>) {
+ for (size_t i = 0; i < value.nbRows(); ++i) {
+ for (size_t j = 0; j < value.nbColumns(); ++j) {
+ fout << ' ' << value(i, j);
+ }
+ }
+ } else {
+ throw UnexpectedError("invalid data type for cell value output: " + demangle<DataType>());
+ }
+}
+
+template <typename DataType, ItemType item_type>
+void
+GnuplotWriter::_writeValue(const ItemValue<DataType, item_type>& item_value,
+ CellId cell_id,
+ NodeId node_id,
+ std::ostream& fout) const
+{
+ if constexpr (item_type == ItemType::cell) {
+ this->_writeCellValue(item_value, cell_id, fout);
+ } else if constexpr (item_type == ItemType::node) {
+ this->_writeNodeValue(item_value, node_id, fout);
+ } else {
+ throw UnexpectedError{"invalid item type"};
+ }
+}
+
+template <typename MeshType>
+void
+GnuplotWriter::_writeValues(const std::shared_ptr<const MeshType>& mesh,
+ const OutputNamedItemValueSet& output_named_item_value_set,
+ std::ostream& fout) const
+{
+ if constexpr (MeshType::Dimension <= 2) {
+ auto cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
+ auto cell_is_owned = mesh->connectivity().cellIsOwned();
+
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ if (cell_is_owned[cell_id]) {
+ const auto& cell_nodes = cell_to_node_matrix[cell_id];
+ for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
+ const NodeId& node_id = cell_nodes[i_node];
+ const TinyVector<MeshType::Dimension>& xr = mesh->xr()[node_id];
+ fout << xr[0];
+ if (MeshType::Dimension == 2) {
+ fout << ' ' << xr[1];
+ }
+ for (auto [name, item_value] : output_named_item_value_set) {
+ std::visit([&](auto&& cell_value) { _writeValue(cell_value, cell_id, node_id, fout); }, item_value);
+ }
+
+ fout << '\n';
+ }
+ fout << '\n';
+ }
+ }
+ } else {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+}
+
+template <typename ValueType>
+void
+GnuplotWriter::_writeNodeValue(const NodeValue<ValueType>& node_value, NodeId node_id, std::ostream& fout) const
+{
+ const auto& value = node_value[node_id];
+ if constexpr (std::is_arithmetic_v<ValueType>) {
+ fout << ' ' << value;
+ } else if constexpr (is_tiny_vector_v<std::decay_t<ValueType>>) {
+ for (size_t i = 0; i < value.dimension(); ++i) {
+ fout << ' ' << value[i];
+ }
+ } else if constexpr (is_tiny_matrix_v<std::decay_t<ValueType>>) {
+ for (size_t i = 0; i < value.nbRows(); ++i) {
+ for (size_t j = 0; j < value.nbColumns(); ++j) {
+ fout << ' ' << value(i, j);
+ }
+ }
+ } else {
+ throw UnexpectedError("invalid data type for cell value output: " + demangle<ValueType>());
+ }
+}
+
+template <typename MeshType>
+void
+GnuplotWriter::_write(const std::shared_ptr<const MeshType>& mesh,
+ const OutputNamedItemValueSet& output_named_item_value_set,
+ double time) const
+{
+ if (parallel::rank() == 0) {
+ std::ofstream fout{_getFilename()};
+ fout.precision(15);
+ fout.setf(std::ios_base::scientific);
+
+ fout << this->_getDateAndVersionComment();
+ fout << "\n# time = " << time << "\n\n";
+
+ this->_writePreamble<MeshType::Dimension>(output_named_item_value_set, fout);
+ }
+
+ for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+ if (i_rank == parallel::rank()) {
+ std::ofstream fout{_getFilename(), std::ios_base::app};
+ fout.precision(15);
+ fout.setf(std::ios_base::scientific);
+
+ this->_writeValues(mesh, output_named_item_value_set, fout);
+ }
+ parallel::barrier();
+ }
+}
+
+void
+GnuplotWriter::writeMesh(const std::shared_ptr<const IMesh>& p_mesh) const
+{
+ OutputNamedItemValueSet output_named_item_value_set{};
+
+ std::ofstream fout(_getFilename());
+ fout.precision(15);
+ fout.setf(std::ios_base::scientific);
+ fout << _getDateAndVersionComment();
+
+ switch (p_mesh->dimension()) {
+ case 1: {
+ this->_writePreamble<1>(output_named_item_value_set, fout);
+ this->_writeValues(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(p_mesh), output_named_item_value_set,
+ fout);
+ break;
+ }
+ case 2: {
+ this->_writePreamble<2>(output_named_item_value_set, fout);
+ this->_writeValues(std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(p_mesh), output_named_item_value_set,
+ fout);
+ break;
+ }
+ default: {
+ throw NormalError("gnuplot format is not available in dimension " + std::to_string(p_mesh->dimension()));
+ }
+ }
+}
+
+void
+GnuplotWriter::write(const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
+ double time) const
+{
+ std::shared_ptr mesh = this->_getMesh(named_discrete_function_list);
+
+ OutputNamedItemValueSet output_named_item_value_set = this->_getOutputNamedItemValueSet(named_discrete_function_list);
+
+ switch (mesh->dimension()) {
+ case 1: {
+ this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(mesh), output_named_item_value_set, time);
+ break;
+ }
+ case 2: {
+ this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(mesh), output_named_item_value_set, time);
+ break;
+ }
+ default: {
+ throw NormalError("gnuplot format is not available in dimension " + std::to_string(mesh->dimension()));
+ }
+ }
+}
diff --git a/src/output/GnuplotWriter.hpp b/src/output/GnuplotWriter.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..bf6487c998c29d8a5c12177a574eea587a1fb6ea
--- /dev/null
+++ b/src/output/GnuplotWriter.hpp
@@ -0,0 +1,57 @@
+#ifndef GNUPLOT_WRITER_HPP
+#define GNUPLOT_WRITER_HPP
+
+#include <output/WriterBase.hpp>
+
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+#include <output/OutputNamedItemValueSet.hpp>
+
+class IMesh;
+
+#include <string>
+
+class GnuplotWriter : public WriterBase
+{
+ private:
+ std::string _getDateAndVersionComment() const;
+
+ std::string _getFilename() const;
+
+ template <typename DataType>
+ void _writeCellValue(const CellValue<DataType>& cell_value, CellId cell_id, std::ostream& fout) const;
+
+ template <typename DataType>
+ void _writeNodeValue(const NodeValue<DataType>& node_value, NodeId cell_id, std::ostream& fout) const;
+
+ template <size_t Dimension>
+ void _writePreamble(const OutputNamedItemValueSet& output_named_item_value_set, std::ostream& fout) const;
+
+ template <typename DataType, ItemType item_type>
+ void _writeValue(const ItemValue<DataType, item_type>& item_value,
+ CellId cell_id,
+ NodeId node_id,
+ std::ostream& fout) const;
+
+ template <typename MeshType>
+ void _writeValues(const std::shared_ptr<const MeshType>& mesh,
+ const OutputNamedItemValueSet& output_named_item_value_set,
+ std::ostream& fout) const;
+
+ template <typename MeshType>
+ void _write(const std::shared_ptr<const MeshType>& mesh,
+ const OutputNamedItemValueSet& output_named_item_value_set,
+ double time) const;
+
+ public:
+ void writeMesh(const std::shared_ptr<const IMesh>& mesh) const final;
+
+ void write(const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
+ double time) const final;
+
+ GnuplotWriter(const std::string& base_filename, const double time_period) : WriterBase(base_filename, time_period) {}
+
+ ~GnuplotWriter() = default;
+};
+
+#endif // GNUPLOT_WRITER_HPP
diff --git a/src/output/IWriter.hpp b/src/output/IWriter.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..02317ad98b055567133f9b969476fedbd2a5c0e1
--- /dev/null
+++ b/src/output/IWriter.hpp
@@ -0,0 +1,28 @@
+#ifndef I_WRITER_HPP
+#define I_WRITER_HPP
+
+#include <output/NamedDiscreteFunction.hpp>
+
+#include <memory>
+#include <vector>
+
+class IMesh;
+
+class IWriter
+{
+ public:
+ virtual void writeMesh(const std::shared_ptr<const IMesh>& mesh) const = 0;
+
+ virtual void writeIfNeeded(
+ const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
+ double time) const = 0;
+
+ virtual void writeForced(
+ const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
+ double time) const = 0;
+
+ IWriter() = default;
+ virtual ~IWriter() = default;
+};
+
+#endif // I_WRITER_HPP
diff --git a/src/output/NamedDiscreteFunction.hpp b/src/output/NamedDiscreteFunction.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..1b0850040a1fd944c5e3d3ba361f7d764bece3dd
--- /dev/null
+++ b/src/output/NamedDiscreteFunction.hpp
@@ -0,0 +1,39 @@
+#ifndef NAMED_DISCRETE_FUNCTION_HPP
+#define NAMED_DISCRETE_FUNCTION_HPP
+
+#include <memory>
+#include <string>
+
+class IDiscreteFunction;
+
+class NamedDiscreteFunction
+{
+ private:
+ std::shared_ptr<const IDiscreteFunction> m_discrete_function;
+ std::string m_name;
+
+ public:
+ const std::string&
+ name() const
+ {
+ return m_name;
+ }
+
+ const std::shared_ptr<const IDiscreteFunction>
+ discreteFunction() const
+ {
+ return m_discrete_function;
+ }
+
+ NamedDiscreteFunction(const std::shared_ptr<const IDiscreteFunction>& discrete_function, const std::string& name)
+ : m_discrete_function{discrete_function}, m_name{name}
+ {}
+
+ NamedDiscreteFunction(const NamedDiscreteFunction&) = default;
+ NamedDiscreteFunction(NamedDiscreteFunction&&) = default;
+
+ NamedDiscreteFunction() = default;
+ ~NamedDiscreteFunction() = default;
+};
+
+#endif // NAMED_DISCRETE_FUNCTION_HPP
diff --git a/src/output/OutputNamedItemValueSet.hpp b/src/output/OutputNamedItemValueSet.hpp
index 7dcc23e3d0aadaa95cb3a4333e38b73d1edad9a5..c7464135e3531c36ba3b1d3a801088cba002fbd7 100644
--- a/src/output/OutputNamedItemValueSet.hpp
+++ b/src/output/OutputNamedItemValueSet.hpp
@@ -54,19 +54,29 @@ class NamedItemValue
class OutputNamedItemValueSet
{
public:
- using ItemValueVariant = std::variant<NodeValue<const int>,
+ using ItemValueVariant = std::variant<NodeValue<const bool>,
+ NodeValue<const int>,
NodeValue<const long int>,
+ NodeValue<const unsigned long int>,
NodeValue<const double>,
NodeValue<const TinyVector<1, double>>,
NodeValue<const TinyVector<2, double>>,
NodeValue<const TinyVector<3, double>>,
+ NodeValue<const TinyMatrix<1, double>>,
+ NodeValue<const TinyMatrix<2, double>>,
+ NodeValue<const TinyMatrix<3, double>>,
+ CellValue<const bool>,
CellValue<const int>,
CellValue<const long int>,
+ CellValue<const unsigned long int>,
CellValue<const double>,
CellValue<const TinyVector<1, double>>,
CellValue<const TinyVector<2, double>>,
- CellValue<const TinyVector<3, double>>>;
+ CellValue<const TinyVector<3, double>>,
+ CellValue<const TinyMatrix<1, double>>,
+ CellValue<const TinyMatrix<2, double>>,
+ CellValue<const TinyMatrix<3, double>>>;
private:
std::map<std::string, ItemValueVariant> m_name_itemvariant_map;
@@ -103,6 +113,13 @@ class OutputNamedItemValueSet
return m_name_itemvariant_map.end();
}
+ template <typename DataType, ItemType item_type>
+ void
+ add(const NamedItemValue<DataType, item_type>& named_itemvalue)
+ {
+ _doInsert(named_itemvalue);
+ }
+
template <typename... DataType, ItemType... item_type>
OutputNamedItemValueSet(NamedItemValue<DataType, item_type>... named_itemvalue)
{
diff --git a/src/output/VTKWriter.cpp b/src/output/VTKWriter.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..48f81dd376bd96c36820a535b2ff2fa20a4e831b
--- /dev/null
+++ b/src/output/VTKWriter.cpp
@@ -0,0 +1,585 @@
+#include <output/VTKWriter.hpp>
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
+#include <utils/Messenger.hpp>
+#include <utils/RevisionInfo.hpp>
+
+#include <ctime>
+#include <fstream>
+#include <iomanip>
+#include <sstream>
+#include <unordered_map>
+
+std::string
+VTKWriter::_getFilenamePVTU() const
+{
+ std::ostringstream sout;
+ sout << m_base_filename;
+ sout << '.' << std::setfill('0') << std::setw(4) << m_saved_times.size() << ".pvtu";
+ return sout.str();
+}
+
+std::string
+VTKWriter::_getDateAndVersionComment() const
+{
+ std::ostringstream os;
+
+ std::time_t now = std::time(nullptr);
+ os << "<!--\n";
+ os << " Generated by pugs: " << std::ctime(&now);
+ os << " version: " << RevisionInfo::version() << '\n';
+ os << " tag: " << RevisionInfo::gitTag() << '\n';
+ os << " HEAD: " << RevisionInfo::gitHead() << '\n';
+ os << " hash: " << RevisionInfo::gitHash() << " (" << ((RevisionInfo::gitIsClean()) ? "clean" : "dirty") << ")\n";
+ os << "-->\n";
+ return os.str();
+}
+
+std::string
+VTKWriter::_getFilenameVTU(int rank_number) const
+{
+ std::ostringstream sout;
+ sout << m_base_filename;
+ if (parallel::size() > 1) {
+ sout << '-' << std::setfill('0') << std::setw(4) << rank_number;
+ }
+ sout << '.' << std::setfill('0') << std::setw(4) << m_saved_times.size() << ".vtu";
+ return sout.str();
+}
+
+template <typename DataType>
+void
+VTKWriter::_write_node_pvtu(std::ofstream& os, const std::string& name, const NodeValue<const DataType>&) const
+{
+ os << "<PDataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\"/>\n";
+}
+
+template <size_t N, typename DataType>
+void
+VTKWriter::_write_node_pvtu(std::ofstream& os,
+ const std::string& name,
+ const NodeValue<const TinyVector<N, DataType>>&) const
+{
+ os << "<PDataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\"" << N
+ << "\"/>\n";
+}
+
+template <size_t N, typename DataType>
+void
+VTKWriter::_write_node_pvtu(std::ofstream& os,
+ const std::string& name,
+ const NodeValue<const TinyMatrix<N, DataType>>&) const
+{
+ os << "<PDataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\"" << N * N
+ << "\"/>\n";
+}
+
+template <typename DataType>
+void
+VTKWriter::_write_node_pvtu(std::ofstream&, const std::string&, const CellValue<const DataType>&) const
+{}
+
+template <typename DataType>
+void
+VTKWriter::_write_cell_pvtu(std::ofstream& os, const std::string& name, const CellValue<const DataType>&) const
+{
+ os << "<PDataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\"/>\n";
+}
+
+template <size_t N, typename DataType>
+void
+VTKWriter::_write_cell_pvtu(std::ofstream& os,
+ const std::string& name,
+ const CellValue<const TinyVector<N, DataType>>&) const
+{
+ os << "<PDataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\"" << N
+ << "\"/>\n";
+}
+
+template <size_t N, typename DataType>
+void
+VTKWriter::_write_cell_pvtu(std::ofstream& os,
+ const std::string& name,
+ const CellValue<const TinyMatrix<N, DataType>>&) const
+{
+ os << "<PDataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\"" << N * N
+ << "\"/>\n";
+}
+
+template <typename DataType>
+void
+VTKWriter::_write_cell_pvtu(std::ofstream&, const std::string&, const NodeValue<const DataType>&) const
+{}
+
+template <typename DataType>
+struct VTKWriter::VTKType
+{
+ inline const static std::string name = [] {
+ static_assert(std::is_arithmetic_v<DataType>, "invalid data type");
+
+ if constexpr (std::is_integral_v<DataType>) {
+ if constexpr (std::is_unsigned_v<DataType>) {
+ return "UInt" + std::to_string(sizeof(DataType) * 8);
+ } else {
+ return "UInt" + std::to_string(sizeof(DataType) * 8);
+ }
+ } else if constexpr (std::is_floating_point_v<DataType>) {
+ return "Float" + std::to_string(sizeof(DataType) * 8);
+ }
+ }();
+};
+
+template <typename DataType>
+void
+VTKWriter::_write_array(std::ofstream& os, const std::string& name, const Array<DataType>& item_value) const
+{
+ os << "<DataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\">\n";
+ for (typename Array<DataType>::index_type i = 0; i < item_value.size(); ++i) {
+ // The following '+' enforces integer output for char types
+ os << +item_value[i] << ' ';
+ }
+ os << "\n</DataArray>\n";
+}
+
+template <size_t N, typename DataType>
+void
+VTKWriter::_write_array(std::ofstream& os,
+ const std::string& name,
+ const Array<TinyVector<N, DataType>>& item_value) const
+{
+ os << "<DataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\"" << N
+ << "\">\n";
+ for (typename Array<DataType>::index_type i = 0; i < item_value.size(); ++i) {
+ for (size_t j = 0; j < N; ++j) {
+ // The following '+' enforces integer output for char types
+ os << +item_value[i][j] << ' ';
+ }
+ }
+ os << "\n</DataArray>\n";
+}
+
+template <typename DataType>
+void
+VTKWriter::_write_node_value(std::ofstream& os,
+ const std::string& name,
+ const NodeValue<const DataType>& item_value) const
+{
+ os << "<DataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\">\n";
+ for (NodeId i = 0; i < item_value.size(); ++i) {
+ // The following '+' enforces integer output for char types
+ os << +item_value[i] << ' ';
+ }
+ os << "\n</DataArray>\n";
+}
+
+template <size_t N, typename DataType>
+void
+VTKWriter::_write_node_value(std::ofstream& os,
+ const std::string& name,
+ const NodeValue<const TinyVector<N, DataType>>& item_value) const
+{
+ os << "<DataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\"" << N
+ << "\">\n";
+ for (NodeId i = 0; i < item_value.size(); ++i) {
+ for (size_t j = 0; j < N; ++j) {
+ // The following '+' enforces integer output for char types
+ os << +item_value[i][j] << ' ';
+ }
+ }
+ os << "\n</DataArray>\n";
+}
+
+template <size_t N, typename DataType>
+void
+VTKWriter::_write_node_value(std::ofstream& os,
+ const std::string& name,
+ const NodeValue<const TinyMatrix<N, DataType>>& item_value) const
+{
+ os << "<DataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\"" << N * N
+ << "\">\n";
+ for (NodeId i = 0; i < item_value.size(); ++i) {
+ for (size_t j = 0; j < N; ++j) {
+ for (size_t k = 0; k < N; ++k) {
+ // The following '+' enforces integer output for char types
+ os << +item_value[i](j, k) << ' ';
+ }
+ }
+ }
+ os << "\n</DataArray>\n";
+}
+
+template <typename DataType>
+void
+VTKWriter::_write_node_value(std::ofstream&, const std::string&, const CellValue<const DataType>&) const
+{}
+
+template <typename DataType>
+void
+VTKWriter::_write_cell_value(std::ofstream& os,
+ const std::string& name,
+ const CellValue<const DataType>& item_value) const
+{
+ os << "<DataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\">\n";
+ for (CellId i = 0; i < item_value.size(); ++i) {
+ // The following '+' enforces integer output for char types
+ os << +item_value[i] << ' ';
+ }
+ os << "\n</DataArray>\n";
+}
+
+template <size_t N, typename DataType>
+void
+VTKWriter::_write_cell_value(std::ofstream& os,
+ const std::string& name,
+ const CellValue<const TinyVector<N, DataType>>& item_value) const
+{
+ os << "<DataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\"" << N
+ << "\">\n";
+ for (CellId i = 0; i < item_value.size(); ++i) {
+ for (size_t j = 0; j < N; ++j) {
+ // The following '+' enforces integer output for char types
+ os << +item_value[i][j] << ' ';
+ }
+ }
+ os << "\n</DataArray>\n";
+}
+
+template <size_t N, typename DataType>
+void
+VTKWriter::_write_cell_value(std::ofstream& os,
+ const std::string& name,
+ const CellValue<const TinyMatrix<N, DataType>>& item_value) const
+{
+ os << "<DataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\"" << N * N
+ << "\">\n";
+ for (CellId i = 0; i < item_value.size(); ++i) {
+ for (size_t j = 0; j < N; ++j) {
+ for (size_t k = 0; k < N; ++k) {
+ // The following '+' enforces integer output for char types
+ os << +item_value[i](j, k) << ' ';
+ }
+ }
+ }
+ os << "\n</DataArray>\n";
+}
+
+template <typename DataType>
+void
+VTKWriter::_write_cell_value(std::ofstream&, const std::string&, const NodeValue<const DataType>&) const
+{}
+
+template <typename MeshType>
+void
+VTKWriter::_write(const std::shared_ptr<const MeshType>& mesh,
+ const OutputNamedItemValueSet& given_output_named_item_value_set,
+ double time) const
+{
+ OutputNamedItemValueSet output_named_item_value_set{given_output_named_item_value_set};
+ // Adding basic mesh information
+ output_named_item_value_set.add(NamedItemValue{"cell_number", mesh->connectivity().cellNumber()});
+ output_named_item_value_set.add(NamedItemValue{"node_number", mesh->connectivity().nodeNumber()});
+
+ if (parallel::rank() == 0) { // write PVTK file
+ std::ofstream fout(_getFilenamePVTU());
+ fout << "<?xml version=\"1.0\"?>\n";
+ fout << _getDateAndVersionComment();
+ fout << "<VTKFile type=\"PUnstructuredGrid\">\n";
+ fout << "<PUnstructuredGrid GhostLevel=\"0\">\n";
+
+ fout << "<PPoints>\n";
+ fout << "<PDataArray Name=\"Positions\" NumberOfComponents=\"3\" "
+ "type=\"Float64\"/>\n";
+ fout << "</PPoints>\n";
+
+ fout << "<PCells>\n";
+ fout << "<PDataArray type=\"Int32\" Name=\"connectivity\" "
+ "NumberOfComponents=\"1\"/>\n";
+ fout << "<PDataArray type=\"UInt32\" Name=\"offsets\" "
+ "NumberOfComponents=\"1\"/>\n";
+ fout << "<PDataArray type=\"Int8\" Name=\"types\" "
+ "NumberOfComponents=\"1\"/>\n";
+ for (const auto& [name, item_value_variant] : output_named_item_value_set) {
+ std::visit([&, name = name](auto&& item_value) { return this->_write_cell_pvtu(fout, name, item_value); },
+ item_value_variant);
+ }
+ if constexpr (MeshType::Dimension == 3) {
+ fout << "<PDataArray type=\"Int64\" IdType=\"1\" Name=\"faces\"/>\n";
+ fout << "<PDataArray type=\"Int64\" IdType=\"1\" Name=\"faceoffsets\"/>\n";
+ }
+ fout << "</PCells>\n";
+
+ fout << "<PPointData>\n";
+ for (const auto& [name, item_value_variant] : output_named_item_value_set) {
+ std::visit([&, name = name](auto&& item_value) { return this->_write_node_pvtu(fout, name, item_value); },
+ item_value_variant);
+ }
+ fout << "</PPointData>\n";
+
+ fout << "<PCellData>\n";
+ for (const auto& [name, item_value_variant] : output_named_item_value_set) {
+ std::visit([&, name = name](auto&& item_value) { return this->_write_cell_pvtu(fout, name, item_value); },
+ item_value_variant);
+ }
+ fout << "</PCellData>\n";
+
+ for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+ fout << "<Piece Source=\"" << _getFilenameVTU(i_rank) << "\"/>\n";
+ }
+ fout << "</PUnstructuredGrid>\n";
+ fout << "</VTKFile>\n";
+ }
+
+ { // write VTK files
+ std::ofstream fout(_getFilenameVTU(parallel::rank()));
+ fout << "<?xml version=\"1.0\"?>\n";
+ fout << _getDateAndVersionComment();
+ fout << "<VTKFile type=\"UnstructuredGrid\">\n";
+ fout << "<UnstructuredGrid>\n";
+ fout << "<Piece NumberOfPoints=\"" << mesh->numberOfNodes() << "\" NumberOfCells=\"" << mesh->numberOfCells()
+ << "\">\n";
+ fout << "<CellData>\n";
+ for (const auto& [name, item_value_variant] : output_named_item_value_set) {
+ std::visit([&, name = name](auto&& item_value) { return this->_write_cell_value(fout, name, item_value); },
+ item_value_variant);
+ }
+ fout << "</CellData>\n";
+ fout << "<PointData>\n";
+ for (const auto& [name, item_value_variant] : output_named_item_value_set) {
+ std::visit([&, name = name](auto&& item_value) { return this->_write_node_value(fout, name, item_value); },
+ item_value_variant);
+ }
+ fout << "</PointData>\n";
+ fout << "<Points>\n";
+ {
+ using Rd = TinyVector<MeshType::Dimension>;
+ const NodeValue<const Rd>& xr = mesh->xr();
+ Array<TinyVector<3>> positions(mesh->numberOfNodes());
+ parallel_for(
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+ for (unsigned short i = 0; i < MeshType::Dimension; ++i) {
+ positions[r][i] = xr[r][i];
+ }
+ for (unsigned short i = MeshType::Dimension; i < 3; ++i) {
+ positions[r][i] = 0;
+ }
+ });
+ _write_array(fout, "Positions", positions);
+ }
+ fout << "</Points>\n";
+
+ fout << "<Cells>\n";
+ {
+ const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
+
+ _write_array(fout, "connectivity", cell_to_node_matrix.entries());
+ }
+
+ {
+ const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
+ Array<unsigned int> offsets(mesh->numberOfCells());
+ unsigned int offset = 0;
+ for (CellId j = 0; j < mesh->numberOfCells(); ++j) {
+ const auto& cell_nodes = cell_to_node_matrix[j];
+ offset += cell_nodes.size();
+ offsets[j] = offset;
+ }
+ _write_array(fout, "offsets", offsets);
+ }
+
+ {
+ Array<int8_t> types(mesh->numberOfCells());
+ const auto& cell_type = mesh->connectivity().cellType();
+ const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
+
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ switch (cell_type[j]) {
+ case CellType::Line: {
+ types[j] = 3;
+ break;
+ }
+ case CellType::Triangle: {
+ types[j] = 5;
+ break;
+ }
+ case CellType::Quadrangle: {
+ types[j] = 9;
+ break;
+ }
+ case CellType::Tetrahedron: {
+ types[j] = 10;
+ break;
+ }
+ case CellType::Pyramid: {
+ if (cell_to_node_matrix[j].size() == 5) {
+ types[j] = 14; // quadrangle basis
+ } else {
+ types[j] = 41; // polygonal basis
+ }
+ break;
+ }
+ case CellType::Prism: {
+ types[j] = 13;
+ break;
+ }
+ case CellType::Hexahedron: {
+ types[j] = 12;
+ break;
+ }
+ case CellType::Diamond: {
+ types[j] = 41;
+ break;
+ }
+ default: {
+ std::ostringstream os;
+ os << __FILE__ << ':' << __LINE__ << ": unknown cell type";
+ throw UnexpectedError(os.str());
+ }
+ }
+ });
+ _write_array(fout, "types", types);
+ if constexpr (MeshType::Dimension == 3) {
+ const bool has_general_polyhedron = [&] {
+ for (size_t i = 0; i < types.size(); ++i) {
+ if (types[i] == 41)
+ return true;
+ }
+ return false;
+ }();
+ if (has_general_polyhedron) {
+ const auto& cell_to_face_matrix = mesh->connectivity().cellToFaceMatrix();
+ const auto& face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
+ const auto& cell_face_is_reversed = mesh->connectivity().cellFaceIsReversed();
+
+ Array<size_t> faces_offsets(mesh->numberOfCells());
+ size_t next_offset = 0;
+ fout << "<DataArray type=\"Int64\" IdType=\"1\" Name=\"faces\">\n";
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ const auto& cell_nodes = cell_to_node_matrix[cell_id];
+ std::unordered_map<size_t, size_t> node_id_to_local_number_in_cell;
+ for (size_t i_cell_node = 0; i_cell_node < cell_nodes.size(); ++i_cell_node) {
+ node_id_to_local_number_in_cell[cell_nodes[i_cell_node]] = i_cell_node;
+ }
+
+ const auto& cell_faces = cell_to_face_matrix[cell_id];
+ fout << cell_faces.size() << '\n';
+ next_offset++;
+ for (size_t i_cell_face = 0; i_cell_face < cell_faces.size(); ++i_cell_face) {
+ const FaceId& face_id = cell_faces[i_cell_face];
+ const auto& face_nodes = face_to_node_matrix[face_id];
+ fout << face_nodes.size();
+ next_offset++;
+ Array<size_t> face_node_in_cell(face_nodes.size());
+
+ for (size_t i_face_node = 0; i_face_node < face_nodes.size(); ++i_face_node) {
+ const NodeId& node_id = face_nodes[i_face_node];
+ auto i_node_id_to_local_number_in_cell = node_id_to_local_number_in_cell.find(node_id);
+ Assert(i_node_id_to_local_number_in_cell != node_id_to_local_number_in_cell.end());
+ face_node_in_cell[i_face_node] = i_node_id_to_local_number_in_cell->second;
+ }
+
+ if (cell_face_is_reversed(cell_id, i_cell_face)) {
+ for (size_t i = 0; i < face_node_in_cell.size(); ++i) {
+ fout << ' ' << face_node_in_cell[face_node_in_cell.size() - 1 - i];
+ }
+ } else {
+ for (size_t i = 0; i < face_node_in_cell.size(); ++i) {
+ fout << ' ' << face_node_in_cell[i];
+ }
+ }
+
+ next_offset += face_nodes.size();
+
+ fout << '\n';
+ }
+ faces_offsets[cell_id] = next_offset;
+ }
+ fout << "</DataArray>\n";
+ fout << "<DataArray type=\"Int64\" IdType=\"1\" Name=\"faceoffsets\">\n";
+ for (size_t i_face_offsets = 0; i_face_offsets < faces_offsets.size(); ++i_face_offsets) {
+ fout << faces_offsets[i_face_offsets] << '\n';
+ }
+ fout << "</DataArray>\n";
+ }
+ }
+ }
+
+ fout << "</Cells>\n";
+ fout << "</Piece>\n";
+ fout << "</UnstructuredGrid>\n";
+ fout << "</VTKFile>\n";
+ }
+
+ if (parallel::rank() == 0) { // write PVD file
+ std::ofstream fout(m_base_filename + ".pvd");
+
+ fout << "<?xml version=\"1.0\"?>\n";
+ fout << "<VTKFile type=\"Collection\" version=\"0.1\">\n";
+
+ fout << "<Collection>\n";
+ for (size_t i_time = 0; i_time < m_saved_times.size(); ++i_time) {
+ std::ostringstream sout;
+ sout << m_base_filename;
+ sout << '.' << std::setfill('0') << std::setw(4) << i_time << ".pvtu";
+
+ fout << "<DataSet timestep=\"" << m_saved_times[i_time] << "\" "
+ << "file=\"" << sout.str() << "\"/>\n";
+ }
+ fout << "<DataSet timestep=\"" << time << "\" group=\"\" part=\"0\" file=\"" << _getFilenamePVTU() << "\"/>\n";
+
+ fout << "</Collection>\n";
+ fout << "</VTKFile>\n";
+ }
+}
+
+void
+VTKWriter::writeMesh(const std::shared_ptr<const IMesh>& mesh) const
+{
+ OutputNamedItemValueSet output_named_item_value_set;
+
+ switch (mesh->dimension()) {
+ case 1: {
+ this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(mesh), output_named_item_value_set, 0);
+ break;
+ }
+ case 2: {
+ this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(mesh), output_named_item_value_set, 0);
+ break;
+ }
+ case 3: {
+ this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<3>>>(mesh), output_named_item_value_set, 0);
+ break;
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ }
+}
+
+void
+VTKWriter::write(const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
+ double time) const
+{
+ std::shared_ptr mesh = this->_getMesh(named_discrete_function_list);
+
+ OutputNamedItemValueSet output_named_item_value_set = this->_getOutputNamedItemValueSet(named_discrete_function_list);
+
+ switch (mesh->dimension()) {
+ case 1: {
+ this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(mesh), output_named_item_value_set, time);
+ break;
+ }
+ case 2: {
+ this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(mesh), output_named_item_value_set, time);
+ break;
+ }
+ case 3: {
+ this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<3>>>(mesh), output_named_item_value_set, time);
+ break;
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ }
+}
diff --git a/src/output/VTKWriter.hpp b/src/output/VTKWriter.hpp
index f3c202ffc51633d291d36150d662fe950023ef58..704ba63f990ead2b6b8e8f64cafdda627c2be64b 100644
--- a/src/output/VTKWriter.hpp
+++ b/src/output/VTKWriter.hpp
@@ -1,456 +1,110 @@
#ifndef VTK_WRITER_HPP
#define VTK_WRITER_HPP
+#include <output/WriterBase.hpp>
+
+#include <algebra/TinyMatrix.hpp>
#include <algebra/TinyVector.hpp>
-#include <mesh/CellType.hpp>
-#include <mesh/IConnectivity.hpp>
-#include <mesh/ItemValue.hpp>
#include <output/OutputNamedItemValueSet.hpp>
-#include <utils/Exceptions.hpp>
-#include <utils/Messenger.hpp>
-#include <utils/PugsAssert.hpp>
-
-#include <fstream>
-#include <iomanip>
-#include <iostream>
-#include <sstream>
+
+class IMesh;
+
#include <string>
-#include <unordered_map>
-class VTKWriter
+class VTKWriter : public WriterBase
{
private:
- const std::string m_base_filename;
- unsigned int m_file_number;
- double m_last_time;
- const double m_time_period;
-
- std::string
- _getFilenamePVTU()
- {
- std::ostringstream sout;
- sout << m_base_filename;
- sout << '.' << std::setfill('0') << std::setw(4) << m_file_number << ".pvtu";
- return sout.str();
- }
- std::string
- _getFilenameVTU(int rank_number) const
- {
- std::ostringstream sout;
- sout << m_base_filename;
- if (parallel::size() > 1) {
- sout << '-' << std::setfill('0') << std::setw(4) << rank_number;
- }
- sout << '.' << std::setfill('0') << std::setw(4) << m_file_number << ".vtu";
- return sout.str();
- }
+ std::string _getDateAndVersionComment() const;
+
+ std::string _getFilenamePVTU() const;
+
+ std::string _getFilenameVTU(int rank_number) const;
template <typename DataType>
- void
- _write_node_pvtu(std::ofstream& os, const std::string& name, const NodeValue<const DataType>&)
- {
- os << "<PDataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\"/>\n";
- }
+ void _write_node_pvtu(std::ofstream& os, const std::string& name, const NodeValue<const DataType>&) const;
template <size_t N, typename DataType>
- void
- _write_node_pvtu(std::ofstream& os, const std::string& name, const NodeValue<const TinyVector<N, DataType>>&)
- {
- os << "<PDataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\"" << N
- << "\"/>\n";
- }
+ void _write_node_pvtu(std::ofstream& os,
+ const std::string& name,
+ const NodeValue<const TinyVector<N, DataType>>&) const;
+
+ template <size_t N, typename DataType>
+ void _write_node_pvtu(std::ofstream& os,
+ const std::string& name,
+ const NodeValue<const TinyMatrix<N, DataType>>&) const;
template <typename DataType>
- void
- _write_node_pvtu(std::ofstream&, const std::string&, const CellValue<const DataType>&)
- {}
+ void _write_node_pvtu(std::ofstream&, const std::string&, const CellValue<const DataType>&) const;
template <typename DataType>
- void
- _write_cell_pvtu(std::ofstream& os, const std::string& name, const CellValue<const DataType>&)
- {
- os << "<PDataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\"/>\n";
- }
+ void _write_cell_pvtu(std::ofstream& os, const std::string& name, const CellValue<const DataType>&) const;
+
+ template <size_t N, typename DataType>
+ void _write_cell_pvtu(std::ofstream& os,
+ const std::string& name,
+ const CellValue<const TinyVector<N, DataType>>&) const;
template <size_t N, typename DataType>
- void
- _write_cell_pvtu(std::ofstream& os, const std::string& name, const CellValue<const TinyVector<N, DataType>>&)
- {
- os << "<PDataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\"" << N
- << "\"/>\n";
- }
+ void _write_cell_pvtu(std::ofstream& os,
+ const std::string& name,
+ const CellValue<const TinyMatrix<N, DataType>>&) const;
template <typename DataType>
- void
- _write_cell_pvtu(std::ofstream&, const std::string&, const NodeValue<const DataType>&)
- {}
+ void _write_cell_pvtu(std::ofstream&, const std::string&, const NodeValue<const DataType>&) const;
template <typename DataType>
- struct VTKType
- {
- inline const static std::string name = [] {
- static_assert(std::is_arithmetic_v<DataType>, "invalid data type");
-
- if constexpr (std::is_integral_v<DataType>) {
- if constexpr (std::is_unsigned_v<DataType>) {
- return "UInt" + std::to_string(sizeof(DataType) * 8);
- } else {
- return "UInt" + std::to_string(sizeof(DataType) * 8);
- }
- } else if constexpr (std::is_floating_point_v<DataType>) {
- return "Float" + std::to_string(sizeof(DataType) * 8);
- }
- }();
- };
+ struct VTKType;
template <typename DataType>
- void
- _write_array(std::ofstream& os, const std::string& name, const Array<DataType>& item_value)
- {
- os << "<DataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\">\n";
- for (typename Array<DataType>::index_type i = 0; i < item_value.size(); ++i) {
- // The following '+' enforces integer output for char types
- os << +item_value[i] << ' ';
- }
- os << "\n</DataArray>\n";
- }
+ void _write_array(std::ofstream& os, const std::string& name, const Array<DataType>& item_value) const;
template <size_t N, typename DataType>
- void
- _write_array(std::ofstream& os, const std::string& name, const Array<TinyVector<N, DataType>>& item_value)
- {
- os << "<DataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\"" << N
- << "\">\n";
- for (typename Array<DataType>::index_type i = 0; i < item_value.size(); ++i) {
- for (size_t j = 0; j < N; ++j) {
- // The following '+' enforces integer output for char types
- os << +item_value[i][j] << ' ';
- }
- }
- os << "\n</DataArray>\n";
- }
+ void _write_array(std::ofstream& os, const std::string& name, const Array<TinyVector<N, DataType>>& item_value) const;
template <typename DataType>
- void
- _write_node_value(std::ofstream& os, const std::string& name, const NodeValue<const DataType>& item_value)
- {
- os << "<DataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\">\n";
- for (NodeId i = 0; i < item_value.size(); ++i) {
- // The following '+' enforces integer output for char types
- os << +item_value[i] << ' ';
- }
- os << "\n</DataArray>\n";
- }
+ void _write_node_value(std::ofstream& os, const std::string& name, const NodeValue<const DataType>& item_value) const;
template <size_t N, typename DataType>
- void
- _write_node_value(std::ofstream& os,
- const std::string& name,
- const NodeValue<const TinyVector<N, DataType>>& item_value)
- {
- os << "<DataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\"" << N
- << "\">\n";
- for (NodeId i = 0; i < item_value.size(); ++i) {
- for (size_t j = 0; j < N; ++j) {
- // The following '+' enforces integer output for char types
- os << +item_value[i][j] << ' ';
- }
- }
- os << "\n</DataArray>\n";
- }
+ void _write_node_value(std::ofstream& os,
+ const std::string& name,
+ const NodeValue<const TinyVector<N, DataType>>& item_value) const;
+
+ template <size_t N, typename DataType>
+ void _write_node_value(std::ofstream& os,
+ const std::string& name,
+ const NodeValue<const TinyMatrix<N, DataType>>& item_value) const;
template <typename DataType>
- void
- _write_node_value(std::ofstream&, const std::string&, const CellValue<const DataType>&)
- {}
+ void _write_node_value(std::ofstream&, const std::string&, const CellValue<const DataType>&) const;
template <typename DataType>
- void
- _write_cell_value(std::ofstream& os, const std::string& name, const CellValue<const DataType>& item_value)
- {
- os << "<DataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\">\n";
- for (CellId i = 0; i < item_value.size(); ++i) {
- // The following '+' enforces integer output for char types
- os << +item_value[i] << ' ';
- }
- os << "\n</DataArray>\n";
- }
+ void _write_cell_value(std::ofstream& os, const std::string& name, const CellValue<const DataType>& item_value) const;
template <size_t N, typename DataType>
- void
- _write_cell_value(std::ofstream& os,
- const std::string& name,
- const CellValue<const TinyVector<N, DataType>>& item_value)
- {
- os << "<DataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\"" << N
- << "\">\n";
- for (CellId i = 0; i < item_value.size(); ++i) {
- for (size_t j = 0; j < N; ++j) {
- // The following '+' enforces integer output for char types
- os << +item_value[i][j] << ' ';
- }
- }
- os << "\n</DataArray>\n";
- }
+ void _write_cell_value(std::ofstream& os,
+ const std::string& name,
+ const CellValue<const TinyVector<N, DataType>>& item_value) const;
+
+ template <size_t N, typename DataType>
+ void _write_cell_value(std::ofstream& os,
+ const std::string& name,
+ const CellValue<const TinyMatrix<N, DataType>>& item_value) const;
template <typename DataType>
- void
- _write_cell_value(std::ofstream&, const std::string&, const NodeValue<const DataType>&)
- {}
+ void _write_cell_value(std::ofstream&, const std::string&, const NodeValue<const DataType>&) const;
- public:
template <typename MeshType>
- void
- write(const std::shared_ptr<const MeshType>& mesh,
- const OutputNamedItemValueSet& output_named_item_value_set,
- double time,
- bool forced_output = false)
- {
- if (time == m_last_time)
- return; // output already performed
- if ((time - m_last_time >= m_time_period) or forced_output) {
- m_last_time = time;
- } else {
- return;
- }
-
- if (parallel::rank() == 0) { // write PVTK file
- std::ofstream fout(_getFilenamePVTU());
- fout << "<?xml version=\"1.0\"?>\n";
- fout << "<VTKFile type=\"PUnstructuredGrid\">\n";
- fout << "<PUnstructuredGrid GhostLevel=\"0\">\n";
-
- fout << "<PPoints>\n";
- fout << "<PDataArray Name=\"Positions\" NumberOfComponents=\"3\" "
- "type=\"Float64\"/>\n";
- fout << "</PPoints>\n";
-
- fout << "<PCells>\n";
- fout << "<PDataArray type=\"Int32\" Name=\"connectivity\" "
- "NumberOfComponents=\"1\"/>\n";
- fout << "<PDataArray type=\"UInt32\" Name=\"offsets\" "
- "NumberOfComponents=\"1\"/>\n";
- fout << "<PDataArray type=\"Int8\" Name=\"types\" "
- "NumberOfComponents=\"1\"/>\n";
- for (const auto& [name, item_value_variant] : output_named_item_value_set) {
- std::visit([&, name = name](auto&& item_value) { return this->_write_cell_pvtu(fout, name, item_value); },
- item_value_variant);
- }
- if constexpr (MeshType::Dimension == 3) {
- fout << "<PDataArray type=\"Int64\" IdType=\"1\" Name=\"faces\"/>\n";
- fout << "<PDataArray type=\"Int64\" IdType=\"1\" Name=\"faceoffsets\"/>\n";
- }
- fout << "</PCells>\n";
-
- fout << "<PPointData>\n";
- for (const auto& [name, item_value_variant] : output_named_item_value_set) {
- std::visit([&, name = name](auto&& item_value) { return this->_write_node_pvtu(fout, name, item_value); },
- item_value_variant);
- }
- fout << "</PPointData>\n";
-
- fout << "<PCellData>\n";
- for (const auto& [name, item_value_variant] : output_named_item_value_set) {
- std::visit([&, name = name](auto&& item_value) { return this->_write_cell_pvtu(fout, name, item_value); },
- item_value_variant);
- }
- fout << "</PCellData>\n";
-
- for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
- fout << "<Piece Source=\"" << _getFilenameVTU(i_rank) << "\"/>\n";
- }
- fout << "</PUnstructuredGrid>\n";
- fout << "</VTKFile>\n";
- }
-
- { // write VTK files
- std::ofstream fout(_getFilenameVTU(parallel::rank()));
- fout << "<?xml version=\"1.0\"?>\n";
- fout << "<VTKFile type=\"UnstructuredGrid\">\n";
- fout << "<UnstructuredGrid>\n";
- fout << "<Piece NumberOfPoints=\"" << mesh->numberOfNodes() << "\" NumberOfCells=\"" << mesh->numberOfCells()
- << "\">\n";
- fout << "<CellData>\n";
- for (const auto& [name, item_value_variant] : output_named_item_value_set) {
- std::visit([&, name = name](auto&& item_value) { return this->_write_cell_value(fout, name, item_value); },
- item_value_variant);
- }
- fout << "</CellData>\n";
- fout << "<PointData>\n";
- for (const auto& [name, item_value_variant] : output_named_item_value_set) {
- std::visit([&, name = name](auto&& item_value) { return this->_write_node_value(fout, name, item_value); },
- item_value_variant);
- }
- fout << "</PointData>\n";
- fout << "<Points>\n";
- {
- using Rd = TinyVector<MeshType::Dimension>;
- const NodeValue<const Rd>& xr = mesh->xr();
- Array<TinyVector<3>> positions(mesh->numberOfNodes());
- parallel_for(
- mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) {
- for (unsigned short i = 0; i < MeshType::Dimension; ++i) {
- positions[r][i] = xr[r][i];
- }
- for (unsigned short i = MeshType::Dimension; i < 3; ++i) {
- positions[r][i] = 0;
- }
- });
- _write_array(fout, "Positions", positions);
- }
- fout << "</Points>\n";
-
- fout << "<Cells>\n";
- {
- const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
-
- _write_array(fout, "connectivity", cell_to_node_matrix.entries());
- }
-
- {
- const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
- Array<unsigned int> offsets(mesh->numberOfCells());
- unsigned int offset = 0;
- for (CellId j = 0; j < mesh->numberOfCells(); ++j) {
- const auto& cell_nodes = cell_to_node_matrix[j];
- offset += cell_nodes.size();
- offsets[j] = offset;
- }
- _write_array(fout, "offsets", offsets);
- }
-
- {
- Array<int8_t> types(mesh->numberOfCells());
- const auto& cell_type = mesh->connectivity().cellType();
- const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
-
- parallel_for(
- mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
- switch (cell_type[j]) {
- case CellType::Line: {
- types[j] = 3;
- break;
- }
- case CellType::Triangle: {
- types[j] = 5;
- break;
- }
- case CellType::Quadrangle: {
- types[j] = 9;
- break;
- }
- case CellType::Tetrahedron: {
- types[j] = 10;
- break;
- }
- case CellType::Pyramid: {
- if (cell_to_node_matrix[j].size() == 5) {
- types[j] = 14; // quadrangle basis
- } else {
- types[j] = 41; // polygonal basis
- }
- break;
- }
- case CellType::Prism: {
- types[j] = 13;
- break;
- }
- case CellType::Hexahedron: {
- types[j] = 12;
- break;
- }
- case CellType::Diamond: {
- types[j] = 41;
- break;
- }
- default: {
- std::ostringstream os;
- os << __FILE__ << ':' << __LINE__ << ": unknown cell type";
- throw UnexpectedError(os.str());
- }
- }
- });
- _write_array(fout, "types", types);
- if constexpr (MeshType::Dimension == 3) {
- const bool has_general_polyhedron = [&] {
- for (size_t i = 0; i < types.size(); ++i) {
- if (types[i] == 41)
- return true;
- }
- return false;
- }();
- if (has_general_polyhedron) {
- const auto& cell_to_face_matrix = mesh->connectivity().cellToFaceMatrix();
- const auto& face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
- const auto& cell_face_is_reversed = mesh->connectivity().cellFaceIsReversed();
-
- Array<size_t> faces_offsets(mesh->numberOfCells());
- size_t next_offset = 0;
- fout << "<DataArray type=\"Int64\" IdType=\"1\" Name=\"faces\">\n";
- for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
- const auto& cell_nodes = cell_to_node_matrix[cell_id];
- std::unordered_map<size_t, size_t> node_id_to_local_number_in_cell;
- for (size_t i_cell_node = 0; i_cell_node < cell_nodes.size(); ++i_cell_node) {
- node_id_to_local_number_in_cell[cell_nodes[i_cell_node]] = i_cell_node;
- }
-
- const auto& cell_faces = cell_to_face_matrix[cell_id];
- fout << cell_faces.size() << '\n';
- next_offset++;
- for (size_t i_cell_face = 0; i_cell_face < cell_faces.size(); ++i_cell_face) {
- const FaceId& face_id = cell_faces[i_cell_face];
- const auto& face_nodes = face_to_node_matrix[face_id];
- fout << face_nodes.size();
- next_offset++;
- Array<size_t> face_node_in_cell(face_nodes.size());
-
- for (size_t i_face_node = 0; i_face_node < face_nodes.size(); ++i_face_node) {
- const NodeId& node_id = face_nodes[i_face_node];
- auto i_node_id_to_local_number_in_cell = node_id_to_local_number_in_cell.find(node_id);
- Assert(i_node_id_to_local_number_in_cell != node_id_to_local_number_in_cell.end());
- face_node_in_cell[i_face_node] = i_node_id_to_local_number_in_cell->second;
- }
-
- if (cell_face_is_reversed(cell_id, i_cell_face)) {
- for (size_t i = 0; i < face_node_in_cell.size(); ++i) {
- fout << ' ' << face_node_in_cell[face_node_in_cell.size() - 1 - i];
- }
- } else {
- for (size_t i = 0; i < face_node_in_cell.size(); ++i) {
- fout << ' ' << face_node_in_cell[i];
- }
- }
-
- next_offset += face_nodes.size();
-
- fout << '\n';
- }
- faces_offsets[cell_id] = next_offset;
- }
- fout << "</DataArray>\n";
- fout << "<DataArray type=\"Int64\" IdType=\"1\" Name=\"faceoffsets\">\n";
- for (size_t i_face_offsets = 0; i_face_offsets < faces_offsets.size(); ++i_face_offsets) {
- fout << faces_offsets[i_face_offsets] << '\n';
- }
- fout << "</DataArray>\n";
- }
- }
- }
-
- fout << "</Cells>\n";
- fout << "</Piece>\n";
- fout << "</UnstructuredGrid>\n";
- fout << "</VTKFile>\n";
- }
- m_file_number++;
- }
-
- VTKWriter(const std::string& base_filename, const double time_period)
- : m_base_filename(base_filename),
- m_file_number(0),
- m_last_time(-std::numeric_limits<double>::max()),
- m_time_period(time_period)
- {}
+ void _write(const std::shared_ptr<const MeshType>& mesh,
+ const OutputNamedItemValueSet& output_named_item_value_set,
+ double time) const;
+
+ public:
+ void writeMesh(const std::shared_ptr<const IMesh>& mesh) const final;
+
+ void write(const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
+ double time) const final;
+
+ VTKWriter(const std::string& base_filename, const double time_period) : WriterBase(base_filename, time_period) {}
~VTKWriter() = default;
};
diff --git a/src/output/WriterBase.cpp b/src/output/WriterBase.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0359cf7b3d33c8aa9fb33f8b62c671425e5958ba
--- /dev/null
+++ b/src/output/WriterBase.cpp
@@ -0,0 +1,201 @@
+#include <output/WriterBase.hpp>
+
+#include <mesh/IMesh.hpp>
+#include <output/OutputNamedItemValueSet.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/IDiscreteFunction.hpp>
+#include <scheme/IDiscreteFunctionDescriptor.hpp>
+#include <utils/Exceptions.hpp>
+
+template <size_t Dimension, typename DataType>
+void
+WriterBase::registerDiscreteFunctionP0(const std::string& name,
+ const IDiscreteFunction& i_discrete_function,
+ OutputNamedItemValueSet& named_item_value_set)
+{
+ const DiscreteFunctionP0<Dimension, DataType>& discrete_function =
+ dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(i_discrete_function);
+ named_item_value_set.add(NamedItemValue{name, discrete_function.cellValues()});
+}
+
+template <size_t Dimension>
+void
+WriterBase::registerDiscreteFunctionP0(const std::string& name,
+ const IDiscreteFunction& i_discrete_function,
+ OutputNamedItemValueSet& named_item_value_set)
+{
+ const ASTNodeDataType& data_type = i_discrete_function.dataType();
+ switch (data_type) {
+ case ASTNodeDataType::bool_t: {
+ registerDiscreteFunctionP0<Dimension, bool>(name, i_discrete_function, named_item_value_set);
+ break;
+ }
+ case ASTNodeDataType::unsigned_int_t: {
+ registerDiscreteFunctionP0<Dimension, uint64_t>(name, i_discrete_function, named_item_value_set);
+ break;
+ }
+ case ASTNodeDataType::int_t: {
+ registerDiscreteFunctionP0<Dimension, int64_t>(name, i_discrete_function, named_item_value_set);
+ break;
+ }
+ case ASTNodeDataType::double_t: {
+ registerDiscreteFunctionP0<Dimension, double>(name, i_discrete_function, named_item_value_set);
+ break;
+ }
+ case ASTNodeDataType::vector_t: {
+ switch (data_type.dimension()) {
+ case 1: {
+ registerDiscreteFunctionP0<Dimension, TinyVector<1, double>>(name, i_discrete_function, named_item_value_set);
+ break;
+ }
+ case 2: {
+ registerDiscreteFunctionP0<Dimension, TinyVector<2, double>>(name, i_discrete_function, named_item_value_set);
+ break;
+ }
+ case 3: {
+ registerDiscreteFunctionP0<Dimension, TinyVector<3, double>>(name, i_discrete_function, named_item_value_set);
+ break;
+ }
+ default: {
+ throw UnexpectedError("invalid vector dimension");
+ }
+ }
+ break;
+ }
+ case ASTNodeDataType::matrix_t: {
+ Assert(data_type.nbRows() == data_type.nbColumns(), "invalid matrix dimensions");
+ switch (data_type.nbRows()) {
+ case 1: {
+ registerDiscreteFunctionP0<Dimension, TinyMatrix<1, double>>(name, i_discrete_function, named_item_value_set);
+ break;
+ }
+ case 2: {
+ registerDiscreteFunctionP0<Dimension, TinyMatrix<2, double>>(name, i_discrete_function, named_item_value_set);
+ break;
+ }
+ case 3: {
+ registerDiscreteFunctionP0<Dimension, TinyMatrix<3, double>>(name, i_discrete_function, named_item_value_set);
+ break;
+ }
+ default: {
+ throw UnexpectedError("invalid matrix dimension");
+ }
+ }
+ break;
+ }
+ default: {
+ throw UnexpectedError("invalid data type " + dataTypeName(data_type));
+ }
+ }
+}
+
+void
+WriterBase::registerDiscreteFunctionP0(const NamedDiscreteFunction& named_discrete_function,
+ OutputNamedItemValueSet& named_item_value_set)
+{
+ const IDiscreteFunction& i_discrete_function = *named_discrete_function.discreteFunction();
+ const std::string& name = named_discrete_function.name();
+ switch (i_discrete_function.mesh()->dimension()) {
+ case 1: {
+ registerDiscreteFunctionP0<1>(name, i_discrete_function, named_item_value_set);
+ break;
+ }
+ case 2: {
+ registerDiscreteFunctionP0<2>(name, i_discrete_function, named_item_value_set);
+ break;
+ }
+ case 3: {
+ registerDiscreteFunctionP0<3>(name, i_discrete_function, named_item_value_set);
+ break;
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ }
+}
+
+std::shared_ptr<const IMesh>
+WriterBase::_getMesh(
+ const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list) const
+{
+ Assert(named_discrete_function_list.size() > 0);
+
+ std::shared_ptr mesh = named_discrete_function_list[0]->discreteFunction()->mesh();
+ for (size_t i = 1; i < named_discrete_function_list.size(); ++i) {
+ if (mesh != named_discrete_function_list[i]->discreteFunction()->mesh()) {
+ std::ostringstream error_msg;
+ error_msg << "discrete functions must be defined on the same mesh!\n"
+ << rang::fgB::yellow << "note:" << rang::style::reset << " cannot write " << rang::fgB::blue
+ << named_discrete_function_list[0]->name() << rang::style::reset << " and " << rang::fgB::blue
+ << named_discrete_function_list[i]->name() << rang::style::reset << " in the same file.";
+ throw NormalError(error_msg.str());
+ }
+ }
+
+ return mesh;
+}
+
+OutputNamedItemValueSet
+WriterBase::_getOutputNamedItemValueSet(
+ const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list) const
+{
+ OutputNamedItemValueSet named_item_value_set;
+
+ for (auto& named_discrete_function : named_discrete_function_list) {
+ const IDiscreteFunction& i_discrete_function = *named_discrete_function->discreteFunction();
+
+ switch (i_discrete_function.descriptor().type()) {
+ case DiscreteFunctionType::P0: {
+ WriterBase::registerDiscreteFunctionP0(*named_discrete_function, named_item_value_set);
+ break;
+ }
+ default: {
+ std::ostringstream error_msg;
+ error_msg << "the type of discrete function of " << rang::fgB::blue << named_discrete_function->name()
+ << rang::style::reset << " is not supported";
+ throw NormalError(error_msg.str());
+ }
+ }
+ }
+
+ return named_item_value_set;
+}
+
+double
+WriterBase::getLastTime() const
+{
+ if (m_saved_times.size() > 0) {
+ return m_saved_times[m_saved_times.size() - 1];
+ } else {
+ return -std::numeric_limits<double>::max();
+ }
+}
+
+WriterBase::WriterBase(const std::string& base_filename, const double& time_period)
+ : m_base_filename{base_filename}, m_time_period{time_period}
+{}
+
+void
+WriterBase::writeIfNeeded(const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
+ double time) const
+{
+ const double last_time = getLastTime();
+ if (time == last_time)
+ return; // output already performed
+
+ if (time >= last_time + m_time_period) {
+ this->write(named_discrete_function_list, time);
+ m_saved_times.push_back(time);
+ }
+}
+
+void
+WriterBase::writeForced(const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
+ double time) const
+{
+ if (time == getLastTime())
+ return; // output already performed
+
+ this->write(named_discrete_function_list, time);
+ m_saved_times.push_back(time);
+}
diff --git a/src/output/WriterBase.hpp b/src/output/WriterBase.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..6e9a08ad123361178327ddb984e5a3da5bdfa909
--- /dev/null
+++ b/src/output/WriterBase.hpp
@@ -0,0 +1,54 @@
+#ifndef WRITER_BASE_HPP
+#define WRITER_BASE_HPP
+
+#include <output/IWriter.hpp>
+
+#include <string>
+
+class IMesh;
+class OutputNamedItemValueSet;
+
+class WriterBase : public IWriter
+{
+ protected:
+ const std::string m_base_filename;
+ const double m_time_period;
+
+ mutable std::vector<double> m_saved_times;
+
+ private:
+ template <size_t Dimension, typename DataType>
+ static void registerDiscreteFunctionP0(const std::string& name, const IDiscreteFunction&, OutputNamedItemValueSet&);
+
+ template <size_t Dimension>
+ static void registerDiscreteFunctionP0(const std::string& name, const IDiscreteFunction&, OutputNamedItemValueSet&);
+
+ static void registerDiscreteFunctionP0(const NamedDiscreteFunction&, OutputNamedItemValueSet&);
+
+ protected:
+ std::shared_ptr<const IMesh> _getMesh(
+ const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list) const;
+
+ OutputNamedItemValueSet _getOutputNamedItemValueSet(
+ const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list) const;
+
+ virtual void write(const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
+ double time) const = 0;
+
+ public:
+ double getLastTime() const;
+
+ void writeIfNeeded(const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
+ double time) const final;
+
+ void writeForced(const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
+ double time) const final;
+
+ WriterBase() = delete;
+
+ WriterBase(const std::string& base_filename, const double& time_period);
+
+ virtual ~WriterBase() = default;
+};
+
+#endif // WRITER_BASE_HPP
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..1467393f2cbcbd9e31dd495e2905e9be478a1d78
--- /dev/null
+++ b/src/scheme/CMakeLists.txt
@@ -0,0 +1,5 @@
+# ------------------- Source files --------------------
+
+add_library(
+ PugsScheme
+ DiscreteFunctionInterpoler.cpp)
diff --git a/src/scheme/DiscreteFunctionDescriptorP0.hpp b/src/scheme/DiscreteFunctionDescriptorP0.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..97bfbbf729d4224dbb7af41bf61d489c5a3080cb
--- /dev/null
+++ b/src/scheme/DiscreteFunctionDescriptorP0.hpp
@@ -0,0 +1,24 @@
+#ifndef DISCRETE_FUNCTION_DESCRIPTOR_P_0HPP
+#define DISCRETE_FUNCTION_DESCRIPTOR_P_0HPP
+
+#include <scheme/IDiscreteFunctionDescriptor.hpp>
+
+class DiscreteFunctionDescriptorP0 : public IDiscreteFunctionDescriptor
+{
+ public:
+ DiscreteFunctionType
+ type() const final
+ {
+ return DiscreteFunctionType::P0;
+ }
+
+ DiscreteFunctionDescriptorP0() noexcept = default;
+
+ DiscreteFunctionDescriptorP0(const DiscreteFunctionDescriptorP0&) = default;
+
+ DiscreteFunctionDescriptorP0(DiscreteFunctionDescriptorP0&&) noexcept = default;
+
+ ~DiscreteFunctionDescriptorP0() noexcept = default;
+};
+
+#endif // DISCRETE_FUNCTION_DESCRIPTOR_P_0HPP
diff --git a/src/scheme/DiscreteFunctionInterpoler.cpp b/src/scheme/DiscreteFunctionInterpoler.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..fda64fdcb82d1eb36b8925e9cdc438989cd92306
--- /dev/null
+++ b/src/scheme/DiscreteFunctionInterpoler.cpp
@@ -0,0 +1,103 @@
+#include <scheme/DiscreteFunctionInterpoler.hpp>
+
+#include <scheme/DiscreteFunctionP0.hpp>
+#include <utils/Exceptions.hpp>
+
+template <size_t Dimension, typename DataType>
+std::shared_ptr<IDiscreteFunction>
+DiscreteFunctionInterpoler::_interpolate() const
+{
+ std::shared_ptr mesh = std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(m_mesh);
+ return std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(mesh, m_function_id);
+}
+
+template <size_t Dimension>
+std::shared_ptr<IDiscreteFunction>
+DiscreteFunctionInterpoler::_interpolate() const
+{
+ const auto& function_descriptor = m_function_id.descriptor();
+ Assert(function_descriptor.domainMappingNode().children[1]->m_data_type == ASTNodeDataType::typename_t);
+
+ const ASTNodeDataType& data_type = function_descriptor.domainMappingNode().children[1]->m_data_type.contentType();
+
+ switch (data_type) {
+ case ASTNodeDataType::bool_t: {
+ return this->_interpolate<Dimension, bool>();
+ }
+ case ASTNodeDataType::unsigned_int_t: {
+ return this->_interpolate<Dimension, uint64_t>();
+ }
+ case ASTNodeDataType::int_t: {
+ return this->_interpolate<Dimension, int64_t>();
+ }
+ case ASTNodeDataType::double_t: {
+ return this->_interpolate<Dimension, double>();
+ }
+ case ASTNodeDataType::vector_t: {
+ switch (data_type.dimension()) {
+ case 1: {
+ return this->_interpolate<Dimension, TinyVector<1>>();
+ }
+ case 2: {
+ return this->_interpolate<Dimension, TinyVector<2>>();
+ }
+ case 3: {
+ return this->_interpolate<Dimension, TinyVector<3>>();
+ }
+ default: {
+ std::ostringstream os;
+ os << "invalid vector dimension " << rang::fgB::red << data_type.dimension() << rang::style::reset;
+
+ throw UnexpectedError(os.str());
+ }
+ }
+ }
+ case ASTNodeDataType::matrix_t: {
+ Assert(data_type.nbColumns() == data_type.nbRows(), "undefined matrix type");
+ switch (data_type.nbColumns()) {
+ case 1: {
+ return this->_interpolate<Dimension, TinyMatrix<1>>();
+ }
+ case 2: {
+ return this->_interpolate<Dimension, TinyMatrix<2>>();
+ }
+ case 3: {
+ return this->_interpolate<Dimension, TinyMatrix<3>>();
+ }
+ default: {
+ std::ostringstream os;
+ os << "invalid vector dimension " << rang::fgB::red << data_type.dimension() << rang::style::reset;
+
+ throw UnexpectedError(os.str());
+ }
+ }
+ }
+ default: {
+ std::ostringstream os;
+ os << "invalid interpolation value type: " << rang::fgB::red << dataTypeName(data_type) << rang::style::reset;
+
+ throw UnexpectedError(os.str());
+ }
+ }
+}
+
+std::shared_ptr<IDiscreteFunction>
+DiscreteFunctionInterpoler::interpolate() const
+{
+ std::shared_ptr<IDiscreteFunction> discrete_function;
+ switch (m_mesh->dimension()) {
+ case 1: {
+ return this->_interpolate<1>();
+ }
+ case 2: {
+ return this->_interpolate<2>();
+ }
+ case 3: {
+ return this->_interpolate<3>();
+ }
+ default: {
+ throw UnexpectedError("invalid dimension");
+ }
+ }
+ return nullptr;
+}
diff --git a/src/scheme/DiscreteFunctionInterpoler.hpp b/src/scheme/DiscreteFunctionInterpoler.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..758ff2454ec673ba2cc46c0c561eca1a4a256b51
--- /dev/null
+++ b/src/scheme/DiscreteFunctionInterpoler.hpp
@@ -0,0 +1,39 @@
+#ifndef DISCRETE_FUNCTION_INTERPOLER_HPP
+#define DISCRETE_FUNCTION_INTERPOLER_HPP
+
+#include <language/utils/FunctionSymbolId.hpp>
+#include <mesh/IMesh.hpp>
+#include <scheme/IDiscreteFunction.hpp>
+#include <scheme/IDiscreteFunctionDescriptor.hpp>
+
+#include <memory>
+
+class DiscreteFunctionInterpoler
+{
+ private:
+ std::shared_ptr<const IMesh> m_mesh;
+ std::shared_ptr<const IDiscreteFunctionDescriptor> m_discrete_function_descriptor;
+ const FunctionSymbolId m_function_id;
+
+ template <size_t Dimension, typename DataType>
+ std::shared_ptr<IDiscreteFunction> _interpolate() const;
+
+ template <size_t Dimension>
+ std::shared_ptr<IDiscreteFunction> _interpolate() const;
+
+ public:
+ std::shared_ptr<IDiscreteFunction> interpolate() const;
+
+ DiscreteFunctionInterpoler(const std::shared_ptr<const IMesh>& mesh,
+ const std::shared_ptr<const IDiscreteFunctionDescriptor>& discrete_function_descriptor,
+ const FunctionSymbolId& function_id)
+ : m_mesh{mesh}, m_discrete_function_descriptor{discrete_function_descriptor}, m_function_id{function_id}
+ {}
+
+ DiscreteFunctionInterpoler(const DiscreteFunctionInterpoler&) = delete;
+ DiscreteFunctionInterpoler(DiscreteFunctionInterpoler&&) = delete;
+
+ ~DiscreteFunctionInterpoler() = default;
+};
+
+#endif // DISCRETE_FUNCTION_INTERPOLER_HPP
diff --git a/src/scheme/DiscreteFunctionP0.hpp b/src/scheme/DiscreteFunctionP0.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..bf97079b776c80fec10e26608cfc6ff5e40a664a
--- /dev/null
+++ b/src/scheme/DiscreteFunctionP0.hpp
@@ -0,0 +1,72 @@
+#ifndef DISCRETE_FUNCTION_P0_HPP
+#define DISCRETE_FUNCTION_P0_HPP
+
+#include <scheme/IDiscreteFunction.hpp>
+
+#include <language/utils/InterpolateItemValue.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <scheme/DiscreteFunctionDescriptorP0.hpp>
+
+template <size_t Dimension, typename DataType>
+class DiscreteFunctionP0 : public IDiscreteFunction
+{
+ private:
+ using MeshType = Mesh<Connectivity<Dimension>>;
+
+ std::shared_ptr<const MeshType> m_mesh;
+ CellValue<DataType> m_cell_values;
+
+ DiscreteFunctionDescriptorP0 m_discrete_function_descriptor;
+
+ public:
+ ASTNodeDataType
+ dataType() const final
+ {
+ return ast_node_data_type_from<DataType>;
+ }
+
+ CellValue<DataType>
+ cellValues() const
+ {
+ return m_cell_values;
+ }
+
+ std::shared_ptr<const IMesh>
+ mesh() const
+ {
+ return m_mesh;
+ }
+
+ const IDiscreteFunctionDescriptor&
+ descriptor() const final
+ {
+ return m_discrete_function_descriptor;
+ }
+
+ PUGS_FORCEINLINE
+ DataType&
+ operator[](const CellId& cell_id) const noexcept(NO_ASSERT)
+ {
+ return m_cell_values[cell_id];
+ }
+
+ DiscreteFunctionP0(const std::shared_ptr<const MeshType>& mesh, const FunctionSymbolId& function_id) : m_mesh(mesh)
+ {
+ using MeshDataType = MeshData<Dimension>;
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+ m_cell_values =
+ InterpolateItemValue<DataType(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(function_id,
+ mesh_data.xj());
+ }
+
+ DiscreteFunctionP0(const DiscreteFunctionP0&) noexcept = default;
+ DiscreteFunctionP0(DiscreteFunctionP0&&) noexcept = default;
+
+ ~DiscreteFunctionP0() = default;
+};
+
+#endif // DISCRETE_FUNCTION_P0_HPP
diff --git a/src/scheme/DiscreteFunctionType.hpp b/src/scheme/DiscreteFunctionType.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..debf863499df64cd161ac0b7c20a1a5ce4657107
--- /dev/null
+++ b/src/scheme/DiscreteFunctionType.hpp
@@ -0,0 +1,9 @@
+#ifndef DISCRETE_FUNCTION_TYPE_HPP
+#define DISCRETE_FUNCTION_TYPE_HPP
+
+enum class DiscreteFunctionType
+{
+ P0
+};
+
+#endif // DISCRETE_FUNCTION_TYPE_HPP
diff --git a/src/scheme/IDiscreteFunction.hpp b/src/scheme/IDiscreteFunction.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a8bf77f43efb48b03a8967eb16d81966fde86570
--- /dev/null
+++ b/src/scheme/IDiscreteFunction.hpp
@@ -0,0 +1,26 @@
+#ifndef I_DISCRETE_FUNCTION_HPP
+#define I_DISCRETE_FUNCTION_HPP
+
+class IMesh;
+class IDiscreteFunctionDescriptor;
+
+#include <language/utils/ASTNodeDataTypeTraits.hpp>
+#include <memory>
+
+class IDiscreteFunction
+{
+ public:
+ virtual std::shared_ptr<const IMesh> mesh() const = 0;
+ virtual const IDiscreteFunctionDescriptor& descriptor() const = 0;
+ virtual ASTNodeDataType dataType() const = 0;
+
+ IDiscreteFunction() = default;
+
+ IDiscreteFunction(const IDiscreteFunction&) = default;
+
+ IDiscreteFunction(IDiscreteFunction&&) noexcept = default;
+
+ virtual ~IDiscreteFunction() noexcept = default;
+};
+
+#endif // I_DISCRETE_FUNCTION_HPP
diff --git a/src/scheme/IDiscreteFunctionDescriptor.hpp b/src/scheme/IDiscreteFunctionDescriptor.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..19d411507e087c0e3262f8e712aa393d96796b20
--- /dev/null
+++ b/src/scheme/IDiscreteFunctionDescriptor.hpp
@@ -0,0 +1,20 @@
+#ifndef I_DISCRETE_FUNCTION_DESCRIPTOR_HPP
+#define I_DISCRETE_FUNCTION_DESCRIPTOR_HPP
+
+#include <scheme/DiscreteFunctionType.hpp>
+
+class IDiscreteFunctionDescriptor
+{
+ public:
+ virtual DiscreteFunctionType type() const = 0;
+
+ IDiscreteFunctionDescriptor() noexcept = default;
+
+ IDiscreteFunctionDescriptor(const IDiscreteFunctionDescriptor&) = default;
+
+ IDiscreteFunctionDescriptor(IDiscreteFunctionDescriptor&&) noexcept = default;
+
+ virtual ~IDiscreteFunctionDescriptor() noexcept = default;
+};
+
+#endif // I_DISCRETE_FUNCTION_DESCRIPTOR_HPP
diff --git a/src/utils/PugsTraits.hpp b/src/utils/PugsTraits.hpp
index f51b455cb5ea4efa678dd980d1fec2721ab1fd2a..e1dc8b48d6b15fa05ac9b9b9378385e351fadb2b 100644
--- a/src/utils/PugsTraits.hpp
+++ b/src/utils/PugsTraits.hpp
@@ -3,6 +3,7 @@
#include <cstddef>
#include <memory>
+#include <tuple>
#include <type_traits>
#include <variant>
#include <vector>
@@ -40,7 +41,7 @@ inline constexpr bool is_shared_ptr_v = false;
template <typename T>
inline constexpr bool is_shared_ptr_v<std::shared_ptr<T>> = true;
-// Traits is_shared_ptr
+// Traits is_unique_ptr
template <typename T>
inline constexpr bool is_unique_ptr_v = false;
@@ -76,6 +77,14 @@ inline constexpr bool is_std_vector_v = false;
template <typename T>
inline constexpr bool is_std_vector_v<std::vector<T>> = true;
+// Traits is_std_tuple
+
+template <typename... T>
+inline constexpr bool is_std_tuple_v = false;
+
+template <typename... T>
+inline constexpr bool is_std_tuple_v<std::tuple<T...>> = true;
+
// Traits is_tiny_vector
template <typename T>
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 2ae89a2790455514249a6097ff4dbe748e7a0839..c24b57517dc2261a930605823d0c1a736d6c6a2d 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -116,6 +116,8 @@ target_link_libraries (unit_tests
PugsLanguage
PugsMesh
PugsAlgebra
+ PugsScheme
+ PugsOutput
PugsUtils
kokkos
${PARMETIS_LIBRARIES}
@@ -137,7 +139,10 @@ target_link_libraries (mpi_unit_tests
PugsMesh
PugsAlgebra
PugsUtils
- PugsLanguageUtils PugsUtils
+ PugsLanguageUtils
+ PugsScheme
+ PugsOutput
+ PugsUtils
PugsAlgebra
PugsMesh
kokkos
diff --git a/tests/test_BuiltinFunctionEmbedder.cpp b/tests/test_BuiltinFunctionEmbedder.cpp
index 2651c0e94659c1bbffc0bd47e6ab946cce18e3c9..ba88ac0b0a0e3d6e7c3a16977e03437a1e14a05c 100644
--- a/tests/test_BuiltinFunctionEmbedder.cpp
+++ b/tests/test_BuiltinFunctionEmbedder.cpp
@@ -286,6 +286,62 @@ TEST_CASE("BuiltinFunctionEmbedder", "[language]")
REQUIRE(i_embedded_c->getReturnDataType() == ASTNodeDataType::double_t);
}
+ SECTION("R*R -> R*R^2*shared_double BuiltinFunctionEmbedder")
+ {
+ std::function c = [](double a, double b) -> std::tuple<double, TinyVector<2>, std::shared_ptr<double>> {
+ return std::make_tuple(a + b, TinyVector<2>{b, -a}, std::make_shared<double>(a - b));
+ };
+
+ std::unique_ptr<IBuiltinFunctionEmbedder> i_embedded_c = std::make_unique<
+ BuiltinFunctionEmbedder<std::tuple<double, TinyVector<2>, std::shared_ptr<double>>(double, double)>>(c);
+
+ const double a = 3.2;
+ const double b = 1.5;
+
+ REQUIRE(a + b == std::get<0>(c(a, b)));
+ REQUIRE(TinyVector<2>{b, -a} == std::get<1>(c(a, b)));
+ REQUIRE(a - b == *std::get<2>(c(a, b)));
+ const AggregateDataVariant value_list =
+ std::get<AggregateDataVariant>(i_embedded_c->apply(std::vector<DataVariant>{a, b}));
+
+ REQUIRE(std::get<double>(value_list[0]) == a + b);
+ REQUIRE(std::get<TinyVector<2>>(value_list[1]) == TinyVector<2>{b, -a});
+ auto data_type = i_embedded_c->getReturnDataType();
+ REQUIRE(data_type == ASTNodeDataType::list_t);
+
+ REQUIRE(*data_type.contentTypeList()[0] == ASTNodeDataType::double_t);
+ REQUIRE(*data_type.contentTypeList()[1] == ASTNodeDataType::build<ASTNodeDataType::vector_t>(2));
+ REQUIRE(*data_type.contentTypeList()[2] == ast_node_data_type_from<std::shared_ptr<double>>);
+ }
+
+ SECTION("void -> N*R*shared_double BuiltinFunctionEmbedder")
+ {
+ std::function c = [](void) -> std::tuple<uint64_t, double, std::shared_ptr<double>> {
+ uint64_t a = 1;
+ double b = 3.5;
+ return std::make_tuple(a, b, std::make_shared<double>(a + b));
+ };
+
+ std::unique_ptr<IBuiltinFunctionEmbedder> i_embedded_c =
+ std::make_unique<BuiltinFunctionEmbedder<std::tuple<uint64_t, double, std::shared_ptr<double>>(void)>>(c);
+
+ REQUIRE(1ul == std::get<0>(c()));
+ REQUIRE(3.5 == std::get<1>(c()));
+ REQUIRE((1ul + 3.5) == *std::get<2>(c()));
+ const AggregateDataVariant value_list =
+ std::get<AggregateDataVariant>(i_embedded_c->apply(std::vector<DataVariant>{}));
+
+ REQUIRE(std::get<uint64_t>(value_list[0]) == 1ul);
+ REQUIRE(std::get<double>(value_list[1]) == 3.5);
+
+ auto data_type = i_embedded_c->getReturnDataType();
+ REQUIRE(data_type == ASTNodeDataType::list_t);
+
+ REQUIRE(*data_type.contentTypeList()[0] == ASTNodeDataType::unsigned_int_t);
+ REQUIRE(*data_type.contentTypeList()[1] == ASTNodeDataType::double_t);
+ REQUIRE(*data_type.contentTypeList()[2] == ast_node_data_type_from<std::shared_ptr<double>>);
+ }
+
SECTION("void(void) BuiltinFunctionEmbedder")
{
std::function c = [](void) -> void {};