/*****************************************************************************
*
* Copyright (c) 2000 - 2007, The Regents of the University of California
* Produced at the Lawrence Livermore National Laboratory
* All rights reserved.
*
* This file is part of VisIt. For details, see http://www.llnl.gov/visit/. 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 materials provided with the distribution.
*  - Neither the name of the UC/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  THE  REGENTS  OF  THE  UNIVERSITY OF
* CALIFORNIA, 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.
*
*****************************************************************************/

// ************************************************************************* //
//                            avtCarpetHDF5FileFormat.C                           //
// ************************************************************************* //

#include <avtCarpetHDF5FileFormat.h>

#include <string>
#include <map>
#include <set>
#include <cassert>
#include <cstdlib>
#include <ctime>
#include <cstring>

// POSIX specific!
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>

#include <vtkIdList.h>
#include <vtkIntArray.h>
#include <vtkFloatArray.h>
#include <vtkRectilinearGrid.h>
#include <vtkStructuredGrid.h>
#include <vtkUnstructuredGrid.h>
#include <vtkPoints.h>
#include <avtGhostData.h>
#include <vtkCellData.h>
#include <vtkUnsignedCharArray.h>
#include <DebugStream.h>

#include <avtDatabaseMetaData.h>
#include <avtStructuredDomainNesting.h>
#include <avtVariableCache.h>
#include <avtStructuredDomainBoundaries.h>
#include <avtIntervalTree.h>

#include <Expression.h>

#include <InvalidVariableException.h>

#include "H5Index.h"
using namespace H5Index;

using namespace std;



// ****************************************************************************
//  Method: avtCarpetHDF5FileFormat constructor
//
//  Programmer: Christian Reisswig -- generated by xml2avt
//  Creation:   Tue Jan 22 15:54:40 PST 2008
//
// ****************************************************************************

std::map<string, avtCarpetHDF5FileFormat::multi_file*> avtCarpetHDF5FileFormat::file_cache;
avtCarpetHDF5FileFormat::avtCarpetHDF5FileFormat(const char *filename)
  : avtMTMDFileFormat(filename), data_file(0), xcoord_file(0), ycoord_file(0), zcoord_file(0)
{
    // INITIALIZE DATA MEMBERS
    open_all_files(filename);
    
}


// ****************************************************************************
//  Method: avtEMSTDFileFormat::GetNTimesteps
//
//  Purpose:
//      Tells the rest of the code how many timesteps there are in this file.
//
//  Programmer: Christian Reisswig -- generated by xml2avt
//  Creation:   Tue Jan 22 15:54:40 PST 2008
//
// ****************************************************************************

int
avtCarpetHDF5FileFormat::GetNTimesteps(void)
{
    return data_file->n_timesteps;
}


// ****************************************************************************
//  Method: avtCarpetHDF5FileFormat::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: Christian Reisswig -- generated by xml2avt
//  Creation:   Tue Jan 22 15:54:40 PST 2008
//
// ****************************************************************************

void
avtCarpetHDF5FileFormat::FreeUpResources(void)
{
   for (int i=0; i < data_file->timesteps.size(); i++)
   {
      data_file->timesteps[i].unset_ranges();
   }
}


// ****************************************************************************
//  Method: avtCarpetHDF5FileFormat::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: Christian Reisswig -- generated by xml2avt
//  Creation:   Tue Jan 22 15:54:40 PST 2008
//
// ****************************************************************************

void
avtCarpetHDF5FileFormat::PopulateDatabaseMetaData(avtDatabaseMetaData *md, int timeState)
{
   md->SetTimes(data_file->times);
   md->SetTimesAreAccurate(true);
   md->SetCycles(data_file->cycles);
   md->SetCyclesAreAccurate(true);

    // Prevent VisIt from sorting the variables.
   md->SetMustAlphabetizeVariables(false);

   md->SetNumStates(data_file->n_timesteps);
   md->SetMustRepopulateOnStateChange(true);
   
   /*cout << "----" << endl;
   cout << data_file->n_timesteps <<  ", "
        << data_file->cycles.size() << ", "
        << data_file->times.size() << endl;
   
   
   for (int i=0; i < data_file->cycles.size(); ++i) {
      cout << data_file->cycles[i] << endl;
      cout << data_file->times[i] << endl;
   }*/
   
   
   // set up Cartesian grids (if any) and Multi-Patch grids (if any)
   for (int i=0; i < 2; ++i)
   {
      bool isCartesian = data_file->timesteps[timeState].cart_comp[0].size() > 0;
      if (!isCartesian && i == 0) continue;
      if (i == 1 && data_file->timesteps[timeState].multi_comp[0].size() > 0) isCartesian = false;
      if (isCartesian && i == 1) continue;
      const string meshname = isCartesian ? "Carpet AMR-grid" : "Carpet Multipatch";
      const avtMeshType meshType = isCartesian ? AVT_AMR_MESH : AVT_CURVILINEAR_MESH;
      const int nblocks = get_nblocks(timeState, isCartesian);

      vector<vector<dataset_entry> > &comp = isCartesian ?
        data_file->timesteps[timeState].cart_comp :
        data_file->timesteps[timeState].multi_comp;
      // find the maximum number of dimensions over all components
      int maxndims = -1;
      for (int i=0; i < comp.size(); i++)
        if (maxndims < comp[i][0].ndims()) maxndims = comp[i][0].ndims();

      // Here's the call that tells the meta-data object that we have a mesh:
      AddMeshToMetaData(md, meshname, meshType, NULL, nblocks, 0, maxndims, maxndims);
          
      avtMeshMetaData *mesh = (avtMeshMetaData *) md->GetMesh(meshname);
      mesh->xUnits = mesh->yUnits = mesh->zUnits = "M";
      mesh->hasSpatialExtents = false;
      mesh->blockTitle = "comps";
      mesh->blockPieceName = "comp";
      mesh->numGroups = get_max_reflevels(timeState, isCartesian);
      mesh->groupTitle = "levels";
      mesh->groupPieceName = "level";
      vector<int> groupIds(nblocks);
      vector<string> blockPieceNames(nblocks);
      for (int i = 0 ; i < nblocks; i++)
      {
          char tmpName[128];
          int level = comp[0][i].rl();
          int local_patch = i;
          groupIds[i] = level; 
          sprintf(tmpName, "ref-level=%d, component=%d", level, local_patch);
          blockPieceNames[i] = tmpName;
      }
      mesh->blockNames = blockPieceNames;
      md->AddGroupInformation(mesh->numGroups, nblocks, groupIds);
          
      for (int i=0; i < comp.size(); i++) {
          string str = comp[i][0].varname() + (isCartesian ? "" : "(MP)");
          if (comp[i][0].comps() == 1) {
            // add scalar variable
            avtScalarMetaData *smd = new avtScalarMetaData(str.replace(str.find("::"), 2, "--"), meshname, AVT_NODECENT);
            smd->centering = AVT_NODECENT;
            md->Add(smd);
          } else if (comp[i][0].comps() > 1) {
            // add vector variables
            int dim = comp[i][0].comps();
            avtVectorMetaData *vmd = new avtVectorMetaData(str.replace(str.find("::"), 2, "--"), meshname, AVT_NODECENT, dim);
            vmd->centering = AVT_NODECENT;
            md->Add(vmd);
          }
      }
   
   }
}


// ****************************************************************************
//  Method: avtCarpetHDF5::BuildDomainAuxiliaryInfo
//
//  Purpose: Build the two data structures needed to support nesting and
//  abutting of AMR subgrids.
//  This is restricted to the Cartesian AMR grid and does not apply to the 
//  curvilinear multiblock grid.
//
//  Note: These are *never* explicitly served up to VisIt like a GetMesh or
//  GetVar call would do. Instead, we essentially 'publish' them to VisIt
//  by sticking the structures we create here into the database cache. VisIt
//  will try to look for them there when it needs them.
//
//  Programmer: Christian Reisswig
//  Creation:   Sun Feb 28 14:31:17 PDT 2010
//
// ****************************************************************************

void
avtCarpetHDF5FileFormat::BuildDomainAuxiliaryInfo(int timeState)
{
#ifdef MDSERVER
    return;
#endif

    const vector<vector<dataset_entry> > &comp = data_file->timesteps[timeState].cart_comp;
    vector<int> dims(3,1);

    int num_dims = dim;   // defined in vtCarpetHDF5FileFormat.h as 'static const int dim = 3;'
    int num_levels = data_file->timesteps[timeState].max_cart_rl;
    int num_patches = data_file->timesteps[timeState].cart_comp[0].size();
    int visit_2d_granularity_fudge = (getenv("VISIT_2D_GRANULARITY_FUDGE") != NULL) &&
                                     (strcmp(getenv("VISIT_2D_GRANULARITY_FUDGE"), "yes") == 0);

    // first, look to see if we don't already have it cached
    void_ref_ptr vrTmp = cache->GetVoidRef("any_mesh",
                                   AUXILIARY_DATA_DOMAIN_NESTING_INFORMATION,
                                  timeState, -1);
    if (*vrTmp == NULL)
    {
        //
        // build the avtDomainNesting object
        //
        avtStructuredDomainNesting *dn =
            new avtStructuredDomainNesting(num_patches, num_levels);

        dn->SetNumDimensions(num_dims);

        //
        // Set refinement level ratio information
        // We hardcode this to 2:1 ratio for each direction!
        //
        vector<int> ratios(3,1);
        dn->SetLevelRefinementRatios(0, ratios);
        for (int i = 1; i < num_levels; i++)
        {
           ratios[0] = 2;
           ratios[1] = 2;
           ratios[2] = 2;
           dn->SetLevelRefinementRatios(i, ratios);
        }
        
        //
        // set each domain's level, children and logical extents
        //
        for (int i = 0; i < num_patches; i++)
        {
            vector<int> childPatches;
            float x0 = comp[0][i].xmin();
            float x1 = comp[0][i].xmax();
            float y0 = comp[0][i].ymin();
            float y1 = comp[0][i].ymax();
            float z0 = comp[0][i].zmin();
            float z1 = comp[0][i].zmax();

            bool ghost = true;
                dims = comp[0][i].npoints(ghost);
            
            for (int j = 0; j < num_patches; j++)
            {
                if (comp[0][j].rl() != comp[0][i].rl()+1)
                    continue;

                float a0 = comp[0][j].xmin();
                float a1 = comp[0][j].xmax();
                float b0 = comp[0][j].ymin();
                float b1 = comp[0][j].ymax();
                float c0 = comp[0][j].zmin();
                float c1 = comp[0][j].zmax();

                float dx = comp[0][j].delta()[0]/2;
                float dy = comp[0][j].delta()[1]/2;
                float dz = comp[0][j].delta()[2]/2;

                // compare with a half-finer level offset to work around floating point roundoff issues
                if (((dims[0] > 1 || visit_2d_granularity_fudge) && (a0 >= x1-dx || x0 >= a1-dx)) ||
                    ((dims[1] > 1 || visit_2d_granularity_fudge) && (b0 >= y1-dy || y0 >= b1-dy)) ||
                    ((dims[2] > 1 || visit_2d_granularity_fudge) && (c0 >= z1-dz || z0 >= c1-dz)))
                    continue;

                childPatches.push_back(j);
            }

            //bool isOutermostX = comp[0][i].is_at_upper_total_rl_domain_boundary(0, data_file->timesteps[timeState].totalCartDomainIExt[comp[0][i].rl()][1]);
                //bool isOutermostY = comp[0][i].is_at_upper_total_rl_domain_boundary(1, data_file->timesteps[timeState].totalCartDomainIExt[comp[0][i].rl()][3]);
                //bool isOutermostZ = comp[0][i].is_at_upper_total_rl_domain_boundary(2, data_file->timesteps[timeState].totalCartDomainIExt[comp[0][i].rl()][5]);

            // fill logical extends. Note that these are CELL-indices!
            vector<int> logExts(6);
            logExts[0] = comp[0][i].iorigin(ghost)[0];
            logExts[1] = comp[0][i].iorigin(ghost)[1];
            logExts[2] = comp[0][i].iorigin(ghost)[2];
            logExts[3] = logExts[0] + dims[0] -2;// + int(!isOutermostX); // we have here -2 because logExts is the cell-index of the last cell.
            logExts[4] = logExts[1] + dims[1] -2;// + int(!isOutermostY);
            //logExts[5] = num_dims == 3 ? logExts[2] + dims[2] -2 /*+ int(!isOutermostZ)*/ : 0;
            logExts[5] = logExts[2] + dims[2] -2 >= 0 ? logExts[2] + dims[2] -2 /*+ int(!isOutermostZ)*/ : 0;

            //debug1 << "logExts[0] = " << logExts[0] << ", logExts[1] = " << logExts[1] << ", logExt[2] = " << logExts[2] << endl;
            //debug1 << "logExts[3] = " << logExts[3] << ", logExts[4] = " << logExts[4] << ", logExt[5] = " << logExts[5] << endl;

            dn->SetNestingForDomain(i, comp[0][i].rl() /*-1*/,
                childPatches, logExts);
        }

        void_ref_ptr vr = void_ref_ptr(dn, avtStructuredDomainNesting::Destruct);

        cache->CacheVoidRef("any_mesh", AUXILIARY_DATA_DOMAIN_NESTING_INFORMATION,
            timeState, -1, vr);
    }
}

// ****************************************************************************
//  Method: avtCarpetHDF5FileFormat::GetMesh
//
//  Purpose:
//      Gets the mesh associated with this file.  The mesh is returned as a
//      derived type of vtkDataSet (ie vtkRectilinearGrid, vtkStructuredGrid,
//      vtkUnstructuredGrid, etc).
//
//  Arguments:
//      timestate   The index of the timestate.  If GetNTimesteps returned
//                  'N' time steps, this is guaranteed to be between 0 and N-1.
//      domain      The index of the domain.  If there are NDomains, this
//                  value is guaranteed to be between 0 and NDomains-1,
//                  regardless of block origin.
//      meshname    The name of the mesh of interest.  This can be ignored if
//                  there is only one mesh.
//
//  Programmer: Christian Reisswig -- generated by xml2avt
//  Creation:   Tue Jan 22 15:54:40 PST 2008
//
// ****************************************************************************

vtkDataSet *
avtCarpetHDF5FileFormat::GetMesh(int timestate, int domain, const char *meshname)
{
   vector<int> dims(3,1);
   
   // Read the ndims and number of X,Y,Z nodes from file.
   int ndims = dim;   // defined in vtCarpetHDF5FileFormat.h as 'static const int dim = 3;'

   if (strcmp(meshname, "Carpet Multipatch") == 0) {
      
      float *xarray = 0;
      float *yarray = 0;
      float *zarray = 0;
      
      
      // check if the iteration is present in coordinate files.
      // If not, we use the last stored iteration.
      int actual_timestate = 0;
      
      vector<vector<dataset_entry> > &comp = data_file->timesteps[timestate].multi_comp;
      
      for (int i=0; i < xcoord_file->timesteps.size(); ++i)
         if (xcoord_file->timesteps[i].multi_comp[0][0].cycle() <= comp[0][0].cycle())
            actual_timestate = i;
      
      xcoord_file->get_data(false, actual_timestate, domain, 0, &xarray);
      if(ndims > 1) {
         ycoord_file->get_data(false, actual_timestate, domain, 0, &yarray);
      }
      
      if(ndims > 2) {
         zcoord_file->get_data(false, actual_timestate, domain, 0, &zarray);
      }
      
      dims = xcoord_file->timesteps[actual_timestate].multi_comp[0][domain].npoints();
      const int nnodes = dims[0]*dims[1]*dims[2];
      
      //
      // Create the vtkStructuredGrid and vtkPoints objects.
      //
      vtkStructuredGrid *sgrid = vtkStructuredGrid::New();
      vtkPoints          *points = vtkPoints::New();
      sgrid->SetPoints(points);
      sgrid->SetDimensions(&dims[0]);
      points->Delete();
      points->SetNumberOfPoints(nnodes);
      //
      // Copy the coordinate values into the vtkPoints object.
      //
      float *pts = (float *) points->GetVoidPointer(0);
      float *xc = xarray;
      float *yc = yarray;
      float *zc = zarray;
      if(ndims == 3)
      {
         for(int k = 0; k < dims[2]; ++k)
         for(int j = 0; j < dims[1]; ++j)
         for(int i = 0; i < dims[0]; ++i)
         {
            *pts++ = *xc++;
            *pts++ = *yc++;
            *pts++ = *zc++;
         }
      }
      else if(ndims == 2)
      {
         for(int j = 0; j < dims[1]; ++j)
         for(int i = 0; i < dims[0]; ++i)
         {
            *pts++ = *xc++;
            *pts++ = *yc++;
            *pts++ = 0.;
         }
      }
      
      // Delete temporary arrays.
      delete [] xarray;
      delete [] yarray;
      delete [] zarray;
      return sgrid;
   }
   else
   {
      vtkFloatArray* coords[3] = {0, 0, 0};
      
      vector<vector<dataset_entry> > &comp = data_file->timesteps[timestate].cart_comp;
      
      // Feed AMR domain abutting to VisIt's cache
      BuildDomainAuxiliaryInfo(timestate);
      
      dims = comp[0][domain].npoints();
      
      for (int dim = 0; dim < 3; dim++) {
        coords[dim] = vtkFloatArray::New();
        if (dim < ndims) {
          coords[dim]->SetNumberOfTuples(dims[dim]);
          float *array = (float *)coords[dim]->GetVoidPointer(0);
          for (int j=0; j < dims[dim]; j++) {
            array[j] = comp[0][domain].origin()[dim] +  j * comp[0][domain].delta()[dim];
          }
        } else {
          coords[dim]->SetNumberOfTuples(1);
          coords[dim]->SetComponent(0, 0, 0.);
        }
      }
      
      //
      // Create the vtkRectilinearGrid object and set its dimensions
      // and coordinates.
      //
      vtkRectilinearGrid *rgrid = vtkRectilinearGrid::New();
      rgrid->SetDimensions(&dims[0]);
      rgrid->SetXCoordinates(coords[0]);
      coords[0]->Delete();
      rgrid->SetYCoordinates(coords[1]);
      coords[1]->Delete();
      rgrid->SetZCoordinates(coords[2]);
      coords[2]->Delete();
      
      // The following is not needed because we cut away ghotszones during reading. 
      // When setting up the following ghostzones, I had problems with the nested domains and ghostzones
      // created by Visit
#if 0 
      // Now that you have your mesh, figure out which cells need
      // to be removed.
      bool isOutermostX = comp[0][domain].is_at_upper_total_rl_domain_boundary(0, data_file->timesteps[timestate].totalCartDomainIExt[comp[0][domain].rl()][1]);
      bool isOutermostY = comp[0][domain].is_at_upper_total_rl_domain_boundary(1, data_file->timesteps[timestate].totalCartDomainIExt[comp[0][domain].rl()][3]);
      bool isOutermostZ = comp[0][domain].is_at_upper_total_rl_domain_boundary(2, data_file->timesteps[timestate].totalCartDomainIExt[comp[0][domain].rl()][5]);
      
      // Initialize as all nodes being ghost nodes
      vector<bool> ghostPoints(dims[0]*dims[1]*dims[2], true);
      
      vector<int> first(3);
      first[0] = comp[0][domain].ghosts(true)[0];
      first[1] = comp[0][domain].ghosts(true)[2];
      first[2] = comp[0][domain].ghosts(true)[4];
      vector<int> last(3);
      last[0] = comp[0][domain].npoints(true)[0]+comp[0][domain].ghosts(true)[0]+int(!isOutermostX);
      last[1] = comp[0][domain].npoints(true)[1]+comp[0][domain].ghosts(true)[2]+int(!isOutermostY);
      last[2] = comp[0][domain].npoints(true)[2]+comp[0][domain].ghosts(true)[4]+int(!isOutermostZ);
      
      // Figure out nominal grid nodes.
      for (int k=first[2]; k < last[2]; ++k)
         for (int j=first[1]; j < last[1]; ++j)
            for (int i=first[0]; i < last[0]; ++i)
            {
               const int ijk = k*dims[1]*dims[0] + j*dims[0] + i;
               ghostPoints[ijk] = false;
            }
      
      //
      //  okay, now we have ghost points, but what we really want
      //  are ghost cells ... convert:  if all points associated with
      //  cell are 'real' then so is the cell.
      //

      unsigned char realVal = 0;
      unsigned char ghostVal = 1;
      avtGhostData::AddGhostZoneType(ghostVal,
                                       DUPLICATED_ZONE_INTERNAL_TO_PROBLEM);
      const int ncells = rgrid->GetNumberOfCells();
      vtkIdList *ptIds = vtkIdList::New();
      vtkUnsignedCharArray *ghostCells = vtkUnsignedCharArray::New();
      ghostCells->SetName("avtGhostZones");
      ghostCells->Allocate(ncells);

      for (int i = 0; i < ncells; i++)
      {
         rgrid->GetCellPoints(i, ptIds);
         bool ghost = false;
         for (int idx = 0; idx < ptIds->GetNumberOfIds(); idx++)
              ghost |= ghostPoints[ptIds->GetId(idx)];

         if (ghost)
            ghostCells->InsertNextValue(ghostVal);
         else
            ghostCells->InsertNextValue(realVal);

      }
      
      rgrid->GetCellData()->AddArray(ghostCells);
      ghostCells->Delete();
      ptIds->Delete();

      vtkIntArray *realDims = vtkIntArray::New();
      realDims->SetName("avtRealDims");
      realDims->SetNumberOfValues(6);
      realDims->SetValue(0, first[0]);
      realDims->SetValue(1, last[0]-1);
      realDims->SetValue(2, first[1]);
      realDims->SetValue(3, last[1]-1);
      realDims->SetValue(4, first[2]);
      realDims->SetValue(5, last[2]-1);
      rgrid->GetFieldData()->AddArray(realDims);
      rgrid->GetFieldData()->CopyFieldOn("avtRealDims");
      realDims->Delete();

      rgrid->SetUpdateGhostLevel(0);
#endif
      
      return rgrid;
   }
   
}


// ****************************************************************************
//  Method: avtCarpetHDF5FileFormat::GetVar
//
//  Purpose:
//      Gets a scalar variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      timestate  The index of the timestate.  If GetNTimesteps returned
//                 'N' time steps, this is guaranteed to be between 0 and N-1.
//      domain     The index of the domain.  If there are NDomains, this
//                 value is guaranteed to be between 0 and NDomains-1,
//                 regardless of block origin.
//      varname    The name of the variable requested.
//
//  Programmer: Christian Reisswig -- generated by xml2avt
//  Creation:   Tue Jan 22 15:54:40 PST 2008
//
// ****************************************************************************

vtkDataArray *
avtCarpetHDF5FileFormat::GetVar(int timestate, int domain, const char *varname)
{
   
#ifdef USE_SELECTIONS
   ProcessDataSelections();
#endif

   

   string str = varname;
   const bool isMultiPatch = str.find("(MP)") != string::npos;
   vector<vector<dataset_entry> > &comp = isMultiPatch ? data_file->timesteps[timestate].multi_comp : data_file->timesteps[timestate].cart_comp;

   if (isMultiPatch) {
     str = str.erase(str.find_last_of("(MP)"));
   }
   str.replace(str.find("--"), 2, "::");
       
   int var_no = 0;
   for (int i=0; i < comp.size(); i++) {
     if (comp[i][0].varname() == str) {
       var_no = i;
     }
   }
   
   int dims[] = {1, 1, 1};
   float* data;
   data_file->get_data(not isMultiPatch, timestate, domain, var_no, &data);
   
   int ntuples = comp[0][domain].dims()[0]*comp[0][domain].dims()[1]*comp[0][domain].dims()[2];
   
   vtkFloatArray *rv = vtkFloatArray::New();
   rv->SetNumberOfTuples(ntuples);
   for (int i = 0 ; i < ntuples ; i++) {
      rv->SetTuple1(i, data[i]);  // you must determine value for ith entry.
   }
   
   delete [] data;
   
   return rv;
}


// ****************************************************************************
//  Method: avtCarpetHDF5FileFormat::GetVectorVar
//
//  Purpose:
//      Gets a vector variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      timestate  The index of the timestate.  If GetNTimesteps returned
//                 'N' time steps, this is guaranteed to be between 0 and N-1.
//      domain     The index of the domain.  If there are NDomains, this
//                 value is guaranteed to be between 0 and NDomains-1,
//                 regardless of block origin.
//      varname    The name of the variable requested.
//
//  Programmer: Christian Reisswig -- generated by xml2avt
//  Creation:   Tue Jan 22 15:54:40 PST 2008
//
// ****************************************************************************

vtkDataArray *
avtCarpetHDF5FileFormat::GetVectorVar(int timestate, int domain,const char *varname)
{
   string str = varname; 
   
   const size_t pos = str.find_last_of("(MP)");
   const bool isMultipatch = pos != string::npos;
   if (isMultipatch) str = str.erase(pos);
   vector<vector<dataset_entry> > &comp = isMultipatch ? data_file->timesteps[timestate].multi_comp : data_file->timesteps[timestate].cart_comp;

   str.replace(str.find("--"), 2, "::");

   int var_no = 0;
   for (int i=0; i < comp.size(); i++) {
      if (comp[i][0].varname() == str)
         var_no = i;
   }
   const int ncomponents = comp[var_no][0].comps();
   const int type = comp[var_no][0].type();
   
   float** data = NULL;
   int dims[3] = {1, 1, 1};
   data_file->get_vector_data(not isMultipatch, timestate, domain, var_no, &data);
   
   const int ntuples = dims[0]*dims[1]*dims[2];  // this is the number of entries in the variable.
   
   vtkFloatArray *rv = vtkFloatArray::New();
   int ucomps = (ncomponents == 2 ? 3 : ncomponents);
   rv->SetNumberOfComponents(ucomps);
   rv->SetNumberOfTuples(ntuples);
   float* one_entry = new float[ucomps];
   for (int i = 0 ; i < ntuples ; i++)
   {
      if (type == 1)
      {
         for (int j = 0 ; j < ncomponents ; j++)
            one_entry[j] = data[i][j];
         for (int j = ncomponents ; j < ucomps ; j++)
            one_entry[j] = 0;
      }         
      if (type == 0)
      { 
         for (int j = 0 ; j < ncomponents ; j++)
            one_entry[j] = data[j][i];
         for (int j = ncomponents ; j < ucomps ; j++)
            one_entry[j] = 0;
      }
      rv->SetTuple(i, one_entry);
   }

   delete [] one_entry;
   if (type == 0)
   {
      for (int i=0; i < ncomponents; ++i)
         delete [] data[i];
   }
   if (type == 1)
   {
      for (int i=0; i < ntuples; ++i)
         delete [] data[i];
   }
   delete [] data;
   
   return rv;

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





void avtCarpetHDF5FileFormat::open_all_files(const char* fname)
{
   data_file = open_cached_file(fname);
   
   bool curvilinear = false;
   
   // determine file prefix
   char cslash;
#ifdef WINDOWS
   cslash = '\\';
#else
   cslash = '/';
#endif
   size_t pos = string(fname).find_last_of(cslash);
   string prefix = string(fname).substr(0,pos+1);
   
   
   for (int i=0; i < data_file->timesteps.size(); i++)
   {
      if (data_file->timesteps[i].multi_comp[0].size() > 0)
        curvilinear = true;
   }
   
   
   if (curvilinear) // only look for coordinate-files if we have curvilinear meshes
   {
     if (data_file->files.size() > 1)
     {
       xcoord_file = open_cached_file(prefix + "x.file_0.h5");
       ycoord_file = open_cached_file(prefix + "y.file_0.h5");
       zcoord_file = open_cached_file(prefix + "z.file_0.h5");
     }
     else
     {
       xcoord_file = open_cached_file(prefix + "x.h5");
       ycoord_file = open_cached_file(prefix + "y.h5");
       zcoord_file = open_cached_file(prefix + "z.h5");
     }
   }
}

avtCarpetHDF5FileFormat::multi_file* avtCarpetHDF5FileFormat::open_cached_file(const string fname)
{
  if(getenv("CARPETHDF5_CACHE_METADATA") && strcasecmp(getenv("CARPETHDF5_CACHE_METADATA"), "no") == 0)
    return new multi_file(fname.c_str());

  // if file already exists in cache, check if it is not outdated
  if(file_cache.count(fname))
  {
    multi_file* file = file_cache[fname];

    for(int i=0; i < file->files.size(); i++)
    {
      if(file->files[i]->file_changed_on_disk())
      {
        // remove outdated file from cache. Memory will only be freeed when the file is closed by VisIt.
        file_cache.erase(fname);
        close_cached_file(file);
        break;
      }
    }
  }

  // open file and add to cache
  if(file_cache.count(fname) == 0)
  {
    multi_file* file = new multi_file(fname.c_str());
    file->refcount += 1; // one reference for the cache

    file_cache[fname] = file;
  }
  else
  {
    file_cache[fname]->refcount++;
  }

  return file_cache[fname];
}

void avtCarpetHDF5FileFormat::close_cached_file(avtCarpetHDF5FileFormat::multi_file* file)
{
  // this logic has two maybe unexpected side effects:
  // a file is only ever closed when it is replaced by a new versin of it
  // the close button in visit has no effect whatsoever
  if(file)
  {
    if(--file->refcount == 0)
       delete file;
    else
       file->close();
  }
}



avtCarpetHDF5FileFormat::multi_file::multi_file(const char* fname) : refcount(1)
{
   hid_t H5f = H5Fopen(haveIndex(fname) ? indexFilename(fname).c_str() : fname, H5F_ACC_RDONLY, H5P_DEFAULT);
   
   // access global parameters and attributes group to get the number of io-processes
   hid_t group = H5Gopen (H5f, "Parameters and Global Attributes", H5P_DEFAULT);
   
   hid_t attribute = H5Aopen(group, "nioprocs", H5P_DEFAULT);
   int n_io_procs = 0;
   hid_t status = H5Aread(attribute, H5T_NATIVE_INT, &n_io_procs);
   status = H5Aclose(attribute);

   H5Gclose(group);
   H5Fclose(H5f);
   
   files.resize(n_io_procs);
   
   // open all data-files
   string str(fname);
   size_t pos = str.find_last_of(".");
   str = str.erase(pos);
   pos = str.find_last_of(".");
   if (pos != string::npos)
      str = str.erase(pos);
   
   for (int i=0; i < n_io_procs; i++)
   {
      if (n_io_procs > 1)
      {
         char myname[1000]; 
         sprintf(myname, "%s.file_%d.h5", str.c_str(), i);
         files[i] = new file_t(myname);
      }
      else
      {
         files[i] = new file_t(fname);
      }
      
      // obtain all valid dataset names found in the current file
      // and add them to the global datasetname list
      for (int j=0; j < files[i]->dsetnames.size(); j++)
      {
         files[i]->dsetnames[j].set_file_id(i);
      }
   }
      
   // get number of times and cycles by browsing all available datasets of ALL files
   std::set<int> cycle_seen;
   for (int f=0; f < files.size(); f++)
   {
      vector<dataset_entry>& dsetnames = files[f]->dsetnames;
      for (int i=0; i < dsetnames.size(); i++)
      {
         bool found = cycle_seen.count(dsetnames[i].cycle());

         if (!found)
         {
            cycles.push_back(dsetnames[i].cycle());
            times.push_back(dsetnames[i].time());
            cycle_seen.insert(dsetnames[i].cycle());
         }
      }
   }
   n_timesteps = cycles.size();
   
   
   sort(cycles.begin(), cycles.end());
   sort(times.begin(), times.end());
   
   // get varnames
   std::set<string> varname_seen;
   for (int f=0; f < files.size(); f++)
   {
      vector<dataset_entry> &dsetnames = files[f]->dsetnames;
      for (int i=0; i < dsetnames.size(); i++)
      {
         // replacable by varnames.assign(varname_seen.begin(), varname_end())
         bool found = varname_seen.count(dsetnames[i].varname());
         if (!found)
         {
            varnames.push_back(dsetnames[i].varname());
            //varnames.back().replace(varnames.back().find("::"), 2, "--");
            varname_seen.insert(dsetnames[i].varname());
         }
      }
   }
   
   
   // set the varnames in each file
   for (int i=0; i < n_io_procs; i++)
   {
      files[i]->set_varnames(varnames, cycles);
   }
   
   // now that all files are opened, we need to collect the
   // global meta-data (dataset-attributes etc.) and sort them in an
   // appropriate way
   timesteps.resize(n_timesteps);
   
   // we could probably also replace these with some binary searches, but the
   // maps are simpler to use and fulfill the same purpose
   map<string,int> varmap;
   map<int,int> cyclemap;
   // for very large files, using dataset_entry* might be required
   vector<vector<vector<const dataset_entry *> > > dsets_per_timestep_per_variable(timesteps.size());
   for (int j=0; j < varnames.size(); j++)
   {
     varmap[varnames[j]] = j;
   }
   for (int it=0; it < cycles.size(); it++)
   {
     cyclemap[cycles[it]] = it;
     dsets_per_timestep_per_variable[it].resize(varnames.size());
   }
   
   for (int i=0; i < files.size(); i++)
   {
     for (int c=0; c < files[i]->dsetnames.size(); c++)
     {
       const int it = cyclemap[files[i]->dsetnames[c].cycle()];
       const int j = varmap[files[i]->dsetnames[c].varname()];
       dsets_per_timestep_per_variable[it][j].push_back(&(files[i]->dsetnames[c]));
     }
   }
   for (int it=0; it < timesteps.size(); it++)
   {
      timesteps[it] = dsets_per_timestep_per_variable[it];
   }
}

/*class DisableAutoPrint
{
public:
  DisableAutoPrint()
  {
    Exception::getAutoPrint(func, &client_data);
    Exception::dontPrint();
  }

  ~DisableAutoPrint()
  {
    Exception::setAutoPrint(func, client_data);
  }

private:
        H5E_auto2_t func;
        void *client_data;
};*/

static herr_t H5iter(hid_t group_id, const char *member_name, void *operator_data)
{
   vector<dataset_entry>* dsetnames = (vector<dataset_entry>*) operator_data;
   char   rootname[1000];
   char   fullname[1000];
   H5G_stat_t object_info;
   int dim = avtCarpetHDF5FileFormat::dim;
   
   // build the full name for the current object to process
   H5Iget_name(group_id, rootname, 256);
   
   // the root name is "/", so does not need an extra "/"
   sprintf(fullname, "%s%s%s", rootname, rootname[strlen(rootname)-1]=='/' ? "" : "/",member_name);
   
   // we are interested in datasets only - skip anything else
   H5Gget_objinfo (group_id, member_name, 0, &object_info);
   if (object_info.type != H5G_DATASET)
   {
      if (object_info.type == H5G_GROUP)
      {
         // iterate over all datasets in this group (if it isn't file metadata)
         if (strcmp (member_name, "Parameters and Global Attributes"))
         {
            H5Giterate (group_id, member_name, NULL, H5iter, operator_data);
         }
      }
      
   
      return 0;
   }
   else
   {
      hid_t dataset = H5Dopen(group_id, member_name, H5P_DEFAULT);
      
      hid_t dataspace = H5Dget_space(dataset);
      hid_t attrib;
      
      hsize_t dims[3] = {1,1,1};
      int ndims  = H5Sget_simple_extent_dims(dataspace, dims, NULL);

      // First try to read the shape attribute which is written to
      // index files.  If this is not available, use the shape given
      // in the dataspace.
      H5E_BEGIN_TRY {
        hid_t attrib = H5Aopen(dataset, "h5shape", H5P_DEFAULT);
         if(attrib >= 0) {
            H5Aread(attrib, H5T_NATIVE_HSIZE, dims);
            H5Aclose(attrib);
         }
      } H5E_END_TRY;

      // shift "dims" around because we always treat all datasets as being 3d
      // and HDF5 is in C-order!
      if (ndims == 2)
      {
         dims[2] = dims[1];
         dims[1] = dims[0];
         dims[0] = 1;
      }
      if (ndims == 1)
      {
         dims[2] = dims[0];
         dims[1] = 1;
         dims[0] = 1;
      }
            
      int cycle = 0;
      int rl = 0;
      int map = 0;
      int tl, c; // dummy targets for scanf to get conversion number right
      const int len = strlen(member_name); // used to check that scanf consumed the full string
      int used_len;
      char varname[1000];
      // if metadata can be parsed from dataset name then we are about 6x faster
      // then using attributes
      if(sscanf(member_name, "%s it=%d tl=%d rl=%d c=%d%n", varname, &cycle, &tl, &rl, &c, &used_len) == 5 &&
         len == used_len) {
         map = 0; // single map run
      } else if(sscanf(member_name, "%s it=%d tl=%d m=%d rl=%d c=%d%n", varname, &cycle, &tl, &map, &rl, &c, &used_len) == 6 &&
                len == used_len) {
          /* no-op */
      } else if(sscanf(member_name, "%s it=%d tl=%d c=%d%n", varname, &cycle, &tl, &c, &used_len) == 4 &&
          len == used_len) {
          rl = 0; // unigrid run
         map = 0; // single map run
      } else if(sscanf(member_name, "%s it=%d tl=%d m=%d c=%d%n", varname, &cycle, &tl, &map, &c, &used_len) == 5 &&
         len == used_len) {
          rl = 0; // unigrid run
      } else {
         static bool have_warned = false;
         if(!have_warned) {
           have_warned = true;
           cout << "Could not parse name '" << member_name << "' which can slow down things considerably\n";
         }

         attrib = H5Aopen(dataset, "timestep", H5P_DEFAULT);
         H5Aread(attrib, H5T_NATIVE_INT, &cycle);
         H5Aclose(attrib);
      
         attrib = H5Aopen(dataset, "level", H5P_DEFAULT);
         H5Aread(attrib, H5T_NATIVE_INT, &rl);
         H5Aclose(attrib);
      
         //attrib = dataset.openAttribute("name");
         //attrib.read(StrType(PredType::C_S1, strlen(fullname)), &varname);
         hid_t atype = H5Tcopy(H5T_C_S1);
         H5Tset_size(atype, strlen(fullname));
         H5Tset_strpad(atype, H5T_STR_NULLTERM);
         attrib = H5Aopen(dataset, "name", H5P_DEFAULT);
         H5Aread(attrib, atype, &varname);
         H5Aclose(attrib);
         H5Tclose(atype);

         // get map-attribute from dataset-name
         map = 0;
         if (strstr(member_name, "m="))
         {
            map = atoi(strstr(member_name, "m=")+2);
         }
      }

      attrib = H5Aopen(dataset, "time", H5P_DEFAULT);
      double time = 0;
      H5Aread(attrib, H5T_NATIVE_DOUBLE, &time);
      H5Aclose(attrib);
      
      
      bool is_Cartesian = true;
      
      H5E_BEGIN_TRY {
         hid_t attrib = H5Aopen(dataset, "MapIsCartesian", H5P_DEFAULT);
         if(attrib >= 0) {
           int MapIsCartesian = 1;
            H5Aread(attrib, H5T_NATIVE_INT, &MapIsCartesian);
            is_Cartesian = MapIsCartesian;
            H5Aclose(attrib);
         }
      } H5E_END_TRY;
      
      double delta[3] = {1.0,1.0,1.0};
      H5E_BEGIN_TRY {
         hid_t attrib = H5Aopen(dataset, "delta", H5P_DEFAULT);
         if(attrib >= 0) {
            H5Aread(attrib, H5T_NATIVE_DOUBLE, delta);
            H5Aclose(attrib);
         }
      } H5E_END_TRY;
      
      attrib = H5Aopen(dataset, "iorigin", H5P_DEFAULT);
      int iorigin[3] = {0,0,0};
      H5Aread(attrib, H5T_NATIVE_INT, iorigin);
      H5Aclose(attrib);
      
      int cctk_ghosts[3] = {0,0,0};
      H5E_BEGIN_TRY {
         hid_t attrib = H5Aopen(dataset, "cctk_nghostzones", H5P_DEFAULT);
         if(attrib >= 0) {
            H5Aread(attrib, H5T_NATIVE_INT, cctk_ghosts);
            H5Aclose(attrib);
         }
      } H5E_END_TRY;
      int ghosts[6] = {0,0,0,0,0,0};
      for (int i=0; i < dim; ++i)
        ghosts[2*i] = ghosts[2*i+1] = cctk_ghosts[i];
      
      double origin[3] = {0,0,0};
      H5E_BEGIN_TRY {
         hid_t attrib = H5Aopen(dataset, "origin", H5P_DEFAULT);
         if(attrib >= 0) {
            H5Aread(attrib, H5T_NATIVE_DOUBLE, origin);
            H5Aclose(attrib);
         } else {
            // if orign attrib is not there we probably have an array
            is_Cartesian = true;
            for (int i=0; i < dim; i++)
               origin[i] = iorigin[i];
            }
      } H5E_END_TRY;
      
      
      // factor set to one for now....will get set once timestep_t is constructed
      int factor = 1;
      
      // find out scalar type of variable
      hid_t dtype = H5Dget_type(dataset);
      H5T_class_t type_class = H5Tget_class(dtype);
      int type = 0;
      int comps = 1;
      if (type_class == H5T_COMPOUND)
      {
        type = 1;
        comps = 2;  // add complex variables as 2d-vectors
      }
      else
        type = 0;
      
      H5Tclose(dtype);
      H5Sclose(dataspace);
      H5Dclose(dataset);
      
      
      (*dsetnames).push_back(dataset_entry(fullname, varname, cycle, time, rl, map, factor, type, comps, origin, delta, iorigin, dim, dims, is_Cartesian, ghosts));
   }
   
   return 0;
}


hid_t avtCarpetHDF5FileFormat::file_t::gethandle()
{
   if (file < 0) {
      struct stat s;

      file = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
      if (stat(filename.c_str(), &s) == 0) {
        mtime = s.st_mtime;
      } else {
        mtime = 0;
      }
   }
   return file;
}

void avtCarpetHDF5FileFormat::file_t::open(const char* fname)
{
   // get hierarchy datasetnames contained only in this file
   if (haveIndex(fname)) {
     struct stat s;

     hid_t index_file = H5Fopen(indexFilename(fname).c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
     if (stat(indexFilename(fname).c_str(), &s) == 0) {
       idx_mtime = s.st_mtime;
     } else {
       idx_mtime = 0;
     }
     //index_file.iterateElems("/", NULL, H5iter, &dsetnames);
     H5Giterate (index_file, "/", NULL, H5iter, &dsetnames);
     H5Fclose(index_file);
   } 
   else {
     //file->iterateElems("/", NULL, H5iter, &dsetnames);
     H5Giterate (gethandle(), "/", NULL, H5iter, &dsetnames);
   }

   // try to free memory allocated during file traversal
   H5garbage_collect();
}

void avtCarpetHDF5FileFormat::file_t::close()
{
   if (file >= 0) {
      H5Fclose(file);
      file = -1;
   }
}

bool avtCarpetHDF5FileFormat::file_t::file_changed_on_disk()
{
    struct stat s;

    if (stat(filename.c_str(), &s) != 0 || s.st_mtime > mtime) {
      return true;
    }
    if (haveIndex(filename) &&
        (stat(indexFilename(filename).c_str(), &s) != 0 || s.st_mtime > idx_mtime)) {
      return true;
    }

    return false;
}



void avtCarpetHDF5FileFormat::file_t::set_varnames(const vector<string>& v, const vector<int>& cycles)
{
   varnames = v;
}






void avtCarpetHDF5FileFormat::file_t::get_data(const dataset_entry& dset, float** data, const bool removeGhosts)
{
   hid_t dataset = H5Dopen(gethandle(), dset.name().c_str(), H5P_DEFAULT);
      
   hid_t dataspace = H5Dget_space(dataset);
   hid_t attrib;
   
   int rank  = H5Sget_simple_extent_ndims(dataspace);
      
   hsize_t dims_out[3] = {1, 1, 1};
   int ndims  = H5Sget_simple_extent_dims(dataspace, dims_out, NULL);
      
   //hsize_t      offset[3] = { dset.ghosts(removeGhosts)[4], dset.ghosts(removeGhosts)[2], dset.ghosts(removeGhosts)[0] };   // hyperslab offset in the file
   //hsize_t      count[3] = { dset.dims(removeGhosts)[0], dset.dims(removeGhosts)[1],  dset.dims(removeGhosts)[2] };    // size of the hyperslab in the file
   
   //DataSpace memspace( rank, &dset.dims(removeGhosts).front() );
   
   hsize_t      offset[3] = { 0, 0, 0 };
   hsize_t      count[3] = { 1, 1, 1 };
   hsize_t      dims[3] = { 1, 1, 1 };
  
   if (ndims == 3)
   {
      offset[0] = dset.ghosts(removeGhosts)[4];
      offset[1] = dset.ghosts(removeGhosts)[2];
      offset[2] = dset.ghosts(removeGhosts)[0];
      count[0] = dset.dims(removeGhosts)[0];
      count[1] = dset.dims(removeGhosts)[1];
      count[2] = dset.dims(removeGhosts)[2];
      dims[0] = dset.dims(removeGhosts)[0];
      dims[1] = dset.dims(removeGhosts)[1];
      dims[2] = dset.dims(removeGhosts)[2];
   }
   if (ndims == 2)
   {
      offset[0] = dset.ghosts(removeGhosts)[2];
      offset[1] = dset.ghosts(removeGhosts)[0];
      count[0] = dset.dims(removeGhosts)[1];
      count[1] = dset.dims(removeGhosts)[2];
      dims[0] = dset.dims(removeGhosts)[1];
      dims[1] = dset.dims(removeGhosts)[2];      
   }
   
   //DataSpace memspace( rank, dims );
   hid_t memspace = H5Screate_simple(rank, dims, NULL);
   
   //*data = new float[memspace.getSimpleExtentNpoints()];
   *data = new float[H5Sget_simple_extent_npoints(memspace)];
   
   //dataspace.selectHyperslab( H5S_SELECT_SET, count, offset );
   H5Sselect_hyperslab(dataspace, H5S_SELECT_SET, offset, NULL, count, NULL);
   //dataset.read( *data, PredType::NATIVE_FLOAT, memspace, dataspace );
   H5Dread(dataset, H5T_NATIVE_FLOAT, memspace, dataspace, H5P_DEFAULT, *data);
   
   H5Sclose(memspace);
   H5Sclose(dataspace);
   H5Dclose(dataset);
}




void avtCarpetHDF5FileFormat::file_t::get_complex_data(const dataset_entry& dset, float*** data, const bool removeGhosts)
{
/*
   DataSet dataset = file->openDataSet(dset.name());
   
   DataSpace dataspace = dataset.getSpace();
   
   int rank = dataspace.getSimpleExtentNdims();
   
   hsize_t dims_out[3] = {1, 1, 1};
   int ndims = dataspace.getSimpleExtentDims( dims_out, NULL);
   
   hsize_t      offset[3] = { dset.ghosts(removeGhosts)[4], dset.ghosts(removeGhosts)[2], dset.ghosts(removeGhosts)[0] };   // hyperslab offset in the file
   hsize_t      count[3] = { dset.dims(removeGhosts)[0], dset.dims(removeGhosts)[1], dset.dims(removeGhosts)[2] };    // size of the hyperslab in the file
   
   DataSpace memspace( rank, &dset.dims(removeGhosts).front() );
   
   *data = new float*[dims_out[0]*dims_out[1]*dims_out[2]];
   for (int i=0; i < dims_out[0]*dims_out[1]*dims_out[2]; ++i)
     (*data)[i] = new float[2];
   
   dataspace.selectHyperslab( H5S_SELECT_SET, count, offset );
   
   CompType type(2*sizeof(float));
   type.insertMember("real", 0, PredType::NATIVE_FLOAT);
   type.insertMember("imag", sizeof(float), PredType::NATIVE_FLOAT);
   
   float* buffer = new float[2*memspace.getSimpleExtentNpoints()];
   
   dataset.read( buffer, type, memspace, dataspace );
   
   for (int i=0; i < memspace.getSimpleExtentNpoints(); ++i)
   {
      (*data)[i][0] = buffer[2*i];
      (*data)[i][1] = buffer[2*i+1];
   }
   
   delete [] buffer;
*/
}







#ifdef USE_SELECTIONS
void avtCarpetHDF5FileFormat::RegisterDataSelections(
    const vector<avtDataSelection_p> &selections,
    vector<bool> *selectionsApplied)
{
  this->selections = selections;
  this->selectionsApplied = selectionsApplied;
}


void avtCarpetHDF5FileFormat::ProcessDataSelections()
{
  assert(selections.size() == (*selectionsApplied).size());

  (*selectionsApplied).assign((*selectionsApplied).size(), false);

//    selList     = selections;
//    selsApplied = selectionsApplied;
  avtSpatialBoxSelection composedSelection;

  for (size_t i = 0; i < selections.size(); i++) {
    if (string(selections[i]->GetType()) == "Spatial Box Data Selection") {
      const avtSpatialBoxSelection *selection = static_cast<const avtSpatialBoxSelection *>(*selections[i]);
      avtSpatialBoxSelection::InclusionMode inclusionMode = selection->GetInclusionMode();
      if (inclusionMode == avtSpatialBoxSelection::Whole or inclusionMode == avtSpatialBoxSelection::Partial) {

#if 0
        // default box ranges from FLT_MIN to FLT_MAX
        double mins[3], maxs[3];
        selection->GetMins(mins);
        selection->GetMaxs(maxs);
        cout << "Spatial Box Data Selection: " << mins[0] << "/" << maxs[0] << ", " << mins[1] << "/" << maxs[1] << ", " << mins[2] << "/" << maxs[2] << endl;
#endif

        composedSelection.Compose(*selection, composedSelection);
        (*selectionsApplied)[i] = true;
      }
    } else {
      cerr << "Unsupported selection type: " << selections[i]->GetType() << endl;
    }
  }

  double mins[3], maxs[3];
  composedSelection.GetMins(mins);
  composedSelection.GetMaxs(maxs);
  cout << "    composed Selection: " << mins[0] << "/" << maxs[0] << ", " << mins[1] << "/" << maxs[1] << ", " << mins[2] << "/" << maxs[2] << endl;
}
#endif






void*
avtCarpetHDF5FileFormat::GetAuxiliaryData(const char *var, int timestep, int domain,
                              const char *type, void *, DestructorFunction &df)
{
    void *rv = NULL;
    if (strcmp(type, AUXILIARY_DATA_SPATIAL_EXTENTS) == 0)
    {
        // Read the number of domains for the mesh.
      if (strcmp(var, "Carpet AMR-grid") == 0) {
        int ndoms = data_file->timesteps[timestep].cart_comp[0].size();
        
        // Read the spatial extents for each domain of the
        // mesh. This information should be in a single
        // and should be available without having to
        // read the real data. The expected format for
        // the data in the spatialextents array is to
        // repeat the following pattern for each domain:
        // xmin, xmax, ymin, ymax, zmin, zmax.
        vector<double> spatialextents(ndoms * 6);
        // READ ndoms*6 DOUBLE VALUES INTO spatialextents ARRAY.
        for (int i=0; i < ndoms; ++i)
        {
           spatialextents[i*6  ] = data_file->timesteps[timestep].cart_comp[0][i].xmin();
           spatialextents[i*6+1] = data_file->timesteps[timestep].cart_comp[0][i].xmax();
           spatialextents[i*6+2] = data_file->timesteps[timestep].cart_comp[0][i].ymin();
           spatialextents[i*6+3] = data_file->timesteps[timestep].cart_comp[0][i].ymax();
           spatialextents[i*6+4] = data_file->timesteps[timestep].cart_comp[0][i].zmin();
           spatialextents[i*6+5] = data_file->timesteps[timestep].cart_comp[0][i].zmax();
        }
        
        // Create an interval tree
        avtIntervalTree *itree = new avtIntervalTree(ndoms, 3);
        vector<double> extents = spatialextents;
        for(int dom = 0; dom < ndoms; ++dom)
        {
           itree->AddElement(dom, &extents[dom*6]);
        }
        itree->Calculate(true);

        // Set return values
        rv = (void *)itree;
        df = avtIntervalTree::Destruct;
      }
    }
/*    else if (strcmp(type, AUXILIARY_DATA_DATA_EXTENTS) == 0)
    {
        if (ignoreDataExtents)
        {
            debug1 << "Read options ignore request for data extents" << endl;
            return 0;
        }
        rv = (void *) GetDataExtents(var);
        df = avtIntervalTree::Destruct;
    }*/


    return rv;
}