Commit 8394083f authored by John Biddiscombe's avatar John Biddiscombe
Browse files

ENH:Major overhaul of the particle tracer code. Improve the caching by using...

ENH:Major overhaul of the particle tracer code. Improve the caching by using CellLocators to find cells in Untructired grids. Improve the caching of individual particle cell IDs from time step to timestep. Fix handling of dynamic/static meshes to optimizae the caching and construction of BSPTrees/Locators. Fix bug causing particle to disappear when passed between processors. Add Ternimation speed threshold back in, prepare the way for ParticleTerminationStrategy and ParticleInjectionStrategy. Optimize a little the handling of static seed reinjection in static meshes. Allow Multiple inputs of multiblock temporal data. Allow multiple inputs for Seed locations. Fix computation of rotation/vorticity etc.
parent d46a7e35
......@@ -16,6 +16,7 @@ ENDIF(VTK_HAS_EXODUS)
SET ( Kit_SRCS
vtkBranchExtentTranslator.cxx
vtkCachingInterpolatedVelocityField.cxx
vtkCollectGraph.cxx
vtkCollectPolyData.cxx
vtkCollectTable.cxx
......
/*=========================================================================
Program: Visualization Toolkit
Module: vtkCachingInterpolatedVelocityField.cxx
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
All rights reserved.
See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
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.
=========================================================================*/
#include "vtkCachingInterpolatedVelocityField.h"
#include "vtkDataArray.h"
#include "vtkDataSet.h"
#include "vtkPointData.h"
#include "vtkGenericCell.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkAbstractCellLocator.h"
#include "vtkSmartPointer.h"
#include "vtkFloatArray.h"
#include "vtkDoubleArray.h"
#include "vtkCellLocator.h"
#define Custom_TreeType vtkCellLocator
#include <vtkstd/vector>
#ifdef JB_BSP_TREE
#include "vtkModifiedBSPTree.h"
#undef Custom_TreeType
#define Custom_TreeType vtkModifiedBSPTree
#endif
//---------------------------------------------------------------------------
vtkCxxRevisionMacro(vtkCachingInterpolatedVelocityField, "1.1");
vtkStandardNewMacro(vtkCachingInterpolatedVelocityField);
//---------------------------------------------------------------------------
const double IVFDataSetInfo::TOLERANCE_SCALE = 1.0E-8;
//---------------------------------------------------------------------------
IVFDataSetInfo::IVFDataSetInfo() {}
//---------------------------------------------------------------------------
void IVFDataSetInfo::SetDataSet(vtkDataSet *data, char *velocity, bool staticdataset, vtkAbstractCellLocator *locator)
{
this->VelocityFloat = NULL;
this->VelocityDouble = NULL;
this->DataSet = data;
this->Cell = vtkSmartPointer<vtkGenericCell>::New();
this->StaticDataSet = staticdataset;
if (locator) {
this->BSPTree = locator;
}
else if (this->DataSet->IsA("vtkUnstructuredGrid")) {
if (!this->BSPTree) {
this->BSPTree = vtkSmartPointer<Custom_TreeType>::New();
}
this->BSPTree->SetLazyEvaluation(1);
this->BSPTree->SetDataSet(this->DataSet);
this->BSPTree->SetUseExistingSearchStructure(this->StaticDataSet);
}
this->Tolerance =
this->DataSet->GetLength() * IVFDataSetInfo::TOLERANCE_SCALE;
//
vtkDataArray *vectors = this->DataSet->GetPointData()->GetArray(velocity);
if (vtkFloatArray::SafeDownCast(vectors)) {
this->VelocityFloat = vtkFloatArray::SafeDownCast(vectors)->GetPointer(0);
}
else if (vtkDoubleArray::SafeDownCast(vectors)) {
this->VelocityDouble = vtkDoubleArray::SafeDownCast(vectors)->GetPointer(0);
}
else {
vtkGenericWarningMacro("We only support float/double velocity vectors at the current time");
}
}
//---------------------------------------------------------------------------
IVFDataSetInfo::IVFDataSetInfo(const IVFDataSetInfo &ivfci)
{
this->VelocityFloat = ivfci.VelocityFloat;
this->VelocityDouble = ivfci.VelocityDouble;
this->DataSet = ivfci.DataSet;
this->Cell = ivfci.Cell;
this->BSPTree = ivfci.BSPTree;
this->Tolerance = ivfci.Tolerance;
this->StaticDataSet = ivfci.StaticDataSet;
}
//---------------------------------------------------------------------------
IVFDataSetInfo &IVFDataSetInfo::operator=(const IVFDataSetInfo &ivfci)
{
this->VelocityFloat = ivfci.VelocityFloat;
this->VelocityDouble = ivfci.VelocityDouble;
this->DataSet = ivfci.DataSet;
this->Cell = ivfci.Cell;
this->BSPTree = ivfci.BSPTree;
this->Tolerance = ivfci.Tolerance;
this->StaticDataSet = ivfci.StaticDataSet;
return *this;
}
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
vtkCachingInterpolatedVelocityField::vtkCachingInterpolatedVelocityField()
{
this->NumFuncs = 3; // u, v, w
this->NumIndepVars = 4; // x, y, z, t
this->VectorsSelection = NULL;
this->TempCell = vtkGenericCell::New();
this->CellCacheHit = 0;
this->DataSetCacheHit = 0;
this->CacheMiss = 0;
this->LastCacheIndex = 0;
this->Cache = NULL;
this->LastCellId = -1;
}
//---------------------------------------------------------------------------
vtkCachingInterpolatedVelocityField::~vtkCachingInterpolatedVelocityField()
{
this->NumFuncs = 0;
this->NumIndepVars = 0;
this->TempCell->Delete();
this->SetVectorsSelection(0);
}
//---------------------------------------------------------------------------
void vtkCachingInterpolatedVelocityField::SetDataSet(int I, vtkDataSet* dataset, bool staticdataset, vtkAbstractCellLocator *locator)
{
int N = vtkstd::max(I+1, static_cast<int>(this->CacheList.size()));
this->CacheList.resize(N);
this->CacheList[I].SetDataSet(dataset, this->VectorsSelection, staticdataset, locator);
//
int maxsize = vtkstd::max(static_cast<int>(this->Weights.size()), dataset->GetMaxCellSize());
this->Weights.assign(maxsize, 0.0);
}
//---------------------------------------------------------------------------
void vtkCachingInterpolatedVelocityField::SetLastCellInfo(vtkIdType c, int datasetindex)
{
if ((this->LastCacheIndex != datasetindex) || (this->LastCellId != c))
{
this->LastCacheIndex = datasetindex;
this->Cache = &this->CacheList[this->LastCacheIndex];
this->LastCellId = c;
// if the dataset changes, then the cached cell is invalidated
// we might as well prefetch the new cached cell - we'll need it on the next test anyway
if (this->LastCellId!=-1)
{
this->Cache->DataSet->GetCell(this->LastCellId, this->Cache->Cell);
}
}
}
//---------------------------------------------------------------------------
void vtkCachingInterpolatedVelocityField::ClearLastCellInfo() {
this->Cache = NULL;
this->LastCellId = -1;
}
//---------------------------------------------------------------------------
vtkGenericCell *vtkCachingInterpolatedVelocityField::GetLastCell()
{
if (this->Cache)
{
return this->Cache->Cell;
}
return NULL;
}
//---------------------------------------------------------------------------
// Evaluate {u,v,w} at {x,y,z,t}
int vtkCachingInterpolatedVelocityField::FunctionValues(double* x, double* f)
{
// Test using whatever cached information we have
if (this->Cache) {
if (this->FunctionValues(this->Cache, x, f)) {
this->DataSetCacheHit++;
return 1;
}
}
// don't reset this->Cache as we don't want to test it again,
// see == comparison below
int oldCacheIndex = this->LastCacheIndex;
// now try each of the datasets in turn
for (this->LastCacheIndex = 0;
this->LastCacheIndex < static_cast<int>(this->CacheList.size());
this->LastCacheIndex++)
{
IVFDataSetInfo *data = &this->CacheList[this->LastCacheIndex];
if (data==this->Cache) continue;
//
this->LastCellId = -1;
if (this->FunctionValues(data, x, f))
{
this->Cache = data;
this->CacheMiss++;
return 1;
}
}
// failed, so clear data and set the cache index to something sensible
this->CacheMiss++;
this->ClearLastCellInfo();
this->LastCacheIndex = oldCacheIndex;
return 0;
}
//---------------------------------------------------------------------------
// same as FunctionValues, but only testing in/out of cells
int vtkCachingInterpolatedVelocityField::InsideTest(double* x)
{
// Test using whatever cached information we have
if (this->Cache) {
// check the last cell
int subId;
if (this->LastCellId!=-1 && this->Cache->Cell->EvaluatePosition(
x, 0, subId, this->Cache->PCoords,
this->Cache->Tolerance, &this->Weights[0])==1) return 1;
// check this dataset
if (this->InsideTest(this->Cache, x)) return 1;
}
// don't reset this->Cache as we don't want to test it again,
// see != comparison below
// now try each of the other datasets in turn
for (this->LastCacheIndex = 0;
this->LastCacheIndex < static_cast<int>(this->CacheList.size());
this->LastCacheIndex++)
{
IVFDataSetInfo *data = &this->CacheList[this->LastCacheIndex];
if (data==this->Cache) continue;
//
this->LastCellId = -1;
if (this->InsideTest(data, x))
{
this->Cache = data;
return 1;
}
}
// failed, so clear data
this->ClearLastCellInfo();
return 0;
}
//---------------------------------------------------------------------------
int vtkCachingInterpolatedVelocityField::InsideTest(IVFDataSetInfo *data, double* x)
{
int cellId = data->BSPTree->FindCell(x, data->Tolerance, data->Cell, data->PCoords, &this->Weights[0]);
if (cellId!=-1)
{
this->LastCellId = cellId;
return 1;
}
return 0;
}
//---------------------------------------------------------------------------
// Evaluate {u,v,w} at {x,y,z,t}
int vtkCachingInterpolatedVelocityField::FunctionValues(
IVFDataSetInfo *data, double *x, double *f)
{
int subId;
double dist2;
if (this->LastCellId>=0)
{
bool inbox = true;
if (data->BSPTree && !data->BSPTree->InsideCellBounds(x, this->LastCellId)) {
inbox = false;
}
if (inbox && data->Cell->EvaluatePosition(
x, 0, subId, data->PCoords, dist2, &this->Weights[0])==1)
{
this->FastCompute(data, f);
this->CellCacheHit++;
return 1;
}
}
// we need to search the whole dataset
if (data->BSPTree)
{
int cellId = data->BSPTree->FindCell(x, data->Tolerance,
data->Cell, data->PCoords, &this->Weights[0]);
this->LastCellId = cellId;
}
else
{
data->DataSet->GetCell(this->LastCellId, this->TempCell);
this->LastCellId =
data->DataSet->FindCell(x, this->TempCell, data->Cell, -1,
data->Tolerance, subId, data->PCoords, &this->Weights[0]);
if (this->LastCellId != -1)
{
data->DataSet->GetCell(this->LastCellId, data->Cell);
}
}
// if the cell is valid
if (this->LastCellId != -1)
{
this->FastCompute(data, f);
}
else
{
return 0;
}
return 1;
}
//---------------------------------------------------------------------------
void vtkCachingInterpolatedVelocityField::FastCompute(
IVFDataSetInfo *data, double f[3])
{
f[0] = f[1] = f[2] = 0.0;
int numPts = data->Cell->GetNumberOfPoints();
// interpolate the vectors
double *dvectors = data->VelocityDouble;
if (dvectors)
{
for (int j=0; j<numPts; j++)
{
int id = data->Cell->PointIds->GetId(j);
f[0] += dvectors[id*3 + 0] * this->Weights[j];
f[1] += dvectors[id*3 + 1] * this->Weights[j];
f[2] += dvectors[id*3 + 2] * this->Weights[j];
}
}
else
{
float *fvectors = data->VelocityFloat;
for (int j=0; j<numPts; j++)
{
int id = data->Cell->PointIds->GetId(j);
f[0] += fvectors[id*3 + 0] * this->Weights[j];
f[1] += fvectors[id*3 + 1] * this->Weights[j];
f[2] += fvectors[id*3 + 2] * this->Weights[j];
}
}
}
//---------------------------------------------------------------------------
bool vtkCachingInterpolatedVelocityField::InterpolatePoint(
vtkPointData *outPD, vtkIdType outIndex)
{
if (!this->Cache || !this->Cache->DataSet)
{
return 0;
}
outPD->InterpolatePoint(
this->Cache->DataSet->GetPointData(), outIndex,
this->Cache->Cell->PointIds, &this->Weights[0]);
return 1;
}
//---------------------------------------------------------------------------
int vtkCachingInterpolatedVelocityField::GetLastWeights(double* w)
{
int j, numPts;
// If last cell is valid, fill w with the interpolation weights
// and return true
if (this->Cache && this->LastCellId >= 0)
{
numPts = this->Cache->Cell->GetNumberOfPoints();
for (j=0; j < numPts; j++)
{
w[j] = this->Weights[j];
}
return 1;
}
// otherwise, return false
else
{
return 0;
}
}
//---------------------------------------------------------------------------
int vtkCachingInterpolatedVelocityField::GetLastLocalCoordinates(double pcoords[3])
{
int j;
// If last cell is valid, fill p with the local coordinates
// and return true
if (this->Cache && this->LastCellId >= 0)
{
for (j=0; j < 3; j++)
{
pcoords[j] = this->Cache->PCoords[j];
}
return 1;
}
// otherwise, return false
else
{
return 0;
}
}
//---------------------------------------------------------------------------
void vtkCachingInterpolatedVelocityField::PrintSelf(ostream& os, vtkIndent indent)
{
this->Superclass::PrintSelf(os, indent);
os << indent << "Weights: " << &this->Weights[0] << endl;
os << indent << "Cell Cache hit: " << this->CellCacheHit << endl;
os << indent << "DataSet Cache hit: " << this->DataSetCacheHit << endl;
os << indent << "Cache miss: " << this->CacheMiss << endl;
os << indent << "VectorsSelection: "
<< (this->VectorsSelection?this->VectorsSelection:"(none)") << endl;
os << indent << "Cache->DataSet : "
<< this->Cache->DataSet << endl;
os << indent << "LastCacheIndex : "
<< this->LastCacheIndex << endl;
}
//---------------------------------------------------------------------------
/*=========================================================================
Program: Visualization Toolkit
Module: vtkCachingInterpolatedVelocityField.h
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
All rights reserved.
See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
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.
=========================================================================*/
// .NAME vtkCachingInterpolatedVelocityField - Interface for obtaining
// interpolated velocity values
// .SECTION Description
// vtkCachingInterpolatedVelocityField acts as a continuous velocity field
// by performing cell interpolation on the underlying vtkDataSet.
// This is a concrete sub-class of vtkFunctionSet with
// NumberOfIndependentVariables = 4 (x,y,z,t) and
// NumberOfFunctions = 3 (u,v,w). Normally, every time an evaluation
// is performed, the cell which contains the point (x,y,z) has to
// be found by calling FindCell. This is a computationally expensive
// operation. In certain cases, the cell search can be avoided or shortened
// by providing a guess for the cell id. For example, in streamline
// integration, the next evaluation is usually in the same or a neighbour
// cell. For this reason, vtkCachingInterpolatedVelocityField stores the last
// cell id. If caching is turned on, it uses this id as the starting point.
// .SECTION Caveats
// vtkCachingInterpolatedVelocityField is not thread safe. A new instance should
// be created by each thread.
// .SECTION See Also
// vtkFunctionSet vtkStreamer
// .SECTION TODO
// Need to clean up style to match vtk/Kitware standards. Please help.
#ifndef __vtkTInterpolatedVelocityField_h
#define __vtkTInterpolatedVelocityField_h
#include "vtkFunctionSet.h"
#include "vtkSmartPointer.h" // this is allowed
//BTX
#include <vtkstd/vector> // we need them
//ETX
class vtkDataSet;
class vtkDataArray;
class vtkPointData;
class vtkGenericCell;
class vtkAbstractCellLocator;
//BTX
//---------------------------------------------------------------------------
class IVFDataSetInfo;
//---------------------------------------------------------------------------
class IVFCacheList : public vtkstd::vector< IVFDataSetInfo > {};
//---------------------------------------------------------------------------
//ETX
class VTK_PARALLEL_EXPORT vtkCachingInterpolatedVelocityField : public vtkFunctionSet
{
public:
vtkTypeRevisionMacro(vtkCachingInterpolatedVelocityField,vtkFunctionSet);
virtual void PrintSelf(ostream& os, vtkIndent indent);
// Description:
// Construct a vtkCachingInterpolatedVelocityField with no initial data set.
// LastCellId is set to -1.
static vtkCachingInterpolatedVelocityField *New();
// Description:
// Evaluate the velocity field, f={u,v,w}, at {x, y, z}.
// returns 1 if valid, 0 if test failed
virtual int FunctionValues(double* x, double* f);
virtual int InsideTest(double* x);
// Description:
// Add a dataset used by the interpolation function evaluation.
virtual void SetDataSet(int I, vtkDataSet* dataset, bool staticdataset, vtkAbstractCellLocator *locator);
// Description:
// If you want to work with an arbitrary vector array, then set its name
// here. By default this in NULL and the filter will use the active vector
// array.
vtkGetStringMacro(VectorsSelection);
void SelectVectors(const char *fieldName)
{this->SetVectorsSelection(fieldName);}
// Description:
// Return the cell id cached from last evaluation.
void SetLastCellInfo(vtkIdType c, int datasetindex);
// Description:
// Set the last cell id to -1 so that the next search does not
// start from the previous cell
void ClearLastCellInfo();
// Description:
// Returns the interpolation weights/pcoords cached from last evaluation
// if the cached cell is valid (returns 1). Otherwise, it does not
// change w and returns 0.
int GetLastWeights(double* w);
int GetLastLocalCoordinates(double pcoords[3]);
// Description:
// Caching statistics.
vtkGetMacro(CellCacheHit, int);
vtkGetMacro(DataSetCacheHit, int);
vtkGetMacro(CacheMiss, int);
protected:
vtkCachingInterpolatedVelocityField();
~vtkCachingInterpolatedVelocityField();
vtkGenericCell *TempCell;
int CellCacheHit;
int DataSetCacheHit;
int CacheMiss;
int LastCacheIndex;
int LastCellId;
IVFDataSetInfo *Cache;
IVFCacheList CacheList;
char *VectorsSelection;
//BTX
vtkstd::vector<double> Weights;
//ETX
vtkSetStringMacro(VectorsSelection);
// private versions which work on the passed dataset/cache
// these do the real computation
int FunctionValues(IVFDataSetInfo *cache, double *x, double *f);
int InsideTest(IVFDataSetInfo *cache, double* x);
//BTX
friend class vtkTemporalInterpolatedVelocityField;
// Description:
// If all weights have been computed (parametric coords etc all valid)
// then we can quickly interpolate a scalar/vector using the known weights
// and the generic cell which has been stored.
// This function is primarily reserved for use by
// vtkTemporalInterpolatedVelocityField
void FastCompute(IVFDataSetInfo *cache, double f[3]);
bool InterpolatePoint(vtkPointData *outPD, vtkIdType outIndex);
vtkGenericCell *GetLastCell();
//ETX
private:
vtkCachingInterpolatedVelocityField(const vtkCachingInterpolatedVelocityField&); // Not implemented.
void operator=(const vtkCachingInterpolatedVelocityField&); // Not implemented.
};
//---------------------------------------------------------------------------
///////////////////////////////////////////////////////////////////////////////
// IVFDataSetInfo
///////////////////////////////////////////////////////////////////////////////
#ifndef DOXYGEN_SHOULD_SKIP_THIS
//
//BTX
//
class IVFDataSetInfo
{
public:
vtkSmartPointer<vtkDataSet> DataSet;
vtkSmartPointer<vtkAbstractCellLocator> BSPTree;
vtkSmartPointer<vtkGenericCell> Cell;
double PCoords[3];
float *VelocityFloat;
double *VelocityDouble;
double Tolerance;
bool StaticDataSet;
IVFDataSetInfo();
IVFDataSetInfo(const IVFDataSetInfo &ivfci);
IVFDataSetInfo &operator=(const IVFDataSetInfo &ivfci);
void SetDataSet(vtkDataSet *data, char *velocity, bool staticdataset, vtkAbstractCellLocator *locator);
//
static const double TOLERANCE_SCALE;
};
//
//ETX