Select Git revision
VTKWriter.cpp
-
Stéphane Del Pino authoredStéphane Del Pino authored
VTKWriter.cpp 39.49 KiB
#include <output/VTKWriter.hpp>
#include <mesh/Connectivity.hpp>
#include <mesh/Mesh.hpp>
#include <mesh/MeshData.hpp>
#include <mesh/MeshDataManager.hpp>
#include <mesh/MeshTraits.hpp>
#include <mesh/MeshVariant.hpp>
#include <mesh/PolynomialMesh.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;
if constexpr (std::is_arithmetic_v<DataType>) {
os << "<PDataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\""
<< item_data.sizeOfArrays() << "\"/>\n";
} else {
throw NormalError("cannot write item array of non arithmetic values");
}
}
}
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()});
if constexpr (is_polygonal_mesh_v<MeshType>) {
output_named_item_data_set.add(NamedItemData{"node_number", mesh.connectivity().nodeNumber()});
}
createDirectoryIfNeeded(m_base_filename);
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) {
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";
if constexpr (is_polynomial_mesh_v<MeshType>) {
fout << "<PDataArray type=\"UInt32\" Name=\"node_number\"/>\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);
}
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";
}
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=\"";
if constexpr (is_polygonal_mesh_v<MeshType>) {
fout << mesh.numberOfNodes();
} else if constexpr (is_polynomial_mesh_v<MeshType>) {
fout << mesh.numberOfNodes() + (mesh.degree() - 1) * mesh.numberOfFaces() +
(mesh.degree() > 2) * mesh.numberOfCells();
} else {
throw NotImplementedError("unexpected mesh type");
}
fout << "\" NumberOfCells=\"" << mesh.numberOfCells() << "\">\n";
fout << "<CellData>\n";
for (const auto& [name, item_value_variant] : output_named_item_data_set) {
std::visit(
[&, var_name = name](auto&& item_value) {
using IVType = std::decay_t<decltype(item_value)>;
using DataType = typename IVType::data_type;
if constexpr (is_item_array_v<IVType> and not std::is_arithmetic_v<DataType>) {
throw NotImplementedError("DiscreteFunctionP0Vector of non arithmetic type");
} else {
return this->_write_cell_data(fout, var_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(
[&, var_name = name](auto&& item_value) {
using IVType = std::decay_t<decltype(item_value)>;
using DataType = typename IVType::data_type;
if constexpr (is_item_array_v<IVType> and not std::is_arithmetic_v<DataType>) {
throw NotImplementedError("DiscreteFunctionP0Vector of non arithmetic type");
} else {
return this->_write_node_data(fout, var_name, item_value, serialize_data_list);
}
},
item_value_variant);
}
if constexpr (is_polynomial_mesh_v<MeshType>) {
size_t number_of_inner_nodes = 0;
if (mesh.degree() > 2) {
const auto cell_type = mesh.connectivity().cellType();
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
switch (cell_type[cell_id]) {
case CellType::Triangle: {
number_of_inner_nodes += (mesh.degree() - 1) * (mesh.degree() - 2) / 2;
break;
}
case CellType::Quadrangle: {
number_of_inner_nodes += (mesh.degree() - 2) * (mesh.degree() - 2);
break;
}
default: {
throw NormalError("invalid cell type for polynomal mesh of degree > 2");
}
}
}
}
Array<int> node_number(mesh.numberOfNodes() + (mesh.degree() - 1) * mesh.numberOfFaces() + number_of_inner_nodes);
node_number.fill(0);
_write_array(fout, "node_number", node_number, serialize_data_list);
}
fout << "</PointData>\n";
fout << "<Points>\n";
using Rd = TinyVector<MeshType::Dimension>;
if constexpr (is_polygonal_mesh_v<MeshType>) {
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);
} else if constexpr (is_polynomial_mesh_v<MeshType>) {
const NodeValue<const Rd>& xr = mesh.xr();
const FaceArray<const Rd>& xl = mesh.xl();
Array<TinyVector<3>> positions(mesh.numberOfNodes() + (mesh.degree() - 1) * mesh.numberOfFaces() +
(mesh.degree() > 2) * mesh.numberOfCells());
parallel_for(
mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
for (unsigned short i = 0; i < MeshType::Dimension; ++i) {
positions[node_id][i] = xr[node_id][i];
}
for (unsigned short i = MeshType::Dimension; i < 3; ++i) {
positions[node_id][i] = 0;
}
});
const size_t number_of_node_per_face = xl.sizeOfArrays();
const size_t number_of_vertices = mesh.numberOfNodes();
parallel_for(
mesh.numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) {
for (size_t i_node = 0; i_node < number_of_node_per_face; ++i_node) {
const size_t index = number_of_vertices + number_of_node_per_face * face_id + i_node;
for (unsigned short i = 0; i < MeshType::Dimension; ++i) {
positions[index][i] = xl[face_id][i_node][i];
}
for (unsigned short i = MeshType::Dimension; i < 3; ++i) {
positions[index][i] = 0;
}
}
});
if (mesh.degree() > 2) {
const size_t cell_begining = mesh.numberOfNodes() + (mesh.degree() - 1) * mesh.numberOfFaces();
auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
for (size_t i_node = 0; i_node < number_of_node_per_face; ++i_node) {
const size_t index = cell_begining + cell_id;
for (unsigned short i = 0; i < MeshType::Dimension; ++i) {
positions[index][i] = xj[cell_id][i];
}
for (unsigned short i = MeshType::Dimension; i < 3; ++i) {
positions[index][i] = 0;
}
}
});
}
_write_array(fout, "Positions", positions, serialize_data_list);
} else {
throw NotImplementedError("unexpected mesh type");
}
fout << "</Points>\n";
fout << "<Cells>\n";
if constexpr (is_polygonal_mesh_v<MeshType>) {
{
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);
}
} else if constexpr (is_polynomial_mesh_v<MeshType>) {
const size_t degree = mesh.degree();
const size_t number_of_vertices = mesh.numberOfNodes();
const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
const auto& cell_type = mesh.connectivity().cellType();
Array<unsigned int> offsets(mesh.numberOfCells());
unsigned int offset = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
const auto& cell_nodes = cell_to_node_matrix[cell_id];
offset += degree * cell_nodes.size();
if (degree > 2) {
if (cell_type[cell_id] == CellType::Triangle) {
offset += (mesh.degree() - 1) * (mesh.degree() - 2) / 2;
} else if (cell_type[cell_id] == CellType::Quadrangle) {
offset += (mesh.degree() - 1) * (mesh.degree() - 1);
}
}
offsets[cell_id] = offset;
}
_write_array(fout, "offsets", offsets, serialize_data_list);
const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
const auto& cell_face_is_reversed_matrix = mesh.connectivity().cellFaceIsReversed();
Array<NodeId::base_type> nodes{offset};
size_t index = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
auto cell_nodes = cell_to_node_matrix[cell_id];
auto cell_faces = cell_to_face_matrix[cell_id];
for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
nodes[index++] = cell_nodes[i_node];
}
if (degree > 1) {
auto cell_face_is_reversed = cell_face_is_reversed_matrix[cell_id];
if (cell_type[cell_id] == CellType::Triangle) {
for (size_t i_face = 0; i_face < cell_faces.size(); ++i_face) {
const size_t i_face_local_number = (i_face + 2) % 3;
if (cell_face_is_reversed[i_face_local_number]) {
for (ssize_t i_face_inner_node = degree - 2; i_face_inner_node >= 0; --i_face_inner_node) {
nodes[index++] =
number_of_vertices + (degree - 1) * cell_faces[i_face_local_number] + i_face_inner_node;
}
} else {
for (size_t i_face_inner_node = 0; i_face_inner_node < degree - 1; ++i_face_inner_node) {
nodes[index++] =
number_of_vertices + (degree - 1) * cell_faces[i_face_local_number] + i_face_inner_node;
}
}
}
if (degree > 2) {
for (size_t i = 0; i < (mesh.degree() - 1) * (mesh.degree() - 2) / 2; ++i) {
nodes[index++] = number_of_vertices + (degree - 1) * mesh.numberOfFaces() + cell_id;
}
}
} else {
for (size_t i_face = 0; i_face < cell_faces.size(); ++i_face) {
if ((i_face >= 2) xor cell_face_is_reversed[i_face]) {
for (ssize_t i_face_inner_node = degree - 2; i_face_inner_node >= 0; --i_face_inner_node) {
nodes[index++] = number_of_vertices + (degree - 1) * cell_faces[i_face] + i_face_inner_node;
}
} else {
for (size_t i_face_inner_node = 0; i_face_inner_node < degree - 1; ++i_face_inner_node) {
nodes[index++] = number_of_vertices + (degree - 1) * cell_faces[i_face] + i_face_inner_node;
}
}
}
if (degree > 2) {
if (cell_type[cell_id] != CellType::Quadrangle) {
std::ostringstream error_msg;
error_msg << "cell type: " << name(cell_type[cell_id])
<< " is not supported for VTK output of polynomial meshes of degree " << mesh.degree()
<< '\n';
throw NormalError(error_msg.str());
}
for (size_t i = 0; i < (mesh.degree() - 1) * (mesh.degree() - 1); ++i) {
nodes[index++] = number_of_vertices + (degree - 1) * mesh.numberOfFaces() + cell_id;
}
}
}
}
}
_write_array(fout, "connectivity", nodes, serialize_data_list);
} else {
throw NotImplementedError("unexpected mesh type");
}
{
Array<int8_t> types(mesh.numberOfCells());
const auto& cell_type = mesh.connectivity().cellType();
const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
const size_t element_degree = [&] {
if constexpr (is_polygonal_mesh_v<MeshType>) {
return 1;
} else if constexpr (is_polynomial_mesh_v<MeshType>) {
return mesh.degree();
} else {
throw NotImplementedError("unexpected mesh type");
}
}();
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
switch (cell_type[j]) {
case CellType::Line: {
types[j] = 3;
break;
}
case CellType::Triangle: {
if constexpr (is_polygonal_mesh_v<MeshType>) {
types[j] = 5;
} else if constexpr (is_polynomial_mesh_v<MeshType>) {
if (element_degree == 1) {
types[j] = 3;
} else if (element_degree == 2) {
types[j] = 22;
} else {
types[j] = 69;
}
}
break;
}
case CellType::Quadrangle: {
if constexpr (is_polygonal_mesh_v<MeshType>) {
types[j] = 9;
} else if constexpr (is_polynomial_mesh_v<MeshType>) {
if (element_degree == 1) {
types[j] = 9;
} else if (element_degree == 2) {
types[j] = 23;
} else {
types[j] = 70;
}
}
break;
}
case CellType::Polygon: {
if constexpr (is_polygonal_mesh_v<MeshType>) {
types[j] = 7;
} else if constexpr (is_polynomial_mesh_v<MeshType>) {
if (element_degree == 1) {
types[j] = 7;
} else if (element_degree == 2) {
types[j] = 36;
} else {
throw NormalError("output for polynomial cells of degree > 2 is not available for polygonal cells");
}
}
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(
[&, var_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
<< var_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
<< var_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(
[&, var_name = name](auto&& item_value) {
using IVType = std::decay_t<decltype(item_value)>;
using DataType = typename IVType::data_type;
if constexpr (is_item_array_v<IVType> and not std::is_arithmetic_v<DataType>) {
throw NotImplementedError("DiscreteFunctionP0Vector of non arithmetic type");
} else {
return this->_write_node_pvtu(fout, var_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(
[&, var_name = name](auto&& item_value) {
using IVType = std::decay_t<decltype(item_value)>;
using DataType = typename IVType::data_type;
if constexpr (is_item_array_v<IVType> and not std::is_arithmetic_v<DataType>) {
throw NotImplementedError("DiscreteFunctionP0Vector of non arithmetic type");
} else {
return this->_write_cell_pvtu(fout, var_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());
}