src/SMESH/SMESH_Block.cxx

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