Commit c59912e6 authored by Cory Quammen's avatar Cory Quammen
Browse files

Added parallel connectivity filter

This new filter is a parallel version of vtkConnectivityFilter. It
first computes the connectivity on each process, then ties together
the different regions across process boundaries.

* ScalarConnectivity mode is not supported.

* Point and cell seeded connectivity is not supported.

* The extraction modes VTK_EXTRACT_POINT_SEEDED_REGIONS,
  VTK_EXTRACT_CELL_SEEDED_REGIONS, and VTK_EXTRACT_SPECIFIED_REGIONS
    are not supported.

Also added a test for vtkPConnectivityFilter that runs on 1 process
and 4 processes.
parent 08d60354
......@@ -11,6 +11,7 @@ set(Module_SRCS
vtkExtractUnstructuredGridPiece.cxx
vtkExtractUserDefinedPiece.cxx
vtkPCellDataToPointData.cxx
vtkPConnectivityFilter.cxx
vtkPExtractArraysOverTime.cxx
vtkPeriodicFilter.cxx
vtkPKdTree.cxx
......
......@@ -17,6 +17,12 @@ vtk_add_test_mpi(${vtk-module}CxxTests-MPI no_data_tests
ParallelResampling.cxx
UnitTestPMaskPoints.cxx
)
set(${vtk-module}CxxTests-MPI_NUMPROCS 1)
vtk_add_test_mpi(${vtk-module}CxxTests-MPI data_tests_1_proc
ParallelConnectivity1,ParallelConnectivity.cxx,TESTING_DATA,NO_VALID
)
# We want 4 processes to test the vtkAggregateDataSetFilter
# properly
set(${vtk-module}CxxTests-MPI_NUMPROCS 4)
......@@ -24,10 +30,16 @@ vtk_add_test_mpi(${vtk-module}CxxTests-MPI no_data_tests_4_procs
AggregateDataSet.cxx
)
vtk_add_test_mpi(${vtk-module}CxxTests-MPI data_tests_4_procs
ParallelConnectivity4,ParallelConnectivity.cxx,TESTING_DATA,NO_VALID
)
set(all_tests
${tests}
${data_tests_1_proc}
${no_data_tests}
${no_data_tests_4_procs}
${data_tests_4_procs}
)
vtk_test_mpi_executable(${vtk-module}CxxTests-MPI all_tests)
vtk_test_cxx_executable(${vtk-module}CxxTests testsStd)
/*=========================================================================
Program: Visualization Toolkit
Module: ParallelConnectivity.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 "vtkPConnectivityFilter.h"
#include "vtkContourFilter.h"
#include "vtkDataSetTriangleFilter.h"
#include "vtkDistributedDataFilter.h"
#include "vtkMPIController.h"
#include "vtkPUnstructuredGridGhostCellsGenerator.h"
#include "vtkPConnectivityFilter.h"
#include "vtkRemoveGhosts.h"
#include "vtkStructuredPoints.h"
#include "vtkStructuredPointsReader.h"
#include "vtkTestUtilities.h"
#include "vtkUnstructuredGrid.h"
#include <mpi.h>
int ParallelConnectivity(int argc, char* argv[])
{
int returnValue = EXIT_SUCCESS;
MPI_Init(&argc, &argv);
// Note that this will create a vtkMPIController if MPI
// is configured, vtkThreadedController otherwise.
vtkMPIController *contr = vtkMPIController::New();
contr->Initialize(&argc, &argv, 1);
vtkMultiProcessController::SetGlobalController(contr);
int me = contr->GetLocalProcessId();
vtkNew<vtkStructuredPointsReader> reader;
vtkDataSet* ds;
vtkSmartPointer<vtkUnstructuredGrid> ug = vtkSmartPointer<vtkUnstructuredGrid>::New();
if (me == 0)
{
char* fname =
vtkTestUtilities::ExpandDataFileName(argc, argv, "Data/ironProt.vtk");
std::cout << fname << std::endl;
reader->SetFileName(fname);
delete[] fname;
reader->Update();
ds = reader->GetOutput();
}
else
{
ds = ug;
}
vtkNew<vtkDistributedDataFilter> dd;
dd->SetInputData(ds);
dd->SetController(contr);
dd->UseMinimalMemoryOff();
dd->SetBoundaryModeToAssignToOneRegion();
vtkNew<vtkContourFilter> contour;
contour->SetInputConnection(dd->GetOutputPort());
contour->SetNumberOfContours(1);
contour->SetValue(0, 240.0);
vtkNew<vtkDataSetTriangleFilter> tetrahedralize;
tetrahedralize->SetInputConnection(contour->GetOutputPort());
vtkNew<vtkPUnstructuredGridGhostCellsGenerator> ghostCells;
ghostCells->SetController(contr);
ghostCells->SetBuildIfRequired(false);
ghostCells->SetMinimumNumberOfGhostLevels(1);
ghostCells->SetInputConnection(tetrahedralize->GetOutputPort());
vtkNew<vtkPConnectivityFilter> connectivity;
connectivity->SetInputConnection(ghostCells->GetOutputPort());
connectivity->Update();
// Remove ghost points/cells so that the cell count is the same regardless
// of the number of processes.
vtkNew<vtkRemoveGhosts> removeGhosts;
removeGhosts->SetInputConnection(connectivity->GetOutputPort());
// Check the number of regions
int numberOfRegions = connectivity->GetNumberOfExtractedRegions();
int expectedNumberOfRegions = 19;
if (numberOfRegions != expectedNumberOfRegions)
{
std::cerr << "Expected " << expectedNumberOfRegions << " regions but got "
<< numberOfRegions << std::endl;
returnValue = EXIT_FAILURE;
}
// Check the number of cells in the largest region when the extraction mode
// is set to largest region.
connectivity->SetExtractionModeToLargestRegion();
removeGhosts->Update();
int numberOfCells =
vtkPointSet::SafeDownCast(removeGhosts->GetOutput())->GetNumberOfCells();
int globalNumberOfCells = 0;
contr->AllReduce(&numberOfCells, &globalNumberOfCells, 1, vtkCommunicator::SUM_OP);
int expectedNumberOfCells = 2124;
if (globalNumberOfCells != expectedNumberOfCells)
{
std::cerr << "Expected " << expectedNumberOfCells << " cells in largest "
<< "region bug got " << globalNumberOfCells << std::endl;
returnValue = EXIT_FAILURE;
}
// Closest point region test
connectivity->SetExtractionModeToClosestPointRegion();
removeGhosts->Update();
numberOfCells =
vtkPointSet::SafeDownCast(removeGhosts->GetOutput())->GetNumberOfCells();
contr->AllReduce(&numberOfCells, &globalNumberOfCells, 1, vtkCommunicator::SUM_OP);
expectedNumberOfCells = 862; // point (0, 0, 0)
if (globalNumberOfCells != expectedNumberOfCells)
{
std::cerr << "Expected " << expectedNumberOfCells << " cells in closest "
<< "point extraction mode but got " << globalNumberOfCells << std::endl;
returnValue = EXIT_FAILURE;
}
contr->Finalize();
contr->Delete();
return returnValue;
}
......@@ -9,6 +9,7 @@ vtk_module(vtkFiltersParallel
vtkIOXML
vtkRendering${VTK_RENDERING_BACKEND}
vtkRenderingParallel
vtkFiltersParallelGeometry
vtkFiltersParallelMPI
vtkFiltersParallelImaging
vtkFiltersFlowPaths
......
/*=========================================================================
Program: Visualization Toolkit
Module: vtkPConnectivityFilter.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 "vtkPConnectivityFilter.h"
#include "vtkCellData.h"
#include "vtkCellIterator.h"
#include "vtkDataSet.h"
#include "vtkDataSetSurfaceFilter.h"
#include "vtkIdTypeArray.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkKdTreePointLocator.h"
#include "vtkMath.h"
#include "vtkMultiProcessController.h"
#include "vtkNew.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkThreshold.h"
#include "vtkUnsignedCharArray.h"
#include "vtkUnstructuredGrid.h"
#include <algorithm>
#include <map>
#include <numeric>
#include <set>
#include <vector>
vtkStandardNewMacro(vtkPConnectivityFilter);
vtkPConnectivityFilter::vtkPConnectivityFilter()
{
}
vtkPConnectivityFilter::~vtkPConnectivityFilter()
{
}
int vtkPConnectivityFilter::RequestData(
vtkInformation *request,
vtkInformationVector **inputVector,
vtkInformationVector *outputVector)
{
if (this->ExtractionMode == VTK_EXTRACT_POINT_SEEDED_REGIONS ||
this->ExtractionMode == VTK_EXTRACT_CELL_SEEDED_REGIONS ||
this->ExtractionMode == VTK_EXTRACT_SPECIFIED_REGIONS)
{
vtkErrorMacro("ExtractionMode " << this->GetExtractionModeAsString()
<< " is not supported in " << this->GetClassName());
return 1;
}
vtkMultiProcessController* controller =
vtkMultiProcessController::GetGlobalController();
int numRanks = controller->GetNumberOfProcesses();
int myRank = controller->GetLocalProcessId();
// Compute local connectivity. If we are running in parallel, we need the full
// connectivity first, and will handle the extraction mode later.
int success = 1;
if (numRanks > 1)
{
int saveScalarConnectivity = this->ScalarConnectivity;
int saveExtractionMode = this->ExtractionMode;
int saveColorRegions = this->ColorRegions;
// Overwrite custom member variables temporarily.
this->ScalarConnectivity = 0;
this->ExtractionMode = VTK_EXTRACT_ALL_REGIONS;
this->ColorRegions = 1;
// Invoke the connectivity algorithm in the superclass.
success = this->Superclass::RequestData(request, inputVector, outputVector);
this->ScalarConnectivity = saveScalarConnectivity;
this->ExtractionMode = saveExtractionMode;
this->ColorRegions = saveColorRegions;
}
else
{
// Only 1 process, just invoke the superclass and return.
return this->Superclass::RequestData(request, inputVector, outputVector);
}
// Get the info objects
vtkInformation *outInfo = outputVector->GetInformationObject(0);
// Get the input and output
vtkPointSet *output = vtkPointSet::SafeDownCast(
outInfo->Get(vtkDataObject::DATA_OBJECT()));
// Check that we have at least one ghost level. If we don't have any ghost
// points, we say that we have ghost points so other ranks can proceed.
bool hasGhostPoints = output->GetNumberOfPoints() == 0 ||
(output->GetNumberOfPoints() > 0 && output->GetPointGhostArray());
success *= hasGhostPoints ? 1 : 0;
// Check that all ranks succeeded in local connectivity connection and
// have ghost cells.
int globalSuccess = 0;
controller->AllReduce(&success, &globalSuccess, 1, vtkCommunicator::MIN_OP);
if (!globalSuccess)
{
if (!hasGhostPoints)
{
vtkErrorMacro("At least one ghost level is required to run this filter in "
"parallel, but no ghost cells are available.");
}
return 0;
}
// Exchange number of regions. We assume the RegionIDs are contiguous.
int numRegions = this->GetNumberOfExtractedRegions();
std::vector<int> regionCounts(numRanks, 0);
std::vector<int> regionStarts(numRanks + 1, 0);
controller->AllGather(&numRegions, &regionCounts[0], 1);
// Compute starting region Ids on each rank
std::partial_sum(regionCounts.begin(), regionCounts.end(),
regionStarts.begin() + 1);
vtkPoints* outputPoints = output->GetPoints();
vtkPointData* outputPD = output->GetPointData();
vtkCellData* outputCD = output->GetCellData();
vtkIdTypeArray* cellRegionIds =
vtkIdTypeArray::SafeDownCast(outputCD->GetArray("RegionId"));
vtkIdTypeArray* pointRegionIds =
vtkIdTypeArray::SafeDownCast(outputPD->GetArray("RegionId"));
// Gather up ghost points in the dataset and send them to all the other ranks.
vtkUnsignedCharArray* pointGhostArray = output->GetPointGhostArray();
// Extract ghost points and their regions.
vtkSmartPointer<vtkDataArray> ghostPoints;
ghostPoints.TakeReference(outputPoints->GetData()->NewInstance());
ghostPoints->SetNumberOfComponents(outputPoints->GetData()->GetNumberOfComponents());
vtkSmartPointer<vtkIdTypeArray> ghostRegionIds =
vtkSmartPointer<vtkIdTypeArray>::New();
for (vtkIdType i = 0; i < output->GetNumberOfPoints(); ++i)
{
if (pointGhostArray->GetTypedComponent(i, 0) &
vtkDataSetAttributes::DUPLICATEPOINT)
{
ghostPoints->InsertNextTuple(i, outputPoints->GetData());
// Add the region id starting value for this rank so it is in terms of the
// global region numbering.
vtkIdType regionId = pointRegionIds->GetTypedComponent(i, 0);
ghostRegionIds->InsertNextValue(regionId + regionStarts[myRank]);
}
}
// Gather up number of point tuples on each rank
std::vector<vtkIdType> remotePointLengths(numRanks, 0);
vtkIdType localPointsLength = ghostPoints->GetNumberOfValues();
controller->AllGather(&localPointsLength, &remotePointLengths[0], 1);
std::vector<vtkIdType> remotePointOffsets(numRanks+1, 0);
std::partial_sum(remotePointLengths.begin(), remotePointLengths.end(),
remotePointOffsets.begin() + 1);
// Gather points
vtkSmartPointer<vtkDataArray> remotePointData;
remotePointData.TakeReference(ghostPoints->NewInstance());
remotePointData->SetNumberOfComponents(3);
remotePointData->SetNumberOfTuples(remotePointOffsets[numRanks] / 3);
controller->AllGatherV(ghostPoints, remotePointData,
&remotePointLengths[0], &remotePointOffsets[0]);
// The remote point lengths and offsets account for the point data being
// stored as 3-tuples. Divide by 3 in preparation for distributing 1-tuple
// point region ids.
for (int i = 0; i < numRanks; ++i)
{
remotePointLengths[i] /= 3;
remotePointOffsets[i+1] /= 3;
}
// Gather the RegionIds associated with the ghost points.
vtkSmartPointer<vtkIdTypeArray> remoteRegionIds =
vtkSmartPointer<vtkIdTypeArray>::New();
remoteRegionIds->SetNumberOfComponents(1);
remoteRegionIds->SetNumberOfTuples(remotePointOffsets[numRanks]);
remoteRegionIds->FillValue(-1); // Invalid region id
controller->AllGatherV(ghostRegionIds, remoteRegionIds,
&remotePointLengths[0], &remotePointOffsets[0]);
ghostRegionIds = nullptr;
// Links from local region ids to remote region ids. Vector index is local
// region id, and the set contains linked remote ids.
typedef std::vector< std::set< vtkIdType > > RegionLinksType;
RegionLinksType links(regionStarts[numRanks]);
if (output->GetNumberOfPoints() > 0)
{
// Now resolve the remote ghost points to local points if possible.
// TODO - use a surface filter on the output to reduce the number of points
// in the point locator.
vtkNew<vtkKdTreePointLocator> localPointLocator;
localPointLocator->SetDataSet(output);
localPointLocator->BuildLocator();
// Map the local and remote ids
for (int rank = 0; rank < numRanks; ++rank)
{
if (rank == myRank)
{
continue;
}
for (vtkIdType remoteIdx = remotePointOffsets[rank];
remoteIdx < remotePointOffsets[rank+1]; ++remoteIdx)
{
double x[3];
remotePointData->GetTuple(remoteIdx, x);
vtkIdType localId = localPointLocator->FindClosestPoint(x);
// Skip local ghost points as we do not need ghost-ghost links.
if (pointGhostArray->GetTypedComponent(localId, 0) &
vtkDataSetAttributes::DUPLICATEPOINT)
{
continue;
}
double y[3];
output->GetPoints()->GetPoint(localId, y);
double dist2 = vtkMath::Distance2BetweenPoints(x, y);
if (dist2 > 1e-6)
{
// Nearest point is too far away, move on.
continue;
}
// Save association between local and remote ids
vtkIdTypeArray* localRegionIds =
vtkIdTypeArray::SafeDownCast(outputPD->GetArray("RegionId"));
vtkIdType localRegionId = localRegionIds->GetTypedComponent(localId, 0) +
regionStarts[myRank];
vtkIdType remoteRegionId = remoteRegionIds->GetTypedComponent(remoteIdx, 0);
if (links[localRegionId].count(remoteRegionId) == 0)
{
// Link not already established. Add it here.
links[localRegionId].insert(remoteRegionId);
}
}
}
}
remoteRegionIds = nullptr;
// Set up storage for gathering all links from all processors. This is an
// interleaved vector containing one regionId and its connected regionId.
std::vector< vtkIdType > localLinks;
for (size_t i = 0; i < links.size(); ++i)
{
for (std::set< vtkIdType >::iterator iter = std::begin(links[i]);
iter != std::end(links[i]); ++iter)
{
localLinks.push_back(static_cast<vtkIdType>(i));
localLinks.push_back(*iter);
}
}
// Gather all the links on each rank. This is possibly suboptimal, but it
// avoids needing a connected components algorithm on a distributed graph.
vtkIdType localNumLinks = static_cast<vtkIdType>(localLinks.size());
std::vector< vtkIdType > linkCounts(numRanks, -1);
std::vector< vtkIdType > linkStarts(numRanks + 1, 0);
controller->AllGather(&localNumLinks, &linkCounts[0], 1);
// Compute starting region IDs on each rank
for (int i = 0; i < numRanks; ++i)
{
linkStarts[i + 1] = linkCounts[i] + linkStarts[i];
}
std::vector< vtkIdType > allLinks(linkStarts[numRanks]);
controller->AllGatherV(&localLinks[0], &allLinks[0],
static_cast<vtkIdType>(localLinks.size()),
&linkCounts[0], &linkStarts[0]);
// Set up a graph of all the region-to-region links.
typedef struct _RegionNode {
// Stored for relabeling step
vtkIdType OriginalRegionId;
// Current local region id
vtkIdType CurrentRegionId;
std::vector< vtkIdType > Links;
} RegionNode;
size_t linkIdx = 0;
std::vector< RegionNode > regionNodes(regionStarts[numRanks]);
for (vtkIdType regionId = 0; regionId < regionStarts[numRanks]; ++regionId)
{
regionNodes[regionId].OriginalRegionId = regionId;
regionNodes[regionId].CurrentRegionId = regionId;
while (linkIdx < allLinks.size() && allLinks[linkIdx] == regionId)
{
regionNodes[regionId].Links.push_back(allLinks[linkIdx + 1]);
linkIdx += 2;
}
}
// Now run connected components on this graph. The algorithm labels all
// connected nodes in the graph with the lowest region id in the connected
// component. This is a breadth-first algorithm. I'm not 100% sure that the
// multiple passes in the do-while loop are required, but I suspect there may
// be graph configurations where a single pass is not sufficient for the
// relabeling to converge.
bool componentChanged = false;
do
{
componentChanged = false;
for (std::vector< RegionNode >::iterator nodeIter = regionNodes.begin();
nodeIter != regionNodes.end(); ++nodeIter)
{
for (std::vector< vtkIdType >::iterator linkIter = nodeIter->Links.begin();
linkIter != nodeIter->Links.end(); ++linkIter)
{
vtkIdType linkedRegionId = *linkIter;
if (nodeIter->CurrentRegionId < regionNodes[linkedRegionId].CurrentRegionId)
{
regionNodes[linkedRegionId].CurrentRegionId = nodeIter->CurrentRegionId;
componentChanged = true;
}
}
}
} while (componentChanged);
// Collect all the current ids remaining after the connected components
// algorithm.
std::set< vtkIdType > currentRegionIds;
for (std::vector< RegionNode >::iterator nodeIter = regionNodes.begin();
nodeIter != regionNodes.end(); ++nodeIter)
{
currentRegionIds.insert(nodeIter->CurrentRegionId);
}
// Create a map from current region id after relabeling to a new, contiguous
// label. Maps current region id -> relabeled array.
std::map< vtkIdType, vtkIdType > relabeledRegionMap;
vtkIdType contiguousLabel = 0;
for (std::set< vtkIdType >::iterator setIter = currentRegionIds.begin();
setIter != currentRegionIds.end(); ++setIter)
{
relabeledRegionMap[*setIter] = contiguousLabel++;
}
// Now do the relabing to the contiguous region id.
std::vector< vtkIdType > regionIdMap(regionNodes.size(), -1);
for (std::vector< RegionNode >::iterator nodeIter = regionNodes.begin();
nodeIter != regionNodes.end(); ++nodeIter)
{
nodeIter->CurrentRegionId = relabeledRegionMap[nodeIter->CurrentRegionId];
}
// Relabel the points and cells according to the contiguous renumbering.
for (vtkIdType i = 0; i < output->GetNumberOfCells(); ++i)
{
// Offset the cellRegionId by the starting region id on this rank.
vtkIdType cellRegionId = cellRegionIds->GetValue(i) + regionStarts[myRank];
cellRegionIds->SetValue(i, regionNodes[cellRegionId].CurrentRegionId);
}
for (vtkIdType i = 0; i < output->GetNumberOfPoints(); ++i)
{
// Offset the pointRegionId by the starting region id on this rank.
vtkIdType pointRegionId = pointRegionIds->GetValue(i) + regionStarts[myRank];
pointRegionIds->SetValue(i, regionNodes[pointRegionId].CurrentRegionId);
}
// Sum up number of cells in each region.
vtkIdType numContiguousLabels = contiguousLabel;
std::vector< vtkIdType > localRegionSizes(numContiguousLabels, 0);
if (cellRegionIds)
{
// Iterate over cells and count how many are in different regions.
for (vtkIdType i = 0; i < cellRegionIds->GetNumberOfValues(); ++i)
{
localRegionSizes[cellRegionIds->GetValue(i)]++;
}
}
// AllReduce to sum up the number of cells in each region on each process.
std::vector< vtkIdType > globalRegionSizes(numContiguousLabels, 0);
controller->AllReduce(&localRegionSizes[0], &globalRegionSizes[0],
numContiguousLabels, vtkCommunicator::SUM_OP);
// Store the region sizes
this->RegionSizes->Reset();