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     if ( !myVolForward )
 527       myCurFace = -1; // previous setFace(0) didn't take myVolForward into account
 528   }
 529   return true;
 530 }
 531
 532 //=======================================================================
 533 //function : Inverse
 534 //purpose  : Inverse volume
 535 //=======================================================================
 536
 537 #define SWAP_NODES(nodes,i1,i2)           \
 538 {                                         \
 539   const SMDS_MeshNode* tmp = nodes[ i1 ]; \
 540   nodes[ i1 ] = nodes[ i2 ];              \
 541   nodes[ i2 ] = tmp;                      \
 542 }
 543 void SMDS_VolumeTool::Inverse ()
 544 {
 545   if ( !myVolume ) return;
 546
 547   if (myVolume->IsPoly()) {
 548     MESSAGE("Warning: attempt to inverse polyhedral volume");
 549     return;
 550   }
 551
 552   myVolForward = !myVolForward;
 553   myCurFace = -1;
 554
 555   // inverse top and bottom faces
 556   switch ( myVolumeNbNodes ) {
 557   case 4:
 558     SWAP_NODES( myVolumeNodes, 1, 2 );
 559     break;
 560   case 5:
 561     SWAP_NODES( myVolumeNodes, 1, 3 );
 562     break;
 563   case 6:
 564     SWAP_NODES( myVolumeNodes, 1, 2 );
 565     SWAP_NODES( myVolumeNodes, 4, 5 );
 566     break;
 567   case 8:
 568     SWAP_NODES( myVolumeNodes, 1, 3 );
 569     SWAP_NODES( myVolumeNodes, 5, 7 );
 570     break;
 571   case 12:
 572     SWAP_NODES( myVolumeNodes, 1, 5 );
 573     SWAP_NODES( myVolumeNodes, 2, 4 );
 574     SWAP_NODES( myVolumeNodes, 7, 11 );
 575     SWAP_NODES( myVolumeNodes, 8, 10 );
 576     break;
 577
 578   case 10:
 579     SWAP_NODES( myVolumeNodes, 1, 2 );
 580     SWAP_NODES( myVolumeNodes, 4, 6 );
 581     SWAP_NODES( myVolumeNodes, 8, 9 );
 582     break;
 583   case 13:
 584     SWAP_NODES( myVolumeNodes, 1, 3 );
 585     SWAP_NODES( myVolumeNodes, 5, 8 );
 586     SWAP_NODES( myVolumeNodes, 6, 7 );
 587     SWAP_NODES( myVolumeNodes, 10, 12 );
 588     break;
 589   case 15:
 590     SWAP_NODES( myVolumeNodes, 1, 2 );
 591     SWAP_NODES( myVolumeNodes, 4, 5 );
 592     SWAP_NODES( myVolumeNodes, 6, 8 );
 593     SWAP_NODES( myVolumeNodes, 9, 11 );
 594     SWAP_NODES( myVolumeNodes, 13, 14 );
 595     break;
 596   case 20:
 597     SWAP_NODES( myVolumeNodes, 1, 3 );
 598     SWAP_NODES( myVolumeNodes, 5, 7 );
 599     SWAP_NODES( myVolumeNodes, 8, 11 );
 600     SWAP_NODES( myVolumeNodes, 9, 10 );
 601     SWAP_NODES( myVolumeNodes, 12, 15 );
 602     SWAP_NODES( myVolumeNodes, 13, 14 );
 603     SWAP_NODES( myVolumeNodes, 17, 19 );
 604     break;
 605   case 27:
 606     SWAP_NODES( myVolumeNodes, 1, 3 );
 607     SWAP_NODES( myVolumeNodes, 5, 7 );
 608     SWAP_NODES( myVolumeNodes, 8, 11 );
 609     SWAP_NODES( myVolumeNodes, 9, 10 );
 610     SWAP_NODES( myVolumeNodes, 12, 15 );
 611     SWAP_NODES( myVolumeNodes, 13, 14 );
 612     SWAP_NODES( myVolumeNodes, 17, 19 );
 613     SWAP_NODES( myVolumeNodes, 21, 24 );
 614     SWAP_NODES( myVolumeNodes, 22, 23 );
 615     break;
 616   default:;
 617   }
 618 }
 619
 620 //=======================================================================
 621 //function : GetVolumeType
 622 //purpose  :
 623 //=======================================================================
 624
 625 SMDS_VolumeTool::VolumeType SMDS_VolumeTool::GetVolumeType() const
 626 {
 627   if ( myPolyedre )
 628     return POLYHEDA;
 629
 630   switch( myVolumeNbNodes ) {
 631   case 4: return TETRA;
 632   case 5: return PYRAM;
 633   case 6: return PENTA;
 634   case 8: return HEXA;
 635   case 12: return HEX_PRISM;
 636   case 10: return QUAD_TETRA;
 637   case 13: return QUAD_PYRAM;
 638   case 15: return QUAD_PENTA;
 639   case 20: return QUAD_HEXA;
 640   case 27: return QUAD_HEXA;
 641   default: break;
 642   }
 643
 644   return UNKNOWN;
 645 }
 646
 647 //=======================================================================
 648 //function : getTetraVolume
 649 //purpose  :
 650 //=======================================================================
 651
 652 static double getTetraVolume(const SMDS_MeshNode* n1,
 653                              const SMDS_MeshNode* n2,
 654                              const SMDS_MeshNode* n3,
 655                              const SMDS_MeshNode* n4)
 656 {
 657   double X1 = n1->X();
 658   double Y1 = n1->Y();
 659   double Z1 = n1->Z();
 660
 661   double X2 = n2->X();
 662   double Y2 = n2->Y();
 663   double Z2 = n2->Z();
 664
 665   double X3 = n3->X();
 666   double Y3 = n3->Y();
 667   double Z3 = n3->Z();
 668
 669   double X4 = n4->X();
 670   double Y4 = n4->Y();
 671   double Z4 = n4->Z();
 672
 673   double Q1 = -(X1-X2)*(Y3*Z4-Y4*Z3);
 674   double Q2 =  (X1-X3)*(Y2*Z4-Y4*Z2);
 675   double R1 = -(X1-X4)*(Y2*Z3-Y3*Z2);
 676   double R2 = -(X2-X3)*(Y1*Z4-Y4*Z1);
 677   double S1 =  (X2-X4)*(Y1*Z3-Y3*Z1);
 678   double S2 = -(X3-X4)*(Y1*Z2-Y2*Z1);
 679
 680   return (Q1+Q2+R1+R2+S1+S2)/6.0;
 681 }
 682
 683 //=======================================================================
 684 //function : GetSize
 685 //purpose  : Return element volume
 686 //=======================================================================
 687
 688 double SMDS_VolumeTool::GetSize() const
 689 {
 690   double V = 0.;
 691   if ( !myVolume )
 692     return 0.;
 693
 694   if ( myVolume->IsPoly() )
 695   {
 696     if ( !myPolyedre )
 697       return 0.;
 698
 699     // split a polyhedron into tetrahedrons
 700
 701     int saveCurFace = myCurFace;
 702     SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* > ( this );
 703     for ( int f = 0; f < NbFaces(); ++f )
 704     {
 705       me->setFace( f );
 706       XYZ area (0,0,0), p1( myFaceNodes[0] );
 707       for ( int n = 0; n < myFaceNbNodes; ++n )
 708       {
 709         XYZ p2( myFaceNodes[ n+1 ]);
 710         area = area + p1.Crossed( p2 );
 711         p1 = p2;
 712       }
 713       V += p1.Dot( area );
 714     }
 715     V /= 6;
 716     if ( saveCurFace > -1 && saveCurFace != myCurFace )
 717       me->setFace( myCurFace );
 718   }
 719   else
 720   {
 721     const static int ind[] = {
 722       0, 1, 3, 6, 11, 23, 31, 44, 58, 78 };
 723     const static int vtab[][4] = { // decomposition into tetra in the order of enum VolumeType
 724       // tetrahedron
 725       { 0, 1, 2, 3 },
 726       // pyramid
 727       { 0, 1, 3, 4 },
 728       { 1, 2, 3, 4 },
 729       // pentahedron
 730       { 0, 1, 2, 3 },
 731       { 1, 5, 3, 4 },
 732       { 1, 5, 2, 3 },
 733       // hexahedron
 734       { 1, 4, 3, 0 },
 735       { 4, 1, 6, 5 },
 736       { 1, 3, 6, 2 },
 737       { 4, 6, 3, 7 },
 738       { 1, 4, 6, 3 },
 739       // hexagonal prism
 740       { 0, 1, 2, 7 },
 741       { 0, 7, 8, 6 },
 742       { 2, 7, 8, 0 },
 743
 744       { 0, 3, 4, 9 },
 745       { 0, 9, 10, 6 },
 746       { 4, 9, 10, 0 },
 747
 748       { 0, 3, 4, 9 },
 749       { 0, 9, 10, 6 },
 750       { 4, 9, 10, 0 },
 751
 752       { 0, 4, 5, 10 },
 753       { 0, 10, 11, 6 },
 754       { 5, 10, 11, 0 },
 755
 756       // quadratic tetrahedron
 757       { 0, 4, 6, 7 },
 758       { 1, 5, 4, 8 },
 759       { 2, 6, 5, 9 },
 760       { 7, 8, 9, 3 },
 761       { 4, 6, 7, 9 },
 762       { 4, 5, 6, 9 },
 763       { 4, 7, 8, 9 },
 764       { 4, 5, 9, 8 },
 765
 766       // quadratic pyramid
 767       { 0, 5, 8, 9 },
 768       { 1, 5,10, 6 },
 769       { 2, 6,11, 7 },
 770       { 3, 7,12, 8 },
 771       { 4, 9,11,10 },
 772       { 4, 9,12,11 },
 773       { 10, 5, 9, 8 },
 774       { 10, 8, 9,12 },
 775       { 10, 8,12, 7 },
 776       { 10, 7,12,11 },
 777       { 10, 7,11, 6 },
 778       { 10, 5, 8, 6 },
 779       { 10, 6, 8, 7 },
 780
 781       // quadratic pentahedron
 782       { 12, 0, 8, 6 },
 783       { 12, 8, 7, 6 },
 784       { 12, 8, 2, 7 },
 785       { 12, 6, 7, 1 },
 786       { 12, 1, 7,13 },
 787       { 12, 7, 2,13 },
 788       { 12, 2,14,13 },
 789
 790       { 12, 3, 9,11 },
 791       { 12,11, 9,10 },
 792       { 12,11,10, 5 },
 793       { 12, 9, 4,10 },
 794       { 12,14, 5,10 },
 795       { 12,14,10, 4 },
 796       { 12,14, 4,13 },
 797
 798       // quadratic hexahedron
 799       { 16, 0,11, 8 },
 800       { 16,11, 9, 8 },
 801       { 16, 8, 9, 1 },
 802       { 16,11, 3,10 },
 803       { 16,11,10, 9 },
 804       { 16,10, 2, 9 },
 805       { 16, 3,19, 2 },
 806       { 16, 2,19,18 },
 807       { 16, 2,18,17 },
 808       { 16, 2,17, 1 },
 809
 810       { 16, 4,12,15 },
 811       { 16,12, 5,13 },
 812       { 16,12,13,15 },
 813       { 16,13, 6,14 },
 814       { 16,13,14,15 },
 815       { 16,14, 7,15 },
 816       { 16, 6, 5,17 },
 817       { 16,18, 6,17 },
 818       { 16,18, 7, 6 },
 819       { 16,18,19, 7 },
 820
 821     };
 822
 823     int type = GetVolumeType();
 824     int n1 = ind[type];
 825     int n2 = ind[type+1];
 826
 827     for (int i = n1; i <  n2; i++) {
 828       V -= getTetraVolume( myVolumeNodes[ vtab[i][0] ],
 829                            myVolumeNodes[ vtab[i][1] ],
 830                            myVolumeNodes[ vtab[i][2] ],
 831                            myVolumeNodes[ vtab[i][3] ]);
 832     }
 833   }
 834   return V;
 835 }
 836
 837 //=======================================================================
 838 //function : GetBaryCenter
 839 //purpose  :
 840 //=======================================================================
 841
 842 bool SMDS_VolumeTool::GetBaryCenter(double & X, double & Y, double & Z) const
 843 {
 844   X = Y = Z = 0.;
 845   if ( !myVolume )
 846     return false;
 847
 848   for ( int i = 0; i < myVolumeNbNodes; i++ ) {
 849     X += myVolumeNodes[ i ]->X();
 850     Y += myVolumeNodes[ i ]->Y();
 851     Z += myVolumeNodes[ i ]->Z();
 852   }
 853   X /= myVolumeNbNodes;
 854   Y /= myVolumeNbNodes;
 855   Z /= myVolumeNbNodes;
 856
 857   return true;
 858 }
 859
 860 //================================================================================
 861 /*!
 862  * \brief Classify a point
 863  *  \param tol - thickness of faces
 864  */
 865 //================================================================================
 866
 867 bool SMDS_VolumeTool::IsOut(double X, double Y, double Z, double tol) const
 868 {
 869   // LIMITATION: for convex volumes only
 870   XYZ p( X,Y,Z );
 871   for ( int iF = 0; iF < myNbFaces; ++iF )
 872   {
 873     XYZ faceNormal;
 874     if ( !GetFaceNormal( iF, faceNormal.x, faceNormal.y, faceNormal.z ))
 875       continue;
 876     if ( !IsFaceExternal( iF ))
 877       faceNormal = XYZ() - faceNormal; // reverse
 878
 879     XYZ face2p( p - XYZ( myFaceNodes[0] ));
 880     if ( face2p.Dot( faceNormal ) > tol )
 881       return true;
 882   }
 883   return false;
 884 }
 885
 886 //=======================================================================
 887 //function : SetExternalNormal
 888 //purpose  : Node order will be so that faces normals are external
 889 //=======================================================================
 890
 891 void SMDS_VolumeTool::SetExternalNormal ()
 892 {
 893   myExternalFaces = true;
 894   myCurFace = -1;
 895 }
 896
 897 //=======================================================================
 898 //function : NbFaceNodes
 899 //purpose  : Return number of nodes in the array of face nodes
 900 //=======================================================================
 901
 902 int SMDS_VolumeTool::NbFaceNodes( int faceIndex ) const
 903 {
 904     if ( !setFace( faceIndex ))
 905       return 0;
 906     return myFaceNbNodes;
 907 }
 908
 909 //=======================================================================
 910 //function : GetFaceNodes
 911 //purpose  : Return pointer to the array of face nodes.
 912 //           To comfort link iteration, the array
 913 //           length == NbFaceNodes( faceIndex ) + 1 and
 914 //           the last node == the first one.
 915 //=======================================================================
 916
 917 const SMDS_MeshNode** SMDS_VolumeTool::GetFaceNodes( int faceIndex ) const
 918 {
 919   if ( !setFace( faceIndex ))
 920     return 0;
 921   return myFaceNodes;
 922 }
 923
 924 //=======================================================================
 925 //function : GetFaceNodesIndices
 926 //purpose  : Return pointer to the array of face nodes indices
 927 //           To comfort link iteration, the array
 928 //           length == NbFaceNodes( faceIndex ) + 1 and
 929 //           the last node index == the first one.
 930 //=======================================================================
 931
 932 const int* SMDS_VolumeTool::GetFaceNodesIndices( int faceIndex ) const
 933 {
 934   if ( !setFace( faceIndex ))
 935     return 0;
 936
 937   if (myPolyedre)
 938   {
 939     SMDS_VolumeTool* me = const_cast< SMDS_VolumeTool* > ( this );
 940     me->myPolyIndices.resize( myFaceNbNodes + 1 );
 941     me->myFaceNodeIndices = & me->myPolyIndices[0];
 942     for ( int i = 0; i <= myFaceNbNodes; ++i )
 943       me->myFaceNodeIndices[i] = myVolume->GetNodeIndex( myFaceNodes[i] );
 944   }
 945   return myFaceNodeIndices;
 946 }
 947
 948 //=======================================================================
 949 //function : GetFaceNodes
 950 //purpose  : Return a set of face nodes.
 951 //=======================================================================
 952
 953 bool SMDS_VolumeTool::GetFaceNodes (int                        faceIndex,
 954                                     set<const SMDS_MeshNode*>& theFaceNodes ) const
 955 {
 956   if ( !setFace( faceIndex ))
 957     return false;
 958
 959   theFaceNodes.clear();
 960   theFaceNodes.insert( myFaceNodes, myFaceNodes + myFaceNbNodes );
 961
 962   return true;
 963 }
 964
 965 //=======================================================================
 966 //function : IsFaceExternal
 967 //purpose  : Check normal orientation of a given face
 968 //=======================================================================
 969
 970 bool SMDS_VolumeTool::IsFaceExternal( int faceIndex ) const
 971 {
 972   if ( myExternalFaces || !myVolume )
 973     return true;
 974
 975   if (myVolume->IsPoly()) {
 976     XYZ aNormal, baryCenter, p0 (myPolyedre->GetFaceNode(faceIndex + 1, 1));
 977     GetFaceNormal(faceIndex, aNormal.x, aNormal.y, aNormal.z);
 978     GetBaryCenter(baryCenter.x, baryCenter.y, baryCenter.z);
 979     XYZ insideVec (baryCenter - p0);
 980     if (insideVec.Dot(aNormal) > 0)
 981       return false;
 982     return true;
 983   }
 984
 985   // switch ( myVolumeNbNodes ) {
 986   // case 4:
 987   // case 5:
 988   // case 10:
 989   // case 13:
 990   //   // only the bottom of a reversed tetrahedron can be internal
 991   //   return ( myVolForward || faceIndex != 0 );
 992   // case 6:
 993   // case 15:
 994   // case 12:
 995   //   // in a forward prism, the top is internal, in a reversed one - bottom
 996   //   return ( myVolForward ? faceIndex != 1 : faceIndex != 0 );
 997   // case 8:
 998   // case 20:
 999   // case 27: {
1000   //   // in a forward hexahedron, even face normal is external, odd - internal
1001   //   bool odd = faceIndex % 2;
1002   //   return ( myVolForward ? !odd : odd );
1003   // }
1004   // default:;
1005   // }
1006   // return false;
1007   return true;
1008 }
1009
1010 //=======================================================================
1011 //function : GetFaceNormal
1012 //purpose  : Return a normal to a face
1013 //=======================================================================
1014
1015 bool SMDS_VolumeTool::GetFaceNormal (int faceIndex, double & X, double & Y, double & Z) const
1016 {
1017   if ( !setFace( faceIndex ))
1018     return false;
1019
1020   const int iQuad = ( myFaceNbNodes > 6 && !myPolyedre ) ? 2 : 1;
1021   XYZ p1 ( myFaceNodes[0*iQuad] );
1022   XYZ p2 ( myFaceNodes[1*iQuad] );
1023   XYZ p3 ( myFaceNodes[2*iQuad] );
1024   XYZ aVec12( p2 - p1 );
1025   XYZ aVec13( p3 - p1 );
1026   XYZ cross = aVec12.Crossed( aVec13 );
1027
1028   if ( myFaceNbNodes >3*iQuad ) {
1029     XYZ p4 ( myFaceNodes[3*iQuad] );
1030     XYZ aVec14( p4 - p1 );
1031     XYZ cross2 = aVec13.Crossed( aVec14 );
1032     cross = cross + cross2;
1033   }
1034
1035   double size = cross.Magnitude();
1036   if ( size <= numeric_limits<double>::min() )
1037     return false;
1038
1039   X = cross.x / size;
1040   Y = cross.y / size;
1041   Z = cross.z / size;
1042
1043   return true;
1044 }
1045
1046 //================================================================================
1047 /*!
1048  * \brief Return barycenter of a face
1049  */
1050 //================================================================================
1051
1052 bool SMDS_VolumeTool::GetFaceBaryCenter (int faceIndex, double & X, double & Y, double & Z) const
1053 {
1054   if ( !setFace( faceIndex ))
1055     return false;
1056
1057   X = Y = Z = 0.0;
1058   for ( int i = 0; i < myFaceNbNodes; ++i )
1059   {
1060     X += myFaceNodes[i]->X() / myFaceNbNodes;
1061     Y += myFaceNodes[i]->Y() / myFaceNbNodes;
1062     Z += myFaceNodes[i]->Z() / myFaceNbNodes;
1063   }
1064   return true;
1065 }
1066
1067 //=======================================================================
1068 //function : GetFaceArea
1069 //purpose  : Return face area
1070 //=======================================================================
1071
1072 double SMDS_VolumeTool::GetFaceArea( int faceIndex ) const
1073 {
1074   if (myVolume->IsPoly()) {
1075     MESSAGE("Warning: attempt to obtain area of a face of polyhedral volume");
1076     return 0;
1077   }
1078
1079   if ( !setFace( faceIndex ))
1080     return 0;
1081
1082   XYZ p1 ( myFaceNodes[0] );
1083   XYZ p2 ( myFaceNodes[1] );
1084   XYZ p3 ( myFaceNodes[2] );
1085   XYZ aVec12( p2 - p1 );
1086   XYZ aVec13( p3 - p1 );
1087   double area = aVec12.Crossed( aVec13 ).Magnitude() * 0.5;
1088
1089   if ( myFaceNbNodes == 4 ) {
1090     XYZ p4 ( myFaceNodes[3] );
1091     XYZ aVec14( p4 - p1 );
1092     area += aVec14.Crossed( aVec13 ).Magnitude() * 0.5;
1093   }
1094   return area;
1095 }
1096
1097 //================================================================================
1098 /*!
1099  * \brief Return index of the node located at face center of a quadratic element like HEX27
1100  */
1101 //================================================================================
1102
1103 int SMDS_VolumeTool::GetCenterNodeIndex( int faceIndex ) const
1104 {
1105   if ( myAllFacesNbNodes && myVolumeNbNodes == 27 ) // classic element with 27 nodes
1106   {
1107     switch ( faceIndex ) {
1108     case 0: return 20;
1109     case 1: return 25;
1110     default:
1111       return faceIndex + 19;
1112     }
1113   }
1114   return -1;
1115 }
1116
1117 //=======================================================================
1118 //function : GetOppFaceIndex
1119 //purpose  : Return index of the opposite face if it exists, else -1.
1120 //=======================================================================
1121
1122 int SMDS_VolumeTool::GetOppFaceIndex( int faceIndex ) const
1123 {
1124   int ind = -1;
1125   if (myPolyedre) {
1126     MESSAGE("Warning: attempt to obtain opposite face on polyhedral volume");
1127     return ind;
1128   }
1129
1130   const int nbHoriFaces = 2;
1131
1132   if ( faceIndex >= 0 && faceIndex < NbFaces() ) {
1133     switch ( myVolumeNbNodes ) {
1134     case 6:
1135     case 15:
1136       if ( faceIndex == 0 || faceIndex == 1 )
1137         ind = 1 - faceIndex;
1138         break;
1139     case 8:
1140     case 12:
1141       if ( faceIndex <= 1 ) // top or bottom
1142         ind = 1 - faceIndex;
1143       else {
1144         const int nbSideFaces = myAllFacesNbNodes[0];
1145         ind = ( faceIndex - nbHoriFaces + nbSideFaces/2 ) % nbSideFaces + nbHoriFaces;
1146       }
1147       break;
1148     case 20:
1149     case 27:
1150       if ( faceIndex <= 1 ) // top or bottom
1151         ind = 1 - faceIndex;
1152       else {
1153         const int nbSideFaces = myAllFacesNbNodes[0] / 2;
1154         ind = ( faceIndex - nbHoriFaces + nbSideFaces/2 ) % nbSideFaces + nbHoriFaces;
1155       }
1156       break;
1157     default:;
1158     }
1159   }
1160   return ind;
1161 }
1162
1163 //=======================================================================
1164 //function : IsLinked
1165 //purpose  : return true if theNode1 is linked with theNode2
1166 // If theIgnoreMediumNodes then corner nodes of quadratic cell are considered linked as well
1167 //=======================================================================
1168
1169 bool SMDS_VolumeTool::IsLinked (const SMDS_MeshNode* theNode1,
1170                                 const SMDS_MeshNode* theNode2,
1171                                 const bool           theIgnoreMediumNodes) const
1172 {
1173   if ( !myVolume )
1174     return false;
1175
1176   if (myVolume->IsPoly()) {
1177     if (!myPolyedre) {
1178       MESSAGE("Warning: bad volumic element");
1179       return false;
1180     }
1181     bool isLinked = false;
1182     int iface;
1183     for (iface = 1; iface <= myNbFaces && !isLinked; iface++) {
1184       int inode, nbFaceNodes = myPolyedre->NbFaceNodes(iface);
1185
1186       for (inode = 1; inode <= nbFaceNodes && !isLinked; inode++) {
1187         const SMDS_MeshNode* curNode = myPolyedre->GetFaceNode(iface, inode);
1188
1189         if (curNode == theNode1 || curNode == theNode2) {
1190           int inextnode = (inode == nbFaceNodes) ? 1 : inode + 1;
1191           const SMDS_MeshNode* nextNode = myPolyedre->GetFaceNode(iface, inextnode);
1192
1193           if ((curNode == theNode1 && nextNode == theNode2) ||
1194               (curNode == theNode2 && nextNode == theNode1)) {
1195             isLinked = true;
1196           }
1197         }
1198       }
1199     }
1200     return isLinked;
1201   }
1202
1203   // find nodes indices
1204   int i1 = -1, i2 = -1, nbFound = 0;
1205   for ( int i = 0; i < myVolumeNbNodes && nbFound < 2; i++ )
1206   {
1207     if ( myVolumeNodes[ i ] == theNode1 )
1208       i1 = i, ++nbFound;
1209     else if ( myVolumeNodes[ i ] == theNode2 )
1210       i2 = i, ++nbFound;
1211   }
1212   return IsLinked( i1, i2 );
1213 }
1214
1215 //=======================================================================
1216 //function : IsLinked
1217 //purpose  : return true if the node with theNode1Index is linked
1218 //           with the node with theNode2Index
1219 // If theIgnoreMediumNodes then corner nodes of quadratic cell are considered linked as well
1220 //=======================================================================
1221
1222 bool SMDS_VolumeTool::IsLinked (const int theNode1Index,
1223                                 const int theNode2Index,
1224                                 bool      theIgnoreMediumNodes) const
1225 {
1226   if ( myVolume->IsPoly() ) {
1227     return IsLinked(myVolumeNodes[theNode1Index], myVolumeNodes[theNode2Index]);
1228   }
1229
1230   int minInd = min( theNode1Index, theNode2Index );
1231   int maxInd = max( theNode1Index, theNode2Index );
1232
1233   if ( minInd < 0 || maxInd > myVolumeNbNodes - 1 || maxInd == minInd )
1234     return false;
1235
1236   VolumeType type = GetVolumeType();
1237   if ( myVolume->IsQuadratic() )
1238   {
1239     int firstMediumInd = myVolume->NbCornerNodes();
1240     if ( minInd >= firstMediumInd )
1241       return false; // both nodes are medium - not linked
1242     if ( maxInd < firstMediumInd ) // both nodes are corners
1243     {
1244       if ( theIgnoreMediumNodes )
1245         type = quadToLinear(type); // to check linkage of corner nodes only
1246       else
1247         return false; // corner nodes are not linked directly in a quadratic cell
1248     }
1249   }
1250
1251   switch ( type ) {
1252   case TETRA:
1253     return true;
1254   case HEXA:
1255     switch ( maxInd - minInd ) {
1256     case 1: return minInd != 3;
1257     case 3: return minInd == 0 || minInd == 4;
1258     case 4: return true;
1259     default:;
1260     }
1261     break;
1262   case PYRAM:
1263     if ( maxInd == 4 )
1264       return true;
1265     switch ( maxInd - minInd ) {
1266     case 1:
1267     case 3: return true;
1268     default:;
1269     }
1270     break;
1271   case PENTA:
1272     switch ( maxInd - minInd ) {
1273     case 1: return minInd != 2;
1274     case 2: return minInd == 0 || minInd == 3;
1275     case 3: return true;
1276     default:;
1277     }
1278     break;
1279   case QUAD_TETRA:
1280     {
1281       switch ( minInd ) {
1282       case 0: if( maxInd==4 ||  maxInd==6 ||  maxInd==7 ) return true;
1283       case 1: if( maxInd==4 ||  maxInd==5 ||  maxInd==8 ) return true;
1284       case 2: if( maxInd==5 ||  maxInd==6 ||  maxInd==9 ) return true;
1285       case 3: if( maxInd==7 ||  maxInd==8 ||  maxInd==9 ) return true;
1286       default:;
1287       }
1288       break;
1289     }
1290   case QUAD_HEXA:
1291     {
1292       switch ( minInd ) {
1293       case 0: if( maxInd==8 ||  maxInd==11 ||  maxInd==16 ) return true;
1294       case 1: if( maxInd==8 ||  maxInd==9 ||  maxInd==17 ) return true;
1295       case 2: if( maxInd==9 ||  maxInd==10 ||  maxInd==18 ) return true;
1296       case 3: if( maxInd==10 ||  maxInd==11 ||  maxInd==19 ) return true;
1297       case 4: if( maxInd==12 ||  maxInd==15 ||  maxInd==16 ) return true;
1298       case 5: if( maxInd==12 ||  maxInd==13 ||  maxInd==17 ) return true;
1299       case 6: if( maxInd==13 ||  maxInd==14 ||  maxInd==18 ) return true;
1300       case 7: if( maxInd==14 ||  maxInd==15 ||  maxInd==19 ) return true;
1301       default:;
1302       }
1303       break;
1304     }
1305   case QUAD_PYRAM:
1306     {
1307       switch ( minInd ) {
1308       case 0: if( maxInd==5 ||  maxInd==8 ||  maxInd==9 ) return true;
1309       case 1: if( maxInd==5 ||  maxInd==6 ||  maxInd==10 ) return true;
1310       case 2: if( maxInd==6 ||  maxInd==7 ||  maxInd==11 ) return true;
1311       case 3: if( maxInd==7 ||  maxInd==8 ||  maxInd==12 ) return true;
1312       case 4: if( maxInd==9 ||  maxInd==10 ||  maxInd==11 ||  maxInd==12 ) return true;
1313       default:;
1314       }
1315       break;
1316     }
1317   case QUAD_PENTA:
1318     {
1319       switch ( minInd ) {
1320       case 0: if( maxInd==6 ||  maxInd==8 ||  maxInd==12 ) return true;
1321       case 1: if( maxInd==6 ||  maxInd==7 ||  maxInd==13 ) return true;
1322       case 2: if( maxInd==7 ||  maxInd==8 ||  maxInd==14 ) return true;
1323       case 3: if( maxInd==9 ||  maxInd==11 ||  maxInd==12 ) return true;
1324       case 4: if( maxInd==9 ||  maxInd==10 ||  maxInd==13 ) return true;
1325       case 5: if( maxInd==10 ||  maxInd==11 ||  maxInd==14 ) return true;
1326       default:;
1327       }
1328       break;
1329     }
1330   case HEX_PRISM:
1331     {
1332       const int diff = maxInd-minInd;
1333       if ( diff > 6  ) return false;// not linked top and bottom
1334       if ( diff == 6 ) return true; // linked top and bottom
1335       return diff == 1 || diff == 7;
1336     }
1337   default:;
1338   }
1339   return false;
1340 }
1341
1342 //=======================================================================
1343 //function : GetNodeIndex
1344 //purpose  : Return an index of theNode
1345 //=======================================================================
1346
1347 int SMDS_VolumeTool::GetNodeIndex(const SMDS_MeshNode* theNode) const
1348 {
1349   if ( myVolume ) {
1350     for ( int i = 0; i < myVolumeNbNodes; i++ ) {
1351       if ( myVolumeNodes[ i ] == theNode )
1352         return i;
1353     }
1354   }
1355   return -1;
1356 }
1357
1358 //================================================================================
1359 /*!
1360  * \brief Fill vector with boundary faces existing in the mesh
1361   * \param faces - vector of found nodes
1362   * \retval int - nb of found faces
1363  */
1364 //================================================================================
1365
1366 int SMDS_VolumeTool::GetAllExistingFaces(vector<const SMDS_MeshElement*> & faces) const
1367 {
1368   faces.clear();
1369   for ( int iF = 0; iF < NbFaces(); ++iF ) {
1370     const SMDS_MeshFace* face = 0;
1371     const SMDS_MeshNode** nodes = GetFaceNodes( iF );
1372     switch ( NbFaceNodes( iF )) {
1373     case 3:
1374       face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2] ); break;
1375     case 4:
1376       face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2], nodes[3] ); break;
1377     case 6:
1378       face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2],
1379                                   nodes[3], nodes[4], nodes[5]); break;
1380     case 8:
1381       face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2], nodes[3],
1382                                   nodes[4], nodes[5], nodes[6], nodes[7]); break;
1383     }
1384     if ( face )
1385       faces.push_back( face );
1386   }
1387   return faces.size();
1388 }
1389
1390
1391 //================================================================================
1392 /*!
1393  * \brief Fill vector with boundary edges existing in the mesh
1394   * \param edges - vector of found edges
1395   * \retval int - nb of found faces
1396  */
1397 //================================================================================
1398
1399 int SMDS_VolumeTool::GetAllExistingEdges(vector<const SMDS_MeshElement*> & edges) const
1400 {
1401   edges.clear();
1402   edges.reserve( myVolumeNbNodes * 2 );
1403   for ( int i = 0; i < myVolumeNbNodes-1; ++i ) {
1404     for ( int j = i + 1; j < myVolumeNbNodes; ++j ) {
1405       if ( IsLinked( i, j )) {
1406         const SMDS_MeshElement* edge =
1407           SMDS_Mesh::FindEdge( myVolumeNodes[i], myVolumeNodes[j] );
1408         if ( edge )
1409           edges.push_back( edge );
1410       }
1411     }
1412   }
1413   return edges.size();
1414 }
1415
1416 //================================================================================
1417 /*!
1418  * \brief Return minimal square distance between connected corner nodes
1419  */
1420 //================================================================================
1421
1422 double SMDS_VolumeTool::MinLinearSize2() const
1423 {
1424   double minSize = 1e+100;
1425   int iQ = myVolume->IsQuadratic() ? 2 : 1;
1426
1427   // store current face data
1428   int curFace = myCurFace, nbN = myFaceNbNodes;
1429   int* ind = myFaceNodeIndices;
1430   myFaceNodeIndices = NULL;
1431   const SMDS_MeshNode** nodes = myFaceNodes;
1432   myFaceNodes = NULL;
1433
1434   // it seems that compute distance twice is faster than organization of a sole computing
1435   myCurFace = -1;
1436   for ( int iF = 0; iF < myNbFaces; ++iF )
1437   {
1438     setFace( iF );
1439     for ( int iN = 0; iN < myFaceNbNodes; iN += iQ )
1440     {
1441       XYZ n1( myFaceNodes[ iN ]);
1442       XYZ n2( myFaceNodes[(iN + iQ) % myFaceNbNodes]);
1443       minSize = std::min( minSize, (n1 - n2).SquareMagnitude());
1444     }
1445   }
1446   // restore current face data
1447   myCurFace = curFace;
1448   myFaceNbNodes = nbN;
1449   myFaceNodeIndices = ind;
1450   delete [] myFaceNodes; myFaceNodes = nodes;
1451
1452   return minSize;
1453 }
1454
1455 //================================================================================
1456 /*!
1457  * \brief Return maximal square distance between connected corner nodes
1458  */
1459 //================================================================================
1460
1461 double SMDS_VolumeTool::MaxLinearSize2() const
1462 {
1463   double maxSize = -1e+100;
1464   int iQ = myVolume->IsQuadratic() ? 2 : 1;
1465
1466   // store current face data
1467   int curFace = myCurFace, nbN = myFaceNbNodes;
1468   int* ind = myFaceNodeIndices;
1469   myFaceNodeIndices = NULL;
1470   const SMDS_MeshNode** nodes = myFaceNodes;
1471   myFaceNodes = NULL;
1472
1473   // it seems that compute distance twice is faster than organization of a sole computing
1474   myCurFace = -1;
1475   for ( int iF = 0; iF < myNbFaces; ++iF )
1476   {
1477     setFace( iF );
1478     for ( int iN = 0; iN < myFaceNbNodes; iN += iQ )
1479     {
1480       XYZ n1( myFaceNodes[ iN ]);
1481       XYZ n2( myFaceNodes[(iN + iQ) % myFaceNbNodes]);
1482       maxSize = std::max( maxSize, (n1 - n2).SquareMagnitude());
1483     }
1484   }
1485   // restore current face data
1486   myCurFace         = curFace;
1487   myFaceNbNodes     = nbN;
1488   myFaceNodeIndices = ind;
1489   delete [] myFaceNodes; myFaceNodes = nodes;
1490
1491   return maxSize;
1492 }
1493
1494 //================================================================================
1495 /*!
1496  * \brief check that only one volume is build on the face nodes
1497  *
1498  * If a face is shared by one of <ignoreVolumes>, it is considered free
1499  */
1500 //================================================================================
1501
1502 bool SMDS_VolumeTool::IsFreeFace( int faceIndex, const SMDS_MeshElement** otherVol/*=0*/ ) const
1503 {
1504   const bool isFree = true;
1505
1506   if (!setFace( faceIndex ))
1507     return !isFree;
1508
1509   const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
1510   const int nbFaceNodes = myFaceNbNodes;
1511
1512   // evaluate nb of face nodes shared by other volumes
1513   int maxNbShared = -1;
1514   typedef map< const SMDS_MeshElement*, int > TElemIntMap;
1515   TElemIntMap volNbShared;
1516   TElemIntMap::iterator vNbIt;
1517   for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) {
1518     const SMDS_MeshNode* n = nodes[ iNode ];
1519     SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( SMDSAbs_Volume );
1520     while ( eIt->more() ) {
1521       const SMDS_MeshElement* elem = eIt->next();
1522       if ( elem != myVolume ) {
1523         vNbIt = volNbShared.insert( make_pair( elem, 0 )).first;
1524         (*vNbIt).second++;
1525         if ( vNbIt->second > maxNbShared )
1526           maxNbShared = vNbIt->second;
1527       }
1528     }
1529   }
1530   if ( maxNbShared < 3 )
1531     return isFree; // is free
1532
1533   // find volumes laying on the opposite side of the face
1534   // and sharing all nodes
1535   XYZ intNormal; // internal normal
1536   GetFaceNormal( faceIndex, intNormal.x, intNormal.y, intNormal.z );
1537   if ( IsFaceExternal( faceIndex ))
1538     intNormal = XYZ( -intNormal.x, -intNormal.y, -intNormal.z );
1539   XYZ p0 ( nodes[0] ), baryCenter;
1540   for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end();  ) {
1541     const int& nbShared = (*vNbIt).second;
1542     if ( nbShared >= 3 ) {
1543       SMDS_VolumeTool volume( (*vNbIt).first );
1544       volume.GetBaryCenter( baryCenter.x, baryCenter.y, baryCenter.z );
1545       XYZ intNormal2( baryCenter - p0 );
1546       if ( intNormal.Dot( intNormal2 ) < 0 ) {
1547         // opposite side
1548         if ( nbShared >= nbFaceNodes )
1549         {
1550           // a volume shares the whole facet
1551           if ( otherVol ) *otherVol = vNbIt->first;
1552           return !isFree;
1553         }
1554         ++vNbIt;
1555         continue;
1556       }
1557     }
1558     // remove a volume from volNbShared map
1559     volNbShared.erase( vNbIt++ );
1560   }
1561
1562   // here volNbShared contains only volumes laying on the opposite side of
1563   // the face and sharing 3 or more but not all face nodes with myVolume
1564   if ( volNbShared.size() < 2 ) {
1565     return isFree; // is free
1566   }
1567
1568   // check if the whole area of a face is shared
1569   for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
1570   {
1571     const SMDS_MeshNode* n = nodes[ iNode ];
1572     // check if n is shared by one of volumes of volNbShared
1573     bool isShared = false;
1574     SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( SMDSAbs_Volume );
1575     while ( eIt->more() && !isShared )
1576       isShared = volNbShared.count( eIt->next() );
1577     if ( !isShared )
1578       return isFree;
1579   }
1580   if ( otherVol ) *otherVol = volNbShared.begin()->first;
1581   return !isFree;
1582
1583 //   if ( !myVolume->IsPoly() )
1584 //   {
1585 //     bool isShared[] = { false, false, false, false }; // 4 triangle parts of a quadrangle
1586 //     for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); vNbIt++ ) {
1587 //       SMDS_VolumeTool volume( (*vNbIt).first );
1588 //       bool prevLinkShared = false;
1589 //       int nbSharedLinks = 0;
1590 //       for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) {
1591 //         bool linkShared = volume.IsLinked( nodes[ iNode ], nodes[ iNode + 1] );
1592 //         if ( linkShared )
1593 //           nbSharedLinks++;
1594 //         if ( linkShared && prevLinkShared &&
1595 //              volume.IsLinked( nodes[ iNode - 1 ], nodes[ iNode + 1] ))
1596 //           isShared[ iNode ] = true;
1597 //         prevLinkShared = linkShared;
1598 //       }
1599 //       if ( nbSharedLinks == nbFaceNodes )
1600 //         return !free; // is not free
1601 //       if ( nbFaceNodes == 4 ) {
1602 //         // check traingle parts 1 & 3
1603 //         if ( isShared[1] && isShared[3] )
1604 //           return !free; // is not free
1605 //         // check triangle parts 0 & 2;
1606 //         // 0 part could not be checked in the loop; check it here
1607 //         if ( isShared[2] && prevLinkShared &&
1608 //              volume.IsLinked( nodes[ 0 ], nodes[ 1 ] ) &&
1609 //              volume.IsLinked( nodes[ 1 ], nodes[ 3 ] ) )
1610 //           return !free; // is not free
1611 //       }
1612 //     }
1613 //   }
1614 //  return free;
1615 }
1616
1617 //=======================================================================
1618 //function : GetFaceIndex
1619 //purpose  : Return index of a face formed by theFaceNodes
1620 //=======================================================================
1621
1622 int SMDS_VolumeTool::GetFaceIndex( const set<const SMDS_MeshNode*>& theFaceNodes ) const
1623 {
1624   for ( int iFace = 0; iFace < myNbFaces; iFace++ ) {
1625     const SMDS_MeshNode** nodes = GetFaceNodes( iFace );
1626     int nbFaceNodes = NbFaceNodes( iFace );
1627     set<const SMDS_MeshNode*> nodeSet;
1628     for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
1629       nodeSet.insert( nodes[ iNode ] );
1630     if ( theFaceNodes == nodeSet )
1631       return iFace;
1632   }
1633   return -1;
1634 }
1635
1636 //=======================================================================
1637 //function : GetFaceIndex
1638 //purpose  : Return index of a face formed by theFaceNodes
1639 //=======================================================================
1640
1641 /*int SMDS_VolumeTool::GetFaceIndex( const set<int>& theFaceNodesIndices )
1642 {
1643   for ( int iFace = 0; iFace < myNbFaces; iFace++ ) {
1644     const int* nodes = GetFaceNodesIndices( iFace );
1645     int nbFaceNodes = NbFaceNodes( iFace );
1646     set<int> nodeSet;
1647     for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
1648       nodeSet.insert( nodes[ iNode ] );
1649     if ( theFaceNodesIndices == nodeSet )
1650       return iFace;
1651   }
1652   return -1;
1653 }*/
1654
1655 //=======================================================================
1656 //function : setFace
1657 //purpose  :
1658 //=======================================================================
1659
1660 bool SMDS_VolumeTool::setFace( int faceIndex ) const
1661 {
1662   if ( !myVolume )
1663     return false;
1664
1665   if ( myCurFace == faceIndex )
1666     return true;
1667
1668   myCurFace = -1;
1669
1670   if ( faceIndex < 0 || faceIndex >= NbFaces() )
1671     return false;
1672
1673   if (myFaceNodes != NULL) {
1674     delete [] myFaceNodes;
1675     myFaceNodes = NULL;
1676   }
1677
1678   if (myVolume->IsPoly())
1679   {
1680     if (!myPolyedre) {
1681       MESSAGE("Warning: bad volumic element");
1682       return false;
1683     }
1684
1685     // set face nodes
1686     int iNode;
1687     myFaceNbNodes = myPolyedre->NbFaceNodes(faceIndex + 1);
1688     myFaceNodes = new const SMDS_MeshNode* [myFaceNbNodes + 1];
1689     for ( iNode = 0; iNode < myFaceNbNodes; iNode++ )
1690       myFaceNodes[ iNode ] = myPolyedre->GetFaceNode(faceIndex + 1, iNode + 1);
1691     myFaceNodes[ myFaceNbNodes ] = myFaceNodes[ 0 ]; // last = first
1692
1693     // check orientation
1694     if (myExternalFaces)
1695     {
1696       myCurFace = faceIndex; // avoid infinite recursion in IsFaceExternal()
1697       myExternalFaces = false; // force normal computation by IsFaceExternal()
1698       if ( !IsFaceExternal( faceIndex ))
1699         for ( int i = 0, j = myFaceNbNodes; i < j; ++i, --j )
1700           std::swap( myFaceNodes[i], myFaceNodes[j] );
1701       myExternalFaces = true;
1702     }
1703   }
1704   else
1705   {
1706     if ( !myAllFacesNodeIndices_F )
1707     {
1708       // choose data for an element type
1709       switch ( myVolumeNbNodes ) {
1710       case 4:
1711         myAllFacesNodeIndices_F  = &Tetra_F [0][0];
1712         //myAllFacesNodeIndices_FE = &Tetra_F [0][0];
1713         myAllFacesNodeIndices_RE = &Tetra_RE[0][0];
1714         myAllFacesNbNodes        = Tetra_nbN;
1715         myMaxFaceNbNodes         = sizeof(Tetra_F[0])/sizeof(Tetra_F[0][0]);
1716         break;
1717       case 5:
1718         myAllFacesNodeIndices_F  = &Pyramid_F [0][0];
1719         //myAllFacesNodeIndices_FE = &Pyramid_F [0][0];
1720         myAllFacesNodeIndices_RE = &Pyramid_RE[0][0];
1721         myAllFacesNbNodes        = Pyramid_nbN;
1722         myMaxFaceNbNodes         = sizeof(Pyramid_F[0])/sizeof(Pyramid_F[0][0]);
1723         break;
1724       case 6:
1725         myAllFacesNodeIndices_F  = &Penta_F [0][0];
1726         //myAllFacesNodeIndices_FE = &Penta_FE[0][0];
1727         myAllFacesNodeIndices_RE = &Penta_RE[0][0];
1728         myAllFacesNbNodes        = Penta_nbN;
1729         myMaxFaceNbNodes         = sizeof(Penta_F[0])/sizeof(Penta_F[0][0]);
1730         break;
1731       case 8:
1732         myAllFacesNodeIndices_F  = &Hexa_F [0][0];
1733         ///myAllFacesNodeIndices_FE = &Hexa_FE[0][0];
1734         myAllFacesNodeIndices_RE = &Hexa_RE[0][0];
1735         myAllFacesNbNodes        = Hexa_nbN;
1736         myMaxFaceNbNodes         = sizeof(Hexa_F[0])/sizeof(Hexa_F[0][0]);
1737         break;
1738       case 10:
1739         myAllFacesNodeIndices_F  = &QuadTetra_F [0][0];
1740         //myAllFacesNodeIndices_FE = &QuadTetra_F [0][0];
1741         myAllFacesNodeIndices_RE = &QuadTetra_RE[0][0];
1742         myAllFacesNbNodes        = QuadTetra_nbN;
1743         myMaxFaceNbNodes         = sizeof(QuadTetra_F[0])/sizeof(QuadTetra_F[0][0]);
1744         break;
1745       case 13:
1746         myAllFacesNodeIndices_F  = &QuadPyram_F [0][0];
1747         //myAllFacesNodeIndices_FE = &QuadPyram_F [0][0];
1748         myAllFacesNodeIndices_RE = &QuadPyram_RE[0][0];
1749         myAllFacesNbNodes        = QuadPyram_nbN;
1750         myMaxFaceNbNodes         = sizeof(QuadPyram_F[0])/sizeof(QuadPyram_F[0][0]);
1751         break;
1752       case 15:
1753         myAllFacesNodeIndices_F  = &QuadPenta_F [0][0];
1754         //myAllFacesNodeIndices_FE = &QuadPenta_FE[0][0];
1755         myAllFacesNodeIndices_RE = &QuadPenta_RE[0][0];
1756         myAllFacesNbNodes        = QuadPenta_nbN;
1757         myMaxFaceNbNodes         = sizeof(QuadPenta_F[0])/sizeof(QuadPenta_F[0][0]);
1758         break;
1759       case 20:
1760       case 27:
1761         myAllFacesNodeIndices_F  = &QuadHexa_F [0][0];
1762         //myAllFacesNodeIndices_FE = &QuadHexa_FE[0][0];
1763         myAllFacesNodeIndices_RE = &QuadHexa_RE[0][0];
1764         myAllFacesNbNodes        = QuadHexa_nbN;
1765         myMaxFaceNbNodes         = sizeof(QuadHexa_F[0])/sizeof(QuadHexa_F[0][0]);
1766         if ( !myIgnoreCentralNodes && myVolumeNbNodes == 27 )
1767         {
1768           myAllFacesNodeIndices_F  = &TriQuadHexa_F [0][0];
1769           //myAllFacesNodeIndices_FE = &TriQuadHexa_FE[0][0];
1770           myAllFacesNodeIndices_RE = &TriQuadHexa_RE[0][0];
1771           myAllFacesNbNodes        = TriQuadHexa_nbN;
1772           myMaxFaceNbNodes         = sizeof(TriQuadHexa_F[0])/sizeof(TriQuadHexa_F[0][0]);
1773         }
1774         break;
1775       case 12:
1776         myAllFacesNodeIndices_F  = &HexPrism_F [0][0];
1777         //myAllFacesNodeIndices_FE = &HexPrism_FE[0][0];
1778         myAllFacesNodeIndices_RE = &HexPrism_RE[0][0];
1779         myAllFacesNbNodes        = HexPrism_nbN;
1780         myMaxFaceNbNodes         = sizeof(HexPrism_F[0])/sizeof(HexPrism_F[0][0]);
1781         break;
1782       default:
1783         return false;
1784       }
1785     }
1786     myFaceNbNodes = myAllFacesNbNodes[ faceIndex ];
1787     // if ( myExternalFaces )
1788     //   myFaceNodeIndices = (int*)( myVolForward ? myAllFacesNodeIndices_FE + faceIndex*myMaxFaceNbNodes : myAllFacesNodeIndices_RE + faceIndex*myMaxFaceNbNodes );
1789     // else
1790     //   myFaceNodeIndices = (int*)( myAllFacesNodeIndices_F + faceIndex*myMaxFaceNbNodes );
1791     myFaceNodeIndices = (int*)( myVolForward ? myAllFacesNodeIndices_F + faceIndex*myMaxFaceNbNodes : myAllFacesNodeIndices_RE + faceIndex*myMaxFaceNbNodes );
1792
1793     // set face nodes
1794     myFaceNodes = new const SMDS_MeshNode* [myFaceNbNodes + 1];
1795     for ( int iNode = 0; iNode < myFaceNbNodes; iNode++ )
1796       myFaceNodes[ iNode ] = myVolumeNodes[ myFaceNodeIndices[ iNode ]];
1797     myFaceNodes[ myFaceNbNodes ] = myFaceNodes[ 0 ];
1798   }
1799
1800   myCurFace = faceIndex;
1801
1802   return true;
1803 }
1804
1805 //=======================================================================
1806 //function : GetType
1807 //purpose  : return VolumeType by nb of nodes in a volume
1808 //=======================================================================
1809
1810 SMDS_VolumeTool::VolumeType SMDS_VolumeTool::GetType(int nbNodes)
1811 {
1812   switch ( nbNodes ) {
1813   case 4: return TETRA;
1814   case 5: return PYRAM;
1815   case 6: return PENTA;
1816   case 8: return HEXA;
1817   case 10: return QUAD_TETRA;
1818   case 13: return QUAD_PYRAM;
1819   case 15: return QUAD_PENTA;
1820   case 20:
1821   case 27: return QUAD_HEXA;
1822   case 12: return HEX_PRISM;
1823   default:return UNKNOWN;
1824   }
1825 }
1826
1827 //=======================================================================
1828 //function : NbFaces
1829 //purpose  : return nb of faces by volume type
1830 //=======================================================================
1831
1832 int SMDS_VolumeTool::NbFaces( VolumeType type )
1833 {
1834   switch ( type ) {
1835   case TETRA     :
1836   case QUAD_TETRA: return 4;
1837   case PYRAM     :
1838   case QUAD_PYRAM: return 5;
1839   case PENTA     :
1840   case QUAD_PENTA: return 5;
1841   case HEXA      :
1842   case QUAD_HEXA : return 6;
1843   case HEX_PRISM : return 8;
1844   default:         return 0;
1845   }
1846 }
1847
1848 //================================================================================
1849 /*!
1850  * \brief Useful to know nb of corner nodes of a quadratic volume
1851   * \param type - volume type
1852   * \retval int - nb of corner nodes
1853  */
1854 //================================================================================
1855
1856 int SMDS_VolumeTool::NbCornerNodes(VolumeType type)
1857 {
1858   switch ( type ) {
1859   case TETRA     :
1860   case QUAD_TETRA: return 4;
1861   case PYRAM     :
1862   case QUAD_PYRAM: return 5;
1863   case PENTA     :
1864   case QUAD_PENTA: return 6;
1865   case HEXA      :
1866   case QUAD_HEXA : return 8;
1867   case HEX_PRISM : return 12;
1868   default:         return 0;
1869   }
1870   return 0;
1871 }
1872   //
1873
1874 //=======================================================================
1875 //function : GetFaceNodesIndices
1876 //purpose  : Return the array of face nodes indices
1877 //           To comfort link iteration, the array
1878 //           length == NbFaceNodes( faceIndex ) + 1 and
1879 //           the last node index == the first one.
1880 //=======================================================================
1881
1882 const int* SMDS_VolumeTool::GetFaceNodesIndices(VolumeType type,
1883                                                 int        faceIndex,
1884                                                 bool       external)
1885 {
1886   switch ( type ) {
1887   case TETRA: return Tetra_F[ faceIndex ];
1888   case PYRAM: return Pyramid_F[ faceIndex ];
1889   case PENTA: return external ? Penta_F[ faceIndex ] : Penta_F[ faceIndex ];
1890   case HEXA:  return external ? Hexa_F[ faceIndex ] : Hexa_F[ faceIndex ];
1891   case QUAD_TETRA: return QuadTetra_F[ faceIndex ];
1892   case QUAD_PYRAM: return QuadPyram_F[ faceIndex ];
1893   case QUAD_PENTA: return external ? QuadPenta_F[ faceIndex ] : QuadPenta_F[ faceIndex ];
1894     // what about SMDSEntity_TriQuad_Hexa?
1895   case QUAD_HEXA:  return external ? QuadHexa_F[ faceIndex ] : QuadHexa_F[ faceIndex ];
1896   case HEX_PRISM:  return external ? HexPrism_F[ faceIndex ] : HexPrism_F[ faceIndex ];
1897   default:;
1898   }
1899   return 0;
1900 }
1901
1902 //=======================================================================
1903 //function : NbFaceNodes
1904 //purpose  : Return number of nodes in the array of face nodes
1905 //=======================================================================
1906
1907 int SMDS_VolumeTool::NbFaceNodes(VolumeType type,
1908                                  int        faceIndex )
1909 {
1910   switch ( type ) {
1911   case TETRA: return Tetra_nbN[ faceIndex ];
1912   case PYRAM: return Pyramid_nbN[ faceIndex ];
1913   case PENTA: return Penta_nbN[ faceIndex ];
1914   case HEXA:  return Hexa_nbN[ faceIndex ];
1915   case QUAD_TETRA: return QuadTetra_nbN[ faceIndex ];
1916   case QUAD_PYRAM: return QuadPyram_nbN[ faceIndex ];
1917   case QUAD_PENTA: return QuadPenta_nbN[ faceIndex ];
1918     // what about SMDSEntity_TriQuad_Hexa?
1919   case QUAD_HEXA:  return QuadHexa_nbN[ faceIndex ];
1920   case HEX_PRISM:  return HexPrism_nbN[ faceIndex ];
1921   default:;
1922   }
1923   return 0;
1924 }
1925
1926 //=======================================================================
1927 //function : Element
1928 //purpose  : return element
1929 //=======================================================================
1930
1931 const SMDS_MeshVolume* SMDS_VolumeTool::Element() const
1932 {
1933   return static_cast<const SMDS_MeshVolume*>( myVolume );
1934 }
1935
1936 //=======================================================================
1937 //function : ID
1938 //purpose  : return element ID
1939 //=======================================================================
1940
1941 int SMDS_VolumeTool::ID() const
1942 {
1943   return myVolume ? myVolume->GetID() : 0;
1944 }