src/SMESH/SMESH_Block.cxx

   1 //  Copyright (C) 2003  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   2 //  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   3 //
   4 //  This library is free software; you can redistribute it and/or
   5 //  modify it under the terms of the GNU Lesser General Public
   6 //  License as published by the Free Software Foundation; either
   7 //  version 2.1 of the License.
   8 //
   9 //  This library is distributed in the hope that it will be useful,
  10 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12 //  Lesser General Public License for more details.
  13 //
  14 //  You should have received a copy of the GNU Lesser General Public
  15 //  License along with this library; if not, write to the Free Software
  16 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  17 //
  18 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  19
  20 // File      : SMESH_Pattern.hxx
  21 // Created   : Mon Aug  2 10:30:00 2004
  22 // Author    : Edward AGAPOV (eap)
  23
  24 #include "SMESH_Block.hxx"
  25
  26 #include <BRepAdaptor_Curve.hxx>
  27 #include <BRepAdaptor_Curve2d.hxx>
  28 #include <BRepAdaptor_Surface.hxx>
  29 #include <BRepTools.hxx>
  30 #include <BRepTools_WireExplorer.hxx>
  31 #include <BRep_Builder.hxx>
  32 #include <BRep_Tool.hxx>
  33 #include <Bnd_Box.hxx>
  34 #include <Extrema_ExtPC.hxx>
  35 #include <Extrema_POnCurv.hxx>
  36 #include <Geom2d_Curve.hxx>
  37 #include <TopExp_Explorer.hxx>
  38 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
  39 #include <TopTools_ListIteratorOfListOfShape.hxx>
  40 #include <TopTools_ListOfShape.hxx>
  41 #include <TopoDS.hxx>
  42 #include <TopoDS_Compound.hxx>
  43 #include <TopoDS_Face.hxx>
  44 #include <TopoDS_Iterator.hxx>
  45 #include <TopoDS_Wire.hxx>
  46 #include <gp_Trsf.hxx>
  47 #include <gp_Vec.hxx>
  48 #include <math_FunctionSetRoot.hxx>
  49 #include <math_Matrix.hxx>
  50 #include <math_Vector.hxx>
  51
  52 #include "SMDS_MeshNode.hxx"
  53 #include "SMDS_MeshVolume.hxx"
  54 #include "SMDS_VolumeTool.hxx"
  55 #include "utilities.h"
  56
  57 #include <list>
  58
  59 using namespace std;
  60
  61 //#define DEBUG_PARAM_COMPUTE
  62
  63 //================================================================================
  64 /*!
  65  * \brief Set edge data
  66   * \param edgeID - block subshape ID
  67   * \param curve - edge geometry
  68   * \param isForward - is curve orientation coincides with edge orientation in the block
  69  */
  70 //================================================================================
  71
  72 void SMESH_Block::TEdge::Set( const int edgeID, Adaptor3d_Curve* curve, const bool isForward )
  73 {
  74   myCoordInd = SMESH_Block::GetCoordIndOnEdge( edgeID );
  75   if ( myC3d ) delete myC3d;
  76   myC3d = curve;
  77   myFirst = curve->FirstParameter();
  78   myLast = curve->LastParameter();
  79   if ( !isForward )
  80     std::swap( myFirst, myLast );
  81 }
  82
  83 //================================================================================
  84 /*!
  85  * \brief Set coordinates of nodes at edge ends to work with mesh block
  86   * \param edgeID - block subshape ID
  87   * \param node1 - coordinates of node with lower ID
  88   * \param node2 - coordinates of node with upper ID
  89  */
  90 //================================================================================
  91
  92 void SMESH_Block::TEdge::Set( const int edgeID, const gp_XYZ& node1, const gp_XYZ& node2 )
  93 {
  94   myCoordInd = SMESH_Block::GetCoordIndOnEdge( edgeID );
  95   myNodes[ 0 ] = node1;
  96   myNodes[ 1 ] = node2;
  97
  98   if ( myC3d ) delete myC3d;
  99   myC3d = 0;
 100 }
 101
 102 //=======================================================================
 103 //function : SMESH_Block::TEdge::GetU
 104 //purpose  :
 105 //=======================================================================
 106
 107 double SMESH_Block::TEdge::GetU( const gp_XYZ& theParams ) const
 108 {
 109   double u = theParams.Coord( myCoordInd );
 110   if ( !myC3d ) // if mesh block
 111     return u;
 112   return ( 1 - u ) * myFirst + u * myLast;
 113 }
 114
 115 //=======================================================================
 116 //function : SMESH_Block::TEdge::Point
 117 //purpose  :
 118 //=======================================================================
 119
 120 gp_XYZ SMESH_Block::TEdge::Point( const gp_XYZ& theParams ) const
 121 {
 122   double u = GetU( theParams );
 123   if ( myC3d ) return myC3d->Value( u ).XYZ();
 124   // mesh block
 125   return myNodes[0] * ( 1 - u ) + myNodes[1] * u;
 126 }
 127
 128 //================================================================================
 129 /*!
 130  * \brief Destructor
 131  */
 132 //================================================================================
 133
 134 SMESH_Block::TEdge::~TEdge()
 135 {
 136   if ( myC3d ) delete myC3d;
 137 }
 138
 139 //================================================================================
 140 /*!
 141  * \brief Set face data
 142   * \param faceID - block subshape ID
 143   * \param S - face surface geometry
 144   * \param c2d - 4 pcurves in the order as returned by GetFaceEdgesIDs(faceID)
 145   * \param isForward - orientation of pcurves comparing with block edge direction
 146  */
 147 //================================================================================
 148
 149 void SMESH_Block::TFace::Set( const int          faceID,
 150                               Adaptor3d_Surface* S,
 151                               Adaptor2d_Curve2d* c2D[4],
 152                               const bool         isForward[4] )
 153 {
 154   if ( myS ) delete myS;
 155   myS = S;
 156   // pcurves
 157   vector< int > edgeIdVec;
 158   GetFaceEdgesIDs( faceID, edgeIdVec );
 159   for ( int iE = 0; iE < edgeIdVec.size(); iE++ ) // loop on 4 edges
 160   {
 161     myCoordInd[ iE ] = GetCoordIndOnEdge( edgeIdVec[ iE ] );
 162     if ( myC2d[ iE ]) delete myC2d[ iE ];
 163     myC2d[ iE ] = c2D[ iE ];
 164     myFirst[ iE ] = myC2d[ iE ]->FirstParameter();
 165     myLast [ iE ] = myC2d[ iE ]->LastParameter();
 166     if ( !isForward[ iE ])
 167       std::swap( myFirst[ iE ], myLast[ iE ] );
 168   }
 169   // 2d corners
 170   myCorner[ 0 ] = myC2d[ 0 ]->Value( myFirst[0] ).XY();
 171   myCorner[ 1 ] = myC2d[ 0 ]->Value( myLast[0] ).XY();
 172   myCorner[ 2 ] = myC2d[ 1 ]->Value( myLast[1] ).XY();
 173   myCorner[ 3 ] = myC2d[ 1 ]->Value( myFirst[1] ).XY();
 174 }
 175
 176 //================================================================================
 177 /*!
 178  * \brief Set face data to work with mesh block
 179   * \param faceID - block subshape ID
 180   * \param edgeU0 - filled data of edge u0 = GetFaceEdgesIDs(faceID)[ 0 ]
 181   * \param edgeU1 - filled data of edge u1 = GetFaceEdgesIDs(faceID)[ 1 ]
 182  */
 183 //================================================================================
 184
 185 void SMESH_Block::TFace::Set( const int faceID, const TEdge& edgeU0, const TEdge& edgeU1 )
 186 {
 187   vector< int > edgeIdVec;
 188   GetFaceEdgesIDs( faceID, edgeIdVec );
 189   myNodes[ 0 ] = edgeU0.NodeXYZ( 1 );
 190   myNodes[ 1 ] = edgeU0.NodeXYZ( 0 );
 191   myNodes[ 2 ] = edgeU1.NodeXYZ( 0 );
 192   myNodes[ 3 ] = edgeU1.NodeXYZ( 1 );
 193   myCoordInd[ 0 ] = GetCoordIndOnEdge( edgeIdVec[ 0 ] );
 194   myCoordInd[ 1 ] = GetCoordIndOnEdge( edgeIdVec[ 1 ] );
 195   myCoordInd[ 2 ] = GetCoordIndOnEdge( edgeIdVec[ 2 ] );
 196   myCoordInd[ 3 ] = GetCoordIndOnEdge( edgeIdVec[ 3 ] );
 197   if ( myS ) delete myS;
 198   myS = 0;
 199 }
 200
 201 //================================================================================
 202 /*!
 203  * \brief Destructor
 204  */
 205 //================================================================================
 206
 207 SMESH_Block::TFace::~TFace()
 208 {
 209   if ( myS ) delete myS;
 210   for ( int i = 0 ; i < 4; ++i )
 211     if ( myC2d[ i ]) delete myC2d[ i ];
 212 }
 213
 214 //=======================================================================
 215 //function : SMESH_Block::TFace::GetCoefs
 216 //purpose  : return coefficients for addition of [0-3]-th edge and vertex
 217 //=======================================================================
 218
 219 void SMESH_Block::TFace::GetCoefs(int           iE,
 220                                   const gp_XYZ& theParams,
 221                                   double&       Ecoef,
 222                                   double&       Vcoef ) const
 223 {
 224   double dU = theParams.Coord( GetUInd() );
 225   double dV = theParams.Coord( GetVInd() );
 226   switch ( iE ) {
 227   case 0:
 228     Ecoef = ( 1 - dV ); // u0
 229     Vcoef = ( 1 - dU ) * ( 1 - dV ); break; // 00
 230   case 1:
 231     Ecoef = dV; // u1
 232     Vcoef = dU * ( 1 - dV ); break; // 10
 233   case 2:
 234     Ecoef = ( 1 - dU ); // 0v
 235     Vcoef = dU * dV  ; break; // 11
 236   case 3:
 237     Ecoef = dU  ; // 1v
 238     Vcoef = ( 1 - dU ) * dV  ; break; // 01
 239   default: ASSERT(0);
 240   }
 241 }
 242
 243 //=======================================================================
 244 //function : SMESH_Block::TFace::GetUV
 245 //purpose  :
 246 //=======================================================================
 247
 248 gp_XY SMESH_Block::TFace::GetUV( const gp_XYZ& theParams ) const
 249 {
 250   gp_XY uv(0.,0.);
 251   for ( int iE = 0; iE < 4; iE++ ) // loop on 4 edges
 252   {
 253     double Ecoef = 0, Vcoef = 0;
 254     GetCoefs( iE, theParams, Ecoef, Vcoef );
 255     // edge addition
 256     double u = theParams.Coord( myCoordInd[ iE ] );
 257     u = ( 1 - u ) * myFirst[ iE ] + u * myLast[ iE ];
 258     uv += Ecoef * myC2d[ iE ]->Value( u ).XY();
 259     // corner addition
 260     uv -= Vcoef * myCorner[ iE ];
 261   }
 262   return uv;
 263 }
 264
 265 //=======================================================================
 266 //function : SMESH_Block::TFace::Point
 267 //purpose  :
 268 //=======================================================================
 269
 270 gp_XYZ SMESH_Block::TFace::Point( const gp_XYZ& theParams ) const
 271 {
 272   gp_XYZ p(0.,0.,0.);
 273   if ( !myS ) // if mesh block
 274   {
 275     for ( int iE = 0; iE < 4; iE++ ) // loop on 4 edges
 276     {
 277       double Ecoef = 0, Vcoef = 0;
 278       GetCoefs( iE, theParams, Ecoef, Vcoef );
 279       // edge addition
 280       double u = theParams.Coord( myCoordInd[ iE ] );
 281       int i1 = 0, i2 = 1;
 282       switch ( iE ) {
 283       case 1: i1 = 3; i2 = 2; break;
 284       case 2: i1 = 1; i2 = 2; break;
 285       case 3: i1 = 0; i2 = 3; break;
 286       }
 287       p += Ecoef * ( myNodes[ i1 ] * ( 1 - u ) + myNodes[ i2 ] * u );
 288       // corner addition
 289       p -= Vcoef * myNodes[ iE ];
 290     }
 291
 292   }
 293   else // shape block
 294   {
 295     gp_XY uv = GetUV( theParams );
 296     p = myS->Value( uv.X(), uv.Y() ).XYZ();
 297   }
 298   return p;
 299 }
 300
 301 //=======================================================================
 302 //function : GetShapeCoef
 303 //purpose  :
 304 //=======================================================================
 305
 306 double* SMESH_Block::GetShapeCoef (const int theShapeID)
 307 {
 308   static double shapeCoef[][3] = {
 309     //    V000,        V100,        V010,         V110
 310     { -1,-1,-1 }, {  1,-1,-1 }, { -1, 1,-1 }, {  1, 1,-1 },
 311     //    V001,        V101,        V011,         V111,
 312     { -1,-1, 1 }, {  1,-1, 1 }, { -1, 1, 1 }, {  1, 1, 1 },
 313     //    Ex00,        Ex10,        Ex01,         Ex11,
 314     {  0,-1,-1 }, {  0, 1,-1 }, {  0,-1, 1 }, {  0, 1, 1 },
 315     //    E0y0,        E1y0,        E0y1,         E1y1,
 316     { -1, 0,-1 }, {  1, 0,-1 }, { -1, 0, 1 }, {  1, 0, 1 },
 317     //    E00z,        E10z,        E01z,         E11z,
 318     { -1,-1, 0 }, {  1,-1, 0 }, { -1, 1, 0 }, {  1, 1, 0 },
 319     //    Fxy0,        Fxy1,        Fx0z,         Fx1z,         F0yz,           F1yz,
 320     {  0, 0,-1 }, {  0, 0, 1 }, {  0,-1, 0 }, {  0, 1, 0 }, { -1, 0, 0 }, {  1, 0, 0 },
 321     // ID_Shell
 322     {  0, 0, 0 }
 323   };
 324   if ( theShapeID < ID_V000 || theShapeID > ID_F1yz )
 325     return shapeCoef[ ID_Shell - 1 ];
 326
 327   return shapeCoef[ theShapeID - 1 ];
 328 }
 329
 330 //=======================================================================
 331 //function : ShellPoint
 332 //purpose  : return coordinates of a point in shell
 333 //=======================================================================
 334
 335 bool SMESH_Block::ShellPoint( const gp_XYZ& theParams, gp_XYZ& thePoint ) const
 336 {
 337   thePoint.SetCoord( 0., 0., 0. );
 338   for ( int shapeID = ID_V000; shapeID < ID_Shell; shapeID++ )
 339   {
 340     // coef
 341     double* coefs = GetShapeCoef( shapeID );
 342     double k = 1;
 343     for ( int iCoef = 0; iCoef < 3; iCoef++ ) {
 344       if ( coefs[ iCoef ] != 0 ) {
 345         if ( coefs[ iCoef ] < 0 )
 346           k *= ( 1. - theParams.Coord( iCoef + 1 ));
 347         else
 348           k *= theParams.Coord( iCoef + 1 );
 349       }
 350     }
 351     // add point on a shape
 352     if ( fabs( k ) > DBL_MIN )
 353     {
 354       gp_XYZ Ps;
 355       if ( shapeID < ID_Ex00 ) // vertex
 356         VertexPoint( shapeID, Ps );
 357       else if ( shapeID < ID_Fxy0 ) { // edge
 358         EdgePoint( shapeID, theParams, Ps );
 359         k = -k;
 360       } else // face
 361         FacePoint( shapeID, theParams, Ps );
 362
 363       thePoint += k * Ps;
 364     }
 365   }
 366   return true;
 367 }
 368
 369 //=======================================================================
 370 //function : ShellPoint
 371 //purpose  : computes coordinates of a point in shell by points on sub-shapes;
 372 //           thePointOnShape[ subShapeID ] must be a point on a subShape
 373 //=======================================================================
 374
 375 bool SMESH_Block::ShellPoint(const gp_XYZ&         theParams,
 376                              const vector<gp_XYZ>& thePointOnShape,
 377                              gp_XYZ&               thePoint )
 378 {
 379   if ( thePointOnShape.size() < ID_F1yz )
 380     return false;
 381
 382   double x = theParams.X(), y = theParams.Y(), z = theParams.Z();
 383   double x1 = 1. - x,       y1 = 1. - y,       z1 = 1. - z;
 384   const vector<gp_XYZ>& p = thePointOnShape;
 385
 386   thePoint =
 387     x1 * p[ID_F0yz] + x * p[ID_F1yz] +
 388     y1 * p[ID_Fx0z] + y * p[ID_Fx1z] +
 389     z1 * p[ID_Fxy0] + z * p[ID_Fxy1] +
 390     x1 * (y1 * (z1 * p[ID_V000] + z * p[ID_V001])  +
 391           y  * (z1 * p[ID_V010] + z * p[ID_V011])) +
 392     x  * (y1 * (z1 * p[ID_V100] + z * p[ID_V101])  +
 393           y  * (z1 * p[ID_V110] + z * p[ID_V111]));
 394   thePoint -=
 395     x1 * (y1 * p[ID_E00z] + y * p[ID_E01z]) +
 396     x  * (y1 * p[ID_E10z] + y * p[ID_E11z]) +
 397     y1 * (z1 * p[ID_Ex00] + z * p[ID_Ex01]) +
 398     y  * (z1 * p[ID_Ex10] + z * p[ID_Ex11]) +
 399     z1 * (x1 * p[ID_E0y0] + x * p[ID_E1y0]) +
 400     z  * (x1 * p[ID_E0y1] + x * p[ID_E1y1]);
 401
 402   return true;
 403 }
 404
 405 //=======================================================================
 406 //function : Constructor
 407 //purpose  :
 408 //=======================================================================
 409
 410 SMESH_Block::SMESH_Block():
 411   myNbIterations(0),
 412   mySumDist(0.),
 413   myTolerance(-1.) // to be re-initialized
 414 {
 415 }
 416
 417
 418 //=======================================================================
 419 //function : NbVariables
 420 //purpose  :
 421 //=======================================================================
 422
 423 Standard_Integer SMESH_Block::NbVariables() const
 424 {
 425   return 3;
 426 }
 427
 428 //=======================================================================
 429 //function : NbEquations
 430 //purpose  :
 431 //=======================================================================
 432
 433 Standard_Integer SMESH_Block::NbEquations() const
 434 {
 435   return 1;
 436 }
 437
 438 //=======================================================================
 439 //function : Value
 440 //purpose  :
 441 //=======================================================================
 442
 443 Standard_Boolean SMESH_Block::Value(const math_Vector& theXYZ, math_Vector& theFxyz)
 444 {
 445   gp_XYZ P, params( theXYZ(1), theXYZ(2), theXYZ(3) );
 446   if ( params.IsEqual( myParam, DBL_MIN )) { // same param
 447     theFxyz( 1 ) = funcValue( myValues[ SQUARE_DIST ]);
 448   }
 449   else {
 450     ShellPoint( params, P );
 451     gp_Vec dP( P - myPoint );
 452     theFxyz(1) = funcValue( dP.SquareMagnitude() );
 453   }
 454   return true;
 455 }
 456
 457 //=======================================================================
 458 //function : Derivatives
 459 //purpose  :
 460 //=======================================================================
 461
 462 Standard_Boolean SMESH_Block::Derivatives(const math_Vector& XYZ,math_Matrix& Df)
 463 {
 464   math_Vector F(1,3);
 465   return Values(XYZ,F,Df);
 466 }
 467
 468 //=======================================================================
 469 //function : GetStateNumber
 470 //purpose  :
 471 //=======================================================================
 472
 473 Standard_Integer SMESH_Block::GetStateNumber ()
 474 {
 475   return 0; //myValues[0] < 1e-1;
 476 }
 477
 478 //=======================================================================
 479 //function : Values
 480 //purpose  :
 481 //=======================================================================
 482
 483 Standard_Boolean SMESH_Block::Values(const math_Vector& theXYZ,
 484                                      math_Vector&       theFxyz,
 485                                      math_Matrix&       theDf)
 486 {
 487   gp_XYZ P, params( theXYZ(1), theXYZ(2), theXYZ(3) );
 488   if ( params.IsEqual( myParam, DBL_MIN )) { // same param
 489     theFxyz( 1 )      = funcValue( myValues[ SQUARE_DIST ] );
 490     theDf( 1, DRV_1 ) = myValues[ DRV_1 ];
 491     theDf( 1, DRV_2 ) = myValues[ DRV_2 ];
 492     theDf( 1, DRV_3 ) = myValues[ DRV_3 ];
 493     return true;
 494   }
 495 #ifdef DEBUG_PARAM_COMPUTE
 496   cout << "PARAM GUESS: " << params.X() << " "<< params.Y() << " "<< params.X() << endl;
 497   myNbIterations++; // how many times call ShellPoint()
 498 #endif
 499   ShellPoint( params, P );
 500
 501   gp_Vec dP( myPoint, P );
 502   double sqDist = dP.SquareMagnitude();
 503   theFxyz(1) = funcValue( sqDist );
 504
 505   if ( sqDist < myTolerance * myTolerance ) { // a solution found
 506     myParam = params;
 507     myValues[ SQUARE_DIST ] = sqDist;
 508     theFxyz(1)  = theDf( 1,1 ) = theDf( 1,2 ) = theDf( 1,3 ) = 0;
 509     return true;
 510   }
 511
 512   if ( sqDist < myValues[ SQUARE_DIST ] ) // a better guess
 513   {
 514     // 3 partial derivatives
 515     gp_Vec drv[ 3 ]; // where we move with a small step in each direction
 516     for ( int iP = 1; iP <= 3; iP++ ) {
 517       if ( iP == myFaceIndex ) {
 518         drv[ iP - 1 ] = gp_Vec(0,0,0);
 519         continue;
 520       }
 521       gp_XYZ Pi;
 522       bool onEdge = ( theXYZ( iP ) + 0.001 > 1. );
 523       if ( onEdge )
 524         params.SetCoord( iP, theXYZ( iP ) - 0.001 );
 525       else
 526         params.SetCoord( iP, theXYZ( iP ) + 0.001 );
 527       ShellPoint( params, Pi );
 528       params.SetCoord( iP, theXYZ( iP ) ); // restore params
 529       gp_Vec dPi ( P, Pi );
 530       if ( onEdge ) dPi *= -1.;
 531       double mag = dPi.Magnitude();
 532       if ( mag > DBL_MIN )
 533         dPi /= mag;
 534       drv[ iP - 1 ] = dPi;
 535     }
 536     for ( int iP = 0; iP < 3; iP++ ) {
 537 #if 1
 538       theDf( 1, iP + 1 ) = dP * drv[iP];
 539 #else
 540       // Distance from P to plane passing through myPoint and defined
 541       // by the 2 other derivative directions:
 542       // like IntAna_IntConicQuad::Perform (const gp_Lin& L, const gp_Pln& P)
 543       // where L is (P -> myPoint), P is defined by the 2 other derivative direction
 544       int iPrev = ( iP ? iP - 1 : 2 );
 545       int iNext = ( iP == 2 ? 0 : iP + 1 );
 546       gp_Vec plnNorm = drv[ iPrev ].Crossed( drv [ iNext ] );
 547       double Direc = plnNorm * drv[ iP ];
 548       if ( Abs(Direc) <= DBL_MIN )
 549         theDf( 1, iP + 1 ) = dP * drv[ iP ];
 550       else {
 551         double Dis = plnNorm * P - plnNorm * myPoint;
 552         theDf( 1, iP + 1 ) = Dis/Direc;
 553       }
 554 #endif
 555     }
 556 #ifdef DEBUG_PARAM_COMPUTE
 557     cout << "F = " << theFxyz(1) <<
 558       " DRV: " << theDf(1,1) << " " << theDf(1,2) << " " << theDf(1,3)  << endl;
 559     myNbIterations +=3; // how many times call ShellPoint()
 560 #endif
 561
 562     // store better values
 563     myParam              = params;
 564     myValues[SQUARE_DIST]= sqDist;
 565     myValues[DRV_1]      = theDf(1,DRV_1);
 566     myValues[DRV_2]      = theDf(1,DRV_2);
 567     myValues[DRV_3]      = theDf(1,DRV_3);
 568   }
 569
 570   return true;
 571 }
 572
 573 //============================================================================
 574 //function : computeParameters
 575 //purpose  : compute point parameters in the block using math_FunctionSetRoot
 576 //============================================================================
 577
 578 bool SMESH_Block::computeParameters(const gp_Pnt& thePoint,
 579                                     gp_XYZ&       theParams,
 580                                     const gp_XYZ& theParamsHint)
 581 {
 582   myPoint = thePoint.XYZ();
 583
 584   myParam.SetCoord( -1,-1,-1 );
 585   myValues[ SQUARE_DIST ] = 1e100;
 586
 587   math_Vector low  ( 1, 3, 0.0 );
 588   math_Vector up   ( 1, 3, 1.0 );
 589   math_Vector tol  ( 1, 3, 1e-4 );
 590   math_Vector start( 1, 3, 0.0 );
 591   start( 1 ) = theParamsHint.X();
 592   start( 2 ) = theParamsHint.Y();
 593   start( 3 ) = theParamsHint.Z();
 594
 595   math_FunctionSetRoot paramSearch( *this, tol );
 596
 597   mySquareFunc = 0; // large approaching steps
 598   //if ( hasHint ) mySquareFunc = 1; // small approaching steps
 599
 600   double loopTol = 10 * myTolerance;
 601   int nbLoops = 0;
 602   while ( distance() > loopTol && nbLoops <= 3 )
 603   {
 604     paramSearch.Perform ( *static_cast<math_FunctionSetWithDerivatives*>(this),
 605                           start, low, up );
 606     start( 1 ) = myParam.X();
 607     start( 2 ) = myParam.Y();
 608     start( 3 ) = myParam.Z();
 609     mySquareFunc = !mySquareFunc;
 610     nbLoops++;
 611   }
 612 #ifdef DEBUG_PARAM_COMPUTE
 613   mySumDist += distance();
 614   cout << " ------ SOLUTION: ( "<< myParam.X() <<" "<< myParam.Y() <<" "<< myParam.Z() <<" )"<<endl
 615        << " ------ DIST : " << distance() << "\t Tol=" << myTolerance << "\t Nb LOOPS=" << nbLoops << endl
 616        << " ------ NB IT: " << myNbIterations << ",  SUM DIST: " << mySumDist << endl;
 617 #endif
 618
 619   theParams = myParam;
 620
 621   if ( myFaceIndex > 0 )
 622     theParams.SetCoord( myFaceIndex, myFaceParam );
 623
 624   return true;
 625 }
 626
 627 //=======================================================================
 628 //function : ComputeParameters
 629 //purpose  : compute point parameters in the block
 630 //=======================================================================
 631
 632 bool SMESH_Block::ComputeParameters(const gp_Pnt& thePoint,
 633                                     gp_XYZ&       theParams,
 634                                     const int     theShapeID,
 635                                     const gp_XYZ& theParamsHint)
 636 {
 637   if ( VertexParameters( theShapeID, theParams ))
 638     return true;
 639
 640   if ( IsEdgeID( theShapeID )) {
 641     TEdge& e = myEdge[ theShapeID - ID_FirstE ];
 642     Adaptor3d_Curve* curve = e.GetCurve();
 643     Extrema_ExtPC anExtPC( thePoint, *curve, curve->FirstParameter(), curve->LastParameter() );
 644     int i, nb = anExtPC.IsDone() ? anExtPC.NbExt() : 0;
 645     for ( i = 1; i <= nb; i++ ) {
 646       if ( anExtPC.IsMin( i ))
 647         return EdgeParameters( theShapeID, anExtPC.Point( i ).Parameter(), theParams );
 648     }
 649     return false;
 650   }
 651
 652   const bool isOnFace = IsFaceID( theShapeID );
 653   double * coef = GetShapeCoef( theShapeID );
 654
 655   // Find the first guess paremeters
 656
 657   gp_XYZ start(0, 0, 0);
 658
 659   bool hasHint = ( 0 <= theParamsHint.X() && theParamsHint.X() <= 1 &&
 660                    0 <= theParamsHint.Y() && theParamsHint.Y() <= 1 &&
 661                    0 <= theParamsHint.Y() && theParamsHint.Y() <= 1 );
 662   if ( !hasHint && !myGridComputed )
 663   {
 664     // define the first guess by thePoint projection on lines
 665     // connecting vertices
 666     bool needGrid = false;
 667     gp_XYZ par000( 0, 0, 0 ), par111( 1, 1, 1 );
 668     double zero = DBL_MIN * DBL_MIN;
 669     for ( int iEdge = 0, iParam = 1; iParam <= 3 && !needGrid; iParam++ )
 670     {
 671       if ( isOnFace && coef[ iParam - 1 ] != 0 ) {
 672         iEdge += 4;
 673         continue;
 674       }
 675       double sumParam = 0;
 676       for ( int iE = 0; iE < 4; iE++, iEdge++ ) { // loop on 4 parallel edges
 677         gp_Pnt p0 = myEdge[ iEdge ].Point( par000 );
 678         gp_Pnt p1 = myEdge[ iEdge ].Point( par111 );
 679         gp_Vec v01( p0, p1 ), v0P( p0, thePoint );
 680         double len2 = v01.SquareMagnitude();
 681         double par = 0;
 682         if ( len2 > zero ) {
 683           par = v0P.Dot( v01 ) / len2;
 684           if ( par < 0 || par > 1 ) { // projection falls out of line ends => needGrid
 685             needGrid = true;
 686             break;
 687           }
 688         }
 689         sumParam += par;
 690       }
 691       start.SetCoord( iParam, sumParam / 4.);
 692     }
 693     if ( needGrid ) {
 694       // compute nodes of 3 x 3 x 3 grid
 695       int iNode = 0;
 696       Bnd_Box box;
 697       for ( double x = 0.25; x < 0.9; x += 0.25 )
 698         for ( double y = 0.25; y < 0.9; y += 0.25 )
 699           for ( double z = 0.25; z < 0.9; z += 0.25 ) {
 700             TxyzPair & prmPtn = my3x3x3GridNodes[ iNode++ ];
 701             prmPtn.first.SetCoord( x, y, z );
 702             ShellPoint( prmPtn.first, prmPtn.second );
 703             box.Add( gp_Pnt( prmPtn.second ));
 704           }
 705       myGridComputed = true;
 706       myTolerance = sqrt( box.SquareExtent() ) * 1e-5;
 707     }
 708   }
 709
 710   if ( hasHint )
 711   {
 712     start = theParamsHint;
 713   }
 714   else if ( myGridComputed )
 715   {
 716     double minDist = DBL_MAX;
 717     gp_XYZ* bestParam = 0;
 718     for ( int iNode = 0; iNode < 27; iNode++ ) {
 719       TxyzPair & prmPtn = my3x3x3GridNodes[ iNode ];
 720       double dist = ( thePoint.XYZ() - prmPtn.second ).SquareModulus();
 721       if ( dist < minDist ) {
 722         minDist = dist;
 723         bestParam = & prmPtn.first;
 724       }
 725     }
 726     start = *bestParam;
 727   }
 728
 729   myFaceIndex = -1;
 730   myFaceParam = 0.;
 731   if ( isOnFace ) {
 732     // put a point on the face
 733     for ( int iCoord = 0; iCoord < 3; iCoord++ )
 734       if ( coef[ iCoord ] ) {
 735         myFaceIndex = iCoord + 1;
 736         myFaceParam = ( coef[ iCoord ] < 0.5 ) ? 0.0 : 1.0;
 737         start.SetCoord( myFaceIndex, myFaceParam );
 738       }
 739   }
 740
 741 #ifdef DEBUG_PARAM_COMPUTE
 742   cout << " #### POINT " <<thePoint.X()<<" "<<thePoint.Y()<<" "<<thePoint.Z()<<" ####"<< endl;
 743 #endif
 744
 745   if ( myTolerance < 0 ) myTolerance = 1e-6;
 746
 747   const double parDelta = 1e-4;
 748   const double sqTolerance = myTolerance * myTolerance;
 749
 750   gp_XYZ solution = start, params = start;
 751   double sqDistance = 1e100;
 752   int nbLoops = 0, nbGetWorst = 0;
 753
 754   while ( nbLoops <= 100 )
 755   {
 756     gp_XYZ P, Pi;
 757     ShellPoint( params, P );
 758
 759     gp_Vec dP( thePoint, P );
 760     double sqDist = dP.SquareMagnitude();
 761
 762     if ( sqDist > sqDistance ) { // solution get worse
 763       if ( ++nbGetWorst > 2 )
 764         return computeParameters( thePoint, theParams, solution );
 765     }
 766 #ifdef DEBUG_PARAM_COMPUTE
 767     cout << "PARAMS: ( " << params.X() <<" "<< params.Y() <<" "<< params.Z() <<" )"<< endl;
 768     cout << "DIST: " << sqrt( sqDist ) << endl;
 769 #endif
 770
 771     if ( sqDist < sqDistance ) { // get better
 772       sqDistance = sqDist;
 773       solution   = params;
 774       nbGetWorst = 0;
 775       if ( sqDistance < sqTolerance ) // a solution found
 776         break;
 777     }
 778
 779         // look for a next better solution
 780     for ( int iP = 1; iP <= 3; iP++ ) {
 781       if ( iP == myFaceIndex )
 782         continue;
 783       // see where we move with a small (=parDelta) step in this direction
 784       gp_XYZ nearParams = params;
 785       bool onEdge = ( params.Coord( iP ) + parDelta > 1. );
 786       if ( onEdge )
 787         nearParams.SetCoord( iP, params.Coord( iP ) - parDelta );
 788       else
 789         nearParams.SetCoord( iP, params.Coord( iP ) + parDelta );
 790       ShellPoint( nearParams, Pi );
 791       gp_Vec dPi ( P, Pi );
 792       if ( onEdge ) dPi *= -1.;
 793       // modify a parameter
 794       double mag = dPi.Magnitude();
 795       if ( mag < DBL_MIN )
 796         continue;
 797       gp_Vec dir = dPi / mag; // dir we move modifying the parameter
 798       double dist = dir * dP; // where we should get to
 799       double dPar = dist / mag * parDelta; // predict parameter change
 800       double curPar = params.Coord( iP );
 801       double par = curPar - dPar; // new parameter value
 802       while ( par > 1 || par < 0 ) {
 803         dPar /= 2.;
 804         par = curPar - dPar;
 805       }
 806       params.SetCoord( iP, par );
 807     }
 808
 809     nbLoops++;
 810   }
 811 #ifdef DEBUG_PARAM_COMPUTE
 812   myNbIterations += nbLoops*4; // how many times ShellPoint called
 813   mySumDist += sqrt( sqDistance );
 814   cout << " ------ SOLUTION: ( "<<solution.X()<<" "<<solution.Y()<<" "<<solution.Z()<<" )"<<endl
 815        << " ------ DIST : " << sqrt( sqDistance ) << "\t Tol=" << myTolerance << "\t Nb LOOPS=" << nbLoops << endl
 816        << " ------ NB IT: " << myNbIterations << ",  SUM DIST: " << mySumDist << endl;
 817 #endif
 818
 819   theParams = solution;
 820
 821   if ( myFaceIndex > 0 )
 822     theParams.SetCoord( myFaceIndex, myFaceParam );
 823
 824   return true;
 825 }
 826
 827 //=======================================================================
 828 //function : VertexParameters
 829 //purpose  : return parameters of a vertex given by TShapeID
 830 //=======================================================================
 831
 832 bool SMESH_Block::VertexParameters(const int theVertexID, gp_XYZ& theParams)
 833 {
 834   switch ( theVertexID ) {
 835   case ID_V000: theParams.SetCoord(0., 0., 0.); return true;
 836   case ID_V100: theParams.SetCoord(1., 0., 0.); return true;
 837   case ID_V110: theParams.SetCoord(1., 1., 0.); return true;
 838   case ID_V010: theParams.SetCoord(0., 1., 0.); return true;
 839   default:;
 840   }
 841   return false;
 842 }
 843
 844 //=======================================================================
 845 //function : EdgeParameters
 846 //purpose  : return parameters of a point given by theU on edge
 847 //=======================================================================
 848
 849 bool SMESH_Block::EdgeParameters(const int theEdgeID, const double theU, gp_XYZ& theParams)
 850 {
 851   if ( IsEdgeID( theEdgeID )) {
 852     vector< int > vertexVec;
 853     GetEdgeVertexIDs( theEdgeID, vertexVec );
 854     VertexParameters( vertexVec[0], theParams );
 855     TEdge& e = myEdge[ theEdgeID - ID_Ex00 ];
 856     double param = ( theU - e.EndParam(0) ) / ( e.EndParam(1) - e.EndParam(0) );
 857     theParams.SetCoord( e.CoordInd(), param );
 858     return true;
 859   }
 860   return false;
 861 }
 862
 863 //=======================================================================
 864 //function : DumpShapeID
 865 //purpose  : debug an id of a block sub-shape
 866 //=======================================================================
 867
 868 #define CASEDUMP(id,strm) case id: strm << #id; break;
 869
 870 ostream& SMESH_Block::DumpShapeID (const int id, ostream& stream)
 871 {
 872   switch ( id ) {
 873   CASEDUMP( ID_V000, stream );
 874   CASEDUMP( ID_V100, stream );
 875   CASEDUMP( ID_V010, stream );
 876   CASEDUMP( ID_V110, stream );
 877   CASEDUMP( ID_V001, stream );
 878   CASEDUMP( ID_V101, stream );
 879   CASEDUMP( ID_V011, stream );
 880   CASEDUMP( ID_V111, stream );
 881   CASEDUMP( ID_Ex00, stream );
 882   CASEDUMP( ID_Ex10, stream );
 883   CASEDUMP( ID_Ex01, stream );
 884   CASEDUMP( ID_Ex11, stream );
 885   CASEDUMP( ID_E0y0, stream );
 886   CASEDUMP( ID_E1y0, stream );
 887   CASEDUMP( ID_E0y1, stream );
 888   CASEDUMP( ID_E1y1, stream );
 889   CASEDUMP( ID_E00z, stream );
 890   CASEDUMP( ID_E10z, stream );
 891   CASEDUMP( ID_E01z, stream );
 892   CASEDUMP( ID_E11z, stream );
 893   CASEDUMP( ID_Fxy0, stream );
 894   CASEDUMP( ID_Fxy1, stream );
 895   CASEDUMP( ID_Fx0z, stream );
 896   CASEDUMP( ID_Fx1z, stream );
 897   CASEDUMP( ID_F0yz, stream );
 898   CASEDUMP( ID_F1yz, stream );
 899   CASEDUMP( ID_Shell, stream );
 900   default: stream << "ID_INVALID";
 901   }
 902   return stream;
 903 }
 904
 905 //=======================================================================
 906 //function : GetShapeIDByParams
 907 //purpose  : define an id of the block sub-shape by normlized point coord
 908 //=======================================================================
 909
 910 int SMESH_Block::GetShapeIDByParams ( const gp_XYZ& theCoord )
 911 {
 912   //   id ( 0 - 26 ) computation:
 913
 914   //   vertex     ( 0 - 7 )  : id = 1*x + 2*y + 4*z
 915
 916   //   edge || X  ( 8 - 11 ) : id = 8   + 1*y + 2*z
 917   //   edge || Y  ( 12 - 15 ): id = 1*x + 12  + 2*z
 918   //   edge || Z  ( 16 - 19 ): id = 1*x + 2*y + 16
 919
 920   //   face || XY ( 20 - 21 ): id = 8   + 12  + 1*z - 0
 921   //   face || XZ ( 22 - 23 ): id = 8   + 1*y + 16  - 2
 922   //   face || YZ ( 24 - 25 ): id = 1*x + 12  + 16  - 4
 923
 924   static int iAddBnd[]    = { 1, 2, 4 };
 925   static int iAddNotBnd[] = { 8, 12, 16 };
 926   static int iFaceSubst[] = { 0, 2, 4 };
 927
 928   int id = 0;
 929   int iOnBoundary = 0;
 930   for ( int iCoord = 0; iCoord < 3; iCoord++ )
 931   {
 932     double val = theCoord.Coord( iCoord + 1 );
 933     if ( val == 0.0 )
 934       iOnBoundary++;
 935     else if ( val == 1.0 )
 936       id += iAddBnd[ iOnBoundary++ ];
 937     else
 938       id += iAddNotBnd[ iCoord ];
 939   }
 940   if ( iOnBoundary == 1 ) // face
 941     id -= iFaceSubst[ (id - 20) / 4 ];
 942   else if ( iOnBoundary == 0 ) // shell
 943     id = 26;
 944
 945   if ( id > 26 || id < 0 ) {
 946     MESSAGE( "GetShapeIDByParams() = " << id
 947             <<" "<< theCoord.X() <<" "<< theCoord.Y() <<" "<< theCoord.Z() );
 948   }
 949
 950   return id + 1; // shape ids start at 1
 951 }
 952
 953 //=======================================================================
 954 //function : GetOrderedEdges
 955 //purpose  : return nb wires and a list of oredered edges
 956 //=======================================================================
 957
 958 int SMESH_Block::GetOrderedEdges (const TopoDS_Face&   theFace,
 959                                   TopoDS_Vertex        theFirstVertex,
 960                                   list< TopoDS_Edge >& theEdges,
 961                                   list< int >  &       theNbVertexInWires)
 962 {
 963   // put wires in a list, so that an outer wire comes first
 964   list<TopoDS_Wire> aWireList;
 965   TopoDS_Wire anOuterWire = BRepTools::OuterWire( theFace );
 966   aWireList.push_back( anOuterWire );
 967   for ( TopoDS_Iterator wIt (theFace); wIt.More(); wIt.Next() )
 968     if ( !anOuterWire.IsSame( wIt.Value() ))
 969       aWireList.push_back( TopoDS::Wire( wIt.Value() ));
 970
 971   // loop on edges of wires
 972   theNbVertexInWires.clear();
 973   list<TopoDS_Wire>::iterator wlIt = aWireList.begin();
 974   for ( ; wlIt != aWireList.end(); wlIt++ )
 975   {
 976     int iE;
 977     BRepTools_WireExplorer wExp( *wlIt, theFace );
 978     for ( iE = 0; wExp.More(); wExp.Next(), iE++ )
 979     {
 980       TopoDS_Edge edge = wExp.Current();
 981       edge = TopoDS::Edge( edge.Oriented( wExp.Orientation() ));
 982       theEdges.push_back( edge );
 983     }
 984     theNbVertexInWires.push_back( iE );
 985     iE = 0;
 986     if ( wlIt == aWireList.begin() && theEdges.size() > 1 ) { // the outer wire
 987       // orient closed edges
 988       list< TopoDS_Edge >::iterator eIt, eIt2;
 989       for ( eIt = theEdges.begin(); eIt != theEdges.end(); eIt++ )
 990       {
 991         TopoDS_Edge& edge = *eIt;
 992         if ( TopExp::FirstVertex( edge ).IsSame( TopExp::LastVertex( edge ) ))
 993         {
 994           eIt2 = eIt;
 995           bool isNext = ( eIt2 == theEdges.begin() );
 996           TopoDS_Edge edge2 = isNext ? *(++eIt2) : *(--eIt2);
 997           double f1,l1,f2,l2;
 998           Handle(Geom2d_Curve) c1 = BRep_Tool::CurveOnSurface( edge, theFace, f1,l1 );
 999           Handle(Geom2d_Curve) c2 = BRep_Tool::CurveOnSurface( edge2, theFace, f2,l2 );
1000           gp_Pnt2d pf = c1->Value( edge.Orientation() == TopAbs_FORWARD ? f1 : l1 );
1001           gp_Pnt2d pl = c1->Value( edge.Orientation() == TopAbs_FORWARD ? l1 : f1 );
1002           bool isFirst = ( edge2.Orientation() == TopAbs_FORWARD ? isNext : !isNext );
1003           gp_Pnt2d p2 = c2->Value( isFirst ? f2 : l2 );
1004           isFirst = ( p2.SquareDistance( pf ) < p2.SquareDistance( pl ));
1005           if ( isNext ? isFirst : !isFirst )
1006             edge.Reverse();
1007           // to make a seam go first
1008           if ( theFirstVertex.IsNull() )
1009             theFirstVertex = TopExp::FirstVertex( edge, true );
1010         }
1011       }
1012       // rotate theEdges until it begins from theFirstVertex
1013       if ( ! theFirstVertex.IsNull() ) {
1014         TopoDS_Vertex vv[2];
1015         TopExp::Vertices( theEdges.front(), vv[0], vv[1], true );
1016         // on closed face, make seam edge the first in the list
1017         while ( !vv[0].IsSame( theFirstVertex ) || vv[0].IsSame( vv[1] ))
1018         {
1019           theEdges.splice(theEdges.end(), theEdges,
1020                           theEdges.begin(), ++theEdges.begin());
1021           TopExp::Vertices( theEdges.front(), vv[0], vv[1], true );
1022           if ( iE++ > theNbVertexInWires.back() ) {
1023 #ifdef _DEBUG_
1024             gp_Pnt p = BRep_Tool::Pnt( theFirstVertex );
1025             cout << " : Warning : vertex "<< theFirstVertex.TShape().operator->()
1026                  << " ( " << p.X() << " " << p.Y() << " " << p.Z() << " )"
1027                  << " not found in outer wire of face "<< theFace.TShape().operator->()
1028                  << " with vertices: " <<  endl;
1029             wExp.Init( *wlIt, theFace );
1030             for ( int i = 0; wExp.More(); wExp.Next(), i++ )
1031             {
1032               TopoDS_Edge edge = wExp.Current();
1033               edge = TopoDS::Edge( edge.Oriented( wExp.Orientation() ));
1034               TopoDS_Vertex v = TopExp::FirstVertex( edge, true );
1035               gp_Pnt p = BRep_Tool::Pnt( v );
1036               cout << i << " " << v.TShape().operator->() << " "
1037                    << p.X() << " " << p.Y() << " " << p.Z() << " " << endl;
1038             }
1039 #endif
1040             break; // break infinite loop
1041           }
1042         }
1043       }
1044     } // end outer wire
1045   }
1046
1047   return aWireList.size();
1048 }
1049 //================================================================================
1050 /*!
1051  * \brief Call it after geometry initialisation
1052  */
1053 //================================================================================
1054
1055 void SMESH_Block::init()
1056 {
1057   myNbIterations = 0;
1058   mySumDist = 0;
1059   myGridComputed = false;
1060 }
1061
1062 //=======================================================================
1063 //function : LoadMeshBlock
1064 //purpose  : prepare to work with theVolume
1065 //=======================================================================
1066
1067 #define gpXYZ(n) gp_XYZ(n->X(),n->Y(),n->Z())
1068
1069 bool SMESH_Block::LoadMeshBlock(const SMDS_MeshVolume*        theVolume,
1070                                 const int                     theNode000Index,
1071                                 const int                     theNode001Index,
1072                                 vector<const SMDS_MeshNode*>& theOrderedNodes)
1073 {
1074   MESSAGE(" ::LoadMeshBlock()");
1075   init();
1076
1077   SMDS_VolumeTool vTool;
1078   if (!vTool.Set( theVolume ) || vTool.NbNodes() != 8 ||
1079       !vTool.IsLinked( theNode000Index, theNode001Index )) {
1080     MESSAGE(" Bad arguments ");
1081     return false;
1082   }
1083   vTool.SetExternalNormal();
1084   // In terms of indices used for access to nodes and faces in SMDS_VolumeTool:
1085   int V000, V100, V010, V110, V001, V101, V011, V111; // 8 vertices
1086   int Fxy0, Fxy1; // bottom and top faces
1087   // vertices of faces
1088   vector<int> vFxy0, vFxy1;
1089
1090   V000 = theNode000Index;
1091   V001 = theNode001Index;
1092
1093   // get faces sharing V000 and V001
1094   list<int> fV000, fV001;
1095   int i, iF, iE, iN;
1096   for ( iF = 0; iF < vTool.NbFaces(); ++iF ) {
1097     const int* nid = vTool.GetFaceNodesIndices( iF );
1098     for ( iN = 0; iN < 4; ++iN )
1099       if ( nid[ iN ] == V000 ) {
1100         fV000.push_back( iF );
1101       } else if ( nid[ iN ] == V001 ) {
1102         fV001.push_back( iF );
1103       }
1104   }
1105
1106   // find the bottom (Fxy0), the top (Fxy1) faces
1107   list<int>::iterator fIt1, fIt2, Fxy0Pos;
1108   for ( fIt1 = fV000.begin(); fIt1 != fV000.end(); fIt1++) {
1109     fIt2 = std::find( fV001.begin(), fV001.end(), *fIt1 );
1110     if ( fIt2 != fV001.end() ) { // *fIt1 is in the both lists
1111       fV001.erase( fIt2 ); // erase Fx0z or F0yz from fV001
1112     } else { // *fIt1 is in fV000 only
1113       Fxy0Pos = fIt1; // points to Fxy0
1114     }
1115   }
1116   Fxy0 = *Fxy0Pos;
1117   Fxy1 = fV001.front();
1118   const SMDS_MeshNode** nn = vTool.GetNodes();
1119
1120   // find bottom veritices, their order is that a face normal is external
1121   vFxy0.resize(4);
1122   const int* nid = vTool.GetFaceNodesIndices( Fxy0 );
1123   for ( i = 0; i < 4; ++i )
1124     if ( nid[ i ] == V000 )
1125       break;
1126   for ( iN = 0; iN < 4; ++iN, ++i ) {
1127     if ( i == 4 ) i = 0;
1128     vFxy0[ iN ] = nid[ i ];
1129   }
1130   // find top veritices, their order is that a face normal is external
1131   vFxy1.resize(4);
1132   nid = vTool.GetFaceNodesIndices( Fxy1 );
1133   for ( i = 0; i < 4; ++i )
1134     if ( nid[ i ] == V001 )
1135       break;
1136   for ( iN = 0; iN < 4; ++iN, ++i ) {
1137     if ( i == 4 ) i = 0;
1138     vFxy1[ iN ] = nid[ i ];
1139   }
1140   // find indices of the rest veritices
1141   V100 = vFxy0[3];
1142   V010 = vFxy0[1];
1143   V110 = vFxy0[2];
1144   V101 = vFxy1[1];
1145   V011 = vFxy1[3];
1146   V111 = vFxy1[2];
1147
1148   // set points coordinates
1149   myPnt[ ID_V000 - 1 ] = gpXYZ( nn[ V000 ] );
1150   myPnt[ ID_V100 - 1 ] = gpXYZ( nn[ V100 ] );
1151   myPnt[ ID_V010 - 1 ] = gpXYZ( nn[ V010 ] );
1152   myPnt[ ID_V110 - 1 ] = gpXYZ( nn[ V110 ] );
1153   myPnt[ ID_V001 - 1 ] = gpXYZ( nn[ V001 ] );
1154   myPnt[ ID_V101 - 1 ] = gpXYZ( nn[ V101 ] );
1155   myPnt[ ID_V011 - 1 ] = gpXYZ( nn[ V011 ] );
1156   myPnt[ ID_V111 - 1 ] = gpXYZ( nn[ V111 ] );
1157
1158   // fill theOrderedNodes
1159   theOrderedNodes.resize( 8 );
1160   theOrderedNodes[ 0 ] = nn[ V000 ];
1161   theOrderedNodes[ 1 ] = nn[ V100 ];
1162   theOrderedNodes[ 2 ] = nn[ V010 ];
1163   theOrderedNodes[ 3 ] = nn[ V110 ];
1164   theOrderedNodes[ 4 ] = nn[ V001 ];
1165   theOrderedNodes[ 5 ] = nn[ V101 ];
1166   theOrderedNodes[ 6 ] = nn[ V011 ];
1167   theOrderedNodes[ 7 ] = nn[ V111 ];
1168
1169   // fill edges
1170   vector< int > vertexVec;
1171   for ( iE = 0; iE < NbEdges(); ++iE ) {
1172     GetEdgeVertexIDs(( iE + ID_FirstE ), vertexVec );
1173     myEdge[ iE ].Set(( iE + ID_FirstE ),
1174                      myPnt[ vertexVec[0] - 1 ],
1175                      myPnt[ vertexVec[1] - 1 ]);
1176   }
1177
1178   // fill faces' corners
1179   for ( iF = ID_Fxy0; iF < ID_Shell; ++iF )
1180   {
1181     TFace& tFace = myFace[ iF - ID_FirstF ];
1182     vector< int > edgeIdVec(4, -1);
1183     GetFaceEdgesIDs( iF, edgeIdVec );
1184     tFace.Set( iF, myEdge[ edgeIdVec [ 0 ] - ID_Ex00], myEdge[ edgeIdVec [ 1 ] - ID_Ex00]);
1185   }
1186
1187   return true;
1188 }
1189
1190 //=======================================================================
1191 //function : LoadBlockShapes
1192 //purpose  : Initialize block geometry with theShell,
1193 //           add sub-shapes of theBlock to theShapeIDMap so that they get
1194 //           IDs acoording to enum TShapeID
1195 //=======================================================================
1196
1197 bool SMESH_Block::LoadBlockShapes(const TopoDS_Shell&         theShell,
1198                                   const TopoDS_Vertex&        theVertex000,
1199                                   const TopoDS_Vertex&        theVertex001,
1200                                   TopTools_IndexedMapOfOrientedShape& theShapeIDMap )
1201 {
1202   MESSAGE(" ::LoadBlockShapes()");
1203   return ( FindBlockShapes( theShell, theVertex000, theVertex001, theShapeIDMap ) &&
1204            LoadBlockShapes( theShapeIDMap ));
1205 }
1206
1207 //=======================================================================
1208 //function : LoadBlockShapes
1209 //purpose  : add sub-shapes of theBlock to theShapeIDMap so that they get
1210 //           IDs acoording to enum TShapeID
1211 //=======================================================================
1212
1213 bool SMESH_Block::FindBlockShapes(const TopoDS_Shell&         theShell,
1214                                   const TopoDS_Vertex&        theVertex000,
1215                                   const TopoDS_Vertex&        theVertex001,
1216                                   TopTools_IndexedMapOfOrientedShape& theShapeIDMap )
1217 {
1218   MESSAGE(" ::FindBlockShapes()");
1219
1220   // 8 vertices
1221   TopoDS_Shape V000, V100, V010, V110, V001, V101, V011, V111;
1222   // 12 edges
1223   TopoDS_Shape Ex00, Ex10, Ex01, Ex11;
1224   TopoDS_Shape E0y0, E1y0, E0y1, E1y1;
1225   TopoDS_Shape E00z, E10z, E01z, E11z;
1226   // 6 faces
1227   TopoDS_Shape Fxy0, Fx0z, F0yz, Fxy1, Fx1z, F1yz;
1228
1229   // nb of faces bound to a vertex in TopTools_IndexedDataMapOfShapeListOfShape
1230   // filled by TopExp::MapShapesAndAncestors()
1231   const int NB_FACES_BY_VERTEX = 6;
1232
1233   TopTools_IndexedDataMapOfShapeListOfShape vfMap;
1234   TopExp::MapShapesAndAncestors( theShell, TopAbs_VERTEX, TopAbs_FACE, vfMap );
1235   if ( vfMap.Extent() != 8 ) {
1236     MESSAGE(" Wrong nb of vertices in the block: " << vfMap.Extent() );
1237     return false;
1238   }
1239
1240   V000 = theVertex000;
1241   V001 = theVertex001;
1242
1243   if ( V000.IsNull() ) {
1244     // find vertex 000 - the one with smallest coordinates
1245     double minVal = DBL_MAX, minX, val;
1246     for ( int i = 1; i <= 8; i++ ) {
1247       const TopoDS_Vertex& v = TopoDS::Vertex( vfMap.FindKey( i ));
1248       gp_Pnt P = BRep_Tool::Pnt( v );
1249       val = P.X() + P.Y() + P.Z();
1250       if ( val < minVal || ( val == minVal && P.X() < minX )) {
1251         V000 = v;
1252         minVal = val;
1253         minX = P.X();
1254       }
1255     }
1256     // find vertex 001 - the one on the most vertical edge passing through V000
1257     TopTools_IndexedDataMapOfShapeListOfShape veMap;
1258     TopExp::MapShapesAndAncestors( theShell, TopAbs_VERTEX, TopAbs_EDGE, veMap );
1259     gp_Vec dir001 = gp::DZ();
1260     gp_Pnt p000 = BRep_Tool::Pnt( TopoDS::Vertex( V000 ));
1261     double maxVal = -DBL_MAX;
1262     TopTools_ListIteratorOfListOfShape eIt ( veMap.FindFromKey( V000 ));
1263     for (  ; eIt.More(); eIt.Next() ) {
1264       const TopoDS_Edge& e = TopoDS::Edge( eIt.Value() );
1265       TopoDS_Vertex v = TopExp::FirstVertex( e );
1266       if ( v.IsSame( V000 ))
1267         v = TopExp::LastVertex( e );
1268       val = dir001 * gp_Vec( p000, BRep_Tool::Pnt( v )).Normalized();
1269       if ( val > maxVal ) {
1270         V001 = v;
1271         maxVal = val;
1272       }
1273     }
1274   }
1275
1276   // find the bottom (Fxy0), Fx0z and F0yz faces
1277
1278   const TopTools_ListOfShape& f000List = vfMap.FindFromKey( V000 );
1279   const TopTools_ListOfShape& f001List = vfMap.FindFromKey( V001 );
1280   if (f000List.Extent() != NB_FACES_BY_VERTEX ||
1281       f001List.Extent() != NB_FACES_BY_VERTEX ) {
1282     MESSAGE(" LoadBlockShapes() " << f000List.Extent() << " " << f001List.Extent());
1283     return false;
1284   }
1285   TopTools_ListIteratorOfListOfShape f001It, f000It ( f000List );
1286   int i, j, iFound1, iFound2;
1287   for ( j = 0; f000It.More(); f000It.Next(), j++ )
1288   {
1289     if ( NB_FACES_BY_VERTEX == 6 && j % 2 ) continue; // each face encounters twice
1290     const TopoDS_Shape& F = f000It.Value();
1291     for ( i = 0, f001It.Initialize( f001List ); f001It.More(); f001It.Next(), i++ ) {
1292       if ( NB_FACES_BY_VERTEX == 6 && i % 2 ) continue; // each face encounters twice
1293       if ( F.IsSame( f001It.Value() ))
1294         break;
1295     }
1296     if ( f001It.More() ) // Fx0z or F0yz found
1297       if ( Fx0z.IsNull() ) {
1298         Fx0z = F;
1299         iFound1 = i;
1300       } else {
1301         F0yz = F;
1302         iFound2 = i;
1303       }
1304     else // F is the bottom face
1305       Fxy0 = F;
1306   }
1307   if ( Fxy0.IsNull() || Fx0z.IsNull() || F0yz.IsNull() ) {
1308     MESSAGE( Fxy0.IsNull() <<" "<< Fx0z.IsNull() <<" "<< F0yz.IsNull() );
1309     return false;
1310   }
1311
1312   // choose the top face (Fxy1)
1313   for ( i = 0, f001It.Initialize( f001List ); f001It.More(); f001It.Next(), i++ ) {
1314     if ( NB_FACES_BY_VERTEX == 6 && i % 2 ) continue; // each face encounters twice
1315     if ( i != iFound1 && i != iFound2 )
1316       break;
1317   }
1318   Fxy1 = f001It.Value();
1319   if ( Fxy1.IsNull() ) {
1320     MESSAGE(" LoadBlockShapes() error ");
1321     return false;
1322   }
1323
1324   // find bottom edges and veritices
1325   list< TopoDS_Edge > eList;
1326   list< int >         nbVertexInWires;
1327   GetOrderedEdges( TopoDS::Face( Fxy0 ), TopoDS::Vertex( V000 ), eList, nbVertexInWires );
1328   if ( nbVertexInWires.size() != 1 || nbVertexInWires.front() != 4 ) {
1329     MESSAGE(" LoadBlockShapes() error ");
1330     return false;
1331   }
1332   list< TopoDS_Edge >::iterator elIt = eList.begin();
1333   for ( i = 0; elIt != eList.end(); elIt++, i++ )
1334     switch ( i ) {
1335     case 0: E0y0 = *elIt; V010 = TopExp::LastVertex( *elIt, true ); break;
1336     case 1: Ex10 = *elIt; V110 = TopExp::LastVertex( *elIt, true ); break;
1337     case 2: E1y0 = *elIt; V100 = TopExp::LastVertex( *elIt, true ); break;
1338     case 3: Ex00 = *elIt; break;
1339     default:;
1340     }
1341   if ( i != 4 || E0y0.IsNull() || Ex10.IsNull() || E1y0.IsNull() || Ex00.IsNull() ) {
1342     MESSAGE(" LoadBlockShapes() error, eList.size()=" << eList.size());
1343     return false;
1344   }
1345
1346
1347   // find top edges and veritices
1348   eList.clear();
1349   GetOrderedEdges( TopoDS::Face( Fxy1 ), TopoDS::Vertex( V001 ), eList, nbVertexInWires );
1350   if ( nbVertexInWires.size() != 1 || nbVertexInWires.front() != 4 ) {
1351     MESSAGE(" LoadBlockShapes() error ");
1352     return false;
1353   }
1354   for ( i = 0, elIt = eList.begin(); elIt != eList.end(); elIt++, i++ )
1355     switch ( i ) {
1356     case 0: Ex01 = *elIt; V101 = TopExp::LastVertex( *elIt, true ); break;
1357     case 1: E1y1 = *elIt; V111 = TopExp::LastVertex( *elIt, true ); break;
1358     case 2: Ex11 = *elIt; V011 = TopExp::LastVertex( *elIt, true ); break;
1359     case 3: E0y1 = *elIt; break;
1360     default:;
1361     }
1362   if ( i != 4 || Ex01.IsNull() || E1y1.IsNull() || Ex11.IsNull() || E0y1.IsNull() ) {
1363     MESSAGE(" LoadBlockShapes() error, eList.size()=" << eList.size());
1364     return false;
1365   }
1366
1367   // swap Fx0z and F0yz if necessary
1368   TopExp_Explorer exp( Fx0z, TopAbs_VERTEX );
1369   for ( ; exp.More(); exp.Next() ) // Fx0z shares V101 and V100
1370     if ( V101.IsSame( exp.Current() ) || V100.IsSame( exp.Current() ))
1371       break; // V101 or V100 found
1372   if ( !exp.More() ) { // not found
1373     std::swap( Fx0z, F0yz);
1374   }
1375
1376   // find Fx1z and F1yz faces
1377   const TopTools_ListOfShape& f111List = vfMap.FindFromKey( V111 );
1378   const TopTools_ListOfShape& f110List = vfMap.FindFromKey( V110 );
1379   if (f111List.Extent() != NB_FACES_BY_VERTEX ||
1380       f110List.Extent() != NB_FACES_BY_VERTEX ) {
1381     MESSAGE(" LoadBlockShapes() " << f111List.Extent() << " " << f110List.Extent());
1382     return false;
1383   }
1384   TopTools_ListIteratorOfListOfShape f111It, f110It ( f110List);
1385   for ( j = 0 ; f110It.More(); f110It.Next(), j++ ) {
1386     if ( NB_FACES_BY_VERTEX == 6 && j % 2 ) continue; // each face encounters twice
1387     const TopoDS_Shape& F = f110It.Value();
1388     for ( i = 0, f111It.Initialize( f111List ); f111It.More(); f111It.Next(), i++ ) {
1389       if ( NB_FACES_BY_VERTEX == 6 && i % 2 ) continue; // each face encounters twice
1390       if ( F.IsSame( f111It.Value() )) { // Fx1z or F1yz found
1391         if ( Fx1z.IsNull() )
1392           Fx1z = F;
1393         else
1394           F1yz = F;
1395       }
1396     }
1397   }
1398   if ( Fx1z.IsNull() || F1yz.IsNull() ) {
1399     MESSAGE(" LoadBlockShapes() error ");
1400     return false;
1401   }
1402
1403   // swap Fx1z and F1yz if necessary
1404   for ( exp.Init( Fx1z, TopAbs_VERTEX ); exp.More(); exp.Next() )
1405     if ( V010.IsSame( exp.Current() ) || V011.IsSame( exp.Current() ))
1406       break;
1407   if ( !exp.More() ) {
1408     std::swap( Fx1z, F1yz);
1409   }
1410
1411   // find vertical edges
1412   for ( exp.Init( Fx0z, TopAbs_EDGE ); exp.More(); exp.Next() ) {
1413     const TopoDS_Edge& edge = TopoDS::Edge( exp.Current() );
1414     const TopoDS_Shape& vFirst = TopExp::FirstVertex( edge, true );
1415     if ( vFirst.IsSame( V001 ))
1416       E00z = edge;
1417     else if ( vFirst.IsSame( V100 ))
1418       E10z = edge;
1419   }
1420   if ( E00z.IsNull() || E10z.IsNull() ) {
1421     MESSAGE(" LoadBlockShapes() error ");
1422     return false;
1423   }
1424   for ( exp.Init( Fx1z, TopAbs_EDGE ); exp.More(); exp.Next() ) {
1425     const TopoDS_Edge& edge = TopoDS::Edge( exp.Current() );
1426     const TopoDS_Shape& vFirst = TopExp::FirstVertex( edge, true );
1427     if ( vFirst.IsSame( V111 ))
1428       E11z = edge;
1429     else if ( vFirst.IsSame( V010 ))
1430       E01z = edge;
1431   }
1432   if ( E01z.IsNull() || E11z.IsNull() ) {
1433     MESSAGE(" LoadBlockShapes() error ");
1434     return false;
1435   }
1436
1437   // load shapes in theShapeIDMap
1438
1439   theShapeIDMap.Clear();
1440
1441   theShapeIDMap.Add(V000.Oriented( TopAbs_FORWARD ));
1442   theShapeIDMap.Add(V100.Oriented( TopAbs_FORWARD ));
1443   theShapeIDMap.Add(V010.Oriented( TopAbs_FORWARD ));
1444   theShapeIDMap.Add(V110.Oriented( TopAbs_FORWARD ));
1445   theShapeIDMap.Add(V001.Oriented( TopAbs_FORWARD ));
1446   theShapeIDMap.Add(V101.Oriented( TopAbs_FORWARD ));
1447   theShapeIDMap.Add(V011.Oriented( TopAbs_FORWARD ));
1448   theShapeIDMap.Add(V111.Oriented( TopAbs_FORWARD ));
1449
1450   theShapeIDMap.Add(Ex00);
1451   theShapeIDMap.Add(Ex10);
1452   theShapeIDMap.Add(Ex01);
1453   theShapeIDMap.Add(Ex11);
1454
1455   theShapeIDMap.Add(E0y0);
1456   theShapeIDMap.Add(E1y0);
1457   theShapeIDMap.Add(E0y1);
1458   theShapeIDMap.Add(E1y1);
1459
1460   theShapeIDMap.Add(E00z);
1461   theShapeIDMap.Add(E10z);
1462   theShapeIDMap.Add(E01z);
1463   theShapeIDMap.Add(E11z);
1464
1465   theShapeIDMap.Add(Fxy0);
1466   theShapeIDMap.Add(Fxy1);
1467   theShapeIDMap.Add(Fx0z);
1468   theShapeIDMap.Add(Fx1z);
1469   theShapeIDMap.Add(F0yz);
1470   theShapeIDMap.Add(F1yz);
1471
1472   theShapeIDMap.Add(theShell);
1473
1474   return true;
1475 }
1476
1477 //================================================================================
1478 /*!
1479  * \brief Initialize block geometry with shapes from theShapeIDMap
1480   * \param theShapeIDMap - map of block subshapes
1481   * \retval bool - is a success
1482  */
1483 //================================================================================
1484
1485 bool SMESH_Block::LoadBlockShapes(const TopTools_IndexedMapOfOrientedShape& theShapeIDMap)
1486 {
1487   init();
1488
1489   // store shapes geometry
1490   for ( int shapeID = 1; shapeID < theShapeIDMap.Extent(); shapeID++ )
1491   {
1492     const TopoDS_Shape& S = theShapeIDMap( shapeID );
1493     switch ( S.ShapeType() )
1494     {
1495     case TopAbs_VERTEX: {
1496
1497       if ( !IsVertexID( ID_V111 )) return false;
1498       myPnt[ shapeID - ID_V000 ] = BRep_Tool::Pnt( TopoDS::Vertex( S )).XYZ();
1499       break;
1500     }
1501     case TopAbs_EDGE: {
1502
1503       if ( !IsEdgeID( shapeID )) return false;
1504       const TopoDS_Edge& edge = TopoDS::Edge( S );
1505       TEdge& tEdge = myEdge[ shapeID - ID_FirstE ];
1506       tEdge.Set( shapeID,
1507                  new BRepAdaptor_Curve( edge ),
1508                  IsForwardEdge( edge, theShapeIDMap ));
1509       break;
1510     }
1511     case TopAbs_FACE: {
1512
1513       if ( !LoadFace( TopoDS::Face( S ), shapeID, theShapeIDMap ))
1514         return false;
1515       break;
1516     }
1517     default: break;
1518     }
1519   } // loop on shapes in theShapeIDMap
1520
1521   return true;
1522 }
1523
1524 //================================================================================
1525 /*!
1526  * \brief Load face geometry
1527   * \param theFace - face
1528   * \param theFaceID - face in-block ID
1529   * \param theShapeIDMap - map of block subshapes
1530   * \retval bool - is a success
1531  *
1532  * It is enough to compute params or coordinates on the face.
1533  * Face subshapes must be loaded into theShapeIDMap before
1534  */
1535 //================================================================================
1536
1537 bool SMESH_Block::LoadFace(const TopoDS_Face& theFace,
1538                            const int          theFaceID,
1539                            const TopTools_IndexedMapOfOrientedShape& theShapeIDMap)
1540 {
1541   if ( !IsFaceID( theFaceID ) ) return false;
1542   // pcurves
1543   Adaptor2d_Curve2d* c2d[4];
1544   bool isForward[4];
1545   vector< int > edgeIdVec;
1546   GetFaceEdgesIDs( theFaceID, edgeIdVec );
1547   for ( int iE = 0; iE < edgeIdVec.size(); iE++ ) // loop on 4 edges
1548   {
1549     if ( edgeIdVec[ iE ] > theShapeIDMap.Extent() )
1550       return false;
1551     const TopoDS_Edge& edge = TopoDS::Edge( theShapeIDMap( edgeIdVec[ iE ]));
1552     c2d[ iE ] = new BRepAdaptor_Curve2d( edge, theFace );
1553     isForward[ iE ] = IsForwardEdge( edge, theShapeIDMap );
1554   }
1555   TFace& tFace = myFace[ theFaceID - ID_FirstF ];
1556   tFace.Set( theFaceID, new BRepAdaptor_Surface( theFace ), c2d, isForward );
1557   return true;
1558 }
1559
1560 //================================================================================
1561 /*!
1562  * \brief/ Insert theShape into theShapeIDMap with theShapeID
1563   * \param theShape - shape to insert
1564   * \param theShapeID - shape in-block ID
1565   * \param theShapeIDMap - map of block subshapes
1566  */
1567 //================================================================================
1568
1569 bool SMESH_Block::Insert(const TopoDS_Shape& theShape,
1570                          const int           theShapeID,
1571                          TopTools_IndexedMapOfOrientedShape& theShapeIDMap)
1572 {
1573   if ( !theShape.IsNull() && theShapeID > 0 )
1574   {
1575     if ( theShapeIDMap.Contains( theShape ))
1576       return ( theShapeIDMap.FindIndex( theShape ) == theShapeID );
1577
1578     if ( theShapeID <= theShapeIDMap.Extent() ) {
1579         theShapeIDMap.Substitute( theShapeID, theShape );
1580     }
1581     else {
1582       while ( theShapeIDMap.Extent() < theShapeID - 1 ) {
1583         TopoDS_Compound comp;
1584         BRep_Builder().MakeCompound( comp );
1585         theShapeIDMap.Add( comp );
1586       }
1587       theShapeIDMap.Add( theShape );
1588     }
1589     return true;
1590   }
1591   return false;
1592 }
1593
1594 //=======================================================================
1595 //function : GetFaceEdgesIDs
1596 //purpose  : return edges IDs in the order u0, u1, 0v, 1v
1597 //           u0 means "|| u, v == 0"
1598 //=======================================================================
1599
1600 void SMESH_Block::GetFaceEdgesIDs (const int faceID, vector< int >& edgeVec )
1601 {
1602   edgeVec.resize( 4 );
1603   switch ( faceID ) {
1604   case ID_Fxy0:
1605     edgeVec[ 0 ] = ID_Ex00;
1606     edgeVec[ 1 ] = ID_Ex10;
1607     edgeVec[ 2 ] = ID_E0y0;
1608     edgeVec[ 3 ] = ID_E1y0;
1609     break;
1610   case ID_Fxy1:
1611     edgeVec[ 0 ] = ID_Ex01;
1612     edgeVec[ 1 ] = ID_Ex11;
1613     edgeVec[ 2 ] = ID_E0y1;
1614     edgeVec[ 3 ] = ID_E1y1;
1615     break;
1616   case ID_Fx0z:
1617     edgeVec[ 0 ] = ID_Ex00;
1618     edgeVec[ 1 ] = ID_Ex01;
1619     edgeVec[ 2 ] = ID_E00z;
1620     edgeVec[ 3 ] = ID_E10z;
1621     break;
1622   case ID_Fx1z:
1623     edgeVec[ 0 ] = ID_Ex10;
1624     edgeVec[ 1 ] = ID_Ex11;
1625     edgeVec[ 2 ] = ID_E01z;
1626     edgeVec[ 3 ] = ID_E11z;
1627     break;
1628   case ID_F0yz:
1629     edgeVec[ 0 ] = ID_E0y0;
1630     edgeVec[ 1 ] = ID_E0y1;
1631     edgeVec[ 2 ] = ID_E00z;
1632     edgeVec[ 3 ] = ID_E01z;
1633     break;
1634   case ID_F1yz:
1635     edgeVec[ 0 ] = ID_E1y0;
1636     edgeVec[ 1 ] = ID_E1y1;
1637     edgeVec[ 2 ] = ID_E10z;
1638     edgeVec[ 3 ] = ID_E11z;
1639     break;
1640   default:
1641     MESSAGE(" GetFaceEdgesIDs(), wrong face ID: " << faceID );
1642   }
1643 }
1644
1645 //=======================================================================
1646 //function : GetEdgeVertexIDs
1647 //purpose  : return vertex IDs of an edge
1648 //=======================================================================
1649
1650 void SMESH_Block::GetEdgeVertexIDs (const int edgeID, vector< int >& vertexVec )
1651 {
1652   vertexVec.resize( 2 );
1653   switch ( edgeID ) {
1654
1655   case ID_Ex00:
1656     vertexVec[ 0 ] = ID_V000;
1657     vertexVec[ 1 ] = ID_V100;
1658     break;
1659   case ID_Ex10:
1660     vertexVec[ 0 ] = ID_V010;
1661     vertexVec[ 1 ] = ID_V110;
1662     break;
1663   case ID_Ex01:
1664     vertexVec[ 0 ] = ID_V001;
1665     vertexVec[ 1 ] = ID_V101;
1666     break;
1667   case ID_Ex11:
1668     vertexVec[ 0 ] = ID_V011;
1669     vertexVec[ 1 ] = ID_V111;
1670     break;
1671
1672   case ID_E0y0:
1673     vertexVec[ 0 ] = ID_V000;
1674     vertexVec[ 1 ] = ID_V010;
1675     break;
1676   case ID_E1y0:
1677     vertexVec[ 0 ] = ID_V100;
1678     vertexVec[ 1 ] = ID_V110;
1679     break;
1680   case ID_E0y1:
1681     vertexVec[ 0 ] = ID_V001;
1682     vertexVec[ 1 ] = ID_V011;
1683     break;
1684   case ID_E1y1:
1685     vertexVec[ 0 ] = ID_V101;
1686     vertexVec[ 1 ] = ID_V111;
1687     break;
1688
1689   case ID_E00z:
1690     vertexVec[ 0 ] = ID_V000;
1691     vertexVec[ 1 ] = ID_V001;
1692     break;
1693   case ID_E10z:
1694     vertexVec[ 0 ] = ID_V100;
1695     vertexVec[ 1 ] = ID_V101;
1696     break;
1697   case ID_E01z:
1698     vertexVec[ 0 ] = ID_V010;
1699     vertexVec[ 1 ] = ID_V011;
1700     break;
1701   case ID_E11z:
1702     vertexVec[ 0 ] = ID_V110;
1703     vertexVec[ 1 ] = ID_V111;
1704     break;
1705   default:
1706     vertexVec.resize(0);
1707     MESSAGE(" GetEdgeVertexIDs(), wrong edge ID: " << edgeID );
1708   }
1709 }