vtkLegacyParticleTracerBase.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 vtkLegacyParticleTracerBase_h
#define vtkLegacyParticleTracerBase_h

#include "vtkDeprecation.h"                            // For VTK_DEPRECATED_IN_9_2_0
#include "vtkPVVTKExtensionsFiltersGeneralMPIModule.h" // For export macro
#include "vtkPolyDataAlgorithm.h"
#include "vtkSmartPointer.h" // For vtkSmartPointer

#include <list>   // STL Header
#include <mutex>  // STL Header
#include <vector> // STL Header

VTK_ABI_NAMESPACE_BEGIN
class vtkAbstractParticleWriter;
class vtkCellArray;
class vtkCompositeDataSet;
class vtkDataArray;
class vtkDataSet;
class vtkDoubleArray;
class vtkFloatArray;
class vtkGenericCell;
class vtkInitialValueProblemSolver;
class vtkIntArray;
class vtkMultiBlockDataSet;
class vtkMultiProcessController;
class vtkPointData;
class vtkPoints;
class vtkPolyData;
class vtkSignedCharArray;
class vtkTemporalInterpolatedVelocityField;
VTK_ABI_NAMESPACE_END

namespace vtkLegacyParticleTracerBaseNamespace
{
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
  int UniqueParticleId;
  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;
  // if PointId is negative then in parallel this particle was just
  // received and we need to get the tuple value from vtkPParticleTracerBase::Tail.
  vtkIdType TailPointId;
};
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 VTKPVVTKEXTENSIONSFILTERSGENERALMPI_EXPORT vtkLegacyParticleTracerBase
  : public vtkPolyDataAlgorithm
{
public:
  friend struct vtkLegacyParticleTracerBaseNamespace::ParticleTracerFunctor;
  enum Solvers
  {
    RUNGE_KUTTA2,
    RUNGE_KUTTA4,
    RUNGE_KUTTA45,
    NONE,
    UNKNOWN
  };

  vtkTypeMacro(vtkLegacyParticleTracerBase, vtkPolyDataAlgorithm);
  void PrintSelf(ostream& os, vtkIndent indent) override;
  void PrintParticleHistories();


  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 SetTerminationTime(double t);
  vtkGetMacro(TerminationTime, double);

  void SetIntegrator(vtkInitialValueProblemSolver*);
  vtkGetObjectMacro(Integrator, vtkInitialValueProblemSolver);

  void SetIntegratorType(int type);
  int GetIntegratorType();


  vtkGetMacro(StartTime, double);
  void SetStartTime(double t);


  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);


  VTK_DEPRECATED_IN_9_2_0("Use SetMeshOverTime instead")
  virtual void SetStaticMesh(vtkTypeBool staticMesh)
  {
    this->SetMeshOverTime(staticMesh ? STATIC : DIFFERENT);
  }
  VTK_DEPRECATED_IN_9_2_0("Use GetMeshOverTime instead")
  virtual vtkTypeBool GetStaticMesh() { return this->MeshOverTime == STATIC; }

  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);


  vtkSetMacro(DisableResetCache, vtkTypeBool);
  vtkGetMacro(DisableResetCache, vtkTypeBool);
  vtkBooleanMacro(DisableResetCache, vtkTypeBool);


  void AddSourceConnection(vtkAlgorithmOutput* input);
  void RemoveAllSources();


  vtkGetMacro(ForceSerialExecution, bool);
  vtkSetMacro(ForceSerialExecution, bool);
  vtkBooleanMacro(ForceSerialExecution, bool);
protected:
  vtkSmartPointer<vtkPolyData> Output; // managed by child classes

  vtkSmartPointer<vtkPointData> ProtoPD;
  vtkIdType UniqueIdCounter; // global Id counter used to give particles a stamp
  vtkLegacyParticleTracerBaseNamespace::ParticleDataList ParticleHistories;
  vtkSmartPointer<vtkPointData> ParticlePointData; // the current particle point data consistent
                                                   // with particle history
  // Everything related to time
  vtkTypeBool IgnorePipelineTime; // whether to use the pipeline time for termination
  vtkTypeBool DisableResetCache;  // whether to enable ResetCache() method

  // Control execution as serial or threaded
  bool ForceSerialExecution;

  vtkLegacyParticleTracerBase();
  ~vtkLegacyParticleTracerBase() override;

  //
  // Make sure the pipeline knows what type we expect as input
  //
  int FillInputPortInformation(int port, vtkInformation* info) override;

  //
  // The usual suspects
  //
  vtkTypeBool ProcessRequest(vtkInformation* request, vtkInformationVector** inputVector,
    vtkInformationVector* outputVector) override;

  //
  // Store any information we need in the output and fetch what we can
  // from the input
  //
  int RequestInformation(vtkInformation* request, vtkInformationVector** inputVector,
    vtkInformationVector* outputVector) override;

  //
  // Compute input time steps given the output step
  //
  int RequestUpdateExtent(vtkInformation* request, vtkInformationVector** inputVector,
    vtkInformationVector* outputVector) override;

  //
  // what the pipeline calls for each time step
  //
  int RequestData(vtkInformation* request, vtkInformationVector** inputVector,
    vtkInformationVector* outputVector) override;

  //
  // these routines are internally called to actually generate the output
  //
  virtual int ProcessInput(vtkInformationVector** inputVector);

  // This is the main part of the algorithm:
  //  * move all the particles one step
  //  * Reinject particles (by adding them to this->ParticleHistories)
  //    either at the beginning or at the end of each step (modulo
  //    this->ForceReinjectionEveryNSteps)
  //  * Output a polydata representing the moved particles
  // Note that if the starting and the ending time coincide, the polydata is still valid.
  virtual vtkPolyData* Execute(vtkInformationVector** inputVector);

  // the RequestData will call these methods in turn
  virtual void Initialize() {}                        // the first iteration
  virtual int OutputParticles(vtkPolyData* poly) = 0; // every iteration
  virtual void Finalize() {}                          // the last iteration

  virtual std::vector<vtkDataSet*> GetSeedSources(vtkInformationVector* inputVector, int timeStep);

  // Initialization of input (vector-field) geometry
  int InitializeInterpolator();
  int UpdateDataCache(vtkDataObject* td);

  void TestParticles(vtkLegacyParticleTracerBaseNamespace::ParticleVector& candidates,
    vtkLegacyParticleTracerBaseNamespace::ParticleVector& passed, int& count);

  void TestParticles(
    vtkLegacyParticleTracerBaseNamespace::ParticleVector& candidates, std::vector<int>& passed);

  virtual void AssignSeedsToProcessors(double time, vtkDataSet* source, int sourceID, int ptId,
    vtkLegacyParticleTracerBaseNamespace::ParticleVector& localSeedPoints, int& localAssignedCount);

  virtual void AssignUniqueIds(
    vtkLegacyParticleTracerBaseNamespace::ParticleVector& localSeedPoints);

  void UpdateParticleList(vtkLegacyParticleTracerBaseNamespace::ParticleVector& candidates);

  virtual bool UpdateParticleListFromOtherProcesses() { return false; }

  void IntegrateParticle(vtkLegacyParticleTracerBaseNamespace::ParticleListIterator& it,
    double currentTime, double targetTime, vtkInitialValueProblemSolver* integrator,
    vtkTemporalInterpolatedVelocityField* interpolator, vtkDoubleArray* cellVectors,
    std::atomic<vtkIdType>& particleCount, std::mutex& eraseMutex, bool sequential);

  // if the particle is added to send list, then returns value is 1,
  // if it is kept on this process after a retry return value is 0
  virtual bool SendParticleToAnotherProcess(
    vtkLegacyParticleTracerBaseNamespace::ParticleInformation&,
    vtkLegacyParticleTracerBaseNamespace::ParticleInformation&, vtkPointData*)
  {
    return true;
  }

  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]);

  //------------------------------------------------------

  double GetCacheDataTime(int i);
  double GetCacheDataTime();

  virtual void ResetCache();
  void SetParticle(vtkLegacyParticleTracerBaseNamespace::ParticleInformation& info,
    double* velocity, vtkTemporalInterpolatedVelocityField* interpolator, vtkIdType particleId,
    vtkDoubleArray* cellVectors);


  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);
  vtkGetMacro(CurrentTimeValue, double);

  void ResizeArrays(vtkIdType numTuples);

  virtual void InitializeExtraPointDataArrays(vtkPointData* vtkNotUsed(outputPD)) {}

  virtual void SetToExtraPointDataArrays(
    vtkIdType, vtkLegacyParticleTracerBaseNamespace::ParticleInformation&)
  {
  }

  vtkTemporalInterpolatedVelocityField* GetInterpolator();

  virtual void AddRestartSeeds(vtkInformationVector** /*inputVector*/) {}

  virtual void RenameGhostArray(vtkPointData* pd);

private:
  bool RetryWithPush(vtkLegacyParticleTracerBaseNamespace::ParticleInformation& info,
    double* point1, double delT, int subSteps, vtkTemporalInterpolatedVelocityField* interpolator);

  bool SetTerminationTimeNoModify(double t);

  // 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;

  std::vector<double> InputTimeValues;
  double StartTime;
  double TerminationTime;
  double CurrentTimeValue;

  int StartTimeStep; // InputTimeValues[StartTimeStep] <= StartTime <=
                     // InputTimeValues[StartTimeStep+1]
  int CurrentTimeStep;
  int TerminationTimeStep; // computed from start time
  bool FirstIteration;

  // Innjection parameters
  int ForceReinjectionEveryNSteps;
  vtkTimeStamp ParticleInjectionTime;
  bool HasCache;

  // Particle writing to disk
  vtkAbstractParticleWriter* ParticleWriter;
  char* ParticleFileName;
  vtkTypeBool EnableParticleWriting;

  // The main lists which are held during operation- between time step updates
  vtkLegacyParticleTracerBaseNamespace::ParticleVector LocalSeeds;

  // The velocity interpolator
  vtkSmartPointer<vtkTemporalInterpolatedVelocityField> Interpolator;

  // Data for time step CurrentTimeStep-1 and CurrentTimeStep
  vtkSmartPointer<vtkMultiBlockDataSet> CachedData[2];

  // 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];

  vtkSmartPointer<vtkPoints> OutputCoordinates;
  vtkSmartPointer<vtkIdTypeArray> ParticleCellsConnectivity;
  vtkSmartPointer<vtkIdTypeArray> ParticleCellsOffsets;
  vtkSmartPointer<vtkCellArray> ParticleCells;

  vtkSmartPointer<vtkFloatArray> ParticleAge;
  vtkSmartPointer<vtkIntArray> ParticleIds;
  vtkSmartPointer<vtkSignedCharArray> ParticleSourceIds;
  vtkSmartPointer<vtkIntArray> InjectedPointIds;
  vtkSmartPointer<vtkIntArray> InjectedStepIds;
  vtkSmartPointer<vtkIntArray> ErrorCodeArray;
  vtkSmartPointer<vtkFloatArray> ParticleVorticity;
  vtkSmartPointer<vtkFloatArray> ParticleRotation;
  vtkSmartPointer<vtkFloatArray> ParticleAngularVel;
  vtkSmartPointer<vtkPointData> OutputPointData;

  // temp array
  vtkSmartPointer<vtkDoubleArray> CellVectors;

  vtkLegacyParticleTracerBase(const vtkLegacyParticleTracerBase&) = delete;
  void operator=(const vtkLegacyParticleTracerBase&) = delete;
  vtkTimeStamp ExecuteTime;

  unsigned int NumberOfParticles();

  friend class vtkLegacyPParticlePathFilter;

  static const double Epsilon;
};

VTK_ABI_NAMESPACE_END
#endif