text: /*=========================================================================
  
  Program:   Visualization Toolkit
  Module:    vtkHyperTreeGrid.cxx
  
  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
  
  This software is distributed WITHOUT ANY WARRANTY; without even
  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  PURPOSE.  See the above copyright notice for more information.
  
  =========================================================================*/
  #include "vtkHyperTreeGrid.h"
  
  #include "vtkBitArray.h"
  #include "vtkCellData.h"
  #include "vtkCellLinks.h"
  #include "vtkCellType.h"
  #include "vtkCollection.h"
  #include "vtkDoubleArray.h"
  #include "vtkDataSetAttributes.h"
  #include "vtkGenericCell.h"
  #include "vtkHyperTree.h"
  #include "vtkHyperTreeCursor.h"
  #include "vtkIdTypeArray.h"
  #include "vtkInformation.h"
  #include "vtkInformationDoubleVectorKey.h"
  #include "vtkInformationIntegerKey.h"
  #include "vtkInformationVector.h"
  #include "vtkLine.h"
  #include "vtkMath.h"
  #include "vtkNew.h"
  #include "vtkObjectFactory.h"
  #include "vtkPixel.h"
  #include "vtkPointData.h"
  #include "vtkPoints.h"
  #include "vtkStructuredData.h"
  #include "vtkTimerLog.h"
  #include "vtkVoxel.h"
  
  #include <cassert>
  
  
  vtkInformationKeyMacro( vtkHyperTreeGrid, LEVELS, Integer );
  vtkInformationKeyMacro( vtkHyperTreeGrid, DIMENSION, Integer );
  vtkInformationKeyRestrictedMacro( vtkHyperTreeGrid, SIZES, DoubleVector, 3 );
  
  vtkStandardNewMacro( vtkHyperTreeGrid );
  
  vtkCxxSetObjectMacro( vtkHyperTreeGrid, MaterialMask, vtkBitArray );
  vtkCxxSetObjectMacro( vtkHyperTreeGrid, MaterialMaskIndex, vtkIdTypeArray );
  vtkCxxSetObjectMacro( vtkHyperTreeGrid, XCoordinates, vtkDataArray );
  vtkCxxSetObjectMacro( vtkHyperTreeGrid, YCoordinates, vtkDataArray );
  vtkCxxSetObjectMacro( vtkHyperTreeGrid, ZCoordinates, vtkDataArray );
  
  // Helpers to quickly fetch a HT at a given index or iterator
  #define GetHTGHyperTreeAtIndexMacro( _obj_, _index_ )             \
    ( static_cast<vtkHyperTree*>(                                   \
    _obj_->HyperTrees.find( _index_ ) != _obj_->HyperTrees.end() ?  \
    _obj_->HyperTrees[ _index_ ] : 0 ) )
  
  #define GetHyperTreeAtIndexMacro( _index_ ) \
    GetHTGHyperTreeAtIndexMacro( this, _index_ )
  
  //-----------------------------------------------------------------------------
  vtkHyperTreeGrid::vtkHyperTreeGrid()
  {
    // Dual grid corners (primal grid leaf centers)
    this->Points = 0;
    this->Connectivity = 0;
  
    // Internal links
    this->Links = 0;
  
    // Grid topology
    this->GridSize[0] = 0;
    this->GridSize[1] = 0;
    this->GridSize[2] = 0;
    this->TransposedRootIndexing = false;
  
    // Grid parameters
    this->BranchFactor = 2;
    this->Dimension =  1;
    this->NumberOfChildren = 2;
  
    // Masked primal leaves
    this->MaterialMask = vtkBitArray::New();
    this->MaterialMaskIndex = 0;
  
    // Grid geometry
    this->XCoordinates = vtkDoubleArray::New();
    this->YCoordinates = vtkDoubleArray::New();
    this->ZCoordinates = vtkDoubleArray::New();
  
    // For dataset API
    this->Voxel = vtkVoxel::New();
    this->Pixel = vtkPixel::New();
    this->Line = vtkLine::New();
    int extent[6];
    extent[0] = 0;
    extent[1] = this->GridSize[0] - 1;
    extent[2] = 0;
    extent[3] = this->GridSize[1] - 1;
    extent[4] = 0;
    extent[5] = this->GridSize[2] - 1;
  
    memcpy( this->Extent, extent, 6 * sizeof(int) );
  
    this->Information->Set( vtkDataObject::DATA_EXTENT_TYPE(), VTK_3D_EXTENT) ;
    this->Information->Set( vtkDataObject::DATA_EXTENT(), this->Extent, 6 );
  }
  
  //-----------------------------------------------------------------------------
  vtkHyperTreeGrid::~vtkHyperTreeGrid()
  {
    if ( this->MaterialMask )
      {
      this->MaterialMask->UnRegister( this );
      }
  
    if ( this->MaterialMaskIndex )
      {
      this->MaterialMaskIndex->UnRegister( this );
      }
  
    if ( this->XCoordinates )
      {
      this->XCoordinates->UnRegister( this );
      }
  
    if ( this->YCoordinates )
      {
      this->YCoordinates->UnRegister( this );
      }
  
    if ( this->ZCoordinates )
      {
      this->ZCoordinates->UnRegister( this );
      }
  
    if ( this->Voxel )
      {
      this->Voxel->UnRegister( this );
      }
  
    if ( this->Pixel )
      {
      this->Pixel->UnRegister( this );
      }
  
    if ( this->Line )
      {
      this->Line->UnRegister( this );
      }
  
    this->DeleteInternalArrays();
    this->DeleteTrees();
  }
  
  //-----------------------------------------------------------------------------
  void vtkHyperTreeGrid::PrintSelf( ostream& os, vtkIndent indent )
  {
    this->Superclass::PrintSelf( os, indent );
    os << indent << "Dimension: " << this->Dimension << endl;
    os << indent << "GridSize: "
       << this->GridSize[0] <<","
       << this->GridSize[1] <<","
       << this->GridSize[2] << endl;
  
    if ( this->XCoordinates )
      {
      this->XCoordinates->PrintSelf( os, indent.GetNextIndent() );
      }
    if ( this->YCoordinates )
      {
      this->YCoordinates->PrintSelf( os, indent.GetNextIndent() );
      }
    if ( this->ZCoordinates )
      {
      this->ZCoordinates->PrintSelf( os, indent.GetNextIndent() );
      }
  }
  
  //-----------------------------------------------------------------------------
  // Description:
  // Return what type of dataset this is.
  int vtkHyperTreeGrid::GetDataObjectType()
  {
    return VTK_HYPER_TREE_GRID;
  }
  
  //-----------------------------------------------------------------------------
  
  void vtkHyperTreeGrid::DeleteTrees()
  {
    if ( this->HyperTrees.size() )
      {
      vtkHyperTreeIterator it;
      this->InitializeTreeIterator( it );
      while ( vtkHyperTree* tree = it.GetNextTree() )
        {
        tree->UnRegister( this );
        }
      this->HyperTrees.clear();
      }
  }
  
  //-----------------------------------------------------------------------------
  // Description:
  // Copy the geometric and topological structure of a hyper tree grid
  // object.
  void vtkHyperTreeGrid::CopyStructure( vtkDataSet* ds )
  {
    assert( "pre: ds_exists" && ds!=0 );
    assert( "pre: same_type" && vtkHyperTreeGrid::SafeDownCast(ds)!=0 );
  
    vtkHyperTreeGrid* htg = vtkHyperTreeGrid::SafeDownCast( ds );
    assert(htg);
  
    // Copy grid parameters
    this->Dimension = htg->Dimension;
    this->BranchFactor = htg->BranchFactor;
    this->NumberOfChildren = htg->NumberOfChildren;
    this->TransposedRootIndexing = htg->TransposedRootIndexing;
    memcpy( this->GridSize, htg->GetGridSize(), 3 * sizeof( int ) );
  
    // Un-register existing trees
    DeleteTrees();
  
    // Shallow copy and register new trees
    this->HyperTrees = htg->HyperTrees;
  
    if ( this->HyperTrees.size() )
      {
      vtkHyperTreeIterator it;
      this->InitializeTreeIterator( it );
      while ( vtkHyperTree* tree = it.GetNextTree() )
        {
        tree->Register( this );
        }
      }
  
    this->DeleteInternalArrays();
  
    if ( htg->Points )
      {
      this->Points = htg->Points;
      this->Points->Register( this );
      }
  
    if ( htg->Connectivity )
      {
      this->Connectivity = htg->Connectivity;
      this->Connectivity->Register( this );
      }
  
    if ( htg->Links )
      {
      this->Links = htg->Links;
      this->Links->Register( this );
      }
  
    // Shallow copy masked leaf IDs
    this->SetMaterialMask( htg->MaterialMask );
  
    // Shallow copy masked leaf IDs
    this->SetMaterialMaskIndex( htg->MaterialMaskIndex );
  
    // Shallow copy coordinates
    this->SetXCoordinates( htg->XCoordinates );
    this->SetYCoordinates( htg->YCoordinates );
    this->SetZCoordinates( htg->ZCoordinates );
  }
  
  //----------------------------------------------------------------------------
  void vtkHyperTreeGrid::SetGridSize( unsigned int dim[3] )
  {
    this->SetGridExtent( 0, dim[0]-1, 0, dim[1]-1, 0, dim[2]-1 );
  }
  
  //----------------------------------------------------------------------------
  void vtkHyperTreeGrid::SetGridSize( unsigned int i, unsigned int j, unsigned int k )
  {
    this->SetGridExtent( 0, i-1, 0, j-1, 0, k-1 );
  }
  
  //----------------------------------------------------------------------------
  void vtkHyperTreeGrid::SetGridExtent( int extent[6] )
  {
    int description;
  
    description = vtkStructuredData::SetExtent( extent, this->Extent );
    if ( description < 0 ) //improperly specified
      {
      vtkErrorMacro ( << "Bad Extent, retaining previous values" );
      return;
      }
  
    if ( description == VTK_UNCHANGED )
      {
      return;
      }
  
    this->GridSize[0] = extent[1] - extent[0] + 1;
    this->GridSize[1] = extent[3] - extent[2] + 1;
    this->GridSize[2] = extent[5] - extent[4] + 1;
  
    this->Modified();
  }
  
  //----------------------------------------------------------------------------
  void vtkHyperTreeGrid::SetGridExtent(int iMin, int iMax, int jMin, int jMax,
                                       int kMin, int kMax)
  {
    int extent[6];
  
    extent[0] = iMin; extent[1] = iMax;
    extent[2] = jMin; extent[3] = jMax;
    extent[4] = kMin; extent[5] = kMax;
  
    this->SetGridExtent( extent );
  }
  
  //-----------------------------------------------------------------------------
  // Set the dimension of the tree with `dim'. See GetDimension() for details.
  // \pre valid_dim: dim >= 1 && dim <= 3
  // \post dimension_is_set: GetDimension()==dim
  void vtkHyperTreeGrid::SetDimension( unsigned int dim )
  {
    assert( "pre: valid_dim" && dim >= 1 && dim <= 3 );
    if ( this->Dimension == dim )
      {
      return;
      }
    this->Dimension = dim;
  
    // Number of children is factor^dimension
    this->NumberOfChildren = this->BranchFactor;
    for ( unsigned int i = 1; i < this->Dimension; ++ i )
      {
      this->NumberOfChildren *= this->BranchFactor;
      }
  
    this->Modified();
  }
  
  //-----------------------------------------------------------------------------
  // \pre valid_dim: factor == 2 or factor == 3;
  // \post dimension_is_set: GetBranchFactor()==dim
  void vtkHyperTreeGrid::SetBranchFactor( unsigned int factor )
  {
    assert( "pre: valid_factor" && factor >= 2 && factor <= 3 );
    if ( this->BranchFactor == factor )
      {
      return;
      }
    this->BranchFactor = factor;
  
    // Number of children is factor^dimension
    this->NumberOfChildren = this->BranchFactor;
    for ( unsigned int i = 1; i < this->Dimension; ++ i )
      {
      this->NumberOfChildren *= this->BranchFactor;
      }
  
    this->Modified();
  }
  
  //-----------------------------------------------------------------------------
  void vtkHyperTreeGrid::GenerateTrees()
  {
    // Clean up existing trees
    this->DeleteTrees();
  
    // Generate concrete instance of hyper tree and append it to list of roots
    vtkIdType nr = this->GetNumberOfTrees();
  
    for ( vtkIdType r = 0; r < nr; ++ r )
      {
      vtkIdType idx = this->MaterialMaskIndex ?
        this->MaterialMaskIndex->GetValue(r) : r;
      this->HyperTrees[ idx ] =
        vtkHyperTree::CreateInstance( this->BranchFactor, this->Dimension );
      }
  
    this->Modified();
    this->DeleteInternalArrays();
  }
  
  //----------------------------------------------------------------------------
  void vtkHyperTreeGrid::ComputeBounds()
  {
    // Retrieve coordinate arrays
    vtkDataArray* coords[3] =
      {
      this->XCoordinates,
      this->YCoordinates,
      this->ZCoordinates
      };
    for ( unsigned int i = 0; i < 3; ++ i )
      {
      if ( ! coords[i] || ! coords[i]->GetNumberOfTuples() )
        {
        return;
        }
      }
  
    // Get bounds from coordinate arrays
    vtkMath::UninitializeBounds( this->Bounds );
    for ( unsigned int i = 0; i < 3; ++ i )
      {
      unsigned int di = 2 * i;
      unsigned int dip = di + 1;
      this->Bounds[di] = coords[i]->GetComponent( 0, 0 );
      this->Bounds[dip] = coords[i]->GetComponent( coords[i]->GetNumberOfTuples() - 1, 0 );
  
      // Ensure that the bounds are increasing
      if ( this->Bounds[di] > this->Bounds[dip] )
        {
        std::swap( this->Bounds[di], this->Bounds[dip] );
        }
      }
  }
  
  //-----------------------------------------------------------------------------
  vtkIdType vtkHyperTreeGrid::GetNumberOfLeaves()
  {
    vtkIdType nLeaves = 0;
    vtkHyperTreeIterator it;
    this->InitializeTreeIterator( it );
    while ( vtkHyperTree* tree = it.GetNextTree() )
      {
      nLeaves += tree->GetNumberOfLeaves();
      }
    return nLeaves;
  }
  
  //-----------------------------------------------------------------------------
  vtkIdType vtkHyperTreeGrid::GetNumberOfLevels( vtkIdType id )
  {
    vtkHyperTree* tree = GetHyperTreeAtIndexMacro( id );
    return tree ? tree->GetNumberOfLevels() : 0;
  }
  
  //-----------------------------------------------------------------------------
  vtkIdType vtkHyperTreeGrid::GetNumberOfTrees()
  {
    return this->MaterialMaskIndex ?
      this->MaterialMaskIndex->GetNumberOfTuples() :
      this->GridSize[0] * this->GridSize[1] * this->GridSize[2];
  }
  
  //-----------------------------------------------------------------------------
  void vtkHyperTreeGrid::InitializeTreeIterator( vtkHyperTreeIterator& it )
  {
    it.Initialize( this );
  }
  
  //-----------------------------------------------------------------------------
  vtkHyperTreeCursor* vtkHyperTreeGrid::NewCursor( vtkIdType id )
  {
    vtkHyperTree* tree = GetHyperTreeAtIndexMacro( id );
    return tree ? tree->NewCursor() : 0;
  }
  
  //-----------------------------------------------------------------------------
  void vtkHyperTreeGrid::SubdivideLeaf( vtkHyperTreeCursor* leaf, vtkIdType id )
  {
    assert( "pre: leaf_exists" && leaf );
    assert( "pre: is_a_leaf" && leaf->IsLeaf() );
    vtkHyperTree* tree = GetHyperTreeAtIndexMacro( id );
    if ( tree )
      {
      tree->SubdivideLeaf( leaf );
      this->DeleteInternalArrays();
      }
  }
  
  //-----------------------------------------------------------------------------
  void vtkHyperTreeGrid::Initialize()
  {
    if ( this->HyperTrees.size() )
      {
      vtkHyperTreeIterator it;
      this->InitializeTreeIterator( it );
      while ( vtkHyperTree* tree = it.GetNextTree() )
        {
        tree->Initialize();
        }
      }
  
    this->DeleteInternalArrays();
  }
  
  //-----------------------------------------------------------------------------
  int vtkHyperTreeGrid::GetMaxCellSize()
  {
    switch( this->Dimension )
      {
      case 3:
        // Hexahedron, 8 vertices
        return 8;
      case 2:
        // Quadrangle, 4 vertices
        return 4;
      case 1:
        // Line segment, 2 vertices
        return 2;
      default:
        // This is useless, just to avoid a warning
        assert( "check: bad grid dimension" && 0 );
        return 0;
      }
  }
  
  //-----------------------------------------------------------------------------
  // Description:
  // Shallow and Deep copy.
  void vtkHyperTreeGrid::ShallowCopy( vtkDataObject* src )
  {
    assert( "src_same_type" && vtkHyperTreeGrid::SafeDownCast( src ) );
    this->CopyStructure(vtkHyperTreeGrid::SafeDownCast( src ) );
  
    // Call superclass
    this->Superclass::ShallowCopy( src );
  }
  
  //-----------------------------------------------------------------------------
  void vtkHyperTreeGrid::DeepCopy( vtkDataObject* src )
  {
    assert( "src_same_type" && vtkHyperTreeGrid::SafeDownCast( src ) );
    this->CopyStructure( vtkHyperTreeGrid::SafeDownCast( src ) );
  
    // Call superclass
    this->Superclass::DeepCopy( src );
  }
  
  //=============================================================================
  // DataSet API that returns dual grid.
  
  //-----------------------------------------------------------------------------
  // Description:
  // Return the number of leaves.
  // \post positive_result: result>=0
  vtkIdType vtkHyperTreeGrid::GetNumberOfCells()
  {
    this->ComputeDualGrid();
    return this->GetConnectivity()->GetNumberOfTuples();
  }
  
  //-----------------------------------------------------------------------------
  vtkIdType vtkHyperTreeGrid::GetNumberOfPoints()
  {
    return this->GetNumberOfLeaves();
  }
  
  //-----------------------------------------------------------------------------
  double* vtkHyperTreeGrid::GetPoint( vtkIdType ptId )
  {
    this->ComputeDualGrid();
    vtkPoints* leafCenters = this->GetPoints();
    assert( "Index out of bounds." &&
            ptId >= 0 && ptId < leafCenters->GetNumberOfPoints() );
    return leafCenters->GetPoint( ptId );
  }
  
  //-----------------------------------------------------------------------------
  void vtkHyperTreeGrid::GetPoint( vtkIdType ptId, double x[3] )
  {
    this->ComputeDualGrid();
    vtkPoints* leafCenters = this->GetPoints();
    assert( "Index out of bounds." &&
            ptId >= 0 && ptId < leafCenters->GetNumberOfPoints() );
    leafCenters->GetPoint( ptId, x );
  }
  
  //-----------------------------------------------------------------------------
  void vtkHyperTreeGrid::GetCell( vtkIdType cellId, vtkCell* cell )
  {
    assert( "Null cell ptr." && cell != 0 );
  
    int numPts = 1 << this->Dimension;
    double x[3];
  
    this->ComputeDualGrid();
    vtkIdTypeArray* cornerLeafIds = this->GetConnectivity();
    assert( "Index out of bounds." &&
            cellId >= 0 && cellId < cornerLeafIds->GetNumberOfTuples() );
    vtkPoints* leafCenters = this->GetPoints();
    vtkIdType* ptr = cornerLeafIds->GetPointer( 0 ) + cellId*numPts;
    for ( int ptIdx = 0; ptIdx < numPts; ++ptIdx, ++ptr )
      {
      cell->PointIds->SetId( ptIdx, *ptr );
      leafCenters->GetPoint( *ptr, x );
      cell->Points->SetPoint( ptIdx, x );
      }
  }
  
  //-----------------------------------------------------------------------------
  vtkCell* vtkHyperTreeGrid::GetCell( vtkIdType cellId )
  {
    vtkCell* cell = 0;
    switch ( this->Dimension )
      {
      case 1:
        cell = this->Line;
        break;
      case 2:
        cell = this->Pixel;
        break;
      case 3:
        cell = this->Voxel;
        break;
     default:
        assert( "post: bad grid dimension" && false );
        return 0;          // impossible case
      }
  
    GetCell( cellId, cell );
    return cell;
  }
  
  //-----------------------------------------------------------------------------
  void vtkHyperTreeGrid::GetCell( vtkIdType cellId, vtkGenericCell* cell )
  {
    assert( "GetCell on null cell." && cell != 0 );
    switch ( this->Dimension )
      {
      case 1:
        cell->SetCellTypeToLine();
        break;
      case 2:
        cell->SetCellTypeToPixel();
        break;
      case 3:
        cell->SetCellTypeToVoxel();
        break;
      default:
        assert( "post: bad grid dimension" && false );
        return;            // impossible case
      }
  
    GetCell( cellId, static_cast<vtkCell*>( cell ) );
  }
  
  //-----------------------------------------------------------------------------
  int vtkHyperTreeGrid::GetCellType( vtkIdType vtkNotUsed(cellId) )
  {
    switch ( this->Dimension )
      {
      case 3:
        return VTK_VOXEL; // hexahedron = 8 points
      case 2:
        return VTK_PIXEL; // quad = 4 points
      case 1:
        return VTK_LINE;  // line = 2 points
      default:
        assert( "post: bad grid dimension" && false );
        return 0;         // impossible case
      }
  }
  
  //-----------------------------------------------------------------------------
  void vtkHyperTreeGrid::GetCellPoints( vtkIdType cellId, vtkIdList* ptIds )
  {
    int numPts = 1 << this->Dimension;
    ptIds->Initialize();
    ptIds->SetNumberOfIds( numPts );
  
    this->ComputeDualGrid();
    vtkIdTypeArray* cornerLeafIds = this->GetConnectivity();
    assert( "Index out of bounds." &&
            cellId >= 0 && cellId < cornerLeafIds->GetNumberOfTuples() );
    vtkIdType* ptr = cornerLeafIds->GetPointer( 0 ) + cellId * numPts;
    memcpy( ptIds->GetPointer(0), ptr, numPts * sizeof( vtkIdType ) );
  }
  
  //----------------------------------------------------------------------------
  void vtkHyperTreeGrid::GetCellPoints( vtkIdType cellId,
                                        vtkIdType& npts,
                                        vtkIdType* &pts )
  {
    this->ComputeDualGrid();
    vtkIdTypeArray* cornerLeafIds = this->GetConnectivity();
    assert( "Index out of bounds." &&
            cellId >= 0 && cellId < cornerLeafIds->GetNumberOfTuples() );
    npts = static_cast<vtkIdType>( 1 << this->Dimension );
    pts = cornerLeafIds->GetPointer( 0 ) + cellId * npts;
  }
  
  //-----------------------------------------------------------------------------
  void vtkHyperTreeGrid::GetPointCells( vtkIdType ptId, vtkIdList* cellIds )
  {
    if ( ! this->Links )
      {
      this->BuildLinks();
      }
    cellIds->Reset();
  
    int numCells = this->Links->GetNcells( ptId );
    cellIds->SetNumberOfIds( numCells );
  
    vtkIdType* cells = this->Links->GetCells( ptId );
    for ( int i = 0; i < numCells; ++ i )
      {
      cellIds->SetId( i, cells[i] );
      }
  }
  
  //-----------------------------------------------------------------------------
  void vtkHyperTreeGrid::BuildLinks()
  {
    this->Links = vtkCellLinks::New();
    this->Links->Allocate( this->GetNumberOfPoints() );
    this->Links->Register( this );
    this->Links->BuildLinks( this );
    this->Links->UnRegister( this );
  }
  
  //-----------------------------------------------------------------------------
  void vtkHyperTreeGrid::GetCellNeighbors( vtkIdType cellId,
                                           vtkIdList* ptIds,
                                           vtkIdList* cellIds )
  {
    if ( ! this->Links )
      {
      this->BuildLinks();
      }
  
    cellIds->Reset();
  
    vtkIdType numPts = ptIds->GetNumberOfIds();
    if (numPts <= 0)
      {
      vtkErrorMacro("input point ids empty.");
      return;
      }
  
    int minNumCells = VTK_INT_MAX;
    vtkIdType* pts = ptIds->GetPointer( 0 );
    vtkIdType* minCells = 0;
    vtkIdType minPtId = 0;
    // Find the point used by the fewest number of cells
    for ( vtkIdType i = 0; i < numPts; i++ )
      {
      vtkIdType ptId = pts[i];
      int numCells = this->Links->GetNcells( ptId );
      if ( numCells < minNumCells )
        {
        minNumCells = numCells;
        minCells = this->Links->GetCells( ptId );
        minPtId = ptId;
        }
      }
  
    cellIds->Allocate( minNumCells );
  
    // For all cells that contNow for each cell, see if it contains all the points
    // in the ptIds list.
    for ( int i = 0; i < minNumCells; i++ )
      {
      // Do not include current cell
      if ( minCells[i] != cellId )
        {
        vtkIdType *cellPts;
        vtkIdType npts;
        this->GetCellPoints( minCells[i], npts, cellPts );
        // Iterate over all points in input cell
        bool match = true;
        for ( vtkIdType j = 0; j < numPts && match; j++ )
          {
          // Skip point with index minPtId which is contained by current cell
          if ( pts[j] != minPtId )
            {
            // Iterate over all points in current cell
            match = false;
            for ( vtkIdType k = 0; k < npts; ++ k )
              {
              if ( pts[j] == cellPts[k] )
                {
                // A match was found
                match = true;
                break;
                }
              } // For all points in current cell
            } // If not guaranteed match
          } // For all points in input cell
        if ( match )
          {
          cellIds->InsertNextId( minCells[i] );
          }
        } // If not the reference cell
      } // For all candidate cells attached to point
  }
  
  //----------------------------------------------------------------------------
  vtkIdType vtkHyperTreeGrid::FindPoint( double x[3] )
  {
    vtkIdType ix = 0;
    vtkIdType nx = this->XCoordinates->GetNumberOfTuples();
    while ( ix < nx && x[0] > this->XCoordinates->GetTuple1( ix ) )
      {
      ++ ix;
      }
    if ( ix )
      {
      -- ix;
      }
  
    vtkIdType iy = 0;
    vtkIdType ny = this->YCoordinates->GetNumberOfTuples();
    while ( iy < ny && x[1] > this->YCoordinates->GetTuple1( iy ) )
      {
      ++ iy;
      }
    if ( iy )
      {
      -- iy;
      }
  
    vtkIdType iz = 0;
    vtkIdType nz = this->ZCoordinates->GetNumberOfTuples();
    while ( iz < nz && x[2] > this->ZCoordinates->GetTuple1( iz ) )
      {
      ++ iz;
      }
    if ( iz )
      {
      -- iz;
      }
  
    int index = ( this->TransposedRootIndexing ) ?
      ( ix * this->GridSize[1] + iy ) * this->GridSize[2] + iz :
      ( iz * this->GridSize[1] + iy ) * this->GridSize[0] + ix;
  
    vtkHyperTreeSimpleCursor cursor;
    int pos[] = { 0, 0, 0 };
    cursor.Initialize( this, index, pos );
  
    // Geometry of the cell
    double origin[3] =
      {
      this->XCoordinates->GetTuple1( ix ),
      this->YCoordinates->GetTuple1( iy ),
      this->ZCoordinates->GetTuple1( iz )
      };
  
    double extreme[3] =
      {
      this->XCoordinates->GetTuple1( ix + 1 ),
      this->YCoordinates->GetTuple1( iy + 1 ),
      this->ZCoordinates->GetTuple1( iz + 1 )
      };
  
    double size[3] =
      {
      extreme[0] - origin[0],
      extreme[1] - origin[1],
      extreme[2] - origin[2]
      };
  
    return this->RecursiveFindPoint( x, &cursor, origin, size );
  }
  
  //----------------------------------------------------------------------------
  vtkIdType vtkHyperTreeGrid::RecursiveFindPoint( double x[3],
                                                  vtkHyperTreeSimpleCursor* cursor,
                                                  double* origin,
                                                  double* size )
  {
    if ( cursor->IsLeaf() )
      {
      return cursor->GetLeafIndex();
      }
  
    vtkHyperTreeSimpleCursor newCursor = *cursor;
    double newSize[3];
    double newOrigin[3];
    int child = 0;
    for ( int i = 0; i < 3; ++ i )
      {
      newSize[i] = size[i] * 0.5;
      newOrigin[i] = origin[i];
      if ( x[i] >= origin[i] + newSize[i] )
        {
        child = child | ( 1 << i );
        newOrigin[i] += newSize[i];
        }
      }
    newCursor.ToChild( child );
  
    return this->RecursiveFindPoint( x, &newCursor, newOrigin, newSize );
  }
  
  //----------------------------------------------------------------------------
  vtkIdType vtkHyperTreeGrid::FindCell( double x[3], vtkCell* cell,
                                        vtkGenericCell* gencell, vtkIdType cellId,
                                        double tol2, int& subId, double pcoords[3],
                                        double* weights )
  {
    vtkIdType ptId = this->FindPoint( x );
    if ( ptId < 0 )
      {
      // Return invalid Id if point is completely outside of data set
      return -1;
      }
  
    vtkNew<vtkIdList> cellIds;
    cellIds->Allocate( 8, 100 );
    this->GetPointCells( ptId, cellIds.GetPointer() );
    if ( cellIds->GetNumberOfIds() <= 0 )
      {
      return -1;
      }
  
    double closestPoint[3];
    double dist2;
    vtkIdType num = cellIds->GetNumberOfIds();
    for ( vtkIdType i = 0; i < num; ++ i )
      {
      cellId = cellIds->GetId( i );
      if ( gencell )
        {
        this->GetCell( cellId, gencell );
        }
      else
        {
        cell = this->GetCell( cellId );
        }
  
      // See whether this cell contains the point
      if ( ( gencell &&
             gencell->EvaluatePosition( x, closestPoint, subId,
                                        pcoords, dist2, weights ) == 1
             && dist2 <= tol2 )  ||
           ( !gencell &&
             cell->EvaluatePosition( x, closestPoint, subId,
                                     pcoords, dist2, weights ) == 1
             && dist2 <= tol2 ) )
        {
        return cellId;
        }
      }
  
    // This should never happen.
    vtkErrorMacro( "Could not find cell." );
    return -1;
  }
  
  //----------------------------------------------------------------------------
  vtkIdType vtkHyperTreeGrid::FindCell( double x[3], vtkCell* cell,
                                        vtkIdType cellId, double tol2,
                                        int& subId, double pcoords[3],
                                        double* weights )
  {
    return this->FindCell( x, cell, NULL, cellId, tol2, subId, pcoords, weights );
  }
  
  //----------------------------------------------------------------------------
  unsigned long vtkHyperTreeGrid::GetActualMemorySize()
  {
    unsigned long size = this->vtkDataSet::GetActualMemorySize();
  
    vtkHyperTreeIterator it;
    this->InitializeTreeIterator( it );
    while ( vtkHyperTree* tree = it.GetNextTree() )
      {
      size += tree->GetActualMemorySize();
      }
  
    // Approximate map memory size
    size += static_cast<unsigned long>(
      ( this->HyperTrees.size() * sizeof(vtkIdType) * 3 ) / 1024 );
  
    if ( this->XCoordinates )
      {
      size += this->XCoordinates->GetActualMemorySize();
      }
  
    if ( this->YCoordinates )
      {
      size += this->YCoordinates->GetActualMemorySize();
      }
  
    if ( this->ZCoordinates )
      {
      size += this->ZCoordinates->GetActualMemorySize();
      }
  
    if ( this->Points )
      {
      size += this->Points->GetActualMemorySize();
      }
  
    if ( this->Connectivity )
      {
      size += this->Connectivity->GetActualMemorySize();
      }
  
    if ( this->MaterialMask )
      {
      size += this->MaterialMask->GetActualMemorySize();
      }
  
    if ( this->MaterialMaskIndex )
      {
      size += this->MaterialMaskIndex->GetActualMemorySize();
      }
  
    return size;
  }
  
  //=============================================================================
  // Internal arrays used to generate dual grid.  Random access to cells
  // requires the cell leaves connectively array which costs memory.
  
  //-----------------------------------------------------------------------------
  vtkPoints* vtkHyperTreeGrid::GetPoints()
  {
    this->ComputeDualGrid();
    return this->Points;
  }
  
  //-----------------------------------------------------------------------------
  vtkIdTypeArray* vtkHyperTreeGrid::GetConnectivity()
  {
    this->ComputeDualGrid();
    return this->Connectivity;
  }
  
  //-----------------------------------------------------------------------------
  void vtkHyperTreeGrid::GetLevelZeroCoordsFromIndex( vtkIdType index,
                                                      unsigned int &i,
                                                      unsigned int &j,
                                                      unsigned int &k )
  {
    if ( !this->TransposedRootIndexing )
      {
      k = index / ( this->GridSize[0] * this->GridSize[1] );
      vtkIdType rk = k * ( this->GridSize[0] * this->GridSize[1] );
      j = ( index - rk ) / this->GridSize[0];
      i = index - ( j * this->GridSize[0] ) - rk;
      }
    else
      {
      i = index / ( this->GridSize[2] * this->GridSize[1] );
      vtkIdType rk = i * ( this->GridSize[2] * this->GridSize[1] );
      j = ( index - rk ) / this->GridSize[2];
      k = index - ( j * this->GridSize[2] ) - rk;
      }
  }
  
  //-----------------------------------------------------------------------------
  void vtkHyperTreeGrid::InitializeSuperCursor( vtkHyperTreeGridSuperCursor* sc,
                                                vtkIdType index )
  {
    unsigned int i, j, k;
    this->GetLevelZeroCoordsFromIndex( index, i, j, k );
    // Initialize center cursor
    this->InitializeSuperCursor( sc, i, j, k, index );
  }
  
  //-----------------------------------------------------------------------------
  void vtkHyperTreeGrid::InitializeSuperCursor( vtkHyperTreeGridSuperCursor* sc,
                                                unsigned int i,
                                                unsigned int j,
                                                unsigned int k,
                                                vtkIdType index )
  {
    // Location and size of the middle cursor/node
    double origin[3] =
      {
      this->XCoordinates->GetTuple1( i ),
      this->YCoordinates->GetTuple1( j ),
      this->ZCoordinates->GetTuple1( k )
      };
  
    double extreme[3] =
      {
      this->XCoordinates->GetTuple1( i + 1 ),
      this->YCoordinates->GetTuple1( j + 1 ),
      this->ZCoordinates->GetTuple1( k + 1 )
      };
  
    memcpy( sc->Origin, origin, 3 * sizeof( double ) );
    sc->Size[0] = extreme[0] - origin[0];
    sc->Size[1] = extreme[1] - origin[1];
    sc->Size[2] = extreme[2] - origin[2];
  
    // Initialize middle cursors and bounds for other cursors
    sc->NumberOfCursors = 3;
    for ( unsigned int c = 1; c < this->Dimension; ++ c )
      {
      sc->NumberOfCursors *= 3;
      }
    sc->MiddleCursorId = sc->NumberOfCursors / 2;
  
    int lowI = -1;
    int highI = 1;
    int lowJ = this->Dimension > 1 ? -1 : 0;
    int highJ = this->Dimension > 1 ? 1 : 0;
    int lowK = this->Dimension > 2 ? -1 : 0;
    int highK = this->Dimension > 2 ? 1 : 0;
    switch ( this->Dimension )
      {
      case 3:
        lowK = ( k == 0 ) ? 0 : -1;
        highK = ( k + 1 < this->GridSize[2] ) ? 1 : 0;
        VTK_FALLTHROUGH;
      case 2:
        lowJ = ( j == 0 ) ? 0 : -1;
        highJ = ( j + 1 < this->GridSize[1] ) ? 1 : 0;
        VTK_FALLTHROUGH;
      case 1:
        lowI = ( i == 0 ) ? 0 : -1;
        highI = ( i + 1 < this->GridSize[0] ) ? 1 : 0;
      }
  
    // Initialize all connectivity cursors
    for ( int ck = lowK; ck <= highK; ++ ck )
      {
      for ( int cj = lowJ; cj <= highJ; ++ cj )
        {
        for ( int ci = lowI; ci <= highI; ++ ci )
          {
          int pos[3] = { ci, cj, ck };
          int d = ci + 3 * cj + 9 * ck;
          sc->GetCursor( d )->Initialize( this, index, pos );
          } // i
        } // j
      } // k
  }
  
  //-----------------------------------------------------------------------------
  void vtkHyperTreeGrid::InitializeSuperCursorChild( vtkHyperTreeGridSuperCursor* parent,
                                                     vtkHyperTreeGridSuperCursor* child,
                                                     unsigned int childIdx )
  {
    // Retrieve child's parameters that are identical to parent's ones
    child->NumberOfCursors = parent->NumberOfCursors;
    child->MiddleCursorId = parent->MiddleCursorId;
  
    // Compute size of child
    child->Size[0] = parent->Size[0] / double( this->BranchFactor );
    child->Size[1] = parent->Size[1] / double( this->BranchFactor );
    child->Size[2] = parent->Size[2] / double( this->BranchFactor );
  
    // Compute origin of child
    unsigned int x, y, z;
    if ( this->BranchFactor == 2 )
      {
      x = childIdx & 1;
      y = ( childIdx & 2 ) >> 1;
      z = ( childIdx & 4 ) >> 2;
      }
    else
      {
      div_t d = div( childIdx, 9 );
      z = d.quot;
      y = d.rem / 3;
      x = childIdx % 3;
      }
    child->Origin[0] = parent->Origin[0] + ( x * child->Size[0] );
    child->Origin[1] = parent->Origin[1] + ( y * child->Size[1] );
    child->Origin[2] = parent->Origin[2] + ( z * child->Size[2] );
  
    // Move each cursor in the superCursor down to a child
    vtkSuperCursorEntry* cursorPtr =
      this->SuperCursorTraversalTable + ( childIdx * 27 );
    for ( int cursorIdx = 0; cursorIdx < child->NumberOfCursors; ++ cursorIdx )
      {
      // Extract the parent and child of the new node from the traversal table
      // Child is encoded in the first three bits for all dimensions
      int tParent = cursorPtr[cursorIdx].Parent;
      child->Cursors[cursorIdx] = parent->Cursors[tParent];
      if ( parent->Cursors[tParent].GetTree()
        && ! parent->Cursors[tParent].IsLeaf() )
        {
        // Move to child
        child->Cursors[cursorIdx] = parent->Cursors[tParent];
        int tChild = cursorPtr[cursorIdx].Child;
        child->Cursors[cursorIdx].ToChild( tChild );
        }
      }
  }
  
  //-----------------------------------------------------------------------------
  // Traverse tree with 3x3x3 super cursor.  Center cursor generates dual point
  // Smallest leaf (highest level) owns corners/dual cell.  Ties are given to
  // smallest index (z,y,x order)
  // post: Generate Points and Connectivity.
  void vtkHyperTreeGrid::ComputeDualGrid()
  {
    // Check if we can break out early
    vtkIdType numPoints = this->GetNumberOfPoints();
    if ( this->Points )
      {
      if ( this->Points->GetNumberOfPoints() == numPoints )
        {
        return;
        }
      this->Points->UnRegister( this );
      this->Connectivity->UnRegister( this );
      }
  
    vtkNew<vtkTimerLog> timer;
    timer->StartTimer();
  
    // Primal cell centers are dual points
    this->Points = vtkPoints::New();
    this->Points->SetNumberOfPoints( numPoints );
  
    this->Connectivity = vtkIdTypeArray::New();
    int numVerts = 1 << this->Dimension;
    this->Connectivity->SetNumberOfComponents( numVerts );
  
    // Create an array of cursors that occupy 1 3x3x3 neighborhood
    // Will traverse the tree as one
    // NB: Lower dimensions will not use them all
    this->GenerateSuperCursorTraversalTable();
  
    // Initialize grid depth
    vtkIdType gridDepth = 0;
  
    // Compute and assign scales of all tree roots
    double scale[] = { 1., 1., 1. };
  
    // Check whether coordinate arrays match grid size
    // If coordinates array are complete, compute all tree scales
    if ( static_cast<int>( this->GridSize[0] ) + 1 == this->XCoordinates->GetNumberOfTuples()
         && static_cast<int>( this->GridSize[1] ) + 1 == this->YCoordinates->GetNumberOfTuples()
         && static_cast<int>( this->GridSize[2] ) + 1 == this->ZCoordinates->GetNumberOfTuples() )
      {
      // Iterate over all hyper trees depending on indexing mode
  
      std::map<vtkIdType, vtkHyperTree*>::iterator it = this->HyperTrees.begin();
      std::map<vtkIdType, vtkHyperTree*>::iterator endit = this->HyperTrees.end();
  
      if ( this->TransposedRootIndexing )
        {
        // I-J-K indexing
        for ( unsigned int i = 0; i < this->GridSize[0] && it != endit; ++ i )
          {
          // Compute scale along x-axis
          scale[0] = this->XCoordinates->GetTuple1( i + 1 ) -
            this->XCoordinates->GetTuple1( i );
  
          for ( unsigned int j = 0; j < this->GridSize[1] && it != endit; ++ j )
            {
            // Compute scale along y-axis
            scale[1] = this->YCoordinates->GetTuple1( j + 1 ) -
              this->YCoordinates->GetTuple1( j );
  
            for ( unsigned int k = 0; k < this->GridSize[2] && it != endit; ++ k )
              {
              // Retrieve hyper tree
              vtkHyperTree* tree = it->second;
              ++ it;
              if ( !tree )
                {
                continue;
                }
  
              // Compute scale along z-axis and set tree scale
              scale[2] = this->ZCoordinates->GetTuple1( k + 1 ) -
                this->ZCoordinates->GetTuple1( k );
              tree->SetScale( scale );
  
              // Update hyper tree grid depth
              vtkIdType treeDepth = tree->GetNumberOfLevels();
              if ( treeDepth > gridDepth )
                {
                gridDepth = treeDepth;
                }
              } // i
            } // j
          } // k
        } // if ( this->TransposedRootIndexing )
      else
        {
        // K-J-I indexing
        for ( unsigned int k = 0; k < this->GridSize[2] && it != endit; ++ k )
          {
          // Compute scale along z-axis
          scale[2] = this->ZCoordinates->GetTuple1( k + 1 ) -
            this->ZCoordinates->GetTuple1( k );
  
          for ( unsigned int j = 0; j < this->GridSize[1] && it != endit; ++ j )
            {
            // Compute scale along y-axis
            scale[1] = this->YCoordinates->GetTuple1( j + 1 ) -
              this->YCoordinates->GetTuple1( j );
  
            for ( unsigned int i = 0; i < this->GridSize[0] && it != endit; ++ i )
              {
              // Retrieve hyper tree
              vtkHyperTree* tree = it->second;
              ++ it;
              if ( !tree )
                {
                continue;
                }
  
              // Compute scale along x-axis and set tree scale
              scale[0] = this->XCoordinates->GetTuple1( i + 1 ) -
                this->XCoordinates->GetTuple1( i );
              tree->SetScale( scale );
  
              // Update hyper tree grid depth
              vtkIdType treeDepth = tree->GetNumberOfLevels();
              if ( treeDepth > gridDepth )
                {
                gridDepth = treeDepth;
                }
              } // i
            } // j
          } // k
        } // else indexing mode
      } // if coordinate cardinalities match
  
    // Compute and store reduction factors for speed
    double factor = 1.;
    for ( unsigned short p = 0; p < gridDepth; ++ p )
      {
      this->ReductionFactors[p] = .5 * factor;
      factor /= this->BranchFactor;
      } // p
  
    // Iterate over all hyper trees
    vtkIdType index;
    vtkHyperTreeIterator it;
    this->InitializeTreeIterator( it );
    while ( it.GetNextTree( index ) )
      {
      // Storage for super cursors
      vtkHyperTreeGridSuperCursor superCursor;
  
      // Initialize center cursor
      this->InitializeSuperCursor( &superCursor, index );
  
      // Traverse and populate dual recursively
      this->TraverseDualRecursively( &superCursor, 0 );
      } // it
  
    // Adjust dual points as needed to fit the primal boundary
    for ( unsigned int d = 0; d < this->Dimension; ++ d )
      {
      // Iterate over all adjustments for current dimension
      for ( std::map<vtkIdType, double>::const_iterator itps =
        this->PointShifts[d].begin();
        itps != this->PointShifts[d].end(); ++ itps )
        {
        double pt[3];
        this->Points->GetPoint( itps->first, pt );
        pt[d] += itps->second;
        this->Points->SetPoint( itps->first, pt );
        } // it
      this->PointShifts[d].clear();
      } // d
    this->PointShifted.clear();
  
    timer->StopTimer();
    vtkDebugMacro(<< "Internal dual update : " << timer->GetElapsedTime());
  }
  
  //-----------------------------------------------------------------------------
  // Iterate over leaves.  Generate dual point.  Highest level (smallest leaf)
  // owns the corner and generates that dual cell.
  void vtkHyperTreeGrid::TraverseDualRecursively(
    vtkHyperTreeGridSuperCursor* superCursor, unsigned int level )
  {
    // Get cursor at super cursor center
    vtkHyperTreeSimpleCursor* cursor0 = superCursor->GetCursor( 0 );
  
    // Retrieve global index of center cursor
    vtkIdType id0 = cursor0->GetGlobalNodeIndex();
  
    // Initialize non leaf point coordinates
    double pt[3] = { 0., 0., 0. };
    for ( unsigned int d = 0; d < this->Dimension; ++ d )
      {
      pt[d] = superCursor->Origin[d] + superCursor->Size[d] / 2.;
      }
    // Retrieve global index of center cursor
    this->Points->SetPoint( id0, pt );
    if ( !cursor0->IsLeaf() )
      {
      // If cursor 0 is not at leaf, recurse to all children
      for ( unsigned int child = 0; child < this->NumberOfChildren; ++ child )
        {
        vtkHyperTreeGridSuperCursor newSuperCursor;
        this->InitializeSuperCursorChild( superCursor, &newSuperCursor, child );
        this->TraverseDualRecursively( &newSuperCursor, level + 1 );
        }
      }
    else
      {
      // Center is a leaf, create a dual point
      if ( this->MaterialMask->GetValue( id0 ) )
        {
        this->TraverseDualMaskedLeaf( superCursor );
        }
      else
        {
        this->TraverseDualLeaf( superCursor );
        }
      }
  }
  
  //-----------------------------------------------------------------------------
  void vtkHyperTreeGrid::TraverseDualMaskedLeaf(
    vtkHyperTreeGridSuperCursor* superCursor )
  {
    // Get cursor at super cursor center
    vtkHyperTreeSimpleCursor* cursor0 = superCursor->GetCursor( 0 );
  
    // Check across D-face neighbors whether point must be adjusted
    int f = 1;
    for ( unsigned int d = 0; d < this->Dimension; ++ d, f *= 3 )
      {
      // For each direction, check both orientations
      for ( int o = -1; o < 2; o += 2 )
        {
        // Retrieve face neighbor cursor
        vtkHyperTreeSimpleCursor* cursor = superCursor->GetCursor( o * f );
  
        // Detect faces shared by an unmasked cell, break ties at same level
        if ( cursor->GetTree()
             && cursor->IsLeaf()
             && cursor->GetLevel() < cursor0->GetLevel() )
          {
          vtkIdType id = cursor->GetGlobalNodeIndex();
          if ( ! this->MaterialMask->GetValue( id ) )
            {
            // Move to corresponding D-face
            this->PointShifted[id] = true;
            this->PointShifts[d][id] = - o
              * this->ReductionFactors[cursor->GetLevel()]
              * cursor->GetTree()->GetScale( d );
            }
          } // if cursor
        } // o
      } // d
  
    switch ( this->Dimension )
      {
      case 2:
        {
        // Check across (D-1)-face neighbors (corners)
        for ( int o2 = -1; o2 < 2; o2 += 2 )
          {
          int c = o2 + 3;
          for ( int o1 = -1; o1 < 2; o1 += 2 )
            {
            vtkHyperTreeSimpleCursor* cursor = superCursor->GetCursor( o1 * c );
            if ( cursor->GetTree()
                 &&
                 cursor->IsLeaf()
                 && cursor->GetLevel() < cursor0->GetLevel() )
              {
              vtkIdType id = cursor->GetGlobalNodeIndex();
              if ( ! this->MaterialMask->GetValue( id )
                   && ! this->PointShifted[id] )
                {
                // Move to corresponding corner
                double halfL[3];
                cursor->GetTree()->GetScale( halfL );
                double fac = this->ReductionFactors[cursor->GetLevel()];
                this->PointShifts[0][id] = - o1 * o2 * fac * halfL[0];
                this->PointShifts[1][id] = - o1 * fac * halfL[1];
                }
              } // if cursor
            } // o1
          } // o2
        break;
        } // case 2
      case 3:
        {
        // Check across (D-1)-face neighbors (edges)
        int tpa1 = 1;
        for ( int a1 = 0; a1 < 2; ++ a1, tpa1 *= 3 )
          {
          int tpa2 = 3 * tpa1;
          for ( int a2 = a1 + 1; a2 < 3; ++ a2, tpa2 *= 3 )
            {
            for ( int o2 = -1; o2 < 2; o2 += 2 )
              {
              int c  = tpa1 * o2 + tpa2;
              for ( int o1 = -1; o1 < 2; o1 += 2 )
                {
                vtkHyperTreeSimpleCursor* cursor = superCursor->GetCursor( o1 * c );
                if ( cursor->GetTree()
                     &&
                     cursor->IsLeaf()
                     && cursor->GetLevel() < cursor0->GetLevel() )
                  {
                  vtkIdType id = cursor->GetGlobalNodeIndex();
                  if ( ! this->MaterialMask->GetValue( id )
                       && ! this->PointShifted[id] )
                    {
                    // Move to corresponding edge
                    double halfL[3];
                    cursor->GetTree()->GetScale( halfL );
                    double fac = this->ReductionFactors[cursor->GetLevel()];
                    this->PointShifts[a1][id] = - o1 * o2 * fac * halfL[a1];
                    this->PointShifts[a2][id] = - o1 * fac * halfL[a2];
                    this->PointShifted[id] = true;
                    }
                  } // if cursor
                } // o1
              } // o2
            } // a2
          } // a1
  
        // Check across (D-2)-face neighbors (corners)
        for ( int o3 = -1; o3 < 2; o3 += 2 )
          {
          for ( int o2 = -1; o2 < 2; o2 += 2 )
            {
            int c = o2 * ( o3 + 3 ) + 9;
            for ( int o1 = -1; o1 < 2; o1 += 2 )
              {
              vtkHyperTreeSimpleCursor* cursor = superCursor->GetCursor( o1 * c );
              if ( cursor->GetTree()
                   &&
                   cursor->IsLeaf()
                   && cursor->GetLevel() < cursor0->GetLevel() )
                {
                vtkIdType id = cursor->GetGlobalNodeIndex();
                if ( ! this->MaterialMask->GetValue( id )
                     && ! this->PointShifted[id] )
                  {
                  // Move to corresponding corner
                  double halfL[3];
                  cursor->GetTree()->GetScale( halfL );
                  double fac = this->ReductionFactors[cursor->GetLevel()];
                  this->PointShifts[0][id] = - o1 * o2 * o3 * fac * halfL[0];
                  this->PointShifts[1][id] = - o1 * o2 * fac * halfL[1];
                  this->PointShifts[2][id] = - o1 * fac * halfL[2];
                  this->PointShifted[id] = true;
                  }
                } // if cursor
              } // o1
            } // o2
          } // o3
        break;
        } // case 3
      } // switch  ( this->Dimension )
  }
  
  //-----------------------------------------------------------------------------
  void vtkHyperTreeGrid::TraverseDualLeaf( vtkHyperTreeGridSuperCursor* superCursor )
  {
    // Get cursor at super cursor center
    vtkHyperTreeSimpleCursor* cursor0 = superCursor->GetCursor( 0 );
  
    // Initialize dual point coordinates and D-face adjustment flag
    double pt[] = { 0., 0., 0. };
    double halfL[] = { 0., 0., 0. };
    bool movedToDFace = false;
  
    // In 1D:
    //   (D-0)-faces are corners, neighbors are +/- 1
    //   (D-1)-faces do not exist
    //   (D-2)-faces do not exist
    // In 2D:
    //   (D-0)-faces are edges, neighbors are +/- 1, 3
    //   (D-1)-faces are corners, neighbors are +/- 2, 4
    //   (D-2)-faces do not exist
    // In 3D:
    //   (D-0)-faces are faces, neighbors are +/- 1, 3, 9
    //   (D-1)-faces are edges, neighbors are +/- 2, 4, 6, 8, 10, 12
    //   (D-2)-faces are corners, neighbors are +/-  5, 7, 11, 13
  
    // Check across D-face neighbors whether point must be adjusted
    int f = 1;
    for ( unsigned int d = 0; d < this->Dimension; ++ d, f *= 3 )
      {
      // Start at center
      halfL[d] = .5 * superCursor->Size[d];
      pt[d] = superCursor->Origin[d] + halfL[d];
  
      // Check
      for ( int o = -1; o < 2; o += 2 )
        {
        vtkHyperTreeSimpleCursor* cursor = superCursor->GetCursor( o * f );
        if ( ! cursor->GetTree()
             ||
             ( cursor->IsLeaf()
             && this->MaterialMask->GetValue( cursor->GetGlobalNodeIndex() ) ) )
          {
          // Move to corresponding D-face
          pt[d] += o * halfL[d];
          movedToDFace = true;
          } // if cursor
        } // o
      } // d
  
    // Only when point was not moved to D-face, check D-1 and D-2 neighbors
    if ( ! movedToDFace )
      {
      switch ( this->Dimension )
        {
        case 2:
          {
          // Check across (D-1)-face neighbors (corners)
          for ( int o2 = -1; o2 < 2; o2 += 2 )
            {
            int c = o2 + 3;
            for ( int o1 = -1; o1 < 2; o1 += 2 )
              {
              vtkHyperTreeSimpleCursor* cursor = superCursor->GetCursor( o1 * c );
              if ( ! cursor->GetTree()
                   ||
                   ( cursor->IsLeaf()
                     && this->MaterialMask->GetValue( cursor->GetGlobalNodeIndex() ) ) )
                {
                // Move to corresponding corner
                pt[0] += o1 * o2 * halfL[0];
                pt[1] += o1 * halfL[1];
                } // if cursor
              } // o1
            } // o2
          break;
          } // case 2
        case 3:
          {
          // Initialize edge adjustment flag
          bool movedToEdge = false;
  
          // Check across (D-1)-face neighbors (edges)
          int tpa1 = 1;
          for ( int a1 = 0; a1 < 2; ++ a1, tpa1 *= 3 )
            {
            int tpa2 = 3 * tpa1;
            for ( int a2 = a1 + 1; a2 < 3; ++ a2, tpa2 *= 3 )
              {
              for ( int o2 = -1; o2 < 2; o2 += 2 )
                {
                int c  = tpa1 * o2 + tpa2;
                for ( int o1 = -1; o1 < 2; o1 += 2 )
                  {
                  vtkHyperTreeSimpleCursor* cursor = superCursor->GetCursor( o1 * c );
                  if ( ! cursor->GetTree()
                       ||
                       ( cursor->IsLeaf()
                         && this->MaterialMask->GetValue( cursor->GetGlobalNodeIndex() ) ) )
                    {
                    // Move to corresponding edge
                    pt[a1] += o1 * o2 * halfL[a1];
                    pt[a2] += o1 * halfL[a2];
                    movedToEdge = true;
                    } // if cursor
                  } // o1
                } // o2
              } // a2
            } // a1
  
          // Only when point was not moved to edge, check across corners
          if ( ! movedToEdge )
            {
            for ( int o3 = -1; o3 < 2; o3 += 2 )
              {
              for ( int o2 = -1; o2 < 2; o2 += 2 )
                {
                int c = o2 * ( o3 + 3 ) + 9;
                for ( int o1 = -1; o1 < 2; o1 += 2 )
                  {
                  vtkHyperTreeSimpleCursor* cursor = superCursor->GetCursor( o1 * c );
                  if ( ! cursor->GetTree()
                       ||
                       ( cursor->IsLeaf()
                         && this->MaterialMask->GetValue( cursor->GetGlobalNodeIndex() ) ) )
                    {
                    // Move to corresponding corner
                    pt[0] += o1 * o2 * o3 * halfL[0];
                    pt[1] += o1 * o2 * halfL[1];
                    pt[2] += o1 * halfL[2];
                    } // if cursor
                  } // o1
                } // o2
              } // o3
            } // if ( ! movedToEdge )
          break;
          } // case 3
        } // switch ( this->Dimension )
      } // if ( ! movedToDFace )
  
    // Retrieve global index of center cursor
    vtkIdType id0 = cursor0->GetGlobalNodeIndex();
  
    // Insert dual point corresponding to current primal cell
    this->Points->SetPoint( id0, pt );
  
    // If cell is masked, terminate recursion, no dual cell will be generated
    if ( this->MaterialMask->GetValue( id0 ) )
      {
      return;
      }
  
    // Now see if the center leaf owns any of the corners
    // If it does, create the dual cell
    // Iterate over the corners around the middle leaf
    int numLeavesCorners = 1 << this->Dimension;
    for ( int cornerIdx = 0; cornerIdx < numLeavesCorners; ++ cornerIdx )
      {
      bool owner = true;
      vtkIdType leaves[8];
      // Iterate over every leaf touching the corner
      for ( int leafIdx = 0; leafIdx < numLeavesCorners && owner; ++ leafIdx )
        {
        // Compute the cursor index into the superCursor
        int cursorIdx = 0;
        switch ( this->Dimension )
          {
          case 3:
            cursorIdx += 9 * ( ( ( cornerIdx >> 2 ) & 1 ) + ( ( leafIdx >> 2 ) & 1 ) );
            VTK_FALLTHROUGH;
          case 2:
            cursorIdx += 3 * ( ( ( cornerIdx >> 1 ) & 1 ) + ( ( leafIdx >> 1 ) & 1 ) );
            VTK_FALLTHROUGH;
          case 1:
            cursorIdx += ( cornerIdx & 1) + ( leafIdx & 1);
          }
  
        // Collect the leaf indices for the dual cell
        leaves[leafIdx] = superCursor->Cursors[cursorIdx].GetGlobalNodeIndex();
  
        // Compute if the mid leaf owns the corner
        if ( cursorIdx != superCursor->MiddleCursorId )
          {
          vtkHyperTreeSimpleCursor* cursor = superCursor->Cursors + cursorIdx;
          if ( ! cursor->GetTree() || ! cursor->IsLeaf() )
            {
            // If neighbor leaf is out of bounds or has not been
            // refined to a leaf, this leaf does not own the corner
            owner = false;
            }
          else if ( this->MaterialMask->GetValue( cursor->GetGlobalNodeIndex() ) )
            {
            owner = false;
            }
          else if ( cursor->GetLevel() == cursor0->GetLevel()
                    && superCursor->MiddleCursorId < cursorIdx )
            {
            // A level tie is broken in favor of the largest index
            // All points are set before defining the cell
            owner = false;
            }
          }
        } // leafIdx
  
      if ( owner )
        {
        this->Connectivity->InsertNextTupleValue( leaves );
        }
      } // cornerIdx
  }
  
  //----------------------------------------------------------------------------
  // This table is used to move a 3x3x3 neighborhood of cursors through the tree.
  void vtkHyperTreeGrid::GenerateSuperCursorTraversalTable()
  {
    int bf = static_cast<int>( this->BranchFactor );
  
    int xChildDim = bf;
    int yChildDim = this->Dimension > 1 ? bf : 1;
    int zChildDim = this->Dimension > 2 ? bf : 1;
    int xCursorDim = 3;
    int yCursorDim = this->Dimension > 1 ? 3 : 1;
    int zCursorDim = this->Dimension > 2 ? 3 : 1;
  
    int childId = 0;
    for ( int zChild = 0; zChild < zChildDim; ++ zChild )
      {
      for ( int yChild = 0; yChild < yChildDim; ++ yChild )
        {
        for ( int xChild = 0; xChild < xChildDim; ++ xChild, ++ childId )
          {
          int cursorId = 0;
          for ( int zCursor = 0; zCursor < zCursorDim; ++ zCursor )
            {
            for ( int yCursor = 0; yCursor < yCursorDim; ++ yCursor )
              {
              for ( int xCursor = 0; xCursor < xCursorDim; ++ xCursor, ++ cursorId )
                {
                // Compute the x, y, z index into the
                // 6x6x6 (9x9x9) neighborhood of children.
                int xNeighbor = xCursor + xChild + xChildDim - 1;
                int yNeighbor = yCursor + yChild + yChildDim - 1;
                int zNeighbor = zCursor + zChild + zChildDim - 1;
  
                // Separate neighbor index into Cursor/Child index.
                div_t dx = div( xNeighbor, bf );
                div_t dy = div( yNeighbor, bf );
                div_t dz = div( zNeighbor, bf );
                int tableId = childId * 27 + cursorId;
                this->SuperCursorTraversalTable[tableId].Parent
                  = dx.quot + 3 * ( dy.quot + 3 * dz.quot );
                this->SuperCursorTraversalTable[tableId].Child
                  = dx.rem + bf * ( dy.rem + bf * dz.rem );
                } // xCursor
              } // yCursor
            } // zCursor
          } // xChild
        } // yChild
      } // zChild
  }
  
  //-----------------------------------------------------------------------------
  void vtkHyperTreeGrid::DeleteInternalArrays()
  {
    if ( this->Points )
      {
      this->Points->UnRegister( this );
      this->Points = 0;
      }
  
    if ( this->Connectivity )
      {
      this->Connectivity->UnRegister( this );
      this->Connectivity = 0;
      }
  
    if ( this->Links )
      {
      this->Links->UnRegister( this );
      this->Links = 0;
      }
  }
  
  //=============================================================================
  // Hyper tree grid iterator
  // Implemented here because it needs access to the internal classes.
  
  //-----------------------------------------------------------------------------
  void vtkHyperTreeGrid::vtkHyperTreeIterator::Initialize( vtkHyperTreeGrid* tree )
  {
    this->Tree = tree;
    this->Iterator = tree->HyperTrees.begin();
  }
  
  //-----------------------------------------------------------------------------
  vtkHyperTree* vtkHyperTreeGrid::vtkHyperTreeIterator::GetNextTree( vtkIdType &index )
  {
    if ( this->Iterator == this->Tree->HyperTrees.end() )
      {
      return 0;
      }
    vtkHyperTree* t = this->Iterator->second;
    index = this->Iterator->first;
    ++ this->Iterator;
    return t;
  }
  
  //-----------------------------------------------------------------------------
  vtkHyperTree* vtkHyperTreeGrid::vtkHyperTreeIterator::GetNextTree()
  {
    vtkIdType index;
    return GetNextTree( index );
  }
  
  //=============================================================================
  // Hyper tree grid cursor
  // Implemented here because it needs access to the internal classes.
  
  //-----------------------------------------------------------------------------
  // Constructor.
  vtkHyperTreeGrid::vtkHyperTreeSimpleCursor::vtkHyperTreeSimpleCursor()
  {
    this->Clear();
  }
  
  //-----------------------------------------------------------------------------
  // Set the state back to the initial constructed state
  void vtkHyperTreeGrid::vtkHyperTreeSimpleCursor::Clear()
  {
    this->Tree = 0;
    this->Index = 0;
    this->Leaf = false;
    this->Level = 0;
  }
  
  //-----------------------------------------------------------------------------
  void vtkHyperTreeGrid::vtkHyperTreeSimpleCursor::Initialize( vtkHyperTreeGrid* grid,
                                                               vtkIdType index,
                                                               int pos[3] )
  {
    // Convert local index into global one
    unsigned int n[3];
    grid->GetGridSize( n );
    vtkIdType globalIndex = grid->GetTransposedRootIndexing() ?
      index + pos[2] +
      pos[1] * static_cast<int>( n[2] ) +
      pos[0] * static_cast<int>( n[2] ) * static_cast<int>( n[1] ) :
      index + pos[0] +
      pos[1] * static_cast<int>( n[0] ) +
      pos[2] * static_cast<int>( n[0] ) * static_cast<int>( n[1] );
  
    // Assign hypertree to this cursor
    this->Tree = GetHTGHyperTreeAtIndexMacro( grid, globalIndex );
    this->ToRoot();
  }
  
  //-----------------------------------------------------------------------------
  bool vtkHyperTreeGrid::vtkHyperTreeSimpleCursor::IsLeaf()
  {
    // Empty cursors appear like a leaf so that recursion stop
    return this->Tree ? this->Leaf : true;
  }
  
  //-----------------------------------------------------------------------------
  void vtkHyperTreeGrid::vtkHyperTreeSimpleCursor::ToRoot()
  {
    if ( ! this->Tree )
      {
      return;
      }
  
    // Return to root level
    this->Level = 0;
    this->Index = 0;
  
    this->Leaf = ( this->Tree->GetNumberOfLeaves() == 1 );
  }
  
  //-----------------------------------------------------------------------------
  void vtkHyperTreeGrid::vtkHyperTreeSimpleCursor::ToChild( int child )
  {
    if ( ! this->Tree || this->Leaf )
      {
      // Leaves do not have children.
      return;
      }
  
    this->Tree->FindChildParameters( child, this->Index, this->Leaf );
    this->Level++;
  
    assert( "Bad index" && this->Index >= 0 );
  }
  
  //-----------------------------------------------------------------------------
  vtkIdType vtkHyperTreeGrid::vtkHyperTreeSimpleCursor::GetGlobalNodeIndex()
  {
    return this->Tree ?
      this->Tree->GetGlobalIndexFromLocal( this->Index ) : 0;
  }
  