vtkParticleTracerBase.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
// SPDX-License-Identifier: BSD-3-Clause
#ifndef vtkParticleTracerBase_h
#define vtkParticleTracerBase_h
 
#include "vtkFiltersFlowPathsModule.h" // For export macro
#include "vtkPolyDataAlgorithm.h"
#include "vtkSmartPointer.h"      // For vtkSmartPointer
#include "vtkTemporalAlgorithm.h" // For vtkTemporalAlgorithm
#include "vtkWeakPointer.h"       // For vtkWeakPointer
 
#include <list>          // STL Header
#include <mutex>         // STL Header
#include <numeric>       // STL Header
#include <unordered_map> // STL Header
#include <vector>        // STL Header
 
#ifndef __VTK_WRAP__
#define vtkPolyDataAlgorithm vtkTemporalAlgorithm<vtkPolyDataAlgorithm>
#endif
 
VTK_ABI_NAMESPACE_BEGIN
class vtkAbstractParticleWriter;
class vtkCompositeDataSet;
class vtkDataArray;
class vtkDataSet;
class vtkDoubleArray;
class vtkFloatArray;
class vtkGenericCell;
class vtkInitialValueProblemSolver;
class vtkIntArray;
class vtkMultiBlockDataSet;
class vtkMultiProcessController;
class vtkPartitionedDataSet;
class vtkPointData;
class vtkPoints;
class vtkPolyData;
class vtkSignedCharArray;
class vtkTemporalInterpolatedVelocityField;
VTK_ABI_NAMESPACE_END
 
namespace vtkParticleTracerBaseNamespace
{
VTK_ABI_NAMESPACE_BEGIN
struct Position_t
{
  double x[4];
};
using Position = struct Position_t;
 
struct ParticleInformation_t
{
  // These are used during iteration
  Position CurrentPosition;
  int CachedDataSetId[2];
  vtkIdType CachedCellId[2];
  int LocationState;
  // These are computed scalars we might display
  int SourceID;
  int TimeStepAge; // amount of time steps the particle has advanced
  int InjectedPointId;
  int InjectedStepId; // time step the particle was injected
  double SimulationTime;
  // These are useful to track for debugging etc
  int ErrorCode;
  float age;
  // these are needed across time steps to compute vorticity
  float rotation;
  float angularVel;
  float time;
  float speed;
  // once the particle is added, PointId is valid and is the tuple location
  // in ProtoPD.
  vtkIdType PointId;
 
  double velocity[3];
};
using ParticleInformation = struct ParticleInformation_t;
 
typedef std::vector<ParticleInformation> ParticleVector;
typedef ParticleVector::iterator ParticleIterator;
typedef std::list<ParticleInformation> ParticleDataList;
typedef ParticleDataList::iterator ParticleListIterator;
struct ParticleTracerFunctor;
VTK_ABI_NAMESPACE_END
}
 
VTK_ABI_NAMESPACE_BEGIN
class VTKFILTERSFLOWPATHS_EXPORT vtkParticleTracerBase : public vtkPolyDataAlgorithm
{
public:
  friend struct vtkParticleTracerBaseNamespace::ParticleTracerFunctor;
  enum Solvers
  {
    RUNGE_KUTTA2,
    RUNGE_KUTTA4,
    RUNGE_KUTTA45,
    NONE,
    UNKNOWN
  };
 
  vtkTypeMacro(vtkParticleTracerBase, vtkPolyDataAlgorithm);
  void PrintParticleHistories();
 
#ifndef __VTK_WRAP__
#undef vtkPolyDataAlgorithm
#endif
  void PrintSelf(ostream& os, vtkIndent indent) override;
 
#if defined(__VTK_WRAP__) || defined(__WRAP_GCCXML)
  vtkCreateWrappedTemporalAlgorithmInterface();
#endif
 
 
  virtual void SetController(vtkMultiProcessController*);
  virtual vtkMultiProcessController* GetController();
 
 
  vtkGetMacro(ComputeVorticity, bool);
  void SetComputeVorticity(bool);
 
 
  vtkGetMacro(TerminalSpeed, double);
  void SetTerminalSpeed(double);
 
 
  vtkGetMacro(RotationScale, double);
  void SetRotationScale(double);
 
 
  vtkSetMacro(IgnorePipelineTime, vtkTypeBool);
  vtkGetMacro(IgnorePipelineTime, vtkTypeBool);
  vtkBooleanMacro(IgnorePipelineTime, vtkTypeBool);
 
 
  vtkGetMacro(ForceReinjectionEveryNSteps, int);
  void SetForceReinjectionEveryNSteps(int);
 
  void SetIntegrator(vtkInitialValueProblemSolver*);
  vtkGetObjectMacro(Integrator, vtkInitialValueProblemSolver);
 
  void SetIntegratorType(int type);
  int GetIntegratorType();
 
 
  vtkSetMacro(StaticSeeds, vtkTypeBool);
  vtkGetMacro(StaticSeeds, vtkTypeBool);
 
  enum MeshOverTimeTypes
  {
    DIFFERENT = 0,
    STATIC = 1,
    LINEAR_TRANSFORMATION = 2,
    SAME_TOPOLOGY = 3
  };
 
  /*
   * Set/Get the type of variance of the mesh over time.
   *
   * DIFFERENT = 0,
   * STATIC = 1,
   * LINEAR_TRANSFORMATION = 2
   * SAME_TOPOLOGY = 3
   */
  virtual void SetMeshOverTime(int meshOverTime);
  virtual int GetMeshOverTimeMinValue() { return DIFFERENT; }
  virtual int GetMeshOverTimeMaxValue() { return SAME_TOPOLOGY; }
  void SetMeshOverTimeToDifferent() { this->SetMeshOverTime(DIFFERENT); }
  void SetMeshOverTimeToStatic() { this->SetMeshOverTime(STATIC); }
  void SetMeshOverTimeToLinearTransformation() { this->SetMeshOverTime(LINEAR_TRANSFORMATION); }
  void SetMeshOverTimeToSameTopology() { this->SetMeshOverTime(SAME_TOPOLOGY); }
  vtkGetMacro(MeshOverTime, int);
 
  enum
  {
    INTERPOLATOR_WITH_DATASET_POINT_LOCATOR,
    INTERPOLATOR_WITH_CELL_LOCATOR
  };
 
  void SetInterpolatorType(int interpolatorType);
 
  void SetInterpolatorTypeToDataSetPointLocator();
 
  void SetInterpolatorTypeToCellLocator();
 
 
  virtual void SetParticleWriter(vtkAbstractParticleWriter* pw);
  vtkGetObjectMacro(ParticleWriter, vtkAbstractParticleWriter);
 
 
  vtkSetFilePathMacro(ParticleFileName);
  vtkGetFilePathMacro(ParticleFileName);
 
 
  vtkSetMacro(EnableParticleWriting, vtkTypeBool);
  vtkGetMacro(EnableParticleWriting, vtkTypeBool);
  vtkBooleanMacro(EnableParticleWriting, vtkTypeBool);
 
 
  void AddSourceConnection(vtkAlgorithmOutput* input);
  void RemoveAllSources();
 
 
  vtkGetMacro(ForceSerialExecution, bool);
  vtkSetMacro(ForceSerialExecution, bool);
  vtkBooleanMacro(ForceSerialExecution, bool);
protected:
 
  vtkSmartPointer<vtkPointData> ProtoPD;
  vtkParticleTracerBaseNamespace::ParticleDataList ParticleHistories;
  vtkTypeBool IgnorePipelineTime; // whether to use the pipeline time for termination
 
  // Control execution as serial or threaded
  bool ForceSerialExecution;
 
  void EnqueueParticleToAnotherProcess(vtkParticleTracerBaseNamespace::ParticleInformation&);
 
  vtkParticleTracerBase();
  ~vtkParticleTracerBase() override;
 
  //
  // Make sure the pipeline knows what type we expect as input
  //
  int FillInputPortInformation(int port, vtkInformation* info) override;
 
  int Initialize(vtkInformation* request, vtkInformationVector** inputVector,
    vtkInformationVector* outputVector) override;
 
  int Execute(vtkInformation* request, vtkInformationVector** inputVector,
    vtkInformationVector* outputVector) override;
 
  int Finalize(vtkInformation*, vtkInformationVector**, vtkInformationVector*) override;
 
  virtual std::vector<vtkDataSet*> GetSeedSources(vtkInformationVector* inputVector);
 
  // Initialization of input (vector-field) geometry
  int InitializeInterpolator();
 
  vtkSmartPointer<vtkDataSet> Seeds;
 
  void TestParticles(vtkParticleTracerBaseNamespace::ParticleVector& candidates,
    vtkParticleTracerBaseNamespace::ParticleVector& passed);
 
  void TestParticles(
    vtkParticleTracerBaseNamespace::ParticleVector& candidates, std::vector<int>& passed);
 
  virtual void AssignSeedsToProcessors(double time, vtkDataSet* source,
    vtkParticleTracerBaseNamespace::ParticleVector& localSeedPoints);
 
  bool SendReceiveParticles(std::vector<vtkIdType>&);
 
  void UpdateParticleList(vtkParticleTracerBaseNamespace::ParticleVector& candidates);
 
  virtual bool UpdateParticleListFromOtherProcesses();
 
  void IntegrateParticle(vtkParticleTracerBaseNamespace::ParticleListIterator& it,
    double currentTime, double targetTime, vtkInitialValueProblemSolver* integrator,
    vtkTemporalInterpolatedVelocityField* interpolator, vtkDoubleArray* cellVectors,
    std::atomic<vtkIdType>& particleCount, std::mutex& eraseMutex, bool sequential);
 
  bool ComputeDomainExitLocation(
    double pos[4], double p2[4], double intersection[4], vtkGenericCell* cell);
 
  //
  // Scalar arrays that are generated as each particle is updated
  //
  void CreateProtoPD(vtkDataObject* input);
 
  vtkFloatArray* GetParticleAge(vtkPointData*);
  vtkIntArray* GetParticleIds(vtkPointData*);
  vtkSignedCharArray* GetParticleSourceIds(vtkPointData*);
  vtkIntArray* GetInjectedPointIds(vtkPointData*);
  vtkIntArray* GetInjectedStepIds(vtkPointData*);
  vtkIntArray* GetErrorCodeArr(vtkPointData*);
  vtkFloatArray* GetParticleVorticity(vtkPointData*);
  vtkFloatArray* GetParticleRotation(vtkPointData*);
  vtkFloatArray* GetParticleAngularVel(vtkPointData*);
 
  // utility function we use to test if a point is inside any of our local datasets
  bool InsideBounds(double point[]);
 
  void CalculateVorticity(
    vtkGenericCell* cell, double pcoords[3], vtkDoubleArray* cellVectors, double vorticity[3]);
 
  //------------------------------------------------------
 
  void SetParticle(vtkParticleTracerBaseNamespace::ParticleInformation& info,
    vtkTemporalInterpolatedVelocityField* interpolator, vtkDoubleArray* cellVectors);
 
  std::unordered_map<vtkIdType, int> InjectedPointIdToProcessId;
 
 
  virtual bool IsPointDataValid(vtkDataObject* input);
  bool IsPointDataValid(vtkCompositeDataSet* input, std::vector<std::string>& arrayNames);
  void GetPointDataArrayNames(vtkDataSet* input, std::vector<std::string>& names);
 
  vtkGetMacro(ReinjectionCounter, int);
 
  void ResizeArrays(vtkIdType numTuples);
 
  virtual void InitializeExtraPointDataArrays(vtkPointData* vtkNotUsed(outputPD)) {}
 
  virtual void SetToExtraPointDataArrays(
    vtkIdType, vtkParticleTracerBaseNamespace::ParticleInformation&)
  {
  }
 
  vtkTemporalInterpolatedVelocityField* GetInterpolator();
 
  virtual void AddRestartSeeds(vtkInformationVector** /*inputVector*/) {}
 
  bool RetryWithPush(vtkParticleTracerBaseNamespace::ParticleInformation& info, double* point1,
    double delT, int subSteps, vtkTemporalInterpolatedVelocityField* interpolator);
 
  bool SetTerminationTimeNoModify(double t);
 
  double CachedTimeStep;
 
  // Parameters of tracing
  vtkInitialValueProblemSolver* Integrator;
  double IntegrationStep;
  double MaximumError;
  bool ComputeVorticity;
  double RotationScale;
  double TerminalSpeed;
 
  // A counter to keep track of how many times we reinjected
  int ReinjectionCounter;
 
  // Important for Caching of Cells/Ids/Weights etc
  vtkTypeBool AllFixedGeometry;
  int MeshOverTime;
  vtkTypeBool StaticSeeds;
 
  // Innjection parameters
  int ForceReinjectionEveryNSteps;
  vtkTimeStamp ParticleInjectionTime;
  bool HasCache;
 
  // Particle writing to disk
  vtkAbstractParticleWriter* ParticleWriter;
  char* ParticleFileName;
  vtkTypeBool EnableParticleWriting;
 
  vtkParticleTracerBaseNamespace::ParticleVector LocalSeeds;
 
  // The velocity interpolator
  vtkSmartPointer<vtkTemporalInterpolatedVelocityField> Interpolator;
 
  // MPI controller needed when running in parallel
  vtkSmartPointer<vtkMultiProcessController> Controller;
 
  vtkParticleTracerBaseNamespace::ParticleVector MPISendList;
 
  std::unordered_map<vtkIdType, vtkParticleTracerBaseNamespace::ParticleInformation> MPIRecvList;
 
  // Cache bounds info for each dataset we will use repeatedly
  struct bounds_t
  {
    double b[6];
  };
  using bounds = struct bounds_t;
  std::vector<bounds> CachedBounds[2];
 
  // variables used by Execute() to produce output
 
  vtkSmartPointer<vtkDataSet> DataReferenceT[2];
 
  vtkNew<vtkPoints> OutputCoordinates;
  vtkNew<vtkIdTypeArray> ParticleCellsConnectivity;
 
  vtkNew<vtkFloatArray> ParticleAge;
  vtkNew<vtkIntArray> ParticleIds;
  vtkNew<vtkSignedCharArray> ParticleSourceIds;
  vtkNew<vtkIdTypeArray> InjectedPointIds;
  vtkNew<vtkIntArray> InjectedStepIds;
  vtkNew<vtkIntArray> ErrorCodeArray;
  vtkNew<vtkFloatArray> ParticleVorticity;
  vtkNew<vtkFloatArray> ParticleRotation;
  vtkNew<vtkFloatArray> ParticleAngularVel;
  vtkNew<vtkPointData> OutputPointData;
 
  // temp array
  vtkNew<vtkDoubleArray> CellVectors;
 
  vtkParticleTracerBase(const vtkParticleTracerBase&) = delete;
  void operator=(const vtkParticleTracerBase&) = delete;
  vtkTimeStamp ExecuteTime;
 
  unsigned int NumberOfParticles();
 
  friend class ParticlePathFilterInternal;
  friend class StreaklineFilterInternal;
 
  static const double Epsilon;
 
private:
  // Internal method to initialize CachedData[1] using the provided input
  void InitializeNextCachedData(vtkDataObject* input);
 
  // Data for time step CurrentTimeStep-1 and CurrentTimeStep
  vtkSmartPointer<vtkPartitionedDataSet> CachedData[2];
};
 
VTK_ABI_NAMESPACE_END
#endif