vtkStructuredImplicitConnectivity.cxx 48.1 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkStructuredImplicitConnectivity.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.

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

// VTK includes
18
#include "vtkDataArray.h"
19
#include "vtkFieldDataSerializer.h"
20
#include "vtkImageData.h"
21 22 23 24 25 26
#include "vtkMPIController.h"
#include "vtkMultiProcessController.h"
#include "vtkMultiProcessStream.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkPoints.h"
27
#include "vtkRectilinearGrid.h"
28 29 30 31 32 33 34 35 36 37 38 39 40 41 42
#include "vtkStructuredData.h"
#include "vtkStructuredExtent.h"
#include "vtkStructuredGrid.h"

// C/C++ includes
#include <algorithm>
#include <cassert>
#include <map>
#include <sstream>
#include <vector>

//==============================================================================
// INTERNAL DATASTRUCTURES & DEFINITIONS
//==============================================================================

Kunda's avatar
Kunda committed
43
// Some useful extent macros
44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87
#define IMIN(ext) ext[0]
#define IMAX(ext) ext[1]
#define JMIN(ext) ext[2]
#define JMAX(ext) ext[3]
#define KMIN(ext) ext[4]
#define KMAX(ext) ext[5]

#define I(ijk) ijk[0]
#define J(ijk) ijk[1]
#define K(ijk) ijk[2]

namespace vtk
{
namespace detail
{

// Given two intervals A=[a1,a2] and B[b1,b2] the IntervalsConnect struct
// enumerates the cases where interval A connects to Interval B.
struct IntervalsConnect
{
  // NOTE: This enum is arranged s.t., negating a value in [-4,4] will yield
  // the mirror inverse
  enum connectivity_t
  {
    IMPLICIT_LO = -4, // Interval A implicitly connects with B on A's low end
    SUBSET      = -3, // Interval A is completely inside interval B
    OVERLAP_LO  = -2, // Interval A intersects with B on A's low end
    LO          = -1, // A's low end touches B's high end A.Low() == B.High()
    ONE_TO_ONE  =  0, // Intervals A,B are exactly the same.
    HI          =  1, // A's high end touches B's low end A.High() == B.Low()
    OVERLAP_HI  =  2, // Interval A intersects with B on A's high end
    SUPERSET    =  3, // Interval A *contains* all of interval B
    IMPLICIT_HI =  4, // Interval A implicitly connects with B on its high end.

    DISJOINT    =  5, // Intervals A,B are completely disjoint.
    UNDEFINED   =  6  // Undefined
  };

  static
  std::string OrientationToString( int orient[3] )
  {
    std::ostringstream oss;
    oss << "(";
    for(int i=0; i < 3; ++i)
88
    {
89
      if(i==1 || i==2)
90
      {
91
        oss << ", ";
92
      }
93
      switch( orient[i] )
94
      {
95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129
        case IMPLICIT_LO:
          oss << "IMPLICIT_LO";
          break;
        case SUBSET:
          oss << "SUBSET";
          break;
        case OVERLAP_LO:
          oss << "OVERLAP_LO";
          break;
        case LO:
          oss << "LO";
          break;
        case ONE_TO_ONE:
          oss << "ONE_TO_ONE";
          break;
        case HI:
          oss << "HI";
          break;
        case OVERLAP_HI:
          oss << "OVERLAP_HI";
          break;
        case SUPERSET:
          oss << "SUPERSET";
          break;
        case IMPLICIT_HI:
          oss << "IMPLICIT_HI";
          break;
        case DISJOINT:
          oss << "DISJOINT";
          break;
        case UNDEFINED:
          oss << "UNDEFINED";
          break;
        default:
          oss << "*UNKNOWN*";
130 131
      } // END switch
    } // END for
132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176
    oss << ")";
    return( oss.str() );
  }

}; // END struct IntervalsConnect

//------------------------------------------------------------------------------
//  Interval class Definition
//------------------------------------------------------------------------------
class Interval
{
public:
  Interval() : lo(0), hi(-1) {};
  Interval(const int l, const int h) : lo(l), hi(h) {};
  ~Interval() {};

  int Low() const { return this->lo; };
  int High() const { return this->hi; };
  int Cardinality() const { return(this->hi-this->lo+1); };
  bool Valid() const { return(this->lo <= this->hi); };
  void Set(const int l, const int h) { this->lo=l; this->hi=h; };
  void Invalidate() {this->Set(0,-1);}
  bool Within(const Interval& B) const
    { return( (this->lo >= B.Low()) && (this->hi <= B.High()) ); };

  bool ImplicitNeighbor(const Interval& B, int& type);
  static bool ImplicitNeighbors(
      const Interval& A, const Interval& B, int& type);

  bool Intersects(const Interval& B, Interval& Overlap, int& type);
  static bool Intersects(const Interval& A, const Interval& B,
                         Interval& Overlap, int& type);
private:
  int lo;
  int hi;
};

//------------------------------------------------------------------------------
bool Interval::ImplicitNeighbors(const Interval& A, const Interval& B, int& t)
{
  assert("pre: interval is not valid!" && A.Valid());
  assert("pre: B interval is not valid!" && B.Valid() );

  bool status = false;
  if( A.High()+1 == B.Low() )
177
  {
178 179
    status = true;
    t = IntervalsConnect::IMPLICIT_HI;
180
  }
181
  else if( B.High()+1 == A.Low() )
182
  {
183 184
    status = true;
    t = IntervalsConnect::IMPLICIT_LO;
185
  }
186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205
  return( status );
}

//------------------------------------------------------------------------------
bool Interval::ImplicitNeighbor(const Interval& B, int& type)
{
  return( Interval::ImplicitNeighbors(*this,B,type) );
}

//------------------------------------------------------------------------------
bool Interval::Intersects(const Interval& A, const Interval& B,
                          Interval& Overlap, int& type)
{
  assert("pre: interval is not valid!" && A.Valid());
  assert("pre: B interval is not valid!" && B.Valid() );

  bool status = false;

  // Disjoint cases
  if( A.High() < B.Low() )
206
  {
207 208 209
    type = IntervalsConnect::DISJOINT;
    Overlap.Invalidate();
    status = false;
210
  }
211
  else if( B.High() < A.Low() )
212
  {
213 214 215
    type = IntervalsConnect::DISJOINT;
    Overlap.Invalidate();
    status = false;
216
  }
217 218 219 220
  // ONE_TO_ONE case
  else if( A.Cardinality()==B.Cardinality() &&
           A.Low()==B.Low() &&
           A.High()==B.High() )
221
  {
222 223 224
    type = IntervalsConnect::ONE_TO_ONE;
    Overlap.Set(A.Low(),A.High());
    status = true;
225
  }
226 227
  // A is a SUBSET of B
  else if( A.Within(B) )
228
  {
229 230 231
    type = IntervalsConnect::SUBSET;
    Overlap.Set(A.Low(),A.High());
    status = true;
232
  }
233 234
  // A is a superset of B
  else if( B.Within(A) )
235
  {
236 237 238
    type = IntervalsConnect::SUPERSET;
    Overlap.Set(B.Low(),B.High());
    status = true;
239
  }
240 241
  // A touches B on the high end
  else if( A.High() == B.Low() )
242
  {
243 244 245
    type = IntervalsConnect::HI;
    Overlap.Set(A.High(),A.High());
    status = true;
246
  }
247 248
  // A touches B on the low end
  else if( A.Low() == B.High() )
249
  {
250 251 252
    type = IntervalsConnect::LO;
    Overlap.Set(A.Low(),A.Low());
    status = true;
253
  }
254 255
  // A intersects B on its low end
  else if( (A.Low() >= B.Low()) && (A.Low() <= B.High()) )
256
  {
257 258 259
    type = IntervalsConnect::OVERLAP_LO;
    Overlap.Set(A.Low(),B.High());
    status = true;
260
  }
261 262
  // A intersects B on its high end
  else if( (A.High() >= B.Low() ) && (A.High() <= B.High()) )
263
  {
264 265 266
    type = IntervalsConnect::OVERLAP_HI;
    Overlap.Set(B.Low(),A.High());
    status = true;
267
  }
268
  else
269
  {
270 271 272 273 274 275
    vtkGenericWarningMacro(
        << "Undefined interval intersection!"
        << "Code should not reach here!!!");
    type = IntervalsConnect::UNDEFINED;
    status = false;
    Overlap.Invalidate();
276
  }
277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294
  return( status );
}

//------------------------------------------------------------------------------
bool Interval::Intersects(const Interval& B, Interval& Overlap, int& type)
{
  return( Interval::Intersects(*this, B, Overlap, type) );
}

//------------------------------------------------------------------------------
struct ImplicitNeighbor
{
  int Rank;           // the rank of the neighbor
  int Extent[6];      // the extent of the neighbor
  int Orientation[3]; // the orientation w.r.t the local extent
  int Overlap[6];     // the overlap extent

  std::string ToString()
295
  {
296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316
    std::ostringstream oss;

    oss << "rank=" << this->Rank << " ";
    oss << "extent=[";
    oss << this->Extent[0] << ", ";
    oss << this->Extent[1] << ", ";
    oss << this->Extent[2] << ", ";
    oss << this->Extent[3] << ", ";
    oss << this->Extent[4] << ", ";
    oss << this->Extent[5] << "] ";
    oss << "overlap=[";
    oss << this->Overlap[0] << ", ";
    oss << this->Overlap[1] << ", ";
    oss << this->Overlap[2] << ", ";
    oss << this->Overlap[3] << ", ";
    oss << this->Overlap[4] << ", ";
    oss << this->Overlap[5] << "] ";
    oss << "orientation=";
    oss << IntervalsConnect::OrientationToString(this->Orientation);

    return( oss.str() );
317
  }
318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352
};

//------------------------------------------------------------------------------
struct DomainMetaData
{
  int WholeExtent[6];     // Extent of the entire domain

  int DataDescription;    // Data-description of the distributed dataset.
  int NDim;               // Number of dimensions according to DataDescription.
  int DimIndex[3];        // Stores the dimensions of the dataset in the
                          // the right order. This essentially allows to
                          // process 2-D (XY,XZ,YZ) and 3-D datasets in a
                          // transparent way.

  int GlobalImplicit[3];  // indicates for each dimension if there is globally
                          // implicit connectivity. Any value > 0 indicates
                          // implicit connectivity in the given direction.

  // Flat list of extents. Extents are organized as follows:
  // [id, imin, imax, jmin, jmax, kmin, kmax]
  std::vector< int > ExtentListInfo;

  /// \brief Checks if a grid with the given extent is within this domain
  /// \param ext the extent of the grid in query
  /// \return status true if the grid is insided, else false.
  bool HasGrid(int ext[6])
    { return( vtkStructuredExtent::Smaller(ext,this->WholeExtent) );  };

  /// \brief Initializes the domain metadata.
  void Initialize(int wholeExt[6])
  {
    memcpy(this->WholeExtent,wholeExt,6*sizeof(int));
    this->DataDescription =
        vtkStructuredData::GetDataDescriptionFromExtent(wholeExt);

353
    if (this->DataDescription == VTK_EMPTY)
354
    {
355
      return;
356
    }
357

358 359 360 361 362 363 364 365 366 367 368
    // Sanity checks!
    assert( "pre: data description is VTK_EMPTY!" &&
             (this->DataDescription != VTK_EMPTY) );
    assert( "pre: dataset must be 2-D or 3-D" &&
            (this->DataDescription >= VTK_XY_PLANE) );

    this->NDim = -1;
    std::fill(this->DimIndex,this->DimIndex+3,-1);
    std::fill(this->GlobalImplicit,this->GlobalImplicit+3,0);

    switch( this->DataDescription )
369
    {
370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394
      case VTK_XY_PLANE:
        this->NDim        = 2;
        this->DimIndex[0] = 0;
        this->DimIndex[1] = 1;
        break;
      case VTK_XZ_PLANE:
        this->NDim        = 2;
        this->DimIndex[0] = 0;
        this->DimIndex[1] = 2;
        break;
      case VTK_YZ_PLANE:
        this->NDim        = 2;
        this->DimIndex[0] = 1;
        this->DimIndex[1] = 2;
        break;
      case VTK_XYZ_GRID:
        this->NDim        = 3;
        this->DimIndex[0] = 0;
        this->DimIndex[1] = 1;
        this->DimIndex[2] = 2;
        break;
      default:
        vtkGenericWarningMacro(
            << "Cannot handle data description: "
            << this->DataDescription << "\n");
395
    } // END switch
396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414

    assert( "post: NDim==2 || NDim==3" && ( this->NDim==2 || this->NDim==3 ) );
  }

};

//------------------------------------------------------------------------------
struct StructuredGrid
{
  int ID;
  int Extent[6];
  int DataDescription;

  int Grow[3];     // indicates if the grid grows to the right along each dim.
  int Implicit[3]; // indicates implicit connectivity alone each dim.

  vtkPoints* Nodes;
  vtkPointData* PointData;

415 416 417 418 419
  // arrays used if the grid is a rectilinear grid
  vtkDataArray* X_Coords;
  vtkDataArray* Y_Coords;
  vtkDataArray* Z_Coords;

420 421
  std::vector< ImplicitNeighbor > Neighbors;

422 423 424
//------------------------------------------------------------------------------
  bool IsRectilinearGrid()
  {
425 426 427
    if( (this->X_Coords != nullptr) &&
        (this->Y_Coords != nullptr) &&
        (this->Z_Coords != nullptr) )
428
    {
429
      return true;
430
    }
431 432 433 434 435 436
    return false;
  }

//------------------------------------------------------------------------------
  void Clear()
  {
437
    if( this->Nodes != nullptr )
438
    {
439
      this->Nodes->Delete();
440
      this->Nodes = nullptr;
441
    }
442
    if( this->PointData != nullptr )
443
    {
444
      this->PointData->Delete();
445
      this->PointData = nullptr;
446
    }
447
    if( this->X_Coords != nullptr )
448
    {
449
      this->X_Coords->Delete();
450
      this->X_Coords = nullptr;
451
    }
452
    if( this->Y_Coords != nullptr )
453
    {
454
      this->Y_Coords->Delete();
455
      this->Y_Coords = nullptr;
456
    }
457
    if( this->Z_Coords != nullptr )
458
    {
459
      this->Z_Coords->Delete();
460
      this->Z_Coords = nullptr;
461
    }
462 463 464
    this->Neighbors.clear();
  }

465 466 467
//------------------------------------------------------------------------------
  void Initialize(StructuredGrid* grid)
  {
468
    assert("pre: input grid is nullptr!" && (grid != nullptr) );
469

470
    this->Initialize(grid->ID,grid->Extent,nullptr,nullptr);
471 472 473

    // Grow the extent in each dimension as needed
    for(int i=0; i < 3; ++i)
474
    {
475
      if( grid->Grow[i]==1 )
476
      {
477
        this->Extent[i*2+1] += 1;
478 479
      } // END if
    } // END for all dimensions
480 481

    // the number of nodes in the grown extent
482
    vtkIdType nnodes = vtkStructuredData::GetNumberOfPoints(
483 484 485
        this->Extent,grid->DataDescription);

   // Allocate coordinates, if needed
486
   if( grid->Nodes != nullptr )
487
   {
488 489 490
     this->Nodes = vtkPoints::New();
     this->Nodes->SetDataType( grid->Nodes->GetDataType() );
     this->Nodes->SetNumberOfPoints( nnodes );
491
   } // END if has points
492
   else
493
   {
494
     this->Nodes = nullptr;
495
   }
496

497
   // Allocate rectilinear grid coordinates, if needed
498 499 500
   if( (grid->X_Coords != nullptr) &&
       (grid->Y_Coords != nullptr) &&
       (grid->Z_Coords != nullptr) )
501
   {
502 503 504 505 506 507 508 509
     int dims[3];
     vtkStructuredData::GetDimensionsFromExtent(
         this->Extent,dims,this->DataDescription);

     this->X_Coords = vtkDataArray::CreateDataArray(
         grid->X_Coords->GetDataType());
     this->X_Coords->SetNumberOfTuples( dims[0] );
     for(vtkIdType idx=0; idx < grid->X_Coords->GetNumberOfTuples(); ++idx)
510
     {
511
       this->X_Coords->SetTuple(idx,idx,grid->X_Coords);
512
     }
513 514 515 516 517

     this->Y_Coords = vtkDataArray::CreateDataArray(
         grid->Y_Coords->GetDataType());
     this->Y_Coords->SetNumberOfTuples( dims[1] );
     for(vtkIdType idx=0; idx < grid->Y_Coords->GetNumberOfTuples(); ++idx)
518
     {
519
        this->Y_Coords->SetTuple(idx,idx,grid->Y_Coords);
520
     }
521 522 523 524 525

     this->Z_Coords = vtkDataArray::CreateDataArray(
         grid->Z_Coords->GetDataType());
     this->Z_Coords->SetNumberOfTuples( dims[2] );
     for(vtkIdType idx=0; idx < grid->Z_Coords->GetNumberOfTuples(); ++idx)
526
     {
527
        this->Z_Coords->SetTuple(idx,idx,grid->Z_Coords);
528 529
     }
   } // END if rectilinear grid
530
   else
531
   {
532 533 534
     grid->X_Coords = nullptr;
     grid->Y_Coords = nullptr;
     grid->Z_Coords = nullptr;
535
   }
536

537
   // Allocate fields, if needed
538
   if( grid->PointData != nullptr )
539
   {
540 541 542 543 544 545
     this->PointData = vtkPointData::New();
     this->PointData->CopyAllocate(grid->PointData,nnodes);

     // NOTE: CopyAllocate, allocates the buffers internally, but, does not
     // set the number of tuples of each array to nnodes.
     for(int array=0; array < this->PointData->GetNumberOfArrays(); ++array)
546
     {
547 548
       vtkDataArray* a = this->PointData->GetArray( array );
       a->SetNumberOfTuples(nnodes);
549
     } // END for all arrays
550

551
   }
552
   else
553
   {
554
     this->PointData = nullptr;
555
   }
556 557 558 559 560 561

   // copy everything from the given grid
   int desc   = grid->DataDescription;
   int ijk[3] = {0,0,0};

   for( I(ijk)=IMIN(grid->Extent); I(ijk) <= IMAX(grid->Extent); ++I(ijk) )
562
   {
563
     for( J(ijk)=JMIN(grid->Extent); J(ijk) <= JMAX(grid->Extent); ++J(ijk) )
564
     {
565
       for( K(ijk)=KMIN(grid->Extent); K(ijk) <= KMAX(grid->Extent); ++K(ijk) )
566
       {
567 568 569 570 571 572 573 574 575
         // Compute the source index
         vtkIdType srcIdx =
             vtkStructuredData::ComputePointIdForExtent(grid->Extent,ijk,desc);

         // Compute the target index
         vtkIdType targetIdx =
             vtkStructuredData::ComputePointIdForExtent(this->Extent,ijk,desc);

         // Copy nodes
576
         if( this->Nodes != nullptr )
577
         {
578
           this->Nodes->SetPoint(targetIdx,grid->Nodes->GetPoint(srcIdx));
579
         }
580 581

         // Copy node-centered fields
582
         if( this->PointData != nullptr )
583
         {
584
           this->PointData->CopyData(grid->PointData,srcIdx,targetIdx);
585
         }
586

587 588 589
       } // END for all k
     } // END for all j
   } // END for all i
590 591 592

  }

593 594 595 596
//------------------------------------------------------------------------------
  void Initialize(int id, int ext[6], vtkDataArray* x_coords,
      vtkDataArray* y_coords, vtkDataArray* z_coords, vtkPointData* fields)
  {
597 598 599
    assert("pre: nullptr x_coords!" && (x_coords != nullptr) );
    assert("pre: nullptr y_coords!" && (y_coords != nullptr) );
    assert("pre: nullptr z_coords!" && (z_coords != nullptr) );
600 601 602 603 604 605 606

    this->ID = id;
    memcpy(this->Extent,ext,6*sizeof(int));
    this->DataDescription = vtkStructuredData::GetDataDescriptionFromExtent(ext);
    std::fill(this->Grow,this->Grow+3,0);
    std::fill(this->Implicit,this->Implicit+3,0);

607
    this->Nodes = nullptr;
608 609 610 611 612 613 614 615 616 617 618 619 620 621 622

    // Effectively, shallow copy the coordinate arrays and maintain ownership
    // of these arrays in the caller.
    this->X_Coords = vtkDataArray::CreateDataArray(x_coords->GetDataType());
    this->X_Coords->SetVoidArray(
        x_coords->GetVoidPointer(0),x_coords->GetNumberOfTuples(),1);

    this->Y_Coords = vtkDataArray::CreateDataArray(y_coords->GetDataType());
    this->Y_Coords->SetVoidArray(
        y_coords->GetVoidPointer(0),y_coords->GetNumberOfTuples(),1);

    this->Z_Coords = vtkDataArray::CreateDataArray(z_coords->GetDataType());
    this->Z_Coords->SetVoidArray(
        z_coords->GetVoidPointer(0),z_coords->GetNumberOfTuples(),1);

623
    if(fields != nullptr)
624
    {
625 626
      this->PointData = vtkPointData::New();
      this->PointData->ShallowCopy(fields);
627
    }
628
    else
629
    {
630
      this->PointData = nullptr;
631
    }
632 633
  }

634 635 636 637 638 639 640 641 642
//------------------------------------------------------------------------------
  void Initialize(int id, int ext[6], vtkPoints* nodes, vtkPointData* fields)
  {
  this->ID = id;
  memcpy(this->Extent,ext,6*sizeof(int));
  this->DataDescription = vtkStructuredData::GetDataDescriptionFromExtent(ext);
  std::fill(this->Grow,this->Grow+3,0);
  std::fill(this->Implicit,this->Implicit+3,0);

643 644 645
  this->X_Coords = nullptr;
  this->Y_Coords = nullptr;
  this->Z_Coords = nullptr;
646

647
  if(nodes != nullptr)
648
  {
649 650
    this->Nodes = vtkPoints::New();
    this->Nodes->ShallowCopy(nodes);
651
  }
652
  else
653
  {
654
    this->Nodes = nullptr;
655
  }
656

657
  if(fields != nullptr)
658
  {
659 660
    this->PointData = vtkPointData::New();
    this->PointData->ShallowCopy(fields);
661
  }
662
  else
663
  {
664
    this->PointData = nullptr;
665
  }
666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709
  }

};

//------------------------------------------------------------------------------
//  CommManager class Definition
//------------------------------------------------------------------------------

class CommunicationManager
{
public:
  CommunicationManager() {};
  ~CommunicationManager() { this->Clear(); };

  unsigned char* GetRcvBuffer(const int fromRank);
  unsigned int GetRcvBufferSize(const int fromRank);

  void EnqueueRcv(const int fromRank);
  void EnqueueSend(const int toRank, unsigned char* data, unsigned int nbytes);
  void Exchange(vtkMPIController* comm);
  int NumMsgs();
  void Clear();

private:
  // map send/rcv buffers based on rank.
  std::map< int, unsigned char* > Send;
  std::map< int, int> SendByteSize;
  std::map< int, unsigned char* > Rcv;
  std::map< int, int> RcvByteSize;
  std::vector< vtkMPICommunicator::Request > Requests;

  // exchanges buffer-sizes
  void AllocateRcvBuffers(vtkMPIController* comm);
};

//------------------------------------------------------------------------------
void CommunicationManager::Clear()
{
  this->Requests.clear();
  this->SendByteSize.clear();
  this->RcvByteSize.clear();

  std::map<int,unsigned char*>::iterator it;
  for(it=this->Send.begin(); it != this->Send.end(); ++it)
710
  {
711
    delete [] it->second;
712
  }
713 714 715
  this->Send.clear();

  for(it=this->Rcv.begin(); it != this->Rcv.end(); ++it)
716
  {
717
    delete [] it->second;
718
  }
719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740
  this->Rcv.clear();
}

//------------------------------------------------------------------------------
unsigned char* CommunicationManager::GetRcvBuffer(const int fromRank)
{
  assert( "pre: cannot find buffer for requested rank!" &&
          (this->Rcv.find( fromRank ) != this->Rcv.end()) );
  return( this->Rcv[ fromRank ] );
}

//------------------------------------------------------------------------------
unsigned int CommunicationManager::GetRcvBufferSize(const int fromRank)
{
  assert( "pre: cannot find bytesize size of requested rank!" &&
          (this->RcvByteSize.find( fromRank ) != this->RcvByteSize.end()) );
  return( this->RcvByteSize[ fromRank ]  );
}

//------------------------------------------------------------------------------
int CommunicationManager::NumMsgs()
{
741
  return static_cast<int>( this->Send.size()+this->Rcv.size() );
742 743 744 745 746 747 748 749
}

//------------------------------------------------------------------------------
void CommunicationManager::EnqueueRcv(const int fromRank)
{
  assert("pre: rcv from rank has already been enqueued!" &&
          (this->Rcv.find(fromRank)==this->Rcv.end()) );

750
  this->Rcv[ fromRank ] = nullptr;
751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775
  this->RcvByteSize[ fromRank ] = 0;
}

//------------------------------------------------------------------------------
void CommunicationManager::EnqueueSend(
      const int toRank, unsigned char* data, unsigned int nbytes)
{
  assert("pre: send to rank has already been enqueued!" &&
          (this->Send.find(toRank)==this->Send.end()));

  this->Send[ toRank ] = data;
  this->SendByteSize[ toRank ] = nbytes;
}

//------------------------------------------------------------------------------
void CommunicationManager::AllocateRcvBuffers(vtkMPIController* comm)
{
  std::map<int,int>::iterator it;

  // STEP 0: Allocate vector to store request objects for non-blocking comm.
  int rqstIdx = 0;
  this->Requests.resize(this->NumMsgs());

  // STEP 1: Post receives
  for(it=this->RcvByteSize.begin(); it != this->RcvByteSize.end(); ++it)
776
  {
777 778 779 780
    int fromRank = it->first;
    int* dataPtr = &(it->second);
    comm->NoBlockReceive(dataPtr,1,fromRank,0,this->Requests[rqstIdx]);
    ++rqstIdx;
781
  }
782 783 784

  // STEP 2: Post Sends
  for(it=this->SendByteSize.begin(); it != this->SendByteSize.end(); ++it)
785
  {
786 787 788 789
    int toRank   = it->first;
    int* dataPtr = &(it->second);
    comm->NoBlockSend(dataPtr,1,toRank,0,this->Requests[rqstIdx]);
    ++rqstIdx;
790
  }
791 792

  // STEP 3: WaitAll
793
  if (!this->Requests.empty())
794
  {
795
    comm->WaitAll(this->NumMsgs(),&this->Requests[0]);
796
  }
797 798 799 800 801
  this->Requests.clear();

  // STEP 4: Allocate rcv buffers
  std::map<int,unsigned char*>::iterator bufferIter = this->Rcv.begin();
  for( ;bufferIter != this->Rcv.end(); ++bufferIter)
802
  {
803
    int fromRank = bufferIter->first;
804
    assert("pre: rcv buffer should be nullptr!" && (this->Rcv[fromRank]==nullptr) );
805
    this->Rcv[ fromRank ] = new unsigned char[ this->RcvByteSize[fromRank] ];
806
  }
807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822
}

//------------------------------------------------------------------------------
void CommunicationManager::Exchange(vtkMPIController* comm)
{
  std::map<int,unsigned char*>::iterator it;

  // STEP 0: exchange & allocate buffer sizes
  this->AllocateRcvBuffers(comm);

  // STEP 1: Allocate vector to store request objects for non-blocking comm.
  int rqstIdx = 0;
  this->Requests.resize(this->NumMsgs());

  // STEP 2: Post Rcvs
  for(it=this->Rcv.begin(); it != this->Rcv.end(); ++it)
823
  {
824 825 826 827 828 829 830 831
    int fromRank          = it->first;
    unsigned char* buffer = it->second;
    assert("pre: rcv buffer size not found!" &&
            this->RcvByteSize.find(fromRank) != this->RcvByteSize.end() );
    unsigned int bytesize = this->RcvByteSize[ fromRank ];

    comm->NoBlockReceive(buffer,bytesize,fromRank,0,this->Requests[rqstIdx]);
    ++rqstIdx;
832
  }
833 834 835

  // STEP 3: Post Sends
  for(it=this->Send.begin(); it != this->Send.end(); ++it)
836
  {
837 838 839 840 841 842 843 844
    int toRank            = it->first;
    unsigned char* buffer = it->second;
    assert("pre: rcv buffer size not found!" &&
            this->SendByteSize.find(toRank) != this->SendByteSize.end() );
    unsigned int bytesize = this->SendByteSize[ toRank ];

    comm->NoBlockSend(buffer,bytesize,toRank,0,this->Requests[rqstIdx]);
    ++rqstIdx;
845
  }
846 847

  // STEP 4: WaitAll
848
  if (!this->Requests.empty())
849
  {
850
    comm->WaitAll(this->NumMsgs(),&this->Requests[0]);
851
  }
852 853 854 855 856 857 858 859 860 861 862 863 864 865
  this->Requests.clear();
}

} // END namespace detail
} // END namespace vtk
//==============================================================================
// END INTERNAL DATASTRUCTURE DEFINITIONS
//==============================================================================

vtkStandardNewMacro(vtkStructuredImplicitConnectivity);

//------------------------------------------------------------------------------
vtkStructuredImplicitConnectivity::vtkStructuredImplicitConnectivity()
{
866 867 868 869
  this->DomainInfo  = nullptr;
  this->InputGrid   = nullptr;
  this->OutputGrid  = nullptr;
  this->CommManager = nullptr;
870 871 872 873 874 875 876
  this->Controller  = vtkMPIController::SafeDownCast(
            vtkMultiProcessController::GetGlobalController());
}

//------------------------------------------------------------------------------
vtkStructuredImplicitConnectivity::~vtkStructuredImplicitConnectivity()
{
877
  delete this->DomainInfo;
878
  this->DomainInfo = nullptr;
879

880
  if( this->InputGrid != nullptr )
881
  {
882
    this->InputGrid->Clear();
883
    delete this->InputGrid;
884
    this->InputGrid = nullptr;
885
  }
886

887
  if( this->OutputGrid != nullptr )
888
  {
889
    this->OutputGrid->Clear();
890
    delete this->OutputGrid;
891
    this->OutputGrid = nullptr;
892
  }
893

894
  if( this->CommManager != nullptr )
895
  {
896 897
    this->CommManager->Clear();
    delete this->CommManager;
898
    this->CommManager = nullptr;
899
  }
900

901
  this->Controller = nullptr;
902 903 904 905 906 907 908
}

//------------------------------------------------------------------------------
void vtkStructuredImplicitConnectivity::PrintSelf(ostream& os,vtkIndent indent)
{
  this->Superclass::PrintSelf(os,indent);
  os << "Controller: "      << this->Controller << std::endl;
909
  if( this->Controller != nullptr )
910
  {
911 912
    os << "Number of Ranks: " << this->Controller->GetNumberOfProcesses();
    os << std::endl;
913
  } // END if Controller != nullptr
914 915

  os << "Input Grid: " << this->InputGrid  << std::endl;
916
  if( this->InputGrid != nullptr )
917
  {
918 919 920 921 922 923 924 925 926 927 928 929 930 931
    os << "Extent: [" << this->InputGrid->Extent[0];
    os << ", " << this->InputGrid->Extent[1];
    os << ", " << this->InputGrid->Extent[2];
    os << ", " << this->InputGrid->Extent[3];
    os << ", " << this->InputGrid->Extent[4];
    os << ", " << this->InputGrid->Extent[5];
    os << "] " << std::endl;

    os << "Grow: [" << this->InputGrid->Grow[0];
    os << ", " << this->InputGrid->Grow[1];
    os << ", " << this->InputGrid->Grow[2];
    os << "] " << std::endl;

    os << "Number of Neighbors: " << this->InputGrid->Neighbors.size();
932
    os << std::endl;
933 934
    size_t N = this->InputGrid->Neighbors.size();
    for(size_t nei=0; nei < N; ++nei)
935
    {
936 937
      os << "\t" << this->InputGrid->Neighbors[ nei ].ToString();
      os << std::endl;
938
    } // END for all neighbors
939
  } // END if InputGrid != nullptr
940 941 942 943 944
}

//------------------------------------------------------------------------------
void vtkStructuredImplicitConnectivity::SetWholeExtent(int wholeExt[6])
{
945
  delete this->DomainInfo;
946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961

  this->DomainInfo = new vtk::detail::DomainMetaData();
  this->DomainInfo->Initialize(wholeExt);

  assert("post: Domain description does not match across ranks!" &&
          this->GlobalDataDescriptionMatch() );
}

//------------------------------------------------------------------------------
void vtkStructuredImplicitConnectivity::RegisterGrid(
      const int gridID,
      int extent[6],
      vtkPoints* gridNodes,
      vtkPointData* pointData)
{
  // Sanity Checks!
962 963 964 965
  assert("pre: nullptr Domain, whole extent is not set!" &&
          (this->DomainInfo != nullptr) );
  assert("pre: input not nullptr in this process!" &&
          (this->InputGrid == nullptr) );
966 967
  assert("pre: input grid ID should be >= 0" && (gridID >= 0) );

968
  delete this->InputGrid;
969
  this->InputGrid = nullptr;
970

971 972 973
  // Only add if the grid falls within the output extent. Processes that do
  // not contain the VOI will fail this test.
  if (this->DomainInfo->HasGrid(extent))
974
  {
975 976
    this->InputGrid = new vtk::detail::StructuredGrid();
    this->InputGrid->Initialize(gridID,extent,gridNodes,pointData);
977
  }
978 979
}

980 981 982 983 984 985 986 987 988 989
//------------------------------------------------------------------------------
void vtkStructuredImplicitConnectivity::RegisterRectilinearGrid(
            const int gridID,
            int extent[6],
            vtkDataArray* xcoords,
            vtkDataArray* ycoords,
            vtkDataArray* zcoords,
            vtkPointData* pointData)
{
  // Sanity Checks!
990 991 992 993
  assert("pre: nullptr Domain, whole extent is not set!" &&
          (this->DomainInfo != nullptr) );
  assert("pre: input not nullptr in this process!" &&
          (this->InputGrid == nullptr) );
994 995
  assert("pre: input grid ID should be >= 0" && (gridID >= 0) );

996
  delete this->InputGrid;
997
  this->InputGrid = nullptr;
998

999 1000 1001
  // Only add if the grid falls within the output extent. Processes that do
  // not contain the VOI will fail this test.
  if (this->DomainInfo->HasGrid(extent))
1002
  {
1003 1004 1005
    this->InputGrid = new vtk::detail::StructuredGrid();
    this->InputGrid->Initialize(gridID, extent, xcoords, ycoords, zcoords,
                                pointData);
1006
  }
1007 1008
}

1009 1010 1011 1012
//------------------------------------------------------------------------------
void vtkStructuredImplicitConnectivity::ExchangeExtents()
{
  // Sanity checks!
1013 1014
  assert("pre: null controller!" && (this->Controller != nullptr) );
  assert("pre: null domain!" && (this->DomainInfo != nullptr) );
1015 1016 1017 1018

  // STEP 0: Construct the extent buffer that will be sent from each process.
  // Each process sends 7 ints: [gridId imin imax jmin jmax kmin kmax]
  int extbuffer[7];
1019
  if( this->InputGrid == nullptr )
1020
  {
1021 1022
    // pad the buffer with -1, indicating that this process has no grid
    std::fill(extbuffer,extbuffer+7,-1);
1023
  }
1024
  else
1025
  {
1026 1027
    extbuffer[0] = this->InputGrid->ID;
    memcpy(&extbuffer[1],this->InputGrid->Extent,6*sizeof(int));
1028
  }
1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041

  // STEP 1: Allocate receive buffer, we receive 7 ints for each rank
  int nranks = this->Controller->GetNumberOfProcesses();
  this->DomainInfo->ExtentListInfo.resize(7*nranks,0);

  // STEP 2: AllGather
  int* rcvbuffer = &(this->DomainInfo->ExtentListInfo)[0];
  this->Controller->AllGather(extbuffer,rcvbuffer,7);
}

//------------------------------------------------------------------------------
void vtkStructuredImplicitConnectivity::ComputeNeighbors()
{
1042
  if (!this->InputGrid)
1043
  {
1044
    return;
1045
  }
1046

1047 1048 1049 1050 1051 1052 1053 1054
  int type;
  vtk::detail::Interval A; // used to store the local interval at each dim
  vtk::detail::Interval B; // used to store the remote interval at each dim
  vtk::detail::Interval Overlap; // used to store the computed overlap
  vtk::detail::ImplicitNeighbor Neighbor; // used to store neighbor information

  int nranks = this->Controller->GetNumberOfProcesses();
  for(int rank=0; rank < nranks; ++rank)
1055
  {
1056 1057
    int rmtID = this->DomainInfo->ExtentListInfo[rank*7];
    if( (rmtID == this->InputGrid->ID) || (rmtID == -1) )
1058
    {
1059 1060
      // skip self or empty remote grid
      continue;
1061
    }
1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075

    int* rmtExtent = &(this->DomainInfo->ExtentListInfo)[rank*7+1];

    // Initialize neighbor data-structure
    Neighbor.Rank = rank;
    memcpy(Neighbor.Extent,rmtExtent,6*sizeof(int));
    memcpy(Neighbor.Overlap,Neighbor.Extent,6*sizeof(int));
    std::fill(Neighbor.Orientation,Neighbor.Orientation+3,
        vtk::detail::IntervalsConnect::UNDEFINED);

    bool disregard = false;
    int nimplicit  = 0;

    for(int dim=0; dim < this->DomainInfo->NDim; ++dim)
1076
    {
1077 1078 1079 1080 1081 1082 1083
      int d = this->DomainInfo->DimIndex[dim];
      assert("pre: invalid dimension!" && (d >= 0) && (d <= 2) );

      A.Set( this->InputGrid->Extent[d*2], this->InputGrid->Extent[d*2+1] );
      B.Set( rmtExtent[d*2], rmtExtent[d*2+1] );

      if( A.ImplicitNeighbor(B,type) )
1084
      {
1085 1086 1087 1088 1089 1090
        this->InputGrid->Implicit[ d ]= 1;
        Neighbor.Orientation[ d ] = type;

        // Compute overlap based on the fact that we are communicating
        // data to the left <=> grow to the right.
        if( type == vtk::detail::IntervalsConnect::IMPLICIT_HI )
1091
        {
1092 1093 1094 1095
          ++nimplicit;
          Neighbor.Overlap[d*2]    =
          Neighbor.Overlap[d*2+1]  = Neighbor.Extent[d*2];
          this->InputGrid->Grow[d] = 1; /* increment by 1 in this dimension */
1096
        } // END if IMPLICIT_HI
1097
        else if( type == vtk::detail::IntervalsConnect::IMPLICIT_LO )
1098
        {
1099 1100 1101
          ++nimplicit;
          Neighbor.Overlap[d*2]   =
          Neighbor.Overlap[d*2+1] = this->InputGrid->Extent[d*2];
1102
        } // END else if IMPLICIT_LO
1103
        else
1104
        {
1105 1106 1107 1108
          vtkGenericWarningMacro(
              << "Invalid implicit connectivity type! "
              << "Code should not reach here!\n"
              );
1109 1110
        } // END else
      } // END if implicit
1111
      else if( A.Intersects(B,Overlap,type) )
1112
      {
1113 1114 1115
        Neighbor.Orientation[ d ] = type;
        Neighbor.Overlap[d*2]     = Overlap.Low();
        Neighbor.Overlap[d*2+1]   = Overlap.High();
1116
      } // END if intersect
1117
      else
1118
      {
1119 1120
        disregard = true;
        Neighbor.Orientation[ d ] = type;
1121 1122
      } // END else
    } // END for all dimensions
1123 1124 1125 1126 1127

    // Determine whether to include the neighbor to the list of neighbors in
    // this rank.

    if( !(nimplicit > 1 || disregard) )
1128
    {
1129
      this->InputGrid->Neighbors.push_back( Neighbor );
1130
    }
1131

1132
  } // END for all ranks
1133 1134 1135 1136 1137 1138 1139 1140 1141 1142
}

//------------------------------------------------------------------------------
bool vtkStructuredImplicitConnectivity::GlobalDataDescriptionMatch()
{
  int sum = -1;
  this->Controller->AllReduce(
      &this->DomainInfo->DataDescription,&sum,1,vtkCommunicator::SUM_OP);
  if( (sum/this->Controller->GetNumberOfProcesses()) ==
      this->DomainInfo->DataDescription)
1143
  {
1144
    return true;
1145
  }
1146 1147 1148 1149 1150 1151
  return false;
}

//------------------------------------------------------------------------------
bool vtkStructuredImplicitConnectivity::HasImplicitConnectivity()
{
1152
  if( this->DomainInfo == nullptr )
1153
  {
1154
    vtkGenericWarningMacro(<< "nullptr domain, WholeExtent not set!");
1155
    return false;
1156
  }
1157 1158 1159 1160

  if( (this->DomainInfo->GlobalImplicit[0] > 0) ||
      (this->DomainInfo->GlobalImplicit[1] > 0) ||
      (this->DomainInfo->GlobalImplicit[2] > 0) )
1161
  {
1162
    return true;
1163
  }
1164 1165 1166 1167 1168 1169 1170
  return false;
}

//------------------------------------------------------------------------------
void vtkStructuredImplicitConnectivity::GetGlobalImplicitConnectivityState()
{
  // Sanity checks!
1171
  assert("pre: null controller!" && (this->Controller != nullptr) );
1172 1173

  int sndbuffer[3];
1174
  if( this->InputGrid == nullptr)
1175
  {
1176
    std::fill(sndbuffer,sndbuffer+3,0);
1177
  }
1178
  else
1179
  {