src/SMDS/SMDS_VolumeTool.cxx

   1 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 // File      : SMDS_VolumeTool.cxx
  24 // Created   : Tue Jul 13 12:22:13 2004
  25 // Author    : Edward AGAPOV (eap)
  26 //
  27 #ifdef _MSC_VER
  28 #pragma warning(disable:4786)
  29 #endif
  30
  31 #include "SMDS_VolumeTool.hxx"
  32
  33 #include "SMDS_MeshElement.hxx"
  34 #include "SMDS_MeshNode.hxx"
  35 #include "SMDS_VtkVolume.hxx"
  36 #include "SMDS_Mesh.hxx"
  37
  38 #include "utilities.h"
  39
  40 #include <map>
  41 #include <limits>
  42 #include <cmath>
  43
  44 using namespace std;
  45
  46 // ======================================================
  47 // Node indices in faces depending on volume orientation
  48 // making most faces normals external
  49 // ======================================================
  50 // For all elements, 0-th face is bottom based on the first nodes.
  51 // For prismatic elements (tetra,hexa,prisms), 1-th face is a top one.
  52 // For all elements, side faces follow order of bottom nodes
  53 // ======================================================
  54
  55 /*
  56 //           N3
  57 //           +
  58 //          /|\
  59 //         / | \
  60 //        /  |  \
  61 //    N0 +---|---+ N1                TETRAHEDRON
  62 //       \   |   /
  63 //        \  |  /
  64 //         \ | /
  65 //          \|/
  66 //           +
  67 //           N2
  68 */
  69 static int Tetra_F [4][4] = { // FORWARD == EXTERNAL
  70   { 0, 1, 2, 0 },              // All faces have external normals
  71   { 0, 3, 1, 0 },
  72   { 1, 3, 2, 1 },
  73   { 0, 2, 3, 0 }};
  74 static int Tetra_RE [4][4] = { // REVERSED -> FORWARD (EXTERNAL)
  75   { 0, 2, 1, 0 },              // All faces have external normals
  76   { 0, 1, 3, 0 },
  77   { 1, 2, 3, 1 },
  78   { 0, 3, 2, 0 }};
  79 static int Tetra_nbN [] = { 3, 3, 3, 3 };
  80
  81 //
  82 //     PYRAMID
  83 //
  84 static int Pyramid_F [5][5] = { // FORWARD == EXTERNAL
  85   { 0, 1, 2, 3, 0 },            // All faces have external normals
  86   { 0, 4, 1, 0, 4 },
  87   { 1, 4, 2, 1, 4 },
  88   { 2, 4, 3, 2, 4 },
  89   { 3, 4, 0, 3, 4 }
  90 };
  91 static int Pyramid_RE [5][5] = { // REVERSED -> FORWARD (EXTERNAL)
  92   { 0, 3, 2, 1, 0 },             // All faces but a bottom have external normals
  93   { 0, 1, 4, 0, 4 },
  94   { 1, 2, 4, 1, 4 },
  95   { 2, 3, 4, 2, 4 },
  96   { 3, 0, 4, 3, 4 }};
  97 static int Pyramid_nbN [] = { 4, 3, 3, 3, 3 };
  98
  99 /*
 100 //            + N4
 101 //           /|\
 102 //          / | \
 103 //         /  |  \
 104 //        /   |   \
 105 //    N3 +---------+ N5
 106 //       |    |    |
 107 //       |    + N1 |
 108 //       |   / \   |                PENTAHEDRON
 109 //       |  /   \  |
 110 //       | /     \ |
 111 //       |/       \|
 112 //    N0 +---------+ N2
 113 */
 114 static int Penta_F [5][5] = { // FORWARD
 115   { 0, 1, 2, 0, 0 },          // All faces have external normals
 116   { 3, 5, 4, 3, 3 },          // 0 is bottom, 1 is top face
 117   { 0, 3, 4, 1, 0 },
 118   { 1, 4, 5, 2, 1 },
 119   { 0, 2, 5, 3, 0 }};
 120 static int Penta_RE [5][5] = { // REVERSED -> EXTERNAL
 121   { 0, 2, 1, 0, 0 },
 122   { 3, 4, 5, 3, 3 },
 123   { 0, 1, 4, 3, 0 },
 124   { 1, 2, 5, 4, 1 },
 125   { 0, 3, 5, 2, 0 }};
 126 static int Penta_nbN [] = { 3, 3, 4, 4, 4 };
 127
 128 /*
 129 //         N5+----------+N6
 130 //          /|         /|
 131 //         / |        / |
 132 //        /  |       /  |
 133 //     N4+----------+N7 |
 134 //       |   |      |   |           HEXAHEDRON
 135 //       |   |      |   |
 136 //       |   |      |   |
 137 //       | N1+------|---+N2
 138 //       |  /       |  /
 139 //       | /        | /
 140 //       |/         |/
 141 //     N0+----------+N3
 142 */
 143 static int Hexa_F [6][5] = { // FORWARD
 144   { 0, 1, 2, 3, 0 },
 145   { 4, 7, 6, 5, 4 },          // all face normals are external
 146   { 0, 4, 5, 1, 0 },
 147   { 1, 5, 6, 2, 1 },
 148   { 3, 2, 6, 7, 3 },
 149   { 0, 3, 7, 4, 0 }};
 150 static int Hexa_RE [6][5] = { // REVERSED -> EXTERNAL
 151   { 0, 3, 2, 1, 0 },
 152   { 4, 5, 6, 7, 4 },          // all face normals are external
 153   { 0, 1, 5, 4, 0 },
 154   { 1, 2, 6, 5, 1 },
 155   { 3, 7, 6, 2, 3 },
 156   { 0, 4, 7, 3, 0 }};
 157 static int Hexa_nbN [] = { 4, 4, 4, 4, 4, 4 };
 158
 159 /*
 160 //      N8 +------+ N9
 161 //        /        \
 162 //       /          \
 163 //   N7 +            + N10
 164 //       \          /
 165 //        \        /
 166 //      N6 +------+ N11
 167 //                             HEXAGONAL PRISM
 168 //      N2 +------+ N3
 169 //        /        \
 170 //       /          \
 171 //   N1 +            + N4
 172 //       \          /
 173 //        \        /
 174 //      N0 +------+ N5
 175 */
 176 static int HexPrism_F [8][7] = { // FORWARD
 177   { 0, 1, 2, 3, 4, 5, 0 },
 178   { 6,11,10, 9, 8, 7, 6 },
 179   { 0, 6, 7, 1, 0, 0, 0 },
 180   { 1, 7, 8, 2, 1, 1, 1 },
 181   { 2, 8, 9, 3, 2, 2, 2 },
 182   { 3, 9,10, 4, 3, 3, 3 },
 183   { 4,10,11, 5, 4, 4, 4 },
 184   { 5,11, 6, 0, 5, 5, 5 }};
 185 static int HexPrism_RE [8][7] = { // REVERSED -> EXTERNAL
 186   { 0, 5, 4, 3, 2, 1, 0 },
 187   { 6,11,10, 9, 8, 7, 6 },
 188   { 0, 6, 7, 1, 0, 0, 0 },
 189   { 1, 7, 8, 2, 1, 1, 1 },
 190   { 2, 8, 9, 3, 2, 2, 2 },
 191   { 3, 9,10, 4, 3, 3, 3 },
 192   { 4,10,11, 5, 4, 4, 4 },
 193   { 5,11, 6, 0, 5, 5, 5 }};
 194 static int HexPrism_nbN [] = { 6, 6, 4, 4, 4, 4, 4, 4 };
 195
 196
 197 /*
 198 //           N3
 199 //           +
 200 //          /|\
 201 //        7/ | \8
 202 //        /  |4 \                    QUADRATIC
 203 //    N0 +---|---+ N1                TETRAHEDRON
 204 //       \   +9  /
 205 //        \  |  /
 206 //        6\ | /5
 207 //          \|/
 208 //           +
 209 //           N2
 210 */
 211 static int QuadTetra_F [4][7] = { // FORWARD
 212   { 0, 4, 1, 5, 2, 6, 0 },        // All faces have external normals
 213   { 0, 7, 3, 8, 1, 4, 0 },
 214   { 1, 8, 3, 9, 2, 5, 1 },
 215   { 0, 6, 2, 9, 3, 7, 0 }};
 216 static int QuadTetra_RE [4][7] = { // REVERSED -> FORWARD (EXTERNAL)
 217   { 0, 6, 2, 5, 1, 4, 0 },         // All faces have external normals
 218   { 0, 4, 1, 8, 3, 7, 0 },
 219   { 1, 5, 2, 9, 3, 8, 1 },
 220   { 0, 7, 3, 9, 2, 6, 0 }};
 221 static int QuadTetra_nbN [] = { 6, 6, 6, 6 };
 222
 223 //
 224 //     QUADRATIC
 225 //     PYRAMID
 226 //
 227 //            +4
 228 //
 229 //
 230 //       10+-----+11
 231 //         |     |        9 - middle point for (0,4) etc.
 232 //         |     |
 233 //        9+-----+12
 234 //
 235 //            6
 236 //      1+----+----+2
 237 //       |         |
 238 //       |         |
 239 //      5+         +7
 240 //       |         |
 241 //       |         |
 242 //      0+----+----+3
 243 //            8
 244 static int QuadPyram_F [5][9] = {  // FORWARD
 245   { 0, 5, 1, 6, 2, 7, 3, 8, 0 },   // All faces have external normals
 246   { 0, 9, 4, 10,1, 5, 0, 4, 4 },
 247   { 1, 10,4, 11,2, 6, 1, 4, 4 },
 248   { 2, 11,4, 12,3, 7, 2, 4, 4 },
 249   { 3, 12,4, 9, 0, 8, 3, 4, 4 }};
 250 static int QuadPyram_RE [5][9] = { // REVERSED -> FORWARD (EXTERNAL)
 251   { 0, 8, 3, 7, 2, 6, 1, 5, 0 },   // All faces but a bottom have external normals
 252   { 0, 5, 1, 10,4, 9, 0, 4, 4 },
 253   { 1, 6, 2, 11,4, 10,1, 4, 4 },
 254   { 2, 7, 3, 12,4, 11,2, 4, 4 },
 255   { 3, 8, 0, 9, 4, 12,3, 4, 4 }};
 256 static int QuadPyram_nbN [] = { 8, 6, 6, 6, 6 };
 257
 258 /*
 259 //            + N4
 260 //           /|\
 261 //         9/ | \10
 262 //         /  |  \
 263 //        /   |   \
 264 //    N3 +----+----+ N5
 265 //       |    |11  |
 266 //       |    |    |
 267 //       |    +13  |                QUADRATIC
 268 //       |    |    |                PENTAHEDRON
 269 //     12+    |    +14
 270 //       |    |    |
 271 //       |    |    |
 272 //       |    + N1 |
 273 //       |   / \   |
 274 //       | 6/   \7 |
 275 //       | /     \ |
 276 //       |/       \|
 277 //    N0 +---------+ N2
 278 //            8
 279 */
 280 static int QuadPenta_F [5][9] = {  // FORWARD
 281   { 0, 6, 1, 7, 2, 8, 0, 0, 0 },
 282   { 3,11, 5, 10,4, 9, 3, 3, 3 },
 283   { 0, 12,3, 9, 4, 13,1, 6, 0 },
 284   { 1, 13,4, 10,5, 14,2, 7, 1 },
 285   { 0, 8, 2, 14,5, 11,3, 12,0 }};
 286 static int QuadPenta_RE [5][9] = { // REVERSED -> EXTERNAL
 287   { 0, 8, 2, 7, 1, 6, 0, 0, 0 },
 288   { 3, 9, 4, 10,5, 11,3, 3, 3 },
 289   { 0, 6, 1, 13,4, 9, 3, 12,0 },
 290   { 1, 7, 2, 14,5, 10,4, 13,1 },
 291   { 0, 12,3, 11,5, 14,2, 8, 0 }};
 292 static int QuadPenta_nbN [] = { 6, 6, 8, 8, 8 };
 293
 294 /*
 295 //                 13
 296 //         N5+-----+-----+N6
 297 //          /|          /|
 298 //       12+ |       14+ |
 299 //        /  |        /  |
 300 //     N4+-----+-----+N7 |           QUADRATIC
 301 //       |   | 15    |   |           HEXAHEDRON
 302 //       |   |       |   |
 303 //       | 17+       |   +18
 304 //       |   |       |   |
 305 //       |   |       |   |
 306 //       |   |       |   |
 307 //     16+   |       +19 |
 308 //       |   |       |   |
 309 //       |   |     9 |   |
 310 //       | N1+-----+-|---+N2
 311 //       |  /        |  /
 312 //       | +8        | +10
 313 //       |/          |/
 314 //     N0+-----+-----+N3
 315 //             11
 316 */
 317 static int QuadHexa_F [6][9] = {  // FORWARD
 318   { 0, 8, 1, 9, 2, 10,3, 11,0 },   // all face normals are external,
 319   { 4, 15,7, 14,6, 13,5, 12,4 },
 320   { 0, 16,4, 12,5, 17,1, 8, 0 },
 321   { 1, 17,5, 13,6, 18,2, 9, 1 },
 322   { 3, 10,2, 18,6, 14,7, 19,3 },
 323   { 0, 11,3, 19,7, 15,4, 16,0 }};
 324 static int QuadHexa_RE [6][9] = {  // REVERSED -> EXTERNAL
 325   { 0, 11,3, 10,2, 9, 1, 8, 0 },   // all face normals are external
 326   { 4, 12,5, 13,6, 14,7, 15,4 },
 327   { 0, 8, 1, 17,5, 12,4, 16,0 },
 328   { 1, 9, 2, 18,6, 13,5, 17,1 },
 329   { 3, 19,7, 14,6, 18,2, 10,3 },
 330   { 0, 16,4, 15,7, 19,3, 11,0 }};
 331 static int QuadHexa_nbN [] = { 8, 8, 8, 8, 8, 8 };
 332
 333 static int TriQuadHexa_F [6][9] = {  // FORWARD
 334   { 0, 8, 1, 9, 2, 10,3, 11, 20 },   // all face normals are external
 335   { 4, 15,7, 14,6, 13,5, 12, 25 },
 336   { 0, 16,4, 12,5, 17,1, 8,  21 },
 337   { 1, 17,5, 13,6, 18,2, 9,  22 },
 338   { 3, 10,2, 18,6, 14,7, 19, 23 },
 339   { 0, 11,3, 19,7, 15,4, 16, 24 }};
 340 static int TriQuadHexa_RE [6][9] = {  // REVERSED -> EXTERNAL
 341   { 0, 11,3, 10,2, 9, 1, 8,  20 },   // opposite faces are neighbouring,
 342   { 4, 12,5, 13,6, 14,7, 15, 25 },   // all face normals are external
 343   { 0, 8, 1, 17,5, 12,4, 16, 21 },
 344   { 1, 9, 2, 18,6, 13,5, 17, 22 },
 345   { 3, 19,7, 14,6, 18,2, 10, 23 },
 346   { 0, 16,4, 15,7, 19,3, 11, 24 }};
 347 static int TriQuadHexa_nbN [] = { 9, 9, 9, 9, 9, 9 };
 348
 349
 350 // ========================================================
 351 // to perform some calculations without linkage to CASCADE
 352 // ========================================================
 353 namespace
 354 {
 355 struct XYZ {
 356   double x;
 357   double y;
 358   double z;
 359   XYZ()                               { x = 0; y = 0; z = 0; }
 360   XYZ( double X, double Y, double Z ) { x = X; y = Y; z = Z; }
 361   XYZ( const XYZ& other )             { x = other.x; y = other.y; z = other.z; }
 362   XYZ( const SMDS_MeshNode* n )       { x = n->X(); y = n->Y(); z = n->Z(); }
 363   inline XYZ operator-( const XYZ& other );
 364   inline XYZ operator+( const XYZ& other );
 365   inline XYZ Crossed( const XYZ& other );
 366   inline double Dot( const XYZ& other );
 367   inline double Magnitude();
 368   inline double SquareMagnitude();
 369 };
 370 inline XYZ XYZ::operator-( const XYZ& Right ) {
 371   return XYZ(x - Right.x, y - Right.y, z - Right.z);
 372 }
 373 inline XYZ XYZ::operator+( const XYZ& Right ) {
 374   return XYZ(x + Right.x, y + Right.y, z + Right.z);
 375 }
 376 inline XYZ XYZ::Crossed( const XYZ& Right ) {
 377   return XYZ (y * Right.z - z * Right.y,
 378               z * Right.x - x * Right.z,
 379               x * Right.y - y * Right.x);
 380 }
 381 inline double XYZ::Dot( const XYZ& Other ) {
 382   return(x * Other.x + y * Other.y + z * Other.z);
 383 }
 384 inline double XYZ::Magnitude() {
 385   return sqrt (x * x + y * y + z * z);
 386 }
 387 inline double XYZ::SquareMagnitude() {
 388   return (x * x + y * y + z * z);
 389 }
 390
 391   //================================================================================
 392   /*!
 393    * \brief Return linear type corresponding to a quadratic one
 394    */
 395   //================================================================================
 396
 397   SMDS_VolumeTool::VolumeType quadToLinear(SMDS_VolumeTool::VolumeType quadType)
 398   {
 399     SMDS_VolumeTool::VolumeType linType = SMDS_VolumeTool::VolumeType( int(quadType)-4 );
 400     const int nbCornersByQuad = SMDS_VolumeTool::NbCornerNodes( quadType );
 401     if ( SMDS_VolumeTool::NbCornerNodes( linType ) == nbCornersByQuad )
 402       return linType;
 403
 404     int iLin = 0;
 405     for ( ; iLin < SMDS_VolumeTool::NB_VOLUME_TYPES; ++iLin )
 406       if ( SMDS_VolumeTool::NbCornerNodes( SMDS_VolumeTool::VolumeType( iLin )) == nbCornersByQuad)
 407         return SMDS_VolumeTool::VolumeType( iLin );
 408
 409     return SMDS_VolumeTool::UNKNOWN;
 410   }
 411
 412 } // namespace
 413
 414 //=======================================================================
 415 //function : SMDS_VolumeTool
 416 //purpose  :
 417 //=======================================================================
 418
 419 SMDS_VolumeTool::SMDS_VolumeTool ()
 420   : myVolumeNodes( NULL ),
 421     myFaceNodes( NULL )
 422 {
 423   Set( 0 );
 424 }
 425
 426 //=======================================================================
 427 //function : SMDS_VolumeTool
 428 //purpose  :
 429 //=======================================================================
 430
 431 SMDS_VolumeTool::SMDS_VolumeTool (const SMDS_MeshElement* theVolume,
 432                                   const bool              ignoreCentralNodes)
 433   : myVolumeNodes( NULL ),
 434     myFaceNodes( NULL )
 435 {
 436   Set( theVolume, ignoreCentralNodes );
 437 }
 438
 439 //=======================================================================
 440 //function : SMDS_VolumeTool
 441 //purpose  :
 442 //=======================================================================
 443
 444 SMDS_VolumeTool::~SMDS_VolumeTool()
 445 {
 446   if ( myVolumeNodes != NULL ) delete [] myVolumeNodes;
 447   if ( myFaceNodes != NULL   ) delete [] myFaceNodes;
 448
 449   myFaceNodeIndices = NULL;
 450   myVolumeNodes = myFaceNodes = NULL;
 451 }
 452
 453 //=======================================================================
 454 //function : SetVolume
 455 //purpose  : Set volume to iterate on
 456 //=======================================================================
 457
 458 bool SMDS_VolumeTool::Set (const SMDS_MeshElement* theVolume,
 459                            const bool              ignoreCentralNodes)
 460 {
 461   // reset fields
 462   myVolume = 0;
 463   myPolyedre = 0;
 464   myIgnoreCentralNodes = ignoreCentralNodes;
 465
 466   myVolForward = true;
 467   myNbFaces = 0;
 468   myVolumeNbNodes = 0;
 469   if (myVolumeNodes != NULL) {
 470     delete [] myVolumeNodes;
 471     myVolumeNodes = NULL;
 472   }
 473   myPolyIndices.clear();
 474
 475   myExternalFaces = false;
 476
 477   myAllFacesNodeIndices_F  = 0;
 478   //myAllFacesNodeIndices_FE = 0;
 479   myAllFacesNodeIndices_RE = 0;
 480   myAllFacesNbNodes        = 0;
 481
 482   myCurFace = -1;
 483   myFaceNbNodes = 0;
 484   myFaceNodeIndices = NULL;
 485   if (myFaceNodes != NULL) {
 486     delete [] myFaceNodes;
 487     myFaceNodes = NULL;
 488   }
 489
 490   // set volume data
 491   if ( !theVolume || theVolume->GetType() != SMDSAbs_Volume )
 492     return false;
 493
 494   myVolume = theVolume;
 495   if (myVolume->IsPoly())
 496     myPolyedre = dynamic_cast<const SMDS_VtkVolume*>( myVolume );
 497
 498   myNbFaces = theVolume->NbFaces();
 499   myVolumeNbNodes = theVolume->NbNodes();
 500
 501   // set nodes
 502   int iNode = 0;
 503   myVolumeNodes = new const SMDS_MeshNode* [myVolumeNbNodes];
 504   SMDS_ElemIteratorPtr nodeIt = myVolume->nodesIterator();
 505   while ( nodeIt->more() )
 506     myVolumeNodes[ iNode++ ] = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
 507
 508   // check validity
 509   if ( !setFace(0) )
 510     return ( myVolume = 0 );
 511
 512   if ( !myPolyedre )
 513   {
 514     // define volume orientation
 515     XYZ botNormal;
 516     GetFaceNormal( 0, botNormal.x, botNormal.y, botNormal.z );
 517     const SMDS_MeshNode* botNode = myVolumeNodes[ 0 ];
 518     int topNodeIndex = myVolume->NbCornerNodes() - 1;
 519     while ( !IsLinked( 0, topNodeIndex, /*ignoreMediumNodes=*/true )) --topNodeIndex;
 520     const SMDS_MeshNode* topNode = myVolumeNodes[ topNodeIndex ];
 521     XYZ upDir (topNode->X() - botNode->X(),
 522                topNode->Y() - botNode->Y(),
 523                topNode->Z() - botNode->Z() );
 524     myVolForward = ( botNormal.Dot( upDir ) < 0 );
 525   }
 526   return true;
 527 }
 528
 529 //=======================================================================
 530 //function : Inverse
 531 //purpose  : Inverse volume
 532 //=======================================================================
 533
 534 #define SWAP_NODES(nodes,i1,i2)           \
 535 {                                         \
 536   const SMDS_MeshNode* tmp = nodes[ i1 ]; \
 537   nodes[ i1 ] = nodes[ i2 ];              \
 538   nodes[ i2 ] = tmp;                      \
 539 }
 540 void SMDS_VolumeTool::Inverse ()
 541 {
 542   if ( !myVolume ) return;
 543
 544   if (myVolume->IsPoly()) {
 545     MESSAGE("Warning: attempt to inverse polyhedral volume");
 546     return;
 547   }
 548
 549   myVolForward = !myVolForward;
 550   myCurFace = -1;
 551
 552   // inverse top and bottom faces
 553   switch ( myVolumeNbNodes ) {
 554   case 4:
 555     SWAP_NODES( myVolumeNodes, 1, 2 );
 556     break;
 557   case 5:
 558     SWAP_NODES( myVolumeNodes, 1, 3 );
 559     break;
 560   case 6:
 561     SWAP_NODES( myVolumeNodes, 1, 2 );
 562     SWAP_NODES( myVolumeNodes, 4, 5 );
 563     break;
 564   case 8:
 565     SWAP_NODES( myVolumeNodes, 1, 3 );
 566     SWAP_NODES( myVolumeNodes, 5, 7 );
 567     break;
 568   case 12:
 569     SWAP_NODES( myVolumeNodes, 1, 5 );
 570     SWAP_NODES( myVolumeNodes, 2, 4 );
 571     SWAP_NODES( myVolumeNodes, 7, 11 );
 572     SWAP_NODES( myVolumeNodes, 8, 10 );
 573     break;
 574
 575   case 10:
 576     SWAP_NODES( myVolumeNodes, 1, 2 );
 577     SWAP_NODES( myVolumeNodes, 4, 6 );
 578     SWAP_NODES( myVolumeNodes, 8, 9 );
 579     break;
 580   case 13:
 581     SWAP_NODES( myVolumeNodes, 1, 3 );
 582     SWAP_NODES( myVolumeNodes, 5, 8 );
 583     SWAP_NODES( myVolumeNodes, 6, 7 );
 584     SWAP_NODES( myVolumeNodes, 10, 12 );
 585     break;
 586   case 15:
 587     SWAP_NODES( myVolumeNodes, 1, 2 );
 588     SWAP_NODES( myVolumeNodes, 4, 5 );
 589     SWAP_NODES( myVolumeNodes, 6, 8 );
 590     SWAP_NODES( myVolumeNodes, 9, 11 );
 591     SWAP_NODES( myVolumeNodes, 13, 14 );
 592     break;
 593   case 20:
 594     SWAP_NODES( myVolumeNodes, 1, 3 );
 595     SWAP_NODES( myVolumeNodes, 5, 7 );
 596     SWAP_NODES( myVolumeNodes, 8, 11 );
 597     SWAP_NODES( myVolumeNodes, 9, 10 );
 598     SWAP_NODES( myVolumeNodes, 12, 15 );
 599     SWAP_NODES( myVolumeNodes, 13, 14 );
 600     SWAP_NODES( myVolumeNodes, 17, 19 );
 601     break;
 602   case 27:
 603     SWAP_NODES( myVolumeNodes, 1, 3 );
 604     SWAP_NODES( myVolumeNodes, 5, 7 );
 605     SWAP_NODES( myVolumeNodes, 8, 11 );
 606     SWAP_NODES( myVolumeNodes, 9, 10 );
 607     SWAP_NODES( myVolumeNodes, 12, 15 );
 608     SWAP_NODES( myVolumeNodes, 13, 14 );
 609     SWAP_NODES( myVolumeNodes, 17, 19 );
 610     SWAP_NODES( myVolumeNodes, 21, 24 );
 611     SWAP_NODES( myVolumeNodes, 22, 23 );
 612     break;
 613   default:;
 614   }
 615 }
 616
 617 //=======================================================================
 618 //function : GetVolumeType
 619 //purpose  :
 620 //=======================================================================
 621
 622 SMDS_VolumeTool::VolumeType SMDS_VolumeTool::GetVolumeType() const
 623 {
 624   if ( myPolyedre )
 625     return POLYHEDA;
 626
 627   switch( myVolumeNbNodes ) {
 628   case 4: return TETRA;
 629   case 5: return PYRAM;
 630   case 6: return PENTA;
 631   case 8: return HEXA;
 632   case 12: return HEX_PRISM;
 633   case 10: return QUAD_TETRA;
 634   case 13: return QUAD_PYRAM;
 635   case 15: return QUAD_PENTA;
 636   case 20: return QUAD_HEXA;
 637   case 27: return QUAD_HEXA;
 638   default: break;
 639   }
 640
 641   return UNKNOWN;
 642 }
 643
 644 //=======================================================================
 645 //function : getTetraVolume
 646 //purpose  :
 647 //=======================================================================
 648
 649 static double getTetraVolume(const SMDS_MeshNode* n1,
 650                              const SMDS_MeshNode* n2,
 651                              const SMDS_MeshNode* n3,
 652                              const SMDS_MeshNode* n4)
 653 {
 654   double X1 = n1->X();
 655   double Y1 = n1->Y();
 656   double Z1 = n1->Z();
 657
 658   double X2 = n2->X();
 659   double Y2 = n2->Y();
 660   double Z2 = n2->Z();
 661
 662   double X3 = n3->X();
 663   double Y3 = n3->Y();
 664   double Z3 = n3->Z();
 665
 666   double X4 = n4->X();
 667   double Y4 = n4->Y();
 668   double Z4 = n4->Z();
 669
 670   double Q1 = -(X1-X2)*(Y3*Z4-Y4*Z3);
 671   double Q2 =  (X1-X3)*(Y2*Z4-Y4*Z2);
 672   double R1 = -(X1-X4)*(Y2*Z3-Y3*Z2);
 673   double R2 = -(X2-X3)*(Y1*Z4-Y4*Z1);
 674   double S1 =  (X2-X4)*(Y1*Z3-Y3*Z1);
 675   double S2 = -(X3-X4)*(Y1*Z2-Y2*Z1);
 676
 677   return (Q1+Q2+R1+R2+S1+S2)/6.0;
 678 }
 679
 680 //=======================================================================
 681 //function : GetSize
 682 //purpose  : Return element volume
 683 //=======================================================================
 684
 685 double SMDS_VolumeTool::GetSize() const
 686 {
 687   double V = 0.;
 688   if ( !myVolume )
 689     return 0.;
 690
 691   if ( myVolume->IsPoly() )
 692   {
 693     if ( !myPolyedre )
 694       return 0.;
 695
 696     // split a polyhedron into tetrahedrons
 697
 698     int saveCurFace = myCurFace;
 699     SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* > ( this );
 700     for ( int f = 0; f < NbFaces(); ++f )
 701     {
 702       me->setFace( f );
 703       XYZ area (0,0,0), p1( myFaceNodes[0] );
 704       for ( int n = 0; n < myFaceNbNodes; ++n )
 705       {
 706         XYZ p2( myFaceNodes[ n+1 ]);
 707         area = area + p1.Crossed( p2 );
 708         p1 = p2;
 709       }
 710       V += p1.Dot( area );
 711     }
 712     V /= 6;
 713     if ( saveCurFace > -1 && saveCurFace != myCurFace )
 714       me->setFace( myCurFace );
 715   }
 716   else
 717   {
 718     const static int ind[] = {
 719       0, 1, 3, 6, 11, 23, 31, 44, 58, 78 };
 720     const static int vtab[][4] = { // decomposition into tetra in the order of enum VolumeType
 721       // tetrahedron
 722       { 0, 1, 2, 3 },
 723       // pyramid
 724       { 0, 1, 3, 4 },
 725       { 1, 2, 3, 4 },
 726       // pentahedron
 727       { 0, 1, 2, 3 },
 728       { 1, 5, 3, 4 },
 729       { 1, 5, 2, 3 },
 730       // hexahedron
 731       { 1, 4, 3, 0 },
 732       { 4, 1, 6, 5 },
 733       { 1, 3, 6, 2 },
 734       { 4, 6, 3, 7 },
 735       { 1, 4, 6, 3 },
 736       // hexagonal prism
 737       { 0, 1, 2, 7 },
 738       { 0, 7, 8, 6 },
 739       { 2, 7, 8, 0 },
 740
 741       { 0, 3, 4, 9 },
 742       { 0, 9, 10, 6 },
 743       { 4, 9, 10, 0 },
 744
 745       { 0, 3, 4, 9 },
 746       { 0, 9, 10, 6 },
 747       { 4, 9, 10, 0 },
 748
 749       { 0, 4, 5, 10 },
 750       { 0, 10, 11, 6 },
 751       { 5, 10, 11, 0 },
 752
 753       // quadratic tetrahedron
 754       { 0, 4, 6, 7 },
 755       { 1, 5, 4, 8 },
 756       { 2, 6, 5, 9 },
 757       { 7, 8, 9, 3 },
 758       { 4, 6, 7, 9 },
 759       { 4, 5, 6, 9 },
 760       { 4, 7, 8, 9 },
 761       { 4, 5, 9, 8 },
 762
 763       // quadratic pyramid
 764       { 0, 5, 8, 9 },
 765       { 1, 5,10, 6 },
 766       { 2, 6,11, 7 },
 767       { 3, 7,12, 8 },
 768       { 4, 9,11,10 },
 769       { 4, 9,12,11 },
 770       { 10, 5, 9, 8 },
 771       { 10, 8, 9,12 },
 772       { 10, 8,12, 7 },
 773       { 10, 7,12,11 },
 774       { 10, 7,11, 6 },
 775       { 10, 5, 8, 6 },
 776       { 10, 6, 8, 7 },
 777
 778       // quadratic pentahedron
 779       { 12, 0, 8, 6 },
 780       { 12, 8, 7, 6 },
 781       { 12, 8, 2, 7 },
 782       { 12, 6, 7, 1 },
 783       { 12, 1, 7,13 },
 784       { 12, 7, 2,13 },
 785       { 12, 2,14,13 },
 786
 787       { 12, 3, 9,11 },
 788       { 12,11, 9,10 },
 789       { 12,11,10, 5 },
 790       { 12, 9, 4,10 },
 791       { 12,14, 5,10 },
 792       { 12,14,10, 4 },
 793       { 12,14, 4,13 },
 794
 795       // quadratic hexahedron
 796       { 16, 0,11, 8 },
 797       { 16,11, 9, 8 },
 798       { 16, 8, 9, 1 },
 799       { 16,11, 3,10 },
 800       { 16,11,10, 9 },
 801       { 16,10, 2, 9 },
 802       { 16, 3,19, 2 },
 803       { 16, 2,19,18 },
 804       { 16, 2,18,17 },
 805       { 16, 2,17, 1 },
 806
 807       { 16, 4,12,15 },
 808       { 16,12, 5,13 },
 809       { 16,12,13,15 },
 810       { 16,13, 6,14 },
 811       { 16,13,14,15 },
 812       { 16,14, 7,15 },
 813       { 16, 6, 5,17 },
 814       { 16,18, 6,17 },
 815       { 16,18, 7, 6 },
 816       { 16,18,19, 7 },
 817
 818     };
 819
 820     int type = GetVolumeType();
 821     int n1 = ind[type];
 822     int n2 = ind[type+1];
 823
 824     for (int i = n1; i <  n2; i++) {
 825       V -= getTetraVolume( myVolumeNodes[ vtab[i][0] ],
 826                            myVolumeNodes[ vtab[i][1] ],
 827                            myVolumeNodes[ vtab[i][2] ],
 828                            myVolumeNodes[ vtab[i][3] ]);
 829     }
 830   }
 831   return V;
 832 }
 833
 834 //=======================================================================
 835 //function : GetBaryCenter
 836 //purpose  :
 837 //=======================================================================
 838
 839 bool SMDS_VolumeTool::GetBaryCenter(double & X, double & Y, double & Z) const
 840 {
 841   X = Y = Z = 0.;
 842   if ( !myVolume )
 843     return false;
 844
 845   for ( int i = 0; i < myVolumeNbNodes; i++ ) {
 846     X += myVolumeNodes[ i ]->X();
 847     Y += myVolumeNodes[ i ]->Y();
 848     Z += myVolumeNodes[ i ]->Z();
 849   }
 850   X /= myVolumeNbNodes;
 851   Y /= myVolumeNbNodes;
 852   Z /= myVolumeNbNodes;
 853
 854   return true;
 855 }
 856
 857 //================================================================================
 858 /*!
 859  * \brief Classify a point
 860  *  \param tol - thickness of faces
 861  */
 862 //================================================================================
 863
 864 bool SMDS_VolumeTool::IsOut(double X, double Y, double Z, double tol) const
 865 {
 866   // LIMITATION: for convex volumes only
 867   XYZ p( X,Y,Z );
 868   for ( int iF = 0; iF < myNbFaces; ++iF )
 869   {
 870     XYZ faceNormal;
 871     if ( !GetFaceNormal( iF, faceNormal.x, faceNormal.y, faceNormal.z ))
 872       continue;
 873     if ( !IsFaceExternal( iF ))
 874       faceNormal = XYZ() - faceNormal; // reverse
 875
 876     XYZ face2p( p - XYZ( myFaceNodes[0] ));
 877     if ( face2p.Dot( faceNormal ) > tol )
 878       return true;
 879   }
 880   return false;
 881 }
 882
 883 //=======================================================================
 884 //function : SetExternalNormal
 885 //purpose  : Node order will be so that faces normals are external
 886 //=======================================================================
 887
 888 void SMDS_VolumeTool::SetExternalNormal ()
 889 {
 890   myExternalFaces = true;
 891   myCurFace = -1;
 892 }
 893
 894 //=======================================================================
 895 //function : NbFaceNodes
 896 //purpose  : Return number of nodes in the array of face nodes
 897 //=======================================================================
 898
 899 int SMDS_VolumeTool::NbFaceNodes( int faceIndex ) const
 900 {
 901     if ( !setFace( faceIndex ))
 902       return 0;
 903     return myFaceNbNodes;
 904 }
 905
 906 //=======================================================================
 907 //function : GetFaceNodes
 908 //purpose  : Return pointer to the array of face nodes.
 909 //           To comfort link iteration, the array
 910 //           length == NbFaceNodes( faceIndex ) + 1 and
 911 //           the last node == the first one.
 912 //=======================================================================
 913
 914 const SMDS_MeshNode** SMDS_VolumeTool::GetFaceNodes( int faceIndex ) const
 915 {
 916   if ( !setFace( faceIndex ))
 917     return 0;
 918   return myFaceNodes;
 919 }
 920
 921 //=======================================================================
 922 //function : GetFaceNodesIndices
 923 //purpose  : Return pointer to the array of face nodes indices
 924 //           To comfort link iteration, the array
 925 //           length == NbFaceNodes( faceIndex ) + 1 and
 926 //           the last node index == the first one.
 927 //=======================================================================
 928
 929 const int* SMDS_VolumeTool::GetFaceNodesIndices( int faceIndex ) const
 930 {
 931   if ( !setFace( faceIndex ))
 932     return 0;
 933
 934   if (myPolyedre)
 935   {
 936     SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* > ( this );
 937     me->myPolyIndices.resize( myFaceNbNodes + 1 );
 938     me->myFaceNodeIndices = & me->myPolyIndices[0];
 939     for ( int i = 0; i <= myFaceNbNodes; ++i )
 940       me->myFaceNodeIndices[i] = myVolume->GetNodeIndex( myFaceNodes[i] );
 941   }
 942   return myFaceNodeIndices;
 943 }
 944
 945 //=======================================================================
 946 //function : GetFaceNodes
 947 //purpose  : Return a set of face nodes.
 948 //=======================================================================
 949
 950 bool SMDS_VolumeTool::GetFaceNodes (int                        faceIndex,
 951                                     set<const SMDS_MeshNode*>& theFaceNodes ) const
 952 {
 953   if ( !setFace( faceIndex ))
 954     return false;
 955
 956   theFaceNodes.clear();
 957   theFaceNodes.insert( myFaceNodes, myFaceNodes + myFaceNbNodes );
 958
 959   return true;
 960 }
 961
 962 //=======================================================================
 963 //function : IsFaceExternal
 964 //purpose  : Check normal orientation of a given face
 965 //=======================================================================
 966
 967 bool SMDS_VolumeTool::IsFaceExternal( int faceIndex ) const
 968 {
 969   if ( myExternalFaces || !myVolume )
 970     return true;
 971
 972   if (myVolume->IsPoly()) {
 973     XYZ aNormal, baryCenter, p0 (myPolyedre->GetFaceNode(faceIndex + 1, 1));
 974     GetFaceNormal(faceIndex, aNormal.x, aNormal.y, aNormal.z);
 975     GetBaryCenter(baryCenter.x, baryCenter.y, baryCenter.z);
 976     XYZ insideVec (baryCenter - p0);
 977     if (insideVec.Dot(aNormal) > 0)
 978       return false;
 979     return true;
 980   }
 981
 982   // switch ( myVolumeNbNodes ) {
 983   // case 4:
 984   // case 5:
 985   // case 10:
 986   // case 13:
 987   //   // only the bottom of a reversed tetrahedron can be internal
 988   //   return ( myVolForward || faceIndex != 0 );
 989   // case 6:
 990   // case 15:
 991   // case 12:
 992   //   // in a forward prism, the top is internal, in a reversed one - bottom
 993   //   return ( myVolForward ? faceIndex != 1 : faceIndex != 0 );
 994   // case 8:
 995   // case 20:
 996   // case 27: {
 997   //   // in a forward hexahedron, even face normal is external, odd - internal
 998   //   bool odd = faceIndex % 2;
 999   //   return ( myVolForward ? !odd : odd );
1000   // }
1001   // default:;
1002   // }
1003   // return false;
1004   return true;
1005 }
1006
1007 //=======================================================================
1008 //function : GetFaceNormal
1009 //purpose  : Return a normal to a face
1010 //=======================================================================
1011
1012 bool SMDS_VolumeTool::GetFaceNormal (int faceIndex, double & X, double & Y, double & Z) const
1013 {
1014   if ( !setFace( faceIndex ))
1015     return false;
1016
1017   XYZ p1 ( myFaceNodes[0] );
1018   XYZ p2 ( myFaceNodes[1] );
1019   XYZ p3 ( myFaceNodes[2] );
1020   XYZ aVec12( p2 - p1 );
1021   XYZ aVec13( p3 - p1 );
1022   XYZ cross = aVec12.Crossed( aVec13 );
1023
1024   //if ( myFaceNbNodes == 4 ) {
1025   if ( myFaceNbNodes >3 ) {
1026     XYZ p4 ( myFaceNodes[3] );
1027     XYZ aVec14( p4 - p1 );
1028     XYZ cross2 = aVec13.Crossed( aVec14 );
1029     cross = cross + cross2;
1030   }
1031
1032   double size = cross.Magnitude();
1033   if ( size <= numeric_limits<double>::min() )
1034     return false;
1035
1036   X = cross.x / size;
1037   Y = cross.y / size;
1038   Z = cross.z / size;
1039
1040   return true;
1041 }
1042
1043 //================================================================================
1044 /*!
1045  * \brief Return barycenter of a face
1046  */
1047 //================================================================================
1048
1049 bool SMDS_VolumeTool::GetFaceBaryCenter (int faceIndex, double & X, double & Y, double & Z) const
1050 {
1051   if ( !setFace( faceIndex ))
1052     return false;
1053
1054   X = Y = Z = 0.0;
1055   for ( int i = 0; i < myFaceNbNodes; ++i )
1056   {
1057     X += myFaceNodes[i]->X() / myFaceNbNodes;
1058     Y += myFaceNodes[i]->Y() / myFaceNbNodes;
1059     Z += myFaceNodes[i]->Z() / myFaceNbNodes;
1060   }
1061   return true;
1062 }
1063
1064 //=======================================================================
1065 //function : GetFaceArea
1066 //purpose  : Return face area
1067 //=======================================================================
1068
1069 double SMDS_VolumeTool::GetFaceArea( int faceIndex ) const
1070 {
1071   if (myVolume->IsPoly()) {
1072     MESSAGE("Warning: attempt to obtain area of a face of polyhedral volume");
1073     return 0;
1074   }
1075
1076   if ( !setFace( faceIndex ))
1077     return 0;
1078
1079   XYZ p1 ( myFaceNodes[0] );
1080   XYZ p2 ( myFaceNodes[1] );
1081   XYZ p3 ( myFaceNodes[2] );
1082   XYZ aVec12( p2 - p1 );
1083   XYZ aVec13( p3 - p1 );
1084   double area = aVec12.Crossed( aVec13 ).Magnitude() * 0.5;
1085
1086   if ( myFaceNbNodes == 4 ) {
1087     XYZ p4 ( myFaceNodes[3] );
1088     XYZ aVec14( p4 - p1 );
1089     area += aVec14.Crossed( aVec13 ).Magnitude() * 0.5;
1090   }
1091   return area;
1092 }
1093
1094 //================================================================================
1095 /*!
1096  * \brief Return index of the node located at face center of a quadratic element like HEX27
1097  */
1098 //================================================================================
1099
1100 int SMDS_VolumeTool::GetCenterNodeIndex( int faceIndex ) const
1101 {
1102   if ( myAllFacesNbNodes && myVolumeNbNodes == 27 ) // classic element with 27 nodes
1103   {
1104     switch ( faceIndex ) {
1105     case 0: return 20;
1106     case 1: return 25;
1107     default:
1108       return faceIndex + 19;
1109     }
1110   }
1111   return -1;
1112 }
1113
1114 //=======================================================================
1115 //function : GetOppFaceIndex
1116 //purpose  : Return index of the opposite face if it exists, else -1.
1117 //=======================================================================
1118
1119 int SMDS_VolumeTool::GetOppFaceIndex( int faceIndex ) const
1120 {
1121   int ind = -1;
1122   if (myPolyedre) {
1123     MESSAGE("Warning: attempt to obtain opposite face on polyhedral volume");
1124     return ind;
1125   }
1126
1127   const int nbHoriFaces = 2;
1128
1129   if ( faceIndex >= 0 && faceIndex < NbFaces() ) {
1130     switch ( myVolumeNbNodes ) {
1131     case 6:
1132     case 15:
1133       if ( faceIndex == 0 || faceIndex == 1 )
1134         ind = 1 - faceIndex;
1135         break;
1136     case 8:
1137     case 12:
1138       if ( faceIndex <= 1 ) // top or bottom
1139         ind = 1 - faceIndex;
1140       else {
1141         const int nbSideFaces = myAllFacesNbNodes[0];
1142         ind = ( faceIndex - nbHoriFaces + nbSideFaces/2 ) % nbSideFaces + nbHoriFaces;
1143       }
1144       break;
1145     case 20:
1146     case 27:
1147       if ( faceIndex <= 1 ) // top or bottom
1148         ind = 1 - faceIndex;
1149       else {
1150         const int nbSideFaces = myAllFacesNbNodes[0] / 2;
1151         ind = ( faceIndex - nbHoriFaces + nbSideFaces/2 ) % nbSideFaces + nbHoriFaces;
1152       }
1153       break;
1154     default:;
1155     }
1156   }
1157   return ind;
1158 }
1159
1160 //=======================================================================
1161 //function : IsLinked
1162 //purpose  : return true if theNode1 is linked with theNode2
1163 // If theIgnoreMediumNodes then corner nodes of quadratic cell are considered linked as well
1164 //=======================================================================
1165
1166 bool SMDS_VolumeTool::IsLinked (const SMDS_MeshNode* theNode1,
1167                                 const SMDS_MeshNode* theNode2,
1168                                 const bool           theIgnoreMediumNodes) const
1169 {
1170   if ( !myVolume )
1171     return false;
1172
1173   if (myVolume->IsPoly()) {
1174     if (!myPolyedre) {
1175       MESSAGE("Warning: bad volumic element");
1176       return false;
1177     }
1178     bool isLinked = false;
1179     int iface;
1180     for (iface = 1; iface <= myNbFaces && !isLinked; iface++) {
1181       int inode, nbFaceNodes = myPolyedre->NbFaceNodes(iface);
1182
1183       for (inode = 1; inode <= nbFaceNodes && !isLinked; inode++) {
1184         const SMDS_MeshNode* curNode = myPolyedre->GetFaceNode(iface, inode);
1185
1186         if (curNode == theNode1 || curNode == theNode2) {
1187           int inextnode = (inode == nbFaceNodes) ? 1 : inode + 1;
1188           const SMDS_MeshNode* nextNode = myPolyedre->GetFaceNode(iface, inextnode);
1189
1190           if ((curNode == theNode1 && nextNode == theNode2) ||
1191               (curNode == theNode2 && nextNode == theNode1)) {
1192             isLinked = true;
1193           }
1194         }
1195       }
1196     }
1197     return isLinked;
1198   }
1199
1200   // find nodes indices
1201   int i1 = -1, i2 = -1, nbFound = 0;
1202   for ( int i = 0; i < myVolumeNbNodes && nbFound < 2; i++ )
1203   {
1204     if ( myVolumeNodes[ i ] == theNode1 )
1205       i1 = i, ++nbFound;
1206     else if ( myVolumeNodes[ i ] == theNode2 )
1207       i2 = i, ++nbFound;
1208   }
1209   return IsLinked( i1, i2 );
1210 }
1211
1212 //=======================================================================
1213 //function : IsLinked
1214 //purpose  : return true if the node with theNode1Index is linked
1215 //           with the node with theNode2Index
1216 // If theIgnoreMediumNodes then corner nodes of quadratic cell are considered linked as well
1217 //=======================================================================
1218
1219 bool SMDS_VolumeTool::IsLinked (const int theNode1Index,
1220                                 const int theNode2Index,
1221                                 bool      theIgnoreMediumNodes) const
1222 {
1223   if ( myVolume->IsPoly() ) {
1224     return IsLinked(myVolumeNodes[theNode1Index], myVolumeNodes[theNode2Index]);
1225   }
1226
1227   int minInd = min( theNode1Index, theNode2Index );
1228   int maxInd = max( theNode1Index, theNode2Index );
1229
1230   if ( minInd < 0 || maxInd > myVolumeNbNodes - 1 || maxInd == minInd )
1231     return false;
1232
1233   VolumeType type = GetVolumeType();
1234   if ( myVolume->IsQuadratic() )
1235   {
1236     int firstMediumInd = myVolume->NbCornerNodes();
1237     if ( minInd >= firstMediumInd )
1238       return false; // both nodes are medium - not linked
1239     if ( maxInd < firstMediumInd ) // both nodes are corners
1240     {
1241       if ( theIgnoreMediumNodes )
1242         type = quadToLinear(type); // to check linkage of corner nodes only
1243       else
1244         return false; // corner nodes are not linked directly in a quadratic cell
1245     }
1246   }
1247
1248   switch ( type ) {
1249   case TETRA:
1250     return true;
1251   case HEXA:
1252     switch ( maxInd - minInd ) {
1253     case 1: return minInd != 3;
1254     case 3: return minInd == 0 || minInd == 4;
1255     case 4: return true;
1256     default:;
1257     }
1258     break;
1259   case PYRAM:
1260     if ( maxInd == 4 )
1261       return true;
1262     switch ( maxInd - minInd ) {
1263     case 1:
1264     case 3: return true;
1265     default:;
1266     }
1267     break;
1268   case PENTA:
1269     switch ( maxInd - minInd ) {
1270     case 1: return minInd != 2;
1271     case 2: return minInd == 0 || minInd == 3;
1272     case 3: return true;
1273     default:;
1274     }
1275     break;
1276   case QUAD_TETRA:
1277     {
1278       switch ( minInd ) {
1279       case 0: if( maxInd==4 ||  maxInd==6 ||  maxInd==7 ) return true;
1280       case 1: if( maxInd==4 ||  maxInd==5 ||  maxInd==8 ) return true;
1281       case 2: if( maxInd==5 ||  maxInd==6 ||  maxInd==9 ) return true;
1282       case 3: if( maxInd==7 ||  maxInd==8 ||  maxInd==9 ) return true;
1283       default:;
1284       }
1285       break;
1286     }
1287   case QUAD_HEXA:
1288     {
1289       switch ( minInd ) {
1290       case 0: if( maxInd==8 ||  maxInd==11 ||  maxInd==16 ) return true;
1291       case 1: if( maxInd==8 ||  maxInd==9 ||  maxInd==17 ) return true;
1292       case 2: if( maxInd==9 ||  maxInd==10 ||  maxInd==18 ) return true;
1293       case 3: if( maxInd==10 ||  maxInd==11 ||  maxInd==19 ) return true;
1294       case 4: if( maxInd==12 ||  maxInd==15 ||  maxInd==16 ) return true;
1295       case 5: if( maxInd==12 ||  maxInd==13 ||  maxInd==17 ) return true;
1296       case 6: if( maxInd==13 ||  maxInd==14 ||  maxInd==18 ) return true;
1297       case 7: if( maxInd==14 ||  maxInd==15 ||  maxInd==19 ) return true;
1298       default:;
1299       }
1300       break;
1301     }
1302   case QUAD_PYRAM:
1303     {
1304       switch ( minInd ) {
1305       case 0: if( maxInd==5 ||  maxInd==8 ||  maxInd==9 ) return true;
1306       case 1: if( maxInd==5 ||  maxInd==6 ||  maxInd==10 ) return true;
1307       case 2: if( maxInd==6 ||  maxInd==7 ||  maxInd==11 ) return true;
1308       case 3: if( maxInd==7 ||  maxInd==8 ||  maxInd==12 ) return true;
1309       case 4: if( maxInd==9 ||  maxInd==10 ||  maxInd==11 ||  maxInd==12 ) return true;
1310       default:;
1311       }
1312       break;
1313     }
1314   case QUAD_PENTA:
1315     {
1316       switch ( minInd ) {
1317       case 0: if( maxInd==6 ||  maxInd==8 ||  maxInd==12 ) return true;
1318       case 1: if( maxInd==6 ||  maxInd==7 ||  maxInd==13 ) return true;
1319       case 2: if( maxInd==7 ||  maxInd==8 ||  maxInd==14 ) return true;
1320       case 3: if( maxInd==9 ||  maxInd==11 ||  maxInd==12 ) return true;
1321       case 4: if( maxInd==9 ||  maxInd==10 ||  maxInd==13 ) return true;
1322       case 5: if( maxInd==10 ||  maxInd==11 ||  maxInd==14 ) return true;
1323       default:;
1324       }
1325       break;
1326     }
1327   case HEX_PRISM:
1328     {
1329       const int diff = maxInd-minInd;
1330       if ( diff > 6  ) return false;// not linked top and bottom
1331       if ( diff == 6 ) return true; // linked top and bottom
1332       return diff == 1 || diff == 7;
1333     }
1334   default:;
1335   }
1336   return false;
1337 }
1338
1339 //=======================================================================
1340 //function : GetNodeIndex
1341 //purpose  : Return an index of theNode
1342 //=======================================================================
1343
1344 int SMDS_VolumeTool::GetNodeIndex(const SMDS_MeshNode* theNode) const
1345 {
1346   if ( myVolume ) {
1347     for ( int i = 0; i < myVolumeNbNodes; i++ ) {
1348       if ( myVolumeNodes[ i ] == theNode )
1349         return i;
1350     }
1351   }
1352   return -1;
1353 }
1354
1355 //================================================================================
1356 /*!
1357  * \brief Fill vector with boundary faces existing in the mesh
1358   * \param faces - vector of found nodes
1359   * \retval int - nb of found faces
1360  */
1361 //================================================================================
1362
1363 int SMDS_VolumeTool::GetAllExistingFaces(vector<const SMDS_MeshElement*> & faces) const
1364 {
1365   faces.clear();
1366   for ( int iF = 0; iF < NbFaces(); ++iF ) {
1367     const SMDS_MeshFace* face = 0;
1368     const SMDS_MeshNode** nodes = GetFaceNodes( iF );
1369     switch ( NbFaceNodes( iF )) {
1370     case 3:
1371       face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2] ); break;
1372     case 4:
1373       face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2], nodes[3] ); break;
1374     case 6:
1375       face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2],
1376                                   nodes[3], nodes[4], nodes[5]); break;
1377     case 8:
1378       face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2], nodes[3],
1379                                   nodes[4], nodes[5], nodes[6], nodes[7]); break;
1380     }
1381     if ( face )
1382       faces.push_back( face );
1383   }
1384   return faces.size();
1385 }
1386
1387
1388 //================================================================================
1389 /*!
1390  * \brief Fill vector with boundary edges existing in the mesh
1391   * \param edges - vector of found edges
1392   * \retval int - nb of found faces
1393  */
1394 //================================================================================
1395
1396 int SMDS_VolumeTool::GetAllExistingEdges(vector<const SMDS_MeshElement*> & edges) const
1397 {
1398   edges.clear();
1399   edges.reserve( myVolumeNbNodes * 2 );
1400   for ( int i = 0; i < myVolumeNbNodes-1; ++i ) {
1401     for ( int j = i + 1; j < myVolumeNbNodes; ++j ) {
1402       if ( IsLinked( i, j )) {
1403         const SMDS_MeshElement* edge =
1404           SMDS_Mesh::FindEdge( myVolumeNodes[i], myVolumeNodes[j] );
1405         if ( edge )
1406           edges.push_back( edge );
1407       }
1408     }
1409   }
1410   return edges.size();
1411 }
1412
1413 //================================================================================
1414 /*!
1415  * \brief Return minimal square distance between connected corner nodes
1416  */
1417 //================================================================================
1418
1419 double SMDS_VolumeTool::MinLinearSize2() const
1420 {
1421   double minSize = 1e+100;
1422   int iQ = myVolume->IsQuadratic() ? 2 : 1;
1423
1424   // store current face data
1425   int curFace = myCurFace, nbN = myFaceNbNodes;
1426   int* ind = myFaceNodeIndices;
1427   myFaceNodeIndices = NULL;
1428   const SMDS_MeshNode** nodes = myFaceNodes;
1429   myFaceNodes = NULL;
1430
1431   // it seems that compute distance twice is faster than organization of a sole computing
1432   myCurFace = -1;
1433   for ( int iF = 0; iF < myNbFaces; ++iF )
1434   {
1435     setFace( iF );
1436     for ( int iN = 0; iN < myFaceNbNodes; iN += iQ )
1437     {
1438       XYZ n1( myFaceNodes[ iN ]);
1439       XYZ n2( myFaceNodes[(iN + iQ) % myFaceNbNodes]);
1440       minSize = std::min( minSize, (n1 - n2).SquareMagnitude());
1441     }
1442   }
1443   // restore current face data
1444   myCurFace = curFace;
1445   myFaceNbNodes = nbN;
1446   myFaceNodeIndices = ind;
1447   delete [] myFaceNodes; myFaceNodes = nodes;
1448
1449   return minSize;
1450 }
1451
1452 //================================================================================
1453 /*!
1454  * \brief check that only one volume is build on the face nodes
1455  *
1456  * If a face is shared by one of <ignoreVolumes>, it is considered free
1457  */
1458 //================================================================================
1459
1460 bool SMDS_VolumeTool::IsFreeFace( int faceIndex, const SMDS_MeshElement** otherVol/*=0*/ ) const
1461 {
1462   const bool isFree = true;
1463
1464   if (!setFace( faceIndex ))
1465     return !isFree;
1466
1467   const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
1468   const int nbFaceNodes = myFaceNbNodes;
1469
1470   // evaluate nb of face nodes shared by other volumes
1471   int maxNbShared = -1;
1472   typedef map< const SMDS_MeshElement*, int > TElemIntMap;
1473   TElemIntMap volNbShared;
1474   TElemIntMap::iterator vNbIt;
1475   for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) {
1476     const SMDS_MeshNode* n = nodes[ iNode ];
1477     SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( SMDSAbs_Volume );
1478     while ( eIt->more() ) {
1479       const SMDS_MeshElement* elem = eIt->next();
1480       if ( elem != myVolume ) {
1481         vNbIt = volNbShared.insert( make_pair( elem, 0 )).first;
1482         (*vNbIt).second++;
1483         if ( vNbIt->second > maxNbShared )
1484           maxNbShared = vNbIt->second;
1485       }
1486     }
1487   }
1488   if ( maxNbShared < 3 )
1489     return isFree; // is free
1490
1491   // find volumes laying on the opposite side of the face
1492   // and sharing all nodes
1493   XYZ intNormal; // internal normal
1494   GetFaceNormal( faceIndex, intNormal.x, intNormal.y, intNormal.z );
1495   if ( IsFaceExternal( faceIndex ))
1496     intNormal = XYZ( -intNormal.x, -intNormal.y, -intNormal.z );
1497   XYZ p0 ( nodes[0] ), baryCenter;
1498   for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end();  ) {
1499     const int& nbShared = (*vNbIt).second;
1500     if ( nbShared >= 3 ) {
1501       SMDS_VolumeTool volume( (*vNbIt).first );
1502       volume.GetBaryCenter( baryCenter.x, baryCenter.y, baryCenter.z );
1503       XYZ intNormal2( baryCenter - p0 );
1504       if ( intNormal.Dot( intNormal2 ) < 0 ) {
1505         // opposite side
1506         if ( nbShared >= nbFaceNodes )
1507         {
1508           // a volume shares the whole facet
1509           if ( otherVol ) *otherVol = vNbIt->first;
1510           return !isFree;
1511         }
1512         ++vNbIt;
1513         continue;
1514       }
1515     }
1516     // remove a volume from volNbShared map
1517     volNbShared.erase( vNbIt++ );
1518   }
1519
1520   // here volNbShared contains only volumes laying on the opposite side of
1521   // the face and sharing 3 or more but not all face nodes with myVolume
1522   if ( volNbShared.size() < 2 ) {
1523     return isFree; // is free
1524   }
1525
1526   // check if the whole area of a face is shared
1527   for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
1528   {
1529     const SMDS_MeshNode* n = nodes[ iNode ];
1530     // check if n is shared by one of volumes of volNbShared
1531     bool isShared = false;
1532     SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( SMDSAbs_Volume );
1533     while ( eIt->more() && !isShared )
1534       isShared = volNbShared.count( eIt->next() );
1535     if ( !isShared )
1536       return isFree;
1537   }
1538   if ( otherVol ) *otherVol = volNbShared.begin()->first;
1539   return !isFree;
1540
1541 //   if ( !myVolume->IsPoly() )
1542 //   {
1543 //     bool isShared[] = { false, false, false, false }; // 4 triangle parts of a quadrangle
1544 //     for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); vNbIt++ ) {
1545 //       SMDS_VolumeTool volume( (*vNbIt).first );
1546 //       bool prevLinkShared = false;
1547 //       int nbSharedLinks = 0;
1548 //       for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) {
1549 //         bool linkShared = volume.IsLinked( nodes[ iNode ], nodes[ iNode + 1] );
1550 //         if ( linkShared )
1551 //           nbSharedLinks++;
1552 //         if ( linkShared && prevLinkShared &&
1553 //              volume.IsLinked( nodes[ iNode - 1 ], nodes[ iNode + 1] ))
1554 //           isShared[ iNode ] = true;
1555 //         prevLinkShared = linkShared;
1556 //       }
1557 //       if ( nbSharedLinks == nbFaceNodes )
1558 //         return !free; // is not free
1559 //       if ( nbFaceNodes == 4 ) {
1560 //         // check traingle parts 1 & 3
1561 //         if ( isShared[1] && isShared[3] )
1562 //           return !free; // is not free
1563 //         // check triangle parts 0 & 2;
1564 //         // 0 part could not be checked in the loop; check it here
1565 //         if ( isShared[2] && prevLinkShared &&
1566 //              volume.IsLinked( nodes[ 0 ], nodes[ 1 ] ) &&
1567 //              volume.IsLinked( nodes[ 1 ], nodes[ 3 ] ) )
1568 //           return !free; // is not free
1569 //       }
1570 //     }
1571 //   }
1572 //  return free;
1573 }
1574
1575 //=======================================================================
1576 //function : GetFaceIndex
1577 //purpose  : Return index of a face formed by theFaceNodes
1578 //=======================================================================
1579
1580 int SMDS_VolumeTool::GetFaceIndex( const set<const SMDS_MeshNode*>& theFaceNodes ) const
1581 {
1582   for ( int iFace = 0; iFace < myNbFaces; iFace++ ) {
1583     const SMDS_MeshNode** nodes = GetFaceNodes( iFace );
1584     int nbFaceNodes = NbFaceNodes( iFace );
1585     set<const SMDS_MeshNode*> nodeSet;
1586     for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
1587       nodeSet.insert( nodes[ iNode ] );
1588     if ( theFaceNodes == nodeSet )
1589       return iFace;
1590   }
1591   return -1;
1592 }
1593
1594 //=======================================================================
1595 //function : GetFaceIndex
1596 //purpose  : Return index of a face formed by theFaceNodes
1597 //=======================================================================
1598
1599 /*int SMDS_VolumeTool::GetFaceIndex( const set<int>& theFaceNodesIndices )
1600 {
1601   for ( int iFace = 0; iFace < myNbFaces; iFace++ ) {
1602     const int* nodes = GetFaceNodesIndices( iFace );
1603     int nbFaceNodes = NbFaceNodes( iFace );
1604     set<int> nodeSet;
1605     for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
1606       nodeSet.insert( nodes[ iNode ] );
1607     if ( theFaceNodesIndices == nodeSet )
1608       return iFace;
1609   }
1610   return -1;
1611 }*/
1612
1613 //=======================================================================
1614 //function : setFace
1615 //purpose  :
1616 //=======================================================================
1617
1618 bool SMDS_VolumeTool::setFace( int faceIndex ) const
1619 {
1620   if ( !myVolume )
1621     return false;
1622
1623   if ( myCurFace == faceIndex )
1624     return true;
1625
1626   myCurFace = -1;
1627
1628   if ( faceIndex < 0 || faceIndex >= NbFaces() )
1629     return false;
1630
1631   if (myFaceNodes != NULL) {
1632     delete [] myFaceNodes;
1633     myFaceNodes = NULL;
1634   }
1635
1636   if (myVolume->IsPoly())
1637   {
1638     if (!myPolyedre) {
1639       MESSAGE("Warning: bad volumic element");
1640       return false;
1641     }
1642
1643     // set face nodes
1644     int iNode;
1645     myFaceNbNodes = myPolyedre->NbFaceNodes(faceIndex + 1);
1646     myFaceNodes = new const SMDS_MeshNode* [myFaceNbNodes + 1];
1647     for ( iNode = 0; iNode < myFaceNbNodes; iNode++ )
1648       myFaceNodes[ iNode ] = myPolyedre->GetFaceNode(faceIndex + 1, iNode + 1);
1649     myFaceNodes[ myFaceNbNodes ] = myFaceNodes[ 0 ]; // last = first
1650
1651     // check orientation
1652     if (myExternalFaces)
1653     {
1654       myCurFace = faceIndex; // avoid infinite recursion in IsFaceExternal()
1655       myExternalFaces = false; // force normal computation by IsFaceExternal()
1656       if ( !IsFaceExternal( faceIndex ))
1657         for ( int i = 0, j = myFaceNbNodes; i < j; ++i, --j )
1658           std::swap( myFaceNodes[i], myFaceNodes[j] );
1659       myExternalFaces = true;
1660     }
1661   }
1662   else
1663   {
1664     if ( !myAllFacesNodeIndices_F )
1665     {
1666       // choose data for an element type
1667       switch ( myVolumeNbNodes ) {
1668       case 4:
1669         myAllFacesNodeIndices_F  = &Tetra_F [0][0];
1670         //myAllFacesNodeIndices_FE = &Tetra_F [0][0];
1671         myAllFacesNodeIndices_RE = &Tetra_RE[0][0];
1672         myAllFacesNbNodes        = Tetra_nbN;
1673         myMaxFaceNbNodes         = sizeof(Tetra_F[0])/sizeof(Tetra_F[0][0]);
1674         break;
1675       case 5:
1676         myAllFacesNodeIndices_F  = &Pyramid_F [0][0];
1677         //myAllFacesNodeIndices_FE = &Pyramid_F [0][0];
1678         myAllFacesNodeIndices_RE = &Pyramid_RE[0][0];
1679         myAllFacesNbNodes        = Pyramid_nbN;
1680         myMaxFaceNbNodes         = sizeof(Pyramid_F[0])/sizeof(Pyramid_F[0][0]);
1681         break;
1682       case 6:
1683         myAllFacesNodeIndices_F  = &Penta_F [0][0];
1684         //myAllFacesNodeIndices_FE = &Penta_FE[0][0];
1685         myAllFacesNodeIndices_RE = &Penta_RE[0][0];
1686         myAllFacesNbNodes        = Penta_nbN;
1687         myMaxFaceNbNodes         = sizeof(Penta_F[0])/sizeof(Penta_F[0][0]);
1688         break;
1689       case 8:
1690         myAllFacesNodeIndices_F  = &Hexa_F [0][0];
1691         ///myAllFacesNodeIndices_FE = &Hexa_FE[0][0];
1692         myAllFacesNodeIndices_RE = &Hexa_RE[0][0];
1693         myAllFacesNbNodes        = Hexa_nbN;
1694         myMaxFaceNbNodes         = sizeof(Hexa_F[0])/sizeof(Hexa_F[0][0]);
1695         break;
1696       case 10:
1697         myAllFacesNodeIndices_F  = &QuadTetra_F [0][0];
1698         //myAllFacesNodeIndices_FE = &QuadTetra_F [0][0];
1699         myAllFacesNodeIndices_RE = &QuadTetra_RE[0][0];
1700         myAllFacesNbNodes        = QuadTetra_nbN;
1701         myMaxFaceNbNodes         = sizeof(QuadTetra_F[0])/sizeof(QuadTetra_F[0][0]);
1702         break;
1703       case 13:
1704         myAllFacesNodeIndices_F  = &QuadPyram_F [0][0];
1705         //myAllFacesNodeIndices_FE = &QuadPyram_F [0][0];
1706         myAllFacesNodeIndices_RE = &QuadPyram_RE[0][0];
1707         myAllFacesNbNodes        = QuadPyram_nbN;
1708         myMaxFaceNbNodes         = sizeof(QuadPyram_F[0])/sizeof(QuadPyram_F[0][0]);
1709         break;
1710       case 15:
1711         myAllFacesNodeIndices_F  = &QuadPenta_F [0][0];
1712         //myAllFacesNodeIndices_FE = &QuadPenta_FE[0][0];
1713         myAllFacesNodeIndices_RE = &QuadPenta_RE[0][0];
1714         myAllFacesNbNodes        = QuadPenta_nbN;
1715         myMaxFaceNbNodes         = sizeof(QuadPenta_F[0])/sizeof(QuadPenta_F[0][0]);
1716         break;
1717       case 20:
1718       case 27:
1719         myAllFacesNodeIndices_F  = &QuadHexa_F [0][0];
1720         //myAllFacesNodeIndices_FE = &QuadHexa_FE[0][0];
1721         myAllFacesNodeIndices_RE = &QuadHexa_RE[0][0];
1722         myAllFacesNbNodes        = QuadHexa_nbN;
1723         myMaxFaceNbNodes         = sizeof(QuadHexa_F[0])/sizeof(QuadHexa_F[0][0]);
1724         if ( !myIgnoreCentralNodes && myVolumeNbNodes == 27 )
1725         {
1726           myAllFacesNodeIndices_F  = &TriQuadHexa_F [0][0];
1727           //myAllFacesNodeIndices_FE = &TriQuadHexa_FE[0][0];
1728           myAllFacesNodeIndices_RE = &TriQuadHexa_RE[0][0];
1729           myAllFacesNbNodes        = TriQuadHexa_nbN;
1730           myMaxFaceNbNodes         = sizeof(TriQuadHexa_F[0])/sizeof(TriQuadHexa_F[0][0]);
1731         }
1732         break;
1733       case 12:
1734         myAllFacesNodeIndices_F  = &HexPrism_F [0][0];
1735         //myAllFacesNodeIndices_FE = &HexPrism_FE[0][0];
1736         myAllFacesNodeIndices_RE = &HexPrism_RE[0][0];
1737         myAllFacesNbNodes        = HexPrism_nbN;
1738         myMaxFaceNbNodes         = sizeof(HexPrism_F[0])/sizeof(HexPrism_F[0][0]);
1739         break;
1740       default:
1741         return false;
1742       }
1743     }
1744     myFaceNbNodes = myAllFacesNbNodes[ faceIndex ];
1745     // if ( myExternalFaces )
1746     //   myFaceNodeIndices = (int*)( myVolForward ? myAllFacesNodeIndices_FE + faceIndex*myMaxFaceNbNodes : myAllFacesNodeIndices_RE + faceIndex*myMaxFaceNbNodes );
1747     // else
1748     //   myFaceNodeIndices = (int*)( myAllFacesNodeIndices_F + faceIndex*myMaxFaceNbNodes );
1749     myFaceNodeIndices = (int*)( myVolForward ? myAllFacesNodeIndices_F + faceIndex*myMaxFaceNbNodes : myAllFacesNodeIndices_RE + faceIndex*myMaxFaceNbNodes );
1750
1751     // set face nodes
1752     myFaceNodes = new const SMDS_MeshNode* [myFaceNbNodes + 1];
1753     for ( int iNode = 0; iNode < myFaceNbNodes; iNode++ )
1754       myFaceNodes[ iNode ] = myVolumeNodes[ myFaceNodeIndices[ iNode ]];
1755     myFaceNodes[ myFaceNbNodes ] = myFaceNodes[ 0 ];
1756   }
1757
1758   myCurFace = faceIndex;
1759
1760   return true;
1761 }
1762
1763 //=======================================================================
1764 //function : GetType
1765 //purpose  : return VolumeType by nb of nodes in a volume
1766 //=======================================================================
1767
1768 SMDS_VolumeTool::VolumeType SMDS_VolumeTool::GetType(int nbNodes)
1769 {
1770   switch ( nbNodes ) {
1771   case 4: return TETRA;
1772   case 5: return PYRAM;
1773   case 6: return PENTA;
1774   case 8: return HEXA;
1775   case 10: return QUAD_TETRA;
1776   case 13: return QUAD_PYRAM;
1777   case 15: return QUAD_PENTA;
1778   case 20:
1779   case 27: return QUAD_HEXA;
1780   case 12: return HEX_PRISM;
1781   default:return UNKNOWN;
1782   }
1783 }
1784
1785 //=======================================================================
1786 //function : NbFaces
1787 //purpose  : return nb of faces by volume type
1788 //=======================================================================
1789
1790 int SMDS_VolumeTool::NbFaces( VolumeType type )
1791 {
1792   switch ( type ) {
1793   case TETRA     :
1794   case QUAD_TETRA: return 4;
1795   case PYRAM     :
1796   case QUAD_PYRAM: return 5;
1797   case PENTA     :
1798   case QUAD_PENTA: return 5;
1799   case HEXA      :
1800   case QUAD_HEXA : return 6;
1801   case HEX_PRISM : return 8;
1802   default:         return 0;
1803   }
1804 }
1805
1806 //================================================================================
1807 /*!
1808  * \brief Useful to know nb of corner nodes of a quadratic volume
1809   * \param type - volume type
1810   * \retval int - nb of corner nodes
1811  */
1812 //================================================================================
1813
1814 int SMDS_VolumeTool::NbCornerNodes(VolumeType type)
1815 {
1816   switch ( type ) {
1817   case TETRA     :
1818   case QUAD_TETRA: return 4;
1819   case PYRAM     :
1820   case QUAD_PYRAM: return 5;
1821   case PENTA     :
1822   case QUAD_PENTA: return 6;
1823   case HEXA      :
1824   case QUAD_HEXA : return 8;
1825   case HEX_PRISM : return 12;
1826   default:         return 0;
1827   }
1828   return 0;
1829 }
1830   //
1831
1832 //=======================================================================
1833 //function : GetFaceNodesIndices
1834 //purpose  : Return the array of face nodes indices
1835 //           To comfort link iteration, the array
1836 //           length == NbFaceNodes( faceIndex ) + 1 and
1837 //           the last node index == the first one.
1838 //=======================================================================
1839
1840 const int* SMDS_VolumeTool::GetFaceNodesIndices(VolumeType type,
1841                                                 int        faceIndex,
1842                                                 bool       external)
1843 {
1844   switch ( type ) {
1845   case TETRA: return Tetra_F[ faceIndex ];
1846   case PYRAM: return Pyramid_F[ faceIndex ];
1847   case PENTA: return external ? Penta_F[ faceIndex ] : Penta_F[ faceIndex ];
1848   case HEXA:  return external ? Hexa_F[ faceIndex ] : Hexa_F[ faceIndex ];
1849   case QUAD_TETRA: return QuadTetra_F[ faceIndex ];
1850   case QUAD_PYRAM: return QuadPyram_F[ faceIndex ];
1851   case QUAD_PENTA: return external ? QuadPenta_F[ faceIndex ] : QuadPenta_F[ faceIndex ];
1852     // what about SMDSEntity_TriQuad_Hexa?
1853   case QUAD_HEXA:  return external ? QuadHexa_F[ faceIndex ] : QuadHexa_F[ faceIndex ];
1854   case HEX_PRISM:  return external ? HexPrism_F[ faceIndex ] : HexPrism_F[ faceIndex ];
1855   default:;
1856   }
1857   return 0;
1858 }
1859
1860 //=======================================================================
1861 //function : NbFaceNodes
1862 //purpose  : Return number of nodes in the array of face nodes
1863 //=======================================================================
1864
1865 int SMDS_VolumeTool::NbFaceNodes(VolumeType type,
1866                                  int        faceIndex )
1867 {
1868   switch ( type ) {
1869   case TETRA: return Tetra_nbN[ faceIndex ];
1870   case PYRAM: return Pyramid_nbN[ faceIndex ];
1871   case PENTA: return Penta_nbN[ faceIndex ];
1872   case HEXA:  return Hexa_nbN[ faceIndex ];
1873   case QUAD_TETRA: return QuadTetra_nbN[ faceIndex ];
1874   case QUAD_PYRAM: return QuadPyram_nbN[ faceIndex ];
1875   case QUAD_PENTA: return QuadPenta_nbN[ faceIndex ];
1876     // what about SMDSEntity_TriQuad_Hexa?
1877   case QUAD_HEXA:  return QuadHexa_nbN[ faceIndex ];
1878   case HEX_PRISM:  return HexPrism_nbN[ faceIndex ];
1879   default:;
1880   }
1881   return 0;
1882 }
1883
1884 //=======================================================================
1885 //function : Element
1886 //purpose  : return element
1887 //=======================================================================
1888
1889 const SMDS_MeshVolume* SMDS_VolumeTool::Element() const
1890 {
1891   return static_cast<const SMDS_MeshVolume*>( myVolume );
1892 }
1893
1894 //=======================================================================
1895 //function : ID
1896 //purpose  : return element ID
1897 //=======================================================================
1898
1899 int SMDS_VolumeTool::ID() const
1900 {
1901   return myVolume ? myVolume->GetID() : 0;
1902 }