#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());
}