#ifndef MEDUSA_BITS_UTILS_WRITEVTK_HPP_
#define MEDUSA_BITS_UTILS_WRITEVTK_HPP_

#include <string>
#include <iomanip>
#include <fstream>
#include <medusa/bits/types/Vec.hpp>
#include <medusa/bits/types/Range.hpp>
#include <cmath>
#include <Eigen/Core>

namespace mm
{
    template <int dim>
    void writePntVTK(const char *path, const Eigen::MatrixXd &pnts)
    {
        // std::cout << "Writing " << dim << "*N points to VTK file..." << std::endl;

        std::ofstream out(path);
        out << "# vtk DataFile Version 3.0\n"
               "Volume Mesh\n"
               "ASCII\n"
               "DATASET UNSTRUCTURED_GRID"
            << std::endl;
        out << "POINTS " << pnts.rows() << " float" << std::endl;
        for (int i = 0; i < pnts.rows(); ++i)
        {
            for (int j = 0; j < dim; ++j)
            {
                out << std::setprecision(4) << pnts(i, j) << " ";
            }
            for (int j = dim; j < 3; ++j)
            {
                out << "0 ";
            }
            out << std::endl;
        }
        out << "CELLS " << pnts.rows() << " " << pnts.rows() * (1 + 1) << std::endl;
        for (int i = 0; i < pnts.rows(); ++i)
        {
            out << "1 " << i << std::endl;
        }
        out << "CELL_TYPES " << pnts.rows() << std::endl;
        for (int i = 0; i < pnts.rows(); ++i)
        {
            out << 1 << std::endl;
        }
    }
}; // namespace mm

#endif