text: /*=========================================================================
  
   Program:   Visualization Toolkit
   Module:    vtkPGenericIOMultiBlockReader.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 "vtkPGenericIOMultiBlockReader.h"
  
  #include "vtkCallbackCommand.h"
  #include "vtkCellArray.h"
  #include "vtkCellType.h"
  #include "vtkCompositeDataPipeline.h"
  #include "vtkDataArraySelection.h"
  #include "vtkGenericIOUtilities.h"
  #include "vtkIdList.h"
  #include "vtkInformation.h"
  #include "vtkInformationKey.h"
  #include "vtkInformationDoubleKey.h"
  #include "vtkInformationVector.h"
  #include "vtkMultiProcessController.h"
  #include "vtkMultiBlockDataSet.h"
  #include "vtkNew.h"
  #include "vtkObjectFactory.h"
  #include "vtkPointData.h"
  #include "vtkSmartPointer.h"
  #include "vtkStringArray.h"
  #include "vtkType.h"
  #include "vtkUnstructuredGrid.h"
  #include "vtkDoubleArray.h"
  #include "vtkTypeUInt64Array.h"
  #include "vtkFieldData.h"
  
  #include "GenericIOReader.h"
  #include "GenericIOUtilities.h"
  
  #include <map>
  #include <set>
  #include <stdexcept>
  #include <utility>
  #include <vector>
  
  #include <unistd.h>
  
  // Uncomment the line below to get debugging information
  //#define DEBUG
  
  namespace {
  struct block_t
  {
    vtkTypeUInt64 GlobalId;
    vtkTypeUInt64 NumberOfElements;
    double bounds[6];
    uint64_t coords[3];
    int ProcessId;
  
    std::map< std::string, void* > RawCache;
    std::map< std::string, bool > VariableStatus;
  };
  
  struct most_of_block_t {
    vtkTypeUInt64 GlobalId;
    vtkTypeUInt64 NumberOfElements;
    double bounds[6];
    uint64_t coords[3];
    int ProcessId;
  };
  
  }
  
  //------------------------------------------------------------------------------
  class vtkPGenericIOMultiBlockReader::vtkGenericIOMultiBlockMetaData
  {
  public:
    int NumberOfBlocks; // Total number of blocks across all ranks
    int TotalNumberOfElements; // Total number of points across all blocks
    std::map < std::string, gio::VariableInfo > VariableInformation;
    std::map < std::string, int > VariableGenericIOType;
    std::map < int, block_t > Blocks;
  
  
    /**
     * @brief Metadata constructor.
     */
    vtkGenericIOMultiBlockMetaData(){}
  
    /**
     * @brief Destructor
     */
    ~vtkGenericIOMultiBlockMetaData() { this->Clear();}
  
    /**
     * @brief Checks if the supplied rank should load data.
     * @param r the rank in query.
     * @return status true or false.
     */
    bool HasBlock( const int blockId )
    {
      return (this->Blocks.find(blockId) != this->Blocks.end());
    }
  
    bool BlockIsOnMyProcess(const int blockId, const int myProcessId)
    {
      return this->Blocks.find(blockId) != this->Blocks.end() &&
          this->Blocks[blockId].ProcessId == myProcessId;
    }
  
    /**
     * @brief Get the raw MPI communicator from a Multi-process controller.
     * @param controller the multi-process controller
     * /
    void InitCommunicator(vtkMultiProcessController *controller)
    {
      assert("pre: controller is NULL!" && (controller != NULL) );
      this->MPICommunicator =
          vtkGenericIOUtilities::GetMPICommunicator(controller);
    } */
  
    /**
     * @brief Performs a quick sanity on the metadata
     * @return status false if the metadata is somehow corrupted.
     */
    bool SanityCheck()
    {
      return  (this->VariableInformation.size() == this->VariableGenericIOType.size());
    }
  
    /**
     * @brief Checks if a variable exists
     * @param varName the name of the variable in query
     * @return status true or false
     */
    bool HasVariable(const std::string &varName)
    {
      return (this->VariableInformation.find(varName) !=
        this->VariableInformation.end());
    }
  
    /**
     * @brief Clears the metadata
     */
    void Clear()
    {
      this->NumberOfBlocks = 0;
      this->VariableGenericIOType.clear();
      this->VariableInformation.clear();
  
      for (std::map< int, block_t >::iterator blockItr = this->Blocks.begin();
            blockItr != this->Blocks.end(); ++blockItr)
        {
        std::map< std::string, void* >::iterator varItr = blockItr->second.RawCache.begin();
        for ( ; varItr != blockItr->second.RawCache.end(); ++varItr)
          {
          if (varItr->second != NULL )
            {
            delete [] static_cast< char* >(varItr->second);
            }
          }
        }
      this->Blocks.clear();
  
    }
  };
  
  vtkStandardNewMacro(vtkPGenericIOMultiBlockReader)
  //------------------------------------------------------------------------------
  vtkPGenericIOMultiBlockReader::vtkPGenericIOMultiBlockReader()
  {
    this->SetNumberOfInputPorts(0);
    this->SetNumberOfOutputPorts(1);
  
    this->Controller        = vtkMultiProcessController::GetGlobalController();
    this->Reader            = NULL;
    this->FileName          = NULL;
    this->XAxisVariableName = NULL;
    this->YAxisVariableName = NULL;
    this->ZAxisVariableName = NULL;
    this->HaloIdVariableName = NULL;
    this->GenericIOType     = IOTYPEMPI;
    this->BlockAssignment   = ROUND_ROBIN;
    this->BuildMetaData     = false;
  
    this->SetXAxisVariableName("x");
    this->SetYAxisVariableName("y");
    this->SetZAxisVariableName("z");
  
    this->MetaData = new vtkGenericIOMultiBlockMetaData();
    //this->MetaData->InitCommunicator(this->Controller);
  
    this->ArrayList = vtkStringArray::New();
    this->PointDataArraySelection = vtkDataArraySelection::New();
    this->HaloList = vtkIdList::New();
    this->SelectionObserver = vtkCallbackCommand::New();
    this->SelectionObserver->SetCallback(
          &vtkPGenericIOMultiBlockReader::SelectionModifiedCallback);
    this->SelectionObserver->SetClientData( this );
    this->PointDataArraySelection->AddObserver(
          vtkCommand::ModifiedEvent,this->SelectionObserver);
  }
  
  //------------------------------------------------------------------------------
  vtkPGenericIOMultiBlockReader::~vtkPGenericIOMultiBlockReader()
  {
    if (this->Reader != NULL)
      {
      this->Reader->Close();
      delete this->Reader;
      }
  
    vtkGenericIOUtilities::SafeDeleteString(this->FileName);
    vtkGenericIOUtilities::SafeDeleteString(this->XAxisVariableName);
    vtkGenericIOUtilities::SafeDeleteString(this->YAxisVariableName);
    vtkGenericIOUtilities::SafeDeleteString(this->ZAxisVariableName);
    vtkGenericIOUtilities::SafeDeleteString(this->HaloIdVariableName);
  
    if (this->MetaData != NULL)
      {
      delete this->MetaData;
      }
  
    this->ArrayList->Delete();
    this->PointDataArraySelection->RemoveObserver(this->SelectionObserver);
    this->HaloList->Delete();
    this->SelectionObserver->Delete();
    this->PointDataArraySelection->Delete();
  
    this->Controller = NULL;
  }
  
  //------------------------------------------------------------------------------
  void vtkPGenericIOMultiBlockReader::PrintSelf(std::ostream &os, vtkIndent indent)
  {
    this->Superclass::PrintSelf(os,indent);
    os << indent << "FileName: " << this->FileName << endl;
    os << indent << "x-axis: " << this->XAxisVariableName << endl;
    os << indent << "y-axis: " << this->YAxisVariableName << endl;
    os << indent << "z-axis: " << this->ZAxisVariableName << endl;
    os << indent << "GenericIOType: " << this->GenericIOType << endl;
    os << indent << "BlockAssignment: " << this->BlockAssignment << endl;
    os << indent << "ArrayList: " << endl;
    this->ArrayList->PrintSelf(os,indent.GetNextIndent());
    os << indent << "PointDataSelection: " << endl;
    this->PointDataArraySelection->PrintSelf(os,indent.GetNextIndent());
    if( Controller != NULL )
      {
      os << indent << "Controller: " << this->Controller << endl;
      }
    else
      {
      os << indent << "Controller: (null)" << endl;
      }
  }
  
  //------------------------------------------------------------------------------
  int vtkPGenericIOMultiBlockReader::GetNumberOfPointArrays()
  {
    return this->PointDataArraySelection->GetNumberOfArrays();
  }
  
  //------------------------------------------------------------------------------
  const char* vtkPGenericIOMultiBlockReader::GetPointArrayName(int i)
  {
    assert("pre: array index is out-of-bounds!" &&
      (i >= 0) && (i < this->GetNumberOfPointArrays()));
    return this->PointDataArraySelection->GetArrayName(i);
  }
  
  //------------------------------------------------------------------------------
  int vtkPGenericIOMultiBlockReader::GetPointArrayStatus(const char* name)
  {
    assert(this->PointDataArraySelection->ArrayExists(name));
    return this->PointDataArraySelection->ArrayIsEnabled(name);
  }
  
  //------------------------------------------------------------------------------
  void vtkPGenericIOMultiBlockReader::SetPointArrayStatus(
    const char* name, int status)
  {
    assert(this->PointDataArraySelection->ArrayExists(name));
    if (status)
      {
      this->PointDataArraySelection->EnableArray(name);
      }
    else
      {
      this->PointDataArraySelection->DisableArray(name);
      }
  }
  
  //------------------------------------------------------------------------------
  vtkIdType vtkPGenericIOMultiBlockReader::GetRequestedHaloId(vtkIdType i)
  {
    assert("pre: array index out of bounds" &&
           (i >= 0 && this->HaloList->GetNumberOfIds() > i));
    return this->HaloList->GetId(i);
  }
  
  //------------------------------------------------------------------------------
  vtkIdType vtkPGenericIOMultiBlockReader::GetNumberOfRequestedHaloIds()
  {
    return this->HaloList->GetNumberOfIds();
  }
  
  //------------------------------------------------------------------------------
  void vtkPGenericIOMultiBlockReader::SetNumberOfRequestedHaloIds(vtkIdType numIds)
  {
    this->HaloList->SetNumberOfIds(numIds);
    this->Modified();
  }
  
  //------------------------------------------------------------------------------
  void vtkPGenericIOMultiBlockReader::AddRequestedHaloId(vtkIdType haloId)
  {
    this->SetRequestedHaloId(this->GetNumberOfRequestedHaloIds(),haloId);
  }
  
  //------------------------------------------------------------------------------
  void vtkPGenericIOMultiBlockReader::ClearRequestedHaloIds()
  {
    this->HaloList->Reset();
    this->Modified();
  }
  
  //------------------------------------------------------------------------------
  void vtkPGenericIOMultiBlockReader::SetRequestedHaloId(vtkIdType i, vtkIdType haloId)
  {
    *this->HaloList->WritePointer(i,1) = haloId;
    this->Modified();
  }
  
  //------------------------------------------------------------------------------
  bool vtkPGenericIOMultiBlockReader::ReaderParametersChanged()
  {
    assert("pre: internal reader is NULL!" && (this->Reader != NULL) );
  
    if(this->Reader->GetFileName() != std::string(this->FileName) )
      {
  #ifdef DEBUG
      std::cout << "\t[INFO]: File name has changed!\n";
      std::cout.flush();
  #endif
      return true;
      }
  
    bool status = false;
    switch(this->Reader->GetIOStrategy())
      {
      case gio::GenericIOBase::FileIOMPI:
        status = (this->GenericIOType!=IOTYPEMPI)? true : false;
  #ifdef DEBUG
        if(status == true)
          {
          std::cout << "\t[INFO]: I/O strategy changed from MPI\n";
          std::cout.flush();
          }
  #endif
        break;
      case gio::GenericIOBase::FileIOPOSIX:
        status = (this->GenericIOType!=IOTYPEPOSIX)? true : false;
  #ifdef DEBUG
        if(status == true)
          {
          std::cout << "\t[INFO]: I/O strategy changed from POSIX\n";
          std::cout.flush();
          }
  #endif
        break;
      default:
        assert("pre: Code should not reach here!" && false);
        throw std::runtime_error("Invalid I/O strategy!\n");
      } // END switch on I/O strategy
  
    if( status == true )
      {
      /* short-circuit here */
      return( status );
      }
  
    switch(this->Reader->GetBlockAssignmentStrategy())
      {
      case gio::RR_BLOCK_ASSIGNMENT:
        status = (this->BlockAssignment!=ROUND_ROBIN)? true : false;
  #ifdef DEBUG
        if(status == true)
          {
          std::cout << "\t[INFO]: I/O block assignment changed to Round-Robin\n";
          std::cout.flush();
          }
  #endif
        break;
      case gio::RCB_BLOCK_ASSIGNMENT:
        status = (this->BlockAssignment!=RCB)? true : false;
  #ifdef DEBUG
        if(status == true)
          {
          std::cout << "\t[INFO]: I/O block assignment changed to RCB\n";
          std::cout.flush();
          }
  #endif
        break;
      default:
        assert("pre: Code should not reach here!" && false);
        throw std::runtime_error("Invalid BlockAssignment strategy!\n");
      } // END switch on BlockAssignment
  
    return( status );
  }
  
  //------------------------------------------------------------------------------
  gio::GenericIOReader* vtkPGenericIOMultiBlockReader::GetInternalReader()
  {
    if( this->Reader != NULL )
      {
      if( this->ReaderParametersChanged() )
        {
  #ifdef DEBUG
        std::cout << "\t[INFO]: Deleting Reader instance...\n";
        std::cout.flush();
  #endif
        this->Reader->Close();
        delete this->Reader;
        this->Reader=NULL;
        } // END if the reader parameters
      else
        {
        return( this->Reader );
        }
      } // END if the reader is not NULL
  
    this->BuildMetaData = true; // signal to re-build metadata
  
    assert("pre: Reader should be NULL" && (this->Reader==NULL));
    gio::GenericIOReader *r = NULL;
    bool posix              = (this->GenericIOType==IOTYPEMPI)? false : true;
    int distribution        = (this->BlockAssignment==RCB)?
                                  gio::RCB_BLOCK_ASSIGNMENT :
                                  gio::RR_BLOCK_ASSIGNMENT;
  
    r = vtkGenericIOUtilities::GetReader(
          vtkGenericIOUtilities::GetMPICommunicator(this->Controller),
          posix,distribution,std::string(this->FileName));
    assert("post: Internal GenericIO reader should not be NULL!" && (r!=NULL) );
  
    return( r );
  }
  
  //------------------------------------------------------------------------------
  void vtkPGenericIOMultiBlockReader::LoadMetaData()
  {
    if( !this->BuildMetaData )
      {
  #ifdef DEBUG
      std::cout << "\t[INFO]: No need to update metadata!\n";
      std::cout.flush();
  #endif
      return;
      }
  
    this->MetaData->Clear();
    this->PointDataArraySelection->RemoveAllArrays();
    this->ArrayList->SetNumberOfValues(0);
  
  #ifdef DEBUG
      std::cout << "\t[INFO]: Reading header to build metadata!\n";
      std::cout << "\t[INFO]: Filename: " << this->FileName << std::endl;
      std::cout.flush();
  #endif
    this->Reader->OpenAndReadHeader();
    this->MetaData->TotalNumberOfElements = this->Reader->GetNumberOfElements();
  
    // load variable information
    for (int i = 0; i < this->Reader->GetNumberOfVariablesInFile(); ++i)
      {
      std::string vname = this->Reader->GetVariableName(i);
      this->ArrayList->InsertNextValue( vname.c_str() );
      this->PointDataArraySelection->AddArray( vname.c_str() );
      this->PointDataArraySelection->DisableArray( vname.c_str() );
  
      this->MetaData->VariableInformation[vname]=
        this->Reader->GetFileVariableInfo(i);
      this->MetaData->VariableGenericIOType[vname]=
        gio::GenericIOUtilities::DetectVariablePrimitiveType(
          this->MetaData->VariableInformation[vname]);
      } // end for all variables in file
  
    // load block information
    this->MetaData->NumberOfBlocks = this->Reader->GetTotalNumberOfBlocks();
  
    std::vector< most_of_block_t > localBlocks;
    localBlocks.resize(this->Reader->GetNumberOfAssignedBlocks());
    double min[3];
    double max[3];
    int myProcessId = this->Controller->GetLocalProcessId();
    for (int i = 0; i < this->Reader->GetNumberOfBlockHeaders(); ++i)
      {
      most_of_block_t& myBlock = localBlocks[i];
      gio::RankHeader block = this->Reader->GetBlockHeader(i);
      assert("pre: loading duplicate block in metadata!" &&
        !this->MetaData->HasBlock(block.GlobalRank));
  
      myBlock.GlobalId = block.GlobalRank;
      myBlock.NumberOfElements = block.NElems;
  
  
      if (this->Reader->IsSpatiallyDecomposed())
        {
        this->Reader->GetBlockBounds(i,min,max);
        this->Reader->GetBlockCoords(i,myBlock.coords);
        myBlock.bounds[0] = min[0];
        myBlock.bounds[1] = max[0];
        myBlock.bounds[2] = min[1];
        myBlock.bounds[3] = max[1];
        myBlock.bounds[4] = min[2];
        myBlock.bounds[5] = max[2];
        }
      else
        {
        for (int j = 0; j < 6; ++j)
          {
          myBlock.bounds[j] = 0;
          }
        myBlock.coords[0] = myBlock.coords[1] = myBlock.coords[2] = 0;
        }
      myBlock.ProcessId = myProcessId;
      }
    std::vector< most_of_block_t > allBlocks;
    allBlocks.resize(this->Reader->GetTotalNumberOfBlocks());
    assert(allBlocks.size() == localBlocks.size() * this->Controller->GetNumberOfProcesses());
    this->Controller->AllGather((const char*)&localBlocks[0],(char*)&allBlocks[0],
        localBlocks.size() * sizeof(most_of_block_t));
  
    for (unsigned i = 0; i < allBlocks.size(); ++i)
      {
      most_of_block_t& blockData = allBlocks[i];
      block_t block;
      block.GlobalId = blockData.GlobalId;
      block.NumberOfElements = blockData.NumberOfElements;
      block.ProcessId = blockData.ProcessId;
      block.coords[0] = blockData.coords[0];
      block.coords[1] = blockData.coords[1];
      block.coords[2] = blockData.coords[2];
      block.bounds[0] = blockData.bounds[0];
      block.bounds[1] = blockData.bounds[1];
      block.bounds[2] = blockData.bounds[2];
      block.bounds[3] = blockData.bounds[3];
      block.bounds[4] = blockData.bounds[4];
      block.bounds[5] = blockData.bounds[5];
      if (block.ProcessId == myProcessId)
        {
        for (int var = 0; var < this->Reader->GetNumberOfVariablesInFile(); ++var)
          {
          std::string vname = this->Reader->GetVariableName(var);
          block.RawCache[vname] = NULL;
          block.VariableStatus[vname] = false;
          }
        }
      this->MetaData->Blocks[ block.GlobalId ] = block;
      }
  
    this->BuildMetaData = false;
  }
  
  //------------------------------------------------------------------------------
  void vtkPGenericIOMultiBlockReader::LoadRawVariableDataForBlock(
    const std::string& varName, int blockId)
  {
  #ifdef DEBUG
    std::cout << "[INFO]: Loading variable: " << varName << std::endl;
    std::cout.flush();
  #endif
    assert("pre: Reader is NULL" && (this->Reader != NULL));
    assert("pre: metadata is NULL" && (this->MetaData != NULL));
    assert("pre: metadata is corrupt" && (this->MetaData->SanityCheck()));
    assert("pre: block is owned by another process" &&
      this->MetaData->HasBlock(blockId));
    assert("pre: metadata does not have requested variable" &&
      this->MetaData->HasVariable(varName));
  
    block_t& dataBlock = this->MetaData->Blocks[blockId];
  
    assert("pre: no variable by this name!" &&
      (dataBlock.RawCache.find(varName) != dataBlock.RawCache.end()) &&
      (dataBlock.VariableStatus.find(varName) != dataBlock.VariableStatus.end()));
  
    if (dataBlock.VariableStatus[varName])
      {
  #ifdef DEBUG
      std::cout << "\t[INFO]: Variable appears to be already loaded "
                << "for BLOCK=" << blockId << std::endl;
      std::cout.flush();
  #endif
      return;
      }
  
    dataBlock.RawCache[varName] =
      gio::GenericIOUtilities::AllocateVariableArray(
        this->MetaData->VariableInformation[varName],dataBlock.NumberOfElements);
  
    this->Reader->AddVariable(
      this->MetaData->VariableInformation[varName],dataBlock.RawCache[varName]);
  
    dataBlock.VariableStatus[varName] = true;
  
  #ifdef DEBUG
    std::cout << "\t[INFO]: Variable [" << varName << "] is now loaded "
              << "for BLOCK=" << blockId << std::endl;
    std::cout.flush();
  #endif
  
  }
  
  //------------------------------------------------------------------------------
  void vtkPGenericIOMultiBlockReader::LoadRawDataForBlock(int blockId)
  {
    assert("pre: Reader is NULL" && (this->Reader != NULL));
    assert("pre: metadata is NULL" && (this->MetaData != NULL));
    assert("pre: metadata is corrupt!" && (this->MetaData->SanityCheck()));
    assert("pre: block is not owned by this process!" &&
      this->MetaData->HasBlock(blockId));
  
    // This method is called for every block, so we must clear any previously
    // registered variables
    this->Reader->ClearVariables();
  
    std::string xaxis = std::string(this->XAxisVariableName);
    xaxis = vtkGenericIOUtilities::trim(xaxis);
  
    std::string yaxis = std::string(this->YAxisVariableName);
    yaxis = vtkGenericIOUtilities::trim(yaxis);
  
    std::string zaxis = std::string(this->ZAxisVariableName);
    zaxis = vtkGenericIOUtilities::trim(zaxis);
  
    this->LoadRawVariableDataForBlock( xaxis, blockId );
    this->LoadRawVariableDataForBlock( yaxis, blockId );
    this->LoadRawVariableDataForBlock( zaxis, blockId );
  
    if (this->HaloList->GetNumberOfIds() > 0)
      {
      std::string haloIds = std::string(this->HaloIdVariableName);
      haloIds = vtkGenericIOUtilities::trim(haloIds);
      this->LoadRawVariableDataForBlock(haloIds,blockId);
      }
  
  #ifdef DEBUG
    std::cout << "\t==========\n";
    std::cout << "\tNUMBER OF ARRAYS: "
        << this->PointDataArraySelection->GetNumberOfArrays() << std::endl;
    std::cout.flush();
  #endif
  
    int arrayIdx = 0;
    for(;arrayIdx < this->PointDataArraySelection->GetNumberOfArrays(); ++arrayIdx)
      {
      const char *name = this->PointDataArraySelection->GetArrayName(arrayIdx);
  #ifdef DEBUG
      std::cout << "\tARRAY " << name << " is ";
  #endif
      if( this->PointDataArraySelection->ArrayIsEnabled(name) )
        {
  #ifdef DEBUG
        std::cout << "ENABLED\n";
        std::cout.flush();
  #endif
        std::string varName = std::string( name );
        this->LoadRawVariableDataForBlock( varName, blockId );
        } // END if the array is enabled
      else
        {
  #ifdef DEBUG
        std::cout << "DISBLED\n";
        std::cout.flush();
  #endif
        }
      } // END for all arrays
  
  #ifdef DEBUG
    std::cout << "\t[INFO]: Reading data...";
  #endif
  
    this->Reader->ReadBlock(blockId);
  
  #ifdef DEBUG
    std::cout << "[DONE]\n";
    std::cout.flush();
  #endif
  }
  
  //------------------------------------------------------------------------------
  void vtkPGenericIOMultiBlockReader::GetPointFromRawData(
          int xType, void* xBuffer, int yType, void* yBuffer, int zType, void* zBuffer,
          vtkIdType idx, double pnt[3])
  {
    void* buffer[3] = {xBuffer, yBuffer, zBuffer};
  
    int type[3] = {xType, yType, zType};
  
    for( int i=0; i < 3; ++i )
      {
  //    type = this->MetaData->VariableGenericIOType[ name[i] ];
  //    void *rawBuffer = this->MetaData->Blocks[ blockId ].RawCache[ name[i] ];
      assert("pre: raw buffer is NULL!" && (buffer[i] != NULL) );
  
      pnt[i] = vtkGenericIOUtilities::GetDoubleFromRawBuffer(type[i],buffer[i],idx);
      } // END for all dimensions
  }
  
  //------------------------------------------------------------------------------
  void vtkPGenericIOMultiBlockReader::LoadCoordinatesForBlock(
      vtkUnstructuredGrid *grid, std::set< vtkIdType >& pointsInSelectedHalos,
      int blockId)
  {
    assert("pre: metadata is NULL!" && (this->MetaData != NULL));
    assert("pre: grid is NULL!" && (grid != NULL) );
    assert("pre: block is not owned by this process!" &&
      this->MetaData->HasBlock(blockId));
  
    std::string xaxis = std::string(this->XAxisVariableName);
    xaxis = vtkGenericIOUtilities::trim(xaxis);
  
    std::string yaxis = std::string(this->YAxisVariableName);
    yaxis = vtkGenericIOUtilities::trim(yaxis);
  
    std::string zaxis = std::string(this->ZAxisVariableName);
    zaxis = vtkGenericIOUtilities::trim(zaxis);
  
    if( !this->MetaData->HasVariable(xaxis) ||
         !this->MetaData->HasVariable(yaxis) ||
         !this->MetaData->HasVariable(zaxis))
      {
      vtkErrorMacro(<< "Don't have one or more coordinate arrays!\n");
      return;
      }
    block_t& dataBlock = this->MetaData->Blocks[blockId];
  
    int xType = this->MetaData->VariableGenericIOType[xaxis];
    void* xBuffer = dataBlock.RawCache[xaxis];
    int yType = this->MetaData->VariableGenericIOType[yaxis];
    void* yBuffer = dataBlock.RawCache[yaxis];
    int zType = this->MetaData->VariableGenericIOType[zaxis];
    void* zBuffer = dataBlock.RawCache[zaxis];
  
  
    int nparticles = dataBlock.NumberOfElements;
  
    vtkSmartPointer< vtkCellArray > cells =
        vtkSmartPointer< vtkCellArray >::New();
    cells->Allocate(cells->EstimateSize(nparticles,1));
  
    vtkSmartPointer< vtkPoints > pnts =
        vtkSmartPointer< vtkPoints >::New();
    pnts->SetDataTypeToDouble();
    pnts->SetNumberOfPoints( nparticles );
  
  
    double pnt[3];
    vtkIdType idx = 0;
    if (this->HaloList->GetNumberOfIds() == 0)
      {
      for( ;idx < nparticles; ++idx)
        {
        this->GetPointFromRawData(xType, xBuffer, yType, yBuffer, zType, zBuffer,
                                  idx, pnt);
        pnts->SetPoint(idx,pnt);
        cells->InsertNextCell(1,&idx);
        } // END for all points
      }
    else
      {
      std::string haloVarName = std::string(this->HaloIdVariableName);
      haloVarName = vtkGenericIOUtilities::trim(haloVarName);
      int haloType = this->MetaData->VariableGenericIOType[haloVarName];
      void* haloBuffer = dataBlock.RawCache[haloVarName];
      for (vtkIdType i = 0; idx < nparticles; ++idx)
        {
        vtkIdType haloId = vtkGenericIOUtilities::GetIdFromRawBuffer(haloType,haloBuffer,idx);
        bool isInRequestedHalo = false;
        for (vtkIdType j = 0; j < this->GetNumberOfRequestedHaloIds(); ++j)
          {
          if (haloId == this->HaloList->GetId(j))
            {
            isInRequestedHalo = true;
            pointsInSelectedHalos.insert(idx);
            break;
            }
          }
        if (isInRequestedHalo)
          {
          this->GetPointFromRawData(xType,xBuffer,yType,yBuffer,zType,zBuffer,idx,pnt);
          pnts->SetPoint(i,pnt);
          cells->InsertNextCell(1,&i);
          ++i;
          }
        }
      pnts->SetNumberOfPoints(pointsInSelectedHalos.size());
      }
  
    grid->SetPoints(pnts);
  
    grid->SetCells(VTK_VERTEX,cells);
  
    grid->Squeeze();
  }
  
  namespace {
  template< typename T >
  void GetOnlyDataInHalo(vtkDataArray* allData, vtkDataArray* haloData, std::set< vtkIdType > pointsInHalo)
  {
    T* data = (T*) allData->GetVoidPointer(0);
    T* filteredData = (T*) haloData->GetVoidPointer(0);
    vtkIdType i = 0;
    for (std::set< vtkIdType >::iterator itr = pointsInHalo.begin();
         itr != pointsInHalo.end(); ++itr)
      {
      filteredData[i++] = data[*itr];
      }
  }
  }
  
  //------------------------------------------------------------------------------
  void vtkPGenericIOMultiBlockReader::LoadDataArraysForBlock(
      vtkUnstructuredGrid *grid, const std::set< vtkIdType >& pointsInSelectedHalos,
      int blockId)
  {
    assert("pre: metadata is NULL!" && (this->MetaData != NULL));
    assert("pre: grid is NULL!" && (grid != NULL) );
    assert("pre: block is not owned by this process!" &&
      this->MetaData->HasBlock(blockId));
  
    block_t& dataBlock = this->MetaData->Blocks[blockId];
  
  //  assert("pre: # points in dataset different from points in block" &&
  //    (static_cast<uint64_t>(grid->GetNumberOfPoints()) == dataBlock.NumberOfElements));
    vtkPointData* PD = grid->GetPointData();
  
    int arrayIdx = 0;
    for(;arrayIdx < this->PointDataArraySelection->GetNumberOfArrays(); ++arrayIdx)
      {
      const char *name = this->PointDataArraySelection->GetArrayName( arrayIdx );
      if( this->PointDataArraySelection->ArrayIsEnabled(name) )
        {
        std::string varName( name );
        vtkSmartPointer< vtkDataArray > dataArray;
        dataArray.TakeReference(
            vtkGenericIOUtilities::GetVtkDataArray(
                varName,
                this->MetaData->VariableGenericIOType[ varName ],
                dataBlock.RawCache[ varName ],
                dataBlock.NumberOfElements
                ));
        if (this->HaloList->GetNumberOfIds() != 0)
          {
          vtkSmartPointer< vtkDataArray > onlyDataInHalo;
          onlyDataInHalo.TakeReference(dataArray->NewInstance());
          onlyDataInHalo->SetNumberOfTuples(grid->GetNumberOfPoints());
          onlyDataInHalo->SetName(dataArray->GetName());
          switch (dataArray->GetDataType()) {
          vtkTemplateMacro(GetOnlyDataInHalo<VTK_TT>(
                             dataArray,onlyDataInHalo,pointsInSelectedHalos));
            }
          dataArray = onlyDataInHalo;
          }
  
  
        PD->AddArray( dataArray );
        } // END if the array is enabled
      } // END for all arrays
  
  }
  
  //------------------------------------------------------------------------------
  vtkUnstructuredGrid* vtkPGenericIOMultiBlockReader::LoadBlock(int blockId)
  {
    // Sanity Check
    assert("pre: metadata is null" && (this->MetaData != NULL));
    assert("pre: block is not owned by this process!" &&
      this->MetaData->HasBlock(blockId));
    // STEP 1: Load raw data
    this->LoadRawDataForBlock(blockId);
  
    vtkUnstructuredGrid* grid = vtkUnstructuredGrid::New();
    std::set< vtkIdType > pointsInSelectedHalos;
  
    // STEP 2: Load coordinates
    this->LoadCoordinatesForBlock(grid,pointsInSelectedHalos,blockId);
  
    // STEP 3: Load data
    this->LoadDataArraysForBlock(grid,pointsInSelectedHalos,blockId);
  
    if (this->Reader->IsSpatiallyDecomposed())
      {
      vtkSmartPointer< vtkTypeUInt64Array > coords =
          vtkSmartPointer< vtkTypeUInt64Array >::New();
      coords->SetNumberOfComponents(3);
      coords->SetNumberOfTuples(1);
      coords->SetName("genericio_block_coords");
      coords->SetTupleValue(0,(vtkTypeUInt64*)this->MetaData->Blocks[blockId].coords);
      grid->GetFieldData()->AddArray(coords);
      }
  
    return grid;
  }
  
  //------------------------------------------------------------------------------
  void vtkPGenericIOMultiBlockReader::SelectionModifiedCallback(
      vtkObject* vtkNotUsed(caller), unsigned long vtkNotUsed(eid),
      void *clientdata, void* vtkNotUsed(calldata))
  {
    assert("pre: clientdata is NULL!" && (clientdata != NULL) );
    static_cast<vtkPGenericIOMultiBlockReader*>(clientdata)->Modified();
  }
  
  //------------------------------------------------------------------------------
  int vtkPGenericIOMultiBlockReader::RequestInformation(
        vtkInformation *vtkNotUsed(request),
        vtkInformationVector **vtkNotUsed(inputVector),
        vtkInformationVector *outputVector)
  {
  #ifdef DEBUG
    std::cout << "[INFO]: Call RequestInformation()\n";
    std::cout.flush();
  #endif
  
  //  ++this->RequestInfoCounter;
  
    // tell the pipeline that this dataset is distributed
    outputVector->GetInformationObject(0)->Set(
        CAN_HANDLE_PIECE_REQUEST(),
        1
        );
    outputVector->GetInformationObject(0)->Set(
        vtkDataObject::DATA_NUMBER_OF_PIECES(),
        this->Controller->GetNumberOfProcesses()
        );
  
    this->Reader = this->GetInternalReader();
    assert("pre: internal reader is NULL!" && (this->Reader != NULL) );
  
    this->LoadMetaData();
  
    vtkSmartPointer< vtkMultiBlockDataSet > outline =
      vtkSmartPointer< vtkMultiBlockDataSet >::New();
    outline->SetNumberOfBlocks(this->MetaData->NumberOfBlocks);
    int myProcessId = this->Controller->GetLocalProcessId();
    for (int i = 0; i < this->MetaData->NumberOfBlocks; ++i)
      {
      outline->SetBlock(i,NULL);
      if (this->MetaData->HasBlock(i))
        {
        vtkInformation* blockInfo = outline->GetMetaData(i);
        assert("pre: block info is NULL!" && (blockInfo != NULL));
        // bounds are meaningless if this is false
        if (this->Reader->IsSpatiallyDecomposed())
          {
          blockInfo->Set(
                vtkStreamingDemandDrivenPipeline::BOUNDS(),
                this->MetaData->Blocks[i].bounds,6
                );
          }
        blockInfo->Set(
          vtkCompositeDataPipeline::BLOCK_AMOUNT_OF_DETAIL(),
              this->MetaData->Blocks[i].NumberOfElements);
        blockInfo->Set(
              vtkCompositeDataSet::CURRENT_PROCESS_CAN_LOAD_BLOCK(),
              this->MetaData->Blocks[i].ProcessId == myProcessId ? 1 : 0);
        }
      }
    outline->GetInformation()->Set(vtkCompositeDataPipeline::BLOCK_AMOUNT_OF_DETAIL(),
                                   this->MetaData->TotalNumberOfElements);
    outputVector->GetInformationObject(0)->Set(
      vtkCompositeDataPipeline::COMPOSITE_DATA_META_DATA(), outline);
    return 1;
  }
  
  //------------------------------------------------------------------------------
  int vtkPGenericIOMultiBlockReader::RequestData(
      vtkInformation *vtkNotUsed(request),
      vtkInformationVector **vtkNotUsed(inputVector),
      vtkInformationVector *outputVector)
  {
  #ifdef DEBUG
    std::cout << "[INFO]: Call RequestData()\n";
    std::cout.flush();
  #endif
  
    // STEP 0: get the output grid
    vtkMultiBlockDataSet* output = vtkMultiBlockDataSet::GetData(outputVector,0);
    assert("pre: output dataset is NULL!" && (output != NULL) );
  
    vtkInformation* outInfo = outputVector->GetInformationObject(0);
  
    output->SetNumberOfBlocks(this->MetaData->NumberOfBlocks);
  
    // Get the global dimensions and physical orign & scale from
    // the genericio file and add them to the dataset
    uint64_t tmpDims[3];
    double tmpDouble[3];
    vtkSmartPointer< vtkTypeUInt64Array > dims =
        vtkSmartPointer< vtkTypeUInt64Array >::New();
    vtkSmartPointer< vtkDoubleArray > origin =
        vtkSmartPointer< vtkDoubleArray >::New();
    vtkSmartPointer< vtkDoubleArray > scale =
        vtkSmartPointer< vtkDoubleArray >::New();
    dims->SetNumberOfComponents(3);
    dims->SetName("genericio_global_dimensions");
    origin->SetNumberOfComponents(3);
    origin->SetName("genericio_phys_origin");
    scale->SetNumberOfComponents(3);
    scale->SetName("genericio_phys_scale");
  
    this->Reader->GetGlobalDimensions(tmpDims);
    dims->InsertNextTupleValue((vtkTypeUInt64*)tmpDims);
    this->Reader->GetPhysOrigin(tmpDouble);
    origin->InsertNextTupleValue(tmpDouble);
    this->Reader->GetPhysScale(tmpDouble);
    scale->InsertNextTupleValue(tmpDouble);
  
    output->GetFieldData()->AddArray(origin);
    output->GetFieldData()->AddArray(scale);
    output->GetFieldData()->AddArray(dims);
  
    int myProcessId = this->Controller->GetLocalProcessId();
  
    if (outInfo->Has(vtkCompositeDataPipeline::LOAD_REQUESTED_BLOCKS()))
      {
      int size = outInfo->Length(vtkCompositeDataPipeline::UPDATE_COMPOSITE_INDICES());
      int* ids = outInfo->Get(vtkCompositeDataPipeline::UPDATE_COMPOSITE_INDICES());
      for (int i = 0; i < size; ++i)
        {
        int blockId = ids[i];
        if (this->MetaData->BlockIsOnMyProcess(blockId,myProcessId))
          {
          vtkSmartPointer< vtkUnstructuredGrid > grid =
            vtkSmartPointer< vtkUnstructuredGrid >::Take(
              this->LoadBlock(blockId));
          output->SetBlock(blockId,grid);
          }
        }
      }
    else
      {
      std::map< int, block_t >::iterator blockItr = this->MetaData->Blocks.begin();
      for ( ; blockItr != this->MetaData->Blocks.end(); ++blockItr)
        {
        if (blockItr->second.ProcessId == myProcessId)
          {
          vtkSmartPointer< vtkUnstructuredGrid > grid =
              vtkSmartPointer< vtkUnstructuredGrid >::Take(
                this->LoadBlock(blockItr->first));
          output->SetBlock(blockItr->first,grid);
          }
        }
      }
  
    return 1;
  }
  