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/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 542690a4ffcfb73ba174f0014b358548208d830f..c23d3794e4c26a95301599b6933fe0d643391010 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -5,10 +5,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>
@@ -18,9 +22,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 96418b1693d377641f9bd096d263d5d17c6b929d..3cfbb254a8f6e2ac36d16c1a78e53f3f995f56fc 100644
--- a/src/language/modules/SchemeModule.hpp
+++ b/src/language/modules/SchemeModule.hpp
@@ -15,6 +15,16 @@ template <>
inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const IBoundaryConditionDescriptor>> =
ASTNodeDataType::build<ASTNodeDataType::type_id_t>("boundary_condition");
+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/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/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/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/tests/CMakeLists.txt b/tests/CMakeLists.txt
index b2c9f5015a442538ec23f8d421fb07fa93153480..c4413ddc7c6f06665b4efb7b091548fb26140fab 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -113,6 +113,8 @@ target_link_libraries (unit_tests
PugsLanguage
PugsMesh
PugsAlgebra
+ PugsScheme
+ PugsOutput
PugsUtils
kokkos
${PARMETIS_LIBRARIES}
@@ -134,7 +136,10 @@ target_link_libraries (mpi_unit_tests
PugsMesh
PugsAlgebra
PugsUtils
- PugsLanguageUtils PugsUtils
+ PugsLanguageUtils
+ PugsScheme
+ PugsOutput
+ PugsUtils
PugsAlgebra
PugsMesh
kokkos