/*****************************************************************************
*
* Copyright (c) 2000 - 2017, Lawrence Livermore National Security, LLC
* Produced at the Lawrence Livermore National Laboratory
* LLNL-CODE-442911
* All rights reserved.
*
* This file is  part of VisIt. For  details, see https://visit.llnl.gov/.  The
* full copyright notice is contained in the file COPYRIGHT located at the root
* of the VisIt distribution or at http://www.llnl.gov/visit/copyright.html.
*
* Redistribution  and  use  in  source  and  binary  forms,  with  or  without
* modification, are permitted provided that the following conditions are met:
*
*  - Redistributions of  source code must  retain the above  copyright notice,
*    this list of conditions and the disclaimer below.
*  - Redistributions in binary form must reproduce the above copyright notice,
*    this  list of  conditions  and  the  disclaimer (as noted below)  in  the
*    documentation and/or other materials provided with the distribution.
*  - Neither the name of  the LLNS/LLNL nor the names of  its contributors may
*    be used to endorse or promote products derived from this software without
*    specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT  HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR  IMPLIED WARRANTIES, INCLUDING,  BUT NOT  LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND  FITNESS FOR A PARTICULAR  PURPOSE
* ARE  DISCLAIMED. IN  NO EVENT  SHALL LAWRENCE  LIVERMORE NATIONAL  SECURITY,
* LLC, THE  U.S.  DEPARTMENT OF  ENERGY  OR  CONTRIBUTORS BE  LIABLE  FOR  ANY
* DIRECT,  INDIRECT,   INCIDENTAL,   SPECIAL,   EXEMPLARY,  OR   CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT  LIMITED TO, PROCUREMENT OF  SUBSTITUTE GOODS OR
* SERVICES; LOSS OF  USE, DATA, OR PROFITS; OR  BUSINESS INTERRUPTION) HOWEVER
* CAUSED  AND  ON  ANY  THEORY  OF  LIABILITY,  WHETHER  IN  CONTRACT,  STRICT
* LIABILITY, OR TORT  (INCLUDING NEGLIGENCE OR OTHERWISE)  ARISING IN ANY  WAY
* OUT OF THE  USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
* DAMAGE.
*
*****************************************************************************/

// ************************************************************************* //
//                            avtUintahFileFormat.h                          //
// ************************************************************************* //

#include <avtUintahFileFormat.h>

#include <cstdio>
#include <string>
#include <cmath>
#include <dlfcn.h>
#include <string>
#include <vector>
#include <climits>

#ifdef PARALLEL
#  include <mpi.h>
#endif

#include <vtkPoints.h>
#include <vtkPolyData.h>
#include <vtkCellData.h>
#include <vtkFloatArray.h>
#include <vtkDoubleArray.h>
#include <vtkLongArray.h>
#include <vtkIntArray.h>
#include <vtkRectilinearGrid.h>
#include <vtkUnstructuredGrid.h>

#include <avtParallel.h>
#include <avtDatabase.h>
#include <avtDatabaseMetaData.h>
#include <avtIntervalTree.h>
#include <avtStructuredDomainBoundaries.h>
#include <avtStructuredDomainNesting.h>
#include <avtVariableCache.h>
#include <avtMaterial.h>
#include <avtCallback.h>

#include <TimingsManager.h>
#include <DebugStream.h>
#include <InvalidVariableException.h>
#include <InvalidDBTypeException.h>
#include <DBOptionsAttributes.h>
#include <InvalidFilesException.h>
#include <InstallationFunctions.h>

using std::max;
using std::min;
using std::vector;
using std::string;

//#ifdef PARALLEL
// we only get the option for  serialized reads if compiling the parallel version
//#  define SERIALIZED_READS
//#endif

const double NAN_REPLACE_VAL=1.0E9;

// Macros to stringize the define value so it can be used in a string context.
#define XSTR(x) #x
#define STR(x) XSTR(x)

// ****************************************************************************
//  Method: avtUintahFileFormat constructor
//
//  Programmer: sshankar -- generated by xml2avt
//  Creation:   Tue May 13 19:02:26 PST 2008
//
// ****************************************************************************
avtUintahFileFormat::avtUintahFileFormat(const char *filename,
                                         DBOptionsAttributes* attrs) :
  avtMTMDFileFormat(filename),
  useExtraCells(true),
  dataVariesOverTime(true),
  nodeCentered(false),
  forceMeshReload(true),
  archive(NULL),
  grid(NULL)
{
  int t1 = visitTimer->StartTimer();

  for (int i=0; attrs!=0 && i<attrs->GetNumberOfOptions(); ++i) {
    if (attrs->GetName(i) == "Load extra cells") {
      useExtraCells = attrs->GetBool("Load extra cells");
    }
    else if (attrs->GetName(i) == "Data varies over time") {
      dataVariesOverTime = attrs->GetBool("Data varies over time");
    }
  }

  // Verify that it is a UDA index.xml file:
  // The 2nd line should look like this <Uintah_DataArchive>.
  FILE * fp = fopen( filename, "r" );
  if( fp == NULL ) {
    string error = string( "Failed to open file: " ) + filename;
    EXCEPTION1( InvalidDBTypeException, error.c_str() );
  }

  char line[1024];
  char * result = fgets( line, 1024, fp );
  if( result ) { 
    result = fgets( line, 1024, fp );
  }

  string lineStr = line;
  if( !result || lineStr.find( "<Uintah_DataArchive>" ) == string::npos ) {
    string error = string( filename ) + " does not appear to be a <Uintah_DataArchive>.";
//    printf("here: %s\n", error.c_str());
    EXCEPTION1( InvalidDBTypeException, error.c_str() );
  }
  fclose( fp );

  // This environment variable prevents Core/Thread/Thread_pthreads.cc
  // from adding it's atexit handler, which will kill visit's process
  // when dlclose is called.
  char *tmp = new char[50];
  strcpy(tmp, "THREAD_NO_ATEXIT=1");
  putenv(tmp);
  
  int dlopen_mode = RTLD_NOW;

#ifdef UINTAH_UDA2VIS_LIB
  const char * lib_name = STR( UINTAH_UDA2VIS_LIB );
#else
#error  "UINTAH_UDA2VIS_LIB has not been defined"
#endif

  setenv("THREAD_NO_CATCH_SIGNALS", "True", 1);

  // First try to open the library without any paths - assumes the
  // user has the Uintah library parth in their *_LIBRARY_PATH.
  libHandle = dlopen(lib_name, dlopen_mode);

  // if( libHandle )
  //   std::cerr << __LINE__ << " Uintah lib " << lib_name << std::endl;

  // Try the visit installation library directory.
  if (!libHandle)
  {
    const char *vlibdir = GetVisItLibraryDirectory().c_str();

    if( vlibdir )
    {
      char *lib = (char *) malloc( strlen(vlibdir) + strlen(lib_name) + 8 );

      sprintf( lib, "%s/%s", vlibdir, lib_name );

      libHandle = dlopen(lib, dlopen_mode);

      // if( libHandle )
      //        std::cerr << __LINE__ << " Uintah lib " << lib << std::endl;
    }

    // Try a relative installed path
    if (!libHandle)
    {
      char *lib = (char *) malloc( strlen(vlibdir) + strlen(lib_name) + 12 );

      sprintf( lib, "%s/uintah/%s", vlibdir, lib_name );

      libHandle = dlopen(lib, dlopen_mode);
      
      // if( libHandle )
      //        std::cerr << __LINE__ << " Uintah lib " << lib << std::endl;
    }
  }

  // Did not open so try to open the library using the path from the
  // calling library which would be in plugins/databases.
#ifdef HAVE_WINDOWS_H
  // Not sure to do here???
#else
  if (!libHandle)
  {
    char *pathname = 0;

    Dl_info info;
    if (dladdr(__builtin_return_address(0), &info)) {
      const char *lastslash = strrchr(info.dli_fname,'/');

      if( lastslash )
      {
        int pathLen = strlen(info.dli_fname) - strlen(lastslash);
        
        pathname = (char *) malloc(pathLen+2);
        strncpy( pathname, info.dli_fname, pathLen );
        pathname[pathLen] = '\0';
      }
    }

    if( pathname )
    {
      char *lib = (char *) malloc( strlen(pathname) + strlen(lib_name) + 8 );

      sprintf( lib, "%s/%s", pathname, lib_name );

      libHandle = dlopen(lib, dlopen_mode);

      // if( libHandle )
      //        std::cerr << __LINE__ << " Uintah lib " << lib << std::endl;

      // Try a relative installed path
      if (!libHandle)
      {
        char *lib = (char *) malloc( strlen(pathname) + strlen(lib_name) + 24 );
        
        sprintf( lib, "%s/../../lib/uintah/%s", pathname, lib_name );

        libHandle = dlopen(lib, dlopen_mode);
        
        // if( libHandle )
        //        std::cerr << __LINE__ << " Uintah lib " << lib << std::endl;
      }
      
      // Try a relative installed path
      if (!libHandle)
      {
        char *lib = (char *) malloc( strlen(pathname) + strlen(lib_name) + 16 );
        
        sprintf( lib, "%s/../../lib/%s", pathname, lib_name );
        
        libHandle = dlopen(lib, dlopen_mode);
        
        // if( libHandle )
        //        std::cerr << __LINE__ << " Uintah lib " << lib << std::endl;
      }

      free( pathname );
    }
  }

#endif

  // Did not open so try to open the library using the library path
  // found where the user installed Uintah.
#if __APPLE__
  if (!libHandle)
  {
    char *lib = (char *) malloc( strlen(lib_name) + 32 );
  
    sprintf( lib, "@executable_path/../lib/uintah/%s", lib_name );

    libHandle = dlopen(lib, dlopen_mode);

    // if( libHandle )
    //   std::cerr << __LINE__ << " Uintah lib " << lib << std::endl;
  }

  if (!libHandle)
  {
    char *lib = (char *) malloc( strlen(lib_name) + 32 );
  
    sprintf( lib, "@executable_path/../lib/%s", lib_name );

    libHandle = dlopen(lib, dlopen_mode);

    // if( libHandle )
    //   std::cerr << __LINE__ << " Uintah lib " << lib << std::endl;
  }

#endif

  // Did not open so try to open the library using the Uintah library path
  // found when VisIt was configured.
  if (!libHandle) {

#ifdef UINTAH_LIBRARY_DIR
  const char * lib_path = STR( UINTAH_LIBRARY_DIR );
#else
#error  "UINTAH_LIBRARY_DIR has not been defined"
#endif

    char *lib = (char *) malloc( strlen(lib_path) + strlen(lib_name) + 8 );

    sprintf( lib, "%s/%s", lib_path, lib_name );

    libHandle = dlopen(lib, dlopen_mode);

    // if( libHandle )
    //   std::cerr << __LINE__ << " Uintah lib " << lib << std::endl;
  }

  // Library was never opened so report back an error.
  if (!libHandle) {
    char* errString = dlerror();
    std::string str = "The library " + std::string(lib_name) +
      " could not be opened. The following error was reported: " +
      std::string( errString );
    EXCEPTION1(InvalidDBTypeException, str.c_str());
  }

  // All possible function calls - check here
  char *error;

  openDataArchive = (DataArchive* (*)(const std::string&)) dlsym(libHandle, "openDataArchive");
  if((error = dlerror()) != NULL) {
    EXCEPTION1(InvalidDBTypeException, "The function openDataArchive could not be located in the library!!!");
  }

  closeDataArchive = (void (*)(DataArchive*)) dlsym(libHandle, "closeDataArchive");
  if((error = dlerror()) != NULL) {
    EXCEPTION1(InvalidDBTypeException, "The function closeDataArchive could not be located in the library!!!");
  }

  getGrid = (GridP* (*)(DataArchive*, int)) dlsym(libHandle, "getGrid");
  if((error = dlerror()) != NULL) {
    EXCEPTION1(InvalidDBTypeException, "The function getGrid could not be located in the library!!!");
  }

  releaseGrid = (void (*)(GridP*)) dlsym(libHandle, "releaseGrid");
  if((error = dlerror()) != NULL) {
    EXCEPTION1(InvalidDBTypeException, "The function releaseGrid could not be located in the library!!!");
  }


  getCycleTimes = (std::vector<double> (*)(DataArchive*)) dlsym(libHandle, "getCycleTimes");
  if((error = dlerror()) != NULL) {
    EXCEPTION1(InvalidDBTypeException, "The function getCycleTimes could not be located in the library!!!");
  }

  getTimeStepInfo = (TimeStepInfo* (*)(DataArchive*, GridP*, int, bool)) dlsym(libHandle, "getTimeStepInfo");
  if((error = dlerror()) != NULL) {
    EXCEPTION1(InvalidDBTypeException, "The function getTimeStepInfo could not be located in the library!!!");
  }

  getGridData = (GridDataRaw* (*)(DataArchive*, GridP*, int, int, std::string, int, int, int[3], int[3])) dlsym(libHandle, "getGridData");
  if((error = dlerror()) != NULL) {
    EXCEPTION1(InvalidDBTypeException, "The function getGridData could not be located in the library!!!");
  }

#if (1 <= UINTAH_MAJOR_VERSION && 7 <= UINTAH_MINOR_VERSION )
  variableExists = (bool (*)(DataArchive*, std::string)) dlsym(libHandle, "variableExists");
  if((error = dlerror()) != NULL) {
    EXCEPTION1(InvalidDBTypeException, "The function variableExists could not be located in the library!!!");
  }
#endif

  getParticleData = (ParticleDataRaw* (*)(DataArchive*, GridP*, int, int, std::string, int, int)) dlsym(libHandle, "getParticleData");
  if((error = dlerror()) != NULL) {
    EXCEPTION1(InvalidDBTypeException, "The function getParticleData could not be located in the library!!!");
  }

  getParticlePositionName = (std::string (*)(DataArchive*)) dlsym(libHandle, "getParticlePositionName");
  if((error = dlerror()) != NULL) {
    EXCEPTION1(InvalidDBTypeException, "The function getParticlePositionName could not be located in the library!!!");
  }

  // use the folder name, not the index.xml file name to open the archive
  string folder(filename);
  size_t found = folder.find_last_of("/");
  folder = folder.substr(0, found);
  archive = (*openDataArchive)(folder);

  // timestep times
  int t2 = visitTimer->StartTimer();

  debug5<<__FILE__<<"  "<<__FUNCTION__<<"  "<< __LINE__<< std::endl;
  cycleTimes = (*getCycleTimes)(archive);
  debug5<<__FILE__<<"  "<<__FUNCTION__<<"  "<< __LINE__<< std::endl;
  
  visitTimer->StopTimer(t2, "avtUintahFileFormat::GetVar getCycleTimes");

  // haven't loaded any timestep data yet
  stepInfo = NULL;  
  currTimeStep = -1;

  visitTimer->StopTimer(t1, "avtUintahFileFormat::avtUintahFileFormat");
}


// Destructor
avtUintahFileFormat::~avtUintahFileFormat()
{
  if (grid)
    (*releaseGrid)(grid);

  if (archive)
    (*closeDataArchive)(archive);

  if (stepInfo)
    delete stepInfo;

  dlclose(libHandle);
}


// ****************************************************************************
//  Method: avtEMSTDFileFormat::GetNTimesteps
//
//  Purpose:
//      Tells the rest of the code how many timesteps there are in this file.
//
//  Programmer: sshankar -- generated by xml2avt
//  Creation:   Tue May 13 19:02:26 PST 2008
//
// ****************************************************************************

int
avtUintahFileFormat::GetNTimesteps(void)
{
  return cycleTimes.size();
}


// ****************************************************************************
// Method: avtUintahFileForma::GetTime
//
// Purpose: 
//   Get the time.
//
// Programmer: sshankar 
// Creation:   Fri Feb 6 15:31 MST 2009 
//
// ****************************************************************************

double 
avtUintahFileFormat::GetTime(int ts)
{
  return cycleTimes[ts];
}


// ****************************************************************************
// Method: avtUintahFileForma::ActivateTimestep
//
// Purpose: 
//   Get ready to read data for the given timestep
//
// ****************************************************************************

void
avtUintahFileFormat::ActivateTimestep(int ts)
{
  if (currTimeStep == ts)
    return;

  int t1 = visitTimer->StartTimer();

  // get the uda grid for the new timestep
  int t2 = visitTimer->StartTimer();

  if (grid)
    (*releaseGrid)(grid);

  debug5<<__FILE__<<"  "<<__FUNCTION__<<"  "<< __LINE__<< std::endl;
  grid = (*getGrid)(archive, ts);
  debug5<<__FILE__<<"  "<<__FUNCTION__<<"  "<< __LINE__<< std::endl;

  visitTimer->StopTimer(t2, "avtUintahFileFormat::ActivateTimestep getGrid");

  // get the time step info for the new timestep
  int t3 = visitTimer->StartTimer();

  if (stepInfo)
    delete stepInfo;

  debug5<<__FILE__<<"  "<<__FUNCTION__<<"  "<< __LINE__<< std::endl;
  stepInfo = (*getTimeStepInfo)(archive, grid, ts, useExtraCells);
  debug5<<__FILE__<<"  "<<__FUNCTION__<<"  "<< __LINE__<< std::endl;

  currTimeStep = ts; 
  forceMeshReload = true;

  visitTimer->StopTimer(t3, "avtUintahFileFormat::ActivateTimestep getTimeStepInfo");

  visitTimer->StopTimer(t1, "avtUintahFileFormat::ActivateTimestep");
}



// ****************************************************************************
//  Method: avtUintahFileFormat::AddExpressionsToMetadata
// ****************************************************************************
void avtUintahFileFormat::AddExpressionsToMetadata(avtDatabaseMetaData *md) {
  string home(getenv("HOME"));
  string fname = home +"/.visit/udaExpressions.txt";
  FILE *f = fopen(fname.c_str(), "r");
  if (!f) {
    debug5<<"couldn't open uda expressions file"<<endl;
    return;
  }


  int line = 0;
  char cline[2048];
  while (fgets(cline, 2048, f)) {
    line++;

    string line(cline);
    if (line[0] == '#' ||
        line.length() == 0 ||
        line == "\n" || line=="\r\n")
      continue;

    // first token is the expression name
    int space = line.find(" ");
    if (space < 0) {
      debug5<<"uda expression syntax error on line "<<line<<endl;
      continue;
    }
    string name = line.substr(0, space);
    line = line.substr(space+1);

    // second token is the expression type
    space = line.find(" ");
    if (space < 0) {
      debug5<<"uda expression syntax error on line "<<line<<endl;
      continue;
    }
    string type = line.substr(0, space);
    line = line.substr(space+1);

    Expression::ExprType etype = Expression::Unknown;
    if (type == "scalar")
      etype = Expression::ScalarMeshVar;
    if (type == "vector")
      etype = Expression::VectorMeshVar;
    if (type == "tensor")
      etype = Expression::TensorMeshVar;
    if (type == "symtensor")
      etype = Expression::SymmetricTensorMeshVar;
    if (type == "array")
      etype = Expression::ArrayMeshVar;
    if (type == "curve")
      etype = Expression::CurveMeshVar;

    if (etype == Expression::Unknown) {
      debug5<<"uda expression - unknown type on line "<<line<<endl;
      continue;
    }


    Expression e;
    e.SetName(name);
    e.SetType(etype);
    e.SetDefinition(line);
    md->AddExpression(&e);
  }


  fclose(f);
}


// ****************************************************************************
//  Method: getBounds
//
//  Purpose:
//   Returns the bounds for the given patch of the specified mesh 
//   based on periodicity and type.
//
//  Node centered data uses the same mesh as cell centered, 
//  but face centered meshes need an extra value for one axis,
//  unless they are periodic on that axis.
//
//  use patch_id=-1 to query all patches.
//
// ****************************************************************************
void getBounds(int low[3], int high[3], const string meshName, const LevelInfo &levelInfo,int patch_id=-1)
{
  levelInfo.getBounds(low,high,meshName,patch_id);
  
  //debug5<<"getBounds("<<meshName<<",id="<<patch_id<<")=["<<low[0]<<","<<low[1]<<","<<low[2]<<"] to ["<<high[0]<<","<<high[1]<<","<<high[2]<<"]\n";
}

// ****************************************************************************
//  Method: avtUintahFileFormat::ReadMetaData
//
//  Purpose:
//      Does the actual work for PopulateMetaData().
//
//     ADDS MESHES
// ****************************************************************************
void
avtUintahFileFormat::ReadMetaData(avtDatabaseMetaData *md, int timeState)
{
  int t1 = visitTimer->StartTimer();

  ActivateTimestep(timeState);

  int numLevels = stepInfo->levelInfo.size();

  int totalPatches = 0;
  for (int i = 0; i < numLevels; i++)
    totalPatches += stepInfo->levelInfo[i].patchInfo.size();

  vector<int> groupIds(totalPatches);
  vector<string> pieceNames(totalPatches);

  for (int i = 0; i < totalPatches; i++) {
    char tmpName[64];
    int level, local_patch;

    GetLevelAndLocalPatchNumber(i, level, local_patch);
    sprintf(tmpName,"level%d, patch%d", level, local_patch);

    groupIds[i] = level;
    pieceNames[i] = tmpName;
  }

  // Compute the bounding box of the mesh from the grid indices of
  // level 0
  LevelInfo &levelInfo = stepInfo->levelInfo[0];

  // Don't add proc id unless CC_Mesh or NC_Mesh exists (some only
  // have SFCk_MESH)
  bool addProcId = false;
  string mesh_for_procid("NC_Mesh");

  // Grid meshes are shared between materials, and particle meshes are
  // shared between variables - keep track of what has been added so
  // they're only added once
  std::set<string> meshes_added;

  // If a variable exists in multiple materials, we don't want to add
  // it more than once to the meta data - it can mess up visit's
  // expressions variable lists.
  std::set<string> mesh_vars_added;

  // Get CC bounds
  int low[3],high[3];
  getBounds(low,high,"CC_Mesh",levelInfo);

  // This can be done once for everything because the spatial range is
  // the same for all meshes
  double box_min[3] = { levelInfo.anchor[0] + low[0] * levelInfo.spacing[0],
                        levelInfo.anchor[1] + low[1] * levelInfo.spacing[1],
                        levelInfo.anchor[2] + low[2] * levelInfo.spacing[2] };
  double box_max[3] = { levelInfo.anchor[0] + high[0] * levelInfo.spacing[0],
                        levelInfo.anchor[1] + high[1] * levelInfo.spacing[1],
                        levelInfo.anchor[2] + high[2] * levelInfo.spacing[2] };
  //debug5<<"box_min/max=["<<box_min[0]<<","<<box_min[1]<<","<<box_min[2]<<"] to ["<<box_max[0]<<","<<box_max[1]<<","<<box_max[2]<<"]\n";

  int logical[3];
  for (int i=0; i<3; ++i)
    logical[i] = high[i]-low[i];
  //debug5 <<"logical: "<< logical[0] << ", "<< logical[1] << ", "<< logical[2]<<endl;

  for (int i=0; i<(int)stepInfo->varInfo.size(); ++i) {
    if (stepInfo->varInfo[i].type.find("ParticleVariable") == string::npos) {
      string varname = stepInfo->varInfo[i].name;
      string vartype = stepInfo->varInfo[i].type;

      string mesh_for_this_var;
      avtCentering cent=AVT_ZONECENT;

      if (vartype.find("NC") != string::npos) {
        cent = AVT_NODECENT;
        mesh_for_this_var.assign("NC_Mesh"); 
        addProcId=true;
      }  
      else if (vartype.find("CC") != string::npos) {  
        cent = AVT_ZONECENT;
        mesh_for_this_var.assign("CC_Mesh");
        addProcId=true;
        mesh_for_procid=mesh_for_this_var;
      }
      else if (vartype.find("SFC") != string::npos) { 
        cent = AVT_ZONECENT;

        if (vartype.find("SFCX") != string::npos)               
          mesh_for_this_var.assign("SFCX_Mesh");
        else if (vartype.find("SFCY") != string::npos)          
          mesh_for_this_var.assign("SFCY_Mesh");
        else if (vartype.find("SFCZ") != string::npos)          
          mesh_for_this_var.assign("SFCZ_Mesh");
      }  
      else
        debug5<<"unknown vartype: "<<vartype<<endl;

      if (meshes_added.find(mesh_for_this_var)==meshes_added.end()) {
        avtMeshMetaData *mesh = new avtMeshMetaData;

        mesh->name = mesh_for_this_var;
        mesh->meshType = AVT_AMR_MESH;
        mesh->topologicalDimension = 3;
        mesh->spatialDimension = 3;

        mesh->numBlocks = totalPatches;
        mesh->blockTitle = "patches";
        mesh->blockPieceName = "patch";
        mesh->numGroups = numLevels;
        mesh->groupTitle = "levels";
        mesh->groupPieceName = "level";
        mesh->blockNames = pieceNames;
        mesh->containsExteriorBoundaryGhosts = false;

        mesh->hasSpatialExtents = true; 
        mesh->minSpatialExtents[0] = box_min[0];
        mesh->maxSpatialExtents[0] = box_max[0];
        mesh->minSpatialExtents[1] = box_min[1];
        mesh->maxSpatialExtents[1] = box_max[1];
        mesh->minSpatialExtents[2] = box_min[2];
        mesh->maxSpatialExtents[2] = box_max[2];

        mesh->hasLogicalBounds = true;
        mesh->logicalBounds[0] = logical[0];
        mesh->logicalBounds[1] = logical[1];
        mesh->logicalBounds[2] = logical[2];

        md->Add(mesh);
        meshes_added.insert(mesh_for_this_var);
      }

      // Add meshvars
      for (int j=0; j<(int)stepInfo->varInfo[i].materials.size(); j++) {
        char buffer[128];
        string newVarname = varname;
        sprintf(buffer, "%d", stepInfo->varInfo[i].materials[j]);
        newVarname.append("/");
        newVarname.append(buffer);

        if (mesh_vars_added.find(mesh_for_this_var+newVarname)==mesh_vars_added.end()) {
          mesh_vars_added.insert(mesh_for_this_var+newVarname);

          if (vartype.find("Vector") != string::npos)
            AddVectorVarToMetaData(md, newVarname, mesh_for_this_var, cent, 3); // 3 -> vector dimension
          else if (vartype.find("Matrix3") != string::npos)
            AddTensorVarToMetaData(md, newVarname, mesh_for_this_var, cent, 9); // 9 -> tensor 
          else if (vartype.find("Stencil7") != string::npos)
            AddVectorVarToMetaData(md, newVarname, mesh_for_this_var, cent, 7); // 7 -> vector
          else if (vartype.find("Stencil4") != string::npos)
            AddVectorVarToMetaData(md, newVarname, mesh_for_this_var, cent, 4); // 4 -> vector
          else 
            AddScalarVarToMetaData(md, newVarname, mesh_for_this_var, cent);
        }
      }
    }   
  }

  // Add a proc id enum variable
  if (addProcId)
  {
    avtScalarMetaData *scalar = new avtScalarMetaData();

    scalar->name = "proc_id";
    scalar->meshName = mesh_for_procid;
    scalar->centering = AVT_ZONECENT;
    scalar->hasDataExtents = false;
    scalar->treatAsASCII = false;
    md->Add(scalar);
  }
  

  // Nothing needs to be modifed for particle data, as they exist only
  // on a single level
  for (int i=0; i<(int)stepInfo->varInfo.size(); ++i) {
    if (stepInfo->varInfo[i].type.find("ParticleVariable") != string::npos) {
      string varname = stepInfo->varInfo[i].name;
      string vartype = stepInfo->varInfo[i].type;

      // j=-1 -> all materials (*)
      for (int j=-1; j<(int)stepInfo->varInfo[i].materials.size(); ++j) {
        string mesh_for_this_var = string("Particle_Mesh/");
        string newVarname = varname+"/";

        if (j >= 0) {
          char buffer[128];
          sprintf(buffer, "%d", stepInfo->varInfo[i].materials[j]);
          mesh_for_this_var.append(buffer);
          newVarname.append(buffer);
        }
        else {
          mesh_for_this_var.append("*");
          newVarname.append("*");
        }

        if (meshes_added.find(mesh_for_this_var)==meshes_added.end()) {

          avtMeshMetaData *mesh = new avtMeshMetaData;

          mesh->name = mesh_for_this_var;
          mesh->meshType = AVT_POINT_MESH;
          mesh->topologicalDimension = 0;
          mesh->spatialDimension = 3;

          mesh->numBlocks = totalPatches;
          mesh->blockTitle = "patches";
          mesh->blockPieceName = "patch";
          mesh->numGroups = numLevels;
          mesh->groupTitle = "levels";
          mesh->groupPieceName = "level";
          mesh->blockNames = pieceNames;

          mesh->hasSpatialExtents = true; 
          mesh->minSpatialExtents[0] = box_min[0];
          mesh->maxSpatialExtents[0] = box_max[0];
          mesh->minSpatialExtents[1] = box_min[1];
          mesh->maxSpatialExtents[1] = box_max[1];
          mesh->minSpatialExtents[2] = box_min[2];
          mesh->maxSpatialExtents[2] = box_max[2];

          mesh->hasLogicalBounds = true;
          mesh->logicalBounds[0] = logical[0];
          mesh->logicalBounds[1] = logical[1];
          mesh->logicalBounds[2] = logical[2];

          md->Add(mesh); 
          meshes_added.insert(mesh_for_this_var);
        }

        if (mesh_vars_added.find(mesh_for_this_var+newVarname)==mesh_vars_added.end()) {
          mesh_vars_added.insert(mesh_for_this_var+newVarname);

          avtCentering cent = AVT_NODECENT;
          // 3 -> vector dimension
          if ((vartype.find("Vector") != string::npos) ||
              (vartype.find("Point") != string::npos))
            AddVectorVarToMetaData(md, newVarname, mesh_for_this_var, cent, 3);
          // 9 -> tensor
          else if (vartype.find("Matrix3") != string::npos)
            AddTensorVarToMetaData(md, newVarname, mesh_for_this_var, cent, 9);
          // 7 -> vector
          else if (vartype.find("Stencil7") != string::npos)
            AddVectorVarToMetaData(md, newVarname, mesh_for_this_var, cent, 7);
          // 4 -> vector
          else if (vartype.find("Stencil4") != string::npos)
            AddVectorVarToMetaData(md, newVarname, mesh_for_this_var, cent, 4);
          else
            AddScalarVarToMetaData(md, newVarname, mesh_for_this_var, cent);
        }
      }
    }   
  }
  
  md->AddGroupInformation(numLevels, totalPatches, groupIds);
  md->AddDefaultSILRestrictionDescription(std::string("!TurnOnAll"));

  // Set the cycles and times.
  md->SetCyclesAreAccurate(true);

  std::vector<int> cycles;

  cycles.resize( cycleTimes.size() );

  for(int i=0; i<(int)cycleTimes.size(); ++i )
    cycles[i] = i;

  md->SetCycles( cycles );

  md->SetTimesAreAccurate(true);
  md->SetTimes( cycleTimes );
  
  AddExpressionsToMetadata(md);

  visitTimer->StopTimer(t1, "avtUintahFileFormat::ActivateTimestep");
}


// ****************************************************************************
//  Method: avtUintahFileFormat::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: sshankar -- generated by xml2avt
//  Creation:   Tue May 13 19:02:26 PST 2008
//
// ****************************************************************************
void
avtUintahFileFormat::PopulateDatabaseMetaData(avtDatabaseMetaData *md, int timeState)
{
#ifdef SERIALIZED_READS

  int numProcs, rank;
  int msg = 128, tag = 256;
  MPI_Status status;

  MPI_Comm_size(VISIT_MPI_COMM, &numProcs);
  MPI_Comm_rank(VISIT_MPI_COMM, &rank);
  //debug5 << "Proc: " << rank << " sent to mdserver" << endl;  

  if (rank == 0) {
    ReadMetaData(md, timeState);
    MPI_Send(&msg, 1, MPI_INT, 1, tag, VISIT_MPI_COMM);
  }
  else {
    MPI_Recv(&msg, 1, MPI_INT, rank - 1, tag, VISIT_MPI_COMM, &status);
    if (msg == 128 && tag == 256) {
      ReadMetaData(md, timeState);
      if (rank < (numProcs - 1))
        MPI_Send(&msg, 1, MPI_INT, rank + 1, tag, VISIT_MPI_COMM);
    }
  }
#else      
  ReadMetaData(md, timeState);
#endif
}


// ****************************************************************************
//  Method: avtUintahFileFormat::GetLevelAndLocalPatchNumber
//
//  Purpose:
//      Translates the global patch identifier to a refinement level and patch
//      number local to that refinement level.
//  
//  Programmer: sshankar, taken from implementation of the plugin, CHOMBO
//  Creation:   May 20, 2008
//
// ****************************************************************************
void
avtUintahFileFormat::GetLevelAndLocalPatchNumber(int global_patch, 
                                                 int &level, int &local_patch)
{
//  int t1 = visitTimer->StartTimer();

  int num_levels = stepInfo->levelInfo.size();
  int num_patches = 0;
  int tmp = global_patch;
  level = 0;
  while (level < num_levels) {
    num_patches = stepInfo->levelInfo[level].patchInfo.size();
    if (tmp < num_patches)
      break;
    tmp -= num_patches;
    level++;
  }
  local_patch = tmp;

//  visitTimer->StopTimer(t1, "avtUintahFileFormat::GetLevelAndLocalPatchNumber");
}


// ****************************************************************************
//  Method: avtUintahFileFormat::GetGlobalDomainNumber
//
//  Purpose:
//      Translates the level and local patch number into a global patch id.
//  
// ****************************************************************************
int avtUintahFileFormat::GetGlobalDomainNumber(int level, int local_patch) {
  int g=0;
  for (int l=0; l<level; l++)
    g += stepInfo->levelInfo[l].patchInfo.size();
  g += local_patch;

  return g;
}


// ****************************************************************************
//  Method: avtUintahFileFormat::CalculateDomainNesting
//
//  Purpose:
//      Calculates two important data structures.  One is the structure domain
//      nesting, which tells VisIt how the AMR patches are nested, which allows
//      VisIt to ghost out coarse zones that are refined by smaller zones.
//      The other structure is the rectilinear domain boundaries, which tells
//      VisIt which patches are next to each other, allowing VisIt to create
//      a layer of ghost zones around each patch.  Note that this only works
//      within a refinement level, not across refinement levels.
//  
//
// NOTE: The cache variable for the mesh MUST be called "any_mesh",
// which is a problem when there are multiple meshes or one of them is
// actually named "any_mesh" (see
// https://visitbugs.ornl.gov/issues/52). Thus, for each mesh we keep
// around our own cache variable and if this function finds it then it
// just uses it again instead of recomputing it.
//
// ****************************************************************************
void
avtUintahFileFormat::CalculateDomainNesting(int timestate, const std::string &meshname)
{
//#warning "FIX ME - IS THIS CODE VALID ??????"

#ifdef MDSERVER

    return;

#else
  int t1 = visitTimer->StartTimer();

  //lookup mesh in our cache and if it's not there, compute it
  if (*this->mesh_domains[meshname]==NULL || forceMeshReload == true)
  {
    //
    // Calculate some info we will need in the rest of the routine.
    //
    int num_levels = stepInfo->levelInfo.size();
    int totalPatches = 0;
    for (int level = 0 ; level < num_levels ; level++)
      totalPatches += stepInfo->levelInfo[level].patchInfo.size();

    //
    // Now set up the data structure for patch boundaries.  The data 
    // does all the work ... it just needs to know the extents of each patch.
    //
    avtRectilinearDomainBoundaries *rdb = new avtRectilinearDomainBoundaries(true);
    //debug5<<"Calculating avtRectilinearDomainBoundaries for "<<meshname<<" mesh ("<<rdb<<").\n";
    rdb->SetNumDomains(totalPatches);

    for (int patch = 0 ; patch < totalPatches ; patch++) {
      int my_level, local_patch;
      GetLevelAndLocalPatchNumber(patch, my_level, local_patch);

      PatchInfo &patchInfo = stepInfo->levelInfo[my_level].patchInfo[local_patch];

      int low[3],high[3];
      patchInfo.getBounds(low,high,meshname);

      int e[6] = { low[0], high[0],
                   low[1], high[1],
                   low[2], high[2] };
      //debug5<<"\trdb->SetIndicesForAMRPatch("<<patch<<","<<my_level<<", <"<<e[0]<<","<<e[2]<<","<<e[4]<<"> to <"<<e[1]<<","<<e[3]<<","<<e[5]<<">)\n";
      rdb->SetIndicesForAMRPatch(patch, my_level, e);
    }

    rdb->CalculateBoundaries();

    this->mesh_boundaries[meshname]=void_ref_ptr(rdb,avtStructuredDomainBoundaries::Destruct);
    
    //
    // Domain Nesting
    //
    avtStructuredDomainNesting *dn = new avtStructuredDomainNesting(totalPatches, num_levels);
    //debug5<<"Calculating avtStructuredDomainNesting for "<<meshname<<" mesh ("<<dn<<").\n";
    dn->SetNumDimensions(3);

    //
    // Calculate what the refinement ratio is from one level to the next.
    //
    for (int level = 0 ; level < num_levels ; level++) {
      // SetLevelRefinementRatios requires data as a vector<int>
      vector<int> rr(3);
      for (int i=0; i<3; i++)
        rr[i] = stepInfo->levelInfo[level].refinementRatio[i];

      //debug5<<"\tdn->SetLevelRefinementRatios("<<level<<", <"<<rr[0]<<","<<rr[1]<<","<<rr[2]<<">)\n";
      dn->SetLevelRefinementRatios(level, rr);
    }

    //
    // Calculating the child patches really needs some better sorting than
    // what I am doing here.  This is likely to become a bottleneck in extreme
    // cases.  Although this routine has performed well for a previous 55K
    // patch run.
    //
    vector< vector<int> > childPatches(totalPatches);
    for (int level = num_levels-1 ; level > 0 ; level--) {
      int prev_level = level-1;
      LevelInfo &levelInfoParent = stepInfo->levelInfo[prev_level];
      LevelInfo &levelInfoChild = stepInfo->levelInfo[level];

      for (int child=0; child<(int)levelInfoChild.patchInfo.size(); child++) {
        PatchInfo &childPatchInfo = levelInfoChild.patchInfo[child];
        int child_low[3],child_high[3];
        childPatchInfo.getBounds(child_low,child_high,meshname);

        for (int parent=0; parent<(int)levelInfoParent.patchInfo.size(); parent++) {
          PatchInfo &parentPatchInfo = levelInfoParent.patchInfo[parent];
          int parent_low[3],parent_high[3];
          parentPatchInfo.getBounds(parent_low,parent_high,meshname);

          int mins[3], maxs[3];
          for (int i=0; i<3; i++) {
            mins[i] = max(child_low[i],  parent_low[i] *levelInfoChild.refinementRatio[i]);
            maxs[i] = min(child_high[i], parent_high[i]*levelInfoChild.refinementRatio[i]);
          }
          bool overlap = (mins[0]<maxs[0] &&
                          mins[1]<maxs[1] &&
                          mins[2]<maxs[2]);

          if (overlap) {
            int child_gpatch = GetGlobalDomainNumber(level, child);
            int parent_gpatch = GetGlobalDomainNumber(prev_level, parent);
            childPatches[parent_gpatch].push_back(child_gpatch);
          }
        }
      }
    }

    //
    // Now that we know the extents for each patch and what its children are,
    // tell the structured domain boundary that information.
    //
    for (int p=0; p<totalPatches ; p++) {
      int my_level, local_patch;
      GetLevelAndLocalPatchNumber(p, my_level, local_patch);

      PatchInfo &patchInfo = stepInfo->levelInfo[my_level].patchInfo[local_patch];
      int low[3],high[3];
      patchInfo.getBounds(low,high,meshname);

      vector<int> e(6);
      for (int i=0; i<3; i++) {
        e[i+0] = low[i];
        e[i+3] = high[i]-1;
      }
      //debug5<<"\tdn->SetNestingForDomain("<<p<<","<<my_level<<", <>, <"<<e[0]<<","<<e[1]<<","<<e[2]<<"> to <"<<e[3]<<","<<e[4]<<","<<e[5]<<">)\n";
      dn->SetNestingForDomain(p, my_level, childPatches[p], e);
    }

    this->mesh_domains[meshname]=void_ref_ptr(dn, avtStructuredDomainNesting::Destruct);
    forceMeshReload = false;
  }

  //
  // Register these structures with the generic database so that it knows
  // to ghost out the right cells.
  //
  cache->CacheVoidRef("any_mesh", // key MUST be called any_mesh
                      AUXILIARY_DATA_DOMAIN_BOUNDARY_INFORMATION,
                      timestate, -1, this->mesh_boundaries[meshname]);
  cache->CacheVoidRef("any_mesh", // key MUST be called any_mesh
                      AUXILIARY_DATA_DOMAIN_NESTING_INFORMATION,
                      timestate, -1, this->mesh_domains[meshname]);

  //VERIFY we got the mesh boundary and domain in there
  void_ref_ptr vrTmp = cache->GetVoidRef("any_mesh", // MUST be called any_mesh
                            AUXILIARY_DATA_DOMAIN_BOUNDARY_INFORMATION,
                            timestate, -1);
  if (*vrTmp == NULL || *vrTmp != *this->mesh_boundaries[meshname])
    throw InvalidFilesException("uda boundary mesh not registered");

  vrTmp = cache->GetVoidRef("any_mesh", // MUST be called any_mesh
                            AUXILIARY_DATA_DOMAIN_NESTING_INFORMATION,
                            timestate, -1);
  if (*vrTmp == NULL || *vrTmp != *this->mesh_domains[meshname])
    throw InvalidFilesException("uda domain mesh not registered");

  visitTimer->StopTimer(t1, "avtUintahFileFormat::CalculateDomainNesting");
#endif
}


// ***************************************************************************
//  Method: avtUintahFileFormat::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: sshankar -- generated by xml2avt
//  Creation:   Tue May 13 19:02:26 PST 2008
//
// ****************************************************************************

vtkIntArray *ConvertToInt(vtkDataArray *input)
{
  vtkIntArray *iarr = vtkIntArray::New();
  iarr->SetNumberOfComponents(input->GetNumberOfComponents());
  iarr->SetNumberOfTuples(input->GetNumberOfTuples());
  for(vtkIdType i = 0; i < input->GetNumberOfTuples(); ++i)
  {
    //iarr->SetTuple(i, input->GetTuple(i));
    double foo=i;
    iarr->SetTuple(i,&foo);
  }
  return iarr;
}

vtkLongArray *ConvertToLong(vtkDataArray *input)
{
  vtkLongArray *iarr = vtkLongArray::New();
  iarr->SetNumberOfComponents(input->GetNumberOfComponents());
  iarr->SetNumberOfTuples(input->GetNumberOfTuples());
  for(vtkIdType i = 0; i < input->GetNumberOfTuples(); ++i)
  {
    //iarr->SetTuple(i, input->GetTuple(i));
    double foo=i;
    iarr->SetTuple(i,&foo);
  }
  return iarr;
}

vtkDataSet *
avtUintahFileFormat::GetMesh(int timestate, int domain, const char *meshname)
{
  int t1 = visitTimer->StartTimer();
  
  //debug5<<"avtUintahFileFormat::GetMesh(timestate="<<timestate<<",domain="<<domain<<",meshname="<<meshname<<std::endl;

  ActivateTimestep(timestate);

  string meshName(meshname);

  int level, local_patch;
  GetLevelAndLocalPatchNumber(domain, level, local_patch);

  // volume data
  if (meshName.find("Particle_Mesh") == string::npos) {

    // make sure we have ghosting info for this mesh
    CalculateDomainNesting(timestate, meshname);

    LevelInfo &levelInfo = stepInfo->levelInfo[level];

    //get global bounds
    int glow[3], ghigh[3];
    getBounds(glow,ghigh,meshName,levelInfo);

    //get patch bounds
    int low[3], high[3];
    getBounds(low,high,meshName,levelInfo,local_patch);

    vtkRectilinearGrid *rgrid = vtkRectilinearGrid::New();
    //debug5<<"Calculating vtkRectilinearGrid mesh for "<<meshName<<" mesh ("<<rgrid<<").\n";

    if (meshName.find("NC_") != string::npos)
      nodeCentered = true;

    int dims[3];
    for (int i=0; i<3; i++) 
    {
      int offset = 1; // always one for non-node-centered
      if (nodeCentered && high[i] == ghigh[i]) 
        offset = 0; // nodeCentered and patch end is on high boundary

      dims[i] = high[i]-low[i]+offset;
    }

    rgrid->SetDimensions(dims);

    // need these to offset grid points in order to preserve face 
    // centered locations on node-centered domain.
    bool sfck[3] = {meshName.find("SFCX") != string::npos,
                    meshName.find("SFCY") != string::npos,
                    meshName.find("SFCZ") != string::npos};

    // Set the coordinates of the grid points in each direction.
    for (int c=0; c<3; c++) {
      vtkFloatArray *coords = vtkFloatArray::New(); 
      coords->SetNumberOfTuples(dims[c]); 
      float *array = (float *)coords->GetVoidPointer(0); 
      for (int i=0; i<dims[c]; i++)
      {
        // Face centered data gets shifted towards -inf by half a cell.
        // Boundary patches are special shifted to preserve global domain.
        // Internal patches are always just shifted.
        float face_offset=0;
        if (sfck[c]) 
        {
          if (i==0)
            if (low[c]==glow[c]) // patch is on low boundary
              face_offset = 0.0;
            else
              face_offset = -0.5;       // patch boundary is internal to the domain
          else if (i==dims[c]-1)
            if (high[c]==ghigh[c]) // patch is on high boundary
              if (levelInfo.periodic[c])  // periodic means one less value in the face-centered direction
                face_offset = 0.0;
              else
                face_offset = -1;
            else                        // patch boundary is internal to the domain
              face_offset = -0.5;
          else
            face_offset = -0.5;
        }
        array[i] = levelInfo.anchor[c] +
          (i + low[c] + face_offset) * levelInfo.spacing[c];
      }

      switch(c) {
      case 0:
        rgrid->SetXCoordinates(coords); break;
      case 1:
        rgrid->SetYCoordinates(coords); break;
      case 2:
        rgrid->SetZCoordinates(coords); break;
      }

      coords->Delete();
    }

    visitTimer->StopTimer(t1, "avtUintahFileFormat::GetMesh() Volume Grid");

    return rgrid;
  }

  // particle data
  else if (meshName.find("Particle_Mesh") != string::npos) {

    size_t found = meshName.find("/");
    string matl = meshName.substr(found + 1);

    int matlNo = -1;
    if (matl.compare("*") != 0)
      matlNo = atoi(matl.c_str());

    int t2 = visitTimer->StartTimer();

    // we always want p.x when setting up the mesh
//    string vars = "p.x";
    debug5<<__FILE__<<"  "<<__FUNCTION__<<"  "<< __LINE__<< std::endl;
    string vars = (*getParticlePositionName)(archive);
    debug5<<__FILE__<<"  "<<__FUNCTION__<<"  "<< __LINE__<< std::endl;

    debug5<<__FILE__<<"  "<<__FUNCTION__<<"  "<< __LINE__<< std::endl;
    ParticleDataRaw *pd = (*getParticleData)(archive, grid, level, local_patch, vars, matlNo, timestate);
    debug5<<__FILE__<<"  "<<__FUNCTION__<<"  "<< __LINE__<< std::endl;

    visitTimer->StopTimer(t2, "avtUintahFileFormat::GetMesh() getParticleData");

    // Create the vtkPoints object and copy points into it.
    vtkDoubleArray *doubleArray = vtkDoubleArray::New();
    doubleArray->SetNumberOfComponents(3);
    doubleArray->SetArray(pd->data, pd->num*pd->components, 0);

    vtkPoints *points = vtkPoints::New();
    points->SetData(doubleArray);
    doubleArray->Delete();

    // 
    // Create a vtkUnstructuredGrid to contain the point cells. 
    // 
    vtkUnstructuredGrid *ugrid = vtkUnstructuredGrid::New(); 
    ugrid->SetPoints(points); 
    points->Delete(); 
    ugrid->Allocate(pd->num); 
    vtkIdType onevertex; 

    for(int i = 0; i < pd->num; ++i) {
      onevertex = i; 
      ugrid->InsertNextCell(VTK_VERTEX, 1, &onevertex); 
    } 

    // don't delete pd->data - vtk owns it now!
    delete pd;

    //try to retrieve existing cache ref
    void_ref_ptr vrTmp = cache->GetVoidRef(meshname,AUXILIARY_DATA_GLOBAL_NODE_IDS,timestate,domain);
    vtkDataArray *pID=NULL;
    if (*vrTmp == NULL)
    {
      //
      // add globel node ids to facilitate point cloud usage
      //
      //basically same as GetVar(timestate, domain, "particleID");
      int level, local_patch;
      //debug5<<"\tGetLevelAndLocalPatchNumber...\n";
      GetLevelAndLocalPatchNumber(domain, level, local_patch);
      int matlNo = -1;
      if (matl.compare("*") != 0)
        matlNo = atoi(matl.c_str());

      ParticleDataRaw *pd = NULL;

      //debug5<<"\t(*getParticleData)...\n";
      //todo: this returns an array of doubles. Need to return expected datatype to avoid unnecessary conversion.

#if (1 <= UINTAH_MAJOR_VERSION && 7 <= UINTAH_MINOR_VERSION )
      if( variableExists(archive, "p.particleID") )
#endif
      {
        int t2 = visitTimer->StartTimer();

        debug5<<__FILE__<<"  "<<__FUNCTION__<<"  "<< __LINE__<< std::endl;
        pd = (*getParticleData)(archive, grid, level, local_patch, "p.particleID", matlNo, timestate);
        debug5<<__FILE__<<"  "<<__FUNCTION__<<"  "<< __LINE__<< std::endl;

        visitTimer->StopTimer(t2, "avtUintahFileFormat::GetMesh() getParticleData");
      }
      
      //debug5 << "got particle data: "<<pd<<"\n";
      if (pd)
      {
        vtkDoubleArray *rv = vtkDoubleArray::New();
        //vtkLongArray *rv = vtkLongArray::New();
        //debug5<<"\tSetNumberOfComponents("<<pd->components<<")...\n";
        rv->SetNumberOfComponents(pd->components);

        //debug5<<"\tSetArray...\n";
        rv->SetArray(pd->data, pd->num*pd->components, 0);

        // don't delete pd->data - vtk owns it now!
        delete pd;
        
        //todo: this is the unnecesary conversion, from long int->double->int, to say nothing of the implicit curtailing that might occur (note also: this is a VisIt bug that uses ints to store particle ids rather than long ints)
        vtkIntArray *iv=ConvertToInt(rv);
        //vtkLongArray *iv=ConvertToLong(rv);
        rv->Delete(); // this should now delete pd->data

        pID=iv;
      }

      //debug5<<"read particleID ("<<pID<<")\n";
      if(pID != NULL)
      {
        //debug5<<"adding global node ids from particleID\n";
        pID->SetName("avtGlobalNodeId");
        void_ref_ptr vr = void_ref_ptr( pID , avtVariableCache::DestructVTKObject );
        cache->CacheVoidRef( meshname , AUXILIARY_DATA_GLOBAL_NODE_IDS , timestate , domain , vr);

        //make sure it worked
        void_ref_ptr vrTmp = cache->GetVoidRef(meshname,AUXILIARY_DATA_GLOBAL_NODE_IDS,timestate,domain);
        if (*vrTmp == NULL || *vrTmp != *vr)
          throw InvalidFilesException("failed to register uda particle global node");
      }
    }

    visitTimer->StopTimer(t1, "avtUintahFileFormat::GetMesh() Particle Grid");

    return ugrid;
  }

  return NULL;
}


// ****************************************************************************
//  Method: avtUintahFileFormat::GetVar
//
//  Purpose:
//      Gets a scalar variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkDoubleArray, 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: sshankar -- generated by xml2avt
//  Creation:   Tue May 13 19:02:26 PST 2008
//
// ****************************************************************************
vtkDataArray *
avtUintahFileFormat::GetVar(int timestate, int domain, const char *varname)
{
  int t1 = visitTimer->StartTimer();

  ActivateTimestep(timestate);

  string varName(varname);
  bool isParticleVar = false;
    
  size_t found = varName.find("/");
  string tmpVarName = varName;
    
  string matl = varName.substr(found + 1);
  varName = varName.substr(0, found);
    
  string varType="CC_Mesh";
  if (strcmp(varname, "proc_id")!=0) {
    for (int k=0; k<(int)stepInfo->varInfo.size(); k++) {
      if (stepInfo->varInfo[k].name == varName) {
        varType = stepInfo->varInfo[k].type;
        if (stepInfo->varInfo[k].type.find("ParticleVariable") != string::npos) {
          isParticleVar = true;
          break;
        }
      }
    }
  }

  vtkDoubleArray *rv = vtkDoubleArray::New();

  int level, local_patch;
  GetLevelAndLocalPatchNumber(domain, level, local_patch);

  // volume data
  if (!isParticleVar) { 

    LevelInfo &levelInfo = stepInfo->levelInfo[level];
    PatchInfo &patchInfo = levelInfo.patchInfo[local_patch];

    int qlow[3], qhigh[3]; // region we're going to ask uintah for (from qlow to qhigh-1)
    patchInfo.getBounds(qlow,qhigh,varType);

    GridDataRaw *gd=NULL;

    if (strcmp(varname, "proc_id")==0) {
      gd = new GridDataRaw;
      for (int i=0; i<3; i++) {
        gd->low[i ] = qlow[i];
        gd->high[i] = qhigh[i];
      }
      gd->components=1;

      int ncells = (qhigh[0]-qlow[0])*(qhigh[1]-qlow[1])*(qhigh[2]-qlow[2]);
      gd->data = new double[ncells];
      for (int i=0; i<ncells; i++) 
        gd->data[i] = patchInfo.getProcId();
    }

    else {
      if (nodeCentered == true){
        int glow[3], ghigh[3];
        getBounds(glow,ghigh,varType, levelInfo);
        patchInfo.getBounds(qlow,qhigh,varType);
                
        for (int j=0; j<3; j++)
          if (qhigh[j] != ghigh[j]) // patch is on low boundary
            qhigh[j] = qhigh[j]+1;
          else
            qhigh[j] = qhigh[j];
      }

#ifdef SERIALIZED_READS
      int numProcs, rank;
      int msg = 128, tag = 256;
      MPI_Status status;

      MPI_Comm_size(VISIT_MPI_COMM, &numProcs);
      MPI_Comm_rank(VISIT_MPI_COMM, &rank);

      int totalPatches = 0;
      for (int i = 0; i < stepInfo->levelInfo.size(); i++)
        totalPatches += stepInfo->levelInfo[i].patchInfo.size();

      // calculate which process we should wait for a message from
      // if we're processing doiman 0 don't wait for anyone else
      int prev = (rank+numProcs-1)%numProcs;
      int next = (rank+1)%numProcs;

      // domain 0 always reads right away
      if (domain==0)
        prev = -1;
      //debug5 << "Proc: " << rank << " sent to GetVar" << endl;

      // wait for previous read to finish
      if (prev>=0)
        MPI_Recv(&msg, 1, MPI_INT, prev, tag, VISIT_MPI_COMM, &status);

      int t2 = visitTimer->StartTimer();

      debug5<<__FILE__<<"  "<<__FUNCTION__<<"  "<< __LINE__<< std::endl;
      gd = (*getGridData)(archive, grid, level, local_patch, varName, atoi(matl.c_str()), timestate, qlow, qhigh);
      debug5<<__FILE__<<"  "<<__FUNCTION__<<"  "<< __LINE__<< std::endl;

      visitTimer->StopTimer(t2, "avtUintahFileFormat::GetMesh() getGridData");
      // let the next read go
      if (next>=0)
        MPI_Send(&msg, 1, MPI_INT, next, tag, VISIT_MPI_COMM);

#else
      int t2 = visitTimer->StartTimer();
      
      debug5<<__FILE__<<"  "<<__FUNCTION__<<"  "<< __LINE__<< std::endl;
      gd = (*getGridData)(archive, grid, level, local_patch, varName, atoi(matl.c_str()), timestate, qlow, qhigh);
      debug5<<__FILE__<<"  "<<__FUNCTION__<<"  "<< __LINE__<< std::endl;

      visitTimer->StopTimer(t2, "avtUintahFileFormat::GetVar getGridData");
#endif
    }

    int n = (qhigh[0]-qlow[0])*(qhigh[1]-qlow[1])*(qhigh[2]-qlow[2]);

    CheckNaNs(gd->data,n*gd->components,varname,level,local_patch);

    rv->SetNumberOfComponents(gd->components);
    rv->SetArray(gd->data, n*gd->components, 0);

    // don't delete gd->data - vtk owns it now!
    delete gd;
  }

  // particle data
  else {
    int matlNo = -1;
    if (matl.compare("*") != 0)
      matlNo = atoi(matl.c_str());

    ParticleDataRaw *pd = NULL;

#ifdef SERIALIZED_READS
    int numProcs, rank;
    int msg = 128, tag = 256;
    MPI_Status status;

    MPI_Comm_size(VISIT_MPI_COMM, &numProcs);
    MPI_Comm_rank(VISIT_MPI_COMM, &rank);

    int totalPatches = 0;
    for (int i = 0; i < stepInfo->levelInfo.size(); i++)
      totalPatches += stepInfo->levelInfo[i].patchInfo.size();

    // calculate which process we should wait for a message from
    // if we're processing doiman 0 don't wait for anyone else
    int prev = (rank+numProcs-1)%numProcs;
    int next = (rank+1)%numProcs;

    // domain 0 always reads right away
    if (domain==0)
      prev = -1;
    //debug5 << "Proc: " << rank << " sent to GetVar" << endl;

    // wait for previous read to finish
    if (prev>=0)
      MPI_Recv(&msg, 1, MPI_INT, prev, tag, VISIT_MPI_COMM, &status);

    int t2 = visitTimer->StartTimer();

    debug5<<__FILE__<<"  "<<__FUNCTION__<<"  "<< __LINE__<< std::endl;
    pd = (*getParticleData)(archive, grid, level, local_patch, varName, matlNo, timestate);
    debug5<<__FILE__<<"  "<<__FUNCTION__<<"  "<< __LINE__<< std::endl;

    visitTimer->StopTimer(t2, "avtUintahFileFormat::GetVar getParticleData");

    // let the next read go
    if (next>=0)
      MPI_Send(&msg, 1, MPI_INT, next, tag, VISIT_MPI_COMM);

#else
    int t2 = visitTimer->StartTimer();

    debug5<<__FILE__<<"  "<<__FUNCTION__<<"  "<< __LINE__<< std::endl;
    pd = (*getParticleData)(archive, grid, level, local_patch, varName, matlNo, timestate);
    debug5<<__FILE__<<"  "<<__FUNCTION__<<"  "<< __LINE__<< std::endl;

    visitTimer->StopTimer(t2, "avtUintahFileFormat::GetVar getParticleData");
#endif

    CheckNaNs(pd->data,pd->num*pd->components,varname,level,local_patch);

    rv->SetNumberOfComponents(pd->components);
    rv->SetArray(pd->data, pd->num*pd->components, 0);

    // don't delete pd->data - vtk owns it now!
    delete pd;
  }

  visitTimer->StopTimer(t1, "avtUintahFileFormat::GetVar");

  return rv;
}


// ****************************************************************************
//  Method: avtUintahFileFormat::GetVectorVar
//
//  Purpose:
//      Gets a vector variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkDoubleArray, 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: sshankar -- generated by xml2avt
//  Creation:   Tue May 13 19:02:26 PST 2008
//
// ****************************************************************************
vtkDataArray *
avtUintahFileFormat::GetVectorVar(int timestate, int domain,const char *varname)
{
  // we handle vector variables exactly the same way as scalars
  return GetVar(timestate, domain, varname);
}


// ****************************************************************************
//  Method: avtBoxlib2DFileFormat::GetAuxiliaryData
//
//  Purpose:
//      Gets the auxiliary data specified.
//
//  Arguments:
//      var        The variable of interest.
//      dom        The domain of interest.
//      type       The type of auxiliary data.
//      <unnamed>  The arguments for that type -- not used.
//      df         Destructor function.
//
//  Returns:    The auxiliary data.
//
//  Programmer: Hank Childs
//  Creation:   January 22, 2006
//
//  Modifications:
//    Kathleen Bonnell, Mon Aug 14 16:40:30 PDT 2006
//    API change for avtIntervalTree.
//
//    Gunther H. Weber, Tue Aug  7 16:01:28 PDT 2007
//    Return material information
// ****************************************************************************
void *
avtUintahFileFormat::GetAuxiliaryData(const char *var, int dom,
                                             const char * type, void *,
                                             DestructorFunction &df)
{
  return NULL;
}

// ****************************************************************************
//  Method: avtUintahFileFormat::CheckNaNs
//
//  Purpose:
//      Check for and warn about NaN values in the file.
//
//  Arguments:
//      num        data size
//      data       data
//      level      level that contains this patch
//      patch      patch that contains these cells
//
//  Returns:    none
//
//  Programmer: cchriste
//  Creation:   06.02.2012
//
//  Modifications:
void
avtUintahFileFormat::CheckNaNs(double *data, const int num,
                               const char* varname,
                               const int level, const int patch)
{
  // replace nan's with a large negative number
  std::vector<int> nanCells;
  for (int i=0; i<num; i++) 
  {
    if (std::isnan(data[i]))
    {
      data[i] = NAN_REPLACE_VAL;
      nanCells.push_back(i);
    }
  }

  if (!nanCells.empty())
  {
    std::stringstream sstr;
    sstr << "NaNs exist for variable " << varname
         << " in patch " << patch << " of level " << level
         << " and " << nanCells.size() << "/" << num
         << " cells have been replaced by the value "
         <<  NAN_REPLACE_VAL << ".";

    // if ((int)nanCells.size()>40)
    // {
    //   sstr<<"\nFirst 20: ";
    //   for (int i=0;i<(int)nanCells.size() && i<20;i++)
    //     sstr<<nanCells[i]<<",";
    //   sstr<<"\nLast 20: ";
    //   for (int i=(int)nanCells.size()-21;i<(int)nanCells.size();i++)
    //     sstr<<nanCells[i]<<",";
    // }
    // else
    // {
    //   for (int i=0;i<(int)nanCells.size();i++)
    //     sstr<<nanCells[i]<<((int)nanCells.size()!=(i+1)?",":".");
    // }

    avtCallback::IssueWarning(sstr.str().c_str());
  }
}

