Commit 4db27d8e authored by Andrew Bauer's avatar Andrew Bauer
Browse files

More particletracer cleanup that shouldn't change behavior.

Change-Id: If4b3a83fcb842bb0a692df10efbd7da71ef29141
parent 90cbb8f8
......@@ -14,36 +14,29 @@ PURPOSE. See the above copyright notice for more information.
=========================================================================*/
#include "vtkParticleTracerBase.h"
#include "vtkAbstractParticleWriter.h"
#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkCharArray.h"
#include "vtkCompositeDataIterator.h"
#include "vtkCompositeDataPipeline.h"
#include "vtkDataSetAttributes.h"
#include "vtkDoubleArray.h"
#include "vtkExecutive.h"
#include "vtkGenericCell.h"
#include "vtkIdList.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkIntArray.h"
#include "vtkFloatArray.h"
#include "vtkDoubleArray.h"
#include "vtkCharArray.h"
#include "vtkMath.h"
#include "vtkMultiBlockDataSet.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkPointSet.h"
#include "vtkPolyData.h"
#include "vtkPolyLine.h"
#include "vtkRungeKutta2.h"
#include "vtkRungeKutta4.h"
#include "vtkRungeKutta45.h"
#include "vtkSmartPointer.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkTemporalInterpolatedVelocityField.h"
#include "vtkOutputWindow.h"
#include "vtkAbstractParticleWriter.h"
#include "vtkToolkits.h" // For VTK_USE_MPI
#include <cassert>
#ifdef WIN32
......@@ -74,18 +67,6 @@ using namespace vtkParticleTracerBaseNamespace;
vtkCxxSetObjectMacro(vtkParticleTracerBase, ParticleWriter, vtkAbstractParticleWriter);
vtkCxxSetObjectMacro(vtkParticleTracerBase,Integrator,vtkInitialValueProblemSolver); //XXX
#define ParticleTracerSetMacro(name,type) \
void vtkParticleTracerBase::Set##name (type _arg) \
{ \
if (this->name == _arg) \
{\
return; \
}\
this->name = _arg; \
this->ResetCache(); \
this->Modified(); \
}
namespace
{
//return the interval i, such that a belongs to the interval (A[i],A[i+1]]
......@@ -111,8 +92,6 @@ namespace
//---------------------------------------------------------------------------
vtkParticleTracerBase::vtkParticleTracerBase()
{
this->SetNumberOfInputPorts(2);
// by default process active point vectors
this->SetInputArrayToProcess(0,0,0,vtkDataObject::FIELD_ASSOCIATION_POINTS,
vtkDataSetAttributes::VECTORS);
......@@ -140,14 +119,9 @@ vtkParticleTracerBase::vtkParticleTracerBase()
this->RotationScale = 1.0;
this->MaximumError = 1.0e-6;
this->TerminalSpeed = vtkParticleTracerBase::Epsilon;
//
this->IntegrationStep = 0.5;
// we are not actually using these for now
//
this->Interpolator = vtkSmartPointer<vtkTemporalInterpolatedVelocityField>::New();
//
this->SetNumberOfInputPorts(2);
#ifdef JB_H5PART_PARTICLE_OUTPUT
......@@ -165,27 +139,23 @@ vtkParticleTracerBase::vtkParticleTracerBase()
this->SetIntegratorType(RUNGE_KUTTA4);
this->DisableResetCache = 0;
}
//---------------------------------------------------------------------------
vtkParticleTracerBase::~vtkParticleTracerBase()
{
this->SetParticleWriter(NULL);
if (this->ParticleFileName)
{
delete []this->ParticleFileName;
this->ParticleFileName = NULL;
}
this->SetParticleFileName(NULL);
this->CachedData[0] = NULL;
this->CachedData[1] = NULL;;
this->CachedData[1] = NULL;
this->SetIntegrator(0);
this->SetInterpolatorPrototype(0);
}
//----------------------------------------------------------------------------
int vtkParticleTracerBase::FillInputPortInformation(
int port,
vtkInformation* info)
int port, vtkInformation* info)
{
// port 0 must be a temporal collection of any type
// the executive should put a temporal collection in when
......@@ -202,16 +172,19 @@ int vtkParticleTracerBase::FillInputPortInformation(
}
return 1;
}
//----------------------------------------------------------------------------
void vtkParticleTracerBase::AddSourceConnection(vtkAlgorithmOutput* input)
{
this->AddInputConnection(1, input);
}
//----------------------------------------------------------------------------
void vtkParticleTracerBase::RemoveAllSources()
{
this->SetInputConnection(1, 0);
}
//----------------------------------------------------------------------------
int vtkParticleTracerBase::ProcessRequest(
vtkInformation* request,
......@@ -235,6 +208,7 @@ int vtkParticleTracerBase::ProcessRequest(
}
return 1;
}
//----------------------------------------------------------------------------
int vtkParticleTracerBase::RequestInformation(
vtkInformation *vtkNotUsed(request),
......@@ -246,7 +220,8 @@ int vtkParticleTracerBase::RequestInformation(
if (inInfo->Has(vtkStreamingDemandDrivenPipeline::TIME_STEPS()) )
{
unsigned int numberOfInputTimeSteps = inInfo->Length(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
unsigned int numberOfInputTimeSteps =
inInfo->Length(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
vtkDebugMacro(<<"vtkParticleTracerBase "
"inputVector TIME_STEPS " << numberOfInputTimeSteps);
//
......@@ -254,8 +229,8 @@ int vtkParticleTracerBase::RequestInformation(
this->InputTimeValues.resize(numberOfInputTimeSteps);
inInfo->Get( vtkStreamingDemandDrivenPipeline::TIME_STEPS(),
&this->InputTimeValues[0] );
if (numberOfInputTimeSteps==1 && this->DisableResetCache==0) //warning would be skipped in coprocessing work flow
{
if (numberOfInputTimeSteps==1 && this->DisableResetCache==0)
{ //warning would be skipped in coprocessing work flow
vtkWarningMacro(<<"Not enough input time steps for particle integration");
}
......@@ -268,7 +243,6 @@ int vtkParticleTracerBase::RequestInformation(
{
this->SetStartTime(this->InputTimeValues.back());
}
}
else
{
......@@ -279,20 +253,22 @@ int vtkParticleTracerBase::RequestInformation(
outInfo->Set(
vtkStreamingDemandDrivenPipeline::MAXIMUM_NUMBER_OF_PIECES(), -1);
return 1;
}
//----------------------------------------------------------------------------
class WithinTolerance: public std::binary_function<double, double, bool>
namespace
{
public:
//----------------------------------------------------------------------------
class WithinTolerance: public std::binary_function<double, double, bool>
{
public:
result_type operator()(first_argument_type a, second_argument_type b) const
{
bool result = (fabs(a-b)<=(a*1E-6));
return (result_type)result;
}
};
};
}
//----------------------------------------------------------------------------
int vtkParticleTracerBase::RequestUpdateExtent(
vtkInformation *vtkNotUsed(request),
......@@ -388,6 +364,7 @@ int vtkParticleTracerBase::RequestUpdateExtent(
return 1;
}
//---------------------------------------------------------------------------
int vtkParticleTracerBase::InitializeInterpolator()
{
......@@ -504,6 +481,8 @@ int vtkParticleTracerBase::InitializeInterpolator()
//
return VTK_OK;
}
//---------------------------------------------------------------------------
int vtkParticleTracerBase::UpdateDataCache(vtkDataObject *data)
{
double dataTime = data->GetInformation()->Get(vtkDataObject::DATA_TIME_STEP());
......@@ -579,7 +558,6 @@ int vtkParticleTracerBase::UpdateDataCache(vtkDataObject *data)
return 1;
}
//---------------------------------------------------------------------------
bool vtkParticleTracerBase::InsideBounds(double point[])
{
......@@ -596,6 +574,7 @@ bool vtkParticleTracerBase::InsideBounds(double point[])
}
return false;
}
//---------------------------------------------------------------------------
void vtkParticleTracerBase::TestParticles(
ParticleVector &candidates, ParticleVector &passed, int &count)
......@@ -611,7 +590,10 @@ void vtkParticleTracerBase::TestParticles(
}
}
void vtkParticleTracerBase::TestParticles(vtkParticleTracerBaseNamespace::ParticleVector &candidates, std::vector<int> &passed)
//---------------------------------------------------------------------------
void vtkParticleTracerBase::TestParticles(
vtkParticleTracerBaseNamespace::ParticleVector &candidates,
std::vector<int> &passed)
{
int i = 0;
for (ParticleIterator it=candidates.begin(); it!=candidates.end(); ++it, ++i)
......@@ -640,8 +622,8 @@ void vtkParticleTracerBase::TestParticles(vtkParticleTracerBaseNamespace::Partic
}
//---------------------------------------------------------------------------
void vtkParticleTracerBase::AssignSeedsToProcessors( double time,
vtkDataSet *source, int sourceID, int ptId,
void vtkParticleTracerBase::AssignSeedsToProcessors(
double time, vtkDataSet *source, int sourceID, int ptId,
ParticleVector &localSeedPoints, int &localAssignedCount)
{
ParticleVector candidates;
......@@ -685,6 +667,7 @@ void vtkParticleTracerBase::AssignSeedsToProcessors( double time,
//
}
//---------------------------------------------------------------------------
void vtkParticleTracerBase::AssignUniqueIds(
vtkParticleTracerBaseNamespace::ParticleVector &localSeedPoints)
......@@ -711,6 +694,7 @@ void vtkParticleTracerBase::UpdateParticleList(ParticleVector &candidates)
vtkDebugMacro(<< "UpdateParticleList completed with " << this->NumberOfParticles() << " particles");
}
//---------------------------------------------------------------------------
int vtkParticleTracerBase::ProcessInput(vtkInformationVector** inputVector)
{
Assert(this->CurrentTimeStep>=StartTimeStep && this->CurrentTimeStep<=TerminationTimeStep);
......@@ -734,7 +718,7 @@ int vtkParticleTracerBase::ProcessInput(vtkInformationVector** inputVector)
return 1;
}
//---------------------------------------------------------------------------
vtkPolyData* vtkParticleTracerBase::Execute(vtkInformationVector** inputVector)
{
Assert(this->CurrentTimeStep>=this->StartTimeStep);
......@@ -1001,6 +985,7 @@ vtkPolyData* vtkParticleTracerBase::Execute(vtkInformationVector** inputVector)
return output;
}
//---------------------------------------------------------------------------
int vtkParticleTracerBase::RequestData(
vtkInformation *request,
vtkInformationVector **inputVector,
......@@ -1076,20 +1061,18 @@ int vtkParticleTracerBase::RequestData(
PRINT("RD: "<<this->CurrentTimeStep<<" "<<this->StartTimeStep<<" "<<this->TerminationTimeStep);
return 1;
}
//---------------------------------------------------------------------------
void vtkParticleTracerBase::IntegrateParticle(
ParticleListIterator &it,
double currenttime, double targettime,
ParticleListIterator &it, double currenttime, double targettime,
vtkInitialValueProblemSolver* integrator)
{
double epsilon = (targettime-currenttime)/100.0;
double velocity[3], point1[4], point2[4] = {0.0, 0.0, 0.0, 0.0};
double minStep=0, maxStep=0;
double stepWanted, stepTaken=0.0;
int substeps = 0;
ParticleInformation &info = (*it);
ParticleInformation previous = (*it);
bool particle_good = true;
......@@ -1132,7 +1115,6 @@ void vtkParticleTracerBase::IntegrateParticle(
maxStep = stepWanted;
}
// Calculate the next step using the integrator provided.
// If the next point is out of bounds, send it to another process
if (integrator->ComputeNextStep(
......@@ -1142,7 +1124,6 @@ void vtkParticleTracerBase::IntegrateParticle(
{
// if the particle is sent, remove it from the list
info.ErrorCode = 1;
if (!this->RetryWithPush(info, point1, delT, substeps))
{
if(previous.PointId <0)
......@@ -1250,7 +1231,7 @@ void vtkParticleTracerBase::IntegrateParticle(
Assert (point1[3]>=(this->GetCacheDataTime(0)-eps) && point1[3]<=(this->GetCacheDataTime(1)+eps));
#endif
}
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
void vtkParticleTracerBase::PrintSelf(ostream& os, vtkIndent indent)
{
......@@ -1267,6 +1248,7 @@ void vtkParticleTracerBase::PrintSelf(ostream& os, vtkIndent indent)
os << indent << "TerminationTime: " << this->TerminationTime << endl;
os << indent << "StaticSeeds: " << this->StaticSeeds << endl;
}
//---------------------------------------------------------------------------
bool vtkParticleTracerBase::ComputeDomainExitLocation(
double pos[4], double p2[4], double intersection[4],
......@@ -1293,6 +1275,7 @@ bool vtkParticleTracerBase::ComputeDomainExitLocation(
return 1;
}
}
//---------------------------------------------------------------------------
void vtkParticleTracerBase::SetIntegratorType(int type)
{
......@@ -1319,7 +1302,7 @@ void vtkParticleTracerBase::SetIntegratorType(int type)
}
}
//---------------------------------------------------------------------------
int vtkParticleTracerBase::GetIntegratorType()
{
if (!this->Integrator)
......@@ -1341,7 +1324,10 @@ int vtkParticleTracerBase::GetIntegratorType()
return UNKNOWN;
}
void vtkParticleTracerBase::CalculateVorticity( vtkGenericCell* cell, double pcoords[3], vtkDoubleArray* cellVectors, double vorticity[3] )
//---------------------------------------------------------------------------
void vtkParticleTracerBase::CalculateVorticity(
vtkGenericCell* cell, double pcoords[3],
vtkDoubleArray* cellVectors, double vorticity[3] )
{
double* cellVel;
double derivs[9];
......@@ -1353,11 +1339,13 @@ void vtkParticleTracerBase::CalculateVorticity( vtkGenericCell* cell, double pco
vorticity[2] = derivs[3] - derivs[1];
}
//---------------------------------------------------------------------------
double vtkParticleTracerBase::GetCacheDataTime(int i)
{
return this->CachedData[i]->GetInformation()->Get(vtkDataObject::DATA_TIME_STEP());
}
//---------------------------------------------------------------------------
double vtkParticleTracerBase::GetCacheDataTime()
{
double currentTime = CachedData[1]? this->GetCacheDataTime(1):
......@@ -1366,11 +1354,13 @@ double vtkParticleTracerBase::GetCacheDataTime()
}
//---------------------------------------------------------------------------
unsigned int vtkParticleTracerBase::NumberOfParticles()
{
return static_cast<unsigned int>(this->ParticleHistories.size());
}
//---------------------------------------------------------------------------
void vtkParticleTracerBase::ResetCache()
{
if(this->DisableResetCache==0)
......@@ -1389,6 +1379,7 @@ void vtkParticleTracerBase::ResetCache()
}
}
//---------------------------------------------------------------------------
void vtkParticleTracerBase::SetTerminationTime(double t)
{
if(t==this->TerminationTime)
......@@ -1411,7 +1402,7 @@ void vtkParticleTracerBase::SetTerminationTime(double t)
this->TerminationTime = t;
}
//---------------------------------------------------------------------------
void vtkParticleTracerBase::CreateProtoPD(vtkDataObject* input)
{
this->ProtoPD = NULL;
......@@ -1438,14 +1429,11 @@ void vtkParticleTracerBase::CreateProtoPD(vtkDataObject* input)
}
//---------------------------------------------------------------------------
bool vtkParticleTracerBase::RetryWithPush(ParticleInformation &info, double* point1,double delT, int substeps)
bool vtkParticleTracerBase::RetryWithPush(
ParticleInformation &info, double* point1,double delT, int substeps)
{
double velocity[3];
this->Interpolator->ClearCache();
if (info.UniqueParticleId==3)
{
vtkDebugMacro(<< "3 is about to be sent");
}
info.LocationState = this->Interpolator->TestPoint(point1);
......@@ -1491,6 +1479,7 @@ bool vtkParticleTracerBase::RetryWithPush(ParticleInformation &info, double* po
info.CurrentPosition.x[3] += delT;
info.LocationState = this->Interpolator->TestPoint(info.CurrentPosition.x);
if (info.LocationState!=ID_OUTSIDE_ALL)
{
// a push helped the particle get back into a dataset,
......@@ -1501,7 +1490,9 @@ bool vtkParticleTracerBase::RetryWithPush(ParticleInformation &info, double* po
return 0;
}
void vtkParticleTracerBase::AddParticle( vtkParticleTracerBaseNamespace::ParticleInformation &info, double* velocity)
//---------------------------------------------------------------------------
void vtkParticleTracerBase::AddParticle(
vtkParticleTracerBaseNamespace::ParticleInformation &info, double* velocity)
{
const double *coord = info.CurrentPosition.x;
vtkIdType tempId = this->OutputCoordinates->InsertNextPoint(coord);
......@@ -1581,6 +1572,7 @@ void vtkParticleTracerBase::AddParticle( vtkParticleTracerBaseNamespace::Particl
}
//---------------------------------------------------------------------------
bool vtkParticleTracerBase::IsPointDataValid(vtkDataObject* input)
{
if(vtkCompositeDataSet* cdInput = vtkCompositeDataSet::SafeDownCast(input))
......@@ -1592,6 +1584,7 @@ bool vtkParticleTracerBase::IsPointDataValid(vtkDataObject* input)
return true;
}
//---------------------------------------------------------------------------
bool vtkParticleTracerBase::IsPointDataValid(vtkCompositeDataSet* input,
std::vector<std::string>& arrayNames)
{
......@@ -1616,6 +1609,7 @@ bool vtkParticleTracerBase::IsPointDataValid(vtkCompositeDataSet* input,
return true;
}
//---------------------------------------------------------------------------
void vtkParticleTracerBase::GetPointDataArrayNames(
vtkDataSet* input, std::vector<std::string>& names)
{
......@@ -1631,66 +1625,69 @@ void vtkParticleTracerBase::GetPointDataArrayNames(
}
}
//---------------------------------------------------------------------------
vtkFloatArray* vtkParticleTracerBase::GetParticleAge(vtkPointData* pd)
{
return vtkFloatArray::SafeDownCast(pd->GetArray("ParticleAge"));
}
//---------------------------------------------------------------------------
vtkIntArray* vtkParticleTracerBase::GetParticleIds(vtkPointData* pd)
{
return vtkIntArray::SafeDownCast(pd->GetArray("ParticleId"));
}
//---------------------------------------------------------------------------
vtkCharArray* vtkParticleTracerBase::GetParticleSourceIds(vtkPointData* pd)
{
return vtkCharArray::SafeDownCast(pd->GetArray("ParticleSourceId"));
}
//---------------------------------------------------------------------------
vtkIntArray* vtkParticleTracerBase::GetInjectedPointIds(vtkPointData* pd)
{
return vtkIntArray::SafeDownCast(pd->GetArray("InjectedPointId"));
}
//---------------------------------------------------------------------------
vtkIntArray* vtkParticleTracerBase::GetInjectedStepIds(vtkPointData* pd)
{
return vtkIntArray::SafeDownCast(pd->GetArray("InjectionStepId"));
}
//---------------------------------------------------------------------------
vtkIntArray* vtkParticleTracerBase::GetErrorCodeArr(vtkPointData* pd)
{
return vtkIntArray::SafeDownCast(pd->GetArray("ErrorCode"));
}
//---------------------------------------------------------------------------
vtkFloatArray* vtkParticleTracerBase::GetParticleVorticity(vtkPointData* pd)
{
return vtkFloatArray::SafeDownCast(pd->GetArray("Vorticity"));
}
//---------------------------------------------------------------------------
vtkFloatArray* vtkParticleTracerBase::GetParticleRotation(vtkPointData* pd)
{
return vtkFloatArray::SafeDownCast(pd->GetArray("Rotation"));
}
//---------------------------------------------------------------------------
vtkFloatArray* vtkParticleTracerBase::GetParticleAngularVel(vtkPointData* pd)
{
return vtkFloatArray::SafeDownCast(pd->GetArray("AngularVelocity"));
}
//---------------------------------------------------------------------------
void vtkParticleTracerBase::PrintParticleHistories()
{
cout<<"Particle id, ages: "<<endl;
for(ParticleListIterator itr = this->ParticleHistories.begin(); itr!=this->ParticleHistories.end();itr++)
for(ParticleListIterator itr = this->ParticleHistories.begin();
itr!=this->ParticleHistories.end();itr++)
{
ParticleInformation& info(*itr);
cout<<info.InjectedPointId<<" "<<info.age<<" "<<endl;
}
cout<<endl;
}
ParticleTracerSetMacro(StartTime, double)
ParticleTracerSetMacro(ComputeVorticity, bool);
ParticleTracerSetMacro(RotationScale, double)
ParticleTracerSetMacro(ForceReinjectionEveryNSteps,int);
ParticleTracerSetMacro(TerminalSpeed, double);
......@@ -108,20 +108,19 @@ public:
// (necessary for generating proper stream-ribbons using the
// vtkRibbonFilter.
vtkGetMacro(ComputeVorticity, bool);
void SetComputeVorticity(bool);
vtkSetMacro(ComputeVorticity, bool);
// Description
// Specify the terminal speed value, below which integration is terminated.
vtkGetMacro(TerminalSpeed, double);
void SetTerminalSpeed(double);
vtkSetMacro(TerminalSpeed, double);
// Description
// This can be used to scale the rate with which the streamribbons
// twist. The default is 1.
void SetRotationScale(double);
vtkSetMacro(RotationScale, double);
vtkGetMacro(RotationScale, double);
// Description:
// To get around problems with the Paraview Animation controls
// we can just animate the time step and ignore the TIME_ requests
......@@ -137,8 +136,8 @@ public:
// Note that if the particle source is also animated, this flag will be
// redundant as the particles will be reinjected whenever the source changes
// anyway
void SetForceReinjectionEveryNSteps(int);
vtkGetMacro(ForceReinjectionEveryNSteps,int);
vtkSetMacro(ForceReinjectionEveryNSteps,int);
// Description:
// Setting TerminationTime to a positive value will cause particles
......@@ -159,10 +158,9 @@ public:
// to terminate when the time is reached. Use a vlue of zero to
// diable termination. The units of time should be consistent with the
// primary time variable.
void SetStartTime(double t);
vtkSetMacro(StartTime, double);
vtkGetMacro(StartTime,double);
// Description:
// if StaticSeeds is set and the mesh is static,
// then every time particles are injected we can re-use the same
......@@ -294,7 +292,6 @@ public:
int InitializeInterpolator();
int UpdateDataCache(vtkDataObject *td);
// Description : Test the list of particles to see if they are
// inside our data. Add good ones to passed list and set count to the
// number that passed
......@@ -313,14 +310,14 @@ public:
// If either are non static, then this step is skipped.
virtual void AssignSeedsToProcessors(double time,
vtkDataSet *source, int sourceID, int ptId,
vtkParticleTracerBaseNamespace::ParticleVector &LocalSeedPoints,
int &LocalAssignedCount);
vtkParticleTracerBaseNamespace::ParticleVector &localSeedPoints,
int &localAssignedCount);
// Description : once seeds have been assigned to a process, we
// give each one a uniqu ID. We need to use MPI to find out
// who is using which numbers.
virtual void AssignUniqueIds(
vtkParticleTracerBaseNamespace::ParticleVector &LocalSeedPoints);
vtkParticleTracerBaseNamespace::ParticleVector &localSeedPoints);
// Description : copy list of particles from a vector used for testing particles
// and sending between processors, into a list, which is used as the master
......
......@@ -127,7 +127,6 @@ protected:
this->BoundingBox[5]=1;
}
void GetSpacing(double dx[3])
{
for(int i=0; i<3; i++)
......@@ -207,11 +206,18 @@ protected:
{
int* uExtent = outInfo->Get(