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