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