src/SMDS/SMDS_VolumeTool.cxx

   1 // Copyright (C) 2007-2015  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 namespace
  47 {
  48 // ======================================================
  49 // Node indices in faces depending on volume orientation
  50 // making most faces normals external
  51 // ======================================================
  52 // For all elements, 0-th face is bottom based on the first nodes.
  53 // For prismatic elements (tetra,hexa,prisms), 1-th face is a top one.
  54 // For all elements, side faces follow order of bottom nodes
  55 // ======================================================
  56
  57 /*
  58 //           N3
  59 //           +
  60 //          /|\
  61 //         / | \
  62 //        /  |  \
  63 //    N0 +---|---+ N1                TETRAHEDRON
  64 //       \   |   /
  65 //        \  |  /
  66 //         \ | /
  67 //          \|/
  68 //           +
  69 //           N2
  70 */
  71 static int Tetra_F [4][4] = { // FORWARD == EXTERNAL
  72   { 0, 1, 2, 0 },              // All faces have external normals
  73   { 0, 3, 1, 0 },
  74   { 1, 3, 2, 1 },
  75   { 0, 2, 3, 0 }};
  76 static int Tetra_RE [4][4] = { // REVERSED -> FORWARD (EXTERNAL)
  77   { 0, 2, 1, 0 },              // All faces have external normals
  78   { 0, 1, 3, 0 },
  79   { 1, 2, 3, 1 },
  80   { 0, 3, 2, 0 }};
  81 static int Tetra_nbN [] = { 3, 3, 3, 3 };
  82
  83 //
  84 //     PYRAMID
  85 //
  86 static int Pyramid_F [5][5] = { // FORWARD == EXTERNAL
  87   { 0, 1, 2, 3, 0 },            // All faces have external normals
  88   { 0, 4, 1, 0, 4 },
  89   { 1, 4, 2, 1, 4 },
  90   { 2, 4, 3, 2, 4 },
  91   { 3, 4, 0, 3, 4 }
  92 };
  93 static int Pyramid_RE [5][5] = { // REVERSED -> FORWARD (EXTERNAL)
  94   { 0, 3, 2, 1, 0 },             // All faces but a bottom have external normals
  95   { 0, 1, 4, 0, 4 },
  96   { 1, 2, 4, 1, 4 },
  97   { 2, 3, 4, 2, 4 },
  98   { 3, 0, 4, 3, 4 }};
  99 static int Pyramid_nbN [] = { 4, 3, 3, 3, 3 };
 100
 101 /*
 102 //            + N4
 103 //           /|\
 104 //          / | \
 105 //         /  |  \
 106 //        /   |   \
 107 //    N3 +---------+ N5
 108 //       |    |    |
 109 //       |    + N1 |
 110 //       |   / \   |                PENTAHEDRON
 111 //       |  /   \  |
 112 //       | /     \ |
 113 //       |/       \|
 114 //    N0 +---------+ N2
 115 */
 116 static int Penta_F [5][5] = { // FORWARD
 117   { 0, 1, 2, 0, 0 },          // All faces have external normals
 118   { 3, 5, 4, 3, 3 },          // 0 is bottom, 1 is top face
 119   { 0, 3, 4, 1, 0 },
 120   { 1, 4, 5, 2, 1 },
 121   { 0, 2, 5, 3, 0 }};
 122 static int Penta_RE [5][5] = { // REVERSED -> EXTERNAL
 123   { 0, 2, 1, 0, 0 },
 124   { 3, 4, 5, 3, 3 },
 125   { 0, 1, 4, 3, 0 },
 126   { 1, 2, 5, 4, 1 },
 127   { 0, 3, 5, 2, 0 }};
 128 static int Penta_nbN [] = { 3, 3, 4, 4, 4 };
 129
 130 /*
 131 //         N5+----------+N6
 132 //          /|         /|
 133 //         / |        / |
 134 //        /  |       /  |
 135 //     N4+----------+N7 |
 136 //       |   |      |   |           HEXAHEDRON
 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 static int Hexa_oppF[] = { 1, 0, 4, 5, 2, 3 }; // oppopsite facet indices
 159
 160 /*
 161 //      N8 +------+ N9
 162 //        /        \
 163 //       /          \
 164 //   N7 +            + N10
 165 //       \          /
 166 //        \        /
 167 //      N6 +------+ N11
 168 //                             HEXAGONAL PRISM
 169 //      N2 +------+ N3
 170 //        /        \
 171 //       /          \
 172 //   N1 +            + N4
 173 //       \          /
 174 //        \        /
 175 //      N0 +------+ N5
 176 */
 177 static int HexPrism_F [8][7] = { // FORWARD
 178   { 0, 1, 2, 3, 4, 5, 0 },
 179   { 6,11,10, 9, 8, 7, 6 },
 180   { 0, 6, 7, 1, 0, 0, 0 },
 181   { 1, 7, 8, 2, 1, 1, 1 },
 182   { 2, 8, 9, 3, 2, 2, 2 },
 183   { 3, 9,10, 4, 3, 3, 3 },
 184   { 4,10,11, 5, 4, 4, 4 },
 185   { 5,11, 6, 0, 5, 5, 5 }};
 186 static int HexPrism_RE [8][7] = { // REVERSED -> EXTERNAL
 187   { 0, 5, 4, 3, 2, 1, 0 },
 188   { 6,11,10, 9, 8, 7, 6 },
 189   { 0, 6, 7, 1, 0, 0, 0 },
 190   { 1, 7, 8, 2, 1, 1, 1 },
 191   { 2, 8, 9, 3, 2, 2, 2 },
 192   { 3, 9,10, 4, 3, 3, 3 },
 193   { 4,10,11, 5, 4, 4, 4 },
 194   { 5,11, 6, 0, 5, 5, 5 }};
 195 static int HexPrism_nbN [] = { 6, 6, 4, 4, 4, 4, 4, 4 };
 196
 197
 198 /*
 199 //           N3
 200 //           +
 201 //          /|\
 202 //        7/ | \8
 203 //        /  |4 \                    QUADRATIC
 204 //    N0 +---|---+ N1                TETRAHEDRON
 205 //       \   +9  /
 206 //        \  |  /
 207 //        6\ | /5
 208 //          \|/
 209 //           +
 210 //           N2
 211 */
 212 static int QuadTetra_F [4][7] = { // FORWARD
 213   { 0, 4, 1, 5, 2, 6, 0 },        // All faces have external normals
 214   { 0, 7, 3, 8, 1, 4, 0 },
 215   { 1, 8, 3, 9, 2, 5, 1 },
 216   { 0, 6, 2, 9, 3, 7, 0 }};
 217 static int QuadTetra_RE [4][7] = { // REVERSED -> FORWARD (EXTERNAL)
 218   { 0, 6, 2, 5, 1, 4, 0 },         // All faces have external normals
 219   { 0, 4, 1, 8, 3, 7, 0 },
 220   { 1, 5, 2, 9, 3, 8, 1 },
 221   { 0, 7, 3, 9, 2, 6, 0 }};
 222 static int QuadTetra_nbN [] = { 6, 6, 6, 6 };
 223
 224 //
 225 //     QUADRATIC
 226 //     PYRAMID
 227 //
 228 //            +4
 229 //
 230 //
 231 //       10+-----+11
 232 //         |     |        9 - middle point for (0,4) etc.
 233 //         |     |
 234 //        9+-----+12
 235 //
 236 //            6
 237 //      1+----+----+2
 238 //       |         |
 239 //       |         |
 240 //      5+         +7
 241 //       |         |
 242 //       |         |
 243 //      0+----+----+3
 244 //            8
 245 static int QuadPyram_F [5][9] = {  // FORWARD
 246   { 0, 5, 1, 6, 2, 7, 3, 8, 0 },   // All faces have external normals
 247   { 0, 9, 4, 10,1, 5, 0, 4, 4 },
 248   { 1, 10,4, 11,2, 6, 1, 4, 4 },
 249   { 2, 11,4, 12,3, 7, 2, 4, 4 },
 250   { 3, 12,4, 9, 0, 8, 3, 4, 4 }};
 251 static int QuadPyram_RE [5][9] = { // REVERSED -> FORWARD (EXTERNAL)
 252   { 0, 8, 3, 7, 2, 6, 1, 5, 0 },   // All faces but a bottom have external normals
 253   { 0, 5, 1, 10,4, 9, 0, 4, 4 },
 254   { 1, 6, 2, 11,4, 10,1, 4, 4 },
 255   { 2, 7, 3, 12,4, 11,2, 4, 4 },
 256   { 3, 8, 0, 9, 4, 12,3, 4, 4 }};
 257 static int QuadPyram_nbN [] = { 8, 6, 6, 6, 6 };
 258
 259 /*
 260 //            + N4
 261 //           /|\
 262 //         9/ | \10
 263 //         /  |  \
 264 //        /   |   \
 265 //    N3 +----+----+ N5
 266 //       |    |11  |
 267 //       |    |    |
 268 //       |    +13  |                QUADRATIC
 269 //       |    |    |                PENTAHEDRON
 270 //     12+    |    +14
 271 //       |    |    |
 272 //       |    |    |
 273 //       |    + N1 |
 274 //       |   / \   |
 275 //       | 6/   \7 |
 276 //       | /     \ |
 277 //       |/       \|
 278 //    N0 +---------+ N2
 279 //            8
 280 */
 281 static int QuadPenta_F [5][9] = {  // FORWARD
 282   { 0, 6, 1, 7, 2, 8, 0, 0, 0 },
 283   { 3, 11,5, 10,4, 9, 3, 3, 3 },
 284   { 0, 12,3, 9, 4, 13,1, 6, 0 },
 285   { 1, 13,4, 10,5, 14,2, 7, 1 },
 286   { 0, 8, 2, 14,5, 11,3, 12,0 }};
 287 static int QuadPenta_RE [5][9] = { // REVERSED -> EXTERNAL
 288   { 0, 8, 2, 7, 1, 6, 0, 0, 0 },
 289   { 3, 9, 4, 10,5, 11,3, 3, 3 },
 290   { 0, 6, 1, 13,4, 9, 3, 12,0 },
 291   { 1, 7, 2, 14,5, 10,4, 13,1 },
 292   { 0, 12,3, 11,5, 14,2, 8, 0 }};
 293 static int QuadPenta_nbN [] = { 6, 6, 8, 8, 8 };
 294
 295 /*
 296 //                 13
 297 //         N5+-----+-----+N6                          +-----+-----+
 298 //          /|          /|                           /|          /|
 299 //       12+ |       14+ |                          + |   +25   + |
 300 //        /  |        /  |                         /  |        /  |
 301 //     N4+-----+-----+N7 |       QUADRATIC        +-----+-----+   |  Central nodes
 302 //       |   | 15    |   |       HEXAHEDRON       |   |       |   |  of tri-quadratic
 303 //       |   |       |   |                        |   |       |   |  HEXAHEDRON
 304 //       | 17+       |   +18                      |   +   22+ |   +
 305 //       |   |       |   |                        |21 |       |   |
 306 //       |   |       |   |                        | + | 26+   | + |
 307 //       |   |       |   |                        |   |       |23 |
 308 //     16+   |       +19 |                        +   | +24   +   |
 309 //       |   |       |   |                        |   |       |   |
 310 //       |   |     9 |   |                        |   |       |   |
 311 //       | N1+-----+-|---+N2                      |   +-----+-|---+
 312 //       |  /        |  /                         |  /        |  /
 313 //       | +8        | +10                        | +   20+   | +
 314 //       |/          |/                           |/          |/
 315 //     N0+-----+-----+N3                          +-----+-----+
 316 //             11
 317 */
 318 static int QuadHexa_F [6][9] = {  // FORWARD
 319   { 0, 8, 1, 9, 2, 10,3, 11,0 },   // all face normals are external,
 320   { 4, 15,7, 14,6, 13,5, 12,4 },
 321   { 0, 16,4, 12,5, 17,1, 8, 0 },
 322   { 1, 17,5, 13,6, 18,2, 9, 1 },
 323   { 3, 10,2, 18,6, 14,7, 19,3 },
 324   { 0, 11,3, 19,7, 15,4, 16,0 }};
 325 static int QuadHexa_RE [6][9] = {  // REVERSED -> EXTERNAL
 326   { 0, 11,3, 10,2, 9, 1, 8, 0 },   // all face normals are external
 327   { 4, 12,5, 13,6, 14,7, 15,4 },
 328   { 0, 8, 1, 17,5, 12,4, 16,0 },
 329   { 1, 9, 2, 18,6, 13,5, 17,1 },
 330   { 3, 19,7, 14,6, 18,2, 10,3 },
 331   { 0, 16,4, 15,7, 19,3, 11,0 }};
 332 static int QuadHexa_nbN [] = { 8, 8, 8, 8, 8, 8 };
 333
 334 static int TriQuadHexa_F [6][9] = {  // FORWARD
 335   { 0, 8, 1, 9, 2, 10,3, 11, 20 },   // all face normals are external
 336   { 4, 15,7, 14,6, 13,5, 12, 25 },
 337   { 0, 16,4, 12,5, 17,1, 8,  21 },
 338   { 1, 17,5, 13,6, 18,2, 9,  22 },
 339   { 3, 10,2, 18,6, 14,7, 19, 23 },
 340   { 0, 11,3, 19,7, 15,4, 16, 24 }};
 341 static int TriQuadHexa_RE [6][9] = {  // REVERSED -> EXTERNAL
 342   { 0, 11,3, 10,2, 9, 1, 8,  20 },   // opposite faces are neighbouring,
 343   { 4, 12,5, 13,6, 14,7, 15, 25 },   // all face normals are external
 344   { 0, 8, 1, 17,5, 12,4, 16, 21 },
 345   { 1, 9, 2, 18,6, 13,5, 17, 22 },
 346   { 3, 19,7, 14,6, 18,2, 10, 23 },
 347   { 0, 16,4, 15,7, 19,3, 11, 24 }};
 348 static int TriQuadHexa_nbN [] = { 9, 9, 9, 9, 9, 9 };
 349
 350
 351 // ========================================================
 352 // to perform some calculations without linkage to CASCADE
 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  *        This check is valid for conformal meshes only
1502  */
1503 //================================================================================
1504
1505 bool SMDS_VolumeTool::IsFreeFace( int faceIndex, const SMDS_MeshElement** otherVol/*=0*/ ) const
1506 {
1507   const bool isFree = true;
1508
1509   if (!setFace( faceIndex ))
1510     return !isFree;
1511
1512   const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
1513
1514   const int  di = myVolume->IsQuadratic() ? 2 : 1;
1515   const int nbN = ( myFaceNbNodes/di <= 4 && !IsPoly()) ? 3 : myFaceNbNodes/di; // nb nodes to check
1516
1517   SMDS_ElemIteratorPtr eIt = nodes[0]->GetInverseElementIterator( SMDSAbs_Volume );
1518   while ( eIt->more() )
1519   {
1520     const SMDS_MeshElement* vol = eIt->next();
1521     if ( vol == myVolume )
1522       continue;
1523     int iN;
1524     for ( iN = 1; iN < nbN; ++iN )
1525       if ( vol->GetNodeIndex( nodes[ iN*di ]) < 0 )
1526         break;
1527     if ( iN == nbN ) // nbN nodes are shared with vol
1528     {
1529       // if ( vol->IsPoly() || vol->NbFaces() > 6 ) // vol is polyhed or hex prism
1530       // {
1531       //   int nb = myFaceNbNodes;
1532       //   if ( myVolume->GetEntityType() != vol->GetEntityType() )
1533       //     nb -= ( GetCenterNodeIndex(0) > 0 );
1534       //   set<const SMDS_MeshNode*> faceNodes( nodes, nodes + nb );
1535       //   if ( SMDS_VolumeTool( vol ).GetFaceIndex( faceNodes ) < 0 )
1536       //     continue;
1537       // }
1538       if ( otherVol ) *otherVol = vol;
1539       return !isFree;
1540     }
1541   }
1542   if ( otherVol ) *otherVol = 0;
1543   return isFree;
1544 }
1545
1546 //================================================================================
1547 /*!
1548  * \brief Thorough check that only one volume is build on the face nodes
1549  */
1550 //================================================================================
1551
1552 bool SMDS_VolumeTool::IsFreeFaceAdv( int faceIndex, const SMDS_MeshElement** otherVol/*=0*/ ) const
1553 {
1554   const bool isFree = true;
1555
1556   if (!setFace( faceIndex ))
1557     return !isFree;
1558
1559   const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
1560   const int nbFaceNodes = myFaceNbNodes;
1561
1562   // evaluate nb of face nodes shared by other volumes
1563   int maxNbShared = -1;
1564   typedef map< const SMDS_MeshElement*, int > TElemIntMap;
1565   TElemIntMap volNbShared;
1566   TElemIntMap::iterator vNbIt;
1567   for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) {
1568     const SMDS_MeshNode* n = nodes[ iNode ];
1569     SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( SMDSAbs_Volume );
1570     while ( eIt->more() ) {
1571       const SMDS_MeshElement* elem = eIt->next();
1572       if ( elem != myVolume ) {
1573         vNbIt = volNbShared.insert( make_pair( elem, 0 )).first;
1574         (*vNbIt).second++;
1575         if ( vNbIt->second > maxNbShared )
1576           maxNbShared = vNbIt->second;
1577       }
1578     }
1579   }
1580   if ( maxNbShared < 3 )
1581     return isFree; // is free
1582
1583   // find volumes laying on the opposite side of the face
1584   // and sharing all nodes
1585   XYZ intNormal; // internal normal
1586   GetFaceNormal( faceIndex, intNormal.x, intNormal.y, intNormal.z );
1587   if ( IsFaceExternal( faceIndex ))
1588     intNormal = XYZ( -intNormal.x, -intNormal.y, -intNormal.z );
1589   XYZ p0 ( nodes[0] ), baryCenter;
1590   for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end();  ) {
1591     const int& nbShared = (*vNbIt).second;
1592     if ( nbShared >= 3 ) {
1593       SMDS_VolumeTool volume( (*vNbIt).first );
1594       volume.GetBaryCenter( baryCenter.x, baryCenter.y, baryCenter.z );
1595       XYZ intNormal2( baryCenter - p0 );
1596       if ( intNormal.Dot( intNormal2 ) < 0 ) {
1597         // opposite side
1598         if ( nbShared >= nbFaceNodes )
1599         {
1600           // a volume shares the whole facet
1601           if ( otherVol ) *otherVol = vNbIt->first;
1602           return !isFree;
1603         }
1604         ++vNbIt;
1605         continue;
1606       }
1607     }
1608     // remove a volume from volNbShared map
1609     volNbShared.erase( vNbIt++ );
1610   }
1611
1612   // here volNbShared contains only volumes laying on the opposite side of
1613   // the face and sharing 3 or more but not all face nodes with myVolume
1614   if ( volNbShared.size() < 2 ) {
1615     return isFree; // is free
1616   }
1617
1618   // check if the whole area of a face is shared
1619   for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
1620   {
1621     const SMDS_MeshNode* n = nodes[ iNode ];
1622     // check if n is shared by one of volumes of volNbShared
1623     bool isShared = false;
1624     SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( SMDSAbs_Volume );
1625     while ( eIt->more() && !isShared )
1626       isShared = volNbShared.count( eIt->next() );
1627     if ( !isShared )
1628       return isFree;
1629   }
1630   if ( otherVol ) *otherVol = volNbShared.begin()->first;
1631   return !isFree;
1632
1633 //   if ( !myVolume->IsPoly() )
1634 //   {
1635 //     bool isShared[] = { false, false, false, false }; // 4 triangle parts of a quadrangle
1636 //     for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); vNbIt++ ) {
1637 //       SMDS_VolumeTool volume( (*vNbIt).first );
1638 //       bool prevLinkShared = false;
1639 //       int nbSharedLinks = 0;
1640 //       for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) {
1641 //         bool linkShared = volume.IsLinked( nodes[ iNode ], nodes[ iNode + 1] );
1642 //         if ( linkShared )
1643 //           nbSharedLinks++;
1644 //         if ( linkShared && prevLinkShared &&
1645 //              volume.IsLinked( nodes[ iNode - 1 ], nodes[ iNode + 1] ))
1646 //           isShared[ iNode ] = true;
1647 //         prevLinkShared = linkShared;
1648 //       }
1649 //       if ( nbSharedLinks == nbFaceNodes )
1650 //         return !free; // is not free
1651 //       if ( nbFaceNodes == 4 ) {
1652 //         // check traingle parts 1 & 3
1653 //         if ( isShared[1] && isShared[3] )
1654 //           return !free; // is not free
1655 //         // check triangle parts 0 & 2;
1656 //         // 0 part could not be checked in the loop; check it here
1657 //         if ( isShared[2] && prevLinkShared &&
1658 //              volume.IsLinked( nodes[ 0 ], nodes[ 1 ] ) &&
1659 //              volume.IsLinked( nodes[ 1 ], nodes[ 3 ] ) )
1660 //           return !free; // is not free
1661 //       }
1662 //     }
1663 //   }
1664 //  return free;
1665 }
1666
1667 //=======================================================================
1668 //function : GetFaceIndex
1669 //purpose  : Return index of a face formed by theFaceNodes
1670 //=======================================================================
1671
1672 int SMDS_VolumeTool::GetFaceIndex( const set<const SMDS_MeshNode*>& theFaceNodes,
1673                                    const int                        theFaceIndexHint ) const
1674 {
1675   if ( theFaceIndexHint >= 0 )
1676   {
1677     int nbNodes = NbFaceNodes( theFaceIndexHint );
1678     if ( nbNodes == (int) theFaceNodes.size() )
1679     {
1680       const SMDS_MeshNode** nodes = GetFaceNodes( theFaceIndexHint );
1681       while ( nbNodes )
1682         if ( theFaceNodes.count( nodes[ nbNodes-1 ]))
1683           --nbNodes;
1684         else
1685           break;
1686       if ( nbNodes == 0 )
1687         return theFaceIndexHint;
1688     }
1689   }
1690   for ( int iFace = 0; iFace < myNbFaces; iFace++ )
1691   {
1692     if ( iFace == theFaceIndexHint )
1693       continue;
1694     int nbNodes = NbFaceNodes( iFace );
1695     if ( nbNodes == (int) theFaceNodes.size() )
1696     {
1697       const SMDS_MeshNode** nodes = GetFaceNodes( iFace );
1698       while ( nbNodes )
1699         if ( theFaceNodes.count( nodes[ nbNodes-1 ]))
1700           --nbNodes;
1701         else
1702           break;
1703       if ( nbNodes == 0 )
1704         return iFace;
1705     }
1706   }
1707   return -1;
1708 }
1709
1710 //=======================================================================
1711 //function : GetFaceIndex
1712 //purpose  : Return index of a face formed by theFaceNodes
1713 //=======================================================================
1714
1715 /*int SMDS_VolumeTool::GetFaceIndex( const set<int>& theFaceNodesIndices )
1716 {
1717   for ( int iFace = 0; iFace < myNbFaces; iFace++ ) {
1718     const int* nodes = GetFaceNodesIndices( iFace );
1719     int nbFaceNodes = NbFaceNodes( iFace );
1720     set<int> nodeSet;
1721     for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
1722       nodeSet.insert( nodes[ iNode ] );
1723     if ( theFaceNodesIndices == nodeSet )
1724       return iFace;
1725   }
1726   return -1;
1727 }*/
1728
1729 //=======================================================================
1730 //function : setFace
1731 //purpose  :
1732 //=======================================================================
1733
1734 bool SMDS_VolumeTool::setFace( int faceIndex ) const
1735 {
1736   if ( !myVolume )
1737     return false;
1738
1739   if ( myCurFace == faceIndex )
1740     return true;
1741
1742   myCurFace = -1;
1743
1744   if ( faceIndex < 0 || faceIndex >= NbFaces() )
1745     return false;
1746
1747   if (myFaceNodes != NULL) {
1748     delete [] myFaceNodes;
1749     myFaceNodes = NULL;
1750   }
1751
1752   if (myVolume->IsPoly())
1753   {
1754     if (!myPolyedre) {
1755       MESSAGE("Warning: bad volumic element");
1756       return false;
1757     }
1758
1759     // set face nodes
1760     int iNode;
1761     myFaceNbNodes = myPolyedre->NbFaceNodes(faceIndex + 1);
1762     myFaceNodes = new const SMDS_MeshNode* [myFaceNbNodes + 1];
1763     for ( iNode = 0; iNode < myFaceNbNodes; iNode++ )
1764       myFaceNodes[ iNode ] = myPolyedre->GetFaceNode(faceIndex + 1, iNode + 1);
1765     myFaceNodes[ myFaceNbNodes ] = myFaceNodes[ 0 ]; // last = first
1766
1767     // check orientation
1768     if (myExternalFaces)
1769     {
1770       myCurFace = faceIndex; // avoid infinite recursion in IsFaceExternal()
1771       myExternalFaces = false; // force normal computation by IsFaceExternal()
1772       if ( !IsFaceExternal( faceIndex ))
1773         for ( int i = 0, j = myFaceNbNodes; i < j; ++i, --j )
1774           std::swap( myFaceNodes[i], myFaceNodes[j] );
1775       myExternalFaces = true;
1776     }
1777   }
1778   else
1779   {
1780     if ( !myAllFacesNodeIndices_F )
1781     {
1782       // choose data for an element type
1783       switch ( myVolumeNbNodes ) {
1784       case 4:
1785         myAllFacesNodeIndices_F  = &Tetra_F [0][0];
1786         //myAllFacesNodeIndices_FE = &Tetra_F [0][0];
1787         myAllFacesNodeIndices_RE = &Tetra_RE[0][0];
1788         myAllFacesNbNodes        = Tetra_nbN;
1789         myMaxFaceNbNodes         = sizeof(Tetra_F[0])/sizeof(Tetra_F[0][0]);
1790         break;
1791       case 5:
1792         myAllFacesNodeIndices_F  = &Pyramid_F [0][0];
1793         //myAllFacesNodeIndices_FE = &Pyramid_F [0][0];
1794         myAllFacesNodeIndices_RE = &Pyramid_RE[0][0];
1795         myAllFacesNbNodes        = Pyramid_nbN;
1796         myMaxFaceNbNodes         = sizeof(Pyramid_F[0])/sizeof(Pyramid_F[0][0]);
1797         break;
1798       case 6:
1799         myAllFacesNodeIndices_F  = &Penta_F [0][0];
1800         //myAllFacesNodeIndices_FE = &Penta_FE[0][0];
1801         myAllFacesNodeIndices_RE = &Penta_RE[0][0];
1802         myAllFacesNbNodes        = Penta_nbN;
1803         myMaxFaceNbNodes         = sizeof(Penta_F[0])/sizeof(Penta_F[0][0]);
1804         break;
1805       case 8:
1806         myAllFacesNodeIndices_F  = &Hexa_F [0][0];
1807         ///myAllFacesNodeIndices_FE = &Hexa_FE[0][0];
1808         myAllFacesNodeIndices_RE = &Hexa_RE[0][0];
1809         myAllFacesNbNodes        = Hexa_nbN;
1810         myMaxFaceNbNodes         = sizeof(Hexa_F[0])/sizeof(Hexa_F[0][0]);
1811         break;
1812       case 10:
1813         myAllFacesNodeIndices_F  = &QuadTetra_F [0][0];
1814         //myAllFacesNodeIndices_FE = &QuadTetra_F [0][0];
1815         myAllFacesNodeIndices_RE = &QuadTetra_RE[0][0];
1816         myAllFacesNbNodes        = QuadTetra_nbN;
1817         myMaxFaceNbNodes         = sizeof(QuadTetra_F[0])/sizeof(QuadTetra_F[0][0]);
1818         break;
1819       case 13:
1820         myAllFacesNodeIndices_F  = &QuadPyram_F [0][0];
1821         //myAllFacesNodeIndices_FE = &QuadPyram_F [0][0];
1822         myAllFacesNodeIndices_RE = &QuadPyram_RE[0][0];
1823         myAllFacesNbNodes        = QuadPyram_nbN;
1824         myMaxFaceNbNodes         = sizeof(QuadPyram_F[0])/sizeof(QuadPyram_F[0][0]);
1825         break;
1826       case 15:
1827         myAllFacesNodeIndices_F  = &QuadPenta_F [0][0];
1828         //myAllFacesNodeIndices_FE = &QuadPenta_FE[0][0];
1829         myAllFacesNodeIndices_RE = &QuadPenta_RE[0][0];
1830         myAllFacesNbNodes        = QuadPenta_nbN;
1831         myMaxFaceNbNodes         = sizeof(QuadPenta_F[0])/sizeof(QuadPenta_F[0][0]);
1832         break;
1833       case 20:
1834       case 27:
1835         myAllFacesNodeIndices_F  = &QuadHexa_F [0][0];
1836         //myAllFacesNodeIndices_FE = &QuadHexa_FE[0][0];
1837         myAllFacesNodeIndices_RE = &QuadHexa_RE[0][0];
1838         myAllFacesNbNodes        = QuadHexa_nbN;
1839         myMaxFaceNbNodes         = sizeof(QuadHexa_F[0])/sizeof(QuadHexa_F[0][0]);
1840         if ( !myIgnoreCentralNodes && myVolumeNbNodes == 27 )
1841         {
1842           myAllFacesNodeIndices_F  = &TriQuadHexa_F [0][0];
1843           //myAllFacesNodeIndices_FE = &TriQuadHexa_FE[0][0];
1844           myAllFacesNodeIndices_RE = &TriQuadHexa_RE[0][0];
1845           myAllFacesNbNodes        = TriQuadHexa_nbN;
1846           myMaxFaceNbNodes         = sizeof(TriQuadHexa_F[0])/sizeof(TriQuadHexa_F[0][0]);
1847         }
1848         break;
1849       case 12:
1850         myAllFacesNodeIndices_F  = &HexPrism_F [0][0];
1851         //myAllFacesNodeIndices_FE = &HexPrism_FE[0][0];
1852         myAllFacesNodeIndices_RE = &HexPrism_RE[0][0];
1853         myAllFacesNbNodes        = HexPrism_nbN;
1854         myMaxFaceNbNodes         = sizeof(HexPrism_F[0])/sizeof(HexPrism_F[0][0]);
1855         break;
1856       default:
1857         return false;
1858       }
1859     }
1860     myFaceNbNodes = myAllFacesNbNodes[ faceIndex ];
1861     // if ( myExternalFaces )
1862     //   myFaceNodeIndices = (int*)( myVolForward ? myAllFacesNodeIndices_FE + faceIndex*myMaxFaceNbNodes : myAllFacesNodeIndices_RE + faceIndex*myMaxFaceNbNodes );
1863     // else
1864     //   myFaceNodeIndices = (int*)( myAllFacesNodeIndices_F + faceIndex*myMaxFaceNbNodes );
1865     myFaceNodeIndices = (int*)( myVolForward ? myAllFacesNodeIndices_F + faceIndex*myMaxFaceNbNodes : myAllFacesNodeIndices_RE + faceIndex*myMaxFaceNbNodes );
1866
1867     // set face nodes
1868     myFaceNodes = new const SMDS_MeshNode* [myFaceNbNodes + 1];
1869     for ( int iNode = 0; iNode < myFaceNbNodes; iNode++ )
1870       myFaceNodes[ iNode ] = myVolumeNodes[ myFaceNodeIndices[ iNode ]];
1871     myFaceNodes[ myFaceNbNodes ] = myFaceNodes[ 0 ];
1872   }
1873
1874   myCurFace = faceIndex;
1875
1876   return true;
1877 }
1878
1879 //=======================================================================
1880 //function : GetType
1881 //purpose  : return VolumeType by nb of nodes in a volume
1882 //=======================================================================
1883
1884 SMDS_VolumeTool::VolumeType SMDS_VolumeTool::GetType(int nbNodes)
1885 {
1886   switch ( nbNodes ) {
1887   case 4: return TETRA;
1888   case 5: return PYRAM;
1889   case 6: return PENTA;
1890   case 8: return HEXA;
1891   case 10: return QUAD_TETRA;
1892   case 13: return QUAD_PYRAM;
1893   case 15: return QUAD_PENTA;
1894   case 20:
1895   case 27: return QUAD_HEXA;
1896   case 12: return HEX_PRISM;
1897   default:return UNKNOWN;
1898   }
1899 }
1900
1901 //=======================================================================
1902 //function : NbFaces
1903 //purpose  : return nb of faces by volume type
1904 //=======================================================================
1905
1906 int SMDS_VolumeTool::NbFaces( VolumeType type )
1907 {
1908   switch ( type ) {
1909   case TETRA     :
1910   case QUAD_TETRA: return 4;
1911   case PYRAM     :
1912   case QUAD_PYRAM: return 5;
1913   case PENTA     :
1914   case QUAD_PENTA: return 5;
1915   case HEXA      :
1916   case QUAD_HEXA : return 6;
1917   case HEX_PRISM : return 8;
1918   default:         return 0;
1919   }
1920 }
1921
1922 //================================================================================
1923 /*!
1924  * \brief Useful to know nb of corner nodes of a quadratic volume
1925   * \param type - volume type
1926   * \retval int - nb of corner nodes
1927  */
1928 //================================================================================
1929
1930 int SMDS_VolumeTool::NbCornerNodes(VolumeType type)
1931 {
1932   switch ( type ) {
1933   case TETRA     :
1934   case QUAD_TETRA: return 4;
1935   case PYRAM     :
1936   case QUAD_PYRAM: return 5;
1937   case PENTA     :
1938   case QUAD_PENTA: return 6;
1939   case HEXA      :
1940   case QUAD_HEXA : return 8;
1941   case HEX_PRISM : return 12;
1942   default:         return 0;
1943   }
1944   return 0;
1945 }
1946   //
1947
1948 //=======================================================================
1949 //function : GetFaceNodesIndices
1950 //purpose  : Return the array of face nodes indices
1951 //           To comfort link iteration, the array
1952 //           length == NbFaceNodes( faceIndex ) + 1 and
1953 //           the last node index == the first one.
1954 //=======================================================================
1955
1956 const int* SMDS_VolumeTool::GetFaceNodesIndices(VolumeType type,
1957                                                 int        faceIndex,
1958                                                 bool       external)
1959 {
1960   switch ( type ) {
1961   case TETRA: return Tetra_F[ faceIndex ];
1962   case PYRAM: return Pyramid_F[ faceIndex ];
1963   case PENTA: return external ? Penta_F[ faceIndex ] : Penta_F[ faceIndex ];
1964   case HEXA:  return external ? Hexa_F[ faceIndex ] : Hexa_F[ faceIndex ];
1965   case QUAD_TETRA: return QuadTetra_F[ faceIndex ];
1966   case QUAD_PYRAM: return QuadPyram_F[ faceIndex ];
1967   case QUAD_PENTA: return external ? QuadPenta_F[ faceIndex ] : QuadPenta_F[ faceIndex ];
1968     // what about SMDSEntity_TriQuad_Hexa?
1969   case QUAD_HEXA:  return external ? QuadHexa_F[ faceIndex ] : QuadHexa_F[ faceIndex ];
1970   case HEX_PRISM:  return external ? HexPrism_F[ faceIndex ] : HexPrism_F[ faceIndex ];
1971   default:;
1972   }
1973   return 0;
1974 }
1975
1976 //=======================================================================
1977 //function : NbFaceNodes
1978 //purpose  : Return number of nodes in the array of face nodes
1979 //=======================================================================
1980
1981 int SMDS_VolumeTool::NbFaceNodes(VolumeType type,
1982                                  int        faceIndex )
1983 {
1984   switch ( type ) {
1985   case TETRA: return Tetra_nbN[ faceIndex ];
1986   case PYRAM: return Pyramid_nbN[ faceIndex ];
1987   case PENTA: return Penta_nbN[ faceIndex ];
1988   case HEXA:  return Hexa_nbN[ faceIndex ];
1989   case QUAD_TETRA: return QuadTetra_nbN[ faceIndex ];
1990   case QUAD_PYRAM: return QuadPyram_nbN[ faceIndex ];
1991   case QUAD_PENTA: return QuadPenta_nbN[ faceIndex ];
1992     // what about SMDSEntity_TriQuad_Hexa?
1993   case QUAD_HEXA:  return QuadHexa_nbN[ faceIndex ];
1994   case HEX_PRISM:  return HexPrism_nbN[ faceIndex ];
1995   default:;
1996   }
1997   return 0;
1998 }
1999
2000 //=======================================================================
2001 //function : Element
2002 //purpose  : return element
2003 //=======================================================================
2004
2005 const SMDS_MeshVolume* SMDS_VolumeTool::Element() const
2006 {
2007   return static_cast<const SMDS_MeshVolume*>( myVolume );
2008 }
2009
2010 //=======================================================================
2011 //function : ID
2012 //purpose  : return element ID
2013 //=======================================================================
2014
2015 int SMDS_VolumeTool::ID() const
2016 {
2017   return myVolume ? myVolume->GetID() : 0;
2018 }