// File:        delaunayVTK.cpp
// Description: Converts ASCII output of QHull's Delaunay tessellation to VTK Unstructured Grid
// Author:      Alexei Razoumov
// Date:        2014-Dec
//-----------------------------------------------------------------------------
#include <stdio.h>
#include <iostream>
#include <sstream>
#include <fstream>
#include <vector>
#include <string>
#include <vtkHexahedron.h>
#include <vtkFloatArray.h>
#include <vtkPoints.h>
#include <vtkVersion.h>
#include <vtkSmartPointer.h>
#include <vtkTetra.h>
#include <vtkDoubleArray.h>
#include <vtkPolyhedron.h>
#include <vtkCellArray.h>
#include <vtkDataSetMapper.h>
#include <vtkActor.h>
#include <vtkXMLUnstructuredGridWriter.h>
#include <vtkUnstructuredGrid.h>
#include <vtkPointData.h>
#include <vtkVertexGlyphFilter.h>

using namespace std;

int main() {
  float xyz[3];
  FILE *file;
  int i, j, ip, pid[8], nvertices, ncells, nlocal, **cell = NULL, vertex, linenum, itemnum;
  double *x = NULL, *y = NULL, *z = NULL, *p = NULL, *rho = NULL;
  string item, line;
  ifstream inFile;

  cout << "reading tetrahedra ..." << endl;
  inFile.open("qhullTetrahedra.txt");
  linenum = 0;
  while (getline(inFile, line)) {
    linenum++;
    if (linenum == 1) {
      istringstream linestream(line);
      while (getline(linestream, item, ' '))
	ncells = atoi(item.c_str());
      cell = new int*[ncells];
    }
    if (linenum > 1) {
      istringstream linestream(line);
      ip = linenum - 2;
      cell[ip] = new int[4];
      itemnum = 0;
      while (getline(linestream, item, ' ')) {
	cell[ip][itemnum] = atoi(item.c_str());
	itemnum++;
      }
    }
  }
  inFile.close();

  cout << "reading particles ..." << endl;
  inFile.open("qhullParticles.txt");
  linenum = 0;
  while (getline(inFile, line)) {
    linenum++;
    if (linenum == 2) {
      istringstream linestream(line);
      while (getline(linestream, item, ' '))
	nvertices = atoi(item.c_str());
      x = new double[nvertices];
      y = new double[nvertices];
      z = new double[nvertices];
    }
    if (linenum > 2) {
      istringstream linestream(line);
      itemnum = 0;
      ip = linenum - 3;
      while (getline(linestream, item, ' ')) {
	if (itemnum == 0) x[ip] = atof(item.c_str());
	if (itemnum == 1) y[ip] = atof(item.c_str());
	if (itemnum == 2) z[ip] = atof(item.c_str());
	itemnum++;
      }
    }
  }
  inFile.close();

  cout << "reading 3D fields ..." << endl;
  inFile.open("qhullFields.txt"); 
  linenum = 0;
  p = new double[nvertices];
  rho = new double[nvertices];
  while (getline(inFile, line)) {
    linenum++;
    if (linenum > 2) {
      istringstream linestream(line);
      itemnum = 0;
      ip = linenum - 3;
      while (getline(linestream, item, ' ')) {
	if (itemnum == 0) p[ip] = atof(item.c_str());
	if (itemnum == 1) rho[ip] = atof(item.c_str());
	itemnum++;
      }
    }
  }
  inFile.close();

  vtkSmartPointer points = vtkSmartPointer::New();
  vtkSmartPointer scalars1 = vtkSmartPointer::New();
  scalars1 -> SetName("pressure");
  vtkSmartPointer scalars2 = vtkSmartPointer::New();
  scalars2 -> SetName("density");

  for (i = 0 ; i < nvertices ; i++) {
    points->InsertNextPoint(x[i],y[i],z[i]);
    scalars1->InsertNextValue(p[i]);
    scalars2->InsertNextValue(rho[i]);
  }

  cout << "writing VTK cells ..." << endl;
  vtkSmartPointer unstructuredGrid = vtkSmartPointer::New();
  unstructuredGrid -> SetPoints(points);
  unstructuredGrid -> GetPointData() -> SetScalars(scalars1);
  unstructuredGrid -> GetPointData() -> AddArray(scalars2);
  for (i = 0 ; i < ncells ; i++) {
    vtkIdType ptIds[] = {cell[i][0], cell[i][1], cell[i][2], cell[i][3]};
    unstructuredGrid->InsertNextCell(VTK_TETRA, 4, ptIds);
    if (i%1000000 == 0) cout << i << " out of " << ncells << endl;
  }

  cout << "writing data to file ..." << endl;
  vtkXMLUnstructuredGridWriter* writer = vtkXMLUnstructuredGridWriter::New();
  writer -> SetInputData(unstructuredGrid);
  writer -> SetFileName("a1.vtu");
  writer -> Write();

  return 0;
}