Commit 49db7bb2 authored by hrchilds's avatar hrchilds
Browse files

Update from February 14, 2006

git-svn-id: http://visit.ilight.com/svn/visit/trunk/src@619 18c085ea-50e0-402c-830e-de6fd14e8384
parent c46cb1d2
......@@ -4,6 +4,8 @@
#include <avtSIL.h>
#include <algorithm>
#include <avtSILCollection.h>
#include <avtSILNamespace.h>
#include <avtSILSet.h>
......
......@@ -5183,6 +5183,12 @@ avtGenericDatabase::ActivateTimestep(int stateIndex)
// Hank Childs, Mon Jun 27 16:24:23 PDT 2005
// Add argument to GetDataset.
//
// Hank Childs, Thu Feb 9 16:24:32 PST 2006
// If we have domain boundary information and we are doing reconstruction
// for mixed material zones, then we need to read in the materials for
// every domain, so the ghost zone communication routines can work
// correctly.
//
// ****************************************************************************
void
......@@ -5267,11 +5273,21 @@ avtGenericDatabase::ReadDataset(avtDatasetCollection &ds, intVector &domains,
src->DatabaseProgress(0, 0, progressString);
int nDomains = domains.size();
avtSILRestrictionTraverser trav(silr);
bool forceMIR = spec->MustDoMaterialInterfaceReconstruction();
if (spec->NeedMixedVariableReconstruction())
{
avtDatasetCollection emptyCollection(0);
intVector emptyDomainList;
avtDomainBoundaries *dbi = GetDomainBoundaryInformation(emptyCollection,
emptyDomainList, spec);
if (dbi != NULL)
forceMIR = true;
}
for (i = 0 ; i < nDomains ; i++)
{
stringVector labels;
stringVector matnames;
bool forceMIR = spec->MustDoMaterialInterfaceReconstruction();
bool doSelect = PrepareMaterialSelect(domains[i], forceMIR, trav,
matnames);
int nmats = matnames.size();
......
......@@ -17,6 +17,9 @@
#include <vtkPointData.h>
#include <vtkPointDataToCellData.h>
#include <vtkRectilinearGrid.h>
#include <vtkStructuredGrid.h>
#include <avtCallback.h>
#include <DebugStream.h>
#include <ExpressionException.h>
......@@ -36,11 +39,14 @@
// Hank Childs, Fri Mar 4 08:21:04 PST 2005
// Removed centering conversion modules.
//
// Hank Childs, Mon Feb 13 15:08:41 PST 2006
// Add support for logical gradients.
//
// ****************************************************************************
avtGradientFilter::avtGradientFilter()
{
;
doLogicalGradients = false;
}
......@@ -66,6 +72,25 @@ avtGradientFilter::~avtGradientFilter()
}
// ****************************************************************************
// Method: avtGradientFilter::PreExecute
//
// Purpose:
// Initializes a flag saying whether or not we've issued a warning.
//
// Programmer: Hank Childs
// Creation: February 13, 2006
//
// ****************************************************************************
void
avtGradientFilter::PreExecute(void)
{
avtSingleInputExpressionFilter::PreExecute();
haveIssuedWarning = false;
}
// ****************************************************************************
// Method: avtGradientFilter::DeriveVariable
//
......@@ -106,6 +131,9 @@ avtGradientFilter::~avtGradientFilter()
// Hank Childs, Fri Mar 11 16:01:21 PST 2005
// Fix memory leak.
//
// Hank Childs, Mon Feb 13 15:08:41 PST 2006
// Add support for logical gradients ['4385].
//
// ****************************************************************************
vtkDataArray *
......@@ -115,6 +143,30 @@ avtGradientFilter::DeriveVariable(vtkDataSet *in_ds)
{
return RectilinearGradient((vtkRectilinearGrid *) in_ds);
}
if (doLogicalGradients)
{
if (in_ds->GetDataObjectType() != VTK_STRUCTURED_GRID)
{
if (!haveIssuedWarning)
avtCallback::IssueWarning("You can only do logical gradients "
"on structured grids.");
haveIssuedWarning = true;
int nvals = 0;
if (in_ds->GetPointData()->GetScalars() != NULL)
nvals = in_ds->GetNumberOfPoints();
else
nvals = in_ds->GetNumberOfCells();
vtkFloatArray *rv = vtkFloatArray::New();
rv->SetNumberOfComponents(3);
rv->SetNumberOfTuples(nvals);
float vals[3] = { 0., 0., 0. };
for (int i = 0 ; i < nvals ; i++)
rv->SetTuple(i, vals);
return rv;
}
return LogicalGradient((vtkStructuredGrid *) in_ds);
}
vtkDataArray *scalarValues = in_ds->GetPointData()->GetScalars();
bool recentered = false;
......@@ -395,6 +447,9 @@ float avtGradientFilter::EvaluateValue(float x, float y, float z,
// Jeremy Meredith, Fri Jul 2 15:58:01 PDT 2004
// Added a check to make sure the scalars existed before proceeding.
//
// Hank Childs, Mon Feb 13 16:12:23 PST 2006
// Removed misleading comment.
//
// ****************************************************************************
vtkDataArray *
......@@ -495,7 +550,6 @@ avtGradientFilter::RectilinearGradient(vtkRectilinearGrid *rg)
const int kskip = dims0*dims1;
if (dims2 <= 1)
{
// Do all the cases where we have valid values on both sides.
for (j = 0 ; j < dims1 ; j++)
{
for (i = 0 ; i < dims0 ; i++)
......@@ -524,7 +578,6 @@ avtGradientFilter::RectilinearGradient(vtkRectilinearGrid *rg)
}
else
{
// Do all the cases where we have valid values on both sides.
for (k = 0 ; k < dims2 ; k++)
{
for (j = 0 ; j < dims1 ; j++)
......@@ -577,3 +630,278 @@ avtGradientFilter::RectilinearGradient(vtkRectilinearGrid *rg)
}
// ****************************************************************************
// Method: avtGradientFilter::LogicalGradient
//
// Purpose:
// Determines the logical gradient of a curvilinear dataset.
//
// Programmer: Hank Childs
// Creation: February 13, 2006
//
// ****************************************************************************
vtkDataArray *
avtGradientFilter::LogicalGradient(vtkStructuredGrid *sg)
{
int i, j, k;
vtkPoints *vtkpts = sg->GetPoints();
float *pts = (float *) vtkpts->GetVoidPointer(0);
bool deletePoints = false;
int dims[3];
sg->GetDimensions(dims);
bool isNodal = true;
vtkDataArray *s = sg->GetPointData()->GetScalars();
if (s == NULL)
{
s = sg->GetCellData()->GetScalars();
if (s == NULL)
{
EXCEPTION1(ExpressionException, "the scalar variable could not"
" be found.");
}
isNodal = false;
dims[0] -= 1;
dims[1] -= 1;
dims[2] -= 1;
float *pts2 = NULL;
if (dims[2] > 1)
{
pts2 = new float[3*dims[0]*dims[1]*dims[2]];
for (k = 0 ; k < dims[2] ; k++)
{
for (j = 0 ; j < dims[1] ; j++)
{
for (i = 0 ; i < dims[0] ; i++)
{
int c_idx = k*(dims[0])*(dims[1]) + j*dims[0] + i;
float *p = pts2 + 3*c_idx;
p[0] = 0.;
p[1] = 0.;
p[2] = 0.;
for (int l = 0 ; l < 8 ; l++)
{
int p_idx = (l&1 ? i+1 : i);
p_idx += (l&2 ? j+1 : j)*(dims[0]+1);
p_idx += (l&4 ? k+1 : k)*(dims[0]+1)*(dims[1]+1);
p[0] += pts[3*p_idx] / 8.;
p[1] += pts[3*p_idx+1] / 8.;
p[2] += pts[3*p_idx+2] / 8.;
}
}
}
}
}
else
{
pts2 = new float[3*dims[0]*dims[1]];
for (j = 0 ; j < dims[1] ; j++)
{
for (i = 0 ; i < dims[0] ; i++)
{
int c_idx = j*dims[0] + i;
float *p = pts2 + 3*c_idx;
p[0] = 0.;
p[1] = 0.;
p[2] = 0.;
for (int l = 0 ; l < 4 ; l++)
{
int p_idx = (l&1 ? i+1 : i);
p_idx += (l&2 ? j+1 : j)*(dims[0]+1);
p[0] += pts[3*p_idx] / 4.;
p[1] += pts[3*p_idx+1] / 4.;
}
}
}
}
pts = pts2;
deletePoints = true;
}
vtkDataArray *out_array = s->NewInstance();
out_array->SetNumberOfComponents(3);
out_array->SetNumberOfTuples(s->GetNumberOfTuples());
float *in = (float *) s->GetVoidPointer(0);
float *out = (float *) out_array->GetVoidPointer(0);
const int dims0 = dims[0];
const int dims1 = dims[1];
const int dims2 = dims[2];
const int iskip = 1;
const int jskip = dims0;
const int kskip = dims0*dims1;
if (dims2 <= 1)
{
for (j = 0 ; j < dims1 ; j++)
{
for (i = 0 ; i < dims0 ; i++)
{
int index = j*jskip + i*iskip;
int vec_index = 3*index;
float *pt1 = pts + 3*(index+iskip);
float *pt2 = pts + 3*(index-iskip);
if ((i > 0) && (i < (dims0-1)))
out[vec_index] = in[index+iskip] - in[index-iskip];
else if (i == 0)
{
pt2 = pts + 3*index;
out[vec_index] = in[index+iskip] - in[index];
}
else // i == dims0-1
{
pt1 = pts + 3*index;
out[vec_index] = in[index] - in[index-iskip];
}
float diff[3];
diff[0] = pt1[0] - pt2[0];
diff[1] = pt1[1] - pt2[1];
diff[2] = pt1[2] - pt2[2];
float dist = sqrt(diff[0]*diff[0] + diff[1]*diff[1] + diff[2]*diff[2]);
if (dist == 0.)
out[vec_index++] = 0.;
else
out[vec_index++] /= dist;
pt1 = pts + 3*(index+jskip);
pt2 = pts + 3*(index-jskip);
if ((j > 0) && (j < (dims1-1)))
out[vec_index] = in[index+jskip] - in[index-jskip];
else if (j == 0)
{
pt2 = pts + 3*index;
out[vec_index] = in[index+jskip] - in[index];
}
else // j == dims1-1
{
pt1 = pts + 3*index;
out[vec_index] = in[index] - in[index-jskip];
}
diff[0] = pt1[0] - pt2[0];
diff[1] = pt1[1] - pt2[1];
diff[2] = pt1[2] - pt2[2];
dist = sqrt(diff[0]*diff[0] + diff[1]*diff[1] + diff[2]*diff[2]);
if (dist == 0.)
out[vec_index++] = 0.;
else
out[vec_index++] /= dist;
out[vec_index++] = 0.;
}
}
}
else
{
for (k = 0 ; k < dims2 ; k++)
{
for (j = 0 ; j < dims1 ; j++)
{
for (i = 0 ; i < dims0 ; i++)
{
int index = k*kskip + j*jskip + i*iskip;
int vec_index = 3*index;
float *pt1 = pts + 3*(index+iskip);
float *pt2 = pts + 3*(index-iskip);
if ((i > 0) && (i < (dims0-1)))
out[vec_index] = in[index+iskip]-in[index-iskip];
else if (i == 0)
{
pt2 = pts + 3*index;
out[vec_index] = in[index+iskip] - in[index];
}
else // i == dims0-1
{
pt1 = pts + 3*index;
out[vec_index] = in[index] - in[index-iskip];
}
float diff[3];
diff[0] = pt1[0] - pt2[0];
diff[1] = pt1[1] - pt2[1];
diff[2] = pt1[2] - pt2[2];
float dist = sqrt(diff[0]*diff[0] + diff[1]*diff[1] + diff[2]*diff[2]);
if (dist == 0.)
out[vec_index++] = 0.;
else
out[vec_index++] /= dist;
pt1 = pts + 3*(index+jskip);
pt2 = pts + 3*(index-jskip);
if ((j > 0) && (j < (dims1-1)))
out[vec_index] = in[index+jskip] - in[index-jskip];
else if (j == 0)
{
pt2 = pts + 3*index;
out[vec_index] = in[index+jskip] - in[index];
}
else // j == dims1-1
{
pt1 = pts + 3*index;
out[vec_index] = in[index] - in[index-jskip];
}
diff[0] = pt1[0] - pt2[0];
diff[1] = pt1[1] - pt2[1];
diff[2] = pt1[2] - pt2[2];
dist = sqrt(diff[0]*diff[0] + diff[1]*diff[1] + diff[2]*diff[2]);
if (dist == 0.)
out[vec_index++] = 0.;
else
out[vec_index++] /= dist;
pt1 = pts + 3*(index+kskip);
pt2 = pts + 3*(index-kskip);
if ((k > 0) && (k < (dims2-1)))
out[vec_index] = in[index+kskip] - in[index-kskip];
else if (k == 0)
{
pt2 = pts + 3*index;
out[vec_index] = in[index+kskip] - in[index];
}
else // k == dims2-1
{
pt1 = pts + 3*index;
out[vec_index] = in[index] - in[index-kskip];
}
diff[0] = pt1[0] - pt2[0];
diff[1] = pt1[1] - pt2[1];
diff[2] = pt1[2] - pt2[2];
dist = sqrt(diff[0]*diff[0] + diff[1]*diff[1] + diff[2]*diff[2]);
if (dist == 0.)
out[vec_index++] = 0.;
else
out[vec_index++] /= dist;
}
}
}
}
return out_array;
}
// ****************************************************************************
// Method: avtGradientFilter::PerformRestriction
//
// Purpose:
// Request ghost zones.
//
// Programmer: Hank Childs
// Creation: February 13, 2006
//
// ****************************************************************************
avtPipelineSpecification_p
avtGradientFilter::PerformRestriction(avtPipelineSpecification_p in_spec)
{
avtPipelineSpecification_p spec2 =
avtSingleInputExpressionFilter::PerformRestriction(in_spec);
spec2->GetDataSpecification()->SetDesiredGhostDataType(GHOST_ZONE_DATA);
return spec2;
}
......@@ -13,6 +13,7 @@ class vtkDataArray;
class vtkDataSet;
class vtkIdList;
class vtkRectilinearGrid;
class vtkStructuredGrid;
// ****************************************************************************
......@@ -34,6 +35,10 @@ class vtkRectilinearGrid;
// Hank Childs, Fri Mar 4 08:21:04 PST 2005
// Removed data centering conversion modules.
//
// Hank Childs, Mon Feb 13 14:45:18 PST 2006
// Add support for logical gradients. Also add perform restriction, so we
// can request ghost zones.
//
// ****************************************************************************
class EXPRESSION_API avtGradientFilter : public avtSingleInputExpressionFilter
......@@ -42,11 +47,18 @@ class EXPRESSION_API avtGradientFilter : public avtSingleInputExpressionFilter
avtGradientFilter();
virtual ~avtGradientFilter();
void SetDoLogicalGradient(bool b)
{ doLogicalGradients = b; };
virtual const char *GetType(void) { return "avtGradientFilter"; };
virtual const char *GetDescription(void)
{ return "Calculating Gradient of Each Node"; };
{ return "Calculating Gradient"; };
virtual void PreExecute(void);
protected:
bool doLogicalGradients;
bool haveIssuedWarning;
virtual vtkDataArray *DeriveVariable(vtkDataSet *);
virtual int GetVariableDimension() { return 3; }
......@@ -57,6 +69,9 @@ class EXPRESSION_API avtGradientFilter : public avtSingleInputExpressionFilter
float EvaluateValue(float, float, float, vtkDataSet *,
vtkDataArray *,vtkIdList *,bool &);
vtkDataArray *RectilinearGradient(vtkRectilinearGrid *);
vtkDataArray *LogicalGradient(vtkStructuredGrid *);
virtual avtPipelineSpecification_p
PerformRestriction(avtPipelineSpecification_p);
};
......
......@@ -431,6 +431,9 @@ avtVectorExpr::CreateFilters(ExprPipelineState *state)
// Hank Childs, Sat Jan 21 14:43:45 PST 2006
// Added symm_eval_transform.
//
// Hank Childs, Tue Feb 14 14:03:47 PST 2006
// Added logical gradient.
//
// ****************************************************************************
void
avtFunctionExpr::CreateFilters(ExprPipelineState *state)
......@@ -523,6 +526,12 @@ avtFunctionExpr::CreateFilters(ExprPipelineState *state)
f = new avtMeshCoordinateFilter();
else if (functionName == "procid")
f = new avtProcessorIdFilter();
else if (functionName == "ijk_gradient" || functionName == "ij_gradient")
{
avtGradientFilter *g = new avtGradientFilter();
g->SetDoLogicalGradient(true);
f = g;
}
else if (functionName == "gradient")
f = new avtGradientFilter();
else if (functionName == "curl")
......
......@@ -211,6 +211,9 @@ avtWholeImageCompositerWithZ::~avtWholeImageCompositerWithZ()
// Jeremy Meredith, October 20, 2004
// Allowed for the use of an allreduce instead of a simple reduce.
//
// Hank Childs, Mon Feb 6 14:55:39 PST 2006
// Fix memory leak ['6829].
//
// ****************************************************************************
void
......@@ -265,7 +268,6 @@ avtWholeImageCompositerWithZ::Execute(void)
iorgb = zeroImageRep.GetRGBBuffer();
}
if (mpiRoot >= 0)
{
// only root allocates output AVT image (for a non-allreduce)
......@@ -291,7 +293,7 @@ avtWholeImageCompositerWithZ::Execute(void)
{
if (shouldOutputZBuffer)
{
avtImageRepresentation theOutput(mergedGlobalImage,rioz);
avtImageRepresentation theOutput(mergedGlobalImage,rioz,true);
SetOutputImage(theOutput);
}
else
......@@ -314,7 +316,7 @@ avtWholeImageCompositerWithZ::Execute(void)
{
if (shouldOutputZBuffer)
{
avtImageRepresentation theOutput(mergedLocalImage,ioz);
avtImageRepresentation theOutput(mergedLocalImage,ioz,true);
SetOutputImage(theOutput);
}
else
......@@ -327,16 +329,8 @@ avtWholeImageCompositerWithZ::Execute(void)
}
else
{
if (shouldOutputZBuffer)
{
avtImageRepresentation theOutput(zeroImageRep.GetImageVTK(),zeroImageRep.GetZBuffer());
SetOutputImage(theOutput);
}
else
{
avtImageRepresentation theOutput(zeroImageRep.GetImageVTK());
SetOutputImage(theOutput);
}
avtImageRepresentation theOutput(zeroImageRep);
SetOutputImage(theOutput);
}
}
}
......
......@@ -68,9 +68,15 @@ avtImageRepresentation::avtImageRepresentation(vtkImageData *d)
// Programmer: Hank Childs
// Creation: February 13, 2001
//
// Modifications:
//
// Hank Childs, Mon Feb 6 14:59:43 PST 2006
// Allow the z-buffer to be directly acquired instead of copied.
//
// ****************************************************************************
avtImageRepresentation::avtImageRepresentation(vtkImageData *d, float *z)
avtImageRepresentation::avtImageRepresentation(vtkImageData *d, float *z,
bool iNowOwn)
{
Initialize();
int width = 0;
......@@ -86,8 +92,13 @@ avtImageRepresentation::avtImageRepresentation(vtkImageData *d, float *z)
if (z != NULL)
{
zbuffer = new float[width * height];
memcpy(zbuffer, z, sizeof(float) * width * height);
if (iNowOwn)
zbuffer = z;
else
{
zbuffer = new float[width * height];