Commit ef1611f7 authored by George Zagaris's avatar George Zagaris

ENH: RectilinearGrid implicit connectivity support

Added support for vtkRectilinearGrid type in the structured
implicity connectivity class. Also, extended the test to
check for correctness of the new functionality. Further,
updated RectilinearGrid partitioner such that it can either leave a
gap between partitions or not. Also, refactored and improved the
implementation.

Change-Id: I658e2798bd9720a7baafd7eb1eed3f1e8befafe0
parent 8a4e1cf1
......@@ -13,16 +13,19 @@
=========================================================================*/
#include "vtkRectilinearGridPartitioner.h"
#include "vtkObjectFactory.h"
#include "vtkIndent.h"
// VTK includes
#include "vtkDoubleArray.h"
#include "vtkExtentRCBPartitioner.h"
#include "vtkRectilinearGrid.h"
#include "vtkStructuredData.h"
#include "vtkStructuredExtent.h"
#include "vtkIndent.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkMultiBlockDataSet.h"
#include "vtkDoubleArray.h"
#include "vtkObjectFactory.h"
#include "vtkRectilinearGrid.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkStructuredData.h"
#include "vtkStructuredExtent.h"
#include <cassert>
......@@ -33,6 +36,7 @@ vtkRectilinearGridPartitioner::vtkRectilinearGridPartitioner()
{
this->NumberOfPartitions = 2;
this->NumberOfGhostLayers = 0;
this->DuplicateNodesOn();
this->SetNumberOfInputPorts(1);
this->SetNumberOfOutputPorts(1);
}
......@@ -80,169 +84,31 @@ void vtkRectilinearGridPartitioner::ExtractGridCoordinates(
int dataDescription = vtkStructuredData::GetDataDescriptionFromExtent(subext);
int dims[3];
vtkStructuredData::GetDimensionsFromExtent(subext,dims,dataDescription );
int ndims[3];
vtkStructuredData::GetDimensionsFromExtent(subext,ndims,dataDescription );
int i,j,k;
int lidx; // local linear index
double p[3];
switch( dataDescription )
vtkDoubleArray* coords[3];
coords[0] = xcoords;
coords[1] = ycoords;
coords[2] = zcoords;
vtkDataArray* src_coords[3];
src_coords[0] = grd->GetXCoordinates();
src_coords[1] = grd->GetYCoordinates();
src_coords[2] = grd->GetZCoordinates();
for(int dim=0; dim < 3; ++dim)
{
case VTK_XY_PLANE:
// Sanity checks
assert("pre: dimension should be greater than 0" && (dims[0] > 0) );
assert("pre: dimension should be greater than 0" && (dims[1] > 0) );
// Allocate coordinate arrays
xcoords->SetNumberOfComponents(1);
xcoords->SetNumberOfTuples(dims[0]);
ycoords->SetNumberOfComponents(1);
ycoords->SetNumberOfTuples(dims[1]);
zcoords->SetNumberOfComponents(1);
zcoords->SetNumberOfTuples(0);
// Copy coordinates from the whole grid
i=subext[0]; k=subext[4]; j=subext[2];
for( ; i <= subext[1]; ++i )
{
lidx = i-subext[0];
assert("pre: local linear index is out-of-bounds!" &&
(lidx >= 0) && (lidx < dims[0]) );
grd->GetPoint( i,j,k,p );
xcoords->SetComponent(lidx,0, p[0] );
} // END for all i
i=subext[0]; k=subext[4]; j=subext[2];
for( ; j <= subext[3]; ++j )
{
lidx = j-subext[2];
assert("pre: local linear index is out-of-bounds!" &&
(lidx >= 0) && (lidx < dims[1]) );
grd->GetPoint( i,j,k,p );
ycoords->SetComponent(lidx,0,p[1]);
} // END for all j
break;
case VTK_XZ_PLANE:
// Sanity checks
assert("pre: dimension should be greater than 0" && (dims[0] > 0) );
assert("pre: dimension should be greater than 0" && (dims[2] > 0) );
// Allocate coordinate arrays
xcoords->SetNumberOfComponents(1);
xcoords->SetNumberOfTuples(dims[0]);
zcoords->SetNumberOfComponents(1);
zcoords->SetNumberOfTuples(dims[2]);
ycoords->SetNumberOfComponents(1);
ycoords->SetNumberOfTuples(0);
// Copy coordinates from the whole grid
i=subext[0]; k=subext[4]; j=subext[2];
for( ; i <= subext[1]; ++i )
{
lidx = i-subext[0];
assert("pre: local linear index is out-of-bounds!" &&
(lidx >= 0) && (lidx < dims[0]) );
grd->GetPoint( i,j,k,p );
xcoords->SetComponent(lidx,0, p[0] );
} // END for all i
i=subext[0]; k=subext[4]; j=subext[2];
for( ; k <= subext[5]; ++k )
{
lidx = k-subext[4];
assert("pre: local linear index is out-of-bounds!" &&
(lidx >= 0) && (lidx < dims[2]) );
grd->GetPoint( i,j,k,p );
zcoords->SetComponent( lidx,0,p[2] );
} // END for all k
break;
case VTK_YZ_PLANE:
// Sanity checks
assert("pre: dimension should be greater than 0" && (dims[1] > 0) );
assert("pre: dimension should be greater than 0" && (dims[2] > 0) );
// Allocate coordinate arrays
ycoords->SetNumberOfComponents(1);
ycoords->SetNumberOfTuples(dims[1]);
zcoords->SetNumberOfComponents(1);
zcoords->SetNumberOfTuples(dims[2]);
xcoords->SetNumberOfComponents(1);
xcoords->SetNumberOfTuples(dims[0]);
// Copy coordinates from the whole grid
i=subext[0]; k=subext[4]; j=subext[2];
for( ; j <= subext[3]; ++j )
{
lidx = j-subext[2];
assert("pre: local linear index is out-of-bounds!" &&
(lidx >= 0) && (lidx < dims[1]) );
grd->GetPoint( i,j,k,p );
ycoords->SetComponent(lidx,0,p[1]);
} // END for all j
i=subext[0]; k=subext[4]; j=subext[2];
for( ; k <= subext[5]; ++k )
{
lidx = k-subext[4];
assert("pre: local linear index is out-of-bounds!" &&
(lidx >= 0) && (lidx < dims[2]) );
grd->GetPoint( i,j,k,p );
zcoords->SetComponent( lidx,0,p[2] );
} // END for all k
break;
case VTK_XYZ_GRID:
// Sanity checks
assert("pre: dimension should be greater than 0" && (dims[0] > 0) );
assert("pre: dimension should be greater than 0" && (dims[1] > 0) );
assert("pre: dimension should be greater than 0" && (dims[2] > 0) );
// Allocate coordinate arrays
xcoords->SetNumberOfComponents(1);
xcoords->SetNumberOfTuples( dims[0] );
ycoords->SetNumberOfComponents(1);
ycoords->SetNumberOfTuples( dims[1] );
zcoords->SetNumberOfComponents(1);
zcoords->SetNumberOfTuples( dims[2] );
// Copy coordinates from the whole grid
i=subext[0]; k=subext[4]; j=subext[2];
for( ; i <= subext[1]; ++i )
{
lidx = i-subext[0];
assert("pre: local linear index is out-of-bounds!" &&
(lidx >= 0) && (lidx < dims[0]) );
grd->GetPoint( i,j,k,p );
xcoords->SetComponent(lidx,0, p[0] );
} // END for all i
i=subext[0]; k=subext[4]; j=subext[2];
for( ; j <= subext[3]; ++j )
{
lidx = j-subext[2];
assert("pre: local linear index is out-of-bounds!" &&
(lidx >= 0) && (lidx < dims[1]) );
grd->GetPoint( i,j,k,p );
ycoords->SetComponent(lidx,0,p[1]);
} // END for all j
i=subext[0]; k=subext[4]; j=subext[2];
for( ; k <= subext[5]; ++k )
{
lidx = k-subext[4];
assert("pre: local linear index is out-of-bounds!" &&
(lidx >= 0) && (lidx < dims[2]) );
grd->GetPoint( i,j,k,p );
zcoords->SetComponent( lidx,0,p[2] );
} // END for all k
break;
default:
vtkErrorMacro("Cannot handle structured data!");
}
coords[ dim ]->SetNumberOfComponents( 1 );
coords[ dim ]->SetNumberOfTuples( ndims[ dim ] );
for(int idx=subext[dim*2]; idx <= subext[dim*2+1]; ++idx)
{
vtkIdType lidx = idx-subext[dim*2];
coords[ dim ]->SetTuple1(lidx,src_coords[dim]->GetTuple1(idx));
} // END for all ids
} // END for all dimensions
}
//------------------------------------------------------------------------------
......@@ -274,6 +140,14 @@ int vtkRectilinearGridPartitioner::RequestData(
extentPartitioner->SetGlobalExtent( extent );
extentPartitioner->SetNumberOfPartitions( this->NumberOfPartitions );
extentPartitioner->SetNumberOfGhostLayers( this->NumberOfGhostLayers );
if( this->DuplicateNodes == 1 )
{
extentPartitioner->DuplicateNodesOn();
}
else
{
extentPartitioner->DuplicateNodesOff();
}
// STEP 4: Partition
extentPartitioner->Partition();
......@@ -281,6 +155,10 @@ int vtkRectilinearGridPartitioner::RequestData(
// STEP 5: Extract partition in a multi-block
multiblock->SetNumberOfBlocks( extentPartitioner->GetNumExtents( ) );
// Set the whole extent of the grid
multiblock->GetInformation()->Set(
vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),extent,6);
int subext[6];
unsigned int blockIdx = 0;
for( ; blockIdx < multiblock->GetNumberOfBlocks(); ++blockIdx )
......@@ -293,12 +171,21 @@ int vtkRectilinearGridPartitioner::RequestData(
vtkDoubleArray *ycoords = vtkDoubleArray::New();
vtkDoubleArray *zcoords = vtkDoubleArray::New();
this->ExtractGridCoordinates(grd,subext,xcoords,ycoords,zcoords);
subgrid->SetXCoordinates( xcoords );
subgrid->SetYCoordinates( ycoords );
subgrid->SetZCoordinates( zcoords );
xcoords->Delete();
ycoords->Delete();
zcoords->Delete();
vtkInformation *metadata = multiblock->GetMetaData( blockIdx );
assert( "pre: metadata is NULL" && (metadata != NULL) );
metadata->Set( vtkDataObject::PIECE_EXTENT(), subext, 6);
multiblock->SetBlock(blockIdx, subgrid);
subgrid->Delete();
} // END for all blocks
extentPartitioner->Delete();
......
......@@ -52,6 +52,11 @@ public:
vtkGetMacro(NumberOfGhostLayers,int);
vtkSetMacro(NumberOfGhostLayers,int);
// Description:
vtkGetMacro(DuplicateNodes,int);
vtkSetMacro(DuplicateNodes,int);
vtkBooleanMacro(DuplicateNodes,int);
protected:
vtkRectilinearGridPartitioner();
virtual ~vtkRectilinearGridPartitioner();
......@@ -72,6 +77,7 @@ protected:
int NumberOfPartitions;
int NumberOfGhostLayers;
int DuplicateNodes;
private:
vtkRectilinearGridPartitioner(const vtkRectilinearGridPartitioner &); // Not implemented
......
......@@ -39,6 +39,9 @@
#include "vtkMathUtilities.h"
#include "vtkMultiBlockDataSet.h"
#include "vtkPointData.h"
#include "vtkRectilinearGrid.h"
#include "vtkRectilinearGridPartitioner.h"
#include "vtkRectilinearGridWriter.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkStructuredGrid.h"
#include "vtkStructuredGridPartitioner.h"
......@@ -48,7 +51,8 @@
#include "vtkXMLPMultiBlockDataWriter.h"
//#define DEBUG_ON
#define DEBUG_ON
//------------------------------------------------------------------------------
// G L O B A L D A T A
//------------------------------------------------------------------------------
......@@ -88,8 +92,8 @@ void AddNodeCenteredXYZField( vtkMultiBlockDataSet *mbds )
for( unsigned int block=0; block < mbds->GetNumberOfBlocks(); ++block )
{
vtkStructuredGrid *grid =
vtkStructuredGrid::SafeDownCast(mbds->GetBlock(block));
vtkDataSet* grid = vtkDataSet::SafeDownCast(mbds->GetBlock(block));
if( grid == NULL )
{
continue;
......@@ -212,6 +216,175 @@ vtkMultiBlockDataSet* GetDataSet(
return( mbds );
}
//------------------------------------------------------------------------------
double exponential_distribution(const int i, const double beta)
{
double xi=0.0;
xi = ( ( exp( i*beta ) - 1 ) /( exp( beta ) - 1 ) );
return( xi );
}
//------------------------------------------------------------------------------
void GenerateRectGrid( vtkRectilinearGrid* grid, int ext[6], double origin[3])
{
grid->Initialize();
grid->SetExtent(ext);
vtkDataArray* coords[3];
int dims[3];
int dataDesc = vtkStructuredData::GetDataDescriptionFromExtent(ext);
vtkStructuredData::GetDimensionsFromExtent(ext,dims,dataDesc);
// compute & populate coordinate vectors
double beta = 0.01; /* controls the intensity of the stretching */
for(int i=0; i < 3; ++i)
{
coords[i] = vtkDataArray::CreateDataArray(VTK_DOUBLE);
if( dims[i] == 0 )
{
continue;
}
coords[i]->SetNumberOfTuples(dims[i]);
double prev = origin[i];
for(int j=0; j < dims[i]; ++j)
{
double val = prev + ( (j==0)? 0.0 : exponential_distribution(j,beta) );
coords[ i ]->SetTuple( j, &val );
prev = val;
} // END for all points along this dimension
} // END for all dimensions
grid->SetXCoordinates( coords[0] );
grid->SetYCoordinates( coords[1] );
grid->SetZCoordinates( coords[2] );
coords[0]->Delete();
coords[1]->Delete();
coords[2]->Delete();
}
//------------------------------------------------------------------------------
// Description:
// Generates a distributed multi-block dataset, each grid is added using
// round-robin assignment.
vtkMultiBlockDataSet* GetRectGridDataSet(
const int numPartitions,
double origin[3],
int wholeExtent[6])
{
int dims[3];
int desc = vtkStructuredData::GetDataDescriptionFromExtent(wholeExtent);
vtkStructuredData::GetDimensionsFromExtent(wholeExtent,dims,desc);
vtkRectilinearGrid* wholeGrid = vtkRectilinearGrid::New();
GenerateRectGrid(wholeGrid,wholeExtent,origin);
//#ifdef DEBUG_ON
// if( Controller->GetLocalProcessId() == 0 )
// {
// vtkRectilinearGridWriter* writer = vtkRectilinearGridWriter::New();
// writer->SetFileName("RectilinearGrid.vtk");
// writer->SetInputData( wholeGrid );
// writer->Write();
// writer->Delete();
// }
//#endif
vtkRectilinearGridPartitioner *gridPartitioner =
vtkRectilinearGridPartitioner::New();
gridPartitioner->SetInputData( wholeGrid );
gridPartitioner->SetNumberOfPartitions( numPartitions );
gridPartitioner->SetNumberOfGhostLayers(0);
gridPartitioner->DuplicateNodesOff(); /* to create a gap */
gridPartitioner->Update();
vtkMultiBlockDataSet *partitionedGrid =
vtkMultiBlockDataSet::SafeDownCast( gridPartitioner->GetOutput() );
assert( "pre: partitionedGrid != NULL" && (partitionedGrid != NULL) );
// Each process has the same number of blocks, i.e., the same structure,
// however some block entries are NULL indicating that the data lives on
// some other process
vtkMultiBlockDataSet *mbds = vtkMultiBlockDataSet::New();
mbds->SetNumberOfBlocks( numPartitions );
mbds->GetInformation()->Set(
vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),
partitionedGrid->GetInformation()->Get(
vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT()),
6 );
// Populate blocks for this process
unsigned int block=0;
for( ; block < partitionedGrid->GetNumberOfBlocks(); ++block )
{
if( Rank == static_cast<int>( block%NumberOfProcessors ) )
{
// Copy the structured grid
vtkRectilinearGrid *grid = vtkRectilinearGrid::New();
grid->DeepCopy( partitionedGrid->GetBlock(block) );
mbds->SetBlock( block, grid );
grid->Delete();
// Copy the global extent for the blockinformation
vtkInformation *info = partitionedGrid->GetMetaData( block );
assert( "pre: null metadata!" && (info != NULL) );
assert( "pre: must have a piece extent!" &&
(info->Has(vtkDataObject::PIECE_EXTENT() ) ) );
vtkInformation *metadata = mbds->GetMetaData( block );
assert( "pre: null metadata!" && (metadata != NULL) );
metadata->Set(
vtkDataObject::PIECE_EXTENT(),
info->Get( vtkDataObject::PIECE_EXTENT() ),
6 );
} // END if we own the block
else
{
mbds->SetBlock( block, NULL );
} // END else we don't own the block
} // END for all blocks
wholeGrid->Delete();
gridPartitioner->Delete();
assert( "pre: mbds is NULL" && (mbds != NULL) );
AddNodeCenteredXYZField( mbds );
Controller->Barrier();
return mbds;
}
//------------------------------------------------------------------------------
void RegisterRectGrid(vtkMultiBlockDataSet* mbds,
vtkStructuredImplicitConnectivity* connectivity)
{
assert( "pre: Multi-block is NULL!" && (mbds != NULL) );
assert( "pre: connectivity is NULL!" && (connectivity != NULL) );
for(unsigned int block=0; block < mbds->GetNumberOfBlocks(); ++block)
{
vtkRectilinearGrid* grid =
vtkRectilinearGrid::SafeDownCast(mbds->GetBlock(block));
if( grid != NULL )
{
vtkInformation* info = mbds->GetMetaData( block );
assert( "pre: metadata should not be NULL" && (info != NULL) );
assert( "pre: must have piece extent!" &&
info->Has(vtkDataObject::PIECE_EXTENT() ) );
connectivity->RegisterRectilinearGrid(
block,
info->Get(vtkDataObject::PIECE_EXTENT()),
grid->GetXCoordinates(),
grid->GetYCoordinates(),
grid->GetZCoordinates(),
grid->GetPointData()
);
} // END if block belongs to this process
} // END for all blocks
}
//------------------------------------------------------------------------------
void RegisterGrid(
vtkMultiBlockDataSet *mbds,
......@@ -242,7 +415,7 @@ void RegisterGrid(
}
//------------------------------------------------------------------------------
int CheckGrid(vtkStructuredGrid* grid)
int CheckGrid(vtkDataSet* grid)
{
assert("pre: grid should not be NULL!" && (grid != NULL) );
int rc = 0;
......@@ -299,18 +472,36 @@ int TestOutput(vtkMultiBlockDataSet* mbds, int wholeExtent[6])
vtkStructuredImplicitConnectivity::New();
gridConnectivity->SetWholeExtent(wholeExtent);
vtkStructuredGrid* grid = NULL;
vtkDataSet* grid = NULL;
for(unsigned int block=0; block < mbds->GetNumberOfBlocks(); ++block)
{
grid = vtkStructuredGrid::SafeDownCast(mbds->GetBlock(block));
grid = vtkDataSet::SafeDownCast( mbds->GetBlock(block) );
if( grid != NULL )
{
gridConnectivity->RegisterGrid(
block,
grid->GetExtent(),
grid->GetPoints(),
grid->GetPointData()
);
if( grid->IsA("vtkStructuredGrid") )
{
vtkStructuredGrid* sGrid = vtkStructuredGrid::SafeDownCast(grid);
gridConnectivity->RegisterGrid(
block,
sGrid->GetExtent(),
sGrid->GetPoints(),
sGrid->GetPointData()
);
}
else
{
assert("pre: expected rectilinear grid!" &&
grid->IsA("vtkRectilinearGrid"));
vtkRectilinearGrid* rGrid = vtkRectilinearGrid::SafeDownCast(grid);
gridConnectivity->RegisterRectilinearGrid(
block,
rGrid->GetExtent(),
rGrid->GetXCoordinates(),
rGrid->GetYCoordinates(),
rGrid->GetZCoordinates(),
rGrid->GetPointData()
);
}
rc += CheckGrid(grid);
} // END if grid != NULL
......@@ -697,6 +888,92 @@ int TestImplicitGridConnectivity3D( )
return( rc );
}
//------------------------------------------------------------------------------
int TestRectGridImplicitConnectivity3D()
{
assert( "pre: MPI Controller is NULL!" && (Controller != NULL) );
vtkMPIUtilities::Printf(vtkMPIController::SafeDownCast(Controller),
"=======================\nTesting 3-D Rectilinear Grid Dataset\n");
int rc = 0;
int wholeExtent[6];
wholeExtent[0] = 0;
wholeExtent[1] = 99;
wholeExtent[2] = 0;
wholeExtent[3] = 99;
wholeExtent[4] = 0;
wholeExtent[5] = 99;
double origin[3] = {0.0,0.0,0.0};
// STEP 0: We generate the same number of partitions as processes
int numPartitions = NumberOfProcessors;
// STEP 1: Acquire the distributed rectilinear grid for this process.
// Each process has the same number of blocks, but not all entries are
// populated. A NULL entry indicates that the block belongs to a different
// process.
vtkMultiBlockDataSet *mbds =
GetRectGridDataSet( numPartitions,origin,wholeExtent );
Controller->Barrier();
assert( "pre: mbds != NULL" && (mbds != NULL) );
assert( "pre: numBlocks mismatch" &&
(static_cast<int>(mbds->GetNumberOfBlocks())==numPartitions) );
WriteDistributedDataSet("INPUT-3D-RECTGRID",mbds);
// STEP 2: Setup the grid connecti