vtkXdmfReader.cxx

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

  Program:   Visualization Toolkit
  Module:    vtkXdmfReader.cxx
  Language:  C++
  Date:      $Date$
  Version:   $Revision$


Copyright (c) 1993-2001 Ken Martin, Will Schroeder, Bill Lorensen  
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

 * Redistributions of source code must retain the above copyright notice,
   this list of conditions and the following disclaimer.

 * Redistributions in binary form must reproduce the above copyright notice,
   this list of conditions and the following disclaimer in the documentation
   and/or other materials provided with the distribution.

 * Neither name of Ken Martin, Will Schroeder, or Bill Lorensen nor the names
   of any contributors may be used to endorse or promote products derived
   from this software without specific prior written permission.

 * Modified source versions must be plainly marked as such, and must not be
   misrepresented as being the original software.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ``AS IS''
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

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

#include "vtkCallbackCommand.h"
#include "vtkDataArraySelection.h"
#include "vtkDataSet.h"
#include "vtkDoubleArray.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkObjectFactory.h"
#include "vtkRectilinearGrid.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkStructuredGrid.h"
#include "vtkUnstructuredGrid.h"
#include "vtkXdmfDataArray.h"
#include "vtkPointData.h"
#include "vtkCellData.h"
#include "vtkCell.h"
#include "vtkCellArray.h"
#include "vtkCharArray.h"
#include "vtkXMLParser.h"
#include "vtkImageData.h"
#include "vtkUniformGrid.h"
#include "vtkMultiGroupDataInformation.h"
#include "vtkMultiGroupDataSet.h"
// #include "vtkHierarchicalDataSet.h"
#include "vtkCompositeDataPipeline.h"
#include "vtkMultiProcessController.h"

#include "XdmfArray.h"
#include "XdmfAttribute.h"
#include "XdmfDOM.h"
#include "XdmfDataDesc.h"
#include "XdmfDataItem.h"
#include "XdmfGrid.h"
#include "XdmfTopology.h"
#include "XdmfGeometry.h"

#include <sys/stat.h>
#include <vtkstd/set>
#include <vtkstd/map>
#include <vtkstd/string>
#include <vtkstd/vector>
#include <assert.h>

#define USE_IMAGE_DATA // otherwise uniformgrid

vtkStandardNewMacro(vtkXdmfReader);
vtkCxxRevisionMacro(vtkXdmfReader, "1.16");

vtkCxxSetObjectMacro(vtkXdmfReader,Controller,vtkMultiProcessController);

#if defined(_WIN32) && (defined(_MSC_VER) || defined(__BORLANDC__))
#  include <direct.h>
#  define GETCWD _getcwd
#else
#include <unistd.h>
#  define GETCWD getcwd
#endif


#define vtkMAX(x, y) (((x)>(y))?(x):(y))
#define vtkMIN(x, y) (((x)<(y))?(x):(y))

#define PRINT_EXTENT(x) "[" << (x)[0] << " " << (x)[1] << " " << (x)[2] << " " << (x)[3] << " " << (x)[4] << " " << (x)[5] << "]" 

//----------------------------------------------------------------------------
class vtkXdmfReaderGrid
{
public:
  vtkXdmfReaderGrid() : XMGrid(0), DataDescription(0) {}
  ~vtkXdmfReaderGrid()
    {
    delete this->XMGrid;
    }

  XdmfGrid       *XMGrid;
  XdmfDataDesc   *DataDescription;
  vtkstd::string Name;
  int Level;
};

//----------------------------------------------------------------------------
class vtkXdmfReaderGridCollection
{
public:
  vtkXdmfReaderGrid* GetXdmfGrid(const char* gridName,
                                 int level);

  typedef vtkstd::map<vtkstd::string, vtkXdmfReaderGrid*> SetOfGrids;
  SetOfGrids Grids;
  // Update number of levels and number of datasets.
  void UpdateCounts();
  
  // Return the last number of levels computed by UpdateCounts.
  int GetNumberOfLevels()
    {
      return this->NumberOfLevels;
    }
  // Return the last number of dataset for a given level
  // computed by UpdateCounts.
  int GetNumberOfDataSets(int level)
    {
      assert("pre: valid_level" && level>=0 && level<GetNumberOfLevels());
      return this->NumberOfDataSets[level];
    }
  
protected:
  int NumberOfLevels;
  vtkstd::vector<int> NumberOfDataSets;
};

//----------------------------------------------------------------------------
vtkXdmfReaderGrid* vtkXdmfReaderGridCollection::GetXdmfGrid(
  const char *gridName,
  int level)
{
  vtkXdmfReaderGridCollection::SetOfGrids::iterator it = this->Grids.find(gridName);
  if ( it == this->Grids.end() || it->second == 0 )
    {
    this->Grids[gridName] = new vtkXdmfReaderGrid;
    }
  this->Grids[gridName]->Level=level;
  return this->Grids[gridName];
}

//----------------------------------------------------------------------------
void vtkXdmfReaderGridCollection::UpdateCounts()
{
  // Update the number of levels.
  vtkXdmfReaderGridCollection::SetOfGrids::iterator it;
  vtkXdmfReaderGridCollection::SetOfGrids::iterator itEnd;
  
  int maxLevel=0;
  it=this->Grids.begin();
  itEnd=this->Grids.end();
  while(it!=itEnd)
    {
    int level=it->second->Level;
    if(level>maxLevel)
      {
      maxLevel=level;
      }
    ++it;
    }
  this->NumberOfLevels=maxLevel+1;
  this->NumberOfDataSets.resize(this->NumberOfLevels);
 
  // Update the number of datasets
  int l=this->NumberOfLevels;
  int i=0;
  while(i<l)
    {
    this->NumberOfDataSets[i]=0;
    ++i;
    }
  
  it=this->Grids.begin();
  while(it!=itEnd)
    {
    int level=it->second->Level;
    ++this->NumberOfDataSets[level];
    ++it;
    }
}

//----------------------------------------------------------------------------
class vtkXdmfReaderActualGrid
{
public:
  vtkXdmfReaderActualGrid() : Enabled(0), Grid(0), Collection(0) {}

  int Enabled;
  vtkXdmfReaderGrid* Grid;
  vtkXdmfReaderGridCollection* Collection;
};

//----------------------------------------------------------------------------
class vtkXdmfReaderInternal
{
public:
  vtkXdmfReaderInternal()
    {
    this->DataItem = 0;
    this->ArrayConverter = vtkXdmfDataArray::New();
    this->NumberOfEnabledActualGrids=0;
    }
  ~vtkXdmfReaderInternal()
    {
    if ( this->DataItem )
      {
      // cout << "....Deleting DataItem " << this->DataItem << endl;
      delete this->DataItem;
      // cout << "....Back from Deleting DataItem " << this->DataItem << endl;
      this->DataItem = 0;
      }
    this->ArrayConverter->Delete();
    this->ArrayConverter = 0;
    }

  typedef vtkstd::vector<vtkstd::string> StringListType;

  int RequestSingleGridData(const char* currentGridName,
                            vtkXdmfReaderGrid *grid,
                            vtkInformation *outInfo,
                            vtkDataObject *output,
                            int isSubBlock);
  
  int RequestActualGridData(const char* currentGridName,
                            vtkXdmfReaderActualGrid* currentActualGrid,
                            int outputGrid,
                            int numberOfGrids,
                            vtkInformationVector *outputVector);
  
  // outInfo is null in the  multi-block case.
  int RequestSingleGridInformation(vtkXdmfReaderGrid *grid,
                                   vtkInformation *outInfo);

  int RequestActualGridInformation(vtkXdmfReaderActualGrid* currentActualGrid,
                                   int outputGrid,
                                    int numberOfGrids,
                                   vtkInformationVector* outputVector);
  
  typedef vtkstd::map<vtkstd::string,vtkXdmfReaderActualGrid> MapOfActualGrids;
  MapOfActualGrids ActualGrids;
  int NumberOfEnabledActualGrids;
  
  vtkXdmfReaderGridCollection* GetCollection(const char* collectionName);
  vtkXdmfReaderActualGrid* GetGrid(const char* gridName);
  vtkXdmfReaderActualGrid* GetGrid(int idx);
  vtkXdmfReaderGrid *GetXdmfGrid(const char *gridName,
                                 const char *collectionName,
                                 const char *levelName);

  StringListType DomainList;
  XdmfXmlNode DomainPtr;
  vtkXdmfReader* Reader;
  XdmfDataItem *DataItem;

  // For converting arrays from XDMF to VTK format
  vtkXdmfDataArray *ArrayConverter;
};

//----------------------------------------------------------------------------
vtkXdmfReaderGridCollection* vtkXdmfReaderInternal::GetCollection(
  const char* collectionName)
{
  if ( !collectionName )
    {
    return 0;
    }
  vtkXdmfReaderActualGrid* actualGrid = &this->ActualGrids[collectionName];
  if ( !actualGrid->Collection )
    {
    if ( actualGrid->Grid )
      {
      cerr << "Trying to create collection with the same name as an existing grid" << endl;
      return 0;
      }
    actualGrid->Collection = new vtkXdmfReaderGridCollection;
    }
  return actualGrid->Collection;
}

//----------------------------------------------------------------------------
vtkXdmfReaderGrid* vtkXdmfReaderInternal::GetXdmfGrid(
  const char *gridName,
  const char *collectionName,
  const char *levelName)
{
  if ( !gridName )
    {
    return 0;
    }

  if ( collectionName )
    {
    vtkXdmfReaderGridCollection *collection=this->GetCollection(collectionName);
    if ( !collection )
      {
      return 0;
      }
    int level;
    if(levelName==0)
      {
      level=0;
      }
    else
      {
      char *tmp=new char[strlen(levelName)+1];
      strcpy(tmp,levelName);
      istrstream s(tmp,strlen(tmp));
      s>>level;
      if(level<0)
        {
        cerr << "Expect a positive Level value, level=" << level <<endl;
        delete[] tmp;
        return 0;
        }
      else
        {
        delete[] tmp;
        }
      }
    return collection->GetXdmfGrid(gridName,level);
    }

  vtkXdmfReaderActualGrid* grid = this->GetGrid(gridName);
  if ( grid->Collection )
    {
    cerr << "Trying to create a grid with the same name as an existing collection" << endl;
    return 0;
    }
  grid->Grid = new vtkXdmfReaderGrid;
  return grid->Grid;
}

//----------------------------------------------------------------------------
vtkXdmfReaderActualGrid* vtkXdmfReaderInternal::GetGrid(const char* gridName)
{
  return &this->ActualGrids[gridName];
}

//----------------------------------------------------------------------------
vtkXdmfReaderActualGrid* vtkXdmfReaderInternal::GetGrid(int idx)
{
  if ( idx < 0 )
    {
    return 0;
    }
  vtkXdmfReaderInternal::MapOfActualGrids::iterator it;
  int cnt = 0;
  for ( it = this->ActualGrids.begin();
    it != this->ActualGrids.end();
    ++ it )
    {
    if ( cnt == idx )
      {
      return &it->second;
      }
    cnt ++;
    }
  return 0;
}

//----------------------------------------------------------------------------
vtkXdmfReader::vtkXdmfReader()
{
//this->SetNumberOfInputPorts(0);
//this->SetNumberOfOutputPorts(0);
  
  this->Internals = new vtkXdmfReaderInternal;
  this->Internals->Reader = this;

  // I have this flag because I do not want to change the initialization
  // of the output to the generic output.  It might break something;
  this->OutputsInitialized = 0;

  this->DOM = 0;

  this->PointDataArraySelection = vtkDataArraySelection::New();
  this->CellDataArraySelection = vtkDataArraySelection::New();

  // Setup the selection callback to modify this object when an array
  // selection is changed.
  this->SelectionObserver = vtkCallbackCommand::New();
  this->SelectionObserver->SetCallback(
    &vtkXdmfReader::SelectionModifiedCallback);
  this->SelectionObserver->SetClientData(this);
  this->PointDataArraySelection->AddObserver(vtkCommand::ModifiedEvent,
                                             this->SelectionObserver);
  this->CellDataArraySelection->AddObserver(vtkCommand::ModifiedEvent,
                                            this->SelectionObserver);

  this->DomainName = 0;
  this->Internals->DomainPtr = 0;
  this->GridName = 0;
  
  for (int i = 0; i < 3; i ++ )
    {
    this->Stride[i] = 1;
    }

  this->GridsModified = 0;
  
  this->NumberOfEnabledActualGrids=0;
  
  this->Controller = 0;
  this->SetController(vtkMultiProcessController::GetGlobalController());
}

//----------------------------------------------------------------------------
vtkXdmfReader::~vtkXdmfReader()
{
  if ( this->DOM )
    {
    delete this->DOM;
    }

  this->CellDataArraySelection->RemoveObserver(this->SelectionObserver);
  this->PointDataArraySelection->RemoveObserver(this->SelectionObserver);
  this->SelectionObserver->Delete();
  this->CellDataArraySelection->Delete();
  this->PointDataArraySelection->Delete();

  this->SetDomainName(0);
  
  vtkXdmfReaderInternal::MapOfActualGrids::iterator actualGridIt;
  for ( actualGridIt = this->Internals->ActualGrids.begin();
    actualGridIt != this->Internals->ActualGrids.end();
    ++ actualGridIt )
    {
    vtkXdmfReaderActualGrid* grid = &actualGridIt->second;
    if ( grid->Grid )
      {
      delete grid->Grid;
      }
    if ( grid->Collection )
      {
      vtkXdmfReaderGridCollection::SetOfGrids::iterator gridIt;
      for ( gridIt = grid->Collection->Grids.begin();
        gridIt != grid->Collection->Grids.end();
        ++ gridIt )
        {
        delete gridIt->second;
        }
      delete grid->Collection;
      }
    }

  delete this->Internals;
  H5garbage_collect();
  
  this->SetController(0);
}

//----------------------------------------------------------------------------
int vtkXdmfReader::ProcessRequest(vtkInformation* request,
                                  vtkInformationVector** inputVector,
                                  vtkInformationVector* outputVector)
{
  // create the output
  if(request->Has(vtkDemandDrivenPipeline::REQUEST_DATA_OBJECT()))
    {
    return this->RequestDataObject(outputVector);
    }
  return this->Superclass::ProcessRequest(request, inputVector, outputVector);
}

//----------------------------------------------------------------------------
int vtkXdmfReader::RequestDataObject(vtkInformationVector *outputVector)
{
  // Set the number of outputs and the create the output data objects
  // needed because we may change the number of output ports
  // inside RequestDataObject. So we may change the output vector.
  vtkInformationVector *newOutputVector=outputVector;
  
  ////////////
  vtkIdType cc;
  XdmfConstString CurrentFileName;

  if ( !this->FileName )
    {
    vtkErrorMacro("File name not set");
    return 1;
    }
  // First make sure the file exists.  This prevents an empty file
  // from being created on older compilers.
  struct stat fs;
  if(stat(this->FileName, &fs) != 0)
    {
    vtkErrorMacro("Error opening file " << this->FileName);
    return 1;
    }
  if ( !this->DOM )
    {
    this->DOM = new XdmfDOM();
    }
    // this->DOM->GlobalDebugOn();
  if ( !this->Internals->DataItem )
    {
    this->Internals->DataItem = new XdmfDataItem();
    this->Internals->DataItem->SetDOM(this->DOM);
    }
  CurrentFileName = this->DOM->GetInputFileName();
  
  // Don't Re-Parse : Could Reset Parameters
  if((CurrentFileName == NULL) ||
     (STRCASECMP(CurrentFileName, this->FileName) != 0 ))
    {
    this->DOM->SetInputFileName(this->FileName);
    vtkDebugMacro( << ".!!............Preparing to Parse " << this->FileName);
    this->DOM->Parse();
    // Added By Jerry Clarke
    this->GridsModified = 1;
    }

  XdmfXmlNode domain = 0;
  int done = 0;
  this->Internals->DomainList.erase(this->Internals->DomainList.begin(),
    this->Internals->DomainList.end());
  for ( cc = 0; !done; cc ++ )
    {
    ostrstream str1, str2;
    domain = this->DOM->FindElement("Domain", cc);
    if ( !domain )
      {
      break;
      }
    XdmfConstString domainName = this->DOM->Get( domain, "Name" );
    ostrstream str;
    if ( !domainName )
      {
      str << "Domain" << cc << ends;
      domainName = str.str();
      }
    this->Internals->DomainList.push_back(domainName);
    str.rdbuf()->freeze(0);
    }
  if ( this->DomainName )
    {
    for ( cc = 0; !done; cc ++ )
      {
      domain = this->DOM->FindElement("Domain", cc);
      if ( !domain )
        {
        break;
        }
      XdmfConstString domainName = this->DOM->Get( domain, "Name" );
      ostrstream str;
      if ( !domainName )
        {
        str << "Domain" << cc << ends;
        domainName = str.str();
        }
      if( domainName && strcmp(domainName, this->DomainName) == 0)
        {
        str.rdbuf()->freeze(0);
        break;
        }      
      str.rdbuf()->freeze(0);
      }
    }

  if ( !domain )
    {
    domain = this->DOM->FindElement("Domain", 0); // 0 - domain index    
    }

  if ( !domain )
    {
    vtkErrorMacro("Cannot find any domain...");
    return 1;
    }

  this->Internals->DomainPtr = domain;
//  this->UpdateGrids();

  char* filename = 0;
  if ( this->FileName )
    {
    filename = new char[ strlen(this->FileName)+ 1];
    strcpy(filename, this->FileName);
    }
  int len = static_cast<int>(strlen(filename));
  for ( cc = len-1; cc >= 0; cc -- )
    {
    if ( filename[cc] != '/' && filename[cc] != '\\' )
      {
      filename[cc] = 0;
      }
    else
      {
      break;
      }
    }
  if ( filename[0] == 0 )
    {
    char buffer[1024];
    if ( GETCWD(buffer, 1023) )
      {
      delete [] filename;
      if ( buffer )
        {
        filename = new char[ strlen(buffer)+ 1];
        strcpy(filename, buffer);
        }
      }
    }
  this->DOM->SetWorkingDirectory(filename);
  delete [] filename;
////////////////
  
  
  this->UpdateGrids();

/*  
  if(1 !=this->GetNumberOfOutputPorts())
    {
    this->SetNumberOfOutputPorts(1);

    // We have to refresh the outputVector with this new number of ports.
    // The one in argument was given by the executive.
    // The only way to do that, is to ask the output vector from the executive.
    newOutputVector=this->GetExecutive()->GetOutputInformation();
    }
*/
  
  // for each output
  vtkXdmfReaderInternal::MapOfActualGrids::iterator currentGridIterator;
  currentGridIterator = this->Internals->ActualGrids.begin();
  int i=0;
  
  int someOutputChanged=0;
  vtkInformation *jinfo=newOutputVector->GetInformationObject(0);
  vtkMultiGroupDataSet *output=vtkMultiGroupDataSet::SafeDownCast(jinfo->Get(vtkDataObject::DATA_OBJECT()));
  if(output==0)
    {
    someOutputChanged=1;
    output=vtkMultiGroupDataSet::New();
    output->SetPipelineInformation(jinfo);
    output->Delete();
    }
    // Collapse on second level of Hierarchy
  output->SetNumberOfGroups(this->NumberOfEnabledActualGrids);
  //
  //while(currentGridIterator != this->Internals->ActualGrids.end())
  if(someOutputChanged)
    {
    this->GetPointDataArraySelection()->RemoveAllArrays();
    this->GetCellDataArraySelection()->RemoveAllArrays();
    }
  
  return 1;
}

//----------------------------------------------------------------------------
void vtkXdmfReader::SetGridName(const char* grid)
{
  this->DisableAllGrids();
  this->EnableGrid(grid);
}

//----------------------------------------------------------------------------
vtkDataObject *vtkXdmfReader::GetOutput()
{
  return this->GetOutput(0);
}

//----------------------------------------------------------------------------
void vtkXdmfReader::SetStride(int x, int y, int z)
{
  if ( x <= 0 || y <= 0 || z <= 0 )
    {
    vtkErrorMacro("Strides have to be greater than 0.");
    return;
    }
  vtkDebugMacro(<< this->GetClassName() << " (" << this 
    << "): setting Stride to (" << x << "," << y << "," << z << ")");
  if ((this->Stride[0] != x)||(this->Stride[1] != y)||(this->Stride[2] != z))
    {
    this->Stride[0] = x;
    this->Stride[1] = y;
    this->Stride[2] = z;
    this->Modified();
    }
}

//----------------------------------------------------------------------------
vtkDataObject *vtkXdmfReader::GetOutput(int idx)
{
  return this->GetExecutive()->GetOutputData(idx);
}

//----------------------------------------------------------------------------
int vtkXdmfReader::RequestData(
  vtkInformation *vtkNotUsed(request),
  vtkInformationVector **vtkNotUsed(inputVector),
  vtkInformationVector *outputVector)
{
  vtkDebugMacro("Execute");
  if ( !this->FileName )
    {
    vtkErrorMacro("File name not set");
    return 0;
    }

  if ( !this->DOM )
    {
    return 0;
    }
  
  int outputGrid = 0;
  vtkXdmfReaderInternal::MapOfActualGrids::iterator currentGridIterator;
  for ( currentGridIterator = this->Internals->ActualGrids.begin();
    currentGridIterator != this->Internals->ActualGrids.end();
    ++currentGridIterator )
    {
    if (currentGridIterator->second.Enabled )
      {
      this->Internals->RequestActualGridData(
        currentGridIterator->first.c_str(),&currentGridIterator->second,
        outputGrid, this->NumberOfEnabledActualGrids, outputVector);
      outputGrid++;
       // cout << "Progress " << 1.0 * outputGrid / this->NumberOfEnabledActualGrids << endl;
      this->UpdateProgress(1.0 * outputGrid / this->NumberOfEnabledActualGrids);
      }
    }
  
  return 1;
}

//----------------------------------------------------------------------------
int vtkXdmfReaderInternal::RequestActualGridData(
  const char* currentGridName,
  vtkXdmfReaderActualGrid* currentActualGrid,
  int outputGrid,
  int numberOfGrids,
  vtkInformationVector *outputVector)
{

  vtkInformation *info=outputVector->GetInformationObject(0);
  int procId=info->Get(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER());
  int nbProcs=info->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES());
  
  vtkInformation *outInfo = outputVector->GetInformationObject(0);
  vtkMultiGroupDataSet *mgd=vtkMultiGroupDataSet::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));

  if (!currentActualGrid->Collection)
    {
    return 0; //work around an undo/redo bug
    }

    unsigned int numberOfDataSets=currentActualGrid->Collection->Grids.size();
    
    currentActualGrid->Collection->UpdateCounts();
    int levels=currentActualGrid->Collection->GetNumberOfLevels();
    // mgd->SetNumberOfLevels(levels);
    // int levels = 1;
    
    int level=0;
    // cout << "mgd->SetNumberOfDataSets(" << outputGrid << "," << currentActualGrid->Collection->GetNumberOfDataSets(level) << ")" << endl;
    mgd->SetNumberOfDataSets(outputGrid, currentActualGrid->Collection->GetNumberOfDataSets(level));
    while(level<levels)
      {
      // mgd->SetNumberOfDataSets(level,currentActualGrid->Collection->GetNumberOfDataSets(level));
      ++level;
      }
 
   // cout << "level = " << level << " levels = " << levels << " numberOfDataSets = " << numberOfDataSets << endl; 
    vtkXdmfReaderGridCollection::SetOfGrids::iterator gridIt;
    vtkXdmfReaderGridCollection::SetOfGrids::iterator gridItEnd;
    
    int numBlocksPerProcess=numberOfDataSets/nbProcs;
    int leftOverBlocks=numberOfDataSets-(numBlocksPerProcess*nbProcs);
    
    int blockStart;
    int blockEnd;
    
    if(procId<leftOverBlocks)
      {
      blockStart=(numBlocksPerProcess+1)*procId;
      blockEnd=blockStart+(numBlocksPerProcess+1)-1;
      }
    else
      {
      blockStart=numBlocksPerProcess*procId+leftOverBlocks;
      blockEnd=blockStart+numBlocksPerProcess-1;
      }
    
    gridIt=currentActualGrid->Collection->Grids.begin();
    gridItEnd=currentActualGrid->Collection->Grids.end();
    int result=1;
    
    int datasetIdx=0;
    
    vtkMultiGroupDataInformation *compInfo=vtkMultiGroupDataInformation::SafeDownCast(info->Get(vtkCompositeDataPipeline::COMPOSITE_DATA_INFORMATION()));
    
    // currentIndex in each level.
    vtkstd::vector<int> currentIndex(levels);
    level=0;
    while(level<levels)
      {
      currentIndex[level]=0;
      ++level;
      }
    
    while(gridIt != gridItEnd && result)
      {
      vtkXdmfReaderGrid *subgrid=gridIt->second;
      level=subgrid->Level;
      int index=currentIndex[level];
      if(datasetIdx<blockStart || datasetIdx>blockEnd)
        {
             // cout << " mgd->SetDataSet(0) for index " << index << " outputGrid " << outputGrid << " datasetIdx " << datasetIdx << endl;
        // mgd->SetDataSet(level,index,0); // empty, on another processor
        mgd->SetDataSet(outputGrid, index, 0); // empty, on another processor
        }
      else
        {
        XdmfGrid *xdmfGrid=subgrid->XMGrid;
        if( xdmfGrid->GetTopology()->GetClass() == XDMF_UNSTRUCTURED ) 
          {
          vtkUnstructuredGrid *ds=vtkUnstructuredGrid::New();
          ds->SetMaximumNumberOfPieces(1);
             // cout << " 1 mgd->SetDataSet(ds) for index " << index << " outputGrid " << outputGrid << " datasetIdx " << datasetIdx << endl;
          // mgd->SetDataSet(level,index,ds);
          mgd->SetDataSet(outputGrid, index, ds);
          ds->Delete();
          } 
        else if( xdmfGrid->GetTopology()->GetTopologyType() == XDMF_2DSMESH ||
                 xdmfGrid->GetTopology()->GetTopologyType() == XDMF_3DSMESH )
          {
          vtkStructuredGrid *ds=vtkStructuredGrid::New();
             // cout << " 2 mgd->SetDataSet(ds) for index " << index << " level " << level << " datasetIdx " << datasetIdx << endl;
          // mgd->SetDataSet(level,index,ds);
          mgd->SetDataSet(outputGrid, index, ds);
          ds->Delete();
          }
        else if ( xdmfGrid->GetTopology()->GetTopologyType() == XDMF_2DCORECTMESH ||
                  xdmfGrid->GetTopology()->GetTopologyType() == XDMF_3DCORECTMESH )
          {
#ifdef USE_IMAGE_DATA
            vtkImageData *ds=vtkImageData::New();
#else
            vtkUniformGrid *ds=vtkUniformGrid::New();
#endif
            mgd->SetDataSet(level,index,ds);
            ds->Delete();
          }
        else if ( xdmfGrid->GetTopology()->GetTopologyType() == XDMF_2DRECTMESH ||
                  xdmfGrid->GetTopology()->GetTopologyType() == XDMF_3DRECTMESH )
          {
            vtkRectilinearGrid *ds=vtkRectilinearGrid::New();
            mgd->SetDataSet(level,index,ds);
            ds->Delete();
          }
        else
          {
          // Unknown type for this sub grid. 
          return 0;
          }
        // vtkDataObject *ds=mgd->GetDataSet(level,index);
        vtkDataObject *ds=mgd->GetDataSet(outputGrid,index);
        vtkInformation *subInfo=compInfo->GetInformation(level,index);
        result=this->RequestSingleGridData("",gridIt->second,subInfo,ds,1);
        // cout << "ds level " << level << " index " << index << " = " << endl;
        // ds->Print(cout);
        }
      ++currentIndex[level];
      ++gridIt;
      ++datasetIdx;
      // cout << "Progress " << 1.0 * datasetIdx / numberOfDataSets << endl;
      this->Reader->UpdateProgress(1.0 * datasetIdx / numberOfDataSets);
      }
      // mgd->Print(cout);
      // mgd->DebugOn();
    return result;
}

//----------------------------------------------------------------------------
int vtkXdmfReaderInternal::RequestSingleGridData(
  const char* currentGridName,
  vtkXdmfReaderGrid *grid,
  vtkInformation *outInfo,
  vtkDataObject *output,
  int isSubBlock)
{
  int *readerStride = this->Reader->GetStride();

  vtkDataArraySelection* pointDataArraySelection = this->Reader->GetPointDataArraySelection();
  vtkDataArraySelection* cellDataArraySelection = this->Reader->GetCellDataArraySelection();

  // Handle single grid
  XdmfGrid* xdmfGrid = grid->XMGrid;
  XdmfDOM* xdmfDOM = xdmfGrid->GetDOM();
  
  if ( !grid->DataDescription )
    {
    grid->DataDescription = xdmfGrid->GetTopology()->GetShapeDesc();
    //continue;
    }
  
  vtkDebugWithObjectMacro( this->Reader, << "Reading Heavy Data for " << xdmfGrid->GetName());
  xdmfGrid->Update();
  
  // True for all 3d datasets except unstructured grids
  int globalrank = 3;
  switch(xdmfGrid->GetTopology()->GetTopologyType())
    {
    case XDMF_2DSMESH: case XDMF_2DRECTMESH: case XDMF_2DCORECTMESH:
      globalrank = 2;
    }
  if ( xdmfGrid->GetTopology()->GetClass() == XDMF_UNSTRUCTURED )
    {
    globalrank = 1;
    }
  
  
  vtkIdType cc;
  // XdmfXmlNode attrNode = this->DOM->FindElement("Attribute");
  // XdmfXmlNode dataNode = this->DOM->FindElement("DataStructure", 0, attrNode);
  XdmfXmlNode attrNode;
  XdmfXmlNode dataNode;
  // XdmfInt64 start[3]  = { 0, 0, 0 };
  // XdmfInt64 stride[3] = { 1, 1, 1 };
  // XdmfInt64 count[3] = { 0, 0, 0 };
  // XdmfInt64 end[3] = { 0, 0, 0 };
  // Must Account for Vectors and Tensors
  XdmfInt64 start[4]  = { 0, 0, 0, 0 };
  XdmfInt64 stride[4] = { 1, 1, 1, 1 };
  XdmfInt64 count[4] = { 0, 0, 0, 0 };
  // XdmfInt64 end[4] = { 0, 0, 0, 0 };
  grid->DataDescription->GetShape(count);

  int upext[6];
  int whext[6];
  
  if( xdmfGrid->GetTopology()->GetClass() != XDMF_UNSTRUCTURED)
    {
    if(isSubBlock)
      {
      // the composite pipeline does not set update extent, so just take
      // the whole extent for as update extent.
      outInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), upext);
      }
    else
      {
      outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), upext);
      }
    
    outInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), whext);
    
    start[2] = vtkMAX(0, upext[0]);
    start[1] = vtkMAX(0, upext[2]);
    start[0] = vtkMAX(0, upext[4]);
    
    count[2] = upext[1] - upext[0];
    count[1] = upext[3] - upext[2];
    count[0] = upext[5] - upext[4];
    }

  XdmfGeometry  *Geometry = xdmfGrid->GetGeometry();
    
    
  // Read Topology for Unstructured Grid
  if( xdmfGrid->GetTopology()->GetClass() == XDMF_UNSTRUCTURED ) 
    {
    vtkUnstructuredGrid  *vGrid = static_cast<vtkUnstructuredGrid *>(output);
    vtkCellArray                *verts;
    XdmfInt32           vType;
    XdmfInt32           NodesPerElement;
    vtkIdType           NumberOfElements;
    vtkIdType           i, j, index;    
    XdmfInt64           Length, *Connections;
    vtkIdType           *connections;
    int                 *cell_types, *ctp;
    
    vtkDebugWithObjectMacro(this->Reader,<< "Unstructured Topology is " 
                            << xdmfGrid->GetTopology()->GetTopologyTypeAsString() );
    switch ( xdmfGrid->GetTopology()->GetTopologyType() )
      {
      case  XDMF_POLYVERTEX :
        vType = VTK_POLY_VERTEX;
        break;
      case  XDMF_POLYLINE :
        vType = VTK_POLY_LINE;
        break;
      case  XDMF_POLYGON :
        vType = VTK_POLYGON;
        break;
      case  XDMF_TRI :
        vType = VTK_TRIANGLE;
        break;
      case  XDMF_QUAD :
        vType = VTK_QUAD;
        break;
      case  XDMF_TET :
        vType = VTK_TETRA;
        break;
      case  XDMF_PYRAMID :
        vType = VTK_PYRAMID;
        break;
      case  XDMF_WEDGE :
        vType = VTK_WEDGE;
        break;
      case  XDMF_HEX :
        vType = VTK_HEXAHEDRON;
        break;
      case  XDMF_EDGE_3 :
        vType = VTK_QUADRATIC_EDGE ;
        break;
      case  XDMF_TRI_6 :
        vType = VTK_QUADRATIC_TRIANGLE ;
        break;
      case  XDMF_QUAD_8 :
        vType = VTK_QUADRATIC_QUAD ;
        break;
      case  XDMF_TET_10 :
        vType = VTK_QUADRATIC_TETRA ;
        break;
      case  XDMF_PYRAMID_13 :
        vType = VTK_QUADRATIC_PYRAMID ;
        break;
      case  XDMF_WEDGE_15 :
        vType = VTK_QUADRATIC_WEDGE ;
        break;
      case  XDMF_HEX_20 :
        vType = VTK_QUADRATIC_HEXAHEDRON ;
        break;
      case XDMF_MIXED :
        vType = -1;
        break;
      default :
        XdmfErrorMessage("Unknown Topology Type = " << xdmfGrid->GetTopology()->GetTopologyType());
        return 1;
      }
      if( xdmfGrid->GetTopology()->GetTopologyType() != XDMF_MIXED){
        NodesPerElement = xdmfGrid->GetTopology()->GetNodesPerElement();
        if ( xdmfGrid->GetTopology()->GetConnectivity()->GetRank() == 2 )
          {
          NodesPerElement = xdmfGrid->GetTopology()->GetConnectivity()->GetDimension(1);
          }
    
        /* Create Cell Type Array */
        Length = xdmfGrid->GetTopology()->GetConnectivity()->GetNumberOfElements();
        Connections = new XdmfInt64[ Length ];
        xdmfGrid->GetTopology()->GetConnectivity()->GetValues( 0, Connections, Length );
    
        NumberOfElements = xdmfGrid->GetTopology()->GetShapeDesc()->GetNumberOfElements();
        ctp = cell_types = new int[ NumberOfElements ];
    
        /* Create Cell Array */
        verts = vtkCellArray::New();
    
        /* Get the pointer */
        connections = verts->WritePointer(
          NumberOfElements,
          NumberOfElements * ( 1 + NodesPerElement ));
    
        /* Connections : N p1 p2 ... pN */
        /* i.e. Triangles : 3 0 1 2    3 3 4 5   3 6 7 8 */
        index = 0;
        for( j = 0 ; j < NumberOfElements; j++ )
          {
          *ctp++ = vType;
          *connections++ = NodesPerElement;
          for( i = 0 ; i < NodesPerElement; i++ )
            {
            *connections++ = Connections[index++];
            }
          }
    } else {
        // Mixed Topology
        /* Create Cell Type Array */
        vtkIdTypeArray *IdArray;
        vtkIdType RealSize;

        Length = xdmfGrid->GetTopology()->GetConnectivity()->GetNumberOfElements();
        Connections = new XdmfInt64[ Length ];
        xdmfGrid->GetTopology()->GetConnectivity()->GetValues( 0, Connections, Length );
        NumberOfElements = xdmfGrid->GetTopology()->GetShapeDesc()->GetNumberOfElements();
        ctp = cell_types = new int[ NumberOfElements ];
    
        /* Create Cell Array */
        verts = vtkCellArray::New();
    
        /* Get the pointer. Make it Big enough ... too big for now */
        // cout << "::::: Length = " << Length << endl;
        connections = verts->WritePointer(
          NumberOfElements,
          Length);
        //   Length);
        /* Connections : N p1 p2 ... pN */
        /* i.e. Triangles : 3 0 1 2    3 3 4 5   3 6 7 8 */
        index = 0;
        for( j = 0 ; j < NumberOfElements; j++ )
          {
            switch ( Connections[index++] )
              {
              case  XDMF_POLYVERTEX :
                vType = VTK_POLY_VERTEX;
                NodesPerElement = Connections[index++];
                break;
              case  XDMF_POLYLINE :
                vType = VTK_POLY_LINE;
                NodesPerElement = Connections[index++];
                break;
              case  XDMF_POLYGON :
                vType = VTK_POLYGON;
                NodesPerElement = Connections[index++];
                break;
              case  XDMF_TRI :
                vType = VTK_TRIANGLE;
                NodesPerElement = 3;
                break;
              case  XDMF_QUAD :
                vType = VTK_QUAD;
                NodesPerElement = 4;
                break;
              case  XDMF_TET :
                vType = VTK_TETRA;
                NodesPerElement = 4;
                break;
              case  XDMF_PYRAMID :
                vType = VTK_PYRAMID;
                NodesPerElement = 5;
                break;
              case  XDMF_WEDGE :
                vType = VTK_WEDGE;
                NodesPerElement = 6;
                break;
              case  XDMF_HEX :
                vType = VTK_HEXAHEDRON;
                NodesPerElement = 8;
                break;
              case  XDMF_EDGE_3 :
                vType = VTK_QUADRATIC_EDGE ;
                NodesPerElement = 3;
                break;
              case  XDMF_TRI_6 :
                vType = VTK_QUADRATIC_TRIANGLE ;
                NodesPerElement = 6;
                break;
              case  XDMF_QUAD_8 :
                vType = VTK_QUADRATIC_QUAD ;
                NodesPerElement = 8;
                break;
              case  XDMF_TET_10 :
                vType = VTK_QUADRATIC_TETRA ;
                NodesPerElement = 10;
                break;
              case  XDMF_PYRAMID_13 :
                vType = VTK_QUADRATIC_PYRAMID ;
                NodesPerElement = 13;
                break;
              case  XDMF_WEDGE_15 :
                vType = VTK_QUADRATIC_WEDGE ;
                NodesPerElement = 15;
                break;
              case  XDMF_HEX_20 :
                vType = VTK_QUADRATIC_HEXAHEDRON ;
                NodesPerElement = 20;
                break;
              default :
                XdmfErrorMessage("Unknown Topology Type");
                return 1;
              }
          *ctp++ = vType;
          *connections++ = NodesPerElement;
          for( i = 0 ; i < NodesPerElement; i++ )
            {
            *connections++ = Connections[index++];
            }
          }
    // Resize the Array to the Proper Size
    IdArray = verts->GetData();
    RealSize = index - 1;
    vtkDebugWithObjectMacro(this->Reader, "Resizing to " << RealSize << " elements");
    IdArray->Resize(RealSize);
    }
    delete [] Connections;
    vGrid->SetCells(cell_types, verts);
    /* OK, because of reference counting */
    verts->Delete();
    delete [] cell_types;
    vGrid->Modified();
    }  // if( xdmfGrid->GetClass() == XDMF_UNSTRUCTURED ) 
  else if( xdmfGrid->GetTopology()->GetTopologyType() == XDMF_2DSMESH ||
           xdmfGrid->GetTopology()->GetTopologyType() == XDMF_3DSMESH )
    {
    vtkDebugWithObjectMacro(this->Reader, << "Setting Extents for vtkStructuredGrid");
    vtkStructuredGrid  *vGrid = vtkStructuredGrid::SafeDownCast(output);
    vGrid->SetExtent(upext);    
    } 
  else if (xdmfGrid->GetTopology()->GetTopologyType() == XDMF_2DCORECTMESH ||
           xdmfGrid->GetTopology()->GetTopologyType() == XDMF_3DCORECTMESH ) 
    {
#ifdef USE_IMAGE_DATA
    vtkImageData *idata = vtkImageData::SafeDownCast(output);
#else
    vtkUniformGrid *idata = vtkUniformGrid::SafeDownCast(output);
#endif
    idata->SetExtent(upext);
    }
  else if ( xdmfGrid->GetTopology()->GetTopologyType() == XDMF_2DRECTMESH ||
            xdmfGrid->GetTopology()->GetTopologyType() == XDMF_3DRECTMESH )
    {
    vtkRectilinearGrid *vGrid = vtkRectilinearGrid::SafeDownCast(output);
    vGrid->SetExtent(upext);    
    }
  else
    {
    vtkErrorWithObjectMacro(this->Reader,"Do not understand topology type: " 
                            << xdmfGrid->GetTopology()->GetTopologyType());
    }
  // Read Geometry
  if( ( Geometry->GetGeometryType() == XDMF_GEOMETRY_X_Y_Z ) ||
      ( Geometry->GetGeometryType() == XDMF_GEOMETRY_XYZ ) ||
      ( Geometry->GetGeometryType() == XDMF_GEOMETRY_X_Y ) ||
      ( Geometry->GetGeometryType() == XDMF_GEOMETRY_XY ) )
    {
    XdmfInt64   Length;
    vtkPoints   *Points;
    vtkPointSet *Pointset = vtkPointSet::SafeDownCast(output);
    
    // Special flag, for structured data
    int structured_grid = 0;
    if ( vtkStructuredGrid::SafeDownCast(output) )
      {
      structured_grid = 1;
      }

    Points = Pointset->GetPoints();
    if( !Points )
      {
      vtkDebugWithObjectMacro(this->Reader,<<  "Creating vtkPoints" );
      Points = vtkPoints::New();
      Pointset->SetPoints( Points );
      // OK Because of Reference Counting
      Points->Delete();
      }

    if( Geometry->GetPoints() )
      {
      if( Points )
        {
        if ( Geometry->GetPoints()->GetNumberType() == XDMF_FLOAT32_TYPE )
          {
          if ( Points->GetData()->GetDataType() != VTK_FLOAT)
            {
            vtkFloatArray* da = vtkFloatArray::New();
            da->SetNumberOfComponents(3);
            Points->SetData(da);
            da->Delete();
            }
          }
        else // means == XDMF_FLOAT64_TYPE
          {
          if ( Points->GetData()->GetDataType() != VTK_DOUBLE )
            {