text: /*=========================================================================
  
    Program:   Visualization Toolkit
    Module:    vtkPStreamTracer.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 "vtkPStreamTracer.h"
  
  #include "vtkAppendPolyData.h"
  #include "vtkCellData.h"
  #include "vtkCompositeDataSet.h"
  #include "vtkFloatArray.h"
  #include "vtkIdList.h"
  #include "vtkInformation.h"
  #include "vtkInformationVector.h"
  #include "vtkIntArray.h"
  #include "vtkAbstractInterpolatedVelocityField.h"
  #include "vtkMultiProcessController.h"
  #include "vtkObjectFactory.h"
  #include "vtkPointData.h"
  #include "vtkPoints.h"
  #include "vtkPolyData.h"
  #include "vtkStreamingDemandDrivenPipeline.h"
  #include "vtkRungeKutta2.h"
  #include "vtkOverlappingAMR.h"
  #include "vtkAMRInterpolatedVelocityField.h"
  #include "vtkUniformGrid.h"
  #include "vtkParallelAMRUtilities.h"
  #include "vtkMath.h"
  #include "vtkCompositeDataIterator.h"
  #include "vtkNew.h"
  #include "vtkMultiProcessStream.h"
  #include "vtkCellArray.h"
  #include "vtkMPIController.h"
  #include "vtkCharArray.h"
  #include "vtkDoubleArray.h"
  #include "vtkTimerLog.h"
  
  #include <list>
  #include <vector>
  #include <cassert>
  
  #ifndef NDEBUG
  // #define DEBUGTRACE
  // #define LOGTRACE
  #else
  #endif
  
  #ifdef DEBUGTRACE
  #define PRINT(x)  \
    {\
    cout<<this->Rank<<")"<<x<<endl;\
    }
  #define ALLPRINT(x)\
    for(int rank=0; rank<NumProcs; rank++){\
      Controller->Barrier();\
      if(rank==Rank) cout<<"("<<this->Rank<<")"<<x<<endl;\
    }
  #define Assert(a,msg)\
  {\
      if(!a) \
      {\
        cerr<<msg<<endl;\
        assert(false);\
      }\
  }
  
  #define AssertEq(a, b) \
  {\
      if(a!=b) \
      {\
        cerr<<a<<" != "<<b<<endl; \
        assert(false);\
      }\
  }
  
  #define AssertNe(a, b) \
  {\
      if(a==b) \
      {\
        cerr<<a<<" == "<<b<<endl; \
        assert(false);\
      }\
  }
  
  #define AssertGe(a, b) \
  {\
      if(a<b) \
      {\
        cerr<<a<<" < "<<b<<endl; \
        assert(false);\
      }\
  }
  
  #define AssertGt(a, b) \
  {\
      if(a<=b) \
      {\
        cerr<<a<<" < "<<b<<endl; \
        assert(false);\
      }\
  }
  //#define PRINT(id, x)
  #else
  #define PRINT(x)
  #define ALLPRINT(x)
  
  #define AssertEq(a,b)
  #define AssertGt(a,b)
  #define AssertGe(a,b)
  #define Assert(a,b)
  #define AssertNe(a,b)
  #endif
  
  namespace
  {
    inline int CNext(int i, int n)
    {
      return (i+1)%n;
    }
  
    class MyStream
    {
    public:
      MyStream(int BufferSize):Size(BufferSize)
      {
        this->Data = new char[Size];
        this->Head = Data;
      }
      ~MyStream()
      {
        delete [] this->Data;
      }
      int GetSize() { return Size;}
  
      template<class T>
      MyStream& operator << (T t)
      {
        unsigned int size =sizeof(T);
        char* value = reinterpret_cast<char*>(&t);
        for(unsigned int i=0; i<size;i++)
        {
          AssertGe(Size,this->Head-this->Data);
          *(this->Head++) = (*(value++));
        }
        return (*this);
      }
  
      template<class T>
      MyStream& operator >> (T& t)
      {
        unsigned int size =sizeof(T);
        AssertGe(Size,this->Head+size-this->Data);
        t = *(reinterpret_cast<T*>(this->Head));
        this->Head+=size;
        return (*this);
      }
  
      char * GetRawData(){ return this->Data;}
      int GetLength(){ return this->Head-this->Data;}
  
      void Reset()
      {
        this->Head = this->Data;
      }
  
    private:
      MyStream(const MyStream& ){};
      char* Data;
      char* Head;
      int Size;
    };
  
    inline void InitBB(double* Bounds)
    {
      Bounds[0] = DBL_MAX;
      Bounds[1] =-DBL_MAX;
      Bounds[2] = DBL_MAX;
      Bounds[3] =-DBL_MAX;
      Bounds[4] = DBL_MAX;
      Bounds[5] =-DBL_MAX;
    }
  
    inline bool InBB(const double* x, const double* bounds)
    {
      return       bounds[0] <= x[0] && x[0] <=bounds[1]
        &&  bounds[2] <= x[1] && x[1] <=bounds[3]
        &&  bounds[4] <= x[2] && x[2] <=bounds[5];
    }
  
    inline void UpdateBB(double* a, const double* b)
    {
      for(int i=0; i<=4; i+=2)
      {
        if(b[i]<a[i])
        {
          a[i] = b[i];
        }
      }
      for(int i=1; i<=5; i+=2)
      {
        if(b[i]>a[i])
        {
          a[i] = b[i];
        }
      }
    }
  }
  
  class PStreamTracerPoint: public vtkObject
  {
  public:
    vtkTypeMacro(PStreamTracerPoint, vtkObject);
  
    static PStreamTracerPoint *New();
  
  
    vtkGetMacro(Id,int)
    vtkGetVector3Macro(Seed,double)
    vtkGetVector3Macro(Normal,double)
    vtkGetMacro(Direction,int)
    vtkGetMacro(NumSteps,int)
    vtkGetMacro(Propagation,double)
    vtkGetMacro(Rank,int)
    vtkGetMacro(IntegrationTime,double);
  
    vtkSetMacro(Id,int)
    vtkSetMacro(Direction,int)
    vtkSetVector3Macro(Seed,double)
    vtkSetMacro(NumSteps,int)
    vtkSetMacro(Propagation,double)
    vtkSetMacro(Rank,int)
    vtkSetMacro(IntegrationTime,double);
  
    void Reseed(double* seed, double* normal, vtkPolyData* poly, int id,
                double propagation, double integrationTime)
    {
      memcpy( this->Seed, seed, 3*sizeof(double));
      memcpy( this->Normal, normal,3*sizeof(double));
  
      this->AllocateTail(poly->GetPointData());
      double* x = poly->GetPoints()->GetPoint(id);
      this->Tail->GetPoints()->SetPoint(0,x);
      this->Tail->GetPointData()->CopyData(poly->GetPointData(),id,0);
      this->Rank = -1; //someone else figure this out
      this->IntegrationTime = integrationTime;
      this->Propagation = propagation;
    }
  
    vtkPolyData* GetTail()
    {
      return this->Tail;
    }
  
    void CopyTail(PStreamTracerPoint* other)
    {
      if(other->Tail)
      {
        vtkPointData* pd = other->Tail->GetPointData();
        if(!this->Tail)
        {
          AllocateTail(pd);
        }
        this->Tail->GetPointData()->DeepCopy(pd);
      }
      else
      {
        Tail = nullptr;
      }
    }
  
  
    //allocate a one point vtkPolyData whose PointData setup matches pd
    void AllocateTail(vtkPointData* pd)
    {
      if(!this->Tail)
      {
        Tail = vtkSmartPointer<vtkPolyData>::New();
        vtkNew<vtkPoints> points;
        {
          points->SetNumberOfPoints(1);
        }
          Tail->SetPoints(points);
      }
  
  
      this->Tail->GetPointData()->CopyAllocate(pd);
    }
  
    virtual int GetSize()
    {
      int size(0);
      vtkPointData* data = this->GetTail()->GetPointData();
      for(int i=0; i<data->GetNumberOfArrays();i++)
      {
        size+= data->GetArray(i)->GetNumberOfComponents();
      }
      return size*sizeof(double) + sizeof(PStreamTracerPoint);
    }
  
  
    virtual void Read(MyStream& stream)
    {
      stream>> this->Id;
      stream>> this->Seed[0];
      stream>> this->Seed[1];
      stream>> this->Seed[2];
      stream>> this->Direction;
      stream>> this->NumSteps;
      stream>> this->Propagation;
      stream>> this->IntegrationTime;
  
      char hasTail(0);
      stream>> hasTail;
      if(hasTail)
      {
        double x[3];
        for(int i=0; i<3; i++)
        {
          stream>>x[i];
        }
        AssertNe(this->Tail,nullptr); //someone should have allocated it by prototype
        this->Tail->SetPoints( vtkSmartPointer<vtkPoints>::New());
        this->Tail->GetPoints()->InsertNextPoint(x);
  
        vtkPointData* pointData = this->Tail->GetPointData();
        for(int i=0; i<pointData->GetNumberOfArrays();i++)
        {
          int numComponents = pointData->GetArray(i)->GetNumberOfComponents();
          std::vector<double> xi(numComponents);
          for(int j=0; j<numComponents; j++)
          {
            double& xj(xi[j]);
            stream>>xj;
          }
          pointData->GetArray(i)->InsertNextTuple(&xi[0]);
        }
      }
      else
      {
        this->Tail = nullptr;
      }
    }
  
    virtual void Write(MyStream& stream)
    {
      stream <<this->Id
             <<this->Seed[0]
             <<this->Seed[1]
             <<this->Seed[2]
             <<this->Direction
             <<this->NumSteps
             <<this->Propagation
             <<this->IntegrationTime;
  
      stream<<(char)(this->Tail!=nullptr);
  
      if(this->Tail)
      {
        double* x = this->Tail->GetPoints()->GetPoint(0);
        for(int i=0; i<3; i++)
        {
          stream<<x[i];
        }
        vtkPointData* pData = this->Tail->GetPointData();
        int numArrays(pData->GetNumberOfArrays());
        for(int i=0; i<numArrays;i++)
        {
          vtkDataArray* arr = pData->GetArray(i);
          int numComponents = arr->GetNumberOfComponents();
          double* y = arr->GetTuple(0);
          for(int j=0; j<numComponents;j++)
            stream<<y[j];
        }
      }
    }
  
  private:
    int Id;
    double Seed[3];
    double Normal[3];
    int Direction;
    int NumSteps;
    double Propagation;
    vtkSmartPointer<vtkPolyData> Tail;
    int Rank;
    double IntegrationTime;
  
  protected:
    PStreamTracerPoint(): Id(-1)
                         ,Direction(0)
                         ,NumSteps(0)
                         ,Propagation(0)
                         ,Rank(-1)
                         ,IntegrationTime(0)
    {
      this->Seed[0] = this->Seed[1] = this->Seed[2] = -999;
    }
  };
  
  vtkStandardNewMacro(PStreamTracerPoint);
  
  
  class AMRPStreamTracerPoint: public PStreamTracerPoint
  {
  public:
    vtkTypeMacro(AMRPStreamTracerPoint, PStreamTracerPoint);
  
    static AMRPStreamTracerPoint *New();
    vtkSetMacro(Level,int)
    vtkGetMacro(Level,int)
  
    vtkSetMacro(GridId,int)
    vtkGetMacro(GridId,int)
  
    virtual int GetSize() override
    {
      return PStreamTracerPoint::GetSize()+2*sizeof(int);
    }
  
    virtual void Read(MyStream& stream) override
    {
      PStreamTracerPoint::Read(stream);
      stream>>Level;
      stream>>GridId;
    }
    virtual void Write(MyStream& stream) override
    {
      PStreamTracerPoint::Write(stream);
      stream<<Level<<GridId;
    }
  
  private:
    AMRPStreamTracerPoint():
      Level(-1),
      GridId(-1)
    {
    }
    int Level;
    int GridId;
  };
  
  typedef std::vector<vtkSmartPointer<PStreamTracerPoint> > PStreamTracerPointArray;
  
  vtkStandardNewMacro(AMRPStreamTracerPoint);
  
  
  class ProcessLocator: public vtkObject
  {
  public:
    vtkTypeMacro(ProcessLocator,vtkObject);
    static ProcessLocator *New();
    void Initialize(vtkCompositeDataSet* data)
    {
      this->Controller = vtkMultiProcessController::GetGlobalController();
      this->Rank = this->Controller->GetLocalProcessId();
      this->NumProcs = this->Controller->GetNumberOfProcesses();
      this->InitBoundingBoxes(this->NumProcs);
  
      double bb[6];
      InitBB(bb);
  
      if(data)
      {
        vtkCompositeDataIterator * iter = data->NewIterator();
        iter->InitTraversal();
        while ( !iter->IsDoneWithTraversal() )
        {
          vtkDataSet* dataSet = vtkDataSet::SafeDownCast(iter->GetCurrentDataObject());
          AssertNe(dataSet,nullptr);
          UpdateBB(bb,dataSet->GetBounds());
          iter->GoToNextItem();
        }
        iter->Delete();
      }
  
      PRINT(bb[0]<<" "<<bb[1]<<" "<<bb[2]<<" "<<bb[3]<<" "<<bb[4]<<" "<<bb[5]);
      this->Controller->AllGather(bb,&this->BoundingBoxes[0],6);
  
  #ifdef DEBUGTRACE
      cout<<"("<<Rank<<") BoundingBoxes: ";
      for(int i=0; i<NumProcs;i++)
      {
        double* box = this->GetBoundingBox(i);
        cout<<box[0]<<" "<<box[1]<<" "<<box[2]<<" "<<box[3]<<" "<<box[4]<<" "<<box[5]<<";  ";
      }
      cout<<endl;
  #endif
    }
  
    bool InCurrentProcess(double* p)
    {
      return InBB(p,GetBoundingBox(Rank));
    }
    int FindNextProcess(double* p)
    {
      for(int rank=CNext(this->Rank,this->NumProcs);
          rank!=Rank;rank=CNext(rank,this->NumProcs))
      {
        if(InBB(p,GetBoundingBox(rank)))
        {
          return rank;
        }
      }
      return -1;
    }
  private:
    ProcessLocator()
    {
    }
    vtkMultiProcessController* Controller;
    int Rank;
    int NumProcs;
  
    double* GetBoundingBox(int i)
    {
      return &BoundingBoxes[6*i];
    }
  
    void InitBoundingBoxes(int num)
    {
      for(int i=0; i<6*num; i++)
      {
        this->BoundingBoxes.push_back(0);
      }
    }
    std::vector<double> BoundingBoxes;
  };
  vtkStandardNewMacro(ProcessLocator);
  
  class AbstractPStreamTracerUtils: public vtkObject
  {
  public:
    vtkTypeMacro(AbstractPStreamTracerUtils,vtkObject);
  
    vtkGetMacro(VecName,char*);
    vtkGetMacro(VecType,int);
    vtkGetMacro(Input0,vtkDataSet*);
  
    virtual ProcessLocator* GetProcessLocator()
    {
      return nullptr;
    }
  
    vtkSmartPointer<PStreamTracerPoint> GetProto()
    {
      return this->Proto;
    }
  
    virtual void InitializeVelocityFunction(PStreamTracerPoint*, vtkAbstractInterpolatedVelocityField*)
    {
      return;
    }
  
    virtual bool PreparePoint(PStreamTracerPoint*, vtkAbstractInterpolatedVelocityField*)
    {
      return true;
    }
  
    void ComputeSeeds(vtkDataSet* source,PStreamTracerPointArray& out, int& maxId)
    {
      vtkDataArray* seeds;
      vtkIdList* seedIds;
      vtkIntArray* integrationDirections;
      this->Tracer->InitializeSeeds(seeds,
                              seedIds,
                              integrationDirections,
                              source);
  
      int numSeeds = seedIds->GetNumberOfIds();
      for (int i = 0; i < numSeeds; i ++ )
      {
        double seed[3];
        seeds->GetTuple(seedIds->GetId(i),seed);
        vtkSmartPointer<PStreamTracerPoint> point = NewPoint(i,seed,integrationDirections->GetValue(i));
        if(this->InBound(point))
        {
          out.push_back(point);
        }
      }
      if(seeds)
      {
        seeds->Delete();
      }
      if(seedIds)
      {
        seedIds->Delete();
      }
      if(integrationDirections)
      {
        integrationDirections->Delete();
      }
  
      maxId = numSeeds-1;
  
    }
  
    virtual void Initialize(vtkPStreamTracer* tracer)
    {
      this->Tracer = tracer;
      this->Controller = tracer->Controller;
      this->Rank = tracer->Rank;
      this->NumProcs = tracer->NumProcs;
      this->InputData = tracer->InputData;
      this->VecType = 0;
      this->VecName = nullptr;
      this->Input0 = 0;
      if(!tracer->EmptyData)
      {
        vtkCompositeDataIterator* iter = tracer->InputData->NewIterator();
        vtkSmartPointer<vtkCompositeDataIterator> iterP(iter);
        iter->Delete();
        iterP->GoToFirstItem();
        if(!iterP->IsDoneWithTraversal())
        {
          Input0 = vtkDataSet::SafeDownCast(iterP->GetCurrentDataObject());
          //iterP->GotoNextitem();
        }
        vtkDataArray* vectors = tracer->GetInputArrayToProcess(0,this->Input0,this->VecType);
        this->VecName = vectors->GetName();
      }
  
      if(!tracer->EmptyData)
      {
        this->CreatePrototype(this->Input0->GetPointData(),VecType,VecName);
      }
    }
  
  protected:
    AbstractPStreamTracerUtils():Tracer(nullptr), Controller(nullptr), VecName(nullptr),Input0(nullptr),InputData(nullptr)
    {
    }
  
    virtual vtkSmartPointer<PStreamTracerPoint> NewPoint(int id, double* x, int dir)=0;
    virtual bool InBound(PStreamTracerPoint* p) = 0;
  
    void CreatePrototype(vtkPointData* pointData, int fieldType, const char* vecName)
    {
      this->Proto  = NewPoint(-1,nullptr,-1);
  
      vtkNew<vtkPointData> protoPD;
      protoPD->InterpolateAllocate(pointData,1);
      vtkSmartPointer<vtkDoubleArray> time = vtkSmartPointer<vtkDoubleArray>::New();
      time->SetName("IntegrationTime");
      protoPD->AddArray(time);
  
      if(fieldType==vtkDataObject::FIELD_ASSOCIATION_CELLS)
      {
        vtkSmartPointer<vtkDoubleArray> velocityVectors = vtkSmartPointer<vtkDoubleArray>::New();
        velocityVectors->SetName(vecName);
        velocityVectors->SetNumberOfComponents(3);
        protoPD->AddArray(velocityVectors);
      }
  
      if(Tracer->GetComputeVorticity())
      {
        vtkSmartPointer<vtkDoubleArray> vorticity = vtkSmartPointer<vtkDoubleArray>::New();
        vorticity->SetName("Vorticity");
        vorticity->SetNumberOfComponents(3);
        protoPD->AddArray(vorticity);
  
        vtkSmartPointer<vtkDoubleArray> rotation = vtkSmartPointer<vtkDoubleArray>::New();
        rotation->SetName("Rotation");
        protoPD->AddArray(rotation);
  
        vtkSmartPointer<vtkDoubleArray> angularVel  = vtkSmartPointer<vtkDoubleArray>::New();
        angularVel->SetName("AngularVelocity");
        protoPD->AddArray(angularVel);
      }
  
      if(Tracer->GenerateNormalsInIntegrate)
      {
        PRINT("Generate normals prototype");
        vtkSmartPointer<vtkDoubleArray> normals = vtkSmartPointer<vtkDoubleArray>::New();
        normals->SetName("Normals");
        normals->SetNumberOfComponents(3);
        protoPD->AddArray(normals);
      }
      AssertEq(this->Proto->GetTail(),nullptr);
      this->Proto->AllocateTail(protoPD);
    }
  
    vtkPStreamTracer* Tracer;
    vtkMultiProcessController* Controller;
    vtkSmartPointer<PStreamTracerPoint> Proto;
    int VecType;
    char *VecName;
    vtkDataSet* Input0;
    vtkCompositeDataSet* InputData;
  
    int Rank;
    int NumProcs;
  
  };
  
  class PStreamTracerUtils: public AbstractPStreamTracerUtils
  {
  public:
    vtkTypeMacro(PStreamTracerUtils,AbstractPStreamTracerUtils);
  
    static PStreamTracerUtils *New();
  
    PStreamTracerUtils()
    {
      this->Locator = nullptr;
    }
  
    virtual void Initialize(vtkPStreamTracer* tracer) override
    {
      this->Superclass::Initialize(tracer);
      this->Locator = vtkSmartPointer<ProcessLocator>::New();
      this->Locator->Initialize(tracer->InputData);
    }
  
    virtual ProcessLocator* GetProcessLocator() override
    {
      return this->Locator;
    }
  
  
    virtual bool InBound(PStreamTracerPoint*) override
    {
      return true;
    }
  
    virtual vtkSmartPointer<PStreamTracerPoint> NewPoint(int id, double* x, int dir) override
    {
      vtkSmartPointer<PStreamTracerPoint>  p = vtkSmartPointer<PStreamTracerPoint>::New();
      p->SetId(id);
      if(x)
      {
        p->SetSeed(x);
      }
      p->SetDirection(dir);
      return p;
    }
  private:
    vtkSmartPointer<ProcessLocator> Locator;
  };
  
  vtkStandardNewMacro(PStreamTracerUtils);
  
  class AMRPStreamTracerUtils: public AbstractPStreamTracerUtils
  {
  public:
    vtkTypeMacro(AMRPStreamTracerUtils,AbstractPStreamTracerUtils);
    static AMRPStreamTracerUtils *New();
    vtkSetMacro(AMR,vtkOverlappingAMR*);
  
    virtual void InitializeVelocityFunction(PStreamTracerPoint* point, vtkAbstractInterpolatedVelocityField* func) override
    {
      AMRPStreamTracerPoint* amrPoint = AMRPStreamTracerPoint::SafeDownCast(point);
      assert(amrPoint);
      vtkAMRInterpolatedVelocityField* amrFunc = vtkAMRInterpolatedVelocityField::SafeDownCast(func);
      assert(amrFunc);
      if(amrPoint->GetLevel()>=0)
      {
        amrFunc->SetLastDataSet(amrPoint->GetLevel(),amrPoint->GetGridId());
  #ifdef DEBUGTRACE
        vtkUniformGrid* grid = this->AMR->GetDataSet(amrPoint->GetLevel(),amrPoint->GetGridId());
        if(!grid || !InBB(amrPoint->GetSeed(),grid->GetBounds()))
        {
          PRINT("WARNING: Bad AMR Point "<<(grid)<<" "<<amrPoint->GetSeed()[0]<<" "<<amrPoint->GetSeed()[1]<<" "<<amrPoint->GetSeed()[2]<<
                " "<<amrPoint->GetLevel()<<" "<<amrPoint->GetGridId());
        }
  #endif
      }
    }
  
  
    virtual bool PreparePoint(PStreamTracerPoint* point, vtkAbstractInterpolatedVelocityField* func) override
    {
      AMRPStreamTracerPoint* amrPoint = AMRPStreamTracerPoint::SafeDownCast(point);
      vtkAMRInterpolatedVelocityField* amrFunc = vtkAMRInterpolatedVelocityField::SafeDownCast(func);
      unsigned int level, id;
      if(amrFunc->GetLastDataSetLocation(level,id))
      {
        amrPoint->SetLevel(level);
        amrPoint->SetId(id);
        int blockIndex = this->AMR->GetCompositeIndex(level,id);
        amrPoint->SetRank(this->BlockProcess[blockIndex]);
        return true;
      }
      else
      {
        PRINT("Invalid AMR : "<<point->GetSeed()[0]<<" "<<point->GetSeed()[1]<<" "<<point->GetSeed()[2]<<" "<<"Probably out of bound");
        amrPoint->SetLevel(-1);
        amrPoint->SetGridId(-1);
        amrPoint->SetRank(-1);
        return false;
      }
    }
  
    //this assume that p's AMR information has been set correctly
    //it makes no attempt to look for it
    virtual bool InBound(PStreamTracerPoint* p) override
    {
      AMRPStreamTracerPoint* amrp = AMRPStreamTracerPoint::SafeDownCast(p);
      if(amrp->GetLevel()<0)
      {
        return false;
      }
      AssertNe(amrp,nullptr);
      vtkUniformGrid* grid= this->AMR->GetDataSet(amrp->GetLevel(),amrp->GetGridId());
      return grid!=nullptr;
    }
  
    virtual vtkSmartPointer<PStreamTracerPoint> NewPoint(int id, double* x, int dir) override
    {
  
      vtkSmartPointer<AMRPStreamTracerPoint>  amrp  = vtkSmartPointer<AMRPStreamTracerPoint>::New();
      vtkSmartPointer<PStreamTracerPoint>  p  = amrp;
      p->SetId(id);
      if(x)
      {
        p->SetSeed(x);
      }
      p->SetDirection(dir);
  
      if(x)
      {
        unsigned int level, gridId;
        if(vtkAMRInterpolatedVelocityField::FindGrid(x,this->AMR,level,gridId))
        {
          amrp->SetLevel((int)level);
          amrp->SetGridId((int)gridId);
          int blockIndex =this->AMR->GetCompositeIndex(level,gridId);
          int process  =this->BlockProcess[blockIndex];
          AssertGe(process,0);
          amrp->SetRank(process);
        }
        else
        {
        }
      }
  
      return p;
    }
    virtual void Initialize(vtkPStreamTracer* tracer) override
    {
      this->Superclass::Initialize(tracer);
      AssertNe(this->InputData,nullptr);
      this->AMR = vtkOverlappingAMR::SafeDownCast(this->InputData);
  
      vtkParallelAMRUtilities::DistributeProcessInformation(this->AMR, this->Controller, BlockProcess);
      this->AMR->GenerateParentChildInformation();
    }
  
  protected:
    AMRPStreamTracerUtils()
    {
      this->AMR = nullptr;
    }
    vtkOverlappingAMR* AMR;
  
    std::vector<int> BlockProcess; //stores block->process information
  };
  vtkStandardNewMacro(AMRPStreamTracerUtils);
  
  
  
  
  namespace
  {
    inline double normvec3(double* x, double* y)
    {
      return sqrt( ( x[0] - y[0] ) * ( x[0] - y[0] )
                   + ( x[1] - y[1] ) * ( x[1] - y[1] )
                   + ( x[2] - y[2] ) * ( x[2] - y[2] ) );
    }
  
    inline vtkIdType LastPointIndex(vtkPolyData* pathPoly)
    {
      vtkCellArray* pathCells = pathPoly->GetLines();
      AssertGt(pathCells->GetNumberOfCells(),0);
      vtkIdType* path(0);
      vtkIdType nPoints(0);
      pathCells->InitTraversal();
      pathCells->GetNextCell(nPoints,path);
      int lastPointIndex = path[nPoints-1];
      return lastPointIndex;
  
    }
  
  #ifdef DEBUGTRACE
    inline double ComputeLength(vtkIdList* poly, vtkPoints* pts)
    {
      int n = poly->GetNumberOfIds();
      if(n==0) return 0;
  
      double s(0);
      double p[3];
      pts->GetPoint(poly->GetId(0),p);
      for(int j=1; j<n;j++)
      {
        int pIndex = poly->GetId(j);
        double q[3];
        pts->GetPoint(pIndex,q);
        s+= sqrt( vtkMath::Distance2BetweenPoints(p,q));
        memcpy(p,q,3*sizeof(double));
      }
      return s;
    }
  
    inline void PrintNames(ostream& out, vtkPointData* a)
    {
      for(int i=0; i<a->GetNumberOfArrays();i++)
      {
        out<< a->GetArray(i)->GetName()<<" ";
      }
      out<<endl;
    }
  
    inline bool SameShape(vtkPointData* a, vtkPointData* b)
    {
      if (!a || !b)
      {
        return false;
      }
  
      if(a->GetNumberOfArrays()!=b->GetNumberOfArrays())
      {
        PrintNames(cerr,a);
        PrintNames(cerr,b);
        return false;
      }
  
      int numArrays(a->GetNumberOfArrays());
      for(int i=0; i<numArrays;i++)
      {
        if(a->GetArray(i)->GetNumberOfComponents()!=b->GetArray(i)->GetNumberOfComponents())
        {
          return false;
        }
      }
  
      return true;
    }
  #endif
  
    class MessageBuffer
    {
    public:
      MessageBuffer(int size)
      {
        this->Stream = new MyStream(size);
      }
  
      ~MessageBuffer()
      {
        delete this->Stream;
      }
  
      vtkMPICommunicator::Request& GetRequest()
      {
        return this->Request;
      }
      MyStream& GetStream()
      {
        return *(this->Stream);
      }
  
    private:
      vtkMPICommunicator::Request Request;
      MyStream* Stream;
      MessageBuffer(const MessageBuffer& ){};
    };
  
    typedef MyStream MessageStream;
  
    class Task: public vtkObject
    {
    public:
  
      static Task *New();
  
      int GetId()
      {
        return this->Point->GetId();
      }
      vtkGetMacro(TraceExtended,bool)
      vtkGetMacro(TraceTerminated,bool)
      vtkSetMacro(TraceExtended,bool)
      vtkSetMacro(TraceTerminated,bool)
  
      PStreamTracerPoint* GetPoint()
      {
        return this->Point;
      }
      void IncHop()
      {
        this->NumHops++;
      }
    private:
      vtkSmartPointer<PStreamTracerPoint> Point;
  
      int NumPeeks;
      int NumHops;
      bool TraceTerminated;
      bool TraceExtended;
  
      Task()
        :NumPeeks(0)
        ,NumHops(0)
        ,TraceTerminated(false)
        ,TraceExtended(false)
      {
  
      }
      friend class TaskManager;
      friend MessageStream& operator<<(MessageStream& stream, const Task& task);
    };
    vtkStandardNewMacro(Task);
  
    MessageStream& operator<<(MessageStream& stream, const Task& t)
    {
      t.Point->Write(stream);
      stream<<t.NumPeeks;
      stream<<t.NumHops;
      return stream;
    }
  
    //Description:
    //Manages the communication of traces between processes
    class TaskManager
    {
    public:
      enum Message
      {
        NewTask,
        NoMoreTasks,
        TaskFinished
      };
  
      TaskManager(ProcessLocator* locator, PStreamTracerPoint* proto):
        Locator(locator),  Proto(proto)
      {
        this->Controller = vtkMPIController::SafeDownCast(vtkMultiProcessController::GetGlobalController());
        AssertNe(this->Controller,nullptr);
        this->NumProcs = this->Controller->GetNumberOfProcesses();
        this->Rank = this->Controller->GetLocalProcessId();
  
        int prototypeSize = Proto==nullptr? 0: Proto->GetSize();
        this->MessageSize = prototypeSize+sizeof(Task);
        this->ReceiveBuffer = nullptr;
  
        this->NumSends =0;
        this->Timer = vtkSmartPointer<vtkTimerLog>::New();
        this->ReceiveTime = 0;
      }
  
      void Initialize(bool hasData, const PStreamTracerPointArray& seeds, int MaxId)
      {
        AssertGe(MaxId,0);
        int numSeeds  = static_cast<int>(seeds.size());
        this->HasData.clear();
        {
          for(int i=0; i<NumProcs;i++)
            this->HasData.push_back(0);
  
          std::vector<int> hasDataIn(NumProcs);
          for(int i=0; i<NumProcs;i++)
          {
            hasDataIn[i] = i==Rank? hasData: 0;
          }
          this->Controller->AllReduce(&hasDataIn[0],&this->HasData[0],NumProcs,vtkCommunicator::MAX_OP);
        }
  
        for(int i=0; i<NumProcs;i++)
        {
          if(this->HasData[i])
          {
            this->Leader=i;
            break;
          }
        }
  
        std::vector<int> processMap0(MaxId+1);
        for(int i=0; i<MaxId+1; i++)
        {
          processMap0[i] = -1;
        }
        for (int i = 0; i < numSeeds; i ++ )
        {
          int rank  = seeds[i]->GetRank();
          int id = seeds[i]->GetId();
          if(rank<0 && this->Locator)
          {
            rank = this->Locator->InCurrentProcess(seeds[i]->GetSeed())? this->Rank : -1;
          }
          processMap0[id] = rank;
        }
  
        std::vector<int> processMap(MaxId+1);
        this->Controller->AllReduce(&processMap0[0], &processMap[0],MaxId+1,vtkCommunicator::MAX_OP);
  
        int totalNumTasks(0);
        for (int id = 0; id <=MaxId; id++)
        {
          if(processMap[id]>=0)
          {
            totalNumTasks++;
          }
        }
        this->TotalNumTasks = Rank==this->Leader? totalNumTasks: INT_MAX; //only the master process knows how many are left
  
        for (int i = 0; i < numSeeds; i++ )
        {
          int id = seeds[i]->GetId();
          if(processMap[id]==Rank)
          {
            vtkNew<Task> task;
            task->Point = seeds[i];
            NTasks.push_back(task);
          }
        }
        ALLPRINT(NTasks.size()<<" initial seeds out of "<<totalNumTasks);
  
      }
  
      Task* NextTask()
      {
        if(!this->HasData[Rank])
        {
          return nullptr;
        }
  
        //---------------------------------------------------------
        // Send messages
        //---------------------------------------------------------
  
        while(!this->PTasks.empty())
        {
          vtkSmartPointer<Task> task  = PTasks.back();
          PTasks.pop_back();
  
          if(task->GetTraceTerminated())
          {
            //send to the master process
            this->Send(TaskFinished,this->Leader,task);
          }
          else
          {
            if(!task->GetTraceExtended())
            {
              //increment the peak
              task->NumPeeks++;
              PRINT("Skip "<<task->GetId()<<" with "<<task->NumPeeks<<" Peeks");
            }
            else
            {
              task->NumPeeks=1;
            }
            int nextProcess = -1;
            if(task->NumPeeks<this->NumProcs)
            {
              nextProcess = NextProcess(task);
              if(nextProcess>=0)
              {
                task->IncHop();
                //send it to the next guy
                this->Send(NewTask,NextProcess(task),task);
              }
            }
  
            if(nextProcess<0)
            {
              this->Send(TaskFinished,this->Leader,task); //no one can do it, norminally finished
              PRINT("Bail on "<<task->GetId());
            }
          }
        }
  
        //---------------------------------------------------------
        // Receive messages
        //---------------------------------------------------------
  
        do
        {
          this->Receive(this->TotalNumTasks!=0 && this->Msgs.empty() && NTasks.empty()); //wait if there is nothing to do
          while(!this->Msgs.empty())
          {
            Message msg = this->Msgs.back();
            this->Msgs.pop_back();
            switch(msg)
            {
              case NewTask: break;
              case TaskFinished:
                AssertEq(Rank,this->Leader);
                this->TotalNumTasks--;
                PRINT(TotalNumTasks<<" tasks left");
                break;
              case NoMoreTasks:
                AssertNe(Rank,this->Leader);
                this->TotalNumTasks=0;
                break;
              default: assert(false);
            }
          }
        }while(this->TotalNumTasks!=0 && NTasks.empty());
  
        vtkSmartPointer<Task> nextTask;
        if(NTasks.empty())
        {
          AssertEq( this->TotalNumTasks,0);
          if(this->Rank==this->Leader)   //let everyeone know
          {
            for(int i=(this->Rank+1)%NumProcs; i!=this->Rank;i=(i+1)%NumProcs)
            {
              if(this->HasData[i])
              {
                this->Send(NoMoreTasks,i,0);
              }
            }
          }
        }
        else
        {
          nextTask = this->NTasks.back();
          this->NTasks.pop_back();
          this->PTasks.push_back(nextTask);
        }
  
        return nextTask;
      }
  
      ~TaskManager()
      {
        for( BufferList::iterator itr=SendBuffers.begin();itr!=SendBuffers.end();++itr)
        {
          MessageBuffer* buf = *itr;
          AssertNe(buf->GetRequest().Test(),0);
          delete buf;
        }
        if(this->ReceiveBuffer)
        {
          this->ReceiveBuffer->GetRequest().Cancel();
          delete ReceiveBuffer;
        }
      }
    private:
      ProcessLocator* Locator;
      vtkSmartPointer<PStreamTracerPoint> Proto;
      vtkMPIController* Controller;
      std::vector<vtkSmartPointer<Task> > NTasks;
      std::vector<vtkSmartPointer<Task> > PTasks;
      std::vector<Message> Msgs;
      int NumProcs;
      int Rank;
      int TotalNumTasks;
      int MessageSize;
      std::vector<int> HasData;
      int Leader;
      typedef  std::list<MessageBuffer*> BufferList;
      BufferList SendBuffers;
      MessageBuffer* ReceiveBuffer;
  
      void Send(int msg, int rank, Task* task)
      {
        if(task && (msg==TaskFinished))
        {
          PRINT("Done in "<<task->Point->GetNumSteps()<<" steps "<<task->NumHops<<" hops");
        }
        if(rank==this->Rank)
        {
          switch(msg)
          {
            case TaskFinished:
              this->TotalNumTasks--;
              PRINT(TotalNumTasks<<" tasks left");
              break;
            default:
              PRINT("Unhandled message "<<msg);
              assert(false);
          }
        }
        else
        {
          MessageBuffer& buf = this->NewSendBuffer();
          MessageStream& outStream(buf.GetStream());
  
          outStream<<msg<<this->Rank;
          AssertNe(this->Rank,rank);
  
          if(task)
          {
            outStream<<(*task);
          }
  
          AssertGe(this->MessageSize,outStream.GetLength());
          this->Controller->NoBlockSend(outStream.GetRawData(),outStream.GetLength(),rank,561,buf.GetRequest());
  
          NumSends++;
          if(task)
          {
            PRINT("Send "<<msg<<"; task "<<task->GetId());//<<" "<<task->Seed[0]<<" "<<task->Seed[1]<<" "<<task->Seed[2]<<" to "<<rank);
          }
          else
          {
            PRINT("Send "<<msg);
          }
        }
  
      }
      int NextProcess(Task* task)
      {
        PStreamTracerPoint* p = task->GetPoint();
        int rank = p->GetRank();
        if(rank>=0)
        {
          return rank;
        }
  
        if(this->Locator)
        {
          rank = this->Locator->FindNextProcess(p->GetSeed());
        }
        AssertNe(rank,Rank);
        return rank;
      }
  
      vtkSmartPointer<Task> NewTaskInstance()
      {
        vtkSmartPointer<Task> task = vtkSmartPointer<Task>::New();
  
        task->Point  = this->Proto->NewInstance();
        task->Point->CopyTail(this->Proto);
        task->Point->Delete();
        return task;
      }
  
      void Read(MessageStream& stream, Task& task)
      {
        task.Point->Read(stream);
        stream>> task.NumPeeks;
        stream>> task.NumHops;
      }
  
      MessageBuffer& NewSendBuffer()
      {
        //remove all empty buffers
        BufferList::iterator itr = SendBuffers.begin();
        while(itr!=SendBuffers.end())
        {
          MessageBuffer* buf(*itr);
          BufferList::iterator next = itr;
          ++next;
          if(buf->GetRequest().Test())
          {
            delete buf;
            SendBuffers.erase(itr);
          }
          itr = next;
        }
  
        MessageBuffer* buf = new MessageBuffer(this->MessageSize);
        SendBuffers.push_back(buf);
        return *buf;
      }
  
      void Receive(bool wait = false)
      {
        int msg=-1;
        int sender(0);
  
  #ifdef DEBUGTRACE
        this->StartTimer();
  #endif
        if(ReceiveBuffer && wait)
        {
          ReceiveBuffer->GetRequest().Wait();
        }
  
        if(ReceiveBuffer && ReceiveBuffer->GetRequest().Test())
        {
          MyStream& inStream(ReceiveBuffer->GetStream());
          inStream >>msg >> sender;
          this->Msgs.push_back( (Message)msg);
          if(msg==NewTask)
          {
            PRINT("Received message "<<msg<<" from "<<sender)
  
            vtkSmartPointer<Task> task = this->NewTaskInstance();
            this->Read(inStream,*task);
            PRINT("Received task "<<task->GetId());//<<" "<<task->Seed[0]<<" "<<task->Seed[1]<<" "<<task->Seed[2]);
            this->NTasks.push_back(task);
          }
          delete ReceiveBuffer;  ReceiveBuffer = nullptr;
        }
        if(ReceiveBuffer==nullptr)
        {
          ReceiveBuffer = new MessageBuffer(this->MessageSize);
          MyStream& inStream(ReceiveBuffer->GetStream());
          this->Controller->NoBlockReceive(inStream.GetRawData(),
                                           inStream.GetSize(),
                                           vtkMultiProcessController::ANY_SOURCE,561,
                                           ReceiveBuffer->GetRequest());
        }
  
  #ifdef DEBUGTRACE
        double time = this->StopTimer();
        if(msg>=0)
        {
          this->ReceiveTime+=time;
        }
  #endif
  
      }
  
  
      int NumSends;
      double ReceiveTime;
      vtkSmartPointer<vtkTimerLog> Timer;
    };
  
  };
  
  vtkCxxSetObjectMacro(vtkPStreamTracer, Controller, vtkMultiProcessController);
  vtkCxxSetObjectMacro(vtkPStreamTracer,
                       Interpolator,
                       vtkAbstractInterpolatedVelocityField);
  vtkStandardNewMacro( vtkPStreamTracer );
  
  vtkPStreamTracer::vtkPStreamTracer()
  {
    this->Controller = vtkMultiProcessController::GetGlobalController();
    if (this->Controller)
    {
      this->Controller->Register(this);
    }
    this->Interpolator = 0;
    this->GenerateNormalsInIntegrate = 0;
  
    this->EmptyData = 0;
  }
  
  vtkPStreamTracer::~vtkPStreamTracer()
  {
    if (this->Controller)
    {
      this->Controller->UnRegister(this);
      this->Controller = 0;
    }
    this->SetInterpolator(0);
  }
  
  int vtkPStreamTracer::RequestUpdateExtent(
    vtkInformation *vtkNotUsed(request),
    vtkInformationVector **inputVector,
    vtkInformationVector *outputVector)
  {
    vtkInformation *outInfo = outputVector->GetInformationObject(0);
    int piece =
      outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER());
    int numPieces =
      outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES());
    int ghostLevel =
      outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS());
  
    int numInputs = this->GetNumberOfInputConnections(0);
    for (int idx = 0; idx < numInputs; ++idx)
    {
      vtkInformation *info = inputVector[0]->GetInformationObject(idx);
      if (info)
      {
        info->Set(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER(),
                  piece);
        info->Set(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES(),
                  numPieces);
        info->Set(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS(),
                  ghostLevel);
      }
    }
  
  
    vtkInformation *sourceInfo = inputVector[1]->GetInformationObject(0);
    if (sourceInfo)
    {
      sourceInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER(),
                      0);
      sourceInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES(),
                      1);
      sourceInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS(),
                      ghostLevel);
    }
  
    return 1;
  }
  
  int vtkPStreamTracer::RequestData(
    vtkInformation *request,
    vtkInformationVector **inputVector,
    vtkInformationVector *outputVector)
  {
    if (!vtkMPIController::SafeDownCast(this->Controller) || this->Controller->GetNumberOfProcesses() == 1)
    {
      this->GenerateNormalsInIntegrate = 1;
      int result = vtkStreamTracer::RequestData(request,inputVector,outputVector);
      this->GenerateNormalsInIntegrate = 0;
      return result;
    }
  
    this->Rank = this->Controller->GetLocalProcessId();
    this->NumProcs = this->Controller->GetNumberOfProcesses();
  
  
    vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
    vtkInformation *outInfo = outputVector->GetInformationObject(0);
    if (!this->SetupOutput(inInfo, outInfo))
    {
      return 0;
    }
  
    vtkInformation *sourceInfo = inputVector[1]->GetInformationObject(0);
    vtkDataSet *source = 0;
    if (sourceInfo)
    {
      source = vtkDataSet::SafeDownCast(
        sourceInfo->Get(vtkDataObject::DATA_OBJECT()));
    }
    vtkPolyData* output = vtkPolyData::SafeDownCast(
      outInfo->Get(vtkDataObject::DATA_OBJECT()));
  
    // init 'func' with nullptr such that we can check it later to determine
    // if we need to deallocate 'func' in case CheckInputs() fails (note
    // that a process may be assigned no any dataset when the number of
    // processes is greater than that of the blocks)
    vtkAbstractInterpolatedVelocityField * func = nullptr;
    int maxCellSize = 0;
    if (this->CheckInputs(func, &maxCellSize) != VTK_OK)
    {
      vtkDebugMacro("No appropriate inputs have been found..");
      this->EmptyData = 1;
      PRINT("Has Empty Data")
  
      // the if-statement below is a MUST since 'func' may be still nullptr
      // when this->InputData is nullptr ---- no any data has been assigned
      // to this process
      if ( func )
      {
        func->Delete();
        func = nullptr;
      }
    }
    else
    {
      func->SetCaching(0);
      this->SetInterpolator(func);
      func->Delete();
    }
  
    if(vtkOverlappingAMR::SafeDownCast(this->InputData))
    {
      this->Utils = vtkSmartPointer<AMRPStreamTracerUtils>::New();
    }
    else
    {
      this->Utils = vtkSmartPointer<PStreamTracerUtils>::New();
    }
    this->Utils->Initialize(this);
    ALLPRINT("Vec Name: "<<this->Utils->GetVecName());
    typedef std::vector< vtkSmartPointer<vtkPolyData> > traceOutputsType;
    traceOutputsType traceOutputs;
  
    TaskManager taskManager(this->Utils->GetProcessLocator(),this->Utils->GetProto());
    PStreamTracerPointArray seedPoints;
  
  
    int maxId;
    this->Utils->ComputeSeeds(source,seedPoints,maxId);
    taskManager.Initialize(this->EmptyData==0,seedPoints,maxId);
  
    Task* task(0);
    std::vector<int> traceIds;
    int iterations = 0;
    while( (task = taskManager.NextTask()))
    {
      iterations++;
      PStreamTracerPoint* point = task->GetPoint();
  
      vtkSmartPointer<vtkPolyData> traceOut;
      this->Trace(this->Utils->GetInput0(),
            this->Utils->GetVecType(),
            this->Utils->GetVecName(),
            point,
            traceOut,
            func,
            maxCellSize);
  
      task->SetTraceExtended(traceOut->GetNumberOfPoints()>0);
  
      if(task->GetTraceExtended() && task->GetPoint()->GetTail())
      {
        // if we got this streamline from another process then this
        // process is responsible for filling in the gap over
        // the subdomain boundary
        this->Prepend(traceOut,task->GetPoint()->GetTail());
      }
  
      int resTerm=vtkStreamTracer::OUT_OF_DOMAIN;
      vtkIntArray* resTermArray = vtkArrayDownCast<vtkIntArray>(
        traceOut->GetCellData()->GetArray("ReasonForTermination"));
      if (resTermArray)
      {
        resTerm = resTermArray->GetValue(0);
      }
  
      //construct a new seed from the last point
      task->SetTraceTerminated(this->Controller->GetNumberOfProcesses()==1
                               || resTerm != vtkStreamTracer::OUT_OF_DOMAIN
                               || point->GetPropagation() > this->MaximumPropagation
                               || point->GetNumSteps()    >=  this->MaximumNumberOfSteps);
      if(task->GetTraceExtended() && !task->GetTraceTerminated())
      {
        task->SetTraceTerminated(!this->TraceOneStep(traceOut,func,point)); //we don't know where to go, just terminate it
      }
      if(!task->GetTraceTerminated())
      {
        task->SetTraceTerminated(!this->Utils->PreparePoint(point,func));
      }
  
      traceIds.push_back(task->GetId());
      traceOutputs.push_back(traceOut);
    }
  
    this->Controller->Barrier();
  
  #ifdef LOGTRACE
    double receiveTime = taskManager.ComputeReceiveTime();
    if(this->Rank==0)
    {
      PRINT("Total receive time: "<<receiveTime)
    }
    this->Controller->Barrier();
  #endif
  
    PRINT("Done");
  
    // The parallel integration adds all streamlines to traceOutputs
    // container. We append them all together here.
    vtkNew<vtkAppendPolyData> append;
    for (traceOutputsType::iterator it = traceOutputs.begin();
         it != traceOutputs.end(); ++it)
    {
      vtkPolyData* inp = it->GetPointer();
      if ( inp->GetNumberOfCells() > 0 )
      {
        append->AddInputData(inp);
      }
    }
    if (append->GetNumberOfInputConnections(0) > 0)
    {
      append->Update();
      vtkPolyData* appoutput = append->GetOutput();
      output->CopyStructure(appoutput);
      output->GetPointData()->PassData(appoutput->GetPointData());
      output->GetCellData()->PassData(appoutput->GetCellData());
    }
  
  
    this->InputData->UnRegister(this);
  
  #ifdef DEBUGTRACE
    int maxSeeds(maxId+1);
    std::vector<double> lengths(maxSeeds);
    for(int i=0; i<maxSeeds;i++)
    {
      lengths[i] = 0;
    }
  
    AssertEq(traceOutputs.size(),traceIds.size());
    for(unsigned int i=0; i<traceOutputs.size();i++)
    {
      vtkPolyData* poly = traceOutputs[i];
      int id = traceIds[i];
      double length(0);
      vtkCellArray* lines = poly->GetLines();
      if(lines)
      {
        lines->InitTraversal();
        vtkNew<vtkIdList> trace;
        lines->GetNextCell(trace);
        length= ComputeLength(trace,poly->GetPoints());
      }
      lengths[id] += length;
    }
    std::vector<double> totalLengths(maxSeeds);
    this->Controller->AllReduce(&lengths[0],&totalLengths[0],maxSeeds,vtkCommunicator::SUM_OP);
  
    int numNonZeros(0);
    double totalLength(0);
    for(int i=0; i<maxSeeds;i++)
    {
      totalLength+=totalLengths[i];
      if(totalLengths[i]>0)
      {
        numNonZeros++;
      }
    }
  
    if(this->Rank==0)
    {
      PRINT("Summary: "<<maxSeeds<<" seeds,"<<numNonZeros<<" traces"<<" total length "<<totalLength);
    }
  
  #endif
    PRINT("Done in "<<iterations<<" iterations");
  
    traceOutputs.erase(traceOutputs.begin(), traceOutputs.end());
    return 1;
  }
  
  void vtkPStreamTracer::PrintSelf(ostream& os, vtkIndent indent)
  {
    this->Superclass::PrintSelf(os,indent);
    os << indent << "Controller: " << this->Controller << endl;
  }
  
  
  void vtkPStreamTracer::Trace( vtkDataSet *input,
                                int vecType,
                                const char* vecName,
                                PStreamTracerPoint* point,
                                vtkSmartPointer<vtkPolyData>& traceOut,
                                vtkAbstractInterpolatedVelocityField* func,
                                int maxCellSize)
  {
    double* seedSource = point->GetSeed();
    int direction = point->GetDirection();
  
    this->Utils->InitializeVelocityFunction(point, func);
  
    double lastPoint[3];
    vtkSmartPointer<vtkFloatArray> seeds = vtkSmartPointer<vtkFloatArray>::New();
    seeds->SetNumberOfComponents(3);
    seeds->InsertNextTuple(seedSource);
  
    vtkNew<vtkIdList> seedIds;
    seedIds->InsertNextId(0);
  
    vtkNew<vtkIntArray> integrationDirections;
    integrationDirections->InsertNextValue(direction);
    traceOut =  vtkSmartPointer<vtkPolyData>::New();
  
    double propagation = point->GetPropagation();
    vtkIdType numSteps = point->GetNumSteps();
    double integrationTime = point->GetIntegrationTime();
  
    vtkStreamTracer::Integrate(input->GetPointData(),
                      traceOut,
                      seeds,
                      seedIds,
                      integrationDirections,
                      lastPoint,
                      func,
                      maxCellSize,
                      vecType,
                      vecName,
                      propagation,
                      numSteps,
                      integrationTime);
    AssertGe(propagation, point->GetPropagation());
    AssertGe(numSteps, point->GetNumSteps());
  
    point->SetPropagation(propagation);
    point->SetNumSteps(numSteps);
    point->SetIntegrationTime(integrationTime);
  
    if(this->GenerateNormalsInIntegrate)
    {
      this->GenerateNormals(traceOut, point->GetNormal(), vecName);
    }
  
    if(traceOut->GetNumberOfPoints()>0)
    {
      if (traceOut->GetLines()->GetNumberOfCells()==0)
      {
        PRINT( "Fix Single Point Path")
        AssertEq(traceOut->GetNumberOfPoints(),1); //fix it
        vtkNew<vtkCellArray> newCells;
        vtkIdType cell;
        cell = 0;
        newCells->InsertNextCell(1,&cell);
        traceOut->SetLines(newCells);
  
        // Don't forget to add the ReasonForTermination cell array.
        vtkNew<vtkIntArray> retVals;
        retVals->SetName("ReasonForTermination");
        retVals->SetNumberOfTuples(1);
        retVals->SetValue(0, vtkStreamTracer::OUT_OF_DOMAIN);
        traceOut->GetCellData()->AddArray(retVals);
      }
  
      vtkNew<vtkIntArray> ids;
      ids->SetName("SeedIds");
      ids->SetNumberOfTuples(1);
      ids->SetValue(0, point->GetId());
      traceOut->GetCellData()->AddArray(ids);
    }
    Assert(SameShape(traceOut->GetPointData(),this->Utils->GetProto()->GetTail()->GetPointData()),"trace data does not match prototype");
  
  }
  bool vtkPStreamTracer::TraceOneStep(vtkPolyData* traceOut,  vtkAbstractInterpolatedVelocityField* func, PStreamTracerPoint* point)
  {
    double outPoint[3],outNormal[3];
  
    vtkIdType lastPointIndex = LastPointIndex(traceOut);
    double lastPoint[3];
    // Continue the integration a bit further to obtain a point
    // outside. The main integration step can not always be used
    // for this, specially if the integration is not 2nd order.
    traceOut->GetPoint(lastPointIndex, lastPoint);
  
    vtkInitialValueProblemSolver* ivp = this->Integrator;
    ivp->Register(this);
  
    vtkNew<vtkRungeKutta2> tmpSolver;
    this->SetIntegrator(tmpSolver);
  
    memcpy(outPoint,lastPoint,sizeof(double)*3);
  
    double timeStepTaken = this->SimpleIntegrate(0, outPoint, this->LastUsedStepSize, func);
    PRINT("Simple Integrate from :"<<lastPoint[0]<<" "<<lastPoint[1]<<" "<<lastPoint[2]<<" to "<<outPoint[0]<<" "<<outPoint[1]<<" "<<outPoint[2]);
    double d = sqrt(vtkMath::Distance2BetweenPoints(lastPoint,outPoint));
  
    this->SetIntegrator(ivp);
    ivp->UnRegister(this);
  
    vtkDataArray* normals = traceOut->GetPointData()->GetArray("Normals");
    if (normals)
    {
      normals->GetTuple(lastPointIndex, outNormal);
    }
  
    bool res = d>0;
    if(res)
    {
      Assert(SameShape(traceOut->GetPointData(),this->Utils->GetProto()->GetTail()->GetPointData()),"Point data mismatch");
      point->Reseed(outPoint,outNormal,traceOut,lastPointIndex,
                    point->GetPropagation()+d, point->GetIntegrationTime()+timeStepTaken);
      AssertEq(point->GetTail()->GetPointData()->GetNumberOfTuples(),1);
    }
  
    return res;
  }
  
  void vtkPStreamTracer::Prepend(vtkPolyData* pathPoly, vtkPolyData* headPoly)
  {
    vtkCellArray* pathCells = pathPoly->GetLines();
    AssertEq(pathCells->GetNumberOfCells(),1);
    AssertEq(headPoly->GetNumberOfPoints(),1);
  
    double* newPoint = headPoly->GetPoint(0);
    AssertEq(headPoly->GetPointData()->GetNumberOfArrays(),pathPoly->GetPointData()->GetNumberOfArrays());
  
    vtkIdType* path(0);
    vtkIdType nPoints(0);
    pathCells->InitTraversal();
    pathCells->GetNextCell(nPoints,path);
    AssertNe(path,nullptr);
    AssertEq(nPoints,pathPoly->GetNumberOfPoints());
  
    vtkIdType newPointId = pathPoly->GetPoints()->InsertNextPoint(newPoint);
  
    vtkPointData* headData = headPoly->GetPointData();
    vtkPointData* pathData = pathPoly->GetPointData();
    Assert(SameShape(headData,pathData),"Prepend failure");
    int numArrays(headData->GetNumberOfArrays());
    for(int i=0; i<numArrays;i++)
    {
      pathData->CopyTuple(headData->GetAbstractArray(i),pathData->GetAbstractArray(i),0,newPointId);
    }
  
    PRINT("Prepend Point "<<newPointId<<" "<<newPoint[0]<<" "<<newPoint[1]<<" "<<newPoint[2]);
    vtkNew<vtkIdList> newPath;
    newPath->InsertNextId(newPointId);
    for(int i=0; i<nPoints;i++)
    {
      newPath->InsertNextId(path[i]);
    }
  
    pathCells->Reset();
    pathCells->InsertNextCell(newPath);
    AssertEq(pathCells->GetNumberOfCells(),1);
    vtkIdType newNumPoints(0);
    pathCells->GetNextCell(newNumPoints,path);
    AssertEq(newNumPoints,nPoints+1);
    AssertEq(newNumPoints,pathPoly->GetNumberOfPoints());
  }
  