src/SMESHUtils/SMESH_Block.cxx

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