Select Git revision
GnuplotWriter1D.cpp 14.10 KiB
#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 <mesh/MeshTraits.hpp>
#include <mesh/MeshVariant.hpp>
#include <utils/Filesystem.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_face_data(const ItemDataType&) const
{
return ItemDataType::item_t == ItemType::face;
}
template <typename ItemDataType>
bool
GnuplotWriter1D::_is_edge_data(const ItemDataType&) const
{
return ItemDataType::item_t == ItemType::edge;
}
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 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 nb_columns = 1;
for (const auto& [name, item_data] : output_named_item_data_set) {
std::visit([&](auto&& value) { nb_columns += _itemDataNbRow(value); }, item_data);
}
return nb_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 (const 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::_write(const MeshType& mesh, const OutputNamedItemDataSet& output_named_item_data_set,
std::optional<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);
}
for (const auto& [name, item_data_variant] : output_named_item_data_set) {
std::visit(
[&, name = name](auto&& item_data) {
if (this->_is_face_data(item_data)) {
std::ostringstream error_msg;
error_msg << "gnuplot_1d_writer does not support face data, cannot save variable \"" << rang::fgB::yellow
<< name << rang::fg::reset << '"';
throw NormalError(error_msg.str());
}
},
item_data_variant);
}
for (const auto& [name, item_data_variant] : output_named_item_data_set) {
std::visit(
[&, name = name](auto&& item_data) {
if (this->_is_edge_data(item_data)) {
std::ostringstream error_msg;
error_msg << "gnuplot_1d_writer does not support edge data, cannot save variable \"" << rang::fgB::yellow
<< name << rang::fg::reset << '"';
throw NormalError(error_msg.str());
}
},
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());
if (not fout) {
std::ostringstream error_msg;
error_msg << "cannot create file \"" << rang::fgB::yellow << _getFilename() << rang::fg::reset << '"';
throw NormalError(error_msg.str());
}
fout.precision(15);
fout.setf(std::ios_base::scientific);
fout << _getDateAndVersionComment();
if (time.has_value()) {
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 MeshVariant&) 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 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) {
using MeshType = mesh_type_t<decltype(p_mesh)>;
if constexpr (MeshType::Dimension == 1) {
this->_write(*p_mesh, output_named_item_data_set, time);
} else if constexpr (MeshType::Dimension == 2) {
std::ostringstream errorMsg;
errorMsg << "gnuplot_1d_writer is not available in dimension " << stringify(MeshType::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());
} else {
throw NormalError("gnuplot format is not available in dimension " + stringify(MeshType::Dimension));
}
},
mesh_v.variant());
}
void
GnuplotWriter1D::_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) {
using MeshType = mesh_type_t<decltype(p_mesh)>;
if constexpr (MeshType::Dimension == 1) {
this->_write(*p_mesh, output_named_item_data_set, {});
} else if constexpr (MeshType::Dimension == 2) {
std::ostringstream errorMsg;
errorMsg << "gnuplot_1d_writer is not available in dimension " << stringify(MeshType::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());
} else {
throw NormalError("gnuplot format is not available in dimension " + stringify(MeshType::Dimension));
}
},
mesh_v.variant());
}