1 /*=========================================================================
3 Program: Visualization Toolkit
4 Module: vtkPUnstructuredGridGhostCellsGenerator.cxx
6 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
8 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
10 This software is distributed WITHOUT ANY WARRANTY; without even
11 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12 PURPOSE. See the above copyright notice for more information.
14 =========================================================================*/
16 #include "vtkPUnstructuredGridGhostCellsGenerator.h"
18 #include "vtkCellArray.h"
19 #include "vtkCharArray.h"
20 #include "vtkDataSetSurfaceFilter.h"
21 #include "vtkExtractCells.h"
22 #include "vtkIdList.h"
23 #include "vtkIdTypeArray.h"
24 #include "vtkInformation.h"
25 #include "vtkInformationVector.h"
26 #include "vtkMergeCells.h"
27 #include "vtkMergePoints.h"
28 #include "vtkMPICommunicator.h"
29 #include "vtkMultiProcessController.h"
31 #include "vtkObjectFactory.h"
32 #include "vtkPointData.h"
33 #include "vtkPoints.h"
34 #include "vtkStreamingDemandDrivenPipeline.h"
35 #include "vtkUnstructuredGrid.h"
41 //----------------------------------------------------------------------------
46 bool AllGatherV(vtkMultiProcessController* controller,
49 std::vector<T>& globalV,
50 std::vector<vtkIdType>& sizes,
51 std::vector<vtkIdType>& offsets)
53 int nbOfRanks = controller->GetNumberOfProcesses();
54 sizes.resize(nbOfRanks);
55 int ret = controller->AllGather(&localSize, &sizes[0], 1);
58 vtkErrorWithObjectMacro(controller, << "Communication error!");
62 offsets.resize(nbOfRanks);
63 for (int i = 0; i < nbOfRanks; i++)
68 globalV.resize(count);
71 controller->AllGatherV(localSize > 0 ? localV : 0,
72 &globalV[0], localSize, &sizes[0], &offsets[0]);
78 //----------------------------------------------------------------------------
79 // Internal data structures
81 // Class to hold asynchronous communication information
85 CommDataInfo() : SendLen(-1), RecvLen(-1), CommStep(0)
87 this->SendBuffer = vtkCharArray::New();
88 this->RecvBuffer = vtkCharArray::New();
91 CommDataInfo(const CommDataInfo& c)
94 if (this->SendBuffer) { this->SendBuffer->Register(0); }
95 if (this->RecvBuffer) { this->RecvBuffer->Register(0); }
100 if (this->SendBuffer) { this->SendBuffer->Delete(); }
101 if (this->RecvBuffer) { this->RecvBuffer->Delete(); }
104 vtkMPICommunicator::Request SendReqs[2];
105 vtkMPICommunicator::Request RecvReqs[2];
106 vtkCharArray *SendBuffer;
107 vtkCharArray *RecvBuffer;
113 // Communication arrays
114 struct vtkPUnstructuredGridGhostCellsGenerator::vtkInternals
117 std::map<vtkIdType, vtkIdType> GlobalToLocalPointIdMap;
118 std::vector<vtkIdType> AllGlobalIdsOfSurfacePoints;
120 // For point coordinates
121 vtkNew<vtkMergePoints> LocalPoints;
122 std::vector<vtkIdType> LocalPointsMap;
123 std::vector<double> AllPointsOfSurfacePoints;
125 std::vector<vtkIdType> AllSizes;
126 std::vector<vtkIdType> AllOffsets;
127 std::map<int, std::set<vtkIdType> > NeighborRanksCells;
128 std::map<int, CommDataInfo> CommData;
129 vtkUnstructuredGridBase* Input;
131 vtkDataArray* InputGlobalPointIds;
132 bool UseGlobalPointIds;
135 static const int UGGCG_SIZE_EXCHANGE_TAG = 9000;
136 static const int UGGCG_DATA_EXCHANGE_TAG = 9001;
137 static const char* UGGCG_GLOBAL_POINT_IDS = "GlobalNodeIds";
139 //----------------------------------------------------------------------------
141 vtkStandardNewMacro(vtkPUnstructuredGridGhostCellsGenerator)
142 vtkSetObjectImplementationMacro(
143 vtkPUnstructuredGridGhostCellsGenerator, Controller, vtkMultiProcessController);
145 //----------------------------------------------------------------------------
146 vtkPUnstructuredGridGhostCellsGenerator::vtkPUnstructuredGridGhostCellsGenerator()
148 this->Controller = NULL;
149 this->SetController(vtkMultiProcessController::GetGlobalController());
151 this->Internals = NULL;
152 this->GlobalPointIdsArrayName = NULL;
153 this->UseGlobalPointIds = true;
154 this->SetGlobalPointIdsArrayName(UGGCG_GLOBAL_POINT_IDS);
157 //----------------------------------------------------------------------------
158 vtkPUnstructuredGridGhostCellsGenerator::~vtkPUnstructuredGridGhostCellsGenerator()
160 this->SetController(NULL);
161 this->SetGlobalPointIdsArrayName(NULL);
163 delete this->Internals;
167 //-----------------------------------------------------------------------------
168 void vtkPUnstructuredGridGhostCellsGenerator::PrintSelf(ostream& os, vtkIndent indent)
170 Superclass::PrintSelf(os, indent);
173 //-----------------------------------------------------------------------------
174 int vtkPUnstructuredGridGhostCellsGenerator::RequestData(
175 vtkInformation *vtkNotUsed(request),
176 vtkInformationVector **inputVector,
177 vtkInformationVector *outputVector)
179 // get the info objects
180 vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
181 vtkInformation *outInfo = outputVector->GetInformationObject(0);
183 // get the input and output. Input may just have the UnstructuredGridBase
184 // interface, but output should be an unstructured grid.
185 vtkUnstructuredGridBase *input = vtkUnstructuredGridBase::SafeDownCast(
186 inInfo->Get(vtkDataObject::DATA_OBJECT()));
187 vtkUnstructuredGrid *output = vtkUnstructuredGrid::SafeDownCast(
188 outInfo->Get(vtkDataObject::DATA_OBJECT()));
191 //ghostLevels = outInfo->Get(
192 // vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS());
194 if (!this->Controller)
196 this->Controller = vtkMultiProcessController::GetGlobalController();
200 this->NumRanks = this->Controller ? this->Controller->GetNumberOfProcesses() : 1;
201 this->RankId = this->Controller ? this->Controller->GetLocalProcessId() : 0;
202 if (ghostLevels == 0 || !this->Controller || this->NumRanks == 1)
204 // Ghost levels are not requested. Nothing to do but pass the dataset.
205 if (this->RankId == 0)
207 vtkWarningMacro(<< "Ghost cells not requested or not needed.");
209 output->ShallowCopy(input);
213 delete this->Internals;
214 this->Internals = new vtkPUnstructuredGridGhostCellsGenerator::vtkInternals();
215 this->Internals->Input = input;
217 vtkPointData *inputPD = input->GetPointData();
218 this->Internals->InputGlobalPointIds = inputPD->GetGlobalIds();
220 if (!this->Internals->InputGlobalPointIds)
222 this->Internals->InputGlobalPointIds =
223 inputPD->GetArray(this->GlobalPointIdsArrayName);
224 inputPD->SetGlobalIds(this->Internals->InputGlobalPointIds);
227 if (!this->UseGlobalPointIds)
229 this->Internals->InputGlobalPointIds = NULL;
232 this->ExtractAndReduceSurfacePoints();
233 this->UpdateProgress(0.3);
235 this->ComputeSharedPoints();
236 this->UpdateProgress(0.6);
238 this->ExtractAndSendGhostCells();
239 this->UpdateProgress(0.8);
241 this->ReceiveAndMergeGhostCells(output);
242 this->UpdateProgress(1.0);
244 output->GetInformation()->Set(vtkDataObject::DATA_NUMBER_OF_GHOST_LEVELS(), 1);
246 this->Controller->Barrier();
248 delete this->Internals;
254 //-----------------------------------------------------------------------------
255 // Step 1: Extract surface geometry and all reduce global ids of surface points
256 void vtkPUnstructuredGridGhostCellsGenerator::ExtractAndReduceSurfacePoints()
258 // Extract boundary cells and points with the surface filter
259 vtkNew<vtkDataSetSurfaceFilter> surfaceFilter;
260 surfaceFilter->SetInputData(this->Internals->Input);
261 surfaceFilter->PassThroughPointIdsOn();
262 surfaceFilter->Update();
264 vtkPolyData* surface = surfaceFilter->GetOutput();
265 vtkIdType nbSurfacePoints = surface->GetNumberOfPoints();
266 vtkCellArray* surfaceCells = surface->GetPolys();
267 surfaceCells->InitTraversal();
268 vtkIdType npts, *pts;
270 vtkIdTypeArray* surfaceOriginalPointIds = vtkIdTypeArray::SafeDownCast(
271 surface->GetPointData()->GetArray(surfaceFilter->GetOriginalPointIdsName()));
273 if (this->Internals->InputGlobalPointIds)
275 std::vector<vtkIdType> globalIdsOfSurfacePoints;
276 globalIdsOfSurfacePoints.reserve(nbSurfacePoints);
278 // Browse surface cells and save global and local ids of cell points
279 while (surfaceCells->GetNextCell(npts, pts))
281 for (vtkIdType i = 0; i < npts; i++)
283 vtkIdType origPtId = surfaceOriginalPointIds->GetValue(pts[i]);
284 vtkIdType globalPtId = static_cast<vtkIdType>(
285 this->Internals->InputGlobalPointIds->GetTuple1(origPtId));
287 if (this->Internals->GlobalToLocalPointIdMap.find(globalPtId) ==
288 this->Internals->GlobalToLocalPointIdMap.end())
290 this->Internals->GlobalToLocalPointIdMap[globalPtId] = origPtId;
291 globalIdsOfSurfacePoints.push_back(globalPtId);
296 // Now reduce surface point global ids on ALL ranks
297 ::AllGatherV(this->Controller, &globalIdsOfSurfacePoints[0],
298 globalIdsOfSurfacePoints.size(),
299 this->Internals->AllGlobalIdsOfSurfacePoints,
300 this->Internals->AllSizes, this->Internals->AllOffsets);
304 // We can't use global ids, so we will process point coordinates instead
305 vtkPoints* inputPoints = this->Internals->Input->GetPoints();
306 vtkNew<vtkPoints> surfacePoints;
307 surfacePoints->SetDataTypeToDouble();
308 surfacePoints->Allocate(nbSurfacePoints);
309 this->Internals->LocalPoints->InitPointInsertion(
310 surfacePoints.Get(), surface->GetPoints()->GetBounds());
311 this->Internals->LocalPointsMap.reserve(nbSurfacePoints);
313 // Browse surface cells and push point coordinates to the locator
314 while (surfaceCells->GetNextCell(npts, pts))
316 for (vtkIdType i = 0; i < npts; i++)
318 vtkIdType origPtId = surfaceOriginalPointIds->GetValue(pts[i]);
320 inputPoints->GetPoint(origPtId, p);
322 if (this->Internals->LocalPoints->InsertUniquePoint(p, sid))
324 // New point, save the id of the original grid point id associated
325 // to this surface point
326 if (static_cast<vtkIdType>(this->Internals->LocalPointsMap.size()) <= sid)
328 this->Internals->LocalPointsMap.resize(sid + 1);
330 this->Internals->LocalPointsMap[sid] = origPtId;
335 // Now reduce surface point coordinates on ALL ranks
336 ::AllGatherV(this->Controller,
337 (double*)surfacePoints->GetVoidPointer(0),
338 surfacePoints->GetNumberOfPoints() * 3,
339 this->Internals->AllPointsOfSurfacePoints,
340 this->Internals->AllSizes, this->Internals->AllOffsets);
344 //---------------------------------------------------------------------------
345 // Step 2: browse global ids/point coordinates of other ranks and check if some
346 // are duplicated locally.
347 // For each neighbor rank, save the ids of the cells adjacent to the surface
348 // points shared, those cells are the ghost cells we will send them.
349 void vtkPUnstructuredGridGhostCellsGenerator::ComputeSharedPoints()
351 vtkNew<vtkIdList> cellIdsList;
352 for (int i = 0; i < this->NumRanks; i++)
354 if (i == this->RankId)
358 for (vtkIdType j = 0, idx = this->Internals->AllOffsets[i];
359 j < this->Internals->AllSizes[i]; j++, idx++)
361 vtkIdType localPointId = -1;
362 if (this->Internals->InputGlobalPointIds)
364 // Check if this point exists locally from its global ids, if so
366 vtkIdType gid = this->Internals->AllGlobalIdsOfSurfacePoints[idx];
367 std::map<vtkIdType, vtkIdType>::iterator iter =
368 this->Internals->GlobalToLocalPointIdMap.find(gid);
369 if (iter != this->Internals->GlobalToLocalPointIdMap.end())
371 localPointId = iter->second;
376 // Check if this point exists locally from its coordinates, if so
378 double *p = &this->Internals->AllPointsOfSurfacePoints[idx];
379 localPointId = this->Internals->LocalPoints->IsInsertedPoint(p);
381 if (localPointId != -1)
383 localPointId = this->Internals->LocalPointsMap[localPointId];
385 idx += 2; // jump to next coordinates
389 if (localPointId != -1)
391 // Current rank also has a copy of this global point
392 cellIdsList->Reset();
393 // Get the cells connected to this point
394 this->Internals->Input->GetPointCells(localPointId, cellIdsList.Get());
395 vtkIdType nbIds = cellIdsList->GetNumberOfIds();
396 // Add those cells to the list of cells to send to this rank
397 for (vtkIdType k = 0; k < nbIds; k++)
399 this->Internals->NeighborRanksCells[i].insert(cellIdsList->GetId(k));
405 // Release memory of all reduced arrays
406 this->Internals->AllGlobalIdsOfSurfacePoints.resize(0);
407 this->Internals->AllPointsOfSurfacePoints.resize(0);
408 this->Internals->AllSizes.resize(0);
409 this->Internals->AllOffsets.resize(0);
410 // Now we know our neighbors and which points we have in common and the
411 // ghost cells to share.
414 //-----------------------------------------------------------------------------
415 // Step 3: extract and send the ghost cells to the neighbor ranks
416 void vtkPUnstructuredGridGhostCellsGenerator::ExtractAndSendGhostCells()
418 vtkNew<vtkIdList> cellIdsList;
419 vtkNew<vtkExtractCells> extractCells;
420 extractCells->SetInputData(this->Internals->Input);
421 std::map<int, std::set<vtkIdType> >::iterator iter =
422 this->Internals->NeighborRanksCells.begin();
424 vtkMPICommunicator* com =
425 vtkMPICommunicator::SafeDownCast(this->Controller->GetCommunicator());
427 // Browse all neighbor ranks and extract the cells connected to the points we share
428 for (; iter != this->Internals->NeighborRanksCells.end(); ++iter)
430 int toRank = iter->first;
431 std::set<vtkIdType>& cellsToShare = iter->second;
432 cellIdsList->SetNumberOfIds(cellsToShare.size());
433 std::set<vtkIdType>::iterator sIter = cellsToShare.begin();
434 for (vtkIdType i = 0; sIter != cellsToShare.end(); ++sIter, i++)
436 cellIdsList->SetId(i, *sIter);
438 extractCells->SetCellList(cellIdsList.Get());
439 extractCells->Update();
440 vtkUnstructuredGrid* extractGrid = extractCells->GetOutput();
442 // Send the extracted grid to the neighbor rank asynchronously
443 CommDataInfo& c = this->Internals->CommData[toRank];
444 if (vtkCommunicator::MarshalDataObject(extractGrid, c.SendBuffer))
446 c.SendLen = c.SendBuffer->GetNumberOfTuples();
448 com->NoBlockSend(&c.SendLen, 1, toRank, UGGCG_SIZE_EXCHANGE_TAG, c.SendReqs[0]);
450 com->NoBlockSend((char*)c.SendBuffer->GetVoidPointer(0), c.SendLen, toRank,
451 UGGCG_DATA_EXCHANGE_TAG, c.SendReqs[1]);
456 //-----------------------------------------------------------------------------
457 // Step 4: receive the ghost cells from the neighbor ranks and merge them
458 // to the local grid.
459 void vtkPUnstructuredGridGhostCellsGenerator::ReceiveAndMergeGhostCells(
460 vtkUnstructuredGrid *output)
462 // We need to compute a rough estimation of the total number of cells and
463 // points for vtkMergeCells
464 vtkIdType totalNbCells = this->Internals->Input->GetNumberOfCells();
465 vtkIdType totalNbPoints = this->Internals->Input->GetNumberOfPoints();
467 vtkMPICommunicator* com =
468 vtkMPICommunicator::SafeDownCast(this->Controller->GetCommunicator());
470 // Browse all neighbor ranks and receive the mesh that contains cells
471 int nbNeighbors = static_cast<int>(this->Internals->NeighborRanksCells.size());
472 std::vector<vtkUnstructuredGridBase*> neighborGrids;
473 neighborGrids.reserve(nbNeighbors);
475 // First create requests to receive the size of the mesh to receive
476 std::map<int, std::set<vtkIdType> >::iterator iter;
477 for (iter = this->Internals->NeighborRanksCells.begin();
478 iter != this->Internals->NeighborRanksCells.end(); ++iter)
480 vtkIdType fromRank = iter->first;
481 CommDataInfo& c = this->Internals->CommData[fromRank];
483 &c.RecvLen, 1, fromRank, UGGCG_SIZE_EXCHANGE_TAG, c.RecvReqs[0]);
486 // Then, once the data length is received, create requests to receive the mesh data
488 while (counter != nbNeighbors)
490 for (iter = this->Internals->NeighborRanksCells.begin();
491 iter != this->Internals->NeighborRanksCells.end(); ++iter)
493 vtkIdType fromRank = iter->first;
494 CommDataInfo& c = this->Internals->CommData[fromRank];
495 if (!c.RecvReqs[0].Test() || c.CommStep != 0)
500 c.RecvBuffer->SetNumberOfValues(c.RecvLen);
502 (char*)c.RecvBuffer->GetVoidPointer(0), c.RecvLen, fromRank,
503 UGGCG_DATA_EXCHANGE_TAG, c.RecvReqs[1]);
508 // Browse all neighbor ranks and receive the mesh that contains cells
509 // that are ghost cells for current rank.
511 while (counter != nbNeighbors)
513 for (iter = this->Internals->NeighborRanksCells.begin();
514 iter != this->Internals->NeighborRanksCells.end(); ++iter)
516 vtkIdType fromRank = iter->first;
517 CommDataInfo& c = this->Internals->CommData[fromRank];
519 if (!c.RecvReqs[1].Test() || c.CommStep != 1)
525 vtkUnstructuredGrid* grid = vtkUnstructuredGrid::New();
526 vtkCommunicator::UnMarshalDataObject(c.RecvBuffer, grid);
527 c.RecvBuffer->Delete();
530 // Flag the received grid elements as ghosts
531 grid->AllocatePointGhostArray();
532 grid->AllocateCellGhostArray();
533 grid->GetPointGhostArray()->FillComponent(0, 1);
534 grid->GetCellGhostArray()->FillComponent(0, 1);
536 // Make sure the global point ids array is tagged accordingly
537 if (this->Internals->InputGlobalPointIds &&
538 !grid->GetPointData()->GetGlobalIds())
540 grid->GetPointData()->SetGlobalIds(grid->GetPointData()->GetArray(
541 this->Internals->InputGlobalPointIds->GetName()));
544 totalNbCells += grid->GetNumberOfCells();
545 totalNbPoints += grid->GetNumberOfPoints();
547 neighborGrids.push_back(grid);
553 // Shallow copy the input grid and initialize the ghost arrays
554 vtkNew<vtkUnstructuredGrid> inputCopy;
555 inputCopy->ShallowCopy(this->Internals->Input);
556 inputCopy->AllocatePointGhostArray();
557 inputCopy->AllocateCellGhostArray();
559 // MergeCells merge input + grids that contains ghost cells to the output grid
560 vtkNew<vtkMergeCells> mergeCells;
561 mergeCells->SetUnstructuredGrid(output);
562 mergeCells->SetTotalNumberOfCells(totalNbCells);
563 mergeCells->SetTotalNumberOfPoints(totalNbPoints);
564 mergeCells->SetTotalNumberOfDataSets(
565 1 + static_cast<int>(this->Internals->NeighborRanksCells.size()));
566 mergeCells->SetUseGlobalIds(this->Internals->InputGlobalPointIds != 0 ? 1 : 0);
567 mergeCells->SetPointMergeTolerance(0.0);
569 // Merge input grid first
570 mergeCells->MergeDataSet(inputCopy.Get());
572 // Then merge ghost grid from neighbor rank
573 for (std::size_t i = 0; i < neighborGrids.size(); i++)
575 mergeCells->MergeDataSet(neighborGrids[i]);
576 neighborGrids[i]->Delete();
579 // Finalize the merged output
580 mergeCells->Finish();