Commit 971b2ba6 authored by Will Schroeder's avatar Will Schroeder
Browse files

ENH: Got tensors working

parent 1a79eef5
......@@ -15,13 +15,13 @@ Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen 1993, 1994
=========================================================================*/
#include "FTensors.hh"
vtkTensor &vtkFloatTensors::GetTensor(int i)
vtkTensor *vtkFloatTensors::GetTensor(int i)
{
static vtkTensor t;
t.SetDimension(this->Dimension);
t.T = this->T.GetPtr(this->Dimension*this->Dimension*i);
return t;
return &t;
}
vtkTensors *vtkFloatTensors::MakeObject(int sze, int d, int ext)
......
......@@ -313,7 +313,7 @@ void vtkPointData::InterpolatePoint(vtkPointData *fromPd, int toId, vtkIdList *p
{
int i, j;
float s, *pv, v[3], *pn, n[3], *ptc, tc[3];
static vtkTensor tensor(3), &pt=tensor;
static vtkTensor tensor(3), *pt;
void *ud;
if ( fromPd->Scalars && this->Scalars && this->CopyScalars )
......@@ -391,9 +391,9 @@ void vtkPointData::InterpolatePoint(vtkPointData *fromPd, int toId, vtkIdList *p
pt = cellTensors->GetTensor(i);
for (j=0; j<cellTensors->GetDimension(); j++)
for (int k=0; k<cellTensors->GetDimension(); k++)
tensor.AddComponent(j,k,pt.GetComponent(j,k)*weights[i]);
tensor.AddComponent(j,k,pt->GetComponent(j,k)*weights[i]);
}
this->Tensors->InsertTensor(toId,tensor);
this->Tensors->InsertTensor(toId,&tensor);
}
if ( fromPd->UserDefined && this->UserDefined && this->CopyUserDefined )
......
......@@ -137,20 +137,21 @@ void vtkPointLoad::Execute()
xP[1] = (this->ModelBounds[2] + this->ModelBounds[3]) / 2.0;
xP[2] = this->ModelBounds[5]; // at top of box
//
// Traverse all points evaluating implicit function at each point
// Traverse all points evaluating implicit function at each point. Note that
// points are evaluated in local coordinate system of applied force.
//
twoPi = 2.0*math.Pi();
P = -this->LoadValue;
for (k=0; k<this->Dimensions[2]; k++)
{
z = this->Origin[2] + k*this->AspectRatio[2];
z = xP[2] - (this->Origin[2] + k*this->AspectRatio[2]);
for (j=0; j<this->Dimensions[1]; j++)
{
y = this->Origin[1] + j*this->AspectRatio[1];
y = xP[1] - (this->Origin[1] + j*this->AspectRatio[1]);
for (i=0; i<this->Dimensions[0]; i++)
{
x = this->Origin[0] + i*this->AspectRatio[0];
x = (this->Origin[0] + i*this->AspectRatio[0]) - xP[0];
rho = sqrt((x-xP[0])*(x-xP[0]) + (y-xP[1])*(y-xP[1]) +
(z-xP[2])*(z-xP[2]));
if ( rho < 1.0e-10 )
......@@ -165,7 +166,7 @@ void vtkPointLoad::Execute()
tensor.SetComponent(1,2,0.0);
tensor.SetComponent(2,0,0.0);
tensor.SetComponent(2,1,0.0);
newTensors->InsertNextTensor(tensor);
newTensors->InsertNextTensor(&tensor);
if ( this->ComputeEffectiveStress )
newScalars->InsertNextScalar(LARGE_FLOAT);
continue;
......@@ -197,13 +198,13 @@ void vtkPointLoad::Execute()
tensor.SetComponent(0,0, sx);
tensor.SetComponent(1,1, sy);
tensor.SetComponent(2,2, sz);
tensor.SetComponent(0,1, txy); // symmetric matrix
tensor.SetComponent(0,1, txy); // real symmetric matrix
tensor.SetComponent(1,0, txy);
tensor.SetComponent(0,2, txz);
tensor.SetComponent(2,0, txz);
tensor.SetComponent(1,2, tyz);
tensor.SetComponent(2,1, tyz);
newTensors->InsertNextTensor(tensor);
newTensors->InsertNextTensor(&tensor);
if ( this->ComputeEffectiveStress )
{
......
......@@ -20,13 +20,17 @@ Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen 1993, 1994
#include "vtkMath.hh"
// Description
// Construct object with scaling on and scale factor 1.0.
// Construct object with scaling on and scale factor 1.0. Eigenvalues are
// extracted, glyphs are colored with input scalar data, and logarithmic
// scaling is turned off.
vtkTensorGlyph::vtkTensorGlyph()
{
this->Source = NULL;
this->Scaling = 1;
this->ScaleFactor = 1.0;
this->ExtractEigenvalues = 1;
this->ColorGlyphs = 1;
this->LogScaling = 1;
}
vtkTensorGlyph::~vtkTensorGlyph()
......@@ -37,24 +41,39 @@ void vtkTensorGlyph::Execute()
{
vtkPointData *pd;
vtkTensors *inTensors;
vtkTensor *tensor;
vtkScalars *inScalars;
int numPts, numSourcePts, numSourceCells;
int inPtId, i;
int inPtId, i, j;
vtkPoints *sourcePts;
vtkNormals *sourceNormals;
vtkCellArray *sourceCells;
vtkFloatPoints *newPts;
float *x;
vtkFloatScalars *newScalars=NULL;
vtkFloatNormals *newNormals=NULL;
float *x, s;
vtkTransform trans;
vtkCell *cell;
vtkIdList *cellPts;
int npts, pts[MAX_CELL_SIZE];
int ptIncr, cellId;
vtkMath math;
vtkMatrix4x4 matrix;
double *m[3], w[3], *e[3];
double m0[3], m1[3], m2[3];
double e0[3], e1[3], e2[3];
float xv[3], yv[3], zv[3];
// set up working matrices
m[0] = m0; m[1] = m1; m[2] = m2;
e[0] = e0; e[1] = e1; e[2] = e2;
vtkDebugMacro(<<"Generating tensor glyphs");
this->Initialize();
pd = this->Input->GetPointData();
inTensors = pd->GetTensors();
inScalars = pd->GetScalars();
numPts = this->Input->GetNumberOfPoints();
if ( !inTensors || numPts < 1 )
......@@ -91,9 +110,16 @@ void vtkTensorGlyph::Execute()
// only copy scalar data through
pd = this->Source->GetPointData();
this->PointData.CopyAllOff();
this->PointData.CopyScalarsOn();
this->PointData.CopyAllocate(pd,numPts*numSourcePts);
if ( inScalars && this->ColorGlyphs )
newScalars = new vtkFloatScalars(numPts*numSourcePts);
else
{
this->PointData.CopyAllOff();
this->PointData.CopyScalarsOn();
this->PointData.CopyAllocate(pd,numPts*numSourcePts);
}
if ( sourceNormals = pd->GetNormals() )
newNormals = new vtkFloatNormals(numPts*numSourcePts);
//
// First copy all topology (transformation independent)
//
......@@ -112,6 +138,8 @@ void vtkTensorGlyph::Execute()
//
// Traverse all Input points, transforming glyph at Source points
//
trans.PreMultiply();
for (inPtId=0; inPtId < numPts; inPtId++)
{
ptIncr = inPtId * numSourcePts;
......@@ -122,23 +150,82 @@ void vtkTensorGlyph::Execute()
x = this->Input->GetPoint(inPtId);
trans.Translate(x[0], x[1], x[2]);
// extract appropriate eigenfunctions
tensor = inTensors->GetTensor(inPtId);
// compute orientation vectors and scale factors from tensor
if ( this->ExtractEigenvalues ) // extract appropriate eigenfunctions
{
for (j=0; j<3; j++)
for (i=0; i<3; i++)
m[i][j] = tensor->GetComponent(i,j);
math.SingularValueDecomposition(m,3,3,w,e);
//copy eigenvectors
xv[0] = e[0][0]; xv[1] = e[1][0]; xv[2] = e[2][0];
yv[0] = e[0][1]; yv[1] = e[1][1]; yv[2] = e[2][1];
zv[0] = e[0][2]; zv[1] = e[1][2]; zv[2] = e[2][2];
}
else //use tensor columns as eigenvectors
{
for (i=0; i<3; i++)
{
w[i] = 1.0;
xv[i] = tensor->GetComponent(i,0);
yv[i] = tensor->GetComponent(i,1);
zv[i] = tensor->GetComponent(i,2);
}
}
// eigenvectors (assumed normalized) rotate object
// normalize eigenvectors / compute scale factors
w[0] = fabs(w[0] * math.Normalize(xv) * this->ScaleFactor);
w[1] = fabs(w[1] * math.Normalize(yv) * this->ScaleFactor);
w[2] = fabs(w[2] * math.Normalize(zv) * this->ScaleFactor);
if ( this->LogScaling )
{
for (i=0; i<3; i++)
if ( w[i] != 0.0 )
if ( (w[i]=log10(w[i])) < 0.0 ) w[i] = -1.0/w[i];
}
// if scaling modify matrix to scale according to eigenvalues
// normalized eigenvectors rotate object
matrix.Element[0][0] = xv[0];
matrix.Element[0][1] = yv[0];
matrix.Element[0][2] = zv[0];
matrix.Element[1][0] = xv[1];
matrix.Element[1][1] = yv[1];
matrix.Element[1][2] = zv[1];
matrix.Element[2][0] = xv[2];
matrix.Element[2][1] = yv[2];
matrix.Element[2][2] = zv[2];
trans.Concatenate(matrix);
trans.Scale(w[0], w[1], w[2]);
// multiply points by resulting matrix
// multiply points (and normals if available) by resulting matrix
trans.MultiplyPoints(sourcePts,newPts);
if ( newNormals )
trans.MultiplyNormals(sourceNormals,newNormals);
// Copy point data from source
for (i=0; i < numSourcePts; i++)
this->PointData.CopyData(pd,i,ptIncr+i);
if ( inScalars && this->ColorGlyphs )
{
s = inScalars->GetScalar(inPtId);
for (i=0; i < numSourcePts; i++)
newScalars->InsertScalar(ptIncr+i, s);
}
else
{
for (i=0; i < numSourcePts; i++)
this->PointData.CopyData(pd,i,ptIncr+i);
}
}
//
// Update ourselves
//
this->SetPoints(newPts);
if ( newScalars ) this->PointData.SetScalars(newScalars);
if ( newNormals ) this->PointData.SetNormals(newNormals);
this->Squeeze();
}
......@@ -184,5 +271,8 @@ void vtkTensorGlyph::PrintSelf(ostream& os, vtkIndent indent)
os << indent << "Source: " << this->Source << "\n";
os << indent << "Scaling: " << (this->Scaling ? "On\n" : "Off\n");
os << indent << "Scale Factor: " << this->ScaleFactor << "\n";
os << indent << "Extract Eigenvalues: " << (this->ExtractEigenvalues ? "On\n" : "Off\n");
os << indent << "Color Glyphs: " << (this->ColorGlyphs ? "On\n" : "Off\n");
os << indent << "Log Scaling: " << (this->LogScaling ? "On\n" : "Off\n");
}
......@@ -22,8 +22,8 @@ vtkTensors::vtkTensors(int dim)
void vtkTensors::GetTensor(int id, vtkTensor &ft)
{
vtkTensor& t = this->GetTensor(id);
ft = t;
vtkTensor *t = this->GetTensor(id);
ft = *t;
}
void vtkTensors::InsertTensor(int id, float t11, float t12, float t13,
......
......@@ -22,6 +22,7 @@ vtkWarpScalar::vtkWarpScalar()
this->Normal[0] = 0.0;
this->Normal[1] = 0.0;
this->Normal[2] = 1.0;
this->XYPlane = 0;
}
float *vtkWarpScalar::DataNormal(int id, vtkNormals *normals)
......@@ -34,6 +35,12 @@ float *vtkWarpScalar::InstanceNormal(int id, vtkNormals *normals)
return this->Normal;
}
float *vtkWarpScalar::ZNormal(int id, vtkNormals *normals)
{
static float zNormal[3]={0.0,0.0,1.0};
return zNormal;
}
void vtkWarpScalar::Execute()
{
vtkPoints *inPts;
......@@ -64,6 +71,11 @@ void vtkWarpScalar::Execute()
this->PointNormal = &vtkWarpScalar::DataNormal;
vtkDebugMacro(<<"Using data normals");
}
else if ( this->XYPlane )
{
this->PointNormal = &vtkWarpScalar::ZNormal;
vtkDebugMacro(<<"Using x-y plane normal");
}
else
{
this->PointNormal = &vtkWarpScalar::InstanceNormal;
......@@ -78,7 +90,8 @@ void vtkWarpScalar::Execute()
{
x = inPts->GetPoint(ptId);
n = (this->*(PointNormal))(ptId,inNormals);
s = inScalars->GetScalar(ptId);
if ( this->XYPlane ) s = x[2];
else s = inScalars->GetScalar(ptId);
for (i=0; i<3; i++)
{
newX[i] = x[i] + this->ScaleFactor * s * n[i];
......@@ -102,4 +115,5 @@ void vtkWarpScalar::PrintSelf(ostream& os, vtkIndent indent)
os << indent << "Use Normal: " << (this->UseNormal ? "On\n" : "Off\n");
os << indent << "Normal: (" << this->Normal[0] << ", "
<< this->Normal[1] << ", " << this->Normal[2] << ")\n";
os << indent << "XY Plane: " << (this->XYPlane ? "On\n" : "Off\n");
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment