src/SMESHUtils/SMESH_Block.cxx

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