src/SMESHUtils/SMESH_Block.cxx

   1 // Copyright (C) 2007-2013  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
  23 // File      : SMESH_Pattern.hxx
  24 // Created   : Mon Aug  2 10:30:00 2004
  25 // Author    : Edward AGAPOV (eap)
  26 //
  27 #include "SMESH_Block.hxx"
  28
  29 #include <BRepAdaptor_Curve.hxx>
  30 #include <BRepAdaptor_Curve2d.hxx>
  31 #include <BRepAdaptor_Surface.hxx>
  32 #include <BRepTools.hxx>
  33 #include <BRepTools_WireExplorer.hxx>
  34 #include <BRep_Builder.hxx>
  35 #include <BRep_Tool.hxx>
  36 #include <Bnd_B2d.hxx>
  37 #include <Bnd_Box.hxx>
  38 #include <Extrema_ExtPC.hxx>
  39 #include <Extrema_ExtPS.hxx>
  40 #include <Extrema_POnCurv.hxx>
  41 #include <Extrema_POnSurf.hxx>
  42 #include <Geom2d_Curve.hxx>
  43 #include <ShapeAnalysis.hxx>
  44 #include <TopExp_Explorer.hxx>
  45 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
  46 #include <TopTools_ListIteratorOfListOfShape.hxx>
  47 #include <TopTools_ListOfShape.hxx>
  48 #include <TopoDS.hxx>
  49 #include <TopoDS_Compound.hxx>
  50 #include <TopoDS_Face.hxx>
  51 #include <TopoDS_Iterator.hxx>
  52 #include <TopoDS_Wire.hxx>
  53 #include <gp_Trsf.hxx>
  54 #include <gp_Vec.hxx>
  55 #include <math_FunctionSetRoot.hxx>
  56 #include <math_Matrix.hxx>
  57 #include <math_Vector.hxx>
  58
  59 #include "SMDS_MeshNode.hxx"
  60 #include "SMDS_MeshVolume.hxx"
  61 #include "SMDS_VolumeTool.hxx"
  62 #include "utilities.h"
  63
  64 #include <list>
  65 #include <limits>
  66
  67 using namespace std;
  68
  69 //#define DEBUG_PARAM_COMPUTE
  70
  71 //================================================================================
  72 /*!
  73  * \brief Set edge data
  74   * \param edgeID - block sub-shape ID
  75   * \param curve - edge geometry
  76   * \param isForward - is curve orientation coincides with edge orientation in the block
  77  */
  78 //================================================================================
  79
  80 void SMESH_Block::TEdge::Set( const int edgeID, Adaptor3d_Curve* curve, const bool isForward )
  81 {
  82   myCoordInd = SMESH_Block::GetCoordIndOnEdge( edgeID );
  83   if ( myC3d ) delete myC3d;
  84   myC3d = curve;
  85   myFirst = curve->FirstParameter();
  86   myLast = curve->LastParameter();
  87   if ( !isForward )
  88     std::swap( myFirst, myLast );
  89 }
  90
  91 //================================================================================
  92 /*!
  93  * \brief Set coordinates of nodes at edge ends to work with mesh block
  94   * \param edgeID - block sub-shape ID
  95   * \param node1 - coordinates of node with lower ID
  96   * \param node2 - coordinates of node with upper ID
  97  */
  98 //================================================================================
  99
 100 void SMESH_Block::TEdge::Set( const int edgeID, const gp_XYZ& node1, const gp_XYZ& node2 )
 101 {
 102   myCoordInd = SMESH_Block::GetCoordIndOnEdge( edgeID );
 103   myNodes[ 0 ] = node1;
 104   myNodes[ 1 ] = node2;
 105
 106   if ( myC3d ) delete myC3d;
 107   myC3d = 0;
 108 }
 109
 110 //=======================================================================
 111 //function : SMESH_Block::TEdge::GetU
 112 //purpose  :
 113 //=======================================================================
 114
 115 double SMESH_Block::TEdge::GetU( const gp_XYZ& theParams ) const
 116 {
 117   double u = theParams.Coord( myCoordInd );
 118   if ( !myC3d ) // if mesh block
 119     return u;
 120   return ( 1 - u ) * myFirst + u * myLast;
 121 }
 122
 123 //=======================================================================
 124 //function : SMESH_Block::TEdge::Point
 125 //purpose  :
 126 //=======================================================================
 127
 128 gp_XYZ SMESH_Block::TEdge::Point( const gp_XYZ& theParams ) const
 129 {
 130   double u = GetU( theParams );
 131   if ( myC3d ) return myC3d->Value( u ).XYZ();
 132   // mesh block
 133   return myNodes[0] * ( 1 - u ) + myNodes[1] * u;
 134 }
 135
 136 //================================================================================
 137 /*!
 138  * \brief Destructor
 139  */
 140 //================================================================================
 141
 142 SMESH_Block::TEdge::~TEdge()
 143 {
 144   if ( myC3d ) delete myC3d;
 145 }
 146
 147 //================================================================================
 148 /*!
 149  * \brief Set face data
 150   * \param faceID - block sub-shape ID
 151   * \param S - face surface geometry
 152   * \param c2d - 4 pcurves in the order as returned by GetFaceEdgesIDs(faceID)
 153   * \param isForward - orientation of pcurves comparing with block edge direction
 154  */
 155 //================================================================================
 156
 157 void SMESH_Block::TFace::Set( const int          faceID,
 158                               Adaptor3d_Surface* S,
 159                               Adaptor2d_Curve2d* c2D[4],
 160                               const bool         isForward[4] )
 161 {
 162   if ( myS ) delete myS;
 163   myS = S;
 164   // pcurves
 165   vector< int > edgeIdVec;
 166   GetFaceEdgesIDs( faceID, edgeIdVec );
 167   for ( int iE = 0; iE < edgeIdVec.size(); iE++ ) // loop on 4 edges
 168   {
 169     myCoordInd[ iE ] = GetCoordIndOnEdge( edgeIdVec[ iE ] );
 170     if ( myC2d[ iE ]) delete myC2d[ iE ];
 171     myC2d[ iE ] = c2D[ iE ];
 172     myFirst[ iE ] = myC2d[ iE ]->FirstParameter();
 173     myLast [ iE ] = myC2d[ iE ]->LastParameter();
 174     if ( !isForward[ iE ])
 175       std::swap( myFirst[ iE ], myLast[ iE ] );
 176   }
 177   // 2d corners
 178   myCorner[ 0 ] = myC2d[ 0 ]->Value( myFirst[0] ).XY();
 179   myCorner[ 1 ] = myC2d[ 0 ]->Value( myLast[0] ).XY();
 180   myCorner[ 2 ] = myC2d[ 1 ]->Value( myLast[1] ).XY();
 181   myCorner[ 3 ] = myC2d[ 1 ]->Value( myFirst[1] ).XY();
 182 }
 183
 184 //================================================================================
 185 /*!
 186  * \brief Set face data to work with mesh block
 187   * \param faceID - block sub-shape ID
 188   * \param edgeU0 - filled data of edge u0 = GetFaceEdgesIDs(faceID)[ 0 ]
 189   * \param edgeU1 - filled data of edge u1 = GetFaceEdgesIDs(faceID)[ 1 ]
 190  */
 191 //================================================================================
 192
 193 void SMESH_Block::TFace::Set( const int faceID, const TEdge& edgeU0, const TEdge& edgeU1 )
 194 {
 195   vector< int > edgeIdVec;
 196   GetFaceEdgesIDs( faceID, edgeIdVec );
 197   myNodes[ 0 ] = edgeU0.NodeXYZ( 1 );
 198   myNodes[ 1 ] = edgeU0.NodeXYZ( 0 );
 199   myNodes[ 2 ] = edgeU1.NodeXYZ( 0 );
 200   myNodes[ 3 ] = edgeU1.NodeXYZ( 1 );
 201   myCoordInd[ 0 ] = GetCoordIndOnEdge( edgeIdVec[ 0 ] );
 202   myCoordInd[ 1 ] = GetCoordIndOnEdge( edgeIdVec[ 1 ] );
 203   myCoordInd[ 2 ] = GetCoordIndOnEdge( edgeIdVec[ 2 ] );
 204   myCoordInd[ 3 ] = GetCoordIndOnEdge( edgeIdVec[ 3 ] );
 205   if ( myS ) delete myS;
 206   myS = 0;
 207 }
 208
 209 //================================================================================
 210 /*!
 211  * \brief Destructor
 212  */
 213 //================================================================================
 214
 215 SMESH_Block::TFace::~TFace()
 216 {
 217   if ( myS ) delete myS;
 218   for ( int i = 0 ; i < 4; ++i )
 219     if ( myC2d[ i ]) delete myC2d[ i ];
 220 }
 221
 222 //=======================================================================
 223 //function : SMESH_Block::TFace::GetCoefs
 224 //purpose  : return coefficients for addition of [0-3]-th edge and vertex
 225 //=======================================================================
 226
 227 void SMESH_Block::TFace::GetCoefs(int           iE,
 228                                   const gp_XYZ& theParams,
 229                                   double&       Ecoef,
 230                                   double&       Vcoef ) const
 231 {
 232   double dU = theParams.Coord( GetUInd() );
 233   double dV = theParams.Coord( GetVInd() );
 234   switch ( iE ) {
 235   case 0:
 236     Ecoef = ( 1 - dV ); // u0
 237     Vcoef = ( 1 - dU ) * ( 1 - dV ); break; // 00
 238   case 1:
 239     Ecoef = dV; // u1
 240     Vcoef = dU * ( 1 - dV ); break; // 10
 241   case 2:
 242     Ecoef = ( 1 - dU ); // 0v
 243     Vcoef = dU * dV  ; break; // 11
 244   case 3:
 245     Ecoef = dU  ; // 1v
 246     Vcoef = ( 1 - dU ) * dV  ; break; // 01
 247   default: ASSERT(0);
 248   }
 249 }
 250
 251 //=======================================================================
 252 //function : SMESH_Block::TFace::GetUV
 253 //purpose  :
 254 //=======================================================================
 255
 256 gp_XY SMESH_Block::TFace::GetUV( const gp_XYZ& theParams ) const
 257 {
 258   gp_XY uv(0.,0.);
 259   for ( int iE = 0; iE < 4; iE++ ) // loop on 4 edges
 260   {
 261     double Ecoef = 0, Vcoef = 0;
 262     GetCoefs( iE, theParams, Ecoef, Vcoef );
 263     // edge addition
 264     double u = theParams.Coord( myCoordInd[ iE ] );
 265     u = ( 1 - u ) * myFirst[ iE ] + u * myLast[ iE ];
 266     uv += Ecoef * myC2d[ iE ]->Value( u ).XY();
 267     // corner addition
 268     uv -= Vcoef * myCorner[ iE ];
 269   }
 270   return uv;
 271 }
 272
 273 //=======================================================================
 274 //function : SMESH_Block::TFace::Point
 275 //purpose  :
 276 //=======================================================================
 277
 278 gp_XYZ SMESH_Block::TFace::Point( const gp_XYZ& theParams ) const
 279 {
 280   gp_XYZ p(0.,0.,0.);
 281   if ( !myS ) // if mesh block
 282   {
 283     for ( int iE = 0; iE < 4; iE++ ) // loop on 4 edges
 284     {
 285       double Ecoef = 0, Vcoef = 0;
 286       GetCoefs( iE, theParams, Ecoef, Vcoef );
 287       // edge addition
 288       double u = theParams.Coord( myCoordInd[ iE ] );
 289       int i1 = 0, i2 = 1;
 290       switch ( iE ) {
 291       case 1: i1 = 3; i2 = 2; break;
 292       case 2: i1 = 1; i2 = 2; break;
 293       case 3: i1 = 0; i2 = 3; break;
 294       }
 295       p += Ecoef * ( myNodes[ i1 ] * ( 1 - u ) + myNodes[ i2 ] * u );
 296       // corner addition
 297       p -= Vcoef * myNodes[ iE ];
 298     }
 299
 300   }
 301   else // shape block
 302   {
 303     gp_XY uv = GetUV( theParams );
 304     p = myS->Value( uv.X(), uv.Y() ).XYZ();
 305   }
 306   return p;
 307 }
 308
 309
 310 namespace
 311 {
 312   bool isPntInTria( const gp_XY& p, const gp_XY& t0, const gp_XY& t1, const gp_XY& t2  )
 313   {
 314     const double // matrix 2x2
 315       T11 = t0.X()-t2.X(), T12 = t1.X()-t2.X(),
 316       T21 = t0.Y()-t2.Y(), T22 = t1.Y()-t2.Y();
 317     const double Tdet = T11*T22 - T12*T21; // matrix determinant
 318     if ( Abs( Tdet ) < std::numeric_limits<double>::min() )
 319       return false;
 320     // matrix inverse
 321     const double t11 = T22, t12 = -T12, t21 = -T21, t22 = T11;
 322     // vector
 323     const double r11 = p.X()-t2.X(), r12 = p.Y()-t2.Y();
 324     // barycentric coordinates: mutiply matrix by vector
 325     const double bc0 = (t11 * r11 + t12 * r12)/Tdet;
 326     const double bc1 = (t21 * r11 + t22 * r12)/Tdet;
 327     return ( bc0 >= 0. && bc1 >= 0. && bc0 + bc1 <= 1. );
 328   }
 329
 330   bool isPntInQuad( const gp_XY& p,
 331                     const gp_XY& q0, const gp_XY& q1, const gp_XY& q2, const gp_XY& q3 )
 332   {
 333     const int in1 = isPntInTria( p, q0, q1, q2 );
 334     const int in2 = isPntInTria( p, q0, q2, q3 );
 335     return in1 + in2 == 1;
 336   }
 337 }
 338
 339 //=======================================================================
 340 //function : IsUVInQuad
 341 //purpose  : Checks if UV is in a quardilateral defined by 4 nornalized points
 342 //=======================================================================
 343
 344 bool SMESH_Block::TFace::IsUVInQuad( const gp_XY& uv,
 345                                      const gp_XYZ& param0, const gp_XYZ& param1,
 346                                      const gp_XYZ& param2, const gp_XYZ& param3 ) const
 347 {
 348   gp_XY q0 = GetUV( param0 );
 349   gp_XY q1 = GetUV( param1 );
 350   gp_XY q2 = GetUV( param2 );
 351   gp_XY q3 = GetUV( param3 );
 352   return isPntInQuad( uv, q0,q1,q2,q3);
 353 }
 354
 355 //=======================================================================
 356 //function : GetUVRange
 357 //purpose  : returns UV range of the face
 358 //=======================================================================
 359
 360 gp_XY SMESH_Block::TFace::GetUVRange() const
 361 {
 362   if ( !myS ) return gp_XY(1.,1.);
 363
 364   Bnd_B2d bb;
 365   for ( int iE = 0; iE < 4; ++iE )
 366   {
 367     //TColStd_Array1OfReal T(1,
 368   }
 369   return bb.CornerMax() - bb.CornerMin();
 370 }
 371
 372 //=======================================================================
 373 //function : GetShapeCoef
 374 //purpose  :
 375 //=======================================================================
 376
 377 double* SMESH_Block::GetShapeCoef (const int theShapeID)
 378 {
 379   static double shapeCoef[][3] = {
 380     //    V000,        V100,        V010,         V110
 381     { -1,-1,-1 }, {  1,-1,-1 }, { -1, 1,-1 }, {  1, 1,-1 },
 382     //    V001,        V101,        V011,         V111,
 383     { -1,-1, 1 }, {  1,-1, 1 }, { -1, 1, 1 }, {  1, 1, 1 },
 384     //    Ex00,        Ex10,        Ex01,         Ex11,
 385     {  0,-1,-1 }, {  0, 1,-1 }, {  0,-1, 1 }, {  0, 1, 1 },
 386     //    E0y0,        E1y0,        E0y1,         E1y1,
 387     { -1, 0,-1 }, {  1, 0,-1 }, { -1, 0, 1 }, {  1, 0, 1 },
 388     //    E00z,        E10z,        E01z,         E11z,
 389     { -1,-1, 0 }, {  1,-1, 0 }, { -1, 1, 0 }, {  1, 1, 0 },
 390     //    Fxy0,        Fxy1,        Fx0z,         Fx1z,         F0yz,           F1yz,
 391     {  0, 0,-1 }, {  0, 0, 1 }, {  0,-1, 0 }, {  0, 1, 0 }, { -1, 0, 0 }, {  1, 0, 0 },
 392     // ID_Shell
 393     {  0, 0, 0 }
 394   };
 395   if ( theShapeID < ID_V000 || theShapeID > ID_F1yz )
 396     return shapeCoef[ ID_Shell - 1 ];
 397
 398   return shapeCoef[ theShapeID - 1 ];
 399 }
 400
 401 //=======================================================================
 402 //function : ShellPoint
 403 //purpose  : return coordinates of a point in shell
 404 //=======================================================================
 405
 406 bool SMESH_Block::ShellPoint( const gp_XYZ& theParams, gp_XYZ& thePoint ) const
 407 {
 408   thePoint.SetCoord( 0., 0., 0. );
 409   for ( int shapeID = ID_V000; shapeID < ID_Shell; shapeID++ )
 410   {
 411     // coef
 412     double* coefs = GetShapeCoef( shapeID );
 413     double k = 1;
 414     for ( int iCoef = 0; iCoef < 3; iCoef++ ) {
 415       if ( coefs[ iCoef ] != 0 ) {
 416         if ( coefs[ iCoef ] < 0 )
 417           k *= ( 1. - theParams.Coord( iCoef + 1 ));
 418         else
 419           k *= theParams.Coord( iCoef + 1 );
 420       }
 421     }
 422     // add point on a shape
 423     if ( fabs( k ) > DBL_MIN )
 424     {
 425       gp_XYZ Ps;
 426       if ( shapeID < ID_Ex00 ) // vertex
 427         VertexPoint( shapeID, Ps );
 428       else if ( shapeID < ID_Fxy0 ) { // edge
 429         EdgePoint( shapeID, theParams, Ps );
 430         k = -k;
 431       } else // face
 432         FacePoint( shapeID, theParams, Ps );
 433
 434       thePoint += k * Ps;
 435     }
 436   }
 437   return true;
 438 }
 439
 440 //=======================================================================
 441 //function : ShellPoint
 442 //purpose  : computes coordinates of a point in shell by points on sub-shapes;
 443 //           thePointOnShape[ subShapeID ] must be a point on a sub-shape
 444 //=======================================================================
 445
 446 bool SMESH_Block::ShellPoint(const gp_XYZ&         theParams,
 447                              const vector<gp_XYZ>& thePointOnShape,
 448                              gp_XYZ&               thePoint )
 449 {
 450   if ( thePointOnShape.size() < ID_F1yz )
 451     return false;
 452
 453   const double x = theParams.X(), y = theParams.Y(), z = theParams.Z();
 454   const double x1 = 1. - x,       y1 = 1. - y,       z1 = 1. - z;
 455   const vector<gp_XYZ>& p = thePointOnShape;
 456
 457   thePoint =
 458     x1 * p[ID_F0yz] + x * p[ID_F1yz] +
 459     y1 * p[ID_Fx0z] + y * p[ID_Fx1z] +
 460     z1 * p[ID_Fxy0] + z * p[ID_Fxy1] +
 461     x1 * (y1 * (z1 * p[ID_V000] + z * p[ID_V001])  +
 462           y  * (z1 * p[ID_V010] + z * p[ID_V011])) +
 463     x  * (y1 * (z1 * p[ID_V100] + z * p[ID_V101])  +
 464           y  * (z1 * p[ID_V110] + z * p[ID_V111]));
 465   thePoint -=
 466     x1 * (y1 * p[ID_E00z] + y * p[ID_E01z]) +
 467     x  * (y1 * p[ID_E10z] + y * p[ID_E11z]) +
 468     y1 * (z1 * p[ID_Ex00] + z * p[ID_Ex01]) +
 469     y  * (z1 * p[ID_Ex10] + z * p[ID_Ex11]) +
 470     z1 * (x1 * p[ID_E0y0] + x * p[ID_E1y0]) +
 471     z  * (x1 * p[ID_E0y1] + x * p[ID_E1y1]);
 472
 473   return true;
 474 }
 475
 476 //=======================================================================
 477 //function : Constructor
 478 //purpose  :
 479 //=======================================================================
 480
 481 SMESH_Block::SMESH_Block():
 482   myNbIterations(0),
 483   mySumDist(0.),
 484   myTolerance(-1.) // to be re-initialized
 485 {
 486 }
 487
 488
 489 //=======================================================================
 490 //function : NbVariables
 491 //purpose  :
 492 //=======================================================================
 493
 494 Standard_Integer SMESH_Block::NbVariables() const
 495 {
 496   return 3;
 497 }
 498
 499 //=======================================================================
 500 //function : NbEquations
 501 //purpose  :
 502 //=======================================================================
 503
 504 Standard_Integer SMESH_Block::NbEquations() const
 505 {
 506   return 1;
 507 }
 508
 509 //=======================================================================
 510 //function : Value
 511 //purpose  :
 512 //=======================================================================
 513
 514 Standard_Boolean SMESH_Block::Value(const math_Vector& theXYZ, math_Vector& theFxyz)
 515 {
 516   gp_XYZ P, params( theXYZ(1), theXYZ(2), theXYZ(3) );
 517   if ( params.IsEqual( myParam, DBL_MIN )) { // same param
 518     theFxyz( 1 ) = funcValue( myValues[ SQUARE_DIST ]);
 519   }
 520   else {
 521     ShellPoint( params, P );
 522     gp_Vec dP( P - myPoint );
 523     theFxyz(1) = funcValue( dP.SquareMagnitude() );
 524   }
 525   return true;
 526 }
 527
 528 //=======================================================================
 529 //function : Derivatives
 530 //purpose  :
 531 //=======================================================================
 532
 533 Standard_Boolean SMESH_Block::Derivatives(const math_Vector& XYZ,math_Matrix& Df)
 534 {
 535   math_Vector F(1,3);
 536   return Values(XYZ,F,Df);
 537 }
 538
 539 //=======================================================================
 540 //function : GetStateNumber
 541 //purpose  :
 542 //=======================================================================
 543
 544 Standard_Integer SMESH_Block::GetStateNumber ()
 545 {
 546   return 0; //myValues[0] < 1e-1;
 547 }
 548
 549 //=======================================================================
 550 //function : Values
 551 //purpose  :
 552 //=======================================================================
 553
 554 Standard_Boolean SMESH_Block::Values(const math_Vector& theXYZ,
 555                                      math_Vector&       theFxyz,
 556                                      math_Matrix&       theDf)
 557 {
 558   gp_XYZ P, params( theXYZ(1), theXYZ(2), theXYZ(3) );
 559   if ( params.IsEqual( myParam, DBL_MIN )) { // same param
 560     theFxyz( 1 )      = funcValue( myValues[ SQUARE_DIST ] );
 561     theDf( 1, DRV_1 ) = myValues[ DRV_1 ];
 562     theDf( 1, DRV_2 ) = myValues[ DRV_2 ];
 563     theDf( 1, DRV_3 ) = myValues[ DRV_3 ];
 564     return true;
 565   }
 566 #ifdef DEBUG_PARAM_COMPUTE
 567   MESSAGE ( "PARAM GUESS: " << params.X() << " "<< params.Y() << " "<< params.X() );
 568   myNbIterations++; // how many times call ShellPoint()
 569 #endif
 570   ShellPoint( params, P );
 571
 572   gp_Vec dP( myPoint, P );
 573   double sqDist = dP.SquareMagnitude();
 574   theFxyz(1) = funcValue( sqDist );
 575
 576   if ( sqDist < myTolerance * myTolerance ) { // a solution found
 577     myParam = params;
 578     myValues[ SQUARE_DIST ] = sqDist;
 579     theFxyz(1)  = theDf( 1,1 ) = theDf( 1,2 ) = theDf( 1,3 ) = 0;
 580     return true;
 581   }
 582
 583   if ( sqDist < myValues[ SQUARE_DIST ] ) // a better guess
 584   {
 585     // 3 partial derivatives
 586     gp_Vec drv[ 3 ]; // where we move with a small step in each direction
 587     for ( int iP = 1; iP <= 3; iP++ ) {
 588       if ( iP == myFaceIndex ) {
 589         drv[ iP - 1 ] = gp_Vec(0,0,0);
 590         continue;
 591       }
 592       gp_XYZ Pi;
 593       bool onEdge = ( theXYZ( iP ) + 0.001 > 1. );
 594       if ( onEdge )
 595         params.SetCoord( iP, theXYZ( iP ) - 0.001 );
 596       else
 597         params.SetCoord( iP, theXYZ( iP ) + 0.001 );
 598       ShellPoint( params, Pi );
 599       params.SetCoord( iP, theXYZ( iP ) ); // restore params
 600       gp_Vec dPi ( P, Pi );
 601       if ( onEdge ) dPi *= -1.;
 602       double mag = dPi.Magnitude();
 603       if ( mag > DBL_MIN )
 604         dPi /= mag;
 605       drv[ iP - 1 ] = dPi;
 606     }
 607     for ( int iP = 0; iP < 3; iP++ ) {
 608 #if 1
 609       theDf( 1, iP + 1 ) = dP * drv[iP];
 610 #else
 611       // Distance from P to plane passing through myPoint and defined
 612       // by the 2 other derivative directions:
 613       // like IntAna_IntConicQuad::Perform (const gp_Lin& L, const gp_Pln& P)
 614       // where L is (P -> myPoint), P is defined by the 2 other derivative direction
 615       int iPrev = ( iP ? iP - 1 : 2 );
 616       int iNext = ( iP == 2 ? 0 : iP + 1 );
 617       gp_Vec plnNorm = drv[ iPrev ].Crossed( drv [ iNext ] );
 618       double Direc = plnNorm * drv[ iP ];
 619       if ( Abs(Direc) <= DBL_MIN )
 620         theDf( 1, iP + 1 ) = dP * drv[ iP ];
 621       else {
 622         double Dis = plnNorm * P - plnNorm * myPoint;
 623         theDf( 1, iP + 1 ) = Dis/Direc;
 624       }
 625 #endif
 626     }
 627 #ifdef DEBUG_PARAM_COMPUTE
 628     MESSAGE ( "F = " << theFxyz(1) << " DRV: " << theDf(1,1) << " " << theDf(1,2) << " " << theDf(1,3) );
 629     myNbIterations +=3; // how many times call ShellPoint()
 630 #endif
 631
 632     // store better values
 633     myParam              = params;
 634     myValues[SQUARE_DIST]= sqDist;
 635     myValues[DRV_1]      = theDf(1,DRV_1);
 636     myValues[DRV_2]      = theDf(1,DRV_2);
 637     myValues[DRV_3]      = theDf(1,DRV_3);
 638   }
 639
 640   return true;
 641 }
 642
 643 //============================================================================
 644 //function : computeParameters
 645 //purpose  : compute point parameters in the block using math_FunctionSetRoot
 646 //============================================================================
 647
 648 bool SMESH_Block::computeParameters(const gp_Pnt& thePoint,
 649                                     gp_XYZ&       theParams,
 650                                     const gp_XYZ& theParamsHint,
 651                                     int           theShapeID)
 652 {
 653   myPoint = thePoint.XYZ();
 654
 655   myParam.SetCoord( -1,-1,-1 );
 656   myValues[ SQUARE_DIST ] = 1e100;
 657
 658   math_Vector low  ( 1, 3, 0.0 );
 659   math_Vector up   ( 1, 3, 1.0 );
 660   math_Vector tol  ( 1, 3, 1e-4 );
 661   math_Vector start( 1, 3, 0.0 );
 662   start( 1 ) = theParamsHint.X();
 663   start( 2 ) = theParamsHint.Y();
 664   start( 3 ) = theParamsHint.Z();
 665
 666   math_FunctionSetRoot paramSearch( *this, tol );
 667
 668   mySquareFunc = 0; // large approaching steps
 669   //if ( hasHint ) mySquareFunc = 1; // small approaching steps
 670
 671   double loopTol = 10 * myTolerance;
 672   int nbLoops = 0;
 673   while ( distance() > loopTol && nbLoops <= 3 )
 674   {
 675     paramSearch.Perform ( *static_cast<math_FunctionSetWithDerivatives*>(this),
 676                           start, low, up );
 677     start( 1 ) = myParam.X();
 678     start( 2 ) = myParam.Y();
 679     start( 3 ) = myParam.Z();
 680     mySquareFunc = !mySquareFunc;
 681     nbLoops++;
 682   }
 683 #ifdef DEBUG_PARAM_COMPUTE
 684   mySumDist += distance();
 685   MESSAGE ( " ------ SOLUTION: ( "<< myParam.X() <<" "<< myParam.Y() <<" "<< myParam.Z() <<" )"<<endl
 686          << " ------ DIST : " << distance() << "\t Tol=" << myTolerance << "\t Nb LOOPS=" << nbLoops << endl
 687          << " ------ NB IT: " << myNbIterations << ",  SUM DIST: " << mySumDist );
 688 #endif
 689
 690   theParams = myParam;
 691
 692   if ( myFaceIndex > 0 )
 693   {
 694     theParams.SetCoord( myFaceIndex, myFaceParam );
 695     if ( distance() > loopTol )
 696       refineParametersOnFace( thePoint, theParams, theShapeID );
 697   }
 698   return true;
 699 }
 700
 701 //=======================================================================
 702 //function : ComputeParameters
 703 //purpose  : compute point parameters in the block
 704 //=======================================================================
 705
 706 bool SMESH_Block::ComputeParameters(const gp_Pnt& thePoint,
 707                                     gp_XYZ&       theParams,
 708                                     const int     theShapeID,
 709                                     const gp_XYZ& theParamsHint)
 710 {
 711   if ( VertexParameters( theShapeID, theParams ))
 712     return true;
 713
 714   if ( IsEdgeID( theShapeID )) {
 715     TEdge& e = myEdge[ theShapeID - ID_FirstE ];
 716     Adaptor3d_Curve* curve = e.GetCurve();
 717     Extrema_ExtPC anExtPC( thePoint, *curve, curve->FirstParameter(), curve->LastParameter() );
 718     int i, nb = anExtPC.IsDone() ? anExtPC.NbExt() : 0;
 719     for ( i = 1; i <= nb; i++ ) {
 720       if ( anExtPC.IsMin( i ))
 721         return EdgeParameters( theShapeID, anExtPC.Point( i ).Parameter(), theParams );
 722     }
 723     return false;
 724   }
 725
 726   const bool isOnFace = IsFaceID( theShapeID );
 727   double * coef = GetShapeCoef( theShapeID );
 728
 729   // Find the first guess paremeters
 730
 731   gp_XYZ start(0, 0, 0);
 732
 733   bool hasHint = ( 0 <= theParamsHint.X() && theParamsHint.X() <= 1 &&
 734                    0 <= theParamsHint.Y() && theParamsHint.Y() <= 1 &&
 735                    0 <= theParamsHint.Y() && theParamsHint.Y() <= 1 );
 736   if ( !hasHint && !myGridComputed )
 737   {
 738     // define the first guess by thePoint projection on lines
 739     // connecting vertices
 740     bool needGrid = false;
 741     gp_XYZ par000( 0, 0, 0 ), par111( 1, 1, 1 );
 742     double zero = DBL_MIN * DBL_MIN;
 743     for ( int iEdge = 0, iParam = 1; iParam <= 3 && !needGrid; iParam++ )
 744     {
 745       if ( isOnFace && coef[ iParam - 1 ] != 0 ) {
 746         iEdge += 4;
 747         continue;
 748       }
 749       double sumParam = 0;
 750       for ( int iE = 0; iE < 4; iE++, iEdge++ ) { // loop on 4 parallel edges
 751         gp_Pnt p0 = myEdge[ iEdge ].Point( par000 );
 752         gp_Pnt p1 = myEdge[ iEdge ].Point( par111 );
 753         gp_Vec v01( p0, p1 ), v0P( p0, thePoint );
 754         double len2 = v01.SquareMagnitude();
 755         double par = 0;
 756         if ( len2 > zero ) {
 757           par = v0P.Dot( v01 ) / len2;
 758           if ( par < 0 || par > 1 ) { // projection falls out of line ends => needGrid
 759             needGrid = true;
 760             break;
 761           }
 762         }
 763         sumParam += par;
 764       }
 765       start.SetCoord( iParam, sumParam / 4.);
 766     }
 767     if ( needGrid ) {
 768       // compute nodes of 10 x 10 x 10 grid
 769       int iNode = 0;
 770       Bnd_Box box;
 771       for ( double x = 0.05; x < 1.; x += 0.1 )
 772         for ( double y = 0.05; y < 1.; y += 0.1 )
 773           for ( double z = 0.05; z < 1.; z += 0.1 ) {
 774             TxyzPair & prmPtn = my3x3x3GridNodes[ iNode++ ];
 775             prmPtn.first.SetCoord( x, y, z );
 776             ShellPoint( prmPtn.first, prmPtn.second );
 777             box.Add( gp_Pnt( prmPtn.second ));
 778           }
 779       myGridComputed = true;
 780       if ( myTolerance < 0 )
 781         myTolerance = sqrt( box.SquareExtent() ) * 1e-5;
 782     }
 783   }
 784
 785   if ( hasHint )
 786   {
 787     start = theParamsHint;
 788   }
 789   if ( myGridComputed )
 790   {
 791     double minDist = DBL_MAX;
 792     if ( hasHint )
 793     {
 794       gp_XYZ p;
 795       if ( ShellPoint( start, p ))
 796         minDist = thePoint.SquareDistance( p );
 797     }
 798     gp_XYZ* bestParam = 0;
 799     for ( int iNode = 0; iNode < 1000; iNode++ ) {
 800       TxyzPair & prmPtn = my3x3x3GridNodes[ iNode ];
 801       double dist = ( thePoint.XYZ() - prmPtn.second ).SquareModulus();
 802       if ( dist < minDist ) {
 803         minDist = dist;
 804         bestParam = & prmPtn.first;
 805       }
 806     }
 807     if ( bestParam )
 808       start = *bestParam;
 809   }
 810
 811   myFaceIndex = -1;
 812   myFaceParam = 0.;
 813   if ( isOnFace ) {
 814     // put a point on the face
 815     for ( int iCoord = 0; iCoord < 3; iCoord++ )
 816       if ( coef[ iCoord ] ) {
 817         myFaceIndex = iCoord + 1;
 818         myFaceParam = ( coef[ iCoord ] < 0.5 ) ? 0.0 : 1.0;
 819         start.SetCoord( myFaceIndex, myFaceParam );
 820       }
 821   }
 822
 823 #ifdef DEBUG_PARAM_COMPUTE
 824   MESSAGE ( " #### POINT " <<thePoint.X()<<" "<<thePoint.Y()<<" "<<thePoint.Z()<<" ####" );
 825 #endif
 826
 827   if ( myTolerance < 0 ) myTolerance = 1e-6;
 828
 829   const double parDelta = 1e-4;
 830   const double sqTolerance = myTolerance * myTolerance;
 831
 832   gp_XYZ solution = start, params = start;
 833   double sqDistance = 1e100;
 834   int nbLoops = 0, nbGetWorst = 0;
 835
 836   while ( nbLoops <= 100 )
 837   {
 838     gp_XYZ P, Pi;
 839     ShellPoint( params, P );
 840
 841     gp_Vec dP( thePoint, P );
 842     double sqDist = dP.SquareMagnitude();
 843
 844     if ( sqDist > sqDistance ) { // solution get worse
 845       if ( ++nbGetWorst > 2 )
 846         return computeParameters( thePoint, theParams, solution, theShapeID );
 847     }
 848 #ifdef DEBUG_PARAM_COMPUTE
 849     MESSAGE ( "PARAMS: ( " << params.X() <<" "<< params.Y() <<" "<< params.Z() <<" )" );
 850     MESSAGE ( "DIST: " << sqrt( sqDist ) );
 851 #endif
 852
 853     if ( sqDist < sqDistance ) { // get better
 854       sqDistance = sqDist;
 855       solution   = params;
 856       nbGetWorst = 0;
 857       if ( sqDistance < sqTolerance ) // a solution found
 858         break;
 859     }
 860
 861         // look for a next better solution
 862     for ( int iP = 1; iP <= 3; iP++ ) {
 863       if ( iP == myFaceIndex )
 864         continue;
 865       // see where we move with a small (=parDelta) step in this direction
 866       gp_XYZ nearParams = params;
 867       bool onEdge = ( params.Coord( iP ) + parDelta > 1. );
 868       if ( onEdge )
 869         nearParams.SetCoord( iP, params.Coord( iP ) - parDelta );
 870       else
 871         nearParams.SetCoord( iP, params.Coord( iP ) + parDelta );
 872       ShellPoint( nearParams, Pi );
 873       gp_Vec dPi ( P, Pi );
 874       if ( onEdge ) dPi *= -1.;
 875       // modify a parameter
 876       double mag = dPi.Magnitude();
 877       if ( mag < DBL_MIN )
 878         continue;
 879       gp_Vec dir = dPi / mag; // dir we move modifying the parameter
 880       double dist = dir * dP; // where we should get to
 881       double dPar = dist / mag * parDelta; // predict parameter change
 882       double curPar = params.Coord( iP );
 883       double par = curPar - dPar; // new parameter value
 884       while ( par > 1 || par < 0 ) {
 885         dPar /= 2.;
 886         par = curPar - dPar;
 887       }
 888       params.SetCoord( iP, par );
 889     }
 890
 891     nbLoops++;
 892   }
 893 #ifdef DEBUG_PARAM_COMPUTE
 894   myNbIterations += nbLoops*4; // how many times ShellPoint called
 895   mySumDist += sqrt( sqDistance );
 896   MESSAGE ( " ------ SOLUTION: ( "<<solution.X()<<" "<<solution.Y()<<" "<<solution.Z()<<" )"<< std::endl
 897          << " ------ DIST : " << sqrt( sqDistance ) << "\t Tol=" << myTolerance << "\t Nb LOOPS=" << nbLoops << std::endl
 898          << " ------ NB IT: " << myNbIterations << ",  SUM DIST: " << mySumDist );
 899 #endif
 900
 901   if ( myFaceIndex > 0 )
 902     theParams.SetCoord( myFaceIndex, myFaceParam );
 903
 904   const double reachedDist = sqrt( sqDistance );
 905   // if ( reachedDist > 1000 * myTolerance &&
 906   //      computeParameters( thePoint, theParams, solution ) &&
 907   //      reachedDist > distance() )
 908   //   return true;
 909
 910   if ( reachedDist > 10 * myTolerance &&
 911        computeParameters( thePoint, theParams, solution, theShapeID ) &&
 912        reachedDist > distance() )
 913     return true;
 914
 915   theParams = solution;
 916   myValues[ SQUARE_DIST ] = sqDistance;
 917
 918   if ( reachedDist > 10 * myTolerance && myFaceIndex > 0 )
 919     refineParametersOnFace( thePoint, theParams, theShapeID );
 920
 921   return true;
 922 }
 923
 924 //================================================================================
 925 /*!
 926  * \brief Find more precise solution
 927  *  \param [in] thePoint - the point
 928  *  \param [in,out] theParams - solution to precise
 929  *  \param [in] theFaceID - FACE ID
 930  */
 931 //================================================================================
 932
 933 void SMESH_Block::refineParametersOnFace( const gp_Pnt& thePoint,
 934                                           gp_XYZ&       theParams,
 935                                           int           theFaceID )
 936 {
 937   // find UV of thePoint on the FACE
 938   Standard_Real U,V;
 939
 940   const TFace& tface = myFace[ theFaceID - ID_FirstF ];
 941   if ( !tface.Surface() ) return;
 942
 943   Extrema_ExtPS extPS( thePoint, *tface.Surface(),
 944                        tface.Surface()->UResolution( myTolerance ),
 945                        tface.Surface()->VResolution( myTolerance ));
 946   if ( !extPS.IsDone() || extPS.NbExt() < 1 )
 947     return;
 948
 949   double minDist = 1e100;
 950   for ( int i = 1; i <= extPS.NbExt(); ++i )
 951     if ( extPS.SquareDistance( i ) < minDist )
 952     {
 953       minDist = extPS.SquareDistance( i );
 954       extPS.Point( i ).Parameter( U,V );
 955     }
 956   if ( minDist > 100 * myTolerance * myTolerance )
 957     return;
 958
 959   gp_XY uv(U,V);
 960   if ( findUVByHalfDivision( thePoint, uv, tface, theParams))
 961     return;
 962
 963   int nbGetWorstLimit = 20;
 964   if ( findUVAround( thePoint, uv, tface, theParams, nbGetWorstLimit ))
 965     return;
 966
 967   double dist2, prevSolDist = distance();
 968   gp_XYZ sol = theParams;
 969   for ( double delta = 1./10; delta > 0.001; delta /= 2.5, nbGetWorstLimit *= 2 )
 970   {
 971     for ( double y = delta; y < 1.; y += delta )
 972     {
 973       sol.SetCoord( tface.GetVInd(), y );
 974       for ( double x = delta; x < 1.; x += delta )
 975       {
 976         sol.SetCoord( tface.GetUInd(), x );
 977         dist2 = thePoint.SquareDistance( tface.Point( sol ));
 978         if ( dist2 < prevSolDist * prevSolDist )
 979         {
 980           if ( findUVAround( thePoint, uv, tface, theParams, nbGetWorstLimit ))
 981             return;
 982           if ( distance() < 1000 * myTolerance )
 983             return;
 984           prevSolDist = distance();
 985         }
 986       }
 987     }
 988   }
 989 }
 990
 991 //================================================================================
 992 /*!
 993  * \brief Finds parameters corresponding to a given UV of a given face using half-division
 994  *  \param [in] theUV - the UV to locate
 995  *  \param [in] tface - the face
 996  *  \param [in,out] theParams - the starting parameters to improve
 997  *  \return bool - \c true if found solution is within myTolerance
 998  */
 999 //================================================================================
1000
1001 bool SMESH_Block::findUVByHalfDivision( const gp_Pnt&             thePoint,
1002                                         const gp_XY&              theUV,
1003                                         const SMESH_Block::TFace& tface,
1004                                         gp_XYZ&                   theParams)
1005 {
1006   int nbGetUV = 0; // just for statistics
1007
1008   // find a range of parameters including the UV
1009
1010   double xMin, xMax, yMin, yMax;
1011   //#define _DEBUG_REFINE_
1012 #ifdef _DEBUG_REFINE_
1013   cout << "SMESH_Block::refineParametersOnFace(): dividing Starts at dist " << distance()<< endl;
1014 #endif
1015   double dx = 0.1, xSol = theParams.Coord( tface.GetUInd() );
1016   double dy = 0.1, ySol = theParams.Coord( tface.GetVInd() );
1017   gp_XYZ xXYZ( 0,0,0 ); xXYZ.SetCoord( tface.GetUInd(), 1 );
1018   gp_XYZ yXYZ( 0,0,0 ); yXYZ.SetCoord( tface.GetVInd(), 1 );
1019   gp_XYZ xy0,xy1,xy2,xy3;
1020   bool isInQuad = false;
1021   while ( !isInQuad )
1022   {
1023     xMin = Max( 0., xSol - 0.5*dx ); xMax = Min( 1.0, xSol + 0.5*dx );
1024     yMin = Max( 0., ySol - 0.5*dy ); yMax = Min( 1.0, ySol + 0.5*dy );
1025     xy0.SetLinearForm( xMin, xXYZ, yMin, yXYZ );
1026     xy1.SetLinearForm( xMax, xXYZ, yMin, yXYZ );
1027     xy2.SetLinearForm( xMax, xXYZ, yMax, yXYZ );
1028     xy3.SetLinearForm( xMin, xXYZ, yMax, yXYZ );
1029     isInQuad = tface.IsUVInQuad( theUV, xy0,xy1,xy2,xy3 );
1030     nbGetUV += 4;
1031     if ( !isInQuad )
1032     {
1033       dx *= 1.2;
1034       dy *= 1.2;
1035       xSol = 0.5 * (xMax + xMin) ;
1036       ySol = 0.5 * (yMax + yMin) ;
1037       if ( xMin == 0. && yMin == 0. && xMax == 1. && yMax == 1. ) // avoid infinit loop
1038       {
1039 #ifdef _DEBUG_REFINE_
1040         cout << "SMESH_Block::refineParametersOnFace(): tface.IsUVInQuad() fails" << endl;
1041       cout << " nbGetUV = " << nbGetUV << endl;
1042 #endif
1043         break;
1044       }
1045     }
1046   }
1047
1048   // refine solution using half-division technic
1049
1050   gp_XYZ sol = theParams;
1051
1052   const double paramTol = 0.001;
1053   while ( dx > paramTol || dy > paramTol )
1054   {
1055     // divide along X
1056     bool xDivided = ( dx > paramTol );
1057     if ( xDivided )
1058     {
1059       double xMid = 0.5 * ( xMin + xMax );
1060       gp_XYZ parMid1 = xMid * xXYZ + yMin * yXYZ;
1061       gp_XYZ parMid2 = xMid * xXYZ + yMax * yXYZ;
1062       nbGetUV += 4;
1063       if ( tface.IsUVInQuad( theUV, xy0,parMid1,parMid2,xy3 ))
1064       {
1065         xMax = xMid;
1066         xy1 = parMid1; xy2 = parMid2;
1067       }
1068       else if ( tface.IsUVInQuad( theUV, parMid1,xy1,xy2,parMid2 ))
1069       {
1070         nbGetUV += 4;
1071         xMin = xMid;
1072         xy0 = parMid1; xy3 = parMid2;
1073       }
1074       else
1075       {
1076         nbGetUV += 8;
1077         xDivided = false;
1078       }
1079       dx = xMax - xMin;
1080     }
1081     // divide along Y
1082     bool yDivided = ( dy > paramTol );
1083     if ( yDivided )
1084     {
1085       double yMid = 0.5 * ( yMin + yMax );
1086       gp_XYZ parMid2 = xMax * xXYZ + yMid * yXYZ;
1087       gp_XYZ parMid3 = xMin * xXYZ + yMid * yXYZ;
1088       nbGetUV += 4;
1089       if ( tface.IsUVInQuad( theUV, xy0,xy1,parMid2,parMid3 ))
1090       {
1091         yMax = yMid;
1092         xy2 = parMid2; xy3 = parMid3;
1093       }
1094       else if ( tface.IsUVInQuad( theUV, parMid3,parMid2,xy2,xy3 ))
1095       {
1096         nbGetUV += 4;
1097         yMin = yMid;
1098         xy0 = parMid3; xy1 = parMid2;
1099       }
1100       else
1101       {
1102         nbGetUV += 8;
1103         yDivided = false;
1104       }
1105       dy = yMax - yMin;
1106     }
1107     if ( !xDivided && !yDivided )
1108     {
1109 #ifdef _DEBUG_REFINE_
1110       cout << "SMESH_Block::refineParametersOnFace(): nothing divided" << endl;
1111       cout << " nbGetUV = " << nbGetUV << endl;
1112 #endif
1113       break;
1114     }
1115
1116     // evaluate reached distance to thePoint
1117     sol.SetCoord( tface.GetUInd(), 0.5 * ( xMin + xMax ));
1118     sol.SetCoord( tface.GetVInd(), 0.5 * ( yMin + yMax ));
1119     if ( saveBetterSolution( sol, theParams, thePoint.SquareDistance( tface.Point( sol ))))
1120     {
1121 #ifdef _DEBUG_REFINE_
1122       cout << "SMESH_Block::refineParametersOnFace(): dividing suceeded" << endl;
1123       cout << " nbGetUV = " << nbGetUV << endl;
1124 #endif
1125         return true;
1126     }
1127   }
1128 #ifdef _DEBUG_REFINE_
1129   cout << "SMESH_Block::refineParametersOnFace(): dividing Ends at dist " << distance()<< endl;
1130   cout << " nbGetUV = " << nbGetUV << endl;
1131 #endif
1132
1133   return false;
1134 }
1135
1136 //================================================================================
1137 /*!
1138  * \brief Finds parameters corresponding to a given UV of a given face by searching
1139  * around the starting solution
1140  *  \param [in] theUV - the UV to locate
1141  *  \param [in] tface - the face
1142  *  \param [in,out] theParams - the starting parameters to improve
1143  *  \param [in] nbGetWorstLimit - nb of steps from the starting solution w/o improvement
1144  *         to stop searching in this direction
1145  *  \return bool - \c true if found solution is within myTolerance
1146  */
1147 //================================================================================
1148
1149 bool SMESH_Block::findUVAround( const gp_Pnt&             thePoint,
1150                                 const gp_XY&              theUV,
1151                                 const SMESH_Block::TFace& tface,
1152                                 gp_XYZ&                   theParams,
1153                                 int                       nbGetWorstLimit )
1154 {
1155 #ifdef _DEBUG_REFINE_
1156   cout << "SMESH_Block::refineParametersOnFace(): walk around Starts at dist " << distance()<< endl;
1157   cout << " nbGetUV = " << (nbGetUV=0) << endl;
1158 #endif
1159   const double paramTol = 0.001;
1160   const double dx = 0.01, dy = 0.01;
1161   double xMin = theParams.Coord( tface.GetUInd() ), xMax;
1162   double yMin = theParams.Coord( tface.GetVInd() ), yMax;
1163   yMax = yMin;
1164   if ( xMin + dx < 1. )
1165     xMax = xMin + dx;
1166   else
1167     xMax = 1, xMin = 1 - dx;
1168   gp_XYZ sol = theParams;
1169   sol.SetCoord( tface.GetUInd(), xMax );
1170   sol.SetCoord( tface.GetVInd(), yMax );
1171   //nbGetUV++;
1172   if ( saveBetterSolution( sol, theParams, thePoint.SquareDistance( tface.Point( sol ))))
1173     return true;
1174
1175   int xMaxNbGetWorst = 0, xMinNbGetWorst = 0, yMaxNbGetWorst = 0, yMinNbGetWorst = 0;
1176   double xMaxBestDist = 1e100, xMinBestDist = 1e100, yMaxBestDist = 1e100, yMinBestDist = 1e100;
1177   double x, y, bestDist, dist;
1178   while ( xMax - xMin < 1 || yMax - yMin < 1 )
1179   {
1180     // walk along X
1181     if ( yMin > 0. )
1182     {
1183       bestDist = 1e100;
1184       for ( x = Max(0.,xMin); x <= xMax+paramTol; x += dx )
1185       {
1186         y = Max( 0., yMin - dy );
1187         sol.SetCoord( tface.GetUInd(), x );
1188         sol.SetCoord( tface.GetVInd(), y );
1189         //nbGetUV++;
1190         dist = thePoint.SquareDistance( tface.Point( sol ));
1191         bestDist = Min( dist, bestDist );
1192         if ( saveBetterSolution( sol, theParams, dist ))
1193           return true;
1194         sol.SetCoord( tface.GetUInd(), Min( 1., x + 0.5*dx ));
1195         sol.SetCoord( tface.GetVInd(),          y + 0.5*dy );
1196         //nbGetUV++;
1197         dist = thePoint.SquareDistance( tface.Point( sol ));
1198         bestDist = Min( dist, bestDist );
1199         if ( saveBetterSolution( sol, theParams, dist ))
1200           return true;
1201       }
1202       yMin = Max(0., yMin-dy );
1203       yMinNbGetWorst += ( yMinBestDist < bestDist );
1204       yMinBestDist = Min( yMinBestDist, bestDist );
1205       if ( yMinNbGetWorst > nbGetWorstLimit )
1206         yMin = 0;
1207     }
1208     if ( yMax < 1. )
1209     {
1210       bestDist = 1e100;
1211       for ( x = Max(0.,xMin); x <= xMax+paramTol; x += dx )
1212       {
1213         y = Min( 1., yMax + dy );
1214         sol.SetCoord( tface.GetUInd(), x );
1215         sol.SetCoord( tface.GetVInd(), y );
1216         //nbGetUV++;
1217         dist = thePoint.SquareDistance( tface.Point( sol ));
1218         bestDist = Min( dist, bestDist );
1219         if ( saveBetterSolution( sol, theParams, dist ))
1220           return true;
1221         sol.SetCoord( tface.GetUInd(), Min( 1., x + 0.5*dx ));
1222         sol.SetCoord( tface.GetVInd(),          y - 0.5*dy );
1223         //nbGetUV++;
1224         dist = thePoint.SquareDistance( tface.Point( sol ));
1225         bestDist = Min( dist, bestDist );
1226         if ( saveBetterSolution( sol, theParams, dist ))
1227           return true;
1228       }
1229       yMax = Min(1., yMax+dy );
1230       yMaxNbGetWorst += ( yMaxBestDist < bestDist );
1231       yMaxBestDist = Min( yMaxBestDist, bestDist );
1232       if ( yMaxNbGetWorst > nbGetWorstLimit )
1233         yMax = 1;
1234     }
1235     // walk along Y
1236     if ( xMin > 0. )
1237     {
1238       bestDist = 1e100;
1239       for ( y = Max(0.,yMin); y <= yMax+paramTol; y += dy )
1240       {
1241         x = Max( 0., xMin - dx );
1242         sol.SetCoord( tface.GetUInd(), x );
1243         sol.SetCoord( tface.GetVInd(), y );
1244         //nbGetUV++;
1245         dist = thePoint.SquareDistance( tface.Point( sol ));
1246         bestDist = Min( dist, bestDist );
1247         if ( saveBetterSolution( sol, theParams, dist ))
1248           return true;
1249         sol.SetCoord( tface.GetUInd(),          x + 0.5*dx );
1250         sol.SetCoord( tface.GetVInd(), Min( 1., y + 0.5*dy ));
1251         //nbGetUV++;
1252         dist = thePoint.SquareDistance( tface.Point( sol ));
1253         bestDist = Min( dist, bestDist );
1254         if ( saveBetterSolution( sol, theParams, dist ))
1255           return true;
1256       }
1257       xMin = Max(0., xMin-dx );
1258       xMinNbGetWorst += ( xMinBestDist < bestDist );
1259       xMinBestDist = Min( xMinBestDist, bestDist );
1260       if ( xMinNbGetWorst > nbGetWorstLimit )
1261         xMin = 0;
1262     }
1263     if ( xMax < 1. )
1264     {
1265       bestDist = 1e100;
1266       for ( y = Max(0.,yMin); y <= yMax+paramTol; y += dy )
1267       {
1268         x = Min( 1., xMax + dx );
1269         sol.SetCoord( tface.GetUInd(), x );
1270         sol.SetCoord( tface.GetVInd(), y );
1271         //nbGetUV++;
1272         dist = thePoint.SquareDistance( tface.Point( sol ));
1273         bestDist = Min( dist, bestDist );
1274         if ( saveBetterSolution( sol, theParams, dist ))
1275           return true;
1276         sol.SetCoord( tface.GetUInd(),          x - 0.5*dx);
1277         sol.SetCoord( tface.GetVInd(), Min( 1., y + 0.5*dy ));
1278         //nbGetUV++;
1279         dist = thePoint.SquareDistance( tface.Point( sol ));
1280         bestDist = Min( dist, bestDist );
1281         if ( saveBetterSolution( sol, theParams, dist ))
1282           return true;
1283       }
1284       xMax = Min(1., xMax+dx );
1285       xMaxNbGetWorst += ( xMaxBestDist < bestDist );
1286       xMaxBestDist = Min( xMaxBestDist, bestDist );
1287       if ( xMaxNbGetWorst > nbGetWorstLimit )
1288         xMax = 1;
1289     }
1290   }
1291 #ifdef _DEBUG_REFINE_
1292   cout << "SMESH_Block::refineParametersOnFace(): walk around failed at dist " << distance()<< endl;
1293   //cout << " nbGetUV = " << nbGetUV << endl;
1294 #endif
1295
1296   return false;
1297 }
1298
1299 //================================================================================
1300 /*!
1301  * \brief Store a solution if it's better than a previous one
1302  *  \param [in] theNewParams - a new solution
1303  *  \param [out] theParams - the parameters to store solution in
1304  *  \param [in] sqDistance - a square distance reached at theNewParams
1305  *  \return bool - true if the reached distance is within the tolerance
1306  */
1307 //================================================================================
1308
1309 bool SMESH_Block::saveBetterSolution( const gp_XYZ& theNewParams,
1310                                       gp_XYZ&       theParams,
1311                                       double        sqDistance )
1312 {
1313   if ( myValues[ SQUARE_DIST ] > sqDistance )
1314   {
1315     myValues[ SQUARE_DIST ] = sqDistance;
1316     theParams = theNewParams;
1317     if ( distance() <= myTolerance )
1318       return true;
1319   }
1320   return false;
1321 }
1322
1323 //=======================================================================
1324 //function : SetTolerance
1325 //purpose  : set tolerance for ComputeParameters()
1326 //=======================================================================
1327
1328 void SMESH_Block::SetTolerance(const double tol)
1329 {
1330   if ( tol > 0 )
1331     myTolerance = tol;
1332 }
1333
1334 //=======================================================================
1335 //function : IsToleranceReached
1336 //purpose  : return true if solution found by ComputeParameters() is within the tolerance
1337 //=======================================================================
1338
1339 bool SMESH_Block::IsToleranceReached() const
1340 {
1341   return distance() < myTolerance;
1342 }
1343
1344 //=======================================================================
1345 //function : VertexParameters
1346 //purpose  : return parameters of a vertex given by TShapeID
1347 //=======================================================================
1348
1349 bool SMESH_Block::VertexParameters(const int theVertexID, gp_XYZ& theParams)
1350 {
1351   switch ( theVertexID ) {
1352   case ID_V000: theParams.SetCoord(0., 0., 0.); return true;
1353   case ID_V100: theParams.SetCoord(1., 0., 0.); return true;
1354   case ID_V110: theParams.SetCoord(1., 1., 0.); return true;
1355   case ID_V010: theParams.SetCoord(0., 1., 0.); return true;
1356   default:;
1357   }
1358   return false;
1359 }
1360
1361 //=======================================================================
1362 //function : EdgeParameters
1363 //purpose  : return parameters of a point given by theU on edge
1364 //=======================================================================
1365
1366 bool SMESH_Block::EdgeParameters(const int theEdgeID, const double theU, gp_XYZ& theParams)
1367 {
1368   if ( IsEdgeID( theEdgeID )) {
1369     vector< int > vertexVec;
1370     GetEdgeVertexIDs( theEdgeID, vertexVec );
1371     VertexParameters( vertexVec[0], theParams );
1372     TEdge& e = myEdge[ theEdgeID - ID_Ex00 ];
1373     double param = ( theU - e.EndParam(0) ) / ( e.EndParam(1) - e.EndParam(0) );
1374     theParams.SetCoord( e.CoordInd(), param );
1375     return true;
1376   }
1377   return false;
1378 }
1379
1380 //=======================================================================
1381 //function : DumpShapeID
1382 //purpose  : debug an id of a block sub-shape
1383 //=======================================================================
1384
1385 #define CASEDUMP(id,strm) case id: strm << #id; break;
1386
1387 ostream& SMESH_Block::DumpShapeID (const int id, ostream& stream)
1388 {
1389   switch ( id ) {
1390   CASEDUMP( ID_V000, stream );
1391   CASEDUMP( ID_V100, stream );
1392   CASEDUMP( ID_V010, stream );
1393   CASEDUMP( ID_V110, stream );
1394   CASEDUMP( ID_V001, stream );
1395   CASEDUMP( ID_V101, stream );
1396   CASEDUMP( ID_V011, stream );
1397   CASEDUMP( ID_V111, stream );
1398   CASEDUMP( ID_Ex00, stream );
1399   CASEDUMP( ID_Ex10, stream );
1400   CASEDUMP( ID_Ex01, stream );
1401   CASEDUMP( ID_Ex11, stream );
1402   CASEDUMP( ID_E0y0, stream );
1403   CASEDUMP( ID_E1y0, stream );
1404   CASEDUMP( ID_E0y1, stream );
1405   CASEDUMP( ID_E1y1, stream );
1406   CASEDUMP( ID_E00z, stream );
1407   CASEDUMP( ID_E10z, stream );
1408   CASEDUMP( ID_E01z, stream );
1409   CASEDUMP( ID_E11z, stream );
1410   CASEDUMP( ID_Fxy0, stream );
1411   CASEDUMP( ID_Fxy1, stream );
1412   CASEDUMP( ID_Fx0z, stream );
1413   CASEDUMP( ID_Fx1z, stream );
1414   CASEDUMP( ID_F0yz, stream );
1415   CASEDUMP( ID_F1yz, stream );
1416   CASEDUMP( ID_Shell, stream );
1417   default: stream << "ID_INVALID";
1418   }
1419   return stream;
1420 }
1421
1422 //=======================================================================
1423 //function : GetShapeIDByParams
1424 //purpose  : define an id of the block sub-shape by normlized point coord
1425 //=======================================================================
1426
1427 int SMESH_Block::GetShapeIDByParams ( const gp_XYZ& theCoord )
1428 {
1429   //   id ( 0 - 26 ) computation:
1430
1431   //   vertex     ( 0 - 7 )  : id = 1*x + 2*y + 4*z
1432
1433   //   edge || X  ( 8 - 11 ) : id = 8   + 1*y + 2*z
1434   //   edge || Y  ( 12 - 15 ): id = 1*x + 12  + 2*z
1435   //   edge || Z  ( 16 - 19 ): id = 1*x + 2*y + 16
1436
1437   //   face || XY ( 20 - 21 ): id = 8   + 12  + 1*z - 0
1438   //   face || XZ ( 22 - 23 ): id = 8   + 1*y + 16  - 2
1439   //   face || YZ ( 24 - 25 ): id = 1*x + 12  + 16  - 4
1440
1441   static int iAddBnd[]    = { 1, 2, 4 };
1442   static int iAddNotBnd[] = { 8, 12, 16 };
1443   static int iFaceSubst[] = { 0, 2, 4 };
1444
1445   int id = 0;
1446   int iOnBoundary = 0;
1447   for ( int iCoord = 0; iCoord < 3; iCoord++ )
1448   {
1449     double val = theCoord.Coord( iCoord + 1 );
1450     if ( val == 0.0 )
1451       iOnBoundary++;
1452     else if ( val == 1.0 )
1453       id += iAddBnd[ iOnBoundary++ ];
1454     else
1455       id += iAddNotBnd[ iCoord ];
1456   }
1457   if ( iOnBoundary == 1 ) // face
1458     id -= iFaceSubst[ (id - 20) / 4 ];
1459   else if ( iOnBoundary == 0 ) // shell
1460     id = 26;
1461
1462   if ( id > 26 || id < 0 ) {
1463     MESSAGE( "GetShapeIDByParams() = " << id
1464             <<" "<< theCoord.X() <<" "<< theCoord.Y() <<" "<< theCoord.Z() );
1465   }
1466
1467   return id + 1; // shape ids start at 1
1468 }
1469
1470 //================================================================================
1471 /*!
1472  * \brief Return number of wires and a list of oredered edges.
1473  *  \param theFace - the face to process
1474  *  \param theEdges - all ordered edges of theFace (outer edges goes first).
1475  *  \param theNbEdgesInWires - nb of edges (== nb of vertices in closed wire) in each wire
1476  *  \param theFirstVertex - the vertex of the outer wire to set first in the returned
1477  *         list ( theFirstVertex may be NULL )
1478  *  \param theShapeAnalysisAlgo - if true, ShapeAnalysis::OuterWire() is used to find
1479  *         the outer wire else BRepTools::OuterWire() is used.
1480  *  \retval int - nb of wires
1481  *
1482  * Always try to set a seam edge first.
1483  * BRepTools::OuterWire() fails e.g. in the case of issue 0020184,
1484  * ShapeAnalysis::OuterWire() fails in the case of issue 0020452
1485  */
1486 //================================================================================
1487
1488 int SMESH_Block::GetOrderedEdges (const TopoDS_Face&   theFace,
1489                                   list< TopoDS_Edge >& theEdges,
1490                                   list< int >  &       theNbEdgesInWires,
1491                                   TopoDS_Vertex        theFirstVertex,
1492                                   const bool           theShapeAnalysisAlgo)
1493 {
1494   // put wires in a list, so that an outer wire comes first
1495   list<TopoDS_Wire> aWireList;
1496   TopoDS_Wire anOuterWire =
1497     theShapeAnalysisAlgo ? ShapeAnalysis::OuterWire( theFace ) : BRepTools::OuterWire( theFace );
1498   for ( TopoDS_Iterator wIt (theFace); wIt.More(); wIt.Next() )
1499     if ( wIt.Value().ShapeType() == TopAbs_WIRE ) // it can be internal vertex!
1500     {
1501       if ( !anOuterWire.IsSame( wIt.Value() ))
1502         aWireList.push_back( TopoDS::Wire( wIt.Value() ));
1503       else
1504         aWireList.push_front( TopoDS::Wire( wIt.Value() ));
1505     }
1506
1507   // loop on edges of wires
1508   theNbEdgesInWires.clear();
1509   list<TopoDS_Wire>::iterator wlIt = aWireList.begin();
1510   for ( ; wlIt != aWireList.end(); wlIt++ )
1511   {
1512     int iE;
1513     BRepTools_WireExplorer wExp( *wlIt, theFace );
1514     for ( iE = 0; wExp.More(); wExp.Next(), iE++ )
1515     {
1516       TopoDS_Edge edge = wExp.Current();
1517       // commented for issue 0020557, other related ones: 0020526, PAL19080
1518       // edge = TopoDS::Edge( edge.Oriented( wExp.Orientation() ));
1519       theEdges.push_back( edge );
1520     }
1521     if ( iE == 0 ) // wExp returns nothing if e.g. the wire contains one internal edge
1522     { // Issue 0020676
1523       for ( TopoDS_Iterator e( *wlIt ); e.More(); e.Next(), ++iE )
1524         theEdges.push_back( TopoDS::Edge( e.Value() ));
1525     }
1526     theNbEdgesInWires.push_back( iE );
1527     iE = 0;
1528     if ( wlIt == aWireList.begin() && theEdges.size() > 1 ) { // the outer wire
1529       // orient closed edges
1530       list< TopoDS_Edge >::iterator eIt, eIt2;
1531       for ( eIt = theEdges.begin(); eIt != theEdges.end(); eIt++ )
1532       {
1533         TopoDS_Edge& edge = *eIt;
1534         if ( TopExp::FirstVertex( edge ).IsSame( TopExp::LastVertex( edge ) ))
1535         {
1536           eIt2 = eIt;
1537           bool isNext = ( eIt2 == theEdges.begin() );
1538           TopoDS_Edge edge2 = isNext ? *(++eIt2) : *(--eIt2);
1539           double f1,l1,f2,l2;
1540           Handle(Geom2d_Curve) c1 = BRep_Tool::CurveOnSurface( edge, theFace, f1,l1 );
1541           Handle(Geom2d_Curve) c2 = BRep_Tool::CurveOnSurface( edge2, theFace, f2,l2 );
1542           gp_Pnt2d pf = c1->Value( edge.Orientation() == TopAbs_FORWARD ? f1 : l1 );
1543           gp_Pnt2d pl = c1->Value( edge.Orientation() == TopAbs_FORWARD ? l1 : f1 );
1544           bool isFirst = ( edge2.Orientation() == TopAbs_FORWARD ? isNext : !isNext );
1545           gp_Pnt2d p2 = c2->Value( isFirst ? f2 : l2 );
1546           isFirst = ( p2.SquareDistance( pf ) < p2.SquareDistance( pl ));
1547           if ( isNext ? isFirst : !isFirst )
1548             edge.Reverse();
1549           // to make a seam go first
1550           if ( theFirstVertex.IsNull() )
1551             theFirstVertex = TopExp::FirstVertex( edge, true );
1552         }
1553       }
1554       // rotate theEdges until it begins from theFirstVertex
1555       if ( ! theFirstVertex.IsNull() ) {
1556         TopoDS_Vertex vv[2];
1557         TopExp::Vertices( theEdges.front(), vv[0], vv[1], true );
1558         // on closed face, make seam edge the first in the list
1559         while ( !vv[0].IsSame( theFirstVertex ) || vv[0].IsSame( vv[1] ))
1560         {
1561           theEdges.splice(theEdges.end(), theEdges,
1562                           theEdges.begin(), ++theEdges.begin());
1563           TopExp::Vertices( theEdges.front(), vv[0], vv[1], true );
1564           if ( iE++ > theNbEdgesInWires.back() ) {
1565 #ifdef _DEBUG_
1566             gp_Pnt p = BRep_Tool::Pnt( theFirstVertex );
1567             MESSAGE ( " : Warning : vertex "<< theFirstVertex.TShape().operator->()
1568                    << " ( " << p.X() << " " << p.Y() << " " << p.Z() << " )"
1569                    << " not found in outer wire of face "<< theFace.TShape().operator->()
1570                    << " with vertices: " );
1571             wExp.Init( *wlIt, theFace );
1572             for ( int i = 0; wExp.More(); wExp.Next(), i++ )
1573             {
1574               TopoDS_Edge edge = wExp.Current();
1575               edge = TopoDS::Edge( edge.Oriented( wExp.Orientation() ));
1576               TopoDS_Vertex v = TopExp::FirstVertex( edge, true );
1577               gp_Pnt p = BRep_Tool::Pnt( v );
1578               MESSAGE_ADD ( i << " " << v.TShape().operator->() << " "
1579                             << p.X() << " " << p.Y() << " " << p.Z() << " " << std::endl );
1580             }
1581 #endif
1582             break; // break infinite loop
1583           }
1584         }
1585       }
1586     } // end outer wire
1587   }
1588
1589   return aWireList.size();
1590 }
1591 //================================================================================
1592 /*!
1593  * \brief Call it after geometry initialisation
1594  */
1595 //================================================================================
1596
1597 void SMESH_Block::init()
1598 {
1599   myNbIterations = 0;
1600   mySumDist = 0;
1601   myGridComputed = false;
1602 }
1603
1604 //=======================================================================
1605 //function : LoadMeshBlock
1606 //purpose  : prepare to work with theVolume
1607 //=======================================================================
1608
1609 #define gpXYZ(n) gp_XYZ(n->X(),n->Y(),n->Z())
1610
1611 bool SMESH_Block::LoadMeshBlock(const SMDS_MeshVolume*        theVolume,
1612                                 const int                     theNode000Index,
1613                                 const int                     theNode001Index,
1614                                 vector<const SMDS_MeshNode*>& theOrderedNodes)
1615 {
1616   MESSAGE(" ::LoadMeshBlock()");
1617   init();
1618
1619   SMDS_VolumeTool vTool;
1620   if (!vTool.Set( theVolume ) || vTool.NbNodes() != 8 ||
1621       !vTool.IsLinked( theNode000Index, theNode001Index )) {
1622     MESSAGE(" Bad arguments ");
1623     return false;
1624   }
1625   vTool.SetExternalNormal();
1626   // In terms of indices used for access to nodes and faces in SMDS_VolumeTool:
1627   int V000, V100, V010, V110, V001, V101, V011, V111; // 8 vertices
1628   int Fxy0, Fxy1; // bottom and top faces
1629   // vertices of faces
1630   vector<int> vFxy0, vFxy1;
1631
1632   V000 = theNode000Index;
1633   V001 = theNode001Index;
1634
1635   // get faces sharing V000 and V001
1636   list<int> fV000, fV001;
1637   int i, iF, iE, iN;
1638   for ( iF = 0; iF < vTool.NbFaces(); ++iF ) {
1639     const int* nid = vTool.GetFaceNodesIndices( iF );
1640     for ( iN = 0; iN < 4; ++iN )
1641       if ( nid[ iN ] == V000 ) {
1642         fV000.push_back( iF );
1643       } else if ( nid[ iN ] == V001 ) {
1644         fV001.push_back( iF );
1645       }
1646   }
1647
1648   // find the bottom (Fxy0), the top (Fxy1) faces
1649   list<int>::iterator fIt1, fIt2, Fxy0Pos;
1650   for ( fIt1 = fV000.begin(); fIt1 != fV000.end(); fIt1++) {
1651     fIt2 = std::find( fV001.begin(), fV001.end(), *fIt1 );
1652     if ( fIt2 != fV001.end() ) { // *fIt1 is in the both lists
1653       fV001.erase( fIt2 ); // erase Fx0z or F0yz from fV001
1654     } else { // *fIt1 is in fV000 only
1655       Fxy0Pos = fIt1; // points to Fxy0
1656     }
1657   }
1658   Fxy0 = *Fxy0Pos;
1659   Fxy1 = fV001.front();
1660   const SMDS_MeshNode** nn = vTool.GetNodes();
1661
1662   // find bottom veritices, their order is that a face normal is external
1663   vFxy0.resize(4);
1664   const int* nid = vTool.GetFaceNodesIndices( Fxy0 );
1665   for ( i = 0; i < 4; ++i )
1666     if ( nid[ i ] == V000 )
1667       break;
1668   for ( iN = 0; iN < 4; ++iN, ++i ) {
1669     if ( i == 4 ) i = 0;
1670     vFxy0[ iN ] = nid[ i ];
1671   }
1672   // find top veritices, their order is that a face normal is external
1673   vFxy1.resize(4);
1674   nid = vTool.GetFaceNodesIndices( Fxy1 );
1675   for ( i = 0; i < 4; ++i )
1676     if ( nid[ i ] == V001 )
1677       break;
1678   for ( iN = 0; iN < 4; ++iN, ++i ) {
1679     if ( i == 4 ) i = 0;
1680     vFxy1[ iN ] = nid[ i ];
1681   }
1682   // find indices of the rest veritices
1683   V100 = vFxy0[3];
1684   V010 = vFxy0[1];
1685   V110 = vFxy0[2];
1686   V101 = vFxy1[1];
1687   V011 = vFxy1[3];
1688   V111 = vFxy1[2];
1689
1690   // set points coordinates
1691   myPnt[ ID_V000 - 1 ] = gpXYZ( nn[ V000 ] );
1692   myPnt[ ID_V100 - 1 ] = gpXYZ( nn[ V100 ] );
1693   myPnt[ ID_V010 - 1 ] = gpXYZ( nn[ V010 ] );
1694   myPnt[ ID_V110 - 1 ] = gpXYZ( nn[ V110 ] );
1695   myPnt[ ID_V001 - 1 ] = gpXYZ( nn[ V001 ] );
1696   myPnt[ ID_V101 - 1 ] = gpXYZ( nn[ V101 ] );
1697   myPnt[ ID_V011 - 1 ] = gpXYZ( nn[ V011 ] );
1698   myPnt[ ID_V111 - 1 ] = gpXYZ( nn[ V111 ] );
1699
1700   // fill theOrderedNodes
1701   theOrderedNodes.resize( 8 );
1702   theOrderedNodes[ 0 ] = nn[ V000 ];
1703   theOrderedNodes[ 1 ] = nn[ V100 ];
1704   theOrderedNodes[ 2 ] = nn[ V010 ];
1705   theOrderedNodes[ 3 ] = nn[ V110 ];
1706   theOrderedNodes[ 4 ] = nn[ V001 ];
1707   theOrderedNodes[ 5 ] = nn[ V101 ];
1708   theOrderedNodes[ 6 ] = nn[ V011 ];
1709   theOrderedNodes[ 7 ] = nn[ V111 ];
1710
1711   // fill edges
1712   vector< int > vertexVec;
1713   for ( iE = 0; iE < NbEdges(); ++iE ) {
1714     GetEdgeVertexIDs(( iE + ID_FirstE ), vertexVec );
1715     myEdge[ iE ].Set(( iE + ID_FirstE ),
1716                      myPnt[ vertexVec[0] - 1 ],
1717                      myPnt[ vertexVec[1] - 1 ]);
1718   }
1719
1720   // fill faces' corners
1721   for ( iF = ID_Fxy0; iF < ID_Shell; ++iF )
1722   {
1723     TFace& tFace = myFace[ iF - ID_FirstF ];
1724     vector< int > edgeIdVec(4, -1);
1725     GetFaceEdgesIDs( iF, edgeIdVec );
1726     tFace.Set( iF, myEdge[ edgeIdVec [ 0 ] - ID_Ex00], myEdge[ edgeIdVec [ 1 ] - ID_Ex00]);
1727   }
1728
1729   return true;
1730 }
1731
1732 //=======================================================================
1733 //function : LoadBlockShapes
1734 //purpose  : Initialize block geometry with theShell,
1735 //           add sub-shapes of theBlock to theShapeIDMap so that they get
1736 //           IDs acoording to enum TShapeID
1737 //=======================================================================
1738
1739 bool SMESH_Block::LoadBlockShapes(const TopoDS_Shell&         theShell,
1740                                   const TopoDS_Vertex&        theVertex000,
1741                                   const TopoDS_Vertex&        theVertex001,
1742                                   TopTools_IndexedMapOfOrientedShape& theShapeIDMap )
1743 {
1744   MESSAGE(" ::LoadBlockShapes()");
1745   return ( FindBlockShapes( theShell, theVertex000, theVertex001, theShapeIDMap ) &&
1746            LoadBlockShapes( theShapeIDMap ));
1747 }
1748
1749 //=======================================================================
1750 //function : LoadBlockShapes
1751 //purpose  : add sub-shapes of theBlock to theShapeIDMap so that they get
1752 //           IDs acoording to enum TShapeID
1753 //=======================================================================
1754
1755 bool SMESH_Block::FindBlockShapes(const TopoDS_Shell&         theShell,
1756                                   const TopoDS_Vertex&        theVertex000,
1757                                   const TopoDS_Vertex&        theVertex001,
1758                                   TopTools_IndexedMapOfOrientedShape& theShapeIDMap )
1759 {
1760   MESSAGE(" ::FindBlockShapes()");
1761
1762   // 8 vertices
1763   TopoDS_Shape V000, V100, V010, V110, V001, V101, V011, V111;
1764   // 12 edges
1765   TopoDS_Shape Ex00, Ex10, Ex01, Ex11;
1766   TopoDS_Shape E0y0, E1y0, E0y1, E1y1;
1767   TopoDS_Shape E00z, E10z, E01z, E11z;
1768   // 6 faces
1769   TopoDS_Shape Fxy0, Fx0z, F0yz, Fxy1, Fx1z, F1yz;
1770
1771   // nb of faces bound to a vertex in TopTools_IndexedDataMapOfShapeListOfShape
1772   // filled by TopExp::MapShapesAndAncestors()
1773   const int NB_FACES_BY_VERTEX = 6;
1774
1775   TopTools_IndexedDataMapOfShapeListOfShape vfMap;
1776   TopExp::MapShapesAndAncestors( theShell, TopAbs_VERTEX, TopAbs_FACE, vfMap );
1777   if ( vfMap.Extent() != 8 ) {
1778     MESSAGE(" Wrong nb of vertices in the block: " << vfMap.Extent() );
1779     return false;
1780   }
1781
1782   V000 = theVertex000;
1783   V001 = theVertex001;
1784
1785   if ( V000.IsNull() ) {
1786     // find vertex 000 - the one with smallest coordinates
1787     double minVal = DBL_MAX, minX, val;
1788     for ( int i = 1; i <= 8; i++ ) {
1789       const TopoDS_Vertex& v = TopoDS::Vertex( vfMap.FindKey( i ));
1790       gp_Pnt P = BRep_Tool::Pnt( v );
1791       val = P.X() + P.Y() + P.Z();
1792       if ( val < minVal || ( val == minVal && P.X() < minX )) {
1793         V000 = v;
1794         minVal = val;
1795         minX = P.X();
1796       }
1797     }
1798     // find vertex 001 - the one on the most vertical edge passing through V000
1799     TopTools_IndexedDataMapOfShapeListOfShape veMap;
1800     TopExp::MapShapesAndAncestors( theShell, TopAbs_VERTEX, TopAbs_EDGE, veMap );
1801     gp_Vec dir001 = gp::DZ();
1802     gp_Pnt p000 = BRep_Tool::Pnt( TopoDS::Vertex( V000 ));
1803     double maxVal = -DBL_MAX;
1804     TopTools_ListIteratorOfListOfShape eIt ( veMap.FindFromKey( V000 ));
1805     for (  ; eIt.More(); eIt.Next() ) {
1806       const TopoDS_Edge& e = TopoDS::Edge( eIt.Value() );
1807       TopoDS_Vertex v = TopExp::FirstVertex( e );
1808       if ( v.IsSame( V000 ))
1809         v = TopExp::LastVertex( e );
1810       val = dir001 * gp_Vec( p000, BRep_Tool::Pnt( v )).Normalized();
1811       if ( val > maxVal ) {
1812         V001 = v;
1813         maxVal = val;
1814       }
1815     }
1816   }
1817
1818   // find the bottom (Fxy0), Fx0z and F0yz faces
1819
1820   const TopTools_ListOfShape& f000List = vfMap.FindFromKey( V000 );
1821   const TopTools_ListOfShape& f001List = vfMap.FindFromKey( V001 );
1822   if (f000List.Extent() != NB_FACES_BY_VERTEX ||
1823       f001List.Extent() != NB_FACES_BY_VERTEX ) {
1824     MESSAGE(" LoadBlockShapes() " << f000List.Extent() << " " << f001List.Extent());
1825     return false;
1826   }
1827   TopTools_ListIteratorOfListOfShape f001It, f000It ( f000List );
1828   int i, j, iFound1, iFound2;
1829   for ( j = 0; f000It.More(); f000It.Next(), j++ )
1830   {
1831     if ( NB_FACES_BY_VERTEX == 6 && j % 2 ) continue; // each face encounters twice
1832     const TopoDS_Shape& F = f000It.Value();
1833     for ( i = 0, f001It.Initialize( f001List ); f001It.More(); f001It.Next(), i++ ) {
1834       if ( NB_FACES_BY_VERTEX == 6 && i % 2 ) continue; // each face encounters twice
1835       if ( F.IsSame( f001It.Value() ))
1836         break;
1837     }
1838     if ( f001It.More() ) // Fx0z or F0yz found
1839       if ( Fx0z.IsNull() ) {
1840         Fx0z = F;
1841         iFound1 = i;
1842       } else {
1843         F0yz = F;
1844         iFound2 = i;
1845       }
1846     else // F is the bottom face
1847       Fxy0 = F;
1848   }
1849   if ( Fxy0.IsNull() || Fx0z.IsNull() || F0yz.IsNull() ) {
1850     MESSAGE( Fxy0.IsNull() <<" "<< Fx0z.IsNull() <<" "<< F0yz.IsNull() );
1851     return false;
1852   }
1853
1854   // choose the top face (Fxy1)
1855   for ( i = 0, f001It.Initialize( f001List ); f001It.More(); f001It.Next(), i++ ) {
1856     if ( NB_FACES_BY_VERTEX == 6 && i % 2 ) continue; // each face encounters twice
1857     if ( i != iFound1 && i != iFound2 )
1858       break;
1859   }
1860   Fxy1 = f001It.Value();
1861   if ( Fxy1.IsNull() ) {
1862     MESSAGE(" LoadBlockShapes() error ");
1863     return false;
1864   }
1865
1866   // find bottom edges and veritices
1867   list< TopoDS_Edge > eList;
1868   list< int >         nbVertexInWires;
1869   GetOrderedEdges( TopoDS::Face( Fxy0 ), eList, nbVertexInWires, TopoDS::Vertex( V000 ) );
1870   if ( nbVertexInWires.size() != 1 || nbVertexInWires.front() != 4 ) {
1871     MESSAGE(" LoadBlockShapes() error ");
1872     return false;
1873   }
1874   list< TopoDS_Edge >::iterator elIt = eList.begin();
1875   for ( i = 0; elIt != eList.end(); elIt++, i++ )
1876     switch ( i ) {
1877     case 0: E0y0 = *elIt; V010 = TopExp::LastVertex( *elIt, true ); break;
1878     case 1: Ex10 = *elIt; V110 = TopExp::LastVertex( *elIt, true ); break;
1879     case 2: E1y0 = *elIt; V100 = TopExp::LastVertex( *elIt, true ); break;
1880     case 3: Ex00 = *elIt; break;
1881     default:;
1882     }
1883   if ( i != 4 || E0y0.IsNull() || Ex10.IsNull() || E1y0.IsNull() || Ex00.IsNull() ) {
1884     MESSAGE(" LoadBlockShapes() error, eList.size()=" << eList.size());
1885     return false;
1886   }
1887
1888
1889   // find top edges and veritices
1890   eList.clear();
1891   GetOrderedEdges( TopoDS::Face( Fxy1 ), eList, nbVertexInWires, TopoDS::Vertex( V001 ) );
1892   if ( nbVertexInWires.size() != 1 || nbVertexInWires.front() != 4 ) {
1893     MESSAGE(" LoadBlockShapes() error ");
1894     return false;
1895   }
1896   for ( i = 0, elIt = eList.begin(); elIt != eList.end(); elIt++, i++ )
1897     switch ( i ) {
1898     case 0: Ex01 = *elIt; V101 = TopExp::LastVertex( *elIt, true ); break;
1899     case 1: E1y1 = *elIt; V111 = TopExp::LastVertex( *elIt, true ); break;
1900     case 2: Ex11 = *elIt; V011 = TopExp::LastVertex( *elIt, true ); break;
1901     case 3: E0y1 = *elIt; break;
1902     default:;
1903     }
1904   if ( i != 4 || Ex01.IsNull() || E1y1.IsNull() || Ex11.IsNull() || E0y1.IsNull() ) {
1905     MESSAGE(" LoadBlockShapes() error, eList.size()=" << eList.size());
1906     return false;
1907   }
1908
1909   // swap Fx0z and F0yz if necessary
1910   TopExp_Explorer exp( Fx0z, TopAbs_VERTEX );
1911   for ( ; exp.More(); exp.Next() ) // Fx0z shares V101 and V100
1912     if ( V101.IsSame( exp.Current() ) || V100.IsSame( exp.Current() ))
1913       break; // V101 or V100 found
1914   if ( !exp.More() ) { // not found
1915     std::swap( Fx0z, F0yz);
1916   }
1917
1918   // find Fx1z and F1yz faces
1919   const TopTools_ListOfShape& f111List = vfMap.FindFromKey( V111 );
1920   const TopTools_ListOfShape& f110List = vfMap.FindFromKey( V110 );
1921   if (f111List.Extent() != NB_FACES_BY_VERTEX ||
1922       f110List.Extent() != NB_FACES_BY_VERTEX ) {
1923     MESSAGE(" LoadBlockShapes() " << f111List.Extent() << " " << f110List.Extent());
1924     return false;
1925   }
1926   TopTools_ListIteratorOfListOfShape f111It, f110It ( f110List);
1927   for ( j = 0 ; f110It.More(); f110It.Next(), j++ ) {
1928     if ( NB_FACES_BY_VERTEX == 6 && j % 2 ) continue; // each face encounters twice
1929     const TopoDS_Shape& F = f110It.Value();
1930     for ( i = 0, f111It.Initialize( f111List ); f111It.More(); f111It.Next(), i++ ) {
1931       if ( NB_FACES_BY_VERTEX == 6 && i % 2 ) continue; // each face encounters twice
1932       if ( F.IsSame( f111It.Value() )) { // Fx1z or F1yz found
1933         if ( Fx1z.IsNull() )
1934           Fx1z = F;
1935         else
1936           F1yz = F;
1937       }
1938     }
1939   }
1940   if ( Fx1z.IsNull() || F1yz.IsNull() ) {
1941     MESSAGE(" LoadBlockShapes() error ");
1942     return false;
1943   }
1944
1945   // swap Fx1z and F1yz if necessary
1946   for ( exp.Init( Fx1z, TopAbs_VERTEX ); exp.More(); exp.Next() )
1947     if ( V010.IsSame( exp.Current() ) || V011.IsSame( exp.Current() ))
1948       break;
1949   if ( !exp.More() ) {
1950     std::swap( Fx1z, F1yz);
1951   }
1952
1953   // find vertical edges
1954   for ( exp.Init( Fx0z, TopAbs_EDGE ); exp.More(); exp.Next() ) {
1955     const TopoDS_Edge& edge = TopoDS::Edge( exp.Current() );
1956     const TopoDS_Shape& vFirst = TopExp::FirstVertex( edge, true );
1957     if ( vFirst.IsSame( V001 ))
1958       E00z = edge;
1959     else if ( vFirst.IsSame( V100 ))
1960       E10z = edge;
1961   }
1962   if ( E00z.IsNull() || E10z.IsNull() ) {
1963     MESSAGE(" LoadBlockShapes() error ");
1964     return false;
1965   }
1966   for ( exp.Init( Fx1z, TopAbs_EDGE ); exp.More(); exp.Next() ) {
1967     const TopoDS_Edge& edge = TopoDS::Edge( exp.Current() );
1968     const TopoDS_Shape& vFirst = TopExp::FirstVertex( edge, true );
1969     if ( vFirst.IsSame( V111 ))
1970       E11z = edge;
1971     else if ( vFirst.IsSame( V010 ))
1972       E01z = edge;
1973   }
1974   if ( E01z.IsNull() || E11z.IsNull() ) {
1975     MESSAGE(" LoadBlockShapes() error ");
1976     return false;
1977   }
1978
1979   // load shapes in theShapeIDMap
1980
1981   theShapeIDMap.Clear();
1982
1983   theShapeIDMap.Add(V000.Oriented( TopAbs_FORWARD ));
1984   theShapeIDMap.Add(V100.Oriented( TopAbs_FORWARD ));
1985   theShapeIDMap.Add(V010.Oriented( TopAbs_FORWARD ));
1986   theShapeIDMap.Add(V110.Oriented( TopAbs_FORWARD ));
1987   theShapeIDMap.Add(V001.Oriented( TopAbs_FORWARD ));
1988   theShapeIDMap.Add(V101.Oriented( TopAbs_FORWARD ));
1989   theShapeIDMap.Add(V011.Oriented( TopAbs_FORWARD ));
1990   theShapeIDMap.Add(V111.Oriented( TopAbs_FORWARD ));
1991
1992   theShapeIDMap.Add(Ex00);
1993   theShapeIDMap.Add(Ex10);
1994   theShapeIDMap.Add(Ex01);
1995   theShapeIDMap.Add(Ex11);
1996
1997   theShapeIDMap.Add(E0y0);
1998   theShapeIDMap.Add(E1y0);
1999   theShapeIDMap.Add(E0y1);
2000   theShapeIDMap.Add(E1y1);
2001
2002   theShapeIDMap.Add(E00z);
2003   theShapeIDMap.Add(E10z);
2004   theShapeIDMap.Add(E01z);
2005   theShapeIDMap.Add(E11z);
2006
2007   theShapeIDMap.Add(Fxy0);
2008   theShapeIDMap.Add(Fxy1);
2009   theShapeIDMap.Add(Fx0z);
2010   theShapeIDMap.Add(Fx1z);
2011   theShapeIDMap.Add(F0yz);
2012   theShapeIDMap.Add(F1yz);
2013
2014   theShapeIDMap.Add(theShell);
2015
2016   return true;
2017 }
2018
2019 //================================================================================
2020 /*!
2021  * \brief Initialize block geometry with shapes from theShapeIDMap
2022   * \param theShapeIDMap - map of block sub-shapes
2023   * \retval bool - is a success
2024  */
2025 //================================================================================
2026
2027 bool SMESH_Block::LoadBlockShapes(const TopTools_IndexedMapOfOrientedShape& theShapeIDMap)
2028 {
2029   init();
2030
2031   // store shapes geometry
2032   for ( int shapeID = 1; shapeID < theShapeIDMap.Extent(); shapeID++ )
2033   {
2034     const TopoDS_Shape& S = theShapeIDMap( shapeID );
2035     switch ( S.ShapeType() )
2036     {
2037     case TopAbs_VERTEX: {
2038
2039       if ( !IsVertexID( ID_V111 )) return false;
2040       myPnt[ shapeID - ID_V000 ] = BRep_Tool::Pnt( TopoDS::Vertex( S )).XYZ();
2041       break;
2042     }
2043     case TopAbs_EDGE: {
2044
2045       if ( !IsEdgeID( shapeID )) return false;
2046       const TopoDS_Edge& edge = TopoDS::Edge( S );
2047       TEdge& tEdge = myEdge[ shapeID - ID_FirstE ];
2048       tEdge.Set( shapeID,
2049                  new BRepAdaptor_Curve( edge ),
2050                  IsForwardEdge( edge, theShapeIDMap ));
2051       break;
2052     }
2053     case TopAbs_FACE: {
2054
2055       if ( !LoadFace( TopoDS::Face( S ), shapeID, theShapeIDMap ))
2056         return false;
2057       break;
2058     }
2059     default: break;
2060     }
2061   } // loop on shapes in theShapeIDMap
2062
2063   return true;
2064 }
2065
2066 //================================================================================
2067 /*!
2068  * \brief Load face geometry
2069   * \param theFace - face
2070   * \param theFaceID - face in-block ID
2071   * \param theShapeIDMap - map of block sub-shapes
2072   * \retval bool - is a success
2073  *
2074  * It is enough to compute params or coordinates on the face.
2075  * Face sub-shapes must be loaded into theShapeIDMap before
2076  */
2077 //================================================================================
2078
2079 bool SMESH_Block::LoadFace(const TopoDS_Face& theFace,
2080                            const int          theFaceID,
2081                            const TopTools_IndexedMapOfOrientedShape& theShapeIDMap)
2082 {
2083   if ( !IsFaceID( theFaceID ) ) return false;
2084   // pcurves
2085   Adaptor2d_Curve2d* c2d[4];
2086   bool isForward[4];
2087   vector< int > edgeIdVec;
2088   GetFaceEdgesIDs( theFaceID, edgeIdVec );
2089   for ( int iE = 0; iE < edgeIdVec.size(); iE++ ) // loop on 4 edges
2090   {
2091     if ( edgeIdVec[ iE ] > theShapeIDMap.Extent() )
2092       return false;
2093     const TopoDS_Edge& edge = TopoDS::Edge( theShapeIDMap( edgeIdVec[ iE ]));
2094     c2d[ iE ] = new BRepAdaptor_Curve2d( edge, theFace );
2095     isForward[ iE ] = IsForwardEdge( edge, theShapeIDMap );
2096   }
2097   TFace& tFace = myFace[ theFaceID - ID_FirstF ];
2098   tFace.Set( theFaceID, new BRepAdaptor_Surface( theFace ), c2d, isForward );
2099   return true;
2100 }
2101
2102 //================================================================================
2103 /*!
2104  * \brief/ Insert theShape into theShapeIDMap with theShapeID
2105   * \param theShape - shape to insert
2106   * \param theShapeID - shape in-block ID
2107   * \param theShapeIDMap - map of block sub-shapes
2108  */
2109 //================================================================================
2110
2111 bool SMESH_Block::Insert(const TopoDS_Shape& theShape,
2112                          const int           theShapeID,
2113                          TopTools_IndexedMapOfOrientedShape& theShapeIDMap)
2114 {
2115   if ( !theShape.IsNull() && theShapeID > 0 )
2116   {
2117     if ( theShapeIDMap.Contains( theShape ))
2118       return ( theShapeIDMap.FindIndex( theShape ) == theShapeID );
2119
2120     if ( theShapeID <= theShapeIDMap.Extent() ) {
2121         theShapeIDMap.Substitute( theShapeID, theShape );
2122     }
2123     else {
2124       while ( theShapeIDMap.Extent() < theShapeID - 1 ) {
2125         TopoDS_Compound comp;
2126         BRep_Builder().MakeCompound( comp );
2127         theShapeIDMap.Add( comp );
2128       }
2129       theShapeIDMap.Add( theShape );
2130     }
2131     return true;
2132   }
2133   return false;
2134 }
2135
2136 //=======================================================================
2137 //function : GetFaceEdgesIDs
2138 //purpose  : return edges IDs in the order u0, u1, 0v, 1v
2139 //           u0 means "|| u, v == 0"
2140 //=======================================================================
2141
2142 void SMESH_Block::GetFaceEdgesIDs (const int faceID, vector< int >& edgeVec )
2143 {
2144   edgeVec.resize( 4 );
2145   switch ( faceID ) {
2146   case ID_Fxy0:
2147     edgeVec[ 0 ] = ID_Ex00;
2148     edgeVec[ 1 ] = ID_Ex10;
2149     edgeVec[ 2 ] = ID_E0y0;
2150     edgeVec[ 3 ] = ID_E1y0;
2151     break;
2152   case ID_Fxy1:
2153     edgeVec[ 0 ] = ID_Ex01;
2154     edgeVec[ 1 ] = ID_Ex11;
2155     edgeVec[ 2 ] = ID_E0y1;
2156     edgeVec[ 3 ] = ID_E1y1;
2157     break;
2158   case ID_Fx0z:
2159     edgeVec[ 0 ] = ID_Ex00;
2160     edgeVec[ 1 ] = ID_Ex01;
2161     edgeVec[ 2 ] = ID_E00z;
2162     edgeVec[ 3 ] = ID_E10z;
2163     break;
2164   case ID_Fx1z:
2165     edgeVec[ 0 ] = ID_Ex10;
2166     edgeVec[ 1 ] = ID_Ex11;
2167     edgeVec[ 2 ] = ID_E01z;
2168     edgeVec[ 3 ] = ID_E11z;
2169     break;
2170   case ID_F0yz:
2171     edgeVec[ 0 ] = ID_E0y0;
2172     edgeVec[ 1 ] = ID_E0y1;
2173     edgeVec[ 2 ] = ID_E00z;
2174     edgeVec[ 3 ] = ID_E01z;
2175     break;
2176   case ID_F1yz:
2177     edgeVec[ 0 ] = ID_E1y0;
2178     edgeVec[ 1 ] = ID_E1y1;
2179     edgeVec[ 2 ] = ID_E10z;
2180     edgeVec[ 3 ] = ID_E11z;
2181     break;
2182   default:
2183     MESSAGE(" GetFaceEdgesIDs(), wrong face ID: " << faceID );
2184   }
2185 }
2186
2187 //=======================================================================
2188 //function : GetEdgeVertexIDs
2189 //purpose  : return vertex IDs of an edge
2190 //=======================================================================
2191
2192 void SMESH_Block::GetEdgeVertexIDs (const int edgeID, vector< int >& vertexVec )
2193 {
2194   vertexVec.resize( 2 );
2195   switch ( edgeID ) {
2196
2197   case ID_Ex00:
2198     vertexVec[ 0 ] = ID_V000;
2199     vertexVec[ 1 ] = ID_V100;
2200     break;
2201   case ID_Ex10:
2202     vertexVec[ 0 ] = ID_V010;
2203     vertexVec[ 1 ] = ID_V110;
2204     break;
2205   case ID_Ex01:
2206     vertexVec[ 0 ] = ID_V001;
2207     vertexVec[ 1 ] = ID_V101;
2208     break;
2209   case ID_Ex11:
2210     vertexVec[ 0 ] = ID_V011;
2211     vertexVec[ 1 ] = ID_V111;
2212     break;
2213
2214   case ID_E0y0:
2215     vertexVec[ 0 ] = ID_V000;
2216     vertexVec[ 1 ] = ID_V010;
2217     break;
2218   case ID_E1y0:
2219     vertexVec[ 0 ] = ID_V100;
2220     vertexVec[ 1 ] = ID_V110;
2221     break;
2222   case ID_E0y1:
2223     vertexVec[ 0 ] = ID_V001;
2224     vertexVec[ 1 ] = ID_V011;
2225     break;
2226   case ID_E1y1:
2227     vertexVec[ 0 ] = ID_V101;
2228     vertexVec[ 1 ] = ID_V111;
2229     break;
2230
2231   case ID_E00z:
2232     vertexVec[ 0 ] = ID_V000;
2233     vertexVec[ 1 ] = ID_V001;
2234     break;
2235   case ID_E10z:
2236     vertexVec[ 0 ] = ID_V100;
2237     vertexVec[ 1 ] = ID_V101;
2238     break;
2239   case ID_E01z:
2240     vertexVec[ 0 ] = ID_V010;
2241     vertexVec[ 1 ] = ID_V011;
2242     break;
2243   case ID_E11z:
2244     vertexVec[ 0 ] = ID_V110;
2245     vertexVec[ 1 ] = ID_V111;
2246     break;
2247   default:
2248     vertexVec.resize(0);
2249     MESSAGE(" GetEdgeVertexIDs(), wrong edge ID: " << edgeID );
2250   }
2251 }