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

  Program:   ParaView
  Module:    vtkAMRConnectivity.cxx

  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.

  Copyright 2013 Sandia Corporation.
  Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
  the U.S. Government retains certain rights in this software.

=========================================================================*/
#include "vtkAMRConnectivity.h"

#include "vtkObject.h"
#include "vtkObjectFactory.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"

#include "vtkAMRDualGridHelper.h"
#include "vtkCellData.h"
#include "vtkCell.h"
#include "vtkCompositeDataIterator.h"
#include "vtkDataArray.h"
#include "vtkDataSet.h"
#include "vtkDataObject.h"
#include "vtkDoubleArray.h"
#include "vtkIdList.h"
#include "vtkIdTypeArray.h"
#include "vtkIntArray.h"
#include "vtkMultiProcessController.h"
#include "vtkNonOverlappingAMR.h"
#include "vtkSmartPointer.h"
#include "vtkTimerLog.h"
#include "vtkUniformGrid.h"
#include "vtkUnsignedCharArray.h"
#include "vtksys/SystemTools.hxx"

#include "vtkPVConfig.h"
#ifdef PARAVIEW_USE_MPI
#include "vtkMPIController.h"
#endif

#include <map>
#include <list>

vtkStandardNewMacro (vtkAMRConnectivity);

class vtkAMRConnectivityEquivalence
{
public:
  vtkAMRConnectivityEquivalence () 
    {
    set_to_min_id = vtkSmartPointer<vtkIntArray>::New ();
    set_to_min_id->SetNumberOfComponents (1);
    set_to_min_id->SetNumberOfTuples (0);
    }

  ~vtkAMRConnectivityEquivalence ()
    {
    }

  int AddEquivalence (int id1, int id2) 
    {
    int min_id;
    if (id1 < id2)
      {
      min_id = id1;
      }
    else
      {
      min_id = id2;
      }

    std::map<int, int>::iterator iter;

    iter = id_to_set.find (id1);
    int set1 = (iter == id_to_set.end () ? -1 : iter->second);
    iter = id_to_set.find (id2);
    int set2 = (iter == id_to_set.end () ? -1 : iter->second);

    if (set1 >= 0 && set_to_min_id->GetValue (set1) < 0)  
      {
      set1 = -1;
      }
    if (set2 >= 0 && set_to_min_id->GetValue (set2) < 0)
      {
      set2 = -1;
      }

    if (set1 >= 0 && set2 >= 0 && set1 == set2) 
      {
      return 0;
      }
    else if (set1 >= 0 && set2 >= 0)
      {
      // merge sets
      int min_set, max_set;
      if (set1 < set2) 
        {
        min_set = set1;
        max_set = set2;
        }
      else
        {
        min_set = set2;
        max_set = set1;
        }
      for (iter = id_to_set.begin (); iter != id_to_set.end (); iter ++) 
        {
        if (iter->second == max_set)
          {
          id_to_set[iter->first] = min_set;
          }
        }
      int max_set_min = set_to_min_id->GetValue (max_set);
      int min_set_min = set_to_min_id->GetValue (min_set);
      // pick the smallest of the two mins to represent the set.
      if (max_set_min < min_set_min) 
        {
        set_to_min_id->SetValue (min_set, max_set_min);
        }
      set_to_min_id->SetValue (max_set, -1);
      }
    else if (set1 >= 0)
      {
      id_to_set[id2] = set1;
      if (id2 < set_to_min_id->GetValue (set1)) 
        {
        set_to_min_id->SetValue (set1, id2);
        }
      }
    else if (set2 >= 0)
      {
      id_to_set[id1] = set2;
      if (id1 < set_to_min_id->GetValue (set2)) 
        {
        set_to_min_id->SetValue (set2, id1);
        }
      }
    else
      {
      // find an empty set.
      int first_empty = -1;
      for (int i = 0; i < set_to_min_id->GetNumberOfTuples (); i ++)
        {
        if (set_to_min_id->GetValue (i) < 0)
          {
          first_empty = i;
          break;
          }
        }
      if (first_empty < 0)
        {
        first_empty = set_to_min_id->InsertNextValue (-1);
        }

      id_to_set[id1] = first_empty;
      id_to_set[id2] = first_empty;
      set_to_min_id->SetValue (first_empty, min_id);
      // return zero here because its not values we've previously cared about.
      }
    return 1;
    }

  int GetMinimumSetId (int id)
    {
    std::map<int, int>::iterator iter = id_to_set.find (id);
    if (iter == id_to_set.end ())
      {
      // vtkErrorWithObjectMacro (id_to_set, << "ID out of range " << id << " (expected 0 <= x < " << id_to_set->GetNumberOfTuples () << ")");
      return -1;
      }
    int set = id_to_set[id];
    return (set >= 0 ? set_to_min_id->GetValue (set) : -1);
    }
    
private:
  std::map<int,int> id_to_set;
  vtkSmartPointer<vtkIntArray> set_to_min_id;
};


#ifdef PARAVIEW_USE_MPI

static const int BOUNDARY_TAG = 857089;
static const int EQUIV_SIZE_TAG = 748957;
static const int EQUIV_TAG = 357345;

//-----------------------------------------------------------------------------
// Simple containers for managing asynchronous communication.
struct vtkAMRConnectivityCommRequest
{
  vtkMPICommunicator::Request Request;
  vtkSmartPointer<vtkIntArray> Buffer;
  int SendProcess;
  int ReceiveProcess;
};

// This class is a STL list of vtkAMRConnectivityCommRequest structs with some
// helper methods added.
class vtkAMRConnectivityCommRequestList
  : public std::list<vtkAMRConnectivityCommRequest>
{
public:
  // Description:
  // Waits for all of the communication to complete.
  void WaitAll()
  {
    for (iterator i = this->begin(); i != this->end(); i++) 
      {
      i->Request.Wait();
      // this->erase (i);
      }
  }
  // Description:
  // Waits for one of the communications to complete, removes it from the list,
  // and returns it.
  value_type WaitAny()
  {
    while (!this->empty())
      {
      for (iterator i = this->begin(); i != this->end(); i++)
        {
        if (i->Request.Test())
          {
          value_type retval = *i;
          this->erase(i);
          return retval;
          }
        }
      vtksys::SystemTools::Delay(1);
      }
    vtkGenericWarningMacro(<< "Nothing to wait for.");
    return value_type();
  }
};
#endif /* PARAVIEW_USE_MPI */


//-----------------------------------------------------------------------------
vtkAMRConnectivity::vtkAMRConnectivity ()
{
  this->VolumeFractionSurfaceValue = 0.5;
  this->Helper = 0;
  this->Equivalence = 0;
  this->ResolveBlocks = 1;
  this->PropagateGhosts = 0;
}

vtkAMRConnectivity::~vtkAMRConnectivity ()
{
}

void vtkAMRConnectivity::PrintSelf (ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf (os, indent);
  os << indent << "VolumeFractionSurfaceValue: " << this->VolumeFractionSurfaceValue << endl;
}

void vtkAMRConnectivity::AddInputVolumeArrayToProcess (const char *name) 
{
  this->VolumeArrays.push_back (name);
  this->Modified();
}

void vtkAMRConnectivity::ClearInputVolumeArrayToProcess ()
{
  this->VolumeArrays.clear ();
  this->Modified();
}

int vtkAMRConnectivity::FillInputPortInformation (int port, vtkInformation *info)
{
  switch (port) 
    {
    case 0:
      info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkNonOverlappingAMR");
      break;
    default:
      return 0;
    }

  return 1;
}

//----------------------------------------------------------------------------
int vtkAMRConnectivity::FillOutputPortInformation(int port, vtkInformation *info)
{
  switch (port)
    {
    case 0:
      info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkNonOverlappingAMR");
      break;
    default:
      return 0;
    }

  return 1;
}

//----------------------------------------------------------------------------
int vtkAMRConnectivity::RequestData (vtkInformation* vtkNotUsed(request),
                                  vtkInformationVector** inputVector,
                                  vtkInformationVector* outputVector)
{
  vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
  vtkNonOverlappingAMR* amrInput=vtkNonOverlappingAMR::SafeDownCast(
    inInfo->Get(vtkDataObject::DATA_OBJECT()));

  vtkInformation *outInfo;
  outInfo = outputVector->GetInformationObject(0);
  vtkNonOverlappingAMR* amrOutput = vtkNonOverlappingAMR::SafeDownCast(
    outInfo->Get(vtkDataObject::DATA_OBJECT()));

  amrOutput->ShallowCopy (amrInput);

  this->Helper = vtkAMRDualGridHelper::New ();
  vtkMultiProcessController* controller = vtkMultiProcessController::GetGlobalController ();
  this->Helper->SetController(controller);
  this->Helper->Initialize(amrInput);

  unsigned int noOfArrays = static_cast<unsigned int>(this->VolumeArrays.size());
  for(unsigned int i = 0; i < noOfArrays; i++)
    {
    if (this->DoRequestData (amrOutput, this->VolumeArrays[i].c_str()) == 0) 
      {
      return 0;
      }
    }
  this->Helper->Delete ();

  return 1;
}


//----------------------------------------------------------------------------
int vtkAMRConnectivity::DoRequestData (vtkNonOverlappingAMR* volume,
                                      const char* volumeName)
{
  vtkMultiProcessController* controller = vtkMultiProcessController::GetGlobalController ();
#ifdef PARAVIEW_USE_MPI
  vtkMPIController *mpiController = vtkMPIController::SafeDownCast(controller);
#endif
  int myProc = controller->GetLocalProcessId ();
  int numProcs = controller->GetNumberOfProcesses ();

  this->RegionName = std::string ("RegionId-");
  this->RegionName += volumeName;

  // initialize with a global unique region id
  this->NextRegionId = myProc+1;

  vtkTimerLog::MarkStartEvent ("Initial fragment seeding");

  // Find the block local fragments
  vtkCompositeDataIterator* iter = volume->NewIterator ();
  for (iter->InitTraversal (); !iter->IsDoneWithTraversal (); iter->GoToNextItem ())
    {
    // Go through each block and create an array RegionId set all zero.
    vtkUniformGrid* grid = vtkUniformGrid::SafeDownCast (iter->GetCurrentDataObject ());
    if (!grid)
      {
      vtkErrorMacro ("NonOverlappingAMR not made up of UniformGrids");
      return 0;
      }
    vtkSmartPointer<vtkIdTypeArray> regionId = vtkSmartPointer<vtkIdTypeArray>::New ();
    regionId->SetName (this->RegionName.c_str());
    regionId->SetNumberOfComponents (1);
    regionId->SetNumberOfTuples (grid->GetNumberOfCells ());
    for (int i = 0; i < grid->GetNumberOfCells (); i ++) 
      {
      regionId->SetTuple1 (i, 0);
      }
    grid->GetCellData ()->AddArray (regionId);

    vtkDataArray* volArray = grid->GetCellData ()->GetArray (volumeName);
    if (!volArray)
      {
      vtkErrorMacro (<< "There is no " << volumeName << " in cell field");
      return 0;
      }

    vtkUnsignedCharArray* ghostArray = grid->GetCellGhostArray();
    if (!ghostArray) 
      {
      vtkErrorMacro ("No ghost array attached to the CTH volume data");
      return 0;
      }

    // within each block find all fragments
    int extents[6];
    grid->GetExtent (extents);

    for (int i = extents[0]; i < extents[1]; i ++) 
      {
      for (int j = extents[2]; j < extents[3]; j ++) 
        {
        for (int k = extents[4]; k < extents[5]; k ++) 
          {
          int ijk[3] = { i, j, k };
          vtkIdType cellId = grid->ComputeCellId (ijk);
          if (regionId->GetTuple1 (cellId) == 0 &&
              volArray->GetTuple1 (cellId) > this->VolumeFractionSurfaceValue &&
              (ghostArray->GetValue (cellId) & 
               vtkDataSetAttributes::DUPLICATECELL) == 0)
            {
            // wave propagation sets the region id as it propagates
            this->WavePropagation (cellId, grid, regionId, volArray, ghostArray);
            // increment this way so the region id remains globally unique
            this->NextRegionId += numProcs;
            }
          }
        }
      }
    }

  vtkTimerLog::MarkEndEvent ("Initial fragment seeding");

  if (this->ResolveBlocks)
    {
    vtkTimerLog::MarkStartEvent ("Computing boundary regions");

    // Determine boundaries at the block that need to be sent to neighbors
    this->BoundaryArrays.resize (numProcs);
    this->ReceiveList.resize (numProcs);
    this->ValidNeighbor.resize (numProcs);
    this->NeighborList.resize (this->Helper->GetNumberOfLevels ());
    for (int level = 0; level < this->Helper->GetNumberOfLevels (); level ++)
      {
      int maxId = 0;
      for (int blockId = 0; blockId < this->Helper->GetNumberOfBlocksInLevel (level); blockId ++) 
        {
        vtkAMRDualGridHelperBlock* block = this->Helper->GetBlock (level, blockId);
        if (block->BlockId > maxId)
          {
          maxId = block->BlockId;
          }
        }
      this->NeighborList[level].resize (maxId + 1);
      for (int blockId = 0; blockId < maxId; blockId ++) 
        {
        this->NeighborList[level][blockId].resize (0);
        }
      }

    for (int level = 0; level < this->Helper->GetNumberOfLevels (); level ++)
      {
      for (int blockId = 0; blockId < this->Helper->GetNumberOfBlocksInLevel (level); blockId ++) 
        {
        vtkAMRDualGridHelperBlock* block = this->Helper->GetBlock (level, blockId);
        for (int dir = 0; dir < 3; dir ++)
          {
          vtkAMRDualGridHelperBlock* neighbor = this->GetBlockNeighbor (block, dir);
          if (neighbor != 0)
            {
            this->ProcessBoundaryAtBlock (volume, block, neighbor, dir);
            if (neighbor && neighbor->Level > block->Level) 
              {
              int n_level = neighbor->Level;
              int index[3] = { neighbor->GridIndex[0], neighbor->GridIndex[1], neighbor->GridIndex[2] };
  
              index[(dir+1) % 3] ++;
              neighbor = this->Helper->GetBlock (n_level, index[0], index[1], index[2]);
              if (neighbor != 0)
                { 
                this->ProcessBoundaryAtBlock (volume, block, neighbor, dir);
                }
  
              index[(dir+2) % 3] ++;
              neighbor = this->Helper->GetBlock (n_level, index[0], index[1], index[2]);
              if (neighbor != 0)
                { 
                this->ProcessBoundaryAtBlock (volume, block, neighbor, dir);
                }
  
              index[(dir+1) % 3] --;
              neighbor = this->Helper->GetBlock (n_level, index[0], index[1], index[2]);
              if (neighbor != 0)
                { 
                this->ProcessBoundaryAtBlock (volume, block, neighbor, dir);
                }
              }
            }
          }
        }
      }

#ifdef PARAVIEW_USE_MPI
    // Exchange all boundaries between processes where block and neighbor are different procs
    if (numProcs > 1  && !this->ExchangeBoundaries (mpiController))
      {
      return 0;
      }
#endif 
    // Process all boundaries at the neighbors to find the equivalence pairs at the boundaries
    this->Equivalence = new vtkAMRConnectivityEquivalence;

    // Initialize equivalence with all regions independent
    for (size_t i = 0; i < this->BoundaryArrays.size (); i ++) 
      {
      for (size_t j = 0; j < this->BoundaryArrays[i].size (); j ++)
        {
        if (i == static_cast<unsigned int> (myProc))
          {
          this->ProcessBoundaryAtNeighbor (volume, this->BoundaryArrays[myProc][j]);
          }
        }
        this->BoundaryArrays[i].clear ();
      }
    vtkTimerLog::MarkEndEvent ("Computing boundary regions");
  
    vtkTimerLog::MarkStartEvent ("Transferring equivalence");
    this->EquivPairs.resize (numProcs);
    while (true)
      {
      int sets_changed = 0;
      // Relabel all fragment IDs with the equivalence set number 
      // (set numbers start with 1 and 0 is considered "no set" or "no fragment")
      for (int level = 0; level < this->Helper->GetNumberOfLevels (); level ++)
        {
        for (int blockId = 0; blockId < this->Helper->GetNumberOfBlocksInLevel (level); blockId ++) 
          {
          vtkAMRDualGridHelperBlock* block = this->Helper->GetBlock (level, blockId);
          if (block->ProcessId != myProc) 
            {
            continue;
            }
          vtkUniformGrid* grid = volume->GetDataSet (block->Level, block->BlockId);
          vtkIdTypeArray* regionIdArray = vtkIdTypeArray::SafeDownCast (
                                            grid->GetCellData ()->GetArray (this->RegionName.c_str()));
          if (regionIdArray == 0)
            {
            vtkErrorMacro ("block Image doesn't not contain the regionId just added");
            return 0;
            }
          for (int i = 0; i < regionIdArray->GetNumberOfTuples (); i ++) 
            {
            vtkIdType regionId = regionIdArray->GetTuple1 (i);
            if (regionId > 0) 
              {
              int setId = this->Equivalence->GetMinimumSetId (regionId);
              if (setId > 0 && setId != regionId) 
                {
                regionIdArray->SetTuple1 (i, setId);
                sets_changed = 1;
                for (size_t p = 0; p < this->NeighborList[block->Level][block->BlockId].size (); p ++)
                  {
                  int n = this->NeighborList[block->Level][block->BlockId][p];
                  if (this->EquivPairs[n] == 0)
                    {
                    this->EquivPairs[n] = vtkSmartPointer<vtkIntArray>::New ();
                    this->EquivPairs[n]->SetNumberOfComponents (1);
                    this->EquivPairs[n]->SetNumberOfTuples (0);
                    }
                  int contained = 0;
                  for (int e = 0; e < this->EquivPairs[n]->GetNumberOfTuples (); e += 2)
                    {
                    int v1 = this->EquivPairs[n]->GetValue (e);
                    int v2 = this->EquivPairs[n]->GetValue (e+1);
                    if ((v1 == regionId && v2 == setId) ||
                        (v2 == regionId && v1 == setId))
                      {
                      contained = 1;
                      }
                    }
                  if (contained == 0) 
                    {
                    this->EquivPairs[n]->InsertNextValue (regionId);
                    this->EquivPairs[n]->InsertNextValue (setId);
                    }
                  }
                }
              }
            else
              {
              regionIdArray->SetTuple1 (i, 0);
              }
            }
          }
        }
#ifdef PARAVIEW_USE_MPI
      if (numProcs > 1)
        {
        int out;
        controller->AllReduce (&sets_changed, &out, 1, vtkCommunicator::MAX_OP);
        sets_changed = out;
        }
#endif
      if (sets_changed == 0)
        {
        break;
        }
#ifdef PARAVIEW_USE_MPI
      if (numProcs > 1 && !this->ExchangeEquivPairs (mpiController)) 
        {
        return 0;
        }
#endif
      // clear out the pairs after sending them.
      for (int i = 0; i < numProcs; i ++)
        {
        if (this->EquivPairs[i] != 0)
          {
          this->EquivPairs[i]->SetNumberOfTuples (0);
          }
        }
      }

    // clean up EquivPairs for the last time this update
    this->EquivPairs.clear ();
    ValidNeighbor.clear ();
    NeighborList.clear ();

    vtkTimerLog::MarkEndEvent ("Transferring equivalence");

    delete this->Equivalence;
    this->Equivalence = 0;
    }

  if (PropagateGhosts) 
    {
    vtkTimerLog::MarkStartEvent ("Propagating ghosts");
    // Propagate the region IDs out to the ghosts
    vtkSmartPointer<vtkIdList> ptIds = vtkIdList::New ();
    vtkSmartPointer<vtkIdList> cellIds = vtkIdList::New ();
    for (iter->InitTraversal (); !iter->IsDoneWithTraversal (); iter->GoToNextItem ())
      {
      vtkUniformGrid* grid = vtkUniformGrid::SafeDownCast (iter->GetCurrentDataObject ());
      vtkIdTypeArray* regionIdArray = vtkIdTypeArray::SafeDownCast (
                                        grid->GetCellData ()->GetArray (this->RegionName.c_str()));
      vtkDataArray* volArray = grid->GetCellData ()->GetArray (volumeName);
      vtkUnsignedCharArray* ghostArray = grid->GetCellGhostArray();
      for (int cellId = 0; cellId < regionIdArray->GetNumberOfTuples (); cellId ++) 
        {
        if ((ghostArray->GetValue (cellId) & 
             vtkDataSetAttributes::DUPLICATECELL) == 0 ||
            volArray->GetTuple1 (cellId) < this->VolumeFractionSurfaceValue)
          {
          continue;
          }
        grid->GetCellPoints (cellId, ptIds);
        bool isSet = false;
        for (int i = 0; i < ptIds->GetNumberOfIds (); i ++)
          {
          grid->GetPointCells (ptIds->GetId (i), cellIds);
          for (int j = 0; j < cellIds->GetNumberOfIds (); j ++)
            {
            vtkIdType neighbor = cellIds->GetId (j);
            if (neighbor != cellId &&
                volArray->GetTuple1 (neighbor) > 
                  this->VolumeFractionSurfaceValue && 
                (ghostArray->GetValue (neighbor) & 
                 vtkDataSetAttributes::DUPLICATECELL) == 0)
              {
              regionIdArray->SetTuple1 (cellId, regionIdArray->GetTuple1 (neighbor));
              isSet = true;
              break;
              }
            }
          if (isSet) 
            {
            break;
            }
          }
        }
      }

    vtkTimerLog::MarkEndEvent ("Propagating ghosts");
    }


  return 1;
}

//----------------------------------------------------------------------------
int vtkAMRConnectivity::WavePropagation (vtkIdType cellIdStart, 
                                         vtkUniformGrid* grid, 
                                         vtkIdTypeArray* regionId,
                                         vtkDataArray* volArray,
                                         vtkUnsignedCharArray* ghostArray)
{
  vtkSmartPointer<vtkIdList> todoList = vtkIdList::New ();
  todoList->SetNumberOfIds (0);
  todoList->InsertNextId (cellIdStart);

  vtkSmartPointer<vtkIdList> ptIds = vtkIdList::New ();
  vtkSmartPointer<vtkIdList> cellIds = vtkIdList::New ();

  while (todoList->GetNumberOfIds () > 0)
    {
    vtkIdType cellId = todoList->GetId (0);
    todoList->DeleteId (cellId);
    regionId->SetTuple1 (cellId, this->NextRegionId);

    grid->GetCellPoints (cellId, ptIds);
    for (int i = 0; i < ptIds->GetNumberOfIds (); i ++)
      {
      cellIds->SetNumberOfIds (0);
      grid->GetPointCells (ptIds->GetId (i), cellIds);
      for (int j = 0; j < cellIds->GetNumberOfIds (); j ++)
        {
        vtkIdType neighbor = cellIds->GetId (j);
        if (neighbor == cellId || todoList->IsId (neighbor) >= 0) { continue; }
        if (regionId->GetTuple1 (neighbor) == 0 &&
            volArray->GetTuple1 (neighbor) > this->VolumeFractionSurfaceValue &&
            (ghostArray->GetValue (neighbor) &
             vtkDataSetAttributes::DUPLICATECELL) == 0)
          {
          todoList->InsertNextId (neighbor);
          }
        }
      }
    }
  return 1;
}

//----------------------------------------------------------------------------
vtkAMRDualGridHelperBlock* vtkAMRConnectivity::GetBlockNeighbor (
                       vtkAMRDualGridHelperBlock* block, 
                       int dir)
{
  int gridIndex[3] = { 
          block->GridIndex[0], 
          block->GridIndex[1],
          block->GridIndex[2] };
  gridIndex[dir] ++;

  vtkAMRDualGridHelperBlock* neighbor = 0;
  neighbor = this->Helper->GetBlock (block->Level, gridIndex[0], gridIndex[1], gridIndex[2]);
  if (neighbor != 0)
    {
    return neighbor;
    }
  neighbor = this->Helper->GetBlock (block->Level-1, gridIndex[0] >> 1, gridIndex[1] >> 1, gridIndex[2] >> 1);
  if (neighbor != 0)
    {
    return neighbor;
    }
  // It's only the 0,0,0 corner block of the set.  
  // We'll need 3 more +1 blocks depending on direction
  // e.g., 0, 0, 1  and 0, 1, 0 and 0, 1, 1  If it's in the x direction.
  neighbor = this->Helper->GetBlock (block->Level+1, gridIndex[0] << 1, gridIndex[1] << 1, gridIndex[2] << 1);
  if (neighbor != 0)
    {
    return neighbor;
    }

  return 0;
}

//----------------------------------------------------------------------------
void vtkAMRConnectivity::ProcessBoundaryAtBlock (
                       vtkNonOverlappingAMR* volume,
                       vtkAMRDualGridHelperBlock* block, 
                       vtkAMRDualGridHelperBlock* neighbor, 
                       int dir)
{
  if (block == 0 || neighbor == 0)
    {
    vtkErrorMacro ("Cannot process undefined blocks");
    return;
    }

  vtkMultiProcessController* controller = vtkMultiProcessController::GetGlobalController ();
  int myProc = controller->GetLocalProcessId ();
  // int numProcs = controller->GetNumberOfProcesses ();

  if (block->ProcessId != myProc && neighbor->ProcessId != myProc)
    {
    return;
    }

  if (block->ProcessId != neighbor->ProcessId) 
    {
    int levelId, blockId, processId;
    if (neighbor->ProcessId != myProc) 
      {
      levelId = block->Level;
      blockId = block->BlockId;
      processId = neighbor->ProcessId;
      }
    else /* if (block->ProcessId != myProc) */
      {
      levelId = neighbor->Level;
      blockId = neighbor->BlockId;
      processId = block->ProcessId;
      }

    // Add this neighbor to this block's neighbor list.
    bool contained = false;
    for (size_t i = 0; i < this->NeighborList[levelId][blockId].size (); i ++)
      {
      if (this->NeighborList[levelId][blockId][i] == processId)
        {
        contained = true;
        break;
        }
      }
    if (!contained)
      {
      this->NeighborList[levelId][blockId].push_back (processId);
      this->ValidNeighbor[processId] = true;
      }
    }

  if (block->ProcessId == myProc) 
    {
    vtkUniformGrid* grid = volume->GetDataSet (block->Level, block->BlockId);
    vtkIdTypeArray* regionId = vtkIdTypeArray::SafeDownCast (
                                  grid->GetCellData ()->GetArray (this->RegionName.c_str()));

    int extent[6];
    grid->GetExtent (extent);

    // move inside the overlap areas.
    extent[dir*2+1] --; 

    // we only want the plane at the edge in the direction
    extent[dir*2] = extent[dir*2+1] - 1;

    if (block->Level < neighbor->Level) 
      {
      // we only want the quarter of the boundary that applies to this neighbor
      int del[3];
      del[0] = (extent[1] - extent[0]) >> 1;
      del[1] = (extent[3] - extent[2]) >> 1;
      del[2] = (extent[5] - extent[4]) >> 1;

      for (int i = 0; i < 3; i ++)
        {
        if (i == dir) 
          { 
          continue; 
          }
        // shrink the orthogonal directions by half depending on which offset the neighbor is
        if ((neighbor->GridIndex[i] % 2) == 0) 
          {
            extent[i*2+1] -= del[i];
          } else {
            extent[i*2+0] += del[i];
          }
        }
      }

    vtkSmartPointer<vtkIdTypeArray> array = vtkSmartPointer<vtkIdTypeArray>::New ();
    array->SetNumberOfComponents (1);
    array->SetNumberOfTuples (9 + 
                              (extent[1] - extent[0]) *
                              (extent[3] - extent[2]) *
                              (extent[5] - extent[4]));

    int index = 0;
    array->SetTuple1 (index, block->Level);
    index ++;
    array->SetTuple1 (index, block->GridIndex[0]);
    index ++;
    array->SetTuple1 (index, block->GridIndex[1]);
    index ++;
    array->SetTuple1 (index, block->GridIndex[2]);
    index ++;
    array->SetTuple1 (index, dir);
    index ++;
    array->SetTuple1 (index, neighbor->Level);
    index ++;
    array->SetTuple1 (index, neighbor->GridIndex[0]);
    index ++;
    array->SetTuple1 (index, neighbor->GridIndex[1]);
    index ++;
    array->SetTuple1 (index, neighbor->GridIndex[2]);
    index ++;

    if (block->Level > neighbor->Level)
      { 
      int ijk[3];
      for (ijk[0] = extent[0]; ijk[0] < extent[1]; ijk[0] += 2)
        {
        for (ijk[1] = extent[2]; ijk[1] < extent[3]; ijk[1] += 2)
          {
          for (ijk[2] = extent[4]; ijk[2] < extent[5]; ijk[2] += 2)
            {
            vtkIdType cellId = grid->ComputeCellId (ijk); 
            int region = regionId->GetTuple1 (cellId);
            int setting = 1;
            // if this cell doesn't have a region
            // check one of the other 3 that align with the neighbor cell
            while (region == 0 && setting <= 8) 
              {
              int ijk_off[3];
              ijk_off[0] = ijk[0];
              if ((setting & 1) == 1 && extent[1] != extent[0]) 
                {
                ijk_off[0] ++;
                } 
              ijk_off[1] = ijk[1];
              if ((setting & 2) == 1 && extent[3] != extent[2]) 
                {
                ijk_off[1] ++;
                } 
              ijk_off[2] = ijk[2];
              if ((setting & 4) == 1 && extent[5] != extent[4]) 
                {
                ijk_off[2] ++;
                } 
              cellId = grid->ComputeCellId (ijk_off); 
              region = regionId->GetTuple1 (cellId);
              setting ++;
              }
            array->SetTuple1 (index, region);
            index ++;
            }
          }
      }
      }
    else
      {
      int ijk[3];
      for (ijk[0] = extent[0]; ijk[0] < extent[1]; ijk[0] ++)
        {
        for (ijk[1] = extent[2]; ijk[1] < extent[3]; ijk[1] ++)
          {
          for (ijk[2] = extent[4]; ijk[2] < extent[5]; ijk[2] ++)
            {
            vtkIdType cellId = grid->ComputeCellId (ijk); 
            array->SetTuple1 (index, regionId->GetTuple1 (cellId));
            index ++;
            }
          }
        }
      }

    this->BoundaryArrays[neighbor->ProcessId].push_back (array);
    }
  else if (neighbor->ProcessId == myProc)
    {
    vtkUniformGrid* grid = volume->GetDataSet (neighbor->Level, neighbor->BlockId);
    // estimate the size of the receive by the extent of this neighbor
    int extent[6];
    grid->GetExtent (extent);
    extent[dir*2 + 1] = extent[dir*2] + 2;

    this->ReceiveList[block->ProcessId].push_back (9 + 
                                                   (extent[1] - extent[0]) *
                                                   (extent[3] - extent[2]) *
                                                   (extent[5] - extent[4]));
    }
}
//----------------------------------------------------------------------------
int vtkAMRConnectivity::ExchangeBoundaries (vtkMPIController *controller)
{
  if (controller == 0)
    {
    vtkErrorMacro ("vtkAMRConnectivity only works parallel in MPI environment");
    return 0;
    }

#ifdef PARAVIEW_USE_MPI
  int myProc = controller->GetLocalProcessId ();

  vtkAMRConnectivityCommRequestList receiveList;
  for (size_t i = 0; i < this->ReceiveList.size (); i ++) 
    {
    if (i == static_cast<unsigned int> (myProc)) { continue; }

    int messageLength = 0; 
    for (size_t j = 0; j < this->ReceiveList[i].size (); j ++)
      {
      messageLength += 1 + this->ReceiveList[i][j];
      }
    if (messageLength == 0) { continue; }
    messageLength ++; // end of message marker

    vtkSmartPointer<vtkIntArray> array = vtkSmartPointer<vtkIntArray>::New ();
    array->SetNumberOfComponents (1);
    array->SetNumberOfTuples (messageLength);

    vtkAMRConnectivityCommRequest request;
    request.SendProcess = static_cast<int>(i);
    request.ReceiveProcess = myProc;
    request.Buffer = array;

    controller->NoBlockReceive(array->GetPointer (0),
                               messageLength,
                               static_cast<int>(i), 
                               BOUNDARY_TAG,
                               request.Request);

    this->ReceiveList[i].clear ();
    receiveList.push_back(request);
    } 

  vtkAMRConnectivityCommRequestList sendList;
  for (size_t i = 0; i < this->BoundaryArrays.size (); i ++)
    {
    if (i == static_cast<unsigned int> (myProc)) { continue; }
    vtkSmartPointer<vtkIntArray> array = vtkSmartPointer<vtkIntArray>::New ();
    array->SetNumberOfComponents (1);
    array->SetNumberOfTuples (0);
    for (size_t j = 0; j < this->BoundaryArrays[i].size (); j ++)
      {
      int tuples = this->BoundaryArrays[i][j]->GetNumberOfTuples ();
      array->InsertNextTuple1 (tuples);
      for (int k = 0; k < tuples; k ++) 
        {
        array->InsertNextTuple1 (this->BoundaryArrays[i][j]->GetTuple1 (k));
        }
      }

    if (array->GetNumberOfTuples () == 0) { continue; }
    
    array->InsertNextTuple1 (-1);

    vtkAMRConnectivityCommRequest request;
    request.SendProcess = myProc;
    request.ReceiveProcess = static_cast<int>(i);
    request.Buffer = array;

    controller->NoBlockSend(array->GetPointer(0),
                            array->GetNumberOfTuples (),
                            static_cast<int>(i), 
                            BOUNDARY_TAG,
                            request.Request);

    sendList.push_back(request);
    }

  while (!receiveList.empty())
    {
    vtkAMRConnectivityCommRequest request = receiveList.WaitAny();
    vtkSmartPointer<vtkIntArray> array = request.Buffer;
    int total = array->GetNumberOfTuples ();
    int index = 0;
    while (index < total)
      {
      int tuples = array->GetTuple1 (index);
      if (tuples < 0) 
        {
        break;
        }
      index ++;
      vtkIdTypeArray* boundary = vtkIdTypeArray::New ();
      boundary->SetNumberOfComponents (1);
      boundary->SetNumberOfTuples (tuples);
      for (int i = 0; i < tuples; i ++) 
        {
        boundary->SetTuple1 (i, array->GetTuple1 (index));
        index ++;
        }
      this->BoundaryArrays[myProc].push_back (boundary);
      }
    }

  sendList.WaitAll();
  sendList.clear ();
#endif /* PARAVIEW_USE_MPI */
  return 1;
}

int vtkAMRConnectivity::ExchangeEquivPairs (vtkMPIController *controller)
{
  if (controller == 0)
    {
    vtkErrorMacro ("vtkAMRConnectivity only works parallel in MPI environment");
    return 0;
    }
#ifdef PARAVIEW_USE_MPI
  int myProc = controller->GetLocalProcessId ();
  int numProcs = controller->GetNumberOfProcesses ();

  vtkAMRConnectivityCommRequestList receiveList;
  for (size_t i = 0; i < static_cast<unsigned int> (numProcs); i ++) 
    {
    if (i == static_cast<unsigned int> (myProc) || !this->ValidNeighbor[i]) 
      {
      continue; 
      }

    vtkSmartPointer<vtkIntArray> array = vtkSmartPointer<vtkIntArray>::New ();
    array->SetNumberOfComponents (1);
    array->SetNumberOfTuples (1);

    vtkAMRConnectivityCommRequest request;
    request.SendProcess = static_cast<int>(i);
    request.ReceiveProcess = myProc;
    request.Buffer = array;

    controller->NoBlockReceive(array->GetPointer (0),
                               1,
                               static_cast<int>(i), 
                               EQUIV_SIZE_TAG,
                               request.Request);

    this->ReceiveList[i].clear ();
    receiveList.push_back(request);
    }

  vtkAMRConnectivityCommRequestList sendList;
  for (size_t i = 0; i < static_cast<unsigned int> (numProcs); i ++)
    {
    if (i == static_cast<unsigned int> (myProc) || !this->ValidNeighbor[i])  
      { 
      continue; 
      }

    vtkSmartPointer<vtkIntArray> array = vtkSmartPointer<vtkIntArray>::New ();
    array->SetNumberOfComponents (1);
    array->SetNumberOfTuples (1);
    int size = this->EquivPairs[i] == 0 ? 0 : this->EquivPairs[i]->GetNumberOfTuples ();
    array->SetTuple1 (0, size);

    vtkAMRConnectivityCommRequest request;
    request.SendProcess = myProc;
    request.ReceiveProcess = static_cast<int>(i);
    request.Buffer = array;

    controller->NoBlockSend(array->GetPointer (0),
                            1,
                            static_cast<int>(i), 
                            EQUIV_SIZE_TAG,
                            request.Request);

    sendList.push_back(request);
    }

  std::vector<int> receive_sizes (numProcs);
  for (int i = 0; i < numProcs; i ++)
    {
    receive_sizes[i] = 0;
    }

  while (!receiveList.empty())
    {
    vtkAMRConnectivityCommRequest request = receiveList.WaitAny();
    vtkSmartPointer<vtkIntArray> array = request.Buffer;
    receive_sizes[request.SendProcess] = array->GetTuple1 (0);
    }
  
  receiveList.clear ();
  sendList.WaitAll();
  sendList.clear ();

  for (size_t i = 0; i < static_cast<unsigned int>(numProcs); i ++) 
    {
    if (i == static_cast<unsigned int> (myProc) || receive_sizes[i] == 0) 
      {
      continue; 
      }

    vtkSmartPointer<vtkIntArray> array = vtkSmartPointer<vtkIntArray>::New ();
    array->SetNumberOfComponents (1);
    array->SetNumberOfTuples (receive_sizes[i]);

    vtkAMRConnectivityCommRequest request;
    request.SendProcess = static_cast<int>(i);
    request.ReceiveProcess = myProc;
    request.Buffer = array;

    controller->NoBlockReceive(array->GetPointer (0),
                               receive_sizes[i],
                               static_cast<int>(i), 
                               EQUIV_TAG,
                               request.Request);

    this->ReceiveList[i].clear ();
    receiveList.push_back(request);
    }

  for (size_t i = 0; i < static_cast<unsigned int>(numProcs); i ++)
    {
    if (i == static_cast<unsigned int> (myProc) 
        || this->EquivPairs[i] == 0
        || this->EquivPairs[i]->GetNumberOfTuples () == 0)  
      { 
      continue; 
      }

    vtkAMRConnectivityCommRequest request;
    request.SendProcess = myProc;
    request.ReceiveProcess = static_cast<int>(i);
    request.Buffer = this->EquivPairs[i];

    controller->NoBlockSend(this->EquivPairs[i]->GetPointer (0),
                            this->EquivPairs[i]->GetNumberOfTuples (),
                            static_cast<int>(i), 
                            EQUIV_TAG,
                            request.Request);

    sendList.push_back(request);
    }

  while (!receiveList.empty())
    {
    vtkAMRConnectivityCommRequest request = receiveList.WaitAny();
    vtkSmartPointer<vtkIntArray> array = request.Buffer;
    for (int i = 0; i < array->GetNumberOfTuples (); i += 2)  
      {
      int v1 = array->GetValue (i);
      int v2 = array->GetValue (i+1);
      this->Equivalence->AddEquivalence (v1, v2);
      }
    }

  receiveList.clear ();
  sendList.WaitAll();
  sendList.clear ();
#endif /* PARAVIEW_USE_MPI */
  return 1;
}

//----------------------------------------------------------------------------
void vtkAMRConnectivity::ProcessBoundaryAtNeighbor (
                vtkNonOverlappingAMR* volume, vtkIdTypeArray *array)
{

  vtkMultiProcessController* controller = vtkMultiProcessController::GetGlobalController ();
  int myProc = controller->GetLocalProcessId ();
  // int numProcs = controller->GetNumberOfProcesses ();

  int index = 0;
  int blockLevel = array->GetTuple1 (index);
  index ++;
  int blockGrid[3];
  blockGrid[0] = array->GetTuple1 (index);
  index ++;
  blockGrid[1] = array->GetTuple1 (index);
  index ++;
  blockGrid[2] = array->GetTuple1 (index);
  index ++;
  int dir = array->GetTuple1 (index);
  index ++;
  int neighborLevel = array->GetTuple1 (index);
  index ++;
  int gridX = array->GetTuple1 (index);
  index ++;
  int gridY = array->GetTuple1 (index);
  index ++;
  int gridZ = array->GetTuple1 (index);
  index ++;

  vtkAMRDualGridHelperBlock* neighbor = 
          this->Helper->GetBlock (neighborLevel, gridX, gridY, gridZ);

  if (neighbor->ProcessId != myProc) 
    {
    vtkErrorMacro ("Shouldn't get here.  Cannot process neighbor on remote process");
    return;
    }

  vtkUniformGrid* grid = volume->GetDataSet (neighbor->Level, neighbor->BlockId);
  vtkIdTypeArray* regionId = vtkIdTypeArray::SafeDownCast (
                                grid->GetCellData ()->GetArray (this->RegionName.c_str()));

  int extent[6];
  grid->GetExtent (extent);

  // move inside the overlap areas.
  extent[dir*2] ++; 

  // we only want the plane at the edge in dir
  extent[dir*2+1] = extent[dir*2] + 1;

  if (blockLevel > neighborLevel) 
    {
    // we only want the quarter of the boundary that applies to this neighbor
    int del[3];
    del[0] = (extent[1] - extent[0]) >> 1;
    del[1] = (extent[3] - extent[2]) >> 1;
    del[2] = (extent[5] - extent[4]) >> 1;

    for (int i = 0; i < 3; i ++)
      {
      if (i == dir) 
        { 
        continue; 
        }
      // shrink the orthogonal directions by half depending on which offset the neighbor is
      // TODO need the blocks offset within the octet
      if ((blockGrid[i] % 2) == 0) 
        {
          extent[i*2+1] -= del[i];
        } else {
          extent[i*2+0] += del[i];
        }
      }
    }

  if (blockLevel < neighborLevel)
    {
    int ijk[3];
    for (ijk[0] = extent[0]; ijk[0] < extent[1]; ijk[0] += 2)
      {
      for (ijk[1] = extent[2]; ijk[1] < extent[3]; ijk[1] += 2)
        {
        for (ijk[2] = extent[4]; ijk[2] < extent[5]; ijk[2] += 2)
          {
          vtkIdType cellId = grid->ComputeCellId (ijk); 
          int neighborRegion = regionId->GetTuple1 (cellId);
          int setting = 1;
          // if this cell doesn't have a region
          // check one of the other 3 that align with the neighbor cell
          while (neighborRegion == 0 && setting <= 8) 
            {
            int ijk_off[3];
            ijk_off[0] = ijk[0];
            if ((setting & 1) == 1 && extent[1] != extent[0]) 
              {
              ijk_off[0] ++;
              } 
            ijk_off[1] = ijk[1];
            if ((setting & 2) == 1 && extent[3] != extent[2]) 
              {
              ijk_off[1] ++;
              } 
            ijk_off[2] = ijk[2];
            if ((setting & 4) == 1 && extent[5] != extent[4]) 
              {
              ijk_off[2] ++;
              } 
            cellId = grid->ComputeCellId (ijk_off); 
            neighborRegion = regionId->GetTuple1 (cellId);
            setting ++;
            }
          int blockRegion = array->GetTuple1 (index);
          if (neighborRegion != 0 && blockRegion != 0) 
            {
            this->Equivalence->AddEquivalence (neighborRegion, blockRegion);
            }
          }
          index ++;
        }
      }
    }
  else /* blockLevel >= neighborLevel */
    {
    int ijk[3];
    for (ijk[0] = extent[0]; ijk[0] < extent[1]; ijk[0] ++)
      {
      for (ijk[1] = extent[2]; ijk[1] < extent[3]; ijk[1] ++)
        {
        for (ijk[2] = extent[4]; ijk[2] < extent[5]; ijk[2] ++)
          {
          vtkIdType cellId = grid->ComputeCellId (ijk); 
          int neighborRegion = regionId->GetTuple1 (cellId);
          int blockRegion = array->GetTuple1 (index);
          if (neighborRegion != 0 && blockRegion != 0) 
            {
            this->Equivalence->AddEquivalence (neighborRegion, blockRegion);
            }
          index ++;
          }
        }
      }
    }
}
