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