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

  Program:   Visualization Toolkit
  Module:    vtkCONVERGECFDReader.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 "vtkCONVERGECFDReader.h"

#include "vtkBuffer.h"
#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkCommand.h"
#include "vtkDataArraySelection.h"
#include "vtkDirectory.h"
#include "vtkFloatArray.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkMultiBlockDataSet.h"
#include "vtkNew.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkPolyData.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkUnstructuredGrid.h"

#include <vtksys/RegularExpression.hxx>

#include <array>
#include <map>
#include <set>
#include <sstream>

#define H5_USE_16_API
#include "vtk_hdf5.h"

vtkStandardNewMacro(vtkCONVERGECFDReader);

//----------------------------------------------------------------------------
class vtkCONVERGECFDReader::vtkInternal
{
public:
  std::vector<std::string> CellDataScalarVariables;
  std::vector<std::string> CellDataVectorVariables;
  std::vector<std::string> ParcelDataScalarVariables;
  std::vector<std::string> ParcelDataVectorVariables;
};

//----------------------------------------------------------------------------
vtkCONVERGECFDReader::vtkCONVERGECFDReader()
  : FileName(nullptr)
  , Internal(new vtkCONVERGECFDReader::vtkInternal())
{
  this->SetNumberOfInputPorts(0);
  this->SetNumberOfOutputPorts(1);

  this->CellDataArraySelection->AddObserver(
    vtkCommand::ModifiedEvent, this, &vtkCONVERGECFDReader::Modified);
  this->ParcelDataArraySelection->AddObserver(
    vtkCommand::ModifiedEvent, this, &vtkCONVERGECFDReader::Modified);
}

//----------------------------------------------------------------------------
vtkCONVERGECFDReader::~vtkCONVERGECFDReader()
{
  delete[] this->FileName;
  this->FileName = nullptr;
  delete this->Internal;
}

//----------------------------------------------------------------------------
void vtkCONVERGECFDReader::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);
}

//----------------------------------------------------------------------------
namespace
{

//----------------------------------------------------------------------------
/**
 * RAII class for automatically closing H5 handles.
 */
#define DefineScopedHandle(name)                                                                   \
  class ScopedH5##name##Handle                                                                     \
  {                                                                                                \
  public:                                                                                          \
    ScopedH5##name##Handle(const ScopedH5##name##Handle& other) { this->Handle = other.Handle; }   \
    ScopedH5##name##Handle(hid_t handle)                                                           \
      : Handle(handle)                                                                             \
    {                                                                                              \
    }                                                                                              \
    virtual ~ScopedH5##name##Handle()                                                              \
    {                                                                                              \
      if (this->Handle >= 0)                                                                       \
      {                                                                                            \
        H5##name##close(this->Handle);                                                             \
      }                                                                                            \
    }                                                                                              \
                                                                                                   \
    operator hid_t() const { return this->Handle; }                                                \
                                                                                                   \
  private:                                                                                         \
    hid_t Handle;                                                                                  \
  };

// Defines ScopedH5AHandle closed with H5Aclose
DefineScopedHandle(A);

// Defines ScopedH5DHandle closed with H5Dclose
DefineScopedHandle(D);

// Defines ScopedH5FHandle closed with H5Fclose
DefineScopedHandle(F);

// Defines ScopedH5GHandle closed with H5Gclose
DefineScopedHandle(G);

// Defines ScopedH5SHandle closed with H5Sclose
DefineScopedHandle(S);

// Defines ScopedH5THandle closed with H5Tclose
DefineScopedHandle(T);

//----------------------------------------------------------------------------
/**
 *  Check existence of array defined by pathName relative to fileId.
 */
bool ArrayExists(hid_t fileId, const char* pathName)
{
  return (H5Lexists(fileId, pathName, H5P_DEFAULT) > 0);
}

//----------------------------------------------------------------------------
/**
 *  Get length of array defined by pathName relative to fileId.
 */
hsize_t GetDataLength(hid_t fileId, const char* pathName)
{
  ScopedH5DHandle arrayId = H5Dopen(fileId, pathName);
  if (arrayId < 0)
  {
    vtkGenericWarningMacro("No array named " << pathName << " available");
    return 0;
  }

  ScopedH5DHandle dataspace = H5Dget_space(arrayId);
  if (H5Sget_simple_extent_ndims(dataspace) != 1)
  {
    vtkGenericWarningMacro("Array " << pathName << " dimensionality is not 1");
    return 0;
  }

  hsize_t length = 0;
  int numDimensions = H5Sget_simple_extent_dims(dataspace, &length, nullptr);
  if (numDimensions < 0)
  {
    vtkGenericWarningMacro("Failed to get length of array");
    return 0;
  }

  return length;
}

//----------------------------------------------------------------------------
/**
 *  Read a typed array and into an array passed in by the caller. Checks that the
 *  number of elements in the array specified by fileId and pathName matches the length
 *  argument n.
 *
 *  Returns true if reading succeeded, false otherwise.
 */
template <typename T>
bool ReadArray(hid_t fileId, const char* pathName, T* data, hsize_t n)
{
  ScopedH5DHandle arrayId = H5Dopen(fileId, pathName);
  if (arrayId < 0)
  {
    vtkGenericWarningMacro("No array named " << pathName << " available");
    return false;
  }

  ScopedH5DHandle rawType = H5Dget_type(arrayId);
  ScopedH5THandle dataType = H5Tget_native_type(rawType, H5T_DIR_ASCEND);
  ScopedH5DHandle dataspace = H5Dget_space(arrayId);
  if (H5Sget_simple_extent_ndims(dataspace) != 1)
  {
    vtkGenericWarningMacro("Array " << pathName << " dimensionality is not 1");
    return false;
  }

  hsize_t length = 0;
  int numDims = H5Sget_simple_extent_dims(dataspace, &length, nullptr);
  if (numDims < 0)
  {
    vtkGenericWarningMacro("Failed to get length of array");
    return false;
  }

  if (n != length)
  {
    vtkGenericWarningMacro(
      "Size of array passed in does not match length of array. Skipping array.");
    return false;
  }

  if (H5Dread(arrayId, dataType, H5S_ALL, H5S_ALL, H5P_DEFAULT, data) < 0)
  {
    vtkGenericWarningMacro("Could not read " << pathName);
    return false;
  }

  return true;
}

//----------------------------------------------------------------------------
/**
 * Get an array of strings from a table defined by pathName relative to fileId.
 * Strings are returned in the vector of strings parameter that was passed in.
 */
bool ReadStrings(hid_t fileId, const char* path, std::vector<std::string>& strings)
{
  ScopedH5DHandle stringsId = H5Dopen(fileId, path);
  if (stringsId < 0)
  {
    vtkGenericWarningMacro("Could not read " << path);
    return false;
  }

  ScopedH5THandle filetype = H5Dget_type(stringsId);
  size_t sdim = H5Tget_size(filetype);
  sdim++; /* Make room for null terminator */

  ScopedH5SHandle space = H5Dget_space(stringsId);
  hsize_t dim;
  int ndims = H5Sget_simple_extent_dims(space, &dim, nullptr);
  if (ndims != 1)
  {
    vtkGenericWarningMacro("String array dimension not 1");
    return false;
  }

  char** rdata = new char*[dim];
  rdata[0] = new char[dim * sdim];
  for (hsize_t i = 1; i < dim; ++i)
  {
    rdata[i] = rdata[0] + i * sdim;
  }

  ScopedH5THandle memtype = H5Tcopy(H5T_C_S1);
  H5Tset_size(memtype, sdim);
  if (H5Dread(stringsId, memtype, H5S_ALL, H5S_ALL, H5P_DEFAULT, rdata[0]) < 0)
  {
    vtkGenericWarningMacro("Could not read " << path);
    return false;
  }

  strings.clear();
  for (hsize_t i = 0; i < dim; ++i)
  {
    strings.emplace_back(std::string(rdata[i]));
  }

  delete[] rdata[0];
  delete[] rdata;

  return true;
}

void SplitScalarAndVectorVariables(
  std::vector<std::string>& allVariables, std::vector<std::string>& vectorVariables)
{
  vectorVariables.clear();
  for (const auto& varName : allVariables)
  {
    // See if variable is an array.
    if (varName.find_last_of('_') == varName.size() - 2)
    {
      const char componentName = varName[varName.size() - 1];
      if (componentName == 'X')
      {
        // Check that components Y and Z exist as well
        std::string baseName = varName.substr(0, varName.size() - 2);
        if (std::find(allVariables.begin(), allVariables.end(), baseName + "_Y") !=
            allVariables.end() &&
          std::find(allVariables.begin(), allVariables.end(), baseName + "_Z") !=
            allVariables.end())
        {
          vectorVariables.emplace_back(baseName);
        }
      }
    }
  }

  // Now remove the vector variables from all variables. At the end, allVariables will contain
  // only scalar array names.
  for (const auto& varName : vectorVariables)
  {
    allVariables.erase(std::find(allVariables.begin(), allVariables.end(), varName + "_X"));
    allVariables.erase(std::find(allVariables.begin(), allVariables.end(), varName + "_Y"));
    allVariables.erase(std::find(allVariables.begin(), allVariables.end(), varName + "_Z"));
  }
}

} // anonymous namespace

//----------------------------------------------------------------------------
int vtkCONVERGECFDReader::RequestInformation(
  vtkInformation*, vtkInformationVector**, vtkInformationVector* outInfos)
{
  if (!this->FileName || this->FileName[0] == '\0')
  {
    return 1;
  }

  ScopedH5FHandle fileId = H5Fopen(this->FileName, H5F_ACC_RDONLY, H5P_DEFAULT);
  if (fileId < 0)
  {
    vtkErrorMacro("Could not open HDF5 file '" << this->FileName << "'");
    return 0;
  }

  // Assume that the variables are the same in all streams
  ScopedH5GHandle streamId = H5Gopen(fileId, "/STREAM_00");
  if (streamId < 0)
  {
    vtkErrorMacro("Cannot open group /STREAM_00");
    return 0;
  }

  if (!ReadStrings(
        streamId, "VARIABLE_NAMES/CELL_VARIABLES", this->Internal->CellDataScalarVariables))
  {
    vtkErrorMacro("Could not read cell variable names");
    return 0;
  }

  bool defaultEnabledState = true;

  // Split cell variables into scalar and vector arrays
  SplitScalarAndVectorVariables(
    this->Internal->CellDataScalarVariables, this->Internal->CellDataVectorVariables);

  for (const auto& varName : this->Internal->CellDataScalarVariables)
  {
    if (!this->CellDataArraySelection->ArrayExists(varName.c_str()))
    {
      this->CellDataArraySelection->AddArray(varName.c_str(), defaultEnabledState);
    }
  }

  for (const auto& varName : this->Internal->CellDataVectorVariables)
  {
    if (!this->CellDataArraySelection->ArrayExists(varName.c_str()))
    {
      this->CellDataArraySelection->AddArray(varName.c_str(), defaultEnabledState);
    }
  }

  if (ArrayExists(streamId, "VARIABLE_NAMES/PARCEL_VARIABLES"))
  {
    if (!ReadStrings(
          streamId, "VARIABLE_NAMES/PARCEL_VARIABLES", this->Internal->ParcelDataScalarVariables))
    {
      vtkErrorMacro("Could not read parcel variable names");
      return 0;
    }

    // Split parcel arrays into scalar and vector variables
    SplitScalarAndVectorVariables(
      this->Internal->ParcelDataScalarVariables, this->Internal->ParcelDataVectorVariables);

    for (const auto& varName : this->Internal->ParcelDataScalarVariables)
    {
      if (!this->ParcelDataArraySelection->ArrayExists(varName.c_str()))
      {
        this->ParcelDataArraySelection->AddArray(varName.c_str(), defaultEnabledState);
      }
    }

    for (const auto& varName : this->Internal->ParcelDataVectorVariables)
    {
      // Skip X, Y, Z points
      if (varName == "PARCEL")
      {
        continue;
      }

      if (!this->ParcelDataArraySelection->ArrayExists(varName.c_str()))
      {
        this->ParcelDataArraySelection->AddArray(varName.c_str(), defaultEnabledState);
      }
    }
  }

  // Get time information
  vtkInformation* outInfo = outInfos->GetInformationObject(0);
  this->ReadTimeSteps(outInfo);

  return 1;
}

//----------------------------------------------------------------------------
int vtkCONVERGECFDReader::RequestData(
  vtkInformation*, vtkInformationVector**, vtkInformationVector* outputVector)
{
  vtkInformation* outInfo = outputVector->GetInformationObject(0);

  size_t fileIndex = this->SelectTimeStepIndex(outInfo);
  std::string fileName = this->FileNames[fileIndex];

  if (fileName.empty())
  {
    vtkErrorMacro("No file sequence found");
    return 0;
  }

  vtkMultiBlockDataSet* outputMB = vtkMultiBlockDataSet::GetData(outInfo);

  if (!outputMB)
  {
    vtkErrorMacro("No output available!");
    return 0;
  }

  ScopedH5FHandle fileId = H5Fopen(fileName.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
  if (fileId < 0)
  {
    vtkErrorMacro("Could not open HDF5 file '" << fileName << "'");
    return 0;
  }

  ScopedH5GHandle boundaryHandle = H5Gopen(fileId, "/BOUNDARIES");
  if (boundaryHandle < 0)
  {
    vtkErrorMacro("Cannot open group/BOUNDARIES");
    return 0;
  }

  // Iterate over stream groups
  hsize_t numObjs = 0;
  herr_t err = H5Gget_num_objs(fileId, &numObjs);
  if (err < 0)
  {
    vtkErrorMacro("Cannot get number of groups from file");
    return 0;
  }

  int streamCount = 0;
  do
  {
    herr_t status = 0;
    std::ostringstream streamName;
    streamName << "/STREAM_" << std::setw(2) << std::setfill('0') << streamCount;
    H5Eset_auto(nullptr, nullptr);
    status = H5Gget_objinfo(fileId, streamName.str().c_str(), false, nullptr);
    if (status < 0)
    {
      break;
    }

    // Open the group
    ScopedH5GHandle streamId = H5Gopen(fileId, streamName.str().c_str());
    if (streamId < 0)
    {
      vtkErrorMacro("Could not open stream " << streamName.str());
      break;
    }

    if (!ArrayExists(streamId, "VERTEX_COORDINATES/X"))
    {
      vtkErrorMacro("Could not find array VERTEX_COORDINATES/X");
      break;
    }

    hsize_t xCoordsLength = GetDataLength(streamId, "VERTEX_COORDINATES/X");

    // Temporary buffer for reading vector array components
    vtkNew<vtkBuffer<float>> floatBuffer;

    vtkNew<vtkFloatArray> pointArray;
    pointArray->SetNumberOfComponents(3);
    pointArray->SetNumberOfTuples(xCoordsLength);

    floatBuffer->Allocate(xCoordsLength);

    std::array<char, 3> dimensionNames = { 'X', 'Y', 'Z' };

    for (size_t c = 0; c < dimensionNames.size(); ++c)
    {
      std::stringstream name;
      name << "VERTEX_COORDINATES/" << dimensionNames[c];
      if (!ReadArray(streamId, name.str().c_str(), floatBuffer->GetBuffer(), xCoordsLength))
      {
        vtkErrorMacro(
          "No coordinate array " << name.str() << " dataset available in " << streamName.str());
        return 0;
      }

      for (hsize_t j = 0; j < xCoordsLength; ++j)
      {
        pointArray->SetTypedComponent(j, static_cast<int>(c), floatBuffer->GetBuffer()[j]);
      }
    }

    // ++++ POLYGON_OFFSET ++++
    hsize_t polygonOffsetsLength = GetDataLength(streamId, "CONNECTIVITY/POLYGON_OFFSET");
    std::vector<int> polygonOffsets(polygonOffsetsLength);
    if (!ReadArray(
          streamId, "CONNECTIVITY/POLYGON_OFFSET", polygonOffsets.data(), polygonOffsetsLength))
    {
      vtkErrorMacro("Could not read CONNECTIVITY/POLYGON_OFFSET");
      return 0;
    }

    // Reduce the number of polygons by one to make up for the fact that the POLYGON_OFFSETS
    // array is longer by one row.
    vtkIdType numPolygons = static_cast<vtkIdType>(polygonOffsets.size()) - 1;

    // ++++ POLYGON_TO_VERTEX ++++
    hsize_t polygonsLength = GetDataLength(streamId, "CONNECTIVITY/POLYGON_TO_VERTEX");
    std::vector<int> polygons(polygonsLength);
    if (!ReadArray(streamId, "CONNECTIVITY/POLYGON_TO_VERTEX", polygons.data(), polygonsLength))
    {
      vtkErrorMacro("Could not read CONNECTIVITY/POLYGON_TO_VERTEX");
      return 0;
    }

    // ++++ CONNECTED_CELLS ++++
    hsize_t connectedCellsLength = GetDataLength(streamId, "CONNECTIVITY/CONNECTED_CELLS");
    std::vector<int> connectedCells(connectedCellsLength);
    if (!ReadArray(
          streamId, "CONNECTIVITY/CONNECTED_CELLS", connectedCells.data(), connectedCellsLength))
    {
      vtkErrorMacro("Could not read CONNECTIVITY/CONNECTED_CELLS");
      return 0;
    }

    // ++++ CREATE VTK DATA SETS ++++
    vtkNew<vtkPoints> points;
    points->SetData(pointArray);

    // boundaryIdToIndex must be size of max id... ids are not guaranteed to be sequential,
    // i.e., 1, 3, 5, 30, 31, 32, 1001 so it's better to use map instead of array lookup
    std::map<int, int> boundaryIdToIndex;

    hsize_t numBoundaryNames = GetDataLength(boundaryHandle, "BOUNDARY_NAMES");
    std::vector<std::string> boundaryNames(numBoundaryNames);
    ReadStrings(boundaryHandle, "BOUNDARY_NAMES", boundaryNames);
    hsize_t numBoundaries = GetDataLength(boundaryHandle, "NUM_ELEMENTS");
    std::vector<int> boundaryNumElements(numBoundaries);
    ReadArray(
      boundaryHandle, "NUM_ELEMENTS", boundaryNumElements.data(), boundaryNumElements.size());
    std::vector<int> boundaryIds(numBoundaries);
    ReadArray(boundaryHandle, "BOUNDARY_IDS", boundaryIds.data(), boundaryIds.size());
    if (numBoundaries != numBoundaryNames)
    {
      vtkErrorMacro("Number of BOUNDARY_NAMES does not match NUM_ELEMENTS");
      return 0;
    }

    // Multiple surfaces can exist in a single file. We create a vtkPolyData for
    // each one and store them under another multiblock dataset.
    vtkNew<vtkMultiBlockDataSet> multiSurfaces;
    multiSurfaces->SetNumberOfBlocks(static_cast<unsigned int>(boundaryIds.size()));
    for (size_t i = 0; i < boundaryIds.size(); ++i)
    {
      // If boundary index 0 has boundary id == 1, index 1 of boundaryIdToIndex will be 0.
      boundaryIdToIndex[boundaryIds[i]] = static_cast<int>(i);

      vtkNew<vtkPolyData> boundarySurface;
      multiSurfaces->SetBlock(static_cast<unsigned int>(i), boundarySurface);
      multiSurfaces->GetMetaData(static_cast<unsigned int>(i))
        ->Set(vtkCompositeDataSet::NAME(), boundaryNames[i]);

      vtkNew<vtkCellArray> polys;
      polys->AllocateEstimate(boundaryNumElements[i], 4);
      boundarySurface->SetPolys(polys);
    }

    // Create maps from surface point IDs for each block
    std::vector<std::set<vtkIdType>> blocksSurfacePointIds(boundaryIds.size());
    for (int polyId = 0; polyId < numPolygons; ++polyId)
    {
      if (connectedCells[2 * polyId + 0] >= 0 && connectedCells[2 * polyId + 1] >= 0)
      {
        // Polygon is not part of a surface, so skip.
        continue;
      }

      int boundaryId = -(connectedCells[2 * polyId + 0] + 1);
      int boundaryIndex = boundaryIdToIndex[boundaryId];
      vtkIdType numCellPts =
        static_cast<vtkIdType>(polygonOffsets[polyId + 1] - polygonOffsets[polyId]);
      for (vtkIdType id = 0; id < numCellPts; ++id)
      {
        vtkIdType ptId = polygons[polygonOffsets[polyId] + id];
        blocksSurfacePointIds[boundaryIndex].insert(ptId);
      }
    }

    // Create maps from original point IDs to surface point IDs for each block
    std::vector<std::map<vtkIdType, vtkIdType>> blocksOriginalToBlockPointId(numBoundaryNames);
    for (hsize_t boundaryIndex = 0; boundaryIndex < numBoundaryNames; ++boundaryIndex)
    {
      // Create a map from original point ID in the global points list
      vtkIdType newIndex = 0;
      for (vtkIdType id : blocksSurfacePointIds[boundaryIndex])
      {
        blocksOriginalToBlockPointId[boundaryIndex][id] = newIndex++;
      }

      // Clear some memory
      blocksSurfacePointIds[boundaryIndex].clear();

      // Create localized points for this block
      vtkNew<vtkPoints> blockPoints;
      blockPoints->SetDataType(points->GetDataType());
      blockPoints->SetNumberOfPoints(newIndex);
      vtkFloatArray* toArray = vtkFloatArray::SafeDownCast(blockPoints->GetData());
      for (auto it = blocksOriginalToBlockPointId[boundaryIndex].begin();
           it != blocksOriginalToBlockPointId[boundaryIndex].end(); ++it)
      {
        vtkIdType from = it->first;
        vtkIdType to = it->second;
        float xyz[3];
        pointArray->GetTypedTuple(from, xyz);
        toArray->SetTypedTuple(to, xyz);
      }

      vtkPolyData* boundarySurface =
        vtkPolyData::SafeDownCast(multiSurfaces->GetBlock(boundaryIndex));
      boundarySurface->SetPoints(blockPoints);
    }

    // Go through polygons again and add them to the polydata blocks
    vtkIdType numSurfacePolys = 0;
    for (int polyId = 0; polyId < numPolygons; ++polyId)
    {
      if (connectedCells[2 * polyId + 0] >= 0 && connectedCells[2 * polyId + 1] >= 0)
      {
        // Polygon is not part of a surface, so skip.
        continue;
      }

      numSurfacePolys++;

      int boundaryId = -(connectedCells[2 * polyId + 0] + 1);
      int boundaryIndex = boundaryIdToIndex[boundaryId];
      vtkPolyData* polyData = vtkPolyData::SafeDownCast(multiSurfaces->GetBlock(boundaryIndex));
      vtkIdType numCellPts =
        static_cast<vtkIdType>(polygonOffsets[polyId + 1] - polygonOffsets[polyId]);
      std::vector<vtkIdType> ptIds(numCellPts);
      for (vtkIdType id = 0; id < numCellPts; ++id)
      {
        vtkIdType ptId = polygons[polygonOffsets[polyId] + id];
        ptIds[id] = blocksOriginalToBlockPointId[boundaryIndex][ptId];
      }
      polyData->GetPolys()->InsertNextCell(numCellPts, &ptIds[0]);
    }

    // Clear some memory
    blocksOriginalToBlockPointId.clear();

    // Create a map from cell to polygons
    std::map<int, std::set<int>> cellToPoly;

    // Create a map from polygon to the volumetric cell to which it is attached
    std::vector<int> polyToCell(numSurfacePolys);

    // Create a map from polygon to boundary
    std::vector<int> polyToBoundary(numSurfacePolys);

    vtkIdType surfacePolyCount = 0;
    for (int polyId = 0; polyId < numPolygons; ++polyId)
    {
      int cell0 = connectedCells[2 * polyId + 0];
      int cell1 = connectedCells[2 * polyId + 1];
      if (cell0 >= 0)
      {
        // Add polyId to cell 0's list of polygons
        cellToPoly[cell0].insert(polyId);
      }
      if (cell1 >= 0)
      {
        // Add polyId to cell 1's list of polygons
        cellToPoly[cell1].insert(polyId);
      }

      if (cell0 < 0 || cell1 < 0)
      {
        assert(polyToBoundary.size() > static_cast<size_t>(surfacePolyCount));
        polyToBoundary[surfacePolyCount] = cell0 >= 0 ? -(cell1 + 1) : -(cell0 + 1);
        polyToCell[surfacePolyCount] = cell0 >= 0 ? cell0 : cell1;
        surfacePolyCount++;
      }
    }

    // Create a new unstructured grid
    vtkNew<vtkUnstructuredGrid> ugrid;
    ugrid->SetPoints(points);

    // Create polyhedra from their faces
    vtkNew<vtkIdList> faces;
    for (const auto& cellIdAndPolys : cellToPoly)
    {
      faces->Reset();
      // Number of faces
      faces->InsertNextId(static_cast<vtkIdType>(cellIdAndPolys.second.size()));
      for (auto polyId : cellIdAndPolys.second)
      {
        // Get polygon
        int numPts = polygonOffsets[polyId + 1] - polygonOffsets[polyId];

        // Number of points in face
        faces->InsertNextId(numPts);
        for (int i = 0; i < numPts; ++i)
        {
          // Polygon vertex
          faces->InsertNextId(polygons[polygonOffsets[polyId] + i]);
        }
      }

      ugrid->InsertNextCell(VTK_POLYHEDRON, faces);
    }

    // ++++ CELL DATA ++++
    for (int i = 0; i < this->CellDataArraySelection->GetNumberOfArrays(); ++i)
    {
      std::string varName(this->CellDataArraySelection->GetArrayName(i));
      if (this->GetCellDataArraySelection()->ArrayIsEnabled(varName.c_str()) == 0)
      {
        continue;
      }

      auto begin = this->Internal->CellDataVectorVariables.begin();
      auto end = this->Internal->CellDataVectorVariables.end();
      bool isVector = std::find(begin, end, varName) != end;

      // This would be a lot simpler using a vtkSOADataArrayTemplate<float>, but
      // until GetVoidPointer() is removed from more of the VTK code base, we
      // will use a vtkFloatArray.
      vtkNew<vtkFloatArray> dataArray;
      bool success = true;
      if (isVector)
      {
        std::string rootPath("CELL_CENTER_DATA/");
        rootPath += varName;
        std::string pathX = rootPath + "_X";
        std::string pathY = rootPath + "_Y";
        std::string pathZ = rootPath + "_Z";

        hsize_t length = GetDataLength(streamId, pathX.c_str());
        dataArray->SetNumberOfComponents(3);
        dataArray->SetNumberOfTuples(length);
        dataArray->SetName(varName.c_str());

        if (floatBuffer->GetSize() != static_cast<vtkIdType>(length))
        {
          floatBuffer->Allocate(length);
        }

        success = success && ReadArray(streamId, pathX.c_str(), floatBuffer->GetBuffer(), length);
        for (hsize_t j = 0; j < length; ++j)
        {
          dataArray->SetTypedComponent(j, 0, floatBuffer->GetBuffer()[j]);
        }
        success = success && ReadArray(streamId, pathY.c_str(), floatBuffer->GetBuffer(), length);
        for (hsize_t j = 0; j < length; ++j)
        {
          dataArray->SetTypedComponent(j, 1, floatBuffer->GetBuffer()[j]);
        }
        success = success && ReadArray(streamId, pathZ.c_str(), floatBuffer->GetBuffer(), length);
        for (hsize_t j = 0; j < length; ++j)
        {
          dataArray->SetTypedComponent(j, 2, floatBuffer->GetBuffer()[j]);
        }
      }
      else
      {
        std::string path("CELL_CENTER_DATA/");
        path += varName;

        hsize_t length = GetDataLength(streamId, path.c_str());
        dataArray->SetNumberOfComponents(1);
        dataArray->SetNumberOfTuples(length);
        dataArray->SetName(varName.c_str());
        success = success && ReadArray(streamId, path.c_str(), dataArray->GetPointer(0), length);
      }

      if (success)
      {
        ugrid->GetCellData()->AddArray(dataArray);

        // Now pull out the values needed for the surface geometry
        for (int boundaryIndex = 0;
             boundaryIndex < static_cast<int>(multiSurfaces->GetNumberOfBlocks()); ++boundaryIndex)
        {
          vtkPolyData* boundarySurface =
            vtkPolyData::SafeDownCast(multiSurfaces->GetBlock(boundaryIndex));
          int numBoundaryPolys = boundarySurface->GetNumberOfCells();
          vtkNew<vtkFloatArray> surfaceDataArray;
          surfaceDataArray->SetNumberOfComponents(dataArray->GetNumberOfComponents());
          surfaceDataArray->SetNumberOfTuples(numBoundaryPolys);
          surfaceDataArray->SetName(varName.c_str());
          vtkIdType localDataCount = 0;
          for (vtkIdType id = 0; id < numSurfacePolys; ++id)
          {
            assert(polyToBoundary.size() > static_cast<size_t>(id));
            if (boundaryIdToIndex.find(polyToBoundary[id]) == boundaryIdToIndex.end())
            {
              vtkErrorMacro(
                "polyToBoundary[id] is not found within boundaryIdToIndex" << polyToBoundary[id]);
              return 0;
            }
            const int polyBoundaryIndex = boundaryIdToIndex[polyToBoundary[id]];
            if (polyBoundaryIndex != boundaryIndex)
            {
              continue;
            }
            for (int c = 0; c < surfaceDataArray->GetNumberOfComponents(); ++c)
            {
              assert(polyToCell.size() > static_cast<size_t>(id));
              surfaceDataArray->SetTypedComponent(
                localDataCount, c, dataArray->GetTypedComponent(polyToCell[id], c));
            }
            ++localDataCount;
          }
          boundarySurface->GetCellData()->AddArray(surfaceDataArray);
        }
      }
    }

    // Set the mesh and surface in the multiblock output
    vtkNew<vtkMultiBlockDataSet> streamMB;
    streamMB->SetNumberOfBlocks(2);
    streamMB->SetBlock(0, ugrid);
    streamMB->GetMetaData(0u)->Set(vtkCompositeDataSet::NAME(), "Mesh");
    streamMB->SetBlock(1, multiSurfaces);
    streamMB->GetMetaData(1u)->Set(vtkCompositeDataSet::NAME(), "Surface");

    // ++++ PARCEL DATA ++++
    bool parcelExists = ArrayExists(streamId, "PARCEL_DATA/PARCEL_X");
    if (parcelExists)
    {
      // Read parcel point locations
      hsize_t parcelLength = GetDataLength(streamId, "PARCEL_DATA/PARCEL_X");

      vtkNew<vtkFloatArray> parcelPointArray;
      parcelPointArray->SetNumberOfComponents(3);
      parcelPointArray->SetNumberOfTuples(parcelLength);

      floatBuffer->Allocate(parcelLength);

      for (size_t c = 0; c < dimensionNames.size(); ++c)
      {
        std::stringstream name;
        name << "PARCEL_DATA/PARCEL_" << dimensionNames[c];
        if (!ReadArray(streamId, name.str().c_str(), floatBuffer->GetBuffer(), parcelLength))
        {
          vtkErrorMacro("No parcel coordinate array " << name.str() << " dataset available in "
                                                      << streamName.str());
          return 0;
        }

        for (hsize_t j = 0; j < parcelLength; ++j)
        {
          parcelPointArray->SetTypedComponent(j, static_cast<int>(c), floatBuffer->GetBuffer()[j]);
        }
      }

      vtkNew<vtkPoints> parcelPoints;
      parcelPoints->SetData(parcelPointArray);

      vtkNew<vtkPolyData> parcels;
      parcels->SetPoints(parcelPoints);

      // Create a vertex for each parcel point
      vtkNew<vtkCellArray> parcelCells;
      parcelCells->AllocateExact(parcelLength, 1);
      for (vtkIdType id = 0; id < static_cast<vtkIdType>(parcelLength); ++id)
      {
        parcelCells->InsertNextCell(1, &id);
      }
      parcels->SetVerts(parcelCells);

      // Read parcel data arrays
      for (int i = 0; i < this->ParcelDataArraySelection->GetNumberOfArrays(); ++i)
      {
        std::string varName(this->ParcelDataArraySelection->GetArrayName(i));
        if (varName == "PARCEL_X" || varName == "PARCEL_Y" || varName == "PARCEL_Z" ||
          this->ParcelDataArraySelection->ArrayIsEnabled(varName.c_str()) == 0)
        {
          continue;
        }

        auto begin = this->Internal->ParcelDataVectorVariables.begin();
        auto end = this->Internal->ParcelDataVectorVariables.end();
        bool isVector = std::find(begin, end, varName) != end;

        // This would be a lot simpler using a vtkSOADataArrayTemplate<float>, but
        // until GetVoidPointer() is removed from more of the VTK code base, we
        // will use a vtkFloatArray.
        vtkNew<vtkFloatArray> dataArray;
        bool success = true;
        if (isVector)
        {
          std::string rootPath("PARCEL_DATA/");
          rootPath += varName;
          std::string pathX = rootPath + "_X";
          std::string pathY = rootPath + "_Y";
          std::string pathZ = rootPath + "_Z";

          // hsize_t length = GetDataLength(streamId, pathX.c_str());
          dataArray->SetNumberOfComponents(3);
          dataArray->SetNumberOfTuples(parcelLength);
          dataArray->SetName(varName.c_str());

          if (static_cast<hsize_t>(floatBuffer->GetSize()) != parcelLength)
          {
            floatBuffer->Allocate(parcelLength);
          }
          success =
            success && ReadArray(streamId, pathX.c_str(), floatBuffer->GetBuffer(), parcelLength);
          for (hsize_t j = 0; j < parcelLength; ++j)
          {
            dataArray->SetTypedComponent(j, 0, floatBuffer->GetBuffer()[j]);
          }
          success =
            success && ReadArray(streamId, pathY.c_str(), floatBuffer->GetBuffer(), parcelLength);
          for (hsize_t j = 0; j < parcelLength; ++j)
          {
            dataArray->SetTypedComponent(j, 1, floatBuffer->GetBuffer()[j]);
          }
          success =
            success && ReadArray(streamId, pathZ.c_str(), floatBuffer->GetBuffer(), parcelLength);
          for (hsize_t j = 0; j < parcelLength; ++j)
          {
            dataArray->SetTypedComponent(j, 2, floatBuffer->GetBuffer()[j]);
          }
        }
        else
        {
          std::string path("PARCEL_DATA/");
          path += varName;

          dataArray->SetNumberOfComponents(1);
          dataArray->SetNumberOfTuples(parcelLength);
          dataArray->SetName(varName.c_str());
          success =
            success && ReadArray(streamId, path.c_str(), dataArray->GetPointer(0), parcelLength);
        }

        if (success)
        {
          parcels->GetPointData()->AddArray(dataArray);
        }
      }

      streamMB->SetNumberOfBlocks(3);
      streamMB->SetBlock(2, parcels);
      streamMB->GetMetaData(2u)->Set(vtkCompositeDataSet::NAME(), "Parcels");
    }

    std::stringstream streamss;
    streamss << "STREAM_" << std::setw(2) << std::setfill('0') << streamCount;
    outputMB->SetNumberOfBlocks(streamCount + 1);
    outputMB->GetMetaData(streamCount)->Set(vtkCompositeDataSet::NAME(), streamss.str().c_str());
    outputMB->SetBlock(streamCount, streamMB);

    streamCount++;
  } while (true);

  // Everything succeeded
  return 1;
}

//----------------------------------------------------------------------------
void vtkCONVERGECFDReader::ReadTimeSteps(vtkInformation* outInfo)
{
  if (!this->FileName || this->FileName[0] == '\0')
  {
    return;
  }

  // Scan for other files with the same naming pattern in the same directory.
  // We are looking for files with the following convention
  //
  // <file><zero-padded index>_<time>.h5
  //
  // We load each file and extract the time from within.
  this->FileNames.clear();
  std::string originalFile = this->FileName;
  std::string path;
  std::string baseName;
  std::string::size_type dirseppos = originalFile.find_last_of("/\\");
  if (dirseppos == std::string::npos)
  {
    path = "./";
    baseName = originalFile;
  }
  else
  {
    path = originalFile.substr(0, dirseppos + 1);
    baseName = originalFile.substr(dirseppos + 1);
  }

  vtksys::RegularExpression regEx("^([^0-9]*)([0-9]*)[_]?(.*).h5$");
  if (!regEx.find(baseName))
  {
    this->FileNames.emplace_back(originalFile);
    return;
  }

  std::string prefix = regEx.match(1);
  std::string indexString = regEx.match(2);

  vtkNew<vtkDirectory> dir;
  if (!dir->Open(path.c_str()))
  {
    vtkWarningMacro(<< "Could not open directory " << originalFile.c_str()
                    << " is supposed to be from (" << path.c_str() << ")");
    this->FileNames.emplace_back(originalFile);
    return;
  }

  for (vtkIdType i = 0; i < dir->GetNumberOfFiles(); i++)
  {
    const char* file = dir->GetFile(i);
    if (!regEx.find(file))
    {
      continue;
    }
    if (regEx.match(1) != prefix)
    {
      continue;
    }
    this->FileNames.emplace_back(path + file);
  }

  std::vector<double> times;
  double timeRange[2] = { VTK_DOUBLE_MAX, VTK_DOUBLE_MIN };
  for (const auto& file : this->FileNames)
  {
    double time = 0.0;
    bool timeRead = this->ReadOutputTime(file, time);
    if (timeRead)
    {
      timeRange[0] = std::min(timeRange[0], time);
      timeRange[1] = std::max(timeRange[1], time);
      times.emplace_back(time);
    }
  }

  if (!times.empty())
  {
    outInfo->Set(vtkStreamingDemandDrivenPipeline::TIME_RANGE(), timeRange, 2);
    outInfo->Set(
      vtkStreamingDemandDrivenPipeline::TIME_STEPS(), &times[0], static_cast<int>(times.size()));
  }
}

//----------------------------------------------------------------------------
bool vtkCONVERGECFDReader::ReadOutputTime(const std::string& filePath, double& time)
{
  if (filePath[0] == '\0')
  {
    return false;
  }

  ScopedH5FHandle fileId = H5Fopen(filePath.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
  if (fileId < 0)
  {
    return false;
  }

  if (H5Aexists(fileId, "OUTPUT_TIME"))
  {
    ScopedH5AHandle outputTimeId =
      H5Aopen_by_name(fileId, ".", "OUTPUT_TIME", H5P_DEFAULT, H5P_DEFAULT);
    ScopedH5THandle rawType = H5Aget_type(outputTimeId);
    ScopedH5THandle dataType = H5Tget_native_type(rawType, H5T_DIR_ASCEND);

    double outputTime = 0.0;
    if (H5Aread(outputTimeId, dataType, &outputTime) >= 0)
    {
      time = outputTime;
      return true;
    }
  }

  return false;
}

//----------------------------------------------------------------------------
size_t vtkCONVERGECFDReader::SelectTimeStepIndex(vtkInformation* info)
{
  if (!info->Has(vtkStreamingDemandDrivenPipeline::TIME_STEPS()) ||
    !info->Has(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP()))
  {
    return 0;
  }

  double* times = info->Get(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
  int nTimes = info->Length(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
  double t = info->Get(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP());

  double resultDiff = VTK_DOUBLE_MAX;
  size_t result = 0;
  for (int i = 0; i < nTimes; ++i)
  {
    double diff = std::fabs(times[i] - t);
    if (diff < resultDiff)
    {
      resultDiff = diff;
      result = static_cast<size_t>(i);
    }
  }

  return result;
}

//----------------------------------------------------------------------------
int vtkCONVERGECFDReader::CanReadFile(const char* fname)
{
  if (H5Fis_hdf5(fname) == 0)
  {
    return 0;
  }

  ScopedH5FHandle fileId = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
  if (fileId < 0)
  {
    return 0;
  }

  // Require a /BOUNDARIES group and at least one STREAM_00 group
  if (H5Lexists(fileId, "/BOUNDARIES", H5P_DEFAULT) == 0 ||
    H5Lexists(fileId, "/STREAM_00", H5P_DEFAULT) == 0)
  {
    return 0;
  }

  // Everything succeeded
  return 1;
}
