1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkTemporalStreamTracer.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
16 // Hide VTK_DEPRECATED_IN_9_0_0() warnings for this class.
17 #define VTK_DEPRECATION_LEVEL 0
18
19 #include "vtkTemporalStreamTracer.h"
20
21 #include "vtkAbstractParticleWriter.h"
22 #include "vtkCellArray.h"
23 #include "vtkCellData.h"
24 #include "vtkCompositeDataIterator.h"
25 #include "vtkCompositeDataPipeline.h"
26 #include "vtkDataSetAttributes.h"
27 #include "vtkDoubleArray.h"
28 #include "vtkExecutive.h"
29 #include "vtkFloatArray.h"
30 #include "vtkGenericCell.h"
31 #include "vtkIdList.h"
32 #include "vtkInformation.h"
33 #include "vtkInformationVector.h"
34 #include "vtkIntArray.h"
35 #include "vtkMath.h"
36 #include "vtkMultiBlockDataSet.h"
37 #include "vtkObjectFactory.h"
38 #include "vtkOutputWindow.h"
39 #include "vtkPointData.h"
40 #include "vtkPointSet.h"
41 #include "vtkPolyData.h"
42 #include "vtkPolyLine.h"
43 #include "vtkRungeKutta2.h"
44 #include "vtkRungeKutta4.h"
45 #include "vtkRungeKutta45.h"
46 #include "vtkSignedCharArray.h"
47 #include "vtkSmartPointer.h"
48 #include "vtkTemporalInterpolatedVelocityField.h"
49 #include <cassert>
50
51 #ifdef _WIN32
52 #undef JB_H5PART_PARTICLE_OUTPUT
53 #else
54 // #define JB_H5PART_PARTICLE_OUTPUT
55 #endif
56
57 #ifdef JB_H5PART_PARTICLE_OUTPUT
58 // #include "vtkH5PartWriter.h"
59 #include "vtkXMLParticleWriter.h"
60 #endif
61
62 #include <algorithm>
63 #include <functional>
64
65 using namespace vtkTemporalStreamTracerNamespace;
66
67 // The 3D cell with the maximum number of points is VTK_LAGRANGE_HEXAHEDRON.
68 // We support up to 6th order hexahedra.
69 #define VTK_MAXIMUM_NUMBER_OF_POINTS 216
70
71 //------------------------------------------------------------------------------
72 //------------------------------------------------------------------------------
73 //#define JB_DEBUG__
74 #if defined JB_DEBUG__
75 #ifdef _WIN32
76 #define OUTPUTTEXT(a) vtkOutputWindowDisplayText(a);
77 #else
78 #endif
79
80 #undef vtkDebugMacro
81 #define vtkDebugMacro(a) \
82 { \
83 vtkOStreamWrapper::EndlType endl; \
84 vtkOStreamWrapper::UseEndl(endl); \
85 vtkOStrStreamWrapper vtkmsg; \
86 vtkmsg << "P(" << this->UpdatePieceId << "): " a << "\n"; \
87 OUTPUTTEXT(vtkmsg.str()); \
88 vtkmsg.rdbuf()->freeze(0); \
89 }
90
91 #undef vtkErrorMacro
92 #define vtkErrorMacro(a) vtkDebugMacro(a)
93
94 #endif
95 //------------------------------------------------------------------------------
96 vtkStandardNewMacro(vtkTemporalStreamTracer);
97 vtkCxxSetObjectMacro(vtkTemporalStreamTracer, ParticleWriter, vtkAbstractParticleWriter);
98 //------------------------------------------------------------------------------
vtkTemporalStreamTracer()99 vtkTemporalStreamTracer::vtkTemporalStreamTracer()
100 {
101 this->IntegrationDirection = FORWARD;
102 this->TimeStep = 0;
103 this->ActualTimeStep = 0;
104 this->NumberOfInputTimeSteps = 0;
105 this->ForceReinjectionEveryNSteps = 1;
106 this->ReinjectionFlag = false;
107 this->ReinjectionCounter = 0;
108 this->UpdatePieceId = 0;
109 this->UpdateNumPieces = 0;
110 this->AllFixedGeometry = 1;
111 this->StaticMesh = 1;
112 this->StaticSeeds = 1;
113 this->ComputeVorticity = true;
114 this->IgnorePipelineTime = 0;
115 this->ParticleWriter = nullptr;
116 this->ParticleFileName = nullptr;
117 this->EnableParticleWriting = false;
118 this->UniqueIdCounter = 0;
119 this->UniqueIdCounterMPI = 0;
120 this->InterpolationCount = 0;
121 //
122 this->NumberOfParticles = 0;
123 this->TimeStepResolution = 1.0;
124 this->TerminationTime = 0.0;
125 this->TerminationTimeUnit = TERMINATION_STEP_UNIT;
126 this->EarliestTime = -1E6;
127 // we are not actually using these for now
128
129 this->MaximumPropagation = 1.0;
130
131 this->IntegrationStepUnit = LENGTH_UNIT;
132 this->MinimumIntegrationStep = 1.0E-2;
133 this->MaximumIntegrationStep = 1.0;
134 this->InitialIntegrationStep = 0.5;
135 //
136 this->Interpolator = vtkSmartPointer<vtkTemporalInterpolatedVelocityField>::New();
137 //
138 this->SetNumberOfInputPorts(2);
139
140 #ifdef JB_H5PART_PARTICLE_OUTPUT
141 #ifdef _WIN32
142 vtkDebugMacro(<< "Setting vtkH5PartWriter");
143 vtkH5PartWriter* writer = vtkH5PartWriter::New();
144 #else
145 vtkDebugMacro(<< "Setting vtkXMLParticleWriter");
146 vtkXMLParticleWriter* writer = vtkXMLParticleWriter::New();
147 #endif
148 this->SetParticleWriter(writer);
149 writer->Delete();
150 #endif
151
152 this->SetIntegratorType(RUNGE_KUTTA4);
153 this->RequestIndex = 0;
154 VTK_LEGACY_BODY(vtkTemporalStreamTracer::vtkTemporalStreamTracer, "VTK 9.0");
155 }
156 //------------------------------------------------------------------------------
~vtkTemporalStreamTracer()157 vtkTemporalStreamTracer::~vtkTemporalStreamTracer()
158 {
159 this->SetParticleWriter(nullptr);
160 delete[] this->ParticleFileName;
161 this->ParticleFileName = nullptr;
162 }
163 //------------------------------------------------------------------------------
FillInputPortInformation(int port,vtkInformation * info)164 int vtkTemporalStreamTracer::FillInputPortInformation(int port, vtkInformation* info)
165 {
166 // port 0 must be a temporal collection of any type
167 // the executive should put a temporal collection in when
168 // we request multiple time steps.
169 if (port == 0)
170 {
171 info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataObject");
172 info->Set(vtkAlgorithm::INPUT_IS_REPEATABLE(), 1);
173 // info->Set(vtkAlgorithm::INPUT_IS_OPTIONAL(), 1);
174 }
175 else if (port == 1)
176 {
177 info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataSet");
178 info->Set(vtkAlgorithm::INPUT_IS_REPEATABLE(), 1);
179 // info->Set(vtkAlgorithm::INPUT_IS_OPTIONAL(), 1);
180 }
181 return 1;
182 }
183 //------------------------------------------------------------------------------
AddSourceConnection(vtkAlgorithmOutput * input)184 void vtkTemporalStreamTracer::AddSourceConnection(vtkAlgorithmOutput* input)
185 {
186 this->AddInputConnection(1, input);
187 }
188 //------------------------------------------------------------------------------
RemoveAllSources()189 void vtkTemporalStreamTracer::RemoveAllSources()
190 {
191 this->SetInputConnection(1, nullptr);
192 }
193 //------------------------------------------------------------------------------
ProcessRequest(vtkInformation * request,vtkInformationVector ** inputVector,vtkInformationVector * outputVector)194 vtkTypeBool vtkTemporalStreamTracer::ProcessRequest(
195 vtkInformation* request, vtkInformationVector** inputVector, vtkInformationVector* outputVector)
196 {
197 if (request->Has(vtkDemandDrivenPipeline::REQUEST_INFORMATION()))
198 {
199 return this->RequestInformation(request, inputVector, outputVector);
200 }
201 if (request->Has(vtkStreamingDemandDrivenPipeline::REQUEST_UPDATE_EXTENT()))
202 {
203 return this->RequestUpdateExtent(request, inputVector, outputVector);
204 }
205 if (request->Has(vtkDemandDrivenPipeline::REQUEST_DATA()))
206 {
207 return this->RequestData(request, inputVector, outputVector);
208 }
209 return this->Superclass::ProcessRequest(request, inputVector, outputVector);
210 }
211 //------------------------------------------------------------------------------
RequestInformation(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)212 int vtkTemporalStreamTracer::RequestInformation(vtkInformation* vtkNotUsed(request),
213 vtkInformationVector** inputVector, vtkInformationVector* outputVector)
214 {
215 vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
216 vtkInformation* outInfo = outputVector->GetInformationObject(0);
217
218 if (inInfo->Has(vtkStreamingDemandDrivenPipeline::TIME_STEPS()))
219 {
220 this->NumberOfInputTimeSteps = inInfo->Length(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
221 vtkDebugMacro(<< "vtkTemporalStreamTracer "
222 "inputVector TIME_STEPS "
223 << this->NumberOfInputTimeSteps);
224 //
225 // Get list of input time step values
226 this->InputTimeValues.resize(this->NumberOfInputTimeSteps);
227 inInfo->Get(vtkStreamingDemandDrivenPipeline::TIME_STEPS(), &this->InputTimeValues[0]);
228 if (this->NumberOfInputTimeSteps == 1)
229 {
230 vtkErrorMacro(<< "Not enough input time steps for particle integration");
231 return 0;
232 }
233 //
234 // We only output T-1 time steps
235 //
236 this->OutputTimeValues.resize(this->NumberOfInputTimeSteps - 1);
237 this->OutputTimeValues.clear();
238 this->OutputTimeValues.insert(this->OutputTimeValues.begin(), this->InputTimeValues.begin() + 1,
239 this->InputTimeValues.end());
240 }
241 else
242 {
243 vtkErrorMacro(<< "Input information has no TIME_STEPS set");
244 return 0;
245 }
246
247 outInfo->Set(vtkStreamingDemandDrivenPipeline::TIME_STEPS(), &this->OutputTimeValues[0],
248 static_cast<int>(this->OutputTimeValues.size()));
249
250 return 1;
251 }
252 //------------------------------------------------------------------------------
253 class WithinTolerance : public std::binary_function<double, double, bool>
254 {
255 public:
operator ()(first_argument_type a,second_argument_type b) const256 result_type operator()(first_argument_type a, second_argument_type b) const
257 {
258 bool result = (fabs(a - b) <= (a * 1E-6));
259 return (result_type)result;
260 }
261 };
262 //------------------------------------------------------------------------------
RequestUpdateExtent(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)263 int vtkTemporalStreamTracer::RequestUpdateExtent(vtkInformation* vtkNotUsed(request),
264 vtkInformationVector** inputVector, vtkInformationVector* outputVector)
265 {
266 int numInputs = inputVector[0]->GetNumberOfInformationObjects();
267 vtkInformation* outInfo = outputVector->GetInformationObject(0);
268
269 //
270 // The output has requested a time value, what times must we ask from our input
271 //
272 if (this->IgnorePipelineTime ||
273 !outInfo->Has(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP()))
274 {
275 double requestedTimeValue;
276 //
277 // ideally we want the output information to be requesting a time step,
278 // but since it isn't we must use the SetTimeStep value as a Time request
279 //
280 if (this->TimeStep < this->OutputTimeValues.size())
281 {
282 requestedTimeValue = this->OutputTimeValues[this->TimeStep];
283 }
284 else
285 {
286 requestedTimeValue = this->OutputTimeValues.back();
287 }
288 this->ActualTimeStep = this->TimeStep;
289
290 vtkDebugMacro(<< "SetTimeStep : requestedTimeValue " << requestedTimeValue
291 << " ActualTimeStep " << this->ActualTimeStep);
292 (void)requestedTimeValue; // silence unused variable warning
293 }
294 else
295 {
296 //
297 // Get the requested time step.
298 //
299 double requestedTimeValue = outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP());
300 this->ActualTimeStep =
301 std::find_if(this->OutputTimeValues.begin(), this->OutputTimeValues.end(),
302 std::bind(WithinTolerance(), std::placeholders::_1, requestedTimeValue)) -
303 this->OutputTimeValues.begin();
304 if (this->ActualTimeStep >= this->OutputTimeValues.size())
305 {
306 this->ActualTimeStep = 0;
307 }
308 vtkDebugMacro(<< "UPDATE_TIME_STEPS : requestedTimeValue " << requestedTimeValue
309 << " ActualTimeStep " << this->ActualTimeStep);
310 }
311
312 if (this->ActualTimeStep < this->OutputTimeValues.size())
313 {
314 for (int i = 0; i < numInputs; i++)
315 {
316 vtkInformation* inInfo = inputVector[0]->GetInformationObject(i);
317 // our output timestep T is timestep T+1 in the source
318 // so output inputTimeSteps[T], inputTimeSteps[T+1]
319 inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP(),
320 this->InputTimeValues[this->ActualTimeStep + this->RequestIndex]);
321 vtkDebugMacro(<< "requested 1 time values : "
322 << this->InputTimeValues[this->ActualTimeStep + this->RequestIndex]);
323 }
324 }
325 else
326 {
327 vtkDebugMacro(<< "UPDATE_TIME_STEPS : Error getting requested time step");
328 return 0;
329 }
330
331 return 1;
332 }
333 //------------------------------------------------------------------------------
InitializeInterpolator()334 int vtkTemporalStreamTracer::InitializeInterpolator()
335 {
336 if (!this->InputDataT[0] || !this->InputDataT[1])
337 {
338 return 0;
339 }
340 //
341 // When Multiblock arrays are processed, some may be empty
342 // if the first is empty, we won't find the correct vector name
343 // so scan until we get one
344 //
345 vtkSmartPointer<vtkCompositeDataIterator> iterP;
346 iterP.TakeReference(this->InputDataT[0]->NewIterator());
347 iterP->GoToFirstItem();
348 char* vecname = nullptr;
349 while (!iterP->IsDoneWithTraversal())
350 {
351 vtkDataArray* vectors = this->GetInputArrayToProcess(0, iterP->GetCurrentDataObject());
352 if (vectors)
353 {
354 vecname = vectors->GetName();
355 break;
356 }
357 iterP->GoToNextItem();
358 }
359 if (!vecname)
360 {
361 vtkDebugMacro(<< "Couldn't find vector array " << vecname);
362 return VTK_ERROR;
363 }
364
365 vtkDebugMacro(<< "Interpolator using array " << vecname);
366 this->Interpolator->SelectVectors(vecname);
367
368 this->AllFixedGeometry = 1;
369
370 int numValidInputBlocks[2] = { 0, 0 };
371 int numTotalInputBlocks[2] = { 0, 0 };
372 this->DataReferenceT[0] = this->DataReferenceT[1] = nullptr;
373 for (int T = 0; T < 2; T++)
374 {
375 this->CachedBounds[T].clear();
376 int index = 0;
377 // iterate over all blocks of input and cache the bounds information
378 // and determine fixed/dynamic mesh status.
379
380 vtkSmartPointer<vtkCompositeDataIterator> anotherIterP;
381 anotherIterP.TakeReference(this->InputDataT[T]->NewIterator());
382 anotherIterP->GoToFirstItem();
383 while (!anotherIterP->IsDoneWithTraversal())
384 {
385 numTotalInputBlocks[T]++;
386 vtkDataSet* inp = vtkDataSet::SafeDownCast(anotherIterP->GetCurrentDataObject());
387 if (inp)
388 {
389 if (inp->GetNumberOfCells() == 0)
390 {
391 vtkDebugMacro("Skipping an empty dataset");
392 }
393 else if (!inp->GetPointData()->GetVectors(vecname) && inp->GetNumberOfPoints() > 0)
394 {
395 vtkDebugMacro("One of the input datasets has no velocity vector.");
396 }
397 else
398 {
399 // vtkDebugMacro("pass " << i << " Found dataset with " << inp->GetNumberOfCells() << "
400 // cells");
401 //
402 // store the bounding boxes of each local dataset for faster 'point-in-dataset' testing
403 //
404 bounds bbox;
405 inp->ComputeBounds();
406 inp->GetBounds(&bbox.b[0]);
407 this->CachedBounds[T].push_back(bbox);
408 bool static_dataset = (this->StaticMesh != 0);
409 this->AllFixedGeometry = this->AllFixedGeometry && static_dataset;
410 // add the dataset to the interpolator
411 this->Interpolator->SetDataSetAtTime(
412 index++, T, this->CurrentTimeSteps[T], inp, static_dataset);
413 if (!this->DataReferenceT[T])
414 {
415 this->DataReferenceT[T] = inp;
416 }
417 //
418 numValidInputBlocks[T]++;
419 }
420 }
421 anotherIterP->GoToNextItem();
422 }
423 }
424 if (numValidInputBlocks[0] == 0 || numValidInputBlocks[1] == 0)
425 {
426 vtkDebugMacro("Not enough inputs have been found. Can not execute."
427 << numValidInputBlocks[0] << " " << numValidInputBlocks[1]);
428 return VTK_ERROR;
429 }
430 if (numValidInputBlocks[0] != numValidInputBlocks[1])
431 {
432 vtkDebugMacro("The number of datasets is different between time steps "
433 << numValidInputBlocks[0] << " " << numValidInputBlocks[1]);
434 return VTK_ERROR;
435 }
436 //
437 vtkDebugMacro("Number of Valid input blocks is " << numValidInputBlocks[0] << " from "
438 << numTotalInputBlocks[0]);
439 vtkDebugMacro("AllFixedGeometry " << this->AllFixedGeometry);
440
441 // force optimizations if StaticMesh is set.
442 if (this->StaticMesh)
443 {
444 vtkDebugMacro("Static Mesh optimizations Forced ON");
445 this->AllFixedGeometry = 1;
446 }
447
448 //
449 return VTK_OK;
450 }
SetTemporalInput(vtkDataObject * data,int i)451 int vtkTemporalStreamTracer::SetTemporalInput(vtkDataObject* data, int i)
452 {
453 // if not set, create a multiblock dataset to hold all input blocks
454 if (!this->InputDataT[i])
455 {
456 this->InputDataT[i] = vtkSmartPointer<vtkMultiBlockDataSet>::New();
457 }
458 // if simple dataset, add to our list, otherwise, add blocks
459 vtkDataSet* dsInput = vtkDataSet::SafeDownCast(data);
460 vtkMultiBlockDataSet* mbInput = vtkMultiBlockDataSet::SafeDownCast(data);
461
462 if (dsInput)
463 {
464 vtkSmartPointer<vtkDataSet> copy;
465 copy.TakeReference(dsInput->NewInstance());
466 copy->ShallowCopy(dsInput);
467 this->InputDataT[i]->SetBlock(this->InputDataT[i]->GetNumberOfBlocks(), copy);
468 }
469 else if (mbInput)
470 {
471 vtkSmartPointer<vtkCompositeDataIterator> iter;
472 iter.TakeReference(mbInput->NewIterator());
473 for (iter->InitTraversal(); !iter->IsDoneWithTraversal(); iter->GoToNextItem())
474 {
475 vtkDataSet* ds = vtkDataSet::SafeDownCast(iter->GetCurrentDataObject());
476 if (ds)
477 {
478 vtkSmartPointer<vtkDataSet> copy;
479 copy.TakeReference(ds->NewInstance());
480 copy->ShallowCopy(ds);
481 this->InputDataT[i]->SetBlock(this->InputDataT[i]->GetNumberOfBlocks(), copy);
482 }
483 }
484 }
485 else
486 {
487 vtkDebugMacro(
488 "This filter cannot handle inputs of type: " << (data ? data->GetClassName() : "(none)"));
489 return 0;
490 }
491
492 return 1;
493 }
494
495 //------------------------------------------------------------------------------
InsideBounds(double point[])496 bool vtkTemporalStreamTracer::InsideBounds(double point[])
497 {
498 double delta[3] = { 0.0, 0.0, 0.0 };
499 for (int t = 0; t < 2; ++t)
500 {
501 for (unsigned int i = 0; i < (this->CachedBounds[t].size()); ++i)
502 {
503 if (vtkMath::PointIsWithinBounds(point, &((this->CachedBounds[t])[i].b[0]), delta))
504 return true;
505 }
506 }
507 return false;
508 }
509 //------------------------------------------------------------------------------
TestParticles(ParticleVector & candidates,ParticleVector & passed,int & count)510 void vtkTemporalStreamTracer::TestParticles(
511 ParticleVector& candidates, ParticleVector& passed, int& count)
512 {
513 int i = 0;
514 int div = static_cast<int>(candidates.size() / 10);
515 if (div == 0)
516 {
517 div = 1;
518 }
519 count = 0;
520 for (ParticleIterator it = candidates.begin(); it != candidates.end(); ++it, ++i)
521 {
522 ParticleInformation& info = (*it);
523 double* pos = &info.CurrentPosition.x[0];
524 // if outside bounds, reject instantly
525 if (this->InsideBounds(pos))
526 {
527 if (info.UniqueParticleId == 602)
528 {
529 vtkDebugMacro(<< "TestParticles got 602");
530 }
531 // since this is first test, avoid bad cache tests
532 this->Interpolator->ClearCache();
533 info.LocationState = this->Interpolator->TestPoint(pos);
534 if (info.LocationState == ID_OUTSIDE_ALL /*|| location==ID_OUTSIDE_T0*/)
535 {
536 // can't really use this particle.
537 vtkDebugMacro(<< "TestParticles rejected particle");
538 }
539 else
540 {
541 // get the cached ids and datasets from the TestPoint call
542 this->Interpolator->GetCachedCellIds(info.CachedCellId, info.CachedDataSetId);
543 passed.push_back(info);
544 count++;
545 }
546 }
547 if (i % div == 0)
548 {
549 // vtkDebugMacro(<< "TestParticles " << i);
550 }
551 }
552 }
553 //------------------------------------------------------------------------------
AssignSeedsToProcessors(vtkDataSet * source,int sourceID,int ptId,ParticleVector & LocalSeedPoints,int & LocalAssignedCount)554 void vtkTemporalStreamTracer::AssignSeedsToProcessors(vtkDataSet* source, int sourceID, int ptId,
555 ParticleVector& LocalSeedPoints, int& LocalAssignedCount)
556 {
557 ParticleVector candidates;
558 //
559 // take points from the source object and create a particle list
560 //
561 vtkIdType numSeeds = source->GetNumberOfPoints();
562 candidates.resize(numSeeds);
563 //
564 for (vtkIdType i = 0; i < numSeeds; i++)
565 {
566 ParticleInformation& info = candidates[i];
567 memcpy(&(info.CurrentPosition.x[0]), source->GetPoint(i), sizeof(double) * 3);
568 info.CurrentPosition.x[3] = this->CurrentTimeSteps[0];
569 info.LocationState = 0;
570 info.CachedCellId[0] = -1;
571 info.CachedCellId[1] = -1;
572 info.CachedDataSetId[0] = 0;
573 info.CachedDataSetId[1] = 0;
574 info.SourceID = sourceID;
575 info.InjectedPointId = i + ptId;
576 info.InjectedStepId = this->ReinjectionCounter;
577 info.TimeStepAge = 0;
578 info.UniqueParticleId = -1;
579 info.rotation = 0.0;
580 info.angularVel = 0.0;
581 info.time = 0.0;
582 info.age = 0.0;
583 info.speed = 0.0;
584 info.ErrorCode = 0;
585 }
586 //
587 // Gather all Seeds to all processors for classification
588 //
589
590 #ifndef NDEBUG
591 int numTested = static_cast<int>(candidates.size());
592 #endif
593 this->TestParticles(candidates, LocalSeedPoints, LocalAssignedCount);
594 int TotalAssigned = LocalAssignedCount;
595 (void)TotalAssigned;
596
597 // Assign unique identifiers taking into account uneven distribution
598 // across processes and seeds which were rejected
599 this->AssignUniqueIds(LocalSeedPoints);
600
601 #ifndef NDEBUG
602 vtkDebugMacro(<< "Tested " << numTested << " LocallyAssigned " << LocalAssignedCount);
603 if (this->UpdatePieceId == 0)
604 {
605 vtkDebugMacro(<< "Total Assigned to all processes " << TotalAssigned);
606 }
607 #endif
608 }
609 //------------------------------------------------------------------------------
AssignUniqueIds(vtkTemporalStreamTracerNamespace::ParticleVector & LocalSeedPoints)610 void vtkTemporalStreamTracer::AssignUniqueIds(
611 vtkTemporalStreamTracerNamespace::ParticleVector& LocalSeedPoints)
612 {
613 vtkIdType ParticleCountOffset = 0;
614 vtkIdType numParticles = static_cast<vtkIdType>(LocalSeedPoints.size());
615 for (vtkIdType i = 0; i < numParticles; i++)
616 {
617 LocalSeedPoints[i].UniqueParticleId = this->UniqueIdCounter + ParticleCountOffset + i;
618 }
619 this->UniqueIdCounter += numParticles;
620 }
621 //------------------------------------------------------------------------------
TransmitReceiveParticles(ParticleVector &,ParticleVector &,bool)622 void vtkTemporalStreamTracer::TransmitReceiveParticles(ParticleVector&, ParticleVector&, bool) {}
623
624 //------------------------------------------------------------------------------
UpdateParticleList(ParticleVector & candidates)625 void vtkTemporalStreamTracer::UpdateParticleList(ParticleVector& candidates)
626 {
627 int numSeedsNew = static_cast<int>(candidates.size());
628 //
629 for (int i = 0; i < numSeedsNew; i++)
630 {
631 // allocate a new particle on the list and get a reference to it
632 this->ParticleHistories.push_back(candidates[i]);
633 }
634 this->NumberOfParticles = static_cast<int>(ParticleHistories.size());
635
636 vtkDebugMacro(<< "UpdateParticleList completed with " << this->NumberOfParticles << " particles");
637 }
638
ProcessInput(vtkInformationVector ** inputVector)639 int vtkTemporalStreamTracer::ProcessInput(vtkInformationVector** inputVector)
640 {
641 assert(this->RequestIndex >= 0 && this->RequestIndex < 2);
642 int numInputs = inputVector[0]->GetNumberOfInformationObjects();
643 if (numInputs != 1)
644 {
645 if (numInputs == 0)
646 {
647 vtkErrorMacro(<< "No input found.");
648 return 0;
649 }
650 vtkWarningMacro(<< "Multiple inputs founds. Use only the first one.");
651 }
652
653 // inherited from streamtracer, make sure it is null
654 this->InputData = nullptr;
655 this->InputDataT[this->RequestIndex] = nullptr;
656
657 vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
658 if (inInfo)
659 {
660 vtkDataObject* input = inInfo->Get(vtkDataObject::DATA_OBJECT());
661 SetTemporalInput(input, this->RequestIndex);
662 //
663 // Get the timestep information for this instant
664 //
665 std::vector<double> timesteps;
666 if (inInfo->Has(vtkDataObject::DATA_TIME_STEP()))
667 {
668 timesteps.resize(1);
669 timesteps[0] = inInfo->Get(vtkDataObject::DATA_TIME_STEP());
670 }
671 else
672 {
673 vtkErrorMacro(<< "No time step info");
674 return 1;
675 }
676 this->CurrentTimeSteps[this->RequestIndex] = timesteps[0] * this->TimeStepResolution;
677 }
678 return 1;
679 }
680
GenerateOutput(vtkInformationVector ** inputVector,vtkInformationVector * outputVector)681 int vtkTemporalStreamTracer::GenerateOutput(
682 vtkInformationVector** inputVector, vtkInformationVector* outputVector)
683
684 {
685
686 //
687 // Parallel/Piece information
688 //
689 vtkInformation* outInfo = outputVector->GetInformationObject(0);
690
691 this->UpdatePieceId = outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER());
692 this->UpdateNumPieces = outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES());
693
694 //
695 // How many Seed point sources are connected?
696 // Copy the sources into a vector for later use
697 //
698
699 int numSources = inputVector[1]->GetNumberOfInformationObjects();
700 std::vector<vtkDataSet*> SeedSources;
701 for (int idx = 0; idx < numSources; ++idx)
702 {
703 vtkDataObject* dobj = nullptr;
704 vtkInformation* inInfo = inputVector[1]->GetInformationObject(idx);
705 if (inInfo)
706 {
707 dobj = inInfo->Get(vtkDataObject::DATA_OBJECT());
708 SeedSources.push_back(vtkDataSet::SafeDownCast(dobj));
709 }
710 }
711
712 if (this->IntegrationDirection != FORWARD)
713 {
714 vtkErrorMacro(<< "We can only handle forward time particle tracking at the moment");
715 return 1;
716 }
717
718 //
719 // Add the datasets to an interpolator object
720 //
721 if (this->InitializeInterpolator() != VTK_OK)
722 {
723 if (this->InputDataT[0])
724 this->InputDataT[0] = nullptr;
725 if (this->InputDataT[1])
726 this->InputDataT[1] = nullptr;
727 vtkErrorMacro(<< "InitializeInterpolator failed");
728 return 1;
729 }
730
731 //
732 // Setup some variables
733 //
734 vtkSmartPointer<vtkInitialValueProblemSolver> integrator;
735 integrator.TakeReference(this->GetIntegrator()->NewInstance());
736 integrator->SetFunctionSet(this->Interpolator);
737
738 //
739 // Make sure the Particle Positions are initialized with Seed particles
740 //
741 this->ReinjectionFlag = false;
742 if (this->ForceReinjectionEveryNSteps > 0)
743 {
744 if ((this->ActualTimeStep % this->ForceReinjectionEveryNSteps) == 0)
745 {
746 this->ReinjectionFlag = true;
747 }
748 }
749 //
750 // If T=0 reset everything to allow us to setup stuff then start an animation
751 // with a clean slate
752 //
753 if (this->ActualTimeStep == 0) // XXX: what if I start from some other time?
754 {
755 this->LocalSeeds.clear();
756 this->ParticleHistories.clear();
757 this->EarliestTime = -1E6;
758 this->ReinjectionFlag = true;
759 this->ReinjectionCounter = 0;
760 this->UniqueIdCounter = 0;
761 this->UniqueIdCounterMPI = 0;
762 }
763 else if (this->CurrentTimeSteps[0] < this->EarliestTime)
764 {
765 //
766 // We don't want to go back in time, so just reuse whatever we have
767 //
768 vtkDebugMacro("skipping particle tracking because we have seen this timestep before");
769 outInfo->Set(vtkDataObject::DATA_TIME_STEP(), this->OutputTimeValues[this->ActualTimeStep]);
770 if (this->InputDataT[0])
771 this->InputDataT[0] = nullptr;
772 if (this->InputDataT[1])
773 this->InputDataT[1] = nullptr;
774 return 1;
775 }
776 this->EarliestTime = (this->CurrentTimeSteps[0] > this->EarliestTime) ? this->CurrentTimeSteps[0]
777 : this->EarliestTime;
778 //
779 //
780 //
781 for (unsigned int i = 0; i < SeedSources.size(); i++)
782 {
783 if (SeedSources[i]->GetMTime() > this->ParticleInjectionTime)
784 {
785 // this->ReinjectionFlag = 1;
786 }
787 }
788
789 //
790 // Lists for seed particles
791 //
792 ParticleVector candidates;
793 ParticleVector received;
794 //
795
796 if (this->ReinjectionFlag)
797 {
798 int seedPointId = 0;
799 if (this->StaticSeeds && this->AllFixedGeometry && this->LocalSeeds.empty())
800 {
801 for (unsigned int i = 0; i < SeedSources.size(); i++)
802 {
803 this->AssignSeedsToProcessors(SeedSources[i], i, 0, this->LocalSeeds, seedPointId);
804 }
805 }
806 else
807 {
808 // wipe the list and reclassify for each injection
809 this->LocalSeeds.clear();
810 for (unsigned int i = 0; i < SeedSources.size(); i++)
811 {
812 this->AssignSeedsToProcessors(SeedSources[i], i, 0, this->LocalSeeds, seedPointId);
813 }
814 }
815 this->ParticleInjectionTime.Modified();
816
817 // Now update our main list with the ones we are keeping
818 vtkDebugMacro(<< "Reinjection about to update candidates (" << this->LocalSeeds.size()
819 << " particles)");
820 this->UpdateParticleList(this->LocalSeeds);
821 this->ReinjectionCounter += 1;
822 }
823
824 //
825 // setup all our output arrays
826 //
827 vtkDebugMacro(<< "About to allocate point arrays ");
828 vtkPolyData* output = vtkPolyData::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));
829 this->ParticleAge = vtkSmartPointer<vtkFloatArray>::New();
830 this->ParticleIds = vtkSmartPointer<vtkIntArray>::New();
831 this->ParticleSourceIds = vtkSmartPointer<vtkSignedCharArray>::New();
832 this->InjectedPointIds = vtkSmartPointer<vtkIntArray>::New();
833 this->InjectedStepIds = vtkSmartPointer<vtkIntArray>::New();
834 this->ErrorCodeArray = vtkSmartPointer<vtkIntArray>::New();
835 this->ParticleVorticity = vtkSmartPointer<vtkFloatArray>::New();
836 this->ParticleRotation = vtkSmartPointer<vtkFloatArray>::New();
837 this->ParticleAngularVel = vtkSmartPointer<vtkFloatArray>::New();
838 this->cellVectors = vtkSmartPointer<vtkDoubleArray>::New();
839 this->ParticleCells = vtkSmartPointer<vtkCellArray>::New();
840 this->OutputCoordinates = vtkSmartPointer<vtkPoints>::New();
841 this->OutputPointData = output->GetPointData();
842 this->OutputPointData->Initialize();
843 this->InterpolationCount = 0;
844 vtkDebugMacro(<< "About to Interpolate allocate space");
845 this->OutputPointData->InterpolateAllocate(this->DataReferenceT[1]->GetPointData());
846 //
847 this->ParticleAge->SetName("ParticleAge");
848 this->ParticleIds->SetName("ParticleId");
849 this->ParticleSourceIds->SetName("ParticleSourceId");
850 this->InjectedPointIds->SetName("InjectedPointId");
851 this->InjectedStepIds->SetName("InjectionStepId");
852 this->ErrorCodeArray->SetName("ErrorCode");
853
854 if (this->ComputeVorticity)
855 {
856 this->cellVectors->SetNumberOfComponents(3);
857 this->cellVectors->Allocate(3 * VTK_CELL_SIZE);
858 this->ParticleVorticity->SetName("Vorticity");
859 this->ParticleRotation->SetName("Rotation");
860 this->ParticleAngularVel->SetName("AngularVelocity");
861 }
862
863 output->SetPoints(this->OutputCoordinates);
864 output->SetVerts(this->ParticleCells);
865 vtkDebugMacro(<< "Finished allocating point arrays ");
866
867 //
868 // Perform 2 passes
869 // Pass 0 : Integration of particles created by a source in this process
870 // or received at start from a source in another process.
871 //
872 // Pass 1 : Particles that were sent in mid integration from another process
873 // are added in and their integration continued here. In actual fact, the process
874 // should be repeated until all particles are finished, but the chances of
875 // a particle stepping inside and out again through a single domain
876 // in one time step are small (hopefully!)
877
878 vtkDebugMacro(<< "Clear MPI send list ");
879 this->MPISendList.clear();
880
881 #ifndef NDEBUG
882 int Number = static_cast<int>(this->ParticleHistories.size());
883 #endif
884 ParticleListIterator it_first = this->ParticleHistories.begin();
885 ParticleListIterator it_last = this->ParticleHistories.end();
886 ParticleListIterator it_next;
887 #define PASSES 2
888 for (int pass = 0; pass < PASSES; pass++)
889 {
890 vtkDebugMacro(<< "Begin Pass " << pass << " with " << Number << " Particles");
891 for (ParticleListIterator it = it_first; it != it_last;)
892 {
893 // Keep the 'next' iterator handy because if a particle is terminated
894 // or leaves the domain, the 'current' iterator will be deleted.
895 it_next = it;
896 it_next++;
897 //
898 // Shall we terminate this particle
899 //
900 double interval = (this->CurrentTimeSteps[1] - this->CurrentTimeSteps[0]);
901 bool terminated = false;
902 if (this->TerminationTime > 0)
903 {
904 if (this->TerminationTimeUnit == TERMINATION_TIME_UNIT &&
905 (it->age + interval) > this->TerminationTime)
906 {
907 terminated = true;
908 }
909 else if (this->TerminationTimeUnit == TERMINATION_STEP_UNIT &&
910 (it->TimeStepAge + 1) > this->TerminationTime)
911 {
912 terminated = true;
913 }
914 }
915 if (terminated)
916 {
917 this->ParticleHistories.erase(it);
918 }
919 else
920 {
921 this->IntegrateParticle(
922 it, this->CurrentTimeSteps[0], this->CurrentTimeSteps[1], integrator);
923 }
924 //
925 if (this->GetAbortExecute())
926 {
927 break;
928 }
929 it = it_next;
930 }
931 // Particles might have been deleted during the first pass as they move
932 // out of domain or age. Before adding any new particles that are sent
933 // to us, we must know the starting point ready for the second pass
934 bool list_valid = (!this->ParticleHistories.empty());
935 if (list_valid)
936 {
937 // point to one before the end
938 it_first = --this->ParticleHistories.end();
939 }
940 // Send and receive any particles which exited/entered the domain
941 if (this->UpdateNumPieces > 1 && pass < (PASSES - 1))
942 {
943 // the Particle lists will grow if any are received
944 // so we must be very careful with our iterators
945 vtkDebugMacro(<< "End of Pass " << pass << " with " << this->ParticleHistories.size() << " "
946 << " about to Transmit/Receive " << this->MPISendList.size());
947 this->TransmitReceiveParticles(this->MPISendList, received, true);
948 // don't want the ones that we sent away
949 this->MPISendList.clear();
950 int assigned;
951 // classify all the ones we received
952 if (!received.empty())
953 {
954 this->TestParticles(received, candidates, assigned);
955 vtkDebugMacro(<< "received " << received.size() << " : assigned locally " << assigned);
956 received.clear();
957 }
958 // Now update our main list with the ones we are keeping
959 this->UpdateParticleList(candidates);
960 // free up unwanted memory
961 #ifndef NDEBUG
962 Number = static_cast<int>(candidates.size());
963 #endif
964 candidates.clear();
965 }
966 it_last = this->ParticleHistories.end();
967 if (list_valid)
968 {
969 // increment to point to first new entry
970 it_first++;
971 }
972 else
973 {
974 it_first = this->ParticleHistories.begin();
975 }
976 }
977 if (!this->MPISendList.empty())
978 {
979 // If a particle went out of domain on the second pass, it should be sent
980 // can it really pass right through a domain in one step?
981 // what about grazing the edge of rotating zone?
982 vtkDebugMacro(<< "MPISendList not empty " << this->MPISendList.size());
983 this->MPISendList.clear();
984 }
985
986 //
987 // We must only add these scalar arrays at the end because the
988 // existing scalars on the input get interpolated during iteration
989 // over the particles
990 //
991 this->OutputPointData->AddArray(this->ParticleIds);
992 this->OutputPointData->AddArray(this->ParticleSourceIds);
993 this->OutputPointData->AddArray(this->InjectedPointIds);
994 this->OutputPointData->AddArray(this->InjectedStepIds);
995 this->OutputPointData->AddArray(this->ErrorCodeArray);
996 this->OutputPointData->AddArray(this->ParticleAge);
997 if (this->ComputeVorticity)
998 {
999 this->OutputPointData->AddArray(this->ParticleVorticity);
1000 this->OutputPointData->AddArray(this->ParticleRotation);
1001 this->OutputPointData->AddArray(this->ParticleAngularVel);
1002 }
1003
1004 if (this->InterpolationCount != this->OutputCoordinates->GetNumberOfPoints())
1005 {
1006 vtkErrorMacro(<< "Mismatch in point array/data counts");
1007 }
1008 //
1009 outInfo->Set(vtkDataObject::DATA_TIME_STEP(), this->OutputTimeValues[this->ActualTimeStep]);
1010
1011 // save some locator building, by re-using them as time progresses
1012 this->Interpolator->AdvanceOneTimeStep();
1013
1014 //
1015 // Let go of inputs
1016 //
1017 if (this->InputDataT[0])
1018 this->InputDataT[0] = nullptr;
1019 if (this->InputDataT[1])
1020 this->InputDataT[1] = nullptr;
1021
1022 //
1023 // Write Particles out if necessary
1024 //
1025 // NB. We don't want our writer to trigger any updates,
1026 // so shallow copy the output
1027 if (this->ParticleWriter && this->EnableParticleWriting)
1028 {
1029 vtkSmartPointer<vtkPolyData> polys = vtkSmartPointer<vtkPolyData>::New();
1030 polys->ShallowCopy(output);
1031 int N = polys->GetNumberOfPoints();
1032 (void)N;
1033 this->ParticleWriter->SetFileName(this->ParticleFileName);
1034 this->ParticleWriter->SetTimeStep(this->ActualTimeStep);
1035 this->ParticleWriter->SetTimeValue(this->CurrentTimeSteps[1]);
1036 this->ParticleWriter->SetInputData(polys);
1037 this->ParticleWriter->Write();
1038 this->ParticleWriter->CloseFile();
1039 this->ParticleWriter->SetInputData(nullptr);
1040 vtkDebugMacro(<< "Written " << N);
1041 }
1042 return 1;
1043 }
1044
RequestData(vtkInformation * request,vtkInformationVector ** inputVector,vtkInformationVector * outputVector)1045 int vtkTemporalStreamTracer::RequestData(
1046 vtkInformation* request, vtkInformationVector** inputVector, vtkInformationVector* outputVector)
1047 {
1048 //
1049 // Inputs information
1050 //
1051 bool result(true);
1052 if (RequestIndex < 2)
1053 {
1054 result = (this->ProcessInput(inputVector) == 1);
1055 if (result && RequestIndex == 1)
1056 {
1057 this->GenerateOutput(inputVector, outputVector);
1058 }
1059 }
1060
1061 this->RequestIndex++;
1062 if (result && this->RequestIndex < 2)
1063 {
1064 request->Set(vtkStreamingDemandDrivenPipeline::CONTINUE_EXECUTING(), 1);
1065 }
1066 else
1067 {
1068 request->Remove(vtkStreamingDemandDrivenPipeline::CONTINUE_EXECUTING());
1069 this->RequestIndex = 0;
1070 }
1071
1072 // this->Interpolator->ShowCacheResults();
1073 // vtkErrorMacro(<<"RequestData done");
1074 return 1;
1075 }
1076 //------------------------------------------------------------------------------
IntegrateParticle(ParticleListIterator & it,double currenttime,double targettime,vtkInitialValueProblemSolver * integrator)1077 void vtkTemporalStreamTracer::IntegrateParticle(ParticleListIterator& it, double currenttime,
1078 double targettime, vtkInitialValueProblemSolver* integrator)
1079 {
1080 double epsilon = (targettime - currenttime) / 100.0;
1081 double velocity[3], point1[4], point2[4] = { 0.0, 0.0, 0.0, 0.0 };
1082 double minStep = 0, maxStep = 0;
1083 double stepWanted, stepTaken = 0.0;
1084 substeps = 0;
1085
1086 ParticleInformation& info = (*it);
1087 // Get the Initial point {x,y,z,t}
1088 memcpy(point1, &info.CurrentPosition, sizeof(Position));
1089
1090 if (point1[3] < (currenttime - epsilon) || point1[3] > (targettime + epsilon))
1091 {
1092 vtkDebugMacro(<< "Bad particle time : expected (" << this->CurrentTimeSteps[0] << "-"
1093 << this->CurrentTimeSteps[1] << ") got " << point1[3]);
1094 }
1095
1096 IntervalInformation delT;
1097 delT.Unit = LENGTH_UNIT;
1098 delT.Interval = (targettime - currenttime) * this->InitialIntegrationStep;
1099 epsilon = delT.Interval * 1E-3;
1100
1101 //
1102 // begin interpolation between available time values, if the particle has
1103 // a cached cell ID and dataset - try to use it,
1104 //
1105 this->Interpolator->SetCachedCellIds(info.CachedCellId, info.CachedDataSetId);
1106
1107 bool particle_good = true;
1108 info.ErrorCode = 0;
1109 while (point1[3] < (targettime - epsilon))
1110 {
1111 //
1112 // Here beginneth the real work
1113 //
1114 double error = 0;
1115
1116 // If, with the next step, propagation will be larger than
1117 // max, reduce it so that it is (approximately) equal to max.
1118 stepWanted = delT.Interval;
1119 if ((point1[3] + stepWanted) > targettime)
1120 {
1121 stepWanted = targettime - point1[3];
1122 maxStep = stepWanted;
1123 }
1124 this->LastUsedStepSize = stepWanted;
1125
1126 // Calculate the next step using the integrator provided.
1127 // If the next point is out of bounds, send it to another process
1128 if (integrator->ComputeNextStep(point1, point2, point1[3], stepWanted, stepTaken, minStep,
1129 maxStep, this->MaximumError, error) != 0)
1130 {
1131 // if the particle is sent, remove it from the list
1132 info.ErrorCode = 1;
1133 if (this->SendParticleToAnotherProcess(info, point1, this->LastUsedStepSize))
1134 {
1135 this->ParticleHistories.erase(it);
1136 particle_good = false;
1137 break;
1138 }
1139 else
1140 {
1141 // particle was not sent, retry saved it, so copy info back
1142 substeps++;
1143 memcpy(point1, &info.CurrentPosition, sizeof(Position));
1144 }
1145 }
1146 else // success, increment position/time
1147 {
1148 substeps++;
1149
1150 // increment the particle time
1151 point2[3] = point1[3] + stepTaken;
1152 info.age += stepTaken;
1153
1154 // Point is valid. Insert it.
1155 memcpy(&info.CurrentPosition, point2, sizeof(Position));
1156 memcpy(point1, point2, sizeof(Position));
1157 }
1158
1159 // If the solver is adaptive and the next time step (delT.Interval)
1160 // that the solver wants to use is smaller than minStep or larger
1161 // than maxStep, re-adjust it. This has to be done every step
1162 // because minStep and maxStep can change depending on the Cell
1163 // size (unless it is specified in time units)
1164 if (integrator->IsAdaptive())
1165 {
1166 // code removed. Put it back when this is stable
1167 }
1168 }
1169 if (particle_good)
1170 {
1171 // The integration succeeded, but check the computed final position
1172 // is actually inside the domain (the intermediate steps taken inside
1173 // the integrator were ok, but the final step may just pass out)
1174 // if it moves out, we can't interpolate scalars, so we must send it away
1175 info.LocationState = this->Interpolator->TestPoint(info.CurrentPosition.x);
1176 if (info.LocationState == ID_OUTSIDE_ALL)
1177 {
1178 info.ErrorCode = 2;
1179 // if the particle is sent, remove it from the list
1180 if (this->SendParticleToAnotherProcess(info, point1, this->LastUsedStepSize))
1181 {
1182 this->ParticleHistories.erase(it);
1183 particle_good = false;
1184 }
1185 }
1186 }
1187
1188 //
1189 // Has this particle stagnated
1190 //
1191 if (particle_good)
1192 {
1193 this->Interpolator->GetLastGoodVelocity(velocity);
1194 info.speed = vtkMath::Norm(velocity);
1195 if (it->speed <= this->TerminalSpeed)
1196 {
1197 this->ParticleHistories.erase(it);
1198 particle_good = false;
1199 }
1200 }
1201
1202 //
1203 // We got this far without error :
1204 // Insert the point into the output
1205 // Create any new scalars and interpolate existing ones
1206 // Cache cell ids and datasets
1207 //
1208 if (particle_good)
1209 {
1210 //
1211 // store the last Cell Ids and dataset indices for next time particle is updated
1212 //
1213 this->Interpolator->GetCachedCellIds(info.CachedCellId, info.CachedDataSetId);
1214 //
1215 info.TimeStepAge += 1;
1216 //
1217 // Now generate the output geometry and scalars
1218 //
1219 double* coord = &info.CurrentPosition.x[0];
1220 vtkIdType tempId = this->OutputCoordinates->InsertNextPoint(coord);
1221 // create the cell
1222 this->ParticleCells->InsertNextCell(1, &tempId);
1223 // set the easy scalars for this particle
1224 this->ParticleIds->InsertNextValue(info.UniqueParticleId);
1225 this->ParticleSourceIds->InsertNextValue(info.SourceID);
1226 this->InjectedPointIds->InsertNextValue(info.InjectedPointId);
1227 this->InjectedStepIds->InsertNextValue(info.InjectedStepId);
1228 this->ErrorCodeArray->InsertNextValue(info.ErrorCode);
1229 this->ParticleAge->InsertNextValue(info.age);
1230 //
1231 // Interpolate all existing point attributes
1232 // In principle we always integrate the particle until it reaches Time2
1233 // - so we don't need to do any interpolation of the scalars
1234 // between T0 and T1, just fetch the values
1235 // of the spatially interpolated scalars from T1.
1236 //
1237 if (info.LocationState == ID_OUTSIDE_T1)
1238 {
1239 this->Interpolator->InterpolatePoint(0, this->OutputPointData, tempId);
1240 }
1241 else
1242 {
1243 this->Interpolator->InterpolatePoint(1, this->OutputPointData, tempId);
1244 }
1245 this->InterpolationCount++;
1246 //
1247 // Compute vorticity
1248 //
1249 if (this->ComputeVorticity)
1250 {
1251 vtkGenericCell* cell;
1252 double pcoords[3], vorticity[3], weights[VTK_MAXIMUM_NUMBER_OF_POINTS];
1253 double rotation, omega;
1254 // have to use T0 if particle is out at T1, otherwise use T1
1255 if (info.LocationState == ID_OUTSIDE_T1)
1256 {
1257 this->Interpolator->GetVorticityData(0, pcoords, weights, cell, this->cellVectors);
1258 }
1259 else
1260 {
1261 this->Interpolator->GetVorticityData(1, pcoords, weights, cell, this->cellVectors);
1262 }
1263 vtkStreamTracer::CalculateVorticity(cell, pcoords, cellVectors, vorticity);
1264 this->ParticleVorticity->InsertNextTuple(vorticity);
1265 // local rotation = vorticity . unit tangent ( i.e. velocity/speed )
1266 if (info.speed != 0.0)
1267 {
1268 omega = vtkMath::Dot(vorticity, velocity);
1269 omega /= info.speed;
1270 omega *= this->RotationScale;
1271 }
1272 else
1273 {
1274 omega = 0.0;
1275 }
1276 vtkIdType index = this->ParticleAngularVel->InsertNextValue(omega);
1277 if (index > 0)
1278 {
1279 rotation =
1280 info.rotation + (info.angularVel + omega) / 2 * (info.CurrentPosition.x[3] - info.time);
1281 }
1282 else
1283 {
1284 rotation = 0.0;
1285 }
1286 this->ParticleRotation->InsertNextValue(rotation);
1287 info.rotation = rotation;
1288 info.angularVel = omega;
1289 info.time = info.CurrentPosition.x[3];
1290 }
1291 }
1292 else
1293 this->Interpolator->ClearCache();
1294
1295 double eps = (this->CurrentTimeSteps[1] - this->CurrentTimeSteps[0]) / 100;
1296 if (point1[3] < (this->CurrentTimeSteps[0] - eps) ||
1297 point1[3] > (this->CurrentTimeSteps[1] + eps))
1298 {
1299 vtkDebugMacro(<< "Unexpected time ending IntegrateParticle - expected ("
1300 << this->CurrentTimeSteps[0] << "-" << this->CurrentTimeSteps[1] << ") got "
1301 << point1[3]);
1302 }
1303 }
1304 //------------------------------------------------------------------------------
RetryWithPush(ParticleInformation & info,double velocity[3],double delT)1305 bool vtkTemporalStreamTracer::RetryWithPush(
1306 ParticleInformation& info, double velocity[3], double delT)
1307 {
1308 // try adding a one increment push to the particle to get over a rotating/moving boundary
1309 for (int v = 0; v < 3; v++)
1310 info.CurrentPosition.x[v] += velocity[v] * delT;
1311 info.CurrentPosition.x[3] += delT;
1312 info.LocationState = this->Interpolator->TestPoint(info.CurrentPosition.x);
1313 if (info.LocationState != ID_OUTSIDE_ALL)
1314 {
1315 // a push helped the particle get back into a dataset,
1316 info.age += delT;
1317 info.ErrorCode = 6;
1318 return true;
1319 }
1320 return false;
1321 }
1322 //------------------------------------------------------------------------------
SendParticleToAnotherProcess(ParticleInformation & info,double point1[4],double delT)1323 bool vtkTemporalStreamTracer::SendParticleToAnotherProcess(
1324 ParticleInformation& info, double point1[4], double delT)
1325 {
1326 // return 1;
1327 double velocity[3];
1328 this->Interpolator->ClearCache();
1329 if (info.UniqueParticleId == 3)
1330 {
1331 vtkDebugMacro(<< "3 is about to be sent");
1332 }
1333 info.LocationState = this->Interpolator->TestPoint(point1);
1334 if (info.LocationState == ID_OUTSIDE_ALL)
1335 {
1336 // something is wrong, the particle has left the building completely
1337 // we can't get the last good velocity as it won't be valid
1338 // send the particle 'as is' and hope it lands in another process
1339 if (substeps > 0)
1340 {
1341 this->Interpolator->GetLastGoodVelocity(velocity);
1342 }
1343 else
1344 {
1345 velocity[0] = velocity[1] = velocity[2] = 0.0;
1346 }
1347 info.ErrorCode = 3;
1348 }
1349 else if (info.LocationState == ID_OUTSIDE_T0)
1350 {
1351 // the particle left the volume but can be tested at T2, so use the velocity at T2
1352 this->Interpolator->GetLastGoodVelocity(velocity);
1353 info.ErrorCode = 4;
1354 }
1355 else if (info.LocationState == ID_OUTSIDE_T1)
1356 {
1357 // the particle left the volume but can be tested at T1, so use the velocity at T1
1358 this->Interpolator->GetLastGoodVelocity(velocity);
1359 info.ErrorCode = 5;
1360 }
1361 else
1362 {
1363 // The test returned INSIDE_ALL, so test failed near start of integration,
1364 this->Interpolator->GetLastGoodVelocity(velocity);
1365 }
1366 if (this->RetryWithPush(info, velocity, delT))
1367 return false;
1368 this->AddParticleToMPISendList(info);
1369 return true;
1370 }
1371 //------------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)1372 void vtkTemporalStreamTracer::PrintSelf(ostream& os, vtkIndent indent)
1373 {
1374 this->Superclass::PrintSelf(os, indent);
1375
1376 os << indent << "TimeStepResolution: " << this->TimeStepResolution << endl;
1377 os << indent << "ParticleWriter: " << this->ParticleWriter << endl;
1378 os << indent << "ParticleFileName: " << (this->ParticleFileName ? this->ParticleFileName : "None")
1379 << endl;
1380 os << indent << "TimeStep: " << this->TimeStep << endl;
1381 os << indent << "ForceReinjectionEveryNSteps: " << this->ForceReinjectionEveryNSteps << endl;
1382 os << indent << "EnableParticleWriting: " << this->EnableParticleWriting << endl;
1383 os << indent << "IgnorePipelineTime: " << this->IgnorePipelineTime << endl;
1384 os << indent << "StaticMesh: " << this->StaticMesh << endl;
1385 os << indent << "TerminationTime: " << this->TerminationTime << endl;
1386 os << indent << "TerminationTimeUnit: " << this->TerminationTimeUnit << endl;
1387 os << indent << "StaticSeeds: " << this->StaticSeeds << endl;
1388 }
1389 //------------------------------------------------------------------------------
ComputeDomainExitLocation(double pos[4],double p2[4],double intersection[4],vtkGenericCell * cell)1390 bool vtkTemporalStreamTracer::ComputeDomainExitLocation(
1391 double pos[4], double p2[4], double intersection[4], vtkGenericCell* cell)
1392 {
1393 double t, pcoords[3];
1394 int subId;
1395 if (cell->IntersectWithLine(pos, p2, 1E-3, t, intersection, pcoords, subId) == 0)
1396 {
1397 vtkDebugMacro(<< "No cell/domain exit was found");
1398 return false;
1399 }
1400 else
1401 {
1402 // We found an intersection on the edge of the cell.
1403 // Shift it by a small amount to ensure that it crosses over the edge
1404 // into the adjoining cell.
1405 for (int i = 0; i < 3; i++)
1406 intersection[i] = pos[i] + (t + 0.01) * (p2[i] - pos[i]);
1407 // intersection stored, compute T for intersection
1408 intersection[3] = pos[3] + (t + 0.01) * (p2[3] - pos[3]);
1409 return true;
1410 }
1411 }
1412 //------------------------------------------------------------------------------
AddParticleToMPISendList(ParticleInformation & info)1413 void vtkTemporalStreamTracer::AddParticleToMPISendList(ParticleInformation& info)
1414 {
1415 double eps = (this->CurrentTimeSteps[1] - this->CurrentTimeSteps[0]) / 100;
1416 if (info.CurrentPosition.x[3] < (this->CurrentTimeSteps[0] - eps) ||
1417 info.CurrentPosition.x[3] > (this->CurrentTimeSteps[1] + eps))
1418 {
1419 vtkDebugMacro(<< "Unexpected time value in MPISendList - expected ("
1420 << this->CurrentTimeSteps[0] << "-" << this->CurrentTimeSteps[1] << ") got "
1421 << info.CurrentPosition.x[3]);
1422 }
1423 }
1424 //------------------------------------------------------------------------------
1425 /*
1426 // Now try to compute the trajectory exit location from the cell on the edge
1427 if (this->Interpolator->GetLastValidCellId(0)!=-1) {
1428 vtkDataSet *hitdata = this->Interpolator->GetLastDataSet(0);
1429 hitdata->GetCell(this->Interpolator->GetLastValidCellId(0), this->GenericCell);
1430 double intersection[4];
1431 if (this->ComputeDomainExitLocation(point1, point2, intersection, this->GenericCell))
1432 {
1433 vtkDebugMacro(<< "SendParticleToAnotherProcess : Sending Particle " << particleId << " Time "
1434 << intersection[3]); this->AddParticleToMPISendList(intersection, velocity); vtkIdType nextPoint =
1435 ParticleCoordinates->InsertNextPoint(point2);
1436 this->ParticleHistories[particleId].push_back(nextPoint);
1437 this->LiveParticleIds[particleId] = 0;
1438 return 1;
1439 }
1440 else {
1441 vtkDebugMacro(<< "Domain-Exit aborted : Domain Intersection failed" << particleId );
1442 this->LiveParticleIds[particleId] = 0;
1443 return 0;
1444 }
1445 }
1446 else {
1447 vtkDebugMacro(<< "Domain-Exit aborted : Couldn't copy cell from earlier test" << particleId );
1448 this->LiveParticleIds[particleId] = 0;
1449 return 0;
1450 }
1451 */
1452