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 : GetVolumeType
 377 //purpose  :
 378 //=======================================================================
 379
 380 SMDS_VolumeTool::VolumeType SMDS_VolumeTool::GetVolumeType() const
 381 {
 382   if ( myPolyedre )
 383     return POLYHEDA;
 384
 385   if ( myVolume ) {
 386     static const VolumeType types[] = {
 387       TETRA,     // myVolumeNbNodes = 4
 388       PYRAM,     // myVolumeNbNodes = 5
 389       PENTA,     // myVolumeNbNodes = 6
 390       UNKNOWN,   // myVolumeNbNodes = 7
 391       HEXA       // myVolumeNbNodes = 8
 392     };
 393     return types[ myVolumeNbNodes - 4 ];
 394   }
 395
 396   return UNKNOWN;
 397 }
 398
 399 //=======================================================================
 400 //function : getTetraVolume
 401 //purpose  :
 402 //=======================================================================
 403
 404 static double getTetraVolume(const SMDS_MeshNode* n1,
 405                              const SMDS_MeshNode* n2,
 406                              const SMDS_MeshNode* n3,
 407                              const SMDS_MeshNode* n4)
 408 {
 409   double X1 = n1->X();
 410   double Y1 = n1->Y();
 411   double Z1 = n1->Z();
 412
 413   double X2 = n2->X();
 414   double Y2 = n2->Y();
 415   double Z2 = n2->Z();
 416
 417   double X3 = n3->X();
 418   double Y3 = n3->Y();
 419   double Z3 = n3->Z();
 420
 421   double X4 = n4->X();
 422   double Y4 = n4->Y();
 423   double Z4 = n4->Z();
 424
 425   double Q1 = -(X1-X2)*(Y3*Z4-Y4*Z3);
 426   double Q2 =  (X1-X3)*(Y2*Z4-Y4*Z2);
 427   double R1 = -(X1-X4)*(Y2*Z3-Y3*Z2);
 428   double R2 = -(X2-X3)*(Y1*Z4-Y4*Z1);
 429   double S1 =  (X2-X4)*(Y1*Z3-Y3*Z1);
 430   double S2 = -(X3-X4)*(Y1*Z2-Y2*Z1);
 431
 432   return (Q1+Q2+R1+R2+S1+S2)/6.0;
 433 }
 434
 435 //=======================================================================
 436 //function : GetSize
 437 //purpose  : Return element volume
 438 //=======================================================================
 439
 440 double SMDS_VolumeTool::GetSize() const
 441 {
 442   double V = 0.;
 443   if ( !myVolume )
 444     return 0.;
 445
 446   if ( myVolume->IsPoly() )
 447   {
 448     if ( !myPolyedre )
 449       return 0.;
 450
 451     // split a polyhedron into tetrahedrons
 452
 453     SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* > ( this );
 454     XYZ baryCenter;
 455     me->GetBaryCenter(baryCenter.x, baryCenter.y, baryCenter.z);
 456     SMDS_MeshNode bcNode ( baryCenter.x, baryCenter.y, baryCenter.z );
 457
 458     for ( int f = 0; f < NbFaces(); ++f )
 459     {
 460       bool externalFace = me->IsFaceExternal( f ); // it calls setFace()
 461       for ( int n = 2; n < myFaceNbNodes; ++n )
 462       {
 463         double Vn = getTetraVolume( myFaceNodes[ 0 ],
 464                                     myFaceNodes[ n-1 ],
 465                                     myFaceNodes[ n ],
 466                                     & bcNode );
 467 ///         cout <<"++++   " << Vn << "   nodes " <<myFaceNodes[ 0 ]->GetID() << " " <<myFaceNodes[ n-1 ]->GetID() << " " <<myFaceNodes[ n ]->GetID() << "        < " << V << endl;
 468         V += externalFace ? -Vn : Vn;
 469       }
 470     }
 471   }
 472   else
 473   {
 474     const static int ind[] = {
 475       0, 1, 3, 6, 11 };
 476     const static int vtab[][4] = {
 477       // tetrahedron
 478       { 0, 1, 2, 3 },
 479       // pyramid
 480       { 0, 1, 3, 4 },
 481       { 1, 2, 3, 4 },
 482       // pentahedron
 483       { 0, 1, 2, 3 },
 484       { 1, 5, 3, 4 },
 485       { 1, 5, 2, 3 },
 486       // hexahedron
 487       { 1, 4, 3, 0 },
 488       { 4, 1, 6, 5 },
 489       { 1, 3, 6, 2 },
 490       { 4, 6, 3, 7 },
 491       { 1, 4, 6, 3 }
 492     };
 493
 494     int type = GetVolumeType();
 495     int n1 = ind[type];
 496     int n2 = ind[type+1];
 497
 498     for (int i = n1; i <  n2; i++)
 499     {
 500       V -= getTetraVolume( myVolumeNodes[ vtab[i][0] ],
 501                            myVolumeNodes[ vtab[i][1] ],
 502                            myVolumeNodes[ vtab[i][2] ],
 503                            myVolumeNodes[ vtab[i][3] ]);
 504     }
 505   }
 506   return V;
 507 }
 508
 509 //=======================================================================
 510 //function : GetBaryCenter
 511 //purpose  :
 512 //=======================================================================
 513
 514 bool SMDS_VolumeTool::GetBaryCenter(double & X, double & Y, double & Z) const
 515 {
 516   X = Y = Z = 0.;
 517   if ( !myVolume )
 518     return false;
 519
 520   for ( int i = 0; i < myVolumeNbNodes; i++ ) {
 521     X += myVolumeNodes[ i ]->X();
 522     Y += myVolumeNodes[ i ]->Y();
 523     Z += myVolumeNodes[ i ]->Z();
 524   }
 525   X /= myVolumeNbNodes;
 526   Y /= myVolumeNbNodes;
 527   Z /= myVolumeNbNodes;
 528
 529   return true;
 530 }
 531
 532 //=======================================================================
 533 //function : SetExternalNormal
 534 //purpose  : Node order will be so that faces normals are external
 535 //=======================================================================
 536
 537 void SMDS_VolumeTool::SetExternalNormal ()
 538 {
 539   myExternalFaces = true;
 540   myCurFace = -1;
 541 }
 542
 543 //=======================================================================
 544 //function : NbFaceNodes
 545 //purpose  : Return number of nodes in the array of face nodes
 546 //=======================================================================
 547
 548 int SMDS_VolumeTool::NbFaceNodes( int faceIndex )
 549 {
 550     if ( !setFace( faceIndex ))
 551       return 0;
 552     return myFaceNbNodes;
 553 }
 554
 555 //=======================================================================
 556 //function : GetFaceNodes
 557 //purpose  : Return pointer to the array of face nodes.
 558 //           To comfort link iteration, the array
 559 //           length == NbFaceNodes( faceIndex ) + 1 and
 560 //           the last node == the first one.
 561 //=======================================================================
 562
 563 const SMDS_MeshNode** SMDS_VolumeTool::GetFaceNodes( int faceIndex )
 564 {
 565   if ( !setFace( faceIndex ))
 566     return 0;
 567   return myFaceNodes;
 568 }
 569
 570 //=======================================================================
 571 //function : GetFaceNodesIndices
 572 //purpose  : Return pointer to the array of face nodes indices
 573 //           To comfort link iteration, the array
 574 //           length == NbFaceNodes( faceIndex ) + 1 and
 575 //           the last node index == the first one.
 576 //=======================================================================
 577
 578 const int* SMDS_VolumeTool::GetFaceNodesIndices( int faceIndex )
 579 {
 580   if (myVolume->IsPoly()) {
 581     MESSAGE("Warning: attempt to obtain FaceNodesIndices of polyhedral volume");
 582     return NULL;
 583   }
 584   if ( !setFace( faceIndex ))
 585     return 0;
 586   return myFaceNodeIndices;
 587 }
 588
 589 //=======================================================================
 590 //function : GetFaceNodes
 591 //purpose  : Return a set of face nodes.
 592 //=======================================================================
 593
 594 bool SMDS_VolumeTool::GetFaceNodes (int                        faceIndex,
 595                                     set<const SMDS_MeshNode*>& theFaceNodes )
 596 {
 597   if ( !setFace( faceIndex ))
 598     return false;
 599
 600   theFaceNodes.clear();
 601   int iNode, nbNode = myFaceNbNodes;
 602   for ( iNode = 0; iNode < nbNode; iNode++ )
 603     theFaceNodes.insert( myFaceNodes[ iNode ]);
 604
 605   return true;
 606 }
 607
 608 //=======================================================================
 609 //function : IsFaceExternal
 610 //purpose  : Check normal orientation of a returned face
 611 //=======================================================================
 612
 613 bool SMDS_VolumeTool::IsFaceExternal( int faceIndex )
 614 {
 615   if ( myExternalFaces || !myVolume )
 616     return true;
 617
 618   if (myVolume->IsPoly()) {
 619     XYZ aNormal, baryCenter, p0 (myPolyedre->GetFaceNode(faceIndex + 1, 1));
 620     GetFaceNormal(faceIndex, aNormal.x, aNormal.y, aNormal.z);
 621     GetBaryCenter(baryCenter.x, baryCenter.y, baryCenter.z);
 622     XYZ insideVec (baryCenter - p0);
 623     if (insideVec.Dot(aNormal) > 0)
 624       return false;
 625     return true;
 626   }
 627
 628   switch ( myVolumeNbNodes ) {
 629   case 4:
 630   case 5:
 631     // only the bottom of a reversed tetrahedron can be internal
 632     return ( myVolForward || faceIndex != 0 );
 633   case 6:
 634     // in a forward pentahedron, the top is internal, in a reversed one - bottom
 635     return ( myVolForward ? faceIndex != 1 : faceIndex != 0 );
 636   case 8: {
 637     // in a forward hexahedron, even face normal is external, odd - internal
 638     bool odd = faceIndex % 2;
 639     return ( myVolForward ? !odd : odd );
 640   }
 641   default:;
 642   }
 643   return false;
 644 }
 645
 646 //=======================================================================
 647 //function : GetFaceNormal
 648 //purpose  : Return a normal to a face
 649 //=======================================================================
 650
 651 bool SMDS_VolumeTool::GetFaceNormal (int faceIndex, double & X, double & Y, double & Z)
 652 {
 653   if ( !setFace( faceIndex ))
 654     return false;
 655
 656   XYZ p1 ( myFaceNodes[0] );
 657   XYZ p2 ( myFaceNodes[1] );
 658   XYZ p3 ( myFaceNodes[2] );
 659   XYZ aVec12( p2 - p1 );
 660   XYZ aVec13( p3 - p1 );
 661   XYZ cross = aVec12.Crossed( aVec13 );
 662
 663   if ( myFaceNbNodes == 4 ) {
 664     XYZ p4 ( myFaceNodes[3] );
 665     XYZ aVec14( p4 - p1 );
 666     XYZ cross2 = aVec13.Crossed( aVec14 );
 667     cross.x += cross2.x;
 668     cross.y += cross2.y;
 669     cross.z += cross2.z;
 670   }
 671
 672   double size = cross.Magnitude();
 673   if ( size <= DBL_MIN )
 674     return false;
 675
 676   X = cross.x / size;
 677   Y = cross.y / size;
 678   Z = cross.z / size;
 679
 680   return true;
 681 }
 682
 683 //=======================================================================
 684 //function : GetFaceArea
 685 //purpose  : Return face area
 686 //=======================================================================
 687
 688 double SMDS_VolumeTool::GetFaceArea( int faceIndex )
 689 {
 690   if (myVolume->IsPoly()) {
 691     MESSAGE("Warning: attempt to obtain area of a face of polyhedral volume");
 692     return 0;
 693   }
 694
 695   if ( !setFace( faceIndex ))
 696     return 0;
 697
 698   XYZ p1 ( myFaceNodes[0] );
 699   XYZ p2 ( myFaceNodes[1] );
 700   XYZ p3 ( myFaceNodes[2] );
 701   XYZ aVec12( p2 - p1 );
 702   XYZ aVec13( p3 - p1 );
 703   double area = aVec12.Crossed( aVec13 ).Magnitude() * 0.5;
 704
 705   if ( myFaceNbNodes == 4 ) {
 706     XYZ p4 ( myFaceNodes[3] );
 707     XYZ aVec14( p4 - p1 );
 708     area += aVec14.Crossed( aVec13 ).Magnitude() * 0.5;
 709   }
 710   return area;
 711 }
 712
 713 //=======================================================================
 714 //function : GetOppFaceIndex
 715 //purpose  : Return index of the opposite face if it exists, else -1.
 716 //=======================================================================
 717
 718 int SMDS_VolumeTool::GetOppFaceIndex( int faceIndex ) const
 719 {
 720   int ind = -1;
 721   if (myVolume->IsPoly()) {
 722     MESSAGE("Warning: attempt to obtain opposite face on polyhedral volume");
 723     return ind;
 724   }
 725
 726   if ( faceIndex >= 0 && faceIndex < NbFaces() ) {
 727     switch ( myVolumeNbNodes ) {
 728     case 6:
 729       if ( faceIndex == 0 || faceIndex == 1 )
 730         ind = 1 - faceIndex;
 731         break;
 732     case 8:
 733       ind = faceIndex + ( faceIndex % 2 ? -1 : 1 );
 734       break;
 735     default:;
 736     }
 737   }
 738   return ind;
 739 }
 740
 741 //=======================================================================
 742 //function : IsLinked
 743 //purpose  : return true if theNode1 is linked with theNode2
 744 //=======================================================================
 745
 746 bool SMDS_VolumeTool::IsLinked (const SMDS_MeshNode* theNode1,
 747                                 const SMDS_MeshNode* theNode2) const
 748 {
 749   if ( !myVolume )
 750     return false;
 751
 752   if (myVolume->IsPoly()) {
 753     if (!myPolyedre) {
 754       MESSAGE("Warning: bad volumic element");
 755       return false;
 756     }
 757     bool isLinked = false;
 758     int iface;
 759     for (iface = 1; iface <= myNbFaces && !isLinked; iface++) {
 760       int inode, nbFaceNodes = myPolyedre->NbFaceNodes(iface);
 761
 762       for (inode = 1; inode <= nbFaceNodes && !isLinked; inode++) {
 763         const SMDS_MeshNode* curNode = myPolyedre->GetFaceNode(iface, inode);
 764
 765         if (curNode == theNode1 || curNode == theNode2) {
 766           int inextnode = (inode == nbFaceNodes) ? 1 : inode + 1;
 767           const SMDS_MeshNode* nextNode = myPolyedre->GetFaceNode(iface, inextnode);
 768
 769           if ((curNode == theNode1 && nextNode == theNode2) ||
 770               (curNode == theNode2 && nextNode == theNode1)) {
 771             isLinked = true;
 772           }
 773         }
 774       }
 775     }
 776     return isLinked;
 777   }
 778
 779   // find nodes indices
 780   int i1 = -1, i2 = -1;
 781   for ( int i = 0; i < myVolumeNbNodes; i++ ) {
 782     if ( myVolumeNodes[ i ] == theNode1 )
 783       i1 = i;
 784     else if ( myVolumeNodes[ i ] == theNode2 )
 785       i2 = i;
 786   }
 787   return IsLinked( i1, i2 );
 788 }
 789
 790 //=======================================================================
 791 //function : IsLinked
 792 //purpose  : return true if the node with theNode1Index is linked
 793 //           with the node with theNode2Index
 794 //=======================================================================
 795
 796 bool SMDS_VolumeTool::IsLinked (const int theNode1Index,
 797                                 const int theNode2Index) const
 798 {
 799   if (myVolume->IsPoly()) {
 800     return IsLinked(myVolumeNodes[theNode1Index], myVolumeNodes[theNode2Index]);
 801   }
 802
 803   int minInd = theNode1Index < theNode2Index ? theNode1Index : theNode2Index;
 804   int maxInd = theNode1Index < theNode2Index ? theNode2Index : theNode1Index;
 805
 806   if ( minInd < 0 || maxInd > myVolumeNbNodes - 1 || maxInd == minInd )
 807     return false;
 808
 809   switch ( myVolumeNbNodes ) {
 810   case 4:
 811     return true;
 812   case 5:
 813     if ( maxInd == 4 )
 814       return true;
 815     switch ( maxInd - minInd ) {
 816     case 1:
 817     case 3: return true;
 818     default:;
 819     }
 820     break;
 821   case 6:
 822     switch ( maxInd - minInd ) {
 823     case 1: return minInd != 2;
 824     case 2: return minInd == 0 || minInd == 3;
 825     case 3: return true;
 826     default:;
 827     }
 828     break;
 829   case 8:
 830     switch ( maxInd - minInd ) {
 831     case 1: return minInd != 3;
 832     case 3: return minInd == 0 || minInd == 4;
 833     case 4: return true;
 834     default:;
 835     }
 836     break;
 837   default:;
 838   }
 839   return false;
 840 }
 841
 842 //=======================================================================
 843 //function : GetNodeIndex
 844 //purpose  : Return an index of theNode
 845 //=======================================================================
 846
 847 int SMDS_VolumeTool::GetNodeIndex(const SMDS_MeshNode* theNode) const
 848 {
 849   if ( myVolume ) {
 850     for ( int i = 0; i < myVolumeNbNodes; i++ ) {
 851       if ( myVolumeNodes[ i ] == theNode )
 852         return i;
 853     }
 854   }
 855   return -1;
 856 }
 857
 858 //=======================================================================
 859 //function : IsFreeFace
 860 //purpose  : check that only one volume is build on the face nodes
 861 //=======================================================================
 862
 863 bool SMDS_VolumeTool::IsFreeFace( int faceIndex )
 864 {
 865   const int free = true;
 866
 867   if (!setFace( faceIndex ))
 868     return !free;
 869
 870   const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
 871   int nbFaceNodes = myFaceNbNodes;
 872
 873   // evaluate nb of face nodes shared by other volume
 874   int maxNbShared = -1;
 875   typedef map< const SMDS_MeshElement*, int > TElemIntMap;
 876   TElemIntMap volNbShared;
 877   TElemIntMap::iterator vNbIt;
 878   for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
 879   {
 880     const SMDS_MeshNode* n = nodes[ iNode ];
 881     SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator();
 882     while ( eIt->more() ) {
 883       const SMDS_MeshElement* elem = eIt->next();
 884       if ( elem != myVolume && elem->GetType() == SMDSAbs_Volume ) {
 885         int nbShared = 1;
 886         vNbIt = volNbShared.find( elem );
 887         if ( vNbIt == volNbShared.end() )
 888           volNbShared.insert ( TElemIntMap::value_type( elem, nbShared ));
 889         else
 890           nbShared = ++(*vNbIt).second;
 891         if ( nbShared > maxNbShared )
 892           maxNbShared = nbShared;
 893       }
 894     }
 895   }
 896   if ( maxNbShared < 3 )
 897     return free; // is free
 898
 899   // find volumes laying on the opposite side of the face
 900   // and sharing all nodes
 901   XYZ intNormal; // internal normal
 902   GetFaceNormal( faceIndex, intNormal.x, intNormal.y, intNormal.z );
 903   if ( IsFaceExternal( faceIndex ))
 904     intNormal = XYZ( -intNormal.x, -intNormal.y, -intNormal.z );
 905   XYZ p0 ( nodes[0] ), baryCenter;
 906   for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); vNbIt++ )
 907   {
 908     int nbShared = (*vNbIt).second;
 909     if ( nbShared >= 3 ) {
 910       SMDS_VolumeTool volume( (*vNbIt).first );
 911       volume.GetBaryCenter( baryCenter.x, baryCenter.y, baryCenter.z );
 912       XYZ intNormal2( baryCenter - p0 );
 913       if ( intNormal.Dot( intNormal2 ) < 0 )
 914         continue; // opposite side
 915     }
 916     // remove a volume from volNbShared map
 917     volNbShared.erase( vNbIt );
 918   }
 919   // here volNbShared contains only volumes laying on the
 920   // opposite side of the face
 921   if ( volNbShared.empty() )
 922     return free; // is free
 923
 924   // check if the whole area of a face is shared
 925   bool isShared[] = { false, false, false, false }; // 4 triangle parts of a quadrangle
 926   for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); vNbIt++ )
 927   {
 928     SMDS_VolumeTool volume( (*vNbIt).first );
 929     bool prevLinkShared = false;
 930     int nbSharedLinks = 0;
 931     for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
 932     {
 933       bool linkShared = volume.IsLinked( nodes[ iNode ], nodes[ iNode + 1] );
 934       if ( linkShared )
 935         nbSharedLinks++;
 936       if ( linkShared && prevLinkShared &&
 937           volume.IsLinked( nodes[ iNode - 1 ], nodes[ iNode + 1] ))
 938         isShared[ iNode ] = true;
 939       prevLinkShared = linkShared;
 940     }
 941     if ( nbSharedLinks == nbFaceNodes )
 942       return !free; // is not free
 943     if ( nbFaceNodes == 4 ) {
 944       // check traingle parts 1 & 3
 945       if ( isShared[1] && isShared[3] )
 946         return !free; // is not free
 947       // check triangle parts 0 & 2;
 948       // 0 part could not be checked in the loop; check it here
 949       if ( isShared[2] && prevLinkShared &&
 950           volume.IsLinked( nodes[ 0 ], nodes[ 1 ] ) &&
 951           volume.IsLinked( nodes[ 1 ], nodes[ 3 ] ) )
 952         return !free; // is not free
 953     }
 954   }
 955   return free;
 956 }
 957
 958 //=======================================================================
 959 //function : GetFaceIndex
 960 //purpose  : Return index of a face formed by theFaceNodes
 961 //=======================================================================
 962
 963 int SMDS_VolumeTool::GetFaceIndex( const set<const SMDS_MeshNode*>& theFaceNodes )
 964 {
 965   for ( int iFace = 0; iFace < myNbFaces; iFace++ ) {
 966     const SMDS_MeshNode** nodes = GetFaceNodes( iFace );
 967     int nbFaceNodes = NbFaceNodes( iFace );
 968     set<const SMDS_MeshNode*> nodeSet;
 969     for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
 970       nodeSet.insert( nodes[ iNode ] );
 971     if ( theFaceNodes == nodeSet )
 972       return iFace;
 973   }
 974   return -1;
 975 }
 976
 977 //=======================================================================
 978 //function : GetFaceIndex
 979 //purpose  : Return index of a face formed by theFaceNodes
 980 //=======================================================================
 981
 982 /*int SMDS_VolumeTool::GetFaceIndex( const set<int>& theFaceNodesIndices )
 983 {
 984   for ( int iFace = 0; iFace < myNbFaces; iFace++ ) {
 985     const int* nodes = GetFaceNodesIndices( iFace );
 986     int nbFaceNodes = NbFaceNodes( iFace );
 987     set<int> nodeSet;
 988     for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
 989       nodeSet.insert( nodes[ iNode ] );
 990     if ( theFaceNodesIndices == nodeSet )
 991       return iFace;
 992   }
 993   return -1;
 994 }*/
 995
 996 //=======================================================================
 997 //function : setFace
 998 //purpose  :
 999 //=======================================================================
1000
1001 bool SMDS_VolumeTool::setFace( int faceIndex )
1002 {
1003   if ( !myVolume )
1004     return false;
1005
1006   if ( myCurFace == faceIndex )
1007     return true;
1008
1009   myCurFace = -1;
1010
1011   if ( faceIndex < 0 || faceIndex >= NbFaces() )
1012     return false;
1013
1014   if (myFaceNodes != NULL) {
1015     delete [] myFaceNodes;
1016     myFaceNodes = NULL;
1017   }
1018
1019   if (myVolume->IsPoly()) {
1020     if (!myPolyedre) {
1021       MESSAGE("Warning: bad volumic element");
1022       return false;
1023     }
1024
1025     // check orientation
1026     bool isGoodOri = true;
1027     if (myExternalFaces)
1028       isGoodOri = IsFaceExternal( faceIndex );
1029
1030     // set face nodes
1031     int iNode;
1032     myFaceNbNodes = myPolyedre->NbFaceNodes(faceIndex + 1);
1033     myFaceNodes = new const SMDS_MeshNode* [myFaceNbNodes + 1];
1034     if (isGoodOri) {
1035       for ( iNode = 0; iNode < myFaceNbNodes; iNode++ )
1036         myFaceNodes[ iNode ] = myPolyedre->GetFaceNode(faceIndex + 1, iNode + 1);
1037     } else {
1038       for ( iNode = 0; iNode < myFaceNbNodes; iNode++ )
1039         myFaceNodes[ iNode ] = myPolyedre->GetFaceNode(faceIndex + 1, myFaceNbNodes - iNode);
1040     }
1041     myFaceNodes[ myFaceNbNodes ] = myFaceNodes[ 0 ]; // last = first
1042
1043   } else {
1044     // choose face node indices
1045     switch ( myVolumeNbNodes ) {
1046     case 4:
1047       myFaceNbNodes = Tetra_nbN[ faceIndex ];
1048       if ( myExternalFaces )
1049         myFaceNodeIndices = myVolForward ? Tetra_F[ faceIndex ] : Tetra_RE[ faceIndex ];
1050       else
1051         myFaceNodeIndices = myVolForward ? Tetra_F[ faceIndex ] : Tetra_R[ faceIndex ];
1052       break;
1053     case 5:
1054       myFaceNbNodes = Pyramid_nbN[ faceIndex ];
1055       if ( myExternalFaces )
1056         myFaceNodeIndices = myVolForward ? Pyramid_F[ faceIndex ] : Pyramid_RE[ faceIndex ];
1057       else
1058         myFaceNodeIndices = myVolForward ? Pyramid_F[ faceIndex ] : Pyramid_R[ faceIndex ];
1059       break;
1060     case 6:
1061       myFaceNbNodes = Penta_nbN[ faceIndex ];
1062       if ( myExternalFaces )
1063         myFaceNodeIndices = myVolForward ? Penta_FE[ faceIndex ] : Penta_RE[ faceIndex ];
1064       else
1065         myFaceNodeIndices = myVolForward ? Penta_F[ faceIndex ] : Penta_R[ faceIndex ];
1066       break;
1067     case 8:
1068       myFaceNbNodes = Hexa_nbN[ faceIndex ];
1069       if ( myExternalFaces )
1070         myFaceNodeIndices = myVolForward ? Hexa_FE[ faceIndex ] : Hexa_RE[ faceIndex ];
1071       else
1072         myFaceNodeIndices = Hexa_F[ faceIndex ];
1073       break;
1074     default:
1075       return false;
1076     }
1077
1078     // set face nodes
1079     myFaceNodes = new const SMDS_MeshNode* [myFaceNbNodes + 1];
1080     for ( int iNode = 0; iNode < myFaceNbNodes; iNode++ )
1081       myFaceNodes[ iNode ] = myVolumeNodes[ myFaceNodeIndices[ iNode ]];
1082     myFaceNodes[ myFaceNbNodes ] = myFaceNodes[ 0 ];
1083   }
1084
1085   myCurFace = faceIndex;
1086
1087   return true;
1088 }
1089
1090 //=======================================================================
1091 //function : GetType
1092 //purpose  : return VolumeType by nb of nodes in a volume
1093 //=======================================================================
1094
1095 SMDS_VolumeTool::VolumeType SMDS_VolumeTool::GetType(int nbNodes)
1096 {
1097   switch ( nbNodes ) {
1098   case 4: return TETRA;
1099   case 5: return PYRAM;
1100   case 6: return PENTA;
1101   case 8: return HEXA;
1102   default:return UNKNOWN;
1103   }
1104 }
1105
1106 //=======================================================================
1107 //function : NbFaces
1108 //purpose  : return nb of faces by volume type
1109 //=======================================================================
1110
1111 int SMDS_VolumeTool::NbFaces( VolumeType type )
1112 {
1113   switch ( type ) {
1114   case TETRA: return 4;
1115   case PYRAM: return 5;
1116   case PENTA: return 5;
1117   case HEXA : return 6;
1118   default:    return 0;
1119   }
1120 }
1121
1122 //=======================================================================
1123 //function : GetFaceNodesIndices
1124 //purpose  : Return the array of face nodes indices
1125 //           To comfort link iteration, the array
1126 //           length == NbFaceNodes( faceIndex ) + 1 and
1127 //           the last node index == the first one.
1128 //=======================================================================
1129
1130 const int* SMDS_VolumeTool::GetFaceNodesIndices(VolumeType type,
1131                                                 int        faceIndex,
1132                                                 bool       external)
1133 {
1134   switch ( type ) {
1135   case TETRA: return Tetra_F[ faceIndex ];
1136   case PYRAM: return Pyramid_F[ faceIndex ];
1137   case PENTA: return external ? Penta_FE[ faceIndex ] : Penta_F[ faceIndex ];
1138   case HEXA:  return external ? Hexa_FE[ faceIndex ] : Hexa_F[ faceIndex ];
1139   default:;
1140   }
1141   return 0;
1142 }
1143
1144 //=======================================================================
1145 //function : NbFaceNodes
1146 //purpose  : Return number of nodes in the array of face nodes
1147 //=======================================================================
1148
1149 int SMDS_VolumeTool::NbFaceNodes(VolumeType type,
1150                                  int        faceIndex )
1151 {
1152   switch ( type ) {
1153   case TETRA: return Tetra_nbN[ faceIndex ];
1154   case PYRAM: return Pyramid_nbN[ faceIndex ];
1155   case PENTA: return Penta_nbN[ faceIndex ];
1156   case HEXA:  return Hexa_nbN[ faceIndex ];
1157   default:;
1158   }
1159   return 0;
1160 }
1161