Commit a8a43d4a authored by Bill Lorensen's avatar Bill Lorensen Committed by Kitware Robot

Merge topic 'UnitTestCells'

8044bf4e ENH: Unit Test for Linear, Quadratic and BiQuadratic Cells
Acked-by: Kitware Robot's avatarKitware Robot <kwrobot@kitware.com>
Reviewed-by: Berk Geveci's avatarBerk Geveci <berk.geveci@kitware.com>
Merge-request: !338
parents b0691ec1 8044bf4e
......@@ -46,6 +46,7 @@ vtk_add_test_cxx(${vtk-module}CxxTests tests
TestStructuredData.cxx
TestDataObjectTypes.cxx
TestPolyDataRemoveDeletedCells.cxx
UnitTestCells.cxx
)
vtk_add_test_cxx(${vtk-module}CxxTests data_tests
TestCellIterators.cxx,NO_VALID,NO_OUTPUT
......
......@@ -77,6 +77,7 @@ int TestOneInterpolationFunction(double eps = VTK_EPSILON)
{
if( fabs(sf[j] - 1) > eps)
{
std::cout << "fabs(sf[" << j << "] - 1): " << fabs(sf[j] - 1) << std::endl;
++r;
}
}
......@@ -84,6 +85,7 @@ int TestOneInterpolationFunction(double eps = VTK_EPSILON)
{
if( fabs(sf[j] - 0) > eps )
{
std::cout << "fabs(sf[" << j << "] - 0): " << fabs(sf[j] - 0) << std::endl;
++r;
}
}
......
This diff is collapsed.
......@@ -414,12 +414,9 @@ vtkBiQuadraticQuad::Derivatives (int vtkNotUsed (subId),
void vtkBiQuadraticQuad::InterpolationFunctions (double pcoords[3],
double weights[9])
{
//VTK needs parametric coordinates to be between (0,1). Isoparametric
//shape functions are normaly formulated between (-1,1). But because we need
//to derivate these functions in x and y direction, we are formulating the
//shape functions in parametric coordinates. Normaly these coordinates are
//named r and s, but I chose x and y, because you can easily mark and
//paste these functions to the gnuplot splot function. :D
//Normally these coordinates are named r and s, but I chose x and y,
//because you can easily mark and paste these functions to the
//gnuplot splot function. :D
double x = pcoords[0];
double y = pcoords[1];
......
......@@ -275,7 +275,7 @@ int vtkBiQuadraticQuadraticHexahedron::EvaluatePosition(double* x,
d=vtkMath::Determinant3x3(rcol,scol,tcol);
if ( fabs(d) < 1.e-20)
{
vtkErrorMacro (<<"Determinant incorrect, iteration " << iteration);
vtkDebugMacro (<<"Determinant incorrect, iteration " << iteration);
return -1;
}
......@@ -296,7 +296,6 @@ int vtkBiQuadraticQuadraticHexahedron::EvaluatePosition(double* x,
(fabs(pcoords[1]) > VTK_DIVERGED) ||
(fabs(pcoords[2]) > VTK_DIVERGED))
{
vtkErrorMacro (<<"Newton did not converged, iteration " << iteration);
return -1;
}
......@@ -313,7 +312,6 @@ int vtkBiQuadraticQuadraticHexahedron::EvaluatePosition(double* x,
// outside of element
if ( !converged )
{
vtkErrorMacro (<<"Newton did not converged, iteration " << iteration);
return -1;
}
......@@ -549,7 +547,6 @@ void vtkBiQuadraticQuadraticHexahedron::JacobianInverse(double pcoords[3],
// now find the inverse
if ( vtkMath::InvertMatrix(m,inverse,3) == 0 )
{
vtkErrorMacro(<<"Jacobian inverse not found");
return;
}
}
......
......@@ -196,6 +196,7 @@ int vtkBiQuadraticQuadraticWedge::EvaluatePosition (double *x,
d = vtkMath::Determinant3x3 (rcol, scol, tcol);
if (fabs (d) < 1.e-20)
{
vtkDebugMacro (<<"Determinant incorrect, iteration " << iteration);
return -1;
}
......
......@@ -151,7 +151,7 @@ int vtkBiQuadraticTriangle::EvaluatePosition(double* x, double* closestPoint,
pcoords[0] = 0.5*pc0;
pcoords[1] = 0.5 + 0.5*pc1;
}
pcoords[2] = 1.0 - pcoords[0] - pcoords[1];
pcoords[2] = 0.0;
this->InterpolationFunctions(pcoords,weights);
}
......
......@@ -133,6 +133,7 @@ int vtkHexagonalPrism::EvaluatePosition(double x[3], double* closestPoint,
d=vtkMath::Determinant3x3(rcol,scol,tcol);
if ( fabs(d) < 1.e-20)
{
vtkDebugMacro (<<"Determinant incorrect, iteration " << iteration);
return -1;
}
......
......@@ -134,6 +134,7 @@ int vtkPentagonalPrism::EvaluatePosition(double x[3], double closestPoint[3],
d=vtkMath::Determinant3x3(rcol,scol,tcol);
if ( fabs(d) < 1.e-20)
{
vtkDebugMacro (<<"Determinant incorrect, iteration " << iteration);
return -1;
}
......
......@@ -618,7 +618,6 @@ int vtkPolygon::PointInPolygon (double x[3], int numPts, double *pts,
int iterNumber;
int maxComp, comps[2];
int deltaVotes;
// do a quick bounds check
if ( x[0] < bounds[0] || x[0] > bounds[1] ||
x[1] < bounds[2] || x[1] > bounds[3] ||
......@@ -735,9 +734,18 @@ int vtkPolygon::PointInPolygon (double x[3], int numPts, double *pts,
// Fire the ray and compute the number of intersections. Be careful
// of degenerate cases (e.g., ray intersects at vertex).
//
if ((status=vtkLine::Intersection(x,xray,x1,x2,u,v)) == VTK_POLYGON_INTERSECTION)
{
if ( (VTK_POLYGON_RAY_TOL < v) && (v < 1.0-VTK_POLYGON_RAY_TOL) )
// This test checks for vertex and edge intersections
// For example
// Vertex intersection
// (u=0 v=0), (u=0 v=1), (u=1 v=0), (u=1 v=0)
// Edge intersection
// (u=0 v!=0 v!=1), (u=1 v!=0 v!=1)
// (u!=0 u!=1 v=0), (u!=0 u!=1 v=1)
if ( (VTK_POLYGON_RAY_TOL < u) && (u < 1.0-VTK_POLYGON_RAY_TOL) &&
(VTK_POLYGON_RAY_TOL < v) && (v < 1.0-VTK_POLYGON_RAY_TOL) )
{
numInts++;
}
......@@ -753,7 +761,7 @@ int vtkPolygon::PointInPolygon (double x[3], int numPts, double *pts,
}
if ( testResult == VTK_POLYGON_CERTAIN )
{
if ( (numInts % 2) == 0)
if ( numInts % 2 == 0)
{
--deltaVotes;
}
......
......@@ -218,4 +218,3 @@ private:
};
#endif
......@@ -1909,7 +1909,23 @@ int vtkPolyhedron::IsInside(double x[3], double tolerance)
this->PolyData->GetCell(this->CellIds->GetId(idx), this->Cell);
if ( this->Cell->IntersectWithLine(x, xray, tol, t, xint, pcoords, subId) )
{
numInts++;
// Check for vertex, edge or face intersections
// count the number of 0 or 1 pcoords
int pcount = 0;
for (int p = 0; p < 3; ++p)
{
if (pcoords[p] == 0.0 || pcoords[p] == 1.0)
{
pcount++;
}
}
// pcount = 1, exact face intersection
// pcount = 2, exact edge intersection
// pcount = 3, exact vertex intersection
if (pcount == 0)
{
numInts++;
}
}
} //for all candidate cells
}
......@@ -1923,13 +1939,29 @@ int vtkPolyhedron::IsInside(double x[3], double tolerance)
this->PolyData->GetCell(idx, this->Cell);
if ( this->Cell->IntersectWithLine(x, xray, tol, t, xint, pcoords, subId) )
{
numInts++;
// Check for vertex, edge or face intersections
// count the number of 0 or 1 pcoords
int pcount = 0;
for (int p = 0; p < 3; ++p)
{
if (pcoords[p] == 0.0 || pcoords[p] == 1.0)
{
pcount++;
}
}
// pcount = 1, exact face intersection
// pcount = 2, exact edge intersection
// pcount = 3, exact vertex intersection
if (pcount == 0)
{
numInts++;
}
}
} //for all candidate cells
}
// Count the result
if ( (numInts % 2) == 0)
if ( numInts != 0 && (numInts % 2) == 0)
{
--deltaVotes;
}
......
......@@ -95,7 +95,7 @@ int vtkPyramid::EvaluatePosition(double x[3], double closestPoint[3],
// square it here because we're looking at dist2^2.
if(dist2 == 0. || ( length2 != 0. && dist2/length2 < 1.e-6) )
{
pcoords[0] = pcoords[1] = .5;
pcoords[0] = pcoords[1] = 0;
pcoords[2] = 1;
this->InterpolationFunctions(pcoords, weights);
if(closestPoint)
......@@ -150,6 +150,7 @@ int vtkPyramid::EvaluatePosition(double x[3], double closestPoint[3],
d=vtkMath::Determinant3x3(rcol,scol,tcol);
if ( fabs(d) < 1.e-20)
{
vtkDebugMacro (<<"Determinant incorrect, iteration " << iteration);
return -1;
}
......
......@@ -199,8 +199,8 @@ void vtkQuadraticHexahedron::Subdivide(vtkPointData *inPd, vtkCellData *inCd,
//----------------------------------------------------------------------------
static const double VTK_DIVERGED = 1.e6;
static const int VTK_HEX_MAX_ITERATION=10;
static const double VTK_HEX_CONVERGED=1.e-03;
static const int VTK_HEX_MAX_ITERATION=20;
static const double VTK_HEX_CONVERGED=1.e-04;
int vtkQuadraticHexahedron::EvaluatePosition(double* x,
double* closestPoint,
......@@ -252,6 +252,7 @@ int vtkQuadraticHexahedron::EvaluatePosition(double* x,
d=vtkMath::Determinant3x3(rcol,scol,tcol);
if ( fabs(d) < 1.e-20)
{
vtkDebugMacro (<<"Determinant incorrect, iteration " << iteration);
return -1;
}
......
......@@ -484,10 +484,9 @@ void vtkQuadraticPolygon::ConvertFromPolygon(vtkIdList *ids)
//----------------------------------------------------------------------------
void vtkQuadraticPolygon::Derivatives(int vtkNotUsed(subId),
double pcoords[3],
double vtkNotUsed(pcoords)[3],
double *vtkNotUsed(values),
int vtkNotUsed(dim),
double *vtkNotUsed(derivs))
{
pcoords[0] = pcoords[1] = pcoords[2] = 0.0;
}
......@@ -160,10 +160,50 @@ int vtkQuadraticPyramid::EvaluatePosition(double* x,
int& subId, double pcoords[3],
double& dist2, double *weights)
{
int i, j;
subId = 0;
// There are problems searching for the apex point so we check if
// we are there first before doing the full parametric inversion.
vtkPoints* points = this->GetPoints();
double apexPoint[3];
points->GetPoint(4, apexPoint);
dist2 = vtkMath::Distance2BetweenPoints(apexPoint, x);
double baseMidpoint[3];
points->GetPoint(0, baseMidpoint);
for(i=1;i<4;i++)
{
double tmp[3];
points->GetPoint(i, tmp);
for(j=0;j<3;j++)
{
baseMidpoint[j] += tmp[j];
}
}
for(i=0;i<3;i++)
{
baseMidpoint[i] /= 4.;
}
double length2 = vtkMath::Distance2BetweenPoints(apexPoint, baseMidpoint);
// we use .001 as the relative tolerance here since that is the same
// that is used for the interior cell check below but we need to
// square it here because we're looking at dist2^2.
if(dist2 == 0. || ( length2 != 0. && dist2/length2 < 1.e-6) )
{
pcoords[0] = pcoords[1] = 0;
pcoords[2] = 1;
this->InterpolationFunctions(pcoords, weights);
if(closestPoint)
{
memcpy(closestPoint, x, 3*sizeof(double));
dist2 = 0.;
}
return 1;
}
int iteration, converged;
double params[3];
double fcol[3], rcol[3], scol[3], tcol[3];
int i, j;
double d, pt[3];
double derivs[3*13];
......@@ -205,6 +245,7 @@ int vtkQuadraticPyramid::EvaluatePosition(double* x,
d=vtkMath::Determinant3x3(rcol,scol,tcol);
if ( fabs(d) < 1.e-20)
{
vtkDebugMacro (<<"Determinant incorrect, iteration " << iteration);
return -1;
}
......@@ -646,39 +687,30 @@ void vtkQuadraticPyramid::InterpolationFunctions(double pcoords[3],
// VTK needs parametric coordinates to be between (0,1). Isoparametric
// shape functions are formulated between (-1,1). Here we do a
// coordinate system conversion from (0,1) to (-1,1).
const double r = 2*pcoords[0] - 1;
const double s = 2*pcoords[1] - 1;
const double t = 2*pcoords[2] - 1;
const double rm = 1.0 - r;
const double rp = 1.0 + r;
const double sm = 1.0 - s;
const double sp = 1.0 + s;
const double tm = 1.0 - t;
const double tp = 1.0 + t;
const double r2 = 1.0 - r*r;
const double s2 = 1.0 - s*s;
const double t2 = 1.0 - t*t;
const double r = 2.0*(pcoords[0] - 0.5);
const double s = 2.0*(pcoords[1] - 0.5);
const double t = 2.0*(pcoords[2] - 0.5);
const double r2 = r*r;
const double s2 = s*s;
const double t2 = t*t;
weights[0] = -(1 - r) * (1 - s) * (1 - t) * (4 + 3*r + 3*s + 2*r*s + 2*t + r*t + s*t + 2*r*s*t) / 16.0;
weights[1] = -(1 + r) * (1 - s) * (1 - t) * (4 - 3*r + 3*s - 2*r*s + 2*t - r*t + s*t - 2*r*s*t) / 16.0;
weights[2] = -(1 + r) * (1 + s) * (1 - t) * (4 - 3*r - 3*s + 2*r*s + 2*t - r*t - s*t + 2*r*s*t) / 16.0;
weights[3] = -(1 - r) * (1 + s) * (1 - t) * (4 + 3*r - 3*s - 2*r*s + 2*t + r*t - s*t - 2*r*s*t) / 16.0;
weights[4] = t * (1 + t) / 2.0;
weights[5] = (1 - r2) * (1 - s) * (1 - t) * (2 + s + s*t) / 8.0;
weights[6] = (1 + r) * (1 - s2) * (1 - t) * (2 - r - r*t) / 8.0;
weights[7] = (1 - r2) * (1 + s) * (1 - t) * (2 - s - s*t) / 8.0;
weights[8] = (1 - r) * (1 - s2) * (1 - t) * (2 + r + r*t) / 8.0;
weights[9] = (1 - r) * (1 - s) * (1 - t2) / 4.0;
weights[10] = (1 + r) * (1 - s) * (1 - t2) / 4.0;
weights[11] = (1 + r) * (1 + s) * (1 - t2) / 4.0;
weights[12] = (1 - r) * (1 + s) * (1 - t2) / 4.0;
// corners
weights[0] = 0.125 * rm * sm * tm * (-r - s - t - 2.0);
weights[1] = 0.125 * rp * sm * tm * ( r - s - t - 2.0);
weights[2] = 0.125 * rp * sp * tm * ( r + s - t - 2.0);
weights[3] = 0.125 * rm * sp * tm * (-r + s - t - 2.0);
weights[4] = 0.5 * t * tp;
// midsides of rectangles
weights[5] = 0.25 * r2 * sm * tm;
weights[6] = 0.25 * s2 * rp * tm;
weights[7] = 0.25 * r2 * sp * tm;
weights[8] = 0.25 * s2 * rm * tm;
// midsides of triangles
weights[9] = 0.25 * ( 1 - r ) * (1 - s ) * t2;
weights[10] = 0.25 * ( 1 + r ) * (1 - s ) * t2;
weights[11] = 0.25 * ( 1 + r ) * (1 + s ) * t2;
weights[12] = 0.25 * ( 1 - r ) * (1 + s ) * t2;
}
//----------------------------------------------------------------------------
......
......@@ -155,8 +155,8 @@ private:
//
inline int vtkQuadraticPyramid::GetParametricCenter(double pcoords[3])
{
pcoords[0] = pcoords[1] = 6./13;
pcoords[2] = 3./13;
pcoords[0] = pcoords[1] = 6.0/13.0;
pcoords[2] = 3.0/13.0;
return 0;
}
......
......@@ -528,9 +528,6 @@ void vtkQuadraticQuad::Derivatives(int vtkNotUsed(subId),
void vtkQuadraticQuad::InterpolationFunctions(double pcoords[3],
double weights[8])
{
//VTK needs parametric coordinates to be between (0,1). Isoparametric
//shape functions are formulated between (-1,1). Here we do a
//coordinate system conversion from (0,1) to (-1,1).
double r = pcoords[0];
double s = pcoords[1];
......
......@@ -116,8 +116,8 @@ vtkCell *vtkQuadraticTetra::GetFace(int faceId)
//----------------------------------------------------------------------------
static const double VTK_DIVERGED = 1.e6;
static const int VTK_TETRA_MAX_ITERATION=10;
static const double VTK_TETRA_CONVERGED=1.e-03;
static const int VTK_TETRA_MAX_ITERATION=20;
static const double VTK_TETRA_CONVERGED=1.e-04;
int vtkQuadraticTetra::EvaluatePosition(double* x,
double* closestPoint,
......
......@@ -127,7 +127,7 @@ int vtkQuadraticTriangle::EvaluatePosition(double* x, double* closestPoint,
pcoords[0] = 0.5 - pcoords[0]/2.0;
pcoords[1] = 0.5 - pcoords[1]/2.0;
}
pcoords[2] = 1.0 - pcoords[0] - pcoords[1];
pcoords[2] = 0.0;
if(closestPoint!=0)
{
// Compute both closestPoint and weights
......@@ -260,12 +260,11 @@ int vtkQuadraticTriangle::Triangulate(int vtkNotUsed(index), vtkIdList *ptIds,
//----------------------------------------------------------------------------
void vtkQuadraticTriangle::Derivatives(int vtkNotUsed(subId),
double pcoords[3],
double vtkNotUsed(pcoords)[3],
double *vtkNotUsed(values),
int vtkNotUsed(dim),
double *vtkNotUsed(derivs))
{
pcoords[0] = pcoords[1] = pcoords[2] = 0.0;
}
......
......@@ -202,6 +202,7 @@ int vtkQuadraticWedge::EvaluatePosition(double* x,
d=vtkMath::Determinant3x3(rcol,scol,tcol);
if ( fabs(d) < 1.e-20)
{
vtkDebugMacro (<<"Determinant incorrect, iteration " << iteration);
return -1;
}
......
......@@ -190,7 +190,7 @@ int vtkTriQuadraticHexahedron::EvaluatePosition (double *x,
d = vtkMath::Determinant3x3 (rcol, scol, tcol);
if (fabs (d) < 1.e-20)
{
vtkErrorMacro (<<"Determinant incorrect, iteration " << iteration);
vtkDebugMacro (<<"Determinant incorrect, iteration " << iteration);
return -1;
}
......@@ -210,7 +210,6 @@ int vtkTriQuadraticHexahedron::EvaluatePosition (double *x,
else if ((fabs (pcoords[0]) > VTK_DIVERGED) ||
(fabs (pcoords[1]) > VTK_DIVERGED) || (fabs (pcoords[2]) > VTK_DIVERGED))
{
vtkErrorMacro (<<"Newton did not converged, iteration " << iteration << " det " << d);
return -1;
}
......@@ -228,7 +227,6 @@ int vtkTriQuadraticHexahedron::EvaluatePosition (double *x,
// outside of element
if (!converged)
{
vtkErrorMacro (<<"Newton did not converged, iteration " << iteration << " det " << d);
return -1;
}
......
......@@ -110,6 +110,7 @@ int vtkWedge::EvaluatePosition(double x[3], double* closestPoint,
d=vtkMath::Determinant3x3(rcol,scol,tcol);
if ( fabs(d) < 1.e-20)
{
vtkDebugMacro (<<"Determinant incorrect, iteration " << iteration);
return -1;
}
......
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