src/SMDS/SMDS_VolumeTool.cxx

   1 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 // File      : SMDS_VolumeTool.cxx
  24 // Created   : Tue Jul 13 12:22:13 2004
  25 // Author    : Edward AGAPOV (eap)
  26 //
  27 #ifdef _MSC_VER
  28 #pragma warning(disable:4786)
  29 #endif
  30
  31 #include "SMDS_VolumeTool.hxx"
  32
  33 #include "SMDS_MeshElement.hxx"
  34 #include "SMDS_MeshNode.hxx"
  35 #include "SMDS_VtkVolume.hxx"
  36 #include "SMDS_Mesh.hxx"
  37
  38 #include "utilities.h"
  39
  40 #include <map>
  41 #include <limits>
  42 #include <cmath>
  43
  44 using namespace std;
  45
  46 // ======================================================
  47 // Node indices in faces depending on volume orientation
  48 // making most faces normals external
  49 // ======================================================
  50 // For all elements, 0-th face is bottom based on the first nodes.
  51 // For prismatic elements (tetra,hexa,prisms), 1-th face is a top one.
  52 // For all elements, side faces follow order of bottom nodes
  53 // ======================================================
  54
  55 /*
  56 //           N3
  57 //           +
  58 //          /|\
  59 //         / | \
  60 //        /  |  \
  61 //    N0 +---|---+ N1                TETRAHEDRON
  62 //       \   |   /
  63 //        \  |  /
  64 //         \ | /
  65 //          \|/
  66 //           +
  67 //           N2
  68 */
  69 static int Tetra_F [4][4] = { // FORWARD == EXTERNAL
  70   { 0, 1, 2, 0 },              // All faces have external normals
  71   { 0, 3, 1, 0 },
  72   { 1, 3, 2, 1 },
  73   { 0, 2, 3, 0 }};
  74 static int Tetra_RE [4][4] = { // REVERSED -> FORWARD (EXTERNAL)
  75   { 0, 2, 1, 0 },              // All faces have external normals
  76   { 0, 1, 3, 0 },
  77   { 1, 2, 3, 1 },
  78   { 0, 3, 2, 0 }};
  79 static int Tetra_nbN [] = { 3, 3, 3, 3 };
  80
  81 //
  82 //     PYRAMID
  83 //
  84 static int Pyramid_F [5][5] = { // FORWARD == EXTERNAL
  85   { 0, 1, 2, 3, 0 },            // All faces have external normals
  86   { 0, 4, 1, 0, 4 },
  87   { 1, 4, 2, 1, 4 },
  88   { 2, 4, 3, 2, 4 },
  89   { 3, 4, 0, 3, 4 }
  90 };
  91 static int Pyramid_RE [5][5] = { // REVERSED -> FORWARD (EXTERNAL)
  92   { 0, 3, 2, 1, 0 },             // All faces but a bottom have external normals
  93   { 0, 1, 4, 0, 4 },
  94   { 1, 2, 4, 1, 4 },
  95   { 2, 3, 4, 2, 4 },
  96   { 3, 0, 4, 3, 4 }};
  97 static int Pyramid_nbN [] = { 4, 3, 3, 3, 3 };
  98
  99 /*
 100 //            + N4
 101 //           /|\
 102 //          / | \
 103 //         /  |  \
 104 //        /   |   \
 105 //    N3 +---------+ N5
 106 //       |    |    |
 107 //       |    + N1 |
 108 //       |   / \   |                PENTAHEDRON
 109 //       |  /   \  |
 110 //       | /     \ |
 111 //       |/       \|
 112 //    N0 +---------+ N2
 113 */
 114 static int Penta_F [5][5] = { // FORWARD
 115   { 0, 1, 2, 0, 0 },          // All faces have external normals
 116   { 3, 5, 4, 3, 3 },          // 0 is bottom, 1 is top face
 117   { 0, 3, 4, 1, 0 },
 118   { 1, 4, 5, 2, 1 },
 119   { 0, 2, 5, 3, 0 }};
 120 static int Penta_RE [5][5] = { // REVERSED -> EXTERNAL
 121   { 0, 2, 1, 0, 0 },
 122   { 3, 4, 5, 3, 3 },
 123   { 0, 1, 4, 3, 0 },
 124   { 1, 2, 5, 4, 1 },
 125   { 0, 3, 5, 2, 0 }};
 126 static int Penta_nbN [] = { 3, 3, 4, 4, 4 };
 127
 128 /*
 129 //         N5+----------+N6
 130 //          /|         /|
 131 //         / |        / |
 132 //        /  |       /  |
 133 //     N4+----------+N7 |
 134 //       |   |      |   |           HEXAHEDRON
 135 //       |   |      |   |
 136 //       |   |      |   |
 137 //       | N1+------|---+N2
 138 //       |  /       |  /
 139 //       | /        | /
 140 //       |/         |/
 141 //     N0+----------+N3
 142 */
 143 static int Hexa_F [6][5] = { // FORWARD
 144   { 0, 1, 2, 3, 0 },
 145   { 4, 7, 6, 5, 4 },          // all face normals are external
 146   { 0, 4, 5, 1, 0 },
 147   { 1, 5, 6, 2, 1 },
 148   { 3, 2, 6, 7, 3 },
 149   { 0, 3, 7, 4, 0 }};
 150 static int Hexa_RE [6][5] = { // REVERSED -> EXTERNAL
 151   { 0, 3, 2, 1, 0 },
 152   { 4, 5, 6, 7, 4 },          // all face normals are external
 153   { 0, 1, 5, 4, 0 },
 154   { 1, 2, 6, 5, 1 },
 155   { 3, 7, 6, 2, 3 },
 156   { 0, 4, 7, 3, 0 }};
 157 static int Hexa_nbN [] = { 4, 4, 4, 4, 4, 4 };
 158
 159 /*
 160 //      N8 +------+ N9
 161 //        /        \
 162 //       /          \
 163 //   N7 +            + N10
 164 //       \          /
 165 //        \        /
 166 //      N6 +------+ N11
 167 //                             HEXAGONAL PRISM
 168 //      N2 +------+ N3
 169 //        /        \
 170 //       /          \
 171 //   N1 +            + N4
 172 //       \          /
 173 //        \        /
 174 //      N0 +------+ N5
 175 */
 176 static int HexPrism_F [8][7] = { // FORWARD
 177   { 0, 1, 2, 3, 4, 5, 0 },
 178   { 6,11,10, 9, 8, 7, 6 },
 179   { 0, 6, 7, 1, 0, 0, 0 },
 180   { 1, 7, 8, 2, 1, 1, 1 },
 181   { 2, 8, 9, 3, 2, 2, 2 },
 182   { 3, 9,10, 4, 3, 3, 3 },
 183   { 4,10,11, 5, 4, 4, 4 },
 184   { 5,11, 6, 0, 5, 5, 5 }};
 185 static int HexPrism_RE [8][7] = { // REVERSED -> EXTERNAL
 186   { 0, 5, 4, 3, 2, 1, 0 },
 187   { 6,11,10, 9, 8, 7, 6 },
 188   { 0, 6, 7, 1, 0, 0, 0 },
 189   { 1, 7, 8, 2, 1, 1, 1 },
 190   { 2, 8, 9, 3, 2, 2, 2 },
 191   { 3, 9,10, 4, 3, 3, 3 },
 192   { 4,10,11, 5, 4, 4, 4 },
 193   { 5,11, 6, 0, 5, 5, 5 }};
 194 static int HexPrism_nbN [] = { 6, 6, 4, 4, 4, 4, 4, 4 };
 195
 196
 197 /*
 198 //           N3
 199 //           +
 200 //          /|\
 201 //        7/ | \8
 202 //        /  |4 \                    QUADRATIC
 203 //    N0 +---|---+ N1                TETRAHEDRON
 204 //       \   +9  /
 205 //        \  |  /
 206 //        6\ | /5
 207 //          \|/
 208 //           +
 209 //           N2
 210 */
 211 static int QuadTetra_F [4][7] = { // FORWARD
 212   { 0, 4, 1, 5, 2, 6, 0 },        // All faces have external normals
 213   { 0, 7, 3, 8, 1, 4, 0 },
 214   { 1, 8, 3, 9, 2, 5, 1 },
 215   { 0, 6, 2, 9, 3, 7, 0 }};
 216 static int QuadTetra_RE [4][7] = { // REVERSED -> FORWARD (EXTERNAL)
 217   { 0, 6, 2, 5, 1, 4, 0 },         // All faces have external normals
 218   { 0, 4, 1, 8, 3, 7, 0 },
 219   { 1, 5, 2, 9, 3, 8, 1 },
 220   { 0, 7, 3, 9, 2, 6, 0 }};
 221 static int QuadTetra_nbN [] = { 6, 6, 6, 6 };
 222
 223 //
 224 //     QUADRATIC
 225 //     PYRAMID
 226 //
 227 //            +4
 228 //
 229 //
 230 //       10+-----+11
 231 //         |     |        9 - middle point for (0,4) etc.
 232 //         |     |
 233 //        9+-----+12
 234 //
 235 //            6
 236 //      1+----+----+2
 237 //       |         |
 238 //       |         |
 239 //      5+         +7
 240 //       |         |
 241 //       |         |
 242 //      0+----+----+3
 243 //            8
 244 static int QuadPyram_F [5][9] = {  // FORWARD
 245   { 0, 5, 1, 6, 2, 7, 3, 8, 0 },   // All faces have external normals
 246   { 0, 9, 4, 10,1, 5, 0, 4, 4 },
 247   { 1, 10,4, 11,2, 6, 1, 4, 4 },
 248   { 2, 11,4, 12,3, 7, 2, 4, 4 },
 249   { 3, 12,4, 9, 0, 8, 3, 4, 4 }};
 250 static int QuadPyram_RE [5][9] = { // REVERSED -> FORWARD (EXTERNAL)
 251   { 0, 8, 3, 7, 2, 6, 1, 5, 0 },   // All faces but a bottom have external normals
 252   { 0, 5, 1, 10,4, 9, 0, 4, 4 },
 253   { 1, 6, 2, 11,4, 10,1, 4, 4 },
 254   { 2, 7, 3, 12,4, 11,2, 4, 4 },
 255   { 3, 8, 0, 9, 4, 12,3, 4, 4 }};
 256 static int QuadPyram_nbN [] = { 8, 6, 6, 6, 6 };
 257
 258 /*
 259 //            + N4
 260 //           /|\
 261 //         9/ | \10
 262 //         /  |  \
 263 //        /   |   \
 264 //    N3 +----+----+ N5
 265 //       |    |11  |
 266 //       |    |    |
 267 //       |    +13  |                QUADRATIC
 268 //       |    |    |                PENTAHEDRON
 269 //     12+    |    +14
 270 //       |    |    |
 271 //       |    |    |
 272 //       |    + N1 |
 273 //       |   / \   |
 274 //       | 6/   \7 |
 275 //       | /     \ |
 276 //       |/       \|
 277 //    N0 +---------+ N2
 278 //            8
 279 */
 280 static int QuadPenta_F [5][9] = {  // FORWARD
 281   { 0, 6, 1, 7, 2, 8, 0, 0, 0 },
 282   { 3,11, 5, 10,4, 9, 3, 3, 3 },
 283   { 0, 12,3, 9, 4, 13,1, 6, 0 },
 284   { 1, 13,4, 10,5, 14,2, 7, 1 },
 285   { 0, 8, 2, 14,5, 11,3, 12,0 }};
 286 static int QuadPenta_RE [5][9] = { // REVERSED -> EXTERNAL
 287   { 0, 8, 2, 7, 1, 6, 0, 0, 0 },
 288   { 3, 9, 4, 10,5, 11,3, 3, 3 },
 289   { 0, 6, 1, 13,4, 9, 3, 12,0 },
 290   { 1, 7, 2, 14,5, 10,4, 13,1 },
 291   { 0, 12,3, 11,5, 14,2, 8, 0 }};
 292 static int QuadPenta_nbN [] = { 6, 6, 8, 8, 8 };
 293
 294 /*
 295 //                 13
 296 //         N5+-----+-----+N6
 297 //          /|          /|
 298 //       12+ |       14+ |
 299 //        /  |        /  |
 300 //     N4+-----+-----+N7 |           QUADRATIC
 301 //       |   | 15    |   |           HEXAHEDRON
 302 //       |   |       |   |
 303 //       | 17+       |   +18
 304 //       |   |       |   |
 305 //       |   |       |   |
 306 //       |   |       |   |
 307 //     16+   |       +19 |
 308 //       |   |       |   |
 309 //       |   |     9 |   |
 310 //       | N1+-----+-|---+N2
 311 //       |  /        |  /
 312 //       | +8        | +10
 313 //       |/          |/
 314 //     N0+-----+-----+N3
 315 //             11
 316 */
 317 static int QuadHexa_F [6][9] = {  // FORWARD
 318   { 0, 8, 1, 9, 2, 10,3, 11,0 },   // all face normals are external,
 319   { 4, 15,7, 14,6, 13,5, 12,4 },
 320   { 0, 16,4, 12,5, 17,1, 8, 0 },
 321   { 1, 17,5, 13,6, 18,2, 9, 1 },
 322   { 3, 10,2, 18,6, 14,7, 19,3 },
 323   { 0, 11,3, 19,7, 15,4, 16,0 }};
 324 static int QuadHexa_RE [6][9] = {  // REVERSED -> EXTERNAL
 325   { 0, 11,3, 10,2, 9, 1, 8, 0 },   // all face normals are external
 326   { 4, 12,5, 13,6, 14,7, 15,4 },
 327   { 0, 8, 1, 17,5, 12,4, 16,0 },
 328   { 1, 9, 2, 18,6, 13,5, 17,1 },
 329   { 3, 19,7, 14,6, 18,2, 10,3 },
 330   { 0, 16,4, 15,7, 19,3, 11,0 }};
 331 static int QuadHexa_nbN [] = { 8, 8, 8, 8, 8, 8 };
 332
 333 static int TriQuadHexa_F [6][9] = {  // FORWARD
 334   { 0, 8, 1, 9, 2, 10,3, 11, 20 },   // all face normals are external
 335   { 4, 15,7, 14,6, 13,5, 12, 25 },
 336   { 0, 16,4, 12,5, 17,1, 8,  21 },
 337   { 1, 17,5, 13,6, 18,2, 9,  22 },
 338   { 3, 10,2, 18,6, 14,7, 19, 23 },
 339   { 0, 11,3, 19,7, 15,4, 16, 24 }};
 340 static int TriQuadHexa_RE [6][9] = {  // REVERSED -> EXTERNAL
 341   { 0, 11,3, 10,2, 9, 1, 8,  20 },   // opposite faces are neighbouring,
 342   { 4, 12,5, 13,6, 14,7, 15, 25 },   // all face normals are external
 343   { 0, 8, 1, 17,5, 12,4, 16, 21 },
 344   { 1, 9, 2, 18,6, 13,5, 17, 22 },
 345   { 3, 19,7, 14,6, 18,2, 10, 23 },
 346   { 0, 16,4, 15,7, 19,3, 11, 24 }};
 347 static int TriQuadHexa_nbN [] = { 9, 9, 9, 9, 9, 9 };
 348
 349
 350 // ========================================================
 351 // to perform some calculations without linkage to CASCADE
 352 // ========================================================
 353 namespace
 354 {
 355 struct XYZ {
 356   double x;
 357   double y;
 358   double z;
 359   XYZ()                               { x = 0; y = 0; z = 0; }
 360   XYZ( double X, double Y, double Z ) { x = X; y = Y; z = Z; }
 361   XYZ( const XYZ& other )             { x = other.x; y = other.y; z = other.z; }
 362   XYZ( const SMDS_MeshNode* n )       { x = n->X(); y = n->Y(); z = n->Z(); }
 363   inline XYZ operator-( const XYZ& other );
 364   inline XYZ operator+( const XYZ& other );
 365   inline XYZ Crossed( const XYZ& other );
 366   inline double Dot( const XYZ& other );
 367   inline double Magnitude();
 368   inline double SquareMagnitude();
 369 };
 370 inline XYZ XYZ::operator-( const XYZ& Right ) {
 371   return XYZ(x - Right.x, y - Right.y, z - Right.z);
 372 }
 373 inline XYZ XYZ::operator+( const XYZ& Right ) {
 374   return XYZ(x + Right.x, y + Right.y, z + Right.z);
 375 }
 376 inline XYZ XYZ::Crossed( const XYZ& Right ) {
 377   return XYZ (y * Right.z - z * Right.y,
 378               z * Right.x - x * Right.z,
 379               x * Right.y - y * Right.x);
 380 }
 381 inline double XYZ::Dot( const XYZ& Other ) {
 382   return(x * Other.x + y * Other.y + z * Other.z);
 383 }
 384 inline double XYZ::Magnitude() {
 385   return sqrt (x * x + y * y + z * z);
 386 }
 387 inline double XYZ::SquareMagnitude() {
 388   return (x * x + y * y + z * z);
 389 }
 390
 391   //================================================================================
 392   /*!
 393    * \brief Return linear type corresponding to a quadratic one
 394    */
 395   //================================================================================
 396
 397   SMDS_VolumeTool::VolumeType quadToLinear(SMDS_VolumeTool::VolumeType quadType)
 398   {
 399     SMDS_VolumeTool::VolumeType linType = SMDS_VolumeTool::VolumeType( int(quadType)-4 );
 400     const int nbCornersByQuad = SMDS_VolumeTool::NbCornerNodes( quadType );
 401     if ( SMDS_VolumeTool::NbCornerNodes( linType ) == nbCornersByQuad )
 402       return linType;
 403
 404     int iLin = 0;
 405     for ( ; iLin < SMDS_VolumeTool::NB_VOLUME_TYPES; ++iLin )
 406       if ( SMDS_VolumeTool::NbCornerNodes( SMDS_VolumeTool::VolumeType( iLin )) == nbCornersByQuad)
 407         return SMDS_VolumeTool::VolumeType( iLin );
 408
 409     return SMDS_VolumeTool::UNKNOWN;
 410   }
 411
 412 } // namespace
 413
 414 //=======================================================================
 415 //function : SMDS_VolumeTool
 416 //purpose  :
 417 //=======================================================================
 418
 419 SMDS_VolumeTool::SMDS_VolumeTool ()
 420   : myVolumeNodes( NULL ),
 421     myFaceNodes( NULL )
 422 {
 423   Set( 0 );
 424 }
 425
 426 //=======================================================================
 427 //function : SMDS_VolumeTool
 428 //purpose  :
 429 //=======================================================================
 430
 431 SMDS_VolumeTool::SMDS_VolumeTool (const SMDS_MeshElement* theVolume,
 432                                   const bool              ignoreCentralNodes)
 433   : myVolumeNodes( NULL ),
 434     myFaceNodes( NULL )
 435 {
 436   Set( theVolume, ignoreCentralNodes );
 437 }
 438
 439 //=======================================================================
 440 //function : SMDS_VolumeTool
 441 //purpose  :
 442 //=======================================================================
 443
 444 SMDS_VolumeTool::~SMDS_VolumeTool()
 445 {
 446   if ( myVolumeNodes != NULL ) delete [] myVolumeNodes;
 447   if ( myFaceNodes != NULL   ) delete [] myFaceNodes;
 448
 449   myFaceNodeIndices = NULL;
 450   myVolumeNodes = myFaceNodes = NULL;
 451 }
 452
 453 //=======================================================================
 454 //function : SetVolume
 455 //purpose  : Set volume to iterate on
 456 //=======================================================================
 457
 458 bool SMDS_VolumeTool::Set (const SMDS_MeshElement* theVolume,
 459                            const bool              ignoreCentralNodes)
 460 {
 461   // reset fields
 462   myVolume = 0;
 463   myPolyedre = 0;
 464   myIgnoreCentralNodes = ignoreCentralNodes;
 465
 466   myVolForward = true;
 467   myNbFaces = 0;
 468   myVolumeNbNodes = 0;
 469   if (myVolumeNodes != NULL) {
 470     delete [] myVolumeNodes;
 471     myVolumeNodes = NULL;
 472   }
 473   myPolyIndices.clear();
 474
 475   myExternalFaces = false;
 476
 477   myAllFacesNodeIndices_F  = 0;
 478   //myAllFacesNodeIndices_FE = 0;
 479   myAllFacesNodeIndices_RE = 0;
 480   myAllFacesNbNodes        = 0;
 481
 482   myCurFace = -1;
 483   myFaceNbNodes = 0;
 484   myFaceNodeIndices = NULL;
 485   if (myFaceNodes != NULL) {
 486     delete [] myFaceNodes;
 487     myFaceNodes = NULL;
 488   }
 489
 490   // set volume data
 491   if ( !theVolume || theVolume->GetType() != SMDSAbs_Volume )
 492     return false;
 493
 494   myVolume = theVolume;
 495   if (myVolume->IsPoly())
 496     myPolyedre = dynamic_cast<const SMDS_VtkVolume*>( myVolume );
 497
 498   myNbFaces = theVolume->NbFaces();
 499   myVolumeNbNodes = theVolume->NbNodes();
 500
 501   // set nodes
 502   int iNode = 0;
 503   myVolumeNodes = new const SMDS_MeshNode* [myVolumeNbNodes];
 504   SMDS_ElemIteratorPtr nodeIt = myVolume->nodesIterator();
 505   while ( nodeIt->more() )
 506     myVolumeNodes[ iNode++ ] = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
 507
 508   // check validity
 509   if ( !setFace(0) )
 510     return ( myVolume = 0 );
 511
 512   if ( !myPolyedre )
 513   {
 514     // define volume orientation
 515     XYZ botNormal;
 516     GetFaceNormal( 0, botNormal.x, botNormal.y, botNormal.z );
 517     const SMDS_MeshNode* botNode = myVolumeNodes[ 0 ];
 518     int topNodeIndex = myVolume->NbCornerNodes() - 1;
 519     while ( !IsLinked( 0, topNodeIndex, /*ignoreMediumNodes=*/true )) --topNodeIndex;
 520     const SMDS_MeshNode* topNode = myVolumeNodes[ topNodeIndex ];
 521     XYZ upDir (topNode->X() - botNode->X(),
 522                topNode->Y() - botNode->Y(),
 523                topNode->Z() - botNode->Z() );
 524     myVolForward = ( botNormal.Dot( upDir ) < 0 );
 525
 526     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   XYZ p1 ( myFaceNodes[0] );
1020   XYZ p2 ( myFaceNodes[1] );
1021   XYZ p3 ( myFaceNodes[2] );
1022   XYZ aVec12( p2 - p1 );
1023   XYZ aVec13( p3 - p1 );
1024   XYZ cross = aVec12.Crossed( aVec13 );
1025
1026   //if ( myFaceNbNodes == 4 ) {
1027   if ( myFaceNbNodes >3 ) {
1028     XYZ p4 ( myFaceNodes[3] );
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       if ( faceIndex <= 1 ) // top or bottom
1150         ind = 1 - faceIndex;
1151       else {
1152         const int nbSideFaces = myAllFacesNbNodes[0] / 2;
1153         ind = ( faceIndex - nbHoriFaces + nbSideFaces/2 ) % nbSideFaces + nbHoriFaces;
1154       }
1155       break;
1156     default:;
1157     }
1158   }
1159   return ind;
1160 }
1161
1162 //=======================================================================
1163 //function : IsLinked
1164 //purpose  : return true if theNode1 is linked with theNode2
1165 // If theIgnoreMediumNodes then corner nodes of quadratic cell are considered linked as well
1166 //=======================================================================
1167
1168 bool SMDS_VolumeTool::IsLinked (const SMDS_MeshNode* theNode1,
1169                                 const SMDS_MeshNode* theNode2,
1170                                 const bool           theIgnoreMediumNodes) const
1171 {
1172   if ( !myVolume )
1173     return false;
1174
1175   if (myVolume->IsPoly()) {
1176     if (!myPolyedre) {
1177       MESSAGE("Warning: bad volumic element");
1178       return false;
1179     }
1180     bool isLinked = false;
1181     int iface;
1182     for (iface = 1; iface <= myNbFaces && !isLinked; iface++) {
1183       int inode, nbFaceNodes = myPolyedre->NbFaceNodes(iface);
1184
1185       for (inode = 1; inode <= nbFaceNodes && !isLinked; inode++) {
1186         const SMDS_MeshNode* curNode = myPolyedre->GetFaceNode(iface, inode);
1187
1188         if (curNode == theNode1 || curNode == theNode2) {
1189           int inextnode = (inode == nbFaceNodes) ? 1 : inode + 1;
1190           const SMDS_MeshNode* nextNode = myPolyedre->GetFaceNode(iface, inextnode);
1191
1192           if ((curNode == theNode1 && nextNode == theNode2) ||
1193               (curNode == theNode2 && nextNode == theNode1)) {
1194             isLinked = true;
1195           }
1196         }
1197       }
1198     }
1199     return isLinked;
1200   }
1201
1202   // find nodes indices
1203   int i1 = -1, i2 = -1, nbFound = 0;
1204   for ( int i = 0; i < myVolumeNbNodes && nbFound < 2; i++ )
1205   {
1206     if ( myVolumeNodes[ i ] == theNode1 )
1207       i1 = i, ++nbFound;
1208     else if ( myVolumeNodes[ i ] == theNode2 )
1209       i2 = i, ++nbFound;
1210   }
1211   return IsLinked( i1, i2 );
1212 }
1213
1214 //=======================================================================
1215 //function : IsLinked
1216 //purpose  : return true if the node with theNode1Index is linked
1217 //           with the node with theNode2Index
1218 // If theIgnoreMediumNodes then corner nodes of quadratic cell are considered linked as well
1219 //=======================================================================
1220
1221 bool SMDS_VolumeTool::IsLinked (const int theNode1Index,
1222                                 const int theNode2Index,
1223                                 bool      theIgnoreMediumNodes) const
1224 {
1225   if ( myVolume->IsPoly() ) {
1226     return IsLinked(myVolumeNodes[theNode1Index], myVolumeNodes[theNode2Index]);
1227   }
1228
1229   int minInd = min( theNode1Index, theNode2Index );
1230   int maxInd = max( theNode1Index, theNode2Index );
1231
1232   if ( minInd < 0 || maxInd > myVolumeNbNodes - 1 || maxInd == minInd )
1233     return false;
1234
1235   VolumeType type = GetVolumeType();
1236   if ( myVolume->IsQuadratic() )
1237   {
1238     int firstMediumInd = myVolume->NbCornerNodes();
1239     if ( minInd >= firstMediumInd )
1240       return false; // both nodes are medium - not linked
1241     if ( maxInd < firstMediumInd ) // both nodes are corners
1242     {
1243       if ( theIgnoreMediumNodes )
1244         type = quadToLinear(type); // to check linkage of corner nodes only
1245       else
1246         return false; // corner nodes are not linked directly in a quadratic cell
1247     }
1248   }
1249
1250   switch ( type ) {
1251   case TETRA:
1252     return true;
1253   case HEXA:
1254     switch ( maxInd - minInd ) {
1255     case 1: return minInd != 3;
1256     case 3: return minInd == 0 || minInd == 4;
1257     case 4: return true;
1258     default:;
1259     }
1260     break;
1261   case PYRAM:
1262     if ( maxInd == 4 )
1263       return true;
1264     switch ( maxInd - minInd ) {
1265     case 1:
1266     case 3: return true;
1267     default:;
1268     }
1269     break;
1270   case PENTA:
1271     switch ( maxInd - minInd ) {
1272     case 1: return minInd != 2;
1273     case 2: return minInd == 0 || minInd == 3;
1274     case 3: return true;
1275     default:;
1276     }
1277     break;
1278   case QUAD_TETRA:
1279     {
1280       switch ( minInd ) {
1281       case 0: if( maxInd==4 ||  maxInd==6 ||  maxInd==7 ) return true;
1282       case 1: if( maxInd==4 ||  maxInd==5 ||  maxInd==8 ) return true;
1283       case 2: if( maxInd==5 ||  maxInd==6 ||  maxInd==9 ) return true;
1284       case 3: if( maxInd==7 ||  maxInd==8 ||  maxInd==9 ) return true;
1285       default:;
1286       }
1287       break;
1288     }
1289   case QUAD_HEXA:
1290     {
1291       switch ( minInd ) {
1292       case 0: if( maxInd==8 ||  maxInd==11 ||  maxInd==16 ) return true;
1293       case 1: if( maxInd==8 ||  maxInd==9 ||  maxInd==17 ) return true;
1294       case 2: if( maxInd==9 ||  maxInd==10 ||  maxInd==18 ) return true;
1295       case 3: if( maxInd==10 ||  maxInd==11 ||  maxInd==19 ) return true;
1296       case 4: if( maxInd==12 ||  maxInd==15 ||  maxInd==16 ) return true;
1297       case 5: if( maxInd==12 ||  maxInd==13 ||  maxInd==17 ) return true;
1298       case 6: if( maxInd==13 ||  maxInd==14 ||  maxInd==18 ) return true;
1299       case 7: if( maxInd==14 ||  maxInd==15 ||  maxInd==19 ) return true;
1300       default:;
1301       }
1302       break;
1303     }
1304   case QUAD_PYRAM:
1305     {
1306       switch ( minInd ) {
1307       case 0: if( maxInd==5 ||  maxInd==8 ||  maxInd==9 ) return true;
1308       case 1: if( maxInd==5 ||  maxInd==6 ||  maxInd==10 ) return true;
1309       case 2: if( maxInd==6 ||  maxInd==7 ||  maxInd==11 ) return true;
1310       case 3: if( maxInd==7 ||  maxInd==8 ||  maxInd==12 ) return true;
1311       case 4: if( maxInd==9 ||  maxInd==10 ||  maxInd==11 ||  maxInd==12 ) return true;
1312       default:;
1313       }
1314       break;
1315     }
1316   case QUAD_PENTA:
1317     {
1318       switch ( minInd ) {
1319       case 0: if( maxInd==6 ||  maxInd==8 ||  maxInd==12 ) return true;
1320       case 1: if( maxInd==6 ||  maxInd==7 ||  maxInd==13 ) return true;
1321       case 2: if( maxInd==7 ||  maxInd==8 ||  maxInd==14 ) return true;
1322       case 3: if( maxInd==9 ||  maxInd==11 ||  maxInd==12 ) return true;
1323       case 4: if( maxInd==9 ||  maxInd==10 ||  maxInd==13 ) return true;
1324       case 5: if( maxInd==10 ||  maxInd==11 ||  maxInd==14 ) return true;
1325       default:;
1326       }
1327       break;
1328     }
1329   case HEX_PRISM:
1330     {
1331       const int diff = maxInd-minInd;
1332       if ( diff > 6  ) return false;// not linked top and bottom
1333       if ( diff == 6 ) return true; // linked top and bottom
1334       return diff == 1 || diff == 7;
1335     }
1336   default:;
1337   }
1338   return false;
1339 }
1340
1341 //=======================================================================
1342 //function : GetNodeIndex
1343 //purpose  : Return an index of theNode
1344 //=======================================================================
1345
1346 int SMDS_VolumeTool::GetNodeIndex(const SMDS_MeshNode* theNode) const
1347 {
1348   if ( myVolume ) {
1349     for ( int i = 0; i < myVolumeNbNodes; i++ ) {
1350       if ( myVolumeNodes[ i ] == theNode )
1351         return i;
1352     }
1353   }
1354   return -1;
1355 }
1356
1357 //================================================================================
1358 /*!
1359  * \brief Fill vector with boundary faces existing in the mesh
1360   * \param faces - vector of found nodes
1361   * \retval int - nb of found faces
1362  */
1363 //================================================================================
1364
1365 int SMDS_VolumeTool::GetAllExistingFaces(vector<const SMDS_MeshElement*> & faces) const
1366 {
1367   faces.clear();
1368   for ( int iF = 0; iF < NbFaces(); ++iF ) {
1369     const SMDS_MeshFace* face = 0;
1370     const SMDS_MeshNode** nodes = GetFaceNodes( iF );
1371     switch ( NbFaceNodes( iF )) {
1372     case 3:
1373       face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2] ); break;
1374     case 4:
1375       face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2], nodes[3] ); break;
1376     case 6:
1377       face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2],
1378                                   nodes[3], nodes[4], nodes[5]); break;
1379     case 8:
1380       face = SMDS_Mesh::FindFace( nodes[0], nodes[1], nodes[2], nodes[3],
1381                                   nodes[4], nodes[5], nodes[6], nodes[7]); break;
1382     }
1383     if ( face )
1384       faces.push_back( face );
1385   }
1386   return faces.size();
1387 }
1388
1389
1390 //================================================================================
1391 /*!
1392  * \brief Fill vector with boundary edges existing in the mesh
1393   * \param edges - vector of found edges
1394   * \retval int - nb of found faces
1395  */
1396 //================================================================================
1397
1398 int SMDS_VolumeTool::GetAllExistingEdges(vector<const SMDS_MeshElement*> & edges) const
1399 {
1400   edges.clear();
1401   edges.reserve( myVolumeNbNodes * 2 );
1402   for ( int i = 0; i < myVolumeNbNodes-1; ++i ) {
1403     for ( int j = i + 1; j < myVolumeNbNodes; ++j ) {
1404       if ( IsLinked( i, j )) {
1405         const SMDS_MeshElement* edge =
1406           SMDS_Mesh::FindEdge( myVolumeNodes[i], myVolumeNodes[j] );
1407         if ( edge )
1408           edges.push_back( edge );
1409       }
1410     }
1411   }
1412   return edges.size();
1413 }
1414
1415 //================================================================================
1416 /*!
1417  * \brief Return minimal square distance between connected corner nodes
1418  */
1419 //================================================================================
1420
1421 double SMDS_VolumeTool::MinLinearSize2() const
1422 {
1423   double minSize = 1e+100;
1424   int iQ = myVolume->IsQuadratic() ? 2 : 1;
1425
1426   // store current face data
1427   int curFace = myCurFace, nbN = myFaceNbNodes;
1428   int* ind = myFaceNodeIndices;
1429   myFaceNodeIndices = NULL;
1430   const SMDS_MeshNode** nodes = myFaceNodes;
1431   myFaceNodes = NULL;
1432
1433   // it seems that compute distance twice is faster than organization of a sole computing
1434   myCurFace = -1;
1435   for ( int iF = 0; iF < myNbFaces; ++iF )
1436   {
1437     setFace( iF );
1438     for ( int iN = 0; iN < myFaceNbNodes; iN += iQ )
1439     {
1440       XYZ n1( myFaceNodes[ iN ]);
1441       XYZ n2( myFaceNodes[(iN + iQ) % myFaceNbNodes]);
1442       minSize = std::min( minSize, (n1 - n2).SquareMagnitude());
1443     }
1444   }
1445   // restore current face data
1446   myCurFace = curFace;
1447   myFaceNbNodes = nbN;
1448   myFaceNodeIndices = ind;
1449   delete [] myFaceNodes; myFaceNodes = nodes;
1450
1451   return minSize;
1452 }
1453
1454 //================================================================================
1455 /*!
1456  * \brief check that only one volume is build on the face nodes
1457  *
1458  * If a face is shared by one of <ignoreVolumes>, it is considered free
1459  */
1460 //================================================================================
1461
1462 bool SMDS_VolumeTool::IsFreeFace( int faceIndex, const SMDS_MeshElement** otherVol/*=0*/ ) const
1463 {
1464   const bool isFree = true;
1465
1466   if (!setFace( faceIndex ))
1467     return !isFree;
1468
1469   const SMDS_MeshNode** nodes = GetFaceNodes( faceIndex );
1470   const int nbFaceNodes = myFaceNbNodes;
1471
1472   // evaluate nb of face nodes shared by other volumes
1473   int maxNbShared = -1;
1474   typedef map< const SMDS_MeshElement*, int > TElemIntMap;
1475   TElemIntMap volNbShared;
1476   TElemIntMap::iterator vNbIt;
1477   for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) {
1478     const SMDS_MeshNode* n = nodes[ iNode ];
1479     SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( SMDSAbs_Volume );
1480     while ( eIt->more() ) {
1481       const SMDS_MeshElement* elem = eIt->next();
1482       if ( elem != myVolume ) {
1483         vNbIt = volNbShared.insert( make_pair( elem, 0 )).first;
1484         (*vNbIt).second++;
1485         if ( vNbIt->second > maxNbShared )
1486           maxNbShared = vNbIt->second;
1487       }
1488     }
1489   }
1490   if ( maxNbShared < 3 )
1491     return isFree; // is free
1492
1493   // find volumes laying on the opposite side of the face
1494   // and sharing all nodes
1495   XYZ intNormal; // internal normal
1496   GetFaceNormal( faceIndex, intNormal.x, intNormal.y, intNormal.z );
1497   if ( IsFaceExternal( faceIndex ))
1498     intNormal = XYZ( -intNormal.x, -intNormal.y, -intNormal.z );
1499   XYZ p0 ( nodes[0] ), baryCenter;
1500   for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end();  ) {
1501     const int& nbShared = (*vNbIt).second;
1502     if ( nbShared >= 3 ) {
1503       SMDS_VolumeTool volume( (*vNbIt).first );
1504       volume.GetBaryCenter( baryCenter.x, baryCenter.y, baryCenter.z );
1505       XYZ intNormal2( baryCenter - p0 );
1506       if ( intNormal.Dot( intNormal2 ) < 0 ) {
1507         // opposite side
1508         if ( nbShared >= nbFaceNodes )
1509         {
1510           // a volume shares the whole facet
1511           if ( otherVol ) *otherVol = vNbIt->first;
1512           return !isFree;
1513         }
1514         ++vNbIt;
1515         continue;
1516       }
1517     }
1518     // remove a volume from volNbShared map
1519     volNbShared.erase( vNbIt++ );
1520   }
1521
1522   // here volNbShared contains only volumes laying on the opposite side of
1523   // the face and sharing 3 or more but not all face nodes with myVolume
1524   if ( volNbShared.size() < 2 ) {
1525     return isFree; // is free
1526   }
1527
1528   // check if the whole area of a face is shared
1529   for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
1530   {
1531     const SMDS_MeshNode* n = nodes[ iNode ];
1532     // check if n is shared by one of volumes of volNbShared
1533     bool isShared = false;
1534     SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( SMDSAbs_Volume );
1535     while ( eIt->more() && !isShared )
1536       isShared = volNbShared.count( eIt->next() );
1537     if ( !isShared )
1538       return isFree;
1539   }
1540   if ( otherVol ) *otherVol = volNbShared.begin()->first;
1541   return !isFree;
1542
1543 //   if ( !myVolume->IsPoly() )
1544 //   {
1545 //     bool isShared[] = { false, false, false, false }; // 4 triangle parts of a quadrangle
1546 //     for ( vNbIt = volNbShared.begin(); vNbIt != volNbShared.end(); vNbIt++ ) {
1547 //       SMDS_VolumeTool volume( (*vNbIt).first );
1548 //       bool prevLinkShared = false;
1549 //       int nbSharedLinks = 0;
1550 //       for ( int iNode = 0; iNode < nbFaceNodes; iNode++ ) {
1551 //         bool linkShared = volume.IsLinked( nodes[ iNode ], nodes[ iNode + 1] );
1552 //         if ( linkShared )
1553 //           nbSharedLinks++;
1554 //         if ( linkShared && prevLinkShared &&
1555 //              volume.IsLinked( nodes[ iNode - 1 ], nodes[ iNode + 1] ))
1556 //           isShared[ iNode ] = true;
1557 //         prevLinkShared = linkShared;
1558 //       }
1559 //       if ( nbSharedLinks == nbFaceNodes )
1560 //         return !free; // is not free
1561 //       if ( nbFaceNodes == 4 ) {
1562 //         // check traingle parts 1 & 3
1563 //         if ( isShared[1] && isShared[3] )
1564 //           return !free; // is not free
1565 //         // check triangle parts 0 & 2;
1566 //         // 0 part could not be checked in the loop; check it here
1567 //         if ( isShared[2] && prevLinkShared &&
1568 //              volume.IsLinked( nodes[ 0 ], nodes[ 1 ] ) &&
1569 //              volume.IsLinked( nodes[ 1 ], nodes[ 3 ] ) )
1570 //           return !free; // is not free
1571 //       }
1572 //     }
1573 //   }
1574 //  return free;
1575 }
1576
1577 //=======================================================================
1578 //function : GetFaceIndex
1579 //purpose  : Return index of a face formed by theFaceNodes
1580 //=======================================================================
1581
1582 int SMDS_VolumeTool::GetFaceIndex( const set<const SMDS_MeshNode*>& theFaceNodes ) const
1583 {
1584   for ( int iFace = 0; iFace < myNbFaces; iFace++ ) {
1585     const SMDS_MeshNode** nodes = GetFaceNodes( iFace );
1586     int nbFaceNodes = NbFaceNodes( iFace );
1587     set<const SMDS_MeshNode*> nodeSet;
1588     for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
1589       nodeSet.insert( nodes[ iNode ] );
1590     if ( theFaceNodes == nodeSet )
1591       return iFace;
1592   }
1593   return -1;
1594 }
1595
1596 //=======================================================================
1597 //function : GetFaceIndex
1598 //purpose  : Return index of a face formed by theFaceNodes
1599 //=======================================================================
1600
1601 /*int SMDS_VolumeTool::GetFaceIndex( const set<int>& theFaceNodesIndices )
1602 {
1603   for ( int iFace = 0; iFace < myNbFaces; iFace++ ) {
1604     const int* nodes = GetFaceNodesIndices( iFace );
1605     int nbFaceNodes = NbFaceNodes( iFace );
1606     set<int> nodeSet;
1607     for ( int iNode = 0; iNode < nbFaceNodes; iNode++ )
1608       nodeSet.insert( nodes[ iNode ] );
1609     if ( theFaceNodesIndices == nodeSet )
1610       return iFace;
1611   }
1612   return -1;
1613 }*/
1614
1615 //=======================================================================
1616 //function : setFace
1617 //purpose  :
1618 //=======================================================================
1619
1620 bool SMDS_VolumeTool::setFace( int faceIndex ) const
1621 {
1622   if ( !myVolume )
1623     return false;
1624
1625   if ( myCurFace == faceIndex )
1626     return true;
1627
1628   myCurFace = -1;
1629
1630   if ( faceIndex < 0 || faceIndex >= NbFaces() )
1631     return false;
1632
1633   if (myFaceNodes != NULL) {
1634     delete [] myFaceNodes;
1635     myFaceNodes = NULL;
1636   }
1637
1638   if (myVolume->IsPoly())
1639   {
1640     if (!myPolyedre) {
1641       MESSAGE("Warning: bad volumic element");
1642       return false;
1643     }
1644
1645     // set face nodes
1646     int iNode;
1647     myFaceNbNodes = myPolyedre->NbFaceNodes(faceIndex + 1);
1648     myFaceNodes = new const SMDS_MeshNode* [myFaceNbNodes + 1];
1649     for ( iNode = 0; iNode < myFaceNbNodes; iNode++ )
1650       myFaceNodes[ iNode ] = myPolyedre->GetFaceNode(faceIndex + 1, iNode + 1);
1651     myFaceNodes[ myFaceNbNodes ] = myFaceNodes[ 0 ]; // last = first
1652
1653     // check orientation
1654     if (myExternalFaces)
1655     {
1656       myCurFace = faceIndex; // avoid infinite recursion in IsFaceExternal()
1657       myExternalFaces = false; // force normal computation by IsFaceExternal()
1658       if ( !IsFaceExternal( faceIndex ))
1659         for ( int i = 0, j = myFaceNbNodes; i < j; ++i, --j )
1660           std::swap( myFaceNodes[i], myFaceNodes[j] );
1661       myExternalFaces = true;
1662     }
1663   }
1664   else
1665   {
1666     if ( !myAllFacesNodeIndices_F )
1667     {
1668       // choose data for an element type
1669       switch ( myVolumeNbNodes ) {
1670       case 4:
1671         myAllFacesNodeIndices_F  = &Tetra_F [0][0];
1672         //myAllFacesNodeIndices_FE = &Tetra_F [0][0];
1673         myAllFacesNodeIndices_RE = &Tetra_RE[0][0];
1674         myAllFacesNbNodes        = Tetra_nbN;
1675         myMaxFaceNbNodes         = sizeof(Tetra_F[0])/sizeof(Tetra_F[0][0]);
1676         break;
1677       case 5:
1678         myAllFacesNodeIndices_F  = &Pyramid_F [0][0];
1679         //myAllFacesNodeIndices_FE = &Pyramid_F [0][0];
1680         myAllFacesNodeIndices_RE = &Pyramid_RE[0][0];
1681         myAllFacesNbNodes        = Pyramid_nbN;
1682         myMaxFaceNbNodes         = sizeof(Pyramid_F[0])/sizeof(Pyramid_F[0][0]);
1683         break;
1684       case 6:
1685         myAllFacesNodeIndices_F  = &Penta_F [0][0];
1686         //myAllFacesNodeIndices_FE = &Penta_FE[0][0];
1687         myAllFacesNodeIndices_RE = &Penta_RE[0][0];
1688         myAllFacesNbNodes        = Penta_nbN;
1689         myMaxFaceNbNodes         = sizeof(Penta_F[0])/sizeof(Penta_F[0][0]);
1690         break;
1691       case 8:
1692         myAllFacesNodeIndices_F  = &Hexa_F [0][0];
1693         ///myAllFacesNodeIndices_FE = &Hexa_FE[0][0];
1694         myAllFacesNodeIndices_RE = &Hexa_RE[0][0];
1695         myAllFacesNbNodes        = Hexa_nbN;
1696         myMaxFaceNbNodes         = sizeof(Hexa_F[0])/sizeof(Hexa_F[0][0]);
1697         break;
1698       case 10:
1699         myAllFacesNodeIndices_F  = &QuadTetra_F [0][0];
1700         //myAllFacesNodeIndices_FE = &QuadTetra_F [0][0];
1701         myAllFacesNodeIndices_RE = &QuadTetra_RE[0][0];
1702         myAllFacesNbNodes        = QuadTetra_nbN;
1703         myMaxFaceNbNodes         = sizeof(QuadTetra_F[0])/sizeof(QuadTetra_F[0][0]);
1704         break;
1705       case 13:
1706         myAllFacesNodeIndices_F  = &QuadPyram_F [0][0];
1707         //myAllFacesNodeIndices_FE = &QuadPyram_F [0][0];
1708         myAllFacesNodeIndices_RE = &QuadPyram_RE[0][0];
1709         myAllFacesNbNodes        = QuadPyram_nbN;
1710         myMaxFaceNbNodes         = sizeof(QuadPyram_F[0])/sizeof(QuadPyram_F[0][0]);
1711         break;
1712       case 15:
1713         myAllFacesNodeIndices_F  = &QuadPenta_F [0][0];
1714         //myAllFacesNodeIndices_FE = &QuadPenta_FE[0][0];
1715         myAllFacesNodeIndices_RE = &QuadPenta_RE[0][0];
1716         myAllFacesNbNodes        = QuadPenta_nbN;
1717         myMaxFaceNbNodes         = sizeof(QuadPenta_F[0])/sizeof(QuadPenta_F[0][0]);
1718         break;
1719       case 20:
1720       case 27:
1721         myAllFacesNodeIndices_F  = &QuadHexa_F [0][0];
1722         //myAllFacesNodeIndices_FE = &QuadHexa_FE[0][0];
1723         myAllFacesNodeIndices_RE = &QuadHexa_RE[0][0];
1724         myAllFacesNbNodes        = QuadHexa_nbN;
1725         myMaxFaceNbNodes         = sizeof(QuadHexa_F[0])/sizeof(QuadHexa_F[0][0]);
1726         if ( !myIgnoreCentralNodes && myVolumeNbNodes == 27 )
1727         {
1728           myAllFacesNodeIndices_F  = &TriQuadHexa_F [0][0];
1729           //myAllFacesNodeIndices_FE = &TriQuadHexa_FE[0][0];
1730           myAllFacesNodeIndices_RE = &TriQuadHexa_RE[0][0];
1731           myAllFacesNbNodes        = TriQuadHexa_nbN;
1732           myMaxFaceNbNodes         = sizeof(TriQuadHexa_F[0])/sizeof(TriQuadHexa_F[0][0]);
1733         }
1734         break;
1735       case 12:
1736         myAllFacesNodeIndices_F  = &HexPrism_F [0][0];
1737         //myAllFacesNodeIndices_FE = &HexPrism_FE[0][0];
1738         myAllFacesNodeIndices_RE = &HexPrism_RE[0][0];
1739         myAllFacesNbNodes        = HexPrism_nbN;
1740         myMaxFaceNbNodes         = sizeof(HexPrism_F[0])/sizeof(HexPrism_F[0][0]);
1741         break;
1742       default:
1743         return false;
1744       }
1745     }
1746     myFaceNbNodes = myAllFacesNbNodes[ faceIndex ];
1747     // if ( myExternalFaces )
1748     //   myFaceNodeIndices = (int*)( myVolForward ? myAllFacesNodeIndices_FE + faceIndex*myMaxFaceNbNodes : myAllFacesNodeIndices_RE + faceIndex*myMaxFaceNbNodes );
1749     // else
1750     //   myFaceNodeIndices = (int*)( myAllFacesNodeIndices_F + faceIndex*myMaxFaceNbNodes );
1751     myFaceNodeIndices = (int*)( myVolForward ? myAllFacesNodeIndices_F + faceIndex*myMaxFaceNbNodes : myAllFacesNodeIndices_RE + faceIndex*myMaxFaceNbNodes );
1752
1753     // set face nodes
1754     myFaceNodes = new const SMDS_MeshNode* [myFaceNbNodes + 1];
1755     for ( int iNode = 0; iNode < myFaceNbNodes; iNode++ )
1756       myFaceNodes[ iNode ] = myVolumeNodes[ myFaceNodeIndices[ iNode ]];
1757     myFaceNodes[ myFaceNbNodes ] = myFaceNodes[ 0 ];
1758   }
1759
1760   myCurFace = faceIndex;
1761
1762   return true;
1763 }
1764
1765 //=======================================================================
1766 //function : GetType
1767 //purpose  : return VolumeType by nb of nodes in a volume
1768 //=======================================================================
1769
1770 SMDS_VolumeTool::VolumeType SMDS_VolumeTool::GetType(int nbNodes)
1771 {
1772   switch ( nbNodes ) {
1773   case 4: return TETRA;
1774   case 5: return PYRAM;
1775   case 6: return PENTA;
1776   case 8: return HEXA;
1777   case 10: return QUAD_TETRA;
1778   case 13: return QUAD_PYRAM;
1779   case 15: return QUAD_PENTA;
1780   case 20:
1781   case 27: return QUAD_HEXA;
1782   case 12: return HEX_PRISM;
1783   default:return UNKNOWN;
1784   }
1785 }
1786
1787 //=======================================================================
1788 //function : NbFaces
1789 //purpose  : return nb of faces by volume type
1790 //=======================================================================
1791
1792 int SMDS_VolumeTool::NbFaces( VolumeType type )
1793 {
1794   switch ( type ) {
1795   case TETRA     :
1796   case QUAD_TETRA: return 4;
1797   case PYRAM     :
1798   case QUAD_PYRAM: return 5;
1799   case PENTA     :
1800   case QUAD_PENTA: return 5;
1801   case HEXA      :
1802   case QUAD_HEXA : return 6;
1803   case HEX_PRISM : return 8;
1804   default:         return 0;
1805   }
1806 }
1807
1808 //================================================================================
1809 /*!
1810  * \brief Useful to know nb of corner nodes of a quadratic volume
1811   * \param type - volume type
1812   * \retval int - nb of corner nodes
1813  */
1814 //================================================================================
1815
1816 int SMDS_VolumeTool::NbCornerNodes(VolumeType type)
1817 {
1818   switch ( type ) {
1819   case TETRA     :
1820   case QUAD_TETRA: return 4;
1821   case PYRAM     :
1822   case QUAD_PYRAM: return 5;
1823   case PENTA     :
1824   case QUAD_PENTA: return 6;
1825   case HEXA      :
1826   case QUAD_HEXA : return 8;
1827   case HEX_PRISM : return 12;
1828   default:         return 0;
1829   }
1830   return 0;
1831 }
1832   //
1833
1834 //=======================================================================
1835 //function : GetFaceNodesIndices
1836 //purpose  : Return the array of face nodes indices
1837 //           To comfort link iteration, the array
1838 //           length == NbFaceNodes( faceIndex ) + 1 and
1839 //           the last node index == the first one.
1840 //=======================================================================
1841
1842 const int* SMDS_VolumeTool::GetFaceNodesIndices(VolumeType type,
1843                                                 int        faceIndex,
1844                                                 bool       external)
1845 {
1846   switch ( type ) {
1847   case TETRA: return Tetra_F[ faceIndex ];
1848   case PYRAM: return Pyramid_F[ faceIndex ];
1849   case PENTA: return external ? Penta_F[ faceIndex ] : Penta_F[ faceIndex ];
1850   case HEXA:  return external ? Hexa_F[ faceIndex ] : Hexa_F[ faceIndex ];
1851   case QUAD_TETRA: return QuadTetra_F[ faceIndex ];
1852   case QUAD_PYRAM: return QuadPyram_F[ faceIndex ];
1853   case QUAD_PENTA: return external ? QuadPenta_F[ faceIndex ] : QuadPenta_F[ faceIndex ];
1854     // what about SMDSEntity_TriQuad_Hexa?
1855   case QUAD_HEXA:  return external ? QuadHexa_F[ faceIndex ] : QuadHexa_F[ faceIndex ];
1856   case HEX_PRISM:  return external ? HexPrism_F[ faceIndex ] : HexPrism_F[ faceIndex ];
1857   default:;
1858   }
1859   return 0;
1860 }
1861
1862 //=======================================================================
1863 //function : NbFaceNodes
1864 //purpose  : Return number of nodes in the array of face nodes
1865 //=======================================================================
1866
1867 int SMDS_VolumeTool::NbFaceNodes(VolumeType type,
1868                                  int        faceIndex )
1869 {
1870   switch ( type ) {
1871   case TETRA: return Tetra_nbN[ faceIndex ];
1872   case PYRAM: return Pyramid_nbN[ faceIndex ];
1873   case PENTA: return Penta_nbN[ faceIndex ];
1874   case HEXA:  return Hexa_nbN[ faceIndex ];
1875   case QUAD_TETRA: return QuadTetra_nbN[ faceIndex ];
1876   case QUAD_PYRAM: return QuadPyram_nbN[ faceIndex ];
1877   case QUAD_PENTA: return QuadPenta_nbN[ faceIndex ];
1878     // what about SMDSEntity_TriQuad_Hexa?
1879   case QUAD_HEXA:  return QuadHexa_nbN[ faceIndex ];
1880   case HEX_PRISM:  return HexPrism_nbN[ faceIndex ];
1881   default:;
1882   }
1883   return 0;
1884 }
1885
1886 //=======================================================================
1887 //function : Element
1888 //purpose  : return element
1889 //=======================================================================
1890
1891 const SMDS_MeshVolume* SMDS_VolumeTool::Element() const
1892 {
1893   return static_cast<const SMDS_MeshVolume*>( myVolume );
1894 }
1895
1896 //=======================================================================
1897 //function : ID
1898 //purpose  : return element ID
1899 //=======================================================================
1900
1901 int SMDS_VolumeTool::ID() const
1902 {
1903   return myVolume ? myVolume->GetID() : 0;
1904 }