// SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
// SPDX-License-Identifier: BSD-3-Clause

#include "vtkAMRResampleFilter.h"
#include "vtkAMRUtilities.h"
#include "vtkBoundingBox.h"
#include "vtkCell.h"
#include "vtkCellData.h"
#include "vtkCompositeDataPipeline.h"
#include "vtkDataArray.h"
#include "vtkExtentRCBPartitioner.h"
#include "vtkFieldData.h"
#include "vtkIndent.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkMath.h"
#include "vtkMultiBlockDataSet.h"
#include "vtkMultiProcessController.h"
#include "vtkObjectFactory.h"
#include "vtkOverlappingAMR.h"
#include "vtkOverlappingAMRMetaData.h"
#include "vtkPointData.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkUniformGrid.h"
#include "vtkUniformGridPartitioner.h"

#include "vtkTimerLog.h"

#include <algorithm>
#include <cassert>
#include <cmath>
#include <iostream>
#include <sstream>

VTK_ABI_NAMESPACE_BEGIN
vtkStandardNewMacro(vtkAMRResampleFilter);
vtkCxxSetObjectMacro(vtkAMRResampleFilter, Controller, vtkMultiProcessController);

//------------------------------------------------------------------------------
vtkAMRResampleFilter::vtkAMRResampleFilter()
{
  this->TransferToNodes = 1;
  this->DemandDrivenMode = 0;
  this->NumberOfPartitions = 1;
  this->LevelOfResolution = 0;
  this->AMRMetaData = nullptr;
  this->NumberOfSamples[0] = this->NumberOfSamples[1] = this->NumberOfSamples[2] = 10;
  this->Controller = nullptr;
  this->SetController(vtkMultiProcessController::GetGlobalController());
  this->ROI = vtkMultiBlockDataSet::New();

  for (int i = 0; i < 3; ++i)
  {
    this->Min[i] = 0.;
    this->Max[i] = 1.;
  }
  this->SetNumberOfInputPorts(1);
  this->SetNumberOfOutputPorts(1);
  this->UseBiasVector = false;
  this->BiasVector[0] = this->BiasVector[1] = this->BiasVector[2] = 0.0;
}

//------------------------------------------------------------------------------
vtkAMRResampleFilter::~vtkAMRResampleFilter()
{
  this->BlocksToLoad.clear();
  this->SetController(nullptr);

  if (this->ROI != nullptr)
  {
    this->ROI->Delete();
  }
  this->ROI = nullptr;

  //  if( this->AMRMetaData != nullptr )
  //    {
  //    this->AMRMetaData->Delete();
  //    }
  //  this->AMRMetaData = nullptr;
}

//------------------------------------------------------------------------------
void vtkAMRResampleFilter::PrintSelf(std::ostream& oss, vtkIndent indent)
{
  this->Superclass::PrintSelf(oss, indent);
}

//------------------------------------------------------------------------------
int vtkAMRResampleFilter::FillInputPortInformation(int vtkNotUsed(port), vtkInformation* info)
{
  assert("pre: information object is nullptr" && (info != nullptr));
  info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkOverlappingAMR");
  return 1;
}

//------------------------------------------------------------------------------
int vtkAMRResampleFilter::FillOutputPortInformation(int vtkNotUsed(port), vtkInformation* info)
{
  assert("pre: information object is nullptr" && (info != nullptr));
  info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkMultiBlockDataSet");
  return 1;
}

//------------------------------------------------------------------------------
int vtkAMRResampleFilter::RequestUpdateExtent(vtkInformation*, vtkInformationVector** inputVector,
  vtkInformationVector* vtkNotUsed(outputVector))
{
  assert("pre: inputVector is nullptr" && (inputVector != nullptr));
  vtkInformation* info = inputVector[0]->GetInformationObject(0);
  assert("pre: info is nullptr" && (info != nullptr));

  if (this->DemandDrivenMode == 1)
  {
    // Tell reader to load all requested blocks.
    info->Set(vtkCompositeDataPipeline::LOAD_REQUESTED_BLOCKS(), 1);

    // Tell reader which blocks this process requires
    info->Set(vtkCompositeDataPipeline::UPDATE_COMPOSITE_INDICES(), this->BlocksToLoad.data(),
      static_cast<int>(this->BlocksToLoad.size()));
  }
  return 1;
}

//------------------------------------------------------------------------------
int vtkAMRResampleFilter::RequestInformation(vtkInformation* vtkNotUsed(rqst),
  vtkInformationVector** inputVector, vtkInformationVector* vtkNotUsed(outputVector))
{

  assert("pre: inputVector is nullptr" && (inputVector != nullptr));

  vtkInformation* input = inputVector[0]->GetInformationObject(0);
  assert("pre: input is nullptr" && (input != nullptr));

  if (this->DemandDrivenMode == 1 &&
    input->Has(vtkCompositeDataPipeline::COMPOSITE_DATA_META_DATA()))
  {
    //    this->AMRMetaData = vtkOverlappingAMR::New();
    this->AMRMetaData = vtkOverlappingAMR::SafeDownCast(
      input->Get(vtkCompositeDataPipeline::COMPOSITE_DATA_META_DATA()));

    // Get Region
    double h[3];
    this->ComputeAndAdjustRegionParameters(this->AMRMetaData, h);
    this->GetRegion(h);

    // Compute which blocks to load
    this->ComputeAMRBlocksToLoad(this->AMRMetaData);
  }

  // Don't we need to call this->Modified() here?
  // this->Modified();
  return 1;
}

//------------------------------------------------------------------------------
int vtkAMRResampleFilter::RequestData(vtkInformation* vtkNotUsed(rqst),
  vtkInformationVector** inputVector, vtkInformationVector* outputVector)
{
  std::cerr << "Running Resampler\n";

  // STEP 0: Get input object
  vtkInformation* input = inputVector[0]->GetInformationObject(0);
  assert("pre: Null information object!" && (input != nullptr));
  vtkOverlappingAMR* amrds =
    vtkOverlappingAMR::SafeDownCast(input->Get(vtkDataObject::DATA_OBJECT()));
  assert("pre: input AMR dataset is nullptr" && (amrds != nullptr));

  // STEP 1: Get output object
  vtkInformation* output = outputVector->GetInformationObject(0);
  assert("pre: Null output information object!" && (output != nullptr));

  vtkMultiBlockDataSet* mbds =
    vtkMultiBlockDataSet::SafeDownCast(output->Get(vtkDataObject::DATA_OBJECT()));
  assert("pre: output grid is nullptr" && (mbds != nullptr));

  // STEP 2: Get Metadata
  if (this->DemandDrivenMode == 1)
  {
    assert("pre: Metadata must have been populated in RqstInfo" && (this->AMRMetaData != nullptr));
    this->ExtractRegion(amrds, mbds, this->AMRMetaData);
  }
  else
  {
    // GetRegion
    double h[3];
    this->ComputeAndAdjustRegionParameters(amrds, h);
    this->GetRegion(h);
    this->ExtractRegion(amrds, mbds, amrds);
  }

  return 1;
}

//------------------------------------------------------------------------------
bool vtkAMRResampleFilter::FoundDonor(double q[3], vtkUniformGrid*& donorGrid, int& cellIdx)
{
  assert("pre: donor grid is nullptr" && (donorGrid != nullptr));
  double gbounds[6];
  // Lets do a trivial spatial check
  this->NumberOfBlocksTested++;
  donorGrid->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;
  }
  int ijk[3];
  double pcoords[3];
  int status = donorGrid->ComputeStructuredCoordinates(q, ijk, pcoords);
  if (status == 1)
  {
    cellIdx = vtkStructuredData::ComputeCellId(donorGrid->GetDimensions(), ijk);
    return true;
  }
  return false;
}

//------------------------------------------------------------------------------
void vtkAMRResampleFilter::InitializeFields(vtkFieldData* f, vtkIdType size, vtkCellData* src)
{
  assert("pre: field data is nullptr!" && (f != nullptr));
  assert("pre: source cell data is nullptr" && (src != nullptr));

  for (int arrayIdx = 0; arrayIdx < src->GetNumberOfArrays(); ++arrayIdx)
  {
    int dataType = src->GetArray(arrayIdx)->GetDataType();
    vtkDataArray* array = vtkDataArray::CreateDataArray(dataType);
    assert("pre: failed to create array!" && (array != nullptr));

    array->SetName(src->GetArray(arrayIdx)->GetName());
    array->SetNumberOfComponents(src->GetArray(arrayIdx)->GetNumberOfComponents());
    array->SetNumberOfTuples(size);
    assert("post: array size mismatch" && (array->GetNumberOfTuples() == size));

    f->AddArray(array);
    array->Delete();

    //      int myIndex = -1;
    //      vtkDataArray *arrayPtr = f->GetArray(
    //          src->GetArray(arrayIdx)->GetName(), myIndex );
    //      assert( "post: array index mismatch" && (myIndex == arrayIdx) );
    //      assert( "post: array size mismatch" &&
    //              (arrayPtr->GetNumberOfTuples()==size) );

    assert("post: array size mismatch" && (f->GetArray(arrayIdx)->GetNumberOfTuples() == size));
  } // END for all arrays
}

//------------------------------------------------------------------------------
void vtkAMRResampleFilter::CopyData(
  vtkFieldData* target, vtkIdType targetIdx, vtkCellData* src, vtkIdType srcIdx)
{
  assert("pre: target field data is nullptr" && (target != nullptr));
  assert("pre: source field data is nullptr" && (src != nullptr));
  assert("pre: number of arrays does not match" &&
    (target->GetNumberOfArrays() == src->GetNumberOfArrays()));

  int arrayIdx = 0;
  for (; arrayIdx < src->GetNumberOfArrays(); ++arrayIdx)
  {
    vtkDataArray* targetArray = target->GetArray(arrayIdx);
    vtkDataArray* srcArray = src->GetArray(arrayIdx);
    assert("pre: target array is nullptr!" && (targetArray != nullptr));
    assert("pre: source array is nullptr!" && (srcArray != nullptr));
    assert("pre: target/source array number of components mismatch!" &&
      (targetArray->GetNumberOfComponents() == srcArray->GetNumberOfComponents()));
    assert("pre: target/source array names mismatch!" &&
      (strcmp(targetArray->GetName(), srcArray->GetName()) == 0));
    assert("pre: source index is out-of-bounds" && (srcIdx >= 0) &&
      (srcIdx < srcArray->GetNumberOfTuples()));
    assert("pre: target index is out-of-bounds" && (targetIdx >= 0) &&
      (targetIdx < targetArray->GetNumberOfTuples()));

    int c = 0;
    for (; c < srcArray->GetNumberOfComponents(); ++c)
    {
      double f = srcArray->GetComponent(srcIdx, c);
      targetArray->SetComponent(targetIdx, c, f);
    } // END for all components

  } // END for all arrays
}

//------------------------------------------------------------------------------
void vtkAMRResampleFilter::ComputeCellCentroid(vtkUniformGrid* g, vtkIdType cellIdx, double c[3])
{
  assert("pre: uniform grid is nullptr" && (g != nullptr));
  assert("pre: centroid is nullptr" && (c != nullptr));
  assert("pre: cell index out-of-bounds" && (cellIdx >= 0) && (cellIdx < g->GetNumberOfCells()));

  vtkCell* myCell = g->GetCell(cellIdx);
  assert("post: cell is nullptr!" && (myCell != nullptr));

  double pc[3]; // the parametric center
  double* weights = new double[myCell->GetNumberOfPoints()];
  assert("post: weights vector is nullptr" && (weights != nullptr));

  int subId = myCell->GetParametricCenter(pc);
  myCell->EvaluateLocation(subId, pc, c, weights);
  delete[] weights;
}

//------------------------------------------------------------------------------
void vtkAMRResampleFilter::TransferToCellCenters(vtkUniformGrid* g, vtkOverlappingAMR* amrds)
{
  assert("pre: uniform grid is nullptr" && (g != nullptr));
  assert("pre: AMR data-structure is nullptr" && (amrds != nullptr));

  // STEP 0: Get the first block so that we know the arrays
  vtkUniformGrid* refGrid = this->GetReferenceGrid(amrds);

  // STEP 1: Get the cell-data of the reference grid
  vtkCellData* CD = refGrid->GetCellData();
  assert("pre: Donor CellData is nullptr!" && (CD != nullptr));

  // STEP 2: Get the cell data of the resampled grid
  vtkCellData* fieldData = g->GetCellData();
  assert("pre: Target PointData is nullptr!" && (fieldData != nullptr));

  // STEP 3: Initialize the fields on the resampled grid
  this->InitializeFields(fieldData, g->GetNumberOfCells(), CD);

  if (fieldData->GetNumberOfArrays() == 0)
  {
    return;
  }

  // TODO: this is a very naive implementation and should be optimized. However,
  // mostly this filter is used to transfer the solution to the grid nodes and
  // not on the cell nodes.
  vtkIdType cellIdx = 0;
  for (; cellIdx < g->GetNumberOfCells(); ++cellIdx)
  {
    double qPoint[3];
    this->ComputeCellCentroid(g, cellIdx, qPoint);

    unsigned int level = 0;
    for (; level < amrds->GetNumberOfBlocks(level); ++level)
    {
      unsigned int dataIdx = 0;
      for (; dataIdx < amrds->GetNumberOfBlocks(level); ++dataIdx)
      {
        int donorCellIdx = -1;
        vtkUniformGrid* donorGrid =
          vtkUniformGrid::SafeDownCast(amrds->GetDataSetAsCartesianGrid(level, dataIdx));
        if ((donorGrid != nullptr) && this->FoundDonor(qPoint, donorGrid, donorCellIdx))
        {
          assert("pre: donorCellIdx is invalid" && (donorCellIdx >= 0) &&
            (donorCellIdx < donorGrid->GetNumberOfCells()));
          CD = donorGrid->GetCellData();
          this->CopyData(fieldData, cellIdx, CD, donorCellIdx);
        } // END if
      }   // END for all datasets
    }     // END for all levels
  }       // END for all cells
}

//------------------------------------------------------------------------------
bool vtkAMRResampleFilter::SearchForDonorGridAtLevel(double q[3], vtkOverlappingAMR* amrds,
  unsigned int level, unsigned int& donorGridId, int& donorCellIdx)
{
  assert("pre: AMR dataset is nullptr" && (amrds != nullptr));
  this->NumberOfBlocksTestedForLevel = 0;
  std::ostringstream oss;
  oss << "SearchLevel-" << level;

  vtkTimerLog::MarkStartEvent(oss.str().c_str());

  for (donorGridId = 0; donorGridId < amrds->GetNumberOfBlocks(level); ++donorGridId)
  {
    donorCellIdx = -1;
    this->NumberOfBlocksTestedForLevel++;
    if (amrds->GetOverlappingAMRMetaData()->FindCell(q, level, donorGridId, donorCellIdx))
    {
      assert("pre: donorCellIdx is invalid" && (donorCellIdx >= 0)); //  &&
      // (donorCellIdx < donorGrid->GetNumberOfCells()) );
      vtkTimerLog::MarkEndEvent(oss.str().c_str());
      return true;
    } // END if

  } // END for all data at level

  // No suitable grid is found at the requested level, set donorGrid to nullptr
  // to indicate that to the caller.
  vtkTimerLog::MarkEndEvent(oss.str().c_str());
  return false;
}

//------------------------------------------------------------------------------
int vtkAMRResampleFilter::ProbeGridPointInAMR(double q[3], unsigned int& donorLevel,
  unsigned int& donorGridId, vtkOverlappingAMR* amrds, unsigned int maxLevel, bool hadDonorGrid)
{
  assert("pre: AMR dataset is nullptr" && amrds != nullptr);

  vtkUniformGrid* currentGrid = nullptr;
  int currentCellIdx = -1;
  int donorCellIdx = -1;
  unsigned int currentLevel = 0;
  unsigned int currentGridId = 0;
  vtkUniformGrid* donorGrid = hadDonorGrid
    ? vtkUniformGrid::SafeDownCast(amrds->GetDataSetAsCartesianGrid(donorLevel, donorGridId))
    : nullptr;

  // STEP 0: Check the previously cached donor-grid
  if (hadDonorGrid)
  {
    this->NumberOfBlocksTested++;
    bool res(true);
    if (!amrds->GetOverlappingAMRMetaData()->FindCell(q, donorLevel, donorGridId, donorCellIdx))
    {
      // Lets see if the point is contained by a grid at the same donar level
      res = this->SearchForDonorGridAtLevel(q, amrds, donorLevel, donorGridId, donorCellIdx);
      donorGrid = res
        ? vtkUniformGrid::SafeDownCast(amrds->GetDataSetAsCartesianGrid(donorLevel, donorGridId))
        : nullptr;
      this->NumberOfBlocksTested += this->NumberOfBlocksTestedForLevel;
    }

    // If donorGrid is still not nullptr then we found the grid and potential starting
    // level
    if (res)
    {
      assert("pre: donorCellIdx is invalid" && (donorCellIdx >= 0) &&
        (donorCellIdx < donorGrid->GetNumberOfCells()));

      this->NumberOfTimesFoundOnDonorLevel++;

      // Initialize values for step 1 s.t. that the search will start from the
      // current donorLevel
      currentGrid = donorGrid;
      currentGridId = donorGridId;
      currentCellIdx = donorCellIdx;
      currentLevel = donorLevel;
      assert(!donorGrid ||
        vtkUniformGrid::SafeDownCast(amrds->GetDataSetAsCartesianGrid(donorLevel, donorGridId)) ==
          donorGrid);
    }
    else if (donorLevel == 0)
    {
      // if we are here then the point is not contained in any of the level 0
      // blocks!
      this->NumberOfFailedPoints++;
      donorGrid = nullptr;
      donorLevel = 0;
      return -1;
    }
    else
    {
      // If we are here then we know the point is not on the donor level
      // and therefore not contained in any of the more refined levels -
      // Base on the assumption of overlapping AMR
      assert("pre:Donor Level is 0" && donorLevel != 0);
      // Initialize values for step 1 s.t. the search will start from level 0.
      donorGrid = nullptr;
      maxLevel = donorLevel;
      donorLevel = 0;
      currentLevel = 0;
    }
  }

  // If we didn't have an initial donor grid or if we still have one
  // we need to test higher res grids
  int startLevel, endLevel;
  int incLevel;
  if (!((donorGrid == nullptr) && hadDonorGrid))
  {
    startLevel = donorGrid == nullptr ? currentLevel : currentLevel + 1;
    endLevel = maxLevel;
    incLevel = 1;
  }
  else
  {
    startLevel = maxLevel - 1;
    endLevel = -1;
    incLevel = -1;
  }
  // STEP 1: Search in the AMR hierarchy for the donor-grid
  for (int level = startLevel; level != endLevel; level += incLevel)
  {
    if (incLevel == 1)
    {
      this->NumberOfTimesLevelUp++;
    }
    else
    {
      this->NumberOfTimesLevelDown++;
    }
    bool res = this->SearchForDonorGridAtLevel(q, amrds, level, donorGridId, donorCellIdx);
    donorGrid = res
      ? vtkUniformGrid::SafeDownCast(amrds->GetDataSetAsCartesianGrid(level, donorGridId))
      : nullptr;

    this->NumberOfBlocksTested += this->NumberOfBlocksTestedForLevel;
    if (res)
    {
      donorLevel = level;
      // if we are going from fine to coarse then we can stop the search
      if (incLevel == -1)
      {
        assert(vtkUniformGrid::SafeDownCast(
                 amrds->GetDataSetAsCartesianGrid(donorLevel, donorGridId)) == donorGrid);
        return donorCellIdx;
      }

      // Lets see if this is the highest resolution grid that contains the
      // point
      if (donorGrid->IsCellVisible(donorCellIdx))
      {
        // return donorCellIdx;
      }
      // we found a grid that contains the point at level l, let's store it
      // here temporatily in case there is a grid at a higher resolution that
      // we need to use.
      currentGrid = donorGrid;
      currentCellIdx = donorCellIdx;
      currentLevel = level;
      currentGridId = donorGridId;
    }
    else if (currentGrid != nullptr)
    {
      // we did not find the point at a higher res, but, we did find at a lower
      // resolution, so we will use the solution we found previously
      // THIS SHOULD NOW NOT HAPPEN!!
      // vtkErrorMacro("Could not find point in an unblanked cell.");
      this->NumberOfBlocksVisSkipped += this->NumberOfBlocksTestedForLevel;
      donorGrid = currentGrid;
      donorCellIdx = currentCellIdx;
      donorLevel = currentLevel;
      donorGridId = currentGridId;
      assert(!donorGrid ||
        vtkUniformGrid::SafeDownCast(amrds->GetDataSetAsCartesianGrid(donorLevel, donorGridId)) ==
          donorGrid);
      break;
    }
    else
    {
      // we are not able to find a grid/cell that contains the query point, in
      // this case we will just return.
      this->NumberOfFailedPoints++;
      donorCellIdx = -1;
      donorGrid = nullptr;
      donorLevel = 0;
      break;
    }
  } // END for all levels
  assert(!donorGrid ||
    vtkUniformGrid::SafeDownCast(amrds->GetDataSetAsCartesianGrid(donorLevel, donorGridId)) ==
      donorGrid);
  return (donorCellIdx);
}

//------------------------------------------------------------------------------
bool vtkAMRResampleFilter::SearchGridAncestors(
  double q[3], vtkOverlappingAMR* amrds, unsigned int& level, unsigned int& gridId, int& cellId)
{
  assert("pre: AMR dataset is nullptr" && (amrds != nullptr));
  unsigned int *parents, plevel;
  for (; level > 0; --level)
  {
    ++this->NumberOfTimesLevelUp;
    // Get the parents of the grid

    unsigned int numParents;
    parents = amrds->GetParents(level, gridId, numParents);
    plevel = level - 1;
    // There should be at least 1 parent
    assert("Found non-level 0 grid with no parents" && (parents != nullptr) && (numParents > 0));
    if (numParents > 1)
    {
      vtkDebugMacro("Number of parents: " << numParents << " - Only processing 1 route");
    }
    gridId = parents[0];
    if (amrds->GetOverlappingAMRMetaData()->FindCell(q, plevel, gridId, cellId))
    {
      level = plevel;
      return true;
    }
  }
  // If we are here then we could not find an ancestor
  cellId = -1;
  return false;
}

//------------------------------------------------------------------------------
void vtkAMRResampleFilter::SearchGridDecendants(double q[3], vtkOverlappingAMR* amrds,
  unsigned int maxLevel, unsigned int& level, unsigned int& gridId, int& cellId)
{
  assert("pre: AMR dataset is nullptr" && (amrds != nullptr));
  unsigned int *children, clevel, n, i;
  for (; level < maxLevel - 1; ++level)
  {
    // Get the children of the grid
    children = amrds->GetChildren(level, gridId, n);
    clevel = level + 1;
    // If there are no children then we found the grid!
    if (children == nullptr)
    {
      return;
    }
    //    assert(n == children[0]);
    for (i = 0; i < n; ++i)
    {
      if (amrds->GetOverlappingAMRMetaData()->FindCell(q, clevel, children[i], cellId))
      {
        // We found a descendant so stop searching the
        // children and can instead search that grid's
        // children
        gridId = children[i];
        ++this->NumberOfTimesLevelDown;
        break;
      }
    }
    if (i >= n)
    {
      // We tested some children that we didn't need to if
      // we had visibility info
      this->NumberOfBlocksVisSkipped += n;
      // If we are here then no child contains the point
      // so don't search any further
      return;
    }
  }
}

//------------------------------------------------------------------------------
int vtkAMRResampleFilter::ProbeGridPointInAMRGraph(double q[3], unsigned int& donorLevel,
  unsigned int& donorGridId, vtkOverlappingAMR* amrds, unsigned int maxLevel, bool useCached)
{
  assert("pre: AMR dataset is nullptr" && amrds != nullptr);

  int donorCellIdx = -1;

  vtkUniformGrid* donorGrid = nullptr;
  // STEP 0: Check the previously cached donor-grid
  if (useCached)
  {
    if (!amrds->GetOverlappingAMRMetaData()->FindCell(q, donorLevel, donorGridId, donorCellIdx))
    {
      // Lets find the grid's ancestor that contains the point
      bool res = this->SearchGridAncestors(q, amrds, donorLevel, donorGridId, donorCellIdx);
      donorGrid = res
        ? vtkUniformGrid::SafeDownCast(amrds->GetDataSetAsCartesianGrid(donorLevel, donorGridId))
        : nullptr;
    }
    else
    {
      donorGrid =
        vtkUniformGrid::SafeDownCast(amrds->GetDataSetAsCartesianGrid(donorLevel, donorGridId));
      ++this->NumberOfTimesFoundOnDonorLevel;
    }
    // if the point is not contained in an ancestor then lets just assume its on level
    // 0 which is the default
  }

  // If there is no initial donor grid then search level 0
  if (donorGrid == nullptr)
  {
    bool res = SearchForDonorGridAtLevel(q, amrds, 0, donorGridId, donorCellIdx);
    // If we still can't find a grid then the point is not contained in the
    // AMR Data
    if (!res)
    {
      this->NumberOfFailedPoints++;
      donorLevel = 0;
      return -1;
    }
  }

  // Now search the descendants of the donor grid
  this->SearchGridDecendants(q, amrds, maxLevel, donorLevel, donorGridId, donorCellIdx);
  return (donorCellIdx);
}

//------------------------------------------------------------------------------
void vtkAMRResampleFilter::TransferToGridNodes(vtkUniformGrid* g, vtkOverlappingAMR* amrds)
{
  this->NumberOfBlocksTested = 0;
  this->NumberOfBlocksVisSkipped = 0;
  this->NumberOfTimesFoundOnDonorLevel = 0;
  this->NumberOfTimesLevelUp = 0;
  this->NumberOfTimesLevelDown = 0;
  this->NumberOfFailedPoints = 0;
  this->AverageLevel = 0.0;
  assert("pre: uniform grid is nullptr" && (g != nullptr));
  assert("pre: AMR data-structure is nullptr" && (amrds != nullptr));

  // STEP 0: Initialize the fields on the grid
  vtkUniformGrid* refGrid = this->GetReferenceGrid(amrds);

  vtkCellData* CD = refGrid->GetCellData();
  assert("pre: Donor CellData is nullptr!" && (CD != nullptr));

  vtkPointData* PD = g->GetPointData();
  assert("pre: Target PointData is nullptr!" && (PD != nullptr));

  // STEP 0: Initialize the fields on the grid
  this->InitializeFields(PD, g->GetNumberOfPoints(), CD);

  // STEP 1: If no arrays are selected, there is no need to interpolate
  // anything on the grid, just return
  if (PD->GetNumberOfArrays() == 0)
  {
    return;
  }

  // STEP 2: Fix the maximum level at which the search algorithm will operate
  unsigned int maxLevelToLoad = 0;
  if (this->LevelOfResolution < static_cast<int>(amrds->GetNumberOfLevels()) &&
    this->DemandDrivenMode == 1)
  {
    maxLevelToLoad = this->LevelOfResolution + 1;
  }
  else
  {
    maxLevelToLoad = amrds->GetNumberOfLevels();
  }

  // STEP 3: Loop through all the points and find the donors.
  unsigned int donorLevel = 0;
  unsigned int donorGridId = 0;
  double qPoint[3];
  vtkIdType pIdx;
  int donorCellIdx;
  bool useCached(false);
  // Do we have parent/child meta information (yes, we always do)
  if (this->AMRMetaData)
  {
    for (pIdx = 0; pIdx < g->GetNumberOfPoints(); ++pIdx)
    {
      g->GetPoint(pIdx, qPoint);
      donorCellIdx = this->ProbeGridPointInAMRGraph(
        qPoint, donorLevel, donorGridId, amrds, maxLevelToLoad, useCached);
      if (donorCellIdx != -1)
      {
        useCached = true;
        vtkUniformGrid* amrGrid =
          vtkUniformGrid::SafeDownCast(amrds->GetDataSetAsCartesianGrid(donorLevel, donorGridId));
        this->AverageLevel += donorLevel;
        CD = amrGrid->GetCellData();
        this->CopyData(PD, pIdx, CD, donorCellIdx);
      }
      else
      {
        useCached = false;
        // Point is outside the domain, blank it
        g->BlankPoint(pIdx);
      }
    } // END for all grid nodes
  }
  else
  {
    for (pIdx = 0; pIdx < g->GetNumberOfPoints(); ++pIdx)
    {
      g->GetPoint(pIdx, qPoint);

      donorCellIdx = this->ProbeGridPointInAMR(
        qPoint, donorLevel, donorGridId, amrds, maxLevelToLoad, useCached);

      if (donorCellIdx != -1)
      {
        useCached = true;
        this->AverageLevel += donorLevel;
        vtkUniformGrid* donorGrid =
          vtkUniformGrid::SafeDownCast(amrds->GetDataSetAsCartesianGrid(donorLevel, donorGridId));
        assert(donorGrid != nullptr);
        CD = donorGrid->GetCellData();
        this->CopyData(PD, pIdx, CD, donorCellIdx);
      }
      else
      {
        useCached = false;
        // Point is outside the domain, blank it
        g->BlankPoint(pIdx);
      }
    } // END for all grid nodes
  }
  std::cerr << "********* Resample Stats *************\n";
  double c = this->NumberOfSamples[0] * this->NumberOfSamples[1] * this->NumberOfSamples[2];
  double b = g->GetNumberOfPoints();
  std::cerr << "Number of Requested Points: " << c << " Number of Actual Points: " << b << "\n";
  std::cerr << " Percentage of Requested Points in Grid: " << 100.0 * b / c << "\n";
  std::cerr << "Total Number of Blocks Tested: " << this->NumberOfBlocksTested << "\n";
  std::cerr << " Number of Blocks that could be skipped by Visibility: "
            << this->NumberOfBlocksVisSkipped << "\n";
  double a = 100.0 * (double)(this->NumberOfBlocksVisSkipped) / (double)this->NumberOfBlocksTested;
  std::cerr << "Percentage of Blocks skipped via Visibility: " << a << "\n";
  a = (double)this->NumberOfBlocksTested / b;
  std::cerr << "Ave Number of Blocks Tested per Point: " << a << "\n";
  a = 100.0 * (double)this->NumberOfTimesFoundOnDonorLevel / b;
  std::cerr << "Percentage of Times we found point on Previous Level: " << a << "\n";
  a = 100.0 * (double)this->NumberOfTimesLevelUp / b;
  std::cerr << "Percentage of Times went to finer level: " << a << "\n";
  a = 100.0 * (double)this->NumberOfTimesLevelDown / b;
  std::cerr << "Percentage of Times went to coarser level: " << a << "\n";
  a = this->AverageLevel / b;
  std::cerr << "Average Level: " << a << "\n";
  std::cerr << "Number Of Failed Points: " << this->NumberOfFailedPoints << "\n";
}

//------------------------------------------------------------------------------
void vtkAMRResampleFilter::TransferSolution(vtkUniformGrid* g, vtkOverlappingAMR* amrds)
{
  assert("pre: uniform grid is nullptr" && (g != nullptr));
  assert("pre: AMR data-structure is nullptr" && (amrds != nullptr));

  if (this->TransferToNodes == 1)
  {
    this->TransferToGridNodes(g, amrds);
  }
  else
  {
    this->TransferToCellCenters(g, amrds);
  }
}

//------------------------------------------------------------------------------
void vtkAMRResampleFilter::ExtractRegion(
  vtkOverlappingAMR* amrds, vtkMultiBlockDataSet* mbds, vtkOverlappingAMR* vtkNotUsed(metadata))
{

  assert("pre: input AMR data-structure is nullptr" && (amrds != nullptr));
  assert("pre: resampled grid should not be nullptr" && (mbds != nullptr));

  //  std::cout << "NumBlocks: " << this->ROI->GetNumberOfBlocks() << std::endl;
  //  std::cout << "NumProcs: "  << this->Controller->GetNumberOfProcesses() << std::endl;
  //  std::cout.flush();

  assert("pre: NumProcs must be less than or equal to NumBlocks" &&
    (!this->Controller ||
      (static_cast<int>(this->ROI->GetNumberOfBlocks()) <=
        this->Controller->GetNumberOfProcesses())));

  mbds->SetNumberOfBlocks(this->ROI->GetNumberOfBlocks());
  for (unsigned int block = 0; block < this->ROI->GetNumberOfBlocks(); ++block)
  {
    if (this->CheckAbort())
    {
      break;
    }
    if (this->IsRegionMine(block))
    {
      vtkUniformGrid* grid = vtkUniformGrid::New();
      grid->ShallowCopy(this->ROI->GetBlock(block));
      this->TransferSolution(grid, amrds);
      mbds->SetBlock(block, grid);
      grid->Delete();
    }
    else
    {
      mbds->SetBlock(block, nullptr);
    }
  } // END for all blocks
}

//------------------------------------------------------------------------------
void vtkAMRResampleFilter::ComputeAMRBlocksToLoad(vtkOverlappingAMR* metadata)
{
  assert("pre: metadata is nullptr" && (metadata != nullptr));

  this->BlocksToLoad.clear();

  unsigned int maxLevelToLoad = 0;
  if (this->LevelOfResolution < static_cast<int>(metadata->GetNumberOfLevels()))
  {
    maxLevelToLoad = this->LevelOfResolution + 1;
  }
  else
  {
    maxLevelToLoad = metadata->GetNumberOfLevels();
  }

  unsigned int level = 0;
  for (; level < maxLevelToLoad; ++level)
  {
    unsigned int dataIdx = 0;
    for (; dataIdx < metadata->GetNumberOfBlocks(level); ++dataIdx)
    {
      double grd[6];
      metadata->GetBounds(level, dataIdx, grd);
      if (this->IsBlockWithinBounds(grd))
      {
        this->BlocksToLoad.push_back(metadata->GetAbsoluteBlockIndex(level, dataIdx));
      } // END check if the block is within the bounds of the ROI
    }   // END for all data
  }     // END for all levels

  std::sort(this->BlocksToLoad.begin(), this->BlocksToLoad.end());
  std::cerr << "Number Levels Loaded = " << maxLevelToLoad
            << " Number of Blocks = " << BlocksToLoad.size() << "\n";
}

//------------------------------------------------------------------------------
void vtkAMRResampleFilter::GetDomainParameters(vtkOverlappingAMR* amr, double domainMin[3],
  double domainMax[3], double h[3], int dims[3], double& rf)
{
  assert("pre: AMR dataset is nullptr!" && (amr != nullptr));

  rf = amr->GetRefinementRatio(1);
  amr->GetAMRBox(0, 0).GetNumberOfNodes(dims);
  amr->GetMin(domainMin);
  amr->GetMax(domainMax);

  amr->GetSpacing(0, h);
}

//------------------------------------------------------------------------------
void vtkAMRResampleFilter::SnapBounds(const double* vtkNotUsed(h0[3]), const double domainMin[3],
  const double domainMax[3], const int* vtkNotUsed(dims[3]), bool outside[6])
{
  int i, j;
  for (i = 0, j = 0; i < 3; ++i)
  {
    // Snap the parts of the bounds that lie outside of the AMR data
    if (this->Min[i] < domainMin[i])
    {
      outside[j++] = true;
      this->GridMin[i] = domainMin[i];
    }
    else
    {
      outside[j++] = false;
      this->GridMin[i] = this->Min[i];
    }

    if (this->Max[i] > domainMax[i])
    {
      outside[j++] = true;
      this->GridMax[i] = domainMax[i];
    }
    else
    {
      outside[j++] = false;
      this->GridMax[i] = this->Max[i];
    }
  }
}

//------------------------------------------------------------------------------
void vtkAMRResampleFilter::ComputeLevelOfResolution(
  const int N[3], const double h0[3], const double L[3], double rf)
{
  this->LevelOfResolution = 0;
  for (int i = 0; i < 3; ++i)
  {
    double c1 = ((N[i] * h0[i]) / L[i]);
    int currentLevel = vtkMath::Floor(0.5 + (log(c1) / log(rf)));
    this->LevelOfResolution = std::max(currentLevel, this->LevelOfResolution);
  } // END for all i
  std::cerr << "Requested Max Level = " << this->LevelOfResolution << "\n";
}

//------------------------------------------------------------------------------
bool vtkAMRResampleFilter::RegionIntersectsWithAMR(
  double domainMin[3], double domainMax[3], double regionMin[3], double regionMax[3])
{
  vtkBoundingBox domain;
  domain.SetMinPoint(domainMin);
  domain.SetMaxPoint(domainMax);

  vtkBoundingBox region;
  region.SetMinPoint(regionMin);
  region.SetMaxPoint(regionMax);

  return domain.Intersects(region);
}

//------------------------------------------------------------------------------
void vtkAMRResampleFilter::AdjustNumberOfSamplesInRegion(
  const double Rh[3], const bool outside[6], int N[3])
{
  int startIndex;
  int endIndex;
  double dx = 0.0;
  for (int i = 0; i < 3; ++i)
  {
    N[i] = this->NumberOfSamples[i];

    // Get ijk of the snapped bounding wrt the requested virtual grid.
    if (outside[i * 2] || outside[i * 2 + 1])
    {
      dx = this->GridMin[i] - this->Min[i];
      if (dx > 0.0)
      {
        startIndex = static_cast<int>(dx / Rh[i] + 1);
      }
      else
      {
        startIndex = 0;
      }

      dx = this->GridMax[i] - this->Min[i];
      endIndex = static_cast<int>(dx / Rh[i] + 1);

      endIndex = std::min(endIndex, N[i]);
      int newN = endIndex - startIndex + 1;
      if (newN <= N[i])
      {
        N[i] = newN;
      }
      else
      {
        assert("ERROR: code should not reach here!" && false);
      }
    }
  }
  std::cerr << "Request Grid Dim : " << this->NumberOfSamples[0] << ", " << this->NumberOfSamples[1]
            << ", " << this->NumberOfSamples[2] << "\n";
  std::cerr << "Computed Grid Dim: " << N[0] << ", " << N[1] << ", " << N[2] << "\n";
  if (this->UseBiasVector)
  {
    double a[3];
    a[0] = fabs(this->BiasVector[0]);
    a[1] = fabs(this->BiasVector[1]);
    a[2] = fabs(this->BiasVector[2]);

    // Find the max component
    // NOLINTNEXTLINE(readability-avoid-nested-conditional-operator)
    int bdir = (a[0] > a[1]) ? ((a[0] > a[2]) ? 0 : 2) : ((a[1] > a[2]) ? 1 : 2);

    if (bdir == 0)
    {
      N[0] = std::min(N[0], std::max(N[1], N[2]));
    }
    else if (bdir == 1)
    {
      N[1] = std::min(N[1], std::max(N[0], N[2]));
    }
    else
    {
      N[2] = std::min(N[2], std::max(N[0], N[1]));
    }
    std::cerr << "Adjusted Grid Dim: " << N[0] << ", " << N[1] << ", " << N[2] << "\n";
  }
}

//------------------------------------------------------------------------------
void vtkAMRResampleFilter::ComputeAndAdjustRegionParameters(vtkOverlappingAMR* amrds, double h[3])
{
  assert("pre: AMR dataset is nullptr" && (amrds != nullptr));

  // STEP 0: Get domain parameters from root level metadata
  int dims[3];
  double h0[3];
  double domainMin[3];
  double domainMax[3];
  double rf;
  this->GetDomainParameters(amrds, domainMin, domainMax, h0, dims, rf);

  // STEP 1: Check to see if the requested region intersects the AMR domain
  if (!this->RegionIntersectsWithAMR(domainMin, domainMax, this->Min, this->Max))
  {
    h[0] = h[1] = h[2] = 0.0;
    return;
  }

  // STEP 3: Get requested region parameters
  double L0[3]; // initial length of each box side
  double Rh[3]; // initial spacing based on the number of samples requested.
  for (int i = 0; i < 3; ++i)
  {
    L0[i] = this->Max[i] - this->Min[i];
    Rh[i] = L0[i] / (this->NumberOfSamples[i] - 1);
  }

  // DEBUG
  //  this->WriteUniformGrid(this->Min,this->NumberOfSamples,Rh,"RequestedGrid" );
  // END DEBUG
  // DEBUG
  //  this->WriteUniformGrid(domainMin,dims,h0,"RootGrid");
  // END DEBUG

  // STEP 4: Snap region to domain bounds
  bool outside[6];
  // Determine the Min/Max of the computed grid
  this->SnapBounds(h0, domainMin, domainMax, dims, outside);

  // STEP 5: Compute grid parameters on the snapped region
  double L[3];
  for (int i = 0; i < 3; ++i)
  {
    L[i] = this->GridMax[i] - this->GridMin[i];
    h[i] = L[i] / (this->NumberOfSamples[i] - 1);
  }

  // DEBUG
  //  this->WriteUniformGrid(min,this->NumberOfSamples,h,"SnappedGrid");
  // END DEBUG

  // STEP 6: Adjust N according to how much of the requested region is cropped
  int N[3];
  this->AdjustNumberOfSamplesInRegion(Rh, outside, N);

  // STEP 7: Adjust region parameters
  for (int i = 0; i < 3; ++i)
  {
    this->GridNumberOfSamples[i] = (N[i] > 1) ? N[i] : 2;
    h[i] = L[i] / (this->GridNumberOfSamples[i] - 1);
  }

  this->ComputeLevelOfResolution(this->GridNumberOfSamples, h0, L, rf);
}

//------------------------------------------------------------------------------
void vtkAMRResampleFilter::GetRegion(double h[3])
{
  assert("pre: Region of interest is nullptr!" && (this->ROI != nullptr));

  unsigned int block = 0;
  for (; block < this->ROI->GetNumberOfBlocks(); ++block)
  {
    this->ROI->RemoveBlock(block);
  }

  if (h[0] == 0.0 && h[1] == 0.0 && h[2] == 0.0)
  {
    return;
  }

  vtkUniformGrid* grd = vtkUniformGrid::New();
  grd->SetOrigin(this->GridMin);
  grd->SetSpacing(h);
  grd->SetDimensions(this->GridNumberOfSamples);
  if (grd->GetNumberOfPoints() == 0)
  {
    vtkErrorMacro("Empty Grid!");
    return;
  }

  vtkUniformGridPartitioner* gridPartitioner = vtkUniformGridPartitioner::New();
  gridPartitioner->SetInputData(grd);
  grd->Delete();

  gridPartitioner->SetNumberOfPartitions(this->NumberOfPartitions);
  gridPartitioner->Update();

  this->ROI->DeepCopy(gridPartitioner->GetOutput());

  gridPartitioner->Delete();
}

//------------------------------------------------------------------------------
bool vtkAMRResampleFilter::GridsIntersect(double* g1, double* g2)
{
  assert("pre: g1 is nullptr" && (g1 != nullptr));
  assert("pre: g2 is nullptr" && (g2 != nullptr));

  vtkBoundingBox b1;
  b1.SetBounds(g1);

  vtkBoundingBox b2;
  b2.SetBounds(g2);

  return b1.IntersectBox(b2);
}

//------------------------------------------------------------------------------
bool vtkAMRResampleFilter::IsBlockWithinBounds(double* grd)
{
  assert("pre: Input AMR grid is nullptr" && (grd != nullptr));

  for (unsigned int block = 0; block < this->ROI->GetNumberOfBlocks(); ++block)
  {
    if (this->IsRegionMine(block))
    {
      vtkUniformGrid* blk = vtkUniformGrid::SafeDownCast(this->ROI->GetBlock(block));
      assert("pre: block is nullptr" && (blk != nullptr));

      if (this->GridsIntersect(grd, blk->GetBounds()))
      {
        return true;
      }
    } // END if region is mine
  }   // END for all blocks

  return false;
}

//------------------------------------------------------------------------------
int vtkAMRResampleFilter::GetRegionProcessId(int regionIdx)
{
  if (!this->IsParallel())
  {
    return 0;
  }

  int N = this->Controller->GetNumberOfProcesses();
  return (regionIdx % N);
}

//------------------------------------------------------------------------------
bool vtkAMRResampleFilter::IsRegionMine(int regionIdx)
{
  if (!this->IsParallel())
  {
    return true;
  }

  int myRank = this->Controller->GetLocalProcessId();
  return myRank == this->GetRegionProcessId(regionIdx);
}

//------------------------------------------------------------------------------
bool vtkAMRResampleFilter::IsParallel()
{
  if (this->Controller == nullptr)
  {
    return false;
  }

  return this->Controller->GetNumberOfProcesses() > 1;
}

//------------------------------------------------------------------------------
vtkUniformGrid* vtkAMRResampleFilter::GetReferenceGrid(vtkOverlappingAMR* amrds)
{
  assert("pre:AMR dataset is nullptr" && (amrds != nullptr));

  unsigned int numLevels = amrds->GetNumberOfLevels();
  for (unsigned int l = 0; l < numLevels; ++l)
  {
    unsigned int numDatasets = amrds->GetNumberOfBlocks(l);
    for (unsigned int dataIdx = 0; dataIdx < numDatasets; ++dataIdx)
    {
      vtkUniformGrid* refGrid =
        vtkUniformGrid::SafeDownCast(amrds->GetDataSetAsCartesianGrid(l, dataIdx));
      if (refGrid != nullptr)
      {
        return (refGrid);
      }
    } // END for all datasets
  }   // END for all number of levels

  // This process has no grids
  return nullptr;
}
VTK_ABI_NAMESPACE_END