1 /*=========================================================================
2 
3   Program:   Visualization Toolkit
4   Module:    vtkLagrangianParticleTracker.cxx
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 #include "vtkLagrangianParticleTracker.h"
16 
17 #include "vtkBoundingBox.h"
18 #include "vtkCellData.h"
19 #include "vtkCompositeDataIterator.h"
20 #include "vtkCompositeDataSet.h"
21 #include "vtkDataSet.h"
22 #include "vtkDataSetSurfaceFilter.h"
23 #include "vtkDoubleArray.h"
24 #include "vtkExecutive.h"
25 #include "vtkIdList.h"
26 #include "vtkInformation.h"
27 #include "vtkInformationVector.h"
28 #include "vtkLagrangianMatidaIntegrationModel.h"
29 #include "vtkLagrangianParticle.h"
30 #include "vtkLongLongArray.h"
31 #include "vtkNew.h"
32 #include "vtkObjectFactory.h"
33 #include "vtkPointData.h"
34 #include "vtkPoints.h"
35 #include "vtkPolyData.h"
36 #include "vtkPolyDataNormals.h"
37 #include "vtkPolyLine.h"
38 #include "vtkPolygon.h"
39 #include "vtkRungeKutta2.h"
40 #include "vtkStreamingDemandDrivenPipeline.h"
41 #include "vtkVoxel.h"
42 
43 #include <algorithm>
44 #include <limits>
45 #include <sstream>
46 
47 vtkObjectFactoryNewMacro(vtkLagrangianParticleTracker);
48 vtkCxxSetObjectMacro(vtkLagrangianParticleTracker, IntegrationModel, vtkLagrangianBasicIntegrationModel);
49 vtkCxxSetObjectMacro(vtkLagrangianParticleTracker, Integrator, vtkInitialValueProblemSolver);
50 
51 //---------------------------------------------------------------------------
vtkLagrangianParticleTracker()52 vtkLagrangianParticleTracker::vtkLagrangianParticleTracker()
53 {
54   this->IntegrationModel = vtkLagrangianMatidaIntegrationModel::New();
55   this->Integrator = vtkRungeKutta2::New();
56 
57   this->SetNumberOfInputPorts(3);
58   this->SetNumberOfOutputPorts(2);
59 
60   this->CellLengthComputationMode = STEP_LAST_CELL_LENGTH;
61   this->AdaptiveStepReintegration = false;
62   this->StepFactor = 1.0;
63   this->StepFactorMin = 0.5;
64   this->StepFactorMax = 1.5;
65   this->MaximumNumberOfSteps = 100;
66   this->MaximumIntegrationTime = -1.0;
67 
68   this->MinimumVelocityMagnitude = 0.001;
69   this->MinimumReductionFactor = 1.1;
70 
71   this->UseParticlePathsRenderingThreshold = false;
72   this->GeneratePolyVertexInteractionOutput = false;
73   this->ParticlePathsRenderingPointsThreshold = 100;
74 
75   this->ParticleCounter = 0;
76 
77   this->FlowCache = nullptr;
78   this->FlowTime = 0;
79   this->SurfacesCache = nullptr;
80   this->SurfacesTime = 0;
81 }
82 
83 //---------------------------------------------------------------------------
~vtkLagrangianParticleTracker()84 vtkLagrangianParticleTracker::~vtkLagrangianParticleTracker()
85 {
86   this->SetIntegrator(nullptr);
87   this->SetIntegrationModel(nullptr);
88 }
89 
90 //---------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)91 void vtkLagrangianParticleTracker::PrintSelf(ostream& os, vtkIndent indent)
92 {
93   this->Superclass::PrintSelf(os, indent);
94   if (this->IntegrationModel)
95   {
96     os << indent << "IntegrationModel: " << endl;
97     this->IntegrationModel->PrintSelf(os, indent.GetNextIndent());
98   }
99   else
100   {
101     os << indent << "IntegrationModel: " << this->IntegrationModel << endl;
102   }
103   if (this->Integrator)
104   {
105     os << indent << "Integrator: " << endl;
106     this->Integrator->PrintSelf(os, indent.GetNextIndent());
107   }
108   else
109   {
110     os << indent << "Integrator: " << this->Integrator << endl;
111   }
112   os << indent << "CellLengthComputationMode: "
113     << this->CellLengthComputationMode << endl;
114   os << indent << "StepFactor: " << this->StepFactor << endl;
115   os << indent << "StepFactorMin: " << this->StepFactorMin << endl;
116   os << indent << "StepFactorMax: " << this->StepFactorMax << endl;
117   os << indent << "MaximumNumberOfSteps: " << this->MaximumNumberOfSteps << endl;
118   os << indent << "MaximumIntegrationTime: " << this->MaximumIntegrationTime << endl;
119   os << indent << "AdaptiveStepReintegration: " << this->AdaptiveStepReintegration << endl;
120   os << indent << "UseParticlePathsRenderingThreshold: "
121     << this->UseParticlePathsRenderingThreshold << endl;
122   os << indent << "ParticlePathsRenderingPointsThreshold: "
123     << this->ParticlePathsRenderingPointsThreshold << endl;
124   os << indent << "MinimumVelocityMagnitude: " << this->MinimumVelocityMagnitude << endl;
125   os << indent << "MinimumReductionFactor: " << this->MinimumReductionFactor << endl;
126   os << indent << "ParticleCounter: " << this->ParticleCounter << endl;
127 }
128 
129 //---------------------------------------------------------------------------
SetSourceConnection(vtkAlgorithmOutput * algInput)130 void vtkLagrangianParticleTracker::SetSourceConnection(vtkAlgorithmOutput* algInput)
131 {
132   this->SetInputConnection(1, algInput);
133 }
134 
135 //---------------------------------------------------------------------------
SetSourceData(vtkDataObject * source)136 void vtkLagrangianParticleTracker::SetSourceData(vtkDataObject *source)
137 {
138   this->SetInputData(1, source);
139 }
140 
141 //---------------------------------------------------------------------------
GetSource()142 vtkDataObject* vtkLagrangianParticleTracker::GetSource()
143 {
144   if (this->GetNumberOfInputConnections(1) < 1)
145   {
146     return nullptr;
147   }
148   return vtkDataObject::SafeDownCast(this->GetExecutive()->GetInputData(1, 0));
149 }
150 
151 //---------------------------------------------------------------------------
SetSurfaceConnection(vtkAlgorithmOutput * algOutput)152 void vtkLagrangianParticleTracker::SetSurfaceConnection(vtkAlgorithmOutput* algOutput)
153 {
154   this->SetInputConnection(2, algOutput);
155 }
156 
157 //---------------------------------------------------------------------------
SetSurfaceData(vtkDataObject * surface)158 void vtkLagrangianParticleTracker::SetSurfaceData(vtkDataObject *surface)
159 {
160   this->SetInputData(2, surface);
161 }
162 
163 //---------------------------------------------------------------------------
GetSurface()164 vtkDataObject* vtkLagrangianParticleTracker::GetSurface()
165 {
166   if (this->GetNumberOfInputConnections(2) < 1)
167   {
168     return nullptr;
169   }
170   return this->GetExecutive()->GetInputData(2, 0);
171 }
172 
173 //---------------------------------------------------------------------------
FillInputPortInformation(int port,vtkInformation * info)174 int vtkLagrangianParticleTracker::FillInputPortInformation(int port,
175   vtkInformation *info)
176 {
177   if (port == 2)
178   {
179     info->Set(vtkAlgorithm::INPUT_IS_OPTIONAL(), 1);
180   }
181   return this->Superclass::FillInputPortInformation(port, info);
182 }
183 
184 //----------------------------------------------------------------------------
FillOutputPortInformation(int port,vtkInformation * info)185 int vtkLagrangianParticleTracker::FillOutputPortInformation(int port,
186   vtkInformation* info)
187 {
188   if (port == 0)
189   {
190     info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkPolyData");
191   }
192   return this->Superclass::FillOutputPortInformation(port, info);
193 }
194 
195 //----------------------------------------------------------------------------
RequestDataObject(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)196 int vtkLagrangianParticleTracker::RequestDataObject(
197   vtkInformation* vtkNotUsed(request),
198   vtkInformationVector** inputVector,
199   vtkInformationVector* outputVector)
200 {
201   // Create a polydata output
202   vtkInformation* info = outputVector->GetInformationObject(0);
203   vtkNew<vtkPolyData> particlePathsOutput;
204   info->Set(vtkDataObject::DATA_OBJECT(), particlePathsOutput);
205 
206   // Create a surface interaction output
207   // First check for composite
208   vtkInformation* inInfo = inputVector[2]->GetInformationObject(0);
209   info = outputVector->GetInformationObject(1);
210   if (inInfo)
211   {
212     vtkDataObject *input = vtkDataObject::SafeDownCast(
213       inInfo->Get(vtkDataObject::DATA_OBJECT()));
214     if (input)
215     {
216       vtkCompositeDataSet *hdInput = vtkCompositeDataSet::SafeDownCast(input);
217       if (hdInput)
218       {
219         vtkDataObject* interactionOutput = input->NewInstance();
220         info->Set(vtkDataObject::DATA_OBJECT(), interactionOutput);
221         interactionOutput->Delete();
222         return 1;
223       }
224     }
225   }
226   // In any other case, create a polydata
227   vtkNew<vtkPolyData> interactionOutput;
228   info->Set(vtkDataObject::DATA_OBJECT(), interactionOutput);
229   return 1;
230 }
231 
232 //---------------------------------------------------------------------------
RequestData(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)233 int vtkLagrangianParticleTracker::RequestData(
234   vtkInformation* vtkNotUsed(request),
235   vtkInformationVector **inputVector,
236   vtkInformationVector *outputVector)
237 {
238   // Initialize inputs
239   vtkDataObject* flow = nullptr;
240   vtkDataObject* seeds = nullptr;
241   vtkDataObject* surfaces = nullptr;
242   std::queue<vtkLagrangianParticle*> particlesQueue;
243 
244   if (!this->IntegrationModel)
245   {
246     vtkErrorMacro(<< "Integration Model is nullptr, cannot integrate");
247     return 0;
248   }
249   this->IntegrationModel->SetTracker(this);
250 
251   vtkNew<vtkPointData> seedData;
252   if (!this->InitializeInputs(inputVector, flow, seeds, surfaces,
253     particlesQueue, seedData))
254   {
255     vtkErrorMacro(<< "Cannot initialize inputs");
256     return 0;
257   }
258 
259   // Initialize outputs
260   vtkPolyData* particlePathsOutput;
261   vtkDataObject* interactionOutput;
262   if (!this->InitializeOutputs(outputVector, seedData,
263     static_cast<vtkIdType>(particlesQueue.size()), surfaces,
264     particlePathsOutput, interactionOutput))
265   {
266     vtkErrorMacro(<< "Cannot initialize outputs");
267     return 0;
268   }
269 
270   // Let model a chance to change the particles or compute things
271   // before integration.
272   this->IntegrationModel->PreIntegrate(particlesQueue);
273 
274   // Integrate each particle
275   while (!this->GetAbortExecute())
276   {
277     // Check for particle feed
278     this->GetParticleFeed(particlesQueue);
279     if (particlesQueue.empty())
280     {
281       break;
282     }
283 
284     // Recover particle
285     vtkLagrangianParticle* particle = particlesQueue.front();
286     particlesQueue.pop();
287 
288     // Create polyLine output cell
289     vtkNew<vtkPolyLine> particlePath;
290 
291     // Integrate
292     this->Integrate(particle, particlesQueue, particlePathsOutput,
293       particlePath->GetPointIds(), interactionOutput);
294 
295     if (particlePath->GetPointIds()->GetNumberOfIds() == 1)
296     {
297       particlePath->GetPointIds()->InsertNextId(particlePath->GetPointId(0));
298     }
299 
300     // Duplicate single point particle paths, to avoid degenerated lines.
301     if (particlePath->GetPointIds()->GetNumberOfIds() > 0)
302     {
303       // Add particle path or vertex to cell array
304       particlePathsOutput->GetLines()->InsertNextCell(particlePath);
305       this->InsertPathData(particle, particlePathsOutput->GetCellData());
306       this->IntegrationModel->InsertModelPathData(particle, particlePathsOutput->GetCellData());
307 
308       // Insert data from seed data only on not yet written arrays
309       this->InsertSeedData(particle, particlePathsOutput->GetCellData());
310     }
311 
312     // Delete integrated particle
313     delete particle;
314   }
315 
316   // Abort if necessary
317   if (this->GetAbortExecute())
318   {
319     // delete all remaining particle
320     while (!particlesQueue.empty())
321     {
322       vtkLagrangianParticle* particle = particlesQueue.front();
323       particlesQueue.pop();
324       delete particle;
325     }
326   }
327   // Finalize outputs
328   else if (!this->FinalizeOutputs(particlePathsOutput, interactionOutput))
329   {
330     vtkErrorMacro(<< "Cannot Finalize outputs");
331     return 0;
332   }
333   return 1;
334 }
335 
336 //---------------------------------------------------------------------------
CheckParticlePathsRenderingThreshold(vtkPolyData * particlePathsOutput)337 bool vtkLagrangianParticleTracker::CheckParticlePathsRenderingThreshold(vtkPolyData* particlePathsOutput)
338 {
339   return this->UseParticlePathsRenderingThreshold &&
340     particlePathsOutput->GetNumberOfPoints() > this->ParticlePathsRenderingPointsThreshold;
341 }
342 
343 //---------------------------------------------------------------------------
GetMTime()344 vtkMTimeType vtkLagrangianParticleTracker::GetMTime()
345 {
346   // Take integrator and integration model MTime into account
347   return std::max(this->Superclass::GetMTime(),
348     std::max(this->IntegrationModel ? this->IntegrationModel->GetMTime() : 0,
349       this->Integrator ? this->Integrator->GetMTime() : 0));
350 }
351 
352 //---------------------------------------------------------------------------
GetNewParticleId()353 vtkIdType vtkLagrangianParticleTracker::GetNewParticleId()
354 {
355   vtkIdType id = this->ParticleCounter;
356   this->ParticleCounter++;
357   return id;
358 }
359 
360 //---------------------------------------------------------------------------
InitializeInputs(vtkInformationVector ** inputVector,vtkDataObject * & flow,vtkDataObject * & seeds,vtkDataObject * & surfaces,std::queue<vtkLagrangianParticle * > & particlesQueue,vtkPointData * seedData)361 bool vtkLagrangianParticleTracker::InitializeInputs(vtkInformationVector **inputVector,
362   vtkDataObject*& flow, vtkDataObject*& seeds, vtkDataObject*& surfaces,
363   std::queue<vtkLagrangianParticle*>& particlesQueue, vtkPointData* seedData)
364 {
365   // Initialize flow
366   vtkInformation* flowInInfo = inputVector[0]->GetInformationObject(0);
367   flow = flowInInfo->Get(vtkDataObject::DATA_OBJECT());
368   vtkBoundingBox bounds;
369   if (!this->InitializeFlow(flow, &bounds))
370   {
371     vtkErrorMacro(<< "Could not initialize flow, aborting.");
372     return false;
373   }
374 
375   // Recover seeds
376   vtkInformation *seedsInInfo = inputVector[1]->GetInformationObject(0);
377   seeds = vtkDataObject::SafeDownCast(seedsInInfo->Get(vtkDataObject::DATA_OBJECT()));
378   if (!seeds)
379   {
380     vtkErrorMacro(<< "Cannot recover seeds, aborting.");
381     return false;
382   }
383 
384   // Configure integrator (required before particle initialization)
385   this->Integrator->SetFunctionSet(this->IntegrationModel);
386 
387   // Initialize Particles
388   if (!this->InitializeParticles(&bounds, seeds, particlesQueue, seedData))
389   {
390     vtkErrorMacro(<< "Could not initialize particles, aborting.");
391     return false;
392   }
393 
394   // Recover surfaces
395   vtkInformation* surfacesInInfo = inputVector[2]->GetInformationObject(0);
396   if (surfacesInInfo != nullptr)
397   {
398     surfaces = surfacesInInfo->Get(vtkDataObject::DATA_OBJECT());
399     if (this->UpdateSurfaceCacheIfNeeded(surfaces))
400     {
401       this->InitializeSurface(surfaces);
402     }
403   }
404   return true;
405 }
406 
407 //---------------------------------------------------------------------------
InitializeOutputs(vtkInformationVector * outputVector,vtkPointData * seedData,vtkIdType numberOfSeeds,vtkDataObject * surfaces,vtkPolyData * & particlePathsOutput,vtkDataObject * & interactionOutput)408 bool vtkLagrangianParticleTracker::InitializeOutputs(
409   vtkInformationVector *outputVector, vtkPointData* seedData,
410   vtkIdType numberOfSeeds, vtkDataObject* surfaces,
411   vtkPolyData*& particlePathsOutput, vtkDataObject*& interactionOutput)
412 {
413   if (!this->InitializePathsOutput(outputVector, seedData, numberOfSeeds, particlePathsOutput))
414   {
415     return false;
416   }
417   if (!this->InitializeInteractionOutput(outputVector, seedData, surfaces, interactionOutput))
418   {
419     return false;
420   }
421   return true;
422 }
423 
424 //---------------------------------------------------------------------------
InitializePathsOutput(vtkInformationVector * outputVector,vtkPointData * seedData,vtkIdType numberOfSeeds,vtkPolyData * & particlePathsOutput)425 bool vtkLagrangianParticleTracker::InitializePathsOutput(vtkInformationVector *outputVector,
426   vtkPointData* seedData, vtkIdType numberOfSeeds, vtkPolyData*& particlePathsOutput)
427 {
428   // Prepare path output
429   vtkInformation* particleOutInfo = outputVector->GetInformationObject(0);
430   particlePathsOutput = vtkPolyData::SafeDownCast(particleOutInfo->Get(
431     vtkPolyData::DATA_OBJECT()));
432   if (!particlePathsOutput)
433   {
434     vtkErrorMacro(<< "Cannot find a vtkPolyData particle paths output. aborting");
435     return false;
436   }
437 
438   // Set information keys
439   particlePathsOutput->GetInformation()->Set(
440     vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES(),
441     particleOutInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES()));
442   particlePathsOutput->GetInformation()->Set(
443     vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER(),
444     particleOutInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER()));
445   particlePathsOutput->GetInformation()->Set(
446     vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS(),
447     particleOutInfo->Get(
448       vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS()));
449 
450   vtkNew<vtkPoints> particlePathsPoints;
451   vtkNew<vtkCellArray> particlePaths;
452   vtkNew<vtkCellArray> particleVerts;
453   particlePathsOutput->SetPoints(particlePathsPoints);
454   particlePathsOutput->SetLines(particlePaths);
455   particlePathsOutput->SetVerts(particleVerts);
456 
457   // Prepare particle paths output point data
458   vtkCellData* particlePathsCellData = particlePathsOutput->GetCellData();
459   particlePathsCellData->CopyStructure(seedData);
460   this->InitializePathData(particlePathsCellData);
461   this->IntegrationModel->InitializeModelPathData(particlePathsCellData);
462 
463   // Initialize Particle Paths Point Data
464   vtkPointData* particlePathsPointData = particlePathsOutput->GetPointData();
465   this->InitializeParticleData(particlePathsPointData, numberOfSeeds);
466 
467   // Initialize particle data from integration model, if any
468   this->IntegrationModel->
469     InitializeVariablesParticleData(particlePathsPointData, numberOfSeeds);
470   return true;
471 }
472 
473 //---------------------------------------------------------------------------
InitializeInteractionOutput(vtkInformationVector * outputVector,vtkPointData * seedData,vtkDataObject * surfaces,vtkDataObject * & interactionOutput)474 bool vtkLagrangianParticleTracker::InitializeInteractionOutput(
475   vtkInformationVector *outputVector, vtkPointData* seedData,
476   vtkDataObject* surfaces, vtkDataObject*& interactionOutput)
477 {
478   // Prepare interaction output
479   vtkInformation* particleOutInfo = outputVector->GetInformationObject(1);
480   interactionOutput = particleOutInfo->Get(vtkDataObject::DATA_OBJECT());
481   if (!interactionOutput)
482   {
483     vtkErrorMacro(<< "Cannot find a vtkDataObject particle interaction output. aborting");
484     return false;
485   }
486 
487   // Set information keys
488   interactionOutput->GetInformation()->Set(
489     vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES(),
490     particleOutInfo->Get(
491       vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES()));
492   interactionOutput->GetInformation()->Set(
493     vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER(),
494     particleOutInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER()));
495   interactionOutput->GetInformation()->Set(
496     vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS(),
497     particleOutInfo->Get(
498       vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS()));
499 
500   // Check surfaces dataset type
501   vtkCompositeDataSet *hdInput = vtkCompositeDataSet::SafeDownCast(surfaces);
502   vtkDataSet* dsInput = vtkDataSet::SafeDownCast(surfaces);
503   if (hdInput)
504   {
505     // Composite data
506     vtkCompositeDataSet* hdOutput = vtkCompositeDataSet::SafeDownCast(interactionOutput);
507     if (!hdOutput)
508     {
509       vtkErrorMacro(<< "Cannot find composite interaction output, aborting");
510       return false;
511     }
512 
513     hdOutput->CopyStructure(hdInput);
514     vtkSmartPointer<vtkCompositeDataIterator> iter;
515     iter.TakeReference(hdInput->NewIterator());
516     for (iter->InitTraversal(); !iter->IsDoneWithTraversal(); iter->GoToNextItem())
517     {
518       vtkNew<vtkPolyData> pd;
519       vtkNew<vtkCellArray> cells;
520       vtkNew<vtkPoints> points;
521       pd->SetPoints(points);
522       pd->GetPointData()->CopyStructure(seedData);
523       this->InitializePathData(pd->GetPointData());
524       this->IntegrationModel->InitializeModelPathData(pd->GetPointData());
525       this->InitializeInteractionData(pd->GetPointData());
526       this->InitializeParticleData(pd->GetPointData());
527       this->IntegrationModel->
528         InitializeVariablesParticleData(pd->GetPointData());
529       hdOutput->SetDataSet(iter, pd);
530     }
531   }
532   else if (dsInput)
533   {
534     vtkPolyData* pd = vtkPolyData::SafeDownCast(interactionOutput);
535     if (!pd)
536     {
537       vtkErrorMacro(<< "Cannot find polydata interaction output, aborting");
538       return false;
539     }
540 
541     vtkNew<vtkPoints> points;
542     vtkNew<vtkCellArray> cells;
543     pd->SetPoints(points);
544     pd->GetPointData()->CopyStructure(seedData);
545     this->InitializePathData(pd->GetPointData());
546     this->IntegrationModel->InitializeModelPathData(pd->GetPointData());
547     this->InitializeInteractionData(pd->GetPointData());
548     this->InitializeParticleData(pd->GetPointData());
549     this->IntegrationModel->
550       InitializeVariablesParticleData(pd->GetPointData());
551   }
552   return true;
553 }
554 
555 //---------------------------------------------------------------------------
InitializeParticleData(vtkFieldData * particleData,int maxTuple)556 void vtkLagrangianParticleTracker::InitializeParticleData(vtkFieldData* particleData, int maxTuple)
557 {
558   vtkNew<vtkIntArray> particleStepNumArray;
559   particleStepNumArray->SetName("StepNumber");
560   particleStepNumArray->SetNumberOfComponents(1);
561   particleStepNumArray->Allocate(maxTuple);
562   particleData->AddArray(particleStepNumArray);
563 
564   vtkNew<vtkDoubleArray> particleVelArray;
565   particleVelArray->SetName("ParticleVelocity");
566   particleVelArray->SetNumberOfComponents(3);
567   particleVelArray->Allocate(maxTuple*3);
568   particleData->AddArray(particleVelArray);
569 
570   vtkNew<vtkDoubleArray> particleIntegrationTimeArray;
571   particleIntegrationTimeArray->SetName("IntegrationTime");
572   particleIntegrationTimeArray->SetNumberOfComponents(1);
573   particleIntegrationTimeArray->Allocate(maxTuple);
574   particleData->AddArray(particleIntegrationTimeArray);
575 }
576 
577 //---------------------------------------------------------------------------
InitializePathData(vtkFieldData * data)578 void vtkLagrangianParticleTracker::InitializePathData(vtkFieldData* data)
579 {
580   vtkNew<vtkLongLongArray> particleIdArray;
581   particleIdArray->SetName("Id");
582   particleIdArray->SetNumberOfComponents(1);
583   data->AddArray(particleIdArray);
584 
585   vtkNew<vtkLongLongArray> particleParentIdArray;
586   particleParentIdArray->SetName("ParentId");
587   particleParentIdArray->SetNumberOfComponents(1);
588   data->AddArray(particleParentIdArray);
589 
590   vtkNew<vtkLongLongArray> particleSeedIdArray;
591   particleSeedIdArray->SetName("SeedId");
592   particleSeedIdArray->SetNumberOfComponents(1);
593   data->AddArray(particleSeedIdArray);
594 
595   vtkNew<vtkIntArray> particleTerminationArray;
596   particleTerminationArray->SetName("Termination");
597   particleTerminationArray->SetNumberOfComponents(1);
598   data->AddArray(particleTerminationArray);
599 }
600 
601 //---------------------------------------------------------------------------
InitializeInteractionData(vtkFieldData * data)602 void vtkLagrangianParticleTracker::InitializeInteractionData(vtkFieldData* data)
603 {
604   vtkNew<vtkIntArray> interactionArray;
605   interactionArray->SetName("Interaction");
606   interactionArray->SetNumberOfComponents(1);
607   data->AddArray(interactionArray);
608 }
609 
610 //---------------------------------------------------------------------------
FinalizeOutputs(vtkPolyData * particlePathsOutput,vtkDataObject * interactionOutput)611 bool vtkLagrangianParticleTracker::FinalizeOutputs(
612   vtkPolyData* particlePathsOutput,
613   vtkDataObject* interactionOutput)
614 {
615   // Recover structures
616   vtkPointData* particlePathsPointData = particlePathsOutput->GetPointData();
617   vtkPoints* particlePathsPoints = particlePathsOutput->GetPoints();
618 
619   // Squeeze and resize point data
620   for (int i = 0; i < particlePathsPointData->GetNumberOfArrays(); i++)
621   {
622     vtkDataArray* array = particlePathsPointData->GetArray(i);
623     array->Resize(particlePathsPoints->GetNumberOfPoints());
624     array->Squeeze();
625   }
626 
627   // Insert interaction poly-vertex cell
628   vtkCompositeDataSet *hd = vtkCompositeDataSet::SafeDownCast(interactionOutput);
629   vtkPolyData* pd = vtkPolyData::SafeDownCast(interactionOutput);
630   if (hd)
631   {
632     vtkSmartPointer<vtkCompositeDataIterator> iter;
633     iter.TakeReference(hd->NewIterator());
634     for (iter->InitTraversal(); !iter->IsDoneWithTraversal(); iter->GoToNextItem())
635     {
636       vtkPolyData* pdBlock = vtkPolyData::SafeDownCast(hd->GetDataSet(iter));
637       if (!pdBlock)
638       {
639         vtkErrorMacro(<< "Cannot recover interaction output, something went wrong");
640         return false;
641       }
642       if (this->GeneratePolyVertexInteractionOutput)
643       {
644         this->InsertPolyVertexCell(pdBlock);
645       }
646       else
647       {
648         this->InsertVertexCells(pdBlock);
649       }
650     }
651   }
652   else if (pd)
653   {
654     if (this->GeneratePolyVertexInteractionOutput)
655     {
656       this->InsertPolyVertexCell(pd);
657     }
658     else
659     {
660       this->InsertVertexCells(pd);
661     }
662   }
663 
664   // Enable model post processing
665   this->IntegrationModel->FinalizeOutputs(particlePathsOutput, interactionOutput);
666 
667   // Optional paths rendering
668   if (this->CheckParticlePathsRenderingThreshold(particlePathsOutput))
669   {
670     particlePathsOutput->Initialize();
671   }
672   return true;
673 }
674 
675 //---------------------------------------------------------------------------
InsertPolyVertexCell(vtkPolyData * polydata)676 void vtkLagrangianParticleTracker::InsertPolyVertexCell(vtkPolyData* polydata)
677 {
678   // Insert a vertex cell for each point
679   vtkIdType nPoint = polydata->GetNumberOfPoints();
680   if (nPoint > 0)
681   {
682     vtkNew<vtkCellArray> polyVertex;
683     polyVertex->Allocate(polyVertex->EstimateSize(1, nPoint));
684     polyVertex->InsertNextCell(nPoint);
685     for (vtkIdType i = 0; i < nPoint; i++)
686     {
687       polyVertex->InsertCellPoint(i);
688     }
689     polydata->SetVerts(polyVertex);
690   }
691 }
692 
693 //---------------------------------------------------------------------------
InsertVertexCells(vtkPolyData * polydata)694 void vtkLagrangianParticleTracker::InsertVertexCells(vtkPolyData* polydata)
695 {
696   // Insert a vertex cell for each point
697   vtkIdType nPoint = polydata->GetNumberOfPoints();
698   if (nPoint > 0)
699   {
700     vtkNew<vtkCellArray> polyVertex;
701     polyVertex->Allocate(polyVertex->EstimateSize(1, nPoint));
702     for (vtkIdType i = 0; i < nPoint; i++)
703     {
704       polyVertex->InsertNextCell(1);
705       polyVertex->InsertCellPoint(i);
706     }
707     polydata->SetVerts(polyVertex);
708   }
709 }
710 
711 //---------------------------------------------------------------------------
InitializeFlow(vtkDataObject * input,vtkBoundingBox * bounds)712 bool vtkLagrangianParticleTracker::InitializeFlow(vtkDataObject* input, vtkBoundingBox* bounds)
713 {
714   // Check for updated cache
715   if (input == this->FlowCache &&
716      input->GetMTime() <= this->FlowTime &&
717      this->IntegrationModel->GetLocatorsBuilt())
718   {
719     bounds->Reset();
720     bounds->AddBox(this->FlowBoundsCache);
721     return true;
722   }
723 
724   // No Cache, do the initialization
725   // Clear previously setup flow
726   this->IntegrationModel->ClearDataSets();
727 
728   // Check flow dataset type
729   vtkCompositeDataSet *hdInput = vtkCompositeDataSet::SafeDownCast(input);
730   vtkDataSet* dsInput = vtkDataSet::SafeDownCast(input);
731   if (hdInput)
732   {
733     // Composite data
734     vtkSmartPointer<vtkCompositeDataIterator> iter;
735     iter.TakeReference(hdInput->NewIterator());
736     for (iter->InitTraversal(); !iter->IsDoneWithTraversal(); iter->GoToNextItem())
737     {
738       vtkDataSet *ds = vtkDataSet::SafeDownCast(iter->GetCurrentDataObject());
739       if (ds)
740       {
741         // Add each leaf to the integration model
742         this->IntegrationModel->AddDataSet(ds);
743         ds->ComputeBounds();
744         bounds->AddBounds(ds->GetBounds());
745       }
746     }
747   }
748   else if (dsInput)
749   {
750     // Add dataset to integration model
751     this->IntegrationModel->AddDataSet(dsInput);
752     dsInput->ComputeBounds();
753     bounds->AddBounds(dsInput->GetBounds());
754   }
755   else
756   {
757     vtkErrorMacro(<< "This filter cannot handle input of type: " <<
758                   (input ? input->GetClassName() : "(none)"));
759     return false;
760   }
761   this->IntegrationModel->SetLocatorsBuilt(true);
762   this->FlowCache = input;
763   this->FlowTime = input->GetMTime();
764   this->FlowBoundsCache.Reset();
765   this->FlowBoundsCache.AddBox(*bounds);
766   return true;
767 }
768 
769 //---------------------------------------------------------------------------
UpdateSurfaceCacheIfNeeded(vtkDataObject * & surfaces)770 bool vtkLagrangianParticleTracker::UpdateSurfaceCacheIfNeeded(vtkDataObject*& surfaces)
771 {
772   if (surfaces != this->SurfacesCache || surfaces->GetMTime() > this->SurfacesTime)
773   {
774     this->SurfacesCache = surfaces;
775     this->SurfacesTime = surfaces->GetMTime();
776     return true;
777   }
778   return false;
779 }
780 
781 //---------------------------------------------------------------------------
InitializeSurface(vtkDataObject * & surfaces)782 void vtkLagrangianParticleTracker::InitializeSurface(vtkDataObject*& surfaces)
783 {
784   // Clear previously setup surfaces
785   this->IntegrationModel->ClearDataSets(/*surface*/ true);
786 
787   // Check surfaces dataset type
788   vtkCompositeDataSet *hdInput = vtkCompositeDataSet::SafeDownCast(surfaces);
789   vtkDataSet* dsInput = vtkDataSet::SafeDownCast(surfaces);
790 
791   if (hdInput)
792   {
793     // Composite data
794     vtkSmartPointer<vtkCompositeDataIterator> iter;
795     iter.TakeReference(hdInput->NewIterator());
796     for (iter->InitTraversal(); !iter->IsDoneWithTraversal(); iter->GoToNextItem())
797     {
798       vtkDataSet* ds = vtkDataSet::SafeDownCast(iter->GetCurrentDataObject());
799       if (ds)
800       {
801         vtkPolyData* pd = vtkPolyData::SafeDownCast(iter->GetCurrentDataObject());
802         vtkNew<vtkDataSetSurfaceFilter> surfaceFilter;
803         if (pd == nullptr)
804         {
805           surfaceFilter->SetInputData(ds);
806           surfaceFilter->Update();
807           pd = surfaceFilter->GetOutput();
808         }
809 
810         // Add each leaf to the integration model surfaces
811         // Compute normals if non-present
812         vtkNew<vtkPolyDataNormals> normals;
813         if (!pd->GetCellData()->GetNormals())
814         {
815           normals->ComputePointNormalsOff();
816           normals->ComputeCellNormalsOn();
817           normals->SetInputData(pd);
818           normals->Update();
819           pd = normals->GetOutput();
820         }
821         if (pd->GetNumberOfCells() > 0)
822         {
823           this->IntegrationModel->AddDataSet(pd, /*surface*/ true, iter->GetCurrentFlatIndex());
824         }
825       }
826     }
827   }
828   else if (dsInput)
829   {
830     vtkPolyData* pd = vtkPolyData::SafeDownCast(dsInput);
831     vtkNew<vtkDataSetSurfaceFilter> surfaceFilter;
832     if (pd == nullptr)
833     {
834       surfaceFilter->SetInputData(dsInput);
835       surfaceFilter->Update();
836       pd = surfaceFilter->GetOutput();
837     }
838 
839     // Add surface to integration model
840     // Compute normals if non-present
841     vtkNew<vtkPolyDataNormals> normals;
842     if (!pd->GetCellData()->GetNormals())
843     {
844       normals->ComputePointNormalsOff();
845       normals->ComputeCellNormalsOn();
846       normals->SetInputData(pd);
847       normals->Update();
848       pd = normals->GetOutput();
849     }
850     if (pd->GetNumberOfCells() > 0)
851     {
852       this->IntegrationModel->AddDataSet(pd, /*surface*/ true);
853     }
854   }
855 }
856 
857 //---------------------------------------------------------------------------
InitializeParticles(const vtkBoundingBox * bounds,vtkDataObject * seeds,std::queue<vtkLagrangianParticle * > & particles,vtkPointData * seedData)858 bool vtkLagrangianParticleTracker::InitializeParticles(
859   const vtkBoundingBox* bounds, vtkDataObject* seeds,
860   std::queue<vtkLagrangianParticle*>& particles, vtkPointData* seedData)
861 {
862   // Sanity check
863   if (seeds == nullptr)
864   {
865     vtkErrorMacro(<< "Cannot generate Particles without seeds");
866     return false;
867   }
868 
869   // Check seed dataset type
870   vtkCompositeDataSet *hdInput = vtkCompositeDataSet::SafeDownCast(seeds);
871   vtkDataSet* actualSeeds = vtkDataSet::SafeDownCast(seeds);
872   if (hdInput)
873   {
874     // Composite data
875     vtkSmartPointer<vtkCompositeDataIterator> iter;
876     iter.TakeReference(hdInput->NewIterator());
877     bool leafFound = false;
878     for (iter->InitTraversal(); !iter->IsDoneWithTraversal(); iter->GoToNextItem())
879     {
880       vtkDataSet *ds = vtkDataSet::SafeDownCast(iter->GetCurrentDataObject());
881       if (ds)
882       {
883         // We show the warning only when the input contains more than one leaf
884         if (leafFound)
885         {
886           vtkWarningMacro("Only the first block of seeds have been used to "
887             "generate seeds, other blocks are ignored");
888           break;
889         }
890         actualSeeds = ds;
891         leafFound = true;
892       }
893     }
894   }
895 
896   if (!actualSeeds)
897   {
898     vtkErrorMacro(<< "This filter cannot handle input of type: " <<
899       (seeds ? seeds->GetClassName() : "(none)"));
900     return false;
901   }
902 
903   // Recover data
904   int nVar = this->IntegrationModel->GetNumberOfIndependentVariables();
905   seedData->DeepCopy(actualSeeds->GetPointData());
906 
907   vtkDataArray* initialVelocities = nullptr;
908   vtkDataArray* initialIntegrationTimes = nullptr;
909 
910   if (actualSeeds->GetNumberOfPoints() > 0)
911   {
912     // Recover initial velocities, index 0
913     initialVelocities = vtkDataArray::SafeDownCast(
914       this->IntegrationModel->GetSeedArray(0, actualSeeds->GetPointData()));
915     if (initialVelocities == nullptr)
916     {
917       vtkErrorMacro(<< "initialVelocity is not set in particle data, "
918         "unable to initialize particles!");
919       return false;
920     }
921 
922     // Recover initial integration time if any, index 1
923     if (this->IntegrationModel->GetUseInitialIntegrationTime())
924     {
925       initialIntegrationTimes = vtkDataArray::SafeDownCast(
926         this->IntegrationModel->GetSeedArray(1, actualSeeds->GetPointData()));
927       if (initialVelocities == nullptr)
928       {
929         vtkWarningMacro("initialIntegrationTimes is not set in particle data, "
930           "initial integration time set to zero!");
931       }
932     }
933   }
934 
935   // Create one particle for each point
936   this->GenerateParticles(bounds, actualSeeds, initialVelocities,
937     initialIntegrationTimes, seedData, nVar, particles);
938   return true;
939 }
940 
941 //---------------------------------------------------------------------------
GenerateParticles(const vtkBoundingBox * vtkNotUsed (bounds),vtkDataSet * seeds,vtkDataArray * initialVelocities,vtkDataArray * initialIntegrationTimes,vtkPointData * seedData,int nVar,std::queue<vtkLagrangianParticle * > & particles)942 void vtkLagrangianParticleTracker::GenerateParticles(
943   const vtkBoundingBox* vtkNotUsed(bounds), vtkDataSet* seeds,
944   vtkDataArray* initialVelocities, vtkDataArray* initialIntegrationTimes,
945   vtkPointData* seedData, int nVar, std::queue<vtkLagrangianParticle*>& particles)
946 {
947   this->ParticleCounter = 0;
948   for (vtkIdType i = 0; i < seeds->GetNumberOfPoints(); i++)
949   {
950     double position[3];
951     seeds->GetPoint(i, position);
952     double initialIntegrationTime = initialIntegrationTimes ?
953       initialIntegrationTimes->GetTuple1(i) : 0;
954     vtkIdType particleId = this->GetNewParticleId();
955     vtkLagrangianParticle* particle = new vtkLagrangianParticle(nVar, particleId,
956       particleId, i, initialIntegrationTime, seedData);
957     memcpy(particle->GetPosition(), position, 3 * sizeof(double));
958     initialVelocities->GetTuple(i, particle->GetVelocity());
959     this->IntegrationModel->InitializeParticle(particle);
960     if (this->IntegrationModel->FindInLocators(particle->GetPosition()))
961     {
962       particles.push(particle);
963     }
964     else
965     {
966       delete particle;
967     }
968   }
969 }
970 
971 //---------------------------------------------------------------------------
GetParticleFeed(std::queue<vtkLagrangianParticle * > & vtkNotUsed (particleQueue))972 void vtkLagrangianParticleTracker::GetParticleFeed(
973   std::queue<vtkLagrangianParticle*>& vtkNotUsed(particleQueue))
974 {
975 }
976 
977 //---------------------------------------------------------------------------
Integrate(vtkLagrangianParticle * particle,std::queue<vtkLagrangianParticle * > & particlesQueue,vtkPolyData * particlePathsOutput,vtkIdList * particlePathPointId,vtkDataObject * interactionOutput)978 int vtkLagrangianParticleTracker::Integrate(vtkLagrangianParticle* particle,
979   std::queue<vtkLagrangianParticle*>& particlesQueue,
980   vtkPolyData* particlePathsOutput, vtkIdList* particlePathPointId,
981   vtkDataObject* interactionOutput)
982 {
983   // Sanity check
984   if (particle == nullptr)
985   {
986     vtkErrorMacro(<< "Cannot integrate nullptr particle");
987     return -1;
988   }
989 
990   // Set the current particle
991   this->IntegrationModel->SetCurrentParticle(particle);
992 
993   // Integrate until MaximumNumberOfSteps or MaximumIntegrationTime is reached or special case stops integration
994   int integrationRes = 0;
995   double stepFactor = this->StepFactor;
996   double reintegrationFactor = 1;
997   double& stepTimeActual = particle->GetStepTimeRef();
998   while (particle->GetTermination() ==
999          vtkLagrangianParticle::PARTICLE_TERMINATION_NOT_TERMINATED)
1000   {
1001     // Update progress
1002     if (particle->GetNumberOfSteps() % 100 == 0 && this->ParticleCounter > 0)
1003     {
1004       double progress = 1.0;
1005       if (this->MaximumNumberOfSteps != -1)
1006       {
1007         progress = static_cast<double>(particle->GetId() +
1008           static_cast<double>(particle->GetNumberOfSteps()) / this->MaximumNumberOfSteps) /
1009           this->ParticleCounter;
1010       }
1011       else if (this->MaximumIntegrationTime >= 0.0)
1012       {
1013         progress = static_cast<double>(particle->GetId() +
1014           particle->GetIntegrationTime() / this->MaximumIntegrationTime) /
1015           this->ParticleCounter;
1016       }
1017       this->UpdateProgress(progress);
1018       if (this->GetAbortExecute())
1019       {
1020         break;
1021       }
1022     }
1023 
1024     // Compute step
1025     double velocityMagnitude = reintegrationFactor * std::max(
1026       this->MinimumVelocityMagnitude,
1027       vtkMath::Norm(particle->GetVelocity()));
1028     double cellLength = this->ComputeCellLength(particle);
1029 
1030     double stepLength    = stepFactor          * cellLength;
1031     double stepLengthMin = this->StepFactorMin * cellLength;
1032     double stepLengthMax = this->StepFactorMax * cellLength;
1033     double stepTime    = stepLength    / (reintegrationFactor * velocityMagnitude);
1034     double stepTimeMin = stepLengthMin / (reintegrationFactor * velocityMagnitude);
1035     double stepTimeMax = stepLengthMax / (reintegrationFactor * velocityMagnitude);
1036 
1037     // Integrate one step
1038     if (!this->ComputeNextStep(particle->GetEquationVariables(),
1039       particle->GetNextEquationVariables(), particle->GetIntegrationTime(),
1040       stepTime, stepTimeActual, stepTimeMin, stepTimeMax, integrationRes))
1041     {
1042       vtkErrorMacro(<< "Integration Error");
1043       break;
1044     }
1045 
1046     bool stagnating =
1047       std::abs(particle->GetPosition()[0] - particle->GetNextPosition()[0]) <
1048       std::numeric_limits<double>::epsilon() &&
1049       std::abs(particle->GetPosition()[1] - particle->GetNextPosition()[1]) <
1050       std::numeric_limits<double>::epsilon() &&
1051       std::abs(particle->GetPosition()[2] - particle->GetNextPosition()[2]) <
1052       std::numeric_limits<double>::epsilon();
1053 
1054     // Only stagnating OUT_OF_DOMAIN are actually out of domain
1055     bool outOfDomain = integrationRes ==
1056       vtkInitialValueProblemSolver::OUT_OF_DOMAIN && stagnating;
1057 
1058     // Simpler Adaptive Step Reintegration code
1059     if (this->AdaptiveStepReintegration &&
1060         this->IntegrationModel->CheckAdaptiveStepReintegration(particle))
1061     {
1062       double stepLengthCurr2 = vtkMath::Distance2BetweenPoints(
1063         particle->GetPosition(), particle->GetNextPosition());
1064       double stepLengthMax2 = stepLengthMax * stepLengthMax;
1065       if (stepLengthCurr2 > stepLengthMax2)
1066       {
1067         reintegrationFactor *= 2;
1068         continue;
1069       }
1070       reintegrationFactor = 1;
1071     }
1072 
1073     if (outOfDomain)
1074     {
1075       // Stop integration
1076       particle->SetTermination(vtkLagrangianParticle::PARTICLE_TERMINATION_OUT_OF_DOMAIN);
1077       break;
1078     }
1079 
1080     // We care only about non-stagnating particle
1081     if (!stagnating)
1082     {
1083       // Surface interaction
1084       vtkLagrangianBasicIntegrationModel::PassThroughParticlesType passThroughParticles;
1085       unsigned int interactedSurfaceFlaxIndex;
1086       vtkLagrangianParticle* interactionParticle =
1087         this->IntegrationModel->ComputeSurfaceInteraction(
1088         particle, particlesQueue, interactedSurfaceFlaxIndex, passThroughParticles);
1089       if (interactionParticle != nullptr)
1090       {
1091         this->InsertInteractionOutputPoint(interactionParticle,
1092           interactedSurfaceFlaxIndex, interactionOutput);
1093         delete interactionParticle;
1094         interactionParticle = nullptr;
1095       }
1096 
1097       // Insert pass through interaction points
1098       // Note: when going out of domain right after going some pass through
1099       // surfaces, the pass through interaction point will not be
1100       // on a particle track, since we do not want to show out of domain particle
1101       // track. The pass through interaction still has occurred and it is not a bug.
1102       while (!passThroughParticles.empty())
1103       {
1104         vtkLagrangianBasicIntegrationModel::PassThroughParticlesItem item =
1105           passThroughParticles.front();
1106         passThroughParticles.pop();
1107         this->InsertInteractionOutputPoint(item.second, item.first, interactionOutput);
1108 
1109         // the pass through particles needs to be deleted
1110         delete item.second;
1111       }
1112 
1113       // Particle has been correctly integrated and interacted, record it
1114       // Insert Current particle as an output point
1115       this->InsertPathOutputPoint(particle, particlePathsOutput, particlePathPointId);
1116 
1117       // Particle has been terminated by surface
1118       if (particle->GetTermination() !=
1119         vtkLagrangianParticle::PARTICLE_TERMINATION_NOT_TERMINATED)
1120       {
1121         // Insert last particle path point on surface
1122         particle->MoveToNextPosition();
1123         this->InsertPathOutputPoint(particle, particlePathsOutput, particlePathPointId);
1124 
1125         // stop integration
1126         break;
1127       }
1128     }
1129 
1130     if (this->IntegrationModel->CheckFreeFlightTermination(particle))
1131     {
1132       particle->SetTermination(
1133         vtkLagrangianParticle::PARTICLE_TERMINATION_FLIGHT_TERMINATED);
1134       break;
1135     }
1136 
1137     // Keep integrating
1138     particle->MoveToNextPosition();
1139 
1140     // Compute now adaptive step
1141     if (this->Integrator->IsAdaptive() || this->AdaptiveStepReintegration)
1142     {
1143       stepFactor = stepTime * reintegrationFactor * velocityMagnitude / cellLength;
1144     }
1145     if (this->MaximumNumberOfSteps > -1 &&
1146         particle->GetNumberOfSteps() == this->MaximumNumberOfSteps &&
1147         particle->GetTermination() ==
1148         vtkLagrangianParticle::PARTICLE_TERMINATION_NOT_TERMINATED)
1149     {
1150       particle->SetTermination(
1151         vtkLagrangianParticle::PARTICLE_TERMINATION_OUT_OF_STEPS);
1152     }
1153     if (this->MaximumIntegrationTime >= 0.0 &&
1154         particle->GetIntegrationTime() >= this->MaximumIntegrationTime &&
1155         particle->GetTermination() ==
1156         vtkLagrangianParticle::PARTICLE_TERMINATION_NOT_TERMINATED)
1157     {
1158       particle->SetTermination(
1159         vtkLagrangianParticle::PARTICLE_TERMINATION_OUT_OF_TIME);
1160     }
1161   }
1162 
1163   this->IntegrationModel->SetCurrentParticle(nullptr);
1164   return integrationRes;
1165 }
1166 
1167 //---------------------------------------------------------------------------
InsertPathOutputPoint(vtkLagrangianParticle * particle,vtkPolyData * particlePathsOutput,vtkIdList * particlePathPointId,bool prev)1168 void vtkLagrangianParticleTracker::InsertPathOutputPoint(
1169   vtkLagrangianParticle* particle, vtkPolyData* particlePathsOutput,
1170   vtkIdList* particlePathPointId, bool prev)
1171 {
1172   // Recover structures
1173   vtkPoints* particlePathsPoints = particlePathsOutput->GetPoints();
1174   vtkPointData* particlePathsPointData = particlePathsOutput->GetPointData();
1175 
1176   // Store previous position
1177   vtkIdType pointId = particlePathsPoints->InsertNextPoint(
1178     prev ? particle->GetPrevPosition() : particle->GetPosition());
1179 
1180   particlePathPointId->InsertNextId(pointId);
1181 
1182   // Insert particle data
1183   this->InsertParticleData(particle, particlePathsPointData,
1184     prev ? vtkLagrangianBasicIntegrationModel::VARIABLE_STEP_PREV :
1185     vtkLagrangianBasicIntegrationModel::VARIABLE_STEP_CURRENT);
1186 
1187   // Add Variables data
1188   this->IntegrationModel->InsertVariablesParticleData(particle,
1189     particlePathsPointData, prev ?
1190     vtkLagrangianBasicIntegrationModel::VARIABLE_STEP_PREV :
1191     vtkLagrangianBasicIntegrationModel::VARIABLE_STEP_CURRENT);
1192 }
1193 
1194 //---------------------------------------------------------------------------
InsertInteractionOutputPoint(vtkLagrangianParticle * particle,unsigned int interactedSurfaceFlatIndex,vtkDataObject * interactionOutput)1195 void vtkLagrangianParticleTracker::InsertInteractionOutputPoint(
1196   vtkLagrangianParticle* particle, unsigned int interactedSurfaceFlatIndex,
1197   vtkDataObject* interactionOutput)
1198 {
1199   // Find the correct output
1200   vtkCompositeDataSet *hdOutput = vtkCompositeDataSet::SafeDownCast(interactionOutput);
1201   vtkPolyData* pdOutput = vtkPolyData::SafeDownCast(interactionOutput);
1202   vtkPolyData* interactionPd = nullptr;
1203   if (hdOutput)
1204   {
1205     vtkSmartPointer<vtkCompositeDataIterator> iter;
1206     iter.TakeReference(hdOutput->NewIterator());
1207     for (iter->InitTraversal(); !iter->IsDoneWithTraversal(); iter->GoToNextItem())
1208     {
1209       if (interactedSurfaceFlatIndex == iter->GetCurrentFlatIndex())
1210       {
1211         interactionPd = vtkPolyData::SafeDownCast(hdOutput->GetDataSet(iter));
1212         break;
1213       }
1214     }
1215   }
1216   else if (pdOutput)
1217   {
1218     interactionPd = pdOutput;
1219   }
1220 
1221   if (!interactionPd)
1222   {
1223     vtkErrorMacro(<< "Something went wrong with interaction output, "
1224       "cannot find correct interaction output polydata");
1225     return;
1226   }
1227 
1228   // "Next" Point
1229   vtkPoints* points = interactionPd->GetPoints();
1230   points->InsertNextPoint(particle->GetNextPosition());
1231 
1232   // Fill up interaction point data
1233   vtkPointData* pointData = interactionPd->GetPointData();
1234   this->InsertPathData(particle, pointData);
1235   this->IntegrationModel->InsertModelPathData(particle, pointData);
1236   this->InsertInteractionData(particle, pointData);
1237   this->InsertParticleData(particle, pointData,
1238     vtkLagrangianBasicIntegrationModel::VARIABLE_STEP_NEXT);
1239 
1240   // Add Variables data
1241   this->IntegrationModel->InsertVariablesParticleData(particle, pointData,
1242     vtkLagrangianBasicIntegrationModel::VARIABLE_STEP_NEXT);
1243 
1244   // Finally, Insert data from seed data only on not yet written arrays
1245   this->InsertSeedData(particle, pointData);
1246 }
1247 
1248 //---------------------------------------------------------------------------
InsertSeedData(vtkLagrangianParticle * particle,vtkFieldData * data)1249 void vtkLagrangianParticleTracker::InsertSeedData(vtkLagrangianParticle* particle,
1250   vtkFieldData* data)
1251 {
1252   // Check for max number of tuples in arrays
1253   vtkIdType maxTuples = 0;
1254   for (int i = 0; i < data->GetNumberOfArrays(); i++)
1255   {
1256     maxTuples = std::max(data->GetArray(i)->GetNumberOfTuples(), maxTuples);
1257   }
1258 
1259   // Copy seed data in not yet written array only
1260   // ie not yet at maxTuple
1261   vtkPointData* seedData = particle->GetSeedData();
1262   for (int i = 0; i < seedData->GetNumberOfArrays(); i++)
1263   {
1264     const char* name = seedData->GetArrayName(i);
1265     vtkDataArray* arr = data->GetArray(name);
1266     if (arr->GetNumberOfTuples() < maxTuples)
1267     {
1268       arr->InsertNextTuple(
1269         seedData->GetArray(i)->GetTuple(particle->GetSeedArrayTupleIndex()));
1270     }
1271   }
1272   // here all arrays from data should have the exact same size
1273 }
1274 
1275 //---------------------------------------------------------------------------
InsertPathData(vtkLagrangianParticle * particle,vtkFieldData * data)1276 void vtkLagrangianParticleTracker::InsertPathData(vtkLagrangianParticle* particle,
1277   vtkFieldData* data)
1278 {
1279   vtkLongLongArray::SafeDownCast(
1280     data->GetArray("Id"))->InsertNextValue(particle->GetId());
1281   vtkLongLongArray::SafeDownCast(
1282     data->GetArray("ParentId"))
1283     ->InsertNextValue(particle->GetParentId());
1284   vtkLongLongArray::SafeDownCast(
1285     data->GetArray("SeedId"))->InsertNextValue(particle->GetSeedId());
1286   vtkIntArray::SafeDownCast(
1287     data->GetArray("Termination"))->InsertNextValue(particle->GetTermination());
1288 }
1289 
1290 //---------------------------------------------------------------------------
InsertInteractionData(vtkLagrangianParticle * particle,vtkFieldData * data)1291 void vtkLagrangianParticleTracker::InsertInteractionData(
1292   vtkLagrangianParticle* particle, vtkFieldData* data)
1293 {
1294   vtkIntArray::SafeDownCast(
1295     data->GetArray("Interaction"))->InsertNextValue(particle->GetInteraction());
1296 }
1297 
1298 //---------------------------------------------------------------------------
InsertParticleData(vtkLagrangianParticle * particle,vtkFieldData * data,int stepEnum)1299 void vtkLagrangianParticleTracker::InsertParticleData(vtkLagrangianParticle* particle,
1300   vtkFieldData* data, int stepEnum)
1301 {
1302   switch (stepEnum)
1303   {
1304     case vtkLagrangianBasicIntegrationModel::VARIABLE_STEP_PREV:
1305       vtkIntArray::SafeDownCast(
1306         data->GetArray("StepNumber"))->InsertNextValue(particle->GetNumberOfSteps() - 1);
1307       data->GetArray("ParticleVelocity")->InsertNextTuple(particle->GetPrevVelocity());
1308       data->GetArray("IntegrationTime")->InsertNextTuple1(particle->GetPrevIntegrationTime());
1309       break;
1310     case vtkLagrangianBasicIntegrationModel::VARIABLE_STEP_CURRENT:
1311       vtkIntArray::SafeDownCast(
1312         data->GetArray("StepNumber"))->InsertNextValue(particle->GetNumberOfSteps());
1313       data->GetArray("ParticleVelocity")->InsertNextTuple(particle->GetVelocity());
1314       data->GetArray("IntegrationTime")->InsertNextTuple1(particle->GetIntegrationTime());
1315       break;
1316     case vtkLagrangianBasicIntegrationModel::VARIABLE_STEP_NEXT:
1317       vtkIntArray::SafeDownCast(
1318         data->GetArray("StepNumber"))->InsertNextValue(particle->GetNumberOfSteps() + 1);
1319       data->GetArray("ParticleVelocity")->InsertNextTuple(particle->GetNextVelocity());
1320       data->GetArray("IntegrationTime")->InsertNextTuple1(particle->GetIntegrationTime() +
1321         particle->GetStepTimeRef());
1322       break;
1323     default:
1324       break;
1325   }
1326 }
1327 
1328 //---------------------------------------------------------------------------
ComputeCellLength(vtkLagrangianParticle * particle)1329 double vtkLagrangianParticleTracker::ComputeCellLength(
1330   vtkLagrangianParticle* particle)
1331 {
1332   double cellLength = 1.0;
1333   vtkDataSet* dataset = nullptr;
1334   vtkCell* cell = nullptr;
1335   bool forceLastCell = false;
1336   if (this->CellLengthComputationMode == STEP_CUR_CELL_LENGTH ||
1337     this->CellLengthComputationMode == STEP_CUR_CELL_VEL_DIR ||
1338     this->CellLengthComputationMode == STEP_CUR_CELL_DIV_THEO)
1339   {
1340     vtkIdType cellId;
1341     if (this->IntegrationModel->FindInLocators(particle->GetPosition(), dataset, cellId))
1342     {
1343       cell = dataset->GetCell(cellId);
1344     }
1345     else
1346     {
1347       forceLastCell = true;
1348     }
1349   }
1350   if (this->CellLengthComputationMode == STEP_LAST_CELL_LENGTH ||
1351     this->CellLengthComputationMode == STEP_LAST_CELL_VEL_DIR ||
1352     this->CellLengthComputationMode == STEP_LAST_CELL_DIV_THEO ||
1353     forceLastCell)
1354   {
1355     dataset = particle->GetLastDataSet();
1356     if (!dataset)
1357     {
1358       return cellLength;
1359     }
1360     cell = dataset->GetCell(particle->GetLastCellId());
1361     if (!cell)
1362     {
1363       return cellLength;
1364     }
1365   }
1366   if (cell == nullptr)
1367   {
1368     vtkWarningMacro("Unsupported Cell Length Computation Mode"
1369       " or could not find a cell to compute cell length with");
1370     return 1.0;
1371   }
1372 
1373   double* vel = particle->GetVelocity();
1374   if ((this->CellLengthComputationMode == STEP_CUR_CELL_VEL_DIR ||
1375     this->CellLengthComputationMode == STEP_LAST_CELL_VEL_DIR) &&
1376       vtkMath::Norm(vel) > 0.0)
1377   {
1378     double velHat[3] = { vel[0], vel[1], vel[2] };
1379     vtkMath::Normalize(velHat);
1380     double tmpCellLength = 0.0;
1381     for (int ne = 0; ne < cell->GetNumberOfEdges(); ++ne)
1382     {
1383       double evect[3], x0[3], x1[3];
1384       vtkCell* edge = cell->GetEdge(ne);
1385       vtkIdType e0 = edge->GetPointId(0);
1386       vtkIdType e1 = edge->GetPointId(1);
1387       dataset->GetPoint(e0, x0);
1388       dataset->GetPoint(e1, x1);
1389       vtkMath::Subtract(x0, x1, evect);
1390       double elength = std::fabs(vtkMath::Dot(evect, velHat));
1391       tmpCellLength = std::max(tmpCellLength, elength);
1392     }
1393     cellLength = tmpCellLength;
1394   }
1395   else if ((this->CellLengthComputationMode == STEP_CUR_CELL_DIV_THEO ||
1396     this->CellLengthComputationMode == STEP_LAST_CELL_DIV_THEO) &&
1397       vtkMath::Norm(vel) > 0.0 && !vtkVoxel::SafeDownCast(cell))
1398   {
1399     double velHat[3] = {vel[0], vel[1], vel[2]};
1400     vtkMath::Normalize(velHat);
1401     double xa = 0.0;  // cell cross-sectional area in velHat direction
1402     double vol = 0.0; // cell volume
1403     for (int nf = 0; nf < cell->GetNumberOfFaces(); ++nf)
1404     {
1405       double norm[3];  // cell face normal
1406       double centroid[3] = {0.0, 0.0, 0.0}; // cell face centroid
1407       vtkCell* face = cell->GetFace(nf);
1408       vtkPoints* pts = face->GetPoints();
1409       vtkIdType nPoints = pts->GetNumberOfPoints();
1410       const double area = vtkPolygon::ComputeArea(pts, nPoints, nullptr, norm);
1411       const double fact = 1.0 / static_cast<double>(nPoints);
1412       for (int np = 0; np < nPoints; ++np)
1413       {
1414         double* x = pts->GetPoint(np);
1415         for (int nc = 0; nc < 3; ++nc)
1416         {
1417           centroid[nc] += x[nc] * fact;
1418         }
1419       }
1420       xa += std::fabs(vtkMath::Dot(norm, velHat) * area) / 2.0;   // sum unsigned areas
1421       vol += vtkMath::Dot(norm, centroid) * area / 3.0;           // using divergence theorem
1422     }
1423     // characteristic length is cell volume / cell cross-sectional area in velocity direction
1424     // Absolute value of volume because of some Fluent cases where all the volumes seem negative
1425     cellLength = std::fabs(vol) / xa;
1426   }
1427   else
1428   {
1429     cellLength = std::sqrt(cell->GetLength2());
1430   }
1431   return cellLength;
1432 }
1433 
1434 //---------------------------------------------------------------------------
ComputeNextStep(double * xprev,double * xnext,double t,double & delT,double & delTActual,double minStep,double maxStep,int & integrationRes)1435 bool vtkLagrangianParticleTracker::ComputeNextStep(
1436   double* xprev, double* xnext,
1437   double t, double& delT, double& delTActual,
1438   double minStep, double maxStep,
1439   int& integrationRes)
1440 {
1441   // Check for potential manual integration
1442   double error;
1443   if (!this->IntegrationModel->ManualIntegration(xprev, xnext, t, delT, delTActual,
1444     minStep, maxStep, this->IntegrationModel->GetTolerance(), error, integrationRes))
1445   {
1446     // integrate one step
1447     integrationRes =
1448       this->Integrator->ComputeNextStep(xprev, xnext, t, delT, delTActual,
1449         minStep, maxStep, this->IntegrationModel->GetTolerance(), error);
1450   }
1451 
1452   // Check failure cases
1453   if (integrationRes == vtkInitialValueProblemSolver::NOT_INITIALIZED)
1454   {
1455     vtkErrorMacro(<< "Integrator is not initialized. Aborting.");
1456     return false;
1457   }
1458   if (integrationRes == vtkInitialValueProblemSolver::UNEXPECTED_VALUE)
1459   {
1460     vtkErrorMacro(<< "Integrator encountered an unexpected value. Dropping particle.");
1461     return false;
1462   }
1463   return true;
1464 }
1465