Commit 6774210d authored by Yuanxin Liu's avatar Yuanxin Liu
Browse files

Fix bug in streakline filter introduced by f329b260

Previous change to streakline filter is not correct when vtkParticleTracerBase::HasCache
is true. This is now fixed. In addition, to simplify the code, particle
(re)injection is done after advection.

Change-Id: I2c2077c386a0c7d7a5edf75eecf4e8367cdc5aac
parent 209f86ab
......@@ -164,7 +164,6 @@ vtkParticleTracerBase::vtkParticleTracerBase()
writer->Delete();
#endif
this->ForceReinjectionAtTermination = false;
this->SetIntegratorType(RUNGE_KUTTA4);
}
//---------------------------------------------------------------------------
......@@ -529,10 +528,6 @@ int vtkParticleTracerBase::UpdateDataCache(vtkDataObject *data)
{
this->CachedData[1] = this->CachedData[0];
}
// if(this->CachedData[1])
// {
// cout<<GetCacheDataTime(0)<<" "<<GetCacheDataTime(1)<<endl;
// }
return 1;
}
......@@ -597,7 +592,7 @@ void vtkParticleTracerBase::TestParticles(vtkParticleTracerBaseNamespace::Partic
}
//---------------------------------------------------------------------------
void vtkParticleTracerBase::AssignSeedsToProcessors(
void vtkParticleTracerBase::AssignSeedsToProcessors( double time,
vtkDataSet *source, int sourceID, int ptId,
ParticleVector &LocalSeedPoints, int &LocalAssignedCount)
{
......@@ -612,7 +607,7 @@ void vtkParticleTracerBase::AssignSeedsToProcessors(
{
ParticleInformation &info = candidates[i];
memcpy(&(info.CurrentPosition.x[0]), source->GetPoint(i), sizeof(double)*3);
info.CurrentPosition.x[3] = this->GetCacheDataTime(0);
info.CurrentPosition.x[3] = time;
info.LocationState = 0;
info.CachedCellId[0] =-1;
info.CachedCellId[1] =-1;
......@@ -792,18 +787,9 @@ vtkPolyData* vtkParticleTracerBase::Execute(vtkInformationVector** inputVector)
//
// Make sure the Particle Positions are initialized with Seed particles
//
bool injectionFlag = !this->HasCache; //if we haven't started tracing, we need to inject
if (this->ForceReinjectionEveryNSteps>0)
{
injectionFlag = (this->CurrentTimeStep%this->ForceReinjectionEveryNSteps)==0;
}
//
// Lists for seed particles
//
if (injectionFlag)
if (this->StartTime==this->CurrentTime)
{
Assert(!this->HasCache); //shouldn't have cache if restarting
int seedPointId=0;
if (!(this->StaticSeeds && this->AllFixedGeometry && this->LocalSeeds.size()==0))
{
......@@ -813,7 +799,7 @@ vtkPolyData* vtkParticleTracerBase::Execute(vtkInformationVector** inputVector)
for (size_t i=0; i<seedSources.size(); i++)
{
this->AssignSeedsToProcessors(seedSources[i], static_cast<int>(i), 0, this->LocalSeeds, seedPointId);
this->AssignSeedsToProcessors(this->CurrentTime,seedSources[i], static_cast<int>(i), 0, this->LocalSeeds, seedPointId);
}
this->ParticleInjectionTime.Modified();
......@@ -822,7 +808,7 @@ vtkPolyData* vtkParticleTracerBase::Execute(vtkInformationVector** inputVector)
vtkDebugMacro(<< "Reinjection about to update candidates (" << this->LocalSeeds.size() << " particles)");
this->UpdateParticleList(this->LocalSeeds);
this->ReinjectionCounter += 1;
}
}
if(this->CurrentTimeStep==this->StartTimeStep) //just add all the particles
{
......@@ -837,7 +823,7 @@ vtkPolyData* vtkParticleTracerBase::Execute(vtkInformationVector** inputVector)
}
}
else
{
{
ParticleListIterator it_first = this->ParticleHistories.begin();
ParticleListIterator it_last = this->ParticleHistories.end();
ParticleListIterator it_next;
......@@ -896,21 +882,31 @@ vtkPolyData* vtkParticleTracerBase::Execute(vtkInformationVector** inputVector)
}//end of pass
}
if(this->ForceReinjectionAtTermination && this->CurrentTime==this->TerminationTime) //reinject again in the last step
bool injectionFlag(false);
if (this->CurrentTime!=this->StartTime && this->ForceReinjectionEveryNSteps>0)
{
injectionFlag = (this->CurrentTimeStep - this->StartTimeStep)%this->ForceReinjectionEveryNSteps==0;
}
if(injectionFlag) //reinject again in the last step
{
this->ReinjectionCounter += 1;
ParticleListIterator lastParticle = --this->ParticleHistories.end();
int seedPointId=0;
this->LocalSeeds.clear();
for (size_t i=0; i<seedSources.size(); i++)
{
this->AssignSeedsToProcessors(seedSources[i], static_cast<int>(i), 0, this->LocalSeeds, seedPointId);
this->AssignSeedsToProcessors(this->CurrentTime,seedSources[i], static_cast<int>(i), 0, this->LocalSeeds, seedPointId);
}
this->ParticleInjectionTime.Modified();
this->ParticleHistories.clear();
this->UpdateParticleList(this->LocalSeeds);
this->ReinjectionCounter += 1;
for(ParticleListIterator itr = this->ParticleHistories.begin(); itr!=this->ParticleHistories.end();itr++)
ParticleListIterator itr = lastParticle;
for(++itr; itr!=this->ParticleHistories.end(); ++itr)
{
ParticleInformation& info(*itr);
ParticleInformation& info(*lastParticle);
this->Interpolator->TestPoint(info.CurrentPosition.x);
double velocity[3];
this->Interpolator->GetLastGoodVelocity(velocity);
......@@ -939,6 +935,12 @@ vtkPolyData* vtkParticleTracerBase::Execute(vtkInformationVector** inputVector)
this->ExecuteTime.Modified();
this->HasCache = true;
PRINT("Output "<<output->GetNumberOfPoints()<<" particles");
//Check post condition
for (ParticleListIterator it=this->ParticleHistories.begin(); it!=this->ParticleHistories.end();it++)
{
Assert( (*it).CurrentPosition.x[3] == this->CurrentTime);
}
return output;
}
......@@ -1320,6 +1322,12 @@ void vtkParticleTracerBase::SetTerminationTime(double t)
this->ResetCache();
}
if( t <= this->StartTime)
{
vtkWarningMacro("Can't go backward");
t = this->StartTime;
}
this->Modified();
this->TerminationTime = t;
}
......
......@@ -260,6 +260,12 @@ public:
//
virtual int ProcessInput(vtkInformationVector** inputVector);
// This is the main part of the algorithm:
// * move all the particles one step
// * Reinject particles (by adding them to this->ParticleHistories)
// either at the beginning or at the end of each step (modulo this->ForceReinjectionEveryNSteps)
// * Output a polydata representing the moved particles
// Note that if the starting and the ending time coincide, the polydata is still valid.
virtual vtkPolyData* Execute(vtkInformationVector** inputVector);
// the RequestData will call these methods in turn
......@@ -290,7 +296,7 @@ public:
// they belong to. This saves us retesting at every injection time
// providing 1) The volumes are static, 2) the seed points are static
// If either are non static, then this step is skipped.
virtual void AssignSeedsToProcessors(
virtual void AssignSeedsToProcessors(double time,
vtkDataSet *source, int sourceID, int ptId,
vtkParticleTracerBaseNamespace::ParticleVector &LocalSeedPoints,
int &LocalAssignedCount);
......@@ -369,8 +375,6 @@ public:
virtual void ResetCache();
void AddParticle(vtkParticleTracerBaseNamespace::ParticleInformation &info, double* velocity);
private:
// Description:
// Hide this because we require a new interpolator type
......@@ -393,7 +397,6 @@ private:
bool ComputeVorticity;
double RotationScale;
double TerminalSpeed;
bool ForceReinjectionAtTermination; //whether to reinject seeds again at the termination step
// Important for Caching of Cells/Ids/Weights etc
int AllFixedGeometry;
......@@ -405,7 +408,7 @@ private:
double TerminationTime;
double CurrentTime;
unsigned int StartTimeStep;
unsigned int StartTimeStep; //InputTimeValues[StartTimeStep] <= StartTime <= InputTimeValues[StartTimeStep+1]
unsigned int CurrentTimeStep;
unsigned int TerminationTimeStep; //computed from start time
bool FirstIteration;
......@@ -429,7 +432,7 @@ private:
vtkSmartPointer<vtkTemporalInterpolatedVelocityField> Interpolator;
vtkAbstractInterpolatedVelocityField * InterpolatorPrototype;
// The input datasets which are stored by time step 0 and 1
// Data for time step CurrentTimeStep-1 and CurrentTimeStep
vtkSmartPointer<vtkMultiBlockDataSet> CachedData[2];
// Cache bounds info for each dataset we will use repeatedly
......
......@@ -64,7 +64,6 @@ StreaklineFilterInternal::StreaklineFilterInternal(vtkParticleTracerBase* filter
{
this->Filter->ForceReinjectionEveryNSteps = 1;
this->Filter->IgnorePipelineTime = 1;
this->Filter->ForceReinjectionAtTermination = true;
}
int StreaklineFilterInternal::OutputParticles(vtkPolyData* particles)
......
......@@ -126,13 +126,13 @@ bool vtkPParticleTracerBase::SendParticleToAnotherProcess(ParticleInformation &i
}
//---------------------------------------------------------------------------
void vtkPParticleTracerBase::AssignSeedsToProcessors(
void vtkPParticleTracerBase::AssignSeedsToProcessors(double t,
vtkDataSet *source, int sourceID, int ptId,
ParticleVector &LocalSeedPoints, int &LocalAssignedCount)
{
if(!this->Controller)
{
return Superclass::AssignSeedsToProcessors(source, sourceID, ptId,
return Superclass::AssignSeedsToProcessors(t, source, sourceID, ptId,
LocalSeedPoints, LocalAssignedCount);
}
ParticleVector candidates;
......@@ -146,7 +146,7 @@ void vtkPParticleTracerBase::AssignSeedsToProcessors(
{
ParticleInformation &info = candidates[i];
memcpy(&(info.CurrentPosition.x[0]), source->GetPoint(i), sizeof(double)*3);
info.CurrentPosition.x[3] = this->GetCacheDataTime(0);
info.CurrentPosition.x[3] = t;
info.LocationState = 0;
info.CachedCellId[0] =-1;
info.CachedCellId[1] =-1;
......
......@@ -103,7 +103,7 @@ public:
// they belong to. This saves us retesting at every injection time
// providing 1) The volumes are static, 2) the seed points are static
// If either are non static, then this step is skipped.
virtual void AssignSeedsToProcessors(
virtual void AssignSeedsToProcessors(double time,
vtkDataSet *source, int sourceID, int ptId,
vtkParticleTracerBaseNamespace::ParticleVector &LocalSeedPoints,
int &LocalAssignedCount);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment