/*****************************************************************************
*
* 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.
*
*****************************************************************************/

// ************************************************************************* //
//                            avtH5PartFileFormat.C                          //
// ************************************************************************* //

#include <visit-config.h>
#include <avtH5PartFileFormat.h>

// VTK
#include <vtkCellData.h>
#include <vtkCellTypes.h>
#include <vtkUnsignedIntArray.h>
#include <vtkFloatArray.h>
#include <vtkDoubleArray.h>
#include <vtkRectilinearGrid.h>
#include <vtkStructuredGrid.h>
#include <vtkUnstructuredGrid.h>

// VisIt
#include <avtDatabaseMetaData.h>

#include <DBOptionsAttributes.h>
#include <Expression.h>

// FastBit
#ifdef HAVE_LIBFASTBIT
  #include <fastbit-config.h>

  #if FASTBIT_IBIS_INT_VERSION < 2000304
    #error "The H5Part plugin requires FastBit version 2.0.3.4 or newer."
  #endif

  #include <avtDataRangeSelection.h>
  #include <avtHistogramSpecification.h>
  #include <avtIdentifierSelection.h>
#endif

// FastQuery
#ifdef HAVE_LIBFASTQUERY
  #include <fastquery-config.h>

  // #if FASTQUERY_VERSION != "0.8.4.10"
  //   #error "The H5Part plugin requires FastQuery version 0.8.4.10 or newer."
  // #endif

  #include <queryProcessor.h>
#endif

#include <avtCallback.h>
#include <InvalidVariableException.h>
#include <InvalidFilesException.h>
#include <NonCompliantFileException.h>
#include <BadIndexException.h>
#include <DebugStream.h>
#include <TimingsManager.h>

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

// ****************************************************************************
//  Method: avtH5PartFileFormat constructor
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Tue eb 9 13:44:50 PST 2010
//
//   Modifications:
//    Gunther H. Weber, Thu Apr  1 18:17:11 PDT 2010
//    Enable domain decomposition per default.
//
//    Gunther H. Weber, Wed Aug  4 15:34:04 PDT 2010
//    Read and store number of time steps in file for use by GetNTimesteps().
//    Use lazy evaluation for queries. Do not run query in
//    RegisterDataSelections(). Instead, use flag "queryResultsValid" to
//    indicate whether current queryResults reflect current query.
//
//    Jeremy Meredith, Wed Mar 30 12:55:29 EDT 2011
//    Close the file if it's no an H5Part file -- HDF5 can be sensitive about
//    not closing files.
//
//    Gunther H. Weber, Wed Aug 17 16:28:46 PDT 2011
//    Write to debug log file if reader was compiled with FastBit support.
//
// ****************************************************************************

avtH5PartFileFormat::avtH5PartFileFormat(const char *filename,
        DBOptionsAttributes *readOpts) : avtMTSDFileFormat(&filename, 1)
{
    // ibis::gVerbose = 1000;
    
#ifdef HAVE_LIBFASTQUERY
    debug5 << "avtH5PartFileFormat compiled with FastBit support." << std::endl;
#endif
    int t1 = visitTimer->StartTimer();

    // Defaults
    enableDomainDecomposition = true;
    variablePathPrefix = std::string("Step#");

#ifdef HAVE_LIBFASTQUERY
    useFastBitIndex = true;
    fastBitIndexPathPrefix = std::string("__H5PartIndex__");
#endif

    // Check options
    if (readOpts != NULL)
    {
        for (int i = 0; i < readOpts->GetNumberOfOptions(); ++i)
        {
            if (readOpts->GetName(i) == "Enable domain decomposition")
                enableDomainDecomposition =
                  readOpts->GetBool("Enable domain decomposition");
            else if (readOpts->GetName(i) == "Variable path prefix")
                variablePathPrefix =
                  readOpts->GetString("Variable path prefix");
#ifdef HAVE_LIBFASTQUERY
            else if (readOpts->GetName(i) == "Use FastBit index")
                useFastBitIndex = readOpts->GetBool("Use FastBit index");
            else if (readOpts->GetName(i) == "FastBit index path prefix")
                fastBitIndexPathPrefix =
                  readOpts->GetString("FastBit index path prefix");
#endif      
        }
    }

    // only open the file briefly to ensure that it is really an
    // H5Part file
    file = H5PartOpenFile(filename, H5PART_READ);

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

    // This function was once supported by H5Part and will be again in the new
    // release. Until it is widely available, comment it out.
    if (H5PartFileIsValid(file) != H5PART_SUCCESS)
    {
        debug1 << "avtH5PartFileFormat::avtH5PartFileFormat(): "
               << "H5PartFileIsValid check failed." << std::endl;

        H5PartCloseFile(file);
        EXCEPTION1(InvalidFilesException, filename);
    }

    // Check for the first time step - if it does not exit bailout as
    // the file is not a H5Part file.
    activeTimeStep = 0;
    if (!H5PartHasStep(file, activeTimeStep))
    {
        H5PartCloseFile(file);
        EXCEPTION1(InvalidFilesException, "Cannot find first time step 0.");
    }

    // Assume a valid H5Part file from here out.
    
    // Get number of time steps in file
    numTimestepsInFile = H5PartGetNumSteps(file);

    // The following should be equivalent to testing the presence of a
    // group named "Group #0".
    if (numTimestepsInFile <= 0)
    {
        debug1 << "avtH5PartFileFormat::avtH5PartFileFormat(): "
               << "File contains " << numTimestepsInFile
               << " <= 0 time steps." << std::endl;

        H5PartCloseFile(file);
        EXCEPTION1(InvalidFilesException, filename);
    }

    for( int i=0; i<numTimestepsInFile; ++i )
    {
      if (H5PartSetStep(file, i) != H5PART_SUCCESS)
      {
          H5PartCloseFile(file);
          EXCEPTION1(InvalidFilesException, "Cannot activate time step.");
      }
      
      cycles.push_back( i );

      h5part_int64_t numStepAttr = H5PartGetNumStepAttribs( file );

      for( h5part_int64_t i=0; i<numStepAttr; ++i )
      {
        char attrib_name[maxVarNameLen];
        h5part_int64_t len_of_attrib_name;
        h5part_int64_t attrib_type;
        h5part_int64_t attrib_nelem;

        H5PartGetStepAttribInfo( file, i, attrib_name, maxVarNameLen,
                                 &attrib_type, &attrib_nelem );
        
        if( strncmp( attrib_name, "Time", 4 ) == 0 ||
            strncmp( attrib_name, "time", 4 ) == 0 )
        {
          h5part_int64_t status = H5PART_SUCCESS;

          if (attrib_type == H5PART_INT64)
          {
            h5part_int64_t time;
            status = H5PartReadStepAttrib(file, attrib_name, &time);
            times.push_back( time );
          }
          else if (attrib_type == H5PART_INT32)
          {
            h5part_int32_t time;
            status = H5PartReadStepAttrib(file, attrib_name, &time);
            times.push_back( time );
          }
          else if (attrib_type == H5PART_FLOAT64)
          {
            h5part_float64_t time;
            status = H5PartReadStepAttrib(file, attrib_name, &time);
            times.push_back( time );
          }
          else if (attrib_type == H5PART_FLOAT32)
          {
            h5part_float32_t time;
            status = H5PartReadStepAttrib(file, attrib_name, &time);
            times.push_back( time );
          }
          else
          {
            EXCEPTION2(NonCompliantFileException, "H5Part Open",
                       "Unsupported value type for coordinates. "
                       "(Supported tyes are FLOAT64, FLOAT32 and INT64.)");
          }
          // FIXME: Possibly handle other data types for coordinates.

          if (status != H5PART_SUCCESS)
          {
            EXCEPTION2(NonCompliantFileException, "H5Part Open",
                       "Could not read time.");
          }
        }
      }
    }
    
    // Activate first time step
    activeTimeStep = 0;
    if (H5PartSetStep(file, activeTimeStep) != H5PART_SUCCESS)
    {
        H5PartCloseFile(file);
        EXCEPTION1(InvalidFilesException, "Cannot activate time step 0.");
    }
  
    // Iterate over particle var names
    h5part_int64_t nPointVars = H5PartGetNumDatasets(file);
    if (nPointVars < 0)
    {
        H5PartCloseFile(file);
        EXCEPTION2(NonCompliantFileException, "H5Part Constructor",
                "Could not obtain the number of particle variables.");
    }

    // Get information about particle variables
    for (int i=0; i<nPointVars; ++i)
    {
        // Get information. We are interested in name and type
        char varName[maxVarNameLen];
        h5part_int64_t type;
        h5part_int64_t nElem;

        if (H5PartGetDatasetInfo(file, i, varName, maxVarNameLen, &type, &nElem)
                != H5PART_SUCCESS)
        {
            H5PartCloseFile(file);
            EXCEPTION2(NonCompliantFileException, "H5Part Constructor",
                       "Could not obtain particle data set information.");
        }

        // Store information about type in map so that we can access
        // it in GetVar()
        if (particleVarNameToTypeMap.find(varName) != particleVarNameToTypeMap.end())
        {
            H5PartCloseFile(file);
            EXCEPTION2(NonCompliantFileException, "H5Part Constructor",
                       "Particle data contains two variables with the same name.");
        }
        else
        {
            particleVarNameToTypeMap[varName] = type;
        }
    }

    particleNSpatialDims = 2;
    unsigned char cartesian = 0;
    unsigned char cylindrical = 0;
    unsigned char spherical = 0;

    // If there is particle data in the file (check needed since a
    // file could contain only block data).
    if (particleVarNameToTypeMap.size())
    {
        // Determine the coordinate system used and whether the file
        // is 2D or 3D
        for (VarNameToInt64Map_t::const_iterator it = particleVarNameToTypeMap.begin();
             it != particleVarNameToTypeMap.end(); ++it)
        {
            std::string currVarName = it->first; 

            std::transform(currVarName.begin(), currVarName.end(),
                           currVarName.begin(), ::tolower);

            if (currVarName == "x" || currVarName == "y" || currVarName == "z" )
                ++cartesian;

            if (currVarName == "r" || currVarName == "phi" || currVarName == "z")
                ++cylindrical;

            if (currVarName == "r" || currVarName == "phi" || currVarName == "theta" )
                ++spherical;

            if (currVarName == "z" || currVarName == "theta" )
                particleNSpatialDims = 3;

            
        }

        // Spherical coordinates always requires three values. 
        if( spherical == 3 && (cylindrical >= 3 || cartesian >= 3) )
        {
          std::stringstream buf;

          buf << "Ambiguous coordinate system: both "
              << "spherical and cylindrical or cartesian coordinate systems "
              << " were found. Can not disambiguate. ";

          debug1 << "avtH5PartFileFormat::avtH5PartFileFormat(): "
                 << buf.str() << std::endl;

          EXCEPTION2(NonCompliantFileException,
                     "H5Part Constructor", buf.str().c_str() );

          // callback is not working so use the exception above
          avtCallback::IssueWarning( buf.str().c_str() );
        }

        // Cylindrical and cartesian coordinates have two or three values.
        else if( cylindrical >= 2 && cartesian >= 2 )
        {
          std::stringstream buf;

          buf << "Ambiguous coordinate system: both "
              << "cylindrical and cartesian coordinate systems were found. "
              << "Can not disambiguate. ";

          debug1 << "avtH5PartFileFormat::avtH5PartFileFormat(): "
                 << buf.str() << std::endl;

          EXCEPTION2(NonCompliantFileException,
                     "H5Part Constructor", buf.str().c_str());

          // callback is not working so use the exception above
          avtCallback::IssueWarning( buf.str().c_str() );
        }

        if( spherical == 3)           { coordType = sphericalCoordSystem; }
        else if( cylindrical >= 2)    { coordType = cylindricalCoordSystem; }
        else /* if( cartesian >= 2)*/ { coordType = cartesianCoordSystem; }

        debug1 << "File contains a " << particleNSpatialDims
               << "D mesh using a "
               << (coordType == cartesianCoordSystem ? "Cartesian" :
                   (coordType == cylindricalCoordSystem ? "Cylindrical" :
                    (coordType == sphericalCoordSystem ? "Spherical" : "Unknown") ) )
               << " coordinate system." << std::endl;
    }

    // Get information about field variables

    // FIXME: Currently we assume that all fields live on the same
    // mesh. This may not be true.
    h5part_int64_t nFieldVars = H5BlockGetNumFields(file);
    if (nFieldVars < 0)
    {
        H5PartCloseFile(file);
        EXCEPTION2(NonCompliantFileException, "H5Part Constructor",
                   "Could not read number of field variables.");
    }
    else if (nFieldVars > 0)
    {
        for (h5part_int64_t idx=0; idx < nFieldVars; ++idx)
        {
            char varName[maxVarNameLen];
            h5part_int64_t gridRank;
            h5part_int64_t gridDims[3];
            h5part_int64_t fieldDims;
            h5part_int64_t type;

            if (H5BlockGetFieldInfo (file, idx, varName, maxVarNameLen,
                                     &gridRank, gridDims, &fieldDims, &type) !=
                H5PART_SUCCESS ) 
            {
                H5PartCloseFile(file);
                EXCEPTION2(NonCompliantFileException, "H5Part Constructor",
                           "Could not read field information.");
            }

            if (fieldDims == 1)
            {
                if (fieldScalarVarNameToTypeMap.find(varName) !=
                    fieldScalarVarNameToTypeMap.end())
                {
                    H5PartCloseFile(file);
                    EXCEPTION2(NonCompliantFileException, "H5Part Constructor",
                               "Field data contains two scalar variables with "
                               "the same name.");
                }
                else 
                {
                    fieldScalarVarNameToTypeMap[varName] = type;
                }
            }
            else
            {
                if (fieldVectorVarNameToTypeMap.find(varName) !=
                    fieldVectorVarNameToTypeMap.end())
                {
                    H5PartCloseFile(file);
                    EXCEPTION2(NonCompliantFileException, "H5Part Constructor",
                               "Field data contains two vector variables with "
                               "the same name.");
                }
                else 
                {
                    fieldVectorVarNameToTypeMap[varName] = type;
                    fieldVectorVarNameToFieldRankMap[varName] = fieldDims;
                }
            }
        }
    }

    // FIXME: This information should be read as attribute
    // from file.
    defaultIdVariableName = "id";

    // This key is really only used for FastBit.
    defaultSortedVariableName = "";
    
    h5part_int64_t numFileAttr = H5PartGetNumFileAttribs( file );
    
    for( h5part_int64_t i=0; i<numFileAttr; ++i )
    {
        char attrib_name[maxVarNameLen];
        h5part_int64_t len_of_attrib_name;
        h5part_int64_t attrib_type;
        h5part_int64_t attrib_nelem;

        H5PartGetFileAttribInfo ( file, i, attrib_name, maxVarNameLen,
                                  &attrib_type, &attrib_nelem );

        if( strncmp( attrib_name, "sortedKey", 9 ) == 0 )
        {
            // If the data was sorted use the sortedkey as the
            // default. Otherwise assume the id var name is "id".
            char sortedKey[maxVarNameLen];
            memset(sortedKey, 0, sizeof(char) * maxVarNameLen);
            
            if (H5PartReadFileAttrib(file, "sortedKey", sortedKey) == H5PART_SUCCESS)
            {
                if( std::string(sortedKey) != std::string("unsorted") )
                {
                    defaultSortedVariableName = sortedKey;
                    // FIXME: This information should be read as
                    // attribute from file.
                    defaultIdVariableName = sortedKey;
                }
            }
        }
    }

    idVariableName = defaultIdVariableName;
    
    // FIXME: Still need to check whether there are duplicate variable
    // names in field, scalar and vector variables and possibly modify
    // the name so that VisIt can distinguish between them.
    H5PartCloseFile(file);
    file = 0; // Mark file as closed

#ifdef HAVE_LIBFASTQUERY
    querySpecified = noQuery;
    queryResultsValid = false;
    dataSelectionActive = false;
    queryString = "";
   
    histoCache.clearCache();
    
    // Flag to see if any FastBit indexing was found.
    int foundFastBitIndexing = 0;

    if (useFastBitIndex)
    {
      FastQuery::FQ fq =
        FastQuery::FQ( filename, FastQuery::FQ_H5Part, filename );

      fq.setVariablePathPrefix( getVariablePathPrefix(activeTimeStep).c_str() );
      fq.setIndexPathPrefix( getFastBitIndexPathPrefix(activeTimeStep).c_str() );

      // Get FastBit information about particle variables
      for ( VarNameToInt64Map_t::const_iterator
             it = particleVarNameToTypeMap.begin();
           it != particleVarNameToTypeMap.end(); ++it)
      // for( const auto &it : particleVarNameToTypeMap )
      {
          // std::string varName = it.first;
          std::string varName = it->first;

          // NOTE: Assumption is that if the first time step has
          // FastBit data that all time steps have that data.
          if( fq.checkForIndex( varName ) == true )
          {
              particleVarNameToFastBitMap[varName.c_str()] = true;
              ++foundFastBitIndexing;
          }
          else
          {
              particleVarNameToFastBitMap[varName.c_str()] = false;
          }
      }

      // If no variables are using fastbit then it can be turned off.
      if( foundFastBitIndexing != particleVarNameToTypeMap.size() )
      {
          useFastBitIndex = false;

          std::stringstream buf;

          buf << "Not all H5Part variables contain FastBit indexing. "
              << "Turn the 'File open option' requesting FastBit off "
              << "and re-open the file.";

          debug1 << "avtH5PartFileFormat::avtH5PartFileFormat(): "
                 << buf.str() << std::endl;

          if( foundFastBitIndexing )
            EXCEPTION2(NonCompliantFileException,
                       "H5Part Constructor", buf.str().c_str() );

          // Callback is not working so use the exception above
          avtCallback::IssueWarning( buf.str().c_str() );
      }
    }
    else
    {
      for (VarNameToInt64Map_t::const_iterator
             it = particleVarNameToTypeMap.begin();
           it != particleVarNameToTypeMap.end(); ++it)
        {
          std::string varName = it->first;
          particleVarNameToFastBitMap[varName.c_str()] = false;
        }
    }
#endif
    
    visitTimer->StopTimer(t1, "H5PartFileFormat::avtH5PartFileFormat()");
}

// ****************************************************************************
//  Method: avtH5PartFileFormat destructor
//
//  Programmer: ghweber
//  Creation:   Tue eb 9 13:44:50 PST 2010
//
//  Modifications:
//    Gunther H. Weber, Wed Aug 17 16:27:07 PDT 2011
//    Close HDF5 file if this has not happened so far.
//
// ****************************************************************************

avtH5PartFileFormat::~avtH5PartFileFormat()
{
    FreeUpResources();
}

// ****************************************************************************
//  Method: avtH5PartFileFormat::GetAuxiliaryData
//
//  Purpose:
//      Gets the auxiliary data, e.g., histograms, from a H5Part file.
//
//  Arguments:
//      var        The variable of interest.
//      type       The type of auxiliary data.
//      <unnamed>  The arguments for that -- not used for any Silo types.
//
//  Returns:    The auxiliary data.  Throws an exception if this is not a
//              supported data type.
//
//  Programmer: Prabhat
//  Creation:   February 26, 2008
//
//  Modifications:
//
//    Hank Childs, Thu Mar  6 08:59:57 PST 2008
//    Add skeleton for creating an avtIdentifierSelection.
//
//    Hank Childs, Thu Jan  7 16:46:10 PST 2010
//    Handle the case where an expression is sent in more gracefully.
// 
//    Gunther H. Weber, Mon Feb 15 20:14:18 PST 2010
//    Merged into new H5Part database plugin.
//
//    Gunther H. Weber, Wed Aug 17 13:10:57 PDT 2011
//    Honor useFastBitIndex option.
//
// ****************************************************************************
#ifdef HAVE_LIBFASTQUERY

void*
avtH5PartFileFormat::GetAuxiliaryData(const char *varName, int ts,
                                      const char *type, void *s,
                                      DestructorFunction &df)
{
    int t1 = visitTimer->StartTimer();

    if (useFastBitIndex && particleVarNameToFastBitMap[varName])
    {
        if (strcmp(type, AUXILIARY_DATA_HISTOGRAM) == 0)
        {
            debug5 << "H5Part trying to get histogram!" << std::endl;

            avtHistogramSpecification *spec = (avtHistogramSpecification *) s;

            TRY
            {
                ConstructHistogram(spec);
            }
            CATCHALL
            {
                debug1 << "Exception thrown during histogram construction."
                       << std::endl
                       << "This is normal when an expression is passed in."
                       << std::endl;
            }
            ENDTRY

            // Don't return the spec ... it was passed in as an input,
            // that was then populated.
            visitTimer->StopTimer(t1, "H5PartFileFormat::GetAuxiliaryData() [Histogram]");
            return NULL;
        }
        else if (strcmp(type, AUXILIARY_DATA_IDENTIFIERS) == 0)
        {
            debug5 << "H5Part trying to get auxiliary data for identifiers "
                   << std::endl;

            std::vector<avtDataSelection *> *ds =
              (std::vector<avtDataSelection *> *) s;

            std::vector<avtDataSelection *> drs;

            for (size_t i = 0; i < ds->size(); ++i)
            {
                if ((strcmp((*ds)[i]->GetType(), "Data Range Selection") == 0) ||
                    ((strcmp((*ds)[i]->GetType(), "Identifier Data Selection") == 0)))
                    drs.push_back((*ds)[i]);
                else
                {
                    visitTimer->StopTimer(t1,
                            "H5PartFileFormat::GetAuxiliaryData() [Data identifiers]");
                    return NULL;
                }
            }

            avtIdentifierSelection *ids =
              ConstructIdentifiersFromDataRangeSelection(drs);
            
            df = avtIdentifierSelection::Destruct;

            visitTimer->StopTimer(t1,
                    "H5PartFileFormat::GetAuxiliaryData() [Data identifiers]");

            return static_cast<void*>(ids);
        }
    }

    visitTimer->StopTimer(t1,
            "H5PartFileFormat::GetAuxiliaryData() [Unimplemented request]");

    return NULL;
}
#endif

// ****************************************************************************
//  Method: avtEMSTDFileFormat::GetNTimesteps
//
//  Purpose:
//      Tells the rest of the code how many timesteps there are in this file.
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Tue Feb 9 13:44:50 PST 2010
//
// Modifications:
//    Gunther H. Weber, Wed Aug  4 15:22:23 PDT 2010
//    Originally this function opened the file and used H5PartGetNumSteps()
//    to obtain the number of time steps. However, this call takes a very
//    long time and was slowing down VisIt considerably. To solve this
//    problem, we read the number of time steps when opening the file and
//    cache it as member variable.
//
// ****************************************************************************

int
avtH5PartFileFormat::GetNTimesteps(void)
{
    return numTimestepsInFile;;
}


// ****************************************************************************
//  Method: avtH5PartFileFormat::GetCycles
//
//  Purpose:
//      Returns the cycles
//
//  Arguments:
//      c          the cycles
//
//  Programmer: allen
//  Creation:   
//
// ****************************************************************************


void avtH5PartFileFormat::GetCycles(std::vector<int> &c)
{
    c = cycles;
}


// ****************************************************************************
//  Method: avtH5PartFileFormat::GetTimes
//
//  Purpose:
//      Returns the times
//
//  Arguments:
//      t          the times
//
//  Programmer: allen
//  Creation:   
//
// ****************************************************************************

void avtH5PartFileFormat::GetTimes(std::vector<double> &t)
{
    t = times;
}


// ****************************************************************************
//  Method: avtH5PartFileFormat::FreeUpResources
//
//  Purpose:
//      When VisIt is done focusing on a particular timestep, it asks that
//      timestep to free up any resources (memory, file descriptors) that
//      it has associated with it. This method is the mechanism for doing
//      that.
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Tue Feb 9 13:44:50 PST 2010
//
// Modifications:
//    Gunther H. Weber, Wed Aug 17 13:14:27 PDT 2011
//    Honor useFastBitIndex option.
//
// ****************************************************************************

void
avtH5PartFileFormat::FreeUpResources(void)
{
    int t1 = visitTimer->StartTimer();

    if (file)
    {
        H5PartCloseFile(file);
        file = 0;
    }

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


// ****************************************************************************
//  Method: avtH5PartFileFormat::DoubleToString
//
//  Purpose:
//      Create string representation of a double variable.
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Tue Feb 9 13:44:50 PST 2010
//
// ****************************************************************************


std::string inline avtH5PartFileFormat::DoubleToString(double x)
{
    std::ostringstream o;
    if (!(o << setprecision(32) << x))
        EXCEPTION1(VisItException, "Error converting double to string.");
    return o.str();
}


// ****************************************************************************
//  Method: avtH5PartFileFormat::RegisterDataSelections
//
//  Purpose:
//      This is the mechanism that allows the database to communicate to this
//      format what data selections are requested. This particular database
//      plugin is interested in data range selections, which are produced
//      by the Threshold operator or identifier data selections, which are
//      produced by named selections.
//
//  Arguments:
//      sels          the data selections
//      selsApplied   a Boolean list stating whether the file format was able
//                    to apply the selections.
//
//  Programmer: Prabhat
//  Creation:   December 19, 2007
//
//  Modifications:
//    Gunther H. Weber, Thu Feb 11 17:02:13 PST 2010
//    Ported/moved from HDF_UC plugin to H5Part plugin
//
//    Gunther H. Weber, Wed Aug  4 16:28:24 PDT 2010
//    Changed to lazy evaluation of query. Just set queryResultsValid to
//    false and perform query when data is read for the first time.
//
//    Gunther H. Weber, Wed Aug 17 13:15:46 PDT 2011
//    Make queries inclusive to match Threshold operator
//
// ****************************************************************************
#ifdef HAVE_LIBFASTQUERY

void
avtH5PartFileFormat::RegisterDataSelections(
        const std::vector<avtDataSelection_p> &selList,
        std::vector<bool> *selsApplied)
{
    int t1 = visitTimer->StartTimer();

    //time_t startTime = time(0);
    debug5 << "avtH5PartFileFormat::RegisterDataSelection() "
           << "useFastBitIndex is " << useFastBitIndex << std::endl;

    dataSelectionActive = false;
    querySpecified = noQuery;
    queryResultsValid = false;
    queryString.clear();
    queryIdList.clear();
    queryResults.clear();

    if (useFastBitIndex)
    {
        idVariableName = defaultIdVariableName;

        if (selList.size() == 1 &&
            std::string(selList[0]->GetType()) ==
            std::string("Identifier Data Selection"))
        {
            // Shortcut: There is only one query and it is a list of
            // identifiers.  Thus, we can skip the step of creating a
            // string and having FastBit parse it again.
            avtIdentifierSelection *ids =
              (avtIdentifierSelection *) *(selList[0]);

            const std::vector<double> &identifiers =
              ids->GetIdentifiers();

            if(!ids->GetIdVariable().empty() &&
               ids->GetIdVariable() != std::string("avtOriginalCellNumbers"))
              idVariableName = ids->GetIdVariable();
            
            bool available = 
              ( (particleVarNameToFastBitMap.find(idVariableName) !=
                 particleVarNameToFastBitMap.end()) &&
                particleVarNameToFastBitMap[idVariableName] );
            
            (*selsApplied)[0] = available;

            debug5 << "avtH5PartFileFormat::RegisterDataSelection() "
                   << "a single Identifier Data Selection was found with "
                   << identifiers.size() << " ids specified "
                   << (available ? "and can" : "but *CANNOT*")
                   << " limit read data based on the variable \""
                   << idVariableName << ".\"" << std::endl;

            if (available && identifiers.size())
            {
                querySpecified = idListQuery;
                queryIdList.resize(identifiers.size());
                
                for (size_t i = 0 ; i < identifiers.size() ; i++)
                {
                    queryIdList[i] = identifiers[i];
                    
                    debug5 << "identifies[ " <<i << "] = " << identifiers[i]
                           << std::endl;
                }
            }
        }
        else
        {
            for (size_t i = 0; i < selList.size(); ++i)
            {
                if (std::string(selList[i]->GetType()) ==
                    std::string("Data Range Selection"))
                {
                    // Check if the variable name is available
                    avtDataRangeSelection *dr =
                      (avtDataRangeSelection *) *(selList[i]);
                    
                    std::string varName = dr->GetVariable();    
 
                    bool available =
                      ( (particleVarNameToFastBitMap.find(varName) !=
                         particleVarNameToFastBitMap.end()) &&
                        particleVarNameToFastBitMap[varName] );

                    (*selsApplied)[i] = available;

                    debug5 << "avtH5PartFileFormat::RegisterDataSelection() "
                           << "for selList[ " << i << "] "
                           << "Data Range Selection was found "
                           << (available ? "and can" : "but *CANNOT*")
                           << " limit read data based on variable \""
                           << varName << ".\"" << std::endl;

                    // If the variable is available then update the
                    // query string accordingly
                    if (available)
                    {
                        double min = dr->GetMin();
                        double max = dr->GetMax();

                        std::string min_cond =
                          varName + ">=" + DoubleToString(min);
                          // DoubleToString(min) + "<=" + varName;
                        std::string max_cond =
                          varName + "<=" + DoubleToString(max);

                        if (queryString.empty())
                            queryString = min_cond + " && " + max_cond;
                        else
                          queryString += " && " + min_cond + " && " + max_cond;

                        querySpecified = stringQuery;
                    }
                }
                else if (std::string(selList[i]->GetType()) ==
                         std::string("Identifier Data Selection"))
                {
                    avtIdentifierSelection *ids =
                        (avtIdentifierSelection *) *(selList[i]);

                    const std::vector<double> &identifiers =
                      ids->GetIdentifiers();

                    if(!ids->GetIdVariable().empty() &&
                       ids->GetIdVariable() != std::string("avtOriginalCellNumbers"))
                      idVariableName = ids->GetIdVariable();
                        
                    bool available =
                      ( (particleVarNameToFastBitMap.find(idVariableName) !=
                         particleVarNameToFastBitMap.end()) &&
                        particleVarNameToFastBitMap[idVariableName] );
                    
                    (*selsApplied)[0] = available;

                    debug5 << "avtH5PartFileFormat::RegisterDataSelection() "
                           << "for selList[ " << i << "] "
                           << "an Identifier Data Selection was found with"
                           << identifiers.size() << " ids specified "
                           << (available ? "and can" : "but *CANNOT*")
                           << " limit read data based on the variable \""
                           << idVariableName << ".\"" << std::endl;

                    if (available && identifiers.size())
                    {
                        // FIXME: Constructing a string is very
                        // inefficient for large selection. However,
                        // there seems no HDF5_FastQuery function to
                        // combine a query for a passed array of
                        // identifiers with a string specifying
                        // thresholds. It may make sense to add this
                        // functionality or add a non-string based API
                        // for range queries and logical combination
                        // of queries.
                        std::string id_string;
                        ConstructIdQueryString(identifiers,
                                               idVariableName, id_string);

                        if (queryString.empty())
                            queryString = id_string;
                        else
                            queryString += " && " + id_string;

                        querySpecified = stringQuery;
                    }
                }
                else
                {
                    (*selsApplied)[i] = false;
                }
            } // over all selList
        }
    }

    visitTimer->StopTimer(t1, "H5PartFileFormat::RegisterDataSelections()");
}
#endif


// ****************************************************************************
//  Method: avtH5PartFileFormat::PopulateDatabaseMetaData
//
//  Purpose:
//      This database meta-data object is like a table of contents for the
//      file.  By populating it, you are telling the rest of VisIt what
//      information it can request from you.
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Tue Feb 9 13:44:50 PST 2010
//
// ****************************************************************************

void
avtH5PartFileFormat::PopulateDatabaseMetaData(avtDatabaseMetaData *md,
                                              int timeState)
{
    int t1 = visitTimer->StartTimer();

    if (particleVarNameToTypeMap.size())
    {
        // AddMeshToMetaData(md, meshname, mt, extents, nblocks, block_origin,
        //                   spatial_dimension, topological_dimension);
        std::string particleMeshName = "particles";
        AddMeshToMetaData(md, particleMeshName, AVT_POINT_MESH,
                          0 /* extents not known */,
                          1 /* single domain */,
                          0, particleNSpatialDims, 1);
        
        for (VarNameToInt64Map_t::const_iterator
               it = particleVarNameToTypeMap.begin();
             it != particleVarNameToTypeMap.end(); ++it)
        {
            // AddScalarVarToMetaData(md, varName, mesh_for_this_var, cent);
            AddScalarVarToMetaData(md, it->first, particleMeshName, AVT_NODECENT);
        }
    }

    // FIXME: Currently we assume that all fields live on the same
    // mesh. This may not be true.
    if (fieldScalarVarNameToTypeMap.size() ||
        fieldVectorVarNameToTypeMap.size())
    {
        std::string fieldMeshName = "fields";

        // AddMeshToMetaData(md, meshname, mt, extents, nblocks, block_origin,
        //                   spatial_dimension, topological_dimension);
        AddMeshToMetaData(md, fieldMeshName, AVT_RECTILINEAR_MESH,
                          0, 1, 0, 3, 3);
        // FIXME: Assume three dimensions for grid data
        // FIXME: Need to determine dimensions of mesh and create a
        // new mesh if dimensions for variables differ

        for (VarNameToInt64Map_t::const_iterator
               it = fieldScalarVarNameToTypeMap.begin();
             it != fieldScalarVarNameToTypeMap.end(); ++it)
        {
            // AddScalarVarToMetaData(md, varName, mesh_for_this_var, cent);
            AddScalarVarToMetaData(md, it->first, fieldMeshName, AVT_NODECENT);
        }

        for (VarNameToInt64Map_t::const_iterator
               it = fieldVectorVarNameToFieldRankMap.begin();
             it != fieldVectorVarNameToFieldRankMap.end(); ++it)
        {
            // AddVectorVarToMetaData(md, varName, mesh_for_this_var, cent,vector_dim);
            AddVectorVarToMetaData(md, it->first, fieldMeshName, AVT_NODECENT, it->second);
        }
    }

    if (enableDomainDecomposition)
        md->SetFormatCanDoDomainDecomposition(true);

    md->SetCyclesAreAccurate(true);
    md->SetCycles( cycles );

    if( times.size() )
    {
      md->SetTimesAreAccurate(true);
      md->SetTimes( times );
    }
    
    //
    // Here's the way to add expressions:
    //Expression momentum_expr;
    //momentum_expr.SetName("momentum");
    //momentum_expr.SetDefinition("{u, v}");
    //momentum_expr.SetType(Expression::VectorMeshVar);
    //md->AddExpression(&momentum_expr);
    //Expression KineticEnergy_expr;
    //KineticEnergy_expr.SetName("KineticEnergy");
    //KineticEnergy_expr.SetDefinition("0.5*(momentum*momentum)/(rho*rho)");
    //KineticEnergy_expr.SetType(Expression::ScalarMeshVar);
    //md->AddExpression(&KineticEnergy_expr);
    //
    visitTimer->StopTimer(t1, "H5PartFileFormat::PopulateDatabaseMetaData()");
}

// ****************************************************************************
//  Method: avtH5PartFileFormat::SelectParticlesToRead
//
//  Purpose:
//      Restrict the "view" in H5Part to those particles that we are
//      interested in: serial, no query -> all particles; serial, query ->
//      query result particles; parallel -> a subset of particles based on
//      query and domain decomposition settings.
//
//  Arguments:
//
//  Programmer: ghweber
//  Creation:   Tue Feb 9 13:44:50 PST 2010
//
//  Modifications:
//    Gunther H. Weber, Thu Apr  1 18:09:09 PDT 2010
//    Fixed computation of index range. Indices passed to H5PartSetView are
//    inclusive!
//
//    Gunther H. Weber, Wed Aug  4 19:18:23 PDT 2010
//    Move re-running of query to this method since it seems more logical
//    here.
//
//    Iuri Prilepov and Gunther H. Weber, Mon Aug 15 15:19:21 PDT 2011
//    Fix parallel decomposition of selections bug.
//
// ****************************************************************************

void
avtH5PartFileFormat::SelectParticlesToRead( const char *varName )
{
    int t1 = visitTimer->StartTimer();

    // Reset the view so all particles are considered.
    H5PartSetView(file, -1, -1);

#ifdef HAVE_LIBFASTQUERY
    // Using FastBit with a query, re-run query if necessary.
    if (useFastBitIndex && particleVarNameToFastBitMap[ varName ] &&
        querySpecified && !queryResultsValid)
    {
        debug5 << "SelectParticlesToRead(): "
               << "Query results invalid: Re-running query for " << varName
               << std::endl;

        // Re-run query
        PerformQuery();
    }
    else if (useFastBitIndex && particleVarNameToFastBitMap[ varName ] &&
             dataSelectionActive)

    {
        debug5 << "SelectParticlesToRead(): "
               << "Query results valid, have " << queryResults.size()
               << " particles" << std::endl;
    }
    
#endif // HAVE_LIBFASTQUERY

#ifdef PARALLEL
    h5part_int64_t nRanks = PAR_Size();
    h5part_int64_t rank   = PAR_Rank();

#ifdef HAVE_LIBFASTQUERY
    if (useFastBitIndex &&  particleVarNameToFastBitMap[ varName ] &&
        dataSelectionActive)
    {
        debug5 << "SelectParticlesToRead(): "
               << "FastBit enabled with a selection active" << std::endl;

        // There is a query/selection. Read only those particles.
        h5part_int64_t nParticles = queryResults.size();

        // Domain decomposition so read only the particles for this rank.
        if (enableDomainDecomposition)
        {
            // Determine what portion of the query will be read by each rank.
            h5part_int64_t myStart, myEnd, totalRead;
            GetDecomp( nParticles, myStart, myEnd, totalRead );
            
            // Restrict the view to the particles read on this rank.
            if (totalRead)
            {
                h5part_int64_t *indices = new h5part_int64_t[totalRead];

                for (h5part_int64_t i=0, j=myStart; j<=myEnd; ++i, ++j)
                    indices[i] = queryResults[j];
                
                H5PartSetViewIndices(file, indices, totalRead);
                delete[] indices;
                
                debug5 << "SelectParticlesToRead(): Rank " << rank
                       << " reads " << totalRead << " of " << nParticles
                       << " particles from " << myStart << " to " << myEnd
                       << " in the query. " << std::endl;
            }
            // This rank does not have any particles to read.
            else
            {
                h5part_int64_t *indices = new h5part_int64_t[0];
                H5PartSetViewIndices(file, indices, 0);
                delete[] indices;

                debug5 << "SelectParticlesToRead(): Rank " << rank
                       << " no particles to read." << std::endl;
            }
        }

        // No domain decomposition so read all the particles rank zero.
        else
        {
            debug5 << "SelectParticlesToRead(): "
                   << "Domain decomposition disabled." << std::endl;
            
            if( PAR_Rank() )
            {
                h5part_int64_t *indices = new h5part_int64_t[0];
                H5PartSetViewIndices(file, indices, 0);
                delete[] indices;
                
                debug5 << "SelectParticlesToRead(): Rank " << rank
                       << " no particles to read." << std::endl;
            }
            else
            {
                // All particles selected are placed on this rank.
                h5part_int64_t *indices = new h5part_int64_t[nParticles];
                std::copy(queryResults.begin(), queryResults.end(), indices);
                H5PartSetViewIndices(file, indices, nParticles);
                delete[] indices;

                debug5 << "SelectParticlesToRead(): Rank " << rank
                       << " All " << nParticles << " particles "
                       << " in the query read." << std::endl;
            }
        }

    }
    else
    {
#endif // HAVE_LIBFASTQUERY

        debug5 << "SelectParticlesToRead(): Rank " << rank
               << " FastBit disabled or no selection active" << std::endl;
        
        // Not using FastBit or no selection. Read the whole file.
        h5part_int64_t nParticles = H5PartGetNumParticles(file);
        
        // Domain decomposition so read only the particles for this rank.
        if (enableDomainDecomposition)
        {
            // Determine what portion of the file will be read by each rank.
            h5part_int64_t idStart, idEnd, totalRead;
            GetDecomp( nParticles, idStart, idEnd, totalRead );

            // Restrict the view to the particles read on this rank.
            if( totalRead )
            {
                H5PartSetView(file, idStart, idEnd);
              
                debug5 << "SelectParticlesToRead(): Rank " << rank
                       << " reads " << totalRead << " of " << nParticles
                       << " particles from " << idStart << " to " << idEnd
                       << " in the file. " << std::endl;
            }
            // This rank does not have any particles to read.
            else
            {
                h5part_int64_t *indices = new h5part_int64_t[0];
                H5PartSetViewIndices(file, indices, 0);
                delete[] indices;
                
                debug5 << "SelectParticlesToRead(): Rank " << rank
                       << " no particles to read." << std::endl;
            }            
        }
        // Read the entire data set on rank zero
        else
        {
            debug5 << "SelectParticlesToRead(): Rank " << rank
                   << " domain decomposition disabled." << std::endl;
            
            if( PAR_Rank() )
            {
                h5part_int64_t *indices = new h5part_int64_t[0];
                H5PartSetViewIndices(file, indices, 0);
                delete[] indices;

                debug5 << "SelectParticlesToRead(): Rank " << rank
                       << " no particles to read." << std::endl;
            }
            else
            {
                H5PartSetView(file, -1, -1);
          
                debug5 << "SelectParticlesToRead(): Rank " << rank
                       << " All " << nParticles << " particles "
                       << " in the file read." << std::endl;
            }
        }

#ifdef HAVE_LIBFASTQUERY
    }
#endif // HAVE_LIBFASTQUERY

#else  // PARALLEL

    // Serial read
#ifdef HAVE_LIBFASTQUERY
    if (useFastBitIndex && particleVarNameToFastBitMap[ varName ] &&
        dataSelectionActive)
    {
        debug5 << "SelectParticlesToRead(): Serial "
               << "FastBit enabled with a selection active" << std::endl;

        // All particles selected are read.
        h5part_int64_t *indices = new h5part_int64_t[queryResults.size()];
        std::copy(queryResults.begin(), queryResults.end(), indices);
        H5PartSetViewIndices(file, indices, queryResults.size());
        delete[] indices;

        debug5 << "SelectParticlesToRead(): Serial "
               << "all particles (" << queryResults.size()
               << ") in the query read." << std::endl;
    }
    else
    {
#endif      
        // Not using FastBit or no selection
        debug5 << "SelectParticlesToRead(): Serial "
               << "FastBit disabled or no selection active" << std::endl;

        // Not using FastBit or no selection. Read the whole file.
        H5PartSetView(file, -1, -1);

        h5part_int64_t nParticles = H5PartGetNumParticles(file);

        debug5 << "SelectParticlesToRead(): Serial "
               << " reading all particles (" << nParticles
               << ") in file (no query)." << std::endl;
        
#ifdef HAVE_LIBFASTQUERY
    }
#endif // HAVE_LIBFASTQUERY

#endif // PARALLEL

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


// ****************************************************************************
//  Method: avtH5PartFileFormat::GetMesh
//
//  Purpose:
//      Gets the mesh associated with this file.  The mesh is returned as a
//      derived type of vtkDataSet (ie vtkRectilinearGrid, vtkStructuredGrid,
//      vtkUnstructuredGrid, etc).
//
//  Arguments:
//      timestate   The index of the timestate.  If GetNTimesteps returned
//                  'N' time steps, this is guaranteed to be between 0 and N-1.
//      meshname    The name of the mesh of interest.  This can be ignored if
//                  there is only one mesh.
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Tue Feb 9 13:44:50 PST 2010
//
//  Modifications:
//
// ****************************************************************************

vtkDataSet *
avtH5PartFileFormat::GetMesh(int timestate, const char *meshname)
{
    if (strcmp(meshname, "particles") == 0)
    {
        return GetParticleMesh(timestate); // There is only one particle mesh
    }
    else
    {
        return GetFieldMesh(timestate, meshname); // In the future there may be
        // more than one mesh for field data
    }
}

// ****************************************************************************
//  Method: avtH5PartFileFormat::GetParticleMesh
//
//  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.
//      meshname    The name of the mesh of interest.  This can be ignored if
//                  there is only one mesh.
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Tue Feb 9 13:44:50 PST 2010
//
//  Modifications:
//    Brad Whitlock, Mon Dec 12 16:06:13 PST 2011
//    Construct avtOriginalCellNumbers here if we're doing a data selection.
//
// ****************************************************************************

vtkDataSet *
avtH5PartFileFormat::GetParticleMesh(int timestate)
{
    int t1 = visitTimer->StartTimer();

    // Switch to appropriate time step (opens file if necessary)
    ActivateTimestep(timestate);

    idVariableName = defaultIdVariableName;

    VarNameToInt64Map_t::const_iterator it =
      particleVarNameToTypeMap.find(idVariableName);
    
    if (it == particleVarNameToTypeMap.end())
      EXCEPTION1(InvalidVariableException, idVariableName);
    
    // Select particles to actually read
    SelectParticlesToRead( idVariableName.c_str() );
    
    // Read data
    h5part_int64_t nPoints = H5PartGetNumParticles(file);

    // If we do domain decomposition or particle selection, this will only
    // return the number of particles in the current view, i.e., the number of
    // particles we actually read
    debug5 << "GetParticleMesh(): Reading " << nPoints
           << " particles." << std::endl;

    if (nPoints ==  0)
    {
        visitTimer->StopTimer(t1, "H5PartFileFormat::GetParticleMesh()");
        return 0;
    }

    const char *coordNames[3][3] = {
        { "x", "y", "z" }, {"r", "phi", "z" } , { "r", "phi", "theta" }
    };

    h5part_int64_t coordValType = particleVarNameToTypeMap[coordNames[coordType][0]]; 
    if (coordValType !=
        (h5part_int64_t) particleVarNameToTypeMap[coordNames[coordType][1]] ||
        ((particleNSpatialDims > 2) &&
         coordValType != (h5part_int64_t) particleVarNameToTypeMap[coordNames[coordType][2]])) /// TODO: check parenthsis placement logic 
    {
        EXCEPTION2(NonCompliantFileException, "H5Part GetParticleMesh",
                   "Coordinate data sets do not have the same value type "
                   "(double, float, ...).");
    }

    void* coordsRead[3] = { 0, 0, 0 };

    for (int i=0; i < particleNSpatialDims; ++i)
    {
        h5part_int64_t status = H5PART_SUCCESS; 

        if (coordValType == H5PART_INT64)
        {
            coordsRead[i] = malloc(sizeof(h5part_int64_t) * nPoints);
            status = H5PartReadDataInt64(file, coordNames[coordType][i],
                    static_cast<h5part_int64_t*>(coordsRead[i]));
        }
        else if (coordValType == H5PART_INT32)
        {
            coordsRead[i] = malloc(sizeof(h5part_int32_t) * nPoints);
            status = H5PartReadDataInt32(file, coordNames[coordType][i],
                    static_cast<h5part_int32_t*>(coordsRead[i]));
        }
        else if (coordValType == H5PART_FLOAT64)
        {
            coordsRead[i] = malloc(sizeof(h5part_float64_t) * nPoints);
            status = H5PartReadDataFloat64(file, coordNames[coordType][i],
                    static_cast<h5part_float64_t*>(coordsRead[i]));

        }
        else if (coordValType == H5PART_FLOAT32)
        {
            coordsRead[i] = malloc(sizeof(h5part_float32_t) * nPoints);
            status = H5PartReadDataFloat32(file, coordNames[coordType][i],
                    static_cast<h5part_float32_t*>(coordsRead[i]));
        }
        else
        {
            EXCEPTION2(NonCompliantFileException, "H5Part GetParticleMesh",
                       "Unsupported value type for coordinates. "
                       "(Supported tyes are FLOAT64, FLOAT32 and INT64.)");
        }
        // FIXME: Possibly handle other data types for coordinates.

        if (status != H5PART_SUCCESS)
        {
            EXCEPTION2(NonCompliantFileException, "H5Part GetParticleMesh",
                       "Could not read coordinates.");
        }
    }

    // Construct VTK data set
    vtkUnstructuredGrid *dataset = vtkUnstructuredGrid::New();
    vtkPoints *vtkpoints = vtkPoints::New();

    // FIXME: Once VisIt supports dobule data, we need to ensure that
    // native precsission is preserved
    vtkpoints->SetNumberOfPoints((vtkIdType) nPoints);

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

    if (coordValType == H5PART_INT64)
    {
        for (h5part_int64_t ptNo = 0; ptNo < nPoints; ++ptNo)
        {
            pts[3*ptNo] =
                float(static_cast<h5part_int64_t*>(coordsRead[0])[ptNo]);
            pts[3*ptNo+1] =
                float(static_cast<h5part_int64_t*>(coordsRead[1])[ptNo]);
            if (particleNSpatialDims > 2)
                pts[3*ptNo+2] =
                    float(static_cast<h5part_int64_t*>(coordsRead[2])[ptNo]);
            else
                pts[3*ptNo+2] = 0.0;
        }
    }
    else if (coordValType == H5PART_INT32)
    {
        for (h5part_int32_t ptNo = 0; ptNo < nPoints; ++ptNo)
        {
            pts[3*ptNo] =
                float(static_cast<h5part_int32_t*>(coordsRead[0])[ptNo]);
            pts[3*ptNo+1] =
                float(static_cast<h5part_int32_t*>(coordsRead[1])[ptNo]);
            if (particleNSpatialDims > 2)
                pts[3*ptNo+2] =
                    float(static_cast<h5part_int32_t*>(coordsRead[2])[ptNo]);
            else
                pts[3*ptNo+2] = 0.0;
        }
    }
    else if (coordValType == H5PART_FLOAT64)
    {
        for (h5part_int64_t ptNo = 0; ptNo < nPoints; ++ptNo)
        {
            pts[3*ptNo] =
                float(static_cast<h5part_float64_t*>(coordsRead[0])[ptNo]);
            pts[3*ptNo+1] =
                float(static_cast<h5part_float64_t*>(coordsRead[1])[ptNo]);
            if (particleNSpatialDims > 2)
                pts[3*ptNo+2] =
                    float(static_cast<h5part_float64_t*>(coordsRead[2])[ptNo]);
            else
                pts[3*ptNo+2] = 0.0;
        }
    }
    else if (coordValType == H5PART_FLOAT32)
    {
        for (h5part_int64_t ptNo = 0; ptNo < nPoints; ptNo++)
        {
            pts[3*ptNo] =
                float(static_cast<h5part_float32_t*>(coordsRead[0])[ptNo]);
            pts[3*ptNo+1] =
                float(static_cast<h5part_float32_t*>(coordsRead[1])[ptNo]);
            if (particleNSpatialDims > 2)
                pts[3*ptNo+2] =
                    float(static_cast<h5part_float32_t*>(coordsRead[2])[ptNo]);
            else
                pts[3*ptNo+2] = 0.0;
        }
    }
    else
    {
        EXCEPTION2(NonCompliantFileException, "H5Part GetParticleMesh",
                   "Unsupported value type for coordinates. "
                   "Supported types are FLOAT64. FLOAT32 and INT64.");
    }
    // FIXME: Possibly handle other data types for the coordinates

    for (int i=0; i <particleNSpatialDims; ++i)
        free(coordsRead[i]);

    dataset->Allocate(nPoints);

    for (h5part_int64_t i=0; i<nPoints; ++i)
    {
        vtkIdType onevertex = (vtkIdType) i;
        dataset->InsertNextCell(VTK_VERTEX, 1, &onevertex);
    }

    dataset->SetPoints(vtkpoints);
    vtkpoints->Delete();

    //
    // If a data selection was processed successfully then dataset
    // with different connectivity from the original full dataset will
    // be produced.  This means that if we're requesting original cell
    // numbers above up in the database then the values calculated
    // there would be wrong. We need to calculate them here and send
    // them along so named selections will work properly.
    //
    debug5 << "avtH5PartFileFormat::GetParticleMesh(): "
           << " Creating avtOriginalCellNumbers early" << std::endl;

    vtkUnsignedIntArray *origZones = vtkUnsignedIntArray::New();
    origZones->SetName("avtOriginalCellNumbers");
    origZones->SetNumberOfComponents(2);
    origZones->SetNumberOfTuples(nPoints);
    unsigned int *iptr = (unsigned int *) origZones->GetVoidPointer(0);

    if (it->second == H5PART_FLOAT64)
    {
        h5part_float64_t *idList = new h5part_float64_t[nPoints];;
        if (H5PartReadDataFloat64(file, idVariableName.c_str(), idList) != H5PART_SUCCESS)
        {
          delete[] idList;
          EXCEPTION1(InvalidVariableException, idVariableName);
        }
        else
        {
          for(size_t i = 0; i < nPoints; ++i)
          {
            // domain 0. This format is single domain so it's probably okay.
            *iptr++ = 0;
            *iptr++ = (unsigned int) idList[i];
        
            debug5 << "Domain " << 0 << " cell " << idList[i]
                   << std::endl;
          }
    
          delete[] idList;
        }
    }
    else if (it->second == H5PART_FLOAT32)
    {
        h5part_float32_t *idList = new h5part_float32_t[nPoints];
        if (H5PartReadDataFloat32(file, idVariableName.c_str(), idList) != H5PART_SUCCESS)
        {
          delete[] idList;
          EXCEPTION1(InvalidVariableException, idVariableName);
        }
        else
        {
          for(size_t i = 0; i < nPoints; ++i)
          {
            // domain 0. This format is single domain so it's probably okay.
            *iptr++ = 0;
            *iptr++ = (unsigned int) idList[i];
        
            debug5 << "Domain " << 0 << " cell " << idList[i]
                   << std::endl;
          }

          delete[] idList;
        }
    }
    else if (it->second == H5PART_INT64)
    {
      h5part_int64_t *idList = new h5part_int64_t[nPoints];
      if (H5PartReadDataInt64(file, idVariableName.c_str(), idList) != H5PART_SUCCESS)
      {
          delete[] idList;
          EXCEPTION1(InvalidVariableException, idVariableName);
      }
      else
      {
          for(size_t i = 0; i < nPoints; ++i)
          {
            // domain 0. This format is single domain so it's probably okay.
            *iptr++ = 0;
            *iptr++ = (unsigned int) idList[i];
        
            debug5 << "Domain " << 0 << " cell " << idList[i]
                   << std::endl;
          }

          delete[] idList;
      }
    }
    else if (it->second == H5PART_INT32)
    {
      h5part_int32_t *idList = new h5part_int32_t[nPoints];
      if (H5PartReadDataInt32(file, idVariableName.c_str(), idList) != H5PART_SUCCESS)
      {
          delete[] idList;
          EXCEPTION1(InvalidVariableException, idVariableName);
      }
      else
      {
          for(size_t i = 0; i < nPoints; ++i)
          {
            // domain 0. This format is single domain so it's probably okay.
            *iptr++ = 0;
            *iptr++ = (unsigned int) idList[i];
        
            debug5 << "Domain " << 0 << " cell " << idList[i]
                   << std::endl;
          }

          delete[] idList;
      }
    }
    else
    {
        EXCEPTION2(NonCompliantFileException, "H5Part GetVar",
                   std::string("Scalar variable ") + idVariableName + 
                   std::string(" uses an unsupported value type. Supported "
                               "types are FLOAT64, FLOAT32 and INT64."));
    }
    // FIXME: Possibly handle other data types for the id variable

    dataset->GetCellData()->AddArray(origZones);
    origZones->Delete();
    
    visitTimer->StopTimer(t1, "H5PartFileFormat::GetParticleMesh()");

    // Reset view
    H5PartSetView(file, -1, -1);

    return dataset;
}

// ****************************************************************************
//  Method: avtH5PartFileFormat::GetFieldMesh
//
//  Purpose:
//      Gets the mesh associated with fields of this file
//
//  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.
//
//  Programmer: kurts
//  Creation:   Tue Aug 28 17:35:50 PDT 2007
//
//  Modifications:
//    Gunther H. Weber, Fri Apr 17 13:55:07 PDT 2009
//    Read block origin from file.
//
//    Gunther H. Weber, Fri Jul 17 11:11:36 PDT 2009
//    Clean up. Read block origin via dedicated API call.
//
//    Gunther H. Weber, Wed Feb 10 11:47:32 PST 2010
//    Conversion from MD plugin to SD plugin supporting domain decomposition
//    and clean-up.
//
//    Gunther H. Weber, Thu Apr  1 18:14:00 PDT 2010
//    Adapted to GetSubBlock paramter change.
//
// ****************************************************************************

vtkDataSet *
avtH5PartFileFormat::GetFieldMesh(int timestate, const char *meshname)
{
    int t1 = visitTimer->StartTimer();

    // Switch to appropriate time step (opens file if necessary)
    ActivateTimestep(timestate);

    // FIXME: assume all fields have the same dimensions
    // get information from first field
    char fieldName[maxVarNameLen];
    h5part_int64_t gridRank;
    h5part_int64_t gridDims[3];
    h5part_int64_t fieldDims;
    h5part_int64_t type;
    if (H5BlockGetFieldInfo (file, 0, fieldName, maxVarNameLen,
                &gridRank, gridDims, &fieldDims, &type) != H5PART_SUCCESS)
        EXCEPTION2(NonCompliantFileException, "H5Part GetFieldMesh",
                   "Could not read field information.");

    // Determine which data portion we need to construct a mesh for
    h5part_int64_t subBlockDims[6];
    GetSubBlock(gridDims, subBlockDims);

    double xOrigin = 0;
    double yOrigin = 0;
    double zOrigin = 0;

    if (H5Block3dGetFieldOrigin(file, fieldName,
                &xOrigin, &yOrigin, &zOrigin) != H5PART_SUCCESS)
    {
        EXCEPTION2(NonCompliantFileException, "H5Part GetFieldMesh",
                   "Could not read field origin.");
    }

    debug5 << "xOrigin: " << xOrigin << "\n";
    debug5 << "yOrigin: " << yOrigin << "\n";
    debug5 << "zOrigin: " << zOrigin << "\n";

    double xSpacing = 0;
    double ySpacing = 0;
    double zSpacing = 0;

    if (H5Block3dGetFieldSpacing (file, fieldName,
                &xSpacing, &ySpacing, &zSpacing) != H5PART_SUCCESS)
    {
        EXCEPTION2(NonCompliantFileException, "H5Part GetFieldMesh", 
                   "Could not read field spacing.");
    }

    debug5 << "xSpacing: " << xSpacing << "\n";
    debug5 << "ySpacing: " << ySpacing << "\n";
    debug5 << "zSpacing: " << zSpacing << "\n";

    vtkFloatArray *coords[3] = {0, 0, 0};
    // set x coordinates
    coords[0] = vtkFloatArray::New();
    coords[0]->SetNumberOfTuples(subBlockDims[1] - subBlockDims[0] + 1);
    float *xarray = (float *) coords[0]->GetVoidPointer(0);
    for (int i=subBlockDims[0]; i <= subBlockDims[1]; i++)
    {
        xarray[i-subBlockDims[0]] = xOrigin + i * xSpacing;
    }

    // set y coordinates
    coords[1] = vtkFloatArray::New();
    coords[1]->SetNumberOfTuples(subBlockDims[3] - subBlockDims[2] + 1);
    float *yarray = (float *) coords[1]->GetVoidPointer(0);
    for (int i=subBlockDims[2]; i <= subBlockDims[3]; i++)
    {
        yarray[i-subBlockDims[2]] = yOrigin + i * ySpacing;
    }

    // set z coordinates
    coords[2] = vtkFloatArray::New();
    coords[2]->SetNumberOfTuples(subBlockDims[5] - subBlockDims[4] + 1);
    float *zarray = (float *) coords[2]->GetVoidPointer(0);
    for (int i=subBlockDims[4]; i <= subBlockDims[5]; i++)
    {
        zarray[i-subBlockDims[4]] = zOrigin + i * zSpacing;
    }

    // create vtkRectilinearGrid objects + set dims and coords
    vtkRectilinearGrid *rgrid = vtkRectilinearGrid::New();
    rgrid->SetDimensions(
            coords[0]->GetNumberOfTuples(),
            coords[1]->GetNumberOfTuples(),
            coords[2]->GetNumberOfTuples());

    rgrid->SetXCoordinates(coords[0]);
    coords[0]->Delete();
    rgrid->SetYCoordinates(coords[1]);
    coords[1]->Delete();
    rgrid->SetZCoordinates(coords[2]);
    coords[2]->Delete();

    visitTimer->StopTimer(t1, "H5PartFileFormat::GetFieldMesh()");

    return rgrid;
}

// ****************************************************************************
//  Method: avtH5PartFileFormat::GetVar
//
//  Purpose:
//      Gets a scalar variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      timestate  The index of the timestate.  If GetNTimesteps returned
//                 'N' time steps, this is guaranteed to be between 0 and N-1.
//      varName    The name of the variable requested.
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Tue Feb 9 13:44:50 PST 2010
//
// ****************************************************************************

vtkDataArray *
avtH5PartFileFormat::GetVar(int timestate, const char *varName)
{
    int t1 = visitTimer->StartTimer();

    VarNameToInt64Map_t::const_iterator it =
      particleVarNameToTypeMap.find(varName);

    if (it == particleVarNameToTypeMap.end())
    {
        visitTimer->StopTimer(t1, "H5PartFileFormat::GetVar() [invoking GetFieldVar]");
        return GetFieldVar(timestate, varName);
    }

    // Switch to appropriate time step (opens file if necessary)
    ActivateTimestep(timestate);

    // Select particles to actually read
    SelectParticlesToRead( varName );

    // Read data
    h5part_int64_t nPoints = H5PartGetNumParticles(file);

    debug5 << "GetVar() reading " << nPoints << " particles." << std::endl;

    if (nPoints == 0)
    {
        visitTimer->StopTimer(t1, "H5PartFileFormat::GetVar()");
        return 0;
    }

    vtkDataArray *scalars = 0;
    if (it->second == H5PART_FLOAT64)
    {
        scalars = vtkDoubleArray::New();
        scalars->SetNumberOfTuples(nPoints);
        if (H5PartReadDataFloat64(file, varName,
                    static_cast<h5part_float64_t*>(scalars->GetVoidPointer(0)))
                != H5PART_SUCCESS)
        {
            EXCEPTION1(InvalidVariableException, varName);
        }
    }
    else if (it->second == H5PART_FLOAT32)
    {
        scalars = vtkFloatArray::New();
        scalars->SetNumberOfTuples(nPoints);
        if (H5PartReadDataFloat32(file, varName,
                    static_cast<h5part_float32_t*>(scalars->GetVoidPointer(0)))
                != H5PART_SUCCESS)
        {
            EXCEPTION1(InvalidVariableException, varName);
        }
    }
    else if (it->second == H5PART_INT64)
    {
        // Unlike FLOAT32 and FLOAT64 which should map to float and
        // double we are not sure what int type is used to achieve 64 bits.
        // To be safe, read into separate array and copy.
        h5part_int64_t *idvar = new h5part_int64_t[nPoints];
        if (H5PartReadDataInt64(file, varName, idvar) != H5PART_SUCCESS)
            EXCEPTION1(InvalidVariableException, varName);

        scalars = vtkFloatArray::New(); // FIXME: Can an integer array be used?
        scalars->SetNumberOfTuples(nPoints);
        float *ptr = static_cast<float*>(scalars->GetVoidPointer(0));
        for (h5part_int64_t i=0; i<nPoints; ++i)
        {
            ptr[i] = float(idvar[i]);
        }
        delete[] idvar;
    }
    else if (it->second == H5PART_INT32)
    {
        // Unlike FLOAT32 and FLOAT64 which should map to float and
        // double we are not sure what int type is used to achieve 32 bits.
        // To be safe, read into separate array and copy.
        h5part_int32_t *idvar = new h5part_int32_t[nPoints];
        if (H5PartReadDataInt32(file, varName, idvar) != H5PART_SUCCESS)
            EXCEPTION1(InvalidVariableException, varName);

        scalars = vtkFloatArray::New(); // FIXME: Can an integer array be used?
        scalars->SetNumberOfTuples(nPoints);
        float *ptr = static_cast<float*>(scalars->GetVoidPointer(0));
        for (h5part_int32_t i=0; i<nPoints; ++i)
        {
            ptr[i] = float(idvar[i]);
        }
        delete[] idvar;
    }
    else
    {
        EXCEPTION2(NonCompliantFileException, "H5Part GetVar",
                   std::string("Scalar variable ") + varName + 
                   std::string(" uses an unsupported value type. Supported "
                               "types are FLOAT64, FLOAT32 and INT64."));
    }
    // FIXME: Possibly handle other data types for scalar variables

    // Reset view
    H5PartSetView(file, -1, -1);

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

    return scalars;
}


// ****************************************************************************
//  Method: avtH5PartFileFormat::GetFieldVar
//
//  Purpose:
//      Get a scalar variable associated with the field mesh of this file
//
//  Arguments:
//      timestate   The index of the timestate.  If GetNTimesteps returned
//                  'N' time steps, this is guaranteed to be between 0 and N-1.
//      varName    The name of the variable requested.
//
//  Programmer: kurts
//  Creation:   Tue Aug 28 17:35:50 PDT 2007
//
//  Modifications:
//    Gunther H. Weber, Fri Jul 17 11:13:39 PDT 2009
//    Clean-up.
//
//    Gunther H. Weber, Wed Feb 10 11:47:32 PST 2010
//    Conversion from MD plugin to SD plugin supporting domain decomposition.
//
//    Gunther H. Weber, Thu Apr  1 18:14:00 PDT 2010
//    Adapted to GetSubBlock paramter change.
//
// ****************************************************************************

vtkDataArray *
avtH5PartFileFormat::GetFieldVar(int timestate, const char* varName)
{
    int t1 = visitTimer->StartTimer();

    // Check if we know about variable and get type
    VarNameToInt64Map_t::const_iterator it =
      fieldScalarVarNameToTypeMap.find(varName);

    if (it == fieldScalarVarNameToTypeMap.end())
        EXCEPTION1(InvalidVariableException, varName);

    // Activate correct time step (opens file if necessary)
    ActivateTimestep(timestate);

    char fieldName[maxVarNameLen];
    h5part_int64_t gridRank;
    h5part_int64_t gridDims[3];
    h5part_int64_t fieldDims;
    h5part_int64_t type;

    // FIXME: This assumes that all blocks have the same mesh
    if (H5BlockGetFieldInfo (file, 0, fieldName, maxVarNameLen,
                &gridRank, gridDims, &fieldDims, &type) != H5PART_SUCCESS)
        EXCEPTION1(InvalidVariableException,
                   "Could not read field information.");

    h5part_int64_t subBlockDims[6];
    GetSubBlock(gridDims, subBlockDims);

    // set field layout
    if(H5BlockDefine3DFieldLayout(file, subBlockDims[0], subBlockDims[1], 
                subBlockDims[2], subBlockDims[3], 
                subBlockDims[4], subBlockDims[5]) != H5PART_SUCCESS)
        EXCEPTION1(VisItException, "Could not set field layout.");

    // read data
    int nValues = (subBlockDims[1] - subBlockDims[0] + 1) *
        (subBlockDims[3] - subBlockDims[2] + 1) *
        (subBlockDims[5] - subBlockDims[4] + 1);

    vtkDataArray *scalars;
    if (it->second == H5PART_FLOAT64)
    {
        scalars = vtkDoubleArray::New();
        scalars->SetNumberOfTuples(nValues);
        if (H5Block3dReadScalarFieldFloat64(file, varName,
                    static_cast<h5part_float64_t*>(scalars->GetVoidPointer(0)))
                != H5PART_SUCCESS)
        {
            EXCEPTION1(InvalidVariableException, varName);
        }
    }
    else if (it->second == H5PART_FLOAT32)
    {
        scalars = vtkFloatArray::New();
        scalars->SetNumberOfTuples(nValues);
        if (H5Block3dReadScalarFieldFloat32(file, varName,
                    static_cast<h5part_float32_t*>(scalars->GetVoidPointer(0)))
                != H5PART_SUCCESS)
        {
            EXCEPTION1(InvalidVariableException, varName);
        }
    }
    else if (it->second == H5PART_INT64)
    {
        // Unlike FLOAT32 and FLOAT64 which should map to float and
        // double we are not sure what int type is used to achieve 64 bits.
        // To be safe, read into separate array and copy.
        h5part_int64_t *idvar = new h5part_int64_t[nValues];
        if (H5Block3dReadScalarFieldInt64(file, varName, idvar) != H5PART_SUCCESS)
            EXCEPTION1(InvalidVariableException, varName);

        scalars = vtkFloatArray::New(); // FIXME: Can an integer array be used?
        scalars->SetNumberOfTuples(nValues);
        float *ptr = static_cast<float*>(scalars->GetVoidPointer(0));
        for (h5part_int64_t i=0; i<nValues; ++i)
        {
            ptr[i] = float(idvar[i]);
        }
        delete[] idvar;
    }
    else if (it->second == H5PART_INT32)
    {
        // Unlike FLOAT32 and FLOAT64 which should map to float and
        // double we are not sure what int type is used to achieve 32 bits.
        // To be safe, read into separate array and copy.
        h5part_int32_t *idvar = new h5part_int32_t[nValues];
        if (H5Block3dReadScalarFieldInt32(file, varName, idvar) != H5PART_SUCCESS)
            EXCEPTION1(InvalidVariableException, varName);

        scalars = vtkFloatArray::New(); // FIXME: Can an integer array be used?
        scalars->SetNumberOfTuples(nValues);
        float *ptr = static_cast<float*>(scalars->GetVoidPointer(0));
        for (h5part_int32_t i=0; i<nValues; ++i)
        {
            ptr[i] = float(idvar[i]);
        }
        delete[] idvar;
    }
    else
    {
        EXCEPTION2(NonCompliantFileException, "H5Part GetVar",
                   std::string("Field variable ") + varName + 
                   std::string(" uses an unsupported value type. Supported "
                               "types are FLOAT64, FLOAT32 and INT64."));
    }
    // FIXME: Possibly handle other data types for field variables

    visitTimer->StopTimer(t1, "H5PartFileFormat::GetFieldVar()");

    return scalars;
}

// ****************************************************************************
//  Method: avtH5PartFileFormat::GetVectorVar
//
//  Purpose:
//      Gets a vector variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      timestate  The index of the timestate.  If GetNTimesteps returned
//                 'N' time steps, this is guaranteed to be between 0 and N-1.
//      varName    The name of the variable requested.
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Tue Feb 9 13:44:50 PST 2010
//
//  Modifications:
//    Kurt Stockinger, Tue Aug 28 17:35:50 PDT 2007
//    Added support for field data
//
//    Gunther H. Weber, Fri Jul 17 11:13:39 PDT 2009
//    Clean-up.
//
//    Gunther H. Weber, Thu Apr  1 18:14:00 PDT 2010
//    Adapted to GetSubBlock paramter change.
//
// ****************************************************************************
vtkDataArray *
avtH5PartFileFormat::GetVectorVar(int timestate, const char *varName)
{
    int t1 = visitTimer->StartTimer();

    // Check if we know about variable and get type
    VarNameToInt64Map_t::const_iterator it =
      fieldVectorVarNameToTypeMap.find(varName);
    
    if (it == fieldVectorVarNameToTypeMap.end())
        EXCEPTION1(InvalidVariableException, varName);

    // Activate correct time step (opens file if necessary)
    ActivateTimestep(timestate);

    char fieldName[maxVarNameLen];
    h5part_int64_t gridRank;
    h5part_int64_t gridDims[3];
    h5part_int64_t fieldDims;
    h5part_int64_t type;

    // FIXME: This assumes that all blocks have the same size
    if (H5BlockGetFieldInfo (file, 0, fieldName, maxVarNameLen,
                &gridRank, gridDims, &fieldDims, &type) != H5PART_SUCCESS)
        EXCEPTION1(InvalidVariableException, "Could not read field information.");

    h5part_int64_t subBlockDims[6];
    GetSubBlock(gridDims, subBlockDims);

    // set field layout
    if (H5BlockDefine3DFieldLayout (file, 
                subBlockDims[0], subBlockDims[1], 
                subBlockDims[2], subBlockDims[3], 
                subBlockDims[4], subBlockDims[5]) != H5PART_SUCCESS)
        EXCEPTION1(InvalidVariableException, "Could not set field layout");

    // Read data
    int nValues = (subBlockDims[1] - subBlockDims[0] + 1) *
        (subBlockDims[3] - subBlockDims[2] + 1) *
        (subBlockDims[5] - subBlockDims[4] + 1);

    void* vecsRead[3] = { 0, 0, 0 };
    h5part_int64_t status = H5PART_SUCCESS; (void) status;
    if (it->second == H5PART_INT64)
    {
        for (int i=0; i<3; ++i)
            vecsRead[i] = malloc(sizeof(h5part_int64_t) * nValues);

        status = H5Block3dRead3dVectorFieldInt64(file, varName,  
                static_cast<h5part_int64_t*>(vecsRead[0]),
                static_cast<h5part_int64_t*>(vecsRead[1]),
                static_cast<h5part_int64_t*>(vecsRead[2]));
    }
    else if (it->second == H5PART_INT32)
    {
        for (int i=0; i<3; ++i)
            vecsRead[i] = malloc(sizeof(h5part_int32_t) * nValues);

        status = H5Block3dRead3dVectorFieldInt32(file, varName,  
                static_cast<h5part_int32_t*>(vecsRead[0]),
                static_cast<h5part_int32_t*>(vecsRead[1]),
                static_cast<h5part_int32_t*>(vecsRead[2]));
    }
    else if (it->second == H5PART_FLOAT64)
    {
        for (int i=0; i<3; ++i)
            vecsRead[i] = malloc(sizeof(h5part_float64_t) * nValues);

        status = H5Block3dRead3dVectorFieldFloat64(file, varName,  
                static_cast<h5part_float64_t*>(vecsRead[0]),
                static_cast<h5part_float64_t*>(vecsRead[1]),
                static_cast<h5part_float64_t*>(vecsRead[2]));


    }
    else if (it->second == H5PART_FLOAT32)
    {
        for (int i=0; i<3; ++i)
            vecsRead[i] = malloc(sizeof(h5part_float32_t) * nValues);

        status = H5Block3dRead3dVectorFieldFloat32(file, varName,  
                static_cast<h5part_float32_t*>(vecsRead[0]),
                static_cast<h5part_float32_t*>(vecsRead[1]),
                static_cast<h5part_float32_t*>(vecsRead[2]));
    }
    else
    {
        EXCEPTION2(NonCompliantFileException, "H5Part GetVectorVar",
                   std::string("Vector variable ") + varName +
                   std::string(" uses an unsupported value type. Supported"
                               "types are FLOAT64, FLOAT32 and INT64."));
    }
    // FIXME: Possibly handle other data types for vector variables

    // Store data in VTK object
    vtkDataArray *array;
    if (it->second == H5PART_FLOAT64)
    {
        array = vtkDoubleArray::New();
    }
    else
    {
        array = vtkFloatArray::New(); //TODO: check fix
        // FIXME: Possibly check for integer and create an integer array
    }

    // Set array size
    array->SetNumberOfComponents(3);
    array->SetNumberOfTuples(nValues);

    // Copy data
    if (it->second == H5PART_INT64)
    {
        // FIXME: Possibly check for integer and create an integer array
        float *vecs = static_cast<float*>(array->GetVoidPointer(0));
        for (h5part_int64_t valNo = 0; valNo < nValues; valNo++)
        {
            vecs[3*valNo] =
                float(static_cast<h5part_int64_t*>(vecsRead[0])[valNo]);
            vecs[3*valNo+1] =
                float(static_cast<h5part_int64_t*>(vecsRead[1])[valNo]);
            vecs[3*valNo+2] =
                float(static_cast<h5part_int64_t*>(vecsRead[2])[valNo]);
        }
    }
    else if (it->second == H5PART_INT32)
    {
        // FIXME: Possibly check for integer and create an integer array
        float *vecs = static_cast<float*>(array->GetVoidPointer(0));
        for (h5part_int32_t valNo = 0; valNo < nValues; valNo++)
        {
            vecs[3*valNo] =
                float(static_cast<h5part_int32_t*>(vecsRead[0])[valNo]);
            vecs[3*valNo+1] =
                float(static_cast<h5part_int32_t*>(vecsRead[1])[valNo]);
            vecs[3*valNo+2] =
                float(static_cast<h5part_int32_t*>(vecsRead[2])[valNo]);
        }
    }
    else if (it->second == H5PART_FLOAT64)
    {
        double *vecs = static_cast<double*>(array->GetVoidPointer(0));
        for (h5part_int64_t valNo = 0; valNo < nValues; valNo++)
        {
            vecs[3*valNo] =
                double(static_cast<h5part_float64_t*>(vecsRead[0])[valNo]);
            vecs[3*valNo+1] =
                double(static_cast<h5part_float64_t*>(vecsRead[1])[valNo]);
            vecs[3*valNo+2] =
                double(static_cast<h5part_float64_t*>(vecsRead[2])[valNo]);
        }
    }
    else if (it->second == H5PART_FLOAT32)
    {
        float *vecs = static_cast<float*>(array->GetVoidPointer(0));
        for (h5part_int64_t valNo = 0; valNo < nValues; valNo++)
        {
            vecs[3*valNo] =
                float(static_cast<h5part_float32_t*>(vecsRead[0])[valNo]);
            vecs[3*valNo+1] =
                float(static_cast<h5part_float32_t*>(vecsRead[1])[valNo]);
            vecs[3*valNo+2] =
                float(static_cast<h5part_float32_t*>(vecsRead[2])[valNo]);
        }
    }
    else
    {
        EXCEPTION2(NonCompliantFileException, "H5Part GetVectorVar",
                   std::string("Vector variable ") + varName +
                   std::string(" uses an unsupported value type. Supported "
                               "types are FLOAT64, FLOAT32 and INT64.)"));
    }
    // FIXME: Possibly handle other data types for vector variables
 
    for (int i=0; i<3; ++i)
      free(vecsRead[i]);

    visitTimer->StopTimer(t1, "H5PartFileFormat::GetVectorVar()");
    return array;
}

// ****************************************************************************
//  Method: avtH5PartFileFormat::ActivateTimestep
//
//  Purpose:
//      Select the current time step
//
//  Arguments:
//      ts      The new current time step
//
//  Output:
//
//  Programmer: ghweber
//  Creation:   Thu Feb 11 15:38:49 PST 2010
//
//  Modifications:
//    Gunther H. Weber, Mon Aug  9 14:58:58 PDT 2010
//    Fixed bug identified by regression tests that sometimes (it seems
//    mainly for field data) H5PartSetStep is not called when a new file
//    is opened.
//
//    Gunther H. Weber, Wed Aug 17 13:44:13 PDT 2011
//    Do not open file with HDF5_FQ if useFastBitIndex is false.
//
// ****************************************************************************

void
avtH5PartFileFormat::ActivateTimestep(int ts)
{
    int t1 = visitTimer->StartTimer();

    // Open file if necessary
    if (!file)
    {
        debug5 << "avtH5PartFileFormat::ActivateTimestep(): Opening file and ";
        debug5 << "activating time step " << ts << std::endl;

        file = H5PartOpenFile(filenames[0], H5PART_READ);

        if (!file)
          EXCEPTION1(InvalidFilesException, "Cannot open file.");
        
        if (H5PartSetStep(file, ts) != H5PART_SUCCESS)
        {
          std::ostringstream msg;
          msg << "Cannot activate time step " << ts << ".";

          debug1 << "avtH5PartFileFormat::ActivateTimestep(): "
                 << msg.str() << std::endl;

          EXCEPTION2(NonCompliantFileException, "H5Part AcitvateTimestep",
                     msg.str());
        }

#ifdef HAVE_LIBFASTQUERY
        // New file is opened -> Need to update query results.
        queryResultsValid = false;
#endif
    } 
    else if (ts != activeTimeStep)
    {
        debug5 << "avtH5PartFileFormat::ActivateTimestep(): Activating time ";
        debug5 << "step " << ts << std::endl;

        if (H5PartSetStep(file, ts) != H5PART_SUCCESS)
        {
          std::ostringstream msg;
          msg << "Cannot activate time step " << ts << ".";

          debug1 << "avtH5PartFileFormat::ActivateTimestep(): "
                 << msg.str() << std::endl;

          EXCEPTION2(NonCompliantFileException, "H5Part AcitvateTimestep",
                     msg.str());
        }

#ifdef HAVE_LIBFASTQUERY
        // New time step -> Need to update query results.
        queryResultsValid = false;
#endif
    }

    // Update global time step information
    activeTimeStep = ts;

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

// ****************************************************************************
//  Method: avtH5PartFileFormat::GetSubBlock
//
//  Purpose:
//      Compute extents of a sub bliock for a given processor id
//
//  Arguments:
//      gridDims      Number of nodes in x-/y-/z-direction
//
//  Output:
//      subBlockDims  Range of samples that this processor should read. Order
//                    i_min, i_max, j_min, j_max, k_min, k_max. Indices are
//                    inclusive, i.e., processor should also read samples i_min
//                    and i_max as well as all samples with indices in between.
//
//  Programmer: kurts and ghweber
//  Creation:   Tue Aug 28 17:35:50 PDT 2007
//
//  Modifications:
//    Gunther H. Weber, Thu Apr  1 15:43:05 PDT 2010
//    Use PAR_Rank() instead of block number as parameter. Fix problem
//    that cuased data to be only read partially if domain
//    decomposition was disabled.
//
// ****************************************************************************

void avtH5PartFileFormat::GetSubBlock(h5part_int64_t gridDims[3],
        h5part_int64_t subBlockDims[6])
{
    int t1 = visitTimer->StartTimer();
#ifdef PARALLEL
    if (enableDomainDecomposition)
    {
        int partitionAxis =0;
        int axis1 = 1;
        int axis2 = 2;

        if (gridDims[1] > gridDims[0] && gridDims[1] > gridDims[2])
        {
            partitionAxis = 1;
            axis1 = 0;
            axis2 = 2;
        }
        else if (gridDims[2] > gridDims[0] && gridDims[2] > gridDims[1])
        {
            axis1 = 0;
            axis2 = 1;
            partitionAxis = 2;
        }
        
        debug1 << "Grid dims are "
               << gridDims[0] << " " << gridDims[1] << " " << gridDims[2]
               << " partition axis is " << partitionAxis
               << " other axes are " << axis1 << " " << axis2 << std::endl;

        subBlockDims[2*axis1+0] = 0;
        subBlockDims[2*axis1+1] = gridDims[axis1] - 1;

        subBlockDims[2*axis2+0] = 0;
        subBlockDims[2*axis2+1] = gridDims[axis2] - 1;

        int numPartsForPartition = std::min(h5part_int64_t(PAR_Size()),
                                            gridDims[partitionAxis]);
        float slicesPerRank =
          float(gridDims[partitionAxis]) / float(numPartsForPartition);
        
        debug1 << "Slices per rank: " << slicesPerRank << std::endl;

        subBlockDims[2*partitionAxis+0] = int(PAR_Rank()*slicesPerRank);
        subBlockDims[2*partitionAxis+1] = int((PAR_Rank()+1)*slicesPerRank);
        
        if (PAR_Rank() == numPartsForPartition - 1)
        {
            // Make sure that the final slice is included and do not
            // include anything beyond.
            subBlockDims[2*partitionAxis+1] = gridDims[partitionAxis] - 1; 
        }
        
        debug1 << "Rank " << PAR_Rank() << " has block "
               << subBlockDims[0] << " " << subBlockDims[1] << " "
               << subBlockDims[2] << " " << subBlockDims[3] << " "
               << subBlockDims[4] << " " << subBlockDims[5] << std::endl;
    }
    else
    {
        debug1 << "Domain decomposition disabled. Returning entire block."
               << std::endl;
#endif
        subBlockDims[2*0+0] = 0;
        subBlockDims[2*0+1] = gridDims[0] - 1;

        subBlockDims[2*1+0] = 0;
        subBlockDims[2*1+1] = gridDims[1] - 1;

        subBlockDims[2*2+0] = 0;
        subBlockDims[2*2+1] = gridDims[2] - 1;

        debug1 << "avtH5PartFileFormat::GetSubBlock() returning block "
               << "comprising entire data set:"
               << subBlockDims[0] << " " << subBlockDims[1] << " "
               << subBlockDims[2] << " " << subBlockDims[3] << " "
               << subBlockDims[4] << " " << subBlockDims[5] << std::endl;
#ifdef PARALLEL
    }
#endif    

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

// ****************************************************************************
//  Method: avtH5PartFileFormat::ConstructHistogram
//
//  Purpose:
//      Given a specification, fill in the histogram data through FastBit calls
//
//  Arguments:
//      spec   The input specification
//
//  Returns:    Specification with count/bounds data filled in. Need to
//              determine appropriate exceptions to throw
//
//  Programmer: Prabhat
//  Creation:   February 26, 2008
//
//  Modifications:
//
//    Gunther H. Weber, Mon Feb 15 20:20:56 PST 2010
//    Merged into new H5Part database plugin.
//
//    Gunther H. Weber, Wed Aug  4 18:45:06 PDT 2010
//    Some clean-up and reformatting.
//
// ****************************************************************************
#ifdef HAVE_LIBFASTQUERY

void
avtH5PartFileFormat::ConstructHistogram(avtHistogramSpecification *spec)
{
    std::string method = "avtH5PartFileFormat::ConstructHistogram(): ";

    // NOTE: Only called from GetAuxiliaryData() if useFastBitIndex is true.

    if (!spec) 
        EXCEPTION1(ImproperUseException, "NULL Histogram Specification");

    // Compute the histogram
    int                       timestep = spec->GetTimestep();
    bool                regularBinning = spec->IsRegularBinning();
    std::vector<std::string> variables = spec->GetVariables();
    std::vector<int>           numBins = spec->GetNumberOfBins();
    bool               boundsSpecified = spec->BoundsSpecified();
    std::string              condition = spec->GetCondition();

    int                     boundsSize = spec->GetBounds().size();
    //VISIT_LONG_LONG *           counts = spec->GetCounts();
    FastQuery::QueryProcessor qp =
      FastQuery::QueryProcessor( filenames[0], FastQuery::FQ_H5Part,
                                 filenames[0] );

    qp.setVariablePathPrefix( getVariablePathPrefix(activeTimeStep).c_str() );
    qp.setIndexPathPrefix( getFastBitIndexPathPrefix(timestep).c_str() );
      
    for (size_t i=0; i<numBins.size(); i++) 
    {
        if (regularBinning)
        {
            if (numBins[i]<1)
            {
                debug5 << method
                       <<  "Need to specify a valid #bins for histogram!"
                       << std::endl;
                EXCEPTION1(ImproperUseException,
                           "numBins cannot be less than one.");
            }
        }
    }

    if (variables.empty())
    {
        debug5 << method
               << "Histogram Spec needs to specify variables to plot!"
               << std::endl;
        EXCEPTION1(ImproperUseException,
                   "No variables for histogram specified.");
    }

    for (size_t i=0; i<variables.size(); i++)
    {
        if( 0 && qp.checkForIndex( variables[i] ) == false )
        {
            debug5 << method
                   << "Histogram Spec lists variable not present in file!"
                   << std::endl;
            EXCEPTION1(InvalidVariableException, variables[i]);
        }
    }

    if (boundsSpecified)
    {
        if (variables.size() != (size_t)boundsSize)
        {
            debug5 << method
                   << " Histogram Spec variable list does not match bounds "
                   << " list size " << std::endl;

            EXCEPTION1(ImproperUseException, "Variable bounds list mismatch.");
        }
    }

    // Try to retrieve the histogram from cache
    bool isCached = histoCache.getCached(spec);

    if (isCached)
    {
        // FIXME: Does it really make sense to override boundsSpecified?
        boundsSpecified = spec->BoundsSpecified();  
        debug5 << method << "Histogram retrieved from cache!" << std::endl;
    }  

    // Define the begin and end 
    std::vector<double> begins;
    std::vector<double> ends;
    begins.resize(variables.size());
    ends.resize(variables.size());

    // If bounds are specified then copy the extends into the histogram
    if (boundsSpecified)
    {
        for (size_t i=0; i<variables.size(); i++)
        {
            begins[i] = spec->GetBounds()[i][0];
            ends[i]   = spec->GetBounds()[i][spec->GetBounds()[i].size()-1];

            debug5 << method
                   << "boundsSpecified to be "<< begins[i] << ", " << ends[i]
                   << " for " << variables[i] 
                   << std::endl;
        }
    }
    // If bounds are not specified than ask the reader for the extents 
    else
    {
        debug5 << method
               << "Detected that bounds are not set. Setting them too."
               << std::endl;

        for(size_t i=0; i<variables.size() ; i++)
        {
            std::string varName = variables[i];
            std::string fullVarName;

            std::vector <uint64_t> dims;
            FastQuery::DataType type;

            qp.getVariableInfo( varName, fullVarName, dims, &type );

            debug5 << "HDF-FQ/FastBit get var info for "
                   << varName << " type is "  << type
                   << std::endl;
            
            if (type == FastQuery::FQT_FLOAT)
            {
                float vals[2];
                bool ret = qp.getActualRange(varName, vals);

                begins[i] = vals[0];
                ends[i]   = vals[1];

                debug5 << "HDF-FQ/FastBit calculated bounds for "
                       << varName << (ret ? " passed" : " failed")
                       << " are "  << begins[i] << ", " << ends[i]
                       << std::endl;
            }
            else if (type == FastQuery::FQT_DOUBLE)
            {
                double vals[2];
                bool ret = qp.getActualRange(varName, vals);

                begins[i] = vals[0];
                ends[i]   = vals[1];

                debug5 << "HDF-FQ/FastBit calculated bounds for "
                       << varName << (ret ? " passed" : " failed")
                       << " are "  << begins[i] << ", " << ends[i]
                       << std::endl;
            }
            else
            {
                debug5 << method
                       << "Warning:: data type min/max not yet supported "
                       << "for this " << variables[i]
                       << "'s histogram computation"
                       << std::endl;
                
                EXCEPTION2(NonCompliantFileException,
                           "H5Part ConstructHistogram",
                           "Data type min/max not yet not supported for histogram computation.");
            }
        }
    }

    bool cacheCurrentHisto = false;

    // We need to check the dimensionality of the histogram to
    // determine what FastBit call to make. We'll assume that the
    // bounds contain junk for now and that num_bins is what will
    // drive the histogram creation.
    if (!isCached)
    {
        if (variables.size() == 1)
        {
            double stride = (ends[0] - begins[0]) / double(numBins[0]);

            // make sure that begin + stride*num_bins > end
            stride = ibis::util::incrDouble(stride);

            // Need to have a dummy query
            if( condition.empty() )
            {
              std::string var = variables[0];
              double min = begins[0];
              double max = ends[0];
              // std::string min_cond = DoubleToString(min) + "<=" + var;
              std::string min_cond = var + ">=" + DoubleToString(min);
              std::string max_cond = var + "<=" + DoubleToString(max);

              // FIXME: Check if all variables are availale and have
              // FastBit index.  Possibly issue warning when
              // encountering a variable for which this
              //  is not true.
              condition = min_cond + " && " + max_cond;
            }
              
            debug5 << method
                   << "Asking reader for 1D histogram for"
                   << " timestep=" << timestep << std::endl
                   << " condition=" << (condition.empty() ? "None" : condition)
                   << " variables=" << variables[0]
                   << " begin=" << begins[0]
                   << " end=" << ends[0]
                   << " stride=" << stride
                   << " bins=" << numBins[0]
                   << std::endl;

            std::vector<uint32_t> counts;

            debug5 << method << __LINE__ << std::endl;          

            qp.get1DHistogram(condition.c_str(),
                              variables[0], begins[0], ends[0], stride,
                              counts);
            
            // TODO: Possibly replace in order to avoid copying of the
            // histograms

            debug5 << method << __LINE__ << std::endl;          

            // Copy the counts into the specification
            spec->SetCounts(counts);

            debug5 << method << __LINE__ << std::endl;          

            // Fill in the bounds array.
            std::vector<double> &bounds = spec->GetBounds()[0];
            bounds.resize(numBins[0]+1);
            
            debug5 << method << __LINE__ << std::endl;          

            for (unsigned int i=0; i<=numBins[0]; ++i)
              bounds[i] = begins[0] + double(i) * stride;

            debug5 << method << __LINE__ << std::endl;          

            // Make sure that the bounds specified flag is defined properly
            spec->SetBoundsSpecified();

            debug5 << method << __LINE__ << std::endl;          

            cacheCurrentHisto = true;
        }
        else if (variables.size() == 2)
        {
            std::vector<uint32_t> counts;

            if (regularBinning)
            {
                double stride0 = (ends[0] - begins[0]) / double(numBins[0]);

                // make sure that begin + stride*num_bins > end
                stride0 = ibis::util::incrDouble(stride0);
                
                double stride1 = (ends[1] - begins[1]) / double(numBins[1]);

                // make sure that begin + stride*num_bins > end
                stride1 = ibis::util::incrDouble(stride1);

                // Need to have a dummy query
                if( condition.empty() )
                {
                  std::string var = variables[0];
                  double min = begins[0];
                  double max = ends[0];
                  // std::string min_cond = DoubleToString(min) + "<=" + var;
                  std::string min_cond = var + ">=" + DoubleToString(min);
                  std::string max_cond = var + "<=" + DoubleToString(max);

                  // FIXME: Check if all variables are availale and have
                  // FastBit index.  Possibly issue warning when
                  // encountering a variable for which this
                  //  is not true.
                  condition = min_cond + " && " + max_cond;

                  var = variables[1];
                  min = begins[1];
                  max = ends[1];
                  // std::string min_cond = DoubleToString(min) + "<=" + var;
                  min_cond = var + ">=" + DoubleToString(min);
                  max_cond = var + "<=" + DoubleToString(max);

                  // FIXME: Check if all variables are availale and have
                  // FastBit index.  Possibly issue warning when
                  // encountering a variable for which this
                  //  is not true.
                  condition += " && " + min_cond + " && " + max_cond;
                }
              
                debug5 << method
                       << "Asking reader for 2D histogram for timestep for "
                       << " timestep " << timestep << std::endl
                       << " condition " << (condition.empty() ? "None" : condition)
                       << " variables " << variables[0]
                       << " and " << variables[1] << std::endl;

                qp.get2DHistogram(condition.c_str(),
                                  variables[0], begins[0], ends[0], stride0, 
                                  variables[1], begins[1], ends[1], stride1,
                                  counts);

                // Fill in the bounds array.
                std::vector<double> &bounds0 = spec->GetBounds()[0];
                bounds0.resize(numBins[0]+1);
            
                for (unsigned int i=0; i<=numBins[0]; ++i)
                  bounds0[i] = begins[0] + double(i) * stride0;

                // Fill in the bounds array.
                std::vector<double> &bounds1 = spec->GetBounds()[1];
                bounds1.resize(numBins[1]+1);
            
                for (unsigned int i=0; i<=numBins[1]; ++i)
                  bounds1[i] = begins[1] + double(i) * stride1;
            }
            else
            {
              // NOTE - ARS/JW - deprecate adaptive binning.

                // Adaptive binning
                // if (condition.empty())
                // {
                //     // No condition.
                //     debug5 << method
                //            << "Asking reader for adaptive 2D histogram for"
                //            << " timestep " << timestep << std::endl
                //            << " variables " << variables[0]
                //            << " and " << variables[1] << std::endl;

                //     qp.get2DHistogram(variables[0],
                //                    variables[1],
                //                    numBins[0],
                //                    numBins[1],
                //                    spec->GetBounds()[0],
                //                    spec->GetBounds()[1],
                //                    counts, 0);
                // }
                // else
                {
                    // Conditional histogram.
                    debug5 << method
                           << "Asking reader for adaptive 2D histogram for "
                           << "timestep " << timestep << std::endl
                           << " condition " << condition
                           << " variables " << variables[0]
                           << " and " << variables[1] << std::endl;

                    qp.get2DHistogram(condition.c_str(),
                                      variables[0],
                                      variables[1],
                                      numBins[0],
                                      numBins[1],
                                      spec->GetBounds()[0],
                                      spec->GetBounds()[1],
                                      counts);
                }
            }

            debug5 << "begin0 = " << begins[0] << ", end0 = " << ends[0]
                   << ", begin1 " << begins[1] << ", end1 = " << ends[1]
                   << std::endl
                   << "numbins0 = " << numBins[0]
                   << ", numBins1 = " << numBins[1]
                   << std::endl;

            // Copy the counts into the specification
            spec->SetCounts(counts);

            // Make sure that the bounds specified flag is defined properly
            spec->SetBoundsSpecified();

            cacheCurrentHisto = true;
        }
        else
        {
            // FIXME: generic histogram, will need to take slices.
            debug5 << method << "Need to implement generic histogram slices."
                   << std::endl;
            EXCEPTION1(VisItException, "Generic histogram slices not implemented.");
        }
    }
    
            debug5 << method << __LINE__ << std::endl;          

    //Add the current histogram to the cache
    if (cacheCurrentHisto)
        histoCache.addToCache(spec);
            debug5 << method << __LINE__ << std::endl;          

}
#endif


// ****************************************************************************
//  Method: avtH5PartFileFormat::ConstructIdentifiersFromDataRangeSelection
//
//  Purpose:
//      Constructs the identifiers that fall within a range selection.
//
//  Programmer: Prabhat (Hank added the stub)
//  Creation:   March 6, 2008
//
//  Modifications:
//
//    Gunther H. Weber, Mon Feb 15 20:48:08 PST 2010
//    Ported/moved to new H5Part plugin.
//
//    Gunther H. Weber, Wed Aug  4 18:01:04 PDT 2010
//    Ensure that variables used for constructing identifiers from query are
//    different from global variables for query. (We do not want a call to
//    this method to change the currently active query!)
//
//    Gunther H. Weber, Wed Aug 17 13:15:46 PDT 2011
//    Make queries inclusive to match Threshold operator
//
//    Brad Whitlock, Thu Mar 15 14:30:44 PDT 2012
//    Set the avtIdentifierSelection's name so it will have the right variable
//    later in it when we later store the results as an 
//    avtFloatingPointIdNamedSelection.
//
// ****************************************************************************
#ifdef HAVE_LIBFASTQUERY

avtIdentifierSelection *
avtH5PartFileFormat::ConstructIdentifiersFromDataRangeSelection(
    std::vector<avtDataSelection *> &drs)
{
    // NOTE: This function is only called from GetAuxiliaryData if
    // useFastBitIndex is true. So no need to check that option here.
    int t1 = visitTimer->StartTimer();

    std::string method =
        "avtH5PartFileFormat::ConstructIdentifiersFromDataRangeSelection(): ";

    std::string selectionQuery;
    bool selectionSpecified = false;

    debug5 << method << "Creating query string from selections." << std::endl;

    // FIXME: This function has some similarity/shared code with
    // RegisterDataSelections() that should probably be put in a
    // common helper function.

    // Iterate over all the selections
    idVariableName = defaultIdVariableName;
    
    for (size_t i=0; i<drs.size(); ++i)
    {
        if (std::string(drs[i]->GetType()) ==
            std::string("Data Range Selection"))
        {
            // Data Range is valid. -> Create the query string
            avtDataRangeSelection *dr = (avtDataRangeSelection*)drs[i];

            std::string var = dr->GetVariable();
            double min = dr->GetMin();
            double max = dr->GetMax();
            // std::string min_cond = DoubleToString(min) + "<=" + var;
            std::string min_cond = var + ">=" + DoubleToString(min);
            std::string max_cond = var + "<=" + DoubleToString(max);

            // FIXME: Check if all variables are availale and have
            // FastBit index.  Possibly issue warning when
            // encountering a variable for which this
            //  is not true.
            if (selectionQuery.empty())
                selectionQuery = min_cond + " && " + max_cond;
            else
                selectionQuery += " && " + min_cond + " && " + max_cond;

            selectionSpecified = true;
        }
        else if (std::string(drs[i]->GetType()) ==
                std::string("Identifier Data Selection"))
        {
            // Identifier selection.
            debug5 << method << "Named Selection found."  << std::endl;

            avtIdentifierSelection *id = (avtIdentifierSelection*)drs[i];
            const std::vector<double>& ids = id->GetIdentifiers();

            if (!ids.empty())
            {
                std::string id_string;
                // FIXME: Constructing a string is very inefficient
                // for large selection. However, there seems no
                // HDF5_FastQuery function to combine a query for a
                // passed array of identifiers with a string
                // specifying thresholds. It may make sense to add
                // this functionality or add a non-string based API
                // for range queries and logical combination of
                // queries.
                if(!id->GetIdVariable().empty())
                    idVariableName = id->GetIdVariable();
                debug5 << method << "Setting idVariableName to "
                       << idVariableName << std::endl;

                ConstructIdQueryString(ids, idVariableName, id_string);

                if (selectionQuery.empty())
                    selectionQuery = id_string;
                else
                    selectionQuery += " && " + id_string;

                selectionSpecified = true;
            }
        }
    }

    debug5 << method << "Query string is " << selectionQuery << std::endl;

    std::vector<double> returnIds;

    if (selectionSpecified)
    {
        if (!file)
        {
            file = H5PartOpenFile(filenames[0], H5PART_READ);

            if (!file)
              EXCEPTION1(InvalidFilesException, "Could not open file.");
        }

        H5PartSetStep(file, activeTimeStep);

        FastQuery::QueryProcessor qp =
          FastQuery::QueryProcessor( filenames[0], FastQuery::FQ_H5Part,
                                     filenames[0] );

        qp.setVariablePathPrefix( getVariablePathPrefix(activeTimeStep).c_str() );
        qp.setIndexPathPrefix( getFastBitIndexPathPrefix(activeTimeStep).c_str() );

        std::vector<uint64_t> selectionQueryResults;
        qp.executeQuery(selectionQuery, selectionQueryResults);
        
        debug5 << method << "Execution of query string (length "
               << selectionQuery.size() << ") resulted in "
               << selectionQueryResults.size() << " entries." << std::endl;

        VarNameToInt64Map_t::const_iterator it =
            particleVarNameToTypeMap.find(idVariableName);
        if (it == particleVarNameToTypeMap.end())
            EXCEPTION1(InvalidVariableException, idVariableName);

        H5PartSetView(file, -1, -1);
        h5part_int64_t *indices = new h5part_int64_t[selectionQueryResults.size()];
        std::copy(selectionQueryResults.begin(), selectionQueryResults.end(), indices);
        H5PartSetViewIndices(file, indices, selectionQueryResults.size());
        delete[] indices;

        if (it->second == H5PART_FLOAT64)
        {
            h5part_float64_t *idList = new h5part_float64_t[selectionQueryResults.size()];;
            if (H5PartReadDataFloat64(file, idVariableName.c_str(), idList) != H5PART_SUCCESS)
            {
                delete[] idList;
                EXCEPTION1(InvalidVariableException, idVariableName);
            }
            else
            {
                returnIds.resize(selectionQueryResults.size());
                for (size_t i=0; i<selectionQueryResults.size(); ++i)
                    returnIds[i] = idList[i];
                delete[] idList;
            }
        }
        else if (it->second == H5PART_FLOAT32)
        {
            h5part_float32_t *idList = new h5part_float32_t[selectionQueryResults.size()];
            if (H5PartReadDataFloat32(file, idVariableName.c_str(), idList) != H5PART_SUCCESS)
            {
                delete[] idList;
                EXCEPTION1(InvalidVariableException, idVariableName);
            }
            else
            {
                returnIds.resize(selectionQueryResults.size());
                for (size_t i=0; i<selectionQueryResults.size(); ++i)
                    returnIds[i] = idList[i];
                delete[] idList;
            }
        }
        else if (it->second == H5PART_INT64)
        {
            h5part_int64_t *idList = new h5part_int64_t[selectionQueryResults.size()];
            if (H5PartReadDataInt64(file, idVariableName.c_str(), idList) != H5PART_SUCCESS)
            {
                delete[] idList;
                EXCEPTION1(InvalidVariableException, idVariableName);
            }
            else
            {
                returnIds.resize(selectionQueryResults.size());
                for (size_t i=0; i<selectionQueryResults.size(); ++i)
                    returnIds[i] = idList[i];
                delete[] idList;
            }
        }
        else if (it->second == H5PART_INT32)
        {
            h5part_int32_t *idList = new h5part_int32_t[selectionQueryResults.size()];
            if (H5PartReadDataInt32(file, idVariableName.c_str(), idList) != H5PART_SUCCESS)
            {
                delete[] idList;
                EXCEPTION1(InvalidVariableException, idVariableName);
            }
            else
            {
                returnIds.resize(selectionQueryResults.size());
                for (size_t i=0; i<selectionQueryResults.size(); ++i)
                    returnIds[i] = idList[i];
                delete[] idList;
            }
        }
        else
        {
            EXCEPTION2(NonCompliantFileException,
                       "H5Part ConstructIdentifiersFromDataRangeSelection",
                       std::string("Id variable ") + idVariableName +
                       std::string(" uses an unsupported value type. Supported "
                                   "types are FLOAT64, FLOAT32 and INT64.)"));
        }
        // FIXME: Possibly handle other data types for the id variable
    }

    avtIdentifierSelection *rv = new avtIdentifierSelection();
    rv->SetIdentifiers(returnIds);
    rv->SetIdVariable(idVariableName);

    visitTimer->StopTimer(t1, "H5PartFileFormat::ConstructIdentifiersFromDataRangeSelection()");
    return rv;
}
#endif


// ****************************************************************************
//  Method: avtH5PartFileFormat::ConstructIdQueryString
//
//  Purpose:
//      Generate a string from a list of identifiers
//
//  Programmer: Prabhat
//  Creation:   March 6, 2008
//
//  Modifications:
//
//    Gunther H. Weber, Mon Jul 26 11:46:15 PDT 2010
//    Construct entire string in stream buffer for speed up.
//
//    Gunther H. Weber, Wed Aug 17 13:36:48 PDT 2011
//    Commented out potentially very long debug message.
//
//    Brad Whitlock, Thu Mar 15 14:26:52 PDT 2012
//    Pass in the id variable name.
//
// ****************************************************************************
#ifdef HAVE_LIBFASTQUERY

void avtH5PartFileFormat::ConstructIdQueryString(
    const std::vector<double>& identifiers,
    const std::string &idVar,
    std::string& id_string)
{
    int t1 = visitTimer->StartTimer();

    //time_t startTime = time(0);

    id_string.clear();

    debug5 << "ConstructIdQueryString(): " << idVar
           << " length = " << identifiers.size() << std::endl;

    if (identifiers.size()>0)
    {
        std::ostringstream queryStream;
        queryStream << "" << idVar << " in ( ";

        // FIXME: Precision always sufficient? In particular with
        // floats this may not be true.
        for (size_t j=0; j<identifiers.size()-1; j++)
            queryStream << setprecision(32) << identifiers[j] << ", ";

        queryStream << setprecision(32) << identifiers[identifiers.size()-1]
                    << " )";

        if (queryStream)
            id_string =  queryStream.str();
        else
            EXCEPTION1(VisItException,
                       "Error constructing identifier query string.")
    }

    //debug5 << "id_string is " << id_string << std::endl;

    visitTimer->StopTimer(t1, "H5PartFileFormat::ConstructIdQueryString()");
}
#endif


// ****************************************************************************
//  Method: avtH5PartFileFormat::PerformQuery
//
//  Purpose:
//      Generate a string from a list of identifiers
//
//  Programmer: Prabhat
//  Creation:   March 6, 2008
//
//  Modifications:
//
// ****************************************************************************
#ifdef HAVE_LIBFASTQUERY

void avtH5PartFileFormat::PerformQuery()
{
    // NOTE: Only called from SelectParticlesToRead if useFastBintIndex is true.

    int t1 = visitTimer->StartTimer();

    debug5 << "avtH5PartFileFormat::PerformQuery(): Running query."
           << std::endl;

    FastQuery::QueryProcessor qp =
      FastQuery::QueryProcessor( filenames[0], FastQuery::FQ_H5Part,
                                 filenames[0] );

    qp.setVariablePathPrefix( getVariablePathPrefix(activeTimeStep).c_str() );
    qp.setIndexPathPrefix( getFastBitIndexPathPrefix(activeTimeStep).c_str() );

    queryResults.clear();

    if (querySpecified == stringQuery)
    {
        debug5 << "String query specified: " << queryString << std::endl;

        qp.executeQuery(queryString, queryResults);
        
        dataSelectionActive = true;
    }
    else if (querySpecified == idListQuery)
    {
        debug5 << "Id list (" << idVariableName
               << ") query with " << queryIdList.size()
               << " ids specified." << std::endl;

        if (queryIdList.size() > 0)
        {
          // Note: Ignore the number of hits that was returned
          // if( defaultSortedVariableName.size() )
          // {
            qp.executeEqualitySelectionQuery(idVariableName,
                                             queryIdList,
                                             queryResults);
          // }
          // else
          // {
          //     std::ostringstream queryStream;
            
          //     for( int i=0; i<queryIdList.size(); ++i)
          //     {
          //         queryStream << idVariableName << "=="
          //                  << setprecision(32) << queryIdList[i];
                  
          //         if( i < queryIdList.size()-1)
          //           queryStream << " || ";            
          //     }
              
          //     debug5 << "String query specified: " << queryStream.str()
          //            << std::endl;
              
          //     qp.executeQuery(queryStream.str(), queryResults);
          // }
          
          dataSelectionActive = true;
        }
    }
    else
    {
        debug5 << "No query specified." << std::endl;
    }

    debug5 << "Query resulted in " << queryResults.size() << " hits."
           << std::endl;

    queryResultsValid = true;

    visitTimer->StopTimer(t1, "H5PartFileFormat::PerformQuery()");
}
#endif

// ****************************************************************************
// Method: avtH5PartFileFormat::GetDecomp
//
// Purpose: 
//   Distribute a set of objects evenly amongst all ranks
//
// Arguments:
//   nObjects   : The number of objects to be distributed
//   firstIndex : The index of the first object on this rank
//   lastIndex  : The index of the lastst object on this rank
//   total      : The total number of ojects on this rank
//
// Programmer: Allen Sanderson
// Creation:   Fri Jan 20 2017
//
// ****************************************************************************
void avtH5PartFileFormat::GetDecomp( h5part_int64_t nObjects,
                                     h5part_int64_t &firstIndex,
                                     h5part_int64_t &lastIndex,
                                     h5part_int64_t &total )
{
#ifdef PARALLEL
  h5part_int64_t rank   = PAR_Rank();
  h5part_int64_t nProcs = PAR_Size();

  h5part_int64_t nObjectsPerProc = (nObjects / nProcs);
  h5part_int64_t oneExtraUntil   = (nObjects % nProcs);
    
  if (rank < oneExtraUntil)
  {
    firstIndex = (rank  ) * (nObjectsPerProc+1);
    lastIndex  = (rank+1) * (nObjectsPerProc+1) - 1;

    total = lastIndex - firstIndex + 1;
  }
  else if( nObjectsPerProc ) 
  {
    firstIndex = (rank  ) * (nObjectsPerProc) + oneExtraUntil;
    lastIndex  = (rank+1) * (nObjectsPerProc) + oneExtraUntil - 1;

    total = lastIndex - firstIndex + 1;
  }
  else
  {
    firstIndex = -1;
    lastIndex  = -2;
    total = 0;
  }
#else   
  firstIndex = 0;
  lastIndex = nObjects-1;
  total = nObjects;
#endif
}

// ****************************************************************************
// Method: avtH5PartFileFormat::getVariablePathPrefix
//
// Purpose: 
//   Get the variablePathPrefix
//
// Programmer: Allen Sanderson
// Creation:   31 May 2017
//
// ****************************************************************************
const std::string
avtH5PartFileFormat::getVariablePathPrefix( int timestep )
{
  std::ostringstream prefix;

  prefix << variablePathPrefix << timestep;

  return prefix.str();
}

// ****************************************************************************
// Method: avtH5PartFileFormat::getFastBitIndexPathPrefix
//
// Purpose: 
//   Get the FastBitIndexPathPrefix
//
// Programmer: Allen Sanderson
// Creation:   31 May 2017
//
// ****************************************************************************
#ifdef HAVE_LIBFASTQUERY
const std::string
avtH5PartFileFormat::getFastBitIndexPathPrefix( int timestep )
{
  std::ostringstream prefix;

  prefix << fastBitIndexPathPrefix << "/"
         << getVariablePathPrefix( timestep );

  return prefix.str();
}
#endif
