1 /*========================================================================= 2 3 Program: Visualization Toolkit 4 Module: vtkEvenlySpacedStreamlines2D.h 5 6 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen 7 All rights reserved. 8 See Copyright.txt or http://www.kitware.com/Copyright.htm for details. 9 10 This software is distributed WITHOUT ANY WARRANTY; without even 11 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 12 PURPOSE. See the above copyright notice for more information. 13 14 =========================================================================*/ 15 /** 16 * @class vtkEvenlySpacedStreamlines2D 17 * @brief Evenly spaced streamline generator for 2D. 18 * 19 * vtkEvenlySpacedStreamlines2D is a filter that integrates a 2D 20 * vector field to generate evenly-spaced streamlines. 21 * 22 * We implement 23 * the algorithm described in: 24 * Jobard, Bruno, and Wilfrid Lefer. "Creating evenly-spaced 25 * streamlines of arbitrary density." Visualization in Scientific 26 * Computing '97. Springer Vienna, 1997. 43-55. 27 * The loop detection is described in: 28 * Liu, Zhanping, Robert Moorhead, and Joe Groner. 29 * "An advanced evenly-spaced streamline placement algorithm." 30 * IEEE Transactions on Visualization and Computer Graphics 12.5 (2006): 965-972. 31 * 32 * The integration is performed using a specified integrator, 33 * by default Runge-Kutta2. 34 * 35 * vtkEvenlySpacedStreamlines2D produces polylines as the output, with 36 * each cell (i.e., polyline) representing a streamline. The attribute 37 * values associated with each streamline are stored in the cell data, 38 * whereas those associated with streamline-points are stored in the 39 * point data. 40 * 41 * vtkEvenlySpacedStreamlines2D integrates streamlines both forward 42 * and backward. The integration for a streamline terminates upon 43 * exiting the flow field domain, or if the particle speed is reduced 44 * to a value less than a specified terminal speed, if the current 45 * streamline gets too close to other streamlines 46 * (vtkStreamTracer::FIXED_REASONS_FOR_TERMINATION_COUNT + 1) or if 47 * the streamline forms a loop 48 * (vtkStreamTracer::FIXED_REASONS_FOR_TERMINATION_COUNT). The 49 * specific reason for the termination is stored in a cell array named 50 * ReasonForTermination. 51 * 52 * Note that normalized vectors are adopted in streamline integration, 53 * which achieves high numerical accuracy/smoothness of flow lines that is 54 * particularly guaranteed for Runge-Kutta45 with adaptive step size and 55 * error control). In support of this feature, the underlying step size is 56 * ALWAYS in arc length unit (LENGTH_UNIT) while the 'real' time interval 57 * (virtual for steady flows) that a particle actually takes to trave in a 58 * single step is obtained by dividing the arc length by the LOCAL speed. 59 * The overall elapsed time (i.e., the life span) of the particle is the 60 * sum of those individual step-wise time intervals. 61 * 62 * The quality of streamline integration can be controlled by setting 63 * the initial integration step (InitialIntegrationStep), particularly 64 * for Runge-Kutta2 and Runge-Kutta4 (with a fixed step size). We do 65 * not support Runge-Kutta45 (with an adaptive step size and error 66 * control) because a requirement of the algorithm is that sample 67 * points along a streamline be evenly spaced. These steps are in 68 * either LENGTH_UNIT or CELL_LENGTH_UNIT. 69 * 70 * The integration time, vorticity, rotation and angular velocity are stored 71 * in point data arrays named "IntegrationTime", "Vorticity", "Rotation" and 72 * "AngularVelocity", respectively (vorticity, rotation and angular velocity 73 * are computed only when ComputeVorticity is on). All point data attributes 74 * in the source dataset are interpolated on the new streamline points. 75 * 76 * vtkEvenlySpacedStreamlines2D supports integration through any type of 2D dataset. 77 * 78 * The starting point, or the so-called 'seed', of the first streamline is set 79 * by setting StartPosition 80 * 81 * @sa 82 * vtkStreamTracer vtkRibbonFilter vtkRuledSurfaceFilter vtkInitialValueProblemSolver 83 * vtkRungeKutta2 vtkRungeKutta4 vtkRungeKutta45 vtkTemporalStreamTracer 84 * vtkAbstractInterpolatedVelocityField vtkInterpolatedVelocityField 85 * vtkCellLocatorInterpolatedVelocityField 86 * 87 */ 88 89 #ifndef vtkEvenlySpacedStreamlines2D_h 90 #define vtkEvenlySpacedStreamlines2D_h 91 92 #include "vtkFiltersFlowPathsModule.h" // For export macro 93 #include "vtkPolyDataAlgorithm.h" 94 95 #include <array> 96 #include <vector> 97 98 99 class vtkAbstractInterpolatedVelocityField; 100 class vtkCompositeDataSet; 101 class vtkDataArray; 102 class vtkDoubleArray; 103 class vtkExecutive; 104 class vtkGenericCell; 105 class vtkIdList; 106 class vtkInitialValueProblemSolver; 107 class vtkImageData; 108 class vtkIntArray; 109 class vtkPolyDataCollection; 110 class vtkPoints; 111 class vtkStreamTracer; 112 113 class VTKFILTERSFLOWPATHS_EXPORT vtkEvenlySpacedStreamlines2D : public vtkPolyDataAlgorithm 114 { 115 public: 116 vtkTypeMacro(vtkEvenlySpacedStreamlines2D,vtkPolyDataAlgorithm); 117 void PrintSelf(ostream& os, vtkIndent indent) override; 118 119 /** 120 * Construct object to start from position (0,0,0), with forward 121 * integration, terminal speed 1.0E-12, vorticity computation on, 122 * integration step size 0.5 (in cell length unit), maximum number 123 * of steps 2000, using Runge-Kutta2, and maximum propagation 1.0 124 * (in arc length unit). 125 */ 126 static vtkEvenlySpacedStreamlines2D *New(); 127 128 //@{ 129 /** 130 * Specify the starting point (seed) of the first streamline in the global 131 * coordinate system. Search must be performed to find the initial cell 132 * from which to start integration. If the seed is not specified a 133 * random position in the input data is chosen. 134 */ 135 vtkSetVector3Macro(StartPosition, double); 136 vtkGetVector3Macro(StartPosition, double); 137 //@} 138 139 //@{ 140 /** 141 * Set/get the integrator type to be used for streamline generation. 142 * The object passed is not actually used but is cloned with 143 * NewInstance in the process of integration (prototype pattern). 144 * The default is Runge-Kutta2. The integrator can also be changed 145 * using SetIntegratorType. The recognized solvers are: 146 * RUNGE_KUTTA2 = 0 147 * RUNGE_KUTTA4 = 1 148 */ 149 void SetIntegrator(vtkInitialValueProblemSolver *); 150 vtkGetObjectMacro ( Integrator, vtkInitialValueProblemSolver ); 151 void SetIntegratorType(int type); 152 int GetIntegratorType(); 153 void SetIntegratorTypeToRungeKutta2(); 154 void SetIntegratorTypeToRungeKutta4(); 155 //@} 156 157 /** 158 * Set the velocity field interpolator type to the one involving 159 * a dataset point locator. 160 */ 161 void SetInterpolatorTypeToDataSetPointLocator(); 162 163 /** 164 * Set the velocity field interpolator type to the one involving 165 * a cell locator. 166 */ 167 void SetInterpolatorTypeToCellLocator(); 168 169 /** 170 * Specify a uniform integration step unit for 171 * InitialIntegrationStep, and SeparatingDistance. Valid units are 172 * LENGTH_UNIT (1) (value is in global coordinates) and CELL_LENGTH_UNIT (2) 173 * (the value is in number of cell lengths) 174 */ 175 void SetIntegrationStepUnit( int unit ); GetIntegrationStepUnit()176 int GetIntegrationStepUnit() { return this->IntegrationStepUnit; } 177 178 //@{ 179 /** 180 * Specify the maximum number of steps for integrating a streamline. 181 */ 182 vtkSetMacro(MaximumNumberOfSteps, vtkIdType); 183 vtkGetMacro(MaximumNumberOfSteps, vtkIdType); 184 //@} 185 186 //@{ 187 /** 188 * We don't try to eliminate loops with fewer points than this. Default value 189 * is 4. 190 */ 191 vtkSetMacro(MinimumNumberOfLoopPoints, vtkIdType); 192 vtkGetMacro(MinimumNumberOfLoopPoints, vtkIdType); 193 //@} 194 195 196 //@{ 197 /** 198 * Specify the Initial step size used for line integration, expressed in 199 * IntegrationStepUnit 200 * 201 * This is the constant / fixed size for non-adaptive integration 202 * methods, i.e., RK2 and RK4 203 */ 204 vtkSetMacro(InitialIntegrationStep, double); 205 vtkGetMacro(InitialIntegrationStep, double); 206 //@} 207 208 //@{ 209 /** 210 * Specify the separation distance between streamlines expressed in 211 * IntegrationStepUnit. 212 */ 213 vtkSetMacro(SeparatingDistance, double); 214 vtkGetMacro(SeparatingDistance, double); 215 //@} 216 217 //@{ 218 /** 219 * Streamline integration is stopped if streamlines are closer than 220 * SeparatingDistance*SeparatingDistanceRatio to other streamlines. 221 */ 222 vtkSetMacro(SeparatingDistanceRatio, double); 223 vtkGetMacro(SeparatingDistanceRatio, double); 224 //@} 225 226 //@{ 227 /** 228 * Loops are considered closed if the have two points at distance less than this. 229 * This is expressed in IntegrationStepUnit. 230 */ 231 vtkSetMacro(ClosedLoopMaximumDistance, double); 232 vtkGetMacro(ClosedLoopMaximumDistance, double); 233 //@} 234 235 //@{ 236 /** 237 * The angle (in radians) between the vector created by p0p1 and the 238 * velocity in the point closing the loop. p0 is the current point 239 * and p1 is the point before that. Default value is 20 degrees in radians. 240 */ 241 vtkSetMacro(LoopAngle, double); 242 vtkGetMacro(LoopAngle, double); 243 //@} 244 245 246 //@{ 247 /** 248 * Specify the terminal speed value, below which integration is terminated. 249 */ 250 vtkSetMacro(TerminalSpeed, double); 251 vtkGetMacro(TerminalSpeed, double); 252 //@} 253 254 //@{ 255 /** 256 * Turn on/off vorticity computation at streamline points 257 * (necessary for generating proper stream-ribbons using the 258 * vtkRibbonFilter. 259 */ 260 vtkSetMacro(ComputeVorticity, bool); 261 vtkGetMacro(ComputeVorticity, bool); 262 //@} 263 264 /** 265 * The object used to interpolate the velocity field during 266 * integration is of the same class as this prototype. 267 */ 268 void SetInterpolatorPrototype( vtkAbstractInterpolatedVelocityField * ivf ); 269 270 /** 271 * Set the type of the velocity field interpolator to determine whether 272 * vtkInterpolatedVelocityField (INTERPOLATOR_WITH_DATASET_POINT_LOCATOR) or 273 * vtkCellLocatorInterpolatedVelocityField (INTERPOLATOR_WITH_CELL_LOCATOR) 274 * is employed for locating cells during streamline integration. The latter 275 * (adopting vtkAbstractCellLocator sub-classes such as vtkCellLocator and 276 * vtkModifiedBSPTree) is more robust then the former (through vtkDataSet / 277 * vtkPointSet::FindCell() coupled with vtkPointLocator). 278 */ 279 void SetInterpolatorType( int interpType ); 280 281 protected: 282 vtkEvenlySpacedStreamlines2D(); 283 ~vtkEvenlySpacedStreamlines2D() override; 284 285 /** 286 * Do we test for separating distance or a ratio of the separating distance. 287 */ 288 enum DistanceType 289 { 290 DISTANCE, 291 DISTANCE_RATIO 292 }; 293 // hide the superclass' AddInput() from the user and the compiler AddInput(vtkDataObject *)294 void AddInput(vtkDataObject *) 295 { 296 vtkErrorMacro(<< "AddInput() must be called with a vtkDataSet not a vtkDataObject."); 297 } 298 299 int RequestData(vtkInformation *, 300 vtkInformationVector **, vtkInformationVector *) override; 301 int FillInputPortInformation(int, vtkInformation *) override; 302 303 int SetupOutput(vtkInformation* inInfo, vtkInformation* outInfo); 304 int CheckInputs(vtkAbstractInterpolatedVelocityField*& func, 305 int* maxCellSize); 306 double ConvertToLength(double interval, int unit, double cellLength ); 307 308 static void GetBounds(vtkCompositeDataSet* cds, double bounds[6]); 309 void InitializeSuperposedGrid(double* bounds); 310 void AddToAllPoints(vtkPolyData* streamline); 311 void AddToCurrentPoints(vtkIdType pointId); 312 template<typename T> void InitializePoints(T& points); 313 void InitializeMinPointIds(); 314 315 static bool IsStreamlineLooping( 316 void* clientdata, 317 vtkPoints* points, vtkDataArray* velocity, int direction); 318 static bool IsStreamlineTooCloseToOthers( 319 void* clientdata, 320 vtkPoints* points, vtkDataArray* velocity, int direction); 321 template<typename CellCheckerType> 322 bool ForEachCell(double* point, CellCheckerType checker, 323 vtkPoints* points = nullptr, 324 vtkDataArray* velocity = nullptr, 325 int direction = 1); 326 template <int distanceType> 327 bool IsTooClose(double* point, vtkIdType cellId, 328 vtkPoints* points, 329 vtkDataArray* velocity, int direction); 330 bool IsLooping(double* point, vtkIdType cellId, 331 vtkPoints* points, vtkDataArray* velocity, int direction); 332 const char* GetInputArrayToProcessName(); 333 int ComputeCellLength(double* cellLength); 334 335 // starting from global x-y-z position 336 double StartPosition[3]; 337 338 double TerminalSpeed; 339 340 double InitialIntegrationStep; 341 double SeparatingDistance; 342 // SeparatingDistance can be in cell length or arc length. This member 343 // stores SeparatingDistance in arc length. It is computed when 344 // the filter executes. 345 double SeparatingDistanceArcLength; 346 double SeparatingDistanceRatio; 347 double ClosedLoopMaximumDistance; 348 // ClosedLoopMaximumDistance can be in cell length or arc length. 349 // This member stores ClosedLoopMaximumDistance in arc length. It is 350 // computed when the filter executes. 351 double ClosedLoopMaximumDistanceArcLength; 352 double LoopAngle; 353 int IntegrationStepUnit; 354 355 vtkIdType MaximumNumberOfSteps; 356 vtkIdType MinimumNumberOfStreamlinePoints; 357 vtkIdType MinimumNumberOfLoopPoints; 358 359 // Prototype showing the integrator type to be set by the user. 360 vtkInitialValueProblemSolver* Integrator; 361 362 bool ComputeVorticity; 363 364 vtkAbstractInterpolatedVelocityField * InterpolatorPrototype; 365 366 vtkCompositeDataSet* InputData; 367 // grid superposed over InputData. The grid cell height and width is 368 // SeparatingDistance 369 vtkImageData* SuperposedGrid; 370 // AllPoints[i][j] is the point for point j on the streamlines that 371 // falls over cell id i in SuperposedGrid. AllPoint[i].size() tell 372 // us how many points fall over cell id i. 373 std::vector<std::vector<std::array<double,3> > > AllPoints; 374 375 // CurrentPoints[i][j] is the point id for point j on the current streamline that 376 // falls over cell id i in SuperposedGrid. CurrentPoints[i].size() tell us 377 // how many points fall over cell id i. 378 std::vector<std::vector<vtkIdType> > CurrentPoints; 379 // Min and Max point ids stored in a cell of SuperposedGrid 380 std::vector<vtkIdType> MinPointIds; 381 // The index of the first point for the current 382 // direction. Note we integrate streamlines both forward and 383 // backward. 384 vtkIdType DirectionStart; 385 // The previous integration direction. 386 int PreviousDirection; 387 388 // queue of streamlines to be processed 389 vtkPolyDataCollection* Streamlines; 390 private: 391 vtkEvenlySpacedStreamlines2D( 392 const vtkEvenlySpacedStreamlines2D&) = delete; 393 void operator=(const vtkEvenlySpacedStreamlines2D&) = delete; 394 }; 395 396 397 #endif 398 399 // VTK-HeaderTest-Exclude: vtkEvenlySpacedStreamlines2D.h 400