Index: vtkBiQuadraticQuadraticHexahedron.h =================================================================== RCS file: /cvsroot/VTK/VTK/Filtering/vtkBiQuadraticQuadraticHexahedron.h,v retrieving revision 1.8 diff -u -r1.8 vtkBiQuadraticQuadraticHexahedron.h --- vtkBiQuadraticQuadraticHexahedron.h 27 Feb 2007 23:34:50 -0000 1.8 +++ vtkBiQuadraticQuadraticHexahedron.h 9 Aug 2007 14:52:52 -0000 @@ -26,7 +26,33 @@ // these midedge nodes correspond lie // on the edges defined by (0,1), (1,2), (2,3), (3,0), (4,5), (5,6), (6,7), // (7,4), (0,4), (1,5), (2,6), (3,7). The center face nodes lieing in quad -// 20-(0,1,5,4), 21-(1,2,6,5), 22-(2,3,7,6) and 23-(3,0,4,7) +// 22-(0,1,5,4), 21-(1,2,6,5), 23-(2,3,7,6) and 22-(3,0,4,7) +// +// \verbatim +// +// top +// 7--14--6 +// | | +// 15 13 +// | | +// 4--12--5 +// +// middle +// 19--23--18 +// | | +// 20 21 +// | | +// 16--22--17 +// +// bottom +// 3--10--2 +// | | +// 11 9 +// | | +// 0-- 8--1 +// +// \endverbatim +// // // .SECTION See Also // vtkQuadraticEdge vtkQuadraticTriangle vtkQuadraticTetra @@ -99,14 +125,14 @@ // Description: // @deprecated Replaced by vtkBiQuadraticQuadraticHexahedron::InterpolateFunctions as of VTK 5.2 - static void InterpolationFunctions(double pcoords[3], double weights[20]); + static void InterpolationFunctions(double pcoords[3], double weights[24]); // Description: // @deprecated Replaced by vtkBiQuadraticQuadraticHexahedron::InterpolateDerivs as of VTK 5.2 static void InterpolationDerivs(double pcoords[3], double derivs[72]); // Description: // Compute the interpolation functions/derivatives // (aka shape functions/derivatives) - virtual void InterpolateFunctions(double pcoords[3], double weights[20]) + virtual void InterpolateFunctions(double pcoords[3], double weights[24]) { vtkBiQuadraticQuadraticHexahedron::InterpolationFunctions(pcoords,weights); } Index: vtkBiQuadraticQuadraticHexahedron.cxx =================================================================== RCS file: /cvsroot/VTK/VTK/Filtering/vtkBiQuadraticQuadraticHexahedron.cxx,v retrieving revision 1.9 diff -u -r1.9 vtkBiQuadraticQuadraticHexahedron.cxx --- vtkBiQuadraticQuadraticHexahedron.cxx 12 May 2007 15:23:44 -0000 1.9 +++ vtkBiQuadraticQuadraticHexahedron.cxx 9 Aug 2007 14:52:52 -0000 @@ -29,7 +29,7 @@ #include "vtkBiQuadraticQuad.h" #include "vtkPoints.h" -vtkCxxRevisionMacro(vtkBiQuadraticQuadraticHexahedron, "$Revision: 1.9 $"); +vtkCxxRevisionMacro(vtkBiQuadraticQuadraticHexahedron, "$Revision: 1.8 $"); vtkStandardNewMacro(vtkBiQuadraticQuadraticHexahedron); //---------------------------------------------------------------------------- @@ -59,7 +59,8 @@ this->CellScalars = vtkDoubleArray::New(); this->CellScalars->SetNumberOfTuples(27); this->Scalars = vtkDoubleArray::New(); - this->Scalars->SetNumberOfTuples(8); + this->Scalars->SetNumberOfTuples(8); // vertices of a linear hexahedron + } //---------------------------------------------------------------------------- @@ -220,10 +221,23 @@ int i, j; double d, pt[3]; double derivs[72]; + double hexweights[8]; // set initial position for Newton's method + pcoords[0] = pcoords[1] = pcoords[2] = params[0] = params[1] = params[2]=0.0; subId = 0; - pcoords[0] = pcoords[1] = pcoords[2] = params[0] = params[1] = params[2]=0.5; + + // Use a tri-linear hexahederon to get good starting values + vtkHexahedron *hex = vtkHexahedron::New(); + for(i = 0; i < 8; i++) + hex->GetPoints()->SetPoint(i, this->Points->GetPoint(i)); + + hex->EvaluatePosition(x, closestPoint, subId, pcoords, dist2, hexweights); + hex->Delete(); + + params[0] = pcoords[0]; + params[1] = pcoords[1]; + params[2] = pcoords[2]; // enter iteration loop for (iteration=converged=0; @@ -259,6 +273,7 @@ d=vtkMath::Determinant3x3(rcol,scol,tcol); if ( fabs(d) < 1.e-20) { + vtkErrorMacro (<<"Determinant incorrect, iteration " << iteration); return -1; } @@ -279,6 +294,7 @@ (fabs(pcoords[1]) > VTK_DIVERGED) || (fabs(pcoords[2]) > VTK_DIVERGED)) { + vtkErrorMacro (<<"Newton did not converged, iteration " << iteration); return -1; } @@ -295,6 +311,7 @@ // outside of element if ( !converged ) { + vtkErrorMacro (<<"Newton did not converged, iteration " << iteration); return -1; } @@ -627,12 +644,13 @@ weights[19] =(-0.25*(x*(1-x))*(y*(1+y)) - 0.25*(1+x)*(1-x)*(1+y)*(1-y)) *((1+z)*(1-z)); //Face center Nodes in xz and yz direction - weights[22] = 0.5*((1+x)*(1-x))*(1-y) *((1+z)*(1-z)); + weights[20] = 0.5*((1+y)*(1-y))*(1-x) *((1+z)*(1-z)); weights[21] = 0.5*((1+y)*(1-y))*(1+x) *((1+z)*(1-z)); + weights[22] = 0.5*((1+x)*(1-x))*(1-y) *((1+z)*(1-z)); weights[23] = 0.5*((1+x)*(1-x))*(1+y) *((1+z)*(1-z)); - weights[20] = 0.5*((1+y)*(1-y))*(1-x) *((1+z)*(1-z)); } + //---------------------------------------------------------------------------- // Derivatives in parametric space. void vtkBiQuadraticQuadraticHexahedron::InterpolationDerivs(double pcoords[3], @@ -645,92 +663,86 @@ double y = 2.0*(pcoords[1]-0.5); double z = 2.0*(pcoords[2]-0.5); - //x-direction - //The eight corner points - derivs[0] =( 0.25*(1-2*x)*(y*(1-y)) - 0.25*(-2*x)*(1+y)*(1-y))*(-0.5*z*(1-z)); - derivs[1] =(-0.25*(1+2*x)*(y*(1-y)) - 0.25*(-2*x)*(1+y)*(1-y))*(-0.5*z*(1-z)); - derivs[2] =( 0.25*(1+2*x)*(y*(1+y)) - 0.25*(-2*x)*(1+y)*(1-y))*(-0.5*z*(1-z)); - derivs[3] =(-0.25*(1-2*x)*(y*(1+y)) - 0.25*(-2*x)*(1+y)*(1-y))*(-0.5*z*(1-z)); - derivs[4] =( 0.25*(1-2*x)*(y*(1-y)) - 0.25*(-2*x)*(1+y)*(1-y))*( 0.5*z*(1+z)); - derivs[5] =(-0.25*(1+2*x)*(y*(1-y)) - 0.25*(-2*x)*(1+y)*(1-y))*( 0.5*z*(1+z)); - derivs[6] =( 0.25*(1+2*x)*(y*(1+y)) - 0.25*(-2*x)*(1+y)*(1-y))*( 0.5*z*(1+z)); - derivs[7] =(-0.25*(1-2*x)*(y*(1+y)) - 0.25*(-2*x)*(1+y)*(1-y))*( 0.5*z*(1+z)); - //The mid-edge nodes - derivs[8] = 0.5*((-2*x))*(1-y) * (-0.5*z*(1-z)); - derivs[9] = 0.5*((1+y)*(1-y)) * (-0.5*z*(1-z)); - derivs[10] = 0.5*((-2*x))*(1+y) * (-0.5*z*(1-z)); - derivs[11] =-0.5*((1+y)*(1-y)) * (-0.5*z*(1-z)); - derivs[12] = 0.5*((-2*x))*(1-y) * ( 0.5*z*(1+z)); - derivs[13] = 0.5*((1+y)*(1-y)) * ( 0.5*z*(1+z)); - derivs[14] = 0.5*((-2*x))*(1+y) * ( 0.5*z*(1+z)); - derivs[15] =-0.5*((1+y)*(1-y)) * ( 0.5*z*(1+z)); - derivs[16] =( 0.25*(1-2*x)*(y*(1-y)) - 0.25*(-2*x)*(1+y)*(1-y)) *((1+z)*(1-z)); - derivs[17] =(-0.25*(1+2*x)*(y*(1-y)) - 0.25*(-2*x)*(1+y)*(1-y)) *((1+z)*(1-z)); - derivs[18] =( 0.25*(1+2*x)*(y*(1+y)) - 0.25*(-2*x)*(1+y)*(1-y)) *((1+z)*(1-z)); - derivs[19] =(-0.25*(1-2*x)*(y*(1+y)) - 0.25*(-2*x)*(1+y)*(1-y)) *((1+z)*(1-z)); - //Face center Nodes in xz and yz direction - derivs[20] = 0.5*((-2*x))*(1-y) *((1+z)*(1-z)); - derivs[21] = 0.5*((1+y)*(1-y)) *((1+z)*(1-z)); - derivs[22] = 0.5*((-2*x))*(1+y) *((1+z)*(1-z)); - derivs[23] =-0.5*((1+y)*(1-y)) *((1+z)*(1-z)); - - //y-direction - //The eight Corner points - derivs[24] =( 0.25*(x*(1-x))*(1-2*y) - 0.25*(1+x)*(1-x)*(-2*y))*(-0.5*z*(1-z)); - derivs[25] =(-0.25*(x*(1+x))*(1-2*y) - 0.25*(1+x)*(1-x)*(-2*y))*(-0.5*z*(1-z)); - derivs[26] =( 0.25*(x*(1+x))*(1+2*y) - 0.25*(1+x)*(1-x)*(-2*y))*(-0.5*z*(1-z)); - derivs[27] =(-0.25*(x*(1-x))*(1+2*y) - 0.25*(1+x)*(1-x)*(-2*y))*(-0.5*z*(1-z)); - derivs[28] =( 0.25*(x*(1-x))*(1-2*y) - 0.25*(1+x)*(1-x)*(-2*y))*( 0.5*z*(1+z)); - derivs[29] =(-0.25*(x*(1+x))*(1-2*y) - 0.25*(1+x)*(1-x)*(-2*y))*( 0.5*z*(1+z)); - derivs[30] =( 0.25*(x*(1+x))*(1+2*y) - 0.25*(1+x)*(1-x)*(-2*y))*( 0.5*z*(1+z)); - derivs[31] =(-0.25*(x*(1-x))*(1+2*y) - 0.25*(1+x)*(1-x)*(-2*y))*( 0.5*z*(1+z)); - //The mid-edge nodes - derivs[32] =-0.5*((1+x)*(1-x)) *(-0.5*z*(1-z)); - derivs[33] = 0.5*((-2*y))*(1+x) *(-0.5*z*(1-z)); - derivs[34] = 0.5*((1+x)*(1-x)) *(-0.5*z*(1-z)); - derivs[35] = 0.5*((-2*y))*(1-x) *(-0.5*z*(1-z)); - derivs[36] =-0.5*((1+x)*(1-x)) *( 0.5*z*(1+z)); - derivs[37] = 0.5*((-2*y))*(1+x) *( 0.5*z*(1+z)); - derivs[38] = 0.5*((1+x)*(1-x)) *( 0.5*z*(1+z)); - derivs[39] = 0.5*((-2*y))*(1-x) *( 0.5*z*(1+z)); - derivs[40] =( 0.25*(x*(1-x))*(1-2*y) - 0.25*(1+x)*(1-x)*(-2*y)) *((1+z)*(1-z)); - derivs[41] =(-0.25*(x*(1+x))*(1-2*y) - 0.25*(1+x)*(1-x)*(-2*y)) *((1+z)*(1-z)); - derivs[42] =( 0.25*(x*(1+x))*(1+2*y) - 0.25*(1+x)*(1-x)*(-2*y)) *((1+z)*(1-z)); - derivs[43] =(-0.25*(x*(1-x))*(1+2*y) - 0.25*(1+x)*(1-x)*(-2*y)) *((1+z)*(1-z)); - //Face center Nodes in xz and yz direction - derivs[44] =-0.5*((1+x)*(1-x)) * ((1+z)*(1-z)); - derivs[45] = 0.5*((-2*y))*(1+x) * ((1+z)*(1-z)); - derivs[46] = 0.5*((1+x)*(1-x)) * ((1+z)*(1-z)); - derivs[47] = 0.5*((-2*y))*(1-x) * ((1+z)*(1-z)); - - //z-direction - //The eight corner points - derivs[48] =( 0.25*(x*(1-x))*(y*(1-y)) - 0.25*(1+x)*(1-x)*(1+y)*(1-y))*(-0.5*(1-2*z)); - derivs[49] =(-0.25*(x*(1+x))*(y*(1-y)) - 0.25*(1+x)*(1-x)*(1+y)*(1-y))*(-0.5*(1-2*z)); - derivs[50] =( 0.25*(x*(1+x))*(y*(1+y)) - 0.25*(1+x)*(1-x)*(1+y)*(1-y))*(-0.5*(1-2*z)); - derivs[51] =(-0.25*(x*(1-x))*(y*(1+y)) - 0.25*(1+x)*(1-x)*(1+y)*(1-y))*(-0.5*(1-2*z)); - derivs[52] =( 0.25*(x*(1-x))*(y*(1-y)) - 0.25*(1+x)*(1-x)*(1+y)*(1-y))*( 0.5*(1+2*z)); - derivs[53] =(-0.25*(x*(1+x))*(y*(1-y)) - 0.25*(1+x)*(1-x)*(1+y)*(1-y))*( 0.5*(1+2*z)); - derivs[54] =( 0.25*(x*(1+x))*(y*(1+y)) - 0.25*(1+x)*(1-x)*(1+y)*(1-y))*( 0.5*(1+2*z)); - derivs[55] =(-0.25*(x*(1-x))*(y*(1+y)) - 0.25*(1+x)*(1-x)*(1+y)*(1-y))*( 0.5*(1+2*z)); - //The mid-edge nodes - derivs[56] = 0.5*((1+x)*(1-x))*(1-y) *(-0.5*(1-2*z)); - derivs[57] = 0.5*((1+y)*(1-y))*(1+x) *(-0.5*(1-2*z)); - derivs[58] = 0.5*((1+x)*(1-x))*(1+y) *(-0.5*(1-2*z)); - derivs[59] = 0.5*((1+y)*(1-y))*(1-x) *(-0.5*(1-2*z)); - derivs[60] = 0.5*((1+x)*(1-x))*(1-y) *( 0.5*(1+2*z)); - derivs[61] = 0.5*((1+y)*(1-y))*(1+x) *( 0.5*(1+2*z)); - derivs[62] = 0.5*((1+x)*(1-x))*(1+y) *( 0.5*(1+2*z)); - derivs[63] = 0.5*((1+y)*(1-y))*(1-x) *( 0.5*(1+2*z)); - derivs[64] =( 0.25*(x*(1-x))*(y*(1-y)) - 0.25*(1+x)*(1-x)*(1+y)*(1-y)) *((-2*z)); - derivs[65] =(-0.25*(x*(1+x))*(y*(1-y)) - 0.25*(1+x)*(1-x)*(1+y)*(1-y)) *((-2*z)); - derivs[66] =( 0.25*(x*(1+x))*(y*(1+y)) - 0.25*(1+x)*(1-x)*(1+y)*(1-y)) *((-2*z)); - derivs[67] =(-0.25*(x*(1-x))*(y*(1+y)) - 0.25*(1+x)*(1-x)*(1+y)*(1-y)) *((-2*z)); - //Face center Nodes in xz and yz direction - derivs[68] =0.5*((1+x)*(1-x))*(1-y) *((-2*z)); - derivs[69] =0.5*((1+y)*(1-y))*(1+x) *((-2*z)); - derivs[70] =0.5*((1+x)*(1-x))*(1+y) *((-2*z)); - derivs[71] =0.5*((1+y)*(1-y))*(1-x) *((-2*z)); + // x direction + derivs[0] =-((y*y+(2*x-1)*y-2*x)*z*z+(-y*y+(1-2*x)*y+2*x)*z)/8; + derivs[1] =((y*y+(-2*x-1)*y+2*x)*z*z+(-y*y+(2*x+1)*y-2*x)*z)/8; + derivs[2] =((y*y+(2*x+1)*y+2*x)*z*z+(-y*y+(-2*x-1)*y-2*x)*z)/8; + derivs[3] =-((y*y+(1-2*x)*y-2*x)*z*z+(-y*y+(2*x-1)*y+2*x)*z)/8; + derivs[4] =-((y*y+(2*x-1)*y-2*x)*z*z+(y*y+(2*x-1)*y-2*x)*z)/8; + derivs[5] =((y*y+(-2*x-1)*y+2*x)*z*z+(y*y+(-2*x-1)*y+2*x)*z)/8; + derivs[6] =((y*y+(2*x+1)*y+2*x)*z*z+(y*y+(2*x+1)*y+2*x)*z)/8; + derivs[7] =-((y*y+(1-2*x)*y-2*x)*z*z+(y*y+(1-2*x)*y-2*x)*z)/8; + derivs[8] =((x*y-x)*z*z+(x-x*y)*z)/2; + derivs[9] =-((y*y-1)*z*z+(1-y*y)*z)/4; + derivs[10] =-((x*y+x)*z*z+(-x*y-x)*z)/2; + derivs[11] =((y*y-1)*z*z+(1-y*y)*z)/4; + derivs[12] =((x*y-x)*z*z+(x*y-x)*z)/2; + derivs[13] =-((y*y-1)*z*z+(y*y-1)*z)/4; + derivs[14] =-((x*y+x)*z*z+(x*y+x)*z)/2; + derivs[15] =((y*y-1)*z*z+(y*y-1)*z)/4; + derivs[16] =((y*y+(2*x-1)*y-2*x)*z*z-y*y+(1-2*x)*y+2*x)/4; + derivs[17] =-((y*y+(-2*x-1)*y+2*x)*z*z-y*y+(2*x+1)*y-2*x)/4; + derivs[18] =-((y*y+(2*x+1)*y+2*x)*z*z-y*y+(-2*x-1)*y-2*x)/4; + derivs[19] =((y*y+(1-2*x)*y-2*x)*z*z-y*y+(2*x-1)*y+2*x)/4; + derivs[20] =-((y*y-1)*z*z-y*y+1)/2; + derivs[21] =((y*y-1)*z*z-y*y+1)/2; + derivs[22] =(x-x*y)*z*z+x*y-x; + derivs[23] =(x*y+x)*z*z-x*y-x; + // y direction + derivs[24] =-(((2*x-2)*y+x*x-x)*z*z+((2-2*x)*y-x*x+x)*z)/8; + derivs[25] =(((2*x+2)*y-x*x-x)*z*z+((-2*x-2)*y+x*x+x)*z)/8; + derivs[26] =(((2*x+2)*y+x*x+x)*z*z+((-2*x-2)*y-x*x-x)*z)/8; + derivs[27] =-(((2*x-2)*y-x*x+x)*z*z+((2-2*x)*y+x*x-x)*z)/8; + derivs[28] =-(((2*x-2)*y+x*x-x)*z*z+((2*x-2)*y+x*x-x)*z)/8; + derivs[29] =(((2*x+2)*y-x*x-x)*z*z+((2*x+2)*y-x*x-x)*z)/8; + derivs[30] =(((2*x+2)*y+x*x+x)*z*z+((2*x+2)*y+x*x+x)*z)/8; + derivs[31] =-(((2*x-2)*y-x*x+x)*z*z+((2*x-2)*y-x*x+x)*z)/8; + derivs[32] =((x*x-1)*z*z+(1-x*x)*z)/4; + derivs[33] =-((x+1)*y*z*z+(-x-1)*y*z)/2; + derivs[34] =-((x*x-1)*z*z+(1-x*x)*z)/4; + derivs[35] =((x-1)*y*z*z+(1-x)*y*z)/2; + derivs[36] =((x*x-1)*z*z+(x*x-1)*z)/4; + derivs[37] =-((x+1)*y*z*z+(x+1)*y*z)/2; + derivs[38] =-((x*x-1)*z*z+(x*x-1)*z)/4; + derivs[39] =((x-1)*y*z*z+(x-1)*y*z)/2; + derivs[40] =(((2*x-2)*y+x*x-x)*z*z+(2-2*x)*y-x*x+x)/4; + derivs[41] =-(((2*x+2)*y-x*x-x)*z*z+(-2*x-2)*y+x*x+x)/4; + derivs[42] =-(((2*x+2)*y+x*x+x)*z*z+(-2*x-2)*y-x*x-x)/4; + derivs[43] =(((2*x-2)*y-x*x+x)*z*z+(2-2*x)*y+x*x-x)/4; + derivs[44] =(1-x)*y*z*z+(x-1)*y; + derivs[45] =(x+1)*y*z*z+(-x-1)*y; + derivs[46] =-((x*x-1)*z*z-x*x+1)/2; + derivs[47] =((x*x-1)*z*z-x*x+1)/2; + // z direction + derivs[48] =-(((2*x-2)*y*y+(2*x*x-2*x)*y-2*x*x+2)*z+(1-x)*y*y+(x-x*x)*y+x*x-1)/8; + derivs[49] =(((2*x+2)*y*y+(-2*x*x-2*x)*y+2*x*x-2)*z+(-x-1)*y*y+(x*x+x)*y-x*x+1)/8; + derivs[50] =(((2*x+2)*y*y+(2*x*x+2*x)*y+2*x*x-2)*z+(-x-1)*y*y+(-x*x-x)*y-x*x+1)/8; + derivs[51] =-(((2*x-2)*y*y+(2*x-2*x*x)*y-2*x*x+2)*z+(1-x)*y*y+(x*x-x)*y+x*x-1)/8; + derivs[52] =-(((2*x-2)*y*y+(2*x*x-2*x)*y-2*x*x+2)*z+(x-1)*y*y+(x*x-x)*y-x*x+1)/8; + derivs[53] =(((2*x+2)*y*y+(-2*x*x-2*x)*y+2*x*x-2)*z+(x+1)*y*y+(-x*x-x)*y+x*x-1)/8; + derivs[54] =(((2*x+2)*y*y+(2*x*x+2*x)*y+2*x*x-2)*z+(x+1)*y*y+(x*x+x)*y+x*x-1)/8; + derivs[55] =-(((2*x-2)*y*y+(2*x-2*x*x)*y-2*x*x+2)*z+(x-1)*y*y+(x-x*x)*y-x*x+1)/8; + derivs[56] =(((2*x*x-2)*y-2*x*x+2)*z+(1-x*x)*y+x*x-1)/4; + derivs[57] =-(((2*x+2)*y*y-2*x-2)*z+(-x-1)*y*y+x+1)/4; + derivs[58] =-(((2*x*x-2)*y+2*x*x-2)*z+(1-x*x)*y-x*x+1)/4; + derivs[59] =(((2*x-2)*y*y-2*x+2)*z+(1-x)*y*y+x-1)/4; + derivs[60] =(((2*x*x-2)*y-2*x*x+2)*z+(x*x-1)*y-x*x+1)/4; + derivs[61] =-(((2*x+2)*y*y-2*x-2)*z+(x+1)*y*y-x-1)/4; + derivs[62] =-(((2*x*x-2)*y+2*x*x-2)*z+(x*x-1)*y+x*x-1)/4; + derivs[63] =(((2*x-2)*y*y-2*x+2)*z+(x-1)*y*y-x+1)/4; + derivs[64] =((x-1)*y*y+(x*x-x)*y-x*x+1)*z/2; + derivs[65] =-((x+1)*y*y+(-x*x-x)*y+x*x-1)*z/2; + derivs[66] =-((x+1)*y*y+(x*x+x)*y+x*x-1)*z/2; + derivs[67] =((x-1)*y*y+(x-x*x)*y-x*x+1)*z/2; + derivs[68] =((1-x)*y*y+x-1)*z; + derivs[69] =((x+1)*y*y-x-1)*z; + derivs[70] =((1-x*x)*y+x*x-1)*z; + derivs[71] =((x*x-1)*y+x*x-1)*z; + + // we compute derivatives in in [-1; 1] but we need them in [ 0; 1] + for(int i = 0; i < 72; i++) + derivs[i] *= 2; + } //----------------------------------------------------------------------------