vtkCellLocatorInterpolatedVelocityField.cxx

/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkCellLocatorInterpolatedVelocityField.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 "vtkCellLocatorInterpolatedVelocityField.h"

#include "vtkMath.h"
#include "vtkDataSet.h"
#include "vtkDataArray.h"
#include "vtkPointData.h"
#include "vtkGenericCell.h"
#include "vtkCellLocator.h"
#include "vtkSmartPointer.h"
#include "vtkObjectFactory.h"
#include "vtkModifiedBSPTree.h"

vtkStandardNewMacro ( vtkCellLocatorInterpolatedVelocityField ); 
vtkCxxSetObjectMacro( vtkCellLocatorInterpolatedVelocityField, CellLocatorPrototype, vtkAbstractCellLocator );

//----------------------------------------------------------------------------
typedef vtkstd::vector< vtkSmartPointer < vtkAbstractCellLocator > > CellLocatorsTypeBase;
class vtkCellLocatorInterpolatedVelocityFieldCellLocatorsType : public CellLocatorsTypeBase { };

//----------------------------------------------------------------------------
vtkCellLocatorInterpolatedVelocityField::vtkCellLocatorInterpolatedVelocityField()
{   
  this->LastCellLocator  = 0;
  this->CellLocatorPrototype = 0;
  this->CellLocators = new vtkCellLocatorInterpolatedVelocityFieldCellLocatorsType;
}

//----------------------------------------------------------------------------
vtkCellLocatorInterpolatedVelocityField::~vtkCellLocatorInterpolatedVelocityField()
{   
  this->LastCellLocator = 0;
  this->SetCellLocatorPrototype( 0 );
  
  if ( this->CellLocators )
    {
    delete this->CellLocators; 
    this->CellLocators = 0; 
    }
}

//----------------------------------------------------------------------------
void vtkCellLocatorInterpolatedVelocityField::SetLastCellId
  ( vtkIdType c, int dataindex ) 
{
  this->LastCellId       = c; 
  this->LastDataSet      = ( *this->DataSets )[dataindex];
  this->LastCellLocator  = ( *this->CellLocators )[dataindex].GetPointer();
  this->LastDataSetIndex = dataindex;
  
  // If the dataset changes, then the cached cell is invalidated. We might as 
  // well prefetch the cached cell either way.
  if ( this->LastCellId != -1 )
  {
    this->LastDataSet->GetCell( this->LastCellId, this->GenCell );
  } 
}

//----------------------------------------------------------------------------
int vtkCellLocatorInterpolatedVelocityField::FunctionValues
  ( double * x, double * f )
{
  vtkDataSet *             vds = NULL;
  vtkAbstractCellLocator * loc = NULL;

  if( !this->LastDataSet && !this->DataSets->empty() )
    {
    vds = ( *this->DataSets )[0];
    loc = ( *this->CellLocators )[0].GetPointer();
    this->LastDataSet      = vds;
    this->LastCellLocator  = loc;
    this->LastDataSetIndex = 0;
    }
  else
    {
    vds = this->LastDataSet;
    loc = this->LastCellLocator;
    }
  
  int retVal;
  if ( loc )
    {
    // resort to vtkAbstractCellLocator::FindCell()
    retVal = this->FunctionValues( vds, loc, x, f );
    }
  else
    {
    // turn to vtkImageData/vtkRectilinearGrid::FindCell()
    retVal = this->FunctionValues( vds, x, f );
    }
  
  if ( !retVal )
    {
    for( this->LastDataSetIndex = 0; 
         this->LastDataSetIndex < static_cast<int>( this->DataSets->size() );
         this->LastDataSetIndex ++ )
      {
      vds = this->DataSets->operator[]( this->LastDataSetIndex );
      loc = this->CellLocators->operator[]( this->LastDataSetIndex ).GetPointer();
      if( vds && vds != this->LastDataSet )
        {
        this->ClearLastCellId();
        
        if ( loc )
          {
          // resort to vtkAbstractCellLocator::FindCell()
          retVal = this->FunctionValues( vds, loc, x, f );
          }
        else
          {
          // turn to vtkImageData/vtkRectilinearGrid::FindCell()
          retVal = this->FunctionValues( vds, x, f );
          }

        if ( retVal ) 
          {
          this->LastDataSet     = vds;
          this->LastCellLocator = loc;
          vds = NULL;
          loc = NULL;
          return retVal;
          }
        }
      }
      
    this->LastCellId       = -1;
    this->LastDataSet      = ( *this->DataSets )[0];
    this->LastCellLocator  = ( *this->CellLocators )[0].GetPointer();
    this->LastDataSetIndex = 0;
    vds = NULL;
    loc = NULL;
    return 0;
    }
  
  vds = NULL;
  loc = NULL;
  return retVal;
}

//----------------------------------------------------------------------------
int vtkCellLocatorInterpolatedVelocityField::FunctionValues
  ( vtkDataSet * dataset, vtkAbstractCellLocator * loc, double * x, double * f )
{ 
  f[0] = f[1] = f[2] = 0.0;
  vtkDataArray * vectors = NULL;
  
  if ( !dataset || !loc || !dataset->IsA( "vtkPointSet" ) ||
       !( vectors = dataset->GetPointData()
                           ->GetVectors( this->VectorsSelection )
        )
     )
    {
    vtkErrorMacro( <<"Can't evaluate dataset!" );
    vectors = NULL;
    return  0;
    }
    
  int    i;
  int    subIdx;
  int    numPts;
  int    pntIdx;
  int    bFound = 0;
  double vector[3];
  double dstns2 = 0.0;
  double toler2 = dataset->GetLength() * 
                  vtkCellLocatorInterpolatedVelocityField::TOLERANCE_SCALE;

  // check if the point is in the cached cell AND can be successfully evaluated
  if ( this->LastCellId != -1 &&
       this->GenCell->EvaluatePosition
             ( x, 0, subIdx, this->LastPCoords, dstns2, this->Weights ) == 1
     ) 
    {
    bFound = 1;
    this->CacheHit ++;
    }
  
  if ( !bFound )
    {
    // cache missing or evaluation failure and then we have to find the cell
    this->CacheMiss += !(  !( this->LastCellId + 1 )  );
    this->LastCellId = loc->FindCell( x, toler2, this->GenCell, 
                                      this->LastPCoords, this->Weights );
    bFound = !(  !( this->LastCellId + 1 )  );
    }
  
  // interpolate vectors if possible                              
  if ( bFound )
    {
    numPts = this->GenCell->GetNumberOfPoints();
    for ( i = 0; i < numPts; i ++ )
      {
      pntIdx = this->GenCell->PointIds->GetId( i );
      vectors->GetTuple( pntIdx, vector );
      f[0] += vector[0] * this->Weights[i];
      f[1] += vector[1] * this->Weights[i];
      f[2] += vector[2] * this->Weights[i];
      }
      
    if ( this->NormalizeVector == true )
      {
      vtkMath::Normalize( f );
      }  
    }
   
  vectors = NULL;
  return  bFound;
}

//----------------------------------------------------------------------------
void vtkCellLocatorInterpolatedVelocityField::AddDataSet( vtkDataSet * dataset )
{
  if ( !dataset )
    {
    vtkErrorMacro( <<"Dataset NULL!" );
    return;
    }

  // insert the dataset (do NOT register the dataset to 'this')
  this->DataSets->push_back( dataset );
  
  // We need to attach a valid vtkAbstractCellLocator to any vtkPointSet for
  // robust cell location as vtkPointSet::FindCell() may incur failures. For 
  // any non-vtkPointSet dataset, either vtkImageData or vtkRectilinearGrid, 
  // we do not need to associate a vtkAbstractCellLocator with it (though a 
  // NULL vtkAbstractCellLocator is still inserted to this->CellLocators to
  // enable proper access to those valid cell locators) since these two kinds
  // of datasets themselves are able to guarantee robust as well as fast cell 
  // location via vtkImageData/vtkRectilinearGrid::FindCell().
  vtkSmartPointer< vtkAbstractCellLocator > locator = 0; // MUST inited with 0
  if (  dataset->IsA( "vtkPointSet" )  )
    {
     
    if ( !this->CellLocatorPrototype )
      {
      locator = vtkSmartPointer < vtkModifiedBSPTree >::New();
      }
    else
      {
      locator.TakeReference( this->CellLocatorPrototype->NewInstance() );
      }
      
    locator->SetLazyEvaluation( 1 );
    locator->SetDataSet( dataset );
    }
  this->CellLocators->push_back( locator );
  
  int  size = dataset->GetMaxCellSize();
  if ( size > this->WeightsSize )
    {
    this->WeightsSize = size;
    if ( this->Weights )
      {
      delete[] this->Weights;
      this->Weights = NULL;
      }
    this->Weights = new double[size]; 
    }
}

//----------------------------------------------------------------------------
void vtkCellLocatorInterpolatedVelocityField::CopyParameters
  ( vtkAbstractInterpolatedVelocityField * from )
{
  this->Superclass::CopyParameters( from );
  
  if (  from->IsA( "vtkCellLocatorInterpolatedVelocityField" )  )
    {
    this->SetCellLocatorPrototype
          (  vtkCellLocatorInterpolatedVelocityField::SafeDownCast( from )
             ->GetCellLocatorPrototype()
          );
    }
}

//----------------------------------------------------------------------------
void vtkCellLocatorInterpolatedVelocityField::PrintSelf( ostream & os, vtkIndent indent )
{
  this->Superclass::PrintSelf( os, indent );
 
  os << indent << "CellLocators: "     << this->CellLocators     << endl;
  if ( this->CellLocators )
    {
    os << indent << "Number of Cell Locators: " << this->CellLocators->size();
    }
  os << indent << "LastCellLocator: "      << this->LastCellLocator      << endl;
  os << indent << "CellLocatorPrototype: " << this->CellLocatorPrototype << endl;
}