X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FSMESH%2FSMESH_Block.cxx;h=dba5e0d32d833575799db38c44281c6a8410568c;hp=6550a81ea3391ed5fc3ef15bb6cb28fbdef2679f;hb=ebbfb1406423887dbaed49719d5327851a0d5841;hpb=ed456586bfb1411c5bff73b221658766689a6253 diff --git a/src/SMESH/SMESH_Block.cxx b/src/SMESH/SMESH_Block.cxx index 6550a81ea..dba5e0d32 100644 --- a/src/SMESH/SMESH_Block.cxx +++ b/src/SMESH/SMESH_Block.cxx @@ -15,7 +15,7 @@ // License along with this library; if not, write to the Free Software // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA // -// See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // File : SMESH_Pattern.hxx // Created : Mon Aug 2 10:30:00 2004 @@ -23,16 +23,23 @@ #include "SMESH_Block.hxx" +#include +#include +#include #include #include +#include #include -#include +#include +#include +#include +#include #include -#include #include #include #include #include +#include #include #include #include @@ -40,13 +47,55 @@ #include #include +#include "SMDS_MeshNode.hxx" +#include "SMDS_MeshVolume.hxx" +#include "SMDS_VolumeTool.hxx" #include "utilities.h" #include using namespace std; -#define SQRT_FUNC 1 +//#define DEBUG_PARAM_COMPUTE + +//================================================================================ +/*! + * \brief Set edge data + * \param edgeID - block subshape ID + * \param curve - edge geometry + * \param isForward - is curve orientation coincides with edge orientation in the block + */ +//================================================================================ + +void SMESH_Block::TEdge::Set( const int edgeID, Adaptor3d_Curve* curve, const bool isForward ) +{ + myCoordInd = SMESH_Block::GetCoordIndOnEdge( edgeID ); + if ( myC3d ) delete myC3d; + myC3d = curve; + myFirst = curve->FirstParameter(); + myLast = curve->LastParameter(); + if ( !isForward ) + std::swap( myFirst, myLast ); +} + +//================================================================================ +/*! + * \brief Set coordinates of nodes at edge ends to work with mesh block + * \param edgeID - block subshape ID + * \param node1 - coordinates of node with lower ID + * \param node2 - coordinates of node with upper ID + */ +//================================================================================ + +void SMESH_Block::TEdge::Set( const int edgeID, const gp_XYZ& node1, const gp_XYZ& node2 ) +{ + myCoordInd = SMESH_Block::GetCoordIndOnEdge( edgeID ); + myNodes[ 0 ] = node1; + myNodes[ 1 ] = node2; + + if ( myC3d ) delete myC3d; + myC3d = 0; +} //======================================================================= //function : SMESH_Block::TEdge::GetU @@ -56,6 +105,8 @@ using namespace std; double SMESH_Block::TEdge::GetU( const gp_XYZ& theParams ) const { double u = theParams.Coord( myCoordInd ); + if ( !myC3d ) // if mesh block + return u; return ( 1 - u ) * myFirst + u * myLast; } @@ -66,10 +117,125 @@ double SMESH_Block::TEdge::GetU( const gp_XYZ& theParams ) const gp_XYZ SMESH_Block::TEdge::Point( const gp_XYZ& theParams ) const { - gp_XYZ p = myC3d->Value( GetU( theParams )).XYZ(); - if ( myTrsf.Form() != gp_Identity ) - myTrsf.Transforms( p ); - return p; + double u = GetU( theParams ); + if ( myC3d ) return myC3d->Value( u ).XYZ(); + // mesh block + return myNodes[0] * ( 1 - u ) + myNodes[1] * u; +} + +//================================================================================ +/*! + * \brief Destructor + */ +//================================================================================ + +SMESH_Block::TEdge::~TEdge() +{ + if ( myC3d ) delete myC3d; +} + +//================================================================================ +/*! + * \brief Set face data + * \param faceID - block subshape ID + * \param S - face surface geometry + * \param c2d - 4 pcurves in the order as returned by GetFaceEdgesIDs(faceID) + * \param isForward - orientation of pcurves comparing with block edge direction + */ +//================================================================================ + +void SMESH_Block::TFace::Set( const int faceID, + Adaptor3d_Surface* S, + Adaptor2d_Curve2d* c2D[4], + const bool isForward[4] ) +{ + if ( myS ) delete myS; + myS = S; + // pcurves + vector< int > edgeIdVec; + GetFaceEdgesIDs( faceID, edgeIdVec ); + for ( int iE = 0; iE < edgeIdVec.size(); iE++ ) // loop on 4 edges + { + myCoordInd[ iE ] = GetCoordIndOnEdge( edgeIdVec[ iE ] ); + if ( myC2d[ iE ]) delete myC2d[ iE ]; + myC2d[ iE ] = c2D[ iE ]; + myFirst[ iE ] = myC2d[ iE ]->FirstParameter(); + myLast [ iE ] = myC2d[ iE ]->LastParameter(); + if ( !isForward[ iE ]) + std::swap( myFirst[ iE ], myLast[ iE ] ); + } + // 2d corners + myCorner[ 0 ] = myC2d[ 0 ]->Value( myFirst[0] ).XY(); + myCorner[ 1 ] = myC2d[ 0 ]->Value( myLast[0] ).XY(); + myCorner[ 2 ] = myC2d[ 1 ]->Value( myLast[1] ).XY(); + myCorner[ 3 ] = myC2d[ 1 ]->Value( myFirst[1] ).XY(); +} + +//================================================================================ +/*! + * \brief Set face data to work with mesh block + * \param faceID - block subshape ID + * \param edgeU0 - filled data of edge u0 = GetFaceEdgesIDs(faceID)[ 0 ] + * \param edgeU1 - filled data of edge u1 = GetFaceEdgesIDs(faceID)[ 1 ] + */ +//================================================================================ + +void SMESH_Block::TFace::Set( const int faceID, const TEdge& edgeU0, const TEdge& edgeU1 ) +{ + vector< int > edgeIdVec; + GetFaceEdgesIDs( faceID, edgeIdVec ); + myNodes[ 0 ] = edgeU0.NodeXYZ( 1 ); + myNodes[ 1 ] = edgeU0.NodeXYZ( 0 ); + myNodes[ 2 ] = edgeU1.NodeXYZ( 0 ); + myNodes[ 3 ] = edgeU1.NodeXYZ( 1 ); + myCoordInd[ 0 ] = GetCoordIndOnEdge( edgeIdVec[ 0 ] ); + myCoordInd[ 1 ] = GetCoordIndOnEdge( edgeIdVec[ 1 ] ); + myCoordInd[ 2 ] = GetCoordIndOnEdge( edgeIdVec[ 2 ] ); + myCoordInd[ 3 ] = GetCoordIndOnEdge( edgeIdVec[ 3 ] ); + if ( myS ) delete myS; + myS = 0; +} + +//================================================================================ +/*! + * \brief Destructor + */ +//================================================================================ + +SMESH_Block::TFace::~TFace() +{ + if ( myS ) delete myS; + for ( int i = 0 ; i < 4; ++i ) + if ( myC2d[ i ]) delete myC2d[ i ]; +} + +//======================================================================= +//function : SMESH_Block::TFace::GetCoefs +//purpose : return coefficients for addition of [0-3]-th edge and vertex +//======================================================================= + +void SMESH_Block::TFace::GetCoefs(int iE, + const gp_XYZ& theParams, + double& Ecoef, + double& Vcoef ) const +{ + double dU = theParams.Coord( GetUInd() ); + double dV = theParams.Coord( GetVInd() ); + switch ( iE ) { + case 0: + Ecoef = ( 1 - dV ); // u0 + Vcoef = ( 1 - dU ) * ( 1 - dV ); break; // 00 + case 1: + Ecoef = dV; // u1 + Vcoef = dU * ( 1 - dV ); break; // 10 + case 2: + Ecoef = ( 1 - dU ); // 0v + Vcoef = dU * dV ; break; // 11 + case 3: + Ecoef = dU ; // 1v + Vcoef = ( 1 - dU ) * dV ; break; // 01 + default: ASSERT(0); + } } //======================================================================= @@ -80,26 +246,10 @@ gp_XYZ SMESH_Block::TEdge::Point( const gp_XYZ& theParams ) const gp_XY SMESH_Block::TFace::GetUV( const gp_XYZ& theParams ) const { gp_XY uv(0.,0.); - double dU = theParams.Coord( GetUInd() ); - double dV = theParams.Coord( GetVInd() ); for ( int iE = 0; iE < 4; iE++ ) // loop on 4 edges { double Ecoef = 0, Vcoef = 0; - switch ( iE ) { - case 0: - Ecoef = ( 1 - dV ); // u0 - Vcoef = ( 1 - dU ) * ( 1 - dV ); break; // 00 - case 1: - Ecoef = dV; // u1 - Vcoef = dU * ( 1 - dV ); break; // 10 - case 2: - Ecoef = ( 1 - dU ); // 0v - Vcoef = dU * dV ; break; // 11 - case 3: - Ecoef = dU ; // 1v - Vcoef = ( 1 - dU ) * dV ; break; // 01 - default:; - } + GetCoefs( iE, theParams, Ecoef, Vcoef ); // edge addition double u = theParams.Coord( myCoordInd[ iE ] ); u = ( 1 - u ) * myFirst[ iE ] + u * myLast[ iE ]; @@ -117,10 +267,32 @@ gp_XY SMESH_Block::TFace::GetUV( const gp_XYZ& theParams ) const gp_XYZ SMESH_Block::TFace::Point( const gp_XYZ& theParams ) const { - gp_XY uv = GetUV( theParams ); - gp_XYZ p = myS->Value( uv.X(), uv.Y() ).XYZ(); - if ( myTrsf.Form() != gp_Identity ) - myTrsf.Transforms( p ); + gp_XYZ p(0.,0.,0.); + if ( !myS ) // if mesh block + { + for ( int iE = 0; iE < 4; iE++ ) // loop on 4 edges + { + double Ecoef = 0, Vcoef = 0; + GetCoefs( iE, theParams, Ecoef, Vcoef ); + // edge addition + double u = theParams.Coord( myCoordInd[ iE ] ); + int i1 = 0, i2 = 1; + switch ( iE ) { + case 1: i1 = 3; i2 = 2; break; + case 2: i1 = 1; i2 = 2; break; + case 3: i1 = 0; i2 = 3; break; + } + p += Ecoef * ( myNodes[ i1 ] * ( 1 - u ) + myNodes[ i2 ] * u ); + // corner addition + p -= Vcoef * myNodes[ iE ]; + } + + } + else // shape block + { + gp_XY uv = GetUV( theParams ); + p = myS->Value( uv.X(), uv.Y() ).XYZ(); + } return p; } @@ -174,155 +346,70 @@ bool SMESH_Block::ShellPoint( const gp_XYZ& theParams, gp_XYZ& thePoint ) const k *= theParams.Coord( iCoef + 1 ); } } - // point on a shape - gp_XYZ Ps; - if ( shapeID < ID_Ex00 ) // vertex - VertexPoint( shapeID, Ps ); - else if ( shapeID < ID_Fxy0 ) { // edge - EdgePoint( shapeID, theParams, Ps ); - k = -k; - } else // face - FacePoint( shapeID, theParams, Ps ); - - thePoint += k * Ps; + // add point on a shape + if ( fabs( k ) > DBL_MIN ) + { + gp_XYZ Ps; + if ( shapeID < ID_Ex00 ) // vertex + VertexPoint( shapeID, Ps ); + else if ( shapeID < ID_Fxy0 ) { // edge + EdgePoint( shapeID, theParams, Ps ); + k = -k; + } else // face + FacePoint( shapeID, theParams, Ps ); + + thePoint += k * Ps; + } } return true; } //======================================================================= -//function : NbVariables -//purpose : -//======================================================================= - -Standard_Integer SMESH_Block::NbVariables() const -{ - return 3; -} - -//======================================================================= -//function : NbEquations -//purpose : +//function : ShellPoint +//purpose : computes coordinates of a point in shell by points on sub-shapes; +// thePointOnShape[ subShapeID ] must be a point on a subShape //======================================================================= -Standard_Integer SMESH_Block::NbEquations() const +bool SMESH_Block::ShellPoint(const gp_XYZ& theParams, + const vector& thePointOnShape, + gp_XYZ& thePoint ) { - return 1; -} + if ( thePointOnShape.size() < ID_F1yz ) + return false; -//======================================================================= -//function : Value -//purpose : -//======================================================================= + double x = theParams.X(), y = theParams.Y(), z = theParams.Z(); + double x1 = 1. - x, y1 = 1. - y, z1 = 1. - z; + 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])); + 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]); -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 ) = myValues[ 0 ]; - } - else { - ShellPoint( params, P ); - gp_Vec dP( P - myPoint ); - theFxyz(1) = SQRT_FUNC ? dP.SquareMagnitude() : dP.Magnitude(); - } return true; } //======================================================================= -//function : Derivatives +//function : Constructor //purpose : //======================================================================= -Standard_Boolean SMESH_Block::Derivatives(const math_Vector& XYZ,math_Matrix& Df) +SMESH_Block::SMESH_Block(): + myNbIterations(0), + mySumDist(0.), + myTolerance(-1.) // to be re-initialized { - MESSAGE( "SMESH_Block::Derivatives()"); - math_Vector F(1,3); - return Values(XYZ,F,Df); -} - -//======================================================================= -//function : Values -//purpose : -//======================================================================= - -Standard_Boolean SMESH_Block::Values(const math_Vector& theXYZ, - math_Vector& theFxyz, - math_Matrix& theDf) -{ -// MESSAGE( endl<<"SMESH_Block::Values( "< DBL_MIN ) - dPi /= mag; - drv[ iP - 1 ] = dPi; - } - for ( int iP = 0; iP < 3; iP++ ) { - if ( iP == myFaceIndex ) - theDf( 1, iP + 1 ) = myFaceParam; - else { - // 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; - } - } - } - //myNbIterations +=3; // how many time call ShellPoint() - - // store better values - myParam = params; - myValues[0]= theFxyz(1); - myValues[1]= theDf(1,1); - myValues[2]= theDf(1,2); - myValues[3]= theDf(1,3); - -// SCRUTE( theFxyz(1) ); -// SCRUTE( theDf( 1,1 )); -// SCRUTE( theDf( 1,2 )); -// SCRUTE( theDf( 1,3 )); - } - - return true; } //======================================================================= @@ -332,22 +419,35 @@ Standard_Boolean SMESH_Block::Values(const math_Vector& theXYZ, bool SMESH_Block::ComputeParameters(const gp_Pnt& thePoint, gp_XYZ& theParams, - const int theShapeID) + const int theShapeID, + const gp_XYZ& theParamsHint) { -// MESSAGE( endl<<"SMESH_Block::ComputeParameters( " -// <FirstParameter(), curve->LastParameter() ); + int i, nb = anExtPC.IsDone() ? anExtPC.NbExt() : 0; + for ( i = 1; i <= nb; i++ ) { + if ( anExtPC.IsMin( i )) + return EdgeParameters( theShapeID, anExtPC.Point( i ).Parameter(), theParams ); + } + return false; + } const bool isOnFace = IsFaceID( theShapeID ); double * coef = GetShapeCoef( theShapeID ); - // the first guess - math_Vector start( 1, 3, 0.0 ); - if ( !myGridComputed ) + // Find the first guess paremeters + + gp_XYZ start(0, 0, 0); + + bool hasHint = ( 0 <= theParamsHint.X() && theParamsHint.X() <= 1 && + 0 <= theParamsHint.Y() && theParamsHint.Y() <= 1 && + 0 <= theParamsHint.Y() && theParamsHint.Y() <= 1 ); + if ( !hasHint && !myGridComputed ) { // define the first guess by thePoint projection on lines // connecting vertices @@ -360,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 ); @@ -368,29 +469,38 @@ bool SMESH_Block::ComputeParameters(const gp_Pnt& thePoint, double par = 0; if ( len2 > zero ) { par = v0P.Dot( v01 ) / len2; - if ( par < 0 || par > 1 ) { + if ( par < 0 || par > 1 ) { // projection falls out of line ends => needGrid needGrid = true; 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 int iNode = 0; + Bnd_Box box; for ( double x = 0.25; x < 0.9; x += 0.25 ) for ( double y = 0.25; y < 0.9; y += 0.25 ) for ( double z = 0.25; z < 0.9; z += 0.25 ) { TxyzPair & prmPtn = my3x3x3GridNodes[ iNode++ ]; prmPtn.first.SetCoord( x, y, z ); ShellPoint( prmPtn.first, prmPtn.second ); + box.Add( gp_Pnt( prmPtn.second )); } myGridComputed = true; + myTolerance = sqrt( box.SquareExtent() ) * 1e-5; } } - if ( myGridComputed ) { + + if ( hasHint ) + { + start = theParamsHint; + } + else if ( myGridComputed ) + { double minDist = DBL_MAX; gp_XYZ* bestParam = 0; for ( int iNode = 0; iNode < 27; iNode++ ) { @@ -401,61 +511,135 @@ 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; } - int 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 ); } } - math_Vector low ( 1, 3, 0.0 ); - math_Vector up ( 1, 3, 1.0 ); - math_Vector tol ( 1, 3, 1e-4 ); - math_FunctionSetRoot paramSearch( *this, tol ); +#ifdef DEBUG_PARAM_COMPUTE + cout << " #### POINT " < 1e-1 && nbLoops++ < 10 ) { - paramSearch.Perform ( *this, start, low, up ); - if ( !paramSearch.IsDone() ) { - //MESSAGE( " !paramSearch.IsDone() " ); + while ( nbLoops <= 100 ) + { + 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; } - else { - //MESSAGE( " NB ITERATIONS: " << paramSearch.NbIterations() ); + 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 ); } - start( 1 ) = myParam.X(); - start( 2 ) = myParam.Y(); - start( 3 ) = myParam.Z(); - //MESSAGE( "Distance: " << ( SQRT_FUNC ? sqrt(myValues[0]) : myValues[0] )); + + nbLoops++; } -// MESSAGE( endl << myParam.X() << " " << myParam.Y() << " " << myParam.Z() << endl); -// mySumDist += myValues[0]; -// MESSAGE( " TOTAL NB ITERATIONS: " << myNbIterations << -// " DIST: " << ( SQRT_FUNC ? sqrt(mySumDist) : mySumDist )); +#ifdef DEBUG_PARAM_COMPUTE + myNbIterations += nbLoops*4; // how many times ShellPoint called + mySumDist += sqrt( sqDistance ); + cout << " ------ SOLUTION: ( "<= 0 ) + theParams.SetCoord( faceIndex + 1, faceParam ); return true; } //======================================================================= -//function : GetStateNumber -//purpose : +//function : VertexParameters +//purpose : return parameters of a vertex given by TShapeID //======================================================================= -Standard_Integer SMESH_Block::GetStateNumber () +bool SMESH_Block::VertexParameters(const int theVertexID, gp_XYZ& theParams) { -// MESSAGE( endl<<"SMESH_Block::GetStateNumber( "< vertexVec; + GetEdgeVertexIDs( theEdgeID, vertexVec ); + VertexParameters( vertexVec[0], theParams ); + TEdge& e = myEdge[ theEdgeID - ID_Ex00 ]; + double param = ( theU - e.EndParam(0) ) / ( e.EndParam(1) - e.EndParam(0) ); + theParams.SetCoord( e.CoordInd(), param ); + return true; + } + return false; } //======================================================================= @@ -548,16 +732,15 @@ int SMESH_Block::GetShapeIDByParams ( const gp_XYZ& theCoord ) return id + 1; // shape ids start at 1 } - //======================================================================= -//function : getOrderedEdges +//function : GetOrderedEdges //purpose : return nb wires and a list of oredered edges //======================================================================= -static int getOrderedEdges (const TopoDS_Face& theFace, - const TopoDS_Vertex& theFirstVertex, - list< TopoDS_Edge >& theEdges, - list< int > & theNbVertexInWires) +int SMESH_Block::GetOrderedEdges (const TopoDS_Face& theFace, + TopoDS_Vertex theFirstVertex, + list< TopoDS_Edge >& theEdges, + list< int > & theNbVertexInWires) { // put wires in a list, so that an outer wire comes first list aWireList; @@ -603,43 +786,220 @@ static int getOrderedEdges (const TopoDS_Face& theFace, isFirst = ( p2.SquareDistance( pf ) < p2.SquareDistance( pl )); if ( isNext ? isFirst : !isFirst ) edge.Reverse(); + // to make a seam go first + if ( theFirstVertex.IsNull() ) + theFirstVertex = TopExp::FirstVertex( edge, true ); } } // rotate theEdges until it begins from theFirstVertex - if ( ! theFirstVertex.IsNull() ) - while ( !theFirstVertex.IsSame( TopExp::FirstVertex( theEdges.front(), true ))) + if ( ! theFirstVertex.IsNull() ) { + TopoDS_Vertex vv[2]; + TopExp::Vertices( theEdges.front(), vv[0], vv[1], true ); + // on closed face, make seam edge the first in the list + while ( !vv[0].IsSame( theFirstVertex ) || vv[0].IsSame( vv[1] )) { theEdges.splice(theEdges.end(), theEdges, - theEdges.begin(), ++ theEdges.begin()); - if ( iE++ > theNbVertexInWires.back() ) + theEdges.begin(), ++theEdges.begin()); + TopExp::Vertices( theEdges.front(), vv[0], vv[1], true ); + if ( iE++ > theNbVertexInWires.back() ) { +#ifdef _DEBUG_ + gp_Pnt p = BRep_Tool::Pnt( theFirstVertex ); + cout << " : Warning : vertex "<< theFirstVertex.TShape().operator->() + << " ( " << p.X() << " " << p.Y() << " " << p.Z() << " )" + << " not found in outer wire of face "<< theFace.TShape().operator->() + << " with vertices: " << endl; + wExp.Init( *wlIt, theFace ); + for ( int i = 0; wExp.More(); wExp.Next(), i++ ) + { + TopoDS_Edge edge = wExp.Current(); + edge = TopoDS::Edge( edge.Oriented( wExp.Orientation() )); + TopoDS_Vertex v = TopExp::FirstVertex( edge, true ); + gp_Pnt p = BRep_Tool::Pnt( v ); + cout << i << " " << v.TShape().operator->() << " " + << p.X() << " " << p.Y() << " " << p.Z() << " " << endl; + } +#endif break; // break infinite loop + } } - } + } + } // end outer wire } return aWireList.size(); } +//================================================================================ +/*! + * \brief Call it after geometry initialisation + */ +//================================================================================ + +void SMESH_Block::init() +{ + myNbIterations = 0; + mySumDist = 0; + myGridComputed = false; +} + +//======================================================================= +//function : LoadMeshBlock +//purpose : prepare to work with theVolume +//======================================================================= + +#define gpXYZ(n) gp_XYZ(n->X(),n->Y(),n->Z()) + +bool SMESH_Block::LoadMeshBlock(const SMDS_MeshVolume* theVolume, + const int theNode000Index, + const int theNode001Index, + vector& theOrderedNodes) +{ + MESSAGE(" ::LoadMeshBlock()"); + init(); + + SMDS_VolumeTool vTool; + if (!vTool.Set( theVolume ) || vTool.NbNodes() != 8 || + !vTool.IsLinked( theNode000Index, theNode001Index )) { + MESSAGE(" Bad arguments "); + return false; + } + vTool.SetExternalNormal(); + // In terms of indices used for access to nodes and faces in SMDS_VolumeTool: + int V000, V100, V010, V110, V001, V101, V011, V111; // 8 vertices + int Fxy0, Fxy1; // bottom and top faces + // vertices of faces + vector vFxy0, vFxy1; + + V000 = theNode000Index; + V001 = theNode001Index; + + // get faces sharing V000 and V001 + list fV000, fV001; + int i, iF, iE, iN; + for ( iF = 0; iF < vTool.NbFaces(); ++iF ) { + const int* nid = vTool.GetFaceNodesIndices( iF ); + for ( iN = 0; iN < 4; ++iN ) + if ( nid[ iN ] == V000 ) { + fV000.push_back( iF ); + } else if ( nid[ iN ] == V001 ) { + fV001.push_back( iF ); + } + } + + // find the bottom (Fxy0), the top (Fxy1) faces + list::iterator fIt1, fIt2, Fxy0Pos; + for ( fIt1 = fV000.begin(); fIt1 != fV000.end(); fIt1++) { + fIt2 = std::find( fV001.begin(), fV001.end(), *fIt1 ); + if ( fIt2 != fV001.end() ) { // *fIt1 is in the both lists + fV001.erase( fIt2 ); // erase Fx0z or F0yz from fV001 + } else { // *fIt1 is in fV000 only + Fxy0Pos = fIt1; // points to Fxy0 + } + } + Fxy0 = *Fxy0Pos; + Fxy1 = fV001.front(); + const SMDS_MeshNode** nn = vTool.GetNodes(); + + // find bottom veritices, their order is that a face normal is external + vFxy0.resize(4); + const int* nid = vTool.GetFaceNodesIndices( Fxy0 ); + for ( i = 0; i < 4; ++i ) + if ( nid[ i ] == V000 ) + break; + for ( iN = 0; iN < 4; ++iN, ++i ) { + if ( i == 4 ) i = 0; + vFxy0[ iN ] = nid[ i ]; + } + // find top veritices, their order is that a face normal is external + vFxy1.resize(4); + nid = vTool.GetFaceNodesIndices( Fxy1 ); + for ( i = 0; i < 4; ++i ) + if ( nid[ i ] == V001 ) + break; + for ( iN = 0; iN < 4; ++iN, ++i ) { + if ( i == 4 ) i = 0; + vFxy1[ iN ] = nid[ i ]; + } + // find indices of the rest veritices + V100 = vFxy0[3]; + V010 = vFxy0[1]; + V110 = vFxy0[2]; + V101 = vFxy1[1]; + V011 = vFxy1[3]; + V111 = vFxy1[2]; + + // set points coordinates + myPnt[ ID_V000 - 1 ] = gpXYZ( nn[ V000 ] ); + myPnt[ ID_V100 - 1 ] = gpXYZ( nn[ V100 ] ); + myPnt[ ID_V010 - 1 ] = gpXYZ( nn[ V010 ] ); + myPnt[ ID_V110 - 1 ] = gpXYZ( nn[ V110 ] ); + myPnt[ ID_V001 - 1 ] = gpXYZ( nn[ V001 ] ); + myPnt[ ID_V101 - 1 ] = gpXYZ( nn[ V101 ] ); + myPnt[ ID_V011 - 1 ] = gpXYZ( nn[ V011 ] ); + myPnt[ ID_V111 - 1 ] = gpXYZ( nn[ V111 ] ); + + // fill theOrderedNodes + theOrderedNodes.resize( 8 ); + theOrderedNodes[ 0 ] = nn[ V000 ]; + theOrderedNodes[ 1 ] = nn[ V100 ]; + theOrderedNodes[ 2 ] = nn[ V010 ]; + theOrderedNodes[ 3 ] = nn[ V110 ]; + theOrderedNodes[ 4 ] = nn[ V001 ]; + theOrderedNodes[ 5 ] = nn[ V101 ]; + theOrderedNodes[ 6 ] = nn[ V011 ]; + theOrderedNodes[ 7 ] = nn[ V111 ]; + + // fill edges + vector< int > vertexVec; + for ( iE = 0; iE < NbEdges(); ++iE ) { + GetEdgeVertexIDs(( iE + ID_FirstE ), vertexVec ); + myEdge[ iE ].Set(( iE + ID_FirstE ), + myPnt[ vertexVec[0] - 1 ], + myPnt[ vertexVec[1] - 1 ]); + } + + // fill faces' corners + for ( iF = ID_Fxy0; iF < ID_Shell; ++iF ) + { + TFace& tFace = myFace[ iF - ID_FirstF ]; + vector< int > edgeIdVec(4, -1); + GetFaceEdgesIDs( iF, edgeIdVec ); + tFace.Set( iF, myEdge[ edgeIdVec [ 0 ] - ID_Ex00], myEdge[ edgeIdVec [ 1 ] - ID_Ex00]); + } + + return true; +} //======================================================================= //function : LoadBlockShapes -//purpose : add sub-shapes of theBlock to theShapeIDMap so that they get +//purpose : Initialize block geometry with theShell, +// add sub-shapes of theBlock to theShapeIDMap so that they get // IDs acoording to enum TShapeID //======================================================================= bool SMESH_Block::LoadBlockShapes(const TopoDS_Shell& theShell, const TopoDS_Vertex& theVertex000, const TopoDS_Vertex& theVertex001, -// TopTools_IndexedMapOfShape& theShapeIDMap TopTools_IndexedMapOfOrientedShape& theShapeIDMap ) { MESSAGE(" ::LoadBlockShapes()"); + return ( FindBlockShapes( theShell, theVertex000, theVertex001, theShapeIDMap ) && + LoadBlockShapes( theShapeIDMap )); +} - myShell = theShell; - myNbIterations = 0; - mySumDist = 0; - myGridComputed = false; +//======================================================================= +//function : LoadBlockShapes +//purpose : add sub-shapes of theBlock to theShapeIDMap so that they get +// IDs acoording to enum TShapeID +//======================================================================= - // 6 vertices +bool SMESH_Block::FindBlockShapes(const TopoDS_Shell& theShell, + const TopoDS_Vertex& theVertex000, + const TopoDS_Vertex& theVertex001, + TopTools_IndexedMapOfOrientedShape& theShapeIDMap ) +{ + MESSAGE(" ::FindBlockShapes()"); + + // 8 vertices TopoDS_Shape V000, V100, V010, V110, V001, V101, V011, V111; // 12 edges TopoDS_Shape Ex00, Ex10, Ex01, Ex11; @@ -653,7 +1013,7 @@ bool SMESH_Block::LoadBlockShapes(const TopoDS_Shell& theShell, const int NB_FACES_BY_VERTEX = 6; TopTools_IndexedDataMapOfShapeListOfShape vfMap; - TopExp::MapShapesAndAncestors( myShell, TopAbs_VERTEX, TopAbs_FACE, vfMap ); + TopExp::MapShapesAndAncestors( theShell, TopAbs_VERTEX, TopAbs_FACE, vfMap ); if ( vfMap.Extent() != 8 ) { MESSAGE(" Wrong nb of vertices in the block: " << vfMap.Extent() ); return false; @@ -677,7 +1037,7 @@ bool SMESH_Block::LoadBlockShapes(const TopoDS_Shell& theShell, } // find vertex 001 - the one on the most vertical edge passing through V000 TopTools_IndexedDataMapOfShapeListOfShape veMap; - TopExp::MapShapesAndAncestors( myShell, TopAbs_VERTEX, TopAbs_EDGE, veMap ); + TopExp::MapShapesAndAncestors( theShell, TopAbs_VERTEX, TopAbs_EDGE, veMap ); gp_Vec dir001 = gp::DZ(); gp_Pnt p000 = BRep_Tool::Pnt( TopoDS::Vertex( V000 )); double maxVal = -DBL_MAX; @@ -746,7 +1106,7 @@ bool SMESH_Block::LoadBlockShapes(const TopoDS_Shell& theShell, // find bottom edges and veritices list< TopoDS_Edge > eList; list< int > nbVertexInWires; - getOrderedEdges( TopoDS::Face( Fxy0 ), TopoDS::Vertex( V000 ), eList, nbVertexInWires ); + GetOrderedEdges( TopoDS::Face( Fxy0 ), TopoDS::Vertex( V000 ), eList, nbVertexInWires ); if ( nbVertexInWires.size() != 1 || nbVertexInWires.front() != 4 ) { MESSAGE(" LoadBlockShapes() error "); return false; @@ -768,7 +1128,7 @@ bool SMESH_Block::LoadBlockShapes(const TopoDS_Shell& theShell, // find top edges and veritices eList.clear(); - getOrderedEdges( TopoDS::Face( Fxy1 ), TopoDS::Vertex( V001 ), eList, nbVertexInWires ); + GetOrderedEdges( TopoDS::Face( Fxy1 ), TopoDS::Vertex( V001 ), eList, nbVertexInWires ); if ( nbVertexInWires.size() != 1 || nbVertexInWires.front() != 4 ) { MESSAGE(" LoadBlockShapes() error "); return false; @@ -792,7 +1152,7 @@ bool SMESH_Block::LoadBlockShapes(const TopoDS_Shell& theShell, if ( V101.IsSame( exp.Current() ) || V100.IsSame( exp.Current() )) break; // V101 or V100 found if ( !exp.More() ) { // not found - TopoDS_Shape f = Fx0z; Fx0z = F0yz; F0yz = f; + std::swap( Fx0z, F0yz); } // find Fx1z and F1yz faces @@ -827,7 +1187,7 @@ bool SMESH_Block::LoadBlockShapes(const TopoDS_Shell& theShell, if ( V010.IsSame( exp.Current() ) || V011.IsSame( exp.Current() )) break; if ( !exp.More() ) { - TopoDS_Shape f = Fx1z; Fx1z = F1yz; F1yz = f; + std::swap( Fx1z, F1yz); } // find vertical edges @@ -891,12 +1251,22 @@ bool SMESH_Block::LoadBlockShapes(const TopoDS_Shell& theShell, theShapeIDMap.Add(F0yz); theShapeIDMap.Add(F1yz); - theShapeIDMap.Add(myShell); + theShapeIDMap.Add(theShell); - if ( theShapeIDMap.Extent() != 27 ) { - MESSAGE("LoadBlockShapes() " << theShapeIDMap.Extent() ); - return false; - } + return true; +} + +//================================================================================ +/*! + * \brief Initialize block geometry with shapes from theShapeIDMap + * \param theShapeIDMap - map of block subshapes + * \retval bool - is a success + */ +//================================================================================ + +bool SMESH_Block::LoadBlockShapes(const TopTools_IndexedMapOfOrientedShape& theShapeIDMap) +{ + init(); // store shapes geometry for ( int shapeID = 1; shapeID < theShapeIDMap.Extent(); shapeID++ ) @@ -906,59 +1276,24 @@ bool SMESH_Block::LoadBlockShapes(const TopoDS_Shell& theShell, { case TopAbs_VERTEX: { - if ( shapeID > ID_V111 ) { - MESSAGE(" shapeID =" << shapeID ); - return false; - } - myPnt[ shapeID - ID_V000 ] = - BRep_Tool::Pnt( TopoDS::Vertex( S )).XYZ(); + if ( !IsVertexID( ID_V111 )) return false; + myPnt[ shapeID - ID_V000 ] = BRep_Tool::Pnt( TopoDS::Vertex( S )).XYZ(); break; } case TopAbs_EDGE: { + if ( !IsEdgeID( shapeID )) return false; const TopoDS_Edge& edge = TopoDS::Edge( S ); - if ( shapeID < ID_Ex00 || shapeID > ID_E11z || edge.IsNull() ) { - MESSAGE(" shapeID =" << shapeID ); - return false; - } - TEdge& tEdge = myEdge[ shapeID - ID_Ex00 ]; - tEdge.myCoordInd = GetCoordIndOnEdge( shapeID ); - TopLoc_Location loc; - tEdge.myC3d = BRep_Tool::Curve( edge, loc, tEdge.myFirst, tEdge.myLast ); - if ( !IsForwardEdge( edge, theShapeIDMap )) - Swap( tEdge.myFirst, tEdge.myLast ); - tEdge.myTrsf = loc; + TEdge& tEdge = myEdge[ shapeID - ID_FirstE ]; + tEdge.Set( shapeID, + new BRepAdaptor_Curve( edge ), + IsForwardEdge( edge, theShapeIDMap )); break; } case TopAbs_FACE: { - const TopoDS_Face& face = TopoDS::Face( S ); - if ( shapeID < ID_Fxy0 || shapeID > ID_F1yz || face.IsNull() ) { - MESSAGE(" shapeID =" << shapeID ); + if ( !LoadFace( TopoDS::Face( S ), shapeID, theShapeIDMap )) return false; - } - TFace& tFace = myFace[ shapeID - ID_Fxy0 ]; - // pcurves - vector< int > edgeIdVec(4, -1); - GetFaceEdgesIDs( shapeID, edgeIdVec ); - for ( int iE = 0; iE < 4; iE++ ) // loop on 4 edges - { - const TopoDS_Edge& edge = TopoDS::Edge( theShapeIDMap( edgeIdVec[ iE ])); - tFace.myCoordInd[ iE ] = GetCoordIndOnEdge( edgeIdVec[ iE ] ); - tFace.myC2d[ iE ] = - BRep_Tool::CurveOnSurface( edge, face, tFace.myFirst[iE], tFace.myLast[iE] ); - if ( !IsForwardEdge( edge, theShapeIDMap )) - Swap( tFace.myFirst[ iE ], tFace.myLast[ iE ] ); - } - // 2d corners - tFace.myCorner[ 0 ] = tFace.myC2d[ 0 ]->Value( tFace.myFirst[0] ).XY(); - tFace.myCorner[ 1 ] = tFace.myC2d[ 0 ]->Value( tFace.myLast[0] ).XY(); - tFace.myCorner[ 2 ] = tFace.myC2d[ 1 ]->Value( tFace.myLast[1] ).XY(); - tFace.myCorner[ 3 ] = tFace.myC2d[ 1 ]->Value( tFace.myFirst[1] ).XY(); - // sufrace - TopLoc_Location loc; - tFace.myS = BRep_Tool::Surface( face, loc ); - tFace.myTrsf = loc; break; } default: break; @@ -968,6 +1303,76 @@ bool SMESH_Block::LoadBlockShapes(const TopoDS_Shell& theShell, return true; } +//================================================================================ +/*! + * \brief Load face geometry + * \param theFace - face + * \param theFaceID - face in-block ID + * \param theShapeIDMap - map of block subshapes + * \retval bool - is a success + * + * It is enough to compute params or coordinates on the face. + * Face subshapes must be loaded into theShapeIDMap before + */ +//================================================================================ + +bool SMESH_Block::LoadFace(const TopoDS_Face& theFace, + const int theFaceID, + const TopTools_IndexedMapOfOrientedShape& theShapeIDMap) +{ + if ( !IsFaceID( theFaceID ) ) return false; + // pcurves + Adaptor2d_Curve2d* c2d[4]; + bool isForward[4]; + vector< int > edgeIdVec; + GetFaceEdgesIDs( theFaceID, edgeIdVec ); + for ( int iE = 0; iE < edgeIdVec.size(); iE++ ) // loop on 4 edges + { + if ( edgeIdVec[ iE ] > theShapeIDMap.Extent() ) + return false; + const TopoDS_Edge& edge = TopoDS::Edge( theShapeIDMap( edgeIdVec[ iE ])); + c2d[ iE ] = new BRepAdaptor_Curve2d( edge, theFace ); + isForward[ iE ] = IsForwardEdge( edge, theShapeIDMap ); + } + TFace& tFace = myFace[ theFaceID - ID_FirstF ]; + tFace.Set( theFaceID, new BRepAdaptor_Surface( theFace ), c2d, isForward ); + return true; +} + +//================================================================================ +/*! + * \brief/ Insert theShape into theShapeIDMap with theShapeID + * \param theShape - shape to insert + * \param theShapeID - shape in-block ID + * \param theShapeIDMap - map of block subshapes + */ +//================================================================================ + +bool SMESH_Block::Insert(const TopoDS_Shape& theShape, + const int theShapeID, + TopTools_IndexedMapOfOrientedShape& theShapeIDMap) +{ + if ( !theShape.IsNull() && theShapeID > 0 ) + { + if ( theShapeIDMap.Contains( theShape )) + return ( theShapeIDMap.FindIndex( theShape ) == theShapeID ); + + if ( theShapeID <= theShapeIDMap.Extent() ) { + theShapeIDMap.Substitute( theShapeID, theShape ); + } + else { + while ( theShapeIDMap.Extent() < theShapeID - 1 ) { + TopoDS_Compound comp; + BRep_Builder().MakeCompound( comp ); + theShapeIDMap.Add( comp ); + } + theShapeIDMap.Add( theShape ); + } + return true; + } + return false; +} + //======================================================================= //function : GetFaceEdgesIDs //purpose : return edges IDs in the order u0, u1, 0v, 1v @@ -976,6 +1381,7 @@ bool SMESH_Block::LoadBlockShapes(const TopoDS_Shell& theShell, void SMESH_Block::GetFaceEdgesIDs (const int faceID, vector< int >& edgeVec ) { + edgeVec.resize( 4 ); switch ( faceID ) { case ID_Fxy0: edgeVec[ 0 ] = ID_Ex00; @@ -1017,3 +1423,70 @@ void SMESH_Block::GetFaceEdgesIDs (const int faceID, vector< int >& edgeVec ) MESSAGE(" GetFaceEdgesIDs(), wrong face ID: " << faceID ); } } + +//======================================================================= +//function : GetEdgeVertexIDs +//purpose : return vertex IDs of an edge +//======================================================================= + +void SMESH_Block::GetEdgeVertexIDs (const int edgeID, vector< int >& vertexVec ) +{ + vertexVec.resize( 2 ); + switch ( edgeID ) { + + case ID_Ex00: + vertexVec[ 0 ] = ID_V000; + vertexVec[ 1 ] = ID_V100; + break; + case ID_Ex10: + vertexVec[ 0 ] = ID_V010; + vertexVec[ 1 ] = ID_V110; + break; + case ID_Ex01: + vertexVec[ 0 ] = ID_V001; + vertexVec[ 1 ] = ID_V101; + break; + case ID_Ex11: + vertexVec[ 0 ] = ID_V011; + vertexVec[ 1 ] = ID_V111; + break; + + case ID_E0y0: + vertexVec[ 0 ] = ID_V000; + vertexVec[ 1 ] = ID_V010; + break; + case ID_E1y0: + vertexVec[ 0 ] = ID_V100; + vertexVec[ 1 ] = ID_V110; + break; + case ID_E0y1: + vertexVec[ 0 ] = ID_V001; + vertexVec[ 1 ] = ID_V011; + break; + case ID_E1y1: + vertexVec[ 0 ] = ID_V101; + vertexVec[ 1 ] = ID_V111; + break; + + case ID_E00z: + vertexVec[ 0 ] = ID_V000; + vertexVec[ 1 ] = ID_V001; + break; + case ID_E10z: + vertexVec[ 0 ] = ID_V100; + vertexVec[ 1 ] = ID_V101; + break; + case ID_E01z: + vertexVec[ 0 ] = ID_V010; + vertexVec[ 1 ] = ID_V011; + break; + case ID_E11z: + vertexVec[ 0 ] = ID_V110; + vertexVec[ 1 ] = ID_V111; + break; + default: + vertexVec.resize(0); + MESSAGE(" GetEdgeVertexIDs(), wrong edge ID: " << edgeID ); + } +} +