// SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
// SPDX-License-Identifier: BSD-3-Clause
#include "vtkStructuredGridConnectivity.h"

// VTK includes
#include "vtkIdList.h"
#include "vtkObjectFactory.h"
#include "vtkStructuredData.h"
#include "vtkStructuredExtent.h"
#include "vtkStructuredNeighbor.h"
#include "vtkUnsignedCharArray.h"

#include <iostream>

VTK_ABI_NAMESPACE_BEGIN
vtkStandardNewMacro(vtkStructuredGridConnectivity);
VTK_ABI_NAMESPACE_END

// Description:
// An enum to define the 6 block faces
namespace BlockFace
{
VTK_ABI_NAMESPACE_BEGIN
enum
{
  FRONT = 0,  // (+k direction)
  BACK = 1,   // (-k direction)
  RIGHT = 2,  // (+i direction)
  LEFT = 3,   // (-i direction)
  TOP = 4,    // (+j direction)
  BOTTOM = 5, // (-j direction)
  NOT_ON_BLOCK_FACE = 6
};
VTK_ABI_NAMESPACE_END
}

VTK_ABI_NAMESPACE_BEGIN
//------------------------------------------------------------------------------
vtkStructuredGridConnectivity::vtkStructuredGridConnectivity()
{
  this->DataDimension = 0;
  this->DataDescription = -1;
  this->NumberOfGrids = 0;
  this->WholeExtent[0] = this->WholeExtent[1] = this->WholeExtent[2] = this->WholeExtent[3] =
    this->WholeExtent[4] = this->WholeExtent[5] = -1;
  this->GridExtents.clear();
}

//------------------------------------------------------------------------------
vtkStructuredGridConnectivity::~vtkStructuredGridConnectivity()
{
  this->GridExtents.clear();
  this->NeighborPair2NeighborListIndex.clear();
}

//------------------------------------------------------------------------------
void vtkStructuredGridConnectivity::PrintSelf(std::ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);

  os << "========================\n";
  os << "DATA DIMENSION: " << this->DataDimension << std::endl;
  os << "WHOLE EXTENT: [ ";
  for (int i = 0; i < 6; i++)
  {
    os << this->WholeExtent[i] << " ";
  }
  os << "]\n";
  os << "CONNECTIVITY INFORMATION: \n";
  for (unsigned int gridID = 0; gridID < this->NumberOfGrids; ++gridID)
  {
    int GridExtent[6];
    int RealExtent[6];
    this->GetGridExtent(gridID, GridExtent);
    this->GetRealExtent(gridID, GridExtent, RealExtent);
    os << "GRID[ " << gridID << "]: ";
    for (int i = 0; i < 6; i += 2)
    {
      os << " [";
      os << GridExtent[i] << ", " << GridExtent[i + 1] << "]";
    }
    os << " REAL EXTENT: ";
    for (int i = 0; i < 6; i += 2)
    {
      os << " [";
      os << RealExtent[i] << ", " << RealExtent[i + 1] << "]";
    }
    os << std::endl;
    os << " Connecting faces: " << this->GetNumberOfConnectingBlockFaces(gridID) << " ";

    os << "[ ";
    if (this->HasBlockConnection(gridID, BlockFace::FRONT))
    {
      os << "FRONT(+k) ";
    }
    if (this->HasBlockConnection(gridID, BlockFace::BACK))
    {
      os << "BACK(-k) ";
    }
    if (this->HasBlockConnection(gridID, BlockFace::RIGHT))
    {
      os << "RIGHT(+i) ";
    }
    if (this->HasBlockConnection(gridID, BlockFace::LEFT))
    {
      os << "LEFT(-i) ";
    }
    if (this->HasBlockConnection(gridID, BlockFace::TOP))
    {
      os << "TOP(+j) ";
    }
    if (this->HasBlockConnection(gridID, BlockFace::BOTTOM))
    {
      os << "BOTTOM(-j) ";
    }
    os << "] ";
    os << std::endl;

    for (unsigned int nei = 0; nei < this->Neighbors[gridID].size(); ++nei)
    {
      int NeiExtent[6];
      this->GetGridExtent(this->Neighbors[gridID][nei].NeighborID, NeiExtent);

      os << "\t N[" << nei << "] GRID ID:" << this->Neighbors[gridID][nei].NeighborID << " ";
      for (int i = 0; i < 6; i += 2)
      {
        os << " [";
        os << NeiExtent[i] << ", " << NeiExtent[i + 1] << "] ";
      }

      os << " overlaps @ ";
      for (int i = 0; i < 6; i += 2)
      {
        os << " [";
        os << this->Neighbors[gridID][nei].OverlapExtent[i] << ", ";
        os << this->Neighbors[gridID][nei].OverlapExtent[i + 1] << "] ";
      }

      os << " orientation: (";
      os << this->Neighbors[gridID][nei].Orientation[0] << ", ";
      os << this->Neighbors[gridID][nei].Orientation[1] << ", ";
      os << this->Neighbors[gridID][nei].Orientation[2] << ")\n ";
      os << std::endl;

      os << "\t RCVEXTENT: ";
      for (int i = 0; i < 6; i += 2)
      {
        os << " [";
        os << this->Neighbors[gridID][nei].RcvExtent[i] << ", ";
        os << this->Neighbors[gridID][nei].RcvExtent[i + 1] << "] ";
      }
      os << std::endl;

      os << "\t SNDEXTENT: ";
      for (int i = 0; i < 6; i += 2)
      {
        os << " [";
        os << this->Neighbors[gridID][nei].SendExtent[i] << ", ";
        os << this->Neighbors[gridID][nei].SendExtent[i + 1] << "] ";
      }
      os << std::endl << std::endl;

    } // END for all neis
  }   // END for all grids
}

//------------------------------------------------------------------------------
void vtkStructuredGridConnectivity::RegisterGrid(int gridID, int ext[6],
  vtkUnsignedCharArray* nodesGhostArray, vtkUnsignedCharArray* cellGhostArray,
  vtkPointData* pointData, vtkCellData* cellData, vtkPoints* gridNodes)
{
  assert("pre: gridID out-of-bounds!" &&
    (gridID >= 0 && gridID < static_cast<int>(this->NumberOfGrids)));

  for (int i = 0; i < 6; ++i)
  {
    this->GridExtents[gridID * 6 + i] = ext[i];
  }

  this->RegisterGridGhostArrays(gridID, nodesGhostArray, cellGhostArray);
  this->RegisterFieldData(gridID, pointData, cellData);
  this->RegisterGridNodes(gridID, gridNodes);
}

//------------------------------------------------------------------------------
void vtkStructuredGridConnectivity::AcquireDataDescription()
{
  if (this->DataDescription != -1)
  {
    return;
  }

  int dims[3];

  vtkStructuredExtent::GetDimensions(this->WholeExtent, dims);

  this->DataDescription = vtkStructuredData::GetDataDescription(dims);
  this->DataDimension = vtkStructuredData::GetDataDimension(this->DataDescription);

  assert("pre: Error acquiring data description" && (this->DataDescription >= 0));
  assert("pre: grid description cannot be empty" &&
    (this->DataDescription != vtkStructuredData::VTK_STRUCTURED_EMPTY));
}

//------------------------------------------------------------------------------
vtkStructuredNeighbor vtkStructuredGridConnectivity::GetGridNeighbor(int gridID, int nei)
{
  assert("pre: gridID out-of-bounds!" &&
    (gridID >= 0 && gridID < static_cast<int>(this->NumberOfGrids)));
  assert(
    "pre: nei index is out-of-bounds!" && (nei >= 0) && (nei < this->GetNumberOfNeighbors(gridID)));

  return (this->Neighbors[gridID][nei]);
}

//------------------------------------------------------------------------------
vtkIdList* vtkStructuredGridConnectivity::GetNeighbors(int gridID, int* extents)
{
  assert("pre: input extents array is nullptr" && (extents != nullptr));

  int N = this->GetNumberOfNeighbors(gridID);
  if (N < 1)
  {
    return nullptr;
  }

  vtkIdList* neiList = vtkIdList::New();
  neiList->SetNumberOfIds(N);

  unsigned int nei = 0;
  for (; nei < this->Neighbors[gridID].size(); ++nei)
  {
    vtkIdType neiId = this->Neighbors[gridID][nei].NeighborID;
    neiList->SetId(nei, neiId);
    for (int i = 0; i < 6; ++i)
    {
      extents[nei * 6 + i] = this->Neighbors[gridID][nei].OverlapExtent[i];
    }
  } // END for all neighbors

  assert("post: N==neiList.size()" && (N == neiList->GetNumberOfIds()));
  return (neiList);
}

//------------------------------------------------------------------------------
void vtkStructuredGridConnectivity::ComputeNeighbors()
{
  // STEP 0: Acquire data description, i.e., determine how the structured data
  // is laid out, e.g., is it volumetric or 2-D along some plane, XY, XZ, or YZ.
  this->AcquireDataDescription();
  if (this->DataDescription == vtkStructuredData::VTK_STRUCTURED_EMPTY ||
    this->DataDescription == vtkStructuredData::VTK_STRUCTURED_SINGLE_POINT)
  {
    return;
  }

  // STEP 1: Establish neighbors based on the structured extents.
  for (unsigned int i = 0; i < this->NumberOfGrids; ++i)
  {
    this->SetBlockTopology(i);
    for (unsigned int j = i + 1; j < this->NumberOfGrids; ++j)
    {
      this->EstablishNeighbors(i, j);
    } // END for all j
  }   // END for all i

  // STEP 2: Fill the ghost arrays
  for (unsigned int i = 0; i < this->NumberOfGrids; ++i)
  {
    // NOTE: typically remote grids have nullptr ghost arrays, by this approach
    // ComputeNeighbors() can be called transparently from
    // vtkPStructuredGridConnectivity without any modification.
    if (this->GridPointGhostArrays[i] != nullptr)
    {
      this->FillGhostArrays(i, this->GridPointGhostArrays[i], this->GridCellGhostArrays[i]);
    }
  } // END for all grids
}

//------------------------------------------------------------------------------
void vtkStructuredGridConnectivity::SearchNeighbors(
  int gridID, int i, int j, int k, vtkIdList* neiList)
{
  assert("pre: neiList should not be nullptr" && (neiList != nullptr));
  assert("pre: gridID is out-of-bounds" &&
    ((gridID >= 0) && (gridID < static_cast<int>(this->NumberOfGrids))));

  for (unsigned int nei = 0; nei < this->Neighbors[gridID].size(); ++nei)
  {
    vtkStructuredNeighbor* myNei = &this->Neighbors[gridID][nei];
    assert("pre: myNei != nullptr" && (myNei != nullptr));

    if (this->IsNodeWithinExtent(i, j, k, myNei->OverlapExtent))
    {
      neiList->InsertNextId(myNei->NeighborID);
    }
  } // END for all neis
}

//------------------------------------------------------------------------------
void vtkStructuredGridConnectivity::MarkCellProperty(
  unsigned char& pfield, unsigned char* nodeGhostFields, int numNodes)
{
  // Sanity check
  assert("pre: node ghostfields should not be nullptr" && (nodeGhostFields != nullptr));

  pfield = 0;

  for (int i = 0; i < numNodes; ++i)
  {
    if (nodeGhostFields[i] & vtkDataSetAttributes::DUPLICATEPOINT)
    {
      pfield |= vtkDataSetAttributes::DUPLICATECELL;
      return;
    }
  } // END for all nodes
}

//------------------------------------------------------------------------------
void vtkStructuredGridConnectivity::MarkNodeProperty(
  int gridID, int i, int j, int k, int ext[6], int realExtent[6], unsigned char& p)
{
  p = 0;

  // Check if the node is an interior a node, i.e., it is not on any boundary
  // shared or real boundary and not in a ghost region. Interior nodes can only
  // be internal nodes!
  if (!this->IsNodeInterior(i, j, k, realExtent))
  {
    // If the node is on the boundary of the computational domain mark it
    // if( this->IsNodeOnBoundary(i,j,k) )
    // {
    //   BOUNDARY might be used in the future
    // }

    // Check if the node is also on a shared boundary or if it is a ghost node
    if (this->IsNodeOnSharedBoundary(gridID, realExtent, i, j, k))
    {
      // SHARED might be used in the future

      // For shared nodes we must check for ownership
      vtkIdList* neiList = vtkIdList::New();
      this->SearchNeighbors(gridID, i, j, k, neiList);

      if (neiList->GetNumberOfIds() > 0)
      {
        int neiRealExtent[6];
        int neiGridExtent[6];

        for (vtkIdType nei = 0; nei < neiList->GetNumberOfIds(); ++nei)
        {
          vtkIdType neiIdx = neiList->GetId(nei);
          this->GetGridExtent(neiIdx, neiGridExtent);
          this->GetRealExtent(neiIdx, neiGridExtent, neiRealExtent);

          // If my gridID is not the smallest gridID that shares the point,then
          // I don't own the point.
          // The convention is that the grid with the smallest gridID will own the
          // point and all other grids should IGNORE it when computing statistics
          // etc.
          if (this->IsNodeWithinExtent(i, j, k, neiRealExtent) && gridID > neiList->GetId(nei))
          {
            p |= vtkDataSetAttributes::HIDDENPOINT;
            break;
          }
        } // END for all neis
      }   // END if neisList isn't empty
      neiList->Delete();
    } // END if node is on a shared boundary
    else if (this->IsGhostNode(ext, realExtent, i, j, k))
    {
      p |= vtkDataSetAttributes::DUPLICATEPOINT;
    }
  }
}

//------------------------------------------------------------------------------
void vtkStructuredGridConnectivity::FillNodesGhostArray(int gridID, int dataDescription,
  int GridExtent[6], int RealExtent[6], vtkUnsignedCharArray* nodesArray)
{
  int ijk[3];
  for (int i = GridExtent[0]; i <= GridExtent[1]; ++i)
  {
    for (int j = GridExtent[2]; j <= GridExtent[3]; ++j)
    {
      for (int k = GridExtent[4]; k <= GridExtent[5]; ++k)
      {
        ijk[0] = i;
        ijk[1] = j;
        ijk[2] = k;
        vtkIdType idx =
          vtkStructuredData::ComputePointIdForExtent(GridExtent, ijk, dataDescription);

        this->MarkNodeProperty(
          gridID, i, j, k, GridExtent, RealExtent, *nodesArray->GetPointer(idx));
      } // END for all k
    }   // END for all j
  }     // END for all i
}

//------------------------------------------------------------------------------
void vtkStructuredGridConnectivity::FillCellsGhostArray(int dataDescription, int numNodesPerCell,
  int dims[3], int CellExtent[6], vtkUnsignedCharArray* nodesArray,
  vtkUnsignedCharArray* cellsArray)
{
  assert("pre: nodes array should not be nullptr" && (nodesArray != nullptr));

  if (cellsArray == nullptr)
  {
    return;
  }

  vtkIdList* cellNodeIds = vtkIdList::New();
  unsigned char* cellNodeGhostFields = new unsigned char[numNodesPerCell];

  int ijk[3];
  for (int i = CellExtent[0]; i <= CellExtent[1]; ++i)
  {
    for (int j = CellExtent[2]; j <= CellExtent[3]; ++j)
    {
      for (int k = CellExtent[4]; k <= CellExtent[5]; ++k)
      {
        ijk[0] = i;
        ijk[1] = j;
        ijk[2] = k;

        // Note: this is really a cell index, since it is computed from the
        // cell extent
        vtkIdType idx =
          vtkStructuredData::ComputePointIdForExtent(CellExtent, ijk, dataDescription);

        cellNodeIds->Reset();
        vtkStructuredData::GetCellPoints(idx, cellNodeIds, dataDescription, dims);
        assert(cellNodeIds->GetNumberOfIds() == numNodesPerCell);
        assert(cellNodeGhostFields != nullptr);

        for (int ii = 0; ii < numNodesPerCell; ++ii)
        {
          vtkIdType nodeIdx = cellNodeIds->GetId(ii);
          cellNodeGhostFields[ii] = *nodesArray->GetPointer(nodeIdx);
        } // END for all nodes

        this->MarkCellProperty(*cellsArray->GetPointer(idx), cellNodeGhostFields, numNodesPerCell);
      } // END for all cells along k
    }   // END for all cells along j
  }     // END for all cells along i

  delete[] cellNodeGhostFields;
  cellNodeIds->Delete();
}

//------------------------------------------------------------------------------
void vtkStructuredGridConnectivity::FillGhostArrays(
  int gridID, vtkUnsignedCharArray* nodesArray, vtkUnsignedCharArray* cellsArray)
{
  if (nodesArray == nullptr)
  {
    return;
  }

  // STEP 0: Get the grid information
  int GridExtent[6];
  this->GetGridExtent(gridID, GridExtent);

  // STEP 1: Real extent
  int RealExtent[6];
  this->GetRealExtent(gridID, GridExtent, RealExtent);

  // STEP 2: Get the data description
  int dataDescription = vtkStructuredData::GetDataDescriptionFromExtent(GridExtent);

  // STEP 3: Get the cell extent
  int CellExtent[6];
  vtkStructuredData::GetCellExtentFromPointExtent(GridExtent, CellExtent, dataDescription);

  // STEP 4: Get the data dimension
  int dim = vtkStructuredData::GetDataDimension(dataDescription);
  assert("pre: data dimensions must be 1, 2 or 3" && (dim >= 1) && (dim <= 3));

  // STEP 5: Get the grid dimensions from the given extent.
  int dims[3];
  vtkStructuredData::GetDimensionsFromExtent(GridExtent, dims);

  // STEP 6: Get the number of nodes per cell
  int numNodes = this->GetNumberOfNodesPerCell(dim);

  // STEP 7: Mark nodes
  this->FillNodesGhostArray(gridID, dataDescription, GridExtent, RealExtent, nodesArray);

  // STEP 8: Mark Cells
  this->FillCellsGhostArray(dataDescription, numNodes, dims, CellExtent, nodesArray, cellsArray);
}

//------------------------------------------------------------------------------
void vtkStructuredGridConnectivity::GetRealExtent(int gridID, int GridExtent[6], int RealExtent[6])
{
  for (int i = 0; i < 6; ++i)
  {
    RealExtent[i] = GridExtent[i];
  }

  if (this->NumberOfGhostLayers == 0)
  {
    return;
  }

  switch (this->DataDescription)
  {
    case vtkStructuredData::VTK_STRUCTURED_X_LINE:
      if (this->HasBlockConnection(gridID, BlockFace::LEFT))
      {
        RealExtent[0] += this->NumberOfGhostLayers; // imin
      }
      if (this->HasBlockConnection(gridID, BlockFace::RIGHT))
      {
        RealExtent[1] -= this->NumberOfGhostLayers; // imax
      }
      break;
    case vtkStructuredData::VTK_STRUCTURED_Y_LINE:
      if (this->HasBlockConnection(gridID, BlockFace::BOTTOM))
      {
        RealExtent[2] += this->NumberOfGhostLayers; // jmin
      }
      if (this->HasBlockConnection(gridID, BlockFace::TOP))
      {
        RealExtent[3] -= this->NumberOfGhostLayers; // jmax
      }
      break;
    case vtkStructuredData::VTK_STRUCTURED_Z_LINE:
      if (this->HasBlockConnection(gridID, BlockFace::BACK))
      {
        RealExtent[4] += this->NumberOfGhostLayers; // kmin
      }
      if (this->HasBlockConnection(gridID, BlockFace::FRONT))
      {
        RealExtent[5] -= this->NumberOfGhostLayers; // kmax
      }
      break;
    case vtkStructuredData::VTK_STRUCTURED_XY_PLANE:
      if (this->HasBlockConnection(gridID, BlockFace::LEFT))
      {
        RealExtent[0] += this->NumberOfGhostLayers; // imin
      }
      if (this->HasBlockConnection(gridID, BlockFace::RIGHT))
      {
        RealExtent[1] -= this->NumberOfGhostLayers; // imax
      }
      if (this->HasBlockConnection(gridID, BlockFace::BOTTOM))
      {
        RealExtent[2] += this->NumberOfGhostLayers; // jmin
      }
      if (this->HasBlockConnection(gridID, BlockFace::TOP))
      {
        RealExtent[3] -= this->NumberOfGhostLayers; // jmax
      }
      break;
    case vtkStructuredData::VTK_STRUCTURED_YZ_PLANE:
      if (this->HasBlockConnection(gridID, BlockFace::BOTTOM))
      {
        RealExtent[2] += this->NumberOfGhostLayers; // jmin
      }
      if (this->HasBlockConnection(gridID, BlockFace::TOP))
      {
        RealExtent[3] -= this->NumberOfGhostLayers; // jmax
      }
      if (this->HasBlockConnection(gridID, BlockFace::BACK))
      {
        RealExtent[4] += this->NumberOfGhostLayers; // kmin
      }
      if (this->HasBlockConnection(gridID, BlockFace::FRONT))
      {
        RealExtent[5] -= this->NumberOfGhostLayers; // kmax
      }
      break;
    case vtkStructuredData::VTK_STRUCTURED_XZ_PLANE:
      if (this->HasBlockConnection(gridID, BlockFace::LEFT))
      {
        RealExtent[0] += this->NumberOfGhostLayers; // imin
      }
      if (this->HasBlockConnection(gridID, BlockFace::RIGHT))
      {
        RealExtent[1] -= this->NumberOfGhostLayers; // imax
      }
      if (this->HasBlockConnection(gridID, BlockFace::BACK))
      {
        RealExtent[4] += this->NumberOfGhostLayers; // kmin
      }
      if (this->HasBlockConnection(gridID, BlockFace::FRONT))
      {
        RealExtent[5] -= this->NumberOfGhostLayers; // kmax
      }
      break;
    case vtkStructuredData::VTK_STRUCTURED_XYZ_GRID:
      if (this->HasBlockConnection(gridID, BlockFace::LEFT))
      {
        RealExtent[0] += this->NumberOfGhostLayers; // imin
      }
      if (this->HasBlockConnection(gridID, BlockFace::RIGHT))
      {
        RealExtent[1] -= this->NumberOfGhostLayers; // imax
      }
      if (this->HasBlockConnection(gridID, BlockFace::BOTTOM))
      {
        RealExtent[2] += this->NumberOfGhostLayers; // jmin
      }
      if (this->HasBlockConnection(gridID, BlockFace::TOP))
      {
        RealExtent[3] -= this->NumberOfGhostLayers; // jmax
      }
      if (this->HasBlockConnection(gridID, BlockFace::BACK))
      {
        RealExtent[4] += this->NumberOfGhostLayers; // kmin
      }
      if (this->HasBlockConnection(gridID, BlockFace::FRONT))
      {
        RealExtent[5] -= this->NumberOfGhostLayers; // kmax
      }
      break;
    default:
      std::cout << "Data description is: " << this->DataDescription << "\n";
      std::cout.flush();
      assert("pre: Undefined data-description!" && false);
  }
  vtkStructuredExtent::Clamp(RealExtent, this->WholeExtent);
}

//------------------------------------------------------------------------------
bool vtkStructuredGridConnectivity::IsNodeOnSharedBoundary(
  int gridID, int RealExtent[6], int i, int j, int k)
{
  if (this->IsNodeOnBoundaryOfExtent(i, j, k, RealExtent))
  {
    int orient[3];
    this->GetIJKBlockOrientation(i, j, k, RealExtent, orient);
    for (int ii = 0; ii < 3; ++ii)
    {
      if ((orient[ii] != BlockFace::NOT_ON_BLOCK_FACE) &&
        this->HasBlockConnection(gridID, orient[ii]))
      {
        return true;
      }
    } // END for all dimensions
    return false;
  }
  else
  {
    return false;
  }
}

//------------------------------------------------------------------------------
bool vtkStructuredGridConnectivity::IsGhostNode(
  int GridExtent[6], int RealExtent[6], int i, int j, int k)
{
  // STEP 0: Check if there are any ghost-layers. Note, if the original data
  // that the user is registering contains ghost-layers, the users must set
  // the number of ghost-layers.
  if (this->NumberOfGhostLayers == 0)
  {
    // Grid has no ghost-layers, so, the node cannot be a ghost node
    return false;
  }

  bool status = false;
  if (!this->IsNodeWithinExtent(i, j, k, RealExtent) &&
    this->IsNodeWithinExtent(i, j, k, GridExtent))
  {
    status = true;
  }
  return (status);
}

//------------------------------------------------------------------------------
bool vtkStructuredGridConnectivity::IsNodeOnBoundary(int i, int j, int k)
{
  return (this->IsNodeOnBoundaryOfExtent(i, j, k, this->WholeExtent));
}

//------------------------------------------------------------------------------
void vtkStructuredGridConnectivity::EstablishNeighbors(int i, int j)
{
  assert("pre: i < j" && (i < j));

  int iGridExtent[6];
  int jGridExtent[6];
  this->GetGridExtent(i, iGridExtent);
  this->GetGridExtent(j, jGridExtent);

  // A 3-tuple that defines the grid orientation of the form {i,j,k} where
  // i=0;, j=1, k=2. For example, let's say that we want to define the
  // orientation to be in the XZ plane, then, orientation array would be
  // constructed as follows: {0,2 -1}, where -1 indicates a NIL value.
  int orientation[3];

  // A place holder for setting up ndim, which store the dimensionality of
  // the data.
  int ndim = 3;

  switch (this->DataDescription)
  {
    case vtkStructuredData::VTK_STRUCTURED_X_LINE:
      ndim = 1;
      orientation[0] = 0;
      orientation[1] = -1;
      orientation[2] = -1;
      break;
    case vtkStructuredData::VTK_STRUCTURED_Y_LINE:
      ndim = 1;
      orientation[0] = 1;
      orientation[1] = -1;
      orientation[2] = -1;
      break;
    case vtkStructuredData::VTK_STRUCTURED_Z_LINE:
      ndim = 1;
      orientation[0] = 2;
      orientation[1] = -1;
      orientation[2] = -1;
      break;
    case vtkStructuredData::VTK_STRUCTURED_XY_PLANE:
      ndim = 2;
      orientation[0] = 0;
      orientation[1] = 1;
      orientation[2] = -1;
      break;
    case vtkStructuredData::VTK_STRUCTURED_YZ_PLANE:
      ndim = 2;
      orientation[0] = 1;
      orientation[1] = 2;
      orientation[2] = -1;
      break;
    case vtkStructuredData::VTK_STRUCTURED_XZ_PLANE:
      ndim = 2;
      orientation[0] = 0;
      orientation[1] = 2;
      orientation[2] = -1;
      break;
    case vtkStructuredData::VTK_STRUCTURED_XYZ_GRID:
      ndim = 3;
      orientation[0] = 0;
      orientation[1] = 1;
      orientation[2] = 2;
      break;
    default:
      std::cout << "Data description is: " << this->DataDescription << "\n";
      std::cout.flush();
      assert("pre: Undefined data-description!" && false);
  } // END switch

  this->DetectNeighbors(i, j, iGridExtent, jGridExtent, orientation, ndim);
}

//------------------------------------------------------------------------------
void vtkStructuredGridConnectivity::DetectNeighbors(
  int i, int j, int ex1[6], int ex2[6], int orientation[3], int ndim)
{
  std::vector<int> status;
  status.resize(ndim);

  int A[2];
  int B[2];
  int overlap[2];
  int iOrientation[3];
  int jOrientation[3];

  int overlapExtent[6];
  for (int ii = 0; ii < 3; ++ii)
  {
    overlapExtent[ii * 2] = overlapExtent[ii * 2 + 1] = 0;
    iOrientation[ii] = vtkStructuredNeighbor::UNDEFINED;
    jOrientation[ii] = vtkStructuredNeighbor::UNDEFINED;
  }

  int dim = 0;
  int idx = 0;
  for (dim = 0; dim < ndim; ++dim)
  {
    idx = orientation[dim];
    A[0] = ex1[idx * 2];
    A[1] = ex1[idx * 2 + 1];
    B[0] = ex2[idx * 2];
    B[1] = ex2[idx * 2 + 1];
    status[idx] = this->IntervalOverlap(A, B, overlap);
    if (status[idx] == VTK_NO_OVERLAP)
    {
      return; /* No neighbors */
    }

    overlapExtent[idx * 2] = overlap[0];
    overlapExtent[idx * 2 + 1] = overlap[1];

    this->DetermineNeighborOrientation(idx, A, B, overlap, iOrientation);
    this->DetermineNeighborOrientation(idx, B, A, overlap, jOrientation);
  } // END for all dimensions

  this->SetNeighbors(i, j, iOrientation, jOrientation, overlapExtent);
}

//------------------------------------------------------------------------------
void vtkStructuredGridConnectivity::SetBlockTopology(int gridID)
{
  int gridExtent[6];
  this->GetGridExtent(gridID, gridExtent);

  // Check in IMIN
  if (gridExtent[0] > this->WholeExtent[0])
  {
    this->AddBlockConnection(gridID, BlockFace::LEFT);
  }

  // Check in IMAX
  if (gridExtent[1] < this->WholeExtent[1])
  {
    this->AddBlockConnection(gridID, BlockFace::RIGHT);
  }

  // Check in JMIN
  if (gridExtent[2] > this->WholeExtent[2])
  {
    this->AddBlockConnection(gridID, BlockFace::BOTTOM);
  }

  // Check in JMAX
  if (gridExtent[3] < this->WholeExtent[3])
  {
    this->AddBlockConnection(gridID, BlockFace::TOP);
  }

  // Check in KMIN
  if (gridExtent[4] > this->WholeExtent[4])
  {
    this->AddBlockConnection(gridID, BlockFace::BACK);
  }

  // Check in KMAX
  if (gridExtent[5] < this->WholeExtent[5])
  {
    this->AddBlockConnection(gridID, BlockFace::FRONT);
  }
}

//------------------------------------------------------------------------------
void vtkStructuredGridConnectivity::SetNeighbors(
  int i, int j, int i2jOrientation[3], int j2iOrientation[3], int overlapExtent[6])
{
  vtkStructuredNeighbor Ni2j(j, overlapExtent, i2jOrientation);
  vtkStructuredNeighbor Nj2i(i, overlapExtent, j2iOrientation);

  // STEP 0: Setup i-to-j
  this->Neighbors[i].push_back(Ni2j);
  int i2jNeiIdx = static_cast<int>(this->Neighbors[i].size()) - 1;
  std::pair<int, int> i2jPair = std::make_pair(i, j);
  assert("ERROR: Duplicate neighboring pair!" &&
    this->NeighborPair2NeighborListIndex.find(i2jPair) ==
      this->NeighborPair2NeighborListIndex.end());
  this->NeighborPair2NeighborListIndex[i2jPair] = i2jNeiIdx;

  // STEP 1: Setup j-to-i
  this->Neighbors[j].push_back(Nj2i);
  int j2iNeiIdx = static_cast<int>(this->Neighbors[j].size()) - 1;
  std::pair<int, int> j2iPair = std::make_pair(j, i);
  assert("ERROR: Duplicate neighboring pair!" &&
    this->NeighborPair2NeighborListIndex.find(j2iPair) ==
      this->NeighborPair2NeighborListIndex.end());
  this->NeighborPair2NeighborListIndex[j2iPair] = j2iNeiIdx;
}

//------------------------------------------------------------------------------
void vtkStructuredGridConnectivity::PrintExtent(int ex[6])
{
  for (int i = 0; i < 3; ++i)
  {
    std::cout << " [" << ex[i * 2] << ", " << ex[i * 2 + 1] << "] ";
  }
  std::cout << std::endl;
  std::cout.flush();
}

//------------------------------------------------------------------------------
int vtkStructuredGridConnectivity::DoPartialOverlap(int s[2], int S[2], int overlap[2])
{
  if (this->InBounds(s[0], S[0], S[1]) && this->InBounds(s[1], S[0], S[1]))
  {
    overlap[0] = s[0];
    overlap[1] = s[1];
    return VTK_PARTIAL_OVERLAP;
  }
  else if (this->InBounds(s[0], S[0], S[1]))
  {
    overlap[0] = s[0];
    overlap[1] = S[1];
    if (overlap[0] == overlap[1])
    {
      return VTK_NODE_OVERLAP;
    }
    else
    {
      return VTK_PARTIAL_OVERLAP;
    }
  }
  else if (this->InBounds(s[1], S[0], S[1]))
  {
    overlap[0] = S[0];
    overlap[1] = s[1];
    if (overlap[0] == overlap[1])
    {
      return VTK_NODE_OVERLAP;
    }
    else
    {
      return VTK_PARTIAL_OVERLAP;
    }
  }
  return VTK_NO_OVERLAP;
}

//------------------------------------------------------------------------------
int vtkStructuredGridConnectivity::PartialOverlap(
  int A[2], int CardinalityOfA, int B[2], int CardinalityOfB, int overlap[2])
{
  if (CardinalityOfA > CardinalityOfB)
  {
    return (this->DoPartialOverlap(B, A, overlap));
  }
  else if (CardinalityOfB > CardinalityOfA)
  {
    return (this->DoPartialOverlap(A, B, overlap));
  }
  else
  {
    return (this->DoPartialOverlap(A, B, overlap));
  }

  // Code should not reach here!
  //  assert( "Hmm...code should not reach here!" && false );
  //  return VTK_NO_OVERLAP;
}

//------------------------------------------------------------------------------
int vtkStructuredGridConnectivity::IntervalOverlap(int A[2], int B[2], int overlap[2])
{
  // STEP 0: Check if we must check for a partial overlap
  int CardinalityOfA = this->Cardinality(A);
  int CardinalityOfB = this->Cardinality(B);
  return (this->PartialOverlap(A, CardinalityOfA, B, CardinalityOfB, overlap));
}

//------------------------------------------------------------------------------
void vtkStructuredGridConnectivity::GetIJKBlockOrientation(
  int i, int j, int k, int ext[6], int orientation[3])
{
  orientation[0] = orientation[1] = orientation[2] = BlockFace::NOT_ON_BLOCK_FACE;

  switch (this->DataDescription)
  {
    case vtkStructuredData::VTK_STRUCTURED_X_LINE:
      orientation[0] = this->Get1DOrientation(
        i, ext[0], ext[1], BlockFace::LEFT, BlockFace::RIGHT, BlockFace::NOT_ON_BLOCK_FACE);
      break;
    case vtkStructuredData::VTK_STRUCTURED_Y_LINE:
      orientation[1] = this->Get1DOrientation(
        j, ext[2], ext[3], BlockFace::BOTTOM, BlockFace::TOP, BlockFace::NOT_ON_BLOCK_FACE);
      break;
    case vtkStructuredData::VTK_STRUCTURED_Z_LINE:
      orientation[2] = this->Get1DOrientation(
        k, ext[4], ext[5], BlockFace::BACK, BlockFace::FRONT, BlockFace::NOT_ON_BLOCK_FACE);
      break;
    case vtkStructuredData::VTK_STRUCTURED_XY_PLANE:
      orientation[0] = this->Get1DOrientation(
        i, ext[0], ext[1], BlockFace::LEFT, BlockFace::RIGHT, BlockFace::NOT_ON_BLOCK_FACE);
      orientation[1] = this->Get1DOrientation(
        j, ext[2], ext[3], BlockFace::BOTTOM, BlockFace::TOP, BlockFace::NOT_ON_BLOCK_FACE);
      break;
    case vtkStructuredData::VTK_STRUCTURED_YZ_PLANE:
      orientation[1] = this->Get1DOrientation(
        j, ext[2], ext[3], BlockFace::BOTTOM, BlockFace::TOP, BlockFace::NOT_ON_BLOCK_FACE);
      orientation[2] = this->Get1DOrientation(
        k, ext[4], ext[5], BlockFace::BACK, BlockFace::FRONT, BlockFace::NOT_ON_BLOCK_FACE);
      break;
    case vtkStructuredData::VTK_STRUCTURED_XZ_PLANE:
      orientation[0] = this->Get1DOrientation(
        i, ext[0], ext[1], BlockFace::LEFT, BlockFace::RIGHT, BlockFace::NOT_ON_BLOCK_FACE);
      orientation[2] = this->Get1DOrientation(
        k, ext[4], ext[5], BlockFace::BACK, BlockFace::FRONT, BlockFace::NOT_ON_BLOCK_FACE);
      break;
    case vtkStructuredData::VTK_STRUCTURED_XYZ_GRID:
      orientation[0] = this->Get1DOrientation(
        i, ext[0], ext[1], BlockFace::LEFT, BlockFace::RIGHT, BlockFace::NOT_ON_BLOCK_FACE);
      orientation[1] = this->Get1DOrientation(
        j, ext[2], ext[3], BlockFace::BOTTOM, BlockFace::TOP, BlockFace::NOT_ON_BLOCK_FACE);
      orientation[2] = this->Get1DOrientation(
        k, ext[4], ext[5], BlockFace::BACK, BlockFace::FRONT, BlockFace::NOT_ON_BLOCK_FACE);
      break;
    default:
      std::cout << "Data description is: " << this->DataDescription << "\n";
      std::cout.flush();
      assert("pre: Undefined data-description!" && false);
  }
}

//------------------------------------------------------------------------------
void vtkStructuredGridConnectivity::CreateGhostedExtent(int gridID, int N)
{
  assert("pre: gridID is out-of-bounds!" && (gridID >= 0) &&
    (gridID < static_cast<int>(this->NumberOfGrids)));
  assert("pre: ghosted-extents vector has not been allocated" &&
    (this->NumberOfGrids == this->GhostedExtents.size() / 6));
  assert(
    "pre: Number of ghost-layers requested should not be 0" && (this->NumberOfGhostLayers > 0));

  int ext[6];
  this->GetGridExtent(gridID, ext);
  this->SetGhostedGridExtent(gridID, ext);

  int* ghostedExtent = &this->GhostedExtents[gridID * 6];

  switch (this->DataDescription)
  {
    case vtkStructuredData::VTK_STRUCTURED_X_LINE:
      this->GetGhostedExtent(ghostedExtent, ext, 0, 1, N);
      break;
    case vtkStructuredData::VTK_STRUCTURED_Y_LINE:
      this->GetGhostedExtent(ghostedExtent, ext, 2, 3, N);
      break;
    case vtkStructuredData::VTK_STRUCTURED_Z_LINE:
      this->GetGhostedExtent(ghostedExtent, ext, 4, 5, N);
      break;
    case vtkStructuredData::VTK_STRUCTURED_XY_PLANE:
      this->GetGhostedExtent(ghostedExtent, ext, 0, 1, N);
      this->GetGhostedExtent(ghostedExtent, ext, 2, 3, N);
      break;
    case vtkStructuredData::VTK_STRUCTURED_YZ_PLANE:
      this->GetGhostedExtent(ghostedExtent, ext, 2, 3, N);
      this->GetGhostedExtent(ghostedExtent, ext, 4, 5, N);
      break;
    case vtkStructuredData::VTK_STRUCTURED_XZ_PLANE:
      this->GetGhostedExtent(ghostedExtent, ext, 0, 1, N);
      this->GetGhostedExtent(ghostedExtent, ext, 4, 5, N);
      break;
    case vtkStructuredData::VTK_STRUCTURED_XYZ_GRID:
      this->GetGhostedExtent(ghostedExtent, ext, 0, 1, N);
      this->GetGhostedExtent(ghostedExtent, ext, 2, 3, N);
      this->GetGhostedExtent(ghostedExtent, ext, 4, 5, N);
      break;
    default:
      std::cout << "Data description is: " << this->DataDescription << "\n";
      std::cout.flush();
      assert("pre: Undefined data-description!" && false);
  } // END switch
}

//------------------------------------------------------------------------------
void vtkStructuredGridConnectivity::CreateGhostedMaskArrays(int gridID)
{
  // Sanity check
  assert("pre: gridID is out-of-bounds!" && (gridID >= 0) &&
    (gridID < static_cast<int>(this->NumberOfGrids)));
  assert("pre: GhostedPointGhostArray has not been allocated" &&
    (this->NumberOfGrids == this->GhostedPointGhostArray.size()));
  assert("pre: GhostedCellGhostArray has not been allocated" &&
    (this->NumberOfGrids == this->GhostedCellGhostArray.size()));

  // STEP 0: Initialize the ghosted node and cell arrays
  if (this->GhostedPointGhostArray[gridID] == nullptr)
  {
    this->GhostedPointGhostArray[gridID] = vtkUnsignedCharArray::New();
  }
  else
  {
    this->GhostedPointGhostArray[gridID]->Reset();
  }

  if (this->GhostedCellGhostArray[gridID] == nullptr)
  {
    this->GhostedCellGhostArray[gridID] = vtkUnsignedCharArray::New();
  }
  else
  {
    this->GhostedCellGhostArray[gridID]->Reset();
  }

  // STEP 1: Get the ghosted extent
  int ghostedExtent[6];
  this->GetGhostedGridExtent(gridID, ghostedExtent);

  // STEP 2: Get the grid extent
  int gridExtent[6];
  this->GetGridExtent(gridID, gridExtent);

  int numNodes = vtkStructuredData::GetNumberOfPoints(ghostedExtent, this->DataDescription);

  int numCells = vtkStructuredData::GetNumberOfCells(ghostedExtent, this->DataDescription);

  // STEP 3: Allocated the ghosted node and cell arrays and initialize them
  this->GhostedPointGhostArray[gridID]->Allocate(numNodes);
  this->GhostedCellGhostArray[gridID]->Allocate(numCells);

  // Initialize the arrays
  unsigned char* pnodes = this->GhostedPointGhostArray[gridID]->WritePointer(0, numNodes);
  memset(pnodes, 0, numNodes);
  unsigned char* pcells = this->GhostedCellGhostArray[gridID]->WritePointer(0, numCells);
  memset(pcells, 0, numCells);

  // STEP 4: Loop through the ghosted extent and mark the nodes in the ghosted
  // extent accordingly. If the node exists in the grown extent
  int ijk[3];
  unsigned char p = 0;
  for (int i = ghostedExtent[0]; i <= ghostedExtent[1]; ++i)
  {
    for (int j = ghostedExtent[2]; j <= ghostedExtent[3]; ++j)
    {
      for (int k = ghostedExtent[4]; k <= ghostedExtent[5]; ++k)
      {
        ijk[0] = i;
        ijk[1] = j;
        ijk[2] = k;

        vtkIdType idx =
          vtkStructuredData::ComputePointIdForExtent(ghostedExtent, ijk, this->DataDescription);

        if (this->IsNodeWithinExtent(i, j, k, gridExtent))
        {
          // Get index w.r.t. the register extent
          vtkIdType srcidx =
            vtkStructuredData::ComputePointIdForExtent(gridExtent, ijk, this->DataDescription);
          if (this->GridPointGhostArrays[gridID])
          {
            p = this->GridPointGhostArrays[gridID]->GetValue(srcidx);
            this->GhostedPointGhostArray[gridID]->SetValue(idx, p);
          }
        }
        else
        {
          p = 0;

          if (this->IsNodeOnBoundary(i, j, k))
          {
            // We might use BOUNDARY in the future
          }
          p |= vtkDataSetAttributes::DUPLICATEPOINT;
          this->GhostedPointGhostArray[gridID]->SetValue(idx, p);
        }
      } // END for all k
    }   // END for all j
  }     // END for all i

  // STEP 5: Fill the cells ghost arrays for the ghosted grid
  int dim = vtkStructuredData::GetDataDimension(this->DataDescription);
  assert("pre: data dimensions must be 1, 2 or 3" && (dim >= 1) && (dim <= 3));

  int dims[3];
  vtkStructuredData::GetDimensionsFromExtent(ghostedExtent, dims);

  int numNodesPerCell = this->GetNumberOfNodesPerCell(dim);

  int CellExtent[6];
  vtkStructuredData::GetCellExtentFromPointExtent(ghostedExtent, CellExtent);

  this->FillCellsGhostArray(this->DataDescription, numNodesPerCell, dims, CellExtent,
    this->GhostedPointGhostArray[gridID], this->GhostedCellGhostArray[gridID]);
}

//------------------------------------------------------------------------------
void vtkStructuredGridConnectivity::AllocatePointData(vtkPointData* RPD, int N, vtkPointData* PD)
{
  assert("pre: Reference point data is nullptr" && (RPD != nullptr));
  assert("pre: point data is nullptr" && (PD != nullptr));
  assert("pre: N > 0" && (N > 0));

  for (int array = 0; array < RPD->GetNumberOfArrays(); ++array)
  {
    int dataType = RPD->GetArray(array)->GetDataType();
    vtkDataArray* dataArray = vtkDataArray::CreateDataArray(dataType);
    assert("Cannot create data array" && (dataArray != nullptr));

    dataArray->SetName(RPD->GetArray(array)->GetName());
    dataArray->SetNumberOfComponents(RPD->GetArray(array)->GetNumberOfComponents());
    dataArray->SetNumberOfTuples(N);

    PD->AddArray(dataArray);
    dataArray->Delete();
  } // END for all arrays
}

//------------------------------------------------------------------------------
void vtkStructuredGridConnectivity::AllocateCellData(vtkCellData* RCD, int N, vtkCellData* CD)
{
  assert("pre: Reference cell data is nullptr" && (RCD != nullptr));
  assert("pre: cell data is nullptr" && (CD != nullptr));
  assert("pre: N > 0" && (N > 0));

  for (int array = 0; array < RCD->GetNumberOfArrays(); ++array)
  {
    int dataType = RCD->GetArray(array)->GetDataType();
    vtkDataArray* dataArray = vtkDataArray::CreateDataArray(dataType);
    assert("Cannot create data array" && (dataArray != nullptr));

    dataArray->SetName(RCD->GetArray(array)->GetName());
    dataArray->SetNumberOfComponents(RCD->GetArray(array)->GetNumberOfComponents());
    dataArray->SetNumberOfTuples(N);

    CD->AddArray(dataArray);
    dataArray->Delete();
  } // END for all arrays
}

//------------------------------------------------------------------------------
void vtkStructuredGridConnectivity::InitializeGhostData(int gridID)
{
  // Sanity check
  assert("pre: gridID is out-of-bounds!" && (gridID >= 0) &&
    (gridID < static_cast<int>(this->NumberOfGrids)));
  assert("pre: GhostedPointData vector has not been properly allocated!" &&
    (this->NumberOfGrids == this->GhostedGridPointData.size()));
  assert("pre: GhostedCellData vector has not been properly allocated!" &&
    (this->NumberOfGrids == this->GhostedGridCellData.size()));
  assert("pre: Grid has no registered point data!" && (this->GridPointData[gridID] != nullptr));
  assert("pre: Grid has no registered cell data!" && (this->GridCellData[gridID] != nullptr));

  // STEP 0: Get the ghosted grid extent
  int GhostedGridExtent[6];
  this->GetGhostedGridExtent(gridID, GhostedGridExtent);

  // STEP 1: Get the number of nodes/cells in the ghosted extent
  int numNodes = vtkStructuredData::GetNumberOfPoints(GhostedGridExtent, this->DataDescription);
  int numCells = vtkStructuredData::GetNumberOfCells(GhostedGridExtent, this->DataDescription);

  // STEP 2: Allocate coordinates if the grid
  if (this->GridPoints[gridID] != nullptr)
  {

    if (this->GhostedGridPoints[gridID] != nullptr)
    {
      this->GhostedGridPoints[gridID]->Delete();
    }

    this->GhostedGridPoints[gridID] = vtkPoints::New();
    this->GhostedGridPoints[gridID]->SetDataTypeToDouble();
    this->GhostedGridPoints[gridID]->SetNumberOfPoints(numNodes);
  }

  // STEP 3: Allocate point & cell data
  this->GhostedGridPointData[gridID] = vtkPointData::New();
  this->GhostedGridCellData[gridID] = vtkCellData::New();

  this->AllocatePointData(
    this->GridPointData[gridID], numNodes, this->GhostedGridPointData[gridID]);
  this->AllocateCellData(this->GridCellData[gridID], numCells, this->GhostedGridCellData[gridID]);
}

//------------------------------------------------------------------------------
void vtkStructuredGridConnectivity::CopyCoordinates(
  vtkPoints* source, vtkIdType sourceIdx, vtkPoints* target, vtkIdType targetIdx)
{
  assert("pre: source points is nullptr" && (source != nullptr));
  assert("pre: target points is nullptr" && (target != nullptr));
  assert("pre: source index is out-of-bounds!" && (sourceIdx >= 0) &&
    (sourceIdx < source->GetNumberOfPoints()));
  assert("pre: target index is out-of-bounds!" && (targetIdx >= 0) &&
    (targetIdx < target->GetNumberOfPoints()));
  target->SetPoint(targetIdx, source->GetPoint(sourceIdx));
}

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

  int arrayIdx = 0;
  for (; arrayIdx < source->GetNumberOfArrays(); ++arrayIdx)
  {
    // Get source array
    vtkDataArray* sourceArray = source->GetArray(arrayIdx);
    assert("ERROR: encountered nullptr source array" && (sourceArray != nullptr));

    // Get target array
    vtkDataArray* targetArray = target->GetArray(arrayIdx);
    assert("ERROR: encountered nullptr target array" && (targetArray != nullptr));

    // Sanity checks
    assert("ERROR: target/source array name mismatch!" &&
      (strcmp(sourceArray->GetName(), targetArray->GetName()) == 0));
    assert("ERROR: target/source array num components mismatch!" &&
      (sourceArray->GetNumberOfComponents() == targetArray->GetNumberOfComponents()));
    assert("ERROR: sourceIdx out-of-bounds!" && (sourceIdx >= 0) &&
      (sourceIdx < sourceArray->GetNumberOfTuples()));
    assert("ERROR: targetIdx out-of-bounds!" && (targetIdx >= 0) &&
      (targetIdx < targetArray->GetNumberOfTuples()));

    // Copy the tuple
    targetArray->SetTuple(targetIdx, sourceIdx, sourceArray);
  } // END for all arrays
}

//------------------------------------------------------------------------------
void vtkStructuredGridConnectivity::TransferRegisteredDataToGhostedData(int gridID)
{
  // Sanity check
  assert("pre: gridID is out-of-bounds!" && (gridID >= 0) &&
    (gridID < static_cast<int>(this->NumberOfGrids)));

  // STEP 0: Get the registered grid extent
  int GridExtent[6];
  this->GetGridExtent(gridID, GridExtent);

  // STEP 1: Get the registered grid cell extent
  int GridCellExtent[6];
  vtkStructuredData::GetCellExtentFromPointExtent(
    GridExtent, GridCellExtent, this->DataDescription);

  // STEP 2: Get the ghosted grid extent
  int GhostedGridExtent[6];
  this->GetGhostedGridExtent(gridID, GhostedGridExtent);

  // STEP 3: Get the ghosted grid cell extent
  int GhostedGridCellExtent[6];
  vtkStructuredData::GetCellExtentFromPointExtent(
    GhostedGridExtent, GhostedGridCellExtent, this->DataDescription);

  // STEP 2: Loop over the registered grid extent
  int ijk[3];
  for (int i = GridExtent[0]; i <= GridExtent[1]; ++i)
  {
    for (int j = GridExtent[2]; j <= GridExtent[3]; ++j)
    {
      for (int k = GridExtent[4]; k <= GridExtent[5]; ++k)
      {

        ijk[0] = i;
        ijk[1] = j;
        ijk[2] = k;

        // Compute the source index to the registered data
        vtkIdType sourceIdx =
          vtkStructuredData::ComputePointIdForExtent(GridExtent, ijk, this->DataDescription);

        // Compute the target index to the ghosted data
        vtkIdType targetIdx =
          vtkStructuredData::ComputePointIdForExtent(GhostedGridExtent, ijk, this->DataDescription);

        if (this->GridPoints[gridID] != nullptr)
        {
          this->CopyCoordinates(
            this->GridPoints[gridID], sourceIdx, this->GhostedGridPoints[gridID], targetIdx);
        } // END if grid points is not nullptr

        // Transfer node data from the registered grid to the ghosted grid
        this->CopyFieldData(
          this->GridPointData[gridID], sourceIdx, this->GhostedGridPointData[gridID], targetIdx);

        // If the node is within the cell extent, copy the cell datta
        if (this->IsNodeWithinExtent(i, j, k, GridCellExtent))
        {
          // Compute the source cell idx. Note, since we are passing to
          // ComputePointIdForExtent a cell extent, this is a cell id, not
          // a point id.
          vtkIdType sourceCellIdx =
            vtkStructuredData::ComputePointIdForExtent(GridCellExtent, ijk, this->DataDescription);

          // Compute the target cell idx. Note, since we are passing to
          // ComputePointIdForExtent a cell extent, this is a cell id, not
          // a point id.
          vtkIdType targetCellIdx = vtkStructuredData::ComputePointIdForExtent(
            GhostedGridCellExtent, ijk, this->DataDescription);

          // Transfer cell data from the registered grid to the ghosted grid
          this->CopyFieldData(this->GridCellData[gridID], sourceCellIdx,
            this->GhostedGridCellData[gridID], targetCellIdx);
        } // END if node is within cell extent

      } // END for all k
    }   // END for all j
  }     // END for all i
}

//------------------------------------------------------------------------------
void vtkStructuredGridConnectivity::TransferGhostDataFromNeighbors(int gridID)
{
  // Sanity check
  assert("pre: gridID is out-of-bounds!" && (gridID >= 0) &&
    (gridID < static_cast<int>(this->NumberOfGrids)));
  assert(
    "pre: Neighbors is not properly allocated" && (this->NumberOfGrids == this->Neighbors.size()));

  int NumNeis = static_cast<int>(this->Neighbors[gridID].size());
  for (int nei = 0; nei < NumNeis; ++nei)
  {
    this->TransferLocalNeighborData(gridID, this->Neighbors[gridID][nei]);
  } // END for all neighbors
}

//------------------------------------------------------------------------------
void vtkStructuredGridConnectivity::TransferLocalNeighborData(
  int gridID, const vtkStructuredNeighbor& Neighbor)
{
  // Sanity check
  assert("pre: gridID is out-of-bounds!" && (gridID >= 0) &&
    (gridID < static_cast<int>(this->NumberOfGrids)));
  assert("pre: Neighbor gridID is out-of-bounds!" && (Neighbor.NeighborID >= 0) &&
    (Neighbor.NeighborID < static_cast<int>(this->NumberOfGrids)));

  // STEP 0: Get ghosted grid (node) extent and corresponding cell extent
  int GhostedGridExtent[6];
  this->GetGhostedGridExtent(gridID, GhostedGridExtent);
  int GhostedGridCellExtent[6];
  vtkStructuredData::GetCellExtentFromPointExtent(GhostedGridExtent, GhostedGridCellExtent);

  // STEP 1: Get the neighbor (node) extent and corresponding cell extent
  int NeighborExtent[6];
  this->GetGridExtent(Neighbor.NeighborID, NeighborExtent);
  int NeighborCellExtent[6];
  vtkStructuredData::GetCellExtentFromPointExtent(NeighborExtent, NeighborCellExtent);

  int RcvCellExtent[6];
  vtkStructuredData::GetCellExtentFromPointExtent(
    const_cast<int*>(Neighbor.RcvExtent), RcvCellExtent);

  // STEP 3: Transfer the RcvExtent to the grid from the Neighbor
  int ijk[3];
  for (int i = Neighbor.RcvExtent[0]; i <= Neighbor.RcvExtent[1]; ++i)
  {
    for (int j = Neighbor.RcvExtent[2]; j <= Neighbor.RcvExtent[3]; ++j)
    {
      for (int k = Neighbor.RcvExtent[4]; k <= Neighbor.RcvExtent[5]; ++k)
      {
        // Sanity check!
        assert("pre: RcvExtent is outside the GhostExtent!" &&
          this->IsNodeWithinExtent(i, j, k, GhostedGridExtent));
        assert("pre: RcvExtent is outside the NeighborExtent" &&
          this->IsNodeWithinExtent(i, j, k, NeighborExtent));

        ijk[0] = i;
        ijk[1] = j;
        ijk[2] = k;

        // Compute the source index to the registered neighbor data
        vtkIdType srcIdx =
          vtkStructuredData::ComputePointIdForExtent(NeighborExtent, ijk, this->DataDescription);

        // Compute the target index into the ghosted data
        vtkIdType targetIdx =
          vtkStructuredData::ComputePointIdForExtent(GhostedGridExtent, ijk, this->DataDescription);

        if (this->GridPoints[Neighbor.NeighborID] != nullptr)
        {
          this->CopyCoordinates(this->GridPoints[Neighbor.NeighborID], srcIdx,
            this->GhostedGridPoints[gridID], targetIdx);
        } // END if this

        // Transfer node data from the registered grid to the ghosted grid
        this->CopyFieldData(this->GridPointData[Neighbor.NeighborID], srcIdx,
          this->GhostedGridPointData[gridID], targetIdx);

        if (this->IsNodeWithinExtent(i, j, k, RcvCellExtent))
        {
          // Compute the source cell idx. Note, since we are passing to
          // ComputePointIdForExtent a cell extent, this is a cell id, not
          // a point id.
          vtkIdType sourceCellIdx = vtkStructuredData::ComputePointIdForExtent(
            NeighborCellExtent, ijk, this->DataDescription);

          // Compute the target cell idx. Note, since we are passing to
          // ComputePointIdForExtent a cell extent, this is a cell id, not
          // a point id.
          vtkIdType targetCellIdx = vtkStructuredData::ComputePointIdForExtent(
            GhostedGridCellExtent, ijk, this->DataDescription);

          // Transfer cell data from the registered grid to the ghosted grid
          this->CopyFieldData(this->GridCellData[Neighbor.NeighborID], sourceCellIdx,
            this->GhostedGridCellData[gridID], targetCellIdx);
        } // END if node is within cell extent

      } // END for all k
    }   // END for all j
  }     // END for all i
}

//------------------------------------------------------------------------------
void vtkStructuredGridConnectivity::ComputeNeighborSendAndRcvExtent(int gridID, int N)
{
  // Sanity check
  assert("pre: gridID is out-of-bounds!" && (gridID >= 0) &&
    (gridID < static_cast<int>(this->NumberOfGrids)));
  assert(
    "pre: Neighbors is not properly allocated" && (this->NumberOfGrids == this->Neighbors.size()));

  int myRealGridExtent[6];
  this->GetGridExtent(gridID, myRealGridExtent);

  int myGhostedGridExtent[6];
  this->GetGhostedGridExtent(gridID, myGhostedGridExtent);

  int NumNeis = static_cast<int>(this->Neighbors[gridID].size());
  for (int nei = 0; nei < NumNeis; ++nei)
  {
    int neiRealExtent[6];
    this->GetGridExtent(this->Neighbors[gridID][nei].NeighborID, neiRealExtent);

    this->Neighbors[gridID][nei].ComputeSendAndReceiveExtent(
      myRealGridExtent, myGhostedGridExtent, neiRealExtent, this->WholeExtent, N);
  }
}
//------------------------------------------------------------------------------
void vtkStructuredGridConnectivity::CreateGhostLayers(int N)
{
  if (N == 0)
  {
    vtkWarningMacro("N=0 ghost layers requested! No ghost layers will be created");
    return;
  }

  this->NumberOfGhostLayers += N;
  this->AllocateInternalDataStructures();
  this->GhostedExtents.resize(this->NumberOfGrids * 6, -1);

  for (unsigned int i = 0; i < this->NumberOfGrids; ++i)
  {
    this->CreateGhostedExtent(i, N);
    this->CreateGhostedMaskArrays(i);
    this->ComputeNeighborSendAndRcvExtent(i, N);
    this->InitializeGhostData(i);
    this->TransferRegisteredDataToGhostedData(i);
    this->TransferGhostDataFromNeighbors(i);
  } // END for all grids
}
VTK_ABI_NAMESPACE_END