Commit d604dce9 authored by hrchilds's avatar hrchilds

(1) Fix problem with VTK cell lookups for non-rectilinear grids

(2) Incorporate fixes from Christoph for advection not getting off
the ground with small stencils and big velocities.


git-svn-id: http://visit.ilight.com/svn/visit/trunk/src@6761 18c085ea-50e0-402c-830e-de6fd14e8384
parent cf4d8df2
......@@ -60,7 +60,6 @@ Consider the leaveDomains SLs and the balancing at the same time.
#include <vtkCellData.h>
#include <vtkDataSet.h>
#include <vtkFloatArray.h>
#include <vtkInterpolatedVelocityField.h>
#include <vtkLineSource.h>
#include <vtkPlaneSource.h>
#include <vtkPoints.h>
......@@ -72,6 +71,7 @@ Consider the leaveDomains SLs and the balancing at the same time.
#include <vtkGlyph3D.h>
#include <vtkVisItCellLocator.h>
#include <vtkVisItInterpolatedVelocityField.h>
#include <avtCallback.h>
#include <avtDatabase.h>
......@@ -1574,6 +1574,9 @@ avtStreamlineFilter::DomainToRank(DomainType &domain)
// Dave Pugmire, Tue Mar 31 17:01:17 EDT 2009
// Fix memory leak.
//
// Hank Childs, Thu Apr 2 17:58:09 CDT 2009
// Do our own interpolation. The previous one we used was too buggy for ugrids.
//
// ****************************************************************************
avtIVPSolver::Result
......@@ -1589,7 +1592,7 @@ avtStreamlineFilter::IntegrateDomain(avtStreamlineWrapper *slSeg,
<<slSeg->domain<<") HGZ = "<<haveGhostZones <<endl;
// prepare streamline integration ingredients
vtkInterpolatedVelocityField* velocity1=vtkInterpolatedVelocityField::New();
vtkVisItInterpolatedVelocityField* velocity1= vtkVisItInterpolatedVelocityField::New();
if (doPathlines)
{
// Our expression will be the active variable, so reset it.
......@@ -1607,20 +1610,18 @@ avtStreamlineFilter::IntegrateDomain(avtStreamlineWrapper *slSeg,
cellToPt1->SetInput(ds);
cellToPt1->Update();
velocity1->AddDataSet(cellToPt1->GetOutput());
velocity1->SetDataSet(cellToPt1->GetOutput());
}
else
velocity1->AddDataSet(ds);
velocity1->SetDataSet(ds);
velocity1->CachingOn();
vtkInterpolatedVelocityField* velocity2=NULL;
vtkVisItInterpolatedVelocityField* velocity2=NULL;
vtkDataSet *ds2 = NULL;
vtkCellDataToPointData *cellToPt2 = NULL;
double t1, t2;
if (doPathlines)
{
velocity2 = vtkInterpolatedVelocityField::New();
velocity2 = vtkVisItInterpolatedVelocityField::New();
ds2 = (vtkDataSet *) ds->NewInstance();
ds2->ShallowCopy(ds);
......@@ -1642,13 +1643,11 @@ avtStreamlineFilter::IntegrateDomain(avtStreamlineWrapper *slSeg,
cellToPt2->SetInput(ds2);
cellToPt2->Update();
velocity2->AddDataSet(cellToPt2->GetOutput());
velocity2->SetDataSet(cellToPt2->GetOutput());
}
else
velocity2->AddDataSet(ds2);
velocity2->SetDataSet(ds2);
velocity2->CachingOn();
std::string db = GetInput()->GetInfo().GetAttributes().GetFullDBName();
ref_ptr<avtDatabase> dbp = avtCallback::GetDatabase(db, 0, NULL);
if (*dbp == NULL)
......
......@@ -340,6 +340,10 @@ avtIVPDopri5::SetTolerances(const double& relt, const double& abst)
// Dave Pugmire, Tue Aug 19, 17:38:03 EDT 2008
// Consider a negative stepsize.
//
// Hank Childs, Thu Apr 2 16:40:08 PDT 2009
// Incorporate fix from Christoph to do a better job guessing the
// initial step size.
//
// ****************************************************************************
double
......@@ -347,62 +351,91 @@ avtIVPDopri5::GuessInitialStep(const avtIVPField* field,
const double& h_max,
const double& t_max)
{
// make a local copy since we may need to modify it
double local_h_max = h_max;
double direction = sign(1.0, t_max - t);
double sk, sqr;
double dnf = 0.0;
double dny = 0.0;
// loop until an estimate succeeds
while( true )
{
try
{
double sk, sqr;
double dnf = 0.0;
double dny = 0.0;
double h;
double h;
for(size_t i=0 ; i < y.dim() ; i++)
{
sk = abstol + reltol * std::abs(y[i]);
sqr = k1[i] / sk;
dnf += sqr * sqr;
sqr = y[i] / sk;
dny += sqr * sqr;
}
for(size_t i=0 ; i < y.dim() ; i++)
{
sk = abstol + reltol * std::abs(y[i]);
sqr = k1[i] / sk;
dnf += sqr * sqr;
sqr = y[i] / sk;
dny += sqr * sqr;
}
if( (dnf <= 1.0e-10) || (dny <= 1.0e-10) )
h = 1.0e-6;
else
h = sqrt( dny/dnf ) * 0.01;
if( (dnf <= 1.0e-10) || (dny <= 1.0e-10) )
h = 1.0e-6;
else
h = sqrt( dny/dnf ) * 0.01;
h = std::min( h, h_max );
h = sign( h, direction );
h = std::min( h, local_h_max );
h = sign( h, direction );
// perform an explicit Euler step
avtVec k3 = y + h * k1;
avtVec k2 = (*field)( t+h, k3 );
// perform an explicit Euler step
avtVec k3 = y + h * k1;
avtVec k2 = (*field)( t+h, k3 );
n_eval++;
n_eval++;
// estimate the second derivative of the solution
double der2 = 0.0;
// estimate the second derivative of the solution
double der2 = 0.0;
for( size_t i=0; i<y.dim(); i++)
{
sk = abstol + reltol * std::abs( y[i] );
sqr = ( k2[i] - k1[i] ) / sk;
der2 += sqr*sqr;
}
for( size_t i=0; i<y.dim(); i++)
{
sk = abstol + reltol * std::abs( y[i] );
sqr = ( k2[i] - k1[i] ) / sk;
der2 += sqr*sqr;
}
der2 = sqrt( der2 ) / h;
der2 = sqrt( der2 ) / h;
// step size is computed such that
// h**(1/5) * max( norm(k1), norm(der2) ) = 0.01
double der12 = std::max( std::abs(der2), sqrt(dnf) );
double h1;
// step size is computed such that
// h**(1/5) * max( norm(k1), norm(der2) ) = 0.01
double der12 = std::max( std::abs(der2), sqrt(dnf) );
double h1;
if( der12 <= 1.0e-15 )
h1 = std::max( 1.0e-6, std::abs(h)*1.0e-3 );
else
h1 = pow( 0.01/der12, 0.2 );
if( der12 <= 1.0e-15 )
h1 = std::max( 1.0e-6, std::abs(h)*1.0e-3 );
else
h1 = pow( 0.01/der12, 0.2 );
h = std::min( fabs(100.0*h), std::min( h1, h_max ) );
h = sign( h, direction );
h = std::min( fabs(100.0*h), std::min( h1, local_h_max ) );
h = sign( h, direction );
return h;
return h;
}
catch( avtIVPField::Undefined& )
{
// Somehow we couldn't evaluate one of the points we need for the
// starting estimate. The above code adheres to the h_max that is
// passed in - let's reduce that and try again.
// (we're really using local_h_max since h_max is const double&)
// (Oh, and if local_h_max is zero, let's set it to unit length (= direction))
if( !local_h_max )
local_h_max = direction;
else
local_h_max /= 2;
// if local_h_max is smaller then epsilon, we stop trying
// return that back to Step() which will then fail with a
// stepsize underflow
if( local_h_max < std::numeric_limits<double>::epsilon() )
return local_h_max;
}
}
}
......
......@@ -44,7 +44,7 @@
#include <iostream>
#include <vtkCell.h>
#include <vtkDataSet.h>
#include <vtkInterpolatedVelocityField.h>
#include <vtkVisItInterpolatedVelocityField.h>
#include <vtkDoubleArray.h>
#include <vtkPointData.h>
#include <DebugStream.h>
......@@ -55,9 +55,14 @@
// Programmer: Christoph Garth
// Creation: February 25, 2008
//
// Modifications:
//
// Hank Childs, Thu Apr 2 16:40:08 PDT 2009
// Use vtkVisItInterpolatedVelocityField, not vtkInterpolatedVelocityField.
//
// ****************************************************************************
avtIVPVTKField::avtIVPVTKField( vtkInterpolatedVelocityField* velocity )
avtIVPVTKField::avtIVPVTKField( vtkVisItInterpolatedVelocityField* velocity )
: iv(velocity)
{
iv->Register( NULL );
......@@ -88,6 +93,11 @@ avtIVPVTKField::~avtIVPVTKField()
// Programmer: Christoph Garth
// Creation: February 25, 2008
//
// Modifications:
//
// Hank Childs, Thu Apr 2 16:40:08 PDT 2009
// Use vtkVisItInterpolatedVelocityField, not vtkInterpolatedVelocityField.
//
// ****************************************************************************
avtVec
......@@ -95,8 +105,8 @@ avtIVPVTKField::operator()(const double& t, const avtVecRef& x) const
{
avtVec y( x.dim() ), param( pad(x,t));
int result = iv->FunctionValues( param.values(), y.values() );
//debug5<<result<<"= iv->FunctionValues("<<param<<") y= "<<y<<endl;
int result = iv->Evaluate( param.values(), y.values() );
//debug5<<result<<"= iv->Evaluate("<<param<<") y= "<<y<<endl;
if( !result )
throw Undefined();
......@@ -126,6 +136,9 @@ avtIVPVTKField::operator()(const double& t, const avtVecRef& x) const
// Increase the size of the "w" (weights) variable to prevent stack
// overwrites.
//
// Hank Childs, Thu Apr 2 16:40:08 PDT 2009
// Use vtkVisItInterpolatedVelocityField, not vtkInterpolatedVelocityField.
//
// ****************************************************************************
double
......@@ -134,13 +147,13 @@ avtIVPVTKField::ComputeVorticity( const double& t, const avtVecRef& x ) const
avtVec y( x.dim() );
avtVec param = pad(x,t);
int result = iv->FunctionValues( param.values(), y.values() );
int result = iv->Evaluate( param.values(), y.values() );
if( !result )
throw Undefined();
vtkDataSet *ds = iv->GetLastDataSet();
vtkIdType cellID = iv->GetLastCellId();
vtkDataSet *ds = iv->GetDataSet();
vtkIdType cellID = iv->GetLastCell();
vtkCell *cell = ds->GetCell( cellID );
vtkDoubleArray *cellVectors;
......@@ -155,9 +168,8 @@ avtIVPVTKField::ComputeVorticity( const double& t, const avtVecRef& x ) const
inVectors->GetTuples( cell->PointIds, cellVectors );
double *cellVel = cellVectors->GetPointer(0);
double pcoords[3], w[100];
iv->GetLastWeights( w );
iv->GetLastLocalCoordinates( pcoords );
double *w = iv->GetLastWeights();
double *pcoords = iv->GetLastPCoords();
cell->Derivatives( 0, pcoords, cellVel, 3, derivs);
//cout<<"pcoords= "<<pcoords[0]<<" "<<pcoords[1]<<" "<<pcoords[2]<<endl;
......@@ -191,13 +203,18 @@ avtIVPVTKField::ComputeVorticity( const double& t, const avtVecRef& x ) const
// Programmer: Christoph Garth
// Creation: February 25, 2008
//
// Modifications:
//
// Hank Childs, Thu Apr 2 16:40:08 PDT 2009
// Use vtkVisItInterpolatedVelocityField, not vtkInterpolatedVelocityField.
//
// ****************************************************************************
bool
avtIVPVTKField::IsInside( const double& t, const avtVecRef& x ) const
{
avtVec y( x.dim() );
return iv->FunctionValues( pad( x, t ).values(), y.values() );
return iv->Evaluate( pad( x, t ).values(), y.values() );
}
......@@ -210,12 +227,18 @@ avtIVPVTKField::IsInside( const double& t, const avtVecRef& x ) const
// Programmer: Christoph Garth
// Creation: February 25, 2008
//
// Modifications:
//
// Hank Childs, Thu Apr 2 16:40:08 PDT 2009
// Use vtkVisItInterpolatedVelocityField, not vtkInterpolatedVelocityField.
// (The old method was just returning 0 anyways...)
//
// ****************************************************************************
unsigned int
avtIVPVTKField::GetDimension() const
{
return iv->GetNumberOfFunctions();
return 3;
}
// ****************************************************************************
......
......@@ -45,7 +45,7 @@
#include <avtIVPField.h>
#include <vtkInterpolatedVelocityField.h>
#include <vtkVisItInterpolatedVelocityField.h>
#include <ivp_exports.h>
......@@ -66,12 +66,15 @@
// Dave Pugmire, Tue Mar 10 12:41:11 EDT 2009
// Added GetValidTimeRange function.
//
// Hank Childs, Thu Apr 2 16:40:08 PDT 2009
// Use vtkVisItInterpolatedVelocityField, not vtktInterpolatedVelocityField.
//
// ****************************************************************************
class IVP_API avtIVPVTKField: public avtIVPField
{
public:
avtIVPVTKField( vtkInterpolatedVelocityField* velocity );
avtIVPVTKField( vtkVisItInterpolatedVelocityField* velocity );
~avtIVPVTKField();
// avtIVPField interface
......@@ -84,7 +87,7 @@ class IVP_API avtIVPVTKField: public avtIVPField
virtual bool GetValidTimeRange(double range[]) const {return false;}
protected:
vtkInterpolatedVelocityField *iv;
vtkVisItInterpolatedVelocityField *iv;
bool normalized;
};
......
......@@ -37,14 +37,14 @@
*****************************************************************************/
// ************************************************************************* //
// avtIVPVTKTimeVaryingField.C //
// avtIVPVTKTimeVaryingField.C //
// ************************************************************************* //
#include <avtIVPVTKTimeVaryingField.h>
#include <iostream>
#include <vtkCell.h>
#include <vtkDataSet.h>
#include <vtkInterpolatedVelocityField.h>
#include <vtkVisItInterpolatedVelocityField.h>
#include <vtkDoubleArray.h>
#include <vtkPointData.h>
#include <DebugStream.h>
......@@ -57,8 +57,8 @@
//
// ****************************************************************************
avtIVPVTKTimeVaryingField::avtIVPVTKTimeVaryingField( vtkInterpolatedVelocityField *velocity1,
vtkInterpolatedVelocityField *velocity2,
avtIVPVTKTimeVaryingField::avtIVPVTKTimeVaryingField( vtkVisItInterpolatedVelocityField *velocity1,
vtkVisItInterpolatedVelocityField *velocity2,
double t1,
double t2)
{
......@@ -96,9 +96,15 @@ avtIVPVTKTimeVaryingField::~avtIVPVTKTimeVaryingField()
// Programmer: Dave Pugmire (on behalf of Hank Childs)
// Creation: Tue Feb 24 09:24:49 EST 2009
//
// Modifications:
//
// Dave Pugmire, Tue Mar 10 12:41:11 EDT 2009
// Check time bounds and throw execption.
//
// Hank Childs, Thu Apr 2 16:36:50 PDT 2009
// Update to use vtkVisItInterpolatedVelocityField, not
// vtkInterpolatedVelocityField.
//
// ****************************************************************************
avtVec
......@@ -115,8 +121,8 @@ avtIVPVTKTimeVaryingField::operator()(const double& t, const avtVecRef& x) const
// Evaluate the field at both timesteps.
avtVec y1(x.dim()), param(pad(x,t)), y2(x.dim());
if ( ! iv1->FunctionValues(param.values(), y1.values()) ||
! iv2->FunctionValues(param.values(), y2.values()) )
if ( ! iv1->Evaluate(param.values(), y1.values()) ||
! iv2->Evaluate(param.values(), y2.values()) )
{
debug5<<" **OUT of BOUNDS**\n";
throw Undefined();
......@@ -151,6 +157,10 @@ avtIVPVTKTimeVaryingField::operator()(const double& t, const avtVecRef& x) const
//
// Modifications:
//
// Hank Childs, Thu Apr 2 16:36:50 PDT 2009
// Update to use vtkVisItInterpolatedVelocityField, not
// vtkInterpolatedVelocityField.
//
// ****************************************************************************
double
......@@ -160,13 +170,13 @@ EXCEPTION0(ImproperUseException); // didn't do this.
avtVec y( x.dim() );
avtVec param = pad(x,t);
int result = iv1->FunctionValues( param.values(), y.values() );
int result = iv1->Evaluate( param.values(), y.values() );
if( !result )
throw Undefined();
vtkDataSet *ds = iv1->GetLastDataSet();
vtkIdType cellID = iv1->GetLastCellId();
vtkDataSet *ds = iv1->GetDataSet();
vtkIdType cellID = iv1->GetLastCell();
vtkCell *cell = ds->GetCell( cellID );
vtkDoubleArray *cellVectors;
......@@ -181,9 +191,8 @@ EXCEPTION0(ImproperUseException); // didn't do this.
inVectors->GetTuples( cell->PointIds, cellVectors );
double *cellVel = cellVectors->GetPointer(0);
double pcoords[3], w[100];
iv1->GetLastWeights( w );
iv1->GetLastLocalCoordinates( pcoords );
double *w = iv1->GetLastWeights();
double *pcoords = iv1->GetLastPCoords();
cell->Derivatives( 0, pcoords, cellVel, 3, derivs);
//cout<<"pcoords= "<<pcoords[0]<<" "<<pcoords[1]<<" "<<pcoords[2]<<endl;
......@@ -217,6 +226,8 @@ EXCEPTION0(ImproperUseException); // didn't do this.
// Programmer: Dave Pugmire (on behalf of Hank Childs)
// Creation: Tue Feb 24 09:24:49 EST 2009
//
// Modifications:
//
// Dave Pugmire, Tue Mar 10 12:41:11 EDT 2009
// Check time bounds.
//
......@@ -229,8 +240,8 @@ avtIVPVTKTimeVaryingField::IsInside( const double& t, const avtVecRef& x ) const
avtVec param = pad(x,t);
return (t >= time1 && t <= time2 &&
iv1->FunctionValues(param.values(), y.values()) &&
iv2->FunctionValues(param.values(), y.values()));
iv1->Evaluate(param.values(), y.values()) &&
iv2->Evaluate(param.values(), y.values()));
}
......@@ -243,12 +254,20 @@ avtIVPVTKTimeVaryingField::IsInside( const double& t, const avtVecRef& x ) const
// Programmer: Dave Pugmire (on behalf of Hank Childs)
// Creation: Tue Feb 24 09:24:49 EST 2009
//
// Modifications:
//
// Hank Childs, Thu Apr 2 16:36:50 PDT 2009
// Update to use vtkVisItInterpolatedVelocityField, not
// vtkInterpolatedVelocityField.
// (Hard code "3", since that is what the old implementation was doing
// anyways...)
//
// ****************************************************************************
unsigned int
avtIVPVTKTimeVaryingField::GetDimension() const
{
return iv1->GetNumberOfFunctions();
return 3;
}
// ****************************************************************************
......
......@@ -37,7 +37,7 @@
*****************************************************************************/
// ************************************************************************* //
// avtIVPVTKTimeVaryingField.h //
// avtIVPVTKTimeVaryingField.h //
// ************************************************************************* //
#ifndef AVT_IVPVTKTimeVaryingFIELD_H
......@@ -45,7 +45,7 @@
#include <avtIVPField.h>
#include <vtkInterpolatedVelocityField.h>
#include <vtkVisItInterpolatedVelocityField.h>
#include <ivp_exports.h>
......@@ -55,7 +55,7 @@
//
// Purpose:
// A wrapper class to allow the use of vtkDataSets as IVP fields for
// streamline integration. Uses vtkInterpolatedVelocityField on top of
// streamline integration. Uses vtkVisItInterpolatedVelocityField on top of
// the supplied vtkDataSet.
//
// Programmer: Dave Pugmire (on behalf of Hank Childs)
......@@ -66,13 +66,16 @@
// Dave Pugmire, Tue Mar 10 12:41:11 EDT 2009
// Add GetValidTimeRange.
//
// Hank Childs, Thu Apr 2 16:40:08 PDT 2009
// Use vtkVisItInterpolatedVelocityField.
//
// ****************************************************************************
class IVP_API avtIVPVTKTimeVaryingField: public avtIVPField
{
public:
avtIVPVTKTimeVaryingField(vtkInterpolatedVelocityField *velocity1,
vtkInterpolatedVelocityField *velocity2,
avtIVPVTKTimeVaryingField(vtkVisItInterpolatedVelocityField *velocity1,
vtkVisItInterpolatedVelocityField *velocity2,
double time1, double time2);
~avtIVPVTKTimeVaryingField();
......@@ -88,8 +91,8 @@ class IVP_API avtIVPVTKTimeVaryingField: public avtIVPField
protected:
vtkInterpolatedVelocityField *iv1;
vtkInterpolatedVelocityField *iv2;
vtkVisItInterpolatedVelocityField *iv1;
vtkVisItInterpolatedVelocityField *iv2;
double time1;
double time2;
bool normalized;
......
......@@ -193,6 +193,11 @@ avtStreamline::Advance(const avtIVPField* field,
// Reworked the termination code. Added a type enum and value. Made num steps
// a termination criterion. Code cleanup: We no longer need fwd/bwd solvers.
//
// Hank Childs, Thu Apr 2 16:40:08 PDT 2009
// Fix problem with stalling out during initialization in case where
// seed point is close to boundary of domain. Done in consultation with
// Christoph.
//
// ****************************************************************************
avtIVPSolver::Result
......@@ -232,7 +237,17 @@ avtStreamline::DoAdvance(avtIVPSolver* ivp,
// if step size is below given minimum, give up
// restore old state to before failed step
double hBeforePush = ivp->GetNextStepSize();
ivp->PutState( state );
if (ivp->GetNextStepSize() == 0.)
{
// This can happen if we try to look a few points out
// for the very first step and one of those points
// is outside the domain. Just set the step size
// back to what it was before so we can try again
// with a smaller step.
ivp->SetNextStepSize(hBeforePush);
}
double h = ivp->GetNextStepSize();
......
../../visit_vtk/full/vtkVisItInterpolatedVelocityField.h
\ No newline at end of file
......@@ -310,6 +310,9 @@
# Moved vtkMultiFontVectorText and vtkVectorFontCharacterData to
# rendering lib.
#
# Hank Childs, Thu Apr 2 17:54:09 CDT 2009
# Added vtkVisItInterpolatedVelocityField.
#
##############################################################################
##
......@@ -355,6 +358,7 @@ SRC= vtkAxisDepthSort.C \
vtkVisItExtractRectilinearGrid.C \
vtkVisItFeatureEdges.C \
vtkVisItGlyph3D.C \
vtkVisItInterpolatedVelocityField.C \
vtkVisItProbeFilter.C \
vtkVisItPolyDataNormals.C \
vtkVisItRectilinearGrid.C \
......
/*****************************************************************************
*
* Copyright (c) 2000 - 2008, Lawrence Livermore National Security, LLC
* Produced at the Lawrence Livermore National Laboratory
* LLNL-CODE-400142
* All rights reserved.
*
* This file is part of VisIt. For details, see https://visit.llnl.gov/. The
* full copyright notice is contained in the file COPYRIGHT located at the root
* of the VisIt distribution or at http://www.llnl.gov/visit/copyright.html.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* - Redistributions of source code must retain the above copyright notice,
* this list of conditions and the disclaimer below.
* - Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the disclaimer (as noted below) in the
* documentation and/or other materials provided with the distribution.
* - Neither the name of the LLNS/LLNL nor the names of its contributors may
* be used to endorse or promote products derived from this software without
* specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL LAWRENCE LIVERMORE NATIONAL SECURITY,
* LLC, THE U.S. DEPARTMENT OF ENERGY OR CONTRIBUTORS BE LIABLE FOR ANY
* DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
* SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
* DAMAGE.
*
*****************************************************************************/
#include <vtkVisItInterpolatedVelocityField.h>
#include <float.h>
#include <vtkCell.h>
#include <vtkDataArray.h>
#include <vtkDataSet.h>
#include <vtkGenericCell.h>
#include <vtkIdList.h>
#include <vtkObjectFactory.h>
#include <vtkPointData.h>
#include <vtkRectilinearGrid.h>
#include <vtkVisItPointLocator.h>
#include <vtkVisItUtility.h>
#include <DebugStream.h>
vtkVisItInterpolatedVelocityField* vtkVisItInterpolatedVelocityField::New()
{
// First try to create the object from the vtkObjectFactory
vtkObject* ret = vtkObjectFactory::CreateInstance("vtkVisItInterpolatedVelocityField");
if(ret)
{
return (vtkVisItInterpolatedVelocityField*)ret;
}
// If the factory was unable to create the object, then create it here.
return new vtkVisItInterpolatedVelocityField;
}
vtkVisItInterpolatedVelocityField::vtkVisItInterpolatedVelocityField()
{
ds = NULL;
lastCell = 0;
pcoords[0] = pcoords[1] = pcoords[2] = 0.;
for (int i = 0 ; i < 1024 ; i++)
weights[i] = 0.;
locator = NULL;
}
vtkVisItInterpolatedVelocityField::~vtkVisItInterpolatedVelocityField()
{
if (ds != NULL)
ds->Delete();
if (locator != NULL)
locator->Delete();
}
void
vtkVisItInterpolatedVelocityField::SetDataSet(vtkDataSet *ds_)
{
if (ds != NULL)
ds->Delete();
if (locator != NULL)
{
locator->Delete();
locator = NULL;
}
ds = ds_;
ds->Register(NULL);
}
bool
vtkVisItInterpolatedVelocityField::Evaluate(double *pt, double *vel)
{
if (ds == NULL)
{
debug1 <<" vtkVisItInterpolatedVelocityField::No data set to evaluate!" << endl;
return false;
}
vtkDataArray *vectors = ds->GetPointData()->GetVectors();