Commit 86b1bff1 authored by Yuanxin Liu's avatar Yuanxin Liu
Browse files

modify stream tracer to allow AMR input

Change-Id: Icf990a1d09e842c12053004158cf9d11604b2a76
parent aa00b571
/*=========================================================================
Program: Visualization Toolkit
Module: vtkCompositeInterpolatedVelocityField.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 "vtkCompositeInterpolatedVelocityField.h"
#include "vtkMath.h"
#include "vtkDataSet.h"
#include "vtkDataArray.h"
#include "vtkPointData.h"
#include "vtkGenericCell.h"
#include "vtkObjectFactory.h"
const double vtkCompositeInterpolatedVelocityField::TOLERANCE_SCALE = 1.0E-8;
vtkCompositeInterpolatedVelocityField::vtkCompositeInterpolatedVelocityField()
{
this->LastDataSetIndex = 0;
this->DataSets = new vtkCompositeInterpolatedVelocityFieldDataSetsType;
}
vtkCompositeInterpolatedVelocityField::~vtkCompositeInterpolatedVelocityField()
{
if ( this->DataSets )
{
delete this->DataSets;
this->DataSets = NULL;
}
}
void vtkCompositeInterpolatedVelocityField::PrintSelf( ostream & os, vtkIndent indent )
{
this->Superclass::PrintSelf( os, indent );
os << indent << "DataSets: " << this->DataSets << endl;
os << indent << "Last Dataset Index: " << this->LastDataSetIndex << endl;
}
/*=========================================================================
Program: Visualization Toolkit
Module: vtkCompositeInterpolatedVelocityField.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 vtkCompositeInterpolatedVelocityField - An abstract class for
// obtaining the interpolated velocity values at a point
//
// .SECTION Description
// vtkCompositeInterpolatedVelocityField acts as a continuous velocity field
// by performing cell interpolation on the underlying vtkDataSet. This is an
// abstract sub-class of vtkFunctionSet, NumberOfIndependentVariables = 4
// (x,y,z,t) and NumberOfFunctions = 3 (u,v,w). With a brute-force scheme,
// every time an evaluation is performed, the target cell containing point
// (x,y,z) needs to be found by calling FindCell(), via either vtkDataSet or
// vtkAbstractCelllocator's sub-classes (vtkCellLocator & vtkModifiedBSPTree).
// As it incurs a large cost, one (for vtkCellLocatorInterpolatedVelocityField
// via vtkAbstractCellLocator) or two (for vtkInterpolatedVelocityField via
// vtkDataSet that involves vtkPointLocator in addressing vtkPointSet) levels
// of cell caching may be exploited to increase the performance.
//
// For vtkInterpolatedVelocityField, level #0 begins with intra-cell caching.
// Specifically if the previous cell is valid and the next point is still in
// it ( i.e., vtkCell::EvaluatePosition() returns 1, coupled with newly created
// parametric coordinates & weights ), the function values can be interpolated
// and only vtkCell::EvaluatePosition() is invoked. If this fails, then level #1
// follows by inter-cell search for the target cell that contains the next point.
// By an inter-cell search, the previous cell provides an important clue or serves
// as an immediate neighbor to aid in locating the target cell via vtkPointSet::
// FindCell(). If this still fails, a global cell location / search is invoked via
// vtkPointSet::FindCell(). Here regardless of either inter-cell or global search,
// vtkPointLocator is in fact employed (for datasets of type vtkPointSet only, note
// vtkImageData and vtkRectilinearGrid are able to provide rapid and robust cell
// location due to the simple mesh topology) as a crucial tool underlying the cell
// locator. However, the use of vtkPointLocator makes vtkInterpolatedVelocityField
// non-robust in cell location for vtkPointSet.
//
// For vtkCellLocatorInterpolatedVelocityField, the only caching (level #0) works
// by intra-cell trial. In case of failure, a global search for the target cell is
// invoked via vtkAbstractCellLocator::FindCell() and the actual work is done by
// either vtkCellLocator or vtkModifiedBSPTree (for datasets of type vtkPointSet
// only, while vtkImageData and vtkRectilinearGrid themselves are able to provide
// fast robust cell location). Without the involvement of vtkPointLocator, robust
// cell location is achieved for vtkPointSet.
//
// .SECTION Caveats
// vtkCompositeInterpolatedVelocityField is not thread safe. A new instance
// should be created by each thread.
// .SECTION See Also
// vtkInterpolatedVelocityField vtkCellLocatorInterpolatedVelocityField
// vtkGenericInterpolatedVelocityField vtkCachingInterpolatedVelocityField
// vtkTemporalInterpolatedVelocityField vtkFunctionSet vtkStreamer vtkStreamTracer
#ifndef __vtkCompositeInterpolatedVelocityField_h
#define __vtkCompositeInterpolatedVelocityField_h
#include "vtkAbstractInterpolatedVelocityField.h"
//BTX
#include <vector> // STL Header; Required for vector
//ETX
class vtkDataSet;
//BTX
class vtkDataArray;
//ETX
class vtkPointData;
class vtkGenericCell;
class vtkCompositeInterpolatedVelocityFieldDataSetsType;
class VTK_FILTERING_EXPORT vtkCompositeInterpolatedVelocityField : public vtkAbstractInterpolatedVelocityField
{
public:
vtkTypeMacro( vtkCompositeInterpolatedVelocityField, vtkAbstractInterpolatedVelocityField);
void PrintSelf( ostream & os, vtkIndent indent );
// Description:
// Get the most recently visited dataset and it id. The dataset is used
// for a guess regarding where the next point will be, without searching
// through all datasets. When setting the last dataset, care is needed as
// no reference counting or checks are performed. This feature is intended
// for custom interpolators only that cache datasets independently.
vtkGetMacro( LastDataSetIndex, int );
vtkGetObjectMacro( LastDataSet, vtkDataSet );
// Description:
// Add a dataset for implicit velocity function evaluation. If more than
// one dataset is added, the evaluation point is searched in all until a
// match is found. THIS FUNCTION DOES NOT CHANGE THE REFERENCE COUNT OF
// dataset FOR THREAD SAFETY REASONS.
virtual void AddDataSet( vtkDataSet * dataset ) = 0;
protected:
vtkCompositeInterpolatedVelocityField();
~vtkCompositeInterpolatedVelocityField();
static const double TOLERANCE_SCALE;
int LastDataSetIndex;
vtkCompositeInterpolatedVelocityFieldDataSetsType * DataSets;
private:
vtkCompositeInterpolatedVelocityField
( const vtkCompositeInterpolatedVelocityField & ); // Not implemented.
void operator = ( const vtkCompositeInterpolatedVelocityField & ); // Not implemented.
};
//BTX
typedef std::vector< vtkDataSet * > DataSetsTypeBase;
class vtkCompositeInterpolatedVelocityFieldDataSetsType: public DataSetsTypeBase { };
//ETX
#endif
......@@ -274,9 +274,11 @@ int vtkDistributedStreamTracer::ProcessTask(double seed[3],
{
input0 = vtkDataSet::SafeDownCast(iterP->GetCurrentDataObject());
}
vtkDataArray *vectors = this->GetInputArrayToProcess(0,input0);
int vecType(0);
vtkDataArray *vectors = this->GetInputArrayToProcess(0,input0,vecType);
const char *vecName = vectors->GetName();
this->Integrate(input0,
this->Integrate(input0->GetPointData(),
tmpOutput,
seeds,
seedIds,
......@@ -284,6 +286,7 @@ int vtkDistributedStreamTracer::ProcessTask(double seed[3],
lastPoint,
func,
maxCellSize,
vecType,
vecName,
propagation,
numSteps);
......
set(Module_SRCS
vtkAbstractInterpolatedVelocityField.cxx
vtkAMRInterpolatedVelocityField.cxx
vtkCompositeInterpolatedVelocityField.cxx
vtkDashedStreamLine.cxx
vtkInterpolatedVelocityField.cxx
vtkStreamer.cxx
......
set(MyTests
TestBSPTree.cxx
TestAMRInterpolatedVelocityField
)
include(${vtkTestingRendering_SOURCE_DIR}/vtkTestingObjectFactory.cmake)
......
/*=========================================================================
Program: Visualization Toolkit
Module: TestvtkAMRInterpolatedVelocityField.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 <vtkAMRBox.h>
#include <vtkAMRInterpolatedVelocityField.h>
#include <vtkAMRGaussianPulseSource.h>
#include <vtkNew.h>
#include <vtkGradientFilter.h>
#include <vtkOverlappingAMR.h>
#include <vtkMath.h>
#include <vtkCompositeDataPipeline.h>
#define RETURNONFALSE(b)\
if(!b) \
{\
vtkAlgorithm::SetDefaultExecutivePrototype(NULL);\
return EXIT_FAILURE;\
}
int TestAMRInterpolatedVelocityField(int, char*[])
{
vtkNew<vtkCompositeDataPipeline> cexec;
vtkAlgorithm::SetDefaultExecutivePrototype(cexec.GetPointer());
char name[100] = "Gaussian-Pulse";
vtkNew<vtkAMRGaussianPulseSource> imageSource;
vtkNew<vtkGradientFilter> gradientFilter;
gradientFilter->SetInputConnection(imageSource->GetOutputPort());
gradientFilter->SetInputScalars( vtkDataObject::FIELD_ASSOCIATION_CELLS,name);
gradientFilter->SetResultArrayName("Gradient");
gradientFilter->Update();
vtkOverlappingAMR* amrGrad = vtkOverlappingAMR::SafeDownCast(gradientFilter->GetOutputDataObject(0));
amrGrad->GenerateParentChildInformation();
vtkNew<vtkAMRInterpolatedVelocityField> func;
func->SetAMRData(amrGrad);
func->SelectVectors(vtkDataObject::FIELD_ASSOCIATION_CELLS,"Gradient");
double Points[4][3] =
{{-2.1,-0.51,1},
{-1.9,-0.51,1},
{-0.9,-0.51,1},
{-0.1,-0.51,1}
};
double v[3];
bool res;
res = func->FunctionValues(Points[0],v);
RETURNONFALSE(!res);
res = func->FunctionValues(Points[1],v);
RETURNONFALSE(res);
RETURNONFALSE(func->GetLastAMRBox().GetLevel()==1)
res = func->FunctionValues(Points[2],v);
RETURNONFALSE(res);
RETURNONFALSE(func->GetLastAMRBox().GetLevel()==0)
res = func->FunctionValues(Points[3],v);
RETURNONFALSE(res);
RETURNONFALSE(func->GetLastAMRBox().GetLevel()==1)
vtkAlgorithm::SetDefaultExecutivePrototype(NULL);
return EXIT_SUCCESS;
}
vtk_module(vtkFiltersTracers
DEPENDS
vtkCommonExecutionModel
vtkAMRCore
vtkFiltersGeneral
vtkFiltersSources
TEST_DEPENDS
......
#include "vtkAMRInterpolatedVelocityField.h"
#include "vtkObjectFactory.h"
#include "vtkUniformGrid.h"
#include "vtkOverlappingAMR.h"
#include <assert.h>
//----------------------------------------------------------------------------
namespace
{
bool Inside(double q[3], const vtkAMRBox& box)
{
double gbounds[6];
box.GetBounds(gbounds);
if ((q[0] < gbounds[0]) || (q[0] > gbounds[1]) ||
(q[1] < gbounds[2]) || (q[1] > gbounds[3]) ||
(q[2] < gbounds[4]) || (q[2] > gbounds[5]))
{
return false;
}
else
{
return true;
}
}
bool FindInLevel(double q[3], vtkOverlappingAMR *amrds, int level, unsigned int& gridId)
{
vtkAMRBox box;
for(unsigned int i = 0; i < amrds->GetNumberOfDataSets(level); i++ )
{
amrds->GetMetaData(level,i,box);
bool inside = Inside(q,box);
if(inside)
{
gridId = i;
return true;
}
}
return false;
}
};
bool vtkAMRInterpolatedVelocityField::FindGrid(double q[3],vtkOverlappingAMR *amrds,
unsigned int& level, unsigned int& gridId)
{
if (!FindInLevel(q, amrds, 0,gridId))
{
level = -1;
return false;
}
unsigned int maxLevels = amrds->GetNumberOfLevels();
for(level=0; level<maxLevels;level++)
{
unsigned int *children = amrds->GetChildren(level, gridId);
if (children == NULL)
{
break;
}
int n = children[0];
int i;
for (i = 1; i <= n; i++)
{
vtkAMRBox box;
amrds->GetMetaData(level+1,children[i],box);
if(Inside(q,box))
{
gridId = children[i];
break;
}
}
if(i>n)
{
break;
}
}
return true;
}
vtkStandardNewMacro(vtkAMRInterpolatedVelocityField);
vtkAMRInterpolatedVelocityField::vtkAMRInterpolatedVelocityField():
LastAMRBox(3)
{
this->Weights = new double[8];
this->AmrDataSet = NULL;
this->LastAMRBox.Invalidate();
}
vtkAMRInterpolatedVelocityField::~vtkAMRInterpolatedVelocityField()
{
delete [] this->Weights;
this->Weights = NULL;
}
void vtkAMRInterpolatedVelocityField::PrintSelf( ostream & os, vtkIndent indent )
{
this->Superclass::PrintSelf( os, indent );
}
void vtkAMRInterpolatedVelocityField::SetAMRData(vtkOverlappingAMR* amrds)
{
this->AmrDataSet = amrds;
}
int vtkAMRInterpolatedVelocityField::FunctionValues( double * x, double * f )
{
if(this->LastDataSet &&
this->FunctionValues(this->LastDataSet,x,f))
{
return 1;
}
//Either we do not know which data set it is, or existing LastDataSet does not contain x
//In any case, set LastDataSet to NULL and try to find a new one
this->LastDataSet = NULL;
this->LastAMRBox.Invalidate();
this->LastCellId = -1;
unsigned int level, gridId;
if(!FindGrid(x,this->AmrDataSet,level,gridId))
{
return 0;
}
vtkDataSet* ds = this->AmrDataSet->GetDataSet(level,gridId,this->LastAMRBox);
assert(!this->LastAMRBox.IsInvalid());
if(!ds)
{
return 0;
}
if(!this->FunctionValues(ds,x,f))
{
return 0;
}
this->LastDataSet = ds;
return 1;
}
bool vtkAMRInterpolatedVelocityField::SetLastDataSet(int level, int id)
{
this->LastDataSet = this->AmrDataSet->GetDataSet(level,id,this->LastAMRBox);
return this->LastDataSet!=NULL;
}
void vtkAMRInterpolatedVelocityField::SetLastCellId(vtkIdType, int )
{
vtkWarningMacro("Calling SetLastCellId has no effect")
}
/*=========================================================================
Program: Visualization Toolkit
Module: vtkAMRInterpolatedVelocityField.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 vtkAMRInterpolatedVelocityField - A concrete class for obtaining
// the interpolated velocity values at a point in AMR data.
//
//
// .SECTION Description
// The main functionality supported here is the point location inside vtkOverlappingAMR data set.
#ifndef __vtkAMRInterpolatedVelocityField_h
#define __vtkAMRInterpolatedVelocityField_h
#include "vtkFiltersTracersModule.h" // For export macro
#include <vtkAbstractInterpolatedVelocityField.h>
#include <vtkAMRBox.h> //needed for caching the last AMRBox
class vtkOverlappingAMR;
class VTKFILTERSTRACERS_EXPORT vtkAMRInterpolatedVelocityField
: public vtkAbstractInterpolatedVelocityField
{
public:
vtkTypeMacro( vtkAMRInterpolatedVelocityField,
vtkAbstractInterpolatedVelocityField );
static vtkAMRInterpolatedVelocityField * New();
vtkGetMacro(AmrDataSet,vtkOverlappingAMR*);
void SetAMRData(vtkOverlappingAMR* amr);
const vtkAMRBox& GetLastAMRBox() const { return LastAMRBox;}
bool SetLastDataSet(int level, int id);
//Description: This function is no op. Do not call
virtual void SetLastCellId( vtkIdType c, int dataindex );
// Description:
// Set the cell id cached by the last evaluation.
virtual void SetLastCellId( vtkIdType c )
{ this->Superclass::SetLastCellId( c ); }
// Description:
// Evaluate the velocity field f at point p
// If it succeeds, then both this->LastDataSet and this->LastAMRBox
// will be set according to where p is found.
// If it fails, either p is out of bound, in which case
// both this->LastDataSet and this->LastAMRBox are invlaid
// or, in a multi-process setting, p is inbound but not on the processor.
// In the last case, this->LastAMRBox is still valid, and points to
// the exact processor and data set on which p can be found.
virtual int FunctionValues( double * x, double * f );
void PrintSelf( ostream & os, vtkIndent indent );
// Descriptino:
// Point location routine.
static bool FindGrid(double q[3],vtkOverlappingAMR *amrds, unsigned int& level, unsigned int& gridId);
protected:
vtkOverlappingAMR* AmrDataSet;
vtkAMRBox LastAMRBox;
vtkAMRInterpolatedVelocityField();
~vtkAMRInterpolatedVelocityField();
virtual int FunctionValues( vtkDataSet * ds, double * x, double * f )
{ return this->Superclass::FunctionValues( ds, x, f ); }
private:
vtkAMRInterpolatedVelocityField(const vtkAMRInterpolatedVelocityField&); //Not implemented
void operator = ( const vtkAMRInterpolatedVelocityField& ); // Not implemented.
};
#endif
......@@ -40,17 +40,16 @@ vtkAbstractInterpolatedVelocityField::vtkAbstractInterpolatedVelocityField()
this->LastCellId = -1;
this->LastDataSet= 0;
this->LastDataSetIndex = 0;
this->LastPCoords[0] = 0.0;
this->LastPCoords[1] = 0.0;
this->LastPCoords[2] = 0.0;
this->VectorsType = 0;
this->VectorsSelection = 0;
this->NormalizeVector = false;
this->Cell = vtkGenericCell::New();
this->GenCell = vtkGenericCell::New();
this->DataSets = new vtkAbstractInterpolatedVelocityFieldDataSetsType;
}
//---------------------------------------------------------------------------
......@@ -60,7 +59,7 @@ vtkAbstractInterpolatedVelocityField::~vtkAbstractInterpolatedVelocityField()
this->NumIndepVars = 0;
this->LastDataSet = 0;
this->SetVectorsSelection( 0 );
this->SetVectorsSelection(0);
if ( this->Weights )
{
......@@ -79,12 +78,6 @@ vtkAbstractInterpolatedVelocityField::~vtkAbstractInterpolatedVelocityField()
this->GenCell->Delete();
this->GenCell = NULL;
}
if ( this->DataSets )
{
delete this->DataSets;
this->DataSets = NULL;
}
}
//---------------------------------------------------------------------------
......@@ -101,19 +94,19 @@ int vtkAbstractInterpolatedVelocityField::FunctionValues
// See if a dataset has been specified and if there are input vectors
if ( !dataset ||
!( vectors =
dataset->GetPointData()->GetVectors( this->VectorsSelection )
!(vectors = dataset->GetAttributesAsFieldData(this->VectorsType)->GetArray(this->VectorsSelection))
)
)
{
vtkErrorMacro( << "Can't evaluate dataset!" );
vectors = NULL;
return 0;
}