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