diff --git a/CMakeLists.txt b/CMakeLists.txt
index 97a94a4eff3e7b956a219ca9399b5221157d6500..114ceecc5d49685a513708dced3dcffea83bff29 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -530,6 +530,7 @@ target_link_libraries(
PugsAlgebra
PugsScheme
PugsUtils
+ PugsOutput
PugsLanguageUtils
kokkos
${PETSC_LIBRARIES}
@@ -552,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 e9c207b2fc2f256c6e17a35afbafcb82521dea75..16904a7be0cddc9aa54206638ab4b5c759cd0de8 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -19,4 +19,4 @@ add_subdirectory(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..a2a456f0ebc91476b2dbb541e9fdc6eecc77c83c 100644
--- a/src/language/algorithms/AcousticSolverAlgorithm.cpp
+++ b/src/language/algorithms/AcousticSolverAlgorithm.cpp
@@ -7,6 +7,8 @@
#include <scheme/IBoundaryDescriptor.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+#include <iomanip>
+
template <size_t Dimension>
AcousticSolverAlgorithm<Dimension>::AcousticSolverAlgorithm(
const AcousticSolverType& solver_type,
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index e691e248adb885410d7a50d8994734517d4fc921..c23d3794e4c26a95301599b6933fe0d643391010 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -5,6 +5,7 @@
#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>
diff --git a/src/language/modules/VTKModule.cpp b/src/language/modules/VTKModule.cpp
index 3144460fb5729ba29db2f20f864d6a5552e74a9b..9c6b1838f641cb2bab7ac36886b01da8ee398643 100644
--- a/src/language/modules/VTKModule.cpp
+++ b/src/language/modules/VTKModule.cpp
@@ -5,118 +5,164 @@
#include <mesh/Connectivity.hpp>
#include <mesh/GmshReader.hpp>
#include <mesh/Mesh.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>
VTKModule::VTKModule()
{
- this->_addBuiltinFunction("writeMeshVTK",
- std::make_shared<
- BuiltinFunctionEmbedder<void(std::shared_ptr<const IMesh>, const std::string&)>>(
+ this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const NamedDiscreteFunction>>);
- [](std::shared_ptr<const IMesh> p_mesh, const std::string& filename) -> void {
- VTKWriter writer(filename, 0.1);
+ this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const IWriter>>);
- static double time = 0;
+ this->_addBuiltinFunction("vtk_writer",
+ std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IWriter>(const std::string&,
+ const double&)>>(
- 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++;
+ [](const std::string& filename, const double& period) {
+ return std::make_shared<VTKWriter>(filename, period);
}
));
- this->_addBuiltinFunction("writeVTK",
+ this->_addBuiltinFunction("name_output",
std::make_shared<BuiltinFunctionEmbedder<
- void(const std::vector<std::shared_ptr<const IDiscreteFunction>>&, const std::string&)>>(
+ std::shared_ptr<const NamedDiscreteFunction>(std::shared_ptr<const IDiscreteFunction>,
+ const std::string&)>>(
- [](const std::vector<std::shared_ptr<const IDiscreteFunction>>& discrete_function_list,
- const std::string& filename) -> void {
- Assert(discrete_function_list.size() > 0);
+ [](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);
+ }
- VTKWriter writer(filename, 0.1);
+ ));
- static double time = 0;
+ this->_addBuiltinFunction("write_mesh",
+ std::make_shared<BuiltinFunctionEmbedder<void(std::shared_ptr<const IWriter>,
+ std::shared_ptr<const IMesh>)>>(
- std::shared_ptr p_mesh = discrete_function_list[0]->mesh();
- for (size_t i = 1; i < discrete_function_list.size(); ++i) {
- if (p_mesh != discrete_function_list[i]->mesh()) {
- throw NormalError("discrete functions must be defined on the same mesh!");
- }
- }
+ [](std::shared_ptr<const IWriter> writer, std::shared_ptr<const IMesh> p_mesh) -> void {
+ static double number = 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);
+ const std::shared_ptr mesh =
+ std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(p_mesh);
+
+ writer->write(mesh,
+ OutputNamedItemValueSet{
+ NamedItemValue{"cell_number", mesh->connectivity().cellNumber()}},
+ number, 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);
+ const std::shared_ptr mesh =
+ std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(p_mesh);
+
+ writer->write(mesh,
+ OutputNamedItemValueSet{
+ NamedItemValue{"cell_number", mesh->connectivity().cellNumber()}},
+ number, 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);
+ const std::shared_ptr mesh =
+ std::dynamic_pointer_cast<const Mesh<Connectivity<3>>>(p_mesh);
+
+ writer->write(mesh,
+ OutputNamedItemValueSet{
+ NamedItemValue{"cell_number", mesh->connectivity().cellNumber()}},
+ number, true);
break;
}
}
- time++;
+ number++;
}
));
+
+ this->_addBuiltinFunction(
+ "writeVTK",
+ 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 {
+ Assert(named_discrete_function_list.size() > 0);
+
+ std::shared_ptr p_mesh = named_discrete_function_list[0]->discreteFunction()->mesh();
+ for (size_t i = 1; i < named_discrete_function_list.size(); ++i) {
+ if (p_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());
+ }
+ }
+
+ 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: {
+ const DiscreteFunctionP0<2, double>& discrete_function =
+ dynamic_cast<const DiscreteFunctionP0<2, double>&>(i_discrete_function);
+
+ named_item_value_set.add(NamedItemValue{named_discrete_function->name(), discrete_function.cellValues()});
+ 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());
+ }
+ }
+ }
+
+ 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);
+
+ named_item_value_set.add(NamedItemValue{"cell_number", mesh->connectivity().cellNumber()});
+ named_item_value_set.add(NamedItemValue{"node_number", mesh->connectivity().nodeNumber()});
+
+ writer->write(mesh, named_item_value_set, 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);
+
+ named_item_value_set.add(NamedItemValue{"cell_number", mesh->connectivity().cellNumber()});
+ named_item_value_set.add(NamedItemValue{"node_number", mesh->connectivity().nodeNumber()});
+
+ writer->write(mesh, named_item_value_set, 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);
+
+ named_item_value_set.add(NamedItemValue{"cell_number", mesh->connectivity().cellNumber()});
+ named_item_value_set.add(NamedItemValue{"node_number", mesh->connectivity().nodeNumber()});
+
+ writer->write(mesh, named_item_value_set, time, true);
+ break;
+ }
+ }
+ }
+
+ ));
}
diff --git a/src/language/modules/VTKModule.hpp b/src/language/modules/VTKModule.hpp
index 1578fad811b7dae1181053317e891dbb865ad035..55d413e888c0edef1ba138d09ad265df2dd7f2bf 100644
--- a/src/language/modules/VTKModule.hpp
+++ b/src/language/modules/VTKModule.hpp
@@ -2,8 +2,19 @@
#define VTK_MODULE_HPP
#include <language/modules/BuiltinModule.hpp>
+#include <language/utils/ASTNodeDataTypeTraits.hpp>
#include <utils/PugsMacros.hpp>
+class NamedDiscreteFunction;
+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 VTKModule : public BuiltinModule
{
public:
diff --git a/src/output/CMakeLists.txt b/src/output/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..e5b860bc0bf2781a16d7453537e6657600d55250
--- /dev/null
+++ b/src/output/CMakeLists.txt
@@ -0,0 +1,14 @@
+# ------------------- Source files --------------------
+
+add_library(
+ PugsOutput
+ VTKWriter.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/IWriter.hpp b/src/output/IWriter.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..98df0858ee121db43119ba085e783ae757c3022d
--- /dev/null
+++ b/src/output/IWriter.hpp
@@ -0,0 +1,19 @@
+#ifndef I_WRITER_HPP
+#define I_WRITER_HPP
+
+#include <output/OutputNamedItemValueSet.hpp>
+class IMesh;
+
+class IWriter
+{
+ public:
+ virtual void write(const std::shared_ptr<const IMesh>& mesh,
+ const OutputNamedItemValueSet& output_named_item_value_set,
+ double time,
+ bool forced_output = false) 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..c474419eeca269644fa99f4b6c45fa596439c5e0 100644
--- a/src/output/OutputNamedItemValueSet.hpp
+++ b/src/output/OutputNamedItemValueSet.hpp
@@ -103,6 +103,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..11a5d4bb5774076429bae9cbedce05a7b90b7e91
--- /dev/null
+++ b/src/output/VTKWriter.cpp
@@ -0,0 +1,468 @@
+#include <output/VTKWriter.hpp>
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
+#include <utils/Messenger.hpp>
+
+#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_file_number << ".pvtu";
+ return sout.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_file_number << ".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 <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 <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 <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 <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& output_named_item_value_set,
+ double time,
+ bool forced_output) const
+{
+ 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++;
+}
+
+void
+VTKWriter::write(const std::shared_ptr<const IMesh>& mesh,
+ const OutputNamedItemValueSet& output_named_item_value_set,
+ double time,
+ bool forced_output) const
+{
+ switch (mesh->dimension()) {
+ case 1: {
+ this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(mesh), output_named_item_value_set, time,
+ forced_output);
+ break;
+ }
+ case 2: {
+ this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(mesh), output_named_item_value_set, time,
+ forced_output);
+ break;
+ }
+ case 3: {
+ this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<3>>>(mesh), output_named_item_value_set, time,
+ forced_output);
+ break;
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ }
+}
diff --git a/src/output/VTKWriter.hpp b/src/output/VTKWriter.hpp
index f3c202ffc51633d291d36150d662fe950023ef58..7661d0a27257a5ae1dbad96726d3e0263f0eaa21 100644
--- a/src/output/VTKWriter.hpp
+++ b/src/output/VTKWriter.hpp
@@ -1,449 +1,91 @@
#ifndef VTK_WRITER_HPP
#define VTK_WRITER_HPP
+#include <output/IWriter.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 IWriter
{
private:
const std::string m_base_filename;
- unsigned int m_file_number;
- double m_last_time;
+ mutable unsigned int m_file_number;
+ mutable 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 _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 <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>>&)
- {
- 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 TinyVector<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 <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 <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";
+ void _write(const std::shared_ptr<const MeshType>& mesh,
+ const OutputNamedItemValueSet& output_named_item_value_set,
+ double time,
+ bool forced_output = false) const;
- 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++;
- }
+ public:
+ void write(const std::shared_ptr<const IMesh>& mesh,
+ const OutputNamedItemValueSet& output_named_item_value_set,
+ double time,
+ bool forced_output = false) const final;
VTKWriter(const std::string& base_filename, const double time_period)
: m_base_filename(base_filename),
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
index 60b481cdeffcbae4571c5a297a5521099d01686d..7ac4919ab77bd0be973925a88b59a94cf7dc57d5 100644
--- a/src/scheme/DiscreteFunctionInterpoler.cpp
+++ b/src/scheme/DiscreteFunctionInterpoler.cpp
@@ -1,50 +1,8 @@
#include <scheme/DiscreteFunctionInterpoler.hpp>
-#include <language/utils/InterpolateItemValue.hpp>
-#include <mesh/Connectivity.hpp>
-#include <mesh/Mesh.hpp>
-#include <mesh/MeshData.hpp>
-#include <mesh/MeshDataManager.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
#include <utils/Exceptions.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;
-
- public:
- std::shared_ptr<const IMesh>
- mesh() const
- {
- return m_mesh;
- }
-
- 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;
-};
-
std::shared_ptr<IDiscreteFunction>
DiscreteFunctionInterpoler::interpolate() const
{
diff --git a/src/scheme/DiscreteFunctionP0.hpp b/src/scheme/DiscreteFunctionP0.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..4fb0fcc980bc09a71c8668ad16e069ab8c6d0a2d
--- /dev/null
+++ b/src/scheme/DiscreteFunctionP0.hpp
@@ -0,0 +1,65 @@
+#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:
+ 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/IDiscreteFunction.hpp b/src/scheme/IDiscreteFunction.hpp
index eafea0e961c03b4acae7032bf96c3334e8048810..7a7deb08747f74b28e65458a54fe5a775db86b43 100644
--- a/src/scheme/IDiscreteFunction.hpp
+++ b/src/scheme/IDiscreteFunction.hpp
@@ -1,16 +1,16 @@
#ifndef I_DISCRETE_FUNCTION_HPP
#define I_DISCRETE_FUNCTION_HPP
-#include <scheme/DiscreteFunctionType.hpp>
-
#include <memory>
class IMesh;
+class IDiscreteFunctionDescriptor;
class IDiscreteFunction
{
public:
- virtual std::shared_ptr<const IMesh> mesh() const = 0;
+ virtual std::shared_ptr<const IMesh> mesh() const = 0;
+ virtual const IDiscreteFunctionDescriptor& descriptor() const = 0;
IDiscreteFunction() = default;
diff --git a/src/scheme/IDiscreteFunctionDescriptor.hpp b/src/scheme/IDiscreteFunctionDescriptor.hpp
index fcd9fdd1aa8e07973aa9d3503424035d6ce9c646..19d411507e087c0e3262f8e712aa393d96796b20 100644
--- a/src/scheme/IDiscreteFunctionDescriptor.hpp
+++ b/src/scheme/IDiscreteFunctionDescriptor.hpp
@@ -1,6 +1,8 @@
#ifndef I_DISCRETE_FUNCTION_DESCRIPTOR_HPP
#define I_DISCRETE_FUNCTION_DESCRIPTOR_HPP
+#include <scheme/DiscreteFunctionType.hpp>
+
class IDiscreteFunctionDescriptor
{
public:
@@ -15,22 +17,4 @@ class IDiscreteFunctionDescriptor
virtual ~IDiscreteFunctionDescriptor() noexcept = default;
};
-class DiscreteFunctionDescriptorP0 : public IDiscreteFunctionDescriptor
-{
- public:
- DiscreteFunctionType
- type() const final
- {
- return DiscreteFunctionType::P0;
- }
-
- DiscreteFunctionDescriptorP0() noexcept = default;
-
- DiscreteFunctionDescriptorP0(const DiscreteFunctionDescriptorP0& other) = default;
-
- DiscreteFunctionDescriptorP0(DiscreteFunctionDescriptorP0&& other) noexcept = default;
-
- ~DiscreteFunctionDescriptorP0() noexcept = default;
-};
-
#endif // I_DISCRETE_FUNCTION_DESCRIPTOR_HPP
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 06518e70d327fc29760ad6ac2c38142bd2afa138..c4413ddc7c6f06665b4efb7b091548fb26140fab 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -114,6 +114,7 @@ target_link_libraries (unit_tests
PugsMesh
PugsAlgebra
PugsScheme
+ PugsOutput
PugsUtils
kokkos
${PARMETIS_LIBRARIES}
@@ -137,6 +138,7 @@ target_link_libraries (mpi_unit_tests
PugsUtils
PugsLanguageUtils
PugsScheme
+ PugsOutput
PugsUtils
PugsAlgebra
PugsMesh