src/SMESHUtils/SMESH_Block.cxx

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