#include <output/GnuplotWriter1D.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/Stringify.hpp> #include <utils/Demangle.hpp> #include <fstream> #include <iomanip> std::string GnuplotWriter1D::_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 GnuplotWriter1D::_getFilename() 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 << ".gnu"; return sout.str(); } template <typename ItemDataType> bool GnuplotWriter1D::_is_cell_data(const ItemDataType&) const { return ItemDataType::item_t == ItemType::cell; } template <typename ItemDataType> bool GnuplotWriter1D::_is_node_data(const ItemDataType&) const { return ItemDataType::item_t == ItemType::node; } void GnuplotWriter1D::_writePreamble(const OutputNamedItemDataSet& output_named_item_data_set, std::ostream& fout) const { fout << "# list of data\n"; fout << "# 1:x"; uint64_t i = 2; for (const auto& i_named_item_data : output_named_item_data_set) { const std::string name = i_named_item_data.first; const auto& item_data_variant = i_named_item_data.second; std::visit( [&](auto&& item_data) { using ItemDataType = std::decay_t<decltype(item_data)>; using DataType = std::decay_t<typename ItemDataType::data_type>; if constexpr (is_item_value_v<ItemDataType>) { 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{}.numberOfRows(); ++j) { for (size_t k = 0; k < DataType{}.numberOfColumns(); ++k) { fout << ' ' << i++ << ':' << name << '(' << j << ',' << k << ')'; } } } else { throw UnexpectedError("invalid data type"); } } else if constexpr (is_item_array_v<ItemDataType>) { if constexpr (std::is_arithmetic_v<DataType>) { for (size_t j = 0; j < item_data.sizeOfArrays(); ++j) { fout << ' ' << i++ << ':' << name << '[' << j << ']'; } } else { throw UnexpectedError("invalid data type"); } } else { throw UnexpectedError("invalid ItemData type"); } }, item_data_variant); } fout << "\n\n"; } template <typename DataType, ItemType item_type> size_t GnuplotWriter1D::_itemDataNbRow(const ItemValue<DataType, item_type>&) const { if constexpr (std::is_arithmetic_v<DataType>) { return 1; } else if constexpr (is_tiny_vector_v<std::decay_t<DataType>>) { return DataType::Dimension; } else if constexpr (is_tiny_matrix_v<std::decay_t<DataType>>) { return DataType{}.dimension(); } else { throw UnexpectedError("invalid data type for cell value output: " + demangle<DataType>()); } } template <typename DataType, ItemType item_type> size_t GnuplotWriter1D::_itemDataNbRow(const ItemArray<DataType, item_type>& item_array) const { return item_array.sizeOfArrays(); } template <typename MeshType, ItemType item_type> void GnuplotWriter1D::_writeItemDatas(const std::shared_ptr<const MeshType>& mesh, const OutputNamedItemDataSet& output_named_item_data_set, std::ostream& fout) const { using ItemId = ItemIdT<item_type>; const size_t& number_of_columns = [&] { size_t number_of_columns = 1; for (auto [name, item_data] : output_named_item_data_set) { std::visit([&](auto&& value) { number_of_columns += _itemDataNbRow(value); }, item_data); } return number_of_columns; }(); auto is_owned = mesh->connectivity().template isOwned<item_type>(); const size_t& number_of_owned_lines = [&]() { if (parallel::size() > 1) { size_t number_of_owned_items = 0; for (ItemId item_id = 0; item_id < mesh->template numberOf<item_type>(); ++item_id) { if (is_owned[item_id]) { ++number_of_owned_items; } } return number_of_owned_items; } else { return mesh->template numberOf<item_type>(); } }(); Array<double> values{number_of_columns * number_of_owned_lines}; if constexpr (item_type == ItemType::cell) { auto& mesh_data = MeshDataManager::instance().getMeshData(*mesh); const auto& cell_center = mesh_data.xj(); size_t index = 0; for (ItemId item_id = 0; item_id < mesh->template numberOf<item_type>(); ++item_id) { if (is_owned[item_id]) { values[number_of_columns * index++] = cell_center[item_id][0]; } } } else if constexpr (item_type == ItemType::node) { const auto& node_position = mesh->xr(); size_t index = 0; for (ItemId item_id = 0; item_id < mesh->template numberOf<item_type>(); ++item_id) { if (is_owned[item_id]) { values[number_of_columns * index++] = node_position[item_id][0]; } } } else { throw UnexpectedError("invalid item type"); } size_t column_number = 1; for (auto [name, output_item_data] : output_named_item_data_set) { std::visit( [&](auto&& item_data) { using ItemDataT = std::decay_t<decltype(item_data)>; if constexpr (ItemDataT::item_t == item_type) { if constexpr (is_item_value_v<ItemDataT>) { using DataT = std::decay_t<typename ItemDataT::data_type>; size_t index = 0; for (ItemId item_id = 0; item_id < item_data.numberOfItems(); ++item_id) { if (is_owned[item_id]) { if constexpr (std::is_arithmetic_v<DataT>) { values[number_of_columns * index + column_number] = item_data[item_id]; } else if constexpr (is_tiny_vector_v<DataT>) { const size_t k = number_of_columns * index + column_number; for (size_t j = 0; j < DataT::Dimension; ++j) { values[k + j] = item_data[item_id][j]; } } else if constexpr (is_tiny_matrix_v<DataT>) { size_t k = number_of_columns * index + column_number; for (size_t i = 0; i < DataT{}.numberOfRows(); ++i) { for (size_t j = 0; j < DataT{}.numberOfColumns(); ++j) { values[k++] = item_data[item_id](i, j); } } } ++index; } } } else { using DataT = std::decay_t<typename ItemDataT::data_type>; size_t index = 0; for (ItemId item_id = 0; item_id < item_data.numberOfItems(); ++item_id) { if (is_owned[item_id]) { if constexpr (std::is_arithmetic_v<DataT>) { const size_t k = number_of_columns * index + column_number; for (size_t j = 0; j < item_data.sizeOfArrays(); ++j) { values[k + j] = item_data[item_id][j]; } } ++index; } } } } column_number += _itemDataNbRow(item_data); }, output_item_data); } if (parallel::size() > 1) { values = parallel::gatherVariable(values, 0); } if (parallel::rank() == 0) { Assert(values.size() % number_of_columns == 0); std::vector<size_t> line_numbers(values.size() / number_of_columns); for (size_t i = 0; i < line_numbers.size(); ++i) { line_numbers[i] = i; } std::sort(line_numbers.begin(), line_numbers.end(), [&](size_t i, size_t j) { return values[i * number_of_columns] < values[j * number_of_columns]; }); for (auto i_line : line_numbers) { fout << values[i_line * number_of_columns]; for (size_t j = 1; j < number_of_columns; ++j) { fout << ' ' << values[i_line * number_of_columns + j]; } fout << '\n'; } } } template <typename MeshType> void GnuplotWriter1D::_writeAtTime(const std::shared_ptr<const MeshType>& mesh, const OutputNamedItemDataSet& output_named_item_data_set, double time) const { bool has_cell_data = false; for (const auto& [name, item_data_variant] : output_named_item_data_set) { has_cell_data |= std::visit([&](auto&& item_data) { return this->_is_cell_data(item_data); }, item_data_variant); } bool has_node_data = false; for (const auto& [name, item_data_variant] : output_named_item_data_set) { has_node_data |= std::visit([&](auto&& item_data) { return this->_is_node_data(item_data); }, item_data_variant); } if (has_cell_data and has_node_data) { throw NormalError("cannot store both node and cell data in the same gnuplot file"); } std::ofstream fout; if (parallel::rank() == 0) { fout.open(_getFilename()); fout.precision(15); fout.setf(std::ios_base::scientific); fout << _getDateAndVersionComment(); fout << "# time = " << time << "\n\n"; _writePreamble(output_named_item_data_set, fout); } if (has_cell_data) { this->_writeItemDatas<MeshType, ItemType::cell>(mesh, output_named_item_data_set, fout); } else { // has_node_value this->_writeItemDatas<MeshType, ItemType::node>(mesh, output_named_item_data_set, fout); } } void GnuplotWriter1D::writeMesh(const std::shared_ptr<const IMesh>&) const { std::ostringstream errorMsg; errorMsg << "gnuplot_1d_writer does not write meshes\n" << rang::style::bold << "note:" << rang::style::reset << " one can use " << rang::fgB::blue << "gnuplot_writer" << rang::fg::reset << " instead"; throw NormalError(errorMsg.str()); } void GnuplotWriter1D::_writeAtTime( 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); OutputNamedItemDataSet output_named_item_data_set = this->_getOutputNamedItemDataSet(named_discrete_function_list); switch (mesh->dimension()) { case 1: { this->_writeAtTime(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(mesh), output_named_item_data_set, time); break; } case 2: { std::ostringstream errorMsg; errorMsg << "gnuplot_1d_writer is not available in dimension " << stringify(mesh->dimension()) << '\n' << rang::style::bold << "note:" << rang::style::reset << " one can use " << rang::fgB::blue << "gnuplot_writer" << rang::fg::reset << " in dimension 2"; throw NormalError(errorMsg.str()); } default: { throw NormalError("gnuplot format is not available in dimension " + stringify(mesh->dimension())); } } }