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