VTKWriter_grids_st.hpp

/*
 * VTKWriter_grids_st.hpp
 *
 *  Created on: Sep 3, 2015
 *      Author: Pietro Incardona
 */
 
#ifndef SRC_VTKWRITER_GRIDS_ST_HPP_
#define SRC_VTKWRITER_GRIDS_ST_HPP_
 
 
#include <boost/mpl/pair.hpp>
#include "VTKWriter_grids_util.hpp"
#include "util/util_debug.hpp"
#include "util/convert.hpp"
 
template <typename Grid>
struct cell_grid
{
    openfpm::vector<const Grid *> grids;
 
    comb<Grid::dims> cmb;
 
    cell_grid() {}
 
    cell_grid(const comb<Grid::dims> & cmb)
    :cmb(cmb)
    {}
 
    cell_grid(const cell_grid<Grid> & cg)
    {
        this->operator=(cg);
    }
 
    cell_grid(cell_grid<Grid> && cg)
    {
        this->operator=(cg);
    }
 
    cell_grid<Grid> & operator=(const cell_grid<Grid> & cg)
    {
        cmb = cg.cmb;
        grids = cg.grids;
 
        return *this;
    }
 
    cell_grid<Grid> & operator=(cell_grid<Grid> && cg)
    {
        cmb = cg.cmb;
        grids = cg.grids;
 
        return *this;
    }
};
 
template <typename Grid, typename St>
class ele_g_st
{
public:
 
    typedef Grid value_type;
 
    ele_g_st(){};
 
    ele_g_st(const Point<Grid::dims,St> & offset,
             const Point<Grid::dims,St> & spacing,
             const Box<Grid::dims,St> & dom)
    :offset(offset),spacing(spacing),dom(dom)
    {}
 
    std::string dataset;
    openfpm::vector<cell_grid<Grid>> g;
    Point<Grid::dims,St> offset;
    Point<Grid::dims,St> spacing;
    Box<Grid::dims,size_t> dom;
 
    inline ele_g_st(const ele_g_st & ele)
    {
        this->operator=(ele);
    }
 
    inline ele_g_st(ele_g_st && ele)
    {
        this->operator=(ele);
    }
 
    ele_g_st<Grid,St> & operator=(const ele_g_st & ele)
    {
        dataset = ele.dataset;
        g = ele.g;
        offset = ele.offset;
        spacing = ele.spacing;
        dom = ele.dom;
 
        return *this;
    }
 
    ele_g_st<Grid,St> & operator=(ele_g_st && ele)
    {
        dataset = ele.dataset;
        g = ele.g;
        offset = ele.offset;
        spacing = ele.spacing;
        dom = ele.dom;
 
        return *this;
    }
};
 
template <typename pair>
class VTKWriter<pair,VECTOR_ST_GRIDS>
{
    openfpm::vector< ele_g_st<typename pair::first,typename pair::second> > vg;
 
    size_t get_total()
    {
        size_t tot = 0;
 
        for (size_t i = 0 ; i < vg.size() ; i++)
        {
            for (size_t j = 0 ; j < vg.get(i).g.size() ; j++)
            {
                if (vg.get(i).g.get(j).grids.size() != 0)
                    tot += vg.get(i).g.get(j).grids.get(0)->size();
            }
        }
        return tot;
    }
 
    std::string get_vertex_properties_list()
    {
        std::string v_out;
 
        // write the number of vertex
        v_out += "VERTICES " + std::to_string(get_total()) + " " + std::to_string(get_total() * 2) + "\n";
 
        // return the vertex properties string
        return v_out;
    }
 
    std::string get_point_properties_list()
    {
        std::string v_out;
 
        // write the number of vertex
        v_out += "POINTS " + std::to_string(get_total()) + " float" + "\n";
 
        // return the vertex properties string
        return v_out;
    }
 
    std::string get_point_list()
    {
        std::stringstream v_out;
 
        for (size_t i = 0 ; i < vg.size() ; i++)
        {
            // For each position in the cell
            for (size_t j = 0 ; j < vg.get(i).g.size() ; j++)
            {
                // If there are no grid in this position
                if (vg.get(i).g.get(j).grids.size() == 0)
                    continue;
 
                auto it = vg.get(i).g.get(j).grids.get(0)->getIterator();
 
                Box<pair::first::dims,typename pair::second> dom;
 
                // Calculate the offset of the grid considering its cell position
                Point<pair::first::dims,typename pair::second> middle = vg.get(i).spacing / 2;
                Point<pair::first::dims,typename pair::second> one;
                one.one();
                one = one + toPoint<pair::first::dims,typename pair::second>::convert(vg.get(i).g.get(j).cmb);
                Point<pair::first::dims,typename pair::second> offset = pmul(middle,one) + vg.get(i).offset;
 
                // if there is the next element
                while (it.isNext())
                {
                    Point<pair::first::dims,typename pair::second> p;
                    p = it.get().toPoint();
                    p = pmul(p,vg.get(i).spacing) + offset;
 
                    if (pair::first::dims == 2)
                        v_out << p.toString() << " 0.0" << "\n";
                    else
                        v_out << p.toString() << "\n";
 
                    // increment the iterator
                    ++it;
                }
            }
        }
 
        // return the vertex list
        return v_out.str();
    }
 
    std::string get_prop_components(size_t k)
    {
        std::stringstream v_out;
 
        for (size_t i = 0 ; i < vg.size() ; i++)
        {
            // For each position in the cell
            for (size_t j = 0 ; j < vg.get(i).g.size() ; j++)
            {
                if (k < vg.get(i).g.get(j).grids.size())
                {
                    // get the combination string
                    v_out << vg.get(i).g.get(j).cmb.to_string();
                }
            }
        }
 
        return v_out.str();
    }
 
    std::string get_properties_output(size_t k, std::string prop_name)
    {
        std::stringstream v_out;
 
        // Check if T is a supported format
        // for now we support only scalar of native type
 
        typedef typename boost::mpl::at<typename pair::first::value_type::type,boost::mpl::int_<0>>::type ctype;
 
        std::string type = getType<ctype>();
 
        // if the type is not supported return
        if (type.size() == 0)
        {
#ifndef DISABLE_ALL_RTTI
            std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " the type " << demangle(typeid(ctype).name()) << " is not supported by vtk\n";
#endif
            return "";
        }
 
        std::string prp_cp = get_prop_components(k);
 
        // Create point data properties
        v_out << "SCALARS " << prop_name << "_" << prp_cp << " " << type + "\n";
 
        // Default lookup table
        v_out << "LOOKUP_TABLE default\n";
 
        for (size_t i = 0 ; i < vg.size() ; i++)
        {
            // For each position in the cell
            for (size_t j = 0 ; j < vg.get(i).g.size() ; j++)
            {
                // If there are no grid in this position
                if (vg.get(i).g.get(j).grids.size() == 0)
                    continue;
 
                if (k < vg.get(i).g.get(j).grids.size())
                {
                    // Grid source
                    auto & g_src = *vg.get(i).g.get(j).grids.get(k);
 
                    auto it = g_src.getIterator();
 
                    Box<pair::first::dims,typename pair::second> dom;
 
                    // if there is the next element
                    while (it.isNext())
                    {
                        auto key = it.get();
 
                        v_out << std::to_string(g_src.template get<0>(key))  << "\n";
 
                        // increment the iterator
                        ++it;
                    }
                }
                else
                {
                    // Grid source
                    auto & g_src = *vg.get(i).g.get(j).grids.get(0);
 
                    auto it = g_src.getIterator();
 
                    Box<pair::first::dims,typename pair::second> dom;
 
                    // if there is the next element
                    while (it.isNext())
                    {
                        v_out << "0\n";
 
                        // increment the iterator
                        ++it;
                    }
                }
            }
        }
 
        // return the vertex list
        return v_out.str();
    }
 
    std::string lastProp()
    {
        std::stringstream v_out;
 
        // Create point data properties
        v_out << "SCALARS domain float\n";
 
        // Default lookup table
        v_out << "LOOKUP_TABLE default\n";
 
        for (size_t i = 0 ; i < vg.size() ; i++)
        {
            // For each position in the cell
            for (size_t j = 0 ; j < vg.get(i).g.size() ; j++)
            {
                // If there are no grid in this position
                if (vg.get(i).g.get(j).grids.size() == 0)
                    continue;
 
                auto it = vg.get(i).g.get(j).grids.get(0)->getIterator();
 
                // if there is the next element
                while (it.isNext())
                {
                    if (vg.get(i).dom.isInside(it.get().toPoint()) == true)
                        v_out << "1.0\n";
                    else
                        v_out << "0.0\n";
 
                    // increment the iterator and counter
                    ++it;
                }
            }
        }
 
        return v_out.str();
    }
 
    size_t getMaxFused()
    {
        size_t max = 0;
 
        for (size_t i = 0 ; i < vg.size() ; i++)
        {
            // For each position in the cell
            for (size_t j = 0 ; j < vg.get(i).g.size() ; j++)
            {
                // If there are no grid in this position
                if (vg.get(i).g.get(j).grids.size() > max)
                        max = vg.get(i).g.get(j).grids.size();
            }
        }
 
        return max;
    }
 
    std::string get_vertex_list()
    {
        std::string v_out;
 
        size_t k = 0;
 
        for (size_t i = 0 ; i < vg.size() ; i++)
        {
            // For each position in the cell
            for (size_t j = 0 ; j < vg.get(i).g.size() ; j++)
            {
                // If there are no grid in this position
                if (vg.get(i).g.get(j).grids.size() == 0)
                        continue;
                auto it = vg.get(i).g.get(j).grids.get(0)->getIterator();
 
                while (it.isNext())
                {
                    v_out += "1 " + std::to_string(k) + "\n";
 
                    ++k;
                    ++it;
                }
            }
        }
        // return the vertex list
        return v_out;
    }
 
    std::string get_point_data_header()
    {
        std::string v_out;
 
        v_out += "POINT_DATA " + std::to_string(get_total()) + "\n";
 
        return v_out;
    }
 
    void append_grid(size_t id, const typename pair::first & g, const comb<pair::first::dims> & cmb)
    {
        for(size_t i = 0 ; i < vg.get(id).g.size() ; i++)
        {
            // for each defined grid if exist the combination fuse
            if (cmb == vg.get(id).g.get(i).cmb)
            {
                vg.get(id).g.get(i).grids.add(&g);
                return;
            }
        }
 
        // if the combination does not exist add the grid
        cell_grid< typename pair::first> cg(cmb);
        vg.get(id).g.add(cg);
        vg.get(id).g.last().grids.add(&g);
    }
 
public:
 
    VTKWriter()
    {}
 
    void add(size_t i,
             const typename pair::first & g,
             const Point<pair::first::dims,typename pair::second> & offset,
             const Point<pair::first::dims,typename pair::second> & spacing,
             const Box<pair::first::dims,typename pair::second> & dom,
             const comb<pair::first::dims> & cmb)
    {
        if (i >= vg.size())
            vg.resize(i+1);
 
        vg.get(i).offset = offset;
        vg.get(i).spacing = spacing;
        vg.get(i).dom = dom;
 
        // append the grid
        append_grid(i,g,cmb);
    }
 
    template<int prp = -1> bool write(std::string file,
                                      std::string g_name = "grids",
                                      file_type ft = file_type::ASCII)
    {
        // Header for the vtk
        std::string vtk_header;
        // Point list of the VTK
        std::string point_list;
        // Vertex list of the VTK
        std::string vertex_list;
        // Graph header
        std::string vtk_binary_or_ascii;
        // vertex properties header
        std::string point_prop_header;
        // edge properties header
        std::string vertex_prop_header;
        // Data point header
        std::string point_data_header;
        // Data point
        std::string point_data;
 
        // VTK header
        vtk_header = "# vtk DataFile Version 3.0\n"
                     + g_name + "\n";
 
        // Choose if binary or ASCII
        if (ft == file_type::ASCII)
        {vtk_header += "ASCII\n";}
        else
        {vtk_header += "BINARY\n";}
 
        // Data type for graph is DATASET POLYDATA
        vtk_header += "DATASET POLYDATA\n";
 
        // point properties header
        point_prop_header = get_point_properties_list();
 
        // Get point list
        point_list = get_point_list();
 
        // vertex properties header
        vertex_prop_header = get_vertex_properties_list();
 
        // Get vertex list
        vertex_list = get_vertex_list();
 
        // Get the point data header
        point_data_header = get_point_data_header();
 
        // Get the maximum number of fused grids
        size_t mf = getMaxFused();
 
        // For each property in the vertex type produce a point data
        for (size_t i = 0 ; i < mf ; i++)
            point_data += get_properties_output(i,g_name);
 
        lastProp();
 
 
        // write the file
        std::ofstream ofs(file);
 
        // Check if the file is open
        if (ofs.is_open() == false)
        {std::cerr << "Error cannot create the VTK file: " + file + "\n";}
 
        ofs << vtk_header << point_prop_header << point_list <<
                vertex_prop_header << vertex_list << point_data_header << point_data;
 
        // Close the file
 
        ofs.close();
 
        // Completed succefully
        return true;
    }
};
 
 
#endif /* SRC_VTKWRITER_GRIDS_ST_HPP_ */