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

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

#include "vtkAMRBox.h"
#include "vtkBoundingBox.h"
#include "vtkByteSwap.h"
#include "vtkCallbackCommand.h"
#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkCompositeDataIterator.h"
#include "vtkCompositeDataPipeline.h"
#include "vtkDataArraySelection.h"
#include "vtkDoubleArray.h"
#include "vtkFloatArray.h"
//#include "vtkHierarchicalBoxDataSet.h"
#include "vtkNonOverlappingAMR.h"
#include "vtkImageData.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkIntArray.h"
#include "vtkMultiBlockDataSet.h"
#include "vtkMultiProcessController.h"
#include "vtkMultiProcessStream.h"
#include "vtkObjectFactory.h"
#include "vtkProcessGroup.h"
#include "vtkPolyData.h"
#include "vtkPointData.h"
#include "vtkPoints.h"
#include "vtkRectilinearGrid.h"
#include "vtkSmartPointer.h"
#include "vtkUniformGrid.h"
#include "vtkUnsignedCharArray.h"

#include "vtkSpyPlotReaderMap.h"
#include "vtkSpyPlotUniReader.h"
#include "vtkSpyPlotBlock.h"
#include "vtkSpyPlotBlockIterator.h"
#include "vtkSpyPlotIStream.h"

#include <map>
#include <set>
#include <vector>
#include <string>
#include <sys/stat.h>
#include <vtksys/SystemTools.hxx>
#include <assert.h>
#include <cctype>

#define vtkMIN(x, y) \
  (\
   ((x)<(y))?(x):(y)\
   )

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

// #define PARAVIEW_ENABLE_SPYPLOT_MARKERS

vtkStandardNewMacro(vtkSpyPlotReader);
vtkCxxSetObjectMacro(vtkSpyPlotReader,GlobalController,vtkMultiProcessController);

static void createSpyPlotLevelArray(vtkCellData *cd, int size, int level);

class vtkSpyPlotReader::VectorOfDoubles : public std::vector<double> {};

//-----------------------------------------------------------------------------
vtkSpyPlotReader::vtkSpyPlotReader()
{
  this->SetNumberOfInputPorts(0);
  this->SetNumberOfOutputPorts(2);
#ifdef PARAVIEW_ENABLE_SPYPLOT_MARKERS
  this->SetNumberOfOutputPorts(3);
#endif // PARAVIEW_ENABLE_SPYPLOT_MARKERS

  this->Map = new vtkSpyPlotReaderMap;
  this->Bounds = new vtkBoundingBox;
  this->BoxSize[0] = -1;
  this->BoxSize[1] = -1;
  this->BoxSize[2] = -1;
  this->FileName = 0;
  this->CellDataArraySelection = vtkDataArraySelection::New();
  this->TimeStep = 0;
  this->TimeStepRange[0] = 0;
  this->TimeStepRange[1] = 0;
  this->ComputeDerivedVariables = 1;
  this->DownConvertVolumeFraction = 1;
  this->MergeXYZComponents = 1;

  // this has all of the processes.
  this->GlobalController = 0;
  this->SetGlobalController(vtkMultiProcessController::GetGlobalController());

  this->DistributeFiles=0; // by default, distribute blocks, not files.
  this->GenerateLevelArray=0; // by default, do not generate level array.
  this->GenerateBlockIdArray=0; // by default, do not generate block id array.
  this->GenerateActiveBlockArray = 0; // by default do not generate active array
  this->GenerateTracerArray = 1; // by default do not generate tracer array
  this->GenerateMarkers = 1; // by default do generate markers
  this->IsAMR = 1;
  this->FileNameChanged = true;
  this->TimeSteps = new vtkSpyPlotReader::VectorOfDoubles();
  this->TimeRequestedFromPipeline = false;
}

//-----------------------------------------------------------------------------
vtkSpyPlotReader::~vtkSpyPlotReader()
{
  this->SetFileName(0);
  this->CellDataArraySelection->Delete();
  this->Map->Clean(0);
  delete this->Map;
  delete this->Bounds;
  this->Map = 0;
  this->SetGlobalController(0);
  delete this->TimeSteps;
}


//-----------------------------------------------------------------------------
void vtkSpyPlotReader::SetTimeStepsInternal(const vtkSpyPlotReader::VectorOfDoubles& doubles)
{
  *this->TimeSteps = doubles;
  this->TimeStepRange[1] = static_cast<int>(doubles.size()-1);
}

//-----------------------------------------------------------------------------
// Create either vtkHierarchicalBoxDataSet or vtkMultiBlockDataSet based on
// whether the dataset is AMR.
void vtkSpyPlotReader::SetFileName(const char* filename)
{
  if ( this->FileName == NULL && filename == NULL) { return;}
  if ( this->FileName && filename && (!strcmp(this->FileName,filename))) { return;}

  // Filename is going to change.
  delete [] this->FileName;
  this->FileName = vtksys::SystemTools::DuplicateString(filename);
  this->FileNameChanged = true;
  this->Modified();
}

//-----------------------------------------------------------------------------
int vtkSpyPlotReader::FillOutputPortInformation(int port, vtkInformation *info)
{
  switch (port)
    {
    case 0:
    case 2:
      info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkCompositeDataSet");
      break;
    case 1:
      info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkPolyData");
      break;
    default:
      return 0;
    }

  return 1;
}

//-----------------------------------------------------------------------------
// Create either vtkHierarchicalBoxDataSet or vtkMultiBlockDataSet based on
// whether the dataset is AMR.
int vtkSpyPlotReader::RequestDataObject(vtkInformation *req,
  vtkInformationVector **vtkNotUsed(inV),
  vtkInformationVector *outV)
{
  vtkInformation *outInfo = outV->GetInformationObject(0);
  vtkCompositeDataSet* outData = NULL;

  // Call UpdateFile (which sets IsAMR) because RequestInformation isn't called
  // before RequestDataObject
  this->UpdateFile (req, outV);

  if (this->IsAMR)
    {
//    outData = vtkHierarchicalBoxDataSet::New();
    outData = vtkNonOverlappingAMR::New();
    }
  else
    {
    outData = vtkMultiBlockDataSet::New();
    }

  outInfo->Set(vtkDataObject::DATA_EXTENT_TYPE(), outData->GetExtentType());
  outInfo->Set(vtkDataObject::DATA_OBJECT(), outData);
  outData->Delete();

  outInfo = outV->GetInformationObject(1);
  vtkPolyData* polyData = vtkPolyData::New ();
  outInfo->Set(vtkDataObject::DATA_EXTENT_TYPE(), polyData->GetExtentType());
  outInfo->Set(vtkDataObject::DATA_OBJECT(), polyData);
  polyData->Delete ();

#ifdef PARAVIEW_ENABLE_SPYPLOT_MARKERS
  outInfo = outV->GetInformationObject(2);
  vtkMultiBlockDataSet* data = vtkMultiBlockDataSet::New ();
  outInfo->Set(vtkDataObject::DATA_EXTENT_TYPE(), data->GetExtentType());
  outInfo->Set(vtkDataObject::DATA_OBJECT(), data);
  data->Delete ();
#endif // PARAVIEW_ENABLE_SPYPLOT_MARKERS

  return 1;
}

//-----------------------------------------------------------------------------
// Read the case file and the first binary file do get meta
// informations (number of files, number of fields, number of timestep).
int vtkSpyPlotReader::RequestInformation(vtkInformation *request,
                                         vtkInformationVector **inputVector,
                                         vtkInformationVector *outputVector)
{
  if(this->GlobalController==0)
    {
    vtkErrorMacro("Controller not specified. This reader requires controller to be set.");
    //return 0;
    }
  if(!this->Superclass::RequestInformation(request,inputVector,outputVector))
    {
    return 0;
    }


  struct stat fs;
  if(stat(this->FileName,&fs)!=0)
    {
    vtkErrorMacro("Cannot find file " << this->FileName);
    return 0;
    }

  // I don't think UpdateFile is needed here since RequestDataObject() would have
  // got UpdateFile to do the necessary work.
  if (!this->UpdateFile (request, outputVector))
    {
    return 0;
    }

  vtkInformation* outInfo0 = outputVector->GetInformationObject(0);
  outInfo0->Remove(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
  outInfo0->Remove(vtkStreamingDemandDrivenPipeline::TIME_RANGE());
  outInfo0->Set(CAN_HANDLE_PIECE_REQUEST(), 1);

  vtkInformation* outInfo1 = outputVector->GetInformationObject(1);
  outInfo1->Remove(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
  outInfo1->Remove(vtkStreamingDemandDrivenPipeline::TIME_RANGE());
  outInfo1->Set(CAN_HANDLE_PIECE_REQUEST(), 1);
#ifdef PARAVIEW_ENABLE_SPYPLOT_MARKERS
  vtkInformation* outInfo2 = 0;
  if ( this->GenerateMarkers )
    {
    outInfo2 = outputVector->GetInformationObject(2);
    outInfo2->Remove(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
    outInfo2->Remove(vtkStreamingDemandDrivenPipeline::TIME_RANGE());
    outInfo2->Set(CAN_HANDLE_PIECE_REQUEST(), 1);
    }
#endif // PARAVIEW_ENABLE_SPYPLOT_MARKERS
  if (this->TimeSteps->size() > 0)
    {
    outInfo0->Set(vtkStreamingDemandDrivenPipeline::TIME_STEPS(),
      &(*this->TimeSteps)[0],
      static_cast<int>(this->TimeSteps->size()));
    double timeRange[2];
    timeRange[0] = this->TimeSteps->front();
    timeRange[1] = this->TimeSteps->back();
    outInfo0->Set(vtkStreamingDemandDrivenPipeline::TIME_RANGE(),
      timeRange, 2);

    outInfo1->Set(vtkStreamingDemandDrivenPipeline::TIME_STEPS(),
      &(*this->TimeSteps)[0],
      static_cast<int>(this->TimeSteps->size()));

    outInfo1->Set(vtkStreamingDemandDrivenPipeline::TIME_RANGE(),
      timeRange, 2);

#ifdef PARAVIEW_ENABLE_SPYPLOT_MARKERS
    if ( this->GenerateMarkers) 
      {
      outInfo2->Set(vtkStreamingDemandDrivenPipeline::TIME_STEPS(),
        &(*this->TimeSteps)[0],
        static_cast<int>(this->TimeSteps->size()));

      outInfo2->Set(vtkStreamingDemandDrivenPipeline::TIME_RANGE(),
        timeRange, 2);
      }
#endif // PARAVIEW_ENABLE_SPYPLOT_MARKERS
    }
  return 1;
}

//-----------------------------------------------------------------------------
int vtkSpyPlotReader::UpdateFile (vtkInformation* request,
                                  vtkInformationVector* outputVector)
{
  if (!this->FileNameChanged)
    {
    // If filename hasn't changed since the last time we read the core
    // meta-data, we don't need to re-read it. If the file was processed
    // successfully, then this->Map->Files won't be empty.
    return (this->Map->Files.size() != 0)? 1 : 0;
    }

  this->FileNameChanged = false;

  const int procId = this->GlobalController?
    this->GlobalController->GetLocalProcessId() : 0;
  const int numProcs = this->GlobalController?
    this->GlobalController->GetNumberOfProcesses() : 1;
  //cout << procId << " : " << __LINE__ << endl;

  // When running in parallel, we need to ensure that all testing and meta-data
  // loading only happens on the root node. (BUG #12720).
  // We fill up the "Map" on root node, and then share the filename with all
  // nodes.
  if (procId == 0)
    {
    // Clean Map and initialize it with the given file.
    this->Map->Initialize(this->FileName);
    }
  if (numProcs > 1)
    {
    vtkMultiProcessStream stream;
    this->Map->Save(stream);
    this->GlobalController->Broadcast(stream, 0);
    if (procId > 0)
      {
      this->Map->Load(stream);
      }
    }

  return this->Map->Files.size() > 0?
    this->UpdateMetaData(request, outputVector) : 0;
}

//-----------------------------------------------------------------------------
int vtkSpyPlotReader::UpdateMetaData(
  vtkInformation* vtkNotUsed(request),
  vtkInformationVector* vtkNotUsed(outputVector))
{
  if (this->Map->Files.size() == 0)
    {
    vtkErrorMacro("The internal file map is empty!");
    return 0;
    }

  // This method needs to collect meta-data i.e.
  // number of timesteps and names for the cell-arrays.
  // For that we use the first file specified. To avoid IO issues in parallel,
  // this meta-data is collected on root node alone and then broadcast to all
  // processes (BUG #12720).
  const int procId = this->GlobalController?
    this->GlobalController->GetLocalProcessId() : 0;
  const int numProcs = this->GlobalController?
    this->GlobalController->GetNumberOfProcesses() : 1;
  //cout << procId << " : " << __LINE__ << endl;

  vtkSpyPlotReader::VectorOfDoubles timesteps;

  if (procId == 0)
    {
    vtkSpyPlotReaderMap::MapOfStringToSPCTH::iterator iter =
      this->Map->Files.begin();
    assert(iter!=this->Map->Files.end());

    vtkSpyPlotUniReader *uniReader = this->Map->GetReader(iter, this);
    // Open the file and get the number time steps and arrays.
    if (uniReader->ReadInformation() == 0) 
      {
      return 0;
      }
    uniReader->GetTimeStepRange(this->TimeStepRange);
    int num_timesteps  = this->TimeStepRange[1] + 1;
    timesteps.insert(timesteps.end(),
      uniReader->GetTimeArray(), uniReader->GetTimeArray() + num_timesteps);
    this->IsAMR = uniReader->IsAMR();
    }

  if (numProcs > 1)
    {
    vtkMultiProcessStream stream;
    stream << this->IsAMR << static_cast<int>(timesteps.size());
    for (size_t cc = 0; cc < timesteps.size(); cc++)
      {
      stream << timesteps[cc];
      }
    this->GlobalController->Broadcast(stream, 0);

    if (procId > 0)
      {
      int size;
      stream  >> this->IsAMR >> size;
      timesteps.resize(size);
      for (int cc=0; cc < size; cc++)
        {
        double val; stream >> val;
        timesteps[cc] = val;
        }
      }
    }
  this->SetTimeStepsInternal(timesteps);

  // I am not bothering syncing array names between processes since all enabled
  // arrays will be marked explicitly on all processes any ways.

  // To minimize the need for interprocess communication when timesteps change,
  // I am removing the code that updated the cell-arrays available per timestep.
  // The code seemed to assume that the arrays available could change over time.
  // That's sounds highly improbable.
  return 1;
}

// Magic number that encode the message ids for parallel communication
enum
{
  VTK_MSG_SPY_READER_HAS_BOUNDS=288302,
  VTK_MSG_SPY_READER_LOCAL_BOUNDS,
  VTK_MSG_SPY_READER_GLOBAL_BOUNDS,
  VTK_MSG_SPY_READER_LOCAL_NUMBER_OF_LEVELS,
  VTK_MSG_SPY_READER_GLOBAL_NUMBER_OF_LEVELS,
  VTK_MSG_SPY_READER_LOCAL_NUMBER_OF_DATASETS,
  VTK_MSG_SPY_READER_GLOBAL_NUMBER_OF_DATASETS,
  VTK_MSG_SPY_READER_GLOBAL_DATASETS_INDEX
};

template<class DataType>
int vtkSpyPlotRemoveBadGhostCells(
        DataType* dataType,
        vtkDataArray* dataArray,
        int realExtents[6],
        int realDims[3],
        int ptDims[3],
        int realPtDims[3])
{
  /*
    vtkDebugMacro( "//////////////////////////////////////////////////////////////////////////////////" );
    vtkDebugMacro( << dataArray << " vtkSpyPlotRemoveBadGhostCells(" << dataArray->GetName() << ")" );
    vtkDebugMacro( "DataArray: ");
    dataArray->Print(cout);
    vtkDebugMacro( "Real Extents: " << coutVector6(realExtents) );
    vtkDebugMacro( "Real Dims:    " << coutVector3(realDims) );
    vtkDebugMacro( "PT Dims:      " << coutVector3(ptDims) );
    vtkDebugMacro( "Real PT Dims: " << coutVector3(realPtDims) );
    vtkDebugMacro( "//////////////////////////////////////////////////////////////////////////////////" );
  */
  (void)* dataType;
  // skip some cell data.

  //Performance analysis has shown that using ComputeCellId to be a bottleneck of the reader
  //so we are going to replace the method with an incremental algorithm that will reduce the total
  //number of multiplications.
  vtkIdType kOffset[2]={0,0},jOffset[2]={0,0};
  vtkIdType realCellId=0,oldCellId=0;

  int xyz[3];
  int destXyz[3];
  DataType* dataPtr = static_cast<DataType*>(dataArray->GetVoidPointer(0));
  for (xyz[2] = realExtents[4], destXyz[2] = 0;
       xyz[2] < realExtents[5]; ++xyz[2], ++destXyz[2])
    {
    kOffset[0] = destXyz[2] * (realPtDims[1]-1);
    kOffset[1] = xyz[2] * (ptDims[1]-1);

    for (xyz[1] = realExtents[2], destXyz[1] = 0;
         xyz[1] < realExtents[3]; ++xyz[1], ++destXyz[1])
      {
      jOffset[0] = (kOffset[0] + destXyz[1]) * (realPtDims[0]-1);
      jOffset[1] = (kOffset[1] + xyz[1]) * (ptDims[0]-1);

      for (xyz[0] = realExtents[0], destXyz[0] = 0;
           xyz[0] < realExtents[1]; ++xyz[0], ++destXyz[0])
        {
        realCellId = jOffset[0] + destXyz[0];
        oldCellId = jOffset[1] + xyz[0];
        dataPtr[realCellId] = dataPtr[oldCellId];
        //old slow way of calculating cell id
        //dataPtr[vtkStructuredData::ComputeCellId(realPtDims,destXyz)] =
        //  dataPtr[vtkStructuredData::ComputeCellId(ptDims,xyz)];
        }
      }
    }
  dataArray->SetNumberOfTuples(realDims[0]*realDims[1]*realDims[2]);
  return 1;
}
//-----------------------------------------------------------------------------
int vtkSpyPlotReader::UpdateTimeStep(vtkInformation *requestInfo,
                                     vtkInformationVector *outputInfoVec,
                                     vtkCompositeDataSet *outputData)
{
  int port = requestInfo->Get (vtkStreamingDemandDrivenPipeline::FROM_OUTPUT_PORT());

  vtkInformation *outputInfo=outputInfoVec->GetInformationObject(port);

  // Update the timestep.  
  int tsLength =
    outputInfo->Length(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
  double *steps = outputInfo->Get(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
  int closestStep = 0;

  if(outputInfo->Has(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP()))
    {
    // Get the requested time step. We only supprt requests of a single time
    // step in this reader right now
    double requestedTimeStep =
      outputInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP());

    int cnt=0;
    double minDist=-1;
    for (cnt=0;cnt<tsLength;cnt++)
      {
      double tdist=(steps[cnt]-requestedTimeStep>requestedTimeStep-steps[cnt])?
        steps[cnt]-requestedTimeStep:
        requestedTimeStep-steps[cnt];
      if (minDist<0 || tdist<minDist)
        {
        minDist=tdist;
        closestStep=cnt;
        }
      }
    }

  this->CurrentTimeStep = closestStep;

  if (outputData != NULL) 
    {
    outputData->GetInformation()->Set(vtkDataObject::DATA_TIME_STEP(),
                            steps[this->CurrentTimeStep]);
    }
  return 1;
}
//-----------------------------------------------------------------------------
int vtkSpyPlotReader::AddBlockIdArray(vtkCompositeDataSet *cds)
{
  int blockId=0;
  vtkSmartPointer<vtkCompositeDataIterator> cdIter;
  cdIter.TakeReference(cds->NewIterator());
  for (cdIter->InitTraversal(); !cdIter->IsDoneWithTraversal();
    cdIter->GoToNextItem(), blockId++)
    {
    vtkDataObject *dataObject = cdIter->GetCurrentDataObject();
    if (dataObject!=0)
      {
      vtkDataSet *ds = vtkDataSet::SafeDownCast(dataObject);
      assert("check: ds_exists" && ds!=0);

      // Add the block id cell data array
      vtkCellData *cd=ds->GetCellData();
      vtkDataArray *array=cd->GetArray("blockId");
      if(array!=0)
        {
        cd->RemoveArray("blockId"); // if this is not the first step,
        // make sure we have a clean array
        }
      array = vtkIntArray::New();
      cd->AddArray(array);
      array->Delete();
      array->SetName("blockId");
      array->SetNumberOfComponents(1);
      vtkIdType numCells = ds->GetNumberOfCells();
      array->SetNumberOfTuples(numCells);
      array->FillComponent(0,blockId);
      }
    }

  return 1;
}
//-----------------------------------------------------------------------------
int vtkSpyPlotReader::AddActiveBlockArray(
      vtkCellData *cd,
      vtkIdType nCells,
      unsigned char status)
{
  vtkUnsignedCharArray* activeArray = vtkUnsignedCharArray::New();
  activeArray->SetName("ActiveBlock");
  //vtkIdType numCells=realDims[0]*realDims[1]*realDims[2];
  activeArray->SetNumberOfTuples(nCells);
  activeArray->FillComponent(0,status);
//   for (vtkIdType myIdx=0; myIdx<numCells; myIdx++)
//     {
//     activeArray->SetValue(myIdx, block->IsActive());
//     }
  cd->AddArray(activeArray);
  activeArray->Delete();

  return 1;
}
//-----------------------------------------------------------------------------
int vtkSpyPlotReader::AddAttributes(vtkNonOverlappingAMR *hbds)
{
  // global bounds
  double b[6];
  this->Bounds->GetBounds(b);
  vtkDoubleArray *da = vtkDoubleArray::New();
  da->SetNumberOfComponents(1);
  da->SetNumberOfTuples(6);
  da->SetName("GlobalBounds");
  for (int q=0; q<6; ++q)
    {
    da->SetValue(q, b[q]);
    }
  hbds->GetFieldData()->AddArray( da );
  da->Delete();
  // global box size
  vtkIntArray *ia = vtkIntArray::New();
  ia->SetNumberOfComponents(1);
  ia->SetNumberOfTuples(3);
  ia->SetName("GlobalBoxSize");
  for (int q=0; q<3; ++q)
    {
    ia->SetValue(q, this->BoxSize[q]);
    }
  hbds->GetFieldData()->AddArray( ia );
  ia->Delete();
  // minimum level in use
  ia = vtkIntArray::New();
  ia->SetNumberOfComponents(1);
  ia->SetNumberOfTuples(1);
  ia->SetName("MinLevel");
  ia->SetValue(0, this->MinLevel);
  hbds->GetFieldData()->AddArray( ia );
  ia->Delete();
  // grid spacing on the min level
  da = vtkDoubleArray::New();
  da->SetNumberOfComponents(1);
  da->SetNumberOfTuples(3);
  da->SetName("MinLevelSpacing");
  for (int q=0; q<3; ++q)
    {
    da->SetValue(q, this->MinLevelSpacing[q]);
    }
  hbds->GetFieldData()->AddArray( da );
  da->Delete();

  return 1;
}

//-----------------------------------------------------------------------------
int vtkSpyPlotReader::RequestData(
  vtkInformation *request,
  vtkInformationVector **vtkNotUsed(inputVector),
  vtkInformationVector *outputVector)
{
  vtkDebugMacro( "--------------------------- Request Data --------------------------------" );
  vtkSpyPlotUniReader* uniReader = 0;

  std::vector<vtkRectilinearGrid*> grids;

  vtkInformation *info=outputVector->GetInformationObject(0);
  vtkDataObject *doOutput=info->Get(vtkDataObject::DATA_OBJECT());
  vtkCompositeDataSet *cds=vtkCompositeDataSet::SafeDownCast(doOutput);
  if(!cds)
    {
    vtkErrorMacro("The output is not a CompositeDataSet");
    return 0;
    }
  if (!info->Has(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER()) ||
      !info->Has(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES()))
    {
    vtkErrorMacro("Expected information not found. "
                  "Cannot provide update extent.");
    return 0;
    }

  vtkPolyData *tracersData = NULL;

  if (this->GenerateTracerArray == 1)
    {
    info=outputVector->GetInformationObject(1);
    tracersData = vtkPolyData::SafeDownCast(
                  info->Get(vtkDataObject::DATA_OBJECT()));
    }

  cds->Initialize(); // remove all previous blocks
  //int numFiles = this->Map->Files.size();

  const int myGlobalProcId
    = this->GlobalController->GetLocalProcessId();
  const int nProcsAll
    = this->GlobalController->GetNumberOfProcesses();

  // Set to the requested timestep
  this->UpdateTimeStep(request,outputVector,cds);

  // Tell all of the unireaders that they need to make to check to see
  // if they are current
  this->Map->TellReadersToCheck(this);

  vtkSpyPlotBlock *block;
  vtkSpyPlotBlockIterator *blockIterator;
  if(this->DistributeFiles)
    {
    vtkDebugMacro("Distribute files");
    blockIterator=new vtkSpyPlotFileDistributionBlockIterator;
    }
  else
    {
    vtkDebugMacro("Distribute blocks");
    blockIterator=new vtkSpyPlotBlockDistributionBlockIterator;
    }

  // Block iterator could use either the global proc id and number of
  // processes from the global control or those from the sub controller
  blockIterator->Init(
                nProcsAll,
                myGlobalProcId,
                this,
                this->Map,
                this->CurrentTimeStep);

  int nBlocks = blockIterator->GetNumberOfBlocksToProcess();
  int progressInterval = nBlocks / 10 + 1;
  int rightHasBounds = 0;
  int leftHasBounds = 0;

  vtkNonOverlappingAMR *hbds
        = vtkNonOverlappingAMR::SafeDownCast(cds);

  // TODO The following three calls compute meta data
  // this could be done in a single pass, and the meta data
  // should be vtkInformationKeys defined in vtkNonOverlappingAMR
  // rather than placed in the field data as it is here.

  // Note that in the process of getting the bounds
  // all of the readers will get updated appropriately
  this->SetGlobalBounds(blockIterator, nBlocks,
                        progressInterval, &rightHasBounds,
                        &leftHasBounds);
  // Determine if the box size is constant
  this->SetGlobalBoxSize( blockIterator );
  // Determine the minimum level in use
  // and its grid spacing
  this->SetGlobalMinLevelAndSpacing( blockIterator );
  // export global bounds, minimum level, spacing, and box size
  // in field data arrays for use by downstream filters
  if ( hbds != NULL )
    {
    this->AddAttributes(hbds);
    assert("FieldData should not be NULL!" && hbds->GetFieldData()!=NULL );

    vtkFieldData *fd = hbds->GetFieldData();
    (void) fd; // only used for asserts
    assert("Must have a GlobalBounds array!" && fd->HasArray("GlobalBounds"));
    assert("Must have a GlobalBoxSize array!" && fd->HasArray("GlobalBoxSize"));
    assert("Must have a MinLevel array!" && fd->HasArray("MinLevel"));
    assert("Must have a MinLevelSpacing!" && fd->HasArray("MinLevelSpacing"));
    }

  int needTracers = 1;

  // read in the data
  if (nBlocks!=0)
    {
    // TODO This seems wrong, shouldn't the bounds key be defined in
    // the hierarchical dataset?
    double b[6];
    this->Bounds->GetBounds(b);
    info->Set(vtkStreamingDemandDrivenPipeline::BOUNDS(), b, 6);

    // Gather blocks per level information.
    // Note that this is a quick fix. Taking this extra pass shouldbe avoided. -Leo
    if(this->IsAMR)
      {
      std::vector<int> blocksPerLevel;
      for(blockIterator->Start(); blockIterator->IsActive();  blockIterator->Next())
        {
        block=blockIterator->GetBlock();
        int level = 0;
        int realExtents[6];
        int realDims[3];
        int extents[6];
        double spacing[3];
        double origin[3];
        // int hasBadGhostCells =
        block->GetAMRInformation(*this->Bounds,
                                 &level, spacing,
                                 origin, extents,
                                 realExtents,
                                 realDims);
        for(int i=static_cast<int>(blocksPerLevel.size()); i<=level; i++)
          {
          blocksPerLevel.push_back(0);
          }
        blocksPerLevel[level]++;
        }
      hbds->Initialize(static_cast<int>(blocksPerLevel.size()), &blocksPerLevel[0]);
      }

    // Read the blocks/files that are assigned to this process
    int current_block_number;
    std::vector<int> blockId;
    if(this->IsAMR)
      {
      blockId.resize(hbds->GetNumberOfLevels(),0);
      }
    for(blockIterator->Start(), current_block_number=1;
        blockIterator->IsActive();
        blockIterator->Next(), current_block_number++)
      {
      if ( !(current_block_number % progressInterval) )
        {
        this->UpdateProgress(
          0.6 +
          0.4 * static_cast<double>(current_block_number)/nBlocks);
        }
      block=blockIterator->GetBlock();
      int numFields=blockIterator->GetNumberOfFields();
      uniReader=blockIterator->GetUniReader();

      if (this->GenerateTracerArray == 1 && needTracers)
        {
        tracersData->Initialize ();

        vtkFloatArray *tracers = uniReader->GetTracers ();
        if (tracers != 0)
          {

          vtkPoints* points = vtkPoints::New ();
          points->SetNumberOfPoints (tracers->GetNumberOfTuples ());
          points->SetData (tracers);
          tracersData->SetPoints (points);
          points->Delete ();

          vtkCellArray* verts = vtkCellArray::New ();
          for (vtkIdType p = 0; p < tracers->GetNumberOfTuples (); p ++) 
            {
            verts->InsertNextCell (1, &p);
            }
          tracersData->SetVerts (verts);
          verts->Delete ();
          }
        needTracers = 0;
        }

      int dims[3];
      int blockID = blockIterator->GetBlockID();
      int level = 0;
      int hasBadGhostCells;
      int realExtents[6];
      int realDims[3];
      int extents[6];
      vtkCellData* cd;

      // Get the dimensions of the block
      block->GetDimensions(dims);

      // read block
      if (this->IsAMR)
        {
        hasBadGhostCells
          = this->PrepareAMRData(hbds, block, &level, &blockId[0],
                                 extents, realExtents,
                                 realDims, &cd);
        }
      else
        {
        vtkRectilinearGrid *rg;
        vtkDebugMacro( "Preparing Block: " << blockID << " "
                      << uniReader->GetFileName() );
        hasBadGhostCells = this->PrepareData(
          vtkMultiBlockDataSet::SafeDownCast(cds), block, &rg, extents, realExtents,
          realDims, &cd);

        //currently only for rectilinear grids do we computed derived variables
        grids.push_back(rg);
        }

      vtkDebugMacro("Executing block: " << blockID);

      // Deal with the field data and map out where the true
      // ghost cells are
      if(!hasBadGhostCells)
        {
        this->UpdateFieldData(numFields, dims, level, blockID,
                              uniReader, cd);
        }
      else // we have some bad ghost cells
        {
        this->UpdateBadGhostFieldData(numFields, dims, realDims,
                                      realExtents, level, blockID,
                                      uniReader, cd);
        }
      if(!this->IsAMR)
      {
        this->ComputeDerivedVars(cd, block, uniReader, blockID);
      }

      // Add active block array, for debugging
      if (this->GenerateActiveBlockArray)
        {
        this->AddActiveBlockArray(cd,
                                  realDims[0]*realDims[1]*realDims[2],
                                  block->IsActive());
        }
      // vectorize
      if (this->MergeXYZComponents)
        {
        this->MergeVectors(cd);
        }
      }
    delete blockIterator;
    }

#ifdef PARAVIEW_ENABLE_SPYPLOT_MARKERS
  if (this->GenerateMarkers)
    {
    info=outputVector->GetInformationObject(2);
    doOutput=info->Get(vtkDataObject::DATA_OBJECT());
    vtkMultiBlockDataSet *mbds=vtkMultiBlockDataSet::SafeDownCast(doOutput);

    mbds->SetNumberOfBlocks (0);

    int maxMat = 0;
    vtkSpyPlotReaderMap::MapOfStringToSPCTH::iterator mapIt;
    for ( mapIt = this->Map->Files.begin();
          mapIt != this->Map->Files.end();
          ++ mapIt )
      {
      uniReader = this->Map->GetReader(mapIt, this);
      if (uniReader->GetMarkersOn () == false) 
        {
        // no markers available in this file, don't bother.
        continue;
        }
      int mat = uniReader->GetNumberOfMaterials ();
      if (mat > maxMat) 
        {
        maxMat = mat;      
        }
      }

    int gMaxMat=0;
    this->GlobalController->AllReduce (&maxMat, &gMaxMat, 1, vtkCommunicator::MAX_OP);
    maxMat = gMaxMat;

    this->PrepareBlocks (mbds, maxMat);

    for ( mapIt = this->Map->Files.begin();
          mapIt != this->Map->Files.end();
          ++ mapIt )
      {
      uniReader = this->Map->GetReader(mapIt, this);
      if (uniReader->GetMarkersOn () == false) 
        {
        // no markers available in this file, don't bother.
        continue;
        }
      this->PrepareMarkers (mbds, uniReader);
      }
    }
#endif // PARAVIEW_ENABLE_SPYPLOT_MARKERS

    // At this point, each processor has its own blocks
    // They have to exchange the blocks they have get a unique id for
    // each block over the all dataset.

    // Update the number of levels.
  if(this->GlobalController && this->GlobalController->GetNumberOfProcesses()>1)
    {
    this->SetGlobalLevels(cds);
    }
  // Set the unique block id cell data
  if(this->GenerateBlockIdArray)
    {
    this->AddBlockIdArray(cds);
    }

  return 1;
}

//-----------------------------------------------------------------------------
void vtkSpyPlotReader::AddGhostLevelArray(int numLevels)
{
  (void)numLevels;
  /*
    int numCells = output->GetNumberOfCells();
    int numBlocks = output->GetNumberOfBlocks();
    vtkUnsignedCharArray* array = vtkUnsignedCharArray::New();
    int blockId;
    int dims[3];
    int i, j, k;
    unsigned char* ptr;
    int iLevel, jLevel, kLevel, tmp;

    output->SetNumberOfGhostLevels(numLevels);
    array->SetNumberOfTuples(numCells);
    ptr = (unsigned char*)(array->GetVoidPointer(0));


    for (blockId = 0; blockId < numBlocks; ++blockId)
    {
    output->GetBlockCellDimensions(blockId, dims);
    for (k = 0; k < dims[2]; ++k)
    {
    kLevel = numLevels - k;
    tmp = k - dims[2] + 1 + numLevels;
    if (tmp > kLevel) { kLevel = tmp;}
    if (dims[2] == 1) {kLevel = 0;}
    for (j = 0; j < dims[1]; ++j)
    {
    jLevel = kLevel;
    tmp = numLevels - j;
    if (tmp > jLevel) { jLevel = tmp;}
    tmp = j - dims[1] + 1 + numLevels;
    if (tmp > jLevel) { jLevel = tmp;}
    if (dims[1] == 1) {jLevel = 0;}
    for (i = 0; i < dims[0]; ++i)
    {
    iLevel = jLevel;
    tmp = numLevels - i;
    if (tmp > iLevel) { iLevel = tmp;}
    tmp = i - dims[0] + 1 + numLevels;
    if (tmp > iLevel) { iLevel = tmp;}
    if (dims[0] == 1) {iLevel = 0;}
    if (iLevel <= 0)
    {
    *ptr = 0;
    }
    else
    {
    *ptr = iLevel;
    }
    ++ptr;
    }
    }
    }
    }

    //array->SetName("Test");
    array->SetName("vtkGhostLevels");
    output->GetCellData()->AddArray(array);
    array->Delete();
  */
}

//-----------------------------------------------------------------------------
int vtkSpyPlotReader::GetNumberOfCellArrays()
{
  return this->CellDataArraySelection->GetNumberOfArrays();
}

//-----------------------------------------------------------------------------
const char* vtkSpyPlotReader::GetCellArrayName(int index)
{
  return this->CellDataArraySelection->GetArrayName(index);
}

//-----------------------------------------------------------------------------
int vtkSpyPlotReader::GetCellArrayStatus(const char* name)
{
  return this->CellDataArraySelection->ArrayIsEnabled(name);
}

//-----------------------------------------------------------------------------
void vtkSpyPlotReader::SetCellArrayStatus(const char* name, int status)
{
  vtkDebugMacro("Set cell array \"" << name << "\" status to: " << status);
  if (this->CellDataArraySelection->ArrayIsEnabled(name) == (status? 1 : 0))
    {
    return;
    }

  if(status)
    {
    this->CellDataArraySelection->EnableArray(name);
    }
  else
    {
    this->CellDataArraySelection->DisableArray(name);
    }

  this->Modified();
}

//-----------------------------------------------------------------------------
void vtkSpyPlotReader::MergeVectors(vtkDataSetAttributes* da)
{
  int numArrays = da->GetNumberOfArrays();
  int idx;
  vtkDataArray *a1, *a2, *a3;
  int flag = 1;

  // Loop merging arrays.
  // Since we are modifying the list of arrays that we are traversing,
  // merge one set of arrays at a time.
  while (flag)
    {
    flag = 0;
    for (idx = 0; idx < numArrays-1 && !flag; ++idx)
      {
      a1 = da->GetArray(idx);
      a2 = da->GetArray(idx+1);
      if (idx+2 < numArrays)
        {
        a3 = da->GetArray(idx+2);
        if (this->MergeVectors(da, a1, a2, a3))
          {
          flag = 1;
          continue;
          }
        if (this->MergeVectors(da, a3, a2, a1))
          {
          flag = 1;
          continue;
          }
        }
      if (this->MergeVectors(da, a1, a2))
        {
        flag = 1;
        continue;
        }
      if (this->MergeVectors(da, a2, a1))
        {
        flag = 1;
        continue;
        }
      }
    }
}

//-----------------------------------------------------------------------------
// The templated execute function handles all the data types.
template <class T>
void vtkMergeVectorComponents(vtkIdType length,
                              T *p1, T *p2, T *p3, T *po)
{
  vtkIdType idx;
  if (p3)
    {
    for (idx = 0; idx < length; ++idx)
      {
      *po++ = *p1++;
      *po++ = *p2++;
      *po++ = *p3++;
      }
    }
  else
    {
    for (idx = 0; idx < length; ++idx)
      {
      *po++ = *p1++;
      *po++ = *p2++;
      *po++ = (T)0;
      }
    }
}

//-----------------------------------------------------------------------------
int vtkSpyPlotReader::MergeVectors(vtkDataSetAttributes* da,
                                   vtkDataArray * a1, vtkDataArray * a2, vtkDataArray * a3)
{
  int prefixFlag = 0;

  if (a1 == 0 || a2 == 0 || a3 == 0)
    {
    return 0;
    }
  if(a1->GetNumberOfTuples() != a2->GetNumberOfTuples() ||
     a1->GetNumberOfTuples() != a3->GetNumberOfTuples())
    { // Sanity check.  Should never happen.
    return 0;
    }
  if(a1->GetDataType()!=a2->GetDataType()||a1->GetDataType()!=a3->GetDataType())
    {
    return 0;
    }
  if(a1->GetNumberOfComponents()!=1 || a2->GetNumberOfComponents()!=1 ||
     a3->GetNumberOfComponents()!=1)
    {
    return 0;
    }
  const char *n1, *n2, *n3;
  size_t e1, e2, e3;
  n1 = a1->GetName();
  n2 = a2->GetName();
  n3 = a3->GetName();
  if (n1 == 0 || n2 == 0 || n3 == 0)
    {
    return 0;
    }
  e1 = strlen(n1)-1;
  e2 = strlen(n2)-1;
  e3 = strlen(n3)-1;
  if(e1!=e2 || e1 != e3)
    {
    return 0;
    }
  if (strncmp(n1+1,n2+1,e1)==0 && strncmp(n1+1,n3+1,e1)==0)
    { // Trailing characters are the same. Check for prefix XYZ.
    if ((n1[0]!='X' || n2[0]!='Y' || n3[0]!='Z') &&
        (n1[0]!='x' || n2[0]!='y' || n3[0]!='z'))
      { // Since last characters are the same, there is no way
      // the names can have postfix XYZ.
      return 0;
      }
    prefixFlag = 1;
    }
  else
    { // Check for postfix case.
    if (strncmp(n1,n2,e1)!=0 || strncmp(n1,n3,e1)!=0)
      { // Not pre or postfix.
      return 0;
      }
    if ((n1[e1]!='X' || n2[e2]!='Y' || n3[e3]!='Z') &&
        (n1[e1]!='x' || n2[e2]!='y' || n3[e3]!='z'))
      { // Tails are the same, but postfix not XYZ.
      return 0;
      }
    }
  // Merge the arrays.
  vtkDataArray* newArray = a1->NewInstance();
  newArray->SetNumberOfComponents(3);
  vtkIdType length = a1->GetNumberOfTuples();
  newArray->SetNumberOfTuples(length);
  void *p1 = a1->GetVoidPointer(0);
  void *p2 = a2->GetVoidPointer(0);
  void *p3 = a3->GetVoidPointer(0);
  void *pn = newArray->GetVoidPointer(0);
  switch (a1->GetDataType())
    {
    vtkTemplateMacro(
      vtkMergeVectorComponents( length,
                                (VTK_TT*)p1,
                                (VTK_TT*)p2,
                                (VTK_TT*)p3,
                                (VTK_TT*)pn));
    default:
      vtkErrorMacro(<< "Execute: Unknown ScalarType");
      return 0;
    }
  if (prefixFlag)
    {
    newArray->SetName(n1+1);
    }
  else
    {
    char* name = new char[e1+2];
    strncpy(name,n1,e1);
    name[e1] = '\0';
    newArray->SetName(name);
    delete [] name;
    }
  da->RemoveArray(n1);
  da->RemoveArray(n2);
  da->RemoveArray(n3);
  da->AddArray(newArray);
  newArray->Delete();
  return 1;
}

//-----------------------------------------------------------------------------
int vtkSpyPlotReader::MergeVectors(vtkDataSetAttributes* da,
                                   vtkDataArray * a1, vtkDataArray * a2)
{
  int prefixFlag = 0;

  if (a1 == 0 || a2 == 0)
    {
    return 0;
    }
  if (a1->GetNumberOfTuples() != a2->GetNumberOfTuples())
    { // Sanity check.  Should never happen.
    return 0;
    }
  if (a1->GetDataType() != a2->GetDataType())
    {
    return 0;
    }
  if (a1->GetNumberOfComponents() != 1 || a2->GetNumberOfComponents() != 1)
    {
    return 0;
    }
  const char *n1, *n2;
  size_t e1, e2;
  n1 = a1->GetName();
  n2 = a2->GetName();
  if (n1 == 0 || n2 == 0)
    {
    return 0;
    }
  e1 = strlen(n1)-1;
  e2 = strlen(n2)-1;
  if (e1 != e2 )
    {
    return 0;
    }
  if ( strncmp(n1+1,n2+1,e1) == 0)
    { // Trailing characters are the same. Check for prefix XYZ.
    if ((n1[0]!='X' || n2[0]!='Y') && (n1[0]!='x' || n2[0]!='y'))
      { // Since last characters are the same, there is no way
      // the names can have postfix XYZ.
      return 0;
      }
    prefixFlag = 1;
    }
  else
    { // Check for postfix case.
    if (strncmp(n1,n2,e1) != 0)
      { // Not pre or postfix.
      return 0;
      }
    if ((n1[e1]!='X' || n2[e2]!='Y') && (n1[e1]!='x' || n2[e2]!='y'))
      { // Tails are the same, but postfix not XYZ.
      return 0;
      }
    }
  // Merge the arrays.
  vtkDataArray* newArray = a1->NewInstance();
  // Creae the third componnt and set to 0.
  newArray->SetNumberOfComponents(3);
  vtkIdType length = a1->GetNumberOfTuples();
  newArray->SetNumberOfTuples(length);
  void *p1 = a1->GetVoidPointer(0);
  void *p2 = a2->GetVoidPointer(0);
  void *pn = newArray->GetVoidPointer(0);
  switch (a1->GetDataType())
    {
    vtkTemplateMacro(
      vtkMergeVectorComponents( length,
                                (VTK_TT*)p1,
                                (VTK_TT*)p2,
                                (VTK_TT*)0,
                                (VTK_TT*)pn));
    default:
      vtkErrorMacro(<< "Execute: Unknown ScalarType");
      return 0;
    }
  if (prefixFlag)
    {
    newArray->SetName(n1+1);
    }
  else
    {
    char* name = new char[e1+2];
    strncpy(name,n1,e1);
    name[e1] = '\0';
    newArray->SetName(name);
    delete [] name;
    }
  da->RemoveArray(n1);
  da->RemoveArray(n2);
  da->AddArray(newArray);
  newArray->Delete();
  return 1;
}

//-----------------------------------------------------------------------------
int vtkSpyPlotReader::CanReadFile(const char* fname)
{
  ifstream ifs(fname, ios::binary|ios::in);
  if ( !ifs )
    {
    return 0;
    }
  vtkSpyPlotIStream spis;
  spis.SetStream(&ifs);
  char magic[8];
  if ( !spis.ReadString(magic, 8) )
    {
    vtkDebugMacro( "Cannot read magic" );
    return 0;
    }
  if ( strncmp(magic, "spydata", 7) != 0 &&
       strncmp(magic, "spycase", 7) != 0 )
    {
    return 0;
    }

  return 1;
}

//-----------------------------------------------------------------------------
void vtkSpyPlotReader::SetGenerateMarkers (int gm)
{
  if (this->GenerateMarkers == gm) 
    {
    return;
    }
#ifdef PARAVIEW_ENABLE_SPYPLOT_MARKERS
  vtkSpyPlotReaderMap::MapOfStringToSPCTH::iterator mapIt;
  for ( mapIt = this->Map->Files.begin();
        mapIt != this->Map->Files.end();
        ++ mapIt )
    {
    this->Map->GetReader(mapIt, this)->SetGenerateMarkers(gm);
    }
  this->GenerateMarkers = gm;
#else
  vtkErrorMacro ("GenerateMarkers is currently disabled.  Please issue a support request to enable.");
#endif // PARAVIEW_ENABLE_SPYPLOT_MARKERS
}

//-----------------------------------------------------------------------------
void vtkSpyPlotReader::SetDownConvertVolumeFraction(int vf)
{
  if ( vf == this->DownConvertVolumeFraction )
    {
    return;
    }
  vtkSpyPlotReaderMap::MapOfStringToSPCTH::iterator mapIt;
  for ( mapIt = this->Map->Files.begin();
        mapIt != this->Map->Files.end();
        ++ mapIt )
    {
    this->Map->GetReader(mapIt, this)->SetDownConvertVolumeFraction(vf);
    }
  this->DownConvertVolumeFraction = vf;
  this->Modified();
}

//-----------------------------------------------------------------------------
void vtkSpyPlotReader::SetMergeXYZComponents(int merge)
{
  if ( merge == this->MergeXYZComponents )
    {
    return;
    }
  this->MergeXYZComponents = merge;
  this->Modified();
}
//-----------------------------------------------------------------------------
void vtkSpyPlotReader::PrintBlockList(vtkNonOverlappingAMR *hbds, int
  vtkNotUsed(myProcId))
{
  unsigned int numberOfLevels=hbds->GetNumberOfLevels();
  unsigned int level;
  //  Display the block list for each level
  numberOfLevels=hbds->GetNumberOfLevels();
  for (level = 0; level<numberOfLevels; level++)
    {
   // cout<<myProcId<<" level="<<level<<"/"<<numberOfLevels<<endl;
    int totalNumberOfDataSets=hbds->GetNumberOfDataSets(level);
    int i;
    for (i = 0; i<totalNumberOfDataSets; i++)
      {
     // cout<<myProcId<<" dataset="<<i<<"/"<<totalNumberOfDataSets;
       if(hbds->GetDataSet(level,i)==0)
        {
       // cout<<" Void"<<endl;
        }
      else
        {
       // cout<<" Exists"<<endl;
        }
      }
    }
  /*
    std::vector<vtkRectilinearGrid*>::iterator it;
    for ( it = grids.begin(); it != grids.end(); ++ it )
    {
    (*it)->Print(cout);
    int cc;
    for ( cc = 0; cc < (*it)->GetCellData()->GetNumberOfArrays(); ++ cc )
    {
    (*it)->GetCellData()->GetArray(cc)->Print(cout);
    }
    }
  */
}
//-----------------------------------------------------------------------------
void vtkSpyPlotReader::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os,indent);

  os << "FileName: " << (this->FileName?this->FileName:"(none)") << endl;
  os << "DistributeFiles: ";
  if(this->DistributeFiles)
    {
    os << "true"<<endl;
    }
  else
    {
    os << "false"<<endl;
    }

  os << "DownConvertVolumeFraction: ";
  if(this->DownConvertVolumeFraction)
    {
    os << "true"<<endl;
    }
  else
    {
    os << "false"<<endl;
    }

  os << "MergeXYZComponents: ";
  if(this->MergeXYZComponents)
    {
    os << "true"<<endl;
    }
  else
    {
    os << "false"<<endl;
    }

  os << "GenerateLevelArray: ";
  if(this->GenerateLevelArray)
    {
    os << "true"<<endl;
    }
  else
    {
    os << "false"<<endl;
    }

  os << "GenerateBlockIdArray: ";
  if(this->GenerateBlockIdArray)
    {
    os << "true"<<endl;
    }
  else
    {
    os << "false"<<endl;
    }

  os << "GenerateActiveBlockArray: ";
  if(this->GenerateActiveBlockArray)
    {
    os << "true"<<endl;
    }
  else
    {
    os << "false"<<endl;
    }

  os << "TimeStep: " << this->TimeStep << endl;
  os << "TimeStepRange: " << this->TimeStepRange[0] << " " << this->TimeStepRange[1] << endl;
  if ( this->CellDataArraySelection )
    {
    os << "CellDataArraySelection:" << endl;
    this->CellDataArraySelection->PrintSelf(os, indent.GetNextIndent());
    }
  if ( this->GlobalController)
    {
    os << "GlobalController:" << endl;
    this->GlobalController->PrintSelf(os, indent.GetNextIndent());
    }
}
// Get the minimum level that actuall has data and
// its grid spacing. If this process has no blocks
// then we return spacing and min level close to infinity.
void vtkSpyPlotReader::GetLocalMinLevelAndSpacing(
                              vtkSpyPlotBlockIterator *biter,
                              int *localMinLevel,
                              double spacing[3]) const
{
  biter->Start();
  // If there are no blocks locally we set
  // the spacing and min level close to infinity
  // and return
  if (!biter->IsActive())
    {
    localMinLevel[0]=VTK_INT_MAX;
    spacing[0]=spacing[1]=spacing[2]=VTK_DOUBLE_MAX;
    return;
    }
  // get first block's level
  biter->GetUniReader()->MakeCurrent();
  vtkSpyPlotBlock *thisBlock= biter->GetBlock();
  *localMinLevel = thisBlock->GetLevel();

  vtkSpyPlotBlock *minLevelBlock=thisBlock;

  // compare to all others
  for ( biter->Next(); biter->IsActive(); biter->Next())
    {
    thisBlock = biter->GetBlock();
    int thisMinLevel=thisBlock->GetLevel();
    if (thisMinLevel < *localMinLevel)
      {
      minLevelBlock=thisBlock;
      *localMinLevel=thisMinLevel;
      }
    }
  // now that we have the block with minimum
  // level, we can get the spacing from him.
  minLevelBlock->GetSpacing(spacing);

  return;
}
//
void vtkSpyPlotReader::SetGlobalMinLevelAndSpacing(
                            vtkSpyPlotBlockIterator *biter)
{
  // get the local min level and its spacing
  // if there are no blocks locally these will be
  // close to infinity.
  int localMinLevel;
  double localMinLevelSpacing[3];
  this->GetLocalMinLevelAndSpacing(biter,&localMinLevel,localMinLevelSpacing);

  // If we are not running in parallel then the local
  // min level and spacing is the global min level and spacing
  if (!this->GlobalController)
    {
    this->MinLevel=localMinLevel;
    for (int q=0; q<3; ++q)
      {
      this->MinLevelSpacing[q]=localMinLevelSpacing[q];
      }
    return;
    }

  // package the level and its spacing because
  // we want to reduce min level but also need
  // its associated spacing to go along for the ride.
  double sendBuf[4]={(double)localMinLevel,
                     localMinLevelSpacing[0],
                     localMinLevelSpacing[1],
                     localMinLevelSpacing[2]};
  // collect the information on proc 0
  int numProcs=this->GlobalController->GetNumberOfProcesses();
  int myRank=this->GlobalController->GetLocalProcessId();

  // proc 0 gets a bigger buffer to collect all the data
  double *recvBuf = myRank==0 ? new double[4*numProcs] : 0;
  this->GlobalController->Gather( sendBuf, recvBuf, 4, 0);

  // reduce min level while preserving its
  // associated spacing
  if (myRank==0)
    {
    // sendBuf will hold the result of the reduction
    // it is already initialized with min level and
    // spacing from proc 0, we start looking in data
    // gathered from proc 1...n
    for (int i=0,j=4; i<numProcs-1; ++i, j+=4)
      {
      if (sendBuf[0]>recvBuf[j])
        {
        // copy
        for (int q=0; q<4; ++q)
          {
          sendBuf[q]=recvBuf[j+q];
          }
        }
      }
    delete [] recvBuf;
    }
  // proc 0 knows the min level and spacing
  // share this with everyone else
  this->GlobalController->Broadcast(sendBuf,4,0);

  // now everyone has the min level and associated grid spacing
  this->MinLevel=(int)sendBuf[0];
  for (int q=0; q<3; ++q)
    {
    this->MinLevelSpacing[q]=sendBuf[1+q];
    }
}


// Determine the local box size if the box size is constant.
// If the box size varies over the local tree then size will be
// -1,-1,-1.
// The case where there are no blocks is treated as if the block size
// is constant and the size will be
// inf, inf, inf
//
// returns true if the box size is constant and false if the box size varies
bool vtkSpyPlotReader::GetLocalBoxSize(vtkSpyPlotBlockIterator *biter,
                                       int localBoxSize[3]) const
{
  vtkSpyPlotBlock *block;
  biter->Start();
  // see if there are any blocks on this process
  // if not then, we'll call block size constant
  // and set the size close to infinite
  if (!biter->IsActive())
    {
    localBoxSize[0]=localBoxSize[1]=localBoxSize[2]=VTK_INT_MAX;
    return true;
    }
  // get first box's size
  biter->GetUniReader()->MakeCurrent();
  block = biter->GetBlock();
  int thisBoxSize[3];
  block->GetDimensions(localBoxSize);
  // compare to all others
  for ( biter->Next(); biter->IsActive(); biter->Next())
    {
    block = biter->GetBlock();
    biter->GetUniReader()->MakeCurrent();
    block->GetDimensions(thisBoxSize);
    for (int q=0;q<3;++q)
      {
      // if the size changes stop
      if (thisBoxSize[q]!=localBoxSize[q])
        {
        localBoxSize[0] = -1;
        localBoxSize[1] = -1;
        localBoxSize[2] = -1;
        return false;
        }
      }
    }
  // box size is a constant locally
  return true;
}

// Determines if the box size is a constant over the entire data set
// if so sets this->GlobalBoxSize to that size, otehrwise sets
// it to -1,-1,-1
void vtkSpyPlotReader::SetGlobalBoxSize(vtkSpyPlotBlockIterator *biter)
{
  // Get the local box size
  int localBoxSize[3] = {0,0,0};
  bool isConstantLocally
      = this->GetLocalBoxSize(biter, localBoxSize);

  // If we are not running in parallel then the local
  // size is the global size
  if (!this->GlobalController)
    {
    if (isConstantLocally)
      {
      for (int q=0; q<3; ++q)
        {
        this->BoxSize[q]=localBoxSize[q];
        }
      }
    else
      {
      for (int q=0; q<3; ++q)
        {
        this->BoxSize[q]=-1;
        }
      }
    return;
    }

  // To decide weather or not box size is constant...
  // 1) get the smallest box size across procs
  int globalBoxSize[3] = {-1,-1,-1};
  this->GlobalController->AllReduce(localBoxSize, globalBoxSize, 3, vtkCommunicator::MIN_OP);

  // 2) check if box size we have is the same as the
  // box size that the other procs have.
  bool isConstantGlobally=true;
  for (int q=0; q<3; ++q)
    {
    // if we have no blocks then we have to
    // get the true size from other procs, and pretend
    // that size is constant locally
    if (localBoxSize[q]==VTK_INT_MAX )
      {
      localBoxSize[q]=globalBoxSize[q];
      }
    // if we have blocks the n we just check to see if
    // what we have is the same as what the others
    // have.
    else if (localBoxSize[q]!=globalBoxSize[q])
      {
      isConstantGlobally=false;
      }
    }
  // 3) send a flag indicating change/no change occured
  int lFlag=!isConstantLocally||!isConstantGlobally ? -1 : 1;
  int gFlag=0;
  this->GlobalController->AllReduce(&lFlag,&gFlag,1, vtkCommunicator::MIN_OP);

  // set the global box size acordingly
  switch( gFlag )
    {
    // box size varies
    case -1:
      for (int q=0; q<3; ++q)
        {
        this->BoxSize[q]=-1;
        }
      break;
    // box size is constant
    case 1:
      for (int q=0; q<3; ++q)
        {
        this->BoxSize[q]=localBoxSize[q];
        }
        break;
    default:
      vtkErrorMacro("Invalid flag value verifying that box size is constant.");
    }
 return;
}

// This functions return 1 if the bounds have been set
void vtkSpyPlotReader::GetLocalBounds(vtkSpyPlotBlockIterator *biter,
                                      int nBlocks, int progressInterval)
{
  int i;
  double bounds[6];
  double progressFactor = 0.4 / static_cast<double>(nBlocks);
  vtkSpyPlotBlock *block;

  // every one has bounds, the procs with no blocks
  // have canonical empty bounds.
//   double *pt;
//   pt=const_cast<double *>(this->Bounds->GetMinPoint());
//   pt[0]=pt[1]=pt[2]= VTK_DOUBLE_MAX;
//   pt=const_cast<double *>(this->Bounds->GetMaxPoint());
//   pt[0]=pt[1]=pt[2]=-VTK_DOUBLE_MAX;

  biter->Start();
  for (i = 0; biter->IsActive(); i++, biter->Next())
    {
    // See if we need to update progress
    if (i && !(i % progressInterval))
      {
      this->UpdateProgress(
        static_cast<double>(1.2 + i) * progressFactor);
      }
    // Make sure that the block is up to date
    biter->GetUniReader()->MakeCurrent();
    block = biter->GetBlock();
    block->GetRealBounds(bounds);
    this->Bounds->AddBounds(bounds);
    }
}

// Returns 1 if the bounds are valid else 0
void vtkSpyPlotReader::SetGlobalBounds(vtkSpyPlotBlockIterator *biter,
                                       int total_num_of_blocks,
                                       int progressInterval,
                                       int *rightHasBounds,
                                       int *leftHasBounds)
{
  // Get the local bounds of this reader
  this->GetLocalBounds(biter, total_num_of_blocks,
                       progressInterval);

  // If we are not running in parallel then the local
  // bounds are the global bounds
  if (!this->GlobalController)
    {
    return;
    }
  vtkCommunicator *comm = this->GlobalController->GetCommunicator();
  if (!comm)
    {
    return;
    }

  int processNumber=this->GlobalController->GetLocalProcessId();
  int numProcessors=this->GlobalController->GetNumberOfProcesses();

  if (!comm->ComputeGlobalBounds(processNumber, numProcessors,
                                 this->Bounds, rightHasBounds,
                                 leftHasBounds,
                                 VTK_MSG_SPY_READER_HAS_BOUNDS,
                                 VTK_MSG_SPY_READER_LOCAL_BOUNDS,
                                 VTK_MSG_SPY_READER_GLOBAL_BOUNDS))
    {
    vtkErrorMacro("Problem occurred getting the global bounds");
    }

 return;
}

int vtkSpyPlotReader::PrepareAMRData(vtkNonOverlappingAMR *hb,
                                     vtkSpyPlotBlock *block,
                                     int *level,
                                     int *blockId,
                                     int extents[6],
                                     int realExtents[6],
                                     int realDims[3],
                                     vtkCellData **cd)
{
  double spacing[3];
  double origin[3];
  int needsFixing;

  needsFixing = block->GetAMRInformation(*(this->Bounds),
                                         level, spacing,
                                         origin, extents,
                                         realExtents,
                                         realDims);

//   double bds[6];
//   this->Bounds->GetBounds(bds);
//   cerr << "{\n";
//   cerr << "level:       [" << *level << "]\n";
//   cerr << "Origin:      [" << origin[0] << "," << origin[1] << "," << origin[2] << "]\n";
//   cerr << "Spacing:     [" << spacing[0] << "," << spacing[1] << "," << spacing[2] << "]\n";
//   cerr << "extents:     [" << extents[0] << "," << extents[1] << "," << extents[2] << "|"
//                            << extents[3] << "," << extents[4] << "," << extents[5] << "]\n";
//   cerr << "realExtents: [" << realExtents[0] << "," << realExtents[1] << "," << realExtents[2] << ","
//                            << realExtents[3] << "," << realExtents[4] << "," << realExtents[5] << "]\n";
//   cerr << "realDims:    [" << realDims[0] << "," << realDims[1] << "," << realDims[2] << "]\n";
//   cerr << "bounds:      [" << bds[0] << "," << bds[1] << "," << bds[2] << ","
//                            << bds[3] << "," << bds[4] << "," << bds[5] << "]\n";
//   cerr << "}\n";
//

//  vtkAMRBox box(realExtents);
//  hb->SetDataSet(*level, hb->GetNumberOfDataSets(*level), box, ug);
  vtkUniformGrid* ug = vtkUniformGrid::New();

  ug->SetSpacing(spacing);
  ug->SetExtent(extents);
  ug->SetOrigin(origin);
  *cd = ug->GetCellData();

  hb->SetDataSet( *level, blockId[*level], ug );
  blockId[*level]++;
  ug->Delete();
  return needsFixing;
}

int vtkSpyPlotReader::PrepareData(vtkMultiBlockDataSet *hb,
                                  vtkSpyPlotBlock *block,
                                  vtkRectilinearGrid **rg,
                                  int extents[6],
                                  int realExtents[6],
                                  int realDims[3],
                                  vtkCellData **cd)
{

  int needsFixing;
  vtkDataArray *coordinates[3];
  needsFixing = block->FixInformation(*(this->Bounds),
                                      extents, realExtents,
                                      realDims,coordinates);
  double bounds[6];
  this->Bounds->GetBounds(bounds);
  vtkDebugMacro( << __LINE__ << " Real dims:    "
                 << coutVector3(realDims) );
  vtkDebugMacro( << __LINE__ << " Real Extents: "
                 << coutVector6(realExtents) );
  vtkDebugMacro( << __LINE__ << " Extents:      "
                 << coutVector6(extents) );
  vtkDebugMacro( << __LINE__ << " Global Bounds:"
                 << coutVector6(bounds) );
  vtkDebugMacro( << " Rectilinear grid pointer: " << rg );
  *rg=vtkRectilinearGrid::New();
  (*rg)->SetExtent(extents);
  hb->SetBlock(hb->GetNumberOfBlocks(), *rg);
  if (coordinates[0])
    {
    (*rg)->SetXCoordinates(coordinates[0]);
    vtkDebugMacro( "NT: " << coordinates[0]->GetNumberOfTuples() );
    }
  if (coordinates[1])
    {
    (*rg)->SetYCoordinates(coordinates[1]);
    vtkDebugMacro( "NT: " << coordinates[1]->GetNumberOfTuples() );
    }
  if (coordinates[2])
    {
    (*rg)->SetZCoordinates(coordinates[2]);
    vtkDebugMacro( "NT: " << coordinates[2]->GetNumberOfTuples() );
    }

  vtkDebugMacro( "*******************" );
  vtkDebugMacro( "Coordinates: " );
  /*
    int cor;
    for ( cor = 0; cor < 3; cor ++ )
    {
    if (coordinates[cor])
    {
    coordinates[cor]->Print(cout);
    }
    else
    {
    vtkDebugMacro("No " << cor << "th coordinates");
    }
    vtkDebugMacro( "*******************" );
  */
  *cd = (*rg)->GetCellData();
  (*rg)->Delete();
  return needsFixing;
}

void vtkSpyPlotReader::UpdateFieldData(int numFields, int dims[3],
                                       int level, int blockID,
                                       vtkSpyPlotUniReader *uniReader,
                                       vtkCellData *cd)
{
  int field;
  int fixed= 0;
  int totalSize = dims[0]*dims[1]*dims[2];
  const char *fname;
  for ( field = 0; field < numFields; field ++ )
    {
    fname = uniReader->GetCellFieldName(field);
    if (this->CellDataArraySelection->ArrayIsEnabled(fname))
      {
      vtkDataArray *array=cd->GetArray(fname);
      if(array!=0)
        {
        cd->RemoveArray(fname); // if this is not the first step,
        // make sure we have a clean array
        }
      array = uniReader->GetCellFieldData(blockID, field, &fixed);

      //vtkDebugMacro( << __LINE__ << " Read data block: "
      // << blockID << " " << field << "  [" << array->GetName() << "]" );
      cd->AddArray(array);
      }
    }

  // Add a level array, for debugging
  if(this->GenerateLevelArray)
    {
    createSpyPlotLevelArray(cd, totalSize, level);
    }

  // Mark the bounding cells as ghost cells
  vtkUnsignedCharArray *ghostArray=vtkUnsignedCharArray::New();
  ghostArray->SetNumberOfTuples(totalSize);
  ghostArray->SetName(vtkDataSetAttributes::GhostArrayName());
  cd->AddArray(ghostArray);
  ghostArray->Delete();
  int planeSize = dims[0]*dims[1];
  int j, k;
  unsigned char *ptr =
    static_cast<unsigned char*>(ghostArray->GetVoidPointer(0));
  for (k = 0; k < dims[2]; k++)
    {
    // Is the entire ij plane a set of ghosts
    if ((dims[2] != 1) && ((!k) || (k == dims[2]-1)))
      {
      memset(ptr, vtkDataSetAttributes::DUPLICATECELL, planeSize);
      ptr += planeSize;
      continue;
      }
    for (j = 0; j < dims[1]; j++)
      {
      // Either the row is all non-ghosts except for the ends or
      // its all ghosts
      // Is the entire row a set of ghosts
      if ((dims[1] != 1) && ((!j) || (j == dims[1] - 1)))
        {
        memset(ptr, vtkDataSetAttributes::DUPLICATECELL, dims[0]);
        ptr+= dims[0];
        continue;
        }
      memset(ptr, 0, dims[0]);
      if (dims[0] > 1)
        {
        ptr[0] = 1;
        ptr[dims[0]-1] = 1;
        }
      ptr+= dims[0];
      }
    }
}

void vtkSpyPlotReader::UpdateBadGhostFieldData(int numFields, int dims[3],
                                               int realDims[3],
                                               int realExtents[6],
                                               int level, int blockID,
                                               vtkSpyPlotUniReader *uniReader,
                                               vtkCellData *cd)
{
  int realPtDims[3];
  int ptDims[3];
  int totalSize = realDims[0]*realDims[1]*realDims[2];
  const char *fname;
  int cc, fixed = 0;;
  int field;
  for (cc = 0; cc < 3; cc++)
    {
    realPtDims[cc]=realDims[cc]+1;
    ptDims[cc]=dims[cc]+1;
    }

  for ( field = 0; field < numFields; field ++ )
    {
    fname = uniReader->GetCellFieldName(field);
    if (this->CellDataArraySelection->ArrayIsEnabled(fname))
      {
      vtkDataArray *array=cd->GetArray(fname);
      if(array!=0)
        {
        cd->RemoveArray(fname); // if this is not the first step,
        // make sure we have a clean array
        }

      array = uniReader->GetCellFieldData(blockID, field, &fixed);
      if (array == 0)
        {
        vtkErrorMacro ("Unable to read array " << fname);
        continue;
        }
      //vtkDebugMacro( << __LINE__ << " Read data block: " << blockID
      // << " " << field << "  [" << array->GetName() << "]" );
      cd->AddArray(array);

      if ( !fixed )
        {
        vtkDebugMacro( " Fix bad ghost cells for the array: " << blockID
                       << " / " << field
                       << " (" << uniReader->GetFileName() << ")" );
        switch(array->GetDataType())
          {
          vtkTemplateMacro(
            ::vtkSpyPlotRemoveBadGhostCells(static_cast<VTK_TT*>(0), array,
                                            realExtents, realDims,
                                            ptDims, realPtDims));
          }
        uniReader->MarkCellFieldDataFixed(blockID, field);

        }
      else
        {
        vtkDebugMacro( " Bad ghost cells already fixed for the array" );
        }
      }
    }

  // Add a level array, for debugging
  if(this->GenerateLevelArray)
    {
    createSpyPlotLevelArray(cd, totalSize, level);
    }

  // Mark the remains ghost cell as real ghost cells of level 1.
  vtkUnsignedCharArray *ghostArray=vtkUnsignedCharArray::New();
  ghostArray->SetNumberOfTuples(totalSize);
  ghostArray->SetName(vtkDataSetAttributes::GhostArrayName());
  cd->AddArray(ghostArray);
  ghostArray->Delete();
  unsigned char *ptr
    =static_cast<unsigned char*>(ghostArray->GetVoidPointer(0));
  int k,j;
  int planeSize = realDims[0] * realDims[1];
  int checkILower = (realExtents[0]==0);
  int checkIUpper = (realExtents[1]==dims[0]);
  int checkJLower = (realExtents[2]==0);
  int checkJUpper = (realExtents[3]==dims[1]);
  int checkKLower = (realExtents[4]==0);
  int checkKUpper = (realExtents[5]==dims[2]);
  for (k = 0; k < realDims[2]; k++)
    {
    // Is the entire ij plane a set of ghosts
    if ((realDims[2] != 1) && ((checkKLower && (!k)) ||
                               (checkKUpper && (k == realDims[2]-1))))
      {
      memset(ptr, vtkDataSetAttributes::DUPLICATECELL, planeSize);
      ptr += planeSize;
      continue;
      }
    for (j = 0; j < realDims[1]; j++)
      {
      // Either the row is all non-ghosts except for the ends or
      // its all ghosts
      // Is the entire row a set of ghosts
      if ((realDims[1] != 1) && ((checkJLower && (!j)) ||
                                 (checkJUpper && (j == realDims[1] - 1))))
        {
        memset(ptr, vtkDataSetAttributes::DUPLICATECELL, realDims[0]);
        ptr+= realDims[0];
        continue;
        }
      memset(ptr, 0, realDims[0]);
      if (dims[0] > 1)
        {
        if (checkILower)
          {
          ptr[0] = 1;
          }
        if (checkIUpper)
          {
          ptr[realDims[0]-1] = 1;
          }
        }
      ptr+= realDims[0];
      }
    }
}

void vtkSpyPlotReader::PrepareBlocks (vtkMultiBlockDataSet* mbds, int numMat) 
{
  mbds->SetNumberOfBlocks (numMat);
  for (int m = 0; m < numMat; m ++) 
    {
    vtkPolyData* poly = vtkPolyData::New ();
    poly->Allocate ();
    vtkPoints* points = vtkPoints::New ();
    points->SetNumberOfPoints (0);
    poly->SetPoints (points);
    points->Delete ();

    vtkIntArray* location = vtkIntArray::New ();
    location->SetName ("Location");
    location->SetNumberOfComponents (3);
    location->SetNumberOfTuples (0);
    poly->GetPointData ()->AddArray (location);
    location->Delete ();

    vtkIntArray* blockId = vtkIntArray::New ();
    blockId->SetName ("BlockId");
    blockId->SetNumberOfComponents (1);
    blockId->SetNumberOfTuples (0);
    poly->GetPointData ()->AddArray (blockId);
    blockId->Delete ();

    mbds->SetBlock (m, poly);
    poly->Delete ();
    }
}

int vtkSpyPlotReader::PrepareMarkers (vtkMultiBlockDataSet* mbds, 
                vtkSpyPlotUniReader* reader)
{
  for (int m = 0; m < reader->GetNumberOfMaterials (); m ++)
    {
    if (reader->Markers[m].NumMarks <= 0)
      {
      continue;
      }
    vtkPolyData* poly = vtkPolyData::SafeDownCast (mbds->GetBlock (m));
    vtkPointData* pd = poly->GetPointData ();

    vtkPoints* oldPoints = poly->GetPoints ();
    vtkIntArray* oldLocation = vtkIntArray::SafeDownCast (pd->GetArray ("Location"));
    vtkIntArray* oldBlockId = vtkIntArray::SafeDownCast (pd->GetArray ("BlockId"));

    int offset = oldPoints->GetNumberOfPoints ();
    int length = offset + reader->Markers[m].NumRealMarks;

    vtkPoints* points = vtkPoints::New ();
    points->SetNumberOfPoints (length);

    vtkIntArray* location = vtkIntArray::New ();
    location->SetName ("Location");
    location->SetNumberOfComponents (3);
    location->SetNumberOfTuples (length);

    vtkIntArray* blockId = vtkIntArray::New ();
    blockId->SetName ("BlockId");
    blockId->SetNumberOfComponents (1);
    blockId->SetNumberOfTuples (length);

    for (int o = 0; o < offset; o ++) 
      {
      points->SetPoint (o, oldPoints->GetPoint (o));
      location->SetTuple (o, oldLocation->GetTuple (o));
      blockId->SetTuple (o, oldBlockId->GetTuple (o));
      }

    poly->SetPoints (points);
    points->Delete ();

    pd->AddArray (location);
    location->Delete ();

    pd->AddArray (blockId);
    blockId->Delete ();

    vtkFloatArray** vars = new vtkFloatArray*[reader->Markers[m].NumVars];
    for (int v = 0; v < reader->Markers[m].NumVars; v ++) 
      {

      vtkDataArray *array = pd->GetArray (reader->Markers[m].Variables[v].Label);
      if (array == 0)
        {
        array = vtkFloatArray::New ();
        array->SetName (reader->Markers[m].Variables[v].Label);
        array->SetNumberOfComponents (1);
        array->SetNumberOfTuples (0);
        poly->GetPointData ()->AddArray (array);
        array->Delete ();
        }

      vars[v] = vtkFloatArray::SafeDownCast (array);
      vars[v]->SetNumberOfTuples (length);
      }

    for (int mark = 0; mark < reader->Markers[m].NumRealMarks; mark ++) 
      {
      float x[3];
      int b[3];

      x[0] = reader->MarkersDumps[m].XLoc->GetValue (mark);
      b[0] = reader->MarkersDumps[m].ILoc->GetValue (mark);
      if (reader->GetNumberOfDimensions () > 1)
        {
        x[1] = reader->MarkersDumps[m].YLoc->GetValue (mark);
        b[1] = reader->MarkersDumps[m].JLoc->GetValue (mark);
        }
      else
        {
        x[1] = 0;
        b[1] = 0;
        }
      if (reader->GetNumberOfDimensions () > 2)
        {
        x[2] = reader->MarkersDumps[m].ZLoc->GetValue (mark);
        b[2] = reader->MarkersDumps[m].KLoc->GetValue (mark);
        }
      else
        {
        x[2] = 0;
        b[2] = 0;
        }
      points->SetPoint (mark + offset, x);
      vtkIdType id = mark + offset;
      poly->InsertNextCell (VTK_VERTEX, 1, &id);

      location->SetTupleValue (mark + offset, b);

      blockId->SetValue (mark + offset, reader->MarkersDumps[m].Block->GetValue (mark));

      for (int v = 0; v < reader->Markers[m].NumVars; v ++)
        {
        vars[v]->SetValue (mark + offset, reader->MarkersDumps[m].Variables[v]->GetValue(mark));
        }
      }
    }
  return 1;
}

/*
int vtkSpyPlotReader::DistributeGlobalStructure(vtkCompositeDataSet *cds)
{
  // two cases hierarchical box data and multi block of rectilinear
  // grids

  vtkMultiBlockDataSet* mbDS =
    vtkMultiBlockDataSet::SafeDownCast(composite);

  assert("check: ds must be hierarchical or multiblock" && (hbDS || mbDS));

  if (this->AMR)
    {
    vtkHierarchicalBoxDataSet* hbdS
    = dynamic_cast<vtkHierarchicalBoxDataSet *>(cds);

    unsigned int nLevelsLocal=hbds->GetNumberOfLevels();
    unsigned int nLevelsGlobal=0;
    this->GlobalController->AllReduce(
                &nLevelsLocal,
                &nLevelsGlobal,
                1,vtkCommunicator::MAX_OP);

    }

  return 1;
}*/

//-----------------------------------------------------------------------------
// synch data set structure
void vtkSpyPlotReader::SetGlobalLevels(vtkCompositeDataSet *composite)
{
  // Force the inclusion of empty datatsets.
  bool rightHasBounds=true; //TODO
  bool leftHasBounds=true;
  int processNumber=this->GlobalController->GetLocalProcessId();
  int numProcessors=this->GlobalController->GetNumberOfProcesses();

  int parent = 0;
  int left= vtkCommunicator::GetLeftChildProcessor(processNumber);
  int right=left+1;
  if(processNumber>0) // not root (nothing to do if root)
    {
    parent= vtkCommunicator::GetParentProcessor(processNumber);
    }

  vtkNonOverlappingAMR* hbDS =
    vtkNonOverlappingAMR::SafeDownCast(composite);
  vtkMultiBlockDataSet* mbDS =
    vtkMultiBlockDataSet::SafeDownCast(composite);

  assert("check: ds must be hierarchical or multiblock" && (hbDS || mbDS));
  unsigned int numberOfLevels=1;

  // If this is an AMR SpyPlot we need to first determine
  // the global number of levels. Otherwise, the number of levels is 1.
  if (this->IsAMR)
    {
    // hbDS is non-null.
    assert("check: ds is vtkNonOverlappingAMR" && hbDS !=NULL);
    numberOfLevels = hbDS->GetNumberOfLevels();
    unsigned long ulintMsgValue;
    // Update it from the children
    if(left<numProcessors)
      {
      if(leftHasBounds)
        {
        // Grab the number of levels from left child
        this->GlobalController->Receive(&ulintMsgValue, 1, left,
                                  VTK_MSG_SPY_READER_LOCAL_NUMBER_OF_LEVELS);
        if(numberOfLevels<ulintMsgValue)
          {
          numberOfLevels=ulintMsgValue;
          }
        }
      if(right<numProcessors)
        {
        if(rightHasBounds)
          {
          // Grab the number of levels from right child
          this->GlobalController->Receive(&ulintMsgValue, 1, right,
                                    VTK_MSG_SPY_READER_LOCAL_NUMBER_OF_LEVELS);
          if(numberOfLevels<ulintMsgValue)
            {
            numberOfLevels=ulintMsgValue;
            }
          }
        }
      }

    ulintMsgValue=numberOfLevels;

    // Send local to parent, Receive global from the parent.
    if(processNumber>0) // not root (nothing to do if root)
      {
//      parent=this->GetParentProcessor(processNumber);
      this->GlobalController->Send(&ulintMsgValue, 1, parent,
                             VTK_MSG_SPY_READER_LOCAL_NUMBER_OF_LEVELS);
      this->GlobalController->Receive(&ulintMsgValue, 1, parent,
                                VTK_MSG_SPY_READER_GLOBAL_NUMBER_OF_LEVELS);
      numberOfLevels=ulintMsgValue;
      }

    // Send it to children.
    if(left<numProcessors)
      {
      if(leftHasBounds)
        {
        this->GlobalController->Send(&ulintMsgValue, 1, left,
                               VTK_MSG_SPY_READER_GLOBAL_NUMBER_OF_LEVELS);
        }
      if(right<numProcessors)
        {
        if(rightHasBounds)
          {
          this->GlobalController->Send(&ulintMsgValue, 1, right,
                                 VTK_MSG_SPY_READER_GLOBAL_NUMBER_OF_LEVELS);
          }
        }
      }
    }
  // At this point, the global number of levels is set in each processor.
  // Update each level
  // i.e. for each level synchronize the number of datasets (or pieces).
  std::vector<int> blocksPerLevel(numberOfLevels,0); //collect this
  std::vector<int> globalIndices(numberOfLevels);
  for (unsigned int level = 0; level<numberOfLevels; level++)
    {
    int intMsgValue;
    int numberOfDataSets(0);
    if(hbDS)
      {
      numberOfDataSets = level>=hbDS->GetNumberOfLevels()? 0 : hbDS->GetNumberOfDataSets(level);
      }
    else
      {
      numberOfDataSets = mbDS->GetNumberOfBlocks();
      }
    int totalNumberOfDataSets=numberOfDataSets;
    int leftNumberOfDataSets=0;
    int rightNumberOfDataSets=0;
    // Get number of dataset of each child
    if(left<numProcessors)
      {
      if(leftHasBounds)
        {
        // Grab info the number of datasets from left child
        this->GlobalController->Receive(&intMsgValue, 1, left,
                                  VTK_MSG_SPY_READER_LOCAL_NUMBER_OF_DATASETS);
        leftNumberOfDataSets=intMsgValue;
        }
      if(right<numProcessors)
        {
        if(rightHasBounds)
          {
          // Grab info the number of datasets from right child
          this->GlobalController->Receive(&intMsgValue, 1, right,
                                    VTK_MSG_SPY_READER_LOCAL_NUMBER_OF_DATASETS);
          rightNumberOfDataSets=intMsgValue;
          }
        }
      }

    int globalIndex;
    if(processNumber==0) // root
      {
      totalNumberOfDataSets=numberOfDataSets+leftNumberOfDataSets
        +rightNumberOfDataSets;
      globalIndex=0;
      }
    else
      {
      // Send local to parent, Receive global from the parent.
      intMsgValue=numberOfDataSets+leftNumberOfDataSets+rightNumberOfDataSets;
      this->GlobalController->Send(&intMsgValue, 1, parent,
                             VTK_MSG_SPY_READER_LOCAL_NUMBER_OF_DATASETS);
      this->GlobalController->Receive(&intMsgValue, 1, parent,
                                VTK_MSG_SPY_READER_GLOBAL_NUMBER_OF_DATASETS);
      totalNumberOfDataSets=intMsgValue;
      this->GlobalController->Receive(&intMsgValue, 1, parent,
                                VTK_MSG_SPY_READER_GLOBAL_DATASETS_INDEX);
      globalIndex=intMsgValue;
      }

    globalIndices[level] = globalIndex;

    // Send it to children.
    if(left<numProcessors)
      {
      if(leftHasBounds)
        {
        intMsgValue=totalNumberOfDataSets;
        this->GlobalController->Send(&intMsgValue, 1, left,
                               VTK_MSG_SPY_READER_GLOBAL_NUMBER_OF_DATASETS);
        intMsgValue=globalIndex+numberOfDataSets;
        this->GlobalController->Send(&intMsgValue, 1, left,
                               VTK_MSG_SPY_READER_GLOBAL_DATASETS_INDEX);
        }
      if(right<numProcessors)
        {
        if(rightHasBounds)
          {
          intMsgValue=totalNumberOfDataSets;
          this->GlobalController->Send(&intMsgValue, 1, right,
                                 VTK_MSG_SPY_READER_GLOBAL_NUMBER_OF_DATASETS);
          intMsgValue=globalIndex+numberOfDataSets+leftNumberOfDataSets;
          this->GlobalController->Send(&intMsgValue, 1, right,
                                 VTK_MSG_SPY_READER_GLOBAL_DATASETS_INDEX);
          }
        }
      }
    blocksPerLevel[level] = totalNumberOfDataSets;

    // Update the level.
    if(totalNumberOfDataSets>numberOfDataSets)
      {
      // save the current datasets
      if (hbDS)
        {
        }
      else // if (mbDS)
        {
        if (globalIndex == 0)
          {
          mbDS->SetNumberOfBlocks(totalNumberOfDataSets);
          }
        else
          {
          int kk;
          std::vector<vtkSmartPointer<vtkDataSet> > datasets;
          for (kk=0; kk < numberOfDataSets; kk++)
            {
            datasets.push_back(vtkDataSet::SafeDownCast(mbDS->GetBlock(kk)));
            }
          mbDS->SetNumberOfBlocks(0); // removes all current blocks.
          mbDS->SetNumberOfBlocks(totalNumberOfDataSets);
          // put the datasets back starting at globalIndex.
          // All other indices are in their initialized state.
          for (kk=0; kk < numberOfDataSets; kk++)
            {
            mbDS->SetBlock(kk+globalIndex, datasets[kk]);
            }
          }
        }
      }
    }

  if(hbDS)
    {
    //save all the existing data sets
    std::vector<std::vector<vtkSmartPointer<vtkUniformGrid> > > datasets;
    for(unsigned int level = 0; level<numberOfLevels; level++)
      {
      std::vector<vtkSmartPointer<vtkUniformGrid> > datasetsAtLevel;
      for (unsigned int kk=0; kk < hbDS->GetNumberOfDataSets(level); kk++)
        {
        vtkUniformGrid* ug = hbDS->GetDataSet(level, kk);
        datasetsAtLevel.push_back(ug);
        }
      datasets.push_back(datasetsAtLevel);
      }

    hbDS->Initialize(numberOfLevels,&blocksPerLevel[0]);
    for(unsigned int level = 0; level<numberOfLevels; level++)
      {
      int globalIndex = globalIndices[level];
      int numberOfDataSets = static_cast<int>(datasets[level].size());
      // put the datasets back starting at globalIndex.
      // All other indices are in their initialized state.
      for (int kk=0; kk < numberOfDataSets; kk++)
        {
        vtkUniformGrid* grid = datasets[level][kk];
        hbDS->SetDataSet(level,kk+globalIndex,grid);
        }
      }
    }

}

//-----------------------------------------------------------------------------
// synch data set structure
int vtkSpyPlotReader::ComputeDerivedVars(vtkCellData* data,
  vtkSpyPlotBlock *block, vtkSpyPlotUniReader *reader, const int& blockID)
{
  if ( this->ComputeDerivedVariables != 1 || this->IsAMR)
    {
    return 0;
    }

  int numberOfMaterials = reader->GetNumberOfMaterials();

  //get the mass and material volume array for each material
  vtkDataArray** materialMasses = new vtkDataArray*[numberOfMaterials];
  vtkDataArray** materialVolumeFractions = new vtkDataArray*[numberOfMaterials];

  //bit mask of which materials we have all the information for
  for ( int i=0; i < numberOfMaterials; i++)
    {
    materialMasses[i] = reader->GetMaterialMassField(blockID, i);
    materialVolumeFractions[i] = reader->GetMaterialVolumeFractionField(blockID,i);
    }

  block->SetCoordinateSystem(reader->GetCoordinateSystem());
  block->ComputeDerivedVariables(data, numberOfMaterials,
    materialMasses, materialVolumeFractions, this->DownConvertVolumeFraction);

  //cleanup memory and leave
  delete[] materialMasses;
  delete[] materialVolumeFractions;

  return 1;
}

static void createSpyPlotLevelArray(vtkCellData *cd, int size, int level)
{
  vtkDataArray *array=cd->GetArray("levels");
  if(array!=0)
    {
    cd->RemoveArray("levels"); // if this is not the first step,
    // make sure we have a clean array
    }

  array = vtkIntArray::New();
  cd->AddArray(array);
  array->Delete();

  array->SetName("levels");
  array->SetNumberOfComponents(1);
  array->SetNumberOfTuples(size);
  int *ptr=static_cast<int *>(array->GetVoidPointer(0));
  int i;
  for (i=0; i < size; i++)
    {
    ptr[i]=level;
    }
}

