src/SMDS/SMDS_VolumeTool.cxx

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