Commit af4717da authored by cyrush's avatar cyrush
Browse files

Added ghost zone neighbors enhancement to parallel connected components

Added an algorithm selection option to the gradient expression
Added a new gradient algorithm "nzqh"
Updated macro expressions that use the gradient expression to all for algorithm selection
Fixed case where duplicate points where created in the rectilinear case of VisIt's clipper. 
Fixed var and mesh naming issue with Silo export




git-svn-id: http://visit.ilight.com/svn/visit/trunk/src@1887 18c085ea-50e0-402c-830e-de6fd14e8384
parent c584a0ed
......@@ -59,6 +59,7 @@
#include <vtkDataSetRemoveGhostCells.h>
#include <vtkIntArray.h>
#include <vtkPointData.h>
#include <vtkUnsignedCharArray.h>
#include <vtkUnstructuredGrid.h>
#include <vtkUnstructuredGridRelevantPointsFilter.h>
#include <vtkUnstructuredGridWriter.h>
......@@ -143,6 +144,9 @@ avtConnComponentsExpression::GetNumberOfComponents()
// Cyrus Harrison, Fri Mar 16 10:46:28 PDT 2007
// Added progress update.
//
// Cyrus Harrison, Fri Mar 16 10:46:28 PDT 2007
// Added extraction of ghost zone neighbors for the parallel case
//
// ****************************************************************************
void
avtConnComponentsExpression::Execute()
......@@ -154,8 +158,6 @@ avtConnComponentsExpression::Execute()
// get input data tree to obtain datasets
avtDataTree_p tree = GetInputDataTree();
// holds number of datasets
int nsets;
......@@ -166,18 +168,50 @@ avtConnComponentsExpression::Execute()
vector<int> domain_ids;
tree->GetAllDomainIds(domain_ids);
// check for ghosts
bool have_ghosts = false;
if(nsets > 0)
have_ghosts = data_sets[0]->GetCellData()->GetArray("avtGhostZones");
// set progress related vars
#ifdef PARALLEL
totalSteps = nsets *4;
if(have_ghosts)
totalSteps+= nsets;
#else
totalSteps = nsets *2;
#endif
currentProgress = 0;
#ifdef PARALLEL
if(have_ghosts)
{
// if we have ghosts, label ghost neighbors for reduced comm in global
// resolve
for( i = 0; i <nsets ; i++)
{
LabelGhostNeighbors(data_sets[i]);
UpdateProgress(currentProgress++,totalSteps);
}
}
#endif
// filter out any ghost cells
vtkDataSetRemoveGhostCells **ghost_filters =
new vtkDataSetRemoveGhostCells*[nsets];
for( i = 0 ; i < nsets ; i++)
vtkDataSetRemoveGhostCells **ghost_filters = NULL;
if(have_ghosts)
{
ghost_filters[i] = vtkDataSetRemoveGhostCells::New();
ghost_filters[i]->SetInput(data_sets[i]);
ghost_filters[i]->Update();
vtkDataSet *ds_filtered = ghost_filters[i]->GetOutput();
ds_filtered->Update();
data_sets[i] = ds_filtered;
ghost_filters = new vtkDataSetRemoveGhostCells*[nsets];
for( i = 0 ; i < nsets ; i++)
{
ghost_filters[i] = vtkDataSetRemoveGhostCells::New();
ghost_filters[i]->SetInput(data_sets[i]);
ghost_filters[i]->Update();
vtkDataSet *ds_filtered = ghost_filters[i]->GetOutput();
ds_filtered->Update();
data_sets[i] = ds_filtered;
}
}
// array to hold output sets
......@@ -196,28 +230,12 @@ avtConnComponentsExpression::Execute()
int num_local_comps=0;
int num_comps = 0;
// set progress related vars
#ifdef PARALLEL
totalSteps = nsets *4;
#else
totalSteps = nsets *2;
#endif
currentProgress = 0;
// process all local sets
for(i = 0; i < nsets ; i++)
{
// get current set
vtkDataSet *curr_set = data_sets[i];
// Make sure the input dataset is an unstructured grid
if(curr_set->GetDataObjectType() != VTK_UNSTRUCTURED_GRID)
{
EXCEPTION1(ExpressionException,
"Connected Components Expression requires unstructured mesh data");
}
// perform connected components labeling on current set
// (this only resolves components within the set)
vtkIntArray *res_array = SingleSetLabel(curr_set,results_num_comps[i]);
......@@ -233,7 +251,6 @@ avtConnComponentsExpression::Execute()
vtkDataSet *res_set = (vtkDataSet *) curr_set->NewInstance();
res_set->ShallowCopy(curr_set);
// add array to dataset
res_set->GetCellData()->AddArray(res_array);
......@@ -317,10 +334,12 @@ avtConnComponentsExpression::Execute()
// dec ref pointer for each set's label array
result_arrays[i]->Delete();
// cleanup ghost filters
ghost_filters[i]->Delete();
if(have_ghosts)
ghost_filters[i]->Delete();
}
// cleanup ghost filters array
delete [] ghost_filters;
if(have_ghosts)
delete [] ghost_filters;
// Set progress to complete
UpdateProgress(totalSteps,totalSteps);
......@@ -329,6 +348,79 @@ avtConnComponentsExpression::Execute()
nFinalComps = num_comps;
}
// ****************************************************************************
// Method: avtConnComponentsExpression::LabelGhostNeighbors
//
// Purpose:
// Identifies cells that have ghost neighbors, storing the info in
// a vtkUnsignedCharArray named "avtGhostZoneNeighbors".
//
// Arguments:
// data_set Input mesh
//
// Programmer: Cyrus Harrison
// Creation: August 11, 2007
//
// ****************************************************************************
void
avtConnComponentsExpression::LabelGhostNeighbors(vtkDataSet *data_set)
{
// loop indices
int i,j,k;
vtkUnsignedCharArray *gz_array = (vtkUnsignedCharArray *) data_set
->GetCellData()->GetArray("avtGhostZones");
// if the data set does not have ghosts, we are done
if (!gz_array)
return;
unsigned char *gz_ptr = (unsigned char *)gz_array->GetPointer(0);
int ncells = data_set->GetNumberOfCells();
vtkUnsignedCharArray *gzn_array = vtkUnsignedCharArray::New();
gzn_array->SetName("avtGhostZoneNeighbors");
gzn_array->SetNumberOfComponents(1);
gzn_array->SetNumberOfTuples(ncells);
unsigned char *gzn_ptr = (unsigned char *)gzn_array->GetPointer(0);
// init the ghost zone neighbors array
memset(gzn_ptr,0,ncells * sizeof(unsigned char));
for ( i=0; i < ncells; i++)
{
// if this cell has ghost zones, label it's neighbors
if(gz_ptr[i])
{
// get cell neighbors
vtkIdList *gcell_pts = data_set->GetCell(i)->GetPointIds();
int ngcell_pts = gcell_pts->GetNumberOfIds();
for( j=0; j < ngcell_pts; j++)
{
// neighbors share points with the current cell
vtkIdList *gpt = vtkIdList::New();
gpt->SetNumberOfIds(1);
gpt->SetId(0,gcell_pts->GetId(j));
vtkIdList *nei_cells = vtkIdList::New();
data_set->GetCellNeighbors(i,gpt,nei_cells);
int nnei = nei_cells->GetNumberOfIds();
// tag neighbors
for ( k = 0; k < nnei; k++)
gzn_array->SetTuple1(nei_cells->GetId(k),1);
gpt->Delete();
nei_cells->Delete();
}
}
}
data_set->GetCellData()->AddArray(gzn_array);
gzn_array->Delete();
}
// ****************************************************************************
// Method: avtConnComponentsExpression::SingleSetLabel
//
......@@ -997,6 +1089,25 @@ avtConnComponentsExpression::ShiftLabels(vtkIntArray *labels, int shift)
}
// ****************************************************************************
// Method: avtGradientFilter::PerformRestriction
//
// Purpose:
// Request ghost zones.
//
// Programmer: Cyrus Harrison
// Creation: August 11, 2007
//
// ****************************************************************************
avtPipelineSpecification_p
avtConnComponentsExpression::PerformRestriction(avtPipelineSpecification_p in_spec)
{
avtPipelineSpecification_p spec =
avtExpressionFilter::PerformRestriction(in_spec);
spec->GetDataSpecification()->SetDesiredGhostDataType(GHOST_ZONE_DATA);
return spec;
}
// ****************************************************************************
// Method: UnionFind constructor
//
......@@ -1516,7 +1627,6 @@ avtConnComponentsExpression::BoundarySet::GetBounds(double *bounds) const
// Creation: February 16, 2007
//
// Modifications:
//
// Cyrus Harrison, Thu Mar 1 07:49:44 PST 2007
// Added case for empty input datasets.
//
......@@ -1781,6 +1891,10 @@ avtConnComponentsExpression::BoundarySet::GetBoundsIntersection
// Programmer: Cyrus Harrison
// Creation: February 2, 2007
//
// Modifications:
// Cyrus Harrison, Sat Aug 11 15:08:03 PDT 2007
// Added ghost neighbors optimization
//
// ****************************************************************************
void
avtConnComponentsExpression::BoundarySet::RelocateUsingPartition
......@@ -1821,7 +1935,14 @@ avtConnComponentsExpression::BoundarySet::RelocateUsingPartition
// find which cells from the current mesh to send
vtkUnstructuredGrid *mesh = (vtkUnstructuredGrid*) meshes[i];
int ncells = mesh->GetNumberOfCells();
vtkIntArray *labels = (vtkIntArray*)mesh->GetCellData()->GetArray("ccl");
// get ghost zone neighbor info if available
vtkUnsignedCharArray *gzn_array = (vtkUnsignedCharArray *) mesh
->GetCellData()->GetArray("avtGhostZoneNeighbors");
unsigned char *gzn_ptr = NULL;
if (gzn_array)
gzn_ptr = (unsigned char *)gzn_array->GetPointer(0);
// initialize temporary mesh and cell count holder
for( j = 0; j < nprocs; j++)
......@@ -1833,6 +1954,11 @@ avtConnComponentsExpression::BoundarySet::RelocateUsingPartition
// for each cell
for( j = 0; j < ncells; j++)
{
// if we have ghost neighbors labeled, we only need to send them
if(gzn_ptr)
if(gzn_ptr[j] != 1)
continue;
// find which processors need this cell
vtkCell *cell = mesh->GetCell(j);
spart.GetProcessorList(cell,proc_list);
......@@ -2386,6 +2512,9 @@ PartitionBoundary::AttemptSplit(PartitionBoundary *&b1, PartitionBoundary *&b2)
// Programmer: Cyrus Harrison
// Creation: February 2, 2007
//
// Modifications:
// Cyrus Harrison, Sat Aug 11 15:08:37 PDT 2007
// Added ghost neighbors optimization
//
// ****************************************************************************
......@@ -2472,8 +2601,23 @@ avtConnComponentsExpression::SpatialPartition::CreatePartition
const int ncells = meshes[i]->GetNumberOfCells();
double bbox[6];
double pt[3];
// get ghost zone neighbor info if available
vtkUnsignedCharArray *gzn_array = (vtkUnsignedCharArray *) meshes[i]
->GetCellData()->GetArray("avtGhostZoneNeighbors");
unsigned char *gzn_ptr = NULL;
if (gzn_array)
gzn_ptr = (unsigned char *)gzn_array->GetPointer(0);
for (j = 0 ; j < ncells ; j++)
{
// if we have ghost neighbors labled, we only need to send them
if(gzn_ptr)
if(gzn_ptr[j] !=1)
continue;
vtkCell *cell = meshes[i]->GetCell(j);
cell->GetBounds(bbox);
pt[0] = (bbox[0] + bbox[1]) / 2.;
......
......@@ -64,6 +64,9 @@ class avtIntervalTree;
// Cyrus Harrison, Fri Mar 16 15:52:47 PDT 2007
// Added variables to track progress.
//
// Cyrus Harrison, Sat Aug 11 14:41:01 PDT 2007
// Added LabelGhostNeighbors and PerformRestriction
//
// ****************************************************************************
class EXPRESSION_API avtConnComponentsExpression : public avtExpressionFilter
......@@ -203,7 +206,12 @@ class EXPRESSION_API avtConnComponentsExpression : public avtExpressionFilter
int totalSteps;
virtual void Execute(void);
virtual avtPipelineSpecification_p
PerformRestriction(avtPipelineSpecification_p);
virtual void LabelGhostNeighbors(vtkDataSet *);
virtual vtkIntArray *SingleSetLabel(vtkDataSet *, int &);
virtual int MultiSetResolve(int,
......
......@@ -41,8 +41,7 @@
#include <avtCurlFilter.h>
#include <stdio.h>
#include <snprintf.h>
#include <ExpressionException.h>
......@@ -95,6 +94,9 @@ avtCurlFilter::~avtCurlFilter()
// If we are creating a scalar, then make sure the expression type is a
// scalar as well.
//
// Cyrus Harrison, Sat Aug 11 18:15:50 PDT 2007
// Add second argument for gradient algorithm selection
//
// ****************************************************************************
void
......@@ -111,20 +113,66 @@ avtCurlFilter::GetMacro(std::vector<std::string> &args, std::string &ne,
}
}
char new_expr[1024];
int nargs = args.size();
char new_expr[2048];
if (do3D)
{
sprintf(new_expr, "{gradient(%s[2])[1]-gradient(%s[1])[2],"
"gradient(%s[0])[2]-gradient(%s[2])[0],"
"gradient(%s[1])[0]-gradient(%s[0])[1]}",
args[0].c_str(), args[0].c_str(), args[0].c_str(),
args[0].c_str(), args[0].c_str(), args[0].c_str());
if(nargs == 1)
{
SNPRINTF(new_expr,2048,
"{gradient(%s[2])[1]-gradient(%s[1])[2],"
"gradient(%s[0])[2]-gradient(%s[2])[0],"
"gradient(%s[1])[0]-gradient(%s[0])[1]}",
args[0].c_str(), args[0].c_str(), args[0].c_str(),
args[0].c_str(), args[0].c_str(), args[0].c_str());
}
else if(nargs > 1)
{
SNPRINTF(new_expr,2048,
"{gradient(%s[2],%s)[1]-gradient(%s[1],%s)[2],"
"gradient(%s[0],%s)[2]-gradient(%s[2],%s)[0],"
"gradient(%s[1],%s)[0]-gradient(%s[0],%s)[1]}",
args[0].c_str(), args[1].c_str(),
args[0].c_str(), args[1].c_str(),
args[0].c_str(), args[1].c_str(),
args[0].c_str(), args[1].c_str(),
args[0].c_str(), args[1].c_str(),
args[0].c_str(), args[1].c_str());
}
else
{
EXCEPTION1(ExpressionException, " invalid curl syntax. "
"Expected arguments: "
"vector_var, gradient_algorithm\n"
"[gradient_algorithm is optional]");
}
type = Expression::VectorMeshVar;
}
else
{
sprintf(new_expr, "gradient(%s[1])[0]-gradient(%s[0])[1]",
args[0].c_str(), args[0].c_str());
if(nargs == 1)
{
SNPRINTF(new_expr,2048,
"gradient(%s[1])[0]-gradient(%s[0])[1]",
args[0].c_str(), args[0].c_str());
}
else if(nargs > 1)
{
SNPRINTF(new_expr,2048,
"gradient(%s[1],%s)[0]-gradient(%s[0],%s)[1]",
args[0].c_str(), args[1].c_str(),
args[0].c_str(), args[1].c_str());
}
else
{
EXCEPTION1(ExpressionException, " invalid curl syntax. "
"Expected arguments: "
"vector_var, gradient_algorithm\n"
"[gradient_algorithm is optional]");
}
type = Expression::ScalarMeshVar;
}
ne = new_expr;
......
......@@ -41,8 +41,8 @@
#include <avtDivergenceFilter.h>
#include <stdio.h>
#include <snprintf.h>
#include <ExpressionException.h>
// ****************************************************************************
// Method: avtDivergenceFilter constructor
......@@ -87,6 +87,9 @@ avtDivergenceFilter::~avtDivergenceFilter()
// Hank Childs, Mon Jun 6 11:21:23 PDT 2005
// Add support for 2D.
//
// Cyrus Harrison, Sat Aug 11 18:34:41 PDT 2007
// Add second argument for gradient algorithm selection
//
// ****************************************************************************
void
......@@ -101,14 +104,59 @@ avtDivergenceFilter::GetMacro(std::vector<std::string> &args, std::string &ne,
do3D = (atts.GetTopologicalDimension() == 3);
}
char new_expr[1024];
int nargs = args.size();
char new_expr[2048];
if (do3D)
sprintf(new_expr, "gradient(%s[0])[0]+gradient(%s[1])[1]+"
"gradient(%s[2])[2]",
args[0].c_str(), args[0].c_str(), args[0].c_str());
{
if(nargs == 1)
{
SNPRINTF(new_expr, 2048,
"gradient(%s[0])[0]+gradient(%s[1])[1]+"
"gradient(%s[2])[2]",
args[0].c_str(), args[0].c_str(), args[0].c_str());
}
else if(nargs > 1)
{
SNPRINTF(new_expr, 2048,
"gradient(%s[0],%s)[0]+gradient(%s[1],%s)[1]+"
"gradient(%s[2],%s)[2]",
args[0].c_str(), args[1].c_str(),
args[0].c_str(), args[1].c_str(),
args[0].c_str(), args[1].c_str());
}
else
{
EXCEPTION1(ExpressionException, " invalid divergence syntax. "
"Expected arguments: "
" vector_var, gradient_algorithm\n"
" gradient_algorithm is optional");
}
}
else
sprintf(new_expr, "gradient(%s[0])[0]+gradient(%s[1])[1]",
args[0].c_str(), args[0].c_str());
{
if(nargs == 1)
{
SNPRINTF(new_expr, 2048,
"gradient(%s[0])[0]+gradient(%s[1])[1]",
args[0].c_str(), args[0].c_str());
}
else if(nargs > 1)
{
SNPRINTF(new_expr, 2048,
"gradient(%s[0],%s)[0]+gradient(%s[1],%s)[1]",
args[0].c_str(), args[1].c_str(),
args[0].c_str(), args[1].c_str());
}
else
{
EXCEPTION1(ExpressionException, " invalid divergence syntax. "
"Expected arguments: "
" vector_var, gradient_algorithm\n"
" gradient_algorithm is optional");
}
}
ne = new_expr;
type = Expression::ScalarMeshVar;
}
......
......@@ -47,7 +47,7 @@
#include <vtkCell.h>
#include <vtkCellData.h>
#include <vtkCellDataToPointData.h>
#include <vtkCellDerivatives.h>
#include <vtkCellDerivatives.h>
#include <vtkDataSet.h>
#include <vtkFloatArray.h>
#include <vtkIdList.h>
......@@ -56,6 +56,7 @@
#include <vtkRectilinearGrid.h>
#include <vtkStructuredGrid.h>
#include <avtExprNode.h>
#include <avtCallback.h>
#include <DebugStream.h>
......@@ -79,11 +80,14 @@
// Hank Childs, Mon Feb 13 15:08:41 PST 2006
// Add support for logical gradients.
//
// Cyrus Harrison, Wed Aug 8 11:19:10 PDT 2007
// Add support for multiple gradient algorithms.
//
// ****************************************************************************
avtGradientFilter::avtGradientFilter()
{
doLogicalGradients = false;
gradientAlgo = SAMPLE;
}
......@@ -109,6 +113,112 @@ avtGradientFilter::~avtGradientFilter()
}
// ****************************************************************************
// Method: avtGradientFilter::ProcessArguments
//
// Purpose:
// Parses optional algorithm argument.
//
// Arguments:
// args Expression arguments
// state Expression pipeline state
//
// Programmer: Cyrus Harrison
// Creation: August 8, 2007
//
// ****************************************************************************
void
avtGradientFilter::ProcessArguments(ArgsExpr *args,
ExprPipelineState *state)
{
// get the argument list and # of arguments
std::vector<ArgExpr*> *arguments = args->GetArgs();
int nargs = arguments->size();
// check for call with no args
if (nargs == 0)
{
EXCEPTION1(ExpressionException,
"gradient() Incorrect syntax.\n"
" usage: gradient(varname,algo)\n"
" The algo parameter is optional "
"and specifies which gradient algorithm is used.\n"
"Valid Options:\n"
" type: 0,1,2 or \"sample\",\"logical\",\"nzqh\"\n"
"(Default: algo = sample)");
}
// first argument is the var name, let it do its own magic
ArgExpr *first_arg = (*arguments)[0];
avtExprNode *first_tree = dynamic_cast<avtExprNode*>(first_arg->GetExpr());
first_tree->CreateFilters(state);
// Check to see if the user passed in the 2nd argument (optional)
// that specifies the gradient algo
// Options:
// sample (0) - calc by sampling in x,y(,z) dir
// logical (1) - calc on logical mesh
// nzqh (2) - nodal to zonal logical gradient for strucuted grids made