src/SMDS/SMDS_VolumeTool.cxx

   1 // File      : SMDS_VolumeTool.cxx
   2 // Created   : Tue Jul 13 12:22:13 2004
   3 // Author    : Edward AGAPOV (eap)
   4 // Copyright : Open CASCADE
   5
   6 #ifdef _MSC_VER
   7 #pragma warning(disable:4786)
   8 #endif
   9
  10 #include "SMDS_VolumeTool.hxx"
  11
  12 #include "SMDS_MeshElement.hxx"
  13 #include "SMDS_MeshNode.hxx"
  14 #include <map>
  15 #include <float.h>
  16 #include <math.h>
  17
  18 using namespace std;
  19
  20 // ======================================================
  21 // Node indices in faces depending on volume orientation
  22 // making most faces normals external
  23 // ======================================================
  24
  25 /*
  26 //           N3
  27 //           +
  28 //          /|\
  29 //         / | \
  30 //        /  |  \
  31 //    N0 +---|---+ N1                TETRAHEDRON
  32 //       \   |   /
  33 //        \  |  /
  34 //         \ | /
  35 //          \|/
  36 //           +
  37 //           N2
  38 */
  39 static int Tetra_F [4][4] = { // FORWARD == EXTERNAL
  40   { 0, 1, 2, 0 },              // All faces have external normals
  41   { 0, 3, 1, 0 },
  42   { 1, 3, 2, 1 },
  43   { 0, 2, 3, 0 }};
  44 static int Tetra_R [4][4] = { // REVERSED
  45   { 0, 1, 2, 0 },             // All faces but a bottom have external normals
  46   { 0, 1, 3, 0 },
  47   { 1, 2, 3, 1 },
  48   { 0, 3, 2, 0 }};
  49 static int Tetra_RE [4][4] = { // REVERSED -> FORWARD (EXTERNAL)
  50   { 0, 2, 1, 0 },              // All faces have external normals
  51   { 0, 1, 3, 0 },
  52   { 1, 2, 3, 1 },
  53   { 0, 3, 2, 0 }};
  54 static int Tetra_nbN [] = { 3, 3, 3, 3 };
  55
  56 //
  57 //     PYRAMID
  58 //
  59 static int Pyramid_F [5][5] = { // FORWARD == EXTERNAL
  60   { 0, 1, 2, 3, 0 },            // All faces have external normals
  61   { 0, 4, 1, 0, 4 },
  62   { 1, 4, 2, 1, 4 },
  63   { 2, 4, 3, 2, 4 },
  64   { 3, 4, 0, 3, 4 }};
  65 static int Pyramid_R [5][5] = { // REVERSED
  66   { 0, 1, 2, 3, 0 },            // All faces but a bottom have external normals
  67   { 0, 1, 4, 0, 4 },
  68   { 1, 2, 4, 1, 4 },
  69   { 2, 3, 4, 2, 4 },
  70   { 3, 0, 4, 3, 4 }};
  71 static int Pyramid_RE [5][5] = { // REVERSED -> FORWARD (EXTERNAL)
  72   { 0, 3, 2, 1, 0 },             // All faces but a bottom have external normals
  73   { 0, 1, 4, 0, 4 },
  74   { 1, 2, 4, 1, 4 },
  75   { 2, 3, 4, 2, 4 },
  76   { 3, 0, 4, 3, 4 }};
  77 static int Pyramid_nbN [] = { 4, 3, 3, 3, 3 };
  78
  79 /*
  80 //            + N4
  81 //           /|\
  82 //          / | \
  83 //         /  |  \
  84 //        /   |   \
  85 //    N3 +---------+ N5
  86 //       |    |    |
  87 //       |    + N1 |
  88 //       |   / \   |                PENTAHEDRON
  89 //       |  /   \  |
  90 //       | /     \ |
  91 //       |/       \|
  92 //    N0 +---------+ N2
  93 */
  94 static int Penta_F [5][5] = { // FORWARD
  95   { 0, 1, 2, 0, 0 },          // Top face has an internal normal, other - external
  96   { 3, 4, 5, 3, 3 },          // 0 is bottom, 1 is top face
  97   { 0, 2, 5, 3, 0 },
  98   { 1, 4, 5, 2, 1 },
  99   { 0, 3, 4, 1, 0 }};
 100 static int Penta_R [5][5] = { // REVERSED
 101   { 0, 1, 2, 0, 0 },          // Bottom face has an internal normal, other - external
 102   { 3, 4, 5, 3, 3 },          // 0 is bottom, 1 is top face
 103   { 0, 3, 5, 2, 0 },
 104   { 1, 2, 5, 4, 1 },
 105   { 0, 1, 4, 3, 0 }};
 106 static int Penta_FE [5][5] = { // FORWARD -> EXTERNAL
 107   { 0, 1, 2, 0, 0 },
 108   { 3, 5, 4, 3, 3 },
 109   { 0, 2, 5, 3, 0 },
 110   { 1, 4, 5, 2, 1 },
 111   { 0, 3, 4, 1, 0 }};
 112 static int Penta_RE [5][5] = { // REVERSED -> EXTERNAL
 113   { 0, 2, 1, 0, 0 },
 114   { 3, 4, 5, 3, 3 },
 115   { 0, 3, 5, 2, 0 },
 116   { 1, 2, 5, 4, 1 },
 117   { 0, 1, 4, 3, 0 }};
 118 static int Penta_nbN [] = { 3, 3, 4, 4, 4 };
 119
 120 /*
 121 //         N5+----------+N6
 122 //          /|         /|
 123 //         / |        / |
 124 //        /  |       /  |
 125 //     N4+----------+N7 |
 126 //       |   |      |   |           HEXAHEDRON
 127 //       |   |      |   |
 128 //       |   |      |   |
 129 //       | N1+------|---+N2
 130 //       |  /       |  /
 131 //       | /        | /
 132 //       |/         |/
 133 //     N0+----------+N3
 134 */
 135 static int Hexa_F [6][5] = { // FORWARD
 136   { 0, 1, 2, 3, 0 },         // opposite faces are neighbouring,
 137   { 4, 5, 6, 7, 4 },         // odd face(1,3,5) normal is internal, even(0,2,4) - external
 138   { 1, 0, 4, 5, 1 },         // same index nodes of opposite faces are linked
 139   { 2, 3, 7, 6, 2 },
 140   { 0, 3, 7, 4, 0 },
 141   { 1, 2, 6, 5, 1 }};
 142 // static int Hexa_R [6][5] = { // REVERSED
 143 //   { 0, 3, 2, 1, 0 },         // opposite faces are neighbouring,
 144 //   { 4, 7, 6, 5, 4 },         // odd face(1,3,5) normal is external, even(0,2,4) - internal
 145 //   { 1, 5, 4, 0, 1 },         // same index nodes of opposite faces are linked
 146 //   { 2, 6, 7, 3, 2 },
 147 //   { 0, 4, 7, 3, 0 },
 148 //   { 1, 5, 6, 2, 1 }};
 149 static int Hexa_FE [6][5] = { // FORWARD -> EXTERNAL
 150   { 0, 1, 2, 3, 0 } ,         // opposite faces are neighbouring,
 151   { 4, 7, 6, 5, 4 },          // all face normals are external,
 152   { 0, 4, 5, 1, 0 },          // links in opposite faces: 0-0, 1-3, 2-2, 3-1
 153   { 3, 2, 6, 7, 3 },
 154   { 0, 3, 7, 4, 0 },
 155   { 1, 5, 6, 2, 1 }};
 156 static int Hexa_RE [6][5] = { // REVERSED -> EXTERNAL
 157   { 0, 3, 2, 1, 0 },          // opposite faces are neighbouring,
 158   { 4, 5, 6, 7, 4 },          // all face normals are external,
 159   { 0, 1, 5, 4, 0 },          // links in opposite faces: 0-0, 1-3, 2-2, 3-1
 160   { 3, 7, 6, 2, 3 },
 161   { 0, 4, 7, 3, 0 },
 162   { 1, 2, 6, 5, 1 }};
 163 static int Hexa_nbN [] = { 4, 4, 4, 4, 4, 4 };
 164
 165 // ========================================================
 166 // to perform some calculations without linkage to CASCADE
 167 // ========================================================
 168 struct XYZ {
 169   double x;
 170   double y;
 171   double z;
 172   XYZ()                               { x = 0; y = 0; z = 0; }
 173   XYZ( double X, double Y, double Z ) { x = X; y = Y; z = Z; }
 174   XYZ( const XYZ& other )             { x = other.x; y = other.y; z = other.z; }
 175   XYZ( const SMDS_MeshNode* n )       { x = n->X(); y = n->Y(); z = n->Z(); }
 176   XYZ operator-( const XYZ& other );
 177   XYZ Crossed( const XYZ& other );
 178   double Dot( const XYZ& other );
 179   double Magnitude();
 180 };
 181 XYZ XYZ::operator-( const XYZ& Right ) {
 182   return XYZ(x - Right.x, y - Right.y, z - Right.z);
 183 }
 184 XYZ XYZ::Crossed( const XYZ& Right ) {
 185   return XYZ (y * Right.z - z * Right.y,
 186               z * Right.x - x * Right.z,
 187               x * Right.y - y * Right.x);
 188 }
 189 double XYZ::Dot( const XYZ& Other ) {
 190   return(x * Other.x + y * Other.y + z * Other.z);
 191 }
 192 double XYZ::Magnitude() {
 193   return sqrt (x * x + y * y + z * z);
 194 }
 195
 196 //=======================================================================
 197 //function : SMDS_VolumeTool
 198 //purpose  :
 199 //=======================================================================
 200
 201 SMDS_VolumeTool::SMDS_VolumeTool ()
 202      : myVolume( 0 ),
 203        myVolForward( true ),
 204        myNbFaces( 0 ),
 205        myVolumeNbNodes( 0 ),
 206        myExternalFaces( false )
 207 {
 208 }
 209 //=======================================================================
 210 //function : SMDS_VolumeTool
 211 //purpose  :
 212 //=======================================================================
 213
 214 SMDS_VolumeTool::SMDS_VolumeTool (const SMDS_MeshElement* theVolume)
 215      : myExternalFaces( false )
 216 {
 217   Set( theVolume );
 218 }
 219
 220 //=======================================================================
 221 //function : SMDS_VolumeTool
 222 //purpose  :
 223 //=======================================================================
 224
 225 SMDS_VolumeTool::~SMDS_VolumeTool()
 226 {
 227 }
 228
 229 //=======================================================================
 230 //function : SetVolume
 231 //purpose  : Set volume to iterate on
 232 //=======================================================================
 233
 234 bool SMDS_VolumeTool::Set (const SMDS_MeshElement* theVolume)
 235 {
 236   myVolume = 0;
 237   myVolForward = true;
 238   myCurFace = -1;
 239   myVolumeNbNodes = 0;
 240   myNbFaces = 0;
 241   if ( theVolume && theVolume->GetType() == SMDSAbs_Volume )
 242   {
 243     myVolumeNbNodes = theVolume->NbNodes();
 244     switch ( myVolumeNbNodes ) {
 245     case 4:
 246     case 5:
 247     case 6:
 248     case 8:
 249       {
 250       myVolume = theVolume;
 251       myNbFaces = theVolume->NbFaces();
 252
 253       // set volume nodes
 254       int iNode = 0;
 255       SMDS_ElemIteratorPtr nodeIt = myVolume->nodesIterator();
 256       while ( nodeIt->more() )
 257         myVolumeNodes[ iNode++ ] = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
 258
 259       // nb nodes in each face
 260       if ( myVolumeNbNodes == 4 )
 261         myFaceNbNodes = Tetra_nbN;
 262       else if ( myVolumeNbNodes == 5 )
 263         myFaceNbNodes = Pyramid_nbN;
 264       else if ( myVolumeNbNodes == 6 )
 265         myFaceNbNodes = Penta_nbN;
 266       else
 267         myFaceNbNodes = Hexa_nbN;
 268
 269       // define volume orientation
 270       XYZ botNormal;
 271       GetFaceNormal( 0, botNormal.x, botNormal.y, botNormal.z );
 272       const SMDS_MeshNode* topNode = myVolumeNodes[ myVolumeNbNodes - 1 ];
 273       const SMDS_MeshNode* botNode = myVolumeNodes[ 0 ];
 274       XYZ upDir (topNode->X() - botNode->X(),
 275                  topNode->Y() - botNode->Y(),
 276                  topNode->Z() - botNode->Z() );
 277       myVolForward = ( botNormal.Dot( upDir ) < 0 );
 278       break;
 279     }
 280     default: myVolume = 0;
 281     }
 282   }
 283   return ( myVolume != 0 );
 284 }
 285
 286 //=======================================================================
 287 //function : GetInverseNodes
 288 //purpose  : Return nodes vector of an inverse volume
 289 //=======================================================================
 290
 291 #define SWAP_NODES(nodes,i1,i2)           \
 292 {                                         \
 293   const SMDS_MeshNode* tmp = nodes[ i1 ]; \
 294   nodes[ i1 ] = nodes[ i2 ];              \
 295   nodes[ i2 ] = tmp;                      \
 296 }
 297 void SMDS_VolumeTool::Inverse ()
 298 {
 299   if ( !myVolume ) return;
 300
 301   myVolForward = !myVolForward;
 302   myCurFace = -1;
 303
 304   // inverse top and bottom faces
 305   switch ( myVolumeNbNodes ) {
 306   case 4:
 307     SWAP_NODES( myVolumeNodes, 1, 2 );
 308     break;
 309   case 5:
 310     SWAP_NODES( myVolumeNodes, 1, 3 );
 311     break;
 312   case 6:
 313     SWAP_NODES( myVolumeNodes, 1, 2 );
 314     SWAP_NODES( myVolumeNodes, 4, 5 );
 315     break;
 316   case 8:
 317     SWAP_NODES( myVolumeNodes, 1, 3 );
 318     SWAP_NODES( myVolumeNodes, 5, 7 );
 319     break;
 320   default:;
 321   }
 322 }
 323
 324 //=======================================================================
 325 //function : GetSize
 326 //purpose  : Return element volume
 327 //=======================================================================
 328
 329 double SMDS_VolumeTool::GetSize() const
 330 {
 331   return 0;
 332 }
 333
 334 //=======================================================================
 335 //function : GetBaryCenter
 336 //purpose  :
 337 //=======================================================================
 338
 339 bool SMDS_VolumeTool::GetBaryCenter(double & X, double & Y, double & Z) const
 340 {
 341   X = Y = Z = 0.;
 342   if ( !myVolume )
 343     return false;
 344
 345   for ( int i = 0; i < myVolumeNbNodes; i++ ) {
 346     X += myVolumeNodes[ i ]->X();
 347     Y += myVolumeNodes[ i ]->Y();
 348     Z += myVolumeNodes[ i ]->Z();
 349   }
 350   X /= myVolumeNbNodes;
 351   Y /= myVolumeNbNodes;
 352   Z /= myVolumeNbNodes;
 353
 354   return true;
 355 }
 356
 357 //=======================================================================
 358 //function : SetExternalNormal
 359 //purpose  : Node order will be so that faces normals are external
 360 //=======================================================================
 361
 362 void SMDS_VolumeTool::SetExternalNormal ()
 363 {
 364   myExternalFaces = true;
 365   myCurFace = -1;
 366 }
 367
 368 //=======================================================================
 369 //function : NbFaceNodes
 370 //purpose  : Return number of nodes in the array of face nodes
 371 //=======================================================================
 372
 373 int SMDS_VolumeTool::NbFaceNodes( int faceIndex )
 374 {
 375   if ( !setFace( faceIndex ))
 376     return 0;
 377   return myFaceNbNodes[ faceIndex ];
 378 }
 379
 380 //=======================================================================
 381 //function : GetFaceNodes
 382 //purpose  : Return pointer to the array of face nodes.
 383 //           To comfort link iteration, the array
 384 //           length == NbFaceNodes( faceIndex ) + 1 and
 385 //           the last node == the first one.
 386 //=======================================================================
 387
 388 const SMDS_MeshNode** SMDS_VolumeTool::GetFaceNodes( int faceIndex )
 389 {
 390   if ( !setFace( faceIndex ))
 391     return 0;
 392   return myFaceNodes;
 393 }
 394
 395 //=======================================================================
 396 //function : GetFaceNodesIndices
 397 //purpose  : Return pointer to the array of face nodes indices
 398 //           To comfort link iteration, the array
 399 //           length == NbFaceNodes( faceIndex ) + 1 and
 400 //           the last node index == the first one.
 401 //=======================================================================
 402
 403 const int* SMDS_VolumeTool::GetFaceNodesIndices( int faceIndex )
 404 {
 405   if ( !setFace( faceIndex ))
 406     return 0;
 407   return myFaceNodeIndices;
 408 }
 409
 410 //=======================================================================
 411 //function : GetFaceNodes
 412 //purpose  : Return a set of face nodes.
 413 //=======================================================================
 414
 415 bool SMDS_VolumeTool::GetFaceNodes (int                        faceIndex,
 416                                     set<const SMDS_MeshNode*>& theFaceNodes )
 417 {
 418   if ( !setFace( faceIndex ))
 419     return false;
 420
 421   theFaceNodes.clear();
 422   int iNode, nbNode = myFaceNbNodes[ faceIndex ];
 423   for ( iNode = 0; iNode < nbNode; iNode++ )
 424     theFaceNodes.insert( myFaceNodes[ iNode ]);
 425
 426   return true;
 427 }
 428
 429 //=======================================================================
 430 //function : IsFaceExternal
 431 //purpose  : Check normal orientation of a returned face
 432 //=======================================================================
 433
 434 bool SMDS_VolumeTool::IsFaceExternal( int faceIndex )
 435 {
 436   if ( myExternalFaces || !myVolume )
 437     return true;
 438
 439   switch ( myVolumeNbNodes ) {
 440   case 4:
 441   case 5:
 442     // only the bottom of a reversed tetrahedron can be internal
 443     return ( myVolForward || faceIndex != 0 );
 444   case 6:
 445     // in a forward pentahedron, the top is internal, in a reversed one - bottom
 446     return ( myVolForward ? faceIndex != 1 : faceIndex != 0 );
 447   case 8: {
 448     // in a forward hexahedron, even face normal is external, odd - internal
 449     bool odd = faceIndex % 2;
 450     return ( myVolForward ? !odd : odd );
 451   }
 452   default:;
 453   }
 454   return false;
 455 }
 456
 457 //=======================================================================
 458 //function : GetFaceNormal
 459 //purpose  : Return a normal to a face
 460 //=======================================================================
 461
 462 bool SMDS_VolumeTool::GetFaceNormal (int faceIndex, double & X, double & Y, double & Z)
 463 {
 464   if ( !setFace( faceIndex ))
 465     return false;
 466
 467   XYZ p1 ( myFaceNodes[0] );
 468   XYZ p2 ( myFaceNodes[1] );
 469   XYZ p3 ( myFaceNodes[2] );
 470   XYZ aVec12( p2 - p1 );
 471   XYZ aVec13( p3 - p1 );
 472   XYZ cross = aVec12.Crossed( aVec13 );
 473
 474   double size = cross.Magnitude();
 475   if ( size <= DBL_MIN )
 476     return false;
 477
 478   X = cross.x / size;
 479   Y = cross.y / size;
 480   Z = cross.z / size;
 481
 482   return true;
 483 }
 484
 485
 486 //=======================================================================
 487 //function : GetFaceArea
 488 //purpose  : Return face area
 489 //=======================================================================
 490
 491 double SMDS_VolumeTool::GetFaceArea( int faceIndex )
 492 {
 493   if ( !setFace( faceIndex ))
 494     return 0;
 495
 496   XYZ p1 ( myFaceNodes[0] );
 497   XYZ p2 ( myFaceNodes[1] );
 498   XYZ p3 ( myFaceNodes[2] );
 499   XYZ aVec12( p2 - p1 );
 500   XYZ aVec13( p3 - p1 );
 501   double area = aVec12.Crossed( aVec13 ).Magnitude() * 0.5;
 502
 503   if ( myFaceNbNodes[ faceIndex ] == 4 ) {
 504     XYZ p4 ( myFaceNodes[3] );
 505     XYZ aVec14( p4 - p1 );
 506     area += aVec14.Crossed( aVec13 ).Magnitude() * 0.5;
 507   }
 508   return area;
 509 }
 510
 511 //=======================================================================
 512 //function : GetOppFaceIndex
 513 //purpose  : Return index of the opposite face if it exists, else -1.
 514 //=======================================================================
 515
 516 int SMDS_VolumeTool::GetOppFaceIndex( int faceIndex ) const
 517 {
 518   int ind = -1;
 519   if ( faceIndex >= 0 && faceIndex < NbFaces() ) {
 520     switch ( myVolumeNbNodes ) {
 521     case 6:
 522       if ( faceIndex == 0 || faceIndex == 1 )
 523         ind = 1 - faceIndex;
 524       break;
 525     case 8:
 526       ind = faceIndex + ( faceIndex % 2 ? -1 : 1 );
 527       break;
 528     default:;
 529     }
 530   }
 531   return ind;
 532 }
 533
 534 //=======================================================================
 535 //function : IsLinked
 536 //purpose  : return true if theNode1 is linked with theNode2
 537 //=======================================================================
 538
 539 bool SMDS_VolumeTool::IsLinked (const SMDS_MeshNode* theNode1,
 540                                 const SMDS_MeshNode* theNode2) const
 541 {
 542   if ( !myVolume )
 543     return false;
 544
 545   // find nodes indices
 546   int i1 = -1, i2 = -1;
 547   for ( int i = 0; i < myVolumeNbNodes; i++ ) {
 548     if ( myVolumeNodes[ i ] == theNode1 )
 549       i1 = i;
 550     else if ( myVolumeNodes[ i ] == theNode2 )
 551       i2 = i;
 552   }
 553   return IsLinked( i1, i2 );
 554 }
 555
 556 //=======================================================================
 557 //function : IsLinked
 558 //purpose  : return true if the node with theNode1Index is linked
 559 //           with the node with theNode2Index
 560 //=======================================================================
 561
 562 bool SMDS_VolumeTool::IsLinked (const int theNode1Index,
 563                                 const int theNode2Index) const
 564 {
 565   int minInd = theNode1Index < theNode2Index ? theNode1Index : theNode2Index;
 566   int maxInd = theNode1Index < theNode2Index ? theNode2Index : theNode1Index;
 567
 568   if ( minInd < 0 || maxInd > myVolumeNbNodes - 1 || maxInd == minInd )
 569     return false;
 570
 571   switch ( myVolumeNbNodes ) {
 572   case 4:
 573     return true;
 574   case 5:
 575     if ( maxInd == 4 )
 576       return true;
 577     switch ( maxInd - minInd ) {
 578     case 1:
 579     case 3: return true;
 580     default:;
 581     }
 582     break;
 583   case 6:
 584     switch ( maxInd - minInd ) {
 585     case 1: return minInd != 2;
 586     case 2: return minInd == 0 || minInd == 3;
 587     case 3: return true;
 588     default:;
 589     }
 590     break;
 591   case 8:
 592     switch ( maxInd - minInd ) {
 593     case 1: return minInd != 3;
 594     case 3: return minInd == 0 || minInd == 4;
 595     case 4: return true;
 596     default:;
 597     }
 598     break;
 599   default:;
 600   }
 601   return false;
 602 }
 603
 604 //=======================================================================
 605 //function : GetNodeIndex
 606 //purpose  : Return an index of theNode
 607 //=======================================================================
 608
 609 int SMDS_VolumeTool::GetNodeIndex(const SMDS_MeshNode* theNode) const
 610 {
 611   if ( myVolume ) {
 612     for ( int i = 0; i < myVolumeNbNodes; i++ ) {
 613       if ( myVolumeNodes[ i ] == theNode )
 614         return i;
 615     }
 616   }
 617   return -1;
 618 }
 619
 620
 621 //=======================================================================
 622 //function : IsFreeFace
 623 //purpose  : check that only one volume is build on the face nodes
 624 //=======================================================================
 625
 626 bool SMDS_VolumeTool::IsFreeFace( int faceIndex )
 627 {
 628   const int free = true;
 629   if ( !setFace( faceIndex ))
 630     return !free;
 631
 632   const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
 633   int nbFaceNodes = NbFaceNodes( faceIndex );
 634
 635   // evaluate nb of face nodes shared by other volume
 636   int maxNbShared = -1;
 637   typedef map< const SMDS_MeshElement*, int > TElemIntMap;
 638   TElemIntMap volNbShared;
 639   TElemIntMap::iterator vNbIt;
 640   for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
 641   {
 642     const SMDS_MeshNode* n = nodes[ iNode ];
 643     SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator();
 644     while ( eIt->more() ) {
 645       const SMDS_MeshElement* elem = eIt->next();
 646       if ( elem != myVolume && elem->GetType() == SMDSAbs_Volume ) {
 647         int nbShared = 1;
 648         vNbIt = volNbShared.find( elem );
 649         if ( vNbIt == volNbShared.end() )
 650           volNbShared.insert ( TElemIntMap::value_type( elem, nbShared ));
 651         else
 652           nbShared = ++(*vNbIt).second;
 653         if ( nbShared > maxNbShared )
 654           maxNbShared = nbShared;
 655       }
 656     }
 657   }
 658   if ( maxNbShared < 3 )
 659     return free; // is free
 660
 661   // find volumes laying on the opposite side of the face
 662   // and sharing all nodes
 663   XYZ intNormal; // internal normal
 664   GetFaceNormal( faceIndex, intNormal.x, intNormal.y, intNormal.z );
 665   if ( IsFaceExternal( faceIndex ))
 666     intNormal = XYZ( -intNormal.x, -intNormal.y, -intNormal.z );
 667   XYZ p0 ( nodes[0] ), baryCenter;
 668   for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); vNbIt++ )
 669   {
 670     int nbShared = (*vNbIt).second;
 671     if ( nbShared >= 3 ) {
 672       SMDS_VolumeTool volume( (*vNbIt).first );
 673       volume.GetBaryCenter( baryCenter.x, baryCenter.y, baryCenter.z );
 674       XYZ intNormal2( baryCenter - p0 );
 675       if ( intNormal.Dot( intNormal2 ) < 0 )
 676         continue; // opposite side
 677     }
 678     // remove a volume from volNbShared map
 679     volNbShared.erase( vNbIt );
 680   }
 681   // here volNbShared contains only volumes laying on the
 682   // opposite side of the face
 683   if ( volNbShared.empty() )
 684     return free; // is free
 685
 686   // check if the whole area of a face is shared
 687   bool isShared[] = { false, false, false, false }; // 4 triangle parts of a quadrangle
 688   for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); vNbIt++ )
 689   {
 690     SMDS_VolumeTool volume( (*vNbIt).first );
 691     bool prevLinkShared = false;
 692     int nbSharedLinks = 0;
 693     for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
 694     {
 695       bool linkShared = volume.IsLinked( nodes[ iNode ], nodes[ iNode + 1] );
 696       if ( linkShared )
 697         nbSharedLinks++;
 698       if ( linkShared && prevLinkShared &&
 699           volume.IsLinked( nodes[ iNode - 1 ], nodes[ iNode + 1] ))
 700         isShared[ iNode ] = true;
 701       prevLinkShared = linkShared;
 702     }
 703     if ( nbSharedLinks == nbFaceNodes )
 704       return !free; // is not free
 705     if ( nbFaceNodes == 4 ) {
 706       // check traingle parts 1 & 3
 707       if ( isShared[1] && isShared[3] )
 708         return !free; // is not free
 709       // check traingle parts 0 & 2;
 710       // 0 part could not be checked in the loop; check it here
 711       if ( isShared[2] && prevLinkShared &&
 712           volume.IsLinked( nodes[ 0 ], nodes[ 1 ] ) &&
 713           volume.IsLinked( nodes[ 1 ], nodes[ 3 ] ) )
 714         return !free; // is not free
 715     }
 716   }
 717   return free;
 718 }
 719
 720 //=======================================================================
 721 //function : GetFaceIndex
 722 //purpose  : Return index of a face formed by theFaceNodes
 723 //=======================================================================
 724
 725 int SMDS_VolumeTool::GetFaceIndex( const set<const SMDS_MeshNode*>& theFaceNodes )
 726 {
 727   for ( int iFace = 0; iFace < myNbFaces; iFace++ ) {
 728     const SMDS_MeshNode** nodes = GetFaceNodes( iFace );
 729     int nbFaceNodes = NbFaceNodes( iFace );
 730     set<const SMDS_MeshNode*> nodeSet;
 731     for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
 732       nodeSet.insert( nodes[ iNode ] );
 733     if ( theFaceNodes == nodeSet )
 734       return iFace;
 735   }
 736   return -1;
 737 }
 738
 739 //=======================================================================
 740 //function : GetFaceIndex
 741 //purpose  : Return index of a face formed by theFaceNodes
 742 //=======================================================================
 743
 744 int SMDS_VolumeTool::GetFaceIndex( const set<int>& theFaceNodesIndices )
 745 {
 746   for ( int iFace = 0; iFace < myNbFaces; iFace++ ) {
 747     const int* nodes = GetFaceNodesIndices( iFace );
 748     int nbFaceNodes = NbFaceNodes( iFace );
 749     set<int> nodeSet;
 750     for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
 751       nodeSet.insert( nodes[ iNode ] );
 752     if ( theFaceNodesIndices == nodeSet )
 753       return iFace;
 754   }
 755   return -1;
 756 }
 757
 758 //=======================================================================
 759 //function : setFace
 760 //purpose  :
 761 //=======================================================================
 762
 763 bool SMDS_VolumeTool::setFace( int faceIndex )
 764 {
 765   if ( !myVolume )
 766     return false;
 767
 768   if ( myCurFace == faceIndex )
 769     return true;
 770
 771   if ( faceIndex < 0 || faceIndex >= NbFaces() )
 772     return false;
 773
 774   // choose face node indices
 775   switch ( myVolumeNbNodes ) {
 776   case 4:
 777     if ( myExternalFaces )
 778       myFaceNodeIndices = myVolForward ? Tetra_F[ faceIndex ] : Tetra_RE[ faceIndex ];
 779     else
 780       myFaceNodeIndices = myVolForward ? Tetra_F[ faceIndex ] : Tetra_R[ faceIndex ];
 781     break;
 782   case 5:
 783     if ( myExternalFaces )
 784       myFaceNodeIndices = myVolForward ? Pyramid_F[ faceIndex ] : Pyramid_RE[ faceIndex ];
 785     else
 786       myFaceNodeIndices = myVolForward ? Pyramid_F[ faceIndex ] : Pyramid_R[ faceIndex ];
 787     break;
 788   case 6:
 789     if ( myExternalFaces )
 790       myFaceNodeIndices = myVolForward ? Penta_FE[ faceIndex ] : Penta_RE[ faceIndex ];
 791     else
 792       myFaceNodeIndices = myVolForward ? Penta_F[ faceIndex ] : Penta_R[ faceIndex ];
 793     break;
 794   case 8:
 795     if ( myExternalFaces )
 796       myFaceNodeIndices = myVolForward ? Hexa_FE[ faceIndex ] : Hexa_RE[ faceIndex ];
 797     else
 798       myFaceNodeIndices = Hexa_F[ faceIndex ];
 799     break;
 800   default: return false;
 801   }
 802
 803   // set face nodes
 804   int iNode, nbNode = myFaceNbNodes[ faceIndex ];
 805   for ( iNode = 0; iNode <= nbNode; iNode++ )
 806     myFaceNodes[ iNode ] = myVolumeNodes[ myFaceNodeIndices[ iNode ]];
 807
 808   myCurFace = faceIndex;
 809
 810   return true;
 811 }