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   double size = cross.Magnitude();
 540   if ( size <= DBL_MIN )
 541     return false;
 542
 543   X = cross.x / size;
 544   Y = cross.y / size;
 545   Z = cross.z / size;
 546
 547   return true;
 548 }
 549
 550 //=======================================================================
 551 //function : GetFaceArea
 552 //purpose  : Return face area
 553 //=======================================================================
 554
 555 double SMDS_VolumeTool::GetFaceArea( int faceIndex )
 556 {
 557   if (myVolume->IsPoly()) {
 558     MESSAGE("Warning: attempt to obtain area of a face of polyhedral volume");
 559     return 0;
 560   }
 561
 562   if ( !setFace( faceIndex ))
 563     return 0;
 564
 565   XYZ p1 ( myFaceNodes[0] );
 566   XYZ p2 ( myFaceNodes[1] );
 567   XYZ p3 ( myFaceNodes[2] );
 568   XYZ aVec12( p2 - p1 );
 569   XYZ aVec13( p3 - p1 );
 570   double area = aVec12.Crossed( aVec13 ).Magnitude() * 0.5;
 571
 572   if ( myFaceNbNodes == 4 ) {
 573     XYZ p4 ( myFaceNodes[3] );
 574     XYZ aVec14( p4 - p1 );
 575     area += aVec14.Crossed( aVec13 ).Magnitude() * 0.5;
 576   }
 577   return area;
 578 }
 579
 580 //=======================================================================
 581 //function : GetOppFaceIndex
 582 //purpose  : Return index of the opposite face if it exists, else -1.
 583 //=======================================================================
 584
 585 int SMDS_VolumeTool::GetOppFaceIndex( int faceIndex ) const
 586 {
 587   int ind = -1;
 588   if (myVolume->IsPoly()) {
 589     MESSAGE("Warning: attempt to obtain opposite face on polyhedral volume");
 590     return ind;
 591   }
 592
 593   if ( faceIndex >= 0 && faceIndex < NbFaces() ) {
 594     switch ( myVolumeNbNodes ) {
 595     case 6:
 596       if ( faceIndex == 0 || faceIndex == 1 )
 597         ind = 1 - faceIndex;
 598         break;
 599     case 8:
 600       ind = faceIndex + ( faceIndex % 2 ? -1 : 1 );
 601       break;
 602     default:;
 603     }
 604   }
 605   return ind;
 606 }
 607
 608 //=======================================================================
 609 //function : IsLinked
 610 //purpose  : return true if theNode1 is linked with theNode2
 611 //=======================================================================
 612
 613 bool SMDS_VolumeTool::IsLinked (const SMDS_MeshNode* theNode1,
 614                                 const SMDS_MeshNode* theNode2) const
 615 {
 616   if ( !myVolume )
 617     return false;
 618
 619   if (myVolume->IsPoly()) {
 620     if (!myPolyedre) {
 621       MESSAGE("Warning: bad volumic element");
 622       return false;
 623     }
 624     bool isLinked = false;
 625     int iface;
 626     for (iface = 1; iface <= myNbFaces && !isLinked; iface++) {
 627       int inode, nbFaceNodes = myPolyedre->NbFaceNodes(iface);
 628
 629       for (inode = 1; inode <= nbFaceNodes && !isLinked; inode++) {
 630         const SMDS_MeshNode* curNode = myPolyedre->GetFaceNode(iface, inode);
 631
 632         if (curNode == theNode1 || curNode == theNode2) {
 633           int inextnode = (inode == nbFaceNodes) ? 1 : inode + 1;
 634           const SMDS_MeshNode* nextNode = myPolyedre->GetFaceNode(iface, inextnode);
 635
 636           if ((curNode == theNode1 && nextNode == theNode2) ||
 637               (curNode == theNode2 && nextNode == theNode1)) {
 638             isLinked = true;
 639           }
 640         }
 641       }
 642     }
 643     return isLinked;
 644   }
 645
 646   // find nodes indices
 647   int i1 = -1, i2 = -1;
 648   for ( int i = 0; i < myVolumeNbNodes; i++ ) {
 649     if ( myVolumeNodes[ i ] == theNode1 )
 650       i1 = i;
 651     else if ( myVolumeNodes[ i ] == theNode2 )
 652       i2 = i;
 653   }
 654   return IsLinked( i1, i2 );
 655 }
 656
 657 //=======================================================================
 658 //function : IsLinked
 659 //purpose  : return true if the node with theNode1Index is linked
 660 //           with the node with theNode2Index
 661 //=======================================================================
 662
 663 bool SMDS_VolumeTool::IsLinked (const int theNode1Index,
 664                                 const int theNode2Index) const
 665 {
 666   if (myVolume->IsPoly()) {
 667     return IsLinked(myVolumeNodes[theNode1Index], myVolumeNodes[theNode2Index]);
 668   }
 669
 670   int minInd = theNode1Index < theNode2Index ? theNode1Index : theNode2Index;
 671   int maxInd = theNode1Index < theNode2Index ? theNode2Index : theNode1Index;
 672
 673   if ( minInd < 0 || maxInd > myVolumeNbNodes - 1 || maxInd == minInd )
 674     return false;
 675
 676   switch ( myVolumeNbNodes ) {
 677   case 4:
 678     return true;
 679   case 5:
 680     if ( maxInd == 4 )
 681       return true;
 682     switch ( maxInd - minInd ) {
 683     case 1:
 684     case 3: return true;
 685     default:;
 686     }
 687     break;
 688   case 6:
 689     switch ( maxInd - minInd ) {
 690     case 1: return minInd != 2;
 691     case 2: return minInd == 0 || minInd == 3;
 692     case 3: return true;
 693     default:;
 694     }
 695     break;
 696   case 8:
 697     switch ( maxInd - minInd ) {
 698     case 1: return minInd != 3;
 699     case 3: return minInd == 0 || minInd == 4;
 700     case 4: return true;
 701     default:;
 702     }
 703     break;
 704   default:;
 705   }
 706   return false;
 707 }
 708
 709 //=======================================================================
 710 //function : GetNodeIndex
 711 //purpose  : Return an index of theNode
 712 //=======================================================================
 713
 714 int SMDS_VolumeTool::GetNodeIndex(const SMDS_MeshNode* theNode) const
 715 {
 716   if ( myVolume ) {
 717     for ( int i = 0; i < myVolumeNbNodes; i++ ) {
 718       if ( myVolumeNodes[ i ] == theNode )
 719         return i;
 720     }
 721   }
 722   return -1;
 723 }
 724
 725 //=======================================================================
 726 //function : IsFreeFace
 727 //purpose  : check that only one volume is build on the face nodes
 728 //=======================================================================
 729
 730 bool SMDS_VolumeTool::IsFreeFace( int faceIndex )
 731 {
 732   const int free = true;
 733
 734   if (!setFace( faceIndex ))
 735     return !free;
 736
 737   const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
 738   int nbFaceNodes = myFaceNbNodes;
 739
 740   // evaluate nb of face nodes shared by other volume
 741   int maxNbShared = -1;
 742   typedef map< const SMDS_MeshElement*, int > TElemIntMap;
 743   TElemIntMap volNbShared;
 744   TElemIntMap::iterator vNbIt;
 745   for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
 746   {
 747     const SMDS_MeshNode* n = nodes[ iNode ];
 748     SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator();
 749     while ( eIt->more() ) {
 750       const SMDS_MeshElement* elem = eIt->next();
 751       if ( elem != myVolume && elem->GetType() == SMDSAbs_Volume ) {
 752         int nbShared = 1;
 753         vNbIt = volNbShared.find( elem );
 754         if ( vNbIt == volNbShared.end() )
 755           volNbShared.insert ( TElemIntMap::value_type( elem, nbShared ));
 756         else
 757           nbShared = ++(*vNbIt).second;
 758         if ( nbShared > maxNbShared )
 759           maxNbShared = nbShared;
 760       }
 761     }
 762   }
 763   if ( maxNbShared < 3 )
 764     return free; // is free
 765
 766   // find volumes laying on the opposite side of the face
 767   // and sharing all nodes
 768   XYZ intNormal; // internal normal
 769   GetFaceNormal( faceIndex, intNormal.x, intNormal.y, intNormal.z );
 770   if ( IsFaceExternal( faceIndex ))
 771     intNormal = XYZ( -intNormal.x, -intNormal.y, -intNormal.z );
 772   XYZ p0 ( nodes[0] ), baryCenter;
 773   for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); vNbIt++ )
 774   {
 775     int nbShared = (*vNbIt).second;
 776     if ( nbShared >= 3 ) {
 777       SMDS_VolumeTool volume( (*vNbIt).first );
 778       volume.GetBaryCenter( baryCenter.x, baryCenter.y, baryCenter.z );
 779       XYZ intNormal2( baryCenter - p0 );
 780       if ( intNormal.Dot( intNormal2 ) < 0 )
 781         continue; // opposite side
 782     }
 783     // remove a volume from volNbShared map
 784     volNbShared.erase( vNbIt );
 785   }
 786   // here volNbShared contains only volumes laying on the
 787   // opposite side of the face
 788   if ( volNbShared.empty() )
 789     return free; // is free
 790
 791   // check if the whole area of a face is shared
 792   bool isShared[] = { false, false, false, false }; // 4 triangle parts of a quadrangle
 793   for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); vNbIt++ )
 794   {
 795     SMDS_VolumeTool volume( (*vNbIt).first );
 796     bool prevLinkShared = false;
 797     int nbSharedLinks = 0;
 798     for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
 799     {
 800       bool linkShared = volume.IsLinked( nodes[ iNode ], nodes[ iNode + 1] );
 801       if ( linkShared )
 802         nbSharedLinks++;
 803       if ( linkShared && prevLinkShared &&
 804           volume.IsLinked( nodes[ iNode - 1 ], nodes[ iNode + 1] ))
 805         isShared[ iNode ] = true;
 806       prevLinkShared = linkShared;
 807     }
 808     if ( nbSharedLinks == nbFaceNodes )
 809       return !free; // is not free
 810     if ( nbFaceNodes == 4 ) {
 811       // check traingle parts 1 & 3
 812       if ( isShared[1] && isShared[3] )
 813         return !free; // is not free
 814       // check triangle parts 0 & 2;
 815       // 0 part could not be checked in the loop; check it here
 816       if ( isShared[2] && prevLinkShared &&
 817           volume.IsLinked( nodes[ 0 ], nodes[ 1 ] ) &&
 818           volume.IsLinked( nodes[ 1 ], nodes[ 3 ] ) )
 819         return !free; // is not free
 820     }
 821   }
 822   return free;
 823 }
 824
 825 //=======================================================================
 826 //function : GetFaceIndex
 827 //purpose  : Return index of a face formed by theFaceNodes
 828 //=======================================================================
 829
 830 int SMDS_VolumeTool::GetFaceIndex( const set<const SMDS_MeshNode*>& theFaceNodes )
 831 {
 832   for ( int iFace = 0; iFace < myNbFaces; iFace++ ) {
 833     const SMDS_MeshNode** nodes = GetFaceNodes( iFace );
 834     int nbFaceNodes = NbFaceNodes( iFace );
 835     set<const SMDS_MeshNode*> nodeSet;
 836     for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
 837       nodeSet.insert( nodes[ iNode ] );
 838     if ( theFaceNodes == nodeSet )
 839       return iFace;
 840   }
 841   return -1;
 842 }
 843
 844 //=======================================================================
 845 //function : GetFaceIndex
 846 //purpose  : Return index of a face formed by theFaceNodes
 847 //=======================================================================
 848
 849 /*int SMDS_VolumeTool::GetFaceIndex( const set<int>& theFaceNodesIndices )
 850 {
 851   for ( int iFace = 0; iFace < myNbFaces; iFace++ ) {
 852     const int* nodes = GetFaceNodesIndices( iFace );
 853     int nbFaceNodes = NbFaceNodes( iFace );
 854     set<int> nodeSet;
 855     for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
 856       nodeSet.insert( nodes[ iNode ] );
 857     if ( theFaceNodesIndices == nodeSet )
 858       return iFace;
 859   }
 860   return -1;
 861 }*/
 862
 863 //=======================================================================
 864 //function : setFace
 865 //purpose  :
 866 //=======================================================================
 867
 868 bool SMDS_VolumeTool::setFace( int faceIndex )
 869 {
 870   if ( !myVolume )
 871     return false;
 872
 873   if ( myCurFace == faceIndex )
 874     return true;
 875
 876   myCurFace = -1;
 877
 878   if ( faceIndex < 0 || faceIndex >= NbFaces() )
 879     return false;
 880
 881   if (myFaceNodes != NULL) {
 882     delete [] myFaceNodes;
 883     myFaceNodes = NULL;
 884   }
 885
 886   if (myVolume->IsPoly()) {
 887     if (!myPolyedre) {
 888       MESSAGE("Warning: bad volumic element");
 889       return false;
 890     }
 891
 892     // check orientation
 893     bool isGoodOri = true;
 894     if (myExternalFaces) {
 895       // get natural orientation
 896       XYZ aNormal, baryCenter, p0 (myPolyedre->GetFaceNode(faceIndex + 1, 1));
 897       SMDS_VolumeTool vTool (myPolyedre);
 898       vTool.GetFaceNormal(faceIndex, aNormal.x, aNormal.y, aNormal.z);
 899       vTool.GetBaryCenter(baryCenter.x, baryCenter.y, baryCenter.z);
 900       XYZ insideVec (baryCenter - p0);
 901       if (insideVec.Dot(aNormal) > 0)
 902         isGoodOri = false;
 903     }
 904
 905     // set face nodes
 906     int iNode;
 907     myFaceNbNodes = myPolyedre->NbFaceNodes(faceIndex + 1);
 908     myFaceNodes = new const SMDS_MeshNode* [myFaceNbNodes + 1];
 909     if (isGoodOri) {
 910       for ( iNode = 0; iNode < myFaceNbNodes; iNode++ )
 911         myFaceNodes[ iNode ] = myPolyedre->GetFaceNode(faceIndex + 1, iNode + 1);
 912     } else {
 913       for ( iNode = 0; iNode < myFaceNbNodes; iNode++ )
 914         myFaceNodes[ iNode ] = myPolyedre->GetFaceNode(faceIndex + 1, myFaceNbNodes - iNode);
 915     }
 916     myFaceNodes[ myFaceNbNodes ] = myFaceNodes[ 0 ]; // last = first
 917
 918   } else {
 919     // choose face node indices
 920     switch ( myVolumeNbNodes ) {
 921     case 4:
 922       myFaceNbNodes = Tetra_nbN[ faceIndex ];
 923       if ( myExternalFaces )
 924         myFaceNodeIndices = myVolForward ? Tetra_F[ faceIndex ] : Tetra_RE[ faceIndex ];
 925       else
 926         myFaceNodeIndices = myVolForward ? Tetra_F[ faceIndex ] : Tetra_R[ faceIndex ];
 927       break;
 928     case 5:
 929       myFaceNbNodes = Pyramid_nbN[ faceIndex ];
 930       if ( myExternalFaces )
 931         myFaceNodeIndices = myVolForward ? Pyramid_F[ faceIndex ] : Pyramid_RE[ faceIndex ];
 932       else
 933         myFaceNodeIndices = myVolForward ? Pyramid_F[ faceIndex ] : Pyramid_R[ faceIndex ];
 934       break;
 935     case 6:
 936       myFaceNbNodes = Penta_nbN[ faceIndex ];
 937       if ( myExternalFaces )
 938         myFaceNodeIndices = myVolForward ? Penta_FE[ faceIndex ] : Penta_RE[ faceIndex ];
 939       else
 940         myFaceNodeIndices = myVolForward ? Penta_F[ faceIndex ] : Penta_R[ faceIndex ];
 941       break;
 942     case 8:
 943       myFaceNbNodes = Hexa_nbN[ faceIndex ];
 944       if ( myExternalFaces )
 945         myFaceNodeIndices = myVolForward ? Hexa_FE[ faceIndex ] : Hexa_RE[ faceIndex ];
 946       else
 947         myFaceNodeIndices = Hexa_F[ faceIndex ];
 948       break;
 949     default:
 950       return false;
 951     }
 952
 953     // set face nodes
 954     myFaceNodes = new const SMDS_MeshNode* [myFaceNbNodes + 1];
 955     for ( int iNode = 0; iNode <= myFaceNbNodes; iNode++ )
 956       myFaceNodes[ iNode ] = myVolumeNodes[ myFaceNodeIndices[ iNode ]];
 957   }
 958
 959   myCurFace = faceIndex;
 960
 961   return true;
 962 }
 963
 964 //=======================================================================
 965 //function : GetType
 966 //purpose  : return VolumeType by nb of nodes in a volume
 967 //=======================================================================
 968
 969 SMDS_VolumeTool::VolumeType SMDS_VolumeTool::GetType(int nbNodes)
 970 {
 971   switch ( nbNodes ) {
 972   case 4: return TETRA;
 973   case 5: return PYRAM;
 974   case 6: return PENTA;
 975   case 8: return HEXA;
 976   default:return UNKNOWN;
 977   }
 978 }
 979
 980 //=======================================================================
 981 //function : NbFaces
 982 //purpose  : return nb of faces by volume type
 983 //=======================================================================
 984
 985 int SMDS_VolumeTool::NbFaces( VolumeType type )
 986 {
 987   switch ( type ) {
 988   case TETRA: return 4;
 989   case PYRAM: return 5;
 990   case PENTA: return 5;
 991   case HEXA : return 6;
 992   default:    return 0;
 993   }
 994 }
 995
 996 //=======================================================================
 997 //function : GetFaceNodesIndices
 998 //purpose  : Return the array of face nodes indices
 999 //           To comfort link iteration, the array
1000 //           length == NbFaceNodes( faceIndex ) + 1 and
1001 //           the last node index == the first one.
1002 //=======================================================================
1003
1004 const int* SMDS_VolumeTool::GetFaceNodesIndices(VolumeType type,
1005                                                 int        faceIndex,
1006                                                 bool       external)
1007 {
1008   switch ( type ) {
1009   case TETRA: return Tetra_F[ faceIndex ];
1010   case PYRAM: return Pyramid_F[ faceIndex ];
1011   case PENTA: return external ? Penta_FE[ faceIndex ] : Penta_F[ faceIndex ];
1012   case HEXA:  return external ? Hexa_FE[ faceIndex ] : Hexa_F[ faceIndex ];
1013   default:;
1014   }
1015   return 0;
1016 }
1017
1018 //=======================================================================
1019 //function : NbFaceNodes
1020 //purpose  : Return number of nodes in the array of face nodes
1021 //=======================================================================
1022
1023 int SMDS_VolumeTool::NbFaceNodes(VolumeType type,
1024                                  int        faceIndex )
1025 {
1026   switch ( type ) {
1027   case TETRA: return Tetra_nbN[ faceIndex ];
1028   case PYRAM: return Pyramid_nbN[ faceIndex ];
1029   case PENTA: return Penta_nbN[ faceIndex ];
1030   case HEXA:  return Hexa_nbN[ faceIndex ];
1031   default:;
1032   }
1033   return 0;
1034 }
1035