Commit e6c9eceb authored by bonnell's avatar bonnell

Double precision support.

git-svn-id: http://visit.ilight.com/svn/visit/trunk/src@19132 18c085ea-50e0-402c-830e-de6fd14e8384
parent 89d36f66
......@@ -251,6 +251,68 @@ private:
bool own;
};
// ****************************************************************************
// Class: vtkDirectAccessor
//
// Purpose:
// Access vtkDataArray using direct memory accesses.
//
// *Copied from avtDirectAccessor *
//
// Programmer: Kathleen Biagas
// Creation: September 4, 2012
//
// Modifications:
//
// ****************************************************************************
template <typename T>
class vtkDirectAccessor
{
public:
vtkDirectAccessor(vtkDataArray *da) : arr(da)
{
N = arr->GetNumberOfComponents();
ptr = start = (T *)arr->GetVoidPointer(0);
end = ptr + (N * arr->GetNumberOfTuples());
}
~vtkDirectAccessor() { }
int GetNumberOfTuples() const { return arr->GetNumberOfTuples(); }
int GetNumberOfComponents() const { return N; }
void SetTuple1(T val) { *ptr = val; }
void SetTuple1(vtkIdType i, T val) { start[i*N] = val; }
void SetTuple(vtkIdType i, const T *val) { memcpy(start + i*N, val, sizeof(T)*N); }
void SetTuple(const T *val) { memcpy(ptr, val, sizeof(T)*N); }
T GetTuple1() const { return *ptr; }
T GetTuple1(vtkIdType i) const { return start[i*N]; }
T *GetTuple(vtkIdType i) const { return start + i*N; }
T *GetTuple() const { return ptr; }
void SetComponent(vtkIdType tuple, int comp, T val) { start[tuple * N + comp] = val; }
void SetComponent(int comp, T val) { ptr[comp] = val; }
T GetComponent(vtkIdType tuple, int comp) const { return start[tuple * N + comp]; }
T GetComponent(int comp) const { return ptr[comp]; }
void InitTraversal() { ptr = start; }
bool Iterating() const { return ptr < end; }
void operator ++() { ptr += N; }
void operator ++(int i) { ptr += N; }
private:
vtkDataArray *arr;
T *start;
T *ptr;
T *end;
int N;
};
// ****************************************************************************
// Class: vtkGeneralAccessor
//
......@@ -263,7 +325,9 @@ private:
// Creation: Thu Mar 29 13:11:41 PDT 2012
//
// Modifications:
//
// Kathleen Biagas, Thu Sep 6 11:16:12 MST 2012
// Added more methods (from avtTupleAccessor).
//
// ****************************************************************************
class vtkGeneralAccessor
......@@ -271,19 +335,40 @@ class vtkGeneralAccessor
public:
vtkGeneralAccessor(vtkDataArray *da) : arr(da)
{
index = 0;
size = arr->GetNumberOfTuples();
}
~vtkGeneralAccessor() { }
int GetNumberOfTuples() const { return arr->GetNumberOfTuples(); }
int GetNumberOfComponents() const { return arr->GetNumberOfComponents(); }
inline double GetTuple1(vtkIdType id) const
{
return arr->GetTuple1(id);
}
void SetTuple1(double val) { arr->SetTuple1(index, val); }
void SetTuple1(vtkIdType i, double val) { arr->SetTuple1(i, val); }
void SetTuple(vtkIdType i, const double *val) { arr->SetTuple(i, val); }
void SetTuple(const double *val) { arr->SetTuple(index, val); }
double GetTuple1() { return arr->GetTuple1(index); }
double GetTuple1(vtkIdType i) { return arr->GetTuple1(i); }
double *GetTuple(vtkIdType i) { return arr->GetTuple(i); }
double *GetTuple() { return arr->GetTuple(index); }
void SetComponent(vtkIdType tuple, int comp, double val) { arr->SetComponent(tuple, comp, val); }
void SetComponent(int comp, double val) { arr->SetComponent(index, comp, val); }
double GetComponent(vtkIdType tuple, int comp) const { return arr->GetComponent(tuple, comp); }
double GetComponent(int comp) const { return arr->GetComponent(index, comp); }
void InitTraversal() { index = 0; }
bool Iterating() const { return index < size; }
void operator ++() { index ++; }
void operator ++(int) { index ++; }
inline void SetTuple1(vtkIdType id, double val)
{
arr->SetTuple1(id, val);
}
private:
vtkDataArray *arr;
vtkIdType index;
vtkIdType size;
};
#endif
......@@ -199,7 +199,7 @@ bool
vtkConnectedTubeFilter::PointSequenceList::Build(vtkPoints *points,
vtkCellArray *lines)
{
pts = (float*)(points->GetVoidPointer(0));
pts = points;
len = points->GetNumberOfPoints();
numneighbors = new int[len];
connectivity[0] = new int[len];
......@@ -322,9 +322,11 @@ vtkConnectedTubeFilter::PointSequenceList::GetNextSequence(PointSequence &seq)
// we must skip any sequential identical points:
// 1) they are useless, and 2) they mess up calculations
if (pts[previous*3+0] != pts[current*3+0] ||
pts[previous*3+1] != pts[current*3+1] ||
pts[previous*3+2] != pts[current*3+2])
double *prePt = pts->GetPoint(previous);
double *curPt = pts->GetPoint(current);
if (prePt[0] != curPt[0] ||
prePt[1] != curPt[1] ||
prePt[2] != curPt[2])
{
seq.Add(current, cellindex[current]);
}
......@@ -467,6 +469,9 @@ bool vtkConnectedTubeFilter::BuildConnectivityArrays()
// degenerate polys between 360 and 0 degrees) by partitioning into
// npts-1 sides, instead of npts sides. Fixed that.
//
// Kathleen Biagas, Thu Sep 6 11:15:29 MST 2012
// Preserve coordinate data type.
//
// ****************************************************************************
void vtkConnectedTubeFilter::Execute()
{
......@@ -497,17 +502,15 @@ void vtkConnectedTubeFilter::Execute()
return;
}
float *pts = (float*)(inPts->GetVoidPointer(0));
// Set up the output arrays
int maxNewCells = numCells * (NumberOfSides + 2);
int maxNewPoints = numCells * NumberOfSides * 2;
vtkPolyData *output = this->GetOutput();
vtkPoints *newPts = vtkPoints::New();
vtkPoints *newPts = vtkPoints::New(inPts->GetDataType());
newPts->Allocate(maxNewPoints);
vtkCellArray *newCells = vtkCellArray::New();
newCells->Allocate(maxNewCells, 4*maxNewCells + 2*NumberOfSides);
vtkFloatArray *newNormals = NULL;
vtkDataArray *newNormals = NULL;
vtkPointData *newPD = NULL;
newPD = output->GetPointData();
......@@ -520,7 +523,7 @@ void vtkConnectedTubeFilter::Execute()
if (CreateNormals)
{
newNormals = vtkFloatArray::New();
newNormals = inPts->GetData()->NewInstance();
newNormals->SetNumberOfComponents(3);
newNormals->SetName("Normals");
newPD->SetNormals(newNormals);
......@@ -547,22 +550,19 @@ void vtkConnectedTubeFilter::Execute()
int ix2 = (lastPoint ? seq.index[i] : seq.index[i+1]);
// Use a centered difference approximation for direction
float dir[3] = {pts[ix2*3+0] - pts[ix1*3+0],
pts[ix2*3+1] - pts[ix1*3+1],
pts[ix2*3+2] - pts[ix1*3+2]};
double dir[3];
vtkMath::Subtract(inPts->GetPoint(ix2),inPts->GetPoint(ix1), dir);
// If our centered difference was zero, do a forward
// difference instead. We ensured no sequential points
// are identical, so this can't fail.
if (dir[0]==0 && dir[1]==0 && dir[2]==0)
{
dir[0] = pts[ix2*3+0] - pts[ix*3+0];
dir[1] = pts[ix2*3+1] - pts[ix*3+1];
dir[2] = pts[ix2*3+2] - pts[ix*3+2];
vtkMath::Subtract(inPts->GetPoint(ix2),inPts->GetPoint(ix),dir);
}
// Get a couple vectors orthogonal to our direction
float v1[3], v2[3];
double v1[3], v2[3];
vtkMath::Perpendiculars(dir, v1,v2, 0.0);
vtkMath::Normalize(v1);
vtkMath::Normalize(v2);
......@@ -570,18 +570,19 @@ void vtkConnectedTubeFilter::Execute()
// Hang on to the first point index we create; we need it
// to create the cells
vtkIdType firstIndex = newPts->GetNumberOfPoints();
double *pt = inPts->GetPoint(ix);
for (int j = 0 ; j < NumberOfSides ; j++)
{
float q = (j * 2. * vtkMath::Pi()) / float(NumberOfSides);
float sq = sin(q);
float cq = cos(q);
float normal[3] = { v1[0]*cq + v2[0]*sq,
double q = (j * 2. * vtkMath::Pi()) / double(NumberOfSides);
double sq = sin(q);
double cq = cos(q);
double normal[3] = { v1[0]*cq + v2[0]*sq,
v1[1]*cq + v2[1]*sq,
v1[2]*cq + v2[2]*sq};
float x = pts[ix*3+0] + Radius * normal[0];
float y = pts[ix*3+1] + Radius * normal[1];
float z = pts[ix*3+2] + Radius * normal[2];
double x = pt[0] + Radius * normal[0];
double y = pt[1] + Radius * normal[1];
double z = pt[2] + Radius * normal[2];
vtkIdType id = newPts->InsertNextPoint(x,y,z);
if (CreateNormals)
newNormals->InsertNextTuple(normal);
......
......@@ -132,7 +132,7 @@ class VISIT_VTK_API vtkConnectedTubeFilter :
int *numneighbors;
int *connectivity[2];
int *cellindex;
const float *pts;
vtkPoints *pts;
// traversal variables
bool *visited;
......
......@@ -38,6 +38,7 @@
#include <vtkRectilinearGridFacelistFilter.h>
#include <vtkAccessors.h>
#include <vtkCellArray.h>
#include <vtkCellData.h>
#include <vtkObjectFactory.h>
......@@ -45,6 +46,7 @@
#include <vtkPolyData.h>
#include <vtkRectilinearGrid.h>
#include <vtkUnsignedCharArray.h>
#include <vtkVisItUtility.h>
#include <ImproperUseException.h>
......@@ -167,6 +169,110 @@ SpecializedIndexer::SpecializedIndexer(int x, int y, int z)
cellZBase = 1;
}
// ****************************************************************************
// Kathleen Biagas, Thu Sep 6 11:10:54 MST 2012
// Templated method to process faces, preserves coordinate data type.
//
// ****************************************************************************
template <class Accessor> inline void
vtkRectilinearGridFacelistFilter_ProcessFaces(int nX, int nY, int nZ,
vtkPointData *outPointData, vtkPointData *inPointData,
Accessor x, Accessor y, Accessor z, Accessor p)
{
int i, j, pointId = 0;
p.InitTraversal();
// Front face.
for (i = 0 ; i < nX ; ++i)
{
for (j = 0 ; j < nY ; ++j)
{
p.SetComponent(0, x.GetComponent(i));
p.SetComponent(1, y.GetComponent(j));
p.SetComponent(2, z.GetComponent(0));
++p;
outPointData->CopyData(inPointData, j*nX + i, pointId++);
}
}
// Back face
if (nZ > 1)
{
for (i = 0 ; i < nX ; ++i)
{
for (j = 0 ; j < nY ; ++j)
{
p.SetComponent(0, x.GetComponent(i));
p.SetComponent(1, y.GetComponent(j));
p.SetComponent(2, z.GetComponent(nZ-1));
++p;
outPointData->CopyData(inPointData, (nZ-1)*nX*nY + j*nX + i,pointId++);
}
}
}
// Bottom face
for (i = 0 ; i < nX ; ++i)
{
for (j = 1 ; j < nZ-1 ; ++j)
{
p.SetComponent(0, x.GetComponent(i));
p.SetComponent(1, y.GetComponent(0));
p.SetComponent(2, z.GetComponent(j));
++p;
outPointData->CopyData(inPointData, j*nX*nY + i, pointId++);
}
}
// Top face
if (nY > 1)
{
for (i = 0 ; i < nX ; ++i)
{
for (j = 1 ; j < nZ-1 ; ++j)
{
p.SetComponent(0, x.GetComponent(i));
p.SetComponent(1, y.GetComponent(nY-1));
p.SetComponent(2, z.GetComponent(j));
++p;
outPointData->CopyData(inPointData, j*nX*nY + (nY-1)*nX + i,pointId++);
}
}
}
// Left face
for (i = 1 ; i < nY-1 ; ++i)
{
for (j = 1 ; j < nZ-1 ; ++j)
{
p.SetComponent(0, x.GetComponent(0));
p.SetComponent(1, y.GetComponent(i));
p.SetComponent(2, z.GetComponent(j));
++p;
outPointData->CopyData(inPointData, j*nX*nY + i*nX, pointId++);
}
}
// Right face
if (nX > 1)
{
for (i = 1 ; i < nY-1 ; i++)
{
for (j = 1 ; j < nZ-1 ; j++)
{
p.SetComponent(0, x.GetComponent(nX-1));
p.SetComponent(1, y.GetComponent(i));
p.SetComponent(2, z.GetComponent(j));
++p;
outPointData->CopyData(inPointData, j*nX*nY + i*nX + nX-1, pointId++);
}
}
}
}
// ****************************************************************************
//
// Hank Childs, Thu Aug 15 21:13:43 PDT 2002
......@@ -208,6 +314,9 @@ SpecializedIndexer::SpecializedIndexer(int x, int y, int z)
// Brad Whitlock, Tue Apr 25 15:13:20 PST 2006
// Pass the field data after the shallow copy or it does not get passed.
//
// Kathleen Biagas, Thu Sep 6 11:12:43 MST 2012
// Use new templated method to process faces, to preserve coordinate type.
//
// ****************************************************************************
void vtkRectilinearGridFacelistFilter::Execute()
......@@ -253,32 +362,24 @@ void vtkRectilinearGridFacelistFilter::Execute()
//
vtkDataArray *xc = input->GetXCoordinates();
int nX = xc->GetNumberOfTuples();
float *x = new float[nX];
for (i = 0 ; i < nX ; i++)
{
x[i] = xc->GetTuple1(i);
}
int tX = xc->GetDataType();
vtkDataArray *yc = input->GetYCoordinates();
int nY = yc->GetNumberOfTuples();
float *y = new float[nY];
for (i = 0 ; i < nY ; i++)
{
y[i] = yc->GetTuple1(i);
}
int tY = yc->GetDataType();
vtkDataArray *zc = input->GetZCoordinates();
int nZ = zc->GetNumberOfTuples();
float *z = new float[nZ];
for (i = 0 ; i < nZ ; i++)
{
z[i] = zc->GetTuple1(i);
}
int tZ = zc->GetDataType();
bool same = (tX == tY && tY == tZ);
int type = (same ? tX : (tX == VTK_DOUBLE ? tX : (tY == VTK_DOUBLE ? tY :
(tZ == VTK_DOUBLE ? tZ : VTK_FLOAT))));
//
// Now create the points. Do this so that the total number of points is
// minimal -- this requires sharing points along edges and corners and leads
// to a nasty indexing scheme. Also be wary of 2D issues.
//
vtkPoints *pts = vtkPoints::New();
vtkPoints *pts = vtkPoints::New(type);
int npts = 0;
if (nX <= 1)
{
......@@ -297,103 +398,40 @@ void vtkRectilinearGridFacelistFilter::Execute()
npts = 2*nX*nY + 2*(nZ-2)*nX + 2*(nZ-2)*(nY-2);
}
pts->SetNumberOfPoints(npts);
float *p = (float *) pts->GetVoidPointer(0);
//
// We will be copying the point data as we go so we need to set this up.
//
outPointData->CopyAllocate(input->GetPointData());
int pointId = 0;
// Front face.
for (i = 0 ; i < nX ; i++)
{
for (j = 0 ; j < nY ; j++)
{
p[0] = x[i];
p[1] = y[j];
p[2] = z[0];
p += 3;
outPointData->CopyData(inPointData, j*nX + i, pointId++);
}
}
// Back face
if (nZ > 1)
{
for (i = 0 ; i < nX ; i++)
{
for (j = 0 ; j < nY ; j++)
{
p[0] = x[i];
p[1] = y[j];
p[2] = z[nZ-1];
p += 3;
outPointData->CopyData(inPointData, (nZ-1)*nX*nY + j*nX + i,pointId++);
}
}
}
// Bottom face
for (i = 0 ; i < nX ; i++)
if (same && type == VTK_FLOAT)
{
for (j = 1 ; j < nZ-1 ; j++)
{
p[0] = x[i];
p[1] = y[0];
p[2] = z[j];
p += 3;
outPointData->CopyData(inPointData, j*nX*nY + i, pointId++);
}
vtkRectilinearGridFacelistFilter_ProcessFaces(nX, nY, nZ,
outPointData, inPointData,
vtkDirectAccessor<float>(xc),
vtkDirectAccessor<float>(yc),
vtkDirectAccessor<float>(zc),
vtkDirectAccessor<float>(pts->GetData()));
}
// Top face
if (nY > 1)
else if (same && type == VTK_DOUBLE)
{
for (i = 0 ; i < nX ; i++)
{
for (j = 1 ; j < nZ-1 ; j++)
{
p[0] = x[i];
p[1] = y[nY-1];
p[2] = z[j];
p += 3;
outPointData->CopyData(inPointData, j*nX*nY + (nY-1)*nX + i,pointId++);
}
}
vtkRectilinearGridFacelistFilter_ProcessFaces(nX, nY, nZ,
outPointData, inPointData,
vtkDirectAccessor<double>(xc),
vtkDirectAccessor<double>(yc),
vtkDirectAccessor<double>(zc),
vtkDirectAccessor<double>(pts->GetData()));
}
// Left face
for (i = 1 ; i < nY-1 ; i++)
else
{
for (j = 1 ; j < nZ-1 ; j++)
{
p[0] = x[0];
p[1] = y[i];
p[2] = z[j];
p += 3;
outPointData->CopyData(inPointData, j*nX*nY + i*nX, pointId++);
}
vtkRectilinearGridFacelistFilter_ProcessFaces(nX, nY, nZ,
outPointData, inPointData,
vtkGeneralAccessor(xc),
vtkGeneralAccessor(yc),
vtkGeneralAccessor(zc),
vtkGeneralAccessor(pts->GetData()));
}
// Right face
if (nX > 1)
{
for (i = 1 ; i < nY-1 ; i++)
{
for (j = 1 ; j < nZ-1 ; j++)
{
p[0] = x[nX-1];
p[1] = y[i];
p[2] = z[j];
p += 3;
outPointData->CopyData(inPointData, j*nX*nY + i*nX + nX-1, pointId++);
}
}
}
delete [] x;
delete [] y;
delete [] z;
output->SetPoints(pts);
pts->Delete();
......@@ -907,7 +945,7 @@ void vtkRectilinearGridFacelistFilter::ConsolidateFacesWithoutGhostZones(void)
outCellData->CopyData(inCellData, 0, i);
outPointData->CopyAllocate(inPointData, numOutPoints);
vtkPoints *pts = vtkPoints::New();
vtkPoints *pts = vtkVisItUtility::NewPoints(input);
pts->SetNumberOfPoints(numOutPoints);
for (i = 0 ; i < numOutPoints ; i++)
{
......
......@@ -38,6 +38,7 @@
#include <vtkRectilinearLinesNoDataFilter.h>
#include <vtkAccessors.h>
#include <vtkCellArray.h>
#include <vtkCellData.h>
#include <vtkObjectFactory.h>
......@@ -73,11 +74,15 @@ vtkRectilinearLinesNoDataFilter::vtkRectilinearLinesNoDataFilter()
#define AddLineToPolyData(ai,aj,ak, bi,bj,bk) \
{ \
p[0]=x[ai]; p[1]=y[aj]; p[2]=z[ak]; \
p+=3; \
p.SetComponent(0, x.GetComponent(ai)); \
p.SetComponent(1, y.GetComponent(aj)); \
p.SetComponent(2, z.GetComponent(ak)); \
++p; \
outPointData->CopyData(inPointData, ((ak*nY) + aj)*nX + ai, pointId++); \
p[0]=x[bi]; p[1]=y[bj]; p[2]=z[bk]; \
p+=3; \
p.SetComponent(0, x.GetComponent(bi)); \
p.SetComponent(1, y.GetComponent(bj)); \
p.SetComponent(2, z.GetComponent(bk)); \
++p; \
outPointData->CopyData(inPointData, ((bk*nY) + bj)*nX + bi, pointId++); \
*nl++ = 2; \
*nl++ = pointId-2; \
......@@ -85,7 +90,102 @@ vtkRectilinearLinesNoDataFilter::vtkRectilinearLinesNoDataFilter()
cellId++; \
}
template <class Accessor> inline void
vtkRectilinearLinesNoDataFilter_AddLines(int nX, int nY, int nZ,
vtkIdType *nl, vtkPointData *outPointData, vtkPointData *inPointData,
Accessor x, Accessor y, Accessor z, Accessor p)
{
//
// And now actually create the points/lines
//
int pointId = 0;
int cellId = 0;
p.InitTraversal();
// This case is mutually exclusive with the other ones below....
if ((nX==1 && nY==1) || (nX==1 && nZ==1) || (nY==1 && nZ==1))
{
AddLineToPolyData(0,0,0, nX-1,nY-1,nZ-1);
}
if (nX>1 && nY>1)
{
// even if nz==1
{
// Front, Top to Bottom
for (int i = 0 ; i < nX ; i++)
AddLineToPolyData(i ,0 ,0 , i ,nY-1,0 );
// Front, Left to Right
for (int j = 0 ; j < nY ; j++)
AddLineToPolyData(0 ,j ,0 , nX-1,j ,0 );
}
if (nZ>1)
{
// Back, Top to Bottom
for (int i = 0 ; i < nX ; i++)
AddLineToPolyData(i ,0 ,nZ-1, i ,nY-1,nZ-1);
// Back, Left to Right
for (int j = 0 ; j < nY ; j++)
AddLineToPolyData(0 ,j ,nZ-1, nX-1,j ,nZ-1);
}
}
if (nX>1 && nZ>1)
{
// even if ny==1
{
// Top, Front to Back
for (int i = 0 ; i < nX ; i++)
AddLineToPolyData(i ,0 ,0 , i ,0 ,nZ-1);
// Top, Left to Right
for (int k = 0 ; k < nZ ; k++)
AddLineToPolyData(0 ,0 ,k , nX-1,0 ,k );
}
if (nY>1)
{
// Bottom, Front to Back
for (int i = 0 ; i < nX ; i++)
AddLineToPolyData(i ,nY-1,0 , i ,nY-1,nZ-1);
// Bottom, Left to Right
for (int k = 0 ; k < nZ ; k++)
AddLineToPolyData(0 ,nY-1,k , nX-1,nY-1,k );
}
}
if (nY>1 && nZ>1)
{
// even if nx==1
{
// Left, Front to Back
for (int j = 0 ; j < nY ; j++)
AddLineToPolyData(0 ,j ,0 , 0 ,j ,nZ-1);
// Left, Top to Bottom
for (int k = 0 ; k < nZ ; k++)
AddLineToPolyData(0 ,0 ,k , 0 ,nY-1,k );
}
if (nX>1)
{
// Right, Front to Back
for (int j = 0 ; j < nY ; j++)
AddLineToPolyData(nX-1,j ,0 , nX-1,j ,nZ-1);
// Right, Top to Bottom
for (int k = 0 ; k < nZ ; k++)
AddLineToPolyData(nX-1,0 ,k , nX-1,nY-1,k );
}
}
}