#ifndef TURING_VTK_HH
#define TURING_VTK_HH

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <string>
#include <vector>
#include <mpi.h>
#include "grid.hh"
#include "parallel_grid.hh"

const int TAG_INTRO = 66;
const int TAG_MSG = 67;

class vtk_writer {
    unsigned int cnt;
    const std::string basename;
    const parallel_structured_grid & grid;
    const structured_grid & proc_grid;
    const std::vector<double> & data;
    const unsigned int rank;
public:
    vtk_writer (std::string _basename,
        const parallel_structured_grid & _grid, const structured_grid & _proc_grid,
        const std::vector<double> &_data, unsigned int _rank) :
        cnt(0), basename(_basename), grid(_grid), proc_grid(_proc_grid), data(_data), rank(_rank)
    {}
        
    void write_data(FILE * out, long int origin, long int start_offset, unsigned int * recv_size, double *recv_data) {
        // write data
        unsigned int i[3];
        for (i[2]=0; i[2]<recv_size[2]; i[2]++)
            for (i[1]=0; i[1]<recv_size[1]; i[1]++) {
                long int offset = start_offset+i[2]*grid.global_size(1)*grid.global_size(0)+i[1]*grid.global_size(0);
                fseek(out,origin+sizeof(float)*offset,SEEK_SET);
                for (i[0]=0; i[0]<recv_size[0]; i[0]++)
                {
                    unsigned int idx = (i[2]*recv_size[1]+i[1])*recv_size[0]+i[0];
                    // placement new
                    float f = recv_data[idx];
                    // byte swap
                    char* c = reinterpret_cast<char*>(&f);
                    static const int n = sizeof(float);
                    for (int i=0; i<n/2; i++)
                    {
                        char t = c[n-1-i];
                        c[n-1-i] = c[i];
                        c[i] = t;
                    }
                    fwrite(&f, 1, n, out);
                }
            }
    }

    void write(int c=-1)
    {
        if (c == -1)
            c = cnt;
        if (rank==0) {
            // create filename
            char fname[256];
            snprintf(fname,255,"%s-%05i.vtk", basename.c_str(), c);

            // create file
            FILE *out=fopen(fname, "w");

            // write header
            fprintf(out,"# vtk DataFile Version 2.0\n");
            fprintf(out,"vtk output\n");
            fprintf(out,"BINARY\n ");
            fprintf(out," \n");
            fprintf(out," \n");
            fprintf(out,"DATASET STRUCTURED_POINTS\n");
            fprintf(out,"DIMENSIONS %d %d %d\n",
                grid.global_size(0), grid.global_size(1), grid.global_size(2));
            fprintf(out,"ORIGIN %f %f %f\n", 0.0, 0.0, 0.0);
            fprintf(out,"SPACING %f %f %f\n",
                grid.meshwidth(), grid.meshwidth(), grid.meshwidth());
            fprintf(out," \n");
            fprintf(out,"POINT_DATA %d\n",
                grid.global_size());
            fprintf(out,"SCALARS %s float 1\n", basename.c_str());
            fprintf(out,"LOOKUP_TABLE default\n");
            long int origin = ftell(out);
            // parallel: Daten an Prozess 0 schicken, damit dieser schreiben kann
            unsigned int received = 0;
            while (received < proc_grid.size()-1) {
                received++;
                unsigned int recv_size[3];
                MPI_Status status;
                MPI_Recv(recv_size,3,MPI_INT,MPI_ANY_SOURCE,TAG_INTRO,MPI_COMM_WORLD,&status);
                int count = recv_size[0]*recv_size[1]*recv_size[2];
                int from = status.MPI_SOURCE;
                unsigned int from_pos[3];
                for (int d=0; d<3; d++)
                    from_pos[d] = proc_grid.index(from,d);
                double recv_data[count];
                MPI_Recv(recv_data,count,MPI_DOUBLE,from,TAG_MSG,MPI_COMM_WORLD,NULL);

                // calculate offset of the first vertex
                long int start_offset = 0;
                unsigned int tproc[3] = {0,0,0};
                for (int k=0; k<from_pos[2]; k++) {
                    tproc[2] = k;
                    start_offset += grid.interior_size(proc_grid.index(tproc),2,proc_grid)*grid.global_size(1)*grid.global_size(0);
                }
                tproc[2] = from_pos[2];
                for (int k=0; k<from_pos[1]; k++) {
                    tproc[1] = k;
                    start_offset += grid.interior_size(proc_grid.index(tproc),1,proc_grid)*grid.global_size(0);
                }
                tproc[1] = from_pos[1];
                for (int k=0; k<from_pos[0]; k++) {
                    tproc[0] = k;
                    start_offset += grid.interior_size(proc_grid.index(tproc),0,proc_grid);
                }
                write_data(out,origin,start_offset,recv_size,recv_data);
            }
            // write own data
            unsigned int recv_size[3];
            recv_size[0] = grid.size(0)-grid.overlap(0,0)-grid.overlap(0,1);
            recv_size[1] = grid.size(1)-grid.overlap(1,0)-grid.overlap(1,1);
            recv_size[2] = grid.size(2)-grid.overlap(2,0)-grid.overlap(2,1);
            double interior[recv_size[0]*recv_size[1]*recv_size[2]];
            unsigned int j[3];
            for (j[2]=grid.overlap(2,0); j[2]<grid.size(2)-grid.overlap(2,1); j[2]++)
                for (j[1]=grid.overlap(1,0); j[1]<grid.size(1)-grid.overlap(1,1); j[1]++)
                    for (j[0]=grid.overlap(0,0); j[0]<grid.size(0)-grid.overlap(0,1); j[0]++)
                        interior[(j[2]-grid.overlap(2,0))*recv_size[1]*recv_size[0]+(j[1]-grid.overlap(1,0))*recv_size[0]+j[0]-grid.overlap(0,0)]=data[grid.index(j)];
            write_data(out,origin,0L,recv_size,interior);
            fclose(out);
        } else {
            //retrieve interior
            std::vector<double> interior;
            unsigned int i[3];
            for (i[2]=grid.overlap(2,0); i[2]<grid.size(2)-grid.overlap(2,1); i[2]++)
                for (i[1]=grid.overlap(1,0); i[1]<grid.size(1)-grid.overlap(1,1); i[1]++)
                    for (i[0]=grid.overlap(0,0); i[0]<grid.size(0)-grid.overlap(0,1); i[0]++)
                        interior.push_back(data[grid.index(i)]);
            i[0] = grid.size(0)-grid.overlap(0,0)-grid.overlap(0,1);
            i[1] = grid.size(1)-grid.overlap(1,0)-grid.overlap(1,1);
            i[2] = grid.size(2)-grid.overlap(2,0)-grid.overlap(2,1);
            //send it to 0
            MPI_Send(i,3,MPI_INT,0,TAG_INTRO,MPI_COMM_WORLD);
            MPI_Send(&interior[0],interior.size(),MPI_DOUBLE, 0, TAG_MSG, MPI_COMM_WORLD);
        }
        // wait until all processes have sent their data and 0 has finished writing
        MPI_Barrier(MPI_COMM_WORLD);
        // update timestep
        cnt = c+1;
    }
};
    
#endif // TURING_VTK_HH