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.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org
  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 <BRepTools.hxx>
  27 #include <BRepTools_WireExplorer.hxx>
  28 #include <BRep_Tool.hxx>
  29 #include <TopAbs_ShapeEnum.hxx>
  30 #include <TopExp_Explorer.hxx>
  31 #include <TopLoc_Location.hxx>
  32 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
  33 #include <TopTools_ListIteratorOfListOfShape.hxx>
  34 #include <TopTools_ListOfShape.hxx>
  35 #include <TopoDS.hxx>
  36 #include <TopoDS_Face.hxx>
  37 #include <TopoDS_Iterator.hxx>
  38 #include <TopoDS_Wire.hxx>
  39 #include <gp_Trsf.hxx>
  40 #include <gp_Vec.hxx>
  41 #include <math_FunctionSetRoot.hxx>
  42
  43 #include "utilities.h"
  44
  45 #include <list>
  46
  47 using namespace std;
  48
  49 #define SQRT_FUNC 1
  50
  51 //=======================================================================
  52 //function : SMESH_Block::TEdge::GetU
  53 //purpose  :
  54 //=======================================================================
  55
  56 double SMESH_Block::TEdge::GetU( const gp_XYZ& theParams ) const
  57 {
  58   double u = theParams.Coord( myCoordInd );
  59   return ( 1 - u ) * myFirst + u * myLast;
  60 }
  61
  62 //=======================================================================
  63 //function : SMESH_Block::TEdge::Point
  64 //purpose  :
  65 //=======================================================================
  66
  67 gp_XYZ SMESH_Block::TEdge::Point( const gp_XYZ& theParams ) const
  68 {
  69   gp_XYZ p = myC3d->Value( GetU( theParams )).XYZ();
  70   if ( myTrsf.Form() != gp_Identity )
  71     myTrsf.Transforms( p );
  72   return p;
  73 }
  74
  75 //=======================================================================
  76 //function : SMESH_Block::TFace::GetUV
  77 //purpose  :
  78 //=======================================================================
  79
  80 gp_XY SMESH_Block::TFace::GetUV( const gp_XYZ& theParams ) const
  81 {
  82   gp_XY uv(0.,0.);
  83   double dU = theParams.Coord( GetUInd() );
  84   double dV = theParams.Coord( GetVInd() );
  85   for ( int iE = 0; iE < 4; iE++ ) // loop on 4 edges
  86   {
  87     double Ecoef = 0, Vcoef = 0;
  88     switch ( iE ) {
  89     case 0:
  90       Ecoef = ( 1 - dV ); // u0
  91       Vcoef = ( 1 - dU ) * ( 1 - dV ); break; // 00
  92     case 1:
  93       Ecoef = dV; // u1
  94       Vcoef = dU * ( 1 - dV ); break; // 10
  95     case 2:
  96       Ecoef = ( 1 - dU ); // 0v
  97       Vcoef = dU * dV  ; break; // 11
  98     case 3:
  99       Ecoef = dU  ; // 1v
 100       Vcoef = ( 1 - dU ) * dV  ; break; // 01
 101     default:;
 102     }
 103     // edge addition
 104     double u = theParams.Coord( myCoordInd[ iE ] );
 105     u = ( 1 - u ) * myFirst[ iE ] + u * myLast[ iE ];
 106     uv += Ecoef * myC2d[ iE ]->Value( u ).XY();
 107     // corner addition
 108     uv -= Vcoef * myCorner[ iE ];
 109   }
 110   return uv;
 111 }
 112
 113 //=======================================================================
 114 //function : SMESH_Block::TFace::Point
 115 //purpose  :
 116 //=======================================================================
 117
 118 gp_XYZ SMESH_Block::TFace::Point( const gp_XYZ& theParams ) const
 119 {
 120   gp_XY uv = GetUV( theParams );
 121   gp_XYZ p = myS->Value( uv.X(), uv.Y() ).XYZ();
 122   if ( myTrsf.Form() != gp_Identity )
 123     myTrsf.Transforms( p );
 124   return p;
 125 }
 126
 127 //=======================================================================
 128 //function : GetShapeCoef
 129 //purpose  :
 130 //=======================================================================
 131
 132 double* SMESH_Block::GetShapeCoef (const int theShapeID)
 133 {
 134   static double shapeCoef[][3] = {
 135     //    V000,        V100,        V010,         V110
 136     { -1,-1,-1 }, {  1,-1,-1 }, { -1, 1,-1 }, {  1, 1,-1 },
 137     //    V001,        V101,        V011,         V111,
 138     { -1,-1, 1 }, {  1,-1, 1 }, { -1, 1, 1 }, {  1, 1, 1 },
 139     //    Ex00,        Ex10,        Ex01,         Ex11,
 140     {  0,-1,-1 }, {  0, 1,-1 }, {  0,-1, 1 }, {  0, 1, 1 },
 141     //    E0y0,        E1y0,        E0y1,         E1y1,
 142     { -1, 0,-1 }, {  1, 0,-1 }, { -1, 0, 1 }, {  1, 0, 1 },
 143     //    E00z,        E10z,        E01z,         E11z,
 144     { -1,-1, 0 }, {  1,-1, 0 }, { -1, 1, 0 }, {  1, 1, 0 },
 145     //    Fxy0,        Fxy1,        Fx0z,         Fx1z,         F0yz,           F1yz,
 146     {  0, 0,-1 }, {  0, 0, 1 }, {  0,-1, 0 }, {  0, 1, 0 }, { -1, 0, 0 }, {  1, 0, 0 },
 147     // ID_Shell
 148     {  0, 0, 0 }
 149   };
 150   if ( theShapeID < ID_V000 || theShapeID > ID_F1yz )
 151     return shapeCoef[ ID_Shell - 1 ];
 152
 153   return shapeCoef[ theShapeID - 1 ];
 154 }
 155
 156 //=======================================================================
 157 //function : ShellPoint
 158 //purpose  : return coordinates of a point in shell
 159 //=======================================================================
 160
 161 bool SMESH_Block::ShellPoint( const gp_XYZ& theParams, gp_XYZ& thePoint ) const
 162 {
 163   thePoint.SetCoord( 0., 0., 0. );
 164   for ( int shapeID = ID_V000; shapeID < ID_Shell; shapeID++ )
 165   {
 166     // coef
 167     double* coefs = GetShapeCoef( shapeID );
 168     double k = 1;
 169     for ( int iCoef = 0; iCoef < 3; iCoef++ ) {
 170       if ( coefs[ iCoef ] != 0 ) {
 171         if ( coefs[ iCoef ] < 0 )
 172           k *= ( 1. - theParams.Coord( iCoef + 1 ));
 173         else
 174           k *= theParams.Coord( iCoef + 1 );
 175       }
 176     }
 177     // point on a shape
 178     gp_XYZ Ps;
 179     if ( shapeID < ID_Ex00 ) // vertex
 180       VertexPoint( shapeID, Ps );
 181     else if ( shapeID < ID_Fxy0 ) { // edge
 182       EdgePoint( shapeID, theParams, Ps );
 183       k = -k;
 184     } else // face
 185       FacePoint( shapeID, theParams, Ps );
 186
 187     thePoint += k * Ps;
 188   }
 189   return true;
 190 }
 191
 192 //=======================================================================
 193 //function : NbVariables
 194 //purpose  :
 195 //=======================================================================
 196
 197 Standard_Integer SMESH_Block::NbVariables() const
 198 {
 199   return 3;
 200 }
 201
 202 //=======================================================================
 203 //function : NbEquations
 204 //purpose  :
 205 //=======================================================================
 206
 207 Standard_Integer SMESH_Block::NbEquations() const
 208 {
 209   return 1;
 210 }
 211
 212 //=======================================================================
 213 //function : Value
 214 //purpose  :
 215 //=======================================================================
 216
 217 Standard_Boolean SMESH_Block::Value(const math_Vector& theXYZ, math_Vector& theFxyz)
 218 {
 219   gp_XYZ P, params( theXYZ(1), theXYZ(2), theXYZ(3) );
 220   if ( params.IsEqual( myParam, DBL_MIN )) { // same param
 221     theFxyz( 1 ) = myValues[ 0 ];
 222   }
 223   else {
 224     ShellPoint( params, P );
 225     gp_Vec dP( P - myPoint );
 226     theFxyz(1) = SQRT_FUNC ? dP.SquareMagnitude() : dP.Magnitude();
 227   }
 228   return true;
 229 }
 230
 231 //=======================================================================
 232 //function : Derivatives
 233 //purpose  :
 234 //=======================================================================
 235
 236 Standard_Boolean SMESH_Block::Derivatives(const math_Vector& XYZ,math_Matrix& Df)
 237 {
 238   MESSAGE( "SMESH_Block::Derivatives()");
 239   math_Vector F(1,3);
 240   return Values(XYZ,F,Df);
 241 }
 242
 243 //=======================================================================
 244 //function : Values
 245 //purpose  :
 246 //=======================================================================
 247
 248 Standard_Boolean SMESH_Block::Values(const math_Vector& theXYZ,
 249                                      math_Vector&       theFxyz,
 250                                      math_Matrix&       theDf)
 251 {
 252 //  MESSAGE( endl<<"SMESH_Block::Values( "<<theXYZ(1)<<", "<<theXYZ(2)<<", "<<theXYZ(3)<<")");
 253
 254   gp_XYZ P, params( theXYZ(1), theXYZ(2), theXYZ(3) );
 255   if ( params.IsEqual( myParam, DBL_MIN )) { // same param
 256     theFxyz( 1 ) = myValues[ 0 ];
 257     theDf( 1,1 ) = myValues[ 1 ];
 258     theDf( 1,2 ) = myValues[ 2 ];
 259     theDf( 1,3 ) = myValues[ 3 ];
 260     return true;
 261   }
 262
 263   ShellPoint( params, P );
 264   //myNbIterations++; // how many time call ShellPoint()
 265
 266   gp_Vec dP( P - myPoint );
 267   theFxyz(1) = SQRT_FUNC ? dP.SquareMagnitude() : dP.Magnitude();
 268   if ( theFxyz(1) < 1e-6 ) {
 269     myParam      = params;
 270     myValues[ 0 ]= 0;
 271     theDf( 1,1 ) = 0;
 272     theDf( 1,2 ) = 0;
 273     theDf( 1,3 ) = 0;
 274     return true;
 275   }
 276
 277   if ( theFxyz(1) < myValues[0] )
 278   {
 279     // 3 partial derivatives
 280     gp_Vec drv[ 3 ];
 281     for ( int iP = 1; iP <= 3; iP++ ) {
 282       gp_XYZ Pi;
 283       params.SetCoord( iP, theXYZ( iP ) + 0.001 );
 284       ShellPoint( params, Pi );
 285       params.SetCoord( iP, theXYZ( iP ) ); // restore params
 286       gp_Vec dPi ( P, Pi );
 287       double mag = dPi.Magnitude();
 288       if ( mag > DBL_MIN )
 289         dPi /= mag;
 290       drv[ iP - 1 ] = dPi;
 291     }
 292     for ( int iP = 0; iP < 3; iP++ ) {
 293       if ( iP == myFaceIndex )
 294         theDf( 1, iP + 1 ) = myFaceParam;
 295       else {
 296         // like IntAna_IntConicQuad::Perform (const gp_Lin& L, const gp_Pln& P)
 297         // where L is (P -> myPoint), P is defined by the 2 other derivative direction
 298         int iPrev = ( iP ? iP - 1 : 2 );
 299         int iNext = ( iP == 2 ? 0 : iP + 1 );
 300         gp_Vec plnNorm = drv[ iPrev ].Crossed( drv [ iNext ] );
 301         double Direc = plnNorm * drv[ iP ];
 302         if ( Abs(Direc) <= DBL_MIN )
 303           theDf( 1, iP + 1 ) = dP * drv[ iP ];
 304         else {
 305           double Dis = plnNorm * P - plnNorm * myPoint;
 306           theDf( 1, iP + 1 ) = Dis/Direc;
 307         }
 308       }
 309     }
 310     //myNbIterations +=3; // how many time call ShellPoint()
 311
 312     // store better values
 313     myParam    = params;
 314     myValues[0]= theFxyz(1);
 315     myValues[1]= theDf(1,1);
 316     myValues[2]= theDf(1,2);
 317     myValues[3]= theDf(1,3);
 318
 319 //     SCRUTE( theFxyz(1)  );
 320 //     SCRUTE( theDf( 1,1 ));
 321 //     SCRUTE( theDf( 1,2 ));
 322 //     SCRUTE( theDf( 1,3 ));
 323   }
 324
 325   return true;
 326 }
 327
 328 //=======================================================================
 329 //function : ComputeParameters
 330 //purpose  : compute point parameters in the block
 331 //=======================================================================
 332
 333 bool SMESH_Block::ComputeParameters(const gp_Pnt& thePoint,
 334                                     gp_XYZ&       theParams,
 335                                     const int     theShapeID)
 336 {
 337 //   MESSAGE( endl<<"SMESH_Block::ComputeParameters( "
 338 //           <<thePoint.X()<<", "<<thePoint.Y()<<", "<<thePoint.Z()<<")");
 339
 340   myPoint = thePoint.XYZ();
 341
 342   myParam.SetCoord( -1,-1,-1 );
 343   myValues[0] = 1e100;
 344
 345   const bool isOnFace = IsFaceID( theShapeID );
 346   double * coef = GetShapeCoef( theShapeID );
 347
 348   // the first guess
 349   math_Vector start( 1, 3, 0.0 );
 350   if ( !myGridComputed )
 351   {
 352     // define the first guess by thePoint projection on lines
 353     // connecting vertices
 354     bool needGrid = false;
 355     gp_XYZ par000( 0, 0, 0 ), par111( 1, 1, 1 );
 356     double zero = DBL_MIN * DBL_MIN;
 357     for ( int iEdge = 0, iParam = 1; iParam <= 3 && !needGrid; iParam++ )
 358     {
 359       if ( isOnFace && coef[ iParam - 1 ] != 0 ) {
 360         iEdge += 4;
 361         continue;
 362       }
 363       for ( int iE = 0; iE < 4; iE++, iEdge++ ) { // loop on 4 parallel edges
 364         gp_Pnt p0 = myEdge[ iEdge ].Point( par000 );
 365         gp_Pnt p1 = myEdge[ iEdge ].Point( par111 );
 366         gp_Vec v01( p0, p1 ), v0P( p0, thePoint );
 367         double len2 = v01.SquareMagnitude();
 368         double par = 0;
 369         if ( len2 > zero ) {
 370           par = v0P.Dot( v01 ) / len2;
 371           if ( par < 0 || par > 1 ) {
 372             needGrid = true;
 373             break;
 374           }
 375         }
 376         start( iParam ) += par;
 377       }
 378       start( iParam ) /= 4.;
 379     }
 380     if ( needGrid ) {
 381       // compute nodes of 3 x 3 x 3 grid
 382       int iNode = 0;
 383       for ( double x = 0.25; x < 0.9; x += 0.25 )
 384         for ( double y = 0.25; y < 0.9; y += 0.25 )
 385           for ( double z = 0.25; z < 0.9; z += 0.25 ) {
 386             TxyzPair & prmPtn = my3x3x3GridNodes[ iNode++ ];
 387             prmPtn.first.SetCoord( x, y, z );
 388             ShellPoint( prmPtn.first, prmPtn.second );
 389           }
 390       myGridComputed = true;
 391     }
 392   }
 393   if ( myGridComputed ) {
 394     double minDist = DBL_MAX;
 395     gp_XYZ* bestParam = 0;
 396     for ( int iNode = 0; iNode < 27; iNode++ ) {
 397       TxyzPair & prmPtn = my3x3x3GridNodes[ iNode ];
 398       double dist = ( thePoint.XYZ() - prmPtn.second ).SquareModulus();
 399       if ( dist < minDist ) {
 400         minDist = dist;
 401         bestParam = & prmPtn.first;
 402       }
 403     }
 404     start( 1 ) = bestParam->X();
 405     start( 2 ) = bestParam->Y();
 406     start( 3 ) = bestParam->Z();
 407   }
 408
 409   int myFaceIndex = -1;
 410   if ( isOnFace ) {
 411     // put a point on the face
 412     for ( int iCoord = 0; iCoord < 3; iCoord++ )
 413       if ( coef[ iCoord ] ) {
 414         myFaceIndex = iCoord;
 415         myFaceParam = ( coef[ myFaceIndex ] < 0.5 ) ? 0.0 : 1.0;
 416         start( iCoord + 1 ) = myFaceParam;
 417       }
 418   }
 419   math_Vector low  ( 1, 3, 0.0 );
 420   math_Vector up   ( 1, 3, 1.0 );
 421   math_Vector tol  ( 1, 3, 1e-4 );
 422   math_FunctionSetRoot paramSearch( *this, tol );
 423
 424   int nbLoops = 0;
 425   while ( myValues[0] > 1e-1 && nbLoops++ < 10 ) {
 426     paramSearch.Perform ( *this, start, low, up );
 427     if ( !paramSearch.IsDone() ) {
 428       //MESSAGE( " !paramSearch.IsDone() " );
 429     }
 430     else {
 431       //MESSAGE( " NB ITERATIONS: " << paramSearch.NbIterations() );
 432     }
 433     start( 1 ) = myParam.X();
 434     start( 2 ) = myParam.Y();
 435     start( 3 ) = myParam.Z();
 436     //MESSAGE( "Distance: " << ( SQRT_FUNC ? sqrt(myValues[0]) : myValues[0] ));
 437   }
 438 //   MESSAGE( endl << myParam.X() << " " << myParam.Y() << " " << myParam.Z() << endl);
 439 //   mySumDist += myValues[0];
 440 //   MESSAGE( " TOTAL NB ITERATIONS: " << myNbIterations <<
 441 //             " DIST: " << ( SQRT_FUNC ? sqrt(mySumDist) : mySumDist ));
 442
 443
 444   theParams = myParam;
 445
 446   return true;
 447 }
 448
 449 //=======================================================================
 450 //function : GetStateNumber
 451 //purpose  :
 452 //=======================================================================
 453
 454 Standard_Integer SMESH_Block::GetStateNumber ()
 455 {
 456 //   MESSAGE( endl<<"SMESH_Block::GetStateNumber( "<<myParam.X()<<", "<<
 457 //           myParam.Y()<<", "<<myParam.Z()<<") DISTANCE: " << myValues[0]);
 458   return myValues[0] < 1e-1;
 459 }
 460
 461 //=======================================================================
 462 //function : DumpShapeID
 463 //purpose  : debug an id of a block sub-shape
 464 //=======================================================================
 465
 466 #define CASEDUMP(id,strm) case id: strm << #id; break;
 467
 468 ostream& SMESH_Block::DumpShapeID (const int id, ostream& stream)
 469 {
 470   switch ( id ) {
 471   CASEDUMP( ID_V000, stream );
 472   CASEDUMP( ID_V100, stream );
 473   CASEDUMP( ID_V010, stream );
 474   CASEDUMP( ID_V110, stream );
 475   CASEDUMP( ID_V001, stream );
 476   CASEDUMP( ID_V101, stream );
 477   CASEDUMP( ID_V011, stream );
 478   CASEDUMP( ID_V111, stream );
 479   CASEDUMP( ID_Ex00, stream );
 480   CASEDUMP( ID_Ex10, stream );
 481   CASEDUMP( ID_Ex01, stream );
 482   CASEDUMP( ID_Ex11, stream );
 483   CASEDUMP( ID_E0y0, stream );
 484   CASEDUMP( ID_E1y0, stream );
 485   CASEDUMP( ID_E0y1, stream );
 486   CASEDUMP( ID_E1y1, stream );
 487   CASEDUMP( ID_E00z, stream );
 488   CASEDUMP( ID_E10z, stream );
 489   CASEDUMP( ID_E01z, stream );
 490   CASEDUMP( ID_E11z, stream );
 491   CASEDUMP( ID_Fxy0, stream );
 492   CASEDUMP( ID_Fxy1, stream );
 493   CASEDUMP( ID_Fx0z, stream );
 494   CASEDUMP( ID_Fx1z, stream );
 495   CASEDUMP( ID_F0yz, stream );
 496   CASEDUMP( ID_F1yz, stream );
 497   CASEDUMP( ID_Shell, stream );
 498   default: stream << "ID_INVALID";
 499   }
 500   return stream;
 501 }
 502
 503 //=======================================================================
 504 //function : GetShapeIDByParams
 505 //purpose  : define an id of the block sub-shape by normlized point coord
 506 //=======================================================================
 507
 508 int SMESH_Block::GetShapeIDByParams ( const gp_XYZ& theCoord )
 509 {
 510   //   id ( 0 - 26 ) computation:
 511
 512   //   vertex     ( 0 - 7 )  : id = 1*x + 2*y + 4*z
 513
 514   //   edge || X  ( 8 - 11 ) : id = 8   + 1*y + 2*z
 515   //   edge || Y  ( 12 - 15 ): id = 1*x + 12  + 2*z
 516   //   edge || Z  ( 16 - 19 ): id = 1*x + 2*y + 16
 517
 518   //   face || XY ( 20 - 21 ): id = 8   + 12  + 1*z - 0
 519   //   face || XZ ( 22 - 23 ): id = 8   + 1*y + 16  - 2
 520   //   face || YZ ( 24 - 25 ): id = 1*x + 12  + 16  - 4
 521
 522   static int iAddBnd[]    = { 1, 2, 4 };
 523   static int iAddNotBnd[] = { 8, 12, 16 };
 524   static int iFaceSubst[] = { 0, 2, 4 };
 525
 526   int id = 0;
 527   int iOnBoundary = 0;
 528   for ( int iCoord = 0; iCoord < 3; iCoord++ )
 529   {
 530     double val = theCoord.Coord( iCoord + 1 );
 531     if ( val == 0.0 )
 532       iOnBoundary++;
 533     else if ( val == 1.0 )
 534       id += iAddBnd[ iOnBoundary++ ];
 535     else
 536       id += iAddNotBnd[ iCoord ];
 537   }
 538   if ( iOnBoundary == 1 ) // face
 539     id -= iFaceSubst[ (id - 20) / 4 ];
 540   else if ( iOnBoundary == 0 ) // shell
 541     id = 26;
 542
 543   if ( id > 26 || id < 0 ) {
 544     MESSAGE( "GetShapeIDByParams() = " << id
 545             <<" "<< theCoord.X() <<" "<< theCoord.Y() <<" "<< theCoord.Z() );
 546   }
 547
 548   return id + 1; // shape ids start at 1
 549 }
 550
 551
 552 //=======================================================================
 553 //function : getOrderedEdges
 554 //purpose  : return nb wires and a list of oredered edges
 555 //=======================================================================
 556
 557 static int getOrderedEdges (const TopoDS_Face&   theFace,
 558                             const TopoDS_Vertex& theFirstVertex,
 559                             list< TopoDS_Edge >& theEdges,
 560                             list< int >  &       theNbVertexInWires)
 561 {
 562   // put wires in a list, so that an outer wire comes first
 563   list<TopoDS_Wire> aWireList;
 564   TopoDS_Wire anOuterWire = BRepTools::OuterWire( theFace );
 565   aWireList.push_back( anOuterWire );
 566   for ( TopoDS_Iterator wIt (theFace); wIt.More(); wIt.Next() )
 567     if ( !anOuterWire.IsSame( wIt.Value() ))
 568       aWireList.push_back( TopoDS::Wire( wIt.Value() ));
 569
 570   // loop on edges of wires
 571   theNbVertexInWires.clear();
 572   list<TopoDS_Wire>::iterator wlIt = aWireList.begin();
 573   for ( ; wlIt != aWireList.end(); wlIt++ )
 574   {
 575     int iE;
 576     BRepTools_WireExplorer wExp( *wlIt, theFace );
 577     for ( iE = 0; wExp.More(); wExp.Next(), iE++ )
 578     {
 579       TopoDS_Edge edge = wExp.Current();
 580       edge = TopoDS::Edge( edge.Oriented( wExp.Orientation() ));
 581       theEdges.push_back( edge );
 582     }
 583     theNbVertexInWires.push_back( iE );
 584     iE = 0;
 585     if ( wlIt == aWireList.begin() && theEdges.size() > 1 ) { // the outer wire
 586       // orient closed edges
 587       list< TopoDS_Edge >::iterator eIt, eIt2;
 588       for ( eIt = theEdges.begin(); eIt != theEdges.end(); eIt++ )
 589       {
 590         TopoDS_Edge& edge = *eIt;
 591         if ( TopExp::FirstVertex( edge ).IsSame( TopExp::LastVertex( edge ) ))
 592         {
 593           eIt2 = eIt;
 594           bool isNext = ( eIt2 == theEdges.begin() );
 595           TopoDS_Edge edge2 = isNext ? *(++eIt2) : *(--eIt2);
 596           double f1,l1,f2,l2;
 597           Handle(Geom2d_Curve) c1 = BRep_Tool::CurveOnSurface( edge, theFace, f1,l1 );
 598           Handle(Geom2d_Curve) c2 = BRep_Tool::CurveOnSurface( edge2, theFace, f2,l2 );
 599           gp_Pnt2d pf = c1->Value( edge.Orientation() == TopAbs_FORWARD ? f1 : l1 );
 600           gp_Pnt2d pl = c1->Value( edge.Orientation() == TopAbs_FORWARD ? l1 : f1 );
 601           bool isFirst = ( edge2.Orientation() == TopAbs_FORWARD ? isNext : !isNext );
 602           gp_Pnt2d p2 = c2->Value( isFirst ? f2 : l2 );
 603           isFirst = ( p2.SquareDistance( pf ) < p2.SquareDistance( pl ));
 604           if ( isNext ? isFirst : !isFirst )
 605             edge.Reverse();
 606         }
 607       }
 608       // rotate theEdges until it begins from theFirstVertex
 609       if ( ! theFirstVertex.IsNull() )
 610         while ( !theFirstVertex.IsSame( TopExp::FirstVertex( theEdges.front(), true )))
 611         {
 612           theEdges.splice(theEdges.end(), theEdges,
 613                           theEdges.begin(), ++ theEdges.begin());
 614           if ( iE++ > theNbVertexInWires.back() )
 615             break; // break infinite loop
 616         }
 617     }
 618   }
 619
 620   return aWireList.size();
 621 }
 622
 623 //=======================================================================
 624 //function : LoadBlockShapes
 625 //purpose  : add sub-shapes of theBlock to theShapeIDMap so that they get
 626 //           IDs acoording to enum TShapeID
 627 //=======================================================================
 628
 629 bool SMESH_Block::LoadBlockShapes(const TopoDS_Shell&         theShell,
 630                                   const TopoDS_Vertex&        theVertex000,
 631                                   const TopoDS_Vertex&        theVertex001,
 632 //                             TopTools_IndexedMapOfShape& theShapeIDMap
 633                                   TopTools_IndexedMapOfOrientedShape& theShapeIDMap )
 634 {
 635   MESSAGE(" ::LoadBlockShapes()");
 636
 637   myShell = theShell;
 638   myNbIterations = 0;
 639   mySumDist = 0;
 640   myGridComputed = false;
 641
 642   // 6 vertices
 643   TopoDS_Shape V000, V100, V010, V110, V001, V101, V011, V111;
 644   // 12 edges
 645   TopoDS_Shape Ex00, Ex10, Ex01, Ex11;
 646   TopoDS_Shape E0y0, E1y0, E0y1, E1y1;
 647   TopoDS_Shape E00z, E10z, E01z, E11z;
 648   // 6 faces
 649   TopoDS_Shape Fxy0, Fx0z, F0yz, Fxy1, Fx1z, F1yz;
 650
 651   // nb of faces bound to a vertex in TopTools_IndexedDataMapOfShapeListOfShape
 652   // filled by TopExp::MapShapesAndAncestors()
 653   const int NB_FACES_BY_VERTEX = 6;
 654
 655   TopTools_IndexedDataMapOfShapeListOfShape vfMap;
 656   TopExp::MapShapesAndAncestors( myShell, TopAbs_VERTEX, TopAbs_FACE, vfMap );
 657   if ( vfMap.Extent() != 8 ) {
 658     MESSAGE(" Wrong nb of vertices in the block: " << vfMap.Extent() );
 659     return false;
 660   }
 661
 662   V000 = theVertex000;
 663   V001 = theVertex001;
 664
 665   if ( V000.IsNull() ) {
 666     // find vertex 000 - the one with smallest coordinates
 667     double minVal = DBL_MAX, minX, val;
 668     for ( int i = 1; i <= 8; i++ ) {
 669       const TopoDS_Vertex& v = TopoDS::Vertex( vfMap.FindKey( i ));
 670       gp_Pnt P = BRep_Tool::Pnt( v );
 671       val = P.X() + P.Y() + P.Z();
 672       if ( val < minVal || ( val == minVal && P.X() < minX )) {
 673         V000 = v;
 674         minVal = val;
 675         minX = P.X();
 676       }
 677     }
 678     // find vertex 001 - the one on the most vertical edge passing through V000
 679     TopTools_IndexedDataMapOfShapeListOfShape veMap;
 680     TopExp::MapShapesAndAncestors( myShell, TopAbs_VERTEX, TopAbs_EDGE, veMap );
 681     gp_Vec dir001 = gp::DZ();
 682     gp_Pnt p000 = BRep_Tool::Pnt( TopoDS::Vertex( V000 ));
 683     double maxVal = -DBL_MAX;
 684     TopTools_ListIteratorOfListOfShape eIt ( veMap.FindFromKey( V000 ));
 685     for (  ; eIt.More(); eIt.Next() ) {
 686       const TopoDS_Edge& e = TopoDS::Edge( eIt.Value() );
 687       TopoDS_Vertex v = TopExp::FirstVertex( e );
 688       if ( v.IsSame( V000 ))
 689         v = TopExp::LastVertex( e );
 690       val = dir001 * gp_Vec( p000, BRep_Tool::Pnt( v )).Normalized();
 691       if ( val > maxVal ) {
 692         V001 = v;
 693         maxVal = val;
 694       }
 695     }
 696   }
 697
 698   // find the bottom (Fxy0), Fx0z and F0yz faces
 699
 700   const TopTools_ListOfShape& f000List = vfMap.FindFromKey( V000 );
 701   const TopTools_ListOfShape& f001List = vfMap.FindFromKey( V001 );
 702   if (f000List.Extent() != NB_FACES_BY_VERTEX ||
 703       f001List.Extent() != NB_FACES_BY_VERTEX ) {
 704     MESSAGE(" LoadBlockShapes() " << f000List.Extent() << " " << f001List.Extent());
 705     return false;
 706   }
 707   TopTools_ListIteratorOfListOfShape f001It, f000It ( f000List );
 708   int i, j, iFound1, iFound2;
 709   for ( j = 0; f000It.More(); f000It.Next(), j++ )
 710   {
 711     if ( NB_FACES_BY_VERTEX == 6 && j % 2 ) continue; // each face encounters twice
 712     const TopoDS_Shape& F = f000It.Value();
 713     for ( i = 0, f001It.Initialize( f001List ); f001It.More(); f001It.Next(), i++ ) {
 714       if ( NB_FACES_BY_VERTEX == 6 && i % 2 ) continue; // each face encounters twice
 715       if ( F.IsSame( f001It.Value() ))
 716         break;
 717     }
 718     if ( f001It.More() ) // Fx0z or F0yz found
 719       if ( Fx0z.IsNull() ) {
 720         Fx0z = F;
 721         iFound1 = i;
 722       } else {
 723         F0yz = F;
 724         iFound2 = i;
 725       }
 726     else // F is the bottom face
 727       Fxy0 = F;
 728   }
 729   if ( Fxy0.IsNull() || Fx0z.IsNull() || F0yz.IsNull() ) {
 730     MESSAGE( Fxy0.IsNull() <<" "<< Fx0z.IsNull() <<" "<< F0yz.IsNull() );
 731     return false;
 732   }
 733
 734   // choose the top face (Fxy1)
 735   for ( i = 0, f001It.Initialize( f001List ); f001It.More(); f001It.Next(), i++ ) {
 736     if ( NB_FACES_BY_VERTEX == 6 && i % 2 ) continue; // each face encounters twice
 737     if ( i != iFound1 && i != iFound2 )
 738       break;
 739   }
 740   Fxy1 = f001It.Value();
 741   if ( Fxy1.IsNull() ) {
 742     MESSAGE(" LoadBlockShapes() error ");
 743     return false;
 744   }
 745
 746   // find bottom edges and veritices
 747   list< TopoDS_Edge > eList;
 748   list< int >         nbVertexInWires;
 749   getOrderedEdges( TopoDS::Face( Fxy0 ), TopoDS::Vertex( V000 ), eList, nbVertexInWires );
 750   if ( nbVertexInWires.size() != 1 || nbVertexInWires.front() != 4 ) {
 751     MESSAGE(" LoadBlockShapes() error ");
 752     return false;
 753   }
 754   list< TopoDS_Edge >::iterator elIt = eList.begin();
 755   for ( i = 0; elIt != eList.end(); elIt++, i++ )
 756     switch ( i ) {
 757     case 0: E0y0 = *elIt; V010 = TopExp::LastVertex( *elIt, true ); break;
 758     case 1: Ex10 = *elIt; V110 = TopExp::LastVertex( *elIt, true ); break;
 759     case 2: E1y0 = *elIt; V100 = TopExp::LastVertex( *elIt, true ); break;
 760     case 3: Ex00 = *elIt; break;
 761     default:;
 762     }
 763   if ( i != 4 || E0y0.IsNull() || Ex10.IsNull() || E1y0.IsNull() || Ex00.IsNull() ) {
 764     MESSAGE(" LoadBlockShapes() error, eList.size()=" << eList.size());
 765     return false;
 766   }
 767
 768
 769   // find top edges and veritices
 770   eList.clear();
 771   getOrderedEdges( TopoDS::Face( Fxy1 ), TopoDS::Vertex( V001 ), eList, nbVertexInWires );
 772   if ( nbVertexInWires.size() != 1 || nbVertexInWires.front() != 4 ) {
 773     MESSAGE(" LoadBlockShapes() error ");
 774     return false;
 775   }
 776   for ( i = 0, elIt = eList.begin(); elIt != eList.end(); elIt++, i++ )
 777     switch ( i ) {
 778     case 0: Ex01 = *elIt; V101 = TopExp::LastVertex( *elIt, true ); break;
 779     case 1: E1y1 = *elIt; V111 = TopExp::LastVertex( *elIt, true ); break;
 780     case 2: Ex11 = *elIt; V011 = TopExp::LastVertex( *elIt, true ); break;
 781     case 3: E0y1 = *elIt; break;
 782     default:;
 783     }
 784   if ( i != 4 || Ex01.IsNull() || E1y1.IsNull() || Ex11.IsNull() || E0y1.IsNull() ) {
 785     MESSAGE(" LoadBlockShapes() error, eList.size()=" << eList.size());
 786     return false;
 787   }
 788
 789   // swap Fx0z and F0yz if necessary
 790   TopExp_Explorer exp( Fx0z, TopAbs_VERTEX );
 791   for ( ; exp.More(); exp.Next() ) // Fx0z shares V101 and V100
 792     if ( V101.IsSame( exp.Current() ) || V100.IsSame( exp.Current() ))
 793       break; // V101 or V100 found
 794   if ( !exp.More() ) { // not found
 795     TopoDS_Shape f = Fx0z; Fx0z = F0yz; F0yz = f;
 796   }
 797
 798   // find Fx1z and F1yz faces
 799   const TopTools_ListOfShape& f111List = vfMap.FindFromKey( V111 );
 800   const TopTools_ListOfShape& f110List = vfMap.FindFromKey( V110 );
 801   if (f111List.Extent() != NB_FACES_BY_VERTEX ||
 802       f110List.Extent() != NB_FACES_BY_VERTEX ) {
 803     MESSAGE(" LoadBlockShapes() " << f111List.Extent() << " " << f110List.Extent());
 804     return false;
 805   }
 806   TopTools_ListIteratorOfListOfShape f111It, f110It ( f110List);
 807   for ( j = 0 ; f110It.More(); f110It.Next(), j++ ) {
 808     if ( NB_FACES_BY_VERTEX == 6 && j % 2 ) continue; // each face encounters twice
 809     const TopoDS_Shape& F = f110It.Value();
 810     for ( i = 0, f111It.Initialize( f111List ); f111It.More(); f111It.Next(), i++ ) {
 811       if ( NB_FACES_BY_VERTEX == 6 && i % 2 ) continue; // each face encounters twice
 812       if ( F.IsSame( f111It.Value() )) { // Fx1z or F1yz found
 813         if ( Fx1z.IsNull() )
 814           Fx1z = F;
 815         else
 816           F1yz = F;
 817       }
 818     }
 819   }
 820   if ( Fx1z.IsNull() || F1yz.IsNull() ) {
 821     MESSAGE(" LoadBlockShapes() error ");
 822     return false;
 823   }
 824
 825   // swap Fx1z and F1yz if necessary
 826   for ( exp.Init( Fx1z, TopAbs_VERTEX ); exp.More(); exp.Next() )
 827     if ( V010.IsSame( exp.Current() ) || V011.IsSame( exp.Current() ))
 828       break;
 829   if ( !exp.More() ) {
 830     TopoDS_Shape f = Fx1z; Fx1z = F1yz; F1yz = f;
 831   }
 832
 833   // find vertical edges
 834   for ( exp.Init( Fx0z, TopAbs_EDGE ); exp.More(); exp.Next() ) {
 835     const TopoDS_Edge& edge = TopoDS::Edge( exp.Current() );
 836     const TopoDS_Shape& vFirst = TopExp::FirstVertex( edge, true );
 837     if ( vFirst.IsSame( V001 ))
 838       E00z = edge;
 839     else if ( vFirst.IsSame( V100 ))
 840       E10z = edge;
 841   }
 842   if ( E00z.IsNull() || E10z.IsNull() ) {
 843     MESSAGE(" LoadBlockShapes() error ");
 844     return false;
 845   }
 846   for ( exp.Init( Fx1z, TopAbs_EDGE ); exp.More(); exp.Next() ) {
 847     const TopoDS_Edge& edge = TopoDS::Edge( exp.Current() );
 848     const TopoDS_Shape& vFirst = TopExp::FirstVertex( edge, true );
 849     if ( vFirst.IsSame( V111 ))
 850       E11z = edge;
 851     else if ( vFirst.IsSame( V010 ))
 852       E01z = edge;
 853   }
 854   if ( E01z.IsNull() || E11z.IsNull() ) {
 855     MESSAGE(" LoadBlockShapes() error ");
 856     return false;
 857   }
 858
 859   // load shapes in theShapeIDMap
 860
 861   theShapeIDMap.Clear();
 862
 863   theShapeIDMap.Add(V000.Oriented( TopAbs_FORWARD ));
 864   theShapeIDMap.Add(V100.Oriented( TopAbs_FORWARD ));
 865   theShapeIDMap.Add(V010.Oriented( TopAbs_FORWARD ));
 866   theShapeIDMap.Add(V110.Oriented( TopAbs_FORWARD ));
 867   theShapeIDMap.Add(V001.Oriented( TopAbs_FORWARD ));
 868   theShapeIDMap.Add(V101.Oriented( TopAbs_FORWARD ));
 869   theShapeIDMap.Add(V011.Oriented( TopAbs_FORWARD ));
 870   theShapeIDMap.Add(V111.Oriented( TopAbs_FORWARD ));
 871
 872   theShapeIDMap.Add(Ex00);
 873   theShapeIDMap.Add(Ex10);
 874   theShapeIDMap.Add(Ex01);
 875   theShapeIDMap.Add(Ex11);
 876
 877   theShapeIDMap.Add(E0y0);
 878   theShapeIDMap.Add(E1y0);
 879   theShapeIDMap.Add(E0y1);
 880   theShapeIDMap.Add(E1y1);
 881
 882   theShapeIDMap.Add(E00z);
 883   theShapeIDMap.Add(E10z);
 884   theShapeIDMap.Add(E01z);
 885   theShapeIDMap.Add(E11z);
 886
 887   theShapeIDMap.Add(Fxy0);
 888   theShapeIDMap.Add(Fxy1);
 889   theShapeIDMap.Add(Fx0z);
 890   theShapeIDMap.Add(Fx1z);
 891   theShapeIDMap.Add(F0yz);
 892   theShapeIDMap.Add(F1yz);
 893
 894   theShapeIDMap.Add(myShell);
 895
 896   if ( theShapeIDMap.Extent() != 27 ) {
 897     MESSAGE("LoadBlockShapes() " << theShapeIDMap.Extent() );
 898     return false;
 899   }
 900
 901   // store shapes geometry
 902   for ( int shapeID = 1; shapeID < theShapeIDMap.Extent(); shapeID++ )
 903   {
 904     const TopoDS_Shape& S = theShapeIDMap( shapeID );
 905     switch ( S.ShapeType() )
 906     {
 907     case TopAbs_VERTEX: {
 908
 909       if ( shapeID > ID_V111 ) {
 910         MESSAGE(" shapeID =" << shapeID );
 911         return false;
 912       }
 913       myPnt[ shapeID - ID_V000 ] =
 914         BRep_Tool::Pnt( TopoDS::Vertex( S )).XYZ();
 915       break;
 916     }
 917     case TopAbs_EDGE: {
 918
 919       const TopoDS_Edge& edge = TopoDS::Edge( S );
 920       if ( shapeID < ID_Ex00 || shapeID > ID_E11z || edge.IsNull() ) {
 921         MESSAGE(" shapeID =" << shapeID );
 922         return false;
 923       }
 924       TEdge& tEdge = myEdge[ shapeID - ID_Ex00 ];
 925       tEdge.myCoordInd = GetCoordIndOnEdge( shapeID );
 926       TopLoc_Location loc;
 927       tEdge.myC3d = BRep_Tool::Curve( edge, loc, tEdge.myFirst, tEdge.myLast );
 928       if ( !IsForwardEdge( edge, theShapeIDMap ))
 929         Swap( tEdge.myFirst, tEdge.myLast );
 930       tEdge.myTrsf = loc;
 931       break;
 932     }
 933     case TopAbs_FACE: {
 934
 935       const TopoDS_Face& face = TopoDS::Face( S );
 936       if ( shapeID < ID_Fxy0 || shapeID > ID_F1yz || face.IsNull() ) {
 937         MESSAGE(" shapeID =" << shapeID );
 938         return false;
 939       }
 940       TFace& tFace = myFace[ shapeID - ID_Fxy0 ];
 941       // pcurves
 942       vector< int > edgeIdVec(4, -1);
 943       GetFaceEdgesIDs( shapeID, edgeIdVec );
 944       for ( int iE = 0; iE < 4; iE++ ) // loop on 4 edges
 945       {
 946         const TopoDS_Edge& edge = TopoDS::Edge( theShapeIDMap( edgeIdVec[ iE ]));
 947         tFace.myCoordInd[ iE ] = GetCoordIndOnEdge( edgeIdVec[ iE ] );
 948         tFace.myC2d[ iE ] =
 949           BRep_Tool::CurveOnSurface( edge, face, tFace.myFirst[iE], tFace.myLast[iE] );
 950         if ( !IsForwardEdge( edge, theShapeIDMap ))
 951           Swap( tFace.myFirst[ iE ], tFace.myLast[ iE ] );
 952       }
 953       // 2d corners
 954       tFace.myCorner[ 0 ] = tFace.myC2d[ 0 ]->Value( tFace.myFirst[0] ).XY();
 955       tFace.myCorner[ 1 ] = tFace.myC2d[ 0 ]->Value( tFace.myLast[0] ).XY();
 956       tFace.myCorner[ 2 ] = tFace.myC2d[ 1 ]->Value( tFace.myLast[1] ).XY();
 957       tFace.myCorner[ 3 ] = tFace.myC2d[ 1 ]->Value( tFace.myFirst[1] ).XY();
 958       // sufrace
 959       TopLoc_Location loc;
 960       tFace.myS = BRep_Tool::Surface( face, loc );
 961       tFace.myTrsf = loc;
 962       break;
 963     }
 964     default: break;
 965     }
 966   } // loop on shapes in theShapeIDMap
 967
 968   return true;
 969 }
 970
 971 //=======================================================================
 972 //function : GetFaceEdgesIDs
 973 //purpose  : return edges IDs in the order u0, u1, 0v, 1v
 974 //           u0 means "|| u, v == 0"
 975 //=======================================================================
 976
 977 void SMESH_Block::GetFaceEdgesIDs (const int faceID, vector< int >& edgeVec )
 978 {
 979   switch ( faceID ) {
 980   case ID_Fxy0:
 981     edgeVec[ 0 ] = ID_Ex00;
 982     edgeVec[ 1 ] = ID_Ex10;
 983     edgeVec[ 2 ] = ID_E0y0;
 984     edgeVec[ 3 ] = ID_E1y0;
 985     break;
 986   case ID_Fxy1:
 987     edgeVec[ 0 ] = ID_Ex01;
 988     edgeVec[ 1 ] = ID_Ex11;
 989     edgeVec[ 2 ] = ID_E0y1;
 990     edgeVec[ 3 ] = ID_E1y1;
 991     break;
 992   case ID_Fx0z:
 993     edgeVec[ 0 ] = ID_Ex00;
 994     edgeVec[ 1 ] = ID_Ex01;
 995     edgeVec[ 2 ] = ID_E00z;
 996     edgeVec[ 3 ] = ID_E10z;
 997     break;
 998   case ID_Fx1z:
 999     edgeVec[ 0 ] = ID_Ex10;
1000     edgeVec[ 1 ] = ID_Ex11;
1001     edgeVec[ 2 ] = ID_E01z;
1002     edgeVec[ 3 ] = ID_E11z;
1003     break;
1004   case ID_F0yz:
1005     edgeVec[ 0 ] = ID_E0y0;
1006     edgeVec[ 1 ] = ID_E0y1;
1007     edgeVec[ 2 ] = ID_E00z;
1008     edgeVec[ 3 ] = ID_E01z;
1009     break;
1010   case ID_F1yz:
1011     edgeVec[ 0 ] = ID_E1y0;
1012     edgeVec[ 1 ] = ID_E1y1;
1013     edgeVec[ 2 ] = ID_E10z;
1014     edgeVec[ 3 ] = ID_E11z;
1015     break;
1016   default:
1017     MESSAGE(" GetFaceEdgesIDs(), wrong face ID: " << faceID );
1018   }
1019 }