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

  Program:   Visualization Toolkit
  Module:    vtkPUnstructuredGridGhostCellsGenerator.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 "vtkPUnstructuredGridGhostCellsGenerator.h"

18
#include "vtkAppendFilter.h"
19
#include "vtkCellArray.h"
20
#include "vtkCellData.h"
21 22 23 24 25 26 27
#include "vtkCharArray.h"
#include "vtkDataSetSurfaceFilter.h"
#include "vtkExtractCells.h"
#include "vtkIdList.h"
#include "vtkIdTypeArray.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
28
#include "vtkKdTree.h"
29 30
#include "vtkMPICommunicator.h"
#include "vtkMPIController.h"
31
#include "vtkMath.h"
32 33 34 35 36 37 38
#include "vtkMergeCells.h"
#include "vtkMergePoints.h"
#include "vtkNew.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkPoints.h"
#include "vtkStreamingDemandDrivenPipeline.h"
39
#include "vtkUnsignedCharArray.h"
40 41
#include "vtkUnstructuredGrid.h"

42
#include <algorithm>
43 44
#include <map>
#include <set>
45
#include <vector>
46 47 48 49 50 51 52 53 54

//----------------------------------------------------------------------------
// Helpers
namespace
{
// Class to hold asynchronous communication information
class CommDataInfo
{
public:
55 56 57 58 59 60 61 62
  CommDataInfo()
    : SendLen(-1)
    , RecvLen(-1)
    , CommStep(0)
  {
    this->SendBuffer = vtkCharArray::New();
    this->RecvBuffer = vtkCharArray::New();
  }
63 64

  CommDataInfo(const CommDataInfo& c)
65 66 67 68 69 70 71
  {
    *this = c;
    if (this->SendBuffer)
    {
      this->SendBuffer->Register(nullptr);
    }
    if (this->RecvBuffer)
72
    {
73
      this->RecvBuffer->Register(nullptr);
74
    }
75
  }
76 77

  ~CommDataInfo()
78 79
  {
    if (this->SendBuffer)
80
    {
81
      this->SendBuffer->Delete();
82
    }
83 84 85 86 87
    if (this->RecvBuffer)
    {
      this->RecvBuffer->Delete();
    }
  }
88 89 90

  vtkMPICommunicator::Request SendReqs[2];
  vtkMPICommunicator::Request RecvReqs[2];
91 92
  vtkCharArray* SendBuffer;
  vtkCharArray* RecvBuffer;
93 94 95
  vtkIdType SendLen;
  vtkIdType RecvLen;
  int CommStep;
96
  int RecvSize;
97
};
98
} // end anonymous namespace
99 100 101

struct vtkPUnstructuredGridGhostCellsGenerator::vtkInternals
{
102
  // SubController only has MPI processes which have cells
103
  vtkMPIController* SubController;
104

105 106
  // For global ids
  std::map<vtkIdType, vtkIdType> GlobalToLocalPointIdMap;
107
  std::map<int, std::vector<vtkIdType> > ProcessIdToSurfacePointIds;
108 109 110
  // Ids to send to a specific process. Only the ids of points in the
  // receive process's bounding box are sent.
  std::map<int, std::vector<vtkIdType> > SendIds;
111 112

  // For point coordinates
113
  std::map<int, std::vector<double> > ProcessIdToSurfacePoints;
114
  vtkSmartPointer<vtkIdTypeArray> LocalPointsMap; // from surface id to 3d grid id
115 116 117 118 119
  // Points to send to a specific process. Only the points in the
  // receive process's bounding box are sent.
  std::map<int, std::vector<double> > SendPoints;
  vtkSmartPointer<vtkDataArray> MyPoints;

120 121
  std::map<int, CommDataInfo> CommData;
  vtkUnstructuredGridBase* Input;
122
  vtkSmartPointer<vtkUnstructuredGrid> CurrentGrid;
123

124
  vtkIdTypeArray* InputGlobalPointIds;
125
  bool UseGlobalPointIds;
126 127 128 129

  // cells that need to be sent to a given proc
  std::map<int, std::set<vtkIdType> > CellsToSend;

130 131
  // cells that have been sent to a given proc over the entire time.
  // used to make sure we only send a cell once to a destination process.
132 133
  std::map<int, std::set<vtkIdType> > SentCells;

134
  // cells that have been received from a given proc over the entire time.
135 136
  // stores global cell id. used this to make sure that we don't send
  // a cell back to a process that already sent it to this rank
137 138 139 140 141 142
  std::map<int, std::set<vtkIdType> > ReceivedCells;

  // mapping from global cell id to local cell id.
  // only stores cells which have been received (aka are ghost cells)
  std::map<vtkIdType, vtkIdType> GlobalToLocalCellIdMap;

143 144 145
  // cells that were sent to a proc during the last round,
  // a "round" is receiving one layer of ghost cells
  std::map<int, std::set<vtkIdType> > SentCellsLastRound;
146 147 148 149 150

  // list of processes which are probably my neighbors. this
  // is based on overlapping local bounding boxes so it is
  // not guaranteed that they really are sharing an interprocess boundary
  std::vector<int> Neighbors;
151 152
};

153 154
namespace
{
155 156
static const int UGGCG_SIZE_EXCHANGE_TAG = 9000;
static const int UGGCG_DATA_EXCHANGE_TAG = 9001;
157
}
158 159 160

//----------------------------------------------------------------------------

161
vtkStandardNewMacro(vtkPUnstructuredGridGhostCellsGenerator)
162 163 164 165 166 167
vtkSetObjectImplementationMacro(
  vtkPUnstructuredGridGhostCellsGenerator, Controller, vtkMultiProcessController);

//----------------------------------------------------------------------------
vtkPUnstructuredGridGhostCellsGenerator::vtkPUnstructuredGridGhostCellsGenerator()
{
168
  this->Controller = nullptr;
169 170
  this->SetController(vtkMultiProcessController::GetGlobalController());

171
  this->Internals = nullptr;
172 173 174 175 176
}

//----------------------------------------------------------------------------
vtkPUnstructuredGridGhostCellsGenerator::~vtkPUnstructuredGridGhostCellsGenerator()
{
177
  this->SetController(nullptr);
178 179

  delete this->Internals;
180
  this->Internals = nullptr;
181 182 183 184 185 186
}

//-----------------------------------------------------------------------------
void vtkPUnstructuredGridGhostCellsGenerator::PrintSelf(ostream& os, vtkIndent indent)
{
  Superclass::PrintSelf(os, indent);
187 188
}

189
//-----------------------------------------------------------------------------
190 191
int vtkPUnstructuredGridGhostCellsGenerator::RequestData(vtkInformation* vtkNotUsed(request),
  vtkInformationVector** inputVector, vtkInformationVector* outputVector)
192 193
{
  // get the info objects
194 195
  vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation* outInfo = outputVector->GetInformationObject(0);
196 197 198

  // get the input and output. Input may just have the UnstructuredGridBase
  // interface, but output should be an unstructured grid.
199 200 201 202
  vtkUnstructuredGridBase* input =
    vtkUnstructuredGridBase::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkUnstructuredGrid* output =
    vtkUnstructuredGrid::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));
203

204
  if (!input)
205
  {
206
    vtkErrorMacro("No input data!");
207
    return 0;
208
  }
209

210
  if (!this->Controller)
211
  {
212
    this->Controller = vtkMultiProcessController::GetGlobalController();
213
  }
214

215 216 217 218
  int reqGhostLevel =
    outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS());
  int maxGhostLevel = this->BuildIfRequired ? reqGhostLevel : std::max(reqGhostLevel,
                                                                this->MinimumNumberOfGhostLevels);
219

220
  if (maxGhostLevel == 0 || this->Controller->GetNumberOfProcesses() == 1)
221
  {
222 223 224 225 226 227 228 229 230 231 232 233 234 235 236
    vtkDebugMacro("Don't need ghost cells or only have a single process. Nothing more to do.");
    output->ShallowCopy(input);
    return 1;
  }

  // if only a single process has cells then we can skip ghost cell computations but
  // otherwise we need to do it from scratch since the ghost information coming in
  // may be wrong (it was for the vtkFiltersParallelCxx-MPI-ParallelConnectivity4 test)
  int needsGhosts = input->GetNumberOfCells() > 0 ? 1 : 0;

  int globalNeedsGhosts = 0;
  this->Controller->AllReduce(&needsGhosts, &globalNeedsGhosts, 1, vtkCommunicator::SUM_OP);
  if (globalNeedsGhosts < 2)
  {
    vtkDebugMacro("At most one process has cells. Nothing more to do.");
237 238 239 240
    output->ShallowCopy(input);
    return 1;
  }

241
  // determine which processes have any non-ghost cells and then create a subcontroller
242 243
  // for just them to use
  int hasCells = input->GetNumberOfCells() > 0 ? 1 : 0;
244 245 246 247 248
  if (hasCells && input->GetCellGhostArray() && input->GetCellGhostArray()->GetRange()[0] != 0)
  {
    hasCells = 0; // all the cells are ghost cells which we don't care about anymore
  }

249 250 251 252 253 254 255
  vtkSmartPointer<vtkMPIController> subController;
  subController.TakeReference(
    vtkMPIController::SafeDownCast(this->Controller)->PartitionController(hasCells, 0));

  if (hasCells == 0 || subController->GetNumberOfProcesses() < 2)
  {
    vtkDebugMacro("No work to do since at most one process has data");
256 257
    output->ShallowCopy(input);
    return 1;
258
  }
259

260 261
  vtkNew<vtkUnstructuredGrid> cleanedInput;
  vtkUnsignedCharArray* cellGhostArray = input->GetCellGhostArray();
262
  if (cellGhostArray == nullptr || cellGhostArray->GetValueRange()[1] == 0)
263 264 265 266 267 268 269 270 271 272 273 274 275 276
  {
    // we either have no ghost cells or do but there are no ghost entities so we just need
    // to remove those arrays and can skip modifying the data set itself
    cleanedInput->ShallowCopy(input);
  }
  else
  {
    cleanedInput->DeepCopy(input);
    cleanedInput->RemoveGhostCells();
  }
  cleanedInput->GetPointData()->RemoveArray(vtkDataSetAttributes::GhostArrayName());
  cleanedInput->GetCellData()->RemoveArray(vtkDataSetAttributes::GhostArrayName());
  input = nullptr; // nullify input to make sure we don't use it after this

277 278
  delete this->Internals;
  this->Internals = new vtkPUnstructuredGridGhostCellsGenerator::vtkInternals();
279 280
  this->Internals->SubController = subController;

281
  this->Internals->Input = cleanedInput;
282

283
  vtkPointData* inputPD = cleanedInput->GetPointData();
284
  this->Internals->InputGlobalPointIds = vtkIdTypeArray::FastDownCast(inputPD->GetGlobalIds());
285 286

  if (!this->Internals->InputGlobalPointIds)
287
  {
288
    inputPD = cleanedInput->GetPointData();
289 290
    this->Internals->InputGlobalPointIds =
      vtkIdTypeArray::FastDownCast(inputPD->GetArray(this->GlobalPointIdsArrayName));
291
    inputPD->SetGlobalIds(this->Internals->InputGlobalPointIds);
292
  }
293 294

  if (!this->UseGlobalPointIds)
295
  {
296
    this->Internals->InputGlobalPointIds = nullptr;
297
  }
298
  else
299
  {
300 301
    int useGlobalPointIds = this->Internals->InputGlobalPointIds != nullptr ? 1 : 0;
    int allUseGlobalPointIds;
302 303
    this->Internals->SubController->AllReduce(
      &useGlobalPointIds, &allUseGlobalPointIds, 1, vtkCommunicator::MIN_OP);
304 305 306 307
    if (!allUseGlobalPointIds)
    {
      this->Internals->InputGlobalPointIds = nullptr;
    }
308
  }
309

310 311 312
  // ensure that global cell ids array is there if specified.
  // only need global cell ids when more than one ghost layer is needed
  if (maxGhostLevel > 1)
313
  {
314
    if (this->HasGlobalCellIds)
315
    {
316
      vtkCellData* inputCD = cleanedInput->GetCellData();
317
      if (!inputCD->GetGlobalIds())
318
      {
319
        vtkDataArray* globalCellIdsArray = inputCD->GetArray(this->GlobalCellIdsArrayName);
320
        if (globalCellIdsArray == nullptr)
321
        {
322
          this->SetHasGlobalCellIds(false);
323
        }
324
        else
325
        {
326 327
          inputCD->SetGlobalIds(globalCellIdsArray);
        }
328
      }
329
      else
330
      {
331 332 333
        // make sure GlobalCellIdsArrayName is correct
        this->SetGlobalCellIdsArrayName(inputCD->GetGlobalIds()->GetName());
      }
334
    }
335 336 337 338

    // ensure that everyone has the same value of HasGlobalCellIds
    int hasGlobalCellIds = this->HasGlobalCellIds != 0 ? 1 : 0;
    int allHasGlobalCellIds;
339 340
    this->Internals->SubController->AllReduce(
      &hasGlobalCellIds, &allHasGlobalCellIds, 1, vtkCommunicator::MIN_OP);
341
    if (!allHasGlobalCellIds)
342
    {
343 344
      this->HasGlobalCellIds = false;
    }
345
  }
346 347 348

  // add global cell ids if necessary
  if (!this->HasGlobalCellIds && maxGhostLevel > 1)
349
  {
350
    this->AddGlobalCellIds();
351
  }
352 353

  // obtain first level of ghost cells
354 355
  this->Internals->CurrentGrid = vtkSmartPointer<vtkUnstructuredGrid>::New();
  this->GetFirstGhostLayer(maxGhostLevel, this->Internals->CurrentGrid);
356

357
  // add additional ghost layers one at a time
358
  for (int i = 1; i < maxGhostLevel; i++)
359
  {
360
    this->AddGhostLayer(i + 1, maxGhostLevel);
361
  }
362

363 364
  // remove global cell ids if they were added internally
  if (!this->HasGlobalCellIds && maxGhostLevel > 1)
365
  {
366
    this->RemoveGlobalCellIds();
367
  }
368

369
  // set the output
370
  output->ShallowCopy(this->Internals->CurrentGrid);
371
  output->GetInformation()->Set(vtkDataObject::DATA_NUMBER_OF_GHOST_LEVELS(), maxGhostLevel);
372 373

  delete this->Internals;
374
  this->Internals = nullptr;
375

376
  vtkDebugMacro("Produced " << maxGhostLevel << " ghost levels.");
377 378 379
  return 1;
}

380
//-----------------------------------------------------------------------------
381
// Get the first layer of ghost cells
382
void vtkPUnstructuredGridGhostCellsGenerator::GetFirstGhostLayer(
383
  int maxGhostLevel, vtkUnstructuredGrid* output)
384
{
385 386 387 388
  std::vector<double> allBounds;
  this->ExchangeBoundsAndDetermineNeighbors(allBounds);
  this->ExtractAndReduceSurfacePointsShareData(allBounds);
  allBounds.clear();
389 390
  this->UpdateProgress(1.0 / (3.0 * maxGhostLevel));
  this->ComputeSharedPoints();
391

392 393 394 395 396 397 398 399 400 401
  this->UpdateProgress(2.0 / (3.0 * maxGhostLevel));

  this->ExtractAndSendGhostCells(this->Internals->Input);
  this->UpdateProgress(2.5 / (3.0 * maxGhostLevel));

  // Shallow copy the input grid and initialize the ghost arrays
  vtkNew<vtkUnstructuredGrid> inputCopy;
  inputCopy->ShallowCopy(this->Internals->Input);
  inputCopy->AllocatePointGhostArray();
  inputCopy->AllocateCellGhostArray();
402 403 404 405 406 407 408
  this->ReceiveAndMergeGhostCells(1, maxGhostLevel, inputCopy.Get(), output);

  this->UpdateProgress(1.0 / maxGhostLevel);
}

//-----------------------------------------------------------------------------
// Step 0: Exchange bounds, and determine your neighbors
409 410
void vtkPUnstructuredGridGhostCellsGenerator::ExchangeBoundsAndDetermineNeighbors(
  std::vector<double>& allBounds)
411 412 413 414 415 416 417 418
{
  // increase bounds by a certain percentage to deal with precision stuff
  double epsilon = 0.01;

  double bounds[6];
  this->Internals->Input->GetBounds(bounds);

  // everyone shares bounds
419
  allBounds.resize(this->Internals->SubController->GetNumberOfProcesses() * 6);
420 421 422 423 424 425
  this->Internals->SubController->AllGather(bounds, &allBounds[0], 6);

  double xlength = bounds[1] - bounds[0];
  double ylength = bounds[3] - bounds[2];
  double zlength = bounds[5] - bounds[4];

426 427 428 429 430 431
  double xmin = bounds[0] - xlength * epsilon;
  double xmax = bounds[1] + xlength * epsilon;
  double ymin = bounds[2] - ylength * epsilon;
  double ymax = bounds[3] + ylength * epsilon;
  double zmin = bounds[4] - zlength * epsilon;
  double zmax = bounds[5] + zlength * epsilon;
432 433 434 435

  // go through bounds, and find the ones which intersect my bounds,
  // which are my possible neighbors
  int rank = this->Internals->SubController->GetLocalProcessId();
436
  for (int p = 0; p < this->Internals->SubController->GetNumberOfProcesses(); p++)
437
  {
438
    if (p == rank)
439 440 441 442
    {
      continue;
    }

443 444 445
    double xlength2 = allBounds[p * 6 + 1] - allBounds[p * 6 + 0];
    double xmin2 = allBounds[p * 6 + 0] - xlength2 * epsilon;
    double xmax2 = allBounds[p * 6 + 1] + xlength2 * epsilon;
446

447
    if (xmin <= xmax2 && xmax >= xmin2)
448
    {
449 450 451 452
      double ylength2 = allBounds[p * 6 + 3] - allBounds[p * 6 + 2];
      double ymin2 = allBounds[p * 6 + 2] - ylength2 * epsilon;
      double ymax2 = allBounds[p * 6 + 3] + ylength2 * epsilon;
      if (ymin <= ymax2 && ymax >= ymin2)
453
      {
454 455 456 457
        double zlength2 = allBounds[p * 6 + 5] - allBounds[p * 6 + 4];
        double zmin2 = allBounds[p * 6 + 4] - zlength2 * epsilon;
        double zmax2 = allBounds[p * 6 + 5] + zlength2 * epsilon;
        if (zmin <= zmax2 && zmax >= zmin2)
458 459 460 461 462 463 464
        {
          // this proc is a neighbor
          this->Internals->Neighbors.push_back(p);
        }
      }
    }
  }
465 466
}

467
//-----------------------------------------------------------------------------
468 469
// Step 1a: Extract surface geometry and send to my neighbors. Receive my
// neighbor's surface points
470 471
void vtkPUnstructuredGridGhostCellsGenerator::ExtractAndReduceSurfacePointsShareData(
  std::vector<double>& allBounds)
472 473 474 475 476 477 478
{
  // Extract boundary cells and points with the surface filter
  vtkNew<vtkDataSetSurfaceFilter> surfaceFilter;
  surfaceFilter->SetInputData(this->Internals->Input);
  surfaceFilter->PassThroughPointIdsOn();
  surfaceFilter->Update();

479
  vtkPolyData* surface = surfaceFilter->GetOutput();
480
  vtkIdType nbSurfacePoints = surface->GetNumberOfPoints();
481 482 483

  double bounds[6];
  surface->GetBounds(bounds);
484 485
  double delta[3] = { .0001 * (bounds[1] - bounds[0]), .0001 * (bounds[3] - bounds[2]),
    .0001 * (bounds[5] - bounds[4]) };
486

487
  vtkIdTypeArray* surfaceOriginalPointIds = vtkArrayDownCast<vtkIdTypeArray>(
488 489
    surface->GetPointData()->GetArray(surfaceFilter->GetOriginalPointIdsName()));

490
  std::vector<vtkMPICommunicator::Request> sendReqs(this->Internals->Neighbors.size() * 2);
491 492

  // reset CommStep
493
  std::map<int, CommDataInfo>::iterator comIter = this->Internals->CommData.begin();
494 495 496 497 498
  for (; comIter != this->Internals->CommData.end(); ++comIter)
  {
    comIter->second.CommStep = 0;
  }

499 500
  // we need sizesToSend to stick around for the noblocksends
  std::vector<int> sizesToSend(this->Internals->Neighbors.size());
501
  if (this->Internals->InputGlobalPointIds)
502
  {
503 504 505
    // get all sizes from neighbors
    // first set up the receives
    std::vector<int>::iterator iter = this->Internals->Neighbors.begin();
506
    for (; iter != this->Internals->Neighbors.end(); ++iter)
507 508
    {
      CommDataInfo& c = this->Internals->CommData[*iter];
509 510
      this->Internals->SubController->NoBlockReceive(
        &c.RecvSize, 1, *iter, UGGCG_SIZE_EXCHANGE_TAG, c.RecvReqs[0]);
511 512 513 514
    }

    // store the global point id arrays unique to each process (based on bounding box
    // of the receiving process) to send
515
    this->Internals->ProcessIdToSurfacePointIds.clear();
516

517 518
    for (iter = this->Internals->Neighbors.begin(); iter != this->Internals->Neighbors.end();
         ++iter)
519
    {
520 521 522
      std::vector<vtkIdType>& sendIds = this->Internals->SendIds[*iter];
      sendIds.clear();
      for (vtkIdType i = 0; i < nbSurfacePoints; i++)
523
      {
524 525
        double coord[3];
        surface->GetPoint(i, coord);
526
        if (vtkMath::PointIsWithinBounds(coord, &allBounds[*iter * 6], delta))
527
        {
528 529
          vtkIdType origPtId = surfaceOriginalPointIds->GetValue(i);
          vtkIdType globalPtId = this->Internals->InputGlobalPointIds->GetTuple1(origPtId);
530
          this->Internals->GlobalToLocalPointIdMap[globalPtId] = origPtId;
531
          sendIds.push_back(globalPtId);
532 533
        }
      }
534
    }
535

536
    // send surface point ids to each neighbor
537
    int reqidx = 0;
538 539
    for (iter = this->Internals->Neighbors.begin(); iter != this->Internals->Neighbors.end();
         ++iter)
540
    {
541
      std::vector<vtkIdType>& sendIds = this->Internals->SendIds[*iter];
542
      // send size of vector
543
      sizesToSend[reqidx] = static_cast<int>(sendIds.size());
544 545
      this->Internals->SubController->NoBlockSend(
        &sizesToSend[reqidx], 1, *iter, UGGCG_SIZE_EXCHANGE_TAG, sendReqs[2 * reqidx]);
546 547

      // send the vector
548 549
      this->Internals->SubController->NoBlockSend(
        &sendIds[0], sizesToSend[reqidx], *iter, UGGCG_DATA_EXCHANGE_TAG, sendReqs[2 * reqidx + 1]);
550
      reqidx++;
551 552 553 554 555
    }

    // loop until all sizes are received
    size_t counter = 0;
    size_t numNeighbors = this->Internals->Neighbors.size();
556
    while (counter != numNeighbors)
557 558
    {
      iter = this->Internals->Neighbors.begin();
559
      for (; iter != this->Internals->Neighbors.end(); ++iter)
560 561 562 563 564 565 566 567 568 569 570 571 572
      {
        CommDataInfo& c = this->Internals->CommData[*iter];
        if (!c.RecvReqs[0].Test() || c.CommStep != 0)
        {
          continue;
        }
        c.CommStep = 1;
        counter++;
      }
    }

    // create receive requests for the ids
    iter = this->Internals->Neighbors.begin();
573
    for (; iter != this->Internals->Neighbors.end(); ++iter)
574 575 576
    {
      CommDataInfo& c = this->Internals->CommData[*iter];
      this->Internals->ProcessIdToSurfacePointIds[*iter].resize(c.RecvSize);
577
      this->Internals->SubController->NoBlockReceive(
578
        &this->Internals->ProcessIdToSurfacePointIds[*iter][0], c.RecvSize, *iter,
579 580 581 582 583
        UGGCG_DATA_EXCHANGE_TAG, c.RecvReqs[1]);
    }

    // wait for receives
    counter = 0;
584
    while (counter != numNeighbors)
585 586
    {
      iter = this->Internals->Neighbors.begin();
587
      for (; iter != this->Internals->Neighbors.end(); ++iter)
588 589 590 591 592 593 594 595 596 597 598
      {
        CommDataInfo& c = this->Internals->CommData[*iter];
        if (!c.RecvReqs[1].Test() || c.CommStep != 1)
        {
          continue;
        }
        c.CommStep = 2;
        counter++;
      }
    }
    // should have all id data by now
599
  }
600
  else
601
  {
602
    // We can't use global ids, so we will process point coordinates instead
603 604
    // send surface points to all neighbors
    // could potentially just send points that are in a neighbor's bounding box
605
    this->Internals->ProcessIdToSurfacePoints.clear();
606
    this->Internals->SendPoints.clear();
607 608 609 610 611 612
    vtkPoints* surfacePoints = surface->GetPoints();
    this->Internals->LocalPointsMap = surfaceOriginalPointIds;

    std::vector<int>::iterator iter = this->Internals->Neighbors.begin();
    // get all sizes from neighbors
    // first set up the receives
613
    for (; iter != this->Internals->Neighbors.end(); ++iter)
614
    {
615
      CommDataInfo& c = this->Internals->CommData[*iter];
616 617
      this->Internals->SubController->NoBlockReceive(
        &c.RecvSize, 1, *iter, UGGCG_SIZE_EXCHANGE_TAG, c.RecvReqs[0]);
618
    }
619

620 621 622 623 624
    // keep my own points
    this->Internals->MyPoints = surfacePoints->GetData();

    // store the global point arrays unique to each process (based on bounding box
    // of the receiving process) to send
625 626
    for (iter = this->Internals->Neighbors.begin(); iter != this->Internals->Neighbors.end();
         ++iter)
627 628 629 630 631 632 633
    {
      std::vector<double>& sendPoints = this->Internals->SendPoints[*iter];
      sendPoints.clear();
      for (vtkIdType i = 0; i < nbSurfacePoints; i++)
      {
        double coord[3];
        surface->GetPoint(i, coord);
634
        if (vtkMath::PointIsWithinBounds(coord, &allBounds[*iter * 6], delta))
635
        {
636
          sendPoints.insert(sendPoints.end(), coord, coord + 3);
637 638 639 640 641 642
        }
      }
    }

    // now go through and send the data
    int reqidx = 0;
643 644
    for (iter = this->Internals->Neighbors.begin(); iter != this->Internals->Neighbors.end();
         ++iter)
645
    {
646
      // Send data length
647
      std::vector<double>& sendPoints = this->Internals->SendPoints[*iter];
648
      sizesToSend[reqidx] = static_cast<int>(sendPoints.size());
649 650
      this->Internals->SubController->NoBlockSend(
        &sizesToSend[reqidx], 1, *iter, UGGCG_SIZE_EXCHANGE_TAG, sendReqs[2 * reqidx]);
651 652

      // Send raw data
653
      this->Internals->SubController->NoBlockSend(&sendPoints[0], sizesToSend[reqidx], *iter,
654
        UGGCG_DATA_EXCHANGE_TAG, sendReqs[2 * reqidx + 1]);
655 656 657 658 659 660
      reqidx++;
    }

    // loop until all sizes are received
    size_t counter = 0;
    size_t numNeighbors = this->Internals->Neighbors.size();
661
    while (counter != numNeighbors)
662
    {
663 664
      for (iter = this->Internals->Neighbors.begin(); iter != this->Internals->Neighbors.end();
           ++iter)
665 666 667
      {
        CommDataInfo& c = this->Internals->CommData[*iter];
        if (!c.RecvReqs[0].Test() || c.CommStep != 0)
668
        {
669
          continue;
670
        }
671 672
        c.CommStep = 1;
        counter++;
673
      }
674
    }
675

676
    // create receive requests for point data
677 678
    for (iter = this->Internals->Neighbors.begin(); iter != this->Internals->Neighbors.end();
         ++iter)
679 680
    {
      CommDataInfo& c = this->Internals->CommData[*iter];
681 682
      std::vector<double>& incomingPoints = this->Internals->ProcessIdToSurfacePoints[*iter];
      incomingPoints.resize(c.RecvSize);
683 684
      this->Internals->SubController->NoBlockReceive(
        &incomingPoints[0], c.RecvSize, *iter, UGGCG_DATA_EXCHANGE_TAG, c.RecvReqs[1]);
685 686
    }

687
    // wait for receives of data
688
    counter = 0;
689
    while (counter != numNeighbors)
690 691
    {
      iter = this->Internals->Neighbors.begin();
692
      for (; iter != this->Internals->Neighbors.end(); ++iter)
693 694 695 696 697 698 699 700 701 702
      {
        CommDataInfo& c = this->Internals->CommData[*iter];
        if (!c.RecvReqs[1].Test() || c.CommStep != 1)
        {
          continue;
        }
        c.CommStep = 2;
        counter++;
      }
    }
703
  }
704 705
  // should have all point data by now
  // wait for all my sends to complete
706
  this->Internals->SubController->WaitAll(static_cast<int>(sendReqs.size()), &sendReqs[0]);
707 708 709
}

//---------------------------------------------------------------------------
710
// Step 2a: browse global ids/point coordinates of other ranks and check if some
711 712 713 714 715
// are duplicated locally.
// For each neighbor rank, save the ids of the cells adjacent to the surface
// points shared, those cells are the ghost cells we will send them.
void vtkPUnstructuredGridGhostCellsGenerator::ComputeSharedPoints()
{
716
  this->Internals->CellsToSend.clear();
717
  vtkNew<vtkIdList> cellIdsList;
718
  if (this->Internals->InputGlobalPointIds)
719
  {
720
    for (std::vector<int>::iterator iter = this->Internals->Neighbors.begin();
721
         iter != this->Internals->Neighbors.end(); ++iter)
722
    {
723 724 725
      std::vector<vtkIdType>& surfaceIds = this->Internals->ProcessIdToSurfacePointIds[*iter];
      for (std::vector<vtkIdType>::const_iterator id_iter = surfaceIds.begin();
           id_iter != surfaceIds.end(); ++id_iter)
726
      {
727
        vtkIdType localPointId = -1;
728 729
        // Check if this point exists locally from its global ids, if so
        // get its local id.
730 731
        vtkIdType gid = *id_iter;
        std::map<vtkIdType, vtkIdType>::iterator miter =
732
          this->Internals->GlobalToLocalPointIdMap.find(gid);
733
        if (miter != this->Internals->GlobalToLocalPointIdMap.end())
734
        {
735 736 737 738 739 740 741 742 743 744 745 746
          localPointId = miter->second;
          if (localPointId != -1)
          {
            // Current rank also has a copy of this global point
            // Get the cells connected to this point
            this->Internals->Input->GetPointCells(localPointId, cellIdsList.Get());
            vtkIdType nbIds = cellIdsList->GetNumberOfIds();
            // Add those cells to the list of cells to send to this rank
            for (vtkIdType k = 0; k < nbIds; k++)
            {
              this->Internals->CellsToSend[*iter].insert(cellIdsList->GetId(k));
              this->Internals->SentCellsLastRound[*iter].insert(cellIdsList->GetId(k));
747
              this->Internals->SentCells[*iter].insert(cellIdsList->GetId(k));
748 749
            }
          }
750
        }
751
      }
752 753 754 755 756 757 758 759
    }
  }
  else
  {
    // build kdtree of local surface points
    vtkNew<vtkKdTree> kdtree;
    vtkNew<vtkPoints> points;
    int myRank = this->Internals->SubController->GetLocalProcessId();
760
    points->SetData(this->Internals->MyPoints);
761 762 763
    kdtree->BuildLocatorFromPoints(points);
    double bounds[6];
    kdtree->GetBounds(bounds);
764 765 766
    double tolerance = 1.e-6 * sqrt((bounds[1] - bounds[0]) * (bounds[1] - bounds[0]) +
                                 (bounds[3] - bounds[2]) * (bounds[3] - bounds[2]) +
                                 (bounds[5] - bounds[4]) * (bounds[5] - bounds[4]));
767

768 769 770
    for (std::map<int, std::vector<double> >::iterator iter =
           this->Internals->ProcessIdToSurfacePoints.begin();
         iter != this->Internals->ProcessIdToSurfacePoints.end(); ++iter)
771 772
    {
      if (iter->first == myRank)
773
      {
774
        continue;
775
      }
776
      std::vector<double>& offProcSurfacePoints = iter->second;
777
      double dist2(0); // result will be distance squared
778
      for (size_t i = 0; i < offProcSurfacePoints.size(); i += 3)
779
      {
780 781
        vtkIdType id =
          kdtree->FindClosestPointWithinRadius(tolerance, &offProcSurfacePoints[i], dist2);
782 783 784
        if (id != -1)
        { // matching point...
          vtkIdType inputId = this->Internals->LocalPointsMap->GetValue(id);
785
          this->Internals->Input->GetPointCells(inputId, cellIdsList);
786 787 788 789 790
          // Add those cells to the list of cells to send to this rank
          for (vtkIdType k = 0; k < cellIdsList->GetNumberOfIds(); k++)
          {
            this->Internals->CellsToSend[iter->first].insert(cellIdsList->GetId(k));
            this->Internals->SentCellsLastRound[iter->first].insert(cellIdsList->GetId(k));
791
            this->Internals->SentCells[iter->first].insert(cellIdsList->GetId(k));
792
          }
793 794 795
        }
      }
    }
796
  }
797 798

  // Release memory of all reduced arrays
799 800 801
  this->Internals->ProcessIdToSurfacePointIds.clear();
  this->Internals->ProcessIdToSurfacePoints.clear();
  this->Internals->LocalPointsMap = nullptr;
802 803
  this->Internals->SendIds.clear();
  this->Internals->MyPoints = nullptr;
804 805 806 807 808 809
  // Now we know our neighbors and which points we have in common and the
  // ghost cells to share.
}

//-----------------------------------------------------------------------------
// Step 3: extract and send the ghost cells to the neighbor ranks
810 811
void vtkPUnstructuredGridGhostCellsGenerator::ExtractAndSendGhostCells(
  vtkUnstructuredGridBase* input)
812 813 814
{
  vtkNew<vtkIdList> cellIdsList;
  vtkNew<vtkExtractCells> extractCells;
815
  extractCells->SetInputData(input);
816

817 818
  for (std::vector<int>::iterator iter = this->Internals->Neighbors.begin();
       iter != this->Internals->Neighbors.end(); ++iter)
819
  {
820 821 822 823 824 825
    int toRank = *iter;
    CommDataInfo& c = this->Internals->CommData[toRank];
    std::map<int, std::set<vtkIdType> >::iterator miter = this->Internals->CellsToSend.find(toRank);
    if (miter == this->Internals->CellsToSend.end())
    { // no data to send
      c.SendLen = 0;
826 827
      this->Internals->SubController->NoBlockSend(
        &c.SendLen, 1, toRank, UGGCG_SIZE_EXCHANGE_TAG, c.SendReqs[0]);
828 829 830
      continue;
    }
    std::set<vtkIdType>& cellsToShare = miter->second;
831
    cellIdsList->SetNumberOfIds(static_cast<vtkIdType>(cellsToShare.size()));
832 833
    std::set<vtkIdType>::iterator sIter = cellsToShare.begin();
    for (vtkIdType i = 0; sIter != cellsToShare.end(); ++sIter, i++)
834
    {
835
      cellIdsList->SetId(i, *sIter);
836
    }
837
    extractCells->SetCellList(cellIdsList);
838
    extractCells->Update();
839
    vtkUnstructuredGrid* extractGrid = extractCells->GetOutput();
840 841 842 843 844 845

    //There might be case where the originalcellids needs to be removed
    //but there are definitely cases where it shouldn't.
    //So if you run into that case, think twice before you uncomment this
    //next line and look carefully at paraview issue #18470
    //extractGrid->GetCellData()->RemoveArray("vtkOriginalCellIds");
846 847 848

    // Send the extracted grid to the neighbor rank asynchronously
    if (vtkCommunicator::MarshalDataObject(extractGrid, c.SendBuffer))
849
    {
850 851
      c.SendLen = c.SendBuffer->GetNumberOfTuples();
      // Send data length
852 853
      this->Internals->SubController->NoBlockSend(
        &c.SendLen, 1, toRank, UGGCG_SIZE_EXCHANGE_TAG, c.SendReqs[0]);
854

855
      // Send raw data
856
      this->Internals->SubController->NoBlockSend((char*)c.SendBuffer->GetVoidPointer(0), c.SendLen,
857
        toRank, UGGCG_DATA_EXCHANGE_TAG, c.SendReqs[1]);
858
    }
859
  }
860 861 862
}

//-----------------------------------------------------------------------------
863
// Step 4: Receive the ghost cells from the neighbor ranks and merge them
864
// to the local grid.
865
// Argument output should be an empty unstructured grid.
866 867
void vtkPUnstructuredGridGhostCellsGenerator::ReceiveAndMergeGhostCells(int ghostLevel,
  int maxGhostLevel, vtkUnstructuredGridBase* currentGrid, vtkUnstructuredGrid* output)
868
{
869
  // reset CommStep
870 871
  assert(this->Internals->Neighbors.size() == this->Internals->CommData.size());
  for (std::map<int, CommDataInfo>::iterator comIter = this->Internals->CommData.begin();
872
       comIter != this->Internals->CommData.end(); ++comIter)
873
  {
874
    comIter->second.CommStep = 0;
875
  }
876

877 878
  // We need to compute a rough estimation of the total number of cells and
  // points for vtkMergeCells
879 880
  vtkIdType totalNbCells = currentGrid->GetNumberOfCells();
  vtkIdType totalNbPoints = currentGrid->GetNumberOfPoints();
881

882
  // Browse all neighbor ranks and receive the mesh that contains cells
883
  size_t nbNeighbors = this->Internals->Neighbors.size();
884 885 886 887
  std::vector<vtkUnstructuredGridBase*> neighborGrids;
  neighborGrids.reserve(nbNeighbors);

  // First create requests to receive the size of the mesh to receive
888 889
  std::vector<int>::iterator iter = this->Internals->Neighbors.begin();
  for (; iter != this->Internals->Neighbors.end(); ++iter)
890
  {
891
    int fromRank = *iter;
892
    CommDataInfo& c = this->Internals->CommData[fromRank];
893
    this->Internals->SubController->NoBlockReceive(
894
      &c.RecvLen, 1, fromRank, UGGCG_SIZE_EXCHANGE_TAG, c.RecvReqs[0]);
895
  }
896

897 898
  // Then, once the data length is received, create requests to receive the
  // mesh data
899 900
  size_t counter = 0;
  size_t nonEmptyNeighborCounter = 0; // some neighbors might not have data to send
901
  while (counter != nbNeighbors)
902
  {
903 904
    iter = this->Internals->Neighbors.begin();
    for (; iter != this->Internals->Neighbors.end(); ++iter)
905
    {
906
      int fromRank = *iter;
907 908
      CommDataInfo& c = this->Internals->CommData[fromRank];
      if (!c.RecvReqs[0].Test() || c.CommStep != 0)
909
      {
910
        continue;
911
      }
912 913 914 915
      if (c.RecvLen > 0)
      {
        c.CommStep = 1; // mark that this comm needs to receive the dataset
        c.RecvBuffer->SetNumberOfValues(c.RecvLen);
916 917
        this->Internals->SubController->NoBlockReceive((char*)c.RecvBuffer->GetVoidPointer(0),
          c.RecvLen, fromRank, UGGCG_DATA_EXCHANGE_TAG, c.RecvReqs[1]);
918 919 920 921 922 923
        nonEmptyNeighborCounter++;
      }
      else
      {
        c.CommStep = 2; // mark that this comm doesn't need to receive the dataset
      }
924 925
      counter++;
    }
926
  }
927 928 929 930

  // Browse all neighbor ranks and receive the mesh that contains cells
  // that are ghost cells for current rank.
  counter = 0;
931
  while (counter != nonEmptyNeighborCounter)
932
  {
933 934
    iter = this->Internals->Neighbors.begin();
    for (; iter != this->Internals->Neighbors.end(); ++iter)
935
    {
936
      int fromRank = *iter;
937 938 939
      CommDataInfo& c = this->Internals->CommData[fromRank];

      if (!c.RecvReqs[1].Test() || c.CommStep != 1)
940
      {
941
        continue;
942
      }
943 944

      c.CommStep = 2;
945
      vtkUnstructuredGrid* grid = vtkUnstructuredGrid::New();
946
      vtkCommunicator::UnMarshalDataObject(c.RecvBuffer, grid);
947 948
      // clear out some memory...
      c.RecvBuffer->SetNumberOfTuples(0);
949 950

      if (!grid->HasAnyGhostCells())
951
      {
952 953
        grid->AllocatePointGhostArray();
        grid->AllocateCellGhostArray();
954
      }
955 956 957 958 959

      // Flag the received grid elements as ghosts
      grid->GetPointGhostArray()->FillComponent(0, 1);
      grid->GetCellGhostArray()->FillComponent(0, 1);

960 961
      // record all cells that i received
      // only needed if we need to calculate more ghost layers
962
      if (ghostLevel < maxGhostLevel)
963
      {
964
        if (grid->GetCellData()->GetGlobalIds())
965
        {
966 967
          vtkIdTypeArray* cellids =
            vtkArrayDownCast<vtkIdTypeArray>(grid->GetCellData()->GetGlobalIds());
968
          for (vtkIdType i = 0; i < grid->GetNumberOfCells(); i++)
969
          {
970
            this->Internals->ReceivedCells[fromRank].insert(cellids->GetValue(i));
971 972 973 974
          }
        }
      }

975
      // Make sure the global point ids array is tagged accordingly
976
      if (this->Internals->InputGlobalPointIds && !grid->GetPointData()->GetGlobalIds())
977
      {
978 979
        grid->GetPointData()->SetGlobalIds(
          grid->GetPointData()->GetArray(this->Internals->InputGlobalPointIds->GetName()));
980
      }
981

982 983 984
      // Checking maxGhostLevel to see if global cell ids are needed.
      // If so, make sure the global cell ids array is tagged accordingly
      if (maxGhostLevel > 1)
985
      {
986 987
        if (!grid->GetCellData()->GetGlobalIds())
        {
988 989
          grid->GetCellData()->SetGlobalIds(
            grid->GetCellData()->GetArray(this->GlobalCellIdsArrayName));
990
        }
991
      }
992

993 994 995 996 997 998 999
      totalNbCells += grid->GetNumberOfCells();
      totalNbPoints += grid->GetNumberOfPoints();

      neighborGrids.push_back(grid);

      counter++;
    }
1000
  }
1001

1002
  if (totalNbCells == 0)
1003
  {
1004
    output->ShallowCopy(currentGrid);
1005
    return;
1006
  }
1007

1008
  // Use MergeCells to merge currentGrid + new grids to the output grid
1009 1010 1011 1012
  vtkNew<vtkMergeCells> mergeCells;
  mergeCells->SetUnstructuredGrid(output);
  mergeCells->SetTotalNumberOfCells(totalNbCells);
  mergeCells->SetTotalNumberOfPoints(totalNbPoints);
1013
  mergeCells->SetTotalNumberOfDataSets(1 + static_cast<int>(neighborGrids.size()));
1014
  mergeCells->SetUseGlobalIds(this->Internals->InputGlobalPointIds != nullptr ? 1 : 0);
1015
  mergeCells->SetPointMergeTolerance(0.0);
1016
  mergeCells->SetUseGlobalCellIds(1);
1017

1018
  // Merge current grid first
1019
  mergeCells->MergeDataSet(currentGrid);
1020

1021
  // Then merge ghost grid from neighbor ranks
1022
  for (std::size_t i = 0; i < neighborGrids.size(); i++)
1023
  {
1024 1025 1026 1027
    if (neighborGrids[i]->GetNumberOfCells())
    {
      mergeCells->MergeDataSet(neighborGrids[i]);
    }
1028
    neighborGrids[i]->Delete();
1029
  }
1030 1031 1032

  // Finalize the merged output
  mergeCells->Finish();
1033

1034 1035 1036 1037
  // for all ghost cells, store the global cell id to local cell id mapping.
  // we need this mapping later when determining if cells we want to send
  // have been received before. only needed if we are calculating more
  // ghost layers.
1038
  if (ghostLevel < maxGhostLevel)
1039
  {
1040 1041
    vtkDataArray* ghost = output->GetCellGhostArray();
    vtkDataArray* gids = output->GetCellData()->GetGlobalIds();
1042
    for (vtkIdType lid = 0; lid < output->GetNumberOfCells(); lid++)
1043
    {
1044
      if (ghost->GetTuple1(lid) > 0)
1045
      {
1046 1047
        vtkIdType gid = static_cast<vtkIdType>(gids->GetTuple1(lid));
        if (this->Internals->GlobalToLocalCellIdMap.find(gid) ==
1048
          this->Internals->GlobalToLocalCellIdMap.end())
1049
        {
1050
          this->Internals->GlobalToLocalCellIdMap[gid] = lid;
1051 1052 1053 1054 1055 1056 1057 1058 1059
        }
      }
    }
  }

  // wait here on the sends to make sure we don't corrupt the data before it's fully sent
  counter = 0;
  while (counter != nbNeighbors)
  {
1060 1061
    for (iter = this->Internals->Neighbors.begin(); iter != this->Internals->Neighbors.end();
         ++iter)
1062 1063 1064 1065 1066
    {
      int toRank = *iter;
      CommDataInfo& c = this->Internals->CommData[toRank];
      std::map<int, std::set<vtkIdType> >::iterator miter =
        this->Internals->CellsToSend.find(toRank);
1067
      if (miter == this->Internals->CellsToSend.end())
1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093
      {
        // this is a process that we don't send cells to so we only need to check
        // that the message with the buffer size doesn't get modified
        if (c.CommStep == 3 || !c.SendReqs[0].Test())
        {
          continue;
        }
        c.CommStep = 3;
        counter++;
      }
      else
      {
        if (c.CommStep == 3 || !c.SendReqs[1].Test())
        {
          continue;
        }
        c.CommStep = 3;
        counter++;
        // clear out some memory...
        if (c.SendBuffer)
        {
          c.SendBuffer->SetNumberOfTuples(0);
        }
      }
    }
  }
1094 1095 1096 1097 1098
}

//-----------------------------------------------------------------------------
// Add another ghost layer. Assumes that at least one layer of ghost cells has
// already been created. Must be called after GetFirstGhostLayer.
1099
void vtkPUnstructuredGridGhostCellsGenerator::AddGhostLayer(int ghostLevel, int maxGhostLevel)
1100 1101 1102
{
  this->Internals->CellsToSend.clear();
  this->FindGhostCells();
1103
  this->UpdateProgress((1.0 + ((ghostLevel - 1) * 3.0)) / (maxGhostLevel * 3.0));
1104

1105
  this->ExtractAndSendGhostCells(this->Internals->CurrentGrid);
1106 1107 1108 1109 1110
  this->UpdateProgress((2.0 + ((ghostLevel - 1) * 3.0)) / (maxGhostLevel * 3.0));
  vtkSmartPointer<vtkUnstructuredGrid> outputGrid = vtkSmartPointer<vtkUnstructuredGrid>::New();
  this->