/*****************************************************************************
*
* 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 <string>
#include <vector>
#include <dlfcn.h>

#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>

//#ifdef PARALLEL
// 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),
  forceMeshReload(true),
  mesh_for_patch_data(""),
  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 )
  {
    std::string error = std::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 );
  }

  std::string lineStr = line;
  if( !result || lineStr.find( "<Uintah_DataArchive>" ) == std::string::npos )
  {
    std::string error = std::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], bool)) dlsym(libHandle, "getGridData");
  if((error = dlerror()) != NULL)
  {
    EXCEPTION1(InvalidDBTypeException, "The function getGridData could not be located in the library!!!");
  }

#if (2 <= UINTAH_MAJOR_VERSION && 0 <= 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
  std::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");
}


// ****************************************************************************
//  Method: avtUintahFileFormat Destructor
//
//  Programmer: sshankar -- generated by xml2avt
//  Creation:   Tue May 13 19:02:26 PST 2008
//
// ****************************************************************************

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: avtEMSTDFileFormat::AddExpressionsToMetadata
//
//  Purpose:
//      Adds the Uintah specific expressions to the meta data.
//
//  Programmer: sshankar -- generated by xml2avt
//  Creation:   Tue May 13 19:02:26 PST 2008
//
// ****************************************************************************
void
avtUintahFileFormat::AddExpressionsToMetadata(avtDatabaseMetaData *md)
{
  std::string home(getenv("HOME"));
  std::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++;

    std::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;
    }

    std::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;
    }
    
    std::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: avtUintahFileFormat::ReadMetaData
//
//  Purpose:
//      Does the actual work for PopulateMetaData().
//
//  Programmer: sshankar -- generated by xml2avt
//  Creation:   Tue May 13 19:02:26 PST 2008
//
// ****************************************************************************
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();

  std::vector<int> groupIds(totalPatches);
  std::vector<std::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 node data unless NC_Mesh exists (some only have CC_MESH
  // or SFCk_MESH)
  bool addNodeData = false;
  
  // Don't add patch data unless CC_Mesh or NC_Mesh exists (some only
  // have SFCk_MESH)
  bool addPatchData = false;
  mesh_for_patch_data = "";
  
  // 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<std::string> meshes_added;

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

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

  // 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") == std::string::npos)
    {
      std::string varname = stepInfo->varInfo[i].name;
      std::string vartype = stepInfo->varInfo[i].type;

      avtCentering cent;
      std::string mesh_for_this_var;

      if (vartype.find("NC") != std::string::npos)
      {
        mesh_for_this_var.assign("NC_Mesh");
        cent = AVT_NODECENT;
        // Use the NC_Mesh for node based data.
        addNodeData = true;
        // Use the NC Mesh if there is no CC_Mesh
        if( addPatchData != true )
        {
          addPatchData = true;
          mesh_for_patch_data = mesh_for_this_var;
        }
      }  
      else if (vartype.find("CC") != std::string::npos)
      {  
        mesh_for_this_var.assign("CC_Mesh");
        cent = AVT_ZONECENT;
        // Use the CC_Mesh for patch based data.
        addPatchData = true;
        mesh_for_patch_data = mesh_for_this_var;
      }
      else if (vartype.find("SFC") != std::string::npos)
      { 
        if (vartype.find("SFCX") != std::string::npos)               
          mesh_for_this_var.assign("SFCX_Mesh");
        else if (vartype.find("SFCY") != std::string::npos)          
          mesh_for_this_var.assign("SFCY_Mesh");
        else if (vartype.find("SFCZ") != std::string::npos)          
          mesh_for_this_var.assign("SFCZ_Mesh");

        cent = AVT_ZONECENT;
      }  
      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 the mesh vars
      for (int j=0; j<(int)stepInfo->varInfo[i].materials.size(); j++)
      {
        char buffer[128];
        std::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);

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

  // Add the node data
  if (addNodeData)
  {
    for (std::set<std::string>::iterator it=meshes_added.begin();
         it!=meshes_added.end(); ++it)
    {
      avtVectorMetaData *vector = new avtVectorMetaData();
      vector->name = "patch/nodes/" + *it;
      vector->meshName = "NC_Mesh";
      vector->centering = AVT_NODECENT;
      vector->hasDataExtents = false;
      vector->varDim = 3;
      md->Add(vector);
    }
  }
  
  // Add the patch data
  if (addPatchData)
  {
    avtCentering cent = (mesh_for_patch_data == "CC_Mesh" ?
                         AVT_ZONECENT : AVT_NODECENT);

    avtScalarMetaData *scalar;

#if (2 <= UINTAH_MAJOR_VERSION && 1 <= UINTAH_MINOR_VERSION )
    scalar = new avtScalarMetaData();
    scalar->name = "patch/id";
    scalar->meshName = mesh_for_patch_data;
    scalar->centering = cent;
    scalar->hasDataExtents = false;
    scalar->treatAsASCII = false;
    md->Add(scalar);
#endif

    scalar = new avtScalarMetaData();
    scalar->name = "patch/processor";
    scalar->meshName = mesh_for_patch_data;
    scalar->centering = cent;
    scalar->hasDataExtents = false;
    scalar->treatAsASCII = false;
    md->Add(scalar);

    for (std::set<std::string>::iterator it=meshes_added.begin();
         it!=meshes_added.end(); ++it)
    {
      avtVectorMetaData *vector = new avtVectorMetaData();
      vector->name = "patch/bounds/low/" + *it;
      vector->meshName = mesh_for_patch_data;
      vector->centering = cent;
      vector->hasDataExtents = false;
      vector->varDim = 3;
      md->Add(vector);
      
      vector = new avtVectorMetaData();
      vector->name = "patch/bounds/high/" + *it;
      vector->meshName = mesh_for_patch_data;
      vector->centering = cent;
      vector->hasDataExtents = false;
      vector->varDim = 3;
      md->Add(vector);
    }
  }

  // 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") != std::string::npos)
    {
      std::string varname = stepInfo->varInfo[i].name;
      std::string vartype = stepInfo->varInfo[i].type;

      // j=-1 -> all materials (*)
      for (int j=-1; j<(int)stepInfo->varInfo[i].materials.size(); ++j)
      {
        std::string mesh_for_this_var = std::string("Particle_Mesh/");
        std::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") != std::string::npos) ||
              (vartype.find("Point") != std::string::npos))
            AddVectorVarToMetaData(md, newVarname, mesh_for_this_var, cent, 3);
          // 9 -> tensor
          else if (vartype.find("Matrix3") != std::string::npos)
            AddTensorVarToMetaData(md, newVarname, mesh_for_this_var, cent, 9);
          // 7 -> vector
          else if (vartype.find("Stencil7") != std::string::npos)
            AddVectorVarToMetaData(md, newVarname, mesh_for_this_var, cent, 7);
          // 4 -> vector
          else if (vartype.find("Stencil4") != std::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 the mesh in the cache and if it's not there, compute it
  if (*this->mesh_domains[meshname]==NULL || forceMeshReload == true)
  {
    // Calculate some info needed 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);

    rdb->SetNumDomains(totalPatches);
    
    //debug5<<"Calculating avtRectilinearDomainBoundaries for "<<meshname<<" mesh ("<<rdb<<").\n";

    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 plow[3], phigh[3];
      patchInfo.getBounds(plow, phigh, meshname);

      // For node based meshes add one if there is a neighbor.
      if( meshname.find("NC_") == 0 )
      {
        int nlow[3], nhigh[3];
        patchInfo.getBounds(nlow, nhigh, "NEIGHBORS");
        
        for (int i=0; i<3; i++)
          phigh[i] += nhigh[i];
      }

      int extents[6] = { plow[0], phigh[0],
                         plow[1], phigh[1],
                         plow[2], phigh[2] };

      rdb->SetIndicesForAMRPatch(patch, my_level, extents);

      //debug5<<"\trdb->SetIndicesForAMRPatch("<<patch<<","<<my_level<<", <"<<e[0]<<","<<e[2]<<","<<e[4]<<"> to <"<<e[1]<<","<<e[3]<<","<<e[5]<<">)\n";
    }

    rdb->CalculateBoundaries();

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

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

      dn->SetLevelRefinementRatios(level, rr);

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

    // Calculating the child patches really needs some better sorting
    // than what is crrently being done.  This is likely to become a
    // bottleneck in extreme cases.  Although this routine has
    // performed well for a previous 55K patch run.
    std::vector< std::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 node based meshes add one if there is a neighbor.
        if( meshname.find("NC_") == 0 )
        {
          int nlow[3], nhigh[3];
          childPatchInfo.getBounds(nlow, nhigh, "NEIGHBORS");
          
          for (int i=0; i<3; i++)
            child_high[i] += nhigh[i];
        }

        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);

          // For node based meshes add one if there is a neighbor.
          if( meshname.find("NC_") == 0 )
          {
            int nlow[3], nhigh[3];
            parentPatchInfo.getBounds(nlow, nhigh, "NEIGHBORS");
            
            for (int i=0; i<3; i++)
              parent_high[i] += nhigh[i];
          }

          int mins[3], maxs[3];
          for (int i=0; i<3; i++)
          {
            mins[i] = std::max( child_low[i], 
                                parent_low[i] *levelInfoChild.refinementRatio[i]);
            maxs[i] = std::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 the extents for each patch is known and what its
    // children are, pass 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 plow[3], phigh[3];
      patchInfo.getBounds(plow, phigh, meshname);

      // For node based meshes add one if there is a neighbor.
      if( meshname.find("NC_") == 0 )
      {
        int nlow[3], nhigh[3];
        patchInfo.getBounds(nlow, nhigh, "NEIGHBORS");
        
        for (int i=0; i<3; i++)
          phigh[i] += nhigh[i];
      }

      std::vector<int> extents(6);
      for (int i=0; i<3; i++)
      {
        extents[i+0] = plow[i];
        extents[i+3] = phigh[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], extents);
    }

    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::ConvertTo*
//
//  Purpose:
//      Converts a data arry to ints or longs.
//
//  Arguments:
//      Input array
//
//  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;
}

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

  std::string meshName(meshname);

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

  // Particle data
  if (meshName.find("Particle_Mesh") != std::string::npos)
  {
    size_t found = meshName.find("/");
    std::string matl = meshName.substr(found + 1);

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

    int t2 = visitTimer->StartTimer();

    // Get the particle position name typically p.x
    std::string vars = (*getParticlePositionName)(archive);

    ParticleDataRaw *pd = (*getParticleData)(archive, grid, level, local_patch,
                                             vars, matlNo, timestate);

    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 global 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 (2 <= UINTAH_MAJOR_VERSION && 0 <= 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;
  }

  // Volume data
  else
  {
    vtkRectilinearGrid *rgrid = vtkRectilinearGrid::New();

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

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

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

    int dims[3];

    // Get the patch bounds
    int plow[3], phigh[3];
    patchInfo.getBounds(plow, phigh, meshName);

    // For node based meshes add one if there is a neighbor.
    if( meshName.find("NC_") == 0 )
    {
      int nlow[3], nhigh[3];
      patchInfo.getBounds(nlow, nhigh, "NEIGHBORS");
      
      for (int i=0; i<3; i++)
      {
        phigh[i] += nhigh[i];
        dims[i] = phigh[i] - plow[i];
      }
    }
    else
    {      
      // For cell and face meshes always add one.
      for (int i=0; i<3; i++) 
      {
        dims[i] = phigh[i] - plow[i] + 1;
      }
    }
    
    rgrid->SetDimensions(dims);    

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

    int nlow[3], nhigh[3];

    if( sfck[0] || sfck[1] || sfck[2] )
      patchInfo.getBounds(nlow, nhigh, "NEIGHBORS");

    // 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;
        
        if (sfck[c])
        {
          if (i==0)
          {
            // No neighbor, so the patch is on low boundary
            if (nlow[c] == 0)
              face_offset = 0.0;
            // Patch boundary is internal to the domain
            else
              face_offset = -0.5;
          }
          else if (i == dims[c]-1)
          {
            // No neighbor, so the patch is on high boundary
            if (nhigh[c] == 0)
            {
              // Periodic means one less value in the face-centered direction
              if (levelInfo.periodic[c])
                face_offset = 0.0;
              else
                face_offset = -1;
            }
            // Patch boundary is internal to the domain
            else
              face_offset = -0.5;
          }
          else
            face_offset = -0.5;
        }
        else
          face_offset = 0;

        array[i] = levelInfo.anchor[c] +
          (i + plow[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;
  }

  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);

  bool isParticleVar = false;

  // Get the var name sans the material. If a patch or processor
  // variable then the var name will be either "patch" or "processor".
  std::string varName(varname);    
  size_t found = varName.find("/");
  std::string matl = varName.substr(found + 1);
  varName = varName.substr(0, found);
    
  // Get the varType except for patch based data which does not come
  // from the data warehouse but instead from internal structures.
  std::string varType;

  if( strncmp(varname, "patch/nodes", 11) == 0 )
  {
    varType = "NC_Mesh";
  }
  else if (strncmp(varname, "patch/", 6) == 0)
  {
    varType = mesh_for_patch_data;
  }
  else
  {
    for (int k=0; k<(int)stepInfo->varInfo.size(); k++)
    {
      if (stepInfo->varInfo[k].name == varName)
      {
        varType = stepInfo->varInfo[k].type;

        // Check for a particle variable
        if (stepInfo->varInfo[k].type.find("ParticleVariable") != std::string::npos) {
          isParticleVar = true;
          break;
        }
      }
    }
  }

  vtkDoubleArray *rv = vtkDoubleArray::New();

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

  // particle data
  if (isParticleVar)
  {
    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 will always reads right away
    if (domain == 0)
      prev = -1;
    //debug5 << "Proc: " << rank << " sent to GetVar" << endl;

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

    int t2 = visitTimer->StartTimer();

    pd = (*getParticleData)(archive, grid, level, local_patch, varName, matlNo, timestate);

    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();

    pd = (*getParticleData)(archive, grid, level, local_patch, varName, matlNo, timestate);

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

  // Volume data
  else
  {
    LevelInfo &levelInfo = stepInfo->levelInfo[level];
    PatchInfo &patchInfo = levelInfo.patchInfo[local_patch];

    bool nodeCentered;
    
    // The region we're going to ask uintah for (from plow to phigh-1)
    int plow[3], phigh[3];
    patchInfo.getBounds(plow, phigh, varType);
    
    // For node based meshes add one if there is a neighbor.
    if( varType.find("NC") != std::string::npos )
    {
      nodeCentered = true;
      
      int nlow[3], nhigh[3];
      patchInfo.getBounds(nlow, nhigh, "NEIGHBORS");
      
      for (int i=0; i<3; i++)
        phigh[i] += nhigh[i];
    }
    else
    {  
      nodeCentered = false;      
    }
    
    GridDataRaw *gd = NULL;

    // The data for these variables does not come from the data
    // warehouse but instead from internal structures.
    if (strncmp(varname, "patch/", 6) == 0)
    {
      gd = new GridDataRaw;

      for (int i=0; i<3; i++)
      {
        gd->low[i]  =  plow[i] + int(nodeCentered == false);
        gd->high[i] = phigh[i] + int(nodeCentered == false);
      }

      gd->num = ((gd->high[0]-gd->low[0]) *
                 (gd->high[1]-gd->low[1]) *
                 (gd->high[2]-gd->low[2]));

      // The patch id and processor are scalar values while the bounds
      // are vector values.
      if (strcmp(varname, "patch/id") == 0 ||
          strcmp(varname, "patch/processor") == 0 )     
        gd->components = 1;
      else // if( strncmp(varname, "patch/nodes",  11) == 0 ||
           //     strncmp(varname, "patch/bounds/low",  16) == 0 ||
           //     strncmp(varname, "patch/bounds/high", 17) == 0)
        gd->components = 3;

      gd->data = new double[gd->num * gd->components];

      // Patch Id
      if (strcmp(varname, "patch/id") == 0 )
      {
#if (2 <= UINTAH_MAJOR_VERSION && 1 <= UINTAH_MINOR_VERSION )
        double value = patchInfo.getPatchId();
        
        for (int i=0; i<gd->num; i++) 
          gd->data[i] = value;
#endif      
      }
      // Patch processor Id
      else if (strcmp(varname, "patch/processor") == 0 )
      {
        double value = patchInfo.getProcId();
      
        for (int i=0; i<gd->num; i++) 
          gd->data[i] = value;
      }
      // Patch node ids
      else if (strncmp(varname, "patch/nodes", 11) == 0 )
      {
        int cc = 0;
        
        for( int k=gd->low[2]; k<gd->high[2]; ++k )
        {
          for( int j=gd->low[1]; j<gd->high[1]; ++j )
          {
            for( int i=gd->low[0]; i<gd->high[0]; ++i )
            {
              gd->data[cc++] = i;
              gd->data[cc++] = j;
              gd->data[cc++] = k;
            }
          }
        }
      }
      // Patch bounds
      else
      {
        // Get the bounds for this mesh as a variable (not for the grid).
        std::string meshname = std::string(varname);
        found = meshname.find_last_of("/");
        meshname = meshname.substr(found + 1);
        
        patchInfo.getBounds(plow, phigh, meshname);
      
        int *value;

        if (strncmp(varname, "patch/bounds/low", 16) == 0 )
          value = &plow[0];
        else if( strncmp(varname, "patch/bounds/high", 17) == 0)
          value = &phigh[0];

        for (int i=0; i<gd->num; i++)
          for (int c=0; c<3; c++)
            gd->data[i*gd->components+c] = value[c];
      }
    }
    // Patch data from the warehouse
    else
    {
#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();

      gd = (*getGridData)(archive, grid, level, local_patch,
                          varName, atoi(matl.c_str()), timestate,
                          plow, phigh, nodeCentered);

      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();
      
      gd = (*getGridData)(archive, grid, level, local_patch,
                          varName, atoi(matl.c_str()), timestate,
                          plow, phigh, nodeCentered);

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

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

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

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

  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)
{
  // Handle vector variables exactly the same way as scalar variables
  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
//
// ****************************************************************************
void *
avtUintahFileFormat::GetAuxiliaryData(const char *var, int dom,
                                      const char * type, void *,
                                      DestructorFunction &df)
{
  return NULL;
}

// ****************************************************************************
//  Method: avtUintahFileFormat::CheckNaNs
//
//  Purpose:
//      Check for and warns about NaN values in the data.
//
//  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
//
// ****************************************************************************
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());
  }
}
