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

  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.

=========================================================================*/
// VTK_DEPRECATED_IN_9_2_0() warnings for this class.
#define VTK_DEPRECATION_LEVEL 0

#include "vtkCachingInterpolatedVelocityField.h"

#include "vtkCellLocator.h"
#include "vtkDataSet.h"
#include "vtkDoubleArray.h"
#include "vtkFloatArray.h"
#include "vtkGenericCell.h"
#include "vtkMath.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkSmartPointer.h"

#include <vector>

//------------------------------------------------------------------------------
VTK_ABI_NAMESPACE_BEGIN
vtkStandardNewMacro(vtkCachingInterpolatedVelocityField);
//------------------------------------------------------------------------------
const double IVFDataSetInfo::TOLERANCE_SCALE = 1.0E-8;
//------------------------------------------------------------------------------
IVFDataSetInfo::IVFDataSetInfo()
{
  this->VelocityFloat = nullptr;
  this->VelocityDouble = nullptr;
  this->DataSet = nullptr;
  this->Cell = nullptr;
  this->BSPTree = nullptr;
  this->Tolerance = 0.0;
  this->StaticDataSet = false;
}
//------------------------------------------------------------------------------
void IVFDataSetInfo::SetDataSet(
  vtkDataSet* data, char* velocity, bool staticdataset, vtkAbstractCellLocator* locator)
{
  this->VelocityFloat = nullptr;
  this->VelocityDouble = nullptr;
  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<vtkCellLocator>::New();
    }
    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 (vtkArrayDownCast<vtkFloatArray>(vectors))
  {
    this->VelocityFloat = vtkArrayDownCast<vtkFloatArray>(vectors)->GetPointer(0);
  }
  else if (vtkArrayDownCast<vtkDoubleArray>(vectors))
  {
    this->VelocityDouble = vtkArrayDownCast<vtkDoubleArray>(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;
  this->PCoords[0] = ivfci.PCoords[0];
  this->PCoords[1] = ivfci.PCoords[1];
  this->PCoords[2] = ivfci.PCoords[2];
}
//------------------------------------------------------------------------------
IVFDataSetInfo& IVFDataSetInfo::operator=(const IVFDataSetInfo& ivfci)
{
  if (this == &ivfci)
  {
    return *this;
  }

  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;
  this->PCoords[0] = ivfci.PCoords[0];
  this->PCoords[1] = ivfci.PCoords[1];
  this->PCoords[2] = ivfci.PCoords[2];
  return *this;
}
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
vtkCachingInterpolatedVelocityField::vtkCachingInterpolatedVelocityField()
{
  this->NumFuncs = 3;     // u, v, w
  this->NumIndepVars = 4; // x, y, z, t
  this->VectorsSelection = nullptr;
  this->TempCell = vtkGenericCell::New();
  this->CellCacheHit = 0;
  this->DataSetCacheHit = 0;
  this->CacheMiss = 0;
  this->LastCacheIndex = 0;
  this->Cache = nullptr;
  this->LastCellId = -1;
}
//------------------------------------------------------------------------------
vtkCachingInterpolatedVelocityField::~vtkCachingInterpolatedVelocityField()
{
  this->NumFuncs = 0;
  this->NumIndepVars = 0;
  this->TempCell->Delete();
  this->SetVectorsSelection(nullptr);
}
//------------------------------------------------------------------------------
void vtkCachingInterpolatedVelocityField::SetDataSet(
  int I, vtkDataSet* dataset, bool staticdataset, vtkAbstractCellLocator* locator)
{
  int N = vtkMath::Max(I + 1, static_cast<int>(this->CacheList.size()));
  this->CacheList.resize(N);
  this->CacheList[I].SetDataSet(dataset, this->VectorsSelection, staticdataset, locator);
  //
  int maxsize = vtkMath::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))
  {
    assert(this->CacheList.size() > static_cast<size_t>(datasetindex));
    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)
    {
      assert(c < this->Cache->DataSet->GetNumberOfCells());
      this->Cache->DataSet->GetCell(this->LastCellId, this->Cache->Cell);
    }
  }
}
//------------------------------------------------------------------------------
void vtkCachingInterpolatedVelocityField::ClearLastCellInfo()
{
  this->Cache = nullptr;
  this->LastCellId = -1;
}
//------------------------------------------------------------------------------
vtkGenericCell* vtkCachingInterpolatedVelocityField::GetLastCell()
{
  if (this->Cache)
  {
    return this->Cache->Cell;
  }
  return nullptr;
}
//------------------------------------------------------------------------------
// 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)
    {
      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, nullptr, subId, this->Cache->PCoords, this->Cache->Tolerance, this->Weights.data()) == 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)
    {
      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.data());
  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, nullptr, subId, data->PCoords, dist2, this->Weights.data()) ==
        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.data());
    this->LastCellId = cellId;
  }
  else
  {
    vtkGenericCell* tmpCell = nullptr;
    if (this->LastCellId >= 0)
    {
      data->DataSet->GetCell(this->LastCellId, this->TempCell);
      tmpCell = this->TempCell;
    }
    this->LastCellId = data->DataSet->FindCell(x, tmpCell, data->Cell, this->LastCellId,
      data->Tolerance, subId, data->PCoords, this->Weights.data());
    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++)
    {
      vtkIdType 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++)
    {
      vtkIdType 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 false;
  }
  outPD->InterpolatePoint(this->Cache->DataSet->GetPointData(), outIndex,
    this->Cache->Cell->PointIds, this->Weights.data());
  return true;
}
//------------------------------------------------------------------------------
bool vtkCachingInterpolatedVelocityField::InterpolatePoint(
  vtkCachingInterpolatedVelocityField* inCIVF, vtkPointData* outPD, vtkIdType outIndex)
{
  if (!this->Cache || !this->Cache->DataSet)
  {
    return false;
  }
  vtkPointData* inPD = inCIVF->Cache->DataSet->GetPointData();
  outPD->InterpolatePoint(inPD, outIndex, this->Cache->Cell->PointIds, this->Weights.data());
  return true;
}
//------------------------------------------------------------------------------
int vtkCachingInterpolatedVelocityField::GetLastWeights(double* w)
{
  // If last cell is valid, fill w with the interpolation weights
  // and return true
  if (this->Cache && this->LastCellId >= 0)
  {
    int numPts = this->Cache->Cell->GetNumberOfPoints();
    for (int j = 0; j < numPts; j++)
    {
      w[j] = this->Weights[j];
    }
    return 1;
  }
  // otherwise, return false
  else
  {
    return 0;
  }
}
//------------------------------------------------------------------------------
int vtkCachingInterpolatedVelocityField::GetLastLocalCoordinates(double pcoords[3])
{
  // If last cell is valid, fill p with the local coordinates
  // and return true
  if (this->Cache && this->LastCellId >= 0)
  {
    for (int 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);

  if (!Weights.empty())
  {
    os << indent << "Weights: " << this->Weights.data() << endl;
  }
  else
  {
    os << indent << "Weights: (none)" << 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;

  if (this->Cache)
  {
    os << indent << "Cache->DataSet : " << this->Cache->DataSet << endl;
  }
  else
  {
    os << indent << "Cache->DataSet : (none)" << endl;
  }

  os << indent << "LastCacheIndex : " << this->LastCacheIndex << endl;
}
//------------------------------------------------------------------------------
VTK_ABI_NAMESPACE_END
