From 96cb48791961d788a63c1334410917a1c63b70c1 Mon Sep 17 00:00:00 2001 From: Joachim Pouderoux Date: Thu, 9 Aug 2018 15:30:01 -0400 Subject: [PATCH] Clean & refactorize --- Domains/ParallelChemistry/module.cmake | 1 + .../vtkPSimpleBondPerceiver.cxx | 4 +- .../vtkDistributedPointCloudFilter.cxx | 319 +++++++++++++++--- .../vtkDistributedPointCloudFilter.h | 19 +- Parallel/MPI/vtkMPIUtilities.cxx | 239 ------------- Parallel/MPI/vtkMPIUtilities.h | 11 - 6 files changed, 275 insertions(+), 318 deletions(-) diff --git a/Domains/ParallelChemistry/module.cmake b/Domains/ParallelChemistry/module.cmake index 570ea55726c..7a6fd8e5a9c 100644 --- a/Domains/ParallelChemistry/module.cmake +++ b/Domains/ParallelChemistry/module.cmake @@ -12,4 +12,5 @@ vtk_module(vtkDomainsParallelChemistry PRIVATE_DEPENDS vtkCommonCore vtkParallelMPI + vtkFiltersParallelMPI ) diff --git a/Domains/ParallelChemistry/vtkPSimpleBondPerceiver.cxx b/Domains/ParallelChemistry/vtkPSimpleBondPerceiver.cxx index 63038e842cc..ae091c537c5 100644 --- a/Domains/ParallelChemistry/vtkPSimpleBondPerceiver.cxx +++ b/Domains/ParallelChemistry/vtkPSimpleBondPerceiver.cxx @@ -16,9 +16,9 @@ #include "vtkPSimpleBondPerceiver.h" #include "vtkDataSetAttributes.h" +#include "vtkDistributedPointCloudFilter.h" #include "vtkFloatArray.h" #include "vtkMPIController.h" -#include "vtkMPIUtilities.h" #include "vtkMolecule.h" #include "vtkNew.h" #include "vtkObjectFactory.h" @@ -75,7 +75,7 @@ bool vtkPSimpleBondPerceiver::CreateGhosts(vtkMolecule* molecule) dataArray->DeepCopy(molecule->GetVertexData()); vtkNew outputPoly; - vtkMPIUtilities::GetPointsInsideBounds( + vtkDistributedPointCloudFilter::GetPointsInsideBounds( controller, inputPoly.Get(), outputPoly.Get(), outterBounds); molecule->Initialize( diff --git a/Filters/ParallelMPI/vtkDistributedPointCloudFilter.cxx b/Filters/ParallelMPI/vtkDistributedPointCloudFilter.cxx index 4b9254520f5..4aeb33a7c3a 100644 --- a/Filters/ParallelMPI/vtkDistributedPointCloudFilter.cxx +++ b/Filters/ParallelMPI/vtkDistributedPointCloudFilter.cxx @@ -15,16 +15,19 @@ #include "vtkDistributedPointCloudFilter.h" +#include "vtkCharArray.h" #include "vtkDataArray.h" #include "vtkInformation.h" #include "vtkMPIController.h" #include "vtkMPIUtilities.h" #include "vtkNew.h" #include "vtkObjectFactory.h" +#include "vtkOctreePointLocator.h" #include "vtkPointData.h" #include "vtkPointSet.h" #include "vtkPoints.h" #include "vtkPolyData.h" +#include "vtkSmartPointer.h" #include @@ -48,7 +51,7 @@ vtkDistributedPointCloudFilter::~vtkDistributedPointCloudFilter() } //---------------------------------------------------------------------------- -int vtkDistributedPointCloudFilter::FillOutputPortInformation(int port, vtkInformation *info) +int vtkDistributedPointCloudFilter::FillOutputPortInformation(int, vtkInformation *info) { info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkPolyData"); return 1; @@ -81,62 +84,51 @@ int vtkDistributedPointCloudFilter::RequestData(vtkInformation *, return 1; } - std::vector kdTreeRounds; - if (this->InitializeKdTree(kdTreeRounds)) + std::vector subControllersTree; + if (this->InitializeKdTree(subControllersTree)) { double bounds[6]; - this->OptimizeBoundingBox(kdTreeRounds, input, bounds); - vtkMPIUtilities::GetPointsInsideBounds(controller, input, output, bounds); + this->OptimizeBoundingBox(subControllersTree, input, bounds); + this->GetPointsInsideBounds(controller, input, output, bounds); } else { - vtkErrorMacro("SubCommunicators are not correctly initialized, no distribution performed"); + vtkErrorMacro("Sub-communicators are not correctly initialized, no distribution performed"); return 0; } // Destroy the kdtree controllers - for (const auto &round : kdTreeRounds) + for (auto *c : subControllersTree) { - round.controller->Delete(); + c->Delete(); } return 1; } //---------------------------------------------------------------------------- -bool vtkDistributedPointCloudFilter::InitializeKdTree(std::vector &kdTreeRounds) +bool vtkDistributedPointCloudFilter::InitializeKdTree(std::vector &kdTreeRounds) { this->Controller->Register(this); - kdTreeRounds.push_back({ vtkMPIController::SafeDownCast(this->Controller), - this->Controller->GetNumberOfProcesses(), - this->Controller->GetLocalProcessId() }); + kdTreeRounds.push_back(vtkMPIController::SafeDownCast(this->Controller)); int index = 0; - while (kdTreeRounds[index].np > 2) + while (kdTreeRounds[index]->GetNumberOfProcesses() > 2) { - int split = kdTreeRounds[index].np / 2; - int color, key; + vtkMPIController* ctrl = kdTreeRounds[index]; + int roundRank = ctrl->GetLocalProcessId(); + int split = ctrl->GetNumberOfProcesses() / 2; - if (kdTreeRounds[index].rank < split) - { - color = 0; - key = kdTreeRounds[index].rank; - } - else - { - color = 1; - key = kdTreeRounds[index].rank - split; - } + int color = (roundRank < split) ? 0 : 1; + int key = roundRank - ((roundRank < split) ? 0 : split); - vtkMPIController* ctrl = - kdTreeRounds[index].controller->PartitionController(color, key); - if (ctrl == nullptr) + vtkMPIController *subCtrl = ctrl->PartitionController(color, key); + if (subCtrl == nullptr) { break; } + kdTreeRounds.push_back(subCtrl); ++index; - kdTreeRounds.push_back( - { ctrl, ctrl->GetNumberOfProcesses(), ctrl->GetLocalProcessId() }); } return true; @@ -144,7 +136,7 @@ bool vtkDistributedPointCloudFilter::InitializeKdTree(std::vector &kdTreeRounds, vtkPointSet* pointCloud, double regionBounds[6]) + std::vector &kdTreeRounds, vtkPointSet *pointCloud, double regionBounds[6]) { if (kdTreeRounds.empty()) { @@ -195,7 +187,7 @@ bool vtkDistributedPointCloudFilter::OptimizeBoundingBox( const int histSize = HISTOGRAM_SIZE; std::vector histogram(histSize, 0); std::vector histsum(histSize, 0); - std::vector pointExchangeCount(kdTreeRounds[0].np, 0); + std::vector pointExchangeCount(kdTreeRounds[0]->GetNumberOfProcesses(), 0); std::vector pts; if (numPts > 0) @@ -209,14 +201,20 @@ bool vtkDistributedPointCloudFilter::OptimizeBoundingBox( for (size_t round = 0; round < kdTreeRounds.size(); round++) { - if (kdTreeRounds[round].np == 1) + vtkMPIController* kdTreeRound = kdTreeRounds[round]; + vtkMPICommunicator* roundComm = + vtkMPICommunicator::SafeDownCast(kdTreeRound->GetCommunicator()); + int roundNumRanks = kdTreeRound->GetNumberOfProcesses(); + int roundRank = kdTreeRound->GetLocalProcessId(); + + if (roundNumRanks == 1) { continue; } - kdTreeRounds[round].controller->GetCommunicator()->AllReduceVoidArray( + roundComm->AllReduceVoidArray( localLowerBound, currentGroupLowerBound, 3, VTK_DOUBLE, vtkCommunicator::MIN_OP); - kdTreeRounds[round].controller->GetCommunicator()->AllReduceVoidArray( + roundComm->AllReduceVoidArray( localUpperBound, currentGroupUpperBound, 3, VTK_DOUBLE, vtkCommunicator::MAX_OP); // ---------------------------------------- @@ -247,21 +245,21 @@ bool vtkDistributedPointCloudFilter::OptimizeBoundingBox( // ---------------------------------------- // 3. reduction across round participants to get global histogram - kdTreeRounds[round].controller->GetCommunicator()->ReduceVoidArray( + roundComm->ReduceVoidArray( histogram.data(), histsum.data(), histSize, VTK_INT, vtkCommunicator::SUM_OP, 0); vtkIdType totalNumpts = numPts; - kdTreeRounds[round].controller->GetCommunicator()->ReduceVoidArray( + roundComm->ReduceVoidArray( &numPts, &totalNumpts, 1, VTK_ID_TYPE, vtkCommunicator::SUM_OP, 0); // ---------------------------------------- // 4. process 0 of sub group computes cut position and broadcast it to others int cutpos = 0; - if (kdTreeRounds[round].rank == 0) + if (roundRank == 0) { // ratio of left part over total. it won't be 2 for un-even number of participant processors - double ratio = static_cast(kdTreeRounds[round].np) / - static_cast(kdTreeRounds[round].np / 2); + double ratio = static_cast(roundNumRanks) / + static_cast(roundNumRanks / 2); int i = 1; for (; i < histSize; i++) { @@ -273,13 +271,12 @@ bool vtkDistributedPointCloudFilter::OptimizeBoundingBox( } cutpos = i; } - kdTreeRounds[round].controller->GetCommunicator()->BroadcastVoidArray( - &cutpos, 1, VTK_INT, 0); + roundComm->BroadcastVoidArray(&cutpos, 1, VTK_INT, 0); // ---------------------------------------- // 5. classify points // which half of the group i belongs to. true = left, false = right - bool side = kdTreeRounds[round].rank < (kdTreeRounds[round].np / 2); + bool side = roundRank < (roundNumRanks / 2); std::vector partnerPts(3 * numPts); vtkIdType newNumPts = 0; vtkIdType partnerNumPts = 0; @@ -303,8 +300,7 @@ bool vtkDistributedPointCloudFilter::OptimizeBoundingBox( } } - kdTreeRounds[round].controller->GetCommunicator()->AllGatherVoidArray( - &partnerNumPts, pointExchangeCount.data(), 1, VTK_ID_TYPE); + roundComm->AllGatherVoidArray(&partnerNumPts, pointExchangeCount.data(), 1, VTK_ID_TYPE); // ---------------------------------------- // 6. exchange points @@ -312,18 +308,18 @@ bool vtkDistributedPointCloudFilter::OptimizeBoundingBox( int partner; if (side) { - partner = kdTreeRounds[round].rank + (kdTreeRounds[round].np / 2); + partner = roundRank + (roundNumRanks / 2); } else { - partner = kdTreeRounds[round].rank - (kdTreeRounds[round].np / 2); + partner = roundRank - (roundNumRanks / 2); } vtkIdType toReceive = pointExchangeCount[partner]; - bool even = (kdTreeRounds[round].np % 2) == 0; + bool even = (roundNumRanks % 2) == 0; // If non even number of processes: last one send to 0 and receive nothing. - if (!even && kdTreeRounds[round].rank == kdTreeRounds[round].np - 1) + if (!even && roundRank == roundNumRanks - 1) { partner = 0; toReceive = 0; @@ -331,30 +327,28 @@ bool vtkDistributedPointCloudFilter::OptimizeBoundingBox( pts.resize(3 * (newNumPts + toReceive)); - vtkMPICommunicator *com = - vtkMPICommunicator::SafeDownCast(kdTreeRounds[round].controller->GetCommunicator()); vtkMPICommunicator::Request request; const int EXCHANGE_POINT_TAG = 524821; if (partnerNumPts > 0) { - com->NoBlockSend(partnerPts.data(), static_cast(3 * partnerNumPts), partner, EXCHANGE_POINT_TAG, request); + roundComm->NoBlockSend(partnerPts.data(), static_cast(3 * partnerNumPts), partner, EXCHANGE_POINT_TAG, request); } if (toReceive > 0) { - com->ReceiveVoidArray(&pts[3 * newNumPts], static_cast(3 * toReceive), VTK_DOUBLE, partner, EXCHANGE_POINT_TAG); + roundComm->ReceiveVoidArray(&pts[3 * newNumPts], static_cast(3 * toReceive), VTK_DOUBLE, partner, EXCHANGE_POINT_TAG); } // If non even number of processes: 0 receive from the last one. - if (!even && kdTreeRounds[round].rank == 0) + if (!even && roundRank == 0) { newNumPts += toReceive; - partner = kdTreeRounds[round].np - 1; + partner = roundNumRanks - 1; toReceive = pointExchangeCount[partner]; pts.resize(3 * (newNumPts + toReceive)); if (toReceive > 0) { - com->ReceiveVoidArray(&pts[3 * newNumPts], 3 * toReceive, VTK_DOUBLE, partner, partner); + roundComm->ReceiveVoidArray(&pts[3 * newNumPts], 3 * toReceive, VTK_DOUBLE, partner, partner); } } @@ -396,6 +390,219 @@ bool vtkDistributedPointCloudFilter::OptimizeBoundingBox( return true; } +//---------------------------------------------------------------------------- +void vtkDistributedPointCloudFilter::GetPointsInsideBounds(vtkMPIController *controller, + vtkPointSet *input, vtkPointSet *output, const double outterBounds[6]) +{ + vtkMPICommunicator *com = vtkMPICommunicator::SafeDownCast(controller->GetCommunicator()); + int np = com ? com->GetNumberOfProcesses() : 1; + int rank = com ? com->GetLocalProcessId() : 0; + + if (!com || np == 1) + { + output->ShallowCopy(input); + return; + } + + // round bounds to the nearest float value because locator use float internally. + // Otherwise, points that are exactly on the bounds may be wrongly considered as outside + // because of the cast. + double localOutterBounds[6]; + for (int i = 0; i < 3; i++) + { + localOutterBounds[2 * i] = + std::nextafter(static_cast(outterBounds[2 * i]), static_cast(outterBounds[2 * i]) - 1); + localOutterBounds[2 * i + 1] = + std::nextafter(static_cast(outterBounds[2 * i + 1]), static_cast(outterBounds[2 * i + 1]) + 1); + } + + bool emptyData = input->GetNumberOfPoints() == 0; + + std::vector allOutterBounds(np * 6); + com->AllGather(localOutterBounds, allOutterBounds.data(), 6); + + // size in bytes of messages to be sent to other processes + std::vector messagesSize(np); + + // number of points in messages to be sent to other processes + std::vector messagePointCount(np); + + // array of point ids + vtkNew idArray; + std::vector> dataToSend; + dataToSend.resize(np); + + // we will need a locator to search points inside each processor assigned regions + vtkNew locator; + + if (!emptyData) + { + vtkNew inputPolyData; + inputPolyData->SetPoints(input->GetPoints()); + locator->SetDataSet(inputPolyData.Get()); + locator->BuildLocator(); + } + + // 1st step: define messages to send to each processor (including itself) + // with polydata containing closest points to local data bounding box + for (int partner = 0; partner < np; partner++) + { + idArray->SetNumberOfTuples(0); + vtkIdType nPoints = 0; + vtkIdType *ids = nullptr; + if (!emptyData) + { + double *nbounds = &allOutterBounds[partner * 6]; + locator->FindPointsInArea(nbounds, idArray.Get()); + nPoints = idArray->GetNumberOfTuples(); + ids = idArray->GetPointer(0); + } + + vtkNew pointCloudToSend; + vtkNew pointsToSend; + pointsToSend->SetNumberOfPoints(nPoints); + pointCloudToSend->SetPoints(pointsToSend.Get()); + + vtkPointData *pointsToSendPointData = pointCloudToSend->GetPointData(); + pointsToSendPointData->CopyAllocate(input->GetPointData(), nPoints); + + for (vtkIdType i = 0; i < nPoints; i++) + { + pointsToSend->SetPoint(i, input->GetPoint(ids[i])); + pointsToSendPointData->CopyData(input->GetPointData(), ids[i], i); + } + + // flatten (marshal) point coordinates & data to a raw byte array + messagePointCount[partner] = nPoints; + dataToSend[partner] = vtkSmartPointer::New(); + vtkCommunicator::MarshalDataObject(pointCloudToSend.Get(), dataToSend[partner]); + messagesSize[partner] = dataToSend[partner]->GetNumberOfValues(); + } + + std::vector> dataToReceive(np); + std::vector receiveRequests(np); + + // Calculate size of messages to receive + std::vector receiveSize(np); + std::vector receivePointCount(np); + + for (int i = 0; i < np; i++) + { + com->Gather(messagesSize.data() + i, receiveSize.data(), 1, i); + com->Gather(messagePointCount.data() + i, receivePointCount.data(), 1, i); + } + + // Starting asynchronous receptions + int nReceive = 0; + vtkIdType totalPointsToReceive = 0; + for (int round = 0; round < np - 1; round++) + { + int partner = (rank + round + 1) % np; + if (receiveSize[partner] > 0) + { + dataToReceive[partner] = vtkSmartPointer::New(); + com->NoBlockReceive(dataToReceive[partner]->WritePointer(0, receiveSize[partner]), + receiveSize[partner], partner, 0, receiveRequests[partner]); + totalPointsToReceive += receivePointCount[partner]; + nReceive++; + } + } + + // local sending/receipt is just a pointer assignment + dataToReceive[rank] = dataToSend[rank]; + dataToSend[rank] = nullptr; + if (receiveSize[rank] > 0) + { + totalPointsToReceive += receivePointCount[rank]; + nReceive++; + } + + // Starting asynchronous sends + std::vector sendRequests(np); + for (int round = 0; round < np - 1; round++) + { + int partner = (rank + round + 1) % np; + if (messagesSize[partner] > 0) + { + com->NoBlockSend(dataToSend[partner]->GetPointer(0), messagesSize[partner], + partner, 0, sendRequests[partner]); + } + } + + // sum of received points from the different processors + vtkIdType totalPoints = 0; + vtkPointData *outputPointData = output->GetPointData(); + outputPointData->SetNumberOfTuples(totalPointsToReceive); + + while (nReceive > 0) + { + for (int round = 0; round < np; round++) + { + int partner = (rank + round) % np; + if ((partner == rank || receiveRequests[partner].Test() == 1) && receiveSize[partner] > 0) + { + vtkNew receivedPointCloud; + vtkCommunicator::UnMarshalDataObject(dataToReceive[partner], receivedPointCloud.Get()); + + dataToReceive[partner] = nullptr; + vtkIdType nbReceivedPoints = receivedPointCloud->GetNumberOfPoints(); + vtkPointData *receivedPointData = receivedPointCloud->GetPointData(); + vtkPoints *receivedPoints = receivedPointCloud->GetPoints(); + vtkPoints *outputPoints = output->GetPoints(); + if (!outputPoints) + { + vtkNew points; + outputPoints = points.Get(); + output->SetPoints(outputPoints); + } + vtkIdType outputNbPts = outputPoints->GetNumberOfPoints(); + outputPoints->Resize(outputNbPts + nbReceivedPoints); + for (vtkIdType i = 0; i < nbReceivedPoints; i++) + { + outputPoints->InsertNextPoint(receivedPoints->GetPoint(i)); + } + int nbArray = receivedPointData->GetNumberOfArrays(); + + for (int a = 0; a < nbArray; a++) + { + vtkAbstractArray *fromArray = receivedPointData->GetAbstractArray(a); + if (fromArray) + { + vtkAbstractArray *toArray = outputPointData->GetAbstractArray(fromArray->GetName()); + if (!toArray) + { + toArray = fromArray->NewInstance(); + toArray->SetName(fromArray->GetName()); + toArray->SetNumberOfComponents(fromArray->GetNumberOfComponents()); + toArray->SetNumberOfTuples(totalPointsToReceive); + outputPointData->AddArray(toArray); + toArray->Delete(); + } + + for (vtkIdType i = 0; i < nbReceivedPoints; i++) + { + toArray->SetTuple(totalPoints + i, i, fromArray); + } + } + } + totalPoints += nbReceivedPoints; + nReceive--; + receiveSize[partner] = 0; + } + } + } + + // we wait for sent messages to be received before deleting them + for (int round = 0; round < np - 1; round++) + { + int partner = (rank + round + 1) % np; + if (messagesSize[partner] > 0) + { + sendRequests[partner].Wait(); + } + } +} + //---------------------------------------------------------------------------- void vtkDistributedPointCloudFilter::PrintSelf(ostream& os, vtkIndent indent) { diff --git a/Filters/ParallelMPI/vtkDistributedPointCloudFilter.h b/Filters/ParallelMPI/vtkDistributedPointCloudFilter.h index 78151069a54..884c23618ef 100644 --- a/Filters/ParallelMPI/vtkDistributedPointCloudFilter.h +++ b/Filters/ParallelMPI/vtkDistributedPointCloudFilter.h @@ -39,7 +39,6 @@ class vtkMPIController; class vtkMultiProcessController; -class vtkPointSet; class VTKFILTERSPARALLELMPI_EXPORT vtkDistributedPointCloudFilter : public vtkPointSetAlgorithm { @@ -56,6 +55,13 @@ class VTKFILTERSPARALLELMPI_EXPORT vtkDistributedPointCloudFilter : public vtkPo vtkGetObjectMacro(Controller, vtkMultiProcessController); //@} + /** + * Get the points that are inside innerBounds and put them in output DataSet. + * Ask other MPI ranks for their corresponding points. + */ + static void GetPointsInsideBounds(vtkMPIController*, + vtkPointSet *input, vtkPointSet *output, const double innerBounds[6]); + protected: vtkDistributedPointCloudFilter(); ~vtkDistributedPointCloudFilter() override; @@ -64,13 +70,6 @@ class VTKFILTERSPARALLELMPI_EXPORT vtkDistributedPointCloudFilter : public vtkPo int RequestData(vtkInformation*, vtkInformationVector**, vtkInformationVector*) override; - struct KdTreeBuildRound - { - vtkMPIController *controller; - int np; - int rank; - }; - /** * Optimize bounding box following this rules: * - no intersection of bounding box of different MPI nodes @@ -79,14 +78,14 @@ class VTKFILTERSPARALLELMPI_EXPORT vtkDistributedPointCloudFilter : public vtkPo * Return false if input pointSet is nullptr or if no communicator was found. * Return true otherwise. */ - bool OptimizeBoundingBox(std::vector&, vtkPointSet*, double bounds[6]); + bool OptimizeBoundingBox(std::vector &, vtkPointSet *, double bounds[6]); /** * Initialize KdTreeRound: creates subControllers from Controller. * Delete old values if any. * Return false if KdTree cannot be initialized. */ - bool InitializeKdTree(std::vector&); + bool InitializeKdTree(std::vector &); private: vtkDistributedPointCloudFilter(const vtkDistributedPointCloudFilter&) = delete; diff --git a/Parallel/MPI/vtkMPIUtilities.cxx b/Parallel/MPI/vtkMPIUtilities.cxx index 75ea2929ccc..1866a5718e9 100644 --- a/Parallel/MPI/vtkMPIUtilities.cxx +++ b/Parallel/MPI/vtkMPIUtilities.cxx @@ -15,21 +15,11 @@ #include "vtkMPIUtilities.h" // VTK includes -#include "vtkCharArray.h" -#include "vtkIdTypeArray.h" #include "vtkMPICommunicator.h" #include "vtkMPIController.h" -#include "vtkNew.h" -#include "vtkOctreePointLocator.h" -#include "vtkPointData.h" -#include "vtkPointSet.h" -#include "vtkPoints.h" -#include "vtkPolyData.h" -#include "vtkSmartPointer.h" // C/C++ includes #include -#include #include #include @@ -122,233 +112,4 @@ void SynchronizedPrintf(vtkMPIController* comm, const char* format, ...) comm->Barrier(); } -void GetPointsInsideBounds(vtkMPIController* controller, - vtkPointSet* input, - vtkPointSet* output, - const double outterBounds[6]) -{ - vtkMPICommunicator* com = vtkMPICommunicator::SafeDownCast(controller->GetCommunicator()); - - if (!com) - { - return; - } - - int np = com->GetNumberOfProcesses(); - int rank = com->GetLocalProcessId(); - - if (np == 1) - { - output->ShallowCopy(input); - return; - } - - double localOutterBounds[6]; - - // round bounds to the nearest float value because locator use float internally. - // Otherwise, points that are exactly on the bounds may be wrongly considered as outside - // because of the cast. - for (int i = 0; i < 3; i++) - { - localOutterBounds[2 * i] = - std::nextafter((float)outterBounds[2 * i], (float)outterBounds[2 * i] - 1); - localOutterBounds[2 * i + 1] = - std::nextafter((float)outterBounds[2 * i + 1], (float)outterBounds[2 * i + 1] + 1); - } - - bool emptyData = input->GetNumberOfPoints() == 0; - - std::vector allOutterBounds(np * 6, 0); - com->AllGather(localOutterBounds, allOutterBounds.data(), 6); - - // size in bytes of messages to be sent to other processes - std::vector messagesSize(np, 0); - - // number of points in messages to be sent to other processes - std::vector messagePointCount(np, 0); - - // array of point ids - vtkNew idArray; - std::vector > dataToSend; - dataToSend.resize(np); - - // we will need a locator to search points inside each processor assigned regions - vtkNew locator; - - if (!emptyData) - { - vtkNew inputPolyData; - inputPolyData->SetPoints(input->GetPoints()); - locator->SetDataSet(inputPolyData.Get()); - locator->BuildLocator(); - } - - // 1st step: define messages to send to each processor (including itself) - // with polydata containing closest points to local data bounding box - for (int partner = 0; partner < np; partner++) - { - idArray->SetNumberOfTuples(0); - vtkIdType nPoints = 0; - vtkIdType* ids = nullptr; - if (!emptyData) - { - double* nbounds = &allOutterBounds[partner * 6]; - locator->FindPointsInArea(nbounds, idArray.Get()); - nPoints = idArray->GetNumberOfTuples(); - ids = idArray->GetPointer(0); - } - - vtkNew pointCloudToSend; - vtkNew pointsToSend; - pointsToSend->SetNumberOfPoints(nPoints); - - vtkPointData* pointsToSendPointData = pointCloudToSend->GetPointData(); - pointsToSendPointData->CopyAllocate(input->GetPointData(), nPoints); - - for (vtkIdType i = 0; i < nPoints; i++) - { - pointsToSend->SetPoint(i, input->GetPoint(ids[i])); - pointsToSendPointData->CopyData(input->GetPointData(), ids[i], i); - } - pointCloudToSend->SetPoints(pointsToSend.Get()); - - // flatten point data to byte array - messagePointCount[partner] = nPoints; - dataToSend[partner] = vtkSmartPointer::New(); - vtkCommunicator::MarshalDataObject(pointCloudToSend.Get(), dataToSend[partner]); - messagesSize[partner] = dataToSend[partner]->GetNumberOfTuples(); - } - - std::vector > dataToReceive; - dataToReceive.resize(np); - - std::vector receiveRequests; - receiveRequests.resize(np); - - // Calculate size of messages to receive - std::vector receiveSize(np, 0); - std::vector receivePointCount(np, 0); - - for (int i = 0; i < np; i++) - { - com->Gather(messagesSize.data() + i, receiveSize.data(), 1, i); - com->Gather(messagePointCount.data() + i, receivePointCount.data(), 1, i); - } - - // Starting asynchronous receptions - int nReceive = 0; - vtkIdType totalPointsToReceive = 0; - for (int round = 0; round < np - 1; round++) - { - int partner = (rank + round + 1) % np; - if (receiveSize[partner] > 0) - { - ++nReceive; - dataToReceive[partner] = vtkSmartPointer::New(); - com->NoBlockReceive(dataToReceive[partner]->WritePointer(0, receiveSize[partner]), - receiveSize[partner], - partner, - 0, - receiveRequests[partner]); - totalPointsToReceive += receivePointCount[partner]; - } - } - - // local sending/receipt is just a pointer assignment - dataToReceive[rank] = dataToSend[rank]; - dataToSend[rank] = nullptr; - if (receiveSize[rank] > 0) - { - ++nReceive; - totalPointsToReceive += receivePointCount[rank]; - } - - // Starting asynchronous sends - std::vector sendRequests; - sendRequests.resize(np); - for (int round = 0; round < np - 1; round++) - { - int partner = (rank + round + 1) % np; - if (messagesSize[partner] > 0) - { - com->NoBlockSend(dataToSend[partner]->GetPointer(0), - messagesSize[partner], - partner, - 0, - sendRequests[partner]); - } - } - - // sum of received points from the different processors - vtkIdType totalPoints = 0; - vtkPointData* outputPointData = output->GetPointData(); - outputPointData->SetNumberOfTuples(totalPointsToReceive); - - while (nReceive > 0) - { - for (int round = 0; round < np; round++) - { - int partner = (rank + round) % np; - if ((partner == rank || receiveRequests[partner].Test() == 1) && receiveSize[partner] > 0) - { - vtkNew receivedPointCloud; - vtkCommunicator::UnMarshalDataObject(dataToReceive[partner], receivedPointCloud.Get()); - - dataToReceive[partner] = nullptr; - vtkIdType nbReceivedPoints = receivedPointCloud->GetNumberOfPoints(); - vtkPointData* receivedPointData = receivedPointCloud->GetPointData(); - vtkPoints* receivedPoints = receivedPointCloud->GetPoints(); - vtkPoints* outputPoints = output->GetPoints(); - if (!outputPoints) - { - vtkNew points; - outputPoints = points.Get(); - output->SetPoints(outputPoints); - } - for (vtkIdType i = 0; i < nbReceivedPoints; i++) - { - outputPoints->InsertNextPoint(receivedPoints->GetPoint(i)); - } - int nbArray = receivedPointData->GetNumberOfArrays(); - - for (int a = 0; a < nbArray; a++) - { - vtkDataArray* fromArray = receivedPointData->GetArray(a); - if (fromArray) - { - vtkDataArray* toArray = outputPointData->GetArray(fromArray->GetName()); - if (!toArray) - { - toArray = fromArray->NewInstance(); - toArray->SetName(fromArray->GetName()); - toArray->SetNumberOfComponents(fromArray->GetNumberOfComponents()); - toArray->SetNumberOfTuples(totalPointsToReceive); - outputPointData->AddArray(toArray); - toArray->Delete(); - } - - for (vtkIdType i = 0; i < nbReceivedPoints; i++) - { - toArray->SetTuple(totalPoints + i, fromArray->GetTuple(i)); - } - } - } - totalPoints += nbReceivedPoints; - --nReceive; - receiveSize[partner] = 0; - } - } - } - - // we wait for sent messages to be received before deleting them - for (int round = 0; round < np - 1; round++) - { - int partner = (rank + round + 1) % np; - if (messagesSize[partner] > 0) - { - sendRequests[partner].Wait(); - } - } -} - } // END namespace vtkMPIUtilities diff --git a/Parallel/MPI/vtkMPIUtilities.h b/Parallel/MPI/vtkMPIUtilities.h index e1d94186b85..6f0e1fcd2f4 100644 --- a/Parallel/MPI/vtkMPIUtilities.h +++ b/Parallel/MPI/vtkMPIUtilities.h @@ -19,7 +19,6 @@ // Forward declarations class vtkMPIController; -class vtkPointSet; namespace vtkMPIUtilities { @@ -41,16 +40,6 @@ void Printf(vtkMPIController* comm, const char* format, ...); VTKPARALLELMPI_EXPORT void SynchronizedPrintf(vtkMPIController* comm, const char* format, ...); -/** - * Get the points that are inside innerBounds and put them in output DataSet. - * Ask other MPI ranks for their corresponding points. - */ -VTKPARALLELMPI_EXPORT -void GetPointsInsideBounds(vtkMPIController* controller, - vtkPointSet* input, - vtkPointSet* output, - const double innerBounds[6]); - } // END namespace vtkMPIUtilities #endif // vtkMPIUtilities_h