src/SMESH/SMESH_Block.cxx

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