Commit 8fa7acc0 authored by hrchilds's avatar hrchilds
Browse files

Update from May 16, 2006

git-svn-id: http://visit.ilight.com/svn/visit/trunk/src@667 18c085ea-50e0-402c-830e-de6fd14e8384
parent 77d957c9
......@@ -83,6 +83,9 @@ avtDDFConstructor::~avtDDFConstructor()
// Hank Childs, Sat Feb 25 15:01:56 PST 2006
// Add support for time.
//
// Hank Childs, Mon May 15 14:25:28 PDT 2006
// Fix memory leak.
//
// ****************************************************************************
avtDDF *
......@@ -346,6 +349,7 @@ avtDDFConstructor::ConstructDDF(ConstructDDFAttributes *atts,
}
}
delete [] arr;
delete [] args;
delete [] leaves;
delete [] isNodal;
}
......
// ************************************************************************* //
// avtCurvatureExpression.C //
// ************************************************************************* //
#include <avtCurvatureExpression.h>
#include <vtkCellData.h>
#include <vtkCurvatures.h>
#include <vtkDataSet.h>
#include <vtkFloatArray.h>
#include <vtkPointData.h>
#include <vtkPolyData.h>
#include <vtkRectilinearGrid.h>
#include <ExpressionException.h>
// ****************************************************************************
// Method: avtCurvatureExpression constructor
//
// Purpose:
// Defines the constructor. Note: this should not be inlined in the
// header because it causes problems for certain compilers.
//
// Programmer: Hank Childs
// Creation: May 11, 2006
//
// ****************************************************************************
avtCurvatureExpression::avtCurvatureExpression()
{
doGauss = false;
}
// ****************************************************************************
// Method: avtCurvatureExpression destructor
//
// Purpose:
// Defines the destructor. Note: this should not be inlined in the header
// because it causes problems for certain compilers.
//
// Programmer: Hank Childs
// Creation: May 11, 2006
//
// ****************************************************************************
avtCurvatureExpression::~avtCurvatureExpression()
{
;
}
// ****************************************************************************
// Method: avtCurvatureExpression::DeriveVariable
//
// Purpose:
// Calculates the localized compactness at every point.
//
// Arguments:
// inDS The input dataset.
//
// Returns: The derived variable. The calling class must free this
// memory.
//
// Programmer: Hank Childs
// Creation: May 11, 2006
//
// ****************************************************************************
vtkDataArray *
avtCurvatureExpression::DeriveVariable(vtkDataSet *in_ds)
{
if (in_ds->GetDataObjectType() != VTK_POLY_DATA)
{
EXCEPTION1(ExpressionException, "The curvature expression "
"can only be calculated on surfaces. VisIt tries to "
"evaluate expressions as soon as they are read from a "
"database. For curvature, the expression is typically "
"desired after some operators have been applied. You "
"can defer the evaluation of the curvature expression "
"using the DeferExpression operator."
" The defer expression operator is available through "
"the plugin manager located under the Options menu");
}
vtkPolyData *pd = (vtkPolyData *) in_ds;
vtkCurvatures *curvatures = vtkCurvatures::New();
curvatures->SetCurvatureType((doGauss ? VTK_CURVATURE_GAUSS
: VTK_CURVATURE_MEAN));
curvatures->SetInput(pd);
curvatures->Update();
vtkPolyData *out = curvatures->GetOutput();
vtkDataArray *curvature = out->GetPointData()->GetArray(
doGauss ? "Gauss_Curvature" : "Mean_Curvature");
bool shouldDelete = false;
vtkFloatArray *flt_curvature = vtkFloatArray::New();
int npts = pd->GetNumberOfPoints();
flt_curvature->SetNumberOfTuples(npts);
if (curvature == NULL)
{
for (int i = 0 ; i < npts ; i++)
flt_curvature->SetTuple1(i, 0.);
}
else
{
for (int i = 0 ; i < npts ; i++)
flt_curvature->SetTuple1(i, curvature->GetTuple1(i));
}
curvatures->Delete();
// Calling function will clean up memory.
return flt_curvature;
}
// ************************************************************************* //
// avtCurvatureExpression.h //
// ************************************************************************* //
#ifndef AVT_CURVATURE_EXPRESSION_H
#define AVT_CURVATURE_EXPRESSION_H
#include <avtSingleInputExpressionFilter.h>
// ****************************************************************************
// Class: avtCurvatureExpression
//
// Purpose:
// Calculates the mean or gaussian curvature using a VTK filter. The
// input must be a surface.
//
// Programmer: Hank Childs
// Creation: May 11, 2006
//
// ****************************************************************************
class EXPRESSION_API avtCurvatureExpression
: public avtSingleInputExpressionFilter
{
public:
avtCurvatureExpression();
virtual ~avtCurvatureExpression();
virtual const char *GetType(void)
{ return "avtCurvatureExpression"; };
virtual const char *GetDescription(void)
{return "Calculating curvature"; };
void DoGaussCurvature(bool val) { doGauss = val; };
protected:
bool doGauss;
virtual vtkDataArray *DeriveVariable(vtkDataSet *);
virtual bool IsPointVariable(void) { return true; };
virtual int GetVariableDimension(void) { return 1; };
};
#endif
......@@ -141,6 +141,9 @@
# Hank Childs, Sat Apr 29 14:40:47 PDT 2006
# Added avtLocalizedCompactnessExpression.
#
# Hank Childs, Thu May 11 12:14:51 PDT 2006
# Added avtCurvatureExpression.
#
##############################################################################
@SET_MAKE@
......@@ -187,6 +190,7 @@ CMFE_src= \
Derivations_src=\
Derivations/avtCurvatureExpression.C \
Derivations/avtEffectiveTensorFilter.C \
Derivations/avtLocalizedCompactnessExpression.C \
Derivations/avtTensorMaximumShearFilter.C \
......
......@@ -29,6 +29,7 @@
#include <avtPrincipalDeviatoricTensorFilter.h>
#include <avtPrincipalTensorFilter.h>
#include <avtEffectiveTensorFilter.h>
#include <avtCurvatureExpression.h>
#include <avtGradientFilter.h>
#include <avtCurlFilter.h>
#include <avtDivergenceFilter.h>
......@@ -450,6 +451,9 @@ avtVectorExpr::CreateFilters(ExprPipelineState *state)
// Hank Childs, Sat Apr 29 14:40:47 PDT 2006
// Added localized compactness expression.
//
// Hank Childs, Thu May 11 12:14:51 PDT 2006
// Added curvature.
//
// ****************************************************************************
void
avtFunctionExpr::CreateFilters(ExprPipelineState *state)
......@@ -544,6 +548,18 @@ avtFunctionExpr::CreateFilters(ExprPipelineState *state)
f = new avtMeshCoordinateFilter();
else if (functionName == "procid")
f = new avtProcessorIdFilter();
else if (functionName == "mean_curvature")
{
avtCurvatureExpression *c = new avtCurvatureExpression;
c->DoGaussCurvature(false);
f = c;
}
else if (functionName == "gauss_curvature")
{
avtCurvatureExpression *c = new avtCurvatureExpression;
c->DoGaussCurvature(true);
f = c;
}
else if (functionName == "ijk_gradient" || functionName == "ij_gradient")
{
avtGradientFilter *g = new avtGradientFilter();
......
......@@ -122,6 +122,9 @@
# Hank Childs, Sat Apr 29 14:40:47 PDT 2006
# Added localized and elliptical compactness factor queries.
#
# Hank Childs, Thu May 11 13:21:18 PDT 2006
# Added average mean curvature query.
#
##############################################################################
@SET_MAKE@
......@@ -162,6 +165,7 @@ Queries_src= \
Queries/avtActualDataNumNodesQuery.C \
Queries/avtActualDataNumZonesQuery.C \
Queries/avtAreaBetweenCurvesQuery.C \
Queries/avtAverageMeanCurvatureQuery.C \
Queries/avtBestFitLineQuery.C \
Queries/avtCentroidQuery.C \
Queries/avtCompactnessQuery.C \
......
......@@ -7,6 +7,7 @@
#include <avtActualDataNumNodesQuery.h>
#include <avtActualDataNumZonesQuery.h>
#include <avtAreaBetweenCurvesQuery.h>
#include <avtAverageMeanCurvatureQuery.h>
#include <avtBestFitLineQuery.h>
#include <avtCentroidQuery.h>
#include <avtCompactnessQuery.h>
......@@ -164,6 +165,9 @@ avtQueryFactory::Instance()
// Hank Childs, Sat Apr 29 14:40:47 PDT 2006
// Added localized and elliptical compactness factor queries.
//
// Hank Childs, Thu May 11 13:21:18 PDT 2006
// Added average mean curvature.
//
// ****************************************************************************
......@@ -354,6 +358,10 @@ avtQueryFactory::CreateQuery(const QueryAttributes *qa)
{
query = new avtOriginalDataSpatialExtentsQuery();
}
else if (qname == "Average Mean Curvature")
{
query = new avtAverageMeanCurvatureQuery();
}
return query;
}
......
// ************************************************************************* //
// avtAverageMeanCurvatureQuery.C //
// ************************************************************************* //
#include <avtAverageMeanCurvatureQuery.h>
#include <avtCurvatureExpression.h>
#include <DebugStream.h>
#include <NonQueryableInputException.h>
using std::string;
// ****************************************************************************
// Method: avtAverageMeanCurvatureQuery constructor
//
// Programmer: Hank Childs
// Creation: May 11, 2006
//
// ****************************************************************************
avtAverageMeanCurvatureQuery::avtAverageMeanCurvatureQuery()
: avtWeightedVariableSummationQuery()
{
curvature = new avtCurvatureExpression;
curvature->DoGaussCurvature(false);
curvature->SetOutputVariableName("curvature");
}
// ****************************************************************************
// Method: avtAverageMeanCurvatureQuery destructor
//
// Programmer: Hank Childs
// Creation: May 11, 2006
//
// ****************************************************************************
avtAverageMeanCurvatureQuery::~avtAverageMeanCurvatureQuery()
{
delete curvature;
}
// ****************************************************************************
// Method: avtAverageMeanCurvatureQuery::CreateVariable
//
// Purpose:
// Creates the variable for the summation.
//
// Programmer: Hank Childs
// Creation: May 11, 2006
//
// ****************************************************************************
avtDataObject_p
avtAverageMeanCurvatureQuery::CreateVariable(avtDataObject_p inData)
{
curvature->SetInput(inData);
return curvature->GetOutput();
}
// ****************************************************************************
// Method: avtAverageMeanCurvatureQuery::VerifyInput
//
// Purpose:
// Make sure we are operating on a surface.
//
// Programmer: Hank Childs
// Creation: May 11, 2006
//
// ****************************************************************************
void
avtAverageMeanCurvatureQuery::VerifyInput(void)
{
//
// We want to do this in addition to what the base class does, so call the
// base class' version of this method as well.
//
avtWeightedVariableSummationQuery::VerifyInput();
if (GetInput()->GetInfo().GetAttributes().GetTopologicalDimension() != 2)
{
EXCEPTION1(NonQueryableInputException, "The average mean curvature "
"query can only operate on surfaces.");
}
}
// ************************************************************************* //
// avtAverageMeanCurvatureQuery.h //
// ************************************************************************* //
#ifndef AVT_AVERAGE_MEAN_CURVATURE_QUERY_H
#define AVT_AVERAGE_MEAN_CURVATURE_QUERY_H
#include <query_exports.h>
#include <avtWeightedVariableSummationQuery.h>
class avtCurvatureExpression;
// ****************************************************************************
// Class: avtAverageMeanCurvatureQuery
//
// Purpose:
// A query that will calculate the mean curvature and then find the
// average of that quantity.
//
// Programmer: Hank Childs
// Creation: May 11, 2006
//
// ****************************************************************************
class QUERY_API avtAverageMeanCurvatureQuery : public avtWeightedVariableSummationQuery
{
public:
avtAverageMeanCurvatureQuery();
virtual ~avtAverageMeanCurvatureQuery();
virtual const char *GetType(void)
{ return "avtAverageMeanCurvatureQuery"; };
protected:
avtCurvatureExpression *curvature;
virtual avtDataObject_p CreateVariable(avtDataObject_p d);
virtual std::string GetVarname(std::string &s)
{ return "curvature"; };
virtual void VerifyInput(void);
virtual bool CalculateAverage(void) { return true; };
};
#endif
......@@ -89,6 +89,13 @@ avtEllipticalCompactnessFactorQuery::PreExecute(void)
for (int i = 0 ; i < 3 ; i++)
centroid[i] = 0.;
total_volume = 0.;
bounds[0] = +FLT_MAX;
bounds[1] = -FLT_MAX;
bounds[2] = +FLT_MAX;
bounds[3] = -FLT_MAX;
bounds[4] = +FLT_MAX;
bounds[5] = -FLT_MAX;
}
......@@ -119,19 +126,70 @@ avtEllipticalCompactnessFactorQuery::MidExecute(void)
centroid[1] = C_tmp[1];
centroid[2] = C_tmp[2];
volume_inside = 0.;
radius = pow(total_volume*0.75/M_PI, 0.3333333);
for (int i = 0 ; i < numGuesses ; i++)
volume_inside[i] = 0.;
//
// Now we have to choose the lengths of the major and minor axes, which
// is hard. One choice would be to let them vary from 0 to infinity,
// but we only have a few steps in between, so the steps would be too
// coarse. So we need to narrow the range. Using the bounds of the
// data seems to be a good place to start. But, we don't want to set
// the length of the axes to be based on these axes, because this
// won't cover the things at the extrema. Restated, if X spans (-a,a)
// and Y spans (0,b) (cylindrical), then choosing a and b as axes lengths
// would give us an ellipsoid with volume 4/3pi(a)(b)(b), but the cylinder
// that volume covers would be 2pi(a)(b)(b). So let's lengthen out a and
// b by (50%)^1/3. Then the ellipsoid would have a matching volume.
//
if (is2D)
{
float Amax = ((bounds[1] - bounds[0]) / 2.)*1.5;
float Bmax = (bounds[3])*0.5;
float Bmin = sqrt(total_volume * 0.75 / (M_PI*Amax));
for (int i = 0 ; i < numGuesses ; i++)
{
y_radius[i] = Bmin + (Bmax-Bmin)*((float)i/(float)numGuesses);
z_radius[i] = y_radius[i];
x_radius[i] = total_volume*0.75 / (M_PI*y_radius[i]*z_radius[i]);
}
}
else
{
float Amax = ((bounds[1] - bounds[0]) / 2.)*1.5;
float Bmax = ((bounds[3] - bounds[2]) / 2.)*1.5;
float Cmax = ((bounds[5] - bounds[4]) / 2.)*1.5;
float Amin = sqrt(total_volume * 0.75 / (M_PI*Bmax*Cmax));
float Bmin = sqrt(total_volume * 0.75 / (M_PI*Amax*Cmax));
float Cmin = sqrt(total_volume * 0.75 / (M_PI*Amax*Bmax));
// Get integer square root.
int dims = (int) ceil(sqrt(numGuesses-0.1));
for (int i = 0 ; i < dims ; i++)
{
float A = Amin + (Amax-Amin)*((float)i/(float)dims);
for (int j = 0 ; j < dims ; j++)
{
float B = Bmin + (Bmax-Bmin)*((float)j/(float)dims);
float C = total_volume*0.75 / (M_PI*A*B);
int index = i*dims+j;
x_radius[index] = A;
y_radius[index] = B;
z_radius[index] = C;
}
}
}
if (is2D)
{
sphere_center[0] = centroid[0];
sphere_center[1] = 0.;
sphere_center[2] = 0.;
ellipse_center[0] = centroid[0];
ellipse_center[1] = 0.;
ellipse_center[2] = 0.;
}
else
{
sphere_center[0] = centroid[0];
sphere_center[1] = centroid[1];
sphere_center[2] = centroid[2];
ellipse_center[0] = centroid[0];
ellipse_center[1] = centroid[1];
ellipse_center[2] = centroid[2];
}
}
......@@ -150,10 +208,17 @@ avtEllipticalCompactnessFactorQuery::MidExecute(void)
void
avtEllipticalCompactnessFactorQuery::PostExecute(void)
{
SumDoubleAcrossAllProcessors(volume_inside);
if (total_volume == 0.)
double vi_tmp[numGuesses];
SumDoubleArrayAcrossAllProcessors(volume_inside, vi_tmp, numGuesses);
int biggest = 0;
double biggestVal = vi_tmp[0];
for (int i = 1 ; i < numGuesses ; i++)
{
total_volume = 1;
if (vi_tmp[i] > biggestVal)
{
biggestVal = vi_tmp[i];
biggest = i;
}
}
//
......@@ -163,13 +228,13 @@ avtEllipticalCompactnessFactorQuery::PostExecute(void)
//
char msg[4096];
SNPRINTF(msg, 4096, "Elliptical Compactness Factor = %f. Using centroid "
"for sphere origin. Centroid used was (%f, %f, %f)."
" Radius was %f",
volume_inside / total_volume,
sphere_center[0], sphere_center[1], sphere_center[2],
radius);
"for ellipse origin. Centroid used was (%f, %f, %f)."
" Best fitting axes were %f,%f,%f.",
biggestVal / total_volume,
ellipse_center[0], ellipse_center[1], ellipse_center[2],
x_radius[biggest], y_radius[biggest], z_radius[biggest]);
SetResultMessage(msg);
SetResultValue(volume_inside / total_volume);
SetResultValue(biggestVal / total_volume);
}
......@@ -207,6 +272,12 @@ avtEllipticalCompactnessFactorQuery::Execute1(vtkDataSet *ds, const int dom)
centroid[0] += volume*center[0];
centroid[1] += volume*center[1];
centroid[2] += volume*center[2];
bounds[0] = (center[0] < bounds[0] ? center[0] : bounds[0]);
bounds[1] = (center[0] > bounds[1] ? center[0] : bounds[1]);
bounds[2] = (center[1] < bounds[2] ? center[1] : bounds[2]);
bounds[3] = (center[1] > bounds[3] ? center[1] : bounds[3]);
bounds[4] = (center[2] < bounds[4] ? center[2] : bounds[4]);
bounds[5] = (center[2] > bounds[5] ? center[2] : bounds[5]);
total_volume += volume;
}
}
......@@ -234,7 +305,7 @@ avtEllipticalCompactnessFactorQuery::Execute2(vtkDataSet *ds, const int dom)
{
EXCEPTION0(ImproperUseException);
}
float rad_squared = radius*radius;
float d[3];
for (int i = 0 ; i < nCells ; i++)
{
if (ghosts != NULL && ghosts->GetTuple1(i) != 0.)
......@@ -242,14 +313,22 @@ avtEllipticalCompactnessFactorQuery::Execute2(vtkDataSet *ds, const int dom)
vtkCell *cell = ds->GetCell(i);
double center[3];
vtkVisItUtility::GetCellCenter(cell, center);
float dist = (center[0]-sphere_center[0])*(center[0]-sphere_center[0])
+ (center[1]-sphere_center[1])*(center[1]-sphere_center[1])
+ (center[2]-sphere_center[2])*(center[2]-sphere_center[2]);
if (dist > rad_squared)
continue; // Not inside.
d[0] = center[0] - ellipse_center[0];