// ************************************************************************* //
//                            avth5partFileFormat.C                           //
// ************************************************************************* //

#include <avth5partFileFormat.h>

#include <string>
#include <vector>

#include <vtkFloatArray.h>
#include <vtkRectilinearGrid.h>
#include <vtkStructuredGrid.h>
#include <vtkUnstructuredGrid.h>

#include <avtDatabaseMetaData.h>

#include <Expression.h>

#include <InvalidVariableException.h>
#include <InvalidFilesException.h>
#include <BadIndexException.h>
#include <vtkCellType.h>
#include <vtkPolyData.h>


//h5part specific
#include <stdio.h>
#include <stdlib.h>

#include <fstream>
#include <iomanip>
#include <iostream>



#ifdef PARALLEL_IO
#include <mpi.h>
#include <avtParallel.h>
#endif

using namespace std;

// ****************************************************************************
//  Method: avth5part constructor
//
//  Programmer: cristina -- generated by xml2avt
//  Creation:   Mon Feb 27 13:47:07 PST 2006
//
// ****************************************************************************

avth5partFileFormat::avth5partFileFormat(const char *filename)
    : avtMTMDFileFormat(filename)
{
    // INITIALIZE DATA MEMBERS


    H5PartFile *file;
    fname = filename;

    file = H5PartOpenFile(filename,H5PART_READ);

    if (!file)
      EXCEPTION1(InvalidFilesException, filename);


    int i, j;
    int npoints, npointvars;
    int nspace = 3;

    H5PartSetStep(file,0);
    //points
    npoints= (int) H5PartGetNumParticles(file);
    if (npoints ==  0)
      EXCEPTION1(VisItException, "npoints is zero");
    points.resize(npoints*nspace);
    cout << "constructor: npoints: " << npoints << "\n";

    //point vars
    npointvars= (int) H5PartGetNumDatasets(file); /* get number of datasets in timestep 0 */
    pointvars.resize(npointvars);
    pointvarnames.resize(npointvars);
    cout << "constructor: nvariables: " << npointvars << "\n";

		char name[128];
		h5part_int64_t status;
    for (j=0; j < npointvars; j++){
      status = H5PartGetDatasetName(file,j, name,128);
			if (status != H5PART_SUCCESS){
				EXCEPTION1(VisItException, "could not read a variable name");
			}
      pointvarnames[j] = name;
		}

		H5PartCloseFile(file);
}


// ****************************************************************************
//  Method: avtEMSTDFileFormat::GetNTimesteps
//
//  Purpose:
//      Tells the rest of the code how many timesteps there are in this file.
//
//  Programmer: cristina -- generated by xml2avt
//  Creation:   Mon Feb 27 13:47:07 PST 2006
//
// ****************************************************************************

int
avth5partFileFormat::GetNTimesteps(void)
{
    h5part_int64_t nt;
    H5PartFile *file;
    file = H5PartOpenFile(fname.c_str(),H5PART_READ);
    H5PartSetStep(file,0);
    nt=H5PartGetNumSteps(file); /* get number of steps in file */
    H5PartCloseFile(file);
    return (int) nt;
}


// ****************************************************************************
//  Method: avth5partFileFormat::FreeUpResources
//
//  Purpose:
//      When VisIt is done focusing on a particular timestep, it asks that
//      timestep to free up any resources (memory, file descriptors) that
//      it has associated with it.  This method is the mechanism for doing
//      that.
//
//  Programmer: cristina -- generated by xml2avt
//  Creation:   Mon Feb 27 13:47:07 PST 2006
//
// ****************************************************************************

void
avth5partFileFormat::FreeUpResources(void)
{
}


// ****************************************************************************
//  Method: avth5partFileFormat::PopulateDatabaseMetaData
//
//  Purpose:
//      This database meta-data object is like a table of contents for the
//      file.  By populating it, you are telling the rest of VisIt what
//      information it can request from you.
//
//  Programmer: cristina -- generated by xml2avt
//  Creation:   Mon Feb 27 13:47:07 PST 2006
//
// ****************************************************************************

void
avth5partFileFormat::PopulateDatabaseMetaData(avtDatabaseMetaData *md, int timeState)
{
    //
    // CODE TO ADD A MESH
    //
    // string meshname = ...
    //
    // AVT_RECTILINEAR_MESH, AVT_CURVILINEAR_MESH, AVT_UNSTRUCTURED_MESH,
    // AVT_POINT_MESH, AVT_SURFACE_MESH, AVT_UNKNOWN_MESH
    // avtMeshType mt = AVT_RECTILINEAR_MESH;
    //
    // int nblocks = YOU_MUST_DECIDE;
    // int block_origin = 0;
    // int spatial_dimension = 2;
    // int topological_dimension = 2;
    // float *extents = NULL;
    //
    // Here's the call that tells the meta-data object that we have a mesh:
    //
    // AddMeshToMetaData(md, meshname, mt, extents, nblocks, block_origin,
    //                   spatial_dimension, topological_dimension);
    //

    //
    // CODE TO ADD A SCALAR VARIABLE
    //
    // string mesh_for_this_var = meshname; // ??? -- could be multiple meshes
    // string varname = ...
    //
    // AVT_NODECENT, AVT_ZONECENT, AVT_UNKNOWN_CENT
    // avtCentering cent = AVT_NODECENT;
    //
    //
    // Here's the call that tells the meta-data object that we have a var:
    //
    // AddScalarVarToMetaData(md, varname, mesh_for_this_var, cent);
    //

    //
    // CODE TO ADD A VECTOR VARIABLE
    //
    // string mesh_for_this_var = meshname; // ??? -- could be multiple meshes
    // string varname = ...
    // int vector_dim = 2;
    //
    // AVT_NODECENT, AVT_ZONECENT, AVT_UNKNOWN_CENT
    // avtCentering cent = AVT_NODECENT;
    //
    //
    // Here's the call that tells the meta-data object that we have a var:
    //
    // AddVectorVarToMetaData(md, varname, mesh_for_this_var, cent,vector_dim);
    //

    //
    // CODE TO ADD A TENSOR VARIABLE
    //
    // string mesh_for_this_var = meshname; // ??? -- could be multiple meshes
    // string varname = ...
    // int tensor_dim = 9;
    //
    // AVT_NODECENT, AVT_ZONECENT, AVT_UNKNOWN_CENT
    // avtCentering cent = AVT_NODECENT;
    //
    //
    // Here's the call that tells the meta-data object that we have a var:
    //
    // AddTensorVarToMetaData(md, varname, mesh_for_this_var, cent,tensor_dim);
    //

    //
    // CODE TO ADD A MATERIAL
    //
    // string mesh_for_mat = meshname; // ??? -- could be multiple meshes
    // string matname = ...
    // int nmats = ...;
    // vector<string> mnames;
    // for (int i = 0 ; i < nmats ; i++)
    // {
    //     char str[32];
    //     sprintf(str, "mat%d", i);
    //     -- or -- 
    //     strcpy(str, "Aluminum");
    //     mnames.push_back(str);
    // }
    // 
    // Here's the call that tells the meta-data object that we have a mat:
    //
    // AddMaterialToMetaData(md, matname, mesh_for_mat, nmats, mnames);
    //
    //
    // Here's the way to add expressions:
    //Expression momentum_expr;
    //momentum_expr.SetName("momentum");
    //momentum_expr.SetDefinition("{u, v}");
    //momentum_expr.SetType(Expression::VectorMeshVar);
    //md->AddExpression(&momentum_expr);
    //Expression KineticEnergy_expr;
    //KineticEnergy_expr.SetName("KineticEnergy");
    //KineticEnergy_expr.SetDefinition("0.5*(momentum*momentum)/(rho*rho)");
    //KineticEnergy_expr.SetType(Expression::ScalarMeshVar);
    //md->AddExpression(&KineticEnergy_expr);
    //
		int size;
		size = 1;
#ifdef PARALLEL_IO
		size = PAR_Size();
#endif
	

    if (!points.size()) {
      EXCEPTION1(InvalidFilesException, "Number of points is zero");
    }

		cout << "Populate: size, : " << size << "\n";

    avtMeshMetaData *pmesh = new avtMeshMetaData;

    int dimension = 3;
    pmesh->name = "particles";
    pmesh->originalName = "particles";
    pmesh->meshType = AVT_POINT_MESH;
    pmesh->topologicalDimension = 0;
    pmesh->spatialDimension = dimension;
		pmesh->numBlocks = size;
		pmesh->blockTitle = "subset";
		pmesh->blockPieceName = "subset";
		pmesh->hasSpatialExtents = false; 

    md->Add(pmesh);

    int i;
    for (i=0; i < pointvarnames.size(); i++){
      AddScalarVarToMetaData(md, pointvarnames[i], "particles", AVT_NODECENT);
    }

}


// ****************************************************************************
//  Method: avth5partFileFormat::GetMesh
//
//  Purpose:
//      Gets the mesh associated with this file.  The mesh is returned as a
//      derived type of vtkDataSet (ie vtkRectilinearGrid, vtkStructuredGrid,
//      vtkUnstructuredGrid, etc).
//
//  Arguments:
//      timestate   The index of the timestate.  If GetNTimesteps returned
//                  'N' time steps, this is guaranteed to be between 0 and N-1.
//      domain      The index of the domain.  If there are NDomains, this
//                  value is guaranteed to be between 0 and NDomains-1,
//                  regardless of block origin.
//      meshname    The name of the mesh of interest.  This can be ignored if
//                  there is only one mesh.
//
//  Programmer: cristina -- generated by xml2avt
//  Creation:   Mon Feb 27 13:47:07 PST 2006
//
// ****************************************************************************

vtkDataSet *
avth5partFileFormat::GetMesh(int timestate, int domain, const char *meshname)
{
	cout << "GetMesh domain: " << domain << "\n";

  H5PartFile *file;
  file = H5PartOpenFile(fname.c_str(),H5PART_READ);

  if (!file)
    EXCEPTION1(InvalidFilesException, fname.c_str());

  long int tnpoints, npoints;
	int npointvars;
  int nspace = 3;
	int nprocs = 1;
#ifdef PARALLEL_IO
 	nprocs = PAR_Size();	
#endif

  H5PartSetStep(file,timestate);

  //points
  tnpoints= (int) H5PartGetNumParticles(file);
  h5part_int64_t idStart = (( h5part_int64_t)(tnpoints/nprocs))*domain;
	h5part_int64_t idEnd;	
	if (domain < nprocs-1) 
  	idEnd   = ((h5part_int64_t)(tnpoints/nprocs))*(domain+1);
	else if (domain == nprocs - 1)
  	idEnd   = tnpoints;


  H5PartSetView(file,idStart,idEnd);

  //points
  npoints= (long int) H5PartGetNumParticles(file);
	cout << "GetMesh: npoints for domain " << domain << ": " << npoints << "\n";

  if (strcmp(meshname, "particles") != 0){
    EXCEPTION1(InvalidVariableException, meshname);
  }
  if (npoints ==  0)
    EXCEPTION1(VisItException, "npoints is zero");

  points.resize(npoints*nspace);
  h5part_float64_t *x, *y, *z;
  x = (h5part_float64_t *) malloc(sizeof(h5part_float64_t)*npoints);
  y = (h5part_float64_t *) malloc(sizeof(h5part_float64_t)*npoints);
  z = (h5part_float64_t *) malloc(sizeof(h5part_float64_t)*npoints);


  h5part_int64_t status = H5PART_SUCCESS; 
  status = H5PartReadDataFloat64(file, "x", x);
  if (status != H5PART_SUCCESS)
  EXCEPTION1(VisItException, "Could not read x coordinates");
  status = H5PartReadDataFloat64(file, "y", y);
  if (status != H5PART_SUCCESS)
  EXCEPTION1(VisItException, "Could not read y coordinates");
  status = H5PartReadDataFloat64(file, "z", z);
  if (status != H5PART_SUCCESS)
  EXCEPTION1(VisItException, "Could not read z coordinates");
  for (long int i = 0; i < npoints; i++){
    points[nspace*i] = (float) x[i];
    points[nspace*i+1] = (float) y[i];
    points[nspace*i+2] = (float) z[i];
  }
  free(x);
  free(y);
  free(z);

  H5PartSetView(file,-1, -1);

  vtkPolyData *dataset = vtkPolyData::New();
  vtkPoints *vtkpoints = vtkPoints::New();
  vtkpoints->SetNumberOfPoints((vtkIdType) npoints);

  float *pts = (float *) vtkpoints->GetVoidPointer(0);

  for (long int i=0; i < npoints*nspace; i++){
    pts[i] = points[i];
  }

  dataset->Allocate(npoints*nspace);
  for (long int i=0; i < npoints; i++){
    vtkIdType onevertex = (vtkIdType) i;
    dataset->InsertNextCell(VTK_VERTEX, 1, &onevertex);
  }
  dataset->SetPoints(vtkpoints);
  vtkpoints->Delete();


  H5PartCloseFile(file);
  fprintf(stderr,"proc[%u]:  done\n", domain);

  return dataset;
}


// ****************************************************************************
//  Method: avth5partFileFormat::GetVar
//
//  Purpose:
//      Gets a scalar variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      timestate  The index of the timestate.  If GetNTimesteps returned
//                 'N' time steps, this is guaranteed to be between 0 and N-1.
//      domain     The index of the domain.  If there are NDomains, this
//                 value is guaranteed to be between 0 and NDomains-1,
//                 regardless of block origin.
//      varname    The name of the variable requested.
//
//  Programmer: cristina -- generated by xml2avt
//  Creation:   Mon Feb 27 13:47:07 PST 2006
//
// ****************************************************************************

vtkDataArray *
avth5partFileFormat::GetVar(int timestate, int domain, const char *varname)
{
    //
    // If you have a file format where variables don't apply (for example a
    // strictly polygonal format like the STL (Stereo Lithography) format,
    // then uncomment the code below.
    //
    // EXCEPTION1(InvalidVariableException, varname);
    //

    //
    // If you do have a scalar variable, here is some code that may be helpful.
    //
    // int ntuples = XXX; // this is the number of entries in the variable.
    // vtkFloatArray *rv = vtkFloatArray::New();
    // rv->SetNumberOfTuples(ntuples);
    // for (int i = 0 ; i < ntuples ; i++)
    // {
    //      rv->SetTuple1(i, VAL);  // you must determine value for ith entry.
    // }
    //
    // return rv;
    //

  H5PartFile *file;

	file = H5PartOpenFile(fname.c_str(),H5PART_READ);

  if (!file)
    EXCEPTION1(InvalidFilesException, fname.c_str());

	h5part_int64_t status;
  h5part_int64_t tnpoints, npoints;
	int npointvars;
  int nspace = 3;
	int nprocs = 1;
#ifdef PARALLEL_IO
  nprocs = PAR_Size(); 
#endif

  H5PartSetStep(file,timestate);
  //points
  tnpoints= H5PartGetNumParticles(file);
  //point vars

  char name[64];
  h5part_int64_t *idvar;
  double *data;
	h5part_int64_t idStart = ((h5part_int64_t)(tnpoints/nprocs))*domain;
  h5part_int64_t idEnd; 
  if (domain < nprocs-1) 
    idEnd   = ((h5part_int64_t)(tnpoints/nprocs))*(domain+1);
  else if (domain == nprocs - 1)
    idEnd   = (h5part_int64_t)tnpoints;

  H5PartSetView(file,idStart,idEnd);
	npoints= H5PartGetNumParticles(file);
	cout << "GetVar: npoints for domain " << domain << ": " << npoints << "\n"; 

  for (size_t j=0; j < (size_t)(pointvarnames.size()); j++){
    status = H5PartGetDatasetName(file,j, name,64);
    if (pointvarnames[j] == name) {
      if (strstr(name, "id") != NULL){
        idvar = (h5part_int64_t *) malloc(sizeof(h5part_int64_t)*npoints);
        status = H5PartReadDataInt64(file, name, idvar);
        if (status != H5PART_SUCCESS)
          EXCEPTION1(VisItException, "Could not read dataset");
        pointvars[j].resize(npoints);
        for (size_t i=0; i < (size_t) npoints; i++){
          pointvars[j][i] = (float) idvar[i];
        }
        if (idvar != NULL)
          free(idvar);
      } else {
        data = (h5part_float64_t *) malloc(sizeof(h5part_float64_t)*npoints);
        status = H5PartReadDataFloat64(file, name, data);
        if (status != H5PART_SUCCESS)
          EXCEPTION1(VisItException, "Could not read dataset");
        pointvars[j].resize(npoints);
        for (size_t i=0; i < (size_t)(npoints); i++){
          pointvars[j][i] = (float) data[i];
        }
        if (data != NULL)
          free(data);
      }
    }
  }
	H5PartSetView(file,-1, -1);

  for (int i=0; i < pointvarnames.size(); i++){
    if (pointvarnames[i] == string(varname)){
      vtkFloatArray *scalars = vtkFloatArray::New();
      scalars->SetNumberOfTuples(npoints);
      float *ptr = (float*) scalars->GetVoidPointer(0);
      memcpy(ptr, &pointvars[i][0], sizeof(float)*npoints);
      return scalars;
    }
  }
  H5PartCloseFile(file);
  EXCEPTION1(InvalidVariableException, varname);


}


// ****************************************************************************
//  Method: avth5partFileFormat::GetVectorVar
//
//  Purpose:
//      Gets a vector variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      timestate  The index of the timestate.  If GetNTimesteps returned
//                 'N' time steps, this is guaranteed to be between 0 and N-1.
//      domain     The index of the domain.  If there are NDomains, this
//                 value is guaranteed to be between 0 and NDomains-1,
//                 regardless of block origin.
//      varname    The name of the variable requested.
//
//  Programmer: cristina -- generated by xml2avt
//  Creation:   Mon Feb 27 13:47:07 PST 2006
//
// ****************************************************************************

vtkDataArray *
avth5partFileFormat::GetVectorVar(int timestate, int domain,const char *varname)
{
    //
    // If you have a file format where variables don't apply (for example a
    // strictly polygonal format like the STL (Stereo Lithography) format,
    // then uncomment the code below.
    //
    // EXCEPTION1(InvalidVariableException, varname);
    //

    //
    // If you do have a vector variable, here is some code that may be helpful.
    //
    // int ncomps = YYY;  // This is the rank of the vector - typically 2 or 3.
    // int ntuples = XXX; // this is the number of entries in the variable.
    // vtkFloatArray *rv = vtkFloatArray::New();
    // int ucomps = (ncomps == 2 ? 3 : ncomps);
    // rv->SetNumberOfComponents(ucomps);
    // rv->SetNumberOfTuples(ntuples);
    // float *one_entry = new float[ucomps];
    // for (int i = 0 ; i < ntuples ; i++)
    // {
    //      int j;
    //      for (j = 0 ; j < ncomps ; j++)
    //           one_entry[j] = ...
    //      for (j = ncomps ; j < ucomps ; j++)
    //           one_entry[j] = 0.;
    //      rv->SetTuple(i, one_entry); 
    // }
    //
    // delete [] one_entry;
    // return rv;
    //
		return NULL;
}