From ebbfb1406423887dbaed49719d5327851a0d5841 Mon Sep 17 00:00:00 2001 From: eap Date: Wed, 20 Jun 2007 14:52:43 +0000 Subject: [PATCH] PAL16231 (3D Extrusion produced crossing hexahedrons) improve copmuting point parameters performance by implementing it without math_FunctionSetWithDerivatives usage --- src/SMESH/SMESH_Block.cxx | 305 +++++++++++--------------------------- src/SMESH/SMESH_Block.hxx | 35 +---- 2 files changed, 87 insertions(+), 253 deletions(-) diff --git a/src/SMESH/SMESH_Block.cxx b/src/SMESH/SMESH_Block.cxx index 6f45ca370..dba5e0d32 100644 --- a/src/SMESH/SMESH_Block.cxx +++ b/src/SMESH/SMESH_Block.cxx @@ -382,75 +382,24 @@ bool SMESH_Block::ShellPoint(const gp_XYZ& theParams, const vector& p = thePointOnShape; thePoint = - x1 * p[ID_F0yz] + x * p[ID_F1yz] - + y1 * p[ID_Fx0z] + y * p[ID_Fx1z] - + z1 * p[ID_Fxy0] + z * p[ID_Fxy1] - + x1 * (y1 * (z1 * p[ID_V000] + z * p[ID_V001]) - + y * (z1 * p[ID_V010] + z * p[ID_V011])) - + x * (y1 * (z1 * p[ID_V100] + z * p[ID_V101]) - + y * (z1 * p[ID_V110] + z * p[ID_V111])); + x1 * p[ID_F0yz] + x * p[ID_F1yz] + + y1 * p[ID_Fx0z] + y * p[ID_Fx1z] + + z1 * p[ID_Fxy0] + z * p[ID_Fxy1] + + x1 * (y1 * (z1 * p[ID_V000] + z * p[ID_V001]) + + y * (z1 * p[ID_V010] + z * p[ID_V011])) + + x * (y1 * (z1 * p[ID_V100] + z * p[ID_V101]) + + y * (z1 * p[ID_V110] + z * p[ID_V111])); thePoint -= - x1 * (y1 * p[ID_E00z] + y * p[ID_E01z]) - + x * (y1 * p[ID_E10z] + y * p[ID_E11z]) - + y1 * (z1 * p[ID_Ex00] + z * p[ID_Ex01]) - + y * (z1 * p[ID_Ex10] + z * p[ID_Ex11]) - + z1 * (x1 * p[ID_E0y0] + x * p[ID_E1y0]) - + z * (x1 * p[ID_E0y1] + x * p[ID_E1y1]); + x1 * (y1 * p[ID_E00z] + y * p[ID_E01z]) + + x * (y1 * p[ID_E10z] + y * p[ID_E11z]) + + y1 * (z1 * p[ID_Ex00] + z * p[ID_Ex01]) + + y * (z1 * p[ID_Ex10] + z * p[ID_Ex11]) + + z1 * (x1 * p[ID_E0y0] + x * p[ID_E1y0]) + + z * (x1 * p[ID_E0y1] + x * p[ID_E1y1]); return true; } -//======================================================================= -//function : NbVariables -//purpose : -//======================================================================= - -Standard_Integer SMESH_Block::NbVariables() const -{ - return 3; -} - -//======================================================================= -//function : NbEquations -//purpose : -//======================================================================= - -Standard_Integer SMESH_Block::NbEquations() const -{ - return 1; -} - -//======================================================================= -//function : Value -//purpose : -//======================================================================= - -Standard_Boolean SMESH_Block::Value(const math_Vector& theXYZ, math_Vector& theFxyz) -{ - gp_XYZ P, params( theXYZ(1), theXYZ(2), theXYZ(3) ); - if ( params.IsEqual( myParam, DBL_MIN )) { // same param - theFxyz( 1 ) = funcValue( myValues[ SQUARE_DIST ]); - } - else { - ShellPoint( params, P ); - gp_Vec dP( P - myPoint ); - theFxyz(1) = funcValue( dP.SquareMagnitude() ); - } - return true; -} - -//======================================================================= -//function : Derivatives -//purpose : -//======================================================================= - -Standard_Boolean SMESH_Block::Derivatives(const math_Vector& XYZ,math_Matrix& Df) -{ - MESSAGE( "SMESH_Block::Derivatives()"); - math_Vector F(1,3); - return Values(XYZ,F,Df); -} - //======================================================================= //function : Constructor //purpose : @@ -463,101 +412,6 @@ SMESH_Block::SMESH_Block(): { } -//======================================================================= -//function : Values -//purpose : -//======================================================================= - -Standard_Boolean SMESH_Block::Values(const math_Vector& theXYZ, - math_Vector& theFxyz, - math_Matrix& theDf) -{ - gp_XYZ P, params( theXYZ(1), theXYZ(2), theXYZ(3) ); - if ( params.IsEqual( myParam, DBL_MIN )) { // same param - theFxyz( 1 ) = funcValue( myValues[ SQUARE_DIST ] ); - theDf( 1, DRV_1 ) = myValues[ DRV_1 ]; - theDf( 1, DRV_2 ) = myValues[ DRV_2 ]; - theDf( 1, DRV_3 ) = myValues[ DRV_3 ]; - return true; - } -#ifdef DEBUG_PARAM_COMPUTE - cout << "PARAM GUESS: " << params.X() << " "<< params.Y() << " "<< params.X() << endl; - myNbIterations++; // how many times call ShellPoint() -#endif - ShellPoint( params, P ); - - gp_Vec dP( myPoint, P ); - double sqDist = dP.SquareMagnitude(); - theFxyz(1) = funcValue( sqDist ); - - if ( sqDist < myTolerance * myTolerance ) { // a solution found - myParam = params; - myValues[ SQUARE_DIST ] = sqDist; - theFxyz(1) = theDf( 1,1 ) = theDf( 1,2 ) = theDf( 1,3 ) = 0; - return true; - } - - if ( sqDist < myValues[ SQUARE_DIST ] ) // a better guess - { - // 3 partial derivatives - gp_Vec drv[ 3 ]; // where we move with a small step in each direction - for ( int iP = 1; iP <= 3; iP++ ) { - if ( iP - 1 == myFaceIndex ) { - drv[ iP - 1 ] = gp_Vec(0,0,0); - continue; - } - gp_XYZ Pi; - bool onEdge = ( theXYZ( iP ) + 0.001 > 1. ); - if ( onEdge ) - params.SetCoord( iP, theXYZ( iP ) - 0.001 ); - else - params.SetCoord( iP, theXYZ( iP ) + 0.001 ); - ShellPoint( params, Pi ); - params.SetCoord( iP, theXYZ( iP ) ); // restore params - gp_Vec dPi ( P, Pi ); - if ( onEdge ) dPi *= -1.; - double mag = dPi.Magnitude(); - if ( mag > DBL_MIN ) - dPi /= mag; - drv[ iP - 1 ] = dPi; - } - for ( int iP = 0; iP < 3; iP++ ) { -#if 1 - theDf( 1, iP + 1 ) = dP * drv[iP]; -#else - // Distance from P to plane passing through myPoint and defined - // by the 2 other derivative directions: - // like IntAna_IntConicQuad::Perform (const gp_Lin& L, const gp_Pln& P) - // where L is (P -> myPoint), P is defined by the 2 other derivative direction - int iPrev = ( iP ? iP - 1 : 2 ); - int iNext = ( iP == 2 ? 0 : iP + 1 ); - gp_Vec plnNorm = drv[ iPrev ].Crossed( drv [ iNext ] ); - double Direc = plnNorm * drv[ iP ]; - if ( Abs(Direc) <= DBL_MIN ) - theDf( 1, iP + 1 ) = dP * drv[ iP ]; - else { - double Dis = plnNorm * P - plnNorm * myPoint; - theDf( 1, iP + 1 ) = Dis/Direc; - } -#endif - } -#ifdef DEBUG_PARAM_COMPUTE - cout << "F = " << theFxyz(1) << - " DRV: " << theDf(1,1) << " " << theDf(1,2) << " " << theDf(1,3) << endl; - myNbIterations +=3; // how many times call ShellPoint() -#endif - - // store better values - myParam = params; - myValues[SQUARE_DIST]= sqDist; - myValues[DRV_1] = theDf(1,DRV_1); - myValues[DRV_2] = theDf(1,DRV_2); - myValues[DRV_3] = theDf(1,DRV_3); - } - - return true; -} - //======================================================================= //function : ComputeParameters //purpose : compute point parameters in the block @@ -583,17 +437,12 @@ bool SMESH_Block::ComputeParameters(const gp_Pnt& thePoint, return false; } - myPoint = thePoint.XYZ(); - - myParam.SetCoord( -1,-1,-1 ); - myValues[ SQUARE_DIST ] = 1e100; - const bool isOnFace = IsFaceID( theShapeID ); double * coef = GetShapeCoef( theShapeID ); // Find the first guess paremeters - math_Vector start( 1, 3, 0.0 ); + gp_XYZ start(0, 0, 0); bool hasHint = ( 0 <= theParamsHint.X() && theParamsHint.X() <= 1 && 0 <= theParamsHint.Y() && theParamsHint.Y() <= 1 && @@ -611,6 +460,7 @@ bool SMESH_Block::ComputeParameters(const gp_Pnt& thePoint, iEdge += 4; continue; } + double sumParam = 0; for ( int iE = 0; iE < 4; iE++, iEdge++ ) { // loop on 4 parallel edges gp_Pnt p0 = myEdge[ iEdge ].Point( par000 ); gp_Pnt p1 = myEdge[ iEdge ].Point( par111 ); @@ -624,9 +474,9 @@ bool SMESH_Block::ComputeParameters(const gp_Pnt& thePoint, break; } } - start( iParam ) += par; + sumParam += par; } - start( iParam ) /= 4.; + start.SetCoord( iParam, sumParam / 4.); } if ( needGrid ) { // compute nodes of 3 x 3 x 3 grid @@ -644,22 +494,10 @@ bool SMESH_Block::ComputeParameters(const gp_Pnt& thePoint, myTolerance = sqrt( box.SquareExtent() ) * 1e-5; } } - if ( myTolerance < 0 ) myTolerance = 1e-6; - double loopTol = 10 * myTolerance; // to stop restarting solution search if ( hasHint ) { - gp_Pnt hintPoint = myPoint; - ShellPoint( theParamsHint, hintPoint.ChangeCoord() ); - double sqDist = hintPoint.SquareDistance( myPoint ); - if ( sqrt( sqDist ) < loopTol ) { - // theParamsHint is good, no need to solve - myValues[ SQUARE_DIST ] = sqDist; - myParam = theParamsHint; - } - start( 1 ) = theParamsHint.X(); - start( 2 ) = theParamsHint.Y(); - start( 3 ) = theParamsHint.Z(); + start = theParamsHint; } else if ( myGridComputed ) { @@ -673,60 +511,97 @@ bool SMESH_Block::ComputeParameters(const gp_Pnt& thePoint, bestParam = & prmPtn.first; } } - start( 1 ) = bestParam->X(); - start( 2 ) = bestParam->Y(); - start( 3 ) = bestParam->Z(); + start = *bestParam; } - myFaceIndex = -1; + int faceIndex = -1; + double faceParam = 0.; if ( isOnFace ) { // put a point on the face for ( int iCoord = 0; iCoord < 3; iCoord++ ) if ( coef[ iCoord ] ) { - myFaceIndex = iCoord; - myFaceParam = ( coef[ myFaceIndex ] < 0.5 ) ? 0.0 : 1.0; - start( iCoord + 1 ) = myFaceParam; + faceIndex = iCoord; + faceParam = ( coef[ faceIndex ] < 0.5 ) ? 0.0 : 1.0; + start.SetCoord( iCoord + 1, faceParam ); } } #ifdef DEBUG_PARAM_COMPUTE cout << " #### POINT " < loopTol && nbLoops <= 3 ) + while ( nbLoops <= 100 ) { - paramSearch.Perform ( *static_cast(this), - start, low, up ); - start( 1 ) = myParam.X(); - start( 2 ) = myParam.Y(); - start( 3 ) = myParam.Z(); - mySquareFunc = !mySquareFunc; + gp_XYZ P; + ShellPoint( params, P ); + + gp_Vec dP( thePoint, P ); + double sqDist = dP.SquareMagnitude(); +#ifdef DEBUG_PARAM_COMPUTE + cout << "PARAMS: ( " << params.X() <<" "<< params.Y() <<" "<< params.Z() <<" )"<< endl; + cout << "DIST: " << sqrt( sqDist ) << endl; +#endif + if ( sqDist > sqDistance ) { // solution get worse + break; + } + sqDistance = sqDist; + solution = params; + if ( sqDist < sqTolerance ) { // a solution found + break; + } + // look for a better solution + for ( int iP = 1; iP <= 3; iP++ ) { + if ( iP - 1 == faceIndex ) + continue; + // see where we move with a small (=parDelta) step in this direction + gp_XYZ nearParams = params; + bool onEdge = ( params.Coord( iP ) + parDelta > 1. ); + if ( onEdge ) + nearParams.SetCoord( iP, params.Coord( iP ) - parDelta ); + else + nearParams.SetCoord( iP, params.Coord( iP ) + parDelta ); + gp_XYZ Pi; + ShellPoint( nearParams, Pi ); + gp_Vec dPi ( P, Pi ); + if ( onEdge ) dPi *= -1.; + // modify a parameter + double mag = dPi.Magnitude(); + if ( mag < DBL_MIN ) + continue; + gp_Vec dir = dPi / mag; // dir we move modifying the parameter + double dist = dir * dP; // where we should get to + double dPar = dist / mag * parDelta; // predict parameter change + double curPar = params.Coord( iP ); + double par = curPar - dPar; // new parameter value + while ( par > 1 || par < 0 ) { + dPar /= 2.; + par = curPar - dPar; + } + params.SetCoord( iP, par ); + } + nbLoops++; } #ifdef DEBUG_PARAM_COMPUTE - mySumDist += distance(); - cout << " ------ SOLUTION: ( "<< myParam.X() <<" "<< myParam.Y() <<" "<< myParam.Z() <<" )"<= 0 ) - theParams.SetCoord( myFaceIndex + 1, myFaceParam ); + if ( faceIndex >= 0 ) + theParams.SetCoord( faceIndex + 1, faceParam ); return true; } @@ -767,16 +642,6 @@ bool SMESH_Block::EdgeParameters(const int theEdgeID, const double theU, gp_XYZ& return false; } -//======================================================================= -//function : GetStateNumber -//purpose : -//======================================================================= - -Standard_Integer SMESH_Block::GetStateNumber () -{ - return 0; //myValues[0] < 1e-1; -} - //======================================================================= //function : DumpShapeID //purpose : debug an id of a block sub-shape diff --git a/src/SMESH/SMESH_Block.hxx b/src/SMESH/SMESH_Block.hxx index 7ebe54292..c4e8b2c28 100644 --- a/src/SMESH/SMESH_Block.hxx +++ b/src/SMESH/SMESH_Block.hxx @@ -25,10 +25,6 @@ #ifndef SMESH_Block_HeaderFile #define SMESH_Block_HeaderFile -// #include -// #include -// #include - #include #include #include @@ -36,12 +32,8 @@ #include #include #include -#include #include #include -#include -#include -#include #include #include @@ -58,7 +50,7 @@ class Adaptor3d_Curve; // parameters inside the block and vice versa // ========================================================= -class SMESH_Block: public math_FunctionSetWithDerivatives +class SMESH_Block { public: enum TShapeID { @@ -282,17 +274,6 @@ public: // theFirstVertex may be NULL. // Always try to set a seam edge first - public: - // ----------------------------------------------------------- - // Methods of math_FunctionSetWithDerivatives used internally - // to define parameters by coordinates - // ----------------------------------------------------------- - Standard_Integer NbVariables() const; - Standard_Integer NbEquations() const; - Standard_Boolean Value(const math_Vector& X,math_Vector& F) ; - Standard_Boolean Derivatives(const math_Vector& X,math_Matrix& D) ; - Standard_Boolean Values(const math_Vector& X,math_Vector& F,math_Matrix& D) ; - Standard_Integer GetStateNumber (); protected: @@ -362,20 +343,8 @@ public: // for param computation - enum { SQUARE_DIST = 0, DRV_1, DRV_2, DRV_3 }; - double distance () const { return sqrt( myValues[ SQUARE_DIST ]); } - double funcValue(double sqDist) const { return mySquareFunc ? sqDist : sqrt(sqDist); } - - int myFaceIndex; - double myFaceParam; int myNbIterations; - double mySumDist; - double myTolerance; - bool mySquareFunc; - - gp_XYZ myPoint; // the given point - gp_XYZ myParam; // the best parameters guess - double myValues[ 4 ]; // values computed at myParam: square distance and 3 derivatives + double mySumDist, myTolerance; typedef pair TxyzPair; TxyzPair my3x3x3GridNodes[ 27 ]; // to compute the first param guess -- 2.30.2