+++ /dev/null
-/*=========================================================================
-
- 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"
-
-#include "vtkCellArray.h"
-#include "vtkCharArray.h"
-#include "vtkDataSetSurfaceFilter.h"
-#include "vtkExtractCells.h"
-#include "vtkIdList.h"
-#include "vtkIdTypeArray.h"
-#include "vtkInformation.h"
-#include "vtkInformationVector.h"
-#include "vtkMergeCells.h"
-#include "vtkMergePoints.h"
-#include "vtkMPICommunicator.h"
-#include "vtkMultiProcessController.h"
-#include "vtkNew.h"
-#include "vtkObjectFactory.h"
-#include "vtkPointData.h"
-#include "vtkPoints.h"
-#include "vtkStreamingDemandDrivenPipeline.h"
-#include "vtkUnstructuredGrid.h"
-
-#include <map>
-#include <set>
-#include <vector>
-
-//----------------------------------------------------------------------------
-// Helpers
-namespace
-{
-template<class T>
-bool AllGatherV(vtkMultiProcessController* controller,
- const T* localV,
- vtkIdType localSize,
- std::vector<T>& globalV,
- std::vector<vtkIdType>& sizes,
- std::vector<vtkIdType>& offsets)
-{
- int nbOfRanks = controller->GetNumberOfProcesses();
- sizes.resize(nbOfRanks);
- int ret = controller->AllGather(&localSize, &sizes[0], 1);
- if (ret == 0)
- {
- vtkErrorWithObjectMacro(controller, << "Communication error!");
- return false;
- }
- vtkIdType count = 0;
- offsets.resize(nbOfRanks);
- for (int i = 0; i < nbOfRanks; i++)
- {
- offsets[i] = count;
- count += sizes[i];
- }
- globalV.resize(count);
- if (count > 0)
- {
- controller->AllGatherV(localSize > 0 ? localV : 0,
- &globalV[0], localSize, &sizes[0], &offsets[0]);
- }
- return true;
-}
-}
-
-//----------------------------------------------------------------------------
-// Internal data structures
-
-// Class to hold asynchronous communication information
-class CommDataInfo
-{
-public:
- CommDataInfo() : SendLen(-1), RecvLen(-1), CommStep(0)
- {
- this->SendBuffer = vtkCharArray::New();
- this->RecvBuffer = vtkCharArray::New();
- }
-
- CommDataInfo(const CommDataInfo& c)
- {
- *this = c;
- if (this->SendBuffer) { this->SendBuffer->Register(0); }
- if (this->RecvBuffer) { this->RecvBuffer->Register(0); }
- }
-
- ~CommDataInfo()
- {
- if (this->SendBuffer) { this->SendBuffer->Delete(); }
- if (this->RecvBuffer) { this->RecvBuffer->Delete(); }
- }
-
- vtkMPICommunicator::Request SendReqs[2];
- vtkMPICommunicator::Request RecvReqs[2];
- vtkCharArray *SendBuffer;
- vtkCharArray *RecvBuffer;
- vtkIdType SendLen;
- vtkIdType RecvLen;
- int CommStep;
-};
-
-// Communication arrays
-struct vtkPUnstructuredGridGhostCellsGenerator::vtkInternals
-{
- // For global ids
- std::map<vtkIdType, vtkIdType> GlobalToLocalPointIdMap;
- std::vector<vtkIdType> AllGlobalIdsOfSurfacePoints;
-
- // For point coordinates
- vtkNew<vtkMergePoints> LocalPoints;
- std::vector<vtkIdType> LocalPointsMap;
- std::vector<double> AllPointsOfSurfacePoints;
-
- std::vector<vtkIdType> AllSizes;
- std::vector<vtkIdType> AllOffsets;
- std::map<int, std::set<vtkIdType> > NeighborRanksCells;
- std::map<int, CommDataInfo> CommData;
- vtkUnstructuredGridBase* Input;
-
- vtkDataArray* InputGlobalPointIds;
- bool UseGlobalPointIds;
-};
-
-static const int UGGCG_SIZE_EXCHANGE_TAG = 9000;
-static const int UGGCG_DATA_EXCHANGE_TAG = 9001;
-static const char* UGGCG_GLOBAL_POINT_IDS = "GlobalNodeIds";
-
-//----------------------------------------------------------------------------
-
-vtkStandardNewMacro(vtkPUnstructuredGridGhostCellsGenerator)
-vtkSetObjectImplementationMacro(
- vtkPUnstructuredGridGhostCellsGenerator, Controller, vtkMultiProcessController);
-
-//----------------------------------------------------------------------------
-vtkPUnstructuredGridGhostCellsGenerator::vtkPUnstructuredGridGhostCellsGenerator()
-{
- this->Controller = NULL;
- this->SetController(vtkMultiProcessController::GetGlobalController());
-
- this->Internals = NULL;
- this->GlobalPointIdsArrayName = NULL;
- this->UseGlobalPointIds = true;
- this->SetGlobalPointIdsArrayName(UGGCG_GLOBAL_POINT_IDS);
-}
-
-//----------------------------------------------------------------------------
-vtkPUnstructuredGridGhostCellsGenerator::~vtkPUnstructuredGridGhostCellsGenerator()
-{
- this->SetController(NULL);
- this->SetGlobalPointIdsArrayName(NULL);
-
- delete this->Internals;
- this->Internals = 0;
-}
-
-//-----------------------------------------------------------------------------
-void vtkPUnstructuredGridGhostCellsGenerator::PrintSelf(ostream& os, vtkIndent indent)
-{
- Superclass::PrintSelf(os, indent);
-}
-
-//-----------------------------------------------------------------------------
-int vtkPUnstructuredGridGhostCellsGenerator::RequestData(
- vtkInformation *vtkNotUsed(request),
- vtkInformationVector **inputVector,
- vtkInformationVector *outputVector)
-{
- // get the info objects
- vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
- vtkInformation *outInfo = outputVector->GetInformationObject(0);
-
- // get the input and output. Input may just have the UnstructuredGridBase
- // interface, but output should be an unstructured grid.
- vtkUnstructuredGridBase *input = vtkUnstructuredGridBase::SafeDownCast(
- inInfo->Get(vtkDataObject::DATA_OBJECT()));
- vtkUnstructuredGrid *output = vtkUnstructuredGrid::SafeDownCast(
- outInfo->Get(vtkDataObject::DATA_OBJECT()));
-
- int ghostLevels = 1;
- //ghostLevels = outInfo->Get(
- // vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS());
-
- if (!this->Controller)
- {
- this->Controller = vtkMultiProcessController::GetGlobalController();
- }
- output->Reset();
-
- this->NumRanks = this->Controller ? this->Controller->GetNumberOfProcesses() : 1;
- this->RankId = this->Controller ? this->Controller->GetLocalProcessId() : 0;
- if (ghostLevels == 0 || !this->Controller || this->NumRanks == 1)
- {
- // Ghost levels are not requested. Nothing to do but pass the dataset.
- if (this->RankId == 0)
- {
- vtkWarningMacro(<< "Ghost cells not requested or not needed.");
- }
- output->ShallowCopy(input);
- return 1;
- }
-
- delete this->Internals;
- this->Internals = new vtkPUnstructuredGridGhostCellsGenerator::vtkInternals();
- this->Internals->Input = input;
-
- vtkPointData *inputPD = input->GetPointData();
- this->Internals->InputGlobalPointIds = inputPD->GetGlobalIds();
-
- if (!this->Internals->InputGlobalPointIds)
- {
- this->Internals->InputGlobalPointIds =
- inputPD->GetArray(this->GlobalPointIdsArrayName);
- inputPD->SetGlobalIds(this->Internals->InputGlobalPointIds);
- }
-
- if (!this->UseGlobalPointIds)
- {
- this->Internals->InputGlobalPointIds = NULL;
- }
-
- this->ExtractAndReduceSurfacePoints();
- this->UpdateProgress(0.3);
-
- this->ComputeSharedPoints();
- this->UpdateProgress(0.6);
-
- this->ExtractAndSendGhostCells();
- this->UpdateProgress(0.8);
-
- this->ReceiveAndMergeGhostCells(output);
- this->UpdateProgress(1.0);
-
- output->GetInformation()->Set(vtkDataObject::DATA_NUMBER_OF_GHOST_LEVELS(), 1);
-
- this->Controller->Barrier();
-
- delete this->Internals;
- this->Internals = 0;
-
- return 1;
-}
-
-//-----------------------------------------------------------------------------
-// Step 1: Extract surface geometry and all reduce global ids of surface points
-void vtkPUnstructuredGridGhostCellsGenerator::ExtractAndReduceSurfacePoints()
-{
- // Extract boundary cells and points with the surface filter
- vtkNew<vtkDataSetSurfaceFilter> surfaceFilter;
- surfaceFilter->SetInputData(this->Internals->Input);
- surfaceFilter->PassThroughPointIdsOn();
- surfaceFilter->Update();
-
- vtkPolyData* surface = surfaceFilter->GetOutput();
- vtkIdType nbSurfacePoints = surface->GetNumberOfPoints();
- vtkCellArray* surfaceCells = surface->GetPolys();
- surfaceCells->InitTraversal();
- vtkIdType npts, *pts;
-
- vtkIdTypeArray* surfaceOriginalPointIds = vtkIdTypeArray::SafeDownCast(
- surface->GetPointData()->GetArray(surfaceFilter->GetOriginalPointIdsName()));
-
- if (this->Internals->InputGlobalPointIds)
- {
- std::vector<vtkIdType> globalIdsOfSurfacePoints;
- globalIdsOfSurfacePoints.reserve(nbSurfacePoints);
-
- // Browse surface cells and save global and local ids of cell points
- while (surfaceCells->GetNextCell(npts, pts))
- {
- for (vtkIdType i = 0; i < npts; i++)
- {
- vtkIdType origPtId = surfaceOriginalPointIds->GetValue(pts[i]);
- vtkIdType globalPtId = static_cast<vtkIdType>(
- this->Internals->InputGlobalPointIds->GetTuple1(origPtId));
-
- if (this->Internals->GlobalToLocalPointIdMap.find(globalPtId) ==
- this->Internals->GlobalToLocalPointIdMap.end())
- {
- this->Internals->GlobalToLocalPointIdMap[globalPtId] = origPtId;
- globalIdsOfSurfacePoints.push_back(globalPtId);
- }
- }
- }
-
- // Now reduce surface point global ids on ALL ranks
- ::AllGatherV(this->Controller, &globalIdsOfSurfacePoints[0],
- globalIdsOfSurfacePoints.size(),
- this->Internals->AllGlobalIdsOfSurfacePoints,
- this->Internals->AllSizes, this->Internals->AllOffsets);
- }
- else
- {
- // We can't use global ids, so we will process point coordinates instead
- vtkPoints* inputPoints = this->Internals->Input->GetPoints();
- vtkNew<vtkPoints> surfacePoints;
- surfacePoints->SetDataTypeToDouble();
- surfacePoints->Allocate(nbSurfacePoints);
- this->Internals->LocalPoints->InitPointInsertion(
- surfacePoints.Get(), surface->GetPoints()->GetBounds());
- this->Internals->LocalPointsMap.reserve(nbSurfacePoints);
-
- // Browse surface cells and push point coordinates to the locator
- while (surfaceCells->GetNextCell(npts, pts))
- {
- for (vtkIdType i = 0; i < npts; i++)
- {
- vtkIdType origPtId = surfaceOriginalPointIds->GetValue(pts[i]);
- double p[3];
- inputPoints->GetPoint(origPtId, p);
- vtkIdType sid;
- if (this->Internals->LocalPoints->InsertUniquePoint(p, sid))
- {
- // New point, save the id of the original grid point id associated
- // to this surface point
- if (static_cast<vtkIdType>(this->Internals->LocalPointsMap.size()) <= sid)
- {
- this->Internals->LocalPointsMap.resize(sid + 1);
- }
- this->Internals->LocalPointsMap[sid] = origPtId;
- }
- }
- }
-
- // Now reduce surface point coordinates on ALL ranks
- ::AllGatherV(this->Controller,
- (double*)surfacePoints->GetVoidPointer(0),
- surfacePoints->GetNumberOfPoints() * 3,
- this->Internals->AllPointsOfSurfacePoints,
- this->Internals->AllSizes, this->Internals->AllOffsets);
- }
-}
-
-//---------------------------------------------------------------------------
-// Step 2: browse global ids/point coordinates of other ranks and check if some
-// 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()
-{
- vtkNew<vtkIdList> cellIdsList;
- for (int i = 0; i < this->NumRanks; i++)
- {
- if (i == this->RankId)
- {
- continue;
- }
- for (vtkIdType j = 0, idx = this->Internals->AllOffsets[i];
- j < this->Internals->AllSizes[i]; j++, idx++)
- {
- vtkIdType localPointId = -1;
- if (this->Internals->InputGlobalPointIds)
- {
- // Check if this point exists locally from its global ids, if so
- // get its local id.
- vtkIdType gid = this->Internals->AllGlobalIdsOfSurfacePoints[idx];
- std::map<vtkIdType, vtkIdType>::iterator iter =
- this->Internals->GlobalToLocalPointIdMap.find(gid);
- if (iter != this->Internals->GlobalToLocalPointIdMap.end())
- {
- localPointId = iter->second;
- }
- }
- else
- {
- // Check if this point exists locally from its coordinates, if so
- // get its local id.
- double *p = &this->Internals->AllPointsOfSurfacePoints[idx];
- localPointId = this->Internals->LocalPoints->IsInsertedPoint(p);
-
- if (localPointId != -1)
- {
- localPointId = this->Internals->LocalPointsMap[localPointId];
- }
- idx += 2; // jump to next coordinates
- j += 2;
- }
-
- if (localPointId != -1)
- {
- // Current rank also has a copy of this global point
- cellIdsList->Reset();
- // 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->NeighborRanksCells[i].insert(cellIdsList->GetId(k));
- }
- }
- }
- }
-
- // Release memory of all reduced arrays
- this->Internals->AllGlobalIdsOfSurfacePoints.resize(0);
- this->Internals->AllPointsOfSurfacePoints.resize(0);
- this->Internals->AllSizes.resize(0);
- this->Internals->AllOffsets.resize(0);
- // 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
-void vtkPUnstructuredGridGhostCellsGenerator::ExtractAndSendGhostCells()
-{
- vtkNew<vtkIdList> cellIdsList;
- vtkNew<vtkExtractCells> extractCells;
- extractCells->SetInputData(this->Internals->Input);
- std::map<int, std::set<vtkIdType> >::iterator iter =
- this->Internals->NeighborRanksCells.begin();
-
- vtkMPICommunicator* com =
- vtkMPICommunicator::SafeDownCast(this->Controller->GetCommunicator());
-
- // Browse all neighbor ranks and extract the cells connected to the points we share
- for (; iter != this->Internals->NeighborRanksCells.end(); ++iter)
- {
- int toRank = iter->first;
- std::set<vtkIdType>& cellsToShare = iter->second;
- cellIdsList->SetNumberOfIds(cellsToShare.size());
- std::set<vtkIdType>::iterator sIter = cellsToShare.begin();
- for (vtkIdType i = 0; sIter != cellsToShare.end(); ++sIter, i++)
- {
- cellIdsList->SetId(i, *sIter);
- }
- extractCells->SetCellList(cellIdsList.Get());
- extractCells->Update();
- vtkUnstructuredGrid* extractGrid = extractCells->GetOutput();
-
- // Send the extracted grid to the neighbor rank asynchronously
- CommDataInfo& c = this->Internals->CommData[toRank];
- if (vtkCommunicator::MarshalDataObject(extractGrid, c.SendBuffer))
- {
- c.SendLen = c.SendBuffer->GetNumberOfTuples();
- // Send data length
- com->NoBlockSend(&c.SendLen, 1, toRank, UGGCG_SIZE_EXCHANGE_TAG, c.SendReqs[0]);
- // Send raw data
- com->NoBlockSend((char*)c.SendBuffer->GetVoidPointer(0), c.SendLen, toRank,
- UGGCG_DATA_EXCHANGE_TAG, c.SendReqs[1]);
- }
- }
-}
-
-//-----------------------------------------------------------------------------
-// Step 4: receive the ghost cells from the neighbor ranks and merge them
-// to the local grid.
-void vtkPUnstructuredGridGhostCellsGenerator::ReceiveAndMergeGhostCells(
- vtkUnstructuredGrid *output)
-{
- // We need to compute a rough estimation of the total number of cells and
- // points for vtkMergeCells
- vtkIdType totalNbCells = this->Internals->Input->GetNumberOfCells();
- vtkIdType totalNbPoints = this->Internals->Input->GetNumberOfPoints();
-
- vtkMPICommunicator* com =
- vtkMPICommunicator::SafeDownCast(this->Controller->GetCommunicator());
-
- // Browse all neighbor ranks and receive the mesh that contains cells
- int nbNeighbors = static_cast<int>(this->Internals->NeighborRanksCells.size());
- std::vector<vtkUnstructuredGridBase*> neighborGrids;
- neighborGrids.reserve(nbNeighbors);
-
- // First create requests to receive the size of the mesh to receive
- std::map<int, std::set<vtkIdType> >::iterator iter;
- for (iter = this->Internals->NeighborRanksCells.begin();
- iter != this->Internals->NeighborRanksCells.end(); ++iter)
- {
- vtkIdType fromRank = iter->first;
- CommDataInfo& c = this->Internals->CommData[fromRank];
- com->NoBlockReceive(
- &c.RecvLen, 1, fromRank, UGGCG_SIZE_EXCHANGE_TAG, c.RecvReqs[0]);
- }
-
- // Then, once the data length is received, create requests to receive the mesh data
- int counter = 0;
- while (counter != nbNeighbors)
- {
- for (iter = this->Internals->NeighborRanksCells.begin();
- iter != this->Internals->NeighborRanksCells.end(); ++iter)
- {
- vtkIdType fromRank = iter->first;
- CommDataInfo& c = this->Internals->CommData[fromRank];
- if (!c.RecvReqs[0].Test() || c.CommStep != 0)
- {
- continue;
- }
- c.CommStep = 1;
- c.RecvBuffer->SetNumberOfValues(c.RecvLen);
- com->NoBlockReceive(
- (char*)c.RecvBuffer->GetVoidPointer(0), c.RecvLen, fromRank,
- UGGCG_DATA_EXCHANGE_TAG, c.RecvReqs[1]);
- counter++;
- }
- }
-
- // Browse all neighbor ranks and receive the mesh that contains cells
- // that are ghost cells for current rank.
- counter = 0;
- while (counter != nbNeighbors)
- {
- for (iter = this->Internals->NeighborRanksCells.begin();
- iter != this->Internals->NeighborRanksCells.end(); ++iter)
- {
- vtkIdType fromRank = iter->first;
- CommDataInfo& c = this->Internals->CommData[fromRank];
-
- if (!c.RecvReqs[1].Test() || c.CommStep != 1)
- {
- continue;
- }
-
- c.CommStep = 2;
- vtkUnstructuredGrid* grid = vtkUnstructuredGrid::New();
- vtkCommunicator::UnMarshalDataObject(c.RecvBuffer, grid);
- c.RecvBuffer->Delete();
- c.RecvBuffer = NULL;
-
- // Flag the received grid elements as ghosts
- grid->AllocatePointGhostArray();
- grid->AllocateCellGhostArray();
- grid->GetPointGhostArray()->FillComponent(0, 1);
- grid->GetCellGhostArray()->FillComponent(0, 1);
-
- // Make sure the global point ids array is tagged accordingly
- if (this->Internals->InputGlobalPointIds &&
- !grid->GetPointData()->GetGlobalIds())
- {
- grid->GetPointData()->SetGlobalIds(grid->GetPointData()->GetArray(
- this->Internals->InputGlobalPointIds->GetName()));
- }
-
- totalNbCells += grid->GetNumberOfCells();
- totalNbPoints += grid->GetNumberOfPoints();
-
- neighborGrids.push_back(grid);
-
- counter++;
- }
- }
-
- // Shallow copy the input grid and initialize the ghost arrays
- vtkNew<vtkUnstructuredGrid> inputCopy;
- inputCopy->ShallowCopy(this->Internals->Input);
- inputCopy->AllocatePointGhostArray();
- inputCopy->AllocateCellGhostArray();
-
- // MergeCells merge input + grids that contains ghost cells to the output grid
- vtkNew<vtkMergeCells> mergeCells;
- mergeCells->SetUnstructuredGrid(output);
- mergeCells->SetTotalNumberOfCells(totalNbCells);
- mergeCells->SetTotalNumberOfPoints(totalNbPoints);
- mergeCells->SetTotalNumberOfDataSets(
- 1 + static_cast<int>(this->Internals->NeighborRanksCells.size()));
- mergeCells->SetUseGlobalIds(this->Internals->InputGlobalPointIds != 0 ? 1 : 0);
- mergeCells->SetPointMergeTolerance(0.0);
-
- // Merge input grid first
- mergeCells->MergeDataSet(inputCopy.Get());
-
- // Then merge ghost grid from neighbor rank
- for (std::size_t i = 0; i < neighborGrids.size(); i++)
- {
- mergeCells->MergeDataSet(neighborGrids[i]);
- neighborGrids[i]->Delete();
- }
-
- // Finalize the merged output
- mergeCells->Finish();
-}
