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 "SMDS_PolyhedralVolumeOfNodes.hxx"
  15
  16 #include "utilities.h"
  17
  18 #include <map>
  19 #include <float.h>
  20 #include <math.h>
  21
  22 using namespace std;
  23
  24 // ======================================================
  25 // Node indices in faces depending on volume orientation
  26 // making most faces normals external
  27 // ======================================================
  28
  29 /*
  30 //           N3
  31 //           +
  32 //          /|\
  33 //         / | \
  34 //        /  |  \
  35 //    N0 +---|---+ N1                TETRAHEDRON
  36 //       \   |   /
  37 //        \  |  /
  38 //         \ | /
  39 //          \|/
  40 //           +
  41 //           N2
  42 */
  43 static int Tetra_F [4][4] = { // FORWARD == EXTERNAL
  44   { 0, 1, 2, 0 },              // All faces have external normals
  45   { 0, 3, 1, 0 },
  46   { 1, 3, 2, 1 },
  47   { 0, 2, 3, 0 }};
  48 static int Tetra_R [4][4] = { // REVERSED
  49   { 0, 1, 2, 0 },             // All faces but a bottom have external normals
  50   { 0, 1, 3, 0 },
  51   { 1, 2, 3, 1 },
  52   { 0, 3, 2, 0 }};
  53 static int Tetra_RE [4][4] = { // REVERSED -> FORWARD (EXTERNAL)
  54   { 0, 2, 1, 0 },              // All faces have external normals
  55   { 0, 1, 3, 0 },
  56   { 1, 2, 3, 1 },
  57   { 0, 3, 2, 0 }};
  58 static int Tetra_nbN [] = { 3, 3, 3, 3 };
  59
  60 //
  61 //     PYRAMID
  62 //
  63 static int Pyramid_F [5][5] = { // FORWARD == EXTERNAL
  64   { 0, 1, 2, 3, 0 },            // All faces have external normals
  65   { 0, 4, 1, 0, 4 },
  66   { 1, 4, 2, 1, 4 },
  67   { 2, 4, 3, 2, 4 },
  68   { 3, 4, 0, 3, 4 }};
  69 static int Pyramid_R [5][5] = { // REVERSED
  70   { 0, 1, 2, 3, 0 },            // All faces but a bottom have external normals
  71   { 0, 1, 4, 0, 4 },
  72   { 1, 2, 4, 1, 4 },
  73   { 2, 3, 4, 2, 4 },
  74   { 3, 0, 4, 3, 4 }};
  75 static int Pyramid_RE [5][5] = { // REVERSED -> FORWARD (EXTERNAL)
  76   { 0, 3, 2, 1, 0 },             // All faces but a bottom have external normals
  77   { 0, 1, 4, 0, 4 },
  78   { 1, 2, 4, 1, 4 },
  79   { 2, 3, 4, 2, 4 },
  80   { 3, 0, 4, 3, 4 }};
  81 static int Pyramid_nbN [] = { 4, 3, 3, 3, 3 };
  82
  83 /*
  84 //            + N4
  85 //           /|\
  86 //          / | \
  87 //         /  |  \
  88 //        /   |   \
  89 //    N3 +---------+ N5
  90 //       |    |    |
  91 //       |    + N1 |
  92 //       |   / \   |                PENTAHEDRON
  93 //       |  /   \  |
  94 //       | /     \ |
  95 //       |/       \|
  96 //    N0 +---------+ N2
  97 */
  98 static int Penta_F [5][5] = { // FORWARD
  99   { 0, 1, 2, 0, 0 },          // Top face has an internal normal, other - external
 100   { 3, 4, 5, 3, 3 },          // 0 is bottom, 1 is top face
 101   { 0, 2, 5, 3, 0 },
 102   { 1, 4, 5, 2, 1 },
 103   { 0, 3, 4, 1, 0 }};
 104 static int Penta_R [5][5] = { // REVERSED
 105   { 0, 1, 2, 0, 0 },          // Bottom face has an internal normal, other - external
 106   { 3, 4, 5, 3, 3 },          // 0 is bottom, 1 is top face
 107   { 0, 3, 5, 2, 0 },
 108   { 1, 2, 5, 4, 1 },
 109   { 0, 1, 4, 3, 0 }};
 110 static int Penta_FE [5][5] = { // FORWARD -> EXTERNAL
 111   { 0, 1, 2, 0, 0 },
 112   { 3, 5, 4, 3, 3 },
 113   { 0, 2, 5, 3, 0 },
 114   { 1, 4, 5, 2, 1 },
 115   { 0, 3, 4, 1, 0 }};
 116 static int Penta_RE [5][5] = { // REVERSED -> EXTERNAL
 117   { 0, 2, 1, 0, 0 },
 118   { 3, 4, 5, 3, 3 },
 119   { 0, 3, 5, 2, 0 },
 120   { 1, 2, 5, 4, 1 },
 121   { 0, 1, 4, 3, 0 }};
 122 static int Penta_nbN [] = { 3, 3, 4, 4, 4 };
 123
 124 /*
 125 //         N5+----------+N6
 126 //          /|         /|
 127 //         / |        / |
 128 //        /  |       /  |
 129 //     N4+----------+N7 |
 130 //       |   |      |   |           HEXAHEDRON
 131 //       |   |      |   |
 132 //       |   |      |   |
 133 //       | N1+------|---+N2
 134 //       |  /       |  /
 135 //       | /        | /
 136 //       |/         |/
 137 //     N0+----------+N3
 138 */
 139 static int Hexa_F [6][5] = { // FORWARD
 140   { 0, 1, 2, 3, 0 },         // opposite faces are neighbouring,
 141   { 4, 5, 6, 7, 4 },         // odd face(1,3,5) normal is internal, even(0,2,4) - external
 142   { 1, 0, 4, 5, 1 },         // same index nodes of opposite faces are linked
 143   { 2, 3, 7, 6, 2 },
 144   { 0, 3, 7, 4, 0 },
 145   { 1, 2, 6, 5, 1 }};
 146 // static int Hexa_R [6][5] = { // REVERSED
 147 //   { 0, 3, 2, 1, 0 },         // opposite faces are neighbouring,
 148 //   { 4, 7, 6, 5, 4 },         // odd face(1,3,5) normal is external, even(0,2,4) - internal
 149 //   { 1, 5, 4, 0, 1 },         // same index nodes of opposite faces are linked
 150 //   { 2, 6, 7, 3, 2 },
 151 //   { 0, 4, 7, 3, 0 },
 152 //   { 1, 5, 6, 2, 1 }};
 153 static int Hexa_FE [6][5] = { // FORWARD -> EXTERNAL
 154   { 0, 1, 2, 3, 0 } ,         // opposite faces are neighbouring,
 155   { 4, 7, 6, 5, 4 },          // all face normals are external,
 156   { 0, 4, 5, 1, 0 },          // links in opposite faces: 0-0, 1-3, 2-2, 3-1
 157   { 3, 2, 6, 7, 3 },
 158   { 0, 3, 7, 4, 0 },
 159   { 1, 5, 6, 2, 1 }};
 160 static int Hexa_RE [6][5] = { // REVERSED -> EXTERNAL
 161   { 0, 3, 2, 1, 0 },          // opposite faces are neighbouring,
 162   { 4, 5, 6, 7, 4 },          // all face normals are external,
 163   { 0, 1, 5, 4, 0 },          // links in opposite faces: 0-0, 1-3, 2-2, 3-1
 164   { 3, 7, 6, 2, 3 },
 165   { 0, 4, 7, 3, 0 },
 166   { 1, 2, 6, 5, 1 }};
 167 static int Hexa_nbN [] = { 4, 4, 4, 4, 4, 4 };
 168
 169 // ========================================================
 170 // to perform some calculations without linkage to CASCADE
 171 // ========================================================
 172 struct XYZ {
 173   double x;
 174   double y;
 175   double z;
 176   XYZ()                               { x = 0; y = 0; z = 0; }
 177   XYZ( double X, double Y, double Z ) { x = X; y = Y; z = Z; }
 178   XYZ( const XYZ& other )             { x = other.x; y = other.y; z = other.z; }
 179   XYZ( const SMDS_MeshNode* n )       { x = n->X(); y = n->Y(); z = n->Z(); }
 180   XYZ operator-( const XYZ& other );
 181   XYZ Crossed( const XYZ& other );
 182   double Dot( const XYZ& other );
 183   double Magnitude();
 184 };
 185 XYZ XYZ::operator-( const XYZ& Right ) {
 186   return XYZ(x - Right.x, y - Right.y, z - Right.z);
 187 }
 188 XYZ XYZ::Crossed( const XYZ& Right ) {
 189   return XYZ (y * Right.z - z * Right.y,
 190               z * Right.x - x * Right.z,
 191               x * Right.y - y * Right.x);
 192 }
 193 double XYZ::Dot( const XYZ& Other ) {
 194   return(x * Other.x + y * Other.y + z * Other.z);
 195 }
 196 double XYZ::Magnitude() {
 197   return sqrt (x * x + y * y + z * z);
 198 }
 199
 200 //=======================================================================
 201 //function : SMDS_VolumeTool
 202 //purpose  :
 203 //=======================================================================
 204
 205 SMDS_VolumeTool::SMDS_VolumeTool ()
 206      : myVolume( 0 ),
 207        myPolyedre( 0 ),
 208        myVolForward( true ),
 209        myNbFaces( 0 ),
 210        myVolumeNbNodes( 0 ),
 211        myVolumeNodes( NULL ),
 212        myExternalFaces( false ),
 213        myFaceNbNodes( 0 ),
 214        myCurFace( -1 ),
 215        myFaceNodeIndices( NULL ),
 216        myFaceNodes( NULL )
 217 {
 218 }
 219
 220 //=======================================================================
 221 //function : SMDS_VolumeTool
 222 //purpose  :
 223 //=======================================================================
 224
 225 SMDS_VolumeTool::SMDS_VolumeTool (const SMDS_MeshElement* theVolume)
 226      : myVolume( 0 ),
 227        myPolyedre( 0 ),
 228        myVolForward( true ),
 229        myNbFaces( 0 ),
 230        myVolumeNbNodes( 0 ),
 231        myVolumeNodes( NULL ),
 232        myExternalFaces( false ),
 233        myFaceNbNodes( 0 ),
 234        myCurFace( -1 ),
 235        myFaceNodeIndices( NULL ),
 236        myFaceNodes( NULL )
 237 {
 238   Set( theVolume );
 239 }
 240
 241 //=======================================================================
 242 //function : SMDS_VolumeTool
 243 //purpose  :
 244 //=======================================================================
 245
 246 SMDS_VolumeTool::~SMDS_VolumeTool()
 247 {
 248   if (myVolumeNodes != NULL) {
 249     delete [] myVolumeNodes;
 250     myVolumeNodes = NULL;
 251   }
 252   if (myFaceNodes != NULL) {
 253     delete [] myFaceNodes;
 254     myFaceNodes = NULL;
 255   }
 256 }
 257
 258 //=======================================================================
 259 //function : SetVolume
 260 //purpose  : Set volume to iterate on
 261 //=======================================================================
 262
 263 bool SMDS_VolumeTool::Set (const SMDS_MeshElement* theVolume)
 264 {
 265   myVolume = 0;
 266   myPolyedre = 0;
 267
 268   myVolForward = true;
 269   myNbFaces = 0;
 270   myVolumeNbNodes = 0;
 271   if (myVolumeNodes != NULL) {
 272     delete [] myVolumeNodes;
 273     myVolumeNodes = NULL;
 274   }
 275
 276   myExternalFaces = false;
 277   myFaceNbNodes = 0;
 278
 279   myCurFace = -1;
 280   myFaceNodeIndices = NULL;
 281   if (myFaceNodes != NULL) {
 282     delete [] myFaceNodes;
 283     myFaceNodes = NULL;
 284   }
 285
 286   if ( theVolume && theVolume->GetType() == SMDSAbs_Volume )
 287   {
 288     myVolume = theVolume;
 289
 290     myNbFaces = theVolume->NbFaces();
 291     myVolumeNbNodes = theVolume->NbNodes();
 292
 293     // set volume nodes
 294     int iNode = 0;
 295     myVolumeNodes = new const SMDS_MeshNode* [myVolumeNbNodes];
 296     SMDS_ElemIteratorPtr nodeIt = myVolume->nodesIterator();
 297     while ( nodeIt->more() ) {
 298       myVolumeNodes[ iNode++ ] = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
 299     }
 300
 301     if (myVolume->IsPoly()) {
 302       myPolyedre = static_cast<const SMDS_PolyhedralVolumeOfNodes*>( myVolume );
 303       if (!myPolyedre) {
 304         MESSAGE("Warning: bad volumic element");
 305         return false;
 306       }
 307     } else {
 308       switch ( myVolumeNbNodes ) {
 309       case 4:
 310       case 5:
 311       case 6:
 312       case 8: {
 313         // define volume orientation
 314         XYZ botNormal;
 315         GetFaceNormal( 0, botNormal.x, botNormal.y, botNormal.z );
 316         const SMDS_MeshNode* topNode = myVolumeNodes[ myVolumeNbNodes - 1 ];
 317         const SMDS_MeshNode* botNode = myVolumeNodes[ 0 ];
 318         XYZ upDir (topNode->X() - botNode->X(),
 319                    topNode->Y() - botNode->Y(),
 320                    topNode->Z() - botNode->Z() );
 321         myVolForward = ( botNormal.Dot( upDir ) < 0 );
 322         break;
 323       }
 324       default:
 325         break;
 326       }
 327     }
 328   }
 329   return ( myVolume != 0 );
 330 }
 331
 332 //=======================================================================
 333 //function : Inverse
 334 //purpose  : Inverse volume
 335 //=======================================================================
 336
 337 #define SWAP_NODES(nodes,i1,i2)           \
 338 {                                         \
 339   const SMDS_MeshNode* tmp = nodes[ i1 ]; \
 340   nodes[ i1 ] = nodes[ i2 ];              \
 341   nodes[ i2 ] = tmp;                      \
 342 }
 343 void SMDS_VolumeTool::Inverse ()
 344 {
 345   if ( !myVolume ) return;
 346
 347   if (myVolume->IsPoly()) {
 348     MESSAGE("Warning: attempt to inverse polyhedral volume");
 349     return;
 350   }
 351
 352   myVolForward = !myVolForward;
 353   myCurFace = -1;
 354
 355   // inverse top and bottom faces
 356   switch ( myVolumeNbNodes ) {
 357   case 4:
 358     SWAP_NODES( myVolumeNodes, 1, 2 );
 359     break;
 360   case 5:
 361     SWAP_NODES( myVolumeNodes, 1, 3 );
 362     break;
 363   case 6:
 364     SWAP_NODES( myVolumeNodes, 1, 2 );
 365     SWAP_NODES( myVolumeNodes, 4, 5 );
 366     break;
 367   case 8:
 368     SWAP_NODES( myVolumeNodes, 1, 3 );
 369     SWAP_NODES( myVolumeNodes, 5, 7 );
 370     break;
 371   default:;
 372   }
 373 }
 374
 375 //=======================================================================
 376 //function : GetSize
 377 //purpose  : Return element volume
 378 //=======================================================================
 379
 380 double SMDS_VolumeTool::GetSize() const
 381 {
 382   return 0;
 383 }
 384
 385 //=======================================================================
 386 //function : GetBaryCenter
 387 //purpose  :
 388 //=======================================================================
 389
 390 bool SMDS_VolumeTool::GetBaryCenter(double & X, double & Y, double & Z) const
 391 {
 392   X = Y = Z = 0.;
 393   if ( !myVolume )
 394     return false;
 395
 396   for ( int i = 0; i < myVolumeNbNodes; i++ ) {
 397     X += myVolumeNodes[ i ]->X();
 398     Y += myVolumeNodes[ i ]->Y();
 399     Z += myVolumeNodes[ i ]->Z();
 400   }
 401   X /= myVolumeNbNodes;
 402   Y /= myVolumeNbNodes;
 403   Z /= myVolumeNbNodes;
 404
 405   return true;
 406 }
 407
 408 //=======================================================================
 409 //function : SetExternalNormal
 410 //purpose  : Node order will be so that faces normals are external
 411 //=======================================================================
 412
 413 void SMDS_VolumeTool::SetExternalNormal ()
 414 {
 415   myExternalFaces = true;
 416   myCurFace = -1;
 417 }
 418
 419 //=======================================================================
 420 //function : NbFaceNodes
 421 //purpose  : Return number of nodes in the array of face nodes
 422 //=======================================================================
 423
 424 int SMDS_VolumeTool::NbFaceNodes( int faceIndex )
 425 {
 426     if ( !setFace( faceIndex ))
 427       return 0;
 428     return myFaceNbNodes;
 429 }
 430
 431 //=======================================================================
 432 //function : GetFaceNodes
 433 //purpose  : Return pointer to the array of face nodes.
 434 //           To comfort link iteration, the array
 435 //           length == NbFaceNodes( faceIndex ) + 1 and
 436 //           the last node == the first one.
 437 //=======================================================================
 438
 439 const SMDS_MeshNode** SMDS_VolumeTool::GetFaceNodes( int faceIndex )
 440 {
 441   if ( !setFace( faceIndex ))
 442     return 0;
 443   return myFaceNodes;
 444 }
 445
 446 //=======================================================================
 447 //function : GetFaceNodesIndices
 448 //purpose  : Return pointer to the array of face nodes indices
 449 //           To comfort link iteration, the array
 450 //           length == NbFaceNodes( faceIndex ) + 1 and
 451 //           the last node index == the first one.
 452 //=======================================================================
 453
 454 const int* SMDS_VolumeTool::GetFaceNodesIndices( int faceIndex )
 455 {
 456   if (myVolume->IsPoly()) {
 457     MESSAGE("Warning: attempt to obtain FaceNodesIndices of polyhedral volume");
 458     return NULL;
 459   }
 460   if ( !setFace( faceIndex ))
 461     return 0;
 462   return myFaceNodeIndices;
 463 }
 464
 465 //=======================================================================
 466 //function : GetFaceNodes
 467 //purpose  : Return a set of face nodes.
 468 //=======================================================================
 469
 470 bool SMDS_VolumeTool::GetFaceNodes (int                        faceIndex,
 471                                     set<const SMDS_MeshNode*>& theFaceNodes )
 472 {
 473   if ( !setFace( faceIndex ))
 474     return false;
 475
 476   theFaceNodes.clear();
 477   int iNode, nbNode = myFaceNbNodes;
 478   for ( iNode = 0; iNode < nbNode; iNode++ )
 479     theFaceNodes.insert( myFaceNodes[ iNode ]);
 480
 481   return true;
 482 }
 483
 484 //=======================================================================
 485 //function : IsFaceExternal
 486 //purpose  : Check normal orientation of a returned face
 487 //=======================================================================
 488
 489 bool SMDS_VolumeTool::IsFaceExternal( int faceIndex )
 490 {
 491   if ( myExternalFaces || !myVolume )
 492     return true;
 493
 494   if (myVolume->IsPoly()) {
 495     XYZ aNormal, baryCenter, p0 (myPolyedre->GetFaceNode(faceIndex + 1, 1));
 496     GetFaceNormal(faceIndex, aNormal.x, aNormal.y, aNormal.z);
 497     GetBaryCenter(baryCenter.x, baryCenter.y, baryCenter.z);
 498     XYZ insideVec (baryCenter - p0);
 499     if (insideVec.Dot(aNormal) > 0)
 500       return false;
 501     return true;
 502   }
 503
 504   switch ( myVolumeNbNodes ) {
 505   case 4:
 506   case 5:
 507     // only the bottom of a reversed tetrahedron can be internal
 508     return ( myVolForward || faceIndex != 0 );
 509   case 6:
 510     // in a forward pentahedron, the top is internal, in a reversed one - bottom
 511     return ( myVolForward ? faceIndex != 1 : faceIndex != 0 );
 512   case 8: {
 513     // in a forward hexahedron, even face normal is external, odd - internal
 514     bool odd = faceIndex % 2;
 515     return ( myVolForward ? !odd : odd );
 516   }
 517   default:;
 518   }
 519   return false;
 520 }
 521
 522 //=======================================================================
 523 //function : GetFaceNormal
 524 //purpose  : Return a normal to a face
 525 //=======================================================================
 526
 527 bool SMDS_VolumeTool::GetFaceNormal (int faceIndex, double & X, double & Y, double & Z)
 528 {
 529   if ( !setFace( faceIndex ))
 530     return false;
 531
 532   XYZ p1 ( myFaceNodes[0] );
 533   XYZ p2 ( myFaceNodes[1] );
 534   XYZ p3 ( myFaceNodes[2] );
 535   XYZ aVec12( p2 - p1 );
 536   XYZ aVec13( p3 - p1 );
 537   XYZ cross = aVec12.Crossed( aVec13 );
 538
 539   if ( myFaceNbNodes[ faceIndex ] == 4 ) {
 540     XYZ p4 ( myFaceNodes[3] );
 541     XYZ aVec14( p4 - p1 );
 542     XYZ cross2 = aVec13.Crossed( aVec14 );
 543     cross.x += cross2.x;
 544     cross.y += cross2.y;
 545     cross.z += cross2.z;
 546   }
 547
 548   double size = cross.Magnitude();
 549   if ( size <= DBL_MIN )
 550     return false;
 551
 552   X = cross.x / size;
 553   Y = cross.y / size;
 554   Z = cross.z / size;
 555
 556   return true;
 557 }
 558
 559 //=======================================================================
 560 //function : GetFaceArea
 561 //purpose  : Return face area
 562 //=======================================================================
 563
 564 double SMDS_VolumeTool::GetFaceArea( int faceIndex )
 565 {
 566   if (myVolume->IsPoly()) {
 567     MESSAGE("Warning: attempt to obtain area of a face of polyhedral volume");
 568     return 0;
 569   }
 570
 571   if ( !setFace( faceIndex ))
 572     return 0;
 573
 574   XYZ p1 ( myFaceNodes[0] );
 575   XYZ p2 ( myFaceNodes[1] );
 576   XYZ p3 ( myFaceNodes[2] );
 577   XYZ aVec12( p2 - p1 );
 578   XYZ aVec13( p3 - p1 );
 579   double area = aVec12.Crossed( aVec13 ).Magnitude() * 0.5;
 580
 581   if ( myFaceNbNodes == 4 ) {
 582     XYZ p4 ( myFaceNodes[3] );
 583     XYZ aVec14( p4 - p1 );
 584     area += aVec14.Crossed( aVec13 ).Magnitude() * 0.5;
 585   }
 586   return area;
 587 }
 588
 589 //=======================================================================
 590 //function : GetOppFaceIndex
 591 //purpose  : Return index of the opposite face if it exists, else -1.
 592 //=======================================================================
 593
 594 int SMDS_VolumeTool::GetOppFaceIndex( int faceIndex ) const
 595 {
 596   int ind = -1;
 597   if (myVolume->IsPoly()) {
 598     MESSAGE("Warning: attempt to obtain opposite face on polyhedral volume");
 599     return ind;
 600   }
 601
 602   if ( faceIndex >= 0 && faceIndex < NbFaces() ) {
 603     switch ( myVolumeNbNodes ) {
 604     case 6:
 605       if ( faceIndex == 0 || faceIndex == 1 )
 606         ind = 1 - faceIndex;
 607         break;
 608     case 8:
 609       ind = faceIndex + ( faceIndex % 2 ? -1 : 1 );
 610       break;
 611     default:;
 612     }
 613   }
 614   return ind;
 615 }
 616
 617 //=======================================================================
 618 //function : IsLinked
 619 //purpose  : return true if theNode1 is linked with theNode2
 620 //=======================================================================
 621
 622 bool SMDS_VolumeTool::IsLinked (const SMDS_MeshNode* theNode1,
 623                                 const SMDS_MeshNode* theNode2) const
 624 {
 625   if ( !myVolume )
 626     return false;
 627
 628   if (myVolume->IsPoly()) {
 629     if (!myPolyedre) {
 630       MESSAGE("Warning: bad volumic element");
 631       return false;
 632     }
 633     bool isLinked = false;
 634     int iface;
 635     for (iface = 1; iface <= myNbFaces && !isLinked; iface++) {
 636       int inode, nbFaceNodes = myPolyedre->NbFaceNodes(iface);
 637
 638       for (inode = 1; inode <= nbFaceNodes && !isLinked; inode++) {
 639         const SMDS_MeshNode* curNode = myPolyedre->GetFaceNode(iface, inode);
 640
 641         if (curNode == theNode1 || curNode == theNode2) {
 642           int inextnode = (inode == nbFaceNodes) ? 1 : inode + 1;
 643           const SMDS_MeshNode* nextNode = myPolyedre->GetFaceNode(iface, inextnode);
 644
 645           if ((curNode == theNode1 && nextNode == theNode2) ||
 646               (curNode == theNode2 && nextNode == theNode1)) {
 647             isLinked = true;
 648           }
 649         }
 650       }
 651     }
 652     return isLinked;
 653   }
 654
 655   // find nodes indices
 656   int i1 = -1, i2 = -1;
 657   for ( int i = 0; i < myVolumeNbNodes; i++ ) {
 658     if ( myVolumeNodes[ i ] == theNode1 )
 659       i1 = i;
 660     else if ( myVolumeNodes[ i ] == theNode2 )
 661       i2 = i;
 662   }
 663   return IsLinked( i1, i2 );
 664 }
 665
 666 //=======================================================================
 667 //function : IsLinked
 668 //purpose  : return true if the node with theNode1Index is linked
 669 //           with the node with theNode2Index
 670 //=======================================================================
 671
 672 bool SMDS_VolumeTool::IsLinked (const int theNode1Index,
 673                                 const int theNode2Index) const
 674 {
 675   if (myVolume->IsPoly()) {
 676     return IsLinked(myVolumeNodes[theNode1Index], myVolumeNodes[theNode2Index]);
 677   }
 678
 679   int minInd = theNode1Index < theNode2Index ? theNode1Index : theNode2Index;
 680   int maxInd = theNode1Index < theNode2Index ? theNode2Index : theNode1Index;
 681
 682   if ( minInd < 0 || maxInd > myVolumeNbNodes - 1 || maxInd == minInd )
 683     return false;
 684
 685   switch ( myVolumeNbNodes ) {
 686   case 4:
 687     return true;
 688   case 5:
 689     if ( maxInd == 4 )
 690       return true;
 691     switch ( maxInd - minInd ) {
 692     case 1:
 693     case 3: return true;
 694     default:;
 695     }
 696     break;
 697   case 6:
 698     switch ( maxInd - minInd ) {
 699     case 1: return minInd != 2;
 700     case 2: return minInd == 0 || minInd == 3;
 701     case 3: return true;
 702     default:;
 703     }
 704     break;
 705   case 8:
 706     switch ( maxInd - minInd ) {
 707     case 1: return minInd != 3;
 708     case 3: return minInd == 0 || minInd == 4;
 709     case 4: return true;
 710     default:;
 711     }
 712     break;
 713   default:;
 714   }
 715   return false;
 716 }
 717
 718 //=======================================================================
 719 //function : GetNodeIndex
 720 //purpose  : Return an index of theNode
 721 //=======================================================================
 722
 723 int SMDS_VolumeTool::GetNodeIndex(const SMDS_MeshNode* theNode) const
 724 {
 725   if ( myVolume ) {
 726     for ( int i = 0; i < myVolumeNbNodes; i++ ) {
 727       if ( myVolumeNodes[ i ] == theNode )
 728         return i;
 729     }
 730   }
 731   return -1;
 732 }
 733
 734 //=======================================================================
 735 //function : IsFreeFace
 736 //purpose  : check that only one volume is build on the face nodes
 737 //=======================================================================
 738
 739 bool SMDS_VolumeTool::IsFreeFace( int faceIndex )
 740 {
 741   const int free = true;
 742
 743   if (!setFace( faceIndex ))
 744     return !free;
 745
 746   const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
 747   int nbFaceNodes = myFaceNbNodes;
 748
 749   // evaluate nb of face nodes shared by other volume
 750   int maxNbShared = -1;
 751   typedef map< const SMDS_MeshElement*, int > TElemIntMap;
 752   TElemIntMap volNbShared;
 753   TElemIntMap::iterator vNbIt;
 754   for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
 755   {
 756     const SMDS_MeshNode* n = nodes[ iNode ];
 757     SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator();
 758     while ( eIt->more() ) {
 759       const SMDS_MeshElement* elem = eIt->next();
 760       if ( elem != myVolume && elem->GetType() == SMDSAbs_Volume ) {
 761         int nbShared = 1;
 762         vNbIt = volNbShared.find( elem );
 763         if ( vNbIt == volNbShared.end() )
 764           volNbShared.insert ( TElemIntMap::value_type( elem, nbShared ));
 765         else
 766           nbShared = ++(*vNbIt).second;
 767         if ( nbShared > maxNbShared )
 768           maxNbShared = nbShared;
 769       }
 770     }
 771   }
 772   if ( maxNbShared < 3 )
 773     return free; // is free
 774
 775   // find volumes laying on the opposite side of the face
 776   // and sharing all nodes
 777   XYZ intNormal; // internal normal
 778   GetFaceNormal( faceIndex, intNormal.x, intNormal.y, intNormal.z );
 779   if ( IsFaceExternal( faceIndex ))
 780     intNormal = XYZ( -intNormal.x, -intNormal.y, -intNormal.z );
 781   XYZ p0 ( nodes[0] ), baryCenter;
 782   for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); vNbIt++ )
 783   {
 784     int nbShared = (*vNbIt).second;
 785     if ( nbShared >= 3 ) {
 786       SMDS_VolumeTool volume( (*vNbIt).first );
 787       volume.GetBaryCenter( baryCenter.x, baryCenter.y, baryCenter.z );
 788       XYZ intNormal2( baryCenter - p0 );
 789       if ( intNormal.Dot( intNormal2 ) < 0 )
 790         continue; // opposite side
 791     }
 792     // remove a volume from volNbShared map
 793     volNbShared.erase( vNbIt );
 794   }
 795   // here volNbShared contains only volumes laying on the
 796   // opposite side of the face
 797   if ( volNbShared.empty() )
 798     return free; // is free
 799
 800   // check if the whole area of a face is shared
 801   bool isShared[] = { false, false, false, false }; // 4 triangle parts of a quadrangle
 802   for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); vNbIt++ )
 803   {
 804     SMDS_VolumeTool volume( (*vNbIt).first );
 805     bool prevLinkShared = false;
 806     int nbSharedLinks = 0;
 807     for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
 808     {
 809       bool linkShared = volume.IsLinked( nodes[ iNode ], nodes[ iNode + 1] );
 810       if ( linkShared )
 811         nbSharedLinks++;
 812       if ( linkShared && prevLinkShared &&
 813           volume.IsLinked( nodes[ iNode - 1 ], nodes[ iNode + 1] ))
 814         isShared[ iNode ] = true;
 815       prevLinkShared = linkShared;
 816     }
 817     if ( nbSharedLinks == nbFaceNodes )
 818       return !free; // is not free
 819     if ( nbFaceNodes == 4 ) {
 820       // check traingle parts 1 & 3
 821       if ( isShared[1] && isShared[3] )
 822         return !free; // is not free
 823       // check triangle parts 0 & 2;
 824       // 0 part could not be checked in the loop; check it here
 825       if ( isShared[2] && prevLinkShared &&
 826           volume.IsLinked( nodes[ 0 ], nodes[ 1 ] ) &&
 827           volume.IsLinked( nodes[ 1 ], nodes[ 3 ] ) )
 828         return !free; // is not free
 829     }
 830   }
 831   return free;
 832 }
 833
 834 //=======================================================================
 835 //function : GetFaceIndex
 836 //purpose  : Return index of a face formed by theFaceNodes
 837 //=======================================================================
 838
 839 int SMDS_VolumeTool::GetFaceIndex( const set<const SMDS_MeshNode*>& theFaceNodes )
 840 {
 841   for ( int iFace = 0; iFace < myNbFaces; iFace++ ) {
 842     const SMDS_MeshNode** nodes = GetFaceNodes( iFace );
 843     int nbFaceNodes = NbFaceNodes( iFace );
 844     set<const SMDS_MeshNode*> nodeSet;
 845     for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
 846       nodeSet.insert( nodes[ iNode ] );
 847     if ( theFaceNodes == nodeSet )
 848       return iFace;
 849   }
 850   return -1;
 851 }
 852
 853 //=======================================================================
 854 //function : GetFaceIndex
 855 //purpose  : Return index of a face formed by theFaceNodes
 856 //=======================================================================
 857
 858 /*int SMDS_VolumeTool::GetFaceIndex( const set<int>& theFaceNodesIndices )
 859 {
 860   for ( int iFace = 0; iFace < myNbFaces; iFace++ ) {
 861     const int* nodes = GetFaceNodesIndices( iFace );
 862     int nbFaceNodes = NbFaceNodes( iFace );
 863     set<int> nodeSet;
 864     for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
 865       nodeSet.insert( nodes[ iNode ] );
 866     if ( theFaceNodesIndices == nodeSet )
 867       return iFace;
 868   }
 869   return -1;
 870 }*/
 871
 872 //=======================================================================
 873 //function : setFace
 874 //purpose  :
 875 //=======================================================================
 876
 877 bool SMDS_VolumeTool::setFace( int faceIndex )
 878 {
 879   if ( !myVolume )
 880     return false;
 881
 882   if ( myCurFace == faceIndex )
 883     return true;
 884
 885   myCurFace = -1;
 886
 887   if ( faceIndex < 0 || faceIndex >= NbFaces() )
 888     return false;
 889
 890   if (myFaceNodes != NULL) {
 891     delete [] myFaceNodes;
 892     myFaceNodes = NULL;
 893   }
 894
 895   if (myVolume->IsPoly()) {
 896     if (!myPolyedre) {
 897       MESSAGE("Warning: bad volumic element");
 898       return false;
 899     }
 900
 901     // check orientation
 902     bool isGoodOri = true;
 903     if (myExternalFaces) {
 904       // get natural orientation
 905       XYZ aNormal, baryCenter, p0 (myPolyedre->GetFaceNode(faceIndex + 1, 1));
 906       SMDS_VolumeTool vTool (myPolyedre);
 907       vTool.GetFaceNormal(faceIndex, aNormal.x, aNormal.y, aNormal.z);
 908       vTool.GetBaryCenter(baryCenter.x, baryCenter.y, baryCenter.z);
 909       XYZ insideVec (baryCenter - p0);
 910       if (insideVec.Dot(aNormal) > 0)
 911         isGoodOri = false;
 912     }
 913
 914     // set face nodes
 915     int iNode;
 916     myFaceNbNodes = myPolyedre->NbFaceNodes(faceIndex + 1);
 917     myFaceNodes = new const SMDS_MeshNode* [myFaceNbNodes + 1];
 918     if (isGoodOri) {
 919       for ( iNode = 0; iNode < myFaceNbNodes; iNode++ )
 920         myFaceNodes[ iNode ] = myPolyedre->GetFaceNode(faceIndex + 1, iNode + 1);
 921     } else {
 922       for ( iNode = 0; iNode < myFaceNbNodes; iNode++ )
 923         myFaceNodes[ iNode ] = myPolyedre->GetFaceNode(faceIndex + 1, myFaceNbNodes - iNode);
 924     }
 925     myFaceNodes[ myFaceNbNodes ] = myFaceNodes[ 0 ]; // last = first
 926
 927   } else {
 928     // choose face node indices
 929     switch ( myVolumeNbNodes ) {
 930     case 4:
 931       myFaceNbNodes = Tetra_nbN[ faceIndex ];
 932       if ( myExternalFaces )
 933         myFaceNodeIndices = myVolForward ? Tetra_F[ faceIndex ] : Tetra_RE[ faceIndex ];
 934       else
 935         myFaceNodeIndices = myVolForward ? Tetra_F[ faceIndex ] : Tetra_R[ faceIndex ];
 936       break;
 937     case 5:
 938       myFaceNbNodes = Pyramid_nbN[ faceIndex ];
 939       if ( myExternalFaces )
 940         myFaceNodeIndices = myVolForward ? Pyramid_F[ faceIndex ] : Pyramid_RE[ faceIndex ];
 941       else
 942         myFaceNodeIndices = myVolForward ? Pyramid_F[ faceIndex ] : Pyramid_R[ faceIndex ];
 943       break;
 944     case 6:
 945       myFaceNbNodes = Penta_nbN[ faceIndex ];
 946       if ( myExternalFaces )
 947         myFaceNodeIndices = myVolForward ? Penta_FE[ faceIndex ] : Penta_RE[ faceIndex ];
 948       else
 949         myFaceNodeIndices = myVolForward ? Penta_F[ faceIndex ] : Penta_R[ faceIndex ];
 950       break;
 951     case 8:
 952       myFaceNbNodes = Hexa_nbN[ faceIndex ];
 953       if ( myExternalFaces )
 954         myFaceNodeIndices = myVolForward ? Hexa_FE[ faceIndex ] : Hexa_RE[ faceIndex ];
 955       else
 956         myFaceNodeIndices = Hexa_F[ faceIndex ];
 957       break;
 958     default:
 959       return false;
 960     }
 961
 962     // set face nodes
 963     myFaceNodes = new const SMDS_MeshNode* [myFaceNbNodes + 1];
 964     for ( int iNode = 0; iNode < myFaceNbNodes; iNode++ )
 965       myFaceNodes[ iNode ] = myVolumeNodes[ myFaceNodeIndices[ iNode ]];
 966     myFaceNodes[ myFaceNbNodes ] = myFaceNodes[ 0 ];
 967   }
 968
 969   myCurFace = faceIndex;
 970
 971   return true;
 972 }
 973
 974 //=======================================================================
 975 //function : GetType
 976 //purpose  : return VolumeType by nb of nodes in a volume
 977 //=======================================================================
 978
 979 SMDS_VolumeTool::VolumeType SMDS_VolumeTool::GetType(int nbNodes)
 980 {
 981   switch ( nbNodes ) {
 982   case 4: return TETRA;
 983   case 5: return PYRAM;
 984   case 6: return PENTA;
 985   case 8: return HEXA;
 986   default:return UNKNOWN;
 987   }
 988 }
 989
 990 //=======================================================================
 991 //function : NbFaces
 992 //purpose  : return nb of faces by volume type
 993 //=======================================================================
 994
 995 int SMDS_VolumeTool::NbFaces( VolumeType type )
 996 {
 997   switch ( type ) {
 998   case TETRA: return 4;
 999   case PYRAM: return 5;
1000   case PENTA: return 5;
1001   case HEXA : return 6;
1002   default:    return 0;
1003   }
1004 }
1005
1006 //=======================================================================
1007 //function : GetFaceNodesIndices
1008 //purpose  : Return the array of face nodes indices
1009 //           To comfort link iteration, the array
1010 //           length == NbFaceNodes( faceIndex ) + 1 and
1011 //           the last node index == the first one.
1012 //=======================================================================
1013
1014 const int* SMDS_VolumeTool::GetFaceNodesIndices(VolumeType type,
1015                                                 int        faceIndex,
1016                                                 bool       external)
1017 {
1018   switch ( type ) {
1019   case TETRA: return Tetra_F[ faceIndex ];
1020   case PYRAM: return Pyramid_F[ faceIndex ];
1021   case PENTA: return external ? Penta_FE[ faceIndex ] : Penta_F[ faceIndex ];
1022   case HEXA:  return external ? Hexa_FE[ faceIndex ] : Hexa_F[ faceIndex ];
1023   default:;
1024   }
1025   return 0;
1026 }
1027
1028 //=======================================================================
1029 //function : NbFaceNodes
1030 //purpose  : Return number of nodes in the array of face nodes
1031 //=======================================================================
1032
1033 int SMDS_VolumeTool::NbFaceNodes(VolumeType type,
1034                                  int        faceIndex )
1035 {
1036   switch ( type ) {
1037   case TETRA: return Tetra_nbN[ faceIndex ];
1038   case PYRAM: return Pyramid_nbN[ faceIndex ];
1039   case PENTA: return Penta_nbN[ faceIndex ];
1040   case HEXA:  return Hexa_nbN[ faceIndex ];
1041   default:;
1042   }
1043   return 0;
1044 }
1045