// ************************************************************************* //
// 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;
}