/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkAMREnzoReaderInternal.cxx

  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/
#include "vtkAMREnzoReaderInternal.h"

#define H5_USE_16_API
#include "vtk_hdf5.h" // for the HDF5 library

#include "vtkCellData.h"
#include "vtkDataArray.h"
#include "vtkDataSet.h"
#include "vtkDoubleArray.h"
#include "vtkFloatArray.h"
#include "vtkIntArray.h"
#include "vtkLongArray.h"
#include "vtkLongLongArray.h"
#include "vtkPointData.h"
#include "vtkPolyData.h"
#include "vtkShortArray.h"
#include "vtkUnsignedCharArray.h"
#include "vtkUnsignedIntArray.h"
#include "vtkUnsignedShortArray.h"
#include "vtksys/FStream.hxx"
#include "vtksys/SystemTools.hxx"

//------------------------------------------------------------------------------
//                       Functions for Parsing File Names
//------------------------------------------------------------------------------

static std::string GetEnzoMajorFileName(const char* path)
{
  return (vtksys::SystemTools::GetFilenameName(std::string(path)));
}

//------------------------------------------------------------------------------
//                       Class vtkEnzoReaderBlock (begin)
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void vtkEnzoReaderBlock::Init()
{
  this->BlockFileName = "";
  this->ParticleFileName = "";

  this->Index = -1;
  this->Level = -1;
  this->ParentId = -1;
  this->ChildrenIds.clear();
  this->NumberOfParticles = 0;
  this->NumberOfDimensions = 0;

  this->MinParentWiseIds[0] = this->MinParentWiseIds[1] = this->MinParentWiseIds[2] =
    this->MaxParentWiseIds[0] = this->MaxParentWiseIds[1] = this->MaxParentWiseIds[2] = -1;

  this->MinLevelBasedIds[0] = this->MinLevelBasedIds[1] = this->MinLevelBasedIds[2] =
    this->MaxLevelBasedIds[0] = this->MaxLevelBasedIds[1] = this->MaxLevelBasedIds[2] = -1;

  this->BlockCellDimensions[0] = this->BlockCellDimensions[1] = this->BlockCellDimensions[2] =
    this->BlockNodeDimensions[0] = this->BlockNodeDimensions[1] = this->BlockNodeDimensions[2] = 0;

  this->MinBounds[0] = this->MinBounds[1] = this->MinBounds[2] = VTK_DOUBLE_MAX;
  this->MaxBounds[0] = this->MaxBounds[1] = this->MaxBounds[2] = -VTK_DOUBLE_MAX;

  this->SubdivisionRatio[0] = this->SubdivisionRatio[1] = this->SubdivisionRatio[2] = 1.0;
}

//------------------------------------------------------------------------------
void vtkEnzoReaderBlock::DeepCopy(const vtkEnzoReaderBlock* other)
{
  this->BlockFileName = other->BlockFileName;
  this->ParticleFileName = other->ParticleFileName;

  this->Index = other->Index;
  this->Level = other->Level;
  this->ParentId = other->ParentId;
  this->ChildrenIds = other->ChildrenIds;
  this->NumberOfParticles = other->NumberOfParticles;
  this->NumberOfDimensions = other->NumberOfDimensions;

  this->MinParentWiseIds[0] = other->MinParentWiseIds[0];
  this->MinParentWiseIds[1] = other->MinParentWiseIds[1];
  this->MinParentWiseIds[2] = other->MinParentWiseIds[2];
  this->MaxParentWiseIds[0] = other->MaxParentWiseIds[0];
  this->MaxParentWiseIds[1] = other->MaxParentWiseIds[1];
  this->MaxParentWiseIds[2] = other->MaxParentWiseIds[2];

  this->MinLevelBasedIds[0] = other->MinLevelBasedIds[0];
  this->MinLevelBasedIds[1] = other->MinLevelBasedIds[1];
  this->MinLevelBasedIds[2] = other->MinLevelBasedIds[2];
  this->MaxLevelBasedIds[0] = other->MaxLevelBasedIds[0];
  this->MaxLevelBasedIds[1] = other->MaxLevelBasedIds[1];
  this->MaxLevelBasedIds[2] = other->MaxLevelBasedIds[2];

  this->BlockCellDimensions[0] = other->BlockCellDimensions[0];
  this->BlockCellDimensions[1] = other->BlockCellDimensions[1];
  this->BlockCellDimensions[2] = other->BlockCellDimensions[2];
  this->BlockNodeDimensions[0] = other->BlockNodeDimensions[0];
  this->BlockNodeDimensions[1] = other->BlockNodeDimensions[1];
  this->BlockNodeDimensions[2] = other->BlockNodeDimensions[2];

  this->MinBounds[0] = other->MinBounds[0];
  this->MinBounds[1] = other->MinBounds[1];
  this->MinBounds[2] = other->MinBounds[2];
  this->MaxBounds[0] = other->MaxBounds[0];
  this->MaxBounds[1] = other->MaxBounds[1];
  this->MaxBounds[2] = other->MaxBounds[2];

  this->SubdivisionRatio[0] = other->SubdivisionRatio[0];
  this->SubdivisionRatio[1] = other->SubdivisionRatio[1];
  this->SubdivisionRatio[2] = other->SubdivisionRatio[2];
}

//------------------------------------------------------------------------------
// get the bounding (cell) Ids of this block in terms of its parent block's
// sub-division resolution (indexing is limited to the scope of the parent)
void vtkEnzoReaderBlock::GetParentWiseIds(std::vector<vtkEnzoReaderBlock>& blocks)
{
  if (this->ParentId != 0)
  {
    // the parent is not the root and then we need to determine the offset
    // (in terms of the number of parent divisions / cells) of the current
    // block's beginning / ending position relative to the parent block's
    // beginning position
    vtkEnzoReaderBlock& parent = blocks[this->ParentId];
    this->MinParentWiseIds[0] = static_cast<int>(0.5 +
      parent.BlockCellDimensions[0] * (this->MinBounds[0] - parent.MinBounds[0]) /
        (parent.MaxBounds[0] - parent.MinBounds[0]));
    this->MaxParentWiseIds[0] = static_cast<int>(0.5 +
      parent.BlockCellDimensions[0] * (this->MaxBounds[0] - parent.MinBounds[0]) /
        (parent.MaxBounds[0] - parent.MinBounds[0]));

    this->MinParentWiseIds[1] = static_cast<int>(0.5 +
      parent.BlockCellDimensions[1] * (this->MinBounds[1] - parent.MinBounds[1]) /
        (parent.MaxBounds[1] - parent.MinBounds[1]));
    this->MaxParentWiseIds[1] = static_cast<int>(0.5 +
      parent.BlockCellDimensions[1] * (this->MaxBounds[1] - parent.MinBounds[1]) /
        (parent.MaxBounds[1] - parent.MinBounds[1]));

    if (this->NumberOfDimensions == 3)
    {
      this->MinParentWiseIds[2] = static_cast<int>(0.5 +
        parent.BlockCellDimensions[2] * (this->MinBounds[2] - parent.MinBounds[2]) /
          (parent.MaxBounds[2] - parent.MinBounds[2]));
      this->MaxParentWiseIds[2] = static_cast<int>(0.5 +
        parent.BlockCellDimensions[2] * (this->MaxBounds[2] - parent.MinBounds[2]) /
          (parent.MaxBounds[2] - parent.MinBounds[2]));
    }
    else
    {
      this->MinParentWiseIds[2] = 0;
      this->MaxParentWiseIds[2] = 0;
    }

    // the ratio for mapping two parent-wise Ids to 0 and
    // this->BlockCellDimension[i],
    // respectively, while the same region is covered
    this->SubdivisionRatio[0] = static_cast<double>(this->BlockCellDimensions[0]) /
      (this->MaxParentWiseIds[0] - this->MinParentWiseIds[0]);
    this->SubdivisionRatio[1] = static_cast<double>(this->BlockCellDimensions[1]) /
      (this->MaxParentWiseIds[1] - this->MinParentWiseIds[1]);

    if (this->NumberOfDimensions == 3)
    {
      this->SubdivisionRatio[2] = static_cast<double>(this->BlockCellDimensions[2]) /
        (this->MaxParentWiseIds[2] - this->MinParentWiseIds[2]);
    }
    else
    {
      this->SubdivisionRatio[2] = 1.0;
    }
  }
  else
  {
    // Now that the parent is the root, it can not provide cell-dimensions
    // information (BlockCellDimensions[0 .. 2]) directly, as the above does.
    // Fortunately we can obtain it according to the spatial ratio of the
    // child block (the current one) to the parent (root) and the child block's
    // cell-dimensions information. This derivation is based on the definition
    // of 'level' that all children blocks at the same level (e.g., the current
    // block and its siblings) share the same sub-division ratio relative to
    // their parent (the root herein).
    vtkEnzoReaderBlock& block0 = blocks[0];

    double xRatio =
      (this->MaxBounds[0] - this->MinBounds[0]) / (block0.MaxBounds[0] - block0.MinBounds[0]);
    this->MinParentWiseIds[0] = static_cast<int>(0.5 +
      (this->BlockCellDimensions[0] / xRatio) // parent's dim
        * (this->MinBounds[0] - block0.MinBounds[0]) / (block0.MaxBounds[0] - block0.MinBounds[0]));
    this->MaxParentWiseIds[0] = static_cast<int>(0.5 +
      (this->BlockCellDimensions[0] / xRatio) * (this->MaxBounds[0] - block0.MinBounds[0]) /
        (block0.MaxBounds[0] - block0.MinBounds[0]));

    double yRatio =
      (this->MaxBounds[1] - this->MinBounds[1]) / (block0.MaxBounds[1] - block0.MinBounds[1]);
    this->MinParentWiseIds[1] = static_cast<int>(0.5 +
      (this->BlockCellDimensions[1] / yRatio) * (this->MinBounds[1] - block0.MinBounds[1]) /
        (block0.MaxBounds[1] - block0.MinBounds[1]));
    this->MaxParentWiseIds[1] = static_cast<int>(0.5 +
      (this->BlockCellDimensions[1] / yRatio) * (this->MaxBounds[1] - block0.MinBounds[1]) /
        (block0.MaxBounds[1] - block0.MinBounds[1]));

    if (this->NumberOfDimensions == 3)
    {
      double zRatio =
        (this->MaxBounds[2] - this->MinBounds[2]) / (block0.MaxBounds[2] - block0.MinBounds[2]);
      this->MinParentWiseIds[2] = static_cast<int>(0.5 +
        (this->BlockCellDimensions[2] / zRatio) * (this->MinBounds[2] - block0.MinBounds[2]) /
          (block0.MaxBounds[2] - block0.MinBounds[2]));
      this->MaxParentWiseIds[2] = static_cast<int>(0.5 +
        (this->BlockCellDimensions[2] / zRatio) * (this->MaxBounds[2] - block0.MinBounds[2]) /
          (block0.MaxBounds[2] - block0.MinBounds[2]));
    }
    else
    {
      this->MinParentWiseIds[2] = 0;
      this->MaxParentWiseIds[2] = 0;
    }

    this->SubdivisionRatio[0] = 1.0;
    this->SubdivisionRatio[1] = 1.0;
    this->SubdivisionRatio[2] = 1.0;
  }
}

//------------------------------------------------------------------------------
void vtkEnzoReaderBlock::GetLevelBasedIds(std::vector<vtkEnzoReaderBlock>& blocks)
{
  // note that this function is invoked from the root in a top-down manner
  // and the parent-wise Ids have been determined in advance

  if (this->ParentId != 0)
  {
    // the parent is not the root and therefore we need to exploit the level-
    // based Ids of the parent, of which the shifted version is multiplied by
    // the refinement ratio

    vtkEnzoReaderBlock& parent = blocks[this->ParentId];
    this->MinLevelBasedIds[0] = static_cast<int>(
      (parent.MinLevelBasedIds[0] + this->MinParentWiseIds[0]) * this->SubdivisionRatio[0]);

    this->MinLevelBasedIds[1] = static_cast<int>(
      (parent.MinLevelBasedIds[1] + this->MinParentWiseIds[1]) * this->SubdivisionRatio[1]);
    this->MinLevelBasedIds[2] = static_cast<int>(
      (parent.MinLevelBasedIds[2] + this->MinParentWiseIds[2]) * this->SubdivisionRatio[2]);

    this->MaxLevelBasedIds[0] = static_cast<int>(
      (parent.MinLevelBasedIds[0] + this->MaxParentWiseIds[0]) * this->SubdivisionRatio[0]);
    this->MaxLevelBasedIds[1] = static_cast<int>(
      (parent.MinLevelBasedIds[1] + this->MaxParentWiseIds[1]) * this->SubdivisionRatio[1]);
    this->MaxLevelBasedIds[2] = static_cast<int>(
      (parent.MinLevelBasedIds[2] + this->MaxParentWiseIds[2]) * this->SubdivisionRatio[2]);
  }
  else
  {
    // now that the parent is the root, the parent-wise Ids
    // are just the level-based Ids
    this->MinLevelBasedIds[0] = this->MinParentWiseIds[0];
    this->MinLevelBasedIds[1] = this->MinParentWiseIds[1];
    this->MinLevelBasedIds[2] = this->MinParentWiseIds[2];

    this->MaxLevelBasedIds[0] = this->MaxParentWiseIds[0];
    this->MaxLevelBasedIds[1] = this->MaxParentWiseIds[1];
    this->MaxLevelBasedIds[2] = this->MaxParentWiseIds[2];
  }
}
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
//                       Class vtkEnzoReaderBlock ( end )
//------------------------------------------------------------------------------

//------------------------------------------------------------------------------
//                     Class  vtkEnzoReaderInternal (begin)
//------------------------------------------------------------------------------

vtkEnzoReaderInternal::vtkEnzoReaderInternal()
{
  this->Init();
}

//------------------------------------------------------------------------------
vtkEnzoReaderInternal::~vtkEnzoReaderInternal()
{
  this->ReleaseDataArray();
  this->Init();
  this->FileName = nullptr;
}

//------------------------------------------------------------------------------
void vtkEnzoReaderInternal::Init()
{
  this->DataTime = 0.0;
  this->FileName = nullptr;
  // this->TheReader  = nullptr;
  this->DataArray = nullptr;
  this->CycleIndex = 0;

  this->ReferenceBlock = 0;
  this->NumberOfBlocks = 0;
  this->NumberOfLevels = 0;
  this->NumberOfDimensions = 0;
  this->NumberOfMultiBlocks = 0;

  this->DirectoryName = "";
  this->MajorFileName = "";
  this->BoundaryFileName = "";
  this->HierarchyFileName = "";

  this->Blocks.clear();
  this->BlockAttributeNames.clear();
  this->ParticleAttributeNames.clear();
  this->TracerParticleAttributeNames.clear();
}

//------------------------------------------------------------------------------
void vtkEnzoReaderInternal::ReleaseDataArray()
{
  if (this->DataArray)
  {
    this->DataArray->Delete();
    this->DataArray = nullptr;
  }
}

//------------------------------------------------------------------------------
int vtkEnzoReaderInternal::GetBlockAttribute(
  const char* attribute, int blockIdx, vtkDataSet* pDataSet)
{

  // this function must be called by GetBlock( ... )
  this->ReadMetaData();

  if (attribute == nullptr || blockIdx < 0 || pDataSet == nullptr ||
    blockIdx >= this->NumberOfBlocks)
  {
    return 0;
  }

  // try obtaining the attribute and attaching it to the grid as a cell data
  // NOTE: the 'equal' comparison below is MUST because in some cases (such
  // as cosmological datasets) not all rectilinear blocks contain an assumably
  // common block attribute. This is the case where such a block attribute
  // may result from a particles-to-cells interpolation process. Thus those
  // blocks with particles (e.g., the reference one with the fewest cells)
  // contain it as a block attribute, whereas others without particles just
  // do not contain this block attribute.
  int succeeded = 0;
  if (this->LoadAttribute(attribute, blockIdx) &&
    (pDataSet->GetNumberOfCells() == this->DataArray->GetNumberOfTuples()))
  {
    succeeded = 1;
    pDataSet->GetCellData()->AddArray(this->DataArray);
    this->ReleaseDataArray();
  }

  return succeeded;
}

//------------------------------------------------------------------------------
int vtkEnzoReaderInternal::LoadAttribute(const char* attribute, int blockIdx)
{
  // TODO: implement this
  // called by GetBlockAttribute( ... ) or GetParticlesAttribute()
  this->ReadMetaData();

  if (attribute == nullptr || blockIdx < 0 || blockIdx >= this->NumberOfBlocks)
  {
    return 0;
  }

  // this->Internal->Blocks includes a pseudo block --- the root as block #0
  blockIdx++;

  // currently only the HDF5 file format is supported
  std::string blckFile = this->Blocks[blockIdx].BlockFileName;
  hid_t fileIndx = H5Fopen(blckFile.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
  if (fileIndx < 0)
  {
    return 0;
  }

  // retrieve the contents of the root directory to look for a group
  // corresponding to the target block, if available, open that group

  int blckIndx;
  char blckName[65];
  int objIndex;
  hsize_t numbObjs;
  hid_t rootIndx = H5Gopen(fileIndx, "/");
  H5Gget_num_objs(rootIndx, &numbObjs);
  for (objIndex = 0; objIndex < static_cast<int>(numbObjs); objIndex++)
  {
    int type = H5Gget_objtype_by_idx(rootIndx, objIndex);
    if (type == H5G_GROUP)
    {
      H5Gget_objname_by_idx(rootIndx, objIndex, blckName, 64);
      if ((sscanf(blckName, "Grid%d", &blckIndx) == 1) &&
        ((blckIndx == blockIdx) || (blckIndx == blockIdx + 1)))
      {
        // located the target block
        rootIndx = H5Gopen(rootIndx, blckName);
        break;
      }
    }
  } // END for all objects

  // disable error messages while looking for the attribute (if any) name
  // and enable error messages when it is done
  void* pContext = nullptr;
  H5E_auto_t erorFunc;
  H5Eget_auto(&erorFunc, &pContext);
  H5Eset_auto(nullptr, nullptr);

  hid_t attrIndx = H5Dopen(rootIndx, attribute);

  H5Eset_auto(erorFunc, pContext);
  pContext = nullptr;

  // check if the data attribute exists
  if (attrIndx < 0)
  {
    vtkGenericWarningMacro(
      "Attribute (" << attribute << ") data does not exist in file " << blckFile.c_str());
    H5Gclose(rootIndx);
    H5Fclose(fileIndx);
    return 0;
  }

  // get cell dimensions and the valid number
  hsize_t cellDims[3];
  hid_t spaceIdx = H5Dget_space(attrIndx);
  H5Sget_simple_extent_dims(spaceIdx, cellDims, nullptr);
  hsize_t numbDims = H5Sget_simple_extent_ndims(spaceIdx);

  // number of attribute tuples = number of cells (or particles)
  int numTupls = 0;
  switch (numbDims)
  {
    case 1:
      numTupls = cellDims[0];
      break;
    case 2:
      numTupls = cellDims[0] * cellDims[1];
      break;
    case 3:
      numTupls = cellDims[0] * cellDims[1] * cellDims[2];
      break;
    default:
      H5Gclose(spaceIdx);
      H5Fclose(attrIndx);
      H5Gclose(rootIndx);
      H5Fclose(fileIndx);
      return 0;
  }

  // determine the data type, load the values, and, if necessary, convert
  // the raw data to double type
  // DOUBLE / FLOAT  / INT  / UINT  / CHAR  / UCHAR
  // SHORT  / USHORT / LONG / ULONG / LLONG / ULLONG / LDOUBLE

  this->ReleaseDataArray(); // data array maintained by internal

  hid_t tRawType = H5Dget_type(attrIndx);
  hid_t dataType = H5Tget_native_type(tRawType, H5T_DIR_ASCEND);
  if (H5Tequal(dataType, H5T_NATIVE_FLOAT))
  {
    this->DataArray = vtkFloatArray::New();
    this->DataArray->SetNumberOfTuples(numTupls);
    float* arrayPtr =
      static_cast<float*>(vtkArrayDownCast<vtkFloatArray>(this->DataArray)->GetPointer(0));
    H5Dread(attrIndx, dataType, H5S_ALL, H5S_ALL, H5P_DEFAULT, arrayPtr);
    arrayPtr = nullptr;
  }
  else if (H5Tequal(dataType, H5T_NATIVE_DOUBLE))
  {
    this->DataArray = vtkDoubleArray::New();
    this->DataArray->SetNumberOfTuples(numTupls);
    double* arrayPtr =
      static_cast<double*>(vtkArrayDownCast<vtkDoubleArray>(this->DataArray)->GetPointer(0));
    H5Dread(attrIndx, dataType, H5S_ALL, H5S_ALL, H5P_DEFAULT, arrayPtr);
    arrayPtr = nullptr;
  }
  else if (H5Tequal(dataType, H5T_NATIVE_INT))
  {
    this->DataArray = vtkIntArray::New();
    this->DataArray->SetNumberOfTuples(numTupls);
    int* arrayPtr =
      static_cast<int*>(vtkArrayDownCast<vtkIntArray>(this->DataArray)->GetPointer(0));
    H5Dread(attrIndx, dataType, H5S_ALL, H5S_ALL, H5P_DEFAULT, arrayPtr);
    arrayPtr = nullptr;
  }
  else if (H5Tequal(dataType, H5T_NATIVE_UINT))
  {
    this->DataArray = vtkUnsignedIntArray::New();
    this->DataArray->SetNumberOfTuples(numTupls);
    unsigned int* arrayPtr = static_cast<unsigned int*>(
      vtkArrayDownCast<vtkUnsignedIntArray>(this->DataArray)->GetPointer(0));
    H5Dread(attrIndx, dataType, H5S_ALL, H5S_ALL, H5P_DEFAULT, arrayPtr);
    arrayPtr = nullptr;
  }
  else if (H5Tequal(dataType, H5T_NATIVE_SHORT))
  {
    this->DataArray = vtkShortArray::New();
    this->DataArray->SetNumberOfTuples(numTupls);
    short* arrayPtr =
      static_cast<short*>(vtkArrayDownCast<vtkShortArray>(this->DataArray)->GetPointer(0));
    H5Dread(attrIndx, dataType, H5S_ALL, H5S_ALL, H5P_DEFAULT, arrayPtr);
    arrayPtr = nullptr;
  }
  else if (H5Tequal(dataType, H5T_NATIVE_USHORT))
  {
    this->DataArray = vtkUnsignedShortArray::New();
    this->DataArray->SetNumberOfTuples(numTupls);
    unsigned short* arrayPtr = static_cast<unsigned short*>(
      vtkArrayDownCast<vtkUnsignedShortArray>(this->DataArray)->GetPointer(0));
    H5Dread(attrIndx, dataType, H5S_ALL, H5S_ALL, H5P_DEFAULT, arrayPtr);
    arrayPtr = nullptr;
  }
  else if (H5Tequal(dataType, H5T_NATIVE_UCHAR))
  {

    this->DataArray = vtkUnsignedCharArray::New();
    this->DataArray->SetNumberOfTuples(numTupls);
    unsigned char* arrayPtr = static_cast<unsigned char*>(
      vtkArrayDownCast<vtkUnsignedCharArray>(this->DataArray)->GetPointer(0));
    H5Dread(attrIndx, dataType, H5S_ALL, H5S_ALL, H5P_DEFAULT, arrayPtr);
    arrayPtr = nullptr;
  }
  else if (H5Tequal(dataType, H5T_NATIVE_LONG))
  {
    this->DataArray = vtkLongArray::New();
    this->DataArray->SetNumberOfTuples(numTupls);
    long* arrayPtr =
      static_cast<long*>(vtkArrayDownCast<vtkLongArray>(this->DataArray)->GetPointer(0));
    H5Dread(attrIndx, dataType, H5S_ALL, H5S_ALL, H5P_DEFAULT, arrayPtr);
    arrayPtr = nullptr;
  }
  else if (H5Tequal(dataType, H5T_NATIVE_LLONG))
  {
    this->DataArray = vtkLongLongArray::New();
    this->DataArray->SetNumberOfTuples(numTupls);
    long long* arrayPtr =
      static_cast<long long*>(vtkArrayDownCast<vtkLongLongArray>(this->DataArray)->GetPointer(0));
    H5Dread(attrIndx, dataType, H5S_ALL, H5S_ALL, H5P_DEFAULT, arrayPtr);
    arrayPtr = nullptr;
  }
  else
  {
    //        vtkErrorMacro( "Unknown HDF5 data type --- it is not FLOAT, "       <<
    //                       "DOUBLE, INT, UNSIGNED INT, SHORT, UNSIGNED SHORT, " <<
    //                       "UNSIGNED CHAR, LONG, or LONG LONG." << endl );
    H5Tclose(dataType);
    H5Tclose(tRawType);
    H5Tclose(spaceIdx);
    H5Dclose(attrIndx);
    H5Gclose(rootIndx);
    H5Fclose(fileIndx);
    return 0;
  }

  // do not forget to provide a name for the array
  this->DataArray->SetName(attribute);

  // These close statements cause a crash!
  //  H5Tclose( dataType );
  //  H5Tclose( tRawType );
  //  H5Tclose( spaceIdx );
  //  H5Dclose( attrIndx );
  //  H5Gclose( rootIndx );
  //  H5Fclose( fileIndx );
  return 1;
}

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
// parse the hierarchy file to create block structures, including the bounding
// box, cell dimensions, grid / node dimensions, number of particles, level Id,
// block file name, and particle file name of each block
void vtkEnzoReaderInternal::ReadBlockStructures()
{
  vtksys::ifstream stream(this->HierarchyFileName.c_str());
  if (!stream)
  {
    vtkGenericWarningMacro(
      "Invalid hierarchy file name: " << this->HierarchyFileName.c_str() << endl);
    return;
  }

  // init the root block, addressing only 4 four fields
  vtkEnzoReaderBlock block0;
  block0.Index = 0;
  block0.Level = -1;
  block0.ParentId = -1;
  block0.NumberOfDimensions = this->NumberOfDimensions;
  this->Blocks.push_back(block0);

  int levlId = 0;
  int parent = 0;
  std::string theStr;

  while (stream)
  {
    while (stream && theStr != "Grid" && theStr != "Time" && theStr != "Pointer:")
    {
      stream >> theStr;
    }

    // block information
    if (theStr == "Grid")
    {
      vtkEnzoReaderBlock tmpBlk;
      tmpBlk.NumberOfDimensions = this->NumberOfDimensions;

      stream >> theStr; // '='
      stream >> tmpBlk.Index;

      // the starting and ending (cell --- not node) Ids of the block
      int minIds[3];
      int maxIds[3];
      while (theStr != "GridStartIndex")
      {
        stream >> theStr;
      }
      stream >> theStr; // '='

      if (this->NumberOfDimensions == 3)
      {
        stream >> minIds[0] >> minIds[1] >> minIds[2];
      }
      else
      {
        stream >> minIds[0] >> minIds[1];
      }

      while (theStr != "GridEndIndex")
      {
        stream >> theStr;
      }
      stream >> theStr; // '='

      if (this->NumberOfDimensions == 3)
      {
        stream >> maxIds[0] >> maxIds[1] >> maxIds[2];
      }
      else
      {
        stream >> maxIds[0] >> maxIds[1];
      }

      // the cell dimensions of the block
      tmpBlk.BlockCellDimensions[0] = maxIds[0] - minIds[0] + 1;
      tmpBlk.BlockCellDimensions[1] = maxIds[1] - minIds[1] + 1;
      if (this->NumberOfDimensions == 3)
      {
        tmpBlk.BlockCellDimensions[2] = maxIds[2] - minIds[2] + 1;
      }
      else
      {
        tmpBlk.BlockCellDimensions[2] = 1;
      }

      // the grid (node --- not means the block) dimensions of the block
      tmpBlk.BlockNodeDimensions[0] = tmpBlk.BlockCellDimensions[0] + 1;
      tmpBlk.BlockNodeDimensions[1] = tmpBlk.BlockCellDimensions[1] + 1;
      if (this->NumberOfDimensions == 3)
      {
        tmpBlk.BlockNodeDimensions[2] = tmpBlk.BlockCellDimensions[2] + 1;
      }
      else
      {
        tmpBlk.BlockNodeDimensions[2] = 1;
      }

      // the min bounding box of the block
      while (theStr != "GridLeftEdge")
      {
        stream >> theStr;
      }
      stream >> theStr; // '='

      if (this->NumberOfDimensions == 3)
      {
        stream >> tmpBlk.MinBounds[0] >> tmpBlk.MinBounds[1] >> tmpBlk.MinBounds[2];
      }
      else
      {
        tmpBlk.MinBounds[2] = 0;
        stream >> tmpBlk.MinBounds[0] >> tmpBlk.MinBounds[1];
      }

      // the max bounding box of the block
      while (theStr != "GridRightEdge")
      {
        stream >> theStr;
      }
      stream >> theStr; // '='

      if (this->NumberOfDimensions == 3)
      {
        stream >> tmpBlk.MaxBounds[0] >> tmpBlk.MaxBounds[1] >> tmpBlk.MaxBounds[2];
      }
      else
      {
        tmpBlk.MaxBounds[2] = 0;
        stream >> tmpBlk.MaxBounds[0] >> tmpBlk.MaxBounds[1];
      }

      // obtain the block file name (szName includes the full path)
      std::string szName;
      while (theStr != "BaryonFileName")
      {
        stream >> theStr;
      }
      stream >> theStr; // '='
      stream >> szName;

      //      std::cout << "szname: " << szName.c_str() << std::endl;
      //      std::cout.flush();
      tmpBlk.BlockFileName = this->DirectoryName + "/" + GetEnzoMajorFileName(szName.c_str());

      // obtain the particle file name (szName includes the full path)
      while (theStr != "NumberOfParticles")
      {
        stream >> theStr;
      }
      stream >> theStr; // '='
      stream >> tmpBlk.NumberOfParticles;

      if (tmpBlk.NumberOfParticles > 0)
      {
        while (theStr != "ParticleFileName")
        {
          stream >> theStr;
        }
        stream >> theStr; // '='
        stream >> szName;
        tmpBlk.ParticleFileName = this->DirectoryName + "/" + GetEnzoMajorFileName(szName.c_str());
      }

      tmpBlk.Level = levlId;
      tmpBlk.ParentId = parent;

      if (static_cast<int>(this->Blocks.size()) != tmpBlk.Index)
      {
        vtkGenericWarningMacro("The blocks in the hierarchy file "
          << this->HierarchyFileName.c_str() << " are currently expected to be "
          << " listed in order." << endl);
        return;
      }

      this->Blocks.push_back(tmpBlk);
      this->Blocks[parent].ChildrenIds.push_back(tmpBlk.Index);
      this->NumberOfBlocks = static_cast<int>(this->Blocks.size()) - 1;
    }
    else if (theStr == "Pointer:")
    {
      theStr = "";
      int tmpInt;
      char tmpChr;
      while ((tmpChr = stream.get()) != '[')
        ;
      while ((tmpChr = stream.get()) != ']')
        theStr += tmpChr;

      int blkIdx = atoi(theStr.c_str());
      stream.get(); // -
      stream.get(); // >
      stream >> theStr;
      if (theStr == "NextGridNextLevel")
      {
        stream >> theStr; // '='
        stream >> tmpInt;
        if (tmpInt != 0)
        {
          levlId = this->Blocks[blkIdx].Level + 1;
          this->NumberOfLevels =
            (levlId + 1 > this->NumberOfLevels) ? (levlId + 1) : this->NumberOfLevels;
          parent = blkIdx;
        }
      }
      else // theStr == "NextGridThisLevel"
      {
        stream >> theStr; // '='
        stream >> tmpInt;
      }
    }
    else if (theStr == "Time")
    {
      stream >> theStr; // '='
      stream >> this->DataTime;
    }

    stream >> theStr;
  }

  if (this->NumberOfLevels == 0 && this->NumberOfBlocks > 0)
  {
    // unigrid dataset, change the number of levels to 1
    this->NumberOfLevels = 1;
  }

  stream.close();
}

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
// obtain the general information of the dataset (number of dimensions)
void vtkEnzoReaderInternal::ReadGeneralParameters()
{
  vtksys::ifstream stream(this->MajorFileName.c_str());
  if (!stream)
  {
    vtkGenericWarningMacro("Invalid parameter file " << this->MajorFileName.c_str() << endl);
    return;
  }

  std::string tmpStr;
  while (stream)
  {
    stream >> tmpStr;

    if (tmpStr == "InitialCycleNumber")
    {
      stream >> tmpStr; // '='
      stream >> this->CycleIndex;
    }
    else if (tmpStr == "InitialTime")
    {
      stream >> tmpStr; // '='
      stream >> this->DataTime;
    }
    else if (tmpStr == "TopGridRank")
    {
      stream >> tmpStr; // '='
      stream >> this->NumberOfDimensions;
    }
  }

  stream.close();
}

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
// get the bounding box of the root block based on those of its descendants
void vtkEnzoReaderInternal::DetermineRootBoundingBox()
{
  vtkEnzoReaderBlock& block0 = this->Blocks[0];

  // now loop over all level zero grids
  for (int blkIdx = 1; blkIdx <= this->NumberOfBlocks && this->Blocks[blkIdx].ParentId == 0;
       blkIdx++)
    for (int dimIdx = 0; dimIdx < this->NumberOfDimensions; dimIdx++)
    {
      block0.MinBounds[dimIdx] = (this->Blocks[blkIdx].MinBounds[dimIdx] < block0.MinBounds[dimIdx])
        ? this->Blocks[blkIdx].MinBounds[dimIdx]
        : block0.MinBounds[dimIdx];

      block0.MaxBounds[dimIdx] = (this->Blocks[blkIdx].MaxBounds[dimIdx] > block0.MaxBounds[dimIdx])
        ? this->Blocks[blkIdx].MaxBounds[dimIdx]
        : block0.MaxBounds[dimIdx];
    }
}

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
// perform an initial collection of attribute names (for block and particles)
void vtkEnzoReaderInternal::GetAttributeNames()
{
  int wasFound = 0;       // any block with particles was found?
  int blkIndex = 0;       // index of the block with fewest cells
                          // (either with or without particles)
  int numCells = INT_MAX; // number of cells of a block
  int numbBlks = static_cast<int>(this->Blocks.size());

  for (int i = 1; i < numbBlks; i++)
  {
    vtkEnzoReaderBlock& tmpBlock = this->Blocks[i];
    if (wasFound && (tmpBlock.NumberOfParticles <= 0))
    {
      continue;
    }

    int tempNumb = tmpBlock.BlockCellDimensions[0] * tmpBlock.BlockCellDimensions[1] *
      tmpBlock.BlockCellDimensions[2];

    if ((tempNumb < numCells) || (!wasFound && tmpBlock.NumberOfParticles > 0))
    {
      if (!wasFound || (tmpBlock.NumberOfParticles > 0))
      {
        numCells = tempNumb;
        blkIndex = tmpBlock.Index;
        wasFound = (tmpBlock.NumberOfParticles > 0) ? 1 : 0;
      }
    }
  }
  this->ReferenceBlock = blkIndex;

  // open the block file
  std::string blckFile = this->Blocks[blkIndex].BlockFileName;
  hid_t fileIndx = H5Fopen(blckFile.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);

  if (fileIndx < 0)
  {
    vtkGenericWarningMacro("Failed to open HDF5 grid file " << blckFile.c_str());
    return;
  }

  // retrieve the contents of the root directory to look for a group
  // corresponding to the specified block name (the one with the fewest
  // cells --- either with or without particles) and, if available, open
  // that group

  int objIndex;
  hsize_t numbObjs;
  hid_t rootIndx = H5Gopen(fileIndx, "/");
  H5Gget_num_objs(rootIndx, &numbObjs);

  for (objIndex = 0; objIndex < static_cast<int>(numbObjs); objIndex++)
  {
    if (H5Gget_objtype_by_idx(rootIndx, objIndex) == H5G_GROUP)
    {
      int blckIndx;
      char blckName[65];
      H5Gget_objname_by_idx(rootIndx, objIndex, blckName, 64);

      if (sscanf(blckName, "Grid%d", &blckIndx) == 1 &&
        (blckIndx == blkIndex) // does this block have the fewest cells?
      )
      {
        rootIndx = H5Gopen(rootIndx, blckName); // located the target block
        break;
      }
    }
  }

  // in case of entering a sub-directory, obtain the number of objects here
  // and proceed with the parsing work (now rootIndx points to the group of
  // of target block)
  H5Gget_num_objs(rootIndx, &numbObjs);

  for (objIndex = 0; objIndex < static_cast<int>(numbObjs); objIndex++)
  {
    if (H5Gget_objtype_by_idx(rootIndx, objIndex) == H5G_DATASET)
    {
      char tempName[65];
      H5Gget_objname_by_idx(rootIndx, objIndex, tempName, 64);

      // NOTE: to do the same diligence as HDF4 here, we should
      // really H5Dopen, H5Dget_space, H5Sget_simple_extent_ndims
      // and make sure it is a 3D (or 2D?) object before assuming
      // it is a mesh variable.  For now, assume away!
      if ((strlen(tempName) > 8) && (strncmp(tempName, "particle", 8) == 0))
      {
        // it's a particle variable and skip over coordinate arrays
        if (strncmp(tempName, "particle_position_", 18) != 0)
        {
          this->ParticleAttributeNames.emplace_back(tempName);
        }
      }
      else if ((strlen(tempName) > 16) && (strncmp(tempName, "tracer_particles", 16) == 0))
      {
        // it's a tracer_particle variable and skip over coordinate arrays
        if (strncmp(tempName, "tracer_particle_position_", 25) != 0)
        {
          this->TracerParticleAttributeNames.emplace_back(tempName);
        }
      }
      else
      {
        this->BlockAttributeNames.emplace_back(tempName);
      }
    }
  }

  H5Gclose(rootIndx);
  H5Fclose(fileIndx);
}

//------------------------------------------------------------------------------
// This function checks the block attributes, of which some might be actually
// particle attributes since a flexible (not standard) attributes naming scheme
// (such as the one adopted in cosmological datasets) causes this Enzo reader,
// specifically function GetAttributeNames(), to take particle attributes as
// block attributes. This function detects and corrects such problems, if any.
// Invalid block attributes are removed, possibly re-considered particle ones.
void vtkEnzoReaderInternal::CheckAttributeNames()
{
  // number of cells of the reference block
  vtkEnzoReaderBlock& theBlock = this->Blocks[this->ReferenceBlock];
  int numCells = theBlock.BlockCellDimensions[0] * theBlock.BlockCellDimensions[1] *
    theBlock.BlockCellDimensions[2];

  // number of particles of the reference block, if any
  //  std::cout << "Reference Block: " << this->ReferenceBlock << std::endl;
  //  std::cout << "BlockIdx: " << this->ReferenceBlock - 1 << std::endl;
  //  std::cout.flush();

  vtkPolyData* polyData = vtkPolyData::New();
  //  this->TheReader->GetParticles( this->ReferenceBlock - 1, polyData, 0, 0 );

  int numbPnts = polyData->GetNumberOfPoints();
  polyData->Delete();
  polyData = nullptr;

  // block attributes to be removed and / or exported
  std::vector<std::string> toRemove;
  std::vector<std::string> toExport;
  toRemove.clear();
  toExport.clear();

  // determine to-be-removed and to-be-exported block attributes
  int i;
  int blockAttrs = static_cast<int>(this->BlockAttributeNames.size());
  for (i = 0; i < blockAttrs; i++)
  {
    // the actual number of tuples of a block attribute loaded from the
    // file for the reference block
    int numTupls = 0;
    if (this->DataArray != nullptr)
    {
      numTupls = this->DataArray->GetNumberOfTuples();
      this->ReleaseDataArray();
    }
    else
    {
      numTupls = 0;
    }

    // compare the three numbers
    if (numTupls != numCells)
    {
      if (numTupls == numbPnts)
      {
        toExport.push_back(this->BlockAttributeNames[i]);
      }
      else
      {
        toRemove.push_back(this->BlockAttributeNames[i]);
      }
    }
  }

  int nRemoves = static_cast<int>(toRemove.size());
  int nExports = static_cast<int>(toExport.size());

  // remove block attributes
  for (i = 0; i < nRemoves; i++)
  {
    for (std::vector<std::string>::iterator stringIt = this->BlockAttributeNames.begin();
         stringIt != this->BlockAttributeNames.end(); ++stringIt)
    {
      if ((*stringIt) == toRemove[i])
      {
        this->BlockAttributeNames.erase(stringIt);
        break;
      }
    }
  }

  // export attributes from blocks to particles
  for (i = 0; i < nExports; i++)
  {
    for (std::vector<std::string>::iterator stringIt = this->BlockAttributeNames.begin();
         stringIt != this->BlockAttributeNames.end(); ++stringIt)
    {
      if ((*stringIt) == toExport[i])
      {
        this->ParticleAttributeNames.push_back(*stringIt);
        this->BlockAttributeNames.erase(stringIt);
        break;
      }
    }
  }

  toRemove.clear();
  toExport.clear();
}

//------------------------------------------------------------------------------
// get the meta data
void vtkEnzoReaderInternal::ReadMetaData()
{
  // Check to see if we have read it
  if (this->NumberOfBlocks > 0)
  {
    return;
  }

  // get the general parameters (number of dimensions)
  this->ReadGeneralParameters();

  // obtain the block structures
  this->ReadBlockStructures();

  // determine the bounding box of the root block
  this->DetermineRootBoundingBox();

  // get the parent-wise and level-based bounding Ids of each block in a
  // top-down manner
  int blocks = static_cast<int>(this->Blocks.size());
  for (int i = 1; i < blocks; i++)
  {
    this->Blocks[i].GetParentWiseIds(this->Blocks);
    this->Blocks[i].GetLevelBasedIds(this->Blocks);
  }

  // locate the block that contains the fewest cells (either with or without
  // particles) and collect the attribute names

  this->GetAttributeNames();

  // verify the initial set of attribute names
  this->CheckAttributeNames();
}