#include <output/VTKWriter.hpp> #include <mesh/Connectivity.hpp> #include <mesh/Mesh.hpp> #include <mesh/MeshData.hpp> #include <mesh/MeshDataManager.hpp> #include <mesh/MeshVariant.hpp> #include <utils/Filesystem.hpp> #include <utils/Messenger.hpp> #include <utils/RevisionInfo.hpp> #include <utils/Stringify.hpp> #include <utils/pugs_config.hpp> #include <ctime> #include <fstream> #include <iomanip> #include <sstream> #include <unordered_map> class ICharArrayEmbedder { public: ICharArrayEmbedder() = default; ICharArrayEmbedder(const ICharArrayEmbedder&) = default; ICharArrayEmbedder(ICharArrayEmbedder&&) = default; virtual size_t size() const = 0; virtual void write(std::ostream&) const = 0; virtual ~ICharArrayEmbedder() = default; }; template <typename InputDataT> class CharArrayEmbedder : public ICharArrayEmbedder { CastArray<InputDataT, const char> m_char_cast_array; public: size_t size() const final { return m_char_cast_array.size(); } void write(std::ostream& os) const final { os.write(&(m_char_cast_array[0]), m_char_cast_array.size()); } CharArrayEmbedder(Array<InputDataT> array) : m_char_cast_array{array} {} CharArrayEmbedder() = default; CharArrayEmbedder(const CharArrayEmbedder&) = default; CharArrayEmbedder(CharArrayEmbedder&&) = default; ~CharArrayEmbedder() = default; }; class VTKWriter::SerializedDataList { private: std::vector<std::shared_ptr<ICharArrayEmbedder>> m_serialized_data_list; size_t m_offset = 0; public: size_t offset() const { return m_offset; } template <typename DataT> void add(Array<DataT> array) { auto array_data = std::make_shared<CharArrayEmbedder<DataT>>(array); auto size_data = std::make_shared<CharArrayEmbedder<int>>([&] { Array<int> size_array(1); size_array[0] = array_data->size(); return size_array; }()); m_serialized_data_list.push_back(size_data); m_serialized_data_list.push_back(array_data); m_offset += size_data->size() + array_data->size(); } template <typename DataT, ItemType item_type, typename ConnectivityT> void add(const ItemValue<DataT, item_type, ConnectivityT>& item_value) { Array<std::remove_const_t<DataT>> array(item_value.numberOfItems()); parallel_for( item_value.numberOfItems(), PUGS_LAMBDA(ItemIdT<item_type> item_id) { array[item_id] = item_value[item_id]; }); this->add(array); } template <typename DataT, ItemType item_type, typename ConnectivityT> void add(const ItemArray<DataT, item_type, ConnectivityT>& item_array) { Array<std::remove_const_t<DataT>> array(item_array.numberOfItems() * item_array.sizeOfArrays()); parallel_for( item_array.numberOfItems(), PUGS_LAMBDA(ItemIdT<item_type> item_id) { const size_t begin = item_id * item_array.sizeOfArrays(); for (size_t j = 0; j < item_array.sizeOfArrays(); ++j) { array[begin + j] = item_array[item_id][j]; } }); this->add(array); } void write(std::ostream& os) const { for (const auto& serialized_data : m_serialized_data_list) { serialized_data->write(os); } } SerializedDataList() = default; ~SerializedDataList() = default; }; std::string VTKWriter::_getFilenamePVTU() const { std::ostringstream sout; sout << m_base_filename; if (m_period_manager.has_value()) { sout << '.' << std::setfill('0') << std::setw(4) << m_period_manager->nbSavedTimes(); } sout << ".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; } if (m_period_manager.has_value()) { sout << '.' << std::setfill('0') << std::setw(4) << m_period_manager->nbSavedTimes(); } sout << ".vtu"; return sout.str(); } template <typename ItemDataT> void VTKWriter::_write_item_pvtu(std::ofstream& os, const std::string& name, const ItemDataT& item_data) const { if constexpr (is_item_value_v<ItemDataT>) { using DataType = std::decay_t<typename ItemDataT::data_type>; if constexpr (is_tiny_vector_v<DataType> or is_tiny_matrix_v<DataType>) { using SubDataType = typename DataType::data_type; os << "<PDataArray type=\"" << VTKType<SubDataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\"" << DataType{}.dimension() << "\"/>\n"; } else { static_assert(std::is_arithmetic_v<DataType>); os << "<PDataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\"/>\n"; } } else { static_assert(is_item_array_v<ItemDataT>); using DataType = typename ItemDataT::data_type; static_assert(std::is_arithmetic_v<DataType>); os << "<PDataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\"" << item_data.sizeOfArrays() << "\"/>\n"; } } template <typename ItemDataT> void VTKWriter::_write_node_pvtu(std::ofstream& os, const std::string& name, const ItemDataT& item_data) const { if constexpr (ItemDataT::item_t == ItemType::node) { _write_item_pvtu(os, name, item_data); } } template <typename ItemDataT> void VTKWriter::_write_cell_pvtu(std::ofstream& os, const std::string& name, const ItemDataT& item_data) const { if constexpr (ItemDataT::item_t == ItemType::cell) { _write_item_pvtu(os, name, item_data); } } 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" + stringify(sizeof(DataType) * 8); } else { return "Int" + stringify(sizeof(DataType) * 8); } } else if constexpr (std::is_floating_point_v<DataType>) { return "Float" + stringify(sizeof(DataType) * 8); } }(); }; template <typename DataType> void VTKWriter::_write_array(std::ofstream& os, const std::string& name, const Array<DataType>& array, SerializedDataList& serialized_data_list) const { os << "<DataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" offset=\"" << serialized_data_list.offset() << "\" format=\"appended\"/>\n"; serialized_data_list.add(array); } template <size_t N, typename DataType> void VTKWriter::_write_array(std::ofstream& os, const std::string& name, const Array<TinyVector<N, DataType>>& array, SerializedDataList& serialized_data_list) const { os << "<DataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\"" << N << "\" offset=\"" << serialized_data_list.offset() << "\" format=\"appended\"/>\n"; serialized_data_list.add(array); } template <typename DataType> void VTKWriter::_write_node_value(std::ofstream& os, const std::string& name, const NodeValue<const DataType>& item_value, SerializedDataList& serialized_data_list) const { os << "<DataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" offset=\"" << serialized_data_list.offset() << "\" format=\"appended\"/>\n"; serialized_data_list.add(item_value); } 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, SerializedDataList& serialized_data_list) const { os << "<DataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\"" << N << "\" offset=\"" << serialized_data_list.offset() << "\" format=\"appended\"/>\n"; serialized_data_list.add(item_value); } template <size_t N, typename DataType> void VTKWriter::_write_node_value(std::ofstream& os, const std::string& name, const NodeValue<const TinyMatrix<N, N, DataType>>& item_value, SerializedDataList& serialized_data_list) const { os << "<DataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\"" << N * N << "\" offset=\"" << serialized_data_list.offset() << "\" format=\"appended\"/>\n"; serialized_data_list.add(item_value); } template <typename DataType> void VTKWriter::_write_node_value(std::ofstream&, const std::string&, const CellValue<const DataType>&, SerializedDataList&) const {} template <typename ItemDataT> void VTKWriter::_write_item_data(std::ofstream& os, const std::string& name, const ItemDataT& item_data, SerializedDataList& serialized_data_list) const { using DataType = std::decay_t<typename ItemDataT::data_type>; if constexpr (is_item_value_v<ItemDataT>) { if constexpr (std::is_arithmetic_v<DataType>) { os << "<DataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" offset=\"" << serialized_data_list.offset() << "\" format=\"appended\"/>\n"; } else if constexpr (is_tiny_vector_v<DataType> or is_tiny_matrix_v<DataType>) { using SubDataType = typename DataType::data_type; os << "<DataArray type=\"" << VTKType<SubDataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\"" << DataType{}.dimension() << "\" offset=\"" << serialized_data_list.offset() << "\" format=\"appended\"/>\n"; } else { static_assert(std::is_arithmetic_v<DataType>, "unexpected data type"); } serialized_data_list.add(item_data); } else { static_assert(is_item_array_v<ItemDataT>); if constexpr (std::is_arithmetic_v<DataType>) { os << "<DataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\"" << item_data.sizeOfArrays() << "\" offset=\"" << serialized_data_list.offset() << "\" format=\"appended\"/>\n"; } else { static_assert(std::is_arithmetic_v<DataType>, "unexpected data type"); } serialized_data_list.add(item_data); } } template <typename ItemDataT> void VTKWriter::_write_cell_data(std::ofstream& os, const std::string& name, const ItemDataT& item_data, SerializedDataList& serialized_data_list) const { if constexpr (ItemDataT::item_t == ItemType::cell) { _write_item_data(os, name, item_data, serialized_data_list); } } template <typename ItemDataT> void VTKWriter::_write_node_data(std::ofstream& os, const std::string& name, const ItemDataT& item_data, SerializedDataList& serialized_data_list) const { if constexpr (ItemDataT::item_t == ItemType::node) { _write_item_data(os, name, item_data, serialized_data_list); } } template <MeshConcept MeshType> void VTKWriter::_write(const MeshType& mesh, const OutputNamedItemDataSet& given_output_named_item_data_set, std::optional<double> time) const { OutputNamedItemDataSet output_named_item_data_set{given_output_named_item_data_set}; // Adding basic mesh information output_named_item_data_set.add(NamedItemData{"cell_number", mesh.connectivity().cellNumber()}); output_named_item_data_set.add(NamedItemData{"node_number", mesh.connectivity().nodeNumber()}); createDirectoryIfNeeded(m_base_filename); bool has_general_polyhedron = false; { SerializedDataList serialize_data_list; // write VTK files std::ofstream fout(_getFilenameVTU(parallel::rank())); if (not fout) { std::ostringstream error_msg; error_msg << "cannot create file \"" << rang::fgB::yellow << _getFilenameVTU(parallel::rank()) << rang::fg::reset << '"'; throw NormalError(error_msg.str()); } fout << "<?xml version=\"1.0\"?>\n"; fout << _getDateAndVersionComment(); fout << "<VTKFile type=\"UnstructuredGrid\""; #if defined(PUGS_LITTLE_ENDIAN) fout << " byte_order=\"LittleEndian\""; #else fout << " byte_order=\"BigEndian\""; #endif fout << ">\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_data_set) { std::visit([&, name = name]( auto&& item_value) { return this->_write_cell_data(fout, name, item_value, serialize_data_list); }, item_value_variant); } if (parallel::size() > 1) { CellValue<uint8_t> vtk_ghost_type{mesh.connectivity()}; auto cell_is_owned = mesh.connectivity().cellIsOwned(); parallel_for( mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { vtk_ghost_type[cell_id] = cell_is_owned[cell_id] ? 0 : 1; }); _write_array(fout, "vtkGhostType", vtk_ghost_type.arrayView(), serialize_data_list); } fout << "</CellData>\n"; fout << "<PointData>\n"; for (const auto& [name, item_value_variant] : output_named_item_data_set) { std::visit([&, name = name]( auto&& item_value) { return this->_write_node_data(fout, name, item_value, serialize_data_list); }, 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, serialize_data_list); } fout << "</Points>\n"; fout << "<Cells>\n"; { const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix(); _write_array(fout, "connectivity", cell_to_node_matrix.values(), serialize_data_list); } { 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, serialize_data_list); } { 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::Polygon: { types[j] = 7; 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, serialize_data_list); if constexpr (MeshType::Dimension == 3) { 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 << "<AppendedData encoding=\"raw\">\n"; fout << "_"; serialize_data_list.write(fout); fout << '\n'; fout << "</AppendedData>\n"; fout << "</VTKFile>\n"; } if (parallel::rank() == 0) { // write PVTK file std::ofstream fout(_getFilenamePVTU()); if (not fout) { std::ostringstream error_msg; error_msg << "cannot create file \"" << rang::fgB::yellow << _getFilenamePVTU() << rang::fg::reset << '"'; throw NormalError(error_msg.str()); } 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"; if constexpr (MeshType::Dimension == 3) { if (has_general_polyhedron) { fout << "<PDataArray type=\"Int64\" IdType=\"1\" Name=\"faces\"/>\n"; fout << "<PDataArray type=\"Int64\" IdType=\"1\" Name=\"faceoffsets\"/>\n"; } } fout << "</PCells>\n"; for (const auto& [name, item_value_variant] : output_named_item_data_set) { std::visit( [&, name = name](auto&& item_value) { using ItemValueType = std::decay_t<decltype(item_value)>; if constexpr (ItemValueType::item_t == ItemType::edge) { std::ostringstream error_msg; error_msg << "VTK format does not support edge data, cannot save variable \"" << rang::fgB::yellow << name << rang::fg::reset << '"'; throw NormalError(error_msg.str()); } if constexpr (ItemValueType::item_t == ItemType::face) { std::ostringstream error_msg; error_msg << "VTK format does not support face data, cannot save variable \"" << rang::fgB::yellow << name << rang::fg::reset << '"'; throw NormalError(error_msg.str()); } }, item_value_variant); } fout << "<PPointData>\n"; for (const auto& [name, item_value_variant] : output_named_item_data_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_data_set) { std::visit([&, name = name](auto&& item_value) { return this->_write_cell_pvtu(fout, name, item_value); }, item_value_variant); } if (parallel::size() > 1) { fout << "<PDataArray type=\"UInt8\" Name=\"vtkGhostType\" NumberOfComponents=\"1\"/>\n"; } fout << "</PCellData>\n"; for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) { fout << "<Piece Source=" << std::filesystem::path{_getFilenameVTU(i_rank)}.filename() << "/>\n"; } fout << "</PUnstructuredGrid>\n"; fout << "</VTKFile>\n"; } if (parallel::rank() == 0) { // write PVD file const std::string pvd_filename = m_base_filename + ".pvd"; std::ofstream fout(pvd_filename); if (not fout) { std::ostringstream error_msg; error_msg << "cannot create file \"" << rang::fgB::yellow << pvd_filename << rang::fg::reset << '"'; throw NormalError(error_msg.str()); } fout << "<?xml version=\"1.0\"?>\n"; fout << _getDateAndVersionComment(); fout << "<VTKFile type=\"Collection\" version=\"2.0\">\n"; fout << "<Collection>\n"; if (time.has_value()) { for (size_t i_time = 0; i_time < m_period_manager->nbSavedTimes(); ++i_time) { std::ostringstream sout; sout << m_base_filename; sout << '.' << std::setfill('0') << std::setw(4) << i_time << ".pvtu"; fout << "<DataSet timestep=\"" << m_period_manager->savedTime(i_time) << "\" file=" << std::filesystem::path{sout.str()}.filename() << "/>\n"; } fout << "<DataSet timestep=\"" << *time << "\" file=" << std::filesystem::path{_getFilenamePVTU()}.filename() << "/>\n"; } else { fout << "<DataSet file=" << std::filesystem::path{_getFilenamePVTU()}.filename() << "/>\n"; } fout << "</Collection>\n"; fout << "</VTKFile>\n"; } } void VTKWriter::_writeMesh(const MeshVariant& mesh_v) const { OutputNamedItemDataSet output_named_item_data_set; std::visit([&](auto&& p_mesh) { this->_write(*p_mesh, output_named_item_data_set, {}); }, mesh_v.variant()); } void VTKWriter::_writeAtTime(const MeshVariant& mesh_v, const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list, double time) const { OutputNamedItemDataSet output_named_item_data_set = this->_getOutputNamedItemDataSet(named_discrete_data_list); std::visit([&](auto&& p_mesh) { this->_write(*p_mesh, output_named_item_data_set, time); }, mesh_v.variant()); } void VTKWriter::_write(const MeshVariant& mesh_v, const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const { OutputNamedItemDataSet output_named_item_data_set = this->_getOutputNamedItemDataSet(named_discrete_data_list); std::visit([&](auto&& p_mesh) { this->_write(*p_mesh, output_named_item_data_set, {}); }, mesh_v.variant()); }