Commit 0f33c6a7 authored by George Zagaris's avatar George Zagaris Committed by Code Review
Browse files

Merge topic 'unstructured-ghost-zones' into master

a07749d1 ENH: Ghost zone generation for unstructured grids
7c35e23a ENH: Add vtkMPIUtilities -- parallel printf utils
0f73b924 ENH: NearlyEqual fuzzy comparator and SafeDivision
55eb011a ENH: Support for VtkIdType & enhancements
01204dd4 ENH: multi-process stream enhancements
parents 09a16733 a07749d1
......@@ -46,6 +46,45 @@ bool FuzzyCompare(A a, A b, A epsilon)
return fabs(a - b) < epsilon;
}
// Description:
// Performs safe division that catches overflow and underflow.
template<class A>
A SafeDivision(A a, A b)
{
// Avoid overflow
if( (b < static_cast<A>(1)) && (a > b*std::numeric_limits<A>::max()) )
{
return std::numeric_limits<A>::max();
}
// Avoid underflow
if( (a == static_cast<A>(0)) ||
((b > static_cast<A>(1)) && (a < b*std::numeric_limits<A>::min())) )
{
return static_cast<A>(0);
}
// safe to do the division
return( a/b );
}
// Description:
// A slightly different fuzzy comparator that checks if two values are
// "nearly" equal based on Knuth, "The Art of Computer Programming (vol II)"
template<class A>
bool NearlyEqual(A a, A b, A tol=std::numeric_limits<A>::epsilon())
{
A absdiff = fabs(a-b);
A d1 = vtkMathUtilities::SafeDivision<A>(absdiff,fabs(a));
A d2 = vtkMathUtilities::SafeDivision<A>(absdiff,fabs(b));
if( (d1 <= tol) || (d2 <= tol) )
{
return true;
}
return false;
}
} // End vtkMathUtilities namespace.
#endif // __vtkMathUtilities_h
......
......@@ -3,6 +3,8 @@ set(Module_SRCS
vtkPStructuredGridConnectivity.cxx
vtkPStructuredGridGhostDataGenerator.cxx
vtkPUniformGridGhostDataGenerator.cxx
vtkPUnstructuredGridConnectivity.cxx
vtkPUnstructuredGridGhostDataGenerator.cxx
)
set_source_files_properties(
......
include(vtkMPI)
vtk_add_test_mpi(TestPStructuredGridConnectivity.cxx)
vtk_mpi_link(TestPStructuredGridConnectivity)
vtk_add_test_mpi(TestPStructuredGridGhostDataGenerator.cxx)
vtk_add_test_mpi(TestPUniformGridGhostDataGenerator.cxx)
vtk_mpi_link(TestPStructuredGridConnectivity)
vtk_add_test_mpi(TestPUnstructuredGridConnectivity.cxx)
vtk_mpi_link(TestPUnstructuredGridConnectivity)
vtk_add_test_mpi(TestPUnstructuredGridGhostDataGenerator.cxx)
vtk_mpi_link(TestPUnstructuredGridGhostDataGenerator)
\ No newline at end of file
/*=========================================================================
Program: Visualization Toolkit
Module: TestPStructuredGridConnectivity.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 "vtkCellData.h"
#include "vtkCellType.h"
#include "vtkDoubleArray.h"
#include "vtkExtentRCBPartitioner.h"
#include "vtkIdList.h"
#include "vtkIdTypeArray.h"
#include "vtkMPIController.h"
#include "vtkMPIUtilities.h"
#include "vtkMultiProcessController.h"
#include "vtkPUnstructuredGridConnectivity.h"
#include "vtkPointData.h"
#include "vtkPointData.h"
#include "vtkPoints.h"
#include "vtkStructuredData.h"
#include "vtkTimerLog.h"
#include "vtkUnstructuredGrid.h"
#include "vtkUnstructuredGridWriter.h"
#include <sstream>
#include <string>
//#define DEBUG
#include "UnstructuredGhostZonesCommon.h"
//------------------------------------------------------------------------------
// Program main
int main(int argc, char** argv)
{
int rc = 0;
double ellapsed = 0.0;
vtkTimerLog* timer = vtkTimerLog::New();
// STEP 0: Initialize
vtkMPIController* cntrl = vtkMPIController::New();
cntrl->Initialize( &argc, &argv, 0 );
vtkMultiProcessController::SetGlobalController( cntrl );
Rank = cntrl->GetLocalProcessId();
NRanks = cntrl->GetNumberOfProcesses();
// STEP 1: Generate grid in parallel in each process
Grid = vtkUnstructuredGrid::New();
GenerateDataSet();
// STEP 2: Generate ghost zones
vtkMPIUtilities::Printf(cntrl,"[INFO]: Building ghost zones...");
vtkPUnstructuredGridConnectivity* ghostGen =
vtkPUnstructuredGridConnectivity::New();
ghostGen->SetController(cntrl);
ghostGen->RegisterGrid(Grid);
timer->StartTimer();
ghostGen->BuildGhostZoneConnectivity();
timer->StopTimer();
ellapsed = timer->GetElapsedTime();
vtkMPIUtilities::Printf(cntrl,"[DONE]\n");
// get some performance statistics
double minBuildGhostZonesTime = 0.0;
double maxBuildGhostZonesTime = 0.0;
double avgBuildGhostZonesTime = 0.0;
cntrl->Reduce(&ellapsed,&minBuildGhostZonesTime,1,vtkCommunicator::MIN_OP,0);
cntrl->Reduce(&ellapsed,&maxBuildGhostZonesTime,1,vtkCommunicator::MAX_OP,0);
cntrl->Reduce(&ellapsed,&avgBuildGhostZonesTime,1,vtkCommunicator::SUM_OP,0);
avgBuildGhostZonesTime /= static_cast<double>(cntrl->GetNumberOfProcesses());
vtkMPIUtilities::Printf(
cntrl,"-- Ellapsed Time: min=%f, avg=%f, max=%f\n",
minBuildGhostZonesTime,avgBuildGhostZonesTime,maxBuildGhostZonesTime);
// STEP 3: Update ghost zones
std::ostringstream grdfname; // input grid name at each iteration for I/O
std::ostringstream ghostfname; // ghosted grid name at each iteration for I/O
for(int i=0; i < 2; ++i)
{
vtkUnstructuredGrid* ghostGrid = vtkUnstructuredGrid::New();
grdfname.clear();
grdfname.str("");
grdfname << "INITIAL-T" << i;
ghostfname.clear();
ghostfname.str("");
ghostfname << "GHOSTED-T" << i;
// update grid in this time-step
UpdateGrid(i);
#ifdef DEBUG
WriteDataSet(Grid,grdfname.str().c_str());
#endif
vtkMPIUtilities::Printf(cntrl,"[INFO]: iteration=%d\n",i);
vtkMPIUtilities::Printf(cntrl,"[INFO]: Update ghost zones...");
timer->StartTimer();
ghostGen->UpdateGhosts();
timer->StopTimer();
ellapsed = timer->GetElapsedTime();
vtkMPIUtilities::Printf(cntrl,"[DONE]\n");
// get some performance statistics
double minGhostUpdateTime = 0.0;
double maxGhostUpdateTime = 0.0;
double avgGhostUpdateTime = 0.0;
cntrl->Reduce(&ellapsed,&minGhostUpdateTime,1,vtkCommunicator::MIN_OP,0);
cntrl->Reduce(&ellapsed,&maxGhostUpdateTime,1,vtkCommunicator::MAX_OP,0);
cntrl->Reduce(&ellapsed,&avgGhostUpdateTime,1,vtkCommunicator::SUM_OP,0);
avgGhostUpdateTime /= static_cast<double>(cntrl->GetNumberOfProcesses());
vtkMPIUtilities::Printf(
cntrl,"-- Ellapsed Time: min=%f, avg=%f, max=%f\n",
minGhostUpdateTime,avgGhostUpdateTime,maxGhostUpdateTime);
ghostGrid->DeepCopy(ghostGen->GetGhostedGrid());
#ifdef DEBUG
assert("pre: ghost gird should not be NULL!" && (ghostGrid != NULL) );
WriteDataSet(ghostGrid,ghostfname.str().c_str());
#endif
rc += CheckGrid(ghostGrid,i);
ghostGrid->Delete();
}
// STEP 5: Delete the ghost generator
timer->Delete();
ghostGen->Delete();
Grid->Delete();
cntrl->Finalize();
cntrl->Delete();
return( rc );
}
/*=========================================================================
Program: Visualization Toolkit
Module: TestPStructuredGridConnectivity.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 "vtkCellData.h"
#include "vtkCellType.h"
#include "vtkDoubleArray.h"
#include "vtkExtentRCBPartitioner.h"
#include "vtkIdList.h"
#include "vtkIdTypeArray.h"
#include "vtkMPIController.h"
#include "vtkMPIUtilities.h"
#include "vtkMultiProcessController.h"
#include "vtkPUnstructuredGridConnectivity.h"
#include "vtkPUnstructuredGridGhostDataGenerator.h"
#include "vtkPointData.h"
#include "vtkPointData.h"
#include "vtkPoints.h"
#include "vtkStructuredData.h"
#include "vtkTimerLog.h"
#include "vtkUnstructuredGrid.h"
#include "vtkUnstructuredGridWriter.h"
#include <sstream>
#include <string>
//#define DEBUG
#include "UnstructuredGhostZonesCommon.h"
//------------------------------------------------------------------------------
// Program main
int main(int argc, char** argv)
{
int rc = 0;
double ellapsed = 0.0;
vtkTimerLog* timer = vtkTimerLog::New();
// STEP 0: Initialize
vtkMPIController* cntrl = vtkMPIController::New();
cntrl->Initialize( &argc, &argv, 0 );
vtkMultiProcessController::SetGlobalController( cntrl );
Rank = cntrl->GetLocalProcessId();
NRanks = cntrl->GetNumberOfProcesses();
// STEP 1: Generate grid in parallel in each process
Grid = vtkUnstructuredGrid::New();
GenerateDataSet();
// STEP 2: Setup ghost data generator
vtkPUnstructuredGridGhostDataGenerator* ghostGenerator =
vtkPUnstructuredGridGhostDataGenerator::New();
ghostGenerator->SetInputData(Grid);
// STEP 3: Update ghost zones
std::ostringstream grdfname; // input grid name at each iteration for I/O
std::ostringstream ghostfname; // ghosted grid name at each iteration for I/O
for(int i=0; i < 2; ++i)
{
grdfname.clear();
grdfname.str("");
grdfname << "INITIAL-T" << i;
ghostfname.clear();
ghostfname.str("");
ghostfname << "GHOSTED-T" << i;
// update grid in this iteration...
UpdateGrid(i);
Grid->Modified();
#ifdef DEBUG
WriteDataSet(Grid,grdfname.str().c_str());
#endif
// update ghost zones in this iteration...
vtkMPIUtilities::Printf(cntrl,"[INFO]: iteration=%d\n",i);
vtkMPIUtilities::Printf(cntrl,"[INFO]: Update ghost zones...");
timer->StartTimer();
ghostGenerator->Update();
timer->StopTimer();
ellapsed = timer->GetElapsedTime();
vtkMPIUtilities::Printf(cntrl,"[DONE]\n");
// get some performance statistics
double minGhostUpdateTime = 0.0;
double maxGhostUpdateTime = 0.0;
double avgGhostUpdateTime = 0.0;
cntrl->Reduce(&ellapsed,&minGhostUpdateTime,1,vtkCommunicator::MIN_OP,0);
cntrl->Reduce(&ellapsed,&maxGhostUpdateTime,1,vtkCommunicator::MAX_OP,0);
cntrl->Reduce(&ellapsed,&avgGhostUpdateTime,1,vtkCommunicator::SUM_OP,0);
avgGhostUpdateTime /= static_cast<double>(cntrl->GetNumberOfProcesses());
vtkMPIUtilities::Printf(
cntrl,"-- Ellapsed Time: min=%f, avg=%f, max=%f\n",
minGhostUpdateTime,avgGhostUpdateTime,maxGhostUpdateTime);
vtkUnstructuredGrid* ghostGrid = vtkUnstructuredGrid::New();
ghostGrid->DeepCopy(ghostGenerator->GetOutput());
#ifdef DEBUG
assert("pre: ghost gird should not be NULL!" && (ghostGrid != NULL) );
WriteDataSet(ghostGrid,ghostfname.str().c_str());
#endif
rc += CheckGrid(ghostGrid,i);
ghostGrid->Delete();
} // END for
// STEP 5: Delete the ghost generator
timer->Delete();
ghostGenerator->Delete();
Grid->Delete();
cntrl->Finalize();
cntrl->Delete();
return( rc );
}
/*=========================================================================
Program: Visualization Toolkit
Module: vtkPStructuredGridConnectivity.h
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.
=========================================================================*/
#ifndef UNSTRUCTUREDGHOSTZONESCOMMON_H_
#define UNSTRUCTUREDGHOSTZONESCOMMON_H_
// VTK includes
#include "vtkMathUtilities.h"
#include "vtkMPIUtilities.h"
// C/C++ includes
#include <cmath>
#include <iomanip>
#include <sstream>
//------------------------------------------------------------------------------
// G L O B A L D A T A
//------------------------------------------------------------------------------
double Origin[3] = {0.0,0.0,0.0};
double Spacing[3] = {0.5,0.5,0.5};
int Dims[3] = {50,50,50};
int Rank = -1;
int NRanks = 0;
vtkUnstructuredGrid* Grid;
int CheckGrid(vtkUnstructuredGrid* ghostGrid, const int iteration)
{
int rc = 0;
int numOfErrors = 0;
std::ostringstream out;
std::ostringstream err;
vtkMPIController* cntrl =
vtkMPIController::SafeDownCast(
vtkMultiProcessController::GetGlobalController());
// check node fields by the iteration number
vtkDoubleArray* nodeXYZ =
vtkDoubleArray::SafeDownCast(
ghostGrid->GetPointData()->GetArray("NodeXYZ"));
assert("pre: nodeXYZ != NULL" && (nodeXYZ != NULL) );
assert("pre: nodeXYZ numtuples mismatch!" &&
(ghostGrid->GetNumberOfPoints()==nodeXYZ->GetNumberOfTuples()));
assert("pre: nodeXYZ numcomponents mismatch!" &&
(nodeXYZ->GetNumberOfComponents()==3));
out.str("");
err.str("");
double pnt[3];
double* ptr = static_cast<double*>(nodeXYZ->GetVoidPointer(0));
for(vtkIdType nodeIdx=0; nodeIdx < ghostGrid->GetNumberOfPoints(); ++nodeIdx)
{
ghostGrid->GetPoint(nodeIdx,pnt);
for(int dim=0; dim < 3; ++dim)
{
double actual = ptr[nodeIdx*3+dim];
double expected = pnt[dim]+static_cast<double>(iteration);
if( ! vtkMathUtilities::NearlyEqual<double>(actual,expected) )
{
++numOfErrors;
++rc;
err << std::setprecision(5)
<< "\t[ERROR]: value mismatch at node=" << nodeIdx
<< " expected=" << expected
<< " actual=" << actual
<< " delta=" << std::fabs(actual-expected)
<< std::endl;
} // END if
} // END for all dimensions
} // END for all nodes
out << "[INFO]: " << numOfErrors << "/" << ghostGrid->GetNumberOfPoints()
<< " nodes appear wrong: " << std::endl;
out << err.str();
vtkMPIUtilities::SynchronizedPrintf(cntrl,"%s",out.str().c_str());
// likewise, check cell-fields
vtkDoubleArray* cellXYZ =
vtkDoubleArray::SafeDownCast(
ghostGrid->GetCellData()->GetArray("CentroidXYZ"));
assert("pre: nodeXYZ numtuples mismatch!" &&
(ghostGrid->GetNumberOfCells()==cellXYZ->GetNumberOfTuples()));
assert("pre: nodeXYZ numcomponents mismatch!" &&
(cellXYZ->GetNumberOfComponents()==3));
numOfErrors = 0;
out.str("");
err.str("");
double centroid[3];
double* cptr = static_cast<double*>(cellXYZ->GetVoidPointer(0));
vtkIdList* ptIds = vtkIdList::New();
for(vtkIdType cellIdx=0; cellIdx < ghostGrid->GetNumberOfCells(); ++cellIdx)
{
ghostGrid->GetCellPoints(cellIdx,ptIds);
assert("pre: numpoints per cell must be 8!" &&
ptIds->GetNumberOfIds()==8);
centroid[0] = 0.0;
centroid[1] = 0.0;
centroid[2] = 0.0;
for(vtkIdType n=0; n < ptIds->GetNumberOfIds(); ++n)
{
vtkIdType idx = ptIds->GetId(n);
ghostGrid->GetPoint(idx,pnt);
centroid[0] += pnt[0];
centroid[1] += pnt[1];
centroid[2] += pnt[2];
} // END for all cell nodes
centroid[0] /= static_cast<double>(ptIds->GetNumberOfIds());
centroid[1] /= static_cast<double>(ptIds->GetNumberOfIds());
centroid[2] /= static_cast<double>(ptIds->GetNumberOfIds());
for(int dim=0; dim < 3; ++dim)
{
double actual = cptr[cellIdx*3+dim];
double expected = centroid[dim]+static_cast<double>(iteration);
if( ! vtkMathUtilities::NearlyEqual<double>(actual,expected) )
{
++numOfErrors;
++rc;
err << std::setprecision(5)
<< "\t[ERROR]: cell value mismatch at cell=" << cellIdx
<< " dimension=" << dim
<< " expected=" << expected
<< " actual=" << actual
<< " delta= " << std::fabs(actual-expected)
<< std::endl;
} // END if
} // END for all dimensions
} // END for all cells
out << "[INFO]: " << numOfErrors << "/" << ghostGrid->GetNumberOfCells()
<< " cells appear wrong: " << std::endl;
out << err.str();
vtkMPIUtilities::SynchronizedPrintf(cntrl,"%s",out.str().c_str());
ptIds->Delete();
return( rc );
}
//------------------------------------------------------------------------------
void UpdateGrid(const int iteration)
{
// increment node fields by the iteration number
vtkDoubleArray* nodeXYZ =
vtkDoubleArray::SafeDownCast(
Grid->GetPointData()->GetArray("NodeXYZ"));
assert("pre: nodeXYZ != NULL" && (nodeXYZ != NULL) );
assert("pre: nodeXYZ numtuples mismatch!" &&
(Grid->GetNumberOfPoints()==nodeXYZ->GetNumberOfTuples()));
assert("pre: nodeXYZ numcomponents mismatch!" &&
(nodeXYZ->GetNumberOfComponents()==3));
double* ptr = static_cast<double*>(nodeXYZ->GetVoidPointer(0));
for(vtkIdType nodeIdx=0; nodeIdx < Grid->GetNumberOfPoints(); ++nodeIdx)
{
ptr[ nodeIdx*3 ] += static_cast<double>(iteration);
ptr[ nodeIdx*3+1 ] += static_cast<double>(iteration);
ptr[ nodeIdx*3+2 ] += static_cast<double>(iteration);
} // END for all nodes
// increment cell fields by the iteration number
vtkDoubleArray* cellXYZ =
vtkDoubleArray::SafeDownCast(
Grid->GetCellData()->GetArray("CentroidXYZ"));
assert("pre: nodeXYZ numtuples mismatch!" &&
(Grid->GetNumberOfCells()==cellXYZ->GetNumberOfTuples()));
assert("pre: nodeXYZ numcomponents mismatch!" &&
(cellXYZ->GetNumberOfComponents()==3));
double* cptr = static_cast<double*>(cellXYZ->GetVoidPointer(0));
for(vtkIdType cellIdx=0; cellIdx < Grid->GetNumberOfCells(); ++cellIdx)
{
cptr[ cellIdx*3 ] += static_cast<double>(iteration);
cptr[ cellIdx*3+1 ] += static_cast<double>(iteration);
cptr[ cellIdx*3+2 ] += static_cast<double>(iteration);
} // END for all cells
}
//------------------------------------------------------------------------------
void SetXYZCellField()
{
vtkDoubleArray* centerXYZ = vtkDoubleArray::New();
centerXYZ->SetName("CentroidXYZ");
centerXYZ->SetNumberOfComponents(3);
centerXYZ->SetNumberOfTuples( Grid->GetNumberOfCells() );
double* ptr = static_cast<double*>(centerXYZ->GetVoidPointer(0));
double centroid[3];
vtkIdList* ptIds = vtkIdList::New();
for(vtkIdType cell=0; cell < Grid->GetNumberOfCells(); ++cell)
{
centroid[0] = centroid[1] = centroid[2] = 0.0;
Grid->GetCellPoints(cell,ptIds);
for(vtkIdType n=0; n < ptIds->GetNumberOfIds(); ++n)
{
centroid[0] += Grid->GetPoint(ptIds->GetId(n))[0];
centroid[1] += Grid->GetPoint(ptIds->GetId(n))[1];
centroid[2] += Grid->GetPoint(ptIds->GetId(n))[2];
} // END for all cell nodes
centroid[0] /= static_cast<double>(ptIds->GetNumberOfIds());
centroid[1] /= static_cast<double>(ptIds->GetNumberOfIds());
centroid[2] /= static_cast<double>(ptIds->GetNumberOfIds());
memcpy(&ptr[cell*3],centroid,3*sizeof(double));
} // END for all cells
Grid->GetCellData()->AddArray( centerXYZ );
centerXYZ->Delete();
ptIds->Delete();
}
//------------------------------------------------------------------------------
void SetXYZNodeField()
{
vtkDoubleArray* nodeXYZ = vtkDoubleArray::New();
nodeXYZ->SetName("NodeXYZ");
nodeXYZ->SetNumberOfComponents(3);
nodeXYZ->SetNumberOfTuples( Grid->GetNumberOfPoints() );