vtkParticleTracerBase.cxx 62.7 KB
Newer Older
1
2
/*=========================================================================

3
4
  Program:   Visualization Toolkit
  Module:    vtkParticleTracerBase.cxx
5

6
7
8
  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9

10
11
12
  This software is distributed WITHOUT ANY WARRANTY; without even
  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  PURPOSE.  See the above copyright notice for more information.
13

14
  =========================================================================*/
15
16
17
// VTK_DEPRECATED_IN_9_2_0() warnings for this class.
#define VTK_DEPRECATION_LEVEL 0

18
19
#include "vtkParticleTracerBase.h"

20
#include "vtkAbstractParticleWriter.h"
21
22
#include "vtkCellArray.h"
#include "vtkCellData.h"
23
24
#include "vtkCellLocatorStrategy.h"
#include "vtkClosestPointStrategy.h"
25
#include "vtkCompositeDataIterator.h"
26
#include "vtkDataObjectTreeRange.h"
27
#include "vtkDoubleArray.h"
28
#include "vtkFloatArray.h"
29
30
31
32
33
34
35
36
37
38
39
40
#include "vtkGenericCell.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkIntArray.h"
#include "vtkMath.h"
#include "vtkMultiBlockDataSet.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkPolyData.h"
#include "vtkRungeKutta2.h"
#include "vtkRungeKutta4.h"
#include "vtkRungeKutta45.h"
41
#include "vtkSMPTools.h"
42
#include "vtkSignedCharArray.h"
43
#include "vtkSmartPointer.h"
44
#include "vtkStreamingDemandDrivenPipeline.h"
45
46
#include "vtkTemporalInterpolatedVelocityField.h"

47
#include <algorithm>
48
#include <atomic>
49
#include <functional>
50
#include <mutex>
51
52
#ifdef DEBUGPARTICLETRACE
#define Assert(x) assert(x)
53
#define PRINT(x) cout << __LINE__ << ": " << x << endl;
54
55
56
57
58
#else
#define PRINT(x)
#define Assert(x)
#endif

florian maurin's avatar
florian maurin committed
59
60
61
62
// The 3D cell with the maximum number of points is VTK_LAGRANGE_HEXAHEDRON.
// We support up to 6th order hexahedra.
#define VTK_MAXIMUM_NUMBER_OF_POINTS 216

63
64
65
const double vtkParticleTracerBase::Epsilon = 1.0E-12;

using namespace vtkParticleTracerBaseNamespace;
66
using IDStates = vtkTemporalInterpolatedVelocityField::IDStates;
67

68
//------------------------------------------------------------------------------
69
vtkCxxSetObjectMacro(vtkParticleTracerBase, ParticleWriter, vtkAbstractParticleWriter);
70
vtkCxxSetObjectMacro(vtkParticleTracerBase, Integrator, vtkInitialValueProblemSolver);
71
72
73

// this SetMacro is different than the regular vtkSetMacro
// because it resets the cache as well.
74
75
76
77
78
79
80
81
82
83
#define ParticleTracerSetMacro(name, type)                                                         \
  void vtkParticleTracerBase::Set##name(type _arg)                                                 \
  {                                                                                                \
    if (this->name == _arg)                                                                        \
    {                                                                                              \
      return;                                                                                      \
    }                                                                                              \
    this->name = _arg;                                                                             \
    this->ResetCache();                                                                            \
    this->Modified();                                                                              \
84
  }
85
ParticleTracerSetMacro(StartTime, double);
86
ParticleTracerSetMacro(ComputeVorticity, bool);
87
ParticleTracerSetMacro(RotationScale, double);
88
ParticleTracerSetMacro(ForceReinjectionEveryNSteps, int);
89
ParticleTracerSetMacro(TerminalSpeed, double);
90
91
92

namespace
{
93
94
95
96
// return the interval i, such that a belongs to the interval (A[i],A[i+1]]
inline int FindInterval(double a, const std::vector<double>& A)
{
  if (A.empty() || a < A[0])
97
  {
98
99
    return -1;
  }
100

101
102
103
  for (size_t i = 0; i < A.size() - 1; i++)
  {
    if (a <= A[i + 1])
104
    {
105
      return static_cast<int>(i);
106
    }
107
  }
108
109
110

  return -1;
}
111
112
};

113
//------------------------------------------------------------------------------
114
115
116
vtkParticleTracerBase::vtkParticleTracerBase()
{
  // by default process active point vectors
117
118
  this->SetInputArrayToProcess(
    0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS, vtkDataSetAttributes::VECTORS);
119

120
121
  this->CurrentTimeStep = 0;
  this->CurrentTimeValue = 0;
122
  this->ForceReinjectionEveryNSteps = 0;
123
124
  this->ReinjectionCounter = 0;
  this->AllFixedGeometry = 1;
125
  this->MeshOverTime = MeshOverTimeTypes::DIFFERENT;
126
  this->StaticSeeds = 0;
127
  this->ComputeVorticity = true;
128
129
130
131
132
133
  this->IgnorePipelineTime = 1;
  this->ParticleWriter = nullptr;
  this->ParticleFileName = nullptr;
  this->EnableParticleWriting = false;
  this->UniqueIdCounter = 0;
  this->Integrator = nullptr;
134
135

  this->StartTime = 0.0;
136
  this->TerminationTime = 0.0;
137
138
139
  this->FirstIteration = true;
  this->HasCache = false;

140
141
  this->RotationScale = 1.0;
  this->MaximumError = 1.0e-6;
142
143
  this->TerminalSpeed = vtkParticleTracerBase::Epsilon;
  this->IntegrationStep = 0.5;
144

145
146
147
148
  this->Interpolator = vtkSmartPointer<vtkTemporalInterpolatedVelocityField>::New();
  this->SetNumberOfInputPorts(2);

#ifdef JB_H5PART_PARTICLE_OUTPUT
149
#ifdef _WIN32
150
151
  vtkDebugMacro(<< "Setting vtkH5PartWriter");
  vtkH5PartWriter* writer = vtkH5PartWriter::New();
152
#else
153
154
  vtkDebugMacro(<< "Setting vtkXMLParticleWriter");
  vtkXMLParticleWriter* writer = vtkXMLParticleWriter::New();
155
156
157
158
159
160
#endif
  this->SetParticleWriter(writer);
  writer->Delete();
#endif

  this->SetIntegratorType(RUNGE_KUTTA4);
161
  this->DisableResetCache = 0;
162
  this->ForceSerialExecution = false;
163
}
164

165
//------------------------------------------------------------------------------
166
167
vtkParticleTracerBase::~vtkParticleTracerBase()
{
168
169
  this->SetParticleWriter(nullptr);
  this->SetParticleFileName(nullptr);
170

171
172
  this->CachedData[0] = nullptr;
  this->CachedData[1] = nullptr;
173

174
  this->SetIntegrator(nullptr);
175
}
176

177
//------------------------------------------------------------------------------
178
int vtkParticleTracerBase::FillInputPortInformation(int port, vtkInformation* info)
179
180
181
182
{
  // port 0 must be a temporal collection of any type
  // the executive should put a temporal collection in when
  // we request multiple time steps.
183
  if (port == 0)
184
  {
185
186
    info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataObject");
    info->Set(vtkAlgorithm::INPUT_IS_REPEATABLE(), 1);
187
  }
188
  else if (port == 1)
189
  {
190
191
192
    info->Remove(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE());
    info->Append(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataObjectTree");
    info->Append(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataSet");
193
    info->Set(vtkAlgorithm::INPUT_IS_REPEATABLE(), 1);
194
  }
195
196
  return 1;
}
197

198
//------------------------------------------------------------------------------
199
200
201
202
void vtkParticleTracerBase::AddSourceConnection(vtkAlgorithmOutput* input)
{
  this->AddInputConnection(1, input);
}
203

204
//------------------------------------------------------------------------------
205
206
void vtkParticleTracerBase::RemoveAllSources()
{
207
  this->SetInputConnection(1, nullptr);
208
}
209

210
//------------------------------------------------------------------------------
211
vtkTypeBool vtkParticleTracerBase::ProcessRequest(
212
  vtkInformation* request, vtkInformationVector** inputVector, vtkInformationVector* outputVector)
213
{
214
  if (request->Has(vtkDemandDrivenPipeline::REQUEST_INFORMATION()))
215
  {
216
    if (this->FirstIteration)
217
    {
218
219
      return this->RequestInformation(request, inputVector, outputVector);
    }
220
  }
221
  if (request->Has(vtkStreamingDemandDrivenPipeline::REQUEST_UPDATE_EXTENT()))
222
  {
223
    return this->RequestUpdateExtent(request, inputVector, outputVector);
224
  }
225
  if (request->Has(vtkDemandDrivenPipeline::REQUEST_DATA()))
226
  {
227
    return this->RequestData(request, inputVector, outputVector);
228
  }
229
230
  return 1;
}
231

232
//------------------------------------------------------------------------------
233
234
int vtkParticleTracerBase::RequestInformation(vtkInformation* vtkNotUsed(request),
  vtkInformationVector** inputVector, vtkInformationVector* vtkNotUsed(outputVector))
235
{
236
  vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
237

238
  if (inInfo->Has(vtkStreamingDemandDrivenPipeline::TIME_STEPS()))
239
  {
240
241
    unsigned int numberOfInputTimeSteps =
      inInfo->Length(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
242
243
244
    vtkDebugMacro(<< "vtkParticleTracerBase "
                     "inputVector TIME_STEPS "
                  << numberOfInputTimeSteps);
245
246
    // Get list of input time step values
    this->InputTimeValues.resize(numberOfInputTimeSteps);
247
    inInfo->Get(vtkStreamingDemandDrivenPipeline::TIME_STEPS(), this->InputTimeValues.data());
248
249
250
    if (numberOfInputTimeSteps == 1 && this->DisableResetCache == 0)
    { // warning would be skipped in coprocessing work flow
      vtkWarningMacro(<< "Not enough input time steps for particle integration");
251
    }
252

253
254
    // clamp the default start time to be within the data time range
    if (this->StartTime < this->InputTimeValues[0])
255
    {
256
      this->StartTime = this->InputTimeValues[0];
257
    }
258
    else if (this->StartTime > this->InputTimeValues.back())
259
    {
260
      this->StartTime = this->InputTimeValues.back();
261
    }
262
  }
263
  else
264
  {
265
    vtkErrorMacro(<< "Input information has no TIME_STEPS set");
266
    return 0;
267
  }
268
269
270

  return 1;
}
271

272
//------------------------------------------------------------------------------
273
274
int vtkParticleTracerBase::RequestUpdateExtent(vtkInformation* vtkNotUsed(request),
  vtkInformationVector** inputVector, vtkInformationVector* outputVector)
275
276
{
  int numInputs = inputVector[0]->GetNumberOfInformationObjects();
277
  vtkInformation* outInfo = outputVector->GetInformationObject(0);
278

279
280
  PRINT("RUE: " << this->HasCache << " " << this->FirstIteration << " " << this->StartTime << " "
                << this->TerminationTime << " " << this->CurrentTimeStep);
281
282
  // The output has requested a time value, what times must we ask from our input
  // do this only for the first time
283
  if (this->FirstIteration)
284
  {
285
    if (this->InputTimeValues.size() == 1)
286
    {
287
      this->StartTimeStep = this->InputTimeValues[0] == this->StartTime ? 0 : -1;
288
    }
289
    else
290
    {
291
      this->StartTimeStep = FindInterval(this->StartTime, this->InputTimeValues);
292
    }
293

294
    if (this->StartTimeStep < 0)
295
    {
296
297
      vtkErrorMacro("Start time not in time range");
      return 0;
298
    }
299

300
301
    if (!this->IgnorePipelineTime &&
      outInfo->Has(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP()))
302
    {
303
      double terminationTime = outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP());
304
      PRINT("Pipeline has set termination time: " << terminationTime);
305
      this->SetTerminationTimeNoModify(terminationTime);
306
    }
307

308
    if (this->TerminationTime > this->InputTimeValues.back())
309
    {
310
      this->TerminationTime = this->InputTimeValues.back();
311
    }
312

313
    if (this->InputTimeValues.size() == 1)
314
    {
315
      this->TerminationTimeStep = this->TerminationTime == this->InputTimeValues[0] ? 0 : -1;
316
    }
317
    else
318
    {
319
      this->TerminationTimeStep = FindInterval(this->TerminationTime, this->InputTimeValues) + 1;
320
    }
321

322
    if (this->TerminationTimeStep < 0)
323
    {
324
325
      vtkErrorMacro("Termination time not in time range");
      return 0;
326
    }
327

328
    for (int i = 0; i < this->GetNumberOfInputPorts(); i++)
329
    {
330
      vtkInformation* info = this->GetInputPortInformation(i);
331
      if (info->Get(vtkAlgorithm::INPUT_IS_OPTIONAL()) && this->GetNumberOfInputConnections(i) == 0)
332
      {
333
        continue;
334
      }
335
336
337
338
      vtkAlgorithm* inputAlgorithm = this->GetInputAlgorithm(i, 0);
      vtkStreamingDemandDrivenPipeline* sddp =
        vtkStreamingDemandDrivenPipeline::SafeDownCast(inputAlgorithm->GetExecutive());
      if (sddp)
339
      {
340
        sddp->UpdatePipelineMTime();
Bill Lorensen's avatar
Bill Lorensen committed
341
        vtkMTimeType pmt = sddp->GetPipelineMTime();
342
        if (pmt > this->ExecuteTime.GetMTime())
343
        {
344
          PRINT("Reset cache of because upstream is newer")
345
          this->ResetCache();
346
347
        }
      }
348
    }
349
    if (!this->HasCache)
350
    {
351
      this->CurrentTimeStep = this->StartTimeStep;
352
      this->CurrentTimeValue = -DBL_MAX;
353
    }
354
  }
355

356
  for (int i = 0; i < numInputs; i++)
357
  {
358
359
    vtkInformation* inInfo = inputVector[0]->GetInformationObject(i);
    if (this->CurrentTimeStep < static_cast<int>(this->InputTimeValues.size()))
360
    {
361
362
      inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP(),
        this->InputTimeValues[this->CurrentTimeStep]);
363
    }
364
    else
365
    {
Andrew Bauer's avatar
Andrew Bauer committed
366
      Assert(this->CurrentTimeValue == this->InputTimeValues.back());
367
    }
368
  }
369
370
371

  return 1;
}
372

373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
//------------------------------------------------------------------------------
void vtkParticleTracerBase::SetMeshOverTime(int meshOverTime)
{
  if (this->MeshOverTime !=
    (meshOverTime < DIFFERENT ? DIFFERENT
                              : (meshOverTime > SAME_TOPOLOGY ? SAME_TOPOLOGY : meshOverTime)))
  {
    this->MeshOverTime =
      (meshOverTime < DIFFERENT ? DIFFERENT
                                : (meshOverTime > SAME_TOPOLOGY ? SAME_TOPOLOGY : meshOverTime));
    this->Modified();
    // Needed since the value needs to be set at the same time.
    this->Interpolator->SetMeshOverTime(this->MeshOverTime);
  }
}

389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
//------------------------------------------------------------------------------
void vtkParticleTracerBase::SetInterpolatorType(int interpolatorType)
{
  if (interpolatorType == INTERPOLATOR_WITH_CELL_LOCATOR)
  {
    // create an interpolator equipped with a cell locator (by default)
    vtkNew<vtkCellLocatorStrategy> strategy;
    this->Interpolator->SetFindCellStrategy(strategy);
  }
  else
  {
    // create an interpolator equipped with a point locator
    auto strategy = vtkSmartPointer<vtkClosestPointStrategy>::New();
    this->Interpolator->SetFindCellStrategy(strategy);
  }
}

//------------------------------------------------------------------------------
void vtkParticleTracerBase::SetInterpolatorTypeToDataSetPointLocator()
{
  this->SetInterpolatorType(static_cast<int>(INTERPOLATOR_WITH_DATASET_POINT_LOCATOR));
}

//------------------------------------------------------------------------------
void vtkParticleTracerBase::SetInterpolatorTypeToCellLocator()
{
  this->SetInterpolatorType(static_cast<int>(INTERPOLATOR_WITH_CELL_LOCATOR));
}

418
//------------------------------------------------------------------------------
419
420
int vtkParticleTracerBase::InitializeInterpolator()
{
421
  if (!this->CachedData[0] || !this->CachedData[1])
422
  {
423
424
    vtkErrorMacro("Missing data set to process.");
    return VTK_ERROR;
425
  }
426
427
428
429
430
431
  // When Multiblock arrays are processed, some may be empty
  // if the first is empty, we won't find the correct vector name
  // so scan until we get one
  vtkSmartPointer<vtkCompositeDataIterator> iterP;
  iterP.TakeReference(this->CachedData[0]->NewIterator());
  iterP->GoToFirstItem();
Sean McBride's avatar
Sean McBride committed
432
  const char* vecname = nullptr;
433
  while (!iterP->IsDoneWithTraversal())
434
  {
435
    vtkDataArray* vectors = this->GetInputArrayToProcess(0, iterP->GetCurrentDataObject());
436
    if (vectors)
437
    {
438
439
440
      vecname = vectors->GetName();
      break;
    }
441
442
    iterP->GoToNextItem();
  }
443
  if (!vecname)
444
  {
445
    vtkErrorMacro(<< "Couldn't find vector array " << vecname);
446
    return VTK_ERROR;
447
  }
448

449
  // set strategy if needed
450
451
452
453
454
  if (this->Interpolator->GetFindCellStrategy() == nullptr)
  {
    // cell locator is the default;
    this->SetInterpolatorTypeToCellLocator();
  }
455
456
  this->Interpolator->SelectVectors(vecname);

457
  vtkDebugMacro(<< "Interpolator using array " << vecname);
458
459
  int numValidInputBlocks[2] = { 0, 0 };
  int numTotalInputBlocks[2] = { 0, 0 };
460
  this->DataReferenceT[0] = this->DataReferenceT[1] = nullptr;
461
  for (int T = 0; T < 2; T++)
462
  {
463
464
465
466
467
    this->CachedBounds[T].clear();
    // iterate over all blocks of input and cache the bounds information
    // and determine fixed/dynamic mesh status.
    vtkSmartPointer<vtkCompositeDataIterator> anotherIterP;
    anotherIterP.TakeReference(this->CachedData[T]->NewIterator());
468
469
    for (anotherIterP->GoToFirstItem(); !anotherIterP->IsDoneWithTraversal();
         anotherIterP->GoToNextItem())
470
    {
471
472
473
      numTotalInputBlocks[T]++;
      vtkDataSet* inp = vtkDataSet::SafeDownCast(anotherIterP->GetCurrentDataObject());
      if (inp)
474
      {
475
        if (inp->GetNumberOfCells() == 0)
476
        {
477
          vtkDebugMacro("Skipping an empty dataset");
478
        }
479
        else if (!inp->GetPointData()->GetVectors(vecname) && inp->GetNumberOfPoints() > 0)
480
        {
481
          vtkDebugMacro("One of the input datasets has no velocity vector.");
482
        }
483
        else
484
        {
485
486
487
488
489
          // store the bounding boxes of each local dataset for faster 'point-in-dataset' testing
          bounds bbox;
          inp->GetBounds(&bbox.b[0]);
          this->CachedBounds[T].push_back(bbox);
          // add the dataset to the interpolator
490
          this->Interpolator->AddDataSetAtTime(T, this->GetCacheDataTime(T), inp);
491
          if (!this->DataReferenceT[T])
492
          {
493
            this->DataReferenceT[T] = inp;
494
          }
495
496
497
          numValidInputBlocks[T]++;
        }
      }
498
    }
499
  }
500
  if (numValidInputBlocks[0] == 0 || numValidInputBlocks[1] == 0)
501
  {
502
    vtkErrorMacro("Not enough inputs have been found. Can not execute."
503
      << numValidInputBlocks[0] << " " << numValidInputBlocks[1]);
504
    return VTK_ERROR;
505
  }
506
507
  if (numValidInputBlocks[0] != numValidInputBlocks[1] &&
    this->MeshOverTime != MeshOverTimeTypes::DIFFERENT)
508
  {
509
    vtkErrorMacro(
510
511
      "MeshOverTime is set to STATIC/LINEAR_INTERPOLATION/SAME_TOPOLOGY but the number of "
      "datasets is different between time steps "
512
      << numValidInputBlocks[0] << " " << numValidInputBlocks[1]);
513
  }
514
515
  vtkDebugMacro("Number of Valid input blocks is " << numValidInputBlocks[0] << " from "
                                                   << numTotalInputBlocks[0]);
516
517
518
  vtkDebugMacro("AllFixedGeometry " << this->AllFixedGeometry);

  // force optimizations if StaticMesh is set.
519
520
  this->AllFixedGeometry = this->MeshOverTime == MeshOverTimeTypes::STATIC;
  if (this->MeshOverTime == MeshOverTimeTypes::STATIC)
521
  {
522
    vtkDebugMacro("Static Mesh over time optimizations Forced ON");
523
  }
524

525
526
  this->Interpolator->Initialize(this->CachedData[0], this->CachedData[1]);

527
528
  return VTK_OK;
}
529

530
//------------------------------------------------------------------------------
531
532
533
534
std::vector<vtkDataSet*> vtkParticleTracerBase::GetSeedSources(
  vtkInformationVector* inputVector, int vtkNotUsed(timeStep))
{
  std::vector<vtkDataSet*> seedSources;
535
  for (int idx = 0, max = inputVector->GetNumberOfInformationObjects(); idx < max; ++idx)
536
  {
537
    if (vtkInformation* inInfo = inputVector->GetInformationObject(idx))
538
    {
539
540
      auto datasets = vtkCompositeDataSet::GetDataSets(vtkDataObject::GetData(inInfo));
      seedSources.insert(seedSources.end(), datasets.begin(), datasets.end());
541
    }
542
  }
543
544
545
  return seedSources;
}

546
//------------------------------------------------------------------------------
547
int vtkParticleTracerBase::UpdateDataCache(vtkDataObject* data)
548
549
550
{
  double dataTime = data->GetInformation()->Get(vtkDataObject::DATA_TIME_STEP());

551
552
  Assert(dataTime >= this->GetCacheDataTime());
  if (dataTime == this->GetCacheDataTime())
553
  {
554
    return 1;
555
  }
556
557

  int i = 0;
558
  if (this->CurrentTimeStep == this->StartTimeStep)
559
  {
560
    i = 0;
561
  }
562
  else if (this->CurrentTimeStep == this->StartTimeStep + 1)
563
  {
564
    i = 1;
565
  }
566
  else
567
  {
568
569
    i = 1;
    this->CachedData[0] = this->CachedData[1];
570
    this->CachedData[1] = nullptr;
571
  }
572
573
574
575

  this->CachedData[i] = vtkSmartPointer<vtkMultiBlockDataSet>::New();

  // if simple dataset, add to our list, otherwise, add blocks
576
577
  vtkDataSet* dsInput = vtkDataSet::SafeDownCast(data);
  vtkMultiBlockDataSet* mbInput = vtkMultiBlockDataSet::SafeDownCast(data);
578
579

  if (dsInput)
580
  {
581
582
583
584
    vtkSmartPointer<vtkDataSet> copy;
    copy.TakeReference(dsInput->NewInstance());
    copy->ShallowCopy(dsInput);
    this->CachedData[i]->SetBlock(this->CachedData[i]->GetNumberOfBlocks(), copy);
585
  }
586
  else if (mbInput)
587
  {
588
589
590
    vtkSmartPointer<vtkCompositeDataIterator> iter;
    iter.TakeReference(mbInput->NewIterator());
    for (iter->InitTraversal(); !iter->IsDoneWithTraversal(); iter->GoToNextItem())
591
    {
592
      vtkDataSet* ds = vtkDataSet::SafeDownCast(iter->GetCurrentDataObject());
593
      if (ds)
594
      {
595
596
597
598
599
600
        vtkSmartPointer<vtkDataSet> copy;
        copy.TakeReference(ds->NewInstance());
        copy->ShallowCopy(ds);
        this->CachedData[i]->SetBlock(this->CachedData[i]->GetNumberOfBlocks(), copy);
      }
    }
601
  }
602
  else
603
  {
604
605
    vtkDebugMacro(
      "This filter cannot handle inputs of type: " << (data ? data->GetClassName() : "(none)"));
606
    return 0;
607
  }
608
609

  this->CachedData[i]->GetInformation()->Set(vtkDataObject::DATA_TIME_STEP(), dataTime);
610
  if (this->CurrentTimeStep == this->StartTimeStep)
611
  {
612
    this->CachedData[1] = this->CachedData[0];
613
  }
614
615
616
  return 1;
}

617
//------------------------------------------------------------------------------
618
619
620
bool vtkParticleTracerBase::InsideBounds(double point[])
{
  double delta[3] = { 0.0, 0.0, 0.0 };
621
  for (int t = 0; t < 2; ++t)
622
  {
623
    for (size_t i = 0; i < (this->CachedBounds[t].size()); ++i)
624
    {
625
      if (vtkMath::PointIsWithinBounds(point, &((this->CachedBounds[t])[i].b[0]), delta))
626
      {
627
628
629
        return true;
      }
    }
630
  }
631
632
  return false;
}
633

634
//------------------------------------------------------------------------------
635
void vtkParticleTracerBase::TestParticles(
636
  ParticleVector& candidates, ParticleVector& passed, int& count)
637
638
639
640
641
{
  std::vector<int> passedIndices;
  this->TestParticles(candidates, passedIndices);
  count = static_cast<int>(passedIndices.size());

642
  for (size_t i = 0; i < passedIndices.size(); i++)
643
  {
644
    passed.push_back(candidates[passedIndices[i]]);
645
  }
646
647
}

648
//------------------------------------------------------------------------------
649
void vtkParticleTracerBase::TestParticles(
650
  vtkParticleTracerBaseNamespace::ParticleVector& candidates, std::vector<int>& passed)
651
652
{
  int i = 0;
653
  for (ParticleIterator it = candidates.begin(); it != candidates.end(); ++it, ++i)
654
  {
655
656
    ParticleInformation& info = (*it);
    double* pos = &info.CurrentPosition.x[0];
657
658
    // if outside bounds, reject instantly
    if (this->InsideBounds(pos))
659
    {
660
      // since this is first test, avoid bad cache tests
661
662
      this->Interpolator->ClearCache();
      info.LocationState = this->Interpolator->TestPoint(pos);
663
      if (info.LocationState == IDStates::OUTSIDE_ALL /*|| location==IDStates::OUTSIDE_T0*/)
664
      {
665
666
        // can't really use this particle.
        vtkDebugMacro(<< "TestParticles rejected particle");
667
      }
668
      else
669
      {
670
671
672
673
674
        // get the cached ids and datasets from the TestPoint call
        this->Interpolator->GetCachedCellIds(info.CachedCellId, info.CachedDataSetId);
        passed.push_back(i);
      }
    }
675
  }
676
677
}

678
//------------------------------------------------------------------------------
679
680
void vtkParticleTracerBase::AssignSeedsToProcessors(double time, vtkDataSet* source, int sourceID,
  int ptId, ParticleVector& localSeedPoints, int& localAssignedCount)
681
682
683
{
  ParticleVector candidates;
  // take points from the source object and create a particle list
684
  vtkIdType numSeeds = source->GetNumberOfPoints();
685
  candidates.resize(numSeeds);
686
  for (vtkIdType i = 0; i < numSeeds; i++)
687
  {
688
689
    ParticleInformation& info = candidates[i];
    memcpy(&(info.CurrentPosition.x[0]), source->GetPoint(i), sizeof(double) * 3);
690
    info.CurrentPosition.x[3] = time;
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
    info.LocationState = 0;
    info.CachedCellId[0] = -1;
    info.CachedCellId[1] = -1;
    info.CachedDataSetId[0] = 0;
    info.CachedDataSetId[1] = 0;
    info.SourceID = sourceID;
    info.InjectedPointId = i + ptId;
    info.InjectedStepId = this->ReinjectionCounter;
    info.TimeStepAge = 0;
    info.UniqueParticleId = -1;
    info.rotation = 0.0;
    info.angularVel = 0.0;
    info.time = 0.0;
    info.age = 0.0;
    info.speed = 0.0;
    info.ErrorCode = 0;
    info.SimulationTime = this->GetCurrentTimeValue();
    info.PointId = -1;
    info.TailPointId = -1;
710
  }
711
  // Gather all Seeds to all processors for classification
712
  this->TestParticles(candidates, localSeedPoints, localAssignedCount);
713
714
715

  // Assign unique identifiers taking into account uneven distribution
  // across processes and seeds which were rejected
716
  this->AssignUniqueIds(localSeedPoints);
717
}
718

719
//------------------------------------------------------------------------------
720
void vtkParticleTracerBase::AssignUniqueIds(
721
  vtkParticleTracerBaseNamespace::ParticleVector& localSeedPoints)
722
{
723
  vtkIdType particleCountOffset = 0;
724
  vtkIdType numParticles = static_cast<vtkIdType>(localSeedPoints.size());
725
  for (vtkIdType i = 0; i < numParticles; i++)
726
  {
727
    localSeedPoints[i].UniqueParticleId = this->UniqueIdCounter + particleCountOffset + i;
728
  }
729
730
731
  this->UniqueIdCounter += numParticles;
}

732
//------------------------------------------------------------------------------
733
void vtkParticleTracerBase::UpdateParticleList(ParticleVector& candidates)
734
735
{
  int numSeedsNew = static_cast<int>(candidates.size());
736
  for (int i = 0; i < numSeedsNew; i++)
737
  {
738
739
    // allocate a new particle on the list and get a reference to it
    this->ParticleHistories.push_back(candidates[i]);
740
  }
741
742
  vtkDebugMacro(<< "UpdateParticleList completed with " << this->NumberOfParticles()
                << " particles");
743
744
}

745
//------------------------------------------------------------------------------
746
747
int vtkParticleTracerBase::ProcessInput(vtkInformationVector** inputVector)
{
748
  Assert(this->CurrentTimeStep >= StartTimeStep && this->CurrentTimeStep <= TerminationTimeStep);
749
  int numInputs = inputVector[0]->GetNumberOfInformationObjects();
750
  if (numInputs != 1)
751
  {
752
    if (numInputs == 0)
753
    {
754
755
756
      vtkErrorMacro(<< "No input found.");
      return 0;
    }
757
758
    vtkWarningMacro(<< "Multiple inputs founds. Use only the first one.");
  }
759

760
  vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
761
  if (inInfo)
762
  {
763
    vtkDataObject* input = inInfo->Get(vtkDataObject::DATA_OBJECT());
764
    this->UpdateDataCache(input);
765
  }
766
767
768
  return 1;
}

769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
//------------------------------------------------------------------------------
namespace vtkParticleTracerBaseNamespace
{
struct ParticleTracerFunctor
{
  vtkParticleTracerBase* PT;
  double FromTime;
  bool Sequential;

  std::vector<std::list<ParticleInformation>::iterator> ParticleHistories;
  std::atomic<vtkIdType> ParticleCount;
  std::mutex EraseMutex;

  vtkSMPThreadLocal<vtkSmartPointer<vtkInitialValueProblemSolver>> TLIntegrator;
  vtkSMPThreadLocal<vtkSmartPointer<vtkTemporalInterpolatedVelocityField>> TLInterpolator;
  vtkSMPThreadLocal<vtkSmartPointer<vtkDoubleArray>> TLCellVectors;

  ParticleTracerFunctor(vtkParticleTracerBase* pt, double fromTime, bool sequential)
    : PT(pt)
    , FromTime(fromTime)
    , Sequential(sequential)
  {
    this->ParticleCount = 0;
    size_t particleSize = pt->ParticleHistories.size();
    // Copy the particle histories into a vector for O(1) access
    this->ParticleHistories.reserve(particleSize);
    for (auto it = pt->ParticleHistories.begin(), end = pt->ParticleHistories.end(); it != end;
         ++it)
    {
      this->ParticleHistories.push_back(it);
    }

    this->PT->ResizeArrays(static_cast<vtkIdType>(particleSize));
  }

  void Initialize()
  {
    // Some data members of the local output require per-thread initialization.
    auto& interpolator = this->TLInterpolator.Local();
    interpolator.TakeReference(this->PT->Interpolator->NewInstance());
    interpolator->CopyParameters(this->PT->Interpolator);
    auto& integrator = this->TLIntegrator.Local();
    integrator.TakeReference(this->PT->GetIntegrator()->NewInstance());
    integrator->SetFunctionSet(interpolator);
    auto& cellVectors = this->TLCellVectors.Local();
    cellVectors.TakeReference(vtkDoubleArray::New());

    if (this->PT->ComputeVorticity)
    {
      cellVectors->SetNumberOfComponents(3);
      cellVectors->Allocate(3 * VTK_CELL_SIZE);
    }
  }

  void operator()(vtkIdType begin, vtkIdType end)
  {
    auto& integrator = this->TLIntegrator.Local();
    auto& interpolator = this->TLInterpolator.Local();
    auto& cellVectors = this->TLCellVectors.Local();
    const double& currentTime = this->FromTime;
    const double& targetTime = this->PT->CurrentTimeValue;

    for (vtkIdType i = begin; i < end; ++i)
    {
      auto it = this->ParticleHistories[i];
      this->PT->IntegrateParticle(it, currentTime, targetTime, integrator, interpolator,
        cellVectors, this->ParticleCount, this->EraseMutex, this->Sequential);
      if (this->PT->GetAbortExecute())
      {
        vtkErrorWithObjectMacro(this->PT, "Execute aborted");
        break;
      }
    }
  }

  void Reduce()
  {
    // squeeze possibly extra space
    this->PT->ResizeArrays(this->ParticleCount);
  }