vtkStructuredGridConnectivity.h

Go to the documentation of this file.

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

 Program:   Visualization Toolkit
 Module:    vtkStructuredGridConnectivity.h

 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
 All rights reserved.
 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

 This software is distributed WITHOUT ANY WARRANTY; without even
 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
 PURPOSE.  See the above copyright notice for more information.

 =========================================================================*/
#ifndef vtkStructuredGridConnectivity_h
#define vtkStructuredGridConnectivity_h

#define VTK_NO_OVERLAP      0
#define VTK_NODE_OVERLAP    1
#define VTK_EDGE_OVERLAP    2
#define VTK_PARTIAL_OVERLAP 3

// VTK include directives
#include "vtkFiltersGeometryModule.h" // For export macro
#include "vtkAbstractGridConnectivity.h"
#include "vtkStructuredNeighbor.h" // For Structured Neighbor object definition
#include "vtkStructuredData.h" // For data description definitions

// C++ include directives
#include <iostream> // For cout
#include <vector>   // For STL vector
#include <map>      // For STL map
#include <utility>  // For STL pair and overloaded relational operators
#include <cassert>  // For assert()

// Forward Declarations
class vtkIdList;
class vtkUnsignedCharArray;
class vtkPointData;
class vtkCellData;
class vtkPoints;

class VTKFILTERSGEOMETRY_EXPORT vtkStructuredGridConnectivity :
  public vtkAbstractGridConnectivity
{
public:
  static vtkStructuredGridConnectivity* New();
  vtkTypeMacro( vtkStructuredGridConnectivity, vtkAbstractGridConnectivity );
  void PrintSelf(ostream& os, vtkIndent  indent ) override;


  vtkSetVector6Macro(WholeExtent,int);
  vtkGetVector6Macro(WholeExtent,int);


  vtkGetMacro(DataDimension,int);

  void SetNumberOfGrids( const unsigned int N ) override;

  virtual void RegisterGrid( const int gridID, int extents[6],
      vtkUnsignedCharArray* nodesGhostArray,
      vtkUnsignedCharArray* cellGhostArray,
      vtkPointData* pointData,
      vtkCellData* cellData,
      vtkPoints* gridNodes );

  void GetGridExtent( const int gridID, int extent[6] );

  void SetGhostedGridExtent( const int gridID, int ext[6] );

  void GetGhostedGridExtent( const int gridID, int ext[6] );

  void ComputeNeighbors() override;

  int GetNumberOfNeighbors( const int gridID )
    { return( static_cast<int>(this->Neighbors[ gridID ].size() )); };

  vtkStructuredNeighbor GetGridNeighbor(const int gridID, const int nei);

  vtkIdList* GetNeighbors( const int gridID, int *extents );

  void FillGhostArrays(
     const int gridID,
     vtkUnsignedCharArray *nodesArray,
     vtkUnsignedCharArray *cellsArray ) override;

  void CreateGhostLayers( const int N=1 ) override;

protected:
  vtkStructuredGridConnectivity();
  ~vtkStructuredGridConnectivity() override;

  bool InBounds( const int idx, const int Lo, const int Hi )
  { return( (idx>=Lo) && (idx<=Hi) ); };

  bool StrictlyInsideBounds( const int idx, const int Lo, const int Hi )
  { return( (idx > Lo) && (idx < Hi) ); };

  bool IsSubset( int A[2], int B[2] )
  { return( this->InBounds(A[0], B[0], B[1]) &&
            this->InBounds(A[1], B[0], B[1]) ); };

  int Cardinality( int S[2] ) { return( S[1]-S[0]+1 ); };


  int GetNumberOfNodesPerCell( const int dim )
  {
      int numNodes = 0;
      switch( dim )
      {
        case 1:
         numNodes = 2; // line cell
         break;
        case 2:
         numNodes = 4; // quad cell
         break;
        case 3:
         numNodes = 8; // hex cell
         break;
        default:
         assert( "ERROR: code should not reach here!" && false );
      } // END switch
      return( numNodes );
  }

  void FillNodesGhostArray(
      const int gridID, const int dataDescription,
      int GridExtent[6], int RealExtent[6], vtkUnsignedCharArray *nodeArray );

  void FillCellsGhostArray(
      const int dataDescription, const int numNodesPerCell,
      int dims[3], int CellExtent[6], vtkUnsignedCharArray *nodesArray,
      vtkUnsignedCharArray *cellsArray );

  void SearchNeighbors(
      const int gridID,
      const int i, const int j, const int k,
      vtkIdList *neiList );

  void MarkNodeProperty(
      const int gridID,
      const int i, const int j, const int k,
      int ext[6], int RealExtent[6], unsigned char &pfield );

  void MarkCellProperty(
      unsigned char &pfield,
      unsigned char *nodeGhostFields, const int numNodes );

  void GetRealExtent( const int gridID, int GridExtent[6],int RealExtent[6] );

  bool IsGhostNode(
      int GridExtent[6], int RealExtent[6],
      const int i, const int j, const int k );

  bool IsNodeOnBoundaryOfExtent(
      const int i, const int j, const int k, int ext[6] );

  bool IsNodeOnSharedBoundary(
      const int gridID, int RealExtent[6],
      const int i, const int j, const int k );

  bool IsNodeOnBoundary( const int i, const int j, const int k );

  bool IsNodeInterior(
      const int i, const int j, const int k,
      int GridExtent[6] );

  bool IsNodeWithinExtent(
      const int i, const int j, const int k,
      int GridExtent[6] )
  {
    bool status = false;

    switch( this->DataDescription )
    {
      case VTK_X_LINE:
        if( (GridExtent[0] <= i) && (i <= GridExtent[1]) )
        {
          status = true;
        }
        break;
      case VTK_Y_LINE:
        if( (GridExtent[2] <= j) && (j <= GridExtent[3] ) )
        {
          status = true;
        }
        break;
      case VTK_Z_LINE:
        if( (GridExtent[4] <= k) && (k <= GridExtent[5] ) )
        {
          status = true;
        }
        break;
      case VTK_XY_PLANE:
        if( (GridExtent[0] <= i) && (i <= GridExtent[1]) &&
            (GridExtent[2] <= j) && (j <= GridExtent[3])  )
        {
          status = true;
        }
        break;
      case VTK_YZ_PLANE:
        if( (GridExtent[2] <= j) && (j <= GridExtent[3] ) &&
            (GridExtent[4] <= k) && (k <= GridExtent[5] ) )
        {
          status = true;
        }
        break;
      case VTK_XZ_PLANE:
        if( (GridExtent[0] <= i) && (i <= GridExtent[1] ) &&
            (GridExtent[4] <= k) && (k <= GridExtent[5] ) )
        {
          status = true;
        }
        break;
      case VTK_XYZ_GRID:
        if( (GridExtent[0] <= i) && (i <= GridExtent[1]) &&
            (GridExtent[2] <= j) && (j <= GridExtent[3]) &&
            (GridExtent[4] <= k) && (k <= GridExtent[5]) )
        {
          status = true;
        }
        break;
      default:
        std::cout << "Data description is: " << this->DataDescription << "\n";
        std::cout.flush();
        assert( "pre: Undefined data-description!" && false );
    } // END switch

    return( status );
  }

  void SetNeighbors(
      const int i, const int j,
      int i2jOrientation[3], int j2iOrientation[3],
      int overlapExtent[6] );

  void DetermineNeighborOrientation(
      const int  idx, int A[2], int B[2], int overlap[2], int orient[3] );

  void DetectNeighbors(
      const int i, const int j, int ex1[6], int ex2[6],
      int orientation[3], int ndim);

  int IntervalOverlap( int A[2], int B[2], int overlap[2] );

  int DoPartialOverlap( int s[2], int S[2], int overlap[2] );

  int PartialOverlap(
      int A[2], const int CofA,
      int B[2], const int CofB,
      int overlap[2] );

  void EstablishNeighbors( const int i, const int j );

  void AcquireDataDescription();

  bool HasBlockConnection( const int gridID, const int blockDirection );

  void RemoveBlockConnection( const int gridID, const int blockDirection );

  void AddBlockConnection( const int gridID, const int blockDirection );

  void ClearBlockConnections( const int gridID );

  int GetNumberOfConnectingBlockFaces( const int gridID );

  void SetBlockTopology( const int gridID );

  void GetIJKBlockOrientation(
      const int i, const int j, const int k, int ext[6], int orientation[3] );

  int Get1DOrientation(
      const int idx, const int ExtentLo, const int ExtentHi,
      const int OnLo, const int OnHi, const int NotOnBoundary );

  void CreateGhostedExtent( const int gridID, const int N );

  void GetGhostedExtent(
      int *ghostedExtent, int GridExtent[6],
      const int minIdx, const int maxIdx, const int N);

  void CreateGhostedMaskArrays(const int gridID);

  void InitializeGhostData( const int gridID );

  void AllocatePointData( vtkPointData *RPD, const int N, vtkPointData *PD );

  void AllocateCellData( vtkCellData *RCD, const int N, vtkCellData *CD );

  void TransferRegisteredDataToGhostedData( const int gridID );

  void ComputeNeighborSendAndRcvExtent( const int gridID, const int N );

  virtual void TransferGhostDataFromNeighbors( const int gridID );

  void TransferLocalNeighborData(
      const int gridID, const vtkStructuredNeighbor& Neighor);

  void CopyCoordinates(
      vtkPoints *source, vtkIdType sourceIdx,
      vtkPoints *target, vtkIdType targetIdx );

  void CopyFieldData(
      vtkFieldData *source, vtkIdType sourceIdx,
      vtkFieldData *target, vtkIdType targetIdx );

  int GetNeighborIndex( const int gridIdx, const int NeighborGridIdx );

  void PrintExtent( int extent[6] );

  int DataDimension;
  int DataDescription;
  int WholeExtent[6];

  std::vector< int > GridExtents;
  std::vector< int > GhostedExtents;
  std::vector< unsigned char  > BlockTopology;
  std::vector< std::vector<vtkStructuredNeighbor> > Neighbors;
  std::map< std::pair< int,int >, int > NeighborPair2NeighborListIndex;

private:
  vtkStructuredGridConnectivity( const vtkStructuredGridConnectivity& ) = delete;
  void operator=(const vtkStructuredGridConnectivity& ) = delete;
};

//=============================================================================
//  INLINE METHODS
//=============================================================================

//------------------------------------------------------------------------------
inline int vtkStructuredGridConnectivity::GetNeighborIndex(
    const int gridIdx, const int NeighborGridIdx )
{
  assert("pre: Grid index is out-of-bounds!" &&
         (gridIdx >= 0) &&
         (gridIdx < static_cast<int>(this->NumberOfGrids)));
  assert("pre: Neighbor grid index is out-of-bounds!" &&
         (NeighborGridIdx >= 0) &&
         (NeighborGridIdx < static_cast<int>(this->NumberOfGrids) ) );

  std::pair<int,int> gridPair = std::make_pair(gridIdx,NeighborGridIdx);
  assert("pre: Neighboring grid pair does not exist in hash!" &&
         (this->NeighborPair2NeighborListIndex.find(gridPair) !=
             this->NeighborPair2NeighborListIndex.end() ) );

  return(this->NeighborPair2NeighborListIndex[gridPair]);
}

//------------------------------------------------------------------------------
inline void vtkStructuredGridConnectivity::GetGhostedExtent(
    int *ghostedExtent, int GridExtent[6],
    const int minIdx, const int maxIdx, const int N )
{
  assert( "pre: Number of ghost layers must be N >= 1" && (N >= 1) );
  assert( "pre: ghosted extent pointer is nullptr" && ghostedExtent != nullptr);

  ghostedExtent[minIdx] = GridExtent[minIdx]-N;
  ghostedExtent[maxIdx] = GridExtent[maxIdx]+N;

  // Clamp the ghosted extent to be within the WholeExtent
  ghostedExtent[minIdx] =
   (ghostedExtent[minIdx] < this->WholeExtent[minIdx] )?
       this->WholeExtent[minIdx] : ghostedExtent[minIdx];
  ghostedExtent[maxIdx] =
   (ghostedExtent[maxIdx] > this->WholeExtent[maxIdx])?
       this->WholeExtent[maxIdx] : ghostedExtent[maxIdx];
}

//------------------------------------------------------------------------------
inline void vtkStructuredGridConnectivity::SetGhostedGridExtent(
    const int gridID, int ext[6] )
{
  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 ) );

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

//------------------------------------------------------------------------------
inline void vtkStructuredGridConnectivity::GetGridExtent(
    const int gridID, int ext[6])
{
  assert( "pre: gridID out-of-bounds!" &&
        (gridID >= 0  && gridID < static_cast<int>(this->NumberOfGrids)));
  for( int i=0; i < 6; ++i )
  {
    ext[i] = this->GridExtents[ gridID*6+i ];
  }
}

//------------------------------------------------------------------------------
inline void vtkStructuredGridConnectivity::GetGhostedGridExtent(
    const int gridID, int ext[6])
{
  assert( "pre: gridID out-of-bounds!" &&
        (gridID >= 0  && gridID < static_cast<int>(this->NumberOfGrids)));

  if( this->GhostedExtents.size() == 0 )
  {
    ext[0] = ext[2] = ext[4] = -1;
    ext[1] = ext[3] = ext[5] = 0;
    vtkErrorMacro( "No ghosted extents found for registered grid extends!!!" );
    return;
  }

  assert( "GhostedExtents are not aligned with registered grid extents" &&
        ( this->GhostedExtents.size() == this->GridExtents.size() ) );
  for( int i=0; i < 6; ++i )
  {
    ext[i] = this->GhostedExtents[ gridID*6+i ];
  }
}

//------------------------------------------------------------------------------
inline bool vtkStructuredGridConnectivity::IsNodeOnBoundaryOfExtent(
    const int i, const int j, const int k, int ext[6] )
{
  if( !this->IsNodeWithinExtent( i,j,k, ext) )
  {
    return false;
  }

  bool status = false;
  switch( this->DataDescription )
  {
    case VTK_X_LINE:
       if( i==ext[0] || i==ext[1] )
       {
         status = true;
       }
       break;
     case VTK_Y_LINE:
       if( j==ext[2] || j==ext[3] )
       {
         status = true;
       }
       break;
     case VTK_Z_LINE:
       if( k==ext[4] || k==ext[5] )
       {
         status = true;
       }
       break;
     case VTK_XY_PLANE:
       if( (i==ext[0] || i==ext[1]) ||
           (j==ext[2] || j==ext[3]) )
       {
         status = true;
       }
       break;
     case VTK_YZ_PLANE:
       if( (j==ext[2] || j==ext[3]) ||
           (k==ext[4] || k==ext[5]) )
       {
         status = true;
       }
       break;
     case VTK_XZ_PLANE:
       if( (i==ext[0] || i==ext[1]) ||
           (k==ext[4] || k==ext[5]) )
       {
         status = true;
       }
       break;
     case VTK_XYZ_GRID:
       if( (i==ext[0] || i==ext[1]) ||
           (j==ext[2] || j==ext[3]) ||
           (k==ext[4] || k==ext[5]) )
       {
         status = true;
       }
       break;
     default:
       std::cout << "Data description is: " << this->DataDescription << "\n";
       std::cout.flush();
       assert( "pre: Undefined data-description!" && false );
  } // END switch

  return( status );
}

//------------------------------------------------------------------------------
inline bool vtkStructuredGridConnectivity::IsNodeInterior(
    const int i, const int j, const int k,
    int GridExtent[6] )
{
  bool status = false;

  switch( this->DataDescription )
  {
    case VTK_X_LINE:
      if( (GridExtent[0] < i) && (i < GridExtent[1]) )
      {
        status = true;
      }
      break;
    case VTK_Y_LINE:
      if( (GridExtent[2] < j) && (j < GridExtent[3] ) )
      {
        status = true;
      }
      break;
    case VTK_Z_LINE:
      if( (GridExtent[4] < k) && (k < GridExtent[5] ) )
      {
        status = true;
      }
      break;
    case VTK_XY_PLANE:
      if( (GridExtent[0] < i) && (i < GridExtent[1]) &&
          (GridExtent[2] < j) && (j < GridExtent[3])  )
      {
        status = true;
      }
      break;
    case VTK_YZ_PLANE:
      if( (GridExtent[2] < j) && (j < GridExtent[3] ) &&
          (GridExtent[4] < k) && (k < GridExtent[5] ) )
      {
        status = true;
      }
      break;
    case VTK_XZ_PLANE:
      if( (GridExtent[0] < i) && (i < GridExtent[1] ) &&
          (GridExtent[4] < k) && (k < GridExtent[5] ) )
      {
        status = true;
      }
      break;
    case VTK_XYZ_GRID:
      if( (GridExtent[0] < i) && (i < GridExtent[1]) &&
          (GridExtent[2] < j) && (j < GridExtent[3]) &&
          (GridExtent[4] < k) && (k < GridExtent[5]) )
      {
        status = true;
      }
      break;
    default:
      std::cout << "Data description is: " << this->DataDescription << "\n";
      std::cout.flush();
      assert( "pre: Undefined data-description!" && false );
  } // END switch

  return( status );
}

//------------------------------------------------------------------------------
inline void vtkStructuredGridConnectivity::DetermineNeighborOrientation(
    const int idx, int A[2], int B[2], int overlap[2], int orient[3] )
{
  assert( "pre: idx is out-of-bounds" && (idx >= 0) && (idx < 3)  );

  // A. Non-overlapping cases
  if( overlap[0] == overlap[1] )
  {
    if( A[1] == B[0] )
    {
      orient[ idx ] = vtkStructuredNeighbor::HI;
    }
    else if( A[0] == B[1] )
    {
      orient[ idx ] = vtkStructuredNeighbor::LO;
    }
    else
    {
      orient[ idx ] = vtkStructuredNeighbor::UNDEFINED;
      assert( "ERROR: Code should not reach here!" && false );
    }
  } // END non-overlapping cases
  // B. Sub-set cases
  else if( this->IsSubset( A, B) )
  {
    if( (A[0] == B[0]) && (A[1] == B[1]) )
    {
      orient[ idx ] = vtkStructuredNeighbor::ONE_TO_ONE;
    }
    else if( this->StrictlyInsideBounds( A[0], B[0], B[1] ) &&
             this->StrictlyInsideBounds( A[1], B[0], B[1] ) )
    {
      orient[ idx ] = vtkStructuredNeighbor::SUBSET_BOTH;
    }
    else if( A[0] == B[0] )
    {
      orient[ idx ] = vtkStructuredNeighbor::SUBSET_HI;
    }
    else if( A[1] == B[1] )
    {
      orient[ idx ] = vtkStructuredNeighbor::SUBSET_LO;
    }
    else
    {
      orient[ idx ] = vtkStructuredNeighbor::UNDEFINED;
      assert( "ERROR: Code should not reach here!" && false );
    }
  }
  // C. Super-set cases
  else if( this->IsSubset( B, A ) )
  {
    orient[ idx ] = vtkStructuredNeighbor::SUPERSET;
  }
  // D. Partially-overlapping (non-subset) cases
  else if( !(this->IsSubset(A,B) || this->IsSubset(A,B)) )
  {
    if( this->InBounds( A[0], B[0], B[1] ) )
    {
      orient[ idx ] = vtkStructuredNeighbor::LO;
    }
    else if( this->InBounds( A[1], B[0], B[1] ) )
    {
      orient[ idx ] = vtkStructuredNeighbor::HI;
    }
    else
    {
      orient[ idx ] = vtkStructuredNeighbor::UNDEFINED;
      assert( "ERROR: Code should not reach here!" && false );
    }
  }
  else
  {
    orient[ idx ] = vtkStructuredNeighbor::UNDEFINED;
    assert( "ERROR: Code should not reach here!" && false );
  }
}

//------------------------------------------------------------------------------
inline int vtkStructuredGridConnectivity::Get1DOrientation(
        const int idx, const int ExtentLo, const int ExtentHi,
        const int OnLo, const int OnHi, const int NotOnBoundary )
{
  if( idx == ExtentLo )
  {
    return OnLo;
  }
  else if( idx == ExtentHi )
  {
    return OnHi;
  }
  return NotOnBoundary;
}

//------------------------------------------------------------------------------
inline bool vtkStructuredGridConnectivity::HasBlockConnection(
    const int gridID, const int blockDirection )
{
  // Sanity check
  assert("pre: gridID is out-of-bounds" &&
        (gridID >=0) && (gridID < static_cast<int>(this->NumberOfGrids)));
  assert("pre: BlockTopology has not been properly allocated" &&
        (this->NumberOfGrids == this->BlockTopology.size()));
  assert("pre: blockDirection is out-of-bounds" &&
        (blockDirection >= 0) && (blockDirection < 6) );
  bool status = false;
  if( this->BlockTopology[ gridID ] & (1 << blockDirection) )
  {
    status = true;
  }
  return( status );
}

//------------------------------------------------------------------------------
inline void vtkStructuredGridConnectivity::RemoveBlockConnection(
    const int gridID, const int blockDirection )
{
  // Sanity check
  assert("pre: gridID is out-of-bounds" &&
        (gridID >=0) && (gridID < static_cast<int>(this->NumberOfGrids)));
  assert("pre: BlockTopology has not been properly allocated" &&
        (this->NumberOfGrids == this->BlockTopology.size()));
  assert("pre: blockDirection is out-of-bounds" &&
        (blockDirection >= 0) && (blockDirection < 6) );

  this->BlockTopology[ gridID ] &= ~(1 << blockDirection);
}

//------------------------------------------------------------------------------
inline void vtkStructuredGridConnectivity::AddBlockConnection(
    const int gridID, const int blockDirection )
{
  // Sanity check
  assert("pre: gridID is out-of-bounds" &&
        (gridID >=0) && (gridID < static_cast<int>(this->NumberOfGrids)));
  assert("pre: BlockTopology has not been properly allocated" &&
        (this->NumberOfGrids == this->BlockTopology.size()));
  assert("pre: blockDirection is out-of-bounds" &&
        (blockDirection >= 0) && (blockDirection < 6) );
  this->BlockTopology[ gridID ] |= (1 << blockDirection);
}

//------------------------------------------------------------------------------
inline void vtkStructuredGridConnectivity::ClearBlockConnections(
    const int gridID )
{
  // Sanity check
  assert("pre: gridID is out-of-bounds" &&
        (gridID >=0) && (gridID < static_cast<int>(this->NumberOfGrids)));
  assert("pre: BlockTopology has not been properly allocated" &&
        (this->NumberOfGrids == this->BlockTopology.size()));
  for( int i=0; i < 6; ++i )
  {
    this->RemoveBlockConnection( gridID, i );
  } // END for all block directions
}

//------------------------------------------------------------------------------
inline int vtkStructuredGridConnectivity::GetNumberOfConnectingBlockFaces(
    const int gridID )
{
  // Sanity check
  assert("pre: gridID is out-of-bounds" &&
        (gridID >=0) && (gridID < static_cast<int>(this->NumberOfGrids)));
  assert("pre: BlockTopology has not been properly allocated" &&
        (this->NumberOfGrids == this->BlockTopology.size()));

  int count = 0;
  for( int i=0; i < 6; ++i )
  {
    if( this->HasBlockConnection( gridID, i ) )
    {
      ++count;
    }
  }
  assert( "post: count must be in [0,5]" && (count >=0 && count <= 6) );
  return( count );
}

//------------------------------------------------------------------------------
inline void vtkStructuredGridConnectivity::SetNumberOfGrids(
    const unsigned int N )
{
  if (N == 0)
  {
    vtkErrorMacro("Number of grids cannot be 0.");
    return;
  }

  this->NumberOfGrids = N;
  this->AllocateUserRegisterDataStructures();

  this->GridExtents.resize( 6*N,-1);
  this->Neighbors.resize( N );
  this->BlockTopology.resize( N );
}
#endif /* vtkStructuredGridConnectivity_h */