text: /*=========================================================================
  
  Program:   Visualization Toolkit
  Module:    vtkTemporalStreamTracer.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 "vtkTemporalStreamTracer.h"
  
  #include "vtkCellArray.h"
  #include "vtkCellData.h"
  #include "vtkCompositeDataIterator.h"
  #include "vtkCompositeDataPipeline.h"
  #include "vtkDataSetAttributes.h"
  #include "vtkDoubleArray.h"
  #include "vtkExecutive.h"
  #include "vtkGenericCell.h"
  #include "vtkIdList.h"
  #include "vtkInformation.h"
  #include "vtkInformationVector.h"
  #include "vtkIntArray.h"
  #include "vtkFloatArray.h"
  #include "vtkDoubleArray.h"
  #include "vtkCharArray.h"
  #include "vtkMath.h"
  #include "vtkMultiBlockDataSet.h"
  #include "vtkMultiProcessController.h"
  #include "vtkObjectFactory.h"
  #include "vtkPointData.h"
  #include "vtkPointSet.h"
  #include "vtkPolyData.h"
  #include "vtkPolyLine.h"
  #include "vtkRungeKutta2.h"
  #include "vtkRungeKutta4.h"
  #include "vtkRungeKutta45.h"
  #include "vtkSmartPointer.h"
  #include "vtkTemporalInterpolatedVelocityField.h"
  #include "vtkTemporalDataSet.h"
  #include "vtkOutputWindow.h"
  #include "vtkAbstractParticleWriter.h"
  #include "vtkToolkits.h" // For VTK_USE_MPI 
  
  #ifdef WIN32
    #undef JB_H5PART_PARTICLE_OUTPUT
  #else
  //  #define JB_H5PART_PARTICLE_OUTPUT
  #endif
  
  #ifdef JB_H5PART_PARTICLE_OUTPUT
  // #include "vtkH5PartWriter.h"
    #include "vtkXMLParticleWriter.h"
  #endif
  
  #ifdef VTK_USE_MPI
    #include "vtkMPIController.h"
  #endif
  
  #include <functional>
  #include <algorithm>
  
  using namespace vtkTemporalStreamTracerNamespace;
  
  //----------------------------------------------------------------------------
  //----------------------------------------------------------------------------
  //#define JB_DEBUG__
  #if defined JB_DEBUG__
    #ifdef WIN32
        #define OUTPUTTEXT(a) vtkOutputWindowDisplayText(a);
    #else
    #endif
  
      #undef vtkDebugMacro
      #define vtkDebugMacro(a)  \
      { \
        vtkOStreamWrapper::EndlType endl; \
        vtkOStreamWrapper::UseEndl(endl); \
        vtkOStrStreamWrapper vtkmsg; \
        vtkmsg << "P(" << this->UpdatePiece << "): " a << "\n"; \
        OUTPUTTEXT(vtkmsg.str()); \
        vtkmsg.rdbuf()->freeze(0); \
      }
  
    #undef vtkErrorMacro
    #define vtkErrorMacro(a) vtkDebugMacro(a)  
  
  #endif
  //---------------------------------------------------------------------------
  vtkStandardNewMacro(vtkTemporalStreamTracer);
  vtkCxxSetObjectMacro(vtkTemporalStreamTracer, Controller, vtkMultiProcessController);
  vtkCxxSetObjectMacro(vtkTemporalStreamTracer, ParticleWriter, vtkAbstractParticleWriter);
  //---------------------------------------------------------------------------
  vtkTemporalStreamTracer::vtkTemporalStreamTracer()
  {
    this->IntegrationDirection        = FORWARD;
    this->TimeStep                    = 0;
    this->ActualTimeStep              = 0;
    this->NumberOfInputTimeSteps      = 0;
    this->ForceReinjectionEveryNSteps = 1;
    this->ReinjectionFlag             = 0;
    this->ReinjectionCounter          = 0;
    this->UpdatePiece                 = 0;
    this->UpdateNumPieces             = 0;
    this->AllFixedGeometry            = 1;
    this->StaticMesh                  = 1;
    this->StaticSeeds                 = 1;
    this->ComputeVorticity            = 1;
    this->IgnorePipelineTime          = 0;
    this->ParticleWriter              = NULL;
    this->ParticleFileName            = NULL;
    this->EnableParticleWriting       = false;
    this->UniqueIdCounter             = 0;
    this->UniqueIdCounterMPI          = 0;
    this->InterpolationCount          = 0;
    //
    this->NumberOfParticles     = 0;
    this->TimeStepResolution    = 1.0;
    this->TerminationTime       = 0.0;
    this->TerminationTimeUnit   = TERMINATION_STEP_UNIT;
    this->EarliestTime          =-1E6;
    // we are not actually using these for now
    
    this->MaximumPropagation = 1.0;
    
    this->IntegrationStepUnit    = LENGTH_UNIT;
    this->MinimumIntegrationStep = 1.0E-2;
    this->MaximumIntegrationStep = 1.0;
    this->InitialIntegrationStep = 0.5;
    //
    this->Interpolator = vtkSmartPointer<vtkTemporalInterpolatedVelocityField>::New();
    //
    this->SetNumberOfInputPorts(2);
    this->Controller = NULL;
    this->SetController(vtkMultiProcessController::GetGlobalController());
  
  #ifdef JB_H5PART_PARTICLE_OUTPUT
  #ifdef WIN32
    vtkDebugMacro(<<"Setting vtkH5PartWriter");
    vtkH5PartWriter *writer = vtkH5PartWriter::New();
  #else
    vtkDebugMacro(<<"Setting vtkXMLParticleWriter");
    vtkXMLParticleWriter *writer = vtkXMLParticleWriter::New();
  #endif
    this->SetParticleWriter(writer);
    writer->Delete();
  #endif
  
    this->SetIntegratorType(RUNGE_KUTTA4);
  }
  //---------------------------------------------------------------------------
  vtkTemporalStreamTracer::~vtkTemporalStreamTracer()
  {
    this->SetController(NULL);
    this->SetParticleWriter(NULL);
    if (this->ParticleFileName)
    {
      delete []this->ParticleFileName;
      this->ParticleFileName = NULL;
    }
  }
  //----------------------------------------------------------------------------
  int vtkTemporalStreamTracer::FillInputPortInformation(
    int port, 
    vtkInformation* info)
  {
    // port 0 must be a temporal collection of any type
    // the executive should put a temporal collection in when
    // we request multiple time steps.
    if (port==0) {
      info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataObject");
      info->Set(vtkAlgorithm::INPUT_IS_REPEATABLE(), 1);
  //    info->Set(vtkAlgorithm::INPUT_IS_OPTIONAL(), 1);
    }
    else if (port==1) {
      info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataSet");
      info->Set(vtkAlgorithm::INPUT_IS_REPEATABLE(), 1);
  //    info->Set(vtkAlgorithm::INPUT_IS_OPTIONAL(), 1);
    }
    return 1;
  }
  //----------------------------------------------------------------------------
  void vtkTemporalStreamTracer::AddSourceConnection(vtkAlgorithmOutput* input)
  {
    this->AddInputConnection(1, input);
  }
  //----------------------------------------------------------------------------
  void vtkTemporalStreamTracer::RemoveAllSources()
  {
    this->SetInputConnection(1, 0);
  }
  //----------------------------------------------------------------------------
  int vtkTemporalStreamTracer::ProcessRequest(
    vtkInformation* request,
    vtkInformationVector** inputVector,
    vtkInformationVector* outputVector)
  {
    if(request->Has(vtkDemandDrivenPipeline::REQUEST_INFORMATION()))
      {
      return this->RequestInformation(request, inputVector, outputVector);       
      }
    if(request->Has(vtkStreamingDemandDrivenPipeline::REQUEST_UPDATE_EXTENT()))
      {
        return this->RequestUpdateExtent(request, inputVector, outputVector);
      }
    if(request->Has(vtkDemandDrivenPipeline::REQUEST_DATA()))
      {
      return this->RequestData(request, inputVector, outputVector);       
      }
    return this->Superclass::ProcessRequest(request, inputVector, outputVector);
  }
  //----------------------------------------------------------------------------
  int vtkTemporalStreamTracer::RequestInformation(
    vtkInformation *vtkNotUsed(request),
    vtkInformationVector **inputVector,
    vtkInformationVector *outputVector)
  {
    vtkInformation *inInfo  = inputVector[0]->GetInformationObject(0);
    vtkInformation *outInfo = outputVector->GetInformationObject(0);
  
    if (inInfo->Has(vtkStreamingDemandDrivenPipeline::TIME_STEPS()) )
      {
      this->NumberOfInputTimeSteps = inInfo->Length( 
        vtkStreamingDemandDrivenPipeline::TIME_STEPS() );
      vtkDebugMacro(<<"vtkTemporalStreamTracer "
        "inputVector TIME_STEPS " << this->NumberOfInputTimeSteps);
      //
      // Get list of input time step values
      this->InputTimeValues.resize(this->NumberOfInputTimeSteps);
      inInfo->Get( vtkStreamingDemandDrivenPipeline::TIME_STEPS(), 
        &this->InputTimeValues[0] );
      if (this->NumberOfInputTimeSteps==1) 
        { 
        vtkErrorMacro(<<"Not enough input time steps for particle integration");
        return 0;
        }
      //
      // We only output T-1 time steps   
      //
      this->OutputTimeValues.resize(this->NumberOfInputTimeSteps-1);
      this->OutputTimeValues.clear();
      this->OutputTimeValues.insert(
        this->OutputTimeValues.begin(), 
        this->InputTimeValues.begin()+1, this->InputTimeValues.end());
    }
    else 
      {
      vtkErrorMacro(<<"Input information has no TIME_STEPS set");
      return 0;
      }
    
    outInfo->Set(
      vtkStreamingDemandDrivenPipeline::MAXIMUM_NUMBER_OF_PIECES(), -1);
  
    outInfo->Set(vtkStreamingDemandDrivenPipeline::TIME_STEPS(), 
                 &this->OutputTimeValues[0],
                 static_cast<int>(this->OutputTimeValues.size()));
  
  return 1;
  }
  //----------------------------------------------------------------------------
  class WithinTolerance: public vtkstd::binary_function<double, double, bool>
  {
  public:
      result_type operator()(first_argument_type a, second_argument_type b) const
      {
        bool result = (fabs(a-b)<=(a*1E-6));
        return (result_type)result;
      }
  };
  //----------------------------------------------------------------------------
  int vtkTemporalStreamTracer::RequestUpdateExtent(
    vtkInformation *vtkNotUsed(request),
    vtkInformationVector **inputVector,
    vtkInformationVector *outputVector)
  {
    int numInputs = inputVector[0]->GetNumberOfInformationObjects();
    vtkInformation *outInfo = outputVector->GetInformationObject(0);
  
    //
    // The output has requested a time value, what times must we ask from our input
    //
    double requestedTimeValue;
  
    if (this->IgnorePipelineTime || !outInfo->Has(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEPS())) 
      {
      //
      // ideally we want the output information to be requesting a time step,
      // but since it isn't we must use the SetTimeStep value as a Time request
      //
      if (this->TimeStep<this->OutputTimeValues.size()) {
        requestedTimeValue = this->OutputTimeValues[this->TimeStep];
      } else {
        requestedTimeValue = this->OutputTimeValues.back();
      }
      this->ActualTimeStep = this->TimeStep;
  
      vtkDebugMacro(<< "SetTimeStep       : requestedTimeValue " 
        << requestedTimeValue << " ActualTimeStep " << this->ActualTimeStep);
      } 
    else 
      {
      //
      // Get the requested time step. 
      //
      double *requestedTimeValues = outInfo->Get(
        vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEPS());
      requestedTimeValue = requestedTimeValues[0];
      this->ActualTimeStep = vtkstd::find_if(
        this->OutputTimeValues.begin(), 
        this->OutputTimeValues.end(), 
        vtkstd::bind2nd( WithinTolerance( ), requestedTimeValue )) 
        - this->OutputTimeValues.begin();
      if (this->ActualTimeStep>=this->OutputTimeValues.size())
      {
        this->ActualTimeStep = 0;
      }
      vtkDebugMacro(<< "UPDATE_TIME_STEPS : requestedTimeValue " 
        << requestedTimeValue << " ActualTimeStep " << this->ActualTimeStep);
      }
  
    if (this->ActualTimeStep<this->OutputTimeValues.size()) 
      {
      for (int i=0; i<numInputs; i++) {
        vtkInformation *inInfo = inputVector[0]->GetInformationObject(i);
        // our output timestep T is timestep T+1 in the source
        // so output inputTimeSteps[T], inputTimeSteps[T+1]
        inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEPS(), 
          &this->InputTimeValues[this->ActualTimeStep], 2);
        vtkDebugMacro(<< "requested 2 time values : " 
          << this->InputTimeValues[this->ActualTimeStep] << " "  
          << this->InputTimeValues[this->ActualTimeStep+1]);
        }
      }
    else
    {
      vtkDebugMacro(<< "UPDATE_TIME_STEPS : Error getting requested time step");
      return 0;
    }
  
    return 1;
  }  
  //---------------------------------------------------------------------------
  int vtkTemporalStreamTracer::InitializeInterpolator()
  {
    if (!this->InputDataT[0] || !this->InputDataT[1]) {
      return 0;  
    }
    //
    // When Multiblock arrays are processed, some may be empty
    // if the first is empty, we won't find the correct vector name
    // so scan until we get one
    //
    vtkSmartPointer<vtkCompositeDataIterator> iterP;
    iterP.TakeReference(this->InputDataT[0]->NewIterator());
    iterP->GoToFirstItem();
    char *vecname = NULL;
    while (!iterP->IsDoneWithTraversal())
      {
      vtkDataArray *vectors = this->GetInputArrayToProcess(
        0,iterP->GetCurrentDataObject());
      if (vectors)
        {
        vecname = vectors->GetName();
        break;
        }
      iterP->GoToNextItem();
      }
    if (!vecname)
      {
      vtkDebugMacro(<< "Couldn't find vector array " << vecname);
      return VTK_ERROR;
      }
  
    vtkDebugMacro(<< "Interpolator using array " << vecname);
    this->Interpolator->SelectVectors(vecname);
   
    this->AllFixedGeometry = 1;
  
    int numValidInputBlocks[2] = {0, 0};
    int numTotalInputBlocks[2] = {0, 0};
    this->DataReferenceT[0] = this->DataReferenceT[1] = 0;
    for (int T=0; T<2; T++) {
      this->CachedBounds[T].clear();
      int index = 0;
      // iterate over all blocks of input and cache the bounds information
      // and determine fixed/dynamic mesh status.
  
      vtkSmartPointer<vtkCompositeDataIterator> anotherIterP;
      anotherIterP.TakeReference(this->InputDataT[T]->NewIterator());
      anotherIterP->GoToFirstItem();
      while (!anotherIterP->IsDoneWithTraversal())
        {
        numTotalInputBlocks[T]++;
        vtkDataSet* inp = vtkDataSet::SafeDownCast(anotherIterP->GetCurrentDataObject());
        if (inp)
          {
          if (inp->GetNumberOfCells()==0)
            {
            vtkDebugMacro("Skipping an empty dataset");
            }
          else if (!inp->GetPointData()->GetVectors(vecname) && inp->GetNumberOfPoints()>0)
            {
            vtkDebugMacro("One of the input datasets has no velocity vector.");
            }
          else 
            {
            // vtkDebugMacro("pass " << i << " Found dataset with " << inp->GetNumberOfCells() << " cells");
            //
            // store the bounding boxes of each local dataset for faster 'point-in-dataset' testing
            //
            bounds bbox;
            inp->ComputeBounds();
            inp->GetBounds(&bbox.b[0]);
            this->CachedBounds[T].push_back(bbox);
            bool static_dataset = this->StaticMesh || inp->GetInformation()->Has(vtkDataObject::DATA_GEOMETRY_UNMODIFIED());
            this->AllFixedGeometry = this->AllFixedGeometry && static_dataset;
            // add the dataset to the interpolator
            this->Interpolator->SetDataSetAtTime(index++, T, this->CurrentTimeSteps[T], inp, static_dataset);
            if (!this->DataReferenceT[T]) {
              this->DataReferenceT[T] = inp;
            }
            //
            numValidInputBlocks[T]++;
            }
          }
        anotherIterP->GoToNextItem();
        }
    }
    if (numValidInputBlocks[0]==0 || numValidInputBlocks[1]==0)
      {
      vtkDebugMacro("Not enough inputs have been found. Can not execute." << numValidInputBlocks[0] << " " << numValidInputBlocks[1]);
      return VTK_ERROR;
      }
    if (numValidInputBlocks[0] != numValidInputBlocks[1])
      {
      vtkDebugMacro("The number of datasets is different between time steps " << numValidInputBlocks[0] << " " << numValidInputBlocks[1]);
      return VTK_ERROR;
      }
    //
    vtkDebugMacro("Number of Valid input blocks is " << numValidInputBlocks[0] << " from " << numTotalInputBlocks[0]);
    vtkDebugMacro("AllFixedGeometry " << this->AllFixedGeometry);
  
    // force optimizations if StaticMesh is set.
    if (this->StaticMesh) {
      vtkDebugMacro("Static Mesh optimizations Forced ON");
      this->AllFixedGeometry = 1;
    }
  
    //
    return VTK_OK;
  }
  //---------------------------------------------------------------------------
  int vtkTemporalStreamTracer::AddTemporalInput(vtkTemporalDataSet *td)
  {
    if (td->GetNumberOfTimeSteps()<2)
    {
      vtkErrorMacro(<<"Temporal dataset did not have 2 time steps");
      return 0;
    }
    vtkDataObject *input[2];
    input[0] = td->GetTimeStep(0);
    input[1] = td->GetTimeStep(1);
    for (int i=0; i<2; i++) 
      {
      // if not set, create a multiblock dataset to hold all input blocks
      if (!this->InputDataT[i]) 
      {
        this->InputDataT[i] = vtkSmartPointer<vtkMultiBlockDataSet>::New();
      }
      // if simple dataset, add to our list, otherwise, add blocks
      vtkDataSet           *dsInput = vtkDataSet::SafeDownCast(input[i]);
      vtkMultiBlockDataSet *mbInput = vtkMultiBlockDataSet::SafeDownCast(input[i]);
      if (dsInput)
        {
        vtkSmartPointer<vtkDataSet> copy;
        copy.TakeReference(dsInput->NewInstance());
        copy->ShallowCopy(dsInput);
        this->InputDataT[i]->SetBlock(this->InputDataT[i]->GetNumberOfBlocks(), copy);
        }
      else if (mbInput) 
        {
        vtkSmartPointer<vtkCompositeDataIterator> iter;
        iter.TakeReference(mbInput->NewIterator());
        for (iter->InitTraversal(); !iter->IsDoneWithTraversal(); iter->GoToNextItem())
          {
          vtkDataSet *ds = vtkDataSet::SafeDownCast(iter->GetCurrentDataObject());
          if (ds) 
            {
            vtkSmartPointer<vtkDataSet> copy;
            copy.TakeReference(ds->NewInstance());
            copy->ShallowCopy(ds);
            if (ds->GetInformation()->Has(vtkDataObject::DATA_GEOMETRY_UNMODIFIED()))
              {
              copy->GetInformation()->Set(vtkDataObject::DATA_GEOMETRY_UNMODIFIED(),1);
              }
            this->InputDataT[i]->SetBlock(this->InputDataT[i]->GetNumberOfBlocks(), copy);
            }
          }
        }
      else
        {
        vtkDebugMacro("This filter cannot handle inputs of type: "
                      << (input[i]?input[i]->GetClassName():"(none)"));
        return 0;
        }
      }
    return 1;
  }
  //---------------------------------------------------------------------------
  bool vtkTemporalStreamTracer::InsideBounds(double point[]) 
  {
    double delta[3] = { 0.0, 0.0, 0.0 };
    for (int t=0; t<2; ++t) {
      for (unsigned int i=0; i<(this->CachedBounds[t].size()); ++i) {
        if (vtkMath::PointIsWithinBounds(point, 
          &((this->CachedBounds[t])[i].b[0]), delta)) return true;
      }
    }
    return false;
  }
  //---------------------------------------------------------------------------
  void vtkTemporalStreamTracer::TestParticles(
    ParticleVector &candidates, ParticleVector &passed, int &count)
  {
    int i = 0;
    int div = static_cast<int>(candidates.size()/10);
    if(div==0)
      {
        div=1;
      }
    count = 0;
    for (ParticleIterator it=candidates.begin(); it!=candidates.end(); ++it, ++i) {
      ParticleInformation &info = (*it);
      double *pos = &info.CurrentPosition.x[0];
      // if outside bounds, reject instantly
      if (this->InsideBounds(pos)) {
        if (info.UniqueParticleId==602) {
          vtkDebugMacro(<< "TestParticles got 602");
        }
        // since this is first test, avoid bad cache tests
        this->Interpolator->ClearCache();
        info.LocationState = this->Interpolator->TestPoint(pos);
        if (info.LocationState==ID_OUTSIDE_ALL /*|| location==ID_OUTSIDE_T0*/) {
          // can't really use this particle.
          vtkDebugMacro(<< "TestParticles rejected particle");
        }
        else {
          // get the cached ids and datasets from the TestPoint call
          this->Interpolator->GetCachedCellIds(info.CachedCellId, info.CachedDataSetId);
          passed.push_back(info);
          count++;
        }
      }
      if (i%div==0) {
  //      vtkDebugMacro(<< "TestParticles " << i);
      }
    }
  }
  //---------------------------------------------------------------------------
  void vtkTemporalStreamTracer::AssignSeedsToProcessors(
    vtkDataSet *source, int sourceID, int ptId, 
    ParticleVector &LocalSeedPoints, int &LocalAssignedCount)
  {
    ParticleVector candidates;
    //
    // take points from the source object and create a particle list
    //
    int numSeeds = source->GetNumberOfPoints();
    int numTested = numSeeds;
    candidates.resize(numSeeds);
    //
    for (int i=0; i<numSeeds; i++) {
      ParticleInformation &info = candidates[i];
      memcpy(&(info.CurrentPosition.x[0]), source->GetPoint(i), sizeof(double)*3);
      info.CurrentPosition.x[3] = this->CurrentTimeSteps[0];
      info.LocationState        = 0;
      info.CachedCellId[0]      =-1;
      info.CachedCellId[1]      =-1;
      info.CachedDataSetId[0]   = 0;
      info.CachedDataSetId[1]   = 0;
      info.SourceID             = sourceID;
      info.InjectedPointId      = i+ptId;
      info.InjectedStepId       = this->ReinjectionCounter;
      info.TimeStepAge          = 0;
      info.UniqueParticleId     =-1;
      info.rotation             = 0.0;
      info.angularVel           = 0.0;
      info.time                 = 0.0;
      info.age                  = 0.0;
      info.speed                = 0.0;
      info.ErrorCode            = 0;
    } 
    //
    // Gather all Seeds to all processors for classification
    //
  #ifdef VTK_USE_MPI
    // TODO : can we just use the same array here for send and receive
    ParticleVector allCandidates; 
    if (this->UpdateNumPieces>1) {
      // Gather all seed particles to all processes
      this->TransmitReceiveParticles(candidates, allCandidates, false);
      numTested = static_cast<int>(allCandidates.size());
      vtkDebugMacro(<< "Local Particles " << numSeeds << " TransmitReceive Total " << numTested);
      // Test to see which ones belong to us
      this->TestParticles(allCandidates, LocalSeedPoints, LocalAssignedCount);
    } 
    else {
      numTested = static_cast<int>(candidates.size());
      this->TestParticles(candidates, LocalSeedPoints, LocalAssignedCount);
    }
    int TotalAssigned = 0; 
    this->Controller->Reduce(&LocalAssignedCount, &TotalAssigned, 1, vtkCommunicator::SUM_OP, 0);
  #else 
    numTested = static_cast<int>(candidates.size());
    this->TestParticles(candidates, LocalSeedPoints, LocalAssignedCount);
    int TotalAssigned = LocalAssignedCount; (void)TotalAssigned;
  #endif
    // Assign unique identifiers taking into account uneven distribution 
    // across processes and seeds which were rejected
    this->AssignUniqueIds(LocalSeedPoints);
    //
    vtkDebugMacro(<< "Tested " << numTested << " LocallyAssigned " << LocalAssignedCount);
    if (this->UpdatePiece==0) {
      vtkDebugMacro(<< "Total Assigned to all processes " << TotalAssigned);
    }
  }
  //---------------------------------------------------------------------------
  void vtkTemporalStreamTracer::AssignUniqueIds(
    vtkTemporalStreamTracerNamespace::ParticleVector &LocalSeedPoints)
  {
    vtkIdType ParticleCountOffset = 0;
    vtkIdType numParticles = LocalSeedPoints.size();
    if (this->UpdateNumPieces>1) {
  #ifdef VTK_USE_MPI
      vtkMPICommunicator* com = vtkMPICommunicator::SafeDownCast(
        this->Controller->GetCommunicator()); 
      if (com == 0) {
        vtkErrorMacro("MPICommunicator needed for this operation.");
        return;
      }
      // everyone starts with the master index
      com->Broadcast(&this->UniqueIdCounter, 1, 0);
  //    vtkErrorMacro("UniqueIdCounter " << this->UniqueIdCounter);
      // setup arrays used by the AllGather call.
      vtkstd::vector<vtkIdType> recvNumParticles(this->UpdateNumPieces, 0);
      // Broadcast and receive count to/from all other processes.
      com->AllGather(&numParticles, &recvNumParticles[0], 1);
      // Each process is allocating a certain number.
      // start our indices from sum[0,this->UpdatePiece](numparticles)
      for (int i=0; i<this->UpdatePiece; ++i) {
        ParticleCountOffset += recvNumParticles[i];
      }
      for (vtkIdType i=0; i<numParticles; i++) {
        LocalSeedPoints[i].UniqueParticleId = 
          this->UniqueIdCounter + ParticleCountOffset + i;
      }
      for (int i=0; i<this->UpdateNumPieces; ++i) {
        this->UniqueIdCounter += recvNumParticles[i];
      }
  #endif
    }
    else {
      for (vtkIdType i=0; i<numParticles; i++) {
        LocalSeedPoints[i].UniqueParticleId = 
          this->UniqueIdCounter + ParticleCountOffset + i;
      }
      this->UniqueIdCounter += numParticles;
    }
  }
  //---------------------------------------------------------------------------
  #ifdef VTK_USE_MPI
  void vtkTemporalStreamTracer::TransmitReceiveParticles(
    ParticleVector &sending, ParticleVector &received, bool removeself)
  {
    vtkMPICommunicator* com = vtkMPICommunicator::SafeDownCast(
      this->Controller->GetCommunicator()); 
    if (com == 0) {
      vtkErrorMacro("MPICommunicator needed for this operation.");
      return;
    }
    // 
    // We must allocate buffers for all processor particles
    // 
    vtkIdType OurParticles = sending.size();
    vtkIdType TotalParticles = 0;
    // setup arrays used by the AllGatherV call.
    vtkstd::vector<vtkIdType> recvLengths(this->UpdateNumPieces, 0);
    vtkstd::vector<vtkIdType> recvOffsets(this->UpdateNumPieces, 0);
    // Broadcast and receive size to/from all other processes.
    com->AllGather(&OurParticles, &recvLengths[0], 1);
    // Compute the displacements.
    const vtkIdType TypeSize = sizeof(ParticleInformation);
    for (int i=0; i<this->UpdateNumPieces; ++i)
    {
      //  << i << ": " << recvLengths[i] << "   ";
      recvOffsets[i] = TotalParticles*TypeSize;
      TotalParticles += recvLengths[i];
      recvLengths[i] *= TypeSize;
    }
    //  << '\n';
    // Allocate the space for all particles
    received.resize(TotalParticles);
    if (TotalParticles==0) return;
    // Gather the data from all procs.
    char *sendbuf = (char*) ((sending.size()>0) ? &(sending[0]) : NULL);
    char *recvbuf = (char*) (&(received[0]));
    com->AllGatherV(sendbuf, recvbuf, 
      OurParticles*TypeSize, &recvLengths[0], &recvOffsets[0]);
    // Now all particles from all processors are in one big array
    // remove any from ourself that we have already tested
    if (removeself) {
      vtkstd::vector<ParticleInformation>::iterator first = 
        received.begin() + recvOffsets[this->UpdatePiece]/TypeSize;
      vtkstd::vector<ParticleInformation>::iterator last = 
        first + recvLengths[this->UpdatePiece]/TypeSize;
      received.erase(first, last);
    }
    if (received.size()>0) {
      TotalParticles=TotalParticles; // brkpnt 
    }
  }
  #else // VTK_USE_MPI
  void vtkTemporalStreamTracer::TransmitReceiveParticles(
    ParticleVector &, ParticleVector &, bool )
  {
  }
  #endif // VTK_USE_MPI
  //---------------------------------------------------------------------------
  void vtkTemporalStreamTracer::UpdateParticleList(ParticleVector &candidates)
  {
    int numSeedsNew = static_cast<int>(candidates.size());
    // 
    for (int i=0; i<numSeedsNew; i++) {
      // allocate a new particle on the list and get a reference to it
      this->ParticleHistories.push_back(candidates[i]);
    }
    this->NumberOfParticles = static_cast<int>(ParticleHistories.size());
  
    vtkDebugMacro(<< "UpdateParticleList completed with " << this->NumberOfParticles << " particles");
  }
  //---------------------------------------------------------------------------
  int vtkTemporalStreamTracer::RequestData(
    vtkInformation *vtkNotUsed(request),
    vtkInformationVector **inputVector,
    vtkInformationVector *outputVector)
  {
    //
    // Parallel/Piece information
    //
    vtkInformation *outInfo = outputVector->GetInformationObject(0);
    this->UpdatePiece =
      outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER());
    this->UpdateNumPieces =
      outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES());
  
    //
    // Inputs information
    //
    int numInputs = inputVector[0]->GetNumberOfInformationObjects();
    
    // inherited from streamtracer, make sure it is null
    this->InputData = NULL;
  
    // how many temporal inputs were we given, 
    // from each - add to our list of blocks
    this->InputDataT[0] = NULL;
    this->InputDataT[1] = NULL;
    for (int idx=0; idx<numInputs; ++idx)
      {
      vtkDataObject      *input = 0;
      vtkTemporalDataSet    *td = 0;
      vtkInformation    *inInfo = inputVector[0]->GetInformationObject(idx);
      if (inInfo)
        {
        input = inInfo->Get(vtkDataObject::DATA_OBJECT());
        td    = vtkTemporalDataSet::SafeDownCast(input);
        if (td && this->AddTemporalInput(td)) {
          {
          //
          // Get the timestep information for this instant
          //
          vtkInformation *doInfo = td->GetInformation();
          vtkstd::vector<double> timesteps;
          if (doInfo->Has(vtkDataObject::DATA_TIME_STEPS()))
          {
            int NumberOfDataTimeSteps = doInfo->Length(vtkDataObject::DATA_TIME_STEPS());
            if (NumberOfDataTimeSteps<2) {
              vtkErrorMacro(<<"Not enough time steps in input data");
              return 1;
            }
            timesteps.resize(NumberOfDataTimeSteps);
            doInfo->Get(vtkDataObject::DATA_TIME_STEPS(), &timesteps[0]);
          }
          else {
            vtkErrorMacro(<<"No time step info");
            return 1;
          }
          this->CurrentTimeSteps[0] = timesteps[0]*this->TimeStepResolution;
          this->CurrentTimeSteps[1] = timesteps[1]*this->TimeStepResolution;
          }
        }
        else 
          {
          vtkDebugMacro(<<"AddTemporalDataSet failed : Input dataset was not temporal or had insufficient timesteps");
          return 1;
          }
        }
      }
  
    //
    // How many Seed point sources are connected?
    // Copy the sources into a vector for later use
    //
  
    int numSources = inputVector[1]->GetNumberOfInformationObjects();
    vtkstd::vector<vtkDataSet*> SeedSources;
    for (int idx=0; idx<numSources; ++idx)
      {
      vtkDataObject     *dobj   = 0;
      vtkInformation    *inInfo = inputVector[1]->GetInformationObject(idx);
      if (inInfo)
        {
        dobj   = inInfo->Get(vtkDataObject::DATA_OBJECT());
        SeedSources.push_back(vtkDataSet::SafeDownCast(dobj));
        }
      }
  
    if (this->IntegrationDirection != FORWARD) 
    {
      vtkErrorMacro(<<"We can only handle forward time particle tracking at the moment");
      return 1;
    }
  
    //
    // Add the datasets to an interpolator object
    //
    if (this->InitializeInterpolator() != VTK_OK)
    {
      if (this->InputDataT[0]) this->InputDataT[0] = NULL;
      if (this->InputDataT[1]) this->InputDataT[1] = NULL;
      vtkErrorMacro(<<"InitializeInterpolator failed");
      return 1;
    }
  
    //
    // Setup some variables
    //
    vtkSmartPointer<vtkInitialValueProblemSolver> integrator;
    integrator.TakeReference(this->GetIntegrator()->NewInstance());
    integrator->SetFunctionSet(this->Interpolator);
  
    //
    // Make sure the Particle Positions are initialized with Seed particles
    //
    this->ReinjectionFlag = 0;
    if (this->ForceReinjectionEveryNSteps>0) {
      if ((this->ActualTimeStep%this->ForceReinjectionEveryNSteps)==0) 
      {
        this->ReinjectionFlag = 1;
      }
    }
    //
    // If T=0 reset everything to allow us to setup stuff then start an animation
    // with a clean slate
    //
    if (this->ActualTimeStep==0) {
      this->LocalSeeds.clear();
      this->ParticleHistories.clear();
      this->EarliestTime       =-1E6;
      this->ReinjectionFlag    = 1;
      this->ReinjectionCounter = 0;
      this->UniqueIdCounter    = 0;
      this->UniqueIdCounterMPI = 0;
    }
    else if (this->CurrentTimeSteps[0]<this->EarliestTime) {
      //
      // We don't want to go back in time, so just reuse whatever we have
      //
      vtkDebugMacro("skipping particle tracking because we have seen this timestep before");
      outInfo->Set(vtkDataObject::DATA_TIME_STEPS(), 
        &this->OutputTimeValues[this->ActualTimeStep], 1);
      if (this->InputDataT[0]) this->InputDataT[0] = NULL;
      if (this->InputDataT[1]) this->InputDataT[1] = NULL;
      return 1;
    }
    this->EarliestTime = (this->CurrentTimeSteps[0]>this->EarliestTime)
      ? this->CurrentTimeSteps[0] : this->EarliestTime;
    //
    //
    //
    for (unsigned int i=0; i<SeedSources.size(); i++) 
    {
      if (SeedSources[i]->GetMTime()>this->ParticleInjectionTime)
      {
    //    this->ReinjectionFlag = 1;
      }
    }
  
    //
    // Lists for seed particles
    //
    ParticleVector candidates;
    ParticleVector received;
    //
  
    if (this->ReinjectionFlag) {
      int seedPointId=0;
      if (this->StaticSeeds && this->AllFixedGeometry && this->LocalSeeds.size()==0) {
        for (unsigned int i=0; i<SeedSources.size(); i++) {
          this->AssignSeedsToProcessors(SeedSources[i], i, 0, this->LocalSeeds, seedPointId);
        }
      }
      else {
        // wipe the list and reclassify for each injection
        this->LocalSeeds.clear();
        for (unsigned int i=0; i<SeedSources.size(); i++) {
          this->AssignSeedsToProcessors(SeedSources[i], i, 0, this->LocalSeeds, seedPointId);
        }
      }
      this->ParticleInjectionTime.Modified();
  
      // Now update our main list with the ones we are keeping
      vtkDebugMacro(<< "Reinjection about to update candidates (" << this->LocalSeeds.size() << " particles)");
      this->UpdateParticleList(this->LocalSeeds);
      this->ReinjectionCounter += 1;
    }
  
    //
    // setup all our output arrays
    //
    vtkDebugMacro(<< "About to allocate point arrays ");
    vtkPolyData *output = vtkPolyData::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));
    this->ParticleAge         = vtkSmartPointer<vtkFloatArray>::New();
    this->ParticleIds         = vtkSmartPointer<vtkIntArray>::New();
    this->ParticleSourceIds   = vtkSmartPointer<vtkCharArray>::New();
    this->InjectedPointIds    = vtkSmartPointer<vtkIntArray>::New();
    this->InjectedStepIds     = vtkSmartPointer<vtkIntArray>::New();
    this->ErrorCode           = vtkSmartPointer<vtkIntArray>::New();
    this->ParticleVorticity   = vtkSmartPointer<vtkFloatArray>::New();
    this->ParticleRotation    = vtkSmartPointer<vtkFloatArray>::New();
    this->ParticleAngularVel  = vtkSmartPointer<vtkFloatArray>::New();
    this->cellVectors         = vtkSmartPointer<vtkDoubleArray>::New();
    this->ParticleCells       = vtkSmartPointer<vtkCellArray>::New();
    this->OutputCoordinates   = vtkSmartPointer<vtkPoints>::New();
    this->OutputPointData     = output->GetPointData();
    this->OutputPointData->Initialize();
    this->InterpolationCount  = 0;
    vtkDebugMacro(<< "About to Interpolate allocate space");
    this->OutputPointData->InterpolateAllocate(this->DataReferenceT[1]->GetPointData());
    //
    this->ParticleAge->SetName("ParticleAge");
    this->ParticleIds->SetName("ParticleId");
    this->ParticleSourceIds->SetName("ParticleSourceId");
    this->InjectedPointIds->SetName("InjectedPointId");
    this->InjectedStepIds->SetName("InjectionStepId");
    this->ErrorCode->SetName("ErrorCode");
  
    if (this->ComputeVorticity)
      {
      this->cellVectors->SetNumberOfComponents(3);
      this->cellVectors->Allocate(3*VTK_CELL_SIZE);   
      this->ParticleVorticity->SetName("Vorticity");
      this->ParticleRotation->SetName("Rotation");
      this->ParticleAngularVel->SetName("AngularVelocity");
      }
  
    output->SetPoints(this->OutputCoordinates);
    output->SetVerts(this->ParticleCells);
    vtkDebugMacro(<< "Finished allocating point arrays ");
  
    //
    // Perform 2 passes
    // Pass 0 : Integration of particles created by a source in this process 
    // or received at start from a source in another process.
    //
    // Pass 1 : Particles that were sent in mid integration from another process
    // are added in and their integration continued here. In actual fact, the process
    // should be repeated until all particles are finished, but the chances of 
    // a particle stepping inside and out again through a single domain 
    // in one time step are small (hopefully!)
  
    vtkDebugMacro(<< "Clear MPI send list ");
    this->MPISendList.clear();
  
    int Number = static_cast<int>(this->ParticleHistories.size());
    ParticleListIterator  it_first = this->ParticleHistories.begin();
    ParticleListIterator  it_last  = this->ParticleHistories.end();
    ParticleListIterator  it_next;
  #define PASSES 2
    for (int pass=0; pass<PASSES; pass++) {
      vtkDebugMacro(<<"Begin Pass " << pass << " with " << Number << " Particles");
      for (ParticleListIterator it=it_first; it!=it_last;)
      {
        // Keep the 'next' iterator handy because if a particle is terminated 
        // or leaves the domain, the 'current' iterator will be deleted.
        it_next = it;
        it_next++;
        //
        // Shall we terminate this particle
        //
        double interval = (this->CurrentTimeSteps[1]-this->CurrentTimeSteps[0]);
        bool terminated = false;
        if (this->TerminationTime>0) 
        {
          if (this->TerminationTimeUnit == TERMINATION_TIME_UNIT &&
            (it->age+interval)>this->TerminationTime) {
            terminated = true;
          }
          else if (this->TerminationTimeUnit == TERMINATION_STEP_UNIT &&
            (it->TimeStepAge+1)>this->TerminationTime) {
            terminated = true;
          }
        }
        if (terminated) {
          this->ParticleHistories.erase(it);
        }
        else {
          this->IntegrateParticle(it, this->CurrentTimeSteps[0], this->CurrentTimeSteps[1], integrator);
        }
        //
        if (this->GetAbortExecute()) {
          break;
        }
        it = it_next;
      }
      // Particles might have been deleted during the first pass as they move 
      // out of domain or age. Before adding any new particles that are sent 
      // to us, we must know the starting point ready for the second pass
      bool list_valid = (this->ParticleHistories.size()>0);
      if (list_valid) {
        // point to one before the end
        it_first = --this->ParticleHistories.end();
      }
      // Send and receive any particles which exited/entered the domain
      if (this->UpdateNumPieces>1 && pass<(PASSES-1)) {
        // the Particle lists will grow if any are received
        // so we must be very careful with our iterators
        vtkDebugMacro(<<"End of Pass " << pass << " with " 
          << this->ParticleHistories.size() << " " 
          << " about to Transmit/Receive " << this->MPISendList.size());
        this->TransmitReceiveParticles(this->MPISendList, received, true);
        // don't want the ones that we sent away
        this->MPISendList.clear();
        int assigned;
        // classify all the ones we received
        if (received.size()>0) {
          this->TestParticles(received, candidates, assigned);
          vtkDebugMacro(<<"received " << received.size() << " : assigned locally " << assigned);
          received.clear();
        }
        // Now update our main list with the ones we are keeping
        this->UpdateParticleList(candidates);
        // free up unwanted memory
        Number = static_cast<int>(candidates.size());
        candidates.clear();
      }
      it_last = this->ParticleHistories.end();
      if (list_valid) {
        // increment to point to first new entry
        it_first++;
      }
      else {
        it_first = this->ParticleHistories.begin();
      }
    }
    if (this->MPISendList.size()>0) {
      // If a particle went out of domain on the second pass, it should be sent
      // can it really pass right through a domain in one step?
      // what about grazing the edge of rotating zone?
      vtkDebugMacro(<< "MPISendList not empty " << this->MPISendList.size());
      this->MPISendList.clear();
    }
  
    //
    // We must only add these scalar arrays at the end because the 
    // existing scalars on the input get interpolated during iteration
    // over the particles
    // 
    this->OutputPointData->AddArray(this->ParticleIds);
    this->OutputPointData->AddArray(this->ParticleSourceIds);
    this->OutputPointData->AddArray(this->InjectedPointIds);
    this->OutputPointData->AddArray(this->InjectedStepIds);
    this->OutputPointData->AddArray(this->ErrorCode);
    this->OutputPointData->AddArray(this->ParticleAge);
    if (this->ComputeVorticity)
      {
      this->OutputPointData->AddArray(this->ParticleVorticity);
      this->OutputPointData->AddArray(this->ParticleRotation);
      this->OutputPointData->AddArray(this->ParticleAngularVel);
      }
  
    if (this->InterpolationCount!=this->OutputCoordinates->GetNumberOfPoints()) {
      vtkErrorMacro(<< "Mismatch in point array/data counts");
    }
    //
    outInfo->Set(vtkDataObject::DATA_TIME_STEPS(), 
      &this->OutputTimeValues[this->ActualTimeStep], 1);
  
    // save some locator building, by re-using them as time progresses
    this->Interpolator->AdvanceOneTimeStep();
  
    //
    // Let go of inputs 
    //
    if (this->InputDataT[0]) this->InputDataT[0] = NULL;;
    if (this->InputDataT[1]) this->InputDataT[1] = NULL;;
  
    //
    // Write Particles out if necessary
    //
    // NB. We don't want our writer to trigger any updates, 
    // so shallow copy the output
    if (this->ParticleWriter && this->EnableParticleWriting) {
      vtkSmartPointer<vtkPolyData> polys = vtkSmartPointer<vtkPolyData>::New();
      polys->ShallowCopy(output);
      int N = polys->GetNumberOfPoints(); (void)N;
      this->ParticleWriter->SetFileName(this->ParticleFileName);
      this->ParticleWriter->SetTimeStep(this->ActualTimeStep);
      this->ParticleWriter->SetTimeValue(this->CurrentTimeSteps[1]);
      this->ParticleWriter->SetInput(polys);
      this->ParticleWriter->Write();
      this->ParticleWriter->CloseFile();
      this->ParticleWriter->SetInput(NULL);
  #ifdef VTK_USE_MPI
      this->Controller->Barrier();
  #endif
      vtkDebugMacro(<< "Written " << N);
    }
  //  this->Interpolator->ShowCacheResults();
  //  vtkErrorMacro(<<"RequestData done");  
    return 1;
  }
  //---------------------------------------------------------------------------
  void vtkTemporalStreamTracer::IntegrateParticle(
    ParticleListIterator &it, 
    double currenttime, double targettime,
    vtkInitialValueProblemSolver* integrator)
  {
    double epsilon = (targettime-currenttime)/100.0;
    double velocity[3], point1[4], point2[4] = {0.0, 0.0, 0.0, 0.0};
    double minStep=0, maxStep=0;
    double stepWanted, stepTaken=0.0;
    substeps = 0;
    
    ParticleInformation &info = (*it);
    // Get the Initial point {x,y,z,t}
    memcpy(point1, &info.CurrentPosition, sizeof(Position));
  
    if (point1[3]<(currenttime-epsilon) || point1[3]>(targettime+epsilon)) {
      vtkDebugMacro(<< "Bad particle time : expected (" 
        << this->CurrentTimeSteps[0] << "-" << this->CurrentTimeSteps[1] 
        << ") got " << point1[3]);
    }
  
    IntervalInformation delT;
    delT.Unit     = LENGTH_UNIT;
    delT.Interval = (targettime-currenttime) * this->InitialIntegrationStep;
    epsilon = delT.Interval*1E-3;
  
    //
    // begin interpolation between available time values, if the particle has 
    // a cached cell ID and dataset - try to use it, 
    //
    this->Interpolator->SetCachedCellIds(info.CachedCellId, info.CachedDataSetId);
  
    bool particle_good = true;
    info.ErrorCode = 0;
    while (point1[3] < (targettime-epsilon)) {
      //
      // Here beginneth the real work
      //
      double error = 0;
  
      // If, with the next step, propagation will be larger than
      // max, reduce it so that it is (approximately) equal to max.
      stepWanted = delT.Interval;
      if ( (point1[3] + stepWanted) > targettime )
        {
        stepWanted = targettime - point1[3];
        maxStep = stepWanted;
        }
      this->LastUsedStepSize = stepWanted;
  
      // Calculate the next step using the integrator provided.
      // If the next point is out of bounds, send it to another process
      if (integrator->ComputeNextStep(
            point1, point2, point1[3], stepWanted, 
            stepTaken, minStep, maxStep, 
            this->MaximumError, error) != 0)
      {
        // if the particle is sent, remove it from the list
        info.ErrorCode = 1;
        if (this->SendParticleToAnotherProcess(info, point1, this->LastUsedStepSize)) {
          this->ParticleHistories.erase(it);
          particle_good = false;
          break;
        } 
        else 
        {
          // particle was not sent, retry saved it, so copy info back
          substeps++;
          memcpy(point1, &info.CurrentPosition, sizeof(Position));  
        }
      }
      else // success, increment position/time
      {
        substeps++;
  
        // increment the particle time
        point2[3] = point1[3] + stepTaken;
        info.age += stepTaken;
  
        // Point is valid. Insert it.
        memcpy(&info.CurrentPosition, point2, sizeof(Position));
        memcpy(point1, point2, sizeof(Position));
      }
  
      // If the solver is adaptive and the next time step (delT.Interval)
      // that the solver wants to use is smaller than minStep or larger 
      // than maxStep, re-adjust it. This has to be done every step
      // because minStep and maxStep can change depending on the Cell
      // size (unless it is specified in time units)
      if (integrator->IsAdaptive())
        {
          // code removed. Put it back when this is stable
        }
    }
    if (particle_good) {
      // The integration succeeded, but check the computed final position 
      // is actually inside the domain (the intermediate steps taken inside 
      // the integrator were ok, but the final step may just pass out)
      // if it moves out, we can't interpolate scalars, so we must send it away
      info.LocationState = this->Interpolator->TestPoint(info.CurrentPosition.x);
      if (info.LocationState==ID_OUTSIDE_ALL)
      {
        info.ErrorCode = 2;
        // if the particle is sent, remove it from the list
        if (this->SendParticleToAnotherProcess(info, point1, this->LastUsedStepSize)) {
          this->ParticleHistories.erase(it);
          particle_good = false;
        }
      }
    }
  
    //
    // Has this particle stagnated
    //
    if (particle_good) {
      this->Interpolator->GetLastGoodVelocity(velocity);
      info.speed = vtkMath::Norm(velocity);
      if (it->speed <= this->TerminalSpeed)
      {
        this->ParticleHistories.erase(it);
        particle_good = false;
      }
    }
  
    //
    // We got this far without error :
    // Insert the point into the output
    // Create any new scalars and interpolate existing ones
    // Cache cell ids and datasets
    //
    if (particle_good) {
      //
      // store the last Cell Ids and dataset indices for next time particle is updated
      //
      this->Interpolator->GetCachedCellIds(info.CachedCellId, info.CachedDataSetId);
      //
      info.TimeStepAge += 1;
      //
      // Now generate the output geometry and scalars
      //
      double    *coord = &info.CurrentPosition.x[0];
      vtkIdType tempId = this->OutputCoordinates->InsertNextPoint(coord);
      // create the cell
      this->ParticleCells->InsertNextCell(1, &tempId);
      // set the easy scalars for this particle
      this->ParticleIds->InsertNextValue(info.UniqueParticleId);
      this->ParticleSourceIds->InsertNextValue(info.SourceID);
      this->InjectedPointIds->InsertNextValue(info.InjectedPointId);
      this->InjectedStepIds->InsertNextValue(info.InjectedStepId);
      this->ErrorCode->InsertNextValue(info.ErrorCode);
      this->ParticleAge->InsertNextValue(info.age);
      //
      // Interpolate all existing point attributes
      // In principle we always integrate the particle until it reaches Time2 
      // - so we don't need to do any interpolation of the scalars 
      // between T0 and T1, just fetch the values
      // of the spatially interpolated scalars from T1.
      //
      if (info.LocationState==ID_OUTSIDE_T1) {
        this->Interpolator->InterpolatePoint(0, this->OutputPointData, tempId);
      }
      else {
        this->Interpolator->InterpolatePoint(1, this->OutputPointData, tempId);
      }
      this->InterpolationCount++;
      //
      // Compute vorticity 
      //
      if (this->ComputeVorticity)
      {
        vtkGenericCell *cell;
        double pcoords[3], vorticity[3], weights[256];
        double rotation, omega;
        // have to use T0 if particle is out at T1, otherwise use T1
        if (info.LocationState==ID_OUTSIDE_T1) {
          this->Interpolator->GetVorticityData(
            0, pcoords, weights, cell, this->cellVectors);
        }
        else {
          this->Interpolator->GetVorticityData(
            1, pcoords, weights, cell, this->cellVectors);
        }
        vtkStreamTracer::CalculateVorticity(cell, pcoords, cellVectors, vorticity);
        this->ParticleVorticity->InsertNextTuple(vorticity);
        // local rotation = vorticity . unit tangent ( i.e. velocity/speed )
        if (info.speed != 0.0) {
          omega = vtkMath::Dot(vorticity, velocity);
          omega /= info.speed;
          omega *= this->RotationScale;
        }
        else {
          omega = 0.0;
        }     
        vtkIdType index = this->ParticleAngularVel->InsertNextValue(omega);
        if (index>0) {     
          rotation     = info.rotation + (info.angularVel + omega)/2 * (info.CurrentPosition.x[3] - info.time);
        } 
        else {
          rotation     = 0.0;
        }
        this->ParticleRotation->InsertNextValue(rotation);
        info.rotation   = rotation;
        info.angularVel = omega;
        info.time       = info.CurrentPosition.x[3];
      }
    }
    else this->Interpolator->ClearCache();
  
    double eps = (this->CurrentTimeSteps[1]-this->CurrentTimeSteps[0])/100;
    if (point1[3]<(this->CurrentTimeSteps[0]-eps) || point1[3]>(this->CurrentTimeSteps[1]+eps)) {
      vtkDebugMacro(<< "Unexpected time ending IntegrateParticle - expected (" 
        << this->CurrentTimeSteps[0] << "-" << this->CurrentTimeSteps[1] << ") got " << point1[3]);  
    }
  }
  //---------------------------------------------------------------------------
  bool vtkTemporalStreamTracer::RetryWithPush(ParticleInformation &info, double velocity[3], double delT)
  {
    // try adding a one increment push to the particle to get over a rotating/moving boundary
    for (int v=0; v<3; v++) info.CurrentPosition.x[v] += velocity[v]*delT;
    info.CurrentPosition.x[3] += delT;
    info.LocationState = this->Interpolator->TestPoint(info.CurrentPosition.x);
    if (info.LocationState!=ID_OUTSIDE_ALL) {
      // a push helped the particle get back into a dataset, 
      info.age += delT;
      info.ErrorCode = 6;
     return 1;
    }
    return 0;
  }
  //---------------------------------------------------------------------------
  bool vtkTemporalStreamTracer::SendParticleToAnotherProcess(ParticleInformation &info, double point1[4], double delT)
  {
  //  return 1;
    double velocity[3];
    this->Interpolator->ClearCache();
    if (info.UniqueParticleId==3) {
      vtkDebugMacro(<< "3 is about to be sent");
    }
    info.LocationState = this->Interpolator->TestPoint(point1);
    if (info.LocationState==ID_OUTSIDE_ALL) {
      // something is wrong, the particle has left the building completely
      // we can't get the last good velocity as it won't be valid
      // send the particle 'as is' and hope it lands in another process
      if (substeps>0) {
        this->Interpolator->GetLastGoodVelocity(velocity);
      } else {
        velocity[0] = velocity[1] = velocity[2] = 0.0;
      }
      info.ErrorCode = 3;
    }
    else if (info.LocationState==ID_OUTSIDE_T0) {
      // the particle left the volume but can be tested at T2, so use the velocity at T2
      this->Interpolator->GetLastGoodVelocity(velocity);
      info.ErrorCode = 4;
    }
    else if (info.LocationState==ID_OUTSIDE_T1) {
      // the particle left the volume but can be tested at T1, so use the velocity at T1
      this->Interpolator->GetLastGoodVelocity(velocity);
      info.ErrorCode = 5;
    }
    else {
      // The test returned INSIDE_ALL, so test failed near start of integration,
      this->Interpolator->GetLastGoodVelocity(velocity);
    }
    if (this->RetryWithPush(info, velocity, delT)) return 0;
    this->AddParticleToMPISendList(info);
    return 1;
  }
  //---------------------------------------------------------------------------
  void vtkTemporalStreamTracer::PrintSelf(ostream& os, vtkIndent indent)
  {
    this->Superclass::PrintSelf(os,indent);
  
    os << indent << "TimeStepResolution: " << this->TimeStepResolution << endl;
    os << indent << "Controller: " << this->Controller << endl;
    os << indent << "ParticleWriter: " << this->ParticleWriter << endl;
    os << indent << "ParticleFileName: " << 
      (this->ParticleFileName ? this->ParticleFileName : "None") << endl;
    os << indent << "TimeStep: " << this->TimeStep << endl;
    os << indent << "ForceReinjectionEveryNSteps: "
       << this->ForceReinjectionEveryNSteps << endl;
    os << indent << "EnableParticleWriting: " << this->EnableParticleWriting << endl;
    os << indent << "IgnorePipelineTime: " << this->IgnorePipelineTime << endl;
    os << indent << "StaticMesh: " << this->StaticMesh << endl;
    os << indent << "TerminationTime: " << this->TerminationTime << endl;
    os << indent << "TerminationTimeUnit: " << this->TerminationTimeUnit << endl;
    os << indent << "StaticSeeds: " << this->StaticSeeds << endl;
  }
  //---------------------------------------------------------------------------
  bool vtkTemporalStreamTracer::ComputeDomainExitLocation(
    double pos[4], double p2[4], double intersection[4],
    vtkGenericCell *cell)
  {
    double t, pcoords[3];
    int subId;
    if (cell->IntersectWithLine(pos, p2, 1E-3, t, intersection, pcoords, subId)==0) {
      vtkDebugMacro(<< "No cell/domain exit was found");
      return 0;
    }
    else {
      // We found an intersection on the edge of the cell. 
      // Shift it by a small amount to ensure that it crosses over the edge 
      // into the adjoining cell.
      for (int i=0; i<3; i++) intersection[i] = pos[i] + (t+0.01)*(p2[i]-pos[i]);
      // intersection stored, compute T for intersection
      intersection[3] = pos[3] + (t+0.01)*(p2[3]-pos[3]);
      return 1;    
    }
  }
  //---------------------------------------------------------------------------
  void vtkTemporalStreamTracer::AddParticleToMPISendList(ParticleInformation &info)
  {
    double eps = (this->CurrentTimeSteps[1]-this->CurrentTimeSteps[0])/100;
    if (info.CurrentPosition.x[3]<(this->CurrentTimeSteps[0]-eps) || 
        info.CurrentPosition.x[3]>(this->CurrentTimeSteps[1]+eps)) {
      vtkDebugMacro(<< "Unexpected time value in MPISendList - expected (" 
        << this->CurrentTimeSteps[0] << "-" << this->CurrentTimeSteps[1] << ") got " 
        << info.CurrentPosition.x[3]);
    }
  #ifdef VTK_USE_MPI
    if (this->MPISendList.capacity()<(this->MPISendList.size()+1)) {
      this->MPISendList.reserve(static_cast<int>(this->MPISendList.size()*1.5));
    }
    this->MPISendList.push_back(info);
  #endif
  }
  //---------------------------------------------------------------------------
  /*
    // Now try to compute the trajectory exit location from the cell on the edge
    if (this->Interpolator->GetLastValidCellId(0)!=-1) {
      vtkDataSet *hitdata = this->Interpolator->GetLastDataSet(0);
      hitdata->GetCell(this->Interpolator->GetLastValidCellId(0), this->GenericCell);
      double intersection[4];
      if (this->ComputeDomainExitLocation(point1, point2, intersection, this->GenericCell))
      {
        vtkDebugMacro(<< "SendParticleToAnotherProcess : Sending Particle " << particleId << " Time " << intersection[3]);
        this->AddParticleToMPISendList(intersection, velocity);
        vtkIdType nextPoint = ParticleCoordinates->InsertNextPoint(point2);
        this->ParticleHistories[particleId].push_back(nextPoint);
        this->LiveParticleIds[particleId] = 0;
        return 1;
      }
      else {
        vtkDebugMacro(<< "Domain-Exit aborted : Domain Intersection failed" << particleId );
        this->LiveParticleIds[particleId] = 0;
        return 0;
      }
    }
    else {
      vtkDebugMacro(<< "Domain-Exit aborted : Couldn't copy cell from earlier test" << particleId );
      this->LiveParticleIds[particleId] = 0;
      return 0;
    }
  */
  